Discontinuous Galerkin scheme for elliptic equations on extremely stretched grids

TL;DR

This work has developed a primal DG scheme that is generically applicable to a large class of elliptic equations, including problems on curved and extremely stretched grids, and is able to generate high-quality initial data for binary black hole problems using a parallelizable DG scheme.

Abstract

Discontinuous Galerkin (DG) methods for solving elliptic equations are gaining popularity in the computational physics community for their high-order spectral convergence and their potential for parallelization on computing clusters. However, problems in numerical relativity with extremely stretched grids, such as initial data problems for binary black holes that impose boundary conditions at large distances from the black holes, have proven challenging for DG methods. To alleviate this problem we have developed a primal DG scheme that is generically applicable to a large class of elliptic equations, including problems on curved and extremely stretched grids. The DG scheme accommodates two widely used initial data formulations in numerical relativity, namely the puncture formulation and the extended conformal thin-sandwich (XCTS) formulation. We find that our DG scheme is able to stretch the grid by a factor of $\sim 10^9$ and hence allows to impose boundary conditions at large distances. The scheme converges exponentially with resolution both for the smooth XCTS problem and for the nonsmooth puncture problem. With this method we are able to generate high-quality initial data for binary black hole problems using a parallelizable DG scheme. The code is publicly available in the open-source SpECTRE numerical relativity code.

Figures (10)

• Figure 1: Geometry of an extremely stretched wedge-shaped element $\Omega_k$ in two dimensions. The coordinate transformation $\bm{x}(\bm{\xi})$ deforms a reference cube $\bm{\xi}\in[-1,1]^2$ to a wedge with inner radius $r_0$ and outer radius $r_1 \gg r_0$. Coordinates are mapped with an $1/r$ inverse radial scaling, as illustrated by the light gray lines. Grid points are chosen here as LGL collocation points on the reference cube (black dots).
• Figure 2: Spherical domain used in the Lorentzian problem, \ref{['sec:poisson_lorentzian']}. Two spherical shells wrap a central cube. Each shell is composed of six wedge-shaped elements (black lines) with a grid of LGL collocation points in the element (gray lines). The inner wedges transition from a cubical inner boundary to a spherical outer boundary with a linear radial mapping. The outer wedges stretch to $R=10^9$ with an inverse radial mapping. For this visualization, the inverse radial mapping in the outer shell is undone to be able to see both shells. Note that while the outermost circle of grid points map to $R=10^9$ (outermost black circle), the next-to-outermost circle of grid points shown here already only map to $r = 400$ (outermost gray circle) due to the extreme grid stretching.
• Figure 3: Convergence of the DG scheme with increasing polynomial order $P=1$ to $P=10$ for two refinement levels $L$ in the Lorentzian problem, \ref{['sec:poisson_lorentzian']}. The error is measured as an $L_2$ norm over all grid points of the difference to the exact solution, \ref{['eq:lorentzian']}. The error decreases exponentially with increasing polynomial order $P$ as expected (solid lines), and the $p$-AMR criterion is able to handle the extreme grid stretching as well (dotted lines).
• Figure 4: Spherical domain used in the puncture problem, \ref{['sec:punctures']}, after six $hp$ refinement steps. The domain is rotated such that the $x$ axis with the punctures at $x=\pm 3$ pierces the central cube diagonally. Shown is a slice through the domain spanned by the $x$ axis and another diagonal of the central cube.
• Figure 5: Top: Solution to the puncture problem on the $x$ axis at the highest resolution. The punctures are placed at $x=\pm 3$ on this axis. The solution is interpolated to the $x$ axis using the Lagrange basis within each element, \ref{['eq:lagr_expansion']}. Bottom: Numerical error of the solution at increasing resolution. The error is computed as the difference to the highest resolution shown in the top panel. The error at the punctures is also plotted in \ref{['fig:punctures_convergence']}.