Note：（磁）流體數值模擬程序Zeus與Athena簡介

02-13

在之前的文章Note：數值模擬在黑洞吸積中的應用中，我曾介紹到，「數值模擬是目前研究天體物理的一個重要手段，在黑洞天體物理、太陽物理、空間物理等研究領域發揮著重要作用。」加之目前我正在上海天文台參加一個輻射磁流體力學數值模擬冬季講習班，對兩個流體力學（Hydrodynamics,以下簡稱HD）和磁流體力學（Magneto Hydrodynamics,以下簡稱MHD）的開源程序Zeus和Athena的基本情況有了一些簡單的了解，於是便打算在這裡介紹一下這兩個程序。

Zeus：

上世紀七十年代末，Mike Norman為了完成自己的論文，和Jim Wilson一起寫了一串基於Fortran的流體力學數值模擬的代碼，若干年後，就讀於新墨西哥大學的博士生David Clarke在自己的博士論文中對這串代碼做了大量的改進與修改，並賦予了它「Zeus」，這個古希臘神話中的第三代眾神之王的名字。

圖一：Zeus

圖二：現任聖瑪麗大學教授的David Clarke

在那之後，也就是上世紀80年代末，普林斯頓大學的教授James Stone為了解決2維或2.5維MHD問題，將Zeus移植到了Unix平台上，並將其命名為Zeus-2D。在這個過程中，Stone進一步對Zeus進行了改良：

加入了協變公式，以便進行不同坐標下的數值模擬工作；
加入了人工粘滯張量；
結合CT（Constraint Transport）演算法和MOC（Method Of Characteristics）演算法創造了一個更精確的MHD演算法來處理Alfve Wave（阿爾芬波）。

值得一提的是，Zeus原本的應用是在流體力學和磁流體力學上，而Stone又將輻射動力學（Radiation Hydrodynamics，以下簡稱RHD）帶入了Zeus-2D。

圖三：現任普林斯頓天體物理系主任的James Stone

不久之後，為了解決3維問題，Clarke寫出了Zeus-3D。1996年，Norman和他的研究團隊開發出了並行計算版本的Zeus程序，即Zeus-MP，這大大加快了數值模擬的效率，也是本次講習班中被重點介紹的Zeus程序。

首先，我們考慮一個理想的MHD方程：

$frac{D ho}{Dt}+ hoigtriangledowncdot v=0\ hofrac{Dv}{Dt}=-igtriangledown p-pigtriangledownphi+frac{1}{4pi}(igtriangledown imes B) imes B\ hofrac{D}{Dt}(frac{e}{ ho})=-pigtriangledowncdot v\ frac{partial B}{partial t}=igtriangledown imes(v imes B)$

式中， $ho$ 為質量密度， $v$ 為等離子體的流速， $e$ 為系統的能量密度， $B$ 為磁場， $phi$ 滿足泊松方程 $igtriangledown^{2}phi=4pi G ho$ 為重力勢，運算元 $frac{D}{Dt}=frac{partial}{partial t}+vcdotigtriangledown$ 為拉格朗日微分， $G$ 為引力常數。

在這個方程中，Zeus會採用運算元分裂法（Split Method）對其進行求解。

利用運算元分裂法可以將方程分為源項（Source Step）和傳輸項（Transport Step），對兩項分別採用時域顯式的有限差分法進行求解。

其中源項方程變為：

$hofrac{partial v}{partial t}=-igtriangledown ho - hoigtriangledownphi-igtriangledowncdot Q-igtriangledown(frac{B^{2}}{2})+(igtriangledowncdotigtriangledown)B\ frac{partial e}{partial t}=-pigtriangledowncdot v-Q:igtriangledown v$

式中 $Q$ 為用來捕捉激波間斷面和穩定演算法的人工粘滯，具體形式在此不表。

傳輸項的方程變為：

$frac{d}{dt}int_{V} ho dV=-int_{dV} ho(v-v_{g})cdot dS\ frac{d}{dt}int_{V} ho vdV=-int_{dV} ho v(v-v_{g})cdot dS\ frac{d}{dt}int_{V}edV=-int_{dV}e(v-v_{g})cdot dS\ frac{d}{dt}int_{S}BdS=-int_{partial S} ho(v-v_{g}) imes Bcdot dl$

式中 $v_{g}$ 為網格移動速度，採用歐拉網格時取 $v_{g}=0$ ，且 $frac{d}{dt}=frac{partial}{partial t}+v$ 。具體求解過程在此不表。

在程序發展過程中，Stone首先總結了HD、MHD、和RHD問題的演算法，隨後Low Mac於1995年在程序中加入了雙極擴散效應來分析原恆星中的磁場演化，在之後，Stone又加入歐姆耗散來分析MHD並拓展了數值演算法，而L. Dessart又利用偽3D方法對星際風的衍射問題進行了數值模擬研究，改進了RHD演算法。

開發至今，Zeus在吸積盤的磁旋轉不穩定性、恆星形成的磁解耦階段、活動星系核反饋和冷卻流、星系團內介質中的溫度分布結構、太陽爆發過程中各種波的研究等領域中取得了顯著的成果。

Athena：

Athena是Stone於2008年基於C語言開發的高階數值模擬程序，相比Zeus，Athena更側重於MHD問題的數值模擬。

與「父親」Zeus不同的是，Athena基於有限體積法並採用結合了PPM（Piecewise Parabolic Method）演算法、CTU（Corner Transport Upwind）演算法和CT演算法的unsplit Godunov演算法，收斂精度大大提高。

Athena默認求解的是可壓縮、絕熱且非粘性的理想MHD方程：

$frac{partial ho}{partial t}+igtriangledowncdot( ho v)=0\ frac{partial( ho v)}{partial t}+igtriangledowncdot( ho vv-BB+P^{*})=0\ frac{partial e}{partial t}+igtriangledowncdot[(e+P^{*})v-B(Bcdot v)]=0\ frac{partial B}{partial t}-igtriangledown imes(v imes B)=0\ P^{*}=P+frac{Bcdot B}{2}\ e=frac{P}{gamma-1}+frac{ ho(vcdot v)}{2}+frac{Bcdot B}{2}$

式中， $ho v$ 為動量密度， $P$ 為氣體壓強。具體求解過程不表。

Athena開發過程中，Thomas Gardiner和Stone系統的給出了求解2維和3維理想MHD問題的unsplit Godunov演算法，並分析了演算法的CFL（Courant-FriedrichsLewy）條件，並驗證了演算法的可靠性。隨後，為了研究吸積盤的動力學性質，他們又在程序中加入了剪切盒模型的近似模擬。不久，Aaron Skinner和Eve Ostriker將程序從直角坐標系拓展到柱坐標系。而最近，Stone、Kengo Tomida（本次輻射磁流體力學數值模擬冬季講習班主講老師之一）和Christopher White一起開發了基於C++的Athena++，相比於Athena，Athena++具有以下優點：

坐標和網格選擇更加靈活；
考慮了廣義相對論效應；
顯著提高了性能；
改進了源代碼的可讀性並使之模塊化。

圖四：本次輻射磁流體力學數值模擬冬季講習班中的Kengo Tomida

08年至今，Athena在原行星盤、漩渦吸積盤的動力學研究、磁熱不穩定性的非線性行為、激波與星系際介質雲相互作用的磁流體力學研究、吸積盤冕環內的磁重聯過程等領域取得了一些重要成果。

Compare：

Zeus與Athena的差異體現在為：

Zeus基於Fortran，Athena基於C語言；
Zeus使用有限差分法，Athena使用有限體積法；
Zeus不可以求解熱傳導，而Athena可以通過經典顯式迭代來求解；
Zeus通過輻射轉移方程來求解輻射效應，而Athena則通過光學薄輻射函數來求解。

通過對比我們可以看出，

Athena的計算效率要高於Zeus（辣雞Fortran），且供二次開發的介面很多；
面對複雜的物理問題，Athena比Zeus更加好用。

但是這並不意味著全盤否定Zeus。

Summary：

目前，國內天體物理的MHD數值模擬仍然非常薄弱，因為我們使用的數值模擬程序都是由國外大學或研究所的團隊多年以來開發出來的相對比較成熟的平台，但實際上在對具體的物理問題進行研究時仍然會有各種BUG出現。因此，發展一套國產的數值模擬程序勢在必行。

最後說句題外話，說到這裡，我想起了我在我們學校所做的量子化學的創新項目，也是數值模擬工作，只不過不同的是，目前在量子化學領域，相對與國際上常用的Gauss View程序，我們已經有了具有國人自主知識產權的Xian-CI程序。

參考文獻：

[1]. Stones ZEUS Code Home Page

[2]. ZEUS-3D, v3.6 (ICA)

[3]. Athena++

[4]. Stone J M. The Athena MHD Code: Extensions, Applications, and Comparisons to ZEUS[J]. 2009, 406:277.

P.S.周五下午開始到周天，將進行分課題數值模擬練習，我想選南京大學陳鵬飛教授的課題，但是不知道能不能選上，選不上的話就選清華大學白雪寧老師（Stone的弟子，參與Athena的開發）的課題。

P.P.S.同時，周六下午六點到周日下午六點講和@蔡家麒、@胡子昂、 @薩塔妮亞以及一位不知道知乎ID的大神一起參加International Olympiad in Theoretical Physics這一娛樂性質的比賽，不破大陸最差紀錄是我們的目標。