一种仿真模型的任务分解与并行求解的方法

1.本发明属于仿真与并行计算交叉技术应用领域，尤其涉及一种面向仿真模型的任务分解与并行求解的方法。


背景技术：

2.在面向建模仿真的很多工程领域中，如物理、化学、生物学等，都涉及大规模常微分方程(ode)初值问题(ivp)求解。由于采用离散化方法，在特定求解精度的约束下，会将局部步长设置较小，在每个仿真全局步长时间段内，进行成百上千次的累积迭代ode求解，消耗了大量的仿真时间。典型情况下，对于航天发动机、模拟燃烧、地下水传输等复杂模型，其求解ode的时间消耗会达到墙上时钟时间的90％以上，构成了仿真推进速度的瓶颈。如何利用现代gpu的多核能力，构建高效的模型并行求解方法，已经成为目前亟待解决的一个仿真并行计算的难题。
3.传统并行求解方法主要集中在三个方向：
4.一种是采用类似传输线模型的方法，将整个仿真模型分解为多个组件，每个组件内部都自带求解器，通过协同仿真的方法，将各个组件的运行和求解分布到不同的线程中，达到并行求解的目的；
5.第二种是在编程语言上进行扩展，引入并行编译指令，在需要并行求解的位置提供并行信息，这些可并行的位置主要针对循环和向量运算；
6.第三种是根据代码间的依赖关系建立任务图，采用自顶向下的调度方法，位于任务图中同一层的各个节点可以同时调度，派发到不同的线程中。
7.在传统仿真并行求解方法中，重点是围绕cpu的多核能力，建立多线程的并行执行环境，而很少考虑gpu的多核能力。这主要归因于gpu的硬件架构特征，只能处理具有单指令多数据集(simd)这种形式的计算，要求所有的求解过程完全是同质的，如同一模型的多个实例同时运行这种情况；而复杂产品中包含的各个子模型通常是各不相同的，如汽车，其内部包含发动机、液压、传动、刹车等多种子模型，子模型内部的计算逻辑各不相同，无法发挥gpu的并行计算能力，使得仿真过程中计算节点的gpu处于空闲状态，无法提供算力支撑，进而导致计算节点内部计算资源的负载不均衡。


技术实现要素：

8.为了解决上述问题，本发明提供一种面向系统仿真模型的gpu并行求解方法。
9.本发明的技术方案总体思想为：构建模型间和模型内两种级别并行相结合的求解方法。首先，根据复杂产品的建模描述xml文件，提取所有子模型的微分代数方程组描述信息，根据这些描述信息，动态生成和编译gpu侧的函数体和数值求解代码，同时生成cpu侧表达式动态求值语法树，建立方程组级别求解任务列表。其次，主线程建立线程池，遍历待求解任务列表，主控算法向线程池派发特定常微分方程组的求解任务。线程池中得到求解任务的线程，利用基于时域分解的多分辨率方法，将特定常微分方程求解计算任务异步推送
到gpu的流处理单元上；gpu的多个流处理器并行求解不同的方程组，每个流处理器内部并行求解特定的常微分方程组，当处理完某个推送任务后，产生任务完成事件，线程池中的某个线程被再次唤醒，对gpu计算结果进行修正；如果修正后的结果不满足精度要求，则再次推送异步求解任务到gpu，通过多次迭代，直至满足求解精度要求，并将结果返回给主控算法，求达到并行快速求解的目的。
10.本发明提出的一种仿真模型任务分解与并行求解的方法，包括：
11.s1、提取常微分方程组描述信息，为gpu动态生成代码，为cpu生成抽象语法树，建立方程组级别求解任务列表；
12.s2、建立线程池，遍历求解任务列表，向线程池派发常微分方程组的求解任务；
13.s3、线程池中获得任务的工作线程，采用基于时域分解的多分辨率方法，将常微分方程组求解计算任务异步推送到gpu的流处理器上；
14.s4、gpu的多个流处理器并行求解不同的方程组，每个流处理器内部并行求解具体的某个常微分方程组，当某个流处理器求解完成后，生成任务完成事件；
15.s5、gpu的任务完成事件被信号处理线程捕获，唤醒线程池中的工作线程；
16.s6、被唤醒的工作线程获取gpu计算结果，进行结果修正，如果修正后结果满足精度要求，返回结果；否则再次向gpu推送异步求解任务，重复步骤s4、s5、s6。
17.进一步地，所述步骤s1中，从复杂产品的建模描述xml文件中以子模型为单位，提取方程组描述信息。
18.进一步地，在步骤s1中，为gpu动态生成代码的步骤包括：
19.解析输入描述信息的名称、方程组表达式、方程组系数部分，生成对应的函数体代码；根据生成的函数体，建立函数指针查找表；将数值求解gpu算法代码与生成代码合并，调用gpu运行时编译指令，生成gpu二进制代码。
20.进一步地，在步骤s1中，为cpu生成抽象语法树的步骤包括：
21.解析输入描述信息的名称、方程组表达式、方程组系数部分，将状态变量及变量加入符号表，再将等号右侧的逻辑运算转换为表达式实例，并将符号表中的符号注册到对应的表达式实例中。然后，再将表达式实例转换为抽象语法树(ast)，并缓存在内存中，以使得对应方程组计算过程可以重复执行。每个方程组对应一个ast语法树，最终得到语法树列表exprlist。
22.进一步地，在步骤s1中，建立方程组级别求解任务列表的步骤包括：
23.解析输入描述信息的状态变量初值、积分时间、全局积分步长、局部积分步长、并行求解并发数量部分，得到求解参数列表solveparamlist。遍历exprlist和solveparamlist，建立求解任务列表tasklist,tasklist中的每个任务由有序对task：＝{expr,solveparam}构成。
24.进一步地，所述步骤s2中，主线程首先建立包含一定数量线程的线程池，然后遍历求解任务列表，将其中的每个任务taski派发到线程池。优选的，线程的数量为cpu硬件总核数的2至3倍。
25.进一步地，所述步骤s3中，某个获取任务的工作线程threadk被唤醒，以taski中的信息为基础，采用基于时域分解的多分辨率方法，调用粗粒度数值求解算法，得到常微分方程odei在的第0次估计值并将其赋值给以此估计值为参数，向gpu推送异步并
行求解任务taski。
26.进一步地，所述步骤s4中，gpu根据其内部执行逻辑，在多个流处理器间并行求解不同的方程组，在每个流处理器内部并行求解特定的常微分方程组；在某一时刻，一个流处理器sp被调度，以为初始值，采用细粒度积分方法执行内部并行求解任务；当任务完成时，产生一个任务完成事件eventi。
27.进一步地，所述步骤s5中，信号处理线程检测到任务完成事件eventi后，唤醒线程池中的工作线程threadj。
28.进一步地，所述步骤s6中，被唤醒的工作线程threadj从gpu取回计算结果以中的第1个元素为初值，调用粗粒度数值求解算法，得到odei在的第1次估计值并对最终能结果做如下修正：
[0029][0030]
如果满足精度要求，则返回
[0031]
否则，以估计值为参数，向gpu推送新的常微分方程异步并行求解任务；
[0032]
跳转到步骤s4，继续迭代，直至满足精度要求并输出结果。
[0033]
本发明提供的基于gpu的异构常微分方程组求解方法，与现有方法相比，具有如下有益效果：
[0034]
(1)本发明提供了两级并行求解能力。将常微分方程组的求解过程映射到gpu不同的流处理器上，使得异构的常微分方程组的求解计算能够在流处理器间并行；同时利用gpu流处理器内部的轻量级线程并发调度能力，在流处理器内以单指令多数据集(simd)的形式完成单个常微分方程组的并行求解。既解决了异构常微分方程组并行求解问题，又可以在单个常微分方程求解上获得较好的加速比。
[0035]
(2)本发明的应用领域更广。面向包含大规模异构常微分方程组的模型，适用于所有涉及常微分方程求解的计算领域，不受限于具体模型，能够加快求解速度，尤其是对那些需要采用隐式数值求解算法，能够有效缩短仿真时间；
[0036]
(3)本发明可提升仿真过程的硬件资源利用率。利用本发明提供的方法，可以根据需要，对模型仿真的求解任务做灵活划分，将原本属于cpu的部分计算任务转移到gpu上，使得计算节点的cpu、gpu负载相对均衡，避免gpu闲置和浪费。
附图说明
[0037]
为了更清楚地说明本公开实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本公开的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些图获得其他的附图。
[0038]
图1为本发明一个实施例的并行求解过程示意图；
[0039]
图2为本发明一个实施例的常微分方程描述示意图；
[0040]
图3为本发明一个实施例的gpu代码动态生成示意图；
[0041]
图4为本发明一个实施例的ode求解任务参数实例示意图；
[0042]
图5为本发明一个实施例的ode求解调度示意图；
[0043]
图6为本发明一个实施例的某个ode计算过程时序图。
具体实施方式
[0044]
为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整的描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动的前提下所获得的所有其他实施例，都属于本发明保护的范围。
[0045]
下面结合附图和具体实施方式对本发明作进一步说明。
[0046]
本发明提供一种仿真模型任务分解与并行求解方法，如图1所示，该方法包括如下步骤：
[0047]
s1、提取仿真模型常微分方程组描述信息，为gpu动态生成代码，为cpu生成抽象语法树，建立方程组级别求解任务列表；
[0048]
s2、建立线程池，遍历任务列表，向线程池派发常微分方程组的求解任务；
[0049]
s3、线程池中获得任务的工作线程，利用基于时域分解的多分辨率方法，将常微分方程求解计算任务异步推送到gpu的流处理单元上；
[0050]
s4、gpu的多个流处理器并行求解不同的方程组，每个流处理器内部并行求解常微分方程组，当某个流处理器求解完成后，生成任务完成事件；
[0051]
s5、gpu的任务完成事件被信号处理线程捕获，唤醒线程池中的工作线程；
[0052]
s6、被唤醒的工作线程获取gpu计算结果，进行结果修正，如果修正后结果满足精度要求，返回结果；否则再次向gpu推送异步求解任务，重复步骤s4、s5、s6。
[0053]
在步骤s1中，提取仿真模型常微分方程描述输入信息，可以从产品的建模描述xml文件中以子模型为单位，提取方程组描述信息，包括名称、方程组表达式、方程组系数、状态变量初值、积分时间、局部积分步长、并行求解并发数量这几个要素。在一个实施实例中，如图2所示，包含了两个不同的常微分方程组，分别是洛伦兹系统和调制波系统。其中，洛伦兹系统包含了输入信息的所有要素，而调制波系统只包含一个方程，可被看作是只包含一个方程的特殊方程组；方程组系数描述部分是可选的，与具体的常微分方程组特征有关，在本实施案例中，调制波系统的描述信息就不包含系数部分。
[0054]
在步骤s1的gpu动态生成代码子步骤中，如图3所示，图2中的两个方程组所对应的函数描述被解析，并动态生成对应的函数体，名称分别为lorenz和modulated_wave。生成函数体后，构建函数指针数组function_list。当数值求解算法需要调用某个常微分方程组函数体时，以参数序号为索引，动态定位到function_list数组中的具体函数体指针，执行具体的数值积分步。
[0055]
除函数体和函数指针数组gpu代码是动态生成的之外，数值求解算法可以是预先编制好的代码，存储在文件中。在编译阶段，动态生成代码和数值求解算法代码被合并，生成gpu二进制代码。
[0056]
在本发明中，gpu的数值求解算法采用高分辨率方法，如4阶龙格-库塔或更高阶算法；算法既可以是显式求解方法，也可以是隐式求解方法，既可以是可变步长，也可以是定长步长。在本实施实例中，求解算法采用定长步长4阶龙格-库塔法。
[0057]
步骤s1的生成cpu抽象语法树子步骤中，根据图2中的描述，解析方程组名称、方程组表达式、方程组系数部分，将状态变量及变量加入符号表，再将等号右侧的逻辑运算转换为表达式实例，并将符号表中的符号注册到对应的表达式实例中。然后，再将表达式实例转换为抽象语法树(ast)，并缓存在内存中，以使得对应方程组计算过程可以重复执行。每个方程组对应一个ast语法树，最终得到包含2个元素的语法树列表exprlist。其中exprlist[0]对应洛伦兹系统，exprlist[1]对应调制波系统。
[0058]
具体的，为了支持面向常微分方程的表达式运行时求值，抽象语法树中支持如下四类运算符：
[0059]
算数运算：加、减、乘、除、求模、求幂；
[0060]
通用运算：绝对值、对数、指数、最大值、最小值、取整、平方根、阈值；
[0061]
三角函数运算：正弦、余弦、正切、余切；
[0062]
布尔运算：逻辑与、逻辑或、逻辑非、逻辑异或、按位与、按位或。
[0063]
这些运算被预先编译为二进制代码，在后序遍历抽象语法树的过程中，可以被直接调用。
[0064]
表达式动态求值也可以通过调用现有的开源表达式求值库实现，如exprtk。
[0065]
在步骤s1的建立方程组级别求解任务列表子步骤中，根据图2中的描述，解析输入描述信息的状态变量初值、积分时间、全局积分步长、局部积分步长、并行求解并发数量部分，得到求解参数列表solveparamlist。遍历exprlist和solveparamlist，建立包含两个元素的求解任务列表tasklist,tasklist中的每个任务由有序对task：＝{expr,solveparam}构成。在本实施实例中，solveparamlist数组中包含两个solveparam类型的元素，每个元素中的具体字段值如图4所示。
[0066]
在步骤s2中，主线程首先建立包含2个数量线程的线程池，然后遍历常微分方程求解任务列表，将其中的每个任务taski派发到线程池。线程的数量通常取当前cpu硬件总核数的2至3倍，所有线程初始时处于等待状态。如果输入的方程组数量较少，则以实际要求解的方程组数量作为线程池的配置。
[0067]
在步骤s3中，如图5所示，某个获取任务的线程threadk,k∈{0,1}被唤醒，以taski,i∈{0,1 中的信息为基础，利用基于时域分解的多分辨率方法，调用粗粒度数值求解算法，得到常微分方程组odei在的第0次估计值并将其赋值给以此估计值为参数，向gpu推送异步并行求解任务。
[0068]
在本实施实例中，如果当前线程处理的是洛伦兹系统的求解任务，则对应task0。首先，根据solveparamlist[0]中字段的具体值，可以得到task0的起始积分时间t0＝0.0，全局积分步长δt0＝0.2，局部积分步长δt＝0.01。选择粗粒度数值求解算法2阶龙格-库塔法。以δt为步长，执行从t0时刻到t0 δt0时刻的累积积分，得到长度为20的一组初值估算向量以作为初始值，向gpu推送异步计算任务。
[0069]
在向gpu推送计算任务时，需要指定并行求解的轻量级线程数量及gpu积分步长。在本实施实例中，solveparamlist[0]中的n＝64。但由于64无法整除20，所以需要进行修正，取n＝20
×
(64
÷
20)＝60，对应的δt＝δt
÷
(64
÷
20)。即gpu的积分步长是cpu积分步长的1/3。实际设置δt、δt及n时，应尽可能同时满足如下联立条件：
[0070][0071]
即n可以同时被(δt
÷
δt)和32整除。其中，第一个条件保持粗粒度和细粒度积分步长之间的整数倍关系，避免结果修正时(详见步骤s6)进行内插或外推运算；第二个条件使得gpu的并发调度单位(wrap)尽可能满负荷运行，减少空转线程。如果无法同时满足这两个条件，应首先满足条件一，并保持n尽可能接近32的整数倍关系。为加快收敛速度和屏蔽gpu硬件调度时间，n取值尽可能大。
[0072]
在步骤s4中，如图5、6所示，gpu的多个流处理器并行求解不同的方程组，每个流处理器内部并行求解特定的常微分方程组。gpu根据其内部执行逻辑，在某一时刻调度一个流处理器sp，采用细粒度算法(4阶龙格库塔法，详见步骤s1)，以向量为初始值，并发数为ni,i∈{0,1}，以积分步长为δti,i∈{0,1}，在δti,i∈{0,1}时间段内，执行特定方程组的并行求解任务；任务完成时，产生一个任务完成事件eventi,i∈{0,1}。
[0073]
在步骤s5中，如图5、6所示，信号处理线程检测到任务完成事件eventi,i∈{0,1}后，唤醒线程池中的线程threadj,j∈{0,1}。
[0074]
在步骤s6中，如图5、6所示，被唤醒线程threadj,j∈{0,1}从gpu取回计算结果向量以中的第1个元素为初值，再次调用粗粒度数值求解算法，得到常微分方程组odei,i∈{0,1}的第1次估计值并对最终能结果做如下修正：
[0075][0076]
判断修正后的结果是否满足精度要求：
[0077][0078]
其中∈为误差限。
[0079]
如果满足精度要求，则返回常微分方程组odei,i∈{0,1}在当前全局积分步δti,i∈{0,1}的结果向量
[0080]
否则，以估计值为参数，向gpu推送新的异步并行求解任务，当前线程threadj,j∈{0,1}进入等待状态。算法运行逻辑转到步骤s4，继续迭代，直至满足精度要求并输出结果。
[0081]
对于同一个常微分方程组odei,i∈{0,1}，算法运行逻辑转到步骤s4后，步骤s4至s6的所有带有上标的向量，其对应上标在运算或求解时递增。即
[0082]
本领域普通技术人员可以理解：以上各实施例仅用以说明本发明的技术方案，而非对其限制；尽管参照前述各实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分或者全部技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本发明权利要求所限定的范围。