feat: Solution scaling for thermal single phase

src/coreComponents/physicsSolvers/fluidFlow/CMakeLists.txt

@@ -20,6 +20,7 @@ Contains:
 set( fluidFlowSolvers_headers
      FlowSolverBase.hpp
      FlowSolverBaseFields.hpp
+     FlowSolverUtilities.hpp
      CompositionalMultiphaseBase.hpp
      CompositionalMultiphaseBaseFields.hpp
      CompositionalMultiphaseStatisticsAggregator.hpp
@@ -67,6 +68,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/CompositionalMultiphaseFVM.cpp

@@ -69,7 +69,7 @@ namespace geos
 using namespace dataRepository;
 using namespace constitutive;
 using namespace fields;
-using namespace compositionalMultiphaseUtilities; // for ScalingType
+using namespace flowSolverUtilities; // for ScalingType
 
 CompositionalMultiphaseFVM::CompositionalMultiphaseFVM( const string & name,
                                                         Group * const parent )
@@ -613,6 +613,14 @@ real64 CompositionalMultiphaseFVM::scalingForSystemSolution( DomainPartition & d
 {
   GEOS_MARK_FUNCTION;
 
+  // In a coupled run the coupled solver's NonlinearSolverParameters overwrite the flow subsolver,
+  // so the lineSearchAction / scalingType compatibility is checked here at run time.
+  // A better fix is to check it at input time in CoupledSolver.
+  GEOS_ERROR_IF( m_scalingType == ScalingType::Local &&
+                 m_nonlinearSolverParameters.m_lineSearchAction != NonlinearSolverParameters::LineSearchAction::None,
+                 "line search is not supported for scalingType = Local; set lineSearchAction to None, including in the coupled solver",
+                 getDataContext() );
+
   if( m_formulationType == CompositionalMultiphaseFormulationType::OverallComposition )
   {
     return scalingForSystemSolutionZFormulation( domain, dofManager, localSolution );
@@ -651,6 +659,7 @@ real64 CompositionalMultiphaseFVM::scalingForSystemSolution( DomainPartition & d
                                      createAndLaunch< parallelDevicePolicy<> >( m_maxRelativePresChange,
                                                                                 m_maxAbsolutePresChange,
                                                                                 m_maxRelativeTempChange,
+                                                                                m_maxAbsoluteTempChange,
                                                                                 m_maxCompFracChange,
                                                                                 m_maxRelativeCompDensChange,
                                                                                 pressure,

src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.hpp

@@ -21,6 +21,7 @@
 #define GEOS_PHYSICSSOLVERS_FLUIDFLOW_COMPOSITIONALMULTIPHASEFVM_HPP_
 
 #include "physicsSolvers/fluidFlow/CompositionalMultiphaseBase.hpp"
+#include "physicsSolvers/fluidFlow/FlowSolverUtilities.hpp"
 
 namespace geos
 {
@@ -229,7 +230,7 @@ class CompositionalMultiphaseFVM : public CompositionalMultiphaseBase
   } m_dbcParams;
 
   /// Solution scaling type
-  compositionalMultiphaseUtilities::ScalingType m_scalingType;
+  flowSolverUtilities::ScalingType m_scalingType;
 
   /// scheme for density treatment in gravity
   GravityDensityScheme m_gravityDensityScheme;

src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseHybridFVM.cpp

@@ -616,7 +616,7 @@ bool CompositionalMultiphaseHybridFVM::checkSystemSolution( DomainPartition & do
           SolutionCheckKernelFactory::
           createAndLaunch< parallelDevicePolicy<> >( m_allowCompDensChopping,
                                                      m_allowNegativePressure,
-                                                     compositionalMultiphaseUtilities::ScalingType::Global,
+                                                     flowSolverUtilities::ScalingType::Global,
                                                      scalingFactor,
                                                      pressure,
                                                      compDens,

src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseUtilities.hpp

@@ -48,19 +48,6 @@ ENUM_STRINGS( GravityDensityScheme,
 namespace compositionalMultiphaseUtilities
 {
 
-/**
- * @brief Solution scaling type, used in CompositionalMultiphaseFVM
- */
-enum class ScalingType : integer
-{
-  Global,         ///< Scale the Newton update with a unique scaling factor
-  Local            ///< Scale the Newton update locally (modifies the Newton direction)
-};
-
-ENUM_STRINGS( ScalingType,
-              "Global",
-              "Local" );
-
 /**
  * @brief In each block, shift the elements from 0 to numRowsToShift-1 one position ahead and replaces the first element
  * with the sum of all elements from 0 to numRowsToShift-1 of the input block one-dimensional array of values.

src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.cpp

@@ -168,6 +168,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"; }
 
@@ -308,6 +309,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/FlowSolverUtilities.hpp

@@ -0,0 +1,49 @@
+/*
+ * ------------------------------------------------------------------------------------------------------------
+ * 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 FlowSolverUtilities.hpp
+ */
+
+#ifndef GEOS_PHYSICSSOLVERS_FLUIDFLOW_FLOWSOLVERUTILITIES_HPP_
+#define GEOS_PHYSICSSOLVERS_FLUIDFLOW_FLOWSOLVERUTILITIES_HPP_
+
+#include "common/DataTypes.hpp"
+#include "common/format/EnumStrings.hpp"
+
+namespace geos
+{
+
+namespace flowSolverUtilities
+{
+
+/**
+ * @brief Solution scaling type
+ */
+enum class ScalingType : integer
+{
+  Global,         ///< Scale the Newton update with a unique scaling factor
+  Local            ///< Scale the Newton update locally (modifies the Newton direction)
+};
+
+ENUM_STRINGS( ScalingType,
+              "Global",
+              "Local" );
+
+} // namespace flowSolverUtilities
+
+} // namespace geos
+
+#endif //GEOS_PHYSICSSOLVERS_FLUIDFLOW_FLOWSOLVERUTILITIES_HPP_

src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp

@@ -44,6 +44,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"
@@ -61,7 +62,8 @@ using namespace singlePhaseBaseKernels;
 
 SinglePhaseBase::SinglePhaseBase( const string & name,
                                   Group * const parent ):
-  FlowSolverBase( name, parent )
+  FlowSolverBase( name, parent ),
+  m_scalingType( flowSolverUtilities::ScalingType::Global )
 {
   this->registerWrapper( viewKeyStruct::inputTemperatureString(), &m_inputTemperature ).
     setApplyDefaultValue( 0.0 ).
@@ -103,9 +105,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() );
       }
     } );
 
@@ -1300,37 +1305,110 @@ real64 SinglePhaseBase::scalingForSystemSolution( DomainPartition & domain,
 {
   GEOS_MARK_FUNCTION;
 
+  // In a coupled run the coupled solver's NonlinearSolverParameters overwrite the flow subsolver,
+  // so the lineSearchAction / scalingType compatibility is checked here at run time.
+  // A better fix is to check it at input time in CoupledSolver.
+  GEOS_ERROR_IF( m_scalingType == flowSolverUtilities::ScalingType::Local &&
+                 m_nonlinearSolverParameters.m_lineSearchAction != NonlinearSolverParameters::LineSearchAction::None,
+                 "line search is not supported for scalingType = Local; set lineSearchAction to None, including in the coupled solver",
+                 getDataContext() );
+
   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 );
 
-      auto const subRegionData =
-        singlePhaseBaseKernels::SolutionScalingKernel::
-          launch< parallelDevicePolicy<> >( localSolution, rankOffset, dofNumber, ghostRank, m_maxAbsolutePresChange );
+        if( m_scalingType == flowSolverUtilities::ScalingType::Global )
+        {
+          scalingFactor = std::min( scalingFactor, subRegionData.scalingFactor );
+        }
+        maxDeltaPres  = std::max( maxDeltaPres, subRegionData.maxDeltaPres );
+        maxDeltaTemp  = std::max( maxDeltaTemp, subRegionData.maxDeltaTemp );
 
-      scalingFactor = std::min( scalingFactor, subRegionData.first );
-      maxDeltaPres  = std::max( maxDeltaPres, subRegionData.second );
+        minPresScalingFactor = std::min( minPresScalingFactor, subRegionData.minPresScalingFactor );
+        minTempScalingFactor = std::min( minTempScalingFactor, subRegionData.minTempScalingFactor );
+      } );
     } );
-  } );
+  }
+  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();
+
+        arrayView1d< real64 > pressureScalingFactor = subRegion.getField< flow::pressureScalingFactor >();
+
+        auto const subRegionData = singlePhaseBaseKernels::
+                                     SolutionScalingKernel::
+                                     launch< parallelDevicePolicy<> >( localSolution, rankOffset, dofNumber, ghostRank,
+                                                                       m_maxAbsolutePresChange, pressureScalingFactor );
+
+        if( m_scalingType == flowSolverUtilities::ScalingType::Global )
+        {
+          scalingFactor = std::min( scalingFactor, subRegionData.scalingFactor );
+        }
+        minPresScalingFactor = std::min( minPresScalingFactor, subRegionData.minPresScalingFactor );
+        maxDeltaPres  = std::max( maxDeltaPres, subRegionData.maxDeltaPres );
+      } );
+    } );
+  }
 
   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;
 }
 
@@ -1358,6 +1436,7 @@ bool SinglePhaseBase::checkSystemSolution( DomainPartition & domain,
       arrayView1d< globalIndex const > const dofNumber = subRegion.getReference< array1d< globalIndex > >( dofKey );
       arrayView1d< integer const > const ghostRank = subRegion.ghostRank();
       arrayView1d< real64 const > const pres = subRegion.getField< flow::pressure >();
+      arrayView1d< real64 const > const pressureScalingFactor = subRegion.getField< flow::pressureScalingFactor >();
       auto const negPresCollector = rankNegPressureIds.createCollector( subRegion.localToGlobalMap().toViewConst() );
 
       auto const statistics = singlePhaseBaseKernels::SolutionCheckKernel::
@@ -1366,7 +1445,9 @@ bool SinglePhaseBase::checkSystemSolution( DomainPartition & domain,
                                                                   dofNumber,
                                                                   ghostRank,
                                                                   pres,
-                                                                  scalingFactor,
+                                                                  m_scalingType,
+                                                                  scalingFactor, // global scaling factor
+                                                                  pressureScalingFactor, // per-element scaling factors
                                                                   negPresCollector );
 
       minPressure = std::min( minPressure, statistics.minNegPres );

src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.hpp

@@ -21,6 +21,7 @@
 #define GEOS_PHYSICSSOLVERS_FLUIDFLOW_SINGLEPHASEBASE_HPP_
 
 #include "physicsSolvers/fluidFlow/FlowSolverBase.hpp"
+#include "physicsSolvers/fluidFlow/FlowSolverUtilities.hpp"
 #include "physicsSolvers/fluidFlow/kernels/singlePhase/AccumulationKernels.hpp"
 #include "physicsSolvers/fluidFlow/kernels/singlePhase/ThermalAccumulationKernels.hpp"
 #include "constitutive/ConstitutiveBase.hpp"
@@ -381,6 +382,9 @@ class SinglePhaseBase : public FlowSolverBase
    */
   virtual void saveConvergedState( ElementSubRegionBase & subRegion ) const override;
 
+  /// Solution scaling type
+  flowSolverUtilities::ScalingType m_scalingType;
+
 private:
 
   virtual void setConstitutiveNames( ElementSubRegionBase & subRegion ) const override;