使用Athena++并行仿真

现在让我们进行更实际的模拟。为了运行更大规模的仿真，需要并行化技术。Athena++提供两种方式并行化：基于OpenMP的共享内存并行化和基于MPI的分布式内存并行化。也可以结合二者实现混合并行化。在本章中，您将学习如何运行并行仿真。我们没有介绍如何为并行仿真配置环境；请咨询您的系统管理员。

MPI并行化

OpenMP并行化是共享内存并行化，通常只在一个计算节点内工作。对于更大的模拟，需要分布式内存并行化。消息传递接口（MPI）正是对应的标准。要使用MPI，您需要在系统上安装MPI。大多数的大学集群系统或超级计算机应该都拥有它，但有时你需要加载模块。请参阅系统文档。MPI的安装不在本教程范围内。

MPI有很多变体，但大多数都可以使用。这里我们假设MPI C++编译器的名称是mpicxx，但请相应地替换它。

这里我们还解释了如何使用HDF5输出，这对于并行仿真非常方便。

1a. 使用-mpi选项配置和编译代码

> cd ~/athena
> python configure.py --prob blast -b --flux hlld -mpi
Your Athena++ distribution has now been configured with the following options:
  Problem generator:       blast
  Coordinate system:       cartesian
  Equation of state:       adiabatic
  Riemann solver:          hlld
  Reconstruction method:   plm
  Hydro integrator:        vl2
  Magnetic fields:         ON
  Special relativity:      OFF
  General relativity:      OFF
  Frame transformations:   OFF
  Viscosity:               OFF
  Compiler and flags:      mpicxx  -O3
  Debug flags:             OFF
  Linker flags:              -lhdf5
  MPI parallelism:         ON
  OpenMP parallelism:      OFF
  HDF5 Output:             ON
> make clean
> make

1b. 使用MPI和HDF5配置和编译代码（可选但是推荐）

如果您的系统上安装了HDF5库（并行版本通常称为pHDF5），则可以使用HDF5输出。要启用此功能，请使用-hdf5选项配置代码。该代码假定库名称为libhdf5。 如果不是，则需要编辑Makefile文件（请参阅[[编译]]）。

> python configure.py --prob blast -b --flux hlld -mpi -hdf5
> make clean
> make

2. 编辑输入文件以设置域分解。让我们以三维MHD爆炸波测试为例。首先，注释掉num_threads参数（这不是必需的，只是为了避免混淆）。

<mesh>
...
#num_threads = 4         # Number of OpenMP threads per process

然后指定分解单元的大小MeshBlock，根网格大小为：

<mesh>
nx1        = 64         # Number of zones in X1-direction
nx2        = 64         # Number of zones in X2-direction
nx3        = 64         # Number of zones in X3-direction

要将该区域剖分为32³的单元块，请在输入文件中的<mesh>块之后（实际上任何位置都可以）添加以下输出块：

<meshblock>
nx1        = 32         # Number of zones per MeshBlock in X1-direction
nx2        = 32         # Number of zones per MeshBlock in X2-direction
nx3        = 32         # Number of zones per MeshBlock in X3-direction

这产生了总共包含2×2×2=8个MeshBlock的网格。显然，根网格必须可以被MeshBlock的大小整除。每个MPI进程都可以在一个或多个MeshBlock上工作，并且每个进程的MeshBlock数量不一定相等，尽管建议相等。因此，在此配置中，您最多可以使用8个进程。如果要运行更大规模的仿真，可以增加根网格大小。例如，如果要将分辨率加倍并使用64个进程，可以按如下方式修改输入文件：

<mesh>
nx1        = 128        # Number of zones in X1-direction
nx2        = 128        # Number of zones in X2-direction
nx3        = 128        # Number of zones in X3-direction

<meshblock>
nx1        = 32         # Number of zones per MeshBlock in X1-direction
nx2        = 32         # Number of zones per MeshBlock in X2-direction
nx3        = 32         # Number of zones per MeshBlock in X3-direction

由于这是教程的一部分，我们希望快速获得结果，所以这里使用的问题规模相对较小。为了获得更好的性能，建议使用更大的MeshBlock。对于具有高级矢量指令（AVX）SIMD支持的特殊英特尔处理器（如2.5GHz Haswell Xeon），建议每个内核（进程或线程）大约64³个单元，但最佳块大小取决于问题和计算机（以及您的容差）。高效地使用代码和计算机是用户的责任！

（可选）如果使用HDF5选项配置代码，请将输出格式更改为HDF5。

<output1>
file_type  = hdf5       # HDF5 data dump
variable   = prim       # variables to be output
dt         = 0.1        # time increment between outputs

3. 使用MPI运行代码

如何启动MPI程序取决于系统，因此请参阅系统的文档。在计算机集群和超级计算机上，通常使用诸如PBS和Slurm之类的作业调度软件，且需要一个作业脚本来提交作业。请注意，我们没有足够的人力为特定机器上的问题提供支持。要启动MPI程序，通常使用mpiexec命令。根据系统的不同，它可能是mpirun，aprun等。为了使用8个核启动并行仿真，请使用

> mpiexec -n 8 ~/athena/bin/athena -i athinput.blast > log

标准输出被重定向到名为“log”的文件。

4. 链接VTK文件（用于VTK输出）

指定VTK输出时，每次输出都会为各个MeshBlock创建文件。为方便起见，vis/vtk目录中提供了join_vtk++程序。该程序将VTK文件组合成一个大的VTK文件。你可以按如下方式编译它

> gcc -o join_vtk++ join_vtk++.cpp -lm

按如下方式使用它

> join_vtk++ -o Blast.00000.vtk Blast.block0.out1.00000.vtk Blast.block1.out1.00000.vtk ...

5. 可视化在最后一节进行说明

OpenMP并行化

OpenMP是基于指令的并行化，它可以在计算节点内并行化。许多编译器都支持它，无需另外安装软件。让我们像上一节一样模拟三维MHD爆炸波。

1. 使用-omp选项配置与编译代码

> cd ~/athena
> python configure.py --prob blast -b --flux hlld -omp
Your Athena++ distribution has now been configured with the following options:
  Problem generator:       blast
  Coordinate system:       cartesian
  Equation of state:       adiabatic
  Riemann solver:          hlld
  Reconstruction method:   plm
  Hydro integrator:        vl2
  Magnetic fields:         ON
  Special relativity:      OFF
  General relativity:      OFF
  Frame transformations:   OFF
  Viscosity:               OFF
  Compiler and flags:      g++  -O3 -fopenmp
  Debug flags:             OFF
  Linker flags:
  MPI parallelism:         OFF
  OpenMP parallelism:      ON
  HDF5 Output:             OFF

2. 移动到工作目录并复制示例输入文件

> cd ~/work
> cp ~/athena/inputs/mhd/athinput.blast .

3. 编辑输入文件来设置域分割和每个进程的OpenMP线程数

我们对采用与MPI案例相同的域分解方式进行OpenMP并行化。您必须将计算域分解成名为MeshBlocks的较小单元。

<mesh>
nx1        = 64         # Number of zones in X1-direction
...
nx2        = 64         # Number of zones in X2-direction
...
nx3        = 64         # Number of zones in X3-direction
...
num_threads = 4         # Number of OpenMP threads per process

<meshblock>
nx1        = 64         # Number of zones in X1-direction per MeshBlock
nx2        = 32         # Number of zones in X2-direction per MeshBlock
nx3        = 32         # Number of zones in X3-direction per MeshBlock

在此示例中，根网格大小为64³，可以创建4个包含64x32x32单元的MeshBlock。另外，此示例还需要指定启动4个线程。请注意，MeshBlock的数量必须大于或等于线程数。OpenMP并行化的最佳线程数取决于系统和问题的规模。请注意，此输入文件中的问题规模相对较小。

4. 实际运行仿真

> ~/athena/bin/athena -i athinput.blast

5. 可视化结果

本次运行中，代码为每个MeshBlock的每个输出时间步创建一个VTK文件。您可以使用HDF5输出将数据组合到单个文件中。3D数据的可视化与MPI案例相同，在本页末尾进行了说明。

6. 改变线程数比较性能

在模拟结束时，代码将显示运行的时间和性能。

cpu time used  = 2.86014892578125e+02
zone-cycles/cpu_second = 2.97875406250000e+05

omp wtime used = 7.46549757500034e+01
zone-cycles/omp_wsecond = 1.14120725000000e+06

“omp wtime”是实际的挂钟时间，“cpu time”时间是总CPU时间。 因此，二者之比是OpenMP并行化实现的加速。您需要比较这些数字并进行调整以获得最佳结果。再次强调，每个进程应该使用多少个线程（MPI rank）取决于系统。我们稍后会再次讨论这个问题。

混合仿真（高级用户）

OpenMP和MPI是不同的并行化方式，但在Athena++中，我们可以同时使用两种同样基于MeshBlocks并行化方法。与常用的基于循环的多线程相比，这种粗粒度多线程需要更少的同步，因此可以提供更好的性能。

您可以混合使用OpenMP和MPI。但是，实现性能最大化并非易事。在典型的英特尔处理器（例如2.5GHz Haswell Xeon）上进行MHD模拟，通常每个物理核心应至少有32³个单元，最好是64³个单元。然后，您可以为每个进程启动2个或4个（甚至可能是8个）线程。每个节点的线程总数应与每个节点的物理核心数相匹配（即，不计入超线程中使用的“虚拟/逻辑核”）。例如，如果每个节点有24个物理核心，则每个节点可以运行24个进程（平面MPI），2个线程×12个进程，4个线程×6个进程，甚至24个线程×1个进程（不推荐）。最佳平衡取决于问题和系统。为简单起见，我们建议从平面MPI并行化开始（即仅使用MPI，不使用OpenMI）。

请注意，如果您使用的是Intel Xeon Phi（Knights Landing）或IBM BlueGene/Q，其CPU设计与Intel Xeon不同。因此，您可以且应当为在每个物理核心上使用4（或2）个线程（或MPI进程）。

例如，假设您要在每个节点具有16个核心的集群上运行仿真，并且您希望使用16个节点（256个核心），每个节点使用4个线程×4个进程（总共64个进程），则可以按如下方式完成：

1. 同时使用-omp和-mpi配置和编译代码

> cd ~/athena
> python configure.py --prob blast -b --flux hlld -omp -mpi -hdf5
> make clean
> make

在大型仿真中高度推荐使用HDF5输出。

2. 编辑输入文件

<mesh>
nx1        = 256        # Number of zones in X1-direction
nx2        = 256        # Number of zones in X2-direction
nx3        = 256        # Number of zones in X3-direction
...
num_threads = 4         # Number of OpenMP threads per process

<meshblock>
nx1        = 64         # Number of zones per MeshBlock in X1-direction
nx2        = 32         # Number of zones per MeshBlock in X2-direction
nx3        = 32         # Number of zones per MeshBlock in X3-direction

这将产生总共包含4×8×8=256个MeshBlock的网格。每个进程将拥有4个MeshBlock，每个线程处理一个MeshBlock。

3. 运行程序

您需要告诉系统如何将进程映射到计算节点。这应该在您的系统文档中有所解释。对于OpenMPI，可以使用-map-by ppr：4：node实现（如果有两个套接字，则应该是ppr：2：socket）：

> mpiexec -n 64 --map-by ppr:4:node ~/athena/bin/athena -i athinput.blast > log

三维MHD数据的可视化

三维数据的可视化稍微复杂一些。 这里我们仍然以VisIt为例。

1. 加载数据

如果您使用的是VTK格式，请打开*.vtk数据库。如果您使用的是HDF5格式，请打开*.athdf.xdmf文件，而不是*.athdf文件。VisIt不知道如何直接读取Athena++的HDF5格式，XDMF文件用作接口。

2. 加载表达式文件（仅针对HDF5）

HDF5格式将所有数据存储为标量变量。使用VisIt的表达式，我们可以构造有用的变量。在控制窗口的菜单上，单击“Controls→Expressions”。然后单击“Load”按钮并打开athena/vis/visit/athdf_MHD_primitive.xml，关闭窗口。

3. 添加伪彩色图和切片操作控件

我们首先查看气体压力分布。单击“Add→Pseudocolor→gas_pressure”，然后点击“Draw”。您只能看到一个蓝色框，这是因为只显示了计算域的表面。要看到内部，有许多不同的方式。您可以添加操作控件来调整图形。单击“Operators→Slicing→ThreeSlice”，然后点击“Draw”，则可以展示出在X，Y和Z平面上的切片伪彩色图。您可以在可视化窗口上拖动来旋转绘图。使用时间滑块查看时间演变。下面是t=0.6的压力解使用默认的“hot”颜色表以对数标度显示在128×128×128单元网格上的显示效果。

要更改操作控件的参数，请单击“Pseudocolor”图标旁边的三角形展开操作控件列表。然后双击“ThreeSlice”操作控件更改参数。您可以通过单击“X”图标来删除操作控件。当有多个操作控件时，您可以通过单击右侧的三角形来更改顺序。

4. 速度向量

您可以使用VisIt叠加多层绘图。基于此功能，我们添加速度矢量来查看流动。默认情况下，VisIt将操作控件应用于所有图，这在某些情况下很方便，但现在通过取消选中“Apply operators to all plots”来禁用此功能。 然后添加“Add → Vectors → gas_velocity”。

接下来双击“Vector”图标打开窗口以设置属性。增加向量数量，并通过更改“Glyphs”选项卡中的比例来调整字形大小。而且矢量原点应该设置为中间。

然后返回到控制窗口单击“Draw”，调整参数得到漂亮的图表。

5. 磁场线

另一个重要的分析工具是可视化磁场线。流线图用于此目的，但在VisIt中创建此类图的精确步骤取决于您使用的版本。注意，必须加载上面提到的athdf_MHD_primitive.xml表达式才能使“bfield”变量和“IntegralCurve”操作控件选项可用。

对于VisIt版本2.11.0（2016年8月）及更高版本，可以使用应用于Pseudocolor图的IntegralCurve操作控件绘制磁场线。单击“Add”，然后单击“operators”，“IntegralCurve”，“bfield”添加绘图。打开操作控件列表，然后双击IntegralCurve图标来修改参数。通过矢量场从起始点进行线积分来计算场线。所以首先我们需要指定起点。例如，将源类型设置为“Plane”，然后将“Normal”设置为“1 0 0”（调整这些参数来调整绘图）。然后将X和Y的“Samples”都增加到10.这意味着它将绘制10x10=100条场线。然后将X和Y“Distance”设置为2.0来覆盖整个计算域。最后，将“Integration direction”（从起点开始积分的方向）更改为“Both”。

对于2.11.0之前的VisIt版本，选择“Add → Streamline → bfield”并双击“Streamline”图标。许多设置需要调整。首先，在“Streamlines”选项卡中指定源（一组积分起点）。选择源类型“Plane”，并将“Normal”设置为“1 0 0”。此设置仅适用于此问题，需要针对不同的问题进行更改。然后将X和Y中的采样数增加到10。将X和Y中的距离设置为2.0来覆盖整个域。同时将积分方向更改为“Both”。

下面的菜单截图来自2.11.0之前的VisIt版本。有关VisIt 2.12.1的类似屏幕截图，请参阅该教程的2018年2月版。

注意：将VisIt中的流线积分器应用于具有多个MeshBlock的Athena++输出时，可能会发生已知错误。如果积分器的源位于MeshBlock边界（此处x=0.0），则积分步可能会“陷入两个MeshBlock之间的间隙”。在这种情况下，产生的流线可能不会从源的某一方向或两个方向延伸。有两种可能的解决方法。首先，您可以稍微扰动源点，例如在此示例中为x=0.00001，以便积分器“跳过”MeshBlock边界间隙。或者，通过设置output1/ghost_zones=1将MeshBlock的重影区域的数据添加到输出文件来“填补空白”。

其次，转到“Appearance”选项卡。 将“Data Value”更改为“Speed”来显示幅值，或将其设置为“Solid”来使用单一颜色。您可以选择自己喜欢的颜色表。最后，取消选中“Show seeds”框。不要忘记单击“Apply”按钮。

返回到控制窗口单击Draw”。

解释流线的时间演变需要谨慎。这些流线的运动不会沿着实际磁场线进行时间演变，因为源位置是固定的。换句话说，它显示了各时间步对应于不同场线的流线。

[[继续教程柱坐标与球坐标|柱坐标与球坐标]]