github/CoolProp
1
1
Fork 0
mirror of https://github.com/CoolProp/CoolProp.git synced 2026-04-01 03:00:13 -04:00
Code Issues Packages Projects Releases Wiki Activity

Implemented Eigen solution for all() function, speed testing next

Browse Source 
Signed-off-by: Ian Bell <ian.h.bell@gmail.com>
This commit is contained in:
Ian Bell
2014-08-17 21:09:27 +02:00
parent 10f23992e7
commit f7464fcaa2
5 changed files with 193 additions and 10 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
dev/codelite/coolprop.project
View File

@@ -190,7 +190,7 @@
</Completion>
</Configuration>
<Configuration Name="Release" CompilerType="MinGW ( CodeLite-4.8.1 )" DebuggerType="GNU gdb debugger" Type="Executable" BuildCmpWithGlobalSettings="append" BuildLnkWithGlobalSettings="append" BuildResWithGlobalSettings="append">
<Compiler Options="-O2;-Wall" C_Options="-O2;-Wall" Assembler="" Required="yes" PreCompiledHeader="" PCHInCommandLine="no" UseDifferentPCHFlags="no" PCHFlags="">
<Compiler Options="-O3;-Wall;-msse2" C_Options="-O2;-Wall" Assembler="" Required="yes" PreCompiledHeader="" PCHInCommandLine="no" UseDifferentPCHFlags="no" PCHFlags="">
<IncludePath Value="."/>
<IncludePath Value="../../include"/>
<IncludePath Value="../../externals/Eigen"/>

57
include/Helmholtz.h
View File

@@ -4,6 +4,7 @@
#include <vector>
#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<long double> 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<double> 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<int> l_int, m_int;
Eigen::ArrayXd uE, du_ddeltaE, du_dtauE, d2u_ddelta2E, d2u_dtau2E, d3u_ddelta3E, d3u_dtau3E;
std::vector<ResidualHelmholtzGeneralizedExponentialElement> 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<long double> &n, const std::vector<long double> &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

3
src/Backends/Helmholtz/Fluids/FluidLibrary.h
View File

@@ -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

111
src/Helmholtz.cpp
View File

@@ -1,6 +1,7 @@
#include <numeric>
#include "Helmholtz.h"
namespace CoolProp{
long double kahanSum(std::vector<long double> &x)
@@ -18,13 +19,119 @@ long double kahanSum(std::vector<long double> &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<Eigen::ArrayXd> nE(n.data(), elements.size());
Eigen::Map<Eigen::ArrayXd> dE(d.data(), elements.size());
Eigen::Map<Eigen::ArrayXd> tE(t.data(), elements.size());
Eigen::Map<Eigen::ArrayXd> cE(c.data(), elements.size());
Eigen::Map<Eigen::ArrayXi> l_intE(l_int.data(), elements.size());
Eigen::Map<Eigen::ArrayXd> l_doubleE(l_double.data(), elements.size());
Eigen::Map<Eigen::ArrayXd> eta1E(eta1.data(), elements.size());
Eigen::Map<Eigen::ArrayXd> eta2E(eta2.data(), elements.size());
Eigen::Map<Eigen::ArrayXd> epsilon1E(epsilon1.data(), elements.size());
Eigen::Map<Eigen::ArrayXd> epsilon2E(epsilon2.data(), elements.size());
Eigen::Map<Eigen::ArrayXd> beta1E(beta1.data(), elements.size());
Eigen::Map<Eigen::ArrayXd> beta2E(beta2.data(), elements.size());
Eigen::Map<Eigen::ArrayXd> gamma1E(gamma1.data(), elements.size());
Eigen::Map<Eigen::ArrayXd> 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;
};

30
src/main.cxx
View File

@@ -506,6 +506,36 @@ int main()
shared_ptr<HelmholtzEOSMixtureBackend> Water(new HelmholtzEOSMixtureBackend(names));
Water->set_mole_fractions(std::vector<long double>(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){