Source-linked AI summary
A Domain Decomposition Approach to Implementing Fault Slip in Finite-Element Models of Quasi-static and Dynamic Crustal Deformation
Brad T. Aagaard, Matthew G. Knepley, Charles A. Williams
TL;DR
Modeling the earthquake cycle requires representing slow deformation and rapid rupture without simplifying away their interactions. This paper implements fault slip in PyLith through domain decomposition with Lagrange multipliers, develops a custom fault-block preconditioner, and benchmarks quasi-static and dynamic cases. The implementation agrees with analytical and other-code benchmarks, while scalability is limited principally by load imbalance and nonlinear-iteration sensitivity to fault behavior.
Problem
Earthquake-cycle studies often simplify slow deformation, fault loading, or rupture propagation, motivating models that capture interseismic deformation and rapid rupture together.
Method
PyLith uses a domain decomposition formulation with Lagrange multipliers to implement prescribed or constitutive fault slip across quasi-static and dynamic finite-element simulations.
Results
The fault implementation matches analytical viscoelastic-relaxation results and other dynamic rupture codes, while the custom preconditioner reduces iterations by 80%.
Takeaways & Limitations
A common PyLith framework supports complex-geometry fault slip and integrated quasi-static and dynamic crustal-deformation modeling.
Takeaways & Limitations
Parallel scalability is chiefly constrained by load imbalance, and spontaneous-rupture iteration can be less efficient when deformation is not confined near the fault.
Abstract
from arXiv · showhide
We employ a domain decomposition approach with Lagrange multipliers to implement fault slip in a finite-element code, PyLith, for use in both quasi-static and dynamic crustal deformation applications. This integrated approach to solving both quasi-static and dynamic simulations leverages common finite-element data structures and implementations of various boundary conditions, discretization schemes, and bulk and fault rheologies. We have developed a custom preconditioner for the Lagrange multiplier portion of the system of equations that provides excellent scalability with problem size compared to conventional additive Schwarz methods. We demonstrate application of this approach using benchmarks for both quasi-static viscoelastic deformation and dynamic spontaneous rupture propagation that verify the numerical implementation in PyLith.
1. Introduction
Earthquake-cycle models must span slow interseismic deformation and rapid rupture across broad spatial and temporal scales. PyLith develops fault-slip capabilities intended to reduce simplifications and support integrated quasi-static and dynamic simulations.
- Earthquake-cycle deformation spans fractions of a meter to thousands of kilometers and fractions of a second to thousands of years.These scales encompass asperities and grains through plate boundaries, and rupture slip through interearthquake strain accumulation.
- Many studies simplify one part of the earthquake cycle, such as replacing dynamic rupture with static coseismic slip in slow-deformation models.Such simplifications are used to make calculations tractable or expedite modeling of another portion of the cycle.
- Some earthquake-cycle simulations resolve both interseismic deformation and dynamic rupture, revealing effects that require both slow and rapid processes.Examples include spatially varying fault properties, bends, low-rigidity layers, and damaged zones.
- Reproducing observed earthquake-cycle behavior remains a challenge despite advances in models spanning broader space-time windows.
- Desirable earthquake-cycle models combine elastic, viscoelastic, and viscoelastoplastic rheologies with prescribed or spontaneous fault slip.The broader modeling framework may also couple elasticity to fluid or heat flow.
- PyLith focuses on fault slip for quasi-static interseismic and coseismic deformation and dynamic rupture, with plans to couple these simulations across the earthquake cycle.The approach addresses potentially nonplanar, intersecting faults and permits either prescribed slip or constitutive-model-controlled slip.
2. Numerical Model of Fault Slip
The numerical model treats the fault as an interior frictional interface and augments finite-element elasticity with Lagrange multipliers for fault tractions and slip constraints. It supports prescribed or constitutive fault slip, while specialized solver and line-search strategies address scale and nonlinear convergence challenges.
- The formulation augments conventional finite-element elasticity with domain decomposition to implement fault slip.
- The fault is an interior boundary whose Lagrange multipliers represent equal-and-opposite tractions enforcing the displacement-jump and slip constraint.Slip is defined as displacement of the positive fault side relative to the negative side.
- Prescribed-slip simulations specify d on the fault, whereas constitutive simulations specify how l depends on slip, slip rate, and state variables.
- The resulting system solves jointly for displacement and Lagrange multipliers, whose SI magnitudes can differ by up to nine orders.Displacements are generally mm to m, while fault tractions are generally MPa.
- The iterative procedure is less efficient when friction varies significantly with slip or deformation extends across the entire domain.Its effectiveness depends on the constitutive model and on deformation remaining concentrated near the fault.
- The spontaneous-rupture iteration uses a logarithmic-space line search to select slip increments that best satisfy the fault constitutive model.This greatly accelerates convergence when friction depends on the logarithm of slip rate.
3. Finite-Element Mesh Processing
PyLith represents faults as interior surfaces in a finite-element mesh and preprocesses them into cohesive-cell structures that duplicate fault-side vertices and support Lagrange-multiplier constraints. Mesh topology is updated to maintain consistent connectivity around the fault.
- PyLith uses PETSc to perform finite-element cell and face integrals, assemble global vectors and matrices, impose Dirichlet conditions, and solve nonlinear systems.
- The Sieve mesh API uses directed acyclic graph covering relations and supports one-, two-, and three-dimensional meshes with multiple element types.
- The fault must align with element faces, and each fault vertex needs displacement degrees of freedom on both sides.
- PyLith replaces each fault face with a zero-volume cohesive cell during preprocessing.The preprocessing constructs consistently oriented fault faces from vertices marked as lying on the fault.
- For each fault vertex, preprocessing introduces positive-side, negative-side, and Lagrange-multiplier vertices.The multiplier vertex lies on an edge between the two fault-side vertices.
- Cells touching the fault through an edge or vertex are reclassified and their positive-side vertices are replaced to preserve the modified fault topology.Classification uses neighboring cells and a breadth-first search when necessary.
4. Solver Customization
PyLith customizes solvers for fault-slip systems by combining Krylov methods, AMG for elasticity, and structure-exploiting treatments for Lagrange multipliers. Diagonal approximations and Schur-complement updates enable efficient quasi-static and dynamic solves while preserving accuracy at resolved length scales.
- Quasi-static simulations: Preconditioned Krylov methods solve the large sparse quasi-static systems, with a preconditioner specialized to fault slip with Lagrange multipliers.Sparse matrix-vector multiplication supports scalable parallel processing, but convergence requires problem-specific preconditioning.
- Quasi-static simulations: Introducing Lagrange multipliers creates a saddle-point system for which traditional additive Schwarz preconditioning converges slowly.PETSc field-splitting tools support effective preconditioners for the separate physical fields.
- Quasi-static simulations: Collocated vertex quadrature makes the fault coupling matrix block diagonal, simplifying the Lagrange-multiplier portion of the preconditioner.For spectral elements, Legendre polynomials with Gauss-Lobatto-Legendre points produce the same block-diagonal structure.
- Quasi-static simulations: The preferred configuration combines AMG for the elasticity submatrix with a custom preconditioner for the Lagrange multiplier submatrix.The custom fault preconditioner exploits the sparse Jacobian structure and supports a simple block treatment through the Schur complement.
- Dynamic simulations: Dynamic solves exploit diagonal Jacobian blocks and Schur-complement updates to compute local corrections for fault degrees of freedom and Lagrange multipliers.The diagonal approximation is accurate when deformation length scales are significantly larger than the discretization size.
5. Performance Benchmark
The benchmark evaluates preconditioners and parallel performance for increasingly large multi-fault finite-element problems. The custom Lagrange-multiplier preconditioner improves convergence and runtime, while scalability is limited mainly by load imbalance, memory bandwidth, and AMG setup costs.
- Preconditioner performance: The benchmark uses three intersecting vertical strike-slip faults and compares hexahedral and tetrahedral meshes with nearly equal problem sizes.Problem sizes range from 7.8×10^4 DOF on one process to 7.1×10^6 DOF on 96 processes.
- Preconditioner performance: An 80% reduction in iterations makes the custom Lagrange-multiplier preconditioner the fastest tested preconditioner.Iteration growth is about 20% for a fourfold DOF increase with Schur-complement or AMG field-split methods, versus 60% for ASM.
- Parallel scaling performance: Load imbalance is the main scalability impediment, reaching up to 20% in matrix-vector multiplication on 96 cores.Its sources include unstructured-mesh partitioning, cell-based partitioning, and equal weighting of cohesive and bulk cells.
- Parallel scaling performance: Single-node scaling plateaus from 2 to 4 processes because memory bandwidth saturates, although scaling improves from 12 to 24 processes.The bandwidth limitation depresses overall performance but is not expected to affect inter-node scaling.
- Parallel scaling performance: AMG setup and application dominate solver time, and setup growth causes most deviation from perfect weak scaling.Reducing coarse-grid work could improve scalability but would increase iterations and overall solver time.
6. Code Verification Benchmarks
PyLith is verified against analytical and cross-code benchmarks spanning quasi-static viscoelastic deformation and dynamic spontaneous rupture. The results show convergence with refinement and close agreement with reference solutions and other modeling codes.
- Quasi-static benchmark: The quasi-static benchmark compares surface displacement with the Savage–Prescott analytical solution after multiple earthquake cycles.The numerical solution approaches steady plate motion only after several cycles.
- Quasi-static benchmark: By the tenth earthquake cycle, the numerical surface displacement matches the analytical solution, whereas the second-cycle far field underpredicts displacement.The mismatch reflects the transient approach to steady plate motion in the numerical model.
- Quasi-static benchmark: 6.7 km discretization produces a close analytical fit, with the hexahedral solution indistinguishable and the tetrahedral solution slightly low late in the cycle.Coarse models show large errors near the fault, within about one elastic thickness.
- Dynamic benchmark: The TPV13 dynamic benchmark shows convergence as fault resolution decreases, while 200 m meshes exhibit high-frequency oscillations from insufficient cohesive-zone resolution.Triangular cells generate less oscillation than quadrilateral cells.
- Dynamic benchmark: PyLith slip-rate histories agree very well with four other dynamic rupture codes, with remaining oscillation differences attributed to numerical damping.The comparison is used because TPV13 has no analytical solution.
- Dynamic benchmark: The 3-D TPV13 benchmark produces similar results to the 2-D case but requires greater computational resources.At 100 m resolution, the simulations took 2.5 hours using 64 processes.
- Dynamic benchmark: Across TPV13, PyLith performs similarly to other finite-element and finite-difference codes and captures supershear rupture with the specified bulk and fault rheologies.The benchmark also tests complex geometry and varying discretization size on a dipping normal fault.
7. Conclusions
The paper presents a scalable domain-decomposition implementation of fault slip in PyLith and verifies it for quasi-static and dynamic applications. The benchmarks support its use for coupled earthquake-cycle modeling.
- Method: PyLith combines algebraic multigrid for elasticity with a custom Lagrange-multiplier fault preconditioner to accelerate Krylov convergence with low sensitivity to problem size.The approach uses Lagrange multipliers to represent fault tractions.
- Verification: The fault-slip implementation agrees with an analytical viscoelastic relaxation solution and with other codes for supershear rupture on a dipping normal fault.These benchmarks cover quasi-static and dynamic deformation.
- Implication: The methodology is proposed as an avenue for coupling quasi-static interseismic and postseismic deformation with dynamic earthquake rupture simulations.
Notation
The notation defines the Jacobian operators, material and fault variables, and preconditioning matrices used in the finite-element formulation.
- Operators: The notation includes the Jacobian operators for the full system, elasticity equation, and fault constraint equation.
- Variables: The material and fault variables include the elastic-constant tensor C, fault slip vector d, body force vector f, and Lagrange multiplier vector l.The Lagrange multiplier corresponds to the fault traction vector.
- Preconditioning: The notation also defines P as the preconditioning matrix, with P_elastic and P_fault associated with elasticity and fault-slip constraints.