Summation-by-parts (SBP) finite difference methods are widely used in scientific applications alongside a special treatment of boundary conditions through the simultaneous-approximate-term (SAT) technique which enables the valuable proof of numerical stability. Our work is motivated by multi-scale earthquake cycle simulations described by partial differential equations (PDEs) whose discretizations lead to huge systems of equations and often rely on iterative schemes and parallel implementations to make the numerical solutions tractable. In this study, we consider 2D, variable coefficient elliptic PDEs in complex geometries discretized with the SBP-SAT method. The multigrid method is a well-known, efficient solver or preconditioner for traditional numerical discretizations, but they have not been well-developed for SBP-SAT methods on HPC platforms. We propose a custom geometric-multigrid preconditioned conjugate-gradient (MGCG) method that applies SBP-preserving interpolations. We then present novel, matrix-free GPU kernels designed specifically for SBP operators whose differences from traditional methods make this task nontrivial but that perform 3 × faster than SpMV while requiring only a fraction of memory. The matrix-free GPU implementation of our MGCG method performs 5 × faster than the SpMV counterpart for the largest problems considered (67 million degrees of freedom). When compared to off-the-shelf solvers in the state-of-the-art libraries PETSc and AmgX, our implementation achieves superior performance in both iterations and overall runtime. The method presented in this work offers an attractive solver for simulations using the SBP-SAT method.

Numerical modeling of earthquake dynamics and derived insight for seismic hazard relies on credible, reproducible model results. The SEAS (Sequences of Earthquakes and Aseismic Slip) initiative has set out to facilitate community code comparisons, and verify and advance the next generation of physics-based earthquake models that reproduce all phases of the seismic cycle. With the goal of advancing SEAS models to robustly incorporate physical and geometrical complexities, here we present code comparison results from two new benchmark problems: BP1-FD considers full elastodynamic effects and BP3-QD considers dipping fault geometries. Eight modeling groups participated in each benchmark, allowing us to explore these physical ingredients across multiple codes and better understand associated numerical considerations. We find that numerical resolution and computational domain size are critical parameters to obtain matching results, with increasing domain-size requirements posing challenges for volume-based codes even in 2D settings. Codes for BP1-FD implemented different criteria for switching between quasi-static and dynamic solvers, which require tuning to obtain matching results. In BP3-QD, proper remote boundaries conditions consistent with specified rigid body translation are required to obtain matching surface displacements. With these numerical and mathematical issues resolved, we obtain good agreement among codes in long-term fault behavior, earthquake recurrence intervals, and rupture features of peak slip rates and stress drops for both benchmarks. Including full inertial effects generates events with larger slip rates and rupture speeds compared to the quasi-dynamic counterpart. For BP3-QD, both dip angle and sense of motion (thrust versus normal faulting) alter ground motion on the hanging and foot walls, and influence event patterns, with some sequences exhibiting similar-sized characteristic earthquakes, and others exhibiting several earthquakes of differing magnitudes. These findings underscore the importance of considering full dynamics and non-vertical dip angles in SEAS models, as both influence short and long-term earthquake behavior, and associated hazards.

24 Numerical modeling of earthquake dynamics and derived insight for seismic hazard relies on credi25 ble, reproducible model results. The SEAS (Sequences of Earthquakes and Aseismic Slip) initiative 26 has set out to facilitate community code comparisons, and verify and advance the next gener27 ation of physics-based earthquake models that reproduce all phases of the seismic cycle. With 28 the goal of advancing SEAS models to robustly incorporate physical and geometrical complexities, 29 here we present code comparison results from two new benchmark problems: BP1-FD considers 30 full elastodynamic effects and BP3-QD considers dipping fault geometries. Eight modeling groups 31 participated in each benchmark, allowing us to explore these physical ingredients across multiple 32 codes and better understand associated numerical considerations. We find that numerical resolution 33 and computational domain size are critical parameters to obtain matching results, with increasing 34 domain-size requirements posing challenges for volume-based codes even in 2D settings. Codes 35 for BP1-FD implemented different criteria for switching between quasi-static and dynamic solvers, 36 which require tuning to obtain matching results. In BP3-QD, proper remote boundaries conditions 37 consistent with specified rigid body translation are required to obtain matching surface displace38 ments. With these numerical and mathematical issues resolved, we obtain good agreement among 39 codes in long-term fault behavior, earthquake recurrence intervals, and rupture features of peak 40 slip rates and stress drops for both benchmarks. Including full inertial effects generates events 41 with larger slip rates and rupture speeds compared to the quasi-dynamic counterpart. For BP342 QD, both dip angle and sense of motion (thrust versus normal faulting) alter ground motion on 43 the hanging and foot walls, and influence event patterns, with some sequences exhibiting similar44 sized characteristic earthquakes, and others exhibiting several earthquakes of differing magnitudes. 45 These findings underscore the importance of considering full dynamics and non-vertical dip angles 46 in SEAS models, as both influence short and long-term earthquake behavior, and associated hazards. 47 48