Feature/byer3/well constraint estimator

src/coreComponents/linearAlgebra/utilities/ComponentMask.hpp

@@ -92,10 +92,10 @@ class ComponentMask
 private:
 
   /// Number of bits in mask storage
-  static constexpr int NUM_BITS = internal::roundToNextPowerOfTwo( MAX_COMP );
+  static constexpr int NUM_BITS = geos::internal::roundToNextPowerOfTwo( MAX_COMP );
 
   /// Type used to represent the bit mask
-  using mask_t = typename internal::ComponentMaskType< NUM_BITS >::type;
+  using mask_t = typename geos::internal::ComponentMaskType< NUM_BITS >::type;
 
 public:
 

src/coreComponents/physicsSolvers/fluidFlow/CMakeLists.txt

@@ -116,11 +116,21 @@ set( fluidFlowSolvers_headers
      wells/SinglePhaseWellFields.hpp
      wells/WellConstants.hpp
      wells/WellControls.hpp
+     wells/WellConstraintsBase.hpp
+     wells/WellInjectionConstraint.hpp
+     wells/WellProductionConstraint.hpp
+     wells/WellBHPConstraints.hpp
+     wells/WellVolumeRateConstraint.hpp
+     wells/WellMassRateConstraint.hpp
+     wells/WellPhaseVolumeRateConstraint.hpp
+     wells/WellLiquidRateConstraint.hpp
      wells/WellSolverBase.hpp
      wells/WellSolverBaseFields.hpp
      wells/LogLevelsInfo.hpp
      wells/kernels/SinglePhaseWellKernels.hpp
      wells/kernels/CompositionalMultiphaseWellKernels.hpp
+     wells/kernels/CompositionalMultiphaseWellConstraintKernels.hpp
+     wells/kernels/SinglePhaseWellConstraintKernels.hpp
      proppantTransport/ProppantTransport.hpp
      proppantTransport/ProppantTransportFields.hpp
      proppantTransport/ProppantTransportKernels.hpp )
@@ -150,6 +160,14 @@ set( fluidFlowSolvers_sources
      wells/SinglePhaseWell.cpp
      wells/kernels/SinglePhaseWellKernels.cpp
      wells/WellControls.cpp
+     wells/WellConstraintsBase.cpp
+     wells/WellInjectionConstraint.cpp
+     wells/WellProductionConstraint.cpp
+     wells/WellBHPConstraints.cpp
+     wells/WellVolumeRateConstraint.cpp
+     wells/WellMassRateConstraint.cpp
+     wells/WellPhaseVolumeRateConstraint.cpp
+     wells/WellLiquidRateConstraint.cpp
      wells/WellSolverBase.cpp
      proppantTransport/ProppantTransport.cpp
      proppantTransport/ProppantTransportKernels.cpp )

src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.hpp

@@ -25,6 +25,8 @@
 #include "physicsSolvers/fluidFlow/wells/WellSolverBase.hpp"
 #include "physicsSolvers/fluidFlow/CompositionalMultiphaseBase.hpp"
 
+#include "physicsSolvers/fluidFlow/wells/WellConstraintsBase.hpp"
+
 namespace geos
 {
 
@@ -92,24 +94,51 @@ class CompositionalMultiphaseWell : public WellSolverBase
   /**@{*/
 
 
+  virtual real64
+  calculateWellResidualNorm( real64 const & time_n,
+                             real64 const & dt,
+                             WellElementSubRegion const & subRegion,
+                             DofManager const & dofManager,
+                             arrayView1d< real64 const > const & localRhs ) override;
+
   virtual real64
   calculateResidualNorm( real64 const & time_n,
                          real64 const & dt,
                          DomainPartition const & domain,
                          DofManager const & dofManager,
                          arrayView1d< real64 const > const & localRhs ) override;
 
+  virtual real64
+  scalingForWellSystemSolution( ElementSubRegionBase & subRegion,
+                                DofManager const & dofManager,
+                                arrayView1d< real64 const > const & localSolution ) override;
+
   virtual real64
   scalingForSystemSolution( DomainPartition & domain,
                             DofManager const & dofManager,
                             arrayView1d< real64 const > const & localSolution ) override;
 
+  virtual bool
+  checkWellSystemSolution( ElementSubRegionBase & subRegion,
+                           DofManager const & dofManager,
+                           arrayView1d< real64 const > const & localSolution,
+                           real64 const scalingFactor ) override;
+
   virtual bool
   checkSystemSolution( DomainPartition & domain,
                        DofManager const & dofManager,
                        arrayView1d< real64 const > const & localSolution,
                        real64 const scalingFactor ) override;
 
+  virtual void
+  applyWellSystemSolution( DofManager const & dofManager,
+                           arrayView1d< real64 const > const & localSolution,
+                           real64 const scalingFactor,
+                           real64 const dt,
+                           DomainPartition & domain,
+                           MeshLevel & mesh,
+                           WellElementSubRegion & subRegion ) override;
+
   virtual void
   applySystemSolution( DofManager const & dofManager,
                        arrayView1d< real64 const > const & localSolution,
@@ -159,6 +188,21 @@ class CompositionalMultiphaseWell : public WellSolverBase
    */
   void updateFluidModel( WellElementSubRegion & subRegion );
 
+  /**
+   * @brief Update well separator using current values of pressure and composition at the reference element
+   * @param subRegion the well subregion containing all the primary and dependent fields
+   * @param targetIndex the targetIndex of the subRegion
+   */
+  void updateSeparator( WellElementSubRegion & subRegion );
+
+  /**
+   * @brief  Calculate well rates at reference element
+   * @param subRegion the well subregion containing all the primary and dependent fields
+   * @param targetIndex the targetIndex of the subRegion
+   */
+
+  void calculateReferenceElementRates( WellElementSubRegion & subRegion );
+
   /**
    * @brief Recompute phase volume fractions (saturations) from constitutive and primary variables
    * @param subRegion the well subregion containing all the primary and dependent fields
@@ -172,6 +216,16 @@ class CompositionalMultiphaseWell : public WellSolverBase
    */
   void updateTotalMassDensity( WellElementSubRegion & subRegion ) const;
 
+  /**
+   * @brief Recompute the perforation rates for all the wells
+   * @param domain the domain containing the mesh and fields
+   */
+  virtual void computeWellPerforationRates( real64 const & time_n,
+                                            real64 const & GEOS_UNUSED_PARAM( dt ),
+                                            ElementRegionManager const & elemManager,
+                                            WellElementSubRegion & subRegion )override;
+
+
   /**
    * @brief Recompute the perforation rates for all the wells
    * @param domain the domain containing the mesh and fields
@@ -183,6 +237,7 @@ class CompositionalMultiphaseWell : public WellSolverBase
    * @brief Recompute all dependent quantities from primary variables (including constitutive models)
    * @param subRegion the well subregion containing all the primary and dependent fields
    */
+  virtual real64 updateWellState( WellElementSubRegion & subRegion ) override;
   virtual void updateState( DomainPartition & domain ) override;
 
   virtual real64 updateSubRegionState( WellElementSubRegion & subRegion ) override;
@@ -197,6 +252,13 @@ class CompositionalMultiphaseWell : public WellSolverBase
 
   integer useTotalMassEquation() const { return m_useTotalMassEquation; }
 
+  virtual void assembleWellFluxTerms( real64 const & time,
+                                      real64 const & dt,
+                                      WellElementSubRegion const & subRegion,
+                                      DofManager const & dofManager,
+                                      CRSMatrixView< real64, globalIndex const > const & localMatrix,
+                                      arrayView1d< real64 > const & localRhs ) override;
+
   /**
    * @brief assembles the flux terms for all connections between well elements
    * @param time_n previous time value
@@ -213,6 +275,14 @@ class CompositionalMultiphaseWell : public WellSolverBase
                                   DofManager const & dofManager,
                                   CRSMatrixView< real64, globalIndex const > const & localMatrix,
                                   arrayView1d< real64 > const & localRhs )override;
+
+  virtual void assembleWellAccumulationTerms( real64 const & time,
+                                              real64 const & dt,
+                                              WellElementSubRegion & subRegion,
+                                              DofManager const & dofManager,
+                                              CRSMatrixView< real64, globalIndex const > const & localMatrix,
+                                              arrayView1d< real64 > const & localRhs ) override;
+
   /**
    * @brief assembles the accumulation term for all the well elements
    * @param domain the physical domain object
@@ -227,6 +297,21 @@ class CompositionalMultiphaseWell : public WellSolverBase
                                           CRSMatrixView< real64, globalIndex const > const & localMatrix,
                                           arrayView1d< real64 > const & localRhs ) override;
 
+  virtual void assembleWellConstraintTerms( real64 const & time_n,
+                                            real64 const & dt,
+                                            WellElementSubRegion const & subRegion,
+                                            DofManager const & dofManager,
+                                            CRSMatrixView< real64, globalIndex const > const & localMatrix,
+                                            arrayView1d< real64 > const & localRhs ) override;
+
+
+  virtual void assembleWellPressureRelations( real64 const & time_n,
+                                              real64 const & dt,
+                                              WellElementSubRegion const & subRegion,
+                                              DofManager const & dofManager,
+                                              CRSMatrixView< real64, globalIndex const > const & localMatrix,
+                                              arrayView1d< real64 > const & localRhs ) override;
+
   /**
    * @brief assembles the pressure relations at all connections between well elements except at the well head
    * @param time_n time at the beginning of the time step
@@ -249,6 +334,8 @@ class CompositionalMultiphaseWell : public WellSolverBase
    */
   void chopNegativeDensities( DomainPartition & domain );
 
+  void chopNegativeDensities( WellElementSubRegion & subRegion );
+
   struct viewKeyStruct : WellSolverBase::viewKeyStruct
   {
     static constexpr char const * dofFieldString() { return "compositionalWellVars"; }
@@ -312,6 +399,7 @@ class CompositionalMultiphaseWell : public WellSolverBase
 
   virtual void initializePostInitialConditionsPreSubGroups() override;
 
+  void saveState( WellElementSubRegion & subRegion );
   virtual void postRestartInitialization() override final;
   /*
    * @brief Utility function that checks the consistency of the constitutive models
@@ -330,11 +418,6 @@ class CompositionalMultiphaseWell : public WellSolverBase
   void validateWellControlsForFluid( WellControls const & wellControls,
                                      constitutive::MultiFluidBase const & fluid ) const;
 
-  /**
-   * @brief Checks injection streams for validity (compositions sum to one)
-   * @param subRegion the well subRegion
-   */
-  void validateInjectionStreams( WellElementSubRegion const & subRegion ) const;
 
   /**
    * @brief Make sure that the well constraints are compatible
@@ -357,6 +440,9 @@ class CompositionalMultiphaseWell : public WellSolverBase
                    real64 const & dt,
                    DomainPartition & domain ) override;
 
+  virtual bool evaluateConstraints( real64 const & time_n,
+                                    WellElementSubRegion & subRegion ) override;
+
 private:
 
   /**
@@ -365,8 +451,23 @@ class CompositionalMultiphaseWell : public WellSolverBase
    */
   void initializeWells( DomainPartition & domain, real64 const & time_n ) override;
 
+  void initializeWell( DomainPartition & domain, MeshLevel & mesh, WellElementSubRegion & subRegion, real64 const & time_n ) override;
+
   virtual void setConstitutiveNames( ElementSubRegionBase & subRegion ) const override;
 
+  template< typename ... GROUPTYPES >
+  void selectLimitingConstraint( real64 const & time_n, integer const coupledIterationNumber, WellElementSubRegion & subRegion );
+
+  void solveConstraint( std::shared_ptr< WellConstraintBase > constraint,
+                        real64 const & time_n,
+                        real64 const & dt,
+                        integer const cycleNumber,
+                        integer const coupledIterationNumber,
+                        DomainPartition & domain,
+                        MeshLevel & mesh,
+                        ElementRegionManager & elemManager,
+                        WellElementSubRegion & subRegion,
+                        DofManager const & dofManager );
 
 
   /// flag indicating whether mass or molar formulation should be used