雛形書庫

An Unmoving Arch-Archive

Computational Fluid Dynamics by Julia

JuliaでCFDに入門します。

Original Paper

Disclaimer

内容の正確性には可能な限り留意しておりますが、必ずしも正確性を保証するものではありません。 内容に基づく一切の結果について、一切の責任を負いかねます。 予めご了承ください。

はじめに:この論文について

ポッドキャストEngimono第2回ではJuliaを学び、また第3回ではサイクリスト第4回ではF1の空力シミュレーションに関する論文を読んできた。 お次は自分で実践してみましょうということで、Juliaで流体シミュレーションをやってみる論文。 タイトルにあるCFDとは差金決済取引のことではなく、Computational Fluid Dynamics (CFD) の略。

アブストラクトを読んでいくと、冒頭にあるとおり、学部学生や修士学生がCFDの基礎を学ぶためのJuliaモジュール開発が本論文のメインテーマとのこと。 そのためか本文の構成からは、research articleというよりは講義用の教材に近い。 本論文がカバーしているのはCFDにおける時空間の離散化手法、陽解法と陰解法、高次精度衝撃波捕獲法、など。対象としている方程式は、線形移流方程式、非粘性バーガース方程式、2次元ポアソン方程式、あとは1次元オイラー方程式、2次元非圧縮性ナビエ–ストークス方程式。 取り扱いは、保存形と非保存形での有限差分法。 というふうに眺めてみると、結構手広く紹介されていることがわかる。

そしてページ数が77ページと、とても長い……ので、論文全体を網羅するのは諦めて、また個々の問題に着目することもせず、イントロと本文に基づいて少し俯瞰的な読みを試みたい。

Introduction: why Julia?

CFDの勉強は偏微分方程式の離散化手法の紹介に始まるものが多いけれども、学生にとってみれば座学だけではなくて、ハンズオンで実際にプログラムを動かしてあれこれやってみるほうが理解が捗るだろう。 この論文が紹介しているCFD Juliaというのは、前節で挙げたような種々の問題を取り扱うためのプログラミングモジュールである。 世の中にはこうしたプログラミングモジュールはすでにいくつかあって、例えばナビエ–ストークス方程式であれば CFD Pythonがある。CFD Pythonは名前の通りPythonで書かれていて、他にもPythonをベースにした学習用教材がいくつかある (Pyclaw, HyperPython)。

そうした背景のもとで、このプログラミングモジュールでは(Pythonではなく)Juliaを使った。 科学技術計算の世界ではtwo-language problem1というのがあって、科学技術を計算したい人々にとっては、まずはお手軽に書けるPythonMatlabやRを使ってプログラムを作りたい。 けれども問題の規模が大きくなってくると、やっぱり性能面での補強が必要だということで、FortranC++でプログラムを書き直す。 この書き直しの手間が結構馬鹿にならなくて、シミュレーションコードを開発するうえで常に障壁2となっている。 Juliaはこの問題を解決しようと開発された言語であって3、その恩恵にあやかろうとしたもの。

Advantages of Julia in Numerical Simulation

Juliaはその高速性能で知られているけれど、ではシミュレーションが高速化することで具体的にどのようなメリットがあるだろうか。 シミュレーションの世界ではターンアラウンドタイムというものがあって、これは解析しようと着手してから結果が得られる4までにかかる時間を指す用語。たとえば:

ミシガン州トロイ、2010年4月26日 ― ミシガン州に本拠を置くAltair Engineering, Inc.は、OEMネイティブのCADデータからフル車両衝突有限要素モデルのメッシュ作成-アセンブル-シミュレーションまでの所要時間をわずか24時間にまで短縮することに成功し、製図台からCADステーションへの移行にも匹敵する、先進的製造システムのための飛躍的な進歩を実現しました。この画期的な成果により、同種の解析で通常必要とされるターンアラウンドタイムが2週間から4週間ほど短縮されることになります。

――Altairがフル車両衝突シミュレーションの限界を打破: CADから結果出力までを24時間以内で実現

シミュレーションに要する時間が短いほどこのターンアラウンドタイムも短くなるから、設計する側としては同じ時間でより多くの試行をこなせるようになり、それは設計の競争力向上へと繋がる。

あるいは、ソフトウェア保守の観点からもJuliaを使うご利益がありそう。 従来使われてきたFortranC++ドキュメンテーションツールといえばDoxygenが有名だけど、JuliaにもDocumenter.jlなるツールが存在する。 それにJulia Continuous Integration (CI)としてCIツールも整備されているので、これらはそのままシミュレーションコード開発にも役立つはずだ。

加えて、グラフや可視化図の描画環境との親和性という観点からも、Juliaには有利な側面があるかもしれない。 従来の性能重視の言語では、C++はまだしもFortranで何か描画しようというのは荷が重くて、たいていはシミュレーション結果として何かしらのファイルを出力するようになっている。 こうして出力したファイルを改めて別の描画ソフトで読み込み、グラフや可視化図を作るというのが一般的な操作手順だと思う。 けれどJuliaであれば、例えばPlotsにあるように言語自体に便利な可視化機能が備わっているから、シミュレーション結果を確認するうえでたやすい。 もちろんシミュレーションの規模が大きくなるにつれて可視化の作業も大変になるけれど、この論文で取り上げられている程度の大きさであれば、Juliaが持っている可視化パッケージで十分対応できるだろう。

Partial Differential Equations

CFDでは偏微分方程式 (Partial Differential Equations) を取り扱う。

偏微分方程式にも3つのタイプがあって、Section 2では放物型偏微分方程式の問題を取り扱う。 Section 3では双曲型偏微分方程式の問題を取り扱う。 Section 6では楕円型偏微分方程式の問題を取り扱う。 放物型と双曲型と楕円型で何が違うかというと、情報の伝わり方が異なる:

f:id:hina747:20210131170216p:plain
偏微分方程式の性質(棚橋 隆彦著『はじめてのCFD - 移流拡散方程式 -』(コロナ社)より引用)

流体シミュレーションで用いるナビエ・ストークス方程式は上記3つの型が混合したケースであって、ゆえにそれぞれの型の性質、および対応する諸問題を学ぶことには意義がある。

この論文で紹介されているCFD JuliaのソースコードGitHubで公開されているので、必要に応じて参考にされたい(ただし、この記事では論文にあるコードスニペットを参照するにとどめ、実際のコードまでは立ち入らない)。 Table 1にある通り、全部で22のセクションから構成されていて、各セクションの内容は表にあるDescriptionの通り。

Implementation

実際の実装については、論文で詳しく解説されているので割愛。 代わりにいくつかの一般的なトピックを抑えておきたい。

Error Evaluation

この論文で取り扱われている問題のうち、厳密解 (exact solution) の存在するものについては誤差評価をすることができて、その場合は式 (11) のように厳密解と解析解 (numerical solution) の差(の絶対値)を取ってやる。この誤差はFigure 2の右側のような感じで、離散化誤差 (discretization error) としてプロットすることができる。 両対数グラフにおいて、この誤差がどれくらいの傾きで減少していくかを離散化の精度と呼んでいて、たとえばFigure 18の右側ではプロットがSlope=2の直線に沿って減っていくので、これは2次精度ですね、となる。 ここで、誤差の値そのものの大きさと精度は別物なので、高次精度の方法だからといって誤差の値そのものも小さいとは限らない。 高次精度とは誤差の減っていく様子を言っているので、たとえばp. 7にある記述

The discretization error is less for Runge-Kutta numerical scheme compared to the FTCS numerical scheme due to a higher-order approximation of temporal term.

は、読みかたによっては間違いとまではいえないにしても、ちょっと誤解を招きそう。

CPU Time in Simulation

論文のTable 1ではアスタリスクのついた項目があって、いわくPythonによるプログラムも同梱されているとのこと。 おそらくこれに対応する形で、論文にはシミュレーションに要するCPU時間の比較表がいくつか登場する。 そのうちのひとつ、Table 4を見るとJuliaはPython (pyFFTW) よりも計算時間が長い……ように見えるけれど、その一因としてpyFFTWパッケージはとても高速なCのライブラリを使っているから、とのこと。 また結論にあるように、論文で示されているコードは教育用途の5ものであって、計算速度の面で最適化されていないという側面もある。 イントロでは一般的な話、Juliaの高速性能が謳われていながらも、それを論文中の問題で実感することはやや難しいかもしれない。

Summary

JuliaでCFDに入門してみて、論文での最終課題である2D incompressible Navier-Stokes equations、2次元空間のシミュレーションまでできるようになった(はず)。これを冒頭で挙げたようなサイクリストやF1の空力シミュレーションに適用するとなると、3次元空間シミュレーションへの拡張、シミュレーションに使う格子の準備など、まだやるべきことはたくさんある。 それでも基本的なところを理解したり、あるいは検証過程で基礎的な問題に立ち戻りたくなったときに、こうした小さなプログラムは役に立ってくれるのだろう。 「(Juliaを)HPC環境に完全に馴染ませるのは難しくとも、簡単なsnippetくらいであれば移行したい」というのはかつてこのブログで書いた文言だけど、今回の論文にあったコードの小片たちを、そうした使い方をする際の考え方の一助にできると良さそう。■


  1. Misreading Chatエピソード27での森田さんの解説が詳しい

  2. 原文ではhindrance(辞書を引いた)

  3. Juliaの生産性と高速性能を両立する設計については、本Podcastのエピソード2で論文を読んでいる

  4. もっとも広義の解釈であり、文脈によってはもっと狭い意味になることもある

  5. 原文ではpedagogical(辞書を引いた)