Efficient and Scalable Finite-Element Magnetotelluric Modeling on High-Order Meshes

Summary Three-dimensional (3-D) forward modeling of magnetotelluric (MT) data remains a computationally challenging task, particularly when accurate broadband MT responses are simulated for real-world problems that often involve complex multi-scale bathymetry and/or topography. To overcome this challenge, we developed a new efficient numerical approach for 3-D MT forward modeling 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) preconditioner, 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 preconditioner significantly outperforms the current state-of-the-art AMS preconditioner 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 × 109 degrees of freedom in less than 2 minutes 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 modeling that is crucial for the inversion of complex MT data sets.