A pseudospectral matrix method for time-dependent tensor fields on a spherical shell

We construct a pseudospectral method for the solution of time-dependent, non-linear partial differential equations on a three-dimensional spherical shell. The problem we address is the treatment of tensor fields on the sphere. As a test case we consider the evolution of a single black hole in numerical general relativity. A natural strategy would be the expansion in tensor spherical harmonics in spherical coordinates. Instead, we consider the simpler and potentially more efficient possibility of a double Fourier expansion on the sphere for tensors in Cartesian coordinates. As usual for the double Fourier method, we employ a filter to address time-step limitations and certain stability issues. We find that a tensor filter based on spin-weighted spherical harmonics is successful, while two simplified, non-spin-weighted filters do not lead to stable evolutions. The derivatives and the filter are implemented by matrix multiplication for efficiency. A key technical point is the construction of a matrix multiplication method for the spin-weighted spherical harmonic filter. As example for the efficient parallelization of the double Fourier, spin-weighted filter method we discuss an implementation on a GPU, which achieves a speed-up of up to a factor of 20 compared to a single core CPU implementation.

💡 Research Summary

The paper presents a novel pseudospectral collocation method for solving time‑dependent, nonlinear partial differential equations (PDEs) on a three‑dimensional spherical shell, with a particular focus on tensor fields. Traditional approaches for tensor fields on a sphere rely on tensor spherical harmonics (spin‑weighted spherical harmonics) which, while mathematically natural, are cumbersome to implement and computationally expensive. The authors propose a different strategy: keep the tensor components expressed in global Cartesian coordinates, but perform spatial discretisation on a spherical grid using a double Fourier (CFF) basis for the angular directions and a Chebyshev basis for the radial direction.

Key elements of the method are:

1. Double Fourier expansion – Both the polar angle θ and the azimuthal angle φ are treated as periodic, enabling the use of Fourier transforms for differentiation. To avoid the pole clustering inherent in a standard latitude‑longitude grid, the authors employ a “double covering” of the θ‑direction, temporarily extending the grid to 2 Nθ points when computing θ‑derivatives, then discarding the redundant half. This allows the derivative matrices to be built once and applied repeatedly via matrix multiplication.

2. Spin‑weighted spherical‑harmonic filter (Yn‑filter) – Because Cartesian tensor components behave like a mixture of different spin‑weights on the sphere, a scalar Y‑filter (projection onto ordinary spherical harmonics) is insufficient and leads to instabilities. The authors construct a filter based on spin‑weighted spherical harmonics that correctly projects each tensor component onto its appropriate spin sector, thereby removing high‑frequency noise, mitigating the pole‑clustering time‑step restriction, and stabilising the nonlinear evolution.

3. Matrix‑multiplication implementation – For the modest problem sizes considered (N ≈ 40–100 points per direction), direct matrix multiplication is faster than FFT‑based approaches. All derivative operators and the Yn‑filter are pre‑computed as dense matrices; during the time‑integration they are applied by simple matrix‑vector products. This approach also simplifies GPU implementation because the same dense‑matrix kernels can be reused without the need for separate FFT libraries.

4. GPU acceleration – Using NVIDIA CUDA, the authors implement the entire spatial operator (derivatives + filter) on the GPU, avoiding costly host‑device memory transfers. The dominant operations are products of small square matrices (e.g., 40 × 40) with larger rectangular matrices (e.g., 40 × 40000). Although cuBLAS is optimised for large square matrices, the authors achieve respectable performance (≈50–100 Gflop/s) for these mixed‑size products, resulting in an overall speed‑up of roughly 10–20× compared with a single‑core CPU implementation.

The method is validated on a demanding test case: the evolution of a single Schwarzschild black hole using the generalized harmonic gauge (GHG) formulation of Einstein’s equations. This system involves 50–58 evolved variables (rank‑1,‑2,‑3 tensors) and exhibits strong non‑linearity but no shocks. The CFF + Yn‑filter scheme produces long‑term stable evolutions comparable to previous work that employed tensor spherical harmonics, while the simpler scalar Y‑filters and other non‑spin‑weighted filters fail to maintain stability. A small residual linear growth is observed over very long times, but it is not explored further.

Performance measurements show that for the chosen resolution (Nθ ≈ 2 k + 1, Nφ = 2 Nθ) the matrix‑based derivative operators are competitive with, and often faster than, FFT‑based derivatives. The pre‑computed filter matrices add negligible overhead, and the overall memory footprint remains modest even with dozens of tensor fields, making the approach suitable for multi‑shell or spectral‑element extensions.

In summary, the paper demonstrates that a double‑Fourier spectral discretisation combined with a spin‑weighted spherical‑harmonic filter, both realised through dense matrix multiplications, offers an efficient and stable framework for evolving tensor PDEs on spherical shells. The technique sidesteps the complexity of tensor spherical‑harmonic expansions while retaining their stabilising benefits, and it scales well on modern GPU hardware. The authors suggest that the method can be extended to multi‑domain spectral‑element schemes and to more complex relativistic simulations involving gravitational waves or other high‑frequency phenomena.