GPU Power Boosts Mantle Flow Simulations Up to Elevenfold, New Solver Shows

Simulating the slow, relentless churn of Earth’s mantle has always been a computational brute-force problem. The solid rock beneath our feet flows over millions of years, driven by heat and density contrasts, and capturing that creeping motion numerically means solving the Stokes equations of incompressible viscous flow over and over again at every step of a simulation. Now a team of researchers from Kangwon National University and Utrecht University has shown that the graphics processing units best known for rendering video games can do this heavy lifting dramatically faster than traditional processors. Their new open-source code, GAUZZ v1.0.0, described in the journal Geoscientific Model Development, couples the finite-element platform FEniCS with the CUDA-enabled PETSc library to run geodynamic Stokes solvers almost entirely on a GPU, delivering throughput six to eleven times higher than a 32-core CPU while reproducing CPU results with the same accuracy.

The heart of the challenge lies in the structure of the equations themselves. Discretizing incompressible Stokes flow with finite elements produces a coupled saddle-point system linking velocity and pressure, and for a three-dimensional problem resolved on an N-by-N-by-N grid the number of unknowns grows as the cube of the resolution. Direct factorization methods, while robust, demand memory that scales even faster, making them impractical for large models. The team instead exploited the Schur complement structure of the system: eliminating the velocity unknown leaves a pressure equation that can be solved with a preconditioned conjugate-gradient iteration, where each outer iteration requires one inner solve of a symmetric, positive-definite velocity-only subproblem. That inner solve, handled by GMRES accelerated with the BoomerAMG algebraic multigrid preconditioner from the HYPRE library, dominates the total cost, and it is precisely the kind of workload that maps beautifully onto GPU hardware.

The reason is memory bandwidth. The dominant operations inside Krylov solvers, sparse matrix-vector products and vector updates, perform only about one multiplication and one addition per value read from memory, a low ratio of arithmetic to memory traffic known as low arithmetic intensity. GPUs are built for exactly this streaming pattern. The NVIDIA RTX 6000 Ada cards used in the study have a theoretical bandwidth of 960 gigabytes per second, compared with 153.6 gigabytes per second for the Intel Xeon Gold 6338 workstation serving as the CPU baseline. Measured with a vector-update benchmark, the GPU sustained 824 gigabytes per second, or 86 percent of peak, while the 32-core CPU managed 113.2 gigabytes per second, or 74 percent of peak. That sustained ratio of roughly 7.3 times turns out to be an excellent predictor of the speedups observed across the Stokes benchmarks.

Keeping data on the GPU was the central engineering challenge. The authors used PETSc’s AIJCUSPARSE sparse matrix type and CUDA vectors so that all operators and field vectors remain resident in the GPU’s 48 gigabytes of video memory, minimizing costly transfers between host and device. Because FEniCS 2019 assembles matrices on the CPU, repeatedly rebuilding operators during each conjugate-gradient cycle would have reintroduced transfer bottlenecks. The team’s solution was a pre-assembly strategy: the pressure mass matrix, the gradient operator, and the divergence operator are assembled once and stored in GPU memory for reuse throughout the iteration, requiring just 510 megabytes for 3.4 million degrees of freedom and 1.3 gigabytes for 8.9 million. The same strategy was applied to the CPU implementation to keep the comparison fair.

Intriguingly, the optimal multigrid settings differed between architectures. On the GPU, the PMIS and HMIS coarsening algorithms built coarse grids faster, and scaled Jacobi smoothing with Extended+i interpolation maximized parallel efficiency. On the CPU, the classic Falgout coarsening scheme with symmetric SOR smoothing gave the strongest per-iteration error reduction. In a manufactured-solution test with a viscosity contrast of one million, the GPU’s HMIS configuration needed 17 GMRES iterations in 0.413 seconds, while the CPU’s Falgout setup needed 13 iterations in 5.33 seconds. The GPU took four more iterations yet finished the accumulated work thirteen times faster, a striking illustration of how much cheaper each iteration becomes when thousands of arithmetic units stream data at high bandwidth.

The benchmarks spanned the standard repertoire of computational geodynamics. In the SolCx problem, which features a sharp viscosity jump of six orders of magnitude, the GPU completed five pressure-correction iterations at the highest resolution of 9.4 million degrees of freedom in 2.25 seconds against 17.82 seconds on the CPU, an approximately eightfold speedup, with velocity errors below one part in a billion in both cases. A 3D thermo-mechanical convection run on a 100-by-40-by-100 grid with the MINI element, totaling 8.9 million degrees of freedom, averaged 1.9 seconds per time step on the GPU versus 17.1 seconds on the CPU, and the root-mean-square velocity converged to 42.848, matching the canonical value of 42.865 from the classic Blankenbach benchmark. The Nusselt number of 4.881 likewise agreed with published references, confirming that acceleration did not come at the price of physical fidelity.

Accuracy held up in the harder tests as well. A 2D Rayleigh-Taylor instability benchmark, tracked with a level-set method for the compositional interface, ran at 512-by-512 resolution in 0.6 seconds per step on the GPU versus 6.3 seconds on the CPU, an elevenfold speedup, and the evolving root-mean-square velocity matched both the CPU run and literature reference values. A sticky-air free-surface benchmark at 10.1 million degrees of freedom produced topography evolution identical on GPU and CPU and consistent with the analytical solution. Most demanding of all was a 2D subduction model with nonlinear Mohr-Coulomb viscosity, solved with Picard iteration. As the slab descended, the stiffness matrix’s condition number grew, forcing more GMRES iterations on both platforms, yet the GPU remained roughly six to seven times faster because its per-iteration cost stayed so much lower.

A single GPU’s 48 gigabytes of memory, however, eventually becomes the binding constraint on resolution. The two-GPU configuration, using MPI-based domain decomposition coordinated through NVIDIA’s Multi-Process Service, extended the maximum problem size to 2.5 times ten to the seventh degrees of freedom, well beyond what one card could hold, and increased throughput by up to 2.1 times at sufficiently large problem sizes. Below a resolution threshold, communication and synchronization overhead ate the parallel gains, but in 3D, where degrees of freedom grow rapidly with refinement, the two-GPU advantage appeared consistently. Pairing a 16-core CPU with two GPUs through MPS delivered 5.19 times the end-to-end throughput of the 32-core CPU in the 3D convection benchmark, and even outpaced the community code ASPECT’s block-preconditioned solver by 3.74 times on the same hardware for that case.

The implications reach well beyond benchmark suites. Faster Stokes solves make high-resolution forward models of subduction zones, rifts, and mantle plumes more tractable, and they directly reduce the cost of adjoint inversions, in which each iteration demands both a forward simulation and a backward adjoint integration over the full assimilation window. The authors note that data-assimilation reconstructions of past mantle flow, such as those reconstructing North American subduction history, can cost an order of magnitude more than simplified models, and GPU acceleration could extend the assimilation windows and parameter-calibration loops that fit within a fixed runtime budget. Extrapolating their memory scaling, eight RTX 6000 Ada cards would support roughly 70 million degrees of freedom, and newer 96-gigabyte cards could push past 140 million. HPC-class GPUs with even greater bandwidth and memory offer a path toward global spherical-shell convection at the tens-of-kilometers grid spacing needed to resolve narrow shear zones and thermal boundary layers.

GAUZZ v1.0.0, released under the GNU General Public License with source code, benchmark input scripts, and a Dockerfile archived on Zenodo, represents a pragmatic middle ground in GPU-accelerated geodynamics. Unlike pseudo-transient approaches that avoid global matrices altogether, it brings the Krylov-plus-algebraic-multigrid machinery standard in finite-element community codes onto the GPU, systematically evaluating a combination that had received surprisingly little quantitative attention. Limitations remain, notably the CPU-side matrix assembly inherited from FEniCS and the incompatibility between FEniCS’s high-level layer and GPU-aware MPI, which the authors circumvented with native PETSc calls. But the message is clear: the same silicon that renders game worlds can now churn through the mathematics of Earth’s deep interior an order of magnitude faster, opening the door to planetary simulations that were previously out of reach.

Subject of Research: GPU-accelerated finite-element solution of the Stokes equations for geodynamic mantle flow simulations

Article Title: High-performance geodynamics on GPUs using the PETSc CUDA backend (GAUZZ v1.0.0)

Article References: Lee, K.-M., Jang, D.-K., Thieulot, C., & So, B.-D. (2026). High-performance geodynamics on GPUs using the PETSc CUDA backend (GAUZZ v1.0.0). Geoscientific Model Development, 19(19), 9781-9805. https://doi.org/10.5194/gmd-19-9781-2026

Image Credits: AI Generated

DOI: 10.5194/gmd-19-9781-2026

Keywords: geodynamics, GPU computing, Stokes equations, PETSc, FEniCS, algebraic multigrid, mantle convection, subduction, high-performance computing, finite elements, Schur complement, numerical modeling