Versions

v1.0.0 (2025-09-30)

We are pleased to announce the first public release of PyQuake3D – a Python-based simulation framework for modeling 3-D earthquake sequences, capturing both seismic and aseismic fault slip.

Key Features

  • Simulation of multi-cycle earthquake behavior on 3-D faults

  • Support for both seismic rupture and aseismic creep

  • GPU acceleration using CuPy

  • MPI parallelization for large-scale computations

  • Modular finite-difference / finite-volume numerical implementation

  • Geometry-based input using .msh files

  • Parameter-driven simulation control (parameter.txt)

  • Diagnostic output and visualization tools via HDF5 and Matplotlib

  • Includes benchmark examples such as SCEC BP5-QD

Example Simulations

v1.0.1 (2025-10-06)

Release Summary

Version 1.0.1 introduces major improvements in physics modeling, MPI scalability, and data handling. This update expands PyQuake3D’s ability to simulate fully coupled earthquake cycles while streamlining the MPI code structure and optimizing outputs for large-scale workflows.

What’s Changed

  1. Dilatancy Law Implementation PyQuake3D now includes a fully coupled dilatancy law. The earthquake-cycle model combines 3-D elasticity, rate-and-state friction, and porosity evolution using a Boundary Integral Equation Method (BIEM) coupled to a Finite Difference Method (FDM). Pore-fluid diffusion parallel to the fault is neglected. The FDM time-marching scheme is fully MPI-parallelized.

  2. Enhanced Modularity (MPI Version) The MPI version has been reorganized for clarity. The main driver now uses only six core functions, simplifying the execution pipeline and improving maintainability.

  3. Comprehensive MPI Acceleration MPI is now applied to all computational components, beyond matrix–vector products and Green’s function evaluation. This delivers improved scaling and faster performance on distributed-memory HPC systems.

  4. Optimized Output Format (VTU) Simulation output is now written in VTU format. File sizes are reduced by approximately threefold compared to the previous VTK format, improving storage efficiency and I/O throughput.

  5. Updated `parameter.txt` Several parameter names related to stress and output configuration have been revised. Users should update existing parameter.txt files to align with the new naming conventions.

v1.0.2 (2025-12-01)

Release Summary

Version 1.0.2 introduces two major advancements to PyQuake3D: (1) a scalable Lattice H-matrix architecture for distributed-memory HPC systems, and (2) the implementation of thermal pressurization, enabling thermo–hydro–mechanical weakening during dynamic rupture. Together, these additions greatly enhance PyQuake3D’s capability to model multi-physics earthquake processes at scale.

What’s Changed

  1. Lattice H-matrix Implementation

    PyQuake3D now includes a Lattice H-matrix framework, a distributed-memory extension of the standard H-matrix format designed for large-scale Boundary Element Method (BEM) simulations.

    Lattice H-matrices improve parallel scalability by organizing matrix blocks into a 2D MPI “lattice” layout, reducing communication overhead and enabling efficient compression and multiplication for extremely large, dense BEM matrices.

    Context and Motivation

    • In BEM, the stiffness matrix is fully dense, and matrix–vector products dominate computational cost.

    • Standard H-matrices compress these matrices but require global MPI reductions at each Runge–Kutta stage, leading to prohibitive communication cost as process count increases.

    • In PyQuake3D’s conventional H-matrix implementation, each time step requires: 1. dot products between local H-matrix subblocks and slip rates, 2. MPI_Reduce to gather contributions across all ranks, 3. another MPI_Reduce to compute global slip rates.

      Communication therefore scales like O(n²) with increasing MPI processes.

    Lattice H-matrix Improvement

    • Lattice H-matrices reorganize data into a 2D MPI grid.

    • Dot-product accumulation and slip-rate updates occur locally within each MPI row, independent of other rows.

    • Only diagonal processes perform global reductions at the end of RK stages.

    • This reduces communication complexity from O(n²) to O(n).

    Impact

    • Greatly improved strong-scaling efficiency for large MPI jobs

    • Lower communication overhead

    • Faster time integration for large 3D earthquake-cycle models

    • Enhanced suitability for Tier-1/Tier-0 supercomputers

  2. Thermal Pressurization

    PyQuake3D now supports thermal pressurization (TP) as a dynamic weakening mechanism during seismic slip.

    During rapid coseismic slip, frictional heating raises pore fluid temperature, producing pore pressure that reduces effective normal stress and dramatically weakens the fault. This process can govern rupture speed, stress drop, and arrest behavior.

    The governing equations for pore fluid pressure \(p\) and temperature \(T\) follow [16, 17]:

    (59)\[\frac{\partial p}{\partial t} = c_{hy} \frac{\partial^2 p}{\partial y^2} + \Lambda \frac{\partial T}{\partial t}\]
    (60)\[\frac{\partial T}{\partial t} = \alpha_{th} \frac{\partial^2 T}{\partial y^2} + \frac{\tau V \exp\left(-y^2 / (2 w^2)\right)}{\rho c \sqrt{2 \pi} w}\]

with boundary conditions:

Thermal pressurization only

(61)\[\left. \frac{\partial p}{\partial y} \right|_{y=0^+} = 0, \qquad \left. \frac{\partial T}{\partial y} \right|_{y=0^+} = 0\]

Coupled dilatancy + thermal pressurization

(62)\[\left. \frac{\partial p}{\partial y} \right|_{y=0^+} = \frac{h \dot{\phi}}{2 \beta c_{hyd}}, \qquad \left. \frac{\partial T}{\partial y} \right|_{y=0^+} = 0\]

where:

  • \(T\) — pore fluid temperature

  • \(p\) — pore pressure

  • \(c_{hy}\) — hydraulic diffusivity

  • \(\alpha_{th}\) — thermal diffusivity

  • \(w\) — shear zone half-width

  • \(\rho c\) — volumetric heat capacity

  • \(\Lambda\) — thermal pressurization coefficient

  • \(\tau V\) — frictional heating source term distributed across the shear zone

Physical Interpretation

  • Dilatancy increases porosity → decreases pore pressure → strengthens the fault.

  • Thermal pressurization increases temperature → raises pore pressure → weakens the fault.

These processes compete, producing a fundamental stability trade-off in fault mechanics.

v1.0.3 (2026-4-28)

## What’s Changed

  1. New implementation of MPI-based multi-GPU backend: We implemented a hybrid computing framework where CPUs orchestrate workload distribution, Rugge-Kutta iteration, and communication, while GPUs are dedicated to accelerating computationally intensive numerical kernels (e.g., large-scale Hmatrix operations). To optimize performance, the H-matrix structure undergoes a series of specialized transformations for GPU compatibility. This module provides native support for both standard and Lattice H-matrices.

  2. Pyquake_tools: The pyquake_tools in the root dir provides tools for visualizing and analyzing PyQuake3D simulation outputs, including:

  • Time series plots of maximum slip rate

  • 2D and 3D animations of fault slip

  • Automatic earthquake event detection

  • Slip statistics and seismic moment calculations