A.Ch49-数值模拟
数值模拟是基于理论与实践的数值求解过程，仿真并与现实相接近，用于复现、诠释物理现象，是目前热门的研究领域。
49.1 ParaView
用于后处理器的 ParaView 以其开源、兼容性高等特点，适用于网格、粒子后处理，受到科学家与工程师的青睐。
我常用的有 h5part、VTK ASCII 格式，分别对应粒子和网格后处理，常用过滤器有 SPH、Clip 等。
49.2 自由面条件
\( \eta \) 是自由面方程，\( \phi \) 是速度势，\( P_c \) 是气垫压强。
49.2.1 运动学条件
\[ \frac{\partial \phi}{\partial y} = \frac{\partial \eta }{\partial t} + \frac{\partial \phi}{\partial x} \frac{\partial \eta}{\partial x} + \frac{\partial \phi}{\partial z}\frac{\partial \eta}{\partial z} \]
49.2.2 动力学条件
\[ \frac{\partial \phi}{\partial t} + \frac{|\Delta\phi|^2}{2} + g\eta=0 \]
如果是气垫船，表面压强不为大气压，则气室下方：
\[ \frac{\partial \phi}{\partial t} + \frac{|\Delta\phi|^2}{2} + \frac{P_c}{\rho} + g\eta=0 \]
49.2.3 线性化自由面条件
微幅波假设：当波面起伏是小量时，波面流体质点运动速度也是小量，\(\frac{A}{\lambda} \ll 1 \)。
- 运动学条件：\( \frac{\partial \phi}{\partial y} = \frac{\partial \eta }{\partial t} \)；
- 动力学条件：\( \frac{\partial \phi}{\partial t} + g\eta=0 \)。
可以得到：
\[ g\frac{\partial \phi}{\partial y} + \frac{ \partial ^2\phi}{\partial t^2}=0 \]
49.3 SPH
SPH 光滑粒子流体动力学。
49.3.1 SPH 验证集
《Geo-disaster Modeling and Analysis: An SPH-based Approach》中存在一些很好的验证集。
49.4 气垫船理论研究
49.4.1 流量系数
流量系数的理论公式推导基于伯努利方程（能量方程），流量系数 \(C_d\) 蕴含了局部能量的损耗，取值范围为 \(0 \leq C_d \leq 1\)：
\[ Q = C_d A \sqrt{\frac{2\Delta P}{\rho_a}} \]
49.5 HPC 高性能计算
我能用到的高性能计算/辅助工具有：
- GCC Fortran 编译器；
- Meson 构建工具；
- OpenMP 并行编程；
- OpenBLAS 矩阵库；
- PETSc 线性代数库；
- HDF5 文件格式；
- Fortran_stdlib 标准库；
- Python 绘图脚本；
- Octave 科学工具箱；
- PowerShell 脚本；
- MSC.Patran 建模和网格划分工具；
- ParaView 后处理器。
49.6 CSR 稀疏矩阵格式
CSR（Compressed Sparse Row）稀疏矩阵存储格式是一种常用的稀疏矩阵存储方法，它能够高效地存储稀疏矩阵中的非零元素。CSR 格式主要包括三个组成部分：indices（列索引）、indptr（行指针）以及 data（非零元素值）。下面我将详细介绍 indptr（行指针）的作用。
49.6.1 CSR 格式简介
CSR 格式的主要优点在于能够快速访问矩阵的行，这对于许多基于行的操作来说是非常重要的。CSR 格式由以下三部分组成：
data- 存储所有非零元素的值，按行优先顺序排列。indices- 存储每个非零元素对应的列索引。indptr- 存储每行第一个非零元素在data数组中的位置。对于最后一行，indptr包含的是所有非零元素的总数加一。
49.6.2 行指针 (indptr) 的作用
indptr 数组是一个整数数组，它的长度等于矩阵的行数加一。indptr 的每个元素表示对应行的第一个非零元素在 data 和 indices 数组中的起始位置。具体来说：
indptr[0]总是 0，表示第一行非零元素的起始位置。indptr[row]表示第row行（从 0 开始计数）的第一个非零元素的位置。indptr[row+1] - indptr[row]表示第row行中的非零元素数量。
举个例子，假设有一个 4x4 的稀疏矩阵如下所示：
| |
该矩阵的 CSR 格式表示为：
data:[2, 1, 3, 4]（非零元素值）indices:[1, 3, 2, 0]（列索引）indptr:[0, 2, 3, 4, 4]（行指针）
这里 indptr 的解释如下：
indptr[0] = 0：第一行的第一个非零元素在data中的位置。indptr[1] = 2：第二行的第一个非零元素在data中的位置。indptr[2] = 3：第三行的第一个非零元素在data中的位置。indptr[3] = 4：第四行的第一个非零元素在data中的位置。indptr[4] = 4：最后一行（第 4 行）没有非零元素，因此indptr[4]表示所有非零元素的总数加一。
通过这种方式，我们可以快速访问每一行的非零元素及其对应的列索引。例如，要获取第一行的非零元素，我们可以通过 indptr[0] 和 indptr[1] 访问 data 和 indices 数组：
- 第一行非零元素的值：
data[indptr[0]:indptr[1]] = [2, 1] - 第一行非零元素的列索引：
indices[indptr[0]:indptr[1]] = [1, 3]
49.6.3 总结
indptr 数组是 CSR 格式的一个关键组成部分，它帮助我们快速定位每行的第一个非零元素在 data 和 indices 数组中的位置。这使得在处理稀疏矩阵时能够高效地遍历每一行的非零元素。
CSR 格式不仅在内存使用上节省空间，在许多矩阵操作中也能提供时间上的优势，特别是在处理大型稀疏矩阵时。对于需要频繁访问行或基于行的操作，CSR 格式可以提供高效的性能。 然而，对于需要频繁访问列或基于列的操作，可能需要考虑使用 CSC（Compressed Sparse Column）格式，因为它在这些场景下表现更好。
C 语言的行优先的存储格式适合 CSR，Fortran 语言的列优先的存储格式适合 CSC。
49.7 BLAS 简单指南
当程序运行的很慢时，一定会有至少一个大循环，并且可能存在多个 \( O(n^2) \) 及以上的算法。做科学计算的时候，难免遇到浮点型的矩阵运算，矩阵运算往往满足了前述的 2 种情况， 采用优化的 BLAS/LAPACK 库是一种绝佳选择，直接法可能不是最完美的，但一般是发挥最稳定的。我们常用到以下三种求解例程：
?gemm：矩阵乘法；?getrf/?getrs：矩阵求逆；?getrf/?getrs：线性方程组求解。
它们的时间复杂度一般都在 \( O(n^2) ~ O(n^3) \) 范围内，比如 OpenBLAS/oneMKL 这类被优化和并行的线性代数库，在中小型矩阵中十分有效。 额外地，使用 Intel VTune Profiler 工具查看程序并行情况，尝试为所有耗时部分的 \( O(n^2) \) 及以上的算法代码至少实施 OpenMP。 只要是线性的时间复杂度的必要功能，则不会影响代码运行效率。