CN111782384A - 一种基于精细中子时空动力学格子Boltzmann方法的GPU加速方法 - Google Patents
一种基于精细中子时空动力学格子Boltzmann方法的GPU加速方法 Download PDFInfo
- 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
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)
- Theoretical Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Mathematical Physics (AREA)
- Software Systems (AREA)
- Mathematical Analysis (AREA)
- Pure & Applied Mathematics (AREA)
- Data Mining & Analysis (AREA)
- Mathematical Optimization (AREA)
- General Engineering & Computer Science (AREA)
- Computational Mathematics (AREA)
- Operations Research (AREA)
- Algebra (AREA)
- Databases & Information Systems (AREA)
- Monitoring And Testing Of Nuclear Reactors (AREA)
Abstract
本发明提出了一种基于精细中子时空动力学格子Boltzmann方法的GPU加速方法,属于核反应堆设计及安全分析技术领域。所述方法为:在多核GPU设备上采用格子Boltzmann方法对非均匀堆芯进行扩散计算得到稳态条件下核反应堆不同空间位置、不同能群中的中子注量率分布以及缓发中子先驱核浓度分布;然后,采用格子Boltzmann方法求解准静态时空动力学方程,得到不同时刻下中子标量注量率的时空分布情况以及缓发中子先驱核浓度在不同时刻不同位置的分布。本发明所述加速方法能降低反应堆大规模并行计算的难度,同时大幅度提高反应堆中子时空动力学问题的计算速度;消除了传统方法实现复杂、成本高、加速效率不高的缺陷。
Description
技术领域
本发明涉及一种基于精细中子时空动力学格子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方法的中子扩散计算,得到不同能群在不同空间位置的中子注量率分布,采用计算公式如下:
[φ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能群内的中子分布函数;表示第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表示缓发中子先驱核组编号;表示第g能群第d缓发中子组的中子能谱;λd表示第d组缓发中子衰变常数;[Cd(r)]k表示第k个GPU子区域中r处第d组缓发中子先驱核浓度;表示格子Boltzmann方法中i方向的方向权重;[Dg(r)]k表示第k个GPU子区域中r处g能群的中子扩散系数;cs表示格子Boltzmann声速;keff表示反应堆堆芯的有效增殖因数;
第三步、将各个GPU内存中对应于子区域内边界节点上,指向交界面外的中子注量率分布函数拷贝回CPU,再通过CPU将所述中子注量率分布函数拷贝到原GPU的相邻子区域内边界向内分布函数中,由此完成不同子区域之间的数据交换,具体公式如下:
其中,g表示能群编号;G表示总的能群数量;k表示GPU设备编号;K表示总的GPU设备数;[fi,g(r)]k表示第k个GPU子区域中r处i方向g能群内的中子分布函数;表示第k个GPU子区域的ei相邻子区域中r+eiε处i方向g能群内的中子分布函数;
第四步、根据第二步中获得的中子注量率分布,在各个GPU不同子区域之中通过计算获得稳态条件下的缓发中子先驱核密度分布,所采用的计算公式如下:
其中,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内存中通过调整反应堆系统相关参数使反应堆达到临界状态,所述调整反应堆系统相关参数的具体公式如下:
其中,g表示能群编号;G表示总的能群数量;k表示GPU设备编号;K表示总的GPU设备数;[νΣf,g(r)]k表示第k个GPU子区域中r处g能群的中子产生截面;keff表示反应堆堆芯的有效增殖因数;
第2步:基于与步骤1相同的参数条件,在各个GPU对应子区域中计算伴随中子扩散方程;计算所述伴随中子扩散方程的具体公式如下:
通过以上方程即可得到伴随中子注量率在反应堆不同空间位置的分布,其中,i表示格子Boltzmann方向;g表示能群编号;G表示总的能群数量;k表示GPU设备编号;K表示总的GPU设备数;r表示空间位置;ei表示格子Boltzmann分布函数速度;ε表示努森数,中子平均自由程与特征尺度的比值;表示第k个GPU子区域中r+eiε处i方向g能群内的伴随中子分布函数;表示第k个GPU子区域中r处i方向g能群内的伴随中子分布函数;[τ(r)]k表示第k个GPU子区域中r处的无量纲松弛时间因子;表示第k个GPU子区域中r处i方向g能群内的伴随中子分布函数;表示第k个GPU子区域中r处i方向g能群内的伴随中子平衡态分布函数;表示第k个GPU子区域中r处i方向g能群内的伴随中子源分布函数;表示第k个GPU子区域中r处g能群内的伴随中子总源项;[Σg-g'(r,t0)]k表示第k个GPU子区域中r处g能群t0时刻到g’能群的宏观中子散射截面;表示第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能群的宏观中子移出截面;表示格子Boltzmann方法方向权重;[Dg(r)]k表示第k个GPU子区域中r处g能群的中子扩散系数;cs表示格子Boltzmann声速;keff表示反应堆堆芯的有效增殖因数;通过上述参量即可表达所述伴随中子扩散方程;
第3步:将所有GPU上的中子伴随注量率拷贝回CPU,通过归一化条件将伴随中子注量率进行重构,进而获得归一化的中子伴随注量率,随后引入相应的反应性使反应堆启动;所述重构的具体公式如下:
由此得到归一化的中子伴随注量率,随后引入相应的反应性使反应堆启动
进一步地,步骤三所述获得不同时刻下中子标量注量率的时空分布情况,以及缓发中子先驱核浓度在不同时刻不同位置的分布的具体过程包括:
步骤1、设置大时间步Δtbig、中时间步Δtmiddle以及小时间步Δtlittle;在所有不同GPU设备的不同子区域中的大时间步Δtbig尺度上预估t+Δtbig时刻的中子注量率并计算Δtbig时间步内的中子注量率形状函数,具体计算公式如下:
[φg(r,t)]k=∑i[fi,g(r,t)]k;
其中,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能群内的中子分布函数;表示第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表示缓发中子先驱核组编号;表示第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’能群的中子标量注量率;表示格子Boltzmann方法方向权重;[Dg(r,t)]k表示第k个GPU子区域中r处t时刻g能群的中子扩散系数;cs表示格子Boltzmann声速;
步骤2、将不同GPU上每一个子区域内边界节点上位于t+Δtbig时刻的向外中子注量率分布函数拷贝回CPU,再将其拷贝回相邻GPU子区域内边界节点上向内中子注量率分布函数中,完成不同子区域之间的数据交换,所述数据交换过程应用的具体公式如下:
其中,g表示能群编号;G表示总的能群数量;k表示GPU设备编号;K表示总的GPU设备数;[fi,g(r)]k表示第k个GPU子区域中r处i方向g能群内的中子分布函数;表示第k个GPU子区域的ei相邻子区域中r+eiε处i方向g能群内的中子分布函数;
步骤3、在所有GPU对应的子区域中计算中子注量率在t~t+Δtbig时间段内的形状函数,所述形状函数的计算公式如下:
其中,g表示能群编号;G表示总的能群数量;k表示GPU设备编号;K表示总的GPU设备数;表示第k个GPU子区域中r处g能群t时刻的中子注量率形状函数;vg表示第g能群中子速度;表示第k个GPU子区域中r处g能群的伴随中子标量注量率;[φg(r)]k表示第k个GPU子区域中r处g能群的中子标量注量率;
由以上步骤可以获得该大时间步内的中子注量率形状函数;
步骤4、将所有GPU设备上不同子区域的变量拷贝回CPU,在CPU上进行点堆方程计算;在每一个中时间步Δtmiddle上线性插值中子注量率形状函数;在每一个中时间步Δtmiddle内,以小时间步Δtlittle计算点堆方程,具体公式如下:
其中,n(t)表示t时刻中子注量率的幅值函数;d表示缓发中子先驱核组编号;g表示能群编号;G表示总能群数;k表示GPU设备编号;K表示总的GPU设备数;D表示总缓发中子先驱核组数;md(t)表示t时刻第d组缓发中子先驱核的幅值函数;ρ(t)表示t时刻的反应性;表示t时刻的点堆方程缓发中子总份额;表示t时刻的点堆方程d组缓发中子份额;Λ(t)表示t时刻的点堆方程中子代时间;λd表示第d组缓发中子衰变常数;表示r处g能群的伴随中子标量注量率;Dg(r,t)表示r处t时刻g能群的中子扩散系数;表示r处g能群t时刻的中子注量率形状函数;Σr,g(r,t)表示r处t时刻g能群的宏观中子移出截面;χg表示g能群的中子能谱;表示第g能群第d缓发中子组的中子能谱;β表示缓发中子份额总数;νΣf,g'(r,t)表示r处t时刻g’能群的中子产生截面;βd表示第d组缓发中子份额;vg表示第g能群中子速度;
步骤5、利用步骤3和步骤4的获得的中子注量率形状函数和幅值函数分别在所有GPU的不同子区域上修正中子注量率在t+Δtbig的分布,并计算缓发中子先驱核在t+Δtbig的分布,进而得到t+Δtbig时刻中子注量率与缓发中子先驱核的时空分布,具体公式如下:
其中,d表示缓发中子先驱核组编号;g表示能群编号;G表示总能群数;k表示GPU设备编号;K表示总的GPU设备数;D表示总缓发中子先驱核组数;表示第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方法的中子扩散计算,得到不同能群在不同空间位置的中子注量率分布,采用计算公式如下:
[φ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能群内的中子分布函数;表示第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表示缓发中子先驱核组编号;表示第g能群第d缓发中子组的中子能谱;λd表示第d组缓发中子衰变常数;[Cd(r)]k表示第k个GPU子区域中r处第d组缓发中子先驱核浓度;表示格子Boltzmann方法中i方向的方向权重;[Dg(r)]k表示第k个GPU子区域中r处g能群的中子扩散系数;cs表示格子Boltzmann声速;keff表示反应堆堆芯的有效增殖因数;
第三步、将各个GPU内存中对应于子区域内边界节点上,指向交界面外的中子注量率分布函数拷贝回CPU,再通过CPU将所述中子注量率分布函数拷贝到原GPU的相邻子区域内边界向内分布函数中,由此完成不同子区域之间的数据交换,具体公式如下:
其中,g表示能群编号;G表示总的能群数量;k表示GPU设备编号;K表示总的GPU设备数;[fi,g(r)]k表示第k个GPU子区域中r处i方向g能群内的中子分布函数;表示第k个GPU子区域的ei相邻子区域中r+eiε处i方向g能群内的中子分布函数;
第四步、根据第二步中获得的中子注量率分布,在各个GPU不同子区域之中通过计算获得稳态条件下的缓发中子先驱核密度分布,所采用的计算公式如下:
其中,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内存中通过调整反应堆系统相关参数使反应堆达到临界状态,所述调整反应堆系统相关参数的具体公式如下:
其中,g表示能群编号;G表示总的能群数量;k表示GPU设备编号;K表示总的GPU设备数;[νΣf,g(r)]k表示第k个GPU子区域中r处g能群的中子产生截面;keff表示反应堆堆芯的有效增殖因数;
第2步:基于与步骤1相同的参数条件,在各个GPU对应子区域中计算伴随中子扩散方程;计算所述伴随中子扩散方程的具体公式如下:
通过以上方程即可得到伴随中子注量率在反应堆不同空间位置的分布,其中,i表示格子Boltzmann方向;g表示能群编号;G表示总的能群数量;k表示GPU设备编号;K表示总的GPU设备数;r表示空间位置;ei表示格子Boltzmann分布函数速度;ε表示努森数,中子平均自由程与特征尺度的比值;表示第k个GPU子区域中r+eiε处i方向g能群内的伴随中子分布函数;表示第k个GPU子区域中r处i方向g能群内的伴随中子分布函数;[τ(r)]k表示第k个GPU子区域中r处的无量纲松弛时间因子;[fi,g(r)]k表示第k个GPU子区域中r处i方向g能群内的伴随中子分布函数;表示第k个GPU子区域中r处i方向g能群内的伴随中子平衡态分布函数;表示第k个GPU子区域中r处i方向g能群内的伴随中子源分布函数;表示第k个GPU子区域中r处g能群内的伴随中子总源项;[Σg-g'(r,t0)]k表示第k个GPU子区域中r处g能群t0时刻到g’能群的宏观中子散射截面;表示第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能群的宏观中子移出截面;表示格子Boltzmann方法方向权重;[Dg(r)]k表示第k个GPU子区域中r处g能群的中子扩散系数;cs表示格子Boltzmann声速;keff表示反应堆堆芯的有效增殖因数;通过上述参量即可表达所述伴随中子扩散方程;
第3步:将所有GPU上的中子伴随注量率拷贝回CPU,通过归一化条件将伴随中子注量率进行重构,进而获得归一化的中子伴随注量率,随后引入相应的反应性使反应堆启动;所述重构的具体公式如下:
其中,g表示能群编号;G表示总的能群数量;vg表示第g能群中子速度;表示r处g能群的伴随中子标量注量率;表示r处g能群t时刻的中子注量率形状函数。由此得到归一化的中子伴随注量率,随后引入相应的反应性使反应堆启动。
步骤三所述获得不同时刻下中子标量注量率的时空分布情况,以及缓发中子先驱核浓度在不同时刻不同位置的分布的具体过程包括:
步骤1、设置大时间步Δtbig、中时间步Δtmiddle以及小时间步Δtlittle;在所有不同GPU设备的不同子区域中的大时间步Δtbig尺度上预估t+Δtbig时刻的中子注量率并计算Δtbig时间步内的中子注量率形状函数,具体计算公式如下:
[φg(r,t)]k=∑i[fi,g(r,t)]k;
其中,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能群内的中子分布函数;表示第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表示缓发中子先驱核组编号;表示第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’能群的中子标量注量率;表示格子Boltzmann方法方向权重;[Dg(r,t)]k表示第k个GPU子区域中r处t时刻g能群的中子扩散系数;cs表示格子Boltzmann声速;
步骤2、将不同GPU上每一个子区域内边界节点上位于t+Δtbig时刻的向外中子注量率分布函数拷贝回CPU,再将其拷贝回相邻GPU子区域内边界节点上向内中子注量率分布函数中,完成不同子区域之间的数据交换,所述数据交换过程应用的具体公式如下:
其中,g表示能群编号;G表示总的能群数量;k表示GPU设备编号;K表示总的GPU设备数;[fi,g(r)]k表示第k个GPU子区域中r处i方向g能群内的中子分布函数;表示第k个GPU子区域的ei相邻子区域中r+eiε处i方向g能群内的中子分布函数;
步骤3、在所有GPU对应的子区域中计算中子注量率在t~t+Δtbig时间段内的形状函数,所述形状函数的计算公式如下:
其中,g表示能群编号;G表示总的能群数量;k表示GPU设备编号;K表示总的GPU设备数;表示第k个GPU子区域中r处g能群t时刻的中子注量率形状函数;vg表示第g能群中子速度;表示第k个GPU子区域中r处g能群的伴随中子标量注量率;[φg(r)]k表示第k个GPU子区域中r处g能群的中子标量注量率;
由以上步骤可以获得该大时间步内的中子注量率形状函数;
步骤4、将所有GPU设备上不同子区域的变量拷贝回CPU,在CPU上进行点堆方程计算;在每一个中时间步Δtmiddle上线性插值中子注量率形状函数;在每一个中时间步Δtmiddle内,以小时间步Δtlittle计算点堆方程,具体公式如下:
其中,n(t)表示t时刻中子注量率的幅值函数;d表示缓发中子先驱核组编号;g表示能群编号;G表示总能群数;k表示GPU设备编号;K表示总的GPU设备数;D表示总缓发中子先驱核组数;md(t)表示t时刻第d组缓发中子先驱核的幅值函数;ρ(t)表示t时刻的反应性;表示t时刻的点堆方程缓发中子总份额;表示t时刻的点堆方程d组缓发中子份额;Λ(t)表示t时刻的点堆方程中子代时间;λd表示第d组缓发中子衰变常数;表示r处g能群的伴随中子标量注量率;Dg(r,t)表示r处t时刻g能群的中子扩散系数;表示r处g能群t时刻的中子注量率形状函数;Σr,g(r,t)表示r处t时刻g能群的宏观中子移出截面;χg表示g能群的中子能谱;表示第g能群第d缓发中子组的中子能谱;β表示缓发中子份额总数;νΣf,g'(r,t)表示r处t时刻g’能群的中子产生截面;βd表示第d组缓发中子份额;vg表示第g能群中子速度;
步骤5、利用步骤3和步骤4的获得的中子注量率形状函数和幅值函数分别在所有GPU的不同子区域上修正中子注量率在t+Δtbig的分布,并计算缓发中子先驱核在t+Δtbig的分布,进而得到t+Δtbig时刻中子注量率与缓发中子先驱核的时空分布,具体公式如下:
其中,d表示缓发中子先驱核组编号;g表示能群编号;G表示总能群数;k表示GPU设备编号;K表示总的GPU设备数;D表示总缓发中子先驱核组数;表示第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方法的中子扩散计算,得到不同能群在不同空间位置的中子注量率分布,采用计算公式如下:
[φ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能群内的中子分布函数;表示第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表示缓发中子先驱核组编号;表示第g能群第d缓发中子组的中子能谱;λd表示第d组缓发中子衰变常数;[Cd(r)]k表示第k个GPU子区域中r处第d组缓发中子先驱核浓度;表示格子Boltzmann方法中i方向的方向权重;[Dg(r)]k表示第k个GPU子区域中r处g能群的中子扩散系数;cs表示格子Boltzmann声速;keff表示反应堆堆芯的有效增殖因数;
第三步、将各个GPU内存中对应于子区域内边界节点上,指向交界面外的中子注量率分布函数拷贝回CPU,再通过CPU将所述中子注量率分布函数拷贝到原GPU的相邻子区域内边界向内分布函数中,由此完成不同子区域之间的数据交换,具体公式如下:
其中,g表示能群编号;G表示总的能群数量;k表示GPU设备编号;K表示总的GPU设备数;[fi,g(r)]k表示第k个GPU子区域中r处i方向g能群内的中子分布函数;表示第k个GPU子区域的ei相邻子区域中r+eiε处i方向g能群内的中子分布函数;
第四步、根据第二步中获得的中子注量率分布,在各个GPU不同子区域之中通过计算获得稳态条件下的缓发中子先驱核密度分布,所采用的计算公式如下:
其中,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内存中通过调整反应堆系统相关参数使反应堆达到临界状态,所述调整反应堆系统相关参数的具体公式如下:
其中,g表示能群编号;G表示总的能群数量;k表示GPU设备编号;K表示总的GPU设备数;[νΣf,g(r)]k表示第k个GPU子区域中r处g能群的中子产生截面;keff表示反应堆堆芯的有效增殖因数;
第2步:基于与步骤1相同的参数条件,在各个GPU对应子区域中计算伴随中子扩散方程;
计算所述伴随中子扩散方程的具体公式如下:
通过以上方程即可得到伴随中子注量率在反应堆不同空间位置的分布,其中,i表示格子Boltzmann方向;g表示能群编号;G表示总的能群数量;k表示GPU设备编号;K表示总的GPU设备数;r表示空间位置;ei表示格子Boltzmann分布函数速度;ε表示努森数,中子平均自由程与特征尺度的比值;表示第k个GPU子区域中r+eiε处i方向g能群内的伴随中子分布函数;表示第k个GPU子区域中r处i方向g能群内的伴随中子分布函数;[τ(r)]k表示第k个GPU子区域中r处的无量纲松弛时间因子;表示第k个GPU子区域中r处i方向g能群内的伴随中子分布函数;表示第k个GPU子区域中r处i方向g能群内的伴随中子平衡态分布函数;表示第k个GPU子区域中r处i方向g能群内的伴随中子源分布函数;表示第k个GPU子区域中r处g能群内的伴随中子总源项;[∑g-g′(r,t0)]k表示第k个GPU子区域中r处g能群t0时刻到g’能群的宏观中子散射截面;表示第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能群的宏观中子移出截面;表示格子Boltzmann方法方向权重;[Dg(r)]k表示第k个GPU子区域中r处g能群的中子扩散系数;cs表示格子Boltzmann声速;keff表示反应堆堆芯的有效增殖因数;通过上述参量即可表达所述伴随中子扩散方程;
第3步:将所有GPU上的中子伴随注量率拷贝回CPU,通过归一化条件将伴随中子注量率进行重构,进而获得归一化的中子伴随注量率,随后引入相应的反应性使反应堆启动;所述重构的具体公式如下:
4.根据权利要求1所述GPU加速方法,其特征在于,步骤三所述获得不同时刻下中子标量注量率的时空分布情况,以及缓发中子先驱核浓度在不同时刻不同位置的分布的具体过程包括:
步骤1、设置大时间步Δtbig、中时间步Δtmiddle以及小时间步Δtlittle;在所有不同GPU设备的不同子区域中的大时间步Δtbig尺度上预估t+Δtbig时刻的中子注量率并计算Δtbig时间步内的中子注量率形状函数,具体计算公式如下:
[φg(r,t)]k=∑i[fi,g(r,t)]k;
其中,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能群内的中子分布函数;表示第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表示缓发中子先驱核组编号;表示第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’能群的中子标量注量率;表示格子Boltzmann方法方向权重;[Dg(r,t)]k表示第k个GPU子区域中r处t时刻g能群的中子扩散系数;cs表示格子Boltzmann声速;
步骤2、将不同GPU上每一个子区域内边界节点上位于t+Δtbig时刻的向外中子注量率分布函数拷贝回CPU,再将其拷贝回相邻GPU子区域内边界节点上向内中子注量率分布函数中,完成不同子区域之间的数据交换,所述数据交换过程应用的具体公式如下:
其中,g表示能群编号;G表示总的能群数量;k表示GPU设备编号;K表示总的GPU设备数;[fi,g(r)]k表示第k个GPU子区域中r处i方向g能群内的中子分布函数;表示第k个GPU子区域的ei相邻子区域中r+eiε处i方向g能群内的中子分布函数;
步骤3、在所有GPU对应的子区域中计算中子注量率在t~t+Δtbig时间段内的形状函数,所述形状函数的计算公式如下:
其中,g表示能群编号;G表示总的能群数量;k表示GPU设备编号;K表示总的GPU设备数;表示第k个GPU子区域中r处g能群t时刻的中子注量率形状函数;vg表示第g能群中子速度;表示第k个GPU子区域中r处g能群的伴随中子标量注量率;[φg(r)]k表示第k个GPU子区域中r处g能群的中子标量注量率;
步骤4、将所有GPU设备上不同子区域的变量拷贝回CPU,在CPU上进行点堆方程计算;在每一个中时间步Δtmiddle上线性插值中子注量率形状函数;在每一个中时间步Δtmiddle内,以小时间步Δtlittle计算点堆方程,具体公式如下:
其中,n(t)表示t时刻中子注量率的幅值函数;d表示缓发中子先驱核组编号;g表示能群编号;G表示总能群数;k表示GPU设备编号;K表示总的GPU设备数;D表示总缓发中子先驱核组数;md(t)表示t时刻第d组缓发中子先驱核的幅值函数;ρ(t)表示t时刻的反应性;表示t时刻的点堆方程缓发中子总份额;表示t时刻的点堆方程d组缓发中子份额;Λ(t)表示t时刻的点堆方程中子代时间;λd表示第d组缓发中子衰变常数;表示r处g能群的伴随中子标量注量率;Dg(r,t)表示r处t时刻g能群的中子扩散系数;表示r处g能群t时刻的中子注量率形状函数;Σr,g(r,t)表示r处t时刻g能群的宏观中子移出截面;χg表示g能群的中子能谱;表示第g能群第d缓发中子组的中子能谱;β表示缓发中子份额总数;νΣf,g'(r,t)表示r处t时刻g’能群的中子产生截面;βd表示第d组缓发中子份额;vg表示第g能群中子速度;
步骤5、利用步骤3和步骤4的获得的中子注量率形状函数和幅值函数分别在所有GPU的不同子区域上修正中子注量率在t+Δtbig的分布,并计算缓发中子先驱核在t+Δtbig的分布,进而得到t+Δtbig时刻中子注量率与缓发中子先驱核的时空分布,具体公式如下:
其中,d表示缓发中子先驱核组编号;g表示能群编号;G表示总能群数;k表示GPU设备编号;K表示总的GPU设备数;D表示总缓发中子先驱核组数;表示第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直到达到需要计算的时空中子动力学时间。
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)
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)
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 | 西安交通大学 | 一种高保真时空中子动力学计算的加速方法 |
-
2019
- 2019-04-03 CN CN201910264849.2A patent/CN111782384B/zh active Active
Patent Citations (5)
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)
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)
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 | |
Ren et al. | Sparse LU factorization for parallel circuit simulation on GPU | |
Shimokawabe et al. | 145 TFlops performance on 3990 GPUs of TSUBAME 2.0 supercomputer for an operational weather prediction | |
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 | |
Yin et al. | Multi-group effective cross section calculation method for Fully Ceramic Micro-encapsulated fuel | |
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 | |
Wu et al. | Parallel artificial neural network using CUDA-enabled GPU for extracting hydraulic domain knowledge of large water distribution systems | |
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 | |
Li et al. | Forced propagation method for Monte Carlo fission source convergence acceleration in the RMC | |
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 |
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 |