feat: Add compositional enthalpy model #3705

New issue

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

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

Already on GitHub? Sign in to your account

Open

dkachuma wants to merge 44 commits into develop from dkachuma/compositional-enthalpy

.integrated_tests.yaml

-Original file line number
+Diff line change
@@ -1,6 +1,6 @@
     baselines:
       bucket: geosx
-      baseline: integratedTests/baseline_integratedTests-pr4062-16784-6d8782e
+      baseline: integratedTests/baseline_integratedTests-pr3705-16862-2b262bf
     allow_fail:
       all: ''
@@ Expand Down @@

BASELINE_NOTES.md

-Original file line number
+Diff line change
@@ Expand Up @@
     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) <https://storage.googleapis.com/geosx/integratedTests/baseline_integratedTests-pr3705-16862-2b262bf.tar.gz>
+    Implement compositional enthalpy model
     PR #4062 (2026-05-26) <https://storage.googleapis.com/geosx/integratedTests/baseline_integratedTests-pr4062-16784-6d8782e.tar.gz>
     Add Porous Solid other than PorousElasticity for ALM solver
@@ Expand Down @@

src/coreComponents/constitutive/CMakeLists.txt

-Original file line number
+Diff line change
@@ Expand Up / @@ -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
@@ Expand All / @@ -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
@@ Expand Down Expand Up / @@ -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
@@ Expand All / @@ -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
@@ Expand Down @@

...reComponents/constitutive/fluid/multifluid/compositional/functions/CubicEOSPhaseModel.hpp

-Original file line number
+Diff line change
@@ Expand Up / @@ -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
@@ Expand Down @@

...ponents/constitutive/fluid/multifluid/compositional/functions/CubicEOSPhaseModel_impl.hpp

-Original file line number
+Diff line change
@@ Expand Up / @@ -316,23 +316,41 @@ 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 )
           {
             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] +
+.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];
@@ Expand All / @@ -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;
             }
           }
         }
@@ Expand Down Expand Up / @@ -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
@@ Expand Down @@

...oreComponents/constitutive/fluid/multifluid/compositional/functions/EOSStackVariables.hpp

            
                      Original file line number
                      Diff line number
                      Diff line change
                  
    @@ -145,22 +145,23 @@ 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] ),

        dbic_dt( m_derivativeData[3] ),

        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;

    };

    /**

...reComponents/constitutive/fluid/multifluid/compositional/models/CompositionalEnthalpy.cpp

-Original file line number
+Diff line change
@@ -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

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

feat: Add compositional enthalpy model #3705

Uh oh!

Diff view

Diff view

There are no files selected for viewing

Uh oh!

Uh oh!

feat: Add compositional enthalpy model #3705

Are you sure you want to change the base?

Uh oh!

feat: Add compositional enthalpy model #3705

Uh oh!

Uh oh!

Diff view

Diff view

There are no files selected for viewing

Uh oh!

Uh oh!