天文学中的数值模拟：从PARAMESH说自适应网格
自适应网格这个概念在数值模拟系列的文章中已经提到过好几次了，不过本人之前只是知道大体思路而已，压根没有考虑过具体实现的方法。最近因为在全力突击学习FLASH的缘故，一并去补了补课。因此与本系列的其他文章相比，本文更接近于个人读书笔记，侧重点在网格而非天文应用上。虽说FLASH3之后已经彻底实现了网格和物理的分离，如果只是为了求解具体问题，大可忽略网格部分，不过天知道今后在工作中会不会涉及底层代码呢？毕竟当年第一次了解到FLASH这个东西的时候，本人也没有料到日后会有机会真正接触它，所以这些相关的背景知识也还是多了解些无妨。
自适应网格（Adaptive Mesh Refinement，简称AMR）可以分为两类，其一是解域中格点总数保持不变，只是根据需要随时调整格点具体的位置；其二是实时将现有网格进一步细化，格点总数不再是常量。就本人的经历而言，前者仅见于John D. Anderson, JR.在Computational Fluid Dynamics教材中的介绍，后者才是天文学研究中广泛使用的。鉴于Anderson的教材成书于20年前，当时自适应网格尚处发展初期，计算机性能也十分有限，人们甚至连自适应网格是否真的优于固定网格都不是很确定，书中这一部分内容在如今能有多大的参考价值也是颇值得怀疑的。
在此先回答Anderson在书中提出的疑问，也就是在实际使用中，自适应网格真的存在很大优势么？这是显然的。以FLASH求解二维Sedov点爆炸问题为例，下图给出了不同网格层数下得出的物理量随半径分布，从上到下依次是密度、压强和速度：
图中浅蓝色三角形、绿色方框、红叉以及深蓝色菱形符号分别对应自适应网格取2、4、6、8层的情况，而黑色实线是解析解。网格层数越多，关键区域如激波波面处的网格也就越密，从图中可见所得结果与解析解也越接近。由于网格在有需要的地方更密，差分过程引入的数值振荡或耗散也会大大减弱，数值模拟的真实性和可靠性由此得以提升。当然从理论上说，出于程序编写方便的考虑，大可将整个流场统统划分成最细密的网格，不过其代价就是以指数增长的计算时间了。对于小问题这可能还算不上什么大不了的事情，但是现今天文背景的大型数值模拟工作动辄需要调用国家级超算求解数月，程序的运行时间连带由此产生的高昂费用都是不能不加以考虑的，省时省钱的自适应网格当然更好。
这里本人还要收回先前在FLASH配置一文所说过的一句话，也就是有同学在计算天文课程上用固定网格搞定Sedov解，却未见激波壳层偏离正圆。重翻课程报告之后才发现，其实大家求解的都是一维问题，仅有一位同学模拟了二维。该人的程序用的是大小和数目都固定的网格，所得结果如下，实在是棱角分明到不忍直视啊，看来自适应网格还是很有必要的……
FLASH当前的版本采用了两种自适应网格算法库，分别是基于格点组（block）的PARAMESH与基于格点块（patch）的Chombo。其中前者是NASA与Drexel大学等机构合作开发的，自从上世纪90年代末FLASH问世以来就充当着划分格点的主力，算是立下了汗马功劳，而且从某种意义上说，在2008年因经费原因停止更新之前它也是与FLASH共同发展的；后者则是在FLASH 4之后才刚刚引入，现在尚处测试调整期，据说所得结果不甚理想，而且至今官方小组发表的所有研究成果也都不是基于Chombo网格得出的，因此普通用户要慎用，嗯，其实干脆理解为不要用就对了。
这里问题就来了：block是什么，patch又是什么？这要涉及自适应网格的实现方法。出于方便并行运算分配任务的考虑，PARAMESH（以及其他一些类似的算法库）都是先将解域划分为多个block，每个block由数目固定的格点组成（如对于FLASH而言，默认设置是8乘8乘8总计512个格点）。如下图所示，在block建立之后，网格的细化都要以整个block为单位进行，而且每次细化都是在x、y、z三个方向上各一分为二。同理，细网格的合并也是以block为单位的，每个方向上的两个block合并为一个，所得格点间距随之加倍。
这就是所谓的“八叉树”结构。在三维网格细化后，之前的1个block分裂为8个，这8个子block所覆盖的总范围与母block相当，而且彼此互不重叠。由此就可以对每个子block内的过程单独求解，而不会影响其他区域。作为比较，patch虽然也是解域的一部分，由若干格点组成，但基于patch的算法就不存在这些死板的限制了，它并不要求母子patch覆盖的范围完全一致，也不要求每个patch包含的格点数目必须相同，而且每次细化也毋须对整个patch中的所有格点进行。这样做的代价则是由于patch的覆盖范围变化不定，一切物理过程都必须针对整个解域求算。
至于网格细化的指标，说来也很简单：只要格点两侧某物理量的修正二阶导数超过给定的阈值（换句话说，这个物理量在该格点内的梯度过大），该处网格就需要被细化；反之细网格可以归并变粗，以求节约计算时间。在具体操作中，这个导数的计算当然是依靠对格点差分来实现的。如下式所示，FLASH默认的细化指标不是借助两侧的格点来差分，而是需要用到两侧的第2个格点。现在个人对这种做法的理解是为了方便网格加粗与细化二者，不知是否正确。
这里的指给定物理量，角标
指第
个格点，
是为了筛除数值误差而引入的小量。在具体操作过程中，只要某个block内有一个格点满足细化的标准，那么PARAMESH就会对该block内所有格点统统加细；而只有当block内所有格点的物理量梯度都处在阈值之下时，相应的网格才会变粗。至于细化标准的设定，一般是取两相邻格点物理量出现变化一二成的变化即可，不过具体操作怕是要在实践中总结了。
PARAMESH在设定细化标准方面的高明之处在于，它可以参考多种物理量来加细网格。仍以二维Sedov解为例，以下是其压强分布图。虽然如本文第一张图片所示，在爆发中心处的压强分布较为平坦，但是鉴于这里的温度和速度梯度都比较大，因此格点还是被细化了：
在网格重新划分的时候，需要面对的一大关键问题是疏密格点交界处（以及调整疏密前后同一区域）的通量守恒。由大化小时，细格点的物理量需要通过对粗格点进行插值而得；反之则是对细格点平均得出粗格点的相应数值。若所用坐标系是曲线形式，此时还要考虑格点形状差异带来的权重。为了配合通量守恒的需要，相应的流体代码部分也必须使用形态看上去不是那么整齐的守恒型方程组。而为了应付总能量、速度等某些非守恒量，FLASH也提供了相应的转化功能。
除此之外，并行程序还必须要考虑负载平衡的问题，不能让某个处理器任务过轻或过重，因此这就涉及了调整网格之后计算任务的再分配。PARAMESH使用的是所谓Morton空间填充曲线来实现这一点，并且默认每个子block工作负载相当于母block的2倍。简单说来，这条曲线的作用是保证尽量将相邻的block分配到同一处理器上，减少处理器之间的通信任务量；而子block因为涉及加细过程的插值，需要更多的计算量，因此负载更大。至于更进一步的原理本人就搞不清楚了，还是对计算机底层的运作了解得太少啊……
在重新划分网格的时候，可能是出于方便保持通量守恒、避免引入过多误差的考虑，PARAMESH还额外给出了一条限制：相邻两个block，疏密程度最多相差一级。因此对于激波面这种存在物理量跳变的区域，需要加细的部分绝不是一个窄带。下图是FLASH模拟激波管问题示例，可见波面处的网格分布：
FLASH小组引入Chombo的出发点就是为了克服PARAMESH使用时的诸多限制，如以整个block为单位的细化、对相邻block疏密程度差异的限制以及对block内部格点数的限制等等。因为引入时间不长，可能还要牵扯跨语言编译的问题（Chombo是C++语言写成的，而PARAMESH与FLASH主体一样，都是Fortran90程序），据说Chombo的表现并没有预期的那样好。
对于FLASH用户来说，两种算法库的调用和设置方法差异倒也不是太大。下面不妨了解一下patch的结构分布，还是以二维Sedov解做例子，Chombo网格如下所示，照样是参考了若干不同的物理量来决定格点细化与否：
不过需要注意的是，就算是在Chombo中，也不能单纯细化一个格点，而最好是对某格点周围的区域也一并加细，否则会导致疏密格点的交界处落在物理量梯度过大的区域内，从而引入不必要的麻烦，如简单的激波管可能就会因此出现不收敛的情形。考虑官方使用手册的建议，本人至少在近期是没有打算真正使用Chombo的，只是在此一提。
末了实在是忍不住要吐上几个槽。第一是block这种东西它是有多么的浪费计算资源：
上图只有中心灰色区域的8乘8等于64个格点才是真正要求解的范围，周围每个方向上的4层格点都是所谓的guardcell，理解成相应区域的边界条件就好。对于解域内部的block来说，guardcell需要用周边block的信息来填充；而在解域边界上的block，其guardcell就是真正的边界条件了。因为FLASH默认的流体算法PPM至少需要4个边界格点才能顺利求解，因此guardcell的默认层数也是4层，实际使用时也绝不可以更少。这样看来，为了计算这64个格点的信息，就至少要动用16乘16等于256个格点，三维情况下这一数字会更加可怕。所以这是怎样的一种结构？还是说为了减少浪费，不妨将block内的格点数增加？
第二，下面这个例子算是FLASH中心各位程序员的恶趣味么？FLASH你好歹是天体核爆燃模拟程序，这是要吹什么风洞搞什么飞机呢？
虽说天下流体是一家吧，虽说数值模拟大牛柴田一成有过“哪里有磁场哪里就有我”（同理，哪里有流体哪里也就有他）的豪言壮语吧，不过工程领域更实用的还是基于拉格朗日描述的随体网格，跟天体背景常见的欧拉空间网格差异有点大，于是有哪个工程师会选择FLASH来测试航空器呢？还什么NACA2412号翼形，真的是彻底服了……或者说这只是FLASH在正式求算天体问题之前所经历的测试之一，只有在充分考察了程序对各种已知现象的表现之后才能放任它去探索未知？于是人民群众喜闻乐见的各路不稳定性、自相似解连带诡异的流体实验纷纷登场嗯，其实这些东西也算是提供了参考蓝本吧，用户不妨以相关代码为基础来进行改写，以应付有待解决的问题。（话说流场内添加障碍物的功能其实是新近引入的，莫非这是正式引入之前的测试……）
第三是跟第二相关的。为什么FLASH这种大型模拟程序要进行繁复的测试？前阵系统地通读并推导了各路激波自相似解的一大摞论文，看来看去，还是对纽约大学某位仁兄的某篇短文印象最为深刻。该文写成于10余年前，只有区区4页，至今尚未正式发表，只是上传到了arXiv.org而已，不过有高人对其思路评价甚是不错。这里的槽点在于，该文文末了附上了对激波某类试探性解进行的试探性数值模拟工作，模拟过程不知是用到了怎样不三不四的那么一套代码，说是结果表明，网格划得越细，波后流体的振荡也就越厉害。要命的是该仁兄自己也不知道这诡异的振荡究竟是源自数值误差，还是流体本身的不稳定性，又或者是他的解析解有问题！
好吧，难怪这论文没能发表，理论部分是否靠谱暂且不论，就数值结果而言，连作者自己都搞不清楚可信与否，亏得还好意思写出来。由此可见，选择一套靠谱的模拟程序是多么的有必要啊，模拟程序在正式投入使用之前进行全面的测试也是多么的有必要啊。若是这代码连必须的检测都没能做充足，数值误差各种不明的，用起来还是小心为妙吧！
