Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Refactor fluid unit tests #2695

Merged
merged 6 commits into from
Oct 28, 2023
Merged

Refactor fluid unit tests #2695

Show file tree
Hide file tree
Changes from 3 commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
0b39ae4
Refactor fluid unit tests
dkachuma Sep 21, 2023
98155bc
Merge branch 'develop' into feature/dkachuma/fluid-unit-test-derivatives
dkachuma Sep 22, 2023
8814a7c
Explicitly name lambda captures
dkachuma Sep 22, 2023
8490cd1
Merge branch 'develop' into feature/dkachuma/fluid-unit-test-derivatives
dkachuma Oct 3, 2023
775065e
Change temperature variable name
dkachuma Oct 17, 2023
93875e2
Merge branch 'develop' into feature/dkachuma/fluid-unit-test-derivatives
dkachuma Oct 17, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
68 changes: 68 additions & 0 deletions src/coreComponents/constitutive/unitTests/TestFluidUtilities.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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_
140 changes: 62 additions & 78 deletions src/coreComponents/constitutive/unitTests/testCompositionalProperties.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand All @@ -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 );
Expand All @@ -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;
dkachuma marked this conversation as resolved.
Show resolved Hide resolved
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" );
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;
} );
}
}

Expand All @@ -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
Expand All @@ -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 );
Expand All @@ -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
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;
} );
}
}

Expand All @@ -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();
Expand All @@ -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;
Expand Down Expand Up @@ -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,
Expand Down Expand Up @@ -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,
Expand All @@ -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;
Expand Down
Loading