今度は，チュートリアルと同じ2層地盤に挑戦してみます，って言うほどたいへんな話じゃないはずなのですが，よくわからんのでいろいろハマります。
1. モデル作成
何がたいへんかといって，メッシュ切りがたいへんです。チュートリアルでは表層を800mメッシュ，基盤層を2400mメッシュで切っています。しかし，Gmshを使ってボリュームごとに異なるサイズを切る方法がよくわかりません。もちろん，四面体要素でよければGmshはそれらしく切ってくれますが，EFISPECのほうがそれに対応していませんからいかんともしがたいものがあります。
難しいことを考えるのはやめて，まずは同じサイズの六面体要素を切ることにしました。先の計算では，地盤モデルである直方体を構成する全ての点の座標を入れて，それをもとに体積を構成する，というもっとも原始的な方法を使いました。しかし，このような簡単なモデルの場合，一つの平面を定義してそれを特定の方向へ引き延ばす(スキャン? Extrudeと呼ぶようです)ことで三次元モデルに拡張することがよく行われます。今度はそのような方法を試してみることにします。
1.1 ポイントの作成
まず，Extrudeするもとになる面を構成します。その面を構成する点を作成します。xz平面上に面を構成してその面をy方向にExtrudeしようという目論見です。そのため，xz平面上の点を定義します。
Point(1) = {-29000, -29000, 0, 800.0};
Point(2) = {29000, -29000, 0, 800.0};
Point(3) = {29000, -29000, -34000.0, 800.0};
Point(4) = {-29000, -29000, -34000.0, 800.0};
Point(5) = {-29000, -29000, -1000, 800.0};
Point(6) = {29000, -29000, -1000, 800.0};
前の計算との違いはPoint 5と6が増えていることです。この二つの点を繋いだラインが表層と基盤面の間の境界面になります。
1.2 ラインの作成
面を構成するときには，反時計回りが基本なので，それを意識してラインをつないでいきます。Line 2が境界面になる予定です。
Line(1) = {1, 5};
Line(2) = {5, 6};
Line(3) = {6, 2};
Line(4) = {2, 1};
Line(5) = {5, 4};
Line(6) = {4, 3};
Line(7) = {3, 6};
1.3 ラインループの作成
反時計回りにラインループを構成します。xz平面だけ考えているので，地盤モデルのxz断面だけ作ればOKです。Line Loop 21が表層，Line Loop 22が基盤層です。Line Loop 22を構成するLine 2がマイナスになっているのはLine Loop 22で反時計回りにループを構成するにはLine 2の向きが反対なので左向きのベクトルとなるようにするためです。
Line Loop(21) = {1, 2, 3, 4};
Line Loop(22) = {5, 6, 7, -2};
1.4 面の作成
作成したラインループを用いて面を作成します。これはラインループをそのまま面として定義すれば終わりです。
Plane Surface(31) = {21};
Plane Surface(32) = {22};
1.5 面から体積へ拡張
ここまででxz平面上の地盤モデルの断面が完成したので，これをy方向に引き伸ばして体積を作成します。このためにExtrudeを使います。
Extrude{0,58000,0}{
Surface{31};
Surface{32};
}
Extrudeのあとのカッコの中はどの方向に面を引き延ばすかを指定しています。この例では，y方向に58000だけ引き延ばす，ということを示しています。
このようにして二つの体積をもつ直方体が完成します。面を引き延ばすことで自動的に体積が作られるため面ループや体積を改めて定義する必要はありません。
これはこれで結構なのですが，Extrudeによって作成されたポイントやライン，面，体積の名前(番号)がよくわかりません。
なんらかの番号が自動的につけられているのですが，どういう規則で番号をつけているかわからないので，これは結構，ハマります。なぜなら，番号がわからないとこの後のステップで必要になるTransfiniteができないからです。
しょうがないのでGUI上で番号をコツコツ調べます。
ここまでで作成したgeoファイルをダブルクリックするとGmshが起動して，何か描かれているのでマウスで向きを適当に変えるとそれらしい直方体になっていることがわかります。この状態で，マウスをポイントやラインの上に載せるとポイントやラインの情報がポップアップで表示されます。数が多くなると覚えきれないので，ノートにメモしておかないと大混乱になります。
ところが，面や体積はマウスを載せても何も表示されません。本当はメニューのどこかを使って面や体積の番号を表示する方法があるのかもしれませんが，よくわからないので場当たり的に対応します。左ペインのメニューでModule -> Physical groups -> Add -> Surfaceと選んでいくと面を表す点線が表示されます。その点線の上にマウスを載せると面番号がポップアップで表示されます。同様にVolumeを選ぶと黄色い球が表示されて(この場合は2つ)，この球の上にマウスを載せると体積番号が表示されます。体積は安直に1, 2という番号が振られていますが，面についてはパッと見ただけではよくわからない脈略のない番号がついています。
このようにして調べて見ると以下のようになっていることがわかりました。丸，四角，三角の数字はそれぞれ，ポイント，ライン，面の番号を表しています。ライン番号の横の数字はラインのベクトルの向きを示しています。こうやって改めて見ても番号の振り方にはなんだか脈略がありません。
1.6 面のtransfiniteの設定
Extrudeによって作成された面をメッシュに切るときに四角形で切られるようにtransfiniteを設定します。面番号は前のステップで調べているのでその番号を指定すればOKです。
Transfinite Surface {31};
Transfinite Surface {32};
Transfinite Surface {49};
Transfinite Surface {71};
Transfinite Surface {54};
Transfinite Surface {76};
Transfinite Surface {41};
Transfinite Surface {63};
Transfinite Surface {53};
Transfinite Surface {45};
Transfinite Surface {67};
さらに，メッシュが四角形となるように，recombineします。
Recombine Surface{31,32,49,71,54,76,41,63,53,45,67};
1.7 体積のtransfiniteの設定
最後に体積をtransfiniteして六面体要素が構成されるようにします。六面体メッシュを切るためには，面がそれぞれ4辺からなっていて，かつ体積が6面からなっていなくてはならない，という制約条件があります。そのうえで，transfiniteするに際に表面を構成するポイントの番号を反時計回りに指定してから裏面を構成するポイントの番号を反時計回りに指定しなくてはなりません。
Transfinite Volume{1} = {1,5,6,2,7,8,12,16};
Transfinite Volume{2} = {5,4,3,6,8,18,22,12};
2. メッシュ切り
geoファイルが正しく書けていればメッシュ切りは簡単です。geoファイルをダブルクリックしてGmshでモデルを表示させ，左ペインのメニューからModules -> Mesh -> 3Dと選べばメッシュを切ってくれます。これをAbaquasのinp形式でエクスポートします。
3. inpファイルの修正
基本的に先の計算の場合と同じです。ただ，表層と基盤層の間の境界面の取り扱いに注意が必要です。どうやらEFISPECは境界面としては，吸収境界と自由表面以外はまったく興味がないようなのでそれ以外の面を定義することはできません。
そのため，層境界であるSurface 45は丸ごと削除します。
また，自由表面はこの場合はSurface 53なのでこの面だけELSET=fsuを指定します。それ以外の面は吸収境界としてELSET=prxを一回だけ指定して重複する*ELEMENT行は削除します。
ELSET=Lineを指定しているELEMENTも全て丸ごと削除します。
ELSET=Volume1はELSET=l01に，ELSET=Volume2はELSET=l02に書き換えます。
typeを大文字のTYPEに書き換えて，面はS4R，体積はC3D8Rに修正して，体積(l01, l02)，吸収境界(prx)，自由表面(fsu)の順にELEMENTの定義を並べ替えます。
あとは，ヘッダを少し書きなおして，最後の行に*END PARTを追加すれば完成です。
4. 計算の実行
あとはチュートリアルと同じ条件で計算をするだけです。チュートリアルのt01に置かれているcfgファイル(設定全般)，fsrファイル(地表面のレシーバーの座標)，dcsファイル(ダブルカップルのソースの設定)をコピーしてきてprefixを適切に書きなおして(例えば，t03とだけ書いておく)，さらに，node.txtに計算に用いる計算ノードのホスト名を以下のように列挙します。
Node01
Node03
Node05
Node06
この場合，4 nodes x 2 ways x 8 coresで最大64コアが使えるので，設定ファイルがあるディレクトリで，
mpirun -np 64 --machinefile node.txt ../bin/efispec3d_1.0_avx.exe
とすれば64コアを使って計算をしてくれます(たぶん)。計算に関する情報は計算を行なっているディレクトリにlstという拡張子のファイルとしてログが残ります。このファイルを見ると，以下のようなことがわかります。
チュートリアルでは25,344要素でしたが，私の計算例ではメッシュを表層のサイズにあわせて全部同じサイズで切ったので，240,944要素もありました。かなりアホっぽいですが，練習なのでよしとします。
計算時間は同じように64コアを使ってチュートリアルはトータルで36.5秒，細かいメッシュを使った私の計算例では250.2秒でした。要素数が10倍弱で計算時間は7倍程度なので，計算時間が短い問題を解く場合，並列化のオーバーヘッドが馬鹿にならない，ということを示しているようです。ですのでむやみにコアをたくさん使ってもあまり賢いとは言えないかもしれません。
5. 計算結果
それでも24万要素の計算が4分あまりで終わってしまうなんていうのは遠い昔を知る者としては驚くばかりです。計算を専門にしている人からすればバカみたいなレベルの話ですが，ユーザーとして使うだけの立場からすれば計算機のありがたみを強く感じます。
スナップショットをparaviewで読み込んでjpegのパラパラアニメを出力して，ImageMagicで
convert -loop 0 2-layer_y.0* 2-layer_y_anim.gif
としてgifアニメを出力したものが下のアニメです。
左がチュートリアルをそのまま計算したもの，右が細かいメッシュだけを使った私の計算例です。y方向の速度を表示しています。先の計算に比べてメッシュが細かいのでパッと見た限りでは両者はよくあっているように見えます。
