Updating Limiters for Finite Volume Formulation

framework/doc/content/bib/moose.bib

@@ -621,3 +621,22 @@ @article{rao2018stopping
   url = {https://www.sciencedirect.com/science/article/pii/S0021999117306939},
   author = {Kaustubh Rao and Paul Malan and J. Blair Perot}
 }
+
+@inproceedings{venkatakrishnan1993,
+  title={On the accuracy of limiters and convergence to steady state solutions},
+  author={Venkatakrishnan, Venkat},
+  booktitle={31st Aerospace Sciences Meeting},
+  pages={880},
+  year={1993}
+}
+
+@article{harten1997,
+  title={High resolution schemes for hyperbolic conservation laws},
+  author={Harten, Ami},
+  journal={Journal of computational physics},
+  volume={135},
+  number={2},
+  pages={260--278},
+  year={1997},
+  publisher={Elsevier}
+}

framework/doc/content/syntax/Limiters/index.md

@@ -20,7 +20,10 @@ and $u_j^n = u(x_j,t^n)$. A numerical method is TVD if
 TV(u^{n+1}) \leq TV(u^n)
 \end{equation}
 
-## Limiting Process
+Different formulations are used for compressible and incompressible/weakly-compressible
+flow.
+
+## Limiting Process for Compressible Flow
 
 Borrowing notation from [!citep](greenshields2010implementation), we will now
 discuss computation of limited quantities, represented by $\bm{\Psi}_{f\pm}$ where
@@ -66,7 +69,7 @@ non-skewed mesh the definition in [eq:weighting] and
 [!citep](greenshields2010implementation) are the same.
 
 The flux limiter function $\beta(r_{\pm})$ takes different forms as shown in
-[limiter_summary]. $r_{\pm}$ is computed as follows
+[limiter_summary_compressible]. $r_{\pm}$ is computed as follows
 
 \begin{equation}
 r_{\pm} = 2 \frac{\bm{d}_{\pm}\cdot\left(\nabla
@@ -81,14 +84,95 @@ The following limiters are available in MOOSE. We have noted the convergence
 orders of each (when considering that the solution is smooth), whether they are
 TVD, and what the functional form of the flux limiting function $\beta(r)$ is.
 
-!table id=limiter_summary caption=Limiter Overview
-|Limiter class name | Convergence Order | TVD | $\beta(r)$ |
-|------------------ | ----------------- | --- | --- |
-| `VanLeer`          | 2                 | Yes | $\frac{r +\text{abs}(r)}{1 + \text{abs}(r)}$ |
-| `Upwind`           | 1                 | Yes | 0 |
-| `CentralDifference` | 2                | No | 1 |
-| `MinMod`           | 2                 | Yes | $\text{max}(0, \text{min}(1, r))$ |
-| `SOU`              | 2                 | No | $r$ |
-| `QUICK`            | 2                 | No | $\frac{3+r}{4}$ |
+!table id=limiter_summary_compressible caption=Compressible Limiter Overview
+| Limiter class name  | Convergence Order | TVD | $\beta(r)$                                   |
+| ------------------- | ----------------- | --- | -------------------------------------------- |
+| `VanLeer`           | 2                 | Yes | $\frac{r +\text{abs}(r)}{1 + \text{abs}(r)}$ |
+| `Upwind`            | 1                 | Yes | 0                                            |
+| `CentralDifference` | 2                 | No  | 1                                            |
+| `MinMod`            | 2                 | Yes | $\text{max}(0, \text{min}(1, r))$            |
+| `SOU`               | 2                 | No  | $r$                                          |
+| `QUICK`             | 2                 | No  | $\frac{3+r}{4}$                              |
+
+## Limiting Process for Incompressible and Weakly-Compressible flow
+
+A full second-order upwind reconstruction is used for incompressible and weakly-compressible solvers. In this reconstruction, the limited quantity at the face is expressed as follows:
+
+\begin{equation}
+\bm{\Psi}_f = \bm{\Psi}_C + \beta(r) ((\nabla \bm{\Psi})_C \cdot \bm{d}_{fC})
+\end{equation}
+
+where:
+
+- $\bm{\Psi}_f$ is the value of the variable at the face
+- $\bm{\Psi}_C$ is the value of the variable at the cell
+- $(\nabla \bm{\Psi})_C$ is the value of the gradient at the cell, which is computed with second-order methods (Green-Gauss without skewness correction and Least-Squares for skewness corrected)
+- $\bm{d}_{fC}$ is the distance vector from the face to the cell used in the interpolation
+- $\beta(r)$ is the limiting function
+
+Two kinds of limiters are supported: slope-limited and face-value limited. These limiters are defined below.
+
+For slope-limiting, the approximate gradient ratio (or flux limiting ratio) $r$ is defined as follows:
+
+\begin{equation}
+r = 2 \frac{\bm{d}_{NC} \cdot (\nabla \bm{\Psi})_C}{\bm{d}_{NC} \cdot (\nabla \bm{\Psi})_f} - 1
+\end{equation}
+
+where:
+
+- $\bm{d}_{NC}$ is the vector between the neighbor and current cell adjacent to the face
+- $(\nabla \bm{\Psi})_f$ is the gradient of the variable at the face, which is computed by linear interpolation of second-order gradients at the adjacent cells to the face
+
+For face-value limiting, the limiting function is defined as follows:
+
+\begin{equation}
+r =
+\begin{cases}
+    \frac{|\Delta_f|}{\Delta_{\text{max}}} & \text{if } \Delta_f > 0 \\
+    \frac{|\Delta_f|}{\Delta_{\text{min}}} & \text{if } \Delta_f \leq 0
+\end{cases}
+\end{equation}
+
+where:
+
+- $\Delta_f = (\nabla \bm{\Psi})_C \cdot \bm{d}_{fC}$ is the increment at the face
+- $\Delta_{\text{max}} = \bm{\Psi}_{\text{max}} - \bm{\Psi}_C$ is the maximum increment
+- $\Delta_{\text{min}} = \bm{\Psi}_{\text{min}} - \bm{\Psi}_C$ is the minimum increment
+
+The maximum and minimum variable values, $\Delta_{\text{max}}$ and $\Delta_{\text{min}}$, respectively, are computed with a two-cell stencil. In this method, the maximum value is determined as the maximum cell value of the two faces adjacent to the face and their neighbors, respectively. Similarly, the minimum value is computed as the minimum cell value for these cells.
+
+Each of the limiters implemented along with the implementation reference, limiting type, whether they are TVD, and the functional form of the flux limiting function $\beta(r)$ is shown in [limiter_summary_incompressible].
+
+!table id=limiter_summary_incompressible caption=Incompressible/Weakly-Compressible Limiter Overview
+| Limiter class name                              | Limiting Type | TVD | $\beta(r)$                                                        |
+| ----------------------------------------------- | ------------- | --- | ----------------------------------------------------------------- |
+| `VanLeer` [!citep](harten1997)                  | Slope         | Yes | $\frac{r +\text{abs}(r)}{1 + \text{abs}(r)}$                      |
+| `MinMod` [!citep](harten1997)                   | Slope         | Yes | $\text{max}(0, \text{min}(1, r))$                                 |
+| `QUICK` [!citep](harten1997)                    | Slope         | Yes | $\text{min}(1,\text{max}(\text{min}(\frac{1 + 3r}{4}, 2r, 2),0))$ |
+| `SOU` [!citep](harten1997)                      | Face-Value    | No  | $\text{min}(1,1/r)$                                               |
+| `Venkatakrishnan` [!citep](venkatakrishnan1993) | Face-Value    | No  | $\frac{2r+1}{r(2r+1)+1}$                                          |
+
+To illustrate the performance of the limiters, a dispersion analysis is developedand presented in [dispersion].
+This consists of the advection of a passive scalar in a Cartesian mesh at 45 degrees.
+The exact solution, without numerical diffusion, is a straight line at 45 degrees
+dividing the regions with a scalar concentration of 1 and 0.
+
+!alert note
+In general, we recomment using `VanLeer` and `MinMod` limiters for most of the
+applications considering that they provide truly bounded solutions.
+
+!media framework/finite_volume/dispersion.png
+      style=display: block;margin-left:auto;margin-right:auto;width:40%;
+      id=dispersion
+      caption=Dispersion problem, advection in a Cartesian mesh at 45 degrees.
+
+The results and performance of each of the limiters are shown in [dispersion_line].
+This image provides an idea of the limiting action and results that
+can be expected for each of the limiters.
+
+!media framework/finite_volume/dispersion_line.png
+      style=display: block;margin-left:auto;margin-right:auto;width:40%;
+      id=dispersion_line
+      caption=Performance of each of the limiters in a line perpendicular to the advection front.
 
 !bibtex bibliography

framework/include/base/MooseFunctorArguments.h

@@ -75,6 +75,94 @@ struct ElemPointArg
   ElemArg makeElem() const { return {elem, correct_skewness}; }
 };
 
+/**
+ * Argument for requesting functor evaluation at quadrature point locations on an element side.
+ * Data in the argument:
+ * - The element
+ * - The element side on which the quadrature points are located
+ * - The quadrature point index, e.g. if there are \p n quadrature points, we are requesting the\n
+ *   evaluation of the i-th point
+ * - The quadrature rule that can be used to initialize the functor on the given element and side
+ */
+struct ElemSideQpArg
+{
+  /// The element
+  const libMesh::Elem * elem;
+
+  /// The local side index
+  unsigned int side;
+
+  /// The quadrature point index
+  unsigned int qp;
+
+  /// The qudrature rule
+  const QBase * qrule;
+
+  /// The physical location of the quadrature point
+  Point point;
+
+  /**
+   * @returns The conceptual physical location of this data structure
+   */
+  Point getPoint() const { return point; }
+};
+
+/**
+ * State argument for evaluating functors. The iteration type indicates whether you want to evaluate
+ * a functor based on some iterate state of a transient calculation, nonlinear solve, etc. The state
+ * indicates which iterate of the iterate type we want to evaluate on. A state of 0 indicates
+ * "current", e.g. the current time or the current nonlinear iteration (which should actually be
+ * equivalent); a state of 1 indicates the most-recent "old" time or the most recent previous
+ * nonlinear iteration, etc.
+ */
+
+struct StateArg
+{
+  /**
+   * Prevent implicit conversions from boolean to avoid users accidentally constructing a time
+   * argument when they meant to construct a skewness argument, etc.
+   */
+  StateArg(bool) = delete;
+
+  StateArg(unsigned int state_in) : state(state_in), iteration_type(SolutionIterationType::Time) {}
+
+  StateArg(unsigned int state_in, SolutionIterationType iteration_type_in)
+    : state(state_in), iteration_type(iteration_type_in)
+  {
+  }
+
+  /// The state. Zero represents the most recent state, so for any kind of iteration type, a zero
+  /// state represents the current state, e.g. current solution
+  /// One may represent the 'old' value (one before, in the iteration_type specified), and two an 'older' or two steps away state
+  unsigned int state;
+
+  /// The solution iteration type, e.g. time or nonlinear
+  SolutionIterationType iteration_type;
+
+private:
+  StateArg() : state(0), iteration_type(SolutionIterationType::Time) {}
+
+  friend StateArg currentState();
+};
+
+inline StateArg
+currentState()
+{
+  return {};
+}
+
+inline StateArg
+oldState()
+{
+  return {(unsigned int)1};
+}
+
+inline StateArg
+previousNonlinearState()
+{
+  return {(unsigned int)1, SolutionIterationType::Nonlinear};
+}
+
 /**
  * A structure defining a "face" evaluation calling argument for Moose functors
  */
@@ -104,6 +192,9 @@ struct FaceArg
   /// on one side of the face.
   const Elem * face_side;
 
+  /// A member that can be used to define the instance in which the limiters are executed
+  const Moose::StateArg * state_limiter;
+
   /**
    * @returns The conceptual physical location of this data structure
    */
@@ -169,92 +260,4 @@ struct ElemQpArg
    */
   Point getPoint() const { return point; }
 };
-
-/**
- * Argument for requesting functor evaluation at quadrature point locations on an element side.
- * Data in the argument:
- * - The element
- * - The element side on which the quadrature points are located
- * - The quadrature point index, e.g. if there are \p n quadrature points, we are requesting the\n
- *   evaluation of the i-th point
- * - The quadrature rule that can be used to initialize the functor on the given element and side
- */
-struct ElemSideQpArg
-{
-  /// The element
-  const libMesh::Elem * elem;
-
-  /// The local side index
-  unsigned int side;
-
-  /// The quadrature point index
-  unsigned int qp;
-
-  /// The qudrature rule
-  const QBase * qrule;
-
-  /// The physical location of the quadrature point
-  Point point;
-
-  /**
-   * @returns The conceptual physical location of this data structure
-   */
-  Point getPoint() const { return point; }
-};
-
-/**
- * State argument for evaluating functors. The iteration type indicates whether you want to evaluate
- * a functor based on some iterate state of a transient calculation, nonlinear solve, etc. The state
- * indicates which iterate of the iterate type we want to evaluate on. A state of 0 indicates
- * "current", e.g. the current time or the current nonlinear iteration (which should actually be
- * equivalent); a state of 1 indicates the most-recent "old" time or the most recent previous
- * nonlinear iteration, etc.
- */
-
-struct StateArg
-{
-  /**
-   * Prevent implicit conversions from boolean to avoid users accidentally constructing a time
-   * argument when they meant to construct a skewness argument, etc.
-   */
-  StateArg(bool) = delete;
-
-  StateArg(unsigned int state_in) : state(state_in), iteration_type(SolutionIterationType::Time) {}
-
-  StateArg(unsigned int state_in, SolutionIterationType iteration_type_in)
-    : state(state_in), iteration_type(iteration_type_in)
-  {
-  }
-
-  /// The state. Zero represents the most recent state, so for any kind of iteration type, a zero
-  /// state represents the current state, e.g. current solution
-  /// One may represent the 'old' value (one before, in the iteration_type specified), and two an 'older' or two steps away state
-  unsigned int state;
-
-  /// The solution iteration type, e.g. time or nonlinear
-  SolutionIterationType iteration_type;
-
-private:
-  StateArg() : state(0), iteration_type(SolutionIterationType::Time) {}
-
-  friend StateArg currentState();
-};
-
-inline StateArg
-currentState()
-{
-  return {};
-}
-
-inline StateArg
-oldState()
-{
-  return {(unsigned int)1};
-}
-
-inline StateArg
-previousNonlinearState()
-{
-  return {(unsigned int)1, SolutionIterationType::Nonlinear};
-}
 }

framework/include/fvbcs/FVBoundaryCondition.h

@@ -91,7 +91,8 @@ class FVBoundaryCondition : public MooseObject,
   Moose::FaceArg singleSidedFaceArg(
       const FaceInfo * fi = nullptr,
       Moose::FV::LimiterType limiter_type = Moose::FV::LimiterType::CentralDifference,
-      bool correct_skewness = false) const;
+      bool correct_skewness = false,
+      const Moose::StateArg * state_limiter = nullptr) const;
 
   MooseVariableFV<Real> & _var;
 

framework/include/fviks/FVInterfaceKernel.h

@@ -160,7 +160,8 @@ class FVInterfaceKernel : public MooseObject,
       const MooseVariableFV<Real> & variable,
       const FaceInfo * fi = nullptr,
       Moose::FV::LimiterType limiter_type = Moose::FV::LimiterType::CentralDifference,
-      bool correct_skewness = false) const;
+      bool correct_skewness = false,
+      const Moose::StateArg * state_limiter = nullptr) const;
 
   /// To be consistent with FE interfaces we introduce this quadrature point member. However, for FV
   /// calculations there should every only be one qudrature point and it should be located at the

framework/include/fvkernels/FVFluxKernel.h

@@ -118,7 +118,8 @@ class FVFluxKernel : public FVKernel,
   Moose::FaceArg singleSidedFaceArg(
       const FaceInfo * fi = nullptr,
       Moose::FV::LimiterType limiter_type = Moose::FV::LimiterType::CentralDifference,
-      bool correct_skewness = false) const;
+      bool correct_skewness = false,
+      const Moose::StateArg * state_limiter = nullptr) const;
 
   /**
    * Returns whether to avoid execution on a boundary

framework/include/interfaces/FaceArgInterface.h

@@ -44,7 +44,8 @@ class FaceArgProducerInterface : public FaceArgInterface
   Moose::FaceArg makeFace(const FaceInfo & fi,
                           const Moose::FV::LimiterType limiter_type,
                           const bool elem_is_upwind,
-                          const bool correct_skewness = false) const;
+                          const bool correct_skewness = false,
+                          const Moose::StateArg * state_limiter = nullptr) const;
 
   /**
    * Make a functor face argument with a central differencing limiter, e.g. compose a face