First commits to the multifluid + AC modeling approach

Seems to be working properly, notebook looks decent

doc/source/models/index.rst

@@ -28,5 +28,6 @@ Models
    GERG
    ECS
    MultifluidPlusAssociation
+   MultifluidPlusAC
 
    IdealGas

include/teqp/models/activity/activity_models.hpp

@@ -0,0 +1,120 @@
+#pragma once
+
+
+namespace teqp::activity::activity_models{
+
+/**
+ A residual Helmholtz term that returns nothing (the empty term)
+ */
+template<typename NumType>
+class NullResidualHelmholtzOverRT {
+public:
+    template<typename TType, typename MoleFractions>
+    auto operator () (const TType& /*T*/, const MoleFractions& molefracs) const {
+        std::common_type_t<TType, decltype(molefracs[0])> val = 0.0;
+        return val;
+    }
+};
+
+/**
+ 
+ \f[
+ \frac{a^{E,\gamma}_{total}}{RT} = -sum_iz_i\ln\left(\sum_jz_jOmega_{ij}(T)\right)
+ \f]
+ 
+ \f[
+ \frac{a^{E,\gamma}_{comb}}{RT} = -sum_iz_i\ln\left(\frac{\Omega_i}{z_i}\right)
+  \f]
+ 
+ \f[
+ \frac{a^{E,\gamma}_{res}}{RT} = \frac{a^{E,\gamma}_{total}}{RT} - \frac{a^{E,\gamma}_{comb}}{RT}
+ \f]
+ 
+ Volume fraction of component \f$i\f$
+\f[
+ \phi_i = \frac{z_iv_i}{\sum_j z_j v_j}
+ \f]
+ with \f$v_i = b_i\f$
+ */
+template<typename NumType>
+class WilsonResidualHelmholtzOverRT {
+
+public:
+    const std::vector<double> b;
+    const Eigen::ArrayXXd m, n;
+    WilsonResidualHelmholtzOverRT(const std::vector<double>& b, const Eigen::ArrayXXd& m, const Eigen::ArrayXXd& n) : b(b), m(m), n(n) {};
+
+    template<typename TType, typename MoleFractions>
+    auto combinatorial(const TType& /*T*/, const MoleFractions& molefracs) const {
+        if (b.size() != static_cast<std::size_t>(molefracs.size())){
+            throw teqp::InvalidArgument("Bad size of molefracs");
+        }
+
+        using TYPE = std::common_type_t<TType, decltype(molefracs[0])>;
+        // The denominator in Phi
+        TYPE Vtot = 0.0;
+        for (auto i = 0U; i < molefracs.size(); ++i){
+            auto v_i = b[i];
+            Vtot += molefracs[i]*v_i;
+        }
+
+        TYPE summer = 0.0;
+        for (auto i = 0U; i < molefracs.size(); ++i){
+            auto v_i = b[i];
+            // The ratio phi_i/z_i is expressed like this to better handle
+            // the case of z_i = 0, which would otherwise be a divide by zero
+            // in the case that the composition of one component is zero
+            auto phi_i_over_z_i = v_i/Vtot;
+            summer += molefracs[i]*log(phi_i_over_z_i);
+        }
+        return summer;
+    }
+
+    template<typename TType>
+    auto get_Aij(std::size_t i, std::size_t j, const TType& T) const{
+        return forceeval(m(i,j)*T + n(i,j));
+    }
+
+    template<typename TType, typename MoleFractions>
+    auto total(const TType& T, const MoleFractions& molefracs) const {
+
+        using TYPE = std::common_type_t<TType, decltype(molefracs[0])>;
+        TYPE summer = 0.0;
+        for (auto i = 0U; i < molefracs.size(); ++i){
+            auto v_i = b[i];
+            TYPE summerj = 0.0;
+            for (auto j = 0U; j < molefracs.size(); ++j){
+                auto v_j = b[j];
+                auto Aij = get_Aij(i,j,T);
+                auto Omega_ji = v_j/v_i*exp(-Aij/T);
+                summerj += molefracs[j]*Omega_ji;
+            }
+            summer += molefracs[i]*log(summerj);
+        }
+        return forceeval(-summer);
+    }
+
+    // Returns ares/RT
+    template<typename TType, typename MoleFractions>
+    auto operator () (const TType& T, const MoleFractions& molefracs) const {
+        return forceeval(total(T, molefracs) - combinatorial(T, molefracs));
+    }
+};
+
+using ResidualHelmholtzOverRTVariant = std::variant<NullResidualHelmholtzOverRT<double>, WilsonResidualHelmholtzOverRT<double>>;
+
+inline ResidualHelmholtzOverRTVariant ares_model_factory(const nlohmann::json& armodel) {
+
+    std::string type = armodel.at("type");
+    if (type == "Wilson"){
+        std::vector<double> b = armodel.at("b");
+        auto mWilson = build_square_matrix(armodel.at("m"));
+        auto nWilson = build_square_matrix(armodel.at("n"));
+        return WilsonResidualHelmholtzOverRT<double>(b, mWilson, nWilson);
+    }
+    else{
+        throw teqp::InvalidArgument("bad type of ares model: " + type);
+    }
+};
+
+}

include/teqp/models/cubics.hpp

@@ -17,6 +17,8 @@ Implementations of the canonical cubic equations of state
 #include "teqp/math/pow_templates.hpp"
 
 #include "nlohmann/json.hpp"
+#include "teqp/models/activity/activity_models.hpp"
+using namespace teqp::activity::activity_models;
 
 #include <Eigen/Dense>
 
@@ -460,106 +462,6 @@ inline void from_json(const json& j, AdvancedPRaEOptions& o) {
     }
 }
 
-/**
- A residual Helmholtz term that returns nothing (the empty term)
- */
-template<typename NumType>
-class NullResidualHelmholtzOverRT {
-public:
-    template<typename TType, typename MoleFractions>
-    auto operator () (const TType& /*T*/, const MoleFractions& molefracs) const {
-        std::common_type_t<TType, decltype(molefracs[0])> val = 0.0;
-        return val;
-    }
-};
-
-/**
- 
- \f[
- \frac{a^{E,\gamma}_{total}}{RT} = -sum_iz_i\ln\left(\sum_jz_jOmega_{ij}(T)\right)
- \f]
- 
- \f[
- \frac{a^{E,\gamma}_{comb}}{RT} = -sum_iz_i\ln\left(\frac{\Omega_i}{z_i}\right)
-  \f]
- 
- \f[
- \frac{a^{E,\gamma}_{res}}{RT} = \frac{a^{E,\gamma}_{total}}{RT} - \frac{a^{E,\gamma}_{comb}}{RT}
- \f]
- 
- Volume fraction of component \f$i\f$
-\f[
- \phi_i = \frac{z_iv_i}{\sum_j z_j v_j}
- \f]
- with \f$v_i = b_i\f$
- */
-template<typename NumType>
-class WilsonResidualHelmholtzOverRT {
-
-public:
-    const std::vector<double> b;
-    const Eigen::ArrayXXd m, n;
-    WilsonResidualHelmholtzOverRT(const std::vector<double>& b, const Eigen::ArrayXXd& m, const Eigen::ArrayXXd& n) : b(b), m(m), n(n) {};
-
-    template<typename TType, typename MoleFractions>
-    auto combinatorial(const TType& /*T*/, const MoleFractions& molefracs) const {
-        if (b.size() != static_cast<std::size_t>(molefracs.size())){
-            throw teqp::InvalidArgument("Bad size of molefracs");
-        }
-
-        using TYPE = std::common_type_t<TType, decltype(molefracs[0])>;
-        // The denominator in Phi
-        TYPE Vtot = 0.0;
-        for (auto i = 0U; i < molefracs.size(); ++i){
-            auto v_i = b[i];
-            Vtot += molefracs[i]*v_i;
-        }
-
-        TYPE summer = 0.0;
-        for (auto i = 0U; i < molefracs.size(); ++i){
-            auto v_i = b[i];
-            // The ratio phi_i/z_i is expressed like this to better handle
-            // the case of z_i = 0, which would otherwise be a divide by zero
-            // in the case that the composition of one component is zero
-            auto phi_i_over_z_i = v_i/Vtot;
-            summer += molefracs[i]*log(phi_i_over_z_i);
-        }
-        return summer;
-    }
-
-    template<typename TType>
-    auto get_Aij(std::size_t i, std::size_t j, const TType& T) const{
-        return forceeval(m(i,j)*T + n(i,j));
-    }
-
-    template<typename TType, typename MoleFractions>
-    auto total(const TType& T, const MoleFractions& molefracs) const {
-
-        using TYPE = std::common_type_t<TType, decltype(molefracs[0])>;
-        TYPE summer = 0.0;
-        for (auto i = 0U; i < molefracs.size(); ++i){
-            auto v_i = b[i];
-            TYPE summerj = 0.0;
-            for (auto j = 0U; j < molefracs.size(); ++j){
-                auto v_j = b[j];
-                auto Aij = get_Aij(i,j,T);
-                auto Omega_ji = v_j/v_i*exp(-Aij/T);
-                summerj += molefracs[j]*Omega_ji;
-            }
-            summer += molefracs[i]*log(summerj);
-        }
-        return forceeval(-summer);
-    }
-
-    // Returns ares/RT
-    template<typename TType, typename MoleFractions>
-    auto operator () (const TType& T, const MoleFractions& molefracs) const {
-        return forceeval(total(T, molefracs) - combinatorial(T, molefracs));
-    }
-};
-
-using ResidualHelmholtzOverRTVariant = std::variant<NullResidualHelmholtzOverRT<double>, WilsonResidualHelmholtzOverRT<double>>;
-
 /**
  Cubic EOS with advanced mixing rules, the EoS/aE method of Jaubert and co-workers
  

include/teqp/models/multifluid/multifluid_activity.hpp

@@ -0,0 +1,72 @@
+#pragma once
+
+#include "teqp/types.hpp"
+
+#include "teqp/models/multifluid.hpp"
+#include "teqp/models/activity/activity_models.hpp"
+
+namespace teqp::multifluid::multifluid_activity {
+
+/**
+ Implementing the general approach of:
+ Jaeger et al.,
+ A theoretically based departure function for multi-fluid mixture models,
+ https://doi.org/10.1016/j.fluid.2018.04.015
+*/
+class MultifluidPlusActivity{
+private:
+    using multifluid_t = decltype(multifluidfactory(nlohmann::json{}));
+    const multifluid_t m_multifluid;
+    const ResidualHelmholtzOverRTVariant m_activity;
+    const std::vector<double> b;
+    const double u;
+public:
+    MultifluidPlusActivity(const nlohmann::json &spec) :
+        m_multifluid(multifluidfactory(spec.at("multifluid"))),
+        m_activity(ares_model_factory(spec.at("activity").at("aresmodel"))),
+        b(spec.at("activity").at("options").at("b")),
+        u(spec.at("activity").at("options").at("u")){}
+
+//    const auto& get_activity() const { return m_activity; }
+
+    template<class VecType>
+    auto R(const VecType& molefrac) const {
+        return get_R_gas<decltype(molefrac[0])>();
+    }
+
+    template <typename TType, typename RhoType, typename MoleFractions>
+    auto alphar_activity(const TType& T, const RhoType& rho, const MoleFractions& molefrac) const {
+        auto gER_over_RT = std::visit([T, &molefrac](const auto& mod){return mod(T, molefrac); }, m_activity); // dimensionless
+        if (static_cast<long>(b.size()) != molefrac.size()){
+            throw teqp::InvalidArgument("Size of mole fractions is incorrect");
+        }
+
+        auto bm = contiguous_dotproduct(b, molefrac);
+
+        const auto& Tcvec = m_multifluid.redfunc.Tc;
+        const auto& vcvec = m_multifluid.redfunc.vc;
+
+        auto rhor = m_multifluid.redfunc.get_rhor(molefrac);
+        auto Tr = m_multifluid.redfunc.get_Tr(molefrac);
+        auto tau = forceeval(Tr/T);
+        auto delta_ref = forceeval(1.0/(u*bm*rhor));
+
+        std::decay_t<std::common_type_t<TType, decltype(molefrac[0])>> summer = 0.0;
+        for (auto i = 0; i < molefrac.size(); ++i){
+             auto delta_i_ref = forceeval(1.0/(u*b[i]/vcvec(i)));
+             auto tau_i = forceeval(Tcvec(i)/T);
+             summer += molefrac(i)*(m_multifluid.alphar_taudeltai(tau, delta_ref, i) - m_multifluid.alphar_taudeltai(tau_i, delta_i_ref, i));
+        }
+        return forceeval(log(1.0+rho*bm)/log(1.0+1.0/u)*(gER_over_RT - summer));
+    }
+
+    template <typename TType, typename RhoType, typename MoleFractions>
+    auto alphar(const TType& T, const RhoType& rho, const MoleFractions& molefrac) const {
+        return forceeval(
+            m_multifluid.alphar(T, rho, molefrac)
+            + alphar_activity(T, rho, molefrac)
+        );
+    }
+};
+
+}; // namespace teqp

interface/CPP/model_factory_multifluid_activity.cpp

@@ -0,0 +1,11 @@
+#include "teqp/cpp/teqpcpp.hpp"
+#include "teqp/cpp/deriv_adapter.hpp"
+#include "teqp/models/multifluid/multifluid_activity.hpp"
+
+namespace teqp{
+    namespace cppinterface{
+        std::unique_ptr<teqp::cppinterface::AbstractModel> make_multifluid_activity(const nlohmann::json &j){
+            return teqp::cppinterface::adapter::make_owned(multifluid::multifluid_activity::MultifluidPlusActivity(j));
+        }
+    }
+}

interface/CPP/teqp_impl_factory.cpp

@@ -24,6 +24,7 @@ namespace teqp {
         std::unique_ptr<teqp::cppinterface::AbstractModel> make_LKP(const nlohmann::json &);
         std::unique_ptr<teqp::cppinterface::AbstractModel> make_multifluid(const nlohmann::json &);
         std::unique_ptr<teqp::cppinterface::AbstractModel> make_multifluid_association(const nlohmann::json &);
+        std::unique_ptr<teqp::cppinterface::AbstractModel> make_multifluid_activity(const nlohmann::json &);
         std::unique_ptr<teqp::cppinterface::AbstractModel> make_multifluid_ECS_HuberEly1994(const nlohmann::json &);
         std::unique_ptr<teqp::cppinterface::AbstractModel> make_AmmoniaWaterTillnerRoth();
         std::unique_ptr<teqp::cppinterface::AbstractModel> make_LJ126_TholJPCRD2016();
@@ -78,6 +79,7 @@ namespace teqp {
             {"multifluid", [](const nlohmann::json& spec){ return make_multifluid(spec);}},
             {"multifluid-ECS-HuberEly1994", [](const nlohmann::json& spec){ return make_multifluid_ECS_HuberEly1994(spec);}},
             {"multifluid-association", [](const nlohmann::json& spec){ return make_multifluid_association(spec);}},
+            {"multifluid-activity", [](const nlohmann::json& spec){ return make_multifluid_activity(spec);}},
             {"AmmoniaWaterTillnerRoth", [](const nlohmann::json& /*spec*/){ return make_AmmoniaWaterTillnerRoth();}},
 
             {"LJ126_TholJPCRD2016", [](const nlohmann::json& /*spec*/){ return make_LJ126_TholJPCRD2016();}},