Merge pull request #79 from usnistgov/GERG

Add residual parts of GERG-2004 and GERG-2008

include/teqp/models/fwd.hpp

@@ -26,6 +26,7 @@
 #include "teqp/models/model_potentials/2center_ljf.hpp"
 #include "teqp/models/mie/lennardjones.hpp"
 #include "teqp/models/mie/mie.hpp"
+#include "teqp/models/GERG/GERG.hpp"
 
 namespace teqp {
 

interface/CPP/teqp_impl_factory.cpp

@@ -42,6 +42,9 @@ namespace teqp {
             {"2CLJF-Quadrupole", [](const nlohmann::json& spec){ return make_owned(twocenterljf::build_two_center_model_quadrupole(spec.at("author"), spec.at("L^*"), spec.at("(Q^*)^2")));}},
             {"IdealHelmholtz", [](const nlohmann::json& spec){ return make_owned(IdealHelmholtz(spec));}},
 
+            {"GERG2004resid", [](const nlohmann::json& spec){ return make_owned(GERG2004::GERG2004ResidualModel(spec.at("names")));}},
+            {"GERG2008resid", [](const nlohmann::json& spec){ return make_owned(GERG2008::GERG2008ResidualModel(spec.at("names")));}},
+
             // Implemented in its own compilation unit to help with compilation time
             {"SAFT-VR-Mie", [](const nlohmann::json& spec){ return make_SAFTVRMie(spec); }}
         };

interface/pybind11_wrapper.cpp

@@ -95,6 +95,18 @@ void attach_multifluid_methods(py::object&obj){
     setattr("set_meta", MethodType(py::cpp_function([](py::object& o, const std::string& s){ return get_mutable_typed<TYPE>(o).set_meta(s); }, "self"_a, "s"_a), obj));
     setattr("get_alpharij", MethodType(py::cpp_function([](py::object& o, const int i, const int j, const double tau, const double delta){ return get_typed<TYPE>(o).dep.get_alpharij(i,j,tau,delta); }, "self"_a, "i"_a, "j"_a, "tau"_a, "delta"_a), obj));
 }
+template<typename TYPE>
+void attach_GERG_methods(py::object&obj){
+    auto setattr = py::getattr(obj, "__setattr__");
+    auto MethodType = py::module_::import("types").attr("MethodType");
+    setattr("get_Tcvec", MethodType(py::cpp_function([](py::object& o){ return get_typed<TYPE>(o).red.get_Tcvec(); }), obj));
+    setattr("get_vcvec", MethodType(py::cpp_function([](py::object& o){ return get_typed<TYPE>(o).red.get_vcvec(); }), obj));
+    setattr("get_Tr", MethodType(py::cpp_function([](py::object& o, REArrayd& molefrac){ return get_typed<TYPE>(o).red.get_Tr(molefrac); }, "self"_a, "molefrac"_a.noconvert()), obj));
+    setattr("get_rhor", MethodType(py::cpp_function([](py::object& o, REArrayd& molefrac){ return get_typed<TYPE>(o).red.get_rhor(molefrac); }, "self"_a, "molefrac"_a.noconvert()), obj));
+//    setattr("get_meta", MethodType(py::cpp_function([](py::object& o){ return get_typed<TYPE>(o).get_meta(); }), obj));
+//    setattr("set_meta", MethodType(py::cpp_function([](py::object& o, const std::string& s){ return get_mutable_typed<TYPE>(o).set_meta(s); }, "self"_a, "s"_a), obj));
+//    setattr("get_alpharij", MethodType(py::cpp_function([](py::object& o, const int i, const int j, const double tau, const double delta){ return get_typed<TYPE>(o).dep.get_alpharij(i,j,tau,delta); }, "self"_a, "i"_a, "j"_a, "tau"_a, "delta"_a), obj));
+}
 
 // Type index variables matching the model types, used for runtime attachment of model-specific methods
 const std::type_index vdWEOS1_i{std::type_index(typeid(vdWEOS1))};
@@ -111,6 +123,8 @@ const std::type_index twocenterLJF_i{std::type_index(typeid(twocenterLJF_t))};
 const std::type_index QuantumPR_i{std::type_index(typeid(QuantumPR_t))};
 const std::type_index advancedPRaEres_i{std::type_index(typeid(advancedPRaEres_t))};
 const std::type_index RKPRCismondi2005_i{std::type_index(typeid(RKPRCismondi2005_t))};
+const std::type_index GERG2004ResidualModel_i{std::type_index(typeid(GERG2004::GERG2004ResidualModel))};
+const std::type_index GERG2008ResidualModel_i{std::type_index(typeid(GERG2008::GERG2008ResidualModel))};
 
 /**
  At runtime we can add additional model-specific methods that only apply for a particular model.  We take in a Python-wrapped
@@ -155,10 +169,12 @@ void attach_model_specific_methods(py::object& obj){
         setattr("get_epsilon_over_k_K", MethodType(py::cpp_function([](py::object& o){ return get_typed<SAFTVRMie_t>(o).get_epsilon_over_k_K(); }), obj));
         setattr("get_lambda_r", MethodType(py::cpp_function([](py::object& o){ return get_typed<SAFTVRMie_t>(o).get_lambda_r(); }), obj));
         setattr("get_lambda_a", MethodType(py::cpp_function([](py::object& o){ return get_typed<SAFTVRMie_t>(o).get_lambda_a(); }), obj));
+
 //        setattr("max_rhoN", MethodType(py::cpp_function([](py::object& o, double T, REArrayd& molefrac){ return get_typed<SAFTVRMie_t>(o).max_rhoN(T, molefrac); }, "self"_a, "T"_a, "molefrac"_a), obj));
         setattr("get_kmat", MethodType(py::cpp_function([](py::object& o){ return get_typed<SAFTVRMie_t>(o).get_kmat(); }), obj));
         setattr("has_polar", MethodType(py::cpp_function([](py::object& o){ return get_typed<SAFTVRMie_t>(o).has_polar(); }), obj));
         setattr("get_core_calcs", MethodType(py::cpp_function([](py::object& o, double T, double rhomolar, REArrayd& molefrac){ return get_typed<SAFTVRMie_t>(o).get_core_calcs(T, rhomolar, molefrac); }, "self"_a, "T"_a, "rhomolar"_a, "molefrac"_a), obj));
+
     }
     else if (index == canonical_cubic_i){
         setattr("get_a", MethodType(py::cpp_function([](py::object& o, double T, REArrayd& molefrac){ return get_typed<canonical_cubic_t>(o).get_a(T, molefrac); }, "self"_a, "T"_a, "molefrac"_a), obj));
@@ -230,6 +246,12 @@ void attach_model_specific_methods(py::object& obj){
             return teqp::MultiFluidVLEAncillaries(jancillaries);
         }, "self"_a, py::arg_v("i", std::nullopt, "None")), obj));
     }
+    else if (index == GERG2004ResidualModel_i){
+        attach_GERG_methods<GERG2004::GERG2004ResidualModel>(obj);
+    }
+    else if (index == GERG2008ResidualModel_i){
+        attach_GERG_methods<GERG2008::GERG2008ResidualModel>(obj);
+    }
     else if (index == idealgas_i){
         // Here X-Macros are used to create functions like get_Aig00, get_Aig01, ....
         #define X(i,j) setattr(stringify(get_Aig ## i ## j), MethodType(py::cpp_function([](py::object& o, double T, double rho, REArrayd& molefrac){ \

src/bench_AbstractModel.cpp

@@ -6,9 +6,23 @@
 #include "teqp/derivs.hpp"
 #include "teqp/cpp/deriv_adapter.hpp"
 #include "teqp/models/multifluid.hpp"
+#include "teqp/models/GERG/GERG.hpp"
 
 using namespace teqp;
 
+TEST_CASE("GERG2008 parts", "[GERG2008]")
+{
+    BENCHMARK("bg"){
+        return GERG2008::get_betasgammas("methane","n-decane");
+    };
+    BENCHMARK("bg backwards"){
+        return GERG2008::get_betasgammas("n-decane", "methane");
+    };
+    BENCHMARK("bg fallback to 2004"){
+        return GERG2008::get_betasgammas("ethane", "methane");
+    };
+}
+
 TEST_CASE("multifluid derivatives", "[mf]")
 {
     nlohmann::json j = {
@@ -25,6 +39,9 @@ TEST_CASE("multifluid derivatives", "[mf]")
     auto z = (Eigen::ArrayXd(1) << 1.0).finished();
     auto rhovec = 300.0* z;
 
+    BENCHMARK("construction"){
+        return teqp::cppinterface::make_model(j);
+    };
     BENCHMARK("alphar") {
         return am->get_Arxy(0, 0, 300, 3.0, z);
     };
@@ -84,3 +101,51 @@ TEST_CASE("multifluid derivatives via DerivativeAdapter", "[mf]")
         return view(model)->get_Bnvir(4, 300, z);
     };
 }
+
+
+TEST_CASE("GERG2008 derivatives", "[GERG2008]")
+{
+    nlohmann::json j = {
+        {"kind", "GERG2008resid"},
+        {"model", {
+            {"names", {"methane","ethane","propane","n-butane"}},
+        }
+    }};
+    //std::cout << j.dump(2);
+
+
+    auto z = (Eigen::ArrayXd(1) << 1.0).finished();
+    auto rhovec = 300.0* z;
+    auto am = teqp::cppinterface::make_model(j);
+
+    BENCHMARK("construction"){
+        return teqp::cppinterface::make_model(j);
+    };
+    BENCHMARK("alphar") {
+        return am->get_Arxy(0, 0, 300, 3.0, z);
+    };
+    BENCHMARK("Ar20") {
+        return am->get_Ar20(300, 3.0, z);
+    };
+    BENCHMARK("get_Ar02n") {
+        return am->get_Ar02n(300, 3.0, z);
+    };
+    BENCHMARK("fugacity coefficients") {
+        return am->get_fugacity_coefficients(300.0, rhovec);
+    };
+    BENCHMARK("cvr/R") {
+        return -1*am->get_Arxy(2, 0, 300, 3.0, z);
+    };
+    BENCHMARK("partial_molar_volumes") {
+        return am->get_partial_molar_volumes(300.0, rhovec);
+    };
+    BENCHMARK("get_deriv_mat2") {
+        return am->get_deriv_mat2(300.0, 3.0, z);
+    };
+    BENCHMARK("build_iteration_Jv") {
+        auto mat = am->get_deriv_mat2(300.0, 3.0, z);
+        auto mat2 = am->get_deriv_mat2(300.0, 3.0, z);
+        const std::vector<char> vars = {'T','D','P','S'};
+        return teqp::cppinterface::build_iteration_Jv(vars, mat, mat2, 8.3144, 300.0, 300.0, z);
+    };
+}