一种三维中子输运特征线法的负载均衡并行方法与流程

本发明属于反应堆中子输运计算技术，具体涉及一种三维中子输运特征线法的负载均衡并行方法。

背景技术：

反应堆堆芯由大量可裂变的核燃料构成，其中的原子核与穿行在堆芯中的高速中子发生核裂变反应。裂变反应产生的具有一定的速度的中子继续与其余原子核反应，使这一链式反应过程得以延续。中子数量增长太快时，链式反应过于剧烈，堆芯中的热量无法及时由冷却剂带走，易引起事故；中子数量持续减少时，链式反应无法继续，会造成停堆。堆芯内燃料的类型、数量、排布都可能影响反应的进行，在反应堆设计时应该严格评估大量方案。稳态中子输运(neutrontransport)计算是模拟核裂变反应的重要手段之一，快速、精确的计算不仅可以加速反应堆设计过程，而能在反应堆运行时评估其运行状况。它假定较短时间内同位素不发生损耗，通过求解中子输运方程来模拟短时间内堆芯中发生的核裂变反应。实际工程设计中，稳态中子输运会根据设计方案的调整而反复运行，其运行时间、结果精确性都对工程进度有重大影响。工程中常用粗网格、少能群或均匀化方法完成输运计算，快速但不够精确，其中的代表性方法如节块法。细网格、多群、非均匀化方法有更高精度，但其开销对于小型集群来说是难以承担的，其中的代表性方法如蒙特卡洛法、特征线法。特征线法(themethodofcharacteristics,moc)是确定论方法的一种，相比于蒙特卡洛方法，它也可以适应不规则几何；相比于其他确定论方法，它可以消除其他方法中普遍存在的射线效应，且在合适的射线追踪策略下还可以消除数值弥散，让计算得到的中子通量和功率分布更加精确。

特征线法的主要过程是：读取几何和核数据，射线追踪(即生成大量轨迹并对轨迹分段)，迭代求解特征方程。特征方程即中子输运方程的特征形式如下：

其中，ψ为中子角通量，s为沿特征线或沿轨迹的距离，ω为特征线的方向，e为中子能量，∑^t为总截面，q为总的角中子源。其中，中子源又分为散射源和裂变源，在各向同性的条件下，它是标通量φ的函数。特征方程的求解是沿着轨迹进行的。当轨迹被几何网格截成不同长度的线段后，每一段的长度、所在网格的材料信息就被用于计算线段端点的中子角通量。线段起点的中子角通量用于计算终点的中子角通量，终点的中子角通量又作为下一线段起点的中子角通量，直到计算了整条轨迹的角通量贡献。轨迹到达几何边界后会根据边界条件的不同连接到下一条轨迹，此时的原轨迹称为出射轨迹，下一条轨迹称为入射轨迹。出射轨迹的中子角通量要传递给入射轨迹，作为入射轨迹的初始中子角通量。

现有的技术、软件大都利用了特征方程求解的特点，对求解任务进行划分，加速计算。模块化的射线追踪方法(filippone,w.,woolf,s.,lavigne,r.:particletransportcalculationswiththemethodofstreamingrays.nuclearscienceandengineering77(2),119-136(1981).)利用了规则形状的几何横截面为四边形、正六边形，通过平移对称性将几何划分为数个模块，每个模块上的轨迹分布都是相同的。这种射线追踪方法极大地降低了算法复杂性，减少了计算和内存的开销，适合大规模并行计算。美国轻水堆特征线软件mpact，快堆特征线法软件proteus，特征线法的开源软件openmoc，以及西安交通大学的堆芯计算软件necp-x都成功应用模块化射线追踪，通过分解几何区域、轨迹、能群到不同mpi(messagepassinginterface)进程，利用超级计算机加速求解特征方程。

当几何区域被分解到不同进程时，被截断轨迹需要传递中子角通量，导致相邻进程的边界上产生通信。随着区域划分得更多，被截断的轨迹也更多。当轨迹被分解到不同的进程时，一个进程中的轨迹到达边界后可能要连接到另一进程中的轨迹，它们之间也有通信。分配给进程的轨迹越多，需要通信的中子角通量也越多。上述划分方式中，中子角通量都只能使用点到点通信来发送和接收，这会使得大规模并行计算的效率变低，让整个稳态中子输运计算时间变长。

技术实现要素：

本发明的目的是提供一种能够减少轨迹间通信的三维中子输运特征线法的负载均衡并行方法，用以解决大规模并行特征线法计算过程中通信量过大的问题。

本发明的技术方案如下：一种三维中子输运特征线法的负载均衡并行方法，包括如下步骤：

(1)读取输入参数，包括几何的大小、方位角的数量、轨迹的间距、进程数量；

(2)根据几何的对称性，计算出需要考虑的角度数量i，并计算每个角度下周期性轨迹或反射性轨迹的数量k1,k2,…,ki，以及周期性轨迹或反射性轨迹长度l1,l2,…,li；

(3)根据周期性或反射性轨迹的长度、进程数量以及轨迹的数量，计算出每个进程的负载；

(4)各进程分别进行射线追踪，遍历所有轨迹并结合几何网格信息对轨迹进行分段；

(5)根据几何网格信息、材料数据、射线追踪完成后的线段信息，迭代求解特征方程。

进一步，如上所述的三维中子输运特征线法的负载均衡并行方法，步骤(1)中所述几何的大小包括几何的边长和高度，方位角的数量为4的倍数，轨迹的间距为正实数，进程数量为正整数。

进一步，如上所述的三维中子输运特征线法的负载均衡并行方法，步骤(2)中所述的周期性轨迹是指轨迹从周期性边界出射后，从另一边界入射，如此反复，这些相连的轨迹合在一起被称为周期性轨迹；反射性轨迹是指轨迹从反射性边界出射后，从另一边界入射，如此反复，这些相连的轨迹合在一起被称为反射性轨迹；反射性轨迹的数量为周期性轨迹数量的一半，反射性轨迹的长度为周期性轨迹长度的两倍。

进一步，如上所述的三维中子输运特征线法的负载均衡并行方法，步骤(2)中，根据周期性轨迹或反射性轨迹的数量、长度，以及进程数量，对轨迹数量和间距进行修正，若进程数量大于周期性轨迹或反射性轨迹数量的一半，为每个角度增加周期性轨迹或反射性轨迹的数量，直到进程数不超过轨迹数量的一半，增加轨迹后重新计算轨迹之间的间距。

进一步，如上所述的三维中子输运特征线法的负载均衡并行方法，步骤(3)中，计算每个进程的负载包括如下步骤：

(3-1)初始化进程编号r＝0，初始化当前进程所负责的每个角度下的周期性轨迹或反射性轨迹的起始编号(sr,1,sr,2,…,sr,i)＝(0,0,…,0)，初始化当前进程所负责的周期性轨迹或反射性轨迹数量(wr,1,wr,2,…,wr,i)＝(0,0,…,0)；

(3-2)计算平均负载平均负载在分配任务的过程中不发生变化，p为进程数量；

(3-3)更新起始编号sr,i＝sr,i+wr,i,i＝1,2,…,i；

(3-4)构造整数线性规划问题，求解目标函数wr为编号为r的进程的负载；

(3-5)保存(sr,1,sr,2,…,sr,i)和(wr,1,wr,2,…,wr,i)到数据库，同时将这两个数组发送给进程r，让进程r可以立即进入射线追踪阶段；

(3-6)从总的周期性轨迹或反射性轨迹数量中减去已分配的数量ki＝ki-wr,i,i＝1,2,…,i；

(3-7)如果不是最后一个进程，回到步骤(3-3)；否则，把剩余的轨迹都分配给最后一个进程，完成负载的计算。

进一步，如上所述的三维中子输运特征线法的负载均衡并行方法，步骤(4)中对轨迹进行分段的方法是：从轨迹的起点开始，逐个求出轨迹与几何网格的所有交点，求出每两个相邻交点间的距离作为线段长度，实现轨迹的分段。

进一步，如上所述的三维中子输运特征线法的负载均衡并行方法，步骤(5)中迭代求解特征方程的方法为：

(5-1)在源迭代开始后，每个进程遍历自己负责的所有轨迹，使用每条线段和网格上的截面数据计算源项，优选的，每个进程使用多线程(openmp接口)遍历它的所有轨迹，加速计算；

(5-2)计算每条线段对标通量的贡献量；

(5-3)进行进程间通信，同步标通量；

(5-4)若源迭代已经收敛，则结束；否则，进入下一次迭代。

本发明的有益效果如下：

1)该方法将周期性轨迹和反射性轨迹完整地分配给线程，可以消除进程间的点到点通信。进程之间的中子角通量不再需要通信，因为在边界上相连的轨迹都属于同一进程，不会分散到不同进程中。

2)该方法容易保证负载平衡。周期性轨迹和反射性轨迹的数量和长度可以计算得到。根据这些信息可以建立起关于并行任务分配的最优化模型，使用最优化方法(如线性规划、动态规划)可以得到相应的并行任务分配策略。

3)该方法可以高效地分配并行任务。任务的分配主要由几何形状、角度数量、轨迹间距、进程数量这几个参数决定，在反应堆设计中，几何形状往往不会频繁变化，角度数量和轨迹间距都有经验值，变化范围不大，进程数量则与机器相关。通过调整这些参数，可以预先生成并行任务分配数据库，程序运行时只需从数据库中读取相应的任务分配策略即可，这能极大节省预处理过程的时间。

4)该方法适应于多种几何建模方法和迭代求解策略。使用组合几何、有限元等方法对反应堆进行网格划分时都可以使用本方法划分并行任务。使用内-外迭代法或单层外迭代法都可以使用本方法划分并行任务，因为本方法不影响在轨迹上所做的计算，只影响轨迹的分配方式。

附图说明

图1为本发明具体实施例中任务分配和负载平衡方法的流程图。

图2为本发明具体实施例中四边形几何上的循环轨迹分布示意图。整个几何的长宽分别为a,b，沿坐标轴方向的方向向量为ex,ey，轨迹的垂直间距为δ⊥，轨迹沿坐标轴的间距分别为δx,δy。若周期性轨迹数量为k>1，轨迹间的间距、周期性轨迹的垂直间距、轨迹数量三者满足图中的周期性轨迹由5条轨迹节组成，满足n＝3,m＝2。每条轨迹到达几何边界后都可以按周期性边界条件折回另一边界，即，这些轨迹总是在边界相连。若边界为反射性边界条件，反射后的轨迹的角度与反射前的角度是互补的，反射前后的轨迹是轴对称的。

图3a-图3d为本发明具体实施例中二维横截面上的周期性轨迹示意图。图中注明了一个角度下的轨迹的编号范围是[0,nx+ny]，其中，锐角和钝角的编号方向不同。根据进程r负责的轨迹的起始编号可以得到图中的轨迹编号t＝sr,i，从而可以计算出轨迹起点。图中端点为实心圆点的线表示出射轨迹，实心圆点为出射点；端点为箭头的线表示入射轨迹，入射点为上一条轨迹的出射点。每当轨迹从边界出射，都按照周期性边界条件从另一边界入射，得到的相应入射轨迹的编号可以由图中的公式计算。

图4a、图4b为本发明具体实施例中二维横截面上的反射性轨迹示意图。图中注明了一个角度下的轨迹的编号范围是[0,nx+ny]，其中，锐角和钝角的编号方向不同。根据进程r负责的轨迹的起始编号可以得到图中的轨迹编号t＝sr,i，从而可以计算出轨迹起点。图中端点为实心圆点的线表示出射轨迹，实心圆点为出射点；端点为箭头的线表示入射轨迹，入射点为上一条轨迹的出射点。每当轨迹从边界出射，都按照反射性边界条件从当前边界入射，得到的相应入射轨迹的编号可以由图中的公式计算，比如(φ^c,t-ny)表示到达边界后，角度应该变为补角，编号应该减去ny。

图5为本发明具体实施例中垂直面上的周期性或反射性轨迹示意图。图中注明了垂直面上轨迹经过周期边界或反射边界时编号的对应关系。垂直面的底边是一条二维横截面上的周期性轨迹或反射性轨迹，侧边是纵轴。两种情况下得到的垂直面都是矩形，因此在垂直面内穿过两种边界的轨迹有相同的模式。图中轨迹在垂直面内的编号为l，经过周期性边界或反射性边界后得到的轨迹编号为l+nl，其中nl为垂直面底边上发出的轨迹数量。经过周期性边界的轨迹不改变角度，经过反射性边界的轨迹的角度应变为原角度的补角。

具体实施方式

下面结合附图和实施例对本发明进行详细的描述。

首先对本发明中涉及到的技术术语进行如下说明：

方位角：方位平面内的角度，取值为[0,2π]。

极角：极平面内的角度，取值为[0,π]。

周期性边界条件：由平移变换带来的边界条件。轨迹从这种边界出射后，会从某个对应的边界入射，且轨迹的方向不变。

反射性边界条件：使轨迹发生全反射的边界条件。轨迹从这种边界出射后，入射点就是出射点，且角度会变成原有角度的补角。

周期性轨迹：轨迹从周期性边界出射后，从另一边界入射，如此反复。这些相连的轨迹合在一起被称为周期性轨迹。

反射性轨迹：轨迹从反射性边界出射后，从另一边界入射，如此反复。这些相连的轨迹合在一起被称为反射性轨迹。

源迭代：特征方程的迭代求解过程中引入的最外层迭代。

截面数据：核截面，衡量原子核与中子发生不同核反应可能性大小的参数，与核素类型、中子能量大小相关。

源项：特征方程中的角中子源和标量源项，包括裂变源和散射源。

反应堆通常符合轴向挤出几何的定义，即每个横截面都有完全类似的几何形状，因此可以在一个横截面内分配所有方位角下的轨迹，再以每条周期性轨迹为轴建立垂直于横截面的垂直面，在垂直面中应用四边形几何的射线追踪方法。轻水堆和快堆的输运计算中，主要涉及反射性边界条件和周期性边界条件，而真空边界条件可以很容易地处理。四边形几何和六边形几何都满足平移不变性，应用于它们的平移变换都构成有两个生成元e1,e2的置换群其中每个元素都可以表示为一组向量。在四边形几何中，生成元可以表示为e1＝aex,e2＝bey，其中a,b为几何长宽，ex,ey为单位向量；在六边形几何中，生成元可以表示为其中a为六边形边长。满足这些条件的轨迹都属于循环的轨迹分布，其中的轨迹总是首尾相接，不需要在边界上对通量进行插值。这种首尾相撞构成回路的轨迹称为周期性轨迹。每条周期性轨迹的角度、长度以及轨迹之间的间距都可以计算得到。对于四边形几何，角度为长度为轨迹间的间距为对于六边形几何，角度为长度为轨迹间的间距为四边形的情况如图2所示。减少轨迹间距可以增加周期性轨迹的数量，若数量为k，则相应的间距

以周期性轨迹为例，由上述计算可知，同一角度下的所有周期性轨迹都具有相同长度。假设(0,2π)内取个角度(角度必须是4的倍数)生成轨迹，实际只考虑(0,π)内的个角度，而让每条轨迹具有正向和反向两个方向。则每个角度的周期性轨迹数量为k1,k2,…,ki，长度的取值分别为l1,l2,…,li，且满足ki＝ki+1-i,li＝li+1-i。总的周期性轨迹数量为k＝(k1+k2+…+ki)。由于计算量是随着轨迹长度而增加的，用w＝(k1l1+k2l2+…kili)代表总计算量。在具有p个mpi进程时，理想情况下每个进程的平均负载为为每个进程分配任务时，第r个进程分得的任务可以由一个向量(wr,1,wr,2,…,wr,i)，其中每个值wt,i表示进程t应该分得的第i个角度下的周期性轨迹的数量。进程的负载可以表示为wr＝(wr,1l1,wr,2l2,…,wr,ili)。负载均衡的情况下，各进程的负载相差不大。单个进程分得的轨迹数量可以由如下整数线性规划问题求得：

约束条件：

反射性轨迹的原理是一样的，对于算法来说没有任何区别。

图1说明了本发明中并行特征线法任务分配和负载平衡方法的流程，包括如下步骤：

1.读取输入参数

首先，读取必要的数据，包括几何的大小(四边形几何为边长a,b，正六边形几何为边长a，几何的高度均为c)、方位角的数量轨迹的间距进程数量p。其中，为4的倍数，a,b,均为正实数，p为正整数。各项数据的意义如图2。

2.初始化各项参数

根据几何的对称性可以计算出需要考虑的角度数量为不再为其余的i个角度生成新的轨迹，而是复用已有的轨迹。根据以上参数，结合上取整运算，可以计算出每个角度下分布在各个轴上的轨迹的数量i＝1,2,…,i。满足此条件的轨迹能够构成周期性轨迹和反射性轨迹。

根据轨迹的数量可计算出周期性轨迹的数量ki＝gcd(nx,ny)，即两者的最大公约数(gcd)，得到所有角度下周期性或反射性轨迹的数量k1,k2,…,ki，其中，反射性轨迹的数量为周期性轨迹数量的一半。同时可以计算出周期性或反射性轨迹的长度l1,l2,…,li，其中，反射性轨迹的长度为周期性轨迹长度的两倍。

根据周期性或反射性轨迹的数量、长度，以及进程数量，可以修正轨迹数量和间距。若进程数大于总的轨迹数量的一半为每个角度增加周期性轨迹或反射性轨迹的数量，直到进程数不超过轨迹数量的一半。每次增加数量时，第1个角度增加1条周期性轨迹或反射性轨迹，第i个角度增加条周期性轨迹。增加轨迹后重新计算轨迹之间的间距。

3.计算每个进程的负载

根据周期性或反射性轨迹的长度、进程数量以及修正后的轨迹数量，计算出每个进程的负载。具体步骤如下：

(3-1)初始化进程编号r＝0，初始化当前进程所负责的每个角度下的周期性轨迹的起始编号(sr,1,sr,2,…,sr,i)＝(0,0,…,0)，初始化当前进程所负责的周期性轨迹数量(wr,1,wr,2,…,wr,i)＝(0,0,…,0)；

(3-2)计算平均负载平均负载在分配任务的过程中不发生变化；

(3-3)更新起始编号sr,i＝sr,i+wr,i,i＝1,2,…,i；

(3-4)构造整数线性规划问题，求解目标函数wr为编号为r的进程的负载；

(3-5)保存(sr,1,sr,2,…,sr,i)和(wr,1,wr,2,…,wr,i)到数据库，同时使用非阻塞发送将这两个数组发送给进程r，让进程r可以立即进入射线追踪阶段；

(3-6)从总的周期轨迹数量中减去已分配的数量ki＝ki-wr,i,i＝1,2,…,i；

(3-7)如果不是最后一个进程，回到步骤(3-3)；否则，把剩余的轨迹都分配给最后一个进程，完成负载的计算。

4.各进程分别进行射线追踪

进行射线追踪之前先读取其他必要数据。包括极角的数量j，每个极角取值在(0,π)内。

每个进程根据取得的(sr,1,sr,2,…,sr,i)和(wr,1,wr,2,…,wr,i)两个数组(可从文件读取或从主进程获取)，分别进行射线追踪，即遍历所有轨迹并结合几何网格信息把轨迹截成线段。进程首先遍历二维横截面内的所有方位角i＝1,2,…,i，对于第i个角度，遍历wr,i条周期性轨迹或反射性轨迹，编号为sr,i,sr,i+1,…,sr,i+wr,i-1。对于每条周期性轨迹，遍历其中首尾相连的所有轨迹，轨迹编号的对应关系按照图3中的公式计算。对于每条反射性轨迹，遍历其中首尾相连的所有轨迹，轨迹编号的对应关系按照图4中的公式计算。进程随后遍历垂直面内的所有极角θj,i＝1,2,…,j，垂直面内轨迹的数量计算公式与二维横截面内相同，此时的几何长为周期性或反射性轨迹长度li，几何宽为c。垂直面内轨迹编号的对应关系按图5计算。

5.迭代求解特征方程

根据几何网格的信息、材料数据、射线追踪完成后的线段信息，迭代求解特征方程。若进程负责的周期性或反射性轨迹数量超过1，则使用多线程(openmp接口)，将轨迹分配给不同线程，加速计算。在源迭代开始后，每次迭代中，每个进程遍历自己负责的所有轨迹，使用每条线段和网格上的截面数据计算源项。进程遍历每条轨迹的每条线段，并按特征方程计算每条线段对标通量的贡献量。完成进行进程间通信，同步标通量。若源迭代已经收敛，则结束；否则，进入下一次迭代。

显然，本领域的技术人员可以对本发明进行各种改动和变型而不脱离本发明的精神和范围。这样，倘若对本发明的这些修改和变型属于本发明权利要求及其同等技术的范围之内，则本发明也意图包含这些改动和变型在内。