From fc2e618fb7f62bb7c9516145d75063ae00ed9799 Mon Sep 17 00:00:00 2001 From: Ian Bell Date: Sat, 16 Aug 2014 17:50:04 +0200 Subject: [PATCH] Speed optimizations for saturation routines Signed-off-by: Ian Bell --- src/Backends/Helmholtz/FlashRoutines.cpp | 2 +- src/Backends/Helmholtz/VLERoutines.cpp | 12 +++--------- src/Helmholtz.cpp | 4 ++-- src/main.cxx | 2 +- 4 files changed, 7 insertions(+), 13 deletions(-) diff --git a/src/Backends/Helmholtz/FlashRoutines.cpp b/src/Backends/Helmholtz/FlashRoutines.cpp index 0835dc1f..dfd9f07d 100644 --- a/src/Backends/Helmholtz/FlashRoutines.cpp +++ b/src/Backends/Helmholtz/FlashRoutines.cpp @@ -167,7 +167,7 @@ void FlashRoutines::PQ_flash(HelmholtzEOSMixtureBackend &HEOS) // Use logarithm of delta as independent variables options.use_logdelta = false; - double increment = 0.2; + double increment = 0.4; try{ for (double omega = 1.0; omega > 0; omega -= increment){ diff --git a/src/Backends/Helmholtz/VLERoutines.cpp b/src/Backends/Helmholtz/VLERoutines.cpp index f41a9894..67578a07 100644 --- a/src/Backends/Helmholtz/VLERoutines.cpp +++ b/src/Backends/Helmholtz/VLERoutines.cpp @@ -136,12 +136,8 @@ void SaturationSolvers::saturation_P_pure_1D_T(HelmholtzEOSMixtureBackend *HEOS, double call(double T){ this->T = T; // Recalculate the densities using the current guess values - rhomolar_liq = HEOS->SatL->solver_rho_Tp(T, p, rhomolar_liq); - rhomolar_vap = HEOS->SatV->solver_rho_Tp(T, p, rhomolar_vap); - - // Set the densities in the saturation classes - HEOS->SatL->update(DmolarT_INPUTS, rhomolar_liq, T); - HEOS->SatV->update(DmolarT_INPUTS, rhomolar_vap, T); + HEOS->SatL->update_TP_guessrho(T, p, rhomolar_liq); + HEOS->SatV->update_TP_guessrho(T, p, rhomolar_vap); // Calculate the Gibbs functions for liquid and vapor gL = HEOS->SatL->gibbsmolar(); @@ -370,9 +366,7 @@ void SaturationSolvers::saturation_PHSU_pure(HelmholtzEOSMixtureBackend *HEOS, l throw SolutionError(format("saturation_PHSU_pure solver did not converge after 25 iterations for %s=%Lg current error is %Lg", info.c_str(), specified_value, error)); } } - while (error > 1e-10); - HEOS->SatL->update(DmolarT_INPUTS, rhoL, T); - HEOS->SatV->update(DmolarT_INPUTS, rhoV, T); + while (error > 1e-9); } void SaturationSolvers::saturation_D_pure(HelmholtzEOSMixtureBackend *HEOS, long double rhomolar, saturation_D_pure_options &options) { diff --git a/src/Helmholtz.cpp b/src/Helmholtz.cpp index 278b5a9a..feef56ac 100644 --- a/src/Helmholtz.cpp +++ b/src/Helmholtz.cpp @@ -219,7 +219,7 @@ void ResidualHelmholtzNonAnalytic::all(const long double &tau, const long double long double PI = 4*Bi*ai*(ai-1)*pow(pow(delta-1,2),ai-2)+2*pow(Ai/betai,2)*pow(pow(delta-1,2),1/betai-2)+4*Ai*theta/betai*(1/(2*betai)-1)*pow(pow(delta-1,2),1/(2*betai)-2); long double dPI_dDelta = -8*Bi*ai*(ai-1)*(ai-2)*pow(pow(delta-1,2),ai-5.0/2.0)-8*pow(Ai/betai,2)*(1/(2*betai)-1)*pow(pow(delta-1,2),1/betai-5.0/2.0)-(8*Ai*theta)/betai*(1/(2*betai)-1)*(1/(2*betai)-2)*pow(pow(delta-1,2),1/(2*betai)-5.0/2.0)+4*Ai/betai*(1/(2*betai)-1)*pow(pow(delta-1,2),1/(2*betai)-2)*dtheta_dDelta; dDELTA3_dDelta3 = 1/(delta-1)*dDELTA2_dDelta2-1/pow(delta-1,2)*dDELTA_dDelta+pow(delta-1,2)*dPI_dDelta+2*(delta-1)*PI; - dDELTAbi2_dDelta2 = bi*(pow(DELTA,bi-1.0)*dDELTA2_dDelta2+(bi-1.0)*pow(DELTA,bi-2.0)*pow(dDELTA_dDelta,2)); + dDELTAbi2_dDelta2 = bi*(pow(DELTA,bi-1)*dDELTA2_dDelta2+(bi-1.0)*pow(DELTA,bi-2.0)*pow(dDELTA_dDelta,2)); dDELTAbi3_dDelta3 = bi*(pow(DELTA,bi-1)*dDELTA3_dDelta3+dDELTA2_dDelta2*(bi-1)*pow(DELTA,bi-2)*dDELTA_dDelta+(bi-1)*(pow(DELTA,bi-2)*2*dDELTA_dDelta*dDELTA2_dDelta2+pow(dDELTA_dDelta,2)*(bi-2)*pow(DELTA,bi-3)*dDELTA_dDelta)); } @@ -227,7 +227,7 @@ void ResidualHelmholtzNonAnalytic::all(const long double &tau, const long double long double dDELTAbi2_dDelta_dTau=-Ai*bi*2.0/betai*pow(DELTA,bi-1.0)*(delta-1.0)*pow(pow(delta-1.0,2),1.0/(2.0*betai)-1.0)-2.0*theta*bi*(bi-1.0)*pow(DELTA,bi-2.0)*dDELTA_dDelta; long double dDELTAbi2_dTau2 = 2.0*bi*pow(DELTA,bi-1.0)+4.0*pow(theta,2)*bi*(bi-1.0)*pow(DELTA,bi-2.0); - long double dDELTAbi3_dTau3 = -12.0*theta*bi*(bi-1.0)*pow(DELTA,bi-2)-8*pow(theta,3)*bi*(bi-1)*(bi-2)*pow(DELTA,bi-3.0); + long double dDELTAbi3_dTau3 = -12.0*theta*bi*(bi-1.0)*pow(DELTA,bi-2)-8*pow(theta,3)*bi*(bi-1)*(bi-2)*pow(DELTA,bi-3); long double dDELTAbi3_dDelta_dTau2 = 2*bi*(bi-1)*pow(DELTA,bi-2)*dDELTA_dDelta+4*pow(theta,2)*bi*(bi-1)*(bi-2)*pow(DELTA,bi-3)*dDELTA_dDelta+8*theta*bi*(bi-1)*pow(DELTA,bi-2)*dtheta_dDelta; long double dDELTAbi3_dDelta2_dTau = bi*((bi-1)*pow(DELTA,bi-2)*dDELTA_dTau*dDELTA2_dDelta2 + pow(DELTA,bi-1)*dDELTA3_dDelta2_dTau+(bi-1)*((bi-2)*pow(DELTA,bi-3)*dDELTA_dTau*pow(dDELTA_dDelta,2)+pow(DELTA,bi-2)*2*dDELTA_dDelta*dDELTA2_dDelta_dTau)); diff --git a/src/main.cxx b/src/main.cxx index 169ad715..c51cf28f 100644 --- a/src/main.cxx +++ b/src/main.cxx @@ -502,7 +502,7 @@ int main() time_t t1,t2; long N = 100000; double ss = 0; - std::vector names(1,"Water"); + std::vector names(1,"Propane"); shared_ptr Water(new HelmholtzEOSMixtureBackend(names)); Water->set_mole_fractions(std::vector(1,1)); t1 = clock();