Fortranのオブジェクト指向を数値解析で使う

Advent Calendar 埋め草記事として、Fortranでのオブジェクト指向の良い所と感じている事を書いてみます。

#Fortranとオブジェクト指向について考えてみる
Fortranは、今(Fortranよどこへ行く？)も昔も数値解析を行うための言語であるという方向性はブレていないようです。何かしらの数値解析手法そのものを研究している方々にとっては、計算速度が最優先の場合もあるでしょうから、そのような場合にも計算にオーバーヘッドのあるオブジェクト指向を勧めるつもりはありません。

一方、数値解析を道具として用いて、工業製品の設計に活用したり、色々な試行錯誤をしなければならない人にとって、アイデアを気軽に実行に移せることは非常に重要だと思います。Fortranでオブジェクト指向を使うことの利点は、コードの再利用性が上がるという点に尽きると思ってます。

Fortranにオブジェクト指向の機能が無い(classが無い）としたら、新たに必要となるパラメータをサブルーチンへ渡す方法として次のようなものがあります。

• サブルーチンの引数を追加する。
• 型を新たに定義してそこにパラメータを含め、サブルーチン引数は追加しない。（ただし、引数の型名を変更したら、宣言を関連サブルーチン全てで書き直す必要がある。）
• モジュールにグローバル変数として定義する。

幸い(?)にもFotranにclassができたことで、このような引数追加や型宣言の書き直しが不要になります。

例えば、以下のような計算をしたいと思い立ったとき、オブジェクト指向を活用すれば非常に少ない手間で結果を得ることができます。

1. 多数のパラメータの組に対して、微分方程式の初期値問題を説いて比較したい。
2. 遺伝的アルゴリズムやPSOなどの最適化計算で、構造物の固有振動数を最大化する形状を求めたい。

このような計算は1パターンの計算に数分〜数時間（時には数日）かかる場合もあるため、できるだけ並列化も考慮したいところです。また、データと解法の組をセットにした型を定義しておけば、「数値解法は計算時間の早い非オブジェクト指向ルーチンで実装してメソッドとして構造型に登録」し、「データの取扱は記述が楽なオブジェクト指向で実装」し、さらに「OpenMPやCoArrayで並列化」するというように、Fortranの特長をいかんなく発揮できると思います。

#微分方程式初期値問題のパラメータスタディをやってみる
以前のQiita記事に書いたRunge Kutta法に関するモジュール(rungekutta)を使い、100通りのパラメータそれぞれの組(a,b,c)について、同一の初期値から微分方程式を解いた時の解を求めます。

##ソースコード
ソースコードではまずrungekuttaモジュールで定義したrk型を拡張し、パラメータa,b,cを型の成分として追加定義します。また、状態変数の初期値を型の成分として定義しています。さらに、出力ファイル名とユニット名も型の定義に含めました。

kparallel.f90

program rkparallel
  use,intrinsic :: iso_fortran_env
  use rungekutta
  !$ use omp_lib
  implicit none
  type,extends(rk) :: rkwithpara
    real(kind=real64) :: a=10_real64
    real(kind=real64) :: b=28_real64
    real(kind=real64) :: c=2.66666_real64
    real(kind=real64),dimension(3) :: x=[0.5_real64,0.5_real64,0.5_real64]
    character(len=:),allocatable :: fname
    integer :: unitnum
  end type
  type(rkwithpara),dimension(:),allocatable :: rkobj
  integer :: i
  character(len=4) :: txt
  real(kind=real64),dimension(3) :: tmp
  allocate(rkwithpara::rkobj(100))
  !$omp parallel
  !$omp do private(txt)
  do i=1,size(rkobj)
    !---
    ! open file and set parameter
    call rkobj(i)%setup(rhside,"rk4",1.0_real64)
    write(txt,'(i4.4)')i
    rkobj(i)%fname="rkout"//trim(txt)//".dat"
    write(output_unit,'(1X,"filename:",1X,A)')rkobj(i)%fname
    open(newunit=rkobj(i)%unitnum,file=rkobj(i)%fname)
    call random_number(tmp)
    rkobj(i)%a=rkobj(i)%a+tmp(1)*10_real64
    rkobj(i)%b=rkobj(i)%b+tmp(2)*10_real64
    rkobj(i)%c=rkobj(i)%c+tmp(3)*5_real64
    !---
    ! calc rungekutta and write data
    !$ write(output_unit,*)"Hello! N =", omp_get_num_threads(), " and I am ", omp_get_thread_num()
    associate(time=>rkobj(i)%time,x=>rkobj(i)%x)
      do
        if(time>50.0_real64)exit
        call rkobj(i)%solver(time+4e-3_real64,x)
        write(rkobj(i)%unitnum,*) x
      enddo
    end associate
    !---
    ! close file 
    close(rkobj(i)%unitnum)
  enddo
  !$omp end do
  !$omp end parallel
  !=============================================
  contains
    !=============================================
    function rhside(self,time,wk)
      class(rk),intent(inout) :: self
      real(kind=real64) :: time
      real(kind=real64),dimension(:) :: wk
      real(kind=real64),dimension(size(wk)) :: rhside
      select type(self)
      class is(rkwithpara)
        associate(a=>self%a,b=>self%b,c=>self%c)
          rhside(1)=a*(wk(2)-wk(1))
          rhside(2)=b*wk(1)-wk(2)-wk(1)*wk(3)
          rhside(3)=wk(1)*wk(2)-c*wk(3)
        end associate
      end select
    end function
end program

上記のメインプログラム(kparallel.f90)とRunge Kutta法で記載したモジュールを併せてコンパイルすると実行可能ファイルが作成され、さらにOpenMPを有効化すれば、時刻歴の計算をスレッド並列化することができます。
このように、Fortranのオブジェクト指向(型の拡張)を活用することによって、rungekuttaモジュールに一切手を加えることなく、計算の本質的な部分のみを考慮してコードを作成することが可能になっています。

##計算結果
100通り全てを表示するのは辛いので、ここでは2通りのパラメータの組から得られた結果を図に示します。同じ初期値であっても、パラメータが異なると描く軌跡が変わることが確認できました。

#おわりに
今回の例では計算が一瞬で終わりますが、もっと複雑な系では、オブジェクト指向を使う有難みが増してきます。数値解析でFotranのオブジェクト指向を使ってみてはいかがでしょうか。