From c1ddb8d8baba6e0dfefcceedd73320ab58a62461 Mon Sep 17 00:00:00 2001
From: Ian Bell <ian.bell@nist.gov>
Date: Thu, 8 Apr 2021 18:40:58 -0400
Subject: [PATCH] Add MultiFluid to python interface

Also add the accessors as methods of the class in the python interface (overload resolution has large runtime overhead)
---
 include/teqp/models/multifluid.hpp |  7 +++---
 interface/pybind11_wrapper.cpp     | 40 +++++++++++++++++++++---------
 2 files changed, 31 insertions(+), 16 deletions(-)

diff --git a/include/teqp/models/multifluid.hpp b/include/teqp/models/multifluid.hpp
index 51dc498..368f59e 100644
--- a/include/teqp/models/multifluid.hpp
+++ b/include/teqp/models/multifluid.hpp
@@ -33,7 +33,7 @@ public:
 
     template<typename TauType, typename DeltaType, typename MoleFractions>
     auto alphar(const TauType& tau, const DeltaType& delta, const MoleFractions& molefracs) const {
-        using resulttype = typename std::remove_const<decltype(forceeval(tau* delta* molefracs[0]))>::type; // Type promotion, without the const-ness
+        using resulttype = std::common_type_t<decltype(tau), decltype(molefracs[0]), decltype(delta)>; // Type promotion, without the const-ness
         resulttype alphar = 0.0;
         auto N = molefracs.size();
         for (auto i = 0; i < N; ++i) {
@@ -54,7 +54,7 @@ public:
 
     template<typename TauType, typename DeltaType, typename MoleFractions>
     auto alphar(const TauType& tau, const DeltaType& delta, const MoleFractions& molefracs) const {
-        using resulttype = typename std::remove_const<decltype(forceeval(tau* delta* molefracs[0]))>::type; // Type promotion, without the const-ness
+        using resulttype = std::common_type_t<decltype(tau), decltype(molefracs[0]), decltype(delta)>; // Type promotion, without the const-ness
         resulttype alphar = 0.0;
         auto N = molefracs.size();
         for (auto i = 0; i < N; ++i) {
@@ -106,7 +106,6 @@ public:
 class MultiFluidReducingFunction {
 private:
     Eigen::MatrixXd betaT, gammaT, betaV, gammaV, YT, Yv;
-    
 
     template <typename Num>
     auto cube(Num x) const {
@@ -287,7 +286,7 @@ public:
 auto get_departure_function_matrix(const std::string& coolprop_root, const nlohmann::json& BIPcollection, const std::vector<std::string>& components) {
 
     // Allocate the matrix with default models
-    std::vector<std::vector<MultiFluidDepartureFunction>> funcs(2); for (auto i = 0; i < funcs.size(); ++i) { funcs[i].resize(funcs.size()); }
+    std::vector<std::vector<MultiFluidDepartureFunction>> funcs(components.size()); for (auto i = 0; i < funcs.size(); ++i) { funcs[i].resize(funcs.size()); }
 
     auto depcollection = nlohmann::json::parse(std::ifstream(coolprop_root + "/dev/mixtures/mixture_departure_functions.json"));
 
diff --git a/interface/pybind11_wrapper.cpp b/interface/pybind11_wrapper.cpp
index c52b413..feac7e2 100644
--- a/interface/pybind11_wrapper.cpp
+++ b/interface/pybind11_wrapper.cpp
@@ -6,6 +6,7 @@
 
 #include "teqp/core.hpp"
 #include "teqp/models/pcsaft.hpp"
+#include "teqp/models/multifluid.hpp"
 
 namespace py = pybind11;
 
@@ -32,8 +33,8 @@ void add_virials(py::module& m) {
     m.def("get_B12vir", &vd::get_B12vir, py::arg("model"), py::arg("T"), py::arg("molefrac"));
 }
 
-template<typename Model>
-void add_derivatives(py::module &m) {
+template<typename Model, typename Wrapper>
+void add_derivatives(py::module &m, Wrapper &cls) {
     using id = IsochoricDerivatives<Model, double, Eigen::Array<double,Eigen::Dynamic,1> >;
     m.def("get_Ar00", &id::get_Ar00, py::arg("model"), py::arg("T"), py::arg("rho")); 
     m.def("get_Ar10", &id::get_Ar10, py::arg("model"), py::arg("T"), py::arg("rho"));
@@ -47,21 +48,30 @@ void add_derivatives(py::module &m) {
 
     add_virials<Model>(m);
     add_TDx_derivatives<Model>(m);
-}
 
+    cls.def("get_Ar01", [](const Model& m, const double T, const Eigen::ArrayXd& rhovec) { return id::get_Ar01(m, T, rhovec); });
+    cls.def("get_Ar10", [](const Model& m, const double T, const Eigen::ArrayXd& rhovec) { return id::get_Ar10(m, T, rhovec); });
+    using tdx = TDXDerivatives<Model, double, Eigen::Array<double, Eigen::Dynamic, 1> >;
+    cls.def("get_Ar01n", [](const Model& m, const double T, const double rho, const Eigen::ArrayXd& molefrac) { return tdx::get_Ar0n<1>(m, T, rho, molefrac); }); 
+    cls.def("get_Ar02n", [](const Model& m, const double T, const double rho, const Eigen::ArrayXd& molefrac) { return tdx::get_Ar0n<2>(m, T, rho, molefrac); });
+    cls.def("get_Ar03n", [](const Model& m, const double T, const double rho, const Eigen::ArrayXd& molefrac) { return tdx::get_Ar0n<3>(m, T, rho, molefrac); });
+    cls.def("get_Ar04n", [](const Model& m, const double T, const double rho, const Eigen::ArrayXd& molefrac) { return tdx::get_Ar0n<4>(m, T, rho, molefrac); });
+    cls.def("get_Ar05n", [](const Model& m, const double T, const double rho, const Eigen::ArrayXd& molefrac) { return tdx::get_Ar0n<5>(m, T, rho, molefrac); });
+    cls.def("get_Ar06n", [](const Model& m, const double T, const double rho, const Eigen::ArrayXd& molefrac) { return tdx::get_Ar0n<6>(m, T, rho, molefrac); });
+}
 
 void init_teqp(py::module& m) {
 
     using vdWEOSd = vdWEOS<double>;
-    py::class_<vdWEOSd>(m, "vdWEOS")
+    auto wvdW = py::class_<vdWEOSd>(m, "vdWEOS")
         .def(py::init<const std::valarray<double>&, const std::valarray<double>&>(),py::arg("Tcrit"), py::arg("pcrit"))
         ;
-    add_derivatives<vdWEOSd>(m);
-
-    py::class_<vdWEOS1>(m, "vdWEOS1")
-        .def(py::init<const double&, const double&>(), py::arg("a"), py::arg("b"))
+    add_derivatives<vdWEOSd>(m, wvdW);
+    
+    auto wvdW1 = py::class_<vdWEOS1>(m, "vdWEOS1")
+        .def(py::init<const double&, const double&>(), py::arg("a"), py::arg("b"))    
         ;
-    add_derivatives<vdWEOS1>(m);
+    add_derivatives<vdWEOS1>(m, wvdW1);
 
     py::class_<SAFTCoeffs>(m, "SAFTCoeffs")
         .def(py::init<>())
@@ -72,7 +82,7 @@ void init_teqp(py::module& m) {
         .def_readwrite("BibTeXKey", &SAFTCoeffs::BibTeXKey)
         ;
 
-    py::class_<PCSAFTMixture>(m, "PCSAFTEOS")
+    auto wPCSAFT = py::class_<PCSAFTMixture>(m, "PCSAFTEOS")
         .def(py::init<const std::vector<std::string> &>(), py::arg("names"))
         .def(py::init<const std::vector<SAFTCoeffs> &>(), py::arg("coeffs"))
         .def("print_info", &PCSAFTMixture::print_info)
@@ -80,9 +90,15 @@ void init_teqp(py::module& m) {
         .def("get_m", &PCSAFTMixture::get_m)
         .def("get_sigma_Angstrom", &PCSAFTMixture::get_sigma_Angstrom)
         .def("get_epsilon_over_k_K", &PCSAFTMixture::get_epsilon_over_k_K)
-
         ;
-    add_derivatives<PCSAFTMixture>(m);
+    add_derivatives<PCSAFTMixture>(m, wPCSAFT);
+
+    // Multifluid model
+    m.def("build_multifluid_model", &build_multifluid_model);
+    using MultiFluid = decltype(build_multifluid_model(std::vector<std::string>{"",""},"",""));
+    using idMF = IsochoricDerivatives<MultiFluid, double, Eigen::Array<double, Eigen::Dynamic, 1> >;
+    auto wMF = py::class_<MultiFluid>(m, "MultiFluid");
+    add_derivatives<MultiFluid>(m, wMF);
 
     // for timing testing
     m.def("mysummer", [](const double &c, const Eigen::ArrayXd &x) { return c*x.sum(); });
-- 
GitLab