学際大規模情報基盤共同利用・共同研究拠点 平成 23 年度共同研究 最終報告書 2012 年 5 月

11-NA03

GPGPU の地震ハザード予測シミュレーションへの適応性評価
青井 真（独立行政法人 防災科学技術研究所）
概要

昨年度までに TSUBAME2.0 を使用した単純なモデルによる弱スケーリングの性能評価で格子数が数億規模の比較的大きなモデルまで高いスケーラビリティを得られることができた。本年度は、より高い精度で地震ハザード予測をめざし、より規模の大きなモデルによる超並列計算のパフォーマンス検証を行うと同時に、大量の計算結果を伴うシミュレーションの出力アルゴリズムに関して基礎的な検討を行った。

1

研究の目的と意義

近年の計算機環境の劇的な進歩の恩恵を受けたと兵庫県南部地震を契機に、地震に関する調査研

はいえ、大きな計算機リソース（CPU パワー及びメ究の成果が国民や防災を担当する機関に十分に伝モリ）を必要とし、実務等で一般的に使用可能な達され活用される体制になっていなかったという計算機では数日以上の計算となる事もしばしばで課題意識の下に、行政施策に直結すべき地震に関ある。このような状況を打開する手段として、多する調査研究の責任体制を明らかにし、これを政くの数値計算分野で広く使われるようになりつつ府として一元的に推進するため、地震防災対策特あ る GPGPU （ General Purpose Computation on

別措置法に基づき総理府に設置(現・文部科学省にGraphics Processing Unit）の利用が考えられる。

設置された)機関として、地震調査研究推進本部が本研究では、GPU コンピューティングをスーパー設置されている。地震調査研究推進本部の地震調コンピュータシステムに積極的に取り入れている

査委員会において、地震に関する観測、測量、調TSUBAME 2.0 で地震波伝播シミュレーションを試行

査または研究を行う関係行政機関、大学等の調査し、高度な地震ハザード・リスク評価への適応性結果を収集、整理、分析し、並びにこれに基づきを評価することを目的とする。

総合的な評価を行っており、この中で地震動予測

2地図の高度化を進めている。日本周辺で発生する当拠点公募型共同研究として実施した意義

大地震に関して、地震のリスク評価の基礎となり地震の被害を軽減するためには、個々人の地震得る精度で地震ハザードを予測できるよう、手への防災意識を高め、地震に対する備えを促すこ法・モデルの高度化を目指している。そのために、とが不可欠である。このため、防災科学技術研究必要な精度、分解能を持つ地盤構造の開発を行う

所では、日本全国で発生する地震を対象として、とともに、高精度かつ汎用性のある地震波伝播（強地震調査研究の成果の集大成である地震動予測地震動）シミュレーション手法の開発に関する研究

図を高度化し、地震ハザード・リスク評価に関すを行っている。

る研究を行うとともに、WebGIS 等の技術を用いて、地震波伝播シミュレーションにおいて詳細な 3

地震ハザード・リスク情報、地下構造データ等の次元地下構造を十分な精度で離散化し、短周期の関連情報を網羅的に提供可能な地震ハザード・リ地震波まで計算するためには細かな格子が必要でスク情報ステーションを構築している。
本研究では東京工業大学との超大規模数値計算あるため、実用的な計算においては格子数が数億から数十億に及ぶ規模のモデルを扱うことになる。 応用分野での共同研究を行なっており、これによ1

学際大規模情報基盤共同利用・共同研究拠点 平成 23 年度共同研究 最終報告書 2012 年 5 月

り地震ハザード・リスク評価に用いられる地震波伝播シミュレーションに対する GPGPU に適用について、実務等に使用される汎用計算機環境からTSUBAME 2.0 のような大規模なスーパーコンピュータシステムまでの様々な規模の計算機環境における評価が可能となり、日本の防災力向上に大きく貢献することが期待される。

3

研究成果の詳細と当初計画の達成状況

3.1

GMS による波動場の計算

本研究では、不連続な食い違い格子（Aoi andFujiwara, 1999）を用いた、空間四次・時間二次

図 1 （中央）計算に用いる不連続格子。（右）精度の差分演算子による実用コードである GMS

不連続格子の垂直断面。領域 I と領域 II の接続（Ground Motion Simulator、青井・他、2004）を部分で、内挿のために格子が重なっている。
（左）

ベースに GPGPU の適用性の評価を行なった。GMS

領域 II の最上面（A 面）と領域 I の最下面（Dは、防災科学技術研究所によってパッケージ化さ面）における不連続格子の水平断面。

れた、3 次元有限差分法(FDM)により地震波伝搬シ3.2

ミュレーションを行うためのツール群でおり、主GPU での差分計算処理

GMS の差分法の計算処理部分については、ほぼに Fortran90 で書かれた差分計算ソルバはソース全て GPU 上で処理するように実装を行った。開発コードも公開されている。
Aoi and Fujiwara（1999）や青井・他（2004）

環境には、NVIDIA 社から提供されている CUDAは、大きさの異なる格子を組み合わせることによ（Compute Unified Device Architecture）を使用り効率的かつ高精度に計算を行うことの出来る不している。CUDA では、演算処理を担当する GPU と、連続格子による差分法の定式化を提案した。差分その演算処理に必要なデータを格納するビデオメ法による地震波伝播シミュレーションを行う場合、 モリ（グローバルメモリ）により構成されるデバ
格子サイズは計算すべき最短波長により決定されイスをアクセラレータとして用いることを前提とるため、モデルのごく一部のみが低速度の媒質でしている(図 2 左)。

ある場合でも計算領域全体を小さな格子に分割せ差分法の計算を行う際は、スレッド（GPU で計算ざるを得ず、大規模なモデル計算の大きな障害と処理を行う際の実行の最小単位）の集合であるブなっていた。軟弱（＝地震波速度が低い）な表層

ロックを 2 次元で構成し、計算の対象となる x-yが存在する浅い部分（領域 I）の格子点間隔は細平面にブロックを敷き詰めるように配置している。

かく、深い部分（領域 II）の格子点間隔は領域 I

ブロック内の各スレッドは z 方向の始点から終点の 3 倍の粗い格子点間隔を有する格子モデルを用まで 1 格子点ずつ計算処理を進め、3 次元領域全体いている（図 1）。2 つの領域はオーバーラップしの計算処理を行う（図 2 右）。差分法のように大量ており、波動場の連続性が保たれるよう内挿されのメモリアクセスを伴う計算ではメモリバンド幅る。典型的な盆地構造モデルの計算において、均がボトルネックになることが多いため、GPU 上に搭質な格子による場合と比較し数倍から十数倍程度載されている高速なレジスタやシェアードメモリ効率がよいことが分かっている。

をキャッシュのように使用して、低速なグローバルメモリへのアクセスを軽減する工夫を行っている。
2

学際大規模情報基盤共同利用・共同研究拠点 平成 23 年度共同研究 最終報告書 2012 年 5 月

ことで袖領域のみの計算を別途行う必要性を回避した（図 4）。なお、この方法は、不連続格子を用いていないアルゴリズムにも拡張する可能で、事前計算などの余分なオーバーヘッドが無いことから、通信を隠蔽する際には有効な手段と言える。

図2(左)GPU のアーキテクチャ、(右) GPU での計算処理

3.3

複数 GPU を用いた並列化

複数 GPU を用いた並列化においては、水平二方

図 3 複数 GPU での並列計算（水平方向の領域分向（x 方向、y 方向）の領域分割法を用いている。

割に並列化）

分割した各部分領域には、それぞれ 1 つの GPU を割り付けて計算処理を行う。各部分領域には、隣接する部分領域のデータを収める袖領域を設け、データが更新された際は、各部分領域間でMPI(Message-Passing Interface)ライブラリを用いてデータを交換する。
GPU の演算速度は極めて高速であるために相対的に通信に要する時間の割合が大きくなり、CPU のみを用いて並列化した際と比べて、並列性能が低くなる傾向がある。本研究では、これを改善させるため GPU での計算処理と通信をオーバーラップさせて、通信に要する時間の隠蔽を行っている。

図 4

一般に通信の隠蔽を行う際には、通信の対象とな（袖領域の事前計算が必要ない）

本研究で開発した通信時間の隠蔽方法

る袖領域の計算を事前に行い、内部領域を計算している間に並行して通信を行う手法がとられるこ

3.4とが多い。しかしこのような方法では、袖領域を並列計算の性能試験計算する際に GPU が苦手とする不連続なメモリア

GPU での実行に対応した差分計算ソルバについクセスが生じるため、袖領域の事前計算そのものて TSUBAME 2.0 を用いて性能試験を行った。使用が足かせとなってしまう場合がある（図 3）。本研したノードには GPU に NVIDIA Tesla M2050 x3、

究では、GMS が格子最適化手法として採用していCPU に Intel Xeon X5670 2.93GHz x2 が搭載される不連続格子が 2 つの異なる格子サイズを持つ領ており、ノード間は InfiniBand 4x QDR にて接続域（領域 I、領域 II）からなることに注目し、一されている。
まず、単体 GPU での性能を評価するため、GPU方の領域を計算する間に他方の領域の通信を行う
3

学際大規模情報基盤共同利用・共同研究拠点 平成 23 年度共同研究 最終報告書 2012 年 5 月

版プログラム、および比較用として CPU 版プログ超える）の性能を得られていることが分かる（図ラムの 1 並列での性能試験を実施した。計算モデ5）。一方、モデルを単位モデルに固定し使用する

ル に は 領 域 I が 420x420x100 、 領 域 II がGPU 数を増加させる強スケーリングに関しては、4140x140x200 の 21,560,000 格子からなるモデル

GPU では 3.2 倍の演算性能が得られるものの、16（Unit420）を使用し、浮動小数点演算は GPU 版プGPU では 7.3 倍の性能となり、効率が落ちてゆくログラム、CPU 版プログラム共に単精度となる。

ことが分かった。これは、並列数の増加に伴い各比較用の CPU 版プログラムの実行には Intel Xeon

GPU が担当する格子数が減少するため、相対的にX5670 の 1CPU コアを使用し、コンパイラには通信の時間が大きくなること、および GPU で生成Intel Compiler 11.1 を用いている。

されるスレッド数が減少し計算の効率が低下した

単体 GPU での性能試験の結果から、CPU 版プロことが要因として考えられる。通常行われるシミグラムに比べて GPU 版プログラムは 20.4 倍の性能ュレーションのタイムステップ数は数万ステップとなり、相当な高速化が図られたことが確認されであり、モデルサイズに応じた計算機資源が確保た（表 1）。また、NVIDIA CUDA に付属する性能解さ れ て い る 場 合 に は 計 算 に 要 す る 時 間 （ turn析ツールを用いて GPU-VRAM（グローバルメモリ）

around time）は数分から数十分程度となるため、間のメモリ転送速度を計測したところ、最大で理過剰な並列度による更なる高速化の必要性は実用論性能の 70%以上となることから、律速の要因が的にはほとんど無いと考えられる。

演算処理ではなくメモリ転送速度であることが確認された。計算の主要部分において、メモリ転送が理論性能に対して十分に高速に行われており、ハードウェアの性能を生かしていることから、更なる高速化の余地はそれほど大きくはないと考えられる。
表 1 単体 GPU での性能評価

図 5 並列数 - 高速化率の関係
続いて、複数 GPU を用いた並列性能の評価を実また、実行時間と浮動小数点演算数を用いて施した。計算モデルには、Unit420 を単位モデル

FLOPS（Floating-point Operations Per Second）とし、これを水平方向に 2x2、3x3、4x4、8x8、10x10、を算出した（図 6）。GPU 版プログラムでは、16x16、32x32 個並べて、全体の格子数がそれぞれUnit6720（約 55 億格子）の 256GPU による性能が4 倍、9 倍、16 倍、64 倍、100 倍、256 倍、102421.1TFLOPS、Unit13440（約 220 億格子）の 1024GPU倍となるようなモデルを使用した。

による性能が 79.7TFLOPS になることが確認され複数 GPU での性能試験の結果の結果から弱スケた。CPU 版プログラムでは Unit13440 の 1024 並列ーリングに関しては、ほぼ線形（実際には線形をでの性能が 2.3TFLOPS であることから、約 34 倍の4学際大規模情報基盤共同利用・共同研究拠点 平成 23 年度共同研究 最終報告書 2012 年 5 月

性能が得られていることが分かる。

する時間の割合が高くなってしまう。そこで、本研究では試験的に単体 GPU での実行に対応したプ
ログラムに対して、GPU での計算処理とファイル出力を同時に行う隠蔽アルゴリズムの開発を行った。GPU 側の処理は基本的に非同期に進行する。
そのため、GPU に計算処理を投入した後、前のタイムステップのファイル出力を行うことで、GPUでの計算処理とファイル出力をオーバーラップさせ、ファイル出力に要する時間を隠蔽することが出来る（図 7）。単体 GPU ではこのようなアプロー
チは有効であるが、並列計算を行う場合には、既に述べた通信時間の隠蔽と重複するため、大量の出力に対しては有効な手法とは言えない。

図 6 並列数 - FLOPS の関係（GPU：GPU 版プログラム、CPU：CPU 版プログラム）

3.5

ファイル出力機能の追加と高速化

GPU の演算速度は極めて高速であるため、各

図 7 ファイル出力の隠蔽方法

タイムステップの Open/Close 処理などオーバーヘッド処理が差分計算ソルバ全体のボトルネック

ここでは、GPU での計算結果を CPU に転送後、となることが予想されたため、CPU のメモリを利メモリを利用して出力のバッファリングを行うこ用して出力のバッファリングを行うことにより、とにより、出力時間を低減する実装を行い評価し出力時間を低減する実装を行い評価した。

た。評価方法としては、複数 GPU を使用した際

GMS は計算結果の出力に、ネットワーク透過で自の GMS ソルバの実行時間を測定し、並列度の変己記述的なバイナリファイルを入出力できる HDF5化に伴う性能評価を実施した。約 2.2 千万格子かライブラリ(http://www.hdfgroup.org/HDF5)を利らなるモデル Unit420 を単位モデルとし、これを用している。GMS の差分計算ソルバは、用途に応じ水平方向に 2x2、4x4、8x8 個並べて、全体の格子ていろいろな出力条件を指定できるが、ここでは数がそれぞれ 4 倍(Unit840)、16 倍(Unit1680)、最も自由度の高い 1D Dump 機能（任意の位置の計64 倍(Unit3360)となるようなモデルを使用して、算結果を出力刷る機能）を GPU 版の差分計算ソル出力のパフォーマンスを評価した。出力対象観測

バに実装した。この出力機能では、毎タイムステ点数は全モデルで一定とし、現在実用上最大規模ップに指定した各観測点での波形を HDF5 でファイのモデル計算で設定している規模の 64 万点としルに追記する。

た。

モデル規模が大きくなると出力のすべきデータも

ま ず 、 デ ィ ス ク I/O が 比 較 的 低 速 な PC

膨大になるが、一般的に出力の並列性能を高める(CPU:Intel Xeon E5520, GPU:NVIDIA Teslaことは困難である。また、GPU での計算処理が極S2050, HDD:SATA300 1TB 7200rpm)を用い、めて高速であるため、相対的にファイル出力に要Unit420 で評価を実施した。図 8 に書き出しに関5

学際大規模情報基盤共同利用・共同研究拠点 平成 23 年度共同研究 最終報告書 2012 年 5 月

4

わるデバイス-ホスト間（D2H）、ホスト-ホスト間

今後の展望（H2H）、ホスト-ファイル間(H2F)、それぞれのハザード評価を行う場合、面的に計算結果を出転送時間を個別に計測して示す。バッファリング

力刷る必要がある。PC 環境のように出力環境が比をしない場合（ステップ数 1）に比べ、バッファ較的低速な場合にはバッファリング等の工夫によ

リングステップ数を増やしていくと、ホスト-ファり出力が高速化されることが分かったが、イル間の転送時間が減少し、バッファリングステ

TSUBAME2.0 のように高速な出力環境を備える

ップ数 100 では約 1/3.5 となった。PC 環境では、システムにおいては、このようなアプローチでは１回の書き込みデータ量を増やすバッファリング不十分であることが分かった。大量の面的出力をを用いた実装が効果的であることがわかった。

得るには、現在のパフォーマンスでは計算時間よりも出力時間の方が長くかかることから、根本的450.0
400.0に出力容量を減らすことが、計算全体のターンア消費時間(s)350.0

ラウンドタイムを短縮には必要であることが分か

300.0
250.0
200.0

H2F

150.0

H2H

100.0

D2H50.0

その他

った。ポスト処理の 1 つであるフィルタ処理を行うためには十分細かい時間間隔での出力が必要であるが、ソルバにおける計算と同時にフィルタ処0.0110100

理を行うことで、出力量を劇的に削減出来る可能

波形出力点数 640000 バッ ファリングステップ数性がある。そのためには、GPU 側にデシメーション処理を逐次フィルタとして実装する必要がある。

図 8 波形出力点数 640,000 の際のバッファリ

今後、ソルバ内外のファイル I/O を最適化するこ

ング効果 (PC ローカルディスク)とにより、シミュレーションだけでなくフィルタ

次に TSUBAME2.0 において、Unit420 を用い処理、可視化処理等のポスト処理を含む全ターンて評価した結果を図 9 に示す。TSUBAME2.0 のアラウンドタイムを短縮し、より実用的なシミュ場合、バッファリングステップ数を増やすとホスレーション環境を構築し、実際の地球をモデル化

ト－ファイル間の転送時間はむしろ増大することした非常に大きな不均質性をもつ地下構造モデルがわかった。これは、TSUBAME のファイルシスや面的広がりをもつ断層モデルを用いた現実的でテムが極めて高スループットで、かつ Open/Close大規模モデルによるシミュレーションを目指す。

処理が高速であるため、I/O よりもバッファリン
グ処理に係るオーバーヘッドが比較的高コストに

5

研究成果リストなっていることが原因と考えられる。

(1) 学術論文(2) 国際会議プロシーディングス
S. Aoi, N. Nishizawa, T. Aoki (2012) Large
scale simulation of seismic wave propagation
using

GPGPU,

15thWorld

ConferenceonEarthquake Engineering, in Lisbon「投稿中」.
(3) 国際会議発表(4) 国内会議発表
前田宜浩・森川信之・青井真・藤原広行，2011，
南海トラフの巨大地震による長周期地震動に関する検討，日本地震学会秋季大会，P2-51．

図 9 ファイル出力時間の評価結果(5) その他（特許，プレス発表，著書等）
6