LaxCG

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

The numerical method in LaxCG is an edge-based finite element method for unstructured meshes, composed of tetrahedral cells, intended for the mathematical modelling of flows at all Mach numbers. The implementation borrows elements of RieCG, which is not efficient at low Mach numbers, and thus LaxCG applies a preconditioning technique to enable computation of flows at all Mach numbers. The preconditioner is detailed in [1]. The implementation uses the 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 [2].

The equations of compressible flow

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

\[ \begin{split} \frac{\partial U}{\partial t} + \frac{\partial F_j}{\partial x_j} = 0,\qquad U = \begin{Bmatrix} \rho \\ \rho u_i \\ \rho E \end{Bmatrix}, \enskip F_j = \begin{Bmatrix} \rho u_j \\ \rho u_iu_j + p\delta_{ij} \\ u_j(\rho E + p) \end{Bmatrix} \end{split} \]

where the summation convention on repeated indices has been applied, $\rho$ is the density, $u_i$ is the velocity vector, $p$ is the pressure. and the specific total energy, $E$ , is related to the internal energy, $e$ , as

\[ E = e + u_iu_i/2 \]

The system is closed with the ideal gas law equation of state

\[ p = \rho e (\gamma - 1) \]

where $\gamma$ 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.

Time-derivative preconditioning for all Mach numbers

Because of the large ratio of acoustic and convective timescales density-based solvers, such as RieCG, are not efficient, frequently fail to converge, or may converge to the wrong solution in low-Mach-number or constant-density flows at low-speed flow regimes. One way to remedy this stiffness and associated convergence problems is to apply a preconditioner to the governing equations. The idea is to write the system as

\[ \frac{\partial U}{\partial t} = \frac{\partial U}{\partial Q} \frac{\partial Q}{\partial t} = \Gamma \frac{\partial Q}{\partial t} = -\frac{\partial F_j}{\partial x_j} \]

which is solved for the Q vector of primitive variables of pressure, velocity, and temperature, T,

\[ Q=\begin{Bmatrix} p \\ u_1 \\ u_2 \\ u_3 \\ T \end{Bmatrix} \]

The preconditioner is given by

\[ \Gamma = \begin{pmatrix} \Theta & 0 & 0 & 0 & \rho_T \\ \Theta u_1 & \rho & 0 & 0 & \rho_T u_1 \\ \Theta u_2 & 0 & \rho & 0 & \rho_T u_2 \\ \Theta u_3 & 0 & 0 & \rho & \rho_T u_3 \\ \Theta H - 1 & \rho u_1 & \rho u_2 & \rho u_3 & \rho_T H + \rho C_p \\ \end{pmatrix} \]

where the total enthalpy, H, is related to the total energy as

\[ H = E + p/\rho = h + u_iu_i/2, \qquad h = C_p T \]

and

\[ \Theta = \frac{1}{V_r^2} - \frac{\rho_T}{\rho C_p}, \qquad \rho_T = \left.\frac{\partial \rho}{\partial T}\right|_p, \qquad \rho_p = \left.\frac{\partial \rho}{\partial p}\right|_T \]

where the reference velocity is specified as

\[ V_r = \mathrm{min}\big[ c, \mathrm{max}\left( \sqrt{u_iu_i}, K |v_\infty| \right) \big] \]

where c is the local speed of sound, $v_\infty$ is a fixed reference velocity, and K is a constant.

The numerical method

In LaxCG the preconditioned system

\[ \frac{\partial Q}{\partial t} = -\Gamma^{-1}\frac{\partial F_j}{\partial x_j} \]

is discretized with an edge-based finite-element method for linear tetrahedra. The unknown solution quantities are the primitive variables Q and are located at the nodes of the mesh. The solution at any point within the computational domain is represented as

\[ \begin{split} Q(\vec{x}) = \sum_{v \in \Omega_h} N^v(\vec{x}) Q^v \end{split} \]

where $\Omega_h$ is the tetrahedron containing the point $\vec{x}$ , $Q^v$ is the solution of primitive variables at vertex $v$ , and $N^v(\vec{x})$ is a linear basis function associated with vertex $v$ and element $\Omega_h$ . Applying a Galerkin, lumped-mass approximation to the preconditioned system gives the following semi-discrete form of the governing equations for a vertex v:

\[ \begin{split} \frac{\text{d}Q^v}{\text{d}t} = - \frac{1}{V^v} \left(\Gamma^v\right)^{-1} \sum_j \left[ \sum_{ vw\in v} D^{vw}_j F^{vw}_j + \sum_{vw\in v} B^{vw}_j \left( F^v_j + F^w_j\right) + B^v_j F^v_j \right] \end{split} \]

where vw represents all edges connected to v, $V^v$ is the volume surrounding v, $F^{vw}_j$ is the numerical flux between v and w. The first term on the right-hand-side above is evaluated for all vertices within the domain. The last two terms are evaluated only at the domain boundary. $F^v_j$ denotes the boundary flux. The term $D^{vw}_j$ is effectively the area-weighted normal to the dual face separating v and w and calculated as

\[ \begin{split} D^{vw}_j = \frac{1}{2} \sum_{\Omega_h \in vw} \int_{\Omega_h} \left( N^v \frac{\partial N^w}{\partial x_j} - N^w \frac{\partial N^v}{\partial x_j} \right) \, \text{d} \Omega \end{split} \]

where $\Omega_h \in vw$ represents all elements attached to edge vw. This numerical scheme is guaranteed to be conservative since $D^{vw}_j$ is antisymmetric. The boundary terms $B^{vw}_j$ and $B^v_j$ are defined as

\[ \begin{split} B^{vw}_j &= \frac{1}{2} \sum_{\Gamma_h \in vw} \int_{\Gamma_h} N^v N^w n_j \, \text{d} \Gamma \\ B^v_j &= \sum_{\Gamma_h \in v} \int_{\Gamma_h} N^v N^v n_j \, \text{d} \Gamma \end{split} \]

where n is the outward-pointing surface normal.

Note that the above 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). For the derivation of the above and other details, see [1] and [2].

Numerical flux

The numerical flux $F^{vw}_j$ is evaluated at the midpoint between v and w. One option is to use Rusanov's flux which approximates the flux between v and w as

\[ \begin{split} D^{vw}_j F^{vw}_j = D^{vw}_j \left( F^v_j + F^w_j \right) + \lambda \left| D^{vw}_j \right| \left( U^w - U^v \right) \end{split} \]

where the wave speed of the preconditioned system for edge vw is computed as

\[ \lambda = \mathrm{max}\left( S^v, S^w \right) \]

from the wave speeds at node v (and analogously at node w) as

\[ \begin{split} S^v &= |u'| + c' \\ u' &= u_n(1-\alpha) \\ u_n &= u_kn_k \\ c' &= \sqrt{\alpha^2u_n^2 + V_r^2} \\ \alpha &= (1-\beta V_r^2)/2 \\ \beta &= \rho_p + \rho_T/\rho/C_p \end{split} \]

with all quantities evaluated in vertex v.

Another option is to use the Harten-Lax-van Leer-Contact (HLLC) flux using the modified signal velocities, described in [1], which in LaxCG is implemented as given in [3].

Solution reconstruction

The inputs to the flux function above are approximated using a piecewise limited solution reconstruction of the variables, q, denoting the scalar components of the primitive variables, Q. The reconstruction is performed component-wise for each edge $vw$ as follows:

\[ \begin{split} \hat{q}^v &= q^v + \frac{1}{4} \left[ (1-k) \phi(r^v) \delta_1 + (1+k)\phi\left(\frac{1}{r^v}\right)\delta_2 \right]\\ \hat{q}^w &= q^w - \frac{1}{4} \left[ (1-k) \phi(r^w) \delta_3 + (1+k)\phi\left(\frac{1}{r^w}\right)\delta_2 \right] \end{split} \]

with

\[ \begin{split} \delta_1 &= 2x_i^{vw} \frac{\partial q^v}{\partial x_i} - \delta_2 \\ \delta_2 &= q^w - q^v \\ \delta_3 &= 2x_i^{vw} \frac{\partial q^w}{\partial x_i} - \delta_2 \\ x_i^{vw} &= x_i^w - x_i^v \\ r^v &= \delta_2/\delta_1 \\ r^w &= \delta_2/\delta_3 \end{split} \]

This scheme corresponds to a piecewise linear reconstruction for k = -1 and a piecewise parabolic reconstruction for k = 1/3. The function $\phi$ 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

\[ V^v\left(\frac{\partial q}{\partial x_j}\right)^v = \sum_{vw \in v} D_j^{vw}\left(q^w+q^v\right) + \sum_{vw \in v} B_j^{vw}\left(q^w+q^v\right) + B_j^vq^v \]

where q are scalar components of the primitive variables, Q. The above summations are over all edges $vw \in \Omega_h$ that contribute to point v.

Time integration

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

\[ \begin{split} U^{v,k} = U^{v,0} + \alpha^k \Delta t r^{v,k-1} \end{split} \]

where r is the right hand side of the discretized governing equations, k is the current stage, m is the number of stages, $\Delta t$ is the time step size, and

\[ \begin{split} \alpha^k = \frac{1} {1+m-k} \end{split} \]

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

\[ \begin{split} \Delta t = \text{C} \min_{\forall v} \frac{(V^v)^{1/3}}{\lambda^{v}} \end{split} \]

where $\text{C} \leq 1$ is the Courant number.

References

[1] H. Luo, J.D. Baum, R. Lohner, **Extension of Harten-Lax-van Leer Scheme for Flows at All Speeds, AIAA Journal, 43:6, 1160-1166, 2005.
[2] R. Löhner, Edges, stars, superedges and chains, Computer Methods in Applied Mechanics and Engineering, 111, 3-4, 255-263, 1994.
[3] J. Waltz, N.R. Morgan, T.R. Canfield, M.R.J. Charest, L.D. Risinger, J.G. Wohlbier, A three-dimensional finite element arbitrary Lagrangian–Eulerian method for shock hydrodynamics on unstructured grids, Computers & Fluids, 92: 172-187, 2014.
[4] R. Löhner, Applied Computational Fluid Dynamics Techniques: An Introduction Based on Finite Element Methods, Wiley, 2008.