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/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/HelmholtzEOSMixtureBackend.cpp b/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp index 8b886d11..2533c0c0 100644 --- a/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp +++ b/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp @@ -2932,34 +2932,46 @@ 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 > J(N, std::vector(N, 0)); - Eigen::MatrixXd J0(N, N), adjL = adjugate(MixtureDerivatives::Lstar(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; + 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); - J0(0,0) = (adjL*dLdTau).trace(); - J0(0,1) = (adjL*dLdDelta).trace(); - //std::cout << J0 << std::endl; - std::vector r0 = call(x); + 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)); + + double epsilon = 0.0001; + // Build the Jacobian by column for (std::size_t i = 0; i < N; ++i) { xp = x; - epsilon = 0.00001; 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; }; }; @@ -2987,13 +2999,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), @@ -3050,7 +3060,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 0fe4fb3a..0e53cd36 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_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--; } + 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) - 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--; } + 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) - 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++) + { + 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) { @@ -837,16 +956,51 @@ 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; } +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 of dalphar/(dxi,dxj,dxk) are zero + double term3 = -2*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 of dalphar)/(dxi,dxj,dxk) are zero + double term3 = -2*d3alphar_dxi_dxj_dDelta(HEOS, j, k, xN_flag); + return term1 + term2 + term3; +} + + } /* namespace CoolProp */ #ifdef ENABLE_CATCH @@ -854,6 +1008,16 @@ CoolPropDbl MixtureDerivatives::d2_ndalphardni_dxj_dxk__constdelta_tau_xi(Helmho 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 std::abs(numeric/analytic-1); + } +} + static bool fluids_set = false; static const std::size_t Ncomp_max = 6; @@ -869,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"; @@ -1009,6 +1173,144 @@ 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(); + + 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(0,1), -Mstar(1,0), 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); + 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", "") + { + 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(0,1), -Lstar(1,0), 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(); + 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-8); + } + 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(); + 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-8); + } + 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(); + 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-8); + } + 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-8); + } // These ones only require the i index for (std::size_t i = 0; i< Ncomp; ++i) @@ -1027,7 +1329,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); @@ -1040,7 +1342,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); @@ -1053,7 +1355,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); @@ -1066,7 +1368,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; @@ -1116,31 +1418,29 @@ 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; 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(); 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; 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(); 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; @@ -1151,19 +1451,18 @@ 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; 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(); 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); @@ -1172,12 +1471,11 @@ 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(); 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); @@ -1186,12 +1484,11 @@ 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(); 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); @@ -1200,12 +1497,11 @@ 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(); 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); @@ -1214,12 +1510,11 @@ 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()); 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); @@ -1228,12 +1523,11 @@ 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()); 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); @@ -1247,7 +1541,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); @@ -1256,12 +1550,11 @@ 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(); 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); @@ -1269,12 +1562,11 @@ 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(); 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); @@ -1282,38 +1574,85 @@ 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(); 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); } + std::ostringstream ss3ff; ss3ff << "d4alphar_dxi_dTau3, i=" << i; + SECTION(ss3ff.str(), "") + { + 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 = mix_deriv_err_func(numeric, analytic); + CAPTURE(numeric); + CAPTURE(analytic); + CHECK(err < 1e-7); + } + 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(); 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); } + std::ostringstream ss3gg; ss3gg << "d4alphar_dxi_dDelta_Tau2, i=" << i; + SECTION(ss3gg.str(), "") + { + 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 = mix_deriv_err_func(numeric, analytic); + CAPTURE(numeric); + CAPTURE(analytic); + CHECK(err < 1e-7); + } + std::ostringstream ss3ggg; ss3ggg << "d4alphar_dxi_dDelta2_Tau, i=" << i; + SECTION(ss3ggg.str(), "") + { + 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 = mix_deriv_err_func(numeric, analytic); + CAPTURE(numeric); + CAPTURE(analytic); + CHECK(err < 1e-7); + } + 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(); 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); + } + std::ostringstream ss3hh; ss3hh << "d3_ndalphardni_dTau3, i=" << i; + SECTION(ss3hh.str(), "") + { + 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 = mix_deriv_err_func(numeric, analytic); CAPTURE(numeric); CAPTURE(analytic); CHECK(err < 1e-7); @@ -1321,12 +1660,35 @@ 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(); 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); + } + std::ostringstream ss3ii; ss3ii << "d3_ndalphardni_dDelta_dTau2, i=" << i; + SECTION(ss3ii.str(), "") + { + 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 = mix_deriv_err_func(numeric, analytic); + CAPTURE(numeric); + CAPTURE(analytic); + CHECK(err < 1e-7); + } + std::ostringstream ss3iii; ss3iii << "d3_ndalphardni_dDelta2_dTau, i=" << i; + SECTION(ss3iii.str(), "") + { + 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 = mix_deriv_err_func(numeric, analytic); CAPTURE(numeric); CAPTURE(analytic); CHECK(err < 1e-7); @@ -1334,12 +1696,23 @@ 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(); 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); + } + std::ostringstream ss3k; ss3k << "d3_ndalphardni_dDelta3, i=" << i; + SECTION(ss3k.str(), "") + { + 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 = mix_deriv_err_func(numeric, analytic); CAPTURE(numeric); CAPTURE(analytic); CHECK(err < 1e-7); @@ -1359,12 +1732,11 @@ 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)); 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); @@ -1373,12 +1745,11 @@ 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); 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,12 +1758,11 @@ 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); 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,12 +1771,11 @@ 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); 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); @@ -1415,20 +1784,12 @@ 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); 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); @@ -1437,13 +1798,11 @@ 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); 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); @@ -1456,8 +1815,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); @@ -1469,13 +1827,12 @@ 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); } - 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; } @@ -1483,12 +1840,25 @@ 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); } - std::ostringstream ss3l; ss3l << "d2_ndalphardni_dxj_dTau, i=" << i; + std::ostringstream ss3kk; ss3kk << "d3_ndalphardni_dxj_dDelta2__consttau_xi, i=" << i << ", j=" << j; + 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 = mix_deriv_err_func(numeric, analytic); + CAPTURE(numeric); + CAPTURE(analytic); + CHECK(err < 1e-7); + } + std::ostringstream ss3l; ss3l << "d2_ndalphardni_dxj_dTau, i=" << i << ", j=" << j; SECTION(ss3l.str(), "") { if (i == Ncomp-1){ break; } @@ -1496,25 +1866,72 @@ 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); } - std::ostringstream ss3m; ss3m << "d3alphar_dxi_dxj_dDelta, 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); + 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 = mix_deriv_err_func(numeric, analytic); + CAPTURE(numeric); + CAPTURE(analytic); + CHECK(err < 1e-7); + } + 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); + 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 = mix_deriv_err_func(numeric, analytic); + CAPTURE(numeric); + CAPTURE(analytic); + CHECK(err < 1e-7); + } + std::ostringstream ss3m; ss3m << "d3alphar_dxi_dxj_dDelta, i=" << i << ", j=" << j; 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); 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 ss3n; ss3n << "d3alphar_dxi_dxj_dTau, 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); + 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 = mix_deriv_err_func(numeric, analytic); + CAPTURE(numeric); + CAPTURE(analytic); + CHECK(err < 1e-7); + } + 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); + 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 = mix_deriv_err_func(numeric, analytic); + CAPTURE(numeric); + CAPTURE(analytic); + CHECK(err < 1e-7); + } + std::ostringstream ss3n; ss3n << "d3alphar_dxi_dxj_dTau, i=" << i << ", j=" << j; SECTION(ss3n.str(), "") { if (i == Ncomp-1){ break; } @@ -1522,12 +1939,26 @@ 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); } - std::ostringstream ss3o; ss3o << "d_nd_ndalphardni_dnj_dDelta__consttau_x, i=" << i; + std::ostringstream ss3nn; ss3nn << "d4alphar_dxi_dxj_dTau2, i=" << i << ", j=" << j; + 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 = mix_deriv_err_func(numeric, analytic); + CAPTURE(numeric); + CAPTURE(analytic); + CHECK(err < 1e-7); + } + + std::ostringstream ss3o; ss3o << "d_nd_ndalphardni_dnj_dDelta__consttau_x, i=" << i << ", j=" << j; SECTION(ss3o.str(), "") { if (i == Ncomp-1){ break; } @@ -1535,12 +1966,12 @@ 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); } - 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; } @@ -1548,12 +1979,12 @@ 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); } - 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; } @@ -1561,13 +1992,13 @@ 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); } - 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; } @@ -1575,12 +2006,12 @@ 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); } - 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; } @@ -1588,12 +2019,12 @@ 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); } - 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; } @@ -1601,7 +2032,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); @@ -1623,8 +2104,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); @@ -1637,7 +2117,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); @@ -1650,7 +2130,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); @@ -1663,7 +2143,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); @@ -1676,7 +2156,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); @@ -1689,7 +2169,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); @@ -1702,7 +2182,33 @@ 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); + } + 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 = mix_deriv_err_func(numeric, 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 = mix_deriv_err_func(numeric, analytic); CAPTURE(numeric); CAPTURE(analytic); CHECK(err < 1e-7); @@ -1716,14 +2222,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); @@ -1738,14 +2237,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); @@ -1759,18 +2251,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 << ", k=" << k; + 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 << ", k=" << k; + 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); + } } } } diff --git a/src/Backends/Helmholtz/MixtureDerivatives.h b/src/Backends/Helmholtz/MixtureDerivatives.h index 07b48adc..e6b25868 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,51 @@ 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(); + 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); + 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); + } + 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){ + 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){ std::size_t N = HEOS.mole_fractions.size(); Eigen::MatrixXd L; @@ -198,16 +247,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 @@ -244,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){ @@ -252,18 +294,39 @@ 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; } + 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 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 + 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(); + } + return dMstar; + } /** \brief Table B4, Kunz, JCED, 2012 for the original term and the subsequent substitutions * @@ -497,6 +560,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 +605,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 +617,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 +637,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 +655,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 +679,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*/