CN107908913A

CN107908913A - 基于并行计算机的地球动力数模算法

Info

Publication number: CN107908913A
Application number: CN201711407088.9A
Authority: CN
Inventors: 董昊; 戴黎明; 李三忠; 马芳芳; 陶建丽; 胡梦颖
Original assignee: Ocean University of China
Current assignee: Ocean University of China
Priority date: 2017-12-22
Filing date: 2017-12-22
Publication date: 2018-04-13
Anticipated expiration: 2037-12-22
Also published as: CN107908913B

Abstract

本发明涉及一种基于并行计算机的地球动力数模算法，属于地质地球动力数值模拟技术领域。其包括如下步骤：预设模型参数信息、读取模型参数信息、地质地球动力数值模型分为人工输入模型和计算机转换模型，判断是否达到迭代完成时间、读取每一时刻的状态信息文件、用状态文件含有的模拟参数绘制参数场图、利用模拟参数计算可能的地球物理地球化学参数图。本发明提高了数值模拟的计算速度，并从初始模型、分辨率两个方面直接提升模拟结果的可信度，还支持数值模拟结果与其他研究成果的综合对比；可以将高精度的大尺度模拟结果与常规小尺度数值模拟、小尺度物理模拟、小尺度自然观察等相结合，减小结合时由于大尺度模拟的精度问题所带来的误差。

Description

基于并行计算机的地球动力数模算法

技术领域

本发明涉及一种基于并行计算机的地球动力数模算法，属于地球动力数值模拟技术领域。

背景技术

地质学与地球动力学所研究的对象之一是地球构造演过过程。由于研究对象的空间尺度大、时间尺度长，所以不能通过观察方法研究目标区域不同历史时期的构造演化过程及其动力学机制，而数值模拟方法能够基于基本的物理、化学方程，利用计算机实验的方式解决此类问题。现有的地球动力数模算法均是针对同构众核计算机、并行工作站机群、高性能工作站、个人电脑等平台。以上几种计算机的中央处理器(CPU)包含一个或数个结构相同的核心，这种核心在逻辑运算上具有一定优势，但在浮点运算上效率略低。从效率与准确度的角度来说，上述平台并非运行地球动力数模软件的最优平台；在以上平台上执行的数模算法存在运行速度过慢、分辨率过低的问题。而异构众核计算机的中央处理器除了包含几个主要核心外还包含类似于图形处理器(GPU)的从核，此计算机更适合进行大量的数值计算。

虽然异构众核计算机平台更适合运行数值模拟程序，但是目前地球动力数模算法领域没有面向异构众核的多级并行数模算法。由于中央处理器的架构以及异构众核计算机的编程环境不同，因此地球动力数模算法也与以往的不同。现有的地球动力数模算法存在以下缺陷：(1)是通过人工手动编写、用简单几何形状描述、计算机程序进行差分的方法制作模拟所需的初始模型的。(2)目前没有建立可用于数模算法的复杂三维模型的方法，这导致数模成果的可信度低、人工编写效率低，且这种数模只能局限的应用于理论探讨与验证。(3)目前算法运算速度较慢，延长了项目周期、降低了工作效率。(4)目前的数模算法不包含将计算出来的数值模拟结果转换成其他领域(如地球物理、地球化学)常用样式的数据或图示的步骤，这在一定程度上降低了模拟结果的应用性与可信度。

发明内容

针对现有技术存在的上述缺陷，本发明提出了一种新的基于并行计算机的地球动力数模算法，提高了数值模拟的计算速度，并从初始模型、分辨率两个方面直接提升模拟结果的可信度，还支持数值模拟结果与其他研究成果的综合对比。

本发明是采用以下的技术方案实现的：本发明所述的基于并行计算机的地球动力数模算法，包括如下步骤：

步骤一：预设模型参数信息：将人工编写好的模拟所需的或依据计算机辅助设计生成的模型参数信息，通过输入或读取两种方式读入异构众核超级计算机内存中；

步骤二：读取模型参数信息：读取并分析输入到内存中的模型参数信息，按照读取的模型参数信息，以较高的分辨率来计算网格节点的数目生成数值的差分网格化的地球动力数值模型，然后规定网格节点的空间位置，并将空间位置以数组的形式存储在计算机内存中；

构建地球动力数值模型，地球动力数值模型分为人工输入模型和多计算数据转换模型两种方式，分为如下情况：

情况一：如果为人工输入模型：人工输入网格节点的分层界线的顶点坐标，以变量和数组的形式确定其空间位置；

情况二：如果为多计算数据转换模型：读入带有几何形态的几何模型，将几何模型分类为成组的网格节点的几何顶点参数，以确定其空间位置；

步骤三：读取物质参数文件：地球动力数值模型的每个相连分层界线都分割出一块区域，给每块区域赋予物性参数文件中的某个对应的物质参数信息；

步骤四：将模型离散成点，将模型参数传递给点：在地球动力数值模型的整个模拟空间范围内，生成大量的粒子点；网格节点与粒子点的空间位置关系，用插值法给粒子点赋予物质参数信息；到此完成生成数值的差分网格化的高精度三维地球动力数值模型，将以上生成的变量与数组按照一定的输出规则存储为文件；

步骤五：判断是否达到迭代完成时间，分为如下情况：

情况一：若达到迭代完成时间，进入步骤六；

情况二：若未达到迭代完成时间：计算机主进程利用MPI或OpenMP接口分配任务给其他进程：利用MPI或OpenMP接口进行异构众核计算机中的主核之间的控制、通信、计算；其他进程帮助主进程生成计算所需的矩阵；其他进程调用OpenAcc加速接口与Pthreads加速接口，进行矩阵运算：对于计算量巨大的语句模块，则调用异构众核计算机平台下的从核OpenACC加速接口与Pthreads加速接口进行加速处理；从核完成矩阵运算：利用从核的计算优势来处理地球动力数值模拟运算中的超大型矩阵运算以及IO操作；运算结果汇总，模型状态更新：通过上述的加速计算方式，计算出当前时刻的模型状态方程组及其对应的系数矩阵和系数向量进而得到解向量，将解向量依据参数类型分配到各个参数数组中，从而计算出模型的温度场、速度场、压力场；保存模型状态信息：得到模型下一时刻的状态信息并存储在硬盘中；

步骤六：读取每一时刻的状态信息文件：将得到的新的状态信息计算再下一时刻的状态信息，直到得到目标时刻的状态信息；

步骤七：用状态文件含有的模拟参数绘制参数场图：利用存储在硬盘中的目标时刻的状态信息，将状态信息转换成其他领域关心的参数，并将以其他领域的图标形式绘制出数据图；

步骤八：利用模拟参数计算非常规的地球物理、地球化学参数图：利用密度可以计算出地震学家所需要的地质体声速分布图，利用温度数据可以计算出地热学所需要的地表热通量图，利用应力应变场数据可以计算出构造学家所需的应力分布图和断层分布图，利用密度、空间位置数据可以计算出地球物理学家所需的重力异常图。

优选地，所述步骤一中，模型参数信息包括模型大小、网格分辨率、迭代初始时间、迭代终止时间。

优选地，所述步骤五中，依据速度场改变模型内部离散点的位置、应力参数，依据温度场、压力场改变模型内部离散点的温度、粘度参数。

优选地，所述步骤五，还包括如下小步：

第一步：令初始时刻t＝0，迭代最大时间为tmax；

第二步：利用MPI或OpenMP库将粒子点数据和网格数据分配给其他进程；

第三步：各个进程使用openACC库或Athread库调用从核计算，将粒子点参数以加权平均的方式传输到网格中；

第四步：已知网格的各项物性参数，通过上述的加速计算方式，计算当前状态所对应的模型状态方程组，将其以对应的系数矩阵和系数向量的形式表示；

第五步：系数矩阵为稀疏矩阵，其存储格式为三元组；将系数矩阵与系数向量作为解线性方程组函数的参数并调用此函数，例如利用广义最小残差法解线性方程组的mgmres函数；

第六步：在函数mgmres运行过程中，函数包含的矩阵乘法、矩阵加减法均通过从核进行加速计算；最终求解出解向量；

第七步：将解向量所包含的速度、压力、温度信息赋值给网格，计算出每个点相应的应力、应变、密度、粘度信息；

第八步：将粒子点分配给各从核，使用插值法与网格数据计算粒子点的速度、压力、温度、应力、应变信息；

第九步：依据粒子点的速度与时间步长△t，移动粒子点并改变其坐标值，从而完成模型状态的更新；

第十步：令t＝t+△t，更新当前迭代时间；

第十一步：将当前时刻的迭代数据存储到计算机的硬盘中。

优选地，所述步骤五的判断方法如下：

如果t>tmax，则继续向下执行；如果条件不成立，则跳转到第九步进行迭代计算，直到迭代时间大于最大迭代时间；

到达最大迭代时间后，令计数器N＝0，令Nnum等于输出的状态文件的个数。

优选地，所述步骤七中的参数场图包括温度场图、速度场图、密度场图。

优选地，所述步骤八中，依据密度计算声速与声阻抗界面、声速体位置，并绘制成图声速图、声阻抗界面图、声速体位置图；依据温度计算地表热流值，并绘制成地表热流值分布图；利用应力应变信息计算断层可能出现的位置并绘制成断层位置图；利用密度、空间位置数据可以计算重力异常图；利用其它信息绘制类似的地球物理、地球化学图件。

优选地，所述步骤八中，将图件保存在计算机的硬盘中，令计数器N＝N+1；如果N大于Nnum则正常结束运行，否则继续读取下一个状态文件。

本发明的有益效果是：本发明所述的基于并行计算机的地球动力数模算法，提高了数值模拟的计算速度，并从初始模型、分辨率两个方面直接提升模拟结果的可信度，还支持数值模拟结果与其他研究成果的综合对比。具体来说，本发明相对以往的数模算法在计算速度、模拟准确度上均有较大突破，在实际工作中可以提升工作效率、缩短项目周期、提升成果可信度。除此之外，使用本发明将高精度的大尺度模拟结果与常规小尺度数值模拟、小尺度物理模拟、小尺度自然观察等相结合，减小结合时由于大尺度模拟的精度问题所带来的误差。

附图说明

图1是本发明的流程示意图。

图2是本发明的实现方式图。

具体实施方式

为了使本发明目的、技术方案更加清楚明白，下面结合实施例，对本发明作进一步详细说明。

实施例一：

如图1所示，本发明所述的基于并行计算机的地球动力数模算法，是一种基于异构众核超级计算机多级并行的高速、密集取样、多数据对比的地球动力三维数值模拟算法；算法步骤如下：

(1)将人工编写好的模拟所需的或依据计算机辅助设计生成的数据文件，通过输入或读取两种方式读入计算机内存中。

(2)读取并分析输入到内存中的数据文件，依据数据文件所携带的简单物质信息、几何形态描述来生成数值的差分网格化的高精度三维地球动力数值模型，并按照一定的输出规则将模拟的必要信息转化成数值模拟计算所需的特殊文件格式。将此模型作为初始模型存储到存储介质(如硬盘、U盘)中。

(3)读取初始地球动力模型文件到异构众核超级计算机中，将模型以数值的方式存储。

(4)利用MPI或OpenMP接口进行异构众核计算机中的主核之间的控制、通信、计算。

(5)对于计算量巨大的语句模块，如矩阵计算、数组间数据传递，则调用异构众核计算机平台下的从核OpenACC*加速接口与Pthreads加速接口进行特殊加速处理。利用从核的计算优势来处理地球动力数值模拟运算中的超大型矩阵运算以及IO操作。

(6)通过上述的加速计算方式，计算出当前时刻的模型状态方程组及其对应的系数矩阵和系数向量进而得到解向量。

(7)将解向量依据参数类型分配到各个参数数组中，从而计算出模型的温度场、速度场、压力场；依据速度场改变模型内部离散点的位置、应力等参数，依据温度场、压力场改变模型内部离散点的温度、粘度等参数，从而得到模型下一时刻的状态信息并存储在硬盘中。

(8)跳转到步骤(6)，继续用步骤(7)中得到的新的状态信息计算再下一时刻的状态信息，直到得到目标时刻的状态信息。

(9)利用存储在硬盘中的目标时刻的状态信息，将状态信息转换成其他领域关心的参数，并将以其他领域的图标形式绘制出数据图；如利用密度可以计算出地震学家所需要的地质体声速分布图，利用温度数据可以计算出地热学所需要的地表热通量图，利用应力应变场数据可以计算出构造学家所需的应力分布图和断层分布图，利用密度、空间位置数据可以计算出地球物理学家所需的重力异常图。

算法的具体实施过程可如下所述：

(1)读取预设的数值模拟模型参数，参数主要包括模型大小、网格分辨率、迭代初始时间、迭代终止时间等参数。

(2)按照读取的模型参数信息，以较高的分辨率(节点密度)来计算网格节点的数目，然后规定网格节点的空间位置，并将空间位置以数组的形式存储在计算机内存中，。

(3)依据模型内容，选择人工输入模型参数或多计算数据转换模型文件。

(4)如果选择人工输入模型，则人工输入分层界线的顶点坐标，以变量和数组的形式存储在计算机内存中。

(5)读取硬盘中的物性参数文件，每个相连分层界线都分割出一块区域，给每块区域赋予以及读取的物性参数文件中的某个对应的物质参数信息。

(6)在整个模拟空间范围内，生成大量的粒子点，将粒子点的位置信息存储在计算机内存中。

(7)依据网格节点与粒子点的空间位置关系，用插值法给粒子点赋予物质参数信息。

(8)到此完成生成数值的差分网格化的高精度三维地球动力数值模型，将以上生成的变量与数组按照一定的输出规则存储为文件。即将此模型作为初始模型存储到存储介质(如硬盘、U盘)中。

(8)令初始时刻t＝0，迭代最大时间为tmax。

(9)利用MPI或OpenMP库将粒子点数据和网格数据分配给其他进程。

(10)各个进程使用openACC*库或Athread库调用从核计算，将粒子点参数以加权平均的方式传输到网格中。

(11)已知网格的各项物性参数，通过上述的加速计算方式，计算当前状态所对应的模型状态方程组，将其以对应的系数矩阵和系数向量的形式表示。

(12)系数矩阵为稀疏矩阵，其存储格式为三元组；将系数矩阵与系数向量作为解线性方程组函数的参数并调用此函数，例如利用广义最小残差法解线性方程组的mgmres函数。

(13)在函数mgmres运行过程中，函数包含的矩阵乘法、矩阵加减法均通过从核进行加速计算；最终求解出解向量。

(14)将解向量所包含的速度、压力、温度信息赋值给网格，计算出每个点相应的应力、应变、密度、粘度等信息。

(15)将粒子点分配给各从核，使用插值法与网格数据计算粒子点的速度、压力、温度、应力、应变等信息。

(16)依据粒子点的速度与时间步长△t，移动粒子点并改变其坐标值，从而完成模型状态的更新。

(17)令t＝t+△t，更新当前迭代时间。

(18)将当前时刻的迭代数据存储到计算机的硬盘中。

(19)如果t>tmax，则继续向下执行；如果条件不成立，则跳转到步骤(9)进行迭代计算，直到迭代时间大于最大迭代时间。

(20)到达最大迭代时间后，令计数器N＝0，令Nnum等于输出的状态文件的个数。

(21)读取第N个状态文件，绘制常规的温度场图、速度场图、密度场图等。

(22)依据密度计算声速与声阻抗界面、声速体位置，并绘制成图声速图、声阻抗界面图、声速体位置图。

(23)依据温度计算地表热流值，并绘制成地表热流值分布图。

(24)利用应力应变信息计算断层可能出现的位置并绘制成断层位置图。

(25)利用密度、空间位置数据可以计算重力异常图。

(26)利用其它信息绘制类似的地球物理、地球化学图件。

(27)将图件保存在计算机的硬盘中，令计数器N＝N+1。

(28)如果N大于Nnum则跳转到步骤(29)，否则跳转到步骤(21)继续读取下一个状态文件。

(29)输出运行时间等信息，正常结束运行。

如果在此算法中进一步模拟有机质质点的运动并加入有机质演化的计算模块，则可以进行高精度的地质过程与石油天然气成藏数值模拟；如果将此算法的地形模拟结果导入到沉积与层序地层模拟程序中，则可进行高精度深部地质过程与高精度浅部地层地貌形态耦合过程的数值模拟工作。

实施例二：

如图2所示，本发明所述的基于并行计算机的地球动力数模算法，本发明是一种基于异构众核超级计算机的高速、密集取样、多数据对比的地球动力三维数值模拟算法；算法步骤如下：

(3)读取初始地球动力模型文件到异构众核超级计算机中，将模型

(5)对于计算量巨大的语句模块，如矩阵计算、数组间数据传递，则调用异构众核计算机平台下的从核OpenACC加速接口与Pthreads加速接口进行特殊加速处理。利用从核的计算优势来处理地球动力数值模拟运算中的超大型矩阵运算以及IO操作。

(8)跳转到步骤(6)，继续用步骤(7)中得到的新的状态信息计算再下一时刻的状态信息直到得到目标时刻的状态信息。

(9)利用存储在硬盘中的目标时刻的状态信息，计算出其他领域所需要的数据图；如利用密度可以计算出地震学家所需要的地质体声速分布图，利用温度数据可以计算出地热学所需要的地表热通量图，利用应力应变场数据可以计算出构造学家所需的应力分布图和断层分布图，利用密度、空间位置数据可以计算出地球物理学家所需的重力异常图。

需要特别保护的地方在于：

A.实现通过计算机辅助设计软件设计模型的复杂几何结构，然后通过读取结构文件内容、添加地质物理化学等参数的方式将其转化成模拟的初始模型文件这一算法；

B.针对异构众核计算机及其配套加速库的地球动力数值模拟计算过程的算法；

C.对模拟结果的后处理过程与转化成其他领域(如地球物理、地球化学)的数据或图示的算法。

采用本发明所述的基于并行计算机的地球动力数模算法，这项发明提高了数值模拟的计算速度，并从初始模型、分辨率两个方面直接提升模拟结果的可信度，还支持数值模拟结果与其他研究成果的综合对比。具体来说，本发明相对以往的数模算法在计算速度、模拟准确度上均有较大突破，在实际工作中可以提升工作效率、缩短项目周期、提升成果可信度。

实施例三：

本发明因实施例描述的需要，现在本实施例规定以下术语：

(1)二级加速——将利用异构众核计算机的主核作为任务并行分配的单元和承担主体模拟任务而异构众核计算机的从核作为进一步加速主核工作效率的加速方法称为二级加速。

(2)结果模仿图组——利用温度、速度、密度等数值模拟结果进行进一步计算，进而得出声速与声阻抗界面、声速体位置、地表热流值、应力应变等

信息，进而绘制成图声速图、声阻抗界面图、声速体位置图、地表热流值分布图、断层位置图以及包括其他的但不仅限于类似地球物理、地球化学处理结果的模仿图件，并整合制作出图组；将这种图组成为各学科结果的结果模仿图组。

(3)稀疏矩阵——数值为0的元素远多于非0元素且非0元素没有规律时，称此矩阵为稀疏矩阵；为了提升计算速度，这类矩阵以“三元组”的格式存储在内存中。

(4)三元组——一种存储稀疏矩阵的方法，分别将稀疏矩阵中非零元素的信息——行、列、值分为三个数组，按顺序存储在这三个数组中；按照三元组存储稀疏矩阵，不需要存储大量的0元素。

(5)mgmres函数——指使用“广义最小残差”这一线性代数方法来实现求解线性方程组的函数，是本算法可选的求解线性方程组的函数之一。

本发明主要通过异构众核计算及其加速库来进行地质与地球动力学数值模拟。模拟所需模型可由人工编写的简单参数文件来生成，也可以由模型转换生成。数据后处理除了常规数据处理绘图外，还包括制作结果模仿图组。

以对太平洋板块向菲律宾板块的俯冲模拟这一实际的科学研究为例，本发明可以但不限于以下的具体实施方式来实现：

(1)利用VPN方法登录控制异构众核计算机的外部服务器。

(2)将数值模拟用的含有模拟区域内各个圈层(水圈、大气圈、岩石圈、地幔等)物性参数、模型运行控制信息、俯冲地区边界条件信息等的模型参数文件“parameter.dat”，由第三方软件制作的表示俯冲区域内各个圈层精细的空间几何形态的ASCII格式的几何模型文件“STL_model.stl”、初始模型转换程序“model_marker.c”、迭代计算与显示程序“main.c”等文件上传到外部服务器的存储空间中。

(3)将程序与数据通过外部服务器上传到异构众核计算机的工作队列中，编译并运行此程序。

(4)主程序开始初始化进程，调用初始模型转换程序并准备转换几何模型文件。

(5)使用标准C中的scanf函数以十进制的方式读取ASCII格式的几何模型文件“STL_model.stl”，将其中的地质体信息以成组的空间坐标的形式存储在数组中。

(6)读取“parameter.dat”文件中预设好的数值模拟模型参数，参数主要包括模型大小、网格分辨率、迭代初始时间、迭代终止时间、边界条件、各个圈层物理性质等参数。

(7)按照读取的模型参数信息，以0.1km×0.1km×0.1km的分辨率(即节点密度)来计算网格节点总数并进一步确定网格节点的空间位置，并将空间位置以数组的形式存储在计算机内存中。

(8)将读取好的各个圈层物理性质等参数，与几何模型的空间坐标一一对应；即按照位置每一组几何体空间坐标对应一种西太平洋研究区域内的地质体，只记录对应的索引值便可将几何模型与物理性质相对应。

(9)在整个模拟空间范围内，生成大量的粒子点，粒子点的密度以每个网格格子内大于千个粒子点为佳，将粒子点的位置信息存储在计算机内存中。

(10)遍历所有粒子点，依据网格节点与粒子点的空间位置关系，将粒子点周围的网格节点用插值的方法给粒子点赋予物质参数信息。

(11)到此完成生成数值的差分网格化的高精度西太平洋俯冲带地球动力数值模型，将以上生成的变量与数组按照一定的存储顺序存储为初始状态文件“0.dat”。即将此模型作为初始模型存储到外部服务器或超级计算机的硬盘中。

(12)令初始时刻t＝0，迭代最大时间为tmax。

(13)第0号进程为主进程，主进程负责利用MPI库将粒子点数据和网格数据分配给其他进程。

(14)各个进程使用openACC*库调用从核进行循环计算：遍历所有粒子点，判断粒子点位置；将某一网格节点周围的粒子点参数及其权重做求和规约，再做加权平均除法，最后将计算出的节点数值利用规约将数据传输到网格数组中。

(15)已知网格的各项物性参数，通过上述的加速计算方式，计算当前状态所对应的模型状态方程组，将其以对应的系数矩阵和系数向量的形式表示。

(16)系数矩阵为稀疏矩阵，其存储格式为三元组；将系数矩阵与系数向量作为解线性方程组函数的参数并调用此函数，例如利用广义最小残差法解线性方程组的mgmres函数。

(17)广义最小残差(GMRES)方法的思路是将求解线性方程组问题被转化为数值的迭代计算。由于系数矩阵是稀疏矩阵，mgmres函数的参数数组格式被确定为三元组。在函数mgmres运行过程中，涉及到数组与向量之间的乘法、向量与向量之间的加减、数组与数组之间的加减等；迭代计算是在主进程中进行的，为了加快运行速度主进程采用二级加速。加速的方法可以是：主进程自身执行迭代循环，但将数组运算工作分配给其他各个处理器；其他各个处理器在接到任务后，利用从核进行浮点运算最后规约给主进程。最终，求解出解向量。

(18)解向量的各个元素即为按一定顺序排列的速度、压力、温度信息，将此数据赋值到对应网格节点的数值中，并计算出每个点相应的应力、应变、密度、粘度等信息。

(19)将粒子点平均分配给各从核，由于粒子点之间状态的计算并不受计算顺序的影响故可以进一步采用并行计算手段，如利用openACC*库里的“PARALLEL LOOP”编译指示将任务合并分配给从核；即使用插值法与网格数据计算出粒子点的速度、压力、温度、应力、应变等信息；最后通过MPI规约函数将计算结果规约到一起。

(20)将当前时刻的迭代数据以“i.dat”为文件名，其中i为迭代次数，存储到外部服务器或异构并行计算机自身的硬盘中。

(21)依据粒子点的速度与时间步长△t，移动粒子点。其移动方式可以是：将不同的区域的粒子点位置及其周围节点的速度分发给不同的处理器；处理器使用openACC*库里的“PARALLEL LOOP”编译指示调用从核，利用周围节点计算出每个粒子点速度，并将粒子点的几个方向上分速度分别乘以△t得出位置变化量；最后将这一变化量加到粒子点的原坐标上。

(22)令t＝t+△t，更新当前迭代时间。

(23)如果当前迭代时间t大于目标迭代时间tmax，则继续向下执行；如果条件不成立，则跳转到步骤(13)进行迭代计算，直到迭代时间大于最大迭代时间。

(24)到达最大迭代时间后至此步骤，主要的计算工作已经结束。

(25)令计数器N＝0，令Nnum等于输出的状态文件的个数。

(26)读取第N个状态文件，先绘制西太平洋俯冲带区域内常规的温度场图、速度场图、密度场图、岩性图等常规图件。

(27)依据密度计算声速与声阻抗界面、声速体位置，并绘制成图声速图、声阻抗界面图、声速体位置图。

(28)依据温度计算地表热流值，计算方法可以为：按照一定的分辨率计算出各个水平位置的岩体顶部坐标与从此处垂直向下200米处的坐标，将临近粒子点的温度插值到这两类作用点上，使用下伏点与上覆点的温度值之差除以200米，即可计算出当地的地表热流值。

(29)将地表热流值按照等高线图或色彩图的方式绘制出来，制成地表热流值分布图。

(30)利用应力应变信息计算断层可能出现的位置，判断方式可以以第二应变不变量的变化幅度为指标，最后绘制成断层位置图。

(31)利用密度、空间位置数据计算模型中的自由空间重力异常、布格重力异常等，进一步绘制成各种图件。

(32)利用其它信息绘制类似的有助于了解西太平洋俯冲带的地球物理、地球化学图件。

(33)将以上图件保存在外部服务器的硬盘中，令计数器N＝N+1。

(34)如果计数器N大于Nnum则跳转到步骤(35)，否则跳转到步骤(26)继续读取下一个状态文件。

(35)输出运行时间等信息保存到外部服务器的硬盘中，即在外部服务器的硬盘中创建文件“end_info.dat”并将输出信息写入。

(36)正常结束运行，程序退出并向计算机返还占用的计算资源与存储资源。

综上所述，本实施例利用了异构众核计算机迭代计算出了西太平洋俯冲研究区的各个时间点的状态及其图示。提升分辨率与计算速度的方法是采用主核与从核多级加速，由于计算速度的提高可以在有限的时间内计算出更高精度的模型结果；与常规串行算法或简单并行算法相比，可以大幅度减少科学研究中对数值模拟运算的等待时间，进而提升工作效率；在项目周期固定时提升工作量，在工作量固定时缩短项目周期。此外，由于计算速度的提高，数据量庞大的三维模型将不再被科学研究放弃；使用三维模型可以得到更科学更客观的模拟结果，进而提高可信度。而且，由于模型制作采用了将实际几何模型转换的方式，与常规的用简单参数(如某一深度)来规定地质体界线相比，从模拟可信度来讲有很大的提高。而后处理阶段的结果模仿图组为其他领域的研究人员提供了参考，不但填补了数值模拟结果可信程度的缺陷，还方便其他领域的研究人员直接使用数值模拟的数据成果而非仅仅使用数值模拟研究人员的理论分析结果。

以上所述仅为本发明的较佳实施例而己，并不以本发明为限制，凡在本发明的精神和原则之内所作的均等修改、等同替换和改进等，均应包含在本发明的专利涵盖范围内。

Claims

1.一种基于并行计算机的地球动力数模算法，其特征在于，包括如下步骤：

步骤五：判断是否达到迭代完成时间，分为如下情况：

情况一：若达到迭代完成时间，进入步骤六；

2.根据权利要求1所述的基于并行计算机的地球动力数模算法，其特征在于，所述步骤一中，模型参数信息包括模型大小、网格分辨率、迭代初始时间、迭代终止时间。

3.根据权利要求1所述的基于并行计算机的地球动力数模算法，其特征在于，所述步骤五中，依据速度场改变模型内部离散点的位置、应力参数，依据温度场、压力场改变模型内部离散点的温度、粘度参数。

4.根据权利要求1或3所述的基于并行计算机的地球动力数模算法，其特征在于，所述步骤五，还包括如下小步：

第一步：令初始时刻t＝0，迭代最大时间为tmax；

第十步：令t＝t+△t，更新当前迭代时间；

第十一步：将当前时刻的迭代数据存储到计算机的硬盘中。

5.根据权利要求4所述的基于并行计算机的地球动力数模算法，其特征在于，所述步骤五的判断方法如下：

6.根据权利要求1所述的基于并行计算机的地球动力数模算法，其特征在于，所述步骤七中的参数场图包括温度场图、速度场图、密度场图。

7.根据权利要求1或6所述的基于并行计算机的地球动力数模算法，其特征在于，所述步骤八中，依据密度计算声速与声阻抗界面、声速体位置，并绘制成图声速图、声阻抗界面图、声速体位置图；依据温度计算地表热流值，并绘制成地表热流值分布图；利用应力应变信息计算断层可能出现的位置并绘制成断层位置图；利用密度、空间位置数据可以计算重力异常图；利用其它信息绘制类似的地球物理、地球化学图件。

8.根据权利要求7所述的基于并行计算机的地球动力数模算法，其特征在于，所述步骤八中，将图件保存在计算机的硬盘中，令计数器N＝N+1；如果N大于Nnum则正常结束运行，否则继续读取下一个状态文件。