初めに
初投稿です.拙い部分しかないと思いますが,まだまだ初学者ですゆえどうかご了承ください.
この記事はfortranでスパコンなどの大規模並列計算機を使うシミュレーションプログラムを作成するためのヒントとなればいいなあという記事です.ですが,まず初めにこれからシミュレーションプログラムを作成しようと思われる方はこの記事を読む前にこれとこれを買って読んでください.読むべきです.話はこれからです.シミュレーションにおける,まさにバイブルとも呼べる本です(個人の感想です).私は本当にこの本と同じ時代に生まれたことに感謝しています.シミュレーションでなにをすべきか,全てがここに書いてあります.
ではなぜいまさら私ごときがこんな記事を書くかというと,確かに本の中には「これはC++でもfortranでも同じである」とは書いてあるのですが,中身はC++で書かれてるんですよねー.しかしオブジェクト指向に関するfortranの記事は少ないという….
というわけで,まずは内容で一番最初に躓いた,構造体とポインタについて書いてみようと思いました.いまのところ以下のやり方で特に不具合はありませんが,もしかしたら問題があるかもしれません.詳しい方は是非コメント欄に指摘ください.全力で勉強させて頂きます.
データモデルの表記法
本の中では有限要素法を例に,メモリレイアウトについて言及されています.以下の図のようなデータ構造を保持したいとします.
米印はポインタであることを表しています.私は有限要素法に明るくないのですが,それでも全てのノードでループを回したいとき,エレメントごとにループを回したいときがそれぞれあると思うのでこのメモリ配置がベストなんだと思います.
fortranのコード
あとは上図を実現するコードを書きましょう.
! クラス名にアンダーバーをつけてモジュール名を定義.
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
他言語で言われているように,1モジュール1ソースコードが良さそうです.また,fortranは言語仕様上1モジュールに複数クラスを定義できますが,管理が破綻しそうなので1モジュール1クラスが良いと思われます.
! 同様にモジュール(クラス)を定義していきます.
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
で定義します.
! 同様に(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
これで図に示す構造体の関係はできました.実際の使い方は以下の様になると思います.
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 クラスごとにファイルを分けるような記載に変更.