add logic for changing between mass rate and bhp constraints

GEOS-DEV · Dec 14, 2024 · b5cec3e · b5cec3e
1 parent b5ea7ce
commit b5cec3e
Show file tree

Hide file tree

Showing 5 changed files with 38 additions and 14 deletions.
diff --git 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,
+                                  &currentMassRate,
                                   &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
@@ -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
@@ -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/...eComponents/physicsSolvers/fluidFlow/wells/kernels/CompositionalMultiphaseWellKernels.cpp b/...eComponents/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/...eComponents/physicsSolvers/fluidFlow/wells/kernels/CompositionalMultiphaseWellKernels.hpp b/...eComponents/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 >