An efficient geometric method for incompressible hydrodynamics on the sphere

TL;DR

This paper introduces a scalable geometric numerical method for simulating incompressible fluid flow on the sphere, combining efficient matrix computations with structure-preserving integrators for long-term accuracy.

Contribution

The authors develop a novel, highly efficient geometric algorithm for 2D ideal fluid dynamics on the sphere, leveraging tridiagonal splitting and scalable parallel computing.

Findings

Achieved linear scaling up to 2500 cores for high-resolution simulations.

Enabled long-time, structure-preserving simulations of Euler's equations.

Demonstrated the method's effectiveness at resolution N=2048.

Abstract

We present an efficient and highly scalable geometric 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 $\mathcal{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$.