sakamurray

32nd International Symposium on Shock Wavesで研究成果を発表しました．

投稿日: 2019年7月16日2019年8月3日 by sakamurray

2019年7月14日〜19日の間シンガポールで開催されている32nd International Symposium on Shock Wavesにおいて，本研究室の研究成果を発表しましたので，ご報告します．

演題：Numerical simulation of the shock-induced motion of a solid body using the overset grid functionality of OpenFOAM
著者：Y. Sakamura (登壇者)，K. Nakayama and H. Oshima

熱力学こぼれ話

投稿日: 2019年7月3日2019年7月10日 by sakamurray

本日の講義「エネルギー変換工学」の中で，以下の熱力学関係式を使いました．

$T = \left( \frac{\partial H}{\partial S} \right)_p \qquad \cdots (1)$

ここで， $T$ は温度， $H$ はエンタルピー， $S$ はエントロピー， $p$ は圧力です．授業後に回収した大福帳（コミュニケーションカード）に質問があったので，式(1)の導出過程を示しておきます．

まず， $H=H(S, p)$ とみなせば，その全微分は

$\mathrm{d}H = \left( \frac{\partial H}{\partial S}\right)_p \mathrm{d}S + \left( \frac{\partial H}{\partial p}\right)_S \mathrm{d}p \qquad \cdots(2)$

と表すことができます．一方，熱力学の第1法則から，

$\mathrm{d}H = T \mathrm{d}S + V \mathrm{d}p \qquad \cdots (3)$

が成り立ちます．ここで， $V$ は体積です．

式(2)と(3)の右辺を見比べてみれば，式(1)ならびに

$V = \left( \frac{\partial H}{\partial p} \right)_S \qquad \cdots (4)$

が成り立つことがわかります．

ところで，式(1)をみると，示量性状態量 $H$ を示量性状態量 $S$ で微分したものは示強性状態量 $T$ になっていることがわかりますね．一方，式(4)からは，示量性状態 $H$ を示強性状態量 $p$ で微分したものが示量性状態量 $V$ になることがわかります．だから何だと言われればその通りなのですが，何かのときに役立つかもしれませんww．

上述した示強性状態量と示量性状態量の関係性は，原島鮮先生の教科書「熱力学・統計力学（改訂版）」培風館（1978），p. 78に紹介されています．この本には他にも多くのユニークな「こぼれ話」が紹介されていて，まるで講義を聴いているように読み進めることができる，お薦めの教科書です．

Juliaに傷心（ハートブレイク）

投稿日: 2019年6月4日2019年6月5日 by sakamurray

2019年5月18日（土）に富山CiCビルで開催された第72回オープンCAE勉強会＠富山で，プログラミング言語Juliaに関する話題提供を行いました．その内容をざっくり纏めてみましたので，ご紹介します．

報告内容

Juliaとは？
Juliaのインストール
Jupyterでの利用
Pythonとの速度比較（暫定版）
まとめ

Juliaとは？

Juliaは2012年に発表された新しいプログラミング言語であり，個人的には次のような特徴が気に入っています．プログラミングの導入教育には最適だと考えられます．

Python同様動的型付け言語だが，静的型付け言語に迫る実行速度を持つらしい．
計算実行時にコンパイルするJIT(Just-In-Time)コンパイラを採用しており，PythonやRuby等のインタプリタと同じように簡単に利用できる．
Python同様，Jupyterをインターフェイスとして利用できる．

Juliaのインストール

公式サイトで各種OS用のインストーラが配布されています．私はmacOS版を手元のMacBookにインストールして使っています．

Mac版Juliaは，通常のアプリケーションとして登録されますので，ダブルクリックするとターミナル上でREPL (Read-Eval-Print Loop)環境が動作します．パッケージの追加・アップデートはこの環境で行うことができます．さらに，”?”をつけてコマンド等を入力すると詳しいヘルプ画面が表示されます．

Jupyterでの利用

私にとって最大のJuliaの魅力は，Jupyter notebookやJupyter lab上で実行できることです（そもそもJupyterの”Ju”はJuliaを指しているらしいので，当たり前と言えば当たり前ですが）．

既にJupyterがインストールされていれば，Jupyterとの連携は簡単です．私の環境(Mac+Anaconda)の場合，JuliaのREPL環境で”]”をタイプしてpkgモードにした後，次のようにしてIJuliaパッケージをインストールすれば，Jupyterのカーネル・リストにJuliaが追加され，利用できるようになりました．

(v1.1) pkg> add IJulia

Pythonとの速度比較（暫定版）

Pythonプログラムと，それをJuliaに移植したプログラムとで計算速度の比較を行ってみました．

計算対象は，Barba先生のCFD-Python Step 5で取り上げられている2次元線形移流方程式

$\displaystyle {\partial u \over \partial t}+c{\partial u \over \partial x} + c{\partial u \over \partial y}=0$

の差分近似

$\displaystyle u^{n+1}_{i, j} = u^{n}_{i, j} - c \Delta t \left( {u^n_{i, j}-u^n_{i-1, j} \over \Delta x} + {u^n_{i, j}-u^n_{i, j-1} \over \Delta y} \right)$

です．これを以下の初期条件，境界条件の下で解いてみました．

初期条件：

$u(x, y, 0) = \left\{ \begin{array}{c} 2 \quad \mbox{for} \quad 0.5 \leq x, y \leq 1 \\ 1 \quad \mbox{for everywhere else} \end{array} \right.$

境界条件：

$u(x, y, t) = 1 \quad \mbox{for} \left\{ \begin{array}{c} x=0, 2 \\ y=0, 2 \end{array} \right.$

空間の分割数はx, y方向ともに80分割( $n_x=n_y=81$ )で，100タイムステップ $n_t=100$ の計算にかかる時間を計測しました．

CFD-Pythonのコードのコア部分は以下のようなものです．

for n in range(nt + 1):
un = u.copy()
for j in range(1, ny):
for i in range(1, nx):
u[j, i] = un[j, i] – c * dt*( (un[j, i] – un[j, i – 1])/dx+ (un[j, i] – un[j – 1, i])/dy)

これをJuliaに移植すると

for n in 1:nt
     un = copy(u)
     for i in 2:nx-1
          for j in 2:ny-1
               u[i, j] = un[i, j] – c * dt*((un[i, j]-un[i-1, j])/dx+(un[i, j]-un[i, j-1])/dy)
               u[1, j] = 1
               u[nx, j] = 1
          end
          u[i, 1] = 1
           u[i, ny] = 1
     end
end

シンプルで，初心者にも分かり易いコードだと思います．このコードで計算した結果，Pythonコードでは2.5 s程度掛かっていた時間がJuliaコードでは0.41 sに短縮され，6倍程度高速になりました．

元になったPythonコードの計算過程も分かり易いのですが，3重のfor-loopで書かれていて，ものすごく遅いため，普段Pythonを使っている人は配列操作(Array Operation)を用いるようです．実際，CFD-Pythonでも配列操作を利用するコードが紹介されていて，こちらで計算すると6 ms程度で計算が終了しました．私が書いたJuliaコード（3重のfor-loopを用いたもの）の60倍以上のスピードが出ています（速い！）．

せっかくJuliaを使ったのに配列操作を使ったPythonに惨敗という結果になりまして，傷心気味だったのですが，@ceptreeさんの記事「Juliaの速さを体感する」を見つけ，メインループの関数化を試みたところ，Juliaでの計算時間を0.14 sにまで短縮することができました．それでもなお，配列操作を使ったPythonコードにスピードでは圧倒的に負けていますので，引き続き，Juliaの高速化について調べてみたいと思っています．