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

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.