Skip to content

Commit

Permalink
Add more knobs for fine tuning termination around critical point
Browse files Browse the repository at this point in the history 
Can still be improved, but already a lot better
  • Loading branch information
@ianhbell
ianhbell committed Oct 20, 2023
1 parent ff03495 commit 561c7ec
Show file tree
Hide file tree
Showing 3 changed files with 54 additions and 20 deletions.
55 changes: 40 additions & 15 deletions include/teqp/algorithms/VLE.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -997,13 +997,25 @@ inline auto trace_VLE_isotherm_binary(const AbstractModel &model, double T, cons
auto rhovecV = Eigen::Map<const Eigen::ArrayXd>(&(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<Eigen::ArrayXd>(&(x0[0]), N);
auto rhovecVview = Eigen::Map<Eigen::ArrayXd>(&(x0[0]) + N, N);
rhovecLview = rhovecLnew;
rhovecVview = rhovecVnew;
//std::cout << "[polish]: " << static_cast<int>(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<Eigen::ArrayXd>(&(x0[0]), N);
auto rhovecVview = Eigen::Map<Eigen::ArrayXd>(&(x0[0]) + N, N);
rhovecLview = rhovecLnew;
rhovecVview = rhovecVnew;
}
}

std::swap(previous_drhodt, last_drhodt);
Expand Down Expand Up @@ -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;
}
}
Expand All @@ -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<Eigen::ArrayXd>(&(x0[1]), N);
auto rhovecVview = Eigen::Map<Eigen::ArrayXd>(&(x0[1]) + N, N);
rhovecLview = rhovecLnew;
rhovecVview = rhovecVnew;
//std::cout << "[polish]: " << static_cast<int>(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<Eigen::ArrayXd>(&(x0[1]), N);
auto rhovecVview = Eigen::Map<Eigen::ArrayXd>(&(x0[1]) + N, N);
rhovecLview = rhovecLnew;
rhovecVview = rhovecVnew;
//std::cout << "[polish]: " << static_cast<int>(return_code) << ": " << rhovecLnew.sum() / rhovecL.sum() << " " << rhovecVnew.sum() / rhovecV.sum() << std::endl;
}
}

std::swap(previous_drhodt, last_drhodt);
Expand Down
12 changes: 7 additions & 5 deletions include/teqp/algorithms/VLE_types.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
};
Expand Down
7 changes: 7 additions & 0 deletions interface/pybind11_wrapper.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
;
Expand All @@ -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)
;
Expand Down

0 comments on commit 561c7ec

Please sign in to comment.