Juliaでサメを泳がせます。
Reference
- Gabriel Weymouth, "WaterLily.jl: Real-time fluid simulation in pure Julia," JuliaCon 2021.
Disclaimer
内容の正確性には可能な限り留意しておりますが、必ずしも正確性を保証するものではありません。 内容に基づく一切の結果について、一切の責任を負いかねます。 予めご了承ください。
はじめに
#interaxion ep.24で取り上げられてたJuliaConにやっと追いついた。ゲストのとみやさんのお話も聴いた。講演アーカイブ残ってるのほんとありがたいhttps://t.co/Wb1riOVJnh
— 雛形 (@hinahypersonica) 2021年9月18日
7月末に開催されていたらしいJuliaCon 2021を覗いているとシミュレーションのお話を見つけた。 タイトルにある通り、流体シミュレーションをJuliaだけでやりました、というもの。
なお同じ著者によるThe Julia Language Blogへの投稿記事では、おそらくこの講演で紹介されているコードを使った応用例であろう、サメの泳ぎのシミュレーションを見ることができる。
本記事はこのサメを泳がせる前の下準備として、講演動画の内容を追いかけていく。 そして今後追加で書く（であろう）予定の後編では、応用事例としてのサメの泳ぎについて見ていくこととしたい。
なぜやるか
Youtubeにあるアーカイブ講演動画は約8分なのでさくっと見られる。 学会のオンライン開催が普及したことで、こうして誰でも講演コンテンツにアクセスできるようになったのは良い傾向だと思うし、リアルタイムで追っていなかった身からすると素直にありがたいなと感じる。
さて筆者調べによれば、この動画では以下の3つのポイントが解説されている。
- コードの紹介
- なぜやるか (Why Julia?)
- コードの構成
まず冒頭でコードの概要とできることの紹介（シミュレーションの動画付き）、そしてエディタの画面を共有したシンプルなデモンストレーションが行われる。 2番目になぜJuliaで取り組んでいるかの説明、そして3番目にパッケージが持っている複合型 (composite types) について解説がなされている。
ここで2番目のなぜやるかの部分に着目すると、いわく機械学習、自動微分、GPU計算といった昨今のトレンドをシミュレーションに取り込んでいくうえで、従来のFortranよりもJuliaのほうがやりやすかろうということ。 シミュレーションをFortranでやる一方、機械学習などの部分はPythonによる既存の成果物を使おうとすれば、必然的に2つの言語を駆使することになり、そこには相応の手間がかかるであろうことは想像に難くない。 今回このテーマでJuliaを深堀りしている個人的な動機もまさにここにあって、シミュレーションにいながらにして別の分野との連携が容易にできそうなことが挙げられる。 そしてFortranとは違ったコーディングの流儀を見てみたい、というのがもうひとつの理由にある。
ちなみにコード名称にあるWaterLilyとは、Water（水）+Lily（百合）ではなく"water lily"の一単語でスイレンの意である（と調べてはじめて知った）。百合要素は含まれないことに注意したい
流体の圧縮性と解法
GitHubのWaterLily.jlにある説明を見ると、このコードが取り扱うのは非定常非圧縮性ナビエ・ストークス方程式 (the unsteady incompressible 2D or 3D Navier-Stokes equations) とある。 非圧縮性の近似においては音速が無限大になると考えても良いとあり、流体中の圧力変動は音速で伝搬するから、圧力変動が無限大の速度で伝搬するということになる。 すなわち以前の論文読みの参考文献として取り上げた「点Pのじょう乱はすべての領域に伝わる」ということに相当し、圧力の偏微分方程式が持つ性質としては楕円形となる。
楕円形の偏微分方程式は反復法で解かれる。 非圧縮性流体を解くコストの多くは圧力の方程式を解くための反復法に費やされるとされ（要出典）、その例として以前の論文読みで取り上げた姫野ベンチマークとは "ポアッソン方程式解法をヤコビの反復法で解く場合に主要なループの処理速度を計るものです。" つまり反復法の部分をホットスポットとみなして、そこに要する時間を測定するベンチマークであった。 したがって非圧縮性流体のシミュレーションではここを高速化できるかどうかが効率的な計算のポイントになってくる。
WaterLily.jlにおいては講演動画 (5:18~) にあるとおり、高速化のためにマルチグリッド法を用いている。 マルチグリッド法は以前の論文読みでも出てきた方法で、CFD Juliaの項目17、論文では6.3節に登場している。 この事例に基づけば、残差を落とし切るのに通常の反復解法ではO(105)の反復回数が必要なところを、マルチグリッド法ではO(1)まで落とすことができるということで、反復回数の削減に絶大な威力を持っている。 この後で実際にWaterLily.jlを動かしてみるとその高速性が実感できるのだけど、コードの速さの秘密はこのあたりにありそうだ。
サンプルコードの実行
講演でのサンプルコードとして出てきたスニペットを以下に引用する。
using WaterLily
n,m = 3*2^6,2^7
center,radius = m/2,m/6
body = AutoBody((x,t)->√sum(abs2,x .- center) - radius)
sim = Simulation((n+2,m+2), [1.,0.], radius; body)
sim_step!(sim,10)
include("TwoD_plots.jl")
sim_gif!(sim;duration=20,step=0.25)
途中のinclude("TwoD_plots.jl")については、WaterLily.jlをチェックアウトして得られるexamples/TwoD_plots.jlを実行ディレクトリに置くことでエラーなく実行できた。
実行すると適当なディレクトリに講演動画にあったようなアニメーションgifが生成される。
実行はできたが、どこかで流体の素性を変えられるはず…と思って調べてみたところ、複合型Simulationがキーワード引数として動粘性係数νを取れることがわかった。
VSCodeで複合型の名前にマウスカーソルを乗せると、引数の説明が出てくる。めちゃくちゃ便利
この値はデフォルトではv=0が入っているらしいので、適切な値を代入してやることで流れの様相を変えることができる。
ν=1.0としたとき。渦は見えない。
ν=0.2としたとき。周期的な渦が見られる。
円柱を過ぎる流れの論文によれば、例えば第11図と見比べるとこれらは「双子渦領域」「純カルマン渦領域」に対応しているように見えて、それっぽい流れが得られていることがわかる。 でもそれっぽい以上の信憑性があるかというと、今のところはまだ確たることはいえない――いつも思う（思い悩む）のはこうしたそれっぽさと、実際に正しい・あるいは実際に使えることとのギャップであって、そうしたものをきちんと掘り下げたうえで使いこなせるようになりたいなあと感じている（なお、この講演で用いられているシミュレーションコードの手法的な部分は既報ですでにValidationがなされている様子なので、大方問題のない流れが得られているのだろうと考えている）。
おわりに
Juliaでサメを泳がせるつもりが前口上が長くなってしまい、講演動画にある小さなチュートリアルを動かしたところで一段落した。 この記事はたぶん後編も執筆される予定なので、そこでは冒頭のJuliaブログ記事にあるようなシミュレーションを動かしてみることをやってみたい。
講演動画でも触れられていたように、かつてFortranで書かれていたようなコードを今ではJuliaで書いてみた、というのが一つの新規性、講演をなすことのjustificationとされていた。 ここでふと思い出したのがRebuild.fmのEp.169、森田さんがオブジェクト指向について語る場面で、今となってはオブジェクト指向で実装することは何ら真新しさを持たない、ということが半分冗談めかして語られている。いわく、当時はそうでもしないと研究予算も取れなかったよと。 しかし今回のJulia、あるいはHPCについても同様に、やがては計算効率と開発効率を両立することは珍しいものではなくなるのかもしれないなとも思う。 かつてのオブジェクト指向と同様に、今から5年後か10年後かわからないけれど、その昔two language problemというものがあってね……と古き日々を回顧できる未来がくることをささやかに期待したい。■
ライセンス
WaterLily.jl / MIT "Expat" License
https://github.com/weymouth/WaterLily.jl
Copyright (c) 2020: Gabriel Weymouth.
Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:
The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
