of 14
Discontinuous Galerkin scheme for elliptic equations on extremely stretched grids
Nils L. Vu
*
Theoretical Astrophysics, Walter Burke Institute for Theoretical Physics,
California Institute of Technology
, Pasadena, California 91125, USA
(Received 9 May 2024; accepted 25 September 2024; published 24 October 2024)
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
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
S
p
ECTRE
numerical
relativity code.
DOI:
10.1103/PhysRevD.110.084062
I. INTRODUCTION
Discontinuous Galerkin (DG) methods are a class of
spectral finite-element numerical methods for solving
partial differential equations. Although most prevalently
applied to hyperbolic conservation laws, DG methods can
also be advantageous for elliptic equations
[1,2]
. Their
spectral elements are able to retain exponential conver-
gence rates
[3]
, and the domain decomposition allows for
efficient parallelization on large computing clusters
[4
6]
.
Together, these properties make DG methods attractive for
computationally demanding elliptic problems in computa-
tional physics, such as binary black hole initial data in
numerical relativity
[7,8]
. High-accuracy initial data is
particularly important to underpin ongoing efforts in the
numerical relativity community to develop a next gener-
ation of massively parallel simulation codes
[9
15]
, spe-
cifically to simulate binary black hole mergers with
sufficient accuracy for the next generation of gravitational
wave observatories
[16
18]
.
However, elliptic problems in the field of numerical
relativity often involve extremely stretched grids to impose
boundary conditions at very large distances, which can be
challenging for DG methods. Specifically, initial data
problems in general relativity with black holes or neutron
stars are typically solved on spherical domains centered on
the objects. Since the solution at the outer boundary of the
sphere is unknown but typically falls off as
1
=r
with
coordinate distance from the center, the boundary condition
is set to empty flat spacetime at a large but finite radius
R
(
approximate asymptotic flatness
)
[19]
. This approach
incurs an error of
O
ð
1
=R
Þ
on the solution, which is reduced
below numerical accuracy by pushing the outer boundary
far enough out, often to
R
10
9
. To resolve this large
computational volume numerically, radial grid points are
also distributed on a
1
=r
scale. This way, a large spherical
shell can be resolved by just a few spectral elements that are
extremely stretched. Figure
1
below illustrates such an
extremely stretched spectral element. While this strategy
has worked well for specialized spectral numerical schemes
in the past
[19]
, generic DG schemes struggle with the
extremely large Jacobian factors that arise from the grid
stretching. Specifically, we found that the DG scheme
developed in Ref.
[5]
fails to converge on the extremely
stretched grids used in numerical relativity. In this article
we construct a primal DG scheme for generic elliptic
equations that is capable of handling such extreme grid
stretching. The scheme is implemented in the open-source
numerical relativity code
S
p
ECTRE
[9]
.
To further alleviate the strain that the extreme grid
stretching can pose on the numerical scheme, we propose
new Robin-type boundary conditions at the outer boundary.
These boundary conditions incur an error of only
O
ð
1
=R
2
Þ
,
so the outer boundary can be placed at a smaller radius,
*
Contact author: nilsvu@caltech.edu
PHYSICAL REVIEW D
110,
084062 (2024)
2470-0010
=
2024
=
110(8)
=
084062(14)
084062-1
© 2024 American Physical Society
R
10
5
, while still maintaining the same level of error. To
this end, we develop new Robin-type boundary conditions
for the extended conformal thin sandwich (XCTS) formu-
lation of the Einstein constraint equations.
Alternative approaches to the extreme grid stretching have
been proposed on the level of the elliptic equations them-
selves. For example, the
LORENE
[20,21]
,
KADATH
[22,23]
,
and
SGRID
[24]
libraries employ multiple computational
domains and a radial coordinate transformation to compac-
tify the outermost domain. This allows them to impose
boundary conditions at spatial infinity rather than at a finite
radius, but requires special treatment of derivative operators
in the different domains to incorporate the coordinate trans-
formations. The
T
wo
P
unctures
code
[25]
employs a similar
compactification on a single domain, rewriting the elliptic
equations in compactified coordinates. Other possible
approaches include conformal compactification methods
and hyperboloidal slicing
[26
28]
.Thesemethodsalso
require special treatment of spatial infinity in the elliptic
equations and have not yet been routinely applied to
construct and evolve binary black hole initial data. With
the numerical scheme developed in this article we avoid such
special transformations of the elliptic equations to deal with
asymptotic boundary conditions at spatial infinity. Instead,
the equations can simply be solved on a sufficiently large
domain such that the error incurred by the finite outer radius
falls below numerical accuracy.
This article is structured as follows. Section
II
formulates
the primal DG scheme for generic elliptic equations.
Section
III
demonstrates the versatility of the DG scheme
by solving a simple Poisson problem and two classic
general-relativistic initial data problems with black holes
on extremely stretched grids. We conclude in Sec.
IV
.
II. DISCONTINUOUS GALERKIN
FORMULATION
The DG formulation developed here applies to any set of
coupled, nonlinear elliptic equations that can be written in
first-order flux form,
i
F
α
i
½
u;
u
;
x
S
α
½
u;
u
;
x
f
α
ð
x
Þ
;
ð
1
Þ
where
F
α
i
are the fluxes,
S
α
are the sources, and
f
α
are
the fixed sources. The index
α
enumerates the set of
equations and variables. The fluxes and sources are
functionals of the variables
u
α
ð
x
Þ
and their first deriva-
tives. The fluxes
F
α
i
must be linear in the variables and
their derivatives, which means they can be written in
the form
F
α
i
¼
A
ij
αβ
j
u
β
þ
B
i
αβ
u
β
;
ð
2
Þ
where
A
ij
αβ
must be positive definite to ensure ellipticity
of the equations. The sources
S
α
can be nonlinear. The
fixed sources
f
α
ð
x
Þ
are independent of the variables.
Equation
(1)
closely resembles formulations of hyper-
bolic conservation laws but allows the fluxes
F
α
i
to be
higher-rank tensor fields. It is a simplification of the first-
order flux form defined in Ref.
[5]
that arises by selecting
the gradients of the variables as first-order auxiliary
variables.
For example, a Poisson equation for the variable
u
ð
x
Þ
on
a curved manifold with metric
g
ij
ð
x
Þ
in Cartesian coor-
dinates is
g
ij
i
j
u
ð
x
Þ¼
f
ð
x
Þ
;
ð
3
Þ
where
i
is the covariant derivative compatible with the
metric
g
ij
. It can be written in first-order flux form
(1)
with
the fluxes and sources
F
i
¼
g
ij
j
u;
S
¼
Γ
i
ij
F
j
;
ð
4
Þ
where
Γ
i
jk
¼
1
2
g
il
ð
j
g
kl
þ
k
g
jl
l
g
jk
Þ
are the Christoffel
symbols associated with the metric
g
ij
. Any Poisson
equation with nonlinear sources, such as the puncture
equation solved in Sec.
III B
for simple binary black hole
initial data, can be written in this form as well.
Also the extended conformal thin-sandwich (XCTS)
formulation of the Einstein constraint equations
[8,29,30]
,
that we solve for binary black hole initial data in Sec.
III C
,
can be written in first-order flux form
(1)
. The XCTS
equations are a set of five coupled, nonlinear elliptic
equations for the variables
f
ψ
;
αψ
;
β
i
g
given by
FIG. 1. Geometry of an extremely stretched wedge-shaped
element
Ω
k
in two dimensions. The coordinate transformation
x
ð
ξ
Þ
deforms a reference cube
ξ
½
1
;
1

2
to a wedge with inner
radius
r
0
and outer radius
r
1
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).
NILS L. VU
PHYS. REV. D
110,
084062 (2024)
084062-2
̄
2
ψ
¼
1
8
ψ
̄
R
þ
1
12
ψ
5
K
2
1
8
ψ
7
̄
A
ij
̄
A
ij
2
πψ
5
ρ
;
ð
5a
Þ
̄
2
ð
αψ
Þ¼
αψ

7
8
ψ
8
̄
A
ij
̄
A
ij
þ
5
12
ψ
4
K
2
þ
1
8
̄
R
þ
2
πψ
4
ð
ρ
þ
2
S
Þ

ψ
5
t
K
þ
ψ
5
β
i
̄
i
K;
ð
5b
Þ
̄
i
ð
̄
L
β
Þ
ij
¼ð
̄
L
β
Þ
ij
̄
i
ln
ð
̄
α
Þþ
̄
α
̄
i
ð
̄
α
1
̄
u
ij
Þ
þ
4
3
̄
αψ
6
̄
j
K
þ
16
π
̄
αψ
10
S
j
;
ð
5c
Þ
with
̄
A
ij
¼
1
2
̄
α
ðð
̄
L
β
Þ
ij
̄
u
ij
Þ
and
̄
α
¼
αψ
6
. The XCTS
equations are formulated on a background geometry defined
by the conformal metric
̄
γ
ij
, which also defines the covariant
derivative
̄
, the Ricci scalar
̄
R
, and the longitudinal
operator
ð
̄
L
β
Þ
ij
¼
̄
i
β
j
þ
̄
j
β
i
2
3
̄
γ
ij
̄
k
β
k
:
ð
6
Þ
The conformal metric is chosen in advance alongside the
trace of the extrinsic curvature
K
, their respective time
derivatives
̄
u
ij
and
t
K
, and the matter sources
ρ
,
S
, and
S
i
.
Appendix lists the fluxes and sources to formulate the XCTS
equations in first-order flux form
(1)
.
Therefore, a generic implementation of the DG dis-
cretization for the first-order flux form
(1)
encompasses a
large class of elliptic equations, including the puncture
and XCTS initial data formulations for black holes. It
also encompasses other (nonlinear) Poisson-type equa-
tions, as well as linear elasticity
[31]
. The versatility of
this generic formulation lowers the barrier of entry for
using the DG method in our code, since new equations
can be added by implementing their fluxes and sources
without having to learn the details of the DG discretiza-
tion scheme. In the following we will derive the DG
discretization scheme without reference to any specific
elliptic equations.
A. Domain decomposition
To formulate the DG scheme we decompose the
computational domain
Ω
into deformed cubes, called
elements
Ω
k
, as described in more detail in Refs.
[4,5]
.
In this article we focus on extremely stretched wedge-
shaped elements as pictured in Fig.
1
from which we build
spherical domains.
To deform the hexahedral elements to wedges we
employ an invertible map from logical coordinates
ξ
¼
ð
ξ
;
η
;
ζ
Þ
½
1
;
1

3
on the reference cube to the coordinates
x
Ω
k
in which the elliptic equations
(1)
are formulated.
The coordinate map also defines the Jacobian matrix
J
i
ˆ
|
x
i
ξ
ˆ
|
ð
7
Þ
with determinant J and inverse
ð
J
1
Þ
ˆ
|
i
¼
ξ
ˆ
|
=
x
i
. The
wedge coordinate map pictured in Fig.
1
is
x
ð
ξ
Þ¼
r
ð
ξ
Þ
ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
1
þ
η
2
þ
ζ
2
p
0
B
@
1
η
ζ
1
C
A
ð
8
Þ
in three dimensions, where
r
ð
ξ
Þ
is the radial coordinate
mapping. We choose either a linear, logarithmic, or inverse
radial mapping,
linear
r
ð
ξ
Þ¼
1
2
ðð
1
ξ
Þ
r
0
þð
1
þ
ξ
Þ
r
1
Þ
;
ð
9a
Þ
logarithmic ln
r
ð
ξ
Þ¼
1
2
ðð
1
ξ
Þ
ln
r
0
þð
1
þ
ξ
Þ
ln
r
1
Þ
;
ð
9b
Þ
inverse
1
r
ð
ξ
Þ
¼
1
2

1
ξ
r
0
þ
1
þ
ξ
r
1

:
ð
9c
Þ
Here,
r
0
and
r
1
denote the inner and outer radii of
the wedge.
For illustration, we can quantify how much the reference
cube is stretched by the wedge map with the inverse radial
mapping, Eq.
(9c)
. First, note that the midpoint of the
reference cube,
ξ
¼
0
, is mapped to only
r
mid
¼
2
r
0
in the
limit
r
1
r
0
. Second, the Jacobian determinant at the outer
radius is J
ð
r
1
Þ
4
=r
0
in the limit
r
1
r
0
. This means that
for
r
0
¼
10
and
r
1
¼
10
9
the midpoint is mapped to only
r
mid
20
and the Jacobian determinant is J
10
35
.OurDG
scheme has to be able to handle such extreme stretching.
Forthe DGdiscretizationwechoosea regulargridof
N
k
¼
Q
3
ˆ
{
¼
1
N
k;
ˆ
{
Legendre-Gauss-Lobatto (LGL) or Legendre-
Gauss (LG) collocation points
ξ
p
¼ð
ξ
p
1
;
η
p
2
;
ζ
p
3
Þ
on the
reference cube of element
Ω
k
, as illustrated in Fig.
1
and
described in more detail in Ref.
[5]
. Fields on the computa-
tional domain are represented by their values at these grid
points. This nodal representation is facilitated by expanding
each field
u
ð
x
Þ
on the grid as
u
ð
x
Þ¼
X
N
k
p
u
p
ψ
p

ξ
ð
x
Þ

;
ð
10
Þ
where
u
p
¼
u
ð
x
ð
ξ
p
ÞÞ
are the nodal coefficients and
ψ
p
ð
ξ
Þ
are the basis functions
ψ
p
ð
ξ
Þ¼
Y
3
ˆ
{
¼
1
l
p
ˆ
{
ð
ξ
ˆ
{
Þ
:
ð
11
Þ
These basis functions are products of the one-dimensional
Lagrange interpolating polynomials
l
p
ˆ
{
ð
ξ
Þ
with respect to
the collocation points
ξ
p
on the reference cube.
DISCONTINUOUS GALERKIN SCHEME FOR ELLIPTIC
...
PHYS. REV. D
110,
084062 (2024)
084062-3
The choice between LGL and LG collocation points is
particularly relevant on extremely stretched grids. On LGL
grids, points are placed on element boundaries, as illus-
trated in Fig.
1
. On LG grids, points are only placed in the
interior of the element. For illustration, in the previous
example of a wedge with an inverse radial mapping
extending from
r
0
¼
10
to
r
1
¼
10
9
, the outermost radial
grid point for LGL collocation is at the outer radius of
r
¼
10
9
, but for LG collocation with six radial grid points it
is only at
r
300
. This can make dense matrix operations
over the grid points within an element numerically better
behaved on LG grids because functions that fall off as
1
=r
vary over fewer orders of magnitude when evaluated at the
LG grid points. However, the lack of grid points at the
element boundaries means that boundary conditions must
be projected into the interior of the element, which can
introduce additional errors. In this article we find that both
LGL and LG collocation points can be used successfully on
extremely stretched grids.
B. DG residuals
To derive the primal DG formulation for the elliptic
equations, we project the flux form, Eq.
(1)
, on the set of
basis functions
ψ
p
to obtain the DG residuals
ð
ψ
p
;
i
F
i
Þ
Ω
k
þð
ψ
p
;
S
Þ
Ω
k
¼ð
ψ
p
;f
Þ
Ω
k
;
ð
12
Þ
where the inner product is defined as
ð
φ
;
π
Þ
Ω
k
Z
Ω
k
φ
ð
x
Þ
π
ð
x
Þ
d
3
x;
ð
13a
Þ
¼
Z
½
1
;
1

3
φ

x
ð
ξ
Þ

π

x
ð
ξ
Þ

Jd
3
ξ
:
ð
13b
Þ
We suppress the index
α
that enumerates equations and
variables in this section. A partial integration of Eq.
(12)
gives the weak form of the DG residuals,
ð
i
ψ
p
;
F
i
Þ
Ω
k

ψ
p
;
ð
n
i
F
i
Þ


Ω
k
þð
ψ
p
;
S
Þ
Ω
k
¼ð
ψ
p
;f
Þ
Ω
k
;
ð
14
Þ
where the inner product on the element boundary is
defined as
ð
φ
;
π
Þ
Ω
k
Z
Ω
k
φ
ð
x
Þ
π
ð
x
Þ
d
2
x;
ð
15a
Þ
¼
Z
½
1
;
1

2
φ

x
ð
ξ
Þ

π

x
ð
ξ
Þ

J
Σ
d
2
ξ
ð
15b
Þ
with the surface Jacobian J
Σ
. We have also introduced the
numerical fluxes
ð
n
i
F
i
Þ

on element boundaries in
Eq.
(14)
, which couple the DG residuals on neighboring
elements as discussed in Sec.
II C
. Compared to Ref.
[5]
we
have simplified the Galerkin projections to be performed
over coordinate volume rather than proper volume because
we found no significant difference in performance between
the two.
Now, to discretize the auxiliary first-order equation for
the fluxes
F
i
we employ a primal formulation rather than
the Schur-complement approach used in Ref.
[5]
. This is
because the Schur-complement approach, while simpler to
implement, fails to handle extremely stretched grids. The
two approaches differ only in the boundary correction term
arising from the auxiliary equation for the fluxes
F
i
.We
project the linear decomposition of the fluxes, Eq.
(2)
,on
the derivatives of the basis functions to find an expression
for the first term in Eq.
(14)
, and integrate by parts, to
obtain
ð
i
ψ
p
;
F
i
Þ
Ω
k
¼ð
i
ψ
p
;
A
ij
j
u
Þ
Ω
k
þð
i
ψ
p
;
B
i
u
Þ
Ω
k
;
ð
16a
Þ
¼
ð
j
A
ij
i
ψ
p
;u
Þ
Ω
k
þð
i
ψ
p
;
B
i
u
Þ
Ω
k
þð
i
ψ
p
;
A
ij
n
j
u

Þ
Ω
k
;
ð
16b
Þ
where
u

is the auxiliary numerical flux discussed in
Sec.
II C
. With another partial integration we find
ð
i
ψ
p
;
F
i
Þ
Ω
k
¼ð
i
ψ
p
;
A
ij
j
u
Þ
Ω
k
þð
i
ψ
p
;
B
i
u
Þ
Ω
k
þ

i
ψ
p
;
A
ij
n
j
ð
u

u
Þ

Ω
k
;
ð
17a
Þ
¼ð
i
ψ
p
;
F
i
Þ
Ω
k
þ

i
ψ
p
;
A
ij
n
j
ð
u

u
Þ

Ω
k
:
ð
17b
Þ
Inserting this result into Eq.
(14)
we find the primal DG
residuals in weak form,
ð
i
ψ
p
;
F
i
Þ
Ω
k
þð
i
ψ
p
;
A
ij
n
j
ð
u

u
ÞÞ
Ω
k
ð
ψ
p
;
ð
n
i
F
i
Þ

Þ
Ω
k
þð
ψ
p
;
S
Þ
Ω
k
¼ð
ψ
p
;f
Þ
Ω
k
;
ð
18
Þ
where
u

and
ð
n
i
F
i
Þ

are the numerical fluxes discussed in
the next section. With a final partial integration we find the
primal DG residuals in strong form,
ð
ψ
p
;
i
F
i
Þ
Ω
k
þ

i
ψ
p
;
A
ij
n
j
ð
u

u
Þ

Ω
k

ψ
p
;
ð
n
i
F
i
Þ

n
i
F
i

Ω
k
þð
ψ
p
;
S
Þ
Ω
k
¼ð
ψ
p
;f
Þ
Ω
k
:
ð
19
Þ
Written in terms of discrete matrices over the element
s grid
points, as detailed in Ref.
[5]
, the primal DG residuals in
weak form are
MD
T
i
·
F
i
þ
MLD
T
i
·
A
ij
n
j
ð
u

u
Þ
ML
·
ð
n
i
F
i
Þ

þ
M
·
S
¼
M
·
f;
ð
20
Þ
NILS L. VU
PHYS. REV. D
110,
084062 (2024)
084062-4
and in strong form they are
MD
i
·
F
i
þ
MLD
T
i
·
A
ij
n
j
ð
u

u
Þ
ML
·
ðð
n
i
F
i
Þ

n
i
F
i
Þþ
M
·
S
¼
M
·
f;
ð
21
Þ
where
F
i
¼
A
ij
ð
J
1
Þ
ˆ
|
j
D
ˆ
|
·
u
þ
B
i
u
. The symbol · denotes
matrix-vector multiplication over grid points. We have
defined here the mass matrix,
M
pq
¼
δ
pq
J
j
p
Y
3
ˆ
|
¼
1
w
p
ˆ
|
;
ð
22
Þ
the weak and strong stiffness matrices in the volume of the
element,
MD
T
i;pq
¼
D
T
ˆ
{;pq
J
j
q
ð
J
1
Þ
ˆ
{
i
j
q
Y
3
ˆ
|
¼
1
w
q
ˆ
|
;
ð
23
Þ
MD
i;pq
¼
D
ˆ
{;pq
J
j
q
ð
J
1
Þ
ˆ
{
i
j
q
Y
3
ˆ
|
¼
1
w
p
ˆ
|
;
ð
24
Þ
as well as the lifting matrices from the element boundary to
the volume of the element,
ML
p
̄
q
¼
I
p
̄
q
J
Σ
j
̄
q
Y
2
ˆ
|
¼
1
w
̄
q
ˆ
|
;
ð
25
Þ
MLD
T
i;p
̄
q
¼
D
T
ˆ
{;pq
I
q
̄
q
J
Σ
j
̄
q
ð
J
1
Þ
ˆ
{
i
j
̄
q
Y
2
ˆ
|
¼
1
w
̄
q
ˆ
|
:
ð
26
Þ
To obtain these matrices we have evaluated the integrals in
Eq.
(18)
over the reference cube using Gauss-Lobatto
quadrature (on LGL grids) or Gauss quadrature (on LG
grids)
[5]
. The coefficients
w
p
ˆ
|
are the quadrature weights
in dimension
ˆ
|
of the reference cube. Indices with a bar, i.e.,
̄
q
, only enumerate grid points on the boundary of the
element. The boundary points are
lifted
to the volume of
the element with the interpolation matrices
I
p
̄
q
¼
l
p
ˆ
{
ð
1
Þ
Y
ˆ
|
ˆ
{
δ
p
ˆ
|
̄
q
ˆ
|
ð
27
Þ
for a face in dimension
ˆ
{
of the reference cube. The
interpolation matrices reduce to
I
p
̄
q
¼
δ
p
̄
q
for LGL grids,
since boundary points
ξ
ˆ
{
¼
1
are also grid points. The
logical differentiation matrix is defined as
D
ˆ
{;pq
¼
l
0
q
ˆ
{
ð
ξ
p
ˆ
{
Þ
Y
ˆ
|
ˆ
{
δ
p
ˆ
|
q
ˆ
|
:
ð
28
Þ
It takes the derivative along dimension
ˆ
{
on the refer-
ence cube.
For the weak form, Eq.
(20)
, only the
MLD
T
i
term is
different from the Schur-complement approach used in
Ref.
[5]
, Eq. (30b). It involves a lifting operation with the
weak differentiation matrix and only Jacobians evaluated
on the element faces. For the strong form, Eq.
(21)
, the
volume divergence term is different as well, specifically
Eq.
(24)
. It has been shown in Ref.
[32]
that the strong and
weak forms are equivalent when using the metric identities
to move the Jacobian into the strong divergence,
i
F
i
¼ð
J
1
Þ
ˆ
{
i
ˆ
{
F
i
;
ð
29
Þ
¼
1
J
ˆ
{
J
ð
J
1
Þ
ˆ
{
i
F
i
:
ð
30
Þ
We call Eq.
(30)
the
strong logical
formoftheDG residuals
when used in Eq.
(19)
because it involves a divergence in
logical coordinates. This form leads to the strong logical
stiffness matrix, Eq.
(24)
, which differs from the strong
stiffness matrix in Ref.
[5]
,Eq.
(41)
, by the ordering of
Jacobian operations. Both the weak form and the strong
logical form work equally well on extremely stretched grids,
as they are numerically equivalent. However, the strong form
without applying the metric identities, Eq.
(29)
, as used in
Ref.
[5]
, does not work on extremely stretched grids. We
suspect that the mixing of derivatives through the Jacobian in
Eq.
(29)
is numerically ill behaved on such grids.
C. Generalized internal penalty flux
For the numerical fluxes in the DG residuals
(18)
we
choose the generalized internal penalty flux introduced in
Ref.
[5]
. Using the simplified first-order flux form, Eq.
(1)
,
the generalized internal penalty flux is
u

¼
1
2
ð
u
int
þ
u
ext
Þ
;
ð
31a
Þ
ð
n
i
F
i
Þ

¼
1
2
n
i
ð
F
int
i
þ
F
ext
i
Þ
σ
n
i
A
ij
n
j
ð
u
int
u
ext
Þ
;
ð
31b
Þ
where
int
denotes values on the interior side of a
boundary shared with another element, and
ext
denotes
values on the exterior side that were sent from the
neighboring element. The penalty factor
σ
is
σ
¼
C
ð
max
ð
p
int
;p
ext
Þþ
1
Þ
2
min
ð
h
int
;h
ext
Þ
;
ð
32
Þ
where
p
is the polynomial order perpendicular to the
element boundary, and
h
¼
2
J
=
J
Σ
is a measure of the
element size perpendicular to the element boundary. We
choose the penalty parameter
C
¼
1
.
5
in this article.
DISCONTINUOUS GALERKIN SCHEME FOR ELLIPTIC
...
PHYS. REV. D
110,
084062 (2024)
084062-5
Boundary conditions are imposed through numerical
fluxes as well. They are imposed on the average of the
boundary fluxes to obtain faster convergence
[5,33]
.
Therefore, for Dirichlet-type boundary conditions
u
¼
u
D
on an element boundary we impose
1
2
ð
u
int
þ
u
ext
Þ¼
u
D
ð
33a
Þ
u
ext
¼
2
u
D
u
int
ð
33b
Þ
and
F
ext
i
¼
F
int
i
in the numerical fluxes, Eq.
(31)
.For
Neumann-type boundary conditions
n
i
F
i
¼ð
n
i
F
i
Þ
N
we
instead impose
1
2
n
i
ð
F
int
i
þ
F
ext
i
Þ¼ð
n
i
F
i
Þ
N
ð
34a
Þ
n
i
F
ext
i
¼
2
ð
n
i
F
i
Þ
N
n
i
F
int
i
ð
34b
Þ
and
u
ext
¼
u
int
in the numerical fluxes, Eq.
(31)
.
III. TEST PROBLEMS
We test the DG scheme on three problems with extremely
stretched grids. All test problems are implemented in the
open-source code
S
p
ECTRE
[9]
and solved with the iterative
multigrid-Schwarz preconditioned Newton-Krylov elliptic
solver detailed in Ref.
[6]
.
In addition to demonstrating that the generic DG
scheme can handle extremely stretched grids, we also
demonstrate that we can obtain a measure of the truncation
error in each element and dimension of these extremely
stretched grids, which is an essential ingredient for
adaptive mesh refinement (AMR). To this end we evaluate
the anisotropic relative truncation error estimate detailed in
Ref.
[34]
, Eq. (57), which has been used in the
S
p
EC
code
for a long time
[35]
. In logical dimension
ˆ
{
of an element
Ω
k
it is given by
T
½
P
q
ˆ
{
log
10
ð
max
ð
P
1
;P
2
ÞÞ
P
N
k;i
q
ˆ
{
¼
1
log
10
ð
P
q
ˆ
{
Þ
W
q
ˆ
{
P
N
k;
ˆ
{
q
ˆ
{
¼
1
W
q
ˆ
{
;
ð
35
Þ
with the exponential weights
W
q
ˆ
{
¼
exp

ð
q
ˆ
{
N
k;
ˆ
{
þ
1
=
2
Þ
2

:
ð
36
Þ
Here,
P
q
ˆ
{
are the
power monitors
in logical dimension
ˆ
{
,
which are the spectral modes marginalized over the two
other logical dimensions
ˆ
|;
ˆ
k
ˆ
{
. They are given by
[34]
P
q
ˆ
{
¼
ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
N
k;
ˆ
{
N
k
X
N
k
p
j
̃
u
q
j
2
δ
p
ˆ
{
q
ˆ
{
v
u
u
t
;
ð
37
Þ
where the spectral modes
̃
u
q
¼
V
1
qp
u
p
are obtained from
the nodal data
u
p
using the Vandermonde matrix
V
pq
¼
Q
ˆ
{
Φ
q
ˆ
{
ð
ξ
p
ˆ
{
Þ
, and
Φ
q
ˆ
{
ð
ξ
Þ
are the normalized Legendre
polynomials. Given the relative truncation error estimate
in each dimension, the
p
-AMR criterion increases the
number of grid points in logical dimension
ˆ
{
of an element
Ω
k
, and hence the polynomial order of the expansion
(10)
,
if any tensor component of the variables
u
satisfies the
condition
10
T
½
P
q
ˆ
{

max
ðj
u
p
>
T
target
;
ð
38
Þ
where
T
target
is a given target absolute truncation error. For
further details see Ref.
[34]
, Sec. V.1. We show that this
p
-AMR criterion works well even on the extremely
stretched DG grids, and that we do not need to use the
compactified grid that was necessary in Ref.
[4]
.
A. Poisson problem with Lorentzian solution
First, we test the DG scheme on a simple Poisson
problem with a solution that falls off as
1
=r
at large
distances. We choose the Lorentzian solution
u
analytic
ð
x
Þ¼ð
1
þ
r
2
Þ
1
=
2
;
ð
39
Þ
where
r
¼
ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
x
2
þ
y
2
þ
z
2
p
is the coordinate distance from
the origin. To solve for this solution we choose the fixed
source
f
ð
x
Þ¼
2
u
analytic
ð
x
Þ¼
3
ð
1
þ
r
2
Þ
5
=
2
and solve
for
u
ð
x
Þ
.
This problem was also solved in Ref.
[4]
. In fact, the
generic weak-primal DG scheme developed here reduces to
the scheme used in Ref.
[4]
for a flat-space Poisson
equation, except for the dealiasing techniques employed
in Ref.
[4]
that we have not found to be necessary so far,
and some small differences in the definition of the penalty
factor
σ
.
We choose a spherical domain with outer radius
R
¼
10
9
and impose the solution, Eq.
(39)
, as Dirichlet boundary
conditions on the outer boundary. The domain is illustrated
in Fig.
2
. It consists of two spherical shells wrapped around
a central cube. The spherical shells are composed of six
wedges, each with a grid of LGL collocation points, as
pictured in Fig.
1
. LG collocation points work equally well
for this problem (not shown). The wedges in the inner shell
transition from a cubical inner boundary to a spherical outer
boundary with a linear radial mapping, and the wedges in
the outer shell employ an inverse radial mapping to stretch
the grid to the outer boundary at radius
10
9
.
Figure
3
shows that the DG scheme achieves the
expected exponential convergence with increasing poly-
nomial order
P
for multiple refinement levels
L
. The
refinement levels are constructed from the domain in
Fig.
2
by splitting each of the 13 elements in two along
each dimension, so at
L
¼
1
it has 104 elements, and at
NILS L. VU
PHYS. REV. D
110,
084062 (2024)
084062-6
L
¼
2
it has 832 elements. At every refinement level we
either increase the polynomial order in each element and
dimension from
P
¼
1
to
P
¼
10
uniformly, or we increase
the polynomial order adaptively and anisotropically using
the
p
-AMR criterion, Eq.
(38)
, with
T
target
¼
10
7
for
L
¼
1
and
T
target
¼
10
9
for
L
¼
2
. The error converges exponen-
tially, and the
p
-AMR criterion is able to handle the extreme
grid stretching without the compactified grid that was
necessary in Ref.
[4]
.
B. Puncture initial data with black holes
Next, we solve a classic puncture problem for binary
black hole initial data
[36]
. This method is widely used in
the numerical relativity community to generate initial data
because of its simplicity. It is characterized by the puncture
equation,
Δ
u
¼
β
ð
α
ð
1
þ
u
Þþ
1
Þ
7
;
ð
40
Þ
where
Δ
is the flat-space Laplacian and
u
ð
x
Þ
is the puncture
field that we solve for. The background fields
α
ð
x
Þ
and
β
ð
x
Þ
are given by
[37]
1
α
¼
X
I
M
I
2
r
I
;
β
¼
1
8
α
7
̄
A
ij
̄
A
ij
;
ð
41a
Þ
̄
A
ij
¼
3
2
X
I
1
r
2
I

2
P
ð
i
I
n
j
Þ
I
ð
δ
ij
n
i
I
n
j
I
Þ
P
k
I
n
k
I
þ
4
r
I
n
ð
i
I
ε
j
Þ
kl
S
k
I
n
l
I

ð
41b
Þ
to represent any number of black holes with masses
M
I
,
positions
C
I
, linear momenta
P
I
, and spins
S
I
. Here,
r
I
¼
k
x
C
I
k
is the Euclidean coordinate distance to the
I
th
black hole and
n
I
¼ð
x
C
I
Þ
=r
I
is the radial unit normal to
the
I
th black hole. Since the puncture equation is just a
Poisson equation with a nonlinear source, it is easily
represented by our generic DG scheme.
The puncture equation is an example of a problem with a
nonsmooth solution, which poses a challenge for spectral
methods. Specifically, the puncture field is only
C
2
con-
tinuous at the punctures
[36]
, so spectral methods are not
expected to converge exponentially in elements that contain
a puncture. Therefore, the most prominent implementation
of the puncture method in the community, the
T
wo
P
unctures
code
[25]
, uses a special coordinate transformation to
achieve rapid (but subexponential) convergence despite
the nonsmooth source. Here we achieve exponential con-
vergence with our generic DG scheme by an
hp
refinement
strategy, which presents an alternative and more flexible
approach to the problem. We have solved a similar problem
FIG. 2. Spherical domain used in the Lorentzian problem,
Sec.
III A
. 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.
FIG. 3. Convergence of the DG scheme with increasing
polynomial order
P
¼
1
to
P
¼
10
for two refinement levels
L
in the Lorentzian problem, Sec.
III A
. The error is measured as
an
L
2
norm over all grid points of the difference to the exact
solution, Eq.
(39)
. 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).
DISCONTINUOUS GALERKIN SCHEME FOR ELLIPTIC
...
PHYS. REV. D
110,
084062 (2024)
084062-7
in Ref.
[4]
with a prototype code that was limited to
Poisson-type equations.
We employ the spherical domain already used in the
Lorentzian problem, Sec.
III A
, which is pictured in Fig.
2
.
For the puncture problem we place the outer radius at only
R
¼
10
5
and employ Robin boundary conditions,
r
ð
ru
Þ¼
0
:
ð
42
Þ
They are implemented as Neumann-type boundary con-
ditions
ð
n
i
F
i
Þ
N
¼
u=r
in Eq.
(34)
.
We can place any number of punctures on the grid
to evaluate the background fields, Eq.
(41)
. Here we
choose the same configuration that is solved in Ref.
[25]
,
Sec. V, which consists of two orbiting punctures placed at
x

¼ð
3
;
0
;
0
Þ
with masses
M

¼
0
.
5
, linear momenta
P

¼ð
0
;

0
.
2
;
0
Þ
, and spins
S

¼
0
. To avoid placing the
punctures on element boundaries, we rotate the domain
such that the
x
axis pierces the central cube diagonally and
set the length of the diagonal to
6
k
x

18
(see Fig.
4
).
This precaution ensures that the
C
2
punctures always lie
within elements, demonstrating that the DG scheme can
indeed handle such nonsmooth solutions. Placing the
punctures on element boundaries works equally well, as
long as the punctures are not placed exactly on a grid point
where the background quantities in Eq.
(41)
diverge.
To demonstrate exponential convergence with our DG
scheme we increase the resolution of the domain iteratively
as follows: in every refinement step we split elements that
contain a puncture in two in each dimension (
h
refinement),
and increase the polynomial degree of all other elements by
one in each dimension (
p
refinement). During this pro-
cedure we always impose a two-to-one balance at element
interfaces, meaning that elements may only have two
neighbors per dimension on every side. We enforce the
two-to-one balance by splitting elements at interfaces
where the condition is violated. This strategy leads to
the domain pictured in Fig.
4
after a few refinement steps,
and to the exponential convergence shown in Fig.
5
. The
error at the location of the punctures is also plotted in
Fig.
6
. It decreases exponentially with the expected con-
vergence rate of
ð
N
points
Þ
1
=
4
[38]
. Small variations in the
convergence rate are expected because enforcing the
two-to-one balance at element interfaces leads to some
h-refinement spilling into the wedge-shaped elements
surrounding the central cube (see Fig.
4
).
In addition to the uniform refinement strategy, we also
employ the
p
-AMR criterion, Eq.
(38)
, with
T
target
¼
10
8
to control the polynomial degree in elements that do not
contain a puncture. As in the Lorentzian problem, Sec.
III A
,
we find that the
p
-AMR criterion performs well on the
stretched grids without the compactified grid that was
necessary in Ref.
[4]
. The full
hp
-AMR convergence study
FIG. 4. Spherical domain used in the puncture problem,
Sec.
III B
, after six
hp
refinement steps. The domain is rotated
such that the
x
axis with the punctures at
x
¼
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.
FIG. 5. Top: solution to the puncture problem on the
x
axis at
the highest resolution. The punctures are placed at
x
¼
3
on this
axis. The solution is interpolated to the
x
axis using the Lagrange
basis within each element,
(10)
. 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 Fig.
6
.
NILS L. VU
PHYS. REV. D
110,
084062 (2024)
084062-8
with the 14 refinement levels depicted in Fig.
6
completed in
7
.
5
min on one node of the
C
altech
HPC
computing cluster at
Caltech, running on 56 Intel Xeon Platinum 8276 CPU cores
clocked at
2
.
2
GHz.
C. Binary black hole XCTS initial data
Finally, we test the DG scheme on a binary black hole
initial data problem in the XCTS formulation, Eq.
(5)
. The
XCTS formulation of the Einstein constraint equations is
more challenging to solve than the puncture formulation
but it has several advantages. For example, it provides a
way to impose a quasiequilibrium condition on the initial
data by a choice of
̄
u
ij
¼
0
and
t
K
¼
0
and then yields
appropriate initial gauge choices for the evolution in form
of the lapse
α
ð
x
Þ
and shift
β
i
ð
x
Þ
[8]
. The XCTS formulation
is also equipped to generate initial data for black holes with
high spins, which is difficult to achieve with the puncture
method
[39
41]
. To this end we follow the superposed
Kerr-Schild formalism
[42]
to set the conformal metric and
the trace of the extrinsic curvature to the superpositions
̄
γ
ij
¼
δ
ij
þ
X
2
I
¼
1
e
r
2
I
=w
2
I
ð
γ
I
ij
δ
ij
Þ
;
ð
43a
Þ
K
¼
X
2
n
¼
1
e
r
2
I
=w
2
I
K
I
;
ð
43b
Þ
where
γ
I
ij
and
K
I
are the spatial metric and the trace of the
extrinsic curvature of two isolated Kerr black holes in Kerr-
Schild coordinates, with
r
I
the Euclidean coordinate dis-
tance from either center. The superpositions are modulated
by two Gaussians with widths
w
I
so that the influence of
either of the two isolated solutions at the position of the
other is strongly damped. The Gaussians also avoid
logarithmic terms in the solution far away from the center
where we employ the inverse radial coordinate map-
ping, Eq.
(9c)
.
Orbital motion is imposed on the two black holes by a
split of the shift in a background and an excess part
[43]
,
β
i
¼
β
i
bg
þ
β
i
ex
:
ð
44
Þ
We choose the background shift
β
i
bg
¼ð
Ω
0
×
x
Þ
i
þ
̇
a
0
x
i
þ
v
i
0
;
ð
45
Þ
where
Ω
0
is the orbital angular velocity,
̇
a
0
is a radial
expansion velocity to control the eccentricity of the binary,
and
v
i
0
is a constant velocity to control its linear momen-
tum
[42,44]
. Then, we insert Eq.
(44)
into the XCTS
equations, Eq.
(5)
, and henceforth solve them for
β
i
ex
instead of
β
i
. This split is necessary on the stretched grid to
avoid numerical error in the derivatives of the shift far
away from the center where
β
i
bg
r
. Instead, derivatives of
the background shift are computed analytically and only
derivatives of the variable
β
i
ex
are computed numerically.
In addition, we find that it is essential on these stretched
grids to solve for
ψ
1
and
αψ
1
instead of
ψ
and
αψ
to
avoid numerical error in the derivatives of the conformal
factor and the lapse. This is because numerical cancellation
errors occur in arithmetic operations between field values at
large distances that differ from one only by a small amount.
This effect is further amplified by the conditioning of the
densespectraldifferentiationmatrices,Eq.
(28)
, which couple
field values across the entire stretched grid of an element.
Given these modifications to the XCTS formulation we
solve the equations on the domain illustrated in Fig.
7
. The
domain transitions from two excised spheres that represent
the black holes to a spherical outer boundary at radius
R
.At
the two excision surfaces in the interior of the domain we
impose apparent-horizon boundary conditions
[45]
with an
optional negative expansion
[46]
,
n
k
̄
k
ψ
¼
ψ
3
8
α
n
i
n
j

ð
̄
L
β
Þ
ij
̄
u
ij

ψ
4
̄
m
ij
̄
i
n
j
1
6
K
ψ
3
ψ
3
4
Θ
I
;
ð
46a
Þ
β
i
¼
α
ψ
2
n
i
þ

Ω
I
r
×
ð
x
x
I
Þ

i
þ
n
i
ψ
2
ð
n
j
β
j
I
þ
α
I
Þ
;
ð
46b
Þ
FIG. 6. Convergence of the discretization error at the position
of the punctures (
x
¼
3
in Fig.
5
). The error converges
exponentially with the expected rate of
ð
N
points
Þ
1
=
4
(solid line),
and the
p
-AMR criterion is able to handle the extreme grid
stretching as well (dotted lines).
DISCONTINUOUS GALERKIN SCHEME FOR ELLIPTIC
...
PHYS. REV. D
110,
084062 (2024)
084062-9
where
n
i
is the unit normal to the excision surface pointing
out of the computational domain towards the center
x
I
of the
excision,
̄
m
ij
¼
̄
γ
ij
n
i
n
j
is the conformal surface metric,
and
Ω
I
r
are the rotation parameters that induce spin on the
black hole. The expansion
Θ
I
can be set to zero to impose
that the excision surface
is
an apparent horizon, or to a
negative value obtained by evaluating the isolated Kerr
solution at the location of the excision surface
[46]
. The last
term in Eq.
(46b)
also corresponds to the negative-expansion
boundary conditions and is obtained by evaluating the lapse
α
I
ð
x
Þ
and shift
β
j
I
ð
x
Þ
of the isolated Kerr solution
[46]
.In
addition to Eq.
(46)
we impose the lapse of the isolated Kerr
solution as Dirichlet boundary condition on the excision
surface.
1
At the outer boundary we impose one of two types of
boundary conditions. Either we impose standard asymp-
totically flat Dirichlet boundary conditions,
ψ
¼
1
;
αψ
¼
1
;
β
i
ex
¼
0
;
ð
47
Þ
or we impose new Robin boundary conditions,
r
ð
r
ψ
Þ¼
0
;
r
ð
r
αψ
Þ¼
0
;
r
ð
r
β
i
ex
Þ¼
0
:
ð
48
Þ
Robin-type boundary conditions for the Einstein con-
straints have been proposed as early as Refs.
[47,48]
,
but not for the more modern XCTS formulation. The Robin
boundary conditions on
ψ
and
αψ
are straightforwardly
implemented as Neumann-type conditions on the fluxes in
Eq.
(34)
,
ð
n
i
̄
i
ψ
Þ
N
¼
ψ
=r
and
ð
n
i
̄
i
αψ
Þ
N
¼
αψ
=r
.To
implement the Robin boundary condition for
β
i
ex
we recast
it as a Neumann-type condition on the flux
n
i
F
i
β
ex
¼
n
i
ð
̄
L
β
ex
Þ
ij
as well. To this end, we define the projection
operator
P
i
j
¼
δ
i
j
n
i
n
j
to set the normal component of the
shift gradient,
ð
̄
i
β
j
ex
Þ
N
¼
P
k
i
̄
k
β
j
ex
n
i
β
j
ex
r
;
ð
49
Þ
and then apply the longitudinal operator, Eq.
(6)
, to obtain
ð
n
i
ð
̄
L
β
ex
Þ
ij
Þ
N
¼
n
i
ð
̄
L
β
ex
Þ
ij

β
j
ex
r
þ
n
k
̄
k
β
j
ex

1
3
n
j
n
i

β
i
ex
r
þ
n
k
̄
k
β
i
ex

:
ð
50
Þ
We have assumed conformal flatness at the outer boundary
in this derivation, which is given approximately by the
exponentially damped superposition in Eq.
(43)
.
We solve for an equal-mass and nonspinning black hole
binary in this article. Our implementation also supports
unequal masses and spins, but neither is particularly
relevant to test the extremely stretched grids and outer
boundary conditions. All parameters are listed in the
Supplemental Material
[49]
. We compute a discretization
error
k
u
u
ref
k
2
as the
L
2
norm over all grid points and all
variables
u
α
¼f
ψ
;
αψ
;
β
i
ex
g
of the difference to a high-
resolution reference simulation. For the reference simula-
tion we use a uniform polynomial degree of
P
¼
18
.To
evaluate the reference solution at the grid points of the
lower-resolution domain we use the Lagrange interpolating
polynomial basis within the elements, Eq.
(10)
.
Figure
8
compares the two types of boundary conditions
at the outer boundary. The error is measured as the
difference to the reference simulation with outer radius
R
¼
10
9
and Dirichlet boundary conditions. As expected,
the error incurred by the Dirichlet boundary conditions is of
order
O
ð
1
=R
Þ
, whereas the error incurred by the Robin
boundary conditions is of order
O
ð
1
=R
2
Þ
. Therefore, in the
following we place the outer boundary at
R
¼
10
9
when
using Dirichlet boundary conditions, and at
R
¼
10
5
when
using Robin boundary conditions.
Figure
9
shows the convergence of the DG scheme for
the XCTS problem with increasing polynomial order
P
in
each element. We either increase the polynomial order in all
elements uniformly by one in each dimension (solid line),
FIG. 7. Binary black hole domain used in the XCTS initial data
problem, Sec.
III C
. The domain transitions from two excised
spheres in the center, to two cubes, to an enveloping sphere, all
the way to a spherical outer boundary at a radius of up to
R
¼
10
9
(not shown). The enveloping sphere is split in half once in radial
direction (black ellipse), so the domain shown here has 54
elements in total.
1
Note that the
S
p
EC
code typically imposes the superposed
lapse on the excision surfaces
[46]
, rather than the lapse of the
isolated Kerr solutions.
NILS L. VU
PHYS. REV. D
110,
084062 (2024)
084062-10
or apply the
p
-AMR criterion, Eq.
(38)
, with
T
target
¼
10
9
(dotted line). The discretization error decreases exponen-
tially on the extremely stretched grids, and the
p
-AMR
criterion is able to handle the extreme grid stretching
as well.
The convergence for Dirichlet and Robin boundary
conditions in Fig.
9
is identical and not shown.
However, the smaller outer radius for the Robin boundary
conditions improves the conditioning of the discretized
problem and hence the convergence rate of the iterative
elliptic solver, reducing its time to solution by about 30%.
The full
p
-AMR convergence study with the 14 refine-
ment levels depicted in Fig.
9
completed in
25
min on
one node of the
C
altech
HPC
computing cluster at Caltech,
running on 56 Intel Xeon Platinum 8276 CPU cores
clocked at 2.2 GHz. This performance matches the results
we have demonstrated in Ref.
[6]
despite the increase in
domain radius by many orders of magnitude. Therefore,
we are now able to generate high-accuracy initial data
for binary black holes in the XCTS formulation with
the scalable and parallelizable DG method. Figure
10
FIG. 8. Error incurred by the Dirichlet and Robin outer
boundary conditions in the XCTS problem, Sec.
III C
, with
increasing radius of the outer boundary
R
. The Dirichlet boundary
conditions incur an error of order
O
ð
1
=R
Þ
, whereas the Robin
boundary conditions only incur an error of order
O
ð
1
=R
2
Þ
, thus
allowing to place the outer boundary closer to the center.
FIG. 9. Convergence of the DG scheme with increasing
polynomial order
P
¼
4
to
P
¼
17
for the XCTS problem,
Sec.
III C
. The error converges exponentially (solid line) and
the
p
-AMR criterion is able to handle the extreme grid stretching
as well (dotted lines).
FIG. 10. Solution to the binary black hole initial data problem in the XCTS formulation, Sec.
III C
.
DISCONTINUOUS GALERKIN SCHEME FOR ELLIPTIC
...
PHYS. REV. D
110,
084062 (2024)
084062-11