RieCG

This page describes the numerical method in RieCG a high level, for more details, see Papers.

The numerical method in RieCG is an edge-based finite element method for unstructured meshes, composed of tetrahedral cells, intended for the mathematical modelling of energetic, high-speed, compressible flows. This is an implementation of the method, originally developed by Luo, Baum, and Lohner, Edge-based finite element scheme for the Euler equations, AIAA Journal,1994 to use an asynchronous-by-default, distributed-memory, task-parallel programming paradigm using the Charm++ runtime system. The super-edge-based implementation of the finite element operators help reduce indirect addressing and increases performance.

Note that the RieCG algorithm can be viewed as either an edge-based finite element method on tetrahedron meshes or a node-centered finite volume method (on the dual mesh, consisting of arbitrary polyhedra).

The equations of compressible flow

The equations solved are the 3D unsteady Euler equations, governing inviscid compressible flow,

where the summation convention on repeated indices has been applied, is the density, is the velocity vector, is the specific total energy, is the specific internal energy, and is the pressure. , , and are source terms that arise from the application of the method of manufactured solutions, used for verification; these source terms are zero when computing practical engineering problems. The system is closed with the ideal gas law equation of state

where is the ratio of specific heats. Although the ideal gas law is used most frequently, any arbitrary (i.e. analytic or tabulated) equation of state can be used to relate density and internal energy to pressure.

The numerical method

The above system is solved using an edge-based finite-element method for linear tetrahedra. Using this particular approach, the unknown solution quantities are the conserved variables and are located at the nodes of the mesh. The solution at any point within the computational domain is represented as

where is the tetrahedron containing the point , is the solution at vertex , and is a linear basis function associated with vertex and element . Applying a Galerkin, lumped-mass approximation to the Euler equations gives the following semi-discrete form of the governing equations for a vertex :

where represent all edges connected to , is the volume surrounding , is the numerical flux between and . The first term on the right-hand-side above is evaluated for all vertices within the domain. The last two terms are only evaluated for vertices and surrounding edges that lie on the domain boundary. denotes the boundary flux. The term is effectively the area-weighted normal to the dual face separating and . It is calculated as

where represents all elements attached to edge . This numerical scheme is guaranteed to be conservative since is antisymmetric. The boundary terms and are defined as

For the derivation of the above and other details, see Papers.

Numerical flux

The numerical flux is evaluated at the midpoint between and . We use Rusanov's flux which approximates the flux between and as

where is the maximum wavespeed, defined as

where is the speed of sound at vertex .

Solution reconstruction

The inputs to the flux function above are approximated using a piecewise limited solution reconstruction of the primitive variables. The reconstruction is performed component-wise for each edge as follows:

with

This scheme corresponds to a piecewise linear reconstruction for and a piecewise parabolic reconstruction for . The function can be any total variation diminishing limiter function, e.g., that of van Leer. The gradients of the primitive variables are computed in all mesh points based on the edge-based formulation as

where are scalar components of the primitive variables . The above summations are over all edges that contribute to point v.

Time integration

The solution is advanced using a multi-stage explicit scheme of the form:

where is the right hand side of the discretized governing equations, is the current stage, is the total number of stages, is the time step, and

The explicit Euler time marching scheme is obtained for and the classical 2nd-order Runge-Kutta method is obtained with . The time step size is adaptive and obtained by finding the minimum value for all mesh points at every time step as

where is the Courant number, and is the nodal volume of the dual-mesh associated to point .