diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CMakeLists.txt b/src/coreComponents/physicsSolvers/fluidFlow/CMakeLists.txt index d97e010921b..a648e7d8dfe 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 @@ -67,6 +68,7 @@ set( fluidFlowSolvers_headers kernels/singlePhase/ThermalAccumulationKernels.hpp kernels/singlePhase/ThermalDirichletFluxComputeKernel.hpp kernels/singlePhase/ThermalFluxComputeKernel.hpp + kernels/singlePhase/ThermalSolutionScalingKernel.hpp kernels/singlePhase/proppant/ProppantBaseKernels.hpp kernels/singlePhase/proppant/ProppantFluxKernels.hpp kernels/singlePhase/reactive/AccumulationKernels.hpp diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.cpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.cpp index 235a129d64c..b6191b9942f 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 ) @@ -613,6 +613,14 @@ real64 CompositionalMultiphaseFVM::scalingForSystemSolution( DomainPartition & d { GEOS_MARK_FUNCTION; + // In a coupled run the coupled solver's NonlinearSolverParameters overwrite the flow subsolver, + // so the lineSearchAction / scalingType compatibility is checked here at run time. + // A better fix is to check it at input time in CoupledSolver. + GEOS_ERROR_IF( m_scalingType == ScalingType::Local && + m_nonlinearSolverParameters.m_lineSearchAction != NonlinearSolverParameters::LineSearchAction::None, + "line search is not supported for scalingType = Local; set lineSearchAction to None, including in the coupled solver", + getDataContext() ); + if( m_formulationType == CompositionalMultiphaseFormulationType::OverallComposition ) { return scalingForSystemSolutionZFormulation( domain, dofManager, localSolution ); @@ -651,6 +659,7 @@ real64 CompositionalMultiphaseFVM::scalingForSystemSolution( DomainPartition & d createAndLaunch< parallelDevicePolicy<> >( m_maxRelativePresChange, m_maxAbsolutePresChange, m_maxRelativeTempChange, + m_maxAbsoluteTempChange, m_maxCompFracChange, m_maxRelativeCompDensChange, pressure, 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 be09af79934..a9332402e9f 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/FlowSolverBase.cpp b/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.cpp index 412e1d4e6b4..f0234d014c7 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.cpp @@ -168,6 +168,12 @@ FlowSolverBase::FlowSolverBase( string const & name, setApplyDefaultValue( -1.0 ). // disabled by default setDescription( "Maximum (absolute) pressure change in a Newton iteration" ); + this->registerWrapper( viewKeyStruct::maxAbsoluteTempChangeString(), &m_maxAbsoluteTempChange ). + setSizedFromParent( 0 ). + setInputFlag( InputFlags::OPTIONAL ). + setApplyDefaultValue( -1.0 ). // disabled by default + setDescription( "Maximum (absolute) temperature change in a Newton iteration" ); + this->registerWrapper( viewKeyStruct::maxSequentialPresChangeString(), &m_maxSequentialPresChange ). setSizedFromParent( 0 ). setInputFlag( InputFlags::OPTIONAL ). diff --git a/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.hpp b/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.hpp index e2b99c482ea..c2a0b0bd984 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"; } @@ -308,6 +309,9 @@ class FlowSolverBase : public PhysicsSolverBase /// maximum (absolute) pressure change in a Newton iteration real64 m_maxAbsolutePresChange; + /// maximum (absolute) temperature change in a Newton iteration + real64 m_maxAbsoluteTempChange; + /// maximum (absolute) pressure change in a sequential iteration real64 m_sequentialPresChange; real64 m_maxSequentialPresChange; 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 cf140d68c07..b620a1eb3f5 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp @@ -44,6 +44,7 @@ #include "physicsSolvers/fluidFlow/kernels/singlePhase/MobilityKernel.hpp" #include "physicsSolvers/fluidFlow/kernels/singlePhase/SolutionCheckKernel.hpp" #include "physicsSolvers/fluidFlow/kernels/singlePhase/SolutionScalingKernel.hpp" +#include "physicsSolvers/fluidFlow/kernels/singlePhase/ThermalSolutionScalingKernel.hpp" #include "physicsSolvers/fluidFlow/kernels/singlePhase/StatisticsKernel.hpp" #include "physicsSolvers/fluidFlow/kernels/singlePhase/HydrostaticPressureKernel.hpp" #include "physicsSolvers/fluidFlow/kernels/singlePhase/ThermalHydrostaticPressureKernel.hpp" @@ -61,7 +62,8 @@ using namespace singlePhaseBaseKernels; SinglePhaseBase::SinglePhaseBase( const string & name, Group * const parent ): - FlowSolverBase( name, parent ) + FlowSolverBase( name, parent ), + m_scalingType( flowSolverUtilities::ScalingType::Global ) { this->registerWrapper( viewKeyStruct::inputTemperatureString(), &m_inputTemperature ). setApplyDefaultValue( 0.0 ). @@ -103,9 +105,12 @@ void SinglePhaseBase::registerDataOnMesh( Group & meshBodies ) subRegion.registerField< flow::mass_n >( getName() ); subRegion.registerField< flow::dMass >( getName() ).reference().resizeDimension< 1 >( m_numDofPerCell ); + subRegion.registerField< flow::pressureScalingFactor >( getName() ); + if( m_isThermal ) { subRegion.registerField< flow::dEnergy >( getName() ).reference().resizeDimension< 1 >( m_numDofPerCell ); + subRegion.registerField< flow::temperatureScalingFactor >( getName() ); } } ); @@ -1300,37 +1305,110 @@ real64 SinglePhaseBase::scalingForSystemSolution( DomainPartition & domain, { GEOS_MARK_FUNCTION; + // In a coupled run the coupled solver's NonlinearSolverParameters overwrite the flow subsolver, + // so the lineSearchAction / scalingType compatibility is checked here at run time. + // A better fix is to check it at input time in CoupledSolver. + GEOS_ERROR_IF( m_scalingType == flowSolverUtilities::ScalingType::Local && + m_nonlinearSolverParameters.m_lineSearchAction != NonlinearSolverParameters::LineSearchAction::None, + "line search is not supported for scalingType = Local; set lineSearchAction to None, including in the coupled solver", + getDataContext() ); + string const dofKey = dofManager.getKey( viewKeyStruct::elemDofFieldString() ); real64 scalingFactor = 1.0; real64 maxDeltaPres = 0.0; + real64 maxDeltaTemp = 0.0; - forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &, - MeshLevel & mesh, - string_array const & regionNames ) + real64 minPresScalingFactor = 1.0, minTempScalingFactor = 1.0; + + if( m_isThermal ) { - mesh.getElemManager().forElementSubRegions( regionNames, - [&]( localIndex const, - ElementSubRegionBase & subRegion ) + forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &, + MeshLevel & mesh, + string_array const & regionNames ) { - globalIndex const rankOffset = dofManager.rankOffset(); - arrayView1d< globalIndex const > const dofNumber = subRegion.getReference< array1d< globalIndex > >( dofKey ); - arrayView1d< integer const > const ghostRank = subRegion.ghostRank(); + mesh.getElemManager().forElementSubRegions( regionNames, + [&]( localIndex const, + ElementSubRegionBase & subRegion ) + { + globalIndex const rankOffset = dofManager.rankOffset(); + arrayView1d< globalIndex const > const dofNumber = subRegion.getReference< array1d< globalIndex > >( dofKey ); + arrayView1d< integer const > const ghostRank = subRegion.ghostRank(); + + arrayView1d< real64 > pressureScalingFactor = subRegion.getField< flow::pressureScalingFactor >(); + arrayView1d< real64 > temperatureScalingFactor = subRegion.getField< flow::temperatureScalingFactor >(); + + integer const tempDofOffset = 1; + auto const subRegionData = thermalSinglePhaseBaseKernels:: + SolutionScalingKernel:: + launch< parallelDevicePolicy<> >( localSolution, rankOffset, tempDofOffset, dofNumber, ghostRank, + m_maxAbsolutePresChange, m_maxAbsoluteTempChange, + pressureScalingFactor, temperatureScalingFactor ); - auto const subRegionData = - singlePhaseBaseKernels::SolutionScalingKernel:: - launch< parallelDevicePolicy<> >( localSolution, rankOffset, dofNumber, ghostRank, m_maxAbsolutePresChange ); + if( m_scalingType == flowSolverUtilities::ScalingType::Global ) + { + scalingFactor = std::min( scalingFactor, subRegionData.scalingFactor ); + } + maxDeltaPres = std::max( maxDeltaPres, subRegionData.maxDeltaPres ); + maxDeltaTemp = std::max( maxDeltaTemp, subRegionData.maxDeltaTemp ); - scalingFactor = std::min( scalingFactor, subRegionData.first ); - maxDeltaPres = std::max( maxDeltaPres, subRegionData.second ); + minPresScalingFactor = std::min( minPresScalingFactor, subRegionData.minPresScalingFactor ); + minTempScalingFactor = std::min( minTempScalingFactor, subRegionData.minTempScalingFactor ); + } ); } ); - } ); + } + else + { + GEOS_UNUSED_VAR( maxDeltaTemp ); + + forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &, + MeshLevel & mesh, + string_array const & regionNames ) + { + mesh.getElemManager().forElementSubRegions( regionNames, + [&]( localIndex const, + ElementSubRegionBase & subRegion ) + { + globalIndex const rankOffset = dofManager.rankOffset(); + arrayView1d< globalIndex const > const dofNumber = subRegion.getReference< array1d< globalIndex > >( dofKey ); + arrayView1d< integer const > const ghostRank = subRegion.ghostRank(); + + arrayView1d< real64 > pressureScalingFactor = subRegion.getField< flow::pressureScalingFactor >(); + + auto const subRegionData = singlePhaseBaseKernels:: + SolutionScalingKernel:: + launch< parallelDevicePolicy<> >( localSolution, rankOffset, dofNumber, ghostRank, + m_maxAbsolutePresChange, pressureScalingFactor ); + + if( m_scalingType == flowSolverUtilities::ScalingType::Global ) + { + scalingFactor = std::min( scalingFactor, subRegionData.scalingFactor ); + } + minPresScalingFactor = std::min( minPresScalingFactor, subRegionData.minPresScalingFactor ); + maxDeltaPres = std::max( maxDeltaPres, subRegionData.maxDeltaPres ); + } ); + } ); + } scalingFactor = MpiWrapper::min( scalingFactor ); + minPresScalingFactor = MpiWrapper::min( minPresScalingFactor ); maxDeltaPres = MpiWrapper::max( maxDeltaPres ); GEOS_LOG_LEVEL_RANK_0( logInfo::Solution, GEOS_FMT( " {}: Max pressure change = {} Pa (before scaling)", getName(), fmt::format( "{:.{}f}", maxDeltaPres, 3 ) ) ); + GEOS_LOG_LEVEL_RANK_0( logInfo::Solution, GEOS_FMT( " {}: Min pressure scaling factor = {}", getName(), minPresScalingFactor ) ); + + if( m_isThermal ) + { + minTempScalingFactor = MpiWrapper::min( minTempScalingFactor ); + maxDeltaTemp = MpiWrapper::max( maxDeltaTemp ); + + GEOS_LOG_LEVEL_RANK_0( logInfo::Solution, GEOS_FMT( " {}: Max temperature change = {} K (before scaling)", + getName(), fmt::format( "{:.{}f}", maxDeltaTemp, 3 ) ) ); + + GEOS_LOG_LEVEL_RANK_0( logInfo::Solution, GEOS_FMT( " {}: Min temperature scaling factor = {}", getName(), minTempScalingFactor ) ); + } + return scalingFactor; } @@ -1358,6 +1436,7 @@ bool SinglePhaseBase::checkSystemSolution( DomainPartition & domain, arrayView1d< globalIndex const > const dofNumber = subRegion.getReference< array1d< globalIndex > >( dofKey ); arrayView1d< integer const > const ghostRank = subRegion.ghostRank(); arrayView1d< real64 const > const pres = subRegion.getField< flow::pressure >(); + arrayView1d< real64 const > const pressureScalingFactor = subRegion.getField< flow::pressureScalingFactor >(); auto const negPresCollector = rankNegPressureIds.createCollector( subRegion.localToGlobalMap().toViewConst() ); auto const statistics = singlePhaseBaseKernels::SolutionCheckKernel:: @@ -1366,7 +1445,9 @@ bool SinglePhaseBase::checkSystemSolution( DomainPartition & domain, dofNumber, ghostRank, pres, - scalingFactor, + m_scalingType, + scalingFactor, // global scaling factor + pressureScalingFactor, // per-element scaling factors negPresCollector ); minPressure = std::min( minPressure, statistics.minNegPres ); 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 ad2a994a364..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(), - scalingFactor, + flow::pressureScalingFactor::key(), pressureMask ); - - dofManager.addVectorToField( localSolution, - BASE::viewKeyStruct::elemDofFieldString(), - flow::temperature::key(), - scalingFactor, - 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/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/ThermalSolutionCheckKernel.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/ThermalSolutionCheckKernel.hpp index 10566b06827..dcf94578b2f 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,9 +109,9 @@ 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]); + real64 const newTemp = m_temperature[ei] + (localScaling ? m_temperatureScalingFactor[ei] : m_scalingFactor) * m_localSolution[stack.localRow + m_temperatureOffset]; if( newTemp < minTemperature ) { stack.localMinVal = 0; @@ -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/compositional/ThermalSolutionScalingKernel.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/ThermalSolutionScalingKernel.hpp index 1a4a2b4c782..a548ba12a63 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; - + // compute temperature scaling factor + real64 tempScalingFactor = 1.0; + // 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 ) + { + tempScalingFactor = m_maxAbsoluteTempChange / absTempChange; + } + } 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 = LvArray::math::min( 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/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/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/kernels/singlePhase/SolutionScalingKernel.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/SolutionScalingKernel.hpp index fe0572da818..5da03150d63 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/SolutionScalingKernel.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/SolutionScalingKernel.hpp @@ -31,14 +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 ) + 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 ); @@ -47,6 +66,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,13 +78,18 @@ struct SolutionScalingKernel if( maxAbsolutePresChange > 0.0 && absPresChange > maxAbsolutePresChange ) { real64 const presScalingFactor = maxAbsolutePresChange / absPresChange; + pressureScalingFactor[ei] = presScalingFactor; scalingFactor.min( presScalingFactor ); } } } ); - 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 new file mode 100644 index 00000000000..5b6a08be140 --- /dev/null +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/ThermalSolutionScalingKernel.hpp @@ -0,0 +1,111 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * 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" +#include "physicsSolvers/fluidFlow/kernels/singlePhase/SolutionScalingKernel.hpp" + +namespace geos +{ + +namespace thermalSinglePhaseBaseKernels +{ + +/******************************** SolutionScalingKernel ********************************/ + +struct SolutionScalingKernel +{ + template< typename POLICY > + 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 ); + 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 ); + } + } + + } ); + + 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; + } + +}; + +} // 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 77bf2821028..ce8cd45c5f7 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.cpp @@ -137,7 +137,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 ) { @@ -164,12 +163,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 ). @@ -1615,6 +1608,7 @@ CompositionalMultiphaseWell::scalingForLocalSystemSolution( WellElementSubRegion createAndLaunch< parallelDevicePolicy<> >( m_maxRelativePresChange, m_maxAbsolutePresChange, m_maxRelativeTempChange, + m_maxAbsoluteTempChange, m_maxCompFracChange, m_maxRelativeCompDensChange, pressure, @@ -1769,7 +1763,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 >(); diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.hpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.hpp index de23fec8bfc..b8da93809c3 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.hpp @@ -357,8 +357,6 @@ class CompositionalMultiphaseWell : public WellControls 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"; } @@ -440,9 +438,6 @@ class CompositionalMultiphaseWell : public WellControls /// 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 7b50cd11d1c..111d1cd43a0 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/SinglePhaseWell.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/SinglePhaseWell.cpp @@ -53,6 +53,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/SinglePhaseStatisticsAggregator.hpp" #include "physicsSolvers/fluidFlow/wells/kernels/SinglePhaseWellConstraintKernels.hpp" #include "physicsSolvers/multiphysics/CoupledReservoirAndSinglePhaseWellKernels.hpp" @@ -162,6 +164,12 @@ void SinglePhaseWell::registerWellDataOnMesh( WellElementSubRegion & subRegion ) perforationData.registerField< well::gravityCoefficient >( getName() ); } + subRegion.registerField< well::pressureScalingFactor >( getName() ); + if( isThermal() ) + { + subRegion.registerField< well::temperatureScalingFactor >( getName() ); + } + // write rates output header if( m_writeCSV > 0 && subRegion.isLocallyOwned()) { @@ -936,17 +944,74 @@ void SinglePhaseWell::computeWellPerforationRates( real64 const & time_n, } } } + +real64 +SinglePhaseWell::scalingForLocalSystemSolution( WellElementSubRegion & subRegion, + DofManager const & dofManager, + real64 & maxDeltaPres, + real64 & maxDeltaTemp, + real64 & minPresScalingFactor, + real64 & minTempScalingFactor, + arrayView1d< real64 const > const & localSolution ) +{ + GEOS_MARK_FUNCTION; + + maxDeltaTemp = 0.0; + minTempScalingFactor = 1.0; + + string const wellDofKey = dofManager.getKey( wellElementDofName() ); + globalIndex const rankOffset = dofManager.rankOffset(); + 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 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 = subRegionData.scalingFactor; + maxDeltaPres = subRegionData.maxDeltaPres; + maxDeltaTemp = subRegionData.maxDeltaTemp; + minPresScalingFactor = subRegionData.minPresScalingFactor; + minTempScalingFactor = subRegionData.minTempScalingFactor; + } + else + { + auto const subRegionData = singlePhaseBaseKernels:: + SolutionScalingKernel:: + launch< parallelDevicePolicy<> >( localSolution, rankOffset, dofNumber, ghostRank, + m_maxAbsolutePresChange, pressureScalingFactor ); + + scalingFactor = subRegionData.scalingFactor; + maxDeltaPres = subRegionData.maxDeltaPres; + minPresScalingFactor = subRegionData.minPresScalingFactor; + } + + return scalingFactor; +} + real64 SinglePhaseWell::scalingForWellSystemSolution( WellElementSubRegion & subRegion, DofManager const & dofManager, arrayView1d< real64 const > const & localSolution ) { - GEOS_MARK_FUNCTION; - GEOS_UNUSED_VAR( subRegion ); - GEOS_UNUSED_VAR( dofManager ); - GEOS_UNUSED_VAR( localSolution ); + real64 maxDeltaPres = 0.0, maxDeltaTemp = 0.0; + real64 minPresScalingFactor = 1.0, minTempScalingFactor = 1.0; - return 1.0; + return scalingForLocalSystemSolution( subRegion, dofManager, + maxDeltaPres, maxDeltaTemp, + minPresScalingFactor, minTempScalingFactor, + localSolution ); } array1d< real64 > @@ -1116,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 ); diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/SinglePhaseWell.hpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/SinglePhaseWell.hpp index 12d1697db8d..c61ed53cd82 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/SinglePhaseWell.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/SinglePhaseWell.hpp @@ -176,6 +176,25 @@ class SinglePhaseWell : public WellControls DofManager const & dofManager, arrayView1d< real64 const > const & localRhs ) override; + /** + * @brief Compute the Newton update scaling factor of one well, with the quantities reported in the log + * @param[in] subRegion the well subregion + * @param[in] dofManager degree-of-freedom manager associated with the linear system + * @param[out] maxDeltaPres the maximum absolute pressure change (before scaling) + * @param[out] maxDeltaTemp the maximum absolute temperature change (before scaling) + * @param[out] minPresScalingFactor the minimum pressure scaling factor + * @param[out] minTempScalingFactor the minimum temperature scaling factor + * @param[in] localSolution the solution vector + * @return the scaling factor, local to this rank + */ + real64 scalingForLocalSystemSolution( WellElementSubRegion & subRegion, + DofManager const & dofManager, + real64 & maxDeltaPres, + real64 & maxDeltaTemp, + real64 & minPresScalingFactor, + real64 & minTempScalingFactor, + arrayView1d< real64 const > const & localSolution ); + virtual real64 scalingForWellSystemSolution( WellElementSubRegion & subRegion, DofManager const & dofManager, arrayView1d< real64 const > const & localSolution ) override; diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/WellControls.cpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/WellControls.cpp index c55bcd2c2b7..d3d85043987 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/WellControls.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/WellControls.cpp @@ -73,6 +73,8 @@ WellControls::WellControls( string const & name, Group * const parent ) m_estimateSolution( 0 ), m_enableIsoThermalEstimator( 0 ), m_thermalEffectsEnabled( false ), + m_maxAbsolutePresChange( -1.0 ), // disabled by default + m_maxAbsoluteTempChange( -1.0 ), // disabled by default /// Nonlinear solver parameters m_wellNewtonSolver( groupKeyStruct::wellNewtonSolverString(), this ), m_estimatorDoFManager( name ), @@ -153,6 +155,18 @@ WellControls::WellControls( string const & name, Group * const parent ) setInputFlag( InputFlags::OPTIONAL ). setDescription( "Flag to enable isothermal estimator prior to coupled reservoir and well solve." ); + 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" ); + registerWrapper( viewKeyStruct::statusTableNameString(), &m_statusTableName ). setRTTypeName( rtTypes::CustomTypes::groupNameRef ). setInputFlag( InputFlags::OPTIONAL ). diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/WellControls.hpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/WellControls.hpp index afadd22a724..cb2ee1a2e0f 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/WellControls.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/WellControls.hpp @@ -739,6 +739,11 @@ class WellControls : public dataRepository::Group /// string key for the enable iso thermal estimator flag static constexpr char const * enableIsoThermalEstimatorString() { return "enableIsoThermalEstimator"; } + /// string key for the maximum absolute pressure change in a Newton iteration + static constexpr char const * maxAbsolutePresChangeString() { return "maxAbsolutePressureChange"; } + /// string key for the maximum absolute temperature change in a Newton iteration + static constexpr char const * maxAbsoluteTempChangeString() { return "maxAbsoluteTemperatureChange"; } + // control data (not registered on the mesh) static constexpr char const * massDensityString() { return "massDensity";} @@ -1012,6 +1017,12 @@ class WellControls : public dataRepository::Group integer m_enableIsoThermalEstimator; bool m_thermalEffectsEnabled; + /// maximum (absolute) change in pressure between two Newton iterations + real64 m_maxAbsolutePresChange; + + /// maximum (absolute) change in temperature between two Newton iterations + real64 m_maxAbsoluteTempChange; + WellNewtonSolver m_wellNewtonSolver; diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/WellManager.cpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/WellManager.cpp index 6bf0ad13d28..af95b0d19f4 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/WellManager.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/WellManager.cpp @@ -717,8 +717,61 @@ WellManager::scalingForSystemSolution( DomainPartition & domain, } else { - // Single phase well scaling- not implemented yet - scalingFactor=1.0; + real64 maxDeltaPres = 0.0, maxDeltaTemp = 0.0; + real64 minPresScalingFactor = 1.0, minTempScalingFactor = 1.0; + real64 localMaxDeltaPres = 0.0, localMaxDeltaTemp = 0.0; + real64 localMinPresScalingFactor = 1.0, localMinTempScalingFactor = 1.0; + + forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &, + MeshLevel & mesh, + string_array const & regionNames ) + { + + ElementRegionManager & elemManager = mesh.getElemManager(); + elemManager.forElementSubRegions< WellElementSubRegion >( regionNames, + [&]( localIndex const, + WellElementSubRegion & subRegion ) + + { + SinglePhaseWell * wellControls = dynamic_cast< SinglePhaseWell * >(&getWellControls ( subRegion )); + localScalingFactor = wellControls->scalingForLocalSystemSolution( subRegion, + dofManager, + localMaxDeltaPres, + localMaxDeltaTemp, + localMinPresScalingFactor, + localMinTempScalingFactor, + localSolution ); + maxDeltaPres = LvArray::math::max( localMaxDeltaPres, maxDeltaPres ); + maxDeltaTemp = LvArray::math::max( localMaxDeltaTemp, maxDeltaTemp ); + minPresScalingFactor = LvArray::math::min( localMinPresScalingFactor, minPresScalingFactor ); + minTempScalingFactor = LvArray::math::min( localMinTempScalingFactor, minTempScalingFactor ); + scalingFactor = LvArray::math::min( localScalingFactor, scalingFactor ); + + } ); + } ); + + scalingFactor = MpiWrapper::min( scalingFactor ); + maxDeltaPres = MpiWrapper::max( maxDeltaPres ); + minPresScalingFactor = MpiWrapper::min( minPresScalingFactor ); + + 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 well pressure scaling factor: {}", + getName(), minPresScalingFactor ) ); + + if( m_isThermal ) + { + maxDeltaTemp = MpiWrapper::max( maxDeltaTemp ); + minTempScalingFactor = MpiWrapper::min( minTempScalingFactor ); + 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 well temperature scaling factor: {}", + getName(), minTempScalingFactor ) ); + } } return LvArray::math::max( scalingFactor, m_minScalingFactor ); diff --git a/src/coreComponents/schema/schema.xsd b/src/coreComponents/schema/schema.xsd index 5be40a17911..655ac66048c 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 + + @@ -3326,7 +3328,7 @@ Information output from lower logLevels is added with the desired log level - + @@ -3379,7 +3381,7 @@ When set to `all` output both convergence & iteration information to a csv.--> - + @@ -3420,6 +3422,8 @@ Information output from lower logLevels is added with the desired log level + + @@ -4138,6 +4142,8 @@ Information output from lower logLevels is added with the desired log level + + @@ -4713,6 +4719,8 @@ Information output from lower logLevels is added with the desired log level + + @@ -4777,6 +4785,8 @@ Information output from lower logLevels is added with the desired log level + + @@ -4888,10 +4898,16 @@ Information output from lower logLevels is added with the desired log level + + + + @@ -4943,6 +4959,8 @@ Information output from lower logLevels is added with the desired log level + + @@ -5296,10 +5314,16 @@ Information output from lower logLevels is added with the desired log level + + + + @@ -5357,6 +5381,8 @@ Information output from lower logLevels is added with the desired log level + + @@ -6223,6 +6249,8 @@ Information output from lower logLevels is added with the desired log level + + @@ -6473,6 +6501,10 @@ Information output from lower logLevels is added with the desired log level 1 - Well control information--> + + + +