CN108460195B - 海啸数值计算模型基于gpu并行的快速执行方法 - Google Patents
海啸数值计算模型基于gpu并行的快速执行方法 Download PDFInfo
- Publication number
- CN108460195B CN108460195B CN201810130177.1A CN201810130177A CN108460195B CN 108460195 B CN108460195 B CN 108460195B CN 201810130177 A CN201810130177 A CN 201810130177A CN 108460195 B CN108460195 B CN 108460195B
- Authority
- CN
- China
- Prior art keywords
- tsunami
- gpu
- cpu
- calculation model
- mathematical calculation
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Expired - Fee Related
Links
- 238000000034 method Methods 0.000 title claims abstract description 54
- 238000004364 calculation method Methods 0.000 title claims abstract description 53
- 230000008569 process Effects 0.000 claims abstract description 22
- 238000012545 processing Methods 0.000 claims abstract description 14
- 238000012546 transfer Methods 0.000 claims abstract description 3
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 claims description 29
- 230000006870 function Effects 0.000 claims description 24
- 230000004907 flux Effects 0.000 claims description 10
- 238000004422 calculation algorithm Methods 0.000 claims description 7
- 230000001133 acceleration Effects 0.000 claims description 6
- 238000009826 distribution Methods 0.000 claims description 5
- 238000013316 zoning Methods 0.000 claims description 4
- UZVHFVZFNXBMQJ-UHFFFAOYSA-N butalbital Chemical compound CC(C)CC1(CC=C)C(=O)NC(=O)NC1=O UZVHFVZFNXBMQJ-UHFFFAOYSA-N 0.000 claims description 3
- 230000005484 gravity Effects 0.000 claims description 3
- 238000004321 preservation Methods 0.000 claims description 3
- 241000208340 Araliaceae Species 0.000 claims description 2
- 235000005035 Panax pseudoginseng ssp. pseudoginseng Nutrition 0.000 claims description 2
- 235000003140 Panax quinquefolius Nutrition 0.000 claims description 2
- 230000005540 biological transmission Effects 0.000 claims description 2
- 230000000694 effects Effects 0.000 claims description 2
- 235000008434 ginseng Nutrition 0.000 claims description 2
- 230000002123 temporal effect Effects 0.000 claims description 2
- 239000011800 void material Substances 0.000 claims description 2
- NAWXUBYGYWOOIX-SFHVURJKSA-N (2s)-2-[[4-[2-(2,4-diaminoquinazolin-6-yl)ethyl]benzoyl]amino]-4-methylidenepentanedioic acid Chemical compound C1=CC2=NC(N)=NC(N)=C2C=C1CCC1=CC=C(C(=O)N[C@@H](CC(=C)C(O)=O)C(O)=O)C=C1 NAWXUBYGYWOOIX-SFHVURJKSA-N 0.000 claims 1
- 101000612856 Dictyostelium discoideum Probable serine/threonine-protein kinase tsuA Proteins 0.000 description 93
- 238000005516 engineering process Methods 0.000 description 8
- 230000008901 benefit Effects 0.000 description 5
- 230000006854 communication Effects 0.000 description 5
- 238000003860 storage Methods 0.000 description 4
- 238000004891 communication Methods 0.000 description 3
- 238000010586 diagram Methods 0.000 description 2
- 238000012986 modification Methods 0.000 description 2
- 230000004048 modification Effects 0.000 description 2
- 238000004088 simulation Methods 0.000 description 2
- 238000003325 tomography Methods 0.000 description 2
- RAFZYSUICBQABU-HMMYKYKNSA-N Phytal Chemical compound CC(C)CCCC(C)CCCC(C)CCC\C(C)=C\C=O RAFZYSUICBQABU-HMMYKYKNSA-N 0.000 description 1
- RAFZYSUICBQABU-QYLFUYDXSA-N Phytal Natural products CC(C)CCC[C@@H](C)CCC[C@@H](C)CCC\C(C)=C/C=O RAFZYSUICBQABU-QYLFUYDXSA-N 0.000 description 1
- 238000010521 absorption reaction Methods 0.000 description 1
- 230000015572 biosynthetic process Effects 0.000 description 1
- 238000005520 cutting process Methods 0.000 description 1
- 238000007405 data analysis Methods 0.000 description 1
- 238000013135 deep learning Methods 0.000 description 1
- 230000001419 dependent effect Effects 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 238000006073 displacement reaction Methods 0.000 description 1
- 230000005489 elastic deformation Effects 0.000 description 1
- 230000002706 hydrostatic effect Effects 0.000 description 1
- 230000006872 improvement Effects 0.000 description 1
- 230000010365 information processing Effects 0.000 description 1
- 230000010354 integration Effects 0.000 description 1
- 238000010801 machine learning Methods 0.000 description 1
- 230000009022 nonlinear effect Effects 0.000 description 1
- 238000005457 optimization Methods 0.000 description 1
- RAFZYSUICBQABU-UHFFFAOYSA-N phytenal Natural products CC(C)CCCC(C)CCCC(C)CCCC(C)=CC=O RAFZYSUICBQABU-UHFFFAOYSA-N 0.000 description 1
- 238000006116 polymerization reaction Methods 0.000 description 1
- 238000009877 rendering Methods 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 239000013535 sea water Substances 0.000 description 1
- 239000004575 stone Substances 0.000 description 1
- 238000006467 substitution reaction Methods 0.000 description 1
- 230000001360 synchronised effect Effects 0.000 description 1
- 238000003786 synthesis reaction Methods 0.000 description 1
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F9/00—Arrangements for program control, e.g. control units
- G06F9/06—Arrangements for program control, e.g. control units using stored programs, i.e. using an internal store of processing equipment to receive or retain programs
- G06F9/30—Arrangements for executing machine instructions, e.g. instruction decode
- G06F9/30003—Arrangements for executing specific machine instructions
- G06F9/30007—Arrangements for executing specific machine instructions to perform operations on data operands
- G06F9/3001—Arithmetic instructions
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F9/00—Arrangements for program control, e.g. control units
- G06F9/06—Arrangements for program control, e.g. control units using stored programs, i.e. using an internal store of processing equipment to receive or retain programs
- G06F9/30—Arrangements for executing machine instructions, e.g. instruction decode
- G06F9/38—Concurrent instruction execution, e.g. pipeline or look ahead
- G06F9/3877—Concurrent instruction execution, e.g. pipeline or look ahead using a slave processor, e.g. coprocessor
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F9/00—Arrangements for program control, e.g. control units
- G06F9/06—Arrangements for program control, e.g. control units using stored programs, i.e. using an internal store of processing equipment to receive or retain programs
- G06F9/46—Multiprogramming arrangements
- G06F9/50—Allocation of resources, e.g. of the central processing unit [CPU]
- G06F9/5005—Allocation of resources, e.g. of the central processing unit [CPU] to service a request
- G06F9/5027—Allocation of resources, e.g. of the central processing unit [CPU] to service a request the resource being a machine, e.g. CPUs, Servers, Terminals
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
- Y02A90/00—Technologies having an indirect contribution to adaptation to climate change
- Y02A90/10—Information and communication technologies [ICT] supporting adaptation to climate change, e.g. for weather forecasting or climate simulation
Landscapes
- Engineering & Computer Science (AREA)
- Theoretical Computer Science (AREA)
- Software Systems (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- General Engineering & Computer Science (AREA)
- Evolutionary Computation (AREA)
- Geometry (AREA)
- Computer Hardware Design (AREA)
- Computational Mathematics (AREA)
- Mathematical Analysis (AREA)
- Mathematical Optimization (AREA)
- Pure & Applied Mathematics (AREA)
- Geophysics And Detection Of Objects (AREA)
Abstract
本发明公开了一种海啸数值计算模型基于GPU并行的快速执行方法，方法包括：CPU接收海啸数值计算模型的启动指令后，获取与所述海啸数值计算模型相关的参数信息，并获取初始海表面形变信息；CPU读取到数据传输指令时，CPU将参数信息和初始海表面形变信息传输至GPU；GPU通过内核函数执行所述海啸数值计算模型的海啸传播计算过程，获取最终的计算结果；GPU读取到数据回传指令时，GPU将海啸计算结果传输至CPU；CPU根据输出指令，对计算结果进行格式化处理并输出。上述方法将计算量小的过程交由CPU执行，计算量大的过程通过内核函数实现在GPU端并行处理，进而可显著提高海啸数值计算模型的执行效率，提高了海啸快速预警能力。
Description
技术领域
本发明属于海啸预警和灾害评估技术，尤其涉及一种海啸数值计算模型基于GPU并行的快速执行方法。
背景技术
海啸数值计算模型是实时海啸预警和海啸研究的有效手段之一，二维浅水动力学方程是模型通常的数学物理基础。
海啸数值计算模型的运行需要硬件支持，如工作站、小型工作站或者大型计算机，运行方式包括串行和并行两种。通常，并行方式的计算效率更高，基于CPU的并行技术包括OpenMP(Open Multi-Processing)和MPI(Message Passing Interface)。OpenMP利用单计算节点内含CPU的多核心和共享内存并行提速，可扩展性差，MPI采用多计算节点和分布式内存，可扩展性好。计算节点指的是大型计算机上由1-2颗CPU组成用于指定任务处理的硬件集成环境，一旦执行某项作业，其他作业便无法介入，作业执行阶段具有独占性。
利用OpenMP技术实现单计算节点多计算核心的并行加速的具体方案是在海啸数值计算模型循环代码中加入OpenMP引导语句，通过切割将计算范围分成若干子区域，每个子区域负责执行与其他子区域间无依赖关系的计算代码。每次循环完成之后，在共享内存内完成数据通讯、分配和聚合。计算核心指的是CPU内部的物理核心数量。
上述两种并行方式的缺点：基于MPI的并行技术需要多个计算节点和高速交换机，硬件成本高，代码学习难度大；基于OpenMP并行技术的代码学习成本较低，但是计算性能则受到单个计算节点内含CPU的计算核心数量限制。目前Intel公司最先进的Core i9处理器也只有18个物理核心。上述两种并行加速方案阻碍了越洋海啸数值计算在海啸预警系统中效能的进一步提升。
当前的海啸计算模型都是通过高性能计算机在其中央处理器CPU上执行。英伟达公司(NVIDIA)在1999年发布GeForce256显卡时，首先提出GPU(Graphic Processing Unit)的概念，最初仅将其用于数字图形处理和动画渲染。随着GPU计算能力的不断提高，GPGPU(General Purpose GPU)的概念被提了出来，其应用场景正在向通用计算领域不断渗透。与传统CPU相比，GPU拥有更强的计算性能和更高的访存带宽，是一个天然并行的、数据间无相互依赖关系的纯净计算环境。
2006年英伟达公司发布第一款基于CUDA(Compute Unified DeviceArchitecture)架构的GeForce 8800GTX显卡，配合2007年推出的CUDA_C语言和可编程性越来越强的SDK(Software Development Kit)，GPU逐渐成为当前高性能计算系统中最重要的加速部件，为开发人员有效利用GPU强大的计算性能提供了有利条件。越来越多的算法被成功移植到GPU芯片上执行。
目前，GPU加速技术已在世界各地的政府、实验室、大学、企业以及中小企业得到广泛的应用。随着信息化社会的飞速发展，云计算、大数据分析、深度学习等新技术的不断涌现，人们对计算机信息处理能力的要求越来越高。基于GPU的高性能计算不仅能够应用在石油勘探、航天国防、天气预报等传统领域，互联网、金融、大数据以及机器学习等新领域对高性能计算的需求也在飞速增长。
GPU在高性能计算领域表现出了巨大加速潜力，如何利用GPU对海啸数值计算模型的核心计算模块进行加速的成为当前需要解决的问题。
发明内容
针对现有技术中的问题，本发明提供一种海啸数值计算模型基于GPU并行的快速执行方法，能够有效提高海啸数值计算模型的整体执行效率。
第一方面，本发明提供一种海啸数值计算模型基于GPU并行的快速执行方法，包括：
CPU接收海啸数值计算模型的启动指令后，获取与所述海啸数值计算模型相关的参数信息，并获取初始海表面形变信息；
所述CPU读取到数据传输指令时，所述CPU将所述参数信息、初始海表面形变信息传输至GPU；
所述GPU接收所述参数信息、初始海表面形变信息后，执行所述海啸数值计算模型的海啸传播计算过程，获取海啸的计算结果；
所述GPU读取到数据回传指令时，所述GPU将所述海啸的计算结果传输至所述CPU；
所述CPU根据输出指令，将所述海啸的计算结果进行格式化处理并输出。
可选地，所述GPU接收所述参数信息、初始海表面形变信息后，执行所述海啸数值计算模型的海啸传播计算过程，获取海啸的计算结果的步骤，包括：
所述CPU在确定需要获取至少一个时间节点的海啸计算结果时，所述CPU向所述GPU发送至少一个时间节点的海啸中间计算结果获取指令；
所述GPU根据所述至少一个时间节点的海啸中间计算结果获取指令，将计算过程中符合所述时间节点的海啸中间计算结果传输至所述CPU。
可选地，还包括：
所述CPU将所有的海啸中间计算结果进行格式化处理并输出。
可选地，所述CPU将所述参数信息、初始海表面形变信息传输至GPU，包括：
所述CPU将所述参数信息、初始海表面形变信息通过PCI-E 3.0接口传至所述GPU的显存。
可选地，所述GPU将所述海啸的计算结果传输至所述CPU，包括：
所述GPU将所述海啸的计算结果通过所述PCI-E 3.0接口回传所述CPU的内存。
可选地，所述方法还包括：
所述CPU向所述GPU传输显存释放指令，所述GPU根据显存释放指令，将所述GPU使用的显存释放。
可选地，所述CPU和所述GPU位于同一设备中，所述GPU通过内核函数执行所述海啸数值计算模型的海啸传播计算过程；
所述CPU通过主机函数执行所述海啸数值计算模型中除海啸传播计算过程之外的处理过程。
可选地，所述主机函数为利用Fortran语言编写的CPU负责执行的代码，所述内核函数(也称设备执行函数)为利用CUDA_C语言编写的GPU负责执行的代码；且所述海啸数值计算模型的执行过程由所述CPU主导完成。
本发明具有的有益效果如下：
第一、本发明将具有超多计算核心的GPU(型号为Tesla K40)作为CPU(型号为XeonE5-2620v4)的协处理器，采用CUDA_C语言编写内核函数(Kernel Function)将海啸计算模型的海啸传播计算模块整体移植到GPU端执行，并且将计算需要的参数信息和初始海表面形变信息一次性传递至GPU显存，从而提高海啸数值计算模型的整体计算效率。
第二、相对于现有技术中OpenMP并行的处理方式，利用GPU线程级并行加速方案，将海啸计算模型的海啸传播计算模块通过CUDA_C编写的内核函数整体移植到GPU端执行，充分利用GPU内部数量远超CPU的计算核心对数值计算模块进行加速。(重点是利用GPU更多的计算核心数量)
第三、海啸数值计算模型需要大量浮点运算，内存访问速度(也称访存带宽)是制约海啸计算模型效率的另一个关键因素。CPU与内存通过主板连接，现有技术条件下的速度约为3～6GB/s，即便是加上集成在CPU上的分级缓存器，其综合访存带宽一般不超过70GB/s。本发明中将海啸数值计算模型中海啸传播计算模块的参数信息和初始海表面形变信息通过内核函数一次性传输至GPU，之后每次内核函数调用只传递指令和少量参数。充分利用GPU独特的三级显存结构-依次是网格(Gird)、块(Block)和线程(Thread)，它们各自具有独享的内存空间，每个线程有自己的寄存器和本地内存，每个块有自己的共享内存，每个网格还有全局内存、常量内存和纹理内存，综合访存带宽可达CPU的5～10倍。(重点利用GPU更快的访存带宽速度)
第四、CPU擅长按串行思路执行代码，GPU擅长大规模的并行计算。利用Fortran语言和CUDA_C语言混合编程，将海啸数值计算模型中流程控制和逻辑判断模块通过Fortran语言编程，代码在CPU端执行，充分利用CPU对复杂问题的处理能力；针对海啸传播计算模块，即求解浅水差分方程，包括质量连续方程和动量方程以及海啸波吸收边界算法，需要的大规模数组运算，则利用CUDA_C语言编程将计算需要的相关参数和变量一次性上传至GPU端，并且通过内核函数求解方程，充分利用GPU的并行计算天赋和访存带宽优势；在需要输出结果的时间节点或者整体计算过程结束时，将海啸计算结果从GPU端回传至CPU，完成格式化输出。本发明极大提高了海啸数值计算模型的整体执行效率。
即，本发明实施例中利用Fortran和CUDA_C语言混合编程，通过内核函数实现了将球坐标下海啸传播线性浅水方程移植到GPU端进行空间差分和时间积分计算；利用Fortran和CUDA_C语言编程，通过内核函数实现了将球坐标下海啸传播非线性浅水方程移植到GPU端进行空间差分和时间积分计算；利用Fortran和CUDA_C语言编程，通过内核函数实现海啸波吸收边界算法在GPU端执行；由此，完成对海啸数值计算模型整体执行效率的提升。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其它的附图。
图1为现有技术中的交错式蛙跃差分法T+1时刻求解的示意图；
图2为本发明实施例提供的方法流程示意图。
具体实施方式
为了更好的解释本发明，以便于理解，下面结合附图，通过具体实施方式，对本发明作详细描述。
在以下的描述中，将描述本发明的多个不同的方面，然而，对于本领域内的普通技术人员而言，可以仅仅利用本发明的一些或者全部结构或者流程来实施本发明。为了解释的明确性而言，阐述了特定的数目、配置和顺序，但是很明显，在没有这些特定细节的情况下也可以实施本发明。在其它情况下，为了不混淆本发明，对于一些众所周知的特征将不再进行详细阐述。
海啸数值计算模型一般采用Fortran语言进行编程，加上模型需要链接的数学函数库，编译之后生成可执行文件，然后提交给计算机执行。
海啸计算模型运行首先要声明变量和分配内存、读取模型配置参数和地震断层参数，然后，根据计算区域生成网格。本模型采用球坐标下曲面正交网格，将高精度的洋底地形和水深格点数据插值到网格上。
随后，根据断层参数计算海底地震引起的初始海平面形变(等价于海底形变)。通常采用基于弹性形变理论的Okada模型来估计海底形变量，也可以直接以文件的形式直接读取初始海表面形变信息
海啸波本质上是长周期小振幅重力波的传播，可以采用二维水动力模型对海啸波进行模拟刻画。在大洋深水区，海啸波长很长，传播过程的能量损失很小，一般采用线性浅水方程模拟海啸波的传播过程，球坐标系下的二维线性浅水方程可表示如下：
式(1)至(3)中，η为相对于平均海平面的自由表面位移，h为静水水深，P、Q分别为沿经度方向和纬度方向单位距离的通量；为纬度，ψ为经度，R为地球半径，f为科氏力系数，g为重力加速。
当海啸波传播至陆架浅水区时，非线性效应逐渐显现，因此，应采用非线性浅水方程对海啸传播过程进行刻画，球坐标下方程如下：
式(5)、(6)中，H为总水深，Fx、Fy分别为经度和纬度方向的底摩擦力，其它变量符号意义同上。
考虑到理论求解线性/非线性方程的难度较大，通常借助计算机利用数值求解的方式去获得一个逼近真值的计算解。本发明采用有限差分法来求解浅水长波动力方程，差分方法为交错式显性蛙跃法。有限差分法是一种求偏微分(或常微分)方程的数值解的方法，在空间区域的边界上需要满足的定解条件称为边值条件，本发明在空间边界采用吸收边界条件。在初始时刻所要满足的定解条件，称为初值条件，本发明初始条件为初始海表面形变信息。交错式显性蛙跃法以球坐标系下的线性差分方程(1-3离散化)为例进行说明：
每个变量含义同方程(1)～(3)，变量下标表示网格中的位置，变量上标表示时间点，如图1所示。η、P和Q三个变量在空间上交错分布，当前时刻的η值由上一时刻该变量和其他变量唯一确定，变量P和Q情况相似。
一个完整海啸数值计算模型现有加速方案是，在CPU端利用OpenMP技术将计算区域根据核心数量平均切割，每个核心在自己的范围内负责一个进程，顺序执行各自代码，在共享内存完成数据通讯，并且必须保持同步执行。
技术缺点主要是两个方面：1)单计算节点内部的CPU计算核心数量有限，海啸传播计算模块的代码并行率不高，降低了计算效率；2)对于需要大量浮点运算的海啸计算模型而言，大量数据的内存吞吐不可避免，内存访问速度是另一制约模型效率的核心指标。
本发明在现有技术的基础上，利用CUDA_C语言编程将最耗机时的海啸传播计算模块整体移植到GPU端执行，充分利用并行程度更高、计算核心更多的GPU芯片对模型进行加速；与此同时，利用CUDA_C编写内核函数，将核心计算模块需要的参数和初始海表面形变量通过PCI-E3.0接口一次性传递至GPU显存，充分利用GPU芯片内部更高的显存带宽，优化计算核心和存储器之间的通讯。在需要输出海啸计算结果的时间节点，通过指令将GPU端的数据通过PCI-E 3.0接口回传至CPU端进行输出。通过上述两方面的技术发明，海啸数值计算模型的整体执行效率有了显著的提升。
总结，本发明利用GPU对一个海啸数值计算模型的海啸传播计算模块进行加速，模型的输入、输出以及流程控制仍由CPU来负责执行。充分利用两种芯片的特长，完成对海啸计算模型整体执行效率的提升。
本发明公开了一个高性能地震海啸数值并行计算模型，该模型利用Fortran语言和CUDA_C语言混编实现了基于中央及处理器(CPU)和图形处理器(GPU)异构的高性能海啸传播并行计算。其中，Fortran编程的部分在CPU端(主机端)执行，包括主机端声明变量和分配内存、模型参数读取，生成计算网格和插值水深、计算初始海面形变量以及计算结果输出；CUDA_C编程部分在GPU端执行，包括声明设备端变量和分配内存、从主机向设备传递数组、循环求解离散化的二维浅水动力方程及从设备向主机回传计算结果。流程图详见图2，具体实现步骤如下：
1)搭建CPU+GPU通用高性能计算服务器和软件环境，软件环境包括C、C++、Fortran编译器和CUDA Toolkit等。方案：直接采购GPU图形服务器和与之匹配的高性能计算环境软件，软件包括icc、ifort、nvcc等，还应至少预装NetCDF库。
2)根据地震发生位置和海啸影响区域确定模型计算区域，设定海啸计算空间范围、空间分辨率Δx、时间分辨率Δt以及计算时长；根据空间范围和空间分辨率计算生成球坐标系下正交曲线网格以及相应的网格参数；读取原始地形水深文件，插值得到地形水深网格数据；计算Courant系数(见公式10)，据此判断模型稳定性。
其中，g表示重力加速度，hmax为计算范围内的最大水深。方案：利用Fortran语言编程，在CPU端声明所有变量并赋初值。编写子程序读入一个固定格式的模型参数文件，根据计算区域和空间分辨率计算正交曲面网格每一点的经纬度并插值水深；编写子程序对Courant系数是否满足稳定性条件进行判断，若不满足，对时间步长进行调整。
3)获取海底断层破裂引起的海底形变信息，即初始海表面形变量。在假设海水不可压缩的条件下，可以根据Okada断层公式，结合地震震级、断层的长度、宽度、走向角、倾角、滑动角以及滑移量计算得到初始海表面形变，也可以以文件的形式直接读取。方案：利用Fortran语言编程，根据Okada断层模型，结合读入的地震和断层参数计算初始海表面形变；或者直接读入文件格式的初始海表面形变量和动量通量信息。
4)将海啸传播计算模块需要的变量和参数通过PCI-E 3.0接口传递至GPU显存。方案：利用CUDA_C语言的cudaMalloc命令，在GPU端声明变量、分配显存，然后再用cudaMemcpy(*，cudaMemcpyHostToDevice)命令，实现从CPU向GPU传递所有参与海啸传播计算的变量和参数，包括海平面初始垂向位移、沿经度和纬度动量通量。
5)海啸波在大洋深水区的传播计算，即数值求解线性浅水方程，包括质量连续方程(详见方程1)、经向和纬向动量方程(详见方程2-3)和定解边界条件。方案：利用CUDA_C语言编写在设备端执行的内核函数__global__void**_kernel(…)，包括离散化的质量连续方程(详见方程7)、经向和纬向动量方程(详见方程8-9)和海啸波吸收边界算法。采用交错式显性蛙越方案求解方程(详见图1)，网格中心点是海啸波高和水深，上下及相邻四边为沿经向和纬向的动量通量。时间上，GPU执行逐时间步长迭代计算；空间上，GPU执行沿经度和纬度两个方向求解；当执行至网格边界时，采用吸收边界算法滤波。同一时间步长，GPU必须遍历所有网格点，然后才能进入下一个时间步长计算，循环往复直至设定的计算时长。方程的解包括海啸波高和两个方向的海啸动量通量。该方案最大优点在于，下一时刻整个模拟范围内所有计算格点值(如海啸波幅，海啸动量通量)取决于上一时刻该格点及其周边格点的数值。当前时刻每个格点的数值之间不存在依赖关系，这是利用GPU并行加速的技术基础。
6)海啸波在浅水区域的传播计算，即数值求解非线性浅水方程，包括质量连续方程(详见方程4)、经向和纬向动量方程(详见方程5-6)和吸收边界算发；动量方程中增加了底摩擦项和非线性对流项。方案：基本同5)，仅在内核函数中的非线性动量方程增加了底摩擦项和非线性对流项。
7)输出保存海啸计算结果，采用NetCDF格式对结果进行输出保存。方案：利用CUDA_C语言的cudaMemcpy(*，cudaMemcpyDeviceToHost)命令，实现当前时刻海啸计算结果从GPU向CPU的回传；利用Fortran语言编程，调用NetCDF库实现标准化格式数据的文件输出，调用命令包括定义并打开文件nf90_create、定义变量维度nf90_def_dim、定义变量nf90_def_var、定义变量属性nf90_put_att、结束文件定义nf90_enddef、输入变量nf90_put_var、关闭文件nf90_close。
8)待全部计算过程结束，释放内存。方案：利用CUDA_C语言的cudaFree命令，释放GPU内存。
其中步骤2-4和7-8不需要逐个时间步长迭代，计算量很小，在CPU端执行。步骤5-6需要对大规模数组进行时间迭代，因此，利用CUDA_C语言编程，通过内核函数实现其在GPU端的高性能并行计算，待计算结束，结果回传至CPU端，完成结果的格式化输出和保存。
所述模型首次将天生具备并行能力的GPU用于海啸数值计算，作为CPU的协处理器对海啸传播计算模块进行加速，并且充分利用GPU更大的访存带宽，极大的提升了模型计算效率。实际表现为，无论太平洋上任何地点发生引发海啸的海底强震，模型都可以在不到一分半钟的时间里计算得到海啸传播至中国沿岸的海啸波幅。本发明将显著提高我国的海啸预警能力。
需要明确的是，本发明并不局限于上文所描述并在图中示出的特定配置和处理。为了简明起见，这里省略了对已知方法的详细描述。在上述实施例中，描述和示出了若干具体的步骤作为示例。但是，本发明的方法过程并不限于所描述和示出的具体步骤，本领域的技术人员可以在领会本发明的精神后，做出各种改变、修改和添加，或者改变步骤之间的顺序。
还需要说明的是，本发明中提及的示例性实施例，基于一系列的步骤或者装置描述一些方法或系统。但是，本发明不局限于上述步骤的顺序，也就是说，可以按照实施例中提及的顺序执行步骤，也可以不同于实施例中的顺序，或者若干步骤同时执行。
最后应说明的是：以上所述的各实施例仅用于说明本发明的技术方案，而非对其限制；尽管参照前述实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述实施例所记载的技术方案进行修改，或者对其中部分或全部技术特征进行等同替换；而这些修改或替换，并不使相应技术方案的本质脱离本发明各实施例技术方案的范围。
Claims (6)
1.一种海啸数值计算模型基于GPU并行的快速执行方法，其特征在于，包括：
CPU接收海啸数值计算模型的启动指令后，获取与所述海啸数值计算模型相关的参数信息，并获取初始海表面形变信息；
所述CPU读取到数据传输指令时，所述CPU将所述参数信息、初始海表面形变信息通过PCI-E3.0接口一次性传输至GPU的显存；
所述GPU接收所述参数信息、初始海表面形变信息后，执行所述海啸数值计算模型的海啸从深海传播至近岸计算过程，其中，所述海啸数值计算模型采用有限差分方法求解球坐标系下二维线性方程和非线性方程，获取海啸的计算结果；所述非线性方程考虑底摩擦效应；
所述GPU读取到数据回传指令时，所述GPU将所述海啸的计算结果传输至所述CPU；
所述CPU根据输出指令，将所述海啸的计算结果进行格式化处理并输出；
其中，用于海啸数值计算的GPU作为CPU的协处理器对海啸传播计算模块进行加速；所述GPU接收所述参数信息、初始海表面形变信息后，执行所述海啸数值计算模型的海啸传播计算过程，获取海啸的计算结果的步骤，包括：
所述CPU在确定需要获取至少一个时间节点的海啸计算结果时，所述CPU向所述GPU发送至少一个时间节点的海啸中间计算结果获取指令；
所述GPU根据所述至少一个时间节点的海啸中间计算结果获取指令，将计算过程中符合所述时间节点的海啸中间计算结果传输至所述CPU；
具体地，第一、根据地震发生位置和潜在海啸影响区域确定海啸数值计算模型计算区域，设定海啸计算空间范围、空间分辨率Δx、时间分辨率Δt以及计算时长；根据空间范围和空间分辨率采用公式(10)计算生成球坐标系下正交曲线网格以及相应的网格参数，即读取原始地形水深文件，插值得到地形水深网格数据；计算Courant系数；
其中，g表示重力加速度，hmax为计算范围内的最大水深；
即利用Fortran语言编程，在CPU端声明所有变量并赋初值，编写子程序读入一个固定格式的模型参数文件，根据计算区域和空间分辨率计算正交曲面网格每一点的经纬度并插值水深；编写子程序对Courant系数是否满足稳定性条件进行判断，若不满足，对时间步长进行调整；
第二、获取海底断层破裂引起的海底形变信息，即初始海表面形变量；即利用Fortran语言编程，根据Okada断层模型，结合读入的地震和断层参数计算初始海表面形变；或者直接读入文件格式的初始海表面形变量和动量通量信息；
第三、将海啸传播计算模块需要的变量和参数通过PCI-3.0接口传递至GPU显存；即利用CUDA_C语言的cudaMalloc命令，在GPU端声明变量、分配显存，然后再用cudaMemcpy(*，cudaMemcpyHostToDevice)命令，实现从CPU向GPU传递所有参与海啸传播计算的变量和参数，包括海平面初始垂向位移、沿经度和纬度动量通量；
第四、海啸波在大洋深水区的传播计算，即数值求解线性浅水方程，包括质量连续方程、经向和纬向动量方程和定解边界条件；
即利用CUDA_C语言编写在设备端执行的内核函数__global__void**_kernel(…)，包括离散化的质量连续方程、经向和纬向动量方程和海啸波吸收边界算法；采用交错式显性蛙越方案求解方程，网格中心点是海啸波高和水深，上下及相邻四边为沿经向和纬向的动量通量；
时间上，GPU执行逐时间步长迭代计算；空间上，GPU执行沿经度和纬度两个方向求解；当执行至网格边界时，采用吸收边界算法滤波；同一时间步长，GPU必须遍历所有网格点，然后才能进入下一个时间步长计算，循环往复直至设定的计算时长；方程的解包括海啸波高和两个方向的海啸动量通量；
第五、海啸波在浅水区域的传播计算，即数值求解非线性浅水方程，包括质量连续方程、经向和纬向动量方程和吸收边界算发；动量方程中增加了底摩擦项和非线性对流项；
第六、输出保存海啸计算结果，采用NetCDF格式对结果进行输出保存。
2.根据权利要求1所述的方法，其特征在于，还包括：
所述CPU将所有的海啸中间计算结果进行格式化处理并输出。
3.根据权利要求1所述的方法，其特征在于，所述GPU将所述海啸的计算结果传输至所述CPU，包括：
所述GPU将所述海啸的计算结果通过所述PCI-3.0接口回传所述CPU的内存。
4.根据权利要求1所述的方法，其特征在于，所述方法还包括：
所述CPU向所述GPU传输显存释放指令，所述GPU根据显存释放指令，将所述GPU使用的显存释放。
5.根据权利要求1至4任一所述的方法，其特征在于，
所述CPU和所述GPU位于同一设备中，所述GPU通过内核函数执行所述海啸数值计算模型的海啸传播计算过程；
所述CPU通过主机函数执行所述海啸数值计算模型中除海啸传播计算过程之外的处理过程。
6.根据权利要求5所述的方法，其特征在于，所述主机函数为利用Fortran语言编写的CPU负责执行的代码，所述内核函数为利用CUDA_C语言编写的GPU负责执行的代码；且所述海啸数值计算模型的处理过程由所述CPU主导完成。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201810130177.1A CN108460195B (zh) | 2018-02-08 | 2018-02-08 | 海啸数值计算模型基于gpu并行的快速执行方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201810130177.1A CN108460195B (zh) | 2018-02-08 | 2018-02-08 | 海啸数值计算模型基于gpu并行的快速执行方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN108460195A CN108460195A (zh) | 2018-08-28 |
| CN108460195B true CN108460195B (zh) | 2019-06-14 |
Family
ID=63239638
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201810130177.1A Expired - Fee Related CN108460195B (zh) | 2018-02-08 | 2018-02-08 | 海啸数值计算模型基于gpu并行的快速执行方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN108460195B (zh) |
Families Citing this family (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN110727026B (zh) * | 2019-10-10 | 2021-05-28 | 太原理工大学 | 一种基于地震-海啸连续作用的非弹性位移比谱模型 |
| CN112906224A (zh) * | 2021-02-23 | 2021-06-04 | 中国地质大学（北京） | 一种基于cuda c的短波辐射模式加速方法 |
| CN114895358B (zh) * | 2022-07-13 | 2022-10-14 | 海南浙江大学研究院 | 地震海啸快速预警方法和系统 |
| CN114966854B (zh) * | 2022-07-28 | 2022-11-18 | 海南浙江大学研究院 | 一种单层非静压地震海啸预警系统和方法 |
| CN116187118A (zh) * | 2022-11-25 | 2023-05-30 | 国家海洋环境预报中心 | 基于多gpu并行的多目标区域精细化海啸传播计算方法 |
| CN117687779B (zh) * | 2023-11-30 | 2024-04-26 | 山东诚泉信息科技有限责任公司 | 基于异构多核计算平台的复杂电波传播预测快速计算方法 |
Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102663208A (zh) * | 2012-04-28 | 2012-09-12 | 浪潮电子信息产业股份有限公司 | 一种利用gpu进行大涡模拟实时渲染的方法 |
| CN103425833A (zh) * | 2013-08-07 | 2013-12-04 | 湖南大学 | 一种基于熵格子波尔兹曼模型的并行流体计算实现方法 |
| CN104951601A (zh) * | 2015-06-04 | 2015-09-30 | 大连理工大学 | 一种海冰-海洋结构相互作用的离散元高性能仿真系统 |
| CN106842320A (zh) * | 2017-01-19 | 2017-06-13 | 北京大学 | Gpu并行三维地震波场生成方法和系统 |
Family Cites Families (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN105974471B (zh) * | 2016-07-19 | 2018-11-23 | 中国地质大学（北京） | 一种基于异步流的地震数据多gpu快速正演计算方法 |
| CN106683184B (zh) * | 2017-01-04 | 2019-09-13 | 朱军 | 网络环境下泥石流灾害过程快速模拟与可视化分析方法 |
| CN107515987A (zh) * | 2017-08-25 | 2017-12-26 | 中国地质大学（北京） | 基于多松弛格子玻尔兹曼模型的地下水流动模拟加速方法 |
-
2018
- 2018-02-08 CN CN201810130177.1A patent/CN108460195B/zh not_active Expired - Fee Related
Patent Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102663208A (zh) * | 2012-04-28 | 2012-09-12 | 浪潮电子信息产业股份有限公司 | 一种利用gpu进行大涡模拟实时渲染的方法 |
| CN103425833A (zh) * | 2013-08-07 | 2013-12-04 | 湖南大学 | 一种基于熵格子波尔兹曼模型的并行流体计算实现方法 |
| CN104951601A (zh) * | 2015-06-04 | 2015-09-30 | 大连理工大学 | 一种海冰-海洋结构相互作用的离散元高性能仿真系统 |
| CN106842320A (zh) * | 2017-01-19 | 2017-06-13 | 北京大学 | Gpu并行三维地震波场生成方法和系统 |
Non-Patent Citations (3)
| Title |
|---|
| 基于GPU并行算法的水动力数学模型建立及其效率分析;赵旭东 等;《大连理工大学学报》;20140331;第54卷(第2期);正文第204-209页 |
| 基于GPU并行计算的浅水波运动数值模拟;许栋 等;《计算力学学报》;20160229;第33卷(第1期);正文第114-121页 |
| 海啸淹没范围的数值研究;侯京明;《中国优秀硕士学位论文全文数据库 基础科学辑》;20130315;正文第8-16页 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN108460195A (zh) | 2018-08-28 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN108460195B (zh) | 海啸数值计算模型基于gpu并行的快速执行方法 | |
| CN103135132B (zh) | Cpu/gpu协同并行计算的混合域全波形反演方法 | |
| Petaccia et al. | OpenMP and CUDA simulations of Sella Zerbino Dam break on unstructured grids | |
| Yuan et al. | FUNWAVE‐GPU: Multiple‐GPU acceleration of a Boussinesq‐type wave model | |
| Satria et al. | GPU acceleration of tsunami propagation model | |
| Xu et al. | Optimizing finite volume method solvers on Nvidia GPUs | |
| CN109490948A (zh) | 地震声学波动方程矢量并行计算方法 | |
| CN105974471B (zh) | 一种基于异步流的地震数据多gpu快速正演计算方法 | |
| Horne et al. | A hardware accelerated unstructured overset method to simulate turbulent fluid flow | |
| Abbasnia et al. | OpenMP parallelism in computations of three-dimensional potential numerical wave tank for fully nonlinear simulation of wave-body interaction using NURBS | |
| Topa et al. | GPGPU implementation of cellular automata model of water flow | |
| Kim et al. | GPU-accelerated boussinesq model using compute unified device architecture FORTRAN | |
| Deng et al. | CPU/GPU computing for an implicit multi-block compressible Navier-Stokes solver on heterogeneous platform | |
| CN106353801A (zh) | 三维Laplace域声波方程数值模拟方法及装置 | |
| Mintu et al. | Application of GPGPU to accelerate CFD simulation | |
| Fang et al. | A GPU accelerated Boussinesq-type model for coastal waves | |
| Amouzgar et al. | A GPU-accelerated shallow flow model for tsunami simulations | |
| Sataporn et al. | Parallel Algorithms of Well‐Balanced and Weighted Average Flux for Shallow Water Model Using CUDA | |
| Zhu et al. | Parallel optimization of underwater acoustic models: A survey | |
| Alawneh et al. | Ice simulation using GPGPU | |
| CN110162804A (zh) | 基于cpu加速的波场正演模拟优化方法 | |
| Degtyarev et al. | Virtual testbed: ship motion simulation for personal workstations | |
| Geist et al. | Lattice-boltzmann water waves | |
| Dawes et al. | Acceleration of an unstructured hybrid mesh RANS solver by porting to GPU architectures | |
| Brechtken | Lattice group models: GPU acceleration and numerics |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant | ||
| CF01 | Termination of patent right due to non-payment of annual fee | ||
| CF01 | Termination of patent right due to non-payment of annual fee |
Granted publication date: 20190614 |