github/CoolProp
1
1
Fork 0
mirror of https://github.com/CoolProp/CoolProp.git synced 2026-01-22 20:38:01 -05:00
Code Issues Packages Projects Releases Wiki Activity

Massive speed increase in tests for mixture derivatives

Browse Source
This commit is contained in:
Ian Bell
2015-05-29 22:34:24 -06:00
parent a58a64658d
commit f71313aff3
1 changed files with 167 additions and 146 deletions
Download Patch File Download Diff File Expand all files Collapse all files

313
src/Backends/Helmholtz/MixtureDerivatives.cpp
View File

@@ -380,8 +380,132 @@ CoolPropDbl MixtureDerivatives::d_ndalphardni_dTau(HelmholtzEOSMixtureBackend &H
#include "catch.hpp"
using namespace CoolProp;
static bool fluids_set = false;
static const std::size_t Ncomp_max = 6;
// These are composition invariant
// ** Levels **
// 0: number of components in the mixture
// 1: component index
static std::vector<std::vector<shared_ptr<HelmholtzEOSMixtureBackend> > > HEOS,
HEOS_plusT_constrho, HEOS_minusT_constrho,
HEOS_plusrho_constT, HEOS_minusrho_constT,
HEOS_plusz_xNindep, HEOS_minusz_xNindep,
HEOS_plusz_xNdep, HEOS_minusz_xNdep,
HEOS_plusz_consttaudelta_xNindep, HEOS_minusz_consttaudelta_xNindep,
HEOS_plusz_consttaudelta_xNdep, HEOS_minusz_consttaudelta_xNdep;
static const double T1 = 300, rho1 = 300, dT = 1e-3, drho = 1e-3, dz = 1e-6;
void setup_state(std::vector<shared_ptr<HelmholtzEOSMixtureBackend> > & HEOS, std::size_t Ncomp, double increment, x_N_dependency_flag xN_flag = XN_INDEPENDENT)
{
std::vector<std::string> names(Ncomp);
std::vector<CoolPropDbl> z(Ncomp);
if (Ncomp == 2){
names[0] = "Ethane"; names[1] = "Propane";
z[0] = 0.3; z[1] = 0.7;
}
else if (Ncomp == 3){
names[0] = "Ethane"; names[1] = "Propane"; names[2] = "Methane";
z[0] = 0.3; z[1] = 0.4; z[2] = 0.3;
}
else if (Ncomp == 4){
names[0] = "Ethane"; names[1] = "Propane"; names[2] = "Methane"; names[3] = "n-Butane";
z[0] = 0.3; z[1] = 0.4; z[2] = 0.2; z[3] = 0.1;
}
for (std::size_t i = 0; i < HEOS.size(); ++i){
std::vector<CoolPropDbl> zn = z;
zn[i] += increment;
if (xN_flag == XN_DEPENDENT){ zn[zn.size()-1] -= increment; }
HEOS[i].reset(new HelmholtzEOSMixtureBackend(names));
HEOS[i]->specify_phase(iphase_gas);
HEOS[i]->set_mole_fractions(zn);
}
}
// Set up all the fluids
void connect_fluids(){
if (!fluids_set){
HEOS.resize(Ncomp_max);
HEOS_plusT_constrho.resize(Ncomp_max);
HEOS_minusT_constrho.resize(Ncomp_max);
HEOS_plusrho_constT.resize(Ncomp_max);
HEOS_minusrho_constT.resize(Ncomp_max);
HEOS_plusz_xNindep.resize(Ncomp_max);
HEOS_minusz_xNindep.resize(Ncomp_max);
HEOS_plusz_consttaudelta_xNindep.resize(Ncomp_max);
HEOS_minusz_consttaudelta_xNindep.resize(Ncomp_max);
HEOS_plusz_xNdep.resize(Ncomp_max);
HEOS_minusz_xNdep.resize(Ncomp_max);
HEOS_plusz_consttaudelta_xNdep.resize(Ncomp_max);
HEOS_minusz_consttaudelta_xNdep.resize(Ncomp_max);
for (std::size_t Ncomp = 2; Ncomp <= 4; ++Ncomp){
HEOS[Ncomp].resize(1);
HEOS_plusT_constrho[Ncomp].resize(1);
HEOS_minusT_constrho[Ncomp].resize(1);
HEOS_plusrho_constT[Ncomp].resize(1);
HEOS_minusrho_constT[Ncomp].resize(1);
HEOS_plusz_xNindep[Ncomp].resize(Ncomp);
HEOS_minusz_xNindep[Ncomp].resize(Ncomp);
HEOS_plusz_consttaudelta_xNindep[Ncomp].resize(Ncomp);
HEOS_minusz_consttaudelta_xNindep[Ncomp].resize(Ncomp);
HEOS_plusz_xNdep[Ncomp].resize(Ncomp);
HEOS_minusz_xNdep[Ncomp].resize(Ncomp);
HEOS_plusz_consttaudelta_xNdep[Ncomp].resize(Ncomp);
HEOS_minusz_consttaudelta_xNdep[Ncomp].resize(Ncomp);
setup_state(HEOS[Ncomp], Ncomp, 0);
setup_state(HEOS_plusT_constrho[Ncomp], Ncomp, 0);
setup_state(HEOS_minusT_constrho[Ncomp], Ncomp, 0);
setup_state(HEOS_plusrho_constT[Ncomp], Ncomp, 0);
setup_state(HEOS_minusrho_constT[Ncomp], Ncomp, 0);
setup_state(HEOS_plusz_xNindep[Ncomp], Ncomp, dz);
setup_state(HEOS_minusz_xNindep[Ncomp], Ncomp, -dz);
setup_state(HEOS_plusz_consttaudelta_xNindep[Ncomp], Ncomp, dz);
setup_state(HEOS_minusz_consttaudelta_xNindep[Ncomp], Ncomp, -dz);
setup_state(HEOS_plusz_xNdep[Ncomp], Ncomp, dz, XN_DEPENDENT);
setup_state(HEOS_minusz_xNdep[Ncomp], Ncomp, -dz, XN_DEPENDENT);
setup_state(HEOS_plusz_consttaudelta_xNdep[Ncomp], Ncomp, dz, XN_DEPENDENT);
setup_state(HEOS_minusz_consttaudelta_xNdep[Ncomp], Ncomp, -dz, XN_DEPENDENT);
HEOS[Ncomp][0]->update(DmolarT_INPUTS, rho1, T1);
HEOS_plusT_constrho[Ncomp][0]->update(DmolarT_INPUTS, rho1, T1 + dT);
HEOS_minusT_constrho[Ncomp][0]->update(DmolarT_INPUTS, rho1, T1 - dT);
HEOS_plusrho_constT[Ncomp][0]->update(DmolarT_INPUTS, rho1 + drho, T1);
HEOS_minusrho_constT[Ncomp][0]->update(DmolarT_INPUTS, rho1 - drho, T1);
for (std::size_t i = 0; i < Ncomp; ++i){
HEOS_plusz_xNindep[Ncomp][i]->update(DmolarT_INPUTS, rho1, T1);
HEOS_minusz_xNindep[Ncomp][i]->update(DmolarT_INPUTS, rho1, T1);
HEOS_plusz_xNdep[Ncomp][i]->update(DmolarT_INPUTS, rho1, T1);
HEOS_minusz_xNdep[Ncomp][i]->update(DmolarT_INPUTS, rho1, T1);
HEOS_plusz_consttaudelta_xNindep[Ncomp][i]->calc_reducing_state();
SimpleState red = HEOS_plusz_consttaudelta_xNindep[Ncomp][i]->get_reducing_state();
HEOS_plusz_consttaudelta_xNindep[Ncomp][i]->update(DmolarT_INPUTS, red.rhomolar*HEOS[Ncomp][0]->delta(), red.T/HEOS[Ncomp][0]->tau());
HEOS_plusz_consttaudelta_xNdep[Ncomp][i]->calc_reducing_state();
red = HEOS_plusz_consttaudelta_xNdep[Ncomp][i]->get_reducing_state();
HEOS_plusz_consttaudelta_xNdep[Ncomp][i]->update(DmolarT_INPUTS, red.rhomolar*HEOS[Ncomp][0]->delta(), red.T/HEOS[Ncomp][0]->tau());
HEOS_minusz_consttaudelta_xNindep[Ncomp][i]->calc_reducing_state();
red = HEOS_minusz_consttaudelta_xNindep[Ncomp][i]->get_reducing_state();
HEOS_minusz_consttaudelta_xNindep[Ncomp][i]->update(DmolarT_INPUTS, red.rhomolar*HEOS[Ncomp][0]->delta(), red.T/HEOS[Ncomp][0]->tau());
HEOS_minusz_consttaudelta_xNdep[Ncomp][i]->calc_reducing_state();
red = HEOS_minusz_consttaudelta_xNdep[Ncomp][i]->get_reducing_state();
HEOS_minusz_consttaudelta_xNdep[Ncomp][i]->update(DmolarT_INPUTS, red.rhomolar*HEOS[Ncomp][0]->delta(), red.T/HEOS[Ncomp][0]->tau());
}
}
fluids_set = true;
}
}
TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]")
{
connect_fluids();
for (std::size_t Ncomp = 2; Ncomp <= 4; Ncomp++)
{
std::ostringstream ss00;
@@ -403,86 +527,24 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]")
ss000 << xN_string;
SECTION(ss000.str(),"")
{
std::vector<std::string> names;
std::vector<CoolPropDbl> z;
shared_ptr<HelmholtzEOSMixtureBackend> HEOS, HEOS_plusT_constrho, HEOS_plusrho_constT, HEOS_minusT_constrho, HEOS_minusrho_constT;
names.resize(Ncomp);
z.resize(Ncomp);
{
/* Set up a test class for the mixture tests */
if (Ncomp == 2)
{
names[0] = "Ethane"; names[1] = "Propane";
z[0] = 0.3; z[1] = 0.7;
}
else if (Ncomp == 3)
{
names[0] = "Ethane"; names[1] = "Propane"; names[2] = "Methane";
z[0] = 0.3; z[1] = 0.4; z[2] = 0.3;
}
else if (Ncomp == 4)
{
names[0] = "Ethane"; names[1] = "Propane"; names[2] = "Methane"; names[3] = "n-Butane";
z[0] = 0.3; z[1] = 0.4; z[2] = 0.2; z[3] = 0.1;
}
double T1 = 300, rho1 = 300, dT = 1e-3, drho = 1e-3, dz = 1e-6;
HEOS.reset(new HelmholtzEOSMixtureBackend(names));
HEOS->specify_phase(iphase_gas);
HelmholtzEOSMixtureBackend &rHEOS = *(HEOS.get());
rHEOS.set_mole_fractions(z);
rHEOS.update(DmolarT_INPUTS, rho1, T1);
HEOS_plusT_constrho.reset(new HelmholtzEOSMixtureBackend(names));
HEOS_plusT_constrho->specify_phase(iphase_gas);
HelmholtzEOSMixtureBackend &rHEOS_plusT_constrho = *(HEOS_plusT_constrho.get());
rHEOS_plusT_constrho.set_mole_fractions(z);
rHEOS_plusT_constrho.update(DmolarT_INPUTS, rho1, T1 + dT);
HEOS_minusT_constrho.reset(new HelmholtzEOSMixtureBackend(names));
HEOS_minusT_constrho->specify_phase(iphase_gas);
HelmholtzEOSMixtureBackend &rHEOS_minusT_constrho = *(HEOS_minusT_constrho.get());
rHEOS_minusT_constrho.set_mole_fractions(z);
rHEOS_minusT_constrho.update(DmolarT_INPUTS, rho1, T1 - dT);
HEOS_plusrho_constT.reset(new HelmholtzEOSMixtureBackend(names));
HEOS_plusrho_constT->specify_phase(iphase_gas);
HelmholtzEOSMixtureBackend &rHEOS_plusrho_constT = *(HEOS_plusrho_constT.get());
rHEOS_plusrho_constT.set_mole_fractions(z);
rHEOS_plusrho_constT.update(DmolarT_INPUTS, rho1 + drho, T1);
HEOS_minusrho_constT.reset(new HelmholtzEOSMixtureBackend(names));
HEOS_minusrho_constT->specify_phase(iphase_gas);
HelmholtzEOSMixtureBackend &rHEOS_minusrho_constT = *(HEOS_minusrho_constT.get());
rHEOS_minusrho_constT.set_mole_fractions(z);
rHEOS_minusrho_constT.update(DmolarT_INPUTS, rho1 - drho, T1);
rHEOS.update(DmolarT_INPUTS, rho1, T1);
HelmholtzEOSMixtureBackend &rHEOS = *(HEOS[Ncomp][0].get());
HelmholtzEOSMixtureBackend &rHEOS_plusT_constrho = *(HEOS_plusT_constrho[Ncomp][0].get());
HelmholtzEOSMixtureBackend &rHEOS_minusT_constrho = *(HEOS_minusT_constrho[Ncomp][0].get());
HelmholtzEOSMixtureBackend &rHEOS_plusrho_constT = *(HEOS_plusrho_constT[Ncomp][0].get());
HelmholtzEOSMixtureBackend &rHEOS_minusrho_constT = *(HEOS_minusrho_constT[Ncomp][0].get());
const std::vector<CoolPropDbl> &z = rHEOS.get_mole_fractions();
// These ones only require the i index
for (std::size_t i = 0; i< z.size();++i)
for (std::size_t i = 0; i< Ncomp; ++i)
{
shared_ptr<HelmholtzEOSMixtureBackend> HEOS, HEOS_pluszi, HEOS_minuszi;
HEOS_pluszi.reset(new HelmholtzEOSMixtureBackend(names));
HEOS_pluszi->specify_phase(iphase_gas);
HelmholtzEOSMixtureBackend &rHEOS_pluszi = *(HEOS_pluszi.get());
std::vector<CoolPropDbl> zp = z; /// Copy base composition
zp[i] += dz;
if (xN_flag == XN_DEPENDENT)
zp[z.size()-1] -= dz;
rHEOS_pluszi.set_mole_fractions(zp);
rHEOS_pluszi.update(DmolarT_INPUTS, rho1, T1);
HEOS_minuszi.reset(new HelmholtzEOSMixtureBackend(names));
HEOS_minuszi->specify_phase(iphase_gas);
HelmholtzEOSMixtureBackend &rHEOS_minuszi = *(HEOS_minuszi.get());
std::vector<CoolPropDbl> zm = z; /// Copy base composition
zm[i] -= dz;
if (xN_flag == XN_DEPENDENT)
zm[z.size()-1] += dz;
rHEOS_minuszi.set_mole_fractions(zm);
rHEOS_minuszi.update(DmolarT_INPUTS, rho1, T1);
HelmholtzEOSMixtureBackend & rHEOS_pluszi = (xN_flag == XN_INDEPENDENT) ? *(HEOS_plusz_xNindep[Ncomp][i].get()) : *(HEOS_plusz_xNdep[Ncomp][i].get());
HelmholtzEOSMixtureBackend & rHEOS_minuszi = (xN_flag == XN_INDEPENDENT) ? *(HEOS_minusz_xNindep[Ncomp][i].get()) : *(HEOS_minusz_xNdep[Ncomp][i].get());
HelmholtzEOSMixtureBackend &rHEOS_pluszi_consttaudelta = (xN_flag == XN_INDEPENDENT) ? *(HEOS_plusz_consttaudelta_xNindep[Ncomp][i].get()) : *(HEOS_plusz_consttaudelta_xNdep[Ncomp][i].get());
HelmholtzEOSMixtureBackend &rHEOS_minuszi_consttaudelta = (xN_flag == XN_INDEPENDENT) ? *(HEOS_minusz_consttaudelta_xNindep[Ncomp][i].get()) : *(HEOS_minusz_consttaudelta_xNdep[Ncomp][i].get());
std::ostringstream ss0;
ss0 << "dln_fugacity_i_dT__constrho_n, i=" << i;
@@ -534,7 +596,7 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]")
double err = std::abs((numeric-analytic)/analytic);
CHECK(err < 1e-8);
}
std::ostringstream ss1;
/*std::ostringstream ss1;
ss1 << "dln_fugacity_coefficient_dT__constp_n, i=" << i;
SECTION(ss1.str(), "")
{
@@ -572,7 +634,7 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]")
double err = std::abs((numeric-analytic)/analytic);
CHECK(err < 1e-6);
}
*/
std::ostringstream ss3;
ss3 << "d_ndalphardni_dDelta, i=" << i;
SECTION(ss3.str(), "")
@@ -588,7 +650,7 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]")
ss3a << "d2alphar_dxi_dDelta, i=" << i;
SECTION(ss3a.str(), "")
{
if (i == z.size()-1){break;}
if (i == Ncomp-1){ break; }
double analytic = MixtureDerivatives::d2alphar_dxi_dDelta(rHEOS, i, xN_flag);
double v1 = MixtureDerivatives::dalphar_dxi(rHEOS_plusrho_constT, i, xN_flag), delta1 = rHEOS_plusrho_constT.delta();
double v2 = MixtureDerivatives::dalphar_dxi(rHEOS_minusrho_constT, i, xN_flag), delta2 = rHEOS_minusrho_constT.delta();
@@ -600,7 +662,7 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]")
ss4a << "d2alphar_dxi_dTau, i=" << i;
SECTION(ss4a.str(), "")
{
if (i == z.size()-1){ break; }
if (i == Ncomp-1){ break; }
double analytic = MixtureDerivatives::d2alphar_dxi_dTau(rHEOS, i, xN_flag);
double v1 = MixtureDerivatives::dalphar_dxi(rHEOS_plusT_constrho, i, xN_flag), tau1 = rHEOS_plusT_constrho.tau();
double v2 = MixtureDerivatives::dalphar_dxi(rHEOS_minusT_constrho, i, xN_flag), tau2 = rHEOS_minusT_constrho.tau();
@@ -623,7 +685,7 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]")
ss5 << "dpdxj__constT_V_xi, i=" << i;
SECTION(ss5.str(), "")
{
if (i==z.size()-1){break;}
if (i==Ncomp-1){ break; }
double analytic = MixtureDerivatives::dpdxj__constT_V_xi(rHEOS, i, xN_flag);
double v1 = rHEOS_pluszi.p();
double v2 = rHEOS_minuszi.p();
@@ -637,7 +699,7 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]")
ss5a << "dtaudxj__constT_V_xi, i=" << i;
SECTION(ss5a.str(), "")
{
if (i==z.size()-1){break;}
if (i==Ncomp-1){ break; }
double analytic = MixtureDerivatives::dtau_dxj__constT_V_xi(rHEOS, i, xN_flag);
double v1 = rHEOS_pluszi.tau();
double v2 = rHEOS_minuszi.tau();
@@ -651,7 +713,7 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]")
ss5b << "ddeltadxj__constT_V_xi, i=" << i;
SECTION(ss5b.str(), "")
{
if (i==z.size()-1){break;}
if (i==Ncomp-1){ break; }
double analytic = MixtureDerivatives::ddelta_dxj__constT_V_xi(rHEOS, i, xN_flag);
double v1 = rHEOS_pluszi.delta();
double v2 = rHEOS_minuszi.delta();
@@ -665,7 +727,7 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]")
ss6 << "d_dalpharddelta_dxj__constT_V_xi, i=" << i;
SECTION(ss6.str(), "")
{
if (i==z.size()-1){break;}
if (i==Ncomp-1){ break; }
double analytic = MixtureDerivatives::d_dalpharddelta_dxj__constT_V_xi(rHEOS, i, xN_flag);
double v1 = rHEOS_pluszi.dalphar_dDelta();
double v2 = rHEOS_minuszi.dalphar_dDelta();
@@ -679,7 +741,7 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]")
ss7 << "dTrdxi__constxj, i=" << i;
SECTION(ss7.str(), "")
{
if (i == z.size()-1){break;}
if (i == Ncomp-1){ break; }
double analytic = rHEOS.Reducing->dTrdxi__constxj(rHEOS.get_mole_fractions(), i, xN_flag);
double v1 = rHEOS_pluszi.Reducing->Tr(rHEOS_pluszi.get_mole_fractions());
double v2 = rHEOS_minuszi.Reducing->Tr(rHEOS_minuszi.get_mole_fractions());
@@ -693,7 +755,7 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]")
ss8 << "drhormolardxi__constxj, i=" << i;
SECTION(ss8.str(), "")
{
if (i == z.size()-1){break;}
if (i == Ncomp-1){ break; }
double analytic = rHEOS.Reducing->drhormolardxi__constxj(rHEOS.get_mole_fractions(), i, xN_flag);
double v1 = rHEOS_pluszi.Reducing->rhormolar(rHEOS_pluszi.get_mole_fractions());
double v2 = rHEOS_minuszi.Reducing->rhormolar(rHEOS_minuszi.get_mole_fractions());
@@ -707,7 +769,7 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]")
ss3c << "d2Trdxi2__constxj, i=" << i;
SECTION(ss3c.str(), "")
{
if (i == z.size()-1){break;}
if (i == Ncomp-1){ break; }
double analytic = rHEOS.Reducing->d2Trdxi2__constxj(z, i, xN_flag);
double v1 = rHEOS_pluszi.Reducing->dTrdxi__constxj(rHEOS_pluszi.get_mole_fractions(), i, xN_flag);
double v2 = rHEOS_minuszi.Reducing->dTrdxi__constxj(rHEOS_minuszi.get_mole_fractions(), i, xN_flag);
@@ -721,23 +783,10 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]")
ss3d << "dalphar_dxi, i=" << i;
SECTION(ss3d.str(), "")
{
if (i == z.size()-1){break;}
if (i == Ncomp-1){ break; }
double analytic = MixtureDerivatives::dalphar_dxi(rHEOS, i, xN_flag);
shared_ptr<HelmholtzEOSMixtureBackend> plus(new HelmholtzEOSMixtureBackend(names));
plus->specify_phase(iphase_gas);
plus->set_mole_fractions(zp);
plus->calc_reducing_state();
SimpleState red = plus->get_reducing_state();
plus->update(DmolarT_INPUTS, red.rhomolar*rHEOS.delta(), red.T/rHEOS.tau());
double v1 = plus->alphar();
shared_ptr<HelmholtzEOSMixtureBackend> minus(new HelmholtzEOSMixtureBackend(names));
minus->specify_phase(iphase_gas);
minus->set_mole_fractions(zm);
minus->calc_reducing_state();
red = minus->get_reducing_state();
minus->update(DmolarT_INPUTS, red.rhomolar*rHEOS.delta(), red.T/rHEOS.tau());
double v2 = minus->alphar();
double v1 = rHEOS_pluszi_consttaudelta.alphar();
double v2 = rHEOS_minuszi_consttaudelta.alphar();
double numeric = (v1 - v2)/(2*dz);
double err = std::abs((numeric-analytic)/analytic);
CAPTURE(numeric);
@@ -746,35 +795,20 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]")
}
// These derivatives depend on both the i and j indices
for (std::size_t j = 0; j < z.size(); ++j){
if (xN_flag == XN_DEPENDENT && j == z.size()){ continue; }
shared_ptr<HelmholtzEOSMixtureBackend> HEOS, HEOS_pluszj, HEOS_minuszj;
HEOS_pluszj.reset(new HelmholtzEOSMixtureBackend(names));
HEOS_pluszj->specify_phase(iphase_gas);
HelmholtzEOSMixtureBackend &rHEOS_pluszj = *(HEOS_pluszj.get());
std::vector<CoolPropDbl> zp = z; /// Copy base composition
zp[j] += dz;
if (xN_flag == XN_DEPENDENT)
zp[z.size()-1] -= dz;
rHEOS_pluszj.set_mole_fractions(zp);
rHEOS_pluszj.update(DmolarT_INPUTS, rho1, T1);
HEOS_minuszj.reset(new HelmholtzEOSMixtureBackend(names));
HEOS_minuszj->specify_phase(iphase_gas);
HelmholtzEOSMixtureBackend &rHEOS_minuszj = *(HEOS_minuszj.get());
std::vector<CoolPropDbl> zm = z; /// Copy base composition
zm[j] -= dz;
if (xN_flag == XN_DEPENDENT)
zm[z.size()-1] += dz;
rHEOS_minuszj.set_mole_fractions(zm);
rHEOS_minuszj.update(DmolarT_INPUTS, rho1, T1);
for (std::size_t j = 0; j < Ncomp; ++j){
if (xN_flag == XN_DEPENDENT && j == Ncomp){ continue; }
HelmholtzEOSMixtureBackend & rHEOS_pluszj = (xN_flag == XN_INDEPENDENT) ? *(HEOS_plusz_xNindep[Ncomp][j].get()) : *(HEOS_plusz_xNdep[Ncomp][j].get());
HelmholtzEOSMixtureBackend & rHEOS_minuszj = (xN_flag == XN_INDEPENDENT) ? *(HEOS_minusz_xNindep[Ncomp][j].get()) : *(HEOS_minusz_xNdep[Ncomp][j].get());
HelmholtzEOSMixtureBackend &rHEOS_pluszj_consttaudelta = (xN_flag == XN_INDEPENDENT) ? *(HEOS_plusz_consttaudelta_xNindep[Ncomp][j].get()) : *(HEOS_plusz_consttaudelta_xNdep[Ncomp][j].get());
HelmholtzEOSMixtureBackend &rHEOS_minuszj_consttaudelta = (xN_flag == XN_INDEPENDENT) ? *(HEOS_minusz_consttaudelta_xNindep[Ncomp][j].get()) : *(HEOS_minusz_consttaudelta_xNdep[Ncomp][j].get());
std::ostringstream ss1a;
ss1a << "dln_fugacity_dxj__constT_rho_xi, i=" << i << ", j=" << j;
SECTION(ss1a.str(), "")
{
if (xN_flag == XN_INDEPENDENT){continue;}
if (j == z.size()-1){break;}
if (j == Ncomp-1){break;}
double analytic = MixtureDerivatives::dln_fugacity_dxj__constT_rho_xi(rHEOS, i, j, xN_flag);
double v1 = log(MixtureDerivatives::fugacity_i(rHEOS_pluszj, i, xN_flag));
double v2 = log(MixtureDerivatives::fugacity_i(rHEOS_minuszj, i, xN_flag));
@@ -788,7 +822,7 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]")
ss2 << "d_ndTrdni_dxj, i=" << i << ", j=" << j;
SECTION(ss2.str(), "")
{
if (j == z.size()-1){break;}
if (j == Ncomp-1){ break; }
double analytic = rHEOS.Reducing->d_ndTrdni_dxj__constxi(rHEOS.get_mole_fractions(), i, j, xN_flag);
double v1 = rHEOS_pluszj.Reducing->ndTrdni__constnj(rHEOS_pluszj.get_mole_fractions(), i, xN_flag);
double v2 = rHEOS_minuszj.Reducing->ndTrdni__constnj(rHEOS_minuszj.get_mole_fractions(), i, xN_flag);
@@ -802,7 +836,7 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]")
ss4 << "d_ndrhomolarrdni_dxj, i=" << i << ", j=" << j;
SECTION(ss4.str(), "")
{
if (j == z.size()-1){break;}
if (j == Ncomp-1){ break; }
double analytic = rHEOS.Reducing->d_ndrhorbardni_dxj__constxi(rHEOS.get_mole_fractions(), i, j, xN_flag);
double v1 = rHEOS_pluszj.Reducing->ndrhorbardni__constnj(rHEOS_pluszj.get_mole_fractions(), i, xN_flag);
double v2 = rHEOS_minuszj.Reducing->ndrhorbardni__constnj(rHEOS_minuszj.get_mole_fractions(), i, xN_flag);
@@ -816,7 +850,7 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]")
ss3 << "d_ndalphardni_dxj__constT_V_xi, i=" << i << ", j=" << j;
SECTION(ss3.str(), "")
{
if (j == z.size()-1){break;}
if (j == Ncomp-1){ break; }
double analytic = MixtureDerivatives::d_ndalphardni_dxj__constT_V_xi(rHEOS, i, j, xN_flag);
double v1 = MixtureDerivatives::ndalphar_dni__constT_V_nj(rHEOS_pluszj, i, xN_flag);
double v2 = MixtureDerivatives::ndalphar_dni__constT_V_nj(rHEOS_minuszj, i, xN_flag);
@@ -830,23 +864,10 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]")
ss3a << "d2alphardxidxj, i=" << i << ", j=" << j;
SECTION(ss3a.str(), "")
{
if (j == z.size()-1){break;}
if (j == Ncomp-1){ break; }
double analytic = MixtureDerivatives::d2alphardxidxj(rHEOS,i,j,xN_flag);
shared_ptr<HelmholtzEOSMixtureBackend> plus(new HelmholtzEOSMixtureBackend(names));
plus->specify_phase(iphase_gas);
plus->set_mole_fractions(zp);
plus->calc_reducing_state();
SimpleState red = plus->get_reducing_state();
plus->update(DmolarT_INPUTS, red.rhomolar*rHEOS.delta(), red.T/rHEOS.tau());
double v1 = MixtureDerivatives::dalphar_dxi(*(plus.get()), i, xN_flag);
shared_ptr<HelmholtzEOSMixtureBackend> minus(new HelmholtzEOSMixtureBackend(names));
minus->specify_phase(iphase_gas);
minus->set_mole_fractions(zm);
minus->calc_reducing_state();
red = minus->get_reducing_state();
minus->update(DmolarT_INPUTS, red.rhomolar*rHEOS.delta(), red.T/rHEOS.tau());
double v2 = MixtureDerivatives::dalphar_dxi(*(minus.get()), i, xN_flag);
double v1 = MixtureDerivatives::dalphar_dxi(rHEOS_pluszj_consttaudelta, i, xN_flag);
double v2 = MixtureDerivatives::dalphar_dxi(rHEOS_minuszj_consttaudelta, i, xN_flag);
double numeric = (v1 - v2)/(2*dz);
if (std::abs(numeric) < DBL_EPSILON && std::abs(analytic) < DBL_EPSILON){break;}
double err;
@@ -865,7 +886,7 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]")
ss3b << "d2Trdxidxj, i=" << i << ", j=" << j;
SECTION(ss3b.str(), "")
{
if (j == z.size()-1 || i == j){break;}
if (j == Ncomp-1 || i == j){ break; }
double analytic = rHEOS.Reducing->d2Trdxidxj(z, i, j, xN_flag);
double v1 = rHEOS.Reducing->dTrdxi__constxj(rHEOS_pluszj.get_mole_fractions(), i, xN_flag);
double v2 = rHEOS.Reducing->dTrdxi__constxj(rHEOS_minuszj.get_mole_fractions(), i, xN_flag);