From c694c12861afff777b9aab3812a438ff50b8ca0c Mon Sep 17 00:00:00 2001 From: frankfeifan Date: Sun, 22 Feb 2026 13:03:08 -0800 Subject: [PATCH 1/5] 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 Date: Sun, 22 Feb 2026 14:35:21 -0800 Subject: [PATCH 2/5] 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 Date: Sun, 22 Feb 2026 16:35:00 -0800 Subject: [PATCH 3/5] 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 Date: Mon, 23 Feb 2026 13:35:25 -0800 Subject: [PATCH 4/5] 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 Date: Mon, 23 Feb 2026 13:54:50 -0800 Subject: [PATCH 5/5] 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 );