CN109542406A - 用于模式开发的并行求解方法和系统 - Google Patents

Abstract

本发明对于模式方程提供了高效的并行数值求解方法和系统，将与待求解的模式方程对应的离散形式的表达式，该表达式的操作数为以三维数组形式表示的物理量，该表达式的运算符包括预定义的模式运算符和被重载为支持三维数组的基本运算符；根据可用的并行进程数量和三维数组的维度来划分每个物理量对应的三维数组，并保存每个三维数组在各进程间的数据分布信息；构建与所述表达式对应的表达式图，并基于该表达式图对表达式进行求解计算。该并行求解方法使用户从模式公式快速构建出并行代码，屏蔽复杂且繁琐的并行程序设计细节，简化了模式程序开发的难度。

Description

用于模式开发的并行求解方法和系统

技术领域

本发明涉及领域特定语言和并行计算技术，尤其涉及用于模式开发的并行求解方法和系统。

背景技术

研究气候变化、地球系统演变、海洋演变等的主要手段是对这些系统进行数值模拟，这个过程也可以称为模式计算。诸如地球系统模式、海洋模式等模式的发展正逐步朝着更高分辨率和更多物理过程的方向发展。分辨率的提高对于模式计算提出了严峻的挑战，串行的模式程序难以满足需求。因此要求领域研究者(也可以称为模式开发人员)开发在超级计算机上运行的并行程序以满足诸如大气、海洋预报等高分辨率、高时效性等要求。以海洋为例，通过海洋模式对海面温度和高度的分布、海流分布以及海洋气候演变过程等进行数值模拟。利用海洋模式进行数值模拟本质上就是对描述海洋相关物理过程的原始模式方程进行数值求解的过程。一般而言，海洋模式的原始方程通常有7个，其中包括3个空间方向的动量守恒方程、1个描述海水质量或体积守恒的方程、1个描述海水密度与温度、盐度、压力之间关系的状态方程、以及描述海水位温和盐度守恒方程。完整的模式求解可包括针对原始模式方程设计相应的数值离散方程、实现该数值离散方程的串行程序和实现对应的并行程序等几个基本阶段。因此模式开发人员不仅需要关注其领域内科学问题的研究，还需要熟悉并行程序开发。然而，并行程序的开发非常复杂，不仅需要考虑节点之间通信等底层细节，还需要面对硬件平台多样化带来的挑战，这大大增加了模式程序开发的难度。

发明内容

因此，本发明的目的在于克服上述现有技术的缺陷，提供一种用于模式开发的并行求解方法和系统来降低模式开发的难度，以使模式开发人员和用户能从海量代码和复杂的“并行沼泽”中脱离出来，专注于领域科学问题的研究。

本发明的目的是通过以下技术方案实现的：

一方面，本发明提供了一种用于模式开发的并行求解方法，包括：

获取与待求解的模式方程对应的离散形式的表达式，该表达式的操作数为以三维数组形式表示的物理量，该表达式的运算符包括预定义的模式运算符和被重载为支持三维数组的基本运算符；

根据可用的并行进程数量和三维数组的维度来划分每个物理量对应的三维数组，并保存每个三维数组在各进程间的数据分布信息；

构建与所述表达式对应的表达式图，该表达式图中的数据节点对应表达式中的操作数，该表达式图中运算节点对应表达式中的运算符；

在每个并行进程中执行下列操作：

对于所构建的表达式图中的每个节点，确定其数据类型；

对于表达式图中具有同一父节点的两个数据节点，根据各三维数组的数据分布信息对这两个数据节点进行数据对齐；

基于该表达式图对表达式进行求解计算。

在一个实施例中，预定义的模式运算符可以包括有限差分、有限插值、梯度、散度、旋度、拉普拉斯、积分、投影中的一个或多个。

在一个实施例中，对于每个全局数组，在所有可用的并行进程中可以保存一份相同的数据分布信息。

在一个实施例中，在每个并行进程中执行的操作还可包括：对于所构建的表达式图中的每个节点，确定其边界区域，在所述边界区域中预先保存相邻进程中部分数据的副本。

在一个实施例中，基于表达式图对表达式进行求解计算可以包括：

生成表达式图的逆波兰式并对其进行哈希；

基于所述哈希的结果提取与该表达式图对应的融核函数来对表达式进行求解；所述融核函数中包括对于表达式中各三维数组的每个数据元素进行相应表达式计算以得到最终求解结果的源代码。

在一个实施例中，该方法还可包括响应于基于所述哈希的结果未找到与该表达式图对应的融核函数，动态生成与所述融核函数，其中所述融核函数是通过下列步骤生成的：

生成表达式图的逆波兰式并对其进行哈希，所述哈希的结果用于标识与该表达式图对应的融核函数；

对于表达式图中的每个节点确定其描述信息，所述描述信息包括数据类型、数据分布及边界区域信息；

生成与该表达式图对应的融核函数的源代码，并进行动态编译，将编译后的融核函数保存至融核函数表中，以所述哈希的结果作为该融核函数的唯一标识。

在一个实施例中，根据各三维数组的数据分布信息对这两个数据节点进行数据对齐可以包括：

根据各三维数组的数据分布信息检测两个数据节点的数据分布是否一致；如果不一致，则将其中一个节点的数据转移至与另一节点具有相同数据分布的临时数组中，以方便后续进行数组的逐元素操作。

在一个实施例中，对于所构建的表达式图中的每个节点确定其边界区域可以包括：

在表达式图中自上而下叠加各运算符节点预先设置的边界区域宽度，以确定表达式图中各个数据节点对应的边界区域宽度；

根据每个数组节点的边界区域宽度从相邻进程中获取对应宽度的数据并保存在该数组节点的边界区域中。

在一个实施例中，所构建的表达式图可以采用二叉树的形式。

又一方面，本发明提供了一种用于模式开发的并行求解系统，包括：

输入模块，用于获取与待求解的模式方程对应的离散形式的表达式，该表达式的操作数为以三维数组形式表示的物理量，该表达式的运算符包括预定义的模式运算符和被重载为支持三维数组的基本运算符；

分布管理模块，用于根据可用的并行进程数量和三维数组的维度来划分每个物理量对应的三维数组，并保存每个三维数组在各进程间的数据分布信息；

表达式图构建模块，用于构建与所述表达式对应的表达式图，该表达式图中的数据节点对应表达式中的操作数，该表达式图中运算节点对应表达式中的运算符；

表达式求解模块，用于在每个并行进程中执行下列操作：

对于所构建的表达式图中的每个节点，确定其数据类型；

对于表达式图中具有同一父节点的两个数据节点，根据各三维数组的数据分布信息对这两个数据节点进行数据对齐；

基于该表达式图对表达式进行求解计算。

与现有技术相比，本发明的优点在于：

对于模式方程提供了高效的并行数值求解方法，使用户能从模式公式快速构建出并行代码，屏蔽复杂且繁琐的并行程序设计细节，简化了模式程序开发的难度。同时，基于本发明开发的模式可以方便快速地在多种硬件平台进行移植部署，且与传统方式开发的模式运算速度相当。

附图说明

以下参照附图对本发明实施例作进一步说明，其中：

图1为根据本发明实施例的用于模式开发的并行求解系统的结构示意图；

图2为根据本发明实施例的三维数组在多个进程上的数据分布示意图；

图3为根据本发明实施例的三维数组在多个进程上的数据分布示意图；

图4(a)-(c)为根据本发明实施例的数组的数据分布及数据对齐示意图；

图5为根据本发明实施例的不同形状的数组的数据分布对齐示意图；

图6为根据本发明实施例的节点数据边界区域更新示意图；

图7为根据本发明实施例的表达式图的示意；

图8为根据本发明实施例的基于融核函数的表达式求解的流程示意图；

图9为根据本发明实施例的数组的数据转移和数据对齐示意图；

图10为根据本发明实施例的表达式各节点边界区域更新示意图；

图11为根据本发明实施例的表达式求解的流程示意图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图通过具体实施例对本发明进一步详细说明。应当理解，此处所描述的具体实施例仅用以解释本发明，并不用于限定本发明。

针对模式方程的求解中发生的计算通常是以空间上的网格点为核心进行的，需要不断对空间点上的物理量进行各种循环、差分和插值等操作，对每一个点上的物理量都要进行精确控制才能完成模式方程的求解。以使用Fortran语言开发为例，以传统方式实现某个模式方程的求解通常需要编写几千乃至几万行串行代码，如果改写成并行代码量通常增加3-4倍。模式方程越复杂，对应的模式代码量越大，开发和维护的成本也越高。

在当前的模式开发过程中采用从底向上的方式设计并行的模式程序比较繁琐和低效。在超级计算机中，数据分布在不同节点和进程内，模式开发人员需要通过消息通信和共享内存的方式精心安排计算进程和数据分布，处理复杂频繁的边界通信和协调各个进程的工作，这对于领域模式开发人员的编程能力提出了很高的要求。另外，硬件计算平台多样化也给并行程序的设计和优化带来了很大的挑战，往往要求程序设计者对这些平台的架构细节进行深入的了解。

针对上述问题，本申请提供了用于模式开发的高效并行求解方法和系统，将领域的具体模式开发工作与并行程序的设计分离，支持模式方程组的高效并行数值求解，达到“代码即方程”的效果，大幅度减少模式程序的代码量和设计复杂度，帮助模式开发人员从繁琐的并行程序设计细节中解脱出来，从而专心于领域科学问题本身的研究。该方法将模式原始方程转换为由算子表述的离散方程形式，为模式开发人员提供诸如差分、插值等基础并行算子，以及梯度、散度、旋度、拉普拉斯、积分、投影等高阶的并行算子。这些并行算子也可以称为模式运算符，以与基本数学运算符相区分，通过这些并行算子屏蔽底层复杂的并行、循环、数据通信等过程。以海洋模式中常用的二维连续方程为例，其连续的形式为下面的公式(1)：

其中D表示表示海表面到海底这个水柱的高度，η表示表示海表面相对海平面的高度，U和V分别为x和y方向上的水平速度，t表示时间步长。该模式的连续方程经数值离散后可表示为如公式(2)的形式：

其中公式(2)为公式(1)的离散形式，描述了上一个时刻和下一个时刻各个变量的关系。其中下标b表示上一个时刻，下标f表示下一个时刻。在现有的模式开发中，该公式(2)可以如下的伪码来实现：

其中eta_f对应公式(2)中的η_f，eta_b对应η_b，dt对应Δt，d对应D，u对应U，v对应V。可以看出在上述伪码中需要利用循环对海洋中每个网格点进行这样的运算。

与现有的模式开发不同，本申请使用全局分布式数组表示每个物理量并提供了可直接应用于该数组的一组模式运算符，从而可使用这些模式运算符将公式(2)对应的伪码实现为如下的形式：

Eta_f＝Eta_b-dt＊(DXF(AXB(D)＊U)+DYF(AYB(D)＊V))

其中DXF(即δ^x)指示在x方向上的前向差分、AXB(D)(即)指示x方向上的后向插值；DYF(即δ^y)指示y方向上的前向差分，AYB(D)(即)表示y方向上的后向插值。可以看出，上述伪码的形式与以数值离散形式表示的公式(2)相一致，从而可以使得领域研究者看到方程就能写出代码，而不需要考虑后台多节点间的消息传递、数据对齐等并行技术细节。

图1给出了根据本发明一个实施例的用于模式开发的并行求解系统的结构示意图。该系统结构可以被分为四层：在用户接口层向模式开发人员提供可使用的数据类型和函数，对模式开发人员屏蔽了底层多节点之间的数据管理和并行计算等细节。在表达层，使用延时求解的方法根据模式开发人员所写的代码构建和优化表达式图并对其进行求解。在分布式数据管理层负责分布式数据的描述和管理及性能优化，以使得用户可直接操作分布在不同计算节点上的数组。数学计算层提供求解表达式所需的内核函数和融合内核函数(融核函数)。下面将结合图1具体介绍上述四层的技术细节。

1.接口层

为了简化编程，在本申请的实施例中，模式方程中所处理的物理量都采用三维分布式数组(在下文中记为GArray)的形式来表达，该GArray可以表示在多个MPI进程之间分布的多维数组，例如一维数组、二维数据、三维数组。与传统的数组不同，在本发明的实施例中的数组GArray中的数据可被分布式地存储在多个节点上。例如物理量“温度”，可以用GArray类型的变量T表示，这个数组中每一个点描述的就是各个位置海水的温度，但T存储为分布式的，从而更有利于利用现有的分布式架构的超级计算机进行运算。这样的数组GArray也可被称为分布式数组或全局数组。该数组GArray的数据类型可以是整数(int)、浮点数(float)或双精度浮点数(double)等。对于这些数据类型的大部分常见的运算符号都被重载以兼容该GArray类型。由于采用了这种全局数组的形式表示模式方程中的物理量，使得模式求解不再是对三维空间中的单个点，而是针对整个流场提供三维变量直接运算和处理，简化了模式方程及其代码的实现和维护成本。

在实际的分布式环境下，所创建的全局数组被剖分成小的数据块分散到各个MPI进程中，分别针对每个数据块进行存储和计算。每个全局数组都具有与其关联的数据分布信息，用于描述该全局数组在多个进程或节点之间的数据分布情况。每个全局数组的数据分布信息是对该全局数组进行操作或运算的基础信息。在下文分布式数据管理部分会更详细地介绍如何根据数据不同维度的比例和进程数据在x、y、z方向上合理剖分全局数组。

由于所操作的数据类型的改变，针对数据进行处理的相关函数，例如运算函数、输入和输出接口操作等都需要进行相应的重载以适应新的数据类型。表1给出了一些主要的函数接口示例，其中包括与数组创建、数组操作、输入/输出操作、数学运算符、模式运算符等相关的函数接口，这些重载后的函数接口都可以直接对分布式数组进行整体操作。以全局数组A和B的加法操作为例(即A+B)，重载的“+”函数中逐个将A中元素与B中相同位置的元素相加。另外重载的“+”函数中还支持加法的并行计算，其向调用该函数接口的上层引用隐藏了底层并行计算的细节，例如下文会详细介绍，检查全局数组A和B的数据分布是否对齐，如果没有对齐，则进行相应进程间数据转移，使二者数据对齐以保证加法操作的正确执行等。

模式开发本质是求解方程组，而这些方程组大多是偏微分方程的形式，方程中包含温度、速度、体积等物理量。进过离散化后可以转换为有限差分或插值计算等。也就是说模式代码中会对温度、速度、体积等各种物理量进行差分或插值等算子运算。为了简化模式开发人员编程的难度，该实施例为有限差分/插值计算提供了12个基本模式运算符(也称为算子)。算子用三个字母组合命名，形式为[A|D][X|Y|Z][F|B]，其中第一个字母A或D表示插值或差分算子。第二个字母X，Y或Z表示操作方向，最后一个字母B或F表示后向或前向操作。每个算子的数学表达形式如表2所示。差分算子还可利用物理变量绑定的网格信息，不仅具有逐差运算，还包括除法运算，实现了求导功能，使用物理变量的差值除以其所对应的空间格距，从而使程序代码与数学公式更加吻合。通过组合这些基本的算子可以构建大部分的模式方程的离散形式。因此，可将待求解的模式方程转换为这样的离散形式的数学表达式：该数学表达式由操作数和运算符组成，其中操作数为以三维数组形式表示的物理量，运算符可包括预先定义的模式运算符和被重载为支持三维数组数据类型的各种基本运算符。应理解，上述模式运算符仅是举例说明而非进行任何限制，可以根据具体领域的不同设置不同的模式运算符，例如梯度、散度、旋度、拉普拉斯、积分、投影等高阶并行的模式运算符。

表1

表2：算子运算符

上述函数接口都可以例如使用大部分模式开发人员都熟悉的Fortran语言的运算符重载功能来实现，与底层相关的并行计算细节被封装在具体函数的实现中。而对于模式开发人员，可遵循其已有的编程习惯，利用上述数据类型和相关函数接口实现对模式方程的直接编码，由此大幅度降低模式程序的开发难度。

2.分布式数据管理

如上文讨论的，在底层计算时会将要处理的数据分布在多个节点上并行进行运算。在本申请的一个实施例中，可以根据可用进程数量和数据不同维度的比例在x、y、z方向上对全局数组进行合理剖分来实现数据的分布。这样的分解信息可保存在与该数组对应的分区对象中。该分区对象在数组创建时创建，随着数组的运算动态更新，其保存的分布信息主要包括MPI进程的形状、数组形状和在x、y、z方向上的分段信息。图2给出了使用12个MPI进程分布9×6×6的数组的3D可视化示意图。如图2所示，全局数组GArray整个是三维立体的，可以根据进程个数把GArray剖分成多个小的立体，每个小立体对应一个进程。所以进程划分也可以看作是三维的。进程划分依据可设定为：尽量使得全局数组划分后的每个小立体接近立方体。图2中将MPI进程的形状设置为3×2×2，每个小分块表示1个MPI进程。表3给出了对于这样的分解的分布信息示意，其中记录MPI进程形状为3×2×2，数组对象形状为9×6×6。以x方向为例，在x方向上，数组大小为9，划分为3个进程，因此这3个进程上分得的数组大小分别为3、3、3(3+3+3＝9)，即对应在x方向上的分段为记为[333]。y、z方向以此类推。在y方向上的分段记为[3 3]，在z方向上的分段记为[3 3]。又例如，在形状为2×3×1的MPI进程上分布6×10×1的数组，其分解结果如图3所示。

表3

在又一个实施例中，对于给定的数组可以提取其子数组，以图3所示的数据对象为例，从其x方向提取第3至第5个元素并从其y方向提取第3至第7个元素，可得到如图4(a)所示的子数组对象B；类似地，从其x方向提取第2至第4个元素并从其y方向提取第6至第10个元素可得到如图4(b)所示的子数组对象C。这样提取的子数组并没有改变MPI进程的形状和底层的数据布局，因此也不会引发MPI通信。

在又一个实施例中，由于数组分布在不同的MPI进程上，当需要对数组进行逐元素操作时，需要确保待操作的数组之间数据对齐以便在每个进程中正确执行相应的数学计算。例如计算A＝B+C，以图4(a)和(b)所示的分区数组B和C为例，由于数组B和C在底层MPI进程上的数据布局不同，因此不能直接对B和C进行相加操作，需要将B和C的数据分布布局对齐才能进行相关运算。例如，可根据B的分布来对数组C重新分区，将数组C的数据进行转移，转移至与B有相同分区布局的数组C1中(如图4(c))，然后执行A＝B+C1。这样的数据对齐检测和数据转移对于用户是透明的，由系统通过不同进程间点对点传输、数据广播等操作来实现进程之间的数据转移。例如，可设置数据对齐检测函数和数据转移函数，当执行数组之间的操作时，隐式地自动调用这些函数进行相关的数据对齐和数据转移操作。

对于不同形状的数组之间的计算可以通过维度扩展来完成数据对齐。在模式计算中经常有很多二维变量与三维变量之间的计算，以图5所示的分布在12个MPI进行上的10×8的A数组和10×8×2的B数组为例，当要执行A与B之间的运算时，A的Z方向上部对应的MPI进程中没有数组A的部分，因此可以根据B的分布在Z方向上扩展数组A(例如通过补0的方式)得到临时扩展的三维数组A1来与B进行运算。

考虑到模式计算中要处理的海量数据可能分散地存储在各个节点上，而参考表2所示的模式运算符，在模式计算中常常需要与周围点的数据进行运算，这意味着在MPI进程间需要传递数据从而保证结果的正确性。因此，在又一个实施例中，为改善计算性能，避免频繁的数据传递，每个MPI进程可预先保存相邻进程的数据的副本，也就是在多个MPI进行上分解数组时，允许被划分的数组边界之间有一定的重复。并且还可提供专用的边界区域管理模块来自动更新和维护这些处于边界区域中的数据以提高计算性能。如图6所示，对于每个节点数据需要更新的边界区域可用两个字母[l|u][x|y|z]的组合来表示，其中第一个字母l或u表示下界或上界，x，y或z表示方向。例如，每个节点可使用两个向量l和u来记录要更新的边界区域的信息，每个向量l和u由三个元素构成，分别代表三个方向上的下限/上限的更新区域厚度。例如，如果要更新厚度为1，x方向上的下界区域，则将数据节点的l和u分别设置为[1,0,0]和[0,0,0]。

3.表达式求解

利用上述提到的数据类型及相关函数接口，模式开发人员可以直接按照模式方程的离散形式来开发相应的模式程序，在模式程序中每个方程以表达式的形式表示。该表达式中操作的数据对象以上文介绍的GArray的形式来表示，该表达式中的运算符可以为上文介绍的重载后支持GArray类型的运算符和新建的模式运算符。在表达式求解时，先针对表达式中的各元素建立节点并构建表达式图，所构建的表达式图例如可以采用二叉树的形式。图7示出了针对表达式C＝2.0*A+B构建的表达式图。对于所构建的表达式图可通过下列步骤来进行求解：

1)针对表达式图的各数据节点中全局数组，执行数据对齐检测和数据转移等操作，以为后续运算做好准备；

2)从下而上递归调用每个运算符节点对应的运算函数，最后得到表达式的计算结果；

3)清除上述求解过程中产生的中间变量并销毁所构建的表达式图。

在模式程序中，大多数解决偏微分方程的数值模拟都是按照时间步长循环的，表达式图的构建和销毁被频繁调用，总时间成本不可忽略。因此为了减少这种性能开销，在一些实施例中，该系统还包括节点池来缓存表达式图中的已创建的节点对象，以便以后可以重复使用。另外在表达式求解过程总会产生大量的临时数组对象，例如以A＝B+C+D为例，B+C的结果保存在临时数组T1中，然后执行计算T1+D并将最终结果传递给A。在每次循环中求解完之后都会销毁临时数组T1。临时的中间变量的频繁创建和销毁也会带来不必要的性能开销。由此优选地，该系统还包括数组池来缓存临时创建的数组对象，并为每个缓存的数组对象建立唯一标识，例如对通过全局数组的数据类型及其分布对象中保存的数据分布信息进行哈希得到该数组对象在数组池中的唯一标识。在需要临时数组时，就可以通过数据类型及所需的数据分布信息，在不分配内存空间的情况下直接从数组池中检索并获取临时数组。

为了进一步提高效率，在本发明的一些实施例中对表达式采用了延时求解的方式，重载了赋值运算符，仅通过赋值运算符来触发表达式的求解过程，而在其他重载的运算符中不设置求解计算功能，而是进行表达式相关节点的建立及数据准备工作。以表达式C＝2.0*A+B为例，通常的表达式求解是先计算C1＝2.0*A，然后再计算C＝C1+B。每次调用相关运算符函数都会执行求解的过程，这样会产生临时变量C1并且分两步计算花费时间更长。而本发明实施例中提供的延时求解的方法针对表达式中的各元素先建立节点并构建表达式图，在调用赋值符“＝”时才求解，这样不需要产生临时量，其一次求解速度较快。也就是说，本发明实施例的表达式求解引擎并非直接求解该表达式，而是以运算符为中间节点，数据为叶节点构建与该表达式对应的表达式图，仅在遇到赋值符“＝”时才开始求解。这样的延时求解功能是通过在重载的赋值运算符中调用与表达式图对应的融核函数来完成的。

在上文介绍的基于表达式图从下而上递归调用每个运算符节点对应的运算函数进行求解的实施例中，频繁调用与每个运算节点对应的内核函数不仅耗时而且产生了大量的临时变量，者可能会带来较大性能开销。针对这个问题，该实施例中提供基于融核函数的表达式求解方法，对于每个表达式，动态生成与整个表达式对应的内核函数并通过执行该内核函数来获取表达式的计算结果。由于这个表达式对应的内核函数融合了表达式中各运算符节点对应的内核函数，因此可以将这个函数称为融核函数。对整个表达式的求解可以通过调用这个融核函数一次就可以完成，由此减少了由于多次内核启动和频繁读取临时变量带来的性能开销。

图8给出了根据本发明实施例的基于融核函数的表达式求解方法的流程示意图。如图8所示，该方法主要包括下列步骤：

S1)构建与表达式对应的表达式图并生成与该表达式对应的哈希值。以图8中表达式A＝B-C+D为例，先构建其对应的表达式图，+、-为运算节点；然后通过递归表达式可以生成与该表达式图对应的逆波兰表达式，作为该表达式图专属的签名字符串，接着对该字符串进行哈希得到的哈希值可用于唯一地标识该表达式。同一表达式具有相同的签名字符串和哈希值，并且该哈希值还可用于关联与该表达式对应的融核函数。不同的表达式因为计算逻辑不同，对应不同的融核函数。

S2)确定表达式图中每个节点的数据类型。为了正确生成融核函数来求解表达式，必须确定每个节点的数据类型。表达式图中叶子节点的类型是已经确定的，根据叶子节点的数据类型与运算结点的类型进行数据类型推断，自下而上依次得到表达式图中每个中间节点的类型以及最后运算结果的类型。例如，单精度加双精度结果为双精度数浮点数，两个双精度浮点数进行比较运算，结果类型为整数。根据B，C的数据类型和运算符减号确定中间结果B-C的数据类型，再依次类推，得到最后结果的数据类型。

S3)确定表达式图中每个节点的分布信息。如上文讨论的，每个数组都有其对应的分区对象来保存相关的分布信息。当要对两个数组进行运算时，需要将这两个数组的数据分布对齐才能保证运算的正确性。以A＝B-C+D为例，如图9所示，如果D具有与B和C不同的数据分布，可以在调用核函数前隐式地调用相关的数据转移函数，将数组D转移为与B和C具有相同分布的D1，最后再进行求解。以此类推，根据子节点的数据分布与运算结点的类型进行数据分布推断，自下而上依次得到表达式图中每个中间节点的分布以及最后运算结果的数据分布信息。

S4)确定表达式图中每个节点的边界区域更新。如上文结合图6所介绍的，每个节点需要更新的边界区域可用两个字母[l|u][x|y|z]的组合来表示，其中第一个字母l或u表示下界或上界，x，y或z表示方向。例如，每个节点可使用两个向量l和u来记录要更新的边界区域的信息，每个向量l，u可由三个元素构成，分别代表三个方向上的下限/上限的更新区域厚度。例如，如果某个数据节点需要更新厚度为1，x方向上的下界区域，则可将该节点对应的边界更新向量l和u分别设置为[1,0,0]和[0,0,0]。对于诸如+、－、*、/之类的基本数据运算符节点，不需要更新边界区域，因此，其对应的边界更新向量l和u分别默认设置为[0,0,0]和[0,0,0]；如果是在表达式图中，则与普通数学运算符对应的节点默认与父节点的更新需求一致。对于上文提到的单个的模式运算符节点，比较容易确定其需要更新的边界区域。以AXB(A)为例，可以确定要更新数组A的x方向上的下界区域且更新厚度为1，因此对应AXB运算符的节点的边界更新向量为l和u分别设置为[1,0,0]和[0,0,0]。表4给出了对于单个模式运算符需要更新的边界区域信息。

表4

但对于复杂的表达式，其中很多运算符都涉及边界更新和通信，需要自上而下叠加各运算节点的通信边界效应，由此确定表达式图中每个数据节点里数组的边界更新区域。以上文提到的公式(1)对应的表达式E＝E-t*(DXF(AXB(D)*U)+DYF(AYB(D)*V))为例，图10中给出了该表达式对应的表达式图及其各节点边界更新情况的示意。通过自上而下从根节点到叶子节点的路径叠加更新需求。每个叶子节点再将多条路径的更新情况进行叠加，则得到每个数据节点中数组边界的更新情况。例如，如果第一层节点需更新u(1,0,0),第二层节点也需更新u(1,0,0)，则第三层叶节点叠加路径更新需求后是u(1,0,0)，叠加方式为取两者的最大值。例如，在图10中每个普通数学运算符(+、－、*等)对应的节点对于边界更新无影响，其边界更新向量l和u与其父节点保持一致；对于从上至下路径中首次出现的模式运算符对应的运算节点，可根据表4设置其边界更新向量l和u；对于其余节点，将其父节点的边界更新向量l和u合并至本节点的边界更新向量l和u中。

S5)针对上述表达式图生成融核函数的源代码，并进行动态编译，将编译后的融核函数保存至融核函数表中，以在步骤S1)中生成的与该表达式图对应的哈希值作为该融核函数的唯一标识。在程序执行过程中，当遇到该表达式时，准备好相应输入数据就可以调用该表达式对应的融核函数对其求解了。所生成的融核函数包括对于表达式中以全局数组形式表示的各变量的每个数据元素进行相应表达式计算以得到最终求解结果的源代码。融核函数中的源代码可以是以C或C++语言编写的源代码。关于所生成的融核函数的编译，可以自由地选择不同的编译器。例如，英特尔编译器、基于LLVM的定制的轻量级JIT编译器等。优选可使用英特尔编译器，因为它通常能比LLVM编译器在英特尔CPU上产生更高性能的代码。如果在计算节点上找不到英特尔编译器，会使用LLVM作为备份。编译后的融核函数存储在融核函数表中，每个融核函数在程序运行的生命周期中只需要生成一次，以后再碰到相同的表达式，只需要根据其哈希值直接从融核函数表中取出对应的融核函数并执行即可。

仍以上文求解海平面高度的公式(1)为例，采用上文介绍的数据类型和运算函数可以将经离散后的公式(1)的对应代码写为下列的表达式：

E＝E-t*(DXF(AXB(D)*U)+DYF(AYB(D)*V))。

图11(a)给出了基于该表达式构建的表达式图。图11(b)示出了自动生成的与该表达式图对应的源代码并由编译器将其编译成可执行的融核函数。如图11(b)所示，所生成的融核函数包括对于以全局数组形式表示的各变量的每个数据元素进行相应表达式计算从而得到最终表达式结果的源代码，其主体结构为对于三维数组中各元素进行操作的三重循环结构，每次循环都是将数组中相应数据元素代入具体表达式进行相应计算，直到遍历完所有的数据元素为止。图11(c)可以看出，与之前递归求解的实施例相比，基于融核函数的表达式求解方法只需要将输入数据提供给可执行的融核函数，进行一次函数调用，就可以得到表达式的最终结果，减少了多次函数调用带来的性能开销。

在本发明一个实施例中还提供了利用上述的用于模式开发的并行求解系统进行并行求解的方法，其包括：S1)获取与待求解的模式方程对应的离散形式的表达式，该表达式的操作数为以三维数组形式表示的物理量，该表达式的运算符包括预定义的模式运算符和被重载为支持三维数组的基本运算符；S2)根据可用的并行进程数量和三维数组的维度来划分每个物理量对应的三维数组，并保存每个三维数组在各进程间的数据分布信息；S3)构建与所述表达式对应的表达式图，该表达式图中的数据节点对应表达式中的操作数，该表达式图中运算节点对应表达式中的运算符；以及S4)在每个并行进程中执行下列操作：对于所构建的表达式图中的每个节点，确定其数据类型；对于表达式图中具有同一父节点的两个数据节点，根据各三维数组的数据分布信息对这两个数据节点进行数据对齐；基于该表达式图对表达式进行求解计算。

虽然本发明已经通过优选实施例进行了描述，然而本发明并非局限于这里所描述的实施例，在不脱离本发明范围的情况下还包括所做出的各种改变以及变化。

Claims (10)

1.一种用于模式开发的并行求解方法，包括：

获取与待求解的模式方程对应的离散形式的表达式，该表达式的操作数为以三维数组形式表示的物理量，该表达式的运算符包括预定义的模式运算符和被重载为支持三维数组的基本运算符；

根据可用的并行进程数量和三维数组的维度来划分每个物理量对应的三维数组，并保存每个三维数组在各进程间的数据分布信息；

构建与所述表达式对应的表达式图，该表达式图中的数据节点对应表达式中的操作数，该表达式图中运算节点对应表达式中的运算符；

在每个并行进程中执行下列操作：

对于所构建的表达式图中的每个节点，确定其数据类型；

对于表达式图中具有同一父节点的两个数据节点，根据各三维数组的数据分布信息对这两个数据节点进行数据对齐；

基于该表达式图对表达式进行求解计算。

2.根据权利要求1所述的方法，其中预定义的模式运算符包括有限差分、有限插值、梯度、散度、旋度、拉普拉斯、积分、投影中的一个或多个。

3.根据权利要求1所述的方法，其中对于每个全局数组，在所有可用的并行进程中保存一份相同的数据分布信息。

4.根据权利要求1所述的方法，其中在每个并行进程中执行的操作还包括：对于所构建的表达式图中的每个节点，确定其边界区域，在所述边界区域中预先保存相邻进程中部分数据的副本。

5.根据权利要求1所述的方法，其中基于表达式图对表达式进行求解计算包括：

生成表达式图的逆波兰式并对其进行哈希；

基于所述哈希的结果提取与该表达式图对应的融核函数来对表达式进行求解；所述融核函数中包括对于表达式中各三维数组的每个数据元素进行相应表达式计算以得到最终求解结果的源代码。

6.根据权利要求5所述的方法，还包括响应于基于所述哈希的结果未找到与该表达式图对应的融核函数，动态生成与所述融核函数，其中所述融核函数是通过下列步骤生成的：

生成表达式图的逆波兰式并对其进行哈希，所述哈希的结果用于标识与该表达式图对应的融核函数；

对于表达式图中的每个节点确定其描述信息，所述描述信息包括数据类型、数据分布及边界区域信息；

生成与该表达式图对应的融核函数的源代码，并进行动态编译，将编译后的融核函数保存至融核函数表中，以所述哈希的结果作为该融核函数的唯一标识。

7.根据权利要求1所述的方法，其中根据各三维数组的数据分布信息对这两个数据节点进行数据对齐包括：

根据各三维数组的数据分布信息检测两个数据节点的数据分布是否一致；如果不一致，则将其中一个节点的数据转移至与另一节点具有相同数据分布的临时数组中，以方便后续进行数组的逐元素操作。

8.根据权利要求4所述的方法，其中对于所构建的表达式图中的每个节点确定其边界区域包括：

在表达式图中自上而下叠加各运算符节点预先设置的边界区域宽度，以确定表达式图中各个数据节点对应的边界区域宽度；

根据每个数组节点的边界区域宽度从相邻进程中获取对应宽度的数据并保存在该数组节点的边界区域中。

9.根据权利要求1-8中任一项所述的方法，所构建的表达式图采用二叉树的形式。

10.一种用于模式开发的并行求解系统，包括：

输入模块，用于获取与待求解的模式方程对应的离散形式的表达式，该表达式的操作数为以三维数组形式表示的物理量，该表达式的运算符包括预定义的模式运算符和被重载为支持三维数组的基本运算符；

分布管理模块，用于根据可用的并行进程数量和三维数组的维度来划分每个物理量对应的三维数组，并保存每个三维数组在各进程间的数据分布信息；

表达式图构建模块，用于构建与所述表达式对应的表达式图，该表达式图中的数据节点对应表达式中的操作数，该表达式图中运算节点对应表达式中的运算符；

表达式求解模块，用于在每个并行进程中执行下列操作：

对于所构建的表达式图中的每个节点，确定其数据类型；

对于表达式图中具有同一父节点的两个数据节点，根据各三维数组的数据分布信息对这两个数据节点进行数据对齐；

基于该表达式图对表达式进行求解计算。