diff --git a/include/teqp/algorithms/VLE.hpp b/include/teqp/algorithms/VLE.hpp
index 3a17d69c6e90fc2ba4f71db5ddabd146c651442f..94e2a271371053d1fae597bd63d9419406ce0246 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,8 +162,10 @@ 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);
}
@@ -182,9 +183,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)>;
@@ -198,6 +200,78 @@ auto dpsatdT_pure(const Model& model, Scalar T, Scalar rhoL, Scalar rhoV) {
return dpsatdT;
}
+/***
+ \brief Starting at the critical point, trace the VLE down to a temperature of interest
+
+ \note This method only works for well-behaved EOS, notably absent from that category are EOS in the multiparameter category with orthobaric scaling exponent not equal to 0.5 at the critical point. Most other analytical EOS work fine
+
+ The JSON data structure defines the variables that need to be specified.
+
+ In the current implementation, there are a few steps:
+ 1. Solve for the true critical point satisfying \f$(\partial p/\partial \rho)_{T}=(\partial^2p/\partial\rho^2)_{T}=0\f$
+ 2. Take a small step away from the critical point (this is where the beta=0.5 assumption is invoked)
+ 3. Integrate from the near critical temperature to the temperature of interest
+ */
+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");
+ double R = model.R(z);
+ bool with_deriv = spec.at("with_deriv");
+ double dT = -(Tclose-T)/(Nstep-1);
+ using tdx = TDXDerivatives<decltype(model)>;
+ for (auto T_: Eigen::ArrayXd::LinSpaced(Nstep, Tclose, T)){
+ rhoLrhoVpolished = pure_VLE_T(model, T_, rhoLrhoVpolished[0], rhoLrhoVpolished[1], spec.value("NVLE", 10), z);
+
+ //auto pL = rhoLrhoVpolished[0]*R*T_*(1 + tdx::get_Ar01(model, T_, rhoLrhoVpolished[0], z));
+ //auto pV = rhoLrhoVpolished[1]*R*T_*(1 + tdx::get_Ar01(model, T_, rhoLrhoVpolished[1], z));
+ //std::cout << pL << " " << pV << " " << pL/pV-1 << std::endl;
+ if (with_deriv){
+ // Get drho/dT for both phases
+ auto dpsatdT = dpsatdT_pure(model, T_, rhoLrhoVpolished[0], rhoLrhoVpolished[1], z);
+ auto get_drhodT = [&z, R, dpsatdT](const Model& model, double T, double rho){
+ auto dpdrho = R*T*(1 + 2*tdx::get_Ar01(model, T, rho, z) + tdx::get_Ar02(model, T, rho, z));
+ auto dpdT = R*rho*(1 + tdx::get_Ar01(model, T, rho, z) - tdx::get_Ar11(model, T, rho, z));
+ return -dpdT/dpdrho + dpsatdT/dpdrho;
+ };
+ auto drhodTL = get_drhodT(model, T_, rhoLrhoVpolished[0]);
+ auto drhodTV = get_drhodT(model, T_, rhoLrhoVpolished[1]);
+ // Use the obtained derivative to calculate the step in rho from deltarho = (drhodT)*dT
+ auto DeltarhoL = dT*drhodTL, DeltarhoV = dT*drhodTV;
+ rhoLrhoVpolished[0] += DeltarhoL;
+ rhoLrhoVpolished[1] += DeltarhoV;
+ }
+
+ // Updated values for densities at new T
+ 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 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
@@ -939,7 +1013,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 +1242,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/cpp/teqpcpp.hpp b/include/teqp/cpp/teqpcpp.hpp
index b1981235274304f7613177a8c260c7a99061499e..ab6ddd5a483dde8d0a9d6eb253a533a89e6ce809 100644
--- a/include/teqp/cpp/teqpcpp.hpp
+++ b/include/teqp/cpp/teqpcpp.hpp
@@ -109,13 +109,13 @@ namespace teqp {
virtual double get_dmBnvirdTm(const int Nderiv, const int NTderiv, const double T, const EArrayd& z) const = 0;
// Derivatives from isochoric thermodynamics (all have the same signature whithin each block)
- #define X(f) virtual double f(const double T, const REArrayd& rhovec) const = 0;
+ #define X(f) virtual double f(const double T, const EArrayd& rhovec) const = 0;
ISOCHORIC_double_args
#undef X
- #define X(f) virtual EArrayd f(const double T, const REArrayd& rhovec) const = 0;
+ #define X(f) virtual EArrayd f(const double T, const EArrayd& rhovec) const = 0;
ISOCHORIC_array_args
#undef X
- #define X(f) virtual EMatrixd f(const double T, const REArrayd& rhovec) const = 0;
+ #define X(f) virtual EMatrixd f(const double T, const EArrayd& rhovec) const = 0;
ISOCHORIC_matrix_args
#undef X
diff --git a/include/teqp/derivs.hpp b/include/teqp/derivs.hpp
index 8c26ad8cd5d414c85d3fa4d93f9ead7ba4d6e5ac..c6df0701b41e14999a5cf2aa93751a4ce9a69acf 100644
--- a/include/teqp/derivs.hpp
+++ b/include/teqp/derivs.hpp
@@ -191,7 +191,9 @@ struct TDXDerivatives {
if constexpr (be == ADBackends::autodiff) {
using adtype = autodiff::HigherOrderDual<iT + iD, double>;
adtype Trecipad = 1.0 / T, rhoad = rho;
- auto f = [&w, &molefrac](const adtype& Trecip, const adtype& rho_) { return eval(w.alpha(eval(1.0/Trecip), rho_, molefrac)); };
+ auto f = [&w, &molefrac](const adtype& Trecip, const adtype& rho_) {
+ adtype T_ = 1.0/Trecip;
+ return eval(w.alpha(T_, rho_, molefrac)); };
auto wrts = std::tuple_cat(build_duplicated_tuple<iT>(std::ref(Trecipad)), build_duplicated_tuple<iD>(std::ref(rhoad)));
auto der = derivatives(f, std::apply(wrt_helper(), wrts), at(Trecipad, rhoad));
return powi(1.0 / T, iT) * powi(rho, iD) * der[der.size() - 1];
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
diff --git a/include/teqp/json_builder.hpp b/include/teqp/json_builder.hpp
index d3d78fa1dce88fd1ffe439bec2bf700476aa9948..079954752d49a5b96e16c05ad31924a53477adaa 100644
--- a/include/teqp/json_builder.hpp
+++ b/include/teqp/json_builder.hpp
@@ -95,6 +95,35 @@ namespace teqp {
throw std::invalid_argument("you must provide names or coeffs, but not both");
}
}
+ else if (kind == "SAFT-VR-Mie") {
+ using namespace SAFTVRMie;
+ std::optional<Eigen::ArrayXXd> kmat;
+ if (spec.contains("kmat")){
+ kmat = build_square_matrix(spec["kmat"]);
+ }
+
+ if (spec.contains("names")){
+ return SAFTVRMieMixture(spec["names"], kmat);
+ }
+ else if (spec.contains("coeffs")){
+ std::vector<SAFTVRMieCoeffs> coeffs;
+ for (auto j : spec["coeffs"]) {
+ SAFTVRMieCoeffs c;
+ c.name = j.at("name");
+ c.m = j.at("m");
+ c.sigma_m = (j.contains("sigma_m")) ? j.at("sigma_m").get<double>() : j.at("sigma_Angstrom").get<double>()/1e10;
+ c.epsilon_over_k = j.at("epsilon_over_k");
+ c.lambda_r = j.at("lambda_r");
+ c.lambda_a = j.at("lambda_a");
+ c.BibTeXKey = j.at("BibTeXKey");
+ coeffs.push_back(c);
+ }
+ return SAFTVRMieMixture(coeffs, kmat);
+ }
+ else{
+ throw std::invalid_argument("you must provide names or coeffs, but not both");
+ }
+ }
else if (kind == "multifluid") {
return multifluidfactory(spec);
}
diff --git a/include/teqp/math/quadrature.hpp b/include/teqp/math/quadrature.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..3b0433c537e1ad7c5610ac190e1c76c07aa1476c
--- /dev/null
+++ b/include/teqp/math/quadrature.hpp
@@ -0,0 +1,65 @@
+/**
+
+ Type-agnostic methods for carrying out fixed quadrature with definite integrals
+*/
+
+#pragma once
+
+namespace teqp{
+
+/**
+ Gauss-Legendre quadrature for a function f(x) in the interval [a,b]
+
+ More coefficients here if needed: https://pomax.github.io/bezierinfo/legendre-gauss.html
+*/
+template<int N, typename T>
+inline auto quad(const std::function<T(double)>& F, const double& a, const double& b){
+
+ // Locations x in [-1,1] where the function is to be evaluated, and the corresponding weight
+ static const std::map<int, std::tuple<std::vector<double>, std::vector<double>>> xw = {
+ {
+ 3,
+ {{ 0, sqrt(3.0/5), -sqrt(3.0/5) },
+ {8.0/9.0, 5.0/9.0, 5.0/9.0}}
+ },
+ {
+ 4,
+ {{sqrt(3.0/7.0-2.0/7.0*sqrt(6.0/5)), -sqrt(3.0/7.0-2.0/7.0*sqrt(6.0/5)), sqrt(3.0/7.0+2.0/7.0*sqrt(6.0/5)), -sqrt(3.0/7.0+2.0/7.0*sqrt(6.0/5))},
+ {(18.0+sqrt(30.0))/36.0, (18.0+sqrt(30.0))/36.0, (18.0-sqrt(30.0))/36.0, (18.0-sqrt(30.0))/36.0}}
+ },
+ {
+ 5,
+ {{0, 1.0/3.0*sqrt(5-2*sqrt(10.0/7.0)), -1.0/3.0*sqrt(5-2*sqrt(10.0/7.0)), 1.0/3.0*sqrt(5+2*sqrt(10.0/7.0)), -1.0/3.0*sqrt(5+2*sqrt(10.0/7.0))},
+ {128/225.0, (322.0+13*sqrt(70))/900.0, (322.0+13*sqrt(70))/900.0, (322.0-13*sqrt(70))/900.0, (322.0-13*sqrt(70))/900.0}}
+ },
+ {
+ 7,
+ {{0.0000000000000000,0.4058451513773972,-0.4058451513773972,-0.7415311855993945,0.7415311855993945,-0.9491079123427585,0.9491079123427585},
+ { 0.4179591836734694,0.3818300505051189,0.3818300505051189,0.2797053914892766,0.2797053914892766,0.1294849661688697,0.1294849661688697}}
+ },
+ {
+ 15,
+ {{0.0000000000000000,-0.2011940939974345,0.2011940939974345,-0.3941513470775634,0.3941513470775634,-0.5709721726085388,0.5709721726085388,-0.7244177313601701,0.7244177313601701,-0.8482065834104272,0.8482065834104272,-0.9372733924007060,0.9372733924007060,-0.9879925180204854,0.9879925180204854},
+ {0.2025782419255613,0.1984314853271116,0.1984314853271116,0.1861610000155622,0.1861610000155622,0.1662692058169939,0.1662692058169939,0.1395706779261543,0.1395706779261543,0.1071592204671719,0.1071592204671719,0.0703660474881081,0.0703660474881081,0.0307532419961173,0.0307532419961173}}
+ },
+ {
+ 30,
+ {{-0.0514718425553177,0.0514718425553177,-0.1538699136085835,0.1538699136085835,-0.2546369261678899,0.2546369261678899,-0.3527047255308781,0.3527047255308781,-0.4470337695380892,0.4470337695380892,-0.5366241481420199,0.5366241481420199,-0.6205261829892429,0.6205261829892429,-0.6978504947933158,0.6978504947933158,-0.7677774321048262,0.7677774321048262,-0.8295657623827684,0.8295657623827684,-0.8825605357920527,0.8825605357920527,-0.9262000474292743,0.9262000474292743,-0.9600218649683075,0.9600218649683075,-0.9836681232797472,0.9836681232797472,-0.9968934840746495,0.9968934840746495
+ },
+ {0.1028526528935588,0.1028526528935588,0.1017623897484055,0.1017623897484055,0.0995934205867953,0.0995934205867953,0.0963687371746443,0.0963687371746443,0.0921225222377861,0.0921225222377861,0.0868997872010830,0.0868997872010830,0.0807558952294202,0.0807558952294202,0.0737559747377052,0.0737559747377052,0.0659742298821805,0.0659742298821805,0.0574931562176191,0.0574931562176191,0.0484026728305941,0.0484026728305941,0.0387991925696271,0.0387991925696271,0.0287847078833234,0.0287847078833234,0.0184664683110910,0.0184664683110910,0.0079681924961666,0.0079681924961666}}
+ }
+ };
+
+ T summer = 0.0;
+ // Lookup the coefficients without making a copy or re-allocating
+ const std::tuple<std::vector<double>, std::vector<double>> &pair = xw.at(N);
+ const std::vector<double> &x = std::get<0>(pair);
+ const std::vector<double> &w = std::get<1>(pair);
+ for (auto i = 0; i < N; ++i){
+ double arg = (b-a)/2.0*x[i] + (a+b)/2.0;
+ summer += w[i]*F(arg);
+ }
+ T retval = (b-a)/2.0*summer; // Forces a flattening if T is an autodiff type
+ return retval;
+}
+}
diff --git a/include/teqp/models/fwd.hpp b/include/teqp/models/fwd.hpp
index 03f9781448d7e73be3dfb201620f26a400775ed5..0ece52c686c6151e8a9b86b01b7cdd7bc4eda9d5 100644
--- a/include/teqp/models/fwd.hpp
+++ b/include/teqp/models/fwd.hpp
@@ -15,6 +15,7 @@
#include "teqp/models/cubics.hpp"
#include "teqp/models/CPA.hpp"
#include "teqp/models/pcsaft.hpp"
+#include "teqp/models/saftvrmie.hpp"
#include "teqp/models/multifluid.hpp"
#include "teqp/models/multifluid_mutant.hpp"
#include "teqp/ideal_eosterms.hpp"
@@ -28,11 +29,13 @@
namespace teqp {
using vad = std::valarray<double>;
+ using vecs = std::vector<std::string>;
// Define the EOS types by interrogating the types returned by the respective
// factory function or by alias of the class name
using canonical_cubic_t = decltype(canonical_PR(vad{}, vad{}, vad{}));
using PCSAFT_t = decltype(PCSAFT::PCSAFTfactory(nlohmann::json{}));
+ using SAFTVRMie_t = decltype(SAFTVRMie::SAFTVRMieMixture(vecs{}));
using CPA_t = decltype(CPA::CPAfactory(nlohmann::json{}));
using multifluid_t = decltype(multifluidfactory(nlohmann::json{}));
using multifluidmutant_t = decltype(build_multifluid_mutant(multifluidfactory(nlohmann::json{}), nlohmann::json{}));
@@ -53,6 +56,7 @@ namespace teqp {
vdWEOS_t,
canonical_cubic_t,
PCSAFT_t,
+ SAFTVRMie_t,
CPA_t,
multifluid_t,
multifluidmutant_t,
diff --git a/include/teqp/models/saftvrmie.hpp b/include/teqp/models/saftvrmie.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..cf6740af8210007e7d90c2c3cd3052384b3609c6
--- /dev/null
+++ b/include/teqp/models/saftvrmie.hpp
@@ -0,0 +1,694 @@
+/***
+
+ \brief This file contains the contributions that can be composed together to form SAFT models
+
+*/
+
+#pragma once
+
+#include "nlohmann/json.hpp"
+#include "teqp/types.hpp"
+#include "teqp/exceptions.hpp"
+#include "teqp/constants.hpp"
+#include "teqp/math/quadrature.hpp"
+#include <optional>
+
+namespace teqp {
+namespace SAFTVRMie {
+
+/// Coefficients for one fluid
+struct SAFTVRMieCoeffs {
+ std::string name; ///< Name of fluid
+ double m = -1, ///< number of segments
+ sigma_m = -1, ///< [m] segment diameter
+ epsilon_over_k = -1, ///< [K] depth of pair potential divided by Boltzman constant
+ lambda_a = -1, ///< The attractive exponent (the 6 in LJ 12-6 potential)
+ lambda_r = -1; ///< The repulsive exponent (the 12 in LJ 12-6 potential)
+ std::string BibTeXKey; ///< The BibTeXKey for the reference for these coefficients
+};
+
+/// Manager class for SAFT-VR-Mie coefficients
+class SAFTVRMieLibrary {
+ std::map<std::string, SAFTVRMieCoeffs> coeffs;
+public:
+ SAFTVRMieLibrary() {
+ insert_normal_fluid("Methane", 1.0000, 3.7412e-10, 153.36, 12.650, 6, "Lafitte-JCP-2001");
+ insert_normal_fluid("Ethane", 1.4373, 3.7257e-10, 206.12, 12.400, 6, "Lafitte-JCP-2001");
+ insert_normal_fluid("Propane", 1.6845, 3.9056e-10, 239.89, 13.006, 6, "Lafitte-JCP-2001");
+ }
+ void insert_normal_fluid(const std::string& name, double m, const double sigma_m, const double epsilon_over_k, const double lambda_a, const double lambda_r, const std::string& BibTeXKey) {
+ SAFTVRMieCoeffs coeff;
+ coeff.name = name;
+ coeff.m = m;
+ coeff.sigma_m = sigma_m;
+ coeff.epsilon_over_k = epsilon_over_k;
+ coeff.lambda_r = lambda_r;
+ coeff.lambda_a = lambda_a;
+ coeff.BibTeXKey = BibTeXKey;
+ coeffs.insert(std::pair<std::string, SAFTVRMieCoeffs>(name, coeff));
+ }
+ const auto& get_normal_fluid(const std::string& name) {
+ auto it = coeffs.find(name);
+ if (it != coeffs.end()) {
+ return it->second;
+ }
+ else {
+ throw std::invalid_argument("Bad name:" + name);
+ }
+ }
+ auto get_coeffs(const std::vector<std::string>& names){
+ std::vector<SAFTVRMieCoeffs> c;
+ for (auto n : names){
+ c.push_back(get_normal_fluid(n));
+ }
+ return c;
+ }
+};
+
+/// Things that only depend on the components themselves, but not on composition, temperature, or density
+struct SAFTVRMieChainContributionTerms{
+ private:
+
+ /// The matrix of coefficients needed to evaluate f_k
+ const Eigen::Matrix<double, 7, 7> phi{(Eigen::Matrix<double, 7, 7>() <<
+ 7.5365557, -359.44, 1550.9, -1.19932, -1911.28, 9236.9, 10,
+ -37.60463, 1825.6, -5070.1, 9.063632, 21390.175, -129430, 10,
+ 71.745953, -3168.0, 6534.6, -17.9482, -51320.7, 357230, 0.57,
+ -46.83552, 1884.2, -3288.7, 11.34027, 37064.54, -315530, -6.7,
+ -2.467982, -0.82376, -2.7171, 20.52142, 1103.742, 1390.2, -8,
+ -0.50272, -3.1935, 2.0883, -56.6377, -3264.61, -4518.2, 0,
+ 8.0956883, 3.7090, 0, 40.53683, 2556.181, 4241.6, 0 ).finished()};
+
+ /// The matrix used to obtain the parameters c_1, c_2, c_3, and c_4 in Eq. A18
+ const Eigen::Matrix<double, 4, 4> A{(Eigen::Matrix<double, 4, 4>() <<
+ 0.81096, 1.7888, -37.578, 92.284,
+ 1.0205, -19.341, 151.26, -463.50,
+ -1.9057, 22.845, -228.14, 973.92,
+ 1.0885, -6.1962, 106.98, -677.64).finished()};
+
+ // Eq. A48
+ auto get_lambda_k_ij(const Eigen::ArrayXd& lambda_k) const{
+ Eigen::ArrayXXd mat(N,N);
+ for (auto i = 0; i < lambda_k.size(); ++i){
+ for (auto j = i; j < lambda_k.size(); ++j){
+ mat(i,j) = 3 + sqrt((lambda_k(i)-3)*(lambda_k(j)-3));
+ mat(j,i) = mat(i,j);
+ }
+ }
+ return mat;
+ }
+
+ /// Eq. A3
+ auto get_C_ij() const{
+ Eigen::ArrayXXd C(N,N);
+ for (auto i = 0; i < N; ++i){
+ for (auto j = i; j < N; ++j){
+ C(i,j) = lambda_r_ij(i,j)/(lambda_r_ij(i,j)-lambda_a_ij(i,j))*pow(lambda_r_ij(i,j)/lambda_a_ij(i,j), lambda_a_ij(i,j)/(lambda_r_ij(i,j)-lambda_a_ij(i,j)));
+ C(j,i) = C(i,j); // symmetric
+ }
+ }
+ return C;
+ }
+
+ // Eq. A26
+ auto get_fkij() const{
+ std::vector<Eigen::ArrayXXd> f_(8); // 0-th element is present, but not initialized
+ for (auto k = 1; k < 8; ++k){
+ f_[k].resize(N,N);
+ };
+ for (auto k = 1; k < 8; ++k){
+ auto phik = phi.col(k-1); // phi matrix is indexed to start at 1, but our matrix starts at 0
+ Eigen::ArrayXXd num(N,N), den(N,N); num.setZero(), den.setZero();
+ for (auto n = 0; n < 4; ++n){
+ num += phik[n]*pow(alpha_ij, n);
+ }
+ for (auto n = 4; n < 7; ++n){
+ den += phik[n]*pow(alpha_ij, n-3);
+ }
+ f_[k] = num/(1 + den);
+ }
+ return f_;
+ }
+
+ /// Eq. A45
+ auto get_sigma_ij() const{
+ Eigen::ArrayXXd sigma(N,N);
+ for (auto i = 0; i < N; ++i){
+ for (auto j = i; j < N; ++j){
+ sigma(i,j) = (sigma_A(i) + sigma_A(j))/2.0;
+ sigma(j,i) = sigma(i,j); // symmetric
+ }
+ }
+ return sigma;
+ }
+ /// Eq. A55
+ auto get_epsilon_ij() const{
+ Eigen::ArrayXXd eps_(N,N);
+ for (auto i = 0; i < N; ++i){
+ for (auto j = i; j < N; ++j){
+ eps_(i,j) = (1.0-kmat(i,j))*sqrt(pow(sigma_ij(i,i),3)*pow(sigma_ij(j,j),3)*epsilon_over_k(i)*epsilon_over_k(j))/pow(sigma_ij(i,j), 3);
+ eps_(j,i) = eps_(i,j); // symmetric
+ }
+ }
+ return eps_;
+ }
+ std::size_t get_N(){
+ auto sizes = (Eigen::ArrayXd(5) << m.size(), epsilon_over_k.size(), sigma_A.size(), lambda_a.size(), lambda_r.size()).finished();
+ if (sizes.mean() != sizes.minCoeff()){
+ throw teqp::InvalidArgument("sizes of pure component arrays are not all the same");
+ }
+ return sizes[0];
+ }
+
+ /// Eq. A18 for the attractive exponents
+ auto get_cij(const Eigen::ArrayXXd& lambdaij) const{
+ std::vector<Eigen::ArrayXXd> cij(4);
+ for (auto n = 0; n < 4; ++n){
+ cij[n].resize(N,N);
+ };
+ for (auto i = 0; i < N; ++i){
+ for (auto j = i; j < N; ++j){
+ using CV = Eigen::Vector<double, 4>;
+ const CV b{(CV() << 1, 1.0/lambdaij(i,j), 1.0/pow(lambdaij(i,j),2), 1.0/pow(lambdaij(i,j),3)).finished()};
+ auto c1234 = (A*b).eval();
+ cij[0](i,j) = c1234(0);
+ cij[1](i,j) = c1234(1);
+ cij[2](i,j) = c1234(2);
+ cij[3](i,j) = c1234(3);
+ }
+ }
+ return cij;
+ }
+
+ /// Eq. A18 for the attractive exponents
+ auto get_canij() const{
+ return get_cij(lambda_a_ij);
+ }
+ /// Eq. A18 for 2x the attractive exponents
+ auto get_c2anij() const{
+ return get_cij(2.0*lambda_a_ij);
+ }
+ /// Eq. A18 for the repulsive exponents
+ auto get_crnij() const{
+ return get_cij(lambda_r_ij);
+ }
+ /// Eq. A18 for the 2x the repulsive exponents
+ auto get_c2rnij() const{
+ return get_cij(2.0*lambda_r_ij);
+ }
+ /// Eq. A18 for the 2x the repulsive exponents
+ auto get_carnij() const{
+ return get_cij(lambda_r_ij + lambda_a_ij);
+ }
+
+ template<typename X>
+ auto POW2(const X& x) const{
+ return forceeval(x*x);
+ };
+ template<typename X>
+ auto POW3(const X& x) const{
+ return forceeval(x*POW2(x));
+ };
+ template<typename X>
+ auto POW4(const X& x) const{
+ return forceeval(POW2(x)*POW2(x));
+ };
+
+ public:
+
+ // One entry per component
+ const Eigen::ArrayXd m, epsilon_over_k, sigma_A, lambda_a, lambda_r;
+ const Eigen::ArrayXXd kmat;
+
+ const std::size_t N;
+
+ // Calculated matrices for the ij pair
+ const Eigen::ArrayXXd lambda_r_ij, lambda_a_ij, C_ij, alpha_ij, sigma_ij, epsilon_ij; // Matrices of parameters
+
+ const std::vector<Eigen::ArrayXXd> canij, crnij, c2anij, c2rnij, carnij;
+ const std::vector<Eigen::ArrayXXd> fkij; // Matrices of parameters
+
+ SAFTVRMieChainContributionTerms(
+ const Eigen::ArrayXd& m,
+ const Eigen::ArrayXd& epsilon_over_k,
+ const Eigen::ArrayXd& sigma_m,
+ const Eigen::ArrayXd& lambda_r,
+ const Eigen::ArrayXd& lambda_a,
+ const Eigen::ArrayXXd& kmat)
+ : m(m), epsilon_over_k(epsilon_over_k), sigma_A(sigma_m*1e10), lambda_a(lambda_a), lambda_r(lambda_r), kmat(kmat),
+ N(get_N()),
+ lambda_r_ij(get_lambda_k_ij(lambda_r)), lambda_a_ij(get_lambda_k_ij(lambda_a)),
+ C_ij(get_C_ij()), alpha_ij(C_ij*(1/(lambda_a_ij-3) - 1/(lambda_r_ij-3))),
+ sigma_ij(get_sigma_ij()), epsilon_ij(get_epsilon_ij()),
+ crnij(get_crnij()), canij(get_canij()),
+ c2rnij(get_c2rnij()), c2anij(get_c2anij()), carnij(get_carnij()),
+ fkij(get_fkij())
+ {}
+
+ /// Eq. A2 from Lafitte
+ double get_uii_over_kB(std::size_t i, const double& r) const {
+ double rstarinv = sigma_A[i]/r;
+ return forceeval(C_ij(i,i)*epsilon_over_k[i]*(::pow(rstarinv, lambda_r[i]) - ::pow(rstarinv, lambda_a[i])));
+ }
+
+ /// Eq. A9 from Lafitte
+ template <typename TType>
+ TType get_dii(std::size_t i, const TType &T) const{
+ std::function<TType(double)> integrand = [this, i, &T](const double& r){
+ return forceeval(1.0-exp(-this->get_uii_over_kB(i, r)/T));
+ };
+ return quad<30, TType>(integrand, 0.0, sigma_A[i]);
+ }
+
+ template <typename TType>
+ auto get_dmat(const TType &T) const{
+ Eigen::Array<TType, Eigen::Dynamic, Eigen::Dynamic> d(N,N);
+ // For the pure components, by integration
+ for (auto i = 0; i < N; ++i){
+ d(i,i) = get_dii(i, T);
+ }
+ // The cross terms, using the linear mixing rule
+ for (auto i = 0; i < N; ++i){
+ for (auto j = i+1; j < N; ++j){
+ d(i,j) = (d(i,i) + d(j,j))/2.0;
+ d(j,i) = d(i,j);
+ }
+ }
+ return d;
+ }
+ // Calculate core parameters that depend on temperature, volume, and composition
+ template <typename TType, typename RhoType, typename VecType>
+ auto get_core_calcs(const TType& T, const RhoType& rhomolar, const VecType& molefracs) const{
+
+ if (molefracs.size() != N){
+ throw teqp::InvalidArgument("Length of molefracs of "+std::to_string(molefracs.size()) + " does not match the model size of"+std::to_string(N));
+ }
+
+ using FracType = std::decay_t<decltype(molefracs[0])>;
+ using NumType = std::common_type_t<TType, RhoType, FracType>;
+
+ // Things that are easy to calculate
+ // ....
+
+ auto dmat = get_dmat(T); // Matrix of diameters of pure and cross terms
+ auto rhoN = forceeval(rhomolar*N_A); // Number density, in molecules/m^3
+ auto mbar = forceeval((molefracs*m).sum()); // Mean number of segments, dimensionless
+ auto rhos = forceeval(rhoN*mbar/1e30); // Mean segment number density, in segments/A^3
+ auto xs = forceeval((m*molefracs/mbar).eval()); // Segment fractions
+
+ constexpr double MY_PI = static_cast<double>(EIGEN_PI);
+ auto pi6 = MY_PI/6;
+
+ using TRHOType = std::common_type_t<std::decay_t<TType>, std::decay_t<RhoType>, std::decay_t<decltype(molefracs[0])>, std::decay_t<decltype(m[0])>>;
+ Eigen::Array<TRHOType, 4, 1> zeta;
+ for (auto l = 0; l < 4; ++l){
+ TRHOType summer = 0.0;
+ for (auto i = 0; i < N; ++i){
+ summer += xs(i)*powi(dmat(i,i), l);
+ }
+ zeta(l) = forceeval(pi6*rhos*summer);
+ }
+
+ NumType summer_xi_x = 0.0;
+ TRHOType summer_xi_x_bar = 0.0;
+ for (auto i = 0; i < N; ++i){
+ for (auto j = 0; j < N; ++j){
+ summer_xi_x += xs(i)*xs(j)*powi(dmat(i,j), 3)*rhos;
+ summer_xi_x_bar += xs(i)*xs(j)*powi(sigma_ij(i,j), 3);
+ }
+ }
+
+ auto xi_x = forceeval(pi6*summer_xi_x); // Eq. A13
+ auto xi_x_bar = forceeval(pi6*rhos*summer_xi_x_bar); // Eq. A23
+ auto xi_x_bar5 = forceeval(POW2(xi_x_bar)*POW3(xi_x_bar)); // (xi_x_bar)^5
+ auto xi_x_bar8 = forceeval(xi_x_bar5*POW3(xi_x_bar)); // (xi_x_bar)^8
+
+ // Coefficients in the gdHSij term, do not depend on component,
+ // so calculate them here
+ auto X = forceeval(POW3(1.0 - xi_x)), X3 = X;
+ auto X2 = forceeval(POW2(1.0 - xi_x));
+ auto k0 = forceeval(-log(1.0-xi_x) + (42.0*xi_x - 39.0*POW2(xi_x) + 9.0*POW3(xi_x) - 2.0*POW4(xi_x))/(6.0*X3)); // Eq. A30
+ auto k1 = forceeval((POW4(xi_x) + 6.0*POW2(xi_x) - 12.0*xi_x)/(2.0*X3));
+ auto k2 = forceeval(-3.0*POW2(xi_x)/(8.0*X2));
+ auto k3 = forceeval((-POW4(xi_x) + 3.0*POW2(xi_x) + 3.0*xi_x)/(6.0*X3));
+
+ // Pre-calculate the cubes of the diameters
+ auto dmat3 = dmat.array().pow(3).eval();
+
+ NumType a1kB = 0.0;
+ NumType a2kB2 = 0.0;
+ NumType a3kB3 = 0.0;
+ NumType alphar_chain = 0.0;
+
+ NumType K_HS = get_KHS(xi_x);
+ NumType rho_dK_HS_drho = get_rhos_dK_HS_drhos(xi_x);
+
+ for (auto i = 0; i < N; ++i){
+ for (auto j = i; j < N; ++j){
+ NumType x_0_ij = sigma_ij(i,j)/dmat(i, j);
+
+ // -----------------------
+ // Calculations for a_1/kB
+ // -----------------------
+
+ auto I = [&x_0_ij](double lambda_ij){
+ return forceeval(-(pow(x_0_ij, 3-lambda_ij)-1.0)/(lambda_ij-3.0)); // Eq. A14
+ };
+ auto J = [&x_0_ij](double lambda_ij){
+ return forceeval(-(pow(x_0_ij, 4-lambda_ij)*(lambda_ij-3.0)-pow(x_0_ij, 3.0-lambda_ij)*(lambda_ij-4.0)-1.0)/((lambda_ij-3.0)*(lambda_ij-4.0))); // Eq. A15
+ };
+ auto Bhatij_a = this->get_Bhatij(xi_x, X, I(lambda_a_ij(i,j)), J(lambda_a_ij(i,j)));
+ auto Bhatij_2a = this->get_Bhatij(xi_x, X, I(2*lambda_a_ij(i,j)), J(2*lambda_a_ij(i,j)));
+ auto Bhatij_r = this->get_Bhatij(xi_x, X, I(lambda_r_ij(i,j)), J(lambda_r_ij(i,j)));
+ auto Bhatij_2r = this->get_Bhatij(xi_x, X, I(2*lambda_r_ij(i,j)), J(2*lambda_r_ij(i,j)));
+ auto Bhatij_ar = this->get_Bhatij(xi_x, X, I(lambda_a_ij(i,j)+lambda_r_ij(i,j)), J(lambda_a_ij(i,j)+lambda_r_ij(i,j)));
+
+ auto one_term = [this, &x_0_ij, &I, &J, &xi_x, &X](double lambda_ij, const NumType& xi_x_eff){
+ return forceeval(
+ pow(x_0_ij, lambda_ij)*(
+ this->get_Bhatij(xi_x, X, I(lambda_ij), J(lambda_ij))
+ + this->get_a1Shatij(xi_x_eff, lambda_ij)
+ )
+ );
+ };
+ NumType xi_x_eff_r = crnij[0](i,j)*xi_x + crnij[1](i,j)*POW2(xi_x) + crnij[2](i,j)*POW3(xi_x) + crnij[3](i,j)*POW4(xi_x);
+ NumType xi_x_eff_a = canij[0](i,j)*xi_x + canij[1](i,j)*POW2(xi_x) + canij[2](i,j)*POW3(xi_x) + canij[3](i,j)*POW4(xi_x);
+ NumType dxi_x_eff_dxix_r = crnij[0](i,j) + crnij[1](i,j)*2*xi_x + crnij[2](i,j)*3*POW2(xi_x) + crnij[3](i,j)*4*POW3(xi_x);
+ NumType dxi_x_eff_dxix_a = canij[0](i,j) + canij[1](i,j)*2*xi_x + canij[2](i,j)*3*POW2(xi_x) + canij[3](i,j)*4*POW3(xi_x);
+
+ NumType a1ij = 2.0*MY_PI*rhos*dmat3(i,j)*epsilon_ij(i,j)*C_ij(i,j)*(
+ one_term(lambda_a_ij(i,j), xi_x_eff_a) - one_term(lambda_r_ij(i,j), xi_x_eff_r)
+ ); // divided by k_B
+
+ NumType contribution = xs(i)*xs(j)*a1ij;
+ double factor = (i == j) ? 1.0 : 2.0; // Off-diagonal terms contribute twice
+ a1kB += contribution*factor;
+
+ // --------------------------
+ // Calculations for a_2/k_B^2
+ // --------------------------
+
+ NumType xi_x_eff_2r = c2rnij[0](i,j)*xi_x + c2rnij[1](i,j)*POW2(xi_x) + c2rnij[2](i,j)*POW3(xi_x) + c2rnij[3](i,j)*POW4(xi_x);
+ NumType xi_x_eff_2a = c2anij[0](i,j)*xi_x + c2anij[1](i,j)*POW2(xi_x) + c2anij[2](i,j)*POW3(xi_x) + c2anij[3](i,j)*POW4(xi_x);
+ NumType xi_x_eff_ar = carnij[0](i,j)*xi_x + carnij[1](i,j)*POW2(xi_x) + carnij[2](i,j)*POW3(xi_x) + carnij[3](i,j)*POW4(xi_x);
+ NumType dxi_x_eff_dxix_2r = c2rnij[0](i,j) + c2rnij[1](i,j)*2*xi_x + c2rnij[2](i,j)*3*POW2(xi_x) + c2rnij[3](i,j)*4*POW3(xi_x);
+ NumType dxi_x_eff_dxix_ar = carnij[0](i,j) + carnij[1](i,j)*2*xi_x + carnij[2](i,j)*3*POW2(xi_x) + carnij[3](i,j)*4*POW3(xi_x);
+ NumType dxi_x_eff_dxix_2a = c2anij[0](i,j) + c2anij[1](i,j)*2*xi_x + c2anij[2](i,j)*3*POW2(xi_x) + c2anij[3](i,j)*4*POW3(xi_x);
+
+ NumType chi_ij = fkij[1](i,j)*xi_x_bar + fkij[2](i,j)*xi_x_bar5 + fkij[3](i,j)*xi_x_bar8;
+ auto a2ij = 0.5*K_HS*(1.0+chi_ij)*epsilon_ij(i,j)*POW2(C_ij(i,j))*(2*MY_PI*rhos*dmat3(i,j)*epsilon_ij(i,j))*(
+ one_term(2.0*lambda_a_ij(i,j), xi_x_eff_2a)
+ -2.0*one_term(lambda_a_ij(i,j)+lambda_r_ij(i,j), xi_x_eff_ar)
+ +one_term(2.0*lambda_r_ij(i,j), xi_x_eff_2r)
+ ); // divided by k_B^2
+
+ NumType contributiona2 = xs(i)*xs(j)*a2ij; // Eq. A19
+ a2kB2 += contributiona2*factor;
+
+ // --------------------------
+ // Calculations for a_3/k_B^3
+ // --------------------------
+ auto a3ij = -POW3(epsilon_ij(i,j))*fkij[4](i,j)*xi_x_bar*exp(
+ fkij[5](i,j)*xi_x_bar + fkij[6](i,j)*POW2(xi_x_bar)
+ ); // divided by k_B^3
+ NumType contributiona3 = xs(i)*xs(j)*a3ij; // Eq. A25
+ a3kB3 += contributiona3*factor;
+
+ if (i == j){
+ // ------------------
+ // Chain contribution
+ // ------------------
+
+ // Eq. A29
+ auto gdHSii = exp(k0 + k1*x_0_ij + k2*POW2(x_0_ij) + k3*POW3(x_0_ij));
+
+ // The g1 terms
+ // ....
+
+ // This is the function for the second part (not the partial) that goes in g_{1,ii},
+ // divided by 2*PI*d_ij^3*epsilon*rhos
+ auto g1_term = [&one_term](double lambda_ij, const NumType& xi_x_eff){
+ return forceeval(lambda_ij*one_term(lambda_ij, xi_x_eff));
+ };
+ auto g1_noderivterm = -C_ij(i,i)*(g1_term(lambda_a_ij(i,i), xi_x_eff_a)-g1_term(lambda_r_ij(i,i), xi_x_eff_r));
+
+ // Bhat = B*rho*kappa; diff(Bhat, rho) = Bhat + rho*dBhat/drho; kappa = 2*pi*eps*d^3
+ // This is the function for the partial derivative rhos*(da1ij/drhos),
+ // divided by 2*PI*d_ij^3*epsilon*rhos
+ auto rhosda1iidrhos_term = [this, &x_0_ij, &I, &J, &xi_x, &X](double lambda_ij, const NumType& xi_x_eff, const NumType& dxixeff_dxix, const NumType& Bhatij){
+ auto I_ = I(lambda_ij);
+ auto J_ = J(lambda_ij);
+ auto rhosda1Sdrhos = this->get_rhoda1Shatijdrho(xi_x, xi_x_eff, dxixeff_dxix, lambda_ij);
+ auto rhosdBdrhos = this->get_rhodBijdrho(xi_x, X, I_, J_, Bhatij);
+ return forceeval(pow(x_0_ij, lambda_ij)*(rhosda1Sdrhos + rhosdBdrhos));
+ };
+ // This is rhos*d(a_1ij)/drhos/(2*pi*d^3*eps*rhos)
+ auto da1iidrhos_term = C_ij(i,j)*(
+ rhosda1iidrhos_term(lambda_a_ij(i,i), xi_x_eff_a, dxi_x_eff_dxix_a, Bhatij_a)
+ -rhosda1iidrhos_term(lambda_r_ij(i,i), xi_x_eff_r, dxi_x_eff_dxix_r, Bhatij_r)
+ );
+ auto g1ii = 3.0*da1iidrhos_term + g1_noderivterm;
+
+ // The g2 terms
+ // ....
+
+ // This is the second part (not the partial deriv.) that goes in g_{2,ii},
+ // divided by 2*PI*d_ij^3*epsilon*rhos
+ auto g2_noderivterm = -POW2(C_ij(i,i))*K_HS*(
+ lambda_a_ij(i,j)*one_term(2*lambda_a_ij(i,j), xi_x_eff_2a)
+ -(lambda_a_ij(i,j)+lambda_r_ij(i,j))*one_term(lambda_a_ij(i,j)+lambda_r_ij(i,j), xi_x_eff_ar)
+ +lambda_r_ij(i,j)*one_term(2*lambda_r_ij(i,j), xi_x_eff_2r)
+ );
+ // This is [rhos*d(a_2ij/(1+chi_ij))/drhos]/(2*pi*d^3*eps*rhos)
+ auto da2iidrhos_term = 0.5*POW2(C_ij(i,j))*(
+ rho_dK_HS_drho*(
+ one_term(2.0*lambda_a_ij(i,j), xi_x_eff_2a)
+ -2.0*one_term(lambda_a_ij(i,j)+lambda_r_ij(i,j), xi_x_eff_ar)
+ +one_term(2.0*lambda_r_ij(i,j), xi_x_eff_2r))
+ +K_HS*(
+ rhosda1iidrhos_term(2.0*lambda_a_ij(i,i), xi_x_eff_2a, dxi_x_eff_dxix_2a, Bhatij_2a)
+ -2.0*rhosda1iidrhos_term(lambda_a_ij(i,i)+lambda_r_ij(i,i), xi_x_eff_ar, dxi_x_eff_dxix_ar, Bhatij_ar)
+ +rhosda1iidrhos_term(2.0*lambda_r_ij(i,i), xi_x_eff_2r, dxi_x_eff_dxix_2r, Bhatij_2r)
+ )
+ );
+ auto g2MCAij = 3.0*da2iidrhos_term + g2_noderivterm;
+
+ auto betaepsilon = epsilon_ij(i,i)/T; // (1/(kB*T))/epsilon
+ auto theta = exp(betaepsilon)-1.0;
+ auto phi7 = phi.col(6);
+ auto gamma_cij = phi7(0)*(-tanh(phi7(1)*(phi7(2)-alpha_ij(i,j)))+1.0)*xi_x_bar*theta*exp(phi7(3)*xi_x_bar + phi7(4)*POW2(xi_x_bar)); // Eq. A37
+ auto g2ii = (1.0+gamma_cij)*g2MCAij;
+
+ NumType giiMie = gdHSii*exp((betaepsilon*g1ii + POW2(betaepsilon)*g2ii)/gdHSii);
+ alphar_chain -= molefracs[i]*(m[i]-1.0)*log(giiMie);
+ }
+ }
+ }
+
+ auto ahs = a_HS(rhos, zeta);
+ // Eq. A5 from Lafitte, multiplied by mbar
+ auto alphar_mono = forceeval(mbar*(ahs + a1kB/T + a2kB2/(T*T) + a3kB3/(T*T*T)));
+
+ struct vals{
+ decltype(dmat) dmat;
+ decltype(rhos) rhos;
+ decltype(rhoN) rhoN;
+ decltype(mbar) mbar;
+ decltype(xs) xs;
+ decltype(zeta) zeta;
+ decltype(xi_x) xi_x;
+ decltype(xi_x_bar) xi_x_bar;
+ decltype(alphar_mono) alphar_mono;
+ decltype(a1kB) a1kB;
+ decltype(a2kB2) a2kB2;
+ decltype(a3kB3) a3kB3;
+ decltype(alphar_chain) alphar_chain;
+ };
+ return vals{dmat, rhos, rhoN, mbar, xs, zeta, xi_x, xi_x_bar, alphar_mono, a1kB, a2kB2, a3kB3, alphar_chain};
+ }
+
+ /// Eq. A21 from Lafitte
+ template<typename RhoType>
+ auto get_KHS(const RhoType& pf) const {
+ return forceeval(pow(1.0-pf,4)/(1.0 + 4.0*pf + 4.0*pf*pf - 4.0*pf*pf*pf + pf*pf*pf*pf));
+ }
+
+ /**
+ \f[
+ \rho_s\frac{\partial K_{HS}}{\partial \rho_s} = \xi\frac{\partial K_{HS}}{\partial \xi}
+ \f]
+ */
+ template<typename RhoType>
+ auto get_rhos_dK_HS_drhos(const RhoType& xi_x) const {
+ auto num = -4.0*POW3(xi_x - 1.0)*(POW2(xi_x) - 5.0*xi_x - 2.0);
+ auto den = POW2(POW4(xi_x) - 4.0*POW3(xi_x) + 4.0*POW2(xi_x) + 4.0*xi_x + 1.0);
+ return forceeval(num/den*xi_x);
+ }
+
+ /// Eq. A6 from Lafitte
+ template<typename RhoType, typename XiType>
+ auto a_HS(const RhoType& rhos, const Eigen::Array<XiType, 4, 1>& xi) const{
+ constexpr double MY_PI = static_cast<double>(EIGEN_PI);
+ return forceeval(6.0/(MY_PI*rhos)*(3.0*xi[1]*xi[2]/(1.0-xi[3]) + POW3(xi[2])/(xi[3]*POW2(1.0-xi[3])) + (POW3(xi[2])/POW2(xi[3])-xi[0])*log(1.0-xi[3])));
+ }
+
+ /**
+ \note Starting from Eq. A12 from Lafitte
+
+ Defining:
+ \f[
+ \hat B_{ij} \equiv \frac{B_{ij}}{2\pi\epsilon_{ij}d^3_{ij}\rho_s} = \frac{1-\xi_x/2}{(1-\xi_x)^3}I-\frac{9\xi_x(1+\xi_x)}{2(1-\xi_x)^3}J
+ \f]
+ */
+ template<typename XiType, typename IJ>
+ auto get_Bhatij(const XiType& xi_x, const XiType& one_minus_xi_x3, const IJ& I, const IJ& J) const{
+ return forceeval(
+ (1.0-xi_x/2.0)/one_minus_xi_x3*I - 9.0*xi_x*(1.0+xi_x)/(2.0*one_minus_xi_x3)*J
+ );
+ }
+
+ /**
+ \f[
+ B = \hat B_{ij}\kappa \rho_s
+ \f]
+ \f[
+ \left(\frac{\partial B_{ij}}{\partial \rho_s}\right)_{T,\vec{z}} = \kappa\left(\hat B + \xi_x \frac{\partial \hat B}{\partial \xi_x}\right)
+ \f]
+ and thus
+ \f[
+ \rho_s \left(\frac{\partial B_{ij}}{\partial \rho_s}\right)_{T,\vec{z}} = \hat B + \xi_x \frac{\partial \hat B}{\partial \xi_x}
+ \f]
+ */
+ template<typename XiType, typename IJ>
+ auto get_rhodBijdrho(const XiType& xi_x, const XiType& one_minus_xi_x3, const IJ& I, const IJ& J, const XiType& Bhatij) const{
+ auto dBhatdxix = (-3.0*I*(xi_x - 2.0) - 27.0*J*xi_x*(xi_x + 1.0) + (xi_x - 1.0)*(I + 9.0*J*xi_x + 9.0*J*(xi_x + 1.0)))/(2.0*POW4(1.0-xi_x));
+ return forceeval(Bhatij + dBhatdxix*xi_x);
+ }
+
+ /**
+ \note Starting from Eq. A16 from Lafitte
+
+ \f[
+ \hat a^S_{1,ii} = \frac{a^S_{1,ii}}{2\pi\epsilon_{ij}d^3_{ij}\rho_s}
+ \f]
+ so
+ \f[
+ a^S_{1,ii} = \kappa\rho_s\hat a^S_{1,ii}
+ \f]
+ */
+ template<typename XiType>
+ auto get_a1Shatij(const XiType& xi_x_eff, double lambda_ij) const{
+ return forceeval(
+ -1.0/(lambda_ij-3.0)*(1.0-xi_x_eff/2.0)/POW3(forceeval(1.0-xi_x_eff))
+ );
+ }
+
+ /**
+ \f[
+ \left(\frac{\partial a^S_{1,ii}}{\partial \rho_s}\right)_{T,\vec{z}} = \kappa\left(\hat a^S_{1,ii} + \rho_s\frac{\partial \hat a^S_{1,ii}}{\partial \rho_s} \right)
+ \f]
+
+ \f[
+ \left(\frac{\partial a^S_{1,ii}}{\partial \rho_s}\right)_{T,\vec{z}} = \kappa\left(\hat a^S_{1,ii} + \rho_s\frac{\partial \hat a^S_{1,ii}}{\partial \xi_{x,eff}}\frac{\partial \xi_{x,eff}}{\partial \xi_x}\frac{\partial \xi_x}{\partial \rho_s} \right)
+ \f]
+
+ since \f$\rho_s\frac{\partial \xi_x}{\partial \rho_s} = \xi_x\f$
+ */
+ template<typename XiType>
+ auto get_rhoda1Shatijdrho(const XiType& xi_x, const XiType& xi_x_eff, const XiType& dxixeffdxix, double lambda_ij) const{
+ auto xixda1Shatdxix = ((2.0*xi_x_eff - 5.0)*dxixeffdxix)/(2.0*(lambda_ij-3)*POW4(xi_x_eff-1.0))*xi_x;
+ return forceeval(get_a1Shatij(xi_x_eff, lambda_ij) + xixda1Shatdxix);
+ }
+};
+
+/**
+ \brief A class used to evaluate mixtures using the SAFT-VR-Mie model
+*/
+class SAFTVRMieMixture {
+private:
+
+ std::vector<std::string> names;
+ const SAFTVRMieChainContributionTerms terms;
+
+ void check_kmat(const Eigen::ArrayXXd& kmat, std::size_t N) {
+ if (kmat.size() == 0){
+ return;
+ }
+ if (kmat.cols() != kmat.rows()) {
+ throw teqp::InvalidArgument("kmat rows and columns are not identical");
+ }
+ if (kmat.cols() != N) {
+ throw teqp::InvalidArgument("kmat needs to be a square matrix the same size as the number of components");
+ }
+ };
+ auto get_coeffs_from_names(const std::vector<std::string> &names){
+ SAFTVRMieLibrary library;
+ return library.get_coeffs(names);
+ }
+ auto build_chain(const std::vector<SAFTVRMieCoeffs> &coeffs, const std::optional<Eigen::ArrayXXd>& kmat){
+ if (kmat){
+ check_kmat(kmat.value(), coeffs.size());
+ }
+ const std::size_t N = coeffs.size();
+ Eigen::ArrayXd m(N), epsilon_over_k(N), sigma_m(N), lambda_r(N), lambda_a(N);
+ auto i = 0;
+ for (const auto &coeff : coeffs) {
+ m[i] = coeff.m;
+ epsilon_over_k[i] = coeff.epsilon_over_k;
+ sigma_m[i] = coeff.sigma_m;
+ lambda_r[i] = coeff.lambda_r;
+ lambda_a[i] = coeff.lambda_a;
+ i++;
+ }
+ if (kmat){
+ return SAFTVRMieChainContributionTerms(m, epsilon_over_k, sigma_m, lambda_r, lambda_a, std::move(kmat.value()));
+ }
+ else{
+ auto mat = Eigen::ArrayXXd::Zero(N,N);
+ return SAFTVRMieChainContributionTerms(m, epsilon_over_k, sigma_m, lambda_r, lambda_a, std::move(mat));
+ }
+ }
+public:
+ SAFTVRMieMixture(const std::vector<std::string> &names, const std::optional<Eigen::ArrayXXd>& kmat = std::nullopt) : SAFTVRMieMixture(get_coeffs_from_names(names), kmat){};
+ SAFTVRMieMixture(const std::vector<SAFTVRMieCoeffs> &coeffs, const std::optional<Eigen::ArrayXXd> &kmat = std::nullopt) : terms(build_chain(coeffs, kmat)) {};
+
+// PCSAFTMixture( const PCSAFTMixture& ) = delete; // non construction-copyable
+ SAFTVRMieMixture& operator=( const SAFTVRMieMixture& ) = delete; // non copyable
+
+ const auto& get_terms() const { return terms; }
+
+ auto get_m() const { return terms.m; }
+ auto get_sigma_Angstrom() const { return (terms.sigma_A).eval(); }
+ auto get_sigma_m() const { return terms.sigma_A/1e10; }
+ auto get_epsilon_over_k_K() const { return terms.epsilon_over_k; }
+ auto get_kmat() const { return terms.kmat; }
+ auto get_lambda_r() const { return terms.lambda_r; }
+ auto get_lambda_a() const { return terms.lambda_a; }
+
+ // template<typename VecType>
+ // double max_rhoN(const double T, const VecType& mole_fractions) const {
+ // auto N = mole_fractions.size();
+ // Eigen::ArrayX<double> d(N);
+ // for (auto i = 0; i < N; ++i) {
+ // d[i] = sigma_Angstrom[i] * (1.0 - 0.12 * exp(-3.0 * epsilon_over_k[i] / T));
+ // }
+ // return 6 * 0.74 / EIGEN_PI / (mole_fractions*m*powvec(d, 3)).sum()*1e30; // particles/m^3
+ // }
+
+ template<class VecType>
+ auto R(const VecType& molefrac) const {
+ return get_R_gas<decltype(molefrac[0])>();
+ }
+
+ template<typename TTYPE, typename RhoType, typename VecType>
+ auto alphar(const TTYPE& T, const RhoType& rhomolar, const VecType& mole_fractions) const {
+ // First values for the Mie chain with dispersion (always included)
+ error_if_expr(T); error_if_expr(rhomolar);
+ auto vals = terms.get_core_calcs(T, rhomolar, mole_fractions);
+ auto alphar = vals.alphar_mono + vals.alphar_chain;
+ return forceeval(alphar);
+ }
+};
+
+} /* namespace SAFTVRMie */
+}; // namespace teqp
diff --git a/include/teqp/types.hpp b/include/teqp/types.hpp
index 2f10256313d3c0cabf1a69ef80cd92ceadd4e3d0..4f66d0ddccd958e5e04d963c754f31a4e0e485d3 100644
--- a/include/teqp/types.hpp
+++ b/include/teqp/types.hpp
@@ -54,6 +54,17 @@ namespace teqp {
return expr;
}
}
+
+ /// A constexpr function for ensuring that an argument to a function is NOT an expr,
+ /// which can have surprising behavior
+ template<typename T>
+ void error_if_expr(T&& expr)
+ {
+ using namespace autodiff::detail;
+ if constexpr (isExpr<T>) {
+ static_assert(true, "Argument to function is an expression, but should not be");
+ }
+ }
// See https://stackoverflow.com/a/41438758
template<typename T> struct is_complex_t : public std::false_type {};
diff --git a/interface/CPP/teqp_impl_Arxy.cpp b/interface/CPP/teqp_impl_Arxy.cpp
index ff2fcd2227831485545bb229df128bf1123d85d8..2326d6a6527eb1ce423aa8d4ed579edf8a59d966 100644
--- a/interface/CPP/teqp_impl_Arxy.cpp
+++ b/interface/CPP/teqp_impl_Arxy.cpp
@@ -5,9 +5,9 @@ using MI = teqp::cppinterface::ModelImplementer;
// Derivatives from isochoric thermodynamics (all have the same signature)
#define X(f) \
-double MI::f(const double T, const REArrayd& rhovec) const { \
+double MI::f(const double T, const EArrayd& rhovec) const { \
return std::visit([&](const auto& model) { \
- using id = IsochoricDerivatives<decltype(model), double, REArrayd>; \
+ using id = IsochoricDerivatives<decltype(model), double, EArrayd>; \
return id::f(model, T, rhovec); \
}, m_model); \
}
@@ -15,9 +15,9 @@ ISOCHORIC_double_args
#undef X
#define X(f) \
-EArrayd MI::f(const double T, const REArrayd& rhovec) const { \
+EArrayd MI::f(const double T, const EArrayd& rhovec) const { \
return std::visit([&](const auto& model) { \
- using id = IsochoricDerivatives<decltype(model), double, REArrayd>; \
+ using id = IsochoricDerivatives<decltype(model), double, EArrayd>; \
return id::f(model, T, rhovec); \
}, m_model); \
}
@@ -25,9 +25,9 @@ ISOCHORIC_array_args
#undef X
#define X(f) \
-EMatrixd MI::f(const double T, const REArrayd& rhovec) const { \
+EMatrixd MI::f(const double T, const EArrayd& rhovec) const { \
return std::visit([&](const auto& model) { \
- using id = IsochoricDerivatives<decltype(model), double, REArrayd>; \
+ using id = IsochoricDerivatives<decltype(model), double, EArrayd>; \
return id::f(model, T, rhovec); \
}, m_model); \
}
diff --git a/interface/CPP/teqpcpp.cpp b/interface/CPP/teqpcpp.cpp
index 11e268b7a293a6fa7b5f7846f9518beab035baaa..8bcd676637595c5c889d7dd19a164c03446538ad 100644
--- a/interface/CPP/teqpcpp.cpp
+++ b/interface/CPP/teqpcpp.cpp
@@ -48,17 +48,17 @@ namespace teqp {
// Derivatives from isochoric thermodynamics (all have the same signature)
#define X(f) \
- double f(const double T, const REArrayd& rhovec) const override ;
+ double f(const double T, const EArrayd& rhovec) const override ;
ISOCHORIC_double_args
#undef X
#define X(f) \
- EArrayd f(const double T, const REArrayd& rhovec) const override ;
+ EArrayd f(const double T, const EArrayd& rhovec) const override ;
ISOCHORIC_array_args
#undef X
#define X(f) \
- EMatrixd f(const double T, const REArrayd& rhovec) const override ;
+ EMatrixd f(const double T, const EArrayd& rhovec) const override ;
ISOCHORIC_matrix_args
#undef X
diff --git a/interface/pybind11_wrapper.cpp b/interface/pybind11_wrapper.cpp
index 40813f3ec9320c1cc04655561bd7e38f91831610..15da6cb429fb988f6750f4ff8d96bc6da284185b 100644
--- a/interface/pybind11_wrapper.cpp
+++ b/interface/pybind11_wrapper.cpp
@@ -104,6 +104,16 @@ void attach_model_specific_methods(py::object& obj){
setattr("max_rhoN", MethodType(py::cpp_function([](py::object& o, double T, REArrayd& molefrac){ return get_typed<PCSAFT_t>(o).max_rhoN(T, molefrac); }, "self"_a, "T"_a, "molefrac"_a), obj));
setattr("get_kmat", MethodType(py::cpp_function([](py::object& o){ return get_typed<PCSAFT_t>(o).get_kmat(); }), obj));
}
+ else if (std::holds_alternative<SAFTVRMie_t>(model)){
+ setattr("get_m", MethodType(py::cpp_function([](py::object& o){ return get_typed<SAFTVRMie_t>(o).get_m(); }), obj));
+ setattr("get_sigma_Angstrom", MethodType(py::cpp_function([](py::object& o){ return get_typed<SAFTVRMie_t>(o).get_sigma_Angstrom(); }), obj));
+ setattr("get_sigma_m", MethodType(py::cpp_function([](py::object& o){ return get_typed<SAFTVRMie_t>(o).get_sigma_m(); }), obj));
+ setattr("get_epsilon_over_k_K", MethodType(py::cpp_function([](py::object& o){ return get_typed<SAFTVRMie_t>(o).get_epsilon_over_k_K(); }), obj));
+ setattr("get_lambda_r", MethodType(py::cpp_function([](py::object& o){ return get_typed<SAFTVRMie_t>(o).get_lambda_r(); }), obj));
+ setattr("get_lambda_a", MethodType(py::cpp_function([](py::object& o){ return get_typed<SAFTVRMie_t>(o).get_lambda_a(); }), obj));
+// setattr("max_rhoN", MethodType(py::cpp_function([](py::object& o, double T, REArrayd& molefrac){ return get_typed<SAFTVRMie_t>(o).max_rhoN(T, molefrac); }, "self"_a, "T"_a, "molefrac"_a), obj));
+ setattr("get_kmat", MethodType(py::cpp_function([](py::object& o){ return get_typed<SAFTVRMie_t>(o).get_kmat(); }), obj));
+ }
else if (std::holds_alternative<canonical_cubic_t>(model)){
setattr("get_a", MethodType(py::cpp_function([](py::object& o, double T, REArrayd& molefrac){ return get_typed<canonical_cubic_t>(o).get_a(T, molefrac); }, "self"_a, "T"_a, "molefrac"_a), obj));
setattr("get_b", MethodType(py::cpp_function([](py::object& o, double T, REArrayd& molefrac){ return get_typed<canonical_cubic_t>(o).get_b(T, molefrac); }, "self"_a, "T"_a, "molefrac"_a), obj));
diff --git a/src/tests/catch_test_LJmodels.cxx b/src/tests/catch_test_LJmodels.cxx
index b0604ec1a9592f973b94ce55b62be9d223119e85..116fc61364a83a2a0d9e3ae748764b500ea46bed 100644
--- a/src/tests/catch_test_LJmodels.cxx
+++ b/src/tests/catch_test_LJmodels.cxx
@@ -21,10 +21,23 @@ TEST_CASE("Test for critical point of Kolafa-Nezbeda", "[LJ126]")
TEST_CASE("Test for critical point of Thol", "[LJ126-Thol]")
{
auto m = build_LJ126_TholJPCRD2016();
- auto soln = solve_pure_critical(m, 1.3, 0.3);
+
auto Tc = 1.32, rhoc = 0.31;
- CHECK(std::get<0>(soln) == Approx(Tc).margin(0.01));
+ auto z = (Eigen::ArrayXd(1) << 1.0).finished();
+ using tdx = TDXDerivatives<decltype(m)>;
+ // Generated with ljfeos.cpp file included with paper
+ double tol = 1e-100;
+ CHECK(tdx::get_Arxy<0,0>(m, Tc, rhoc, z) == Approx(-0.848655028421735).margin(tol));
+ CHECK(tdx::get_Arxy<1,0>(m, Tc, rhoc, z) == Approx(-1.72121486947745 ).margin(tol));
+ CHECK(tdx::get_Arxy<0,1>(m, Tc, rhoc, z) == Approx(-0.682159784922984).margin(tol));
+ CHECK(tdx::get_Arxy<2,0>(m, Tc, rhoc, z) == Approx(-3.06580477714764).margin(tol));
+ CHECK(tdx::get_Arxy<1,1>(m, Tc, rhoc, z) == Approx(-1.43011286241234).margin(tol));
+ CHECK(tdx::get_Arxy<0,2>(m, Tc, rhoc, z) == Approx(0.364319525724386).margin(tol));
+
+ auto soln = solve_pure_critical(m, 1.3, 0.3);
+ CHECK(std::get<0>(soln) == Approx(Tc).margin(0.2));
CHECK(std::get<1>(soln) == Approx(rhoc).margin(0.01));
+
}
TEST_CASE("Test for critical point of Johnson", "[LJ126]")
{
diff --git a/src/tests/catch_test_SAFTVRMie.cxx b/src/tests/catch_test_SAFTVRMie.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..aa7d9ec9cc37c622531c235668d63997e3e354a4
--- /dev/null
+++ b/src/tests/catch_test_SAFTVRMie.cxx
@@ -0,0 +1,240 @@
+#include <catch2/catch_test_macros.hpp>
+#include <catch2/catch_approx.hpp>
+
+using Catch::Approx;
+
+#include <fstream>
+
+#include "teqp/core.hpp"
+#include "teqp/derivs.hpp"
+#include "teqp/models/saftvrmie.hpp"
+#include "teqp/models/pcsaft.hpp"
+#include "teqp/math/finite_derivs.hpp"
+#include "teqp/json_builder.hpp"
+#include "teqp/cpp/teqpcpp.hpp"
+#include "teqp/constants.hpp"
+#include "teqp/algorithms/critical_pure.hpp"
+#include "teqp/algorithms/VLE.hpp"
+using namespace teqp::SAFTVRMie;
+using namespace teqp;
+
+#include "boost/multiprecision/cpp_bin_float.hpp"
+#include "boost/multiprecision/cpp_complex.hpp"
+#include "teqp/math/quadrature.hpp"
+
+TEST_CASE("Check integration", "[SAFTVRMIE]"){
+ std::function<double(double)> f = [](const double&x){ return x*sin(x); };
+ auto exact = -2*cos(1) + 2*sin(1);
+ auto deg3 = quad<3, double>(f, -1, 1);
+ auto deg4 = quad<4, double>(f, -1, 1);
+ auto deg5 = quad<5, double>(f, -1, 1);
+ CHECK(deg4 == Approx(exact).margin(1e-12));
+ CHECK(deg5 == Approx(exact).margin(1e-12));
+}
+
+TEST_CASE("Check integration for d", "[SAFTVRMIE]"){
+ double epskB = 206.12;
+ double sigma_m = 3.7257e-10;
+ double lambda_r = 12.4;
+ double lambda_a = 6.0;
+ double C = lambda_r/(lambda_r-lambda_a)*::pow(lambda_r/lambda_a,lambda_a/(lambda_r-lambda_a));
+ double T = 300.0;
+ std::function<double(double)> integrand = [&epskB, &sigma_m, &C, &T, &lambda_a, &lambda_r](const double& r_m){
+ auto u = C*epskB*(::pow(sigma_m/r_m, lambda_r) - ::pow(sigma_m/r_m, lambda_a));
+ return -expm1(-u/T);
+ };
+ auto d30 = quad<30, double>(integrand, 0.0, sigma_m);
+ auto d15 = quad<15, double>(integrand, 0.0, sigma_m);
+ auto d7 = quad<7, double>(integrand, 0.0, sigma_m);
+ auto d5 = quad<5, double>(integrand, 0.0, sigma_m);
+ auto d3 = quad<3, double>(integrand, 0.0, sigma_m);
+ CHECK(d30 == Approx(3.597838581533227e-10).margin(1e-12));
+}
+
+TEST_CASE("Single alphar check value", "[SAFTVRMie]")
+{
+ auto m = (Eigen::ArrayXd(1) << 1.4373).finished();
+ auto eps = (Eigen::ArrayXd(1) << 206.12).finished();
+ auto sigma = (Eigen::ArrayXd(1) << 3.7257e-10).finished();
+ auto lambda_r = (Eigen::ArrayXd(1) << 12.4).finished();
+ auto lambda_a = (Eigen::ArrayXd(1) << 6.0).finished();
+ Eigen::ArrayXXd kmat = Eigen::ArrayXXd::Zero(1,1);
+ SAFTVRMieChainContributionTerms terms(m, eps, sigma, lambda_r, lambda_a, kmat);
+
+ auto z = (Eigen::ArrayXd(1) << 1.0).finished();
+ auto core = terms.get_core_calcs(300.0, 12000.0, z);
+
+ // Check values from Clapeyron.jl
+ CHECK(core.a1kB == Approx(-694.2608061145818));
+ CHECK(core.a2kB2 == Approx(-6741.101051705587));
+ CHECK(core.a3kB3 == Approx(-81372.77460911816));
+ CHECK(core.alphar_mono == Approx(-1.322739797471788));
+ CHECK(core.alphar_chain == Approx(-0.0950261207746853));
+}
+
+template<int i, int j, typename Model, typename TTYPE, typename RhoType, typename VecType>
+auto ijcheck(const Model& model, const TTYPE& T, const RhoType& rho, const VecType& z, double margin=1e-103){
+ using tdx = TDXDerivatives<decltype(model)>;
+ auto Arxy = tdx::template get_Arxy<i, j, ADBackends::autodiff>(model, T, rho, z);
+ auto Arxymcx = tdx::template get_Arxy<i, j, ADBackends::multicomplex>(model, T, rho, z);
+ CAPTURE(i);
+ CAPTURE(j);
+ CHECK(Arxymcx == Approx(Arxy).margin(margin));
+ CHECK(std::isfinite(Arxymcx));
+}
+
+TEST_CASE("Check all xy derivs", "[SAFTVRMie]")
+{
+ Eigen::ArrayXXd kmat = Eigen::ArrayXXd::Zero(1,1);
+ std::vector<std::string> names = {"Ethane"};
+
+ SAFTVRMieMixture model{names, kmat};
+
+ double T = 300.0, rho = 10000.0;
+ auto z = (Eigen::ArrayXd(1) << 1.0).finished();
+
+ ijcheck<0,0>(model, T, rho, z);
+ ijcheck<0,1>(model, T, rho, z);
+ ijcheck<0,2>(model, T, rho, z);
+ ijcheck<0,3>(model, T, rho, z);
+ ijcheck<1,0>(model, T, rho, z);
+ ijcheck<1,1>(model, T, rho, z);
+ ijcheck<1,2>(model, T, rho, z);
+ ijcheck<1,3>(model, T, rho, z);
+ ijcheck<2,0>(model, T, rho, z);
+ ijcheck<2,1>(model, T, rho, z);
+ ijcheck<2,2>(model, T, rho, z);
+ ijcheck<2,3>(model, T, rho, z);
+ int rr = 0;
+}
+
+TEST_CASE("Solve for critical point with three interface approaches", "[SAFTVRMie]")
+{
+ Eigen::ArrayXXd kmat = Eigen::ArrayXXd::Zero(1,1);
+ std::vector<std::string> names = {"Ethane"};
+ SAFTVRMieMixture model_{names, kmat};
+ double T = 300.0, rho = 10000.0;
+ auto z = (Eigen::ArrayXd(1) << 1.0).finished();
+ auto crit1 = solve_pure_critical(model_, 300.0, 10000.0);
+
+ nlohmann::json jcoeffs = nlohmann::json::array();
+ jcoeffs.push_back({
+ {"name", "Ethane"},
+ { "m", 1.4373 },
+ { "sigma_m", 3.7257e-10},
+ {"epsilon_over_k", 206.12},
+ {"lambda_r", 12.4},
+ {"lambda_a", 6.0},
+ {"BibTeXKey", "Lafitte-JCP"}
+ });
+ nlohmann::json model = {
+ {"coeffs", jcoeffs}
+ };
+ nlohmann::json j = {
+ {"kind", "SAFT-VR-Mie"},
+ {"model", model}
+ };
+ auto modelc = cppinterface::make_model(j);
+ auto crit2 = modelc->solve_pure_critical(300.0, 10000.0, {});
+ CHECK(std::get<0>(crit1) == Approx(std::get<0>(crit2)));
+
+ nlohmann::json model2 = {
+ {"names", {"Ethane"}}
+ };
+ nlohmann::json j2 = {
+ {"kind", "SAFT-VR-Mie"},
+ {"model", model2}
+ };
+ auto model3 = cppinterface::make_model(j2);
+ auto crit3 = model3->solve_pure_critical(300.0, 10000.0, {});
+ CHECK(std::get<0>(crit1) == Approx(std::get<0>(crit3)));
+
+ auto z3 = (Eigen::ArrayXd(1) << 1.0).finished();
+ CHECK(std::isfinite(model3->get_Ar11(300.0, 300.0, z3)));
+
+ int rr = 0;
+}
+
+TEST_CASE("A mixture calculation", "[SAFTVRMie]"){
+ std::vector<std::string> names = {"Methane","Ethane","Propane"};
+ SAFTVRMieMixture model_{names};
+ double T = 300.0, rho = 1000.0;
+ nlohmann::json model2 = {
+ {"names", names}
+ };
+ nlohmann::json j2 = {
+ {"kind", "PCSAFT"},
+ {"model", model2}
+ };
+ auto model3 = cppinterface::make_model(j2);
+ auto z = (Eigen::ArrayXd(3) << 0.3, 0.4, 0.3).finished();
+ auto rhovec = (rho*z).eval();
+ auto fugcoeff = model3->get_fugacity_coefficients(T, rhovec);
+}
+
+TEST_CASE("Trace critical locus", "[SAFTVRMie]"){
+ std::vector<std::string> names = {"Methane", "Ethane"};
+ nlohmann::json model2 = {
+ {"names", names}
+ };
+ nlohmann::json j2 = {
+ {"kind", "SAFT-VR-Mie"},
+ {"model", model2}
+ };
+ auto model3 = cppinterface::make_model(j2);
+ auto [T0, rho0] = model3->solve_pure_critical(300, 10000, nlohmann::json{{"alternative_pure_index", 0},{"alternative_length", names.size()}});
+ auto rhovec0 = (Eigen::ArrayXd(2) << rho0, 0).finished();
+ auto al = model3->trace_critical_arclength_binary(T0, rhovec0, "aa.txt");
+ int rr = 0;
+}
+
+TEST_CASE("VLE pure tracing", "[SAFTVRMieVLE]"){
+ std::vector<std::string> names1 = {"Ethane"};
+ SAFTVRMieMixture pure{names1};
+ nlohmann::json spec1{
+ {"Tcguess", 300.0},
+ {"rhocguess", 10000.0},
+ {"Tred", 0.999},
+ {"Nstep", 100},
+ {"with_deriv", true}
+ };
+ auto o1 = pure_trace_VLE(pure, 300, spec1);
+
+ std::vector<std::string> names = {"Methane", "Ethane"};
+ SAFTVRMieMixture model{names};
+ nlohmann::json pure_spec{{"alternative_pure_index", 1}, {"alternative_length", names.size()}};
+ nlohmann::json spec{
+ {"Tcguess", 300.0},
+ {"rhocguess", 10000.0},
+ {"pure_spec", pure_spec},
+ {"Tred", 0.999},
+ {"Nstep", 100},
+ {"with_deriv", true}
+ };
+ auto o = pure_trace_VLE(model, 300, spec);
+ CHECK(o[0] == o1[0]);
+}
+
+TEST_CASE("VLE isotherm tracing", "[SAFTVRMieVLE]"){
+
+ std::vector<std::string> names = {"Ethane", "Propane"};
+ SAFTVRMieMixture model{names};
+ nlohmann::json pure_spec{{"alternative_pure_index", 0}, {"alternative_length", names.size()}};
+ nlohmann::json spec{
+ {"Tcguess", 300.0},
+ {"rhocguess", 10000.0},
+ {"pure_spec", pure_spec},
+ {"Tred", 0.999},
+ {"Nstep", 100},
+ {"with_deriv", true}
+ };
+ double T = 280.0;
+ auto purestart = pure_trace_VLE(model, T, spec);
+ Eigen::ArrayXd rhovecL0(2), rhovecV0(2);
+ int ipure = pure_spec.at("alternative_pure_index");
+ rhovecL0(ipure) = purestart[0]; rhovecV0(ipure) = purestart[1];
+ auto der = get_drhovecdp_Tsat(model, T, rhovecL0, rhovecV0);
+ auto opt = TVLEOptions(); opt.revision = 2; opt.polish = true; opt.init_c = -1;
+ auto iso = trace_VLE_isotherm_binary(model, T, rhovecL0, rhovecV0, opt);
+// std::cout << iso.dump(2) << std::endl;
+}