From d0ee80af5fc63dce2fd6347a831f3fffe65fa073 Mon Sep 17 00:00:00 2001 From: Ian Bell Date: Wed, 24 Aug 2016 22:53:45 -0600 Subject: [PATCH] Update MixtureDerivatives.cpp for VTPR; closes #1201 --- src/Backends/Helmholtz/MixtureDerivatives.cpp | 33 +++++++++++++++++-- src/Backends/Helmholtz/MixtureDerivatives.h | 10 ++++++ 2 files changed, 40 insertions(+), 3 deletions(-) diff --git a/src/Backends/Helmholtz/MixtureDerivatives.cpp b/src/Backends/Helmholtz/MixtureDerivatives.cpp index 8fc7532a..5e23a197 100644 --- a/src/Backends/Helmholtz/MixtureDerivatives.cpp +++ b/src/Backends/Helmholtz/MixtureDerivatives.cpp @@ -765,6 +765,7 @@ CoolPropDbl MixtureDerivatives::dalpha0_dxi(HelmholtzEOSMixtureBackend &HEOS, st #include "catch.hpp" #include "Backends/Cubics/CubicBackend.h" +#include "Backends/Cubics/VTPRBackend.h" using namespace CoolProp; @@ -800,7 +801,7 @@ public: double dtau, ddelta, dz, dn, tol, dT, drho, dp; DerivativeFixture() : xN(XN_INDEPENDENT) { dtau = 1e-6; ddelta = 1e-6; dz = 1e-6; dn = 1e-6; dT = 1e-3; drho = 1e-3; dp = 1; tol = 5e-6; - std::vector names; names.push_back("Methane"); names.push_back("Ethane"); names.push_back("Propane"); names.push_back("n-Butane"); + std::vector names; names.push_back("n-Pentane"); names.push_back("Ethane"); names.push_back("n-Propane"); names.push_back("n-Butane"); std::vector mole_fractions; mole_fractions.push_back(0.1); mole_fractions.push_back(0.2); mole_fractions.push_back(0.3); mole_fractions.push_back(0.4); HEOS.reset(new backend(names)); HEOS->set_mole_fractions(mole_fractions); @@ -873,6 +874,23 @@ public: other->set_mole_fractions(HEOS->get_mole_fractions()); other->specify_phase(CoolProp::iphase_gas); // Something homogeneous } + void zero(const std::string &name, zero_mole_fraction_pointer f, zero_mole_fraction_pointer g = NULL, derivative wrt = NO_DERIV){ + double analytic = f(*HEOS, xN); + double numeric = 0; + if (wrt == TAU){ + numeric = (g(*HEOS_plus_tau, xN) - g(*HEOS_minus_tau, xN))/(2*dtau); + } + else if (wrt == DELTA){ + numeric = (g(*HEOS_plus_delta, xN) - g(*HEOS_minus_delta, xN))/(2*ddelta); + } + CAPTURE(name); + CAPTURE(analytic) + CAPTURE(numeric); + CAPTURE(xN); + double error = mix_deriv_err_func(numeric, analytic); + CAPTURE(error); + CHECK(error < tol); + } void one(const std::string &name, one_mole_fraction_pointer f, one_mole_fraction_pointer g = NULL, derivative wrt = NO_DERIV){ for (int i = 0; i < 4; ++i){ double analytic = f(*HEOS, i, xN); @@ -1076,6 +1094,11 @@ public: one("d2_ndalphardni_dDelta_dTau", MD::d2_ndalphardni_dDelta_dTau, MD::d_ndalphardni_dDelta, TAU); one("d3_ndalphardni_dDelta2_dTau", MD::d3_ndalphardni_dDelta2_dTau, MD::d2_ndalphardni_dDelta2, TAU); one("d3_ndalphardni_dDelta_dTau2", MD::d3_ndalphardni_dDelta_dTau2, MD::d2_ndalphardni_dDelta_dTau, TAU); + + zero("dalphar_dDelta", MD::dalphar_dDelta, MD::alphar, DELTA); + zero("d2alphar_dDelta2", MD::d2alphar_dDelta2, MD::dalphar_dDelta, DELTA); + zero("dalphar_dTau", MD::dalphar_dTau, MD::alphar, TAU); + zero("d2alphar_dTau2", MD::d2alphar_dTau2, MD::dalphar_dTau, TAU); //two_comp("d_ndalphardni_dxj__constT_V_xi", MD::d_ndalphardni_dxj__constT_V_xi, MD::ndalphar_dni__constT_V_nj); @@ -1174,17 +1197,21 @@ TEST_CASE_METHOD(DerivativeFixture, "Check derivativ { run_checks(); }; - TEST_CASE_METHOD(DerivativeFixture, "Check derivatives for Peng-Robinson", "[mixture_derivs2]") { tol = 1e-4; // Relax the tolerance a bit run_checks(); }; - TEST_CASE_METHOD(DerivativeFixture, "Check derivatives for SRK", "[mixture_derivs2]") { tol = 1e-4; // Relax the tolerance a bit run_checks(); }; +// Make sure you set the VTPR UNIFAQ path with something like set_config_string(VTPR_UNIFAQ_PATH, "/Users/ian/Code/CUBAC/dev/unifaq/"); +//TEST_CASE_METHOD(DerivativeFixture, "Check derivatives for VTPR", "[mixture_derivs2]") +//{ +// tol = 1e-4; // Relax the tolerance a bit +// run_checks(); +//}; #endif diff --git a/src/Backends/Helmholtz/MixtureDerivatives.h b/src/Backends/Helmholtz/MixtureDerivatives.h index 075eab25..548ff575 100644 --- a/src/Backends/Helmholtz/MixtureDerivatives.h +++ b/src/Backends/Helmholtz/MixtureDerivatives.h @@ -816,6 +816,16 @@ class MixtureDerivatives{ static CoolPropDbl dalphar_dDelta(CoolProp::HelmholtzEOSMixtureBackend &HEOS, CoolProp::x_N_dependency_flag xN_flag){ return HEOS.dalphar_dDelta(); } + static CoolPropDbl d2alphar_dDelta2(CoolProp::HelmholtzEOSMixtureBackend &HEOS, CoolProp::x_N_dependency_flag xN_flag){ + return HEOS.d2alphar_dDelta2(); + } + static CoolPropDbl dalphar_dTau(CoolProp::HelmholtzEOSMixtureBackend &HEOS, CoolProp::x_N_dependency_flag xN_flag){ + return HEOS.dalphar_dTau(); + } + static CoolPropDbl d2alphar_dTau2(CoolProp::HelmholtzEOSMixtureBackend &HEOS, CoolProp::x_N_dependency_flag xN_flag){ + return HEOS.d2alphar_dTau2(); + } + static CoolPropDbl ln_fugacity(CoolProp::HelmholtzEOSMixtureBackend &HEOS, std::size_t i, CoolProp::x_N_dependency_flag xN_flag){ return log(MixtureDerivatives::fugacity_i(HEOS, i, xN_flag)); }