diff --git a/include/teqp/cpp/teqpcpp.hpp b/include/teqp/cpp/teqpcpp.hpp
index 2478c5cc5ec64ab590a1bff229fda48c9384b095..d5f1bce70fa14e5fb99699ababf8918a37ac1f5a 100644
--- a/include/teqp/cpp/teqpcpp.hpp
+++ b/include/teqp/cpp/teqpcpp.hpp
@@ -176,8 +176,6 @@ namespace teqp {
const std::string& departurepath = {}
);
- std::unique_ptr<AbstractModel> emplace_model(AllowedModels&& model);
-
std::unique_ptr<AbstractModel> build_model_ptr(const nlohmann::json& json);
}
}
diff --git a/include/teqp/json_builder.hpp b/include/teqp/json_builder.hpp
index 23502cced6fc7dd635918389abf398b5bb5a977c..0b2d73f3cda6509d94a26f846d33719d738980fa 100644
--- a/include/teqp/json_builder.hpp
+++ b/include/teqp/json_builder.hpp
@@ -44,46 +44,7 @@ namespace teqp {
return CPA::CPAfactory(spec);
}
else if (kind == "PCSAFT") {
- using namespace PCSAFT;
- std::optional<Eigen::ArrayXXd> kmat;
- if (spec.contains("kmat") && spec.at("kmat").is_array() && spec.at("kmat").size() > 0){
- kmat = build_square_matrix(spec["kmat"]);
- }
-
- if (spec.contains("names")){
- std::vector<std::string> names = spec["names"];
- if (kmat && kmat.value().rows() != names.size()){
- throw teqp::InvalidArgument("Provided length of names of " + std::to_string(names.size()) + " does not match the dimension of the kmat of " + std::to_string(kmat.value().rows()));
- }
- return PCSAFTMixture(names, kmat.value_or(Eigen::ArrayXXd{}));
- }
- else if (spec.contains("coeffs")){
- std::vector<SAFTCoeffs> coeffs;
- for (auto j : spec["coeffs"]) {
- SAFTCoeffs c;
- c.name = j.at("name");
- c.m = j.at("m");
- c.sigma_Angstrom = j.at("sigma_Angstrom");
- c.epsilon_over_k = j.at("epsilon_over_k");
- c.BibTeXKey = j.at("BibTeXKey");
- if (j.contains("(mu^*)^2") && j.contains("nmu")){
- c.mustar2 = j.at("(mu^*)^2");
- c.nmu = j.at("nmu");
- }
- if (j.contains("(Q^*)^2") && j.contains("nQ")){
- c.Qstar2 = j.at("(Q^*)^2");
- c.nQ = j.at("nQ");
- }
- coeffs.push_back(c);
- }
- if (kmat && kmat.value().rows() != coeffs.size()){
- throw teqp::InvalidArgument("Provided length of coeffs of " + std::to_string(coeffs.size()) + " does not match the dimension of the kmat of " + std::to_string(kmat.value().rows()));
- }
- return PCSAFTMixture(coeffs, kmat.value_or(Eigen::ArrayXXd{}));
- }
- else{
- throw std::invalid_argument("you must provide names or coeffs, but not both");
- }
+ return PCSAFT::PCSAFTfactory(spec);
}
else if (kind == "SAFT-VR-Mie") {
return SAFTVRMie::SAFTVRMiefactory(spec);
diff --git a/include/teqp/models/cubics.hpp b/include/teqp/models/cubics.hpp
index 6d8dd256b940fb823097e4c681e4fc2b1724f56a..cd1e16d398af4ed5c3b3214e5cc8c743aac7c91c 100644
--- a/include/teqp/models/cubics.hpp
+++ b/include/teqp/models/cubics.hpp
@@ -12,6 +12,7 @@ Implementations of the canonical cubic equations of state
#include "teqp/constants.hpp"
#include "teqp/exceptions.hpp"
#include "cubicsuperancillary.hpp"
+#include "teqp/json_tools.hpp"
#include "nlohmann/json.hpp"
@@ -20,16 +21,16 @@ Implementations of the canonical cubic equations of state
namespace teqp {
/**
-* \brief The standard alpha function used by Peng-Robinson and SRK
-*/
+ * \brief The standard alpha function used by Peng-Robinson and SRK
+ */
template<typename NumType>
class BasicAlphaFunction {
private:
NumType Tci, ///< The critical temperature
- mi; ///< The "m" parameter
+ mi; ///< The "m" parameter
public:
BasicAlphaFunction(NumType Tci, NumType mi) : Tci(Tci), mi(mi) {};
-
+
template<typename TType>
auto operator () (const TType& T) const {
return forceeval(pow2(forceeval(1.0 + mi * (1.0 - sqrt(T / Tci)))));
@@ -44,18 +45,18 @@ class GenericCubic {
protected:
std::valarray<NumType> ai, bi;
const NumType Delta1, Delta2, OmegaA, OmegaB;
- int superanc_index;
+ int superanc_index;
const AlphaFunctions alphas;
Eigen::ArrayXXd kmat;
-
+
nlohmann::json meta;
-
+
template<typename TType, typename IndexType>
auto get_ai(TType T, IndexType i) const { return ai[i]; }
-
+
template<typename TType, typename IndexType>
auto get_bi(TType T, IndexType i) const { return bi[i]; }
-
+
template<typename IndexType>
void check_kmat(IndexType N) {
if (kmat.cols() != kmat.rows()) {
@@ -68,10 +69,10 @@ protected:
throw teqp::InvalidArgument("kmat needs to be a square matrix the same size as the number of components [" + std::to_string(N) + "]");
}
};
-
+
public:
GenericCubic(NumType Delta1, NumType Delta2, NumType OmegaA, NumType OmegaB, int superanc_index, const std::valarray<NumType>& Tc_K, const std::valarray<NumType>& pc_Pa, const AlphaFunctions& alphas, const Eigen::ArrayXXd& kmat)
- : Delta1(Delta1), Delta2(Delta2), OmegaA(OmegaA), OmegaB(OmegaB), superanc_index(superanc_index), alphas(alphas), kmat(kmat)
+ : Delta1(Delta1), Delta2(Delta2), OmegaA(OmegaA), OmegaB(OmegaB), superanc_index(superanc_index), alphas(alphas), kmat(kmat)
{
ai.resize(Tc_K.size());
bi.resize(Tc_K.size());
@@ -81,13 +82,13 @@ public:
}
check_kmat(ai.size());
};
-
+
void set_meta(const nlohmann::json& j) { meta = j; }
auto get_meta() const { return meta; }
auto get_kmat() const { return kmat; }
-
+
/// Return a tuple of saturated liquid and vapor densities for the EOS given the temperature
- /// Uses the superancillary equations from Bell and Deiters:
+ /// Uses the superancillary equations from Bell and Deiters:
auto superanc_rhoLV(double T) const {
if (ai.size() != 1) {
throw std::invalid_argument("function only available for pure species");
@@ -96,18 +97,18 @@ public:
auto b = get_b(T, z);
auto Ttilde = R(z)*T*b/get_a(T,z);
return std::make_tuple(
- CubicSuperAncillary::supercubic(superanc_index, CubicSuperAncillary::RHOL_CODE, Ttilde)/b,
- CubicSuperAncillary::supercubic(superanc_index, CubicSuperAncillary::RHOV_CODE, Ttilde)/b
- );
+ CubicSuperAncillary::supercubic(superanc_index, CubicSuperAncillary::RHOL_CODE, Ttilde)/b,
+ CubicSuperAncillary::supercubic(superanc_index, CubicSuperAncillary::RHOV_CODE, Ttilde)/b
+ );
}
-
+
const NumType Ru = get_R_gas<double>(); /// Universal gas constant, exact number
-
+
template<class VecType>
auto R(const VecType& molefrac) const {
return Ru;
}
-
+
template<typename TType, typename CompType>
auto get_a(TType T, const CompType& molefracs) const {
std::common_type_t<TType, decltype(molefracs[0])> a_ = 0.0;
@@ -124,7 +125,7 @@ public:
}
return forceeval(a_);
}
-
+
template<typename TType, typename CompType>
auto get_b(TType /*T*/, const CompType& molefracs) const {
std::common_type_t<TType, decltype(molefracs[0])> b_ = 0.0;
@@ -133,11 +134,11 @@ public:
}
return forceeval(b_);
}
-
+
template<typename TType, typename RhoType, typename MoleFracType>
auto alphar(const TType& T,
- const RhoType& rho,
- const MoleFracType& molefrac) const
+ const RhoType& rho,
+ const MoleFracType& molefrac) const
{
if (molefrac.size() != alphas.size()) {
throw std::invalid_argument("Sizes do not match");
@@ -155,16 +156,16 @@ auto canonical_SRK(TCType Tc_K, PCType pc_K, AcentricType acentric, const Eigen:
double Delta1 = 1;
double Delta2 = 0;
AcentricType m = 0.48 + 1.574 * acentric - 0.176 * acentric * acentric;
-
+
std::vector<AlphaFunctionOptions> alphas;
for (auto i = 0; i < Tc_K.size(); ++i) {
alphas.emplace_back(BasicAlphaFunction(Tc_K[i], m[i]));
}
-
+
// See https://doi.org/10.1021/acs.iecr.1c00847
double OmegaA = 1.0 / (9.0 * (cbrt(2) - 1));
double OmegaB = (cbrt(2) - 1) / 3;
-
+
nlohmann::json meta = {
{"Delta1", Delta1},
{"Delta2", Delta2},
@@ -172,12 +173,22 @@ auto canonical_SRK(TCType Tc_K, PCType pc_K, AcentricType acentric, const Eigen:
{"OmegaB", OmegaB},
{"kind", "Soave-Redlich-Kwong"}
};
-
+
auto cub = GenericCubic(Delta1, Delta2, OmegaA, OmegaB, CubicSuperAncillary::SRK_CODE, Tc_K, pc_K, std::move(alphas), kmat);
cub.set_meta(meta);
return cub;
}
+/// A JSON-based factory function for the canonical SRK model
+inline auto make_canonicalSRK(const nlohmann::json& spec){
+ std::valarray<double> Tc_K = spec.at("Tcrit / K"), pc_Pa = spec.at("pcrit / Pa"), acentric = spec.at("acentric");
+ Eigen::ArrayXXd kmat(0, 0);
+ if (spec.contains("kmat")){
+ kmat = build_square_matrix(spec.at("kmat"));
+ }
+ return canonical_SRK(Tc_K, pc_Pa, acentric, kmat);
+}
+
template <typename TCType, typename PCType, typename AcentricType>
auto canonical_PR(TCType Tc_K, PCType pc_K, AcentricType acentric, const Eigen::ArrayXXd& kmat = {}) {
double Delta1 = 1+sqrt(2.0);
@@ -211,4 +222,14 @@ auto canonical_PR(TCType Tc_K, PCType pc_K, AcentricType acentric, const Eigen::
return cub;
}
+/// A JSON-based factory function for the canonical SRK model
+inline auto make_canonicalPR(const nlohmann::json& spec){
+ std::valarray<double> Tc_K = spec.at("Tcrit / K"), pc_Pa = spec.at("pcrit / Pa"), acentric = spec.at("acentric");
+ Eigen::ArrayXXd kmat(0, 0);
+ if (spec.contains("kmat")){
+ kmat = build_square_matrix(spec.at("kmat"));
+ }
+ return canonical_PR(Tc_K, pc_Pa, acentric, kmat);
+}
+
}; // namespace teqp
diff --git a/include/teqp/models/pcsaft.hpp b/include/teqp/models/pcsaft.hpp
index 324c134cea6ef2ba8e7a866bcb7f52674f50a789..04b534b0476045b4273655086b2e64f28b4cc49a 100644
--- a/include/teqp/models/pcsaft.hpp
+++ b/include/teqp/models/pcsaft.hpp
@@ -410,19 +410,48 @@ public:
}
};
-inline auto PCSAFTfactory(const nlohmann::json& json) {
- std::vector<SAFTCoeffs> coeffs;
- for (auto j : json) {
- SAFTCoeffs c;
- c.name = j.at("name");
- c.m = j.at("m");
- c.sigma_Angstrom = j.at("sigma_Angstrom");
- c.epsilon_over_k = j.at("epsilon_over_k");
- c.BibTeXKey = j.at("BibTeXKey");
- coeffs.push_back(c);
+/// A JSON-based factory function for the PC-SAFT model
+inline auto PCSAFTfactory(const nlohmann::json& spec) {
+ std::optional<Eigen::ArrayXXd> kmat;
+ if (spec.contains("kmat") && spec.at("kmat").is_array() && spec.at("kmat").size() > 0){
+ kmat = build_square_matrix(spec["kmat"]);
}
- return PCSAFTMixture(coeffs);
-};
+
+ if (spec.contains("names")){
+ std::vector<std::string> names = spec["names"];
+ if (kmat && kmat.value().rows() != names.size()){
+ throw teqp::InvalidArgument("Provided length of names of " + std::to_string(names.size()) + " does not match the dimension of the kmat of " + std::to_string(kmat.value().rows()));
+ }
+ return PCSAFTMixture(names, kmat.value_or(Eigen::ArrayXXd{}));
+ }
+ else if (spec.contains("coeffs")){
+ std::vector<SAFTCoeffs> coeffs;
+ for (auto j : spec["coeffs"]) {
+ SAFTCoeffs c;
+ c.name = j.at("name");
+ c.m = j.at("m");
+ c.sigma_Angstrom = j.at("sigma_Angstrom");
+ c.epsilon_over_k = j.at("epsilon_over_k");
+ c.BibTeXKey = j.at("BibTeXKey");
+ if (j.contains("(mu^*)^2") && j.contains("nmu")){
+ c.mustar2 = j.at("(mu^*)^2");
+ c.nmu = j.at("nmu");
+ }
+ if (j.contains("(Q^*)^2") && j.contains("nQ")){
+ c.Qstar2 = j.at("(Q^*)^2");
+ c.nQ = j.at("nQ");
+ }
+ coeffs.push_back(c);
+ }
+ if (kmat && kmat.value().rows() != coeffs.size()){
+ throw teqp::InvalidArgument("Provided length of coeffs of " + std::to_string(coeffs.size()) + " does not match the dimension of the kmat of " + std::to_string(kmat.value().rows()));
+ }
+ return PCSAFTMixture(coeffs, kmat.value_or(Eigen::ArrayXXd{}));
+ }
+ else{
+ throw std::invalid_argument("you must provide names or coeffs, but not both");
+ }
+}
} /* namespace PCSAFT */
}; // namespace teqp
diff --git a/interface/CPP/teqp_impl_factory.cpp b/interface/CPP/teqp_impl_factory.cpp
index 1f46fb0123123f2d1bd4552605379c3e99f65b5b..bdecf6f3d38958af8c2786b5f16225af7b4a1c7d 100644
--- a/interface/CPP/teqp_impl_factory.cpp
+++ b/interface/CPP/teqp_impl_factory.cpp
@@ -4,15 +4,6 @@
namespace teqp {
namespace cppinterface {
-
- // std::unique_ptr<AbstractModel> emplace_model(AllowedModels&& model){
- // return make(model);
- // }
-
-#warning PR is not yet implemented
-#warning SRK is not yet implemented
-#warning CPA is not yet implemented
-#warning PCSAFT is not yet implemented
std::unique_ptr<teqp::cppinterface::AbstractModel> make_SAFTVRMie(const nlohmann::json &j);
@@ -21,7 +12,12 @@ namespace teqp {
static std::unordered_map<std::string, makefunc> pointer_factory = {
{"vdW1", [](const nlohmann::json& spec){ return make_owned(vdWEOS1(spec.at("a"), spec.at("b"))); }},
- {"SAFT-VR-Mie", [](const nlohmann::json& spec){ return make_SAFTVRMie(spec); }},
+
+ {"PR", [](const nlohmann::json& spec){ return make_owned(make_canonicalPR(spec));}},
+ {"SRK", [](const nlohmann::json& spec){ return make_owned(make_canonicalSRK(spec));}},
+ {"CPA", [](const nlohmann::json& spec){ return make_owned(CPA::CPAfactory(spec));}},
+ {"PCSAFT", [](const nlohmann::json& spec){ return make_owned(PCSAFT::PCSAFTfactory(spec));}},
+
{"multifluid", [](const nlohmann::json& spec){ return make_owned(multifluidfactory(spec));}},
{"SW_EspindolaHeredia2009", [](const nlohmann::json& spec){ return make_owned(squarewell::EspindolaHeredia2009(spec.at("lambda")));}},
{"EXP6_Kataoka1992", [](const nlohmann::json& spec){ return make_owned(exp6::Kataoka1992(spec.at("alpha")));}},
@@ -32,7 +28,10 @@ namespace teqp {
{"Mie_Pohl2023", [](const nlohmann::json& spec){ return make_owned(Mie::Mie6Pohl2023(spec.at("lambda_a")));}},
{"2CLJF-Dipole", [](const nlohmann::json& spec){ return make_owned(twocenterljf::build_two_center_model_dipole(spec.at("author"), spec.at("L^*"), spec.at("(mu^*)^2")));}},
{"2CLJF-Quadrupole", [](const nlohmann::json& spec){ return make_owned(twocenterljf::build_two_center_model_quadrupole(spec.at("author"), spec.at("L^*"), spec.at("(mu^*)^2")));}},
- {"IdealHelmholtz", [](const nlohmann::json& spec){ return make_owned(IdealHelmholtz(spec));}}
+ {"IdealHelmholtz", [](const nlohmann::json& spec){ return make_owned(IdealHelmholtz(spec));}},
+
+ // Implemented in its own compilation unit to help with compilation time
+ {"SAFT-VR-Mie", [](const nlohmann::json& spec){ return make_SAFTVRMie(spec); }}
};
std::unique_ptr<teqp::cppinterface::AbstractModel> build_model_ptr(const nlohmann::json& json) {