Switch to new methods for pressure calculation

include/teqp/derivs.hpp

@@ -29,6 +29,29 @@ derivrhoi(const FuncType& f, TType T, const ContainerType& rho, Integer i) {
     return f(T, rhocom).imag() / h;
 }
 
+template <typename TType, typename ContainerType, typename Model>
+typename std::enable_if<is_container<ContainerType>::value, typename ContainerType::value_type>::type
+get_Psir(const Model& model, const TType T, const ContainerType& rhovec) {
+    using container = decltype(rhovec);
+    auto rhotot_ = std::accumulate(std::begin(rhovec), std::end(rhovec), (decltype(rhovec[0]))0.0);
+    return model.alphar(T, rhovec)*model.R*T*rhotot_;
+}
+
+/**
+// Calculate the residual pressure from derivatives of alphar
+*/
+template <typename Model, typename TType, typename ContainerType>
+typename std::enable_if<is_container<ContainerType>::value, typename ContainerType::value_type>::type
+get_pr(const Model& model, const TType T, const ContainerType& rhovec) {
+    using container = decltype(rhovec);
+    auto rhotot_ = std::accumulate(std::begin(rhovec), std::end(rhovec), (decltype(rhovec[0]))0.0);
+    decltype(rhovec[0] * T) pr = 0.0;
+    for (auto i = 0; i < rhovec.size(); ++i) {
+        pr += rhovec[i]*derivrhoi([&model](const auto& T, const auto& rhovec){ return model.alphar(T, rhovec); }, T, rhovec, i);
+    }
+    return pr*rhotot_*model.R*T;
+}
+
 template<typename Model, typename TType, typename RhoType>
 auto build_Psir_Hessian(const Model& model, const TType T, const RhoType& rho) {
     // Double derivatives in each component's concentration

src/main.cpp

@@ -7,6 +7,7 @@
 #include <complex>
 #include <chrono>
 #include <optional>
+#include <iomanip>
 
 #include "teqp/core.hpp"
 
@@ -29,22 +30,12 @@ void test_vdW() {
     std::cout << std::chrono::duration<double>(t3 - t2).count() << " from p(T,v)" << std::endl;
 
     const std::valarray<double> rhovec = { rho, 0.0 };
-
     auto t21 = std::chrono::steady_clock::now();
-    
-    auto fPsir = [&vdW](const auto& T, const auto& rhovec) {
-        using container = decltype(rhovec);
-        auto rhotot_ = std::accumulate(std::begin(rhovec), std::end(rhovec), (decltype(rhovec[0]))0.0);
-        return vdW.alphar(T, rhovec) * vdW.R * T * rhotot_;
-    };
-    auto Psir = fPsir(T, rhovec);
-    auto dPsirdrho0 = rhovec[0] * derivrhoi(fPsir, T, rhovec, 0);
-    auto dPsirdrho1 = rhovec[1] * derivrhoi(fPsir, T, rhovec, 1);
-    auto pfromderiv = rho * R * T - Psir + dPsirdrho0 + dPsirdrho1;
-
+    auto pfromderiv = rho*R*T + get_pr(vdW, T, rhovec);
     auto t31 = std::chrono::steady_clock::now();
     std::cout << std::chrono::duration<double>(t31 - t21).count() << " from isochoric" << std::endl;
-    std::cout << pfromderiv / pp - 1 << std::endl;
+    auto err = pfromderiv / pp - 1.0;
+    std::cout << std::setprecision(20) << "Error (fractional): " << err << std::endl;
 }
 
 void test_vdwMix() {
@@ -53,15 +44,13 @@ void test_vdwMix() {
     std::valarray<double> pc_Pa = { 4863000.0, 5842000.0 };
     vdWEOS<double> vdW(Tc_K, pc_Pa);
 
-    double T = 298.15;
+    volatile double T = 298.15;
     auto rho = 3.0;
     auto R = get_R_gas<double>();
     auto rhotot = rho;
 
     const std::valarray<double> rhovec = { rho/2, rho/2 };
 
-    auto t2 = std::chrono::steady_clock::now();
-
     auto fPsir = [&vdW](const auto& T, const auto& rhovec) {
         using container = decltype(rhovec);
         auto rhotot_ = std::accumulate(std::begin(rhovec), std::end(rhovec), (decltype(rhovec[0]))0.0);
@@ -71,21 +60,15 @@ void test_vdwMix() {
     auto dPsirdrho0 = rhovec[0]*derivrhoi(fPsir, T, rhovec, 0);
     auto dPsirdrho1 = rhovec[1]*derivrhoi(fPsir, T, rhovec, 1);
     auto pfromderiv = rho*R*T - Psir + dPsirdrho0 + dPsirdrho1;
-    {
-        auto term0 = rhovec[0] * derivrhoi([&vdW, rhotot](const auto& T, const auto& rhovec) { return vdW.alphar(T, rhovec); }, T, rhovec, 0);
-        auto term1 = rhovec[1] * derivrhoi([&vdW, rhotot](const auto& T, const auto& rhovec) { return vdW.alphar(T, rhovec); }, T, rhovec, 1);
-        auto pr = (term0 + term1)*rhotot*R*T;
-        auto pfromderiv2 = rho*R*T + pr;
-        auto err2 = pfromderiv / pfromderiv2 - 1;
-        int err = 0;
-    }
 
+    auto t2 = std::chrono::steady_clock::now();
+    auto pfromderiv3 = rhotot*R*T + get_pr(vdW, T, rhovec);
     auto t3 = std::chrono::steady_clock::now();
     std::cout << std::chrono::duration<double>(t3 - t2).count() << " from isochoric (mix) " << std::endl;
 }
 
 int main(){
-    //test_vdW();
+    test_vdW();
     test_vdwMix();
     return EXIT_SUCCESS;
 }
\ No newline at end of file