CN104156606A - 基于混合计算架构的用于核素迁移控制方程的求解方法 - Google Patents

Abstract

本发明公开了一种基于混合计算架构的用于核素迁移控制方程的求解方法，其包括以下步骤：读取参数阶段，读入计算参数、边界条件和计算网格文件；生成矩阵阶段：生成刚度稀疏矩阵、时间演化稀疏矩阵和载荷向量，并生成对应的稀疏矩阵，之后刚度稀疏矩阵、时间演化稀疏矩阵；循环计算求解阶段，使用有限差分法循环计算，循环过程中采用Jacobi迭代法求解AX＝B的矩阵方程，得到的解向量X即为该时间点的核素浓度场数据。本发明不仅可以存储更高网格体数规模下的数据矩阵，求解速度更快。

Description

基于混合计算架构的用于核素迁移控制方程的求解方法

技术领域

本发明涉及一种方法，特别是涉及一种基于混合计算架构的用于核素迁移控制方程的求解方法。

背景技术

截止2012年，我国已经有15个核反应堆投入运行，另有26个在建，为解决日益临近的核废物处置问题，我国已经明确提出在2020年建成核废物处置地下实验室，2050年建成核废物地下处置库。核废物处置的发展将直接影响和制约我国核事业和国防力量的规划和发展。

如何有效安全的处置核废物，不仅仅影响国家的发展和未来，对于个人和社会也有着直接的影响，最显著的例子是前苏联的切尔诺贝利核电站事故和今年的日本福岛核电站事故，这些鲜明的灾难性事故都说明：核废物一旦无法安全处置，将给人类社会带来灾难性的后果。

目前，人类社会一般采用深埋(地下500-1000米左右)的方式来处置核废物(主要针对中高放废物)，其核心思想是通过各种工程手段、自然屏障来阻止延迟核素从固化体向周围环境(生物圈)迁移。

目前核素迁移研究常见的研究手段包括：实验室模拟研究、天然类比研究、示踪研究、计算机模拟研究。

然而，在传统的求解计算方法中，使用了全矩阵的方式来存储和计算，这种方式在一定网格规模下(例如：小于10万个网格体数)还是可行的，一旦网格规模超过这个规模，一般的计算机已经非常困难。

现有的技术是读取计算参数和网格数据，采用全矩阵记录和存储矩阵数据，采用LU分解求解AX＝B的矩阵方程，采用基于CPU的串行方式计算。这种技术的缺点是：这种方法在一定规模的网格体数据之下(10万以下)，在计算机上求解核素在缓冲/回填材料中扩散数学模型(即控制方程)的时间消耗和内存消耗还可以接受，一旦超出这个范围，由于内存的物理容量限制，会导致在整个求解过程中反复在内存和硬盘之间交换数据，降低求解速度，且由于CPU计算能力和设计目的的限制，浮点运算耗时会让用户难以接受。

发明内容

本发明所要解决的技术问题是提供一种基于混合计算架构的用于核素迁移控制方程的求解方法，其不仅可以存储更高网格体数规模下的数据矩阵，节省了内存开销，求解速度更快。

本发明是通过下述技术方案来解决上述技术问题的：一种基于混合计算架构的用于核素迁移控制方程的求解方法，其特征在于，其包括以下步骤：

S1：读取参数阶段：读入计算参数、边界条件和计算网格文件；

S2：生成矩阵阶段：生成刚度稀疏矩阵、时间演化稀疏矩阵和载荷向量，并生成对应的稀疏矩阵，之后刚度稀疏矩阵、时间演化稀疏矩阵；

S3：循环计算求解阶段：使用有限差分法循环计算，循环过程中采用Jacobi迭代法求解AX＝B的矩阵方程，得到的解向量X即为该时间点的核素浓度场数据。

优选地，所述核素在缓冲/回填材料中迁移控制方程为：

其中C为核素浓度，D_ij为核素在多孔介质中的动力扩散系数：D_ij＝Di_j'+D_ij″，

α_T，α_L分别为横向和纵向弥散系数，为渗流场速度大小，V_j为j方向的渗流速度分量，D_ij″＝D_dTδ_ij，D_d为分子扩散系数，T为多孔介质的弯曲率；

这里D_d＝D₀exp(σΔT)，D₀为参考温度下的分子扩散系数，σ为温度修正系数，ΔT为相对于参考温度的温度差，F为等效衰减系数，Q为源汇项。

优选地，所述S2包括以下步骤：

S21：估算稀疏矩阵大小，若是从1到计算网格中四面体数量的循环，则进入S22，若不是则进入S210；

S22：获取一个四面体的四个顶点坐标，得到单元节点坐标矩阵；

S23：计算雅克比矩阵；

S24：计算雅克比矩阵行列式值，计算雅克比矩阵逆矩阵；

S25：计算控制方程中衰变系数De对于刚度稀疏矩阵的影响，并写入到刚度稀疏矩阵中；

S26：计算控制方程中渗流系数V_i对刚度稀疏矩阵的影响，并写入到刚度稀疏矩阵中；

S27：计算控制方程中扩散系数D_ij对刚度稀疏矩阵的影响，并写入到刚度稀疏矩阵中；

S28：计算控制方程中源汇项Q对刚度稀疏矩阵的影响，并写入到刚度稀疏矩阵中；

S29：计算时间演化稀疏矩阵，若是从1到计算网格中四面体数量的循环，则进入S22，若不是则进入S210；

S210：若是从1到计算网格中三角形数量的循环则进入S211，若不是则进入S214；

S211：获取三角形点坐标，得到三角形节点坐标矩阵；

S212：计算三角形面积，以及通量值；

S213：计算通量对载荷向量的影响，若是从1到计算网格中三角形数量的循环则进入S211，若不是则进入S214；

S214：装载矩阵结束。

优选地，所述生成矩阵阶段步骤S21为估算稀疏矩阵大小，该步骤的实际含义为：在一定的网格体数环境下，估算稀疏矩阵的元素个数，方便申请相应的稀疏矩阵内存，在网格体数为1800左右的时候，即需要生成的全矩阵元素为1800*1800时，其非零元素的个数占全部元素比例不会超过1％，而在日常计算中，网格数一般在2000个以上，所以，此处可以大致取非零元素比例为1％，即当网格体数为N的时候，单一稀疏矩阵元素的总个数为N_non-zero＝N*N*0.01，在此过程中，稀疏矩阵的存储采用行、列、值的格式存储。

优选地，所述求解AX＝B包括以下步骤：

S31：根据衰变系数，计算第一类边界条件点的浓度值；

S32：去除刚度稀疏矩阵、时间演化稀疏矩阵和载荷向量中的第一类边界条件点，得到新的K_new稀疏矩阵、C_new稀疏矩阵、F_new向量；

S33：计算得到temp1和temp2，temp1＝C_new+K_new*inv*0.5，temp2＝K_new*inv*0.5-C_new；

S34：计算B向量，B＝F_new*inv-temp2*B′，B′为上一次的B向量，如果为第一次，则B为0；

S35：A＝temp1；

S36：使用Jacobi迭代法求解AX＝B；

S37：得到解向量X，将解向量X从GPU迁移至CPU，与第一类边界条件点浓度合并，形成完整核素浓度场数据，写入数据文件；

S38：计算矩阵结束。

优选地，所述步骤S32至步骤S36均在GPU中处理。

本发明的积极进步效果在于：本发明采用了稀疏矩阵来替代全矩阵，其最直接的效果就是极大地节省了内存开销。本发明使用到了有限元法(FEM)、有限差分法(FDM)来求解核素迁移方程，在得到AX＝B的矩阵方程之后，采用迭代方法来求解矩阵方程的数值解，得到指定时间点的3维核素浓度场。

附图说明

图1为本发明基于混合计算架构的用于核素迁移控制方程的求解方法的流程图。

图2为本发明中生成稀疏矩阵的流程图。

图3为本发明中单次求解AX＝B的流程图。

具体实施方式

下面结合附图给出本发明较佳实施例，以详细说明本发明的技术方案。

本发明对于不同的网格体数，非零元素占全矩阵全部元素的比例如下表1所示。

表1 非零元素占全矩阵全部元素的比例表

矩阵尺寸	零元素数量	非零元素数量	非零元素比例
				260X260	64580	3020	4.45％
321X321	98956	4085	3.96％
				1836X1836	3346304	24592	0.73％
2130X2130	4508074	28826	0.64％
				11506X11506	132228540	159496	0.12％

由表1可以看出，随着矩阵尺寸的增加，非零元素的比例呈下降趋势，因此，采用稀疏矩阵来存储数据是合理的，可以节省大量的内存。

除此之外，采用CPU/GPU混合计算架构来替代传统的CPU串行计算，将大量的浮点运算从CPU上移至GPU，实现细颗粒度的并行计算，将成倍地提高运算速度，从而降低求解的时间消耗。采用Jacobi迭代法替代LU分解，获得满足精度要求的数值解。当然，求解出来的值实际上为各个时间点的核素浓度场值，将其写入相应的时间数据文件中，供后续可视化模块调用。

CPU(Central Processing Unit)是中央处理器，是计算机的运算核心和控制核心，主要包括运算器(ALU，Arithmetic and Logic Unit)和控制器(CU，Control Unit)两大部件。

GPU(Graphic Processing Unit，图形处理器)是相对于CPU的一个概念，GPU是显示卡的“大脑”，它决定了该显卡的档次和大部分性能，同时也是2D显示卡和3D显示卡的区别依据。2D显示芯片在处理3D图像和特效时主要依赖CPU的处理能力，称为“软加速”。3D显示芯片是将三维图像和特效处理功能集中在显示芯片内，也即所谓的“硬件加速”功能。GPU通用计算技术发展已经引起业界不少的关注，事实也证明在浮点运算、并行计算等部分计算方面，GPU可以提供数十倍乃至于上百倍于CPU的性能。

混合计算是在进行大规模运算时，在把大量并行浮点运算交给GPU的同时，CPU除负责程序流程控制之外，也承担少量运算任务，通过控制计算任务的比例分配，减少任务计算过程中CPU和GPU的资源闲置，从而达到不浪费计算资源的目的。

稀疏矩阵中非零元素的个数远远小于矩阵元素的总数，并且非零元素的分布没有规律，则称该矩阵为稀疏矩阵(sparse matrix)；在本发明中，对于稀疏矩阵采用特殊的格式来保存。矩阵中无论零还是非零元素全部存在，即为全矩阵。

如图1所示，本发明基于混合计算架构的用于核素迁移控制方程的求解方法包括以下步骤：

S1：读取参数阶段：读入计算参数、边界条件和计算网格文件；

S2：生成矩阵阶段：生成刚度稀疏矩阵(K矩阵)、时间演化稀疏矩阵(C矩阵)和载荷向量(F)，并生成对应的稀疏矩阵，之后销毁刚度稀疏矩阵和间演化稀疏矩阵；

S3：循环计算求解阶段：使用有限差分法循环计算，循环过程中采用Jacobi迭代法求解AX＝B(A代表该求解时刻的刚度矩阵，B代表该求解时刻的载荷向量)的矩阵方程，得到的解向量X即为该时间点的核素浓度场数据。这三个阶段的计算(含逻辑计算)的特点和具体的划分处理单元如表2(三个阶段的特点及划分的处理单元表)所示。

表2 三个阶段的特点及划分的处理单元表

从表2可以看出，三个阶段中，前两个阶段以CPU运算为主，主要完成相关的逻辑控制、IO读写和一部分运算；最后一个阶段则以CPU逻辑控制循环，核心的浮点运算均在GPU上完成，以此替代单纯的CPU串行运算，加速整个控制方程的求解。

本发明采用有限元法、有限差分法，以及Jacobi迭代法求解特定核素扩散方程，核素在缓冲/回填材料中迁移控制方程为式(1)：

其中C为核素浓度，D_ij为核素在多孔介质中的动力扩散系数式(2)和式(3)：

D_ij＝Di_j'+D_ij″…………………………………………………………(2)

α_T，α_L分别为横向和纵向弥散系数，为渗流场速度大小，V_j为j方向的渗流速度分量，D_ij″＝D_dTδ_ij，D_d为分子扩散系数，T为多孔介质的弯曲率；

这里D_d＝D₀exp(σΔT)，D₀为参考温度下的分子扩散系数，σ为温度修正系数，ΔT为相对于参考温度的温度差，F为等效衰减系数，Q为源汇项。

其中，步骤S1在求解方程之前，需要针对特定的迁移环境构造相应的计算算例，如该计算算例存在、且各个计算数据正确，则在求解过程开始之前，需要读入计算算例的边界条件和计算网格数据。文件如表3所列，表中mesh子目录包含了网格数据文件的子目录(该子目录下包含了网格点数据文件、网格面数据文件、网格体数据文件、边界条件文件和计算参数文件)。

表3 计算之前需要读入的数据及文件说明表

读入数据后，使用专门的数据结构(数组和向量表)来分别存储，供计算调用。

如图2所示，步骤S2在读入计算网格和边界条件之后，通过对计算网格的四面体进行循环、以及对计算网格面三角形的循环，计算得到刚度稀疏矩阵(K矩阵)、时间演化稀疏矩阵(C矩阵)和载荷向量(F)。步骤S2包括以下步骤：

S21：估算稀疏矩阵大小，若是从1到计算网格中四面体数量的循环(fori＝1totetnum)，则进入S22，若不是则进入S210；

S22：获取一个四面体的四个顶点坐标，得到单元节点坐标矩阵；

S23：计算雅克比矩阵jacobi；

S24：计算雅克比矩阵行列式值det，计算雅克比矩阵逆矩阵invjacobi；

S25：计算控制方程中衰变系数De对于刚度稀疏矩阵(K矩阵)的影响，并写入到刚度稀疏矩阵(K矩阵)中，刚度稀疏矩阵(K矩阵)即刚度矩阵，刚度是表示物质形变能力的一个量，也就是说物体抵抗变形的能力，或者说，是物体产生单位的唯一所需要加载的载荷量。刚度矩阵则是把刚度变到了多维，考虑了在多维的情况下各个维度的相关性。

S26：计算控制方程中渗流系数V_i对刚度稀疏矩阵(K矩阵)的影响，并写入到刚度稀疏矩阵(K矩阵)中；

S27：计算控制方程中扩散系数D_ij对刚度稀疏矩阵(K矩阵)的影响，并写入到刚度稀疏矩阵(K矩阵)中；

S28：计算控制方程中源汇项Q对刚度稀疏矩阵(K矩阵)的影响，并写入到刚度稀疏矩阵(K矩阵)中；

S29：计算时间演化稀疏矩阵(C矩阵)，若是从1到计算网格中四面体数量的循环，则进入S22，若不是则进入S210，时间演化稀疏矩阵(C矩阵)即时间演化矩阵，表示非平衡过程，落实到扩散的话，是浓度梯度导致的传播速度；

S210：若是从1到计算网格中三角形数量的循环则进入S211，若不是则进入S214；

S211：获取三角形点坐标，得到三角形节点坐标矩阵；

S212：计算三角形面积，以及通量值；

S213：计算通量对F向量的影响，若是从1到计算网格中三角形数量的循环则进入S211，若不是则进入S214；

S214：装载矩阵结束。

生成矩阵阶段的第一步S21为“估算稀疏矩阵大小”，该步骤的实际含义为：在一定的网格体数环境下，估算稀疏矩阵的元素个数，方便申请相应的稀疏矩阵内存，从表1的数据可以看出，在网格体数为1800左右的时候，即需要生成的全矩阵元素为1800*1800时，其非零元素的个数占全部元素比例不会超过1％，而在日常计算中，网格数一般在2000个以上，所以，此处可以大致取非零元素比例为1％，即当网格体数为N的时候，单一稀疏矩阵元素的总个数为N_non-zero＝N*N*0.01，在此过程中，稀疏矩阵的存储采用行、列、值的格式存储。图2中的tetNum是指计算网格数据中四面体的数量；triNum是指计算网格中面三角形的数量。“For i＝1 to tetNum”表示的含义为从1到计算网格中四面体数量(tetNum)的循环，“For i＝1 to triNum”则表示从1到计算网格中面三角形数量(triNum)的循环。

如图3所示，步骤S3(求解AX＝B)包括以下步骤：

S31：根据衰变系数，计算第一类边界条件点的浓度值；

S32：去除刚度稀疏矩阵(K矩阵)、时间演化稀疏矩阵(C矩阵)和载荷向量(F)中的第一类边界条件点，得到新的Knew稀疏矩阵、Cnew稀疏矩阵、Fnew向量，F向量即载荷向量，表示外界条件；

S33：计算得到temp1和temp2，temp1＝C_new+K_new*inv*0.5，temp2＝K_new*inv*0.5-C_new；

S34：计算B向量，B＝F_new*inv-temp2*B′，B′为上一次的B向量，如果为第一次，则B为0；

S35：A＝temp1；

S36：使用Jacobi迭代法求解AX＝B；

S37：得到解向量X，将解向量X从GPU迁移至CPU，与第一类边界条件点浓度合并，形成完整核素浓度场数据，写入数据文件；

S38：计算矩阵结束。

步骤S3实质为一个循环，每次循环均采用有限差分法计算新的时间演化稀疏矩阵和刚度稀疏矩阵，以及载荷向量，即每次循环都需要求解AX＝B的矩阵方程，以求出指定时间点的核素浓度场数据。在本步骤中，在GPU上求解得到的解向量即核素浓度场数据需要从GPU迁移到CPU上，然后将其写入到相关时间点下的核素浓度场数据文件中，供可视化模块调用，显示三维核素浓度场。

步骤S32至步骤S36均在GPU中处理，GPU即图形处理器。GPU是相对于CPU的一个概念，GPU是显示卡的“大脑”，它决定了该显卡的档次和大部分性能，同时也是2D显示卡和3D显示卡的区别依据。2D显示芯片在处理3D图像和特效时主要依赖CPU的处理能力，称为“软加速”。3D显示芯片是将三维图像和特效处理功能集中在显示芯片内，也即所谓的“硬件加速”功能。GPU通用计算技术发展已经引起业界不少的关注，事实也证明在浮点运算、并行计算等部分计算方面，GPU可以提供数十倍乃至于上百倍于CPU的性能。将各个计算任务划分至具体的处理单元(CPU或GPU)，以最大化地加速整个计算流程。

根据以上技术方案提出的本发明，与目前多种方法比较具有以下优点：

一、采用稀疏矩阵替代全矩阵，从前面的分析可以看出，所节省的内存存储空间非常大，在11506X11506的网格体数规模下，采用稀疏矩阵存储所消耗的内存仅为采用全矩阵所消耗空间的0.1％左右，如此可以节省大量的内存空间，或者在相同的内存容量下，可以存储更高网格体数规模下的数据矩阵，矩阵中无论零还是非零元素全部存在，即为全矩阵，矩阵中非零元素的个数远远小于矩阵元素的总数，并且非零元素的分布没有规律，则称该矩阵为稀疏矩阵。

二、采用Jacobi迭代法替代LU分解求解AX＝B的矩阵方程，在大规模网格体数背景下，在一定的工程应用精度要求下，该方法求解速度更快。

三、采用CPU/GPU混合计算架构替代传统的CPU串行计算方式来进行求解计算，由于GPU本身的特点(逻辑控制弱、并行浮点运算能力远超CPU)，将其引入到如公式所示的核素扩散方程求解中，能够利用其强大的并行浮点运算能力，极大地加速求解运算过程；又由于CPU的逻辑控制/运算能力非常强大，在整个求解过程中，将偏重于逻辑运算/控制的部分交给CPU完成，则能够利用两者各自的优点来优化、加速整个求解过程。

以上所述的具体实施例，对本发明的解决的技术问题、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本发明的具体实施例而已，并不用于限制本发明，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (6)

1.一种基于混合计算架构的用于核素迁移控制方程的求解方法，其特征在于，其包括以下步骤：

S1：读取参数阶段：读入计算参数、边界条件和计算网格文件；

S2：生成矩阵阶段：生成刚度稀疏矩阵、时间演化稀疏矩阵和载荷向量，并生成对应的稀疏矩阵，之后刚度稀疏矩阵、时间演化稀疏矩阵；

S3：循环计算求解阶段：使用有限差分法循环计算，循环过程中采用Jacobi迭代法求解AX＝B的矩阵方程，得到的解向量X即为该时间点的核素浓度场数据。

2.如权利要求1所述的基于混合计算架构的用于核素迁移控制方程的求解方法，其特征在于，所述核素在缓冲/回填材料中迁移控制方程为：

$A\frac{\∂C}{\∂t}=\frac{\∂}{\∂x_i}\left(D_{\mathrm{ij}}\frac{\∂C}{\∂x_j}\right)-V_i\frac{\∂C}{\∂\mathrm{xi}}-\mathrm{FC}+Q$

其中C为核素浓度，D_ij为核素在多孔介质中的动力扩散系数：D_ij＝Di_j'+D_ij″， ${D_{\mathrm{ij}}}^{\′}={\mathrm{\α}}_T\stackrel{\‾}{V}{\mathrm{\δ}}_{\mathrm{ij}}+({\mathrm{\α}}_L-{\mathrm{\α}}_T)\frac{V_i{V}_j}{\stackrel{\‾}{V}},$

α_T，α_L分别为横向和纵向弥散系数，为渗流场速度大小，V_j为j方向的渗流速度分量，D_ij″＝D_dTδ_ij，D_d为分子扩散系数，T为多孔介质的弯曲率；

这里D_d＝D₀exp(σΔT)，D₀为参考温度下的分子扩散系数，σ为温度修正系数，ΔT为相对于参考温度的温度差，F为等效衰减系数，Q为源汇项。

3.如权利要求1所述的基于混合计算架构的用于核素迁移控制方程的求解方法，其特征在于，所述S2包括以下步骤：

S21：估算稀疏矩阵大小，若是从1到计算网格中四面体数量的循环，则进入S22，若不是则进入S210；

S22：获取一个四面体的四个顶点坐标，得到单元节点坐标矩阵；

S23：计算雅克比矩阵；

S24：计算雅克比矩阵行列式值，计算雅克比矩阵逆矩阵；

S25：计算控制方程中衰变系数De对于刚度稀疏矩阵的影响，并写入到刚度稀疏矩阵中；

S26：计算控制方程中渗流系数V_i对刚度稀疏矩阵的影响，并写入到刚度稀疏矩阵中；

S27：计算控制方程中扩散系数D_ij对刚度稀疏矩阵的影响，并写入到刚度稀疏矩阵中；

S28：计算控制方程中源汇项Q对刚度稀疏矩阵的影响，并写入到刚度稀疏矩阵中；

S29：计算时间演化稀疏矩阵，若是从1到计算网格中四面体数量的循环，则进入S22，若不是则进入S210；

S210：若是从1到计算网格中三角形数量的循环则进入S211，若不是则进入S214；

S211：获取三角形点坐标，得到三角形节点坐标矩阵；

S212：计算三角形面积，以及通量值；

S213：计算通量对载荷向量的影响，若是从1到计算网格中三角形数量的循环则进入S211，若不是则进入S214；

S214：装载矩阵结束。

4.如权利要求1所述的基于混合计算架构的用于核素迁移控制方程的求解方法，其特征在于，所述生成矩阵阶段步骤S21为估算稀疏矩阵大小，该步骤的实际含义为：在一定的网格体数环境下，估算稀疏矩阵的元素个数，方便申请相应的稀疏矩阵内存，在网格体数为1800左右的时候，即需要生成的全矩阵元素为1800*1800时，其非零元素的个数占全部元素比例不会超过1％，而在日常计算中，网格数一般在2000个以上，所以，此处可以大致取非零元素比例为1％，即当网格体数为N的时候，单一稀疏矩阵元素的总个数为N_non-zero＝N*N*0.01，在此过程中，稀疏矩阵的存储采用行、列、值的格式存储。

5.如权利要求1所述的基于混合计算架构的用于核素迁移控制方程的求解方法，其特征在于，所述S3包括以下步骤：

S31：根据衰变系数，计算第一类边界条件点的浓度值；

S32：去除刚度稀疏矩阵、时间演化稀疏矩阵和载荷向量中的第一类边界条件点，得到新的K_new稀疏矩阵、C_new稀疏矩阵、F_new向量；

S33：计算得到temp1和temp2，temp1＝C_new+K_new*inv*0.5，temp2＝K_new*inv*0.5-C_new；

S34：计算B向量，B＝F_new*inv-temp2*B′，B′为上一次的B向量，如果为第一次，则B为0；

S35：A＝temp1；

S36：使用Jacobi迭代法求解AX＝B；

S37：得到解向量X，将解向量X从GPU迁移至CPU，与第一类边界条件点浓度合并，形成完整核素浓度场数据，写入数据文件；

S38：计算矩阵结束。

6.如权利要求5所述的基于混合计算架构的用于核素迁移控制方程的求解方法，其特征在于，所述步骤S32至步骤S36均在GPU中处理。