（様式 20）成果報告書

TSUBAME 共同利用

平成 26 年度

学術利用

成果報告書

不均質地球構造における地震波・津波伝播シミュレーション
Simulation of seismic/tsunami wave propagation in the heterogeneous earth

1,2

1

古村 孝志 ・前田 拓人
Takashi Furumura and Takuto Maeda
1

2

東京大学地震研究所 ・情報学環総合防災情報研究センター
1
2
Earthquake Research Institute /Center for Integrated Disaster Information Research ,The University of Tokyo
http://www.eri.u-tokyo.ac.jp/people/furumura
http://www.eri.u-tokyo.ac.jp/people/maeda

邦文抄録（300 字程度）
差分法による並列地震動シミュレーションコード SEISM3D に対し，OpenACC を用いた GPU 計算を実現し，その性能評価を行った．TSUBEME システムにおける複数の実装タイプに対する性能評価の結果，これまで一般に行われてきた CUDA を用いたコードの大幅な書き直しを必要とすることなく，OpenACC のみで最大 2.92 倍の加速を確認した．データ出力部分を除けばウィークスケーリング性能も良好であり，今後の大規模シミュレーションへの適用が期待される．
英文抄録（100 words 程度）
We implemented OpenACC-based GPU computing in SEISM3D, a parallelized
finite-difference-method simulation software of seismic wave propagation. Throughperformance measurements of different types of implementation on TSUBEME system, x2.92
speed up at maximum has been achieved only by introducing OpenACC and without significant
rewriting of code.

Keywords: 地震動, 波動伝播, 差分法, OpenACC

背景と目的における地震波動伝播過程そのものも重要な研究対

地震学・固体地球物理学においては地震波形記録象である．しかし，観測網の充実に比して，日本列が地球内部構造ならびに震源での断層破壊過程を知島スケールの空間規模において地震動の主たる周波るための基礎的なデータである．日本国内には数帯域（0­15 Hz 程度まで）をカバーするにはまだ1000 点を越す稠密な基盤的連続観測網が敷かれ，計算機能力は十分とは言えず，継続的な地震動シミ日々地震動記録を蓄積している．さらにこの観測網ュレーション手法ならびにコードの継続的な高度化は海域に向かって拡大を続けており，海底におけるが求められている．

巨大な地震・津波観測網も構築されつつある．大型

本課題で主として取り扱う地震動シミュレーショ

計算機による地震動シミュレーションは，これら観ンソフトウェアは隣接ノード通信を伴う等間隔格子測量を解析・説明するための基礎的なツールであり，の 差 分 法 コ ー ド (e.g., Maeda and Furumura,観測波形と計算波形との比較や同化によってより詳2013)であり，地球シミュレータ(Furumura and細な地球内部構造や震源の破壊過程像が得られるよChen, 2005) や 京 コ ン ピ ュ ー タ (Maeda et al.,

うになると期待されている．加えて，地震波動は不2013)をはじめとした種々のプラットホーム上での均質媒質中における波動伝播の物理学そのもの重要動作実績を持つ．コードは各アーキテクチャ上でチな一フィールドでもあり，不均質な地球内部構造中ューニングされており，京コンピュータではピーク

（様式 20）成果報告書
性能比 20％を超えるなど大規模並列計算機上で十GPU 環境双方で統一的なコードを用いることが出分な性能を達成してきた（井上・他, 2013）．しか来るようになる可能性がある．

し，シミュレーションが密媒質中のベクトル弾性波

そこで，本報告では現在利用されている地震動コである以上必然的に多くの 3 次元変数を必要とし，ードの OpenACC 化を試み，OpenACC による GPUかつ 3 次元的に不均質な地球内部構造も記憶参照さを用いた地震動シミュレーションの実現可能性につれる必要があるため，アプリケーションは全体としいての検討を行った．また，TSUBAME をはじめとnodesroc_xnodes

node (1,4)

node (1,3)

node (1,2)

node (2, 1)

node (1, 1)

node (4, 1)

node (3, 1)

node (Ip, 1)

(b)

np

node (Ip, Jp)

nproc_ynode (1, Jp)(a)node (idx+1, idy)

MPI

MPI

node (idx, idy)
MPI

xyz

com

mu

nic

node (idx+1,idy+1)

MPI

MPI

node (idx, idy+1)

xatio

ns

ns

MPI

nicatiocommuy図 1. SEISM3D の並列化レイアウト (Modified from Maeda et al., 2013)てメモリインテンシブなコードである．今後は計算した複数環境における CPU からの性能向上の比較

能力に対するメモリからのデータ供給能力が低下しと将来的な大規模 GPU 計算の可能性の検討を行っていくと想定されるが，GPU 等の加速器を利用するた．

ことで大きく効率を改善できる可能性がある．実際に，すでに TSUBAME を含む CUDA 環境を用いて概要

大規模な地震動シミュレーションを実現した例が複本課題で用いる地震動シミュレーションコード

数あり（e.g., Okamoto et al., 2010）地震動シミ

SEISM3D (after Furumura and Chen, 2005)は，ュレーションへの GPU 利用の有効性が示されてい地球内部構造を粘弾性体近似のもとに差分法で解くる．しかし，アプリケーションコードが CUDA を並列コードである． 粘弾性緩和は地震動の周波数用いた事実上の書き直しになることが GPU を利用

帯域では地震波エネルギーの散逸（内部減衰）をするうえでの大きな困難となっていた．

表現する．地震動シミュレーションにおいては簡

近年，ハード面では GPU 内蔵のメモリ量が地便のため完全弾性体に一定のダンピングを用い震動シミュレーションに十分耐えうる程度に増加しることが広く行われてきたが，本コードは特徴的たこと，ソフト面では OpenACC が標準化され普及時間の異なる複数の粘弾性体モデルを並列させされつつあることによって，既存コードを GPU 対て扱うことで，広い周波数帯域において地球内部応化するための敷居が大きく下がってきた．

構造の内部減衰特性を良く再現したシミュレー

OpenACC は対応環境以外では無視される指示子ションが実現可能になった[Maeda et al., 2013]．

によって記述できるため，既存コードからの移植が解くべき運動方程式は，たとえば x 方向についてきわめて容易であると期待されるほか，CPU と（様式 20）成果報告書∂vx ( x, t ) ∂σ xx ∂σ xy ∂σ xz=++∂t∂x∂y∂zと表される．ここで ρ は質量密度， vx は弾性体変位速度のｘ成分， σ ij は応力テンソルであり，等ρ (x)また，計算対象として地球構造の一部分を取り出してシミュレーションするため，計算範囲の外周周辺には Perfectly Matched Layer 法（Chewand Liu, 1996）の吸収境界条件をおいた．これ

方粘弾性体では 6 つの独立な成分を持つ．応力テは時間発展方程式を空間微分の方向毎に人工的ンソルは構成方程式により変位あるいは速度とに分解し，境界面の方向に直交する方向のみ地震結びつけられる．粘弾性体の場合には，この構成波振幅を強減衰させることで，吸収境界からの人工反射波をきわめて小さくする手法である．その一方，内部領域とは異なる時間発展方程式を解き，かつ要求メモリが大きくなるという欠点を持つ．
並列化は 3 次元空間を鉛直 2 次元に分割し，分割された領域の境界周辺グリッドを時間積分サイクル毎に隣接ノードと交換することによって行う（図１）．この計算は多くの変数が 3 次元配列であり演算量の割にメモリ入出力が多く，メモリバンド幅が性能を律速するコードである．

結果および考察

図 2. NVVP によるプロファイル結果の例．

これらの問題を解く既存のコードにOpenACC の指示子を追加することで GPU 用コ方程式は粘弾性の緩和時間関数時間方向のたたードの開発を行った．SEISM3D は地球シミュレ

み込みになるが，そのままでは数値計算上不利なータや「京」コンピュータ等に代表される様々なため，メモリ変数(Robertsson, et al., 1994)とアーキテクチャで開発ならびに実行されてきたいう手法により，これを複数の時間発展方程式にため，マシンの特徴にあわせた様々な実装がある．

分割する．本コードで採用されている粘弾性モデ

そこで，それらの実装の特徴を抽出した 2 種類の

ル Generalized Zener Body の場合について，xy実装を OpenACC 化して比較した．

成分は
NN
*
∂σ xy ( x, t )1 M # τ εSm &-# ∂vx ∂vy & M= 2µ R ,1−+%1−(/%( + ∑ rxym (t )∑∂t+ N M m=1 $ τ σ m '.$ ∂y ∂x ' m=1

2 種類のコードは時間発展方程式右辺の取り扱いが異なっている．一つ目のコード（コード A）
では，独立な 3 重ループで前述の方程式の右辺に# τ S &# ∂v ∂v &∂rxym ( x, t )1 1=µ R %1− εm (% x + y (∂tN M τ σ m $ τ σ m '$ ∂y ∂x 'ある空間微分を評価し，それらを 3 次元配列に一のように表される．ここで， N M は並列に用いる重ループで時間積分を計算する．もう一つ（コー粘弾性体モデルの個数，µ R は粘弾性体の relaxedド B）では，空間微分の計算と時間積分を融合しmoduli で媒質による定数， τ は粘弾性の緩和時ている．この場合，右辺を格納する 3 次元の一時間とクリープ時間で，粘弾性による地震波の内部配列が必要なくなり，利用メモリ量が減る．この減衰とその周波数特性を支配する定数である．こことによりメモリ空間の小さなクラスタマシン

れらの方程式の右辺を食い違い格子配置の空間 4でも比較的大規模なシミュレーションが実現可次精度の中心差分で評価し，時間 2 次精度の時間能であり，かつ Intel CPU での実行効率はコード積分を行うことで，速度ならびに応力の時間発展A よりも高い．しかし，メモリ削減の代償として問題を交互に解く．

ループ内の演算が非常に複雑となり，レジスタ溢度セーブする．その結果を再度ロードして別の 3

（様式 20）成果報告書
れを起こしやすい．結果的に「京」コンピュータ
や地球シミュレータ 2 においてはコード A よりも性能が落ちるという特徴がある．

タビリティには課題を残す結果となった．
計算実行性能の比較のため，まず最初に TSBAMEの 2 ノードを用いた試験を行った．CPU での実行

OpenACC 化は Fortran90 コードのループひとは 4CPU24 コアを Flat-MPI で用い，GPU の場合つひとつについて OpenACC の指示子を追加するには 4GPU 一つ一つに CPU の 1 コアを割り当てた．

ことで行った．本コードで必要な通信は入出力をのモデルサイズは両者共通に 384x256x256 グリッ

ぞけば隣接ノード間のみであり，ほとんどのデータドを用い，1000 時間ステップの計算を行った．両は CUDA 上に保持したまま，通信に必要な袖領域者の計算実行時間の比較を図３に示す．
いずれも CPU に比べて相応の性能向上が確認された．しかし，性能向上率には実装の仕方による顕著な差異が見られた．3 次元一時配列を用いたコード A では 2.92 倍の加速率を達成したのに比べ，省メモリのコード B では 2.14 版にとどまった．コー
ド B ではループが融合し複雑化した結果，利用レジスタ数が多くなったため加速率に限界が表れたものと考えられる．コード B は当初はより加速率低い状態であったが，がコンパイルオプションの調整による最大レジスタサイズの指定，再内レベルでのループの分割，OpenACC の vector 指示子のサイズ調

図 3. CPU 実行および OpenACC 利用の GPU実行の計算時間の比較．

整などのチューニングを施した結果，2 倍を超える加速率を達成した．
次に，コード B を対象として，GPU 計算のウィークスケーリングテストを実施した．ノードあたりのみを CPU に転送して MPI 通信を行うことができた．このことにより，GPU とホストとの間の通信を最小化し，効率的な計算が可能になった．GPU-CPU間の通信が MPI 通信およびデータ出力のタイミン
グでのみ発生し，その他の時間帯は GPU のみで計算が進んでいることを NVIDIA Visual Profiler から確認した（図 2）．
ただし，OpenACC 化されたコードは指示子を除いて完全に CPU 版と同一にすることはできなかった ． MPI 通 信 デ ー タ の パ ッ キ ン グ に 用 い て い た
Fortran90 の一部の配列関数が OpenACC コンパイラで非対応であり，該当部分の全面的な書き換えが必要であったほか，性能向上のため，配列のインデックスの順序変更，一部変数の事前計算と 3 次元配列へのストアが必要であった．特に後者のチュー
ニングは CPU 環境での高効率計算と必ずしも整合せず， OpenACC 利用の利点であるコードのポーのグリッド数を 64x64x512，計算時間ステップを1000 ステップに固定し，3GPU から 384GPU までの評価を実施した．結果を図 4 に示す．
出力がない場合スケーリング性能は良好であり，GPU が特に少ない場合ではスケーリングからのずれが見られるが，これは吸収境界条件が外周部周辺で常に一定のグリッド数を占めるため，特に問題規模が小さすぎる場合には計算量がノード数と完全には比例しないためと考えられる．それ以上，24GPU以上の場合にほぼ一定の計算時間で計算規模を拡大することができている．しかし，データの出力がスケーリング性能に顕著な影響を及ぼすことも明らかになった．
SEISM3D は一定時間ステップ毎に，3 次元空間中の 2 次元断面の変数列をバイナリデータとしてファイル出力している．出力されるデータは複数のMPI ノードにまたがるため，MPI の集合通信を用い（様式 20）成果報告書てこれを特定ノードにまず集約した上でファイルを時間の隠蔽を行う必要がある．また，MPI 通信に出力する．GPU 計算の場合には，まず出力データバついても隠蔽を行うことでさらなる高速化が可ッファを GPU から CPU に転送し，さらに CPU 間能となるであろう．本報告ではもっぱら地震動シで MPI 通信がなされるため，通信時間が大きなオーミュレーションを対象にしたが，津波計算の基礎バーヘッドになっている物と考えられる．これは今方程式も地震動シミュレーションと似た構造を

後 MPI 通信ならびにデータ出力の GPU 計算とのオしており，同様の OpenACC 化を今後行うことでーバラップすることによる隠蔽を実装することで大さらなる加速が期待される．

きく改善する可能性がある．
謝辞
まとめ，今後の課題

OpenACC 化に当たり，NVIDIA 社の成瀬彰氏に貴重なアドバイスをいただきました．記して謝意を表OpenACC Weak Scaling, 64x64x512 grids/nodeします．

100

Computation Time [s]

8060

4020
without output
with output

0100101102
Number of GPUs103

図 4. ウィークスケーリングテストの結果．青四角がデータ出力無し，赤丸がデータ出力ありの場合の計算時間をそれぞれ示す．
地震動シミュレーションコード SEISM3D をOpenACC を用いた GPU 化を行った．CPU とGPU で完全に同一コードとすることまでは実現できなかったが，従来一般に行われてきたようなCUDA を用いた大幅な書き直しをすることなく，OpenACC のみで最大 2.92 倍のコード加速を実現できた．多少の性能の差異はあるが，OpenACC化されたコードはそのまま他の CPU 環境でも実行可能であり，さまざまな計算資源を有効に活用したシミュレーションとそれを用いた固体地球物理学の研究への活用が期待される．また，ウィークスケーリング性能も良好であり，今後の大規模シミュレーションへの適用が充分に可能であると考えられる．ただし，そのためには今後出力