From c3508c3f55cfeb764b87d63de4fe1064bf70615b Mon Sep 17 00:00:00 2001
From: Ian Bell <ian.h.bell@gmail.com>
Date: Fri, 15 Aug 2014 15:42:47 +0200
Subject: [PATCH] Added parallel evaluation of residual helmholtz terms

Should in the end only slow down by 1.2x, but calculate at least up to second order derivatives all at once, so net speedup of 5 if calculating all first and second order partials !!

Signed-off-by: Ian Bell <ian.h.bell@gmail.com>
---
 include/Helmholtz.h |  9 ++++++++
 src/Helmholtz.cpp   | 54 +++++++++++++++++++++++++++++++++++++++++++++
 2 files changed, 63 insertions(+)

diff --git a/include/Helmholtz.h b/include/Helmholtz.h
index 96f3d11b..b869ecd8 100644
--- a/include/Helmholtz.h
+++ b/include/Helmholtz.h
@@ -4,6 +4,7 @@
 
 #include <vector>
 #include "rapidjson/rapidjson_include.h"
+#include "time.h"
 
 namespace CoolProp{
 
@@ -97,6 +98,11 @@ public:
 // #############################################################################
 // #############################################################################
 
+struct Derivatives
+{
+    long double alphar, dalphar_ddelta, dalphar_dtau, d2alphar_ddelta2, d2alphar_dtau2, d2alphar_ddelta_dtau;
+};
+
 struct ResidualHelmholtzPowerElement
 {
     long double n,d,t,ld;
@@ -126,6 +132,7 @@ public:
     {
         N = n.size();
         s.resize(N);
+        
         for (std::size_t i = 0; i < n.size(); ++i)
         {
             ResidualHelmholtzPowerElement el;
@@ -153,6 +160,8 @@ public:
     long double dDelta2_dTau(const long double &tau, const long double &delta) throw();
     long double dDelta_dTau2(const long double &tau, const long double &delta) throw();
     long double dTau3(const long double &tau, const long double &delta) throw();
+    
+    void all(const long double &tau, const long double &delta, Derivatives &derivs) throw();
 };
 
 struct ResidualHelmholtzExponentialElement
diff --git a/src/Helmholtz.cpp b/src/Helmholtz.cpp
index 9ed9ea24..5a4642af 100644
--- a/src/Helmholtz.cpp
+++ b/src/Helmholtz.cpp
@@ -34,6 +34,60 @@ long double kahanSum(std::vector<long double> &x)
     return sum;
 }
 bool wayToSort(long double i, long double j) { return std::abs(i) > std::abs(j); }
+
+void ResidualHelmholtzPower::all(const long double &tau, const long double &delta, Derivatives &derivs) throw()
+{
+    long double log_tau = log(tau), log_delta = log(delta), u, du_ddelta, du_dtau, d2u_ddelta2, d2u_dtau2, ndteu, 
+                one_over_delta = 1/delta, one_over_tau = 1/tau, // division is much slower than multiplication, so do one division here
+                B_delta, B_tau, B_delta2, B_tau2;
+    derivs.alphar = 0.0;
+    derivs.dalphar_ddelta = 0.0;
+    derivs.dalphar_dtau = 0.0;
+    derivs.d2alphar_ddelta2 = 0.0;
+    derivs.d2alphar_dtau2 = 0.0;
+    derivs.d2alphar_ddelta_dtau = 0.0;
+    
+    for (std::size_t i = 0; i < N; ++i)
+    {
+        ResidualHelmholtzPowerElement &el = elements[i];
+        long double ni = el.n, di = el.d, ti = el.t;
+        int li = el.l;
+        if (li > 0){
+            u = -pow(delta, li);
+            du_ddelta = li*u*one_over_delta;
+            du_dtau = 0;
+            d2u_ddelta2 = (li-1)*du_ddelta*one_over_delta;
+            d2u_dtau2 = 0;
+        }
+        else{
+            u = 0;
+            du_ddelta = 0;
+            du_dtau = 0;
+            d2u_ddelta2 = 0;
+            d2u_dtau2 = 0;
+        }
+        ndteu = ni*exp(ti*log_tau+di*log_delta+u);
+        
+        B_delta = (delta*du_ddelta + di);
+        B_tau = (tau*du_dtau + ti);
+        B_delta2 = (POW2(delta)*(d2u_ddelta2 + POW2(du_ddelta)) + 2*di*delta*du_ddelta + di*(di-1));
+        B_tau2 = (POW2(tau)*(d2u_dtau2 + POW2(du_dtau)) + 2*ti*tau*du_dtau + ti*(ti-1));
+        
+        derivs.alphar += ndteu;
+        derivs.dalphar_ddelta += ndteu*B_delta;
+        derivs.dalphar_dtau += ndteu*B_tau;
+        derivs.d2alphar_ddelta2 += ndteu*B_delta2;
+        derivs.d2alphar_dtau2 += ndteu*B_tau2;
+        derivs.d2alphar_ddelta_dtau += ndteu*B_delta*B_tau;
+    }
+    derivs.dalphar_ddelta       *= one_over_delta;
+    derivs.dalphar_dtau         *= one_over_tau;
+    derivs.d2alphar_ddelta2     *= POW2(one_over_delta);
+    derivs.d2alphar_dtau2       *= POW2(one_over_tau);
+    derivs.d2alphar_ddelta_dtau *= one_over_delta*one_over_tau;
+    return;
+};
+
 long double ResidualHelmholtzPower::base(const long double &tau, const long double &delta) throw()
 {
     if (N==0){return 0.0;}