diff --git a/.github/workflows/runcatch.yml b/.github/workflows/runcatch.yml
index ef2e281ade3b9c876e68836096dc6b1f5505e98e..2f2a6281910378b838581bd2cf2c2cbbd5215c29 100644
--- a/.github/workflows/runcatch.yml
+++ b/.github/workflows/runcatch.yml
@@ -19,6 +19,8 @@ jobs:
run: cmake --build build --target catch_tests
- name: build multifluid
run: cmake --build build --target multifluid
+ - name: build teqp
+ run: cmake --build build --target teqp
# run tests
- name: run Catch tests
run: build/catch_tests
\ No newline at end of file
diff --git a/include/teqp/models/CPA.hpp b/include/teqp/models/CPA.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..a0419aaef58feee3a2d323ca0b78f473a69d6f65
--- /dev/null
+++ b/include/teqp/models/CPA.hpp
@@ -0,0 +1,279 @@
+#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
+template<typename BType, typename TType, typename RhoType, typename VecType>
+auto get_DeltaAB_pure(radial_dist dist, double epsABi, double betaABi, BType b_cubic, TType RT, RhoType rhomolar, const VecType& molefrac) {
+
+ using eta_type = std::common_type_t<decltype(rhomolar), decltype(b_cubic)>;
+ eta_type eta;
+ eta_type g_vm_ref;
+
+ // Calculate the contact value of the radial distribution function g(v)
+ switch (dist) {
+ case radial_dist::CS: {
+ // Carnahan - Starling EOS, given by Kontogeorgis et al., Ind.Eng.Chem.Res. 2006, 45, 4855 - 4868, Eq. 4a and 4b:
+ eta = (rhomolar / 4.0) * b_cubic;
+ g_vm_ref = (2.0 - eta) / (2.0 * POW3(1.0 - eta));
+ break;
+ }
+ case radial_dist::KG: {
+ // Function proposed by Kontogeorgis, G.M.; Yakoumis, I.V.; Meijer, H.; Hendriks, E.M.; Moorwood, T., Fluid Phase Equilib. 1999, 158 - 160, 201.
+ eta = (rhomolar / 4.0) * b_cubic;
+ g_vm_ref = 1.0 / (1.0 - 1.9 * eta);
+ break;
+ }
+ case radial_dist::OT: {
+ g_vm_ref = 1.0 / (1.0 - 0.475 * rhomolar * b_cubic);
+ break;
+ }
+ default: {
+ throw std::invalid_argument("Bad radial_dist");
+ }
+ }
+
+ // Calculate the association strength between site Ai and Bi for a pure compent
+ auto DeltaAiBj = forceeval(g_vm_ref*(exp(epsABi/RT) - 1.0)*b_cubic* betaABi);
+
+ return DeltaAiBj;
+};
+
+/// Routine that calculates the fractions of sites Ai not bound to other active sites for pure fluids
+/// Some association schemes are explicitly solvable for self-associating compounds, see Huang and Radosz, Ind. Eng. Chem. Res., 29 (11), 1990
+/// So far implemented association schemes : 1A, 2B, 3B, 4C (see Kontogeorgis et al., Ind. Eng. Chem. Res. 2006, 45, 4855 - 4868)
+///
+
+template<typename BType, typename TType, typename RhoType, typename VecType>
+auto XA_calc_pure(int N_sites, association_classes scheme, double epsABi, double betaABi, const BType b_cubic, const TType RT, const RhoType rhomolar, const VecType& molefrac) {
+
+ // Matrix XA(A, j) that contains all of the fractions of sites A not bonded to other active sites for each molecule i
+ // Start values for the iteration(set all sites to non - bonded, = 1)
+ using result_type = std::common_type_t<decltype(RT), decltype(rhomolar), decltype(molefrac[0])>;
+ Eigen::Array<result_type, Eigen::Dynamic, Eigen::Dynamic> XA; // A maximum of 4 association sites(A, B, C, D)
+ XA.resize(N_sites, 1);
+ XA.setOnes();
+
+ // Get the association strength between the associating sites
+ auto dist = radial_dist::KG; // TODO: pass this in
+ auto DeltaAiBj = get_DeltaAB_pure(dist, epsABi, betaABi, b_cubic, RT, rhomolar, molefrac);
+
+ if (scheme == association_classes::a1A) { // Acids
+ // Only one association site "A" (OH - group with C = O - group)
+ XA(0, 0) = forceeval((-1.0 + sqrt(1.0 + 4.0 * rhomolar * DeltaAiBj)) / (2.0 * rhomolar * DeltaAiBj));
+ }
+ else if (scheme == association_classes::a2B) { // Alcohols
+ // Two association sites "A" and "B"
+ XA(0, 0) = forceeval((-1.0 + sqrt(1.0 + 4.0 * rhomolar * DeltaAiBj)) / (2.0 * rhomolar * DeltaAiBj));
+ XA(1, 0) = XA(0, 0); // XB = XA;
+ }
+ else if (scheme == association_classes::a3B) { // Glycols
+ // Three association sites "A", "B", "C"
+ XA(0, 0) = forceeval((-(1.0 - rhomolar * DeltaAiBj) + sqrt(POW2(1.0 + rhomolar * DeltaAiBj) + 4.0 * rhomolar * DeltaAiBj)) / (4.0 * rhomolar * DeltaAiBj));
+ XA(1, 0) = XA(0, 0); // XB = XA
+ XA(2, 0) = 2.0*XA(0, 0) - 1.0; // XC = 2XA - 1
+ }
+ else if (scheme == association_classes::a4C) { // Water
+ // Four association sites "A", "B", "C", "D"
+ XA(0, 0) = forceeval((-1.0 + sqrt(1.0 + 8.0 * rhomolar * DeltaAiBj)) / (4.0 * rhomolar * DeltaAiBj));
+ XA(1, 0) = XA(0, 0); // XB = XA
+ XA(2, 0) = XA(0, 0); // XC = XA
+ XA(3, 0) = XA(0, 0); // XD = XA
+ }
+ else if (scheme == association_classes::not_associating) { // non - associating compound
+ XA(0, 0) = 1;
+ XA(1, 0) = 1;
+ XA(2, 0) = 1;
+ XA(3, 0) = 1;
+ }
+ else {
+ throw std::invalid_argument("Bad scheme");
+ }
+ return XA;
+};
+
+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:
+
+ std::valarray<double> a0, bi, c1, Tc;
+ double delta_1, delta_2;
+ 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) {
+ switch (flag) {
+ case cubic_flag::PR:
+ { delta_1 = 1 + sqrt(2); delta_2 = 1 - sqrt(2); break; }
+ case cubic_flag::SRK:
+ { delta_1 = 0; delta_2 = 1; break; }
+ default:
+ throw std::invalid_argument("Bad cubic flag");
+ }
+ k_ij.resize(Tc.size()); for (auto i = 0; i < k_ij.size(); ++i) { k_ij[i].resize(Tc.size()); }
+ };
+
+ template<typename TType>
+ auto get_ai(TType T, int i) const {
+ return a0[i] * POW2(1.0 + c1[i]*(1.0 - sqrt(T / Tc[i])));
+ }
+
+ template<typename TType, typename VecType>
+ auto get_ab(const TType T, const VecType& molefrac) const {
+ using return_type = std::common_type_t<decltype(T), decltype(molefrac[0])>;
+ return_type asummer = 0.0, bsummer = 0.0;
+ for (auto i = 0; i < molefrac.size(); ++i) {
+ bsummer += molefrac[i] * bi[i];
+ auto ai = get_ai(T, i);
+ for (auto j = 0; j < molefrac.size(); ++j) {
+ auto aj = get_ai(T, j);
+ auto a_ij = (1.0 - k_ij[i][j]) * sqrt(ai * aj);
+ asummer += molefrac[i] * molefrac[j] * a_ij;
+ }
+ }
+ return std::make_tuple(asummer, bsummer);
+ }
+
+ template<typename TType, typename RhoType, typename VecType, typename RType>
+ auto alphar(const TType T, const RhoType rhomolar, const VecType& molefrac, const RType& R_gas) const {
+ auto [a_cubic, b_cubic] = get_ab(T, molefrac);
+ return forceeval(-log(1.0 - b_cubic * rhomolar) - a_cubic / R_gas / T * log((delta_1*b_cubic*rhomolar + 1.0) / (delta_2*b_cubic*rhomolar + 1.0)) / b_cubic / (delta_1 - delta_2));
+ }
+};
+
+template<typename Cubic>
+class CPAAssociation {
+private:
+ const Cubic cubic;
+ const std::vector<association_classes> classes;
+ const std::vector<int> N_sites;
+ const std::valarray<double> epsABi, betaABi;
+
+ auto get_N_sites(const std::vector<association_classes> &classes) {
+ std::vector<int> N_sites_out;
+ auto get_N = [](auto cl) {
+ switch (cl) {
+ case association_classes::a1A: return 1;
+ case association_classes::a2B: return 2;
+ case association_classes::a3B: return 3;
+ case association_classes::a4C: return 4;
+ default: throw std::invalid_argument("Bad association class");
+ }
+ };
+ for (auto cl : classes) {
+ N_sites_out.push_back(get_N(cl));
+ }
+ return N_sites_out;
+ }
+
+public:
+ 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, typename RType>
+ auto alphar(const TType& T, const RhoType& rhomolar, const VecType& molefrac, const RType &R_gas) const {
+ // Calculate a and b of the mixture
+ auto [a_cubic, b_cubic] = cubic.get_ab(T, molefrac);
+
+ // Calculate the fraction of sites not bonded with other active sites
+ auto RT = forceeval(R_gas * T); // R times T
+ auto XA = XA_calc_pure(N_sites[0], classes[0], epsABi[0], betaABi[0], b_cubic, RT, rhomolar, molefrac);
+
+ using return_type = std::common_type_t<decltype(T), decltype(rhomolar), decltype(molefrac[0])>;
+ return_type alpha_r_asso = 0.0;
+ auto i = 0;
+ for (auto xi : molefrac){ // loop over all components
+ auto XAi = XA.col(i);
+ alpha_r_asso += forceeval(xi * (log(XAi) - XAi / 2).sum());
+ alpha_r_asso += xi*static_cast<double>(N_sites[i])/2;
+ i++;
+ }
+ return forceeval(alpha_r_asso);
+ }
+};
+
+template <typename Cubic, typename Assoc>
+class CPAEOS {
+public:
+ const Cubic cubic;
+ const Assoc assoc;
+
+ const double R = 1.380649e-23 * 6.02214076e23; ///< Exact value, given by k_B*N_A
+
+ CPAEOS(Cubic &&cubic, Assoc &&assoc) : cubic(cubic), assoc(assoc) {
+ }
+
+ /// Residual dimensionless Helmholtz energy from the SRK or PR core and contribution due to association
+ /// alphar = a/(R*T) where a and R are both molar quantities
+ template<typename TType, typename RhoType, typename VecType>
+ auto alphar(const TType& T, const RhoType& rhomolar, const VecType& molefrac) const {
+
+ // Calculate the contribution to alphar from the conventional cubic EOS
+ auto alpha_r_cubic = cubic.alphar(T, rhomolar, molefrac, R);
+
+ // Calculate the contribution to alphar from association
+ auto alpha_r_assoc = assoc.alphar(T, rhomolar, molefrac, R);
+
+ return forceeval(alpha_r_cubic + alpha_r_assoc);
+ }
+};
+
+/// 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 CPAEOS(build_cubic(j), build_assoc(build_cubic(j), j));
+}
+
+}; /* namespace CPA */
\ No newline at end of file
diff --git a/interface/pybind11_wrapper.cpp b/interface/pybind11_wrapper.cpp
index 7b13eb231910d258541395e3a62c85095c43dc9b..0d2186255892b1af1fdb5d9d7dcc982170fc8d93 100644
--- a/interface/pybind11_wrapper.cpp
+++ b/interface/pybind11_wrapper.cpp
@@ -1,5 +1,6 @@
#define USE_AUTODIFF
+#include "nlohmann/json.hpp"
#include "pybind11_json/pybind11_json.hpp"
#include <pybind11/pybind11.h>
@@ -12,6 +13,7 @@
#include "teqp/algorithms/critical_tracing.hpp"
#include "teqp/models/pcsaft.hpp"
+#include "teqp/models/CPA.hpp"
#include "teqp/models/multifluid.hpp"
#include "teqp/algorithms/VLE.hpp"
@@ -96,13 +98,18 @@ void init_teqp(py::module& m) {
// Multifluid model
m.def("build_multifluid_model", &build_multifluid_model);
using MultiFluid = decltype(build_multifluid_model(std::vector<std::string>{"",""},"",""));
- using idMF = IsochoricDerivatives<MultiFluid, double, Eigen::Array<double, Eigen::Dynamic, 1> >;
auto wMF = py::class_<MultiFluid>(m, "MultiFluid")
.def("get_Tcvec", [](const MultiFluid& c) { return c.redfunc.Tc; })
.def("get_vcvec", [](const MultiFluid& c) { return c.redfunc.vc; })
;
add_derivatives<MultiFluid>(m, wMF);
+ // CPA model
+ using CPAEOS_ = decltype(CPA::CPAfactory(nlohmann::json()));
+ m.def("CPAfactory", &CPA::CPAfactory);
+ auto wCPA = py::class_<CPAEOS_>(m, "CPAEOS");
+ add_derivatives<CPAEOS_>(m, wCPA);
+
// Some functions for timing overhead of interface
m.def("___mysummer", [](const double &c, const Eigen::ArrayXd &x) { return c*x.sum(); });
using RAX = Eigen::Ref<Eigen::ArrayXd>;
diff --git a/src/tests/catch_tests.cxx b/src/tests/catch_tests.cxx
index af9ea154b2001821e9ec5694d602e1519ccb2235..69045270746fcb918cae5cfedf56c74635f4eff6 100644
--- a/src/tests/catch_tests.cxx
+++ b/src/tests/catch_tests.cxx
@@ -4,6 +4,7 @@
#include "teqp/core.hpp"
#include "teqp/models/pcsaft.hpp"
#include "teqp/models/cubicsuperancillary.hpp"
+#include "teqp/models/CPA.hpp"
#include "teqp/algorithms/VLE.hpp"
auto build_vdW_argon() {
@@ -305,4 +306,53 @@ TEST_CASE("Test pure VLE with non-unity R0/Rr", "") {
auto r1 = residspecial.call(solnspecial);
auto rr = 0;
+}
+
+TEST_CASE("Test water", "") {
+ std::valarray<double> a0 = {0.12277}, bi = {0.000014515}, c1 = {0.67359}, Tc = {647.096},
+ molefrac = {1.0};
+ CPA::CPACubic cub(CPA::cubic_flag::SRK, a0, bi, c1, Tc);
+ double T = 400, rhomolar = 100;
+
+ auto R = 8.3144598;
+ auto z = (Eigen::ArrayXd(1) << 1).finished();
+
+ using tdx = TDXDerivatives<decltype(cub)>;
+ auto alphar = cub.alphar(T, rhomolar, molefrac);
+ double p_noassoc = T*rhomolar*R*(1+tdx::get_Ar01(cub, T, rhomolar, z));
+ CAPTURE(p_noassoc);
+
+ std::vector<CPA::association_classes> schemes = { CPA::association_classes::a4C };
+ std::valarray<double> epsAB = { 16655 }, betaAB = { 0.0692 };
+ CPA::CPAAssociation cpaa(std::move(cub), schemes, epsAB, betaAB);
+
+ CPA::CPAEOS 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