CN111782384A - 一种基于精细中子时空动力学格子Boltzmann方法的GPU加速方法 - Google Patents

一种基于精细中子时空动力学格子Boltzmann方法的GPU加速方法 Download PDF

Info

Publication number
CN111782384A
CN111782384A CN201910264849.2A CN201910264849A CN111782384A CN 111782384 A CN111782384 A CN 111782384A CN 201910264849 A CN201910264849 A CN 201910264849A CN 111782384 A CN111782384 A CN 111782384A
Authority
CN
China
Prior art keywords
neutron
gpu
representing
energy group
kth
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN201910264849.2A
Other languages
English (en)
Other versions
CN111782384B (zh
Inventor
马宇
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
National Sun Yat Sen University
Original Assignee
National Sun Yat Sen University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by National Sun Yat Sen University filed Critical National Sun Yat Sen University
Priority to CN201910264849.2A priority Critical patent/CN111782384B/zh
Publication of CN111782384A publication Critical patent/CN111782384A/zh
Application granted granted Critical
Publication of CN111782384B publication Critical patent/CN111782384B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T1/00—General purpose image data processing
    • G06T1/20—Processor architectures; Processor configuration, e.g. pipelining
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F9/00—Arrangements for program control, e.g. control units
    • G06F9/06—Arrangements for program control, e.g. control units using stored programs, i.e. using an internal store of processing equipment to receive or retain programs
    • G06F9/46—Multiprogramming arrangements
    • G06F9/50—Allocation of resources, e.g. of the central processing unit [CPU]
    • G06F9/5005—Allocation of resources, e.g. of the central processing unit [CPU] to service a request
    • G06F9/5027—Allocation of resources, e.g. of the central processing unit [CPU] to service a request the resource being a machine, e.g. CPUs, Servers, Terminals
    • G06F9/505—Allocation of resources, e.g. of the central processing unit [CPU] to service a request the resource being a machine, e.g. CPUs, Servers, Terminals considering the load
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02E—REDUCTION OF GREENHOUSE GAS [GHG] EMISSIONS, RELATED TO ENERGY GENERATION, TRANSMISSION OR DISTRIBUTION
    • Y02E30/00—Energy generation of nuclear origin
    • Y02E30/30—Nuclear fission reactors

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Software Systems (AREA)
  • Mathematical Physics (AREA)
  • Computational Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • General Engineering & Computer Science (AREA)
  • Data Mining & Analysis (AREA)
  • Operations Research (AREA)
  • Databases & Information Systems (AREA)
  • Algebra (AREA)
  • Monitoring And Testing Of Nuclear Reactors (AREA)

Abstract

本发明提出了一种基于精细中子时空动力学格子Boltzmann方法的GPU加速方法，属于核反应堆设计及安全分析技术领域。所述方法为：在多核GPU设备上采用格子Boltzmann方法对非均匀堆芯进行扩散计算得到稳态条件下核反应堆不同空间位置、不同能群中的中子注量率分布以及缓发中子先驱核浓度分布；然后，采用格子Boltzmann方法求解准静态时空动力学方程，得到不同时刻下中子标量注量率的时空分布情况以及缓发中子先驱核浓度在不同时刻不同位置的分布。本发明所述加速方法能降低反应堆大规模并行计算的难度，同时大幅度提高反应堆中子时空动力学问题的计算速度；消除了传统方法实现复杂、成本高、加速效率不高的缺陷。

Description

一种基于精细中子时空动力学格子Boltzmann方法的GPU加速 方法
技术领域
本发明涉及一种基于精细中子时空动力学格子Boltzmann方法的GPU加速方法，属于核反应堆设计及安全分析技术领域。
背景技术
在当前核反应堆设计及安全分析中，常常需要研究如控制棒移动等外部扰动条件下核反应堆精细功率分布的变化，以验证反应堆设计在不同瞬态条件下是否能满足安全指标，基于此，需要对反应堆运行状态进行精细瞬态时空动力学分析。随着计算机技术及数值模拟手段的不断发展，全堆芯大规模中子扩散精细时空动力学计算成为了主流的研究方法之一。精细中子扩散计算能够获得堆芯不同时刻不同空间位置的中子注量率分布，随着网格数的增加，计算精度不断提高，最终接近真解，满足现代反应堆高分辨率和高精度的计算需求，实现高保真反应堆中子物理模拟。
但是对于高保真精细时空动力学计算而言，由于高计算精度对于网格密度的强烈需求，庞大的网格数量以及巨额的计算时间限制了大规模全堆芯中子时空动力学计算在核反应堆中的进一步应用与发展。为了减少时空动力学循环计算次数，目前国际上采用准静态方法来处理时空动力学的时间演化过程，降低了待求时间步，减少了循环计算的次数，但是即使如此，为了满足计算空间精度，空间网格数需求依旧很高，单步计算时间以久很长。针对这种情况，国际上目前采用大规模并行计算手段来缩短单步计算时间，将整体空间区域划分为多个子区域同步并行计算，降低了单步计算所需时间。但是当前的CPU并行计算对成本需求极高，多核计算机配置价格随着核数的增加价格提升也十分明显。由于利用了GPU的多线程特性，GPU并行计算能实现成千上万线程级的并行计算，同时相对于CPU计算只需要较为低廉的价格成本。虽然GPU计算是一种很好的大规模并行计算手段，但是当前的计算方法对于GPU并行计算的实现并不容易，尤其是有限元方法，由于其实现较为困难且基于矩阵求解，其GPU实现困难且效率不是很高。
发明内容
本发明为了解决大规模GPU并行计算实现难度大，并行效率低的问题，提出了一种精细中子时空动力学格子Boltzmann方法的GPU加速技术，该方法将格子Boltzmann方法应用于中子时空动力学问题并将其应用于GPU加速计算，在不改变网格密度等条件的前提下，能够极大提高计算速度，快速模拟核反应堆内中子注量率的时空分布规律，克服传统数值方法在GPU加速计算中难于实现且加速效果不明显的缺点，使反应堆精细中子时空动力学计算在核反应堆设计及安全分析领域发挥更为重要的作用。为了达到以上需求，本发明采用如下技术方案：
一种基于精细中子时空动力学格子Boltzmann方法的GPU加速方法，如图1至3所示，所述GPU加速方法包括：
步骤一、核反应堆稳态多群扩散方程计算，采用格子Boltzmann方法对非均匀堆芯进行扩散计算，获得稳态条件下核反应堆不同空间位置、不同能群中的中子注量率分布以及缓发中子先驱核浓度分布；
步骤二、将反应堆系统调整到临界状态，并通过临界状态相关参数计算反应堆伴随中子扩散方程，然后将所述反应堆伴随中子扩散方程进行归一化计算，获取临界初始状态下反应堆内伴随中子注量率在空间不同位置上的归一化分布情况；
步骤三、求解准静态中子时空动力学方程，获得不同时刻下中子标量注量率的时空分布情况，以及缓发中子先驱核浓度在不同时刻不同位置的分布。
进一步地，步骤一所述获得稳态条件下核反应堆不同空间位置、不同能群中的中子注量率分布以及缓发中子先驱核浓度分布的具体过程包括：
第一步、从几何卡和材料卡读入核反应堆几何参数、材料配置、不同材料的中子宏观截面参数以及动力学参数；对计算区域采用有限的离散网格进行离散化处理；按照GPU设备数量及性能，将整体计算区域划分为多个子区域，每一个子区域由一个独立的GPU设备进行计算；在各个GPU设备中分配内存并初始化，所述计算项目获得的相关参数包括裂变源项、中子分布函数、中子注量率、缓发中子先驱核、中子宏观截面以及有效增殖因数keff；
第二步、根据第一步计算获得的相关参数和网络信息，采用OpenMP在CPU中开启多个线程同步调用所有的GPU设备并在所有GPU设备中同步进行基于格子Boltzmann方法的中子扩散计算，得到不同能群在不同空间位置的中子注量率分布，采用计算公式如下：
Figure BDA0002016460590000021
Figure BDA0002016460590000022
Figure BDA0002016460590000023
Figure BDA0002016460590000024
[φg(r)]k＝∑i[fi,g(r)]k；
通过上述计算即可获得各个GPU上不同子区域空间各个节点的中子注量率分布；其中，i表示格子Boltzmann方向；g表示能群编号；G表示总的能群数量；k表示GPU设备编号；K表示总的GPU设备数；r表示空间位置；ei表示格子Boltzmann分布函数速度；ε表示努森数，是中子平均自由程与特征尺度的比值；[fi,g(r+εei)]k表示第k个GPU子区域中r+eiε处i方向g能群内的中子分布函数；[fi,g(r)]k表示第k个GPU子区域中r处i方向g能群内的中子分布函数；[τ(r)]k表示第k个GPU子区域中r处的无量纲松弛时间因子；[fi,g(r)]k表示第k个GPU子区域中r处i方向g能群内的中子分布函数；
Figure BDA0002016460590000031
表示第k个GPU子区域中r处i方向g能群内的中子平衡态分布函数；[Ri,g(r)]k表示第k个GPU子区域中r处i方向g能群内的中子源分布函数；[Rg(r)]k表示第k个GPU子区域中r处g能群内的中子总源项；[Σg'-g(r)]k表示第k个GPU子区域中r处g’能群到g能群的宏观中子散射截面；[φg'(r)]k表示第k个GPU子区域中r处g’能群的中子标量注量率；β表示缓发中子份额总数；χg表示g能群的中子能谱；[νΣf,g'(r)]k表示第k个GPU子区域中r处g’能群的中子产生截面；[Σr,g(r)]k表示第k个GPU子区域中r处g能群的宏观中子移出截面；d表示缓发中子先驱核组编号；
Figure BDA0002016460590000032
表示第g能群第d缓发中子组的中子能谱；λd表示第d组缓发中子衰变常数；[Cd(r)]k表示第k个GPU子区域中r处第d组缓发中子先驱核浓度；
Figure BDA0002016460590000033
表示格子Boltzmann方法中i方向的方向权重；[Dg(r)]k表示第k个GPU子区域中r处g能群的中子扩散系数；cs表示格子Boltzmann声速；keff表示反应堆堆芯的有效增殖因数；
第三步、将各个GPU内存中对应于子区域内边界节点上，指向交界面外的中子注量率分布函数拷贝回CPU，再通过CPU将所述中子注量率分布函数拷贝到原GPU的相邻子区域内边界向内分布函数中，由此完成不同子区域之间的数据交换，具体公式如下：
Figure BDA0002016460590000034
其中，g表示能群编号；G表示总的能群数量；k表示GPU设备编号；K表示总的GPU设备数；[fi,g(r)]k表示第k个GPU子区域中r处i方向g能群内的中子分布函数；
Figure BDA0002016460590000035
表示第k个GPU子区域的ei相邻子区域中r+eiε处i方向g能群内的中子分布函数；
第四步、根据第二步中获得的中子注量率分布，在各个GPU不同子区域之中通过计算获得稳态条件下的缓发中子先驱核密度分布，所采用的计算公式如下：
Figure BDA0002016460590000041
其中，g表示能群编号；G表示总的能群数量；k表示GPU设备编号；K表示总的GPU设备数；d表示缓发中子先驱核组编号；D表示缓发中子先驱核总的组数；[Cd(r)]k表示第k个GPU子区域中r处第d组缓发中子先驱核浓度；βd表示第d组缓发中子份额；λd表示第d组缓发中子衰变常数；[νΣf,g'(r)]k表示第k个GPU子区域中r处g’能群的中子产生截面；keff表示反应堆堆芯的有效增殖因数；
通过所述计算公式即可获得稳态下的缓发中子先驱核浓度在空间不同节点的分布情况，由于缓发中子先驱核和中子注量率是线性局部关系，因此不同子区域之间不需要进行缓发中子先驱核的数据交换。
进一步地，步骤二所述获取临界初始状态下反应堆内伴随中子注量率在空间不同位置上的归一化分布情况的具体过程包括：
第1步：在各个GPU内存中通过调整反应堆系统相关参数使反应堆达到临界状态，所述调整反应堆系统相关参数的具体公式如下：
Figure BDA0002016460590000042
其中，g表示能群编号；G表示总的能群数量；k表示GPU设备编号；K表示总的GPU设备数；[νΣf,g(r)]k表示第k个GPU子区域中r处g能群的中子产生截面；keff表示反应堆堆芯的有效增殖因数；
第2步：基于与步骤1相同的参数条件，在各个GPU对应子区域中计算伴随中子扩散方程；计算所述伴随中子扩散方程的具体公式如下：
Figure BDA0002016460590000043
Figure BDA0002016460590000051
Figure BDA0002016460590000052
Figure BDA0002016460590000053
Figure BDA0002016460590000054
通过以上方程即可得到伴随中子注量率在反应堆不同空间位置的分布，其中，i表示格子Boltzmann方向；g表示能群编号；G表示总的能群数量；k表示GPU设备编号；K表示总的GPU设备数；r表示空间位置；ei表示格子Boltzmann分布函数速度；ε表示努森数，中子平均自由程与特征尺度的比值；
Figure BDA0002016460590000055
表示第k个GPU子区域中r+eiε处i方向g能群内的伴随中子分布函数；
Figure BDA0002016460590000056
表示第k个GPU子区域中r处i方向g能群内的伴随中子分布函数；[τ(r)]k表示第k个GPU子区域中r处的无量纲松弛时间因子；
Figure BDA0002016460590000057
表示第k个GPU子区域中r处i方向g能群内的伴随中子分布函数；
Figure BDA0002016460590000058
表示第k个GPU子区域中r处i方向g能群内的伴随中子平衡态分布函数；
Figure BDA0002016460590000059
表示第k个GPU子区域中r处i方向g能群内的伴随中子源分布函数；
Figure BDA00020164605900000510
表示第k个GPU子区域中r处g能群内的伴随中子总源项；[Σg-g'(r,t0)]k表示第k个GPU子区域中r处g能群t0时刻到g’能群的宏观中子散射截面；
Figure BDA00020164605900000511
表示第k个GPU子区域中r处g能群的伴随中子标量注量率；χg’表示g’能群的中子能谱；[νΣf,g(r,t0)]k表示第k个GPU子区域中r处t0时刻g能群的中子产生截面；[Σr,g(r,t0)]k表示第k个GPU子区域中r处t0时刻g能群的宏观中子移出截面；
Figure BDA00020164605900000512
表示格子Boltzmann方法方向权重；[Dg(r)]k表示第k个GPU子区域中r处g能群的中子扩散系数；cs表示格子Boltzmann声速；keff表示反应堆堆芯的有效增殖因数；通过上述参量即可表达所述伴随中子扩散方程；
第3步：将所有GPU上的中子伴随注量率拷贝回CPU，通过归一化条件将伴随中子注量率进行重构，进而获得归一化的中子伴随注量率，随后引入相应的反应性使反应堆启动；所述重构的具体公式如下：
Figure BDA00020164605900000513
其中，g表示能群编号；G表示总的能群数量；vg表示第g能群中子速度；
Figure BDA0002016460590000061
表示r处g能群的伴随中子标量注量率；
Figure BDA0002016460590000062
表示r处g能群t时刻的中子注量率形状函数。
由此得到归一化的中子伴随注量率，随后引入相应的反应性使反应堆启动
进一步地，步骤三所述获得不同时刻下中子标量注量率的时空分布情况，以及缓发中子先驱核浓度在不同时刻不同位置的分布的具体过程包括：
步骤1、设置大时间步Δtbig、中时间步Δtmiddle以及小时间步Δtlittle；在所有不同GPU设备的不同子区域中的大时间步Δtbig尺度上预估t+Δtbig时刻的中子注量率并计算Δtbig时间步内的中子注量率形状函数，具体计算公式如下：
Figure BDA0002016460590000063
Figure BDA0002016460590000064
Figure BDA0002016460590000065
Figure BDA0002016460590000066
[φg(r,t)]k＝∑i[fi,g(r,t)]k；
Figure BDA0002016460590000067
其中，i表示格子Boltzmann方向；g表示能群编号；G表示总的能群数量；k表示GPU设备编号；K表示总的GPU设备数；D表示总缓发中子先驱核组数；r表示空间位置；ei表示格子Boltzmann分布函数速度；Δtbig表示大时间步长；vg表示第g能群中子速度；[fi,g(r+eiΔtbig,t+Δtbig)]k表示第k个GPU子区域中r+eiΔtbig处t+Δtbig时刻i方向g能群内的中子分布函数；[fi,g(r,t)]k表示第k个GPU子区域中r处t时刻i方向g能群内的中子分布函数；[τ(r,t)]k表示第k个GPU子区域中r处t时刻的无量纲松弛时间因子；[fi,g(r,t)]k表示第k个GPU子区域中r处t时刻i方向g能群内的中子分布函数；
Figure BDA0002016460590000071
表示第k个GPU子区域中r处t时刻i方向g能群内的中子平衡态分布函数；[Ri,g(r,t)]k表示第k个GPU子区域中r处t时刻i方向g能群内的中子源分布函数；[Rg(r,t)]k表示第k个GPU子区域中r处t时刻g能群内的中子总源项；[Σg'-g(r,t)]k表示第k个GPU子区域中r处t时刻g’能群到g能群的宏观中子散射截面；[φg'(r,t)]k表示第k个GPU子区域中r处t时刻g’能群的中子标量注量率；β表示缓发中子份额总数；βd表示第d组缓发中子份额；χg表示g能群的中子能谱；[νΣf,g'(r,t)]k表示第k个GPU子区域中r处t时刻g’能群的中子产生截面；[Σr,g(r,t)]k表示第k个GPU子区域中r处t时刻g能群的宏观中子移出截面；d表示缓发中子先驱核组编号；
Figure BDA0002016460590000072
表示第g能群第d缓发中子组的中子能谱；λd表示第d组缓发中子衰变常数；[Cd(r,t)]k表示第k个GPU子区域中r处t时刻第d组缓发中子先驱核浓度；[Cd(r,t+Δtbig)]k表示第k个GPU子区域中r处t+Δtbig时刻第d组缓发中子先驱核浓度；[νΣf,g'(r,t+Δtbig)]k表示第k个GPU子区域中r处t+Δtbig时刻g’能群的中子产生截面；[φg'(r,t+Δtbig)]k表示第k个GPU子区域中r处t+Δtbig时刻g’能群的中子标量注量率；
Figure BDA0002016460590000073
表示格子Boltzmann方法方向权重；[Dg(r,t)]k表示第k个GPU子区域中r处t时刻g能群的中子扩散系数；cs表示格子Boltzmann声速；
步骤2、将不同GPU上每一个子区域内边界节点上位于t+Δtbig时刻的向外中子注量率分布函数拷贝回CPU，再将其拷贝回相邻GPU子区域内边界节点上向内中子注量率分布函数中，完成不同子区域之间的数据交换，所述数据交换过程应用的具体公式如下：
Figure BDA0002016460590000074
其中，g表示能群编号；G表示总的能群数量；k表示GPU设备编号；K表示总的GPU设备数；[fi,g(r)]k表示第k个GPU子区域中r处i方向g能群内的中子分布函数；
Figure BDA0002016460590000075
表示第k个GPU子区域的ei相邻子区域中r+eiε处i方向g能群内的中子分布函数；
步骤3、在所有GPU对应的子区域中计算中子注量率在t～t+Δtbig时间段内的形状函数，所述形状函数的计算公式如下：
Figure BDA0002016460590000081
其中，g表示能群编号；G表示总的能群数量；k表示GPU设备编号；K表示总的GPU设备数；
Figure BDA0002016460590000082
表示第k个GPU子区域中r处g能群t时刻的中子注量率形状函数；vg表示第g能群中子速度；
Figure BDA0002016460590000083
表示第k个GPU子区域中r处g能群的伴随中子标量注量率；[φg(r)]k表示第k个GPU子区域中r处g能群的中子标量注量率；
由以上步骤可以获得该大时间步内的中子注量率形状函数；
步骤4、将所有GPU设备上不同子区域的变量拷贝回CPU，在CPU上进行点堆方程计算；在每一个中时间步Δtmiddle上线性插值中子注量率形状函数；在每一个中时间步Δtmiddle内，以小时间步Δtlittle计算点堆方程，具体公式如下：
Figure BDA0002016460590000084
Figure BDA0002016460590000085
Figure BDA0002016460590000086
Figure BDA0002016460590000087
Figure BDA0002016460590000088
Figure BDA0002016460590000089
其中，n(t)表示t时刻中子注量率的幅值函数；d表示缓发中子先驱核组编号；g表示能群编号；G表示总能群数；k表示GPU设备编号；K表示总的GPU设备数；D表示总缓发中子先驱核组数；md(t)表示t时刻第d组缓发中子先驱核的幅值函数；ρ(t)表示t时刻的反应性；
Figure BDA0002016460590000091
表示t时刻的点堆方程缓发中子总份额；
Figure BDA0002016460590000092
表示t时刻的点堆方程d组缓发中子份额；Λ(t)表示t时刻的点堆方程中子代时间；λd表示第d组缓发中子衰变常数；
Figure BDA0002016460590000093
表示r处g能群的伴随中子标量注量率；Dg(r,t)表示r处t时刻g能群的中子扩散系数；
Figure BDA0002016460590000094
表示r处g能群t时刻的中子注量率形状函数；Σr,g(r,t)表示r处t时刻g能群的宏观中子移出截面；χg表示g能群的中子能谱；
Figure BDA0002016460590000095
表示第g能群第d缓发中子组的中子能谱；β表示缓发中子份额总数；νΣf,g'(r,t)表示r处t时刻g’能群的中子产生截面；βd表示第d组缓发中子份额；vg表示第g能群中子速度；
步骤5、利用步骤3和步骤4的获得的中子注量率形状函数和幅值函数分别在所有GPU的不同子区域上修正中子注量率在t+Δtbig的分布，并计算缓发中子先驱核在t+Δtbig的分布，进而得到t+Δtbig时刻中子注量率与缓发中子先驱核的时空分布，具体公式如下：
Figure BDA0002016460590000096
Figure BDA0002016460590000097
其中，d表示缓发中子先驱核组编号；g表示能群编号；G表示总能群数；k表示GPU设备编号；K表示总的GPU设备数；D表示总缓发中子先驱核组数；
Figure BDA0002016460590000098
表示第k个GPU子区域中r处g能群t时刻的中子注量率形状函数；[φg(r,t)]k表示第k个GPU子区域中r处t时刻g能群的中子标量注量率；n(t)表示t时刻中子注量率的幅值函数；λd表示第d组缓发中子衰变常数；[Cd(r,t)]k表示第k个GPU子区域中r处t时刻第d组缓发中子先驱核浓度；[Cd(r,t+Δtbig)]k表示第k个GPU子区域中r处t+Δtbig时刻第d组缓发中子先驱核浓度；[νΣf,g'(r,t+Δtbig)]k表示第k个GPU子区域中r处t+Δtbig时刻g’能群的中子产生截面；[φg'(r,t+Δtbig)]k表示第k个GPU子区域中r处t+Δtbig时刻g’能群的中子标量注量率；βd表示第d组缓发中子份额；
由此可以得到t+Δtbig时刻中子注量率与缓发中子先驱核的时空分布；
步骤6、重复步骤步骤1至步骤5直到达到需要计算的时空中子动力学时间。
本发明有益效果：
本发明提出的一种基于精细中子时空动力学格子Boltzmann方法的GPU加速方法消除了传统方法实现复杂、成本高、加速效率不高的缺陷，该方法在极低的计算成本与极简单的实现步骤下大幅度提高计算速度，对核反应堆大规模中子时空动力学数值模拟的发展及应用有重要作用；同时，所述基于精细中子时空动力学格子Boltzmann方法的GPU加速方法能够实现任意数量GPU显卡设备的加速策略，能够实现在不同性能显卡之间负载平衡，能够灵活充分利用计算机资源，降低计算成本，提高计算效率与经济性。
附图说明
图1是精细中子时空动力学瞬态计算流程。
图2是堆芯精细稳态计算流程。
图3是精细中子时空动力学瞬态计算的GPU加速流程图。
图4是TWIGL基准题稳态快群中子注量率对比图。
图5是TWIGL基准题稳态热群中子注量率对比图。
图6是TWIGL基准题功率变化对比图。
图7是TWIGL基准题不同GPU配置的加速比。
具体实施方式
下面结合具体实施例对本发明做进一步说明，但本发明不受实施例的限制。
实施例1：
一种基于精细中子时空动力学格子Boltzmann方法的GPU加速方法，所述GPU加速方法包括：
步骤一、核反应堆稳态多群扩散方程计算，采用格子Boltzmann方法对非均匀堆芯进行扩散计算，获得稳态条件下核反应堆不同空间位置、不同能群中的中子注量率分布以及缓发中子先驱核浓度分布；
步骤二、将反应堆系统调整到临界状态，并通过临界状态相关参数计算反应堆伴随中子扩散方程，然后将所述反应堆伴随中子扩散方程进行归一化计算，获取临界初始状态下反应堆内伴随中子注量率在空间不同位置上的归一化分布情况；
步骤三、求解准静态中子时空动力学方程，获得不同时刻下中子标量注量率的时空分布情况，以及缓发中子先驱核浓度在不同时刻不同位置的分布。
其中，上述步骤一至步骤三的具体处理过程如下：
步骤一所述获得稳态条件下核反应堆不同空间位置、不同能群中的中子注量率分布以及缓发中子先驱核浓度分布的具体过程包括：
第一步、从几何卡和材料卡读入核反应堆几何参数、材料配置、不同材料的中子宏观截面参数以及动力学参数；对计算区域采用有限的离散网格进行离散化处理；按照GPU设备数量及性能，将整体计算区域划分为多个子区域，每一个子区域由一个独立的GPU设备进行计算；在各个GPU设备中分配内存并初始化，所述计算项目获得的相关参数包括裂变源项、中子分布函数、中子注量率、缓发中子先驱核、中子宏观截面以及有效增殖因数keff；
第二步、根据第一步计算获得的相关参数和网络信息，采用OpenMP在CPU中开启多个线程同步调用所有的GPU设备并在所有GPU设备中同步进行基于格子Boltzmann方法的中子扩散计算，得到不同能群在不同空间位置的中子注量率分布，采用计算公式如下：
Figure BDA0002016460590000111
Figure BDA0002016460590000112
Figure BDA0002016460590000113
Figure BDA0002016460590000114
[φg(r)]k＝∑i[fi,g(r)]k；
通过上述计算即可获得各个GPU上不同子区域空间各个节点的中子注量率分布；其中，i表示格子Boltzmann方向；g表示能群编号；G表示总的能群数量；k表示GPU设备编号；K表示总的GPU设备数；r表示空间位置；ei表示格子Boltzmann分布函数速度；ε表示努森数，是中子平均自由程与特征尺度的比值；[fi,g(r+εei)]k表示第k个GPU子区域中r+eiε处i方向g能群内的中子分布函数；[fi,g(r)]k表示第k个GPU子区域中r处i方向g能群内的中子分布函数；[τ(r)]k表示第k个GPU子区域中r处的无量纲松弛时间因子；[fi,g(r)]k表示第k个GPU子区域中r处i方向g能群内的中子分布函数；
Figure BDA0002016460590000115
表示第k个GPU子区域中r处i方向g能群内的中子平衡态分布函数；[Ri,g(r)]k表示第k个GPU子区域中r处i方向g能群内的中子源分布函数；[Rg(r)]k表示第k个GPU子区域中r处g能群内的中子总源项；[Σg'-g(r)]k表示第k个GPU子区域中r处g’能群到g能群的宏观中子散射截面；[φg'(r)]k表示第k个GPU子区域中r处g’能群的中子标量注量率；β表示缓发中子份额总数；χg表示g能群的中子能谱；[νΣf,g'(r)]k表示第k个GPU子区域中r处g’能群的中子产生截面；[Σr,g(r)]k表示第k个GPU子区域中r处g能群的宏观中子移出截面；d表示缓发中子先驱核组编号；
Figure BDA0002016460590000121
表示第g能群第d缓发中子组的中子能谱；λd表示第d组缓发中子衰变常数；[Cd(r)]k表示第k个GPU子区域中r处第d组缓发中子先驱核浓度；
Figure BDA0002016460590000122
表示格子Boltzmann方法中i方向的方向权重；[Dg(r)]k表示第k个GPU子区域中r处g能群的中子扩散系数；cs表示格子Boltzmann声速；keff表示反应堆堆芯的有效增殖因数；
第三步、将各个GPU内存中对应于子区域内边界节点上，指向交界面外的中子注量率分布函数拷贝回CPU，再通过CPU将所述中子注量率分布函数拷贝到原GPU的相邻子区域内边界向内分布函数中，由此完成不同子区域之间的数据交换，具体公式如下：
Figure BDA0002016460590000123
其中，g表示能群编号；G表示总的能群数量；k表示GPU设备编号；K表示总的GPU设备数；[fi,g(r)]k表示第k个GPU子区域中r处i方向g能群内的中子分布函数；
Figure BDA0002016460590000124
表示第k个GPU子区域的ei相邻子区域中r+eiε处i方向g能群内的中子分布函数；
第四步、根据第二步中获得的中子注量率分布，在各个GPU不同子区域之中通过计算获得稳态条件下的缓发中子先驱核密度分布，所采用的计算公式如下：
Figure BDA0002016460590000125
其中，g表示能群编号；G表示总的能群数量；k表示GPU设备编号；K表示总的GPU设备数；d表示缓发中子先驱核组编号；D表示缓发中子先驱核总的组数；[Cd(r)]k表示第k个GPU子区域中r处第d组缓发中子先驱核浓度；βd表示第d组缓发中子份额；λd表示第d组缓发中子衰变常数；[νΣf,g'(r)]k表示第k个GPU子区域中r处g’能群的中子产生截面；keff表示反应堆堆芯的有效增殖因数；
通过所述计算公式即可获得稳态下的缓发中子先驱核浓度在空间不同节点的分布情况，由于缓发中子先驱核和中子注量率是线性局部关系，因此不同子区域之间不需要进行缓发中子先驱核的数据交换。
步骤二所述获取临界初始状态下反应堆内伴随中子注量率在空间不同位置上的归一化分布情况的具体过程包括：
第1步：在各个GPU内存中通过调整反应堆系统相关参数使反应堆达到临界状态，所述调整反应堆系统相关参数的具体公式如下：
Figure BDA0002016460590000131
其中，g表示能群编号；G表示总的能群数量；k表示GPU设备编号；K表示总的GPU设备数；[νΣf,g(r)]k表示第k个GPU子区域中r处g能群的中子产生截面；keff表示反应堆堆芯的有效增殖因数；
第2步：基于与步骤1相同的参数条件，在各个GPU对应子区域中计算伴随中子扩散方程；计算所述伴随中子扩散方程的具体公式如下：
Figure BDA0002016460590000132
Figure BDA0002016460590000133
Figure BDA0002016460590000134
Figure BDA0002016460590000135
Figure BDA0002016460590000136
通过以上方程即可得到伴随中子注量率在反应堆不同空间位置的分布，其中，i表示格子Boltzmann方向；g表示能群编号；G表示总的能群数量；k表示GPU设备编号；K表示总的GPU设备数；r表示空间位置；ei表示格子Boltzmann分布函数速度；ε表示努森数，中子平均自由程与特征尺度的比值；
Figure BDA0002016460590000137
表示第k个GPU子区域中r+eiε处i方向g能群内的伴随中子分布函数；
Figure BDA0002016460590000141
表示第k个GPU子区域中r处i方向g能群内的伴随中子分布函数；[τ(r)]k表示第k个GPU子区域中r处的无量纲松弛时间因子；[fi,g(r)]k表示第k个GPU子区域中r处i方向g能群内的伴随中子分布函数；
Figure BDA0002016460590000142
表示第k个GPU子区域中r处i方向g能群内的伴随中子平衡态分布函数；
Figure BDA0002016460590000143
表示第k个GPU子区域中r处i方向g能群内的伴随中子源分布函数；
Figure BDA0002016460590000144
表示第k个GPU子区域中r处g能群内的伴随中子总源项；[Σg-g'(r,t0)]k表示第k个GPU子区域中r处g能群t0时刻到g’能群的宏观中子散射截面；
Figure BDA0002016460590000145
表示第k个GPU子区域中r处g能群的伴随中子标量注量率；χg’表示g’能群的中子能谱；[νΣf,g(r,t0)]k表示第k个GPU子区域中r处t0时刻g能群的中子产生截面；[Σr,g(r,t0)]k表示第k个GPU子区域中r处t0时刻g能群的宏观中子移出截面；
Figure BDA0002016460590000146
表示格子Boltzmann方法方向权重；[Dg(r)]k表示第k个GPU子区域中r处g能群的中子扩散系数；cs表示格子Boltzmann声速；keff表示反应堆堆芯的有效增殖因数；通过上述参量即可表达所述伴随中子扩散方程；
第3步：将所有GPU上的中子伴随注量率拷贝回CPU，通过归一化条件将伴随中子注量率进行重构，进而获得归一化的中子伴随注量率，随后引入相应的反应性使反应堆启动；所述重构的具体公式如下：
Figure BDA0002016460590000147
其中，g表示能群编号；G表示总的能群数量；vg表示第g能群中子速度；
Figure BDA0002016460590000148
表示r处g能群的伴随中子标量注量率；
Figure BDA0002016460590000149
表示r处g能群t时刻的中子注量率形状函数。由此得到归一化的中子伴随注量率，随后引入相应的反应性使反应堆启动。
步骤三所述获得不同时刻下中子标量注量率的时空分布情况，以及缓发中子先驱核浓度在不同时刻不同位置的分布的具体过程包括：
步骤1、设置大时间步Δtbig、中时间步Δtmiddle以及小时间步Δtlittle；在所有不同GPU设备的不同子区域中的大时间步Δtbig尺度上预估t+Δtbig时刻的中子注量率并计算Δtbig时间步内的中子注量率形状函数，具体计算公式如下：
Figure BDA0002016460590000151
Figure BDA0002016460590000152
Figure BDA0002016460590000153
Figure BDA0002016460590000154
[φg(r,t)]k＝∑i[fi,g(r,t)]k；
Figure BDA0002016460590000155
其中，i表示格子Boltzmann方向；g表示能群编号；G表示总的能群数量；k表示GPU设备编号；K表示总的GPU设备数；D表示总缓发中子先驱核组数；r表示空间位置；ei表示格子Boltzmann分布函数速度；Δtbig表示大时间步长；vg表示第g能群中子速度；[fi,g(r+eiΔtbig,t+Δtbig)]k表示第k个GPU子区域中r+eiΔtbig处t+Δtbig时刻i方向g能群内的中子分布函数；[fi,g(r,t)]k表示第k个GPU子区域中r处t时刻i方向g能群内的中子分布函数；[τ(r,t)]k表示第k个GPU子区域中r处t时刻的无量纲松弛时间因子；[fi,g(r,t)]k表示第k个GPU子区域中r处t时刻i方向g能群内的中子分布函数；
Figure BDA0002016460590000156
表示第k个GPU子区域中r处t时刻i方向g能群内的中子平衡态分布函数；[Ri,g(r,t)]k表示第k个GPU子区域中r处t时刻i方向g能群内的中子源分布函数；[Rg(r,t)]k表示第k个GPU子区域中r处t时刻g能群内的中子总源项；[Σg'-g(r,t)]k表示第k个GPU子区域中r处t时刻g’能群到g能群的宏观中子散射截面；[φg'(r,t)]k表示第k个GPU子区域中r处t时刻g’能群的中子标量注量率；β表示缓发中子份额总数；βd表示第d组缓发中子份额；χg表示g能群的中子能谱；[νΣf,g'(r,t)]k表示第k个GPU子区域中r处t时刻g’能群的中子产生截面；[Σr,g(r,t)]k表示第k个GPU子区域中r处t时刻g能群的宏观中子移出截面；d表示缓发中子先驱核组编号；
Figure BDA0002016460590000161
表示第g能群第d缓发中子组的中子能谱；λd表示第d组缓发中子衰变常数；[Cd(r,t)]k表示第k个GPU子区域中r处t时刻第d组缓发中子先驱核浓度；[Cd(r,t+Δtbig)]k表示第k个GPU子区域中r处t+Δtbig时刻第d组缓发中子先驱核浓度；[νΣf,g'(r,t+Δtbig)]k表示第k个GPU子区域中r处t+Δtbig时刻g’能群的中子产生截面；[φg'(r,t+Δtbig)]k表示第k个GPU子区域中r处t+Δtbig时刻g’能群的中子标量注量率；
Figure BDA0002016460590000162
表示格子Boltzmann方法方向权重；[Dg(r,t)]k表示第k个GPU子区域中r处t时刻g能群的中子扩散系数；cs表示格子Boltzmann声速；
步骤2、将不同GPU上每一个子区域内边界节点上位于t+Δtbig时刻的向外中子注量率分布函数拷贝回CPU，再将其拷贝回相邻GPU子区域内边界节点上向内中子注量率分布函数中，完成不同子区域之间的数据交换，所述数据交换过程应用的具体公式如下：
Figure BDA0002016460590000163
其中，g表示能群编号；G表示总的能群数量；k表示GPU设备编号；K表示总的GPU设备数；[fi,g(r)]k表示第k个GPU子区域中r处i方向g能群内的中子分布函数；
Figure BDA0002016460590000164
表示第k个GPU子区域的ei相邻子区域中r+eiε处i方向g能群内的中子分布函数；
步骤3、在所有GPU对应的子区域中计算中子注量率在t～t+Δtbig时间段内的形状函数，所述形状函数的计算公式如下：
Figure BDA0002016460590000165
其中，g表示能群编号；G表示总的能群数量；k表示GPU设备编号；K表示总的GPU设备数；
Figure BDA0002016460590000166
表示第k个GPU子区域中r处g能群t时刻的中子注量率形状函数；vg表示第g能群中子速度；
Figure BDA0002016460590000167
表示第k个GPU子区域中r处g能群的伴随中子标量注量率；[φg(r)]k表示第k个GPU子区域中r处g能群的中子标量注量率；
由以上步骤可以获得该大时间步内的中子注量率形状函数；
步骤4、将所有GPU设备上不同子区域的变量拷贝回CPU，在CPU上进行点堆方程计算；在每一个中时间步Δtmiddle上线性插值中子注量率形状函数；在每一个中时间步Δtmiddle内，以小时间步Δtlittle计算点堆方程，具体公式如下：
Figure BDA0002016460590000171
Figure BDA0002016460590000172
Figure BDA0002016460590000173
Figure BDA0002016460590000174
Figure BDA0002016460590000175
Figure BDA0002016460590000176
其中，n(t)表示t时刻中子注量率的幅值函数；d表示缓发中子先驱核组编号；g表示能群编号；G表示总能群数；k表示GPU设备编号；K表示总的GPU设备数；D表示总缓发中子先驱核组数；md(t)表示t时刻第d组缓发中子先驱核的幅值函数；ρ(t)表示t时刻的反应性；
Figure BDA0002016460590000177
表示t时刻的点堆方程缓发中子总份额；
Figure BDA0002016460590000178
表示t时刻的点堆方程d组缓发中子份额；Λ(t)表示t时刻的点堆方程中子代时间；λd表示第d组缓发中子衰变常数；
Figure BDA0002016460590000179
表示r处g能群的伴随中子标量注量率；Dg(r,t)表示r处t时刻g能群的中子扩散系数；
Figure BDA00020164605900001710
表示r处g能群t时刻的中子注量率形状函数；Σr,g(r,t)表示r处t时刻g能群的宏观中子移出截面；χg表示g能群的中子能谱；
Figure BDA00020164605900001711
表示第g能群第d缓发中子组的中子能谱；β表示缓发中子份额总数；νΣf,g'(r,t)表示r处t时刻g’能群的中子产生截面；βd表示第d组缓发中子份额；vg表示第g能群中子速度；
步骤5、利用步骤3和步骤4的获得的中子注量率形状函数和幅值函数分别在所有GPU的不同子区域上修正中子注量率在t+Δtbig的分布，并计算缓发中子先驱核在t+Δtbig的分布，进而得到t+Δtbig时刻中子注量率与缓发中子先驱核的时空分布，具体公式如下：
Figure BDA0002016460590000181
Figure BDA0002016460590000182
其中，d表示缓发中子先驱核组编号；g表示能群编号；G表示总能群数；k表示GPU设备编号；K表示总的GPU设备数；D表示总缓发中子先驱核组数；
Figure BDA0002016460590000183
表示第k个GPU子区域中r处g能群t时刻的中子注量率形状函数；[φg(r,t)]k表示第k个GPU子区域中r处t时刻g能群的中子标量注量率；n(t)表示t时刻中子注量率的幅值函数；λd表示第d组缓发中子衰变常数；[Cd(r,t)]k表示第k个GPU子区域中r处t时刻第d组缓发中子先驱核浓度；[Cd(r,t+Δtbig)]k表示第k个GPU子区域中r处t+Δtbig时刻第d组缓发中子先驱核浓度；[νΣf,g'(r,t+Δtbig)]k表示第k个GPU子区域中r处t+Δtbig时刻g’能群的中子产生截面；[φg'(r,t+Δtbig)]k表示第k个GPU子区域中r处t+Δtbig时刻g’能群的中子标量注量率；βd表示第d组缓发中子份额；
由此可以得到t+Δtbig时刻中子注量率与缓发中子先驱核的时空分布；
步骤6、重复步骤步骤1至步骤5直到达到需要计算的时空中子动力学时间。
下面以时空动力学基准问题TWIGL系列的稳态及瞬态计算结果为例说明本发明的效果：TWIGL基准问题是一个二维双群非均匀布置堆芯结构；稳态和瞬态基准题均采用D2Q4格子Boltzmann模型计算，稳态下快群及热群中子注量率与基准解的对比如图4和图5所示，瞬态归一化功率辩护额如图6所示。图7所示为GPU加速比。CPU采用Intel Core i7-3370CPU，单GPU采用NVIDIA GTX 1050Ti，双GPU计算采用GTX1050Ti GPU+GTX760双显卡。从计算效率对比可以看出，本发明在不增加实现难度的前提下极大的提高了计算效率，多个GPU计算的格子Boltzmann方法相比传统方法能实现更大的加速效果，测试方法简单可靠易于实现，测试结果准确性高，此计算结果证明本发明具有相当高的实用价值及创新性。
虽然本发明已以较佳的实施例公开如上，但其并非用以限定本发明，任何熟悉此技术的人，在不脱离本发明的精神和范围内，都可以做各种改动和修饰，因此本发明的保护范围应该以权利要求书所界定的为准。

Claims (4)

1.一种基于精细中子时空动力学格子Boltzmann方法的GPU加速方法，其特征在于，所述GPU加速方法包括：
步骤一、核反应堆稳态多群扩散方程计算，采用格子Boltzmann方法对非均匀堆芯进行扩散计算，获得稳态条件下核反应堆不同空间位置、不同能群中的中子注量率分布以及缓发中子先驱核浓度分布；
步骤二、将反应堆系统调整到临界状态，并通过临界状态相关参数计算反应堆伴随中子扩散方程，然后将所述反应堆伴随中子扩散方程进行归一化计算，获取临界初始状态下反应堆内伴随中子注量率在空间不同位置上的归一化分布情况；
步骤三、求解准静态中子时空动力学方程，获得不同时刻下中子标量注量率的时空分布情况，以及缓发中子先驱核浓度在不同时刻不同位置的分布。
2.根据权利要求1所述GPU加速方法，其特征在于，步骤一所述获得稳态条件下核反应堆不同空间位置、不同能群中的中子注量率分布以及缓发中子先驱核浓度分布的具体过程包括：
第一步、从几何卡和材料卡读入核反应堆几何参数、材料配置、不同材料的中子宏观截面参数以及动力学参数；对计算区域采用离散网格进行离散化处理；按照GPU设备数量及性能，将整体计算区域划分为多个子区域，每一个子区域由一个独立的GPU设备进行计算；在各个GPU设备中分配内存并初始化，所述计算项目获得的相关参数包括裂变源项、中子分布函数、中子注量率、缓发中子先驱核、中子宏观截面以及有效增殖因数keff；
第二步、根据第一步计算获得的相关参数和网络信息，采用OpenMP在CPU中开启多个线程同步调用所有的GPU设备并在所有GPU设备中同步进行基于格子Boltzmann方法的中子扩散计算，得到不同能群在不同空间位置的中子注量率分布，采用计算公式如下：
Figure FDA0002016460580000011
Figure FDA0002016460580000012
Figure FDA0002016460580000013
Figure FDA0002016460580000014
[φg(r)]k＝∑i[fi，g(r)]k；
通过上述计算即可获得各个GPU上不同子区域空间各个节点的中子注量率分布；其中，i表示格子Boltzmann方向；g表示能群编号；G表示总的能群数量；k表示GPU设备编号；K表示总的GPU设备数；r表示空间位置；ei表示格子Boltzmann分布函数速度；ε表示努森数，是中子平均自由程与特征尺度的比值；[fi，g(r+εei)]k表示第k个GPU子区域中r+eiε处i方向g能群内的中子分布函数；[fi，g(r)]k表示第k个GPU子区域中r处i方向g能群内的中子分布函数；[τ(r)]k表示第k个GPU子区域中r处的无量纲松弛时间因子；[fi，g(r)]k表示第k个GPU子区域中r处i方向g能群内的中子分布函数；
Figure FDA0002016460580000021
表示第k个GPU子区域中r处i方向g能群内的中子平衡态分布函数；[Ri，g(r)]k表示第k个GPU子区域中r处i方向g能群内的中子源分布函数；[Rg(r)]k表示第k个GPU子区域中r处g能群内的中子总源项；[∑g′-g(r)]k表示第k个GPU子区域中r处g’能群到g能群的宏观中子散射截面；[φg′(r)]k表示第k个GPU子区域中r处g’能群的中子标量注量率；β表示缓发中子份额总数；χg表示g能群的中子能谱；[ν∑f，g′(r)]k表示第k个GPU子区域中r处g’能群的中子产生截面；[∑r，g(r)]k表示第k个GPU子区域中r处g能群的宏观中子移出截面；d表示缓发中子先驱核组编号；
Figure FDA0002016460580000022
表示第g能群第d缓发中子组的中子能谱；λd表示第d组缓发中子衰变常数；[Cd(r)]k表示第k个GPU子区域中r处第d组缓发中子先驱核浓度；
Figure FDA0002016460580000025
表示格子Boltzmann方法中i方向的方向权重；[Dg(r)]k表示第k个GPU子区域中r处g能群的中子扩散系数；cs表示格子Boltzmann声速；keff表示反应堆堆芯的有效增殖因数；
第三步、将各个GPU内存中对应于子区域内边界节点上，指向交界面外的中子注量率分布函数拷贝回CPU，再通过CPU将所述中子注量率分布函数拷贝到原GPU的相邻子区域内边界向内分布函数中，由此完成不同子区域之间的数据交换，具体公式如下：
Figure FDA0002016460580000023
其中，g表示能群编号；G表示总的能群数量；k表示GPU设备编号；K表示总的GPU设备数；[fi，g(r)]k表示第k个GPU子区域中r处i方向g能群内的中子分布函数；
Figure FDA0002016460580000024
表示第k个GPU子区域的ei相邻子区域中r+eiε处i方向g能群内的中子分布函数；
第四步、根据第二步中获得的中子注量率分布，在各个GPU不同子区域之中通过计算获得稳态条件下的缓发中子先驱核密度分布，所采用的计算公式如下：
Figure FDA0002016460580000031
其中，g表示能群编号；G表示总的能群数量；k表示GPU设备编号；K表示总的GPU设备数；d表示缓发中子先驱核组编号；D表示缓发中子先驱核总的组数；[Cd(r)]k表示第k个GPU子区域中r处第d组缓发中子先驱核浓度；βd表示第d组缓发中子份额；λd表示第d组缓发中子衰变常数；[νΣf,g'(r)]k表示第k个GPU子区域中r处g’能群的中子产生截面；keff表示反应堆堆芯的有效增殖因数；
通过所述计算公式即可获得稳态下的缓发中子先驱核浓度在空间不同节点的分布情况，由于缓发中子先驱核和中子注量率是线性局部关系，因此不同子区域之间不需要进行缓发中子先驱核的数据交换。
3.根据权利要求1所述GPU加速方法，其特征在于，步骤二所述获取临界初始状态下反应堆内伴随中子注量率在空间不同位置上的归一化分布情况的具体过程包括：
第1步：在各个GPU内存中通过调整反应堆系统相关参数使反应堆达到临界状态，所述调整反应堆系统相关参数的具体公式如下：
Figure FDA0002016460580000032
其中，g表示能群编号；G表示总的能群数量；k表示GPU设备编号；K表示总的GPU设备数；[νΣf,g(r)]k表示第k个GPU子区域中r处g能群的中子产生截面；keff表示反应堆堆芯的有效增殖因数；
第2步：基于与步骤1相同的参数条件，在各个GPU对应子区域中计算伴随中子扩散方程；
计算所述伴随中子扩散方程的具体公式如下：
Figure FDA0002016460580000033
Figure FDA0002016460580000041
Figure FDA0002016460580000042
Figure FDA0002016460580000043
Figure FDA0002016460580000044
通过以上方程即可得到伴随中子注量率在反应堆不同空间位置的分布，其中，i表示格子Boltzmann方向；g表示能群编号；G表示总的能群数量；k表示GPU设备编号；K表示总的GPU设备数；r表示空间位置；ei表示格子Boltzmann分布函数速度；ε表示努森数，中子平均自由程与特征尺度的比值；
Figure FDA0002016460580000045
表示第k个GPU子区域中r+eiε处i方向g能群内的伴随中子分布函数；
Figure FDA0002016460580000046
表示第k个GPU子区域中r处i方向g能群内的伴随中子分布函数；[τ(r)]k表示第k个GPU子区域中r处的无量纲松弛时间因子；
Figure FDA0002016460580000047
表示第k个GPU子区域中r处i方向g能群内的伴随中子分布函数；
Figure FDA0002016460580000048
表示第k个GPU子区域中r处i方向g能群内的伴随中子平衡态分布函数；
Figure FDA0002016460580000049
表示第k个GPU子区域中r处i方向g能群内的伴随中子源分布函数；
Figure FDA00020164605800000410
表示第k个GPU子区域中r处g能群内的伴随中子总源项；[∑g-g′(r，t0)]k表示第k个GPU子区域中r处g能群t0时刻到g’能群的宏观中子散射截面；
Figure FDA00020164605800000411
表示第k个GPU子区域中r处g能群的伴随中子标量注量率；χg′表示g’能群的中子能谱；[v∑f,g(r，t0)]k表示第k个GPU子区域中r处t0时刻g能群的中子产生截面；[∑r，g(r，t0)]k表示第k个GPU子区域中r处t0时刻g能群的宏观中子移出截面；
Figure FDA00020164605800000412
表示格子Boltzmann方法方向权重；[Dg(r)]k表示第k个GPU子区域中r处g能群的中子扩散系数；cs表示格子Boltzmann声速；keff表示反应堆堆芯的有效增殖因数；通过上述参量即可表达所述伴随中子扩散方程；
第3步：将所有GPU上的中子伴随注量率拷贝回CPU，通过归一化条件将伴随中子注量率进行重构，进而获得归一化的中子伴随注量率，随后引入相应的反应性使反应堆启动；所述重构的具体公式如下：
Figure FDA0002016460580000051
其中，g表示能群编号；G表示总的能群数量；vg表示第g能群中子速度；
Figure FDA0002016460580000052
表示r处g能群的伴随中子标量注量率；
Figure FDA0002016460580000053
表示r处g能群t时刻的中子注量率形状函数。
4.根据权利要求1所述GPU加速方法，其特征在于，步骤三所述获得不同时刻下中子标量注量率的时空分布情况，以及缓发中子先驱核浓度在不同时刻不同位置的分布的具体过程包括：
步骤1、设置大时间步Δtbig、中时间步Δtmiddle以及小时间步Δtlittle；在所有不同GPU设备的不同子区域中的大时间步Δtbig尺度上预估t+Δtbig时刻的中子注量率并计算Δtbig时间步内的中子注量率形状函数，具体计算公式如下：
Figure FDA0002016460580000054
Figure FDA0002016460580000055
Figure FDA0002016460580000056
Figure FDA0002016460580000057
[φg(r,t)]k＝∑i[fi,g(r,t)]k；
Figure FDA0002016460580000058
其中，i表示格子Boltzmann方向；g表示能群编号；G表示总的能群数量；k表示GPU设备编号；K表示总的GPU设备数；D表示总缓发中子先驱核组数；r表示空间位置；ei表示格子Boltzmann分布函数速度；Δtbig表示大时间步长；vg表示第g能群中子速度；[fi，g(r+eiΔtbig，t+Δtbig)]k表示第k个GPU子区域中r+eiΔtbig处t+Δtbig时刻i方向g能群内的中子分布函数；[fi，g(r，t)]k表示第k个GPU子区域中r处t时刻i方向g能群内的中子分布函数；[τ(r，t)]k表示第k个GPU子区域中r处t时刻的无量纲松弛时间因子；[fi，g(r，t)]k表示第k个GPU子区域中r处t时刻i方向g能群内的中子分布函数；
Figure FDA0002016460580000061
表示第k个GPU子区域中r处t时刻i方向g能群内的中子平衡态分布函数；[Ri，g(r，t)]k表示第k个GPU子区域中r处t时刻i方向g能群内的中子源分布函数；[Rg(r，t)]k表示第k个GPU子区域中r处t时刻g能群内的中子总源项；[∑g′-g(r，t)]k表示第k个GPU子区域中r处t时刻g’能群到g能群的宏观中子散射截面；[φg′(r，t)]k表示第k个GPU子区域中r处t时刻g’能群的中子标量注量率；β表示缓发中子份额总数；βd表示第d组缓发中子份额；χg表示g能群的中子能谱；[ν∑f，g′(r，t)]k表示第k个GPU子区域中r处t时刻g’能群的中子产生截面；[∑r，g(r，t)]k表示第k个GPU子区域中r处t时刻g能群的宏观中子移出截面；d表示缓发中子先驱核组编号；
Figure FDA0002016460580000062
表示第g能群第d缓发中子组的中子能谱；λd表示第d组缓发中子衰变常数；[Cd(r，t)]k表示第k个GPU子区域中r处t时刻第d组缓发中子先驱核浓度；[Cd(r，t+Δtbig)]k表示第k个GPU子区域中r处t+Δtbig时刻第d组缓发中子先驱核浓度；[v∑f，g′(r，t+Δtbig)]k表示第k个GPU子区域中r处t+Δtbig时刻g’能群的中子产生截面；[φg′(r，t+Δtbig)]k表示第k个GPU子区域中r处t+Δtbig时刻g’能群的中子标量注量率；
Figure FDA0002016460580000065
表示格子Boltzmann方法方向权重；[Dg(r，t)]k表示第k个GPU子区域中r处t时刻g能群的中子扩散系数；cs表示格子Boltzmann声速；
步骤2、将不同GPU上每一个子区域内边界节点上位于t+Δtbig时刻的向外中子注量率分布函数拷贝回CPU，再将其拷贝回相邻GPU子区域内边界节点上向内中子注量率分布函数中，完成不同子区域之间的数据交换，所述数据交换过程应用的具体公式如下：
Figure FDA0002016460580000063
其中，g表示能群编号；G表示总的能群数量；k表示GPU设备编号；K表示总的GPU设备数；[fi，g(r)]k表示第k个GPU子区域中r处i方向g能群内的中子分布函数；
Figure FDA0002016460580000064
表示第k个GPU子区域的ei相邻子区域中r+eiε处i方向g能群内的中子分布函数；
步骤3、在所有GPU对应的子区域中计算中子注量率在t～t+Δtbig时间段内的形状函数，所述形状函数的计算公式如下：
Figure FDA0002016460580000071
其中，g表示能群编号；G表示总的能群数量；k表示GPU设备编号；K表示总的GPU设备数；
Figure FDA0002016460580000072
表示第k个GPU子区域中r处g能群t时刻的中子注量率形状函数；vg表示第g能群中子速度；
Figure FDA0002016460580000073
表示第k个GPU子区域中r处g能群的伴随中子标量注量率；[φg(r)]k表示第k个GPU子区域中r处g能群的中子标量注量率；
步骤4、将所有GPU设备上不同子区域的变量拷贝回CPU，在CPU上进行点堆方程计算；在每一个中时间步Δtmiddle上线性插值中子注量率形状函数；在每一个中时间步Δtmiddle内，以小时间步Δtlittle计算点堆方程，具体公式如下：
Figure FDA0002016460580000074
Figure FDA0002016460580000075
Figure FDA0002016460580000076
Figure FDA0002016460580000077
Figure FDA0002016460580000078
Figure FDA0002016460580000081
其中，n(t)表示t时刻中子注量率的幅值函数；d表示缓发中子先驱核组编号；g表示能群编号；G表示总能群数；k表示GPU设备编号；K表示总的GPU设备数；D表示总缓发中子先驱核组数；md(t)表示t时刻第d组缓发中子先驱核的幅值函数；ρ(t)表示t时刻的反应性；
Figure FDA0002016460580000082
表示t时刻的点堆方程缓发中子总份额；
Figure FDA0002016460580000083
表示t时刻的点堆方程d组缓发中子份额；Λ(t)表示t时刻的点堆方程中子代时间；λd表示第d组缓发中子衰变常数；
Figure FDA0002016460580000084
表示r处g能群的伴随中子标量注量率；Dg(r,t)表示r处t时刻g能群的中子扩散系数；
Figure FDA0002016460580000085
表示r处g能群t时刻的中子注量率形状函数；Σr,g(r,t)表示r处t时刻g能群的宏观中子移出截面；χg表示g能群的中子能谱；
Figure FDA0002016460580000086
表示第g能群第d缓发中子组的中子能谱；β表示缓发中子份额总数；νΣf,g'(r,t)表示r处t时刻g’能群的中子产生截面；βd表示第d组缓发中子份额；vg表示第g能群中子速度；
步骤5、利用步骤3和步骤4的获得的中子注量率形状函数和幅值函数分别在所有GPU的不同子区域上修正中子注量率在t+Δtbig的分布，并计算缓发中子先驱核在t+Δtbig的分布，进而得到t+Δtbig时刻中子注量率与缓发中子先驱核的时空分布，具体公式如下：
Figure FDA0002016460580000087
Figure FDA0002016460580000088
其中，d表示缓发中子先驱核组编号；g表示能群编号；G表示总能群数；k表示GPU设备编号；K表示总的GPU设备数；D表示总缓发中子先驱核组数；
Figure FDA0002016460580000089
表示第k个GPU子区域中r处g能群t时刻的中子注量率形状函数；[φg(r,t)]k表示第k个GPU子区域中r处t时刻g能群的中子标量注量率；n(t)表示t时刻中子注量率的幅值函数；λd表示第d组缓发中子衰变常数；[Cd(r,t)]k表示第k个GPU子区域中r处t时刻第d组缓发中子先驱核浓度；[Cd(r,t+Δtbig)]k表示第k个GPU子区域中r处t+Δtbig时刻第d组缓发中子先驱核浓度；[νΣf,g'(r,t+Δtbig)]k表示第k个GPU子区域中r处t+Δtbig时刻g’能群的中子产生截面；[φg'(r,t+Δtbig)]k表示第k个GPU子区域中r处t+Δtbig时刻g’能群的中子标量注量率；βd表示第d组缓发中子份额；
步骤6、重复步骤步骤1至步骤5直到达到需要计算的时空中子动力学时间。
CN201910264849.2A 2019-04-03 2019-04-03 一种基于精细中子时空动力学格子Boltzmann方法的GPU加速方法 Active CN111782384B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201910264849.2A CN111782384B (zh) 2019-04-03 2019-04-03 一种基于精细中子时空动力学格子Boltzmann方法的GPU加速方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201910264849.2A CN111782384B (zh) 2019-04-03 2019-04-03 一种基于精细中子时空动力学格子Boltzmann方法的GPU加速方法

Publications (2)

Publication Number Publication Date
CN111782384A true CN111782384A (zh) 2020-10-16
CN111782384B CN111782384B (zh) 2022-08-19

Family

ID=72754716

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201910264849.2A Active CN111782384B (zh) 2019-04-03 2019-04-03 一种基于精细中子时空动力学格子Boltzmann方法的GPU加速方法

Country Status (1)

Country Link
CN (1) CN111782384B (zh)

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114707394A (zh) * 2022-06-07 2022-07-05 西安交通大学 固定网格下考虑变形效应的反应堆瞬态中子通量模拟方法
CN115329250A (zh) * 2022-10-13 2022-11-11 中国空气动力研究与发展中心计算空气动力研究所 基于dg处理数据的方法、装置、设备及可读存储介质

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102411658A (zh) * 2011-11-25 2012-04-11 中国人民解放军国防科学技术大学 一种基于cpu和gpu协作的分子动力学加速方法
CN102681972A (zh) * 2012-04-28 2012-09-19 浪潮电子信息产业股份有限公司 一种利用GPU加速格子-Boltzmann的方法
US20130028364A1 (en) * 2010-03-29 2013-01-31 Jacobs E&C Limited Accelerator-Driven Nuclear System with Control of Effective Neutron Multiplication Coefficent
CN104318063A (zh) * 2014-09-28 2015-01-28 南华大学 基于qs/mc方法的ads次临界嬗变装置中子时空动力学研究方法
CN107145731A (zh) * 2017-04-27 2017-09-08 西安交通大学 一种高保真时空中子动力学计算的加速方法

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20130028364A1 (en) * 2010-03-29 2013-01-31 Jacobs E&C Limited Accelerator-Driven Nuclear System with Control of Effective Neutron Multiplication Coefficent
CN102411658A (zh) * 2011-11-25 2012-04-11 中国人民解放军国防科学技术大学 一种基于cpu和gpu协作的分子动力学加速方法
CN102681972A (zh) * 2012-04-28 2012-09-19 浪潮电子信息产业股份有限公司 一种利用GPU加速格子-Boltzmann的方法
CN104318063A (zh) * 2014-09-28 2015-01-28 南华大学 基于qs/mc方法的ads次临界嬗变装置中子时空动力学研究方法
CN107145731A (zh) * 2017-04-27 2017-09-08 西安交通大学 一种高保真时空中子动力学计算的加速方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
NORIHISA FUJITA ET AL: "Nuclear Fusion Simulation Code Optimization on GPU Clusters", 《2013 INTERNATIONAL CONFERENCE ON PARALLEL AND DISTRIBUTED SYSTEMS》 *
马宇等: "GPU加速的中子输运稳态格子Boltzmann方法", 《核动力工程》 *

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114707394A (zh) * 2022-06-07 2022-07-05 西安交通大学 固定网格下考虑变形效应的反应堆瞬态中子通量模拟方法
CN114707394B (zh) * 2022-06-07 2022-08-26 西安交通大学 固定网格下考虑变形效应的反应堆瞬态中子通量模拟方法
CN115329250A (zh) * 2022-10-13 2022-11-11 中国空气动力研究与发展中心计算空气动力研究所 基于dg处理数据的方法、装置、设备及可读存储介质
CN115329250B (zh) * 2022-10-13 2023-03-10 中国空气动力研究与发展中心计算空气动力研究所 基于dg处理数据的方法、装置、设备及可读存储介质

Also Published As

Publication number Publication date
CN111782384B (zh) 2022-08-19

Similar Documents

Publication Publication Date Title
Wang et al. RMC–A Monte Carlo code for reactor core analysis
Tramm et al. Memory bottlenecks and memory contention in multi-core Monte Carlo transport codes
Ren et al. Sparse LU factorization for parallel circuit simulation on GPU
CN111782384B (zh) 一种基于精细中子时空动力学格子Boltzmann方法的GPU加速方法
Khoshahval et al. A new hybrid method for multi-objective fuel management optimization using parallel PSO-SA
Li et al. OpenKMC: a KMC design for hundred-billion-atom simulation using millions of cores on Sunway Taihulight
Li et al. Simulation of the full-core pin-model by JMCT Monte Carlo neutron-photon transport code
García et al. A Collision-based Domain Decomposition scheme for large-scale depletion with the Serpent 2 Monte Carlo code
Wang et al. GPU accelerated lattice Boltzmann method in neutron kinetics problems II: Neutron transport calculation
Akbari et al. An investigation for an optimized neutron energy-group structure in thermal lattices using Particle Swarm Optimization
Zhang et al. Development of a GPU-based three-dimensional neutron transport code
Jain et al. DIAMOND–A Method of Characteristics neutron transport code using unstructured meshing
Guo et al. Neutronics/Thermal-hydraulics coupling with RMC and CTF for BEAVRS benchmark calculation
Du et al. Parallel 3D finite difference time domain simulations on graphics processors with CUDA
Legrady et al. Full core pin-level VVER-440 simulation of a rod drop experiment with the GPU-based Monte Carlo code GUARDYAN
Li et al. Improved subgroup method coupled with particle swarm optimization algorithm for intra-pellet non-uniform temperature distribution problem
Zhang et al. Heterogeneous programming and optimization of gyrokinetic toroidal code using directives
Gajurel et al. GPU acceleration of sparse neural networks
Pinheiro et al. GPU-based parallel computation in real-time modeling of atmospheric radionuclide dispersion
Du et al. Parallel 3D finite-difference time-domain method on multi-GPU systems
Li et al. Comparative analysis of coupling schemes of Monte Carlo burnup calculation in RMC
Li et al. Investigation of the efficiency optimization for the improved subgroup resonance self-shielding treatment on the GPU platform
Liang et al. Virtual lattice method for efficient Monte Carlo transport simulation of dispersion nuclear fuels
Shafii et al. Neutron flux distribution of slab reactor core using one-dimensional multi-group diffusion equation in the thermal energy region
Wang et al. Effect of Double Heterogeneity Treatment on Neutronics Modeling of HTGR Unit Cell

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant