ソース項vs流束の発散

式(1)と式(7)は数値計算上の取り扱いが異なる。 式(1)は、

$$\frac{\partial \rho\bm{v}}{\partial t}+\nabla\cdot[\rho\bm{v}\otimes\bm{v}+p\mathtt{I}]=0$$ (17)

と

$$\frac{\partial \rho\bm{v}}{\partial t}=\rho \bm{g}$$ (18)

に分解して解くことになるため、厳密な運動量保存は成り立たなくなる。 一方、式(7)を解く立場では、

$$(\rho v_x)^{n+1}_{i,j,k}=(\rho v_x)^{n}_{i,j,k} -\frac{\Delta t}{\Delta x}\left(\mathtt{T_g}_{xx,i+1/2,j,k}^{n+1/2} -\mathtt{T_g}_{xx,i-1/2,j,k}^{n+1/2}\right)$$

$$-\frac{\Delta t}{\Delta y}\left(\mathtt{T_g}_{yx,i,j+1/2,k}^{n+1/2} -\mathtt{T_g}_{yx,i,j-1/2,k}^{n+1/2}\right)$$

$$-\frac{\Delta t}{\Delta z}\left(\mathtt{T_g}_{zx,i,j,k+1/2}^{n+1/2} -\mathtt{T_g}_{zx,i,j,k-1/2}^{n+1/2}\right)$$ (19)

のように時間推進するが、時間発展が流束の発散の形で書かれているから、 自己重力流体の問題に対して、外部境界での流束がなければ、厳密に運動量保存が成り立つ。 Athena(Stone+2008)およびver.19までのAthena++(Stone+2020)では、 式(7)の形式を採用しており、 直ぐ下で見る方法で評価した重力テンソルを運動量流束テンソルに加えて 移流する（つまり式(7)を解く）ことが行われていた。 しかし、密度分布に強い偏りが生じた場合などで、 重力場に回転（rot/curl）およびそれを起源とする誤差を生じる問題がMHG21で指摘された (MHG21の図6）。

Stone+2008,2020の重力テンソル(式(8))の評価法は、

$$\mathtt{\tilde{T}_g}{}_{xx,i+1/2,j,k}= \frac{(g_{x,i+1/2,j,k})^2}{8\pi G} -\frac{(g_{y,i+1,j+1/2,k}+g_{y,i,j+1/2,k}+g_{y,i+1,j-1/2,k}+g_{y,i,j-1/2,k})^2}{4^2\,8\pi G}$$

$$-\frac{(g_{z,i+1,j,k+1/2}+g_{z,i,j,k+1/2}+g_{z,i+1,j,k-1/2}+g_{z,i,j,k-1/2})^2}{4^2\,8\pi G}$$ (20)

$$\mathtt{\tilde{T}_g}{}_{yx,i,j+1/2,k}= \frac{g_{y,i,j+1/2,k}(g_{x,i+1/2,j+1,k}+g_{x,i+1/2,j,k}+g_{x,i-1/2,j+1,k}+g_{x,i-1/2,j,k})}{4\,4\pi G}$$ (21)

$$\mathtt{\tilde{T}_g}{}_{zx,i,j,k+1/2}= \frac{g_{z,i,j,k+1/2}(g_{x,i+1/2,j,k+1}+g_{x,i+1/2,j,k}+g_{x,i-1/2,j,k+1}+g_{x,i-1/2,j,k})}{4\,4\pi G}$$ (22)

$$\ldots$$

である。 ここで、整数の添え字はセル中心を、半整数の添え字はセル表面での値を表す。

MHG21は、ソース項に重力を残す式(1)の方法で、 重力テンソルの発散と同じ性質を持つ、 運動量に誤差を生じないスキームを提案した。 MHG21は、重力テンソルの評価を、Stone+2008,2020のものから、対角成分だけ変更し、 以下のような物にした。

$$\mathtt{T_g}{}_{xx,i+1/2,j,k}= \frac{(g_{x,i+1/2,j,k})^2}{8\pi G} -\frac{g_{y,i+1,j+1/2,k}g_{y,i,j+1/2,k}+g_{y,i+1,j-1/2,k}g_{y,i,j-1/2,k}}{2\,8\pi G}$$

$$-\frac{g_{z,i+1,j,k+1/2}g_{z,i,j,k+1/2}+g_{z,i+1,j,k-1/2}g_{z,i,j,k-1/2}}{2\,8\pi G}$$ (23)

$$\mathtt{T_g}{}_{yx,i,j+1/2,k}= \mathtt{\tilde{T}_g}{}_{yx,i,j+1/2,k}$$ (24)

$$\mathtt{T_g}{}_{zx,i,j,k+1/2}= \mathtt{\tilde{T}_g}{}_{zx,i,j,k+1/2}$$ (25)

$$\ldots$$

（代表的対角成分である式(20)を式(23)に変更。） こう取ると、MHG21の補遺Aで示されるように、このテンソルの発散（の$x$成分）は

$$-(\nabla\cdot \mathtt{T_g})_{x,i,j,k}= \frac{g_{x,i+1/2,j,k}+g_{x,i-1/2,j,k}}{2}$$

$$\times \frac{1}{4\pi G}\left(-\frac{g_{x,i+1/2,j,k}-g_{x,i-1/2,j,k}}{\Delta x}\right.$$

$$\left. -\frac{g_{y,i,j+1/2,k}-g_{y,i,j-1/2,k}}{\Delta y}\right.$$

$$\left. -\frac{g_{z,i,j,k+1/2}-g_{z,i,j,k-1/2}}{\Delta z}\right)$$ (26)

となるが、右辺は$g_x\rho$の自然な差分表現になっており、最初の項 $(g_{x,i+1/2,j,k}+g_{x,i-1/2,j,k})/2$は、セル中心での$g_x$が$x$方向に面したセル表面での 加速度 $g_{x,i\pm 1/2,j,k}$の算術平均で与えられること、 2項目は、セル表面の $\bm{g}=(g_x,g_y,g_z)^t$から計算される発散がセル中心の $\nabla\cdot\bm{g}$を 与えることを示している。 つまり、そのようにとれば、右辺のソース項は重力テンソルTの発散で表される ことがわかる。