Fortranを使って流体シミュレーションで大石泉すき

概要

「大石泉すき」Advent Calendar 2019 19日目の記事です．
Fortranを使って「大石泉すき」をやります．

Fortran

Fortranとは，最古の高級言語であるFORTRANを祖先にもつ，数値計算向けのプログラミング言語です．どのくらい古いかというと，OSが誕生するよりも古くから存在していました．また，アメリカの標準規格団体誕生以前から存在しているので，Fortranをさわっていると，その名残を垣間見ることができます．

大石泉はプログラミングが得意です．現在のプログラミングパラダイムとFortranの設計思想には乖離があるものの，大石泉はいい子なので，Fortranを絶滅危惧種とディスったりはしないでしょう．「なんでこんな書き方なんだろう？」と疑問に思っているところに，八神マキノがレクチャーする場面がありそうな気がします．

Fortranを介して展開されるマキいずとかよくないですか？？？？？

Fortranのコンパイラ

フリーのFortranコンパイラはいくつかあります．

• gfortranは，gccなどと同じGNUのコンパイラです．OSによってはパッケージ管理システムからコマンド一発でインストールできます．Windowsでは，MinGWやMSYS2など，いくつかのバイナリ版が用意されています．gnuのWikiにまとめられています．
• GPUで動くプログラムを作りたい場合は，PGI Fortranが利用できます．PGI FortranにはCommunity editionが存在しており，少々制約はありますが，無料で利用できます．コンパイラは英語ですが，日本の代理店が日本語の記事を公開してくれています．
• Fortranの速さを体感したいのであれば，Intel Fortranがよいでしょう．登録は必要ですが，学生およびオープンソースソフトウェア開発者が無料で利用できるライセンスがあります．詳しくはIntelの公式サイトを見てください．

同人誌ですが，拙著 モダンFortran入門(環境構築編)にて上記3コンパイラの環境構築について解説しています．

Fortranで大石泉すき

まずは，標準出力に「大石泉すき」を出力してみます．
標準出力をする場合は，print文を用います．*は書式を指定するところですが，＊(Asterisk)にしておくと，書式をほどよく判断して決めてくれます．

izumi.f90

program main
    implicit none

    print *, "大石泉すき"
end program main

> gfortran izumi.f90
> a
 大石泉すき

gfortranでは，出力ファイル名を指定せずにコンパイルすると，a.exe（Unix系ではa.out）という名前の実行ファイルを出力します．aというプログラムを実行すると「大石泉すき」と出力されるあたり，語彙力がなくなっていっている感じがします．

次に変数を使ってみます．文字列に「大石泉すき」を代入して，標準出力に文字列の長さと文字列を出力しています．

program main
    implicit none

    character(:),allocatable :: my_heart

    print *,len(my_heart)

    my_heart = "大石泉すき"
    print *,len(my_heart), my_heart
end program main

     8429456
          10 大石泉すき

大石泉に出会うまでは，my_heartは正常ではありませんが，大石泉に出会うと，my_heartは「大石泉すき」で満たされていることがわかります．大石泉すき

Fortranでidol++

Fortranにはインクリメント演算子++が存在しません．そのため，残念ながらidol++はできません．

また，Fortranは関数の戻り値を捨てることや，単項演算子でオペランドの値を変更することができないので，どう頑張ってもidol++はできません．

大石泉はきっと諦めないので，それっぽいことをやってみましょう．

module mod_Idol
    use,intrinsic :: iso_fortran_env
    implicit none
    private
    public :: plusplus
    public :: operator(.pp.)

    interface operator(.pp.)
        module procedure incrementInt
    end interface

    contains

    subroutine plusplus(val)
        implicit none
        integer(int32),intent(inout) :: val
        val = val+1
    end subroutine plusplus

    function incrementInt(val) result(inc)
        implicit none
        integer(int32),intent(in) :: val
        integer(int32) :: inc
        inc = val + 1
    end function incrementInt
end module mod_Idol

program main
    use :: mod_Idol
    use,intrinsic :: iso_fortran_env
    implicit none

    integer(int32) :: idol

    rankup : block
        integer(int32) :: rank
        idol = 0
        do rank = 1, 10
            call plusplus(idol) ! サブルーチンを使う場合
            ! idol = .pp.idol   ! ユーザ定義演算子を使う場合
            print *,idol
        end do
    end block rankup
end program main

moduleは，Fortranにおけるパッケージ管理のような仕組みです．モジュールmod_Idol内で，サブルーチンplusplusと関数incrementIntを定義していますが，incrementIntは外部に公開されていません．サブルーチンは，戻り値のない関数です．サブルーチンと関数をあわせて，procedure（日本語では手続）とよびます．

Fortranでは，ユーザが独自の演算子を作る事ができます．ここでは，引数に1足した値を返す関数を.pp.という演算子で参照できるようにしています．実体は関数incrementIntですが，関数を公開する代わりに演算子.pp.を公開しています．

実行するとランクが上がっていくことを確認できます．

           1
           2
           3
           4
           5
           6
           7
           8
           9
          10

流体シミュレーションで大石泉すき

Fortranはかつて一世を風靡したプログラミング言語ですが，現在は引く手あまたというわけではなく，主な用途である数値計算の分野で力を発揮しています．

京の後継となるスーパーコンピュータ 富岳でも，重点的に取り組むべき社会的・科学的課題（重点課題）のほとんどでFortranが使われる予定です．

流体運動の支配方程式

非圧縮性粘性流れの支配方程式は，次式で表されます．

\begin{eqnarray}
&&\nabla\cdot\boldsymbol{u} = 0 &&(1)\\
&&\frac{\partial \boldsymbol{u}}{\partial t} + (\boldsymbol{u}\cdot \nabla)\boldsymbol{u} = -\nabla p+\frac{1}{Re}\nabla^2\boldsymbol{u}&&(2)
\end{eqnarray}

$\boldsymbol{u}$は速度，$p$は圧力です．物理的には，式(1)は質量の保存，式(2)は運動量（と運動エネルギ）の保存を意味しています．

わかりやすく？説明すると，$\boldsymbol{u}$は各人の心で，内なる大石泉すき活力と見なせるでしょう．$p$は周囲の影響です．式(2)の左辺第2項は，大石泉を好きなことによってさらに好きになっていく効果を表しています．式(2)右辺第1項は，周囲の声（止めろと妨害するか，もっとやれと応援するか）を表しています．式(2)右辺第2項は，大石泉が好きなことによって疲弊する効果を表しています．

式(2)右辺第2項にかかっている係数$\frac{1}{Re}$の分母$Re$は，レイノルズ数とよばれるパラメータですが，ここでは「大石泉すき数」と称することにします．$\frac{1}{Re}$は，大石泉を好きなことによって疲弊する効果に係数としてかかっているので，大石泉すき数（大石泉をすきな度合い）が大きくなると，疲弊しなくなると解釈できます．

レイノルズ数$Re\rightarrow\infty$の極限では，流体はその粘り気の影響を無視することができる理想流体となります．大石泉すきの極限では，疲れも周囲の声も気にならなくなります．流体にとって理想的な状態は，大石泉Pにとっても理想的というわけです．

物体表現

流れの中に大石泉すきを置いて，その周りの流れをシミュレーションします．

物体を表現するために，Navier-Stokes方程式の右辺に，もう1項追加します．

\frac{\partial \boldsymbol{u}}{\partial t} + (\boldsymbol{u}\cdot \nabla)\boldsymbol{u} = -\nabla p+\frac{1}{Re}\nabla^2\boldsymbol{u} +\lambda\chi(\boldsymbol{u_s}-\boldsymbol{u})

右辺第3項は，penalization項とよばれます．$\lambda$はペナルティパラメータとよばれる値，$\chi$は，物体中で1，流体中で0となるマスク関数です．$\boldsymbol{u_s}$は物体の移動速度で，物体が静止していれば$\boldsymbol{u_s}=\boldsymbol{0}$です．

penalization項を追加した式は，浸透流の支配方程式（Brinkman-Navier-Stokes方程式）との類似性があります．その類似性に着目すると，ペナルティパラメータは物体の浸透係数の逆数と解釈できます．つまり，スポンジのような多孔質の物体の目詰まり具合を表すのが$\lambda$，形状を表すのが$\chi$というわけです．

マスク関数$\chi$を画像から生成します．その画像は，Pythonを使って文字列から生成します．

グレイスケールの画像は，黒が0，白が255なので，0~255を0~1の範囲に変換し，白黒を反転することでマスク関数として用いる事ができるようになります．

流体計算では，値が不連続的に変化すると計算が上手くいかなくなる場合があるので，PILのぼかしフィルタPIL.ImageFilter.GaussianBlur(0)で文字の縁を若干ぼかしています．

縁をぼかし，0~1の範囲に変換された画像の各ピクセル値を出力し，Fortranのプログラムから読み込んでマスク関数値としました．

計算方法

時間発展にはFractional Step法を用います．

まず，前進Euler法で時間微分項を離散化し，計算に用いる各物理量の時刻を定めます．

\begin{eqnarray}
&&\nabla\cdot\boldsymbol{u}^{n+1} = 0\\
&&\frac{\boldsymbol{u}^{n+1}-\boldsymbol{u}^{n}}{\varDelta t} + (\boldsymbol{u}^{n}\cdot \nabla)\boldsymbol{u}^{n} = -\nabla p^{n+1}+\frac{1}{Re}\nabla^2\boldsymbol{u}^{n}
\end{eqnarray}

次に，圧力勾配を無視した式を使い，仮速度$u^*$を求めます．

\boldsymbol{u}^{*} = \boldsymbol{u}^{n} - \varDelta t (\boldsymbol{u}^{n}\cdot \nabla)\boldsymbol{u}^{n} + \varDelta t  \frac{1}{Re}\nabla^2\boldsymbol{u}^{n}

仮速度の$\nabla\cdot\boldsymbol{u}^*$から，圧力Poisson方程式を解いて圧力を得ます．

\nabla^2p^{n+1} = \frac{1}{\varDelta t}\nabla\cdot\boldsymbol{u}^*

最後に，圧力を用いて仮速度を修正し，次の時刻の時間を得ます．

\boldsymbol{u}^{n+1} = \boldsymbol{u}^{*}-\varDelta t\nabla p^{n+1}

この計算手順は，最初，周囲の声を無視して活動を行い，その後周囲の影響を反映して最終的な状態を決定していると捉えることができます．

これらの式を有限差分法を用いて離散化しました．

計算結果

計算は，大石泉すき数$Re=1000$として行いました．この条件での計算の難易度は，有償スタージュエル1000個を入手するのと同程度だと思います．

つまり，ジュエル購入の流体計算の経験が少ないと少し難しいですが，息をするようにジュエルを購入するようになるある程度慣れてくると，大したことはないといったところです．

計算結果をアニメーションgifにしてみました．矢印は速度$\boldsymbol{u}$，背景の色は圧力$p$（赤が高圧，青が定圧）です．また，マスク関数$\chi=0.5$で等高線を引き，文字を表示しています．

最初は「ふーん、あんたが大石泉？まぁ、すきかな」くらいですが，時間の経過と共に流れが複雑化していき，大石泉すきを通過するころには心に大きなうねりが生じています．大石泉すき

ある文字から剥離した渦が次の文字にぶつかり，それらの文字の間で定在的な渦が発達します．それらの渦が何らかの拍子で下流に放出されると，同様に発達した渦と合体し，大きく成長していく様子が確認できます．

まとめ

Fortranを使って「大石泉すき」をやりました．

Fortranは，現在多数派のプログラミング言語ほどのシェアはもうありませんが，根強い支持者がいます．こういう書き方をしてみると，アイドルっぽいですね．

Fortranを習得するのは簡単ですし，流体のシミュレーションも，今回程度の規模であれば，さほど難しくはありません．慣れれば，0から作っても1時間ほどしかかかりません（計算には一晩かかりますが）．もし気になった方がいれば，ぜひさわってみてください．川島瑞樹さんも「始めるのに遅すぎることはない」とおっしゃっています．