[论文解读] Defmod - Parallel multiphysics finite element code for modeling crustal deformation during the earthquake/rifting cycle

Shehzad Ali|arXiv (Cornell University)|Feb 3, 2014

Advanced Numerical Methods in Computational Mathematics参考文献 11被引用 11

一句话总结

Defmod 是一个用 Fortran 95 编写的开源、并行、非结构化有限元代码，通过多物理场模拟对地震和裂谷循环中的地壳变形进行建模。它基于 PETSc 求解器，支持动态和准静态过程，包括同震滑动、岩浆侵入、孔弹性回弹和黏弹性松弛，可在数千个核心上实现可扩展的模拟，具有较高的并行效率，并采用稳定化的线性单元以实现精确的孔弹性建模。

ABSTRACT

In this article, we present Defmod, an open source, fully unstructured, two or three dimensional, parallel finite element code for modeling crustal deformation over time scales ranging from milliseconds to thousands of years. Unlike existing public domain numerical codes, Defmod can simulate deformation due to all major processes that make up the earthquake/rifting cycle, in non-homogeneous media. Specifically, it can be used to model deformation due to dynamic and quasistatic processes such as co-seismic slip or dike intrusion(s), poroelastic rebound due to fluid flow and post-seismic or post-rifting viscoelastic relaxation. It can also be used to model deformation due to processes such as post-glacial rebound, hydrological (un)loading, injection and/or withdrawal of fluids from subsurface reservoirs etc. Defmod is written in Fortran 95 and uses PETSc's parallel sparse data structures and implicit solvers. Problems can be solved using (stabilized) linear triangular, quadrilateral, tetrahedral or hexahedral elements on shared or distributed memory machines with hundreds or even thousands of processor cores. In the current version of the code, prescribed loading is supported. Results are written in ASCII VTK format for easy visualization. The source code is released under the terms of GNU General Public License (v3.0) and is freely available from https://bitbucket.org/stali/defmod/.

研究动机与目标

• 开发一种灵活的开源有限元代码，能够在多种时间尺度上模拟复杂的地壳变形过程。
• 在单一统一框架内实现多种相互关联的地球物理过程（如地震滑动、流体引起的孔弹性回弹和黏弹性松弛）的模拟。
• 通过非结构化网格和稳定化有限元，支持具有可变流变性质和复杂几何形状的非均质介质。
• 利用 PETSc 的并行稀疏求解器，在共享内存或分布式内存架构上实现高性能、可扩展的模拟。
• 创建一个不仅适用于研究，也适用于教学和教育用途的工具，因其简洁的 ~2000 行 Fortran 95 实现而具备优势。

提出的方法

• Defmod 使用有限元法求解运动的柯西方程，将弱形式离散化为半离散系统： Mü + Cü + Ku = f。
• 它采用拉格朗日乘子通过增强系统方程强制实现断层滑动、位移和速度约束： Mü + Cu̇ + Ku + Gᵀλ = f 和 Gu = l。
• 对于动态波传播，采用显式中心差分格式计算加速度，使用后向差分格式计算速度，并采用集中质量矩阵以实现高效求逆。
• 采用瑞利阻尼（C = αM + βK）控制数值耗散，时间步长受 CFL 条件限制：Δt ≤ L/c，其中 c = √(E/ρ)。
• 对于准静态问题，使用隐式时间积分，结合 PETSc 的迭代或直接求解器，求解孔弹性的耦合动量和连续性方程。
• 实现局部压力投影方案（Bochev & Dohrmann, 2006），以满足线性单元在孔弹性问题中的 inf-sup 条件。

实验结果

研究问题

• RQ1 一个单一的有限元代码能否高效模拟地震/裂谷循环中的多种相互关联过程，包括动态破裂、流体流动和黏弹性松弛？
• RQ2 Defmod 在数千个处理器核心上的性能如何扩展，适用于显式弹性和隐式准静态问题？
• RQ3 与高阶方法相比，Defmod 在孔弹性问题中采用稳定化线性单元的实现，在精度和稳定性方面表现如何？
• RQ4 一个紧凑的 ~2000 行 Fortran 95 代码在大规模地球物理模拟中，能否实现高性能和可扩展性？
• RQ5 基于 Lysmer 和 Kuhlemeyer 的阻尼器的波吸收方案在多大程度上能有效减少域边界处的虚假反射？

主要发现

• Defmod 在最多 1024 个核心上对显式弹性动力学和隐式准静态问题均实现了 0.90 或以上的并行效率，且在高核心数量下性能保持稳定。
• 对于显式弹性动力学问题，当核心数从 512 增加到 1024 时，效率从 0.93 降低到 0.80，主要原因是每个核心的计算工作量减少（降至 24.4K DOFs/core）。
• 由于每个核心额外的局部计算（如应力恢复和载荷向量组装），隐式准静态求解器的效率更高（在 1024 核时为 0.87）。
• 该代码成功模拟了采用稳定化线性单元的孔弹性回弹，避免了锁定现象，且无需使用高阶单元即可保证精度。
• 在 1024 核上，隐式求解器在第 500 次迭代时的最终预条件残差范数为 4.965×10⁻⁵，表明收敛性稳健。
• Defmod 在包含约 2500 万个自由度的 3D Lamb 问题上表现出强可扩展性，证实其具备在现代 HPC 集群上处理大规模地球物理模拟的能力。