基于gpu的多粒子输运仿真方法

基于gpu的多粒子输运仿真方法
【专利摘要】公开了一种基于GPU的粒子输运仿真方法，包括：初始化；根据粒子的数目划分GPU的线程数目；在GPU并行计算粒子的下一个碰撞点，判断是否有粒子的权重分裂或者死亡，给出碰撞点将要发生的事件；将计算结果返回给CPU，CPU根据事件数目及类型重新划分GPU的线程数目；在GPU中并行处理这些事件，给出事件发生后粒子新的能量，方位角等信息以及产生的二次粒子的信息，并返回给CPU；CPU汇总这些信息并返回到程序开头，依此类推循环，最后统计输出粒子输运的仿真结果。根据本发明实施例，能够实现高效的粒子输运实验仿真，并在降低能耗的同时提高运算效率。
【专利说明】基于GPU的多粒子输运仿真方法
【技术领域】
[0001]本发明示例性实施例总体上涉及多粒子输运仿真【技术领域】，具体地，涉及一种可应用于中子、光子和电子混合输运研究的基于GPU的蒙特卡罗输运方法。
【背景技术】
[0002]粒子输运模拟在科学与工程计算中具有重要意义，其求解算法发展至今，可以分为:确定论方法和蒙特卡洛(Monte Carlo)方法。由于蒙特卡洛方法具有收敛速度与问题的维数无关、受问题的条件限制影响小且程序结构简单等特点，在处理粒子输运相关时，蒙特卡洛方法对于处理截面与能量、散射各向异性、介质非均匀、复杂几何、事件相关等问题时具有得天独厚的优势。但其最大的缺点在于收敛速度慢，误差的概率性质使得需要大量的抽样才能保证准确度。目前粒子输运模拟计算的实现以CPU为主，由于CPU计算能力不足直接影响到其计算效率，因此在可以接受的机时内只能计算很小的空间尺寸和时间尺寸；或者为了提高计算量可以建设大规模甚至超大规模的CPU计算机集群，但这样又会出现耗电量过大、使用和维护成本极高等问题。

【发明内容】

[0003]本发明示例性实施例提出了一种基于GPU的多粒子输运仿真方法，通过在GPU上并行执行蒙特卡洛(MC)算法来进行多粒子输运仿真，在降低能耗的同时提高运算效率，该方法包括以下步骤:
[0004]a、使用蒙特卡洛算法对粒子输运建模，建立包括多种几何体的空间，向CPU和GPU分别分配内存空间，在CPU中进行数据初始化，将初始化后的数据从CPU的内存空间拷贝到GPU的内存空间中；
[0005]b、在CPU中从源产生多个粒子，所述多个粒子包括光子、中子和电子，每个粒子具有粒子信息，包括粒子的种类、位置、速度、权重及所在几何体的编号，将粒子信息连同粒子数目储存在相应数组中；
[0006]C、判断光子和中子数目之和是否大于设定值，如果数目之和大于设定值，则根据该数目之和确定要使用的线程数目，并将每个中子或光子分配给GPU中的一个线程进行处理，计算前进距离以及判断中子或光子是否死亡；如果数目之和不大于设定值，则由CPU处理每个中子或光子，计算前进距离以及判断中子或光子是否死亡；
[0007]d、由CPU依据中子、光子和电子的碰撞事件数目重新划分线程数目，并且GPU中每个线程处理一个碰撞事件；
[0008]e、判断粒子是否死亡，如果继续存活，将粒子的信息更新；
[0009]f、判断是否产生二次粒子，如果产生，将二次粒子的信息存储；
[0010]g、由CPU重新统计所有存活的粒子信息，统计需要处理的粒子的数目；
[0011]h、判断运行时间是否大于设定值，如果大于设定值，则输出统计量并释放CPU和GPU的内存空间，所述方法结束；[0012]1、如果运行时间不大于设定值，则判断要处理的粒子的数目是否小于设定的最小值；如果小于最小值，则返回步骤b，否则返回步骤c ；
[0013]其中GPU包括多个流处理器，所述多个流处理器通过多个线程对粒子输运仿真进行并行处理，并且所述多个流处理器对粒子的输运仿真和碰撞事件的处理是同步进行的。
[0014]根据示例实施例,在步骤c还包括:
[0015]对于处理的中子或光子，获得时间截断距离、到达边界距离、平均自由程和碰撞距离；
[0016]计算以上距离之中的最小距离；
[0017]如果最小距离是时间截断距离，则确定该中子或光子死亡；
[0018]如果最小距离是达到边界距离，则更新该中子或光子的所在几何体的编号；
[0019]如果最小距离是平均自由程，则用其余三个距离减掉平均自由程，并返回到步骤c的开始；
[0020]如果最小距离是碰撞距离，则判断碰撞是中子碰撞还是光子碰撞。
[0021]根据示例实施例，在步骤e中如下处理碰撞事件:
[0022]如果碰撞事件为电子碰撞事件，则计算电子在上一个最小距离中的能量损失，将产生的二次粒子的数目和状态进行记录，并返回给CPU ；
[0023]如果碰撞事件为光子碰撞事件,则抽样碰撞核素，处理光子反应，将产生的二次粒子的数目和状态进行记录，并返回给CPU ；
[0024]如果碰撞事件为中子碰撞事件，则抽样碰撞核素,处理中子反应,将产生的二次粒子的数目和状态进行记录，并返回给CPU。
[0025]本发明实现了基于GPU的多粒子输运蒙特卡洛仿真系统，充分利用GPU的超强浮点运算能力、高带宽及多轻量计算核心的特点，及GPU内众多流处理器，将蒙特卡洛算法合理引入粒子输运仿真中，使得蒙特卡洛粒子输运仿真适合于GPU硬件架构。在本发明中，GPU和CPU充分结合，发挥各自优势，极大地提高了运算效率。在可以接受的能耗和时间条件下，取得了非常好的计算效果，达到了提高运算效率的效果。
【专利附图】

【附图说明】
[0026]下面结合附图对本发明的【具体实施方式】作进一步详细的说明，其中:
[0027]图1示出了根据本发明示例实施例的基于GPU的多粒子输运仿真方法的流程图。
【具体实施方式】
[0028]下面结合附图对本发明的示例实施例进行详述。以下描述包括各种具体细节以辅助理解，但这些具体细节应仅被示为示例性的。因此，本领域普通技术人员将认识到，可以在不脱离本公开范围和精神的情况下对这里描述的各个实施例进行各种改变和修改。此夕卜，为了清楚和简明起见，省略了公知功能和结构的描述。
[0029]以下描述和权利要求中使用的术语和词语不限于其字面含义，而是仅由发明人用于实现本发明的清楚一致的理解。因此，本领域技术人员应当清楚，对本发明各个示例实施例的以下描述仅被提供用于说明目的，而不意在限制由所附权利要求及其等同物限定的本发明。[0030]本发明实施例提供了。以下结合附图详细说明本发明实施例。
[0031]图1示出了根据本发明示例实施例的基于GPU的多粒子输运仿真方法的流程图。如图1所示，该多粒子输运仿真方法包括如下步骤。
[0032]步骤1、使用蒙特卡洛方法对粒子输运建模，建立由各种几何体(后面简称为体，体由面围成,每个体和面均有对应的编号)组成的空间，并在CPU中进行初始化数据并拷入GPU内存。这里，粒子可以包括中子、光子和电子。
[0033]步骤2、在CPU中从源产生多个粒子，其中每个粒子均包含如下信息:该粒子的存活标志live,该粒子编号idf,该粒子的重复数目npa,该粒子所在体的编号icl,种类ipt,截面库的编号iex,绝对坐标XXX, yyy, zzz,方位角uuu, vvv, www,能量erg,权重wgt等。
[0034]在步骤2中，如果粒子是电子，则粒子信息还包括电子当前能量erg在能量网格中对应的编号ngp,小步步数ns (将电子走的路径分为若干大步,大步内的能量损失率不变，每大步又细分为ns个小步)。判断标识El (如El为I则在处理电子时，需要计算电子在当前体的能量损失率qs及小步步数ns。计算完毕后，就恢复El为0，避免下次跳转到该步骤时重复计算qs及ns);如果粒子是中子或者光子，则粒子信息还包括碰撞标识NPl (I代表中子碰撞，2代表光子碰撞)。
[0035]步骤3、如果中子和光子数目超过设定值，则将每个中子和/或光子分配给GPU中一个线程进行处理(电子不参与该过程)，若数目小于该设定值，则由CPU处理每个中子和/或光子，包括计算前进距离，判断粒子是否死亡。
[0036]在步骤3中如果处理的是中子或者光子，则还包含以下子步骤:
[0037]3.1:计算粒子到体边界的的距离(即沿当前运动方向穿出当前所在体需要走的距离dls)、时间截断距离(即粒子在剩余寿命内按照当前的速度所能通过的距离dtc (每种粒子均设置了所能存活的最长寿命))、平均自由程dw (权重窗存在情况下为平均自由程，否则为无穷大，其中权重窗为空间-能量或空间-时间构成的相空间)、及到下次碰撞点的距离pmf ；
[0038]3.2:计算步骤3.1中四种距离中哪种距离具有最小值，并更新粒子坐标XXX，yyy，ZZZ和已存活时间tme ；
[0039]3.3:如果步骤3.1中的最小距离为时间截断距离dtc，则判断粒子死亡，更新其存活标志live为O,跳转至步骤4 ；
[0040]3.4:如果步骤3.1中的最小距离为平均自由程dw，进行粒子在权重窗中的分裂和赌(分裂和赌为方差减少技巧，当粒子进入重要的相空间区域，则进行分裂。虽然粒子权重变小，但数目增多，使抽样结果更加精确，反之进入不重要的相空间区域，它们以一定的概率被杀死，以节省运算时间)，判断结果若为粒子死亡，则更新存活标志live为0，跳转至步骤4，否则其余三个距离均减去平均自由程dw，并回到步骤3.2 ；
[0041]3.5:如果步骤3.1中的最小距离为到达边界的距离dls，更新所在体的编号icl，跳转至步骤4 ；
[0042]3.6:如果步骤3.1中的最小距离为碰撞距离pmf，则更新碰撞标识NPl为I或2 (I代表中子碰撞，2代表光子碰撞)，跳转至步骤4。
[0043]步骤4:在CPU端统计死亡的粒子，对存活的粒子进行汇总，统计有多少个碰撞事件要发生(其中电子碰撞事件的数目等于电子的个数)，重新划分并行的线程数目，GPU每个线程处理一个碰撞事件；
[0044]在步骤4对碰撞事件的处理中，如果粒子是电子碰撞事件，则包含以下子步骤:
[0045]4.1.1:判断标识El是否为I (El初始设置为I)，如为I则计算电子在当前体的能量损失率qs及小步步数ns，否则略过本步至下一步骤4.1.2 ;
[0046]4.1.2:针对电子，分别计算时间截断距离dtc，小步步长pmf，到体边界的距离dls ；
[0047]4.1.3:找到步骤4.1.2中计算的三个距离中的最小距离，计算电子在通过该最小距离后新的方位角uuu, vvv, WWW,更新当前电子的能量erg和能群编号ngp,如有二次粒子产生，得到二次粒子的能量、方位角、权重等信息；
[0048]4.1.5:根据步骤4.1.3，如果
[0049]I).最小距离为时间截断距离dtc，粒子死亡，更新其存活标志live为0，跳转到步骤5 ;
[0050]2).最小距离为到达边界的距离dls，判断粒子偏转后是否还在该体内，如果不在，则粒子必然到达新的体，判断该新体的编号icl，更新判断标识El为1，跳转到步骤5 ;
[0051]3).最小距离为小步步长pmf,小步剩余步数rns减I (rns初始为步数ns),若小步剩余步数rns为0，更新判断标识El为1，跳转到步骤5 ；
[0052]在步骤4中，如果粒子是中子碰撞事件，则包含以下子步骤:
[0053]4.2.1:采用离散型随机变量抽样法抽样中子与哪种核素碰撞，判断是否采用热力学方法进行处理；
[0054]4.2.2:如果满足以下四个条件，则产生光子；
[0055]I).处理的粒子中包含光子；
[0056]2).该体中光子的重要性(Particle cell importance)不为O,而且没有设置光子产生的偏置处理(即按照正常的反应抽样光子的产生)或者进行了光子产生的偏置处理；
[0057]3).产生光子的截面库存在；
[0058]4).当前时间不超过光子的截断时间。
[0059]其中产生光子的数目和权重依据如下步骤:
[0060]1.)如果没有进行光子产生的偏置设置，则:
[0061]1.1).首先得到产生光子总权重Wp:当前中子权重乘以光子产生截面与总截面之比；
[0062]1.2).光子数目Np = Wp / 51+1，Wi为光子权重临界值(Np最大
[0063]值设为10)；
[0064]1.3).单个光子权重为Wp / Np。
[0065]2.)如果进行了光子产生的偏置设置，则:
[0066]1.1).产生的光子数目设置为I ;
[0067]1.2).光子权重为:当前中子权重乘以偏置反应总截面与所有反应的总截面之比；
[0068]1.3).抽样产生光子的偏置反应。
[0069]依据产生每个光子的反应类型，并据此抽样得到每个光子的出射角；[0070]4.2.3:如果属于裂变源，则处理裂变部分，并记录裂变点的位置及其他信息；
[0071]4.2.4:如果满足以下两个条件，则进行模拟俘获处理的判断:
[0072]I).第一权重截断为0，即隐式俘获关闭，并且只在粒子越过界面时进行权重窗处理；
[0073]2).中子能量小于模拟俘获的截断能量；
[0074]如果抽样确实发生模拟俘获(如果是kcode模式，裂变反应也当作模拟俘获进行处理)，中子死亡，跳转至步骤5 ；
[0075]4.2.5:如果没有进行模拟俘获的处理，则进行隐式俘获处理，即把粒子权重按照是否发生俘获的概率分为俘获部分和非俘获部分，俘获部分被强制杀死，非俘获部分继续参加其它反应；
[0076]4.2.6:如果经4.2.4步骤和4.2.5步骤后粒子依然存活，则处理中子的弹性或非弹性碰撞；步骤4.2. 6包含以下子步骤:
[0077]I).判断靶核的原子量是否小于1.5或者中子能量是否小于400kT，如果两个条件至少满足一个，则抽样靶核的速度，得到中子相对于靶核的相对能量；
[0078]2).处理弹性或者非弹性碰撞反应，得到反应后的中子数cmult、质心C系中的能量 colout (I, k)和极角 colout (2, k)、延迟时间 colout (3, k)等；
[0079]3).将中子的速度和方位角转换为实验室L系；
[0080]在步骤4中，如果粒子是光子碰撞事件，包含以下子步骤:
[0081]4.3.1:抽样碰撞核素；
[0082]4.3.2:如果光子能量大于简单处理方式的最低临界值，则采用简单物理处理方式，否则跳转步骤4.3.8，采用详尽处理方式；
[0083]4.3.3:如果允许光子产生电子，进行光电效应的偏置处理，光子的一部分权重产生光电子(不考虑荧光效应)，另一部分继续下一步的输运；
[0084]光电效应的具体处理步骤为:
[0085]I).光电子的权重为光子权重乘以光电效应截面与总截面的比；
[0086]2).将光子的能量赋值为光电子的能量(忽略电子对结合能)；
[0087]3).计算光电子的极角；
[0088]4.3.5:进行模拟俘获处理，如果发生俘获，则此次光子碰撞处理完毕；
[0089]4.3.6:进行隐式俘获；
[0090]4.3.7:抽样碰撞类型进行处理，分别为非相干散射和电子对效应；
[0091]非相干散射的具体处理方式为:微分散射截面Κ(α, μ)由klein_nishina公式给出，
[0092]Κ(α, μ)?μ = 3?02(α ' / α2[α' / α+α / α ' +μ 2-1]?μ
[0093]抽样散射后的光子方向和能量,其中μ为偏转角余弦，α和α '分别为入射和出射的光子能量，α ' =α / [1+α (1-μ)]，(ι为经典电子半径2.817938X10_13cm。
[0094]如果康普顿反弹电子能量大于电子截断能量，并且程序允许电子输运，记录反弹电子信息，否则采用TTB (thick-tatget bremsstrahlung)即厚祀韧致福射模型近似进行处理。
[0095]电子对效应的具体处理方式为:[0096]如果允许电子输运，获得正负电子的能量和极角记录电子信息，否则采用TTB近似进行处理，特别地如果不允许正电子的输运，或者正电子的能量小于截断能量，该正电子转换为一对方向相反的光子。
[0097]4.3.8:采用详尽处理方式，直接抽样碰撞类型，分别为光电效应，电子对效应，非相干散射(康普顿散射)和相干散射。
[0098]光电效应的处理方式为:相较于简单处理方式增加了对荧光光子的处理，具体方法为:当核电荷数小于12时,无突光光子,大于12小于30时,一个突光光子,大于30时两个突光光子。
[0099]电子对效应的处理与简单处理方式相同。
[0100]非相干散射(康普顿散射)的处理方式是:在简单处理方式的基础上增加了对散射微分截面的修正。
[0101]相干散射(汤姆逊散射)的处理方式为:光子能量不变，按照汤姆逊散射截面抽样光子散射角。
[0102]步骤5、判断粒子是否死亡，如果存活，更新粒子的信息，否则跳转至步骤6。
[0103]步骤6、判断是否有二次粒子产生，如有，将二次粒子的整数类型的信息(包括二次粒子的编号，分身数目npa、所在体中的编号icl、种类ipt、截面库编号iex等)存储在多维数组GPBLCM中；实数类型的信息(包括坐标XXX, yyy, zzz，方位角uuu, vvv, www,能量erg,及权重wgt等)存储在数组JPBLCM中。
[0104]步骤7、在CPU端汇总所有粒子的信息，统计需要进行处理的粒子数目。
[0105]步骤8、判断运行时间是否大于指定时间，如果运行时间大于指定时间，输出统计量，释放CPU和GPU的内存空间，方法结束。
[0106]步骤9、如果运动时间不大于指定时间，判断要处理的粒子数目是否小于设定的最小值，如果小于该值，回到步骤2，如果大于该值，回到步骤3。
[0107]以上的详细描述通过使用示意图、流程图和/或示例，已经阐述了根据本发明的基于GPU的多粒子输运仿真方法的实施例。在这种示意图、流程图和/或示例包含一个或多个功能和/或操作的情况下，本领域技术人员应理解，这种示意图、流程图或示例中的每一功能和/或操作可以通过各种结构、硬件、软件、固件或实质上它们的任意组合来单独和/或共同实现。本领域技术人员应认识到，这里所公开的实施例的一些方面在整体上或部分地可以等同地实现在集成电路中，实现为在一台或多台计算机上运行的一个或多个计算机程序(例如，实现为在一台或多台计算机系统上运行的一个或多个程序)，实现为固件，或者实质上实现为上述方式的任意组合，并且本领域技术人员根据本公开，将具备设计电路和/或写入软件和/或固件代码的能力。此外，本领域技术人员将认识到，本公开所述主题的机制能够作为多种形式的程序产品进行分发，并且无论实际用来执行分发的信号承载介质的具体类型如何，本公开所述主题的示例性实施例均适用。信号承载介质的示例包括但不限于:可记录型介质，如软盘、硬盘驱动器、紧致盘(⑶)、数字通用盘(DVD)、数字磁带、计算机存储器等；以及传输型介质，如数字和/或模拟通信介质(例如，光纤光缆、波导、有线通信链路、无线通信链路等)。
[0108]虽然已参照几个典型实施例描述了本发明，但应当理解，所用的术语是说明和示例性、而非限制性的术语。由于本发明能够以多种形式具体实施而不脱离发明的精神或实质，所以应当理解，上述实施例不限于任何前述的细节，而应在随附权利要求所限定的精神和范围内广泛地解释，因此落入权利要求或其等效范围内的全部变化和改型都应为随附权利要求所涵盖。
【权利要求】
1.一种基于GPU的多粒子输运仿真方法，包括以下步骤: a、使用蒙特卡洛算法对粒子输运建模，建立包括多种几何体的空间，向CPU和GPU分别分配内存空间，在CPU中进行数据初始化，将初始化后的数据从CPU的内存空间拷贝到GPU的内存空间中； b、在CPU中从源产生多个粒子，所述多个粒子包括光子、中子和电子，每个粒子具有粒子信息，包括粒子的种类、位置、速度、权重及所在几何体的编号，将粒子信息连同粒子数目储存在相应数组中； C、判断光子和中子数目之和是否大于设定值，如果数目之和大于设定值，则根据该数目之和确定要使用的线程数目，并将每个中子或光子分配给GPU中的一个线程进行处理，计算前进距离以及判断粒子是否死亡；如果数目之和不大于设定值，则由CPU处理每个中子或光子，计算前进距离以及判断粒子是否死亡； d、由CPU依据中子、光子和电子的碰撞事件数目重新划分线程数目，并且GPU中每个线程处理一个碰撞事件； e、判断粒子是否死亡，如果继续存活，将粒子的信息更新； f、判断是否产生二次粒子，如果产生，将二次粒子的信息存储； g、由CPU重新统计所有存活的粒子信息，统计需要处理的粒子的数目； h、判断运行时间是否大于设定值，如果大于设定值，则输出统计量并释放CPU和GPU的内存空间，所述方法结束； 1、如果运行时间不大于设定值，则判断要处理的粒子的数目是否小于设定的最小值；如果小于最小值，则返回步骤b，否则返回步骤c ； 其中GPU包括多个流处理器，所述多个流处理器通过多个线程对粒子输运仿真进行并行处理，并且所述多个流处理器对粒子的输运仿真和碰撞事件的处理是同步进行的。
2.根据权利要求1所述的方法，其中步骤c还包括: 对于处理的中子或光子，获得时间截断距离、到达边界距离、平均自由程以及碰撞距离； 计算以上距离之中的最小距离； 如果最小距离是时间截断距离，则确定该中子或光子死亡； 如果最小距离是达到边界距离，则更新该中子或光子的所在几何体的编号； 如果最小距离是平均自由程，则用其余三个距离减掉平均自由程，并返回到步骤c的开始； 如果最小距离是碰撞距离，则判断碰撞是中子碰撞还是光子碰撞。
3.根据权利要求1所述的方法，其中在步骤e中如下处理碰撞事件: 如果碰撞事件为电子碰撞事件，则计算电子在上一个最小距离中的能量损失，将产生的二次粒子的数目和状态进行记录，并返回给CPU ； 如果碰撞事件为光子碰撞事件，则抽样碰撞核素，处理光子反应，将产生的二次粒子的数目和状态进行记录，并返回给CPU ； 如果碰撞事件为中子碰撞事件，则抽样碰撞核素，处理中子反应，将产生的二次粒子的数目和状态进行记录，并返回给CPU。
【文档编号】G06F17/50GK103955567SQ201410146948
【公开日】2014年7月30日 申请日期:2014年4月10日 优先权日:2014年4月10日
【发明者】杨磊, 张勋超, 高笑菲, 齐记, 付芬, 张雅玲, 张智磊 申请人:中国科学院近代物理研究所