From af8691d0c5efeaf47f6754f1fbc6ed484c03283f Mon Sep 17 00:00:00 2001 From: Ian Bell Date: Fri, 12 Sep 2014 20:59:41 +0200 Subject: [PATCH] Both xN independent and dependent are now tested in mixture derivative code. Signed-off-by: Ian Bell --- src/Backends/Helmholtz/MixtureDerivatives.cpp | 912 +++++++++--------- 1 file changed, 466 insertions(+), 446 deletions(-) diff --git a/src/Backends/Helmholtz/MixtureDerivatives.cpp b/src/Backends/Helmholtz/MixtureDerivatives.cpp index 1311ae34..7750f877 100644 --- a/src/Backends/Helmholtz/MixtureDerivatives.cpp +++ b/src/Backends/Helmholtz/MixtureDerivatives.cpp @@ -385,462 +385,482 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") ss00 << "Mixture with " << Ncomp << " components"; SECTION(ss00.str(),"") { - x_N_dependency_flag xN_flag = XN_INDEPENDENT; - - std::vector names; - std::vector z; - shared_ptr 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; - } - 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); - - // These ones only require the i index - for (std::size_t i = 0; i< z.size();++i) - { - shared_ptr 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 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 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); - - std::ostringstream ss0; - ss0 << "dln_fugacity_i_dT__constrho_n, i=" << i; - SECTION(ss0.str(),"") - { - double analytic = MixtureDerivatives::dln_fugacity_i_dT__constrho_n(rHEOS, i, xN_flag); - double v1 = log(MixtureDerivatives::fugacity_i(rHEOS_plusT_constrho, i, xN_flag)); - double v2 = log(MixtureDerivatives::fugacity_i(rHEOS_minusT_constrho, i, xN_flag)); - double numeric = (v1 - v2)/(2*dT); - double err = std::abs((numeric-analytic)/analytic); - CAPTURE(numeric); - CAPTURE(analytic); - CHECK(err < 1e-8); + x_N_dependency_flag xN_flag; + for (std::size_t xNxN = 0; xNxN <=1; xNxN++){ + std::ostringstream ss000; + std::string xN_string; + if (xNxN == 0){ + xN_flag = XN_DEPENDENT; + xN_string = "Mole fractions dependent"; } - std::ostringstream ss0a; - ss0a << "dln_fugacity_i_dp__constT_n, i=" << i; - SECTION(ss0a.str(),"") - { - double analytic = MixtureDerivatives::dln_fugacity_i_dp__constT_n(rHEOS, i, xN_flag); - double v1 = log(MixtureDerivatives::fugacity_i(rHEOS_plusrho_constT, i, xN_flag)), p1 = rHEOS_plusrho_constT.p(); - double v2 = log(MixtureDerivatives::fugacity_i(rHEOS_minusrho_constT, i, xN_flag)), p2 = rHEOS_minusrho_constT.p(); - double numeric = (v1 - v2)/(p1 - p2); - double err = std::abs((numeric-analytic)/analytic); - CAPTURE(numeric); - CAPTURE(analytic); - CHECK(err < 1e-8); + else{ + xN_flag = XN_INDEPENDENT; + xN_string = "Mole fractions independent"; } - std::ostringstream ss0b; - ss0b << "dln_fugacity_i_drho__constT_n, i=" << i; - SECTION(ss0b.str(), "") + + ss000 << xN_string; + SECTION(ss000.str(),"") { - double analytic = MixtureDerivatives::dln_fugacity_i_drho__constT_n(rHEOS, i, xN_flag); - double v1 = log(MixtureDerivatives::fugacity_i(rHEOS_plusrho_constT, i, xN_flag)); - double v2 = log(MixtureDerivatives::fugacity_i(rHEOS_minusrho_constT, i, xN_flag)); - double numeric = (v1 - v2)/(2*dT); - double err = std::abs((numeric-analytic)/analytic); - CAPTURE(numeric); - CAPTURE(analytic); - CHECK(err < 1e-8); - } - std::ostringstream ss2; - ss2 << "dln_fugacity_coefficient_dp__constT_n, i=" << i; - SECTION(ss2.str(), "") - { - double analytic = MixtureDerivatives::dln_fugacity_coefficient_dp__constT_n(rHEOS, i, xN_flag); - double v1 = MixtureDerivatives::ln_fugacity_coefficient(rHEOS_plusrho_constT, i, xN_flag), p1 = rHEOS_plusrho_constT.p(); - double v2 = MixtureDerivatives::ln_fugacity_coefficient(rHEOS_minusrho_constT, i, xN_flag), p2 = rHEOS_minusrho_constT.p(); - double numeric = (v1 - v2)/(p1 - p2); - double err = std::abs((numeric-analytic)/analytic); - CHECK(err < 1e-8); - } - std::ostringstream ss1; - ss1 << "dln_fugacity_coefficient_dT__constp_n, i=" << i; - SECTION(ss1.str(), "") - { - double T1 = 300, dT = 1e-3; - rHEOS.specify_phase(iphase_gas); + std::vector names; + std::vector z; + shared_ptr HEOS, HEOS_plusT_constrho, HEOS_plusrho_constT, HEOS_minusT_constrho, HEOS_minusrho_constT; + names.resize(Ncomp); + z.resize(Ncomp); - rHEOS.update(PT_INPUTS, 101325, T1); - double analytic = MixtureDerivatives::dln_fugacity_coefficient_dT__constp_n(rHEOS, i, xN_flag); - - rHEOS.update(PT_INPUTS, 101325, T1 + dT); - double v1 = MixtureDerivatives::ln_fugacity_coefficient(rHEOS, i, xN_flag); - rHEOS.update(PT_INPUTS, 101325, T1 - dT); - double v2 = MixtureDerivatives::ln_fugacity_coefficient(rHEOS, i, xN_flag); - - double numeric = (v1 - v2)/(2*dT); - double err = std::abs((numeric-analytic)/analytic); - CHECK(err < 1e-6); - } - std::ostringstream ss1a; - ss1a << "dln_fugacity_i_dT__constp_n, i=" << i; - SECTION(ss1a.str(), "") - { - double T1 = 300, dT = 1e-3; - rHEOS.specify_phase(iphase_gas); - - rHEOS.update(PT_INPUTS, 101325, T1); - double analytic = MixtureDerivatives::dln_fugacity_i_dT__constp_n(rHEOS, i, xN_flag); - - rHEOS.update(PT_INPUTS, 101325, T1 + dT); - double v1 = log(MixtureDerivatives::fugacity_i(rHEOS, i, xN_flag)); - rHEOS.update(PT_INPUTS, 101325, T1 - dT); - double v2 = log(MixtureDerivatives::fugacity_i(rHEOS, i, xN_flag)); - - double numeric = (v1 - v2)/(2*dT); - double err = std::abs((numeric-analytic)/analytic); - CHECK(err < 1e-6); - } - - std::ostringstream ss3; - ss3 << "d_ndalphardni_dDelta, i=" << i; - SECTION(ss3.str(), "") - { - double analytic = MixtureDerivatives::d_ndalphardni_dDelta(rHEOS, i, xN_flag); - double v1 = MixtureDerivatives::ndalphar_dni__constT_V_nj(rHEOS_plusrho_constT, i, xN_flag), delta1 = rHEOS_plusrho_constT.delta(); - double v2 = MixtureDerivatives::ndalphar_dni__constT_V_nj(rHEOS_minusrho_constT, i, xN_flag), delta2 = rHEOS_minusrho_constT.delta(); - double numeric = (v1 - v2)/(delta1 - delta2); - double err = std::abs((numeric-analytic)/analytic); - CHECK(err < 1e-8); - } - std::ostringstream ss3a; - ss3a << "d2alphar_dxi_dDelta, i=" << i; - SECTION(ss3a.str(), "") - { - if (i == z.size()-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(); - double numeric = (v1 - v2)/(delta1 - delta2); - double err = std::abs((numeric-analytic)/analytic); - CHECK(err < 1e-8); - } - std::ostringstream ss4a; - ss4a << "d2alphar_dxi_dTau, i=" << i; - SECTION(ss4a.str(), "") - { - if (i == z.size()-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(); - double numeric = (v1 - v2)/(tau1 - tau2); - double err = std::abs((numeric-analytic)/analytic); - CHECK(err < 1e-8); - } - std::ostringstream ss4; - ss4 << "d_ndalphardni_dTau, i=" << i; - SECTION(ss4.str(), "") - { - double analytic = MixtureDerivatives::d_ndalphardni_dTau(rHEOS, i, xN_flag); - double v1 = MixtureDerivatives::ndalphar_dni__constT_V_nj(rHEOS_plusT_constrho, i, xN_flag), tau1 = rHEOS_plusT_constrho.tau(); - double v2 = MixtureDerivatives::ndalphar_dni__constT_V_nj(rHEOS_minusT_constrho, i, xN_flag), tau2 = rHEOS_minusT_constrho.tau(); - double numeric = (v1 - v2)/(tau1 - tau2); - double err = std::abs((numeric-analytic)/analytic); - CHECK(err < 1e-8); - } - std::ostringstream ss5; - ss5 << "dpdxj__constT_V_xi, i=" << i; - SECTION(ss5.str(), "") - { - if (i==z.size()-1){break;} - double analytic = MixtureDerivatives::dpdxj__constT_V_xi(rHEOS, i, xN_flag); - double v1 = rHEOS_pluszi.p(); - double v2 = rHEOS_minuszi.p(); - double numeric = (v1 - v2)/(2*dz); - double err = std::abs((numeric-analytic)/analytic); - CAPTURE(numeric); - CAPTURE(analytic); - CHECK(err < 1e-8); - } - std::ostringstream ss5a; - ss5a << "dtaudxj__constT_V_xi, i=" << i; - SECTION(ss5a.str(), "") - { - if (i==z.size()-1){break;} - double analytic = MixtureDerivatives::dtau_dxj__constT_V_xi(rHEOS, i, xN_flag); - double v1 = rHEOS_pluszi.tau(); - double v2 = rHEOS_minuszi.tau(); - double numeric = (v1 - v2)/(2*dz); - double err = std::abs((numeric-analytic)/analytic); - CAPTURE(numeric); - CAPTURE(analytic); - CHECK(err < 1e-8); - } - std::ostringstream ss5b; - ss5b << "ddeltadxj__constT_V_xi, i=" << i; - SECTION(ss5b.str(), "") - { - if (i==z.size()-1){break;} - double analytic = MixtureDerivatives::ddelta_dxj__constT_V_xi(rHEOS, i, xN_flag); - double v1 = rHEOS_pluszi.delta(); - double v2 = rHEOS_minuszi.delta(); - double numeric = (v1 - v2)/(2*dz); - double err = std::abs((numeric-analytic)/analytic); - CAPTURE(numeric); - CAPTURE(analytic); - CHECK(err < 1e-8); - } - std::ostringstream ss6; - ss6 << "d_dalpharddelta_dxj__constT_V_xi, i=" << i; - SECTION(ss6.str(), "") - { - if (i==z.size()-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(); - double numeric = (v1 - v2)/(2*dz); - double err = std::abs((numeric-analytic)/analytic); - CAPTURE(numeric); - CAPTURE(analytic); - CHECK(err < 1e-8); - } - std::ostringstream ss7; - ss7 << "dTrdxi__constxj, i=" << i; - SECTION(ss7.str(), "") - { - if (i == z.size()-1){break;} - double analytic = rHEOS.Reducing->dTrdxi__constxj(rHEOS.get_const_mole_fractions(), i, xN_flag); - double v1 = rHEOS_pluszi.Reducing->Tr(rHEOS_pluszi.get_const_mole_fractions()); - double v2 = rHEOS_minuszi.Reducing->Tr(rHEOS_minuszi.get_const_mole_fractions()); - double numeric = (v1 - v2)/(2*dz); - double err = std::abs((numeric-analytic)/analytic); - CAPTURE(numeric); - CAPTURE(analytic); - CHECK(err < 1e-8); - } - std::ostringstream ss8; - ss8 << "drhormolardxi__constxj, i=" << i; - SECTION(ss8.str(), "") - { - if (i == z.size()-1){break;} - double analytic = rHEOS.Reducing->drhormolardxi__constxj(rHEOS.get_const_mole_fractions(), i, xN_flag); - double v1 = rHEOS_pluszi.Reducing->rhormolar(rHEOS_pluszi.get_const_mole_fractions()); - double v2 = rHEOS_minuszi.Reducing->rhormolar(rHEOS_minuszi.get_const_mole_fractions()); - double numeric = (v1 - v2)/(2*dz); - double err = std::abs((numeric-analytic)/analytic); - CAPTURE(numeric); - CAPTURE(analytic); - CHECK(err < 1e-8); - } - std::ostringstream ss3c; - ss3c << "d2Trdxi2__constxj, i=" << i; - SECTION(ss3c.str(), "") - { - if (i == z.size()-1){break;} - double analytic = rHEOS.Reducing->d2Trdxi2__constxj(z, i, xN_flag); - double v1 = rHEOS_pluszi.Reducing->dTrdxi__constxj(rHEOS_pluszi.get_const_mole_fractions(), i, xN_flag); - double v2 = rHEOS_minuszi.Reducing->dTrdxi__constxj(rHEOS_minuszi.get_const_mole_fractions(), i, xN_flag); - double numeric = (v1 - v2)/(2*dz); - double err = std::abs((numeric-analytic)/analytic); - CAPTURE(numeric); - CAPTURE(analytic); - CHECK(err < 1e-8); - } - std::ostringstream ss3d; - ss3d << "dalphar_dxi, i=" << i; - SECTION(ss3d.str(), "") - { - if (i == z.size()-1){break;} - double analytic = MixtureDerivatives::dalphar_dxi(rHEOS, i, xN_flag); - - shared_ptr 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 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 numeric = (v1 - v2)/(2*dz); - double err = std::abs((numeric-analytic)/analytic); - CAPTURE(numeric); - CAPTURE(analytic); - CHECK(err < 1e-8); - } - - // These derivatives depend on both the i and j indices - for (std::size_t j = 0; j < z.size(); ++j){ - - shared_ptr 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 zp = z; /// Copy base composition - zp[j] += dz; 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 zm = z; /// Copy base composition - zm[j] -= dz; zm[z.size()-1] += dz; - rHEOS_minuszj.set_mole_fractions(zm); - rHEOS_minuszj.update(DmolarT_INPUTS, rho1, T1); - - std::ostringstream ss1a; - ss1a << "dln_fugacity_dxj__constT_rho_xi, i=" << i << ", j=" << j; - SECTION(ss1a.str(), "") + /* Set up a test class for the mixture tests */ + if (Ncomp == 2) { - if (j == z.size()-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)); - double numeric = (v1 - v2)/(2*dz); - double err = std::abs((numeric-analytic)/analytic); - CAPTURE(numeric); - CAPTURE(analytic); - CHECK(err < 1e-8); + names[0] = "Ethane"; names[1] = "Propane"; + z[0] = 0.3; z[1] = 0.7; } - std::ostringstream ss2; - ss2 << "d_ndTrdni_dxj, i=" << i << ", j=" << j; - SECTION(ss2.str(), "") + else if (Ncomp == 3) { - if (j == z.size()-1){break;} - double analytic = rHEOS.Reducing->d_ndTrdni_dxj__constxi(rHEOS.get_const_mole_fractions(), i, j, xN_flag); - double v1 = rHEOS_pluszj.Reducing->ndTrdni__constnj(rHEOS_pluszj.get_const_mole_fractions(), i, xN_flag); - double v2 = rHEOS_minuszj.Reducing->ndTrdni__constnj(rHEOS_minuszj.get_const_mole_fractions(), i, xN_flag); - double numeric = (v1 - v2)/(2*dz); - double err = std::abs((numeric-analytic)/analytic); - CAPTURE(numeric); - CAPTURE(analytic); - CHECK(err < 1e-8); + names[0] = "Ethane"; names[1] = "Propane"; names[2] = "Methane"; + z[0] = 0.3; z[1] = 0.4; z[2] = 0.3; } - std::ostringstream ss4; - ss4 << "d_ndrhomolarrdni_dxj, i=" << i << ", j=" << j; - SECTION(ss4.str(), "") + 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); + + // These ones only require the i index + for (std::size_t i = 0; i< z.size();++i) { - if (j == z.size()-1){break;} - double analytic = rHEOS.Reducing->d_ndrhorbardni_dxj__constxi(rHEOS.get_const_mole_fractions(), i, j, xN_flag); - double v1 = rHEOS_pluszj.Reducing->ndrhorbardni__constnj(rHEOS_pluszj.get_const_mole_fractions(), i, xN_flag); - double v2 = rHEOS_minuszj.Reducing->ndrhorbardni__constnj(rHEOS_minuszj.get_const_mole_fractions(), i, xN_flag); - double numeric = (v1 - v2)/(2*dz); - double err = std::abs((numeric-analytic)/analytic); - CAPTURE(numeric); - CAPTURE(analytic); - CHECK(err < 1e-8); - } - std::ostringstream ss3; - ss3 << "d_ndalphardni_dxj__constT_V_xi, i=" << i << ", j=" << j; - SECTION(ss3.str(), "") - { - if (j == z.size()-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); - double numeric = (v1 - v2)/(2*dz); - double err = std::abs((numeric-analytic)/analytic); - CAPTURE(numeric); - CAPTURE(analytic); - CHECK(err < 1e-8); - } - std::ostringstream ss3a; - ss3a << "d2alphardxidxj, i=" << i << ", j=" << j; - SECTION(ss3a.str(), "") - { - if (j == z.size()-1){break;} - double analytic = MixtureDerivatives::d2alphardxidxj(rHEOS,i,j,xN_flag); + shared_ptr 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 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 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); - shared_ptr 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 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 numeric = (v1 - v2)/(2*dz); - double err = std::abs((numeric-analytic)/analytic); - if (std::abs(numeric) < DBL_EPSILON && std::abs(analytic) < DBL_EPSILON){break;} - CAPTURE(numeric); - CAPTURE(analytic); - CHECK(err < 1e-8); - } - std::ostringstream ss3b; - ss3b << "d2Trdxidxj, i=" << i << ", j=" << j; - SECTION(ss3b.str(), "") - { - if (j == z.size()-1 || i == j){break;} - double analytic = rHEOS.Reducing->d2Trdxidxj(z, i, j, xN_flag); - double v1 = rHEOS.Reducing->dTrdxi__constxj(rHEOS_pluszj.get_const_mole_fractions(), i, xN_flag); - double v2 = rHEOS.Reducing->dTrdxi__constxj(rHEOS_minuszj.get_const_mole_fractions(), i, xN_flag); - double numeric = (v1 - v2)/(2*dz); - double err = std::abs((numeric-analytic)/analytic); - if (std::abs(numeric) < DBL_EPSILON && std::abs(analytic) < DBL_EPSILON){break;} - CAPTURE(numeric); - CAPTURE(analytic); - CHECK(err < 1e-8); + std::ostringstream ss0; + ss0 << "dln_fugacity_i_dT__constrho_n, i=" << i; + SECTION(ss0.str(),"") + { + double analytic = MixtureDerivatives::dln_fugacity_i_dT__constrho_n(rHEOS, i, xN_flag); + double v1 = log(MixtureDerivatives::fugacity_i(rHEOS_plusT_constrho, i, xN_flag)); + double v2 = log(MixtureDerivatives::fugacity_i(rHEOS_minusT_constrho, i, xN_flag)); + double numeric = (v1 - v2)/(2*dT); + double err = std::abs((numeric-analytic)/analytic); + CAPTURE(numeric); + CAPTURE(analytic); + CHECK(err < 1e-8); + } + std::ostringstream ss0a; + ss0a << "dln_fugacity_i_dp__constT_n, i=" << i; + SECTION(ss0a.str(),"") + { + double analytic = MixtureDerivatives::dln_fugacity_i_dp__constT_n(rHEOS, i, xN_flag); + double v1 = log(MixtureDerivatives::fugacity_i(rHEOS_plusrho_constT, i, xN_flag)), p1 = rHEOS_plusrho_constT.p(); + double v2 = log(MixtureDerivatives::fugacity_i(rHEOS_minusrho_constT, i, xN_flag)), p2 = rHEOS_minusrho_constT.p(); + double numeric = (v1 - v2)/(p1 - p2); + double err = std::abs((numeric-analytic)/analytic); + CAPTURE(numeric); + CAPTURE(analytic); + CHECK(err < 1e-8); + } + std::ostringstream ss0b; + ss0b << "dln_fugacity_i_drho__constT_n, i=" << i; + SECTION(ss0b.str(), "") + { + double analytic = MixtureDerivatives::dln_fugacity_i_drho__constT_n(rHEOS, i, xN_flag); + double v1 = log(MixtureDerivatives::fugacity_i(rHEOS_plusrho_constT, i, xN_flag)); + double v2 = log(MixtureDerivatives::fugacity_i(rHEOS_minusrho_constT, i, xN_flag)); + double numeric = (v1 - v2)/(2*dT); + double err = std::abs((numeric-analytic)/analytic); + CAPTURE(numeric); + CAPTURE(analytic); + CHECK(err < 1e-8); + } + std::ostringstream ss2; + ss2 << "dln_fugacity_coefficient_dp__constT_n, i=" << i; + SECTION(ss2.str(), "") + { + double analytic = MixtureDerivatives::dln_fugacity_coefficient_dp__constT_n(rHEOS, i, xN_flag); + double v1 = MixtureDerivatives::ln_fugacity_coefficient(rHEOS_plusrho_constT, i, xN_flag), p1 = rHEOS_plusrho_constT.p(); + double v2 = MixtureDerivatives::ln_fugacity_coefficient(rHEOS_minusrho_constT, i, xN_flag), p2 = rHEOS_minusrho_constT.p(); + double numeric = (v1 - v2)/(p1 - p2); + double err = std::abs((numeric-analytic)/analytic); + CHECK(err < 1e-8); + } + std::ostringstream ss1; + ss1 << "dln_fugacity_coefficient_dT__constp_n, i=" << i; + SECTION(ss1.str(), "") + { + double T1 = 300, dT = 1e-3; + rHEOS.specify_phase(iphase_gas); + + rHEOS.update(PT_INPUTS, 101325, T1); + double analytic = MixtureDerivatives::dln_fugacity_coefficient_dT__constp_n(rHEOS, i, xN_flag); + + rHEOS.update(PT_INPUTS, 101325, T1 + dT); + double v1 = MixtureDerivatives::ln_fugacity_coefficient(rHEOS, i, xN_flag); + rHEOS.update(PT_INPUTS, 101325, T1 - dT); + double v2 = MixtureDerivatives::ln_fugacity_coefficient(rHEOS, i, xN_flag); + + double numeric = (v1 - v2)/(2*dT); + double err = std::abs((numeric-analytic)/analytic); + CHECK(err < 1e-6); + } + std::ostringstream ss1a; + ss1a << "dln_fugacity_i_dT__constp_n, i=" << i; + SECTION(ss1a.str(), "") + { + double T1 = 300, dT = 1e-3; + rHEOS.specify_phase(iphase_gas); + + rHEOS.update(PT_INPUTS, 101325, T1); + double analytic = MixtureDerivatives::dln_fugacity_i_dT__constp_n(rHEOS, i, xN_flag); + + rHEOS.update(PT_INPUTS, 101325, T1 + dT); + double v1 = log(MixtureDerivatives::fugacity_i(rHEOS, i, xN_flag)); + rHEOS.update(PT_INPUTS, 101325, T1 - dT); + double v2 = log(MixtureDerivatives::fugacity_i(rHEOS, i, xN_flag)); + + double numeric = (v1 - v2)/(2*dT); + double err = std::abs((numeric-analytic)/analytic); + CHECK(err < 1e-6); + } + + std::ostringstream ss3; + ss3 << "d_ndalphardni_dDelta, i=" << i; + SECTION(ss3.str(), "") + { + double analytic = MixtureDerivatives::d_ndalphardni_dDelta(rHEOS, i, xN_flag); + double v1 = MixtureDerivatives::ndalphar_dni__constT_V_nj(rHEOS_plusrho_constT, i, xN_flag), delta1 = rHEOS_plusrho_constT.delta(); + double v2 = MixtureDerivatives::ndalphar_dni__constT_V_nj(rHEOS_minusrho_constT, i, xN_flag), delta2 = rHEOS_minusrho_constT.delta(); + double numeric = (v1 - v2)/(delta1 - delta2); + double err = std::abs((numeric-analytic)/analytic); + CHECK(err < 1e-8); + } + std::ostringstream ss3a; + ss3a << "d2alphar_dxi_dDelta, i=" << i; + SECTION(ss3a.str(), "") + { + if (i == z.size()-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(); + double numeric = (v1 - v2)/(delta1 - delta2); + double err = std::abs((numeric-analytic)/analytic); + CHECK(err < 1e-8); + } + std::ostringstream ss4a; + ss4a << "d2alphar_dxi_dTau, i=" << i; + SECTION(ss4a.str(), "") + { + if (i == z.size()-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(); + double numeric = (v1 - v2)/(tau1 - tau2); + double err = std::abs((numeric-analytic)/analytic); + CHECK(err < 1e-8); + } + std::ostringstream ss4; + ss4 << "d_ndalphardni_dTau, i=" << i; + SECTION(ss4.str(), "") + { + double analytic = MixtureDerivatives::d_ndalphardni_dTau(rHEOS, i, xN_flag); + double v1 = MixtureDerivatives::ndalphar_dni__constT_V_nj(rHEOS_plusT_constrho, i, xN_flag), tau1 = rHEOS_plusT_constrho.tau(); + double v2 = MixtureDerivatives::ndalphar_dni__constT_V_nj(rHEOS_minusT_constrho, i, xN_flag), tau2 = rHEOS_minusT_constrho.tau(); + double numeric = (v1 - v2)/(tau1 - tau2); + double err = std::abs((numeric-analytic)/analytic); + CHECK(err < 1e-8); + } + std::ostringstream ss5; + ss5 << "dpdxj__constT_V_xi, i=" << i; + SECTION(ss5.str(), "") + { + if (i==z.size()-1){break;} + double analytic = MixtureDerivatives::dpdxj__constT_V_xi(rHEOS, i, xN_flag); + double v1 = rHEOS_pluszi.p(); + double v2 = rHEOS_minuszi.p(); + double numeric = (v1 - v2)/(2*dz); + double err = std::abs((numeric-analytic)/analytic); + CAPTURE(numeric); + CAPTURE(analytic); + CHECK(err < 1e-8); + } + std::ostringstream ss5a; + ss5a << "dtaudxj__constT_V_xi, i=" << i; + SECTION(ss5a.str(), "") + { + if (i==z.size()-1){break;} + double analytic = MixtureDerivatives::dtau_dxj__constT_V_xi(rHEOS, i, xN_flag); + double v1 = rHEOS_pluszi.tau(); + double v2 = rHEOS_minuszi.tau(); + double numeric = (v1 - v2)/(2*dz); + double err = std::abs((numeric-analytic)/analytic); + CAPTURE(numeric); + CAPTURE(analytic); + CHECK(err < 1e-8); + } + std::ostringstream ss5b; + ss5b << "ddeltadxj__constT_V_xi, i=" << i; + SECTION(ss5b.str(), "") + { + if (i==z.size()-1){break;} + double analytic = MixtureDerivatives::ddelta_dxj__constT_V_xi(rHEOS, i, xN_flag); + double v1 = rHEOS_pluszi.delta(); + double v2 = rHEOS_minuszi.delta(); + double numeric = (v1 - v2)/(2*dz); + double err = std::abs((numeric-analytic)/analytic); + CAPTURE(numeric); + CAPTURE(analytic); + CHECK(err < 1e-8); + } + std::ostringstream ss6; + ss6 << "d_dalpharddelta_dxj__constT_V_xi, i=" << i; + SECTION(ss6.str(), "") + { + if (i==z.size()-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(); + double numeric = (v1 - v2)/(2*dz); + double err = std::abs((numeric-analytic)/analytic); + CAPTURE(numeric); + CAPTURE(analytic); + CHECK(err < 1e-8); + } + std::ostringstream ss7; + ss7 << "dTrdxi__constxj, i=" << i; + SECTION(ss7.str(), "") + { + if (i == z.size()-1){break;} + double analytic = rHEOS.Reducing->dTrdxi__constxj(rHEOS.get_const_mole_fractions(), i, xN_flag); + double v1 = rHEOS_pluszi.Reducing->Tr(rHEOS_pluszi.get_const_mole_fractions()); + double v2 = rHEOS_minuszi.Reducing->Tr(rHEOS_minuszi.get_const_mole_fractions()); + double numeric = (v1 - v2)/(2*dz); + double err = std::abs((numeric-analytic)/analytic); + CAPTURE(numeric); + CAPTURE(analytic); + CHECK(err < 1e-8); + } + std::ostringstream ss8; + ss8 << "drhormolardxi__constxj, i=" << i; + SECTION(ss8.str(), "") + { + if (i == z.size()-1){break;} + double analytic = rHEOS.Reducing->drhormolardxi__constxj(rHEOS.get_const_mole_fractions(), i, xN_flag); + double v1 = rHEOS_pluszi.Reducing->rhormolar(rHEOS_pluszi.get_const_mole_fractions()); + double v2 = rHEOS_minuszi.Reducing->rhormolar(rHEOS_minuszi.get_const_mole_fractions()); + double numeric = (v1 - v2)/(2*dz); + double err = std::abs((numeric-analytic)/analytic); + CAPTURE(numeric); + CAPTURE(analytic); + CHECK(err < 1e-8); + } + std::ostringstream ss3c; + ss3c << "d2Trdxi2__constxj, i=" << i; + SECTION(ss3c.str(), "") + { + if (i == z.size()-1){break;} + double analytic = rHEOS.Reducing->d2Trdxi2__constxj(z, i, xN_flag); + double v1 = rHEOS_pluszi.Reducing->dTrdxi__constxj(rHEOS_pluszi.get_const_mole_fractions(), i, xN_flag); + double v2 = rHEOS_minuszi.Reducing->dTrdxi__constxj(rHEOS_minuszi.get_const_mole_fractions(), i, xN_flag); + double numeric = (v1 - v2)/(2*dz); + double err = std::abs((numeric-analytic)/analytic); + CAPTURE(numeric); + CAPTURE(analytic); + CHECK(err < 1e-8); + } + std::ostringstream ss3d; + ss3d << "dalphar_dxi, i=" << i; + SECTION(ss3d.str(), "") + { + if (i == z.size()-1){break;} + double analytic = MixtureDerivatives::dalphar_dxi(rHEOS, i, xN_flag); + + shared_ptr 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 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 numeric = (v1 - v2)/(2*dz); + double err = std::abs((numeric-analytic)/analytic); + CAPTURE(numeric); + CAPTURE(analytic); + CHECK(err < 1e-8); + } + + // These derivatives depend on both the i and j indices + for (std::size_t j = 0; j < z.size(); ++j){ + + shared_ptr 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 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 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); + + std::ostringstream ss1a; + ss1a << "dln_fugacity_dxj__constT_rho_xi, i=" << i << ", j=" << j; + SECTION(ss1a.str(), "") + { + if (j == z.size()-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)); + double numeric = (v1 - v2)/(2*dz); + double err = std::abs((numeric-analytic)/analytic); + CAPTURE(numeric); + CAPTURE(analytic); + CHECK(err < 1e-8); + } + std::ostringstream ss2; + ss2 << "d_ndTrdni_dxj, i=" << i << ", j=" << j; + SECTION(ss2.str(), "") + { + if (j == z.size()-1){break;} + double analytic = rHEOS.Reducing->d_ndTrdni_dxj__constxi(rHEOS.get_const_mole_fractions(), i, j, xN_flag); + double v1 = rHEOS_pluszj.Reducing->ndTrdni__constnj(rHEOS_pluszj.get_const_mole_fractions(), i, xN_flag); + double v2 = rHEOS_minuszj.Reducing->ndTrdni__constnj(rHEOS_minuszj.get_const_mole_fractions(), i, xN_flag); + double numeric = (v1 - v2)/(2*dz); + double err = std::abs((numeric-analytic)/analytic); + CAPTURE(numeric); + CAPTURE(analytic); + CHECK(err < 1e-8); + } + std::ostringstream ss4; + ss4 << "d_ndrhomolarrdni_dxj, i=" << i << ", j=" << j; + SECTION(ss4.str(), "") + { + if (j == z.size()-1){break;} + double analytic = rHEOS.Reducing->d_ndrhorbardni_dxj__constxi(rHEOS.get_const_mole_fractions(), i, j, xN_flag); + double v1 = rHEOS_pluszj.Reducing->ndrhorbardni__constnj(rHEOS_pluszj.get_const_mole_fractions(), i, xN_flag); + double v2 = rHEOS_minuszj.Reducing->ndrhorbardni__constnj(rHEOS_minuszj.get_const_mole_fractions(), i, xN_flag); + double numeric = (v1 - v2)/(2*dz); + double err = std::abs((numeric-analytic)/analytic); + CAPTURE(numeric); + CAPTURE(analytic); + CHECK(err < 1e-8); + } + std::ostringstream ss3; + ss3 << "d_ndalphardni_dxj__constT_V_xi, i=" << i << ", j=" << j; + SECTION(ss3.str(), "") + { + if (j == z.size()-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); + double numeric = (v1 - v2)/(2*dz); + double err = std::abs((numeric-analytic)/analytic); + CAPTURE(numeric); + CAPTURE(analytic); + CHECK(err < 1e-8); + } + std::ostringstream ss3a; + ss3a << "d2alphardxidxj, i=" << i << ", j=" << j; + SECTION(ss3a.str(), "") + { + if (j == z.size()-1){break;} + double analytic = MixtureDerivatives::d2alphardxidxj(rHEOS,i,j,xN_flag); + + shared_ptr 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 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 numeric = (v1 - v2)/(2*dz); + double err = std::abs((numeric-analytic)/analytic); + if (std::abs(numeric) < DBL_EPSILON && std::abs(analytic) < DBL_EPSILON){break;} + CAPTURE(numeric); + CAPTURE(analytic); + CHECK(err < 1e-8); + } + std::ostringstream ss3b; + ss3b << "d2Trdxidxj, i=" << i << ", j=" << j; + SECTION(ss3b.str(), "") + { + if (j == z.size()-1 || i == j){break;} + double analytic = rHEOS.Reducing->d2Trdxidxj(z, i, j, xN_flag); + double v1 = rHEOS.Reducing->dTrdxi__constxj(rHEOS_pluszj.get_const_mole_fractions(), i, xN_flag); + double v2 = rHEOS.Reducing->dTrdxi__constxj(rHEOS_minuszj.get_const_mole_fractions(), i, xN_flag); + double numeric = (v1 - v2)/(2*dz); + double err = std::abs((numeric-analytic)/analytic); + if (std::abs(numeric) < DBL_EPSILON && std::abs(analytic) < DBL_EPSILON){break;} + CAPTURE(numeric); + CAPTURE(analytic); + CHECK(err < 1e-8); + } + } } } }