| 20 May 2026
High-Performance Geodynamics on GPUs Using the PETSc CUDA Backend (GAUZZ v1.0.0)
Abstract. The GPU is a computing architecture designed for parallelism in vector and matrix operations. To solve the Stokes equations in geodynamics, we apply a preconditioned conjugate-gradient (PCG) method to the pressure Schur complement system. The algorithm is dominated by sparse matrix-vector products, vector inner products, vector updates, and repeated velocity subproblem solves, all of which map onto GPU architectures. We implement this solver on GPUs by coupling CUDA-enabled PETSc with FEniCS. PETSc AIJCUSPARSE matrices with CUDA vectors keep operators and field vectors resident in VRAM, which limits CPU-GPU transfers during iterations. The P1 mass matrix, the gradient operator, and the divergence operator are reused in the pressure-correction L2-projection step. HYPRE BoomerAMG was adopted for preconditioning the velocity-only subproblem. Accuracy and performance are evaluated on manufactured solutions, the SolCx benchmark, the sticky-air benchmark, 2D and 3D Rayleigh-Taylor instabilities, 3D thermal convection, and a 2D subduction model with nonlinear viscosity. We confirm that the GPU-based implementation reproduces CPU solutions. We distinguish between the execution time, which includes only the pressure correction step, and the wall time, which includes the workflow from mesh generation to Stokes solves. A single-GPU environment achieves a 5–11× reduction in execution time relative to a 32-core CPU. For nonlinear viscosity cases with per-step matrix updates, the wall time is reduced by 1.14–3.46× on a single GPU relative to a 32-core CPU. Using Multi-Process Service to coordinate two-GPU with a 16-core CPU reduces the wall time by 15 5.83× compared to a 32-core CPU.
The manuscript presents GPU-accelerated solvers for Stokes flow in geodynamic modelling. Matrices are assembled on the CPU, while PETSc’s GPU backend and HYPRE perform the solver operations on the device (PCG on the Schur complement with inner GMRES+HYPRE).
The study is interesting and novel, provides a detailed comparison across many examples, and is strengthened by its emphasis on reproducibility. Despite some major concerns, the study is worthy of publication in my opinion:
1. Too many and often repetitive results. The manuscript can be shortened and provide the same information and conclusions:
- exhaustive comparisons of 1, 2, 4, 8, 16, and 32 CPU cores in various sections doesn't make sense if the focus is GPU-accelerated computations. After the first section, please only show 32 cores.
- The authors spend a lot of figures and text on "verification" of CPU vs GPU results even though the same solver, the same matrices, and the same tolerances are used. In fact, it is actually the identical PETSc algorithm that runs (right?). Unless tolerances are not strict enough, the results should not depend on whether they are solved on the GPU (even with different AMG configuration). If they are, it is equally likely to get different results with different solver configurations for two different CPU configurations. Please shorten this!
- My understanding of the focus of this manuscript is CPU vs GPU performance. Maybe a selection of the examples presented is enough to establish the relative performance (especially if the performance behavior is very similar in different examples)
2. Solver and performance studies typically include strong scaling and weak scaling plots: while timings for different mesh sizes are reported, comparison with the optimal runtime (time/DoF) based on one data point is missing. GPU publications running on a single GPU typically report throughput like DoFs/s (as you cannot run on different numbers of cores). I think this would be relatively easy to add (most of the data is already collected) and would add a lot to the presentation.
3. You are implementing a Schur complement solver (PCG on Schur complement with inner GMRES) but I am surprised that you are not clearly establishing the robustness of the outer solver with mesh refinement (number of PCG iterations required until convergence). Can you add a table?
4. You claim a big disadvantage of block solvers is the fact that (F)GMRES does not have a short recurrence and requires restarts and a large number of temporary vectors. I think this is misleading:
- You can use methods with a short recurrence (BiCGStab, IDR(s), ...). In fact, Clevenger, Heister 2021 uses IDR(2) for improved performance for the large simulations.
- You end up using GMRES anyways for the inner solver (and the velocity problem is not much smaller than the whole system, especially for 3d Taylor-Hood). Why not use something else?
Minor comments:
- Assembly on CPU is not strictly required; this is a limitation of the software (see MFEM for example). Please make this clearer.
- Abstract: "The GPU is a computing architecture designed for parallelism in vector and matrix operations." - I find this sentence confusing. Are you referring to vectors/matrices used for FEM problems, are you referring to scalar and tensor cores in the GPU, or how parallelism in the GPU works (GPUs are much more generic than that)? Either way, please clarify and update the abstract.
- line 105: You report a Xeon Gold 5320 and a Xeon w5-2465X here. Which one is used for the CPU timings in the paper?
- line 165: I am not sure about this approach: your GMRES needs to solve very accurately for the outer CG to work reliably, right? A block preconditioner can get by with a spectrally equivalent operator (a single V-cycle, for example). Your GMRES also needs to store a potentially large basis and requires restarts. Why not use some Krylov method that does not require this?
- line 264 (figure 1):
The way the data is presented is somewhat unconventional (accumulating time over time). Some of this information might be better represented in a table instead of a figure. I don't see the value of showing CPU time for many different core numbers. The focus of the paper is GPU vs CPU comparison, right? A table with required number of iterations and time to solution would be helpful.
- line 275 (figure 2):
please remove <32 cores and rescale plots
- line 310 (figure 4):
The plot a) and its conclusion: Of course they become more expensive with tighter inner tolerance but you are ignoring that the required number of PCG iterations likely changes with the inner tolerance (or PCG doesn't converge if inner is too loose?). I don't understand the point of this figure. You should probably pick an outer tolerance for PCG, vary GMRES inner tolerance to find the best time to solution (and report that). All the scaling numbers between CPU (32 cores) and GPU are nearly identical for all data points, which is not surprising. Even if there are trends, your presentation makes it impossible to see. Why do you need to measure the cost for each PCG iteration (b)? Linear increase is expected and probably not worth showing here.
For c): is the number of required PCG iterations resolution independent (do you have an optimal method)? I find the comparison of GPU vs CPU (32 cores) for each problem size the most interesting here (a single or fixed number of PCG iterations is fine for that).
- figure 14:
What are these lines connecting to one corner of the bars? That looks wrong and confusing. Again, showing throughput might be more helpful here.
- Comparison on p30:
This section is very interesting. Could you provide more information about ASPECT? version, GMG or AMG used, cheap vs expensive iterations used. I am confused about S4: Why does your code not spend any time assembling matrices? Can you measure AMG setup separately (ASPECT should report that as well)?