High-performance geodynamics on GPUs using the PETSc CUDA backend (GAUZZ v1.0.0)
- 1Kangwon National University (Chuncheon, Korea, Republic of)
Abstract
GPUs provide much higher memory bandwidth than CPUs. Data-parallel computations with regular memory access patterns exploit the GPU architecture efficiently. To solve the Stokes equations, we apply a preconditioned conjugate-gradient (PCG) method to the 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 the solver in GAUZZ v1.0.0 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. The P ₁ mass matrix, the gradient operator, and the divergence operator are reused in the pressure-correction L ₂ -projection step. HYPRE BoomerAMG was adopted for preconditioning the velocity-only subproblem. Accuracy and performance are evaluated on a manufactured solution, the SolCx benchmark, the sticky-air benchmark, 2D Rayleigh–Taylor instabilities, 3D thermo-mechanical convection, and a 2D subduction model with nonlinear viscosity. We confirm that the GPU-based implementation reproduces CPU solutions. Performance is quantified by the execution time of the pressure-correction step and by throughput, defined as the number of Stokes degrees of freedom divided by the execution time of one Stokes solve. A single GPU delivers 6–11× the throughput of a 32-core CPU across the Stokes benchmarks. Two GPUs increase the throughput at sufficiently large problem sizes, reaching up to 2.1× the single-GPU throughput. The two-GPU configuration also extends the maximum problem size to 2.5×10 ⁷ DOFs, beyond the 48 GB VRAM limit of a single GPU.
Versions
Related
- GeoSim: Benchmarking Multi-Field Neural Operators for Geodynamics
- GeoAgent: An Agent-Based QGIS Plugin for Natural-Language-Driven Geospatial Analysis with Multi-Backend Large Language Model Support
- Response-surface optimization of high-performance GO-PET-PEG forward osmosis membranes for bentazone removal
- Visible light -driven NiCo2O4 nano-urchins: A high-performance photocatalyst for ciprofloxacin degradation via peroxymonosulfate enhancement
- Speeding up the GEOS-Chem High Performance CTM (GCHP version 14.7.1) through MPI load balancing
BibTeX
@article{lee2026highperformance,
title = {{High-performance geodynamics on GPUs using the PETSc CUDA backend (GAUZZ v1.0.0)}},
author = {Kyeong-Min Lee and Deok-Kyu Jang and Cedric Thieulot and Byung-Dal So},
journal = {Geoscientific Model Development},
year = {2026},
month = oct,
doi = {10.5194/gmd-19-9781-2026},
url = {https://gmd.copernicus.org/articles/19/9781/2026/},
note = {CC BY 4.0},
}
RIS
TY - JOUR
TI - High-performance geodynamics on GPUs using the PETSc CUDA backend (GAUZZ v1.0.0)
AU - Kyeong-Min Lee
AU - Deok-Kyu Jang
AU - Cedric Thieulot
AU - Byung-Dal So
PY - 2026
DA - 2026/10/13
JO - Geoscientific Model Development
DO - 10.5194/gmd-19-9781-2026
UR - https://gmd.copernicus.org/articles/19/9781/2026/
AB - GPUs provide much higher memory bandwidth than CPUs. Data-parallel computations with regular memory access patterns exploit the GPU architecture efficiently. To solve the Stokes equations, we apply a preconditioned conjugate-gradient (PCG) method to the 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 the solver in GAUZZ v1.0.0 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. The P ₁ mass matrix, the gradient operator, and the divergence operator are reused in the pressure-correction L ₂ -projection step. HYPRE BoomerAMG was adopted for preconditioning the velocity-only subproblem. Accuracy and performance are evaluated on a manufactured solution, the SolCx benchmark, the sticky-air benchmark, 2D Rayleigh–Taylor instabilities, 3D thermo-mechanical convection, and a 2D subduction model with nonlinear viscosity. We confirm that the GPU-based implementation reproduces CPU solutions. Performance is quantified by the execution time of the pressure-correction step and by throughput, defined as the number of Stokes degrees of freedom divided by the execution time of one Stokes solve. A single GPU delivers 6–11× the throughput of a 32-core CPU across the Stokes benchmarks. Two GPUs increase the throughput at sufficiently large problem sizes, reaching up to 2.1× the single-GPU throughput. The two-GPU configuration also extends the maximum problem size to 2.5×10 ⁷ DOFs, beyond the 48 GB VRAM limit of a single GPU.
SN - 1991-9603
ER -