diff --git a/include/Helmholtz.h b/include/Helmholtz.h index 124ba32a..10c81572 100644 --- a/include/Helmholtz.h +++ b/include/Helmholtz.h @@ -377,16 +377,18 @@ public: void to_json(rapidjson::Value &el, rapidjson::Document &doc); - long double base(const long double &tau, const long double &delta) throw(); - long double dDelta(const long double &tau, const long double &delta) throw(); - long double dTau(const long double &tau, const long double &delta) throw(); - long double dDelta2(const long double &tau, const long double &delta) throw(); - long double dDelta_dTau(const long double &tau, const long double &delta) throw(); - long double dTau2(const long double &tau, const long double &delta) throw(); - long double dDelta3(const long double &tau, const long double &delta) throw(); - 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(); + long double base(const long double &tau, const long double &delta) throw(){Derivatives deriv; all(tau,delta,deriv); return deriv.alphar;}; + long double dDelta(const long double &tau, const long double &delta) throw(){Derivatives deriv; all(tau,delta,deriv); return deriv.dalphar_ddelta;}; + long double dTau(const long double &tau, const long double &delta) throw(){Derivatives deriv; all(tau,delta,deriv); return deriv.dalphar_dtau;}; + long double dDelta2(const long double &tau, const long double &delta) throw(){Derivatives deriv; all(tau,delta,deriv); return deriv.d2alphar_ddelta2;}; + long double dDelta_dTau(const long double &tau, const long double &delta) throw(){Derivatives deriv; all(tau,delta,deriv); return deriv.d2alphar_ddelta_dtau;}; + long double dTau2(const long double &tau, const long double &delta) throw(){Derivatives deriv; all(tau,delta,deriv); return deriv.d2alphar_dtau2;}; + long double dDelta3(const long double &tau, const long double &delta) throw(){Derivatives deriv; all(tau,delta,deriv); return deriv.d3alphar_ddelta3;}; + long double dDelta2_dTau(const long double &tau, const long double &delta) throw(){Derivatives deriv; all(tau,delta,deriv); return deriv.d3alphar_ddelta2_dtau;}; + long double dDelta_dTau2(const long double &tau, const long double &delta) throw(){Derivatives deriv; all(tau,delta,deriv); return deriv.d3alphar_ddelta_dtau2;}; + long double dTau3(const long double &tau, const long double &delta) throw(){Derivatives deriv; all(tau,delta,deriv); return deriv.d3alphar_dtau3;}; + + void all(const long double &tau, const long double &delta, Derivatives &deriv) throw(); }; class ResidualHelmholtzContainer diff --git a/src/Helmholtz.cpp b/src/Helmholtz.cpp index 8ccaa3ad..edfc6b29 100644 --- a/src/Helmholtz.cpp +++ b/src/Helmholtz.cpp @@ -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) { diff --git a/src/main.cxx b/src/main.cxx index 4031c4d1..80b0d302 100644 --- a/src/main.cxx +++ b/src/main.cxx @@ -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 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); }