From ad643aff7e4e20729d37851f2272e29c676b07fd Mon Sep 17 00:00:00 2001 From: Ian Bell Date: Tue, 11 Aug 2015 23:58:45 -0600 Subject: [PATCH 1/8] Add the rest of the mixture derivatives for the critical point calculations - implemented but something is not quite right. --- .../Helmholtz/HelmholtzEOSMixtureBackend.cpp | 19 ++- src/Backends/Helmholtz/MixtureDerivatives.cpp | 158 ++++++++++++++++++ src/Backends/Helmholtz/MixtureDerivatives.h | 92 ++++++++-- 3 files changed, 252 insertions(+), 17 deletions(-) diff --git a/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp b/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp index 1c72e7f6..b752443d 100644 --- a/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp +++ b/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp @@ -2918,16 +2918,22 @@ void HelmholtzEOSMixtureBackend::calc_critical_point(double rho0, double T0) double epsilon; std::size_t N = x.size(); std::vector r, xp; - std::vector > J(N, std::vector(N, 0)); + std::vector > J(N, std::vector(N, 0)), JJ(N, std::vector(N, 0)); Eigen::MatrixXd J0(N, N), adjL = adjugate(MixtureDerivatives::Lstar(HEOS, XN_INDEPENDENT)), + adjM = adjugate(MixtureDerivatives::Mstar(HEOS, XN_INDEPENDENT)), dLdTau = MixtureDerivatives::dLstar_dX(HEOS, XN_INDEPENDENT, iTau), dLdDelta = MixtureDerivatives::dLstar_dX(HEOS, XN_INDEPENDENT, iDelta), - adjM = adjugate(MixtureDerivatives::Mstar(HEOS, XN_INDEPENDENT)), - dMdTau = dLdTau, dMdDelta = dLdDelta; + dMdTau = MixtureDerivatives::dMstar_dX(HEOS, XN_INDEPENDENT, iTau), + dMdDelta = MixtureDerivatives::dMstar_dX(HEOS, XN_INDEPENDENT, iDelta); + + JJ[0][0] = (adjL*dLdTau).trace(); + JJ[0][1] = (adjL*dLdDelta).trace(); + JJ[1][0] = (adjM*dMdTau).trace(); + JJ[1][1] = (adjM*dMdDelta).trace(); - J0(0,0) = (adjL*dLdTau).trace(); - J0(0,1) = (adjL*dLdDelta).trace(); //std::cout << J0 << std::endl; + std::cout << JJ[0][0] << " " << JJ[0][1] << std::endl; + std::cout << JJ[1][0] << " " << JJ[1][1] << std::endl; std::vector r0 = call(x); // Build the Jacobian by column for (std::size_t i = 0; i < N; ++i) @@ -2942,7 +2948,8 @@ void HelmholtzEOSMixtureBackend::calc_critical_point(double rho0, double T0) J[j][i] = (r[j]-r0[j])/epsilon; } } - + std::cout << J[0][0] << " " << J[0][1] << std::endl; + std::cout << J[1][0] << " " << J[1][1] << std::endl; return J; }; }; diff --git a/src/Backends/Helmholtz/MixtureDerivatives.cpp b/src/Backends/Helmholtz/MixtureDerivatives.cpp index 0fe4fb3a..4db8220c 100644 --- a/src/Backends/Helmholtz/MixtureDerivatives.cpp +++ b/src/Backends/Helmholtz/MixtureDerivatives.cpp @@ -403,6 +403,13 @@ CoolPropDbl MixtureDerivatives::d_ndln_fugacity_i_dnj_dtau__constdelta_x(Helmhol CoolPropDbl MixtureDerivatives::d2_ndln_fugacity_i_dnj_dtau2__constdelta_x(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag){ return d2_ndalphardni_dTau2(HEOS, j, xN_flag) + d2_nd_ndalphardni_dnj_dTau2__constdelta_x(HEOS, i, j, xN_flag); } +CoolPropDbl MixtureDerivatives::d2_ndln_fugacity_i_dnj_ddelta2__consttau_x(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag){ + return d2_ndalphardni_dDelta2(HEOS, j, xN_flag) + d2_nd_ndalphardni_dnj_dDelta2__consttau_x(HEOS, i, j, xN_flag); +} +CoolPropDbl MixtureDerivatives::d2_ndln_fugacity_i_dnj_ddelta_dtau__constx(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag){ + return d2_ndalphardni_dDelta_dTau(HEOS, j, xN_flag) + d2_nd_ndalphardni_dnj_dDelta_dTau__constx(HEOS, i, j, xN_flag); +} + CoolPropDbl MixtureDerivatives::d_ndln_fugacity_i_dnj_ddelta__consttau_x(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag){ return d_ndalphardni_dDelta(HEOS, j, xN_flag) + d_nd_ndalphardni_dnj_dDelta__consttau_x(HEOS, i, j, xN_flag); } @@ -410,6 +417,12 @@ CoolPropDbl MixtureDerivatives::d_ndln_fugacity_i_dnj_ddxk__consttau_delta(Helmh CoolPropDbl s = -Kronecker_delta(i, j)*Kronecker_delta(i, k)/pow(HEOS.mole_fractions[i], 2); return s + d_ndalphardni_dxj__constdelta_tau_xi(HEOS, j, k, xN_flag) + d_nd_ndalphardni_dnj_dxk__consttau_delta(HEOS, i, j, k, xN_flag); } +CoolPropDbl MixtureDerivatives::d2_ndln_fugacity_i_dnj_dxk_dTau__constdelta(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, std::size_t k, x_N_dependency_flag xN_flag){ + return d2_ndalphardni_dxj_dTau__constdelta_xi(HEOS, j, k, xN_flag) + d2_nd_ndalphardni_dnj_dxk_dTau__constdelta(HEOS, i, j, k, xN_flag); +} +CoolPropDbl MixtureDerivatives::d2_ndln_fugacity_i_dnj_dxk_dDelta__consttau(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, std::size_t k, x_N_dependency_flag xN_flag){ + return d2_ndalphardni_dxj_dDelta__consttau_xi(HEOS, j, k, xN_flag) + d2_nd_ndalphardni_dnj_dxk_dDelta__consttau(HEOS, i, j, k, xN_flag); +} CoolPropDbl MixtureDerivatives::nd_ndln_fugacity_i_dnj_dnk__constT_V_xi(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, std::size_t k, x_N_dependency_flag xN_flag){ double sum = d_ndln_fugacity_i_dnj_dtau__constdelta_x(HEOS, i, j, xN_flag)*ndtaudni__constT_V_nj(HEOS, k, xN_flag) + d_ndln_fugacity_i_dnj_ddelta__consttau_x(HEOS, i, j, xN_flag)*nddeltadni__constT_V_nj(HEOS, k, xN_flag) @@ -534,6 +547,10 @@ CoolPropDbl MixtureDerivatives::d_nddeltadni_dxj__constdelta_tau(HelmholtzEOSMix return -HEOS.delta()/rhor*(HEOS.Reducing->d_ndrhorbardni_dxj__constxi(HEOS.mole_fractions, i, j, xN_flag) -1/rhor*HEOS.Reducing->drhormolardxi__constxj(HEOS.mole_fractions, j, xN_flag)*HEOS.Reducing->ndrhorbardni__constnj(HEOS.mole_fractions, i, xN_flag)); } +CoolPropDbl MixtureDerivatives::d2_nddeltadni_dxj_dDelta__consttau(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag) +{ + return d_nddeltadni_dxj__constdelta_tau(HEOS, i, j, xN_flag)/HEOS.delta(); +} CoolPropDbl MixtureDerivatives::ndtaudni__constT_V_nj(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, x_N_dependency_flag xN_flag) { @@ -550,6 +567,10 @@ CoolPropDbl MixtureDerivatives::d_ndtaudni_dxj__constdelta_tau(HelmholtzEOSMixtu return HEOS.tau()/Tr*(HEOS.Reducing->d_ndTrdni_dxj__constxi(HEOS.mole_fractions, i, j, xN_flag) -1/Tr*HEOS.Reducing->dTrdxi__constxj(HEOS.mole_fractions, j, xN_flag)*HEOS.Reducing->ndTrdni__constnj(HEOS.mole_fractions, i, xN_flag)); } +CoolPropDbl MixtureDerivatives::d2_ndtaudni_dxj_dTau__constdelta(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag) +{ + return d_ndtaudni_dxj__constdelta_tau(HEOS, i, j, xN_flag)/HEOS.tau(); +} CoolPropDbl MixtureDerivatives::d_ndalphardni_dxj__constdelta_tau_xi(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag) { double line1 = HEOS._delta.pt()*d2alphar_dxi_dDelta(HEOS, j, xN_flag)*(1-1/HEOS._reducing.rhomolar*HEOS.Reducing->ndrhorbardni__constnj(HEOS.mole_fractions, i, xN_flag)); @@ -617,6 +638,37 @@ CoolPropDbl MixtureDerivatives::d2_nd_ndalphardni_dnj_dTau2__constdelta_x(Helmho double line4 = d3_ndalphardni_dxj_dTau2__constdelta_xi(HEOS, i, j, xN_flag)-summer; return line1 + line2 + line3 + line4; } +CoolPropDbl MixtureDerivatives::d2_nd_ndalphardni_dnj_dDelta2__consttau_x(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag) +{ + double line1 = d3_ndalphardni_dDelta3(HEOS, i, xN_flag)*nddeltadni__constT_V_nj(HEOS, j, xN_flag); + double line2 = 2*d2_ndalphardni_dDelta2(HEOS, i, xN_flag)*d_nddeltadni_dDelta(HEOS, j, xN_flag); + double line3 = d3_ndalphardni_dDelta2_dTau(HEOS, i, xN_flag)*ndtaudni__constT_V_nj(HEOS, j, xN_flag); + double summer = 0; + std::size_t kmax = HEOS.mole_fractions.size(); + if (xN_flag == XN_DEPENDENT){ kmax--; } + for (unsigned int k = 0; k < kmax; k++) + { + summer += HEOS.mole_fractions[k]*d3_ndalphardni_dxj_dDelta2__consttau_xi(HEOS, i, k, xN_flag); + } + double line4 = d3_ndalphardni_dxj_dDelta2__consttau_xi(HEOS, i, j, xN_flag)-summer; + return line1 + line2 + line3 + line4; +} +CoolPropDbl MixtureDerivatives::d2_nd_ndalphardni_dnj_dDelta_dTau__constx(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag) +{ + double line1 = d3_ndalphardni_dDelta2_dTau(HEOS, i, xN_flag)*nddeltadni__constT_V_nj(HEOS, j, xN_flag); + double line2 = d2_ndalphardni_dDelta_dTau(HEOS, i, xN_flag)*d_nddeltadni_dDelta(HEOS, j, xN_flag); + double line3 = d2_ndalphardni_dDelta_dTau(HEOS, i, xN_flag)*d_ndtaudni_dTau(HEOS, j, xN_flag); + double line4 = d3_ndalphardni_dDelta2_dTau(HEOS, i, xN_flag)*ndtaudni__constT_V_nj(HEOS, j, xN_flag); + double summer = 0; + std::size_t kmax = HEOS.mole_fractions.size(); + if (xN_flag == XN_DEPENDENT){ kmax--; } + for (unsigned int k = 0; k < kmax; k++) + { + summer += HEOS.mole_fractions[k]*d3_ndalphardni_dxj_dDelta_dTau__constxi(HEOS, i, k, xN_flag); + } + double line5 = d3_ndalphardni_dxj_dDelta_dTau__constxi(HEOS, i, j, xN_flag)-summer; + return line1 + line2 + line3 + line4 + line5; +} CoolPropDbl MixtureDerivatives::d_nd_ndalphardni_dnj_dDelta__consttau_x(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag) { double line1 = d2_ndalphardni_dDelta2(HEOS, i, xN_flag)*nddeltadni__constT_V_nj(HEOS, j, xN_flag) @@ -651,6 +703,45 @@ CoolPropDbl MixtureDerivatives::d_nd_ndalphardni_dnj_dxk__consttau_delta(Helmhol } return line1 + line2 + line3; } +CoolPropDbl MixtureDerivatives::d2_nd_ndalphardni_dnj_dxk_dTau__constdelta(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, std::size_t k, x_N_dependency_flag xN_flag) +{ + double line1 = d2_ndalphardni_dDelta_dTau(HEOS, i, xN_flag)*d_nddeltadni_dxj__constdelta_tau(HEOS, j, k, xN_flag) + + d3_ndalphardni_dxj_dDelta_dTau__constxi(HEOS, i, k, xN_flag)*nddeltadni__constT_V_nj(HEOS, j, xN_flag); + + double line2 = d_ndalphardni_dTau(HEOS, i, xN_flag)*d2_ndtaudni_dxj_dTau__constdelta(HEOS, j, k, xN_flag) + + d2_ndalphardni_dTau2(HEOS, i, xN_flag)*d_ndtaudni_dxj__constdelta_tau(HEOS, j, k, xN_flag) + + d2_ndalphardni_dxj_dTau__constdelta_xi(HEOS, i, k, xN_flag)*d_ndtaudni_dTau(HEOS, j, xN_flag) + + d3_ndalphardni_dxj_dTau2__constdelta_xi(HEOS, i, k, xN_flag)*ndtaudni__constT_V_nj(HEOS, j, xN_flag); + + double line3 = d3_ndalphardni_dxj_dxk_dTau__constdelta_xi(HEOS, i, j, k, xN_flag) ; + + std::size_t mmax = HEOS.mole_fractions.size(); + if (xN_flag == XN_DEPENDENT){ mmax--; } + for (unsigned int m = 0; m < mmax; m++) + { + line3 -= HEOS.mole_fractions[m]*d3_ndalphardni_dxj_dxk_dTau__constdelta_xi(HEOS, i, m, k, xN_flag); + } + return line1 + line2 + line3; +} +CoolPropDbl MixtureDerivatives::d2_nd_ndalphardni_dnj_dxk_dDelta__consttau(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, std::size_t k, x_N_dependency_flag xN_flag) +{ + double line1 = d_ndalphardni_dDelta(HEOS, i, xN_flag)*d2_nddeltadni_dxj_dDelta__consttau(HEOS, j, k, xN_flag) + + d2_ndalphardni_dDelta2(HEOS, i, xN_flag)*d_nddeltadni_dxj__constdelta_tau(HEOS, j, k, xN_flag) + + d2_ndalphardni_dxj_dDelta__consttau_xi(HEOS, i, k, xN_flag)*d_nddeltadni_dDelta(HEOS, j, xN_flag) + + d3_ndalphardni_dxj_dDelta2__consttau_xi(HEOS, i, k, xN_flag)*nddeltadni__constT_V_nj(HEOS, j, xN_flag); + + double line2 = d2_ndalphardni_dDelta_dTau(HEOS, i, xN_flag)*d_ndtaudni_dxj__constdelta_tau(HEOS, j, k, xN_flag) + + d3_ndalphardni_dxj_dDelta_dTau__constxi(HEOS, i, k, xN_flag)*ndtaudni__constT_V_nj(HEOS, j, xN_flag); + + double line3 = d3_ndalphardni_dxj_dxk_dDelta__consttau_xi(HEOS, i, j, k, xN_flag); + std::size_t mmax = HEOS.mole_fractions.size(); + if (xN_flag == XN_DEPENDENT){ mmax--; } + for (unsigned int m = 0; m < mmax; m++) + { + line3 -= HEOS.mole_fractions[m]*d3_ndalphardni_dxj_dxk_dDelta__consttau_xi(HEOS, i, m, k, xN_flag); + } + return line1 + line2 + line3; +} CoolPropDbl MixtureDerivatives::d_ndalphardni_dDelta(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, x_N_dependency_flag xN_flag) { // The first line @@ -806,6 +897,34 @@ CoolPropDbl MixtureDerivatives::d3_ndalphardni_dxj_dTau2__constdelta_xi(Helmholt } return term1 + term2 + term3 + term4 + term5; } +CoolPropDbl MixtureDerivatives::d3_ndalphardni_dxj_dDelta2__consttau_xi(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag){ + double term1 = (2*HEOS.d2alphar_dDelta2() + HEOS.delta()*HEOS.d3alphar_dDelta3())*HEOS.Reducing->d_PSI_rho_dxj(HEOS.mole_fractions, i, j, xN_flag); + double term2 = (2*d3alphar_dxi_dDelta2(HEOS, j, xN_flag)+HEOS.delta()*d4alphar_dxi_dDelta3(HEOS, j, xN_flag))*HEOS.Reducing->PSI_rho(HEOS.mole_fractions, i, xN_flag); + double term3 = HEOS.tau()*HEOS.d3alphar_dDelta2_dTau()*HEOS.Reducing->d_PSI_T_dxj(HEOS.mole_fractions, i, j, xN_flag); + double term4 = HEOS.tau()*d4alphar_dxi_dDelta2_dTau(HEOS, j, xN_flag)*HEOS.Reducing->PSI_T(HEOS.mole_fractions, i, xN_flag); + double term5 = d4alphar_dxi_dxj_dDelta2(HEOS, i, j, xN_flag); + std::size_t kmax = HEOS.mole_fractions.size(); + if (xN_flag == XN_DEPENDENT){ kmax--; } + for (unsigned int k = 0; k < kmax; k++) + { + term5 -= HEOS.mole_fractions[k]*d4alphar_dxi_dxj_dDelta2(HEOS, k, j, xN_flag) + Kronecker_delta(k, j)*d3alphar_dxi_dDelta2(HEOS, k, xN_flag); + } + return term1 + term2 + term3 + term4 + term5; +} +CoolPropDbl MixtureDerivatives::d3_ndalphardni_dxj_dDelta_dTau__constxi(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag){ + double term1 = (HEOS.d2alphar_dDelta_dTau() + HEOS.delta()*HEOS.d3alphar_dDelta2_dTau())*HEOS.Reducing->d_PSI_rho_dxj(HEOS.mole_fractions, i, j, xN_flag); + double term2 = (d3alphar_dxi_dDelta_dTau(HEOS, j, xN_flag)+HEOS.delta()*d4alphar_dxi_dDelta2_dTau(HEOS, j, xN_flag))*HEOS.Reducing->PSI_rho(HEOS.mole_fractions, i, xN_flag); + double term3 = (HEOS.tau()*HEOS.d3alphar_dDelta_dTau2()+HEOS.d2alphar_dDelta_dTau())*HEOS.Reducing->d_PSI_T_dxj(HEOS.mole_fractions, i, j, xN_flag); + double term4 = (HEOS.tau()*d4alphar_dxi_dDelta_dTau2(HEOS, j, xN_flag)+d3alphar_dxi_dDelta_dTau(HEOS,j,xN_flag))*HEOS.Reducing->PSI_T(HEOS.mole_fractions, i, xN_flag); + double term5 = d4alphar_dxi_dxj_dDelta_dTau(HEOS, i, j, xN_flag); + std::size_t kmax = HEOS.mole_fractions.size(); + if (xN_flag == XN_DEPENDENT){ kmax--; } + for (unsigned int k = 0; k < kmax; k++) + { + term5 -= HEOS.mole_fractions[k]*d4alphar_dxi_dxj_dDelta_dTau(HEOS, k, j, xN_flag) + Kronecker_delta(k, j)*d3alphar_dxi_dDelta_dTau(HEOS, k, xN_flag); + } + return term1 + term2 + term3 + term4 + term5; +} CoolPropDbl MixtureDerivatives::d_ndalphardni_dTau(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, x_N_dependency_flag xN_flag) { @@ -847,6 +966,45 @@ CoolPropDbl MixtureDerivatives::d2_ndalphardni_dxj_dxk__constdelta_tau_xi(Helmho return term1 + term2 + term3 + term4 + term5 + term6 + term7; } +CoolPropDbl MixtureDerivatives::d3_ndalphardni_dxj_dxk_dTau__constdelta_xi(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, std::size_t k, x_N_dependency_flag xN_flag) +{ + double term1a = HEOS.delta()*d3alphar_dxi_dDelta_dTau(HEOS, j, xN_flag)*HEOS.Reducing->d_PSI_rho_dxj(HEOS.mole_fractions, i, k, xN_flag); + double term1b = HEOS.delta()*d4alphar_dxi_dxj_dDelta_dTau(HEOS, j, k, xN_flag)*HEOS.Reducing->PSI_rho(HEOS.mole_fractions, i, xN_flag); + double term1c = HEOS.delta()*HEOS.d2alphar_dDelta_dTau()*HEOS.Reducing->d2_PSI_rho_dxj_dxk(HEOS.mole_fractions, i, j, k, xN_flag); + double term1d = HEOS.delta()*d3alphar_dxi_dDelta_dTau(HEOS, k, xN_flag)*HEOS.Reducing->d_PSI_rho_dxj(HEOS.mole_fractions, i, j, xN_flag); + double term1 = term1a + term1b + term1c + term1d; + + double term2a = (HEOS.tau()*d3alphar_dxi_dTau2(HEOS, j, xN_flag) + d2alphar_dxi_dTau(HEOS, j,xN_flag))*HEOS.Reducing->d_PSI_T_dxj(HEOS.mole_fractions, i, k, xN_flag); + double term2b = (HEOS.tau()*d4alphar_dxi_dxj_dTau2(HEOS, j, k, xN_flag) + d3alphar_dxi_dxj_dTau(HEOS, j, k, xN_flag))*HEOS.Reducing->PSI_T(HEOS.mole_fractions, i, xN_flag); + double term2c = (HEOS.tau()*HEOS.d2alphar_dTau2() + HEOS.dalphar_dTau())*HEOS.Reducing->d2_PSI_T_dxj_dxk(HEOS.mole_fractions, i, j, k, xN_flag); + double term2d = (HEOS.tau()*d3alphar_dxi_dTau2(HEOS, k, xN_flag) + d2alphar_dxi_dTau(HEOS, k, xN_flag))*HEOS.Reducing->d_PSI_T_dxj(HEOS.mole_fractions, i, j, xN_flag); + double term2 = term2a + term2b + term2c + term2d; + + /// All derivatives with dalphar)/(dxi,dxj,dxk) are zero + double term3 = -d3alphar_dxi_dxj_dTau(HEOS, j, k, xN_flag); + return term1 + term2 + term3; +} + +CoolPropDbl MixtureDerivatives::d3_ndalphardni_dxj_dxk_dDelta__consttau_xi(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, std::size_t k, x_N_dependency_flag xN_flag) +{ + double term1a = (HEOS.delta()*d3alphar_dxi_dDelta2(HEOS, j, xN_flag) + d2alphar_dxi_dDelta(HEOS, j, xN_flag))*HEOS.Reducing->d_PSI_rho_dxj(HEOS.mole_fractions, i, k, xN_flag); + double term1b = (HEOS.delta()*d4alphar_dxi_dxj_dDelta2(HEOS, j, k, xN_flag) + d3alphar_dxi_dxj_dDelta(HEOS, j, k, xN_flag))*HEOS.Reducing->PSI_rho(HEOS.mole_fractions, i, xN_flag); + double term1c = (HEOS.delta()*HEOS.d2alphar_dDelta2() + HEOS.dalphar_dDelta())*HEOS.Reducing->d2_PSI_rho_dxj_dxk(HEOS.mole_fractions, i, j, k, xN_flag); + double term1d = (HEOS.delta()*d3alphar_dxi_dDelta2(HEOS, k, xN_flag) + d2alphar_dxi_dDelta(HEOS, k, xN_flag))*HEOS.Reducing->d_PSI_rho_dxj(HEOS.mole_fractions, i, j, xN_flag); + double term1 = term1a + term1b + term1c + term1d; + + double term2a = HEOS.tau()*d3alphar_dxi_dDelta_dTau(HEOS, j, xN_flag)*HEOS.Reducing->d_PSI_T_dxj(HEOS.mole_fractions, i, k, xN_flag); + double term2b = HEOS.tau()*d4alphar_dxi_dxj_dDelta_dTau(HEOS, j, k, xN_flag)*HEOS.Reducing->PSI_T(HEOS.mole_fractions, i, xN_flag); + double term2c = HEOS.tau()*HEOS.d2alphar_dDelta_dTau()*HEOS.Reducing->d2_PSI_T_dxj_dxk(HEOS.mole_fractions, i, j, k, xN_flag); + double term2d = HEOS.tau()*d3alphar_dxi_dDelta_dTau(HEOS, k, xN_flag)*HEOS.Reducing->d_PSI_T_dxj(HEOS.mole_fractions, i, j, xN_flag); + double term2 = term2a + term2b + term2c + term2d; + + /// All derivatives with dalphar)/(dxi,dxj,dxk) are zero + double term3 = -d3alphar_dxi_dxj_dTau(HEOS, j, k, xN_flag); + return term1 + term2 + term3; +} + + } /* namespace CoolProp */ #ifdef ENABLE_CATCH diff --git a/src/Backends/Helmholtz/MixtureDerivatives.h b/src/Backends/Helmholtz/MixtureDerivatives.h index 07b48adc..d20d47bf 100644 --- a/src/Backends/Helmholtz/MixtureDerivatives.h +++ b/src/Backends/Helmholtz/MixtureDerivatives.h @@ -71,7 +71,6 @@ class MixtureDerivatives{ static CoolPropDbl d4alphar_dxi_dxj_dDelta2(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag); static CoolPropDbl d4alphar_dxi_dxj_dDelta_dTau(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag); static CoolPropDbl d4alphar_dxi_dxj_dTau2(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag); - static CoolPropDbl d4alphardxidxjdxk(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, std::size_t k, x_N_dependency_flag xN_flag); /** \brief GERG 2004 Monograph equation 7.63 @@ -162,9 +161,14 @@ class MixtureDerivatives{ static CoolPropDbl ndln_fugacity_i_dnj__constT_V_xi(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag); static CoolPropDbl d_ndln_fugacity_i_dnj_dtau__constdelta_x(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag); static CoolPropDbl d2_ndln_fugacity_i_dnj_dtau2__constdelta_x(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag); + static CoolPropDbl d2_ndln_fugacity_i_dnj_ddelta_dtau__constx(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag); + static CoolPropDbl d2_ndln_fugacity_i_dnj_ddelta2__consttau_x(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag); static CoolPropDbl d_ndln_fugacity_i_dnj_ddelta__consttau_x(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag); static CoolPropDbl d_ndln_fugacity_i_dnj_ddxk__consttau_delta(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, std::size_t k, x_N_dependency_flag xN_flag); + static CoolPropDbl d2_ndln_fugacity_i_dnj_dxk_dTau__constdelta(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, std::size_t k, x_N_dependency_flag xN_flag); + static CoolPropDbl d2_ndln_fugacity_i_dnj_dxk_dDelta__consttau(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, std::size_t k, x_N_dependency_flag xN_flag); + static CoolPropDbl nd_ndln_fugacity_i_dnj_dnk__constT_V_xi(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, std::size_t k, x_N_dependency_flag xN_flag); static CoolPropDbl nAij(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag){ @@ -175,6 +179,38 @@ class MixtureDerivatives{ CoolPropDbl RT = HEOS.gas_constant()*HEOS.T(); return 1/RT*nd_ndln_fugacity_i_dnj_dnk__constT_V_xi(HEOS, i, j, k, xN_flag) - nAij(HEOS, i, j, xN_flag); } + static CoolPropDbl d_nAij_dX(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag, parameters WRT) + { + if (WRT == iTau){ + return 1/(HEOS.gas_constant()*HEOS.T_reducing())*MixtureDerivatives::ndln_fugacity_i_dnj__constT_V_xi(HEOS, i, j, xN_flag) + + 1/(HEOS.gas_constant()*HEOS.T())*MixtureDerivatives::d_ndln_fugacity_i_dnj_dtau__constdelta_x(HEOS, i, j, xN_flag); + } + else if (WRT == iDelta){ + return 1/(HEOS.gas_constant()*HEOS.T())*MixtureDerivatives::d_ndln_fugacity_i_dnj_ddelta__consttau_x(HEOS, i, j, xN_flag); + } + else{ + throw ValueError(format("wrong WRT")); + } + } + static CoolPropDbl d_n2Aijk_dX(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, std::size_t k, x_N_dependency_flag xN_flag, parameters WRT){ + CoolPropDbl RT = HEOS.gas_constant()*HEOS.T(); + CoolPropDbl summer = -d_nAij_dX(HEOS, i, j, xN_flag, WRT); + if (WRT == iTau){ + summer += d2_ndln_fugacity_i_dnj_dtau2__constdelta_x(HEOS, i, j, xN_flag)*ndtaudni__constT_V_nj(HEOS, k, xN_flag); + summer += d_ndln_fugacity_i_dnj_dtau__constdelta_x(HEOS, i, j, xN_flag)*d_ndtaudni_dTau(HEOS, k, xN_flag); + summer += d2_ndln_fugacity_i_dnj_ddelta_dtau__constx(HEOS, i, j, xN_flag)*nddeltadni__constT_V_nj(HEOS, k, xN_flag); + summer += d2_ndln_fugacity_i_dnj_dxk_dTau__constdelta(HEOS, i, j, k, xN_flag); + std::size_t mmax = HEOS.mole_fractions.size(); + if (xN_flag == XN_DEPENDENT){ mmax--; } + for (std::size_t m = 0; m < mmax; ++m){ + summer -= HEOS.mole_fractions[m]*d2_ndln_fugacity_i_dnj_dxk_dTau__constdelta(HEOS, i, j, m, xN_flag); + } + } + else if (WRT== iDelta){ + return -1000; + } + return summer; + } static Eigen::MatrixXd Lstar(HelmholtzEOSMixtureBackend &HEOS, x_N_dependency_flag xN_flag){ std::size_t N = HEOS.mole_fractions.size(); Eigen::MatrixXd L; @@ -198,16 +234,7 @@ class MixtureDerivatives{ Eigen::MatrixXd dLstar_dX(N, N); for (std::size_t i = 0; i < N; ++i){ for (std::size_t j = i; j < N; ++j){ - if (WRT == iTau){ - dLstar_dX(i, j) = 1/(HEOS.gas_constant()*HEOS.T_reducing())*MixtureDerivatives::ndln_fugacity_i_dnj__constT_V_xi(HEOS, i, j, xN_flag) - + 1/(HEOS.gas_constant()*HEOS.T())*MixtureDerivatives::d_ndln_fugacity_i_dnj_dtau__constdelta_x(HEOS, i, j, xN_flag); - } - else if (WRT == iDelta){ - dLstar_dX(i, j) = 1/(HEOS.gas_constant()*HEOS.T())*MixtureDerivatives::d_ndln_fugacity_i_dnj_ddelta__consttau_x(HEOS, i, j, xN_flag); - } - else{ - throw ValueError(format("dLstar_dX invalid WRT")); - } + dLstar_dX(i, j) = d_nAij_dX(HEOS, i, j, xN_flag, WRT); } } // Fill in the symmetric elements @@ -264,6 +291,35 @@ class MixtureDerivatives{ } return M; } + static Eigen::MatrixXd dMstar_dX(HelmholtzEOSMixtureBackend &HEOS, x_N_dependency_flag xN_flag, parameters WRT){ + + std::size_t N = HEOS.mole_fractions.size(); + Eigen::MatrixXd L = Lstar(HEOS, xN_flag), + dL_dX = dLstar_dX(HEOS, xN_flag, WRT), + dMstar = dL_dX, + adjL = adjugate(L), + d_adjL_dX = adjugate_derivative(L, dL_dX); + + // Last row + for (std::size_t i = 0; i < N; ++i){ + Eigen::MatrixXd n2dLdni(N, N), d_n2dLdni_dX(N, N); + for (std::size_t j = 0; j < N; ++j){ + for (std::size_t k = j; k < N; ++k){ + n2dLdni(j, k) = n2Aijk(HEOS, j, k, i, xN_flag); + d_n2dLdni_dX(j, k) = d_n2Aijk_dX(HEOS, j, k, i, xN_flag, WRT); + } + } + // Fill in the symmetric elements + for (std::size_t j = 0; j < N; ++j){ + for (std::size_t k = 0; k < j; ++k){ + n2dLdni(j, k) = n2dLdni(k, j); + d_n2dLdni_dX(j, k) = d_n2dLdni_dX(k, j); + } + } + dMstar(N-1, i) = (n2dLdni*d_adjL_dX + adjL*d_n2dLdni_dX).trace(); + } + return dMstar; + } /** \brief Table B4, Kunz, JCED, 2012 for the original term and the subsequent substitutions * @@ -497,6 +553,8 @@ class MixtureDerivatives{ */ static CoolPropDbl d2_ndalphardni_dxj_dTau__constdelta_xi(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag); static CoolPropDbl d3_ndalphardni_dxj_dTau2__constdelta_xi(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag); + static CoolPropDbl d3_ndalphardni_dxj_dDelta2__consttau_xi(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag); + static CoolPropDbl d3_ndalphardni_dxj_dDelta_dTau__constxi(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag); /** \brief GERG 2004 Monograph equation 7.41 * @@ -540,6 +598,8 @@ class MixtureDerivatives{ */ static CoolPropDbl d_nd_ndalphardni_dnj_dTau__constdelta_x(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag); static CoolPropDbl d2_nd_ndalphardni_dnj_dTau2__constdelta_x(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag); + static CoolPropDbl d2_nd_ndalphardni_dnj_dDelta2__consttau_x(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag); + static CoolPropDbl d2_nd_ndalphardni_dnj_dDelta_dTau__constx(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag); /* \brief \f$\delta\f$ derivative of GERG 2004 7.47 * @@ -550,6 +610,9 @@ class MixtureDerivatives{ * */ static CoolPropDbl d_nd_ndalphardni_dnj_dxk__consttau_delta(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, std::size_t k, x_N_dependency_flag xN_flag); + + static CoolPropDbl d2_nd_ndalphardni_dnj_dxk_dTau__constdelta(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, std::size_t k, x_N_dependency_flag xN_flag); + static CoolPropDbl d2_nd_ndalphardni_dnj_dxk_dDelta__consttau(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, std::size_t k, x_N_dependency_flag xN_flag); /** \brief GERG 2004 Monograph equation 7.48 * @@ -567,6 +630,8 @@ class MixtureDerivatives{ static CoolPropDbl d_nddeltadni_dxj__constdelta_tau(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag); + static CoolPropDbl d2_nddeltadni_dxj_dDelta__consttau(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag); + /** \brief GERG 2004 Monograph equation 7.49 * * The derivative term @@ -583,6 +648,8 @@ class MixtureDerivatives{ static CoolPropDbl d_ndtaudni_dxj__constdelta_tau(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag); + static CoolPropDbl d2_ndtaudni_dxj_dTau__constdelta(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag); + /** \brief GERG 2004 Monograph equation 7.52 * * The derivative term @@ -605,6 +672,9 @@ class MixtureDerivatives{ */ static CoolPropDbl d2_ndalphardni_dxj_dxk__constdelta_tau_xi(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, std::size_t k, x_N_dependency_flag xN_flag); + static CoolPropDbl d3_ndalphardni_dxj_dxk_dTau__constdelta_xi(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, std::size_t k, x_N_dependency_flag xN_flag); + static CoolPropDbl d3_ndalphardni_dxj_dxk_dDelta__consttau_xi(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, std::size_t k, x_N_dependency_flag xN_flag); + }; /* class MixtureDerivatives */ } /* namepsace CoolProp*/ From 14b1bb829ea890e7d3cb0901f71a8dc030379e66 Mon Sep 17 00:00:00 2001 From: Ian Bell Date: Wed, 12 Aug 2015 23:42:46 -0600 Subject: [PATCH 2/8] Fix tau derivative of n2Aijk --- src/Backends/Helmholtz/MixtureDerivatives.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/Backends/Helmholtz/MixtureDerivatives.h b/src/Backends/Helmholtz/MixtureDerivatives.h index d20d47bf..d3caafb7 100644 --- a/src/Backends/Helmholtz/MixtureDerivatives.h +++ b/src/Backends/Helmholtz/MixtureDerivatives.h @@ -194,8 +194,8 @@ class MixtureDerivatives{ } static CoolPropDbl d_n2Aijk_dX(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, std::size_t k, x_N_dependency_flag xN_flag, parameters WRT){ CoolPropDbl RT = HEOS.gas_constant()*HEOS.T(); - CoolPropDbl summer = -d_nAij_dX(HEOS, i, j, xN_flag, WRT); if (WRT == iTau){ + double summer = 0; summer += d2_ndln_fugacity_i_dnj_dtau2__constdelta_x(HEOS, i, j, xN_flag)*ndtaudni__constT_V_nj(HEOS, k, xN_flag); summer += d_ndln_fugacity_i_dnj_dtau__constdelta_x(HEOS, i, j, xN_flag)*d_ndtaudni_dTau(HEOS, k, xN_flag); summer += d2_ndln_fugacity_i_dnj_ddelta_dtau__constx(HEOS, i, j, xN_flag)*nddeltadni__constT_V_nj(HEOS, k, xN_flag); @@ -205,11 +205,11 @@ class MixtureDerivatives{ for (std::size_t m = 0; m < mmax; ++m){ summer -= HEOS.mole_fractions[m]*d2_ndln_fugacity_i_dnj_dxk_dTau__constdelta(HEOS, i, j, m, xN_flag); } + return 1/RT*summer + 1/(HEOS.gas_constant()*HEOS.T_reducing())*nd_ndln_fugacity_i_dnj_dnk__constT_V_xi(HEOS,i,j,k,xN_flag) - d_nAij_dX(HEOS, i, j, xN_flag, WRT); } else if (WRT== iDelta){ return -1000; } - return summer; } static Eigen::MatrixXd Lstar(HelmholtzEOSMixtureBackend &HEOS, x_N_dependency_flag xN_flag){ std::size_t N = HEOS.mole_fractions.size(); From a714857bb02c6895486c5a2a2f84fc61127a4b9c Mon Sep 17 00:00:00 2001 From: Ian Bell Date: Fri, 14 Aug 2015 08:30:18 -0600 Subject: [PATCH 3/8] Fix derivatives for xxdelta --- src/Backends/Helmholtz/MixtureDerivatives.cpp | 213 +++++++++++++++++- 1 file changed, 201 insertions(+), 12 deletions(-) diff --git a/src/Backends/Helmholtz/MixtureDerivatives.cpp b/src/Backends/Helmholtz/MixtureDerivatives.cpp index 4db8220c..043da17f 100644 --- a/src/Backends/Helmholtz/MixtureDerivatives.cpp +++ b/src/Backends/Helmholtz/MixtureDerivatives.cpp @@ -956,13 +956,9 @@ CoolPropDbl MixtureDerivatives::d2_ndalphardni_dxj_dxk__constdelta_tau_xi(Helmho double term5 = HEOS.tau()*d3alphar_dxi_dxj_dTau(HEOS, j, k, xN_flag)*HEOS.Reducing->PSI_T(HEOS.mole_fractions, i, xN_flag); double term6 = HEOS.tau()*HEOS.dalphar_dTau()*HEOS.Reducing->d2_PSI_T_dxj_dxk(HEOS.mole_fractions, i, j, k, xN_flag); - double term7 = d3alphardxidxjdxk(HEOS, i, j, k, xN_flag) - 2*d2alphardxidxj(HEOS, j, k, xN_flag); - std::size_t mmax = HEOS.mole_fractions.size(); - if (xN_flag == XN_DEPENDENT){ mmax--; } - for (unsigned int m = 0; m < mmax; m++) - { - term7 -= HEOS.mole_fractions[m]*d3alphardxidxjdxk(HEOS, m, j, k, xN_flag); - } + /// All derivatives with dalphar/(dxi,dxj,dxk) are zero + double term7 = - 2*d2alphardxidxj(HEOS, j, k, xN_flag); + return term1 + term2 + term3 + term4 + term5 + term6 + term7; } @@ -980,8 +976,8 @@ CoolPropDbl MixtureDerivatives::d3_ndalphardni_dxj_dxk_dTau__constdelta_xi(Helmh double term2d = (HEOS.tau()*d3alphar_dxi_dTau2(HEOS, k, xN_flag) + d2alphar_dxi_dTau(HEOS, k, xN_flag))*HEOS.Reducing->d_PSI_T_dxj(HEOS.mole_fractions, i, j, xN_flag); double term2 = term2a + term2b + term2c + term2d; - /// All derivatives with dalphar)/(dxi,dxj,dxk) are zero - double term3 = -d3alphar_dxi_dxj_dTau(HEOS, j, k, xN_flag); + /// All derivatives of dalphar/(dxi,dxj,dxk) are zero + double term3 = -2*d3alphar_dxi_dxj_dTau(HEOS, j, k, xN_flag); return term1 + term2 + term3; } @@ -999,8 +995,8 @@ CoolPropDbl MixtureDerivatives::d3_ndalphardni_dxj_dxk_dDelta__consttau_xi(Helmh double term2d = HEOS.tau()*d3alphar_dxi_dDelta_dTau(HEOS, k, xN_flag)*HEOS.Reducing->d_PSI_T_dxj(HEOS.mole_fractions, i, j, xN_flag); double term2 = term2a + term2b + term2c + term2d; - /// All derivatives with dalphar)/(dxi,dxj,dxk) are zero - double term3 = -d3alphar_dxi_dxj_dTau(HEOS, j, k, xN_flag); + /// All derivatives of dalphar)/(dxi,dxj,dxk) are zero + double term3 = -2*d3alphar_dxi_dxj_dDelta(HEOS, j, k, xN_flag); return term1 + term2 + term3; } @@ -1450,6 +1446,20 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") CAPTURE(analytic); CHECK(err < 1e-7); } + std::ostringstream ss3ff; ss3ff << "d4alphar_dxi_dTau3, i=" << i; + SECTION(ss3ff.str(), "") + { + if (i == Ncomp-1){ break; } + double analytic = MixtureDerivatives::d4alphar_dxi_dTau3(rHEOS, i, xN_flag); + double v1 = MixtureDerivatives::d3alphar_dxi_dTau2(rHEOS_plusT_constrho, i, xN_flag), tau1 = rHEOS_plusT_constrho.tau(); + double v2 = MixtureDerivatives::d3alphar_dxi_dTau2(rHEOS_minusT_constrho, i, xN_flag), tau2 = rHEOS_minusT_constrho.tau(); + double numeric = (v1 - v2)/(tau1 - tau2); + double err = std::abs((numeric-analytic)/analytic); + CAPTURE(numeric); + CAPTURE(analytic); + CHECK(err < 1e-7); + } + std::ostringstream ss3g; ss3g << "d3alphar_dxi_dDelta_Tau, i=" << i; SECTION(ss3g.str(), "") { @@ -1463,6 +1473,33 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") CAPTURE(analytic); CHECK(err < 1e-7); } + std::ostringstream ss3gg; ss3gg << "d4alphar_dxi_dDelta_Tau2, i=" << i; + SECTION(ss3gg.str(), "") + { + if (i == Ncomp-1){ break; } + double analytic = MixtureDerivatives::d4alphar_dxi_dDelta_dTau2(rHEOS, i, xN_flag); + double v1 = MixtureDerivatives::d3alphar_dxi_dDelta_dTau(rHEOS_plusT_constrho, i, xN_flag), tau1 = rHEOS_plusT_constrho.tau(); + double v2 = MixtureDerivatives::d3alphar_dxi_dDelta_dTau(rHEOS_minusT_constrho, i, xN_flag), tau2 = rHEOS_minusT_constrho.tau(); + double numeric = (v1 - v2)/(tau1 - tau2); + double err = std::abs((numeric-analytic)/analytic); + CAPTURE(numeric); + CAPTURE(analytic); + CHECK(err < 1e-7); + } + std::ostringstream ss3ggg; ss3ggg << "d4alphar_dxi_dDelta2_Tau, i=" << i; + SECTION(ss3ggg.str(), "") + { + if (i == Ncomp-1){ break; } + double analytic = MixtureDerivatives::d4alphar_dxi_dDelta2_dTau(rHEOS, i, xN_flag); + double v1 = MixtureDerivatives::d3alphar_dxi_dDelta2(rHEOS_plusT_constrho, i, xN_flag), tau1 = rHEOS_plusT_constrho.tau(); + double v2 = MixtureDerivatives::d3alphar_dxi_dDelta2(rHEOS_minusT_constrho, i, xN_flag), tau2 = rHEOS_minusT_constrho.tau(); + double numeric = (v1 - v2)/(tau1 - tau2); + double err = std::abs((numeric-analytic)/analytic); + CAPTURE(numeric); + CAPTURE(analytic); + CHECK(err < 1e-7); + } + std::ostringstream ss3h; ss3h << "d2_ndalphardni_dTau2, i=" << i; SECTION(ss3h.str(), "") { @@ -1476,6 +1513,19 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") CAPTURE(analytic); CHECK(err < 1e-7); } + std::ostringstream ss3hh; ss3hh << "d3_ndalphardni_dTau3, i=" << i; + SECTION(ss3hh.str(), "") + { + if (i == Ncomp-1){ break; } + double analytic = MixtureDerivatives::d3_ndalphardni_dTau3(rHEOS, i, xN_flag); + double v1 = MixtureDerivatives::d2_ndalphardni_dTau2(rHEOS_plusT_constrho, i, xN_flag), tau1 = rHEOS_plusT_constrho.tau(); + double v2 = MixtureDerivatives::d2_ndalphardni_dTau2(rHEOS_minusT_constrho, i, xN_flag), tau2 = rHEOS_minusT_constrho.tau(); + double numeric = (v1 - v2)/(tau1 - tau2); + double err = std::abs((numeric-analytic)/analytic); + CAPTURE(numeric); + CAPTURE(analytic); + CHECK(err < 1e-7); + } std::ostringstream ss3i; ss3i << "d2_ndalphardni_dDelta_dTau, i=" << i; SECTION(ss3i.str(), "") { @@ -1489,6 +1539,32 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") CAPTURE(analytic); CHECK(err < 1e-7); } + std::ostringstream ss3ii; ss3ii << "d3_ndalphardni_dDelta_dTau2, i=" << i; + SECTION(ss3ii.str(), "") + { + if (i == Ncomp-1){ break; } + double analytic = MixtureDerivatives::d3_ndalphardni_dDelta_dTau2(rHEOS, i, xN_flag); + double v1 = MixtureDerivatives::d2_ndalphardni_dTau2(rHEOS_plusrho_constT, i, xN_flag), delta1 = rHEOS_plusrho_constT.delta(); + double v2 = MixtureDerivatives::d2_ndalphardni_dTau2(rHEOS_minusrho_constT, i, xN_flag), delta2 = rHEOS_minusrho_constT.delta(); + double numeric = (v1 - v2)/(delta1 - delta2); + double err = std::abs((numeric-analytic)/analytic); + CAPTURE(numeric); + CAPTURE(analytic); + CHECK(err < 1e-7); + } + std::ostringstream ss3iii; ss3iii << "d3_ndalphardni_dDelta2_dTau, i=" << i; + SECTION(ss3iii.str(), "") + { + if (i == Ncomp-1){ break; } + double analytic = MixtureDerivatives::d3_ndalphardni_dDelta2_dTau(rHEOS, i, xN_flag); + double v1 = MixtureDerivatives::d2_ndalphardni_dDelta_dTau(rHEOS_plusrho_constT, i, xN_flag), delta1 = rHEOS_plusrho_constT.delta(); + double v2 = MixtureDerivatives::d2_ndalphardni_dDelta_dTau(rHEOS_minusrho_constT, i, xN_flag), delta2 = rHEOS_minusrho_constT.delta(); + double numeric = (v1 - v2)/(delta1 - delta2); + double err = std::abs((numeric-analytic)/analytic); + CAPTURE(numeric); + CAPTURE(analytic); + CHECK(err < 1e-7); + } std::ostringstream ss3j; ss3j << "d2_ndalphardni_dDelta2, i=" << i; SECTION(ss3j.str(), "") { @@ -1502,6 +1578,19 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") CAPTURE(analytic); CHECK(err < 1e-7); } + std::ostringstream ss3k; ss3k << "d3_ndalphardni_dDelta3, i=" << i; + SECTION(ss3k.str(), "") + { + if (i == Ncomp-1){ break; } + double analytic = MixtureDerivatives::d3_ndalphardni_dDelta3(rHEOS, i, xN_flag); + double v1 = MixtureDerivatives::d2_ndalphardni_dDelta2(rHEOS_plusrho_constT, i, xN_flag), delta1 = rHEOS_plusrho_constT.delta(); + double v2 = MixtureDerivatives::d2_ndalphardni_dDelta2(rHEOS_minusrho_constT, i, xN_flag), delta2 = rHEOS_minusrho_constT.delta(); + double numeric = (v1 - v2)/(delta1 - delta2); + double err = std::abs((numeric-analytic)/analytic); + CAPTURE(numeric); + CAPTURE(analytic); + CHECK(err < 1e-7); + } // These derivatives depend on both the i and j indices for (std::size_t j = 0; j < Ncomp; ++j){ @@ -1646,6 +1735,19 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") CAPTURE(analytic); CHECK(err < 1e-7); } + std::ostringstream ss3kk; ss3kk << "d3_ndalphardni_dxj_dDelta2__consttau_xi, i=" << i; + SECTION(ss3kk.str(), "") + { + if (i == Ncomp-1){ break; } + double analytic = MixtureDerivatives::d3_ndalphardni_dxj_dDelta2__consttau_xi(rHEOS, i, j, xN_flag); + double v1 = MixtureDerivatives::d2_ndalphardni_dDelta2(rHEOS_pluszj_consttaudelta, i, xN_flag); + double v2 = MixtureDerivatives::d2_ndalphardni_dDelta2(rHEOS_minuszj_consttaudelta, i, xN_flag); + double numeric = (v1 - v2)/(2*dz); + double err = std::abs((numeric-analytic)/analytic); + CAPTURE(numeric); + CAPTURE(analytic); + CHECK(err < 1e-7); + } std::ostringstream ss3l; ss3l << "d2_ndalphardni_dxj_dTau, i=" << i; SECTION(ss3l.str(), "") { @@ -1659,10 +1761,33 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") CAPTURE(analytic); CHECK(err < 1e-7); } + std::ostringstream ss3ll; ss3ll << "d3_ndalphardni_dxj_dTau2__constdelta_xi, i=" << i; + SECTION(ss3ll.str(), "") + { + double analytic = MixtureDerivatives::d3_ndalphardni_dxj_dTau2__constdelta_xi(rHEOS, i, j, xN_flag); + double v1 = MixtureDerivatives::d2_ndalphardni_dTau2(rHEOS_pluszj_consttaudelta, i, xN_flag); + double v2 = MixtureDerivatives::d2_ndalphardni_dTau2(rHEOS_minuszj_consttaudelta, i, xN_flag); + double numeric = (v1 - v2)/(2*dz); + double err = std::abs((numeric-analytic)/analytic); + CAPTURE(numeric); + CAPTURE(analytic); + CHECK(err < 1e-7); + } + std::ostringstream ss3lll; ss3lll << "d3_ndalphardni_dxj_dDelta_dTau__constxi, i=" << i; + SECTION(ss3lll.str(), "") + { + double analytic = MixtureDerivatives::d3_ndalphardni_dxj_dDelta_dTau__constxi(rHEOS, i, j, xN_flag); + double v1 = MixtureDerivatives::d2_ndalphardni_dDelta_dTau(rHEOS_pluszj_consttaudelta, i, xN_flag); + double v2 = MixtureDerivatives::d2_ndalphardni_dDelta_dTau(rHEOS_minuszj_consttaudelta, i, xN_flag); + double numeric = (v1 - v2)/(2*dz); + double err = std::abs((numeric-analytic)/analytic); + CAPTURE(numeric); + CAPTURE(analytic); + CHECK(err < 1e-7); + } std::ostringstream ss3m; ss3m << "d3alphar_dxi_dxj_dDelta, i=" << i; SECTION(ss3m.str(), "") { - if (i == Ncomp-1){ break; } double analytic = MixtureDerivatives::d3alphar_dxi_dxj_dDelta(rHEOS, i, j, xN_flag); double v1 = MixtureDerivatives::d2alphar_dxi_dDelta(rHEOS_pluszj_consttaudelta, i, xN_flag); double v2 = MixtureDerivatives::d2alphar_dxi_dDelta(rHEOS_minuszj_consttaudelta, i, xN_flag); @@ -1672,6 +1797,30 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") CAPTURE(analytic); CHECK(err < 1e-7); } + std::ostringstream ss3mm; ss3mm << "d4alphar_dxi_dxj_dDelta2, i=" << i; + SECTION(ss3mm.str(), "") + { + double analytic = MixtureDerivatives::d4alphar_dxi_dxj_dDelta2(rHEOS, i, j, xN_flag); + double v1 = MixtureDerivatives::d3alphar_dxi_dDelta2(rHEOS_pluszj_consttaudelta, i, xN_flag); + double v2 = MixtureDerivatives::d3alphar_dxi_dDelta2(rHEOS_minuszj_consttaudelta, i, xN_flag); + double numeric = (v1 - v2)/(2*dz); + double err = std::abs((numeric-analytic)/analytic); + CAPTURE(numeric); + CAPTURE(analytic); + CHECK(err < 1e-7); + } + std::ostringstream ss3mmm; ss3mmm << "d4alphar_dxi_dxj_dDelta_dTau, i=" << i; + SECTION(ss3mmm.str(), "") + { + double analytic = MixtureDerivatives::d4alphar_dxi_dxj_dDelta_dTau(rHEOS, i, j, xN_flag); + double v1 = MixtureDerivatives::d3alphar_dxi_dDelta_dTau(rHEOS_pluszj_consttaudelta, i, xN_flag); + double v2 = MixtureDerivatives::d3alphar_dxi_dDelta_dTau(rHEOS_minuszj_consttaudelta, i, xN_flag); + double numeric = (v1 - v2)/(2*dz); + double err = std::abs((numeric-analytic)/analytic); + CAPTURE(numeric); + CAPTURE(analytic); + CHECK(err < 1e-7); + } std::ostringstream ss3n; ss3n << "d3alphar_dxi_dxj_dTau, i=" << i; SECTION(ss3n.str(), "") { @@ -1685,6 +1834,20 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") CAPTURE(analytic); CHECK(err < 1e-7); } + std::ostringstream ss3nn; ss3nn << "d4alphar_dxi_dxj_dTau2, i=" << i; + SECTION(ss3nn.str(), "") + { + if (i == Ncomp-1){ break; } + double analytic = MixtureDerivatives::d4alphar_dxi_dxj_dTau2(rHEOS, i, j, xN_flag); + double v1 = MixtureDerivatives::d3alphar_dxi_dTau2(rHEOS_pluszj_consttaudelta, i, xN_flag); + double v2 = MixtureDerivatives::d3alphar_dxi_dTau2(rHEOS_minuszj_consttaudelta, i, xN_flag); + double numeric = (v1 - v2)/(2*dz); + double err = std::abs((numeric-analytic)/analytic); + CAPTURE(numeric); + CAPTURE(analytic); + CHECK(err < 1e-7); + } + std::ostringstream ss3o; ss3o << "d_nd_ndalphardni_dnj_dDelta__consttau_x, i=" << i; SECTION(ss3o.str(), "") { @@ -1865,6 +2028,32 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") CAPTURE(analytic); CHECK(err < 1e-7); } + std::ostringstream ss7a; ss7a << "d3_ndalphardni_dxj_dxk_dDelta__consttau_xi, i=" << i << ", j=" << j << ", k=" << k; + SECTION(ss7a.str(), "") + { + if ((xN_flag == XN_DEPENDENT) && (i == Ncomp-1 || j == Ncomp-1 || k == Ncomp-1)){ break; } + double analytic = MixtureDerivatives::d3_ndalphardni_dxj_dxk_dDelta__consttau_xi(rHEOS, i, j, k, xN_flag); + double v1 = MixtureDerivatives::d2_ndalphardni_dxj_dDelta__consttau_xi(rHEOS_pluszk_consttaudelta, i, j, xN_flag); + double v2 = MixtureDerivatives::d2_ndalphardni_dxj_dDelta__consttau_xi(rHEOS_minuszk_consttaudelta, i, j, xN_flag); + double numeric = (v1 - v2)/(2*dz); + double err = std::abs((numeric-analytic)/analytic); + CAPTURE(numeric); + CAPTURE(analytic); + CHECK(err < 1e-7); + } + std::ostringstream ss7b; ss7b << "d3_ndalphardni_dxj_dxk_dTau__constdelta_xi, i=" << i << ", j=" << j << ", k=" << k; + SECTION(ss7b.str(), "") + { + if ((xN_flag == XN_DEPENDENT) && (i == Ncomp-1 || j == Ncomp-1 || k == Ncomp-1)){ break; } + double analytic = MixtureDerivatives::d3_ndalphardni_dxj_dxk_dTau__constdelta_xi(rHEOS, i, j, k, xN_flag); + double v1 = MixtureDerivatives::d2_ndalphardni_dxj_dTau__constdelta_xi(rHEOS_pluszk_consttaudelta, i, j, xN_flag); + double v2 = MixtureDerivatives::d2_ndalphardni_dxj_dTau__constdelta_xi(rHEOS_minuszk_consttaudelta, i, j, xN_flag); + double numeric = (v1 - v2)/(2*dz); + double err = std::abs((numeric-analytic)/analytic); + CAPTURE(numeric); + CAPTURE(analytic); + CHECK(err < 1e-7); + } std::ostringstream ss8; ss8 << "d3alphardxidxjdxk, i=" << i << ", j=" << j << ", k=" << k; SECTION(ss8.str(), "") { From 045ef05c9de142a02086b941379ae2c291379b82 Mon Sep 17 00:00:00 2001 From: Ian Bell Date: Fri, 14 Aug 2015 18:24:05 -0600 Subject: [PATCH 4/8] Fix a bunch of little things with mixture derivatives; most tests pass now --- src/Backends/Helmholtz/MixtureDerivatives.cpp | 326 ++++++++++++------ 1 file changed, 226 insertions(+), 100 deletions(-) diff --git a/src/Backends/Helmholtz/MixtureDerivatives.cpp b/src/Backends/Helmholtz/MixtureDerivatives.cpp index 043da17f..62ce9094 100644 --- a/src/Backends/Helmholtz/MixtureDerivatives.cpp +++ b/src/Backends/Helmholtz/MixtureDerivatives.cpp @@ -658,7 +658,7 @@ CoolPropDbl MixtureDerivatives::d2_nd_ndalphardni_dnj_dDelta_dTau__constx(Helmho double line1 = d3_ndalphardni_dDelta2_dTau(HEOS, i, xN_flag)*nddeltadni__constT_V_nj(HEOS, j, xN_flag); double line2 = d2_ndalphardni_dDelta_dTau(HEOS, i, xN_flag)*d_nddeltadni_dDelta(HEOS, j, xN_flag); double line3 = d2_ndalphardni_dDelta_dTau(HEOS, i, xN_flag)*d_ndtaudni_dTau(HEOS, j, xN_flag); - double line4 = d3_ndalphardni_dDelta2_dTau(HEOS, i, xN_flag)*ndtaudni__constT_V_nj(HEOS, j, xN_flag); + double line4 = d3_ndalphardni_dDelta_dTau2(HEOS, i, xN_flag)*ndtaudni__constT_V_nj(HEOS, j, xN_flag); double summer = 0; std::size_t kmax = HEOS.mole_fractions.size(); if (xN_flag == XN_DEPENDENT){ kmax--; } @@ -713,7 +713,7 @@ CoolPropDbl MixtureDerivatives::d2_nd_ndalphardni_dnj_dxk_dTau__constdelta(Helmh + d2_ndalphardni_dxj_dTau__constdelta_xi(HEOS, i, k, xN_flag)*d_ndtaudni_dTau(HEOS, j, xN_flag) + d3_ndalphardni_dxj_dTau2__constdelta_xi(HEOS, i, k, xN_flag)*ndtaudni__constT_V_nj(HEOS, j, xN_flag); - double line3 = d3_ndalphardni_dxj_dxk_dTau__constdelta_xi(HEOS, i, j, k, xN_flag) ; + double line3 = d3_ndalphardni_dxj_dxk_dTau__constdelta_xi(HEOS, i, j, k, xN_flag) - d2_ndalphardni_dxj_dTau__constdelta_xi(HEOS, i, k, xN_flag); std::size_t mmax = HEOS.mole_fractions.size(); if (xN_flag == XN_DEPENDENT){ mmax--; } @@ -733,7 +733,7 @@ CoolPropDbl MixtureDerivatives::d2_nd_ndalphardni_dnj_dxk_dDelta__consttau(Helmh double line2 = d2_ndalphardni_dDelta_dTau(HEOS, i, xN_flag)*d_ndtaudni_dxj__constdelta_tau(HEOS, j, k, xN_flag) + d3_ndalphardni_dxj_dDelta_dTau__constxi(HEOS, i, k, xN_flag)*ndtaudni__constT_V_nj(HEOS, j, xN_flag); - double line3 = d3_ndalphardni_dxj_dxk_dDelta__consttau_xi(HEOS, i, j, k, xN_flag); + double line3 = d3_ndalphardni_dxj_dxk_dDelta__consttau_xi(HEOS, i, j, k, xN_flag) - d2_ndalphardni_dxj_dDelta__consttau_xi(HEOS, i, k, xN_flag); std::size_t mmax = HEOS.mole_fractions.size(); if (xN_flag == XN_DEPENDENT){ mmax--; } for (unsigned int m = 0; m < mmax; m++) @@ -1008,6 +1008,16 @@ CoolPropDbl MixtureDerivatives::d3_ndalphardni_dxj_dxk_dDelta__consttau_xi(Helmh using namespace CoolProp; +double mix_deriv_err_func(double numeric, double analytic) +{ + if (std::abs(analytic) < DBL_EPSILON){ + return std::abs(numeric - analytic); + } + else{ + return (numeric/analytic-1); + } +} + static bool fluids_set = false; static const std::size_t Ncomp_max = 6; @@ -1163,6 +1173,54 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") HelmholtzEOSMixtureBackend &rHEOS_minusrho_constT = *(HEOS_minusrho_constT[Ncomp][0].get()); const std::vector &z = rHEOS.get_mole_fractions(); + + std::ostringstream ss3o55; ss3o55 << "dLstar_Tau"; + SECTION(ss3o55.str(), "") + { + Eigen::MatrixXd analytic = MixtureDerivatives::dLstar_dX(rHEOS, xN_flag, CoolProp::iTau); + Eigen::MatrixXd m1 = MixtureDerivatives::Lstar(rHEOS_plusT_constrho, xN_flag); double tau1 = rHEOS_plusT_constrho.tau(); + Eigen::MatrixXd m2 = MixtureDerivatives::Lstar(rHEOS_minusT_constrho, xN_flag); double tau2 = rHEOS_minusT_constrho.tau(); + Eigen::MatrixXd numeric = (m1 - m2)/(tau1 - tau2); + double err = ((analytic-numeric).array()/analytic.array()).cwiseAbs().sum()/Ncomp/Ncomp; + CAPTURE(numeric); + CAPTURE(analytic); + CHECK(err < 1e-12); + } + std::ostringstream ss3o5; ss3o5 << "dLstar_dDelta"; + SECTION(ss3o5.str(), "") { + Eigen::MatrixXd analytic = MixtureDerivatives::dLstar_dX(rHEOS, xN_flag, CoolProp::iDelta); + Eigen::MatrixXd m1 = MixtureDerivatives::Lstar(rHEOS_plusrho_constT, xN_flag); double delta1 = rHEOS_plusrho_constT.delta(); + Eigen::MatrixXd m2 = MixtureDerivatives::Lstar(rHEOS_minusrho_constT, xN_flag); double delta2 = rHEOS_minusrho_constT.delta(); + Eigen::MatrixXd numeric = (m1 - m2)/(delta1 - delta2); + double err = ((analytic-numeric).array()/analytic.array()).cwiseAbs().sum()/Ncomp/Ncomp; + CAPTURE(numeric); + CAPTURE(analytic); + CHECK(err < 1e-12); + } + std::ostringstream ss3o6; ss3o6 << "dMstar_dTau"; + SECTION(ss3o6.str(), "") + { + Eigen::MatrixXd analytic = MixtureDerivatives::dMstar_dX(rHEOS, xN_flag, CoolProp::iTau); + Eigen::MatrixXd m1 = MixtureDerivatives::Mstar(rHEOS_plusT_constrho, xN_flag); double tau1 = rHEOS_plusT_constrho.tau(); + Eigen::MatrixXd m2 = MixtureDerivatives::Mstar(rHEOS_minusT_constrho, xN_flag); double tau2 = rHEOS_minusT_constrho.tau(); + Eigen::MatrixXd numeric = (m1 - m2)/(tau1 - tau2); + double err = ((analytic-numeric).array()/analytic.array()).cwiseAbs().sum()/Ncomp/Ncomp; + CAPTURE(numeric); + CAPTURE(analytic); + CHECK(err < 1e-12); + } + std::ostringstream ss3o7; ss3o7 << "dMstar_dDelta"; + SECTION(ss3o7.str(), "") + { + Eigen::MatrixXd analytic = MixtureDerivatives::dMstar_dX(rHEOS, xN_flag, CoolProp::iDelta); + Eigen::MatrixXd m1 = MixtureDerivatives::Mstar(rHEOS_plusrho_constT, xN_flag); double delta1 = rHEOS_plusrho_constT.delta(); + Eigen::MatrixXd m2 = MixtureDerivatives::Mstar(rHEOS_minusrho_constT, xN_flag); double delta2 = rHEOS_minusrho_constT.delta(); + Eigen::MatrixXd numeric = (m1 - m2)/(delta1 - delta2); + double err = ((analytic-numeric).array()/analytic.array()).cwiseAbs().sum()/Ncomp/Ncomp; + CAPTURE(numeric); + CAPTURE(analytic); + CHECK(err < 1e-12); + } // These ones only require the i index for (std::size_t i = 0; i< Ncomp; ++i) @@ -1181,7 +1239,7 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") 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); + double err = mix_deriv_err_func(numeric, analytic); CAPTURE(numeric); CAPTURE(analytic); CHECK(err < 1e-8); @@ -1194,7 +1252,7 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") 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); + double err = mix_deriv_err_func(numeric, analytic); CAPTURE(numeric); CAPTURE(analytic); CHECK(err < 1e-8); @@ -1207,7 +1265,7 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") 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); + double err = mix_deriv_err_func(numeric, analytic); CAPTURE(numeric); CAPTURE(analytic); CHECK(err < 1e-8); @@ -1220,7 +1278,7 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") 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); + double err = mix_deriv_err_func(numeric, analytic); CHECK(err < 1e-8); } /*std::ostringstream ss1; @@ -1270,7 +1328,7 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") 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); + double err = mix_deriv_err_func(numeric, analytic); CHECK(err < 1e-8); } std::ostringstream ss3a; @@ -1282,7 +1340,7 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") 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); + double err = mix_deriv_err_func(numeric, analytic); CHECK(err < 1e-8); } std::ostringstream ss4a; @@ -1294,7 +1352,7 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") 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); + double err = mix_deriv_err_func(numeric, analytic); CHECK(err < 1e-8); } std::ostringstream ss4; @@ -1305,7 +1363,7 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") 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); + double err = mix_deriv_err_func(numeric, analytic); CHECK(err < 1e-8); } std::ostringstream ss5; @@ -1317,7 +1375,7 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") double v1 = rHEOS_pluszi.p(); double v2 = rHEOS_minuszi.p(); double numeric = (v1 - v2)/(2*dz); - double err = std::abs((numeric-analytic)/analytic); + double err = mix_deriv_err_func(numeric, analytic); CAPTURE(numeric); CAPTURE(analytic); CHECK(err < 1e-8); @@ -1331,7 +1389,7 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") double v1 = rHEOS_pluszi.tau(); double v2 = rHEOS_minuszi.tau(); double numeric = (v1 - v2)/(2*dz); - double err = std::abs((numeric-analytic)/analytic); + double err = mix_deriv_err_func(numeric, analytic); CAPTURE(numeric); CAPTURE(analytic); CHECK(err < 1e-8); @@ -1345,7 +1403,7 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") double v1 = rHEOS_pluszi.delta(); double v2 = rHEOS_minuszi.delta(); double numeric = (v1 - v2)/(2*dz); - double err = std::abs((numeric-analytic)/analytic); + double err = mix_deriv_err_func(numeric, analytic); CAPTURE(numeric); CAPTURE(analytic); CHECK(err < 1e-8); @@ -1359,7 +1417,7 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") 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); + double err = mix_deriv_err_func(numeric, analytic); CAPTURE(numeric); CAPTURE(analytic); CHECK(err < 1e-8); @@ -1373,7 +1431,7 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") double v1 = rHEOS_pluszi.Reducing->Tr(rHEOS_pluszi.get_mole_fractions()); double v2 = rHEOS_minuszi.Reducing->Tr(rHEOS_minuszi.get_mole_fractions()); double numeric = (v1 - v2)/(2*dz); - double err = std::abs((numeric-analytic)/analytic); + double err = mix_deriv_err_func(numeric, analytic); CAPTURE(numeric); CAPTURE(analytic); CHECK(err < 1e-8); @@ -1387,7 +1445,7 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") double v1 = rHEOS_pluszi.Reducing->rhormolar(rHEOS_pluszi.get_mole_fractions()); double v2 = rHEOS_minuszi.Reducing->rhormolar(rHEOS_minuszi.get_mole_fractions()); double numeric = (v1 - v2)/(2*dz); - double err = std::abs((numeric-analytic)/analytic); + double err = mix_deriv_err_func(numeric, analytic); CAPTURE(numeric); CAPTURE(analytic); CHECK(err < 1e-8); @@ -1401,7 +1459,7 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") 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); double numeric = (v1 - v2)/(2*dz); - double err = std::abs((numeric-analytic)/analytic); + double err = mix_deriv_err_func(numeric, analytic); CAPTURE(numeric); CAPTURE(analytic); CHECK(err < 1e-8); @@ -1415,7 +1473,7 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") 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); + double err = mix_deriv_err_func(numeric, analytic); CAPTURE(numeric); CAPTURE(analytic); CHECK(err < 1e-7); @@ -1428,7 +1486,7 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") double v1 = MixtureDerivatives::d2alphar_dxi_dDelta(rHEOS_plusrho_constT, i, xN_flag), delta1 = rHEOS_plusrho_constT.delta(); double v2 = MixtureDerivatives::d2alphar_dxi_dDelta(rHEOS_minusrho_constT, i, xN_flag), delta2 = rHEOS_minusrho_constT.delta(); double numeric = (v1 - v2)/(delta1 - delta2); - double err = std::abs((numeric-analytic)/analytic); + double err = mix_deriv_err_func(numeric, analytic); CAPTURE(numeric); CAPTURE(analytic); CHECK(err < 1e-7); @@ -1441,7 +1499,7 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") double v1 = MixtureDerivatives::d2alphar_dxi_dTau(rHEOS_plusT_constrho, i, xN_flag), tau1 = rHEOS_plusT_constrho.tau(); double v2 = MixtureDerivatives::d2alphar_dxi_dTau(rHEOS_minusT_constrho, i, xN_flag), tau2 = rHEOS_minusT_constrho.tau(); double numeric = (v1 - v2)/(tau1 - tau2); - double err = std::abs((numeric-analytic)/analytic); + double err = mix_deriv_err_func(numeric, analytic); CAPTURE(numeric); CAPTURE(analytic); CHECK(err < 1e-7); @@ -1454,7 +1512,7 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") double v1 = MixtureDerivatives::d3alphar_dxi_dTau2(rHEOS_plusT_constrho, i, xN_flag), tau1 = rHEOS_plusT_constrho.tau(); double v2 = MixtureDerivatives::d3alphar_dxi_dTau2(rHEOS_minusT_constrho, i, xN_flag), tau2 = rHEOS_minusT_constrho.tau(); double numeric = (v1 - v2)/(tau1 - tau2); - double err = std::abs((numeric-analytic)/analytic); + double err = mix_deriv_err_func(numeric, analytic); CAPTURE(numeric); CAPTURE(analytic); CHECK(err < 1e-7); @@ -1468,7 +1526,7 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") double v1 = MixtureDerivatives::d2alphar_dxi_dDelta(rHEOS_plusT_constrho, i, xN_flag), tau1 = rHEOS_plusT_constrho.tau(); double v2 = MixtureDerivatives::d2alphar_dxi_dDelta(rHEOS_minusT_constrho, i, xN_flag), tau2 = rHEOS_minusT_constrho.tau(); double numeric = (v1 - v2)/(tau1 - tau2); - double err = std::abs((numeric-analytic)/analytic); + double err = mix_deriv_err_func(numeric, analytic); CAPTURE(numeric); CAPTURE(analytic); CHECK(err < 1e-7); @@ -1481,7 +1539,7 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") double v1 = MixtureDerivatives::d3alphar_dxi_dDelta_dTau(rHEOS_plusT_constrho, i, xN_flag), tau1 = rHEOS_plusT_constrho.tau(); double v2 = MixtureDerivatives::d3alphar_dxi_dDelta_dTau(rHEOS_minusT_constrho, i, xN_flag), tau2 = rHEOS_minusT_constrho.tau(); double numeric = (v1 - v2)/(tau1 - tau2); - double err = std::abs((numeric-analytic)/analytic); + double err = mix_deriv_err_func(numeric, analytic); CAPTURE(numeric); CAPTURE(analytic); CHECK(err < 1e-7); @@ -1494,7 +1552,7 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") double v1 = MixtureDerivatives::d3alphar_dxi_dDelta2(rHEOS_plusT_constrho, i, xN_flag), tau1 = rHEOS_plusT_constrho.tau(); double v2 = MixtureDerivatives::d3alphar_dxi_dDelta2(rHEOS_minusT_constrho, i, xN_flag), tau2 = rHEOS_minusT_constrho.tau(); double numeric = (v1 - v2)/(tau1 - tau2); - double err = std::abs((numeric-analytic)/analytic); + double err = mix_deriv_err_func(numeric, analytic); CAPTURE(numeric); CAPTURE(analytic); CHECK(err < 1e-7); @@ -1508,7 +1566,7 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") double v1 = MixtureDerivatives::d_ndalphardni_dTau(rHEOS_plusT_constrho, i, xN_flag), tau1 = rHEOS_plusT_constrho.tau(); double v2 = MixtureDerivatives::d_ndalphardni_dTau(rHEOS_minusT_constrho, i, xN_flag), tau2 = rHEOS_minusT_constrho.tau(); double numeric = (v1 - v2)/(tau1 - tau2); - double err = std::abs((numeric-analytic)/analytic); + double err = mix_deriv_err_func(numeric, analytic); CAPTURE(numeric); CAPTURE(analytic); CHECK(err < 1e-7); @@ -1521,7 +1579,7 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") double v1 = MixtureDerivatives::d2_ndalphardni_dTau2(rHEOS_plusT_constrho, i, xN_flag), tau1 = rHEOS_plusT_constrho.tau(); double v2 = MixtureDerivatives::d2_ndalphardni_dTau2(rHEOS_minusT_constrho, i, xN_flag), tau2 = rHEOS_minusT_constrho.tau(); double numeric = (v1 - v2)/(tau1 - tau2); - double err = std::abs((numeric-analytic)/analytic); + double err = mix_deriv_err_func(numeric, analytic); CAPTURE(numeric); CAPTURE(analytic); CHECK(err < 1e-7); @@ -1534,7 +1592,7 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") double v1 = MixtureDerivatives::d_ndalphardni_dTau(rHEOS_plusrho_constT, i, xN_flag), delta1 = rHEOS_plusrho_constT.delta(); double v2 = MixtureDerivatives::d_ndalphardni_dTau(rHEOS_minusrho_constT, i, xN_flag), delta2 = rHEOS_minusrho_constT.delta(); double numeric = (v1 - v2)/(delta1 - delta2); - double err = std::abs((numeric-analytic)/analytic); + double err = mix_deriv_err_func(numeric, analytic); CAPTURE(numeric); CAPTURE(analytic); CHECK(err < 1e-7); @@ -1547,7 +1605,7 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") double v1 = MixtureDerivatives::d2_ndalphardni_dTau2(rHEOS_plusrho_constT, i, xN_flag), delta1 = rHEOS_plusrho_constT.delta(); double v2 = MixtureDerivatives::d2_ndalphardni_dTau2(rHEOS_minusrho_constT, i, xN_flag), delta2 = rHEOS_minusrho_constT.delta(); double numeric = (v1 - v2)/(delta1 - delta2); - double err = std::abs((numeric-analytic)/analytic); + double err = mix_deriv_err_func(numeric, analytic); CAPTURE(numeric); CAPTURE(analytic); CHECK(err < 1e-7); @@ -1560,7 +1618,7 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") double v1 = MixtureDerivatives::d2_ndalphardni_dDelta_dTau(rHEOS_plusrho_constT, i, xN_flag), delta1 = rHEOS_plusrho_constT.delta(); double v2 = MixtureDerivatives::d2_ndalphardni_dDelta_dTau(rHEOS_minusrho_constT, i, xN_flag), delta2 = rHEOS_minusrho_constT.delta(); double numeric = (v1 - v2)/(delta1 - delta2); - double err = std::abs((numeric-analytic)/analytic); + double err = mix_deriv_err_func(numeric, analytic); CAPTURE(numeric); CAPTURE(analytic); CHECK(err < 1e-7); @@ -1573,7 +1631,7 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") double v1 = MixtureDerivatives::d_ndalphardni_dDelta(rHEOS_plusrho_constT, i, xN_flag), delta1 = rHEOS_plusrho_constT.delta(); double v2 = MixtureDerivatives::d_ndalphardni_dDelta(rHEOS_minusrho_constT, i, xN_flag), delta2 = rHEOS_minusrho_constT.delta(); double numeric = (v1 - v2)/(delta1 - delta2); - double err = std::abs((numeric-analytic)/analytic); + double err = mix_deriv_err_func(numeric, analytic); CAPTURE(numeric); CAPTURE(analytic); CHECK(err < 1e-7); @@ -1586,7 +1644,7 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") double v1 = MixtureDerivatives::d2_ndalphardni_dDelta2(rHEOS_plusrho_constT, i, xN_flag), delta1 = rHEOS_plusrho_constT.delta(); double v2 = MixtureDerivatives::d2_ndalphardni_dDelta2(rHEOS_minusrho_constT, i, xN_flag), delta2 = rHEOS_minusrho_constT.delta(); double numeric = (v1 - v2)/(delta1 - delta2); - double err = std::abs((numeric-analytic)/analytic); + double err = mix_deriv_err_func(numeric, analytic); CAPTURE(numeric); CAPTURE(analytic); CHECK(err < 1e-7); @@ -1611,7 +1669,7 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") 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); + double err = mix_deriv_err_func(numeric, analytic); CAPTURE(numeric); CAPTURE(analytic); CHECK(err < 1e-7); @@ -1625,7 +1683,7 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") 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); double numeric = (v1 - v2)/(2*dz); - double err = std::abs((numeric-analytic)/analytic); + double err = mix_deriv_err_func(numeric, analytic); CAPTURE(numeric); CAPTURE(analytic); CHECK(err < 1e-8); @@ -1639,7 +1697,7 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") 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); double numeric = (v1 - v2)/(2*dz); - double err = std::abs((numeric-analytic)/analytic); + double err = mix_deriv_err_func(numeric, analytic); CAPTURE(numeric); CAPTURE(analytic); CHECK(err < 1e-8); @@ -1653,7 +1711,7 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") 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); + double err = mix_deriv_err_func(numeric, analytic); CAPTURE(numeric); CAPTURE(analytic); CHECK(err < 1e-7); @@ -1668,14 +1726,7 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") 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; - if (std::abs(analytic) > DBL_EPSILON){ - err = std::abs((numeric-analytic)/analytic); - } - else{ - err = numeric-analytic; - } - + double err = mix_deriv_err_func(numeric, analytic); CAPTURE(numeric); CAPTURE(analytic); CHECK(err < 1e-8); @@ -1689,8 +1740,7 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") 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); 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;} + double err = mix_deriv_err_func(numeric, analytic); CAPTURE(numeric); CAPTURE(analytic); CHECK(err < 1e-8); @@ -1703,8 +1753,7 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") double v1 = rHEOS.Reducing->PSI_rho(rHEOS_pluszj.get_mole_fractions(), i, xN_flag); double v2 = rHEOS.Reducing->PSI_rho(rHEOS_minuszj.get_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; } + double err = mix_deriv_err_func(numeric, analytic); CAPTURE(numeric); CAPTURE(analytic); CHECK(err < 1e-8); @@ -1716,8 +1765,7 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") double v1 = rHEOS.Reducing->PSI_T(rHEOS_pluszj.get_mole_fractions(), i, xN_flag); double v2 = rHEOS.Reducing->PSI_T(rHEOS_minuszj.get_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; } + double err = mix_deriv_err_func(numeric, analytic); CAPTURE(numeric); CAPTURE(analytic); CHECK(err < 1e-8); @@ -1730,7 +1778,7 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") double v1 = MixtureDerivatives::d_ndalphardni_dDelta(rHEOS_pluszj_consttaudelta, i, xN_flag); double v2 = MixtureDerivatives::d_ndalphardni_dDelta(rHEOS_minuszj_consttaudelta, i, xN_flag); double numeric = (v1 - v2)/(2*dz); - double err = std::abs((numeric-analytic)/analytic); + double err = mix_deriv_err_func(numeric, analytic); CAPTURE(numeric); CAPTURE(analytic); CHECK(err < 1e-7); @@ -1743,7 +1791,7 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") double v1 = MixtureDerivatives::d2_ndalphardni_dDelta2(rHEOS_pluszj_consttaudelta, i, xN_flag); double v2 = MixtureDerivatives::d2_ndalphardni_dDelta2(rHEOS_minuszj_consttaudelta, i, xN_flag); double numeric = (v1 - v2)/(2*dz); - double err = std::abs((numeric-analytic)/analytic); + double err = mix_deriv_err_func(numeric, analytic); CAPTURE(numeric); CAPTURE(analytic); CHECK(err < 1e-7); @@ -1756,7 +1804,7 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") double v1 = MixtureDerivatives::d_ndalphardni_dTau(rHEOS_pluszj_consttaudelta, i, xN_flag); double v2 = MixtureDerivatives::d_ndalphardni_dTau(rHEOS_minuszj_consttaudelta, i, xN_flag); double numeric = (v1 - v2)/(2*dz); - double err = std::abs((numeric-analytic)/analytic); + double err = mix_deriv_err_func(numeric, analytic); CAPTURE(numeric); CAPTURE(analytic); CHECK(err < 1e-7); @@ -1768,7 +1816,7 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") double v1 = MixtureDerivatives::d2_ndalphardni_dTau2(rHEOS_pluszj_consttaudelta, i, xN_flag); double v2 = MixtureDerivatives::d2_ndalphardni_dTau2(rHEOS_minuszj_consttaudelta, i, xN_flag); double numeric = (v1 - v2)/(2*dz); - double err = std::abs((numeric-analytic)/analytic); + double err = mix_deriv_err_func(numeric, analytic); CAPTURE(numeric); CAPTURE(analytic); CHECK(err < 1e-7); @@ -1780,7 +1828,7 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") double v1 = MixtureDerivatives::d2_ndalphardni_dDelta_dTau(rHEOS_pluszj_consttaudelta, i, xN_flag); double v2 = MixtureDerivatives::d2_ndalphardni_dDelta_dTau(rHEOS_minuszj_consttaudelta, i, xN_flag); double numeric = (v1 - v2)/(2*dz); - double err = std::abs((numeric-analytic)/analytic); + double err = mix_deriv_err_func(numeric, analytic); CAPTURE(numeric); CAPTURE(analytic); CHECK(err < 1e-7); @@ -1792,7 +1840,7 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") double v1 = MixtureDerivatives::d2alphar_dxi_dDelta(rHEOS_pluszj_consttaudelta, i, xN_flag); double v2 = MixtureDerivatives::d2alphar_dxi_dDelta(rHEOS_minuszj_consttaudelta, i, xN_flag); double numeric = (v1 - v2)/(2*dz); - double err = std::abs((numeric-analytic)/analytic); + double err = mix_deriv_err_func(numeric, analytic); CAPTURE(numeric); CAPTURE(analytic); CHECK(err < 1e-7); @@ -1804,7 +1852,7 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") double v1 = MixtureDerivatives::d3alphar_dxi_dDelta2(rHEOS_pluszj_consttaudelta, i, xN_flag); double v2 = MixtureDerivatives::d3alphar_dxi_dDelta2(rHEOS_minuszj_consttaudelta, i, xN_flag); double numeric = (v1 - v2)/(2*dz); - double err = std::abs((numeric-analytic)/analytic); + double err = mix_deriv_err_func(numeric, analytic); CAPTURE(numeric); CAPTURE(analytic); CHECK(err < 1e-7); @@ -1816,7 +1864,7 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") double v1 = MixtureDerivatives::d3alphar_dxi_dDelta_dTau(rHEOS_pluszj_consttaudelta, i, xN_flag); double v2 = MixtureDerivatives::d3alphar_dxi_dDelta_dTau(rHEOS_minuszj_consttaudelta, i, xN_flag); double numeric = (v1 - v2)/(2*dz); - double err = std::abs((numeric-analytic)/analytic); + double err = mix_deriv_err_func(numeric, analytic); CAPTURE(numeric); CAPTURE(analytic); CHECK(err < 1e-7); @@ -1829,7 +1877,7 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") double v1 = MixtureDerivatives::d2alphar_dxi_dTau(rHEOS_pluszj_consttaudelta, i, xN_flag); double v2 = MixtureDerivatives::d2alphar_dxi_dTau(rHEOS_minuszj_consttaudelta, i, xN_flag); double numeric = (v1 - v2)/(2*dz); - double err = std::abs((numeric-analytic)/analytic); + double err = mix_deriv_err_func(numeric, analytic); CAPTURE(numeric); CAPTURE(analytic); CHECK(err < 1e-7); @@ -1842,7 +1890,7 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") double v1 = MixtureDerivatives::d3alphar_dxi_dTau2(rHEOS_pluszj_consttaudelta, i, xN_flag); double v2 = MixtureDerivatives::d3alphar_dxi_dTau2(rHEOS_minuszj_consttaudelta, i, xN_flag); double numeric = (v1 - v2)/(2*dz); - double err = std::abs((numeric-analytic)/analytic); + double err = mix_deriv_err_func(numeric, analytic); CAPTURE(numeric); CAPTURE(analytic); CHECK(err < 1e-7); @@ -1856,7 +1904,7 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") double v1 = MixtureDerivatives::nd_ndalphardni_dnj__constT_V(rHEOS_plusrho_constT, i, j, xN_flag), delta1 = rHEOS_plusrho_constT.delta(); double v2 = MixtureDerivatives::nd_ndalphardni_dnj__constT_V(rHEOS_minusrho_constT, i, j, xN_flag), delta2 = rHEOS_minusrho_constT.delta(); double numeric = (v1 - v2)/(delta1 - delta2); - double err = std::abs((numeric-analytic)/analytic); + double err = mix_deriv_err_func(numeric, analytic); CAPTURE(numeric); CAPTURE(analytic); CHECK(err < 1e-7); @@ -1869,7 +1917,7 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") double v1 = MixtureDerivatives::ndln_fugacity_i_dnj__constT_V_xi(rHEOS_plusrho_constT, i, j, xN_flag), delta1 = rHEOS_plusrho_constT.delta(); double v2 = MixtureDerivatives::ndln_fugacity_i_dnj__constT_V_xi(rHEOS_minusrho_constT, i, j, xN_flag), delta2 = rHEOS_minusrho_constT.delta(); double numeric = (v1 - v2)/(delta1 - delta2); - double err = std::abs((numeric-analytic)/analytic); + double err = mix_deriv_err_func(numeric, analytic); CAPTURE(numeric); CAPTURE(analytic); CHECK(err < 1e-7); @@ -1882,7 +1930,7 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") double v1 = MixtureDerivatives::ndln_fugacity_i_dnj__constT_V_xi(rHEOS_plusT_constrho, i, j, xN_flag), tau1 = rHEOS_plusT_constrho.tau(); double v2 = MixtureDerivatives::ndln_fugacity_i_dnj__constT_V_xi(rHEOS_minusT_constrho, i, j, xN_flag), tau2 = rHEOS_minusT_constrho.tau(); double numeric = (v1 - v2)/(tau1 - tau2); - double err = std::abs((numeric-analytic)/analytic); + double err = mix_deriv_err_func(numeric, analytic); CAPTURE(numeric); CAPTURE(analytic); CHECK(err < 1e-7); @@ -1896,7 +1944,7 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") double v1 = MixtureDerivatives::nd_ndalphardni_dnj__constT_V(rHEOS_plusT_constrho, i, j, xN_flag), tau1 = rHEOS_plusT_constrho.tau(); double v2 = MixtureDerivatives::nd_ndalphardni_dnj__constT_V(rHEOS_minusT_constrho, i, j, xN_flag), tau2 = rHEOS_minusT_constrho.tau(); double numeric = (v1 - v2)/(tau1 - tau2); - double err = std::abs((numeric-analytic)/analytic); + double err = mix_deriv_err_func(numeric, analytic); CAPTURE(numeric); CAPTURE(analytic); CHECK(err < 1e-7); @@ -1909,7 +1957,7 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") double v1 = MixtureDerivatives::nddeltadni__constT_V_nj(rHEOS_pluszj_consttaudelta, i, xN_flag); double v2 = MixtureDerivatives::nddeltadni__constT_V_nj(rHEOS_minuszj_consttaudelta, i, xN_flag); double numeric = (v1 - v2)/(2*dz); - double err = std::abs((numeric-analytic)/analytic); + double err = mix_deriv_err_func(numeric, analytic); CAPTURE(numeric); CAPTURE(analytic); CHECK(err < 1e-7); @@ -1922,7 +1970,57 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") double v1 = MixtureDerivatives::ndtaudni__constT_V_nj(rHEOS_pluszj_consttaudelta, i, xN_flag); double v2 = MixtureDerivatives::ndtaudni__constT_V_nj(rHEOS_minuszj_consttaudelta, i, xN_flag); double numeric = (v1 - v2)/(2*dz); - double err = std::abs((numeric-analytic)/analytic); + double err = mix_deriv_err_func(numeric, analytic); + CAPTURE(numeric); + CAPTURE(analytic); + CHECK(err < 1e-7); + } + std::ostringstream ss3s; ss3s << "d2_ndtaudni_dxj_dTau__constdelta, i=" << i << ", j=" << j; + SECTION(ss3s.str(), "") + { + double analytic = MixtureDerivatives::d2_ndtaudni_dxj_dTau__constdelta(rHEOS, i, j, xN_flag); + double v1 = MixtureDerivatives::d_ndtaudni_dxj__constdelta_tau(rHEOS_plusT_constrho, i, j, xN_flag), tau1 = rHEOS_plusT_constrho.tau(); + double v2 = MixtureDerivatives::d_ndtaudni_dxj__constdelta_tau(rHEOS_minusT_constrho, i, j, xN_flag), tau2 = rHEOS_minusT_constrho.tau(); + double numeric = (v1 - v2)/(tau1 - tau2); + double err = mix_deriv_err_func(numeric, analytic); + CAPTURE(numeric); + CAPTURE(analytic); + CHECK(err < 1e-7); + } + std::ostringstream ss3t; ss3t << "d_nAij_dX, i=" << i << ", j=" << j; + SECTION(ss3t.str(), "") + { + double analytic = MixtureDerivatives::d_nAij_dX(rHEOS, i, j, xN_flag, CoolProp::iTau); + double v1 = MixtureDerivatives::nAij(rHEOS_plusT_constrho, i, j, xN_flag), tau1 = rHEOS_plusT_constrho.tau(); + double v2 = MixtureDerivatives::nAij(rHEOS_minusT_constrho, i, j, xN_flag), tau2 = rHEOS_minusT_constrho.tau(); + double numeric = (v1 - v2)/(tau1 - tau2); + double err = mix_deriv_err_func(numeric, analytic); + CAPTURE(numeric); + CAPTURE(analytic); + CHECK(err < 1e-7); + } + std::ostringstream ss3o3; ss3o3 << "d_ndln_fugacity_i_dnj_dtau__constdelta_x, i=" << i << ", j=" << j; + SECTION(ss3o3.str(), "") + { + if (i == Ncomp-1){ break; } + double analytic = MixtureDerivatives::d2_ndln_fugacity_i_dnj_dtau2__constdelta_x(rHEOS, i, j, xN_flag); + double v1 = MixtureDerivatives::d_ndln_fugacity_i_dnj_dtau__constdelta_x(rHEOS_plusT_constrho, i, j, xN_flag), tau1 = rHEOS_plusT_constrho.tau(); + double v2 = MixtureDerivatives::d_ndln_fugacity_i_dnj_dtau__constdelta_x(rHEOS_minusT_constrho, i, j, xN_flag), tau2 = rHEOS_minusT_constrho.tau(); + double numeric = (v1 - v2)/(tau1 - tau2); + double err = mix_deriv_err_func(numeric, analytic); + CAPTURE(numeric); + CAPTURE(analytic); + CHECK(err < 1e-7); + } + std::ostringstream ss3o4; ss3o4 << "d2_ndln_fugacity_i_dnj_ddelta_dtau__constx, i=" << i << ", j=" << j; + SECTION(ss3o4.str(), "") + { + if (i == Ncomp-1){ break; } + double analytic = MixtureDerivatives::d2_ndln_fugacity_i_dnj_ddelta_dtau__constx(rHEOS, i, j, xN_flag); + double v1 = MixtureDerivatives::d_ndln_fugacity_i_dnj_ddelta__consttau_x(rHEOS_plusT_constrho, i, j, xN_flag), tau1 = rHEOS_plusT_constrho.tau(); + double v2 = MixtureDerivatives::d_ndln_fugacity_i_dnj_ddelta__consttau_x(rHEOS_minusT_constrho, i, j, xN_flag), tau2 = rHEOS_minusT_constrho.tau(); + double numeric = (v1 - v2)/(tau1 - tau2); + double err = mix_deriv_err_func(numeric, analytic); CAPTURE(numeric); CAPTURE(analytic); CHECK(err < 1e-7); @@ -1944,8 +2042,7 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") double v1 = rHEOS.Reducing->d2Trdxidxj(rHEOS_pluszk.get_mole_fractions(), i, j, xN_flag); double v2 = rHEOS.Reducing->d2Trdxidxj(rHEOS_minuszk.get_mole_fractions(), i, j, 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){ err = 0; } + double err = mix_deriv_err_func(numeric, analytic); CAPTURE(numeric); CAPTURE(analytic); CHECK(err < 1e-7); @@ -1958,7 +2055,7 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") double v1 = rHEOS.Reducing->d2rhormolardxidxj(rHEOS_pluszk.get_mole_fractions(), i, j, xN_flag); double v2 = rHEOS.Reducing->d2rhormolardxidxj(rHEOS_minuszk.get_mole_fractions(), i, j, xN_flag); double numeric = (v1 - v2)/(2*dz); - double err = std::abs((numeric-analytic)/analytic); + double err = mix_deriv_err_func(numeric, analytic); CAPTURE(numeric); CAPTURE(analytic); CHECK(err < 1e-7); @@ -1971,7 +2068,7 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") double v1 = rHEOS.Reducing->d_ndTrdni_dxj__constxi(rHEOS_pluszk.get_mole_fractions(), i, j, xN_flag); double v2 = rHEOS.Reducing->d_ndTrdni_dxj__constxi(rHEOS_minuszk.get_mole_fractions(), i, j, xN_flag); double numeric = (v1 - v2)/(2*dz); - double err = std::abs((numeric-analytic)/analytic); + double err = mix_deriv_err_func(numeric, analytic); CAPTURE(numeric); CAPTURE(analytic); CHECK(err < 1e-7); @@ -1984,7 +2081,7 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") double v1 = rHEOS.Reducing->d_ndrhorbardni_dxj__constxi(rHEOS_pluszk.get_mole_fractions(), i, j, xN_flag); double v2 = rHEOS.Reducing->d_ndrhorbardni_dxj__constxi(rHEOS_minuszk.get_mole_fractions(), i, j, xN_flag); double numeric = (v1 - v2)/(2*dz); - double err = std::abs((numeric-analytic)/analytic); + double err = mix_deriv_err_func(numeric, analytic); CAPTURE(numeric); CAPTURE(analytic); CHECK(err < 1e-7); @@ -1997,7 +2094,7 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") double v1 = rHEOS.Reducing->d_PSI_T_dxj(rHEOS_pluszk.get_mole_fractions(), i, j, xN_flag); double v2 = rHEOS.Reducing->d_PSI_T_dxj(rHEOS_minuszk.get_mole_fractions(), i, j, xN_flag); double numeric = (v1 - v2)/(2*dz); - double err = std::abs((numeric-analytic)/analytic); + double err = mix_deriv_err_func(numeric, analytic); CAPTURE(numeric); CAPTURE(analytic); CHECK(err < 1e-7); @@ -2010,7 +2107,7 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") double v1 = rHEOS.Reducing->d_PSI_rho_dxj(rHEOS_pluszk.get_mole_fractions(), i, j, xN_flag); double v2 = rHEOS.Reducing->d_PSI_rho_dxj(rHEOS_minuszk.get_mole_fractions(), i, j, xN_flag); double numeric = (v1 - v2)/(2*dz); - double err = std::abs((numeric-analytic)/analytic); + double err = mix_deriv_err_func(numeric, analytic); CAPTURE(numeric); CAPTURE(analytic); CHECK(err < 1e-7); @@ -2023,7 +2120,7 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") double v1 = MixtureDerivatives::d_ndalphardni_dxj__constdelta_tau_xi(rHEOS_pluszk_consttaudelta, i, j, xN_flag); double v2 = MixtureDerivatives::d_ndalphardni_dxj__constdelta_tau_xi(rHEOS_minuszk_consttaudelta, i, j, xN_flag); double numeric = (v1 - v2)/(2*dz); - double err = std::abs((numeric-analytic)/analytic); + double err = mix_deriv_err_func(numeric, analytic); CAPTURE(numeric); CAPTURE(analytic); CHECK(err < 1e-7); @@ -2036,7 +2133,7 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") double v1 = MixtureDerivatives::d2_ndalphardni_dxj_dDelta__consttau_xi(rHEOS_pluszk_consttaudelta, i, j, xN_flag); double v2 = MixtureDerivatives::d2_ndalphardni_dxj_dDelta__consttau_xi(rHEOS_minuszk_consttaudelta, i, j, xN_flag); double numeric = (v1 - v2)/(2*dz); - double err = std::abs((numeric-analytic)/analytic); + double err = mix_deriv_err_func(numeric, analytic); CAPTURE(numeric); CAPTURE(analytic); CHECK(err < 1e-7); @@ -2049,7 +2146,7 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") double v1 = MixtureDerivatives::d2_ndalphardni_dxj_dTau__constdelta_xi(rHEOS_pluszk_consttaudelta, i, j, xN_flag); double v2 = MixtureDerivatives::d2_ndalphardni_dxj_dTau__constdelta_xi(rHEOS_minuszk_consttaudelta, i, j, xN_flag); double numeric = (v1 - v2)/(2*dz); - double err = std::abs((numeric-analytic)/analytic); + double err = mix_deriv_err_func(numeric, analytic); CAPTURE(numeric); CAPTURE(analytic); CHECK(err < 1e-7); @@ -2063,14 +2160,7 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") double v2 = MixtureDerivatives::d2alphardxidxj(rHEOS_minuszk_consttaudelta, i, j, xN_flag); double numeric = (v1 - v2)/(2*dz); if (std::abs(numeric) < DBL_EPSILON && std::abs(analytic) < DBL_EPSILON){ break; } - double err; - if (std::abs(analytic) > DBL_EPSILON){ - err = std::abs((numeric-analytic)/analytic); - } - else{ - err = numeric-analytic; - } - + double err = mix_deriv_err_func(numeric, analytic); CAPTURE(numeric); CAPTURE(analytic); CHECK(err < 1e-8); @@ -2085,14 +2175,7 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") double v2 = MixtureDerivatives::nd_ndalphardni_dnj__constT_V(rHEOS_minuszk_consttaudelta, i, j, xN_flag); double numeric = (v1 - v2)/(2*dz); if (std::abs(numeric) < DBL_EPSILON && std::abs(analytic) < DBL_EPSILON){ break; } - double err; - if (std::abs(analytic) > DBL_EPSILON){ - err = std::abs((numeric-analytic)/analytic); - } - else{ - err = numeric-analytic; - } - + double err = mix_deriv_err_func(numeric, analytic); CAPTURE(numeric); CAPTURE(analytic); CHECK(err < 1e-8); @@ -2106,18 +2189,61 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") double v2 = MixtureDerivatives::ndln_fugacity_i_dnj__constT_V_xi(rHEOS_minuszk_consttaudelta, i, j, xN_flag); double numeric = (v1 - v2)/(2*dz); if (std::abs(numeric) < DBL_EPSILON && std::abs(analytic) < DBL_EPSILON){ break; } - double err; - if (std::abs(analytic) > DBL_EPSILON){ - err = std::abs((numeric-analytic)/analytic); - } - else{ - err = numeric-analytic; - } - + double err = mix_deriv_err_func(numeric, analytic); CAPTURE(numeric); CAPTURE(analytic); CHECK(err < 1e-8); } + std::ostringstream ss3o4; ss3o4 << "d2_ndln_fugacity_i_dnj_dxk_dTau__constdelta, i=" << i << ", j=" << j << ", k=" << k; + SECTION(ss3o4.str(), "") + { + if (i == Ncomp-1){ break; } + double analytic = MixtureDerivatives::d2_ndln_fugacity_i_dnj_dxk_dTau__constdelta(rHEOS, i, j, k, xN_flag); + double v1 = MixtureDerivatives::d_ndln_fugacity_i_dnj_ddxk__consttau_delta(rHEOS_plusT_constrho, i, j, k, xN_flag), tau1 = rHEOS_plusT_constrho.tau(); + double v2 = MixtureDerivatives::d_ndln_fugacity_i_dnj_ddxk__consttau_delta(rHEOS_minusT_constrho, i, j, k, xN_flag), tau2 = rHEOS_minusT_constrho.tau(); + double numeric = (v1 - v2)/(tau1 - tau2); + double err = mix_deriv_err_func(numeric, analytic); + CAPTURE(numeric); + CAPTURE(analytic); + CHECK(err < 1e-7); + } + std::ostringstream ss3oo4; ss3oo4 << "d2_ndln_fugacity_i_dnj_dxk_dDelta__consttau, i=" << i << ", j=" << j << ", k=" << k; + SECTION(ss3oo4.str(), "") + { + if (i == Ncomp-1){ break; } + double analytic = MixtureDerivatives::d2_ndln_fugacity_i_dnj_dxk_dDelta__consttau(rHEOS, i, j, k, xN_flag); + double v1 = MixtureDerivatives::d_ndln_fugacity_i_dnj_ddxk__consttau_delta(rHEOS_plusrho_constT, i, j, k, xN_flag), delta1 = rHEOS_plusrho_constT.delta(); + double v2 = MixtureDerivatives::d_ndln_fugacity_i_dnj_ddxk__consttau_delta(rHEOS_minusrho_constT, i, j, k, xN_flag), delta2 = rHEOS_minusrho_constT.delta(); + double numeric = (v1 - v2)/(delta1 - delta2); + double err = mix_deriv_err_func(numeric, analytic); + CAPTURE(numeric); + CAPTURE(analytic); + CHECK(err < 1e-7); + } + std::ostringstream ss3o5; ss3o5 << "d_n2Aijk_dTau, i=" << i << ", j=" << j; + SECTION(ss3o5.str(), "") + { + double analytic = MixtureDerivatives::d_n2Aijk_dX(rHEOS, i, j, k, xN_flag, CoolProp::iTau); + double v1 = MixtureDerivatives::n2Aijk(rHEOS_plusT_constrho, i, j, k, xN_flag), tau1 = rHEOS_plusT_constrho.tau(); + double v2 = MixtureDerivatives::n2Aijk(rHEOS_minusT_constrho, i, j, k, xN_flag), tau2 = rHEOS_minusT_constrho.tau(); + double numeric = (v1 - v2)/(tau1 - tau2); + double err = mix_deriv_err_func(numeric, analytic); + CAPTURE(numeric); + CAPTURE(analytic); + CHECK(err < 1e-7); + } + std::ostringstream ss3o6; ss3o6 << "d_n2Aijk_dDelta, i=" << i << ", j=" << j; + SECTION(ss3o6.str(), "") + { + double analytic = MixtureDerivatives::d_n2Aijk_dX(rHEOS, i, j, k, xN_flag, CoolProp::iDelta); + double v1 = MixtureDerivatives::n2Aijk(rHEOS_plusrho_constT, i, j, k, xN_flag), delta1 = rHEOS_plusrho_constT.delta(); + double v2 = MixtureDerivatives::n2Aijk(rHEOS_minusrho_constT, i, j, k, xN_flag), delta2 = rHEOS_minusrho_constT.delta(); + double numeric = (v1 - v2)/(delta1 - delta2); + double err = mix_deriv_err_func(numeric, analytic); + CAPTURE(numeric); + CAPTURE(analytic); + CHECK(err < 1e-7); + } } } } From 78d799ee75cc6123c4af20b2d4da8fd01a369e06 Mon Sep 17 00:00:00 2001 From: Ian Bell Date: Fri, 14 Aug 2015 18:24:29 -0600 Subject: [PATCH 5/8] And the changes to the header --- src/Backends/Helmholtz/MixtureDerivatives.h | 15 ++++++++++++++- 1 file changed, 14 insertions(+), 1 deletion(-) diff --git a/src/Backends/Helmholtz/MixtureDerivatives.h b/src/Backends/Helmholtz/MixtureDerivatives.h index d3caafb7..4d7721d0 100644 --- a/src/Backends/Helmholtz/MixtureDerivatives.h +++ b/src/Backends/Helmholtz/MixtureDerivatives.h @@ -208,7 +208,20 @@ class MixtureDerivatives{ return 1/RT*summer + 1/(HEOS.gas_constant()*HEOS.T_reducing())*nd_ndln_fugacity_i_dnj_dnk__constT_V_xi(HEOS,i,j,k,xN_flag) - d_nAij_dX(HEOS, i, j, xN_flag, WRT); } else if (WRT== iDelta){ - return -1000; + double summer = 0; + summer += d2_ndln_fugacity_i_dnj_ddelta_dtau__constx(HEOS, i, j, xN_flag)*ndtaudni__constT_V_nj(HEOS, k, xN_flag); + summer += d2_ndln_fugacity_i_dnj_ddelta2__consttau_x(HEOS, i, j, xN_flag)*nddeltadni__constT_V_nj(HEOS, k, xN_flag); + summer += d_ndln_fugacity_i_dnj_ddelta__consttau_x(HEOS, i, j, xN_flag)*d_nddeltadni_dDelta(HEOS, k, xN_flag); + summer += d2_ndln_fugacity_i_dnj_dxk_dDelta__consttau(HEOS, i, j, k, xN_flag); + std::size_t mmax = HEOS.mole_fractions.size(); + if (xN_flag == XN_DEPENDENT){ mmax--; } + for (std::size_t m = 0; m < mmax; ++m){ + summer -= HEOS.mole_fractions[m]*d2_ndln_fugacity_i_dnj_dxk_dDelta__consttau(HEOS, i, j, m, xN_flag); + } + return 1/RT*summer - d_nAij_dX(HEOS, i, j, xN_flag, iDelta); + } + else{ + return _HUGE; } } static Eigen::MatrixXd Lstar(HelmholtzEOSMixtureBackend &HEOS, x_N_dependency_flag xN_flag){ From ffbdbb9ce9528466375dfd7c4e685630830d298c Mon Sep 17 00:00:00 2001 From: Ian Bell Date: Fri, 14 Aug 2015 22:06:41 -0600 Subject: [PATCH 6/8] Small updates to calculation of critical points; analytic derivatives not fully working --- .../Helmholtz/HelmholtzEOSMixtureBackend.cpp | 21 ++- src/Backends/Helmholtz/MixtureDerivatives.cpp | 175 ++++++++++++------ src/Backends/Helmholtz/MixtureDerivatives.h | 26 +-- 3 files changed, 139 insertions(+), 83 deletions(-) diff --git a/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp b/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp index b752443d..837fc22a 100644 --- a/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp +++ b/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp @@ -2915,9 +2915,9 @@ void HelmholtzEOSMixtureBackend::calc_critical_point(double rho0, double T0) }; std::vector > Jacobian(const std::vector &x) { - double epsilon; + std::size_t N = x.size(); - std::vector r, xp; + std::vector rplus, rminus, xp; std::vector > J(N, std::vector(N, 0)), JJ(N, std::vector(N, 0)); Eigen::MatrixXd J0(N, N), adjL = adjugate(MixtureDerivatives::Lstar(HEOS, XN_INDEPENDENT)), adjM = adjugate(MixtureDerivatives::Mstar(HEOS, XN_INDEPENDENT)), @@ -2931,26 +2931,32 @@ void HelmholtzEOSMixtureBackend::calc_critical_point(double rho0, double T0) JJ[1][0] = (adjM*dMdTau).trace(); JJ[1][1] = (adjM*dMdDelta).trace(); +// return JJ; + //std::cout << J0 << std::endl; std::cout << JJ[0][0] << " " << JJ[0][1] << std::endl; std::cout << JJ[1][0] << " " << JJ[1][1] << std::endl; - std::vector r0 = call(x); + double epsilon; + // Build the Jacobian by column for (std::size_t i = 0; i < N; ++i) { xp = x; - epsilon = 0.00001; + epsilon = 0.0001; xp[i] += epsilon; - r = call(xp); + rplus = call(xp); + xp[i] -= 2*epsilon; + rminus = call(xp); for (std::size_t j = 0; j < N; ++j) { - J[j][i] = (r[j]-r0[j])/epsilon; + J[j][i] = (rplus[j]-rminus[j])/(2*epsilon); } } std::cout << J[0][0] << " " << J[0][1] << std::endl; std::cout << J[1][0] << " " << J[1][1] << std::endl; return J; + }; }; Resid resid(*this); @@ -2976,13 +2982,11 @@ public: return L1; } double deriv(double tau){ - 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); return (adjL*dLdTau).trace(); }; double second_deriv(double tau){ - std::size_t N = HEOS.get_mole_fractions().size(); Eigen::MatrixXd Lstar = MixtureDerivatives::Lstar(HEOS, XN_INDEPENDENT), dLstardTau = MixtureDerivatives::dLstar_dX(HEOS, XN_INDEPENDENT, iTau), d2LstardTau2 = MixtureDerivatives::d2Lstar_dX2(HEOS, XN_INDEPENDENT, iTau, iTau), @@ -3039,7 +3043,6 @@ public: @param theta The angle */ double deriv(double theta){ - std::size_t N = HEOS.get_mole_fractions().size(); double dL1_dtau = (adjLstar*dLstardTau).trace(), dL1_ddelta = (adjLstar*dLstardDelta).trace(); return -R_tau*sin(theta)*dL1_dtau + R_delta*cos(theta)*dL1_ddelta; }; diff --git a/src/Backends/Helmholtz/MixtureDerivatives.cpp b/src/Backends/Helmholtz/MixtureDerivatives.cpp index 62ce9094..bde093ff 100644 --- a/src/Backends/Helmholtz/MixtureDerivatives.cpp +++ b/src/Backends/Helmholtz/MixtureDerivatives.cpp @@ -1014,7 +1014,7 @@ double mix_deriv_err_func(double numeric, double analytic) return std::abs(numeric - analytic); } else{ - return (numeric/analytic-1); + return std::abs(numeric/analytic-1); } } @@ -1174,8 +1174,96 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") const std::vector &z = rHEOS.get_mole_fractions(); - std::ostringstream ss3o55; ss3o55 << "dLstar_Tau"; - SECTION(ss3o55.str(), "") + SECTION("adj(Mstar)", "") + { + if (Ncomp == 2){ + Eigen::MatrixXd Mstar = MixtureDerivatives::Mstar(rHEOS, xN_flag); + Eigen::MatrixXd analytic = adjugate(Mstar); + Eigen::MatrixXd numeric(2,2); + numeric << Mstar(1,1), -Mstar(1,0), -Mstar(0,1), Mstar(0,0); + double err = ((analytic-numeric).array()/analytic.array()).cwiseAbs().sum()/Ncomp/Ncomp; + CAPTURE(numeric); + CAPTURE(analytic); + CAPTURE(Mstar); + CHECK(err < 1e-8); + } + } + SECTION("d(adj(Lstar))/dDelta", "") + { + Eigen::MatrixXd Lstar = MixtureDerivatives::Lstar(rHEOS, xN_flag); + Eigen::MatrixXd dLstar_dDelta = MixtureDerivatives::dLstar_dX(rHEOS, xN_flag, CoolProp::iDelta); + Eigen::MatrixXd analytic = adjugate_derivative(Lstar, dLstar_dDelta); + + Eigen::MatrixXd adj_Lstar_plus = adjugate(MixtureDerivatives::Lstar(rHEOS_plusrho_constT, xN_flag)); double delta1 = rHEOS_plusrho_constT.delta(); + Eigen::MatrixXd adj_Lstar_minus = adjugate(MixtureDerivatives::Lstar(rHEOS_minusrho_constT, xN_flag)); double delta2 = rHEOS_minusrho_constT.delta(); + + Eigen::MatrixXd numeric = (adj_Lstar_plus - adj_Lstar_minus)/(delta1-delta2); + double err = ((analytic-numeric).array()/analytic.array()).cwiseAbs().sum()/Ncomp/Ncomp; + CAPTURE(numeric); + CAPTURE(analytic); + CHECK(err < 1e-8); + } + SECTION("d(M1)/dTau", "") + { + Eigen::MatrixXd Mstar = MixtureDerivatives::Mstar(rHEOS, xN_flag); + Eigen::MatrixXd dMstar_dTau = MixtureDerivatives::dMstar_dX(rHEOS, xN_flag, CoolProp::iTau); + double analytic = (adjugate(Mstar)*dMstar_dTau).trace(); + + double detMstar_plus = MixtureDerivatives::Mstar(rHEOS_plusT_constrho, xN_flag).determinant(); double tau1 = rHEOS_plusT_constrho.tau(); + double detMstar_minus = MixtureDerivatives::Mstar(rHEOS_minusT_constrho, xN_flag).determinant(); double tau2 = rHEOS_minusT_constrho.tau(); + + double numeric = (detMstar_plus - detMstar_minus)/(tau1-tau2); + double err = mix_deriv_err_func(numeric, analytic); + CAPTURE(numeric); + CAPTURE(analytic); + CHECK(err < 1e-8); + } + SECTION("d(M1)/dDelta", "") + { + Eigen::MatrixXd Mstar = MixtureDerivatives::Mstar(rHEOS, xN_flag); + Eigen::MatrixXd dMstar_dDelta = MixtureDerivatives::dMstar_dX(rHEOS, xN_flag, CoolProp::iDelta); + double analytic = (adjugate(Mstar)*dMstar_dDelta).trace(); + + double detMstar_plus = MixtureDerivatives::Mstar(rHEOS_plusrho_constT, xN_flag).determinant(); double delta1 = rHEOS_plusrho_constT.delta(); + double detMstar_minus = MixtureDerivatives::Mstar(rHEOS_minusrho_constT, xN_flag).determinant(); double delta2 = rHEOS_minusrho_constT.delta(); + + double numeric = (detMstar_plus - detMstar_minus)/(delta1-delta2); + double err = mix_deriv_err_func(numeric, analytic); + CAPTURE(numeric); + CAPTURE(analytic); + CHECK(err < -1e-8); + } + SECTION("d(L1)/dDelta", "") + { + Eigen::MatrixXd Lstar = MixtureDerivatives::Lstar(rHEOS, xN_flag); + Eigen::MatrixXd dLstar_dDelta = MixtureDerivatives::dLstar_dX(rHEOS, xN_flag, CoolProp::iDelta); + double analytic = (adjugate(Lstar)*dLstar_dDelta).trace(); + + double detLstar_plus = MixtureDerivatives::Lstar(rHEOS_plusrho_constT, xN_flag).determinant(); double delta1 = rHEOS_plusrho_constT.delta(); + double detLstar_minus = MixtureDerivatives::Lstar(rHEOS_minusrho_constT, xN_flag).determinant(); double delta2 = rHEOS_minusrho_constT.delta(); + + double numeric = (detLstar_plus - detLstar_minus)/(delta1-delta2); + double err = mix_deriv_err_func(numeric, analytic); + CAPTURE(numeric); + CAPTURE(analytic); + CHECK(err < -1e-8); + } + + SECTION("adj(Lstar)", "") + { + if (Ncomp == 2){ + Eigen::MatrixXd Lstar = MixtureDerivatives::Lstar(rHEOS, xN_flag); + Eigen::MatrixXd analytic = adjugate(Lstar); + Eigen::MatrixXd numeric(2,2); + numeric << Lstar(1,1), -Lstar(1,0), -Lstar(0,1), Lstar(0,0); + double err = ((analytic-numeric).array()/analytic.array()).cwiseAbs().sum()/Ncomp/Ncomp; + CAPTURE(numeric); + CAPTURE(analytic); + CAPTURE(Lstar); + CHECK(err < 1e-8); + } + } + SECTION("dLstar_Tau", "") { Eigen::MatrixXd analytic = MixtureDerivatives::dLstar_dX(rHEOS, xN_flag, CoolProp::iTau); Eigen::MatrixXd m1 = MixtureDerivatives::Lstar(rHEOS_plusT_constrho, xN_flag); double tau1 = rHEOS_plusT_constrho.tau(); @@ -1183,11 +1271,11 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") Eigen::MatrixXd numeric = (m1 - m2)/(tau1 - tau2); double err = ((analytic-numeric).array()/analytic.array()).cwiseAbs().sum()/Ncomp/Ncomp; CAPTURE(numeric); - CAPTURE(analytic); - CHECK(err < 1e-12); + CAPTURE(analytic); + CHECK(err < 1e-8); } - std::ostringstream ss3o5; ss3o5 << "dLstar_dDelta"; - SECTION(ss3o5.str(), "") { + SECTION("dLstar_dDelta", "") + { Eigen::MatrixXd analytic = MixtureDerivatives::dLstar_dX(rHEOS, xN_flag, CoolProp::iDelta); Eigen::MatrixXd m1 = MixtureDerivatives::Lstar(rHEOS_plusrho_constT, xN_flag); double delta1 = rHEOS_plusrho_constT.delta(); Eigen::MatrixXd m2 = MixtureDerivatives::Lstar(rHEOS_minusrho_constT, xN_flag); double delta2 = rHEOS_minusrho_constT.delta(); @@ -1195,10 +1283,9 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") double err = ((analytic-numeric).array()/analytic.array()).cwiseAbs().sum()/Ncomp/Ncomp; CAPTURE(numeric); CAPTURE(analytic); - CHECK(err < 1e-12); + CHECK(err < 1e-8); } - std::ostringstream ss3o6; ss3o6 << "dMstar_dTau"; - SECTION(ss3o6.str(), "") + SECTION("dMstar_dTau", "") { Eigen::MatrixXd analytic = MixtureDerivatives::dMstar_dX(rHEOS, xN_flag, CoolProp::iTau); Eigen::MatrixXd m1 = MixtureDerivatives::Mstar(rHEOS_plusT_constrho, xN_flag); double tau1 = rHEOS_plusT_constrho.tau(); @@ -1207,7 +1294,7 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") double err = ((analytic-numeric).array()/analytic.array()).cwiseAbs().sum()/Ncomp/Ncomp; CAPTURE(numeric); CAPTURE(analytic); - CHECK(err < 1e-12); + CHECK(err < 1e-8); } std::ostringstream ss3o7; ss3o7 << "dMstar_dDelta"; SECTION(ss3o7.str(), "") @@ -1219,7 +1306,7 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") double err = ((analytic-numeric).array()/analytic.array()).cwiseAbs().sum()/Ncomp/Ncomp; CAPTURE(numeric); CAPTURE(analytic); - CHECK(err < 1e-12); + CHECK(err < 1e-8); } // These ones only require the i index @@ -1335,7 +1422,6 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") ss3a << "d2alphar_dxi_dDelta, i=" << i; SECTION(ss3a.str(), "") { - 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(); @@ -1347,7 +1433,6 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") ss4a << "d2alphar_dxi_dTau, i=" << i; SECTION(ss4a.str(), "") { - 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(); @@ -1370,7 +1455,6 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") ss5 << "dpdxj__constT_V_xi, i=" << i; SECTION(ss5.str(), "") { - 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(); @@ -1384,7 +1468,6 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") ss5a << "dtaudxj__constT_V_xi, i=" << i; SECTION(ss5a.str(), "") { - 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(); @@ -1398,7 +1481,6 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") ss5b << "ddeltadxj__constT_V_xi, i=" << i; SECTION(ss5b.str(), "") { - 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(); @@ -1412,7 +1494,6 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") ss6 << "d_dalpharddelta_dxj__constT_V_xi, i=" << i; SECTION(ss6.str(), "") { - 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(); @@ -1426,7 +1507,6 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") ss7 << "dTrdxi__constxj, i=" << i; SECTION(ss7.str(), "") { - 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()); @@ -1440,7 +1520,6 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") ss8 << "drhormolardxi__constxj, i=" << i; SECTION(ss8.str(), "") { - 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()); @@ -1468,7 +1547,6 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") ss3d << "dalphar_dxi, i=" << i; SECTION(ss3d.str(), "") { - if (i == Ncomp-1){ break; } double analytic = MixtureDerivatives::dalphar_dxi(rHEOS, i, xN_flag); double v1 = rHEOS_pluszi_consttaudelta.alphar(); double v2 = rHEOS_minuszi_consttaudelta.alphar(); @@ -1481,7 +1559,6 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") std::ostringstream ss3e; ss3e << "d3alphar_dxi_dDelta2, i=" << i; SECTION(ss3e.str(), "") { - if (i == Ncomp-1){ break; } double analytic = MixtureDerivatives::d3alphar_dxi_dDelta2(rHEOS, i, xN_flag); double v1 = MixtureDerivatives::d2alphar_dxi_dDelta(rHEOS_plusrho_constT, i, xN_flag), delta1 = rHEOS_plusrho_constT.delta(); double v2 = MixtureDerivatives::d2alphar_dxi_dDelta(rHEOS_minusrho_constT, i, xN_flag), delta2 = rHEOS_minusrho_constT.delta(); @@ -1494,7 +1571,6 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") std::ostringstream ss3f; ss3f << "d3alphar_dxi_dTau2, i=" << i; SECTION(ss3f.str(), "") { - if (i == Ncomp-1){ break; } double analytic = MixtureDerivatives::d3alphar_dxi_dTau2(rHEOS, i, xN_flag); double v1 = MixtureDerivatives::d2alphar_dxi_dTau(rHEOS_plusT_constrho, i, xN_flag), tau1 = rHEOS_plusT_constrho.tau(); double v2 = MixtureDerivatives::d2alphar_dxi_dTau(rHEOS_minusT_constrho, i, xN_flag), tau2 = rHEOS_minusT_constrho.tau(); @@ -1507,7 +1583,6 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") std::ostringstream ss3ff; ss3ff << "d4alphar_dxi_dTau3, i=" << i; SECTION(ss3ff.str(), "") { - if (i == Ncomp-1){ break; } double analytic = MixtureDerivatives::d4alphar_dxi_dTau3(rHEOS, i, xN_flag); double v1 = MixtureDerivatives::d3alphar_dxi_dTau2(rHEOS_plusT_constrho, i, xN_flag), tau1 = rHEOS_plusT_constrho.tau(); double v2 = MixtureDerivatives::d3alphar_dxi_dTau2(rHEOS_minusT_constrho, i, xN_flag), tau2 = rHEOS_minusT_constrho.tau(); @@ -1521,7 +1596,6 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") std::ostringstream ss3g; ss3g << "d3alphar_dxi_dDelta_Tau, i=" << i; SECTION(ss3g.str(), "") { - if (i == Ncomp-1){ break; } double analytic = MixtureDerivatives::d3alphar_dxi_dDelta_dTau(rHEOS, i, xN_flag); double v1 = MixtureDerivatives::d2alphar_dxi_dDelta(rHEOS_plusT_constrho, i, xN_flag), tau1 = rHEOS_plusT_constrho.tau(); double v2 = MixtureDerivatives::d2alphar_dxi_dDelta(rHEOS_minusT_constrho, i, xN_flag), tau2 = rHEOS_minusT_constrho.tau(); @@ -1534,7 +1608,6 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") std::ostringstream ss3gg; ss3gg << "d4alphar_dxi_dDelta_Tau2, i=" << i; SECTION(ss3gg.str(), "") { - if (i == Ncomp-1){ break; } double analytic = MixtureDerivatives::d4alphar_dxi_dDelta_dTau2(rHEOS, i, xN_flag); double v1 = MixtureDerivatives::d3alphar_dxi_dDelta_dTau(rHEOS_plusT_constrho, i, xN_flag), tau1 = rHEOS_plusT_constrho.tau(); double v2 = MixtureDerivatives::d3alphar_dxi_dDelta_dTau(rHEOS_minusT_constrho, i, xN_flag), tau2 = rHEOS_minusT_constrho.tau(); @@ -1547,7 +1620,6 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") std::ostringstream ss3ggg; ss3ggg << "d4alphar_dxi_dDelta2_Tau, i=" << i; SECTION(ss3ggg.str(), "") { - if (i == Ncomp-1){ break; } double analytic = MixtureDerivatives::d4alphar_dxi_dDelta2_dTau(rHEOS, i, xN_flag); double v1 = MixtureDerivatives::d3alphar_dxi_dDelta2(rHEOS_plusT_constrho, i, xN_flag), tau1 = rHEOS_plusT_constrho.tau(); double v2 = MixtureDerivatives::d3alphar_dxi_dDelta2(rHEOS_minusT_constrho, i, xN_flag), tau2 = rHEOS_minusT_constrho.tau(); @@ -1561,7 +1633,6 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") std::ostringstream ss3h; ss3h << "d2_ndalphardni_dTau2, i=" << i; SECTION(ss3h.str(), "") { - if (i == Ncomp-1){ break; } double analytic = MixtureDerivatives::d2_ndalphardni_dTau2(rHEOS, i, xN_flag); double v1 = MixtureDerivatives::d_ndalphardni_dTau(rHEOS_plusT_constrho, i, xN_flag), tau1 = rHEOS_plusT_constrho.tau(); double v2 = MixtureDerivatives::d_ndalphardni_dTau(rHEOS_minusT_constrho, i, xN_flag), tau2 = rHEOS_minusT_constrho.tau(); @@ -1574,7 +1645,6 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") std::ostringstream ss3hh; ss3hh << "d3_ndalphardni_dTau3, i=" << i; SECTION(ss3hh.str(), "") { - if (i == Ncomp-1){ break; } double analytic = MixtureDerivatives::d3_ndalphardni_dTau3(rHEOS, i, xN_flag); double v1 = MixtureDerivatives::d2_ndalphardni_dTau2(rHEOS_plusT_constrho, i, xN_flag), tau1 = rHEOS_plusT_constrho.tau(); double v2 = MixtureDerivatives::d2_ndalphardni_dTau2(rHEOS_minusT_constrho, i, xN_flag), tau2 = rHEOS_minusT_constrho.tau(); @@ -1587,7 +1657,6 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") std::ostringstream ss3i; ss3i << "d2_ndalphardni_dDelta_dTau, i=" << i; SECTION(ss3i.str(), "") { - if (i == Ncomp-1){ break; } double analytic = MixtureDerivatives::d2_ndalphardni_dDelta_dTau(rHEOS, i, xN_flag); double v1 = MixtureDerivatives::d_ndalphardni_dTau(rHEOS_plusrho_constT, i, xN_flag), delta1 = rHEOS_plusrho_constT.delta(); double v2 = MixtureDerivatives::d_ndalphardni_dTau(rHEOS_minusrho_constT, i, xN_flag), delta2 = rHEOS_minusrho_constT.delta(); @@ -1600,7 +1669,6 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") std::ostringstream ss3ii; ss3ii << "d3_ndalphardni_dDelta_dTau2, i=" << i; SECTION(ss3ii.str(), "") { - if (i == Ncomp-1){ break; } double analytic = MixtureDerivatives::d3_ndalphardni_dDelta_dTau2(rHEOS, i, xN_flag); double v1 = MixtureDerivatives::d2_ndalphardni_dTau2(rHEOS_plusrho_constT, i, xN_flag), delta1 = rHEOS_plusrho_constT.delta(); double v2 = MixtureDerivatives::d2_ndalphardni_dTau2(rHEOS_minusrho_constT, i, xN_flag), delta2 = rHEOS_minusrho_constT.delta(); @@ -1613,7 +1681,6 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") std::ostringstream ss3iii; ss3iii << "d3_ndalphardni_dDelta2_dTau, i=" << i; SECTION(ss3iii.str(), "") { - if (i == Ncomp-1){ break; } double analytic = MixtureDerivatives::d3_ndalphardni_dDelta2_dTau(rHEOS, i, xN_flag); double v1 = MixtureDerivatives::d2_ndalphardni_dDelta_dTau(rHEOS_plusrho_constT, i, xN_flag), delta1 = rHEOS_plusrho_constT.delta(); double v2 = MixtureDerivatives::d2_ndalphardni_dDelta_dTau(rHEOS_minusrho_constT, i, xN_flag), delta2 = rHEOS_minusrho_constT.delta(); @@ -1626,7 +1693,6 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") std::ostringstream ss3j; ss3j << "d2_ndalphardni_dDelta2, i=" << i; SECTION(ss3j.str(), "") { - if (i == Ncomp-1){ break; } double analytic = MixtureDerivatives::d2_ndalphardni_dDelta2(rHEOS, i, xN_flag); double v1 = MixtureDerivatives::d_ndalphardni_dDelta(rHEOS_plusrho_constT, i, xN_flag), delta1 = rHEOS_plusrho_constT.delta(); double v2 = MixtureDerivatives::d_ndalphardni_dDelta(rHEOS_minusrho_constT, i, xN_flag), delta2 = rHEOS_minusrho_constT.delta(); @@ -1639,7 +1705,6 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") std::ostringstream ss3k; ss3k << "d3_ndalphardni_dDelta3, i=" << i; SECTION(ss3k.str(), "") { - if (i == Ncomp-1){ break; } double analytic = MixtureDerivatives::d3_ndalphardni_dDelta3(rHEOS, i, xN_flag); double v1 = MixtureDerivatives::d2_ndalphardni_dDelta2(rHEOS_plusrho_constT, i, xN_flag), delta1 = rHEOS_plusrho_constT.delta(); double v2 = MixtureDerivatives::d2_ndalphardni_dDelta2(rHEOS_minusrho_constT, i, xN_flag), delta2 = rHEOS_minusrho_constT.delta(); @@ -1664,7 +1729,6 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") SECTION(ss1a.str(), "") { if (xN_flag == XN_INDEPENDENT){continue;} - 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)); @@ -1678,7 +1742,6 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") ss2 << "d_ndTrdni_dxj, i=" << i << ", j=" << j; SECTION(ss2.str(), "") { - 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); @@ -1692,7 +1755,6 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") ss4 << "d_ndrhomolarrdni_dxj, i=" << i << ", j=" << j; SECTION(ss4.str(), "") { - 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); @@ -1706,7 +1768,6 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") ss3 << "d_ndalphardni_dxj__constT_V_xi, i=" << i << ", j=" << j; SECTION(ss3.str(), "") { - 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); @@ -1720,7 +1781,6 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") ss3a << "d2alphardxidxj, i=" << i << ", j=" << j; SECTION(ss3a.str(), "") { - if (j == Ncomp-1){ break; } double analytic = MixtureDerivatives::d2alphardxidxj(rHEOS,i,j,xN_flag); double v1 = MixtureDerivatives::dalphar_dxi(rHEOS_pluszj_consttaudelta, i, xN_flag); double v2 = MixtureDerivatives::dalphar_dxi(rHEOS_minuszj_consttaudelta, i, xN_flag); @@ -1735,7 +1795,6 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") ss3b << "d2Trdxidxj, i=" << i << ", j=" << j; SECTION(ss3b.str(), "") { - 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); @@ -1770,7 +1829,7 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") CAPTURE(analytic); CHECK(err < 1e-8); } - std::ostringstream ss3k; ss3k << "d2_ndalphardni_dxj_dDelta, i=" << i; + std::ostringstream ss3k; ss3k << "d2_ndalphardni_dxj_dDelta, i=" << i << ", j=" << j; SECTION(ss3k.str(), "") { if (i == Ncomp-1){ break; } @@ -1783,7 +1842,7 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") CAPTURE(analytic); CHECK(err < 1e-7); } - std::ostringstream ss3kk; ss3kk << "d3_ndalphardni_dxj_dDelta2__consttau_xi, i=" << i; + std::ostringstream ss3kk; ss3kk << "d3_ndalphardni_dxj_dDelta2__consttau_xi, i=" << i << ", j=" << j; SECTION(ss3kk.str(), "") { if (i == Ncomp-1){ break; } @@ -1796,7 +1855,7 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") CAPTURE(analytic); CHECK(err < 1e-7); } - std::ostringstream ss3l; ss3l << "d2_ndalphardni_dxj_dTau, i=" << i; + std::ostringstream ss3l; ss3l << "d2_ndalphardni_dxj_dTau, i=" << i << ", j=" << j; SECTION(ss3l.str(), "") { if (i == Ncomp-1){ break; } @@ -1809,7 +1868,7 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") CAPTURE(analytic); CHECK(err < 1e-7); } - std::ostringstream ss3ll; ss3ll << "d3_ndalphardni_dxj_dTau2__constdelta_xi, i=" << i; + std::ostringstream ss3ll; ss3ll << "d3_ndalphardni_dxj_dTau2__constdelta_xi, i=" << i << ", j=" << j; SECTION(ss3ll.str(), "") { double analytic = MixtureDerivatives::d3_ndalphardni_dxj_dTau2__constdelta_xi(rHEOS, i, j, xN_flag); @@ -1821,7 +1880,7 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") CAPTURE(analytic); CHECK(err < 1e-7); } - std::ostringstream ss3lll; ss3lll << "d3_ndalphardni_dxj_dDelta_dTau__constxi, i=" << i; + std::ostringstream ss3lll; ss3lll << "d3_ndalphardni_dxj_dDelta_dTau__constxi, i=" << i << ", j=" << j; SECTION(ss3lll.str(), "") { double analytic = MixtureDerivatives::d3_ndalphardni_dxj_dDelta_dTau__constxi(rHEOS, i, j, xN_flag); @@ -1833,7 +1892,7 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") CAPTURE(analytic); CHECK(err < 1e-7); } - std::ostringstream ss3m; ss3m << "d3alphar_dxi_dxj_dDelta, i=" << i; + std::ostringstream ss3m; ss3m << "d3alphar_dxi_dxj_dDelta, i=" << i << ", j=" << j; SECTION(ss3m.str(), "") { double analytic = MixtureDerivatives::d3alphar_dxi_dxj_dDelta(rHEOS, i, j, xN_flag); @@ -1845,7 +1904,7 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") CAPTURE(analytic); CHECK(err < 1e-7); } - std::ostringstream ss3mm; ss3mm << "d4alphar_dxi_dxj_dDelta2, i=" << i; + std::ostringstream ss3mm; ss3mm << "d4alphar_dxi_dxj_dDelta2, i=" << i << ", j=" << j; SECTION(ss3mm.str(), "") { double analytic = MixtureDerivatives::d4alphar_dxi_dxj_dDelta2(rHEOS, i, j, xN_flag); @@ -1857,7 +1916,7 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") CAPTURE(analytic); CHECK(err < 1e-7); } - std::ostringstream ss3mmm; ss3mmm << "d4alphar_dxi_dxj_dDelta_dTau, i=" << i; + std::ostringstream ss3mmm; ss3mmm << "d4alphar_dxi_dxj_dDelta_dTau, i=" << i << ", j=" << j; SECTION(ss3mmm.str(), "") { double analytic = MixtureDerivatives::d4alphar_dxi_dxj_dDelta_dTau(rHEOS, i, j, xN_flag); @@ -1869,7 +1928,7 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") CAPTURE(analytic); CHECK(err < 1e-7); } - std::ostringstream ss3n; ss3n << "d3alphar_dxi_dxj_dTau, i=" << i; + std::ostringstream ss3n; ss3n << "d3alphar_dxi_dxj_dTau, i=" << i << ", j=" << j; SECTION(ss3n.str(), "") { if (i == Ncomp-1){ break; } @@ -1882,7 +1941,7 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") CAPTURE(analytic); CHECK(err < 1e-7); } - std::ostringstream ss3nn; ss3nn << "d4alphar_dxi_dxj_dTau2, i=" << i; + std::ostringstream ss3nn; ss3nn << "d4alphar_dxi_dxj_dTau2, i=" << i << ", j=" << j; SECTION(ss3nn.str(), "") { if (i == Ncomp-1){ break; } @@ -1896,7 +1955,7 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") CHECK(err < 1e-7); } - std::ostringstream ss3o; ss3o << "d_nd_ndalphardni_dnj_dDelta__consttau_x, i=" << i; + std::ostringstream ss3o; ss3o << "d_nd_ndalphardni_dnj_dDelta__consttau_x, i=" << i << ", j=" << j; SECTION(ss3o.str(), "") { if (i == Ncomp-1){ break; } @@ -1909,7 +1968,7 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") CAPTURE(analytic); CHECK(err < 1e-7); } - std::ostringstream ss3o1; ss3o1 << "d_ndln_fugacity_i_dnj_ddelta__consttau_x, i=" << i; + std::ostringstream ss3o1; ss3o1 << "d_ndln_fugacity_i_dnj_ddelta__consttau_x, i=" << i << ", j=" << j; SECTION(ss3o1.str(), "") { if (i == Ncomp-1){ break; } @@ -1922,7 +1981,7 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") CAPTURE(analytic); CHECK(err < 1e-7); } - std::ostringstream ss3o2; ss3o2 << "d_ndln_fugacity_i_dnj_dtau__constdelta_x, i=" << i; + std::ostringstream ss3o2; ss3o2 << "d_ndln_fugacity_i_dnj_dtau__constdelta_x, i=" << i << ", j=" << j; SECTION(ss3o2.str(), "") { if (i == Ncomp-1){ break; } @@ -1936,7 +1995,7 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") CHECK(err < 1e-7); } - std::ostringstream ss3p; ss3p << "d_nd_ndalphardni_dnj_dTau__constdelta_x, i=" << i; + std::ostringstream ss3p; ss3p << "d_nd_ndalphardni_dnj_dTau__constdelta_x, i=" << i << ", j=" << j; SECTION(ss3p.str(), "") { if (i == Ncomp-1){ break; } @@ -1949,7 +2008,7 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") CAPTURE(analytic); CHECK(err < 1e-7); } - std::ostringstream ss3q; ss3q << "d_nddeltadni_dxj__constdelta_tau, i=" << i; + std::ostringstream ss3q; ss3q << "d_nddeltadni_dxj__constdelta_tau, i=" << i << ", j=" << j; SECTION(ss3q.str(), "") { if (i == Ncomp-1){ break; } @@ -1962,7 +2021,7 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") CAPTURE(analytic); CHECK(err < 1e-7); } - std::ostringstream ss3r; ss3r << "d_ndtaudni_dxj__constdelta_tau, i=" << i; + std::ostringstream ss3r; ss3r << "d_ndtaudni_dxj__constdelta_tau, i=" << i << ", j=" << j; SECTION(ss3r.str(), "") { if (i == Ncomp-1){ break; } @@ -2220,7 +2279,7 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") CAPTURE(analytic); CHECK(err < 1e-7); } - std::ostringstream ss3o5; ss3o5 << "d_n2Aijk_dTau, i=" << i << ", j=" << j; + std::ostringstream ss3o5; ss3o5 << "d_n2Aijk_dTau, i=" << i << ", j=" << j << ", k=" << k; SECTION(ss3o5.str(), "") { double analytic = MixtureDerivatives::d_n2Aijk_dX(rHEOS, i, j, k, xN_flag, CoolProp::iTau); @@ -2232,7 +2291,7 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") CAPTURE(analytic); CHECK(err < 1e-7); } - std::ostringstream ss3o6; ss3o6 << "d_n2Aijk_dDelta, i=" << i << ", j=" << j; + std::ostringstream ss3o6; ss3o6 << "d_n2Aijk_dDelta, i=" << i << ", j=" << j << ", k=" << k; SECTION(ss3o6.str(), "") { double analytic = MixtureDerivatives::d_n2Aijk_dX(rHEOS, i, j, k, xN_flag, CoolProp::iDelta); diff --git a/src/Backends/Helmholtz/MixtureDerivatives.h b/src/Backends/Helmholtz/MixtureDerivatives.h index 4d7721d0..e6b25868 100644 --- a/src/Backends/Helmholtz/MixtureDerivatives.h +++ b/src/Backends/Helmholtz/MixtureDerivatives.h @@ -284,7 +284,9 @@ class MixtureDerivatives{ static Eigen::MatrixXd Mstar(HelmholtzEOSMixtureBackend &HEOS, x_N_dependency_flag xN_flag){ std::size_t N = HEOS.mole_fractions.size(); - Eigen::MatrixXd L = Lstar(HEOS, xN_flag), M = L; + Eigen::MatrixXd L = Lstar(HEOS, xN_flag), + M = L, + adjL = adjugate(L); // Last row for (std::size_t i = 0; i < N; ++i){ @@ -292,15 +294,11 @@ class MixtureDerivatives{ for (std::size_t j = 0; j < N; ++j){ for (std::size_t k = j; k < N; ++k){ n2dLdni(j, k) = n2Aijk(HEOS, j, k, i, xN_flag); + // Fill in the symmetric elements + n2dLdni(k, j) = n2dLdni(j, k); } } - // Fill in the symmetric elements - for (std::size_t j = 0; j < N; ++j){ - for (std::size_t k = 0; k < j; ++k){ - n2dLdni(j, k) = n2dLdni(k, j); - } - } - M(N-1, i) = (adjugate(L)*n2dLdni).trace(); + M(N-1, i) = (adjL*n2dLdni).trace(); } return M; } @@ -313,20 +311,16 @@ class MixtureDerivatives{ adjL = adjugate(L), d_adjL_dX = adjugate_derivative(L, dL_dX); - // Last row + // Last row in the d(Mstar)/d(X) requires derivatives of for (std::size_t i = 0; i < N; ++i){ Eigen::MatrixXd n2dLdni(N, N), d_n2dLdni_dX(N, N); for (std::size_t j = 0; j < N; ++j){ for (std::size_t k = j; k < N; ++k){ n2dLdni(j, k) = n2Aijk(HEOS, j, k, i, xN_flag); d_n2dLdni_dX(j, k) = d_n2Aijk_dX(HEOS, j, k, i, xN_flag, WRT); - } - } - // Fill in the symmetric elements - for (std::size_t j = 0; j < N; ++j){ - for (std::size_t k = 0; k < j; ++k){ - n2dLdni(j, k) = n2dLdni(k, j); - d_n2dLdni_dX(j, k) = d_n2dLdni_dX(k, j); + // Fill in the symmetric elements + n2dLdni(k, j) = n2dLdni(j, k); + d_n2dLdni_dX(k, j) = d_n2dLdni_dX(j, k); } } dMstar(N-1, i) = (n2dLdni*d_adjL_dX + adjL*d_n2dLdni_dX).trace(); From 1283143e649b240b35fcc545c98592f803e5f034 Mon Sep 17 00:00:00 2001 From: Ian Bell Date: Sun, 16 Aug 2015 16:56:27 -0600 Subject: [PATCH 7/8] Fix bug in adjugate function and fix tests thereof --- include/MatrixMath.h | 2 +- src/Backends/Helmholtz/MixtureDerivatives.cpp | 21 +++++++++++-------- 2 files changed, 13 insertions(+), 10 deletions(-) diff --git a/include/MatrixMath.h b/include/MatrixMath.h index eed5354d..30087dfa 100644 --- a/include/MatrixMath.h +++ b/include/MatrixMath.h @@ -908,7 +908,7 @@ static Eigen::MatrixXd adjugate(const Eigen::MatrixBase& A) for (std::size_t i = 0; i < N; ++i){ for (std::size_t j = 0; j < N; ++j){ int negative_1_to_the_i_plus_j = ((i+j)%2==0) ? 1 : -1; - Aadj(i, j) = negative_1_to_the_i_plus_j*minor_matrix(A, i, j).determinant(); + Aadj(i, j) = negative_1_to_the_i_plus_j*minor_matrix(A, j, i).determinant(); } } return Aadj; diff --git a/src/Backends/Helmholtz/MixtureDerivatives.cpp b/src/Backends/Helmholtz/MixtureDerivatives.cpp index bde093ff..0e53cd36 100644 --- a/src/Backends/Helmholtz/MixtureDerivatives.cpp +++ b/src/Backends/Helmholtz/MixtureDerivatives.cpp @@ -1033,15 +1033,15 @@ static std::vector > > HEOS, 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; +static const double T1 = 319.325, rho1 = 13246.6, dT = 1e-3, drho = 1e-3, dz = 1e-6; void setup_state(std::vector > & HEOS, std::size_t Ncomp, double increment, x_N_dependency_flag xN_flag = XN_INDEPENDENT) { std::vector names(Ncomp); std::vector z(Ncomp); if (Ncomp == 2){ - names[0] = "Ethane"; names[1] = "Propane"; - z[0] = 0.3; z[1] = 0.7; + names[0] = "Methane"; names[1] = "H2S"; + z[0] = 0.4; z[1] = 0.6; } else if (Ncomp == 3){ names[0] = "Ethane"; names[1] = "Propane"; names[2] = "Methane"; @@ -1180,7 +1180,7 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") Eigen::MatrixXd Mstar = MixtureDerivatives::Mstar(rHEOS, xN_flag); Eigen::MatrixXd analytic = adjugate(Mstar); Eigen::MatrixXd numeric(2,2); - numeric << Mstar(1,1), -Mstar(1,0), -Mstar(0,1), Mstar(0,0); + numeric << Mstar(1,1), -Mstar(0,1), -Mstar(1,0), Mstar(0,0); double err = ((analytic-numeric).array()/analytic.array()).cwiseAbs().sum()/Ncomp/Ncomp; CAPTURE(numeric); CAPTURE(analytic); @@ -1207,15 +1207,18 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") { Eigen::MatrixXd Mstar = MixtureDerivatives::Mstar(rHEOS, xN_flag); Eigen::MatrixXd dMstar_dTau = MixtureDerivatives::dMstar_dX(rHEOS, xN_flag, CoolProp::iTau); - double analytic = (adjugate(Mstar)*dMstar_dTau).trace(); + Eigen::MatrixXd adjM = adjugate(Mstar); + double analytic = (adjM*dMstar_dTau).trace(); double detMstar_plus = MixtureDerivatives::Mstar(rHEOS_plusT_constrho, xN_flag).determinant(); double tau1 = rHEOS_plusT_constrho.tau(); double detMstar_minus = MixtureDerivatives::Mstar(rHEOS_minusT_constrho, xN_flag).determinant(); double tau2 = rHEOS_minusT_constrho.tau(); - + double numeric = (detMstar_plus - detMstar_minus)/(tau1-tau2); double err = mix_deriv_err_func(numeric, analytic); CAPTURE(numeric); CAPTURE(analytic); + CAPTURE(dMstar_dTau); + CAPTURE(adjM); CHECK(err < 1e-8); } SECTION("d(M1)/dDelta", "") @@ -1231,7 +1234,7 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") double err = mix_deriv_err_func(numeric, analytic); CAPTURE(numeric); CAPTURE(analytic); - CHECK(err < -1e-8); + CHECK(err < 1e-8); } SECTION("d(L1)/dDelta", "") { @@ -1246,7 +1249,7 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") double err = mix_deriv_err_func(numeric, analytic); CAPTURE(numeric); CAPTURE(analytic); - CHECK(err < -1e-8); + CHECK(err < 1e-8); } SECTION("adj(Lstar)", "") @@ -1255,7 +1258,7 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") Eigen::MatrixXd Lstar = MixtureDerivatives::Lstar(rHEOS, xN_flag); Eigen::MatrixXd analytic = adjugate(Lstar); Eigen::MatrixXd numeric(2,2); - numeric << Lstar(1,1), -Lstar(1,0), -Lstar(0,1), Lstar(0,0); + numeric << Lstar(1,1), -Lstar(0,1), -Lstar(1,0), Lstar(0,0); double err = ((analytic-numeric).array()/analytic.array()).cwiseAbs().sum()/Ncomp/Ncomp; CAPTURE(numeric); CAPTURE(analytic); From 30c5abc35cac75dc1d729adc37ab5c3dc9025b6c Mon Sep 17 00:00:00 2001 From: Ian Bell Date: Sun, 16 Aug 2015 17:04:01 -0600 Subject: [PATCH 8/8] Fully switch over to analytic derivatives; add an IPython notebook to check derivative of determinant --- ...k critical point derivative matrices.ipynb | 162 ++++++++++++++++++ .../Helmholtz/HelmholtzEOSMixtureBackend.cpp | 41 +++-- 2 files changed, 182 insertions(+), 21 deletions(-) create mode 100644 doc/notebooks/Check critical point derivative matrices.ipynb diff --git a/doc/notebooks/Check critical point derivative matrices.ipynb b/doc/notebooks/Check critical point derivative matrices.ipynb new file mode 100644 index 00000000..1779cfa8 --- /dev/null +++ b/doc/notebooks/Check critical point derivative matrices.ipynb @@ -0,0 +1,162 @@ +{ + "metadata": { + "name": "", + "signature": "sha256:a09a62c5ca301d3f63e20ae1ddea45013e338af6b73280fca44334e814ae310d" + }, + "nbformat": 3, + "nbformat_minor": 0, + "worksheets": [ + { + "cells": [ + { + "cell_type": "code", + "collapsed": false, + "input": [ + "from __future__ import division\n", + "from sympy import *\n", + "from IPython.display import display, Math, Latex\n", + "from IPython.core.display import display_html\n", + "init_session(quiet=True, use_latex='mathjax')\n", + "init_printing()" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "output_type": "stream", + "stream": "stdout", + "text": [ + "IPython console for SymPy 0.7.6 (Python 2.7.6-32-bit) (ground types: python)\n" + ] + } + ], + "prompt_number": 1 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "x1,x2,tau,delta = symbols('x1,x2,tau,delta')\n", + "L11 = symbols('L11', cls=Function)(x1,x2,tau,delta)\n", + "L12 = symbols('L12', cls=Function)(x1,x2,tau,delta)\n", + "L21 = symbols('L21', cls=Function)(x1,x2,tau,delta)\n", + "L22 = symbols('L22', cls=Function)(x1,x2,tau,delta)" + ], + "language": "python", + "metadata": {}, + "outputs": [], + "prompt_number": 2 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "Lstar = Matrix([[L11,L12],[L21,L22]])\n", + "L1star = Lstar.det()\n", + "deriv1 = L1star.diff(t)\n", + "deriv2 = (Lstar.adjugate()*Lstar.diff(t)).trace()\n", + "\n", + "Mstar = Matrix([[L11,L12],[Lstar.det().diff(x1),Lstar.det().diff(x2)]])\n", + "\n", + "deriv1 = Mstar.det().diff(tau)\n", + "deriv2 = (Mstar.adjugate()*Mstar.diff(tau)).trace()\n", + "simplify(deriv1-deriv2)" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "latex": [ + "$$0$$" + ], + "metadata": {}, + "output_type": "pyout", + "png": "iVBORw0KGgoAAAANSUhEUgAAAAoAAAAOBAMAAADkjZCYAAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAEJmJdjLNVN0iZu+7\nq0QgoRR7AAAAVklEQVQIHWNgEDJRZWBgSGeQmMDAtYGBOYGB5wID+0cG/gsMfN8Z5BUY+L4wzDdg\nYP0MJeUNQCL8Cgzs3xk4DjBwfWRg2cDAlMDA0M4gHcDAIOxylQEA9FISlFfRJtkAAAAASUVORK5C\nYII=\n", + "prompt_number": 3, + "text": [ + "0" + ] + } + ], + "prompt_number": 3 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "Mstar = Matrix([[0.00112865, 8.76232e-006],[9.57021e-007, 7.42578e-009]])\n", + "dMstar_dTau = Matrix([[-0.000245724,-0.00118232],[3.20921e-006, -7.171e-008]])\n", + "adjM = Matrix([[7.42578e-009,-9.57021e-007],[-8.76232e-006, 0.00112865]])\n", + "(adjM*dMstar_dTau).trace()\n", + "print Mstar.adjugate()\n", + "print adjM" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "output_type": "stream", + "stream": "stdout", + "text": [ + "Matrix([[7.42578000000000e-9, -8.76232000000000e-6], [-9.57021000000000e-7, 0.00112865000000000]])\n", + "Matrix([[7.42578000000000e-9, -9.57021000000000e-7], [-8.76232000000000e-6, 0.00112865000000000]])\n" + ] + } + ], + "prompt_number": 7 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "a,b,c,d = symbols('a,b,c,d')\n", + "M = Matrix([[a,b],[c,d]])\n", + "display(M)\n", + "M.adjugate()" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "latex": [ + "$$\\left[\\begin{matrix}a & b\\\\c & d\\end{matrix}\\right]$$" + ], + "metadata": {}, + "output_type": "display_data", + "png": "iVBORw0KGgoAAAANSUhEUgAAADcAAAAyBAMAAAAKOF7GAAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMA74lUMhAiq3bdRLtm\nmc0lg57xAAABSUlEQVQ4Ee3UPUsDMRwG8OdSr9DT1oIuouANDgoO9w0qWPdzFrlu7rp0cOhQEARB\nEHHVT2D9AIqjo4Pgpi7iVBGlDlI4E/LSZ7hmdjDLJfnxz11yD8Fi/omiJvK8jpnmRpEhbK7XMVtI\nanKiGLfmPIg3H355MBp4UBwTlne7p2poW619lLivvWhNflhQz+wV3xbDBVR7jI0U8tj0PqsDVDqM\nTwgdBj1ka4x9RHIzujKL0UBKOoRczGKCm1KL8AebtxZrnakzQYZ9XMuhXjZsL+0cMm4fyEKDPE99\nXUkT3P1HeRpjDyG6uh+PmE48eJl6cMXzTvRHKJa7/LNkyF3AUD7BO6MMuQsYghh7hCrkL27ZZ84l\nVC4rscN5KpPdQIb83KLKPjcTcnN8Q6BEmumQG3xE+EBoQm5QrN6RwYTcIMuo/xfRe6X6LuNfjWlM\nFpMM9N8AAAAASUVORK5CYII=\n", + "text": [ + "\u23a1a b\u23a4\n", + "\u23a2 \u23a5\n", + "\u23a3c d\u23a6" + ] + }, + { + "latex": [ + "$$\\left[\\begin{matrix}d & - b\\\\- c & a\\end{matrix}\\right]$$" + ], + "metadata": {}, + "output_type": "pyout", + "png": "iVBORw0KGgoAAAANSUhEUgAAAFMAAAAyBAMAAADSNPrMAAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMA74lUMhBEu5ndzSKr\ndmb+gm2XAAABbElEQVRIDWOQ//+JgTBg+v9fgEHYxZWwSgZWF2cBBhEiFIKUsCArTZfArou1cAK6\nUobP2JUycB1AV8r2FYdSfgN0pUwge7CB/ACgKIpbuRqwqQOKaYLEkZTu0fUH2YMNrLnbjKyUTZqh\nHmQPFsD6K4B/A5KpgQ8YurAoAwkxf2XgV0BS2h/AMAlNKevKmUAwy4FjAUM8slJxBtZvaEphXKB3\n8w0QprJ+YmD+wAqTRKX5HzAA7YSHAOtnBvYF1qhKYDz+BFY55BBYy/DygAJMEpXmecDegKw05pL3\n2QRUJTAea+85IBPuAJgwbnpU6WgIDLoQULrrgJZi9+h6QooHNLfuNOBagKoUVJJARNCUCjHwoBVG\nwJKkHVkpvGzALDWBuXouslKYnYwFMBacBpYkX7ApBRY3aACpJEF1K7AYYUB1K1JJgqqU4wBDhAOq\nuYiSBFUpg9K7C6gqGRAlCZpSNHUo3OGrlIRKnvimAwCJ/VUvfMvpJAAAAABJRU5ErkJggg==\n", + "prompt_number": 11, + "text": [ + "\u23a1d -b\u23a4\n", + "\u23a2 \u23a5\n", + "\u23a3-c a \u23a6" + ] + } + ], + "prompt_number": 11 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [], + "language": "python", + "metadata": {}, + "outputs": [] + } + ], + "metadata": {} + } + ] +} \ No newline at end of file diff --git a/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp b/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp index 837fc22a..16d1cb9b 100644 --- a/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp +++ b/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp @@ -2915,34 +2915,34 @@ void HelmholtzEOSMixtureBackend::calc_critical_point(double rho0, double T0) }; std::vector > Jacobian(const std::vector &x) { - + std::size_t N = x.size(); + std::vector > J(N, std::vector(N, 0)); + Eigen::MatrixXd adjL = adjugate(MixtureDerivatives::Lstar(HEOS, XN_INDEPENDENT)), + adjM = adjugate(MixtureDerivatives::Mstar(HEOS, XN_INDEPENDENT)), + dLdTau = MixtureDerivatives::dLstar_dX(HEOS, XN_INDEPENDENT, iTau), + dLdDelta = MixtureDerivatives::dLstar_dX(HEOS, XN_INDEPENDENT, iDelta), + dMdTau = MixtureDerivatives::dMstar_dX(HEOS, XN_INDEPENDENT, iTau), + dMdDelta = MixtureDerivatives::dMstar_dX(HEOS, XN_INDEPENDENT, iDelta); + + J[0][0] = (adjL*dLdTau).trace(); + J[0][1] = (adjL*dLdDelta).trace(); + J[1][0] = (adjM*dMdTau).trace(); + J[1][1] = (adjM*dMdDelta).trace(); + return J; + } + /// Not used, for testing purposes + std::vector > numerical_Jacobian(const std::vector &x) + { std::size_t N = x.size(); std::vector rplus, rminus, xp; - std::vector > J(N, std::vector(N, 0)), JJ(N, std::vector(N, 0)); - Eigen::MatrixXd J0(N, N), adjL = adjugate(MixtureDerivatives::Lstar(HEOS, XN_INDEPENDENT)), - adjM = adjugate(MixtureDerivatives::Mstar(HEOS, XN_INDEPENDENT)), - dLdTau = MixtureDerivatives::dLstar_dX(HEOS, XN_INDEPENDENT, iTau), - dLdDelta = MixtureDerivatives::dLstar_dX(HEOS, XN_INDEPENDENT, iDelta), - dMdTau = MixtureDerivatives::dMstar_dX(HEOS, XN_INDEPENDENT, iTau), - dMdDelta = MixtureDerivatives::dMstar_dX(HEOS, XN_INDEPENDENT, iDelta); - - JJ[0][0] = (adjL*dLdTau).trace(); - JJ[0][1] = (adjL*dLdDelta).trace(); - JJ[1][0] = (adjM*dMdTau).trace(); - JJ[1][1] = (adjM*dMdDelta).trace(); - -// return JJ; + std::vector > J(N, std::vector(N, 0)); - //std::cout << J0 << std::endl; - std::cout << JJ[0][0] << " " << JJ[0][1] << std::endl; - std::cout << JJ[1][0] << " " << JJ[1][1] << std::endl; - double epsilon; + double epsilon = 0.0001; // Build the Jacobian by column for (std::size_t i = 0; i < N; ++i) { xp = x; - epsilon = 0.0001; xp[i] += epsilon; rplus = call(xp); xp[i] -= 2*epsilon; @@ -2956,7 +2956,6 @@ void HelmholtzEOSMixtureBackend::calc_critical_point(double rho0, double T0) std::cout << J[0][0] << " " << J[0][1] << std::endl; std::cout << J[1][0] << " " << J[1][1] << std::endl; return J; - }; }; Resid resid(*this);