feat: feature/byer3/thermal well

src/coreComponents/common/KernelLaunchSelectors.hpp

@@ -0,0 +1,174 @@
+/*
+ * ------------------------------------------------------------------------------------------------------------
+ * 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 S, typename T, typename LAMBDA >
+void invokePhaseDispatchLambda ( S val, T numPhases, LAMBDA && lambda )
+{
+  if( numPhases == 1 )
+  {
+    lambda( val, std::integral_constant< T, 1 >());
+    return;
+  }
+  else if( numPhases == 2 )
+  {
+    lambda( val, std::integral_constant< T, 2 >());
+    return;
+  }
+  else if( numPhases == 3 )
+  {
+    lambda( val, std::integral_constant< T, 3 >());
+    return;
+  }
+  else
+  {
+    GEOS_ERROR( "Unsupported state: " << numPhases );
+  }
+};
+
+template< typename S, typename T, typename LAMBDA >
+void invokeThermalDispatchLambda ( S val, T isThermal, LAMBDA && lambda )
+{
+  if( isThermal == 1 )
+  {
+    lambda( val, std::integral_constant< T, 1 >());
+    return;
+  }
+  else if( isThermal == 0 )
+  {
+    lambda( val, std::integral_constant< T, 0 >());
+    return;
+  }
+  else
+  {
+    GEOS_ERROR( "Unsupported state: " << isThermal );
+  }
+};
+
+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" );
+
+
+  switch( value )
+  {
+    case 1:
+    {
+      invokeThermalDispatchLambda( std::integral_constant< T, 1 >(), isThermal, lambda );  return;
+    }
+    case 2:
+    {
+      invokeThermalDispatchLambda( std::integral_constant< T, 2 >(), isThermal, lambda );
+      return;
+    }
+    case 3:
+    {
+      invokeThermalDispatchLambda( std::integral_constant< T, 3 >(), isThermal, lambda );
+      return;
+    }
+    case 4:
+    {
+      invokeThermalDispatchLambda( std::integral_constant< T, 4 >(), isThermal, lambda );
+      return;
+    }
+    case 5:
+    {
+      invokeThermalDispatchLambda( std::integral_constant< T, 5 >(), isThermal, lambda );
+      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" );
+  switch( value )
+  {
+    case 1:
+    {
+      invokePhaseDispatchLambda( std::integral_constant< T, 1 >(), n_phase, lambda );
+      return;
+    }
+    case 2:
+    {
+      invokePhaseDispatchLambda( std::integral_constant< T, 2 >(), n_phase, lambda );
+      return;
+    }
+    case 3:
+    {
+      invokePhaseDispatchLambda( std::integral_constant< T, 3 >(), n_phase, lambda );
+      return;
+    }
+    case 4:
+    {
+      invokePhaseDispatchLambda( std::integral_constant< T, 4 >(), n_phase, lambda );
+      return;
+    }
+    case 5:
+    {
+      invokePhaseDispatchLambda( std::integral_constant< T, 5 >(), n_phase, lambda );
+      return;
+    }
+    default:
+    { GEOS_ERROR( "Unsupported number of components: " << value ); }
+  }
+}
+
+
+} // 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

@@ -455,7 +455,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 +595,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

@@ -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 );