水力瞬变数值模拟方法

摘要

本发明提供一种水力瞬变数值模拟方法，能够有效地提高计算效率，其特征在于：主机端CPU负责逻辑运算；设备端GPU负责计算各断面的水力参数计算和更新；内核函数都以线程栅格的形式组织，组织形式为一维线性结构，每个线程对应一个计算断面，各断面的索引tx由线程所在块和块内位置唯一确定；边界函数由GPU或CPU计算；将内核函数运算结果先存储在寄存器中，循环结束后再将结果写入全局存储器中，最后复制回主机端；采用统一寻址方式，对压力和流量数组进行内存分配；当内核函数调用时，通过线程索引在数组中定位断面号，以获取相应断面压力、流量值进行计算，更新的水力参数值分别存储在已分配储存空间的另外两个数组中。

权利要求

1.一种水力瞬变数值模拟方法，其特征在于，包括：

(一)并行实现策略

主机端CPU负责逻辑运算，准备初始压力和流量数据，进行边界处理，将初始压力和流量数据复制到显存中、调用内核函数；

设备端GPU负责计算各断面的水力参数计算和更新，内核函数调用完成后将处理好的流量fq[i]和水头fh[i]数据传回到主机端内存中，作为下一时间步的输入数据；

(二)GPU计算内核函数

内核函数都以线程栅格的形式组织，内核函数线程结构的组织形式为一维线性结构，与管道的线形结构相匹配，计算时将管道划分为若干个计算断面，每个线程对应一个计算断面，每个线程的索引号与一个计算断面的编号相一一对应，各断面的索引tx由线程所在块和块内位置唯一确定，内核函数按照如下公式1至5编写：Cp＝HA+BQA BP＝B+R|QA|，    (公式3)CM＝HB-BQB BM＝B+R|QB|，     (公式4)式中，下标P代表需要求解的计算断面的值，下标A代表断面P左边相邻计算断面的上一时刻的值，下标B代表断面P右边相邻计算断面的上一时刻的值，H、Q分别为测压管压力和流量，B为管道特性阻抗，R为管道阻抗系数，a为水锤波速，f为摩阻系数，Δx为空间分段长度，D为管道直径，A为管道面积，g为重力加速度，系数CP、BP、CM、BM为中间变量；断面A和B的压力流量为已知，断面P的压力流量值通过公式1至5求解；

(三)边界函数

边界函数由GPU或CPU计算：

当边界由GPU计算时，将系统所有管道当作一条长管道，统一分配连续的线程索引，所有节点的计算均在GPU上进行；

当边界由CPU计算时，由边界将管道系统划分成多根管段，每个管段单独设置线程索引，节点均在CPU上计算，然后再传递到GPU上；

(四)程序优化策略

将内核函数运算结果先存储在寄存器中，循环结束后再将结果写入全局存储器中，最后复制回主机端；

(五)编程步骤

采用统一寻址方式，对压力和流量数组进行内存分配；当内核函数调用时，通过线程索引在数组中定位断面号，以获取相应断面压力、流量值进行计算，更新的水力参数值分别存储在已分配储存空间的另外两个数组中，计算完成后，将计算结果复制回主机内存中，作为下一时步的输入数据。

2.根据权利要求1所述的水力瞬变数值模拟方法，其特征在于：其中，在(三)边界函数中，设系统有i管段，分段数分别为N1、N2、N3……Ni，则它们的索引范围分别是0到N1-1、0到N2-1、0到N3-1、……、0到Ni-1；管段上第0、N1-1、N2-1、N3-1、Ni-

1节点均在CPU上计算，然后再传递到GPU上。

3.根据权利要求1所述的水力瞬变数值模拟方法，其特征在于：其中，在(五)编程步骤中，是采用Synchronize()函数实现计算的同步。

说明书

水力瞬变数值模拟方法

技术领域

[0001] 本发明属于计算流体动力学领域，具体涉及一种水力瞬变数值模拟方法。技术背景

[0002] 瞬变过程是管道输送系统中经常发生的压力、流量波动现象，相关计算和分析对系统的设计、运行、管理有重要指导意义。目前常用的瞬变过程数值求解方法中，一维特征线法和显式有限差分法都受库朗条件限制，当时间步长和空间段满足一定关系时才稳定，时间步长很小时才能保证较高精度；而隐式有限差分法 (如Pressimann法)虽不受库朗条件限制而无条件稳定，但其迭代求解非线性方程组的过程很复杂，且有较大的数值耗散。而系统瞬变过程计算往往需要较高的精度和计算效率。首先，长距离输水(油、气)系统管线一般长达数百公里，城市供配水管网通常包含非常多且复杂的回路和支路，对它们进行瞬变模拟须划分很细空间段和时间步长，计算十分耗时。另外，在复杂管道系统设计优化中，需要进行大量方案与工况的对比计算，要求很高的计算效率。此外，一些管路系统要求实时监测并进行瞬变过程仿真模拟，以对系统操作提供支持，这种情况对计算的实时性有很高要求。因此，非常有必要针对效率问题开发一种新的水力瞬变模拟方法，能够在不以降低计算精度为代价的前提下，缩短计算时间，为实际工程提供可靠的指导建议。

[0003] 目前，对提高瞬变过程计算效率的研究有不少，主要集中于两种途径，第一种是算法上采用更加高效的方法，第二种是在CPU上进行并行计算。Wood[1-3]、Jung[4]、Boulos[5]等人通过将波特性法(WCM)应用于大规模配水管网瞬变流模拟，由于WCM可取比MOC更大的时间步长，且无需每步都计算管道中间网格上的值，在一定程度上提高了计算效率。Izquierdo[6,7]编制通用水力元素边界处理子程序，将之在主程序中直接调用,避免过细的分段，减少计算资源占用，提高了数值模拟效率。Li[8]使用并行遗传算法(PGA)进行伴随空化的水力过渡过程模拟，相比于串行遗传算法(GA)，计算效率能提升大约10倍。Fan[9]建立了水力机械装置系统瞬变流并行计算的模型，实现了两台计算机对某大型水利水电工程全系统瞬变流的并行模拟，将计算效率提高了1.442倍。上述方法对计算效率的提高十分有限，远难满足现实需求。

[0004] 图形处理器(Graphic Processing Unit,GPU)在硬件设计上采用了超长流水线和大规模线程并行，具有超强的浮点运算能力，最初用于处理大量图像数据。本世纪初科研工作者开始尝试用它来处理大规模科学计算问题，与同等价位CPU 相比，加速比达到数十倍以上。但初始图形处理器的使用非常困难，因为比如像 OpenGL和Direct3D技术需要直接针对芯片编程，程序员必须通过应用程序编程接口来访问其核心，这些应用程序编程接口限制了能够为芯片编写的应用程序的类型。所以只有小部分程序员掌握了利用图形处理器运行特定程序的必要技术，进而限制了该项技术的广泛传播。一切随着NVIDIA公司在2007年推出CUDA 技术而改变。因为规避了传统的应用程序编程接口，通过CUDA，现在程序员们能够用他们所熟悉的类C/C++语言来直接针对芯片编程，从而使得并行编程技术即便对于初学者来说也非常容易上手，大大地推进了这项技术在科学领域的应用。

[0005] 一维特征线法(Method of Characteristics，MOC)由于计算效率高，精度高，编程简单等特点，广泛应用于工程水力瞬变模拟中。该方法在模拟包含复杂边界条件的系统上优势明显，但缺点是受稳定性条件限定而计算时间步长较小，计算时间跨度大的瞬变过程时较费时。

发明内容

[0006] 本发明是为了解决上述问题而进行的，目的在于提供一种基于一维特征线法的水力瞬变数值模拟方法，能够有效地提高计算效率。

[0007] 本发明为了实现上述目的，采用了以下方案：

[0008] 本发明提供一种水力瞬变数值模拟方法，其特征在于，包括：

[0009] (一)并行实现策略

[0010] 主机端CPU负责逻辑运算，包括变量的申明和赋值等，设备端GPU负责计算，以充分发挥CPU的超强逻辑能力和GPU的超高运算速度的优势。

[0011] 如图1所示，在CUDA编程架构中，一个程序分为两个部分：主机端(Host) 和设备端(Device)。主机端是指CPU上执行的部分，包括准备初始数据，进行边界处理，将数据复制到显存中、调用内核函数等；设备端则是在GPU上执行的部分，又称内核函数(Kernel)，它是整个程序的核心部分，负责各断面的水力参数计算和更新，内核函数调用完成后将处理好的流量fq[i]和水头fh[i]数据传回到主机端内存中，作为下一时间步的输入数据。

[0012] (二)GPU计算内核函数

[0013] 内核函数都以线程栅格(Grid)的形式组织，内核函数线程结构的组织形式为一维线性结构，如图1所示，与管道的线形结构相匹配。将管道划分为若干个计算断面，每个线程对应一个计算断面，由于采用线性排列，线程索引号与计算断面的编号是一致的，给定任何一个线程的索引号就能定位到相应的计算断面。给定了执行配置参数，各断面的索引tx由线程所在块和块内位置唯一确定，内核函数按照如下公式1至5编写代码：

[0014]

[0015]

[0016] Cp＝HA+BQA BP＝B+R|QA|，   (公式3)

[0017] CM＝HB-BQB BM＝B+R|QB|，   (公式4)

[0018]

[0019] 式中，下标P代表需要求解的计算断面的值，下标A代表断面P左边相邻计算断面的上一时刻的值，下标B代表断面P右边相邻计算断面的上一时刻的值，H、 Q分别为测压管压力(m)和流量(m3/s)，B为管道特性阻抗，R为管道阻抗系数， a为水锤波速(m/s)，f为摩阻系数，Δx为空间分段长度(m)，D为管道直径(m)， A为管道面积(m2)，g为重力加速度(m/s2)，系数CP、BP、CM、BM为中间变量；断面A和B的压力流量为已知，断面P的压力流量值通过公式1至5求解；t+Δt 时刻的其他断面亦是如此求解。以上求解的格式体现了MOC的显式性和局域性。显式性是指各离散点的值只与上一时刻的值有关，而与同时刻其他点的值无关，即、图1中P点的值与P1、P2、P3、P5…PN无关；局域性是指各离散点的值只与其相邻点的值有关，与相邻点以外的点无关，即、图1中P点的值只与其前后两点A、B(上一时刻)的值有关。

[0020] (三)边界函数

[0021] 边界函数可由GPU计算，也可由CPU计算。

[0022] 当边界由GPU计算时，将系统所有管道当作一条长管道，统一分配连续的线程索引，如系统总共分为N个计算断面，则索引的范围从0到N-1，所有节点的计算均在GPU上进行。

[0023] 当边界由CPU计算时，由边界将管道系统划分成多根管段，每个管段单独设置线程索引，例如系统有三根管段，分段数分别为N1、N2、N3，则它们的索引范围分别是0到N1-1、0到N2-1、0到N3-1。每根管段的计算过程如图1所示，第0、N1-1、N2-1、N3-1节点均在CPU上计算，然后再传递到GPU上。

[0024] (四)程序优化策略

[0025] MOC算法中内核函数代码执行15次浮点运算就需要访问全局存储器9次，浮点运算与全局存储器访问操作之间的比值CGMA(Compute to global memory access)为15/9，明显小于一般设备的理论值，无法充分发挥设备的计算能力。通过将内核函数运算结果先存储在寄存器中，循环结束后再将结果写入全局存储器中，最后复制回主机端，减少了对全局存储器的3次访问，CGMA值变成15/6，提高了访问的速率。

[0026] (五)编程步骤

[0027] 采用统一寻址(Unified Memory)方式，对压力和流量数组进行内存分配；当内核函数调用时，通过线程索引在数组中定位断面号，以获取相应断面压力、流量值进行计算，更新的水力参数值分别存储在已分配储存空间的另外两个数组中，计算完成后，将计算结果复制回主机内存中，作为下一时步的输入数据。

[0028] 本发明提供的水力瞬变数值模拟方法，还可以具有以下特征：在(三)边界函数中，设系统有i管段，分段数分别为N1、N2、N3……Ni，则它们的索引范围分别是0到N1-1、0到N2-1、0到N3-1、……、0到Ni-1；管段上第0、N1-1、 N2-1、N3-1、Ni-1节点均在CPU上计算，然后再传递到GPU上。

[0029] 本发明提供的水力瞬变数值模拟方法，还可以具有以下特征：在(五)编程步骤中，由于输入数据和输出数据是交替存放在两套数组中的，因此需要以奇偶次数交替调用内核函数，最后是采用Synchronize()函数实现计算的同步。

[0030] 基于以上方案，本发明提供的水力瞬变数值模拟方法，具体操作步骤为：

[0031] 1、初始化数据；

[0032] 2、将初始压力和流量数据从主机端内存拷贝到设备显存；

[0033] 3、调用内核函数进行并行计算，更新各断面的水力参数；

[0034] 4、同步所有线程；

[0035] 5、将各断面更新后的压力流量数据从设备拷贝回主机端，作为下一时步的输入数据；

[0036] 6、重复步骤2-5，直至达到计算终点。

[0037] 发明的作用与效果

[0038] 与现有技术相比，本发明具有以下优点和有益效果：

[0039] 1、本方法具有很高的并行性，能够获得数百倍于传统一维特征线法水力过渡过程模拟的计算速度，极大的提高了计算效率；

[0040] 2、本方法能够直接在并行机上实现，多处理器之间的通信极易优化，受计算机容量和速度的限制较小，能模拟计算密度大的的复杂输水工程的水力过渡过程；

[0041] 3、在本方法中，边界计算既可在GPU上也可在CPU上，灵活的边界处理能够带来更为有效的并行加速效果；

[0042] 4、在本方法算法稳定性较好，对大规模输水工程进行过渡过程模拟可取很小的时间步长，易满足稳定条件，适合于模拟实际问题；

[0043] 5、本方法算法简单，编程容易，便于推广。

附图说明

[0044] 图1为一维特征线法在CUDA架构上实现并行计算的示意图；

[0045] 图2为示例的水库-管道-阀门系统的结构示意图；

[0046] 图3为CPU-MOC和GPU-MOC方法模拟的各种水击现象阀端压力变化过程曲线，图3(a)是直接水击，图3(b)是一相水击，图3(c)是极限水击，图3(d)是负水击；

[0047] 图4为Karney文献中案例简单水力管道系统示意图；

[0048] 图5(a)至(f)分别为节点3、6、7处GPU-MOC程序、CPU-MOC程序和文献中的压力和流量过程曲线。

具体实施方式

[0049] 以下结合附图对本发明涉及的水力瞬变数值模拟方法的具体实施方案进行详细地说明。

[0050] <实施例一>

[0051] (一)模拟问题描述

[0052] 以一个简单的水库-管道-阀门系统为对象(图2)，其中水库水位恒定，瞬变过程由阀门开启或关闭引起。对直接水击、一相水击、极限水击和负水击四种典型工况进行计算，各工况计算条件见下表1(各工况均未考虑管道的摩阻损失)。由于系统只有一根管道，只含两个边界条件，将边界放在GPU上计算，采用单精度计算，结果见图3，由图3可知两种方法模拟计算的结果基本完全一致，且符合各种水击现象的压力变化规律。

[0053] 表1四种典型水击计算参数设置

[0054]

[0055] (二)并行计算加速效果

[0056] 1、GPU与CPU性能对比

[0057] 假定进行简单直接水击模拟的管道系统足够长，能够分成足够多的计算断面来体现GPU-MOC模型的并行计算能力。管道划分的计算断面范围为214～220，计算步数设为1000，分别测试GPU程序和CPU程序计算直接水锤的平均耗时，其中GPU程序在每个范围下又分别测试了线程块中分别含有32、64、128、256 个线程的平均耗时。

[0058] 表2显示了不同断面密度情况下单块GPU与单块CPU进行单精度运算的消耗时间和相应断面密度下GPU-MOC相对于CPU-MOC的加速比(Speedup)。表中第一列表示GPU-MOC模型划分的断面数，GPU-MOC列中行表头表示每个线程块中给定的线程数量。单元格中括号前的数据为计算1000步所需要的时间，单位为毫秒，括号中数据为加速比，即程序串行执行时间与并行执行时间比值。

[0059] 表2不同计算密度下GPU和CPU运行时间及加速比 单位：毫秒

[0060]

[0061] 2、GPU并行计算的设备利用率分析

[0062] 计算设置的断面密度范围是214～224，计算时间步数设为1000，每个线程块中设置256个线程。分别测试初始方法case 1和提高CGMA后的改进方法case 2 的运行时间，对比两种方法在不同格点密度下单块GPU每秒可更新百万个格点数(Million MOC reach Updates Per Second，MMUPS)和相应格点密度下带宽利用率(实际计算能力与理论计算能力的比值)。从下表3中可知，两种方法中GPU 设备的运算能力和带宽利用率均随着断面密度的增大而增加，且case 2带宽利用率一直高于case 1。

[0063] 表3两种方法内核函数的利用率

[0064]

[0065] <实施例二>

[0066] (一)问题描述

[0067] 模拟计算了一个Karney et al[10]文章中出现的Case 1(图4)，工程布置和参数设置均与文献中一致。管路总长大约9公里，管线中设有水库、调压室、解压阀、控制阀，瞬变过程由管路末端控制阀引起，其开度在10秒内由0.6线性关到0.2，解压阀在压力超过210m时会在30秒内线性打开，全开后60秒内线性关闭，其他参数可参考文献[11]。计算结果见图5，由图5可知CPU-MOC程序和 CPU-MOC程序的压力流量值完全一致，曲线完全重合，且它们与文献中的结果差别非常小。

[0068] (二)并行计算加速效果分析GPU和CPU计算时间对比

[0069] GPU-MOC程序计算时将管道空间分段长度分别设置成0.1m，0.05m，0.01m， 0.005m，0.001m，相应断面数为1821，90076，180152，900766，1801526，9007600，计算时间步数设为6500，线程块中线程数目设为256。分别测试了GPU程序和 CPU程序的平均耗时，其中GPU程序分别测试了两种边界处理方法，其中 Strategy 1中边界函数在CPU上计算，Strategy 2的边界函数在GPU上计算。

[0070] 下表4显示了不同断面密度下单块GPU与单块CPU进行单精度运算时间和相应断面密度下GPU的两种边界处理方法相对于CPU的加速比。从表中可以看出，在当前断面密度范围Strategy 1可实现63.9～126.4倍于CPU的加速， Strategy 2可实现12.7～316.7倍于CPU的加速，加速比随着线程总数的增加而增加。需要特别说明的是表格中Strategy 2的运行时间包含内核函数运行时间和 CPU计算边界的时间。很明显当断面密度小于900766时，Strategy 1的计算耗时更小，加速效果较好，断面密度大于900766后，Strategy 2的加速效果较Strategy 1有很大的提升，最大达到了Strategy 1的2.5倍。

[0071] 表4 GPU和CPU程序在不同网格密度下的运行时间和加速比 单位：秒

[0072]

[0073] 注：括号前是运行时间，单位秒，括号中是加速比。

[0074] (三)GPU并行计算的设备性能利用率分析

[0075] 下表5比较了Strategy 1和Strategy 2在不同断面密度下单块GPU每秒可更新百万个格点数和带宽利用率。其中Strategy 2中的时间是内核函数运行时间，不包含边界处理时间，计算时间由NVIDIA visual profiler分析软件获得。可以看出，Strategy 2的带宽利用率高于Strategy 1。

[0076] 表5两种策略的每秒网格更新数和设备利用率

[0077]

[0078] 以上实施例仅仅是对本发明技术方案所做的举例说明。本发明所涉及的水力瞬变数值模拟方法并不仅仅限定于在以上实施例中所描述的内容，而是以权利要求所限定的范围为准。本发明所属领域技术人员在该实施例的基础上所做的任何修改或补充或等效替换，都在本发明的权利要求所要求保护的范围内。

[0079] 本申请引证文件：

[0080] [1]Wood D J,Boulos P F,Lingireddy S.Pressure wave analysis of transient flow in pipe distribution systems[M].MWH SOFT,2005.

[0081] [2]Wood D J,Lingireddy S,Boulos P F,et al.Numerical methods for modeling transient flow in distribution systems[J].Journal(American Water Works Association), 2005,97(7):104-115.

[0082] [3]Wood D J,Dorsch R G,Lightner C.Wave-plan analysis of unsteady flow in closed conduits[J].Journal of the Hydraulics Division,1966,92(2):83-110.[0083] [4]Jung  B S,Boulos  P  F,Wood D  J,et  al.A  Lagrangian  wave characteristic method for simulating transient water column separation[J].American Water Works Association.Journal,2009,101(6):64.

[0084] [5]Boulos P F,Karney B W,Wood D J,et al.Hydraulic Transient Guidelines for Protecting Water Distribution Systems[J].Journal,2005,97(5):111-124.

[0085] [6]Izquierdo J,Iglesias P L.Mathematical modelling of hydraulic transients in simple systems[J].Mathematical and Computer Modelling,2002,35(7-8):801-812.

[0086] [7]Izquierdo J,Iglesias P L.Mathematical modelling of hydraulic transients in complex systems[J].Mathematical and computer modelling,2004,39(4-5):529-540.

[0087] [8]Li S,Yang C,Jiang D.Modeling of hydraulic pipeline transients accompanied with cavitation and gas bubbles using parallel genetic algorithms[J].Journal of Applied Mechanics,2008,75(4):041012.

[0088] [9]樊红刚,陈乃祥,杨琳.大型水利水电工程全系统瞬变流并行计算[J]. 清华大学学报:自然科学版,2006,46(5):696-699.

[0089] [10]Karney B W,McInnis D.Efficient calculation of transient flow in simple pipe networks[J].Journal of Hydraulic Engineering,1992,118(7):1014-1030。