Source-linked AI summary
Efficient Spherical Harmonic Transforms aimed at pseudo-spectral numerical simulations
Nathanaël Schaeffer
TL;DR
Spherical harmonic transforms need efficient, accurate high-resolution computation for pseudo-spectral simulations, but existing implementations face performance, memory, or flexibility constraints. The paper develops explicitly vectorized Gauss-Legendre transforms with on-the-fly Legendre-function computation in SHTns. SHTns is between two and ten times faster than the best alternative, while its lower memory requirements support high-resolution simulations.
Problem
Existing SHT implementations can be limited by overhead, instability, inflexible grids, or excessive memory and computation demands at high resolution.
Method
The paper implements explicitly vectorized Gauss-Legendre scalar and vector transforms using on-the-fly Legendre-function computation, FFTs, and multithreading.
Results
SHTns is between two and ten times faster than the best alternative and achieves efficiency very close to 1 for N ≥ 511.
Takeaways & Limitations
The lower memory requirement and performance gain are expected to enable spectral geodynamo simulations at N = 1023.
Takeaways & Limitations
The tested comparisons used a machine without AVX support, while some alternative libraries could not operate at N = 2047 because of memory limitations.
Abstract
from arXiv · showhide
In this paper, we report on very efficient algorithms for the spherical harmonic transform (SHT). Explicitly vectorized variations of the algorithm based on the Gauss-Legendre quadrature are discussed and implemented in the SHTns library which includes scalar and vector transforms. The main breakthrough is to achieve very efficient on-the-fly computations of the Legendre associated functions, even for very high resolutions, by taking advantage of the specific properties of the SHT and the advanced capabilities of current and future computers. This allows us to simultaneously and significantly reduce memory usage and computation time of the SHT. We measure the performance and accuracy of our algorithms. Even though the complexity of the algorithms implemented in SHTns are in $O(N^3)$ (where N is the maximum harmonic degree of the transform), they perform much better than any third party implementation, including lower complexity algorithms, even for truncations as high as N=1023. SHTns is available at https://bitbucket.org/nschaeff/shtns as open source software.
1 Introduction
Spherical harmonic transforms provide a spectral representation for spherical data and are used in pseudo-spectral simulations. The paper develops an explicitly vectorized Gauss-Legendre implementation to improve transform performance.
- Spherical harmonics form a basis for consistently representing data on the sphere in spectral space.SHTs transform spatial data to the spectral domain and back, analogous to Fourier transforms.
- Fast SHT algorithms can have substantial overhead, numerical instability, or inflexible grid-size requirements.Some methods are not claimed to be effectively faster below N < 512, while others require N + 1 to be a power of 2.
- SpharmonicKit requires 2(N + 1) latitudinal points, whereas Gauss-Legendre quadrature requires only N + 1.The larger grid can make the overall numerical simulation slower even if transform timings were otherwise comparable.
- The authors implement a highly optimized, explicitly vectorized Gauss-Legendre SHT algorithm in the SHTns library.The work targets faster numerical simulations and evaluates the transform against other implementations.
2 Spherical Harmonic Transform (SHT)
Spherical harmonics form an orthonormal basis and Laplace-operator eigenfunctions on the sphere. SHT synthesis and analysis combine Fourier transforms in longitude with Legendre transforms in latitude, using Gauss-Legendre quadrature for exact discrete evaluation at sufficient resolution.
- 2.1 Definitions and properties: Spherical harmonics are indexed by degree n and order m, with −n ≤ m ≤ n, and use normalized associated Legendre polynomials.They are constructed as functions on the sphere and form an orthonormal basis.
- 2.1 Definitions and properties: Spherical harmonics are eigenfunctions of the Laplace operator on the unit sphere, making them useful for physical problems involving that operator.The paper identifies this property as appealing for solving problems in spherical geometry.
- 2.2 Synthesis or inverse transform: Spherical harmonic synthesis evaluates a sum of spherical-harmonic basis functions up to degree n = N from complex coefficients.After exchanging sums, the longitude summation becomes a Fourier transform, leaving a Legendre-transform task.
- 2.3 Analysis or forward transform: SHT analysis computes spherical-harmonic coefficients from spatial values using a Fourier transform in longitude and a remaining Legendre transform.The integral over longitude is obtained with the Fourier transform before the Legendre step.
- 2.3 Analysis or forward transform: Gauss-Legendre quadrature replaces the latitudinal integral with a weighted sum over Gauss nodes and weights.Because the relevant integrand is a polynomial of degree at most 2N, the quadrature is exact for Nθ ≥ N + 1.
3 Optimization of the Gauss-Legendre algorithm
SHTns combines standard Gauss-Legendre optimizations with explicitly vectorized on-the-fly Legendre computations, reducing memory use while improving performance. The resulting implementation scales effectively across processors and retains high accuracy at large truncation degrees.
- Standard optimizations: SHTns applies FFTs, Hermitian symmetry, equatorial mirror symmetry, and polar localization to reduce transform work.These optimizations improve speed and can halve the operation count for forward and inverse transforms.
- On-the-fly algorithms and vectorization: On-the-fly Legendre computation avoids storing O(N^3) values, requiring about 8(N + 1)^2 bytes instead of at least 2(N + 1)^3 bytes.At N = 1023, the stated storage is about 8Mb rather than 2Gb.
- On-the-fly algorithms and vectorization: Explicit vectorization makes on-the-fly computation fastest with 4-value AVX vectors, while precomputed matrices are fastest with 2-value vectors for N ≲200.Vectorized recursion computes multiple x values simultaneously and can outperform loading precomputed values from memory.
- On-the-fly algorithms and vectorization: Runtime tuning selects the fastest among precomputed and on-the-fly variants, with observed performance gains between 10 and 30%.The startup overhead can be skipped for applications performing only a few transforms.
- Multi-threaded transform: Parallel speedup is good for N = 255 or above, excellent up to 8 threads for N ≥511, and useful up to 16 threads for N ≥2047.For N below 128, the paper reports no benefit from parallel transforms.
- Performance comparisons: SHTns is 2.8 times faster than Driscoll-Healy at N = 1023 and between two and ten times faster than the best alternative overall.Its efficiency is close to 1 for N ≥511, corresponding to nearly one computed element per clock cycle.
4 Conclusion and perspectives
SHTns combines vectorized on-the-fly computation with flexible transform routines to deliver high performance while reducing memory requirements. These capabilities support demanding spherical-geometry simulations, including prospective geodynamo runs at N=1023.
- SHTns uses explicitly vectorized on-the-fly algorithms that perform better than publicly available implementations, including at N=1023.The approach also reduces sensitivity to memory-bandwidth limitations in multi-threaded transforms.
- The library provides versatile truncation, normalization conventions, and scalar and vector routines for C/C++, Fortran, and Python.These interfaces target high-performance computing involving partial differential equations in spherical geometry.
- N=1023 spectral geodynamo simulations are presented as feasible in the next few years because of the performance gain and lower memory requirements.The paper connects this resolution with simulations closer to the dynamics of Earth’s core.