面向大規(guī)模并行的全隱式托卡馬克MHD數值模擬*

蔣子超，江俊揚，孫哲，姚清河

1. 中山大學航空航天學院，廣東廣州 510006

2. 日本理化學研究所，日本東京197-0804

等離子體的MHD 不穩(wěn)定性是限制托卡馬克等磁約束核聚變裝置長期穩(wěn)定運行的主要因素之一。撕裂模等MHD不穩(wěn)定性的發(fā)展容易導致磁面破壞、輸運增強，從而導致等離子體的約束退化。因此，探索等離子體不穩(wěn)定性的抑制手段已經成為目前磁約束核聚變研究中的關鍵問題之一［1］。

針對托卡馬克中的磁流體不穩(wěn)定性，目前學界多采用基于MHD方程的數值模擬方法進行研究。20 世紀90 年代， Barmin［2］和Nessyahu 等［3］已分別采用逆風格式和二階精度高分辨率無振蕩的中心差分格式求解一維/二維MHD 方程［4］。此后，學界開始嘗試對三維托卡馬克位形下的等離子體進行數值模擬，但是受制于有限的計算能力，針對托卡馬克等離子體的仿真求解器多基于全顯或半隱格式，且問題規(guī)模和模型精度也較為受限。Harned 和Kerner［5］以及Harned 和Schnack 等［6-7］首先提出了簡化代數算子的方法以半隱格式進行求解，近年來經過多次改進與優(yōu)化［5-7］，目前已被廣泛應用于MHD仿真計算領域求解器的開發(fā)，例如：XTOR［8］和NIMROD［9］。此外，較為成熟的托卡馬克磁流體求解器還有M3D［10］、BOUT++［11］等；我國比較成熟的托卡馬克等離子體求解器包括浙江大學馬志為教授團隊開發(fā)的CLT［12］/CLT-K［13］求解器。

雖然基于全顯式和半隱式的求解器已經取得了許多計算成果，但是基于此開發(fā)的求解算法的時間步長嚴格受限于CFL（courant-friedrichs-lewy）條件。撕裂模等MHD 不穩(wěn)定性的形成及飽和過程的模擬需要較長的時間跨度（通常大于104個Alfvén 時間），因此全顯式和半隱式將導致過多的計算步數，而采用全隱格式則可解除CFL 條件的限制，實現較大時間步長的全局求解。此外，捕捉非理想等離子體不穩(wěn)定性需要高分辨率網格，但對計算機內存以及計算機求解能力的要求就越嚴苛，一般的主機和集群無法滿足要求。

基于上述因素，本文提出了一種面向大規(guī)模并行計算的MHD 方程組的全隱式數值求解方法，對托卡馬克等離子體的撕裂模不穩(wěn)定性進行了數值模擬。本研究以浙江大學馬志為教授團隊開發(fā)的顯格式求解器CLT 為參考，采用了與之相同的穩(wěn)態(tài)等離子體物理場數據與幾何模型，并基于MPI和并行開源框架PetSc 開發(fā)了新型全隱式求解器。本研究在日本理化學研究所（RIKEN）的HOKUSAI 高性能計算平臺上完成了對所開發(fā)求解器的并行可擴展性測試，并得到了托卡馬克等離子體撕裂模的模擬結果。

1 數值模型及算法實現

針對托卡馬克等離子體的撕裂模不穩(wěn)定性，采用帶Hall 項的非理想單流體MHD 方程［14-15］，并假定等離子體為冷等離子體，即電子壓力與流體壓力相等，可得speed）vα=Bm(μ0a)1/2，上述參量均可由模型確定。

進一步，將方程（1）中電場和電流密度方程分別代入速度方程與磁場方程，即可消去待解變量E、J得到方程（2），從而減少待求解變量的數目，進而節(jié)約存儲開銷。

為了研究托卡馬克中磁流體動力學效應，可取如圖1 所示的幾何模型，其中a為截面圓半徑，R為截面圓圓心所在的大圓半徑，簡稱環(huán)半徑。

圖1 托卡馬克模型Fig.1 Tokamak model

對應的擾動模式，使體系產生不穩(wěn)定性效應。等離子體處于平衡態(tài)時，方程（2）忽略了霍爾效應，且時間導數項為0，可以得到平衡態(tài)時的方程。借助專門針對平衡場求解的NOVA 程序［16］，對圖1所示的模型進行平衡態(tài)的求解，其結果可作為模擬等離子體不穩(wěn)定性的初始分布場。對于圖1所示模型中的撕裂模不穩(wěn)定性，擾動模式如式（3）所示。式中φ為該點環(huán)向角，θ為該點的極向角；η0是背景電阻率，在計算過程中為固定值；ψ是磁通量，由具體模型確定；系數-ψ為磁軸位置，在本文算例中均取-0.095 83，系數δ為常數，取0.03。系數m與n確定了擾動模式的類型，如撕裂模對應m/n=2/1，扭曲模對應m/n= 1/1。

由式（3）可知，擾動模式與時間無關，且服從系數為-ψ與δ的高斯分布。在對方程（2）的實際數值求解過程中，需根據Maxwell 方程組計算出相應的磁場擾動Bp以引入計算，磁場擾動為

結合方程（2）的數值特性與磁流體不穩(wěn)定性模擬中對求解精度的要求，本算法采用了全隱式的4階Crank-Nicholson 格式進行有限差分離散以保證求解的數值精度與穩(wěn)定性，對于離散后的非線性系統(tǒng)，本文采用了Newton-Krylov 法進行迭代求解，其流程如圖2所示。

由圖2可知，求解過程以平衡態(tài)下的物理場作為求解的初始條件，算法采用了兩層嵌套的迭代計算：在每個時間步中，首先根據差分方程計算出雅克比矩陣，然后利用基于并行Krylov子空間的迭代法對所得線性系統(tǒng)進行求解，求得各牛頓步的增量，若該增量小于收斂判據，則推進至下一個時間步的計算。特別地，在計算開始的一定數量的時間步內（在本文中，取計算時間＜1τa作為判據），在每個牛頓步收斂后將由式（4）計算的擾動磁場添加至計算結果中，隨后進入下一個時間步的循環(huán)。

圖2 求解器算法流程Fig.2 Algorithm flow of the solver

2 結果分析

為驗證本文提出的數值模擬方法與求解器的可靠性，選擇m/n= 2/1 的撕裂模算例，取η0=10-5，與CLT 求解器在同模型下的計算結果進行對比，如圖3 所示。由圖3 可明顯觀察到磁場的撕裂模不穩(wěn)定性，本文提出的求解器與CLT 的計算結果在物理場分布上具有較好的一致性。

圖3 70τa時環(huán)向磁場強度計算結果對比Fig.3 Calculation results of circumferential magnetic field strength when 70τa

為進一步驗證本文提出算法的精度與結果的網格獨立性，分別對同模型下取不同空間步長時的計算結果的相對誤差進行校驗。由于本物理模型不存在精確的解析解，在如圖4所示的精度校驗中，本文選用了具有最高空間網格密度的算例作為計算誤差的基準算例，并以式（5）進行相對精度階數的計算，具體誤差數據如表1所示。計算誤差數據選取環(huán)向磁場與環(huán)向速度場的計算結果，并以環(huán)向空間步長為1/8 時的算例為參考，計算其相對精度階數。

表1 不同環(huán)向空間步長下環(huán)向磁場與速度場的相對誤差與相對精度階數Table 1 Relative error and relative accuracy order of Bφ，Vφ under different Δφ

相對精度階數

式中u*表示網格密度最高的近似精確解。本文以環(huán)向空間步長Δφ作為空間網格密度的控制變量，并假設Δφ= 1/64 時的計算結果為近似精確解；可以證明，當趨近于解析解時，由式（5）計算得到的相對精度階數趨向于實際精度階數。

由圖4 與表1 可知，在該托卡馬克磁流體撕裂模不穩(wěn)定性的算例中，算法表現出良好的網格獨立性，且具有近似4階精度，與算法所采用的差分離散格式的理論精度保持一致。

圖4 不同環(huán)向空間步長下各待解變量的相對誤差Fig.4 The relative error of variables under different Δφ

在并行效率與可擴展性方面，使用基于本文提出的求解器進行了測試，測試結果如圖5 所示，其中相對加速比根據式（6）計算得到，直接反映了求解器對于計算資源的利用效率，即弱可擴展性。

圖5 并行可擴展性測試結果Fig.5 Parallel scalability test results

式中t*與n*分別為基準點核心數與計算時間（本測試中，選用160 核心作為基準測試點），tn為核心數為n時的計算時間。

本文采用的測試平臺是日本理化學研究所（RIKEN）的HOKUSAI超級計算機，測試中分別對160、320、640、1 280核心下求解器的計算效率進行了測試。表2 中展示了單步求解的平均耗時數據，相對加速比參照式（6）計算。由圖5 可知，對于一千核以上級別的大規(guī)模并行計算，基于本文提出的算法開發(fā)的求解器依然可以保持60%以上的并行效率，由此可見本求解器在實際的大規(guī)模磁流體仿真求解中具有較高的使用可行性與參考意義。

表2 100τa內單步求解的平均耗時及相對加速度Table 2 Average time consuming and relative acceleration of single step solution in 100τa

3 結 論

本文基于全隱格式有限差分法及大規(guī)模并行模擬提出了一種新型的并行求解MHD方程的算法，并基于該算法對托卡馬克裝置中的磁流體撕裂模不穩(wěn)定性進行了驗證性的仿真模擬，結果與基于全顯格式開發(fā)的CLT求解器吻合良好。

在并行效率方面，求解器在HOKUSAI 高性能計算平臺上進行了可擴展性測試，結果表明：該算法對于跨節(jié)點的并行平臺具有良好的并行效率與可擴展性，可適用于大規(guī)模并行求解。這對于進一步研究更高精度格式與更接近真實托卡馬克裝置模型的仿真計算具有重要的參考意義。