CN101866381A - 一种基于逐元技术的Lengendre谱元法弹性波传播并行模拟方法 - Google Patents

Abstract

本发明公开一种基于逐元技术的Lengendre谱元法弹性波传播并行模拟方法，该方法基于Legendre展开的谱元法并采用逐元技术，实现主节点控制若干独立计算节点的并行模拟弹性波传播的方法，包含如下步骤：(1)首先，采用Legendre谱元法改写弹性波传播方程；将整个计算区域空间离散后按均衡原则进行各个计算节点的任务分配，各计算节点并行计算各自任务包含的单元质量矩阵和单元刚度矩阵，并存储在各个计算节点的存储器中；(2)然后，在主节点采用时域迭代方法计算整个迭代时间段的全局格点向量，同时所述的各个并行的计算节点采用逐元技术将均分得到的任务的全局矩阵和向量积转化为单元矩阵和向量积。本发明还提出了如何将逐元技术应用在Lengendre谱元法中，为实现并行方法提供基础。

Description

一种基于逐元技术的Lengendre谱元法弹性波传播并行模拟方法

技术领域

本发明涉及弹性波传播模拟领域，具体涉及一种基于逐元技术的Lengendre谱元法弹性波传播并行模拟方法。

背景技术

弹性波的传播模拟在地球物理，工程地震学和计算声学等领域具有重要地位。如何找到一种更精确、更快速、更适宜并行化的数值模拟方法也一直是国、内外相关领域研究者的工作目标和重点。

当前在弹性波方程模拟中普遍使用的数值方法，主要有有限差分法，有限元法和谱方法.它们有各自的优缺点并在不同的领域得到应用.

有限差分法可以被看作是一种配置点方法，它的误差由配置点的数目和Taylor级数的截断误差决定。有限差分法由于编程简捷和计算速度较快被广泛应用于计算地球物理和工程地震中。但对于低阶有限差分法，它需要大量的网格点来保证精度。而对于高阶有限差分法，又很难高效地处理自由界面和复杂结构问题。

有限元方法是一种加权余量法，引进了一组可展开的函数(尝试函数或形函数)作为基函数用于对偏微分方程的解作截断展开。为了保证由有限项基函数截断展开定义的近似函数能最佳地逼近偏微分方程的精确解，测试函数被引进用于验证近似解带来的余量是否达到最小。有限元法适合处理边界和内部复杂结构，在结构分析和瞬态模拟方面获得了广泛的应用，但它在大型的弹性波方程模拟领域一直没有得到很好应用。其原因一是有限元法需要耗用大量计算资源，二是低阶有限元存在频散现象，而通常的高阶有限元存在伪波。

谱方法，它将尝试函数取为无限可微的全局函数(如Fourier级数)，将问题从物理空间变换到谱空间进行处理，具有良好的收敛性。理论上谱方法可以在一个最短波长中取最少的网格点而达到所需要的精度。为了处理通常的边界条件，在空间展开时，正交多项式(Chebychev or Legendre多项式)经常被用来取代Fourier级数。但由于正交多项式的配置点分布很不均匀，给时间步长的选择带来很大的限制。同时，如同有限差分法一样，谱方法不能自如地处理弯曲的自由界面和复杂的结构问题。很多方法如弯曲坐标和区域分解被用于解决这类问题，但代价是计算量的增加。

谱元法由Patera在1984年提出，早期主要应用于流体动力学.它把有限元法和谱方法相结合，兼具了有限元处理边界和结构的灵活性和谱方法的快速收敛特性。在达到相同的精度前提下，可以采用比传统有限元更稀疏的单元划分，减少了计算时间和内存需求。

谱元法的基本方法是在每一个单元上使用谱方法。同时选取以截断的正交多项式表示的基函数，在各个单元上通过利用配置点插值，以提高级数表示的解的收敛速度。其主要步骤是：(1)首先把计算的区域分成许多子域(单元)，每个子域由若干节点(配置点)组成。(2)在每个子域中把近似解表示成截断的正交多项式展开。(3)用Galerkin方法求解正交问题的变分格式，得到全局的近似解。

Seriani等在1991年首次将Chebychev谱元法引入到声波方程的数值模拟中，各国研究者在其后做了大量的相关研究。但在处理大型计算问题时，谱元法仍然受到内存和计算量的限制。

传统的谱元法或有限元法解波动方程时，需要形成全局刚度矩阵和质量矩阵，考虑衰减时还需引进阻尼矩阵，不管是用一维等带宽存贮还是变带宽存贮，都需要耗费大量的存贮空间，尤其在谱元法中，随着阶次的提高，相应的带宽成倍增加。注意到在波动方程的时间迭代过程中，需要做多次矩阵和向量乘积，如S·p_k，这里S和p_k分别为全局刚度(或质量)矩阵和全局向量。考虑到全局矩阵的稀疏特性，全局的矩阵和相量相乘耗费了大量的计算时间。这些都极大地阻碍了谱元法在具体问题中的应用。

另外，传统的谱元法由于存在全局刚度矩阵和质量矩阵，在算法上限制了其并行化的效率。而这对于谱元法在地球物理，工程地震学，计算声学等领域的大规模实际应用是很大的缺陷。

发明内容

本发明的目的在于，为克服针对传统的谱元法在解波动方程时，需要形成全局刚度矩阵和质量矩阵，一方面需要耗费大量的存贮空间，另一方面耗费了大量的计算时间的缺陷；同时为了实现为地球物理，工程地震学，计算声学等领域的大规模高性能计算服务，从而提出一种基于逐元技术的Lengendre谱元法弹性波传播并行模拟方法。

本发明所提出的一种基于逐元技术的Lengendre谱元法弹性波传播并行模拟方法，该方法基于Legendre展开的谱元法并采用逐元技术，实现主节点控制若干独立计算节点的并行模拟弹性波传播的方法，所述的方法包含如下步骤：

(1)首先，采用Legendre谱元法改写弹性波传播方程；将整个计算区域空间离散后按均衡原则进行各个计算节点的任务分配，使各计算节点间达到内存均衡和计算均衡，各计算节点并行计算各自任务包含的单元质量矩阵和单元刚度矩阵，并存储在各个计算节点的存储器中；

(2)然后，在主节点采用时域迭代方法计算整个迭代时间段的全局格点向量，同时所述的各个并行的计算节点采用逐元技术将均分得到的任务的全局矩阵和向量积转化为单元矩阵和向量积，该步骤具体包含如下子步骤：

(2-1)在各计算节点内部，采用逐元技术按扩散原则完成该节点内部的单元矩阵和单元格点向量积，再在该节点内聚合成节点临时全局格点向量，具体步骤为：

首先，单元向量

通过特定的关系从全局向量p_k扩散得到；然后，在各个计算节点上计算单元矩阵向量积，得到单元结果向量

最后，各个计算节点的结果向量

再通过特定的关系聚合成该计算节点的临时全局向量；

(2-2)通过消息传递接口，将各计算节点的临时全局向量收集到主节点，并聚合成主全局向量，并在主节点上完成时间迭代计算；

(2-3)通过消息传递接口，将主节点上时间迭代成的最新时刻的主全局向量U扩散给各个计算节点，重复上述步骤，完成整个迭代时间段的全局格点向量的运算，得到模拟的弹性波传播结果。

上述技术方案，步骤(1)所述的任务分配具体步骤如下：首先采用谱元法对整个计算模型进行网格划分，得到M×M个谱元单元；然后将M×M个谱元单元均分给P个计算节点，并由其独立计算各自所包含单元的单元质量矩阵和单元刚度矩阵，并存储在各个计算节点的存储器中。

上述技术方案，步骤(2-1)所述的扩散原则为各个计算节点将全局格点向量按统一的全局格点编号和单元格点编号对应规则分发给单元格点向量；所述的聚合是指将各单元格点向量按统一的全局格点编号和单元格点编号对应规则汇集相加得到全局格点向量；

其中，计算区域空间离散后得到的谱元单元，每个谱元单元由若干个网格点组成，由于几个谱元单元间有可能共用一些网格点；每个网格点在各自单元内有一个单元格点编号，同时在整个计算区域内有一个唯一统一的全局格点编号且每个网格点的单元格点编号和全局格点编号间有明确的对应关系，所述的扩散和聚合就是根据这些对应关系在全局格点向量和单元格点向量间进行分发和收集。

上述技术方案，所述的主节点进程系统包含：

扩散聚合模块，用于输入整个计算域的参数；按扩散方式将各子域的参数分派给P个子计算节点，按聚合方式收集P个子进程计算的各个单元质量矩阵，最终形成全局质量阵；其中，所述的参数包括：计算节点个数P，全局格点个数，格点编号，计算域弹性参数和时间迭代长度。

时间迭代模块，用于将计算区域中各格点的时间迭代结果扩散给各个计算节点子进程，聚合各子进程计算得到的临时全局向量得到全局向量，再进行时间迭代运算；

和判断存储模块，用于判断主节点进程的时间迭代是否结束，并保存相关文件，得到弹性波传播的模拟结果；

其中，在各个计算节点进行一指定子域的子进程，主节点进行整个计算域的主进程。

上述技术方案，所述的计算节点子进程系统包含：

接受模块，用于接收本子进程负责子域的参数；

单元矩阵生成模块，用于将接受的子域参数形成本子域内各单元的单元质量阵和单元刚度阵，返回主进程单元质量阵；

时间迭代模块，用于接收本子域内各格点的当前时间的计算结果，计算本子域内各单元矩阵向量积，聚合成该计算节点的临时全局向量，并向主进程返回本子域得到的临时全局向量；

和判断模块，用于判断各个计算节点的时间迭代是否结束；

其中，所述的参数包含：该计算节点编号，该计算节点单元个数，格点编号，该计算域弹性参数和时间迭代长度。

本发明的优点在于，首次将逐元技术引入到Lengendre谱元法中，成千倍的缩减了存贮空间和计算量。同时结合Lengendre谱元法能形成对角型全局质量矩阵和逐元技术的优势，创新性的提出了基于逐元技术的Lengendre谱元法并行方法，并将之应用于弹性波的传播模拟中，为地球物理，工程地震学，计算声学等领域的大规模高性能计算服务。

附图说明

图1本发明的一种基于逐元技术的Lengendre谱元法弹性波传播并行模拟方法的主进程流程图；

图2本发明的的一种基于逐元技术的Lengendre谱元法弹性波传播并行模拟方法子进程流程图。

具体实施方式

下面结合附图针对本发明进行详细说明。

首先，阐述逐元技术在Lengendre谱元法中应用的技术方案。

在谱元法中，全局矩阵是单元矩阵按特定的结构聚合而成的。因此从单元的观点来看，S·p_k的结果全局向量可以通过一系列单元层次的计算得到，

$S\·p_k=\left(\underset{e}{\mathrm{\Σ}}{\tilde{S}}^e\right)p_k=\underset{e}{\mathrm{\Σ}}({\tilde{S}}^e\·{\tilde{p}}_k^e)=\underset{e}{\mathrm{\Σ}}{\tilde{v}}_k^e={\mathrm{\ν}}_k,---\left(1\right)$

其中，

和

分别为构成全局矩阵和全局向量的单元矩阵和向量，而

代表聚合.在上述公式中，全局的矩阵向量乘积是在局部单元上完成的。

具体步骤如下：

(1)首先，计算单元矩阵，而单元向量通过全局和局部的关系获得，也就是说，单元向量

通过特定的关系从全局向量p_k扩散得到；

(2)然后，在局部单元上计算单元矩阵向量积，得到局部单元结果向量 ${\tilde{v}}_k^e={\tilde{S}}^e\·{\tilde{p}}_k^e;$

(3)最后，局部的单元结果向量

再通过特定的关系聚合成全局结果向量。就这样，凭借在全局和局部之间正确的扩散和聚合，我们没有使用全局矩阵S就得到了需要的全局向量。

因此，在程序中就没有必要集成并存贮全局矩阵，而只需计算并存贮单元矩阵。单元矩阵是稠密的，这极大地缩减了内存需求和提高了计算效率。

然后，阐述基于逐元技术的Lengendre谱元法并行化技术方案：

由于采用了逐元技术，在基于Legendre展开的谱元法程序中，绝大部分的工作可以局限在局部单元上实施，各单元间很容易获得负载平衡，并可通过聚合和扩散完成整个计算，因此采用逐元技术的谱元程序从结构上就具有优良的并行效率，具备大规模并行计算的特殊基础。

假定我们在拥有P个计算节点的并行集群计算机系统上进行弹性波传播模拟计算，并采用MPI(消息传递接口)作为并行计算开发环境。

进行网格划分后，整个计算模型具有M×M个谱元单元。只需把这M×M个谱元单元均分给P个计算节点，并由其独立计算单元的质量矩阵和单元刚度矩阵，并存储在自己的存储器中。如此在谱元法的空间离散阶段，通过节点的任务分配，达到了内存均衡和计算均衡的目的。

弹性波传播模拟中大量的计算存在于时间离散中的迭代计算中。由于采用了逐元技术，我们将全局矩阵和向量积转化为单元矩阵和向量积。在并行系统中，拟在各个计算节点存储节点临时全局向量，其主要计算步骤如下：

(1)在各计算节点内部，按扩散原则完成该节点内部的单元矩阵和向量积，再在该节点内聚合成节点临时全局向量。

(2)通过MPI，将各节点的节点临时全局向量收集到主节点，并聚合成主全局向量。并在主节点上完成时间迭代计算。由于主全局质量矩阵为对角阵，且在整个迭代计算过程中保持不变，时间迭代计算可以在主节点上轻松完成。

(3)通过MPI，将主节点上迭代完成的最新时刻的主全局向量U广播扩散给各个计算节点，重复上述步骤，完成整个计算。

实施例

假定我们在拥有一个主节点和P个计算节点的并行集群计算机系统上进行弹性波传播模拟计算，并采用MPI(消息传递接口)作为并行计算开发环境。

进行网格划分后，假定整个计算模型具有M×M个谱元单元。只需把这M×M个谱元单元均分给P个计算节点，并由其独立计算单元的质量矩阵和单元刚度矩阵，并存储在自己的存储器中。如此在谱元法的空间离散阶段，通过节点的任务分配，达到了内存均衡和计算均衡的目的。

弹性波传播模拟中大量的计算存在于时间离散中的迭代计算中。由于采用了逐元技术，我们将全局矩阵和向量积转化为单元矩阵和向量积。在并行系统中，拟在各个计算节点存储节点临时全局向量，如附图1所示，其主节点(主进程)的主要步骤如下：

(1)开始

(2)输入整个计算域的参数，包括计算节点个数P和时间迭代长度。

(3)按扩散方式将各子域的参数分派给P个子计算节点(子进程)。

(4)按聚合方式收集P个子进程计算的各个单元质量矩阵，并形成全局质量阵。

(5)开始时间迭代。

(6)扩散计算区域中各节点的时间迭代结果给各个子进程。

(7)聚合各子进程计算得到的节点临时全局向量。

(8)时间迭代运算。

(9)判断迭代是否结束，否转(6)，是转(10)。

(10)保存相关文件。

(11)结束

如附图2所示各子进程的主要步骤如下：

(1)开始

(2)接收本子进程负责子域的参数。

(3)形成本子域内各单元质量阵和刚度阵。

(4)向主进程返回本子域内各单元质量阵。

(5)开始时间迭代。

(6)接收本子域内各节点的当前时间的计算结果。

(7)计算本子域内各单元矩阵向量积。

(8)聚合各单元矩阵向量积，得到本子域内的节点临时全局向量并返回给主进程。

(9)判断迭代是否结束，否转(6)，是转(10)。

(10)结束。

本发明技术效果：

将逐元技术引入到Lengendre谱元法中，可以成千倍的缩减了存贮空间和计算量。以一个二维正方形计算区域Ω为例。

在基于Legendre展开的谱元法中，Ω被离散为M×M个互不重叠的四边形单元。假定Legendre展开的阶次为N，计算区域具有的总节点数为(MN+1)×(MN+1)。这样，全局矩阵S的大小为S_g＝[2×(MN+1)×(MN+1)]²。而假如采用逐元技术，每个单元矩阵的大小为[2×(N+1)×(N+1)]²，M×M个单元的总需求为S_e＝M×M×[2×(N+1)×(N+1)]²。两种方式的内存需求比如下式所示，

$\frac{S_g}{S_e}\≈M^2\×{\left(\frac{N}{N+1}\right)}^4,---\left(2\right)$

假定M＝100，N＝6，则

也就是说，在此例中，通过采用逐元技术，内存需求至少降低了5000倍。

再关注计算需求，在基于Legendre展开的谱元法中，最大的计算负荷是计算全局的矩阵向量乘积S·p_k，这里只关注耗时的乘法运算，在全局方法中需要进行的乘法总次数为T_g＝[2×(MN+1)×(MN+1)]²。采用逐元技术后，一个全局矩阵向量相乘被转换成M×M次单元矩阵向量相乘，它们需要的总的乘法次数T_e＝M×M×[2×(N+1)×(N+1)]²。因此，两种方法的计算量之比大约为

$\frac{T_g}{T_e}\≈M^2\×{\left(\frac{N}{N+1}\right)}^4,---\left(3\right)$

也就是说，通过采用逐元技术，计算需求的减少和内存需求的减少同样巨大。

我们提出的基于逐元技术的Lengendre谱元法并行方法，在采用MPI(消息传递接口)作为并行计算开发环境下运行，其并行效率可以接近100％。

最后所应说明的是，以上实施例仅用以说明本发明的技术方案而非限制。尽管参照实施例对本发明进行了详细说明，本领域的普通技术人员应当理解，对本发明的技术方案进行修改或者等同替换，都不脱离本发明技术方案的精神和范围，其均应涵盖在本发明的权利要求范围当中。

Claims (6)

1.一种基于逐元技术的Lengendre谱元法弹性波传播并行模拟方法，该方法基于Legendre展开的谱元法并采用逐元技术，实现主节点控制若干独立计算节点的并行模拟弹性波传播的方法，所述的方法包含如下步骤：

(1)首先，采用Legendre谱元法改写弹性波传播方程；将整个计算区域空间离散后按均衡原则进行各个计算节点的任务分配，使各计算节点间达到内存均衡和计算均衡，各计算节点并行计算各自任务包含的单元质量矩阵和单元刚度矩阵，并存储在各个计算节点的存储器中；

(2)然后，在主节点采用时域迭代方法计算整个迭代时间段的全局格点向量，同时所述的各个并行的计算节点采用逐元技术将均分得到的任务的全局矩阵和向量积转化为单元矩阵和向量积，该步骤具体包含如下子步骤：

(2-1)在各计算节点内部，采用逐元技术按扩散原则完成该节点内部的单元矩阵和单元格点向量积，再在该节点内聚合成节点临时全局格点向量，具体步骤为：

首先，单元向量通过特定的关系从全局向量p_k扩散得到；然后，在各个计算节点上计算单元矩阵向量积，得到单元结果向量

最后，各个计算节点的结果向量

再通过特定的关系聚合成该计算节点的临时全局向量；

(2-2)通过消息传递接口，将各计算节点的临时全局向量收集到主节点，并聚合成主全局向量，并在主节点上完成时间迭代计算；

(2-3)通过消息传递接口，将主节点上时间迭代成的最新时刻的主全局向量U扩散给各个计算节点，重复上述步骤，完成整个迭代时间段的全局格点向量的运算，得到模拟的弹性波传播结果。

2.根据权利要求1所述的基于逐元技术的Lengendre谱元法弹性波传播并行模拟方法，其特征在于，所述的任务分配具体步骤如下：首先采用谱元法对整个计算模型进行网格划分，得到M×M个谱元单元；然后将M×M个谱元单元均分给P个计算节点，并由其独立计算各自所包含单元的单元质量矩阵和单元刚度矩阵，并存储在各个计算节点的存储器中。

3.根据权利要求1或2所述的基于逐元技术的Lengendre谱元法弹性波传播并行模拟方法，其特征在于，

步骤(2-1)所述的扩散原则为各个计算节点将全局格点向量按统一的全局格点编号和单元格点编号对应规则分发给单元格点向量；所述的聚合是指将各单元格点向量按统一的全局格点编号和单元格点编号对应规则汇集相加得到全局格点向量。

4.根据权利要求1所述的基于逐元技术的Lengendre谱元法弹性波传播并行模拟方法，其特征在于，所述的主节点进程系统包含：

扩散聚合模块，用于输入整个计算域的参数；按扩散方式将各子域的参数分派给P个子计算节点，按聚合方式收集P个子进程计算的各个单元质量矩阵，最终形成全局质量阵；

时间迭代模块，用于将计算区域中各格点的时间迭代结果扩散给各个计算节点子进程，聚合各子进程计算得到的临时全局向量得到全局向量，再进行时间迭代运算；

和判断存储模块，用于判断主节点进程的时间迭代是否结束，并保存相关文件，得到弹性波传播的模拟结果；

其中，在各个计算节点进行一指定子域的子进程，主节点进行整个计算域的主进程。

5.根据权利要求4所述的基于逐元技术的Lengendre谱元法弹性波传播并行模拟方法，其特征在于，所述的参数包括：计算节点个数P，全局格点个数，格点编号，计算域弹性参数和时间迭代长度。

6.根据权利要求1或4所述的基于逐元技术的Lengendre谱元法弹性波传播并行模拟方法，其特征在于，所述的计算节点子进程系统包含：

接受模块，用于接收本子进程负责子域的参数；

单元矩阵生成模块，用于将接受的子域参数形成本子域内各单元的单元质量阵和单元刚度阵，返回主进程单元质量阵；

时间迭代模块，用于接收本子域内各格点的当前时间的计算结果，计算本子域内各单元矩阵向量积，聚合成该计算节点的临时全局向量，并向主进程返回本子域得到的临时全局向量；

和判断模块，用于判断各个计算节点的时间迭代是否结束；

其中，所述的参数包含：该计算节点编号，该计算节点单元个数，格点编号，该计算域弹性参数和时间迭代长度。

Images