diff --git a/include/teqp/algorithms/VLE.hpp b/include/teqp/algorithms/VLE.hpp
index 346a05b021e4b20e5c31f60ec588c402ca39a1b9..25f9f690550051cdd2a40ae4e25a7cc8c2858794 100644
--- a/include/teqp/algorithms/VLE.hpp
+++ b/include/teqp/algorithms/VLE.hpp
@@ -1,8 +1,14 @@
#pragma once
+#include <optional>
#include "teqp/derivs.hpp"
#include <Eigen/Dense>
+// Imports from boost for numerical integration
+#include <boost/numeric/odeint/stepper/controlled_runge_kutta.hpp>
+#include <boost/numeric/odeint/stepper/runge_kutta_cash_karp54.hpp>
+#include <boost/numeric/odeint/stepper/euler.hpp>
+
template<typename Model, typename TYPE = double>
class IsothermPureVLEResiduals {
typedef Eigen::Array<TYPE, 2, 1> EigenArray;
@@ -207,4 +213,189 @@ auto get_drhovecdp_Tsat(const Model& model, const Scalar &T, const VecType& rhov
drhodp_vap = linsolve(PSIVstar, PSILstar*drhodp_liq);
}
return std::make_tuple(drhodp_liq, drhodp_vap);
+}
+
+struct TVLEOptions {
+ double init_dt = 1e-5, abs_err = 1e-8, rel_err = 1e-8, max_dt = 100000;
+ int max_steps = 1000, integration_order = 5;
+ bool polish = false;
+};
+
+/***
+* \brief Trace an isotherm with parametric tracing
+*/
+template<typename Model, typename Scalar, typename VecType>
+auto trace_VLE_isotherm_binary(const Model &model, Scalar T, VecType rhovecL0, VecType rhovecV0, const std::optional<TVLEOptions>& options = std::nullopt)
+{
+ // Get the options, or the default values if not provided
+ TVLEOptions opt = options.value_or(TVLEOptions{});
+ auto N = rhovecL0.size();
+ if (N != 2) {
+ std::invalid_argument("Size must be 2");
+ }
+ if (rhovecL0.size() != rhovecV0.size()) {
+ std::invalid_argument("Both molar concentration arrays must be of the same size");
+ }
+
+ auto norm = [](const auto& v) { return sqrt((v * v).sum()); };
+
+ // Define datatypes and functions for tracing tools
+ auto JSONdata = nlohmann::json::array();
+
+ // Typedefs for the types
+ using namespace boost::numeric::odeint;
+ using state_type = std::vector<double>;
+
+ // Class for simple Euler integration
+ euler<state_type> eul;
+ typedef runge_kutta_cash_karp54< state_type > error_stepper_type;
+ typedef controlled_runge_kutta< error_stepper_type > controlled_stepper_type;
+
+ // Define the tolerances
+ double abs_err = opt.abs_err, rel_err = opt.rel_err, a_x = 1.0, a_dxdt = 1.0;
+ controlled_stepper_type controlled_stepper(default_error_checker< double, range_algebra, default_operations >(abs_err, rel_err, a_x, a_dxdt));
+
+
+ double c = 1.0;
+
+ // The function to be integrated by odeint
+ auto xprime = [&](const state_type& X, state_type& Xprime, double /*t*/) {
+ // Memory maps into the state vector for inputs and their derivatives
+ // These are views, not copies!
+ auto rhovecL = Eigen::Map<const Eigen::ArrayXd>(&(X[0]), N);
+ auto rhovecV = Eigen::Map<const Eigen::ArrayXd>(&(X[0]) + N, N);
+ auto drhovecdtL = Eigen::Map<Eigen::ArrayXd>(&(Xprime[0]), N);
+ auto drhovecdtV = Eigen::Map<Eigen::ArrayXd>(&(Xprime[0]) + N, N);
+ // Get the derivatives with respect to pressure along the isotherm of the phase envelope
+ auto [drhovecdpL, drhovecdpV] = get_drhovecdp_Tsat(model, T, rhovecL, rhovecV);
+ // Get the derivative of p w.r.t. parameter
+ auto dpdt = 1.0/sqrt(norm(drhovecdpL.array()) + norm(drhovecdpV.array()));
+ // And finally the derivatives with respect to the tracing variable
+ drhovecdtL = c*drhovecdpL*dpdt;
+ drhovecdtV = c*drhovecdpV*dpdt;
+ };
+
+ // Set up the initial state vector
+ state_type x0(2*N), last_drhodt(2*N);
+ auto set_init_state = [&](state_type& X) {
+ auto rhovecL = Eigen::Map<Eigen::ArrayXd>(&(X[0]), N);
+ auto rhovecV = Eigen::Map<Eigen::ArrayXd>(&(X[0]) + N, N);
+ rhovecL = rhovecL0;
+ rhovecV = rhovecV0;
+ };
+ set_init_state(x0);
+
+ // Figure out which direction to trace initially
+ double t = 0, dt = opt.init_dt;
+ {
+ auto dxdt = x0;
+ xprime(x0, dxdt, -1.0);
+ const auto dXdt = Eigen::Map<const Eigen::ArrayXd>(&(dxdt[0]), dxdt.size());
+ const auto X = Eigen::Map<const Eigen::ArrayXd>(&(x0[0]), x0.size());
+
+ const Eigen::ArrayXd step = X + dXdt * dt;
+ Eigen::ArrayX<bool> negativestepvals = (step < 0).eval();
+ // Flip the sign if the first step would yield any negative concentrations
+ if (negativestepvals.any()) {
+ c *= -1;
+ }
+ }
+ std::string termination_reason;
+
+ // Then trace...
+ int retry_count = 0;
+ for (auto istep = 0; istep < opt.max_steps; ++istep) {
+
+ auto store_point = [&]() {
+ //// Calculate some other parameters, for debugging
+ auto N = x0.size() / 2;
+ auto rhovecL = Eigen::Map<const Eigen::ArrayXd>(&(x0[0]), N);
+ auto rhovecV = Eigen::Map<const Eigen::ArrayXd>(&(x0[0]) + N, N);
+ auto rhotot = rhovecL.sum();
+ using id = IsochoricDerivatives<decltype(model), Scalar, VecType>;
+ double pL = rhotot * model.R(rhovecL / rhovecL.sum())*T + id::get_pr(model, T, rhovecL);
+
+ // Store the data in a JSON structure
+ nlohmann::json point = {
+ {"t", t},
+ {"T / K", T},
+ {"pL / Pa", pL},
+ {"c", c},
+ {"rhoL / mol/m^3", rhovecL},
+ {"rhoV / mol/m^3", rhovecV},
+
+ };
+ JSONdata.push_back(point);
+ };
+ if (istep == 0 && retry_count == 0) {
+ store_point();
+ }
+
+ double dtold = dt;
+ auto x0_previous = x0;
+
+ // Call the derivative function to get drho/dt in order
+ // to cache it at the *beginning* of the step
+ auto dxdt = x0;
+ try {
+ xprime(x0, dxdt, -1.0);
+ }
+ catch (...) {
+ break;
+ }
+ last_drhodt = dxdt;
+
+ if (opt.integration_order == 5) {
+ controlled_step_result res = controlled_step_result::fail;
+ try {
+ res = controlled_stepper.try_step(xprime, x0, t, dt);
+ }
+ catch (...) {
+ break;
+ }
+
+ if (res != controlled_step_result::success) {
+ // Try again, with a smaller step size
+ istep--;
+ retry_count++;
+ continue;
+ }
+ else {
+ retry_count = 0;
+ }
+ // Reduce step size if greater than the specified max step size
+ dt = std::min(dt, opt.max_dt);
+ }
+ else if (opt.integration_order == 1) {
+ try {
+ eul.do_step(xprime, x0, t, dt);
+ }
+ catch (...) {
+ break;
+ }
+ }
+ else {
+ throw std::invalid_argument("integration order is invalid:" + std::to_string(opt.integration_order));
+ }
+ auto stop_requested = [&]() {
+ //// Calculate some other parameters, for debugging
+ auto N = x0.size() / 2;
+ auto rhovecL = Eigen::Map<const Eigen::ArrayXd>(&(x0[0]), N);
+ auto rhovecV = Eigen::Map<const Eigen::ArrayXd>(&(x0[0]) + N, N);
+ auto x = rhovecL / rhovecL.sum();
+ auto y = rhovecV / rhovecV.sum();
+ if ((x < 0).any() || (x > 1).any() || (y < 0).any() || (y > 1).any()) {
+ return true;
+ }
+ else {
+ return false;
+ }
+ };
+ if (stop_requested()) {
+ break;
+ }
+
+ store_point();
+ }
+ return JSONdata;
}
\ No newline at end of file
diff --git a/interface/pybind11_wrapper.cpp b/interface/pybind11_wrapper.cpp
index 0750c84900c7e70283dc8dc57e8e62cbcc1f12b3..910b0b94cd77c4c50d0c7dd1f8e1e75d9a50700c 100644
--- a/interface/pybind11_wrapper.cpp
+++ b/interface/pybind11_wrapper.cpp
@@ -38,6 +38,17 @@ void init_teqp(py::module& m) {
.def_readwrite("integration_order", &TCABOptions::integration_order)
;
+ // The options class for isotherm tracer, not tied to a particular model
+ py::class_<TVLEOptions>(m, "TVLEOptions")
+ .def(py::init<>())
+ .def_readwrite("abs_err", &TVLEOptions::abs_err)
+ .def_readwrite("rel_err", &TVLEOptions::rel_err)
+ .def_readwrite("init_dt", &TVLEOptions::init_dt)
+ .def_readwrite("max_dt", &TVLEOptions::max_dt)
+ .def_readwrite("max_steps", &TVLEOptions::max_steps)
+ .def_readwrite("integration_order", &TVLEOptions::integration_order)
+ ;
+
// Some functions for timing overhead of interface
m.def("___mysummer", [](const double &c, const Eigen::ArrayXd &x) { return c*x.sum(); });
using RAX = Eigen::Ref<Eigen::ArrayXd>;
diff --git a/interface/pybind11_wrapper.hpp b/interface/pybind11_wrapper.hpp
index dff601d79b6ed1d02f3b900cc4c4c210bbad95af..cc11df5f76b4daa8073e245d208626b98774c421 100644
--- a/interface/pybind11_wrapper.hpp
+++ b/interface/pybind11_wrapper.hpp
@@ -51,6 +51,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_drhovecdp_Tsat", &get_drhovecdp_Tsat<Model, double, RAX>, py::arg("model"), py::arg("T"), py::arg("rhovecL").noconvert(), py::arg("rhovecV").noconvert());
+ m.def("trace_VLE_isotherm_binary", &trace_VLE_isotherm_binary<Model, double, Eigen::ArrayXd>, py::arg("model"), py::arg("T"), py::arg("rhovecL0").noconvert(), py::arg("rhovecV0").noconvert(), py::arg_v("options", std::nullopt, "None"));
//cls.def("get_Ar01", [](const Model& m, const double T, const Eigen::ArrayXd& rhovec) { return id::get_Ar01(m, T, rhovec); });
//cls.def("get_Ar10", [](const Model& m, const double T, const Eigen::ArrayXd& rhovec) { return id::get_Ar10(m, T, rhovec); });
diff --git a/src/tests/catch_test_cubics.cxx b/src/tests/catch_test_cubics.cxx
index 7edd6a1a15728eba6f64a0934a5e5eebe8cfcb1d..737f9208dae1b4f3db69cb7cea42a106ad1fd357 100644
--- a/src/tests/catch_test_cubics.cxx
+++ b/src/tests/catch_test_cubics.cxx
@@ -1,10 +1,12 @@
+
+#include <fstream>
+
#include "catch/catch.hpp"
#include "teqp/models/cubics.hpp"
#include "teqp/derivs.hpp"
#include "teqp/algorithms/VLE.hpp"
-
#include <boost/numeric/odeint/stepper/euler.hpp>
#include <boost/numeric/odeint/stepper/runge_kutta_cash_karp54.hpp>
@@ -88,32 +90,52 @@ TEST_CASE("Check manual integration of subcritical VLE isotherm for binary mixtu
auto pfromderiv = rho * model.R(molefrac) * T + id::get_pr(model, T, rhovecL);
return pfromderiv;
};
-
- for (int i : { 0 }) {
- state_type X0 = get_start(T, i); // Starting point; liquid, then vapor
- double p0 = get_p(X0);
- state_type Xfinal = get_start(T, 1-i); // Ending point; liquid, then vapor
- double pfinal = get_p(Xfinal);
- //euler<state_type> integrator;
- runge_kutta_cash_karp54< state_type > integrator;
- int Nstep = 10000;
- double p = p0, pmax = pfinal, dp = (pmax - p0) / (Nstep - 1);
+ SECTION("Manual integration") {
+
+ for (int i : { 0 }) {
+ state_type X0 = get_start(T, i); // Starting point; liquid, then vapor
+ double p0 = get_p(X0);
+ state_type Xfinal = get_start(T, 1 - i); // Ending point; liquid, then vapor
+ double pfinal = get_p(Xfinal);
+
+ //euler<state_type> integrator;
+ runge_kutta_cash_karp54< state_type > integrator;
+ int Nstep = 10000;
+ double p = p0, pmax = pfinal, dp = (pmax - p0) / (Nstep - 1);
- auto write = [&]() {
- //std::cout << p << " " << X0[0] << "," << X0[1] << std::endl;
- };
- for (auto i = 0; p < pmax; ++i) {
- if (p + dp > pmax) { break; }
+ auto write = [&]() {
+ //std::cout << p << " " << X0[0] << "," << X0[1] << std::endl;
+ };
+ for (auto i = 0; p < pmax; ++i) {
+ if (p + dp > pmax) { break; }
+ write();
+ integrator.do_step(xprime, X0, p, dp);
+ p += dp;
+ }
+ double diffs = 0;
+ for (auto i = 0; i < X0.size(); ++i) {
+ diffs += std::abs(X0[i] - Xfinal[i]);
+ }
+ CHECK(diffs < 0.1);
write();
- integrator.do_step(xprime, X0, p, dp);
- p += dp;
}
- double diffs = 0;
- for (auto i = 0; i < X0.size(); ++i) {
- diffs += std::abs(X0[i] - Xfinal[i]);
- }
- CHECK(diffs < 0.1);
- write();
+ }
+ SECTION("Parametric integration of isotherm") {
+ int i = 0;
+ auto X = get_start(T, 0);
+ state_type Xfinal = get_start(T, 1 - i); // Ending point; liquid, then vapor
+ double pfinal_goal = get_p(Xfinal);
+
+ auto N = X.size() / 2;
+ Eigen::ArrayXd rhovecL0 = Eigen::Map<const Eigen::ArrayXd>(&(X[0]), N);
+ Eigen::ArrayXd rhovecV0 = Eigen::Map<const Eigen::ArrayXd>(&(X[0]) + N, N);
+ auto J = trace_VLE_isotherm_binary(model, T, rhovecL0, rhovecV0);
+ auto Nstep = J.size();
+
+ std::ofstream file("isoT.json"); file << J;
+
+ double pfinal = J.back().at("pL / Pa").back();
+ CHECK(std::abs(pfinal / pfinal_goal-1) < 1e-6);
}
}
\ No newline at end of file