CN106250102A - 交错网格有限差分正演模拟优化的方法 - Google Patents

Abstract

提供一种交错网格有限差分正演模拟优化的方法，包括：获取三维空间的模型区域节点的内存资源，并根据高阶交错网格有限差分正演模拟的计算参数确定进程的网格规模，其中每个进程具有独立的内存空间；根据每个进程内的模型区域，以及x、y、z三个方向的网格规模，利用卷帘式分配建立多线程的数据模型；以及每个线程计算自己区域内的数据模型，以及通过同步操作进行每个时间片下每个数据模型区域立方体六个面上数据的共享。通过上述技术方案，每个进程具有独立的内存空间，并在每个时间片下每个数据模型区域立方体六个面上进行数据共享，提高了交错网格有限差分进行三维弹性波正演数值模拟的数据交换及通讯的效率，从而提高了该方法的运行效率。

Description

交错网格有限差分正演模拟优化的方法

技术领域

本发明属于计算机并行领域，具体涉及一种交错网格有限差分正演模拟优化的方法。

背景技术

波动方程数值模拟不仅能保持地震波的运动学特征，而且还能保持地震波的动力学特征，可用于检验地震数据处理的效果，目前有限差分法在波动方程数值模拟中具有较大的优越性及广泛性，交错网格有限差分法综合其在占用内存、计算精度、计算效率和并行算法实现等方面的特点，与有限元法、虚谱法等比较还是具有相当的优势。但是，如图1所示，目前的交错网格有限差分三维弹性波正演数值模拟的数据交换量及通讯量都很大，传统的方法只实现了进程级的并行，达不到内存共享，不能减少通信时间，从而降低了交错网格有限差分的三维波动方程的求解效率。

针对上述技术问题，现有技术中尚无良好解决方案。

发明内容

本发明的目的在于解决上述现有技术中存在的难题，提供一种交错网格有限差分正演模拟优化的方法，通过采用多线程和多进程的并行算法，将数据量和运算量分散到多个进程、线程上。进程内的计算由于共享了内存，可以有效降低边界数据交换的通信时间，提高交错网格有限差分进行三维弹性波正演数值模拟的数据交换及通讯的效率，从而提高交错网格有限差分三维正演模拟的运行效率。

本发明是通过以下技术方案实现的：

本发明提供了一种交错网格有限差分正演模拟优化的方法，包括：获取三维空间的模型区域节点的内存资源，并根据高阶交错网格有限差分正演模拟的计算参数确定进程的网格规模，其中每个所述进程具有独立的内存空间；根据每个所述进程内的模型区域，以及x、y、z三个方向的网格规模，利用卷帘式分配建立多线程的数据模型；以及每个所述线程计算自己区域内的数据模型，以及通过同步操作进行每个时间片下每个数据模型区域立方体六个面上数据的共享。

进一步地，该方法包括：将所述三维空间的模型区域划分成多个子区域，以及针对每个所述子区域分配一个所述进程。

进一步地，每个所述子区域具有至少一个所述节点。

进一步地，所述计算参数包括：x、y、z三个方向中每个方向上的网格点数、进程数以及有限差分的阶数。

进一步地，将x、y、z三个方向中每个方向上的网格点数的积除以x、y、z三个方向中每个方向上的进程数的积所得到的商作为所述进程的网格规模。

进一步地，利用所述卷帘式分配将多个所述线程分配到每个所述进程，其中，依次循环将多个所述线程分配到所述每个进程。

进一步地，每个循环的所述线程的最大数量小于或等于计算所采用的微处理器的核数。

进一步地，每个时间片下，每个所述进程内的各线程并发地计算自己区域内的数据模型，当一个进程内所有计算线程都执行结束时，交换当前时间片下该进程数据模型区域立方体六个面上的边界数据。

进一步地，所述计算参数包括x、y、z三个方向中每个方向上的网格间距，以及该方法还包括根据所述网格间距和波的传播时长以及时间步长间隔确定总的步长数目，以划分时间片。

与现有技术相比，本发明的有益效果是：本发明中每个进程具有独立的内存空间，并在每个时间片下每个数据模型区域立方体六个面上进行数据共享，在一定程度上提高了交错网格有限差分进行三维弹性波正演数值模拟的数据交换及通讯的效率，从而提高了该方法的运行效率，是一种快速有效的交错网格有限差分三维正演模拟方法。

附图说明

附图是用来提供对本发明的进一步理解，并且构成说明书的一部分，与下面的具体实施方式一起用于解释本发明，但并不构成对本发明的限制。在附图中：

图1示出交错网格有限差分三维正演模拟传统算法；

图2示出本发明实施方式提供的交错网格有限差分正演模拟优化的方法流程图；

图3示出示例的三维空间划分示意图；

图4示出数据模型区域立方体数据共享示意图；

图5示出交错网格的网格几何关系示意图；以及

图6示出优化后的交错网格有限差分三维正演模拟方法的流程图。

具体实施方式

以下结合附图对本发明的具体实施方式进行详细说明。应当理解的是，此处所描述的具体实施方式仅用于说明和解释本发明，并不用于限制本发明。

下面结合附图对本发明作进一步详细描述：

本发明提供了一种交错网格有限差分正演模拟优化的方法，如图2所示，包括：S201，获取三维空间的模型区域节点的内存资源，并根据高阶交错网格有限差分正演模拟的计算参数确定进程的网格规模，其中每个所述进程具有独立的内存空间；S202，根据每个所述进程内的模型区域，以及x、y、z三个方向的网格规模，利用卷帘式分配建立多线程的数据模型；以及S203，每个所述线程计算自己区域内的数据模型，以及通过同步操作进行每个时间片下每个数据模型区域立方体六个面上数据的共享。

通过本发明提供的上述技术方案，每个进程具有独立的内存空间，并在每个时间片下每个数据模型区域立方体六个面上进行数据共享，在一定程度上提高了交错网格有限差分进行三维弹性波正演数值模拟的数据交换及通讯的效率，从而提高了该方法的运行效率，是一种快速有效的交错网格有限差分三维正演模拟方法。

在实施方式中，上述方法可以包括：将三维空间的模型区域划分成多个子区域，以及针对每个所述子区域分配一个所述进程。其中，每个所述子区域可以具有至少一个所述节点。将三维空间的模型区域划分后的每个子区域可以被称为基本处理单元，每个基本处理单元可以被分配一个进程。在实施方式中，每个进程具有独立的内存空间，数据(例如，地震数据)在整个计算空间中不重复占用内存，因而提高了内存的利用率。每个进程各自独立地完成计算，最终由多个进程共同实现三维正演模拟。

在实施方式中，计算参数可以包括：x、y、z三个方向中每个方向上的网格点数、进程数以及有限差分的阶数。可以将x、y、z三个方向中每个方向上的网格点数的积除以x、y、z三个方向中每个方向上的进程数的积所得到的商作为所述进程的网格规模(或称为网格数目)。假设输入的总的模型规模为N＝NX*NY*NZ，其中NX，NY，NZ分别为x、y和z方向的网格点数。将模型沿着x、y和z方向分别使用NPROCX、NPROCY和NPROCZ个进程进行划分，则总的进程数为NP＝NPROCX*NPROCY*NPROCZ，而每个子进程中的网格规模为N/NP，即为进程的最大模型区域。

在实施方式中，利用所述卷帘式分配将多个所述线程分配到每个所述进程，其中，依次循环将多个所述线程分配到所述每个进程。每个循环的线程的数量受到提供线程的微处理器的核数的限制，在实施方式中，每个循环的所述线程的最大数量小于或等于计算所采用的微处理器的核数。假设输入的总的模型规模为N＝NX*NY*NZ，其中NX，NY，NZ分别为x、y和z方向的网格点数。NPROCX、NPROCY和NPROCZ为模型沿着x、y和z方向使用的进程数，沿着x方向的网格规模为NX/NPROCX，沿着y方向的网格规模为NY/NPROCY，沿着z方向的网格规模为NZ/NPROCZ，每个进程内的模型区域需要处理的总的网格点数为NGRIDPOINTS＝(NX/NPROCX)*(NY/NPROCY)*(NZ/NPROCZ)。首先获取计算节点的处理器核数信息NCORES，依据CPU核数，在进程内启动多线程的个数NTHREADS＝NCORES。采用卷帘式分配方式即将NGRIDPOINTS中的第一个网格点分配给0号线程，第二个网格点分配给1号线程，第NTHREADS个网格点分配给NTHREADS-1号线程，第NTHREADS+1个网格点重新分配给0号线程，依此类推，网格点数编号INDEX-GRID对应的线程编号INDEX-THREAD等于INDEX-GRID％NTHREADS-1。依据上述方式可以建立起卷帘式分配的多线程的数据模型(％表示求余数)。

在实施方式中，计算参数可以包括x、y、z三个方向中每个方向上的网格间距，以及该方法还包括根据所述网格间距和波的传播时长以及时间步长间隔确定总的步长数目，以划分时间片。每个时间片下，每个所述进程内的各线程并发地计算自己区域内的数据模型，当一个进程内所有计算线程都执行结束时，交换当前时间片下该进程数据模型区域立方体六个面上的边界数据，如图3和图4所示。

图3中，以NPROCX、NPROCY和NPROCZ都为3的模型划分方式为例说明数据模型立方体中如何进行数据交换，图中示出了该模型的三维空间进程划分方式。

在图3所示的三维空间中，对于进程标识为13的进程，有6个面的数据需要进行进程间的交换，设六个面分别称为顶、底、左、右、前、后。图4中示出了图3中进程标识为13的进程与其6个方向上需进行数据交换的邻近通信进程的关系示意图。

图6详细说明了本发明的步骤，具体说明如下：

(1)MPI主进程读取网格数等各种模型计算参数，包括：模型沿着x、y和z方向分别使用的进程数NPROCX、NPROCY、NPROCZ；x、y和z方向的网格点数NX、NY、NZ；x、y和z方向的网格间距DX、DY、DZ；有限差分的阶数(取值为：2，4，6，8，10，12)；需保存的数据类型；源信号子波类型；波的传播时长及时间步长间隔，并广播给每个进程。本发明的输入数据包括：计算参数配置文件，在该文件中需要详细列出上述计算所需的模型计算参数，还有接收点位置信息，源信号类型定义及模型文件定义。

(2)各个节点进程获取和分析节点资源。例如，通过Linux系统的文件读取操作读取Linux下的处理器相关信息文件/proc/cpuinfo和系统中关于当前内存状况信息文件/proc/meminfo这两个虚拟文件信息获取节点的资源信息。根据选择的有限差分的阶数(即有限差分的阶数)从计算参数配置文件中读取，(取值为：2，4，6，8，10，12)。保存的数据类型(质点速度、应力、旋度和散度)。

在举例的计算参数配置文件中对数据类型进行了定义，例如，根据编码类型选择，1为保存质点速度，2为保存应力，3为保存旋度和散度，4为保存速度、应力和旋度、散度能量。

数据交换消耗的内存：每个时间片下一个进程与其邻近进程需要交换的数据量所需最大内存量(单位为MB)满足如下公式：

(FDORDER)*4.0*6.0*(max((NX/NPROCX*NZ/NPROCZ),max((NY/NPROCY*NZ/NPROCZ),(NX/NPROCX*NY/NPROCY))))*sizeof(MPI_FLOAT)*pow(2.0,-20.0)，其中FDORDER为交错网格有限差分系数的阶数(取值为：2，4，6，8，10，12)。

波场快照及地震记录大小：波场快照及地震记录大小占用内存满足公式(单位为MB)：2*nseismograms*ntr*ns*sizeof(float)*pow(2.0,-20.0)，其中nseismograms值与计算参数配置文件有关，只保存质点速度时nseismograms值为3，只保存应力时nseismograms值为1，保存旋度和散度时nseismograms值为2，保存速度、应力和旋度、散度能量时nseismograms值为6；ntr为接收点个数；ns为一道的样点数。

程序状态控制数据量：通过操作系统线程获取程序接口(在举例的程序中是使用pthread库的程序接口pthread_attr_init()获取的)，把当前进程中的每个线程所占用的线程堆栈获取相加的总数，计算每个节点在x，y，z三个方向上可以处理的网格数，将上述所需内存量求和与节点的物理内存量进行比较，如果总和超过物理内存量则进行参数的调整。

(3)把这些数据规约到主进程。例如，利用MPI的MPI_Reduce函数进行规约。主进程根据每个节点进程处理的网格能力。网格能力指主进程根据获取到的节点物理资源及x，y，z三个方向上的处理网格数与对应的计算节点所需内存量的关系，分配给该物理节点资源所能承受的x，y，z三个方向上可以处理的网格数，即总的模型计算需要的内存量。计算需要的内存量即反映该网格能力。广播各个节点被分到的模型区域。其中，对于模型区域的划分，例如，假设输入的总的模型规模为N＝NX*NY*NZ，其中NX，NY，NZ分别为x、y和z方向的网格点数。NPROCX、NPROCY和NPROCZ为模型沿着x、y和z方向使用的进程数，沿着x方向的网格规模为NX/NPROCX，沿着y方向的网格规模为NY/NPROCY，沿着z方向的网格规模为NZ/NPROCZ，每个进程内的模型区域需要处理的总的网格点数为NGRIDPOINTS＝(NX/NPROCX)*(NY/NPROCY)*(NZ/NPROCZ)。

(4)根据模型区域，每个节点进程分配计算所需的内存量。

(5)每个节点进程根据模型区域和获取、分析的节点计算能力，确定计算线程数，并把自己的模型区域再次均衡分配给计算线程，即线程模型区域，启动计算线程。

(6)判断是否还有炮。前面(1)到(5)是对输入计算参数及物理资源的分析，从而构建进程模型区域的过程，从步骤(6)开始的步骤是针对每一炮的操作，这里是在每一炮操作开始前可以判断是否还有炮数据，作为流程判断条件。需要计算，是则继续下一步执行，否则结束。

(7)计算当前炮的位置及获取子波。在实施方式中可以这样获取子波：依据计算参数配置文件中源信号的子波类型编码选择子波类型：1＝ricker子波类型，2＝fumue子波类型，3＝来自外部信号文件，4＝sin³子波类型。

第一种对应的Ricker子波计算公式：

r(τ)＝(1-2τ²)exp(-τ²) 公式(1)

其中， $\mathrm{\τ}=\frac{\mathrm{\π}(t-1.5/f_c-t_d)}{1.0/f_c};$

第二种对应的Fuchs-M¨uller子波计算公式：

f_m(t)＝sin(2π(t-t_d)f_c)-0.5sin(4π(t-t_d)f_c) 公式(2)

其中，如果t∈|t_d,t_d+1/f_c|，则f_m(t)＝0；

第四种对应的sin³子波计算公式：

s3(t)＝0.75πf_csin(π(t+t_d)f_c)³ 公式(3)

其中，如果t∈|t_d,t_d+1/f_c|，则s3(t)＝0。

(8)判断是否当前炮的正演模拟时间结束，即判断每一炮下的时间片循环是否结束。在计算参数配置文件中预先给出波的传播时长T及时间步长间隔DT，T/DT即总的步长数目，当步长数目全部循环完则当前炮的时间片循环结束，从而当前炮的正演模拟时间结束。是则跳到(6)执行，否则继续下一步执行。

(9)计算当前时间片(dt)在当前线程模型区域的质点速度。其中在(1)中，读取网格参数后就可以根据波的传播时长T及时间步长间隔DT得到总的步长数目，从而进行时间片划分。

基于有限差分的弹性波波动方程的理论背景及求解过程如下，参考附图5：

描述和分析地震波通过岩石时的运动情况首要遵守的即是物理学上的动量守恒定律，基于该定律可推导出如下的一阶弹性波的微分方程：

$\mathrm{\ρ}\frac{{\mathrm{dv}}_i}{dt}+v_i\mathrm{\ρ}\▿v_i=-\frac{\∂p_{ij}}{\∂x_j}+f_i---\left(1\right)$

式中：ρ代表介质密度；v_i，x_j，f_i分别代表i方向的质点速度、坐标和体力(i,j∈x,y,z)；p_ij代表第ij个应力(外力)，代表i方向的质点速度梯度。

对于三维同性弹性介质在应力作用下，同应变之间存在如下的线性关系：

p_ij＝λθδ_ij+2με_ij (2)

式中：λ和μ称为拉梅常数，为应变张量，

θ＝ε₁₁+ε₂₂+ε₃₃为体积膨胀系数，

${\mathrm{\δ}}_{ij}=\left\{\begin{array}{ccc}1& for& i=j\\ 0& else& \end{array}\right\}$ 称为Kronecker记号。

为了求解弹性介质方程(1)和(2)，可以把速度、应力、拉梅常数的计算，放在交错网格下的离散的笛卡尔空间坐标系(x＝i×dx,y＝j×dy,z＝k×dz)和离散的时间(t＝n×dt)中，其几何关系如图5所示。其中dx、dy和dz表示在x、y和z方向上的空间网格点之间的距离，dt表示时间间隔；NX、NY、NZ和NT分别表示空间网格点数和时间间隔数。最后，把微分算子使用有限差分算子代替。为了适应交错网格有限差分具有2N(N＝1,2,…)阶精度，的4阶有限差分具有以下向前、向后差分算子：

$\frac{\∂f}{\∂x}{|}_{i+1/2}=\frac{1}{dh}\[\frac{9}{8}(f_{i+1}-f_i)+\frac{1}{24}(f_{i+2}-f_{i-1})\]---\left(3\right)$

$\frac{\∂f}{\∂x}{|}_{i-1/2}=\frac{1}{dh}\[\frac{9}{8}(f_i-f_{i-1})+\frac{1}{24}(f_{i+1}-f_{i-2})\]---\left(4\right)$

具有以下中心差分算子^[8]：

$\frac{\∂f}{\∂x}{|}_{i-1/2}=\frac{1}{dh}\[\frac{9}{8}(f_i-f_{i-1})+\frac{1}{24}(f_{i+1}-f_{i-2})\]---\left(5\right)$

把上述方程(2)～(5)代入方程(1)，即可解出一阶弹性波微分方程。

(10)等待进程内所有线程计算结束，以及把进程的模型区域的边界质点速度，使用MPI-3.0的邻近集合通信在新构建的三维笛卡尔拓扑结构中进行数据交换。首先依据模型沿着x、y和z方向分配的进程数NPROCX、NPROCY和NPROCZ，利用MPI的创建三维笛卡尔拓扑结构的函数进行拓扑结构的创建，再利用MPI-3.0的邻近集合通信函数一次完成当前进程与其他相邻6个进程的交界面数据的交换。

(11)计算当前时间片(dt)在当前线程模型区域的应力张量。方程及求解如上第(9)步所述。

(12)等待进程内所有线程计算结束，以及把进程的模型区域的边界应力张量，使用邻近集合通信在新构建的三维笛卡尔拓扑结构中进行数据交换。例如，首先依据模型沿着x、y和z方向分配的进程数NPROCX、NPROCY和NPROCZ，利用MPI的创建三维笛卡尔拓扑结构的函数进行拓扑结构的创建，再利用MPI-3.0的邻近集合通信函数一次完成当前进程与其他相邻6个进程的交界面数据的交换。

(13)根据需要，在实施方式中还可以存储波场快照和地震记录等相关数据然后跳转到(6)或(8)继续执行。存储数据可以包括：波场快照数据，x、y、z三个方向质点速度的波场快照及地震记录数据，应力张量的波场快照及地震记录数据，旋度的波场快照及地震记录数据，散度的波场快照及地震记录数据。

(14)计算结束，释放各节点进程的资源。

本发明描述了一种利用计算机的并行优化技术来达到提高计算效率的方法，通过采用多线程和多进程的并行技术，对数据域进行了进程内的线程并行计算，在一定程度上优化了交错网格有限差分正演模拟方法的并行效率，提高了三维地震波动方程数值求解的计算效率。

上述技术方案只是本发明的一种实施方式，对于本领域内的技术人员而言，在本发明公开了应用方法和原理的基础上，很容易做出各种类型的改进或变形，而不仅限于本发明上述具体实施方式所描述的方法，因此前面描述的方式只是优选的，而并不具有限制性的意义。

以上结合附图详细描述了本发明的优选实施方式，但是，本发明并不限于上述实施方式中的具体细节，在本发明的技术构思范围内，可以对本发明的技术方案进行多种简单变型，这些简单变型均属于本发明的保护范围。

另外需要说明的是，在上述具体实施方式中所描述的各个具体技术特征，在不矛盾的情况下，可以通过任何合适的方式进行组合。为了避免不必要的重复，本发明对各种可能的组合方式不再另行说明。

此外，本发明的各种不同的实施方式之间也可以进行任意组合，只要其不违背本发明的思想，其同样应当视为本发明所公开的内容。

Claims (9)

1.一种交错网格有限差分正演模拟优化的方法，其特征在于，该方法包括：

获取三维空间的模型区域节点的内存资源，并根据高阶交错网格有限差分正演模拟的计算参数确定进程的网格规模，其中每个所述进程具有独立的内存空间；

根据每个所述进程内的模型区域，以及x、y、z三个方向的网格规模，利用卷帘式分配建立多线程的数据模型；以及

每个所述线程计算自己区域内的数据模型，以及通过同步操作进行每个时间片下每个数据模型区域立方体六个面上数据的共享。

2.根据权利要求1所述的方法，其特征在于，该方法包括：将所述三维空间的模型区域划分成多个子区域，以及针对每个所述子区域分配一个所述进程。

3.根据权利要求2所述的方法，其特征在于，每个所述子区域具有至少一个所述节点。

4.根据权利要求1所述的方法，其特征在于，所述计算参数包括：x、y、z三个方向中每个方向上的网格点数、进程数以及有限差分的阶数。

5.根据权利要求4所述的方法，其特征在于，将x、y、z三个方向中每个方向上的网格点数的积除以x、y、z三个方向中每个方向上的进程数的积所得到的商作为所述进程的网格规模。

6.根据权利要求1所述的方法，其特征在于，利用所述卷帘式分配将多个所述线程分配到每个所述进程，其中，依次循环将多个所述线程分配到所述每个进程。

7.根据权利要求6所述的方法，其特征在于，每个循环的所述线程的最大数量小于或等于计算所采用的微处理器的核数。

8.根据权利要求1所述的方法，其特征在于，每个时间片下，每个所述进程内的各线程并发地计算自己区域内的数据模型，当一个进程内所有计算线程都执行结束时，交换当前时间片下该进程数据模型区域立方体六个面上的边界数据。

9.根据权利要求1所述的方法，其特征在于，所述计算参数包括x、y、z三个方向中每个方向上的网格间距，以及该方法还包括根据所述网格间距和波的传播时长以及时间步长间隔确定总的步长数目，以划分时间片。