From b5cec3ef0234e06e9ae8b66b279aaa688f5100cb Mon Sep 17 00:00:00 2001 From: Tom Byer <149726499+tjb-ltk@users.noreply.github.com> Date: Fri, 13 Dec 2024 16:06:14 -0800 Subject: [PATCH] add logic for changing between mass rate and bhp constraints --- .../wells/CompositionalMultiphaseWell.cpp | 16 +++++++++++- .../wells/CompositionalMultiphaseWell.hpp | 2 ++ .../fluidFlow/wells/WellControls.hpp | 6 +++++ .../CompositionalMultiphaseWellKernels.cpp | 26 +++++++++---------- .../CompositionalMultiphaseWellKernels.hpp | 2 ++ 5 files changed, 38 insertions(+), 14 deletions(-) diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.cpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.cpp index 1f13044d87..dab39097ae 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.cpp @@ -260,6 +260,8 @@ void CompositionalMultiphaseWell::registerDataOnMesh( Group & meshBodies ) wellControls.registerWrapper< real64 >( viewKeyStruct::massDensityString() ); + wellControls.registerWrapper< real64 >( viewKeyStruct::currentMassRateString() ); + // write rates output header // the rank that owns the reference well element is responsible if( m_writeCSV > 0 && subRegion.isLocallyOwned() ) @@ -686,6 +688,10 @@ void CompositionalMultiphaseWell::updateVolRatesForConstraint( WellElementSubReg real64 & currentTotalVolRate = wellControls.getReference< real64 >( CompositionalMultiphaseWell::viewKeyStruct::currentTotalVolRateString() ); + + real64 & currentMassRate = + wellControls.getReference< real64 >( CompositionalMultiphaseWell::viewKeyStruct::currentMassRateString() ); + arrayView1d< real64 > const & dCurrentTotalVolRate = wellControls.getReference< array1d< real64 > >( CompositionalMultiphaseWell::viewKeyStruct::dCurrentTotalVolRateString() ); @@ -721,6 +727,7 @@ void CompositionalMultiphaseWell::updateVolRatesForConstraint( WellElementSubReg dCurrentTotalVolRate, currentPhaseVolRate, dCurrentPhaseVolRate, + ¤tMassRate, &iwelemRef, &logLevel, &wellControlsName, @@ -757,7 +764,8 @@ void CompositionalMultiphaseWell::updateVolRatesForConstraint( WellElementSubReg // Step 2: update the total volume rate real64 const currentTotalRate = connRate[iwelemRef]; - + // Assumes useMass is true + currentMassRate = currentTotalRate; // Step 2.1: compute the inverse of the total density and derivatives massDensity = totalDens[iwelemRef][0]; real64 const totalDensInv = 1.0 / totalDens[iwelemRef][0]; @@ -1951,6 +1959,12 @@ void CompositionalMultiphaseWell::assemblePressureRelations( real64 const & time GEOS_LOG_LEVEL_INFO_RANK_0( logInfo::WellControl, GEOS_FMT( "Control switch for well {} from BHP constraint to phase volumetric rate constraint", subRegion.getName() ) ); } + else if( wellControls.getInputControl() == WellControls::Control::MASSRATE ) + { + wellControls.switchToMassRateControl( wellControls.getTargetMassRate( timeAtEndOfStep ) ); + GEOS_LOG_LEVEL_INFO_RANK_0( logInfo::WellControl, + GEOS_FMT( "Control switch for well {} from BHP constraint to mass rate constraint", subRegion.getName()) ); + } else { wellControls.switchToTotalRateControl( wellControls.getTargetTotalRate( timeAtEndOfStep ) ); diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.hpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.hpp index aa70d4f450..7261459637 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.hpp @@ -296,6 +296,8 @@ class CompositionalMultiphaseWell : public WellSolverBase static constexpr char const * currentTotalVolRateString() { return "currentTotalVolumetricRate"; } static constexpr char const * dCurrentTotalVolRateString() { return "dCurrentTotalVolumetricRate"; } + static constexpr char const * currentMassRateString() { return "currentMassRate"; } + static constexpr char const * dCurrentTotalVolRate_dPresString() { return "dCurrentTotalVolumetricRate_dPres"; } static constexpr char const * dCurrentTotalVolRate_dCompDensString() { return "dCurrentTotalVolumetricRate_dCompDens"; } diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/WellControls.hpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/WellControls.hpp index 3335b7d98d..9e20209638 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/WellControls.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/WellControls.hpp @@ -149,6 +149,12 @@ class WellControls : public dataRepository::Group */ Control getControl() const { return m_currentControl; } + /** + * @brief Get the input control type for the well. + * @return the Control enum enforced at the well + */ + Control getInputControl() const { return m_inputControl; } + /** * @brief Getter for the reference elevation where the BHP control is enforced * @return the reference elevation diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/kernels/CompositionalMultiphaseWellKernels.cpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/kernels/CompositionalMultiphaseWellKernels.cpp index 42ec691108..01fb94d2e7 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/kernels/CompositionalMultiphaseWellKernels.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/kernels/CompositionalMultiphaseWellKernels.cpp @@ -37,6 +37,7 @@ inline void ControlEquationHelper:: switchControl( bool const isProducer, + WellControls::Control const & inputControl, WellControls::Control const & currentControl, integer const phasePhaseIndex, real64 const & targetBHP, @@ -46,6 +47,7 @@ ControlEquationHelper:: real64 const & currentBHP, arrayView1d< real64 const > const & currentPhaseVolRate, real64 const & currentTotalVolRate, + real64 const & currentMassRate, WellControls::Control & newControl ) { // if isViable is true at the end of the following checks, no need to switch @@ -58,7 +60,7 @@ ControlEquationHelper:: // Currently, the available constraints are: // - Producer: BHP, PHASEVOLRATE - // - Injector: BHP, TOTALVOLRATE + // - Injector: BHP, TOTALVOLRATE, MASSRATE // TODO: support GRAT, WRAT, LIQUID for producers and check if any of the active constraint is violated @@ -71,6 +73,10 @@ ControlEquationHelper:: controlIsViable = ( LvArray::math::abs( currentPhaseVolRate[phasePhaseIndex] ) <= LvArray::math::abs( targetPhaseRate ) ); } // the control is viable if the reference total rate is below the max rate for injectors + else if( inputControl == WellControls::Control::MASSRATE ) + { + controlIsViable = ( LvArray::math::abs( currentMassRate ) <= LvArray::math::abs( targetMassRate ) ); + } else { controlIsViable = ( LvArray::math::abs( currentTotalVolRate ) <= LvArray::math::abs( targetTotalRate ) ); @@ -219,17 +225,6 @@ ControlEquationHelper:: if constexpr ( IS_THERMAL ) dControlEqn[COFFSET_WJ::dT] = massDensity*dCurrentTotalVolRate[COFFSET_WJ::dT]; } - // Total mass rate control - else if( currentControl == WellControls::Control::MASSRATE ) - { - controlEqn = massDensity*currentTotalVolRate - targetMassRate; - dControlEqn[COFFSET_WJ::dP] = massDensity*dCurrentTotalVolRate[COFFSET_WJ::dP]; - dControlEqn[COFFSET_WJ::dQ] = massDensity*dCurrentTotalVolRate[COFFSET_WJ::dQ]; - for( integer ic = 0; ic < NC; ++ic ) - { - dControlEqn[COFFSET_WJ::dC+ic] = massDensity*dCurrentTotalVolRate[COFFSET_WJ::dC+ic]; - } - } else { GEOS_ERROR( "This constraint is not supported in CompositionalMultiphaseWell" ); @@ -323,6 +318,7 @@ PressureRelationKernel:: // static well control data bool const isProducer = wellControls.isProducer(); WellControls::Control const currentControl = wellControls.getControl(); + WellControls::Control const inputControl = wellControls.getInputControl(); real64 const targetBHP = wellControls.getTargetBHP( timeAtEndOfStep ); real64 const targetTotalRate = wellControls.getTargetTotalRate( timeAtEndOfStep ); real64 const targetPhaseRate = wellControls.getTargetPhaseRate( timeAtEndOfStep ); @@ -344,6 +340,9 @@ PressureRelationKernel:: arrayView1d< real64 const > const & dCurrentTotalVolRate = wellControls.getReference< array1d< real64 > >( CompositionalMultiphaseWell::viewKeyStruct::dCurrentTotalVolRateString() ); + real64 const & currentMassRate = + wellControls.getReference< real64 >( CompositionalMultiphaseWell::viewKeyStruct::currentMassRateString() ); + real64 const & massDensity = wellControls.getReference< real64 >( CompositionalMultiphaseWell::viewKeyStruct::massDensityString() ); @@ -358,6 +357,7 @@ PressureRelationKernel:: { WellControls::Control newControl = currentControl; ControlEquationHelper::switchControl( isProducer, + inputControl, currentControl, targetPhaseIndex, targetBHP, @@ -367,6 +367,7 @@ PressureRelationKernel:: currentBHP, currentPhaseVolRate, currentTotalVolRate, + currentMassRate, newControl ); if( currentControl != newControl ) { @@ -737,7 +738,6 @@ RateInitializationKernel:: else if( control == WellControls::Control::MASSRATE ) { connRate[iwelem] = targetMassRate; - connRate[iwelem] = targetMassRate* totalDens[iwelem][0]; } else { diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/kernels/CompositionalMultiphaseWellKernels.hpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/kernels/CompositionalMultiphaseWellKernels.hpp index 41d63148bf..825f58bff2 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/kernels/CompositionalMultiphaseWellKernels.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/kernels/CompositionalMultiphaseWellKernels.hpp @@ -134,6 +134,7 @@ struct ControlEquationHelper static void switchControl( bool const isProducer, + WellControls::Control const & inputControl, WellControls::Control const & currentControl, integer const phasePhaseIndex, real64 const & targetBHP, @@ -143,6 +144,7 @@ struct ControlEquationHelper real64 const & currentBHP, arrayView1d< real64 const > const & currentPhaseVolRate, real64 const & currentTotalVolRate, + real64 const & currentMassRate, WellControls::Control & newControl ); template< integer NC, integer IS_THERMAL >