Efficient and scalable finite-element magnetotelluric modelling on high-order meshes

SUMMARY 3-D Forward modelling of magnetotelluric (MT) data remains a computationally challenging task, particularly when accurate broad-band MT responses are simulated for real-world problems that often involve complex multiscale bathymetry and/or topography. To overcome this challenge, we developed a new efficient numerical approach for 3-D MT forward modelling that combines high-order Nédélec-type finite elements and high-order meshes, allowing us to obtain superior accuracy and account for complex material boundaries and interfaces. Despite gains in accuracy, higher order FE solvers are often considered impractical owing to higher memory consumption and a more ill-conditioned system. To overcome these limitations, we use an iterative solver accelerated by the Low-Order-Refined (LOR) pre-conditioner, which uses spectrally equivalent low-order operators, rendering the complexity independent of the polynomial degree. Another key novelty is a matrix-free implementation, where the action of the high-order operator is computed efficiently without explicit matrix assembly. The low-order system is solved using an Auxiliary Space Maxwell (AMS) solver based on a multigrid solver. We demonstrate the efficiency in a series of numerical experiments. Scalability analysis on a 3-D benchmark model demonstrates that the LOR pre-conditioner significantly outperforms the current state-of-the-art AMS pre-conditioner in terms of CPU time and memory usage, especially for higher polynomial degrees. Excellent scalability is confirmed by solving a problem with up to $1.5\times 10^9$ degrees of freedom in less than 2 min using 16 384 CPU cores, which is, to our knowledge, the largest 3-D MT problem reported to date. We also illustrate that high-order hexahedral meshes allow for accurate discretization of complex geometries, such as topography, with substantially fewer elements than conventional linear meshes. Finally, the capability of the integrated approach is demonstrated on a real 3-D model crossing the ocean trench in the Aleutian subduction zone. The proposed methods pave the way for more efficient and accurate 3-D MT modelling that is crucial for the inversion of complex MT data sets.