52、数值计算与并行算法在不同问题中的应用研究

数值计算与并行算法在不同问题中的应用研究

1. 浅水波方程观测数据同化问题

在解决浅水波方程（SWE）的边值问题数值解时，采用了C语言的MPI库。开发该问题的并行软件时，由于缺乏特定方法解决特定问题的效率信息，遇到了一些困难。因此进行了相关研究，结果不仅与问题本身有关，更在很大程度上涉及到解决问题的工具。具体比较了MPI标准的两种流行实现的效率，并研究了使用不同内存分配方式时软件的行为。

1.1 问题的微分形式

将水域中长波传播问题表述如下：
设 $\Omega_{R_E}$ 是半径为 $R_E$ 的球面上的给定区域，边界 $\Gamma = \Gamma_1 \cup \Gamma_2$，其中 $\Gamma_1$ 是沿海岸的边界部分，$\Gamma_2 = \Gamma \setminus \Gamma_1$ 是穿过水域的边界部分。分别用 $m_1$ 和 $m_2$ 表示这些边界部分的特征函数。不失一般性，可假设 $\phi = 0$ 和 $\phi = \pi$（极点）不在 $\Omega_{R_E}$ 内。设 $\Omega = {(\lambda, \phi) \in [0, 2\pi] \times (0, \pi) : (R_E, \lambda, \phi) \in \Omega_{R_E}}$，其中 $\lambda$ 表示经度，$\phi$ 表示纬度。对于 $\Omega_{R_E} \times (0, T)$ 中的未知函数 $u = u(t, \lambda, \phi)$，$v = v(t, \lambda, \phi)$ 和 $\xi = \xi(t, \lambda, \phi)$，有脉冲平衡方程和连续性方程：
[
\begin{cases}
\frac{\partial u}{\partial t} = lv + mg\frac{\partial \xi}{\partial \lambda} - R_fu + f_1 \
\frac{\partial v}{\partial t} = -lu + ng\frac{\partial \xi}{\partial \phi} - R_fv + f_2 \
\frac{\partial \xi}{\partial t} = m\left(\frac{\partial}{\partial \lambda}(Hu) + \frac{\partial}{\partial \phi}\left(\frac{n}{m}Hv\right)\right) + f_3
\end{cases}
]
其中，$u$ 和 $v$ 分别是速度向量 $U$ 在 $\lambda$ 和 $\phi$ 方向的分量；$\xi$ 是自由表面相对于未扰动水平的偏差；$H(\lambda, \phi) > 0$ 是点 $(\lambda, \phi)$ 处水域的深度；函数 $R_f = r^ |U|/H$ 考虑了底部摩擦力，$r^ $ 是摩擦系数；$l = -2\omega \cos \phi$ 是科里奥利参数；$m = 1/(R_E \sin \phi)$；$n = 1/R_E$；$g$ 是重力加速度；$f_1 = f_1(t, \lambda, \phi)$，$f_2 = f_2(t, \lambda, \phi)$ 和 $f_3 = f_3(t, \lambda, \phi)$ 是给定的外力函数。

边界条件形式如下：
[
H U_n + \beta m_2\sqrt{gH}\xi = m_2\sqrt{gH}d \quad \text{在} \quad \Gamma \times (0, T)
]
其中，$U_n = U \cdot n$，$n = (n_1, n_{mn2})$ 是边界的外法向量；$\beta \in [0, 1]$ 是给定参数，$d = d(t, \lambda, \phi)$ 是定义在边界 $\Gamma_2$ 上的函数。

初始条件为：
[
u(0, \lambda, \phi) = u_0(\lambda, \phi), \quad v(0, \lambda, \phi) = v_0(\lambda, \phi), \quad \xi(0, \lambda, \phi) = \xi_0(\lambda, \phi)
]

对时间进行离散化，将区间 $[0, T]$ 按点 $0 = t_0 < t_1 < \cdots < t_K = T$ 划分为 $K$ 个子区间，步长 $\tau = T/K$。用向后差分近似时间导数，并考虑每个时间区间 $(t_k, t_{k + 1})$ 上的方程组。通过对前一时间层的摩擦项进行线性化，得到半离散椭圆型系统。

一般情况下，函数 $d$ 未知，为了封闭问题，考虑在边界的一部分 $\Gamma_0$ 上定义的封闭方程，特征函数为 $m_0$：
[
m_0\xi = \xi_{obs}
]
其中 $\xi_{obs} \in L^2(\Gamma_0)$ 是给定函数（例如来自观测数据）。

因此，在每个时间步，微分问题可以表述为观测数据同化问题：假设 $\xi_{obs}$ 在 $\Gamma_0$ 上定义，函数 $d$ 在 $\Gamma_2$ 上未知且在 $\Gamma_1$ 上为零。求 $u$，$v$，$\xi$，$d$ 满足方程组、边界条件和封闭条件。

为解决这个不适定的逆问题，采用了基于最优控制方法和伴随方程理论的方法。考虑两个带正则化的最优控制问题族，以确定数值自由表面水平 $\xi_h$ 和观测自由表面水平 $\xi_{obs}$ 之间相对于某种特殊范数的计算误差的最小值。提出了一种迭代数值方法来恢复边界函数，该方法通过交替数值求解正问题和伴随问题来迭代改进边界函数。

正问题和伴随问题的数值解基于有限元方法。对区域 $\Omega$ 进行一致三角剖分 $T = {\omega_i} {i = 1}^{N {el}}$，使用Bubnov - Galerkin方法对问题进行空间离散化，以三角形有限元上的线性函数作为试验函数和测试函数，得到高维线性代数方程组。为离散模拟推导了先验稳定估计，并证明了内部节点的二阶近似。

1.2 并行实现及潜在加速估计

求解线性代数方程组时，使用了雅可比迭代过程。该过程可以有效地并行化，并且通过选择时间步长 $\tau$ 很容易满足其收敛所需的对角占优条件。

有限元方法实现算法有以下特点：全局刚度矩阵与时间有关，每个时间步都必须重新计算。但在有限元上实现雅可比方法时，无需显式存储全局刚度矩阵，程序中仅计算局部刚度矩阵的元素，且只有其对角元素与时间有关，每个时间步都要重新计算。使用局部刚度矩阵的元素在三角形上组装残差。

基于数据的显式并行化，将原始计算域划分为多个部分重叠的子域。在雅可比迭代步框架内，每个子域的计算相互独立。每次雅可比迭代步后，需要对重叠部分的数据进行调整。采用无阴影线的分解方式，将原始域划分为仅在切割边界处相交的子域。对于子域的每个边界点，仅在该子域内的三角形上部分计算残差。每次雅可比迭代步后交换数据时，需要对与相邻处理器计算的向量 $(u, v, \xi)$ 的三个分量的残差部分进行额外求和。

并行程序使用C语言和MPI库函数实现。在数据分布方案中，所有进程执行相同的计算，但针对不同的子域。除了第一个和最后一个进程外，交换结构是均匀的。每次雅可比迭代步后，一个进程与其所有邻居交换数据，邻居数量由分解方式决定，与参与计算的进程数量无关。

算法的潜在加速比估计为单处理器计算时间 $T_1$ 与 $p$ 个处理器计算时间 $T_p$ 的比值：$S_p = T_1/T_p$。

设 $T_{op}$ 为执行一次算术运算的时间，$T_{comm}$ 为传输一个值的时间。假设 $N_{nd}$ 是计算域的总节点数，$s$ 是每个雅可比迭代步在一个节点上执行的操作数，$k$ 是时间步数，$\nu$ 是每个时间步的平均雅可比迭代步数。则单处理器执行算法的时间估计为：$T_1 \sim k\nu sN_{nd}T_{op}$。

假设分解域时能将所有计算均匀分配给处理器，则 $p$ 个处理器执行算法的时间为：
[
T_p \sim \frac{T_1}{p} + T_{over} + T_{comm}
]
其中 $T_{over}$ 是与域分解相关的额外计算时间，$T_{comm}$ 是交换时间。

从数据分布方案可知，每个雅可比迭代步必须执行以下操作：
1. 用于计算迭代过程停止准则的全局归约操作：$T_{comm}^1 \sim (T_{op} + T_{comm})\log_2 p$；
2. 切割处每个点的部分残差向量值的交换：$T_{comm}^2 = \mu k\nu mN_{bnd}T_{comm}$，其中 $m$ 是切割处一个点必须传输给相邻进程的值的数量，$\mu$ 是参与交换的相邻进程数量；
3. 与相邻处理器计算的向量 $(u, v, \xi)$ 的三个分量的残差部分的额外求和：$T_{over} \sim k\nu gN_{bnd}T_{op}$，其中 $g = 3$。

因此，在非阻塞两点交换情况下，并行算法的潜在加速比估计为：
[
S_p = \frac{1}{\frac{1}{p} + \frac{g}{s}R + \frac{\log_2 p}{sN_{nd}}(1 + \kappa) + \frac{\mu m}{s}R\kappa}
]
其中 $R = N_{bnd}/N_{nd}$ 表征计算域的分解情况，$\kappa = T_{comm}/T_{op}$ 表征通信环境。对于相对精细的网格，在相当大的进程数量范围内，潜在加速比接近线性。加速比的大小由两个参数决定，较小的 $R$ 值可提供合理的加速比，因此在构建复杂计算域的分解时，除了要求每个进程的计算负载相等外，还需为每个进程提供最小的子域边界长度。较小的 $\kappa$ 值也有利于加速比，因此应选择 $\kappa$ 值小的集群架构。

1.3 MPI实现和内存管理策略比较

在ICM SB RUS的集群系统上分析了内存管理策略。比较了MPI的两种流行实现：MPICH2 v.1.2.1p1和OpenMPI v.1.4.1。测试了问题的两种修改版本：静态内存分配和动态内存分配（使用calloc - free）。

静态分配版本在两个软件包中没有明显优势，所有测试配置的运行时间和数据交换时间差异太小，可忽略不计。相反，动态内存分配版本显示出对所使用软件包及其设置的依赖性。实际上，问题行为的差异与动态内存管理的区别有关，而非两个软件包中MPI函数实现的特点。

OpenMPI包使用ptmalloc内存管理器进行动态内存分配，用于控制内存碎片化并提高应用程序性能。有一个管理内存分配/释放策略的设置“mpi leave pinned”，默认开启。MPICH2包没有此设置，但可以通过调用 $mallopt()$ 函数并传递相应参数来管理glibc系统库处理内存分配请求时使用的策略。

不同内存管理策略下程序执行时间与节点数量的关系图显示，管理静态内存的程序版本在运行时间和与理论估计的一致性（曲线平滑度）方面表现最佳。“动态分配 + mallopt”策略（MPICH2包）在运行时间和与理论的一致性方面排名第二，该策略动态分配内存，但向系统库发出不使用mmap机制（内存页面映射）进行内存分配的命令。在应用程序以独占模式使用Linux操作系统资源（启动时一次性分配大内存空间，仅在工作结束时释放）的情况下，内存碎片化不是典型问题，从堆中分配的内存的管理速率可能更高。

OpenMPI包中“mpi leave pinned”关闭时的结果排名第三。动态内存分配且不修正策略的版本在速率和时间消耗曲线的平滑度方面表现最差。通信时间图显示，所有策略的交换时间差异不大，在 $p = 4$ 处有预期的跳跃，因为此时处理器3和4之间的数据传输网络开始运行；在 $p = 8$ 处有不太明显但可解释的跳跃，因为承载进程4 - 7的主机开始通过同一网络与两个外部邻居交换数据。之后外部邻居数量不再增加，交换时间图上呈现轻微平滑增长，这与需要交换离散均匀范数的总值有关，该过程的成本与进程数量有关，但依赖性较弱。

2. 格中最短向量问题的并行GaussSieve算法

格是 $\mathbb{R}^n$ 的离散加法子群，格的元素称为向量，格可以由其基（一组线性无关的向量 $b_i \in \mathbb{R}^n$）表示。最短向量问题（SVP）表述为：给定格的基，输出格中最短的非零元素。

在密码学中，格中的困难问题如最短向量问题是数字签名、加密方案、哈希函数等许多密码原语的基础。因此，研究最短向量问题有助于评估基于格的密码系统的难度。

解决SVP基本上有三种方法：
1. 枚举算法：对指定搜索区域内的所有系数向量进行穷举搜索，在搜索树中使用分支限界策略。
2. 计算Voronoi单元的算法。
3. 概率筛算法：以高概率输出最短向量。GaussSieve算法属于第三类，是目前最快的筛算法，因此选择它作为并行筛算法的基础。

提出了格中最短向量问题的并行GaussSieve算法。对于最多5个并行线程，并行版本几乎呈线性扩展；对于更多线程，效率会降低。该并行GaussSieve算法在多核CPU上实现，其思想也可在不同的并行平台上实现。

下面用mermaid图展示浅水波方程问题的求解流程：

graph TD;
    A[问题定义] --> B[时间离散化];
    B --> C[正问题和伴随问题求解];
    C --> D[迭代改进边界函数];
    D --> E[判断是否收敛];
    E -- 否 --> C;
    E -- 是 --> F[输出结果];

解决问题类型	采用方法	关键特点
浅水波方程观测数据同化问题	有限元方法、雅可比迭代、最优控制和伴随方程理论	并行化、时间相关矩阵计算、内存管理影响大
格中最短向量问题	并行GaussSieve算法	线程数量影响效率

数值计算与并行算法在不同问题中的应用研究

3. 两种问题的综合分析与对比

为了更清晰地理解浅水波方程观测数据同化问题和格中最短向量问题的求解方法，下面对它们进行综合分析与对比。

对比项目	浅水波方程观测数据同化问题	格中最短向量问题
问题本质	偏微分方程的不适定逆问题，涉及观测数据同化	离散数学中的搜索问题，用于评估密码系统安全性
求解方法基础	有限元方法、最优控制和伴随方程理论	GaussSieve算法
并行化方式	雅可比迭代过程并行化，基于数据划分计算域	并行GaussSieve算法，线程并行处理
影响因素	时间步长、计算域分解、内存管理策略	线程数量
加速比特点	相对精细网格在一定进程数范围内接近线性加速	最多5个线程时接近线性扩展，更多线程效率降低

从上述对比可以看出，两个问题虽然都采用了并行算法来提高求解效率，但由于问题本质不同，所采用的具体方法和影响因素也有很大差异。浅水波方程问题更侧重于数值计算和边界条件处理，而格中最短向量问题则更关注搜索算法的优化和并行性。

下面用mermaid图展示两个问题求解过程的对比：

graph LR;
    A[浅水波方程问题] --> B[时间离散化];
    B --> C[正问题和伴随问题求解];
    C --> D[迭代改进边界函数];
    D --> E[判断收敛];
    E -- 是 --> F[输出结果];
    G[格中最短向量问题] --> H[初始化格基];
    H --> I[并行GaussSieve算法];
    I --> J[判断找到最短向量];
    J -- 是 --> K[输出最短向量];

4. 实际应用与未来展望

在实际应用中，浅水波方程观测数据同化问题的研究有助于海洋潮汐模拟、海洋环境监测等领域，能够更准确地预测海洋中的长波传播和水位变化。而格中最短向量问题的研究对于基于格的密码学至关重要，能够评估密码系统的安全性，为密码算法的设计和优化提供理论支持。

未来，对于这两个问题的研究可以从以下几个方面展开：
1. 算法优化 ：继续探索更高效的数值算法和并行算法，进一步提高求解效率和精度。例如，对于浅水波方程问题，可以研究更先进的有限元方法或迭代算法；对于格中最短向量问题，可以改进GaussSieve算法的并行化策略。
2. 硬件适配 ：随着计算机硬件技术的不断发展，如多核CPU、GPU和分布式计算集群的广泛应用，需要将算法更好地适配到不同的硬件平台上，充分发挥硬件的性能优势。
3. 跨领域应用 ：尝试将两个问题的研究成果应用到其他相关领域，实现知识的迁移和拓展。例如，将浅水波方程的数值计算方法应用到其他流体力学问题中，将格中最短向量问题的算法应用到机器学习中的特征选择等问题。

总之，数值计算和并行算法在解决复杂问题中具有重要作用，通过不断的研究和创新，有望为各个领域带来更高效、更准确的解决方案。

未来研究方向	具体内容
算法优化	探索更高效的数值算法和并行算法
硬件适配	将算法适配到不同硬件平台
跨领域应用	实现研究成果在其他领域的迁移

下面用mermaid图展示未来研究方向的关系：

graph LR;
    A[算法优化] --> B[提高求解效率和精度];
    C[硬件适配] --> B;
    D[跨领域应用] --> E[拓展应用领域];
    B --> E;