rDG Module

The MOOSE rDG module is a library for the implementation of simulation tools that solve convection-dominated problems using the class of so-called reconstructed discontinuous Galerkin (rDG) methods. The specific rDG method implemented in this module is rDG(P0P1), which is equivalent to the second-order cell-centered finite volume method (FVM). Cell-centered FVMs are regarded as a subset of rDG methods in the case when the baseline polynomial solution in each element is a constant monomial. The FVMs are the most widely used numerical methods in areas such as computational fluid dynamics (CFD) and heat transfer, computational acoustics, and magnetohydrodynamics (MHD). This module provides a systematic solution for implementing all required components in a second-order FVM such as slope reconstruction, slope limiting, numerical flux, and proper boundary conditions.


To avoid ambiguity, the term "rDG" should be understood as "cell-centered FVM" in this module, and these two terms are used interchangeably in the following text.

The Advection Equation

To best help developers and/or users better understand the rDG method so that they can try to develop new applications or use existing rDG-based codes, we use the advection equation (the simplest convection-dominated system that we can think of) as an example to describe the implementation of cell-centered FVM in MOOSE.

The advection equation for a conserved quantity described by a scalar field is expressed mathematically by a continuity equation as

(1) where is the divergence operator, and is the velocity field. If the flow is assumed to be incompressible, that is, the velocity field satisfies the above equation can be rewritten as

Finite Volume Method

Eq. (1) can be discretized in space using a cell-centered FVM. In an FVM, the computational domain is divided by a set of non-overlapping control volumes , which can be one or a combination of the most common element types, e.g. line segment in 1D, triangles and quadrilaterals in 2D, and tetrahedra, prisms, pyramids, and hexahedra in 3D. On each control volume, the integral form of the governing equations is required to be satisfied, (2) The cell-averaged conservative variable, , is taken to be the unknown and defined by where is the volume of the control volume . The following equation can then be derived using the divergence theorem, (3) where denotes an interior common face between cell and , denotes a face on the boundary of domain ; and and are the unit vectors normal to face and , respectively. For each cell , represents a set of neighboring cells, , sharing a common face, .

Because the numerical solution is discontinuous between cell interfaces, the interface fluxes are not uniquely defined. The flux, , appearing in the second term of Eq. (3) should be replaced by a numerical flux function , i.e., where and are the conservative variable at the "left" and "right" side of the cell interface (). In the case of first-order FVM, the solution in each cell is assumed to be constant in space. Then on any interior face, , the two states are simply and . In order to guarantee consistency and conservation, is required to satisfy


Similarly, the flux function on the domain boundary, , should be determined by , i.e., with the use of appropriate boundary conditions satisfying the characteristic theory.

Finally, the boundary integration in Eq. (3) is approximated using one point quadrature at the midpoint of the face, and the semi-discrete form of the equations may be written as (4) where is the length of cell edge in 2D, and area of cell face in 3D.

By assembling all the elemental contributions, a system of ordinary differential equations governing the evolution of the discrete solution in time can be written as where denotes the mass matrix, is the global vector of the degrees of freedom, and is the residual vector. has a block diagonal structure that couples the degrees of freedom of the unknown vector associated to only within . As a result, the inverse of can be easily computed in advance considering one cell at a time.

Numerical Flux Scheme

In the example of advection equation, the flux function is approximated using the upwind scheme, i.e., where and where and

TVD Slope Limiters

First-order FVMs are in general stable on arbitrary grids. However, second-order FVMs based on piecewise linear reconstruction suffer from non-physical oscillations in the vicinity of strong discontinuities for convection-dominant flows. One common approach to address this issue is an appropriate slope limiter. Slope limiters are widely used in FVMs to modify the piecewise linearly reconstructed gradients of solution variables, and thus to satisfy the total-variational diminishing (TVD) condition.

In MOOSE, implementation of the TVD slope limiters is possible on 1D unstructured grids, because it is trivial to know the indices of the "left" and "right" neighboring elements of the -th element, i.e., and , during a loop over the elements. The three classical slope limiters implemented in the example of advection equation are described below.

Minmod Slope Limiter

One choice of slope that gives second-order accuracy for smooth solutions while still satisfying the TVD property is the minmod slope where the minmod function of two arguments is defined by If and have the same sign, then this selects the one that is smaller in modulus, else it returns zero.

Rather than defining the slope on the -th cell by always using the downwind difference (which would give the Lax-Wendroff method), or by always using the upwind difference (which would give the Beam-Warming method), the minmod method compares the two slopes and chooses the one that is smaller in magnitude. If the two slopes have different sign, then the value must be a local maximum or minimum, and it is easy to check in this case that we must set in order to satisfy the TVD condition. The minmod method does a fairly good job of maintaining good accuracy in the smooth hump and also sharp discontinuities in the square wave, with no oscillations. Sharper resolution of discontinuities can be achieved with other limiters that do not reduce the slope as severely as minmod near a discontinuity.

Superbee Slope Limiter

One choice of limiter that gives the sharper reconstruction, while still giving second order accuracy for smooth solutions, is the so-called superbee limiter introduced by Roe (1985): (5) where and \begin{equation} \phi^{(2)}_i = {\rm minmod} \left( 2\frac{Q_{i+1}-Q_i}{\Delta x}, \frac{Q_i-Q_{i-1}}{\Delta x} \right). \begin{equation} Each one-sided slope is compared with twice the opposite one-sided slope. Then the maxmod function in Eq. (5) selects the argument with larger modulus. In regions where the solution is smooth this will tend to return the larger of the two one-sided slopes, but will still be giving an approximation, and hence we expect second-order accuracy. The superbee limiter is also TVD in general.

With the superbee method, the discontinuity stays considerably sharper than with the minmod method. On the other hand, there is a tendency of the smooth hump to become steeper and squared off. This is sometimes a problem with superbee — by choosing the larger of the neighboring slopes it tends to steepen smooth transitions near inflection points.

MC Slope Limiter

Another popular choice is the monotonized central-difference limiter (MC limiter), which was proposed by Van Leer (1977): (6) This compares the central difference of Fromm method with twice the one-sided slope to either side. In smooth regions this reduces to the centered slope of Fromm method and hence does not tend to artificially steepen smooth slopes to the extent that superbee does. The MC limiter appears to be a good default choice for a wide class of problems.

Example Problem

Problem Description

In a test case validating the rDG implementation of the advection equation, the three slope limiters introduced above are used for simulating a right-going square-shaped wave in 1D. The initial condition (I.C.) at contains contact discontinuities at and . For simplicity, the wave speed is set to .

Figure 1: Time evolution of rDG solution for simulating a right-going square-shaped wave in 1D.g

Numerical Results

To demonstrate the oscillation-free solution quality during wave propagation, an animation is presented below. In addition, the numerical results at are presented in Figure 1.

Input File

The content of some input file blocks is described in detail for clarity.


  1. It is mandatory to declare order = CONSTANT and family = MONOMIAL, which specifies the piecewise cell-average solution variable for the cell-centered finite volume method.

  2. It is convenient to provide some parameters for rDG-related objects used in an input file, such as slope_reconstruction = rslope and slope_limiting = lslope for slope reconstruction and slope limiting.

  3. If an explicit time integration method is used, it is convenient to declare implicit = false here, so that Jacobian matrices will not be computed.

  order = CONSTANT
  family = MONOMIAL
  u = u
  slope_limiting = lslope
  implicit = false


  type = GeneratedMesh
  dim = 1
  xmin = 0
  xmax = 1
  nx = 100


  1. A one-dimensional computational domain ranging between and is defined, with 100 elements equally distributed in the domain.


    type = PiecewiseConstant
    axis = x
    direction = right
    xy_data = '0.1 0.5
               0.6 1.0
               1.0 0.5'


  1. In this case, a piecewise constant function is used to specify the initial condition of the square-shaped wave profile.


    type = AEFVSlopeLimitingOneD
    execute_on = 'linear'
    scheme = 'none' #none | minmod | mc | superbee

    type = AEFVUpwindInternalSideFlux
    execute_on = 'linear'

    type = AEFVFreeOutflowBoundaryFlux
    execute_on = 'linear'


  1. The prefix AEFV represents "Advection Equation Finite Volume" — an abbreviation to uniquely name these classes.

  2. AEFVSlopeReconstructionOneD does not do any work, but has to be in place for code consistency. In 1D the slope reconstruction and limiting can be accomplished in one user object.

  3. AEFVSlopeLimitingOneD calculates the limited slope for each element in 1D.

  4. AEFVUpwindInternalSideFlux calculates the internal side flux using a simple upwind scheme.

  5. AEFVFreeOutflowBoundaryFlux calculates the boundary side flux using a free outflow BC.




  1. Declare the nonlinear variable as u. The type and family of the variable have been declared in GlobalParams for convenience.


    implicit = true
    type = TimeDerivative
    variable = u


  1. Always set implicit = true for time derivative kernels when using explicit time integration.

  2. In FVMs, there is no volumetric integration for flux terms.

  3. Other possible kernels in this block include source terms. In this example, we do not have any.


    type = AEFVKernel
    variable = u
    component = 'concentration'
    flux = internal_side_flux


  1. Internal side flux terms should be declared in this block.


    type = AEFVBC
    boundary = 'left right'
    variable = u
    component = 'concentration'
    flux = free_outflow_bc


  1. Boundary side flux terms should be declared in this block.


    type = AEFVMaterial
    block = 0


  1. This block does not calculate actual material properties. It is used to trigger the calculation of slopes in every element, and then interpolate variable values at side centers.



  1. Philip L Roe. Some contributions to the modelling of discontinuous flows. In Large-Scale Computations in Fluid Mechanics, 163–193. 1985.[BibTeX]
  2. Bram Van Leer. Towards the ultimate conservative difference scheme. IV. A new approach to numerical convection. Journal of computational physics, 23(3):276–299, 1977.[BibTeX]