fortranで構造体とオブジェクト指向とシミュレーション

初めに

初投稿です．拙い部分しかないと思いますが，まだまだ初学者ですゆえどうかご了承ください．

この記事はfortranでスパコンなどの大規模並列計算機を使うシミュレーションプログラムを作成するためのヒントとなればいいなあという記事です．ですが，まず初めにこれからシミュレーションプログラムを作成しようと思われる方はこの記事を読む前にこれとこれを買って読んでください．読むべきです．話はこれからです．シミュレーションにおける，まさにバイブルとも呼べる本です（個人の感想です）．私は本当にこの本と同じ時代に生まれたことに感謝しています．シミュレーションでなにをすべきか，全てがここに書いてあります．

ではなぜいまさら私ごときがこんな記事を書くかというと，確かに本の中には「これはC++でもfortranでも同じである」とは書いてあるのですが，中身はC++で書かれてるんですよねー．しかしオブジェクト指向に関するfortranの記事は少ないという…．

というわけで，まずは内容で一番最初に躓いた，構造体とポインタについて書いてみようと思いました．いまのところ以下のやり方で特に不具合はありませんが，もしかしたら問題があるかもしれません．詳しい方は是非コメント欄に指摘ください．全力で勉強させて頂きます．

データモデルの表記法

本の中では有限要素法を例に，メモリレイアウトについて言及されています．以下の図のようなデータ構造を保持したいとします．

米印はポインタであることを表しています．私は有限要素法に明るくないのですが，それでも全てのノードでループを回したいとき，エレメントごとにループを回したいときがそれぞれあると思うのでこのメモリ配置がベストなんだと思います．

fortranのコード

あとは上図を実現するコードを書きましょう．

Node.f90

! クラス名にアンダーバーをつけてモジュール名を定義． 
module Node_
! モジュール名＝クラス名の場合，インデントを下げない方が見やすそう．（エディタが許すかは分かりませんが・・・） 
implicit none
public :: Node

! ここからいわゆるヘッダーに相当する部分
type Node
    private
    real(8) :: x(2), u(2)

contains
    procedure :: init => Node_init
end type

! ここからいわゆるソースに相当する部分
contains

subroutine Node_init(this, x, u)
    class(Node) :: this
    real(8), intent(in) :: x(2), u(2)

    this%x = x
    this%u = u
end subroutine

end module

他言語で言われているように，１モジュール１ソースコードが良さそうです．また，fortranは言語仕様上１モジュールに複数クラスを定義できますが，管理が破綻しそうなので１モジュール１クラスが良いと思われます．

Element.f90

! 同様にモジュール（クラス）を定義していきます． 
module Element_
use Node_
implicit none
public :: Element

type Element
    private
    real(8) :: p, m(4, 4)
    type(Node), pointer :: node1, node2, node3, node4

contains
    procedure :: init => Element_init
    procedure :: link => Element_link
end type

contains

subroutine Element_init(this, p, m)
    class(Element) :: this
    real(8),intent(in) :: p, m(4, 4)

    this%p = p
    this%m = m
end subroutine

subroutine Element_link(this, node1, node2, node3, node4)
    class(Element) :: this
    type(Node),target  :: node1, node2, node3, node4

    this%node1 => node1
    this%node2 => node2
    this%node3 => node3
    this%node4 => node4
end subroutine

end module

値の初期化（この例ではsampleという意味のない値を用いていますが，実際にはちゃんと意味のある値を入れましょう）と，構造体の互いの関係をそれぞれサブルーチンinitとlinkで定義します．

Area.f90

! 同様に(ry
module Area_
use Node_
use Element_
implicit none

type Area
    integer :: num_nodes, num_elems
    type(Node),    allocatable :: nodes(:)
    type(Element), allocatable :: elems(:)
contains
    procedure :: init => Area_init
    procedure :: link => Area_link
end type

contains

subroutine Area_init(this, num_nodes, num_elems)
    class(Area) :: this
    integer, intent(in) :: num_nodes, num_elems
    integer :: i
    real(8) :: sample1 = 1d0, sample2(2) = 2d0, sample44(4, 4) = 4d0

    this%num_nodes = num_nodes
    this%num_elems = num_elems

    allocate(this%nodes(num_nodes))
    do i = 1, this%num_nodes
        call this%nodes(i)%init(sample2, sample2)
    end do

    allocate(this%elems(num_elems))
    do i = 1, this%num_elems
        call this%elems(i)%init(sample1, sample44)
    end do
end subroutine

subroutine Area_link(this)
    class(Area) :: this
    integer     :: i, j

    do i = 1, this%num_elems
        j = 3 * (i - 1)
        call this%elems(i)%link(this%nodes(j+1), this%nodes(j+2), this%nodes(j+3), this%nodes(j+4))
    end do
end subroutine

end module

これで図に示す構造体の関係はできました．実際の使い方は以下の様になると思います．

main.f90

program main
    use Area_
    implicit none

    integer :: num_nodes = 7
    integer :: num_elems = 2

    type(Area) :: this_area

    call this_area%init(num_nodes, num_elems)
    call this_area%link()

end program

これで出来たと思います．どうでしょうか．

更新

2020/02/29　クラスごとにファイルを分けるような記載に変更．