diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.cpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.cpp index 0914b977dfd..fc992716ad2 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() ). @@ -360,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 ); @@ -390,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 ) @@ -768,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(); @@ -1122,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() ) { @@ -1135,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 a2172b71792..aa0284422f4 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" @@ -336,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/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..38d2457e4e7 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/WellControls.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/WellControls.cpp @@ -17,8 +17,21 @@ * @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 "physicsSolvers/fluidFlow/FlowSolverBase.hpp" +#include "physicsSolvers/fluidFlow/CompositionalMultiphaseStatistics.hpp" +#include "physicsSolvers/fluidFlow/SinglePhaseStatistics.hpp" + +#include "LogLevelsInfo.hpp" #include "WellConstants.hpp" #include "dataRepository/InputFlags.hpp" #include "functions/FunctionManager.hpp" @@ -27,6 +40,7 @@ #include "physicsSolvers/fluidFlow/wells/WellFields.hpp" #include "functions/FunctionManager.hpp" + namespace geos { @@ -41,7 +55,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 +71,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 +156,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 >(); } @@ -155,70 +167,16 @@ 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; + 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 ConstitutiveBase here - for( auto & catalogIter: WellConstraintBase::getCatalog()) + // During schema generation, register one of each type derived from WellConstraintBase here + for( auto & catalogIter : WellConstraintBase::getCatalog()) { createChild( catalogIter.first, catalogIter.first ); } @@ -247,6 +205,7 @@ TableFunction * createWellTable( string const & tableName, } } + void WellControls::registerWellDataOnMesh( WellElementSubRegion & subRegion ) { std::string const & regionName = subRegion.getName(); @@ -276,6 +235,7 @@ void WellControls::registerWellDataOnMesh( WellElementSubRegion & subRegion ) } } } + void WellControls::postInputInitialization() { Group::postInputInitialization(); @@ -345,48 +305,145 @@ void WellControls::postInputInitialization() InputError, getDataContext() ); } + // 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 {}", + (isProducerWell ? "producer" : "injector"), getName() ), + InputError, getDataContext() ); +} + +void WellControls::initializePreSubGroups() +{ + // Validate the reference region + validateReferenceRegion(); } + void WellControls::postRestartInitialization( ) {} + +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 ) + { + 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 ) ) ); + } +} + 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 ) { - std::vector< WellConstraintBase * > const constraints = getProdRateConstraints(); - for( auto const & constraint : constraints ) - { - if( isZero( constraint->getConstraintValue( currentTime ) ) ) - { - m_wellStatus = WellControls::Status::CLOSED; - m_currentConstraint=nullptr; - break; - } - } + hasZeroRate = true; } - else + + if( !hasZeroRate ) { - std::vector< WellConstraintBase * > const constraints = getInjRateConstraints(); - for( auto const & constraint : constraints ) + for( auto const * constraint : getRateConstraints() ) { 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 +470,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; @@ -430,27 +488,12 @@ bool WellControls::getWellState() const void WellControls::setNextDtFromTables( real64 const & currentTime, real64 & nextDt ) { - if( isProducer() ) + for( auto const * constraint : getAllConstraints() ) { - if( getMinBHPConstraint() != nullptr ) - { - getMinBHPConstraint()->setNextDtFromTables( currentTime, nextDt ); - } - for( auto const & constraint : m_productionRateConstraintList ) - { - constraint->setNextDtFromTables( currentTime, nextDt ); - } - } - else - { - getMaxBHPConstraint()->setNextDtFromTables( currentTime, nextDt ); - for( auto const & constraint : m_injectionRateConstraintList ) - { - 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 ) @@ -467,81 +510,157 @@ void WellControls::setNextDtFromTable( TableFunction const * table, real64 const } } -real64 WellControls::getTargetBHP( real64 const & targetTime ) 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()) + if( isProducer ) + { + group.template forSubGroups< ProductionConstraint< MassRateConstraint >, + ProductionConstraint< VolumeRateConstraint >, + ProductionConstraint< PhaseVolumeRateConstraint >, + ProductionConstraint< LiquidRateConstraint > >( [&]( CONSTRAINT & constraint ) + { + constraints.push_back( &constraint ); + } ); + } + else { - return m_minBHPConstraint->getConstraintValue( targetTime ); + group.template forSubGroups< InjectionConstraint< MassRateConstraint >, + InjectionConstraint< VolumeRateConstraint >, + InjectionConstraint< PhaseVolumeRateConstraint >, + InjectionConstraint< LiquidRateConstraint > >( [&]( CONSTRAINT & constraint ) + { + constraints.push_back( &constraint ); + } ); } - return m_maxBHPConstraint->getConstraintValue( targetTime ); +} +} + +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 ); +} real64 WellControls::getInjectionTemperature() const { real64 injectionTemperature = 0.0; - this->forInjectionConstraints< InjectionConstraint< PhaseVolumeRateConstraint >, InjectionConstraint< VolumeRateConstraint >, InjectionConstraint< MassRateConstraint > >( [&] ( auto & 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( constraint.isConstraintActive()) + if( firstIndex < 0 && constraint.isConstraintActive() ) { - injectionTemperature = constraint.getInjectionTemperature(); - return; + injectionTemperature = constraint.getInjectionTemperature(); + firstIndex = index; } } ); return injectionTemperature; } - arrayView1d< real64 const > WellControls::getInjectionStream() const { arrayView1d< real64 const > injectionStream; - forInjectionConstraints< InjectionConstraint< PhaseVolumeRateConstraint >, InjectionConstraint< VolumeRateConstraint >, InjectionConstraint< MassRateConstraint > >( [&] ( auto & 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( constraint.isConstraintActive() ) + if( firstIndex < 0 && constraint.isConstraintActive() ) { injectionStream = constraint.getInjectionStream(); - return; + firstIndex = index; } } ); - return injectionStream; } integer WellControls::getConstraintPhaseIndex() const { integer phaseIndex = -1; - - 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< ProductionConstraint< PhaseVolumeRateConstraint >, + InjectionConstraint< PhaseVolumeRateConstraint > >( [&] ( auto & constraint ) { - forProductionConstraints< ProductionConstraint< PhaseVolumeRateConstraint > >( [&] ( auto & constraint ) - { - if( constraint.isConstraintActive() ) - { - phaseIndex = constraint.getPhaseIndex(); - } - } ); - } - else - { - forInjectionConstraints< InjectionConstraint< PhaseVolumeRateConstraint > >( [&] ( auto & constraint ) + if( constraint.isConstraintActive() ) { - if( constraint.isConstraintActive() ) - { - phaseIndex = constraint.getPhaseIndex(); - } - } ); - } - + phaseIndex = constraint.getPhaseIndex(); + } + } ); return phaseIndex; } real64 WellControls::getReferenceElevation() const { - if( isProducer () ) + real64 referenceElevation = 0.0; + // 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 ) { - return getMinBHPConstraint()->getReferenceElevation(); - } - return getMaxBHPConstraint()->getReferenceElevation(); + referenceElevation = constraint.getReferenceElevation(); + } ); + return referenceElevation; } void WellControls::implicitStepSetup( real64 const & time_n, @@ -554,6 +673,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(); @@ -681,8 +801,6 @@ void WellControls::selectWellConstraint( real64 const & time_n, bool useEstimator = coupledIterationNumber < estimateSolution(); - - if( getWellState()) { if( useEstimator ) @@ -718,8 +836,6 @@ void WellControls::selectWellConstraint( real64 const & time_n, MpiWrapper::broadcast( wellControl, owner ); setControl( wellControl ); } - - } bool WellControls::evaluateConstraints( real64 const & time_n, @@ -727,30 +843,16 @@ bool WellControls::evaluateConstraints( real64 const & time_n, { // create list of all constraints to process - std::vector< WellConstraintBase * > constraintList; - if( isProducer() ) - { - constraintList = getProdRateConstraints(); - // Solve minimum bhp constraint first - constraintList.insert( constraintList.begin(), getMinBHPConstraint() ); - } - else - { - constraintList = getInjRateConstraints(); - // Solve maximum bhp constraint first; - constraintList.insert( constraintList.begin(), getMaxBHPConstraint() ); - } + stdVector< WellConstraintBase * > constraintList = getAllConstraints(); + // 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 +873,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 +888,13 @@ 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 = getRateConstraints(); + WellConstraintBase * limitingConstraint = getCurrentConstraint(); if( isProducer() ) { - constraintList = getProdRateConstraints(); if( limitingConstraint->getControl() != ConstraintTypeId::BHP ) { { // remove from list and add BHP constraint @@ -807,38 +903,24 @@ 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 + if( limitingConstraint->getControl() != ConstraintTypeId::BHP ) { - 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( getBHPConstraint() ); } } if( isoThermalEstimatorEnabled() ) @@ -870,7 +952,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 @@ -899,9 +981,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; } @@ -939,7 +1019,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 " << @@ -1000,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 segment 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 3efe3dfcaee..2bebb6ec67a 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/WellControls.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/WellControls.hpp @@ -17,22 +17,15 @@ * @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 "constitutive/fluid/multifluid/MultiFluidBase.hpp" #include "constitutive/fluid/singlefluid/SingleFluidBase.hpp" #include "physicsSolvers/fluidFlow/wells/WellNewtonSolver.hpp" @@ -48,12 +41,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,15 +132,16 @@ 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 + * @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 + /// Expand catalog for schema generation virtual void expandObjectCatalogs() override; virtual void registerWellDataOnMesh( WellElementSubRegion & subRegion ); @@ -170,11 +165,11 @@ class WellControls : public dataRepository::Group //public PhysicsSolverBase 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; /** @@ -381,117 +376,13 @@ class WellControls : public dataRepository::Group //public PhysicsSolverBase virtual void postInputInitialization() override; + virtual void initializePreSubGroups() override; + virtual void initializeWellPostInitialConditionsPreSubGroups( WellElementSubRegion & subRegion ) = 0; virtual void printRates( real64 const & time_n, real64 const & dt, WellElementSubRegion const & subRegion ) = 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 ); - } ); - } - } /** * @name Getters / Setters */ @@ -577,7 +468,6 @@ class WellControls : public dataRepository::Group //public PhysicsSolverBase */ real64 getTargetBHP( real64 const & targetTime ) const; - /** * @brief Const accessor for the temperature of the injection stream * @return the temperature of the injection stream @@ -602,7 +492,6 @@ class WellControls : public dataRepository::Group //public PhysicsSolverBase */ 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 @@ -676,7 +565,6 @@ class WellControls : public dataRepository::Group //public PhysicsSolverBase */ bool getWellState() const; - /** * @brief Set the current consrtaint * @param[in] currentConstraint pointer to constraint @@ -686,9 +574,8 @@ class WellControls : public dataRepository::Group //public PhysicsSolverBase * @brief Get the current consrtaint * @return pointer to constraint */ - WellConstraintBase * getCurrentConstraint() { return m_currentConstraint; } - WellConstraintBase const * getCurrentConstraint() const { return m_currentConstraint; } - + WellConstraintBase * getCurrentConstraint() { return m_currentConstraint; } + WellConstraintBase const * getCurrentConstraint() const { return m_currentConstraint; } /** * @brief Getter for the flag to enable crossflow @@ -719,7 +606,6 @@ class WellControls : public dataRepository::Group //public PhysicsSolverBase void setKeepVariablesConstantDuringInitStep( bool const keepVariablesConstantDuringInitStep ) { m_keepVariablesConstantDuringInitStep = keepVariablesConstantDuringInitStep; } - /** * @brief setter for multi fluid separator * @param[in] fluidSeparatorPtr single or multiphase separator @@ -773,12 +659,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 +681,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 +693,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,34 +722,12 @@ 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 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"; } - /// string key for the minimum BHP presssure for a producer - 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"; } @@ -877,10 +738,17 @@ class WellControls : public dataRepository::Group //public PhysicsSolverBase static constexpr char const * currentTotalVolRateString() { return "currentTotalVolumetricRate"; } static constexpr char const * currentMassRateString() { return "currentMassRate"; } + }; + + /** + * @brief Structure to hold scoped key names + */ + struct groupKeyStruct + { + /// string key for the well Newton solver + static constexpr char const * wellNewtonSolverString() { return "WellNewtonSolver"; } + }; - } - /// ViewKey struct for the WellControls class - viewKeysWellControls; void setPerforationStatus( real64 const & time_n, WellElementSubRegion & subRegion ); void setGravCoef( WellElementSubRegion & subRegion, R1Tensor const & gravVector ); /** @@ -899,22 +767,33 @@ class WellControls : public dataRepository::Group //public PhysicsSolverBase */ 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 Getters for constraints + * @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`. */ - 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; }; + stdVector< WellConstraintBase const * > getRateConstraints() const; + stdVector< WellConstraintBase * > getRateConstraints(); /** - * @brief Getters for constraint lists + * @brief Gets a list of all constraints constraints + * @details Returns a list of all constraints for the WellControl including rate and BHP + * constraints. */ - 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; } + stdVector< WellConstraintBase const * > getAllConstraints() const; + stdVector< WellConstraintBase * > getAllConstraints(); /** * @brief Set thermal effects enable @@ -948,8 +827,58 @@ class WellControls : public dataRepository::Group //public PhysicsSolverBase DofManager const & dofManager ); protected: - virtual void postRestartInitialization( )override; + + /** + * @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 & elementManager, + real64 & averagePressure, + real64 & averageTemperature ) const; + private: /// List of names of regions the solver will be applied to string_array m_targetRegionNames; @@ -987,6 +916,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; @@ -1027,15 +957,7 @@ class WellControls : public dataRepository::Group //public PhysicsSolverBase real64 m_initialPressureCoefficient; // Current constrint - WellConstraintBase * m_currentConstraint; - - // 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; + WellConstraintBase * m_currentConstraint{}; /// Well status WellControls::Status m_wellStatus; 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 >; diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/kernels/CompositionalMultiphaseWellKernels.cpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/kernels/CompositionalMultiphaseWellKernels.cpp index 010a9db3e7f..673b4d4dd84 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.getRateConstraints().front(); // 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.getRateConstraints().front(); // 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..2087eae9c9d 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.getRateConstraints().front(); + 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..ab9e50c7646 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.getRateConstraints().front(); + 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..29e97edaaf2 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.getRateConstraints().front(); 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,