diff --git a/include/teqp/algorithms/VLE.hpp b/include/teqp/algorithms/VLE.hpp
index 3a17d69c6e90fc2ba4f71db5ddabd146c651442f..141347122c1f245be1e110125f15fde498108031 100644
--- a/include/teqp/algorithms/VLE.hpp
+++ b/include/teqp/algorithms/VLE.hpp
@@ -63,18 +63,18 @@ public:
using EigenMatrix = Eigen::Array<TYPE, 2, 2>;
private:
const Model& m_model;
- TYPE m_T;
+ const TYPE m_T;
+ const Eigen::ArrayXd molefracs;
EigenMatrix J;
EigenArray y;
-
public:
std::size_t icall = 0;
double Rr, R0;
- IsothermPureVLEResiduals(const Model& model, TYPE T) : m_model(model), m_T(T) {
- std::valarray<double> molefrac = { 1.0 };
- Rr = m_model.R(molefrac);
- R0 = m_model.R(molefrac);
+ IsothermPureVLEResiduals(const Model& model, const TYPE& T, const std::optional<Eigen::ArrayXd>& molefracs_ = std::nullopt) : m_model(model), m_T(T),
+ molefracs( (molefracs_) ? molefracs_.value() : Eigen::ArrayXd::Ones(1,1)) {
+ Rr = m_model.R(molefracs);
+ R0 = m_model.R(molefracs);
};
const auto& get_errors() { return y; };
@@ -85,9 +85,8 @@ public:
const EigenArray1 rhovecL = rhovec.head(1);
const EigenArray1 rhovecV = rhovec.tail(1);
const auto rhomolarL = rhovecL.sum(), rhomolarV = rhovecV.sum();
- const auto molefracs = (EigenArray1() << 1.0).finished();
- using tdx = TDXDerivatives<Model,TYPE,EigenArray1>;
+ using tdx = TDXDerivatives<Model,TYPE,Eigen::ArrayXd>;
const TYPE &T = m_T;
//const TYPE R = m_model.R(molefracs);
@@ -163,11 +162,48 @@ auto do_pure_VLE_T(Residual &resid, Scalar rhoL, Scalar rhoV, int maxiter) {
}
template<typename Model, typename Scalar, ADBackends backend = ADBackends::autodiff>
-auto pure_VLE_T(const Model& model, Scalar T, Scalar rhoL, Scalar rhoV, int maxiter) {
- auto res = IsothermPureVLEResiduals<Model, Scalar, backend>(model, T);
+auto pure_VLE_T(const Model& model, Scalar T, Scalar rhoL, Scalar rhoV, int maxiter, const std::optional<Eigen::ArrayXd>& molefracs = std::nullopt) {
+ Eigen::ArrayXd molefracs_{Eigen::ArrayXd::Ones(1,1)};
+ if (molefracs){ molefracs_ = molefracs.value(); }
+ auto res = IsothermPureVLEResiduals<Model, Scalar, backend>(model, T, molefracs_);
return do_pure_VLE_T(res, rhoL, rhoV, maxiter);
}
+template<typename Model>
+auto pure_trace_VLE(const Model& model, const double T, const nlohmann::json &spec){
+ // Start at the true critical point, from the specified guess value
+ nlohmann::json pure_spec;
+ Eigen::ArrayXd z{Eigen::ArrayXd::Ones(1,1)};
+ if (spec.contains("pure_spec")){
+ pure_spec = spec.at("pure_spec");
+ z = Eigen::ArrayXd(pure_spec.at("alternative_length").get<int>()); z.setZero();
+ z(pure_spec.at("alternative_pure_index").get<int>()) = 1;
+ }
+
+ auto [Tc, rhoc] = solve_pure_critical(model, spec.at("Tcguess").get<double>(), spec.at("rhocguess").get<double>(), pure_spec);
+
+ // Small step towards lower temperature close to critical temperature
+ double Tclose = spec.at("Tred").get<double>()*Tc;
+ auto rhoLrhoV = extrapolate_from_critical(model, Tc, rhoc, Tclose, z);
+ auto rhoLrhoVpolished = pure_VLE_T(model, Tclose, rhoLrhoV[0], rhoLrhoV[1], spec.value("NVLE", 10), z);
+ if (rhoLrhoVpolished[0] == rhoLrhoVpolished[1]){
+ throw teqp::IterationError("Converged to trivial solution");
+ }
+
+ // Integrate down to temperature of interest
+ int Nstep = spec.at("Nstep");
+ for (auto T_: Eigen::ArrayXd::LinSpaced(Nstep, Tclose, T)){
+ rhoLrhoVpolished = pure_VLE_T(model, T_, rhoLrhoVpolished[0], rhoLrhoVpolished[1], spec.value("NVLE", 10), z);
+ if (!std::isfinite(rhoLrhoVpolished[0])){
+ throw teqp::IterationError("The density is no longer valid; try increasing Nstep");
+ }
+ if (rhoLrhoVpolished[0] == rhoLrhoVpolished[1]){
+ throw teqp::IterationError("Converged to trivial solution; try increasing Nstep");
+ }
+ }
+ return rhoLrhoVpolished;
+}
+
/***
* \brief Calculate the derivative of vapor pressure with respect to temperature
* \param model The model to operate on
@@ -182,9 +218,10 @@ auto pure_VLE_T(const Model& model, Scalar T, Scalar rhoL, Scalar rhoV, int maxi
* where the \f$h''-h'\f$ is given by the difference in residual enthalpy \f$h''-h' = h^r''-h^r'\f$ because the ideal-gas parts cancel
*/
template<typename Model, typename Scalar, ADBackends backend = ADBackends::autodiff>
-auto dpsatdT_pure(const Model& model, Scalar T, Scalar rhoL, Scalar rhoV) {
+auto dpsatdT_pure(const Model& model, Scalar T, Scalar rhoL, Scalar rhoV, const std::optional<Eigen::ArrayXd>& molefracs = std::nullopt) {
auto molefrac = (Eigen::ArrayXd(1) << 1.0).finished();
+ if (molefracs){ molefrac = molefracs.value(); }
using tdx = TDXDerivatives<decltype(model), double, decltype(molefrac)>;
using iso = IsochoricDerivatives<decltype(model), double, decltype(molefrac)>;
@@ -939,7 +976,8 @@ auto get_dpsat_dTsat_isopleth(const Model& model, const Scalar& T, const VecType
}
/***
-* \brief Trace an isotherm with parametric tracing
+ * \brief Trace an isotherm with parametric tracing
+ * \ note If options.revision is 2, the data will be returned in the "data" field, otherwise the data will be returned as root array
*/
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)
@@ -1167,7 +1205,22 @@ auto trace_VLE_isotherm_binary(const Model &model, Scalar T, VecType rhovecL0, V
store_point(); // last_drhodt is updated;
}
- return JSONdata;
+ if (opt.revision == 1){
+ return JSONdata;
+ }
+ else if (opt.revision == 2){
+ nlohmann::json meta{
+ {"termination_reason", termination_reason}
+ };
+ return nlohmann::json{
+ {"meta", meta},
+ {"data", JSONdata}
+ };
+ }
+ else
+ {
+ throw teqp::InvalidArgument("revision is not valid");
+ }
}
/***
diff --git a/include/teqp/algorithms/VLE_types.hpp b/include/teqp/algorithms/VLE_types.hpp
index 90bf24026e9db526734db47638f5c2e6afec2b05..d9dd1666f165c66bb71655e8a0b3102a2805d6fa 100644
--- a/include/teqp/algorithms/VLE_types.hpp
+++ b/include/teqp/algorithms/VLE_types.hpp
@@ -4,7 +4,7 @@ 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;
+ int max_steps = 1000, integration_order = 5, revision = 1;
bool polish = true;
bool calc_criticality = false;
bool terminate_unstable = false;
diff --git a/include/teqp/algorithms/critical_pure.hpp b/include/teqp/algorithms/critical_pure.hpp
index df6afb6876329e5dbe7bb213e1f81f17c86b1798..04f9055fc1f8ec0712d5842d6bfae4aaf383eeef 100644
--- a/include/teqp/algorithms/critical_pure.hpp
+++ b/include/teqp/algorithms/critical_pure.hpp
@@ -77,20 +77,25 @@ namespace teqp {
}
template<typename Model, typename Scalar, ADBackends backend = ADBackends::autodiff>
- auto solve_pure_critical(const Model& model, const Scalar T0, const Scalar rho0, const nlohmann::json& flags = {}) {
+ auto solve_pure_critical(const Model& model, const Scalar T0, const Scalar rho0, const std::optional<nlohmann::json>& flags = std::nullopt) {
auto x = (Eigen::ArrayXd(2) << T0, rho0).finished();
- int maxsteps = (flags.contains("maxsteps")) ? flags.at("maxsteps").get<int>() : 10;
+ int maxsteps = 10;
std::optional<std::size_t> alternative_pure_index;
- if (flags.contains("alternative_pure_index")){
- auto i = flags.at("alternative_pure_index").get<int>();
- if (i < 0){ throw teqp::InvalidArgument("alternative_pure_index cannot be less than 0"); }
- alternative_pure_index = i;
- }
std::optional<std::size_t> alternative_length;
- if (flags.contains("alternative_length")){
- auto i = flags.at("alternative_length").get<int>();
- if (i < 2){ throw teqp::InvalidArgument("alternative_length cannot be less than 2"); }
- alternative_length = i;
+ if (flags){
+ if (flags.value().contains("maxsteps")){
+ maxsteps = flags.value().at("maxsteps");
+ }
+ if (flags.value().contains("alternative_pure_index")){
+ auto i = flags.value().at("alternative_pure_index").get<int>();
+ if (i < 0){ throw teqp::InvalidArgument("alternative_pure_index cannot be less than 0"); }
+ alternative_pure_index = i;
+ }
+ if (flags.value().contains("alternative_length")){
+ auto i = flags.value().at("alternative_length").get<int>();
+ if (i < 2){ throw teqp::InvalidArgument("alternative_length cannot be less than 2"); }
+ alternative_length = i;
+ }
}
// A convenience method to make linear system solving more concise with Eigen datatypes
// All arguments are converted to matrices, the solve is done, and an array is returned
@@ -106,19 +111,29 @@ namespace teqp {
}
template<typename Model, typename Scalar>
- Eigen::ArrayXd extrapolate_from_critical(const Model& model, const Scalar& Tc, const Scalar& rhoc, const Scalar& T) {
+ Scalar get_Brho_critical_extrap(const Model& model, const Scalar& Tc, const Scalar& rhoc, const std::optional<Eigen::ArrayXd>& z = std::nullopt) {
using tdx = TDXDerivatives<Model, Scalar>;
- auto z = (Eigen::ArrayXd(1) << 1.0).finished();
- auto R = model.R(z);
- auto ders = tdx::template get_Ar0n<4>(model, Tc, rhoc, z);
+ auto z_ = (Eigen::ArrayXd(1) << 1.0).finished();
+ if (z){
+ z_ = z.value();
+ }
+ auto R = model.R(z_);
+ auto ders = tdx::template get_Ar0n<4>(model, Tc, rhoc, z_);
//auto dpdrho = R*Tc*(1 + 2 * ders[1] + ders[2]); // Should be zero
//auto d2pdrho2 = R*Tc/rhoc*(2 * ders[1] + 4 * ders[2] + ders[3]); // Should be zero
auto d3pdrho3 = R * Tc / (rhoc * rhoc) * (6 * ders[2] + 6 * ders[3] + ders[4]);
- auto Ar11 = tdx::template get_Ar11(model, Tc, rhoc, z);
- auto Ar12 = tdx::template get_Ar12(model, Tc, rhoc, z);
+ auto Ar11 = tdx::template get_Ar11(model, Tc, rhoc, z_);
+ auto Ar12 = tdx::template get_Ar12(model, Tc, rhoc, z_);
auto d2pdrhodT = R * (1 + 2 * ders[1] + ders[2] - 2 * Ar11 - Ar12);
auto Brho = sqrt(6 * d2pdrhodT * Tc / d3pdrho3);
+ return Brho;
+ }
+
+ template<typename Model, typename Scalar>
+Eigen::Array<double, 2, 1> extrapolate_from_critical(const Model& model, const Scalar& Tc, const Scalar& rhoc, const Scalar& T, const std::optional<Eigen::ArrayXd>& z = std::nullopt) {
+
+ auto Brho = get_Brho_critical_extrap(model, Tc, rhoc, z);
auto drhohat_dT = Brho / Tc;
auto dT = T - Tc;
diff --git a/include/teqp/exceptions.hpp b/include/teqp/exceptions.hpp
index a17fc2ac7c2addeecb092678e0ababa3d60db1ee..90080a39376a2a7a5a5ed12cb969225fd108ded2 100644
--- a/include/teqp/exceptions.hpp
+++ b/include/teqp/exceptions.hpp
@@ -36,5 +36,6 @@ namespace teqp {
public:
IterationFailure(const std::string& msg) : teqpException(100, msg) {};
};
+ using IterationError = IterationFailure;
-}; // namespace teqp
\ No newline at end of file
+}; // namespace teqp