feat: feature/byer3/thermal well

.gitmodules

@@ -11,5 +11,6 @@
 	path = src/coreComponents/fileIO/coupling/hdf5_interface
 	url = ../../GEOS-DEV/hdf5_interface.git
 [submodule "scripts/uberenv"]
-	path = scripts/uberenv
-	url = ../../LLNL/uberenv.git
+        path = scripts/uberenv
+        url = ../../LLNL/uberenv.git
+

src/coreComponents/common/KernelLaunchSelectors.hpp

@@ -0,0 +1,169 @@
+/*
+ * ------------------------------------------------------------------------------------------------------------
+ * 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 KernelLaunchSelectors.hpp
+ */
+
+#ifndef GEOS_COMMON_KERNELLAUNCHSELECTORS_HPP
+#define GEOS_COMMON_KERNELLAUNCHSELECTORS_HPP
+
+namespace geos
+{
+namespace internal
+{
+
+template< typename T, typename LAMBDA >
+void kernelLaunchSelectorThermalSwitch( T value, LAMBDA && lambda )
+{
+  static_assert( std::is_integral< T >::value, "kernelLaunchSelectorThermalSwitch: type should be integral" );
+
+  switch( value )
+  {
+    case 0:
+    { lambda( std::integral_constant< T, 0 >() ); return; }
+    case 1:
+    { lambda( std::integral_constant< T, 1 >() ); return; }
+
+    default:
+    { GEOS_ERROR( "Unsupported thermal state: " << value ); }
+  }
+}
+
+template< typename T, typename LAMBDA >
+void kernelLaunchSelectorCompThermSwitch( T value, bool const isThermal, LAMBDA && lambda )
+{
+  static_assert( std::is_integral< T >::value, "kernelLaunchSelectorCompSwitch: value type should be integral" );
+
+  //constexpr T a = isThermal ? std::integral_constant< T, 1 >() : std::integral_constant< T, 0 >();
+  if( isThermal )
+  {
+    switch( value )
+    {
+      case 1:
+      {
+        lambda( std::integral_constant< T, 1 >(), std::integral_constant< T, 1 >() ); return;
+      }
+      case 2:
+      { lambda( std::integral_constant< T, 2 >(), std::integral_constant< T, 1 >() ); return; }
+      case 3:
+      { lambda( std::integral_constant< T, 3 >(), std::integral_constant< T, 1 >() ); return; }
+      case 4:
+      { lambda( std::integral_constant< T, 4 >(), std::integral_constant< T, 1 >() ); return; }
+      case 5:
+      { lambda( std::integral_constant< T, 5 >(), std::integral_constant< T, 1 >()); return; }
+      default:
+      { GEOS_ERROR( "Unsupported number of components: " << value ); }
+    }
+  }
+  else
+  {
+    switch( value )
+    {
+      case 1:
+      {
+        lambda( std::integral_constant< T, 1 >(), std::integral_constant< T, 0 >() ); return;
+      }
+      case 2:
+      { lambda( std::integral_constant< T, 2 >(), std::integral_constant< T, 0 >() ); return; }
+      case 3:
+      { lambda( std::integral_constant< T, 3 >(), std::integral_constant< T, 0 >() ); return; }
+      case 4:
+      { lambda( std::integral_constant< T, 4 >(), std::integral_constant< T, 0 >() ); return; }
+      case 5:
+      { lambda( std::integral_constant< T, 5 >(), std::integral_constant< T, 0 >() ); return; }
+      default:
+      { GEOS_ERROR( "Unsupported number of components: " << value ); }
+    }
+  }
+}
+
+
+template< typename T, typename LAMBDA >
+void kernelLaunchSelectorCompPhaseSwitch( T value, T n_phase, LAMBDA && lambda )
+{
+  static_assert( std::is_integral< T >::value, "kernelLaunchSelectorCompSwitch: value type should be integral" );
+
+  //constexpr T a = isThermal ? std::integral_constant< T, 1 >() : std::integral_constant< T, 0 >();
+  if( n_phase == 1 )
+  {
+    switch( value )
+    {
+      case 1:
+      {
+        lambda( std::integral_constant< T, 1 >(), std::integral_constant< T, 1 >() ); return;
+      }
+      case 2:
+      { lambda( std::integral_constant< T, 2 >(), std::integral_constant< T, 1 >() ); return; }
+      case 3:
+      { lambda( std::integral_constant< T, 3 >(), std::integral_constant< T, 1 >() ); return; }
+      case 4:
+      { lambda( std::integral_constant< T, 4 >(), std::integral_constant< T, 1 >() ); return; }
+      case 5:
+      { lambda( std::integral_constant< T, 5 >(), std::integral_constant< T, 1 >() ); return; }
+      default:
+      { GEOS_ERROR( "Unsupported number of components: " << value ); }
+    }
+  }
+  else if( n_phase == 2 )
+  {
+    switch( value )
+    {
+      case 1:
+      {
+        lambda( std::integral_constant< T, 1 >(), std::integral_constant< T, 2 >() ); return;
+      }
+      case 2:
+      { lambda( std::integral_constant< T, 2 >(), std::integral_constant< T, 2 >() ); return; }
+      case 3:
+      { lambda( std::integral_constant< T, 3 >(), std::integral_constant< T, 2 >() ); return; }
+      case 4:
+      { lambda( std::integral_constant< T, 4 >(), std::integral_constant< T, 2 >() ); return; }
+      case 5:
+      { lambda( std::integral_constant< T, 5 >(), std::integral_constant< T, 2 >() ); return; }
+      default:
+      { GEOS_ERROR( "Unsupported number of components: " << value ); }
+    }
+
+  }
+  else if( n_phase == 3 )
+  {
+    switch( value )
+    {
+      case 1:
+      {
+        lambda( std::integral_constant< T, 1 >(), std::integral_constant< T, 3 >()  ); return;
+      }
+      case 2:
+      { lambda( std::integral_constant< T, 2 >(), std::integral_constant< T, 3 >()  ); return; }
+      case 3:
+      { lambda( std::integral_constant< T, 3 >(), std::integral_constant< T, 3 >()  ); return; }
+      case 4:
+      { lambda( std::integral_constant< T, 4 >(), std::integral_constant< T, 3 >() ); return; }
+      case 5:
+      { lambda( std::integral_constant< T, 5 >(), std::integral_constant< T, 3 >() ); return; }
+      default:
+      { GEOS_ERROR( "Unsupported number of components: " << value ); }
+    }
+  }
+  else
+  {
+    { GEOS_ERROR( "Unsupported number of phases: " << n_phase ); }
+  }
+}
+} // end namspace internal
+} // end namespace geos
+
+
+#endif // GEOS_COMMON_KERNELLAUNCHSELECTORS_HPP

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

@@ -380,7 +380,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 @@ class SingleFluidBase : public ConstitutiveBase
   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

@@ -30,6 +30,7 @@
 
 
 #include <iostream>
+#include <typeinfo>
 
 #ifndef NOCHARTOSTRING_KEYLOOKUP
 /// macro definition to enable/disable char * lookups
@@ -455,7 +456,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 );
@@ -596,6 +596,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 )
@@ -879,7 +880,6 @@ class Group
 
     for( auto const & subGroupIter : m_subGroups )
     {
-      std::cout << indent << subGroupIter.second->getName() << std::endl;
       subGroupIter.second->generateDataStructureSkeleton( level + 1 );
     }
   }

src/coreComponents/finiteVolume/TwoPointFluxApproximation.cpp

@@ -130,6 +130,7 @@ void TwoPointFluxApproximation::computeCellStencil( MeshLevel & mesh ) const
                                                                    localIndex const ei,
                                                                    CellElementRegion const & )
   {
+    std::cout << "computeCellStencil " << ei << std::endl;
     regionFilter.insert( ei );
   } );
 

src/coreComponents/linearAlgebra/CMakeLists.txt

@@ -109,6 +109,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 @@ void hypre::mgr::createMGR( LinearSolverParameters const & params,
       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 );