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

The Poisson equation governing a planet’s gravitational field is posed on the unbounded domain, $\mathbb{R}^3$, whereas finite-element computations require bounded meshes. We implement and compare three strategies for handling the infinite exterior in the finite-element method: (i) naive domain truncation; (ii) Dirichlet-to-Neumann (DtN) map on a truncated boundary; (iii) multipole expansion on a truncated boundary. While all these methods are known within the geophysical literature, we discuss their parallel implementations within modern open-source finite-element codes, focusing specifically on the widely-used MFEM package. We consider both calculating the gravitational potential for a static density structure and computing the linearised perturbation to the potential caused by a displacement field - a necessary step for coupling self-gravitation into planetary dynamics. In contrast to some earlier studies, we find that the domain truncation method can provide accurate solutions at an acceptable cost, with suitable coarsening of the mesh within the exterior domain. Nevertheless, the DtN and multipole methods provide superior accuracy at a lower cost within large-scale parallel geophysical simulations despite their need for non-local communication associated with spherical harmonic expansions. The DtN method, in particular, admits an efficient parallel implementation based on an MPI-communicator limited to processors that contain part of the mesh’s outer boundary. A series of further illustrative calculations are provided to show the potential of the DtN and multipole methods within realistic geophysical modelling.

💡 Research Summary

The paper addresses a fundamental challenge in planetary gravity modelling: solving Poisson’s equation for the gravitational potential on an unbounded domain using finite‑element methods (FEM), which inherently require a bounded mesh. Three strategies for handling the infinite exterior are implemented, benchmarked, and compared within the modern open‑source FEM library MFEM, with extensions to FEniCS and Firedrake.

1. Large‑domain truncation – The simplest approach expands the computational domain B far enough that homogeneous Dirichlet (φ=0) or Neumann (∂ₙφ=0) conditions can be imposed on its outer surface. The weak form reduces to a standard symmetric positive‑definite system that can be assembled and solved in parallel with PETSc or Hypre. Accuracy improves only slowly with the radius of B, but the authors show that aggressive mesh coarsening in the exterior can keep the number of degrees of freedom modest. Compatibility of the density field (∫ρ dx=0) is required for Neumann conditions; otherwise Dirichlet must be used.

2. Dirichlet‑to‑Neumann (DtN) map – For a spherical outer boundary, the exact relation between the potential and its normal derivative can be expressed in a spherical‑harmonic basis. The DtN operator is diagonal in this basis, and truncating the expansion at degree ℓmax controls the error. Implementing DtN in MFEM requires (i) computing spherical‑harmonic coefficients of the boundary solution, (ii) coupling these coefficients to the FEM degrees of freedom, and (iii) performing a reduced‑communicator MPI reduction that involves only the processes owning the outer boundary. This design limits global communication and yields excellent strong‑scaling up to thousands of cores. The same DtN formulation can be reused for the linearised perturbation problem that arises when coupling gravity to viscoelastic deformation.

3. Multipole expansion – This classic potential‑theory technique evaluates volume integrals of the interior density to obtain multipole moments Qℓm. The truncated series provides non‑homogeneous Dirichlet data on the spherical boundary. Implementation mirrors the DtN workflow: local contributions to Qℓm are accumulated on each rank, a global MPI_Reduce forms the full set of moments, and the boundary values are imposed on the FEM system. Although it requires slightly more global communication than DtN, the method attains comparable accuracy for the same ℓmax.

Benchmark results – Four test cases are examined: (a) a simple offset‑sphere with an analytical solution, (b) a non‑spherical prism‑like interior, (c) the Earth model PREM, and (d) a Phobos‑type moon model. For static gravity, ℓmax≈12 yields relative errors below 10⁻⁶ for both DtN and multipole methods, while the truncation method needs a domain radius at least five times larger to reach the same error. In the linearised displacement‑gravity problem, DtN is the cheapest (≈60 % of the runtime of the truncation method) and uses the least memory; the multipole method costs about 80 % of the truncation runtime. Strong‑scaling tests up to 1024 cores show DtN maintaining >90 % parallel efficiency because only the boundary ranks participate in the global spherical‑harmonic reductions. The truncation method, by contrast, suffers from increased communication as the whole mesh is distributed.

Conclusions and implications – The study demonstrates that while a naïve large‑domain approach is easy to implement and can be acceptable for modest problem sizes, DtN provides the best combination of accuracy, computational cost, and scalability for large‑scale planetary gravity simulations. The multipole expansion offers a viable alternative with similar accuracy and only modestly higher communication overhead. Both methods are readily integrated into MFEM and, with some additional bookkeeping, into higher‑level packages such as FEniCS and Firedrake. Their ability to handle the linearised perturbation problem makes them especially valuable for coupling self‑gravity to viscoelastic or thermo‑mechanical Earth‑system models, such as glacial isostatic adjustment, mantle convection, or tidal deformation studies. The paper therefore supplies the geophysical community with practical, high‑performance tools for realistic, laterally heterogeneous planetary gravity modelling.