ChoCG

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

ChoCG 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 uses a projection scheme, in which a single time step consists of the following steps

  • advective/diffusive prediction: $\mbox{\boldmath$v$}^n \rightarrow \mbox{\boldmath$v$}^*$

    \[ \left[\frac{1}{\Delta t} - \theta\nabla\mu\nabla\right](\mbox{\boldmath$v$}^*-\mbox{\boldmath$v$}^n) + \mbox{\boldmath$v$}^n\cdot\nabla\mbox{\boldmath$v$}^n + \nabla p^n = \nabla\mu\nabla\mbox{\boldmath$v$}^n \]

  • pressure correction: $p^n \rightarrow p^{n+1}$

    \[ \nabla^2(p^{n+1}-p^n) = \frac{\nabla\cdot\mbox{\boldmath$v$}^*}{\Delta t} \\ \]

  • velocity projection: $\mbox{\boldmath$v$}^* \rightarrow \mbox{\boldmath$v$}^{n+1}$

    \[ \mbox{\boldmath$v$}^{n+1} = \mbox{\boldmath$v$}^* - \Delta t\nabla(p^{n+1}-p^n) \]

Here $\theta$ denotes the implicitness factor for the viscous term: $\theta$ =1: implicit, $\theta$ =0.5: Crank-Nicholson.

Discretization of advection and diffusion

The above advective/diffusive step 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 $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 gives the following semi-discrete form for a vertex $i$ :

\[ \begin{split} \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) \end{split} \]

where $ij\in i$ represent all edges $ij$ connected to vertex $i$ . $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 see [3] and [4]. The numerical flux $F^{ij}_{kl}$ above is evaluated on each edge connecting mesh points $i$ and $j$ . This is given by a Galerkin approximation stabilized by edge-based upwinding:

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

where

\[ \begin{split} F^i_{kl} &= v^i_k v^i_l + p^i\delta_{kl} \\ v^{ij} &= \frac{1}{2}\frac{D^{ij}_k}{\left|D^{ij}_k\right|}(v^i_k + v^j_k) \end{split} \]

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

Discretization of pressure correction

The pressure correction step is carried out by solving the Poisson equation for the pressure increment in time, discretized using Galerkin finite elements and linear shapefunctions:

\[ -\sum_{\Omega_h \in i} \int_{\Omega_h} \frac{\partial N^i}{\partial x_k} \frac{\partial N^j}{\partial x_k} \mathrm{d}\Omega \left(p^{n+1}-p^n\right)^j = \frac{1}{\Delta t} \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] \]

where $\Omega_h \in i$ represents all elements attached to vertex $i$ . The numerical flux $F^{ij}_k$ is evaluated on each edge connecting mesh points $i$ and $j$ , given by a Galerkin approximation stabilized as

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

where

\[ \begin{split} F^i_k &= v^i_k \\ \lambda^{ij} &= \frac{\Delta t}{\left|l^{ij}_k\right|} \\ l^{ij}_k &= x^i_k - x^j_k \end{split} \]

The above linear system of equations for the pressure increment in time is solved using a preconditioned conjugate gradient method.

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 advective/diffusive step is advanced in time using the $\theta$ -scheme outlined above. $\theta$ =0 yields forward Euler, $\theta$ =1 results in backward Euler, and $\theta$ =0.5 in Crank-Nicholson time integration for the diffusion term. For $\theta>0$ the resulting linear system is solved numerically using a preconditioned conjugate gradient method.

In addition to integrating the diffusion term implicit in time, the advective/diffusive predictor step can be integrated using a multi-stage Runge-Kutta scheme:

\[ \begin{split} \mbox{\boldmath$v$}^i = \mbox{\boldmath$v$}^n - \alpha^i\Delta t\left( \mbox{\boldmath$v$}^{i-1}\cdot\nabla\mbox{\boldmath$v$}^{i-1} + \nabla p^n -\nabla\mu\nabla\mbox{\boldmath$v$}^{i-1} \right); \quad i=1,\dots,k-1 \\ \left[\frac{1}{\Delta t} - \theta\nabla\mu\nabla\right](\mbox{\boldmath$v$}^k-\mbox{\boldmath$v$}^n) + \mbox{\boldmath$v$}^{k-1}\cdot\nabla\mbox{\boldmath$v$}^{k-1} + \nabla p^n = \nabla\mu\nabla\mbox{\boldmath$v$}^{k-1} \qquad \end{split} \]

where $\alpha^i = 1/(1+k-i)$ are the Runge-Kutta coefficients for a $k$ -step scheme. This improves accuracy and allows larger time steps.

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.