From f7464fcaa23d2efd46c4399d43c04b2631bda6fe Mon Sep 17 00:00:00 2001 From: Ian Bell Date: Sun, 17 Aug 2014 21:09:27 +0200 Subject: [PATCH] Implemented Eigen solution for all() function, speed testing next Signed-off-by: Ian Bell --- dev/codelite/coolprop.project | 2 +- include/Helmholtz.h | 57 ++++++++-- src/Backends/Helmholtz/Fluids/FluidLibrary.h | 3 + src/Helmholtz.cpp | 111 ++++++++++++++++++- src/main.cxx | 30 +++++ 5 files changed, 193 insertions(+), 10 deletions(-) diff --git a/dev/codelite/coolprop.project b/dev/codelite/coolprop.project index 4fd7fcbd..874dc2b0 100644 --- a/dev/codelite/coolprop.project +++ b/dev/codelite/coolprop.project @@ -190,7 +190,7 @@ - + diff --git a/include/Helmholtz.h b/include/Helmholtz.h index 55039bee..026475d3 100644 --- a/include/Helmholtz.h +++ b/include/Helmholtz.h @@ -4,6 +4,7 @@ #include #include "rapidjson/rapidjson_include.h" +#include "Eigen/Core" #include "time.h" namespace CoolProp{ @@ -101,7 +102,7 @@ public: struct HelmholtzDerivatives { - long double alphar, dalphar_ddelta, dalphar_dtau, d2alphar_ddelta2, d2alphar_dtau2, d2alphar_ddelta_dtau, + double alphar, dalphar_ddelta, dalphar_dtau, d2alphar_ddelta2, d2alphar_dtau2, d2alphar_ddelta_dtau, d3alphar_ddelta3, d3alphar_ddelta_dtau2, d3alphar_ddelta2_dtau, d3alphar_dtau3; void reset(){alphar = 0; dalphar_ddelta = 0; dalphar_dtau = 0; d2alphar_ddelta2 = 0; d2alphar_dtau2 = 0; d2alphar_ddelta_dtau = 0; d3alphar_ddelta3 = 0; d3alphar_ddelta_dtau2 = 0; d3alphar_ddelta2_dtau = 0; d3alphar_dtau3 = 0; @@ -121,20 +122,28 @@ struct ResidualHelmholtzGeneralizedExponentialElement ResidualHelmholtzGeneralizedExponentialElement() { - n = _HUGE; d = _HUGE; t = _HUGE; - c = _HUGE; l_double = _HUGE; omega = _HUGE; m_double = _HUGE; - eta1 = _HUGE; epsilon1 = _HUGE; eta2 = _HUGE; epsilon2 = _HUGE; - beta1 = _HUGE; gamma1 = _HUGE; beta2 = _HUGE; gamma2 = _HUGE; - l_int = -100000; - m_int = -100000; + n = 0; d = 0; t = 0; c = 0; + l_double = 0; omega = 0; m_double = 0; + eta1 = 0; epsilon1 = 0; eta2 = 0; epsilon2 = 0; + beta1 = 0; gamma1 = 0; beta2 = 0; gamma2 = 0; + l_int = 0; m_int = 0; } }; class ResidualHelmholtzGeneralizedExponential : public BaseHelmholtzTerm{ public: - bool delta_li_in_u, tau_mi_in_u, eta1_in_u, eta2_in_u, beta1_in_u, beta2_in_u; + bool delta_li_in_u, tau_mi_in_u, eta1_in_u, eta2_in_u, beta1_in_u, beta2_in_u, finished; std::vector s; std::size_t N; + + /// These variables are for the exp(u) part + /// u is given by -c*delta^l_i-omega*tau^m_i-eta1*(delta-epsilon1)-eta2*(delta-epsilon2)^2-beta1*(tau-gamma1)-beta2*(tau-gamma2)^2 + std::vector n,d,t,c, l_double, omega, m_double, eta1, epsilon1, eta2, epsilon2, beta1, gamma1, beta2, gamma2; + /// If l_i or m_i are integers, we will store them as integers in order to call pow(double, int) rather than pow(double, double) + std::vector l_int, m_int; + + Eigen::ArrayXd uE, du_ddeltaE, du_dtauE, d2u_ddelta2E, d2u_dtau2E, d3u_ddelta3E, d3u_dtau3E; + std::vector elements; // Default Constructor ResidualHelmholtzGeneralizedExponential(){N = 0; @@ -144,6 +153,7 @@ public: eta2_in_u = false; beta1_in_u = false; beta2_in_u = false; + finished = false; }; void add_Power(const std::vector &n, const std::vector &d, @@ -252,6 +262,36 @@ public: delta_li_in_u = true; tau_mi_in_u = true; }; + + void finish(){ + n.resize(elements.size()); d.resize(elements.size()); + t.resize(elements.size()); c.resize(elements.size()); + l_double.resize(elements.size()); l_int.resize(elements.size()); + epsilon2.resize(elements.size()); eta2.resize(elements.size()); + gamma2.resize(elements.size()); beta2.resize(elements.size()); + + for (std::size_t i = 0; i < elements.size(); ++i){ + n[i] = elements[i].n; + d[i] = elements[i].d; + t[i] = elements[i].t; + c[i] = elements[i].c; + l_double[i] = elements[i].l_double; + l_int[i] = elements[i].l_int; + epsilon2[i] = elements[i].epsilon2; + eta2[i] = elements[i].eta2; + gamma2[i] = elements[i].gamma2; + beta2[i] = elements[i].beta2; + } + uE.resize(elements.size()); + du_ddeltaE.resize(elements.size()); + du_dtauE.resize(elements.size()); + d2u_ddelta2E.resize(elements.size()); + d2u_dtau2E.resize(elements.size()); + d3u_ddelta3E.resize(elements.size()); + d3u_dtau3E.resize(elements.size()); + + finished = true; + } ///< Destructor for the class. No implementation ~ResidualHelmholtzGeneralizedExponential(){}; @@ -270,6 +310,7 @@ public: long double dTau3(const long double &tau, const long double &delta) throw(){return 0;}; void all(const long double &tau, const long double &delta, HelmholtzDerivatives &derivs) throw(); + void allEigen(const long double &tau, const long double &delta, HelmholtzDerivatives &derivs) throw(); }; struct ResidualHelmholtzNonAnalyticElement diff --git a/src/Backends/Helmholtz/Fluids/FluidLibrary.h b/src/Backends/Helmholtz/Fluids/FluidLibrary.h index 66fe69a2..74ebd984 100644 --- a/src/Backends/Helmholtz/Fluids/FluidLibrary.h +++ b/src/Backends/Helmholtz/Fluids/FluidLibrary.h @@ -129,6 +129,9 @@ protected: throw ValueError(format("Unsupported Residual helmholtz type: %s",type.c_str())); } } + + // Finish adding parts to the Generalized Exponential term, build other vectors + EOS.alphar.GenExp.finish(); }; /// Parse the contributions to the ideal-gas Helmholtz energy diff --git a/src/Helmholtz.cpp b/src/Helmholtz.cpp index 1745db10..d5484939 100644 --- a/src/Helmholtz.cpp +++ b/src/Helmholtz.cpp @@ -1,6 +1,7 @@ #include #include "Helmholtz.h" + namespace CoolProp{ long double kahanSum(std::vector &x) @@ -18,13 +19,119 @@ long double kahanSum(std::vector &x) } bool wayToSort(long double i, long double j) { return std::abs(i) > std::abs(j); } +// define function to be applied coefficient-wise +double ramp(double x) +{ + if (x > 0) + return x; + else + return 0; +} + +void ResidualHelmholtzGeneralizedExponential::allEigen(const long double &tau, const long double &delta, HelmholtzDerivatives &derivs) throw() +{ + double log_tau = log(tau), log_delta = log(delta), + one_over_delta = 1/delta, one_over_tau = 1/tau; // division is much slower than multiplication, so do one division here + + Eigen::Map nE(n.data(), elements.size()); + Eigen::Map dE(d.data(), elements.size()); + Eigen::Map tE(t.data(), elements.size()); + Eigen::Map cE(c.data(), elements.size()); + Eigen::Map l_intE(l_int.data(), elements.size()); + Eigen::Map l_doubleE(l_double.data(), elements.size()); + Eigen::Map eta1E(eta1.data(), elements.size()); + Eigen::Map eta2E(eta2.data(), elements.size()); + Eigen::Map epsilon1E(epsilon1.data(), elements.size()); + Eigen::Map epsilon2E(epsilon2.data(), elements.size()); + Eigen::Map beta1E(beta1.data(), elements.size()); + Eigen::Map beta2E(beta2.data(), elements.size()); + Eigen::Map gamma1E(gamma1.data(), elements.size()); + Eigen::Map gamma2E(gamma2.data(), elements.size()); + + // **************************************** + // The u part in exp(u) and its derivatives + // **************************************** + + // Set the u part of exp(u) to zero + uE.fill(0); + du_ddeltaE.fill(0); + du_dtauE.fill(0); + d2u_ddelta2E.fill(0); + d2u_dtau2E.fill(0); + d3u_ddelta3E.fill(0); + d3u_dtau3E.fill(0); + + if (delta_li_in_u){ + Eigen::ArrayXd u_increment = -cE*Eigen::exp(log_delta*l_doubleE); //pow(delta,L) -> exp(L*log(delta)) + uE += u_increment; + du_ddeltaE += l_doubleE*u_increment*one_over_delta; + d2u_ddelta2E += (l_doubleE-1)*l_doubleE*u_increment*one_over_delta*one_over_delta; + d3u_ddelta3E += (l_doubleE-2)*(l_doubleE-1)*l_doubleE*u_increment*one_over_delta*one_over_delta*one_over_delta; + } + /* + if (tau_mi_in_u){ + long double omegai = el.omega, m_double = el.m_double; + if (std::abs(m_double) > 0){ + long double u_increment = -omegai*pow(tau, m_double); + long double du_dtau_increment = m_double*u_increment*one_over_tau; + long double d2u_dtau2_increment = (m_double-1)*du_dtau_increment*one_over_tau; + long double d3u_dtau3_increment = (m_double-2)*d2u_dtau2_increment*one_over_tau; + u += u_increment; + du_dtau += du_dtau_increment; + d2u_dtau2 += d2u_dtau2_increment; + d3u_dtau3 += d3u_dtau3_increment; + } + }*/ + if (eta1_in_u){ + uE += -eta1E*(delta-epsilon1E); + du_ddeltaE += -eta1E; + } + if (eta2_in_u){ + uE += -eta2E*POW2(delta-epsilon2E); + du_ddeltaE += -2*eta2E*(delta-epsilon2E); + d2u_ddelta2E += -2*eta2E; + } + if (beta1_in_u){ + uE += -beta1E*(tau-gamma1E); + du_dtauE += -beta1E; + } + if (beta2_in_u){ + uE += -beta2E*POW2(tau-gamma2E); + du_dtauE += -2*beta2E*(tau-gamma2E); + d2u_dtau2E += -2*beta2E; + } + + Eigen::ArrayXd ndteuE = nE*exp(tE*log_tau + dE*log_delta + uE); + Eigen::ArrayXd B_deltaE = delta*du_ddeltaE + dE; + Eigen::ArrayXd B_tauE = tau*du_dtauE + tE; + Eigen::ArrayXd B_delta2E = POW2(delta)*(d2u_ddelta2E + du_ddeltaE.square()) + 2*dE*delta*du_ddeltaE + dE*(dE-1); + Eigen::ArrayXd B_tau2E = POW2(tau)*(d2u_dtau2E + du_dtauE.square()) + 2*tE*tau*du_dtauE + tE*(tE-1); + Eigen::ArrayXd B_delta3E = POW3(delta)*d3u_ddelta3E + 3*dE*POW2(delta)*d2u_ddelta2E+3*POW3(delta)*d2u_ddelta2E*du_ddeltaE+3*dE*POW2(delta*du_ddeltaE)+3*dE*(dE-1)*delta*du_ddeltaE+dE*(dE-1)*(dE-2)+POW3(delta*du_ddeltaE); + Eigen::ArrayXd B_tau3E = POW3(tau)*d3u_dtau3E + 3*tE*POW2(tau)*d2u_dtau2E+3*POW3(tau)*d2u_dtau2E*du_dtauE+3*tE*POW2(tau*du_dtauE)+3*tE*(tE-1)*tau*du_dtauE+tE*(tE-1)*(tE-2)+POW3(tau*du_dtauE); + + derivs.alphar += ndteuE.sum(); + derivs.dalphar_ddelta += (ndteuE*B_deltaE).sum()*one_over_delta; + derivs.dalphar_dtau += (ndteuE*B_tauE).sum()*one_over_tau; + derivs.d2alphar_ddelta2 += (ndteuE*B_delta2E).sum()*POW2(one_over_delta); + derivs.d2alphar_dtau2 += (ndteuE*B_tau2E).sum()*POW2(one_over_tau); + derivs.d2alphar_ddelta_dtau += (ndteuE*B_deltaE*B_tauE).sum()*one_over_delta*one_over_tau; + + derivs.d3alphar_ddelta3 += (ndteuE*B_delta3E).sum()*POW3(one_over_delta); + derivs.d3alphar_dtau3 += (ndteuE*B_tau3E).sum()*POW3(one_over_tau); + derivs.d3alphar_ddelta2_dtau += (ndteuE*B_delta2E*B_tauE).sum()*POW2(one_over_delta)*one_over_tau; + derivs.d3alphar_ddelta_dtau2 += (ndteuE*B_deltaE*B_tau2E).sum()*one_over_delta*POW2(one_over_tau); + + return; +}; + void ResidualHelmholtzGeneralizedExponential::all(const long double &tau, const long double &delta, HelmholtzDerivatives &derivs) throw() { long double log_tau = log(tau), log_delta = log(delta), ndteu, one_over_delta = 1/delta, one_over_tau = 1/tau; // division is much slower than multiplication, so do one division here + // Maybe split the construction of u and other parts into two separate loops? + // If both loops can get vectorized, could be worth it. const std::size_t N = elements.size(); - for (std::size_t i = 0; i < N; ++i) { ResidualHelmholtzGeneralizedExponentialElement &el = elements[i]; @@ -131,6 +238,8 @@ void ResidualHelmholtzGeneralizedExponential::all(const long double &tau, const derivs.d3alphar_ddelta2_dtau *= POW2(one_over_delta)*one_over_tau; derivs.d3alphar_ddelta_dtau2 *= one_over_delta*POW2(one_over_tau); + + return; }; diff --git a/src/main.cxx b/src/main.cxx index ab886c18..ad142724 100644 --- a/src/main.cxx +++ b/src/main.cxx @@ -506,6 +506,36 @@ int main() shared_ptr Water(new HelmholtzEOSMixtureBackend(names)); Water->set_mole_fractions(std::vector(1,1)); + ResidualHelmholtzGeneralizedExponential GenExp = Water->get_components()[0]->pEOS->alphar.GenExp; + + HelmholtzDerivatives derivs1, derivs2; + double tau = 0.8, delta = 1.1; + + ss = 0; + t1 = clock(); + for (long i = 0; i < N; ++i){ + derivs1.reset(); + GenExp.all(tau, delta+i*1e-18, derivs1); + ss += derivs1.alphar; + } + t2 = clock(); + std::cout << format("value(all): %0.13g, %0.13g, %g us/call\n", ss, derivs1.alphar, ((double)(t2-t1))/CLOCKS_PER_SEC/double(N)*1e6); + + ss = 0; + t1 = clock(); + for (long i = 0; i < N; ++i){ + derivs2.reset(); + GenExp.allEigen(tau, delta+i*1e-18, derivs2); + ss += derivs2.alphar; + } + t2 = clock(); + std::cout << format("value(allEigen): %0.13g, %0.13g, %g us/call\n", ss, derivs2.alphar, ((double)(t2-t1))/CLOCKS_PER_SEC/double(N)*1e6); + + int r44 =0; + } + #endif + #if 0 + { t1 = clock(); for (long i = 0; i < N; ++i){