2020-08-18

出身研究室の同期が書いた論文の査読依頼がきた

雑記

論文の査読依頼が来た場合、よっぽど専門と外れた内容の論文以外は引き受けるようにしています。存じ上げなかったオープンアクセス誌である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分以内に返信が来て、他の査読者に回してくれることになりました。面倒な奴だなと思われたかもしれませんが、結果としては誠実に対応頂いたように思います。
真夏の珍事でした。

2020-08-17

微分方程式の時間発展を有限要素法で解く（不連続Galerkin法）

fem

はじめに

偏微分方程式の数値解法として有限要素法は一般的な方法であり、多種多様な物理現象の計算に応用されています。

ただし、有限要素法が適用されるのは空間微分項に対してであり、時間微分項、すなわち時間発展（もしくは積分）の方法については有限差分法を用いることが（自分が知る限り）一般的です。複数の文献でも、時間積分については差分を導入することが特に注釈なく当たり前のように記述されています（例外もあるのですが、長くなるので後述します）。

有限要素法による時間積分の手法自体は30-40年前にはすでに提案され、数学者がまとめた有限要素法の書籍に短くはありますが記載されています¹,²。当時の計算技術では、この手法の高い計算コストが致命的な問題となり、多数派は差分法を採用した結果、今日の状況になったのかもしれません。

そこで、軽い気持ちで有限要素法での時間積分の方法について調べてみることにしました。有限要素法の初歩的な内容がわかっていれば理解できる程度までの記述とし、数学的に怪しい表現をしているので、気になる方は下に載せた文献を参照してください。

　弱形式の導出

まず、適当な例題を用意して、弱形式化された方程式を導出する。空間微分項の扱いは通常の有限要素法の作法で扱えるので、本記事では時間微分項の扱いのみを考慮する。その場合、時間に関する常微分方程式の解法を考えればよい。

問題設定

常微分方程式の典型的な例として、下記のスカラー $u$ の時間 $t$ に対する一次遅れ系を考える。

$\frac{\partial u}{\partial t}+Au=f, \text{ for } t\in[0,T], \text{ with } u(0)=v.$

ここで便宜的に適当な時間の最大値 $T$ を用意した。 $A$ については、ここではざっくりと正定値であることのみ定める。

　重み付き残差法

有限要素法の定石通り、まずは重み付き残差法を用いて上の式を弱形式で表現する。重み関数を $w, w(T)=0$ とし、

$\int_{0}^{T}{\left(\frac{\partial u}{\partial t}w+Auw-fw\right)}dt=0,$

が成り立つとする。このとき積分内の第一項は

$\int_{0}^T{\frac{\partial u}{\partial t}w dt}=\int_0^T{\left[\frac{\partial uw}{\partial t}-\frac{\partial w}{\partial t}u \right]dt}, \\=-[uw]_0^T-\int_0^T{\frac{\partial w}{\partial t}u dt}=-vw(0)-\int_0^T{\frac{\partial w}{\partial t}u dt}$

となり、基本境界条件が設定できる。以上をまとめると下記のようになる。

$\int_{0}^{T}{\left(-\frac{\partial w}{\partial t}u+Auw \right)}dt=vw(0)+\int_{0}^{T}{fw}dt.　\tag{I}$

離散化

有限要素法の通常の手順と同様、積分区間を有限の時間間隔で分割することを考え、分割区間内で変数 $u,w$ を多項式近似する。このとき $u,w$ に同じ近似関数を用いれば通常のGalerkin法となるが、全ての時間刻みの変数の値を自由度に持つことになり、本来は空間に対しても自由度を与えたいことを考えれば、自由度が膨大になりすぎて明らかに現実的ではない。解が時間に対して周期解になっていれば別だが、その場合は解を周波数空間に変換して解くほうが効率的であり、スペクトル要素法として整備されている。

不連続Galerkin法

上の問題を回避するための方法として、時間刻みの各点 $t_n$ において、前後の多項式が連続でないことを許容する（不連続Galerkin法）。

適当な時間 $t_n$ , $t_{n-1}$ を考え， $t_n-t_{n-1}=k_n$ として時間刻みを設定する。ここで $u,w$ ともに刻み幅の内部では連続とし、端点のみ不連続性を許容する。

すなわち、 $u^{n-1}$ と $u_+^{n-1}$ を用意し、下記の図のように定義する。

f:id:bigvalley:20200816215819p:plain

要するに

$\lim_{s \to 0+} u(t_{n-1}+s)=u_+^{n-1}$

を定義しておく。なお、Galerkin法の定義通り、 $w$ についても $u$ と同様の不連続性を与える。このとき、式(I)の左辺は

$\int_{0}^{T}{-\frac{\partial w}{\partial t}u}dt=\sum_{n=1}^{N-1}\left[-[uw]_{t_{n-1+0}}^{t_n}+\int_{t_{n-1}}^{t_n}{\frac{\partial u}{\partial t}w}dt \right],$

積分区間の間の不連続性を考慮し、

$[u]^n=u_{+}^{n}-u^n$

を導入すると

$\text{rhs}=\int_{0}^{T}{\frac{\partial u}{\partial t}w}dt+\sum_{n=1}^{N-1}\left([u]^nw \right)+(u_+^0 w^0),$

と表せる。式(I)に上の式を代入して整理すると、 $u$ $^n$ を除いて、各時間刻みで独立であることがわかる。そこで、適当な区間 $I_n=(t_{n-1},t_n$ ]について抜き出してみると、

$\int_{I_n}{\left(\frac{\partial u}{\partial t}w+Auw\right)}dt+\left([u]^{n-1} w_+^{n-1}\right)=\int_{I_n}{fw}dt,\tag{II}$

となる。最初の積分系の式の一部の積分区間のみを抜き出し、隣接区間との間の非独立成分を付け足した形になっている。隣接区間として右端を考慮していないが、未来（右端）から過去は変えられないと解釈すると、直感的に理解できるかもしれない。

多項式近似

最後に、 $u,w$ について区間 $I_n$ で多項式近似する。多項式の形はいろいろ考えられるが、後の計算の都合上、下記の $q-1$ 次の関数として表現する。

$u(t)=\sum_{j=0}^{q-1}\tilde{u}_j\left(\frac{t-t_{n-1}}{k_n}\right)^j.$

ここで $\tilde{u}_j$ は節点 $j$ に割り振られた値であり、括弧の項は通常の有限要素法ならば形状関数と呼ぶ（時間の多項式でも"形状"関数なのだろうか？）。0次近似ならば

$u(t)=\tilde{u}_0=u_+^{n-1}=u^n,$

1次近似ならば

$u(t)=\tilde{u}_0+\tilde{u}_1\left(\frac{t-t_{n-1}}{k_n}\right), u_+^{n-1}=\tilde{u}_0, \text{ }u^n=\tilde{u}_0+\tilde{u}_1.$

である．なお $w$ についても同じ形の近似関数を定義する．

通常の有限要素法の手順ではあまりやらないが、計算の簡略化のため、先に $w$ のみ展開してみる。 $w_+^{n-1}=\tilde{w}_0$ に注意して、1次の場合について書き下すと

$\int_{I_n}{\left(\frac{\partial u}{\partial t}+Au\right)\left(\tilde{w}_0+\tilde{w}_1\frac{t-t_{n-1}}{k_n}\right)}dt+(u_+^{n-1}-u^{n-1})\tilde{w}_0\\ =\int_{I_n}{f\left(\tilde{w}_0+\tilde{w}_1\frac{t-t_{n-1}}{k_n}\right)}dt,$

となる。これが任意の $\tilde{w}_0,\tilde{w}_1$ で成立する条件は下記の通り。

$\int_{I_n}{\left(\frac{\partial u}{\partial t}+Au\right)}dt+(u_+^{n-1}-u^{n-1})=\int_{I_n}{f}dt \tag{III.1}$

$\int_{I_n}{\left(\frac{\partial u}{\partial t}+Au\right)\left(\frac{t-t_{n-1}}{k_n}\right)}dt=\int_{I_n}{f\left(\frac{t-t_{n-1}}{k_n}\right)}dt, \tag{III.2}$

なお、0次の場合は式(III.1)のみが成立する。2次以上を考える場合も同様の手順でいけそう。

離散式

0次近似、1次近似の場合について，離散式を導出してみる。

まず0次近似について、式(III.1)に上の多項式を代入して整理する。積分区間での $u$ の時間微分が0であることに注意すると，下記のようになる。

$Au^nk_n+(u^n-u^{n-1})=\int_{I_n}{f}dt.　 \tag{IV}$

左辺は陰的Euler法の形、右辺は区間内平均の形になっている。

1次近似の場合について、同様に式(III.1),(III.2)に多項式を代入すると、

$\tilde{u}_1+\tilde{u}_0+k_nA\tilde{u}_1+\frac{1}{2}k_nA\tilde{u}_1=u^{n-1}+\int_{I_n}{f}dt \tag{V.1}$

$\frac{1}{2}k_nA\tilde{u}_0+\frac{1}{2}k_n\tilde{u}_1+\frac{1}{3}k_nA\tilde{u}_1=\frac{1}{k_n}\int_{I_n}{(t-t_{n-1})f}dt, \tag{V.2}$

となる。外力項 $f$ は各自の問題設定に従って連続関数として与えればいいので省略する。以上で、最終的に時間発展に必要な離散式（代数方程式）を実装レベルまで書き下せた。

解析解との比較

有限要素法による時間積分を考え、積分区間内での多項式を0次または1次とした場合について検討した。この表現がどの程度適切か、解析解との比較によって考えてみる。 $f=0$ とすれば、もとの常微分方程式の解析解は、漸化式の形で

$u(t_n)=r(Ak_n)u(t_{n-1})$

ここで $r(x)=\exp(-x)$ である。この形で式(IV)を整理すると、

$r_0(x)=\frac{1}{1+x},$

式(V)について

$r_1(x)=\frac{1-\frac{x}{3}}{1+\frac{2}{3}x+\frac{1}{6}x^2},$

となる（この式は $\exp(-x)$ のPadé近似らしい[1])（未確認）。

せっかくなので差分法の場合も考える。 Euler陽解法の場合、

$\frac{u^n-u^{n-1}}{k_n}+Au^{n-1}=0$

だから、 $r_\text{Ex}(x)=(1-x)$ となる。

次に、Crank-Nicolson法を用いると、

$\frac{u^n-u^{n-1}}{k_n}+\frac{1}{2}A(u^n+u^{n-1})=0$

だから、

$r_c(x)=\frac{1-\frac{x}{2}}{1+\frac{x}{2}}$

となる。

近似関数のプロット

解析解とそれぞれの解法で得られた $r(x)$ を描画してみると下記の通り。 f:id:bigvalley:20200817093444p:plain Euler陽解法が解析解から大きく外れていくのはご愛敬。明らかに1次多項式の不連続Galerkin法の結果が最も解析解に近くなった。0次多項式での不連続Galerkin法と等価なEuler陰解法だと、予想通り数値粘性が強く出てきて、解析解より値の減少が少ない。 Crank-Nicolson法は安定かつ実装がシンプル、それでいて二次精度と、使い勝手がよい印象だったが、1次の不連続Galerkin法には勝てなかった。このことは $r(x)$ から予想がついていて、各 $r(x)$ について $x\rightarrow\infty$ とすると、解析解、陰的Euler, 不連続Galerkin法は0に漸近するのに対して、Crank-Nicolson法だと-1に収束する。熱流体分野ではよく使われている印象だったが、気を抜いて刻み幅を大きくすると発散方向に解の誤差が大きくなるということか。

まとめ

有限要素法による時間積分について調べてみました。1次の不連続Galerkin法でさえ、計算の自由度は差分による陰解法の2倍となりますが、低次の差分解法と比較して精度向上が見込めそうです。とはいえ、理屈はかなり明快なので、それなりに使われていてもよさそうなものですが…。

この方法の自然な拡張の一つに有限要素法による流体構造連成計算手法の代表格ともいえるspace-time法(c.f., ³)があり、こちらに研究の動向が移り、わざわざ時間項のみの扱いに着目する研究者や技術者がいなくなったのかもしれません。

数十年前に確立された手法を興味本位で調べて時間を費やしたことに少し後悔していなくもないですが、誰かの興味や好奇心に少しでも刺されば幸いです。

Galerkin Finite Element Methods for Parabolic Problems (Springer Series in Computational Mathematics (25))

作者:Thomee, Vidar
発売日: 2006/09/15
メディア: ハードカバー

Numerical Solution of Partial Differential Equations by the Finite Element Method (Dover Books on Mathematics)

作者:Johnson, Claes
発売日: 2009/01/15
メディア: ペーパーバック

Computational Fluid-Structure Interaction: Methods and Applications (Wiley Series in Computational Mechanics)

作者:Bazilevs, Yuri,Takizawa, Kenji,Tezduyar, Tayfun E.
発売日: 2013/02/11
メディア: ハードカバー

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↩

2019-02-21

VS codeに移行

雑記ライブラリ

Sublime text3をここ数年くらい使用してきたけど、日本語環境が使用できなくなったり、
Latex環境構築があまり楽でなかったりして、だんだん微妙に感じてきたのでVS codeに徐々に移行していきます。

インストールして、特に何か難しい操作もなく直感的にいろいろ遊べそう。ただLatex設定は少し手間取った。
下のまとめで全て解決したけど。
qiita.com

Latex環境としてOverleaf v2を使ってきたけど、特に人と共有して作業することもないので、
オフラインでこっちで作業することにしても問題なさそう。
microsoft社だけあってgitとの相性もいいし、早々にこっちに本格的に移行しよう…

2018-09-06

Qiitaに移行することにしました

雑記ライブラリ

ちょっと前までQiitaはニッチな情報系のエンジニア向けのサービスでしたが、最近は科学技術計算のエンジニアも結構Qiitaにいてて，まじめにFEMやマルチボディダイナミクスの定式化を紹介してたりします。

せっかくなので自分もQiitaに移ろうと思います。練習がてら、すげー短い記事を書いてみました。
qiita.com

コードの埋め込みも簡単そうだし、数式も書けそうだし、悪くなさそうです。

2017-05-14

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などは行っていたのだが、何やらゴミが残っていたらしい。

2017-05-02

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を記入

2017-05-01

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のバージョンは多少古いのかもしれない。