カーネル補間法 - YuriOku/1D_SPH GitHub Wiki

ここでは、カーネル関数の空間微分 $\nabla_i W$ を計算する2つの方法について述べます。

Standard gradient

1つ目は単純に微分を計算する方法で、

$\nabla_i W\left(r_{ij}, h \right) = \frac{\partial r_{ij}}{\partial \mathbf{r}i} \frac{\partial W(r{ij}, h)}{\partial r_{ij}} = \mathbf{e}{ij} \frac{\partial W(r{ij}, h)}{\partial r_{ij}}$

となります。ここで、 $r_{ij} = \left| \mathbf{r}_i - \mathbf{r}j \right|,\ \mathbf{e}{ij} = \frac{\mathbf{r}_i - \mathbf{r}_j}{\left|\mathbf{r}_i - \mathbf{r}_j \right|}$ です。 $\frac{\partial W}{\partial r}$ についてはカーネル関数をご覧ください。

Integral Approach

Integral Approach (IA)は、García-Senz et al. (2012)によって提案された手法で、従来のSPH法が苦手とする Rayleigh-Taylor 不安定性や Kelvin-Helmholtz 不安定性をより精度よく扱うことができます。まず、積分

$I(\mathbf{r}) = \int \left(f(\mathbf{r}^') - f(\mathbf{r}) \right) \left(\mathbf{r}^' - \mathbf{r} \right) W\left(\left|\mathbf{r}^' - \mathbf{r} \right|, h \right) d{r^'}^d \cdot\cdot\cdot\cdot\cdot (1)$

を考えます。ここで $d$ は空間の次元です。 $f(\mathbf{r}^')$ を点 $\mathbf{r}$ 周りで Taylor 展開すると、

$f(\mathbf{r}^') - f(\mathbf{r}) = \nabla f (\mathbf{r}) \cdot (\mathbf{r}^' - \mathbf{r}) + \mathcal{O}({\mathbf{r}^'}^2) \cdot\cdot\cdot\cdot\cdot (2)$

となるので、高次項は無視して式(2)を式(1)に代入し、 $\nabla f(\mathbf{r})$ について解くことを考えます。

1次元

1次元では、式(1)は

$I(r) = \frac{\partial f(r)}{\partial r} \int \left(r^' - r \right)^2 W\left(\left|r^' - r \right|, h \right) dr^' \cdot\cdot\cdot\cdot\cdot (3)$

となるので、式(3)を $\frac{\partial f(r)}{\partial r}$ について解き $I(r)$ に式(1)を代入すると、

$\frac{\partial f(r)}{\partial r} = \frac{\int \left(f(r^') - f(r) \right) \left(r^' - r \right) W\left(\left|r^' - r \right|, h \right) dr^'}{\int \left(r^' - r \right)^2 W\left(\left|r^' - r \right|, h \right) dr^'} \cdot\cdot\cdot\cdot\cdot (4)$

が得られます。カーネル近似によってこの積分をカーネル和に置き換えると

$\frac{\partial f(r_i)}{\partial r} = \frac{\sum_j dV_j \left(f(r_j) - f(r_i) \right) \left(r_j - r_i \right) W\left(\left|r_j - r_i \right|, h \right)}{\sum_j dV_j \left(r_j - r_i \right)^2 W\left(\left|r_j - r_i \right|, h \right)} \cdot\cdot\cdot\cdot\cdot (5)$

となります。この微分は、場が線形であるとき正しい値となります。例えば、

$\rho_j = \rho_i + a(r_j - r_i)$

という密度場の微分は

$\begin{align*} \frac{\partial \rho_i}{\partial r} &= \frac{\sum_j dV_j \left(\rho_j - \rho_i \right) \left(r_j - r_i \right) W\left(\left|r_j - r_i \right|, h \right)}{\sum_j dV_j \left(r_j - r_i \right)^2 W\left(\left|r_j - r_i \right|, h \right)}\ &= \frac{\sum_j dV_j a \left(r_j - r_i \right)^2 W\left(\left|r_j - r_i \right|, h \right)}{\sum_j dV_j \left(r_j - r_i \right)^2 W\left(\left|r_j - r_i \right|, h \right)}\ &= a \end{align*}$

となります。この手法は Linear-exact gradient (Price 2004) と同値のものです (Rosswog 2015)。しかしこの手法では運動方程式が $i,,j$ の交換に対して反対称とならないので、運動量が保存されません。反対称とするためには、

$\begin{align*} &\frac{\partial f_i}{\partial r} = \sum_j f_j G_{ij}\ &(G_{ij} = -G_{ji}) \end{align*}$

という形になっている必要があります。そこで Integral Approach では

$\sum_j dV_j \left(\mathbf{r}_j - \mathbf{r}_i \right) W\left(\left|\mathbf{r}_j - \mathbf{r}_i \right|, h \right) = 0$

を仮定します。これは、粒子の分布に偏りがなく一様であるという仮定です。この仮定を用いて式(1)を

$I(\mathbf{r}) \sim \sum_j dV_j f(\mathbf{r}_j) \left(\mathbf{r}_j - \mathbf{r}_i \right) W\left(\left|\mathbf{r}_j - \mathbf{r}_i \right|, h \right) \cdot\cdot\cdot\cdot\cdot (7)$

と近似します。結果として、微分(5)は

$\frac{\partial f(r_i)}{\partial r} = \frac{\sum_j dV_j f(r_j) (r_j - r_i ) W\left(\left|r_j - r_i \right|, h \right)}{\sum_j dV_j (r_j - r_i )^2 W\left(\left|r_j - r_i \right|, h \right)}$

となります。これを通常のSPHでの微分

$\nabla_i f(r_i) = \sum_j f(r_j) dV_j \nabla_i W(| r_i - r_j|, h)$

と見比べると、Integral Approach でのカーネル関数の空間微分

$\nabla_i W\left(\left| r_i - r_j \right|, h\right) = \frac{(r_j - r_i ) W\left(\left|r_j - r_i \right|, h \right)}{\sum_j dV_j (r_j - r_i )^2 W\left(\left|r_j - r_i \right|, h \right)}$

が得られます。

3次元

3次元でも基本的には同じですが、式(1)は

$\left[ \begin{array}{c} I_1 (\mathbf{r})\ I_2 (\mathbf{r})\ I_3 (\mathbf{r}) \end{array} \right] = \int \left( \left[ \begin{array}{c} \partial f(\mathbf{r})/\partial x_1\ \partial f(\mathbf{r})/\partial x_2\ \partial f(\mathbf{r})/\partial x_3 \end{array} \right] \cdot \left[ \begin{array}{c} x_1^' - x_1\ x_2^' - x_2\ x_3^' - x_3 \end{array} \right] \right) \left[ \begin{array}{c} x_1^' - x_1\ x_2^' - x_2\ x_3^' - x_3 \end{array} \right] W\left(\left| \mathbf{r}^' - \mathbf{r} \right|, h \right) d{r^'}^3$

となります。ただし、 $\mathbf{r} = x_1 \mathbf{i} + x_2 \mathbf{j} + x_3 \mathbf{k}$ 、 $I(\mathbf{r}) = I_1 (\mathbf{r})\mathbf{i} + I_2 (\mathbf{r})\mathbf{j} + I_3 (\mathbf{r})\mathbf{k}$ です。これを $\nabla f(\mathbf{r})$ について解くと、

$\left[ \begin{array}{c} \partial f(\mathbf{r})/\partial x_1\ \partial f(\mathbf{r})/\partial x_2\ \partial f(\mathbf{r})/\partial x_3 \end{array} \right] = \left[ \begin{array}{ccc} \tau_{11} & \tau_{12} & \tau_{13}\ \tau_{21} & \tau_{22} & \tau_{23}\ \tau_{31} & \tau_{32} & \tau_{33} \end{array} \right]^{-1} \left[ \begin{array}{c} I_1 (\mathbf{r})\ I_2 (\mathbf{r})\ I_3 (\mathbf{r}) \end{array} \right]$

となります。ここで

$\tau_{ij} = \int (x_i^' - x_i)(x_j^' - x_j) W\left(\left| \mathbf{r}^' - \mathbf{r}\right|, h\right) d{r^'}^3$

で、これをカーネル和で表すと

$\tau_{ij} = \sum_b dV_b (x_{i,b} - x_{i,a})(x_{j,b} - x_{j,a}) W\left(\left| \mathbf{r}_b - \mathbf{r}_a\right|, h\right)$

となります。式(12)の右辺第一項の逆行列は

$C = \frac{1}{\tau_{11}\tau_{22}\tau_{33} - \tau_{11}\tau_{23}\tau_{32} + \tau_{12}\tau_{21}\tau_{33} -\tau_{12}\tau_{23}\tau_{31} + \tau_{13}\tau_{21}\tau_{32} - \tau_{13}\tau_{22}\tau_{31}} \left[ \begin{array}{ccc} \tau_{22}\tau_{33} - \tau_{23}\tau_{32} & \tau_{13}\tau_{32} - \tau_{12}\tau_{33} & \tau_{12}\tau_{23} - \tau_{13}\tau_{22}\ \tau_{23}\tau_{31} - \tau_{21}\tau_{33} & \tau_{11}\tau_{33} - \tau_{13}\tau_{31} & \tau_{13}\tau_{21} - \tau_{11}\tau_{23}\ \tau_{21}\tau_{32} - \tau_{22}\tau_{31} & \tau_{12}\tau_{31} - \tau_{11}\tau_{32} & \tau_{11}\tau_{22} - \tau_{12}\tau_{21} \end{array} \right]$

と表されます。また、積分(1)の i 成分

$I_i(\mathbf{r}) = \int \left(f(\mathbf{r}^') - f(\mathbf{r}) \right) \left(x_i^' - x_i \right) W\left(\left|\mathbf{r}^' - \mathbf{r} \right|, h \right) d{r^'}^3$

から $f(\mathbf{r})$ を消去してカーネル近似すると、

$I_i(\mathbf{r}a) = \sum_b dV_b f(\mathbf{r}b) \left(x{i,b} - x{i,a} \right) W\left(\left|\mathbf{r}_b - \mathbf{r}_a \right|, h \right)$

となります。関数 $f(\mathbf{r})$ の微分は式(12)から

$\frac{\partial f(\mathbf{r}a)}{\partial x_i} = \sum{l = 1}^3 C^{il}I_l(\mathbf{r}_a)$

と表され、これに式(16)を代入して通常のSPHの微分と比較すると、カーネル関数の微分は

$\frac{\partial W\left(\left|\mathbf{r}b - \mathbf{r}a\right|, h\right)}{\partial x_i} = \sum{l = 1}^3 C^{il} (x{l,b} - x_{l,a}) W\left(\left|\mathbf{r}_b - \mathbf{r}_a\right|, h\right)$

となります。

参考文献

García-Senz, D., Cabezón, R. M., and Escartín, J. A., “Improving smoothed particle hydrodynamics with an integral approach to calculating gradients”, Astronomy and Astrophysics, vol. 538, 2012. doi:10.1051/0004-6361/201117939.
Price, D. J., “Magnetic fields in Astrophysics”, PhDT, 2004.
Rosswog, S., “Boosting the accuracy of SPH techniques: Newtonian and special-relativistic tests”, Monthly Notices of the Royal Astronomical Society, vol. 448, no. 4, pp. 3628–3664, 2015. doi:10.1093/mnras/stv225.