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)?

The paper describes potential speed gains from combining FEniCS with a CUDA enabled solver backend. Such comparisons are always helpful to judge if it is worth the effort of combining such large frameworks. 

A few points from my side:

- We are dealing with memory bandwidth limited processing, please add the theoretical memory bandwidth of the CPU system and the GPUs in use. Scaling should be somewhere in that range.

- The color descriptions in the sticky air figure are incorrect

- Eq 4 only advects T not C_i

- For the sticky air T has to be 0 and the scenario cannot be mixed with thermal convection. Eq 4 and eq 1 cannot be used in a general sense. Please state that.

- Inventing little own benchmarks (3D Taylor) is nice but since you did not compare it to anything else their only point is to show the speedups again, but they all solve the same eqns more or less, so what do we gain from that? If you want to establish it as a benchmark, you need to show the solutions of other available codes.