Lots of small fixes to snippets

Closes #54

CMakeLists.txt

@@ -231,6 +231,7 @@ if (TEQP_SNIPPETS)
       target_link_libraries (${snippet_exe} PRIVATE ${CMAKE_DL_LIBS})
     endif()
     target_compile_definitions(${snippet_exe} PRIVATE -DTEQP_MULTICOMPLEX_ENABLED)
+    target_compile_definitions(${snippet_exe} PRIVATE -DTEQP_MULTIPRECISION_ENABLED)
     target_compile_definitions(${snippet_exe} PRIVATE -DUSE_AUTODIFF)
 
     if(TEQP_JAVASCRIPT_HTML)

include/teqp/models/multifluid.hpp

@@ -99,7 +99,7 @@ class DepartureContribution {
     /// Call a single departure term at i,j 
     template<typename TauType, typename DeltaType>
     auto get_alpharij(const int i, const int j,     const TauType& tau, const DeltaType& delta) const {
-        int N = funcs.size();
+        std::size_t N = funcs.size();
         if (i < 0 || j < 0){
             throw teqp::InvalidArgument("i or j is negative");
         }

src/bench_AbstractModel.cpp

@@ -5,6 +5,7 @@
 #include "teqp/cpp/derivs.hpp"
 #include "teqp/derivs.hpp"
 #include "teqp/cpp/deriv_adapter.hpp"
+#include "teqp/models/multifluid.hpp"
 
 using namespace teqp;
 
@@ -63,8 +64,9 @@ TEST_CASE("multifluid derivatives via DerivativeAdapter", "[mf]")
             {"departure", "../mycp/dev/mixtures/mixture_departure_functions.json"}
         }
     }};
-    auto am = teqp::cppinterface::make_model(j);
-    auto model = multifluidfactory(j["model"]);
+    const auto am = teqp::cppinterface::make_model(j);
+    using multifluid_t = decltype(multifluidfactory(nlohmann::json{}));
+    auto model = teqp::cppinterface::adapter::get_model_cref<multifluid_t>(am.get());
 
     auto z = (Eigen::ArrayXd(1) << 1.0).finished();
 

src/bench_derivbox.cpp

@@ -24,14 +24,15 @@ TEST_CASE("multifluid derivatives", "[mf]")
     nlohmann::json jigs = nlohmann::json::array(); jigs.push_back(jig);
     auto aig = teqp::IdealHelmholtz(jigs);
 
-    auto amm = teqp::cppinterface::make_multifluid_model(names, "../mycp");
-    auto aigg = teqp::cppinterface::make_model({{"kind","IdealHelmholtz"}, {"model", jigs}});
+    using namespace teqp::cppinterface;
+    std::shared_ptr<AbstractModel> amm = teqp::cppinterface::make_multifluid_model(names, "../mycp");
+    std::shared_ptr<AbstractModel> aigg = teqp::cppinterface::make_model({{"kind","IdealHelmholtz"}, {"model", jigs}});
 
     std::vector<char> vars = {'P', 'S'};
     const auto vals = (Eigen::ArrayXd(2) << 300.0, 400.0).finished();
 
     Eigen::Ref<const Eigen::ArrayXd> rvals = vals, rz = z;
-    teqp::iteration::NRIterator NR(amm, aigg, vars, rvals, T, rho, rz);
+    teqp::iteration::NRIterator NR(amm.get(), aigg.get(), vars, rvals, T, rho, rz);
 
     BENCHMARK("alphar") {
         return model.alphar(T, rho, z);
@@ -53,7 +54,7 @@ TEST_CASE("multifluid derivatives", "[mf]")
     };
 
     BENCHMARK("Newton-Raphson construction") {
-        return teqp::iteration::NRIterator(amm, aigg, vars, rvals, T, rho, rz);
+        return teqp::iteration::NRIterator(amm.get(), aigg.get(), vars, rvals, T, rho, rz);
     };
     BENCHMARK("Newton-Raphson calc_step") {
         return NR.calc_step(T, rho);

src/multifluid_crit.cpp

@@ -44,7 +44,7 @@ void trace_critical_loci(const std::string &coolprop_root, const nlohmann::json
                double zi = rhovec[i] / rhovec.sum();
                double T = zi * model.redfunc.Tc[i] + (1 - zi) * model.redfunc.Tc[1L - i];
                double z0 = (i == 0) ? zi : 1-zi;
-               auto [Tnew, rhonew] = ct::critical_polish_molefrac(model, T, rhovec, z0);
+               auto [Tnew, rhonew] = ct::critical_polish_fixedmolefrac(model, T, rhovec, z0);
                T0 = Tnew;
                rhoc0 = rhovec.sum();
            }
@@ -62,4 +62,4 @@ int main(){
     Timer t(1);
     trace_critical_loci(coolprop_root, BIPcollection);
     return EXIT_SUCCESS;
-}
+}

src/sat_Z_accuracy.cpp

@@ -161,10 +161,15 @@ auto with_teqp_and_boost(const Model &model, double T, double rho, const VECTOR
 
     // Now do the third-order derivative of alphar, as a further test
     // Define a generic lambda function taking rho
-    auto ff = [&](const auto& rho){ return model.alphar(T, rho, z); };
+    auto ff = [&](const my_float& rho) -> my_float {
+        auto tau = forceeval(Tc/T);
+        auto delta = forceeval(rho/rhoc);
+        return alphar_Lemmon2009<my_float>(tau, delta);
+    };
     my_float drho = 1e-30*rho;
 
-    o.Ar02exact = static_cast<double>(centered_diff<2,6>(ff,static_cast<my_float>(rho),drho)*pow(rho, 2));
+    auto Ar02_extended = centered_diff<2,6>(ff,static_cast<my_float>(rho),drho)*pow(rho, 2);
+    o.Ar02exact = static_cast<double>(Ar02_extended);
     o.Ar02teqp = tdx::template get_Ar0n<2>(model, T, rho, z)[2];
 
     o.Ar03exact = static_cast<double>(centered_diff<3,6>(ff,static_cast<my_float>(rho),drho)*pow(rho, 3));
@@ -251,4 +256,4 @@ int main()
     do_one("PROPANE", "n-Propane"); 
     do_one("WATER", "Water");
     return EXIT_SUCCESS;
-}
+}