From c454cf851c921a1406d7e3f305a7fba51999102c Mon Sep 17 00:00:00 2001
From: Ian Bell <ian.bell@nist.gov>
Date: Wed, 2 Nov 2022 20:26:18 -0400
Subject: [PATCH] Add all the Ar methods that were already exposed

Closing in
---
 include/teqp/derivs.hpp     |  4 ++--
 interface/CPP/teqpcpp.cpp   | 23 +++++++++++++++++++++
 interface/CPP/teqpcpp.hpp   | 41 ++++++++++++++++++++++++++++++++++++-
 src/bench_AbstractModel.cpp |  8 +++++++-
 4 files changed, 72 insertions(+), 4 deletions(-)

diff --git a/include/teqp/derivs.hpp b/include/teqp/derivs.hpp
index 0b23860..a19ba50 100644
--- a/include/teqp/derivs.hpp
+++ b/include/teqp/derivs.hpp
@@ -211,7 +211,7 @@ struct TDXDerivatives {
     *
     * Note: none of the intermediate derivatives are returned, although they are calculated
     */
-    template<int iT, int iD, ADBackends be>
+    template<int iT, int iD, ADBackends be = ADBackends::autodiff>
     static auto get_Arxy(const Model& model, const Scalar& T, const Scalar& rho, const VectorType& molefrac) {
         auto wrapper = AlphaCallWrapper<AlphaWrapperOption::residual, decltype(model)>(model);
         if constexpr (iT == 0 && iD == 0) {
@@ -230,7 +230,7 @@ struct TDXDerivatives {
     *
     * Note: none of the intermediate derivatives are returned, although they are calculated
     */
-    template<int iT, int iD, ADBackends be>
+    template<int iT, int iD, ADBackends be = ADBackends::autodiff>
     static auto get_Aigxy(const Model& model, const Scalar& T, const Scalar& rho, const VectorType& molefrac) {
         auto wrapper = AlphaCallWrapper<AlphaWrapperOption::idealgas, decltype(model)>(model);
         if constexpr (iT == 0 && iD == 0) {
diff --git a/interface/CPP/teqpcpp.cpp b/interface/CPP/teqpcpp.cpp
index 2d4ee1c..4d4e535 100644
--- a/interface/CPP/teqpcpp.cpp
+++ b/interface/CPP/teqpcpp.cpp
@@ -91,6 +91,29 @@ namespace teqp {
                 }, m_model);
             }
             
+            // Here XMacros are used to create functions like get_Ar00, get_Ar01, ....
+            #define X(i,j) \
+            double get_Ar ## i ## j(const double T, const double rho, const EArrayd& molefracs) const override { \
+                return std::visit([&](const auto& model) { \
+                    using tdx = teqp::TDXDerivatives<decltype(model), double, EArrayd>; \
+                    return tdx::template get_Arxy<i,j>(model, T, rho, molefracs); \
+                }, m_model); \
+            }
+            ARXY_args
+            #undef X
+            
+            // Here XMacros are used to create functions like get_Ar01n, get_Ar02n, ....
+            #define X(i) \
+            EArrayd get_Ar0 ## i ## n(const double T, const double rho, const EArrayd& molefracs) const override { \
+                return std::visit([&](const auto& model) { \
+                    using tdx = teqp::TDXDerivatives<decltype(model), double, EArrayd>; \
+                    auto vals = tdx::template get_Ar0n<i>(model, T, rho, molefracs); \
+                    return Eigen::Map<Eigen::ArrayXd>(&(vals[0]), vals.size());\
+                }, m_model); \
+            }
+            AR0N_args
+            #undef X
+            
             // Methods only available for PC-SAFT
             EArrayd get_m() const override {
                 return get_or_fail<PCSAFT_t>("PCSAFT").get_m();
diff --git a/interface/CPP/teqpcpp.hpp b/interface/CPP/teqpcpp.hpp
index 874c0bc..fb11ce9 100644
--- a/interface/CPP/teqpcpp.hpp
+++ b/interface/CPP/teqpcpp.hpp
@@ -12,6 +12,33 @@ using EArray2 = Eigen::Array<double, 2, 1>;
 using EArrayd = Eigen::ArrayX<double>;
 using EArray33d = Eigen::Array<double, 3, 3>;
 
+#define ARXY_args \
+    X(0,0) \
+    X(0,1) \
+    X(0,2) \
+    X(0,3) \
+    X(0,4) \
+    X(1,0) \
+    X(1,1) \
+    X(1,2) \
+    X(1,3) \
+    X(1,4) \
+    X(2,0) \
+    X(2,1) \
+    X(2,2) \
+    X(2,3) \
+    X(2,4)
+
+#define AR0N_args \
+    X(0) \
+    X(1) \
+    X(2) \
+    X(3) \
+    X(4) \
+    X(5) \
+    X(6)
+
+
 namespace teqp {
     namespace cppinterface {
 
@@ -29,13 +56,25 @@ namespace teqp {
         */
         class AbstractModel {
         public:
-            virtual double get_Arxy(const int, const int, const double, const double, const EArrayd&) const = 0;
+            
             virtual nlohmann::json trace_critical_arclength_binary(const double T0, const EArrayd& rhovec0, const std::optional<std::string>&, const std::optional<TCABOptions> &) const = 0;
             virtual EArray2 pure_VLE_T(const double T, const double rhoL, const double rhoV, int maxiter) const = 0;
             virtual EArrayd get_fugacity_coefficients(const double T, const EArrayd& rhovec) const = 0;
             virtual EArrayd get_partial_molar_volumes(const double T, const EArrayd& rhovec) const = 0;
             virtual EArray33d get_deriv_mat2(const double T, double rho, const EArrayd& z) const = 0;
             
+            virtual double get_Arxy(const int, const int, const double, const double, const EArrayd&) const = 0;
+            
+            // Here XMacros are used to create functions like get_Ar00, get_Ar01, ....
+            #define X(i,j) virtual double get_Ar ## i ## j(const double T, const double rho, const EArrayd& molefrac) const = 0;
+                ARXY_args
+            #undef X
+            // And like get_Ar01n, get_Ar02n, ....
+            #define X(i) virtual EArrayd get_Ar0 ## i ## n(const double T, const double rho, const EArrayd& molefrac) const = 0;
+                AR0N_args
+            #undef X
+            
+            
             // 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/src/bench_AbstractModel.cpp b/src/bench_AbstractModel.cpp
index 89896ac..71701ac 100644
--- a/src/bench_AbstractModel.cpp
+++ b/src/bench_AbstractModel.cpp
@@ -21,11 +21,17 @@ TEST_CASE("multifluid derivatives", "[mf]")
     auto am2 = teqp::cppinterface::make_vdW1(2, 3);
 
     auto z = (Eigen::ArrayXd(1) << 1.0).finished();
-    auto rhovec = 300.0*z;
+    auto rhovec = 300.0* z;
     
     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);
     };
-- 
GitLab