diff --git a/include/teqp/algorithms/critical_tracing.hpp b/include/teqp/algorithms/critical_tracing.hpp
index 72c4e330dec18986c84a173884c8053c73fba9eb..bae1a547962e7f264f0b9d8e85014238c1514a57 100644
--- a/include/teqp/algorithms/critical_tracing.hpp
+++ b/include/teqp/algorithms/critical_tracing.hpp
@@ -36,12 +36,7 @@ struct CriticalTracing {
return std::make_tuple(es.eigenvalues(), es.eigenvectors());
}
- struct EigenData {
- Eigen::ArrayXd v0, v1, eigenvalues;
- Eigen::MatrixXd eigenvectorscols;
- };
-
- static auto eigen_problem(const Model& model, const Scalar T, const VecType& rhovec, const VecType& alignment_v0 = {}) {
+ static auto eigen_problem(const Model& model, const Scalar T, const VecType& rhovec, const std::optional<VecType>& alignment_v0 = std::nullopt) {
EigenData ed;
@@ -114,7 +109,7 @@ struct CriticalTracing {
else {
throw std::invalid_argument("More than one non-zero concentration value found; not currently supported");
}
- if (alignment_v0.size() > 0 && ed.eigenvectorscols.col(0).matrix().dot(alignment_v0.matrix()) < 0) {
+ if (alignment_v0 && ed.eigenvectorscols.col(0).matrix().dot(alignment_v0.value().matrix()) < 0) {
ed.eigenvectorscols.col(0) *= -1;
}
@@ -132,7 +127,7 @@ struct CriticalTracing {
return eigen_problem(model, T, rhovec).eigenvalues[0];
}
- static auto get_derivs(const Model& model, const Scalar T, const VecType& rhovec, const VecType& alignment_v0 = {}) {
+ static auto get_derivs(const Model& model, const Scalar T, const VecType& rhovec, const std::optional<VecType>& alignment_v0 = std::nullopt) {
auto molefrac = rhovec / rhovec.sum();
auto R = model.R(molefrac);
@@ -163,7 +158,7 @@ struct CriticalTracing {
return eval(model.alphar(T, rhotot, molefrac) * model.R(molefrac) * T * rhotot);
};
auto psir_derivs_ = derivatives(wrapper, wrt(varsigma), at(varsigma));
- VecType psir_derivs; psir_derivs.resize(5);
+ Eigen::ArrayXd psir_derivs; psir_derivs.resize(5);
for (auto i = 0; i < 5; ++i) { psir_derivs[i] = psir_derivs_[i]; }
#else
@@ -179,7 +174,7 @@ struct CriticalTracing {
return model.alphar(T, rhotot, molefrac) * model.R(molefrac) * T * rhotot;
};
auto psir_derivs_ = diff_mcx1(wrapper, 0.0, 4, true);
- VecType psir_derivs; psir_derivs.resize(5);
+ Eigen::ArrayXd psir_derivs; psir_derivs.resize(5);
for (auto i = 0; i < 5; ++i) { psir_derivs[i] = psir_derivs_[i]; }
#endif
@@ -211,7 +206,7 @@ struct CriticalTracing {
// The derivatives of total Psi w.r.t.sigma_1 (numerical for residual, analytic for ideal)
// Returns a tuple, with residual, ideal, total dicts with of number of derivatives, value of derivative
- auto all_derivs = get_derivs(model, T, rhovec);
+ auto all_derivs = get_derivs(model, T, rhovec, Eigen::ArrayXd());
auto derivs = all_derivs.tot;
// The temperature derivative of total Psi w.r.t.T from a centered finite difference in T
@@ -284,7 +279,7 @@ struct CriticalTracing {
}
static auto get_criticality_conditions(const Model& model, const Scalar T, const VecType& rhovec) {
- auto derivs = get_derivs(model, T, rhovec);
+ auto derivs = get_derivs(model, T, rhovec, Eigen::ArrayXd());
return (Eigen::ArrayXd(2) << derivs.tot[2], derivs.tot[3]).finished();
}
@@ -341,7 +336,7 @@ struct CriticalTracing {
auto T = x[0];
Eigen::ArrayXd rhovec(2); rhovec << z0*x[1], (1-z0)*x[1];
//auto z0new = rhovec[0] / rhovec.sum();
- auto derivs = get_derivs(model, T, rhovec);
+ auto derivs = get_derivs(model, T, rhovec, {});
// First two are residuals on critical point
return (Eigen::ArrayXd(2) << derivs.tot[2], derivs.tot[3]).finished();
};
@@ -746,11 +741,16 @@ struct CriticalTracing {
/**
* \brief Calculate dp/dT along the critical locus at given T, rhovec
+ * \f[
+ * \rmderivsub{Y}{T}{\CRL} = \deriv{Y}{T}{\vec\rho} + \sum_j \deriv{Y}{\rho_j}{T,\rho_{k\neq j}} \rmderivsub{\rho_j}{T}{\CRL}
+ * \f]
+ * where the derivatives on the RHS without the subscript $\CRL$ are homogeneous derivatives taken at the mixture statepoint, and are NOT along the critical curve.
*/
static auto get_dp_dT_crit(const Model& model, const Scalar& T, const VecType& rhovec) {
using id = IsochoricDerivatives<Model, Scalar, VecType>;
auto dpdrhovec = id::get_dpdrhovec_constT(model, T, rhovec);
- return id::get_dpdT_constrhovec(model, T, rhovec) + (dpdrhovec.matrix() * get_drhovec_dT_crit(model, T, rhovec).matrix()).array();
+ Scalar v = id::get_dpdT_constrhovec(model, T, rhovec) + (dpdrhovec.array()*get_drhovec_dT_crit(model, T, rhovec).array()).sum();
+ return v;
}
}; // namespace VecType
diff --git a/include/teqp/algorithms/critical_tracing_types.hpp b/include/teqp/algorithms/critical_tracing_types.hpp
index 48703d78533c36e75ef3d4e918357f414e601b8a..21d9192ba25221578d0548cebd6c3a3e3429e88d 100644
--- a/include/teqp/algorithms/critical_tracing_types.hpp
+++ b/include/teqp/algorithms/critical_tracing_types.hpp
@@ -24,4 +24,9 @@ struct TCABOptions {
bool pure_endpoint_polish = false; ///< If true, if the last step crossed into negative concentrations, try to interpolate to find the pure fluid endpoint hiding in the data
};
+struct EigenData {
+ Eigen::ArrayXd v0, v1, eigenvalues;
+ Eigen::MatrixXd eigenvectorscols;
+};
+
}
diff --git a/include/teqp/cpp/teqpcpp.hpp b/include/teqp/cpp/teqpcpp.hpp
index e94cc37956256cf3bfd82e1f038179c00d987c12..334bb86c6ce07f04c15d2ecc1a050861da93d0af 100644
--- a/include/teqp/cpp/teqpcpp.hpp
+++ b/include/teqp/cpp/teqpcpp.hpp
@@ -87,10 +87,6 @@ namespace teqp {
virtual double get_R(const EArrayd&) 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 EArray33d get_deriv_mat2(const double T, double rho, const EArrayd& z) const = 0;
-
virtual double get_Arxy(const int, const int, const double, const double, const EArrayd&) const = 0;
// Here X-Macros are used to create functions like get_Ar00, get_Ar01, ....
@@ -120,6 +116,8 @@ namespace teqp {
ISOCHORIC_matrix_args
#undef X
+ virtual EArray33d get_deriv_mat2(const double T, double rho, const EArrayd& z ) const = 0;
+
virtual std::tuple<double, double> solve_pure_critical(const double T, const double rho, const std::optional<nlohmann::json>&) const = 0;
virtual std::tuple<double, double> extrapolate_from_critical(const double Tc, const double rhoc, const double Tgiven) const = 0;
virtual EArray2 pure_VLE_T(const double T, const double rhoL, const double rhoV, int maxiter) const = 0;
@@ -128,6 +126,12 @@ namespace teqp {
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_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;
+ virtual double get_dp_dT_crit(const double T, const REArrayd& rhovec) const = 0;
+ virtual EArray2 get_criticality_conditions(const double T, const REArrayd& rhovec) const = 0;
+ virtual EigenData eigen_problem(const double T, const REArrayd& rhovec, const std::optional<REArrayd>&) const = 0;
+ virtual double get_minimum_eigenvalue_Psi_Hessian(const double T, const REArrayd& rhovec) const = 0;
};
diff --git a/interface/CPP/teqpcpp.cpp b/interface/CPP/teqpcpp.cpp
index 59f4c637de6d212969aad5872effbbae178eadc8..aeb5b343c0920e23378a92172409bcb93e7e7757 100644
--- a/interface/CPP/teqpcpp.cpp
+++ b/interface/CPP/teqpcpp.cpp
@@ -44,6 +44,37 @@ namespace teqp {
return crit::trace_critical_arclength_binary(model, T0, rhovec0, "");
}, m_model);
}
+ EArrayd get_drhovec_dT_crit(const double T, const REArrayd& rhovec) const override {
+ return std::visit([&](const auto& model) {
+ using crit = teqp::CriticalTracing<decltype(model), double, std::decay_t<decltype(rhovec)>>;
+ return crit::get_drhovec_dT_crit(model, T, rhovec);
+ }, m_model);
+ }
+ double get_dp_dT_crit(const double T, const REArrayd& rhovec) const override {
+ return std::visit([&](const auto& model) {
+ using crit = teqp::CriticalTracing<decltype(model), double, std::decay_t<decltype(rhovec)>>;
+ return crit::get_dp_dT_crit(model, T, rhovec);
+ }, m_model);
+ }
+ EArray2 get_criticality_conditions(const double T, const REArrayd& rhovec) const override {
+ return std::visit([&](const auto& model) {
+ using crit = teqp::CriticalTracing<decltype(model), double, std::decay_t<decltype(rhovec)>>;
+ return crit::get_criticality_conditions(model, T, rhovec);
+ }, m_model);
+ }
+ EigenData eigen_problem(const double T, const REArrayd& rhovec, const std::optional<REArrayd>& alignment_v0) const override {
+ return std::visit([&](const auto& model) {
+ using crit = teqp::CriticalTracing<decltype(model), double, std::decay_t<decltype(rhovec)>>;
+ return crit::eigen_problem(model, T, rhovec, alignment_v0.value_or(Eigen::ArrayXd()));
+ }, m_model);
+ }
+ double get_minimum_eigenvalue_Psi_Hessian(const double T, const REArrayd& rhovec) const override {
+ return std::visit([&](const auto& model) {
+ using crit = teqp::CriticalTracing<decltype(model), double, std::decay_t<decltype(rhovec)>>;
+ return crit::get_minimum_eigenvalue_Psi_Hessian(model, T, rhovec);
+ }, m_model);
+ }
+
EArray2 pure_VLE_T(const double T, const double rhoL, const double rhoV, int maxiter) const override {
return std::visit([&](const auto& model) {
return teqp::pure_VLE_T(model, T, rhoL, rhoV, maxiter);
diff --git a/interface/pybind11_wrapper.cpp b/interface/pybind11_wrapper.cpp
index 7746032220d09d42fc4c5716fac4c3375cf2a7cf..99a637ac754a784a9083ac302cc41e3f7545a18c 100644
--- a/interface/pybind11_wrapper.cpp
+++ b/interface/pybind11_wrapper.cpp
@@ -227,25 +227,25 @@ void init_teqp(py::module& m) {
// Routines related to binary mixture critical curve tracing
.def("trace_critical_arclength_binary", &am::trace_critical_arclength_binary, "T0"_a, "rhovec0"_a, py::arg_v("path", std::nullopt, "None"), py::arg_v("options", std::nullopt, "None"))
-// .def("get_criticality_conditions", &am::get_criticality_conditions)
-// .def("eigen_problem", &am::eigen_problem)
-// .def("get_minimum_eigenvalue_Psi_Hessian", &am::get_minimum_eigenvalue_Psi_Hessian)
-// .def("get_drhovec_dT_crit", &am::get_drhovec_dT_crit)
-// .def("get_dp_dT_crit", &am::get_dp_dT_crit)
+ .def("get_criticality_conditions", &am::get_criticality_conditions, "T"_a, "rhovec"_a.noconvert())
+ .def("eigen_problem", &am::eigen_problem, "T"_a, "rhovec"_a, py::arg_v("alignment_v0", std::nullopt, "None"))
+ .def("get_minimum_eigenvalue_Psi_Hessian", &am::get_minimum_eigenvalue_Psi_Hessian, "T"_a, "rhovec"_a.noconvert())
+ .def("get_drhovec_dT_crit", &am::get_drhovec_dT_crit, "T"_a, "rhovec"_a.noconvert())
+ .def("get_dp_dT_crit", &am::get_dp_dT_crit, "T"_a, "rhovec"_a.noconvert())
.def("pure_VLE_T", &am::pure_VLE_T, "T"_a, "rhoL"_a, "rhoV"_a, "max_iter"_a)
.def("get_drhovecdp_Tsat", &am::get_drhovecdp_Tsat, "T"_a, "rhovecL"_a.noconvert(), "rhovecV"_a.noconvert())
.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("get_drhovecdp_Tsat", &get_drhovecdp_Tsat<Model, double, RAX>, "T"_a, "rhovecL"_a.noconvert(), "rhovecV"_a.noconvert())
-// .def("trace_VLE_isotherm_binary", &trace_VLE_isotherm_binary<Model, double, Eigen::ArrayXd>, "T"_a, "rhovecL0"_a.noconvert(), "rhovecV0"_a.noconvert(), py::arg_v("options", std::nullopt, "None"))
-// .def("get_drhovecdT_psat", &get_drhovecdT_psat<Model, double, RAX>, "T"_a, "rhovecL"_a.noconvert(), "rhovecV"_a.noconvert())
-// .def("trace_VLE_isobar_binary", &trace_VLE_isobar_binary<Model, double, Eigen::ArrayXd>, "p"_a, "T0"_a, "rhovecL0"_a.noconvert(), "rhovecV0"_a.noconvert(), py::arg_v("options", std::nullopt, "None"))
-// .def("get_dpsat_dTsat_isopleth", &get_dpsat_dTsat_isopleth<Model, double, RAX>, "T"_a, "rhovecL"_a.noconvert(), "rhovecV"_a.noconvert())
-// .def("mix_VLLE_T", &mix_VLLE_T<Model, double, RAX>);
+
+// .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("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"));
-
;
m.def("_make_model", &teqp::cppinterface::make_model);