粒子法(MPS法)の勾配モデルの式に空間の次元数dが出てくる気持ちを考える

はじめに

粒子法入門 ¹ によると, MPS法 (Moving Particle Semi-implicit) において $i$ 番目の粒子に関する物理量 $\phi_{i}$ の近似値 $\left<\nabla\phi\right>_{i}$ は

\left<\nabla\phi\right>_{i} = \frac{d}{n^{0}} \sum_{j \neq i} \left[\frac{\phi_{j} 
- \phi_{i}}{\left| \vec{r}_{j} - \vec{r}_{i} \right|^{2}} \left( \vec{r}_{j} - \vec{r}_{i} \right) w \left( \left| \vec{r}_{j} - \vec{r}_{i} \right| \right) \right] \tag{1}

で表され, 圧力勾配 $\left<\nabla P\right>_{i}$ を求める際などに使用される. ただし, $n^{0}$ は粒子数密度の基準値, $w$ は重み関数, $d$ は空間の次元数である.

この近似式に空間の次元数 $d$ が含まれるのを不思議に感じました. 「直交した軸方向も考慮する必要があるから」という記述は書籍にありましたが, 僕にはその記述をすんなり受け入れられるだけの物理的直感がありません. 正方格子状の粒子の配置を仮定した場合の近似式であるという記述を見つけたのですが², どういう場合にどの程度の精度で成り立つ近似なのでしょう.

書籍にも参照されている論文「MPS法における勾配計算の高精度化とその応用」² の展開に従い式を追った後に, 正方格子状に配置された粒子を考え, 具体的な計算から $d$ と計算結果との関係を調べてみました.

※Safari で閲覧すると数式を含む文章が崩れる可能性があるので, Chrome での閲覧を推奨します.

勾配についての d を含まない表式を求める

位置に依存する物理量と近傍点でのテイラー展開

位置に依存する物理量 $\phi$ を考えます. ある点 $\vec{r}_{i}$ に注目し, その近傍点 $\vec{r}_{j}$ における物理量 $\phi$ のテイラー展開を考えます.

\begin{eqnarray}
\phi(\vec{r}_{j}) &=& \phi(\vec{r}_{i}) + \nabla \phi(\vec{r}_{i}) \cdot \left(\vec{r}_{j} - \vec{r}_{i}\right) + o\left(\left|\vec{r}_{j} - \vec{r}_{i}\right|^{2}\right) \\
  &\simeq& \phi(\vec{r}_{i}) + \nabla \phi(\vec{r}_{i}) \cdot \left(\vec{r}_{j} - \vec{r}_{i}\right) \tag{2}
\end{eqnarray}

まずここで2次以降の項を捨て, 1次の近似をしています.

テイラー展開からの式変形

(2)より,

\nabla \phi(\vec{r}_{i}) \cdot \left(\vec{r}_{j} - \vec{r}_{i}\right) \simeq \phi(\vec{r}_{j}) - \phi(\vec{r}_{i}) ~. 

両辺に $(\vec{r}_{j} - \vec{r}_{i})~/~|\vec{r}_{j} - \vec{r}_{i}|^{2}$ を掛けて,

\left(\nabla \phi(\vec{r}_{i}) \cdot \frac{\left(\vec{r}_{j} - \vec{r}_{i}\right)}{\left|\vec{r}_{j} - \vec{r}_{i}\right|}\right) \frac{\left(\vec{r}_{j} - \vec{r}_{i}\right)}{\left|\vec{r}_{j} - \vec{r}_{i}\right|} \simeq \frac{\phi(\vec{r}_{j}) - \phi(\vec{r}_{i})}{\left|\vec{r}_{j} - \vec{r}_{i}\right|} \frac{\left(\vec{r}_{j} - \vec{r}_{i}\right)}{\left|\vec{r}_{j} - \vec{r}_{i}\right|} ~.

テンソル積の定義 ³ から,

\left[ \frac{\left(\vec{r}_{j} - \vec{r}_{i}\right)}{\left|\vec{r}_{j} - \vec{r}_{i}\right|} \otimes \frac{\left(\vec{r}_{j} - \vec{r}_{i}\right)}{\left|\vec{r}_{j} - \vec{r}_{i}\right|} \right] \nabla \phi(\vec{r}_{i}) \simeq \frac{\phi(\vec{r}_{j}) - \phi(\vec{r}_{i})}{\left|\vec{r}_{j} - \vec{r}_{i}\right|} \frac{\left(\vec{r}_{j} - \vec{r}_{i}\right)}{\left|\vec{r}_{j} - \vec{r}_{i}\right|} ~.

ここで両辺に $ 1~/~n^{0} $ を掛け, 重み関数を掛けた上で $j~( \neq i )$ について和をとります.

\left[ \frac{1}{n^{0}} \sum_{j\neq i}w\left(\left|\vec{r}_{j} - \vec{r}_{i}\right|\right) \frac{\left(\vec{r}_{j} - \vec{r}_{i}\right)}{\left|\vec{r}_{j} - \vec{r}_{i}\right|} \otimes \frac{\left(\vec{r}_{j} - \vec{r}_{i}\right)}{\left|\vec{r}_{j} - \vec{r}_{i}\right|} \right] \nabla \phi(\vec{r}_{i}) \simeq \frac{1}{n^{0}} \sum_{j\neq i}w\left(\left|\vec{r}_{j} - \vec{r}_{i}\right|\right) \frac{\phi(\vec{r}_{j}) - \phi(\vec{r}_{i})}{\left|\vec{r}_{j} - \vec{r}_{i}\right|} \frac{\left(\vec{r}_{j} - \vec{r}_{i}\right)}{\left|\vec{r}_{j} - \vec{r}_{i}\right|} ~ . \tag{2}

ここで

M \equiv \frac{1}{n^{0}} \sum_{j\neq i}w\left(\left|\vec{r}_{j} - \vec{r}_{i}\right|\right) \frac{\left(\vec{r}_{j} - \vec{r}_{i}\right)}{\left|\vec{r}_{j} - \vec{r}_{i}\right|} \otimes \frac{\left(\vec{r}_{j} - \vec{r}_{i}\right)}{\left|\vec{r}_{j} - \vec{r}_{i}\right|} 

とすると,

M \nabla \phi(\vec{r}_{i}) \simeq \frac{1}{n^{0}} \sum_{j\neq i}w\left(\left|\vec{r}_{j} - \vec{r}_{i}\right|\right) \frac{\phi(\vec{r}_{j}) - \phi(\vec{r}_{i})}{\left|\vec{r}_{j} - \vec{r}_{i}\right|} \frac{\left(\vec{r}_{j} - \vec{r}_{i}\right)}{\left|\vec{r}_{j} - \vec{r}_{i}\right|} ~.

両辺に左から $M^{-1}$ を掛けて,

\nabla \phi(\vec{r}_{i}) \simeq M^{-1} \frac{1}{n^{0}} \sum_{j\neq i}w\left(\left|\vec{r}_{j} - \vec{r}_{i}\right|\right) \frac{\phi(\vec{r}_{j}) - \phi(\vec{r}_{i})}{\left|\vec{r}_{j} - \vec{r}_{i}\right|} \frac{\left(\vec{r}_{j} - \vec{r}_{i}\right)}{\left|\vec{r}_{j} - \vec{r}_{i}\right|} ~. \tag{3}

目的の表式が得られました. ここで正方格子状の粒子配置を考えた時に $M^{-1}$ が $dI$ と近似できる, つまり $M$ が $(1/d)I$ と近似できる, というのが(1)式の説明でした. 次に $M$ を具体的に計算すると, どのような値になるのかを見てみました.

2次元の正方格子状の粒子配置から M の値を計算する

2次元の正方格子状の粒子配置を考え, $M$ の具体的な数値を計算します. 2次元で考えているので, $d = 2$ です. 重み関数 $w$, 粒子数密度の基準値 $n^{0}$ には粒子法入門 ¹ に記載されているものを用います.

w(r) = \left\{ \begin{array}{ll}
  \left(\frac{r_{e}}{r}\right) - 1 & \left(r < r_{e}\right) \\
  0 & \left(r \geq r_{e}\right)
\end{array} \right. ~ ,

n^{0} = \sum_{j \neq i}w\left(\left|\vec{r}_{j} - \vec{r}_{i}\right|\right) ~ .

考慮する粒子 $i~(=0)$ は原点 $(0, 0)$ に位置し, 正方格子状に配置された粒子同士の間隔は $1$ として計算を行います.

• $1 < r_{e} < \sqrt{2}$ の場合（影響範囲内に最近接の粒子のみが含まれる場合）

この場合には最近接の粒子のみが計算に寄与し, それより遠い粒子は重み関数によって影響がゼロになります. すべての $j$ について $\left|\vec{r}_{j} - \vec{r}_{i}\right| = 1$ なので,

n^{0} = \sum_{j = 1,2,3,4}w\left(\left|\vec{r}_{j} - \vec{r}_{0}\right|\right) = 4w(1) ~ ,

\begin{eqnarray}
M &=& \frac{1}{n^{0}} \sum_{j = 1,2,3,4}w\left(\left|\vec{r}_{j} - \vec{r}_{0}\right|\right) \frac{\left(\vec{r}_{j} - \vec{r}_{0}\right)}{\left|\vec{r}_{j} - \vec{r}_{0}\right|} \otimes \frac{\left(\vec{r}_{j} - \vec{r}_{0}\right)}{\left|\vec{r}_{j} - \vec{r}_{0}\right|} \\
  &=& \frac{w(1)}{4w(1)} \left((1,0)\otimes(1,0) + (0,-1)\otimes(0,-1) + (-1,0)\otimes(-1,0) + (0,1)\otimes(0,1) \right) \\
  &=& \frac{w(1)}{4w(1)} \left( \left(\begin{array}{cc} 1 & 0 \\ 0 & 0 \end{array}
 \right) + \left(\begin{array}{cc} 0 & 0 \\ 0 & 1 \end{array}
 \right) + \left(\begin{array}{cc} 1 & 0 \\ 0 & 0 \end{array}
 \right) + \left(\begin{array}{cc} 0 & 0 \\ 0 & 1 \end{array}
 \right) \right) \\
  &=& \frac{w(1)}{4w(1)} \left(\begin{array}{cc} 2 & 0 \\ 0 & 2 \end{array}
 \right) \\
  &=& \frac{1}{2} \left(\begin{array}{cc} 1 & 0 \\ 0 & 1 \end{array}
 \right) \\
  &=& \frac{1}{2} I \\
  &=& \frac{1}{d} I ~~ , ~~ \therefore M^{-1} = dI ~ . \tag{4}
\end{eqnarray}

$d$ が出た！ (4)式を(3)式に代入することで, (1)式が得られます.

粒子間の間隔を一般の $L$ とおいたとしても, 変わらず $M^{-1} = dI$ となります. 3次元を考えた場合にも, 計算を想像すると $M^{-1} = dI$ が成立しそうです.

• $r_{e} = 2.1$ の場合

$r_{e} = 2.1$ という数値は, 粒子法入門¹の中で筆者が用いていると記載されている数値です. この場合には最近接に加えて, 次近接, 次々近接の粒子の全部で12個が影響範囲内に含まれます. 対象の粒子からの距離が同じ粒子は同じ色にしていて, 距離は近い順に $1,~\sqrt{2},~2$ であることはすぐに分かると思います. まず $n^{0}$ を計算すると,

n^{0} = \sum_{j = 1,...,12}w\left(\left|\vec{r}_{j} - \vec{r}_{0}\right|\right) = 4\left(w(1) + w(\sqrt{2}) + w(2)\right) \simeq 4\left(1.1 + 0.48 + 0.05\right) = 6.52

となります. $M$ については,（テンソル積の計算は分量が多いので省略します.）

\begin{eqnarray}
M &=& \frac{1}{n^{0}} \sum_{j = 1,...,12}w\left(\left|\vec{r}_{j} - \vec{r}_{0}\right|\right) \frac{\left(\vec{r}_{j} - \vec{r}_{0}\right)}{\left|\vec{r}_{j} - \vec{r}_{0}\right|} \otimes \frac{\left(\vec{r}_{j} - \vec{r}_{0}\right)}{\left|\vec{r}_{j} - \vec{r}_{0}\right|} \\
  &\simeq& \frac{1}{6.52} \left( 2w(1) + 2.0w(\sqrt{2}) + 2w(2) \right) I \\
  &=& \frac{3.26}{6.52} I \\
  &\simeq& \frac{1}{2} I ~~ , ~~ \therefore M^{-1} \simeq d I ~ .
\end{eqnarray}

$d$ が出た！！

勾配モデルと空間の次元数 d

正方格子状に粒子が配置されている場合, 少なくとも2次元では影響範囲の設定によらず, 係数として空間の次元数 $d$ かもしくはそれに非常に近い数値が出てくるということが分かりました. 検算をしていないのですが, 3次元でも同じ結論は得られそうです.

あとがき

$d$ の気持ちがまったく分からなかった状態から, 計算して起源が分かったことで少しは歩み寄れたような気がします. ただやはり物理的な直感が生まれたかどうかというのはまた別の話で, これは流体力学の勉強不足からきているものだと思うので, 気長に調べていこうかなと思います.

記事中に間違いがありましたら, 指摘いただけると幸いです.

1. 粒子法入門 流体シミュレーションの基礎から並列計算と可視化まで ↩ ↩² ↩³

2. MPS法における勾配計算の高精度化とその応用 / この論文では, 勾配モデルとして $d$ を含んだ式を用いずに,テンソル積を直接計算し $M$ に相当する行列とその逆行列を求めることで勾配の計算を行っています. 論文中には勾配の理論値と計算値のずれが視覚的に分かるような図が挿入されています. ↩ ↩²

3. テンソル - http://www.ms.t.kanazawa-u.ac.jp/~design/hojo/zairiki/text/07Complex/02tensor.htm / 「§テンソルの定義，テンソルの成分と座標変換公式」に式変形に含まれる等式が記載されていて, 参考にさせていただきました. ↩