Pyroclast: A Modular High-Performance Python Solver for Geodynamics

TL;DR

Pyroclast is a modular, high-performance Python framework designed for large-scale geodynamic simulations, combining modern numerical methods, GPU acceleration, and a flexible architecture to modernize legacy codes and enable advanced modeling.

Contribution

It introduces a scalable, GPU-accelerated Python-based geodynamics solver with a modular design, improving performance and extensibility over traditional monolithic codes.

Findings

Achieves 5-10x speedup on NVIDIA A100 GPUs

Demonstrates near-ideal weak scaling up to 896 CPU cores

Provides a flexible, accessible platform for large-scale geodynamics simulations

Abstract

This monograph presents the design, implementation, and evaluation of Pyroclast, a modular high-performance Python framework for large-scale geodynamic simulations. Pyroclast addresses limitations of legacy geodynamics solvers, often implemented in monolithic Fortran, C++, or C codebases with limited GPU support and extensibility, by combining modern numerical methods, hardware-accelerated execution, and a flexible object-oriented architecture. Designed for distributed and GPU-accelerated environments, Pyroclast provides an accessible and efficient platform for simulating mantle convection and lithospheric deformation using the marker-in-cell method and a matrix-free finite difference discretization. The work focuses on a scalable two-dimensional viscous mechanical solver that forms the computational core for future visco-elasto-plastic models. The solver includes a stress-conservative staggered grid discretization of the incompressible Stokes equations, a matrix-free geometric multigrid solver, Krylov and quasi-Newton methods, and MPI-based domain decomposition for distributed execution. Benchmarks evaluate performance and scalability. Shared-memory tests show strong scaling of the Stokes solver and demonstrate a 5-10x speedup on NVIDIA A100 GPUs compared to a multi-core CPU baseline. Distributed advection benchmarks show near-ideal weak scaling up to 896 CPU cores across seven compute nodes. These results demonstrate that Pyroclast achieves high performance while remaining accessible through a high-level Python interface. The framework also provides a blueprint for modernizing legacy geodynamics codes. Its modular architecture and Python-native implementation lower the barrier to entry while enabling interoperability with modern machine learning libraries, enabling hybrid physics-based and data-driven workflows.