Single-material DG hydrodynamics

Contents

Inciter supports multiple hydrodynamics schemes. This page describes the DG method for single-material flows.

The discontinuous Galerkin (DG) method implemented in the code is a high order finite element method developed for solving the conservation laws. Similar to classical finite element method, the DG method can achieve high order accuracy by approximating the numerical solution within the element as high order polynomials and admitting discontinuties at the cell interfaces.

Governing equations of the compressible flows

The governing equations used in the DG algorithm can be represented as

\[ \begin{split} \frac{\partial{\boldsymbol{U}}}{\partial t} + \frac{\partial{\boldsymbol{F_k}}}{\partial{\boldsymbol{x_k}}} = \boldsymbol{S} \end{split} \]

For the scalar advection equation, the scalar variable $\boldsymbol{U}$ and $\boldsymbol{F}$ are defined as

\[ \begin{split} \boldsymbol{U} = \phi, \;\boldsymbol{F} = \boldsymbol{c}\phi \end{split} \]

Here $\boldsymbol{c}$ represents the wave speed.

While for the compressible Euler equations, the conservative vector $\boldsymbol{U}$ and the flux vector $\boldsymbol{F}$ are defined as

\[ \begin{split} \boldsymbol{U} = \begin{bmatrix} \rho \\ \rho u_i \\ \rho E \end{bmatrix}, \quad \boldsymbol{F} = \begin{bmatrix} \rho u \\ \rho u_i u_j + p \delta_{ij} \\ (\rho E+p) u_j \end{bmatrix} \end{split} \]

Here $\rho$ , $p$ , and $E$ denote the density, pressure, and specific total energy of the fluid respectively, and $u_i$ is the velocity of the flow in the coordinate direction $x_i$ .The pressure for ideal gas can be computed from the equation of state

\[ \begin{split} p = (\gamma -1) \rho (E - \frac{1}{2} u_j u_i) \end{split} \]

where $\gamma$ is the ratio of specific heats.

Discontinuous Galerkin discretization

In order to apply the discontinuous Galerkin discretization, a weak formulation of the governing equations over the computational domain $\Omega$ is obtained as

\[ \begin{split} \int_{\Omega} \frac{\partial{\boldsymbol{U}}}{\partial t} \boldsymbol{W} d\Omega + \int_{\Gamma}\boldsymbol{F}_k(\boldsymbol{U}) \boldsymbol{n}_k \boldsymbol{W} d\Gamma - \int_{\Omega} \boldsymbol{F}_k(\boldsymbol{U}) \frac{\partial{\boldsymbol{W}}}{\partial x_k} d\Omega = \int_{\Omega} \boldsymbol{S} \boldsymbol{W} d\Omega \end{split} \]

where $\Gamma=\partial{\Omega}$ denotes the boundary of $\Omega$ , $\boldsymbol{n}_k$ is the unit outward vector to the boundary and $\boldsymbol{W}$ is the test function.

By subdividing the domain $\Omega$ into a collection of non-overlapping elements $\Omega_e$ and considering functions U and W defined within each element, we obtain the following semi-discrete formulation,

\[ \begin{split} \int_{\Omega_e} \frac{\partial{\boldsymbol{U}_h}}{\partial t} \boldsymbol{W}_h d\Omega + \int_{\Gamma_e}\boldsymbol{F}_k(\boldsymbol{U}_h) \boldsymbol{n}_k \boldsymbol{W}_h d\Gamma - \int_{\Omega_e} \boldsymbol{F}_k(\boldsymbol{U}_h) \frac{\partial{\boldsymbol{W}_h}}{\partial x_k} d\Omega = \int_{\Omega} \boldsymbol{S} \boldsymbol{W}_h d\Omega \end{split} \]

where $\Gamma_e=\partial{\Omega_e}$ denotes the boundary of $\Omega_e$ , $\boldsymbol{U}_h$ and $\boldsymbol{W}_h$ represent the piecewise polynomial approximations to the analytical solution and test function. The Galerkin method assumes the test function to be equal to the basis function. Then the above equation becomes the following system with $n$ equations:

\[ \begin{split} \frac{d}{dt}\int_{\Omega_e} \boldsymbol{U}_h B_i d\Omega + \int_{\Gamma_e}\boldsymbol{F}_k(\boldsymbol{U}_h) \boldsymbol{n}_k B_i d\Gamma - \int_{\Omega_e} \boldsymbol{F}_k(\boldsymbol{U}_h) \frac{\partial B_i}{\partial x_k} d\Omega = \int_{\Omega} \boldsymbol{S} B_i d\Omega, \quad i=1,2,\dots,n \end{split} \]

The unknows of the conservative variables on each element are defined as

\[ \begin{split} \boldsymbol{U_h} = \sum_{i = 1}^{n} \boldsymbol{U}_j(t)B_j(\boldsymbol{x}) \end{split} \]

where $B_j$ is the basis functions.

Basis function

The Dubiner basis function designed for tetrahedron grid is implemented in the DG algortithm. This set of basis function is defined in the reference tetrahdral domain, which helps to achieve the orthogonality features of the basis function. The projection procedure between the physical domain and reference domain is linear. Therefore, each point in the either domain can be found its corresponding point using linear transport relations. Let us consider a third-order approximation in three-dimensional space. In such case, the basis functions are shown as,

\[ \begin{split} &B_1=1 \\ &B_2=2\xi+\eta+\zeta-1 \\ &B_3=3\eta+\zeta-1 \\ &B_4=4\zeta-1 \\ &B_5=6\xi^2+\eta^2+\zeta^2+6\xi\eta+6\xi\zeta+2\eta\zeta-6\xi-2\eta-2\zeta+1\\ &B_6=5\eta^2+\zeta^2+10\xi\eta+2\xi\zeta+6\eta\zeta-2\xi-6\eta-2\zeta+1 \\ &B_7=6\zeta^2+12\xi\zeta+6\eta\zeta-2\xi-\eta-7\zeta+1 \\ &B_8=10\eta^2+\zeta^2+8\eta\zeta-8\eta-2\zeta+1 \\ &B_9=6\zeta^2+18\eta\zeta-3\eta-7\zeta+1 \\ &B_{10}=15\zeta^2-10\zeta+ 1 \end{split} \]

Numerical flux

Due to the discontinuous function approximations, flux terms are not uniquely defined at element interfaces. In this respect, the flux function $\boldsymbol{F}_k(\boldsymbol{U}_h) $ is replaced by a numerical flux function $\boldsymbol{h}(\boldsymbol{u}^-,\boldsymbol{u}^+)$ . In the current code, the approximate Lax-Friedrich, HLL and HLLC type of Riemann solvers are implemented.

Temporal integration

The semi-discrete form of the governing equations is able to construct a system of ODE equations.

\[ \begin{split} \boldsymbol{M}_{ij}\frac{du_j}{dt} &= \boldsymbol{R}(\boldsymbol{U}_h) \\ \boldsymbol{M}_{ij} &= \int_{\Omega_e}B_i B_j d\Omega \\ \boldsymbol{R}(\boldsymbol{U}_h) &= -\int_{\Gamma_e}\boldsymbol{F}_k (\boldsymbol{U}_h) \boldsymbol{n}_k B_i d \Gamma + \int_{\Omega_e} \boldsymbol{F}_k(\boldsymbol{U}_h) \frac{\partial B_i}{\partial x_k} d\Omega + \int_{\Omega_e} \boldsymbol{S} B_i d\Omega \end{split} \]

The above system is solved using the TVD-RK3 method with the form of

\[ \begin{split} \boldsymbol{U}^{(1)} &= \boldsymbol{U}^n + \Delta t \boldsymbol{M}^{-1} \boldsymbol{R}(\boldsymbol{U}^n) \\ \boldsymbol{U}^{(2)} &= \frac{3}{4} \boldsymbol{U}^n + \frac{1}{4} [\boldsymbol{U}^{(1)} + \Delta t \boldsymbol{M}^{-1} \boldsymbol{R}(\boldsymbol{U}^{(1)})] \\ \boldsymbol{U}^{n+1} &= \frac{1}{3} \boldsymbol{U}^n + \frac{2}{3} [\boldsymbol{U}^{(2)} + \Delta t \boldsymbol{M}^{-1} \boldsymbol{R}(\boldsymbol{U}^{(2)})] \end{split} \]

The TVD-RK3 method is linearly stable for a CFL number less than or equal to $\frac{1}{2p+1}$ ( $p$ is the order of the polynomial solution) and the global time step is obtained by finding the minimum value for all the elements in the computational domain.