【粒子法 MPS法】アルゴリズムの改善：GC法 勾配モデルの高精度化

上記のnoteで解説したHS法に加えて、下記の書籍で詳しく解説してある標準MPSの高精度化に着手しているよ。目次は次のとおりだよ。

書籍：粒子法:連続体・混相流・粒状体のための計算科学

P.83 HS法 ポアソン方程式生成項の高精度化
P.85 HL法 ラプラシアンモデルの高精度化
P.87 ECS法 ポアソン方程式生成項の誤差補償
P.92 GC法　勾配モデルの高精度化
P109 DS法　人工斥力導入による安定化
P118 SPP法　表面粒子の境界条件適正化
P124 WPP法　壁面における境界条件適正化

粒子法:連続体・混相流・粒状体のための計算科学

書籍の順番に従って、勾配モデルの高精度化であるGC法の導入を行うよ。粒子法は格子点が空間的に不均一であるため、偏微分方程式を粒子法の計算アルゴリズムに書き換えるときのグラディエントモデル並びにラプラシアンモデルは、なんとテーラー展開の１次の精度すら保証していないことが知られているね。そこで登場した、テーラー展開の１次精度を満たすグラディエントモデルやラプラシアンモデルの改良版がGC法と呼ばれるよ。

２つの粒子番号$${i, j}$$の位置ベクトルを$${\boldsymbol{r}_i, \boldsymbol{r}_j}$$、その地点における任意のスカラー関数の値を$${\phi_i, \phi_j }$$と表したときに$${\phi_i}$$と$${\phi_j}$$の関係はテーラー展開を用いて次のように表すことができるね。

$$
\phi_j = \phi_i +\nabla\phi_i\cdot \boldsymbol{r}_{ij} +{\cal O}(h^2)
$$

$${ \boldsymbol{r}_{ij} = \boldsymbol{r}_{j} - \boldsymbol{r}_{i} }$$だよ。２次以降を無視して次のように変形するよ。

$$
\frac{ \phi_j -\phi_i}{|\boldsymbol{r}_{ij} |}= \nabla\phi_i\cdot \frac{ \boldsymbol{r}_{ij}}{|\boldsymbol{r}_{ij} |}
$$

グラディエントモデルの形に近づけるために両辺に$${\boldsymbol{r}_{ij}/|\boldsymbol{r}_{ij} |\omega_{ij}}$$を掛けて、右辺を変形するよ。

$$
\frac{ \phi_j -\phi_i}{|\boldsymbol{r}_{ij} |} \,\frac{ \boldsymbol{r}_{ij}}{|\boldsymbol{r}_{ij} |}\,\omega_{ij} = \left( \nabla\phi_i\cdot \frac{ \boldsymbol{r}_{ij}}{|\boldsymbol{r}_{ij} |} \right) \frac{ \boldsymbol{r}_{ij}}{|\boldsymbol{r}_{ij} |}\,\omega_{ij}\\
                                                = \left(\frac {\boldsymbol r_{ij}\otimes\boldsymbol r_{ij}}{|\boldsymbol r_{ij}|^2}\right) \nabla \phi _i\, \omega _{ij}
$$

ちなみに、$${\otimes }$$は以下に示すテンソル積だよ。

$$
(\boldsymbol{a}\otimes \boldsymbol{b})\boldsymbol{c}:= \boldsymbol{a} ( \boldsymbol{b}\cdot \boldsymbol{c} ) = ( \boldsymbol{c}\cdot \boldsymbol{b} )\boldsymbol{a}
$$

さらにグラディエントモデルに近づけるために両辺を$${n^{(0)}}$$で割って$${j}$$インデックスで和をとると

$$
\sum\limits_{j\ne i}\frac{1}{n^{(0)}}\,\frac{ \phi_j -\phi_i}{|\boldsymbol{r}_{ij} |} \,\frac{ \boldsymbol{r}_{ij}}{|\boldsymbol{r}_{ij} |}\,\omega_{ij} = \sum\limits_{j\ne i}\frac{1}{n^{(0)}}\left(\frac {\boldsymbol r_{ij}\otimes\boldsymbol r_{ij}}{|\boldsymbol r_{ij}|^2}\right) \omega _{ij} \,\nabla \phi _i
$$

となり、$${\nabla \phi _i }$$について解くと次のようになるね。

$$
\nabla \phi _i= \left[ \sum\limits_{j\ne i}\frac{1}{n^{(0)}}\left(\frac {\boldsymbol r_{ij}\otimes\boldsymbol r_{ij}}{|\boldsymbol r_{ij}|^2}\right) \omega _{ij} \right]^{-1} \sum\limits_{j\ne i}\frac{1}{n^{(0)}}\,\frac{ \phi_j -\phi_i}{|\boldsymbol{r}_{ij} |} \,\frac{ \boldsymbol{r}_{ij}}{|\boldsymbol{r}_{ij} |}\,\omega_{ij}
$$

この粒子$${i}$$の近傍粒子位置における平均をとった値に次元数を掛けて改めて粒子法におけるグラディエント$${\langle \nabla \phi \rangle_i}$$と考えるんだね。

CG法グラディエントモデル

$$
\langle \nabla \phi \rangle_i =\frac{\rm dim}{n^{(0)}} \sum\limits_{j\ne i}\frac{ \phi_j -\phi_i}{|\boldsymbol{r}_{ij} |} \,\frac{ \boldsymbol{C}_{i}\cdot \boldsymbol{r}_{ij}}{|\boldsymbol{r}_{ij} |}\,\omega_{ij}
$$

$$
\boldsymbol{C}_{i}= \frac{1}{\rm dim}\left[ \sum\limits_{j\ne i}\frac{1}{n^{(0)}}\left(\frac {\boldsymbol r_{ij}\otimes\boldsymbol r_{ij}}{|\boldsymbol r_{ij}|^2}\right) \omega _{ij} \right]^{-1}
$$

$${\boldsymbol{C}_{i}}$$は粒子配置の不揃いに起因する誤差を補正するテンソルって感じだね。粒子配置が当方的な場合にこのテンソルは単位テンソルになって、元のMPS法のグラディエントモデル

$$
\langle \nabla \phi \rangle_i =\frac{{\rm dim}}{n^{(0)}} \sum\limits_{j\ne i}\frac{ \phi_j -\phi_i}{|\boldsymbol{r}_{ij} |} \,\frac{ \boldsymbol{r}_{ij}}{|\boldsymbol{r}_{ij} |}\,\omega_{ij}
$$

と一致するわけだね。

続いて、CG法ラプラシアンモデルを行うための前段として、ベクトル量の勾配を導出するよ。位置に依存する任意のベクトル量$${\boldsymbol{u}}$$のテーラー展開を次のように表すよ。

$$
\boldsymbol{u}_j = \boldsymbol{u}_i + \nabla \boldsymbol{u}_i \cdot \boldsymbol{r}_{ij} + {\cal O}(h^2)
$$

ここからスカラー関数のときと同様の手続きで、最終的に $${\nabla \boldsymbol{u}_i }$$ について解くよ。

$$
\frac{ \boldsymbol{u}_j - \boldsymbol{u}_i }{|\boldsymbol{r}_{ij} |}= \nabla \boldsymbol{u}_i \cdot \frac{ \boldsymbol{r}_{ij}}{|\boldsymbol{r}_{ij} |}
$$

$$
\frac{ \boldsymbol{u}_j - \boldsymbol{u}_i }{|\boldsymbol{r}_{ij} |} \otimes\frac{ \boldsymbol{r}_{ij}}{|\boldsymbol{r}_{ij} |}\,\omega_{ij} = \nabla \boldsymbol{u}_i \cdot\left(\frac {\boldsymbol r_{ij}\otimes\boldsymbol r_{ij}}{|\boldsymbol r_{ij}|^2}\right) \omega _{ij}
$$

$$
\frac{1}{n^{(0)}}\sum\limits_{j\ne i}\,\frac{ \boldsymbol{u}_j - \boldsymbol{u}_i }{|\boldsymbol{r}_{ij} |} \otimes\frac{ \boldsymbol{r}_{ij}}{|\boldsymbol{r}_{ij} |}\,\omega_{ij} = \nabla \boldsymbol{u}_i \cdot \left[\frac{1}{n^{(0)}}\sum\limits_{j\ne i} \left(\frac {\boldsymbol r_{ij}\otimes\boldsymbol r_{ij}}{|\boldsymbol r_{ij}|^2}\right) \omega _{ij} \right]
$$

$$
\langle \nabla \boldsymbol{u}\rangle_i = \left[ \frac{1}{n^{(0)}}\sum\limits_{j\ne i}\,\frac{ \boldsymbol{u}_{ij} }{|\boldsymbol{r}_{ij} |} \otimes\frac{ \boldsymbol{r}_{ij}}{|\boldsymbol{r}_{ij} |}\,\omega_{ij} \right] \cdot \boldsymbol{B}_i
$$

$$
\boldsymbol{B}_i = \left[\frac{1}{n^{(0)}}\sum\limits_{j\ne i} \left(\frac {\boldsymbol r_{ij}\otimes\boldsymbol r_{ij}}{|\boldsymbol r_{ij}|^2}\right) \omega _{ij} \right]^{-1}
$$

このベクトルの勾配から、次の一般的な関係式を考慮するとベクトルの発散を直ちに計算することができるよ。

$$
\nabla\cdot \boldsymbol{u}={\rm tr}(\nabla \boldsymbol{u})= (\nabla \boldsymbol{u}) :\boldsymbol{ I}
$$

「$${:}$$」はテンソル同士の複内積を表す記号だよ。この関係式を用いるとベクトルの発散は

$$
\langle\nabla\cdot \boldsymbol{u}\rangle_i= \left\{ \left[ \frac{1}{n^{(0)}}\sum\limits_{j\ne i}\,\frac{ \boldsymbol{u}_{ij} }{|\boldsymbol{r}_{ij} |} \otimes\frac{ \boldsymbol{r}_{ij}}{|\boldsymbol{r}_{ij} |}\,\omega_{ij} \right] \cdot \boldsymbol{B}_i\right\} : \boldsymbol{I}
$$

と表すことができるね。目的のスカラー関数のラプラシアンは勾配の発散で得られるので、

$$
\langle \nabla^2 \phi\rangle_i = \langle \nabla\cdot (\nabla \phi) \rangle_i =\left\{ \left[ \frac{1}{n^{(0)}}\sum\limits_{j\ne i}\,\frac{ (\nabla \phi)_{ij} }{|\boldsymbol{r}_{ij} |} \otimes\frac{ \boldsymbol{r}_{ij}}{|\boldsymbol{r}_{ij} |}\,\omega_{ij} \right] \cdot \boldsymbol{B}_i\right\} : \boldsymbol{I}
$$

と表すことができるね。$${(\nabla \phi)_{ij}}$$は$${i}$$番目の粒子を基準として周囲の粒子の位置における$${\nabla \phi }$$の値なので、先に導出したスカラー関数の勾配は和をとる前の

$$
\nabla \phi_i = \left[ \left(\frac {\boldsymbol r_{ij}\otimes\boldsymbol r_{ij}}{|\boldsymbol r_{ij}|^2}\right) \omega _{ij} \right]^{-1}\frac{ \phi_{ij}\boldsymbol{r}_{ij}}{|\boldsymbol{r}_{ij} |^2}\,\omega_{ij}
$$

が$${(\nabla \phi)_{ij}}$$を意味するので代入して整理していくね。

$$
\langle \nabla^2 \phi\rangle_i =\left\{ \left[ \frac{1}{n^{(0)}}\sum\limits_{j\ne i}\,\left[ \left(\frac {\boldsymbol r_{ij}\otimes\boldsymbol r_{ij}}{|\boldsymbol r_{ij}|^2}\right) \omega _{ij} \right]^{-1}\left[ \left(\frac {\boldsymbol r_{ij}\otimes\boldsymbol r_{ij}}{|\boldsymbol r_{ij}|^2}\right) \omega _{ij} \right]\frac{ \phi_{ij}}{|\boldsymbol{r}_{ij} |^2}\,\omega_{ij} \right] \cdot \boldsymbol{B}_i\right\} : \boldsymbol{I}\\
= \left[ \frac{1}{n^{(0)}}\sum\limits_{j\ne i}\frac{ \phi_{ij}}{|\boldsymbol{r}_{ij} |^2}\,\omega_{ij} \right] \boldsymbol{B}_i : \boldsymbol{I} = {\rm tr}\left( \boldsymbol{B}_i\right) \left[ \frac{1}{n^{(0)}}\sum\limits_{j\ne i}\frac{ \phi_{ij}}{|\boldsymbol{r}_{ij} |^2}\,\omega_{ij} \right]                                
$$

CG法ラプラシアンモデル

$$
\langle \nabla^2 \phi\rangle_i ={\rm tr}\left( \boldsymbol{B}_i\right) \left[ \frac{1}{n^{(0)}}\sum\limits_{j\ne i}\frac{ \phi_{ij}}{|\boldsymbol{r}_{ij} |^2}\,\omega_{ij} \right]
$$

なお粒子が等方的な場合、$${{\rm tr}\left( \boldsymbol{B}_i\right) = {\rm dim}^2}$$となるので、

$$
\langle \nabla^2 \phi\rangle_i =\frac{{\rm dim}^2}{n^{(0)}}\sum\limits_{j\ne i}\frac{ \phi_{ij}}{|\boldsymbol{r}_{ij} |^2}\,\omega_{ij}
$$

となるね。これは元のラプラシアンモデルと係数が異なるよ。

実装まとめ

MPS法のグラディエントモデルとラプラシアンモデルをそのままCG法に置き換えるよ。

１．仮速度を計算する際の速度ベクトル（ラプラシアンモデルの差し替え）

$$
\boldsymbol{v}^{(*)}_i = \boldsymbol{v}^{(k)}_i + \left[\nu \left\langle \nabla^2\boldsymbol{v}\, \right\rangle^{(k)}_i + \boldsymbol{g} \right] \Delta t
$$

速度ベクトルのラプラシアンを次のように置き換えるよ。

$$
\left\langle \nabla^2\boldsymbol{v}\, \right\rangle^{(k)}_i = \frac{2 \nu{\rm dim}}{ \lambda^{(0)} n^{(0)} } \sum\limits_{j=1,\, j\ne i}^N (\boldsymbol{v}_j^{(k)} - \boldsymbol{v}_i^{(k)})\, \omega(|\boldsymbol{r}_j^{(k)} - \boldsymbol{r}_i^{(k)}|)
$$

$$
↓
$$

$$
\left\langle \nabla^2\boldsymbol{v}\, \right\rangle^{(k)}_i = {\rm tr}\left( \boldsymbol{B}_i\right)\left[ \frac{\nu}{n^{(0)}}\sum\limits_{j\ne i}\frac{\boldsymbol{v}_{j}^{(k)} -\boldsymbol{v}_{i}^{(k)} }{|\boldsymbol{r}_j^{(k)} - \boldsymbol{r}_i^{(k)} |^2}\,\omega(|\boldsymbol{r}_j^{(k)} - \boldsymbol{r}_i^{(k)}|) \right]
$$

２．圧力ポアソン方程式の左辺と右辺（時刻：仮時刻）

$$
\left\langle \nabla^2 P\, \right\rangle^{(k+1)}_i = \frac{\rho_0}{ \Delta t} \,\left\langle\nabla\cdot\boldsymbol{v}\right\rangle^{(*)}_i
$$

左辺の圧力ラプラシアンを次のように置き換えるよ。

$$
\langle \nabla^2 P\, \rangle_i^{(k+1)} =\frac{2{\rm dim}}{\lambda^{(0)}n^{(0)}} \sum\limits_{j=1,\, j\ne i}^N (P^{(k+1)}_j-P^{(k+1)}_i)\, \omega(|\boldsymbol{r}^{(*)}_j - \boldsymbol{r}^{(*)}_i|)
$$

$$
↓
$$

$$
\langle \nabla^2 P\,\rangle_i^{(k+1)} = {\rm tr}\left( \boldsymbol{B}_i\right)\left[ \frac{1}{n^{(0)}}\sum\limits_{j=1,\, j\ne i}^N(P^{(k+1)}_j-P^{(k+1)}_i)\,\frac{ \omega(|\boldsymbol{r}^{(*)}_{j}-\boldsymbol{r}^{(*)}_{i}|) }{|\boldsymbol{r}^{(*)}_{j}-\boldsymbol{r}^{(*)}_{i} |^2} \right]
$$

続いて、右辺の速度ベクトルの発散を次のように置き換えるよ。

$$
\left\langle\nabla\cdot\boldsymbol{v}\right\rangle^{(*)}_i =\frac{1}{n^{(0)}} \sum\limits_{j\ne i} \frac{r_e}{|\boldsymbol{r}_j^{(*)}-\boldsymbol{r}^{(*)}_i|^3}\, (\boldsymbol{r}_j^{(*)}-\boldsymbol{r}^{(*)}_i) \cdot (\boldsymbol{v}_j^{(*)}-\boldsymbol{v}^{(*)}_i)
$$

$$
↓
$$

$$
\langle\nabla\cdot \boldsymbol{v}\rangle^{(*)}_i=\frac{1}{n^{(0)}} \sum\limits_{j\ne i}\,\frac{ \omega_{ij} }{|\boldsymbol{r}_{ij} |^2} \, \left[ \left( \boldsymbol{v}_{ij} \otimes \boldsymbol{r}_{ij} \right) \cdot \boldsymbol{B}_i\right] : \boldsymbol{I}
$$

３．時刻$${k+1}$$の速度を計算する際の圧力に関するグラディエントモデルの差し替え。

$$
\boldsymbol{v}^{(k+1)}_i = \boldsymbol{v}^{(*)}_i -\frac{1}{\rho_0} \left\langle \nabla P\, \right\rangle^{(k+1)}_i \Delta t
$$

圧力のグラディエントを次のように置き換えるよ。

$$
\langle \nabla P\, \rangle_i^{(k+1)} =\frac{{\rm dim}}{n^{(0)}} \sum\limits_{j=1,\, j\ne i}^N \frac{ P^{(k+1)}_j-P^{(k+1)}_i}{|\boldsymbol{r}^{(*)}_j - \boldsymbol{r}^{(*)}_i|^2}\,(\boldsymbol{r}^{(*)}_j - \boldsymbol{r}^{(*)}_i)\, \omega(|\boldsymbol{r}^{(*)}_j - \boldsymbol{r}^{(*)}_i|)
$$

$$
↓
$$

$$
\langle \nabla P\, \rangle_i^{(k+1)} =\frac{{\rm dim}}{n^{(0)}} \sum\limits_{j=1,\, j\ne i}^N \frac{ P^{(k+1)}_j-P^{(k+1)}_i}{|\boldsymbol{r}^{(*)}_j - \boldsymbol{r}^{(*)}_i|^2} \,\boldsymbol{C}^{(*)}_{i} \cdot (\boldsymbol{r}^{(*)}_j - \boldsymbol{r}^{(*)}_i)\, \omega(|\boldsymbol{r}^{(*)}_j - \boldsymbol{r}^{(*)}_i|)
$$

$$
\boldsymbol{C}^{(*)}_{i}= \frac{1}{{\rm dim}}\left[ \frac{1}{n^{(0)}}\sum\limits_{j=1,\, j\ne i}^N\frac {\omega(|\boldsymbol{r}^{(*)}_j - \boldsymbol{r}^{(*)}_i|)}{|\boldsymbol{r}^{(*)}_j - \boldsymbol{r}^{(*)}_i|^2}\left( \boldsymbol r^{(*)}_{ij}\otimes\boldsymbol r^{(*)}_{ij} \right) \right]^{-1}
$$

今回はアルゴリズムの導出のみだよ。次回、Pythonによる実装と検証を行っていくね。これからも応援よろしくお願いしまーす。