Add AbstractModel methods to access reducing functions

Some additional multi-fluid hybrid models are not yet covered. The pattern just needs to be followed to add them. Closes #158
usnistgov · Dec 12, 2024 · 5b3ec43 · 5b3ec43
1 parent ea91f3a
commit 5b3ec43
Show file tree

Hide file tree

Showing 7 changed files with 142 additions and 0 deletions.
diff --git a/include/teqp/cpp/deriv_adapter.hpp b/include/teqp/cpp/deriv_adapter.hpp
@@ -86,6 +86,16 @@ namespace internal{
     template<class T>struct tag{using type=T;};
 }
 
+template<typename T, typename U>
+concept CallableReducingDensity = requires(T t, U u) {
+    { t.get_reducing_density(u) };
+};
+
+template<typename T, typename U>
+concept CallableReducingTemperature = requires(T t, U u) {
+    { t.get_reducing_temperature(u) };
+};
+
 /**
  This class holds a const reference to a class, and exposes an interface that matches that used in AbstractModel
  
@@ -160,6 +170,25 @@ class DerivativeAdapter : public teqp::cppinterface::AbstractModel{
         return rho*rho*rho*static_cast<double>(centered_diff<3,4>(f, static_cast<my_float_t>(rho), 1e-16*static_cast<my_float_t>(rho)));
     }
 
+    virtual double get_reducing_density(const EArrayd& molefrac) const  override {
+        using Model = std::decay_t<decltype(mp.get_cref())>;
+        if constexpr(CallableReducingDensity<Model, EArrayd>){
+            return mp.get_cref().get_reducing_density(molefrac);
+        }
+        else{
+            throw teqp::NotImplementedError("Cannot call get_reducing_density of a class that doesn't define it");
+        }
+    }
+    virtual double get_reducing_temperature(const EArrayd& molefrac) const  override {
+        using Model = std::decay_t<decltype(mp.get_cref())>;
+        if constexpr(CallableReducingTemperature<Model, EArrayd>){
+            return mp.get_cref().get_reducing_temperature(molefrac);
+        }
+        else{
+            throw teqp::NotImplementedError("Cannot call get_reducing_temperature of a class that doesn't define it");
+        }
+    }
+
     // Virial derivatives
     virtual double get_B2vir(const double T, const EArrayd& z) const override {
         return VirialDerivatives<decltype(mp.get_cref()), double, EArrayd>::get_B2vir(mp.get_cref(), T, z);

diff --git a/include/teqp/cpp/teqpcpp.hpp b/include/teqp/cpp/teqpcpp.hpp
@@ -121,6 +121,11 @@ namespace teqp {
             virtual double get_Ar02ep(const double, const double, const EArrayd&) const = 0;
             virtual double get_Ar03ep(const double, const double, const EArrayd&) const = 0;
 
+            // Pass-through functions to give access to reducing function evaluations
+            // for multi-fluid models in the corresponding-states formulation
+            virtual double get_reducing_density(const EArrayd&) const = 0;
+            virtual double get_reducing_temperature(const EArrayd&) const = 0;
+
             // Virial derivatives
             virtual double get_B2vir(const double T, const EArrayd& z) const = 0;
             virtual std::map<int, double> get_Bnvir(const int Nderiv, const double T, const EArrayd& z) const = 0;

diff --git a/include/teqp/models/GERG/GERG.hpp b/include/teqp/models/GERG/GERG.hpp
@@ -922,6 +922,17 @@ class GERG2004ResidualModel{
         auto val = forceeval(corr.alphar(tau, delta, molefrac) + dep.alphar(tau, delta, molefrac));
         return val;
     }
+
+    template<typename MoleFracType>
+    auto get_reducing_temperature(const MoleFracType& molefrac) const {
+        return forceeval(red.get_Tr(molefrac));
+    }
+
+    template<typename MoleFracType>
+    auto get_reducing_density(const MoleFracType& molefrac) const {
+        return forceeval(red.get_rhor(molefrac));
+    }
+
 };
 
 
@@ -1178,6 +1189,16 @@ class GERG2008ResidualModel{
         auto val = forceeval(corr.alphar(tau, delta, molefrac) + dep.alphar(tau, delta, molefrac));
         return val;
     }
+
+    template<typename MoleFracType>
+    auto get_reducing_temperature(const MoleFracType& molefrac) const {
+        return forceeval(red.get_Tr(molefrac));
+    }
+
+    template<typename MoleFracType>
+    auto get_reducing_density(const MoleFracType& molefrac) const {
+        return forceeval(red.get_rhor(molefrac));
+    }
 };
 
 class GERG2008IdealGasModel{

diff --git a/include/teqp/models/ammonia_water.hpp b/include/teqp/models/ammonia_water.hpp
@@ -57,6 +57,10 @@ namespace teqp{
             auto Tred = forceeval(TcNH3*xNH3*xNH3 + TcH2O*(1-xNH3)*(1-xNH3) + 2.0*xNH3*(1.0-pow(xNH3, alpha))*k_T/2.0*(TcNH3 + TcH2O));
             return Tred;
         }
+        template<typename MoleFracType>
+        auto get_reducing_temperature(const MoleFracType& molefrac) const {
+            return get_Treducing(molefrac);
+        }
 
         template<typename MoleFracType>
         auto get_rhoreducing(const MoleFracType& molefrac) const {
@@ -71,6 +75,10 @@ namespace teqp{
             auto vred = forceeval(vcNH3*xNH3*xNH3 + vcH2O*(1-xNH3)*(1-xNH3) + 2.0*xNH3*(1.0-pow(xNH3, beta))*k_V/2.0*(vcNH3 + vcH2O));
             return forceeval(1/vred);
         }
+        template<typename MoleFracType>
+        auto get_reducing_density(const MoleFracType& molefrac) const {
+            return get_rhoreducing(molefrac);
+        }
 
         template<typename TType, typename RhoType, typename MoleFracType>
         auto alphar(const TType& T,

diff --git a/include/teqp/models/multifluid.hpp b/include/teqp/models/multifluid.hpp
@@ -195,6 +195,16 @@ class MultiFluid {
     {
         return corr.alphari(tau, delta, i);
     }
+
+    template<typename MoleFracType>
+    auto get_reducing_temperature(const MoleFracType& molefrac) const {
+        return forceeval(redfunc.get_Tr(molefrac));
+    }
+
+    template<typename MoleFracType>
+    auto get_reducing_density(const MoleFracType& molefrac) const {
+        return forceeval(redfunc.get_rhor(molefrac));
+    }
 };
 
 

diff --git a/interface/pybind11_wrapper.cpp b/interface/pybind11_wrapper.cpp
@@ -497,6 +497,9 @@ void init_teqp(py::module& m) {
 
         .def("get_R", &am::get_R, "molefrac"_a.noconvert())
 
+        .def("get_reducing_density", &am::get_reducing_density,  "molefrac"_a.noconvert())
+        .def("get_reducing_temperature", &am::get_reducing_temperature,  "molefrac"_a.noconvert())
+
         .def("get_B2vir", &am::get_B2vir, "T"_a, "molefrac"_a.noconvert())
         .def("get_Bnvir", &am::get_Bnvir, "Nderiv"_a, "T"_a, "molefrac"_a.noconvert())
         .def("get_dmBnvirdTm", &am::get_dmBnvirdTm, "Nderiv"_a, "NTderiv"_a, "T"_a, "molefrac"_a.noconvert())

diff --git a/src/tests/catch_test_modelset.cxx b/src/tests/catch_test_modelset.cxx
@@ -11,6 +11,8 @@ using Catch::Matchers::WithinAbs;
 
 #include "catch_fixtures.hpp"
 
+#include "test_common.in"
+
 auto LKPmethane = [](){
     // methane, check values from TREND
     std::vector<double> Tc_K = {190.564};
@@ -146,3 +148,67 @@ TEST_CASE("virials", "[virials]"){
         }
     }
 }
+
+auto GERG2004metheth_ = [](){
+    return R"({"kind":"GERG2004resid", "model":{"names": ["methane", "ethane"]}} )"_json;
+};
+auto GERG2008metheth_ = [](){
+    return R"({"kind":"GERG2008resid", "model":{"names": ["methane", "ethane"]}} )"_json;
+};
+auto PCSAFTmetheth_ = [](){
+    return  R"({
+      "kind": "PCSAFT",
+      "model": {
+          "names": ["Methane", "Ethane"]
+      }
+    })"_json;
+};
+
+auto multifluidmetheth_ = [](){
+    auto j = R"({
+      "kind": "multifluid",
+      "model": {
+          "components": ["Methane", "Ethane"]
+      }
+    })"_json;
+    j["model"]["root"] = FLUIDDATAPATH;
+    j["model"]["BIP"] = FLUIDDATAPATH+"/dev/mixtures/mixture_binary_pairs.json";
+    j["model"]["departure"] = FLUIDDATAPATH+"/dev/mixtures/mixture_departure_functions.json";
+    return j;
+};
+
+auto AmmoniaWaterTillerRoth_ = [](){
+    return R"({ "kind": "AmmoniaWaterTillnerRoth", "model": {} })"_json;
+};
+
+
+// Add a little static_assert to make sure the concept is working properly
+#include "teqp/cpp/deriv_adapter.hpp"
+#include "teqp/models/GERG/GERG.hpp"
+static_assert(teqp::cppinterface::adapter::CallableReducingTemperature<teqp::GERG2008::GERG2008ResidualModel, Eigen::ArrayXd>);
+static_assert(teqp::cppinterface::adapter::CallableReducingDensity<teqp::GERG2008::GERG2008ResidualModel, Eigen::ArrayXd>);
+
+// Binary multi-fluid models
+std::map<std::string, std::pair<nlohmann::json, std::vector<evalpoint>>> MultifluidBinaryTestSet = {
+    {"GERG2004", {GERG2004metheth_(), {}}},
+    {"GERG2008", {GERG2008metheth_(), {}}},
+    {"multifluid", {multifluidmetheth_(), {}}},
+    {"AmmoniaWater", {AmmoniaWaterTillerRoth_(), {}}},
+};
+
+TEST_CASE("reducing", "[MFreducing]"){
+    Eigen::ArrayXd z(2); z(0) = 0.4; z(1) = 0.6;
+    for (const auto& [kind, specdata] : MultifluidBinaryTestSet){
+        auto [spec, points] = specdata;
+        auto model = teqp::cppinterface::make_model(spec);
+
+        CAPTURE(kind);
+        CHECK_NOTHROW(model->get_reducing_density(z));
+        CHECK_NOTHROW(model->get_reducing_temperature(z));
+    }
+
+    SECTION("PCSAFT"){ // this is test that SHOULD throw, to prove the negative
+        auto model = teqp::cppinterface::make_model(PCSAFTmetheth_());
+        CHECK_THROWS(model->get_reducing_density(z));
+    }
+}