MPIを用いた時間領域差分法による大規模3次元地震波動場の並列計算
MPIを用いた時間領域差分法による大規模3次元
地震波動場の並列計算
竹中博士* 藤井雄士郎**
1 はじめに
最近の計算機の進歩のお蔭で,現実的な地下構造モデルに対する地震波伝播の数値シミュ
レーションが身近なものとなりつつあり,高度な地震動予測手法として期待されている.地震
波動の数値シミュレーション手法には,大きく有限要素法や差分法などの領域型解法と境界要
素法等の境界型解法の2種類がある.領域型解法は計算領域を有限化しその領域全体を維散化
して解く手法で,境界型解法に比べて精度が劣る反面,複雑な地下構造モデルにも対応できる
という利点がある.中でも有限要素法や差分法は既にある程度確立された手法として広く用
いられている.一方境界型解法は,境界積分方程式を構成してそれを離散化して解く手法で,
無限に広がった領域を簡単に扱えるという利点を持っており,例えば地震発生場に特有の非常
に多くのクッラクが分布したような媒質中での地震波の散乱問題等を解くのに適している.し
かし,積分方程式を離散化した行列は密行列となり,大規模な問題を扱う際のネックとなって
いる.このような地震波動の数値計算手法について関心のある方は,例えばレビュー[1]. [2],
同をご覧いただきたい.実際的な地下構造モデルのもとで行う地震動予測は地震災害軽減の
ための1つの有効な手段である.予測は全国各地で行う必要があるため,その計算は高速の
pcまたはワークステーションクラスの計算機で行えることが望ましい.しかし,詳細な3次
元地下構造モデルを用いて地震波動シミュレーションを行うためは,莫大な計算メモリーと計
算時間が必要であり,現在のワークステーション単体(単一CPU)の能力では非常に難しい.
そこで我々は,ワークステーションクラスの計算機をクラスターとして用いることによって,
大規模な3次元の計算を行うことができる並列計算コードの開発を目指している.本稿では,
最近我々が開発した差分法の並列計算コードについて紹介する.
2 差分法コード
2.1 支配方程式
地震波動の支配方程式は弾性体の運動方程式と構成方程式(Hookeの法則)である.この支
配方程式のフォーミュレーションには,速度一応力型と変位一応力型の2種類がある.後者
は,差分計算の際に応力をストアする必要がないので,前者に比べて計算メモリーを軽減でき
*九州大学大学院理学研究院地球惑星科学部n **九州大学大学院理学研究科地球惑星科学専攻
E-mail: {takenaka.fuj ii}￠geo.kyushu-u. ac. jp
-85-
九州大学情報基盤センター広報
Vol. 2 No. 1 2002
研 究 開 発
る.本稿の差分コードは,次の変位-応力型の支配方程式に基づいている.
<
いⅣm肌u
+
砦 普
+
+
II
'0-∼.
+
nⅣmmり
用いlmmり
5i
ん
+
,.Ii:I,LJ.
+
+
i
L
・
i
;
,
+
Tr. r. -
(1)
(A+2p砦+A(普+砦-Mr
(A+2,,^+A 有 ・砦-M,yy
(A+2M)^+A 有 ・箸 M,
p(箸+普-M.xy>
p(普+ ) -Mx
du3
Tyy -
ou3
T7, =
'xy -
∂uz
TT, =
Tyz =
dx
+
∂uz
ay
)
- M,yz'
(2)
(1)式が弾性体の運動方程式,
(2)式が等方弾性体の構成方程式である・ここで, uiと7b.
(i,j-x,y,z)は変位と応力, fiとMr ま体積力とモーメント・テンソル, pは密度,A, IL
はLamァの定数である.地震は断層が滑ることによって生じると考えられるが,その滑りと等
価な力源をftまたはMiの項に代入することによって地震波を励起することができる.
2.2 スタガード格子差分法
解法には空間4次精度(一部2次)時間2次精度のスタガード(くいちがい)格子差分法
(時間領城陽解法)を使用する.図1が計算に使用するスタガード格子である.変位と応力の
位置が空間と時間について半格子ずつずらして配置されている.このとき,例えば(1)の最
初の方程式は以下のように離散化できる.
t・,霊 - u霊
+ △t・b [Dxtxx +DyTxy+Dztxz+fx]￨芸3,k
u豊 - ux- +△t・血
(3)
ここで,密度の代わりその逆数である比容H-l/p)を用いている.また転は速度を表す.式
中の添字は, u豊Lを例にとって説明すると,空間座標aAx,j△y,kAz),時刻(n+1)△tを
九州大学情報基盤センター広報
Vol. 2 No. 1 2002
-86-
MPIを用いた時間領域差分法による大規模3次元地震波動場の並列計算
y=G-i/2)A i=1 i=2 i=nx
△ △ A " A
k=1 ▼
△
k=2 ▼
。 △
群:
▼uy,b,QロT , n
≡=
A T ,^
■
。 △
k=nz▼
PE'
i=1 i=2 1=nX
y=jAy
● ∇ ● ∇ ● ● ∇
k=1 °
●
k=2 0
. ●
拝 三=
ux, b,Q"
uz, b, Q"
0 Tn,T r,X;
九,汁
∵
. ●
∇ T M^
k=nz°
Gi-∵鄭
l l 日・-トt
n=1 n=2 n=nt o ux,uy,u2
^xx>・yy>TzzITxyIT粉T yz
図1:スタガード格子の配置
表している.また式中のD。,Dy,Dzは中心差分演算子で,例えば規則格子の場合
Dxux時,j,k
1
Ax
[a(<霊-?/n+1 1
・ C2<n+1 -ォ)]
(4)
ただし,
cl-芸tc2--義(4次)
( C¥-¥,C2-0 (2次)
である.我々のコードではx,y,z方向に不規則格子([4])を採用しており,計算メモリーの
軽減と計算時間の短縮を図っている.境界条件には,上端面に自由表面条件(【5]) ,底面及
び側面に吸収境界条件([6】, [7])を用いている.
-87-
九州大学情報基盤センター広報
Vol.2 No. 1 2002
研 究 開 発
図2:計算領域分割とデータ通信
3 並列化
大規模な3次元シミュレーションの実用化において計算コードの並列化は非常に重要な課題
である.並列計算では分割した計算領域を複数のCPUに割り当て,各小領域の計算を独立に
進行させる.そのため計算メモリーの分散と計算時間の短縮において非常に有効である.今回
の差分法の計算では,図2に示すように計算領域をZ方向に分割し,それぞれの小領域を異な
るプロセッサーに割り当てる.この領域分割によって分断されたデータ(値)を用いて計算を
行う必要もあるため,それぞれの小領域には袖領域を用意し,そこに小領域をまたぐ計算に必
要なデータをプロセッサ一間通信によって転送する.プロセッサー間通信にはMPI (Message
Passing Interface)通信ライブラリ[8])を用いる.使用言語はFORTRAN90である.
4 計算例とパフォーマンス
我々は,このコードを用いて1995年兵庫県南部地震の地震動シミュレーションを実施してい
る([9コ).本稿では,このシミュレーションを情報基盤センターのスカラー並列マシンGS320
で行い,コードの並列化効率について調べた.
まず, 1995年兵庫県南部地震のシミュレーションについて簡単に紹介する.明石海峡付近
を震源とする気象庁マグニチュード7.3のこの地震は,死者6千数百人,全半壊家屋2 0万
棟に達する戦後最悪の地震災害「阪神・淡路大震災」を引き起こした.被害が特に著しかった
地域は神戸市須磨区から西宮市にかけての幅1 km長さ20 kmに渡る狭い帯のように集中して
おり,「震災の帯」と呼ばれている.我々は,この被害に最も大きく寄与していると考えられて
九州大学情報基盤センター広報
Vol. 2 No. 1 2002
-88-
MPIを用いた時間領域差分法による大規模3次元地震波動場の並列計算
135-00'E 135 1O'E 135 20'E
34。 50'N
34。 40'N
34。 30'N
図3:計算領域
いる周期1秒前後の地面の動きを上記の差分コードを用いて再現した[9]).
図3の四角で囲んだ領域が解析対象領域で,これについて深さ24 kmまでを計算領域とし
た.図中の星印が震央で,破壊の開始点(震源)の位置を地表に投影した点である.断層の
すべりを時空間的に記述する震源過程のモデルには,周期1秒前後の地震動を対象にした松
島・川瀬のモデル [10]を採用した.神戸地域の地盤構造モデルは, 2000年10月の時点
で我々の手に入る資料を参考にして作成した.
図4は震源の破壊開始後7秒, 10秒と13秒のスナップショットである.対象領域の上側で
左から右に延びているカラー(左から赤,緑,育)の直線は,鉛直な震源断層を地表に投影し
た線で,その右側の紫色の四角は断層面が若干傾斜している部分で震源断層の最東部分に当
たる.破壊は断層の左端から右に向かって進行しいる.左から右に向かって進行している赤や
緑の帯が地震波の波面である.図5は,計算で求まった周期1秒前後の地動(断層に直交する
速度成分)の最大値をカラーのコンターで表示したものである,赤色が濃い程揺れの強さが大
きかったことを表している.濃い赤色が帯状に分布しており,ちょうど震度7の震災の帯(図
中の赤線で囲った島状の部分)に対応していることがわかる.
-89-
九州大学情報基盤センター広報
Vol. 2 No. 1 2002
研 究 開 発
0.0
図4:波動伝播の様子
九州大学情報基盤センター広報
Vol. 2 No. 1 2002
90-
MPIを用いた時間領域差分法による大規模3次元地震波動場の並列計算
50 1 00 1 50 200
図5:最大速度分布と震度7の領域
表1計算パラメータ
この計算に用いたパラメーターを表1に示す.全体の計算時間はGS320の12PE(CPU)の場
合で約38時間であった.計算メモリーは12PEの場合(袖領域を含めて)で約4.5GBである。
このコードの並列実行の状況を見るために13PE利用時のUpshot (MPIプロファイリングツー
ル)の出力を図6に示す.図の(a)はジョブの開始直後からのUpshot, (b)はMPIによるデー
タ通信部分の拡大図である.黄色い部分はpEが同期を取るために待機している時間,緑色の
部分がMP工によるデータ通信(本コードの場合はsendrecv)を行っている時間を表示してい
る.ジョブの開始直後pE5以降のPEがしばらく待機しているのは PE4までが担当する地盤
浅部の物性値をファイルから入力しているからである.
図7は並列化効率を見るために,利用するpEの数を変えてジョブを実行した際の実行時間
をプロットしたものである.この測定のために実行した時間ステップ数は1500である.利用
-91-
九州大学情報基盤センター広報
Vol. 2 No. 1 2002
研 究 開 発
ih)
図6: Upshot
九州大学情報基盤センター広報
Vol. 2 No. 1 2002
- 92-
」
MPIを用いた時間領域差分法による大規模3次元地震波動場の並列計算
60000
55000
50000
巳iコ
鞄5000
◆
q)
.∈40000
ト
◆
7535000
=I
o30000
E≠
25000
◆◆◆
20000
◆
1 5000
2 4 6 8 10 12 14 16 18
Number of PEs
図7:プロセッサ数と計算時間
するpEの数が増えると当然計算時間が減る傾向にあるが,減少する割合はpEの数にしたがっ
て鈍っている. 9PEのとき8PE利用の場合よりも計算時間が増加する逆転現象が見られるが,
これはGS320特有のアーキテクチャ(利用するpEの配置)と関係しているかもしれない.図
7からこのシミュレーションをGS320で実行する場合, 12PEを超えて使用してもほとんど効
果がないことがわかる.
ただ,この結果はあくまで今回行ったシミュレーション程度のサイズの問題を解く場合にの
みあてはまることに注意して頂きたい.問題のサイズがこれより小さい場合には12PEとい
う上限値はさらに小さくなるであろうし,問題のサイズが大きくなると12PEを超えても十分
な効果が期待できるであろう.
謝辞
本研究は,林田智宏氏(現在,長崎県立佐世保南高等学校) ,川瀬 博教授(九州大学大学院
人間環境学研究院) ,岩田知孝博士(京都大学防災研究所)との共同研究の成果の一部である.
差分の計算には情報基盤センターの無料アカウントを利用させて頂きました.その際,藤野
清次教授と南里豪志助教授にはたい-んお世話になりました.一部図の作成にはGMTマッピ
ング・ツール([11])を使用しました.
参考文献
[1]緒楳-起・竹中博士:近地地震波の伝播に関する理論,地震, 42, 391-403, 1989.
[2]竹中博士:不整形地盤における波動伝播の数値計算法,地震, 46, 191-205, 1993.
93-
九州大学情報基盤センター広報
Vol. 2 No. 1 2002
研 究 開 発
[3]物理探査学会編「物理探査ハンドブック」 (1999年物理探査学会出版)第16章シミュ
レーション.
[4] Pitarka, A. :3D elastic finite-difference modeling of seismic motion using
staggered grids with nonuniform spacing. Bull. Seism. Soc. Am., 89, 54-68,
1999.
[5]林田智宏・竹中博士・岡元太郎:速度-応力型スタガード格子差分法を用いた2次元及
び3次元地震波動計算コードの作成,九大理学部研究報告(地球惑星) , 20, 99--110,
1999.
[6] Clayton, R. and Engquist, B. : Absorbing boundary conditions for acoustic
and elastic wave equations, Bull. Seism. Soc. Am., 67, 1529-1540, 1977.
[7] Cerjan, G., Kosloff, D., Kosloff, R. and Reshef, M.: A non-reflecting bouncU
condition for discrete acoustic and elastic wave equations, Geophysics, 50,
705-708, 1985.
[8]情報基盤センターのWebページhttp: //www. cc.kyushu-u.ac.jp/scp/system/library/
MPL/MPI.htmlを参照.
[9]林田智宏・竹中博士・藤井雄士郎・川瀬 博・岩田知孝 1995年兵庫県南部地震「や
や短周期」地震動の再現,地球惑星科学関連学会2001年合同大会予稿集(CD-ROM) ,
S5-POO2, 2001.
[10]松島信一,川瀬 博:1995年兵庫県南部地震の複数アスペリティモデルの提案とそれに
よる強振動シミュレーション,日本建築学会構造系論文集,第534号 33-40, 2000.
[11] Wessel, P. and Smith, W.H.F.: New, improved version of the Generic Mapping
Tools released. EOS Trans. AGU, 79, 579, 1998.
九州大学情報基盤センター広報
Vo.2 No. 1 2002
-94-
© Copyright 2024 ExpyDoc
