基于CAD鄰居列表和包圍盒的蒙特卡羅粒子輸運幾何跟蹤加速方法研究

陳珍平,宋 婧,吳 斌,郝麗娟,胡麗琴,孫光耀

(1.中國科學技術大學,安徽 合肥230027;2.中國科學院核能安全技術研究所,中國科學院中子輸運理論與輻射安全重點實驗室,安徽 合肥230031)

基于CAD鄰居列表和包圍盒的蒙特卡羅粒子輸運幾何跟蹤加速方法研究

陳珍平1,2,宋 婧2,吳 斌2,郝麗娟2,胡麗琴1,2,孫光耀2

(1.中國科學技術大學,安徽 合肥230027;2.中國科學院核能安全技術研究所,中國科學院中子輸運理論與輻射安全重點實驗室,安徽 合肥230031)

幾何跟蹤主要進行蒙特卡羅粒子輸運計算中粒子位置和徑跡長度的計算,它是蒙特卡羅粒子輸運計算的關鍵技術之一。由于聚變堆幾何結構極其復雜,使得幾何跟蹤在整個蒙特卡羅粒子輸運計算中占據30%～80%的計算時間,因此幾何跟蹤方法的效率是決定聚變堆蒙特卡羅粒子輸運計算效率的重要因素之一。本文提出了基于CAD的鄰居列表和包圍盒加速方法,并基于FDS團隊自主研發的超級蒙特卡羅核計算仿真軟件系統Super MC進行實現。利用聚變堆FDS-Ⅱ和ITER模型對本文方法進行了數值驗證,測試結果表明本文方法不影響計算結果,并能使程序計算效率提高50%～60%,證明了本文方法的正確性和有效性。

蒙特卡羅;粒子輸運;幾何跟蹤;加速方法;聚變堆

蒙特卡羅粒子輸運方法由于其模擬精度高、幾何適應性強等優勢而被廣泛應用于反應堆物理、輻射屏蔽、核探測和醫學物理等領域。然而,蒙特卡羅粒子輸運方法的最大缺點在于,其收斂速度慢,計算耗時,在模擬全堆芯大規模復雜幾何模型[1]時問題尤其顯著。據統計,蒙特卡羅粒子輸運計算約30%～80%的計算時間用于粒子輸運時的幾何跟蹤[2],所以幾何跟蹤方法的效率是決定蒙特卡羅粒子輸運計算效率的重要因素之一。因此,為了提高蒙特卡羅粒子輸運計算效率,幾何跟蹤加速方法研究將具有重要意義。

在蒙特卡羅粒子輸運計算中,幾何跟蹤主要完成兩大功能,根據粒子的空間坐標和運動方向:1)查找粒子所在幾何體;2)計算粒子沿給定方向到幾何邊界的最短距離。由于粒子輸運計算時會頻繁大量使用幾何跟蹤,因此在處理大規模復雜幾何模型(如聚變堆)時,蒙特卡羅粒子輸運計算效率一定程度受到幾何跟蹤效率的限制。本文基于FDS團隊自主研發的超級蒙特卡羅核計算仿真軟件系統Super MC[3-8],在已有(Constructive Solid/Surface Geometry)CSG幾何描述框架基礎上,發展基于CAD的鄰居列表和包圍盒加速方法,以提高蒙特卡羅粒子輸運計算效率。

1 鄰居列表加速方法

在蒙特卡羅粒子輸運計算幾何跟蹤中,當模擬粒子從當前幾何體穿出時,為了查找粒子下次再進入的幾何體,傳統查找方法是一個一個地遍歷查找整個模型中的所有幾何體,直到找到粒子進入的幾何體為止,其算法復雜度為O(N)。由此可知,當處理大規模復雜幾何模型,即N值很大時,該方法的查找效率非常低,這嚴重制約了蒙特卡羅粒子輸運的計算效率。

為了解決傳統查找方法帶來的效率制約問題,鄰居列表方法[9]被提出并被廣泛應用,以提高幾何跟蹤時幾何體的查找效率。鄰居列表方法的加速原理在于,通過減少每次查找時遍歷的候選幾何體數目,即通過降低N值來達到加速的目的。傳統鄰居列表方法的基本思想:1)給模型中的每一個幾何體創建一個鄰居列表,初始列表為空;2)當粒子從當前幾何體穿出時,先查找當前幾何體的鄰居列表,若查找成功則直接返回查找到的幾何體,否則進行第3)步;3)逐個遍歷查找整個模型中的所有幾何體,直到找到為止,并把查找成功的幾何體添加到上一步穿出幾何體的鄰居列表中。在粒子輸運計算時,隨著模擬粒子數目增多,經過第2)和第3)步反復迭代,最終整個模型中的每個幾何體都會構建一個完善的鄰居列表。

然而,傳統鄰居列表方法是一種依賴模擬粒子歷史的方法,其最大缺點在于一個幾何體的鄰居列表需要模擬一定的粒子數后才能構建完成,而這個構建過程所耗費的時間隨著模型復雜度的不同而不同。更重要的是,對于一些復雜模型,這樣鄰居列表的構建過程是非常困難的。因此,傳統鄰居列表方法需要耗費一定計算時間來為每個幾何體構建鄰居列表,而在鄰居列表構建完成之前是不能完全發揮其加速效果的,從而降低了鄰居列表方法的加速效果。

因此,本文提出了一種基于CAD的鄰居列表方法來解決傳統方法中的缺點與不足。在進行蒙特卡羅粒子輸運計算之前,通過分析CAD模型幾何體之間的拓撲信息來方便地對復雜模型中的每個幾何體構建鄰居列表?；贑AD預先構建鄰居列表,一方面簡化了鄰居列表的構建過程,增強了鄰居列表方法對于復雜幾何模型的適應性;另一方面粒子輸運計算在一開始幾何跟蹤時便能利用預先生成的鄰居列表信息,充分發揮鄰居列表的優勢,提高幾何跟蹤效率。

基于CAD的鄰居列表方法如圖1所示。首先,對于任意一個給定CAD模型,自動建模模塊對模型經過簡化和修復之后,會將模型分解成凸實體。其次,將實體模型轉化為面描述模型,即用半空間面來描述凸實體的每個表面。再次,判斷幾何體間的鄰居關系,若兩個幾何體間沒有共面情況,則為非鄰居關系;反之,再判斷兩個幾何體的包圍盒是否相交,若不相交則為非鄰居關系,否則為鄰居關系。最后,根據幾何體間的鄰居關系為各幾何體構建完善的鄰居列表,用于蒙特卡羅粒子輸運計算時加速幾何體查找,提高幾何跟蹤效率。

圖1 基于CAD的鄰居列表方法Fig.1 The scheme of CAD-based neighbor search method

2 包圍盒加速方法

計算機圖形學中,包圍盒方法[10]被廣泛運用于加速射線跟蹤和碰撞檢測等領域。其加速原理在于,通過給場景中的每個物體構建一個包圍盒,在進行射線跟蹤或碰撞檢測時,首先進行射線與包圍盒干涉檢測,若射線與物體的包圍盒不發生干涉,則無須進行射線與真實物體間的檢測,從而避免因射線與真實物體進行不必要干涉檢測而帶來的額外計算。進行蒙特卡羅粒子輸運計算幾何跟蹤時,計算粒子沿給定方向到幾何邊界的最短距離這一操作同射線跟蹤具有相似性,因此可以將包圍盒方法應用于幾何跟蹤,提升蒙特卡羅粒子輸運計算效率。

進行蒙特卡羅粒子輸運模擬時,當粒子穿越幾何邊界時需要特殊處理,因此粒子每模擬一步都需要計算粒子沿其運動方向到幾何邊界的最短距離,以判斷粒子是否到達幾何邊界。傳統方法是依次求解粒子軌跡與模型中所有幾何體聯合建立的數值方程,并從所有有效解中選取最小解作為最短距離。然而,事實上模型中大部分幾何體與粒子軌跡方程并沒有有效解,所以沒有必要對這部分幾何體與粒子軌跡方程進行求解。因此,本文利用包圍盒方法以有效地將這部分幾何體加以篩選排除,減少不必要的額外計算,以加速幾何跟蹤的效率。

然而,對傳統蒙特卡羅粒子輸運計算軟件而言,其最大難點在于無法為模型中的幾何體構建包圍盒,因此本文提出了基于CAD的包圍盒方法。圖2所示為包圍盒加速方法示意圖,整個模型E中填充了A、B、C、D四個幾何體。在進行蒙特卡羅粒子輸運計算之前,通過分析模型中幾何體的結構和尺寸信息,基于CAD預先生成包圍盒數據,然后將數據導入Super MC中對模型中的幾何體A、B、C、D構建包圍盒。粒子輸運計算時,假設粒子當前位置為P,飛行軌跡(射線)為L,需要計算粒子沿軌跡L到幾何邊界的最短距離。在不使用包圍盒加速方法情況下,為了計算粒子沿飛行方向到幾何邊界的最短距離,要依次求解射線L與幾何體A、B、C、D的聯合方程組,然后從有效解中選取最小解作為最短距離。顯而易見,射線L與幾何體A、D的聯合方程組無有效解,即射線L與幾何體A、D無交點,所以這兩組無效解的求解是完全不必要的,增加了額外計算時間,降低了幾何跟蹤的效率。為了解決這類問題,減少不必要的額外計算,在使用包圍盒加速方法時,首先進行射線L與幾何體包圍盒的求交計算,若射線L與某一幾何體的包圍盒不相交,則不必進行射線L與該幾何體的聯合方程組求解,從而節省不必要的計算開銷。以圖2為例,射線L與幾何體A、D的包圍盒不相交,而只與幾何體B、C的包圍盒相交,因此只需求解射線L與幾何體B、C的聯合方程組,再從計算得到的有效解中選取最小解作為最短距離,最終得到即為粒子沿軌跡L到幾何體C邊界的最短距離。

圖2 包圍盒加速方法示意圖Fig.2 An illustration of bounding box method

3 數值驗證

基于FDS團隊自主研發的超級蒙特卡羅核計算仿真軟件系統Super MC,本文實現了基于CAD的鄰居列表和包圍盒加速方法的快速幾何跟蹤功能,并利用聚變堆FDS-Ⅱ[11-12]和ITER[13-15]模型對兩種加速方法進行了數值驗證。

3.1 聚變動力堆—FDS-Ⅱ模型

FDS-Ⅱ模型是FDS團隊發展的一個先進的聚變動力堆概念設計模型,用于驗證和校核聚變能在未來商業應用和發展中的可行性和有效性。其中,FDS-Ⅱ模型的環向主半徑達6 m,軸徑比例為3。整個模型由超過100個布爾幾何體和超過10 000個CSG基本幾何體組成,其幾何結構具有一定的復雜性。模型的具體幾何結構尺寸和材料數據參見文獻[11-12]。進行蒙特卡羅粒子輸運模擬時,其粒子源為一個內半徑4 m、外半徑8 m和高度7.6 m的圓柱殼中子源,中子能量為14.06 MeV,出射方向為各向同性分布。FDS-Ⅱ模型的具體幾何結構如圖3所示。

圖3 聚變堆FDS-Ⅱ模型Fig.3 The nuclear fusion reactor FDS-Ⅱmodel

為了驗證本文幾何跟蹤加速方法對蒙特卡羅粒子輸運計算效率的影響,對FDS-Ⅱ模型采用200 000源粒子進行固定源模擬,并每模擬完成10 000源粒子實時統計當前總計算時間,測試結果如圖4所示。由圖可知,每條測試曲線的斜率(ΔT/ΔN)表征了相應加速方法的計算效率大小,即模擬單位源粒子數所花的計算時間。因此,測試曲線斜率值越小,模擬單位源粒子數所花的計算時間也越少,那么表明相應幾何跟蹤加速方法的效率越高。由圖4可知,針對FDS-Ⅱ模型,鄰居列表方法相比包圍盒方法具有更好的加速效果,這可能是由于對FDS-Ⅱ 模型進行蒙特卡羅粒子輸運計算時,其幾何跟蹤大部分時間主要耗費在查找粒子所在幾何體功能上,此時采用更具針對性的鄰居列表加速方法能夠取得更好的效率提升。對FDS-Ⅱ模型以模擬200 000源粒子所耗費總計算時間為例(見表1),同時采用包圍盒和鄰居列表加速方法能夠獲得約55%的效率提升,其加速效果比較明顯。

圖4 FDS-Ⅱ模型蒙特卡羅粒子輸運計算效率測試Fig.4 Monte Carlo simulation performance test for FDS-Ⅱmodel

表1 FDS-Ⅱ模型20萬粒子蒙特卡羅輸運計算時間對比Table 1 Comparison of Monte Carlo simulation runtime for FDS-Ⅱmodel with 200 000 particle histories

3.2 國際熱核聚變實驗堆—ITER基準模型

國際熱核聚變反應實驗堆模型[13-15]是由國際ITER機構發布,主要用于驗證和校核各種粒子輸運計算程序的計算功能和性能。ITER基準模型如圖5所示,它主要包層、偏濾器、第一壁、CS/TF/PF線圈、上下窗口以及生物屏蔽層等結構組成。整個模型幾何結構極其復雜,由超過3 100個布爾幾何體和22 000個CSG基本幾何體構成。

圖5 國際熱核聚變實驗堆ITER基準模型Fig.5 The nuclear fusion reactor ITER benchmark model

為了進一步驗證本文方法對于處理大尺度復雜幾何模型的加速效果,采用ITER模型對其計算效率進行了測試分析,其測試條件與FDS-Ⅱ測試一致,計算結果見圖6和表2。顯而易見,從圖6可知,對ITER模型進行蒙特卡羅粒子輸運計算時,包圍盒方法比鄰居列表方法具有更好的加速效果,這主要是因為針對幾何結構極其復雜模型,其幾何跟蹤時間主要花在計算粒子沿其運動方向到幾何邊界的最短距離上。而包圍盒方法的加速原理就在于,在計算最短距離時,通過包圍盒與粒子軌跡方程的求交計算,篩選出與粒子軌跡不相交的幾何體,從而減少不必要的額外復雜計算。如表2所示,使用包圍盒加速方法能夠使Super MC的計算效率提升54%。若同時使用包圍盒方法和鄰居列表方法,可以使Super MC的計算效率提高約62%,進一步驗證了本文加速方法對處理復雜幾何模型在效率提升方面的可行性和有效性。

圖6 ITER模型蒙特卡羅粒子輸運計算效率測試Fig.6 Monte Carlo simulation performance test for ITER model

表2 ITER模型20萬粒子蒙特卡羅輸運計算時間對比Table 2 Comparison of Monte Carlo simulation runtime for ITER model with 200 000 particle histories

總結

幾何跟蹤是影響蒙特卡羅粒子輸運計算效率的重要因素之一。針對這個問題,本文發展了基于CAD的鄰居列表和包圍盒方法對幾何跟蹤進行加速處理。采用聚變堆FDS-Ⅱ和ITER基準模型對本文加速方法進行了測試,結果表明本文方法不影響計算結果,并能使蒙特卡羅粒子輸運計算效率提升50%～60%,證明了本文方法的正確性、可行性和有效性。

致謝

本文開展研究工作中,得到了FDS團隊其他成員的大力幫助和支持,在此深表感謝!

[1] Kord Smith,Benoit Forget.Challenges in the Development of High-fidelity LWR Core Neutronics Tools[C]．International Conference on Mathematics and Computational Methods Applied to Nuclear Science and Engineering,Idaho,USA,May 5-9,2013.

[2] Brian R Nease,et al.Geometric Templates for Improved Tracking Performance in Monte Carlo Codes[C]．Joint International Conference on Supercomputing in Nuclear Applications and Monte Carlo,Paris,France,October 27-31,2013.

[3] Wu Y.,et al.CAD-Based Monte Carlo Program for Integrated Simulation of Nuclear System Super MC[J].Annals of Nuclear Energy,2015,82:162－168.

[4] Wu Y.,FDS Team.CAD-Based Interface Programs for Fusion Neutron Transport Simulation[J]．Fusion Engineering and Design,2009,84(7-11):1987-1992.

[5] Wu Y.,et al.A Discrete Ordinates Nodal Method for One-dimensional Neutron Transport Calculation in Curvilinear Geometries[J]．Nuclear Science and Engineering,1999,133(3):350-357.

[6] 吳宜燦等.大型集成多功能中子學計算與分析系統VisualBUS的研究與發展[J]．核科學與工程,2007,27(4):365-373.

[7] Li Y,et al.Benchmarking of MCAM4.0 with the ITER 3D Model[J]．Fusion Engineering and Design,2007,82:2861-2866.

[8] Hu H.,et al.Benchmarking of SNAM with the ITER 3D Model[J]．Fusion Engineering and Design,2007,82:2867-2871.

[9] Forest B.Brown.Fundamentals of Monte Carlo Transport[R]．Los Alamos National Laboratory,USA,2004.

[10] Wald,I.,et al.Ray Tracing Deformable Scenes Using Dynamic Bounding Volume Hierarchies[J]．ACM Transactions on Graphics,26,2007.

[11] Wu Y.,FDS Team.Conceptual Design of the China Fusion Power Plant FDS-Ⅱ [J]．Fusion Engineering and Design,2008,83(10-12):1683-1689.

[12] Wu Y.,FDS Team.Conceptual Design Activities of FDS Series Fusion Power Plants in China[J]．Fusion Engineering and Design,2006,81(23-24):2713-2718.

[13] Wu Y.,FDS Team.Conceptual Design and Testing Strategy of a Dual Functional Lithium-Lead Test Blanket Module in ITER and EAST[J]．Nuclear Fusion,2007,47(11):1533-1539.

[14] Wu Y.,FDS Team.Design Analysis of the China Dual-Functional Lithium Lead(DFLL)Test Blanket Module in ITER[J]．Fusion Engineering and Design,2007,82:1893-1903.

[15] Song J.,Sun G.,Chen Z.,et al.Benchmarking of CAD-based Super MC with ITER benchmark model[J],Fusion Engineering and Design,2014,89(11),2499-2503.

CAD-based neighbor search and bounding box algorithms for geometry navigation acceleration in Monte Carlo particle transport simulation

CHEN Zhen-ping1,2,SONG Jing2,WU Bin2,HAO Li-juan2,HU Li-qin1,2,SUN Guang-yao2
(1.University of Science and Technology of China,Hefei of Anhui Prov.230027,China;2.Key Laboratory of Neutronics and Radiation Safety,Institute of Nuclear Energy Safety Technology,Chinese Academy of Sciences,Hefei of Anhui Prov.230031,China)

Geometry navigation plays the most fundamental role in Monte Carlo particle transport simulation.It's mainly responsible for locating a particle inside which geometry volume it is and computing the distance to the volume boundary along the certain particle trajectory during each particle history.Geometry navigation directlyaffects the run-time performance of the Monte Carlo particle transport simulation,especially for complicated fusion reactor models.Thus,two CAD-based geometry acceleration algorithms,the neighbor search and the bounding box,are presented for improving geometry navigation performance.The algorithms have been implemented in the Super Monte Carlo Simulation Program for Nuclear and Radiation Process(Super MC).The fusion reactors of FDS-Ⅱand ITER benchmark models have been tested to highlight the efficiency gains that can be achieved by using the acceleration algorithms.Testing results showed that efficiency of Monte Carlo simulation can be considerably enhanced by 50%to 60%with the acceleration algorithms.

Monte Carlo;particle transport;Geometry navigation;Acceleration algorithm;Fusion reactor

TL329＋.2

A

0258-0918(2016)02-0212-06

2016-02-21

本工作得到中國科學院戰略性先導科技專項(XDA03040000),國家自然科學基金(91026004,11305205,11305203,11405204),國家ITER 973計劃(2014GB112001)和安徽省自然科學基金(1308085QH138)的資助。

陳珍平(1987—),男,重慶人,博士,研究生,主要從事蒙特卡羅粒子輸運相關研究工作

孫光耀:guangyao.sun＠fds．org．cn