From c532558861d6cdcbc6783ac5d713f200db6a926b Mon Sep 17 00:00:00 2001
From: Ian Bell <ian.bell@nist.gov>
Date: Fri, 7 May 2021 10:29:13 -0400
Subject: [PATCH] Add factory function for CPA
---
include/teqp/models/CPA.hpp | 69 +++++++++++++++++++++++++++++++++----
src/tests/catch_tests.cxx | 30 +++++++++++++---
2 files changed, 88 insertions(+), 11 deletions(-)
diff --git a/include/teqp/models/CPA.hpp b/include/teqp/models/CPA.hpp
index 84fadc5..e9cccfe 100644
--- a/include/teqp/models/CPA.hpp
+++ b/include/teqp/models/CPA.hpp
@@ -1,9 +1,23 @@
#pragma once
+namespace CPA {
+
template<typename X> auto POW2(X x) { return x * x; };
template<typename X> auto POW3(X x) { return x * POW2(x); };
enum class association_classes {not_set, a1A, a2B, a3B, a4C, not_associating};
+
+auto get_association_classes(const std::string& s) {
+ if (s == "1A") { return association_classes::a1A; }
+ else if (s == "2B") { return association_classes::a2B; }
+ else if (s == "2B") { return association_classes::a2B; }
+ else if (s == "3B") { return association_classes::a3B; }
+ else if (s == "4C") { return association_classes::a4C; }
+ else {
+ throw std::invalid_argument("bad association flag:" + s);
+ }
+}
+
enum class radial_dist { CS, KG, OT };
/// Function that calculates the association binding strength between site A of molecule i and site B on molecule j
@@ -97,6 +111,13 @@ auto XA_calc_pure(int N_sites, association_classes scheme, double epsABi, double
};
enum class cubic_flag {not_set, PR, SRK};
+auto get_cubic_flag(const std::string& s) {
+ if (s == "PR") { return cubic_flag::PR; }
+ else if (s == "SRK") { return cubic_flag::SRK; }
+ else {
+ throw std::invalid_argument("bad cubic flag:" + s);
+ }
+}
class CPACubic {
private:
@@ -106,7 +127,7 @@ private:
std::valarray<std::valarray<double>> k_ij;
public:
- CPACubic(cubic_flag flag, const std::valarray<double> a0, const std::valarray<double> bi, const std::valarray<double> c1, const std::valarray<double> Tc) : a0(a0), bi(bi), c1(c1), Tc(Tc) {
+ CPACubic(cubic_flag flag, const std::valarray<double> &a0, const std::valarray<double> &bi, const std::valarray<double> &c1, const std::valarray<double> &Tc) : a0(a0), bi(bi), c1(c1), Tc(Tc) {
switch (flag) {
case cubic_flag::PR:
{ delta_1 = 1 + sqrt(2); delta_2 = 1 - sqrt(2); break; }
@@ -150,7 +171,7 @@ public:
template<typename Cubic>
class CPAAssociation {
private:
- const Cubic& cubic;
+ const Cubic cubic;
const std::vector<association_classes> classes;
const std::vector<int> N_sites;
const std::valarray<double> epsABi, betaABi;
@@ -173,7 +194,7 @@ private:
}
public:
- CPAAssociation(const Cubic &cubic, const std::vector<association_classes>& classes, const std::valarray<double> &epsABi, const std::valarray<double> &betaABi)
+ CPAAssociation(const Cubic &&cubic, const std::vector<association_classes>& classes, const std::valarray<double> &epsABi, const std::valarray<double> &betaABi)
: cubic(cubic), classes(classes), epsABi(epsABi), betaABi(betaABi), N_sites(get_N_sites(classes)) {};
template<typename TType, typename RhoType, typename VecType>
@@ -200,10 +221,10 @@ public:
template <typename Cubic, typename Assoc>
class CPA {
public:
- Cubic cubic;
- Assoc assoc;
+ const Cubic cubic;
+ const Assoc assoc;
- CPA(Cubic cubic, Assoc assoc) : cubic(cubic), assoc(assoc) {
+ CPA(Cubic &&cubic, Assoc &&assoc) : cubic(cubic), assoc(assoc) {
}
/// Residual dimensionless Helmholtz energy from the SRK or PR core and contribution due to association
@@ -219,4 +240,38 @@ public:
return forceeval(alpha_r_cubic + alpha_r_assoc);
}
-};
\ No newline at end of file
+};
+
+/// A factory function to return an instantiated CPA instance given
+/// the JSON representation of the model
+auto CPAfactory(const nlohmann::json &j){
+ auto build_cubic = [](const auto& j) {
+ auto N = j["pures"].size();
+ std::valarray<double> a0i(N), bi(N), c1(N), Tc(N);
+ std::size_t i = 0;
+ for (auto p : j["pures"]) {
+ a0i[i] = p["a0i / Pa m^6/mol^2"];
+ bi[i] = p["bi / m^3/mol"];
+ c1[i] = p["c1"];
+ Tc[i] = p["Tc / K"];
+ i++;
+ }
+ return CPACubic(get_cubic_flag(j["cubic"]), a0i, bi, c1, Tc);
+ };
+ auto build_assoc = [](const auto &&cubic, const auto& j) {
+ auto N = j["pures"].size();
+ std::vector<association_classes> classes;
+ std::valarray<double> epsABi(N), betaABi(N);
+ std::size_t i = 0;
+ for (auto p : j["pures"]) {
+ epsABi[i] = p["epsABi / J/mol"];
+ betaABi[i] = p["betaABi"];
+ classes.push_back(get_association_classes(p["class"]));
+ i++;
+ }
+ return CPAAssociation(std::move(cubic), classes, epsABi, betaABi);
+ };
+ return CPA(build_cubic(j), build_assoc(build_cubic(j), j));
+}
+
+}; /* namespace CPA */
\ No newline at end of file
diff --git a/src/tests/catch_tests.cxx b/src/tests/catch_tests.cxx
index efb6e0f..87603ee 100644
--- a/src/tests/catch_tests.cxx
+++ b/src/tests/catch_tests.cxx
@@ -311,7 +311,7 @@ TEST_CASE("Test pure VLE with non-unity R0/Rr", "") {
TEST_CASE("Test water", "") {
std::valarray<double> a0 = {0.12277}, bi = {0.000014515}, c1 = {0.67359}, Tc = {647.096},
molefrac = {1.0};
- CPACubic cub(cubic_flag::SRK, a0, bi, c1, Tc);
+ CPA::CPACubic cub(CPA::cubic_flag::SRK, a0, bi, c1, Tc);
double T = 400, rhomolar = 100;
auto R = 8.3144598;
@@ -322,15 +322,37 @@ TEST_CASE("Test water", "") {
double p_noassoc = T*rhomolar*R*(1+tdx::get_Ar01(cub, T, rhomolar, z));
CAPTURE(p_noassoc);
- std::vector<association_classes> schemes = { association_classes::a4C };
+ std::vector<CPA::association_classes> schemes = { CPA::association_classes::a4C };
std::valarray<double> epsAB = { 16655 }, betaAB = { 0.0692 };
- CPAAssociation cpaa(cub, schemes, epsAB, betaAB);
+ CPA::CPAAssociation cpaa(std::move(cub), schemes, epsAB, betaAB);
- CPA cpa(cub, cpaa);
+ CPA::CPA cpa(std::move(cub), std::move(cpaa));
using tdc = TDXDerivatives<decltype(cpa)>;
auto Ar01 = tdc::get_Ar01(cpa, T, rhomolar, z);
double p_withassoc = T*rhomolar*R*(1 + Ar01);
CAPTURE(p_withassoc);
REQUIRE(p_withassoc == 3.14);
+}
+
+TEST_CASE("Test water w/ factory", "") {
+ using namespace CPA;
+ nlohmann::json water = {
+ {"a0i / Pa m^6/mol^2",0.12277 }, {"bi / m^3/mol", 0.000014515}, {"c1", 0.67359}, {"Tc / K", 647.096},
+ {"epsABi / J/mol", 16655.0}, {"betaABi",0.0692}, {"class","4C"}
+ };
+ nlohmann::json j = {
+ {"cubic","SRK"},
+ {"pures", {water}}
+ };
+ auto cpa = CPAfactory(j);
+
+ double T = 400, rhomolar = 100, R = 8.3144598;
+ auto z = (Eigen::ArrayXd(1) << 1).finished();
+ using tdc = TDXDerivatives<decltype(cpa)>;
+ auto Ar01 = tdc::get_Ar01(cpa, T, rhomolar, z);
+ double p_withassoc = T * rhomolar * R * (1 + Ar01);
+ CAPTURE(p_withassoc);
+
+ REQUIRE(p_withassoc == 3.14);
}
\ No newline at end of file
--
GitLab