github/CoolProp
1
1
Fork 0
mirror of https://github.com/CoolProp/CoolProp.git synced 2026-02-09 05:15:45 -05:00
Code Issues Packages Projects Releases Wiki Activity

SAFT has transitioned to parallel evaluation of derivatives, though it isn't truly parallel since there is some duplication

Browse Source 
Signed-off-by: Ian Bell <ian.h.bell@gmail.com>
This commit is contained in:
Ian Bell
2014-08-15 22:34:24 +02:00
parent 29632c72e6
commit eee8139306
3 changed files with 46 additions and 113 deletions
Download Patch File Download Diff File Expand all files Collapse all files

108
src/Helmholtz.cpp
View File

@@ -721,90 +721,32 @@ long double ResidualHelmholtzSAFTAssociating::d3g_deta3(const long double &eta)
long double ResidualHelmholtzSAFTAssociating::eta(const long double &delta){
return this->vbarn*delta;
}
long double ResidualHelmholtzSAFTAssociating::base(const long double &tau, const long double &delta) throw()
{
if (disabled){return 0;}
long double X = this->X(delta, this->Deltabar(tau, delta));
return this->m*this->a*((log(X)-X/2.0+0.5));
}
long double ResidualHelmholtzSAFTAssociating::dDelta(const long double &tau, const long double &delta) throw()
{
if (disabled){return 0;}
long double X = this->X(delta, this->Deltabar(tau, delta));
return this->m*this->a*(1/X-0.5)*this->dX_ddelta(tau, delta);
}
long double ResidualHelmholtzSAFTAssociating::dTau(const long double &tau, const long double &delta) throw()
{
if (disabled){return 0;}
long double X = this->X(delta, this->Deltabar(tau, delta));
return this->m*this->a*(1/X-0.5)*this->dX_dtau(tau, delta);
}
long double ResidualHelmholtzSAFTAssociating::dTau2(const long double &tau, const long double &delta) throw()
{
if (disabled){return 0;}
long double X = this->X(delta, this->Deltabar(tau, delta));
long double X_tau = this->dX_dtau(tau, delta);
long double X_tautau = this->d2X_dtau2(tau, delta);
return this->m*this->a*((1/X-0.5)*X_tautau-pow(X_tau/X, 2));
}
long double ResidualHelmholtzSAFTAssociating::dDelta2(const long double &tau, const long double &delta) throw()
{
if (disabled){return 0;}
long double X = this->X(delta, this->Deltabar(tau, delta));
long double X_delta = this->dX_ddelta(tau, delta);
long double X_deltadelta = this->d2X_ddelta2(tau, delta);
return this->m*this->a*((1/X-0.5)*X_deltadelta-pow(X_delta/X,2));
}
long double ResidualHelmholtzSAFTAssociating::dDelta_dTau(const long double &tau, const long double &delta) throw()
{
if (disabled){return 0;}
long double X = this->X(delta, this->Deltabar(tau, delta));
long double X_delta = this->dX_ddelta(tau, delta);
long double X_tau = this->dX_dtau(tau, delta);
long double X_deltatau = this->d2X_ddeltadtau(tau, delta);
return this->m*this->a*((-X_tau/X/X)*X_delta+X_deltatau*(1/X-0.5));
}
long double ResidualHelmholtzSAFTAssociating::dTau3(const long double &tau, const long double &delta) throw()
{
if (disabled){return 0;}
long double X = this->X(delta, this->Deltabar(tau, delta));
long double X_t = this->dX_dtau(tau, delta);
long double X_tt = this->d2X_dtau2(tau, delta);
long double X_ttt = this->d3X_dtau3(tau, delta);
return this->m*this->a*((1/X-1.0/2.0)*X_ttt+(-X_t/pow(X,(int)2))*X_tt-2*(pow(X,(int)2)*(X_t*X_tt)-pow(X_t,(int)2)*(X*X_t))/pow(X,(int)4));
}
long double ResidualHelmholtzSAFTAssociating::dDelta_dTau2(const long double &tau, const long double &delta) throw()
{
if (disabled){return 0;}
long double X = this->X(delta, this->Deltabar(tau, delta));
long double X_t = this->dX_dtau(tau, delta);
long double X_d = this->dX_ddelta(tau, delta);
long double X_tt = this->d2X_dtau2(tau, delta);
long double X_dt = this->d2X_ddeltadtau(tau, delta);
long double X_dtt = this->d3X_ddeltadtau2(tau, delta);
return this->m*this->a*((1/X-1.0/2.0)*X_dtt-X_d/pow(X,(int)2)*X_tt-2*(pow(X,(int)2)*(X_t*X_dt)-pow(X_t,(int)2)*(X*X_d))/pow(X,(int)4));
}
long double ResidualHelmholtzSAFTAssociating::dDelta2_dTau(const long double &tau, const long double &delta) throw()
{
if (disabled){return 0;}
long double X = this->X(delta, this->Deltabar(tau, delta));
long double X_t = this->dX_dtau(tau, delta);
long double X_d = this->dX_ddelta(tau, delta);
long double X_dd = this->d2X_ddelta2(tau, delta);
long double X_dt = this->d2X_ddeltadtau(tau, delta);
long double X_ddt = this->d3X_ddelta2dtau(tau, delta);
return this->m*this->a*((1/X-1.0/2.0)*X_ddt-X_t/pow(X,(int)2)*X_dd-2*(pow(X,(int)2)*(X_d*X_dt)-pow(X_d,(int)2)*(X*X_t))/pow(X,(int)4));
}
long double ResidualHelmholtzSAFTAssociating::dDelta3(const long double &tau, const long double &delta) throw()
{
if (disabled){return 0;}
long double X = this->X(delta, this->Deltabar(tau, delta));
long double X_d = this->dX_ddelta(tau, delta);
long double X_dd = this->d2X_ddelta2(tau, delta);
long double X_ddd = this->d3X_ddelta3(tau, delta);
return this->m*this->a*((1/X-1.0/2.0)*X_ddd-X_d/pow(X,(int)2)*X_dd-2*(pow(X,(int)2)*(X_d*X_dd)-pow(X_d,(int)2)*(X*X_d))/pow(X,(int)4));
}
void ResidualHelmholtzSAFTAssociating::all(const long double &tau, const long double &delta, Derivatives &deriv) throw()
{
if (disabled){return;}
long double X = this->X(delta, this->Deltabar(tau, delta));
long double X_t = this->dX_dtau(tau, delta);
long double X_d = this->dX_ddelta(tau, delta);
long double X_tt = this->d2X_dtau2(tau, delta);
long double X_dd = this->d2X_ddelta2(tau, delta);
long double X_dt = this->d2X_ddeltadtau(tau, delta);
long double X_ttt = this->d3X_dtau3(tau, delta);
long double X_dtt = this->d3X_ddeltadtau2(tau, delta);
long double X_ddt = this->d3X_ddelta2dtau(tau, delta);
long double X_ddd = this->d3X_ddelta3(tau, delta);
deriv.alphar += this->m*this->a*((log(X)-X/2.0+0.5));
deriv.dalphar_ddelta += this->m*this->a*(1/X-0.5)*this->dX_ddelta(tau, delta);
deriv.dalphar_dtau += this->m*this->a*(1/X-0.5)*this->dX_dtau(tau, delta);
deriv.d2alphar_dtau2 += this->m*this->a*((1/X-0.5)*X_tt-pow(X_t/X, 2));
deriv.d2alphar_ddelta2 += this->m*this->a*((1/X-0.5)*X_dd-pow(X_d/X,2));
deriv.d2alphar_ddelta_dtau += this->m*this->a*((-X_t/X/X)*X_d + X_dt*(1/X-0.5));
deriv.d3alphar_dtau3 += this->m*this->a*((1/X-1.0/2.0)*X_ttt+(-X_t/pow(X,(int)2))*X_tt-2*(pow(X,(int)2)*(X_t*X_tt)-pow(X_t,(int)2)*(X*X_t))/pow(X,(int)4));
deriv.d3alphar_ddelta_dtau2 += this->m*this->a*((1/X-1.0/2.0)*X_dtt-X_d/pow(X,(int)2)*X_tt-2*(pow(X,(int)2)*(X_t*X_dt)-pow(X_t,(int)2)*(X*X_d))/pow(X,(int)4));
deriv.d3alphar_ddelta2_dtau += this->m*this->a*((1/X-1.0/2.0)*X_ddt-X_t/pow(X,(int)2)*X_dd-2*(pow(X,(int)2)*(X_d*X_dt)-pow(X_d,(int)2)*(X*X_t))/pow(X,(int)4));
deriv.d3alphar_ddelta3 += this->m*this->a*((1/X-1.0/2.0)*X_ddd-X_d/pow(X,(int)2)*X_dd-2*(pow(X,(int)2)*(X_d*X_dd)-pow(X_d,(int)2)*(X*X_d))/pow(X,(int)4));
}
void IdealHelmholtzCP0PolyT::to_json(rapidjson::Value &el, rapidjson::Document &doc)
{

29
src/main.cxx
View File

@@ -163,7 +163,7 @@ int main()
double rrv = CoolProp::PropsSI("C","T",300,"D",1e-10,"REFPROP::R1234ze");
std::cout << get_global_param_string("errstring");
double rr0 = HumidAir::HAPropsSI("B","T",473.15,"W",0.5,"P",101325);
double rr0 = HumidAir::HAPropsSI("S","T",473.15,"W",0,"P",1e6);
//CoolProp::set_reference_stateS("Air","RESET");
double rr1 = HumidAir::HAPropsSI("B","T",473.15,"W",0.5,"P",101325);
int r = 1;
@@ -505,34 +505,23 @@ int main()
Derivatives derivs;
time_t t1,t2;
long N = 10000;
long N = 100000;
double ss = 0;
std::vector<CoolPropFluid*> components = Water->get_components();
ResidualHelmholtzPower Power = components[0]->pEOS->alphar.Power;
ResidualHelmholtzGeneralizedExponential GenExp = components[0]->pEOS->alphar.GenExp;
long double tau = 0.8, delta = 2.2;
GenExp.all(tau, delta, derivs);
t1 = clock();
for (long i = 0; i < N; ++i){
Power.all(tau, delta+i*1e-10, derivs);
GenExp.all(tau, delta+i*1e-10, derivs);
ss += derivs.alphar+derivs.dalphar_ddelta+derivs.dalphar_dtau+derivs.d2alphar_ddelta2+derivs.d2alphar_ddelta_dtau+derivs.d2alphar_dtau2;
}
t2 = clock();
std::cout << format("value: %0.13g, %g us/call", ss, ((double)(t2-t1))/CLOCKS_PER_SEC/double(N)*1e6);
ss = 0;
t1 = clock();
for (long i = 0; i < N; ++i){
long double alphar = Power.base(tau, delta+i*1e-10);
long double dalphar_ddelta = Power.dDelta(tau, delta+i*1e-10);
long double dalphar_dtau = Power.dTau(tau, delta+i*1e-10);
long double d2alphar_ddelta2 = Power.dDelta2(tau, delta+i*1e-10);
long double d2alphar_ddelta_dtau = Power.dDelta_dTau(tau, delta+i*1e-10);
long double d2alphar_dtau2 = Power.dTau2(tau, delta+i*1e-10);
ss += alphar+dalphar_ddelta+dalphar_dtau+d2alphar_ddelta2+d2alphar_ddelta_dtau+d2alphar_dtau2;
}
t2 = clock();
std::cout << format("value: %0.13g, %g us/call", ss, ((double)(t2-t1))/CLOCKS_PER_SEC/double(N)*1e6);
std::cout << format("value: %0.13g, %g us/call\n", ss, ((double)(t2-t1))/CLOCKS_PER_SEC/double(N)*1e6);
exit(0);
}