Source: Schlottke-Lakemper2019
The discontinuous Galerkin spectral element method (DGSEM) has advantageous properties for problems dominated by wave propagation, in particular the support for high order of accuracy while maintaining a compact numerical stencil. In the following, the general DGSEM is first briefly introduced for cube-shaped elements. Then, the application of the DGSEM to non-conforming meshes, which occur in hierarchical Cartesian grids, is shown in (Extension to non-conforming meshes). A more detailed derivation of the DGSEM in two dimensions is available, e.g., in Schlottke-Lakemper2017.
Derivation for cube-shaped elements
A general system of hyperbolic conservation equations reads
\begin{equation} \mv{U}_t + \mv{\nabla} \cdot \mv{f}(\mv{U}) - \mv{S}(\mv{U}) = 0,\qquad\forall \mv{x} \in \Omega.\label{eqn:conservation_equation} \end{equation}
\(\mv{U}\) is the vector of state quantities, \(\mv{f}\) represents the flux terms, and \(\mv{S}\) denotes the source vector. The computational domain \(\Omega\) is subdivided into elements, which for hierarchical Cartesian grids are cubes. In each element, the differential equation is mapped to the reference element \(E = [-1,1]^3\) for improved efficiency. The transformed system reads \begin{equation} J\mv{U}_t + \mv{\nabla}_{\mv{\xi}} \cdot \mv{f} - J\mv{S} = 0,\label{eqn:dg_transformed} \end{equation} with the divergence operator in reference coordinates \(\mv{\nabla}_{\mv{\xi}} = (\partial/\partial \xi, \partial/\partial \eta, \partial/\partial \zeta)^\top\). The pseudo-Jacobian determinant is given as \(J = h/2\), where \(h\) is the element edge length in physical coordinates.
Next, the differential equation is rewritten in its weak form. \(\eqref{eqn:dg_transformed}\) is multiplied by a test function \(\phi = \phi(\mv{\xi})\) and integrated over the reference element \(E\). Using integration by parts, the weak formulation reads
\begin{equation} \int\limits_E J\mv{U}_t\phi \,d\mv{\xi} + \int\limits_{\partial E} (\mv{f} \cdot \mv{n})^* \phi\,d{\mv{s}} - \int\limits_E \mv{f} \cdot \mv{\nabla_\xi} \phi \,d\mv{\xi} - \int\limits_E J\mv{S}\phi \,d\mv{\xi} = 0,\label{eqn:dg_weak} \end{equation}
where \(\mv{n}\) is the outward facing surface normal vector in the reference system and \(d{\mv{s}}\) is the surface differential on the element boundary \(\partial E\). Since the solutions from adjacent elements are discontinuous at the interface, the normal flux \(\mv{f} \cdot \mv{n}\) on \(\partial E\) is not uniquely defined, as indicated by the superscript \((\cdot)^*\). To obtain a single flux value at element interfaces, e.g. the local Lax-Friedrichs flux formulation can be used.
Inside each element, the solution \(\mv{U}\) is approximated by a polynomial tensor product ansatz with degree \(N\) in each spatial direction \begin{equation} \mv{U}(\mv{\xi}, t) \approx \mv{U}_h(\mv{\xi}, t) = \sum_{i,j,k=0}^N \mv{\hat{U}}_{ijk}(t) \psi_{ijk}(\mv{\xi}).\label{eqn:approx_u} \end{equation} The nodal coefficients \(\mv{\hat{U}}_{ijk}(t)\) are the degrees of freedom that represent the time-dependent solution state, while the basis functions \(\psi_{ijk} \in [-1, 1]^3\) are the products of one-dimensional Lagrange polynomials and are defined by \begin{equation} \psi_{ijk}(\mv{\xi}) = l_i(\xi) l_j(\eta) l_k(\zeta).\label{eqn:dg_basis_function} \end{equation} The Lagrange basis functions \(l_i(\xi)\) of degree \(N\) are defined on the reference interval \(\xi \in [-1, 1]\) by a set of interpolation points \(\{\xi_i\}_{i=0}^N\), which are the Legendre-Gauss nodes. They have the Lagrange property \(l_i(\xi_j) = \delta_{ij}\). The fluxes inside the elements and the source terms are approximated analogously by polynomials defined on the same set of interpolation points.
Following the Galerkin ansatz, the test function \(\phi\) in Eq. \(\eqref{eqn:dg_weak}\) is selected from the same function space as the polynomial basis \(\phi \in \mathcal{P}_N(E)\), i.e., the vector space of polynomials up to degree \(N\) that are defined on the reference element. Choosing \(\phi = \psi_{ijk}\) for \(i,j,k = 0,\ldots,N\) to obtain \((N + 1)^3\) test functions and further replacing the solution, the fluxes and the source term by their approximations, \(\eqref{eqn:dg_weak}\) becomes \begin{equation} \int\limits_E J\mv{U}_{h,t}\psi_{ijk} \,d\mv{\xi} + \int\limits_{\partial E} (\mv{f} \cdot \mv{n})^* \psi_{ijk}\,d \mv{s} - \int\limits_E \mv{f}_h \cdot \mv{\nabla_\xi} \psi_{ijk} \,d\mv{\xi} - \int\limits_E J\mv{S}_{h}\psi_{ijk} \,d\mv{\xi} = 0,\qquad i,j,k=0,\ldots,N.\label{eqn:dg_weak_discretized} \end{equation} To numerically evaluate the integrals, Gauss quadrature rules are used. Integration and interpolation are collocated, i.e., the same Legendre-Gauss nodes are used for the integration as well as for the discretization of the solution and the fluxes.
\(\eqref{eqn:dg_weak_discretized}\) yields one equation per unknown degree of freedom. Inserting the definitions for the approximation of the solution, the fluxes, and the source terms, the semi-discrete DG operator \(\mathcal{L}(\mv{U}, t) = \mv{U}_t\) is obtained. To advance the solution to the next time step, \(\mathcal{L}\) is integrated in time by a low-storage implementation of an explicit five-stage fourth-order Runge-Kutta scheme, which allows large time steps for wave propagation problems. The time step \(\Delta t = \text{CFL}/(2N + 1) \cdot h/\Lambda_\text{max}\) depends on the Courant-Friedrichs-Lewy (CFL) condition and for, e.g. the the acoustic perturbation equations the maximum wave propagation speed \(\Lambda_\text{max} = \max(|\bar{u}|, |\bar{v}|, |\bar{w}|) + c_0\).}
Extension to non-conforming meshes
To support non-conforming elements, i.e., interfaces where a coarse element is adjacent to multiple refined cells, the mortar element method is employed. In the scheme, conformity is restored by adding an intermediate interface between non-conforming cells. Solution data from neighboring element faces is projected to the mortar, where unique surface flux values are computed. Afterwards, the fluxes are projected back to the individual elements.
Since the refinement level difference between neighboring cells in the Cartesian grid is limited to one the non-conforming case is restricted to interfaces with one element on one and four elements on the other side, such as shown in the figure below.
Without loss of generality for such meshes, the following discussion will be limited to a single coarse element \(E^c\) on the left side and four refined elements \(E^{f,st}\) with \(s, t = 0, 1\) on the right side of the non-conforming interface. Four mortars \(\Xi^{st}\) are defined, which coincide with the corresponding faces \(\partial E^{f,st}\) of the refined elements. The superscripts \((\cdot)^s\) and \((\cdot)^t\) specify the relative position in the \(\xi/\eta\)-direction. A value of \(0/1\) indicates the entity located towards the negative/positive coordinate direction. On the element faces and on the mortar elements, the solution states and fluxes are approximated by polynomials of degree \(N\) in each spatial direction. The same nodal basis functions \(\psi_{ij} = l_i l_j\) with \(i,j = 0,\ldots,N\) as in (Derivation for cube-shaped elements) are used. They are defined on the reference interval \([-1, 1]^2\) by the set of Legendre-Gauss nodes \(\{\mv{\xi}_{ij}\}_{i,j=0}^N\) using the face-local coordinate systems \(\mv{\xi}^c = (\xi^c, \eta^c)\) and \(\mv{\xi}^{st} = (\xi^s, \eta^t)\).
In the forward projection step, the solutions on the coarse element face \(\mv{U}^c_h(\mv{\xi}^c)\) and on the fine element faces \(\mv{U}^{f,st}_{h}(\mv{\xi}^{st})\) are projected to the mortars by using \(L^2\) projection. For the coarse element, the four projections are given by \begin{equation} \int\limits_{-1}^1\int\limits_{-1}^1\left(\mv{U}^{\Xi,c,st}_h - \mv{U}^c_h\right)\phi\,d \xi^s d \eta^t = 0,\qquad s,t = 0,1,\label{eqn:forward_proj} \end{equation} where \(\mv{U}^{\Xi,c,st}_h(\mv{\xi}^{st})\) are the sought solution states on the mortars. As in the previous section, the test function \(\phi\) is replaced by the \((N+1)^2\) Lagrange basis functions \(\psi_{ij}\) and interpolation is substituted by Gauss quadrature. After simplifying \eqnref{forward_proj}, the projections read \begin{equation}\label{eqn:forward_proj_tensor} \mv{\hat{U}}^{\Xi,c,st}_{ij} = \sum_{l,m=0}^N \mv{\hat{U}}^c_{lm} l_l\left(\frac{1}{2}(\xi^s_i \pm 1)\right) l_m\left(\frac{1}{2}(\eta^t_j \pm 1)\right),\qquad \begin{aligned} i,j &= 0,\ldots,N,\ s,t &=0,1, \end{aligned} \end{equation} where the hat \(\hat{(\cdot)}\) again denotes the coefficients for the nodal basis functions. The minus sign corresponds to \(s = 0\) or \(t = 0\), while the plus sign corresponds to \(s = 1\) or \(t = 1\). Since the mortars coincide with the fine element faces, the solution states on the right side of the mortars are equal to those on the fine element faces, i.e., \(\mv{\hat{U}}^{\Xi,f,st}_{ij} = \mv{\hat{U}}^{f,st}_{ij}\).
After obtaining the solutions on the mortar elements, the normal flux is calculated at the Gauss nodes using the numerical flux formulation introduced in (Derivation for cube-shaped elements). The mortar fluxes \(\mv{f}^{*,\Xi,st}_h(\mv{\xi}^{st})\) are then projected back onto the element faces using \(L^2\) projection. For the coarse element, the projection is given by
\begin{align} &\phantom{+}\int\limits_{-1}^0\int\limits_{-1}^0 \left(\mv{f}^{*,c}_h - \mv{f}^{*,\Xi,00}_h\right) \phi \,d\xi^c d\eta^c + \int\limits_{ 0}^1\int\limits_{-1}^0 \left(\mv{f}^{*,c}_h - \mv{f}^{*,\Xi,01}_h\right) \phi \,d\xi^c d\eta^c \nonumber \\ &+\int\limits_{-1}^0\int\limits_{ 0}^1 \left(\mv{f}^{*,c}_h - \mv{f}^{*,\Xi,10}_h\right) \phi \,d\xi^c d\eta^c + \int\limits_{ 0}^1\int\limits_{ 0}^1 \left(\mv{f}^{*,c}_h - \mv{f}^{*,\Xi,11}_h\right) \phi \,d\xi^c d\eta^c = 0, \label{eqn:reverse_proj} \end{align}
where \(\mv{f}^{*,c}_h(\mv{\xi}^c)\) is the flux on the coarse element face. Following the same procedure as before, \(\eqref{eqn:reverse_proj}\) is simplified to \begin{equation} \mv{\hat{f}}^{*,c}_{ij} = \sum_{s,t=0}^1\sum_{p,q=0}^N \mv{\hat{f}}^{*,\Xi,st}_{pq} \frac{1}{2}l_i\left(\frac{1}{2}(\xi^s_p \pm 1)\right)\frac{\omega_p}{\omega_i} \frac{1}{2}l_j\left(\frac{1}{2}(\eta^t_q \pm 1)\right)\frac{\omega_q}{\omega_j}, \qquad i,j=0,\ldots,N,\label{eqn:reverse_proj_tensor} \end{equation} where \(\omega_i\) are the \(N+1\) integration weights for Legendre-Gauss quadrature. For the refined elements, the reverse projection operator is again just the identity, i.e., \(\mv{\hat{f}}^{*,f,st}_{ij} = \mv{\hat{f}}^{*,\Xi,st}_{ij}\).
References
- M. Schlottke-Lakemper, A. Niemoeller, M. Meinke, W. Schroeder, Efficient parallelization for volume-coupled multiphysics simulations on hierarchical Cartesian grids, Computer Methods in Applied Mechanics and Engineering, Volume 352, pp 461-487, 2019, 10.1016/j.cma.2019.04.032.
- M. Schlottke-Lakemper, H. Yu, Sven Berger, M. Meinke, W. Schroeder, A fully coupled hybrid computational aeroacoustics method on hierarchical Cartesian meshes, Comput. Fluids 144 (2017) 137–153, 10.1016/j.compfluid.2016.12.001.
