30、高性能计算在地球物理与油藏模拟中的优化策略

高性能计算在地球物理与油藏模拟中的优化策略

在科学计算领域，高性能计算对于解决复杂的地球物理和油藏模拟问题至关重要。本文将探讨两个重要方面的内容：一是地球物理代码的效率优化研究，二是考虑非等温性的裂缝 - 多孔油藏流体动力学过程的并行计算。

地球物理代码的效率优化

在地球物理模拟中，弹性波传播代码的优化是提高计算效率的关键。通过对不同多核系统的优化，我们可以显著提升代码的性能。

图形加速器内存优化

• 全局内存高效利用 ：充分利用图形加速器的全局内存可以将程序的工作速度提高数倍。对于新架构，由于其更完善的内存结构，这种优化的效率会有所降低。
• GPU 内存类型优化
  • 三维数组存储 ：程序使用三维数组，这些数组在内存中作为一个连续区域存储。需要合理分配计算任务，使一个 warp 线程处理相邻存储的元素。
  • 避免重复计算 ：为消除每个时间步长中相同计算的重复，GPU 内存中存储的不是原始网格系数，而是其修改后的版本。
  • 常量内存使用 ：将每个时间步长使用的所有单个常量（如差分方案的大小、每个时间步长源函数的预计算值）写入设备的常量内存，虽然这种优化只能带来约 4% 的性能提升。
  • 缓存和内存访问优化 ：通过一级和二级缓存提高 GPU 全局内存的使用效率。使用 CudaMalloc3D 函数分配内存时对齐所有数组，自动扩展数组以确保每行对齐正确。组织访问模式，使一个 warp 的线程尽可能访问相同的缓存行，将 32 个连续的 4 字节元素的 32 个线程的查询合并为一个。对于 Tesla M2090 和 K40，优先选择 L1 缓存（即 48 KB L1 缓存和 16 KB 共享内存）。
  • 共享内存使用 ：将相应区域的每个块复制到具有每个方向一个附加层（翼区）的共享内存中可以实现更高效的情况。然而，对于 Fermi 和 Kepler 架构的 GPU，由于加载到共享内存中的数据复用率较低，使用共享内存的内核执行时间与不使用共享内存的实现相当。只有 Kepler 架构更高效的内存系统在使用共享内存时具有优势。不同块大小的性能提升结果如下表所示：
    | 线程块大小 | Fermi | Kepler | Pascal |
    | — | — | — | — |
    | 16×8×8 | – | 1.2 | 1.3 |
    | 16×4×4 | 1.0 | 1.2 | 1.3 |
    | 32×4×4 | 0.9 | 1.0 | 1.3 |
    | 64×2×2 | 0.8 | 0.6 | 1.3 |
    | 64×4×2 | 0.9 | – | 1.3 |
    | 64×4×4 | – | 0.98 | 1.5 |
    | 128×2×2 | 0.9 | 0.6 | 1.2 |

多核架构性能比较

弹性波传播代码针对具有非常不同架构特征的不同多核系统进行了优化。不同架构的最佳版本（具有最佳循环序列、流数量、调度器类型和块大小、内存模式、线程大小等）的性能（以 GFLOPS 为单位）如下：
- KNL ：在成熟处理器中脱颖而出，是核心最多且能处理最长向量指令的处理器。
- Tesla P100 ：在 GPU 中处于领先地位，基于 Pascal 架构。其实际性能是其他考虑的 GPU 的 3.2 倍以上，比 KNL 的实际性能领先 1.6 倍以上，尽管其峰值性能是 KNL 峰值性能的 3 倍以上。

优化和并行化的一般建议

• CPU 开发优化规则
  1. 对于二维和三维问题代码，算法的内循环应进行向量化，空间外循环应使用 OpenMP 进行并行化。
  2. 为提高内存访问效率，所有中心数组必须在内存中对齐，例如使用 posix_memalign 或 _mm_malloc 函数。
  3. 使用向量化指令标记代码或使用内在函数进行高效向量化，以适应处理器的向量寄存器大小。
  4. 为了在三维代码中实现更好的缓存和负载平衡，最好选择 (zy)x 循环的嵌套序列，其中 z 和 y 循环合并为一个。
  5. 使用 OpenMP 并行化多核系统的有限差分代码时，最好使用具有最大块大小（默认）的静态调度；对于多核系统，使用具有短块大小的引导调度。
  6. 使用 SMT 并行化有限差分代码时，对于 Intel 处理器，每个核心使用 1 个线程更好；对于 IBM 处理器，每个核心使用最大线程数更好。
  7. 使用 KNL 执行有限差分代码时，最好使用扁平内存模式，将所有主要数组放置在 MCDRAM 内存中。
  8. 尽可能选择多核计算加速器。
• GPU 开发优化规则
  1. 尽可能将整个问题加载到 GPU 的内存中。
  2. 块网格的维度必须与问题的维度匹配。
  3. 对于三维有限差分问题，x 分量的块大小最好使用 32 或 64，y 和 z 分量使用 4 或 2。最大可能的块大小不一定更好。
  4. 使用常量内存存储频繁重用的常量。常量越多，这种存储类型越可取。
  5. 对于从 Pascal 架构开始的较新 GPU，在数据复用率较低（内存受限问题）的情况下，使用共享内存是合理的。
  6. 尽可能选择多核计算加速器。

考虑非等温性的裂缝 - 多孔油藏流体动力学过程的并行计算

在油藏开发中，准确模拟裂缝 - 多孔油藏中的流体动力学过程对于提高采收率至关重要。考虑非等温性的情况下，问题变得更加复杂，需要采用并行计算技术来提高计算效率。

研究背景

随着油气生产难度的增加，大部分油气储量属于难采类型，特别是碳酸盐沉积物中的油藏。碳酸盐油藏具有复杂多变的孔隙结构，其开发既具有挑战性又具有潜力。目前，该领域的研究主要集中在构建整个油田的复杂水动力模型，但这种方法对于解决现场监测中的实际生产问题并不总是有效。因此，需要一种快速计算工具，包括模型、有效算法和软件包。

并行算法的必要性

传统的单孔油藏过滤模型已经有了成熟的算法，但对于存在孔隙油藏和裂缝两种介质的情况，问题变得更加复杂。考虑复杂的热质传递过程只能通过数值实验进行，而全尺寸实验会导致巨大的油损失。数值建模面临着大量未知数和小时间步长的困难，因此需要开发和实现并行数值算法，并在大型超级计算机架构上运行，以显著减少问题的求解时间并降低对运行内存的负载。

数学模型

在一个无限制的油藏中，考虑一口开采裂缝 - 多孔型碳酸盐油藏的生产井。在井开始工作时，地层中会立即形成一个压降漏斗，导致流体从油藏流入井底。基于此问题，建立了一个非等温二元介质模型，描述了热质传递分布，该模型基于组分质量守恒定律，考虑了相的相互作用。
- 质量守恒方程
- 油相：$\frac{\partial(\phi_{\alpha}\rho_{\alpha o}S_{\alpha o})}{\partial t}+\nabla(\rho_{\alpha o}\overrightarrow{U} {\alpha o}) + q {\alpha o} = 0$，其中$q_{m o} = -q_{f o} = -\rho_{m o}\sigma\lambda_{m o}(P_f - P_m)$
- 水相：$\frac{\partial(\phi_{\alpha}\rho_{\alpha w}S_{\alpha w})}{\partial t}+\nabla(\rho_{\alpha w}\overrightarrow{U} {\alpha w}) - q {\alpha w} = 0$，其中$q_{m w} = -q_{f w} = -\rho_{m w}\sigma\lambda_{m w}(P_f - P_m)$
- 能量守恒方程
[
\begin{align }
&\frac{\partial}{\partial t}[(\phi_f\rho_{f o}S_{f o}\varepsilon_{f o} + \phi_m\rho_{m o}S_{m o}\varepsilon_{m o} + \phi_f\rho_{f w}S_{f w}\varepsilon_{f w} + \phi_m\rho_{m w}S_{m w}\varepsilon_{m w}) + (1 - \phi_f - \phi_m)\rho_s\varepsilon_s]\
&+ div[\rho_{f o}\varepsilon_{f o}\overrightarrow{U} {f o} + \rho {f w}\varepsilon_{f w}\overrightarrow{U} {f w}] + div[P_f(\overrightarrow{U} {f o} + \overrightarrow{U}_{f w})] + div[\overrightarrow{W}_f + \overrightarrow{W}_m + \overrightarrow{W}_s] = 0
\end{align }
]
其中，$\lambda_{m o} = \frac{k_mk_{r o}(S_{m o})}{\mu_o}$，$\lambda_{m w} = \frac{k_mk_{r w}(S_{m w})}{\mu_w}$。
- 相关参数定义
- $\overrightarrow{W} f = -(\phi_f[S {f w}\eta_{f w} + (1 - S_{f w})\eta_{f o}])\nabla T$
- $\overrightarrow{W} m = -(\phi_m[S {m w}\eta_{m w} + (1 - S_{m w})\eta_{m o}])\nabla T$
- $\overrightarrow{W} s = -[1 - \phi_f - \phi_m]\eta_s\nabla T$
- $\overrightarrow{W} = \overrightarrow{W}_f + \overrightarrow{W}_m + \overrightarrow{W}_s$
- 油和水的流速：$\overrightarrow{U} {\alpha o} = -\frac{k_{\alpha}k_{r o}(S_{\alpha o})}{\mu_o}gradP_{\alpha}$，$\overrightarrow{U} {\alpha w} = -\frac{k {\alpha}k_{r w}(S_{\alpha w})}{\mu_w}gradP_{\alpha}$（$\alpha = f$，$\overrightarrow{U} {m o} = \overrightarrow{U} {m w} = 0$）
- 密度变化：$\rho_{\alpha w} = \rho_{\alpha w0} + C_w(P_{\alpha} - P_0)$，$\rho_{\alpha o} = \rho_{\alpha o0} + C_o(P_{\alpha} - P_0)$

初始条件为孔隙部分和天然裂缝在初始时刻的原始油藏压力和温度。右边界条件由边界处的恒定压力和温度设定，左边界条件由井底的恒定井底压力和温度设定，且与初始油藏条件不同。考虑到井不受其他作业井的影响，得到的初始边界值问题是一个复杂的准线性混合类型的数学物理方程组。为了解决该系统，使用有限体积法。在某些情况下，线性化系统不再具有自伴性，因此对于初始方程，使用物理过程分裂方法。

并行计算实现

为了求解上述数学模型，提出了一种基于矩阵扫描法的并行算法。计算代数问题的扫描系数需要大量的计算成本，因此将矩阵扫描与方程的系数（矩阵和向量）进行并行化，以加快计算速度。软件实现使用并行运行方法和 MPI 标准，考虑到不仅要在混合计算集群上运行，还要在可能没有 CUDA 的个人计算机上运行。根据计算结果，分析了并行化的效率，并确定了根据使用的网格大小求解方程组所需的最佳进程数。

综上所述，通过对地球物理代码的优化和油藏模拟中的并行计算技术的应用，可以显著提高计算效率，为解决复杂的科学计算问题提供有力支持。在未来的研究中，可以进一步探索更高效的优化策略和并行算法，以满足不断增长的计算需求。

高性能计算在地球物理与油藏模拟中的优化策略（续）

优化策略的技术细节与操作步骤

地球物理代码优化操作步骤

在地球物理代码优化过程中，各项优化措施都有其具体的操作步骤，下面为大家详细介绍：
- 图形加速器全局内存优化
1. 评估全局内存使用情况 ：对程序进行分析，确定哪些部分频繁访问全局内存，以及内存访问的模式。
2. 调整算法逻辑 ：根据图形加速器的内存结构，重新组织算法，确保数据的访问更加连续和高效。例如，将数据存储为连续的块，避免随机访问。
3. 测试与优化 ：对优化后的程序进行测试，通过性能分析工具监测全局内存的使用效率，根据测试结果进一步调整算法。
- GPU 内存类型优化
- 三维数组存储优化
1. 定义数组结构 ：在代码中使用三维数组，并确保其在内存中作为一个连续区域存储。
2. 分配内存 ：使用合适的内存分配函数，如 CudaMalloc3D 进行内存分配，确保数组的对齐。
3. 计算任务分配 ：编写代码逻辑，使得一个 warp 线程处理相邻存储的元素。例如，可以通过循环索引的方式进行任务分配。
- 常量内存使用
1. 确定常量 ：找出每个时间步长中使用的所有单个常量，如差分方案的大小、源函数的预计算值等。
2. 存储到常量内存 ：使用 CUDA 的常量内存存储函数，将这些常量写入设备的常量内存。
3. 代码修改 ：在程序中修改相应的代码，从常量内存中读取这些常量，而不是从全局内存。
- 共享内存使用
1. 数据复制 ：将相应区域的每个块复制到具有每个方向一个附加层（翼区）的共享内存中。可以使用 CUDA 的内存复制函数实现。
2. 数据重用检查 ：在代码中添加逻辑，检查共享内存中数据的重用情况。如果重用率较低，考虑调整块大小或数据布局。
3. 性能测试 ：对使用共享内存的实现进行性能测试，与不使用共享内存的实现进行对比，根据结果进行优化。

油藏模拟并行计算操作步骤

在油藏模拟的并行计算中，主要涉及到数学模型的求解和并行算法的实现，以下是具体步骤：
- 数学模型建立
1. 问题分析 ：明确油藏的物理特性，如孔隙度、渗透率、流体性质等，以及边界条件和初始条件。
2. 方程推导 ：根据质量守恒定律、能量守恒定律等物理原理，推导出描述油藏中热质传递过程的数学方程，如上述的质量守恒方程和能量守恒方程。
3. 参数定义 ：定义方程中的各种参数，如密度、速度、能量等，并根据实际情况进行赋值。
- 并行算法实现
1. 算法选择 ：选择合适的并行算法，如矩阵扫描法。分析算法的复杂度和并行性，确定是否适合当前的问题。
2. 代码编写 ：使用并行编程标准，如 MPI，编写并行算法的代码。将矩阵扫描与方程的系数（矩阵和向量）进行并行化，实现计算任务的分配和通信。
3. 调试与优化 ：对并行代码进行调试，确保其正确性。使用性能分析工具监测并行计算的效率，根据结果调整算法和代码，如调整进程数、块大小等。

优化效果分析

通过对地球物理代码和油藏模拟的优化，我们可以从多个方面分析其效果。

地球物理代码优化效果

• 性能提升 ：从不同架构的性能比较中可以看出，经过优化后，弹性波传播代码的性能得到了显著提升。例如，Tesla P100 的实际性能是其他考虑的 GPU 的 3.2 倍以上，比 KNL 的实际性能领先 1.6 倍以上。
• 内存使用效率 ：通过内存优化措施，如全局内存高效利用、常量内存使用和共享内存使用，减少了内存访问次数，提高了内存使用效率，从而降低了对运行内存的负载。
• 代码可维护性 ：优化过程中对代码进行了重新组织和优化，使得代码的结构更加清晰，逻辑更加合理，提高了代码的可维护性。

油藏模拟并行计算效果

• 计算时间缩短 ：通过并行计算技术，将矩阵扫描与方程的系数进行并行化，显著减少了问题的求解时间。根据计算结果，分析了并行化的效率，并确定了根据使用的网格大小求解方程组所需的最佳进程数，进一步提高了计算效率。
• 资源利用率提高 ：并行计算使得多个处理器核心可以同时参与计算，提高了计算资源的利用率，充分发挥了计算机的性能。

未来研究方向

虽然我们已经取得了一定的成果，但在高性能计算在地球物理与油藏模拟中的应用方面，仍然有许多值得深入研究的方向。
- 更高效的优化策略 ：继续探索新的优化策略，如针对不同架构的自适应优化策略，根据硬件的特点自动调整代码的优化方式，以进一步提高计算效率。
- 新型并行算法 ：研究和开发新型的并行算法，如基于深度学习的并行算法，利用深度学习的强大计算能力和自适应能力，提高油藏模拟的精度和效率。
- 多物理场耦合模拟 ：考虑更多的物理因素，如化学反应、岩石变形等，实现多物理场耦合模拟，以更准确地描述地球物理和油藏中的实际过程。

综上所述，高性能计算在地球物理与油藏模拟中的优化策略具有重要的意义。通过对代码的优化和并行计算技术的应用，我们可以显著提高计算效率，为解决复杂的科学计算问题提供有力支持。未来，我们需要不断探索和创新，以满足不断增长的计算需求。

以下是一个简单的 mermaid 流程图，展示地球物理代码优化的主要步骤：

graph LR
    A[评估全局内存使用情况] --> B[调整算法逻辑]
    B --> C[测试与优化]
    D[定义数组结构] --> E[分配内存]
    E --> F[计算任务分配]
    G[确定常量] --> H[存储到常量内存]
    H --> I[代码修改]
    J[数据复制] --> K[数据重用检查]
    K --> L[性能测试]
    C & F & I & L --> M[优化完成]

通过以上的介绍，相信大家对高性能计算在地球物理与油藏模拟中的优化策略有了更深入的了解。在实际应用中，我们可以根据具体的问题和需求，选择合适的优化策略和并行算法，以提高计算效率和精度。