出身研究室の同期が書いた論文の査読依頼がきた
論文の査読依頼が来た場合、よっぽど専門と外れた内容の論文以外は引き受けるようにしています。存じ上げなかったオープンアクセス誌であるCFD lettersから査読依頼が来て、タイトルだけ見て、まあいいかと思って軽い気持ちで引き受けました。
学生時代に同じ研究室にいた留学生がcorresponding authorでした。
彼も頑張ってるんだなあとほほえましい気持ちになりました。
査読は断りました。近すぎるし。
論文誌のウェブサイト上ですでに承諾してしまったため、editorに下記のメールを直接送ってみました。
Dear Prof. **,
Thank you for giving me an opportunity to contribute to the journal "**". Although I accepted to review this manuscript, I am afraid that I may not be an appropriate reviewer. The corresponding author and I belonged to the same lab. when we were PhD students, and then have a personal relationship. Please would you reconsider more appropriate reviewers for the corresponding author? Thank you.
Sincerely,
30分以内に返信が来て、他の査読者に回してくれることになりました。面倒な奴だなと思われたかもしれませんが、結果としては誠実に対応頂いたように思います。
真夏の珍事でした。
微分方程式の時間発展を有限要素法で解く(不連続Galerkin法)
はじめに
偏微分方程式の数値解法として有限要素法は一般的な方法であり、多種多様な物理現象の計算に応用されています。
ただし、有限要素法が適用されるのは空間微分項に対してであり、時間微分項、すなわち時間発展(もしくは積分)の方法については有限差分法を用いることが(自分が知る限り)一般的です。複数の文献でも、時間積分については差分を導入することが特に注釈なく当たり前のように記述されています(例外もあるのですが、長くなるので後述します)。
有限要素法による時間積分の手法自体は30-40年前にはすでに提案され、数学者がまとめた有限要素法の書籍に短くはありますが記載されています1,2。当時の計算技術では、この手法の高い計算コストが致命的な問題となり、多数派は差分法を採用した結果、今日の状況になったのかもしれません。
そこで、軽い気持ちで有限要素法での時間積分の方法について調べてみることにしました。有限要素法の初歩的な内容がわかっていれば理解できる程度までの記述とし、数学的に怪しい表現をしているので、気になる方は下に載せた文献を参照してください。
弱形式の導出
まず、適当な例題を用意して、弱形式化された方程式を導出する。空間微分項の扱いは通常の有限要素法の作法で扱えるので、本記事では時間微分項の扱いのみを考慮する。その場合、時間に関する常微分方程式の解法を考えればよい。
問題設定
常微分方程式の典型的な例として、下記のスカラーの時間に対する一次遅れ系を考える。
ここで便宜的に適当な時間の最大値を用意した。については、ここではざっくりと正定値であることのみ定める。
重み付き残差法
有限要素法の定石通り、まずは重み付き残差法を用いて上の式を弱形式で表現する。重み関数をとし、
が成り立つとする。このとき積分内の第一項は
となり、基本境界条件が設定できる。以上をまとめると下記のようになる。
離散化
有限要素法の通常の手順と同様、積分区間を有限の時間間隔で分割することを考え、分割区間内で変数を多項式近似する。このときに同じ近似関数を用いれば通常のGalerkin法となるが、全ての時間刻みの変数の値を自由度に持つことになり、本来は空間に対しても自由度を与えたいことを考えれば、自由度が膨大になりすぎて明らかに現実的ではない。 解が時間に対して周期解になっていれば別だが、その場合は解を周波数空間に変換して解くほうが効率的であり、スペクトル要素法として整備されている。
不連続Galerkin法
上の問題を回避するための方法として、時間刻みの各点において、前後の多項式が連続でないことを許容する(不連続Galerkin法)。
適当な時間, を考え,として時間刻みを設定する。ここでともに刻み幅の内部では連続とし、端点のみ不連続性を許容する。
すなわち、とを用意し、下記の図のように定義する。
要するに
を定義しておく。なお、Galerkin法の定義通り、についてもと同様の不連続性を与える。このとき、式(I)の左辺は
を導入すると
と表せる。式(I)に上の式を代入して整理すると、を除いて、各時間刻みで独立であることがわかる。そこで、適当な区間]について抜き出してみると、
となる。最初の積分系の式の一部の積分区間のみを抜き出し、隣接区間との間の非独立成分を付け足した形になっている。隣接区間として右端を考慮していないが、未来(右端)から過去は変えられないと解釈すると、直感的に理解できるかもしれない。
多項式近似
最後に、について区間で多項式近似する。多項式の形はいろいろ考えられるが、後の計算の都合上、下記の次の関数として表現する。
ここでは節点に割り振られた値であり、括弧の項は通常の有限要素法ならば形状関数と呼ぶ(時間の多項式でも"形状"関数なのだろうか?)。0次近似ならば
1次近似ならば
である.なおについても同じ形の近似関数を定義する.
通常の有限要素法の手順ではあまりやらないが、計算の簡略化のため、先にのみ展開してみる。に注意して、1次の場合について書き下すと
となる。これが任意ので成立する条件は下記の通り。
なお、0次の場合は式(III.1)のみが成立する。2次以上を考える場合も同様の手順でいけそう。
離散式
0次近似、1次近似の場合について,離散式を導出してみる。
まず0次近似について、式(III.1)に上の多項式を代入して整理する。積分区間でのの時間微分が0であることに注意すると,下記のようになる。
左辺は陰的Euler法の形、右辺は区間内平均の形になっている。
1次近似の場合について、同様に式(III.1),(III.2)に多項式を代入すると、
となる。外力項は各自の問題設定に従って連続関数として与えればいいので省略する。 以上で、最終的に時間発展に必要な離散式(代数方程式)を実装レベルまで書き下せた。
解析解との比較
有限要素法による時間積分を考え、積分区間内での多項式を0次または1次とした場合について検討した。この表現がどの程度適切か、解析解との比較によって考えてみる。とすれば、もとの常微分方程式の解析解は、漸化式の形で
ここでである。この形で式(IV)を整理すると、
式(V)について
となる(この式はのPadé近似らしい[1])(未確認)。
せっかくなので差分法の場合も考える。 Euler陽解法の場合、
だから、となる。
次に、Crank-Nicolson法を用いると、
だから、
となる。
近似関数のプロット
解析解とそれぞれの解法で得られたを描画してみると下記の通り。 Euler陽解法が解析解から大きく外れていくのはご愛敬。明らかに1次多項式の不連続Galerkin法の結果が最も解析解に近くなった。0次多項式での不連続Galerkin法と等価なEuler陰解法だと、予想通り数値粘性が強く出てきて、解析解より値の減少が少ない。 Crank-Nicolson法は安定かつ実装がシンプル、それでいて二次精度と、使い勝手がよい印象だったが、1次の不連続Galerkin法には勝てなかった。このことはから予想がついていて、各についてとすると、解析解、陰的Euler, 不連続Galerkin法は0に漸近するのに対して、Crank-Nicolson法だと-1に収束する。熱流体分野ではよく使われている印象だったが、気を抜いて刻み幅を大きくすると発散方向に解の誤差が大きくなるということか。
まとめ
有限要素法による時間積分について調べてみました。1次の不連続Galerkin法でさえ、計算の自由度は差分による陰解法の2倍となりますが、低次の差分解法と比較して精度向上が見込めそうです。とはいえ、理屈はかなり明快なので、それなりに使われていてもよさそうなものですが…。
この方法の自然な拡張の一つに有限要素法による流体構造連成計算手法の代表格ともいえるspace-time法(c.f., 3)があり、こちらに研究の動向が移り、わざわざ時間項のみの扱いに着目する研究者や技術者がいなくなったのかもしれません。
数十年前に確立された手法を興味本位で調べて時間を費やしたことに少し後悔していなくもないですが、誰かの興味や好奇心に少しでも刺されば幸いです。
- 作者:Thomee, Vidar
- 発売日: 2006/09/15
- メディア: ハードカバー
- 作者:Johnson, Claes
- 発売日: 2009/01/15
- メディア: ペーパーバック
-
V. Thomée, Galerkin Finite Element Methods for Paraboloc Problems, Springer, 1984. ↩
-
C. Johnson, Numerical solutions of partial differential equations by the finite element method, Dover, 2009↩
-
Y. Bazilevs, T. Takizawa, and T. E. Tezduyar. Computational fluid‐structure interaction: methods and applications. Wiley, 2013↩
VS codeに移行
Sublime text3をここ数年くらい使用してきたけど、日本語環境が使用できなくなったり、
Latex環境構築があまり楽でなかったりして、だんだん微妙に感じてきたのでVS codeに徐々に移行していきます。
インストールして、特に何か難しい操作もなく直感的にいろいろ遊べそう。ただLatex設定は少し手間取った。
下のまとめで全て解決したけど。
qiita.com
Latex環境としてOverleaf v2を使ってきたけど、特に人と共有して作業することもないので、
オフラインでこっちで作業することにしても問題なさそう。
microsoft社だけあってgitとの相性もいいし、早々にこっちに本格的に移行しよう…
netcdfのインストール(C,C++,Fortran)
netcdfは、日本のスパコンには標準装備でインストールされているようだ。汎用的なデータ保存用フォーマットとして便利そうなので、手元の環境にインストールした。とりあえずconfigureのオプションだけまとめる。C版が基本になっており、C++とFortran版はC板のラッピングとなっている。Fortranのmake checkが中々通らず苦労した。
./configure --prefix=/hogehoge/netcdf-4.4.1.1 CC=icc LDFLAGS=-L/hogehoge/hdf5-1.8.18/lib CPPFLAGS="-I/hogehoge/hdf5-1.8.18/include -DNDEBUG -DpgiFortran" CXX=icpc FC=ifort --enable-netcdf-4
./configure --prefix=/hogehoge/netcdf-cxx4-4.3.0 CC=icc CXX=icpc LDFLAGS=-L/hogehoge/netcdf-4.4.1.1/lib CPPFLAGS=-I/hogehoge/netcdf-4.4.1.1/include
./configure --prefix=/hogehoge/libs/netcdf-fortran-4.4.4 CC=icc FC=ifort LDFLAGS="-L/hogehoge/libs/netcdf-4.4.1.1/lib -lnetcdf" CPPFLAGS="-I/hogehoge/libs/netcdf-4.4.1.1/include" --disable-fortran-type-check
何回かコンパイルに失敗し、ダメ元でフォルダを一度消して、tarファイルの解凍からやり直したらうまくいった。distcleanなどは行っていたのだが、何やらゴミが残っていたらしい。
Ubuntu16.04へのVTK7.1のインストール
言わずと知れた可視化ライブラリVTKであるが、その機能はフリーソフトparaviewにより簡単にGUIで使用することができる。物好きでもない限りわざわざVTKを触る必要もない…と思っていた。ただ、paraviewを使用した単純作業を最低でも数百回やることになったら話は別である。最初はparaviewのpythonスクリプトでのり切ろうかと思ったが、どうもいろいろ制約が生じて面倒な予感。保険のつもりでVTKの最新版をインストールした。
面倒なので、やった作業だけまとめておく。
libncurses-devのインストール(apt-getでよい)
cmakeの最新版のインストール(ソースコードから)
端末を再起動させ、cmakeとccmakeのバージョンが一致している(3.8.0以上)であることを確認する。
Qt5のインストール(optional)
vtkの最新版をソースからインストール(ccmakeによりインテルコンパイラを選択)
bashrcにVTK_DIRを記入
Intel Math Kernel Libraryでの大規模疎行列の直接解法<実践>
MKLは普通にLAPACKが使えるのだが、なぜか直接解法ソルバであるPARDISOがデフォルトで使用できる。LAPACKの普通のLU分解とは違い、PARDISOはマルチフロンタル法による高速化、グラフ理論に基づくfill inの削減(metisを使用している?)が行われている。数十万次元くらいの定常解析なら十分実用可能ではないかと期待している。
PARDISOは下記の公式ページにマニュアルがおいてある。サンプルコードを見ればサブルーチンpardisoの使い方がわかるはず。
PARDISO 5.0.0 Solver Project
サブルーチンの引数などは下記をみればだいたいわかるので、正直余裕すぎて書くことがない…と、思っていたが、なぜかプログラムが動かない。サンプルコードは動くが、それをモジュール化して自分のコードに組み込むとセグメンテーション違反が生じる。
仕方ないので調べたところ、/opt/intel/mkl/exampleにpardisoのサンプルコードが入っているらしい。
PARDISO 使用時のヒント | iSUS
確認したところ、公式のサンプルコードと内容が微妙に違う!頭を抱えつつMKLのサンプルを自分のコードにコピペして実行すると、正常に動作することを確認。理不尽。intel MKLのバージョンは2016だけど、PARDISOのバージョンは多少古いのかもしれない。