周生昌刘卫国宋振亚杨晓丹

(1.山东大学 软件学院,山东 济南250101;2.自然资源部 第一海洋研究所,山东 青岛266061;3.青岛海洋科学与技术试点国家实验室 区域海洋动力学与数值模拟功能实验室,山东 青岛266237;4.海洋环境科学和数值模拟自然资源部重点实验室,山东 青岛266061)

研究海洋的手段主要有理论研究、观测试验和数值模拟,其中观测试验是数值模拟和理论研究的基础,理论研究为数值模拟和观测试验设计提供指导,而数值模式则集成了人们从观测数据和理论研究中所获得的知识和理解,不仅可以用来验证理论假说和弥补观测资料的不足,而且是能够最终用于预测未来变化的唯一工具[1]。随着海洋和气候变化研究的不断深入,高性能计算机的发展,海洋数值模式正逐步朝着更高分辨率、更多物理过程和更快计算速度的方向发展[2]。一般来说,水平分辨率每提高1倍,计算量会增大到原来的8～10倍[3]。同时,随着分辨率的提高,会有越来越多的原次网格不能描述的物理过程被包含进来,而新的次网格过程的参数化也需要通过将多组数值试验结果与实测数据进行比较来进行评估和调整,这显然也提高了对计算量和计算能力的要求[4]。因此,充分利用和发挥高性能计算机的性能来提高海洋数值模式计算速度已成为模式发展的一个必要条件。

目前高性能海洋数值模拟的研究主要采用以下几种方案提高模式计算能力:1)基于分布式内存并行化程序。主要采用MPI(Message Passing Interface)消息传递接口使计算任务合理分布在多个独立进程的计算单元内,并通过高速网络进行数据通信。该方法实现复杂,需考虑边界交换、负载均衡等问题,虽然可扩展性好,但节点之间交换会带来通讯开销。2)基于共享式内存并行化程序。主要采用Open MP等接口将计算任务分配给多个线程,并通过共享内存地址的方式实现数据交换。该方法实现相对容易,无通讯开销,但一般来说无法跨节点,这限制了其并行规模。3)基于进程/线程并行化程序。该方法是以上2种并行方式的结合,实现较为复杂,在异构计算机(如CPU+GPU、国产申威处理器等)上效果特别显著。

上述3种方案主要从并行扩展方面提高模式计算性能,这需要大量的计算资源,且程序在设计和实现中并未考虑代码现代化。特别是随着高性能计算机的发展,这些方案已经无法充分利用和发挥现代计算机的优势。杨晓丹等[3]在海浪模式MASNUM上开展了代码现代化优化方面的尝试,通过检查代码的规范性、合理性,并针对语言特性对代码进行改进,使之更好地适应计算机基础架构,充分发挥了现代计算机的计算特点。但是在此海浪模式中,计算量最大的源函数与相邻点无关,无分区间数据交换通讯,并且由于分支判断少,cache命中率高[5],这与海洋环流模式NEMO(Nucleus for European Modeling of the Ocean)[6]计算特点有较大的差异。NEMO不但存在大量的数据交换,而且分支判断较多。因此如何使用代码现代化方法优化和提高海洋环流模式的计算性能,非常具有代表性和必要性。

为了充分利用和发挥现代高性能计算机的性能,进一步促进海洋环流数值模式的应用和发展,本文基于代码现代化的概念,提出了海洋环流模式代码现代化优化的方案,并在此基础上对海洋环流模式NEMO进行了优化,分析和探讨了代码现代化优化方案在海洋环流模式应用中的前景和局限性,为今后海洋环流数值模拟的研究和发展提供参考和建议。

1 优化方案设计

1.1 优化策略

现代高性能计算机是由多核和众核处理器、高速缓存和内存、高带宽处理器以及它们之间的通信结构和高速I/O组合而成的[7]。为了充分适应平台的现代化结构,充分发挥平台的性能,需要发展高性能和现代化的软件,包括对原始代码进行改进,以及为现代计算机重新设计应用程序以获得最大性能,但这些均需要充分利用高性能计算机硬件资源,这是代码现代化的基础。

本文基于代码现代化优化策略,并结合NEMO海洋环流模式的特点以及现有平台架构,制定了6个阶段的优化策略,如图1所示。为了保证结果的准确性,排除平台因素引起的误差,每次均取2次试验结果的平均值作为最终结果。

1.2 试验算例

试验算例为全球海洋环流真实算例ORCA1_LIM,该试验算例只包含海洋动力学、热力学(OPA模块)以及海冰动力学(LIM模块)。模式共积分120个时间步长,水平网格数目为362×292,垂向网格数目为75。为了减小I/O对扩展性的影响,本研究中的测试是针对模式计算性能的无I/O测试。

图1 代码现代化优化策略Fig.1 Code modernization optimization strategy

1.3 测试环境

试验测试环境如表1所示,单节点内共有两颗Intel Xeon Gold 6252处理器,每颗处理器有24核心,因此单节点共有48物理核心,96逻辑核心,内存为384 GB DDR4。

表1 测试环境Table 1 The testing environment

2 NEMO优化

2.1 定位热点函数

为了更准确地获取影响NEMO模式计算效率的瓶颈,从而更有针对性地对NEMO模式进行优化,首先使用Intel Vtune Amplifier工具定位出模式的热点函数,即耗时较多的函数,如图2所示。结果显示NEMO的热点函数极为分散,即热点函数较多,耗时分布较为均衡,本研究选取前五个热点函数(tra_ldf_iso,lim_rhg,tra_adv_tvd,ldf_slp、nonosc)进行重点优化。

图2 串行模式各函数计算时间Fig.2 Serial mode function calculates time

2.2 编译选项优化

在性能分析的基础上,本研究首先使用编译选项对模式进行优化。使用x HOST选项使编译器生成处理器支持的最高指令集,该选项在本文测试环境下等效于x MIC-AVX512;使用ipo启用文件之间的过程间优化,使编译器对单独文件中定义的函数执行内联函数扩展;使用no-prec-div选项对模式的浮点数除法进行优化,将除法改为乘倒数的形式进行计算,从而提高运算效率;此外,还使用了fp-model fast=2,qopt-dynamic-align,qopt-prefetch等优化选项,具体如表2所示。在经过编译选项优化后,NEMO的整体加速比可以达到1.21倍。通过对输出结果进行分析,以海表面温度为例,得出O1和O3两种编译选项所引起的误差绝对值不超过9.3×10-6℃,占仅占平均海表面温度的0.00005%,且误差绝对值在1.0×10-8℃以下的网格点占比超过99.8%,满足精度要求。同样,对其它变量的检查也满足精度要求。因此编译选项优化是有一种有效的优化手段。

表2 编译选项优化列表Table 2 The list of compilation option optimization

2.3 标量串行优化

标量串行优化的目的是通过对热点函数代码进行调整修改,如去除重复计算、减少条件分支、降低循环嵌套等,确保代码使用最少的计算量和合适的精度获得正确的结果。下面以traldf_iso热点函数为例介绍此优化过程。

在如下例子中(图3a),zdk1t和zdkt计算公式是相同的,即zdk1t(ji,jj,jk-1)=zdkt(ji,jj,jk),因此通过将zdk1t从循环中移除,在计算完zdkt后再将结果赋值给zdk1t,可以减少重复性计算,如图3b所示:

图3 重复计算优化前和优化后Fig.3 Repeat calculation before optimization after optimization

当从内存中访问数组元素时,两次访问的寻址间隔会影响Cache命中率,降低计算效率。由于Fortran中的数组是按照列优先的方式存储,因此调整数组元素访问顺序,使得相邻两次访问的寻址间隔变小,可以有效提高Cache命中率。调整数组元素访问顺序前后分别如图4a和4b所示。

另外,减少循环嵌套也能有效提高模式计算效率。在图5a的循环中,z2d是临时数组,目的是存储zftu在三维空间沿K轴方向的积分。由于三重循环计算效率较低,因此利用Fortran提供的forall命令将该循环减少为一重循环,如图5b所示。经过上述串行与标量优化后,NEMO模式整体计算效率提升1.22倍。

图4 调整数组元素访问顺序优化前和优化后Fig.4 Adjust access order of array element before optimization after optimization

图5 循环嵌套优化前和优化后Fig.5 Loop nesting before optimization after optimization

2.4 SIMD优化

SIMD(Single Instruction Multiple Data)即单指令多数据流,以同步方式,在同一时间内对多个数据执行同一条指令。矢量化是现代计算机的一个标志,利用自动矢量化,可以显著提高软件的计算性能。在本测试环境下,利用x HOST选项编译时,已包含了自动矢量化功能。然而自动矢量化需要遵循矢量化规则,如循环中不能存在数据前后依赖以及跳转语句等。以下是自动矢量化失败的例子及采用的优化策略。

在图6a所示的代码中,由于zdit和zdjt都是由指针方式实现的动态数组,Fortran指针可以作为别名指向一块内存地址,因此编译器不能判断zdit和zdjt指针之间是否存在数据依赖,从而无法实现自动矢量化。通过将假定依赖的语句拆分成2个循环,使每个循环满足自动矢量化的条件,可以实现自动矢量化(图6b)。

图6 SIMD数据依赖优化前和优化后Fig.6 SIMD data dependency before optimization after optimization

另外,由于代码中存在大量的指针作为临时数组,因此在不影响物理过程的前提下,通过将指针修改为动态数组,可以使代码安全地进行自动矢量化,如图7所示。经过SIMD优化后,NEMO模式整体计算效率提升至1.23倍。

图7 SIMD指针优化前优化后Fig.7 SIMD pointer before optimization after optimization

2.5 内存带宽优化

当CPU需要读取内存中的数据时,CPU将首先发出一个由内存控制器执行的请求,然后再将数据由内存返回中央处理器,此过程的时间称为延时周期。为了缩短延时周期,中央处理器和内存之间设置了L1和L2缓存。由于L1和L2均采用静态随机访问的原理,可以将其理解为内存带宽[8]。对内存带宽进行优化,提高缓存命中率,也是提升软件计算效率的有效方式。

由于海洋环流模式复杂的物理过程,模式中使用大量的临时数组保存计算过程的中间结果,这就造成了运算时参与计算的数组过多的问题。当CPU计算时,如果一条语句需要同时读取大量数组,由于内存带宽限制,则会减少缓存命中率,降低计算效率。因此减少或避免使用临时数组,是缓解内存带宽占用过高的一种优化手段。

如图8a所示,数组A,B为临时数组,运算时首先通过计算为A,B数组赋值,并在之后的循环中访问该数组。优化策略为将临时数组A、B去掉,而把计算赋值的功能重构为子程序cal_a,cal_b,并在循环中需要此中间结果时调用。

图8 内存带宽优化优化前和优化后Fig.8 Memory bandwidth before optimization after optimization

对热点函数进行分析后,发现traldf_iso、nonosc和lim_rhg中均存在大量临时数组,可以使用上述策略进行优化。图9显示了优化前后各个函数的计算时间以及其加速比,可以看出优化后单个热点函数的计算效率有效提升10%～40%,效果显著。另外,由于NEMO的热点函数呈现分散的特点,此优化使得模式整体计算效率加速至1.31倍。

图9 内存带宽优化前后函数时间对比Fig.9 Comparison of function time before and after memory bandwidth optimization

2.6 MPI多进程扩展

NEMO模式总体是一个循环过程,其中时间是循环的主序列。循环从数据读取开始,以结果输出为结束。循环部分包含模式中的所有核心计算,总计算量占总循环过程的95%及以上[9]。本文使用的NEMO 3.6版本基于MPI消息传递接口编写,使用MPI分布式内存并行计算可以有效的进行任务划分,从而提高计算效率。在经过以上串行优化之后,接下来利用Intel Trace Analyzer Collector工具对NEMO进行多进程扩展测试,以观察模式的负载均衡性。

图10显示了模式在48进程时的负载均衡性,从结果可以看出,计算(蓝色)和通信(红色)时间占比比较合理,且分配给每个进程的计算任务基本一致,其中进程间的计算时间占比最大相差不超过8%,有40个进程相差在1%以内,基本实现了负载均衡,不需要进行进一步优化。

图10 48进程下负载均衡分析Fig.10 Load balancing analysis of 48 processes

为了更全面地对模式的并行能力进行检验,我们进一步测试了模式的扩展性。图11为海洋环流模式在不同进程数下的加速比。在单节点48进程内,随着进程数的增加,加速效率逐渐降低,这是由于增加进程引起的代价不可消除,如进程间通信量增加,带宽限制等。另外,程序中必定有不能并行的串行部分,即使进程无限增多,通过并行所产生的加速比都是受限的,即阿姆达尔定律[10]。因此模式无法随着进程的增加获得线性加速比。此外,当进程数超过48时,由于测试环境没有开启超线程,进程可能会被挂起等待,同时进程间通信增加,加速效率会继续降低。

在经过编译选项优化、标量串行优化、SIMD优化、内存带宽优化几个步骤后,模式整体性能得到了有效提升,如图12所示。选取合适的编译选项可以充分发挥计算机的性能,使得计算效率提升了21%;通过规范代码书写,去除重复计算,减少循环嵌套等方法,实现串行代码优化,以及通过修改代码,使之符合编译器对矢量化的要求,实现编译器自动矢量化,这两步优化共实现了2%的性能提升;另外,针对模式占用内存带宽过高的问题,提出了去掉临时数组的解决方案,使得优化后的单个热点函数加速比最高可以达到1.4倍,模式整体加速了8%。经过上述优化,最终模式的整体计算效率共提升了31%,节省了接近1/3的计算资源。同时,48进程下的MPI并行模式计算效率相较于串行模式提升了15.8倍,虽然单节点内扩展性较弱,但负载均衡性较好。

图11 不同进程下的加速比Fig.11 Acceleration ratio of different processes

图12 各优化步骤加速比Fig.12 Acceleration ratio of each optimization step

3 试验结果分析

3.1 大量指针被假定数据依赖

矢量化有许多限制,例如存在数据依赖、函数调用(数学库调用除外)、在同一循环中混合可矢量化类型以及存在与数据相关的循环退出条件等,则不能进行矢量化。其中数据依赖分为写后读依赖和读后写依赖。写后读依赖即计算一个数组元素时需要依赖该数组已经计算的元素,如果强制将其矢量化,可能后面需要读取的依赖元素并没有被计算完成,因此造成读取数据错误。除了存在不可避免的计算中的前后依赖之外,NEMO的矢量化报告显示模式源代码中存在大量的指针,而由于Fortran中的指针可以作为别名指向内存地址,因此编译器通常无法判断包含指针的代码是否存在数据依赖,从而无法自动矢量化。本文采取了两种方法来避免因指针造成的无法矢量化问题,包括拆分循环和将指针修改为动态数组。但需指出的是,这只是针对热点函数的局部优化,因此效果有限。若要全局优化此问题,则需要将代码重构。

3.2 循环嵌套过多

NEMO的核心计算都包含在循环嵌套中,源代码中的函数多是四重循环,包括一维时间循环以及三维空间循环。当算例规模达到一定程度时,会不可避免地出现效率低下。另外,由于编译器仅会对最内层循环进行自动矢量化,因此循环嵌套过多也会造成SIMD优化效率低下。针对此问题,一种优化思路是遵循内大外小的原则,即将迭代次数多的循环放在最内层,提高SIMD优化效率。但更有效的解决方案则是减少循环嵌套层数,但是需要从根本上改变算法问题,这将是非常复杂的工程。

3.3 内存带宽过高

本研究中的串行与标量优化以及SIMD优化效果有限,而在经过内存带宽分析及优化后,热点函数的性能提升十分明显,这说明NEMO海洋环流模式内存带宽占用过高是个比较显著的问题。这是由于NEMO计算时涉及大量数组运算,因此会出现大量缓存未命中的情况,需要频繁从内存中读取数据。而当内存带宽占用过高时,处理器从内存获取数据的速度会出现瓶颈。去除临时数组是一种解决方案,但却会增加大量的重复计算,增大计算压力,这与串行代码优化中的去除重复计算优化策略恰恰相反。但针对具体的NEMO内存带宽瓶颈问题,去除临时数组后减少的时间远远大于增加重复计算带来的计算时间,因此去除临时数组是非常有效的。

4 结 语

在本文中,首先介绍了NEMO模式特征,并根据代码现代化指导步骤制定了NEMO的优化策略。随后给出了试验算例以及测试环境,并详细介绍了NEMO代码现代化的优化过程。结果表明,通过对选取的热点函数进行编译选项优化、标量串行优化、SIMD优化以及内存带宽优化,串行的海洋环流模式的整体加速比可以达到1.31倍,其中热点函数加速效果非常明显,例如traldf_iso相较于原始性能提升了2.3倍。但因模式热点函数分散,选取的热点函数并不是模式的全部计算瓶颈,因此整体加速比有限。最后对NEMO在MPI多进程下进行了扩展性测试,测试结果表明在48进程(测试环境最大核心数)的加速比为15.8倍,计算和通信时间占比合理,负载均衡。但可能是由于内存带宽占用过高,导致单节点内扩展性较差,这需要进行进一步分析验证。另外,我们还对优化后的模式结果进行了验证,以确保优化后结果的准确性。需要注意的是,虽然对于海洋模式NEMO来说,我们的优化对模拟结果影响不大,在误差允许范围内,但是在海气耦合模式中,由于存在强的非线性海气相互作用,误差可能会被放大,因此在优化时需要注意并检验。总的来说,代码现代化在不增加任何硬件成本的前提下是一种行之有效的优化方案。

此外,在对NEMO进行代码现代化优化的过程中,分析出其主要存在以下几点影响计算效率的问题,如模式大量使用指针、循环嵌套过多、内存带宽占用过高等。针对这些问题,本文给出了相应的解决方案,如将指针修改为动态数组,减少循环嵌套次数以及去除临时数组等,可以为今后海洋模式的设计和改进提供参考。

另外,MPI并行化仅在二维方向上进行分解(沿着jpi和jpj)[11]。当MPI进程数增多时,每个进程间的通信量会成倍增加,通信开销会越来越大,导致模式性能下降。为了减少通信时间,通过引入MPI/Open MP混合并行将会是一个切实可行的方法。Open MP可以实现共享内存层面的任务并行,它是一种简单有效的并行方式[11-12]。这种混合并行方法通过减少MPI进程数,可以减小进程间通信开销,同时利用Open MP打开超线程又能充分利用计算资源,使得并行结构更加合理[13]。由于部分NEMO算例的特征是在三维数组上沿着jpi、jpj和jpk执行的操作,因此可以将Open MP共享式内存并行化应用在jpk(此试验中为75)方向上,提高进程内的计算能力。在热点函数中引入MPI/Open MP并行化后,由于进程间通信数量减少,将有利于提高模式的扩展性。这是接下来研究的方向之一。