Backend-agnostic Julia framework for 3D modeling and inversion of gravity data

This paper presents a high-performance framework for three-dimensional gravity modeling and inversion implemented in Julia, addressing key challenges in geophysical modeling such as computational complexity, ill-posedness, and the non-uniqueness inherent to gravity inversion. The framework adopts a data-space inversion formulation to reduce the dimensionality of the problem, leading to significantly lower memory requirements and improved computational efficiency while maintaining inversion accuracy. Forward modeling and inversion operators are implemented within a backend-agnostic kernel abstraction, enabling execution on both multicore CPUs and GPU accelerators from a single code base. Performance analyses conducted on NVIDIA CUDA GPUs demonstrate substantial reductions in runtime relative to CPU execution, particularly for large-scale datasets involving up to approximately 3.3 million rectangular prisms, highlighting the scalability of the proposed approach. The inversion incorporates implicit model constraints through the data-space formulation and depth-weighted sensitivity, which mitigate depth-related amplitude decay and yield geologically coherent, high-resolution subsurface density models. Validation using synthetic models confirms the ability of the framework to accurately reconstruct complex subsurface structures such as vertical and dipping dykes. Application to field gravity data further demonstrates the robustness and practical utility of the GPU-accelerated framework, with the recovered models showing strong consistency with independent geological constraints and prior interpretations. Overall, this work underscores the potential of GPU-enabled computing in Julia to transform large-scale gravity inversion workflows, providing an efficient, extensible, and accurate computational solution for high-resolution geophysical studies.

💡 Research Summary

The paper introduces a high‑performance, backend‑agnostic framework for three‑dimensional gravity forward modeling and inversion, implemented entirely in the Julia programming language. The authors address two fundamental challenges in gravity inversion: the explosive growth of model parameters (M) relative to the number of observations (N) and the ill‑posed, non‑unique nature of the inverse problem. To mitigate these issues, they adopt a data‑space inversion formulation that operates in the observation space rather than the full model space. By applying a matrix identity, the original normal‑equation system (CᵀW_d⁻¹C + W_m) m = CᵀW_d⁻¹d is transformed into a reduced system of size N × N: δm = W_m Cᵀ(C W_m Cᵀ + W_d)⁻¹ δd. This reduction dramatically lowers memory consumption and computational cost because N (the number of gravity observations) is typically orders of magnitude smaller than M (the number of prisms). The framework also incorporates depth‑weighted sensitivity and implicit model constraints to counteract the depth‑decay of gravity amplitudes, thereby improving the recovery of deeper structures without explicit regularization beyond a standard L₂ term.

A central technical contribution is the use of KernelAbstractions.jl, a Julia package that provides a hardware‑agnostic abstraction for writing GPU‑style kernels. With a single kernel definition, the same code can be dispatched to a multicore CPU, an NVIDIA CUDA GPU, an AMD GPU, or Apple’s Metal backend simply by changing the runtime backend selection. This eliminates the need for separate CPU and device‑specific source files, reduces development effort, and ensures algorithmic consistency across platforms. The authors detail memory‑hierarchy optimizations: the sensitivity matrix resides in global memory, cell parameters are cached in shared memory, physical constants are placed in constant memory, and intermediate results are kept in registers. Such layout maximizes bandwidth utilization and minimizes latency, which is crucial for the large matrix‑vector products and factorizations that dominate gravity inversion.

Performance benchmarks focus on the construction of the sensitivity matrix (the most expensive step). For small models (< 10³ prisms) the CPU is marginally faster due to GPU kernel launch and host‑device transfer overheads. However, once the model size exceeds roughly 10⁴ prisms, the GPU’s massive parallelism yields a sub‑linear increase in runtime, while CPU time grows rapidly. For the largest test case (≈ 3.3 × 10⁶ rectangular prisms), the CPU requires over 1 000 seconds, whereas the CUDA implementation completes in 20–30 seconds, representing a speed‑up of one to two orders of magnitude. The experiment also reveals a memory ceiling: models larger than 3.3 × 10⁶ prisms could not be processed on the tested GPU (Mahti system) due to insufficient device memory, suggesting future work on multi‑GPU distribution or out‑of‑core algorithms.

Synthetic experiments validate the inversion quality. In the first test, two isolated high‑density blocks with different lateral positions and depths are recovered with accurate geometry; minor smoothing near edges is observed, which is expected given the L₂ regularization and the inherent non‑uniqueness of gravity data. The second synthetic case involves a dipping dyke, a more complex geological feature. The data‑space inversion successfully reconstructs the dip angle and thickness, demonstrating robustness to structural complexity.

A field‑data case study further confirms practical applicability. Real gravity measurements are inverted using the same framework, and the resulting density model aligns closely with independent geological interpretations and prior studies. The GPU‑accelerated workflow enables rapid iteration, making the approach suitable for operational exploration settings where timely results are essential.

In summary, the paper delivers three major advances: (1) a mathematically rigorous data‑space inversion that reduces problem dimensionality and memory demand; (2) a Julia‑based, backend‑agnostic GPU acceleration strategy that unifies CPU and GPU implementations under a single code base; and (3) a demonstration of scalable performance on realistic, large‑scale gravity problems, achieving near‑real‑time inversion for models with millions of cells. The authors suggest future extensions such as incorporating sparsity‑promoting (L₁) regularization, joint inversion with seismic or magnetic data, and scaling to multi‑GPU or distributed environments. This work showcases how modern high‑level languages combined with portable GPU abstractions can transform computational geophysics workflows.