A multidomain spectral method for solving elliptic equations

We present a new solver for coupled nonlinear elliptic partial differential equations (PDEs). The solver is based on pseudo-spectral collocation with domain decomposition and can handle one- to three-dimensional problems. It has three distinct features. First, the combined problem of solving the PDE, satisfying the boundary conditions, and matching between different subdomains is cast into one set of equations readily accessible to standard linear and nonlinear solvers. Second, touching as well as overlapping subdomains are supported; both rectangular blocks with Chebyshev basis functions as well as spherical shells with an expansion in spherical harmonics are implemented. Third, the code is very flexible: The domain decomposition as well as the distribution of collocation points in each domain can be chosen at run time, and the solver is easily adaptable to new PDEs. The code has been used to solve the equations of the initial value problem of general relativity and should be useful in many other problems. We compare the new method to finite difference codes and find it superior in both runtime and accuracy, at least for the smooth problems considered here.Comment: 31 pages, 8 figure

"Mariage des Maillages": A new numerical approach for 3D relativistic core collapse simulations

We present a new 3D general relativistic hydrodynamics code for simulations of stellar core collapse to a neutron star, as well as pulsations and instabilities of rotating relativistic stars. It uses spectral methods for solving the metric equations, assuming the conformal flatness approximation for the three-metric. The matter equations are solved by high-resolution shock-capturing schemes. We demonstrate that the combination of a finite difference grid and a spectral grid can be successfully accomplished. This "Mariage des Maillages" (French for grid wedding) approach results in high accuracy of the metric solver and allows for fully 3D applications using computationally affordable resources, and ensures long term numerical stability of the evolution. We compare our new approach to two other, finite difference based, methods to solve the metric equations. A variety of tests in 2D and 3D is presented, involving highly perturbed neutron star spacetimes and (axisymmetric) stellar core collapse, demonstrating the ability to handle spacetimes with and without symmetries in strong gravity. These tests are also employed to assess gravitational waveform extraction, which is based on the quadrupole formula.Comment: 29 pages, 16 figures; added more information about convergence tests and grid setu

hp-adaptive discontinuous Galerkin solver for elliptic equations in numerical relativity

A considerable amount of attention has been given to discontinuous Galerkin methods for hyperbolic problems in numerical relativity, showing potential advantages of the methods in dealing with hydrodynamical shocks and other discontinuities. This paper investigates discontinuous Galerkin methods for the solution of elliptic problems in numerical relativity. We present a novel hp-adaptive numerical scheme for curvilinear and non-conforming meshes. It uses a multigrid preconditioner with a Chebyshev or Schwarz smoother to create a very scalable discontinuous Galerkin code on generic domains. The code employs compactification to move the outer boundary near spatial infinity. We explore the properties of the code on some test problems, including one mimicking Neutron stars with phase transitions. We also apply it to construct initial data for two or three black holes

Matrix-free GPU implementation of a preconditioned conjugate gradient solver for anisotropic elliptic PDEs

Many problems in geophysical and atmospheric modelling require the fast solution of elliptic partial differential equations (PDEs) in "flat" three dimensional geometries. In particular, an anisotropic elliptic PDE for the pressure correction has to be solved at every time step in the dynamical core of many numerical weather prediction models, and equations of a very similar structure arise in global ocean models, subsurface flow simulations and gas and oil reservoir modelling. The elliptic solve is often the bottleneck of the forecast, and an algorithmically optimal method has to be used and implemented efficiently. Graphics Processing Units have been shown to be highly efficient for a wide range of applications in scientific computing, and recently iterative solvers have been parallelised on these architectures. We describe the GPU implementation and optimisation of a Preconditioned Conjugate Gradient (PCG) algorithm for the solution of a three dimensional anisotropic elliptic PDE for the pressure correction in NWP. Our implementation exploits the strong vertical anisotropy of the elliptic operator in the construction of a suitable preconditioner. As the algorithm is memory bound, performance can be improved significantly by reducing the amount of global memory access. We achieve this by using a matrix-free implementation which does not require explicit storage of the matrix and instead recalculates the local stencil. Global memory access can also be reduced by rewriting the algorithm using loop fusion and we show that this further reduces the runtime on the GPU. We demonstrate the performance of our matrix-free GPU code by comparing it to a sequential CPU implementation and to a matrix-explicit GPU code which uses existing libraries. The absolute performance of the algorithm for different problem sizes is quantified in terms of floating point throughput and global memory bandwidth.Comment: 18 pages, 7 figure

IMEX evolution of scalar fields on curved backgrounds

Inspiral of binary black holes occurs over a time-scale of many orbits, far longer than the dynamical time-scale of the individual black holes. Explicit evolutions of a binary system therefore require excessively many time steps to capture interesting dynamics. We present a strategy to overcome the Courant-Friedrichs-Lewy condition in such evolutions, one relying on modern implicit-explicit ODE solvers and multidomain spectral methods for elliptic equations. Our analysis considers the model problem of a forced scalar field propagating on a generic curved background. Nevertheless, we encounter and address a number of issues pertinent to the binary black hole problem in full general relativity. Specializing to the Schwarzschild geometry in Kerr-Schild coordinates, we document the results of several numerical experiments testing our strategy.Comment: 28 pages, uses revtex4. Revised in response to referee's report. One numerical experiment added which incorporates perturbed initial data and adaptive time-steppin

Efficient Multigrid Preconditioners for Atmospheric Flow Simulations at High Aspect Ratio

Many problems in fluid modelling require the efficient solution of highly anisotropic elliptic partial differential equations (PDEs) in "flat" domains. For example, in numerical weather- and climate-prediction an elliptic PDE for the pressure correction has to be solved at every time step in a thin spherical shell representing the global atmosphere. This elliptic solve can be one of the computationally most demanding components in semi-implicit semi-Lagrangian time stepping methods which are very popular as they allow for larger model time steps and better overall performance. With increasing model resolution, algorithmically efficient and scalable algorithms are essential to run the code under tight operational time constraints. We discuss the theory and practical application of bespoke geometric multigrid preconditioners for equations of this type. The algorithms deal with the strong anisotropy in the vertical direction by using the tensor-product approach originally analysed by B\"{o}rm and Hiptmair [Numer. Algorithms, 26/3 (2001), pp. 219-234]. We extend the analysis to three dimensions under slightly weakened assumptions, and numerically demonstrate its efficiency for the solution of the elliptic PDE for the global pressure correction in atmospheric forecast models. For this we compare the performance of different multigrid preconditioners on a tensor-product grid with a semi-structured and quasi-uniform horizontal mesh and a one dimensional vertical grid. The code is implemented in the Distributed and Unified Numerics Environment (DUNE), which provides an easy-to-use and scalable environment for algorithms operating on tensor-product grids. Parallel scalability of our solvers on up to 20,480 cores is demonstrated on the HECToR supercomputer.Comment: 22 pages, 6 Figures, 2 Table