Efficient parallel finite-element methods for planetary gravitation: DtN and multipole expansions

This paper evaluates and implements three strategies—naive domain truncation, Dirichlet-to-Neumann (DtN) maps, and multipole expansions—for handling unbounded domains in parallel finite-element simulations of planetary gravitation using the MFEM package, demonstrating that while truncated meshes with coarsening are viable, DtN and multipole methods offer superior accuracy and cost-efficiency for large-scale geophysical modeling.

The Gist

Imagine you are trying to calculate the gravitational pull of a planet, like Earth or a moon. In the real world, gravity doesn't just stop; it stretches out into the infinite void of space forever. However, when scientists use computers to solve these problems, they hit a wall: computers can't handle "infinity." They need a box with a defined size to do their math.

This paper is about how to build the best possible "box" to simulate a planet's gravity without getting the answer wrong. The authors, working at the University of Cambridge, tested three different ways to handle the "empty space" outside the planet in their computer models.

Here is a breakdown of their findings using simple analogies:

The Problem: The Infinite Ocean

Think of the planet as a lighthouse sitting in the middle of an infinite ocean. You want to calculate the waves (gravity) hitting the lighthouse.

• The Computer's Limit: Your computer is a small boat. It can only measure the water in a small, finite circle around the lighthouse.
• The Challenge: If you just draw a circle around the lighthouse and say, "The water stops here and is flat," you get a wrong answer. The waves don't just stop; they ripple out forever. If your circle is too small, the ripples hit the edge of your boat and bounce back, messing up your calculation.

The Three Strategies Tested

The authors compared three ways to fix this "edge of the boat" problem:

1. The "Big Boat" Method (Naive Domain Truncation)

The Idea: Just make your boat (the computer mesh) huge. If you make the circle around the lighthouse big enough, the waves will be so tiny by the time they hit the edge that you can pretend they are zero.
The Result: This works, but it's inefficient. To get a super-accurate answer, you have to make the boat massive. This means the computer has to do a lot of extra math on empty water just to get the edges right. It's like buying a giant yacht just to measure the waves near the shore.
Verdict: It works, but it's expensive and slow.

2. The "Magic Mirror" Method (Dirichlet-to-Neumann or DtN)

The Idea: Instead of making the boat huge, you put a special "magic mirror" on the edge of your small circle. This mirror knows exactly how the waves should behave if the ocean were infinite. It tells the computer, "Don't stop the wave here; imagine it continues perfectly into infinity."
How it works: The mirror uses a mathematical recipe (spherical harmonics) to predict the future of the waves.
The Result: This is a winner. You can keep your boat small, but because the mirror is so smart, the computer thinks the ocean is infinite. It gives you a highly accurate answer with much less computing power.
The Catch: The mirror needs to "talk" to all the processors in the computer at once to share its prediction. In a supercomputer with thousands of processors, this communication can get a bit noisy, but the authors found a clever way to keep it efficient.

3. The "Ghost Sum" Method (Multipole Expansion)

The Idea: Instead of looking at the waves on the edge, you calculate the "ghost" of the planet's mass. You sum up all the tiny bits of the planet's weight and turn them into a list of numbers (multipole moments) that describe how the gravity looks from far away. You then use this list to set the conditions on the edge of your small circle.
The Result: Like the Magic Mirror, this is incredibly accurate and fast. It's great for static problems (where the planet isn't moving). However, if the planet is wobbling or deforming (like during an earthquake or ice melting), the "ghost sum" has to be recalculated every time, which adds a little extra work.

The Big Winner: The Magic Mirror (DtN)

The paper concludes that while the "Big Boat" method is okay, the Dirichlet-to-Neumann (DtN) method is the champion for modern, large-scale simulations.

• Why? It offers the best balance of speed and accuracy.
• The Analogy: Imagine you are painting a picture of a landscape.
  • Method 1 paints the whole infinite horizon, which takes forever.
  • Method 2 (DtN) paints a small frame but uses a special lens that makes the edges look like they stretch to infinity. It's fast, cheap, and looks perfect.

Why Does This Matter?

This isn't just about math; it's about understanding our planet and others.

• Glacial Rebound: When massive ice sheets melt, the land underneath bounces back up. This changes the planet's shape and gravity. To predict sea-level rise accurately, we need to model this "bouncing" with extreme precision.
• Planetary Science: Whether it's the Moon, Mars, or the weirdly shaped moon Phobos, these methods allow scientists to calculate gravity on irregular shapes without needing supercomputers to run for weeks.

The Takeaway

The authors have successfully built a "smart edge" for their computer models. They proved that you don't need a massive computer to simulate the infinite universe; you just need a smart way to tell the computer how the universe behaves at the edge of its view. This makes future simulations of Earth's climate, earthquakes, and planetary dynamics faster, cheaper, and more accurate.

Technical Summary

1. Problem Statement

The paper addresses a fundamental challenge in geophysical and planetary modeling: solving the Poisson equation for gravitational potential ($\nabla^2 \phi = 4\pi G \rho$) on an unbounded domain ($\mathbb{R}^3$) using Finite Element Methods (FEM).

• The Conflict: FEM requires bounded computational meshes, whereas gravitational fields extend to infinity.
• The Difficulty: Elliptic equations like Poisson's equation have solutions that decay slowly ($1/r$). Truncating the domain introduces significant errors unless the boundary is placed extremely far away, which drastically increases computational cost (degrees of freedom).
• Context: This is critical for Glacial Isostatic Adjustment (GIA) modeling and planetary dynamics, where self-gravitation must be coupled with viscoelastic deformation. While spectral methods work well for spherically symmetric models, they struggle with the lateral heterogeneity found in realistic Earth models, necessitating FEM.

2. Methodology

The authors implement and compare three strategies for handling the infinite exterior domain within the open-source MFEM library (with notes on FEniCS and Firedrake):

A. Naive Domain Truncation

• Approach: The computational domain is extended to a large finite sphere ($B$) with homogeneous Dirichlet ($\phi=0$) or Neumann ($\partial_n \phi=0$) boundary conditions.
• Implementation: Standard FEM discretization.
• Optimization: The mesh is coarsened in the exterior "buffer" layer to prevent an explosion in degrees of freedom as the domain size increases.

B. Dirichlet-to-Neumann (DtN) Map

• Theory: On a spherical boundary $\partial B$, the exterior potential is harmonic and can be represented by a spherical harmonic expansion. The DtN map provides an exact linear relationship between the potential ($\phi$) and its normal derivative ($\partial_n \phi$) on the boundary: $\partial_n \phi = -\sum \frac{\ell+1}{b} \phi_{\ell m} Y_{\ell m}$.
• Implementation:
  • A custom mfem::Operator class is created to handle the non-local boundary term.
  • Parallelization: Uses a split MPI communicator restricted to processors containing the boundary. It employs MPI_Allreduce to sum spherical harmonic coefficients, minimizing communication overhead.
  • Linear System: The DtN term is added as a low-rank perturbation to the stiffness matrix. A preconditioner is built using only the standard diffusion term ($A$) to maintain efficiency.

C. Multipole Expansion

• Theory: The exterior potential is expressed as a sum of multipole moments generated by the interior density distribution. The normal derivative on the boundary is calculated via volume integrals over the interior domain $M$.
• Implementation:
  • Similar to DtN, it uses a custom operator.
  • Static vs. Linearized: For static problems, it modifies the right-hand side (RHS) of the linear system. For linearized problems (coupling with displacement fields), it modifies the RHS based on the displacement vector $\mathbf{u}$.
  • Parallelization: Requires global reduction of multipole coefficients from all processors (interior and boundary), but since it only affects the RHS assembly, the iterative solver time remains comparable to standard Neumann problems.

3. Key Contributions

1. Parallel Implementation in Open-Source Codes: The paper provides the first robust, parallel implementations of DtN and Multipole methods within modern FEM libraries (specifically MFEM), overcoming the non-local communication challenges inherent in spherical harmonic expansions.
2. Efficient Linear Algebra: The authors demonstrate how to implement these non-local operators using matrix-free techniques and low-rank updates, avoiding the assembly of dense matrices that would destroy parallel scalability.
3. Extension to Linearized Problems: The methods are successfully extended to the linearized Poisson equation ($\nabla^2 \phi = -4\pi G \nabla \cdot (\rho \mathbf{u})$), which is essential for coupling self-gravitation with elastic/viscoelastic deformation in GIA modeling.
4. Benchmarking: Rigorous comparison of the three methods using an "offset-sphere" benchmark problem with a known analytical solution.

4. Results

• Accuracy vs. Cost:
  • Domain Truncation: Can achieve high accuracy ($\epsilon \approx 10^{-6}$) but requires a domain 50 times larger than the planet ($b/a=50$), introducing significant mesh complexity.
  • DtN and Multipole: Achieve superior accuracy ($\epsilon \approx 10^{-9}$) with a much smaller domain ($b/a \approx 1.4$) by increasing the spherical harmonic truncation order ($\ell_{max}$).
• Computational Efficiency:
  • The DtN and Multipole methods are 2–3 times faster than the large-domain truncation method for equivalent accuracy.
  • Solver Time: The DtN method adds negligible solver time because the DtN term acts as a low-rank perturbation; the number of conjugate gradient iterations remains similar to standard Neumann problems.
  • Assembly Time: Assembly time increases with $\ell_{max}$ but remains a small fraction of the total solver time for realistic applications.
• Parallel Scaling:
  • Both methods exhibit favorable weak scaling.
  • Strong scaling tests show that communication overhead (via MPI_Allreduce for spherical coefficients) does not become a bottleneck for up to hundreds of processors, provided $\ell_{max}$ is not excessively large.
• Real-World Applications:
  • PREM Earth Model: The DtN method reproduced the radial gravity profile of the Preliminary Reference Earth Model (PREM) with high fidelity compared to spectral methods.
  • Phobos (Moon): Successfully modeled the gravitational potential of the irregular, asymmetric shape of Phobos, demonstrating the method's robustness for non-spherical geometries.

5. Significance

• Enabling High-Fidelity GIA Modeling: By providing an efficient way to handle self-gravitation in laterally heterogeneous Earth models, this work removes a major bottleneck in Glacial Isostatic Adjustment simulations.
• Scalability: The methods are proven to scale well on parallel architectures, making them suitable for the next generation of high-resolution planetary dynamics simulations.
• Accessibility: By implementing these in MFEM (and providing code for FEniCS/Firedrake), the authors make advanced boundary treatment techniques accessible to the broader geophysical community, moving beyond specialized spectral codes.
• Validation: The study confirms that while naive truncation is possible with mesh coarsening, DtN and Multipole expansions are the superior choice for balancing accuracy and computational cost in large-scale parallel geophysical simulations.