diff --git a/include/teqp/algorithms/VLE.hpp b/include/teqp/algorithms/VLE.hpp index ec495583..f855d657 100644 --- a/include/teqp/algorithms/VLE.hpp +++ b/include/teqp/algorithms/VLE.hpp @@ -997,13 +997,25 @@ inline auto trace_VLE_isotherm_binary(const AbstractModel &model, double T, cons auto rhovecV = Eigen::Map(&(x0[0 + N]), N).eval(); auto x = (Eigen::ArrayXd(2) << rhovecL(0) / rhovecL.sum(), rhovecL(1) / rhovecL.sum()).finished(); // Mole fractions in the liquid phase (to be kept constant) auto [return_code, rhovecLnew, rhovecVnew] = model.mix_VLE_Tx(T, rhovecL, rhovecV, x, 1e-10, 1e-8, 1e-10, 1e-8, 10); - - // If the step is accepted, copy into x again ... - auto rhovecLview = Eigen::Map(&(x0[0]), N); - auto rhovecVview = Eigen::Map(&(x0[0]) + N, N); - rhovecLview = rhovecLnew; - rhovecVview = rhovecVnew; - //std::cout << "[polish]: " << static_cast(return_code) << ": " << rhovecLnew.sum() / rhovecL.sum() << " " << rhovecVnew.sum() / rhovecV.sum() << std::endl; + + if (((rhovecL-rhovecLnew).cwiseAbs() > opt.polish_reltol_rho*rhovecL).any()){ + std::string msg; + if (opt.polish_exception_on_fail){ + throw IterationFailure(msg); + } + else{ + if (opt.verbosity > 0){ + std::cout << msg << std::endl; + } + } + } + else{ + // If the step is accepted, copy into x again ... + auto rhovecLview = Eigen::Map(&(x0[0]), N); + auto rhovecVview = Eigen::Map(&(x0[0]) + N, N); + rhovecLview = rhovecLnew; + rhovecVview = rhovecVnew; + } } std::swap(previous_drhodt, last_drhodt); @@ -1220,7 +1232,7 @@ auto trace_VLE_isobar_binary(const Model& model, double p, double T0, const Eige if (opt.calc_criticality) { auto condsL = model.get_criticality_conditions(T, rhovecL); auto condsV = model.get_criticality_conditions(T, rhovecV); - if (condsL[0] < 1e-12 || condsV[0] < 1e-12) { + if (condsL[0] < opt.crit_termination || condsV[0] < opt.crit_termination) { return true; } } @@ -1242,13 +1254,26 @@ auto trace_VLE_isobar_binary(const Model& model, double p, double T0, const Eige auto x = (Eigen::ArrayXd(2) << rhovecL(0) / rhovecL.sum(), rhovecL(1) / rhovecL.sum()).finished(); // Mole fractions in the liquid phase (to be kept constant) auto [return_code, Tnew, rhovecLnew, rhovecVnew] = model.mixture_VLE_px(p, x, T, rhovecL, rhovecV); - // If the step is accepted, copy into x again ... - x0[0] = Tnew; - auto rhovecLview = Eigen::Map(&(x0[1]), N); - auto rhovecVview = Eigen::Map(&(x0[1]) + N, N); - rhovecLview = rhovecLnew; - rhovecVview = rhovecVnew; - //std::cout << "[polish]: " << static_cast(return_code) << ": " << rhovecLnew.sum() / rhovecL.sum() << " " << rhovecVnew.sum() / rhovecV.sum() << std::endl; + if (((rhovecL-rhovecLnew).cwiseAbs() > opt.polish_reltol_rho*rhovecL).any()){ + std::string msg; + if (opt.polish_exception_on_fail){ + throw IterationFailure(msg); + } + else{ + if (opt.verbosity > 0){ + std::cout << msg << std::endl; + } + } + } + else{ + // If the step is accepted, copy into x again ... + x0[0] = Tnew; + auto rhovecLview = Eigen::Map(&(x0[1]), N); + auto rhovecVview = Eigen::Map(&(x0[1]) + N, N); + rhovecLview = rhovecLnew; + rhovecVview = rhovecVnew; + //std::cout << "[polish]: " << static_cast(return_code) << ": " << rhovecLnew.sum() / rhovecL.sum() << " " << rhovecVnew.sum() / rhovecV.sum() << std::endl; + } } std::swap(previous_drhodt, last_drhodt); diff --git a/include/teqp/algorithms/VLE_types.hpp b/include/teqp/algorithms/VLE_types.hpp index d9dd1666..aba0f468 100644 --- a/include/teqp/algorithms/VLE_types.hpp +++ b/include/teqp/algorithms/VLE_types.hpp @@ -4,16 +4,18 @@ namespace teqp{ struct TVLEOptions { double init_dt = 1e-5, abs_err = 1e-8, rel_err = 1e-8, max_dt = 100000, init_c = 1.0, p_termination = 1e15, crit_termination = 1e-12; - int max_steps = 1000, integration_order = 5, revision = 1; - bool polish = true; + int max_steps = 1000, integration_order = 5, revision = 1, verbosity = 0; + bool polish = true, polish_exception_on_fail = false; + double polish_reltol_rho = 0.05; bool calc_criticality = false; bool terminate_unstable = false; }; struct PVLEOptions { - double init_dt = 1e-5, abs_err = 1e-8, rel_err = 1e-8, max_dt = 100000, init_c = 1.0; - int max_steps = 1000, integration_order = 5; - bool polish = true; + double init_dt = 1e-5, abs_err = 1e-8, rel_err = 1e-8, max_dt = 100000, init_c = 1.0, crit_termination = 1e-12; + int max_steps = 1000, integration_order = 5, verbosity = 0; + bool polish = true, polish_exception_on_fail = false; + double polish_reltol_rho = 0.05; bool calc_criticality = false; bool terminate_unstable = false; }; diff --git a/interface/pybind11_wrapper.cpp b/interface/pybind11_wrapper.cpp index 3ada3225..e15d4f93 100644 --- a/interface/pybind11_wrapper.cpp +++ b/interface/pybind11_wrapper.cpp @@ -281,6 +281,9 @@ void init_teqp(py::module& m) { .def_readwrite("max_steps", &TVLEOptions::max_steps) .def_readwrite("integration_order", &TVLEOptions::integration_order) .def_readwrite("polish", &TVLEOptions::polish) + .def_readwrite("polish_reltol_rho", &TVLEOptions::polish_reltol_rho) + .def_readwrite("polish_exception_on_fail", &TVLEOptions::polish_exception_on_fail) + .def_readwrite("verbosity", &TVLEOptions::verbosity) .def_readwrite("calc_criticality", &TVLEOptions::calc_criticality) .def_readwrite("terminate_unstable", &TVLEOptions::terminate_unstable) ; @@ -293,9 +296,13 @@ void init_teqp(py::module& m) { .def_readwrite("init_dt", &PVLEOptions::init_dt) .def_readwrite("init_c", &PVLEOptions::init_c) .def_readwrite("max_dt", &PVLEOptions::max_dt) + .def_readwrite("crit_termination", &PVLEOptions::crit_termination) .def_readwrite("max_steps", &PVLEOptions::max_steps) .def_readwrite("integration_order", &PVLEOptions::integration_order) .def_readwrite("polish", &PVLEOptions::polish) + .def_readwrite("polish_reltol_rho", &PVLEOptions::polish_reltol_rho) + .def_readwrite("polish_exception_on_fail", &PVLEOptions::polish_exception_on_fail) + .def_readwrite("verbosity", &PVLEOptions::verbosity) .def_readwrite("calc_criticality", &PVLEOptions::calc_criticality) .def_readwrite("terminate_unstable", &PVLEOptions::terminate_unstable) ;