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..5fbf1cdbfc5 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,34 +1292,90 @@ real64 SinglePhaseBase::scalingForSystemSolution( DomainPartition & domain, string const dofKey = dofManager.getKey( viewKeyStruct::elemDofFieldString() ); real64 scalingFactor = 1.0; real64 maxDeltaPres = 0.0; + real64 maxDeltaTemp = 0.0; - forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &, - MeshLevel & mesh, - string_array const & regionNames ) + real64 minPresScalingFactor = 1.0, minTempScalingFactor = 1.0; + + if( m_isThermal ) { - mesh.getElemManager().forElementSubRegions( regionNames, - [&]( localIndex const, - ElementSubRegionBase & subRegion ) + forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &, + MeshLevel & mesh, + string_array const & regionNames ) { - globalIndex const rankOffset = dofManager.rankOffset(); - arrayView1d< globalIndex const > const dofNumber = subRegion.getReference< array1d< globalIndex > >( dofKey ); - arrayView1d< integer const > const ghostRank = subRegion.ghostRank(); + mesh.getElemManager().forElementSubRegions( regionNames, + [&]( localIndex const, + ElementSubRegionBase & subRegion ) + { + globalIndex const rankOffset = dofManager.rankOffset(); + arrayView1d< globalIndex const > const dofNumber = subRegion.getReference< array1d< globalIndex > >( dofKey ); + arrayView1d< integer const > const ghostRank = subRegion.ghostRank(); + + arrayView1d< real64 > pressureScalingFactor = subRegion.getField< flow::pressureScalingFactor >(); + arrayView1d< real64 > temperatureScalingFactor = subRegion.getField< flow::temperatureScalingFactor >(); + + integer const tempDofOffset = 1; + auto const subRegionData = thermalSinglePhaseBaseKernels:: + SolutionScalingKernel:: + launch< parallelDevicePolicy<> >( localSolution, rankOffset, tempDofOffset, dofNumber, ghostRank, + m_maxAbsolutePresChange, m_maxAbsoluteTempChange, + pressureScalingFactor, temperatureScalingFactor ); + + scalingFactor = std::min( scalingFactor, std::get< 0 >( subRegionData ) ); + maxDeltaPres = std::max( maxDeltaPres, std::get< 1 >( subRegionData ) ); + maxDeltaTemp = std::max( maxDeltaTemp, std::get< 2 >( subRegionData ) ); + + minPresScalingFactor = std::min( minPresScalingFactor, std::get< 3 >( subRegionData ) ); + minTempScalingFactor = std::min( minTempScalingFactor, std::get< 4 >( subRegionData ) ); + } ); + } ); + } + else + { + GEOS_UNUSED_VAR( maxDeltaTemp ); + + forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &, + MeshLevel & mesh, + string_array const & regionNames ) + { + mesh.getElemManager().forElementSubRegions( regionNames, + [&]( localIndex const, + ElementSubRegionBase & subRegion ) + { + globalIndex const rankOffset = dofManager.rankOffset(); + arrayView1d< globalIndex const > const dofNumber = subRegion.getReference< array1d< globalIndex > >( dofKey ); + arrayView1d< integer const > const ghostRank = subRegion.ghostRank(); - auto const subRegionData = - singlePhaseBaseKernels::SolutionScalingKernel:: - launch< parallelDevicePolicy<> >( localSolution, rankOffset, dofNumber, ghostRank, m_maxAbsolutePresChange ); + auto const subRegionData = singlePhaseBaseKernels:: + SolutionScalingKernel:: + launch< parallelDevicePolicy<> >( localSolution, rankOffset, dofNumber, ghostRank, m_maxAbsolutePresChange ); - scalingFactor = std::min( scalingFactor, subRegionData.first ); - maxDeltaPres = std::max( maxDeltaPres, subRegionData.second ); + scalingFactor = std::min( scalingFactor, subRegionData.first ); + minPresScalingFactor = std::min( minPresScalingFactor, subRegionData.first ); + maxDeltaPres = std::max( maxDeltaPres, subRegionData.second ); + } ); } ); - } ); + } scalingFactor = MpiWrapper::min( scalingFactor ); + minPresScalingFactor = MpiWrapper::min( minPresScalingFactor ); maxDeltaPres = MpiWrapper::max( maxDeltaPres ); GEOS_LOG_LEVEL_RANK_0( logInfo::Solution, GEOS_FMT( " {}: Max pressure change = {} Pa (before scaling)", getName(), fmt::format( "{:.{}f}", maxDeltaPres, 3 ) ) ); + GEOS_LOG_LEVEL_RANK_0( logInfo::Solution, GEOS_FMT( " {}: Min pressure scaling factor = {}", getName(), minPresScalingFactor ) ); + + if( m_isThermal ) + { + minTempScalingFactor = MpiWrapper::min( minTempScalingFactor ); + maxDeltaTemp = MpiWrapper::max( maxDeltaTemp ); + + GEOS_LOG_LEVEL_RANK_0( logInfo::Solution, GEOS_FMT( " {}: Max temperature change = {} K (before scaling)", + getName(), fmt::format( "{:.{}f}", maxDeltaTemp, 3 ) ) ); + + GEOS_LOG_LEVEL_RANK_0( logInfo::Solution, GEOS_FMT( " {}: Min temperature scaling factor = {}", getName(), minTempScalingFactor ) ); + } + return scalingFactor; } 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 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..5ffc3ff7ed7 --- /dev/null +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/ThermalSolutionScalingKernel.hpp @@ -0,0 +1,103 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2024 TotalEnergies + * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2023-2024 Chevron + * Copyright (c) 2019- GEOS/GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +/** + * @file ThermalSolutionScalingKernel.hpp + */ + +#ifndef GEOS_PHYSICSSOLVERS_FLUIDFLOW_SINGLEPHASE_THERMALSOLUTIONSCALINGKERNEL_HPP +#define GEOS_PHYSICSSOLVERS_FLUIDFLOW_SINGLEPHASE_THERMALSOLUTIONSCALINGKERNEL_HPP + +#include "common/DataTypes.hpp" +#include "common/GEOS_RAJA_Interface.hpp" + +namespace geos +{ + +namespace thermalSinglePhaseBaseKernels +{ + +/******************************** SolutionScalingKernel ********************************/ + +struct SolutionScalingKernel +{ + template< typename POLICY > + static std::tuple< real64, real64, real64, real64, real64 > launch( arrayView1d< real64 const > const & localSolution, + globalIndex const rankOffset, + globalIndex const temperatureOffset, + arrayView1d< globalIndex const > const & dofNumber, + arrayView1d< integer const > const & ghostRank, + real64 const maxAbsolutePresChange, + real64 const maxAbsoluteTempChange, + arrayView1d< real64 > pressureScalingFactor, + arrayView1d< real64 > temperatureScalingFactor ) + { + RAJA::ReduceMin< ReducePolicy< POLICY >, real64 > scalingFactor( 1.0 ); + RAJA::ReduceMax< ReducePolicy< POLICY >, real64 > maxDeltaPres( 0.0 ); + RAJA::ReduceMax< ReducePolicy< POLICY >, real64 > maxDeltaTemp( 0.0 ); + + RAJA::ReduceMin< ReducePolicy< POLICY >, real64 > localMinPresScalingFactor( 1.0 ); + RAJA::ReduceMin< ReducePolicy< POLICY >, real64 > localMinTempScalingFactor( 1.0 ); + + forAll< POLICY >( dofNumber.size(), [=] GEOS_HOST_DEVICE ( localIndex const ei ) mutable + { + if( ghostRank[ei] < 0 && dofNumber[ei] >= 0 ) + { + pressureScalingFactor[ei] = 1.0; + temperatureScalingFactor[ei] = 1.0; + + localIndex const lid = dofNumber[ei] - rankOffset; + + // compute the change in pressure + real64 const absPresChange = LvArray::math::abs( localSolution[lid] ); + maxDeltaPres.max( absPresChange ); + + // compute the change in temperature + real64 const absTempChange = LvArray::math::abs( localSolution[lid + temperatureOffset] ); + maxDeltaTemp.max( absTempChange ); + + // maxAbsolutePresChange <= 0.0 means that scaling is disabled, and we are only collecting maxDeltaPres in this kernel + if( maxAbsolutePresChange > 0.0 && absPresChange > maxAbsolutePresChange ) + { + real64 const presScalingFactor = maxAbsolutePresChange / absPresChange; + pressureScalingFactor[ei] = presScalingFactor; + scalingFactor.min( presScalingFactor ); + + localMinPresScalingFactor.min( presScalingFactor ); + } + + // maxAbsoluteTempChange <= 0.0 means that scaling is disabled, and we are only collecting maxDeltaTemps in this kernel + if( maxAbsoluteTempChange > 0.0 && absTempChange > maxAbsoluteTempChange ) + { + real64 const tempScalingFactor = maxAbsoluteTempChange / absTempChange; + temperatureScalingFactor[ei] = tempScalingFactor; + scalingFactor.min( tempScalingFactor ); + + localMinTempScalingFactor.min( tempScalingFactor ); + } + } + + } ); + + return { scalingFactor.get(), maxDeltaPres.get(), maxDeltaTemp.get(), localMinPresScalingFactor.get(), localMinTempScalingFactor.get() }; + } + +}; + +} // namespace thermalSinglePhaseBaseKernels + +} // namespace geos + +#endif //GEOS_PHYSICSSOLVERS_FLUIDFLOW_SINGLEPHASE_THERMALSOLUTIONSCALINGKERNEL_HPP 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..417e11caa0c 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,103 @@ 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 >(); + + integer const tempDofOffset = 2; + auto const subRegionData = thermalSinglePhaseBaseKernels:: + SolutionScalingKernel:: + launch< parallelDevicePolicy<> >( localSolution, rankOffset, tempDofOffset, dofNumber, ghostRank, + m_maxAbsolutePresChange, m_maxAbsoluteTempChange, + pressureScalingFactor, temperatureScalingFactor ); + + scalingFactor = std::min( scalingFactor, std::get< 0 >( subRegionData ) ); + maxDeltaPres = std::max( maxDeltaPres, std::get< 1 >( subRegionData ) ); + maxDeltaTemp = std::max( maxDeltaTemp, std::get< 2 >( subRegionData ) ); + + minPresScalingFactor = std::min( minPresScalingFactor, std::get< 3 >( subRegionData ) ); + minTempScalingFactor = std::min( minTempScalingFactor, std::get< 4 >( subRegionData ) ); + } ); + } ); + } + else + { + 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, @@ -1075,7 +1178,7 @@ SinglePhaseWell::applySystemSolution( DofManager const & dofManager, dofManager.addVectorToField( localSolution, wellElementDofName(), well::pressure::key(), - scalingFactor, + well::pressureScalingFactor::key(), pressureMask ); dofManager.addVectorToField( localSolution, @@ -1091,7 +1194,7 @@ SinglePhaseWell::applySystemSolution( DofManager const & dofManager, dofManager.addVectorToField( localSolution, wellElementDofName(), fields::well::temperature::key(), - scalingFactor, + well::temperatureScalingFactor::key(), temperatureMask ); } 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; }; }