diff --git a/include/teqp/models/CPA.hpp b/include/teqp/models/CPA.hpp index f025dba8..1fca8adf 100644 --- a/include/teqp/models/CPA.hpp +++ b/include/teqp/models/CPA.hpp @@ -311,7 +311,7 @@ inline auto CPAfactory(const nlohmann::json &j){ auto build_assoc_pure = [](const auto& j) -> AssociationVariantWrapper{ auto N = j["pures"].size(); - if (N == 1 && j.at("pures")[0].contains("class") ){ + if (N == 1 && j.at("pures").at(0).contains("class") ){ // This is the backwards compatible approach // with the old style of defining the association class {1,2B...} std::vector classes; diff --git a/include/teqp/models/association/association.hpp b/include/teqp/models/association/association.hpp index 7872a59a..28849eac 100644 --- a/include/teqp/models/association/association.hpp +++ b/include/teqp/models/association/association.hpp @@ -28,7 +28,7 @@ struct AssociationOptions{ association::radial_dist radial_dist; double alpha = 0.5; double rtol = 1e-12, atol = 1e-12; - int max_iters = 10; + int max_iters = 100; }; inline void from_json(const nlohmann::json& j, AssociationOptions& o) { if (j.contains("alpha")){ j.at("alpha").get_to(o.alpha); } diff --git a/src/tests/catch_test_association.cxx b/src/tests/catch_test_association.cxx index a4466dc5..1ce6af37 100644 --- a/src/tests/catch_test_association.cxx +++ b/src/tests/catch_test_association.cxx @@ -1,7 +1,7 @@ #include -#include +#include #include -using Catch::Approx; +using Catch::Matchers::WithinRel; #include @@ -29,6 +29,8 @@ TEST_CASE("Test ethanol + water association", "[association]"){ std::vector> molecule_sites = {{"e", "H"}, {"e", "e", "H", "H"}}; association::AssociationOptions opt; + opt.radial_dist = association::radial_dist::CS; + opt.max_iters = 1000; opt.interaction_partners = {{"e", {"H",}}, {"H", {"e",}}}; association::Association a(b_m3mol, beta, epsilon_Jmol, molecule_sites, opt); @@ -37,20 +39,67 @@ TEST_CASE("Test ethanol + water association", "[association]"){ auto Ngroups = a.mapper.to_siteid.size(); Eigen::ArrayXd X_init = Eigen::ArrayXd::Ones(Ngroups); - double v = 3.0680691201961814e-5; - Eigen::ArrayXd X_A = a.successive_substitution(303.15, 1/v, molefracs, X_init); - CHECK(X_A[0] == Approx(0.06258400385436955)); - CHECK(X_A[3] == Approx(0.10938445109190545)); + double v = 3.0680691201961814e-5, T = 303.15; + auto Delta = a.get_Delta(T, 1/v, molefracs); + CAPTURE(Delta); + CHECK_THAT(Delta(0,0), WithinRel(5.85623687e-27, 1e-8)); + CHECK_THAT(Delta(0,3), WithinRel(4.26510827e-27, 1e-8)); + CHECK_THAT(Delta(3,3), WithinRel(2.18581242e-27, 1e-8)); + Eigen::ArrayXd X_A = a.successive_substitution(T, 1/v, molefracs, X_init); + CHECK_THAT(X_A[0], WithinRel(0.06258400385436955, 1e-8)); + CHECK_THAT(X_A[3], WithinRel(0.10938445109190545, 1e-8)); BENCHMARK("SS"){ - return a.successive_substitution(303.15, 1/v, molefracs, X_init); + return a.successive_substitution(T, 1/v, molefracs, X_init); }; BENCHMARK("alphar"){ - return a.alphar(303.15, 1/v, molefracs); + return a.alphar(T, 1/v, molefracs); }; } -TEST_CASE("Ethanol + water with CPA", "[association]"){ +TEST_CASE("Ethanol with CPA and old class names", "[association]"){ + nlohmann::json ethanol = { + {"a0i / Pa m^6/mol^2", 0.85164}, {"bi / m^3/mol", 0.0491e-3}, {"c1", 0.7502}, {"Tc / K", 513.92}, + {"epsABi / J/mol", 21500.0}, {"betaABi", 0.008}, {"class", "2B"} + }; + nlohmann::json jCPA = { + {"cubic", "SRK"}, + {"radial_dist", "CS"}, + {"pures", {ethanol}}, + {"R_gas / J/mol/K", 8.31446261815324} + }; + nlohmann::json j = { + {"kind", "CPA"}, + {"model", jCPA} + }; + CHECK_NOTHROW(teqp::cppinterface::make_model(j, false)); +} + +TEST_CASE("Ethanol + water with CPA and old class names", "[association]"){ + nlohmann::json water = { + {"a0i / Pa m^6/mol^2", 0.12277}, {"bi / m^3/mol", 0.0000145}, {"c1", 0.6736}, {"Tc / K", 647.13}, + {"epsABi / J/mol", 16655.0}, {"betaABi", 0.0692}, {"sites", {"e","e","H","H"}} + }; + nlohmann::json ethanol = { + {"a0i / Pa m^6/mol^2", 0.85164}, {"bi / m^3/mol", 0.0491e-3}, {"c1", 0.7502}, {"Tc / K", 513.92}, + {"epsABi / J/mol", 21500.0}, {"betaABi", 0.008}, {"sites", {"e","H"}} + }; + nlohmann::json jCPA = { + {"cubic", "SRK"}, + {"radial_dist", "CS"}, +// {"combining", "CR1"}, + {"pures", {ethanol, water}}, + {"R_gas / J/mol/K", 8.31446261815324} + }; + nlohmann::json j = { + {"kind", "CPA"}, + {"validate", false}, + {"model", jCPA} + }; + CHECK_NOTHROW(teqp::cppinterface::make_model(j, false)); +} + +TEST_CASE("Ethanol + water with CPA against Clapeyron.jl", "[association]"){ nlohmann::json water = { {"a0i / Pa m^6/mol^2", 0.12277}, {"bi / m^3/mol", 0.0000145}, {"c1", 0.6736}, {"Tc / K", 647.13}, {"epsABi / J/mol", 16655.0}, {"betaABi", 0.0692}, {"sites", {"e","e","H","H"}} @@ -79,10 +128,10 @@ TEST_CASE("Ethanol + water with CPA", "[association]"){ double R = model->get_R(molefracs); auto alphar = model->get_Ar00(303.15, rhomolar, molefracs); - CHECK(alphar == Approx(-8.333844120879878)); + CHECK_THAT(alphar, WithinRel(-8.333844120879878, 1e-8)); double p = T*R*rhomolar*(1+model->get_Ar01(T, rhomolar, molefracs)); - CHECK(p == Approx(1e5)); + CHECK_THAT(p, WithinRel(1e5, 1e-6)); BENCHMARK("p(T,rho)"){ return T*R*rhomolar*(1+model->get_Ar01(T, rhomolar, molefracs)); };