Add implementation of polishing of pure fluid endpoint in critical tr…

…acing

See #9

include/teqp/algorithms/critical_tracing.hpp

@@ -7,6 +7,7 @@
 
 #include <Eigen/Dense>
 #include "teqp/algorithms/rootfinding.hpp"
+#include "teqp/algorithms/critical_pure.hpp"
 #include "teqp/exceptions.hpp"
 
 // Imports from boost
@@ -33,6 +34,7 @@ struct TCABOptions {
     double stability_rel_drho = 0.001; ///< The relative size of the step (relative to the sum of the molar concentration vector) to be used when taking the step in the direction of \f$\sigma_1\f$ when assessing local stability
     int verbosity = 0; ///< The greater the verbosity, the more output you will get, especially about polishing failures
     bool polish_exception_on_fail = false; ///< If true, when polishing fails, throw an exception, otherwise, terminate tracing
+    bool pure_endpoint_polish = true; ///< If true, if the last step crossed into negative concentrations, try to interpolate to find the pure fluid endpoint hiding in the data
 };
 
 template<typename Model, typename Scalar = double, typename VecType = Eigen::ArrayXd>
@@ -467,6 +469,10 @@ struct CriticalTracing {
         auto extract_drhodt = [](const state_type& dxdt) -> Eigen::ArrayXd {
             return Eigen::Map<const Eigen::ArrayXd>(&(dxdt[0]) + 1, dxdt.size() - 1);
         };
+        // Pull out the dTdt from dxdt, just a convenience function
+        auto extract_dTdt = [](const state_type& dxdt) -> double {
+            return dxdt[0];
+        };
 
         // Define the tolerances
         double abs_err = options.abs_err, rel_err = options.rel_err, a_x = 1.0, a_dxdt = 1.0;
@@ -567,7 +573,9 @@ struct CriticalTracing {
                     res = controlled_stepper.try_step(xprime, x0, t, dt);
                 }
                 catch (const std::exception &e) {
-                    std::cout << e.what() << std::endl;
+                    if (options.verbosity > 0) {
+                        std::cout << e.what() << std::endl;
+                    }
                     break;
                 }
 
@@ -701,6 +709,43 @@ struct CriticalTracing {
                 break;
             }
         }
+        // If the last step crosses a zero concentration, see if it corresponds to a pure fluid
+        // and if so, iterate to find the pure fluid endpoint
+        if (options.pure_endpoint_polish) {
+            // Simple Euler step t
+            auto dxdt = get_dxdt(x0);
+            auto drhodt = extract_drhodt(dxdt);
+            const auto step = (rhovec + drhodt*dt).eval();
+            if ((step * rhovec > 0).any()) {
+                // Calculate step sizes to take concentrations to zero
+                auto step_sizes = ((-rhovec) / drhodt).eval();
+                Eigen::Index ipure;
+                rhovec.maxCoeff(&ipure);
+                // Find new step size based on the pure we are heading towards
+                auto new_step_size = step_sizes(ipure);
+                // Take an Euler step of the new size
+                const auto new_rhovec = (rhovec + drhodt*new_step_size).eval();
+                const double new_T = T + extract_dTdt(dxdt)*new_step_size;
+
+                // Solve for pure fluid critical point
+                // 
+                // A bit trickier here because we need to hack the pure fluid model to put the pure fluid 
+                // composition in an index other than the first if ipure != 0 because we don't want 
+                // to require that a pure fluid model is also provided since zero mole fractions
+                // should be fine
+                nlohmann::json flags = { {"alternative_pure_index", ipure}, {"alternative_length", 2} };
+                auto [TT, rhorho] = solve_pure_critical(model, new_T, new_rhovec.sum(), flags);
+
+                // Replace the values in T and rhovec
+                T = TT;
+                rhovec[ipure] = rhorho;
+                rhovec[1 - ipure] = 0;
+
+                // And store the polished values
+                if (!filename.empty()) { write_line(); }
+                store_point();
+            }
+        }
         //auto N = JSONdata.size();
         return JSONdata;
     }

interface/pybind11_wrapper.cpp

@@ -61,6 +61,7 @@ void init_teqp(py::module& m) {
         .def_readwrite("stability_rel_drho", &TCABOptions::stability_rel_drho)
         .def_readwrite("verbosity", &TCABOptions::verbosity)
         .def_readwrite("polish", &TCABOptions::polish)
+        .def_readwrite("pure_endpoint_polish", &TCABOptions::pure_endpoint_polish)
         .def_readwrite("polish_exception_on_fail", &TCABOptions::polish_exception_on_fail)
         ;
 

interface/pybind11_wrapper.hpp

@@ -74,7 +74,7 @@ void add_derivatives(py::module &m, Wrapper &cls) {
     m.def("extrapolate_from_critical", &extrapolate_from_critical<Model, double>);
     m.def("pure_VLE_T", &pure_VLE_T<Model, double>);
 
-    m.def("get_pure_critical_conditions_Jacobian", &get_pure_critical_conditions_Jacobian<Model, double, ADBackends::autodiff>, py::arg("model"), py::arg("T"), py::arg("rho"));
+    m.def("get_pure_critical_conditions_Jacobian", &get_pure_critical_conditions_Jacobian<Model, double, ADBackends::autodiff>, py::arg("model"), py::arg("T"), py::arg("rho"), py::arg_v("alternative_pure_index", -1), py::arg_v("alternative_length", 2));
     m.def("solve_pure_critical", &solve_pure_critical<Model, double, ADBackends::autodiff>, py::arg("model"), py::arg("T"), py::arg("rho"), py::arg_v("flags", std::nullopt, "None"));
     m.def("mix_VLE_Tx", &mix_VLE_Tx<Model, double, Eigen::ArrayXd>);
     m.def("mixture_VLE_px", &mixture_VLE_px<Model, double, Eigen::ArrayXd>, py::arg("model"), py::arg("p_spec"), py::arg("xmolar_spec").noconvert(), py::arg("T0"), py::arg("rhovecL0").noconvert(), py::arg("rhovecV0").noconvert(), py::arg_v("flags", std::nullopt, "None"));

interface/teqpversion.hpp

@@ -1,2 +1,2 @@
 #include <string>
-const std::string TEQPVERSION = "0.10.0";
+const std::string TEQPVERSION = "0.11.0.dev0";

src/tests/catch_tests.cxx

@@ -247,6 +247,9 @@ TEST_CASE("Trace critical locus for vdW", "[vdW][crit]")
     vdWEOS<double> vdW(Tc_K, pc_Pa);
     auto Zc = 3.0/8.0;
     std::valarray<double> max_spluses(Tc_K.size());
+    auto rhoc0 = pc_Pa[0] / (vdW.R(molefrac) * Tc_K[0]) / Zc;
+    auto rhoc1 = pc_Pa[1] / (vdW.R(molefrac) * Tc_K[1]) / Zc;
+
     for (auto ifluid = 0; ifluid < Tc_K.size(); ++ifluid) {
         auto rhoc0 = pc_Pa[ifluid] / (vdW.R(molefrac) * Tc_K[ifluid]) / Zc;
         double T0 = Tc_K[ifluid];
@@ -258,7 +261,7 @@ TEST_CASE("Trace critical locus for vdW", "[vdW][crit]")
         REQUIRE(splus == Approx(-log(1 - 1.0 / 3.0)));
 
         auto tic0 = std::chrono::steady_clock::now();
-        std::string filename = "";
+        std::string filename = "ArNe";
         using ct = CriticalTracing<decltype(vdW), double, Eigen::ArrayXd>;
         TCABOptions opt;
         opt.polish = true;
@@ -270,6 +273,8 @@ TEST_CASE("Trace critical locus for vdW", "[vdW][crit]")
             max_splus = std::max(max_splus, splus);
         }
         max_spluses[ifluid] = max_splus;
+
+        CHECK(trace.back().at("T / K") == Approx(Tc_K[1 - ifluid]));
     }
     CHECK(max_spluses.min() == Approx(max_spluses.max()).epsilon(0.01));
     CHECK(max_spluses.min() > -log(1 - 1.0 / 3.0));