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

TL;DR

The paper addresses solving the Poisson equation for planetary gravity on unbounded domains by comparing naive domain truncation, Dirichlet-to-Neumann (DtN) maps, and multipole expansions within MFEM. It develops and benchmarks parallel implementations for both the static gravitational field and the linearised perturbation driven by displacement, confirming that DtN and multipole methods achieve markedly higher accuracy at lower cost than large-domain truncation, with DtN offering particularly scalable parallel performance. The work provides detailed theory, numerics, and MFEM-based code for implementing these exterior-domain treatments, and validates them on realistic models such as PREM and Phobos. The results have practical impact for large-scale geophysical modelling, providing robust, scalable tools for self-gravitating simulations in planetary science and glaciation studies, with open-source availability for community use.

Abstract

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.

Figures (8)

• Figure 1: Representative meshes for the offset-sphere configuration. (a) Large-domain truncation. (b) Buffer-layer design for the DtN and multipole methods.
• Figure 2: Large-domain truncation with the zero Dirichlet boundary condition within the offset-sphere configuration. (a) Number of tetrahedral elements normalised by that at the reference case $b/a=10/7$, $N_{e,\mathrm{ref}}=498525$. (b) Relative $L_2$ error over $M$, $\varepsilon_{L^2(M)}$, versus $b/a$.
• Figure 3: DtN method for the offset-sphere configuration with ${b/a=10/7}$. (a) Signed relative error distributions over $M$ for six spherical harmonic orders ($\ell_{\max}=\{0,2,4,8,16,32\}$). (b) Assembly time (normalised by the time for assembling the diffusion term with the zero Dirichlet condition, $T_{ref}=3.81\mathrm{s}$) and solver time (normalised by that achieving $\varepsilon_{L^2(M)}=1.06\times 10^{-6}$ with a large domain of $b/a=50$, $T_{ref}=56.01\mathrm{s}$) versus $\ell_{\max}$. All times are averaged over 20 runs with 8 CPU processors. (c) Relative $L_2$ error over $M$, $\varepsilon_{L^2(M)}$, versus $\ell_{\max}$.
• Figure 4: Multipole expansion method for the offset-sphere configuration with ${a/b=10/7}$. (a) Signed relative error distributions over $M$ for six spherical harmonic orders ($\ell_{\max}=\{0,2,4,8,16,32\}$). (b) Assembly time (normalised by the time for assembling the diffusion term with the zero Dirichlet condition, $T_{ref}=3.81\mathrm{s}$) and solver time (normalised by that achieving $\varepsilon_{L^2(M)}=1.06\times 10^{-6}$ with a large domain of $b/a=50$, $T_{ref}=56.01\mathrm{s}$) versus $\ell_{\max}$. All times are averaged over 20 runs with 8 CPU processors. (c) Relative $L_2$ error over $M$, $\varepsilon_{L^2(M)}$, versus $\ell_{\max}$.
• Figure 5: Linearised gravity problem on the offset-sphere configuration with presumed deformation $\mathbf{u}=(1,1,1)b$. (a) Signed relative error distributions over the offset sphere $M$ for the multipole method at six spherical harmonic orders ($\ell_{\max}=\{0,2,4,8,16,32\}$). (b) Assembly time (dashed, normalised by that for assembling the diffusion term with the zero Dirichlet condition, $T_{ref}=3.81\mathrm{s}$) and solver time (normalised by that achieving $\varepsilon_{L^2(M)}=1.06\times 10^{-6}$ with a large domain of $b/a=50$, $T_{ref}=56.01\mathrm{s}$) versus $\ell_{\max}$ comparing the DtN (blue) and multipole (red) methods. (c) Relative $L_2$ error over $M$ versus $\ell_{\max}$ for both methods.