Elliptic partial differential equations must be solved numerically for many
problems in numerical relativity, such as initial data for every simulation of
merging black holes and neutron stars. Existing elliptic solvers can take
multiple days to solve these problems at high resolution and when matter is
involved, because they are either hard to parallelize or require a large amount
of computational resources. Here we present a new solver for linear and
non-linear elliptic problems that is designed to scale with resolution and to
parallelize on computing clusters. To achieve this we employ a discontinuous
Galerkin discretization, an iterative multigrid-Schwarz preconditioned
Newton-Krylov algorithm, and a task-based parallelism paradigm. To accelerate
convergence of the elliptic solver we have developed novel
subdomain-preconditioning techniques. We find that our multigrid-Schwarz
preconditioned elliptic solves achieve iteration counts that are independent of
resolution, and our task-based parallel programs scale over 200 million degrees
of freedom to at least a few thousand cores. Our new code solves a classic
black-hole binary initial-data problem faster than the spectral code SpEC when
distributed to only eight cores, and in a fraction of the time on more cores.
It is publicly accessible in the next-generation SpECTRE numerical relativity
code. Our results pave the way for highly-parallel elliptic solves in numerical
relativity and beyond.Comment: 24 pages, 18 figures. Results are reproducible with the ancillary
input file
