Cache derivative values of helmholtz energy for individual components; see #726

This commit is contained in:
Ian Bell
2015-08-06 12:58:22 -04:00
parent b585bc744b
commit a9af0df612
2 changed files with 87 additions and 21 deletions

View File

@@ -6,7 +6,7 @@
#include "rapidjson/rapidjson_include.h"
#include "Eigen/Core"
#include "time.h"
#include "CachedElement.h"
namespace CoolProp{
@@ -497,36 +497,99 @@ public:
class ResidualHelmholtzContainer
{
private:
CachedElement _base, _dDelta, _dTau, _dDelta2, _dTau2, _dDelta_dTau, _dDelta3, _dDelta2_dTau, _dDelta_dTau2, _dTau3;
CachedElement _dDelta4, _dDelta3_dTau, _dDelta2_dTau2, _dDelta_dTau3, _dTau4;
public:
ResidualHelmholtzNonAnalytic NonAnalytic;
ResidualHelmholtzSAFTAssociating SAFT;
ResidualHelmholtzGeneralizedExponential GenExp;
ResidualHelmholtzGeneralizedExponential GenExp;
HelmholtzDerivatives all(const CoolPropDbl tau, const CoolPropDbl delta)
HelmholtzDerivatives all(const CoolPropDbl tau, const CoolPropDbl delta, bool cache_values = false)
{
HelmholtzDerivatives derivs; // zeros out the elements
GenExp.all(tau, delta, derivs);
NonAnalytic.all(tau, delta, derivs);
SAFT.all(tau, delta, derivs);
if (cache_values){
_base = derivs.alphar;
_dDelta = derivs.dalphar_ddelta;
_dTau = derivs.dalphar_dtau;
_dDelta2 = derivs.d2alphar_ddelta2;
_dTau2 = derivs.d2alphar_dtau2;
_dDelta_dTau = derivs.d2alphar_ddelta_dtau;
_dDelta3 = derivs.d3alphar_ddelta3;
_dTau3 = derivs.d3alphar_dtau3;
_dDelta2_dTau = derivs.d3alphar_ddelta2_dtau;
_dDelta_dTau2 = derivs.d3alphar_ddelta_dtau2;
}
return derivs;
};
CoolPropDbl base(CoolPropDbl tau, CoolPropDbl delta) { return all(tau,delta).alphar; };
CoolPropDbl dDelta(CoolPropDbl tau, CoolPropDbl delta) { return all(tau,delta).dalphar_ddelta; };
CoolPropDbl dTau(CoolPropDbl tau, CoolPropDbl delta) { return all(tau, delta).dalphar_dtau; };
CoolPropDbl dDelta2(CoolPropDbl tau, CoolPropDbl delta) { return all(tau, delta).d2alphar_ddelta2; };
CoolPropDbl dDelta_dTau(CoolPropDbl tau, CoolPropDbl delta) { return all(tau, delta).d2alphar_ddelta_dtau; };
CoolPropDbl dTau2(CoolPropDbl tau, CoolPropDbl delta) { return all(tau, delta).d2alphar_dtau2; };
CoolPropDbl dDelta3(CoolPropDbl tau, CoolPropDbl delta) { return all(tau, delta).d3alphar_ddelta3; };
CoolPropDbl dDelta2_dTau(CoolPropDbl tau, CoolPropDbl delta) { return all(tau, delta).d3alphar_ddelta2_dtau; };
CoolPropDbl dDelta_dTau2(CoolPropDbl tau, CoolPropDbl delta) { return all(tau, delta).d3alphar_ddelta_dtau2; };
CoolPropDbl dTau3(CoolPropDbl tau, CoolPropDbl delta) { return all(tau, delta).d3alphar_dtau3; };
CoolPropDbl base(CoolPropDbl tau, CoolPropDbl delta) {
if (!_base)
return all(tau, delta).alphar;
else
return _base;
};
CoolPropDbl dDelta(CoolPropDbl tau, CoolPropDbl delta) {
if (!_dDelta)
return all(tau, delta).dalphar_ddelta;
else
return _dDelta;
};
CoolPropDbl dTau(CoolPropDbl tau, CoolPropDbl delta) {
if (!_dTau)
return all(tau, delta).dalphar_dtau;
else
return _dTau;
};
CoolPropDbl dDelta2(CoolPropDbl tau, CoolPropDbl delta) {
if (!_dDelta2)
return all(tau, delta).d2alphar_ddelta2;
else
return _dDelta2;
};
CoolPropDbl dDelta_dTau(CoolPropDbl tau, CoolPropDbl delta) {
if (!_dDelta_dTau)
return all(tau, delta).d2alphar_ddelta_dtau;
else
return _dDelta_dTau;
};
CoolPropDbl dTau2(CoolPropDbl tau, CoolPropDbl delta) {
if (!_dTau2)
return all(tau, delta).d2alphar_dtau2;
else
return _dTau2;
};
CoolPropDbl dDelta3(CoolPropDbl tau, CoolPropDbl delta) {
if (!_dDelta3)
return all(tau, delta).d3alphar_ddelta3;
else
return _dDelta3;
};
CoolPropDbl dDelta2_dTau(CoolPropDbl tau, CoolPropDbl delta) {
if (!_dDelta2_dTau)
return all(tau, delta).d3alphar_ddelta2_dtau;
else
return _dDelta2_dTau;
};
CoolPropDbl dDelta_dTau2(CoolPropDbl tau, CoolPropDbl delta) {
if (!_dDelta_dTau2)
return all(tau, delta).d3alphar_ddelta_dtau2;
else
return _dDelta_dTau2;
};
CoolPropDbl dTau3(CoolPropDbl tau, CoolPropDbl delta) {
if (!_dTau3)
return all(tau, delta).d3alphar_dtau3;
else
return _dTau3;
};
CoolPropDbl dDelta4(CoolPropDbl tau, CoolPropDbl delta) { return all(tau, delta).d4alphar_ddelta4; };
CoolPropDbl dDelta3_dTau(CoolPropDbl tau, CoolPropDbl delta) { return all(tau, delta).d4alphar_ddelta3_dtau; };
CoolPropDbl dDelta2_dTau2(CoolPropDbl tau, CoolPropDbl delta) { return all(tau, delta).d4alphar_ddelta2_dtau2; };
CoolPropDbl dDelta_dTau3(CoolPropDbl tau, CoolPropDbl delta) { return all(tau, delta).d4alphar_ddelta_dtau3; };
CoolPropDbl dTau4(CoolPropDbl tau, CoolPropDbl delta) { return all(tau, delta).d4alphar_dtau4; };
};
// #############################################################################

View File

@@ -2306,7 +2306,8 @@ void HelmholtzEOSMixtureBackend::calc_all_alphar_deriv_cache(const std::vector<C
summer_dTau3 = 0, summer_dDelta3 = 0, summer_dDelta2_dTau = 0, summer_dDelta_dTau2 = 0,
summer_dTau4 = 0, summer_dDelta4 = 0, summer_dDelta3_dTau = 0, summer_dDelta_dTau3 = 0, summer_dDelta2_dTau2 = 0;
for (std::size_t i = 0; i < N; ++i){
HelmholtzDerivatives derivs = components[i].EOS().alphar.all(tau, delta);
bool cache_values = true;
HelmholtzDerivatives derivs = components[i].EOS().alphar.all(tau, delta, cache_values);
CoolPropDbl xi = mole_fractions[i];
summer_base += xi*derivs.alphar;
@@ -2992,6 +2993,7 @@ public:
double delta, tau, M1_last, R_tau, R_delta;
std::vector<CoolProp::SimpleState> critical_points;
int N_critical_points;
Eigen::MatrixXd Lstar, adjLstar, dLstardTau, d2LstardTau2, dLstardDelta;
L0CurveTracer(HelmholtzEOSMixtureBackend &HEOS, double tau0, double delta0) : HEOS(HEOS), delta(delta0), tau(tau0), N_critical_points(0)
{
R_delta = 0.1; R_tau = 0.1;
@@ -3017,7 +3019,11 @@ public:
this->get_tau_delta(theta, tau, delta, tau_new, delta_new);
double rhomolar = HEOS.rhomolar_reducing()*delta_new, T = HEOS.T_reducing()/tau_new;
HEOS.update_DmolarT_direct(rhomolar, T);
double L1 = MixtureDerivatives::Lstar(HEOS, XN_INDEPENDENT).determinant();
Lstar = MixtureDerivatives::Lstar(HEOS, XN_INDEPENDENT);
adjLstar = adjugate(Lstar);
dLstardTau = MixtureDerivatives::dLstar_dX(HEOS, XN_INDEPENDENT, iTau);
dLstardDelta = MixtureDerivatives::dLstar_dX(HEOS, XN_INDEPENDENT, iDelta);
double L1 = Lstar.determinant();
//std::cout << "call: " << theta << " " << L1 << std::endl;
return L1;
};
@@ -3027,10 +3033,7 @@ public:
*/
double deriv(double theta){
std::size_t N = HEOS.get_mole_fractions().size();
Eigen::MatrixXd adjL = adjugate(MixtureDerivatives::Lstar(HEOS, XN_INDEPENDENT)),
dLdTau = MixtureDerivatives::dLstar_dX(HEOS, XN_INDEPENDENT, iTau),
dLdDelta = MixtureDerivatives::dLstar_dX(HEOS, XN_INDEPENDENT, iDelta);
double dL1_dtau = (adjL*dLdTau).trace(), dL1_ddelta = (adjL*dLdDelta).trace();
double dL1_dtau = (adjLstar*dLstardTau).trace(), dL1_ddelta = (adjLstar*dLstardDelta).trace();
return -R_tau*sin(theta)*dL1_dtau + R_delta*cos(theta)*dL1_ddelta;
};