diff --git a/.integrated_tests.yaml b/.integrated_tests.yaml index e1ea844bc05..d729374a0b1 100644 --- a/.integrated_tests.yaml +++ b/.integrated_tests.yaml @@ -1,6 +1,6 @@ baselines: bucket: geosx - baseline: integratedTests/baseline_integratedTests-pr4062-16784-6d8782e + baseline: integratedTests/baseline_integratedTests-pr3705-16862-2b262bf allow_fail: all: '' diff --git a/BASELINE_NOTES.md b/BASELINE_NOTES.md index 2e69f9d42db..fd5b84453a3 100644 --- a/BASELINE_NOTES.md +++ b/BASELINE_NOTES.md @@ -5,6 +5,9 @@ This file is designed to track changes to the integrated test baselines. Any developer who updates the baseline ID in the .integrated_tests.yaml file is expected to create an entry in this file with the pull request number, date, and their justification for rebaselining. These notes should be in reverse-chronological order, and use the following time format: (YYYY-MM-DD). +PR #3705 (2026-06-03) +Implement compositional enthalpy model + PR #4062 (2026-05-26) Add Porous Solid other than PorousElasticity for ALM solver diff --git a/src/coreComponents/constitutive/CMakeLists.txt b/src/coreComponents/constitutive/CMakeLists.txt index 88fd5217b24..7606cfe25f6 100644 --- a/src/coreComponents/constitutive/CMakeLists.txt +++ b/src/coreComponents/constitutive/CMakeLists.txt @@ -102,6 +102,7 @@ set( constitutive_headers fluid/multifluid/compositional/functions/StabilityTest.hpp fluid/multifluid/compositional/functions/Utilities.hpp fluid/multifluid/compositional/models/CompositionalDensity.hpp + fluid/multifluid/compositional/models/CompositionalEnthalpy.hpp fluid/multifluid/compositional/models/ConstantViscosity.hpp fluid/multifluid/compositional/models/FunctionBase.hpp fluid/multifluid/compositional/models/ImmiscibleWaterDensity.hpp @@ -121,6 +122,7 @@ set( constitutive_headers fluid/multifluid/compositional/parameters/CriticalVolume.hpp fluid/multifluid/compositional/parameters/EquationOfState.hpp fluid/multifluid/compositional/parameters/FlashParameters.hpp + fluid/multifluid/compositional/parameters/HeatCapacityCoefficients.hpp fluid/multifluid/compositional/parameters/ImmiscibleWaterParameters.hpp fluid/multifluid/compositional/parameters/KValueFlashParameters.hpp fluid/multifluid/compositional/parameters/ModelParameters.hpp @@ -267,6 +269,7 @@ set( constitutive_sources fluid/multifluid/CO2Brine/functions/PureWaterProperties.cpp fluid/multifluid/CO2Brine/functions/WaterDensity.cpp fluid/multifluid/compositional/models/CompositionalDensity.cpp + fluid/multifluid/compositional/models/CompositionalEnthalpy.cpp fluid/multifluid/compositional/models/ConstantViscosity.cpp fluid/multifluid/compositional/models/ImmiscibleWaterDensity.cpp fluid/multifluid/compositional/models/ImmiscibleWaterFlashModel.cpp @@ -280,6 +283,7 @@ set( constitutive_sources fluid/multifluid/compositional/parameters/ComponentProperties.cpp fluid/multifluid/compositional/parameters/CriticalVolume.cpp fluid/multifluid/compositional/parameters/FlashParameters.cpp + fluid/multifluid/compositional/parameters/HeatCapacityCoefficients.cpp fluid/multifluid/compositional/parameters/ImmiscibleWaterParameters.cpp fluid/multifluid/compositional/parameters/KValueFlashParameters.cpp fluid/multifluid/compositional/parameters/PhaseType.cpp diff --git a/src/coreComponents/constitutive/fluid/multifluid/compositional/functions/CubicEOSPhaseModel.hpp b/src/coreComponents/constitutive/fluid/multifluid/compositional/functions/CubicEOSPhaseModel.hpp index 5728941198a..3c6bae25567 100644 --- a/src/coreComponents/constitutive/fluid/multifluid/compositional/functions/CubicEOSPhaseModel.hpp +++ b/src/coreComponents/constitutive/fluid/multifluid/compositional/functions/CubicEOSPhaseModel.hpp @@ -279,6 +279,38 @@ struct CubicEOSPhaseModel arraySlice1d< real64 > const & logFugacityCoefficients, StackDerivativeType< 2, DERIVATIVES, USD2 > const & logFugacityCoefficientDerivs ); + /** + * @brief Compute the dimensionless enthalpy departure + * @details Computes the dimensionless enthalpy departure (H_dep / RT) alongside its analytical + * derivatives with respect to pressure, temperature, and compositions. Pressure, temperature, + * and componentProperties are required parameters as calculating mixture parameter first and + * second temperature derivatives dictates access to properties like critical Temperature. + * @tparam USD Length descriptor for un-strided array access + * @tparam DERIVATIVES a flag to indicate if derivatives w.r.t p, t, and composition should be evaluated + * @param[in] numComps number of components + * @param[in] pressure pressure of the phase + * @param[in] temperature temperature of the phase + * @param[in] composition composition of the phase + * @param[in] componentProperties The compositional component properties + * @param[in] data The component mixture properties evaluated into an EOSStackVariables object. This should + * always have the derivatives because enthalpy requires the temperature derivative of Z. + * @param[out] enthalpy Calculated dimensionless enthalpy departure H_dep + * @param[out] enthalpyDerivs Slice to hold analytical derivatives H_dep' evaluated if DERIVATIVES is true + * @param[in] selectedRoot Explicit indicator to track designated root selection (AUTO by default) + */ + template< integer USD, bool DERIVATIVES = false > + GEOS_HOST_DEVICE + static void + computeEnthalpy( integer const numComps, + real64 const & pressure, + real64 const & temperature, + arraySlice1d< real64 const, USD > const & composition, + ComponentProperties::KernelWrapper const & componentProperties, + StackVariables< true > const & data, + real64 & enthalpy, + StackDerivativeType< 1, DERIVATIVES > const & enthalpyDerivs, + SelectedRoot const selectedRoot = SelectedRoot::AUTO ); + /** * @brief Helper functions solving a cubic equation using trigonometry * m3 * x^3 + m2 * x^2 + m1 *x + m0 = 0 diff --git a/src/coreComponents/constitutive/fluid/multifluid/compositional/functions/CubicEOSPhaseModel_impl.hpp b/src/coreComponents/constitutive/fluid/multifluid/compositional/functions/CubicEOSPhaseModel_impl.hpp index 07feb4d5a7a..4195e62a9dc 100644 --- a/src/coreComponents/constitutive/fluid/multifluid/compositional/functions/CubicEOSPhaseModel_impl.hpp +++ b/src/coreComponents/constitutive/fluid/multifluid/compositional/functions/CubicEOSPhaseModel_impl.hpp @@ -316,6 +316,8 @@ computeMixtureCoefficients( integer const numComps, { LvArray::forValuesInSlice( stack.daMixture, setZero ); LvArray::forValuesInSlice( stack.dbMixture, setZero ); + LvArray::forValuesInSlice( stack.ddaMixture, setZero ); + for( integer ic = 0; ic < numComps; ++ic ) { for( integer jc = 0; jc < numComps; ++jc ) @@ -323,16 +325,32 @@ computeMixtureCoefficients( integer const numComps, real64 const sqrt_aiaj = LvArray::math::sqrt( stack.aic[ic] * stack.aic[jc] ); real64 const kij_term = 1.0 - kij( ic, jc ); real64 const aij = kij_term * sqrt_aiaj; - real64 const coeff = 0.5 * kij_term / sqrt_aiaj; - real64 const daij_dT = coeff * (stack.aic[jc] * stack.daic_dt[ic] + stack.aic[ic] * stack.daic_dt[jc]); + real64 const f_T = 0.5 / sqrt_aiaj * (stack.aic[jc] * stack.daic_dt[ic] + stack.aic[ic] * stack.daic_dt[jc]); + real64 const f_P = 0.5 / sqrt_aiaj * (stack.aic[jc] * stack.daic_dp[ic] + stack.aic[ic] * stack.daic_dp[jc]); - real64 const daij_dP = coeff * (stack.aic[jc] * stack.daic_dp[ic] + stack.aic[ic] * stack.daic_dp[jc]); + real64 const daij_dT = kij_term * f_T; + real64 const daij_dP = kij_term * f_P; stack.daMixture[Deriv::dP] += composition[ic] * composition[jc] * daij_dP; stack.daMixture[Deriv::dT] += composition[ic] * composition[jc] * daij_dT; stack.daMixture[Deriv::dC+ic] += composition[jc] * aij; stack.daMixture[Deriv::dC+jc] += composition[ic] * aij; + + // Computation for ddaMixture (Derivatives of daMixture[Deriv::dT] wrt to all variables) + stack.ddaMixture[Deriv::dC+ic] += composition[jc] * daij_dT; + stack.ddaMixture[Deriv::dC+jc] += composition[ic] * daij_dT; + + real64 const d2aij_dTdP = kij_term * (f_T * f_P / sqrt_aiaj); + stack.ddaMixture[Deriv::dP] += composition[ic] * composition[jc] * d2aij_dTdP; + + real64 const f_TT = 0.5 / sqrt_aiaj * ( + stack.aic[jc] * stack.d2aic_dt2[ic] + + 2.0 * stack.daic_dt[ic] * stack.daic_dt[jc] + + stack.aic[ic] * stack.d2aic_dt2[jc] + ) - f_T * f_T / sqrt_aiaj; + real64 const d2aij_dT2 = kij_term * f_TT; + stack.ddaMixture[Deriv::dT] += composition[ic] * composition[jc] * d2aij_dT2; } stack.dbMixture[Deriv::dP] += composition[ic] * stack.dbic_dp[ic]; stack.dbMixture[Deriv::dT] += composition[ic] * stack.dbic_dt[ic]; @@ -347,6 +365,16 @@ computeMixtureCoefficients( integer const numComps, real64 const sqrt_aiaj = LvArray::math::sqrt( stack.aic[ic] * stack.aic[jc] ); real64 const dkij_term_dT = -stack.dkij_dT( ic, jc ); stack.daMixture[Deriv::dT] += composition[ic] * composition[jc] * dkij_term_dT * sqrt_aiaj; + + // Updating ddaMixture for the temperature-dependent BICs part + stack.ddaMixture[Deriv::dC+ic] += composition[jc] * dkij_term_dT * sqrt_aiaj; + stack.ddaMixture[Deriv::dC+jc] += composition[ic] * dkij_term_dT * sqrt_aiaj; + + real64 const f_P = 0.5 / sqrt_aiaj * (stack.aic[jc] * stack.daic_dp[ic] + stack.aic[ic] * stack.daic_dp[jc]); + stack.ddaMixture[Deriv::dP] += composition[ic] * composition[jc] * dkij_term_dT * f_P; + + real64 const f_T = 0.5 / sqrt_aiaj * (stack.aic[jc] * stack.daic_dt[ic] + stack.aic[ic] * stack.daic_dt[jc]); + stack.ddaMixture[Deriv::dT] += composition[ic] * composition[jc] * 2.0 * dkij_term_dT * f_T; } } } @@ -608,6 +636,91 @@ computeLogFugacityCoefficients( integer const numComps, } } +template< typename EOS_TYPE > +template< integer USD, bool DERIVATIVES > +GEOS_HOST_DEVICE +void +CubicEOSPhaseModel< EOS_TYPE >:: +computeEnthalpy( integer const numComps, + real64 const & pressure, + real64 const & temperature, + arraySlice1d< real64 const, USD > const & composition, + ComponentProperties::KernelWrapper const & componentProperties, + StackVariables< true > const & data, + real64 & enthalpy, + StackDerivativeType< 1, DERIVATIVES > const & enthalpyDerivs, + SelectedRoot const selectedRoot ) +{ + GEOS_UNUSED_VAR( pressure ); + GEOS_UNUSED_VAR( componentProperties ); + + real64 Z = 0.0; + auto const & compressibilityDerivs = enthalpyDerivs; + + // Dynamically resolve compressibility derivative generation based on the DERIVATIVES boolean flag + computeCompressibilityFactor< USD, DERIVATIVES >( numComps, + composition, + data, + Z, + compressibilityDerivs, + selectedRoot ); + + real64 const A = data.aMixture; + real64 const B = data.bMixture; + real64 const dA_dT = data.daMixture[Deriv::dT]; + + // M represents generalized term inside the logarithm coefficient bracket: A + T * (dA/dT) + real64 const M = A + temperature * dA_dT; + + real64 const expE = ( Z + EOS_TYPE::delta1 * B ) / ( Z + EOS_TYPE::delta2 * B ); +#if defined(GEOS_DEVICE_COMPILE) + GEOS_ERROR_IF( expE < MultiFluidConstants::epsilon, + "Cubic EOS Enthalpy failed: exp(E) is below epsilon." ); +#else + GEOS_ERROR_IF( expE < MultiFluidConstants::epsilon, + GEOS_FMT( "Cubic EOS Enthalpy failed with exp(E)={}", expE ) ); +#endif + real64 const E = log( expE ); + real64 const G = 1.0 / ( ( EOS_TYPE::delta1 - EOS_TYPE::delta2 ) * B ); + + // Dimensionless enthalpy departure: H_dep / RT = Z - 1.0 + (A + T*A_T) * G * E + real64 const dimensionlessEnthalpy = Z - 1.0 + M * G * E; + + // Return scaled dimensional enthalpy H_dep = R * T * (H_dep / RT) + real64 const RT = constants::gasConstant * temperature; + enthalpy = RT * dimensionlessEnthalpy; + + if constexpr ( DERIVATIVES ) + { + integer const numDofs = numComps + 2; + for( integer idof = 0; idof < numDofs; ++idof ) + { + real64 const dZ = compressibilityDerivs[idof]; + real64 const dA = data.daMixture[idof]; + real64 const dB = data.dbMixture[idof]; + real64 const ddA = data.ddaMixture[idof]; + + // Derivatives of M using standard product rules + real64 dM = dA + temperature * ddA; + if( idof == Deriv::dT ) + { + dM += dA_dT; + } + + real64 const dE = ( dZ + EOS_TYPE::delta1 * dB ) / ( Z + EOS_TYPE::delta1 * B ) - + ( dZ + EOS_TYPE::delta2 * dB ) / ( Z + EOS_TYPE::delta2 * B ); + + real64 const dG = -G * dB / B; + + real64 const dDimensionlessEnthalpy = dZ + ( dM * G * E + M * dG * E + M * G * dE ); + + // Apply product rule to dimensional scaling: d(RT * H_dimless) + enthalpyDerivs[idof] = RT * dDimensionlessEnthalpy; + } + enthalpyDerivs[Deriv::dT] += constants::gasConstant * dimensionlessEnthalpy; + } +} + template< typename EOS_TYPE > GEOS_HOST_DEVICE void diff --git a/src/coreComponents/constitutive/fluid/multifluid/compositional/functions/EOSStackVariables.hpp b/src/coreComponents/constitutive/fluid/multifluid/compositional/functions/EOSStackVariables.hpp index 80587bb972a..b92d566d0e1 100644 --- a/src/coreComponents/constitutive/fluid/multifluid/compositional/functions/EOSStackVariables.hpp +++ b/src/coreComponents/constitutive/fluid/multifluid/compositional/functions/EOSStackVariables.hpp @@ -145,7 +145,7 @@ struct EOSStackVariables_Impl< T, true > : public EOSStackVariables_Impl< T, fal arraySlice2d< real64 const > const dbip_dT ): EOSStackVariables_Impl< T, false >( numComps, bip ), dkij_dT( dbip_dT ), - m_derivativeData( 8, numComps+2 ), + m_derivativeData( 9, numComps+2 ), daic_dp( m_derivativeData[0] ), dbic_dp( m_derivativeData[1] ), daic_dt( m_derivativeData[2] ), @@ -153,14 +153,15 @@ struct EOSStackVariables_Impl< T, true > : public EOSStackVariables_Impl< T, fal d2aic_dt2( m_derivativeData[4] ), d2bic_dt2( m_derivativeData[5] ), daMixture( m_derivativeData[6] ), - dbMixture( m_derivativeData[7] ) + dbMixture( m_derivativeData[7] ), + ddaMixture( m_derivativeData[8] ) {} /// Temperature derivatives of binary interaction parameters (dk_ij/dT). arraySlice2d< real64 const > const dkij_dT; /// Internal storage for all derivative-related temporary data. - StackArray< real64, 2, 8 * maxNumDof > m_derivativeData; + StackArray< real64, 2, 9 * maxNumDof > m_derivativeData; /// Pressure derivatives of component-specific 'a_i' coefficients (da_i/dp). arraySlice1d< real64 > const daic_dp; @@ -185,6 +186,9 @@ struct EOSStackVariables_Impl< T, true > : public EOSStackVariables_Impl< T, fal /// Derivatives of mixture parameter 'b' with respect to all degrees of freedom. DerivativeType<> const dbMixture; + + /// Derivatives of temperature derivatives of mixture parameter 'a' with respect to all degrees of freedom. + DerivativeType<> const ddaMixture; }; /** diff --git a/src/coreComponents/constitutive/fluid/multifluid/compositional/models/CompositionalEnthalpy.cpp b/src/coreComponents/constitutive/fluid/multifluid/compositional/models/CompositionalEnthalpy.cpp new file mode 100644 index 00000000000..96e51de6ba1 --- /dev/null +++ b/src/coreComponents/constitutive/fluid/multifluid/compositional/models/CompositionalEnthalpy.cpp @@ -0,0 +1,93 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2024 TotalEnergies + * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2023-2024 Chevron + * Copyright (c) 2019- GEOS/GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +/** + * @file CompositionalEnthalpy.cpp + */ + +#include "CompositionalEnthalpy.hpp" +#include "constitutive/fluid/multifluid/compositional/parameters/HeatCapacityCoefficients.hpp" +#include "common/PhysicsConstants.hpp" + +namespace geos +{ +namespace constitutive +{ +namespace compositional +{ + +CompositionalEnthalpyUpdate::CompositionalEnthalpyUpdate( integer const phaseIndex, + integer const vapourIndex, + PhaseType const phaseType, + EquationOfStateType const equationOfState, + real64 const refTemperature, + arrayView2d< real64 const > const & referenceEnthalpy, + arrayView2d< real64 const > const & coefficients ): + m_phaseIndex( phaseIndex ), + m_vapourIndex( vapourIndex ), + m_phaseType( phaseType ), + m_equationOfState( equationOfState ), + m_refTemperature( refTemperature ), + m_referenceEnthalpy( referenceEnthalpy ), + m_coefficients( coefficients ) +{} + +CompositionalEnthalpy::CompositionalEnthalpy( string const & name, + ComponentProperties const & componentProperties, + integer const phaseIndex, + ModelParameters const & modelParameters ): + FunctionBase( name, componentProperties ), + m_phaseIndex( phaseIndex ), + m_heatCapacityCoefficients ( modelParameters.get< HeatCapacityCoefficients >() ) +{ + EquationOfState const * equationOfState = modelParameters.get< EquationOfState >(); + string const eosName = equationOfState->m_equationsOfStateNames[phaseIndex]; + m_equationOfState = EnumStrings< EquationOfStateType >::fromString( eosName ); + + integer const numPhases = m_heatCapacityCoefficients->m_phaseTypes.size(); + + m_phaseType = m_heatCapacityCoefficients->m_phaseTypes[phaseIndex]; + for( integer ip = 0; ip < numPhases; ip++ ) + { + if( m_heatCapacityCoefficients->m_phaseTypes[ip] == PhaseType::VAPOUR ) + { + m_vapourIndex = ip; + } + } +} + +CompositionalEnthalpy::KernelWrapper +CompositionalEnthalpy::createKernelWrapper() const +{ + return KernelWrapper( m_phaseIndex, + m_vapourIndex, + m_phaseType, + m_equationOfState, + m_heatCapacityCoefficients->m_referenceTemperature, + m_heatCapacityCoefficients->m_referenceEnthalpy.toViewConst(), + m_heatCapacityCoefficients->m_coefficients.toViewConst() ); +} + +std::unique_ptr< ModelParameters > +CompositionalEnthalpy::createParameters( std::unique_ptr< ModelParameters > parameters ) +{ + auto params = EquationOfState::create( std::move( parameters ) ); + params = HeatCapacityCoefficients::create( std::move( params ) ); + return params; +} + +} // namespace compositional +} // namespace constitutive +} // namespace geos diff --git a/src/coreComponents/constitutive/fluid/multifluid/compositional/models/CompositionalEnthalpy.hpp b/src/coreComponents/constitutive/fluid/multifluid/compositional/models/CompositionalEnthalpy.hpp new file mode 100644 index 00000000000..24ac6fe88c8 --- /dev/null +++ b/src/coreComponents/constitutive/fluid/multifluid/compositional/models/CompositionalEnthalpy.hpp @@ -0,0 +1,305 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2024 TotalEnergies + * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2023-2024 Chevron + * Copyright (c) 2019- GEOS/GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +/** + * @file CompositionalEnthalpy.hpp + */ + +#ifndef GEOS_CONSTITUTIVE_FLUID_MULTIFLUID_COMPOSITIONAL_MODELS_COMPOSITIONALENTHALPY_HPP_ +#define GEOS_CONSTITUTIVE_FLUID_MULTIFLUID_COMPOSITIONAL_MODELS_COMPOSITIONALENTHALPY_HPP_ + +#include "FunctionBase.hpp" +#include "constitutive/fluid/multifluid/compositional/parameters/EquationOfState.hpp" +#include "constitutive/fluid/multifluid/compositional/parameters/PhaseType.hpp" + +#include "constitutive/fluid/multifluid/MultiFluidUtils.hpp" +#include "constitutive/fluid/multifluid/MultiFluidConstants.hpp" +#include "constitutive/fluid/multifluid/compositional/functions/CompositionalProperties.hpp" +#include "constitutive/fluid/multifluid/compositional/functions/CubicEOSPhaseModel.hpp" +#include "constitutive/fluid/multifluid/compositional/functions/SoreideWhitsonEOSPhaseModel.hpp" + +namespace geos +{ + +namespace constitutive +{ + +namespace compositional +{ + +class HeatCapacityCoefficients; + +class CompositionalEnthalpyUpdate final : public FunctionBaseUpdate +{ + using Deriv = multifluid::DerivativeOffset; +public: + CompositionalEnthalpyUpdate( integer const phaseIndex, + integer const vapourIndex, + PhaseType const phaseType, + EquationOfStateType const equationOfState, + real64 const refTemperature, + arrayView2d< real64 const > const & referenceEnthalpy, + arrayView2d< real64 const > const & coefficients ); + + template< integer USD1, integer USD2 > + GEOS_HOST_DEVICE + void compute( ComponentProperties::KernelWrapper const & componentProperties, + real64 const & pressure, + real64 const & temperature, + arraySlice1d< real64 const, USD1 > const & phaseComposition, + real64 & enthalpy, + arraySlice1d< real64, USD2 > const & dEnthalpy, + bool useMass ) const; + + /** + * @brief Evaluates the Poling polynomial at a given temperature given a list of coefficients + * @param[in] dT - the temperature difference between the temperature and the reference + * @param[in] a - the coefficients + * @param[out] enthalpy - the enthalpy + * @param[out] heatCapacity - the enthalpy derivative wrt temperature + */ + GEOS_FORCE_INLINE + GEOS_HOST_DEVICE + static void evaluatePolynomial( real64 const & dT, + arraySlice1d< real64 const > const & a, + real64 & enthalpy, + real64 & heatCapacity ) + { + constexpr real64 r2 = 1.0 / 2.0; + constexpr real64 r3 = 1.0 / 3.0; + constexpr real64 r4 = 1.0 / 4.0; + constexpr real64 r5 = 1.0 / 5.0; + + // c(T) = a[0] + a[1]*dT + a[2]*dT^2 + a[3]*dT^3 + a[4]*dT^4 + heatCapacity = a[0] + dT * (a[1] + dT * (a[2] + dT * (a[3] + dT * a[4]))); + + // Evaluate enthalpy using the integral of c(s) ds from Tref to T + // H(T) = a[0]*dT + (a[1]/2)*dT^2 + (a[2]/3)*dT^3 + (a[3]/4)*dT^4 + (a[4]/5)*dT^5 + enthalpy = dT * (a[0] + dT * (a[1] * r2 + dT * (a[2] * r3 + dT * (a[3] * r4 + dT * (a[4] * r5))))); + } + + template< typename LAMBDA > + GEOS_HOST_DEVICE + static void selectEquationOfState( EquationOfStateType equationOfState, LAMBDA && lambda ); + + template< typename EOS, integer USD1, integer USD2 > + GEOS_HOST_DEVICE + void calculateEquationOfStateEnthalpy( integer const numComps, + ComponentProperties::KernelWrapper const & componentProperties, + real64 const & pressure, + real64 const & temperature, + arraySlice1d< real64 const, USD1 > const & phaseComposition, + real64 & enthalpy, + arraySlice1d< real64, USD2 > const & dEnthalpy, + PhaseType const & phaseType ) const; + +private: + integer const m_phaseIndex{-1}; + integer const m_vapourIndex{-1}; + PhaseType const m_phaseType; + EquationOfStateType const m_equationOfState; + real64 const m_refTemperature; + arrayView2d< real64 const > const m_referenceEnthalpy; + arrayView2d< real64 const > const m_coefficients; +}; + +class CompositionalEnthalpy : public FunctionBase +{ +public: + CompositionalEnthalpy( string const & name, + ComponentProperties const & componentProperties, + integer const phaseIndex, + ModelParameters const & modelParameters ); + + static string catalogName() { return "CompositionalEnthalpy"; } + + virtual FunctionType functionType() const override + { + return FunctionType::ENTHALPY; + } + + /// Type of kernel wrapper for in-kernel update + using KernelWrapper = CompositionalEnthalpyUpdate; + + /** + * @brief Create an update kernel wrapper. + * @return the wrapper + */ + KernelWrapper createKernelWrapper() const; + + // Create parameters unique to this model + static std::unique_ptr< ModelParameters > createParameters( std::unique_ptr< ModelParameters > parameters ); + +private: + integer m_phaseIndex{-1}; + integer m_vapourIndex{-1}; + PhaseType m_phaseType; + EquationOfStateType m_equationOfState; + arrayView1d< real64 const > m_criticalTemperature; + HeatCapacityCoefficients const * const m_heatCapacityCoefficients{nullptr}; +}; + +template< integer USD1, integer USD2 > +GEOS_HOST_DEVICE +void CompositionalEnthalpyUpdate::compute( + ComponentProperties::KernelWrapper const & componentProperties, + real64 const & pressure, + real64 const & temperature, + arraySlice1d< real64 const, USD1 > const & phaseComposition, + real64 & enthalpy, + arraySlice1d< real64, USD2 > const & dEnthalpy, + bool useMass ) const +{ + constexpr integer maxNumDof = MultiFluidConstants::MAX_NUM_COMPONENTS + 2; + integer const numComps = componentProperties.m_componentMolarWeight.size(); + integer const numDofs = numComps + 2; + + auto const & componentMolarWeight = componentProperties.m_componentMolarWeight; + auto const & criticalTemperature = componentProperties.m_componentCriticalTemperature; + + // 1. Calculate the polynomial (ideal) enthalpy + real64 const dT = temperature - m_refTemperature; + real64 & hIdealEnthalpy = enthalpy; + auto const & dhIdealEnthalpy = dEnthalpy; + hIdealEnthalpy = 0.0; + dhIdealEnthalpy[Deriv::dT] = 0.0; + dhIdealEnthalpy[Deriv::dP] = 0.0; + for( integer ic = 0; ic < numComps; ++ic ) + { + real64 enthalpyI = 0.0; + real64 heatCapacityI = 0.0; + evaluatePolynomial( dT, m_coefficients[ic], enthalpyI, heatCapacityI ); + // If temperature is greater that critical temperature, use the gas enthalpy if the gas phase exists + if( criticalTemperature[ic] < temperature && 0 <= m_vapourIndex ) + { + enthalpyI += m_referenceEnthalpy( m_vapourIndex, ic ); + } + else + { + enthalpyI += m_referenceEnthalpy( m_phaseIndex, ic ); + } + hIdealEnthalpy += phaseComposition[ic] * enthalpyI; + dhIdealEnthalpy[Deriv::dT] += phaseComposition[ic] * heatCapacityI; + dhIdealEnthalpy[Deriv::dC+ic] = enthalpyI; + } + + // 2. Calculate the departure enthalpy from the EoS + real64 hEosEnthalpy = 0.0; + stackArray1d< real64, maxNumDof > dhEosEnthalpy( numDofs ); + LvArray::forValuesInSlice( dhEosEnthalpy.toSlice(), setZero ); + selectEquationOfState( m_equationOfState, [&]( auto eos ) + { + using CubicModel = CubicEOSPhaseModel< decltype(eos) >; + calculateEquationOfStateEnthalpy< CubicModel >( numComps, + componentProperties, + pressure, + temperature, + phaseComposition, + hEosEnthalpy, + dhEosEnthalpy.toSlice(), + m_phaseType ); + } ); + + // 3. Combine the 2 values for the actual fluid enthalpy + enthalpy += hEosEnthalpy; + for( integer idof = 0; idof < numDofs; idof++ ) + { + dEnthalpy[idof] += dhEosEnthalpy[idof]; + } + + // In the mass formulation, the enthalpy will be required to be (annoyingly) in J/kg. The + // composition however is preconverted to molar composition before we get to this point. + // The derivatives will also be output in molar composition and conversion is done at a + // higher level within the fluid model itself. The only change here is to the "unit" of + // enthalpy. + // Convert from J/mol -> J/kg + if( useMass ) + { + real64 phaseMolarWeight = 0.0; + for( integer ic = 0; ic < numComps; ++ic ) + { + phaseMolarWeight += phaseComposition[ic] * componentMolarWeight[ic]; + } + real64 const invPhaseMolarWeight = 1.0 / phaseMolarWeight; + enthalpy *= invPhaseMolarWeight; + dEnthalpy[Deriv::dP] *= invPhaseMolarWeight; + dEnthalpy[Deriv::dT] *= invPhaseMolarWeight; + for( integer ic = 0; ic < numComps; ++ic ) + { + dEnthalpy[Deriv::dC+ic] = (dEnthalpy[Deriv::dC+ic] - componentMolarWeight[ic] * enthalpy)*invPhaseMolarWeight; + } + } +} + +template< typename LAMBDA > +GEOS_HOST_DEVICE +void CompositionalEnthalpyUpdate::selectEquationOfState( EquationOfStateType equationOfState, LAMBDA && lambda ) +{ + if( equationOfState == EquationOfStateType::SoaveRedlichKwong ) + { + std::forward< LAMBDA >( lambda )( SoaveRedlichKwongEOS{} ); + } + else if( equationOfState == EquationOfStateType::SoreideWhitson ) + { + // Soreide-Whitson uses underying Peng-Robinson formulation + std::forward< LAMBDA >( lambda )( PengRobinsonEOS{} ); + } + else + { + std::forward< LAMBDA >( lambda )( PengRobinsonEOS{} ); + } +} + +template< typename EOS, integer USD1, integer USD2 > +GEOS_HOST_DEVICE +void CompositionalEnthalpyUpdate::calculateEquationOfStateEnthalpy( + integer const numComps, + ComponentProperties::KernelWrapper const & componentProperties, + real64 const & pressure, + real64 const & temperature, + arraySlice1d< real64 const, USD1 > const & phaseComposition, + real64 & enthalpy, + arraySlice1d< real64, USD2 > const & dEnthalpy, + PhaseType const & phaseType ) const +{ + integer sizes[2] = {0, 0}; + arraySlice2d< real64 const > derivs( nullptr, sizes, sizes ); + typename EOS::template StackVariables< true > stack( numComps, + componentProperties.m_componentBinaryCoeff.toSliceConst(), + derivs.toSliceConst() ); + + typename EOS::SelectedRoot const root = phaseType == PhaseType::LIQUID ? + EOS::SelectedRoot::MINIMUM : EOS::SelectedRoot::MAXIMUM; + + EOS::initialiseStack( numComps, pressure, temperature, componentProperties, stack ); + EOS::computeMixtureCoefficients( numComps, phaseComposition, stack ); + EOS::template computeEnthalpy< USD1, true >( numComps, + pressure, + temperature, + phaseComposition, + componentProperties, + stack, + enthalpy, + dEnthalpy, + root ); +} + +} // end namespace compositional + +} // end namespace constitutive + +} // end namespace geos + +#endif //GEOS_CONSTITUTIVE_FLUID_MULTIFLUID_COMPOSITIONAL_MODELS_COMPOSITIONALENTHALPY_HPP_ diff --git a/src/coreComponents/constitutive/fluid/multifluid/compositional/parameters/HeatCapacityCoefficients.cpp b/src/coreComponents/constitutive/fluid/multifluid/compositional/parameters/HeatCapacityCoefficients.cpp new file mode 100644 index 00000000000..72ba65a001a --- /dev/null +++ b/src/coreComponents/constitutive/fluid/multifluid/compositional/parameters/HeatCapacityCoefficients.cpp @@ -0,0 +1,342 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2024 TotalEnergies + * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2023-2024 Chevron + * Copyright (c) 2019- GEOS/GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +/** + * @file HeatCapacityCoefficients.cpp + */ + +#include "HeatCapacityCoefficients.hpp" +#include "ComponentProperties.hpp" +#include "EquationOfState.hpp" +#include "constitutive/fluid/multifluid/MultiFluidBase.hpp" +#include "constitutive/fluid/multifluid/MultiFluidConstants.hpp" +#include "dataRepository/InputFlags.hpp" +#include "common/PhysicsConstants.hpp" + +namespace geos +{ + +namespace constitutive +{ + +namespace compositional +{ + +HeatCapacityCoefficients::HeatCapacityCoefficients( std::unique_ptr< ModelParameters > parameters ): + ModelParameters( std::move( parameters ) ) +{} + +std::unique_ptr< ModelParameters > +HeatCapacityCoefficients::create( std::unique_ptr< ModelParameters > parameters ) +{ + if( parameters && parameters->get< HeatCapacityCoefficients >() != nullptr ) + { + return parameters; + } + parameters = EquationOfState::create( std::move( parameters ) ); + return std::make_unique< HeatCapacityCoefficients >( std::move( parameters ) ); +} + +void HeatCapacityCoefficients::registerParametersImpl( MultiFluidBase * fluid ) +{ + fluid->registerWrapper( viewKeyStruct::enthalpyReferenceTemperatureString(), &m_referenceTemperature ). + setInputFlag( dataRepository::InputFlags::OPTIONAL ). + setDefaultValue( m_referenceTemperature ). + setDescription( "The reference temperature for enthalpy calculation" ); + + fluid->registerWrapper( viewKeyStruct::referenceEnthalpyString(), &m_referenceEnthalpy ). + setInputFlag( dataRepository::InputFlags::OPTIONAL ). + setDescription( "The enthalpy of each component in each of the phases at the reference temperature." ); + + fluid->registerWrapper( viewKeyStruct::componentHeatCapacityCoefficientsString(), &m_coefficients ). + setInputFlag( dataRepository::InputFlags::REQUIRED ). + setDescription( "The polynomial coefficients for the specific heat capacity of each component." ); +} + +void HeatCapacityCoefficients::postInputInitializationImpl( MultiFluidBase const * fluid, + ComponentProperties const & componentProperties ) +{ + integer const numPhases = fluid->numFluidPhases(); + integer const numComps = fluid->numFluidComponents(); + + // If the reference enthalpies are given, then there must be as many as there are components + if( m_referenceEnthalpy.empty() ) + { + m_referenceEnthalpy.resize( numPhases, numComps ); + m_referenceEnthalpy.zero(); + } + else + { + integer const dim0 = m_referenceEnthalpy.size( 0 ); + integer const dim1 = m_referenceEnthalpy.size( 1 ); + + GEOS_THROW_IF_NE_MSG( dim0, numPhases, + GEOS_FMT( "{}: '{}' the first dimension must be equal to the number of phases {}", + fluid->getFullName(), + viewKeyStruct::referenceEnthalpyString(), + numPhases ), + InputError ); + GEOS_THROW_IF_NE_MSG( dim1, numComps, + GEOS_FMT( "{}: '{}' the second dimension must be equal to the number of components {}", + fluid->getFullName(), + viewKeyStruct::referenceEnthalpyString(), + numComps ), + InputError ); + } + + // Reference enthalpies must be ordered + m_phaseTypes.resize( numPhases ); + string_array const & phaseNames = fluid->phaseNames(); + string_array const & componentNames = fluid->componentNames(); + + integer ipL = -1, ipV = -1, ipA = -1; + for( integer ip = 0; ip < numPhases; ip++ ) + { + PhaseType const phaseType = getPhaseTypeFromName( phaseNames[ip] ); + m_phaseTypes[ip] = phaseType; + switch( phaseType ) + { + case PhaseType::LIQUID: + ipL = ip; + break; + case PhaseType::VAPOUR: + ipV = ip; + break; + case PhaseType::AQUEOUS: + ipA = ip; + break; + } + } + // Enthalpy must be Oil < Water < Gas + stdVector< std::pair< integer, integer > > pairs; + if( 0 <= ipL && 0 <= ipV ) + { + pairs.emplace_back( ipL, ipV ); + } + if( 0 <= ipL && 0 <= ipA ) + { + pairs.emplace_back( ipL, ipA ); + } + if( 0 <= ipA && 0 <= ipV ) + { + pairs.emplace_back( ipA, ipV ); + } + for( auto const & [ip1, ip2] : pairs ) + { + for( integer ic = 0; ic < numComps; ++ic ) + { + real64 const phase1Enthalpy = m_referenceEnthalpy( ip1, ic ); + real64 const phase2Enthalpy = m_referenceEnthalpy( ip2, ic ); + GEOS_THROW_IF_GT_MSG( phase1Enthalpy, phase2Enthalpy, + GEOS_FMT( "{}: '{}' for component {}, the {} reference enthalpy {} must not " + "be greater than the {} reference enthalpy {}.", + fluid->getFullName(), + viewKeyStruct::referenceEnthalpyString(), + componentNames[ic], phaseNames[ip1], phase1Enthalpy, phaseNames[ip2], phase2Enthalpy ), + InputError ); + } + } + + { + integer const dim0 = m_coefficients.size( 0 ); + integer const dim1 = m_coefficients.size( 1 ); + + // First dimension must be equal to number of components + GEOS_THROW_IF_NE_MSG( dim0, numComps, + GEOS_FMT( "{}: '{}' the first dimension must be equal to the number of components {}", + fluid->getFullName(), + viewKeyStruct::componentHeatCapacityCoefficientsString(), + numComps ), + InputError ); + // Second dimension must be equal to 5 + GEOS_THROW_IF_NE_MSG( dim1, 5, + GEOS_FMT( "{}: '{}' the second dimension must be equal 5", + fluid->getFullName(), + viewKeyStruct::componentHeatCapacityCoefficientsString() ), + InputError ); + } + + // Reference temperature must not be negative + GEOS_THROW_IF_LE_MSG( m_referenceTemperature, 0.0, + GEOS_FMT( "{}: '{}' the reference temperature {} must not be zero or negative.", + fluid->getFullName(), + viewKeyStruct::enthalpyReferenceTemperatureString(), + m_referenceTemperature ), + InputError ); + + // Determine temperature range + constexpr real64 zeroC = constants::zeroDegreesCelsiusInKelvin; + real64 minTemperature = LvArray::math::min( zeroC, m_referenceTemperature ); + real64 maxTemperature = LvArray::math::max( zeroC, m_referenceTemperature ); + auto const criticalTemperature = componentProperties.getComponentCriticalTemperature(); + for( integer ic = 0; ic < numComps; ++ic ) + { + minTemperature = LvArray::math::min( minTemperature, criticalTemperature[ic] ); + maxTemperature = LvArray::math::max( maxTemperature, criticalTemperature[ic] ); + } + + // Extend interval by 10% in each direction + real64 const dt = LvArray::math::max( maxTemperature - minTemperature, 100.0 ); + minTemperature -= 0.1 * dt; + maxTemperature += 0.1 * dt; + + // Transform to reference temperature space + minTemperature -= m_referenceTemperature; + maxTemperature -= m_referenceTemperature; + + real64 negT = 0.0; + real64 negHT = 0.0; + + for( integer ic = 0; ic < numComps; ++ic ) + { + bool isPositive = isPolynomialPositive( m_coefficients[ic], minTemperature, maxTemperature, negT, negHT ); + GEOS_THROW_IF( !isPositive, + GEOS_FMT( "{}: '{}' coefficients for component {} ({}) give a zero or negative heat " + "capacity of {} at temperature difference {} corresponding to a temperature {}.", + fluid->getFullName(), + viewKeyStruct::componentHeatCapacityCoefficientsString(), + componentNames[ic], m_coefficients[ic], negHT, negT, negT + m_referenceTemperature ), + InputError ); + } +} + +/* + * @brief Evaluates whether a fourth-order polynomial remains non-negative across a specified interval. + * + * @details + * This method verifies that the polynomial h(T) does not fall below a safety threshold + * (epsilon) for all T in the range [T0, T1]. + * The algorithm combines an integration-like stepping mechanism with an adaptively computed + * local Lipschitz constant to determine safe step sizes. + * At any current temperature, the maximum possible rate of decrease (the local Lipschitz constant) + * is calculated by finding the absolute maximum of the polynomial's derivative on the remaining + * interval [currentT, T1]. Because the derivative is a cubic function, its maximums can only + * exist at the boundaries of the interval or at the roots of the second derivative (which is a + * quadratic equation). + * Once this maximum absolute derivative is found, a safe forward step is calculated as + * deltaT = (h(currentT) - epsilon) / localLipschitz. This step mathematically guarantees + * that the function cannot drop below epsilon within the step. + * If the function evaluates to less than epsilon at any point, the iteration halts and + * updates the reference parameters with the failing coordinates. + * + * @param a Array of 5 coefficients representing the polynomial, ordered from T^0 to T^4. + * @param T0 Start of the evaluation range for temperature T. + * @param T1 End of the evaluation range for temperature T. + * @param[out] T Reference output parameter that will store the failing temperature if a drop occurs. + * @param[out] hT Reference output parameter that will store the failing polynomial value. + * @return true if the polynomial remains strictly >= epsilon, false otherwise. + */ +bool HeatCapacityCoefficients::isPolynomialPositive( arraySlice1d< real64 const > const a, + real64 const T0, real64 const T1, + real64 & T, real64 & hT ) +{ + constexpr real64 epsilon = MultiFluidConstants::epsilon; + constexpr real64 minStep = 1.0e-13; + + real64 currentT = T0; + real64 currentH; + + // Precalculate the coefficients of the second derivative h''(T) + real64 const cA = 12.0 * a[4]; + real64 const cB = 6.0 * a[3]; + real64 const cC = 2.0 * a[2]; + + // Precalculate the roots of the second derivative + real64 roots[2]{}; + integer numRoots = 0; + + if( LvArray::math::abs( cA ) < epsilon ) + { + if( !(LvArray::math::abs( cB ) < epsilon)) + { + roots[numRoots++] = -cC / cB; + } + } + else + { + real64 const discriminant = cB * cB - 4.0 * cA * cC; + if( discriminant >= 0.0 ) + { + real64 const sqrtDisc = LvArray::math::sqrt( discriminant ); + roots[numRoots++] = (-cB - sqrtDisc) / (2.0 * cA); + roots[numRoots++] = (-cB + sqrtDisc) / (2.0 * cA); + } + } + + // Lambda function to compute the absolute value of the first derivative + auto const getAbsDeriv = [&]( real64 const t ) -> real64 { + real64 const deriv = a[1] + t * (2.0 * a[2] + t * (3.0 * a[3] + t * 4.0 * a[4])); + return LvArray::math::abs( deriv ); + }; + + while( true ) + { + // Evaluate the polynomial using Horner's method for stability and performance + currentH = a[0] + currentT * (a[1] + currentT * (a[2] + currentT * (a[3] + currentT * a[4]))); + + // If the function falls below the designated epsilon threshold, fail immediately + if( currentH < epsilon ) + { + T = currentT; + hT = currentH; + return false; + } + + // If we successfully reached or surpassed the end of the interval, return success + if( currentT >= T1 ) + { + T = currentT; + hT = currentH; + return true; + } + + // Compute the local Lipschitz constant (max absolute derivative) over [currentT, T1] + real64 const valCur = getAbsDeriv( currentT ); + real64 const valEnd = getAbsDeriv( T1 ); + + real64 maxAbsDeriv = LvArray::math::max( valCur, valEnd ); + + for( integer i = 0; i < numRoots; ++i ) + { + if( roots[i] > currentT && roots[i] < T1 ) + { + maxAbsDeriv = LvArray::math::max( maxAbsDeriv, getAbsDeriv( roots[i] )); + } + } + + // If the derivative is zero across the entire remainder, the function is flat and safe + if( LvArray::math::abs( maxAbsDeriv ) < epsilon ) + { + T = T1; + hT = currentH; + return true; + } + + // Calculate the maximally safe local step size + real64 const calculatedStep = (currentH - epsilon) / maxAbsDeriv; + + // Ensure forward progress even when the derivative is extremely large near epsilon + real64 const actualStep = LvArray::math::max( calculatedStep, minStep ); + + currentT += actualStep; + + // Ensure we do not overshoot the boundary on the final iteration + currentT = LvArray::math::min( currentT, T1 ); + } +} + +} // end namespace compositional +} // end namespace constitutive +} // end namespace geos diff --git a/src/coreComponents/constitutive/fluid/multifluid/compositional/parameters/HeatCapacityCoefficients.hpp b/src/coreComponents/constitutive/fluid/multifluid/compositional/parameters/HeatCapacityCoefficients.hpp new file mode 100644 index 00000000000..9d4375c700b --- /dev/null +++ b/src/coreComponents/constitutive/fluid/multifluid/compositional/parameters/HeatCapacityCoefficients.hpp @@ -0,0 +1,71 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2024 TotalEnergies + * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2023-2024 Chevron + * Copyright (c) 2019- GEOS/GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +/** + * @file HeatCapacityCoefficients.hpp + */ + +#ifndef GEOS_CONSTITUTIVE_FLUID_MULTIFLUID_COMPOSITIONAL_PARAMETERS_HEATCAPACITYCOEFFICIENTS_HPP_ +#define GEOS_CONSTITUTIVE_FLUID_MULTIFLUID_COMPOSITIONAL_PARAMETERS_HEATCAPACITYCOEFFICIENTS_HPP_ + +#include "ModelParameters.hpp" +#include "PhaseType.hpp" +#include "common/DataTypes.hpp" + +namespace geos +{ + +namespace constitutive +{ + +namespace compositional +{ + +class HeatCapacityCoefficients : public ModelParameters +{ +public: + HeatCapacityCoefficients( std::unique_ptr< ModelParameters > parameters ); + ~HeatCapacityCoefficients() override = default; + + static std::unique_ptr< ModelParameters > create( std::unique_ptr< ModelParameters > parameters ); + + struct viewKeyStruct + { + static constexpr char const * enthalpyReferenceTemperatureString() { return "enthalpyReferenceTemperature"; } + static constexpr char const * referenceEnthalpyString() { return "referenceEnthalpy"; } + static constexpr char const * componentHeatCapacityCoefficientsString() { return "componentHeatCapacityCoefficients"; } + }; + + real64 m_referenceTemperature{298.15}; + array2d< real64 > m_referenceEnthalpy; + array2d< real64 > m_coefficients; + array1d< PhaseType > m_phaseTypes; + +protected: + void registerParametersImpl( MultiFluidBase * fluid ) override; + + void postInputInitializationImpl( MultiFluidBase const * fluid, ComponentProperties const & componentProperties ) override; + +private: + static bool isPolynomialPositive( arraySlice1d< real64 const > const a, real64 const T0, real64 const T1, real64 & T, real64 & hT ); +}; + +} // end namespace compositional + +} // end namespace constitutive + +} // end namespace geos + +#endif //GEOS_CONSTITUTIVE_FLUID_MULTIFLUID_COMPOSITIONAL_PARAMETERS_HEATCAPACITYCOEFFICIENTS_HPP_ diff --git a/src/coreComponents/constitutive/unitTests/CMakeLists.txt b/src/coreComponents/constitutive/unitTests/CMakeLists.txt index 9dfbfc2c8a3..03be5f7d139 100644 --- a/src/coreComponents/constitutive/unitTests/CMakeLists.txt +++ b/src/coreComponents/constitutive/unitTests/CMakeLists.txt @@ -4,6 +4,7 @@ set( gtest_geosx_tests testCO2SpycherPruessModels.cpp testCO2BrinePVTModels.cpp testCompositionalDensity.cpp + testCompositionalEnthalpy.cpp testCompositionalPhillipsBrineDensity.cpp testCompositionalPhillipsBrineViscosity.cpp testCompositionalProperties.cpp @@ -12,6 +13,7 @@ set( gtest_geosx_tests testDamageUtilities.cpp testDruckerPrager.cpp testElasticIsotropic.cpp + testHeatCapacityCoefficients.cpp testInvariantImmiscibleFluid.cpp testKValueFlashModel.cpp testKValueInitialization.cpp diff --git a/src/coreComponents/constitutive/unitTests/testCompositionalEnthalpy.cpp b/src/coreComponents/constitutive/unitTests/testCompositionalEnthalpy.cpp new file mode 100644 index 00000000000..bf52664e2ac --- /dev/null +++ b/src/coreComponents/constitutive/unitTests/testCompositionalEnthalpy.cpp @@ -0,0 +1,422 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2024 TotalEnergies + * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2023-2024 Chevron + * Copyright (c) 2019- GEOS/GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +// Source includes +#include "codingUtilities/UnitTestUtilities.hpp" +#include "constitutive/fluid/multifluid/compositional/models/CompositionalEnthalpy.hpp" +#include "constitutive/fluid/multifluid/compositional/parameters/EquationOfState.hpp" +#include "constitutive/fluid/multifluid/compositional/parameters/HeatCapacityCoefficients.hpp" + +#include "TestFluid.hpp" +#include "TestFluidUtilities.hpp" + +using namespace geos::constitutive::compositional; + +namespace geos +{ +namespace testing +{ + +template< int NC > +using EnthalpyData = std::tuple< + real64 const, // pressure + real64 const, // temperature + Feed< NC > const, // phase composition + integer const, // phase + real64 const, // expected enthalpy + real64 const // expected Joule-Thomson coefficient (K/Pa) + >; + +template< int NC > +struct FluidData {}; + +template<> +struct FluidData< 4 > +{ + static constexpr integer numPhases = 2; + static constexpr integer numComps = 4; + static constexpr integer numCoeffs = 5; + + static std::unique_ptr< TestFluid< numComps > > createFluid() + { + return TestFluid< numComps >::create( {Fluid::CO2, Fluid::H2, Fluid::CH4, Fluid::C8H18} ); + } + + static void populateCoefficients( HeatCapacityCoefficients * coefficients ) + { + coefficients->m_referenceEnthalpy.resize( numPhases, numComps ); + coefficients->m_referenceEnthalpy.zero(); + + coefficients->m_coefficients.resize( numComps, numCoeffs ); + std::array< real64, numComps * numCoeffs > coefficientsData { + 3.45314594e+01, 7.07764122e-02, 1.69716907e-04, -1.40256545e-06, 1.83892934e-09, + 2.87501790e+01, 1.01246159e-02, -5.85965230e-05, 1.30968549e-07, -9.36636594e-11, + 3.55287791e+01, 4.22215766e-02, 1.34958295e-04, -4.13365988e-07, 4.05670019e-10, + 2.01824087e+02, 3.07940012e-01, 6.82371480e-04, -1.07019612e-06, -3.48458322e-11 + }; + for( int ic = 0; ic < numComps; ++ic ) + { + for( int j = 0; j < 5; ++j ) + { + coefficients->m_coefficients( ic, j ) = coefficientsData[ic*numCoeffs + j]; + } + } + coefficients->m_phaseTypes.resize( numPhases ); + coefficients->m_phaseTypes[static_cast< integer >(PhaseType::LIQUID)] = PhaseType::LIQUID; + coefficients->m_phaseTypes[static_cast< integer >(PhaseType::VAPOUR)] = PhaseType::VAPOUR; + } +}; + +template<> +struct FluidData< 6 > +{ + static constexpr integer numPhases = 2; + static constexpr integer numComps = 6; + static constexpr integer numCoeffs = 5; + + static std::unique_ptr< TestFluid< numComps > > createFluid() + { + return TestFluid< numComps >::create( {Fluid::CO2, Fluid::N2, Fluid::H2O, Fluid::CH4, Fluid::H2S, Fluid::C5H12} ); + } + + static void populateCoefficients( HeatCapacityCoefficients * coefficients ) + { + coefficients->m_referenceEnthalpy.resize( numPhases, numComps ); + coefficients->m_referenceEnthalpy.zero(); + + coefficients->m_coefficients.resize( numComps, numCoeffs ); + std::array< real64, numComps * numCoeffs > coefficientsData { + 3.40933377e+01, 7.52742132e-02, 2.00271751e-04, -1.61900936e-06, 2.12199565e-09, + 2.89614125e+01, 2.03166820e-03, 2.47038198e-05, -9.54674193e-08, 1.15955044e-10, + 2.69032982e+01, 5.51993739e-02, 6.31405927e-04, -3.66166336e-06, 4.56291585e-09, + 3.54137115e+01, 4.34392612e-02, 1.42852598e-04, -4.70393972e-07, 4.80614694e-10, + 3.18873193e+01, 3.25765666e-02, 2.20667710e-04, -1.29890488e-06, 1.63371768e-09, + 1.21004526e+02, 2.47982404e-01, 7.27115529e-04, -2.64999262e-06, 2.52476069e-09 + }; + for( int ic = 0; ic < numComps; ++ic ) + { + for( int j = 0; j < 5; ++j ) + { + coefficients->m_coefficients( ic, j ) = coefficientsData[ic*numCoeffs + j]; + } + } + coefficients->m_phaseTypes.resize( numPhases ); + coefficients->m_phaseTypes[static_cast< integer >(PhaseType::LIQUID)] = PhaseType::LIQUID; + coefficients->m_phaseTypes[static_cast< integer >(PhaseType::VAPOUR)] = PhaseType::VAPOUR; + } +}; + +template< EquationOfStateType EOS_TYPE, int NUM_COMP > +class CompositionalEnthalpyTestFixture : public ::testing::TestWithParam< EnthalpyData< NUM_COMP > > +{ + static constexpr real64 relTol = 1.0e-5; + static constexpr real64 absTol = 1.0e-7; + static constexpr int numComps = NUM_COMP; + static constexpr int numDofs = NUM_COMP + 2; + using Deriv = geos::constitutive::multifluid::DerivativeOffset; + +public: + using DataType = EnthalpyData< numComps >; + +public: + CompositionalEnthalpyTestFixture() + : m_fluid( FluidData< numComps >::createFluid() ) + { + ComponentProperties const & componentProperties = this->m_fluid->getComponentProperties(); + m_parameters = CompositionalEnthalpy::createParameters( std::make_unique< ModelParameters >() ); + + auto equationOfState = const_cast< EquationOfState * >(m_parameters->get< EquationOfState >()); + string const eosName = EnumStrings< EquationOfStateType >::toString( EOS_TYPE ); + equationOfState->m_equationsOfStateNames.emplace_back( eosName ); + equationOfState->m_equationsOfStateNames.emplace_back( eosName ); + + auto heatCapacityCoefficients = const_cast< HeatCapacityCoefficients * >(m_parameters->get< HeatCapacityCoefficients >()); + FluidData< numComps >::populateCoefficients( heatCapacityCoefficients ); + + string const name = GEOS_FMT( "PhaseEnthalpy{}{}", eosName, numComps ); + m_enthalpy = std::make_unique< CompositionalEnthalpy >( name, componentProperties, 0, *m_parameters ); + } + + ~CompositionalEnthalpyTestFixture() = default; + + void testEnthalpyValues( DataType const & data, bool useMass = false ) + { + real64 const pressure = std::get< 0 >( data ); + real64 const temperature = std::get< 1 >( data ); + stackArray1d< real64, numComps > phaseComposition; + TestFluid< numComps >::createArray( phaseComposition, std::get< 2 >( data )); + setPhase( std::get< 3 >( data ) ); + real64 const & expectedEnthalpy = std::get< 4 >( data ); + real64 const & expectedJouleThomson = std::get< 5 >( data ); + + auto componentProperties = m_fluid->createKernelWrapper(); + auto kernelWrapper = m_enthalpy->createKernelWrapper(); + + real64 enthalpy = 0.0; + stackArray1d< real64, numDofs > enthalpyDerivs( numDofs ); + + kernelWrapper.compute( componentProperties, + pressure, + temperature, + phaseComposition.toSliceConst(), + enthalpy, + enthalpyDerivs.toSlice(), + useMass ); + + // Approximate JT coefficient + real64 const jouleThomson = -enthalpyDerivs[Deriv::dP] / enthalpyDerivs[Deriv::dT]; + + checkRelativeError( enthalpy, expectedEnthalpy, relTol, absTol ); + checkRelativeError( jouleThomson, expectedJouleThomson, relTol, absTol ); + } + + void testEnthalpyDerivatives( DataType const & data, bool useMass ) + { + real64 const pressure = std::get< 0 >( data ); + real64 const temperature = std::get< 1 >( data ); + stackArray1d< real64, numComps > phaseComposition; + TestFluid< numComps >::createArray( phaseComposition, std::get< 2 >( data )); + setPhase( std::get< 3 >( data ) ); + + auto componentProperties = m_fluid->createKernelWrapper(); + auto kernelWrapper = m_enthalpy->createKernelWrapper(); + + real64 enthalpy = 0.0; + stackArray2d< real64, 2*numDofs > derivSpace( 2, numDofs ); + arraySlice1d< real64 > enthalpyDerivs = derivSpace[0]; + arraySlice1d< real64 > tempDerivs = derivSpace[1]; + + kernelWrapper.compute( componentProperties, + pressure, + temperature, + phaseComposition.toSliceConst(), + enthalpy, + enthalpyDerivs, + useMass ); + // Compare against numerical derivatives + // -- Pressure derivative + real64 const dp = 1.0e-4 * pressure; + internal::testNumericalDerivative( + pressure, dp, enthalpyDerivs[Deriv::dP], + [&]( real64 const p ){ + real64 displacedEnthalpy = 0.0; + kernelWrapper.compute( componentProperties, + p, + temperature, + phaseComposition.toSliceConst(), + displacedEnthalpy, + tempDerivs, + useMass ); + return displacedEnthalpy; + } ); + + // -- Temperature derivative + real64 const dT = 1.0e-5 * temperature; + internal::testNumericalDerivative( + temperature, dT, enthalpyDerivs[Deriv::dT], + [&]( real64 const t ){ + real64 displacedEnthalpy = 0.0; + kernelWrapper.compute( componentProperties, + pressure, + t, + phaseComposition.toSliceConst(), + displacedEnthalpy, + tempDerivs, + useMass ); + return displacedEnthalpy; + } ); + + // -- Composition derivatives + real64 constexpr dz = 1.0e-6; + for( integer ic = 0; ic < numComps; ++ic ) + { + internal::testNumericalDerivative( + 0.0, dz, enthalpyDerivs[Deriv::dC+ic], + [&]( real64 const z ){ + real64 const z_old = phaseComposition[ic]; + phaseComposition[ic] += z; + real64 displacedEnthalpy = 0.0; + kernelWrapper.compute( componentProperties, + pressure, + temperature, + phaseComposition.toSliceConst(), + displacedEnthalpy, + tempDerivs, + useMass ); + phaseComposition[ic] = z_old; + return displacedEnthalpy; + } ); + } + } + +protected: + void setPhase( integer const phaseIndex ) + { + ComponentProperties const & componentProperties = m_fluid->getComponentProperties(); + string const name = m_enthalpy->functionName(); + m_enthalpy = std::make_unique< CompositionalEnthalpy >( name, componentProperties, phaseIndex, *m_parameters ); + } + +protected: + std::unique_ptr< CompositionalEnthalpy > m_enthalpy{}; + std::unique_ptr< TestFluid< numComps > > m_fluid{}; + std::unique_ptr< ModelParameters > m_parameters{}; +}; + +using PengRobinson_4 = CompositionalEnthalpyTestFixture< EquationOfStateType::PengRobinson, 4 >; +using SoaveRedlichKwong_6 = CompositionalEnthalpyTestFixture< EquationOfStateType::SoaveRedlichKwong, 6 >; +using SoreideWhitson_6 = CompositionalEnthalpyTestFixture< EquationOfStateType::SoreideWhitson, 6 >; + +TEST_P( PengRobinson_4, testEnthalpyValues ) +{ + testEnthalpyValues( GetParam() ); +} +TEST_P( SoaveRedlichKwong_6, testEnthalpyValues ) +{ + testEnthalpyValues( GetParam() ); +} +TEST_P( SoreideWhitson_6, testEnthalpyValues ) +{ + testEnthalpyValues( GetParam(), true ); +} + +TEST_P( PengRobinson_4, testEnthalpyDerivatives ) +{ + testEnthalpyDerivatives( GetParam(), false ); + testEnthalpyDerivatives( GetParam(), true ); +} +TEST_P( SoaveRedlichKwong_6, testEnthalpyDerivatives ) +{ + testEnthalpyDerivatives( GetParam(), false ); + testEnthalpyDerivatives( GetParam(), true ); +} +TEST_P( SoreideWhitson_6, testEnthalpyDerivatives ) +{ + testEnthalpyDerivatives( GetParam(), false ); + testEnthalpyDerivatives( GetParam(), true ); +} + +/* UNCRUSTIFY-OFF */ + +// Test data + +constexpr integer LIQ = static_cast(PhaseType::LIQUID); +constexpr integer VAP = static_cast(PhaseType::VAPOUR); + +INSTANTIATE_TEST_SUITE_P( + CompositionalEnthalpyTest, PengRobinson_4, + ::testing::ValuesIn< typename PengRobinson_4::DataType >( { + {1.01325e+05, 283.15, {1.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00}, VAP, -5.554727e+02, 1.329174e-05}, + {1.00000e+06, 283.15, {1.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00}, VAP, -9.777036e+02, 1.354551e-05}, + {1.00000e+07, 283.15, {1.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00}, LIQ, -1.258949e+04, 6.097076e-07}, + {1.01325e+05, 298.15, {1.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00}, VAP, -4.172063e+01, 1.189416e-05}, + {1.00000e+06, 298.15, {1.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00}, VAP, -4.271945e+02, 1.203740e-05}, + {1.00000e+07, 298.15, {1.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00}, LIQ, -1.075604e+04, 1.246819e-06}, + {1.01325e+05, 550.00, {1.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00}, VAP, 1.079339e+04, 2.795143e-06}, + {1.00000e+06, 550.00, {1.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00}, VAP, 1.067299e+04, 2.740228e-06}, + {1.00000e+07, 550.00, {1.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00}, VAP, 9.544159e+03, 2.175503e-06}, + {1.01325e+05, 283.15, {0.000000e+00, 1.000000e+00, 0.000000e+00, 0.000000e+00}, VAP, -4.305421e+02, 1.699485e-07}, + {1.00000e+06, 283.15, {0.000000e+00, 1.000000e+00, 0.000000e+00, 0.000000e+00}, VAP, -4.346280e+02, 1.478298e-07}, + {1.00000e+07, 283.15, {0.000000e+00, 1.000000e+00, 0.000000e+00, 0.000000e+00}, LIQ, -4.479027e+02, -3.415949e-08}, + {1.01325e+05, 298.15, {0.000000e+00, 1.000000e+00, 0.000000e+00, 0.000000e+00}, VAP, -3.809923e-01, 1.296216e-07}, + {1.00000e+06, 298.15, {0.000000e+00, 1.000000e+00, 0.000000e+00, 0.000000e+00}, VAP, -3.474631e+00, 1.096517e-07}, + {1.00000e+07, 298.15, {0.000000e+00, 1.000000e+00, 0.000000e+00, 0.000000e+00}, LIQ, -9.208397e+00, -5.550701e-08}, + {1.01325e+05, 550.00, {0.000000e+00, 1.000000e+00, 0.000000e+00, 0.000000e+00}, LIQ, 7.363162e+03, -2.039014e-07}, + {1.00000e+06, 550.00, {0.000000e+00, 1.000000e+00, 0.000000e+00, 0.000000e+00}, LIQ, 7.368612e+03, -2.098955e-07}, + {1.00000e+07, 550.00, {0.000000e+00, 1.000000e+00, 0.000000e+00, 0.000000e+00}, LIQ, 7.431415e+03, -2.626841e-07}, + {1.01325e+05, 283.15, {0.000000e+00, 0.000000e+00, 1.000000e+00, 0.000000e+00}, VAP, -5.482831e+02, 5.621645e-06}, + {1.00000e+06, 283.15, {0.000000e+00, 0.000000e+00, 1.000000e+00, 0.000000e+00}, VAP, -7.262985e+02, 5.522977e-06}, + {1.00000e+07, 283.15, {0.000000e+00, 0.000000e+00, 1.000000e+00, 0.000000e+00}, VAP, -2.519542e+03, 3.642260e-06}, + {1.01325e+05, 298.15, {0.000000e+00, 0.000000e+00, 1.000000e+00, 0.000000e+00}, VAP, -1.838859e+01, 5.095639e-06}, + {1.00000e+06, 298.15, {0.000000e+00, 0.000000e+00, 1.000000e+00, 0.000000e+00}, VAP, -1.819651e+02, 4.994247e-06}, + {1.00000e+07, 298.15, {0.000000e+00, 0.000000e+00, 1.000000e+00, 0.000000e+00}, VAP, -1.781165e+03, 3.367126e-06}, + {1.01325e+05, 550.00, {0.000000e+00, 0.000000e+00, 1.000000e+00, 0.000000e+00}, VAP, 1.066624e+04, 1.144149e-06}, + {1.00000e+06, 550.00, {0.000000e+00, 0.000000e+00, 1.000000e+00, 0.000000e+00}, VAP, 1.061570e+04, 1.110647e-06}, + {1.00000e+07, 550.00, {0.000000e+00, 0.000000e+00, 1.000000e+00, 0.000000e+00}, VAP, 1.017885e+04, 8.033311e-07}, + {1.01325e+05, 283.15, {0.000000e+00, 0.000000e+00, 0.000000e+00, 1.000000e+00}, LIQ, -4.438247e+04, -5.187075e-07}, + {1.00000e+06, 283.15, {0.000000e+00, 0.000000e+00, 0.000000e+00, 1.000000e+00}, LIQ, -4.426586e+04, -5.208575e-07}, + {1.00000e+07, 283.15, {0.000000e+00, 0.000000e+00, 0.000000e+00, 1.000000e+00}, LIQ, -4.307920e+04, -5.386561e-07}, + {1.01325e+05, 298.15, {0.000000e+00, 0.000000e+00, 0.000000e+00, 1.000000e+00}, LIQ, -4.060294e+04, -4.948567e-07}, + {1.00000e+06, 298.15, {0.000000e+00, 0.000000e+00, 0.000000e+00, 1.000000e+00}, LIQ, -4.048960e+04, -4.975484e-07}, + {1.00000e+07, 298.15, {0.000000e+00, 0.000000e+00, 0.000000e+00, 1.000000e+00}, LIQ, -3.933020e+04, -5.194386e-07}, + {1.01325e+05, 550.00, {0.000000e+00, 0.000000e+00, 0.000000e+00, 1.000000e+00}, VAP, 6.288284e+04, 8.574199e-06}, + {1.00000e+06, 550.00, {0.000000e+00, 0.000000e+00, 0.000000e+00, 1.000000e+00}, VAP, 6.015485e+04, 1.126202e-05}, + {1.00000e+07, 550.00, {0.000000e+00, 0.000000e+00, 0.000000e+00, 1.000000e+00}, LIQ, 3.752013e+04, 2.732726e-07}, + {1.01325e+05, 283.15, {3.286000e-01, 3.319000e-01, 3.313000e-01, 8.300000e-03}, VAP, -5.286490e+02, 5.366763e-06}, + {1.01325e+05, 283.15, {1.040000e-02, 3.000000e-04, 2.100000e-03, 9.871000e-01}, LIQ, -4.395667e+04, -5.167514e-07}, + {1.00000e+06, 283.15, {3.121000e-01, 3.468000e-01, 3.400000e-01, 1.000000e-03}, VAP, -6.521012e+02, 4.842468e-06}, + {1.00000e+06, 283.15, {9.190000e-02, 3.600000e-03, 2.110000e-02, 8.835000e-01}, LIQ, -4.041707e+04, -5.013302e-07}, + {1.00000e+07, 283.15, {1.850000e-01, 4.662000e-01, 3.480000e-01, 7.000000e-04}, VAP, -1.481792e+03, 2.276903e-06}, + {1.00000e+07, 283.15, {3.097000e-01, 5.140000e-02, 1.599000e-01, 4.790000e-01}, LIQ, -2.534339e+04, -3.983758e-07}, + {1.01325e+05, 298.15, {3.253000e-01, 3.276000e-01, 3.272000e-01, 1.990000e-02}, VAP, -1.910796e+01, 5.184416e-06}, + {1.01325e+05, 298.15, {8.000000e-03, 4.000000e-04, 1.900000e-03, 9.898000e-01}, LIQ, -4.029301e+04, -4.930559e-07}, + {1.00000e+06, 298.15, {3.171000e-01, 3.431000e-01, 3.374000e-01, 2.500000e-03}, VAP, -1.524379e+02, 4.467927e-06}, + {1.00000e+06, 298.15, {7.240000e-02, 3.800000e-03, 1.880000e-02, 9.050000e-01}, LIQ, -3.758874e+04, -4.793543e-07}, + {1.00000e+07, 298.15, {2.117000e-01, 4.398000e-01, 3.470000e-01, 1.500000e-03}, VAP, -1.021400e+03, 2.300945e-06}, + {1.00000e+07, 298.15, {2.901000e-01, 5.140000e-02, 1.485000e-01, 5.101000e-01}, LIQ, -2.399549e+04, -3.749308e-07}, + {1.01325e+05, 550.00, {2.500000e-01, 2.500000e-01, 2.500000e-01, 2.500000e-01}, VAP, 2.296457e+04, 2.935798e-06}, + {1.00000e+06, 550.00, {2.500000e-01, 2.500000e-01, 2.500000e-01, 2.500000e-01}, VAP, 2.268036e+04, 2.878424e-06}, + {1.00000e+07, 550.00, {2.500000e-01, 2.500000e-01, 2.500000e-01, 2.500000e-01}, VAP, 2.014437e+04, 1.974289e-06}, + } ) +); + +INSTANTIATE_TEST_SUITE_P( + CompositionalEnthalpyTest, SoaveRedlichKwong_6, + ::testing::ValuesIn< typename SoaveRedlichKwong_6::DataType >( { + {1.01325e+05, 283.15, {2.478590e-01, 1.239310e-01, 8.588000e-03, 2.478630e-01, 1.858620e-01, 1.858970e-01}, VAP, -7.928696e+02, 1.221540e-05}, + {1.01325e+05, 283.15, {1.400000e-05, 0.000000e+00, 9.998400e-01, 0.000000e+00, 1.450000e-04, 0.000000e+00}, LIQ, -4.821482e+04, -2.408303e-07}, + {1.00000e+06, 283.15, {2.498020e-01, 1.249180e-01, 1.023000e-03, 2.498360e-01, 1.870430e-01, 1.873770e-01}, VAP, -1.364451e+03, 1.298410e-05}, + {1.00000e+06, 283.15, {1.350000e-04, 0.000000e+00, 9.985240e-01, 1.000000e-06, 1.340000e-03, 0.000000e+00}, LIQ, -4.813521e+04, -2.411439e-07}, + {1.00000e+07, 283.15, {1.606730e-01, 3.112760e-01, 3.509000e-03, 4.560230e-01, 5.605000e-02, 1.247000e-02}, VAP, -2.718582e+03, 4.043744e-06}, + {1.00000e+07, 283.15, {2.139420e-01, 2.509800e-02, 2.696610e-01, 1.092340e-01, 1.833080e-01, 1.987580e-01}, LIQ, -2.059418e+04, -2.359039e-07}, + {1.01325e+05, 298.15, {2.441010e-01, 1.220520e-01, 2.361500e-02, 2.441050e-01, 1.830470e-01, 1.830790e-01}, VAP, -5.567760e+01, 1.118477e-05}, + {1.01325e+05, 298.15, {1.500000e-05, 0.000000e+00, 9.998430e-01, 0.000000e+00, 1.410000e-04, 0.000000e+00}, LIQ, -4.701695e+04, -2.356138e-07}, + {1.00000e+06, 298.15, {2.493660e-01, 1.247020e-01, 2.754000e-03, 2.494030e-01, 1.867230e-01, 1.870520e-01}, VAP, -5.727211e+02, 1.151918e-05}, + {1.00000e+06, 298.15, {1.500000e-04, 0.000000e+00, 9.985160e-01, 2.000000e-06, 1.333000e-03, 0.000000e+00}, LIQ, -4.693758e+04, -2.359289e-07}, + {1.00000e+07, 298.15, {1.885930e-01, 2.808340e-01, 6.328000e-03, 4.344900e-01, 7.191700e-02, 1.783800e-02}, VAP, -2.248771e+03, 4.125190e-06}, + {1.00000e+07, 298.15, {2.048090e-01, 2.376400e-02, 2.816490e-01, 1.011430e-01, 1.829180e-01, 2.057170e-01}, LIQ, -1.967401e+04, -1.819335e-07}, + {1.01325e+05, 550.00, {2.000000e-01, 1.000000e-01, 2.000000e-01, 2.000000e-01, 1.500000e-01, 1.500000e-01}, VAP, 1.424903e+04, 3.202080e-06}, + {1.00000e+06, 550.00, {2.000000e-01, 1.000000e-01, 2.000000e-01, 2.000000e-01, 1.500000e-01, 1.500000e-01}, VAP, 1.405756e+04, 3.173047e-06}, + {1.00000e+07, 550.00, {2.000000e-01, 1.000000e-01, 2.000000e-01, 2.000000e-01, 1.500000e-01, 1.500000e-01}, VAP, 1.216120e+04, 2.686232e-06}, + } ) +); + +// This data is in mass units +INSTANTIATE_TEST_SUITE_P( + CompositionalEnthalpyTest, SoreideWhitson_6, + ::testing::ValuesIn< typename SoreideWhitson_6::DataType >( { + {1.01325e+05, 283.15, {2.478590e-01, 1.239310e-01, 8.588000e-03, 2.478630e-01, 1.858620e-01, 1.858970e-01}, VAP, -2.076518e+04, 1.253196e-05}, + {1.01325e+05, 283.15, {1.400000e-05, 0.000000e+00, 9.998400e-01, 0.000000e+00, 1.450000e-04, 0.000000e+00}, LIQ, -2.598884e+06, -2.346874e-07}, + {1.00000e+06, 283.15, {2.498020e-01, 1.249180e-01, 1.023000e-03, 2.498360e-01, 1.870430e-01, 1.873770e-01}, VAP, -3.590796e+04, 1.327629e-05}, + {1.00000e+06, 283.15, {1.350000e-04, 0.000000e+00, 9.985240e-01, 1.000000e-06, 1.340000e-03, 0.000000e+00}, LIQ, -2.591487e+06, -2.349740e-07}, + {1.00000e+07, 283.15, {1.606730e-01, 3.112760e-01, 3.509000e-03, 4.560230e-01, 5.605000e-02, 1.247000e-02}, VAP, -1.089422e+05, 4.162895e-06}, + {1.00000e+07, 283.15, {2.139420e-01, 2.509800e-02, 2.696610e-01, 1.092340e-01, 1.833080e-01, 1.987580e-01}, LIQ, -5.441170e+05, -2.191882e-07}, + {1.01325e+05, 298.15, {2.441010e-01, 1.220520e-01, 2.361500e-02, 2.441050e-01, 1.830470e-01, 1.830790e-01}, VAP, -1.510413e+03, 1.151375e-05}, + {1.01325e+05, 298.15, {1.500000e-05, 0.000000e+00, 9.998430e-01, 0.000000e+00, 1.410000e-04, 0.000000e+00}, LIQ, -2.537116e+06, -2.293975e-07}, + {1.00000e+06, 298.15, {2.493660e-01, 1.247020e-01, 2.754000e-03, 2.494030e-01, 1.867230e-01, 1.870520e-01}, VAP, -1.534557e+04, 1.183007e-05}, + {1.00000e+06, 298.15, {1.500000e-04, 0.000000e+00, 9.985160e-01, 2.000000e-06, 1.333000e-03, 0.000000e+00}, LIQ, -2.529758e+06, -2.296823e-07}, + {1.00000e+07, 298.15, {1.885930e-01, 2.808340e-01, 6.328000e-03, 4.344900e-01, 7.191700e-02, 1.783800e-02}, VAP, -8.762014e+04, 4.252175e-06}, + {1.00000e+07, 298.15, {2.048090e-01, 2.376400e-02, 2.816490e-01, 1.011430e-01, 1.829180e-01, 2.057170e-01}, LIQ, -5.183653e+05, -1.665260e-07}, + {1.01325e+05, 550.00, {2.000000e-01, 1.000000e-01, 2.000000e-01, 2.000000e-01, 1.500000e-01, 1.500000e-01}, VAP, 4.147675e+05, 3.509469e-06}, + {1.00000e+06, 550.00, {2.000000e-01, 1.000000e-01, 2.000000e-01, 2.000000e-01, 1.500000e-01, 1.500000e-01}, VAP, 4.086661e+05, 3.469178e-06}, + {1.00000e+07, 550.00, {2.000000e-01, 1.000000e-01, 2.000000e-01, 2.000000e-01, 1.500000e-01, 1.500000e-01}, VAP, 3.491882e+05, 2.852633e-06}, + } ) +); + +/* UNCRUSTIFY-ON */ + +} // testing + +} // geos diff --git a/src/coreComponents/constitutive/unitTests/testCubicEOS.cpp b/src/coreComponents/constitutive/unitTests/testCubicEOS.cpp index 74b0d52b115..5b34cb8bdb0 100644 --- a/src/coreComponents/constitutive/unitTests/testCubicEOS.cpp +++ b/src/coreComponents/constitutive/unitTests/testCubicEOS.cpp @@ -93,9 +93,11 @@ class CubicEOSPhaseModelTestFixture : public ::testing::TestWithParam< TestData< void testPureCoefficients( ParamType const & testData ); void testMixtureCoefficients( ParamType const & testData ); + void testATemperatureCoefficients( ParamType const & testData ); void testCompressibilityFactor( ParamType const & testData ); void testCompressibilityFactorValue( ParamType const & testData ); void testLogFugacityCoefficients( ParamType const & testData ); + void testEnthalpy( ParamType const & testData ); protected: std::unique_ptr< TestFluid< NC > > m_fluid{}; @@ -238,6 +240,81 @@ CubicEOSPhaseModelTestFixture< NC, EOS_TYPE >::testMixtureCoefficients( ParamTyp } } +template< integer NC, typename EOS_TYPE > +void +CubicEOSPhaseModelTestFixture< NC, EOS_TYPE >::testATemperatureCoefficients( ParamType const & testData ) +{ + auto componentProperties = this->m_fluid->createKernelWrapper(); + real64 const pressure = std::get< 0 >( testData ); + real64 const temperature = std::get< 1 >( testData ); + stackArray1d< real64, numComps > composition; + TestFluid< numComps >::createArray( composition, std::get< 2 >( testData )); + + auto const binaryInteractionCoefficients = componentProperties.m_componentBinaryCoeff.toSliceConst(); + integer sizes[2] = {0, 0}; + arraySlice2d< real64 const > derivs( nullptr, sizes, sizes ); + typename EOS::template StackVariables< true > stack( numComps, binaryInteractionCoefficients, derivs ); + + EOS::template initialiseStack< true >( + numComps, + pressure, + temperature, + componentProperties, + stack ); + EOS::template computeMixtureCoefficients< 0, true >( numComps, composition.toSliceConst(), stack ); + + integer constexpr numValues = 1; + stackArray1d< real64, numValues > derivatives( numValues ); + + auto concatValues = []( auto const & s, auto & v, real64 const scale = 1.0 ){ + v[0] = scale*s.daMixture[Deriv::dT]; + }; + auto concatDerivatives = []( auto const & s, integer const dof, auto & v, real64 const scale = 1.0 ){ + v[0] = scale*s.ddaMixture[dof]; + }; + + // Pressure derivatives + real64 constexpr pressureScale = 1.0e6; + real64 const dp = 1.0e-4 * pressure; + concatDerivatives( stack, Deriv::dP, derivatives, pressureScale ); + internal::testNumericalDerivative< numValues >( pressure, dp, derivatives.toSliceConst(), [&]( real64 const p, auto & values ) + { + typename EOS::template StackVariables< true > valueStack( numComps, binaryInteractionCoefficients, derivs ); + EOS::template initialiseStack< true >( numComps, p, temperature, componentProperties, valueStack ); + EOS::template computeMixtureCoefficients< 0, true >( numComps, composition.toSliceConst(), valueStack ); + concatValues( valueStack, values, pressureScale ); + }, absTol, relTol ); + + // Temperature derivatives + real64 const dT = 1.0e-4 * temperature; + concatDerivatives( stack, Deriv::dT, derivatives ); + internal::testNumericalDerivative< numValues >( temperature, dT, derivatives.toSliceConst(), [&]( real64 const t, auto & values ) + { + typename EOS::template StackVariables< true > valueStack( numComps, binaryInteractionCoefficients, derivs ); + EOS::template initialiseStack< true >( numComps, pressure, t, componentProperties, valueStack ); + EOS::template computeMixtureCoefficients< 0, true >( numComps, composition.toSliceConst(), valueStack ); + concatValues( valueStack, values ); + }, absTol, relTol ); + + // Composition derivatives + real64 const dz = 1.0e-6; + for( integer ic = 0; ic < numComps; ++ic ) + { + integer const idof = Deriv::dC + ic; + concatDerivatives( stack, idof, derivatives ); + internal::testNumericalDerivative< numValues >( 0, dz, derivatives.toSliceConst(), [&]( real64 const z, auto & values ) + { + real64 const z_orig = composition[ic]; + composition[ic] += z; + typename EOS::template StackVariables< true > valueStack( numComps, binaryInteractionCoefficients, derivs ); + EOS::template initialiseStack< true >( numComps, pressure, temperature, componentProperties, valueStack ); + EOS::template computeMixtureCoefficients< 0, true >( numComps, composition.toSliceConst(), valueStack ); + concatValues( valueStack, values ); + composition[ic] = z_orig; + }, absTol, relTol ); + } +} + template< integer NC, typename EOS_TYPE > void CubicEOSPhaseModelTestFixture< NC, EOS_TYPE >::testCompressibilityFactor( ParamType const & testData ) @@ -410,6 +487,114 @@ CubicEOSPhaseModelTestFixture< NC, EOS_TYPE >::testCompressibilityFactorValue( P checkRelativeError( zFactor, expectedZFactor, relTol, absTol ); } +template< integer NC, typename EOS_TYPE > +void +CubicEOSPhaseModelTestFixture< NC, EOS_TYPE >::testEnthalpy( ParamType const & testData ) +{ + auto componentProperties = this->m_fluid->createKernelWrapper(); + real64 const pressure = std::get< 0 >( testData ); + real64 const temperature = std::get< 1 >( testData ); + stackArray1d< real64, numComps > composition; + TestFluid< numComps >::createArray( composition, std::get< 2 >( testData )); + + constexpr real64 enthalpyScale = 1.0e-2; + + real64 enthalpy = 0.0; + stackArray1d< real64, numDofs > enthalpyDerivs( numDofs ); + + auto const binaryInteractionCoefficients = componentProperties.m_componentBinaryCoeff.toSliceConst(); + integer sizes[2] = {0, 0}; + arraySlice2d< real64 const > derivs( nullptr, sizes, sizes ); + typename EOS::template StackVariables< true > stack( numComps, binaryInteractionCoefficients, derivs ); + + EOS::template initialiseStack< true >( numComps, + pressure, + temperature, + componentProperties, + stack ); + EOS::template computeMixtureCoefficients< 0, true >( numComps, composition.toSliceConst(), stack ); + + EOS::template computeEnthalpy< 0, true >( numComps, + pressure, + temperature, + composition.toSliceConst(), + componentProperties, + stack, + enthalpy, + enthalpyDerivs.toSlice() ); + + // Pressure derivative + real64 const dp = 1.0e-4 * pressure; + internal::testNumericalDerivative( pressure, dp, enthalpyDerivs[Deriv::dP], + [&]( real64 const p ) -> real64 + { + typename EOS::template StackVariables< true > valueStack( numComps, binaryInteractionCoefficients, derivs ); + EOS::template initialiseStack< true >( numComps, p, temperature, componentProperties, valueStack ); + EOS::template computeMixtureCoefficients< 0, true >( numComps, composition.toSliceConst(), valueStack ); + + real64 h = 0.0; + EOS::template computeEnthalpy< 0, false >( numComps, + p, + temperature, + composition.toSliceConst(), + componentProperties, + valueStack, + h, + nullptr ); + return h; + }, absTol, relTol ); + + // Temperature derivative + real64 const dT = 1.0e-4 * temperature; + internal::testNumericalDerivative( temperature, dT, enthalpyScale * enthalpyDerivs[Deriv::dT], + [&]( real64 const t ) -> real64 + { + typename EOS::template StackVariables< true > valueStack( numComps, binaryInteractionCoefficients, derivs ); + EOS::template initialiseStack< true >( numComps, pressure, t, componentProperties, valueStack ); + EOS::template computeMixtureCoefficients< 0, true >( numComps, composition.toSliceConst(), valueStack ); + + real64 h = 0.0; + EOS::template computeEnthalpy< 0, false >( numComps, + pressure, + t, + composition.toSliceConst(), + componentProperties, + valueStack, + h, + nullptr ); + return enthalpyScale * h; + }, absTol, relTol ); + + // Composition derivatives + real64 const dz = 1.0e-6; + for( integer ic = 0; ic < numComps; ++ic ) + { + integer const idof = Deriv::dC + ic; + internal::testNumericalDerivative( 0.0, dz, enthalpyDerivs[idof], + [&]( real64 const z ) -> real64 + { + real64 const z_orig = composition[ic]; + composition[ic] += z; + + typename EOS::template StackVariables< true > valueStack( numComps, binaryInteractionCoefficients, derivs ); + EOS::template initialiseStack< true >( numComps, pressure, temperature, componentProperties, valueStack ); + EOS::template computeMixtureCoefficients< 0, true >( numComps, composition.toSliceConst(), valueStack ); + + real64 h = 0.0; + EOS::template computeEnthalpy< 0, false >( numComps, + pressure, + temperature, + composition.toSliceConst(), + componentProperties, + valueStack, + h, + nullptr ); + composition[ic] = z_orig; + return h; + }, absTol, relTol ); + } +} + using PengRobinson4 = CubicEOSPhaseModelTestFixture< 4, PengRobinsonEOS >; using SoaveRedlichKwong2 = CubicEOSPhaseModelTestFixture< 2, SoaveRedlichKwongEOS >; @@ -418,9 +603,11 @@ TEST_P( PengRobinson4, testCubicModel ) auto const testParam = GetParam(); testPureCoefficients( testParam ); testMixtureCoefficients( testParam ); + testATemperatureCoefficients( testParam ); testCompressibilityFactor( testParam ); testCompressibilityFactorValue( testParam ); testLogFugacityCoefficients( testParam ); + testEnthalpy( testParam ); } TEST_P( SoaveRedlichKwong2, testCubicModel ) @@ -428,6 +615,7 @@ TEST_P( SoaveRedlichKwong2, testCubicModel ) auto const testParam = GetParam(); testPureCoefficients( testParam ); testMixtureCoefficients( testParam ); + testATemperatureCoefficients( testParam ); testCompressibilityFactor( testParam ); testCompressibilityFactorValue( testParam ); testLogFugacityCoefficients( testParam ); diff --git a/src/coreComponents/constitutive/unitTests/testHeatCapacityCoefficients.cpp b/src/coreComponents/constitutive/unitTests/testHeatCapacityCoefficients.cpp new file mode 100644 index 00000000000..fb04fc5314d --- /dev/null +++ b/src/coreComponents/constitutive/unitTests/testHeatCapacityCoefficients.cpp @@ -0,0 +1,215 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2024 TotalEnergies + * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2023-2024 Chevron + * Copyright (c) 2019- GEOS/GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +// Source includes +#include "constitutive/fluid/multifluid/compositional/parameters/HeatCapacityCoefficients.hpp" +#include "constitutive/fluid/multifluid/compositional/parameters/EquationOfState.hpp" +#include "constitutive/fluid/multifluid/MultiFluidBase.hpp" +#include "TestFluid.hpp" + +#include +#include + +using namespace geos::constitutive; +using namespace geos::constitutive::compositional; + +namespace geos +{ +namespace testing +{ + +class MockFluid : public MultiFluidBase +{ +public: + MockFluid( string const & name, Group * const parent ): + MultiFluidBase( name, parent ) + {} + ~MockFluid() override = default; + + string getCatalogName() const override { return "FluidModel"; } + integer getWaterPhaseIndex() const override { return 1; } + void checkTablesParameters( real64, real64 ) const override {} +}; + +class HeatCapacityCoefficientsTestFixture : public ::testing::Test +{ +public: + using Keys = HeatCapacityCoefficients::viewKeyStruct; + static constexpr integer numComps = 4; + static constexpr integer numCoeffs = 5; + static constexpr real64 absTol = 1.0e-10; + +public: + HeatCapacityCoefficientsTestFixture(); + ~HeatCapacityCoefficientsTestFixture() override = default; + +protected: + void createFluid(); +protected: + conduit::Node m_node; + dataRepository::Group m_parent; + std::unique_ptr< MockFluid > m_fluid{}; + std::unique_ptr< TestFluid< numComps > > m_testFluid{}; + std::unique_ptr< ModelParameters > m_parameters{}; + HeatCapacityCoefficients * m_heatCapacityCoefficients{}; +}; + +HeatCapacityCoefficientsTestFixture::HeatCapacityCoefficientsTestFixture(): + m_parent( "parent", m_node ), + m_fluid( std::make_unique< MockFluid >( "fluid", &m_parent )), + m_testFluid( TestFluid< numComps >::create( {Fluid::CO2, Fluid::H2O, Fluid::CH4, Fluid::C5H12} )) +{ + m_parameters = HeatCapacityCoefficients::create( std::move( m_parameters )); + m_heatCapacityCoefficients = dynamic_cast< HeatCapacityCoefficients * >(m_parameters.get()); + GEOS_ERROR_IF( m_heatCapacityCoefficients == nullptr, "Cannot create HeatCapacityCoefficients" ); +} + +void HeatCapacityCoefficientsTestFixture::createFluid() +{ + using FluidKeys = MultiFluidBase::viewKeyStruct; + auto & phaseNames = m_fluid->getWrapper< string_array >( FluidKeys::phaseNamesString() ).reference(); + phaseNames.emplace_back( "gas" ); + phaseNames.emplace_back( "water" ); + phaseNames.emplace_back( "oil" ); + integer const numPhases = static_cast< integer >(phaseNames.size()); + + auto & componentNames = m_fluid->getWrapper< string_array >( FluidKeys::componentNamesString() ).reference(); + for( integer ic = 0; ic < numComps; ++ic ) + { + componentNames.emplace_back( m_testFluid->componentNames[ic] ); + } + + auto & equationsOfState = m_fluid->getWrapper< string_array >( EquationOfState::viewKeyStruct::equationsOfStateString() ).reference(); + string const eosName = EnumStrings< EquationOfStateType >::toString( EquationOfStateType::PengRobinson ); + for( integer ip = 0; ip < numPhases; ++ip ) + { + equationsOfState.emplace_back( eosName ); + } + + m_heatCapacityCoefficients->m_referenceEnthalpy.resize( numPhases, numComps ); + m_heatCapacityCoefficients->m_referenceEnthalpy.zero(); + m_heatCapacityCoefficients->m_coefficients.resize( numComps, numCoeffs ); + m_heatCapacityCoefficients->m_coefficients.zero(); + for( integer ic = 0; ic < numComps; ic++ ) + { + m_heatCapacityCoefficients->m_coefficients( ic, 0 ) = 1.0; + } +} + +TEST_F( HeatCapacityCoefficientsTestFixture, testReferenceTemperature ) +{ + auto componentProperties = m_testFluid->getComponentProperties(); + m_parameters->registerParameters( m_fluid.get()); + createFluid(); + + EXPECT_TRUE( m_fluid->hasWrapper( Keys::enthalpyReferenceTemperatureString()) ); + + real64 & referenceTemperature = m_fluid->getWrapper< real64 >( Keys::enthalpyReferenceTemperatureString() ).reference(); + referenceTemperature = 288.706; + EXPECT_NEAR( m_heatCapacityCoefficients->m_referenceTemperature, referenceTemperature, absTol ); + + referenceTemperature = -1.0; + EXPECT_THROW( m_parameters->postInputInitialization( m_fluid.get(), componentProperties ), InputError ); +} + +TEST_F( HeatCapacityCoefficientsTestFixture, testReferenceEnthalpy ) +{ + auto componentProperties = m_testFluid->getComponentProperties(); + m_parameters->registerParameters( m_fluid.get()); + createFluid(); + + integer const numPhases = m_fluid->numFluidPhases(); + + EXPECT_TRUE( m_fluid->hasWrapper( Keys::referenceEnthalpyString()) ); + + array2d< real64 > & referenceEnthalpy = m_fluid->getWrapper< array2d< real64 > >( Keys::referenceEnthalpyString() ).reference(); + referenceEnthalpy.resize( numPhases - 1, numComps ); + EXPECT_THROW( m_parameters->postInputInitialization( m_fluid.get(), componentProperties ), InputError ); + referenceEnthalpy.resize( numPhases + 1, numComps ); + EXPECT_THROW( m_parameters->postInputInitialization( m_fluid.get(), componentProperties ), InputError ); + referenceEnthalpy.resize( numPhases, numComps - 1 ); + EXPECT_THROW( m_parameters->postInputInitialization( m_fluid.get(), componentProperties ), InputError ); + referenceEnthalpy.resize( numPhases, numComps + 1 ); + EXPECT_THROW( m_parameters->postInputInitialization( m_fluid.get(), componentProperties ), InputError ); + referenceEnthalpy.resize( numPhases, numComps ); + EXPECT_NO_THROW( m_parameters->postInputInitialization( m_fluid.get(), componentProperties ) ); + + constexpr integer vapour = 0; + constexpr integer liquid = 2; + constexpr integer aqueous = 1; + + referenceEnthalpy( vapour, 2 ) = 100.0; + referenceEnthalpy( liquid, 2 ) = 100.0; + referenceEnthalpy( aqueous, 2 ) = 100.0; + + // Reference enthalpy of gas must be greater than oil + referenceEnthalpy( liquid, 2 ) = 110.0; + EXPECT_THROW( m_parameters->postInputInitialization( m_fluid.get(), componentProperties ), InputError ); + referenceEnthalpy( liquid, 2 ) = 100.0; + EXPECT_NO_THROW( m_parameters->postInputInitialization( m_fluid.get(), componentProperties ) ); + + // Reference enthalpy of gas must be greater than water + referenceEnthalpy( aqueous, 2 ) = 120.0; + EXPECT_THROW( m_parameters->postInputInitialization( m_fluid.get(), componentProperties ), InputError ); + referenceEnthalpy( aqueous, 2 ) = 100.0; + EXPECT_NO_THROW( m_parameters->postInputInitialization( m_fluid.get(), componentProperties ) ); + + // Reference enthalpy of water must be greater than oil + referenceEnthalpy( aqueous, 2 ) = 90.0; + EXPECT_THROW( m_parameters->postInputInitialization( m_fluid.get(), componentProperties ), InputError ); + referenceEnthalpy( aqueous, 2 ) = 100.0; + EXPECT_NO_THROW( m_parameters->postInputInitialization( m_fluid.get(), componentProperties ) ); +} + +TEST_F( HeatCapacityCoefficientsTestFixture, testHeatCapacityCoefficients ) +{ + auto componentProperties = m_testFluid->getComponentProperties(); + m_parameters->registerParameters( m_fluid.get()); + createFluid(); + + EXPECT_TRUE( m_fluid->hasWrapper( Keys::componentHeatCapacityCoefficientsString()) ); + + array2d< real64 > & coefficients = m_fluid->getWrapper< array2d< real64 > >( Keys::componentHeatCapacityCoefficientsString() ).reference(); + + auto const setValues = [&]( std::array< real64 const, numCoeffs > const & data ) + { + for( integer ic = 0; ic < numCoeffs; ++ic ) + { + coefficients( 0, ic ) = data[ic]; + } + }; + + setValues( {0.0, 0.0, 0.0, 0.0, 0.0} ); + EXPECT_THROW( m_parameters->postInputInitialization( m_fluid.get(), componentProperties ), InputError ); + setValues( {1.0, 0.0, 0.0, 0.0, 0.0} ); + EXPECT_NO_THROW( m_parameters->postInputInitialization( m_fluid.get(), componentProperties ) ); + + setValues( {9.7529018563e+00, -1.7379753331e-01, 1.3828200497e-03, -4.4470070713e-06, 4.3687436501e-09} ); + EXPECT_THROW( m_parameters->postInputInitialization( m_fluid.get(), componentProperties ), InputError ); + setValues( {1.0, 0.0, 0.0, 0.0, 0.0} ); + EXPECT_NO_THROW( m_parameters->postInputInitialization( m_fluid.get(), componentProperties ) ); + + setValues( {8.440483444351720e-02, -8.156515860472481e-04, 2.921178619177309e-06, -4.649781358999949e-09, 2.775470116151693e-12} ); + EXPECT_THROW( m_parameters->postInputInitialization( m_fluid.get(), componentProperties ), InputError ); + setValues( {1.0, 0.0, 0.0, 0.0, 0.0} ); + EXPECT_NO_THROW( m_parameters->postInputInitialization( m_fluid.get(), componentProperties ) ); + + setValues( {1.469119322002621e-02, -1.822650742125867e-04, 8.537827827668615e-07, -1.777519183416209e-09, 1.387740262103405e-12} ); + EXPECT_NO_THROW( m_parameters->postInputInitialization( m_fluid.get(), componentProperties ) ); +} + +} // namespace testing + +} //namespace geos diff --git a/src/coreComponents/constitutive/unitTests/testSoreideWhitsonEOSPhaseModel.cpp b/src/coreComponents/constitutive/unitTests/testSoreideWhitsonEOSPhaseModel.cpp index 8c56a577404..77e6dd1797c 100644 --- a/src/coreComponents/constitutive/unitTests/testSoreideWhitsonEOSPhaseModel.cpp +++ b/src/coreComponents/constitutive/unitTests/testSoreideWhitsonEOSPhaseModel.cpp @@ -108,9 +108,11 @@ class SoreideWhitsonEOSPhaseModelTestFixture : public ::testing::TestWithParam< void testPureCoefficients( ParamType const & testData ); void testBinaryInteractionCoefficients( ParamType const & testData ); void testMixtureCoefficients( ParamType const & testData ); + void testATemperatureCoefficients( ParamType const & testData ); void testCompressibilityFactor( ParamType const & testData ); void testCompressibilityFactorValue( ParamType const & testData ); void testLogFugacityCoefficients( ParamType const & testData ); + void testEnthalpy( ParamType const & testData ); protected: std::unique_ptr< TestFluid< NC > > m_fluid{}; @@ -325,6 +327,78 @@ SoreideWhitsonEOSPhaseModelTestFixture< NC, EOS_TYPE >::testMixtureCoefficients( } } +template< integer NC, typename EOS_TYPE > +void +SoreideWhitsonEOSPhaseModelTestFixture< NC, EOS_TYPE >::testATemperatureCoefficients( ParamType const & testData ) +{ + auto componentProperties = this->m_fluid->createKernelWrapper(); + real64 const pressure = std::get< 0 >( testData ); + real64 const temperature = std::get< 1 >( testData ); + real64 const salinity = std::get< 2 >( testData ); + stackArray1d< real64, numComps > composition; + TestFluid< numComps >::createArray( composition, std::get< 3 >( testData )); + + typename EOS::template StackVariables< true > stack( numComps ); + + EOS::initialiseStack( numComps, + pressure, + temperature, + componentProperties, + salinity, + stack ); + CubicModel::template computeMixtureCoefficients< 0, true >( numComps, composition.toSliceConst(), stack ); + + integer constexpr numValues = 1; + stackArray1d< real64, numValues > derivatives( numValues ); + + auto concatValues = []( auto const & s, auto & v, real64 const scale = 1.0 ){ + v[0] = scale*s.daMixture[Deriv::dT]; + }; + auto concatDerivatives = []( auto const & s, integer const dof, auto & v, real64 const scale = 1.0 ){ + v[0] = scale*s.ddaMixture[dof]; + }; + + // Pressure derivatives + real64 const dp = 1.0e-4 * pressure; + concatDerivatives( stack, Deriv::dP, derivatives ); + internal::testNumericalDerivative< numValues >( pressure, dp, derivatives.toSliceConst(), [&]( real64 const p, auto & values ) + { + typename EOS::template StackVariables< true > valueStack( numComps ); + EOS::initialiseStack( numComps, p, temperature, componentProperties, salinity, valueStack ); + CubicModel::template computeMixtureCoefficients< 0, true >( numComps, composition.toSliceConst(), valueStack ); + concatValues( valueStack, values ); + }, absTol, relTol ); + + // Temperature derivatives + real64 const dT = 1.0e-4 * temperature; + concatDerivatives( stack, Deriv::dT, derivatives ); + internal::testNumericalDerivative< numValues >( temperature, dT, derivatives.toSliceConst(), [&]( real64 const t, auto & values ) + { + typename EOS::template StackVariables< true > valueStack( numComps ); + EOS::initialiseStack( numComps, pressure, t, componentProperties, salinity, valueStack ); + CubicModel::template computeMixtureCoefficients< 0, true >( numComps, composition.toSliceConst(), valueStack ); + concatValues( valueStack, values ); + }, absTol, relTol ); + + // Composition derivatives + real64 const dz = 1.0e-6; + for( integer ic = 0; ic < numComps; ++ic ) + { + integer const idof = Deriv::dC + ic; + concatDerivatives( stack, idof, derivatives ); + internal::testNumericalDerivative< numValues >( 0, dz, derivatives.toSliceConst(), [&]( real64 const z, auto & values ) + { + real64 const z_orig = composition[ic]; + composition[ic] += z; + typename EOS::template StackVariables< true > valueStack( numComps ); + EOS::initialiseStack( numComps, pressure, temperature, componentProperties, salinity, valueStack ); + CubicModel::template computeMixtureCoefficients< 0, true >( numComps, composition.toSliceConst(), valueStack ); + concatValues( valueStack, values ); + composition[ic] = z_orig; + }, absTol, relTol ); + } +} + template< integer NC, typename EOS_TYPE > void SoreideWhitsonEOSPhaseModelTestFixture< NC, EOS_TYPE >::testCompressibilityFactor( ParamType const & testData ) @@ -498,6 +572,113 @@ SoreideWhitsonEOSPhaseModelTestFixture< NC, EOS_TYPE >::testLogFugacityCoefficie } } +template< integer NC, typename EOS_TYPE > +void +SoreideWhitsonEOSPhaseModelTestFixture< NC, EOS_TYPE >::testEnthalpy( ParamType const & testData ) +{ + auto componentProperties = this->m_fluid->createKernelWrapper(); + real64 const pressure = std::get< 0 >( testData ); + real64 const temperature = std::get< 1 >( testData ); + real64 const salinity = std::get< 2 >( testData ); + StackArray< real64, 1, numComps > composition; + TestFluid< NC >::createArray( composition, std::get< 3 >( testData )); + + constexpr real64 enthalpyScale = 5.0e-4; + + real64 enthalpy = 0.0; + stackArray1d< real64, numDofs > enthalpyDerivs( numDofs ); + + typename EOS::template StackVariables< true > stack( numComps ); + + EOS::initialiseStack( numComps, + pressure, + temperature, + componentProperties, + salinity, + stack ); + CubicModel::template computeMixtureCoefficients< 0, true >( numComps, composition.toSliceConst(), stack ); + + CubicModel::template computeEnthalpy< 0, true >( numComps, + pressure, + temperature, + composition.toSliceConst(), + componentProperties, + stack, + enthalpy, + enthalpyDerivs.toSlice() ); + + // Pressure derivative + real64 const dp = 1.0e-4 * pressure; + internal::testNumericalDerivative( pressure, dp, enthalpyDerivs[Deriv::dP], + [&]( real64 const p ) -> real64 + { + typename EOS::template StackVariables< true > valueStack( numComps ); + EOS::initialiseStack( numComps, p, temperature, componentProperties, salinity, valueStack ); + CubicModel::template computeMixtureCoefficients< 0, true >( numComps, composition.toSliceConst(), valueStack ); + + real64 h = 0.0; + CubicModel::template computeEnthalpy< 0, false >( numComps, + p, + temperature, + composition.toSliceConst(), + componentProperties, + valueStack, + h, + nullptr ); + return h; + }, absTol, relTol ); + + // Temperature derivative + real64 const dT = 1.0e-4 * temperature; + internal::testNumericalDerivative( temperature, dT, enthalpyScale * enthalpyDerivs[Deriv::dT], + [&]( real64 const t ) -> real64 + { + typename EOS::template StackVariables< true > valueStack( numComps ); + EOS::initialiseStack( numComps, pressure, t, componentProperties, salinity, valueStack ); + CubicModel::template computeMixtureCoefficients< 0, true >( numComps, composition.toSliceConst(), valueStack ); + + real64 h = 0.0; + CubicModel::template computeEnthalpy< 0, false >( numComps, + pressure, + t, + composition.toSliceConst(), + componentProperties, + valueStack, + h, + nullptr ); + return enthalpyScale * h; + }, 10*absTol, 10*relTol ); + + // Composition derivatives + real64 const dz = 1.0e-7; + for( integer ic = 0; ic < numComps; ++ic ) + { + integer const idof = Deriv::dC + ic; + internal::testNumericalDerivative( 0.0, dz, enthalpyScale * enthalpyDerivs[idof], + [&]( real64 const z ) -> real64 + { + real64 const z_orig = composition[ic]; + composition[ic] += z; + + typename EOS::template StackVariables< true > valueStack( numComps ); + EOS::initialiseStack( numComps, pressure, temperature, componentProperties, salinity, valueStack ); + CubicModel::template computeMixtureCoefficients< 0, true >( numComps, composition.toSliceConst(), valueStack ); + + real64 h = 0.0; + CubicModel::template computeEnthalpy< 0, false >( numComps, + pressure, + temperature, + composition.toSliceConst(), + componentProperties, + valueStack, + h, + nullptr ); + composition[ic] = z_orig; + return enthalpyScale * h; + }, 10*absTol, relTol ); + } +} + using PengRobinson2 = SoreideWhitsonEOSPhaseModelTestFixture< 2, PengRobinsonEOS >; using PengRobinson4 = SoreideWhitsonEOSPhaseModelTestFixture< 4, PengRobinsonEOS >; using SoaveRedlichKwong3 = SoreideWhitsonEOSPhaseModelTestFixture< 3, SoaveRedlichKwongEOS >; @@ -508,9 +689,11 @@ TEST_P( PengRobinson2, testSWModel ) testPureCoefficients( testParam ); testBinaryInteractionCoefficients( testParam ); testMixtureCoefficients( testParam ); + testATemperatureCoefficients( testParam ); testCompressibilityFactor( testParam ); testCompressibilityFactorValue( testParam ); testLogFugacityCoefficients( testParam ); + testEnthalpy( testParam ); } TEST_P( PengRobinson4, testSWModel ) @@ -535,10 +718,6 @@ TEST_P( SoaveRedlichKwong3, testSWModel ) testLogFugacityCoefficients( testParam ); } -//INSTANTIATE_TEST_SUITE_P( SoreideWhitsonEOSPhaseModelTest, PengRobinson2, ::testing::ValuesIn( generateTestData< 2 >()) ); -//INSTANTIATE_TEST_SUITE_P( SoreideWhitsonEOSPhaseModelTest, PengRobinson4, ::testing::ValuesIn( generateTestData< 4 >()) ); -//INSTANTIATE_TEST_SUITE_P( SoreideWhitsonEOSPhaseModelTest, SoaveRedlichKwong3, ::testing::ValuesIn( generateTestData< 3 >()) ); - /* UNCRUSTIFY-OFF */ INSTANTIATE_TEST_SUITE_P( SoreideWhitsonEOSPhaseModelTest, PengRobinson2,