Note：（磁）流体数值模拟程序Zeus与Athena简介

在之前的文章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进行了改良：

1. 加入了协变公式，以便进行不同坐标下的数值模拟工作；
2. 加入了人工粘滞张量；
3. 结合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方程：

式中， 为质量密度， 为等离子体的流速， 为系统的能量密度， 为磁场， 满足泊松方程 为重力势，算子 为拉格朗日微分， 为引力常数。

在这个方程中，Zeus会采用算子分裂法（Split Method）对其进行求解。

利用算子分裂法可以将方程分为源项（Source Step）和传输项（Transport Step），对两项分别采用时域显式的有限差分法进行求解。

其中源项方程变为：

式中 为用来捕捉激波间断面和稳定算法的人工粘滞，具体形式在此不表。

传输项的方程变为：

式中 为网格移动速度，采用欧拉网格时取 ，且 。具体求解过程在此不表。

在程序发展过程中，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方程：

式中， 为动量密度， 为气体压强。具体求解过程不表。

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++具有以下优点：

1. 坐标和网格选择更加灵活；
2. 考虑了广义相对论效应；
3. 显著提高了性能；
4. 改进了源代码的可读性并使之模块化。

图四：本次辐射磁流体力学数值模拟冬季讲习班中的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程序，我们已经有了具有国人自主知识产权的Xi'an-CI程序。

参考文献：

[1]. Stone's 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这一娱乐性质的比赛，不破大陆最差纪录是我们的目标。