SEMは1998年にKomatitsch and Vilotteによって提案された手法で(BSSA, Vol.88, No.2, pp.368-392)、FEMの一種だろう、と勝手に想像していますが、いかんせんまったく勉強していないのでよくわかりません。勉強はプログラムが使えそうだということがわかってからちゃんとやることにします。公開されているSEMのコードは、de Martinによってかかれたもので(BSSA, Vol.101, No.6, pp.2855-2865)、ここで公開されています。Downloadのページをみると2016年8月にリリースされた1.322というバージョンのソースが置かれています。
これをありがたくダウンロードさせていただいて、コンパイルをしてみます。
インストール方法などについても丁寧なマニュアルがあるのでそのとおりにやればできるのですが、それなりにハマりどころがあるため、忘れないようにメモをしておきます。
実行するために必要な環境は、
1. intelのparallel studio for Linux
2. intelのMPI Library
3. CUBIT mesh generator
4. METIS-5.1.0 library
ということになっていますが、2のMPIの環境はOpenMPIでも問題なく動くようです(私の使っている環境にはintelのMPIはインストールされていないはずだから)。
私が普段使っているクラスタマシンには1の環境が整っているはずなので、使えないことはなさそうです。
1. ファイルの展開
ソースファイルをダウンロードしたらファイルを展開して、そのディレクトリに移動します。
tar zxvf EFISPEC3D.tgz
cd EFISPEC3D/
2. make.configの編集
マニュアルではいきなりmakeをしていますが、コンパイラの指定が異なっているので、config/make.configをまず編集します。
vi config/make.config
として、コンパイラを指定している変数、FC, CCの中身を書き換えます。
#FC = mpiifort
FC = mpif90
#CC = mpiicc
CC = mpicc
もともと書かれていた内容を#でコメントアウトして環境にあわせて書き直します。私の所では上記のようにします。書き換えたら保存して、viを終了します。
3. METIS-5.1.0のコンパイル
EFISPEC3Dをコンパイルする前にMETISをコンパイルします。METISのサイトからソースをもらってきて展開します。
gunzip metis-5.1.0.tar.gz
tar xvf metis-5.1.0
cd metis-5.1.0
intelコンパイラを使うように環境変数を設定してコンパイルします。
export CC=icc
export CXX=icpc
make config CC=icc CXX=icpc
make
./build/Linux-x86_64/libmetis/ディレクトリの下にlibmetis.aができているので、これをEFISPEC3Dのlibディレクトリにコピーします。
cp ./build/Linux-x86_64/libmetis/libmetis.a (somewhere)/EFISPEC3D/lib/
もともと、EFISPEC3D/lib/にlibmetis.aが入っているので64bit環境であればわざわざコンパイルしてコピーする必要はないのかもしれません。しかし、私の環境で作ったlibmetis.aのほうが少しファイルサイズが大きくなっているので何かが違うのかもしれません。実際のところよくわかりません。
4. make
EFISPEC3Dをコンパイルします。基本的にはデフォルト設定のまま
make efispec
とやるだけでよいはずですが、私の環境では、SIMDのオプションをSSEの後継のAVXにしたほうがよいかもしれません。手元のクラスタは既に年季がはいっていてXeon E5-2670の2-way 10 nodeで構成されています。はるか昔のSandy Bridge世代のCPUです。あとから5 nodeを買い足したのでもう少し新しいCPUのものもあるのですが、MPIを使うのに異なる実行ファイルは使えませんから古いものにあわせるしかありません。
EFISPEC3Dのマニュアルではavx2も指定できると書いてありますが、私のところではCPUが対応していないので(Haswell世代のCPUからの対応)、avxを指定します。
make efispec version=1.0 simd=avx
コンパイルに成功すると、bin/efispec3d_1.0_avx.exeができています。実行ファイルはこれだけです。
libmetis.aは先にコンパイルしましたが、他のライブラリ類はlibディレクトリにアーカイバによる静的リンクのためのライブラリ(lib*.a)がパッケージについています。このライブラリそのものをコンパイルするときのオプションがどういうものかまったくわからないので自分でこのライブラリを用意するのはちょっと無理な予感です。兎に角、ついてきたものを黙って使うのがよさそうです。
5. 動作確認
実行ファイルができればあとは実行するだけなのですが、設定ファイルの作り方などはかなりたいへんです。ちゃんと動作するかどうかをtutorialを使って確認します。
cp -r docs/tutorials/t01 test
cd test
としてチュートリアルのt01の設定ファイルを丸ごとコピーしてきます。実行のためのコマンドは環境に依存しますが、私のところでは以下のようにしました。
cat > node.txt
Node01
Node02
^D
としてMPIで使うノードをリストアップしたファイルを作成しておきます。そのうえでmpirunで実行します。
mpirun -np N --machinefile node.txt ../bin/efispec3d_1.0_avx.exe
-np NのNの値は計算に用いるCPUコアの数を書けばよいようです。ノードを二つ使うなら私の環境では最大で8 cores x 2 ways x 2 nodes = 32 coresまで指定できることになります。オーバーヘッドなどを考えるとそれが本当に適切なのかどうかはよくわかりませんが...。
ところがこれを実行するとエラーで死にます。理由は粘性減衰を考慮する、という設定でコンパイルしたのに、設定ファイルはそれを考慮していないからのようです。そこで、粘性減衰を考慮しないで線形の応力-ひずみ関係を使うように設定を変更します。また、出力されるスナップショットをGMTのグリッドファイルではなくpara viewのvtk形式で出力するように変更しておきます。
そのためには、module_global_variables.f90のなかの変数の設定を変えます。
cd ../
vi src/efispec3d_1.0/module_global_variables.f90
として、
!>if .true. --> activate viscoelasticity in the solver module: see mod_solver
!! logical, parameter :: LG_VISCO = .true.
logical, parameter :: LG_VISCO = .false.
!>if .true. --> write snapshots in VTK format
!! logical, parameter :: LG_SNAPSHOT_VTK = .false.
logical, parameter :: LG_SNAPSHOT_VTK = .true.
!>if .true. --> write snapshots in gmt format
!! logical, parameter :: LG_SNAPSHOT_GMT = .true.
logical, parameter :: LG_SNAPSHOT_GMT = .false.
のように書き換えます。これを保存して、改めてコンパイルしてから実行します。
make efispec version=1.0 simd=avx
cd test
mpirun -np N --machinefile node.txt ../bin/efispec3d_1.0_avx.exe
今度は問題なく実行できます。
6. 出力
主な出力ファイルは以下のとおりです。
t01.lst：種々のログ
t01.fsr.*.gpl：レシーバーの一での波形です。gnuplotのバイナリ形式ということですが、どんな形式なのかちょっと謎です。
t01.collection.surface.snapshot.vxyz.pvd：para viewのためのスナップショットのファイルをひとまとめにしたxmlファイルです。しかし、なぜかvtsファイルのファイル名を間違って出力するのでpara viewで読み込んでもエラーになります。次に示すvtsファイルをまとめて読み込んでしまえばアニメーションを見ることはできるのでこのファイルが使えなくてもあまり問題はありません。
t01.snapshot.vxyz.*.vts：para viewのためのスナップショット(速度の三成分)が出力されています。
GMTのグリッドファイルでスナップショットを出力する場合はgrdという拡張子になります。
t01.snapshot.PGXx.grd：XはA (加速度), V (速度), D (変位)のいずれか、xはx, y, zのいずれかで座標軸です。GMTのグリッドファイル形式です。(たぶん)地表面における各ノードの最大値を出力したものです。
t01.snapshot.PGXxyz.vts：上記と同じですが、x, y, z成分をまとめて最大値を出力したもので、para viewのvtk形式です。
7. 図化
7.1 スナップショット
拡張子がvtsのスナップショットはpara viewで見ることができます。番号がついているファイルはグループとしてまとめて読み込むことができるのでアニメーションを見ることができます。
para viewで読み込んで，View -> Animation Viewと選ぶと別窓が開くのでModeを選びます。Snap To TimeStepsにしておくと，全部のスナップショットを使ってアニメーションが出力されます。さらに，Filg -> Save Animation...としてjpgで保存すると連番のファイル名でスナップショットの数だけjpgファイルが保存されます。
出力されたjpgファイルをanimation gifに変換します。Imagemagickを使って，
convert -loop 0 hoge.*.jpeg hoge_anim.gif
というふうにします。*の部分がファイル名についている連番を全部使うということを意味しています。
出来上がりは，こんな感じです。ファイル名とは裏腹にY軸方向の速度のアニメーションです。
拡張子がgplの波形は図化するのがちょっとたいへんです。バイナリファイルですが、バイナリの形式がよくわからないのです。tutorialを読むと
a=2
plot "t03.fsr.000001.gpl" binary format="%10f" u 1:a with line
というようにすればよいことになっています。aの値は、
a=2,3,4のとき、x, y, z方向の変位
a=5,6,7のとき、x, y, z方向の速度
a=8,9,10のとき、x, y, z方向の加速度
です。
ところが、gnuplotはこんなformatは知らない、といって絵を描いてくれません。かわりに、
plot "t03.fsr.000001.gpl" binary format="%float%float%float%float%float%float%float%float%float%float" u 1:a with line
というふうに%floatを10回繰り返して書けばそれらしい絵を書いてくれます。「10回繰り返す」ことを簡単に書く方法がよくわからないので格好が悪いですがとりあえず、これで済ませてしまいます。
計算結果をみてみると、t03.fsr.00000*.gplのファイルサイズが26680 byteでした。t01.cfgに計算の設定がかかれていますが、
シミュレーション波形の継続時間：20秒
シミュレーションの時間刻み：0.003秒
波形の出力の時間刻み：0.03秒
となっています。したがって、1つの波形は2668 byteで一つのデータは4 byteで時間ステップは667ステップ分、ということのようです。どうやら単精度(4 byte = 16 bit)で出力されているようです。あとで波形をみるときにデータを直接読めないと困ってしまいますが、たんに4byteでデータが並んでいるということならば簡単に読み出せそうです。
結構、親切なマニュアルがついている割には意外にもこういうところの情報をマニュアルで見つけられなかったりしてはまります。
いずれにしても、コンパイルして計算をして、結果を見る、ということはなんとかできそうだ、ということがわかりました。