やりたいこと

SPECFEM3Dのプログラムを利用し、地震波動シミュレーションをする。
より詳しくいうと、
自分で作った弾性体内でCMT解がわかっている地震波を発生させて、伝播する様子をみる。
ということです。

基本的には、ドキュメント通りですが、いろいろ他に参照しないとわからないことも多いので、ここでまとめておきます。

本節は、ドキュメントを一読し、シミュレーション方法を整理するフェーズにある方を対象としています。

あると良いツール

推奨されているツールは、CUBIT,ParaViewです。
それぞれの役割は、以下です。
（＊）CUBITで、弾性体と物性値、配置、位置を設定し、メッシュを作成。
（＊）ParaViewで、シミュレーション結果を見る。

Getting Started

ドキュメント2章にあたる。
ここは割愛。クローンして環境を整えるのみ。

open-mpi, gfortranなどをインストールしておくこと。

(2024/1/24 追記)古い設定のままにしている場合、

brew upgrade gfortran

をしないとmakeできない。

操作の流れ

基本的にやることは決まっていますが、簡単なシミュレーションを実行するための手順は以下です。

1. CUBITで、弾性体と物性値、配置、位置を設定し、メッシュを作成する
2.Par_fileを設定
3.CMTSOLUTION,STATIONSを設定
4.run_this_example.shを実行
5.ParaViewで結果を可視化する

CUBITでメッシュを作成

CUBITは、3Dモデルをつくったりメッシュを描いたりするツールです。ここでは超最低限のコマンドをまとめます。

CUBITにはコマンドラインが各種ありますが、基本は"Command"と、"Python"のコマンドラインを使います。（ErrorとHistoryは、エラーログと、コマンド履歴をみられます。）

PythonでCUBITのコマンドを使うときは、cubit.cmd("コマンド")とします。たとえば、resetコマンドをpythonで実行するときは、

cubit.cmd("reset")

とします。また、pythonでcubitコマンドを使うときでも、cubitをimportしなくても使えます！

コマンド	意味	補足
reset	設定の全てをリセット
delete vol 1	volumeID 1のボリュームを除去
brick x 10	一辺の長さが10の、立方体を作る
brick x 10 y 11 z 12	x方向に長さ10,yに11,zに12の直方体を作る
vol 1 scale x 1.5 y 2	volumeID 1のボリュームをx軸方向に1.5倍、y軸方向に2倍する
move vol 1 x 3	volumeID 1の立体を、x方向に3動かす
webcut vol 1 with plane zplane rotate 45 about x	volumeID 1の立体を、z方向を向く平面(y-z平面)をx軸を中心に45度回転(x軸に関して右手系)させた面で切断
imprint all	隣接した2つのボリューム間の共有面を作成	(ⅰ)
merge all	隣接した2つのボリュームの共有面を結合	(ⅰ)
mesh vol all	全ボリュームのメッシュを作成	(ⅰ)
reset vol all	全ボリュームの情報をリセット（meshなど）
block 1 hex in volume 1	volume 1の六面体要素を"block 1"に含める
block 1 name 'elastic 1'	"block 1"の名前を、'elastic 1'に変更
block 1 attribute count 6	"block 1"に、6つの属性枠を作成
block 1 attribute index 3 1500	"block 1"の3つめの属性値を1500に設定	(ⅱ)
block 2 face with z_coord < -1	z座標が−１より小さい部分を、"block 2"とする	(ⅲ)
quality volume all shape global draw mesh	メッシュのクオリティを可視化	(ⅳ)

コマンドの補足

(ⅰ)
２つのボリュームを作成し、imprint,mergeをせずにmesh vol allをすると、２つのボリュームそれぞれのメッシュを作成してしまいます。つまり、１つのボリュームとしてではなく、２つの独立したボリュームがそれぞれでメッシュを作ってしまいます。
また、接地面（共有面）を定義しないままだと、2枚の面がくっついていることになり、期待通りにいきません。複数の弾性体をつくって多層構造にしたいときは、imprint,mergeを行ってからmeshするといいでしょう。

(ⅱ)
このコマンドにより、blockに密度や剛性率、Vp,Vsを設定します。CUBIT上では密度などの意味を持ちませんが、このあとSPECFEM3Dに持ち込んだとき、その意味を持つようになりす。ここで、属性値の順番に注意します。
6つの属性値は順に、material_id,Vp,Vs rho(密度),Q_k, Q_muです。
Q_k, Q_muは、bulk attenuation, shear attenuationで、減衰係数です。

ちなみに、blockはbrickとは異なります。blockはCUBIT上の固有名詞的なもので、ひとまとまりのボリュームを意味しています。特定の形をもつ立体のことではないです。brickは、"レンガ"などの意味をもつ英単語で、CUBITでは直方体を指します。

(ⅲ)
このコマンドにより、弾性体の吸収境界の部分を定義します。6面体の弾性体を扱うなら、6面の吸収境界を作っておきましょう。（さもないと、吸収境界を使ったシミュレーションが実行できないはずです。）

(ⅳ)
これは結構重要で、メッシュがどのような品質になっているか見られます。
例を見ましょう。

sphere radius 1100 #半径1100の球体を作成
mesh vol all #メッシュを作成

実行すると次のようになります。
（Fig1）

ここで、メッシュのクオリティを見てみます。わかりやすいように、断面を見てみます。

quality volume all shape global draw mesh

(Fig2)

各メッシュの体積の大きさが違うことがわかります。メッシュの大きさが大きく違うと、シミュレーションの計算時間に大きな影響があることがあり得ます。できるだけ等しい体積でメッシュできるようにするといいでしょう。CUBITにはメッシュの方法が様々ありますから、適宜使い分けると良いでしょう。

cubitのコマンドについては、やりたいことを英語で調べるとなんとなく出てきますが、情報が少なく、公式ドキュメントが最も参考になります。(ChatGPTに聞いても、正確な情報が出てこない 2024/1/24)

メッシュのエクスポート

作成したメッシュはエクスポートし、シミュレーションに使えるようにします。
CUBITのコマンドラインで、MESHディレクトリを作成します

#結果を入れるディレクトリを作成
mkdir "MESH"

その後、pythonで、次を実行します。
ただし実行する前に、CUBITのメニューから、set directoryを選択し、
/Applications/CUBITWorkingDirectory
などとして、CUBITのディレクトリがどこにあるか、設定してください。
そして、SPECFEM3Dのディレクトリから、geocubitlibを探してコピーし、/Applications/CUBITWorkingDirectory以下において、pythonで呼び出せるようにしてください。

sys.path.append("/Applications/CUBITWorkingDirectory/geocubitlib")

import geocubitlib.cubit2specfem3d
geocubitlib.cubit2specfem3d.export2SPECFEM3D("MESH")

その後、MESHディレクトリに複数のファイルができると思いますので、MESHファイルごとSPECFEM3Dのルートディレクトリに移動します。これでメッシュは完了です。

Par_fileの設定

基本的には、パラメータの意味は、ドキュメントのcreating_databasesの章(4章)の通りですので、詳しくはそちらを参照してください。
ただし、次の設定に注意しましょう。

CREATE_SHAKEMAP                 = .true. #shakemapを作成するか
MOVIE_SURFACE                   = .false. #地表での波動伝播の動画を作成するか
MOVIE_TYPE                      = 2　#弾性体の面全体(2に設定)、topのみ（１に設定）
MOVIE_VOLUME                    = .true.　#内部での波動伝播動画を作成するか
SAVE_DISPLACEMENT               = .true.　#変位を保存するか
MOVIE_VOLUME_STRESS             = .true.　#応力(圧力)を保存するか
...
USE_FORCE_POINT_SOURCE          = .false. #点力源を使うか
...
SAVE_SEISMOGRAMS_DISPLACEMENT   = .true.　#観測点での変位のグラフをとるか

などです。

点力源を使うかどうかは、CMT解を使うかどうかに依存します。基本的にはCMT解を使うと思います。モーメントテンソルを任意に設定して実体波を放射できますので、爆発震源や逆断層型などを入力するといいでしょう。

英単語、surfaceですが、普通は”表面”という意味で使われます。Earth's surfaceは地表ですので、このような場合、単なるsurfaceは、モデルの側面やbottomも含まれていることに注意しましょう。

CMTSOLUTION, STATIONSを設定

CMT解は、好きなように設定できますが、実際のCMT解を使ってみるのが、おそらくいいでしょう。
このフォーマットは、Global CMT Web Pageからコピーできます。
モーメントテンソルや震源位置は好きに設定できますが、time shiftは0でなければなりません。
（ただし、サブイベントを設定するときは、複数のイベントのうち、少なくとも一つのtime shiftは0であればいいようです。）

MLI 1976  1  1  1 29 39.60 -28.6100 -177.6400  59.0 6.2 0.0 KERMADEC ISLANDS REGION      
event name:       010176A  
time shift:       0.0000
half duration:    5 
latorUTM:       67000  
longorUTM:      67000   
depth:          30 
Mrr:      -7.600000e+27
Mtt:       7.700000e+27
Mpp:     -2.000000e+26
Mrt:      -2.500000e+28
Mrp:      4.000000e+26
Mtp:      -2.500000e+27

STATIONSは、観測点の位置です。
X31やDBは文字数制限がありますが、好きに設定できます。観測点に名前とIDをつけておくだけなので、シミュレーション結果には影響しません。
その他の値は、詳細はドキュメントを参照してください。

X31 DB 67000.0 35892.86 0 0.0 0.0
X55 DB 67000.0 64607.14 0 0.0 0.0

ドキュメントにはありませんが、CMTSOLUTIONの震源位置と、STATIONSにある観測点位置はモデルの中になければなりません。震源や観測点を、作成したボリュームの外におくことはできません。

run_this_example.shを実行

./run_this_example.sh

かならず、ルートディレクトリで実行してください。

場合によっては、

chmod +x run_this_example.sh

で実行権限を与えてください。

実行すると、シミュレーションが動くと思います。

おすすめの追加ファイル

シミュレーションが動くと、/OUTPUT_FILES/DATABASES_MPI 以下に、"proc000001_displ_X_it000200.bin"というようなファイルが大量にできます。
これをひとまとめにして動画を作成するわけですが、これをいちいちやるのが結構面倒なので、以下のようなbashファイルを、ルートディレクトリに置いておき、シミュレーション後にこれを実行するといいです。
このファイルを実行すると、一気に全部の動画を作れます。
（本来は、ドキュメントにあるように、個別に./bin/xcombine_vol_data_vtu　を使うのですが、これが結構面倒。）

#!/bin/bash

# ./DATA/Par_fileからNTSTEP_BETWEEN_FRAMESとNSTEPの値を取得
INCREMENT=$(grep -i "NTSTEP_BETWEEN_FRAMES" ./DATA/Par_file | awk '{print $3}')
END=$(grep -i "NSTEP" ./DATA/Par_file | awk '{print $3}')

# 開始値を設定
START=$INCREMENT

# ループでコマンドを実行

for (( i=$START; i<=$END; i+=INCREMENT )); do
  # ゼロ埋めされた値を生成
  NUMBER=$(printf "%06d" $i)
  # コマンドを実行
  ./bin/xcombine_vol_data_vtu 0 3 displ_X_it${NUMBER} ./OUTPUT_FILES/DATABASES_MPI ./OUTPUT_FILES 0
done

for (( i=$START; i<=$END; i+=INCREMENT )); do
  # ゼロ埋めされた値を生成
  NUMBER=$(printf "%06d" $i)
  # コマンドを実行
  ./bin/xcombine_vol_data_vtu 0 3 displ_Y_it${NUMBER} ./OUTPUT_FILES/DATABASES_MPI ./OUTPUT_FILES 0
done

for (( i=$START; i<=$END; i+=INCREMENT )); do
  # ゼロ埋めされた値を生成
  NUMBER=$(printf "%06d" $i)
  # コマンドを実行
  ./bin/xcombine_vol_data_vtu 0 3 displ_Z_it${NUMBER} ./OUTPUT_FILES/DATABASES_MPI ./OUTPUT_FILES 0
done

これは、変位の動画displ_...がある場合にのみ有効です。また、Par_fileからNTSTEP_BETWEEN_FRAMESとNSTEPの値を取得するので、ファイルのあるディレクトリの場所などに注意してください。

ParaViewで、シミュレーション結果を見る

ParaViewを用いて、vtuファイルを時系列順に並べて動画にすることができます。shakemapもfile typeを適宜調べてみることができます。
動画の保存はavi形式などが一般的ですが、mp4は選択できません。各々ツールを使って変換することをおすすめします。

私のシミュレーション結果は以下です。Fig3、Fig4はそれぞれ、displ_Yの結果と、shakemapです。
CMT解のモーメントテンソルは、20110311の、東北地方太平洋沖地震のものを使っています。
CMT解を調べるには、Global CMT Web Pageから検索することが、ドキュメントの方法と一致します。

（Fig3）y方向の変位

（Fig4）shakemap(norm of displacement)

ところで、Fig3のブロック内部をよくみると、仕切りのような壁が見えています。これはプログラムを実行する際に、mpiで並列計算させた時に4つのプロセッサで計算したことを暗示しています。これはCUBITで与えたものではないし、計算にまったく影響しないものです。

ParaViewのAnimationをgifに出力する

imagemagickが一番はやそうです。以下でインストールが可能です。

brew install imagemagick

一応、アニメーション結果を保存するだけなら、ParaViewからavi形式などで動画にすることができますし、結果の情報を別の形式で保存できますが、軽量なままわかりやすく保存するためにはgifがいいと思ったので、メモしておきます。
わざわざメモするのは、ParaViewから直接gifで結果を保存できないからです。

Paraviewで、"file" → "Save Animation"を選び、保存先のディレクトリを作り、pngを選択します。ここでクオリティ選択などができますが、デフォルトでいいとおもいます。
inpやvtuなどは画像ファイルの連続なので、その枚数分の画像が保存先のファイルに出力されます。（このため、デスクトップにそのまま出力しないほうがいいです）

その後、imagemagickのコマンド

convert -delay 20 -loop 0 *.png animation.gif

でgifにできますが、コマンドの意味は、

convert <1フレームの表示時間> <時間> <ループ回数> <回数> <画像(群)> <gif名>

です。

-delay: 1フレームの表示時間
-loop: ループ回数. 0は無限ループ
*.pngは、カレントディレクトリ内のすべてのpngを対象

となっています。

*.pngをコマンド内に用いる場合は、カレントディレクトリ内にその画像があることと、対象の画像が正しい順番に並んでいることを確認してください。必要に応じて、cdコマンドで保存したディレクトリに移動してから実行してください。
ParaViewでpngで出力すると、アニメーション画像名が連番になるので、ディレクトリ内で画像を名前順に並べた後、上記のコマンドをそのまま実行すればいいと思います。

諸注意

CMTSOLUTIONやFORCESOLUTIONのパラメータなどを設定するとき、座標系のとり方には注意です。
参考書や論文などで、z軸の正の方向が下なのか上なのかは、一意(暗黙の了解)ではありません。
このSPECFEM3Dでは、

x軸正方向：東
y軸正方向：北
z軸正方向：上

方向とされています。Aki&Richardsの定義方法とは違うので、注意せよとのことです。

とりあえず、以上！
また追記する予定です。

（補足や間違いなど指摘していただける方は、ぜひ教えてもらえるとうれしいです。）