🔢 mathematics

A Cubed Sphere Fast Multipole Method

This paper presents a new kernel-independent Fast Multipole Method for efficiently summing pairwise particle interactions on a sphere by utilizing barycentric Lagrange interpolation on a cubed sphere grid, demonstrating its effectiveness across various geophysical and tidal applications in both serial and parallel implementations.

Original authors: Anthony Chen, Robert Krasny

Published 2026-04-01

📖 5 min read🧠 Deep dive

CC BY 4.0

Original authors: Anthony Chen, Robert Krasny

Original paper licensed under CC BY 4.0 (http://creativecommons.org/licenses/by/4.0/). ✨ This is an AI-generated explanation of the paper below. It is not written or endorsed by the authors. For technical accuracy, refer to the original paper. Read full disclaimer

Imagine you are standing on a giant, spinning beach ball (the Earth) covered in millions of tiny, glowing fireflies. Each firefly represents a piece of data—maybe the temperature of the ocean, the speed of the wind, or the height of the tide.

Now, imagine you want to know the "total influence" of every single firefly on every other firefly. For example, how much does the heat from a firefly in Tokyo affect the temperature felt by a firefly in New York?

In the world of computer science, calculating this for millions of fireflies is a nightmare. If you have 1 million fireflies, a standard computer has to check every single pair against each other. That's 1 trillion comparisons. It's like trying to introduce every person at a stadium of 1 million people to every other person individually. It would take forever, and your computer would likely overheat and give up.

This paper introduces a clever new trick called the Cubed Sphere Fast Multipole Method (CSFMM) to solve this problem quickly and accurately.

Here is how it works, broken down into simple concepts:

1. The Problem: The "Too Many People" Party

The authors are dealing with complex math problems on a sphere (like predicting weather or ocean tides). These problems require summing up interactions between millions of points. Doing this directly is too slow ( $O(N^2)$ complexity). It's like trying to count every grain of sand on a beach by picking them up one by one.

2. The Solution: The "Group Chat" Strategy

Instead of talking to every single firefly individually, the new method groups them into neighborhoods.

The Cubed Sphere: Imagine wrapping the beach ball in a cube. You can flatten the faces of the cube to make a grid. This creates a neat, organized map of the sphere without the messy distortions you get with traditional latitude/longitude lines (which bunch up at the poles like a crumpled umbrella).
The Tree: The computer builds a family tree of these neighborhoods. It starts with the whole sphere, then splits it into 6 big squares, then splits those into smaller squares, and so on, until the squares are tiny.

3. The Magic Trick: "Proxy" Representatives

This is the most important part. When a firefly in New York wants to know the influence of a huge group of fireflies in Tokyo, it doesn't need to talk to all of them.

The Old Way: The computer calculates the distance and influence from New York to every single Tokyo firefly.
The New Way (CSFMM): The computer picks a few "representatives" (called proxy particles) from the Tokyo group. It asks New York to talk only to these representatives.
The Secret Sauce: The authors use a mathematical tool called Barycentric Lagrange Interpolation. Think of this as a super-smart translator. It allows the computer to say, "If I know how the representatives feel, I can mathematically guess exactly how the whole group feels, without asking them individually."

If the groups are far apart, the computer uses these representatives. If the groups are right next to each other (like fireflies in the same neighborhood), it just calculates them directly because there aren't many of them.

4. Why This is a Big Deal

Speed: This method turns a task that takes days into one that takes minutes. For a massive simulation, it goes from being impossible to being routine.
Versatility: It works for almost any type of "influence" (gravity, heat, fluid flow) as long as the math happens on a sphere.
Accuracy: It doesn't just guess; it uses high-precision math to ensure the "representatives" tell the truth about the whole group.

Real-World Examples from the Paper

The authors tested this on four major problems:

Weather and Fluids: Simulating how wind swirls around the Earth (the Barotropic Vorticity Equation). They watched a "storm" move across the globe and saw that their fast method predicted its path just as well as the slow method, but much faster.
Ocean Tides: Calculating how the ocean's own weight pulls on the Earth (Self-Attraction and Loading). This is crucial for accurate tide predictions for ships and coastal cities.
Earth's Shape: Solving equations that describe how the Earth's surface bends under pressure (Biharmonic equation).

The Bottom Line

Imagine you are at a massive concert. You want to know the vibe of the whole crowd.

The Old Way: You walk up to every single person and ask, "How are you feeling?" (Takes 10 hours).
The New Way: You pick 5 people from the front row, 5 from the back, and 5 from the sides. You ask them, and based on their answers, you instantly know the vibe of the entire stadium with 99.9% accuracy. (Takes 10 seconds).

This paper gives scientists a new, super-fast "group chat" system for the Earth, allowing them to simulate weather, oceans, and gravity with a level of detail and speed that was previously out of reach.

1. Problem Statement

The paper addresses the computational challenge of evaluating integral transforms and convolutions on the sphere ( $S \subset \mathbb{R}^3$ ). These transforms arise in various geophysical and astrophysical applications, such as solving the Poisson and biharmonic equations, modeling fluid flow via the barotropic vorticity equation, and calculating self-attraction and loading (SAL) potentials for ocean tides.

The core mathematical problem is approximating the integral:
$\phi(x) = \int_S K(x, y)f(y)dS(y)$
where $K$ is a kernel (e.g., Green's function) and $f$ is an input data field. When discretized using $N$ particles (grid cells), this becomes an $N$ -body sum:
$\phi(x_i) \approx \sum_{j=1}^N K(x_i, x_j)w_j$
Direct summation of this equation requires $O(N^2)$ operations, which becomes prohibitively expensive for large $N$ (e.g., $N > 10^6$ ). Existing methods often rely on Spherical Harmonic Transforms (SHT), which face challenges such as the need for interpolation on non-lat/lon grids, spurious oscillations from truncation, and high computational costs for localized features.

2. Methodology

The authors propose a new Cubed Sphere Fast Multipole Method (CSFMM) and a corresponding Cubed Sphere Tree Code (CSTC). The methodology adapts the Fast Multipole Method (FMM) framework to the geometry of the sphere using the following key components:

Cubed Sphere Grid: The sphere is partitioned using a cubed sphere grid (projecting a cube's faces onto the sphere). This grid is used for two purposes:
1. Discretizing the integral transform.
2. Constructing a hierarchical quadtree of particle clusters (cells) to accelerate the $N$ -body sum.
Barycentric Lagrange Interpolation: Instead of using analytic multipole expansions (which require specific kernel forms) or evaluating kernels at off-sphere points (which can be undefined for spherical kernels), the method uses barycentric Lagrange interpolation on the cubed sphere cells.
- Kernel Independence: The scheme is "kernel-independent," requiring kernel evaluations only at points on the sphere. This is crucial for kernels like the approximate SAL kernel, which may be undefined for points where the dot product $x \cdot y > 1$ (off-sphere).
- Proxy Particles: Within each cell, Chebyshev points serve as "proxy particles." Interactions between distant clusters are approximated by interpolating the kernel between these proxy points.
Interaction Types: The method categorizes interactions into four types based on the separation of target and source clusters:
1. Particle-Particle (PP): Direct summation for nearby particles.
2. Particle-Cluster (PC): Target particle interacts with a distant source cluster via source interpolation.
3. Cluster-Particle (CP): Target cluster interacts with a distant source particle via target interpolation.
4. Cluster-Cluster (CC): Both target and source are distant clusters; interactions are approximated via interpolation in both variables.
Algorithm Flow:
- Upward Pass: Computes "proxy weights" for source clusters by aggregating particle data up the tree.
- Dual Tree Traversal: Uses a dual tree traversal (target and source trees) to identify well-separated clusters and apply PC, CP, or CC approximations based on a Multipole Acceptance Criterion (MAC).
- Downward Pass: Interpolates accumulated proxy potentials from parent clusters down to leaf particles to reconstruct the final field.

3. Key Contributions

Novel FMM for Spheres: Development of the first FMM specifically designed for the cubed sphere geometry using barycentric Lagrange interpolation, avoiding the need for analytic series expansions.
Kernel Independence: The method works for general kernels defined only on the sphere, overcoming limitations of previous methods (like BLDTT) that require off-sphere evaluations.
Efficiency: Achieves $O(N)$ complexity for the CSFMM and $O(N \log N)$ for the CSTC, a significant improvement over the $O(N^2)$ direct summation.
Application to Real-World Problems: Successfully applied the method to four distinct geophysical problems:
1. Poisson and Biharmonic equations.
2. Barotropic Vorticity Equation (BVE) on a rotating sphere.
3. Self-Attraction and Loading (SAL) potential for tidal calculations.

4. Results

The authors validated the method using the NCAR Derecho supercomputer (AMD EPYC processors) with up to $N \approx 2.6 \times 10^6$ particles.

Accuracy:
- The method exhibits 2nd order convergence with respect to angular grid spacing ( $\Delta \phi$ ).
- The approximation error introduced by the FMM is negligible compared to the discretization error of the midpoint rule when parameters (MAC and interpolation degree $n$ ) are tuned (e.g., MAC=0.7, $n=6$ ).
- Errors are uncorrelated with the exact solution, indicating the error is dominated by the grid discretization rather than the FMM approximation.
Performance (Serial):
- Speedup: For $N = 655,362$ , the CSFMM is approximately 30 times faster than the CSTC and ~180 times faster than direct summation for the Laplace kernel.
- Scalability: The runtime scales linearly ( $O(N)$ ) for CSFMM and $O(N \log N)$ for CSTC.
Performance (Parallel):
- The methods show near-linear speedup in parallel MPI implementations.
- Parallel efficiency drops at high core counts primarily due to the serial nature of the upward and downward passes (which consume a small fraction of total time).
Specific Applications:
- BVE: Simulated Rossby-Haurwitz waves and Gaussian vortices. The CSFMM preserved vorticity dynamics with negligible error compared to direct summation while being 36x faster.
- SAL: Computed the self-attraction potential for a variable-resolution ocean grid ( $N \approx 2.3$ million). The convolution/CSFMM result matched the spherical harmonics approach but was smoother and computationally efficient (90s serial runtime).

5. Significance

This work provides a robust, scalable, and kernel-independent framework for solving integral equations on the sphere, which is critical for modern geophysical fluid dynamics.

Overcoming SHT Limitations: It offers a viable alternative to Spherical Harmonic Transforms, particularly for problems requiring high resolution of localized features or operating on unstructured/variable-resolution grids (like the MPAS-Ocean model).
Enabling High-Fidelity Simulations: By reducing the computational cost of $N$ -body sums from $O(N^2)$ to $O(N)$ , the method enables the inclusion of complex forcing terms (like SAL) in ocean models at every time step, rather than just at prescribed intervals.
Future Impact: The authors plan to integrate this scheme into the Modular Ocean Model 6 (MOM6) and extend it to the Shallow Water Equations, potentially revolutionizing how global and regional tidal and weather models handle particle interactions.