Scalable parallel 3-D TEM inversion via rational approximation of the matrix exponential

TL;DR

This paper introduces a scalable parallel 3-D electromagnetic inversion method using rational approximation of the matrix exponential, enabling efficient large-scale transient simulations on shared-memory systems.

Contribution

It presents a novel parallel implementation combining rational approximation with direct solvers, optimized for shared-memory architectures, and released as open-source in Julia.

Findings

Successfully recovered a synthetic 3-D conductivity structure with 700,000 degrees of freedom.

Demonstrated efficient parallel execution and independence of response computation from observation times.

Discussed computational bottlenecks and potential improvements for high-performance computing.

Abstract

We present a novel parallel implementation for large-scale three-dimensional electromagnetic inversion based on a Gauss-Newton framework combined with a rational near-best approximation of the matrix exponential for transient simulations. The method employs parallel direct solvers for the shifted linear systems arising from the partial fraction representation of the rational approximation and demonstrates efficient parallel execution on a shared-memory architecture using MPI. A key property of the approach is that the time dependence is entirely contained in the residuals of the employed rational functions, such that the computation of forward responses and sensitivities becomes effectively independent of the number of desired observation times. Model regularization is done with smoothness constraints, formulated with Raviart-Thomas elements. The linearized inverse problems are solved using LSQR, using an implicit parallel Jacobian operator. Numerical experiments demonstrate the successful recovery of a synthetic 3-D conductivity structure with approximately 700,000 degrees of freedom. The study further discusses computational bottlenecks related to memory consumption and shared-memory scalability arising from the simultaneous storage of multiple sparse matrix factorizations. Possible improvements based on preconditioned iterative solvers and distributed high-performance computing architectures are outlined. The implementation in the Julia programming language is released as open-source software to support reproducible research and further development by the geophysical inversion community.