From a956e69893366d69ce831a3b3d1d1da2af9c9c89 Mon Sep 17 00:00:00 2001 From: Ian Bell Date: Wed, 10 Aug 2022 12:49:08 -0400 Subject: [PATCH] Add implementation of polishing of pure fluid endpoint in critical tracing See https://github.com/usnistgov/teqp/issues/9 --- include/teqp/algorithms/critical_tracing.hpp | 47 +++++++++++++++++++- interface/pybind11_wrapper.cpp | 1 + interface/pybind11_wrapper.hpp | 2 +- interface/teqpversion.hpp | 2 +- src/tests/catch_tests.cxx | 7 ++- 5 files changed, 55 insertions(+), 4 deletions(-) diff --git a/include/teqp/algorithms/critical_tracing.hpp b/include/teqp/algorithms/critical_tracing.hpp index 2d8c96d5..e776d559 100644 --- a/include/teqp/algorithms/critical_tracing.hpp +++ b/include/teqp/algorithms/critical_tracing.hpp @@ -7,6 +7,7 @@ #include #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 @@ -467,6 +469,10 @@ struct CriticalTracing { auto extract_drhodt = [](const state_type& dxdt) -> Eigen::ArrayXd { return Eigen::Map(&(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; } diff --git a/interface/pybind11_wrapper.cpp b/interface/pybind11_wrapper.cpp index c9bafa43..64b7c3fd 100644 --- a/interface/pybind11_wrapper.cpp +++ b/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) ; diff --git a/interface/pybind11_wrapper.hpp b/interface/pybind11_wrapper.hpp index ed2d21d6..057ab8c4 100644 --- a/interface/pybind11_wrapper.hpp +++ b/interface/pybind11_wrapper.hpp @@ -74,7 +74,7 @@ void add_derivatives(py::module &m, Wrapper &cls) { m.def("extrapolate_from_critical", &extrapolate_from_critical); m.def("pure_VLE_T", &pure_VLE_T); - m.def("get_pure_critical_conditions_Jacobian", &get_pure_critical_conditions_Jacobian, py::arg("model"), py::arg("T"), py::arg("rho")); + m.def("get_pure_critical_conditions_Jacobian", &get_pure_critical_conditions_Jacobian, 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, py::arg("model"), py::arg("T"), py::arg("rho"), py::arg_v("flags", std::nullopt, "None")); m.def("mix_VLE_Tx", &mix_VLE_Tx); m.def("mixture_VLE_px", &mixture_VLE_px, 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")); diff --git a/interface/teqpversion.hpp b/interface/teqpversion.hpp index 0f357825..d1148d8b 100644 --- a/interface/teqpversion.hpp +++ b/interface/teqpversion.hpp @@ -1,2 +1,2 @@ #include -const std::string TEQPVERSION = "0.10.0"; \ No newline at end of file +const std::string TEQPVERSION = "0.11.0.dev0"; \ No newline at end of file diff --git a/src/tests/catch_tests.cxx b/src/tests/catch_tests.cxx index 34fb1185..3c8bc7b2 100644 --- a/src/tests/catch_tests.cxx +++ b/src/tests/catch_tests.cxx @@ -247,6 +247,9 @@ TEST_CASE("Trace critical locus for vdW", "[vdW][crit]") vdWEOS vdW(Tc_K, pc_Pa); auto Zc = 3.0/8.0; std::valarray 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; 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));