Accelerating and scaling up SEAS simulations using GPUs in Julia

Sequences of Earthquakes and Aseismic Slip (SEAS) simulations focus on km scale models and consider all phases of faulting from aseismic slip to earthquake nucleation, propagation and termination. Many codes exist that use different approaches to tackle SEAS simulations; however, solutions designed to leverage the potential of GPUs to parallelize and speed up the simulation are limited (although recent examples are emerging such as PyQuake3D). In this work, we propose a GPU parallelized SEAS quasi-dynamic solver written in Julia adopting the Spectral Boundary Integral Method (SBIM). The SBIM approach is optimal for GPUs as it is more memory efficient in respect to other mesh-based solvers, thus enabling to efficiently run high resolution simulations with around 10 million nodes on the fault plane. We rearrange the rate-and-state equations to solve for the slip rate and adopt a slightly modified Newton-Raphson algorithm for root finding. We introduce elastic bulk by using an analytical stress-slip relationship in the Fourier/wavenumber domain. Most of the operations that are carried out in the solver are element wise and thus can be run in parallel on GPUs, significantly cutting down on computation time as the domain resolution increases. While FFT is inherently not fully parallelizable, GPU kernels are available to efficiently perform Fourier transforms on GPUs. Moreover, by using the FFT algorithm, the numerical complexity for calculating the stress is reduced from O(N²) to O(NlogN). To verify the correctness of our solver, we use the BP4-QD benchmark and show comparable results with other outputs hosted by SCEC. We then measure the runtime of the solver on CPU and 2 NVIDIA GPUs, the RTX4060 8GB and the A100 40GB, and show a x5 to x16 speedup for simulations depending on the GPU. Finally, we run the BP4-QD problem on the A100 GPU, decreasing the indicated node spacing and Dc values by an order of magnitude. This simulation yielded 36 events of Mw > 7 and 181 events of Mw between 5 and 6, showing emergence of complexity. Moreover, we observe that the earthquake’s nucleation points are distributed along the edges of the rate-weakening patch. The smaller events are mostly concentrated on the four corners and the two sides parallel to the slipping direction while the bigger events are distributed more uniformly all around the border.