ALECG hydrodynamics

Contents

Inciter supports multiple hydrodynamics schemes. This page describes the ALECG method at a high level, for more details, see Inciter papers.

The basic ALECG algorithm 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 developed by Waltz and coworkers, see Inciter papers, to solve the governing equations of compressible flow to use an asynchronous-by-default, distributed-memory, task-parallel programming paradigm using the Charm++ runtime system.

Note that the ALECG 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,

\[ \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}, \enskip S = \begin{Bmatrix} S_\rho \\ S_{u,i} \\ S_E \end{Bmatrix}, \end{split} \]

where the summation convention on repeated indices has been applied, $\rho$ is the density, $u_i$ is 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. $S_\rho$ , $S_{u,i}$ , and $S_E$ 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

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

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.

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 Euler equations gives the following semi-discrete form of the governing equations for a vertex $v$ :

\[ \begin{split} \frac{\text{d}U^v}{\text{d}t} = - \frac{1}{V^v} \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$ represent all edges connected to $v$ , $V$ is the volume surrounding $v$ , $F^{vw}_j$ 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. As such, $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 an edge $vw$ . The described numerical scheme is guaranteed conservative since $D^{vw}_j$ is antisymmetric.

The boundary terms $B^{vw}_j$ and $B^v_j$ are defined by the following:

\[ \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 and other details, see Inciter papers.

Numerical flux

The numerical flux $F^{vw}_j$ is evaluated at the midpoint between $v$ and $w$ . Here we used the Rusanov 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) + \left| D^{vw}_j \right| \lambda^{vw} \left( U^w - U^v \right) \end{split} \]

where $\lambda^{vw}=\max (\lambda_v, \lambda_w)$ is the maximum wavespeed, defined as

\[ \begin{split} \lambda_v = \left| u_j \frac{D^{vw}_j}{\left|D^{vw}_j\right|} \right| + c^v \end{split} \]

where $c^v$ is the speed of sound at vertex $v$ .

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 $vw$ as follows:

\[ \begin{split} \hat{u}^v &= u^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{u}^w &= u^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 u^v}{\partial x_i} - \delta_2 \\ \delta_2 &= u^w - u^v \\ \delta_3 &= 2x_i^{vw} \frac{\partial u^w}{\partial x_i} - \delta_2 \end{split} \]

where $x_i^{vw}=x_i^w - x_i^v$ and

\[ \begin{split} r^v &= \frac{\delta_2}{\delta_1} \\ r^w &= \frac{\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. We use the Van Leer and minmod functions.

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 \cdot \left. \frac{\text{d}U^v}{\text{d}t}\right|^{k-1} \end{split} \]

where $k$ is the current stage, $m$ is the total number of stages, $\Delta t$ is the time step, 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 global time step is obtained by finding the minimum value for all edges in the domain. The time step for an edge $vw$ is calculated as

\[ \begin{split} \Delta t^{vw} = \frac{\text{C}l^{vw}}{\lambda^{vw}} \end{split} \]

where $\text{C} \leq 1$ is the Courant number, and $l^{vw}$ is the legnth of an edge $vw$ . The global time step is thus

\[ \begin{split} \Delta t = \min_{\forall vw} \left( \Delta t^{vw} \right). \end{split} \]

Example setup, run, and solution

For an example test set up see the page on steady-state solution of transconic flow around an F16. For other examples, see Inciter papers.