一种核反应堆热工水力模拟软件CFD并行处理方法与流程

本发明涉及高性能计算和核科学技术与工程领域，特别是指一种核反应堆热工水力模拟软件cfd并行处理方法。

背景技术：

当前，以通用cpu加专用加速核心构成的异构新型处理器已成为超级计算机的发展趋势。神威·太湖之光是典型的基于异构新型处理器的超级计算机，它由40960块自主研发的申威26010众核处理器组成，每个处理器包含4个核组(codegroup，cg)，260个异构核心。如图1所示，单个核组由1个控制核心/主核(manageprocessingelement，mpe)和运算核心/从核(computingprocessingelement，cpe)阵列(8×8网格结构)组成，每个从核有16kb的一级指令缓存和64kb的局存，mpe和cpe共享8gb主存。也就是说：主存是主核和从核共享的8g内存，局存是从核独享的内存，只有64kb。

核反应堆热工水力学是研究核反应堆及其回路系统中冷却液的流动和热量传输特性的工程性学科，其内容涉及反应堆实际和假设的稳态和瞬态、事故出现的复杂热工水力现象。由于全尺度实验的限制，需要利用数值计算方法来分析所发生的复杂情况。计算流体力学(computationalfluiddynamics，cfd)数值模拟方法克服了传统实验测量方法的条件限制，被广泛用于核反应堆热工水力分析中，比如堆芯棒束通道内的流动换热过程研究和蒸汽发生器关键部件的流动传热研究。采用cfd进行热工水力分析，计算和存储需求大，比如采用直接数值模拟，计算量为雷诺数的三次方量级因此需要超级计算机的支撑。

采用cfd进行热工水力模拟要处理的基础是纳维斯托克斯(naiver-stokes，n-s)方程，对于非定常不可压缩n-s方程，可以采用谱元法进行空间离散化，谱元法是一种结合有限元法通用性与谱方法的高精度的高阶加权残差技术，基本思想将计算区域划分成若干谱单元，在每个求解单元中利用拉格朗日-洛巴托-勒让德配置点(gauss-legendre-lobatto，gll)点进行展开，利用投影算子技术和高斯数值积分得到离散方程。接着，将离散后的n-s方程拆分为对流项、压力项、速度项，对流项可以用显式外推法求解，而压力项和速度项需要隐式迭代求解，压力项和速度项都可以直接归结为公式(1)：

hu＝(h1k+h2m)u＝f(1)

其中，h1、h2为系数，h为亥姆霍兹算子，k是刚度矩阵，m是质量矩阵，u为待求解的变量域，f是拆解后压力泊松方程或速度亥姆霍兹方程的总合的右边项，当h1＝1，h2＝0即为压力的泊松方程形式，其他h1、h2系数为速度的亥姆霍兹方程形式。

基于谱元法离散后的cfd模拟软件以小矩阵乘为主要计算核心，但数学库和编译器不可能提供最好内核性能，使用通用的优化技术(如openacc、openmpi等)或优化库(如blas、gemm等)也会导致系统性能大幅度下降，所以要针对超算体系结构提供专门优化方案。

技术实现要素：

本发明要解决的技术问题是提供一种核反应堆热工水力模拟软件cfd并行处理方法，以解决现有技术所存在的系统性能下降的问题。

为解决上述技术问题，本发明实施例提供一种核反应堆热工水力模拟软件cfd并行处理方法，包括：

判断热工水力模拟软件cfd中矩阵乘中n2的取值，其中，所述矩阵乘表示为：a、b、c都表示双精度浮点矩阵，n1、n2、n3都表示矩阵维度的大小，n＝n2，n表示谱单元阶数；

若12≤n2≤24，则判断n1是否等于n²，若n1＝n²，则对矩阵a按照从核数m进行数据切分，将切分后的数据分配到m个从核的局存中，将矩阵b完整读入到m个从核的局存中，每个从核根据局存中的数据执行矩阵乘计算任务，任务完成后，将计算结果分配回矩阵c对应主存的地址中。

进一步地，所述对矩阵a按照从核数m进行数据切分，将切分后的数据分配到m个从核的局存中，将矩阵b完整读入到m个从核的局存中，每个从核根据局存中的数据执行矩阵乘计算任务，任务完成后，将计算结果分配回矩阵c对应主存的地址中包括：

将矩阵a矩阵转置为a^t；

确定用于执行矩阵乘计算的m个从核；

调用从核线程组，每个从核按照步骤a1-a5执行矩阵乘计算任务；其中，步骤a1-a5为：

a1，获取当前从核的编号id；

a2，根据公式确定平均每个从核分配的行数，其中，m表示从核的数目，ra为矩阵a的行数；

a3，判断当前从核id是否满足：id<(leave＝ra％m)，若满足，则当前从核需要获取矩阵a的数据量为偏移量为当前从核申请大小的局存空间以为基址取数据放入局存空间中，其中，leave表示剩余量，la为矩阵a的列数；

a4，当前从核申请slabb大小的局存空间spaceb，以b+offsetb为基址取slabb大小数据放入局存空间spaceb中，其中，slabb表示每个从核需要获取的矩阵b的数据量，slabb＝8rblbbyte，rb表示矩阵b的行数，lb表示矩阵b的列数；offsetb表示矩阵b的偏移量，offsetb＝0byte；

a5，当前从核申请大小为的局存空间spacec，当前从核根据局存空间spaceb中的数据进行矩阵乘计算，将计算结果存入局存空间spacec中。

进一步地，在当前从核根据局存空间spaceb中的数据进行矩阵乘计算时，若n2＝12,16,20,24，则利用从核的单指令多数据流扩展部件，采用对界单指令多数据流方式将访问连续的局存空间数据装入向量寄存器进行矩阵乘计算，将计算结果存入局存空间spacec中。

进一步地，所述任务完成后，将计算结果分配回矩阵c对应主存的地址中包括：

任务完成后，根据id将局存空间spacec中的结果以直接内存存取方式向矩阵c对应的主存的基地址c+offsetc跨stridec步存入数据块bsize，直到spacec回传完毕，释放所有局存；其中，offsetc表示矩阵c的偏移量，stridec表示回传跨步大小，bsize表示跨步向量块大小，

进一步地，所述方法还包括：

若12≤n2≤24且n1≠n²，则判断n3是否等于n²，若n3＝n²，则将矩阵a完整读入到m个从核的局存中，对矩阵b按照从核数m进行数据切分，将切分后的数据分配到m个从核的局存中，每个从核根据局存中的数据执行矩阵乘计算任务，任务完成后，将计算结果分配回矩阵c对应主存的地址中。

进一步地，所述方法还包括：

当12≤n2≤15时，n1≠n²且n3≠n²，则采用循环展开方式进行矩阵乘计算，得到矩阵乘计算结果。

进一步地，所述方法还包括：

当16≤n2≤24时，n1≠n²且n3≠n²，则采用非对界单指令多数据流方式在主核中进行向量化计算，得到矩阵乘计算结果。

进一步地，所述方法还包括：

若1≤n2≤11，则采用循环展开方式进行矩阵乘计算。

本发明的上述技术方案的有益效果如下：

上述方案中，基于申威众核架构，提供一种以矩阵乘为计算核心的cfd软件的并行优化处理技术，以充分利用申威众核处理器强大的计算能力，对模拟软件cfd进行并行优化处理，能够大大降低cfd模拟时间，提高程序运行的性能，从而提高热工水力模拟效率。

附图说明

图1为本发明实施例提供的申威众核处理器架构的结构示意图；

图2为本发明实施例提供的向量化的工作原理示意图；

图3为本发明实施例提供的核反应堆热工水力模拟软件cfd并行处理方法的流程示意图；

图4为本发明实施例提供的循环展开矩阵乘的工作原理示意图；

图5为本发明实施例提供的第一种从核优化方案的工作原理示意图；

图6为本发明实施例提供的第二种从核优化方案的工作原理示意图；

图7为本发明实施例提供的核反应堆热工水力模拟软件cfd并行处理方法的详细流程示意图；

图8为本发明实施例提供的性能效果对比示意图。

具体实施方式

为使本发明要解决的技术问题、技术方案和优点更加清楚，下面将结合附图及具体实施例进行详细描述。

本发明针对现有的系统性能下降的问题，提供一种核反应堆热工水力模拟软件cfd并行处理方法。

为了更好地理解本发明实施例提供的核反应堆热工水力模拟软件cfd并行处理方法，对术语：循环展开、向量化、列优先、athread库、字节对齐、数据切分和消息传递接口库进行简要说明：

1)循环展开

由于矩阵乘中的n2的值是可以确定的，传统的矩阵乘法采用三层循环计算，将确定值的一层展开为连续的乘加操作，就是循环展开(见algorithm1)：

该操作能充分使用mpe的流水线性能。

2)向量化

如图2所示，向量化计算是一种特殊的并行计算，可以在同一时间对不同数据执行同样的一个指令，或者说指令应用于一个数组中，在mpe或cpe提供有256位的单指令多数据流(singleinstructionmultipledata，simd)扩展部件，所以在mpe或cpe中可同时计算8个整型(32位)或4个双精度浮点数(64位)。

3)列优先

矩阵(或数组)是按列在内存中连续存储的，如一个二维数组t(10,10),若按列优先就是按照t(1,j)、t(2,j)、t(3,j)…t(10,j)连续在内存中存储，而行优先则按照t(i,1)、t(i,2)、t(i,3)…t(i,10)连续在内存中存储。

4)athread库

athread库是调用从核操作相关的函数库，athread库中的函数以athread_开头。

5)字节对齐

现代计算机中，内存空间按照字节划分，在访问特定类型变量时经常在特定的内存地址访问，这就需要各种类型数据按照一定的规则在空间上排列，而不是顺序一个接一个地存放，这就是对齐。一般双精度浮点数按照32字节地址对齐(变量所在地址可以整除32)，整型按照4字节地址对齐(变量所在地址可以整除4)等。

6)数据切分

数据切分是将一个矩阵分成多个数据片段(slab)，每一个slab有一定的数据量大小。

7)mpi(messagepassinginterface，消息传递接口)库

mpi库是一种消息传递编程模型函数库，函数以mpi_开头。

如图3所示，本发明实施例提供的核反应堆热工水力模拟软件cfd并行处理方法，包括：

判断热工水力模拟软件cfd中矩阵乘中n2的取值，其中，所述矩阵乘表示为：a、b、c都表示双精度浮点矩阵，n1、n2、n3都表示矩阵维度的大小，n＝n2，n表示谱单元阶数；

若12≤n2≤24，则判断n1是否等于n²，若n1＝n²，则对矩阵a按照从核数m进行数据切分，将切分后的数据分配到m个从核的局存中，将矩阵b完整读入到m个从核的局存中，每个从核根据局存中的数据执行矩阵乘计算任务，任务完成后，将计算结果分配回矩阵c对应主存的地址中。

本发明实施例所述的核反应堆热工水力模拟软件cfd并行处理方法，基于申威众核架构，提供一种以矩阵乘为计算核心的cfd软件的并行优化处理技术，以充分利用申威众核处理器强大的计算能力，对模拟软件cfd进行并行优化处理，能够大大降低cfd模拟时间，提高程序运行的性能，从而提高热工水力模拟效率。

本实施例中，以基于谱元法的cfd模拟软件nek5000为例，在申威26010众核处理器上进行优化处理。本实施例中，谱元法离散求解n-s方程，主要是迭代计算隐式项(即式(1))，在nek5000中，每次迭代都需要通过下面例程构造hu算子：

该例程主要由矩阵乘组成，矩阵乘大小对应上面例程的matrix-multiplysize(行维度，列维度，内部维度)，p为谱单元阶数，p的范围为[1,24]，其中h、k、h1、h2、m、u的意义同公式(1)，g是由空间方向的几何项和下标导出的常数对角矩阵。

采用algorithm1得到的nek5000的各子例程的运行时间如表1所示，其中64.98％的时间花费在计算矩阵乘的例程mxfn，n的范围为[1,24]，矩阵乘具有以下形式：

cn1×n3＝an1×n2bn2×n3(2)

其中，a、b、c都表示双精度浮点矩阵，n1、n2、n3都表示矩阵维度的大小，一般情况下：n1、n2、n3的大小不超过524，而且最主要出现下面三种的矩阵大小情况：

1)n1＝n，n2＝n，n3＝n

2)n1＝n²，n2＝n，n3＝n

3)n1＝n，n2＝n，n3＝n²

表1例程的运行时间表

表1中，name表示例程名称，time(％)表示例程占用总运行时间的百分比；cumulativeseconds表示累积时间，即表1中，selfseconds累积和；selfseconds表示运行相应例程时的耗时；calls表示例程调用次数；selfs/call调用例程一次的平均耗时；totals/call表示调用例程一次的平均耗时，包括例程中的子函数。

本实施例中，判断热工水力模拟软件cfd中矩阵乘中n2的取值具体可以包括以下步骤：

首先调用mpi_init函数启动后，调用athread_init函数完成加速线程库初始化；

读入矩阵乘的矩阵维度大小n1、n2、n3，双精度浮点矩阵a、b、c；

判断n2的大小，根据n2的值，选择优化后的mxfn例程(mxfn中的n＝n2)。

本实施例中，若1≤n2≤11，直接采用图4所示的循环展开方式进行矩阵乘计算，提高计算速度，无需采用辅助优化技术。

在前述核反应堆热工水力模拟软件cfd并行处理方法的具体实施方式中，进一步地，所述对矩阵a按照从核数m进行数据切分，将切分后的数据分配到m个从核的局存中，将矩阵b完整读入到m个从核的局存中，每个从核根据局存中的数据执行矩阵乘计算任务，任务完成后，将计算结果分配回矩阵c对应主存的地址中包括：

将矩阵a矩阵转置为a^t；

确定用于执行矩阵乘计算的m个从核；

调用从核线程组，每个从核按照步骤a1-a5执行矩阵乘计算任务；其中，步骤a1-a5为：

a1，获取当前从核的编号id；

a2，根据公式确定平均每个从核分配的行数，其中，m表示从核的数目，ra为矩阵a的行数；

a3，判断当前从核id是否满足：id<(leave＝ra％m)，若满足，则当前从核需要获取矩阵a的数据量为偏移量为当前从核申请大小的局存空间以为基址取数据放入局存空间中，其中，leave表示剩余量，la为矩阵a的列数；

a4，当前从核申请slabb大小的局存空间spaceb，以b+offsetb为基址取slabb大小数据放入局存空间spaceb中，其中，slabb表示每个从核需要获取的矩阵b的数据量，slabb＝8rblbbyte，rb表示矩阵b的行数，lb表示矩阵b的列数；offsetb表示矩阵b的偏移量，offsetb＝0byte；

a5，当前从核申请大小为的局存空间spacec，当前从核根据局存空间spaceb中的数据进行矩阵乘计算，将计算结果存入局存空间spacec中。

本实施例中，对于n1＝n²，12≤n2≤24的矩阵乘，假设从核数为m，采用图5所示的从核优化方案进行优化，具体包括以下步骤：

首先，由于模拟软件cfd所采用的语言是按列优先存储，先将矩阵a矩阵转置为a^t；

然后，调用athread_spawn函数指定用于执行矩阵乘计算的m个从核；

最后，调用athread_join函数显式阻塞调用从核线程组，执行矩阵乘计算任务，具体的：

调用athread_get_id函数获取当前从核的编号，假设编号为id；

假设矩阵a的行为ra，列为la，为平均每个从核分配的行数，若当前从核id<(leave＝ra％m)，leave表示剩余量，那么则当前从核需要获取矩阵a的数据量为偏移量为首先需要通过调用ldm_malloc函数申请大小的局存空间然后调用athread_get函数使用单运算核心模式，运算核心私有连续段以为基址取数据放入局存空间中；

假设矩阵b的行为rb，列为lb，获取的矩阵b的数据量为slabb＝8rblbbyte，偏移量为offsetb＝0byte，当前从核需要完整读入矩阵b，通过调用ldm_malloc函数申请slabb大小的局存空间spaceb，调用athread_get函数以b+offsetb为基址取slabb大小数据放入局存空间spaceb中；

假设矩阵c的行为rc，列为lc，同样当前从核需要申请大小为的局存空间spacec，当前从核根据局存空间spaceb中的数据进行矩阵乘计算，将计算结果存入局存空间spacec中，具体的：

在当前从核根据局存空间spaceb中的数据进行矩阵乘计算时，可以按照图4所示进行循环展开计算(n2≠12,16,20,24)，并将计算结果存入spacec中；但是，若n2＝12,16,20,24，由于在从核申请内存后，局存空间以32字节对齐，则可以利用从核提供的256位单指令多数据流(simd)扩展部件，采用对界单指令多数据流方式将访问连续的内存区域装入向量寄存器进行矩阵乘计算，提高计算效率，并将计算结果存入spacec中；

任务完成后，调用athread_put函数使用单运算核心模式，根据id将局存空间spacec中的结果以直接内存存取(directmemoryaccess，dma)方式向矩阵c对应的主存的基地址c+offsetc跨stridec步存入数据块bsize，直到spacec回传完毕，调用ldm_free函数释放所有局存；其中，offsetc表示矩阵c的偏移量，stridec表示回传跨步大小，bsize表示跨步向量块大小，

在前述核反应堆热工水力模拟软件cfd并行处理方法的具体实施方式中，进一步地，所述方法还包括：

若12≤n2≤24且n1≠n²，则判断n3是否等于n²，若n3＝n²，则将矩阵a完整读入到m个从核的局存中，对矩阵b按照从核数m进行数据切分，将切分后的数据分配到m个从核的局存中，每个从核根据局存中的数据执行矩阵乘计算任务，任务完成后，将计算结果分配回矩阵c对应主存的地址中。

本实施例中，对于12≤n2≤24，n3＝n²的矩阵乘，假设从核数为m，采用图6所示的从核优化方案进行优化，具体包括以下步骤：

首先，由于模拟软件cfd所采用的语言是按列优先存储，先将矩阵a矩阵转置为a^t；

然后，调用athread_spawn函数指定用于执行矩阵乘计算的m个从核；

最后，调用athread_join函数显式阻塞调用从核线程组，执行矩阵乘计算任务，具体的：

调用athread_get_id函数获取当前从核的编号，假设编号为id；

假设矩阵b的行为rb，列为lb，为平均每个从核分配的列数，若当前从核id<(leave＝lb％m)，leave表示剩余量，那么则当前从核需要获取b的数据量为偏移量为首先需要通过调用ldm_malloc函数申请大小的局存空间spaceb，然后通过调用athread_get函数以为基址取数据放入局存空间spaceb中；

假设矩阵a的行为ra，列为la，获取的矩阵a的数据量为slaba＝8ralabyte，偏移量为offseta＝0byte，当前从核需要完整读入矩阵a，通过调用ldm_malloc函数申请slaba大小的局存空间spacea，通过调用athread_get函数以a^t+offsetb为基址取slaba数据放入spacea中；

假设矩阵c的行为rc，列为lc，需要当前从核申请的局存空间spacec；当前从核根据局存空间spaceb中的数据进行矩阵乘计算，将计算结果存入局存空间spacec中，具体的：

在当前从核根据局存空间spaceb中的数据进行矩阵乘计算时，可以按照图4所示进行循环展开计算，并将计算结果存入spacec中；但是，若n2＝12,16,20,24，同样使用从核的simd扩展部件优化加速计算，提高计算效率，并将计算结果存入spacec；

任务完成后，需要将spacec发送回矩阵c对应主存的地址中，这次就不需要跨步回传，假设矩阵c的偏移量为因为回传数据和对应位置是连续的，这时直接以dma方式往主存的基址c+offsetc回传spacec的所有数据，spacec回传完毕，调用ldm_free函数释放所有局存。

本实施例中，每个从核仅有64kb的局存，所以spacea+spaceb+spacec<64kb为宜。

在前述核反应堆热工水力模拟软件cfd并行处理方法的具体实施方式中，进一步地，所述方法还包括：

当12≤n2≤15时，n1≠n²且n3≠n²，则采用循环展开方式进行矩阵乘计算，得到矩阵乘计算结果。

本实施例中，若12≤n2≤15，n1≠n²且n3≠n²，则不采取任何优化措施，直接使用图4所示的循环展开方式进行矩阵乘计算。

在前述核反应堆热工水力模拟软件cfd并行处理方法的具体实施方式中，进一步地，所述方法还包括：

当16≤n2≤24时，n1≠n²且n3≠n²，则采用非对界单指令多数据流方式在主核中进行向量化计算，得到矩阵乘计算结果。

本实施例中，当16≤n2≤24时，n1≠n²且n3≠n²，由于涉及的双精度浮点数并不是按照32位字节对齐，无法采用对界simd方法，所以只能使用非对界simd方式在主核(mpe)中进行向量化计算，得到矩阵乘计算结果。

为了更好地理解本发明实施例提供的核反应堆热工水力模拟软件cfd并行处理方法，结合图7对所述方法的工作流程进行整体说明：

读入矩阵乘的矩阵维度大小n1、n2、n3，双精度浮点矩阵a、b、c；判断矩阵乘中n2的取值；若1≤n2≤11，则不进行优化，直接采用循环展开方式进行矩阵乘计算；若12≤n2≤24，则判断n1是否等于n²，若n1＝n²，则执行图5所示的第一种从核(cpe)优化方案；若n1≠n²，则继续判断n3是否等于n²，若n3＝n²，则执行图6所示的第二种从核(cpe)优化方案；否则，即：当16≤n2≤24时，n1≠n²且n3≠n²，使用主核(mpe)向量化的方法进行优化。

本实施例提供的面向申威众核架构的核反应堆热工水力模拟软件cfd并行处理方法，对于12≤n2≤24的矩阵乘按照相应情况，可以采取cpe上的从核优化方案，或mpe上的向量化优化技术，如图8所示，本实施例提供的核反应堆热工水力模拟软件cfd并行处理方法在申威众核架构上有近30％的性能提升(如图8)，大大缩减了运行时间。另需要说明的是：本实施例提供的面向申威众核架构的核反应堆热工水力模拟软件cfd并行处理方法，比采用纯mpe上的向量化优化技术相比，约有20％性能提升。

以上所述是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明所述原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。