From c694c12861afff777b9aab3812a438ff50b8ca0c Mon Sep 17 00:00:00 2001
From: frankfeifan <fei2@llnl.gov>
Date: Sun, 22 Feb 2026 13:03:08 -0800
Subject: [PATCH 01/15] add scaling output for thermal single phase

---
 .../ThermalCompressibleSinglePhaseFluid.hpp   |   2 +-
 .../physicsSolvers/fluidFlow/CMakeLists.txt   |   1 +
 .../fluidFlow/FlowSolverBase.cpp              |   6 ++
 .../fluidFlow/FlowSolverBase.hpp              |   4 +
 .../fluidFlow/SinglePhaseBase.cpp             |  89 ++++++++++++---
 .../ThermalSolutionScalingKernel.hpp          | 102 ++++++++++++++++++
 .../wells/CompositionalMultiphaseWell.cpp     |   7 --
 .../wells/CompositionalMultiphaseWell.hpp     |   5 -
 .../fluidFlow/wells/SinglePhaseWell.cpp       | 102 ++++++++++++++++++
 .../fluidFlow/wells/SinglePhaseWell.hpp       |   5 +
 .../fluidFlow/wells/WellSolverBase.cpp        |  12 +++
 .../fluidFlow/wells/WellSolverBase.hpp        |   9 ++
 12 files changed, 316 insertions(+), 28 deletions(-)
 create mode 100644 src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/ThermalSolutionScalingKernel.hpp

diff --git a/src/coreComponents/constitutive/fluid/singlefluid/ThermalCompressibleSinglePhaseFluid.hpp b/src/coreComponents/constitutive/fluid/singlefluid/ThermalCompressibleSinglePhaseFluid.hpp
index 2f110c295ef..4051580aeb1 100644
--- a/src/coreComponents/constitutive/fluid/singlefluid/ThermalCompressibleSinglePhaseFluid.hpp
+++ b/src/coreComponents/constitutive/fluid/singlefluid/ThermalCompressibleSinglePhaseFluid.hpp
@@ -229,7 +229,7 @@ class ThermalCompressibleSinglePhaseFluid : public CompressibleSinglePhaseFluid
   using CompressibleSinglePhaseFluid::m_densityModelType;
 
   /// Type of kernel wrapper for in-kernel update (TODO: support multiple EAT combinations, not just this combination)
-  using KernelWrapper = ThermalCompressibleSinglePhaseUpdate< ExponentApproximationType::Full, ExponentApproximationType::Linear, ExponentApproximationType::Linear >;
+  using KernelWrapper = ThermalCompressibleSinglePhaseUpdate< ExponentApproximationType::Linear, ExponentApproximationType::Linear, ExponentApproximationType::Linear >;
 
   /**
    * @brief Create an update kernel wrapper.
diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CMakeLists.txt b/src/coreComponents/physicsSolvers/fluidFlow/CMakeLists.txt
index 6eb55e9775a..063b5e5359e 100644
--- a/src/coreComponents/physicsSolvers/fluidFlow/CMakeLists.txt
+++ b/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
diff --git a/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.cpp b/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.cpp
index bbc5a8466c3..5183914d873 100644
--- a/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.cpp
+++ b/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 ).
diff --git a/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.hpp b/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.hpp
index 8b68a5fad08..10943be607f 100644
--- a/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.hpp
+++ b/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;
diff --git a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp
index 713367c05f0..6800a1f9597 100644
--- a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp
+++ b/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,33 +1292,88 @@ 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();
 
-      auto const subRegionData =
-        singlePhaseBaseKernels::SolutionScalingKernel::
-          launch< parallelDevicePolicy<> >( localSolution, rankOffset, dofNumber, ghostRank, m_maxAbsolutePresChange );
+        arrayView1d< real64 > pressureScalingFactor = subRegion.getField< flow::pressureScalingFactor >();
+        arrayView1d< real64 > temperatureScalingFactor = subRegion.getField< flow::temperatureScalingFactor >();
 
-      scalingFactor = std::min( scalingFactor, subRegionData.first );
-      maxDeltaPres  = std::max( maxDeltaPres, subRegionData.second );
+        auto const subRegionData = thermalSinglePhaseBaseKernels::
+                                     SolutionScalingKernel::
+                                     launch< parallelDevicePolicy<> >( localSolution, rankOffset, 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 );
+
+        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;
 }
diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/ThermalSolutionScalingKernel.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/ThermalSolutionScalingKernel.hpp
new file mode 100644
index 00000000000..4f73ed51e22
--- /dev/null
+++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/ThermalSolutionScalingKernel.hpp
@@ -0,0 +1,102 @@
+/*
+ * ------------------------------------------------------------------------------------------------------------
+ * 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,
+                                                                       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 + 1] );
+        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
diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.cpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.cpp
index 9f18e9741da..9cdb88b6fc4 100644
--- a/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.cpp
+++ b/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 ).
diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.hpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.hpp
index cb8ed9f1a6f..c4b7ea6d328 100644
--- a/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.hpp
+++ b/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;
 
diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/SinglePhaseWell.cpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/SinglePhaseWell.cpp
index 6414693bf22..d4ae96e9677 100644
--- a/src/coreComponents/physicsSolvers/fluidFlow/wells/SinglePhaseWell.cpp
+++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/SinglePhaseWell.cpp
@@ -43,6 +43,8 @@
 #include "physicsSolvers/fluidFlow/wells/kernels/SinglePhasePerforationFluxKernels.hpp"
 #include "physicsSolvers/fluidFlow/kernels/singlePhase/FluidUpdateKernel.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/SinglePhaseStatistics.hpp"
 
 namespace geos
@@ -87,12 +89,16 @@ void SinglePhaseWell::registerDataOnMesh( Group & meshBodies )
       subRegion.registerField< well::connectionRate_n >( getName() );
       subRegion.registerField< well::connectionRate >( getName() );
 
+      subRegion.registerField< well::pressureScalingFactor >( getName() );
+
       PerforationData & perforationData = *subRegion.getPerforationData();
       perforationData.registerField< well::perforationRate >( getName() );
       perforationData.registerField< well::dPerforationRate >( getName() ).
         reference().resizeDimension< 1, 2 >( 2, 2 );
       if( isThermal() )
       {
+        subRegion.registerField< well::temperatureScalingFactor >( getName() );
+        
         perforationData.registerField< well::energyPerforationFlux >( getName() );
         perforationData.registerField< well::dEnergyPerforationFlux >( getName() ).
           reference().resizeDimension< 1, 2 >( 2, 2 );
@@ -1008,6 +1014,102 @@ SinglePhaseWell::calculateResidualNorm( real64 const & time_n,
   return resNorm;
 }
 
+real64
+SinglePhaseWell::scalingForSystemSolution( DomainPartition & domain,
+                                           DofManager const & dofManager,
+                                           arrayView1d< real64 const > const & localSolution )
+{
+  GEOS_MARK_FUNCTION;
+
+  string const wellDofKey = dofManager.getKey( wellElementDofName() );
+
+  real64 scalingFactor = 1.0;
+  real64 maxDeltaPres = 0.0, maxDeltaTemp = 0.0;
+
+  real64 minPresScalingFactor = 1.0, minTempScalingFactor = 1.0;
+
+  if( m_isThermal )
+  {
+    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 > >( wellDofKey );
+        arrayView1d< integer const > const ghostRank = subRegion.ghostRank();
+
+        arrayView1d< real64 > pressureScalingFactor = subRegion.getField< well::pressureScalingFactor >();
+        arrayView1d< real64 > temperatureScalingFactor = subRegion.getField< well::temperatureScalingFactor >();
+
+        auto const subRegionData = thermalSinglePhaseBaseKernels::
+                                     SolutionScalingKernel::
+                                     launch< parallelDevicePolicy<> >( localSolution, rankOffset, 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
+  {
+    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 > >( wellDofKey );
+        arrayView1d< integer const > const ghostRank = subRegion.ghostRank();
+
+        auto const subRegionData = singlePhaseBaseKernels::
+                                     SolutionScalingKernel::
+                                     launch< parallelDevicePolicy<> >( localSolution, rankOffset, dofNumber, ghostRank, m_maxAbsolutePresChange );
+
+        scalingFactor = std::min( subRegionData.first, scalingFactor );
+        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 well pressure change: {} Pa (before scaling)",
+                                   getName(), GEOS_FMT( "{:.{}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 well temperature change: {} K (before scaling)",
+                                     getName(), GEOS_FMT( "{:.{}f}", maxDeltaTemp, 3 ) ) );
+
+    GEOS_LOG_LEVEL_RANK_0( logInfo::Solution, GEOS_FMT( "        {}: Min temperature scaling factor = {}", getName(), minTempScalingFactor ) );
+  }
+
+  return scalingFactor;
+
+}
+
 bool SinglePhaseWell::checkSystemSolution( DomainPartition & domain,
                                            DofManager const & dofManager,
                                            arrayView1d< real64 const > const & localSolution,
diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/SinglePhaseWell.hpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/SinglePhaseWell.hpp
index 8c07e223fb7..e9cb1df4699 100644
--- a/src/coreComponents/physicsSolvers/fluidFlow/wells/SinglePhaseWell.hpp
+++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/SinglePhaseWell.hpp
@@ -102,6 +102,11 @@ class SinglePhaseWell : public WellSolverBase
                          DofManager const & dofManager,
                          arrayView1d< real64 const > const & localRhs ) override;
 
+  virtual real64
+  scalingForSystemSolution( DomainPartition & domain,
+                            DofManager const & dofManager,
+                            arrayView1d< real64 const > const & localSolution ) override;
+
   virtual bool
   checkSystemSolution( DomainPartition & domain,
                        DofManager const & dofManager,
diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/WellSolverBase.cpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/WellSolverBase.cpp
index 2e98c2fbaf2..3bdbd075b6a 100644
--- a/src/coreComponents/physicsSolvers/fluidFlow/wells/WellSolverBase.cpp
+++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/WellSolverBase.cpp
@@ -63,6 +63,18 @@ WellSolverBase::WellSolverBase( string const & name,
     setInputFlag( dataRepository::InputFlags::OPTIONAL ).
     setDescription( "Choose time step to honor rates/bhp tables time intervals" );
 
+  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::maxAbsoluteTempChangeString(), &m_maxAbsoluteTempChange ).
+    setSizedFromParent( 0 ).
+    setInputFlag( InputFlags::OPTIONAL ).
+    setApplyDefaultValue( -1.0 ).       // disabled by default
+    setDescription( "Maximum (absolute) temperature change in a Newton iteration" );
+
   addLogLevel< logInfo::WellControl >();
 }
 
diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/WellSolverBase.hpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/WellSolverBase.hpp
index 04fe58112b4..e9536fc429c 100644
--- a/src/coreComponents/physicsSolvers/fluidFlow/wells/WellSolverBase.hpp
+++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/WellSolverBase.hpp
@@ -289,6 +289,9 @@ class WellSolverBase : public PhysicsSolverBase
     static constexpr char const * timeStepFromTablesFlagString() { return "timeStepFromTables"; }
 
     static constexpr char const * fluidNamesString() { return "fluidNames"; }
+
+    static constexpr char const * maxAbsolutePresChangeString() { return "maxAbsolutePressureChange"; }
+    static constexpr char const * maxAbsoluteTempChangeString() { return "maxAbsoluteTemperatureChange"; }
   };
 
 private:
@@ -357,6 +360,12 @@ class WellSolverBase : public PhysicsSolverBase
 
   /// name of the fluid constitutive model used as a reference for component/phase description
   string m_referenceFluidModelName;
+
+  /// maximum (absolute) change in pressure between two Newton iterations
+  real64 m_maxAbsolutePresChange;
+
+  /// maximum (absolute) change in temperature between two Newton iterations
+  real64 m_maxAbsoluteTempChange;
 };
 
 }

From 7049f106ae66c9ef905ec65aeab3ceb7f541c4a2 Mon Sep 17 00:00:00 2001
From: frankfeifan <fei2@llnl.gov>
Date: Sun, 22 Feb 2026 14:35:21 -0800
Subject: [PATCH 02/15] apply scaling

---
 .../physicsSolvers/fluidFlow/SinglePhaseFVM.cpp               | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseFVM.cpp b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseFVM.cpp
index bff2fbd0f78..127f1a3a40f 100644
--- a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseFVM.cpp
+++ b/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

From 5b9042b92c6a2c5dd4186c241054d1874e03eb4b Mon Sep 17 00:00:00 2001
From: frankfeifan <fei2@llnl.gov>
Date: Sun, 22 Feb 2026 16:35:00 -0800
Subject: [PATCH 03/15] temp dof offset fix

---
 .../physicsSolvers/fluidFlow/SinglePhaseBase.cpp               | 2 +-
 .../kernels/singlePhase/ThermalSolutionScalingKernel.hpp       | 3 ++-
 .../physicsSolvers/fluidFlow/wells/SinglePhaseWell.cpp         | 2 +-
 3 files changed, 4 insertions(+), 3 deletions(-)

diff --git a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp
index 6800a1f9597..a24225e1e00 100644
--- a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp
+++ b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp
@@ -1315,7 +1315,7 @@ real64 SinglePhaseBase::scalingForSystemSolution( DomainPartition & domain,
 
         auto const subRegionData = thermalSinglePhaseBaseKernels::
                                      SolutionScalingKernel::
-                                     launch< parallelDevicePolicy<> >( localSolution, rankOffset, dofNumber, ghostRank,
+                                     launch< parallelDevicePolicy<> >( localSolution, rankOffset, 1, dofNumber, ghostRank,
                                                                        m_maxAbsolutePresChange, m_maxAbsoluteTempChange,
                                                                        pressureScalingFactor, temperatureScalingFactor );
 
diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/ThermalSolutionScalingKernel.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/ThermalSolutionScalingKernel.hpp
index 4f73ed51e22..5ffc3ff7ed7 100644
--- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/ThermalSolutionScalingKernel.hpp
+++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/ThermalSolutionScalingKernel.hpp
@@ -36,6 +36,7 @@ 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,
@@ -64,7 +65,7 @@ struct SolutionScalingKernel
         maxDeltaPres.max( absPresChange );
 
         // compute the change in temperature
-        real64 const absTempChange = LvArray::math::abs( localSolution[lid + 1] );
+        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
diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/SinglePhaseWell.cpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/SinglePhaseWell.cpp
index d4ae96e9677..d7869034822 100644
--- a/src/coreComponents/physicsSolvers/fluidFlow/wells/SinglePhaseWell.cpp
+++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/SinglePhaseWell.cpp
@@ -1047,7 +1047,7 @@ SinglePhaseWell::scalingForSystemSolution( DomainPartition & domain,
 
         auto const subRegionData = thermalSinglePhaseBaseKernels::
                                      SolutionScalingKernel::
-                                     launch< parallelDevicePolicy<> >( localSolution, rankOffset, dofNumber, ghostRank,
+                                     launch< parallelDevicePolicy<> >( localSolution, rankOffset, 2, dofNumber, ghostRank,
                                                                        m_maxAbsolutePresChange, m_maxAbsoluteTempChange,
                                                                        pressureScalingFactor, temperatureScalingFactor );
 

From fff67fd23e21dda394aab8574a0c4fe374519481 Mon Sep 17 00:00:00 2001
From: frankfeifan <fei2@llnl.gov>
Date: Mon, 23 Feb 2026 13:35:25 -0800
Subject: [PATCH 04/15] apply scalingFactor to well solutions

---
 .../physicsSolvers/fluidFlow/wells/SinglePhaseWell.cpp        | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/SinglePhaseWell.cpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/SinglePhaseWell.cpp
index d7869034822..f90384dc90b 100644
--- a/src/coreComponents/physicsSolvers/fluidFlow/wells/SinglePhaseWell.cpp
+++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/SinglePhaseWell.cpp
@@ -1177,7 +1177,7 @@ SinglePhaseWell::applySystemSolution( DofManager const & dofManager,
   dofManager.addVectorToField( localSolution,
                                wellElementDofName(),
                                well::pressure::key(),
-                               scalingFactor,
+                               well::pressureScalingFactor::key(),
                                pressureMask );
 
   dofManager.addVectorToField( localSolution,
@@ -1193,7 +1193,7 @@ SinglePhaseWell::applySystemSolution( DofManager const & dofManager,
     dofManager.addVectorToField( localSolution,
                                  wellElementDofName(),
                                  fields::well::temperature::key(),
-                                 scalingFactor,
+                                 well::temperatureScalingFactor::key(),
                                  temperatureMask );
 
   }

From d20f61ff742f38113a91f553fc916b3347753ab1 Mon Sep 17 00:00:00 2001
From: frankfeifan <fei2@llnl.gov>
Date: Mon, 23 Feb 2026 13:54:50 -0800
Subject: [PATCH 05/15] uncrustify, and restore the acccidentally changed fluid
 model option

---
 .../singlefluid/ThermalCompressibleSinglePhaseFluid.hpp    | 2 +-
 .../physicsSolvers/fluidFlow/SinglePhaseBase.cpp           | 7 ++++---
 .../physicsSolvers/fluidFlow/wells/SinglePhaseWell.cpp     | 5 +++--
 3 files changed, 8 insertions(+), 6 deletions(-)

diff --git a/src/coreComponents/constitutive/fluid/singlefluid/ThermalCompressibleSinglePhaseFluid.hpp b/src/coreComponents/constitutive/fluid/singlefluid/ThermalCompressibleSinglePhaseFluid.hpp
index 4051580aeb1..2f110c295ef 100644
--- a/src/coreComponents/constitutive/fluid/singlefluid/ThermalCompressibleSinglePhaseFluid.hpp
+++ b/src/coreComponents/constitutive/fluid/singlefluid/ThermalCompressibleSinglePhaseFluid.hpp
@@ -229,7 +229,7 @@ class ThermalCompressibleSinglePhaseFluid : public CompressibleSinglePhaseFluid
   using CompressibleSinglePhaseFluid::m_densityModelType;
 
   /// Type of kernel wrapper for in-kernel update (TODO: support multiple EAT combinations, not just this combination)
-  using KernelWrapper = ThermalCompressibleSinglePhaseUpdate< ExponentApproximationType::Linear, ExponentApproximationType::Linear, ExponentApproximationType::Linear >;
+  using KernelWrapper = ThermalCompressibleSinglePhaseUpdate< ExponentApproximationType::Full, ExponentApproximationType::Linear, ExponentApproximationType::Linear >;
 
   /**
    * @brief Create an update kernel wrapper.
diff --git a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp
index a24225e1e00..5fbf1cdbfc5 100644
--- a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp
+++ b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp
@@ -1313,9 +1313,10 @@ real64 SinglePhaseBase::scalingForSystemSolution( DomainPartition & domain,
         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, 1, dofNumber, ghostRank,
+                                     launch< parallelDevicePolicy<> >( localSolution, rankOffset, tempDofOffset, dofNumber, ghostRank,
                                                                        m_maxAbsolutePresChange, m_maxAbsoluteTempChange,
                                                                        pressureScalingFactor, temperatureScalingFactor );
 
@@ -1361,7 +1362,7 @@ real64 SinglePhaseBase::scalingForSystemSolution( DomainPartition & domain,
 
   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 )
@@ -1371,7 +1372,7 @@ real64 SinglePhaseBase::scalingForSystemSolution( DomainPartition & domain,
 
     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 ) );
   }
 
diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/SinglePhaseWell.cpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/SinglePhaseWell.cpp
index f90384dc90b..417e11caa0c 100644
--- a/src/coreComponents/physicsSolvers/fluidFlow/wells/SinglePhaseWell.cpp
+++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/SinglePhaseWell.cpp
@@ -98,7 +98,7 @@ void SinglePhaseWell::registerDataOnMesh( Group & meshBodies )
       if( isThermal() )
       {
         subRegion.registerField< well::temperatureScalingFactor >( getName() );
-        
+
         perforationData.registerField< well::energyPerforationFlux >( getName() );
         perforationData.registerField< well::dEnergyPerforationFlux >( getName() ).
           reference().resizeDimension< 1, 2 >( 2, 2 );
@@ -1045,9 +1045,10 @@ SinglePhaseWell::scalingForSystemSolution( DomainPartition & domain,
         arrayView1d< real64 > pressureScalingFactor = subRegion.getField< well::pressureScalingFactor >();
         arrayView1d< real64 > temperatureScalingFactor = subRegion.getField< well::temperatureScalingFactor >();
 
+        integer const tempDofOffset = 2;
         auto const subRegionData = thermalSinglePhaseBaseKernels::
                                      SolutionScalingKernel::
-                                     launch< parallelDevicePolicy<> >( localSolution, rankOffset, 2, dofNumber, ghostRank,
+                                     launch< parallelDevicePolicy<> >( localSolution, rankOffset, tempDofOffset, dofNumber, ghostRank,
                                                                        m_maxAbsolutePresChange, m_maxAbsoluteTempChange,
                                                                        pressureScalingFactor, temperatureScalingFactor );
 

From 51a189a5f4639116258963769e9f2c4a9bd33a82 Mon Sep 17 00:00:00 2001
From: frankfeifan <frankfei0929@gmail.com>
Date: Fri, 9 Oct 2026 22:17:20 -0700
Subject: [PATCH 06/15] schema update

---
 src/coreComponents/schema/schema.xsd | 28 ++++++++++++++++++++++++++--
 1 file changed, 26 insertions(+), 2 deletions(-)

diff --git a/src/coreComponents/schema/schema.xsd b/src/coreComponents/schema/schema.xsd
index c719a937ad1..1ad66a96106 100644
--- a/src/coreComponents/schema/schema.xsd
+++ b/src/coreComponents/schema/schema.xsd
@@ -3299,6 +3299,8 @@ Information output from lower logLevels is added with the desired log level
 		<xsd:attribute name="logLevel" type="integer" default="0" />
 		<!--maxAbsolutePressureChange => Maximum (absolute) pressure change in a Newton iteration-->
 		<xsd:attribute name="maxAbsolutePressureChange" type="real64" default="-1" />
+		<!--maxAbsoluteTemperatureChange => Maximum (absolute) temperature change in a Newton iteration-->
+		<xsd:attribute name="maxAbsoluteTemperatureChange" type="real64" default="-1" />
 		<!--maxCompFractionChange => Maximum (absolute) allowed change in the composition fraction of any component within a single time step, in a Newton iteration-->
 		<xsd:attribute name="maxCompFractionChange" type="real64" default="0.5" />
 		<!--maxRelativeCompDensChange => Maximum (relative) change in a component density in a Newton iteration-->
@@ -3420,6 +3422,8 @@ Information output from lower logLevels is added with the desired log level
 		<xsd:attribute name="logLevel" type="integer" default="0" />
 		<!--maxAbsolutePressureChange => Maximum (absolute) pressure change in a Newton iteration-->
 		<xsd:attribute name="maxAbsolutePressureChange" type="real64" default="-1" />
+		<!--maxAbsoluteTemperatureChange => Maximum (absolute) temperature change in a Newton iteration-->
+		<xsd:attribute name="maxAbsoluteTemperatureChange" type="real64" default="-1" />
 		<!--maxCompFractionChange => Maximum (absolute) allowed change in the composition fraction of any component within a single time step, in a Newton iteration-->
 		<xsd:attribute name="maxCompFractionChange" type="real64" default="0.5" />
 		<!--maxRelativeCompDensChange => Maximum (relative) change in a component density in a Newton iteration-->
@@ -4138,6 +4142,8 @@ Information output from lower logLevels is added with the desired log level
 		<xsd:attribute name="logLevel" type="integer" default="0" />
 		<!--maxAbsolutePressureChange => Maximum (absolute) pressure change in a Newton iteration-->
 		<xsd:attribute name="maxAbsolutePressureChange" type="real64" default="-1" />
+		<!--maxAbsoluteTemperatureChange => Maximum (absolute) temperature change in a Newton iteration-->
+		<xsd:attribute name="maxAbsoluteTemperatureChange" type="real64" default="-1" />
 		<!--maxSequentialPressureChange => Maximum (absolute) pressure change in a sequential iteration, used for outer loop convergence check-->
 		<xsd:attribute name="maxSequentialPressureChange" type="real64" default="100000" />
 		<!--maxSequentialTemperatureChange => Maximum (absolute) temperature change in a sequential iteration, used for outer loop convergence check-->
@@ -4713,6 +4719,8 @@ Information output from lower logLevels is added with the desired log level
 		<xsd:attribute name="logLevel" type="integer" default="0" />
 		<!--maxAbsolutePressureChange => Maximum (absolute) pressure change in a Newton iteration-->
 		<xsd:attribute name="maxAbsolutePressureChange" type="real64" default="-1" />
+		<!--maxAbsoluteTemperatureChange => Maximum (absolute) temperature change in a Newton iteration-->
+		<xsd:attribute name="maxAbsoluteTemperatureChange" type="real64" default="-1" />
 		<!--maxProppantConcentration => Maximum proppant concentration-->
 		<xsd:attribute name="maxProppantConcentration" type="real64" default="0.6" />
 		<!--maxSequentialPressureChange => Maximum (absolute) pressure change in a sequential iteration, used for outer loop convergence check-->
@@ -4777,6 +4785,8 @@ Information output from lower logLevels is added with the desired log level
 		<xsd:attribute name="logLevel" type="integer" default="0" />
 		<!--maxAbsolutePressureChange => Maximum (absolute) pressure change in a Newton iteration-->
 		<xsd:attribute name="maxAbsolutePressureChange" type="real64" default="-1" />
+		<!--maxAbsoluteTemperatureChange => Maximum (absolute) temperature change in a Newton iteration-->
+		<xsd:attribute name="maxAbsoluteTemperatureChange" type="real64" default="-1" />
 		<!--maxCompFractionChange => Maximum (absolute) change in the composition fraction of any component within a single time step, in a Newton iteration-->
 		<xsd:attribute name="maxCompFractionChange" type="real64" default="1" />
 		<!--maxSequentialPressureChange => Maximum (absolute) pressure change in a sequential iteration, used for outer loop convergence check-->
@@ -4888,6 +4898,8 @@ Information output from lower logLevels is added with the desired log level
 		<xsd:attribute name="logLevel" type="integer" default="0" />
 		<!--maxAbsolutePressureChange => Maximum (absolute) pressure change in a Newton iteration-->
 		<xsd:attribute name="maxAbsolutePressureChange" type="real64" default="-1" />
+		<!--maxAbsoluteTemperatureChange => Maximum (absolute) temperature change in a Newton iteration-->
+		<xsd:attribute name="maxAbsoluteTemperatureChange" type="real64" default="-1" />
 		<!--maxSequentialPressureChange => Maximum (absolute) pressure change in a sequential iteration, used for outer loop convergence check-->
 		<xsd:attribute name="maxSequentialPressureChange" type="real64" default="100000" />
 		<!--maxSequentialTemperatureChange => Maximum (absolute) temperature change in a sequential iteration, used for outer loop convergence check-->
@@ -4943,6 +4955,8 @@ Information output from lower logLevels is added with the desired log level
 		<xsd:attribute name="logLevel" type="integer" default="0" />
 		<!--maxAbsolutePressureChange => Maximum (absolute) pressure change in a Newton iteration-->
 		<xsd:attribute name="maxAbsolutePressureChange" type="real64" default="-1" />
+		<!--maxAbsoluteTemperatureChange => Maximum (absolute) temperature change in a Newton iteration-->
+		<xsd:attribute name="maxAbsoluteTemperatureChange" type="real64" default="-1" />
 		<!--maxSequentialPressureChange => Maximum (absolute) pressure change in a sequential iteration, used for outer loop convergence check-->
 		<xsd:attribute name="maxSequentialPressureChange" type="real64" default="100000" />
 		<!--maxSequentialTemperatureChange => Maximum (absolute) temperature change in a sequential iteration, used for outer loop convergence check-->
@@ -5296,6 +5310,8 @@ Information output from lower logLevels is added with the desired log level
 		<xsd:attribute name="logLevel" type="integer" default="0" />
 		<!--maxAbsolutePressureChange => Maximum (absolute) pressure change in a Newton iteration-->
 		<xsd:attribute name="maxAbsolutePressureChange" type="real64" default="-1" />
+		<!--maxAbsoluteTemperatureChange => Maximum (absolute) temperature change in a Newton iteration-->
+		<xsd:attribute name="maxAbsoluteTemperatureChange" type="real64" default="-1" />
 		<!--maxSequentialPressureChange => Maximum (absolute) pressure change in a sequential iteration, used for outer loop convergence check-->
 		<xsd:attribute name="maxSequentialPressureChange" type="real64" default="100000" />
 		<!--maxSequentialTemperatureChange => Maximum (absolute) temperature change in a sequential iteration, used for outer loop convergence check-->
@@ -5357,6 +5373,8 @@ Information output from lower logLevels is added with the desired log level
 		<xsd:attribute name="logLevel" type="integer" default="0" />
 		<!--maxAbsolutePressureChange => Maximum (absolute) pressure change in a Newton iteration-->
 		<xsd:attribute name="maxAbsolutePressureChange" type="real64" default="-1" />
+		<!--maxAbsoluteTemperatureChange => Maximum (absolute) temperature change in a Newton iteration-->
+		<xsd:attribute name="maxAbsoluteTemperatureChange" type="real64" default="-1" />
 		<!--maxSequentialPressureChange => Maximum (absolute) pressure change in a sequential iteration, used for outer loop convergence check-->
 		<xsd:attribute name="maxSequentialPressureChange" type="real64" default="100000" />
 		<!--maxSequentialTemperatureChange => Maximum (absolute) temperature change in a sequential iteration, used for outer loop convergence check-->
@@ -6223,6 +6241,8 @@ Information output from lower logLevels is added with the desired log level
 		<xsd:attribute name="logLevel" type="integer" default="0" />
 		<!--maxAbsolutePressureChange => Maximum (absolute) pressure change in a Newton iteration-->
 		<xsd:attribute name="maxAbsolutePressureChange" type="real64" default="-1" />
+		<!--maxAbsoluteTemperatureChange => Maximum (absolute) temperature change in a Newton iteration-->
+		<xsd:attribute name="maxAbsoluteTemperatureChange" type="real64" default="-1" />
 		<!--maxCompFractionChange => Maximum (absolute) allowed change in the composition fraction of any component within a single time step, in a Newton iteration-->
 		<xsd:attribute name="maxCompFractionChange" type="real64" default="1" />
 		<!--maxRelativeCompDensChange => Maximum (relative) change in a component density between two Newton iterations-->
@@ -6473,6 +6493,10 @@ Information output from lower logLevels is added with the desired log level
 1
  - Well control information-->
 		<xsd:attribute name="logLevel" type="integer" default="0" />
+		<!--maxAbsolutePressureChange => Maximum (absolute) pressure change in a Newton iteration-->
+		<xsd:attribute name="maxAbsolutePressureChange" type="real64" default="-1" />
+		<!--maxAbsoluteTemperatureChange => Maximum (absolute) temperature change in a Newton iteration-->
+		<xsd:attribute name="maxAbsoluteTemperatureChange" type="real64" default="-1" />
 		<!--referenceReservoirRegion => Name of reservoir region used for obtaining average region pressure used in volume rate constraint calculations.
 Frequency of pressure update is set in SinglePhase/CompositionalMultiphaseStatistics definition.
 Setting cycleFrequency='1' will update the pressure every timestep, note that is a lagged property in constraint propertiesNote the event associated with the statists task must be entered before the solver event.
@@ -8830,8 +8854,8 @@ If you want to do a three-phase simulation, please use instead wettingIntermedia
 		<xsd:attribute name="referencePressure" type="real64" default="0" />
 		<!--referenceViscosity => Reference fluid viscosity-->
 		<xsd:attribute name="referenceViscosity" type="real64" default="0.001" />
-		<!--solventMassPerSolutionVolume => Mass of solvent per unit volume of solution [kg/m^3], used to convert species molality [mol/kg solvent] to molarity [mol/m^3 solution]. The default of 1000 approximates an aqueous solution by the density of pure water.-->
-		<xsd:attribute name="solventMassPerSolutionVolume" type="real64" default="1000" />
+		<!--solventMassFraction => Mass fraction of solvent in the solution [-]. Species molality [mol/kg solvent] times this fraction is the amount per kg of solution, and times the fluid density the amount per m^3 of solution. The default of 1 is a dilute aqueous solution.-->
+		<xsd:attribute name="solventMassFraction" type="real64" default="1" />
 		<!--viscosibility => Fluid viscosity exponential coefficient-->
 		<xsd:attribute name="viscosibility" type="real64" default="0" />
 		<!--viscosityModelType => Type of viscosity model. Valid options:

From 2bcdfd7a6c5d43e1ed5fe8e5ccfab0a0f731817d Mon Sep 17 00:00:00 2001
From: frankfeifan <frankfei0929@gmail.com>
Date: Sun, 11 Oct 2026 20:14:08 -0700
Subject: [PATCH 07/15] add per element pressure scaling for the isothermal

---
 .../fluidFlow/SinglePhaseBase.cpp             |  5 +++-
 .../singlePhase/SolutionScalingKernel.hpp     |  6 ++++-
 .../fluidFlow/wells/SinglePhaseWell.cpp       | 26 ++++++++++++-------
 3 files changed, 25 insertions(+), 12 deletions(-)

diff --git a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp
index 22bfc687bf5..c686dc677cc 100644
--- a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp
+++ b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp
@@ -1360,9 +1360,12 @@ real64 SinglePhaseBase::scalingForSystemSolution( DomainPartition & domain,
         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 );
+                                     launch< parallelDevicePolicy<> >( localSolution, rankOffset, dofNumber, ghostRank,
+                                                                       m_maxAbsolutePresChange, pressureScalingFactor );
 
         scalingFactor = std::min( scalingFactor, subRegionData.first );
         minPresScalingFactor = std::min( minPresScalingFactor, subRegionData.first );
diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/SolutionScalingKernel.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/SolutionScalingKernel.hpp
index fe0572da818..2b61d39d579 100644
--- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/SolutionScalingKernel.hpp
+++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/SolutionScalingKernel.hpp
@@ -38,7 +38,8 @@ struct SolutionScalingKernel
                                              globalIndex const rankOffset,
                                              arrayView1d< globalIndex const > const & dofNumber,
                                              arrayView1d< integer const > const & ghostRank,
-                                             real64 const maxAbsolutePresChange )
+                                             real64 const maxAbsolutePresChange,
+                                             arrayView1d< real64 > pressureScalingFactor )
   {
     RAJA::ReduceMin< ReducePolicy< POLICY >, real64 > scalingFactor( 1.0 );
     RAJA::ReduceMax< ReducePolicy< POLICY >, real64 > maxDeltaPres( 0.0 );
@@ -47,6 +48,8 @@ struct SolutionScalingKernel
     {
       if( ghostRank[ei] < 0 && dofNumber[ei] >= 0 )
       {
+        pressureScalingFactor[ei] = 1.0;
+
         localIndex const lid = dofNumber[ei] - rankOffset;
 
         // compute the change in pressure
@@ -57,6 +60,7 @@ struct SolutionScalingKernel
         if( maxAbsolutePresChange > 0.0 && absPresChange > maxAbsolutePresChange )
         {
           real64 const presScalingFactor = maxAbsolutePresChange / absPresChange;
+          pressureScalingFactor[ei] = presScalingFactor;
           scalingFactor.min( presScalingFactor );
         }
       }
diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/SinglePhaseWell.cpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/SinglePhaseWell.cpp
index 15173c61ab2..db40d36c530 100644
--- a/src/coreComponents/physicsSolvers/fluidFlow/wells/SinglePhaseWell.cpp
+++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/SinglePhaseWell.cpp
@@ -964,9 +964,12 @@ SinglePhaseWell::scalingForLocalSystemSolution( WellElementSubRegion & subRegion
   arrayView1d< globalIndex const > const dofNumber = subRegion.getReference< array1d< globalIndex > >( wellDofKey );
   arrayView1d< integer const > const ghostRank = subRegion.ghostRank();
 
+  arrayView1d< real64 > const pressureScalingFactor = subRegion.getField< well::pressureScalingFactor >();
+
+  real64 scalingFactor = 1.0;
+
   if( isThermal() )
   {
-    arrayView1d< real64 > const pressureScalingFactor = subRegion.getField< well::pressureScalingFactor >();
     arrayView1d< real64 > const temperatureScalingFactor = subRegion.getField< well::temperatureScalingFactor >();
 
     integer const tempDofOffset = 2;
@@ -976,22 +979,25 @@ SinglePhaseWell::scalingForLocalSystemSolution( WellElementSubRegion & subRegion
                                                                    m_maxAbsolutePresChange, m_maxAbsoluteTempChange,
                                                                    pressureScalingFactor, temperatureScalingFactor );
 
+    scalingFactor = std::get< 0 >( subRegionData );
     maxDeltaPres = std::get< 1 >( subRegionData );
     maxDeltaTemp = std::get< 2 >( subRegionData );
     minPresScalingFactor = std::get< 3 >( subRegionData );
     minTempScalingFactor = std::get< 4 >( subRegionData );
-
-    return std::get< 0 >( subRegionData );
   }
+  else
+  {
+    auto const subRegionData = singlePhaseBaseKernels::
+                                 SolutionScalingKernel::
+                                 launch< parallelDevicePolicy<> >( localSolution, rankOffset, dofNumber, ghostRank,
+                                                                   m_maxAbsolutePresChange, pressureScalingFactor );
 
-  auto const subRegionData = singlePhaseBaseKernels::
-                               SolutionScalingKernel::
-                               launch< parallelDevicePolicy<> >( localSolution, rankOffset, dofNumber, ghostRank, m_maxAbsolutePresChange );
-
-  maxDeltaPres = subRegionData.second;
-  minPresScalingFactor = subRegionData.first;
+    scalingFactor = subRegionData.first;
+    maxDeltaPres = subRegionData.second;
+    minPresScalingFactor = subRegionData.first;
+  }
 
-  return subRegionData.first;
+  return scalingFactor;
 }
 
 real64

From 84a18285b1a2c340236187f1a9d989913a74a4f2 Mon Sep 17 00:00:00 2001
From: frankfeifan <frankfei0929@gmail.com>
Date: Sun, 11 Oct 2026 21:27:11 -0700
Subject: [PATCH 08/15] moved ScalingType to all flow solvers and enable local
 and global scaling for single phase FVM

---
 .../physicsSolvers/fluidFlow/CMakeLists.txt   |  1 +
 .../fluidFlow/CompositionalMultiphaseFVM.cpp  |  2 +-
 .../fluidFlow/CompositionalMultiphaseFVM.hpp  |  3 +-
 .../CompositionalMultiphaseHybridFVM.cpp      |  2 +-
 .../CompositionalMultiphaseUtilities.hpp      | 13 -----
 .../fluidFlow/FlowSolverUtilities.hpp         | 49 ++++++++++++++++
 .../fluidFlow/SinglePhaseBase.cpp             |  8 ++-
 .../fluidFlow/SinglePhaseBase.hpp             |  4 ++
 .../fluidFlow/SinglePhaseFVM.cpp              | 57 +++++++++++++++----
 .../fluidFlow/SinglePhaseFVM.hpp              | 10 +++-
 .../compositional/SolutionCheckKernel.hpp     |  9 +--
 .../ThermalSolutionCheckKernel.hpp            |  6 +-
 .../singlePhase/SolutionCheckKernel.hpp       |  7 ++-
 .../wells/CompositionalMultiphaseWell.cpp     |  2 +-
 14 files changed, 133 insertions(+), 40 deletions(-)
 create mode 100644 src/coreComponents/physicsSolvers/fluidFlow/FlowSolverUtilities.hpp

diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CMakeLists.txt b/src/coreComponents/physicsSolvers/fluidFlow/CMakeLists.txt
index 16e641b146e..31a13282532 100644
--- a/src/coreComponents/physicsSolvers/fluidFlow/CMakeLists.txt
+++ b/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
diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.cpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.cpp
index 235a129d64c..33ca6367136 100644
--- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.cpp
+++ b/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 )
diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.hpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.hpp
index 264b6fefec5..57c425a9718 100644
--- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.hpp
+++ b/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;
diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseHybridFVM.cpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseHybridFVM.cpp
index 6fd5cba3a5c..d525ac150ff 100644
--- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseHybridFVM.cpp
+++ b/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,
diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseUtilities.hpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseUtilities.hpp
index 536a52f5822..5b1286248e7 100644
--- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseUtilities.hpp
+++ b/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.
diff --git a/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverUtilities.hpp b/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverUtilities.hpp
new file mode 100644
index 00000000000..e66833b39b6
--- /dev/null
+++ b/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_
diff --git a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp
index c686dc677cc..f2f5165ddc9 100644
--- a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp
+++ b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp
@@ -62,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 ).
@@ -1421,6 +1422,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::
@@ -1429,7 +1431,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 );
diff --git a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.hpp b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.hpp
index 4d01e0a5370..c12e308f95f 100644
--- a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.hpp
+++ b/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;
diff --git a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseFVM.cpp b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseFVM.cpp
index 7ff30c34ea4..bd7457d63ea 100644
--- a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseFVM.cpp
+++ b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseFVM.cpp
@@ -57,12 +57,19 @@ using namespace constitutive;
 using namespace fields;
 using namespace singlePhaseBaseKernels;
 using namespace singlePhaseFVMKernels;
+using namespace flowSolverUtilities; // for ScalingType
 
 template< typename BASE >
 SinglePhaseFVM< BASE >::SinglePhaseFVM( const string & name,
                                         Group * const parent ):
   BASE( name, parent )
 {
+  this->registerWrapper( viewKeyStruct::scalingTypeString(), &m_scalingType ).
+    setInputFlag( dataRepository::InputFlags::OPTIONAL ).
+    setApplyDefaultValue( ScalingType::Global ).
+    setDescription( "Solution scaling type."
+                    "Valid options:\n* " + EnumStrings< ScalingType >::concat( "\n* " ) );
+
   LinearSolverParameters & linParams = m_linearSolverParameters.get();
   linParams.multiscale.fieldName = BASE::viewKeyStruct::elemDofFieldString();
   linParams.multiscale.label = "flow";
@@ -75,6 +82,15 @@ void SinglePhaseFVM< BASE >::initializePreSubGroups()
 
   this->checkDiscretizationName();
 
+  if( m_scalingType == ScalingType::Local &&
+      m_nonlinearSolverParameters.m_lineSearchAction != NonlinearSolverParameters::LineSearchAction::None )
+  {
+    GEOS_ERROR( GEOS_FMT( "line search is not supported for {} = {}",
+                          viewKeyStruct::scalingTypeString(),
+                          EnumStrings< ScalingType >::toString( ScalingType::Local )),
+                this->getDataContext());
+  }
+
   if( m_isThermal )
   {
     // For thermal simulations 2 pdes are considered so we let AMG know.
@@ -270,29 +286,48 @@ void SinglePhaseFVM< BASE >::applySystemSolution( DofManager const & dofManager,
                                                   DomainPartition & domain )
 {
   GEOS_UNUSED_VAR( dt );
-  if( m_isThermal )
-  {
-    DofManager::CompMask pressureMask( m_numDofPerCell, 0, 1 );
-    DofManager::CompMask temperatureMask( m_numDofPerCell, 1, 2 );
 
+  bool const localScaling = m_scalingType == ScalingType::Local;
+
+  DofManager::CompMask pressureMask( m_numDofPerCell, 0, 1 );
+
+  if( localScaling )
+  {
     dofManager.addVectorToField( localSolution,
                                  BASE::viewKeyStruct::elemDofFieldString(),
                                  flow::pressure::key(),
                                  flow::pressureScalingFactor::key(),
                                  pressureMask );
-
-    dofManager.addVectorToField( localSolution,
-                                 BASE::viewKeyStruct::elemDofFieldString(),
-                                 flow::temperature::key(),
-                                 flow::temperatureScalingFactor::key(),
-                                 temperatureMask );
   }
   else
   {
     dofManager.addVectorToField( localSolution,
                                  BASE::viewKeyStruct::elemDofFieldString(),
                                  flow::pressure::key(),
-                                 scalingFactor );
+                                 scalingFactor,
+                                 pressureMask );
+  }
+
+  if( m_isThermal )
+  {
+    DofManager::CompMask temperatureMask( m_numDofPerCell, 1, 2 );
+
+    if( localScaling )
+    {
+      dofManager.addVectorToField( localSolution,
+                                   BASE::viewKeyStruct::elemDofFieldString(),
+                                   flow::temperature::key(),
+                                   flow::temperatureScalingFactor::key(),
+                                   temperatureMask );
+    }
+    else
+    {
+      dofManager.addVectorToField( localSolution,
+                                   BASE::viewKeyStruct::elemDofFieldString(),
+                                   flow::temperature::key(),
+                                   scalingFactor,
+                                   temperatureMask );
+    }
   }
 
   this->forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &,
diff --git a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseFVM.hpp b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseFVM.hpp
index a0adca9dec5..9a32607b84e 100644
--- a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseFVM.hpp
+++ b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseFVM.hpp
@@ -20,6 +20,7 @@
 #ifndef GEOS_PHYSICSSOLVERS_FLUIDFLOW_SINGLEPHASEFVM_HPP_
 #define GEOS_PHYSICSSOLVERS_FLUIDFLOW_SINGLEPHASEFVM_HPP_
 
+#include "physicsSolvers/fluidFlow/FlowSolverUtilities.hpp"
 #include "physicsSolvers/fluidFlow/SinglePhaseBase.hpp"
 #include "physicsSolvers/fluidFlow/SinglePhaseProppantBase.hpp"
 
@@ -51,6 +52,7 @@ class SinglePhaseFVM : public BASE
   // have to use this->member etc.
   using BASE::m_numDofPerCell;
   using BASE::m_isThermal;
+  using BASE::m_scalingType;
 
   /**
    * @brief main constructor for Group Objects
@@ -193,6 +195,12 @@ class SinglePhaseFVM : public BASE
 
   virtual void initializePreSubGroups() override;
 
+  struct viewKeyStruct : BASE::viewKeyStruct
+  {
+    // nonlinear solver parameters
+    static constexpr char const * scalingTypeString() { return "scalingType"; }
+  };
+
 private:
 
   /**
@@ -211,8 +219,6 @@ class SinglePhaseFVM : public BASE
                              CRSMatrixView< real64, globalIndex const > const & localMatrix,
                              arrayView1d< real64 > const & localRhs );
 
-  // no data needed here, see SinglePhaseBase
-
 };
 
 } /* namespace geos */
diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/SolutionCheckKernel.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/SolutionCheckKernel.hpp
index 639decdb418..6ba2f9d9fcb 100644
--- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/SolutionCheckKernel.hpp
+++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/SolutionCheckKernel.hpp
@@ -21,6 +21,7 @@
 #define GEOS_PHYSICSSOLVERS_FLUIDFLOW_COMPOSITIONAL_SOLUTIONCHECKKERNEL_HPP
 
 #include "LvArray/src/math.hpp"
+#include "physicsSolvers/fluidFlow/FlowSolverUtilities.hpp"
 #include "physicsSolvers/fluidFlow/kernels/compositional/SolutionScalingAndCheckingKernelBase.hpp"
 #include "physicsSolvers/fluidFlow/kernels/SolutionCheckKernelsHelpers.hpp"
 
@@ -63,7 +64,7 @@ class SolutionCheckKernel : public SolutionScalingAndCheckingKernelBase< integer
    */
   SolutionCheckKernel( integer const allowCompDensChopping,
                        integer const allowNegativePressure,
-                       compositionalMultiphaseUtilities::ScalingType const scalingType,
+                       flowSolverUtilities::ScalingType const scalingType,
                        real64 const scalingFactor,
                        arrayView1d< real64 const > const pressure,
                        arrayView2d< real64 const, compflow::USD_COMP > const compDens,
@@ -232,7 +233,7 @@ class SolutionCheckKernel : public SolutionScalingAndCheckingKernelBase< integer
                              StackVariables & stack,
                              FUNC && kernelOp = NoOpFunc{} ) const
   {
-    bool const localScaling = m_scalingType == compositionalMultiphaseUtilities::ScalingType::Local;
+    bool const localScaling = m_scalingType == flowSolverUtilities::ScalingType::Local;
 
     real64 const newPres = m_pressure[ei] + (localScaling ? m_pressureScalingFactor[ei] : m_scalingFactor) * m_localSolution[stack.localRow];
     if( newPres < 0.0 && !m_allowNegativePressure )
@@ -306,7 +307,7 @@ class SolutionCheckKernel : public SolutionScalingAndCheckingKernelBase< integer
   real64 const m_scalingFactor;
 
   /// scaling type (global or local)
-  compositionalMultiphaseUtilities::ScalingType const m_scalingType;
+  flowSolverUtilities::ScalingType const m_scalingType;
 
   ElementsReporterCollector const m_negPressureIds;
 
@@ -338,7 +339,7 @@ class SolutionCheckKernelFactory
   static SolutionCheckKernel::KernelStats
   createAndLaunch( integer const allowCompDensChopping,
                    integer const allowNegativePressure,
-                   compositionalMultiphaseUtilities::ScalingType const scalingType,
+                   flowSolverUtilities::ScalingType const scalingType,
                    real64 const scalingFactor,
                    arrayView1d< real64 const > const pressure,
                    arrayView2d< real64 const, compflow::USD_COMP > const compDens,
diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/ThermalSolutionCheckKernel.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/ThermalSolutionCheckKernel.hpp
index 10566b06827..04f7b2e37b8 100644
--- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/ThermalSolutionCheckKernel.hpp
+++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/ThermalSolutionCheckKernel.hpp
@@ -60,7 +60,7 @@ class SolutionCheckKernel : public isothermalCompositionalMultiphaseBaseKernels:
    */
   SolutionCheckKernel( integer const allowCompDensChopping,
                        integer const allowNegativePressure,
-                       compositionalMultiphaseUtilities::ScalingType const scalingType,
+                       flowSolverUtilities::ScalingType const scalingType,
                        real64 const scalingFactor,
                        arrayView1d< real64 const > const pressure,
                        arrayView1d< real64 const > const temperature,
@@ -109,7 +109,7 @@ class SolutionCheckKernel : public isothermalCompositionalMultiphaseBaseKernels:
   {
     Base::computeSolutionCheck( ei, stack, [&] ()
     {
-      bool const localScaling = m_scalingType == compositionalMultiphaseUtilities::ScalingType::Local;
+      bool const localScaling = m_scalingType == flowSolverUtilities::ScalingType::Local;
       // compute the change in temperature
       real64 const newTemp = m_temperature[ei] + (localScaling ? m_temperatureScalingFactor[ei] : m_scalingFactor * m_localSolution[stack.localRow + m_temperatureOffset]);
       if( newTemp < minTemperature )
@@ -155,7 +155,7 @@ class SolutionCheckKernelFactory
   static SolutionCheckKernel::KernelStats
   createAndLaunch( integer const allowCompDensChopping,
                    integer const allowNegativePressure,
-                   compositionalMultiphaseUtilities::ScalingType const scalingType,
+                   flowSolverUtilities::ScalingType const scalingType,
                    real64 const scalingFactor,
                    arrayView1d< real64 const > const pressure,
                    arrayView1d< real64 const > const temperature,
diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/SolutionCheckKernel.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/SolutionCheckKernel.hpp
index 16bd95ea27d..6892f7adc7f 100644
--- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/SolutionCheckKernel.hpp
+++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/SolutionCheckKernel.hpp
@@ -22,6 +22,7 @@
 
 #include "common/DataTypes.hpp"
 #include "common/GEOS_RAJA_Interface.hpp"
+#include "physicsSolvers/fluidFlow/FlowSolverUtilities.hpp"
 #include "physicsSolvers/fluidFlow/kernels/SolutionCheckKernelsHelpers.hpp"
 
 namespace geos
@@ -46,19 +47,23 @@ struct SolutionCheckKernel
                              arrayView1d< globalIndex const > const & dofNumber,
                              arrayView1d< integer const > const & ghostRank,
                              arrayView1d< real64 const > const & pres,
+                             flowSolverUtilities::ScalingType const scalingType,
                              real64 const scalingFactor,
+                             arrayView1d< real64 const > const & perElementScaling,
                              ElementsReporterCollector const & negPressureIds )
   {
     using reducePolicy = ReducePolicy< POLICY >;
     using atomicPolicy = AtomicPolicy< POLICY >;
     RAJA::ReduceMin< reducePolicy, real64 > minNegPres( 0.0 );
 
+    bool const localScaling = scalingType == flowSolverUtilities::ScalingType::Local;
+
     forAll< POLICY >( dofNumber.size(), [=] GEOS_HOST_DEVICE ( localIndex const ei )
     {
       if( ghostRank[ei] < 0 && dofNumber[ei] >= 0 )
       {
         localIndex const lid = dofNumber[ei] - rankOffset;
-        real64 const newPres = pres[ei] + scalingFactor * localSolution[lid];
+        real64 const newPres = pres[ei] + ( localScaling ? perElementScaling[ei] : scalingFactor ) * localSolution[lid];
 
         if( newPres <= 0.0 )
         {
diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.cpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.cpp
index b5edc3d7643..28fc71ba42d 100644
--- a/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.cpp
+++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.cpp
@@ -1757,7 +1757,7 @@ CompositionalMultiphaseWell::checkWellSystemSolution( WellElementSubRegion & sub
   integer localCheck = 1;
 
   integer constexpr allowNegativePressure = 0;
-  compositionalMultiphaseUtilities::ScalingType constexpr scalingType = compositionalMultiphaseUtilities::ScalingType::Global;
+  flowSolverUtilities::ScalingType constexpr scalingType = flowSolverUtilities::ScalingType::Global;
 
   arrayView1d< real64 const > const pressure = subRegion.getField< well::pressure >();
   arrayView1d< real64 const > const temperature = subRegion.getField< well::temperature >();

From 46a8978467851c37337523faabac22489e90e217 Mon Sep 17 00:00:00 2001
From: frankfeifan <frankfei0929@gmail.com>
Date: Mon, 12 Oct 2026 10:23:21 -0700
Subject: [PATCH 09/15] applied global scaling as always for now for
 singlePhaseWell

---
 .../physicsSolvers/fluidFlow/wells/SinglePhaseWell.cpp    | 8 ++++++--
 1 file changed, 6 insertions(+), 2 deletions(-)

diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/SinglePhaseWell.cpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/SinglePhaseWell.cpp
index db40d36c530..200e6f68459 100644
--- a/src/coreComponents/physicsSolvers/fluidFlow/wells/SinglePhaseWell.cpp
+++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/SinglePhaseWell.cpp
@@ -1181,16 +1181,20 @@ bool SinglePhaseWell::checkWellSystemSolution( WellElementSubRegion & subRegion,
 
   // get a reference to the primary variables on well elements
   arrayView1d< real64 const > const & pressure = subRegion.getField< well::pressure >();
+  arrayView1d< real64 const > const & pressureScalingFactor = subRegion.getField< well::pressureScalingFactor >();
   auto const negPresCollector = negPressureIds.createCollector( subRegion.localToGlobalMap().toViewConst() );
 
   using Kernel = singlePhaseBaseKernels::SolutionCheckKernel;
 
+  // the well solution is always scaled with the global factor
   auto const results = Kernel::launch< parallelDevicePolicy<> >( localSolution,
                                                                  rankOffset,
                                                                  dofNumber,
                                                                  ghostRank,
                                                                  pressure,
+                                                                 flowSolverUtilities::ScalingType::Global,
                                                                  scalingFactor,
+                                                                 pressureScalingFactor,
                                                                  negPresCollector );
 
   minPressure = std::min( minPressure, results.minNegPres );
@@ -1214,7 +1218,7 @@ SinglePhaseWell::applyWellSystemSolution( DofManager const & dofManager,
   dofManager.addVectorToField( localSolution,
                                wellElementDofName(),
                                well::pressure::key(),
-                               well::pressureScalingFactor::key(),
+                               scalingFactor,
                                pressureMask );
 
   dofManager.addVectorToField( localSolution,
@@ -1230,7 +1234,7 @@ SinglePhaseWell::applyWellSystemSolution( DofManager const & dofManager,
     dofManager.addVectorToField( localSolution,
                                  wellElementDofName(),
                                  fields::well::temperature::key(),
-                                 well::temperatureScalingFactor::key(),
+                                 scalingFactor,
                                  temperatureMask );
 
   }

From f271abc45a69a1604d941f658b145d3bd47bf7e3 Mon Sep 17 00:00:00 2001
From: frankfeifan <frankfei0929@gmail.com>
Date: Mon, 12 Oct 2026 10:26:22 -0700
Subject: [PATCH 10/15] schema update

---
 src/coreComponents/schema/schema.xsd | 16 ++++++++++++++--
 1 file changed, 14 insertions(+), 2 deletions(-)

diff --git a/src/coreComponents/schema/schema.xsd b/src/coreComponents/schema/schema.xsd
index 1ad66a96106..7ee089d492e 100644
--- a/src/coreComponents/schema/schema.xsd
+++ b/src/coreComponents/schema/schema.xsd
@@ -3328,7 +3328,7 @@ Information output from lower logLevels is added with the desired log level
 		<!--scalingType => Solution scaling type.Valid options:
 * Global
 * Local-->
-		<xsd:attribute name="scalingType" type="geos_compositionalMultiphaseUtilities_ScalingType" default="Global" />
+		<xsd:attribute name="scalingType" type="geos_flowSolverUtilities_ScalingType" default="Global" />
 		<!--solutionChangeScalingFactor => Damping factor for solution change targets-->
 		<xsd:attribute name="solutionChangeScalingFactor" type="real64" default="0.5" />
 		<!--targetCompFracChangeInTimeStep => Target change in component fraction in a time step-->
@@ -3381,7 +3381,7 @@ When set to `all` output both convergence & iteration information to a csv.-->
 			<xsd:pattern value=".*[\[\]`$].*|ArithmeticAverage|PhasePresence" />
 		</xsd:restriction>
 	</xsd:simpleType>
-	<xsd:simpleType name="geos_compositionalMultiphaseUtilities_ScalingType">
+	<xsd:simpleType name="geos_flowSolverUtilities_ScalingType">
 		<xsd:restriction base="xsd:string">
 			<xsd:pattern value=".*[\[\]`$].*|Global|Local" />
 		</xsd:restriction>
@@ -4904,6 +4904,10 @@ Information output from lower logLevels is added with the desired log level
 		<xsd:attribute name="maxSequentialPressureChange" type="real64" default="100000" />
 		<!--maxSequentialTemperatureChange => Maximum (absolute) temperature change in a sequential iteration, used for outer loop convergence check-->
 		<xsd:attribute name="maxSequentialTemperatureChange" type="real64" default="0.1" />
+		<!--scalingType => Solution scaling type.Valid options:
+* Global
+* Local-->
+		<xsd:attribute name="scalingType" type="geos_flowSolverUtilities_ScalingType" default="Global" />
 		<!--targetRegions => Allowable regions that the solver may be applied to. Note that this does not indicate that the solver will be applied to these regions, only that allocation will occur such that the solver may be applied to these regions. The decision about what regions this solver will beapplied to rests in the EventManager.-->
 		<xsd:attribute name="targetRegions" type="groupNameRef_array" use="required" />
 		<!--temperature => Temperature-->
@@ -5316,6 +5320,10 @@ Information output from lower logLevels is added with the desired log level
 		<xsd:attribute name="maxSequentialPressureChange" type="real64" default="100000" />
 		<!--maxSequentialTemperatureChange => Maximum (absolute) temperature change in a sequential iteration, used for outer loop convergence check-->
 		<xsd:attribute name="maxSequentialTemperatureChange" type="real64" default="0.1" />
+		<!--scalingType => Solution scaling type.Valid options:
+* Global
+* Local-->
+		<xsd:attribute name="scalingType" type="geos_flowSolverUtilities_ScalingType" default="Global" />
 		<!--targetRegions => Allowable regions that the solver may be applied to. Note that this does not indicate that the solver will be applied to these regions, only that allocation will occur such that the solver may be applied to these regions. The decision about what regions this solver will beapplied to rests in the EventManager.-->
 		<xsd:attribute name="targetRegions" type="groupNameRef_array" use="required" />
 		<!--temperature => Temperature-->
@@ -6182,6 +6190,10 @@ Information output from lower logLevels is added with the desired log level
 2
  - The summary of declared fields and coupling-->
 		<xsd:attribute name="logLevel" type="integer" default="0" />
+		<!--scalingType => Solution scaling type (Local is only supported for single-phase wells).Valid options:
+* Global
+* Local-->
+		<xsd:attribute name="scalingType" type="geos_flowSolverUtilities_ScalingType" default="Global" />
 		<!--targetRegions => Allowable regions that the solver may be applied to. Note that this does not indicate that the solver will be applied to these regions, only that allocation will occur such that the solver may be applied to these regions. The decision about what regions this solver will beapplied to rests in the EventManager.-->
 		<xsd:attribute name="targetRegions" type="groupNameRef_array" use="required" />
 		<!--timeStepFromTables => Choose time step to honor rates/bhp tables time intervals-->

From afb61102e9de31bb0155aeb5aa357ba5716f14d0 Mon Sep 17 00:00:00 2001
From: frankfeifan <frankfei0929@gmail.com>
Date: Wed, 4 Nov 2026 14:26:31 -0800
Subject: [PATCH 11/15] refactor the output type for solution scaling kernel

---
 .../fluidFlow/SinglePhaseBase.cpp             | 16 ++++-----
 .../singlePhase/SolutionScalingKernel.hpp     | 36 +++++++++++++++----
 .../ThermalSolutionScalingKernel.hpp          | 28 +++++++++------
 .../fluidFlow/wells/SinglePhaseWell.cpp       | 16 ++++-----
 4 files changed, 63 insertions(+), 33 deletions(-)

diff --git a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp
index f2f5165ddc9..6c9a98617e8 100644
--- a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp
+++ b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp
@@ -1336,12 +1336,12 @@ real64 SinglePhaseBase::scalingForSystemSolution( DomainPartition & domain,
                                                                        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 ) );
+        scalingFactor = std::min( scalingFactor, subRegionData.scalingFactor );
+        maxDeltaPres  = std::max( maxDeltaPres, subRegionData.maxDeltaPres );
+        maxDeltaTemp  = std::max( maxDeltaTemp, subRegionData.maxDeltaTemp );
 
-        minPresScalingFactor = std::min( minPresScalingFactor, std::get< 3 >( subRegionData ) );
-        minTempScalingFactor = std::min( minTempScalingFactor, std::get< 4 >( subRegionData ) );
+        minPresScalingFactor = std::min( minPresScalingFactor, subRegionData.minPresScalingFactor );
+        minTempScalingFactor = std::min( minTempScalingFactor, subRegionData.minTempScalingFactor );
       } );
     } );
   }
@@ -1368,9 +1368,9 @@ real64 SinglePhaseBase::scalingForSystemSolution( DomainPartition & domain,
                                      launch< parallelDevicePolicy<> >( localSolution, rankOffset, dofNumber, ghostRank,
                                                                        m_maxAbsolutePresChange, pressureScalingFactor );
 
-        scalingFactor = std::min( scalingFactor, subRegionData.first );
-        minPresScalingFactor = std::min( minPresScalingFactor, subRegionData.first );
-        maxDeltaPres  = std::max( maxDeltaPres, subRegionData.second );
+        scalingFactor = std::min( scalingFactor, subRegionData.scalingFactor );
+        minPresScalingFactor = std::min( minPresScalingFactor, subRegionData.minPresScalingFactor );
+        maxDeltaPres  = std::max( maxDeltaPres, subRegionData.maxDeltaPres );
       } );
     } );
   }
diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/SolutionScalingKernel.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/SolutionScalingKernel.hpp
index 2b61d39d579..5da03150d63 100644
--- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/SolutionScalingKernel.hpp
+++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/SolutionScalingKernel.hpp
@@ -31,15 +31,33 @@ namespace singlePhaseBaseKernels
 
 /******************************** SolutionScalingKernel ********************************/
 
+/**
+ * @struct SolutionScalingData
+ * @brief Rank-local reductions from the single-phase solution scaling kernels
+ */
+struct SolutionScalingData
+{
+  /// the Newton update scaling factor
+  real64 scalingFactor = 1.0;
+  /// the maximum absolute pressure change (before scaling)
+  real64 maxDeltaPres = 0.0;
+  /// the maximum absolute temperature change (before scaling)
+  real64 maxDeltaTemp = 0.0;
+  /// the minimum pressure scaling factor
+  real64 minPresScalingFactor = 1.0;
+  /// the minimum temperature scaling factor
+  real64 minTempScalingFactor = 1.0;
+};
+
 struct SolutionScalingKernel
 {
   template< typename POLICY >
-  static std::pair< real64, real64 > launch( arrayView1d< real64 const > const & localSolution,
-                                             globalIndex const rankOffset,
-                                             arrayView1d< globalIndex const > const & dofNumber,
-                                             arrayView1d< integer const > const & ghostRank,
-                                             real64 const maxAbsolutePresChange,
-                                             arrayView1d< real64 > pressureScalingFactor )
+  static SolutionScalingData launch( arrayView1d< real64 const > const & localSolution,
+                                     globalIndex const rankOffset,
+                                     arrayView1d< globalIndex const > const & dofNumber,
+                                     arrayView1d< integer const > const & ghostRank,
+                                     real64 const maxAbsolutePresChange,
+                                     arrayView1d< real64 > pressureScalingFactor )
   {
     RAJA::ReduceMin< ReducePolicy< POLICY >, real64 > scalingFactor( 1.0 );
     RAJA::ReduceMax< ReducePolicy< POLICY >, real64 > maxDeltaPres( 0.0 );
@@ -67,7 +85,11 @@ struct SolutionScalingKernel
 
     } );
 
-    return { scalingFactor.get(), maxDeltaPres.get() };
+    SolutionScalingData scalingData;
+    scalingData.scalingFactor = scalingFactor.get();
+    scalingData.maxDeltaPres = maxDeltaPres.get();
+    scalingData.minPresScalingFactor = scalingFactor.get();
+    return scalingData;
   }
 
 };
diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/ThermalSolutionScalingKernel.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/ThermalSolutionScalingKernel.hpp
index 5ffc3ff7ed7..5b6a08be140 100644
--- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/ThermalSolutionScalingKernel.hpp
+++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/ThermalSolutionScalingKernel.hpp
@@ -22,6 +22,7 @@
 
 #include "common/DataTypes.hpp"
 #include "common/GEOS_RAJA_Interface.hpp"
+#include "physicsSolvers/fluidFlow/kernels/singlePhase/SolutionScalingKernel.hpp"
 
 namespace geos
 {
@@ -34,15 +35,16 @@ namespace thermalSinglePhaseBaseKernels
 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 )
+  static singlePhaseBaseKernels::SolutionScalingData
+  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 );
@@ -91,7 +93,13 @@ struct SolutionScalingKernel
 
     } );
 
-    return { scalingFactor.get(), maxDeltaPres.get(), maxDeltaTemp.get(), localMinPresScalingFactor.get(), localMinTempScalingFactor.get() };
+    singlePhaseBaseKernels::SolutionScalingData scalingData;
+    scalingData.scalingFactor = scalingFactor.get();
+    scalingData.maxDeltaPres = maxDeltaPres.get();
+    scalingData.maxDeltaTemp = maxDeltaTemp.get();
+    scalingData.minPresScalingFactor = localMinPresScalingFactor.get();
+    scalingData.minTempScalingFactor = localMinTempScalingFactor.get();
+    return scalingData;
   }
 
 };
diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/SinglePhaseWell.cpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/SinglePhaseWell.cpp
index 200e6f68459..111d1cd43a0 100644
--- a/src/coreComponents/physicsSolvers/fluidFlow/wells/SinglePhaseWell.cpp
+++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/SinglePhaseWell.cpp
@@ -979,11 +979,11 @@ SinglePhaseWell::scalingForLocalSystemSolution( WellElementSubRegion & subRegion
                                                                    m_maxAbsolutePresChange, m_maxAbsoluteTempChange,
                                                                    pressureScalingFactor, temperatureScalingFactor );
 
-    scalingFactor = std::get< 0 >( subRegionData );
-    maxDeltaPres = std::get< 1 >( subRegionData );
-    maxDeltaTemp = std::get< 2 >( subRegionData );
-    minPresScalingFactor = std::get< 3 >( subRegionData );
-    minTempScalingFactor = std::get< 4 >( subRegionData );
+    scalingFactor = subRegionData.scalingFactor;
+    maxDeltaPres = subRegionData.maxDeltaPres;
+    maxDeltaTemp = subRegionData.maxDeltaTemp;
+    minPresScalingFactor = subRegionData.minPresScalingFactor;
+    minTempScalingFactor = subRegionData.minTempScalingFactor;
   }
   else
   {
@@ -992,9 +992,9 @@ SinglePhaseWell::scalingForLocalSystemSolution( WellElementSubRegion & subRegion
                                  launch< parallelDevicePolicy<> >( localSolution, rankOffset, dofNumber, ghostRank,
                                                                    m_maxAbsolutePresChange, pressureScalingFactor );
 
-    scalingFactor = subRegionData.first;
-    maxDeltaPres = subRegionData.second;
-    minPresScalingFactor = subRegionData.first;
+    scalingFactor = subRegionData.scalingFactor;
+    maxDeltaPres = subRegionData.maxDeltaPres;
+    minPresScalingFactor = subRegionData.minPresScalingFactor;
   }
 
   return scalingFactor;

From 0e86face73c799bc73a0b2268c49df3ac91ebf0c Mon Sep 17 00:00:00 2001
From: frankfeifan <frankfei0929@gmail.com>
Date: Wed, 4 Nov 2026 19:38:33 -0800
Subject: [PATCH 12/15] add maxAbsoluteTemperatureChange in compositional
 solver

---
 .../fluidFlow/CompositionalMultiphaseFVM.cpp  |  1 +
 .../fluidFlow/SinglePhaseBase.cpp             | 10 ++++-
 .../ThermalSolutionScalingKernel.hpp          | 42 ++++++++++++-------
 .../wells/CompositionalMultiphaseWell.cpp     |  1 +
 4 files changed, 38 insertions(+), 16 deletions(-)

diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.cpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.cpp
index 33ca6367136..59bbf02873d 100644
--- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.cpp
+++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.cpp
@@ -651,6 +651,7 @@ real64 CompositionalMultiphaseFVM::scalingForSystemSolution( DomainPartition & d
                                      createAndLaunch< parallelDevicePolicy<> >( m_maxRelativePresChange,
                                                                                 m_maxAbsolutePresChange,
                                                                                 m_maxRelativeTempChange,
+                                                                                m_maxAbsoluteTempChange,
                                                                                 m_maxCompFracChange,
                                                                                 m_maxRelativeCompDensChange,
                                                                                 pressure,
diff --git a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp
index 6c9a98617e8..0e58c465d1f 100644
--- a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp
+++ b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp
@@ -1336,7 +1336,10 @@ real64 SinglePhaseBase::scalingForSystemSolution( DomainPartition & domain,
                                                                        m_maxAbsolutePresChange, m_maxAbsoluteTempChange,
                                                                        pressureScalingFactor, temperatureScalingFactor );
 
-        scalingFactor = std::min( scalingFactor, subRegionData.scalingFactor );
+        if( m_scalingType == flowSolverUtilities::ScalingType::Global )
+        {
+          scalingFactor = std::min( scalingFactor, subRegionData.scalingFactor );
+        }
         maxDeltaPres  = std::max( maxDeltaPres, subRegionData.maxDeltaPres );
         maxDeltaTemp  = std::max( maxDeltaTemp, subRegionData.maxDeltaTemp );
 
@@ -1368,7 +1371,10 @@ real64 SinglePhaseBase::scalingForSystemSolution( DomainPartition & domain,
                                      launch< parallelDevicePolicy<> >( localSolution, rankOffset, dofNumber, ghostRank,
                                                                        m_maxAbsolutePresChange, pressureScalingFactor );
 
-        scalingFactor = std::min( scalingFactor, subRegionData.scalingFactor );
+        if( m_scalingType == flowSolverUtilities::ScalingType::Global )
+        {
+          scalingFactor = std::min( scalingFactor, subRegionData.scalingFactor );
+        }
         minPresScalingFactor = std::min( minPresScalingFactor, subRegionData.minPresScalingFactor );
         maxDeltaPres  = std::max( maxDeltaPres, subRegionData.maxDeltaPres );
       } );
diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/ThermalSolutionScalingKernel.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/ThermalSolutionScalingKernel.hpp
index 1a4a2b4c782..5f29291ddf1 100644
--- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/ThermalSolutionScalingKernel.hpp
+++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/ThermalSolutionScalingKernel.hpp
@@ -47,6 +47,7 @@ class SolutionScalingKernel : public isothermalCompositionalMultiphaseBaseKernel
    * @param[in] maxRelativePresChange the max allowed relative pressure change
    * @param[in] maxAbsolutePresChange the max allowed absolute pressure change
    * @param[in] maxRelativeTempChange the max allowed relative temperature change
+   * @param[in] maxAbsoluteTempChange the max allowed absolute temperature change
    * @param[in] maxCompFracChange the max allowed comp fraction change
    * @param[in] maxRelativeCompDensChange the max allowed comp density change
    * @param[in] rankOffset the rank offset
@@ -64,6 +65,7 @@ class SolutionScalingKernel : public isothermalCompositionalMultiphaseBaseKernel
   SolutionScalingKernel( real64 const maxRelativePresChange,
                          real64 const maxAbsolutePresChange,
                          real64 const maxRelativeTempChange,
+                         real64 const maxAbsoluteTempChange,
                          real64 const maxCompFracChange,
                          real64 const maxRelativeCompDensChange,
                          globalIndex const rankOffset,
@@ -92,6 +94,7 @@ class SolutionScalingKernel : public isothermalCompositionalMultiphaseBaseKernel
             pressureScalingFactor,
             compDensScalingFactor ),
     m_maxRelativeTempChange( maxRelativeTempChange ),
+    m_maxAbsoluteTempChange( maxAbsoluteTempChange ),
     m_temperature( temperature ),
     m_temperatureScalingFactor( temperatureScalingFactor ),
     m_temperatureOffset( temperatureOffset )
@@ -129,25 +132,33 @@ class SolutionScalingKernel : public isothermalCompositionalMultiphaseBaseKernel
         stack.localMaxDeltaTemp = absTempChange;
       }
 
-      m_temperatureScalingFactor[ei] = 1.0;
-
-      if( temp > eps )
+      // compute temperature scaling factor
+      real64 tempScalingFactor = 1.0;
+      // when enabled, absolute change scaling has a priority over relative change
+      if( m_maxAbsoluteTempChange > 0.0 ) // maxAbsoluteTempChange <= 0.0 means that absolute scaling is disabled
+      {
+        if( absTempChange > m_maxAbsoluteTempChange )
+        {
+          tempScalingFactor = m_maxAbsoluteTempChange / absTempChange;
+        }
+      }
+      else if( temp > eps )
       {
         real64 const relativeTempChange = absTempChange / temp;
         if( relativeTempChange > m_maxRelativeTempChange )
         {
-          real64 const tempScalingFactor = m_maxRelativeTempChange / relativeTempChange;
-          m_temperatureScalingFactor[ei] = tempScalingFactor;
-          if( stack.localMinVal > tempScalingFactor )
-          {
-            stack.localMinVal = tempScalingFactor;
-          }
-          if( stack.localMinTempScalingFactor > tempScalingFactor )
-          {
-            stack.localMinTempScalingFactor = tempScalingFactor;
-          }
+          tempScalingFactor = m_maxRelativeTempChange / relativeTempChange;
         }
       }
+      m_temperatureScalingFactor[ei] = tempScalingFactor;
+      if( stack.localMinVal > tempScalingFactor )
+      {
+        stack.localMinVal = tempScalingFactor;
+      }
+      if( stack.localMinTempScalingFactor > tempScalingFactor )
+      {
+        stack.localMinTempScalingFactor = tempScalingFactor;
+      }
     } );
   }
 
@@ -155,6 +166,7 @@ class SolutionScalingKernel : public isothermalCompositionalMultiphaseBaseKernel
 
   /// Max allowed changes in primary variables
   real64 const m_maxRelativeTempChange;
+  real64 const m_maxAbsoluteTempChange;
 
   /// View on the primary variables
   arrayView1d< real64 const > const m_temperature;
@@ -179,6 +191,7 @@ class SolutionScalingKernelFactory
    * @param[in] maxRelativePresChange the max allowed relative pressure change
    * @param[in] maxAbsolutePresChange the max allowed absolute pressure change
    * @param[in] maxRelativeTempChange the max allowed relative temperature change
+   * @param[in] maxAbsoluteTempChange the max allowed absolute temperature change
    * @param[in] maxCompFracChange the max allowed comp fraction change
    * @param[in] maxRelativeCompdensChange the max allowed relative component density change
    * @param[in] rankOffset the rank offset
@@ -192,6 +205,7 @@ class SolutionScalingKernelFactory
   createAndLaunch( real64 const maxRelativePresChange,
                    real64 const maxAbsolutePresChange,
                    real64 const maxRelativeTempChange,
+                   real64 const maxAbsoluteTempChange,
                    real64 const maxCompFracChange,
                    real64 const maxRelativeCompDensChange,
                    arrayView1d< real64 const > const pressure,
@@ -207,7 +221,7 @@ class SolutionScalingKernelFactory
                    arrayView1d< real64 const > const localSolution,
                    integer const temperatureOffset )
   {
-    SolutionScalingKernel kernel( maxRelativePresChange, maxAbsolutePresChange, maxRelativeTempChange,
+    SolutionScalingKernel kernel( maxRelativePresChange, maxAbsolutePresChange, maxRelativeTempChange, maxAbsoluteTempChange,
                                   maxCompFracChange, maxRelativeCompDensChange,
                                   rankOffset, numComp, dofKey, subRegion, localSolution,
                                   pressure, temperature, compDens, pressureScalingFactor,
diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.cpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.cpp
index d35f68979b6..ce8cd45c5f7 100644
--- a/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.cpp
+++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.cpp
@@ -1608,6 +1608,7 @@ CompositionalMultiphaseWell::scalingForLocalSystemSolution( WellElementSubRegion
       createAndLaunch< parallelDevicePolicy<> >( m_maxRelativePresChange,
                                                  m_maxAbsolutePresChange,
                                                  m_maxRelativeTempChange,
+                                                 m_maxAbsoluteTempChange,
                                                  m_maxCompFracChange,
                                                  m_maxRelativeCompDensChange,
                                                  pressure,

From 1c4422ed247f4fec722331dcf4aa7a1bb38b89b9 Mon Sep 17 00:00:00 2001
From: frankfeifan <frankfei0929@gmail.com>
Date: Thu, 5 Nov 2026 10:36:45 -0800
Subject: [PATCH 13/15] ensure Local scaling compatible with lineSearchAction
 in the coupled solver

---
 .../fluidFlow/CompositionalMultiphaseFVM.cpp              | 8 ++++++++
 .../physicsSolvers/fluidFlow/SinglePhaseBase.cpp          | 8 ++++++++
 2 files changed, 16 insertions(+)

diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.cpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.cpp
index 59bbf02873d..b6191b9942f 100644
--- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.cpp
+++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.cpp
@@ -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 );
diff --git a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp
index 0e58c465d1f..b620a1eb3f5 100644
--- a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp
+++ b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp
@@ -1305,6 +1305,14 @@ 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;

From f326e8401cc1bdf2b6d024850deae3516da56e38 Mon Sep 17 00:00:00 2001
From: frankfeifan <frankfei0929@gmail.com>
Date: Thu, 5 Nov 2026 10:42:05 -0800
Subject: [PATCH 14/15] pick the stricter one from abs and rel limits

---
 .../kernels/compositional/SolutionScalingKernel.hpp         | 6 +++---
 .../kernels/compositional/ThermalSolutionScalingKernel.hpp  | 6 +++---
 .../zFormulation/SolutionScalingZFormulationKernel.hpp      | 6 +++---
 src/coreComponents/schema/schema.xsd                        | 4 ----
 4 files changed, 9 insertions(+), 13 deletions(-)

diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/SolutionScalingKernel.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/SolutionScalingKernel.hpp
index 316859e9de7..da29b394a73 100644
--- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/SolutionScalingKernel.hpp
+++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/SolutionScalingKernel.hpp
@@ -252,7 +252,7 @@ class SolutionScalingKernel : public SolutionScalingAndCheckingKernelBase< real6
 
     // compute pressure scaling factor
     real64 presScalingFactor = 1.0;
-    // when enabled, absolute change scaling has a priority over relative change
+    // the stricter of the absolute and relative change limits is applied
     if( m_maxAbsolutePresChange > 0.0 ) // maxAbsolutePresChange <= 0.0 means that absolute scaling is disabled
     {
       if( absPresChange > m_maxAbsolutePresChange )
@@ -260,12 +260,12 @@ class SolutionScalingKernel : public SolutionScalingAndCheckingKernelBase< real6
         presScalingFactor = m_maxAbsolutePresChange / absPresChange;
       }
     }
-    else if( pres > eps )
+    if( pres > eps )
     {
       real64 const relativePresChange = absPresChange / pres;
       if( relativePresChange > m_maxRelativePresChange )
       {
-        presScalingFactor = m_maxRelativePresChange / relativePresChange;
+        presScalingFactor = LvArray::math::min( presScalingFactor, m_maxRelativePresChange / relativePresChange );
       }
     }
     m_pressureScalingFactor[ei] = presScalingFactor;
diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/ThermalSolutionScalingKernel.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/ThermalSolutionScalingKernel.hpp
index 5f29291ddf1..a548ba12a63 100644
--- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/ThermalSolutionScalingKernel.hpp
+++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/ThermalSolutionScalingKernel.hpp
@@ -134,7 +134,7 @@ class SolutionScalingKernel : public isothermalCompositionalMultiphaseBaseKernel
 
       // compute temperature scaling factor
       real64 tempScalingFactor = 1.0;
-      // when enabled, absolute change scaling has a priority over relative change
+      // the stricter of the absolute and relative change limits is applied
       if( m_maxAbsoluteTempChange > 0.0 ) // maxAbsoluteTempChange <= 0.0 means that absolute scaling is disabled
       {
         if( absTempChange > m_maxAbsoluteTempChange )
@@ -142,12 +142,12 @@ class SolutionScalingKernel : public isothermalCompositionalMultiphaseBaseKernel
           tempScalingFactor = m_maxAbsoluteTempChange / absTempChange;
         }
       }
-      else if( temp > eps )
+      if( temp > eps )
       {
         real64 const relativeTempChange = absTempChange / temp;
         if( relativeTempChange > m_maxRelativeTempChange )
         {
-          tempScalingFactor = m_maxRelativeTempChange / relativeTempChange;
+          tempScalingFactor = LvArray::math::min( tempScalingFactor, m_maxRelativeTempChange / relativeTempChange );
         }
       }
       m_temperatureScalingFactor[ei] = tempScalingFactor;
diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/zFormulation/SolutionScalingZFormulationKernel.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/zFormulation/SolutionScalingZFormulationKernel.hpp
index 47868719506..7130e977577 100644
--- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/zFormulation/SolutionScalingZFormulationKernel.hpp
+++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/zFormulation/SolutionScalingZFormulationKernel.hpp
@@ -233,7 +233,7 @@ class SolutionScalingZFormulationKernel : public SolutionScalingAndCheckingZForm
 
     // compute pressure scaling factor
     real64 presScalingFactor = 1.0;
-    // when enabled, absolute change scaling has a priority over relative change
+    // the stricter of the absolute and relative change limits is applied
     if( m_maxAbsolutePresChange > 0.0 ) // maxAbsolutePresChange <= 0.0 means that absolute scaling is disabled
     {
       if( absPresChange > m_maxAbsolutePresChange )
@@ -241,12 +241,12 @@ class SolutionScalingZFormulationKernel : public SolutionScalingAndCheckingZForm
         presScalingFactor = m_maxAbsolutePresChange / absPresChange;
       }
     }
-    else if( pres > eps )
+    if( pres > eps )
     {
       real64 const relativePresChange = absPresChange / pres;
       if( relativePresChange > m_maxRelativePresChange )
       {
-        presScalingFactor = m_maxRelativePresChange / relativePresChange;
+        presScalingFactor = LvArray::math::min( presScalingFactor, m_maxRelativePresChange / relativePresChange );
       }
     }
     m_pressureScalingFactor[ei] = presScalingFactor;
diff --git a/src/coreComponents/schema/schema.xsd b/src/coreComponents/schema/schema.xsd
index 96b38473f5d..655ac66048c 100644
--- a/src/coreComponents/schema/schema.xsd
+++ b/src/coreComponents/schema/schema.xsd
@@ -6190,10 +6190,6 @@ Information output from lower logLevels is added with the desired log level
 2
  - The summary of declared fields and coupling-->
 		<xsd:attribute name="logLevel" type="integer" default="0" />
-		<!--scalingType => Solution scaling type (Local is only supported for single-phase wells).Valid options:
-* Global
-* Local-->
-		<xsd:attribute name="scalingType" type="geos_flowSolverUtilities_ScalingType" default="Global" />
 		<!--targetRegions => Allowable regions that the solver may be applied to. Note that this does not indicate that the solver will be applied to these regions, only that allocation will occur such that the solver may be applied to these regions. The decision about what regions this solver will beapplied to rests in the EventManager.-->
 		<xsd:attribute name="targetRegions" type="groupNameRef_array" use="required" />
 		<!--timeStepFromTables => Choose time step to honor rates/bhp tables time intervals-->

From a9bcd4406894cce533eac2741fbe787fa15504c7 Mon Sep 17 00:00:00 2001
From: Frank Fei <frankfei0929@gmail.com>
Date: Thu, 5 Nov 2026 15:52:24 -0800
Subject: [PATCH 15/15] Apply suggestion from @jhuang2601

Co-authored-by: Jian Huang <53012159+jhuang2601@users.noreply.github.com>
---
 .../kernels/compositional/ThermalSolutionCheckKernel.hpp        | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/ThermalSolutionCheckKernel.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/ThermalSolutionCheckKernel.hpp
index 04f7b2e37b8..dcf94578b2f 100644
--- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/ThermalSolutionCheckKernel.hpp
+++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/ThermalSolutionCheckKernel.hpp
@@ -111,7 +111,7 @@ class SolutionCheckKernel : public isothermalCompositionalMultiphaseBaseKernels:
     {
       bool const localScaling = m_scalingType == flowSolverUtilities::ScalingType::Local;
       // compute the change in temperature
-      real64 const newTemp = m_temperature[ei] + (localScaling ? m_temperatureScalingFactor[ei] : m_scalingFactor * m_localSolution[stack.localRow + m_temperatureOffset]);
+      real64 const newTemp = m_temperature[ei] + (localScaling ? m_temperatureScalingFactor[ei] : m_scalingFactor) * m_localSolution[stack.localRow + m_temperatureOffset];
       if( newTemp < minTemperature )
       {
         stack.localMinVal = 0;