From cf85ee074b93f957e16616349e1b83852a81d0dd Mon Sep 17 00:00:00 2001
From: Ian Bell <ian.bell@nist.gov>
Date: Wed, 2 Jun 2021 19:32:48 -0400
Subject: [PATCH] Possible now also to modify Fij (though you can't fiddle with
 the departure function)

---
 include/teqp/models/multifluid.hpp | 59 ++++++++++++++++++++++++++++++
 interface/multifluid_mutant.cpp    | 20 +++++++---
 src/multifluid.cpp                 |  5 ++-
 3 files changed, 78 insertions(+), 6 deletions(-)

diff --git a/include/teqp/models/multifluid.hpp b/include/teqp/models/multifluid.hpp
index b003566..23eb0dd 100644
--- a/include/teqp/models/multifluid.hpp
+++ b/include/teqp/models/multifluid.hpp
@@ -704,12 +704,71 @@ auto build_BIPmodified(Model& model, const nlohmann::json& j) {
     return MultiFluidReducingFunctionAdapter(model, std::move(newred));
 }
 
+/**
+This class holds a lightweight reference to the core parts of the model
+
+The reducing and departure functions are moved into this class, while the donor class is used for the corresponding states portion
+*/
+template<typename ReducingFunction, typename DepartureFunction, typename BaseClass>
+class MultiFluidAdapter {
+
+public:
+    const BaseClass& base;
+    const ReducingFunction redfunc;
+    const DepartureFunction depfunc;
 
+    template<class VecType>
+    auto R(const VecType& molefrac) const { return base.R(molefrac); }
 
+    MultiFluidAdapter(const BaseClass& base, ReducingFunction&& redfunc, DepartureFunction &&depfunc) : base(base), redfunc(redfunc), depfunc(depfunc) {};
+
+    template<typename TType, typename RhoType, typename MoleFracType>
+    auto alphar(const TType& T,
+        const RhoType& rho,
+        const MoleFracType& molefrac) const
+    {
+        auto Tred = forceeval(redfunc.get_Tr(molefrac));
+        auto rhored = forceeval(redfunc.get_rhor(molefrac));
+        auto delta = forceeval(rho / rhored);
+        auto tau = forceeval(Tred / T);
+        auto val = base.corr.alphar(tau, delta, molefrac) + depfunc.alphar(tau, delta, molefrac);
+        return forceeval(val);
+    }
+};
 
+template<class Model>
+auto build_multifluid_mutant(Model& model, const nlohmann::json& j) {
 
+    auto red = model.redfunc;
+    auto betaT = red.betaT;
+    betaT(0, 1) = j["betaT"];
+    betaT(1, 0) = 1 / betaT(0, 1);
+    auto betaV = red.betaV;
+    betaV(0, 1) = j["betaV"];
+    betaV(1, 0) = 1 / betaV(0, 1);
+    auto gammaT = red.gammaT, gammaV = red.gammaV;
+    gammaT(0, 1) = j["gammaT"]; gammaT(1, 0) = gammaT(0, 1);
+    gammaV(0, 1) = j["gammaV"]; gammaV(1, 0) = gammaV(0, 1);
+    auto Tc = red.Tc, vc = red.vc;
+    auto newred = MultiFluidReducingFunction(betaT, gammaT, betaV, gammaV, Tc, vc);
 
+    if (j.contains("Fij") && j["Fij"] != 0.0) {
+        throw std::invalid_argument("Don't support Fij != 0 for now");
+    }
+    auto N = 2;
+    // Allocate the matrix with default models
+    Eigen::MatrixXd F(2, 2); F.setZero();
+    std::vector<std::vector<MultiFluidDepartureFunction>> funcs(N); 
+    for (auto i = 0; i < funcs.size(); ++i) {
+        funcs[i].resize(funcs.size());
+        for (auto j = 0; j < N; ++j) {
+            funcs[i][j].set_type("none");
+        }
+    }
+    auto newdep = DepartureContribution(std::move(F), std::move(funcs));
 
+    return MultiFluidAdapter(model, std::move(newred), std::move(newdep));
+}
 
 
 class DummyEOS {
diff --git a/interface/multifluid_mutant.cpp b/interface/multifluid_mutant.cpp
index 26ce267..b665bc1 100644
--- a/interface/multifluid_mutant.cpp
+++ b/interface/multifluid_mutant.cpp
@@ -6,9 +6,19 @@
 void add_multifluid_mutant(py::module& m) {
 
     using MultiFluid = decltype(build_multifluid_model(std::vector<std::string>{"", ""}, "", ""));
-    m.def("build_BIPmodified", &build_BIPmodified<MultiFluid>);
-    using RedType = std::decay_t<decltype(MultiFluid::redfunc)>; 
-    using BIPmod = MultiFluidReducingFunctionAdapter<RedType, MultiFluid>;
-    auto wMFBIP = py::class_<BIPmod>(m, "MultiFluidBIP");
-    add_derivatives<BIPmod>(m, wMFBIP);
+    {
+        m.def("build_BIPmodified", &build_BIPmodified<MultiFluid>);
+        using RedType = std::decay_t<decltype(MultiFluid::redfunc)>;
+        using BIPmod = MultiFluidReducingFunctionAdapter<RedType, MultiFluid>;
+        auto wMFBIP = py::class_<BIPmod>(m, "MultiFluidBIP");
+        add_derivatives<BIPmod>(m, wMFBIP);
+    }
+    {
+        m.def("build_multifluid_mutant", &build_multifluid_mutant<MultiFluid>);
+        using RedType = std::decay_t<decltype(MultiFluid::redfunc)>;
+        using DepType = std::decay_t<decltype(MultiFluid::dep)>;
+        using BIPmod = MultiFluidAdapter<RedType, DepType, MultiFluid>;
+        auto wMFBIP = py::class_<BIPmod>(m, "MultiFluidMutant");
+        add_derivatives<BIPmod>(m, wMFBIP);
+    }
 }
\ No newline at end of file
diff --git a/src/multifluid.cpp b/src/multifluid.cpp
index f532071..c92684b 100644
--- a/src/multifluid.cpp
+++ b/src/multifluid.cpp
@@ -159,8 +159,11 @@ int main(){
 
     //time_calls(coolprop_root, BIPcollection);
     {
-        nlohmann::json flags = { {"estimate", true} };
+        nlohmann::json flags = { {"estimate", true},{"another","key"} };
         auto model = build_multifluid_model({ "Ethane", "R1234ze(E)" }, coolprop_root, BIPcollection, flags);
+
+        nlohmann::json j = { {"betaT", 1.0},{"gammaT", 1.0},{"betaV", 1.0},{"gammaV", 1.0},{"Fij", 0.0} };
+        auto mutant = build_mutant(model, j);
     }
 {
     auto model = build_multifluid_model({ "methane", "ethane" }, coolprop_root, BIPcollection);
-- 
GitLab