From 7ee5a32db41e7101d29f1d95f2cdbf2a2b57f0ec Mon Sep 17 00:00:00 2001
From: Ian Bell <ian.bell@nist.gov>
Date: Fri, 4 Nov 2022 11:03:40 -0400
Subject: [PATCH] The remaining VLE routines
So inconsistent in naming...
---
include/teqp/algorithms/VLE.hpp | 57 ++++-----------------------
include/teqp/algorithms/VLE_types.hpp | 50 +++++++++++++++++++++++
include/teqp/cpp/teqpcpp.hpp | 6 +++
interface/CPP/teqpcpp.cpp | 26 ++++++++++++
interface/pybind11_wrapper.cpp | 11 +++---
5 files changed, 96 insertions(+), 54 deletions(-)
create mode 100644 include/teqp/algorithms/VLE_types.hpp
diff --git a/include/teqp/algorithms/VLE.hpp b/include/teqp/algorithms/VLE.hpp
index 530939b..ddf84e6 100644
--- a/include/teqp/algorithms/VLE.hpp
+++ b/include/teqp/algorithms/VLE.hpp
@@ -5,6 +5,7 @@
#include "teqp/exceptions.hpp"
#include "teqp/algorithms/critical_tracing.hpp"
#include "teqp/algorithms/critical_pure.hpp"
+#include "teqp/algorithms/VLE_types.hpp"
#include <Eigen/Dense>
// Imports from boost for numerical integration
@@ -167,8 +168,6 @@ auto pure_VLE_T(const Model& model, Scalar T, Scalar rhoL, Scalar rhoV, int maxi
return do_pure_VLE_T(res, rhoL, rhoV, maxiter);
}
-enum class VLE_return_code { unset, xtol_satisfied, functol_satisfied, maxfev_met, maxiter_met, notfinite_step };
-
/***
* \brief Do a vapor-liquid phase equilibrium problem for a mixture (binary only for now) with mole fractions specified in the liquid phase
* \param model The model to operate on
@@ -306,25 +305,6 @@ auto mix_VLE_Tx(const Model& model, Scalar T, const Vector& rhovecL0, const Vect
return std::make_tuple(return_code, rhovecLfinal, rhovecVfinal);
}
-struct MixVLETpFlags {
- double atol = 1e-10,
- reltol = 1e-10,
- axtol = 1e-10,
- relxtol = 1e-10,
- relaxation = 1.0;
- int maxiter = 10;
-};
-
-struct MixVLEReturn {
- bool success = false;
- std::string message = "";
- Eigen::ArrayXd rhovecL, rhovecV;
- VLE_return_code return_code;
- int num_iter=-1, num_fev=-1;
- double T=-1;
- Eigen::ArrayXd r, initial_r;
-};
-
template<typename Model>
struct hybrj_functor__mix_VLE_Tp : Functor<double>
{
@@ -419,7 +399,9 @@ struct hybrj_functor__mix_VLE_Tp : Functor<double>
* \param flags Flags controlling the iteration and stopping conditions
*/
template<typename Model, typename Scalar, typename Vector>
-auto mix_VLE_Tp(const Model& model, Scalar T, Scalar pgiven, const Vector& rhovecL0, const Vector& rhovecV0, const MixVLETpFlags& flags = {}) {
+auto mix_VLE_Tp(const Model& model, Scalar T, Scalar pgiven, const Vector& rhovecL0, const Vector& rhovecV0, const std::optional<MixVLETpFlags>& flags_ = std::nullopt) {
+
+ auto flags = flags_.value_or(MixVLETpFlags{});
const Eigen::Index N = rhovecL0.size();
auto lengths = (Eigen::ArrayXi(2) << rhovecL0.size(), rhovecV0.size()).finished();
@@ -518,14 +500,6 @@ auto mix_VLE_Tp(const Model& model, Scalar T, Scalar pgiven, const Vector& rhove
return r;
}
-struct MixVLEpxFlags {
- double atol = 1e-10,
- reltol = 1e-10,
- axtol = 1e-10,
- relxtol = 1e-10;
- int maxiter = 10;
-};
-
/***
* \brief Do vapor-liquid phase equilibrium problem at specified pressure and mole fractions in the bulk phase
* \param model The model to operate on
@@ -538,7 +512,9 @@ struct MixVLEpxFlags {
* \param flags Additional flags
*/
template<typename Model, typename Scalar, typename Vector>
-auto mixture_VLE_px(const Model& model, Scalar p_spec, const Vector& xmolar_spec, Scalar T0, const Vector& rhovecL0, const Vector& rhovecV0, const MixVLEpxFlags& flags = {}) {
+auto mixture_VLE_px(const Model& model, Scalar p_spec, const Vector& xmolar_spec, Scalar T0, const Vector& rhovecL0, const Vector& rhovecV0, const std::optional<MixVLEpxFlags>& flags_ = std::nullopt) {
+
+ auto flags = flags_.value_or(MixVLEpxFlags{});
const Eigen::Index N = rhovecL0.size();
auto lengths = (Eigen::ArrayXi(3) << rhovecL0.size(), rhovecV0.size(), xmolar_spec.size()).finished();
@@ -932,14 +908,6 @@ auto get_dpsat_dTsat_isopleth(const Model& model, const Scalar& T, const VecType
//return der;
}
-struct TVLEOptions {
- 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;
- bool calc_criticality = false;
- bool terminate_unstable = false;
-};
-
/***
* \brief Trace an isotherm with parametric tracing
*/
@@ -1166,15 +1134,6 @@ auto trace_VLE_isotherm_binary(const Model &model, Scalar T, VecType rhovecL0, V
return JSONdata;
}
-
-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;
- bool calc_criticality = false;
- bool terminate_unstable = false;
-};
-
/***
* \brief Trace an isobar with parametric tracing
*/
@@ -1408,4 +1367,4 @@ auto trace_VLE_isobar_binary(const Model& model, Scalar p, Scalar T0, VecType rh
return JSONdata;
}
-}; /* namespace teqp*/
\ No newline at end of file
+}; /* namespace teqp*/
diff --git a/include/teqp/algorithms/VLE_types.hpp b/include/teqp/algorithms/VLE_types.hpp
new file mode 100644
index 0000000..8eb0964
--- /dev/null
+++ b/include/teqp/algorithms/VLE_types.hpp
@@ -0,0 +1,50 @@
+#pragma once
+
+namespace teqp{
+
+struct TVLEOptions {
+ 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;
+ 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;
+ bool calc_criticality = false;
+ bool terminate_unstable = false;
+};
+
+struct MixVLEpxFlags {
+ double atol = 1e-10,
+ reltol = 1e-10,
+ axtol = 1e-10,
+ relxtol = 1e-10;
+ int maxiter = 10;
+};
+
+struct MixVLETpFlags {
+ double atol = 1e-10,
+ reltol = 1e-10,
+ axtol = 1e-10,
+ relxtol = 1e-10,
+ relaxation = 1.0;
+ int maxiter = 10;
+};
+
+enum class VLE_return_code { unset, xtol_satisfied, functol_satisfied, maxfev_met, maxiter_met, notfinite_step };
+
+struct MixVLEReturn {
+ bool success = false;
+ std::string message = "";
+ Eigen::ArrayXd rhovecL, rhovecV;
+ VLE_return_code return_code;
+ int num_iter=-1, num_fev=-1;
+ double T=-1;
+ Eigen::ArrayXd r, initial_r;
+};
+
+}
diff --git a/include/teqp/cpp/teqpcpp.hpp b/include/teqp/cpp/teqpcpp.hpp
index 8c3d435..07b7850 100644
--- a/include/teqp/cpp/teqpcpp.hpp
+++ b/include/teqp/cpp/teqpcpp.hpp
@@ -8,6 +8,7 @@
// The only headers that can be included here are
// ones that define and use POD (plain ole' data) types
#include "teqp/algorithms/critical_tracing_types.hpp"
+#include "teqp/algorithms/VLE_types.hpp"
using EArray2 = Eigen::Array<double, 2, 1>;
using EArrayd = Eigen::ArrayX<double>;
@@ -126,6 +127,11 @@ namespace teqp {
virtual std::tuple<EArrayd, EArrayd> get_drhovecdp_Tsat(const double T, const REArrayd& rhovecL, const REArrayd& rhovecV) const = 0;
virtual std::tuple<EArrayd, EArrayd> get_drhovecdT_psat(const double T, const REArrayd& rhovecL, const REArrayd& rhovecV) const = 0;
virtual double get_dpsat_dTsat_isopleth(const double T, const REArrayd& rhovecL, const REArrayd& rhovecV) const = 0;
+ virtual nlohmann::json trace_VLE_isotherm_binary(const double T0, const EArrayd& rhovec0, const EArrayd& rhovecV0, const std::optional<TVLEOptions> &) const = 0;
+ virtual nlohmann::json trace_VLE_isobar_binary(const double p, const double T0, const EArrayd& rhovecL0, const EArrayd& rhovecV0, const std::optional<PVLEOptions> &) const = 0;
+ virtual std::tuple<VLE_return_code,EArrayd,EArrayd> mix_VLE_Tx(const double T, const REArrayd& rhovecL0, const REArrayd& rhovecV0, const REArrayd& xspec, const double atol, const double reltol, const double axtol, const double relxtol, const int maxiter) const = 0;
+ virtual MixVLEReturn mix_VLE_Tp(const double T, const double pgiven, const REArrayd& rhovecL0, const REArrayd& rhovecV0, const std::optional<MixVLETpFlags> &flags) const = 0;
+ virtual std::tuple<VLE_return_code,double,EArrayd,EArrayd> mixture_VLE_px(const double p_spec, const REArrayd& xmolar_spec, const double T0, const REArrayd& rhovecL0, const REArrayd& rhovecV0, const std::optional<MixVLEpxFlags>& flags) const = 0;
virtual nlohmann::json trace_critical_arclength_binary(const double T0, const EArrayd& rhovec0, const std::optional<std::string>&, const std::optional<TCABOptions> &) const = 0;
virtual EArrayd get_drhovec_dT_crit(const double T, const REArrayd& rhovec) const = 0;
diff --git a/interface/CPP/teqpcpp.cpp b/interface/CPP/teqpcpp.cpp
index 9be1aba..eb4293b 100644
--- a/interface/CPP/teqpcpp.cpp
+++ b/interface/CPP/teqpcpp.cpp
@@ -213,6 +213,32 @@ namespace teqp {
return teqp::get_dpsat_dTsat_isopleth(model, T, rhovecL, rhovecV);
}, m_model);
}
+
+ nlohmann::json trace_VLE_isotherm_binary(const double T0, const EArrayd& rhovecL0, const EArrayd& rhovecV0, const std::optional<TVLEOptions> &options) const override{
+ return std::visit([&](const auto& model) {
+ return teqp::trace_VLE_isotherm_binary(model, T0, rhovecL0, rhovecV0, options);
+ }, m_model);
+ }
+ nlohmann::json trace_VLE_isobar_binary(const double p, const double T0, const EArrayd& rhovecL0, const EArrayd& rhovecV0, const std::optional<PVLEOptions> &options) const override{
+ return std::visit([&](const auto& model) {
+ return teqp::trace_VLE_isobar_binary(model, p, T0, rhovecL0, rhovecV0, options);
+ }, m_model);
+ }
+ std::tuple<VLE_return_code,EArrayd,EArrayd> mix_VLE_Tx(const double T, const REArrayd& rhovecL0, const REArrayd& rhovecV0, const REArrayd& xspec, const double atol, const double reltol, const double axtol, const double relxtol, const int maxiter) const override{
+ return std::visit([&](const auto& model) {
+ return teqp::mix_VLE_Tx(model, T, rhovecL0, rhovecV0, xspec, atol, reltol, axtol, relxtol, maxiter);
+ }, m_model);
+ }
+ MixVLEReturn mix_VLE_Tp(const double T, const double pgiven, const REArrayd& rhovecL0, const REArrayd& rhovecV0, const std::optional<MixVLETpFlags> &flags) const override{
+ return std::visit([&](const auto& model) {
+ return teqp::mix_VLE_Tp(model, T, pgiven, rhovecL0, rhovecV0, flags);
+ }, m_model);
+ }
+ std::tuple<VLE_return_code,double,EArrayd,EArrayd> mixture_VLE_px(const double p_spec, const REArrayd& xmolar_spec, const double T0, const REArrayd& rhovecL0, const REArrayd& rhovecV0, const std::optional<MixVLEpxFlags>& flags ) const override{
+ return std::visit([&](const auto& model) {
+ return teqp::mixture_VLE_px(model, p_spec, xmolar_spec, T0, rhovecL0, rhovecV0, flags);
+ }, m_model);
+ }
};
std::shared_ptr<AbstractModel> make_model(const nlohmann::json& j) {
diff --git a/interface/pybind11_wrapper.cpp b/interface/pybind11_wrapper.cpp
index 017cc23..41a27c5 100644
--- a/interface/pybind11_wrapper.cpp
+++ b/interface/pybind11_wrapper.cpp
@@ -239,11 +239,12 @@ void init_teqp(py::module& m) {
.def("get_drhovecdT_psat", &am::get_drhovecdT_psat, "T"_a, "rhovecL"_a.noconvert(), "rhovecV"_a.noconvert())
.def("get_dpsat_dTsat_isopleth", &am::get_dpsat_dTsat_isopleth, "T"_a, "rhovecL"_a.noconvert(), "rhovecV"_a.noconvert())
-// .def("mix_VLE_Tx", &mix_VLE_Tx<Model, double, RAX>)
-// .def("mix_VLE_Tp", &mix_VLE_Tp<Model, double, RAX>, "T"_a, "p_given"_a, "rhovecL0"_a.noconvert(), "rhovecV0"_a.noconvert(), py::arg_v("flags", MixVLETpFlags{}, "None"))
-// .def("mixture_VLE_px", &mixture_VLE_px<Model, double, RAX>, "p_spec"_a, "xmolar_spec"_a.noconvert(), "T0"_a, "rhovecL0"_a.noconvert(), "rhovecV0"_a.noconvert(), py::arg_v("flags", MixVLEpxFlags{}, "None"))
-// .def("trace_VLE_isotherm_binary", &trace_VLE_isotherm_binary, "T"_a, "rhovecL0"_a.noconvert(), "rhovecV0"_a.noconvert(), py::arg_v("options", std::nullopt, "None"))
-// .def("trace_VLE_isobar_binary", &trace_VLE_isobar_binary, "p"_a, "T0"_a, "rhovecL0"_a.noconvert(), "rhovecV0"_a.noconvert(), py::arg_v("options", std::nullopt, "None"))
+ .def("trace_VLE_isotherm_binary", &am::trace_VLE_isotherm_binary, "T"_a, "rhovecL0"_a.noconvert(), "rhovecV0"_a.noconvert(), py::arg_v("options", std::nullopt, "None"))
+ .def("trace_VLE_isobar_binary", &am::trace_VLE_isobar_binary, "p"_a, "T0"_a, "rhovecL0"_a.noconvert(), "rhovecV0"_a.noconvert(), py::arg_v("options", std::nullopt, "None"))
+ .def("mix_VLE_Tx", &am::mix_VLE_Tx, "T"_a, "rhovecL0"_a.noconvert(), "rhovecV0"_a.noconvert(), "xspec"_a.noconvert(), "atol"_a, "reltol"_a, "axtol"_a, "relxtol"_a, "maxiter"_a)
+ .def("mix_VLE_Tp", &am::mix_VLE_Tp, "T"_a, "p_given"_a, "rhovecL0"_a.noconvert(), "rhovecV0"_a.noconvert(), py::arg_v("options", std::nullopt, "None"))
+ .def("mixture_VLE_px", &am::mixture_VLE_px, "p_spec"_a, "xmolar_spec"_a.noconvert(), "T0"_a, "rhovecL0"_a.noconvert(), "rhovecV0"_a.noconvert(), py::arg_v("options", std::nullopt, "None"))
+
// .def("mix_VLLE_T", &mix_VLLE_T<Model, double, RAX>);
// .def("find_VLLE_T_binary", &am::find_VLLE_T_binary, "traces"_a, py::arg_v("options", std::nullopt, "None"));
;
--
GitLab