xmesh

面向地震动力学模拟的并行三维网格生成程序

概述

xmesh 生成用于 SPECFEM3D 地震波传播求解器的六面体谱元网格。它在 SPECFEM3D 的 xmeshfem3D 网格生成器基础上进行了扩展，支持通过 MPI 实现并行，并通过 Kokkos 实现 GPU 加速。

主要特性：

• 支持层状介质与网格粗化的 Superbrick 网格生成
• SFC-3D 分区算法（Z 加权切分 + Hilbert 曲线水平排序），完全确定性、无外部依赖
• MPI 三维并行分解（列主序约定，xi 最外层）
• 自动物理边界（面 1..6）与 MPI 界面检测
• 兼容 SPECFEM3D 的 Fortran 无格式二进制数据库输出（proc******_Database）
• 通过 Kokkos 实现 GPU 加速的 GLL 节点计算与 Jacobian 求值
• 基于 YAML 的参数输入

依赖

• C++17 编译器（GCC 9+、Clang 10+）
• CMake >= 3.20
• Kokkos（Serial、OpenMP、CUDA 或 HIP 后端）
• MPI（OpenMPI、MPICH 或兼容实现）
• yaml-cpp

可选依赖：

• HDF5（用于替代输出格式）
• ADIOS2（用于可扩展的并行 I/O）

构建

# 使用 Kokkos Serial 后端进行配置
cmake -B build \
    -DKokkos_ROOT=/path/to/kokkos/install \
    -DCMAKE_BUILD_TYPE=Release

# 构建
cmake --build build -j

# 启用单元测试
cmake -B build -DXMESH_BUILD_TESTS=ON
cmake --build build -j

选择 Kokkos 后端

# OpenMP 后端
-DKokkos_ENABLE_OPENMP=ON

# CUDA 后端
-DKokkos_ENABLE_CUDA=ON -DKokkos_ARCH_VOLTA70=ON

# HIP 后端
-DKokkos_ENABLE_HIP=ON -DKokkos_ARCH_VEGA906=ON

用法

# 串行运行
./build/xmesh examples/layered_halfspace/mesh_parameters.yaml

# MPI 并行（4 个进程：2x2 笛卡尔网格）
mpirun -np 4 ./build/xmesh examples/layered_halfspace/mesh_parameters.yaml

# MPI 三维并行（8 个进程：2x2x2 笛卡尔网格）
# 注意：ner_total 必须能被 nproc_z 整除，否则启动时报错
mpirun -np 8 ./build/xmesh examples/layered_halfspace/mesh_parameters.yaml

分区拓扑约束

维度	约束
XI（东西）	nex_xi 必须是 nproc_xi 的整数倍
ETA（南北）	nex_eta 必须是 nproc_eta 的整数倍
Z（垂直）	ner_total 必须是 nproc_z 的整数倍

测试

# 构建并运行所有测试
cmake -B build -DXMESH_BUILD_TESTS=ON
cmake --build build -j
cd build && ctest

# 单独运行单元测试
./test/test_hilbert
./test/test_weight

# MPI 集成测试（不同进程数）
mpirun -np 1 ./test/test_partition_3d    # 1x1x1 拓扑
mpirun -np 2 ./test/test_partition_3d    # 2x1x1 拓扑
mpirun -np 4 ./test/test_partition_3d    # 2x2x1 拓扑
mpirun -np 8 ./test/test_partition_3d    # 2x2x2 拓扑

输入格式

网格参数通过 YAML 文件指定（参见 examples/layered_halfspace/）。

model:
  utm_zone: 11
  x_min: 0.0
  x_max: 50000.0       # XI 方向 50 km
  y_min: 0.0
  y_max: 50000.0       # ETA 方向 50 km
  depth: 30000.0       # 深度 30 km

mesh:
  nex_xi: 64            # XI 方向 64 个谱元
  nex_eta: 64           # ETA 方向 64 个谱元
  ngll: 5               # N=4（5 个 GLL 节点）
  regular: true         # 规则（均匀）网格
  ndoublings: 0         # 无网格粗化

partition:
  nproc_xi: 2           # XI 方向 2 个进程
  nproc_eta: 2          # ETA 方向 2 个进程
  nproc_z: 1            # Z 方向 1 个进程（二维分区）

layers:
  - depth: 15000.0       # 上层 15 km
    material_id: 1
  - depth: 30000.0       # 下层 15 km
    material_id: 2

materials:
  - id: 1
    rho: 2600.0         # kg/m³
    vp: 5800.0          # m/s
    vs: 3400.0          # m/s
    Q: 100.0            # 品质因子

  - id: 2
    rho: 2800.0
    vp: 6500.0
    vs: 3800.0
    Q: 200.0

output:
  format: binary
  path: "./OUTPUT_FILES/DATABASES_MPI/"

架构

src/
├── main.cpp               # 入口点、MPI/Kokkos 初始化、流水线编排
├── parameters.{h,cpp}     # YAML 参数解析与验证
├── geometry.{h,cpp}       # 界面、曲面与层结构
├── meshgen.{h,cpp}        # 网格生成（SFC-3D global_ids 驱动）
├── gll_library.cpp        # GLL 求积节点与 Jacobian 计算
│
├── partition.{h,cpp}      # MPI 笛卡尔拓扑、广播、gather_mesh
├── partition_sfc3d.{h,cpp}# SFC-3D 主控（Z 切 + Hilbert 水平切）
├── partition_z.{h,cpp}    # Z 方向按层界面加权切分
├── partition_weight.{h,cpp}# 单元权值模型（材料/加倍/地形/震源）
├── hilbert.{h,cpp}        # 2D Hilbert 曲线索引
│
├── boundaries.{h,cpp}     # 物理边界（面 1..6）+ MPI 界面检测
├── database_io.{h,cpp}    # SPECFEM3D Fortran 无格式二进制 I/O
├── quality.{h,cpp}        # 网格质量检查与 DT 估计

include/
├── xmesh_config.h         # 编译期常量（默认 NGLL=5、六面体角点数=8）
├── xmesh_types.h          # 核心数据结构（LocalMesh、ProcessInfo 等）
└── xmesh_kernels.h        # Kokkos 内核函数声明

test/
├── unit/
│   ├── test_hilbert.cpp   # HilbertIndex 单元测试
│   └── test_weight.cpp    # WeightModel 单元测试
└── integration/
    └── partition_3d.cpp   # 3D 分区集成测试（4 种 MPI 拓扑）

分区算法（SFC-3D）

MeshParameters → ZAxisPartitioner ─→ Z 切片 {iz_min, iz_max}
                     │
                     ▼
              WeightModel::weight(iex, iey, iez)
                     │
         ┌───────────┴───────────┐
         ▼                       ▼
  水平退路（nex≤1）          Hilbert 路径
  (iex, iey) 行主序均匀切    HilbertIndex::hilbert → 排序
                                     │
                                     ▼
                              加权累计切分
                                     │
                                     ▼
                          Sfc3dPartitioner::compute
                                     │
                                     ▼
                          ElementOwnership.global_ids
                                     │
                                     ▼
                            generate_mesh（只生成 owned 单元）