refactor: clarify new gravity, clean up kernel flags usage, move CFL computations to FVM solver

src/coreComponents/finiteVolume/docs/FiniteVolume.rst.txt

@@ -16,17 +16,27 @@ The numerical flux is obtained using the following expression for the mass flux
 .. math::
   F_{KL} = \Upsilon_{KL} \frac{\rho^{upw}}{\mu^{upw}} \big( p_K - p_L - \rho^{avg} g ( d_K - d_L ) \big),
 
-where :math:`p_K` is the pressure of cell :math:`K`, :math:`d_K` is the depth of cell :math:`K`, and :math:`\Upsilon_{KL}` is the standard TPFA transmissibility coefficient at the interface.
+
+where :math:`p_K` is the pressure of cell :math:`K`, :math:`\rho^{avg}` is the average fluid density, :math:`d_K` is the depth of cell :math:`K`, and :math:`\Upsilon_{KL}` is the standard TPFA transmissibility coefficient at the interface.
 The fluid density, :math:`\rho^{upw}`, and the fluid viscosity, :math:`\mu^{upw}`, are upwinded using the sign of the potential difference at the interface. 
 
-This is currently the only available discretization in the :ref:`CompositionalMultiphaseFlow`. 
+For :ref:`CompositionalMultiphaseFlow` there are different options how :math:`\rho^{avg}` can be computed, which are controlled by `gravityDensityScheme` parameter:
+
+#. `ArithmeticAverage`: :math:`\rho^{avg}` is computed using simple arithmetic average: :math:`\rho^{avg} = 0.5 \cdot (rho_K + rho_L)`, where :math:`rho_K` and :math:`rho_K` are densities in the two cells.
+
+#. `PhasePresence`: average phase density is computed using checking for phase presence:
+
+    * :math:`\rho^{avg} = 0.5 \cdot (\rho_K + \rho_L)` if phase is present in both cells :math:`K` and :math:`L`
+
+    * :math:`\rho^{avg} = \rho_K` if phase is present only in cell :math:`K`
+
+    * :math:`\rho^{avg} = \rho_L` if phase is present only in cell :math:`L`
 
 Hybrid FVM
 ~~~~~~~~~~
 
 This discretization scheme overcomes the limitations of the standard TPFA on non K-orthogonal meshes.
 The hybrid finite-volume scheme--equivalent to the well-known hybrid Mimetic Finite Difference (MFD) scheme--remains consistent with the pressure equation even when the mesh does not satisfy the K-orthogonality condition.
-This numerical scheme is currently implemented in the `SinglePhaseHybridFVM` solver.
 
 The hybrid FVM scheme uses both cell-centered and face-centered pressure degrees of freedom.
 The one-sided face flux, :math:`F_{K,f}`, at face :math:`f` of cell :math:`K` is computed as:
@@ -60,5 +70,3 @@ For a given interior face :math:`f` between two neighboring cells :math:`K` and
 We obtain a numerical scheme with :math:`n_{\textit{cells}}` cell-centered degrees of freedom and :math:`n_{\textit{faces}}` face-centered pressure degrees of freedom.
 The system involves :math:`n_{\textit{cells}}` mass conservation equations and :math:`n_{\textit{faces}}` face-based constraints.
 The linear systems can be efficiently solved using the MultiGrid Reduction (MGR) preconditioner implemented in the Hypre linear algebra package.  
-
-The implementation of the hybrid FVM scheme for :ref:`CompositionalMultiphaseFlow` is in progress.

src/coreComponents/physicsSolvers/PhysicsSolverBase.cpp

@@ -458,15 +458,6 @@ real64 PhysicsSolverBase::setNextDtBasedOnNewtonIter( real64 const & currentDt )
   return nextDt;
 }
 
-
-real64 PhysicsSolverBase::setNextDtBasedOnCFL( const geos::real64 & currentDt, geos::DomainPartition & domain )
-{
-  GEOS_UNUSED_VAR( currentDt, domain );
-  return LvArray::NumericLimits< real64 >::max;       // i.e., not implemented
-}
-
-
-
 real64 PhysicsSolverBase::linearImplicitStep( real64 const & time_n,
                                               real64 const & dt,
                                               integer const GEOS_UNUSED_PARAM( cycleNumber ),

src/coreComponents/physicsSolvers/PhysicsSolverBase.hpp

@@ -239,17 +239,6 @@ class PhysicsSolverBase : public ExecutableGroup
   virtual real64 setNextDtBasedOnStateChange( real64 const & currentDt,
                                               DomainPartition & domain );
 
-  /**
-   * @brief function to set the next dt based on state change
-   * @param [in]  currentDt the current time step size
-   * @param[in] domain the domain object
-   * @return the prescribed time step size
-   */
-  virtual real64 setNextDtBasedOnCFL( real64 const & currentDt,
-                                      DomainPartition & domain );
-
-
-
   /**
    * @brief Entry function for an explicit time integration step
    * @param time_n time at the beginning of the step

src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.hpp

@@ -68,6 +68,8 @@ class CompositionalMultiphaseBase : public FlowSolverBase
 
   virtual void registerDataOnMesh( Group & meshBodies ) override;
 
+  virtual void registerDataForCFL( Group & meshBodies ) { GEOS_UNUSED_VAR( meshBodies ); }
+
   /**
    * @defgroup Solver Interface Functions
    *
@@ -268,7 +270,6 @@ class CompositionalMultiphaseBase : public FlowSolverBase
     static constexpr char const * allowLocalCompDensChoppingString() { return "allowLocalCompDensityChopping"; }
     static constexpr char const * useTotalMassEquationString() { return "useTotalMassEquation"; }
     static constexpr char const * useSimpleAccumulationString() { return "useSimpleAccumulation"; }
-    static constexpr char const * useNewGravityString() { return "useNewGravity"; }
     static constexpr char const * minCompDensString() { return "minCompDens"; }
     static constexpr char const * maxSequentialCompDensChangeString() { return "maxSequentialCompDensChange"; }
     static constexpr char const * minScalingFactorString() { return "minScalingFactor"; }
@@ -370,19 +371,12 @@ class CompositionalMultiphaseBase : public FlowSolverBase
   virtual real64 setNextDtBasedOnStateChange( real64 const & currentDt,
                                               DomainPartition & domain ) override;
 
-  void computeCFLNumbers( DomainPartition & domain, real64 const & dt, real64 & maxPhaseCFL, real64 & maxCompCFL );
-
-  /**
-   * @brief function to set the next time step size
-   * @param[in] currentDt the current time step size
-   * @param[in] domain the domain object
-   * @return the prescribed time step size
-   */
-  real64 setNextDt( real64 const & currentDt,
-                    DomainPartition & domain ) override;
 
-  virtual real64 setNextDtBasedOnCFL( real64 const & currentDt,
-                                      DomainPartition & domain ) override;
+  virtual void computeCFLNumbers( DomainPartition & domain, real64 const & dt, real64 & maxPhaseCFL, real64 & maxCompCFL )
+  {
+    GEOS_UNUSED_VAR( domain, dt, maxPhaseCFL, maxCompCFL );
+    GEOS_ERROR( GEOS_FMT( "{}: computeCFLNumbers is not implemented for {}", getDataContext(), getCatalogName() ) );
+  }
 
   virtual void initializePostInitialConditionsPreSubGroups() override;
 
@@ -487,9 +481,6 @@ class CompositionalMultiphaseBase : public FlowSolverBase
   /// flag indicating whether simple accumulation form is used
   integer m_useSimpleAccumulation;
 
-  /// flag indicating whether new gravity treatment is used
-  integer m_useNewGravity;
-
   /// minimum allowed global component density
   real64 m_minCompDens;
 
@@ -500,9 +491,6 @@ class CompositionalMultiphaseBase : public FlowSolverBase
   real64 m_sequentialCompDensChange;
   real64 m_maxSequentialCompDensChange;
 
-  /// the targeted CFL for timestep
-  real64 m_targetFlowCFL;
-
 private:
 
   /**