并行化洪水演进模拟研究综述

李 健，张大伟，姜晓明，向立云

1.中国地质大学（武汉）海洋学院，武汉430074

2.中国水利水电科学研究院防洪减灾所，北京100038

3.水利部防洪抗旱减灾工程技术研究中心，北京100038

洪水风险图编制是防汛减灾非工程措施的重要手段，对降低洪水灾害至关重要。但在实际应用中，洪水、降雨和地形条件复杂多变，实际汛情与计算方案情景差别很大，制约了洪水风险图的实际应用效果[1]。实施可针对多种情景、任意边界条件、实时演示和人机交互式计算的洪水演进模拟，为防汛决策等提供动态、实时的洪水淹没信息，具有重要的现实意义。但是，基于水动力计算的洪水演进模拟计算量较大时，模拟耗时延长，将妨碍洪水淹没的实时演示和防洪决策，并行化计算可解决此类问题[1]。

洪水演进模拟是计算水力学领域的一个重要研究方向，涉及繁杂的数值算法、研究手段和编程技术等，目前已几乎不可能对其进行全面详尽的综述。然而，实际应用中多采用基于Godunov 型有限体积算法和网格计算的洪水演进模型，涉及到复杂地形处理、干湿前锋追踪、通量格式和谐性、模型鲁棒性等问题[2]，对此国内外学者从不同角度撰写了大量综述性论文，包括：洪水模拟的不同物理机制[3]、洪水演进模型的复杂度分析[4]、洪水淹没的不确定性分析[5]、城市洪水淹没模拟[6]。近年来计算机硬件和软件的发展，尤其异构并行计算技术的发展，极大地推动了洪水演进模型的研发。但是，目前仅有覃金帛等[7]介绍GPU并行优化在水文计算和水库调度研究中的应用，未见有系统性介绍并行化洪水演进模拟的综述性论文。因此，有必要及时总结洪水演进的并行化研究成果，促进今后的洪水演进过程的并行化模拟研究。

考虑到洪水演进模拟研究的主流性、前沿性和实际应用价值，本文将介绍基于网求解浅水方程的并行化洪水演进模型，尤其是非结构网格模式。由于粒子法应用于实际河道或流域的洪水模拟时，存在边界条件难以定义和计算量大等问题，以及计算机图形学中的流体模拟，侧重计算机可视化和高性能制图算法，而流体模拟多采用无条件稳定的快速傅里叶变换（FFT）等算法，并不是求解流体方程，因此，本文综述不包括上述的两种数学模型，粒子法或计算机图形学领域的并行化洪水演进模拟研究综述可参考文献[8]。本文将结合近年洪水演进模拟的最新数值算法和并行化技术，进行文献综述、讨论和展望，期望为今后快速洪水预报研究起到抛砖引玉的作用，为防洪减灾的实时会商和洪水风险图的编制提供技术支撑。

1 洪水演进模拟的并行化需求

近年来洪水演进模型被广泛应用于评估洪水造成的经济损失，洪水演进模型包括1D、2D 和3D 模型，在复杂度和预报能力等方面存在较大差异。1D模型可提供的信息有限，一般仅用于模拟复杂河网或与2D 模型耦合，实现区域洪水模拟，但存在人为处理边界条件的问题[9]。2D 模型采用结构网格或非结构网格求解浅水方程，精度高但计算量增大[10-12]，可采用自适应网格加密（AMR）技术[13]和3D 模型[14]的方法来克服，但实际应用中的计算量任然较大。因此，基于网格计算的并行化洪水演进模型的研发和应用势在必行。

2 洪水演进模型的并行化模式

常见的并行模式见表1，区别包括：并行层次（数据集、线程和进程）、并行粒度（细粒度和粗粒度）、硬件支持、软件实现和编写并行化代码的编程语言。本文涉及的硬件、软件和并行计算的计算机术语，仅列出对应的中文全称和英文缩略词，具体含义可参考文献[15-16]，后文中介绍的洪水演进模型使用的并行模式与表1 所列相对应。如果计算完全在中央处理单元（CPU）中完成，称为同构并行计算，如果计算部分执行于CPU，部分执行于其他硬件设备，如图形处理单元（GPU），称为异构并行计算。从并行层级的角度，表1列出了从最低层级的数据并行，到线程级和进程级并行。

表1 常见的并行技术列表

3 并行化洪水演进模拟研究现状

3.1 同构并行化洪水演进模型

洪水演进模型根据空间网格模式，可分为结构网格和非结构网格；根据时间离散格式，可分为显格式和隐格式；根据对控制方程的不同简化，可分为求解储水单元、扩散波、浅水方程和静水压力假设的Navier-Stokes方程等类型的模型；根据数值方法，可分为有限差分法（FDM）、有限体积法（FVM）和有限单元法（FEM）等模型。针对不同类型的洪水演进模型，需要采用不同的并行模式。本文将根据以上分类，从并行计算的角度，选择有代表性的模型列于表2，其中不包括未得到广泛应用或测试或未并行化的洪水演进模型[11-12，17-29]。另外，表2所列模型还考虑了代码的可获取性，方便读者选择使用，检验洪水淹没模拟结果的可重复性及模型对比分析。

表2 一些有名的同构并行化洪水演进模型列表

较早的洪水模型多使用结构网格和FDM，如LISFLOOD-FP和FloodMAP模型，结构网格模型数据结构简单，易于并行化编程，其中LISFLOOD-FP模型采用了OpenMP、MPI 和Clearspeed 三种并行机制，而Flood-MAP 采用多线程技术，LISFLOOD-FP 和FloodMAP 两种模型采用显格式近似求解黎曼问题的数值通量的方法，允许使用的计算时间步长较小，必须并行化才能应用于精细模拟区域洪水淹没。需要指出的是，不同的并行化洪水演进模型，求解的控制方程，使用的数值方法和并行化机制均不同，笼统地比较不同模型的并行计算效率并无意义，但并行化结构网格模型研究为之后的并行化模型研发奠定了基础。

近年来，洪水演进模拟已由结构网格计算转向非结构网格计算的模式开发，而FVM在单元界面上近似求解黎曼问题，具有较强的洪水激波捕捉能力，如MIKE21[17]、CalTWiMS[18]、ANUGA[19]、TUFLOW[20]和ParBreZo[21]等模型，但缺点是：

（1）对网格质量要求高，网格的正交性严重影响数值格式的计算精度，进而影响计算稳定性。

（2）采用显格式，当网格尺寸很小时，受CFL条件限制，允许使用的计算时间步长很小，严重影响计算效率。

（3）数值通量格式涉及大量数组循环计算。

（4）采用静态网格，可能在干湿地形剧烈变化的区域网格密度不足，影响计算精度。

采用FEM离散方法可克服第（1）点问题，并行化计算可克服第（2）、（3）点问题，其中基于OpenMP 的细粒度线程并行可加速计算数组循环计算，结合区域分解后的MPI 进程并行，可取得更好的加速计算效果。基于AMR动态网格技术可克服第（4）点问题，实现自动加密局部网格密度，跟踪干湿前锋位置。但改进的方法也存在以下局限性：

（1）AMR网格可能会产生正交性较差的网格，一般结合FEM模型使用[22]，已广泛应用于洪水或海啸传播过程的模拟，但增大了并行化的难度，特别是采用MPI 并行时，静态非结构网格可使用METIS[23]或SCOTCH[24]等库实现荷载均衡的网格区域分解，并可降低分区间的数据通信量，提高并行计算效率，而AMR动态网格密度调整会引起荷载不均衡，需要进行分区网格的动态调整，可使用Zoltan库[25]实现，但编程的难度会增大。

（2）AMR非结构网格还会破坏数据结构的局部性，这对GPU异构并行计算效率的影响尤其显著。

（3）局部高密度网格会引起计算不稳定问题。

（4）在非结构网格上难以构建高阶的空间离散格式。

（5）非结构网格隐格式模型，最终形成的具有稀疏系数矩阵的线性方程组求解量巨大，成为计算瓶颈，且系数矩阵的数据局部性对求解效率影响显著。

近年来，基于FEM 的洪水演进模型研究得到广泛开展（见表2），尽管FEM较FVM在编程上难度更大，但FEM 适应低质量网格的能力更高，因此计算稳定性较FVM 更好。FEM 可分为连续迦辽金（CG）法和间断迦辽金（DG）法。CG法结合MPI并行化，一般采用隐格式离散，可高效地模拟洪水演进，但CG法求解对流占优的不连续问题时，存在质量不守恒和格式的耗散性过强的问题，一般采用流线迎风Petrov-Galerkin（SUPG）法解决，如RMA2[26]和Telemac2D[27]模型，也有基于求解通用波动连续性方程（GWCE）的CG法来提高离散格式的质量守恒性，如ADCIRC模型[28]。DG法不要求单元界面处满足连续性条件，即具有FEM的灵活性又具有FVM的激波捕捉能力，同时可方便构建高阶格式，如Fluidity[22]和DGSWE[29]模型。值得注意的是，FEM 模型的数学原理抽象，编程难度较大，模型计算量大，需要并行化才能高精度地模拟洪水演进过程。FEM模型结合基于网格质量指标判断的AMR 动态网格技术，实现洪水演进过程中的干湿前锋的准确捕捉，可明显提高城市洪水和局部冲刷[30]等问题的模拟精度。求解大型线性方程组时，可有效利用多种并行模式的线性方程组求解器来提高计算效率，如PETSc 库[31]，同时可使用SCOTCH 库中的Cuthill-McKee 算法或Gibbs-Poole-Stockmeyer 算法，实施对网格单元和节点编号重排序，进一步改善非结构网格模型的求解效率。需要指出的是，高阶DG法中的高斯求积运算较大，需要合理设计网格分辨率，使计算结果在精度和计算效率上达到平衡。

3.2 异构并行化洪水演进模型

GPU 异构并行计算是利用显卡中的众多核心实现大量浮点运算的线程并行，非常适合用于加速FVM 洪水演进模型中数值通量涉及的数组循环计算，GPU异构并行的研究文献在近年来快速增长。下面选择有代表性的GPU异构并行洪水演进模型列于表3[17，20，32-41]。

表3 一些有名的异构并行化洪水演进模型列表

早期的GPU异构并行洪水演进模型，均采用CUDA C语言编程，且GPU的硬件架构可有效加速结构网格的FDM 模型的计算效率，如JFlow[32]和Flood2D[33]模型。由于FVM 模型在求解不连续问题方面的优势，异构并行的结构网格FVM 洪水演进模型也大量开展，如TUFLOW[20]和SWsolver[34]模型。非结构网格模型的数据结构较为复杂，异构并行编程的难度更大，且数据结构也会影响异构并行计算效率。因此，异构并行非结构洪水演进模型较异构并行结构网格模型的研究要晚，但已成为前沿研究课题，此类模型包括商业软件MIKE21、免费软件但不公开代码的BASEMENT[35]以及开源代码的VOLNA[36]模型。

异构并行非结构网格洪水演进模型存在的局限性有：（1）CUDA 编程复杂，会破坏串行程序的代码结构，代码维护成本高；（2）CUDA SDK版本进化和GPU硬件设备更新速度快，代码更新难度大；（3）CUDA程序只能在Nvidia 公司的显卡上执行，兼容性差；（4）非结构网格模式下，网格单元、节点和边的连接关系必须显式给出，这将导致大量的间接寻址访问，对数据传输带宽要求高，这在同构并行计算中问题不大，因为CPU内存的带宽为Gb/s 量级，而GPU 显存的带宽一般为Mb/s 量级，因此必须使用有利于提高间接寻址访问效率的数组形式。

针对以上问题，研究者们给出了解决方案，例如Zhang等[37]使用OpenACC指令集的异构并行方法，该方法类似于OpenMP指令的同构线程并行方法，对代码结构几乎没有破坏，但是只能用PGI编译器编译。还可采用OpenCL 语言编程，提高模型对硬件的兼容性，如HiPIMS 模型[38]。还有研究者使用特定域语言（DSL）的程序库，如OP2库[39]，可将串行版本的程序，快速转换为不同的同构和异构并行化的程序，如BASEMENT 和VOLNA 模型，DSL 转换可让科学家和工程人员专注本专业领域的研究，而不必关注并行化编程本身，同时可应对软件和硬件架构快速更替的问题。数组的排列形式分为Structure of Array（SoA）排列和Array of Structure（AoS）排列两种，SoA排列和AoS排列对同构并行计算效率的影响不大，但在异构并行计算时，SoA 排列具有显著优势，可提高间接寻址访问效率，同时实施非结构网格编号的重排序。

隐格式的洪水演进模型求解的瓶颈就是求解线性方程组，可调用异构并行的求解器，这样可利用具有异构并行计算资源的集群系统，如Fluidity[40]和CCHE2D[41]模型。值得关注的是，采用异构并行化的系数矩阵预处理的Krylov 空间求解器（如PETSc 库）和多重网格求解器（如AMGX 库[42]），可进一步提高求解效率，已应用于一些计算流体力学商业软件，但尚未见到用于隐格式洪水演进模拟的研究。

3.3 基于其他并行机制的洪水演进模型

除了基于MPI、OpenMP的同构并行和基于CUDA的异构并行洪水演进模型外，还有基于其他加速硬件设备和并行模式的洪水演进模型，此类并行化模式并未成为主流，但涉及最底层的数据并行或是对主流硬件加速研究的补充，也列于表4[11，34，43-44]。基于特殊处理器（如CELL和Clearspeed处理器）的并行化洪水模型，如SWSolver[34]和LISMIN-CS[11]模型。还有采用单指令多数据（SIMD）指令集（如SSE 和AVX）的数据并行，如CITEEC[43]和NUFSAW2D[44]模型。值得关注的是，NUFSAW2D 模型采用AVX 数据矢量化、MPI 和OpenMP 的混合并行，达到近1 000 倍的加速比，混合并行充分利用了并行计算硬件资源，在今后的洪水演进模型研发中应重视混合并行技术。由表4还可以看出：数据并行的洪水演进模型都采用结构网格，在非结构网格模型的数据并行和混合并行研发方面的加速计算效率还有待进一步研究。

表4 基于其他并行机制的洪水演进模型列表

4 讨论

首先，不同并行模式的洪水演进模型的研发技术难度、开发周期和维护成本存在较大差异，见表5，同构线程并行（OpenMP）和异构线程并行（OpenACC）仅需在串行代码中添加一定的指令即可，编程难度最低，代码的修改量最小，且不需要第三方程序库，采用具备一定编译功能的编译器即可编译和运行代码，其开发周期较短。集群并行模式需要使用第三方的MPI 库实施计算节点间数据通信，并实施网格分区间的数据收发、进程管理等，非结构网格区域分解还需要实施荷载均衡的第三方程序库，编程难度中等，且需要在集群系统上运行。GPU异构并行和数据并行，需要特定的SDK，对串行代码的修改量大，甚至要重新采用CUDA语言编程移植，编程难度最大，代码维护成本最高。

表5 不同并行化洪水演进模型的总结

其次，采用MPI+OpenMP 和MPI+CUDA 的混合并行技术，充分挖掘洪水演进模型的细粒度和粗粒度并行潜力[45]，可更进一步提高并行计算效率，见表5。MPI+CUDA 的异构并行称为GPU 集群并行，涉及多个加速器之间的数据交换，如何实现荷载均衡条件下多个设备间的数据交换速度，是近年来的研究热点，实现方法有两种：（1）通过CPU 主机实现多GPU 设备间的数据通信，称为三路数据通信；（2）多GPU 设备间直接数据通信，称为双路数据通信。后一种方式的数据交换速度最快，可使用Nvidia公司开发的GPUDirect技术来实现。

最后，使用类似FORTRAN 和C/C++这样的静态编程语言开发洪水演进模型，尽管具有很高的计算效率，但需要大量的代码编写，且容易出错，开发效率低；尽管FORTRAN或C++语言的面向对象编程可大大降低编程量[46]，但仍无法实现像FEM这样复杂模型的快速建模。使用动态编程语言（如Python），可快速实现FDM、FVM和FEM模型的建立[47-48]，结合Python语言的MPI、CUDA和OpenCL 等库的接口程序[49-50]，实现并行化计算。另外，由于FEM 算法易于形成抽象化范式，例如FEniCS项目建立了使用FEM，求解弱形式偏微分方程的统一形式语言（UFL）[51]，见表5，应用UFL建立了Python语言编程的Firedrake模型[52]，Firedrake模型实现了具有Python 语言接口的MPI、CUDA 和OP2 程序库[53]、具有SSE和AVX指令集优化的COFFEE库[54]、有限元自动组建的FIAT 库[55]、用于线性和非线性问题求解的PETSc库和非结构网格数据管理和编号的DMPlex库[56]之间的无缝耦合，使用FEM 建立了模拟洪水演进的Fluids 模型[57]和模拟海啸波传播的Thetis模型[58]，包括前处理、计算主程序和后处理，仅需要300 行Python 代码。使用DSL 和UFL 可快速实现对串行代码的并行化以及快速实施新的算法，上述不同并行洪水演进模型的优缺点及适用场景的具体总结见表5。

5 结论及展望

本文对具有代表性和实用性的并行化洪水演进模型进行了综述，获得以下几点结论：

（1）基于线程并行（OpenMP）、基于进程并行（MPI）和基于GPU 异构并行（CUDA）的洪水演进模型具有很好的实际应用价值，不同的混合并行模式（如MPI+OpenMP以及MPI+CUDA）将更能充分利用不同层级的计算资源，达到更高的并行加速计算效率。

（2）并行化非结构网格模式的洪水演进模型开发是发展趋势，但存在数据结构复杂、荷载均衡的区域分解以及数据结构对异构并行计算效率的影响较大等问题，可采用相关算法针对不同问题给出解决方案。

（3）使用特定域语言（DSL）和动态编程语言（如Python）开发并行化的洪水演进模型，可快速实施不同并行模式的洪水模拟加速效率的科学比较和新的洪水演进模拟算法，应该是今后的发展方向。

基于对并行化洪水演进模型的综述，对今后并行化洪水演进模型研发提出几点展望：

（1）非结构网格模式下的并行化洪水演进模型研发应成为着力点，应重点研究混合并行机制（如MPI+OpenMP+CUDA）、多GPU 间的数据直接交换技术（如GPUDirect技术）、解决并行化中自适应非结构网格加密的荷载不均衡问题（如使用Zoltan库）。

（2）针对GPU架构特点，研发适合于GPU异构并行的网格形式，如Block Uniform Quadtree（BUG）网格[59]，BUG网格的数据结构具有连续存储、结构紧致的特点，特别是结合BUG 网格和变时间步长技术[60]，可高速求解具有内边界[61]和复杂河网的洪水演进过程。探讨适于GPU加速的网格形式和数据结构仍有很大空间。

（3）基于计算机图形库OpenGL的实时洪水演进可视化和交互式计算是未来的发展趋势[62]，但如何实现并行计算（包括CPU集群并行和GPU集群并行）环境下的实时可视化和交互式计算，仍是有待解决的问题。