An efficient geometric method for incompressible hydrodynamics on the sphere

We present an efficient and highly scalable geometric numerical method for two-dimensional ideal fluid dynamics on the sphere. The starting point is Zeitlin's finite-dimensional model of hydrodynamics. The efficiency stems from exploiting a tridiagonal splitting of the discrete spherical Laplacian combined with highly optimized, scalable numerical algorithms. For time-stepping, we adopt a recently developed isospectral integrator able to preserve the geometric structure of Euler's equations, in particular conservation of the Casimir functions. To overcome previous computational bottlenecks, we formulate the matrix Lie algebra basis through a sequence of tridiagonal eigenvalue problems, efficiently solved by well-established linear algebra libraries. The same tridiagonal splitting allows for computation of the stream matrix, involving the inverse Laplacian, for which we design an efficient parallel implementation on distributed memory systems. The resulting overall computational complexity is O(N^3) per time-step for N^2 spatial degrees of freedom. The dominating computational cost is matrix-matrix multiplication, carried out via the parallel library ScaLAPACK. Scaling tests show approximately linear scaling up to around 2500 cores for the matrix size N=4096 with a computational time per time-step of about 0.55 seconds. These results allow for long-time simulations and the gathering of statistical quantities while simultaneously conserving the Casimir functions. We illustrate the developed algorithm for Euler's equations at the resolution N=2048.