天文学中的数值模拟：简介篇
关于天文中的数值模拟，其实是很早就订下的一个打算，本站会有一系列的文章。不过后面的具体安排还没有考虑好，先写个简要介绍，然后慢慢来。本人初次正面接触数值模拟的概念是在大一或是大二，具体时间已经记不太清楚了，总之是听了本系陈鹏飞老师的一次简介性质的讲座，算是沾了一点皮毛。之后自己看过一些相关的教材，直到大四，才旁听了计算天文的课程，对数值模拟的基本方法有了个大概的了解，求解过几个简单问题，当是初步入门，也才刚刚有了着手写下去的底气，那么就在之前的仪器与历法两个系列临近尾声的时候开始动笔吧。
先说个引子。差不多是去年年底，在arXiv.org看到了一篇关于伽玛暴瞬时辐射的论文，说是已经被《天体物理学杂志》接收。其实这篇文章数值模拟的味道不重，但问题是，作者之一是京都大学的柴田一成教授，陈鹏飞老师博士后期间的导师，一大主攻方向是太阳物理数值模拟研究，也算是圈子里的名人。当时还纳闷半天，心说这也差得太远了些，难道是此君什么时候改行伽玛暴了？或者此柴田非彼柴田乎？后来就此事还专门问过陈老师，答复是，此君曾经口出豪言壮语，曰，哪里有磁场，哪里就有我！所以，他做些什么都属正常，云云。当即听得本人差点崩溃……
之后陈老师开始讲，对于数值模拟来说，柴田的那句话并不算太过分，至少在天文常用的磁流体模拟上，事实就是如此。因为出发点无外乎磁流体基本方程组，对于不同的天文现象，主要的不同只是初始条件和边界条件。倘或真的把基本方法折腾得很熟练，想达到“哪里有磁场，哪里就有我”的境界也并非没有可能。由此也可以一窥当代数值模拟的威力了。
还有一种天文常用的主要方法是粒子模拟，说穿了就是依据牛顿定律来求解粒子运动，如模拟N个粒子在自引力下的演化，规模可以做到大得惊人，出名的成果也不少。有个搞非线性动力学的同学一度就是整日研习相关算法，别人与他见面的寒暄也变成了“你的N-body看得如何了”。话说他所在的小组也是系内使用计算机最狠的一个：随便一个程序，动辄占用几十个CPU的资源，再跑上几个月，众人皆是早已见怪不怪。
不必多说也可以料到，这豪言绝对是说来容易做来麻烦到极点的。天文现象，特别是需要用模拟来解决的天文现象，往往有着极端的环境和恼人的边界条件，要不就是方程组极其庞大而且高度非线性，更可能的情况则是二者兼而有之，比如耀斑，比如宇宙大尺度结构的演化，比如活动天体的喷流，比如恒星或行星系的形成，比如超新星爆发。当然了，倘或方程简单，设法得出解析解岂不是更好，退而求其次的话也可以做展开做近似求个大概么，也就没有多大的必要让计算机搀和进来了。（PS：其实这话更狂，现在看看20世纪早期到中期老一代天体物理学家的工作，比如史瓦西，比如钱德拉塞卡，比如Zel’dovich和他的弟子们，用到的那些技巧，进行的那种计算量，实在是让人倒吸凉气啊……）
好吧，先贴点模拟图吊吊眼球，歇口气接着讲：
马克斯·普朗克研究所的千年模拟，这是迄今规模最大的N体宇宙学模拟。（图片来源：Springel et al. 2005）
马普所的星系碰撞模拟。（图片来源：Tiziana Di Matteo, Volker Springel and Lars Hernquist 2005）
马普所的双中子星对撞模拟，涉及伽玛暴短暴的形成。（图片来源：R. Oechslin and H.-Th. Janka 2005）
Ia型超新星的爆发模拟，由芝加哥大学利用FLASH完成。（图片提供：芝加哥大学）
太阳风与地球磁场相互作用的模拟。（图片提供：Finnish Meteorological Institute）
原行星盘模拟，传说此图让相关研究者一举成名。（图片来源：Department of Astronomy, University of Virginia）
20世纪60年代兴起的数值模拟的出现显然与计算机的发展脱不开干系。当年这一方法初问世的时候，老一辈置疑不断。这也容易理解，数值计算本身就是一种近似，何况计算过程中的舍入误差、算法缺陷和其他因素还会给出非物理解干扰结论，举个例子，对最简单的一阶波动方程，采用迎风格式，对计算步长的绝大部分取值，结果都会出现本来并不存在耗散，明显与波的实际行为不符；而对于另一些问题，情况正好相反：如果不引入人为耗散，数值解也会出现反常行为。
由于一开始计算机性能较差，早年的数值模拟规模不大，结果也是相当粗糙，更不会有人过多考虑“花架子”一般的可视化，以上列出的绚烂图片那是不可能的事情。陈鹏飞老师曾经给大家看过柴田70年代末80年代初读研究生的时候做过的工作，十分怀疑该图是根据计算结果人工绘制的。现在如果谁的模拟还搞成这个样子，估计八成会被退回重做的：
但正如当年某些人的预见，数值模拟的确兴起了，并且日渐成为热点，这不仅仅针对天文，对气象、地理、工程等学科也一样。到如今，方便好用的程序包和可视化软件层出不穷，一流研究机构纷纷配备超级计算机做模拟用，算法和后续数据处理的研究方兴未艾，很多工作从当年的一维二维逐步扩展到三维立体模拟，甚至还有GRAPE这样的模拟专用计算机问世。除了计算机硬件性能的提升和并性计算技术的发展，学术上的理由也同样简单：自然界中，能解析求解（甚至包括可以进行小量展开近似）的东西实在是极少数。涉及非线性方程的问题让专家纷纷头大不已。N大于等于3的多体问题是不可积的。湍流基本是不可解的。就连简化到极点的恒星模型——多方球，描述它的Lane-Emden方程在绝大多数情况下也是不存在解析解的。但研究还是需要进行的，问题也是有待阐明的，哪怕只有数值解，不也是比一无所知好得多吗？何况，经验在时刻积累，算法在逐渐改进，硬件在迅速发展，进而数值解的精度一天天地提高，越来越多的模拟结果也为观测证实，所以现今基本没有谁会无视模拟工作了。
谈到算法问题，其实就是求解微分方程组的方法，不打算介绍太多，因为目前可以买到多本不错的教材，何况这也不是三句两句就可以说完的，再说也须观者配合，静心推导方可知其一二。对于粒子模拟，实话实说，本人没有真正做过，一切从简。直接计算粒子之间的作用力是最直观但也相当耗时的，至于其他，只知道Tree算法的大概，貌似是将作用力分为近场与远场，近场逐一计算，远场按区域算总数，还有Top-Down与Bottom-Up之分。此外又有什么粒子网格之类，实在是没有写下去的胆量，就不再多说了，想有所了解的可以参考这里的简要介绍。
简单说说磁流体模拟。方程从实际情况出发建立，一般为方便起见还要先经过无量纲化处理。一个关键步骤是连续量的离散化，如天文常用的差分法是将解域空间划分为网格，以格点数值代入差分方程，代替微分的求解。化差分的具体方法很多，有单步与多步之别，可分为显式与隐式两大类，区别在于是否可以直接通过前N步的已知数据求出第N+1步的结果。对于每个格点，显式可以由单个表达式求出数值，基本格式如迎风格式、蛙跳格式、Lax-Wendroff格式、Rubin-Burstein格式、MacCormack格式等，它们有各自的精度、适用方程种类和稳定性条件，尽可按需选择；隐式需要求解线性方程组，其优势在于稳定且精度较高，但对多维情况运算量猛增甚至几近不可解，因而往往只是化为局部隐式处理。近年又开发出了一系列的高精度算法，在很大程度上摆脱了传统差分的缺陷。这里网格的划分亦是颇有讲究，除了满足稳定性条件之外，既要保证足够的分辨率显现细节，还需兼顾运算速度，格点间距在解域内也应平稳变化。目前已开发有专用工具完成网格划分工作，一些大规模模拟还会采用自适应网格技术，加上种种选择条件，如果某区域有精细结构出现，自动将网格进一步细分，其他区域则一带而过，以达到速度与精度兼顾的目的。至于边界条件的设定更是麻烦，论直观物理条件当然是首选，追求稳定的话则可以考虑周期性条件（不过不一定要求模拟对象具有周期性），也许这还与研究者的经验或习惯有关吧。
对于磁流体模拟的有限元法、谱方法，以及基于Vlasov方程的方法，连带天体物理背景下往往要加以考虑的辐射转移、对流传导、核反应等因素，也因为没有接触太多而不敢多说。以拉格朗日型方程为基础的的有限元法在天文上用得似乎要比网格少些，可能是由于天体环境下不需要考虑过于复杂的边界的缘故，至少学习计算天文的时候没怎么接触过；光滑粒子流体模拟当下本人只了解个皮毛，总的来说是要用有分布的粒子来对模拟对象进行近似，毋须网格，对于物理条件变化较大或者存在穿越的场合比较合适。好在这一方法有教材可看，对不住想对其有所了解的同行同好了，细节自学吧……
除了FLUENT等通用的商业软件之外，天文专用的模拟程序包也不少，前面提到的FLASH就是一例，可以进行超新星爆发的三维模拟，但程序似乎由于种种非天文的原因不对非美国公民开放，只有零星的子程序放出，不细说也罢。此外还有HYDRA之类，可以在这里找到全面的软件包列表和简介附带相关链接。自己曾经啃过一个名叫PMFAST的N体粒子网格模拟程序代码，只是搞得一知半解，后来不了了之。不过么，既然连专用计算机都有了，专用程序包的出现也就显得很正常了，况且大型程序的开发也不是依靠一两个天文学家的力量在短时间内就可以搞得掂的。
需要说的套话而且也是实话是，数值模拟并非万能。除了算法和计算机本身的误差，模拟也不是可以将一切影响因素都包括在内的，这不仅是出于简便的考虑，也有对某些过程的了解不够充分的问题。如大尺度模拟一旦涉及恒星形成，立即变得复杂无比；超新星爆发的具体机制至今仍属悬而未决。何况当前基于模拟发表的论文往往不会也不可能给出源程序，他人进行结果验证的难度很大。所以，再次使用套话就是，应该让理论、模拟与观测结合，更重要的是，千万不能忽视背后的物理本质。正因为如此，当下在许多一流的研究团队中，对模拟结果分析投入的精力并不亚于模拟过程本身。
末了放个小软件，可以模拟吸积盘的FITDisk，用户要做的只是设置初始条件，运行速度也还算可以，有兴趣就玩玩吧~~~
