From d58c86281579601baf471570e475f0dd8a5600d7 Mon Sep 17 00:00:00 2001 From: Ian Bell Date: Mon, 26 Feb 2024 19:23:42 -0500 Subject: [PATCH 1/9] Implement the LKP EOS --- include/teqp/models/LKP.hpp | 74 +++++++++++++++++++++++++++++++++++++ 1 file changed, 74 insertions(+) create mode 100644 include/teqp/models/LKP.hpp diff --git a/include/teqp/models/LKP.hpp b/include/teqp/models/LKP.hpp new file mode 100644 index 00000000..bf9f71cc --- /dev/null +++ b/include/teqp/models/LKP.hpp @@ -0,0 +1,74 @@ +#pragma once // Only include header once in a compilation unit if it is included multiple times + +#include "teqp/constants.hpp" // used for R +#include "teqp/types.hpp" // needed for forceeval + +namespace teqp{ +namespace LKP{ + +struct LKPFluidParameters { + std::vector b, c, d; + double beta, gamma_, omega; +}; + +class LKPMix{ +public: + const LKPFluidParameters simple{{0.0, 0.1181193, 0.265728, 0.154790, 0.303230e-1}, // b, with pad to match indices + {0.0, 0.236744e-1, 0.186984e-1, 0, 0.4217240e-1}, // c, with pad to match indices + {0, 0.155428e-4, 0.623689e-4}, // d, with pad to match indices + 0.653920, 0.601670e-1, 0.0}, // beta, gamma, omega + ref{{0, 0.2026579, 0.331511, 0.276550e-1, 0.203488}, // b, with pad to match indices + {0, 0.313385e-1, 0.503618e-1, 0.169010e-1,0.41577e-1}, // c, with pad to match indices + {0, 0.487360e-4, 0.740336e-5}, // d, with pad to match indices + 1.226, 0.03754, 0.3978}; // beta, gamma, omega + const std::vector Tcrit, pcrit, acentric; + const std::vector> kmat; + + LKPMix(const std::vector& Tcrit, const std::vector& pcrit, const std::vector& acentric, const std::vector>& kmat) : Tcrit(Tcrit), pcrit(pcrit), acentric(acentric), kmat(kmat){} + + template + auto R(const VecType& molefrac) const { return teqp::constants::R_CODATA2017; } + + /// Calculate the contribution for one of the fluids, depending on the parameter set passed in + template + auto alphar_func(const TTYPE& tau, const RhoType& delta, const ZcType& Zc, const LKPFluidParameters& params) const { + auto theta = 1/tau; + auto B = params.b[1] - params.b[2]*tau - params.b[3]*powi(tau, 2) - params.b[4]*powi(tau, 3); + auto C = params.c[1] - params.c[2]*tau + params.c[3]*powi(tau, 2); + auto D = params.d[1] + params.d[2]*tau; + return forceeval(B/Zc*delta + 1.0/2.0*C/powi(delta/Zc, 2) + 1.0/5.0*D*powi(delta/Zc, 5) - params.c[4]*powi(tau, 3)/(2*params.gamma_)*(params.gamma_*powi(delta/Zc, 2)+params.beta+1.0)*exp(-params.gamma_*powi(delta/Zc, 2)) + params.c[4]*powi(tau,3)/(2*params.gamma_)*(params.beta+1.0)); + } + + template + auto alphar(const TTYPE& T, const RhoType& rhomolar, const VecType& mole_fractions) const { + + const VecType& x = mole_fractions; // just an alias to save typing, no copy is invoked + std::decay_t summer_omega = 0.0, summer_vcmix = 0.0, summer_Tcmix = 0.0; + double Ru = teqp::constants::R_CODATA2017; + + for (auto i = 0; i < mole_fractions.size(); ++i){ + summer_omega += mole_fractions[i]*acentric[i]; + auto v_ci = (0.2905-0.085*acentric[i])*Ru*Tcrit[i]/pcrit[i]; + for (auto j = 0; j < mole_fractions.size(); ++j){ + auto v_cj = (0.2905-0.085*acentric[j])*Ru*Tcrit[j]/pcrit[j]; + auto v_c_ij = 1.0/8.0*(cbrt(v_ci) + cbrt(v_cj)); + auto T_c_ij = kmat[i][j]*sqrt(Tcrit[i]*Tcrit[j]); + summer_vcmix += x[i]*x[j]*v_c_ij; + summer_Tcmix += x[i]*x[j]*pow(v_c_ij, 0.25)*T_c_ij; + } + } + auto omega_mix = summer_omega; + auto vc_mix = summer_vcmix; + auto Tc_mix = 1.0/pow(summer_vcmix, 0.25)*summer_Tcmix; + auto pc_mix = (0.2905-0.085*omega_mix)*Ru*Tc_mix/vc_mix; + auto Zc = pc_mix*vc_mix/(Ru*Tc_mix); + auto tau = Tc_mix/T; + auto delta = vc_mix*rhomolar; + + auto retval = (1.0-omega_mix/ref.omega)*alphar_func(tau, delta, Zc, simple) + (omega_mix/ref.omega)*alphar_func(tau, delta, Zc, ref); + return forceeval(retval); + } +}; + +} +} From 23b4251df9ce010b76635e83a17e25ce16cf8660 Mon Sep 17 00:00:00 2001 From: Ian Bell Date: Mon, 26 Feb 2024 21:40:54 -0500 Subject: [PATCH 2/9] Add LKP test --- src/tests/catch_tests.cxx | 14 ++++++++++++++ 1 file changed, 14 insertions(+) diff --git a/src/tests/catch_tests.cxx b/src/tests/catch_tests.cxx index d6151356..17b65dc9 100644 --- a/src/tests/catch_tests.cxx +++ b/src/tests/catch_tests.cxx @@ -14,6 +14,7 @@ using Catch::Approx; #include using namespace boost::multiprecision; #include "teqp/finite_derivs.hpp" +#include "teqp/models/LKP.hpp" using namespace teqp; @@ -97,6 +98,19 @@ TEST_CASE("Check virial coefficients for vdW", "[virial]") } } +TEST_CASE("Check LKP", "[LKP]"){ + using namespace teqp::LKP; + // THF + water + std::vector Tc_K = { 540.2, 647.096 }; + std::vector pc_Pa = { 5304.44e3, 22064.0e3 }; + std::vector acentric = { 0.234, 0.3443 }; + std::vector> kmat = {{1,1}, {1,1}}; + auto model = LKPMix(Tc_K, pc_Pa, acentric, kmat); + auto T = 303.15, rhomolar = 15000.0; + auto z = (Eigen::ArrayXd(2) << 0.3, 0.7).finished(); + auto Ar01 = TDXDerivatives::get_Ar01(model, T, rhomolar, z); +} + TEST_CASE("Check Hessian of Psir", "[virial]") { From 50fc206fd725d0a3c777e525d262ab93cf0556a5 Mon Sep 17 00:00:00 2001 From: Ian Bell Date: Tue, 27 Feb 2024 14:53:06 -0500 Subject: [PATCH 3/9] Move LKP test into its own file --- src/tests/catch_test_LKP.cxx | 14 ++++++++++++++ src/tests/catch_tests.cxx | 14 -------------- 2 files changed, 14 insertions(+), 14 deletions(-) create mode 100644 src/tests/catch_test_LKP.cxx diff --git a/src/tests/catch_test_LKP.cxx b/src/tests/catch_test_LKP.cxx new file mode 100644 index 00000000..b5155c4e --- /dev/null +++ b/src/tests/catch_test_LKP.cxx @@ -0,0 +1,14 @@ +#include "teqp/models/LKP.hpp" + +TEST_CASE("Check LKP", "[LKP]"){ + using namespace teqp::LKP; + // THF + water + std::vector Tc_K = { 540.2, 647.096 }; + std::vector pc_Pa = { 5304.44e3, 22064.0e3 }; + std::vector acentric = { 0.234, 0.3443 }; + std::vector> kmat = {{1,1}, {1,1}}; + auto model = LKPMix(Tc_K, pc_Pa, acentric, kmat); + auto T = 303.15, rhomolar = 15000.0; + auto z = (Eigen::ArrayXd(2) << 0.3, 0.7).finished(); + auto Ar01 = TDXDerivatives::get_Ar01(model, T, rhomolar, z); +} diff --git a/src/tests/catch_tests.cxx b/src/tests/catch_tests.cxx index 17b65dc9..d6151356 100644 --- a/src/tests/catch_tests.cxx +++ b/src/tests/catch_tests.cxx @@ -14,7 +14,6 @@ using Catch::Approx; #include using namespace boost::multiprecision; #include "teqp/finite_derivs.hpp" -#include "teqp/models/LKP.hpp" using namespace teqp; @@ -98,19 +97,6 @@ TEST_CASE("Check virial coefficients for vdW", "[virial]") } } -TEST_CASE("Check LKP", "[LKP]"){ - using namespace teqp::LKP; - // THF + water - std::vector Tc_K = { 540.2, 647.096 }; - std::vector pc_Pa = { 5304.44e3, 22064.0e3 }; - std::vector acentric = { 0.234, 0.3443 }; - std::vector> kmat = {{1,1}, {1,1}}; - auto model = LKPMix(Tc_K, pc_Pa, acentric, kmat); - auto T = 303.15, rhomolar = 15000.0; - auto z = (Eigen::ArrayXd(2) << 0.3, 0.7).finished(); - auto Ar01 = TDXDerivatives::get_Ar01(model, T, rhomolar, z); -} - TEST_CASE("Check Hessian of Psir", "[virial]") { From d4a42a82e9fe183438d54daab49c924b2754406b Mon Sep 17 00:00:00 2001 From: Ian Bell Date: Tue, 27 Feb 2024 15:51:07 -0500 Subject: [PATCH 4/9] Fix three single-character typo in LKP --- include/teqp/models/LKP.hpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/include/teqp/models/LKP.hpp b/include/teqp/models/LKP.hpp index bf9f71cc..4cda8d2d 100644 --- a/include/teqp/models/LKP.hpp +++ b/include/teqp/models/LKP.hpp @@ -14,7 +14,7 @@ struct LKPFluidParameters { class LKPMix{ public: const LKPFluidParameters simple{{0.0, 0.1181193, 0.265728, 0.154790, 0.303230e-1}, // b, with pad to match indices - {0.0, 0.236744e-1, 0.186984e-1, 0, 0.4217240e-1}, // c, with pad to match indices + {0.0, 0.236744e-1, 0.186984e-1, 0, 0.427240e-1}, // c, with pad to match indices {0, 0.155428e-4, 0.623689e-4}, // d, with pad to match indices 0.653920, 0.601670e-1, 0.0}, // beta, gamma, omega ref{{0, 0.2026579, 0.331511, 0.276550e-1, 0.203488}, // b, with pad to match indices @@ -36,7 +36,7 @@ class LKPMix{ auto B = params.b[1] - params.b[2]*tau - params.b[3]*powi(tau, 2) - params.b[4]*powi(tau, 3); auto C = params.c[1] - params.c[2]*tau + params.c[3]*powi(tau, 2); auto D = params.d[1] + params.d[2]*tau; - return forceeval(B/Zc*delta + 1.0/2.0*C/powi(delta/Zc, 2) + 1.0/5.0*D*powi(delta/Zc, 5) - params.c[4]*powi(tau, 3)/(2*params.gamma_)*(params.gamma_*powi(delta/Zc, 2)+params.beta+1.0)*exp(-params.gamma_*powi(delta/Zc, 2)) + params.c[4]*powi(tau,3)/(2*params.gamma_)*(params.beta+1.0)); + return forceeval(B/Zc*delta + 1.0/2.0*C*powi(delta/Zc, 2) + 1.0/5.0*D*powi(delta/Zc, 5) - params.c[4]*powi(tau, 3)/(2*params.gamma_)*(params.gamma_*powi(delta/Zc, 2)+params.beta+1.0)*exp(-params.gamma_*powi(delta/Zc, 2)) + params.c[4]*powi(tau,3)/(2*params.gamma_)*(params.beta+1.0)); } template @@ -51,7 +51,7 @@ class LKPMix{ auto v_ci = (0.2905-0.085*acentric[i])*Ru*Tcrit[i]/pcrit[i]; for (auto j = 0; j < mole_fractions.size(); ++j){ auto v_cj = (0.2905-0.085*acentric[j])*Ru*Tcrit[j]/pcrit[j]; - auto v_c_ij = 1.0/8.0*(cbrt(v_ci) + cbrt(v_cj)); + auto v_c_ij = 1.0/8.0*powi(cbrt(v_ci) + cbrt(v_cj), 3); auto T_c_ij = kmat[i][j]*sqrt(Tcrit[i]*Tcrit[j]); summer_vcmix += x[i]*x[j]*v_c_ij; summer_Tcmix += x[i]*x[j]*pow(v_c_ij, 0.25)*T_c_ij; From 48192fe610d7d36fb1bd4d74018b88aaa13f578c Mon Sep 17 00:00:00 2001 From: Ian Bell Date: Tue, 27 Feb 2024 15:52:18 -0500 Subject: [PATCH 5/9] Check sizes of vectors of parameters --- include/teqp/models/LKP.hpp | 19 ++++++++++++++++--- 1 file changed, 16 insertions(+), 3 deletions(-) diff --git a/include/teqp/models/LKP.hpp b/include/teqp/models/LKP.hpp index 4cda8d2d..2f4466d7 100644 --- a/include/teqp/models/LKP.hpp +++ b/include/teqp/models/LKP.hpp @@ -1,7 +1,10 @@ #pragma once // Only include header once in a compilation unit if it is included multiple times +#include + #include "teqp/constants.hpp" // used for R #include "teqp/types.hpp" // needed for forceeval +#include "teqp/exceptions.hpp" // to return teqp error messages namespace teqp{ namespace LKP{ @@ -22,12 +25,18 @@ class LKPMix{ {0, 0.487360e-4, 0.740336e-5}, // d, with pad to match indices 1.226, 0.03754, 0.3978}; // beta, gamma, omega const std::vector Tcrit, pcrit, acentric; + const double m_R; ///< molar gas constant to be used in this model, in J/mol/K const std::vector> kmat; - LKPMix(const std::vector& Tcrit, const std::vector& pcrit, const std::vector& acentric, const std::vector>& kmat) : Tcrit(Tcrit), pcrit(pcrit), acentric(acentric), kmat(kmat){} + LKPMix(const std::vector& Tcrit, const std::vector& pcrit, const std::vector& acentric, double R, const std::vector>& kmat) : Tcrit(Tcrit), pcrit(pcrit), acentric(acentric), m_R(R), kmat(kmat){ + if (std::set{Tcrit.size(), pcrit.size(), acentric.size()}.size() > 1){ + throw teqp::InvalidArgument("The arrays should all be the same size."); + } + + } template - auto R(const VecType& molefrac) const { return teqp::constants::R_CODATA2017; } + auto R(const VecType& molefrac) const { return m_R; } /// Calculate the contribution for one of the fluids, depending on the parameter set passed in template @@ -42,9 +51,13 @@ class LKPMix{ template auto alphar(const TTYPE& T, const RhoType& rhomolar, const VecType& mole_fractions) const { + if (mole_fractions.size() != acentric.size()){ + throw teqp::InvalidArgument("The mole fractions should be of of size "+ std::to_string(acentric.size())); + } + const VecType& x = mole_fractions; // just an alias to save typing, no copy is invoked std::decay_t summer_omega = 0.0, summer_vcmix = 0.0, summer_Tcmix = 0.0; - double Ru = teqp::constants::R_CODATA2017; + double Ru = m_R; for (auto i = 0; i < mole_fractions.size(); ++i){ summer_omega += mole_fractions[i]*acentric[i]; From d8c637beca1d0715873f8015782faf4d77bdf470 Mon Sep 17 00:00:00 2001 From: Ian Bell Date: Tue, 27 Feb 2024 15:52:25 -0500 Subject: [PATCH 6/9] One more typo --- include/teqp/models/LKP.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/include/teqp/models/LKP.hpp b/include/teqp/models/LKP.hpp index 2f4466d7..f607ed60 100644 --- a/include/teqp/models/LKP.hpp +++ b/include/teqp/models/LKP.hpp @@ -33,7 +33,7 @@ class LKPMix{ throw teqp::InvalidArgument("The arrays should all be the same size."); } - } + } template auto R(const VecType& molefrac) const { return m_R; } @@ -43,7 +43,7 @@ class LKPMix{ auto alphar_func(const TTYPE& tau, const RhoType& delta, const ZcType& Zc, const LKPFluidParameters& params) const { auto theta = 1/tau; auto B = params.b[1] - params.b[2]*tau - params.b[3]*powi(tau, 2) - params.b[4]*powi(tau, 3); - auto C = params.c[1] - params.c[2]*tau + params.c[3]*powi(tau, 2); + auto C = params.c[1] - params.c[2]*tau + params.c[3]*powi(tau, 3); auto D = params.d[1] + params.d[2]*tau; return forceeval(B/Zc*delta + 1.0/2.0*C*powi(delta/Zc, 2) + 1.0/5.0*D*powi(delta/Zc, 5) - params.c[4]*powi(tau, 3)/(2*params.gamma_)*(params.gamma_*powi(delta/Zc, 2)+params.beta+1.0)*exp(-params.gamma_*powi(delta/Zc, 2)) + params.c[4]*powi(tau,3)/(2*params.gamma_)*(params.beta+1.0)); } From 00144b80f992ae4a08184a208c4decc2bf952555 Mon Sep 17 00:00:00 2001 From: Ian Bell Date: Tue, 27 Feb 2024 15:53:02 -0500 Subject: [PATCH 7/9] Implement check values from TREND 5.0 --- src/tests/catch_test_LKP.cxx | 114 ++++++++++++++++++++++++++++++++--- 1 file changed, 105 insertions(+), 9 deletions(-) diff --git a/src/tests/catch_test_LKP.cxx b/src/tests/catch_test_LKP.cxx index b5155c4e..0ec81c42 100644 --- a/src/tests/catch_test_LKP.cxx +++ b/src/tests/catch_test_LKP.cxx @@ -1,14 +1,110 @@ +#include +#include +using Catch::Approx; + +#include "teqp/derivs.hpp" +#include "teqp/constants.hpp" #include "teqp/models/LKP.hpp" TEST_CASE("Check LKP", "[LKP]"){ using namespace teqp::LKP; - // THF + water - std::vector Tc_K = { 540.2, 647.096 }; - std::vector pc_Pa = { 5304.44e3, 22064.0e3 }; - std::vector acentric = { 0.234, 0.3443 }; - std::vector> kmat = {{1,1}, {1,1}}; - auto model = LKPMix(Tc_K, pc_Pa, acentric, kmat); - auto T = 303.15, rhomolar = 15000.0; - auto z = (Eigen::ArrayXd(2) << 0.3, 0.7).finished(); - auto Ar01 = TDXDerivatives::get_Ar01(model, T, rhomolar, z); + + SECTION("pure"){ + // methane, check values from TREND + std::vector Tc_K = {190.564}; + std::vector pc_Pa = {4.5992e6}; + std::vector acentric = {0.011}; + std::vector> kmat{{1.0}}; + double R = 8.3144598; + auto model = LKPMix(Tc_K, pc_Pa, acentric, R, kmat); + auto z = (Eigen::ArrayXd(1) << 1.0).finished(); + + teqp::TDXDerivatives::get_Ar00(model, 300.0, 8000.0, z); + + struct Point{ double T, rhomolar, alphar_expected; }; + for (auto& pt : std::vector{ + {2.5555555555556E+02,1.7778555555556E+03,-1.0483673775469E-01}, + {2.5555555555556E+02,2.6667333333333E+03,-1.5420416099894E-01}, + {2.5555555555556E+02,3.5556111111111E+03,-2.0148143793072E-01}, + {2.5555555555556E+02,4.4444888888889E+03,-2.4660690326310E-01}, + {2.5555555555556E+02,5.3333666666667E+03,-2.8951201335592E-01}, + {2.5555555555556E+02,6.2222444444444E+03,-3.3012829826090E-01}, + {2.5555555555556E+02,7.1111222222222E+03,-3.6839324530079E-01}, + {2.5555555555556E+02,8.0000000000000E+03,-4.0425381376141E-01}, + {2.7777777777778E+02,1.0000000000000E-01,-5.0166505188877E-06}, + {2.7777777777778E+02,8.8897777777778E+02,-4.3666195153684E-02}, + {2.7777777777778E+02,1.7778555555556E+03,-8.5458261689186E-02}, + {2.7777777777778E+02,2.6667333333333E+03,-1.2535703634585E-01}, + {2.7777777777778E+02,3.5556111111111E+03,-1.6332285849984E-01}, + {2.7777777777778E+02,4.4444888888889E+03,-1.9930364880452E-01}, + {2.7777777777778E+02,5.3333666666667E+03,-2.3323937266306E-01}, + {2.7777777777778E+02,6.2222444444444E+03,-2.6506680393110E-01}, + {2.7777777777778E+02,7.1111222222222E+03,-2.9472345661764E-01}, + {2.7777777777778E+02,8.0000000000000E+03,-3.2214967196215E-01}, + {3.0000000000000E+02,1.0000000000000E-01,-4.1178433147697E-06}, + {3.0000000000000E+02,8.8897777777778E+02,-3.5720133653505E-02}, + {3.0000000000000E+02,1.7778555555556E+03,-6.9656790541858E-02}, + {3.0000000000000E+02,2.6667333333333E+03,-1.0179417406048E-01}, + {3.0000000000000E+02,3.5556111111111E+03,-1.3209892931558E-01}, + {3.0000000000000E+02,4.4444888888889E+03,-1.6052611602765E-01}, + {3.0000000000000E+02,5.3333666666667E+03,-1.8702217613376E-01}, + {3.0000000000000E+02,6.2222444444444E+03,-2.1152813972213E-01}, + {3.0000000000000E+02,7.1111222222222E+03,-2.3398217049700E-01}, + {3.0000000000000E+02,8.0000000000000E+03,-2.5432064692151E-01} + }){ + auto alphar_actual = teqp::TDXDerivatives::get_Ar00(model, pt.T, pt.rhomolar, z); + CHECK(alphar_actual == Approx(pt.alphar_expected).margin(1e-12)); + } + } + + SECTION("methane + nitrogen mix"){ + // methane + nitrogen, check values from TREND + std::vector Tc_K = {190.564, 126.192}; + std::vector pc_Pa = {4.5992e6, 3.3958e6}; + std::vector acentric = {0.011, 0.037}; + std::vector> kmat{{1.0, 0.977},{0.977, 1.0}}; + double R = 8.3144598; + auto model = LKPMix(Tc_K, pc_Pa, acentric, R, kmat); + auto z = (Eigen::ArrayXd(2) << 0.8, 0.2).finished(); + + auto zbad = (Eigen::ArrayXd(3) << 0.3, 0.3, 0.4).finished(); + CHECK_THROWS(model.alphar(300.0, 8000.0, zbad)); + + struct Point{ double T, rhomolar, alphar_expected; }; + for (auto& pt : std::vector{ + {2.55555555555555E+02, 1.00000000000000E-01, -4.91536626760729E-06}, + {2.55555555555555E+02, 8.88977777777777E+02, -4.28036717751216E-02}, + {2.55555555555555E+02, 1.77785555555555E+03, -8.38097840422122E-02}, + {2.55555555555555E+02, 2.66673333333333E+03, -1.23001289158271E-01}, + {2.55555555555555E+02, 3.55561111111111E+03, -1.60342053977241E-01}, + {2.55555555555555E+02, 4.44448888888888E+03, -1.95784417210887E-01}, + {2.55555555555555E+02, 5.33336666666666E+03, -2.29273119192599E-01}, + {2.55555555555555E+02, 6.22224444444444E+03, -2.60749579654642E-01}, + {2.55555555555555E+02, 7.11112222222222E+03, -2.90155546674414E-01}, + {2.55555555555555E+02, 8.00000000000000E+03, -3.17435227813972E-01}, + {2.77777777777777E+02, 1.00000000000000E-01, -3.96079440378770E-06}, + {2.77777777777777E+02, 8.88977777777777E+02, -3.43639521841312E-02}, + {2.77777777777777E+02, 1.77785555555555E+03, -6.70251681046093E-02}, + {2.77777777777777E+02, 2.66673333333333E+03, -9.79698096970028E-02}, + {2.77777777777777E+02, 3.55561111111111E+03, -1.27167929631439E-01}, + {2.77777777777777E+02, 4.44448888888888E+03, -1.54578932250068E-01}, + {2.77777777777777E+02, 5.33336666666666E+03, -1.80154076480535E-01}, + {2.77777777777777E+02, 6.22224444444444E+03, -2.03839249743354E-01}, + {2.77777777777777E+02, 7.11112222222222E+03, -2.25577252196130E-01}, + {2.77777777777777E+02, 8.00000000000000E+03, -2.45308899079892E-01}, + {3.00000000000000E+02, 1.00000000000000E-01, -3.18316290879055E-06}, + {3.00000000000000E+02, 8.88977777777777E+02, -2.74782561342433E-02}, + {3.00000000000000E+02, 1.77785555555555E+03, -5.33101587343266E-02}, + {3.00000000000000E+02, 2.66673333333333E+03, -7.74841431899544E-02}, + {3.00000000000000E+02, 3.55561111111111E+03, -9.99748284135246E-02}, + {3.00000000000000E+02, 4.44448888888888E+03, -1.20746896159584E-01}, + {3.00000000000000E+02, 5.33336666666666E+03, -1.39756583212869E-01}, + {3.00000000000000E+02, 6.22224444444444E+03, -1.56953388700025E-01}, + {3.00000000000000E+02, 7.11112222222222E+03, -1.72281392690949E-01}, + {3.00000000000000E+02, 8.00000000000000E+03, -1.85679636571464E-01} + }){ + auto alphar_actual = teqp::TDXDerivatives::get_Ar00(model, pt.T, pt.rhomolar, z); + CHECK(alphar_actual == Approx(pt.alphar_expected).margin(1e-12)); + } + } } From 0cf634c8c0ab272aa5ddd957ae3d696e4fb4c2c5 Mon Sep 17 00:00:00 2001 From: Ian Bell Date: Tue, 27 Feb 2024 16:04:25 -0500 Subject: [PATCH 8/9] Add factory function and test it --- include/teqp/models/LKP.hpp | 8 +++++++- src/tests/catch_test_LKP.cxx | 15 +++++++++++++++ 2 files changed, 22 insertions(+), 1 deletion(-) diff --git a/include/teqp/models/LKP.hpp b/include/teqp/models/LKP.hpp index f607ed60..f61ae72e 100644 --- a/include/teqp/models/LKP.hpp +++ b/include/teqp/models/LKP.hpp @@ -6,6 +6,8 @@ #include "teqp/types.hpp" // needed for forceeval #include "teqp/exceptions.hpp" // to return teqp error messages +#include "nlohmann/json.hpp" + namespace teqp{ namespace LKP{ @@ -32,7 +34,6 @@ class LKPMix{ if (std::set{Tcrit.size(), pcrit.size(), acentric.size()}.size() > 1){ throw teqp::InvalidArgument("The arrays should all be the same size."); } - } template @@ -83,5 +84,10 @@ class LKPMix{ } }; +// Factory function takes in JSON data and returns an instance of the LKPMix class +auto make_LKPMix(const nlohmann::json& j){ + return LKPMix(j.at("Tcrit / K"), j.at("pcrit / Pa"), j.at("acentric"), j.at("R / J/mol/K"), j.at("kmat")); +} + } } diff --git a/src/tests/catch_test_LKP.cxx b/src/tests/catch_test_LKP.cxx index 0ec81c42..5042e51c 100644 --- a/src/tests/catch_test_LKP.cxx +++ b/src/tests/catch_test_LKP.cxx @@ -55,6 +55,21 @@ TEST_CASE("Check LKP", "[LKP]"){ auto alphar_actual = teqp::TDXDerivatives::get_Ar00(model, pt.T, pt.rhomolar, z); CHECK(alphar_actual == Approx(pt.alphar_expected).margin(1e-12)); } + + nlohmann::json spec{ + {"Tcrit / K", Tc_K}, + {"pcrit / Pa", pc_Pa}, + {"acentric", acentric}, + {"R / J/mol/K", R}, + {"kmat", kmat} + }; +// std::cout << spec.dump(2) << std::endl; + + CHECK_NOTHROW(make_LKPMix(spec)); + nlohmann::json badspec = spec; + badspec["kmat"] = 4.7; + CHECK_THROWS(make_LKPMix(badspec)); + } SECTION("methane + nitrogen mix"){ From 838054c6759c670007baa48a25cc223d35b9293f Mon Sep 17 00:00:00 2001 From: Ian Bell Date: Tue, 27 Feb 2024 16:28:43 -0500 Subject: [PATCH 9/9] Add and test C++ interface --- include/teqp/models/LKP.hpp | 18 ++++++------- include/teqp/models/fwd.hpp | 1 + interface/CPP/teqp_impl_factory.cpp | 2 ++ src/tests/catch_test_LKP.cxx | 42 ++++++++++++++++++----------- 4 files changed, 38 insertions(+), 25 deletions(-) diff --git a/include/teqp/models/LKP.hpp b/include/teqp/models/LKP.hpp index f61ae72e..4acef93f 100644 --- a/include/teqp/models/LKP.hpp +++ b/include/teqp/models/LKP.hpp @@ -37,22 +37,22 @@ class LKPMix{ } template - auto R(const VecType& molefrac) const { return m_R; } + auto R(const VecType& /*molefrac*/) const { return m_R; } /// Calculate the contribution for one of the fluids, depending on the parameter set passed in template auto alphar_func(const TTYPE& tau, const RhoType& delta, const ZcType& Zc, const LKPFluidParameters& params) const { - auto theta = 1/tau; auto B = params.b[1] - params.b[2]*tau - params.b[3]*powi(tau, 2) - params.b[4]*powi(tau, 3); auto C = params.c[1] - params.c[2]*tau + params.c[3]*powi(tau, 3); auto D = params.d[1] + params.d[2]*tau; - return forceeval(B/Zc*delta + 1.0/2.0*C*powi(delta/Zc, 2) + 1.0/5.0*D*powi(delta/Zc, 5) - params.c[4]*powi(tau, 3)/(2*params.gamma_)*(params.gamma_*powi(delta/Zc, 2)+params.beta+1.0)*exp(-params.gamma_*powi(delta/Zc, 2)) + params.c[4]*powi(tau,3)/(2*params.gamma_)*(params.beta+1.0)); + auto deltaZc = forceeval(delta/Zc); + return forceeval(B/Zc*delta + 1.0/2.0*C*powi(deltaZc, 2) + 1.0/5.0*D*powi(deltaZc, 5) - params.c[4]*powi(tau, 3)/(2*params.gamma_)*(params.gamma_*powi(deltaZc, 2)+params.beta+1.0)*exp(-params.gamma_*powi(deltaZc, 2)) + params.c[4]*powi(tau,3)/(2*params.gamma_)*(params.beta+1.0)); } template auto alphar(const TTYPE& T, const RhoType& rhomolar, const VecType& mole_fractions) const { - if (mole_fractions.size() != acentric.size()){ + if (static_cast(mole_fractions.size()) != acentric.size()){ throw teqp::InvalidArgument("The mole fractions should be of of size "+ std::to_string(acentric.size())); } @@ -74,10 +74,10 @@ class LKPMix{ auto omega_mix = summer_omega; auto vc_mix = summer_vcmix; auto Tc_mix = 1.0/pow(summer_vcmix, 0.25)*summer_Tcmix; - auto pc_mix = (0.2905-0.085*omega_mix)*Ru*Tc_mix/vc_mix; - auto Zc = pc_mix*vc_mix/(Ru*Tc_mix); - auto tau = Tc_mix/T; - auto delta = vc_mix*rhomolar; +// auto pc_mix = (0.2905-0.085*omega_mix)*Ru*Tc_mix/vc_mix; + auto Zc = forceeval(0.2905-0.085*omega_mix); + auto tau = forceeval(Tc_mix/T); + auto delta = forceeval(vc_mix*rhomolar); auto retval = (1.0-omega_mix/ref.omega)*alphar_func(tau, delta, Zc, simple) + (omega_mix/ref.omega)*alphar_func(tau, delta, Zc, ref); return forceeval(retval); @@ -85,7 +85,7 @@ class LKPMix{ }; // Factory function takes in JSON data and returns an instance of the LKPMix class -auto make_LKPMix(const nlohmann::json& j){ +inline auto make_LKPMix(const nlohmann::json& j){ return LKPMix(j.at("Tcrit / K"), j.at("pcrit / Pa"), j.at("acentric"), j.at("R / J/mol/K"), j.at("kmat")); } diff --git a/include/teqp/models/fwd.hpp b/include/teqp/models/fwd.hpp index b2b5119f..a72c8308 100644 --- a/include/teqp/models/fwd.hpp +++ b/include/teqp/models/fwd.hpp @@ -27,6 +27,7 @@ #include "teqp/models/mie/lennardjones.hpp" #include "teqp/models/mie/mie.hpp" #include "teqp/models/GERG/GERG.hpp" +#include "teqp/models/LKP.hpp" namespace teqp { diff --git a/interface/CPP/teqp_impl_factory.cpp b/interface/CPP/teqp_impl_factory.cpp index 98141a14..57d5d7f4 100644 --- a/interface/CPP/teqp_impl_factory.cpp +++ b/interface/CPP/teqp_impl_factory.cpp @@ -29,6 +29,8 @@ namespace teqp { {"CPA", [](const nlohmann::json& spec){ return make_owned(CPA::CPAfactory(spec));}}, {"PCSAFT", [](const nlohmann::json& spec){ return make_owned(PCSAFT::PCSAFTfactory(spec));}}, + {"LKP", [](const nlohmann::json& spec){ return make_owned(LKP::make_LKPMix(spec));}}, + {"multifluid", [](const nlohmann::json& spec){ return make_owned(multifluidfactory(spec));}}, {"multifluid-ECS-HuberEly1994", [](const nlohmann::json& spec){ return make_owned(ECSHuberEly::ECSHuberEly1994(spec));}}, {"SW_EspindolaHeredia2009", [](const nlohmann::json& spec){ return make_owned(squarewell::EspindolaHeredia2009(spec.at("lambda")));}}, diff --git a/src/tests/catch_test_LKP.cxx b/src/tests/catch_test_LKP.cxx index 5042e51c..864f8639 100644 --- a/src/tests/catch_test_LKP.cxx +++ b/src/tests/catch_test_LKP.cxx @@ -5,6 +5,7 @@ using Catch::Approx; #include "teqp/derivs.hpp" #include "teqp/constants.hpp" #include "teqp/models/LKP.hpp" +#include "teqp/cpp/teqpcpp.hpp" TEST_CASE("Check LKP", "[LKP]"){ using namespace teqp::LKP; @@ -19,7 +20,28 @@ TEST_CASE("Check LKP", "[LKP]"){ auto model = LKPMix(Tc_K, pc_Pa, acentric, R, kmat); auto z = (Eigen::ArrayXd(1) << 1.0).finished(); - teqp::TDXDerivatives::get_Ar00(model, 300.0, 8000.0, z); +// teqp::TDXDerivatives::get_Ar00(model, 300.0, 8000.0, z); + + nlohmann::json modelspec{ + {"Tcrit / K", Tc_K}, + {"pcrit / Pa", pc_Pa}, + {"acentric", acentric}, + {"R / J/mol/K", R}, + {"kmat", kmat} + }; +// std::cout << spec.dump(2) << std::endl; + + CHECK_NOTHROW(make_LKPMix(modelspec)); + nlohmann::json badspec = modelspec; + badspec["kmat"] = 4.7; + CHECK_THROWS(make_LKPMix(badspec)); + + nlohmann::json spec{ + {"kind", "LKP"}, + {"model", modelspec} + }; + CHECK_NOTHROW(teqp::cppinterface::make_model(spec)); + auto ptr = teqp::cppinterface::make_model(spec); struct Point{ double T, rhomolar, alphar_expected; }; for (auto& pt : std::vector{ @@ -54,22 +76,10 @@ TEST_CASE("Check LKP", "[LKP]"){ }){ auto alphar_actual = teqp::TDXDerivatives::get_Ar00(model, pt.T, pt.rhomolar, z); CHECK(alphar_actual == Approx(pt.alphar_expected).margin(1e-12)); + + auto alphar_actual_ptr = ptr->get_Ar00(pt.T, pt.rhomolar, z); + CHECK(alphar_actual_ptr == Approx(pt.alphar_expected).margin(1e-12)); } - - nlohmann::json spec{ - {"Tcrit / K", Tc_K}, - {"pcrit / Pa", pc_Pa}, - {"acentric", acentric}, - {"R / J/mol/K", R}, - {"kmat", kmat} - }; -// std::cout << spec.dump(2) << std::endl; - - CHECK_NOTHROW(make_LKPMix(spec)); - nlohmann::json badspec = spec; - badspec["kmat"] = 4.7; - CHECK_THROWS(make_LKPMix(badspec)); - } SECTION("methane + nitrogen mix"){