ZalCG

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

The numerical method in ZalCG is an edge-based implementation of a finite element method for unstructured meshes, composed of tetrahedral cells. The solver is intended for the computing of high-speed compressible flows. A good reference for the finite element method itself is [1] and on flux-corrected transport is [2]. The implementation of the method uses an asynchronous-by-default, distributed-memory, task-parallel programming paradigm using the Charm++ runtime system.

The equations of compressible flow

The governing equations solved are the 3D unsteady Euler equations, governing inviscid compressible flow, written here in the flux-conservative form

\[ \begin{split} \frac{\partial U}{\partial t} + \frac{\partial F_j}{\partial x_j} = S,\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, and $\rho$ is the density, $u_i$ the velocity vector, $E=u_iu_i/2+e$ is the specific total energy, $e$ is the specific internal energy, and $p$ is the pressure. The system is closed with the ideal gas law equation of state

\[ \begin{split} p = \rho e (\gamma - 1), \end{split} \]

where $\gamma$ is the ratio of specific heats. The exact form of the equation of state is of secondary importance, 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 $U$ and are located at the nodes of the mesh. The solution at any point within the computational domain is represented as

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

where $\Omega_h$ is the tetrahedron containing the point $\vec{x}$ , $U^v$ is the solution 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 Euler equations gives the following semi-discrete form of the governing equations for a vertex $v$ :

\[ \begin{split} V^v\frac{\text{d}U^v}{\text{d}t} = - \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 \end{split} \]

where $vw$ represent all edges connected to $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 only evaluated for vertices $v$ and surrounding edges $vw$ that lie on 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$ . It is 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} \]

For the derivation of the above and other details, see [1] and [3].

Numerical flux

The numerical flux $F^{vw}_j$ between mesh points $v$ and $w$ is computed as

\[ \begin{split} F^{vw}_j & = 2 F_j( U_{vw}^{n+1/2} ) \\ U_{vw}^{n+1/2} & = \frac{1}{2}\left( U^v + U^w \right) - \frac{\Delta t}{2} \left.\frac{\partial F_j}{\partial x_j}\right|^{vw} \\ \left.\frac{\partial F_j}{\partial x_j}\right|^{vw} & = \frac{l^{vw}_j}{l_kl_k}( F^v_j - F^w_j ) \end{split} \]

where ${l_j}$ denotes the edge difference vector

\[ l^{vw}_j=x^v_j - x^w_j \]

The above flux prescribes a Lax-Wendroff (or Taylor-Galerkin) scheme for the edge-based formulation. Monotonicity is ensured by flux-corrected transport (FCT).

Flux-corrected transport

Flux-corrected transport (FCT) is a way to circumvent the consequence of Godunov's theorem, which states that no linear scheme of order greater than 1 will yield monotonic solutions. Accordingly, FCT is a nonlinear scheme that combines a high-, and a low-order scheme using limiting.

The high-order scheme in the FCT algorithm in ZalCG is given by the Taylor-Galerkin flux above with the Euler equations discretized in time with the forward Euler method as

\[ V \Delta U^h = \Delta t r(U^n) \]

where $U^n$ is the solution at the previous time step $n$ and the notation of the numerical solution represented at mesh node $v$ is dropped for brevity. The low-order scheme is constructed from the high-order scheme in by adding mass diffusion

\[ V \Delta U^l = \Delta t r(U^n) - k(U^n) \]

with the diffusion operator

\[ \begin{split} k(U) & = \sum_{vw \in v} K^{vw} (U^v-U^w) \\ K^{vw} & = c_\tau \sum_{\Omega_h \in vw} \int_{\Omega_h}N^vN^v - N^vN^w\mathrm{d}\Omega \end{split} \]

where $c_\tau$ denotes a diffusion coefficient. The steps to arrive at the diffusion term for edges are analogous to those of the Laplacian operator, invoking the symmetry property of the diffusion operator for linear tetrahedra as

\[ K^{ik}\!=\!-\sum_{j\ne i}K^{jk} \]

In the FCT procedure $c_\tau\!=\!1$ guarantees monotonic solutions. The above specifications lead to the following antidiffusive edge contributions (AEC)

\[ \mathrm{AEC}^{vw} = \Delta U^h - \Delta U^l = k(U^n) \]

which is then used to update the low order solution after performing the limiting procedure to arrive at the solution at the new time step $n+1$

\[ \Delta U^{n+1} = U^l + \sum_{vw \in v} \lim[\mathrm{AEC}^{vw}] \]

Computing the limit coefficients follows the same procedure, detailed in [2], implemented over the edges, given below.

The limiting procedure

The procedure first aggregates to points i the signed contributions of AECs:

\[ \begin{split} P_i^{\pm} = \sum_e \begin{Bmatrix} \max \\ \min \end{Bmatrix} (0,\mathrm{AEC}_e) \end{split} \]

The limiting procedure requires the maximum and minimum nodal values of the low-order solution and the previous solution,

\[ \begin{split} {\mbox{\boldmath$U$}}_i^* = \begin{Bmatrix} \max \\ \min \end{Bmatrix} ({\mbox{\boldmath$U$}}^l_i,{\mbox{\boldmath$U$}}^n). \end{split} \]

Another alternative is to only consider the low order solution, when computing the allowed solution bounds, which leads to the so-called 'clipping limiter' in place of the equation above. This is followed by computing the maximum and minimum nodal values of all edges connected to a node

\[ \begin{split} {\mbox{\boldmath$U$}}_e^* = \begin{Bmatrix} \max \\ \min \end{Bmatrix} ({\mbox{\boldmath$U$}}^*_A,{\mbox{\boldmath$U$}}^*_B), \end{split} \]

then computing the maximum and minimum unknowns of the edges surrounding each node,

\[ \begin{split} {\mbox{\boldmath$U$}}_i^{\tiny\begin{matrix} \max \\ \min \end{matrix}} = \begin{Bmatrix} \max \\ \min \end{Bmatrix} ({\mbox{\boldmath$U$}}^*_1,{\mbox{\boldmath$U$}}^*_2,\dots,{\mbox{\boldmath$U$}}^*_m). \end{split} \]

The limit coefficients will be computed (see below) based on $P_i^{\pm}$ and the maximum and minimum increments and decrements the nodal solution values are allowed to achieve,

\[ \begin{split} Q_i^{\pm} = {\mbox{\boldmath$U$}}_i^{\tiny\begin{matrix} \max \\ \min \end{matrix}} - {\mbox{\boldmath$U$}}^l. \end{split} \]

Defining the ratios of positive and negative edge contributions for each node $i$ that ensure monotonicity as

\[ \begin{split} R_i^{\pm} = \left\{\begin{matrix} \min(1,Q_i^{\pm}/P_i^{\pm}) & P_i^+ > 0 > P_i^- \\ 1 & \mathrm{otherwise} \end{matrix} \right., \end{split} \]

the limit coefficient for each edge is taken as the most conservative ratio

\[ \begin{split} C_e = \min_{i \in \Omega_e} \left\{\begin{matrix} R_i^+ & \mathrm{AEC}>0 \\ R_i^- & \mathrm{AEC}<0 \end{matrix}\right.. \end{split} \]

The limited AEC is then scatter-added to nodes of the low-order solution:

\[ \begin{split} U_i^{n+1} = U_i^l + \sum_{i \in \Omega_e} C_e \cdot \mathrm{AEC}. \end{split} \]

The above procedure is general, works on the numerical solution (instead on fluxes or slopes), and does not require a Riemann solver. The same procedure can be used for each scalar in a system of equations. This works well for independent scalars, but for coupled system of equations additional techniques have been developed to reflect the coupled nature of the equations in the limiting procedure.

References

[1] R. Löhner, Applied Computational Fluid Dynamics Techniques: An Introduction Based on Finite Element Methods, Wiley, 2008.
[2] Lohner, Morgan, Peraire, Vahdati, Finite element flux-corrected transport (FEM-FCT) for the Euler and Navier–Stokes equations, International Journal for Numerical Methods in Fluids, 7, 1093-1109, 1987.
[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.