From 90260a1509d65eb73a698f292aafc7a66d2f537b Mon Sep 17 00:00:00 2001 From: dkachuma Date: Wed, 20 May 2026 19:47:02 -0500 Subject: [PATCH 01/12] Remove constraint containers --- .../wells/CompositionalMultiphaseWell.cpp | 10 +- .../fluidFlow/wells/SinglePhaseWell.cpp | 9 + .../fluidFlow/wells/WellBHPConstraints.cpp | 5 +- .../fluidFlow/wells/WellBHPConstraints.hpp | 4 +- .../fluidFlow/wells/WellConstraintsBase.cpp | 2 +- .../fluidFlow/wells/WellConstraintsBase.hpp | 4 +- .../fluidFlow/wells/WellControls.cpp | 332 +++++++++--------- .../fluidFlow/wells/WellControls.hpp | 134 ++++--- .../fluidFlow/wells/WellManager.hpp | 106 ------ .../CompositionalMultiphaseWellKernels.cpp | 9 +- .../CompositionalMultiphaseWellKernels.hpp | 30 +- .../wells/kernels/SinglePhaseWellKernels.cpp | 20 +- .../wells/kernels/SinglePhaseWellKernels.hpp | 13 +- ...rmalCompositionalMultiphaseWellKernels.hpp | 19 +- .../kernels/ThermalSinglePhaseWellKernels.hpp | 14 +- 15 files changed, 308 insertions(+), 403 deletions(-) diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.cpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.cpp index 0914b977dfd..ba8834b02ce 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.cpp @@ -38,6 +38,15 @@ #include "physicsSolvers/fluidFlow/wells/WellFields.hpp" #include "physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWellFields.hpp" #include "physicsSolvers/fluidFlow/wells/WellControls.hpp" + +#include "physicsSolvers/fluidFlow/wells/WellInjectionConstraint.hpp" +#include "physicsSolvers/fluidFlow/wells/WellProductionConstraint.hpp" +#include "physicsSolvers/fluidFlow/wells/WellBHPConstraints.hpp" +#include "physicsSolvers/fluidFlow/wells/WellVolumeRateConstraint.hpp" +#include "physicsSolvers/fluidFlow/wells/WellPhaseVolumeRateConstraint.hpp" +#include "physicsSolvers/fluidFlow/wells/WellMassRateConstraint.hpp" +#include "physicsSolvers/fluidFlow/wells/WellLiquidRateConstraint.hpp" + #include "physicsSolvers/fluidFlow/wells/kernels/CompositionalMultiphaseWellKernels.hpp" #include "physicsSolvers/fluidFlow/wells/kernels/ThermalCompositionalMultiphaseWellKernels.hpp" #include "physicsSolvers/fluidFlow/wells/kernels/PerforationFluxKernels.hpp" @@ -192,7 +201,6 @@ void CompositionalMultiphaseWell::registerWellDataOnMesh( WellElementSubRegion & subRegion.registerField< well::gravityCoefficient >( getName() ); - subRegion.registerField< well::globalCompDensity >( getName() ). reference().resizeDimension< 1 >( m_numComponents ); subRegion.registerField< well::globalCompDensity_n >( getName() ). diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/SinglePhaseWell.cpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/SinglePhaseWell.cpp index a2172b71792..72b95c20ae2 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/SinglePhaseWell.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/SinglePhaseWell.cpp @@ -38,6 +38,15 @@ #include "physicsSolvers/fluidFlow/wells/WellFields.hpp" #include "physicsSolvers/fluidFlow/wells/WellSolverBaseFields.hpp" #include "physicsSolvers/fluidFlow/wells/SinglePhaseWellFields.hpp" + +#include "physicsSolvers/fluidFlow/wells/WellInjectionConstraint.hpp" +#include "physicsSolvers/fluidFlow/wells/WellProductionConstraint.hpp" +#include "physicsSolvers/fluidFlow/wells/WellBHPConstraints.hpp" +#include "physicsSolvers/fluidFlow/wells/WellVolumeRateConstraint.hpp" +#include "physicsSolvers/fluidFlow/wells/WellPhaseVolumeRateConstraint.hpp" +#include "physicsSolvers/fluidFlow/wells/WellMassRateConstraint.hpp" +#include "physicsSolvers/fluidFlow/wells/WellLiquidRateConstraint.hpp" + #include "physicsSolvers/fluidFlow/wells/kernels/SinglePhaseWellKernels.hpp" #include "physicsSolvers/fluidFlow/wells/kernels/ThermalSinglePhaseWellKernels.hpp" #include "physicsSolvers/fluidFlow/wells/kernels/SinglePhasePerforationFluxKernels.hpp" diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/WellBHPConstraints.cpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/WellBHPConstraints.cpp index c6fbca6e40c..ac03db749b4 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/WellBHPConstraints.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/WellBHPConstraints.cpp @@ -84,14 +84,13 @@ bool BHPConstraint< T >::checkViolation( WellConstraintBase const & currentConst } } - template class BHPConstraint< BHPConstraintTypeId::MIN >; template class BHPConstraint< BHPConstraintTypeId::MAX >; + namespace { -typedef BHPConstraint< BHPConstraintTypeId::MIN > MinimumBHPConstraint; -typedef BHPConstraint< BHPConstraintTypeId::MAX > MaximumBHPConstraint; REGISTER_CATALOG_ENTRY( WellConstraintBase, MinimumBHPConstraint, string const &, Group * const ) REGISTER_CATALOG_ENTRY( WellConstraintBase, MaximumBHPConstraint, string const &, Group * const ) } + } //namespace geos diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/WellBHPConstraints.hpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/WellBHPConstraints.hpp index 4ce9312d019..28f356435c3 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/WellBHPConstraints.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/WellBHPConstraints.hpp @@ -24,6 +24,7 @@ #include "dataRepository/Group.hpp" #include "functions/TableFunction.hpp" #include "WellConstraintsBase.hpp" + namespace geos { @@ -108,8 +109,7 @@ class BHPConstraint : public WellConstraintBase static constexpr char const * targetBHPString() { return "targetBHP"; } /// String key for the well reference elevation (for BHP control) static constexpr char const * refElevString() { return "referenceElevation"; } - } - viewKeysWellBHPConstraint; + }; virtual bool checkViolation( WellConstraintBase const & currentConstraint, real64 const & currentTime ) const override; diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/WellConstraintsBase.cpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/WellConstraintsBase.cpp index 9017e5e36bd..ee77440b105 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/WellConstraintsBase.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/WellConstraintsBase.cpp @@ -123,7 +123,7 @@ void WellConstraintBase::postInputInitialization() } -void WellConstraintBase::setNextDtFromTables( real64 const currentTime, real64 & nextDt ) +void WellConstraintBase::setNextDtFromTables( real64 const currentTime, real64 & nextDt ) const { setNextDtFromTable( m_constraintScheduleTable, currentTime, nextDt ); } diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/WellConstraintsBase.hpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/WellConstraintsBase.hpp index 52708d52ba3..0e732d02d5a 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/WellConstraintsBase.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/WellConstraintsBase.hpp @@ -233,7 +233,7 @@ class WellConstraintBase : public dataRepository::Group * @param[in] currentTime the current time * @param[inout] nextDt the time step */ - void setNextDtFromTables( real64 const currentTime, real64 & nextDt ); + void setNextDtFromTables( real64 const currentTime, real64 & nextDt ) const; protected: @@ -247,7 +247,7 @@ class WellConstraintBase : public dataRepository::Group /// Constraint value real64 m_constraintValue; - void setNextDtFromTable( TableFunction const * table, real64 const currentTime, real64 & nextDt ); + static void setNextDtFromTable( TableFunction const * table, real64 const currentTime, real64 & nextDt ); /// Constraint schedule table name string m_constraintScheduleTableName; diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/WellControls.cpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/WellControls.cpp index f4275a6fa53..45aa94574bc 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/WellControls.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/WellControls.cpp @@ -17,8 +17,17 @@ * @file WellControls.cpp */ -#include "LogLevelsInfo.hpp" #include "WellControls.hpp" + +#include "physicsSolvers/fluidFlow/wells/WellInjectionConstraint.hpp" +#include "physicsSolvers/fluidFlow/wells/WellProductionConstraint.hpp" +#include "physicsSolvers/fluidFlow/wells/WellBHPConstraints.hpp" +#include "physicsSolvers/fluidFlow/wells/WellVolumeRateConstraint.hpp" +#include "physicsSolvers/fluidFlow/wells/WellPhaseVolumeRateConstraint.hpp" +#include "physicsSolvers/fluidFlow/wells/WellMassRateConstraint.hpp" +#include "physicsSolvers/fluidFlow/wells/WellLiquidRateConstraint.hpp" + +#include "LogLevelsInfo.hpp" #include "WellConstants.hpp" #include "dataRepository/InputFlags.hpp" #include "functions/FunctionManager.hpp" @@ -41,7 +50,6 @@ WellControls::WellControls( string const & name, Group * const parent ) m_numDofPerResElement( 0 ), m_isThermal( 0 ), m_keepVariablesConstantDuringInitStep( false ), - //m_ratesOutputDir( joinPath( OutputBase::getOutputDirectory(), "_rates" ) ), m_ratesOutputDir( joinPath( OutputBase::getOutputDirectory(), parent->getName() + "_rates" ) ), m_inputControl( Control::UNINITIALIZED ), m_currentControl( Control::UNINITIALIZED ), @@ -58,7 +66,7 @@ WellControls::WellControls( string const & name, Group * const parent ) m_estimateSolution( 0 ), m_enableIsoThermalEstimator( 0 ), /// Nonlinear solver parameters - m_wellNewtonSolver( viewKeyStruct::wellNewtonSolverString(), this ), + m_wellNewtonSolver( groupKeyStruct::wellNewtonSolverString(), this ), m_estimatorDoFManager( name ), m_dofManagerInitialized( false ) { @@ -143,8 +151,7 @@ WellControls::WellControls( string const & name, Group * const parent ) setDescription( "Name of the well status table when the status of the well is a time dependent function. \n" "If the status function evaluates to a positive value at the current time, the well will be open otherwise the well will be shut." ); - - registerGroup( viewKeyStruct::wellNewtonSolverString(), &m_wellNewtonSolver ); + registerGroup( groupKeyStruct::wellNewtonSolverString(), &m_wellNewtonSolver ); addLogLevel< logInfo::WellControl >(); } @@ -153,78 +160,6 @@ WellControls::WellControls( string const & name, Group * const parent ) WellControls::~WellControls() {} -Group * WellControls::createChild( string const & childKey, string const & childName ) -{ - - Group * child = nullptr; - if( childKey == viewKeyStruct::minimumBHPConstraintString() ) - { - BHPConstraint< BHPConstraintTypeId::MIN > & bhpConstraint = registerGroup< BHPConstraint< BHPConstraintTypeId::MIN > >( childName ); - m_minBHPConstraint = &bhpConstraint; - child = &bhpConstraint; - } - else if( childKey == viewKeyStruct::maximumBHPConstraintString() ) - { - BHPConstraint< BHPConstraintTypeId::MAX > & bhpConstraint = registerGroup< BHPConstraint< BHPConstraintTypeId::MAX > >( childName ); - m_maxBHPConstraint = &bhpConstraint; - child = &bhpConstraint; - } - else if( childKey == viewKeyStruct::productionPhaseVolumeRateConstraintString() ) - { - ProductionConstraint< PhaseVolumeRateConstraint > & phaseConstraint = registerGroup< ProductionConstraint< PhaseVolumeRateConstraint > >( childName ); - m_productionRateConstraintList.emplace_back( &phaseConstraint ); - child = &phaseConstraint; - } - else if( childKey == viewKeyStruct::injectionPhaseVolumeRateConstraint() ) - { - InjectionConstraint< PhaseVolumeRateConstraint > & phaseConstraint = registerGroup< InjectionConstraint< PhaseVolumeRateConstraint > >( childName ); - m_injectionRateConstraintList.emplace_back( &phaseConstraint ); - child = &phaseConstraint; - } - else if( childKey == viewKeyStruct::productionVolumeRateConstraint() ) - { - ProductionConstraint< VolumeRateConstraint > & volConstraint = registerGroup< ProductionConstraint< VolumeRateConstraint > >( childName ); - m_productionRateConstraintList.emplace_back( &volConstraint ); - child = &volConstraint; - } - else if( childKey == viewKeyStruct::injectionVolumeRateConstraint() ) - { - InjectionConstraint< VolumeRateConstraint > & volConstraint = registerGroup< InjectionConstraint< VolumeRateConstraint > >( childName ); - m_injectionRateConstraintList.emplace_back( &volConstraint ); - child = &volConstraint; - } - else if( childKey == viewKeyStruct::productionMassRateConstraint() ) - { - ProductionConstraint< MassRateConstraint > & massConstraint = registerGroup< ProductionConstraint< MassRateConstraint > >( childName ); - m_productionRateConstraintList.emplace_back( &massConstraint ); - child = &massConstraint; - - } - else if( childKey == viewKeyStruct::injectionMassRateConstraint() ) - { - InjectionConstraint< MassRateConstraint > & massConstraint = registerGroup< InjectionConstraint< MassRateConstraint > >( childName ); - m_injectionRateConstraintList.emplace_back( &massConstraint ); - child = &massConstraint; - } - else if( childKey == viewKeyStruct::productionLiquidRateConstraint() ) - { - ProductionConstraint< LiquidRateConstraint > & liquidConstraint = registerGroup< ProductionConstraint< LiquidRateConstraint > >( childName ); - m_productionRateConstraintList.emplace_back( &liquidConstraint ); - child = &liquidConstraint; - } - return child; -} - -void WellControls::expandObjectCatalogs() -{ - // During schema generation, register one of each type derived from ConstitutiveBase here - for( auto & catalogIter: WellConstraintBase::getCatalog()) - { - createChild( catalogIter.first, catalogIter.first ); - } - // tjbcreateChild( WellNewtonSolver::catalogName(), WellNewtonSolver::catalogName() ); -} - namespace { @@ -247,6 +182,7 @@ TableFunction * createWellTable( string const & tableName, } } + void WellControls::registerWellDataOnMesh( WellElementSubRegion & subRegion ) { std::string const & regionName = subRegion.getName(); @@ -276,6 +212,7 @@ void WellControls::registerWellDataOnMesh( WellElementSubRegion & subRegion ) } } } + void WellControls::postInputInitialization() { Group::postInputInitialization(); @@ -348,45 +285,89 @@ void WellControls::postInputInitialization() } void WellControls::postRestartInitialization( ) {} + +void WellControls::logConstraint( WellConstraintBase const * constraint, WellElementSubRegion const & region, real64 time, bool isLimiting ) +{ + bool const needsLog = (constraint != nullptr) && (getLogLevel() > 4) && region.isLocallyOwned(); + if( isLimiting ) + { + GEOS_LOG_RANK_IF ( needsLog, + GEOS_FMT( " Well {}: Limiting Constraint {} - BHP {}, Volume rates {}, Total rate {}, Mass rate {}", + region.getName(), constraint->getName(), constraint->bottomHolePressure(), + constraint->phaseVolumeRates(), constraint->totalVolumeRate(), constraint->massRate()) ); + } + else + { + GEOS_LOG_RANK_IF ( needsLog, + GEOS_FMT( " Well {}: Constraint {} - active {}, value {}", + region.getName(), constraint->getName(), constraint->isConstraintActive(), + constraint->getConstraintValue( time ) ) ); + } +} + +//template< typename T > +//using isConstraint = std::disjunction< +// std::is_same< T, MinimumBHPConstraint >, +//std::is_same< T, MaximumBHPConstraint >, +//std::is_same< T, ProductionConstraint< PhaseVolumeRateConstraint > >, +//std::is_same, +//std::is_same< T, ProductionConstraint< MassRateConstraint > >, +//std::is_same< T, InjectionConstraint< PhaseVolumeRateConstraint > >, +//std::is_same< T, InjectionConstraint< VolumeRateConstraint >, +// std::is_same< T, InjectionConstraint< MassRateConstraint > > +// > +//template< typename ConstraintType > +//ConstraintType const * WellControls::getFirstConstraintOfType() const +//{ +// static_assert( isConstraint< ConstraintType >::value, "ConstraintType must be a constraint type" ); +// return nullptr; +//} + void WellControls::setWellStatus( real64 const & currentTime, WellControls::Status status ) { m_wellStatus = status; if( m_wellStatus == WellControls::Status::OPEN ) { - if( isProducer()) + bool hasZeroRate = false; + + if( m_statusTable->evaluate( ¤tTime ) < LvArray::NumericLimits< real64 >::epsilon ) + { + hasZeroRate = true; + } + + if( !hasZeroRate && isProducer()) { - std::vector< WellConstraintBase * > const constraints = getProdRateConstraints(); - for( auto const & constraint : constraints ) + forSubGroups< ProductionConstraint< PhaseVolumeRateConstraint >, + ProductionConstraint< VolumeRateConstraint >, + ProductionConstraint< MassRateConstraint > >( [&] ( auto const & constraint ) { - if( isZero( constraint->getConstraintValue( currentTime ) ) ) + if( isZero( constraint.getConstraintValue( currentTime ) ) ) { - m_wellStatus = WellControls::Status::CLOSED; - m_currentConstraint=nullptr; - break; + hasZeroRate = true; } - } + } ); } - else + else if( !hasZeroRate ) { - std::vector< WellConstraintBase * > const constraints = getInjRateConstraints(); - for( auto const & constraint : constraints ) + forSubGroups< InjectionConstraint< PhaseVolumeRateConstraint >, + InjectionConstraint< VolumeRateConstraint >, + InjectionConstraint< MassRateConstraint > >( [&] ( auto const & constraint ) { - if( isZero( constraint->getConstraintValue( currentTime ) ) ) + if( isZero( constraint.getConstraintValue( currentTime ) ) ) { - m_wellStatus = WellControls::Status::CLOSED; - m_currentConstraint=nullptr; - break; + hasZeroRate = true; } - } + } ); } - if( m_statusTable->evaluate( ¤tTime ) < LvArray::NumericLimits< real64 >::epsilon ) + if( hasZeroRate ) { - m_wellStatus = WellControls::Status::CLOSED; - m_currentConstraint=nullptr; + m_wellStatus = WellControls::Status::CLOSED; + m_currentConstraint = nullptr; } } } + real64 WellControls::setNextDt( real64 const & currentTime, real64 const & currentDt, WellElementSubRegion & subRegion ) @@ -413,6 +394,7 @@ real64 WellControls::setNextDt( real64 const & currentTime, GEOS_LOG_RANK_0( GEOS_FMT( "{}: next time step based on tables coordinates = {}", getName(), nextDt )); return nextDt; } + bool WellControls::isWellOpen() const { return getWellStatus() == WellControls::Status::OPEN; @@ -432,25 +414,26 @@ void WellControls::setNextDtFromTables( real64 const & currentTime, real64 & nex { if( isProducer() ) { - if( getMinBHPConstraint() != nullptr ) - { - getMinBHPConstraint()->setNextDtFromTables( currentTime, nextDt ); - } - for( auto const & constraint : m_productionRateConstraintList ) + forSubGroups< MinimumBHPConstraint, + ProductionConstraint< PhaseVolumeRateConstraint >, + ProductionConstraint< VolumeRateConstraint >, + ProductionConstraint< MassRateConstraint > >( [&] ( auto const & constraint ) { - constraint->setNextDtFromTables( currentTime, nextDt ); - } + constraint.setNextDtFromTables( currentTime, nextDt ); + } ); } else { - getMaxBHPConstraint()->setNextDtFromTables( currentTime, nextDt ); - for( auto const & constraint : m_injectionRateConstraintList ) + forSubGroups< MaximumBHPConstraint, + InjectionConstraint< PhaseVolumeRateConstraint >, + InjectionConstraint< VolumeRateConstraint >, + InjectionConstraint< MassRateConstraint > >( [&] ( auto const & constraint ) { - constraint->setNextDtFromTables( currentTime, nextDt ); - } + constraint.setNextDtFromTables( currentTime, nextDt ); + } ); } - WellControls::setNextDtFromTable( m_statusTable, currentTime, nextDt ); + setNextDtFromTable( m_statusTable, currentTime, nextDt ); } void WellControls::setNextDtFromTable( TableFunction const * table, real64 const currentTime, real64 & nextDt ) @@ -469,18 +452,30 @@ void WellControls::setNextDtFromTable( TableFunction const * table, real64 const real64 WellControls::getTargetBHP( real64 const & targetTime ) const { + WellConstraintBase const * bhpConstraint = nullptr; if( isProducer()) { - return m_minBHPConstraint->getConstraintValue( targetTime ); + forSubGroups< MinimumBHPConstraint >( [&] ( auto const & constraint ) + { + bhpConstraint = &constraint; + } ); + } + else + { + forSubGroups< MaximumBHPConstraint >( [&] ( auto const & constraint ) + { + bhpConstraint = &constraint; + } ); } - return m_maxBHPConstraint->getConstraintValue( targetTime ); + return bhpConstraint->getConstraintValue( targetTime ); } - real64 WellControls::getInjectionTemperature() const { real64 injectionTemperature = 0.0; - this->forInjectionConstraints< InjectionConstraint< PhaseVolumeRateConstraint >, InjectionConstraint< VolumeRateConstraint >, InjectionConstraint< MassRateConstraint > >( [&] ( auto & constraint ) + forSubGroups< InjectionConstraint< PhaseVolumeRateConstraint >, + InjectionConstraint< VolumeRateConstraint >, + InjectionConstraint< MassRateConstraint > >( [&] ( auto & constraint ) { if( constraint.isConstraintActive()) { @@ -491,11 +486,12 @@ real64 WellControls::getInjectionTemperature() const return injectionTemperature; } - arrayView1d< real64 const > WellControls::getInjectionStream() const { arrayView1d< real64 const > injectionStream; - forInjectionConstraints< InjectionConstraint< PhaseVolumeRateConstraint >, InjectionConstraint< VolumeRateConstraint >, InjectionConstraint< MassRateConstraint > >( [&] ( auto & constraint ) + forSubGroups< InjectionConstraint< PhaseVolumeRateConstraint >, + InjectionConstraint< VolumeRateConstraint >, + InjectionConstraint< MassRateConstraint > >( [&] ( auto & constraint ) { if( constraint.isConstraintActive() ) { @@ -513,7 +509,7 @@ integer WellControls::getConstraintPhaseIndex() const if( isProducer() ) { - forProductionConstraints< ProductionConstraint< PhaseVolumeRateConstraint > >( [&] ( auto & constraint ) + forSubGroups< ProductionConstraint< PhaseVolumeRateConstraint > >( [&] ( auto & constraint ) { if( constraint.isConstraintActive() ) { @@ -523,7 +519,7 @@ integer WellControls::getConstraintPhaseIndex() const } else { - forInjectionConstraints< InjectionConstraint< PhaseVolumeRateConstraint > >( [&] ( auto & constraint ) + forSubGroups< InjectionConstraint< PhaseVolumeRateConstraint > >( [&] ( auto & constraint ) { if( constraint.isConstraintActive() ) { @@ -537,11 +533,22 @@ integer WellControls::getConstraintPhaseIndex() const real64 WellControls::getReferenceElevation() const { + real64 referenceElevation = 0.0; if( isProducer () ) { - return getMinBHPConstraint()->getReferenceElevation(); + forSubGroups< MinimumBHPConstraint >( [&] ( auto & constraint ) + { + referenceElevation = constraint.getReferenceElevation(); + } ); } - return getMaxBHPConstraint()->getReferenceElevation(); + else + { + forSubGroups< MaximumBHPConstraint >( [&] ( auto & constraint ) + { + referenceElevation = constraint.getReferenceElevation(); + } ); + } + return referenceElevation; } void WellControls::implicitStepSetup( real64 const & time_n, @@ -554,6 +561,7 @@ void WellControls::implicitStepSetup( real64 const & time_n, // Set perforation status setPerforationStatus( time_n, subRegion ); } + void WellControls::setPerforationStatus( real64 const & time_n, WellElementSubRegion & subRegion ) { FunctionManager & functionManager = FunctionManager::getInstance(); @@ -718,8 +726,6 @@ void WellControls::selectWellConstraint( real64 const & time_n, MpiWrapper::broadcast( wellControl, owner ); setControl( wellControl ); } - - } bool WellControls::evaluateConstraints( real64 const & time_n, @@ -727,30 +733,35 @@ bool WellControls::evaluateConstraints( real64 const & time_n, { // create list of all constraints to process - std::vector< WellConstraintBase * > constraintList; + stdVector< WellConstraintBase * > constraintList; if( isProducer() ) { - constraintList = getProdRateConstraints(); - // Solve minimum bhp constraint first - constraintList.insert( constraintList.begin(), getMinBHPConstraint() ); + forSubGroups< MinimumBHPConstraint, + ProductionConstraint< PhaseVolumeRateConstraint >, + ProductionConstraint< VolumeRateConstraint >, + ProductionConstraint< MassRateConstraint > >( [&] ( auto & constraint ) + { + constraintList.emplace_back( &constraint ); + } ); } else { - constraintList = getInjRateConstraints(); - // Solve maximum bhp constraint first; - constraintList.insert( constraintList.begin(), getMaxBHPConstraint() ); + forSubGroups< MaximumBHPConstraint, + InjectionConstraint< PhaseVolumeRateConstraint >, + InjectionConstraint< VolumeRateConstraint >, + InjectionConstraint< MassRateConstraint > >( [&] ( auto & constraint ) + { + constraintList.emplace_back( &constraint ); + } ); } // Get current constraint - WellConstraintBase * limitingConstraint = nullptr; + WellConstraintBase * limitingConstraint = nullptr; for( auto & constraint : constraintList ) { if( constraint->getName() == getCurrentConstraint()->getName()) { - limitingConstraint = constraint; - GEOS_LOG_RANK_IF ( getLogLevel() > 4 && subRegion.isLocallyOwned(), - " Well " << subRegion.getName() << " Limiting Constraint " << limitingConstraint->getName() << " " << limitingConstraint->bottomHolePressure() << " " << - limitingConstraint->phaseVolumeRates() << " " << - limitingConstraint->totalVolumeRate() << " " << limitingConstraint->massRate()); + limitingConstraint = constraint; + logConstraint( limitingConstraint, subRegion, time_n, true ); } } constraintList.erase( std::find( constraintList.begin(), constraintList.end(), limitingConstraint ) ); @@ -771,10 +782,7 @@ bool WellControls::evaluateConstraints( real64 const & time_n, } } } - GEOS_LOG_RANK_IF ( getLogLevel() > 4 && subRegion.isLocallyOwned(), - " Well " << subRegion.getName() << " Limiting Constraint " << limitingConstraint->getName() << " " << limitingConstraint->bottomHolePressure() << " " << - limitingConstraint->phaseVolumeRates() << " " << - limitingConstraint->totalVolumeRate() << " " << limitingConstraint->massRate()); + logConstraint( limitingConstraint, subRegion, time_n, true ); return true; } @@ -789,16 +797,19 @@ bool WellControls::evaluateConstraints( real64 const & time_n, WellElementSubRegion & subRegion, DofManager const & dofManager ) { - - // create list of all constraints to solve // note that initializeWells sets the initial constraint // tjb reactive control schema to allow use to set if needed - std::vector< WellConstraintBase * > constraintList; - WellConstraintBase * limitingConstraint = getCurrentConstraint(); + stdVector< WellConstraintBase * > constraintList; + WellConstraintBase * limitingConstraint = getCurrentConstraint(); if( isProducer() ) { - constraintList = getProdRateConstraints(); + forSubGroups< ProductionConstraint< PhaseVolumeRateConstraint >, + ProductionConstraint< VolumeRateConstraint >, + ProductionConstraint< MassRateConstraint > >( [&] ( auto & constraint ) + { + constraintList.emplace_back( &constraint ); + } ); if( limitingConstraint->getControl() != ConstraintTypeId::BHP ) { { // remove from list and add BHP constraint @@ -807,38 +818,33 @@ bool WellControls::evaluateConstraints( real64 const & time_n, { constraintList.erase( it ); } - if( getMinBHPConstraint()->isConstraintActive() ) + forSubGroups< MinimumBHPConstraint >( [&] ( auto & constraint ) { - constraintList.push_back( getMinBHPConstraint() ); - } + constraintList.emplace_back( &constraint ); + } ); constraintList.insert( constraintList.begin(), limitingConstraint ); } } // Solve minimum bhp constraint first - if( false && getMinBHPConstraint()->isConstraintActive() ) - { - // this is related to WHP option which introduces a new BHP constraint - limitingConstraint = getMinBHPConstraint(); - } - else if( limitingConstraint == nullptr ) + if( limitingConstraint == nullptr ) { limitingConstraint = constraintList[0]; } } else { - constraintList = getInjRateConstraints(); - // remove the limiting constraint from the list if present + forSubGroups< InjectionConstraint< PhaseVolumeRateConstraint >, + InjectionConstraint< VolumeRateConstraint >, + InjectionConstraint< MassRateConstraint > >( [&] ( auto & constraint ) { - if( limitingConstraint->getControl() != ConstraintTypeId::BHP ) - { // remove from list and add BHP constraint - //auto it = std::find( constraintList.begin(), constraintList.end(), limitingConstraint ); - //if( it != constraintList.end() ) - //{ - // constraintList.erase( it ); - //} - constraintList.push_back( getMaxBHPConstraint() ); - } + constraintList.emplace_back( &constraint ); + } ); + if( limitingConstraint->getControl() != ConstraintTypeId::BHP ) + { + forSubGroups< MaximumBHPConstraint >( [&] ( auto & constraint ) + { + constraintList.emplace_back( &constraint ); + } ); } } if( isoThermalEstimatorEnabled() ) @@ -870,7 +876,7 @@ bool WellControls::evaluateConstraints( real64 const & time_n, GEOS_LOG_RANK_IF ( getLogLevel() > 4 && subRegion.isLocallyOwned(), " Well " << subRegion.getName() << " Constraint " << constraint->getName() << " active " << constraint->isConstraintActive() << " value " << constraint->getConstraintValue( time_n ) ); - if( constraint->isConstraintActive() && constraint->checkViolation( *limitingConstraint, time_n )) + if( constraint->isConstraintActive() && constraint->checkViolation( *limitingConstraint, time_n )) { limitingConstraint=constraint; setControl( static_cast< WellControls::Control >(constraint->getControl()) ); // tjb old diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/WellControls.hpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/WellControls.hpp index 3efe3dfcaee..3a8343917c6 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/WellControls.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/WellControls.hpp @@ -17,22 +17,22 @@ * @file WellControls.hpp */ - #ifndef GEOS_PHYSICSSOLVERS_FLUIDFLOW_WELLS_WELLCONTROLS_HPP #define GEOS_PHYSICSSOLVERS_FLUIDFLOW_WELLS_WELLCONTROLS_HPP + #include "physicsSolvers/PhysicsSolverBase.hpp" #include "common/format/EnumStrings.hpp" #include "dataRepository/Group.hpp" #include "functions/TableFunction.hpp" #include "constitutive/fluid/multifluid/MultiFluidBase.hpp" -#include "physicsSolvers/fluidFlow/wells/WellInjectionConstraint.hpp" -#include "physicsSolvers/fluidFlow/wells/WellProductionConstraint.hpp" -#include "physicsSolvers/fluidFlow/wells/WellBHPConstraints.hpp" -#include "physicsSolvers/fluidFlow/wells/WellVolumeRateConstraint.hpp" -#include "physicsSolvers/fluidFlow/wells/WellPhaseVolumeRateConstraint.hpp" -#include "physicsSolvers/fluidFlow/wells/WellMassRateConstraint.hpp" -#include "physicsSolvers/fluidFlow/wells/WellLiquidRateConstraint.hpp" +//#include "physicsSolvers/fluidFlow/wells/WellInjectionConstraint.hpp" +//#include "physicsSolvers/fluidFlow/wells/WellProductionConstraint.hpp" +//#include "physicsSolvers/fluidFlow/wells/WellBHPConstraints.hpp" +//#include "physicsSolvers/fluidFlow/wells/WellVolumeRateConstraint.hpp" +//#include "physicsSolvers/fluidFlow/wells/WellPhaseVolumeRateConstraint.hpp" +//#include "physicsSolvers/fluidFlow/wells/WellMassRateConstraint.hpp" +//#include "physicsSolvers/fluidFlow/wells/WellLiquidRateConstraint.hpp" #include "constitutive/fluid/multifluid/MultiFluidBase.hpp" #include "constitutive/fluid/singlefluid/SingleFluidBase.hpp" #include "physicsSolvers/fluidFlow/wells/WellNewtonSolver.hpp" @@ -48,12 +48,13 @@ static constexpr auto wellControls = "WellControls"; } class ElementsReporterBuffer; +class WellConstraintBase; /** * @class WellControls * @brief This class describes the controls used to operate a well. */ -class WellControls : public dataRepository::Group //public PhysicsSolverBase +class WellControls : public dataRepository::Group { public: @@ -138,16 +139,6 @@ class WellControls : public dataRepository::Group //public PhysicsSolverBase /// String used to form the solverName used to register single-physics solvers in CoupledSolver static string coupledSolverAttributePrefix() { return "well"; } - /** - * @brief Create a new geometric object (box, plane, etc) as a child of this group. - * @param childKey the catalog key of the new geometric object to create - * @param childName the name of the new geometric object in the repository - * @return the group child - */ - virtual Group * createChild( string const & childKey, string const & childName ) override; - /// Expand catalog for schema generation - - virtual void expandObjectCatalogs() override; virtual void registerWellDataOnMesh( WellElementSubRegion & subRegion ); virtual void setConstitutiveNames( ElementSubRegionBase & subRegion ) const = 0; @@ -385,6 +376,7 @@ class WellControls : public dataRepository::Group //public PhysicsSolverBase virtual void printRates( real64 const & time_n, real64 const & dt, WellElementSubRegion const & subRegion ) = 0; +#if 0 /**@}*/ /** * @brief Apply a given functor to a container if the container can be @@ -492,6 +484,7 @@ class WellControls : public dataRepository::Group //public PhysicsSolverBase } ); } } +#endif /** * @name Getters / Setters */ @@ -689,6 +682,26 @@ class WellControls : public dataRepository::Group //public PhysicsSolverBase WellConstraintBase * getCurrentConstraint() { return m_currentConstraint; } WellConstraintBase const * getCurrentConstraint() const { return m_currentConstraint; } + /** + * @brief Get the BHP constraint. + * @details Returns the first BHP constraint attached to the well controls. This will be the minimum BHP + * constraint for a producer well and the maximum BHP constraint for an injector well. If there is no BHP + * constraint assigned to the well controls this will return a null value. + * @return The constraint object of one exists nullptr otherwise + */ + //WellConstraintBase const * getBHPConstraint() const; + //WellConstraintBase * getBHPConstraint(); + + /** + * @brief Get the rate constraint. + * @details Returns the first rate constraint attached to the well controls. This will be a prodcution rate + * constraint for a producer well and an injection rate constraint for an injector well. This will return the + * first constraint assigned regardless of type (phase volume, liquid volume or mass). If there is no rate + * constraint assigned to the well controls this will return a null value. + * @return The constraint object of one exists nullptr otherwise + */ + WellConstraintBase const * getRateConstraint() const; + //WellConstraintBase * getRateConstraint(); /** * @brief Getter for the flag to enable crossflow @@ -773,12 +786,8 @@ class WellControls : public dataRepository::Group //public PhysicsSolverBase * @return a Status */ WellControls::Status getWellStatus () const { return m_wellStatus; } - - - ///@} - virtual string wellElementDofName() const = 0; virtual string resElementDofName() const = 0; @@ -799,6 +808,7 @@ class WellControls : public dataRepository::Group //public PhysicsSolverBase virtual localIndex numFluidComponents() const = 0; virtual localIndex numFluidPhases() const = 0; + /** * @brief Struct to serve as a container for variable strings and keys. * @struct viewKeyStruct @@ -810,7 +820,7 @@ class WellControls : public dataRepository::Group //public PhysicsSolverBase /// String key for the write CSV flag static constexpr char const * writeCSVFlagString() { return "writeCSV"; } static constexpr char const * timeStepFromTablesFlagString() { return "timeStepFromTables"; } -/// @return string for the targetRegions wrapper + /// String for the targetRegions wrapper static constexpr char const * targetRegionsString() { return "targetRegions"; } /// String key for the well reference elevation (for BHP control) @@ -839,6 +849,30 @@ class WellControls : public dataRepository::Group //public PhysicsSolverBase /// string key for the initial pressure coefficient static constexpr char const * initialPressureCoefficientString() { return "initialPressureCoefficient"; } + /// string key for the esitmate well solution flag + static constexpr char const * estimateWellSolutionString() { return "estimateWellSolution"; } + /// string key for the enable iso thermal estimator flag + static constexpr char const * enableIsoThermalEstimatorString() { return "enableIsoThermalEstimator"; } + + // control data (not registered on the mesh) + static constexpr char const * massDensityString() { return "massDensity";} + + static constexpr char const * currentBHPString() { return "currentBHP"; } + + static constexpr char const * currentPhaseVolRateString() { return "currentPhaseVolumetricRate"; } + static constexpr char const * currentVolRateString() { return "currentVolRate"; } + + static constexpr char const * currentTotalVolRateString() { return "currentTotalVolumetricRate"; } + + static constexpr char const * currentMassRateString() { return "currentMassRate"; } + }; + + /** + * @brief Structure to hold scoped key names + */ + struct groupKeyStruct + { +#if 0 /// string key for the minimum BHP presssure for a producer static constexpr char const * minimumBHPConstraintString() { return "MinimumBHPConstraint"; } /// string key for the maximum BHP presssure for a injection @@ -857,30 +891,11 @@ class WellControls : public dataRepository::Group //public PhysicsSolverBase static constexpr char const * injectionMassRateConstraint() { return "InjectionMassRateConstraint"; } /// string key for the liquid rate for a producer static constexpr char const * productionLiquidRateConstraint() { return "ProductionLiquidRateConstraint"; } - /// string key for the minimum BHP presssure for a producer +#endif + /// string key for the well Newton solver static constexpr char const * wellNewtonSolverString() { return "WellNewtonSolver"; } + }; - /// string key for the esitmate well solution flag - static constexpr char const * estimateWellSolutionString() { return "estimateWellSolution"; } - /// string key for the enable iso thermal estimator flag - static constexpr char const * enableIsoThermalEstimatorString() { return "enableIsoThermalEstimator"; } - - // control data (not registered on the mesh) - - static constexpr char const * massDensityString() { return "massDensity";} - - static constexpr char const * currentBHPString() { return "currentBHP"; } - - static constexpr char const * currentPhaseVolRateString() { return "currentPhaseVolumetricRate"; } - static constexpr char const * currentVolRateString() { return "currentVolRate"; } - - static constexpr char const * currentTotalVolRateString() { return "currentTotalVolumetricRate"; } - - static constexpr char const * currentMassRateString() { return "currentMassRate"; } - - } - /// ViewKey struct for the WellControls class - viewKeysWellControls; void setPerforationStatus( real64 const & time_n, WellElementSubRegion & subRegion ); void setGravCoef( WellElementSubRegion & subRegion, R1Tensor const & gravVector ); /** @@ -899,22 +914,23 @@ class WellControls : public dataRepository::Group //public PhysicsSolverBase */ template< typename ConstraintType > void createConstraint ( string const & constraintName ); - +#if 0 /** * @brief Getters for constraints */ - BHPConstraint< BHPConstraintTypeId::MIN > * getMinBHPConstraint() { return m_minBHPConstraint; }; - BHPConstraint< BHPConstraintTypeId::MIN > * getMinBHPConstraint() const { return m_minBHPConstraint; }; - BHPConstraint< BHPConstraintTypeId::MAX > * getMaxBHPConstraint() { return m_maxBHPConstraint; }; - BHPConstraint< BHPConstraintTypeId::MAX > * getMaxBHPConstraint() const { return m_maxBHPConstraint; }; + BHPConstraint< BHPConstraintTypeId::MIN > * getMinBHPConstraint() { return nullptr; }; + BHPConstraint< BHPConstraintTypeId::MIN > * getMinBHPConstraint() const { return nullptr; }; + BHPConstraint< BHPConstraintTypeId::MAX > * getMaxBHPConstraint() { return nullptr; }; + BHPConstraint< BHPConstraintTypeId::MAX > * getMaxBHPConstraint() const { return nullptr; }; /** * @brief Getters for constraint lists */ - std::vector< WellConstraintBase * > getProdRateConstraints() { return m_productionRateConstraintList; }; - std::vector< WellConstraintBase * > getProdRateConstraints() const { return m_productionRateConstraintList; }; - std::vector< WellConstraintBase * > getInjRateConstraints() { return m_injectionRateConstraintList; } - std::vector< WellConstraintBase * > getInjRateConstraints() const { return m_injectionRateConstraintList; } + std::vector< WellConstraintBase * > getProdRateConstraints() { return {}; }; + std::vector< WellConstraintBase * > getProdRateConstraints() const { return {}; }; + std::vector< WellConstraintBase * > getInjRateConstraints() { return {}; } + std::vector< WellConstraintBase * > getInjRateConstraints() const { return {}; } +#endif /** * @brief Set thermal effects enable @@ -948,8 +964,10 @@ class WellControls : public dataRepository::Group //public PhysicsSolverBase DofManager const & dofManager ); protected: - virtual void postRestartInitialization( )override; + + void logConstraint( WellConstraintBase const * constraint, WellElementSubRegion const & region, real64 time, bool isLimiting = false ); + private: /// List of names of regions the solver will be applied to string_array m_targetRegionNames; @@ -987,6 +1005,7 @@ class WellControls : public dataRepository::Group //public PhysicsSolverBase // flag to enable time step selection base on rates/bhp tables coordinates integer m_timeStepFromTables; + /// Reference elevation real64 m_refElevation; @@ -1028,7 +1047,7 @@ class WellControls : public dataRepository::Group //public PhysicsSolverBase // Current constrint WellConstraintBase * m_currentConstraint; - +#if 0 // Minimum and maximum BHP and WHP constraints BHPConstraint< BHPConstraintTypeId::MIN > * m_minBHPConstraint; BHPConstraint< BHPConstraintTypeId::MAX > * m_maxBHPConstraint; @@ -1036,6 +1055,7 @@ class WellControls : public dataRepository::Group //public PhysicsSolverBase // Lists of rate constraints std::vector< WellConstraintBase * > m_productionRateConstraintList; std::vector< WellConstraintBase * > m_injectionRateConstraintList; +#endif /// Well status WellControls::Status m_wellStatus; diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/WellManager.hpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/WellManager.hpp index db764b24072..9a3a2e53260 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/WellManager.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/WellManager.hpp @@ -407,109 +407,6 @@ class WellManager : public PhysicsSolverBase * @param domain the physical domain object */ void chopNegativeDensities( DomainPartition & domain ); -#if 0 - /** - * @brief calculate the norm of the global system residual - * @param time the time at the beginning of the step - * @param dt the desired timestep - * @param domain the domain partition - * @param dofManager degree-of-freedom manager associated with the linear system - * @param localRhs the system right-hand side vector - * @return norm of the residual - * - * This function returns the norm of global residual vector, which is suitable for comparison with - * a tolerance. - */ - virtual real64 - calculateResidualNorm( real64 const & time, - real64 const & dt, - DomainPartition const & domain, - DofManager const & dofManager, - arrayView1d< real64 const > const & localRhs ); - /** - * @brief Function to check system solution for physical consistency and constraint violation - * @param domain the domain partition - * @param dofManager degree-of-freedom manager associated with the linear system - * @param localSolution the solution vector - * @param scalingFactor factor to scale the solution prior to application - * @return true if solution can be safely applied without violating physical constraints, false otherwise - * - * @note This function must be overridden in the derived physics solver in order to use an implict - * solution method such as LinearImplicitStep() or NonlinearImplicitStep(). - * - */ - virtual bool - checkSystemSolution( DomainPartition & domain, - DofManager const & dofManager, - arrayView1d< real64 const > const & localSolution, - real64 const scalingFactor ); - - /** - * @brief Function to determine if the solution vector should be scaled back in order to maintain a known constraint. - * @param[in] domain The domain partition. - * @param[in] dofManager degree-of-freedom manager associated with the linear system - * @param[in] localSolution the solution vector - * @return The factor that should be used to scale the solution vector values when they are being applied. - */ - virtual real64 - scalingForSystemSolution( DomainPartition & domain, - DofManager const & dofManager, - arrayView1d< real64 const > const & localSolution ); - - /** - * @brief Function to apply the solution vector to the state - * @param dofManager degree-of-freedom manager associated with the linear system - * @param localSolution the solution vector - * @param scalingFactor factor to scale the solution prior to application - * @param dt the timestep - * @param domain the domain partition - * - * This function performs 2 operations: - * 1) extract the solution vector for the "blockSystem" parameter, and applies the - * contents of the solution vector to the primary variable field data, - * 2) perform a synchronization of the primary field variable such that all ghosts are updated, - * - * The "scalingFactor" parameter allows for the scaled application of the solution vector. For - * instance, a line search may apply a negative scaling factor to remove part of the previously - * applied solution. - * - * @note This function must be overridden in the derived physics solver in order to use an implict - * solution method such as LinearImplicitStep() or NonlinearImplicitStep(). - * - */ - virtual void - applySystemSolution( DofManager const & dofManager, - arrayView1d< real64 const > const & localSolution, - real64 const scalingFactor, - real64 const dt, - DomainPartition & domain ); - - - /** - * @brief Recompute all dependent quantities from primary variables (including constitutive models) - * @param domain the domain containing the mesh and fields - */ - virtual void updateState( DomainPartition & domain ); - - /** - * @brief perform cleanup for implicit timestep - * @param time the time at the beginning of the step - * @param dt the desired timestep - * @param domain the domain partition - * - * This function performs whatever tasks are required to complete an implicit timestep. For - * example, the acceptance of the solution will occur during this step, and deallocation of - * temporaries will be be performed in this function. - * - * @note This function must be overridden in the derived physics solver in order to use an implict - * solution method such as LinearImplicitStep() or NonlinearImplicitStep(). - */ - virtual void - implicitStepComplete( real64 const & time, - real64 const & dt, - DomainPartition & domain ) override; - -#endif /** * @brief Utility function to keep the well variables during a time step (used in @@ -523,15 +420,12 @@ class WellManager : public PhysicsSolverBase protected: - //virtual void postInputInitialization() override; - virtual void initializePostSubGroups() override; virtual void initializePostInitialConditionsPreSubGroups() override; virtual void postRestartInitialization() override final; - private: /// name of the flow solver diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/kernels/CompositionalMultiphaseWellKernels.cpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/kernels/CompositionalMultiphaseWellKernels.cpp index 010a9db3e7f..69ad0227e06 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/kernels/CompositionalMultiphaseWellKernels.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/kernels/CompositionalMultiphaseWellKernels.cpp @@ -466,12 +466,11 @@ RateInitializationKernel:: { // this assumes phase control present integer const targetPhaseIndex = wellControls.getConstraintPhaseIndex(); - std::vector< WellConstraintBase * > const constraints = wellControls.getProdRateConstraints(); + auto const * rateConstraint = wellControls.getRateConstraint(); // Use first rate constraint to set initial connection rates - real64 const rateVal = constraints[0]->getConstraintValue( time ); + real64 const rateVal = rateConstraint->getConstraintValue( time ); forAll< parallelDevicePolicy<> >( subRegionSize, [=] GEOS_HOST_DEVICE ( localIndex const iwelem ) { - connRate[iwelem] = LvArray::math::max( 0.1 * rateVal * phaseDens[iwelem][0][targetPhaseIndex], -1e3 ); } ); } @@ -507,9 +506,9 @@ RateInitializationKernel:: } else if( controlType == ConstraintTypeId::BHP ) { - std::vector< WellConstraintBase * > const constraints = wellControls.getInjRateConstraints(); + auto const * rateConstraint = wellControls.getRateConstraint(); // Use first rate constraint to set initial connection rates - real64 const rateVal = constraints[0]->getConstraintValue( time ); + real64 const rateVal = rateConstraint->getConstraintValue( time ); forAll< parallelDevicePolicy<> >( subRegionSize, [=] GEOS_HOST_DEVICE ( localIndex const iwelem ) { connRate[iwelem] = LvArray::math::min( 0.1 * rateVal * totalDens[iwelem][0], 1e3 ); diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/kernels/CompositionalMultiphaseWellKernels.hpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/kernels/CompositionalMultiphaseWellKernels.hpp index c81c2582537..2da2725cf39 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/kernels/CompositionalMultiphaseWellKernels.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/kernels/CompositionalMultiphaseWellKernels.hpp @@ -519,36 +519,28 @@ class ResidualNormKernel : public physicsSolverBaseKernels::ResidualNormKernelBa m_iwelemControl( subRegion.getTopWellElementIndex() ), m_isProducer( wellControls.isProducer() ), m_currentControl( wellControls.getControl() ), + m_targetBHP( wellControls.getTargetBHP( time ) ), m_volume( subRegion.getElementVolume() ), m_phaseDens_n( fluid.phaseDensity_n() ), m_totalDens_n( fluid.totalDensity_n() ) { + // tjbNote this assumes that there is only one rate constraint + // This is a normalizer for the balance equations. The normalizaer should be the current rate not the constraint value!! + // This is one of the reasons for restricting constraint type for a production well + // another pr will remove fix this (so the cause for difference results is isolated to one change) + auto const * rateConstraint = wellControls.getRateConstraint(); + if( rateConstraint != nullptr ) + { + m_constraintValue = rateConstraint->getConstraintValue( time ); + } if( m_isProducer ) { - if( wellControls.getMinBHPConstraint()->isConstraintActive()) - m_targetBHP = wellControls.getMinBHPConstraint()->getConstraintValue( time ); - // Note this assumes that there is only one rate constraint - // This is a normalizer for the balance equations. The normalizaer should be the current rate not the constraint value!! - // This is one of the reasons for restricting constraint type for a production well - // another pr will remove fix this (so the cause for difference results is isolated to one change) - m_targetPhaseIndex = wellControls.getConstraintPhaseIndex( ); - m_constraintValue = wellControls.getProdRateConstraints()[0]->getConstraintValue( time ); - + m_targetPhaseIndex = wellControls.getConstraintPhaseIndex(); } else { - m_targetBHP = wellControls.getMaxBHPConstraint()->getConstraintValue( time ); - - // tjbNote this assumes that there is only one rate constraint - // This is a normalizer for the balance equations. The normalizaer should be the current rate not the constraint value!! - // This is one of the reasons for restricting constraint type for a production well - // another pr will remove fix this (so the cause for difference results is isolated to one change) m_targetPhaseIndex = -1; - m_constraintValue = wellControls.getInjRateConstraints()[0]->getConstraintValue( time ); - } - - } GEOS_HOST_DEVICE diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/kernels/SinglePhaseWellKernels.cpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/kernels/SinglePhaseWellKernels.cpp index 58a7f34ccf7..47b18cf5657 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/kernels/SinglePhaseWellKernels.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/kernels/SinglePhaseWellKernels.cpp @@ -550,19 +550,9 @@ RateInitializationKernel:: WellControls::Control const control = wellControls.getControl(); bool const isProducer = wellControls.isProducer(); - real64 constraintVal; - if( isProducer ) - { - std::vector< WellConstraintBase * > const constraints = wellControls.getProdRateConstraints(); - // Use first rate constraint to set initial connection rates - constraintVal = constraints[0]->getConstraintValue( currentTime ); - } - else - { - std::vector< WellConstraintBase * > const constraints = wellControls.getInjRateConstraints(); - // Use first rate constraint to set initial connection rates - constraintVal = constraints[0]->getConstraintValue( currentTime );; - } + auto const * rateConstraint = wellControls.getRateConstraint(); + real64 const constraintVal = rateConstraint->getConstraintValue( currentTime ); + // Estimate the connection rates forAll< parallelDevicePolicy<> >( subRegionSize, [=] GEOS_HOST_DEVICE ( localIndex const iwelem ) { @@ -579,6 +569,10 @@ RateInitializationKernel:: connRate[iwelem] = LvArray::math::min( 0.1 * constraintVal * wellElemDens[iwelem][0], 1e3 ); } } + else if( control == WellControls::Control::MASSRATE ) + { + connRate[iwelem] = constraintVal; + } else { connRate[iwelem] = constraintVal * wellElemDens[iwelem][0]; diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/kernels/SinglePhaseWellKernels.hpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/kernels/SinglePhaseWellKernels.hpp index 8365d25c946..3de515df5fc 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/kernels/SinglePhaseWellKernels.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/kernels/SinglePhaseWellKernels.hpp @@ -35,6 +35,7 @@ #include "physicsSolvers/fluidFlow/wells/WellControls.hpp" #include "physicsSolvers/fluidFlow/wells/WellFields.hpp" #include "physicsSolvers/fluidFlow/wells/SinglePhaseWellFields.hpp" +#include "physicsSolvers/fluidFlow/wells/WellConstraintsBase.hpp" #include "physicsSolvers/KernelLaunchSelectors.hpp" @@ -812,18 +813,10 @@ class ResidualNormKernel : public physicsSolverBaseKernels::ResidualNormKernelBa m_iwelemControl( subRegion.getTopWellElementIndex() ), m_currentControl( wellControls.getControl() ), m_constraintValue ( wellControls.getCurrentConstraint()->getConstraintValue( time )), + m_targetBHP( wellControls.getTargetBHP( time ) ), m_volume( subRegion.getElementVolume() ), m_density_n( fluid.density_n() ) - { - if( wellControls.isProducer() ) - { - m_targetBHP = wellControls.getMinBHPConstraint()->getConstraintValue( time ); - } - else - { - m_targetBHP = wellControls.getMaxBHPConstraint()->getConstraintValue( time ); - } - } + {} GEOS_HOST_DEVICE virtual void computeLinf( localIndex const iwelem, diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/kernels/ThermalCompositionalMultiphaseWellKernels.hpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/kernels/ThermalCompositionalMultiphaseWellKernels.hpp index a8b5b91032f..1c6d39c290e 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/kernels/ThermalCompositionalMultiphaseWellKernels.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/kernels/ThermalCompositionalMultiphaseWellKernels.hpp @@ -181,33 +181,34 @@ class ResidualNormKernel : public physicsSolverBaseKernels::ResidualNormKernelBa m_iwelemControl( subRegion.getTopWellElementIndex() ), m_isProducer( wellControls.isProducer() ), m_currentControl( wellControls.getControl() ), + m_targetBHP( wellControls.getTargetBHP( time ) ), m_volume( subRegion.getElementVolume() ), m_phaseDens_n( fluid.phaseDensity_n() ), m_totalDens_n( fluid.totalDensity_n() ), m_phaseVolFraction_n( subRegion.getField< fields::well::phaseVolumeFraction_n >()), m_phaseInternalEnergy_n( fluid.phaseInternalEnergy_n() ) { + auto const * rateConstraint = wellControls.getRateConstraint(); if( m_isProducer ) { // tjb This needs to be fixed should use current constraint rate for normalization - if( wellControls.getMinBHPConstraint()->isConstraintActive() ) - { - m_targetBHP = wellControls.getMinBHPConstraint()->getConstraintValue( time ); - } if( m_currentControl == WellControls::Control::PHASEVOLRATE ) { m_targetPhaseIndex = wellControls.getConstraintPhaseIndex(); - m_constraintValue = wellControls.getProdRateConstraints()[0]->getConstraintValue( time ); + if( rateConstraint != nullptr ) + { + m_constraintValue = rateConstraint->getConstraintValue( time ); + } } } else { - m_targetBHP = wellControls.getMaxBHPConstraint()->getConstraintValue( time ); m_targetPhaseIndex = -1; - m_constraintValue = wellControls.getInjRateConstraints()[0]->getConstraintValue( time ); - + if( rateConstraint != nullptr ) + { + m_constraintValue = rateConstraint->getConstraintValue( time ); + } } - } diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/kernels/ThermalSinglePhaseWellKernels.hpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/kernels/ThermalSinglePhaseWellKernels.hpp index b19ca90845d..420b9f3b76a 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/kernels/ThermalSinglePhaseWellKernels.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/kernels/ThermalSinglePhaseWellKernels.hpp @@ -523,21 +523,11 @@ class ResidualNormKernel : public physicsSolverBaseKernels::ResidualNormKernelBa m_isProducer( wellControls.isProducer() ), m_currentControl( wellControls.getControl() ), m_constraintValue ( wellControls.getCurrentConstraint()->getConstraintValue( time )), + m_targetBHP( wellControls.getTargetBHP( time ) ), m_volume( subRegion.getElementVolume() ), m_density_n( fluid.density_n() ), - m_internalEnergy_n( fluid.internalEnergy_n() ) - { - if( wellControls.isProducer() ) - { - m_targetBHP = wellControls.getMinBHPConstraint()->getConstraintValue( time ); - } - else - { - m_targetBHP = wellControls.getMaxBHPConstraint()->getConstraintValue( time ); - } - } - + {} GEOS_HOST_DEVICE void computeMassEnergyNormalizers( localIndex const iwelem, From a0e310897997021b9bf836313a17c00befd05874 Mon Sep 17 00:00:00 2001 From: dkachuma Date: Thu, 21 May 2026 13:58:05 -0500 Subject: [PATCH 02/12] Well controls --- .../fluidFlow/wells/WellControls.cpp | 218 ++++++++--------- .../fluidFlow/wells/WellControls.hpp | 221 +++--------------- 2 files changed, 149 insertions(+), 290 deletions(-) diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/WellControls.cpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/WellControls.cpp index 45aa94574bc..59e6d9b7db1 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/WellControls.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/WellControls.cpp @@ -213,6 +213,91 @@ void WellControls::registerWellDataOnMesh( WellElementSubRegion & subRegion ) } } +template< typename > +struct isInjectionRateConstraint : std::false_type {}; +template< typename T > +struct isInjectionRateConstraint< InjectionConstraint< T > > : std::true_type {}; +template< bool IS_PRODUCER, typename T > +using RateConstraint = std::conditional_t< IS_PRODUCER, ProductionConstraint< T >, InjectionConstraint< T > >; + +namespace +{ + +// Unpacks a camp::list of types and forwards the lambda to the group's forSubGroups method. +template< typename GROUP, typename LAMBDA, typename ... Ts > +void forConstraints( GROUP & group, camp::list< Ts... >, LAMBDA && lambda ) +{ + group.template forSubGroups< Ts... >( std::forward< LAMBDA >( lambda )); +} + +// Constructs a compile-time list of constraint types based on boolean template parameters, +// then invokes the unpacked forConstraints function above. +template< bool IS_PRODUCER, bool INCLUDE_BHP, typename GROUP, typename LAMBDA > +void forConstraints( GROUP & group, LAMBDA && lambda ) +{ + using BHPConstraint = std::conditional_t< IS_PRODUCER, MinimumBHPConstraint, MaximumBHPConstraint >; + using RateTypes = camp::list< RateConstraint< IS_PRODUCER, PhaseVolumeRateConstraint >, + RateConstraint< IS_PRODUCER, VolumeRateConstraint >, + RateConstraint< IS_PRODUCER, MassRateConstraint > >; + using AllTypes = std::conditional_t< INCLUDE_BHP, + typename camp::prepend< RateTypes, BHPConstraint >::type, + RateTypes >; + forConstraints( group, AllTypes{}, std::forward< LAMBDA >( lambda ) ); +} +} + +template< typename LAMBDA > +void WellControls::forRateConstraints( LAMBDA && lambda ) const +{ + if( isProducer()) + { + forConstraints< true, false >( *this, std::forward< LAMBDA >( lambda ) ); + } + else + { + forConstraints< false, false >( *this, std::forward< LAMBDA >( lambda ) ); + } +} + +template< typename LAMBDA > +void WellControls::forRateConstraints( LAMBDA && lambda ) +{ + if( isProducer()) + { + forConstraints< true, false >( *this, std::forward< LAMBDA >( lambda ) ); + } + else + { + forConstraints< false, false >( *this, std::forward< LAMBDA >( lambda ) ); + } +} + +template< typename LAMBDA > +void WellControls::forAllConstraints( LAMBDA && lambda ) const +{ + if( isProducer()) + { + forConstraints< true, true >( *this, std::forward< LAMBDA >( lambda ) ); + } + else + { + forConstraints< false, true >( *this, std::forward< LAMBDA >( lambda ) ); + } +} + +template< typename LAMBDA > +void WellControls::forAllConstraints( LAMBDA && lambda ) +{ + if( isProducer()) + { + forConstraints< true, true >( *this, std::forward< LAMBDA >( lambda ) ); + } + else + { + forConstraints< false, true >( *this, std::forward< LAMBDA >( lambda ) ); + } +} + void WellControls::postInputInitialization() { Group::postInputInitialization(); @@ -286,7 +371,7 @@ void WellControls::postInputInitialization() void WellControls::postRestartInitialization( ) {} -void WellControls::logConstraint( WellConstraintBase const * constraint, WellElementSubRegion const & region, real64 time, bool isLimiting ) +void WellControls::logConstraint( WellConstraintBase const * constraint, WellElementSubRegion const & region, real64 time, bool isLimiting ) const { bool const needsLog = (constraint != nullptr) && (getLogLevel() > 4) && region.isLocallyOwned(); if( isLimiting ) @@ -305,24 +390,6 @@ void WellControls::logConstraint( WellConstraintBase const * constraint, WellEle } } -//template< typename T > -//using isConstraint = std::disjunction< -// std::is_same< T, MinimumBHPConstraint >, -//std::is_same< T, MaximumBHPConstraint >, -//std::is_same< T, ProductionConstraint< PhaseVolumeRateConstraint > >, -//std::is_same, -//std::is_same< T, ProductionConstraint< MassRateConstraint > >, -//std::is_same< T, InjectionConstraint< PhaseVolumeRateConstraint > >, -//std::is_same< T, InjectionConstraint< VolumeRateConstraint >, -// std::is_same< T, InjectionConstraint< MassRateConstraint > > -// > -//template< typename ConstraintType > -//ConstraintType const * WellControls::getFirstConstraintOfType() const -//{ -// static_assert( isConstraint< ConstraintType >::value, "ConstraintType must be a constraint type" ); -// return nullptr; -//} - void WellControls::setWellStatus( real64 const & currentTime, WellControls::Status status ) { m_wellStatus = status; @@ -335,23 +402,9 @@ void WellControls::setWellStatus( real64 const & currentTime, WellControls::Stat hasZeroRate = true; } - if( !hasZeroRate && isProducer()) - { - forSubGroups< ProductionConstraint< PhaseVolumeRateConstraint >, - ProductionConstraint< VolumeRateConstraint >, - ProductionConstraint< MassRateConstraint > >( [&] ( auto const & constraint ) - { - if( isZero( constraint.getConstraintValue( currentTime ) ) ) - { - hasZeroRate = true; - } - } ); - } - else if( !hasZeroRate ) + if( !hasZeroRate ) { - forSubGroups< InjectionConstraint< PhaseVolumeRateConstraint >, - InjectionConstraint< VolumeRateConstraint >, - InjectionConstraint< MassRateConstraint > >( [&] ( auto const & constraint ) + forRateConstraints( [&] ( auto const & constraint ) { if( isZero( constraint.getConstraintValue( currentTime ) ) ) { @@ -412,26 +465,10 @@ bool WellControls::getWellState() const void WellControls::setNextDtFromTables( real64 const & currentTime, real64 & nextDt ) { - if( isProducer() ) + forAllConstraints( [&] ( auto const & constraint ) { - forSubGroups< MinimumBHPConstraint, - ProductionConstraint< PhaseVolumeRateConstraint >, - ProductionConstraint< VolumeRateConstraint >, - ProductionConstraint< MassRateConstraint > >( [&] ( auto const & constraint ) - { - constraint.setNextDtFromTables( currentTime, nextDt ); - } ); - } - else - { - forSubGroups< MaximumBHPConstraint, - InjectionConstraint< PhaseVolumeRateConstraint >, - InjectionConstraint< VolumeRateConstraint >, - InjectionConstraint< MassRateConstraint > >( [&] ( auto const & constraint ) - { - constraint.setNextDtFromTables( currentTime, nextDt ); - } ); - } + constraint.setNextDtFromTables( currentTime, nextDt ); + } ); setNextDtFromTable( m_statusTable, currentTime, nextDt ); } @@ -473,14 +510,15 @@ real64 WellControls::getTargetBHP( real64 const & targetTime ) const real64 WellControls::getInjectionTemperature() const { real64 injectionTemperature = 0.0; - forSubGroups< InjectionConstraint< PhaseVolumeRateConstraint >, - InjectionConstraint< VolumeRateConstraint >, - InjectionConstraint< MassRateConstraint > >( [&] ( auto & constraint ) + forRateConstraints( [&] ( auto const & constraint ) { - if( constraint.isConstraintActive()) + if constexpr ( isInjectionRateConstraint< decltype(constraint) >::value ) { - injectionTemperature = constraint.getInjectionTemperature(); - return; + if( constraint.isConstraintActive()) + { + injectionTemperature = constraint.getInjectionTemperature(); + return; + } } } ); return injectionTemperature; @@ -489,17 +527,17 @@ real64 WellControls::getInjectionTemperature() const arrayView1d< real64 const > WellControls::getInjectionStream() const { arrayView1d< real64 const > injectionStream; - forSubGroups< InjectionConstraint< PhaseVolumeRateConstraint >, - InjectionConstraint< VolumeRateConstraint >, - InjectionConstraint< MassRateConstraint > >( [&] ( auto & constraint ) + forRateConstraints( [&] ( auto const & constraint ) { - if( constraint.isConstraintActive() ) + if constexpr ( isInjectionRateConstraint< decltype(constraint) >::value ) { - injectionStream = constraint.getInjectionStream(); - return; + if( constraint.isConstraintActive()) + { + injectionStream = constraint.getInjectionStream(); + return; + } } } ); - return injectionStream; } @@ -689,8 +727,6 @@ void WellControls::selectWellConstraint( real64 const & time_n, bool useEstimator = coupledIterationNumber < estimateSolution(); - - if( getWellState()) { if( useEstimator ) @@ -734,26 +770,11 @@ bool WellControls::evaluateConstraints( real64 const & time_n, // create list of all constraints to process stdVector< WellConstraintBase * > constraintList; - if( isProducer() ) - { - forSubGroups< MinimumBHPConstraint, - ProductionConstraint< PhaseVolumeRateConstraint >, - ProductionConstraint< VolumeRateConstraint >, - ProductionConstraint< MassRateConstraint > >( [&] ( auto & constraint ) - { - constraintList.emplace_back( &constraint ); - } ); - } - else + forAllConstraints( [&] ( auto & constraint ) { - forSubGroups< MaximumBHPConstraint, - InjectionConstraint< PhaseVolumeRateConstraint >, - InjectionConstraint< VolumeRateConstraint >, - InjectionConstraint< MassRateConstraint > >( [&] ( auto & constraint ) - { - constraintList.emplace_back( &constraint ); - } ); - } + constraintList.emplace_back( &constraint ); + } ); + // Get current constraint WellConstraintBase * limitingConstraint = nullptr; for( auto & constraint : constraintList ) @@ -801,15 +822,13 @@ bool WellControls::evaluateConstraints( real64 const & time_n, // note that initializeWells sets the initial constraint // tjb reactive control schema to allow use to set if needed stdVector< WellConstraintBase * > constraintList; + forRateConstraints( [&] ( auto & constraint ) + { + constraintList.emplace_back( &constraint ); + } ); WellConstraintBase * limitingConstraint = getCurrentConstraint(); if( isProducer() ) { - forSubGroups< ProductionConstraint< PhaseVolumeRateConstraint >, - ProductionConstraint< VolumeRateConstraint >, - ProductionConstraint< MassRateConstraint > >( [&] ( auto & constraint ) - { - constraintList.emplace_back( &constraint ); - } ); if( limitingConstraint->getControl() != ConstraintTypeId::BHP ) { { // remove from list and add BHP constraint @@ -833,12 +852,6 @@ bool WellControls::evaluateConstraints( real64 const & time_n, } else { - forSubGroups< InjectionConstraint< PhaseVolumeRateConstraint >, - InjectionConstraint< VolumeRateConstraint >, - InjectionConstraint< MassRateConstraint > >( [&] ( auto & constraint ) - { - constraintList.emplace_back( &constraint ); - } ); if( limitingConstraint->getControl() != ConstraintTypeId::BHP ) { forSubGroups< MaximumBHPConstraint >( [&] ( auto & constraint ) @@ -905,9 +918,7 @@ bool WellControls::evaluateConstraints( real64 const & time_n, subRegion, dofManager ); - GEOS_LOG_RANK_IF ( getLogLevel() > 4 && subRegion.isLocallyOwned(), - " Well " << subRegion.getName() << " Limiting Constraint " << limitingConstraint->getName() << " " << limitingConstraint->bottomHolePressure() << " " << limitingConstraint->phaseVolumeRates() << " " << - limitingConstraint->totalVolumeRate() << " " << limitingConstraint->massRate()); + logConstraint( limitingConstraint, subRegion, time_n, true ); return true; } @@ -945,7 +956,6 @@ void WellControls::solveConstraint( WellConstraintBase *constraint, bool useEstimator = coupledIterationNumber < estimateSolution(); if( useEstimator ) { - if( getLogLevel() > 4 ) { GEOS_LOG_RANK_0( "Well " <getName() << " value " << constraint->getConstraintValue( time_n ) << " active " << diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/WellControls.hpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/WellControls.hpp index 3a8343917c6..a4c21c38af7 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/WellControls.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/WellControls.hpp @@ -26,13 +26,6 @@ #include "functions/TableFunction.hpp" #include "constitutive/fluid/multifluid/MultiFluidBase.hpp" -//#include "physicsSolvers/fluidFlow/wells/WellInjectionConstraint.hpp" -//#include "physicsSolvers/fluidFlow/wells/WellProductionConstraint.hpp" -//#include "physicsSolvers/fluidFlow/wells/WellBHPConstraints.hpp" -//#include "physicsSolvers/fluidFlow/wells/WellVolumeRateConstraint.hpp" -//#include "physicsSolvers/fluidFlow/wells/WellPhaseVolumeRateConstraint.hpp" -//#include "physicsSolvers/fluidFlow/wells/WellMassRateConstraint.hpp" -//#include "physicsSolvers/fluidFlow/wells/WellLiquidRateConstraint.hpp" #include "constitutive/fluid/multifluid/MultiFluidBase.hpp" #include "constitutive/fluid/singlefluid/SingleFluidBase.hpp" #include "physicsSolvers/fluidFlow/wells/WellNewtonSolver.hpp" @@ -161,11 +154,11 @@ class WellControls : public dataRepository::Group virtual bool isCompositional() const = 0; /** - * * @brief Initialize well for the beginning of a simulation or restart - * @param domain the domain - * @param mesh the mesh level - * @param subRegion the well subRegion - * @param time_n the current time + * @brief Initialize well for the beginning of a simulation or restart + * @param domain the domain + * @param mesh the mesh level + * @param subRegion the well subRegion + * @param time_n the current time */ virtual void initializeWell( DomainPartition & domain, MeshLevel & mesh, WellElementSubRegion & subRegion, real64 const & time_n ) = 0; /** @@ -376,115 +369,7 @@ class WellControls : public dataRepository::Group virtual void printRates( real64 const & time_n, real64 const & dt, WellElementSubRegion const & subRegion ) = 0; -#if 0 - /**@}*/ - /** - * @brief Apply a given functor to a container if the container can be - * cast to one of the specified types. - * @tparam CASTTYPE the first type that will be used in the attempted casting of container - * @tparam CASTTYPES a variadic list of types that will be used in the attempted casting of container - * @tparam CONTAINERTYPE the type of container - * @tparam LAMBDA the type of lambda function to call in the function - * @param[in] container a pointer to the container which will be passed to the lambda function - * @param[in] lambda the lambda function to call in the function - * @return a boolean to indicate whether the lambda was successfully applied to the container. - */ - template< typename T0, typename T1, typename ... CASTTYPES, typename CONTAINERTYPE, typename LAMBDA > - static bool applyLambdaToContainer( CONTAINERTYPE container, LAMBDA && lambda ) - { - using Pointee = std::remove_pointer_t< std::remove_reference_t< CONTAINERTYPE > >; - using T = std::conditional_t< std::is_const< Pointee >::value, T0 const, T0 >; - T * const castedContainer = dynamic_cast< T * >( container ); - - if( castedContainer != nullptr ) - { - lambda( *castedContainer ); - return true; - } - - return applyLambdaToContainer< T1, CASTTYPES... >( container, std::forward< LAMBDA >( lambda ) ); - } - - // Base case: no more types to try - template< typename CONTAINERTYPE, typename LAMBDA > - static bool applyLambdaToContainer( CONTAINERTYPE /*container*/, LAMBDA && /*lambda*/ ) - { - return false; - } - - // Single-type overload: try only T0 and stop - template< typename T0, typename CONTAINERTYPE, typename LAMBDA > - static bool applyLambdaToContainer( CONTAINERTYPE container, LAMBDA && lambda ) - { - using Pointee = std::remove_pointer_t< std::remove_reference_t< CONTAINERTYPE > >; - using T = std::conditional_t< std::is_const< Pointee >::value, T0 const, T0 >; - T * const castedContainer = dynamic_cast< T * >( container ); - if( castedContainer != nullptr ) - { - lambda( *castedContainer ); - return true; - } - - return false; - } - - - /** - * @copydoc forInjectionConstraints(LAMBDA &&) - */ - template< typename GROUPTYPE = Group, typename ... GROUPTYPES, typename LAMBDA > - void forInjectionConstraints( LAMBDA && lambda ) const - { - for( auto const * constraintIter : m_injectionRateConstraintList ) - { - applyLambdaToContainer< GROUPTYPE, GROUPTYPES... >( constraintIter, [&]( auto const & castedSubGroup ) - { - lambda( castedSubGroup ); - } ); - } - } - - // non-const overload - template< typename GROUPTYPE = Group, typename ... GROUPTYPES, typename LAMBDA > - void forInjectionConstraints( LAMBDA && lambda ) - { - for( auto * constraintIter : m_injectionRateConstraintList ) - { - applyLambdaToContainer< GROUPTYPE, GROUPTYPES... >( constraintIter, [&]( auto & castedSubGroup ) - { - lambda( castedSubGroup ); - } ); - } - } - /** - * @copydoc forProductionConstraints(LAMBDA &&) - */ - template< typename GROUPTYPE = Group, typename ... GROUPTYPES, typename LAMBDA > - void forProductionConstraints( LAMBDA && lambda ) const - { - for( auto const * constraintIter : m_productionRateConstraintList ) - { - applyLambdaToContainer< GROUPTYPE, GROUPTYPES... >( constraintIter, [&]( auto const & castedSubGroup ) - { - lambda( castedSubGroup ); - } ); - } - } - - // non-const overload - template< typename GROUPTYPE = Group, typename ... GROUPTYPES, typename LAMBDA > - void forProductionConstraints( LAMBDA && lambda ) - { - for( auto * constraintIter : m_productionRateConstraintList ) - { - applyLambdaToContainer< GROUPTYPE, GROUPTYPES... >( constraintIter, [&]( auto & castedSubGroup ) - { - lambda( castedSubGroup ); - } ); - } - } -#endif /** * @name Getters / Setters */ @@ -570,7 +455,6 @@ class WellControls : public dataRepository::Group */ real64 getTargetBHP( real64 const & targetTime ) const; - /** * @brief Const accessor for the temperature of the injection stream * @return the temperature of the injection stream @@ -595,7 +479,6 @@ class WellControls : public dataRepository::Group */ real64 getReferenceElevation() const; - /** * @brief Getter for the flag specifying whether we check rates at surface or reservoir conditions * @return 1 if we use surface conditions, and 0 otherwise @@ -669,7 +552,6 @@ class WellControls : public dataRepository::Group */ bool getWellState() const; - /** * @brief Set the current consrtaint * @param[in] currentConstraint pointer to constraint @@ -679,18 +561,8 @@ class WellControls : public dataRepository::Group * @brief Get the current consrtaint * @return pointer to constraint */ - WellConstraintBase * getCurrentConstraint() { return m_currentConstraint; } - WellConstraintBase const * getCurrentConstraint() const { return m_currentConstraint; } - - /** - * @brief Get the BHP constraint. - * @details Returns the first BHP constraint attached to the well controls. This will be the minimum BHP - * constraint for a producer well and the maximum BHP constraint for an injector well. If there is no BHP - * constraint assigned to the well controls this will return a null value. - * @return The constraint object of one exists nullptr otherwise - */ - //WellConstraintBase const * getBHPConstraint() const; - //WellConstraintBase * getBHPConstraint(); + WellConstraintBase * getCurrentConstraint() { return m_currentConstraint; } + WellConstraintBase const * getCurrentConstraint() const { return m_currentConstraint; } /** * @brief Get the rate constraint. @@ -701,7 +573,6 @@ class WellControls : public dataRepository::Group * @return The constraint object of one exists nullptr otherwise */ WellConstraintBase const * getRateConstraint() const; - //WellConstraintBase * getRateConstraint(); /** * @brief Getter for the flag to enable crossflow @@ -732,7 +603,6 @@ class WellControls : public dataRepository::Group void setKeepVariablesConstantDuringInitStep( bool const keepVariablesConstantDuringInitStep ) { m_keepVariablesConstantDuringInitStep = keepVariablesConstantDuringInitStep; } - /** * @brief setter for multi fluid separator * @param[in] fluidSeparatorPtr single or multiphase separator @@ -872,26 +742,6 @@ class WellControls : public dataRepository::Group */ struct groupKeyStruct { -#if 0 - /// string key for the minimum BHP presssure for a producer - static constexpr char const * minimumBHPConstraintString() { return "MinimumBHPConstraint"; } - /// string key for the maximum BHP presssure for a injection - static constexpr char const * maximumBHPConstraintString() { return "MaximumBHPConstraint"; } - /// string key for the maximum phase rate for a producer - static constexpr char const * productionPhaseVolumeRateConstraintString() { return "ProductionPhaseVolumeRateConstraint"; } - /// string key for the maximum phase rate for a injection - static constexpr char const * injectionPhaseVolumeRateConstraint() { return "InjectionPhaseVolumeRateConstraint"; } - /// string key for the maximum volume rate for a producer - static constexpr char const * productionVolumeRateConstraint() { return "ProductionVolumeRateConstraint"; } - /// string key for the maximum volume rate for a injector - static constexpr char const * injectionVolumeRateConstraint() { return "InjectionVolumeRateConstraint"; } - /// string key for the maximum mass rate for a producer - static constexpr char const * productionMassRateConstraint() { return "ProductionMassRateConstraint"; } - /// string key for the maximum mass rate for a injector - static constexpr char const * injectionMassRateConstraint() { return "InjectionMassRateConstraint"; } - /// string key for the liquid rate for a producer - static constexpr char const * productionLiquidRateConstraint() { return "ProductionLiquidRateConstraint"; } -#endif /// string key for the well Newton solver static constexpr char const * wellNewtonSolverString() { return "WellNewtonSolver"; } }; @@ -914,24 +764,6 @@ class WellControls : public dataRepository::Group */ template< typename ConstraintType > void createConstraint ( string const & constraintName ); -#if 0 - /** - * @brief Getters for constraints - */ - BHPConstraint< BHPConstraintTypeId::MIN > * getMinBHPConstraint() { return nullptr; }; - BHPConstraint< BHPConstraintTypeId::MIN > * getMinBHPConstraint() const { return nullptr; }; - BHPConstraint< BHPConstraintTypeId::MAX > * getMaxBHPConstraint() { return nullptr; }; - BHPConstraint< BHPConstraintTypeId::MAX > * getMaxBHPConstraint() const { return nullptr; }; - - /** - * @brief Getters for constraint lists - */ - std::vector< WellConstraintBase * > getProdRateConstraints() { return {}; }; - std::vector< WellConstraintBase * > getProdRateConstraints() const { return {}; }; - std::vector< WellConstraintBase * > getInjRateConstraints() { return {}; } - std::vector< WellConstraintBase * > getInjRateConstraints() const { return {}; } -#endif - /** * @brief Set thermal effects enable * @param[in] true/false @@ -966,7 +798,33 @@ class WellControls : public dataRepository::Group protected: virtual void postRestartInitialization( )override; - void logConstraint( WellConstraintBase const * constraint, WellElementSubRegion const & region, real64 time, bool isLimiting = false ); + /** + * @brief Iterates over rate constraints and applies a lambda function. + * @details This constant member function checks whether the well is a producer or an injector. + * It then dispatches the provided lambda to the underlying rate-specific constraint subgroups + * (PhaseVolumeRateConstraint, VolumeRateConstraint, MassRateConstraint) while explicitly + * excluding BHP constraints. + * @param lambda The callable object to apply to each rate constraint subgroup. + */ + template< typename LAMBDA > + void forRateConstraints( LAMBDA && lambda ) const; + template< typename LAMBDA > + void forRateConstraints( LAMBDA && lambda ); + + /** + * @brief Iterates over all constraints, including BHP, and applies a lambda function. + * @details This constant member function evaluates if the well is a producer or an injector. + * It then applies the provided lambda to all constraint subgroups, which includes both the rate + * constraints and the appropriate Bottom Hole Pressure constraint (MinimumBHPConstraint for + * producers and MaximumBHPConstraint for injectors). + * @param lambda The callable object to apply to each constraint subgroup. + */ + template< typename LAMBDA > + void forAllConstraints( LAMBDA && lambda ) const; + template< typename LAMBDA > + void forAllConstraints( LAMBDA && lambda ); + + void logConstraint( WellConstraintBase const * constraint, WellElementSubRegion const & region, real64 time, bool isLimiting = false ) const; private: /// List of names of regions the solver will be applied to @@ -1046,16 +904,7 @@ class WellControls : public dataRepository::Group real64 m_initialPressureCoefficient; // Current constrint - WellConstraintBase * m_currentConstraint; -#if 0 - // Minimum and maximum BHP and WHP constraints - BHPConstraint< BHPConstraintTypeId::MIN > * m_minBHPConstraint; - BHPConstraint< BHPConstraintTypeId::MAX > * m_maxBHPConstraint; - - // Lists of rate constraints - std::vector< WellConstraintBase * > m_productionRateConstraintList; - std::vector< WellConstraintBase * > m_injectionRateConstraintList; -#endif + WellConstraintBase * m_currentConstraint{}; /// Well status WellControls::Status m_wellStatus; From 410ed120de4509080be9590f6673cf26406f16ff Mon Sep 17 00:00:00 2001 From: dkachuma Date: Thu, 21 May 2026 14:27:15 -0500 Subject: [PATCH 03/12] implement rate constraint getter --- .../fluidFlow/wells/WellControls.cpp | 15 ++++++++++++++- 1 file changed, 14 insertions(+), 1 deletion(-) diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/WellControls.cpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/WellControls.cpp index 59e6d9b7db1..39964e0b5eb 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/WellControls.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/WellControls.cpp @@ -366,11 +366,24 @@ void WellControls::postInputInitialization() m_statusTable->getName() ), InputError, getDataContext() ); } - } + void WellControls::postRestartInitialization( ) {} +WellConstraintBase const * WellControls::getRateConstraint() const +{ + WellConstraintBase const * rateConstraint = nullptr; + forRateConstraints( [&]( auto const & constraint ) + { + if( rateConstraint == nullptr ) + { + rateConstraint = &constraint; + } + } ); + return rateConstraint; +} + void WellControls::logConstraint( WellConstraintBase const * constraint, WellElementSubRegion const & region, real64 time, bool isLimiting ) const { bool const needsLog = (constraint != nullptr) && (getLogLevel() > 4) && region.isLocallyOwned(); From 3ace55bec4d06aa009f9d8a2631a7cef3638de91 Mon Sep 17 00:00:00 2001 From: dkachuma Date: Thu, 21 May 2026 14:34:54 -0500 Subject: [PATCH 04/12] Undo some changes --- .../fluidFlow/wells/WellBHPConstraints.cpp | 5 +- .../fluidFlow/wells/WellBHPConstraints.hpp | 4 +- .../fluidFlow/wells/WellManager.hpp | 106 ++++++++++++++++++ 3 files changed, 111 insertions(+), 4 deletions(-) diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/WellBHPConstraints.cpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/WellBHPConstraints.cpp index ac03db749b4..c6fbca6e40c 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/WellBHPConstraints.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/WellBHPConstraints.cpp @@ -84,13 +84,14 @@ bool BHPConstraint< T >::checkViolation( WellConstraintBase const & currentConst } } + template class BHPConstraint< BHPConstraintTypeId::MIN >; template class BHPConstraint< BHPConstraintTypeId::MAX >; - namespace { +typedef BHPConstraint< BHPConstraintTypeId::MIN > MinimumBHPConstraint; +typedef BHPConstraint< BHPConstraintTypeId::MAX > MaximumBHPConstraint; REGISTER_CATALOG_ENTRY( WellConstraintBase, MinimumBHPConstraint, string const &, Group * const ) REGISTER_CATALOG_ENTRY( WellConstraintBase, MaximumBHPConstraint, string const &, Group * const ) } - } //namespace geos diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/WellBHPConstraints.hpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/WellBHPConstraints.hpp index 28f356435c3..4ce9312d019 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/WellBHPConstraints.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/WellBHPConstraints.hpp @@ -24,7 +24,6 @@ #include "dataRepository/Group.hpp" #include "functions/TableFunction.hpp" #include "WellConstraintsBase.hpp" - namespace geos { @@ -109,7 +108,8 @@ class BHPConstraint : public WellConstraintBase static constexpr char const * targetBHPString() { return "targetBHP"; } /// String key for the well reference elevation (for BHP control) static constexpr char const * refElevString() { return "referenceElevation"; } - }; + } + viewKeysWellBHPConstraint; virtual bool checkViolation( WellConstraintBase const & currentConstraint, real64 const & currentTime ) const override; diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/WellManager.hpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/WellManager.hpp index 9a3a2e53260..db764b24072 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/WellManager.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/WellManager.hpp @@ -407,6 +407,109 @@ class WellManager : public PhysicsSolverBase * @param domain the physical domain object */ void chopNegativeDensities( DomainPartition & domain ); +#if 0 + /** + * @brief calculate the norm of the global system residual + * @param time the time at the beginning of the step + * @param dt the desired timestep + * @param domain the domain partition + * @param dofManager degree-of-freedom manager associated with the linear system + * @param localRhs the system right-hand side vector + * @return norm of the residual + * + * This function returns the norm of global residual vector, which is suitable for comparison with + * a tolerance. + */ + virtual real64 + calculateResidualNorm( real64 const & time, + real64 const & dt, + DomainPartition const & domain, + DofManager const & dofManager, + arrayView1d< real64 const > const & localRhs ); + /** + * @brief Function to check system solution for physical consistency and constraint violation + * @param domain the domain partition + * @param dofManager degree-of-freedom manager associated with the linear system + * @param localSolution the solution vector + * @param scalingFactor factor to scale the solution prior to application + * @return true if solution can be safely applied without violating physical constraints, false otherwise + * + * @note This function must be overridden in the derived physics solver in order to use an implict + * solution method such as LinearImplicitStep() or NonlinearImplicitStep(). + * + */ + virtual bool + checkSystemSolution( DomainPartition & domain, + DofManager const & dofManager, + arrayView1d< real64 const > const & localSolution, + real64 const scalingFactor ); + + /** + * @brief Function to determine if the solution vector should be scaled back in order to maintain a known constraint. + * @param[in] domain The domain partition. + * @param[in] dofManager degree-of-freedom manager associated with the linear system + * @param[in] localSolution the solution vector + * @return The factor that should be used to scale the solution vector values when they are being applied. + */ + virtual real64 + scalingForSystemSolution( DomainPartition & domain, + DofManager const & dofManager, + arrayView1d< real64 const > const & localSolution ); + + /** + * @brief Function to apply the solution vector to the state + * @param dofManager degree-of-freedom manager associated with the linear system + * @param localSolution the solution vector + * @param scalingFactor factor to scale the solution prior to application + * @param dt the timestep + * @param domain the domain partition + * + * This function performs 2 operations: + * 1) extract the solution vector for the "blockSystem" parameter, and applies the + * contents of the solution vector to the primary variable field data, + * 2) perform a synchronization of the primary field variable such that all ghosts are updated, + * + * The "scalingFactor" parameter allows for the scaled application of the solution vector. For + * instance, a line search may apply a negative scaling factor to remove part of the previously + * applied solution. + * + * @note This function must be overridden in the derived physics solver in order to use an implict + * solution method such as LinearImplicitStep() or NonlinearImplicitStep(). + * + */ + virtual void + applySystemSolution( DofManager const & dofManager, + arrayView1d< real64 const > const & localSolution, + real64 const scalingFactor, + real64 const dt, + DomainPartition & domain ); + + + /** + * @brief Recompute all dependent quantities from primary variables (including constitutive models) + * @param domain the domain containing the mesh and fields + */ + virtual void updateState( DomainPartition & domain ); + + /** + * @brief perform cleanup for implicit timestep + * @param time the time at the beginning of the step + * @param dt the desired timestep + * @param domain the domain partition + * + * This function performs whatever tasks are required to complete an implicit timestep. For + * example, the acceptance of the solution will occur during this step, and deallocation of + * temporaries will be be performed in this function. + * + * @note This function must be overridden in the derived physics solver in order to use an implict + * solution method such as LinearImplicitStep() or NonlinearImplicitStep(). + */ + virtual void + implicitStepComplete( real64 const & time, + real64 const & dt, + DomainPartition & domain ) override; + +#endif /** * @brief Utility function to keep the well variables during a time step (used in @@ -420,12 +523,15 @@ class WellManager : public PhysicsSolverBase protected: + //virtual void postInputInitialization() override; + virtual void initializePostSubGroups() override; virtual void initializePostInitialConditionsPreSubGroups() override; virtual void postRestartInitialization() override final; + private: /// name of the flow solver From c85527f7946d5588e870bdc211d64eeced0db47d Mon Sep 17 00:00:00 2001 From: dkachuma Date: Thu, 21 May 2026 15:00:05 -0500 Subject: [PATCH 05/12] Add children --- .../fluidFlow/wells/WellControls.cpp | 18 ++++++++++++++++++ .../fluidFlow/wells/WellControls.hpp | 11 +++++++++++ 2 files changed, 29 insertions(+) diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/WellControls.cpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/WellControls.cpp index 39964e0b5eb..e4cb96a4e1c 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/WellControls.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/WellControls.cpp @@ -160,6 +160,24 @@ WellControls::WellControls( string const & name, Group * const parent ) WellControls::~WellControls() {} +Group * WellControls::createChild( string const & childKey, string const & childName ) +{ + GEOS_LOG_RANK_0( GEOS_FMT( "{}: adding {} {}", getName(), childKey, childName ) ); + std::unique_ptr< WellConstraintBase > constraint = + WellConstraintBase::CatalogInterface::factory( childKey, getDataContext(), childName, this ); + return ®isterGroup< WellConstraintBase >( childName, std::move( constraint ) ); +} + +void WellControls::expandObjectCatalogs() +{ + // During schema generation, register one of each type derived from WellConstraintBase here + for( auto & catalogIter : WellConstraintBase::getCatalog()) + { + createChild( catalogIter.first, catalogIter.first ); + } + // tjbcreateChild( WellNewtonSolver::catalogName(), WellNewtonSolver::catalogName() ); +} + namespace { diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/WellControls.hpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/WellControls.hpp index a4c21c38af7..4d54f6a5bfc 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/WellControls.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/WellControls.hpp @@ -133,6 +133,17 @@ class WellControls : public dataRepository::Group /// String used to form the solverName used to register single-physics solvers in CoupledSolver static string coupledSolverAttributePrefix() { return "well"; } + /** + * @brief Create a new constraint object as a child of this group. + * @param childKey the catalog key of the new constraint object to create + * @param childName the name of the new constraint object in the repository + * @return the group child + */ + virtual Group * createChild( string const & childKey, string const & childName ) override; + + /// Expand catalog for schema generation + virtual void expandObjectCatalogs() override; + virtual void registerWellDataOnMesh( WellElementSubRegion & subRegion ); virtual void setConstitutiveNames( ElementSubRegionBase & subRegion ) const = 0; /** From 3f7a47a7cd37829aee4e8fc435fad5a579b6594c Mon Sep 17 00:00:00 2001 From: dkachuma Date: Thu, 21 May 2026 15:53:31 -0500 Subject: [PATCH 06/12] Fix bypass --- .../physicsSolvers/fluidFlow/wells/WellControls.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/WellControls.cpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/WellControls.cpp index e4cb96a4e1c..df586745ed3 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/WellControls.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/WellControls.cpp @@ -543,7 +543,7 @@ real64 WellControls::getInjectionTemperature() const real64 injectionTemperature = 0.0; forRateConstraints( [&] ( auto const & constraint ) { - if constexpr ( isInjectionRateConstraint< decltype(constraint) >::value ) + if constexpr ( isInjectionRateConstraint< std::decay_t< decltype(constraint) > >::value ) { if( constraint.isConstraintActive()) { @@ -560,7 +560,7 @@ arrayView1d< real64 const > WellControls::getInjectionStream() const arrayView1d< real64 const > injectionStream; forRateConstraints( [&] ( auto const & constraint ) { - if constexpr ( isInjectionRateConstraint< decltype(constraint) >::value ) + if constexpr ( isInjectionRateConstraint< std::decay_t< decltype(constraint) > >::value ) { if( constraint.isConstraintActive()) { From d7f3582cd2268ac3e7f1c8804b7062ca3575903d Mon Sep 17 00:00:00 2001 From: dkachuma Date: Thu, 21 May 2026 21:15:05 -0500 Subject: [PATCH 07/12] Add vector getters --- .../fluidFlow/wells/WellControls.cpp | 230 +++++------------- .../fluidFlow/wells/WellControls.hpp | 45 +--- .../CompositionalMultiphaseWellKernels.cpp | 4 +- .../CompositionalMultiphaseWellKernels.hpp | 2 +- .../wells/kernels/SinglePhaseWellKernels.cpp | 2 +- ...rmalCompositionalMultiphaseWellKernels.hpp | 2 +- 6 files changed, 70 insertions(+), 215 deletions(-) diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/WellControls.cpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/WellControls.cpp index df586745ed3..bb9dda95d69 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/WellControls.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/WellControls.cpp @@ -36,6 +36,7 @@ #include "physicsSolvers/fluidFlow/wells/WellFields.hpp" #include "functions/FunctionManager.hpp" + namespace geos { @@ -231,91 +232,6 @@ void WellControls::registerWellDataOnMesh( WellElementSubRegion & subRegion ) } } -template< typename > -struct isInjectionRateConstraint : std::false_type {}; -template< typename T > -struct isInjectionRateConstraint< InjectionConstraint< T > > : std::true_type {}; -template< bool IS_PRODUCER, typename T > -using RateConstraint = std::conditional_t< IS_PRODUCER, ProductionConstraint< T >, InjectionConstraint< T > >; - -namespace -{ - -// Unpacks a camp::list of types and forwards the lambda to the group's forSubGroups method. -template< typename GROUP, typename LAMBDA, typename ... Ts > -void forConstraints( GROUP & group, camp::list< Ts... >, LAMBDA && lambda ) -{ - group.template forSubGroups< Ts... >( std::forward< LAMBDA >( lambda )); -} - -// Constructs a compile-time list of constraint types based on boolean template parameters, -// then invokes the unpacked forConstraints function above. -template< bool IS_PRODUCER, bool INCLUDE_BHP, typename GROUP, typename LAMBDA > -void forConstraints( GROUP & group, LAMBDA && lambda ) -{ - using BHPConstraint = std::conditional_t< IS_PRODUCER, MinimumBHPConstraint, MaximumBHPConstraint >; - using RateTypes = camp::list< RateConstraint< IS_PRODUCER, PhaseVolumeRateConstraint >, - RateConstraint< IS_PRODUCER, VolumeRateConstraint >, - RateConstraint< IS_PRODUCER, MassRateConstraint > >; - using AllTypes = std::conditional_t< INCLUDE_BHP, - typename camp::prepend< RateTypes, BHPConstraint >::type, - RateTypes >; - forConstraints( group, AllTypes{}, std::forward< LAMBDA >( lambda ) ); -} -} - -template< typename LAMBDA > -void WellControls::forRateConstraints( LAMBDA && lambda ) const -{ - if( isProducer()) - { - forConstraints< true, false >( *this, std::forward< LAMBDA >( lambda ) ); - } - else - { - forConstraints< false, false >( *this, std::forward< LAMBDA >( lambda ) ); - } -} - -template< typename LAMBDA > -void WellControls::forRateConstraints( LAMBDA && lambda ) -{ - if( isProducer()) - { - forConstraints< true, false >( *this, std::forward< LAMBDA >( lambda ) ); - } - else - { - forConstraints< false, false >( *this, std::forward< LAMBDA >( lambda ) ); - } -} - -template< typename LAMBDA > -void WellControls::forAllConstraints( LAMBDA && lambda ) const -{ - if( isProducer()) - { - forConstraints< true, true >( *this, std::forward< LAMBDA >( lambda ) ); - } - else - { - forConstraints< false, true >( *this, std::forward< LAMBDA >( lambda ) ); - } -} - -template< typename LAMBDA > -void WellControls::forAllConstraints( LAMBDA && lambda ) -{ - if( isProducer()) - { - forConstraints< true, true >( *this, std::forward< LAMBDA >( lambda ) ); - } - else - { - forConstraints< false, true >( *this, std::forward< LAMBDA >( lambda ) ); - } -} - void WellControls::postInputInitialization() { Group::postInputInitialization(); @@ -384,24 +300,27 @@ void WellControls::postInputInitialization() m_statusTable->getName() ), InputError, getDataContext() ); } + +//## // 13) Validate constraints +//## bool const producer = isProducer(); +//## ((void)producer); +//## mapBase typeCount; +//## for( auto & catalog : WellConstraintBase::getCatalog()) +//## { +//## typeCount[catalog.first] = 0; +//## } +//## forSubGroups< +//## >([&](auto const & group) +//## { +//## string const catalogName = group.getCatalogName(); +//## ((void)catalogName); +//##//if (typeCount. in ) +//## }); } void WellControls::postRestartInitialization( ) {} -WellConstraintBase const * WellControls::getRateConstraint() const -{ - WellConstraintBase const * rateConstraint = nullptr; - forRateConstraints( [&]( auto const & constraint ) - { - if( rateConstraint == nullptr ) - { - rateConstraint = &constraint; - } - } ); - return rateConstraint; -} - void WellControls::logConstraint( WellConstraintBase const * constraint, WellElementSubRegion const & region, real64 time, bool isLimiting ) const { bool const needsLog = (constraint != nullptr) && (getLogLevel() > 4) && region.isLocallyOwned(); @@ -435,13 +354,13 @@ void WellControls::setWellStatus( real64 const & currentTime, WellControls::Stat if( !hasZeroRate ) { - forRateConstraints( [&] ( auto const & constraint ) + for( auto const * constraint : getRateConstraints() ) { - if( isZero( constraint.getConstraintValue( currentTime ) ) ) + if( isZero( constraint->getConstraintValue( currentTime ) ) ) { hasZeroRate = true; } - } ); + } } if( hasZeroRate ) @@ -496,10 +415,10 @@ bool WellControls::getWellState() const void WellControls::setNextDtFromTables( real64 const & currentTime, real64 & nextDt ) { - forAllConstraints( [&] ( auto const & constraint ) + for( auto const * constraint : getAllConstraints() ) { - constraint.setNextDtFromTables( currentTime, nextDt ); - } ); + constraint->setNextDtFromTables( currentTime, nextDt ); + } setNextDtFromTable( m_statusTable, currentTime, nextDt ); } @@ -520,35 +439,24 @@ void WellControls::setNextDtFromTable( TableFunction const * table, real64 const real64 WellControls::getTargetBHP( real64 const & targetTime ) const { - WellConstraintBase const * bhpConstraint = nullptr; - if( isProducer()) - { - forSubGroups< MinimumBHPConstraint >( [&] ( auto const & constraint ) - { - bhpConstraint = &constraint; - } ); - } - else - { - forSubGroups< MaximumBHPConstraint >( [&] ( auto const & constraint ) - { - bhpConstraint = &constraint; - } ); - } - return bhpConstraint->getConstraintValue( targetTime ); + return getBHPConstraint()->getConstraintValue( targetTime ); } real64 WellControls::getInjectionTemperature() const { real64 injectionTemperature = 0.0; - forRateConstraints( [&] ( auto const & constraint ) + localIndex firstIndex = -1; // Used to "capture" the first one + forSubGroupsIndex< InjectionConstraint< MassRateConstraint >, + InjectionConstraint< VolumeRateConstraint >, + InjectionConstraint< PhaseVolumeRateConstraint >, + InjectionConstraint< LiquidRateConstraint > >( [&] ( localIndex index, auto const & constraint ) { - if constexpr ( isInjectionRateConstraint< std::decay_t< decltype(constraint) > >::value ) + if( firstIndex < 0 ) { - if( constraint.isConstraintActive()) + if( constraint.isConstraintActive() ) { injectionTemperature = constraint.getInjectionTemperature(); - return; + firstIndex = index; } } } ); @@ -558,14 +466,18 @@ real64 WellControls::getInjectionTemperature() const arrayView1d< real64 const > WellControls::getInjectionStream() const { arrayView1d< real64 const > injectionStream; - forRateConstraints( [&] ( auto const & constraint ) + localIndex firstIndex = -1; // Used to "capture" the first one + forSubGroupsIndex< InjectionConstraint< MassRateConstraint >, + InjectionConstraint< VolumeRateConstraint >, + InjectionConstraint< PhaseVolumeRateConstraint >, + InjectionConstraint< LiquidRateConstraint > >( [&] ( localIndex index, auto const & constraint ) { - if constexpr ( isInjectionRateConstraint< std::decay_t< decltype(constraint) > >::value ) + if( firstIndex < 0 ) { if( constraint.isConstraintActive()) { injectionStream = constraint.getInjectionStream(); - return; + firstIndex = index; } } } ); @@ -575,48 +487,29 @@ arrayView1d< real64 const > WellControls::getInjectionStream() const integer WellControls::getConstraintPhaseIndex() const { integer phaseIndex = -1; - - if( isProducer() ) - { - forSubGroups< ProductionConstraint< PhaseVolumeRateConstraint > >( [&] ( auto & constraint ) - { - if( constraint.isConstraintActive() ) - { - phaseIndex = constraint.getPhaseIndex(); - } - } ); - } - else + // Validation should make sure we are not mixing constraints. + // Here we assume that we have zero or one or the other but not both + forSubGroups< ProductionConstraint< PhaseVolumeRateConstraint >, + InjectionConstraint< PhaseVolumeRateConstraint > >( [&] ( auto & constraint ) { - forSubGroups< InjectionConstraint< PhaseVolumeRateConstraint > >( [&] ( auto & constraint ) + if( constraint.isConstraintActive() ) { - if( constraint.isConstraintActive() ) - { - phaseIndex = constraint.getPhaseIndex(); - } - } ); - } - + phaseIndex = constraint.getPhaseIndex(); + } + } ); return phaseIndex; } real64 WellControls::getReferenceElevation() const { real64 referenceElevation = 0.0; - if( isProducer () ) + // Validation should make sure we are not mixing constraints. + // Here we assume that we have zero or one or the other but not both + forSubGroups< MinimumBHPConstraint, + MaximumBHPConstraint >( [&] ( auto & constraint ) { - forSubGroups< MinimumBHPConstraint >( [&] ( auto & constraint ) - { - referenceElevation = constraint.getReferenceElevation(); - } ); - } - else - { - forSubGroups< MaximumBHPConstraint >( [&] ( auto & constraint ) - { - referenceElevation = constraint.getReferenceElevation(); - } ); - } + referenceElevation = constraint.getReferenceElevation(); + } ); return referenceElevation; } @@ -800,11 +693,7 @@ bool WellControls::evaluateConstraints( real64 const & time_n, { // create list of all constraints to process - stdVector< WellConstraintBase * > constraintList; - forAllConstraints( [&] ( auto & constraint ) - { - constraintList.emplace_back( &constraint ); - } ); + stdVector< WellConstraintBase * > constraintList = getAllConstraints(); // Get current constraint WellConstraintBase * limitingConstraint = nullptr; @@ -852,11 +741,7 @@ bool WellControls::evaluateConstraints( real64 const & time_n, // create list of all constraints to solve // note that initializeWells sets the initial constraint // tjb reactive control schema to allow use to set if needed - stdVector< WellConstraintBase * > constraintList; - forRateConstraints( [&] ( auto & constraint ) - { - constraintList.emplace_back( &constraint ); - } ); + stdVector< WellConstraintBase * > constraintList = getRateConstraints(); WellConstraintBase * limitingConstraint = getCurrentConstraint(); if( isProducer() ) { @@ -885,10 +770,7 @@ bool WellControls::evaluateConstraints( real64 const & time_n, { if( limitingConstraint->getControl() != ConstraintTypeId::BHP ) { - forSubGroups< MaximumBHPConstraint >( [&] ( auto & constraint ) - { - constraintList.emplace_back( &constraint ); - } ); + constraintList.emplace_back( getBHPConstraint() ); } } if( isoThermalEstimatorEnabled() ) diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/WellControls.hpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/WellControls.hpp index 4d54f6a5bfc..f5e529862c1 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/WellControls.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/WellControls.hpp @@ -575,16 +575,6 @@ class WellControls : public dataRepository::Group WellConstraintBase * getCurrentConstraint() { return m_currentConstraint; } WellConstraintBase const * getCurrentConstraint() const { return m_currentConstraint; } - /** - * @brief Get the rate constraint. - * @details Returns the first rate constraint attached to the well controls. This will be a prodcution rate - * constraint for a producer well and an injection rate constraint for an injector well. This will return the - * first constraint assigned regardless of type (phase volume, liquid volume or mass). If there is no rate - * constraint assigned to the well controls this will return a null value. - * @return The constraint object of one exists nullptr otherwise - */ - WellConstraintBase const * getRateConstraint() const; - /** * @brief Getter for the flag to enable crossflow * @return the flag deciding whether crossflow is allowed or not @@ -775,6 +765,15 @@ class WellControls : public dataRepository::Group */ template< typename ConstraintType > void createConstraint ( string const & constraintName ); + WellConstraintBase const * getBHPConstraint() const; + WellConstraintBase * getBHPConstraint(); + + stdVector< WellConstraintBase const * > getRateConstraints() const; + stdVector< WellConstraintBase * > getRateConstraints(); + + stdVector< WellConstraintBase const * > getAllConstraints() const; + stdVector< WellConstraintBase * > getAllConstraints(); + /** * @brief Set thermal effects enable * @param[in] true/false @@ -809,32 +808,6 @@ class WellControls : public dataRepository::Group protected: virtual void postRestartInitialization( )override; - /** - * @brief Iterates over rate constraints and applies a lambda function. - * @details This constant member function checks whether the well is a producer or an injector. - * It then dispatches the provided lambda to the underlying rate-specific constraint subgroups - * (PhaseVolumeRateConstraint, VolumeRateConstraint, MassRateConstraint) while explicitly - * excluding BHP constraints. - * @param lambda The callable object to apply to each rate constraint subgroup. - */ - template< typename LAMBDA > - void forRateConstraints( LAMBDA && lambda ) const; - template< typename LAMBDA > - void forRateConstraints( LAMBDA && lambda ); - - /** - * @brief Iterates over all constraints, including BHP, and applies a lambda function. - * @details This constant member function evaluates if the well is a producer or an injector. - * It then applies the provided lambda to all constraint subgroups, which includes both the rate - * constraints and the appropriate Bottom Hole Pressure constraint (MinimumBHPConstraint for - * producers and MaximumBHPConstraint for injectors). - * @param lambda The callable object to apply to each constraint subgroup. - */ - template< typename LAMBDA > - void forAllConstraints( LAMBDA && lambda ) const; - template< typename LAMBDA > - void forAllConstraints( LAMBDA && lambda ); - void logConstraint( WellConstraintBase const * constraint, WellElementSubRegion const & region, real64 time, bool isLimiting = false ) const; private: diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/kernels/CompositionalMultiphaseWellKernels.cpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/kernels/CompositionalMultiphaseWellKernels.cpp index 69ad0227e06..673b4d4dd84 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/kernels/CompositionalMultiphaseWellKernels.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/kernels/CompositionalMultiphaseWellKernels.cpp @@ -466,7 +466,7 @@ RateInitializationKernel:: { // this assumes phase control present integer const targetPhaseIndex = wellControls.getConstraintPhaseIndex(); - auto const * rateConstraint = wellControls.getRateConstraint(); + auto const * rateConstraint = wellControls.getRateConstraints().front(); // Use first rate constraint to set initial connection rates real64 const rateVal = rateConstraint->getConstraintValue( time ); forAll< parallelDevicePolicy<> >( subRegionSize, [=] GEOS_HOST_DEVICE ( localIndex const iwelem ) @@ -506,7 +506,7 @@ RateInitializationKernel:: } else if( controlType == ConstraintTypeId::BHP ) { - auto const * rateConstraint = wellControls.getRateConstraint(); + auto const * rateConstraint = wellControls.getRateConstraints().front(); // Use first rate constraint to set initial connection rates real64 const rateVal = rateConstraint->getConstraintValue( time ); forAll< parallelDevicePolicy<> >( subRegionSize, [=] GEOS_HOST_DEVICE ( localIndex const iwelem ) diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/kernels/CompositionalMultiphaseWellKernels.hpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/kernels/CompositionalMultiphaseWellKernels.hpp index 2da2725cf39..2087eae9c9d 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/kernels/CompositionalMultiphaseWellKernels.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/kernels/CompositionalMultiphaseWellKernels.hpp @@ -528,7 +528,7 @@ class ResidualNormKernel : public physicsSolverBaseKernels::ResidualNormKernelBa // This is a normalizer for the balance equations. The normalizaer should be the current rate not the constraint value!! // This is one of the reasons for restricting constraint type for a production well // another pr will remove fix this (so the cause for difference results is isolated to one change) - auto const * rateConstraint = wellControls.getRateConstraint(); + auto const * rateConstraint = wellControls.getRateConstraints().front(); if( rateConstraint != nullptr ) { m_constraintValue = rateConstraint->getConstraintValue( time ); diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/kernels/SinglePhaseWellKernels.cpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/kernels/SinglePhaseWellKernels.cpp index 47b18cf5657..ab9e50c7646 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/kernels/SinglePhaseWellKernels.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/kernels/SinglePhaseWellKernels.cpp @@ -550,7 +550,7 @@ RateInitializationKernel:: WellControls::Control const control = wellControls.getControl(); bool const isProducer = wellControls.isProducer(); - auto const * rateConstraint = wellControls.getRateConstraint(); + auto const * rateConstraint = wellControls.getRateConstraints().front(); real64 const constraintVal = rateConstraint->getConstraintValue( currentTime ); // Estimate the connection rates diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/kernels/ThermalCompositionalMultiphaseWellKernels.hpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/kernels/ThermalCompositionalMultiphaseWellKernels.hpp index 1c6d39c290e..29e97edaaf2 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/kernels/ThermalCompositionalMultiphaseWellKernels.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/kernels/ThermalCompositionalMultiphaseWellKernels.hpp @@ -188,7 +188,7 @@ class ResidualNormKernel : public physicsSolverBaseKernels::ResidualNormKernelBa m_phaseVolFraction_n( subRegion.getField< fields::well::phaseVolumeFraction_n >()), m_phaseInternalEnergy_n( fluid.phaseInternalEnergy_n() ) { - auto const * rateConstraint = wellControls.getRateConstraint(); + auto const * rateConstraint = wellControls.getRateConstraints().front(); if( m_isProducer ) { // tjb This needs to be fixed should use current constraint rate for normalization From 50f991b8d1e4ddfc0c57c7cc6390d3a2cd777050 Mon Sep 17 00:00:00 2001 From: dkachuma Date: Thu, 21 May 2026 21:40:07 -0500 Subject: [PATCH 08/12] Implement constraint getters --- .../fluidFlow/wells/WellControls.cpp | 101 +++++++++++++++--- .../wells/WellInjectionConstraint.cpp | 6 +- 2 files changed, 92 insertions(+), 15 deletions(-) diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/WellControls.cpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/WellControls.cpp index bb9dda95d69..b399eb04877 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/WellControls.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/WellControls.cpp @@ -437,6 +437,89 @@ void WellControls::setNextDtFromTable( TableFunction const * table, real64 const } } +namespace +{ +template< typename GROUP, + typename CONSTRAINT = std::conditional_t< + std::is_const_v< std::remove_reference_t< GROUP > >, + WellConstraintBase const, + WellConstraintBase > > +void populateConstraints( GROUP & group, bool isProducer, stdVector< CONSTRAINT * > & constraints ) +{ + if( isProducer ) + { + group.template forSubGroups< ProductionConstraint< MassRateConstraint >, + ProductionConstraint< VolumeRateConstraint >, + ProductionConstraint< PhaseVolumeRateConstraint >, + ProductionConstraint< LiquidRateConstraint > >( [&]( CONSTRAINT & constraint ) + { + constraints.push_back( &constraint ); + } ); + } + else + { + group.template forSubGroups< InjectionConstraint< MassRateConstraint >, + InjectionConstraint< VolumeRateConstraint >, + InjectionConstraint< PhaseVolumeRateConstraint >, + InjectionConstraint< LiquidRateConstraint > >( [&]( CONSTRAINT & constraint ) + { + constraints.push_back( &constraint ); + } ); + } +} +} + +WellConstraintBase const * WellControls::getBHPConstraint() const +{ + WellConstraintBase const * bhpConstraint = nullptr; + // Rely on validation here. We assume that there aren't both constraints listed + forSubGroups< MinimumBHPConstraint, + MaximumBHPConstraint >( [&]( WellConstraintBase const & constraint ) + { + bhpConstraint = &constraint; + } ); + return bhpConstraint; +} + +WellConstraintBase * WellControls::getBHPConstraint() +{ + WellConstraintBase * bhpConstraint = nullptr; + // Rely on validation here. We assume that there aren't both constraints listed + forSubGroups< MinimumBHPConstraint, + MaximumBHPConstraint >( [&]( WellConstraintBase & constraint ) + { + bhpConstraint = &constraint; + } ); + return bhpConstraint; +} + +stdVector< WellConstraintBase const * > WellControls::getRateConstraints() const +{ + stdVector< WellConstraintBase const * > constraints; + populateConstraints( *this, isProducer(), constraints ); + return constraints; +} + +stdVector< WellConstraintBase * > WellControls::getRateConstraints() +{ + stdVector< WellConstraintBase * > constraints; + populateConstraints( *this, isProducer(), constraints ); + return constraints; +} + +stdVector< WellConstraintBase const * > WellControls::getAllConstraints() const +{ + stdVector< WellConstraintBase const * > constraints = getRateConstraints(); + constraints.insert( constraints.begin(), getBHPConstraint() ); + return constraints; +} +stdVector< WellConstraintBase * > WellControls::getAllConstraints() +{ + stdVector< WellConstraintBase * > constraints = getRateConstraints(); + constraints.insert( constraints.begin(), getBHPConstraint() ); + return constraints; +} + real64 WellControls::getTargetBHP( real64 const & targetTime ) const { return getBHPConstraint()->getConstraintValue( targetTime ); @@ -451,13 +534,10 @@ real64 WellControls::getInjectionTemperature() const InjectionConstraint< PhaseVolumeRateConstraint >, InjectionConstraint< LiquidRateConstraint > >( [&] ( localIndex index, auto const & constraint ) { - if( firstIndex < 0 ) + if( firstIndex < 0 && constraint.isConstraintActive() ) { - if( constraint.isConstraintActive() ) - { - injectionTemperature = constraint.getInjectionTemperature(); - firstIndex = index; - } + injectionTemperature = constraint.getInjectionTemperature(); + firstIndex = index; } } ); return injectionTemperature; @@ -472,13 +552,10 @@ arrayView1d< real64 const > WellControls::getInjectionStream() const InjectionConstraint< PhaseVolumeRateConstraint >, InjectionConstraint< LiquidRateConstraint > >( [&] ( localIndex index, auto const & constraint ) { - if( firstIndex < 0 ) + if( firstIndex < 0 && constraint.isConstraintActive() ) { - if( constraint.isConstraintActive()) - { - injectionStream = constraint.getInjectionStream(); - firstIndex = index; - } + injectionStream = constraint.getInjectionStream(); + firstIndex = index; } } ); return injectionStream; diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/WellInjectionConstraint.cpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/WellInjectionConstraint.cpp index 57145a752dc..12381d37fee 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/WellInjectionConstraint.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/WellInjectionConstraint.cpp @@ -88,9 +88,9 @@ void InjectionConstraint< ConstraintRateType >::validateInjectionStream( ) } // Register concrete wrapper constraint types and instantiate templates. -template class InjectionConstraint< LiquidRateConstraint >; -using InjectionLiquidRateConstraint = InjectionConstraint< LiquidRateConstraint >; -REGISTER_CATALOG_ENTRY( WellConstraintBase, InjectionLiquidRateConstraint, string const &, Group * const ) +//template class InjectionConstraint< LiquidRateConstraint >; +//using InjectionLiquidRateConstraint = InjectionConstraint< LiquidRateConstraint >; +//REGISTER_CATALOG_ENTRY( WellConstraintBase, InjectionLiquidRateConstraint, string const &, Group * const ) template class InjectionConstraint< MassRateConstraint >; using InjectionMassRateConstraint = InjectionConstraint< MassRateConstraint >; From 5d5df2d9d706801afc045e7e3634f15cd1fcd2d8 Mon Sep 17 00:00:00 2001 From: dkachuma Date: Fri, 22 May 2026 10:41:01 -0500 Subject: [PATCH 09/12] Reference regions --- .../wells/CompositionalMultiphaseWell.cpp | 56 ++---- .../fluidFlow/wells/SinglePhaseWell.cpp | 21 +-- .../fluidFlow/wells/WellControls.cpp | 170 ++++++++++++++++-- .../fluidFlow/wells/WellControls.hpp | 9 + 4 files changed, 182 insertions(+), 74 deletions(-) diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.cpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.cpp index ba8834b02ce..fc992716ad2 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.cpp @@ -368,28 +368,11 @@ void CompositionalMultiphaseWell::validateWellConstraints( real64 const & time_n if( !useSurfaceConditions() ) { bool const useSeg = getReferenceReservoirRegion().empty(); - GEOS_WARNING_IF( useSeg, - GEOS_FMT( "WellControls {} not set and well constraint fluid property calculations will use " - "top segement pressure and temp ", - WellControls::viewKeyStruct::referenceReservoirRegionString() ) ); if( useSeg ) { setRegionAveragePressure( -1 ); setRegionAverageTemperature( -1 ); } - else - { - // Check if region name exists in list of Reservoir's target regions - string const regionName = getReferenceReservoirRegion(); - CompositionalMultiphaseBase const & flowSolver = getParent().getParent().getGroup< CompositionalMultiphaseBase >( getFlowSolverName() ); - string_array const & targetRegionsNames = flowSolver.getTargetRegionNames(); - auto const pos = std::find( targetRegionsNames.begin(), targetRegionsNames.end(), regionName ); - GEOS_ERROR_IF( pos == targetRegionsNames.end(), - GEOS_FMT( "{}: Region {} is not a target of the reservoir solver and cannot be used for referenceReservoirRegion in WellControl {}.", - getDataContext(), regionName, getName() ) ); - - - } } string const & fluidName = subRegion.getReference< string >( viewKeyStruct::fluidNamesString()); MultiFluidBase const & fluid = subRegion.getConstitutiveModel< MultiFluidBase >( fluidName ); @@ -398,21 +381,16 @@ void CompositionalMultiphaseWell::validateWellConstraints( real64 const & time_n { constraint.validatePhaseType( fluid ); } ); - - - } void CompositionalMultiphaseWell::initializePostSubGroups() { WellControls::initializePostSubGroups(); - } void CompositionalMultiphaseWell::initializePostInitialConditionsPreSubGroups() { WellControls::initializePostInitialConditionsPreSubGroups(); - } void CompositionalMultiphaseWell::initializeWellPostInitialConditionsPreSubGroups( WellElementSubRegion & subRegion ) @@ -776,19 +754,11 @@ void CompositionalMultiphaseWell::updateSeparator( ElementRegionManager const & { if( !getReferenceReservoirRegion().empty() ) { - ElementRegionBase const & region = elemManager.getRegion( getReferenceReservoirRegion() ); - GEOS_ERROR_IF ( !region.hasWrapper( CompositionalMultiphaseStatistics::regionStatisticsName()), - GEOS_FMT( "{}: WellControl {} referenceReservoirRegion field requires CompositionalMultiphaseStatistics to be configured for region {} ", - getDataContext(), getName(), getReferenceReservoirRegion() ) ); - - CompositionalMultiphaseStatistics::RegionStatistics const & stats = region.getReference< CompositionalMultiphaseStatistics::RegionStatistics >( - CompositionalMultiphaseStatistics::regionStatisticsName() ); - GEOS_ERROR_IF( stats.averagePressure <= 0.0, - GEOS_FMT( - "{}: No region average quantities computed. WellControl {} referenceReservoirRegion field requires CompositionalMultiphaseStatistics to be configured for region {} ", - getDataContext(), getName(), getReferenceReservoirRegion() )); - setRegionAveragePressure( stats.averagePressure ); - setRegionAverageTemperature( stats.averageTemperature ); + real64 averagePressure = 0.0; + real64 averageTemperature = 0.0; + validateReferenceRegionStatistics< CompositionalMultiphaseStatistics >( elemManager, averagePressure, averageTemperature ); + setRegionAveragePressure( averagePressure ); + setRegionAverageTemperature( averageTemperature ); } // If flashPressure is not set by region the value is defaulted to -1 and indicates to use top segment conditions flashPressure = getRegionAveragePressure(); @@ -1130,10 +1100,10 @@ void CompositionalMultiphaseWell::initializeWell( DomainPartition & domain, Mesh { if( isProducer() ) { - forSubGroups< MinimumBHPConstraint, ProductionConstraint< VolumeRateConstraint >, ProductionConstraint< MassRateConstraint >, - ProductionConstraint< PhaseVolumeRateConstraint > >( [&]( - auto - & constraint ) + forSubGroups< MinimumBHPConstraint, + ProductionConstraint< VolumeRateConstraint >, + ProductionConstraint< MassRateConstraint >, + ProductionConstraint< PhaseVolumeRateConstraint > >( [&]( auto & constraint ) { if( ConstraintTypeId( getControl()) == constraint.getControl() ) { @@ -1143,10 +1113,10 @@ void CompositionalMultiphaseWell::initializeWell( DomainPartition & domain, Mesh } else { - forSubGroups< MaximumBHPConstraint, InjectionConstraint< VolumeRateConstraint >, InjectionConstraint< MassRateConstraint >, InjectionConstraint< PhaseVolumeRateConstraint > >( [&]( - auto - & - constraint ) + forSubGroups< MaximumBHPConstraint, + InjectionConstraint< VolumeRateConstraint >, + InjectionConstraint< MassRateConstraint >, + InjectionConstraint< PhaseVolumeRateConstraint > >( [&]( auto & constraint ) { if( ConstraintTypeId( getControl()) == constraint.getControl() ) { diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/SinglePhaseWell.cpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/SinglePhaseWell.cpp index 72b95c20ae2..aa0284422f4 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/SinglePhaseWell.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/SinglePhaseWell.cpp @@ -345,22 +345,11 @@ void SinglePhaseWell::updateSeparator( ElementRegionManager const & elemManager, { if( !getReferenceReservoirRegion().empty() ) { - ElementRegionBase const & region = elemManager.getRegion( getReferenceReservoirRegion() ); - GEOS_ERROR_IF ( !region.hasWrapper( SinglePhaseStatistics::regionStatisticsName()), - GEOS_FMT( "WellControl {} referenceReservoirRegion field requires SinglePhaseStatistics to be configured for region {} ", - getName(), getReferenceReservoirRegion() ), - getDataContext() ); - - SinglePhaseStatistics::RegionStatistics const & stats = region.getReference< SinglePhaseStatistics::RegionStatistics >( - SinglePhaseStatistics::regionStatisticsName() ); - - GEOS_ERROR_IF( stats.averagePressure <= 0.0, - GEOS_FMT( - "No region average quantities computed. WellControl {} referenceReservoirRegion field requires SinglePhaseStatistics to be configured for region {} ", - getName(), getReferenceReservoirRegion() ), - getDataContext()); - setRegionAveragePressure( stats.averagePressure ); - setRegionAverageTemperature( stats.averageTemperature ); + real64 averagePressure = 0.0; + real64 averageTemperature = 0.0; + validateReferenceRegionStatistics< SinglePhaseStatistics >( elemManager, averagePressure, averageTemperature ); + setRegionAveragePressure( averagePressure ); + setRegionAverageTemperature( averageTemperature ); } // use region conditions flashPressure = getRegionAveragePressure(); diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/WellControls.cpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/WellControls.cpp index b399eb04877..a2d8a39f6c0 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/WellControls.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/WellControls.cpp @@ -27,6 +27,10 @@ #include "physicsSolvers/fluidFlow/wells/WellMassRateConstraint.hpp" #include "physicsSolvers/fluidFlow/wells/WellLiquidRateConstraint.hpp" +#include "physicsSolvers/fluidFlow/FlowSolverBase.hpp" +#include "physicsSolvers/fluidFlow/CompositionalMultiphaseStatistics.hpp" +#include "physicsSolvers/fluidFlow/SinglePhaseStatistics.hpp" + #include "LogLevelsInfo.hpp" #include "WellConstants.hpp" #include "dataRepository/InputFlags.hpp" @@ -301,21 +305,90 @@ void WellControls::postInputInitialization() InputError, getDataContext() ); } -//## // 13) Validate constraints -//## bool const producer = isProducer(); -//## ((void)producer); -//## mapBase typeCount; -//## for( auto & catalog : WellConstraintBase::getCatalog()) -//## { -//## typeCount[catalog.first] = 0; -//## } -//## forSubGroups< -//## >([&](auto const & group) -//## { -//## string const catalogName = group.getCatalogName(); -//## ((void)catalogName); -//##//if (typeCount. in ) -//## }); + // 13) Validate constraints + bool const isProducerWell = isProducer(); + stdVector< std::tuple< string, string, WellConstraintBase const * > > constraints; + forSubGroups< MaximumBHPConstraint, + InjectionConstraint< MassRateConstraint >, + InjectionConstraint< VolumeRateConstraint >, + InjectionConstraint< PhaseVolumeRateConstraint >, + InjectionConstraint< LiquidRateConstraint >, + MinimumBHPConstraint, + ProductionConstraint< MassRateConstraint >, + ProductionConstraint< VolumeRateConstraint >, + ProductionConstraint< PhaseVolumeRateConstraint >, + ProductionConstraint< LiquidRateConstraint > >( [&]( auto const & constraint ) + { + using ConstraintType = std::decay_t< decltype(constraint) >; + constraints.emplace_back( constraint.getName(), ConstraintType::catalogName(), &constraint ); + } ); + + // 13.1) Make sure a producer does not have injector constraints and vice versa + const std::set< string > types = [isProducerWell]() -> std::set< string > + { + if( isProducerWell ) + { + return { MaximumBHPConstraint::catalogName(), + InjectionConstraint< MassRateConstraint >::catalogName(), + InjectionConstraint< VolumeRateConstraint >::catalogName(), + InjectionConstraint< PhaseVolumeRateConstraint >::catalogName(), + InjectionConstraint< LiquidRateConstraint >::catalogName() }; + } + else + { + return { MinimumBHPConstraint::catalogName(), + ProductionConstraint< MassRateConstraint >::catalogName(), + ProductionConstraint< VolumeRateConstraint >::catalogName(), + ProductionConstraint< PhaseVolumeRateConstraint >::catalogName(), + ProductionConstraint< LiquidRateConstraint >::catalogName() }; + } + }(); + for( const auto & [name, type, constraint] : constraints ) + { + GEOS_THROW_IF( types.find( type ) != types.end(), + GEOS_FMT( "Constraint {} of type {} is not allowed for {} wells", + name, type, (isProducerWell ? "producer" : "injector")), + InputError, constraint->getDataContext() ); + } + + // 13.2) Make sure we don't have multiple constraints of the same type + // Track the types we have already seen + mapBase< string, WellConstraintBase const *, std::false_type > seen_types; + for( const auto & [name, type, constraint] : constraints ) + { + auto [it, inserted] = seen_types.insert( {type, constraint} ); + GEOS_THROW_IF( !inserted, + GEOS_FMT( "Constraint of type {} is duplicated by {} and {}", + type, name, it->second->getName() ), + InputError, constraint->getDataContext() ); + } + + // 13.3) Make sure there is a BHP constraint + string const bhp_type = isProducerWell ? MinimumBHPConstraint::catalogName() : MaximumBHPConstraint::catalogName(); + bool const no_match_found = std::none_of( constraints.begin(), constraints.end(), [&bhp_type]( const auto & constraint_tuple ) + { + return std::get< 1 >( constraint_tuple ) == bhp_type; + } ); + GEOS_THROW_IF( no_match_found, + GEOS_FMT( "Constraint of type {} is missing and is required for a {} well", + bhp_type, (isProducerWell ? "producer" : "injector") ), + InputError, getDataContext() ); + + // 13.4) Make sure there is at least one non-BHP constraint + bool const rate_match_found = std::any_of( constraints.begin(), constraints.end(), [&bhp_type]( const auto & constraint_tuple ) + { + return std::get< 1 >( constraint_tuple ) != bhp_type; + } ); + GEOS_THROW_IF( !rate_match_found, + GEOS_FMT( "Missing rate constraint for well {}", + getName(), (isProducerWell ? "producer" : "injector") ), + InputError, getDataContext() ); +} + +void WellControls::initializePreSubGroups() +{ + // Validate the reference region + validateReferenceRegion(); } void WellControls::postRestartInitialization( ) @@ -1006,6 +1079,73 @@ WellControls::assembleSystem( real64 const & time_n, assembleWellPressureRelations( time_n, dt, subRegion, dofManager, localMatrix.toViewConstSizes(), localRhs ); computeWellPerforationRates( time_n, dt, elemManager, subRegion ); assembleWellFluxTerms( time_n, dt, subRegion, dofManager, localMatrix.toViewConstSizes(), localRhs ); +} + +bool WellControls::validateReferenceRegion() const +{ + // If using surface conditions then there is nothing to validate + if( useSurfaceConditions()) + { + return true; + } + bool const isRoot = MpiWrapper::commRank() == 0; + string const regionName = getReferenceReservoirRegion(); + if( regionName.empty() ) + { + GEOS_WARNING_IF( isRoot, + GEOS_FMT( "WellControls {} referenceReservoirRegion not set and well constraint " + "fluid property calculations will use top segement pressure and temperature.", + getName()) ); + } + else + { + FlowSolverBase const * flowSolver = getParent().getParent().getGroupPointer< FlowSolverBase >( getFlowSolverName() ); + if( flowSolver == nullptr ) + { + return true; + } + string_array const & targetRegionsNames = flowSolver->getTargetRegionNames(); + auto const pos = std::find( targetRegionsNames.begin(), targetRegionsNames.end(), regionName ); + GEOS_THROW_IF( pos == targetRegionsNames.end(), + GEOS_FMT( "Region {} is not a target of the reservoir solver {} and cannot " + "be used for referenceReservoirRegion in WellControl {}.", + regionName, flowSolver->getName(), getName() ), + InputError, getDataContext() ); + + return pos != targetRegionsNames.end(); + } + return true; +} +template< typename STATISTICS > +bool WellControls::validateReferenceRegionStatistics( ElementRegionManager const & elemManager, + real64 & averagePressure, + real64 & averageTemperature ) const +{ + averagePressure = 0.0; + averageTemperature = 0.0; + string const regionName = getReferenceReservoirRegion(); + if( !regionName.empty()) + { + ElementRegionBase const & region = elemManager.getRegion( regionName ); + GEOS_THROW_IF( !region.hasWrapper( STATISTICS::regionStatisticsName() ), + GEOS_FMT( "WellControl {} referenceReservoirRegion field requires {} to be configured for region {}", + getName(), STATISTICS::catalogName(), regionName ), + RuntimeError, getDataContext() ); + + auto const & stats = region.getReference< typename STATISTICS::RegionStatistics >( STATISTICS::regionStatisticsName() ); + GEOS_THROW_IF( stats.averagePressure <= 0.0, + GEOS_FMT( "No region average quantities computed. WellControl {} referenceReservoirRegion field requires " + "{} to be configured for region {} ", + getName(), STATISTICS::catalogName(), regionName ), + RuntimeError, getDataContext()); + averagePressure = stats.averagePressure; + averageTemperature = stats.averageTemperature; + } + return true; } + +template bool WellControls::validateReferenceRegionStatistics< SinglePhaseStatistics >( ElementRegionManager const &, real64 &, real64 & ) const; +template bool WellControls::validateReferenceRegionStatistics< CompositionalMultiphaseStatistics >( ElementRegionManager const &, real64 &, real64 & ) const; + } //namespace geos diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/WellControls.hpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/WellControls.hpp index f5e529862c1..adad3a1ebba 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/WellControls.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/WellControls.hpp @@ -376,6 +376,8 @@ class WellControls : public dataRepository::Group virtual void postInputInitialization() override; + virtual void initializePreSubGroups() override; + virtual void initializeWellPostInitialConditionsPreSubGroups( WellElementSubRegion & subRegion ) = 0; virtual void printRates( real64 const & time_n, real64 const & dt, @@ -810,6 +812,13 @@ class WellControls : public dataRepository::Group void logConstraint( WellConstraintBase const * constraint, WellElementSubRegion const & region, real64 time, bool isLimiting = false ) const; + bool validateReferenceRegion() const; + + template< typename STATISTICS > + bool validateReferenceRegionStatistics( ElementRegionManager const & elemManager, + real64 & averagePressure, + real64 & averageTemperature ) const; + private: /// List of names of regions the solver will be applied to string_array m_targetRegionNames; From c5f79e75f3201e9eded32ae4b9623303fe57128f Mon Sep 17 00:00:00 2001 From: dkachuma Date: Fri, 22 May 2026 11:10:58 -0500 Subject: [PATCH 10/12] Some comments --- .../fluidFlow/wells/WellControls.hpp | 64 ++++++++++++++++++- 1 file changed, 62 insertions(+), 2 deletions(-) diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/WellControls.hpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/WellControls.hpp index adad3a1ebba..2bebb6ec67a 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/WellControls.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/WellControls.hpp @@ -767,12 +767,31 @@ class WellControls : public dataRepository::Group */ template< typename ConstraintType > void createConstraint ( string const & constraintName ); + /** + * @brief Gets the defined BHP constraint + * @details Returns the BHP constraint if one is defined for the WellControl. For a producer + * well this will be a minimum BHP constraint and for an injector well this will be a maximum + * BHP constraint. This will possibly return null if no BHP constraint is set. Validation is + * in place to enforce the setting of at least one BHp constraint. + * @return A BHP constraint object of one is defined + */ WellConstraintBase const * getBHPConstraint() const; WellConstraintBase * getBHPConstraint(); + /** + * @brief Gets a list of rate constraints + * @details Returns a list of rate constraints for the WellControl. For a producer + * well these will be a production rate constraints `ProductionConstraint` and for an + * injector well these will be injection rate constraints `InjectionConstraint`. + */ stdVector< WellConstraintBase const * > getRateConstraints() const; stdVector< WellConstraintBase * > getRateConstraints(); + /** + * @brief Gets a list of all constraints constraints + * @details Returns a list of all constraints for the WellControl including rate and BHP + * constraints. + */ stdVector< WellConstraintBase const * > getAllConstraints() const; stdVector< WellConstraintBase * > getAllConstraints(); @@ -810,12 +829,53 @@ class WellControls : public dataRepository::Group protected: virtual void postRestartInitialization( )override; - void logConstraint( WellConstraintBase const * constraint, WellElementSubRegion const & region, real64 time, bool isLimiting = false ) const; + /** + * @brief Logs the state and values of a specific well constraint. + * + * @details This method evaluates whether the provided constraint requires + * logging based on its validity, the current log level being strictly greater + * than 4, and the region being locally owned. When the constraint is flagged + * as the limiting constraint, it logs extensive operational data such as the + * bottom hole pressure, individual phase volume rates, the total volume rate, + * and the mass rate. Otherwise, it logs general information, specifically + * whether the constraint is active and its target value at the current time. + * + * @param constraint Pointer to the base constraint object to evaluate and log. + * @param region Reference to the well element sub-region associated with it. + * @param time The current simulation time used to get the constraint value. + * @param isLimiting Boolean indicating if this is the active limiting constraint. + */ + void logConstraint( WellConstraintBase const * constraint, + WellElementSubRegion const & region, + real64 time, + bool isLimiting = false ) const; + /** + * @brief Validates the reference region + * @details Validates the reference region by ensuring that it is defined for cases that do not + * use surface conditions. If the reference region is provided, it is checked against the flow + * solves regions. + * @return @c true if the region is valid + */ bool validateReferenceRegion() const; + /** + * @brief Validates and retrieves the reference region statistics for average pressure and temperature. + * + * @details This template method checks if a reference reservoir region is configured for the well control. + * If a region is specified, it retrieves the region from the provided element manager and verifies + * that the required statistics wrapper exists. It then extracts the average pressure and temperature, + * throwing an exception if the average pressure has not been properly computed. + * + * @tparam STATISTICS Type providing static methods and types for region statistics (SinglePhaseStatistics or + * CompositionalMultiphaseStatistics). + * @param elementManager Reference to the ElementRegionManager used to look up the reservoir region. + * @param[out] averagePressure Reference to a real64 variable where the retrieved average pressure is stored. + * @param[out] averageTemperature Reference to a real64 variable where the retrieved average temperature is stored. + * @return Boolean value, always returning true upon successful validation. + */ template< typename STATISTICS > - bool validateReferenceRegionStatistics( ElementRegionManager const & elemManager, + bool validateReferenceRegionStatistics( ElementRegionManager const & elementManager, real64 & averagePressure, real64 & averageTemperature ) const; From 44b30d4819996c27b84cbd2a28a7f64fb9170065 Mon Sep 17 00:00:00 2001 From: Dickson Kachuma <81433670+dkachuma@users.noreply.github.com> Date: Fri, 22 May 2026 12:08:45 -0500 Subject: [PATCH 11/12] Potential fix for pull request finding Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com> --- .../physicsSolvers/fluidFlow/wells/WellControls.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/WellControls.cpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/WellControls.cpp index a2d8a39f6c0..eaab2a4736f 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/WellControls.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/WellControls.cpp @@ -380,8 +380,8 @@ void WellControls::postInputInitialization() return std::get< 1 >( constraint_tuple ) != bhp_type; } ); GEOS_THROW_IF( !rate_match_found, - GEOS_FMT( "Missing rate constraint for well {}", - getName(), (isProducerWell ? "producer" : "injector") ), + GEOS_FMT( "Missing rate constraint for {} well {}", + (isProducerWell ? "producer" : "injector"), getName() ), InputError, getDataContext() ); } From 6e653a7b2d28c7df00d8dda1a6e9921947f85b32 Mon Sep 17 00:00:00 2001 From: Dickson Kachuma <81433670+dkachuma@users.noreply.github.com> Date: Fri, 22 May 2026 12:09:17 -0500 Subject: [PATCH 12/12] Potential fix for pull request finding Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com> --- .../physicsSolvers/fluidFlow/wells/WellControls.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/WellControls.cpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/WellControls.cpp index eaab2a4736f..38d2457e4e7 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/WellControls.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/WellControls.cpp @@ -1094,7 +1094,7 @@ bool WellControls::validateReferenceRegion() const { GEOS_WARNING_IF( isRoot, GEOS_FMT( "WellControls {} referenceReservoirRegion not set and well constraint " - "fluid property calculations will use top segement pressure and temperature.", + "fluid property calculations will use top segment pressure and temperature.", getName()) ); } else