diff --git a/src/coreComponents/constitutive/unitTests/TestFluidUtilities.hpp b/src/coreComponents/constitutive/unitTests/TestFluidUtilities.hpp new file mode 100644 index 00000000000..ec2cea463ca --- /dev/null +++ b/src/coreComponents/constitutive/unitTests/TestFluidUtilities.hpp @@ -0,0 +1,68 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2018-2020 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2020 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2018-2020 TotalEnergies + * Copyright (c) 2019- GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +/** + * @file TestFluidUtilities.hpp + */ + +#ifndef GEOS_CONSTITUTIVE_UNITTESTS_TESTFLUIDUTILITIES_HPP_ +#define GEOS_CONSTITUTIVE_UNITTESTS_TESTFLUIDUTILITIES_HPP_ + +#include "codingUtilities/UnitTestUtilities.hpp" + +namespace geos +{ + +namespace testing +{ + +namespace internal +{ + +static constexpr real64 relTol = 1.0e-4; +static constexpr real64 absTol = 1.0e-7; + +/** + * @brief Tests a function against a derivative + * @details Will calculate the left-sided and the right-sided numerical derivatives of a function + * and compare this against a analytically calculated value provided. + * @param x The value at which the function should be evaluated + * @param dx The value to use to perturb @c x in the calculation of the numerical derivatives + * @param derivative The value of the analytically calculated derivative to use for comparison + * @param function The function which is being tested. This should be a function that returns + a single real value. + * @param absTolerance The absolute tolerance to use for the comparison + * @param relTolerance The relative tolerance to use for the comparison + */ +template< typename FUNCTION > +void testNumericalDerivative( real64 const x, real64 const dx, real64 const derivative, + FUNCTION && function, + real64 const absTolerance = absTol, real64 const relTolerance = relTol ) +{ + real64 const leftValue = function( x-dx ); + real64 const centreValue = function( x ); + real64 const rightValue = function( x+dx ); + real64 const leftDerivative = (centreValue - leftValue) / dx; + real64 const rightDerivative = (rightValue - centreValue) / dx; + checkRelativeError( derivative, leftDerivative, relTolerance, absTolerance, "Left derivative" ); + checkRelativeError( derivative, rightDerivative, relTolerance, absTolerance, "Right derivative" ); +} + +}// namespace internal + +}// namespace testing + +}// namespace geos + +#endif //GEOS_CONSTITUTIVE_UNITTESTS_TESTFLUIDUTILITIES_HPP_ diff --git a/src/coreComponents/constitutive/unitTests/testCompositionalProperties.cpp b/src/coreComponents/constitutive/unitTests/testCompositionalProperties.cpp index 9866f027326..3680bfe236c 100644 --- a/src/coreComponents/constitutive/unitTests/testCompositionalProperties.cpp +++ b/src/coreComponents/constitutive/unitTests/testCompositionalProperties.cpp @@ -13,19 +13,16 @@ */ // Source includes -#include "codingUtilities/UnitTestUtilities.hpp" #include "constitutive/fluid/multifluid/compositional/functions/CompositionalProperties.hpp" #include "constitutive/fluid/multifluid/compositional/functions/CubicEOSPhaseModel.hpp" #include "TestFluid.hpp" +#include "TestFluidUtilities.hpp" namespace geos { namespace testing { -static constexpr real64 relTol = 1.0e-4; -static constexpr real64 absTol = 1.0e-7; - template< integer NC > using CompositionalPropertiesTestData = std::tuple< real64 const, // 0 - pressure @@ -64,12 +61,10 @@ class CompositionalPropertiesTestDataTestFixture : public ::testing::TestWithPar const auto [pressure, temperature, composition] = getInputData( data ); real64 const expectedMolarDensity = std::get< 3 >( data ); - real64 molarDensity = 0.0; - computeMolarDensity( pressure, - temperature, - composition, - molarDensity ); - checkRelativeError( molarDensity, expectedMolarDensity, relTol, absTol ); + real64 const molarDensity = computeMolarDensity( pressure, + temperature, + composition ); + checkRelativeError( molarDensity, expectedMolarDensity, internal::relTol, internal::absTol ); } // Compares the calculated molar density derivatives against numerical calculated @@ -79,8 +74,6 @@ class CompositionalPropertiesTestDataTestFixture : public ::testing::TestWithPar const auto [pressure, temperature, composition] = getInputData( data ); real64 molarDensity = 0.0; - real64 currentMolarDensity = 0.0; - real64 fdDerivative = 0.0; real64 dMolarDensity_dp = 0.0; real64 dMolarDensity_dt = 0.0; array1d< real64 > dMolarDensity_dz( numComps ); @@ -96,33 +89,32 @@ class CompositionalPropertiesTestDataTestFixture : public ::testing::TestWithPar // Compare against numerical derivatives // -- Pressure derivative real64 const dp = 1.0e-4 * pressure; - computeMolarDensity( pressure-dp, temperature, composition, currentMolarDensity ); - fdDerivative = -(currentMolarDensity - molarDensity) / dp; - checkDerivative( dMolarDensity_dp, fdDerivative, "Molar density left pressure derivative" ); - computeMolarDensity( pressure+dp, temperature, composition, currentMolarDensity ); - fdDerivative = (currentMolarDensity - molarDensity) / dp; - checkDerivative( dMolarDensity_dp, fdDerivative, "Molar density right pressure derivative" ); + internal::testNumericalDerivative( + pressure, dp, dMolarDensity_dp, + [this, & t=temperature, & zmf=composition]( real64 const p ) -> real64 { + return computeMolarDensity( p, t, zmf ); + } ); + // -- Temperature derivative - real64 const dt = 1.0e-6 * temperature; - computeMolarDensity( pressure, temperature-dt, composition, currentMolarDensity ); - fdDerivative = -(currentMolarDensity - molarDensity) / dt; - checkDerivative( dMolarDensity_dt, fdDerivative, "Molar density left temperature derivative" ); - computeMolarDensity( pressure, temperature+dt, composition, currentMolarDensity ); - fdDerivative = (currentMolarDensity - molarDensity) / dt; - checkDerivative( dMolarDensity_dt, fdDerivative, "Molar density right temperature derivative" ); + real64 const dT = 1.0e-6 * temperature; + internal::testNumericalDerivative( + temperature, dT, dMolarDensity_dt, + [this, & p=pressure, & zmf=composition]( real64 const t ) -> real64 { + return computeMolarDensity( p, t, zmf ); + } ); + // -- Composition derivatives derivative real64 const dz = 1.0e-7; for( integer ic = 0; ic < numComps; ++ic ) { - composition[ic] -= dz; - computeMolarDensity( pressure, temperature, composition, currentMolarDensity ); - fdDerivative = -(currentMolarDensity - molarDensity) / dz; - checkDerivative( dMolarDensity_dz[ic], fdDerivative, "Molar density left composition derivative" ); - composition[ic] += 2.0*dz; - computeMolarDensity( pressure, temperature, composition, currentMolarDensity ); - fdDerivative = (currentMolarDensity - molarDensity) / dz; - checkDerivative( dMolarDensity_dz[ic], fdDerivative, "Molar density right composition derivative" ); - composition[ic] -= dz; + internal::testNumericalDerivative( + 0.0, dz, dMolarDensity_dz[ic], + [this, & p=pressure, & t=temperature, zmf=composition, ic]( real64 const z ) -> real64 { + zmf[ic] += z; + real64 const density = computeMolarDensity( p, t, zmf ); + zmf[ic] -= z; + return density; + } ); } } @@ -132,12 +124,10 @@ class CompositionalPropertiesTestDataTestFixture : public ::testing::TestWithPar const auto [pressure, temperature, composition] = getInputData( data ); real64 const expectedMassDensity = std::get< 5 >( data ); - real64 massDensity = 0.0; - computeMassDensity( pressure, - temperature, - composition, - massDensity ); - checkRelativeError( massDensity, expectedMassDensity, relTol, absTol ); + real64 const massDensity = computeMassDensity( pressure, + temperature, + composition ); + checkRelativeError( massDensity, expectedMassDensity, internal::relTol, internal::absTol ); } // Compares the calculated mass density derivatives against numerical calculated @@ -147,8 +137,6 @@ class CompositionalPropertiesTestDataTestFixture : public ::testing::TestWithPar const auto [pressure, temperature, composition] = getInputData( data ); real64 massDensity = 0.0; - real64 currentMassDensity = 0.0; - real64 fdDerivative = 0.0; real64 dMassDensity_dp = 0.0; real64 dMassDensity_dt = 0.0; array1d< real64 > dMassDensity_dz( numComps ); @@ -164,33 +152,32 @@ class CompositionalPropertiesTestDataTestFixture : public ::testing::TestWithPar // Compare against numerical derivatives // -- Pressure derivative real64 const dp = 1.0e-4 * pressure; - computeMassDensity( pressure-dp, temperature, composition, currentMassDensity ); - fdDerivative = -(currentMassDensity - massDensity) / dp; - checkDerivative( dMassDensity_dp, fdDerivative, "Mass density left pressure derivative" ); - computeMassDensity( pressure+dp, temperature, composition, currentMassDensity ); - fdDerivative = (currentMassDensity - massDensity) / dp; - checkDerivative( dMassDensity_dp, fdDerivative, "Mass density right pressure derivative" ); + internal::testNumericalDerivative( + pressure, dp, dMassDensity_dp, + [this, & t=temperature, & zmf=composition]( real64 const p ) -> real64 { + return computeMassDensity( p, t, zmf ); + } ); + // -- Temperature derivative - real64 const dt = 1.0e-6 * temperature; - computeMassDensity( pressure, temperature-dt, composition, currentMassDensity ); - fdDerivative = -(currentMassDensity - massDensity) / dt; - checkDerivative( dMassDensity_dt, fdDerivative, "Mass density left temperature derivative" ); - computeMassDensity( pressure, temperature+dt, composition, currentMassDensity ); - fdDerivative = (currentMassDensity - massDensity) / dt; - checkDerivative( dMassDensity_dt, fdDerivative, "Mass density right temperature derivative" ); - // -- Composition derivatives derivative + real64 const dT = 1.0e-6 * temperature; + internal::testNumericalDerivative( + temperature, dT, dMassDensity_dt, + [this, & p=pressure, & zmf=composition]( real64 const t ) -> real64 { + return computeMassDensity( p, t, zmf ); + } ); + + // -- Composition derivatives real64 const dz = 1.0e-7; for( integer ic = 0; ic < numComps; ++ic ) { - composition[ic] -= dz; - computeMassDensity( pressure, temperature, composition, currentMassDensity ); - fdDerivative = -(currentMassDensity - massDensity) / dz; - checkDerivative( dMassDensity_dz[ic], fdDerivative, "Mass density left composition derivative" ); - composition[ic] += 2.0*dz; - computeMassDensity( pressure, temperature, composition, currentMassDensity ); - fdDerivative = (currentMassDensity - massDensity) / dz; - checkDerivative( dMassDensity_dz[ic], fdDerivative, "Mass density right composition derivative" ); - composition[ic] -= dz; + internal::testNumericalDerivative( + 0.0, dz, dMassDensity_dz[ic], + [this, & p=pressure, & t=temperature, & zmf=composition, ic]( real64 const z ) -> real64 { + zmf[ic] += z; + real64 const density = computeMassDensity( p, t, zmf ); + zmf[ic] -= z; + return density; + } ); } } @@ -205,14 +192,8 @@ class CompositionalPropertiesTestDataTestFixture : public ::testing::TestWithPar return {pressure, temperature, composition}; } - void checkDerivative( real64 const a, real64 const b, string const & name ) const - { - checkRelativeError( a, b, relTol, absTol, name ); - } - - void computeMolarDensity( real64 const pressure, real64 const temperature, - arrayView1d< real64 const > const & composition, - real64 & molarDensity ) const + real64 computeMolarDensity( real64 const pressure, real64 const temperature, + arrayView1d< real64 const > const & composition ) const { auto criticalPressure = m_fluid->getCriticalPressure(); auto criticalTemperature = m_fluid->getCriticalTemperature(); @@ -221,6 +202,7 @@ class CompositionalPropertiesTestDataTestFixture : public ::testing::TestWithPar real64 const binaryInteractionCoefficients = 0.0; real64 compressibilityFactor = 0.0; + real64 molarDensity = 0.0; array1d< real64 > aPureCoefficient( numComps ); array1d< real64 > bPureCoefficient( numComps ); real64 aMixtureCoefficient = 0.0; @@ -257,7 +239,9 @@ class CompositionalPropertiesTestDataTestFixture : public ::testing::TestWithPar volumeShift, compressibilityFactor, molarDensity ); + return molarDensity; } + void computeMolarDensity( real64 const pressure, real64 const temperature, arrayView1d< real64 const > const & composition, real64 & molarDensity, @@ -371,20 +355,21 @@ class CompositionalPropertiesTestDataTestFixture : public ::testing::TestWithPar dMolarDensity_dt, dMolarDensity_dz ); } - void computeMassDensity( real64 const pressure, real64 const temperature, - arrayView1d< real64 const > const & composition, - real64 & massDensity ) const + + real64 computeMassDensity( real64 const pressure, real64 const temperature, + arrayView1d< real64 const > const & composition ) const { auto molecularWeight = m_fluid->getMolecularWeight(); - real64 molarDensity = 0.0; - computeMolarDensity( pressure, temperature, composition, molarDensity ); + real64 const molarDensity = computeMolarDensity( pressure, temperature, composition ); + real64 massDensity = 0.0; constitutive::CompositionalProperties:: computeMassDensity( numComps, composition, molecularWeight, molarDensity, massDensity ); + return massDensity; } void computeMassDensity( real64 const pressure, real64 const temperature, @@ -395,8 +380,7 @@ class CompositionalPropertiesTestDataTestFixture : public ::testing::TestWithPar arraySlice1d< real64 > const dMassDensity_dz ) const { auto molecularWeight = m_fluid->getMolecularWeight(); - computeMassDensity( pressure, temperature, composition, - massDensity ); + massDensity = computeMassDensity( pressure, temperature, composition ); real64 molarDensity = 0.0; real64 dMolarDensity_dp = 0.0; diff --git a/src/coreComponents/constitutive/unitTests/testCubicEOS.cpp b/src/coreComponents/constitutive/unitTests/testCubicEOS.cpp index 4a1b08b3d9b..a2b2ffdf29c 100644 --- a/src/coreComponents/constitutive/unitTests/testCubicEOS.cpp +++ b/src/coreComponents/constitutive/unitTests/testCubicEOS.cpp @@ -13,10 +13,10 @@ */ // Source includes -#include "codingUtilities/UnitTestUtilities.hpp" #include "common/DataTypes.hpp" #include "constitutive/fluid/multifluid/compositional/functions/CubicEOSPhaseModel.hpp" #include "TestFluid.hpp" +#include "TestFluidUtilities.hpp" // TPL includes #include @@ -368,11 +368,6 @@ class DerivativeTestFixture : public ::testing::TestWithParam< TestData< NC > > ~DerivativeTestFixture() = default; protected: - void checkDerivative( real64 const a, real64 const b, string const & name ) const - { - checkRelativeError( a, b, relTol, absTol, name ); - } - std::unique_ptr< TestFluid< NC > > m_fluid; }; @@ -403,18 +398,14 @@ class MixCoeffDerivativeTestFixture : public DerivativeTestFixture< EOS, NC > public: void testNumericalDerivatives( ParamType const & testData ) const { - auto criticalPressure = this->m_fluid->getCriticalPressure(); - auto criticalTemperature = this->m_fluid->getCriticalTemperature(); - auto omega = this->m_fluid->getAcentricFactor(); - real64 binaryInteractionCoefficients = 0.0; // not implemented yet + auto const & fluid = *this->m_fluid; + auto criticalPressure = fluid.getCriticalPressure(); + auto criticalTemperature = fluid.getCriticalTemperature(); + auto omega = fluid.getAcentricFactor(); + real64 constexpr binaryInteractionCoefficients = 0.0; // not implemented yet array1d< real64 > aPureCoefficient( numComps ); array1d< real64 > bPureCoefficient( numComps ); - real64 aMixtureCoefficient = 0.0; - real64 bMixtureCoefficient = 0.0; - real64 currentAMixtureCoefficient = 0.0; - real64 currentBMixtureCoefficient = 0.0; - real64 fdDerivative = 0.0; real64 daMixtureCoefficient_dp = 0.0; real64 dbMixtureCoefficient_dp = 0.0; @@ -428,7 +419,9 @@ class MixCoeffDerivativeTestFixture : public DerivativeTestFixture< EOS, NC > real64 const temperature = std::get< 1 >( testData ); TestFluid< NC >::createArray( composition, std::get< 2 >( testData )); - auto computeCoefficients = [&]( real64 const p, real64 const t, auto const & zmf, real64 & a, real64 & b ){ + auto computeCoefficients = [&]( real64 const p, real64 const t, auto const & zmf ) -> std::pair< real64 const, real64 const > { + real64 a = 0.0; + real64 b = 0.0; CubicEOSPhaseModel< EOS >::computeMixtureCoefficients( numComps, p, t, zmf, @@ -438,10 +431,12 @@ class MixCoeffDerivativeTestFixture : public DerivativeTestFixture< EOS, NC > bPureCoefficient, a, b ); + return {a, b}; }; // Calculate values - computeCoefficients( pressure, temperature, composition, aMixtureCoefficient, bMixtureCoefficient ); + auto [aMixtureCoefficient, bMixtureCoefficient] = computeCoefficients( pressure, temperature, composition ); + // Calculate derivatives CubicEOSPhaseModel< EOS >::computeMixtureCoefficients( numComps, @@ -462,48 +457,54 @@ class MixCoeffDerivativeTestFixture : public DerivativeTestFixture< EOS, NC > dbMixtureCoefficient_dt, daMixtureCoefficient_dz, dbMixtureCoefficient_dz ); + // Compare against numerical derivatives // -- Pressure derivative real64 const dp = 1.0e-4 * pressure; - computeCoefficients( pressure-dp, temperature, composition, currentAMixtureCoefficient, currentBMixtureCoefficient ); - fdDerivative = -(currentAMixtureCoefficient - aMixtureCoefficient) / dp; - this->checkDerivative( daMixtureCoefficient_dp, fdDerivative, "Mixing Coeff A left pressure derivative" ); - fdDerivative = -(currentBMixtureCoefficient - bMixtureCoefficient) / dp; - this->checkDerivative( dbMixtureCoefficient_dp, fdDerivative, "Mixing Coeff B left pressure derivative" ); - computeCoefficients( pressure+dp, temperature, composition, currentAMixtureCoefficient, currentBMixtureCoefficient ); - fdDerivative = (currentAMixtureCoefficient - aMixtureCoefficient) / dp; - this->checkDerivative( daMixtureCoefficient_dp, fdDerivative, "Mixing Coeff A right pressure derivative" ); - fdDerivative = (currentBMixtureCoefficient - bMixtureCoefficient) / dp; - this->checkDerivative( dbMixtureCoefficient_dp, fdDerivative, "Mixing Coeff B right pressure derivative" ); + geos::testing::internal::testNumericalDerivative( + pressure, dp, daMixtureCoefficient_dp, + [&]( real64 const p ) -> real64 { + return computeCoefficients( p, temperature, composition ).first; + } ); + geos::testing::internal::testNumericalDerivative( + pressure, dp, dbMixtureCoefficient_dp, + [&]( real64 const p ) -> real64 { + return computeCoefficients( p, temperature, composition ).second; + } ); + // -- Temperature derivative - real64 const dt = 1.0e-6 * temperature; - computeCoefficients( pressure, temperature-dt, composition, currentAMixtureCoefficient, currentBMixtureCoefficient ); - fdDerivative = -(currentAMixtureCoefficient - aMixtureCoefficient) / dt; - this->checkDerivative( daMixtureCoefficient_dt, fdDerivative, "Mixing Coeff A left temperature derivative" ); - fdDerivative = -(currentBMixtureCoefficient - bMixtureCoefficient) / dt; - this->checkDerivative( dbMixtureCoefficient_dt, fdDerivative, "Mixing Coeff B left temperature derivative" ); - computeCoefficients( pressure, temperature+dt, composition, currentAMixtureCoefficient, currentBMixtureCoefficient ); - fdDerivative = (currentAMixtureCoefficient - aMixtureCoefficient) / dt; - this->checkDerivative( daMixtureCoefficient_dt, fdDerivative, "Mixing Coeff A right temperature derivative" ); - fdDerivative = (currentBMixtureCoefficient - bMixtureCoefficient) / dt; - this->checkDerivative( dbMixtureCoefficient_dt, fdDerivative, "Mixing Coeff B right temperature derivative" ); + real64 const dT = 1.0e-6 * temperature; + geos::testing::internal::testNumericalDerivative( + temperature, dT, daMixtureCoefficient_dt, + [&]( real64 const t ) -> real64 { + return computeCoefficients( pressure, t, composition ).first; + } ); + geos::testing::internal::testNumericalDerivative( + temperature, dT, dbMixtureCoefficient_dt, + [&]( real64 const t ) -> real64 { + return computeCoefficients( pressure, t, composition ).second; + } ); + // -- Composition derivatives derivative real64 const dz = 1.0e-7; for( integer ic = 0; ic < numComps; ++ic ) { - composition[ic] -= dz; - computeCoefficients( pressure, temperature, composition, currentAMixtureCoefficient, currentBMixtureCoefficient ); - fdDerivative = -(currentAMixtureCoefficient - aMixtureCoefficient) / dz; - this->checkDerivative( daMixtureCoefficient_dz[ic], fdDerivative, "Mixing Coeff A left composition derivative" ); - fdDerivative = -(currentBMixtureCoefficient - bMixtureCoefficient) / dz; - this->checkDerivative( dbMixtureCoefficient_dz[ic], fdDerivative, "Mixing Coeff B left composition derivative" ); - composition[ic] += 2.0*dz; - computeCoefficients( pressure, temperature, composition, currentAMixtureCoefficient, currentBMixtureCoefficient ); - fdDerivative = (currentAMixtureCoefficient - aMixtureCoefficient) / dz; - this->checkDerivative( daMixtureCoefficient_dz[ic], fdDerivative, "Mixing Coeff A right composition derivative" ); - fdDerivative = (currentBMixtureCoefficient - bMixtureCoefficient) / dz; - this->checkDerivative( dbMixtureCoefficient_dz[ic], fdDerivative, "Mixing Coeff B right composition derivative" ); - composition[ic] -= dz; + auto computeComponentCoefficients = [&]( real64 const z ) { + composition[ic] += z; + auto const coefficients = computeCoefficients( pressure, temperature, composition ); + composition[ic] -= z; + return coefficients; + }; + geos::testing::internal::testNumericalDerivative( + 0.0, dz, daMixtureCoefficient_dz[ic], + [&]( real64 const z ) -> real64 { + return computeComponentCoefficients( z ).first; + } ); + geos::testing::internal::testNumericalDerivative( + 0.0, dz, dbMixtureCoefficient_dz[ic], + [&]( real64 const z ) -> real64 { + return computeComponentCoefficients( z ).second; + } ); } } }; @@ -564,10 +565,11 @@ class CompressibilityDerivativeTestFixture : public DerivativeTestFixture< EOS, public: void testNumericalDerivatives( ParamType const & testData ) const { - auto criticalPressure = this->m_fluid->getCriticalPressure(); - auto criticalTemperature = this->m_fluid->getCriticalTemperature(); - auto omega = this->m_fluid->getAcentricFactor(); - real64 binaryInteractionCoefficients = 0.0; // not implemented yet + auto const & fluid = *this->m_fluid; + auto criticalPressure = fluid.getCriticalPressure(); + auto criticalTemperature = fluid.getCriticalTemperature(); + auto omega = fluid.getAcentricFactor(); + real64 constexpr binaryInteractionCoefficients = 0.0; // not implemented yet array1d< real64 > aPureCoefficient( numComps ); array1d< real64 > bPureCoefficient( numComps ); @@ -580,20 +582,17 @@ class CompressibilityDerivativeTestFixture : public DerivativeTestFixture< EOS, array1d< real64 > daMixtureCoefficient_dz( numComps ); array1d< real64 > dbMixtureCoefficient_dz( numComps ); - real64 compressibilityFactor = 0.0; real64 dCompressibilityFactor_dp = 0.0; real64 dCompressibilityFactor_dt = 0.0; array1d< real64 > dCompressibilityFactor_dz( numComps ); - real64 currentCompressibilityFactor = 0.0; - real64 fdDerivative = 0.0; - array1d< real64 > composition; real64 const pressure = std::get< 0 >( testData ); real64 const temperature = std::get< 1 >( testData ); TestFluid< NC >::createArray( composition, std::get< 2 >( testData )); - auto computeCompressibilityFactor = [&]( real64 const p, real64 const t, auto const & zmf, real64 & z ){ + auto computeCompressibilityFactor = [&]( real64 const p, real64 const t, auto const & zmf ) -> real64 { + real64 z = 0.0; CubicEOSPhaseModel< EOS >::computeMixtureCoefficients( numComps, p, t, zmf, @@ -612,10 +611,12 @@ class CompressibilityDerivativeTestFixture : public DerivativeTestFixture< EOS, aMixtureCoefficient, bMixtureCoefficient, z ); + return z; }; // Calculate values - computeCompressibilityFactor( pressure, temperature, composition, compressibilityFactor ); + real64 const compressibilityFactor = computeCompressibilityFactor( pressure, temperature, composition ); + // Calculate derivatives CubicEOSPhaseModel< EOS >::computeMixtureCoefficients( numComps, @@ -650,36 +651,36 @@ class CompressibilityDerivativeTestFixture : public DerivativeTestFixture< EOS, dCompressibilityFactor_dp, dCompressibilityFactor_dt, dCompressibilityFactor_dz ); + // Compare against numerical derivatives // -- Pressure derivative real64 const dp = 1.0e-4 * pressure; - computeCompressibilityFactor( pressure-dp, temperature, composition, currentCompressibilityFactor ); - fdDerivative = -(currentCompressibilityFactor - compressibilityFactor) / dp; - this->checkDerivative( dCompressibilityFactor_dp, fdDerivative, "Compressibility factor left pressure derivative" ); - computeCompressibilityFactor( pressure+dp, temperature, composition, currentCompressibilityFactor ); - fdDerivative = (currentCompressibilityFactor - compressibilityFactor) / dp; - this->checkDerivative( dCompressibilityFactor_dp, fdDerivative, "Compressibility factor right pressure derivative" ); + geos::testing::internal::testNumericalDerivative( + pressure, dp, dCompressibilityFactor_dp, + [&]( real64 const p ) -> real64 { + return computeCompressibilityFactor( p, temperature, composition ); + } ); + // -- Temperature derivative - real64 const dt = 1.0e-6 * temperature; - computeCompressibilityFactor( pressure, temperature-dt, composition, currentCompressibilityFactor ); - fdDerivative = -(currentCompressibilityFactor - compressibilityFactor) / dt; - this->checkDerivative( dCompressibilityFactor_dt, fdDerivative, "Compressibility factor left temperature derivative" ); - computeCompressibilityFactor( pressure, temperature+dt, composition, currentCompressibilityFactor ); - fdDerivative = (currentCompressibilityFactor - compressibilityFactor) / dt; - this->checkDerivative( dCompressibilityFactor_dt, fdDerivative, "Compressibility factor right temperature derivative" ); + real64 const dT = 1.0e-6 * temperature; + geos::testing::internal::testNumericalDerivative( + temperature, dT, dCompressibilityFactor_dt, + [&]( real64 const t ) -> real64 { + return computeCompressibilityFactor( pressure, t, composition ); + } ); + // -- Composition derivatives derivative real64 const dz = 1.0e-7; for( integer ic = 0; ic < numComps; ++ic ) { - composition[ic] -= dz; - computeCompressibilityFactor( pressure, temperature, composition, currentCompressibilityFactor ); - fdDerivative = -(currentCompressibilityFactor - compressibilityFactor) / dz; - this->checkDerivative( dCompressibilityFactor_dz[ic], fdDerivative, "Compressibility factor left composition derivative" ); - composition[ic] += 2.0*dz; - computeCompressibilityFactor( pressure, temperature, composition, currentCompressibilityFactor ); - fdDerivative = (currentCompressibilityFactor - compressibilityFactor) / dz; - this->checkDerivative( dCompressibilityFactor_dz[ic], fdDerivative, "Compressibility factor right composition derivative" ); - composition[ic] -= dz; + geos::testing::internal::testNumericalDerivative( + 0.0, dz, dCompressibilityFactor_dz[ic], + [&]( real64 const z ) -> real64 { + composition[ic] += z; + real64 const compressibility = computeCompressibilityFactor( pressure, temperature, composition ); + composition[ic] -= z; + return compressibility; + } ); } } };