aokomoriuta's blog

青子守歌のブログ

ラックスの等価定理と粒子法の勾配微分演算子モデル

この記事は数値計算 Advent Calendar 2018の18日目の記事です。

はじめに

ご存知の方もいらっしゃると嬉しいのですが、私はOpenMPSというオープンソースの粒子法実装を開発・公開しています。

みなさん、粒子法ってご存知ですか？この記事を読んでいる人で知らない方はいないと思いますが、粒子法のよくある誤解として「でも格子法に比べて精度がでないでしょ？」「数学的な誤差評価ができないんでしょ？」というのがあります。もちろん、研究者の数が少ないという意味で、格子法に比べてまだまだ解析・改善できていないところが多いというのは否定しません。また特に、半陰解法であるMPS法については、越塚先生らが最初に提案した原始MPS法（標準MPS法とも呼ばれる）では圧力擾乱がかなり最悪だったことからこの手の誤解が未だに存在するというのは悲しい事実ではあります。

しかし、現代においてはかなり粒子法についての研究が進み、後藤先生・カイヤー先生らの高精度粒子法と呼ばれる改良スキーム群を始め、改善手法が多く開発されています。これら近年の粒子法の改善については、Gotoh, H. & Khayyer, A. (2016): "Current achievements and future perspectives for projection-based particle methods with applications in ocean engineering", J. Ocean Eng. Mar. Energy, 2:251がとても詳しいですので、（特に水などの非圧縮性流れに対する）粒子法を使う・研究する方は、必読です（最近フリーアクセスになったようです）。OpenMPSでもかなり参考にさせていただきました。最近発売された後藤先生の教科書『粒子法』にも詳しく載っていますので、「英語はちょっと・・・」という方はそちらも合わせてご覧ください。

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

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

本題

さて、前置きが少し長くなってしまいましたが、この記事では、その改善手法の代表として、私が特に重要（効果が顕著）だと経験的に思っている、勾配修正(GC; Gradient Correction)法を紹介したいと思います。この手法は2011年に発表された論文Khayyer A, Gotoh H (2011): "Enhancement of stability and accuracy of the moving particle semi-implicit method", J Comput Phys, 230:3093–3118.によるものです。ここでは厳密さはあまり求めず、私の理解の元で整理し考え方を紹介することに注力しますので、ちゃんとした厳密な証明は元論文等をご覧ください。

ラックスの等価定理

まず勾配修正法の話を始める前に、重要な話として、ラックスの等価定理はご存知でしょうか？数値計算・計算力学のアドベントカレンダーをご覧になるような強い方々はもちろんご存知とは思いますが、この後の話に深く関わってくるので軽くおさらいしておくと、安定かつ適合ならば、その時に限り解は収束するというものです。つまり

  • 安定性：演算によって混入した誤差（丸め誤差、打切り誤差など。離散化誤差は含まない）が、時間発展によって増大しない
  • 適合性：刻み幅を0極限した時に、離散化方程式自体と元の微分方程式が一致する
  • 収束性：刻み幅を0極限した時に、離散化方程式の解と元の微分方程式の解が一致する

というもので、離散化方程式の計算で得られた解の精度を評価するための根本的な定理です。ここで大雑把に言えば、安定性は時間方向の離散化、適合性は空間方向の離散化に対する制約です。今回の勾配修正法は、その名の通り「（ベクトルの微分演算子である）勾配」に関わるものですので、この適合性についての評価になります。

粒子法におけるテイラー展開

さて、適合性を評価する一番簡単な方法は、離散化する時に使ったテイラー展開と比較することです。ということで、ある物理量\phiについて、粒子iの近傍にある粒子jの値をテイラー展開で推定しましょう。その式は以下のようになりますね。

\displaystyle \phi_j - \phi_i = \boldsymbol{\nabla} \phi_i \cdot \boldsymbol{r}_{ij}+ O\left( h^2 \right)

ここで、\boldsymbol{r}_{ij} = \boldsymbol{r}_j - \boldsymbol{r}_iは2粒子間の相対位置ベクトルです。とりあえずベクトルの長さが違うと議論が面倒なので、長さ r_{ij} = \left| \boldsymbol{r}_{ij} \right|で割って正規化（単位ベクトル化）しましょう。ついでに2次以降を無視し、1次精度になるようなスキームを取ります。

\displaystyle \boldsymbol{\nabla} \phi_i \cdot \frac{\boldsymbol{r}_{ij}}{r_{ij}} = \frac{\phi_j - \phi_i}{r_{ij}}

さて、通常のよくあるテイラー展開と違ってこの単位相対位置ベクトルが邪魔で、このままでは\boldsymbol{\nabla} \phiが推定できません。そこで、粒子法では、積分補間子と呼ばれる量で重み付けした値を積分すると、「この相対位置ベクトルが消せる」ことにするというモデルを取ります。積分補間子はMPS法とSPH法で違いますが（むしろその違いがこの両者の違いの本質）、ここでは元論文に合わせてMPS法の流儀として「大きさが重み関数の相対位置ベクトル」すなわち w_{ij} \frac{d\boldsymbol{r}}{r_{ij}}とします。ここでw_{ij}はi-j粒子間の重みを表します。SPHの場合は「粒子の体積とカーネル関数の微分」になるだけで似たような議論ができますのでSPHが好みな人は変わりに \boldsymbol{\nabla} w_{ij} dVを使って計算ください。

さて、「積分補間子を使って積分」 \int w_{ij} \frac{d\boldsymbol{r}}{r_{ij}}しましょう。左辺は先述の通り、この積分補間子による積分で単位相対位置ベクトルが消えますが、重み関数の積分の分が増えているので、

\displaystyle \boldsymbol{\nabla} \phi_i \sum_{j \neq i} w_{ij} = \boldsymbol{\nabla} \phi_i n_0

になります。最後は厳密にはn_iですが、非圧縮性流れであれば基準粒子数密度n_0を使っても問題ないはずです。では右辺はというと

\displaystyle \int \frac{\phi_j - \phi_i}{r_{ij}} \frac{d\boldsymbol{r}}{r_{ij}} w_{ij}

となりますね。これを普通に離散化すると、つまり \intが総和になり、d\boldsymbol{r}=\boldsymbol{r}_{ij}となるので、最終的に

\displaystyle \boldsymbol{\nabla} \phi_i = \frac{1}{n_0} \sum_{j \neq i} \frac{\phi_j - \phi_i}{r_{ij}} \frac{\boldsymbol{r}_{ij}}{r_{ij}} w_{ij}

というモデルができあがります。これが原始MPS法の勾配モデルの正体です（原論文では違う説明がされていますが）。

勾配修正法

さて、気づいた方もいらっしゃると思いますが、これの何が問題かと言うと、このモデルは本当に1次精度なのか？ということです。せっかくテイラー展開を使って1次精度から出発したのに、途中で「左辺は積分補間子によって消せることにする」という別のモデルを導入したりしてしまっていて、ずれてしまっていますよね。

なので、ちゃんと左辺も積分補間子で評価してみましょう。愚直に代入すれば

\displaystyle \int \boldsymbol{\nabla} \phi_i \cdot \frac{\boldsymbol{r}_{ij}}{r_{ij}} w_{ij} \frac{d\boldsymbol{r}}{r_{ij}}

となりますが、\boldsymbol{\nabla} \phi_iはdrに関係ないことと、テンソル代数の基礎によって \left(\boldsymbol{a} \cdot \boldsymbol{b} \right) \boldsymbol{c} = \left(\boldsymbol{c} \otimes \boldsymbol{b} \right) \boldsymbol{a}と交換できることを思いだせば、

\displaystyle \boldsymbol{\nabla} \phi_i \int \frac{\boldsymbol{r}_{ij}}{r_{ij}} \otimes \frac{d\boldsymbol{r}}{r_{ij}} w_{ij}

と変形できることが容易にわかります。はい、全然消えていません。誰ですか、 \boldsymbol{\nabla} \phi_i n_0とみなすモデルだとか言った人？

ということで、勾配修正法では、この右側の逆行列すなわち

\displaystyle C_i = \left( \int \frac{\boldsymbol{r}_{ij}}{r_{ij}} \otimes \frac{d\boldsymbol{r}}{r_{ij}} w_{ij} \right)^{-1}

という修正行列を用いて、従来モデルの勾配モデルにこの修正行列を乗じた

\displaystyle \boldsymbol{\nabla} \phi_i = \sum_{j \neq i} \frac{\phi_j - \phi_i}{r_{ij}} \frac{C_i \boldsymbol{r}_{ij}}{r_{ij}} w_{ij}

というモデルを用います。これで1次精度が保証されますね！この修正行列は粒子jから受ける勾配の方向が修正されるという物理的意味を表しています。

では逆に、このような修正が不要、すなわち従来モデルが成立する時というのはどんな時でしょうか？それはすなわち、修正行列（の逆行列）が不要だった時ですから、

\displaystyle {C_i}^{-1} = \int \frac{\boldsymbol{r}_{ij}}{r_{ij}} \otimes \frac{d\boldsymbol{r}}{r_{ij}} w_{ij} = n_0 I

であればよいですね。ここでIは単位行列です。これが成立する時とはどんな時でしょうか？これは実は近傍粒子が均等に配置されている場合になります。これを確かめるには実際に分解して検算してみるのがいいですが、例えば、xおよびyの対角項は等しく

\displaystyle {{C_i}^{-1}}_{x, x} = {{C_i}^{-1}}_{y, y}

である、つまり

\displaystyle \int_X \frac{x dx}{r^2} w = \int_Y \frac{y dy}{r^2} w

となってほしいわけですが、これが成立するのは、それぞれの積分区間であるX, Yの範囲が同等、つまり、離散化状態で近傍粒子の分布がどの方向にも均等であることを意味します。

一方で、非対角項を見れば、

\displaystyle {{C_i}^{-1}}_{x, y} = {{C_i}^{-1}}_{y, x} = 0

つまり、

\displaystyle \int_Y \frac{x dy}{r^2} w = \int_X \frac{y dx}{r^2} w = 0

ですが、これが成立するためには、x,yはそれぞれ直交するY, X方向の範囲で偶関数であること、つまり、離散化状態で近傍粒子の正・負の分布が均等であるべきということです。

まとめ

以上が、（MPS法を例にした）粒子法の高精度化に関する定式化の一例でした。最初に紹介した論文や教科書には他にもまだまだたくさんの方法が解説されています。ぜひ勉強してみて、実際のコードに導入してみてください。