SPECFEM手册的一些翻译（Chapter 3）

scouttttt

已于 2024-09-23 16:17:50 修改

阅读量426

点赞数

文章标签： github

于 2024-09-23 10:55:03 首次发布

原文链接：https://github.com/SPECFEM

版权

Chapter 3 Mesh Generation

3.1 How to use SPECFEM2D

Notes about DATA/Par_file parameters

3.2 如何使用 Gmsh 生成外部网格

3.3 如何使用 Cubit/Trelis 生成外部网格

3.3.1 关于 Cubit/Trelis 内置 Python 的注意事项

3.4 关于吸收 PML（完美匹配层）的注意事项

Chapter 3 Mesh Generation

3.1 How to use SPECFEM2D

运行谱元模拟的第一步是为研究区域构建一个高质量的网格。我们提供了两种方法来实现这一点：(1) 依赖外部的六面体网格生成器，如 Gmsh 和 CUBIT，或 (2) 使用提供的内部网格生成功能。要运行网格生成器，请编辑主输入文件 DATA/Par_file，该文件描述了模拟的设置。然后在终端中输入以下命令：

./bin/xmeshfem2D

这将创建网格，并将网格文件存储在 OUTPUT_FILES/ 目录中。xmeshfem2D 会输出若干个名为 Database??????.bin 的文件，每个文件对应一个进程。这些文件是后续运行求解器时所必需的。

Notes about DATA/Par_file parameters

代码根目录中提供的默认 DATA/Par_file 包含详细的注释，并且几乎是自解释的（请注意，EXAMPLES 目录中提供的一些较旧的 DATA/Par_file 文件虽然可以正常工作，但其中的某些注释可能已经过时甚至不正确，因此，请参考默认的 DATA/Par_file 以获得可靠的解释）。

如果你需要更多详细信息，我们没有为 2D 版本提供所有参数的详细说明，但你可以在 3D 版本的手册中找到有用的信息，因为许多参数和整体设计理念是相似的。你可以在 SPECFEM3D 用户手册 中找到相关内容。

NPROC：

如果你使用了支持 MPI（消息传递接口）编译的代码，可以为模拟指定多个进程（>1）。否则，必须将其设置为 1 以进行串行运行。网格生成器将使用由参数 PARTITIONING_TYPE 指定的分区算法来对多个进程进行负载均衡。

NGNOD：

关于网格点的编号，在网格生成器创建的文件中，我们使用经典的 4 节点和 9 节点有限元的约定（分别为 NGNOD = 4 或 NGNOD = 9）。

MODEL：

当选择 MODEL = default 时，可以在 Par_file 中的 nbmodels 部分定义多种简单的速度和密度模型。

材料类型可以通过以下格式之一来指定：

I: model_number 1 rho Vp Vs 0 0 QKappa Qmu 0 0 0 0 0 0
II: model_number 2 rho c11 c13 c15 c33 c35 c55 c12 c23 c25 0 QKappa Qmu
III: model_number 3 rhos rhof phi c kxx kxz kzz Ks Kf Kfr etaf mufr Qmu
IV: model_number -1 0 0 A 0 0 0 0 0 0 0 0 0 0

材料类型定义说明：

I: 要将某个区域设置为声学模型，使用 (I) 并将 Vs 设置为 0。要将某个区域设置为各向同性弹性体，使用 (I) 并将 Vs 设置为非零值。更多细节参见第 4.1 节。因此，要创建声学（流体）区域，只需将 S 波速度设为 0，程序将自动识别这些区域为流体，并切换到相应的方程，并与固体区域自动匹配。
II: 要将某个区域设置为各向异性材料，使用 (II)。更多细节参见第 4.3 节。
III: 要将某个区域设置为孔弹性材料，使用 (III)。更多细节参见第 4.4 节。
IV: 要将某个区域设置为层析模型，使用 (IV)。目前仅允许一个层析文件和区域。注意，model_number 必须为正数，但对于层析模型，接下来的区域编号必须为负数（-1）。密度（rho）和 P 波速度（Vp）的值将被忽略，但剪切波速度（Vs，即 A）的值必须为 0.0 才能被识别为声学区域，或者为正的非零值（如 1.0）才能识别为弹性区域。层析模型的值将覆盖该区域，并在 TOMOGRAPHY_FILE 参数指定的文件中定义。

当开启粘弹性模型时，从此处读取的 P 波速度 (Vp) 和 S 波速度 (Vs) 值为未弛豫值，即在无限频率下的值，除非参数 READ_VELOCITIES_AT_f0 设置为 true，此时读取的将是频率 f0 下的值。请注意，Qmu 始终等于 Qs，但 Qkappa 通常不等于 Qp。要进行相互转换，请参考文档 doc/note_on_Qkappa_versus_Qp.pdf 和脚本 utils/attenuation/conversion_from_Qkappa_Qmu_to_Qp_Qs_from_Dahlen_Tromp_959_960.f90。

当选择 MODEL = default 时，可以通过 Par_file 文件中的 nbregions 部分指定各种简单的分层模型配置。请查看 Par_file 中的注释，了解如何为内部网格指定不同的材料区域。对于由外部网格生成器定义的网格，必须在单独的文件中指定这些信息。关于 read_external_mesh 参数的更多解释请见下文。

MODEL 参数还可以设置为除 default 之外的其他值，可能的选项包括：
ascii：用于读取单个 ASCII 模型文件（DATA/proc***_rho_vp_vs.dat）。
binary 或 gll：用于读取二进制模型文件（DATA/proc***_rho.bin, DATA/proc***_vp.bin, DATA/proc***_vs.bin）。如果模拟启用了衰减，也会读取衰减文件（DATA/proc***_Qmu.bin, DATA/proc***_Qkappa.bin）。
binary_voigt：用于读取由二进制文件指定的各向异性模型（DATA/proc***_rho.bin, DATA/proc***_c11.bin, DATA/proc***_c13.bin, DATA/proc***_c15.bin, DATA/proc***_c33.bin, DATA/proc***_c35.bin, DATA/proc***_c55.bin）。
external：从源代码例程 src/specfem2d/define_external_model.f90 中的 define_external_model() 读取外部定义的速度模型。用户可以按照源文件中的说明在此处实现自己的例程。默认情况下，该例程实现的是球对称各向同性的 AK135 模型。
legacy：用于读取非常旧的 SPECFEM2D 版本的单个 ASCII 模型文件（DATA/proc***_model_velocity.dat_input）。

Read_external_mesh：

如果你使用外部网格生成器（如 Gmsh、CUBIT/Trelis 或 GiD），应将 read_external_mesh 参数设置为 .true.，并定义以下文件：

mesh_file：描述网格的文件：第一行为元素的数量，接着每行列出构成每个四边形单元的 4 个节点。
nodes_coords_file：包含每个节点坐标的文件：第一行为节点数，接着每行列出每个节点的 x 和 z 坐标。
materials_file：描述每个单元的材料编号：每行一个整数，范围为 1 到 nbmodels。
nummaterial_velocity_file：指定 nbmodels 中每种材料的材料属性。其格式略有不同，但与 Par_file 中的 nbmodels 部分相似，格式与 SPECFEM3D_Cartesian 版本相同。该文件的示例可在 EXAMPLES/ 文件夹中找到，文件头中有详细说明。
free_surface_file：描述声学自由表面边缘的文件：第一行为边缘数量，接着每行列出：单元编号、自由表面的节点数（1 表示点，2 表示边）、形成自由表面的节点。如果不需要自由表面，第一行写 0，此时将获得一个刚性表面。
axial_elements_file：在轴对称模拟情况下，描述轴向单元的文件。参见第 4.2 节。
absorbing_surface_file：描述吸收边界边缘的文件：第一行为边缘数量，接着每行列出：单元编号、组成吸收边的节点数（始终为 2），形成吸收边的两个节点，然后是吸收边的类型：1 表示底部 (BOTTOM)，2 表示右侧 (RIGHT)，3 表示顶部 (TOP)，4 表示左侧 (LEFT)。每行的第二个参数必须为 2，如果某个单元沿给定吸收轮廓有多个边，应将其列出多次。避免重复列出相同元素及其吸收边，以确保吸收正确进行。如果某个单元仅在吸收轮廓上有一个点，不应列出它。如果使用 9 节点单元，只需列出边缘的第一个和最后一个点，代码会自动恢复正确的 9 节点曲率。
tangential_detection_curve_file：包含用于 source_normal_to_surface 和 rec_normal_to_surface 的包络曲线点的文件。应精细排列，并按顺时针方向排序。第一行为点的数量，接着每行列出 (x, z) 坐标。

Read_velocities_at_f0`：`

参数用于调整从输入文件中读取的速度值，以考虑平均物理色散的影响。也就是说，如果需要，可以更改代码内部定义速度的参考频率：

默认情况下，代码读取的速度值是未弛豫值，即无限频率下的速度值。
如果将该标志设置为 .true.，则从输入文件中读取的速度值将对应于一个特定频率，称为 ATTENUATION_f0_REFERENCE。

请注意，这个标志仅在启用了衰减功能的模拟中生效，即当 ATTENUATION_VISCOELASTIC 或 ATTENUATION_VISCOACOUSTIC 设置为 .true. 时。否则，如果模拟是一个纯弹性或纯声学介质，参考频率的概念就不适用。

3.2 如何使用 Gmsh 生成外部网格

Gmsh 是一个 3D 有限元网格生成器，可用于生成四边形和六面体网格。因此，它是生成可由 SPECFEM2D 处理的网格的良好选择。对于 SPECFEM2D 用户，Gmsh 的两个模块是相关的：几何模块和网格模块。在 EXAMPLES/Gmsh_example 目录中给出了一个示例，展示了如何使用这两个模块生成外部网格。所考虑的模型是一个包含两个由不同材料填充的圆的同质正方形。

步骤：

该操作将生成一个 SqrCirc.msh 文件，该文件必须经过处理以获取 SPECFEM2D 所需的所有文件（参见前一节）。这可以通过运行位于 utils/Gmsh 目录中的 Python 脚本 LibGmsh2Specfem.py 来完成：

此外，还会生成与模型各边界相对应的四个类似于 free_surface_file 的文件。

生成几何结构：可以通过加载文件 SqrCirc.geo 到 Gmsh 中来生成几何结构。在 .geo 文件的末尾有几行代码，用于定义盒子的边界和介质。这是按照以下约定完成的：

Physical Line("Top") = {1};   // 对应盒子的顶部边界
Physical Line("Left") = {2};  // 对应盒子的左边界
Physical Line("Bottom") = {3}; // 对应盒子的底部边界
Physical Line("Right") = {4}; // 对应盒子的右边界
Physical Surface("M1") = {10}; // 外围介质
Physical Surface("M2") = {11,12}; // 两个圆的内部

为不同材料填充圆：如果想为两个圆分别填充不同的材料，需要如下编写代码：

Physical Surface("M1") = {10}; // 外围介质
Physical Surface("M2") = {11}; // 大圆的内部
Physical Surface("M3") = {12}; // 小圆的内部

生成网格：在 Gmsh 中选择适当的选项生成并保存 2D 网格：

细分算法：选择四边形单元(All quads)。

单元顺序：选择 1 表示 4 节点网格，选择 2 表示 9 节点网格。

```
python LibGmsh2Specfem.py SqrCirc -t A -b A -r A -l A
```
脚本选项解释：
- -t：模型顶部
- -b：模型底部
- -r：模型右侧
- -l：模型左侧
- 选项值 A 表示吸收边界，F 表示自由边界（所有边界默认都是吸收边界）。
文件关联：生成的文件与前一节中提到的文件名的关联如下：
- Mesh_SqrCirc 对应于 mesh_file
- Material_SqrCirc 对应于 material_file
- Nodes_SqrCirc 对应于 nodes_coords_file
- Surf_abs_SqrCirc 对应于 absorbing_surface_file
- Surf_free_SqrCirc 对应于 free_surface_file

3.3 如何使用 Cubit/Trelis 生成外部网格

Trelis（以前称为 Cubit）是由 Csimsoft 分发的 2D/3D 有限元网格生成器，可以用于生成四边形和六面体网格。Trelis 提供了一个便捷的 Python 接口（模块 cubit），可通过 Python 脚本创建网格。建议你参考 Csimsoft 网站上的逐步教程以快速入门：http://www.csimsoft.com/tutorials.jsp。Trelis 提供了许多强大的图形工具，虽然非常实用，但这里我们主要关注命令行功能，尤其是 Python 接口，这是 Cubit/Trelis 的真正优势。

快速开始步骤：

打开 Cubit/Trelis：打开软件后，点击相应的图标，选择 Python 文件类型（*.py）并运行以下脚本：
- 对于生成带 PML 层的简单 2D 网格，运行：
  simplest2DexampleWithPmls.py（文件路径：utils/cubit2specfem2d/）
- 如果你需要进行轴对称模拟，运行：
  simpleAxisym2dMesh.py（文件路径：utils/cubit2specfem2d/）
生成网格文件：这些脚本将生成一个带有 PML（完美匹配层）的简单网格。
之后，再次点击图标并运行以下脚本：
cubit2specfem2d.py（文件路径：utils/cubit2specfem2d/）

此脚本将在当前目录中生成所有必要的网格文件，以用于 SPECFEM2D 模拟。
进一步学习：你还可以参考其他带有注释的示例，特别推荐查看 EXAMPLES/paper_axisymmetry_example 文件夹，首先阅读该文件夹中的 README 文档。这些脚本中的注释非常有帮助，建议仔细阅读。
使用 Script 标签：另一种将 Python 与 Cubit/Trelis 结合使用的方法是利用 Script 标签。该标签是一个真实的 Python 终端，可用于通过 cubit 模块将命令行 Python 指令传递给 Cubit/Trelis。

如果命令行面板中没有显示 Script 标签，请依次执行以下操作：
Tools -> Options... -> Layout -> Cubit Layout -> Show script tab

这个标签可以让你直接在 Cubit/Trelis 中逐行执行脚本，从而帮助你理解如何创建网格并将其导出为 SPECFEM2D 格式。

通过这些步骤，你可以快速学习如何使用 Cubit/Trelis 创建网格，并将网格文件导出为 SPECFEM2D 可读格式。

3.3.1 关于 Cubit/Trelis 内置 Python 的注意事项

请注意，Cubit 内置的 Python 与标准 Python 语言有一些令人困扰的差异：

字符串连接问题：使用双引号和单引号连接字符串时，可能会出现问题，即使字符串已经存储。例如：
```
a = "aString"
b = a + 'anotherString'
```
下面这个例子不能正常工作：
```
pathToMeshDir = pathToSpecfem + 'EXAMPLES/paper_axisymmetry_example/MESH'
cubit.cmd('cd "' + pathToMeshDir + '"')
```
冒号后不能加注释：例如，下面的代码不能正常工作：
```
if True: # Just a dummy comment
print "Ok!"
```
删除注释后代码可以正常运行。
os.makedirs("~/aDirectory/") 的问题：该命令会创建一个名为 ~ 的目录。要删除该目录，请执行
```
rm -R ./~
```
，但切记永远不要执行 rm -rf ~，这将删除你的主目录。
不能使用 sys.argv：Cubit 内置 Python 不支持使用 sys.argv。
不能使用多行注释：例如，不能在脚本开头使用 """ """ 形式的多行注释。

此外，可能还有其他限制。编写脚本时要注意这些问题，避免受到 Cubit 内置 Python 版本的限制困扰。

3.4 关于吸收 PML（完美匹配层）的注意事项

如果你使用 CPML（卷积完美匹配层），需要一个额外的文件列出 CPML 元素。文件的第一行是 CPML 元素的总数，接下来每行包含一个 CPML 元素，格式如下：

第一列：CPML 元素编号
第二列：标志（如下表所示）：

表 3.1：CPML 元素标志定义

Flag	Meaning
1	元素仅属于 X 方向的 CPML 层（位于 Xmin 或 Xmax）
2	元素仅属于 Y 方向的 CPML 层（位于 Ymin 或 Ymax）
3	元素同时属于 X 和 Y 方向的 CPML 层（即 CPML 的角落元素）

要了解如何将 PML 层添加到使用外部网格生成器（如 Gmsh）创建的网格或模型中，可以参考 EXAMPLES/CPML_absorbing_layers 目录中的示例。

关于 PML 层中的网格元素

PML 层中的网格元素可以是声学或弹性介质，但不能是粘弹性或孔弹性介质。当定义模型时，你需要将这些吸收元素定义为声学或弹性介质。如果忘记这样做，代码会自动将粘弹性或孔弹性 PML 元素转换为弹性元素。这会导致 PML 层不再完全匹配，因为在进入 PML 层时，物理模型从粘弹性或孔弹性变为弹性介质，但这种改变通常只会导致微小的、可忽略不计的反射。
如果使用 PML 和外部网格（如 Gmsh、CUBIT/TRELIS 等生成的网格），应尽量确保 PML 层内的网格元素尽可能规则，最好是通过“挤压”计算域外边缘的网格元素生成的非变形矩形网格。这样可以使 PML 在时间上更加稳定，因为从数学角度来看，PML 是弱不稳定的，PML 内的严重变形网格更容易触发不稳定性。

自动标记 CPML 元素

如果你有一个以 SPECFEM2D 格式存储的 CUBIT（或类似工具）生成的网格，但不知道如何为其分配 CPML 标志，我们提供了一个小型的 Fortran 程序，可以自动完成这项工作。该程序位于 utils/CPML/convert_external_layers_of_a_given_mesh_to_CPML_layers2D.f90。

使用该脚本创建 PML 层时，你不需要在网格生成器中手动标记这些外部层。但你仍然需要像在吸收边界条件下那样指定网格的边界。脚本会自动提取位于 PML 的元素，并提示你输入 PML 层的厚度。例如，如果你创建了一个 1 米尺寸的区域，脚本提示输入 PML 厚度时可以输入 3.1，生成厚度为 3 个元素的 PML。建议输入略大一些（5-10%）的值，以应对可能的元素偏移，特别是在没有通过 CUBIT/TRELIS 的挤压功能创建 PML 区域的情况下。

稳定 PML 的建议

为了增强 PML 的稳定性，建议添加几何规则的非 PML 过渡层，在该层中也关闭衰减（即设置 Qκ = Qμ = 9999），类似于图 3.4 中的红色过渡层。未来，utils/CPML 目录中的工具将自动实现该过渡层。

更精确的处理方式

如果想使用 PML 层，不需要按照 Python 脚本标记这些层，因为 xmeshfem2d 不识别 CPML 标志。如果开发该脚本的人员修复了这个问题，将对用户大有帮助；目前，这些层不需要物理标识符。
然而，模型的“顶部”、“底部”、“左侧”和“右侧”边界需要重新分配给模型的外边界（例如，左边界为最左侧的 PML，右边界为最右侧的 PML，顶部为最顶层的 PML，底部为底层的边界）。这些边界线需要具备相应的标识符。
如果希望顶部是反射性边界，则不需要创建顶部 PML，因为分配标志的 Fortran 脚本会忽略位于 PML 层厚度范围内的顶部元素。
Fortran 程序 utils/CPML/convert_external_layers_of_a_given_mesh_to_CPML_layers2D.f90 不会创建额外的元素，它仅将边界内的一定距离范围内的元素标记为吸收边界。

注意 PML 中的速度和密度模型

如果使用 PML 且使用外部速度和密度模型（例如，将 MODEL 标志设置为 external），请注意，数学上 PML 不能处理 PML 层内沿法线方向的异质性。这是因为定义的阻尼剖面假定法线方向上的速度和密度模型是常数。因此，你需要修改速度和密度模型，使其在 PML 层内为 1D 模型（如图 3.5 所示）。同样的原则适用于底层，你应确保模型沿垂直方向恒定。

总结

在物理区域内使用 2D 速度和密度模型，而在所有 PML 层中，通过从 PML 内边缘的值进行延续，使速度和密度模型沿法线方向保持恒定。

3.5 控制外部网格的质量

要检查外部生成的网格中各个元素的质量，可以在终端中输入：

./bin/xcheck_quality_external_mesh

在运行该命令时，回答第一个问题时选择“3”。此程序将告诉你整个网格中哪个元素质量最差（最大偏斜度，即元素角度的最大变形）。你可以在使用的网格生成软件中修改这些质量较差的元素，然后重复这一过程，直到整个网格的最大偏斜度小于或等于 0.75（大于 0.75 会有风险，0.75 到 0.80 之间可能还能工作，但如果有单个元素的偏斜度超过 0.80，网格应被改进）。

程序还会显示一个包含 20 个偏斜度等级的直方图，显示有多少元素的偏斜度大于 0.75，以及这些元素占总元素的百分比。要查看此直方图，可以输入以下命令：