feat: feature/byer3/thermal well

.integrated_tests.yaml

@@ -1,6 +1,6 @@
 baselines:
   bucket: geosx
-  baseline: integratedTests/baseline_integratedTests-pr2878-8188-ed2ded2
+  baseline: integratedTests/baseline_integratedTests-pr3156-8300-32d1d28
 
 allow_fail:
   all: ''

BASELINE_NOTES.md

@@ -6,6 +6,10 @@ 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 #3156 (2024-10-29)
+====================
+Restart check errors due to 1) schema node added to enable thermal option in well model and 2) arrays removed/added for option.  Max difference errors due treatment of shutin wells.  Previously non-zero rate value reported for shutin well, new code will set rate arrays to zero. 
+
 PR #2878 (2024-10-17)
 =====================
 Sorted region cellBlocks names alphabetically. Therefore affected ordering of: faceManager/elemSubRegionList, nodeManager/elemList, nodeManager/elemSubRegionList, SurfaceElementSubRegion::fractureElementsToCellSubRegions, field::perforation::reservoirElementSubregion.

src/coreComponents/constitutive/fluid/multifluid/CO2Brine/CO2BrineFluid.hpp

@@ -386,7 +386,6 @@ CO2BrineFluid< PHASE1, PHASE2, FLASH >::KernelWrapper::
 
   if( m_isThermal )
   {
-
     m_phase1.enthalpy.compute( pressure,
                                temperatureInCelsius,
                                phaseCompFraction.value[ip1].toSliceConst(), phaseCompFraction.derivs[ip1].toSliceConst(),

src/coreComponents/constitutive/fluid/multifluid/Layouts.hpp

@@ -30,10 +30,23 @@ namespace geos
 {
 namespace constitutive
 {
+
+namespace singlefluid
+{
+struct DerivativeOffset
+{
+  /// index of derivative wrt pressure
+  static integer constexpr dP = 0;
+  /// index of derivative wrt temperature
+  static integer constexpr dT = 1;
+
+};
+}
 namespace multifluid
 {
 
 /// indices of pressure, temperature, and composition derivatives
+// Fix me - if the order is changed the code crashes
 struct DerivativeOffset
 {
   /// index of derivative wrt pressure
@@ -44,6 +57,33 @@ struct DerivativeOffset
   static integer constexpr dC = 2;
 };
 
+/// indices of pressure, temperature, and composition derivatives
+template< integer NC, integer IS_THERMAL >
+struct DerivativeOffsetC {};
+
+template< integer NC >
+struct DerivativeOffsetC< NC, 1 >
+{
+  /// index of derivative wrt pressure
+  static integer constexpr dP = 0;
+  /// index of derivative wrt temperature
+  static integer constexpr dT = dP + 1;
+  /// index of first derivative wrt compositions
+  static integer constexpr dC = dP+2;
+  /// number of derivatives
+  static integer constexpr nDer =  NC + 2;
+};
+template< integer NC >
+struct DerivativeOffsetC< NC, 0 >
+{
+  /// index of derivative wrt pressure
+  static integer constexpr dP = 0;
+  /// index of first derivative wrt compositions
+  static integer constexpr dC = dP+1;
+  /// number of derivatives
+  static integer constexpr nDer =  NC + 1;
+};
+
 #if defined( GEOS_USE_DEVICE )
 
 /// Constitutive model phase property array layout

src/coreComponents/constitutive/fluid/singlefluid/SingleFluidBase.hpp

@@ -276,6 +276,14 @@
   virtual real64 defaultDensity() const = 0;
   virtual real64 defaultViscosity() const = 0;
 
+/**
+ * @brief Get the thermal flag.
+ * @return boolean value indicating whether the model can be used to assemble the energy balance equation or not
+ * @detail if isThermal is true, the constitutive model compute the enthalpy and internal energy of the phase.
+ *         This can be used to check the compatibility of the constitutive model with the solver
+ */
+  virtual bool isThermal() const { return false; }
+
 protected:
 
   virtual void postInputInitialization() override;

src/coreComponents/dataRepository/Group.hpp

@@ -457,7 +457,6 @@ class Group
   {
     using T = std::conditional_t< std::is_const< CONTAINERTYPE >::value, CASTTYPE const, CASTTYPE >;
     T * const castedContainer = dynamic_cast< T * >( &container );
-
     if( castedContainer != nullptr )
     {
       lambda( *castedContainer );
@@ -598,6 +597,7 @@ class Group
   void forSubGroups( LOOKUP_CONTAINER const & subGroupKeys, LAMBDA && lambda )
   {
     localIndex counter = 0;
+
     for( auto const & subgroup : subGroupKeys )
     {
       applyLambdaToContainer< GROUPTYPE, GROUPTYPES... >( getGroup( subgroup ), [&]( auto & castedSubGroup )

src/coreComponents/linearAlgebra/CMakeLists.txt

@@ -110,6 +110,7 @@ if( ENABLE_HYPRE )
           interfaces/hypre/mgrStrategies/SinglePhaseReservoirHybridFVM.hpp
           interfaces/hypre/mgrStrategies/SolidMechanicsEmbeddedFractures.hpp
           interfaces/hypre/mgrStrategies/ThermalCompositionalMultiphaseFVM.hpp
+          interfaces/hypre/mgrStrategies/ThermalCompositionalMultiphaseReservoirFVM.hpp
           interfaces/hypre/mgrStrategies/ThermalMultiphasePoromechanics.hpp
           interfaces/hypre/mgrStrategies/ThermalSinglePhasePoromechanics.hpp )
     list( APPEND linearAlgebra_sources

src/coreComponents/linearAlgebra/interfaces/hypre/HypreMGR.cpp

@@ -38,6 +38,7 @@
 #include "linearAlgebra/interfaces/hypre/mgrStrategies/SinglePhaseReservoirFVM.hpp"
 #include "linearAlgebra/interfaces/hypre/mgrStrategies/SinglePhaseReservoirHybridFVM.hpp"
 #include "linearAlgebra/interfaces/hypre/mgrStrategies/ThermalCompositionalMultiphaseFVM.hpp"
+#include "linearAlgebra/interfaces/hypre/mgrStrategies/ThermalCompositionalMultiphaseReservoirFVM.hpp"
 #include "linearAlgebra/interfaces/hypre/mgrStrategies/ThermalSinglePhasePoromechanics.hpp"
 #include "linearAlgebra/interfaces/hypre/mgrStrategies/ThermalMultiphasePoromechanics.hpp"
 #include "linearAlgebra/interfaces/hypre/mgrStrategies/SolidMechanicsEmbeddedFractures.hpp"
@@ -112,6 +113,11 @@
       setStrategy< ThermalCompositionalMultiphaseFVM >( params.mgr, numComponentsPerField, precond, mgrData );
       break;
     }
+    case LinearSolverParameters::MGR::StrategyType::thermalCompositionalMultiphaseReservoirFVM:
+    {
+      setStrategy< ThermalCompositionalMultiphaseReservoirFVM >( params.mgr, numComponentsPerField, precond, mgrData );
+      break;
+    }
     case LinearSolverParameters::MGR::StrategyType::hybridSinglePhasePoromechanics:
     {
       setStrategy< HybridSinglePhasePoromechanics >( params.mgr, numComponentsPerField, precond, mgrData );

src/coreComponents/linearAlgebra/interfaces/hypre/mgrStrategies/ThermalCompositionalMultiphaseReservoirFVM.hpp

@@ -0,0 +1,131 @@
+/*
+ * ------------------------------------------------------------------------------------------------------------
+ * SPDX-License-Identifier: LGPL-2.1-only
+ *
+ * Copyright (c) 2018-2020 Lawrence Livermore National Security LLC
+ * Copyright (c) 2018-2020 The Board of Trustees of the Leland Stanford Junior University
+ * Copyright (c) 2018-2020 TotalEnergies
+ * Copyright (c) 2019-     GEOSX Contributors
+ * All rights reserved
+ *
+ * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details.
+ * ------------------------------------------------------------------------------------------------------------
+ */
+
+/**
+ * @file ThermalCompositionalMultiphaseFVM.hpp
+ */
+
+#ifndef GEOS_LINEARALGEBRA_INTERFACES_HYPREMGRTHERMALCOMPOSITIONALMULTIPHASERESERVOIRFVM_HPP_
+#define GEOS_LINEARALGEBRA_INTERFACES_HYPREMGRTHERMALCOMPOSITIONALMULTIPHASERESERVOIRFVM_HPP_
+
+#include "linearAlgebra/interfaces/hypre/HypreMGR.hpp"
+
+namespace geos
+{
+
+namespace hypre
+{
+
+namespace mgr
+{
+
+/**
+ * @brief ThermalCompositionalMultiphaseReservoirFVM strategy.
+ *
+ * Labels description stored in point_marker_array
+ *               0 = reservoir pressure
+ *               1 = reservoir density (component 1 )
+ *             ... = reservoir densities (# components , NC)
+ *          NC + 1 = reservoir temperature
+ *
+ * numResLabels - 1 = reservoir temperature
+ *     numResLabels = well pressure
+ * numResLabels + 1 = well density
+ *              ... = ... (well densities)
+ * numResLabels + numWellLabels - 3 = last well density
+ * numResLabels + numWellLabels - 2 = well rate
+ * numResLabels + numWellLabels - 1 = well temperature
+ *
+ * 3-level MGR reduction strategy
+ *   - 1st level: eliminate the well block
+ *   - 2nd level: eliminate the reservoir density associated with the volume constraint
+ *   - 3rd level: eliminate the remaining the reservoir densities
+ *   - The coarse grid (pressure and temperature system) is solved with BoomerAMG with numFunctions==2.
+ *
+ */
+
+class ThermalCompositionalMultiphaseReservoirFVM : public MGRStrategyBase< 3 >
+{
+public:
+  /**
+   * @brief Constructor.
+   * @param numComponentsPerField array with number of components for each field
+   */
+  explicit ThermalCompositionalMultiphaseReservoirFVM( arrayView1d< int const > const & numComponentsPerField )
+    : MGRStrategyBase( LvArray::integerConversion< HYPRE_Int >( numComponentsPerField[0] + numComponentsPerField[1] ) )
+  {
+    HYPRE_Int const numResLabels = LvArray::integerConversion< HYPRE_Int >( numComponentsPerField[0] );
+
+    // Level 0: eliminate the well block
+    m_labels[0].resize( numResLabels );
+    std::iota( m_labels[0].begin(), m_labels[0].end(), 0 );
+
+    // Level 1: eliminate last density which corresponds to the volume constraint equation
+    m_labels[1].resize( numResLabels - 2 );
+    std::iota( m_labels[1].begin(), m_labels[1].end(), 0 );
+    m_labels[1].push_back( numResLabels-1 ); // keep temperature
+    // Level 2: eliminate the other densities
+    m_labels[2].push_back( 0 ); // keep pressure
+    m_labels[2].push_back( numResLabels-1 ); // keep temperature
+
+    setupLabels();
+
+    // level 0
+    m_levelFRelaxType[0]          = MGRFRelaxationType::gsElimWInverse;
+    m_levelFRelaxIters[0]         = 1;
+    m_levelInterpType[0]          = MGRInterpolationType::blockJacobi;
+    m_levelRestrictType[0]        = MGRRestrictionType::injection;
+    m_levelCoarseGridMethod[0]    = MGRCoarseGridMethod::galerkin;
+    m_levelGlobalSmootherType[0]  = MGRGlobalSmootherType::none;
+
+    m_levelFRelaxType[1]          = MGRFRelaxationType::jacobi;
+    m_levelFRelaxIters[1]         = 1;
+    m_levelInterpType[1]          = MGRInterpolationType::jacobi; // Diagonal scaling (Jacobi)
+    m_levelRestrictType[1]        = MGRRestrictionType::injection;
+    m_levelCoarseGridMethod[1]    = MGRCoarseGridMethod::galerkin; // Standard Galerkin
+    m_levelGlobalSmootherType[1]  = MGRGlobalSmootherType::blockGaussSeidel;
+    m_levelGlobalSmootherIters[1] = 1;
+
+    m_levelFRelaxType[2]          = MGRFRelaxationType::jacobi;
+    m_levelFRelaxIters[2]         = 1;
+    m_levelInterpType[2]          = MGRInterpolationType::injection; // Injection
+    m_levelRestrictType[2]        = MGRRestrictionType::injection;
+    m_levelCoarseGridMethod[2]    = MGRCoarseGridMethod::cprLikeBlockDiag; // Non-Galerkin Quasi-IMPES CPR
+    m_levelGlobalSmootherType[2]  = MGRGlobalSmootherType::ilu0;
+    m_levelGlobalSmootherIters[2] = 1;
+  }
+
+  /**
+   * @brief Setup the MGR strategy.
+   * @param precond preconditioner wrapper
+   * @param mgrData auxiliary MGR data
+   */
+  void setup( LinearSolverParameters::MGR const &,
+              HyprePrecWrapper & precond,
+              HypreMGRData & mgrData )
+  {
+    setReduction( precond, mgrData );
+
+    // Configure the BoomerAMG solver used as mgr coarse solver for the pressure/temperature reduced system
+    setPressureTemperatureAMG( mgrData.coarseSolver );
+  }
+};
+
+} // namespace mgr
+
+} // namespace hypre
+
+} // namespace geos
+
+#endif /*GEOS_LINEARALGEBRA_INTERFACES_HYPREMGRCOMPOSITIONALMULTIPHASERESERVOIRFVM_HPP_*/

src/coreComponents/linearAlgebra/utilities/LinearSolverParameters.hpp

@@ -271,28 +271,29 @@ struct LinearSolverParameters
      */
     enum class StrategyType : integer
     {
-      invalid,                                     ///< default value, to ensure solver sets something
-      singlePhaseReservoirFVM,                     ///< finite volume single-phase flow with wells
-      singlePhaseHybridFVM,                        ///< hybrid finite volume single-phase flow
-      singlePhaseReservoirHybridFVM,               ///< hybrid finite volume single-phase flow with wells
-      singlePhasePoromechanics,                    ///< single phase poromechanics with finite volume single phase flow
-      thermalSinglePhasePoromechanics,             ///< thermal single phase poromechanics with finite volume single phase flow
-      hybridSinglePhasePoromechanics,              ///< single phase poromechanics with hybrid finite volume single phase flow
-      singlePhasePoromechanicsEmbeddedFractures,   ///< single phase poromechanics with FV embedded fractures
-      singlePhasePoromechanicsConformingFractures, ///< single phase poromechanics with FV conforming  fractures
-      singlePhasePoromechanicsReservoirFVM,        ///< single phase poromechanics with finite volume single phase flow with wells
-      compositionalMultiphaseFVM,                  ///< finite volume compositional multiphase flow
-      compositionalMultiphaseHybridFVM,            ///< hybrid finite volume compositional multiphase flow
-      compositionalMultiphaseReservoirFVM,         ///< finite volume compositional multiphase flow with wells
-      compositionalMultiphaseReservoirHybridFVM,   ///< hybrid finite volume compositional multiphase flow with wells
-      reactiveCompositionalMultiphaseOBL,          ///< finite volume reactive compositional flow with OBL
-      thermalCompositionalMultiphaseFVM,           ///< finite volume thermal compositional multiphase flow
-      multiphasePoromechanics,                     ///< multiphase poromechanics with finite volume compositional multiphase flow
-      multiphasePoromechanicsReservoirFVM,         ///< multiphase poromechanics with finite volume compositional multiphase flow with wells
-      thermalMultiphasePoromechanics,              ///< thermal multiphase poromechanics with finite volume compositional multiphase flow
-      hydrofracture,                               ///< hydrofracture
-      lagrangianContactMechanics,                  ///< Lagrangian contact mechanics
-      solidMechanicsEmbeddedFractures              ///< Embedded fractures mechanics
+      invalid,                                   ///< default value, to ensure solver sets something
+      singlePhaseReservoirFVM,                   ///< finite volume single-phase flow with wells
+      singlePhaseHybridFVM,                      ///< hybrid finite volume single-phase flow
+      singlePhaseReservoirHybridFVM,             ///< hybrid finite volume single-phase flow with wells
+      singlePhasePoromechanics,                  ///< single phase poromechanics with finite volume single phase flow
+      thermalSinglePhasePoromechanics,           ///< thermal single phase poromechanics with finite volume single phase flow
+      hybridSinglePhasePoromechanics,            ///< single phase poromechanics with hybrid finite volume single phase flow
+      singlePhasePoromechanicsEmbeddedFractures, ///< single phase poromechanics with FV embedded fractures
+      singlePhasePoromechanicsConformingFractures, ///< single phase poromechanics with FV embedded fractures
+      singlePhasePoromechanicsReservoirFVM,      ///< single phase poromechanics with finite volume single phase flow with wells
+      compositionalMultiphaseFVM,                ///< finite volume compositional multiphase flow
+      compositionalMultiphaseHybridFVM,          ///< hybrid finite volume compositional multiphase flow
+      compositionalMultiphaseReservoirFVM,       ///< finite volume compositional multiphase flow with wells
+      compositionalMultiphaseReservoirHybridFVM, ///< hybrid finite volume compositional multiphase flow with wells
+      reactiveCompositionalMultiphaseOBL,        ///< finite volume reactive compositional flow with OBL
+      thermalCompositionalMultiphaseFVM,         ///< finite volume thermal compositional multiphase flow
+      thermalCompositionalMultiphaseReservoirFVM,///< finite volume thermal compositional multiphase flow
+      multiphasePoromechanics,                   ///< multiphase poromechanics with finite volume compositional multiphase flow
+      multiphasePoromechanicsReservoirFVM,       ///< multiphase poromechanics with finite volume compositional multiphase flow with wells
+      thermalMultiphasePoromechanics,            ///< thermal multiphase poromechanics with finite volume compositional multiphase flow
+      hydrofracture,                             ///< hydrofracture
+      lagrangianContactMechanics,                ///< Lagrangian contact mechanics
+      solidMechanicsEmbeddedFractures            ///< Embedded fractures mechanics
     };
 
     StrategyType strategy = StrategyType::invalid; ///< Predefined MGR solution strategy (solver specific)
@@ -379,6 +380,7 @@ ENUM_STRINGS( LinearSolverParameters::MGR::StrategyType,
               "compositionalMultiphaseReservoirHybridFVM",
               "reactiveCompositionalMultiphaseOBL",
               "thermalCompositionalMultiphaseFVM",
+              "thermalCompositionalMultiphaseReservoirFVM",
               "multiphasePoromechanics",
               "multiphasePoromechanicsReservoirFVM",
               "thermalMultiphasePoromechanics",

src/coreComponents/mesh/PerforationData.cpp

@@ -39,6 +39,8 @@ PerforationData::PerforationData( string const & name, Group * const parent )
   registerField( fields::perforation::reservoirElementRegion{}, &m_toMeshElements.m_toElementRegion );
   registerField( fields::perforation::reservoirElementSubRegion{}, &m_toMeshElements.m_toElementSubRegion );
   registerField( fields::perforation::reservoirElementIndex{}, &m_toMeshElements.m_toElementIndex );
+  registerField( fields::perforation::reservoirElementGlobalIndex{}, &m_reservoirElementGlobalIndex );
+
   registerField( fields::perforation::wellElementIndex{}, &m_wellElementIndex );
   registerField( fields::perforation::location{}, &m_location );
   registerField( fields::perforation::wellTransmissibility{}, &m_wellTransmissibility );

src/coreComponents/mesh/PerforationData.hpp

@@ -150,6 +150,20 @@ class PerforationData : public ObjectManagerBase
    */
   arrayView1d< localIndex const > getWellElements() const { return m_wellElementIndex; }
 
+  /**
+   * @brief Get perforation-to-reservoir-element connectivity.
+   * @return list of global reservoir element index connected to each perforation
+   */
+  arrayView1d< globalIndex > getReservoirElementGlobalIndex() { return m_reservoirElementGlobalIndex; }
+
+
+
+  /**
+   * @brief Provide an immutable accessor to a const perforation-to-reservoir-element connectivity.
+   * @return list of well element index connected to each perforation
+   */
+  arrayView1d< globalIndex const > getReservoirElementGlobalIndex() const { return m_reservoirElementGlobalIndex; }
+
 
   /**
    * @brief Get perforation locations.
@@ -275,6 +289,9 @@ class PerforationData : public ObjectManagerBase
   /// Indices of the well elements to which perforations are attached
   array1d< localIndex > m_wellElementIndex;
 
+  /// Global indices of reservoir cell containing perforation
+  array1d< globalIndex > m_reservoirElementGlobalIndex;
+
   /// Location of the perforations
   array2d< real64 > m_location;
 

src/coreComponents/mesh/PerforationFields.hpp

@@ -58,6 +58,14 @@ DECLARE_FIELD( reservoirElementIndex,
                WRITE_AND_READ,
                "For each perforation, element index of the perforated element" );
 
+DECLARE_FIELD( reservoirElementGlobalIndex,
+               "reservoirElementGlobalIndex",
+               array1d< globalIndex >,
+               0,
+               NOPLOT,
+               WRITE_AND_READ,
+               "For each perforation, global element index of the perforated element" );
+
 DECLARE_FIELD( wellElementIndex,
                "wellElementIndex",
                array1d< localIndex >,