LohCG

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

LohCG solves the constant-density (incompressible) Navier-Stokes equation modeling inviscid or viscous flows at low Mach numbers. The numerical method employs continuous Galerkin finite elements with linear shapefunctions and stores the unknowns at nodes of a tetrahedron-only mesh. The finite element operators are formulated over edges [1] and implemented using tetrahedron and triangle super-edges [2]. The algorithm is implemented using an asynchronous-by-default, distributed-memory, task-parallel programming paradigm on top of the Charm++ runtime system.

The equations of constant-density flow

The equations solved are the 3D unsteady Navier-Stokes equation, governing Newtonian incompressible flows,

\[\begin{split} \frac{\partial\mbox{\boldmath$v$}}{\partial t} + \mbox{\boldmath$v$}\cdot\nabla\mbox{\boldmath$v$} + \nabla p & = \nabla\mu\nabla\mbox{\boldmath$v$} \\ \nabla\cdot\mbox{\boldmath$v$} & = 0 \end{split}\]

where $\mbox{\boldmath$v$}$ is the velocity vector, $p$ is the pressure, and both the pressure and the viscosity $\mu$ have been normalized by the (constant) density.

The numerical method

The method casts the above system of governing equations into a hyperbolic/parabolic system as

\[\begin{split} \frac{\partial\mbox{\boldmath$v$}}{\partial t} + \mbox{\boldmath$v$}\cdot\nabla\mbox{\boldmath$v$} + \nabla p & = \nabla\mu\nabla\mbox{\boldmath$v$} \\ \frac{1}{c^2}\frac{\partial p}{\partial t} + \nabla\cdot\mbox{\boldmath$v$} &=0 \end{split}\]

where $c$ is the (constant) speed of sound.

Discretization of momentum transport

The above system is discretized in space using an edge-based finite-element method with linear tetrahedra. Using this approach, the velocity vector and pressure are located at the nodes of the mesh. The numerical solution at any point within the computational domain is represented as

\[ v_l(\mbox{\boldmath$x$}) = \sum_{i \in \Omega_h} N^i(\mbox{\boldmath$x$})v^i_l \quad\mbox{and}\quad p(\mbox{\boldmath$x$}) = \sum_{i \in \Omega_h} N^i(\mbox{\boldmath$x$}) p^i \]

where $\Omega_h$ is a tetrahedron containing the spatial point $\mbox{\boldmath$x$}$ , superscript $i$ denotes the discrete numerical solution at vertex $i$ , and $N^i(\mbox{\boldmath$x$})$ is a linear basis function associated with vertex $i$ and element $\Omega_h$ . Applying a Galerkin lumped-mass approximation to the momentum equation yields the following semi-discrete form for a vertex $i$ :

\[ \frac{\text{d}v^i_l}{\text{d}t} = - \frac{1}{V^i} \left[ \sum_{ ij\in i} D^{ij}_k F^{ij}_{kl} + \sum_{ij\in i} B^{ij}_k \left( F^i_{kl} + F^j_{kl}\right) + B^i_k F^i_{kl} \right] - \frac{1}{V^i} \sum_{ ij\in i} K^{ij} \left( v^i_l - v^j_l \right) \]

where $ij\in i$ represent all edges $ij$ connected to vertex $i$ and $V^i$ is the volume surrounding $i$ . $F^{ij}_{kl}$ denotes the numerical flux between vertices $i$ and $j$ while $F^i_{kl}$ the boundary flux. Note that summation is implied on the repeated indices $k$ and $l$ representing spatial directions. The terms $D^{ij}_k$ , $B^{ij}_k$ , $B^i_k$ , and $K^{ij}$ are defined as

\[ \begin{split} D^{ij}_k &= \frac{1}{2} \sum_{\Omega_h \in ij} \int_{\Omega_h} \left( N^i \frac{\partial N^j}{\partial x_k} - N^j \frac{\partial N^i}{\partial x_k} \right) \, \text{d} \Omega \\ B^{ij}_k &= \frac{1}{2} \sum_{\Gamma_h \in ij} \int_{\Gamma_h} N^i N^j n_k \, \text{d} \Gamma \\ B^i_k &= \sum_{\Gamma_h \in i} \int_{\Gamma_h} N^i N^i n_k \, \text{d} \Gamma \\ K^{ij} &= \sum_{\Omega_h \in ij} \int_{\Omega_h} \mu \frac{\partial N^i}{\partial x_k} \frac{\partial N^j}{\partial x_k} \mathrm{d}\Omega \end{split} \]

where $\Omega_h \in ij$ represents all domain elements attached to edge $ij$ , while $\Gamma_h \in ij$ all boundary elements attached to edge $ij$ , and $\Gamma_h \in i$ all boundary elements surrounding boundary point $i$ . For the derivation of the above and other details see [3] and [4]. The numerical flux $F^{ij}_{kl}$ above is evaluated on each edge connecting mesh points $i$ and $j$ . In the LohCG solver this is given by a Galerkin approximation stabilized by artificial viscosity:

\[ D^{ij}_k F^{ij}_{kl} = D^{ij}_k \left( F^i_{kl} + F^j_{kl} \right) + c_d \lambda^{ij} \left| D^{ij}_k \right| \left( v^j_l - v^i_l \right) \]

where the flux at mesh point $i$ is

\[ F^i_{kl} = v^i_k v^i_l + p^i\delta_{kl} \]

while $c_d$ is the artificial viscosity (damping) coefficient, and $\lambda^{ij}=\max(\lambda^i, \lambda^j)$ is the maximum wavespeed across edge $ij$ with

\[ \lambda^i = \left| v^i_k \frac{D^{ij}_k}{\left|D^{ij}_k\right|} \right| + c \]

Instead of the above second-order damping term, another, less diffusive, option in LohCG is to stabilize the advection term via solution reconstruction and limited extrapolation, similar to that of implemented in RieCG, applied to the velocity components.

Discretization of pressure transport

The evolution equation governing the pressure is discretized in a similar fashion to the above using the Galerkin finite element operator over edges:

\[ \frac{\text{d}p^i}{\text{d}t} = - \frac{c^2}{V^i} \left[ \sum_{ ij\in i} D^{ij}_k F^{ij}_k + \sum_{ij\in i} B^{ij}_k \left( F^i_k + F^j_k\right) + B^i_k F^i_k \right] \]

with the flux specifying the velocity divergence. The numerical flux $F^{ij}_k$ is evaluated on each edge connecting mesh points $i$ and $j$ . Similar to the advection term the Galerkin approximation is stabilized by artificial viscosity:

\[ D^{ij}_k F^{ij}_k = D^{ij}_k \left( F^i_k + F^j_k \right) + c_d \lambda^{ij} \left| D^{ij}_k \right| \left( p^j - p^i \right) \]

where the flux at mesh point $i$ is

\[ F^i_k = v^i_k \]

Startup procedure

The boundary and initial conditions for the velocity, prescribed by the user for a given problem, are not assumed to yield a divergence-free vector field and the initial pressure with its boundary conditions is not assumed to be consistent with the initial velocity, i.e., the velocity and pressure are allowed to be specified independently. To ensure consistency of the initial and boundary conditions with the governing equations, the user-given velocity field, $\mbox{\boldmath$v$}^*$ , is decomposed into div-free and curl-free parts as $ \mbox{\boldmath$v$}^* = \mbox{\boldmath$v$} + \nabla\phi $ which, taking the divergence, yields the Poisson equation

\[ \nabla^2 \phi = \nabla \cdot \mbox{\boldmath$v$}^* \]

This is solved for $\phi$ which is then used to project the velocity field onto a divergence-free subspace as

\[ \mbox{\boldmath$v$} = \mbox{\boldmath$v$}^* - \nabla \phi \]

This divergence-free velocity field is then used to compute the initial pressure, satisfying the initial and boundary conditions, by solving the Poisson equation

\[ \nabla^2p = -\nabla\cdot\mbox{\boldmath$v$}\nabla\mbox{\boldmath$v$} \]

Time stepping starts after the above startup procedure.

Time integration

The solution is advanced using an explicit $m$ -stage Runge-Kutta scheme:

\[ \mbox{\boldmath$u$}^k = \mbox{\boldmath$u$}^n + \alpha^k \Delta t \mbox{\boldmath$r$}^{k-1} \]

where $\mbox{\boldmath$u$}=(\mbox{\boldmath$v$},p)$ , $r$ is the right hand side, $k$ is the current stage, $\alpha^k = 1/(1+m-k)$ , and the time step size is obtained by finding the minimum value over all mesh points as

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

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

References

[1] H. Luo, J. Baum, and R. Lohner, Edge-based finite element scheme for the Euler equations, AIAA Journal, 32(6), 1183-1190, 1994.
[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.