Skip to content

Commit

Permalink
Merge pull request #28892 from tanoret/lagging_rc
Browse files Browse the repository at this point in the history
Updating Limiters for Finite Volume Formulation
  • Loading branch information
GiudGiud authored Dec 2, 2024
2 parents 507f39e + fe27fcd commit 4a6a925
Show file tree
Hide file tree
Showing 89 changed files with 1,766 additions and 271 deletions.
19 changes: 19 additions & 0 deletions framework/doc/content/bib/moose.bib
Original file line number Diff line number Diff line change
Expand Up @@ -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}
}
106 changes: 95 additions & 11 deletions framework/doc/content/syntax/Limiters/index.md
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand All @@ -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
179 changes: 91 additions & 88 deletions framework/include/base/MooseFunctorArguments.h
Original file line number Diff line number Diff line change
Expand Up @@ -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
*/
Expand Down Expand Up @@ -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
*/
Expand Down Expand Up @@ -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};
}
}
3 changes: 2 additions & 1 deletion framework/include/fvbcs/FVBoundaryCondition.h
Original file line number Diff line number Diff line change
Expand Up @@ -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;

Expand Down
3 changes: 2 additions & 1 deletion framework/include/fviks/FVInterfaceKernel.h
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
3 changes: 2 additions & 1 deletion framework/include/fvkernels/FVFluxKernel.h
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
3 changes: 2 additions & 1 deletion framework/include/interfaces/FaceArgInterface.h
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
Loading

0 comments on commit 4a6a925

Please sign in to comment.