Spectral Element Method (SEM)を使って計算してみる(その1)

ちょっとした波動場の計算をやってみたいという欲望があるのですが、今更いちから差分法のコードをかく気力もなくて横着はできないか、と考えていたところ、某所にてSEMを使った計算結果を見せてもらって、ソースコードも公開されている、という話を聞きつけました。それで、理屈はまったく勉強せずに、とりあえず使えるかどうかを試してみよう、と思い立ちました。

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軸方向の速度のアニメーションです。

x-axis.anim.gif

拡張子が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でデータが並んでいるということならば簡単に読み出せそうです。

結構、親切なマニュアルがついている割には意外にもこういうところの情報をマニュアルで見つけられなかったりしてはまります。

いずれにしても、コンパイルして計算をして、結果を見る、ということはなんとかできそうだ、ということがわかりました。

posted by MOR at 18:00Comment(0)Computer

FreeBSD 11.2Rでアプリケーションサーバー (3)

アプリケーションサーバーをFreeBSd 9.1Rから完全に移行するつもりではなかったので,Xクライアントからのログインの設定をする必要はないだろうと考えていました。しかし,いったん新しいOSが使えて新しいバージョンのアプリケーションが使えるということがわかれば,ユーザーは誰も9.1Rなんて使わないだろう,と予想されるので念のためにxdmも動くようにしておきます。

xdmをpkgからインストールしてあとは設定です。また,サーバー側にもちゃんとwindow managerをインストールするのを忘れないようにしなくてはなりません。これを忘れてログインできるのにすぐにログアウトしてしまう,というワナにはまりまくってしまいました。これまでずっとicewmを使っていてかつ,標準には含まれていないテーマを使い続けているので,pkgでicewmをインストールしたあと,/usr/local/share/icewm/themes/の下にテーマ(私の場合はsortofaqua)の設定ファイル一式をコピーしておきます。

/usr/local/etc/X11/xdm/にある設定ファイルに手を入れます。*.sampleというファイルがありますが,これはデフォルトの設定が書かれています。失敗したら,このファイルをコピーすれば少なくとも初期状態に戻すことはできますので,安心して設定ファイルを書き換えることができます。

xdm-configの一番下の行のDisplayManager.requestPort: 0というのを!でコメントアウト。

Xserversの一番下の行の:0 local /usr/local/bin/X :0を#でコメントアウト。

Xaccessの*だけの行と* CHOOSE BROADCASTと書かれた行、一番下のLISTEN * ff01:0....と書かれている行のコメントをすべてはずす。

/usr/local/etc/rc.d/xdm.shというファイルを作って、
#!/bin/sh
/usr/local/bin/xdm
と書いてchmod +xで実行可能にしておくと起動時にxdmが起動してリモートからX -query HOST_NAMEでログインできるようになります。
posted by MOR at 19:00Comment(0)Computer

FreeBSD 11.2Rでアプリケーションサーバー (2)

OSが動くようになったので次はアプリケーションのインストールです。ユーザーのリクエストにできるだけ応えてかつ,手間をかけない,という方針で何も考えずにpkgからバイナリをガンガンインストールしてしまいます。portsからコンパイルなんて面倒なことはもうやめちゃいます。

最初にOSのバージョンを最新のものにしておきます。

freebsd-update fetch
freebsd-update install

最初にpkgを使うとインストールされてないけどインストールするかい?と聞かれるのでyesと答えてpkgをインストールします。あとは

pkg install XXXXX

として延々とアプリケーションをインストールしていきます。

インストールした主なアプリケーションは以下の通りです。

converter:
dosunix, base64, iconv, mpack, unix2dos

archivers:
gtar, gzip, lhasa, untar, unzip, zip, sharutils

ftp:
wput, wget

japanese:
ja-a2ps, ja-kon2-14dot, ja-less, mtools, ja-nkf, ja-w3m, ja-vim, xpdf

mail:
postfix

dns:
bind-tools

x11:
xorg, xdm

x11-wm:
icewm

www:
opera, firefox, chromium

editors:
emacs

print:
yatex, texlive-full, gv

graphics:
gimp, ipe, gmt5, pgplot, xfig, epdfview

devel:
fpp, gmake, pycharm-ce, py36-virtualenv, py36-pip, py36-pip-tools, ruby24-gems

lang:
ruby24, python3, gcc8, perl5

net:
rsync

math:
R, graph, gnuplot, scilab, octave, octave-forge, maxima, wxMaxima, scalapack, linpack, laspack, paraview

multimedia:
xanim

ローカルでX windowを使わないのであればwindow managerはいらないだろう,と思っていたのですが,リモートホストからxdmを使ってdirect queryでログインしようとするとサーバー側にもwindow managerがいるので,icewmをインストールしています。これに気がつかなくて,リモートからログインできるにもかかわらず,すぐに強制的にログアウトしてしまうというわけのわからんワナにはまりました。考えれば当たり前のことなのですが。

ざっとこんなところです。pythonのanacondaはportsもpkgもなくてインストールしていません。pythonにはいろいろな仮想環境があるようですが,とりあえずvirtualenvでいいんじゃないかと勝手に決め打ちでインストールできるもので済ませてしまっています。


posted by MOR at 05:00Comment(0)Computer