feat: Solution scaling for thermal single phase

src/coreComponents/physicsSolvers/fluidFlow/CMakeLists.txt

@@ -63,6 +63,7 @@ set( fluidFlowSolvers_headers
      kernels/singlePhase/ThermalAccumulationKernels.hpp
      kernels/singlePhase/ThermalDirichletFluxComputeKernel.hpp
      kernels/singlePhase/ThermalFluxComputeKernel.hpp
+     kernels/singlePhase/ThermalSolutionScalingKernel.hpp
      kernels/singlePhase/proppant/ProppantBaseKernels.hpp
      kernels/singlePhase/proppant/ProppantFluxKernels.hpp
      kernels/singlePhase/reactive/AccumulationKernels.hpp

src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.cpp

@@ -130,6 +130,12 @@ FlowSolverBase::FlowSolverBase( string const & name,
     setApplyDefaultValue( -1.0 ).       // disabled by default
     setDescription( "Maximum (absolute) pressure change in a Newton iteration" );
 
+  this->registerWrapper( viewKeyStruct::maxAbsoluteTempChangeString(), &m_maxAbsoluteTempChange ).
+    setSizedFromParent( 0 ).
+    setInputFlag( InputFlags::OPTIONAL ).
+    setApplyDefaultValue( -1.0 ).       // disabled by default
+    setDescription( "Maximum (absolute) temperature change in a Newton iteration" );
+
   this->registerWrapper( viewKeyStruct::maxSequentialPresChangeString(), &m_maxSequentialPresChange ).
     setSizedFromParent( 0 ).
     setInputFlag( InputFlags::OPTIONAL ).

src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.hpp

@@ -77,6 +77,7 @@ class FlowSolverBase : public PhysicsSolverBase
     static constexpr char const * inputTemperatureString() { return "temperature"; }
     static constexpr char const * allowNegativePressureString() { return "allowNegativePressure"; }
     static constexpr char const * maxAbsolutePresChangeString() { return "maxAbsolutePressureChange"; }
+    static constexpr char const * maxAbsoluteTempChangeString() { return "maxAbsoluteTemperatureChange"; }
     static constexpr char const * maxSequentialPresChangeString() { return "maxSequentialPressureChange"; }
     static constexpr char const * maxSequentialTempChangeString() { return "maxSequentialTemperatureChange"; }
 
@@ -276,6 +277,9 @@ class FlowSolverBase : public PhysicsSolverBase
   /// maximum (absolute) pressure change in a Newton iteration
   real64 m_maxAbsolutePresChange;
 
+  /// maximum (absolute) temperature change in a Newton iteration
+  real64 m_maxAbsoluteTempChange;
+
   /// maximum (absolute) pressure change in a sequential iteration
   real64 m_sequentialPresChange;
   real64 m_maxSequentialPresChange;

src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp

@@ -42,6 +42,7 @@
 #include "physicsSolvers/fluidFlow/kernels/singlePhase/MobilityKernel.hpp"
 #include "physicsSolvers/fluidFlow/kernels/singlePhase/SolutionCheckKernel.hpp"
 #include "physicsSolvers/fluidFlow/kernels/singlePhase/SolutionScalingKernel.hpp"
+#include "physicsSolvers/fluidFlow/kernels/singlePhase/ThermalSolutionScalingKernel.hpp"
 #include "physicsSolvers/fluidFlow/kernels/singlePhase/StatisticsKernel.hpp"
 #include "physicsSolvers/fluidFlow/kernels/singlePhase/HydrostaticPressureKernel.hpp"
 #include "physicsSolvers/fluidFlow/kernels/singlePhase/ThermalHydrostaticPressureKernel.hpp"
@@ -101,9 +102,12 @@ void SinglePhaseBase::registerDataOnMesh( Group & meshBodies )
       subRegion.registerField< flow::mass_n >( getName() );
       subRegion.registerField< flow::dMass >( getName() ).reference().resizeDimension< 1 >( m_numDofPerCell );
 
+      subRegion.registerField< flow::pressureScalingFactor >( getName() );
+
       if( m_isThermal )
       {
         subRegion.registerField< flow::dEnergy >( getName() ).reference().resizeDimension< 1 >( m_numDofPerCell );
+        subRegion.registerField< flow::temperatureScalingFactor >( getName() );
       }
     } );
 
@@ -1288,34 +1292,90 @@ real64 SinglePhaseBase::scalingForSystemSolution( DomainPartition & domain,
   string const dofKey = dofManager.getKey( viewKeyStruct::elemDofFieldString() );
   real64 scalingFactor = 1.0;
   real64 maxDeltaPres = 0.0;
+  real64 maxDeltaTemp = 0.0;
 
-  forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &,
-                                                               MeshLevel & mesh,
-                                                               string_array const & regionNames )
+  real64 minPresScalingFactor = 1.0, minTempScalingFactor = 1.0;
+
+  if( m_isThermal )
   {
-    mesh.getElemManager().forElementSubRegions( regionNames,
-                                                [&]( localIndex const,
-                                                     ElementSubRegionBase & subRegion )
+    forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &,
+                                                                 MeshLevel & mesh,
+                                                                 string_array const & regionNames )
     {
-      globalIndex const rankOffset = dofManager.rankOffset();
-      arrayView1d< globalIndex const > const dofNumber = subRegion.getReference< array1d< globalIndex > >( dofKey );
-      arrayView1d< integer const > const ghostRank = subRegion.ghostRank();
+      mesh.getElemManager().forElementSubRegions( regionNames,
+                                                  [&]( localIndex const,
+                                                       ElementSubRegionBase & subRegion )
+      {
+        globalIndex const rankOffset = dofManager.rankOffset();
+        arrayView1d< globalIndex const > const dofNumber = subRegion.getReference< array1d< globalIndex > >( dofKey );
+        arrayView1d< integer const > const ghostRank = subRegion.ghostRank();
+
+        arrayView1d< real64 > pressureScalingFactor = subRegion.getField< flow::pressureScalingFactor >();
+        arrayView1d< real64 > temperatureScalingFactor = subRegion.getField< flow::temperatureScalingFactor >();
+
+        integer const tempDofOffset = 1;
+        auto const subRegionData = thermalSinglePhaseBaseKernels::
+                                     SolutionScalingKernel::
+                                     launch< parallelDevicePolicy<> >( localSolution, rankOffset, tempDofOffset, dofNumber, ghostRank,
+                                                                       m_maxAbsolutePresChange, m_maxAbsoluteTempChange,
+                                                                       pressureScalingFactor, temperatureScalingFactor );
+
+        scalingFactor = std::min( scalingFactor, std::get< 0 >( subRegionData ) );
+        maxDeltaPres  = std::max( maxDeltaPres, std::get< 1 >( subRegionData ) );
+        maxDeltaTemp  = std::max( maxDeltaTemp, std::get< 2 >( subRegionData ) );
+
+        minPresScalingFactor = std::min( minPresScalingFactor, std::get< 3 >( subRegionData ) );
+        minTempScalingFactor = std::min( minTempScalingFactor, std::get< 4 >( subRegionData ) );
+      } );
+    } );
+  }
+  else
+  {
+    GEOS_UNUSED_VAR( maxDeltaTemp );
+
+    forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &,
+                                                                 MeshLevel & mesh,
+                                                                 string_array const & regionNames )
+    {
+      mesh.getElemManager().forElementSubRegions( regionNames,
+                                                  [&]( localIndex const,
+                                                       ElementSubRegionBase & subRegion )
+      {
+        globalIndex const rankOffset = dofManager.rankOffset();
+        arrayView1d< globalIndex const > const dofNumber = subRegion.getReference< array1d< globalIndex > >( dofKey );
+        arrayView1d< integer const > const ghostRank = subRegion.ghostRank();
 
-      auto const subRegionData =
-        singlePhaseBaseKernels::SolutionScalingKernel::
-          launch< parallelDevicePolicy<> >( localSolution, rankOffset, dofNumber, ghostRank, m_maxAbsolutePresChange );
+        auto const subRegionData = singlePhaseBaseKernels::
+                                     SolutionScalingKernel::
+                                     launch< parallelDevicePolicy<> >( localSolution, rankOffset, dofNumber, ghostRank, m_maxAbsolutePresChange );
 
-      scalingFactor = std::min( scalingFactor, subRegionData.first );
-      maxDeltaPres  = std::max( maxDeltaPres, subRegionData.second );
+        scalingFactor = std::min( scalingFactor, subRegionData.first );
+        minPresScalingFactor = std::min( minPresScalingFactor, subRegionData.first );
+        maxDeltaPres  = std::max( maxDeltaPres, subRegionData.second );
+      } );
     } );
-  } );
+  }
 
   scalingFactor = MpiWrapper::min( scalingFactor );
+  minPresScalingFactor = MpiWrapper::min( minPresScalingFactor );
   maxDeltaPres  = MpiWrapper::max( maxDeltaPres );
 
   GEOS_LOG_LEVEL_RANK_0( logInfo::Solution, GEOS_FMT( "        {}: Max pressure change = {} Pa (before scaling)",
                                                       getName(), fmt::format( "{:.{}f}", maxDeltaPres, 3 ) ) );
 
+  GEOS_LOG_LEVEL_RANK_0( logInfo::Solution, GEOS_FMT( "        {}: Min pressure scaling factor = {}", getName(), minPresScalingFactor ) );
+
+  if( m_isThermal )
+  {
+    minTempScalingFactor = MpiWrapper::min( minTempScalingFactor );
+    maxDeltaTemp = MpiWrapper::max( maxDeltaTemp );
+
+    GEOS_LOG_LEVEL_RANK_0( logInfo::Solution, GEOS_FMT( "        {}: Max temperature change = {} K (before scaling)",
+                                                        getName(), fmt::format( "{:.{}f}", maxDeltaTemp, 3 ) ) );
+
+    GEOS_LOG_LEVEL_RANK_0( logInfo::Solution, GEOS_FMT( "        {}: Min temperature scaling factor = {}", getName(), minTempScalingFactor ) );
+  }
+
   return scalingFactor;
 }
 

src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseFVM.cpp

@@ -277,13 +277,13 @@ void SinglePhaseFVM< BASE >::applySystemSolution( DofManager const & dofManager,
     dofManager.addVectorToField( localSolution,
                                  BASE::viewKeyStruct::elemDofFieldString(),
                                  flow::pressure::key(),
-                                 scalingFactor,
+                                 flow::pressureScalingFactor::key(),
                                  pressureMask );
 
     dofManager.addVectorToField( localSolution,
                                  BASE::viewKeyStruct::elemDofFieldString(),
                                  flow::temperature::key(),
-                                 scalingFactor,
+                                 flow::temperatureScalingFactor::key(),
                                  temperatureMask );
   }
   else

src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/ThermalSolutionScalingKernel.hpp

@@ -0,0 +1,103 @@
+/*
+ * ------------------------------------------------------------------------------------------------------------
+ * 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 ThermalSolutionScalingKernel.hpp
+ */
+
+#ifndef GEOS_PHYSICSSOLVERS_FLUIDFLOW_SINGLEPHASE_THERMALSOLUTIONSCALINGKERNEL_HPP
+#define GEOS_PHYSICSSOLVERS_FLUIDFLOW_SINGLEPHASE_THERMALSOLUTIONSCALINGKERNEL_HPP
+
+#include "common/DataTypes.hpp"
+#include "common/GEOS_RAJA_Interface.hpp"
+
+namespace geos
+{
+
+namespace thermalSinglePhaseBaseKernels
+{
+
+/******************************** SolutionScalingKernel ********************************/
+
+struct SolutionScalingKernel
+{
+  template< typename POLICY >
+  static std::tuple< real64, real64, real64, real64, real64 >  launch( arrayView1d< real64 const > const & localSolution,
+                                                                       globalIndex const rankOffset,
+                                                                       globalIndex const temperatureOffset,
+                                                                       arrayView1d< globalIndex const > const & dofNumber,
+                                                                       arrayView1d< integer const > const & ghostRank,
+                                                                       real64 const maxAbsolutePresChange,
+                                                                       real64 const maxAbsoluteTempChange,
+                                                                       arrayView1d< real64 > pressureScalingFactor,
+                                                                       arrayView1d< real64 > temperatureScalingFactor )
+  {
+    RAJA::ReduceMin< ReducePolicy< POLICY >, real64 > scalingFactor( 1.0 );
+    RAJA::ReduceMax< ReducePolicy< POLICY >, real64 > maxDeltaPres( 0.0 );
+    RAJA::ReduceMax< ReducePolicy< POLICY >, real64 > maxDeltaTemp( 0.0 );
+
+    RAJA::ReduceMin< ReducePolicy< POLICY >, real64 > localMinPresScalingFactor( 1.0 );
+    RAJA::ReduceMin< ReducePolicy< POLICY >, real64 > localMinTempScalingFactor( 1.0 );
+
+    forAll< POLICY >( dofNumber.size(), [=] GEOS_HOST_DEVICE ( localIndex const ei ) mutable
+    {
+      if( ghostRank[ei] < 0 && dofNumber[ei] >= 0 )
+      {
+        pressureScalingFactor[ei] = 1.0;
+        temperatureScalingFactor[ei] = 1.0;
+
+        localIndex const lid = dofNumber[ei] - rankOffset;
+
+        // compute the change in pressure
+        real64 const absPresChange = LvArray::math::abs( localSolution[lid] );
+        maxDeltaPres.max( absPresChange );
+
+        // compute the change in temperature
+        real64 const absTempChange = LvArray::math::abs( localSolution[lid + temperatureOffset] );
+        maxDeltaTemp.max( absTempChange );
+
+        // maxAbsolutePresChange <= 0.0 means that scaling is disabled, and we are only collecting maxDeltaPres in this kernel
+        if( maxAbsolutePresChange > 0.0 && absPresChange > maxAbsolutePresChange )
+        {
+          real64 const presScalingFactor = maxAbsolutePresChange / absPresChange;
+          pressureScalingFactor[ei] = presScalingFactor;
+          scalingFactor.min( presScalingFactor );
+
+          localMinPresScalingFactor.min( presScalingFactor );
+        }
+
+        // maxAbsoluteTempChange <= 0.0 means that scaling is disabled, and we are only collecting maxDeltaTemps in this kernel
+        if( maxAbsoluteTempChange > 0.0 && absTempChange > maxAbsoluteTempChange )
+        {
+          real64 const tempScalingFactor = maxAbsoluteTempChange / absTempChange;
+          temperatureScalingFactor[ei] = tempScalingFactor;
+          scalingFactor.min( tempScalingFactor );
+
+          localMinTempScalingFactor.min( tempScalingFactor );
+        }
+      }
+
+    } );
+
+    return { scalingFactor.get(), maxDeltaPres.get(), maxDeltaTemp.get(), localMinPresScalingFactor.get(), localMinTempScalingFactor.get() };
+  }
+
+};
+
+} // namespace thermalSinglePhaseBaseKernels
+
+} // namespace geos
+
+#endif //GEOS_PHYSICSSOLVERS_FLUIDFLOW_SINGLEPHASE_THERMALSOLUTIONSCALINGKERNEL_HPP

src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.cpp

@@ -69,7 +69,6 @@ CompositionalMultiphaseWell::CompositionalMultiphaseWell( const string & name,
   m_useTotalMassEquation( 1 ),
   m_maxCompFracChange( 1.0 ),
   m_maxRelativePresChange( 0.2 ),
-  m_maxAbsolutePresChange( -1 ), // disabled by default
   m_minScalingFactor( 0.01 ),
   m_allowCompDensChopping( 1 ),
   m_targetPhaseIndex( -1 )
@@ -102,12 +101,6 @@ CompositionalMultiphaseWell::CompositionalMultiphaseWell( const string & name,
     setApplyDefaultValue( 1.0 ).
     setDescription( "Maximum (relative) change in pressure between two Newton iterations (recommended with rate control)" );
 
-  this->registerWrapper( viewKeyStruct::maxAbsolutePresChangeString(), &m_maxAbsolutePresChange ).
-    setSizedFromParent( 0 ).
-    setInputFlag( InputFlags::OPTIONAL ).
-    setApplyDefaultValue( -1.0 ).       // disabled by default
-    setDescription( "Maximum (absolute) pressure change in a Newton iteration" );
-
   this->registerWrapper( viewKeyStruct::maxRelativeTempChangeString(), &m_maxRelativeTempChange ).
     setSizedFromParent( 0 ).
     setInputFlag( InputFlags::OPTIONAL ).

src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.hpp

@@ -263,8 +263,6 @@ class CompositionalMultiphaseWell : public WellSolverBase
 
     static constexpr char const * maxRelativePresChangeString() { return "maxRelativePressureChange"; }
 
-    static constexpr char const * maxAbsolutePresChangeString() { return "maxAbsolutePressureChange"; }
-
     static constexpr char const * maxRelativeCompDensChangeString() { return "maxRelativeCompDensChange"; }
 
     static constexpr char const * maxRelativeTempChangeString() { return "maxRelativeTemperatureChange"; }
@@ -380,9 +378,6 @@ class CompositionalMultiphaseWell : public WellSolverBase
   /// maximum (relative) change in pressure between two Newton iterations
   real64 m_maxRelativePresChange;
 
-  /// maximum (absolute) change in pressure between two Newton iterations
-  real64 m_maxAbsolutePresChange;
-
   /// maximum (relative) change in component density between two Newton iterations
   real64 m_maxRelativeCompDensChange;