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

jh130031-NA18

階層分割型数値計算フレームワークを用いた波源から地上構造物までの実地形津波解析
室谷浩平（東京大学大学院 工学系研究科 システム創成学専攻）
概要 本研究では，津波が地上構造物へ及ぼす影響を解析するためのマルチフィジックス解析を 3 つに分けて行う．第 1 の解析は，震源で発生する波源から沿岸部までの津波伝播計算である．第 2 の解析は，沿岸部に押し寄せた津波が地上へ遡上し浮遊物を運搬する計算である．第 3 の解析は，構造物が受ける水圧と浮遊物の衝突による構造解析である．津波が地上構造物へ及ぼす影響を解析するために，それぞれの解析ステージで適した解析手法を用い，それぞれの解析ステージにおいて高性能で高精度な解析を実現する．

1.

研究の目的と意義析手法を用いる．これらの異なる 3 つの解析手法

2011 年 3 月 11 日の東日本大震災によって引き起を効率的に組み合わせるために，我々は「階層分こされた津波により，東日本は甚大な被害にみま割型数値計算フレームワーク」を開発し，そのフわれた．構造物への被害は，水圧によるだけでなレームワークを用いて，それぞれの解析ステージ

く，浮遊する瓦礫が構造物にぶつかり構造物に被において高性能で高精度な解析を実現する．

害を及ぼすケースが多くあった．我々は，津波による構造物への被害を検討する際には，津波の遡2. 当拠点公募型共同研究として実施した意義

上計算を行うだけでは不十分と考えており，浮遊(1) 共同研究を実施した拠点名および役割分担

物の衝突を含んだ構造解析が必要であると考えて「階層分割型数値計算フレームワーク」の開発にいる．

関係する役割分担は次の通りである．


本研究では，津波が地上構造物へ及ぼす影響を荻野正雄(名古屋大学情報基盤センター)，塩解析するためのマルチフィジックス解析を 3 つに谷隆二(東洋大学総合情報学部総合情報学

分けて行う．第 1 の解析は，震源で発生する波源科)：階層領域分割型ソルバーライブラリの開発から沿岸部までの津波伝播計算である．この解析（大規模並列計算向け連立 1 次方程式反復解は，数十～数百キロ四方の広範囲な解析が必要と

法ライブラリ）なるため，粘性浅水波方程式を有限要素法で解く

塩谷隆二，河合浩志(諏訪東京理科大学シスことにする．第 2 の解析は，沿岸部に押し寄せた

テム工学部)，石川格(東洋大学計算力学研究

津波が地上へ遡上し浮遊物を運搬する計算である．

センター)：連続体力学向け DSL（Domain

この解析は，遡上計算と浮遊物の取り扱いが容易

Specific Language，問題領域専用言語）の開な粒子法を用いて計算を行う．粒子法は，解析領発（連続体物理モデル向け行列・テンソル積ラ域に比例する粒子数を必要とするため，数 km 四方

イブラリとトランスレータ）


程度の中規模な範囲の解析しか行うことができな

和田義孝(近畿大学 理工学部 機械工学科)，

い．第 3 の解析は，構造物が受ける水圧と浮遊物

荻野正雄：VSCG ライブラリ（Versatile Scientificの衝突による構造解析である．この解析は，複数Computing Graphics library）の開発（大規模並材料の弾塑性非線形解析であるため，有限要素法列オフラインレンダリングライブラリ）を用いて計算を行う．
我々は，津波が地上構造物へ及ぼす影響を解析「各解析の物理モデル」の開発に関する役割分担するために，それぞれの解析ステージで適した解は次の通りである．
1

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



金山寛(名古屋大学情報基盤センター)，檀啓

きるようになる．これらを，複数のアーキテクチャ上で(九州大学大学院工学研究院機械工学部門)：最適化し，ポストペタコンに向けたアプリ側の課題を津波を粘性浅水波方程式で解くための物理モ明確にすることで，計算科学への貢献を行うことがでデルの構築きる．また，研究成果を幅広く還元することができるよ越塚誠一(東京大学大学院工学系研究科システうになる．

ム創成学専攻)，室谷浩平(東京大学大学院工本研究課題では，荻野正雄，塩谷隆二，河合浩志，学系研究科システム創成学専攻)：津波の沿岸

石川格，和田義孝らが開発した「階層分割型数値計から遡上までを浮遊物を含めて粒子法で解くた算フレームワーク」を利用して，金山寛と檀啓は「粘性めの物理モデルの構築浅水波方程式による津波解析」，越塚誠一と室谷浩平

塩谷隆二，河合浩志：地上構造物への水圧や浮が「MPS 法による津波遡上解析」を行っている．この様遊物の衝突による構造解析のための物理モデに，各種計算機に応じたチューニングと各解析の物理ルの構築モデルの構築を分業することで，開発コストを大幅に「連続体力学向け DSL」の開発と「各解析の物理モ向上させることに成功した．今回は，高並列かつ長時デル」の開発は連携しながら開発を進めている．

間利用における課題を明らかにするため，九大CX400(Xeon + K20)と名大 CX400(Xeon + Xeon Phi)

(2) 共同研究分野を対象にして研究開発を進めている．

大規模数値計算，高精度高速計算，大規模可視化，粘性浅水波方程式による津波解析，MPS(Moving

3.

Particle Simulation)法による津波遡上解析，地上構(1) 研究成果の詳細について

造物の衝突構造解析．

【中間報告書までの成果】

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

本研究課題では，3 段階の津波計算を行っている．
(3) 当公募型共同研究ならではという事項など

第 1 の解析では震源で発生する波源から沿岸部まで連続体力学系シミュレーションは，従来からの有限の津波伝播計算(数十～数百キロ四方程度)，第 2 の要素法，有限体積法，境界要素法など非構造格子ア

解析では沿岸部に押し寄せた津波が地上へ遡上す

プローチに加え，近年ではメッシュフリーや粒子ベーる解析(数 km～10km 四方程度)，第 3 の解析では，ス手法による大規模計算の台頭が目覚しい．非構造

市街地に浮遊物が衝突しながら浮遊する市街地浸水格子および節点，粒子ベースモデルを統一的に扱え解析(500m 四方程度)である．現在，第 1 解析と第 2る連続体力学系シミュレーション向けフレームワークの解析が完了したので，その研究成果を述べる．

が整備されれば，連続体力学系シミュレーションでの本研究課題では，石巻市街地の津波解析をターゲスーパーコンピュータ利用の拡大に大きく貢献するこットとしている．図 1 は東アジアでの石巻市の位置，図とができる．さらに，「連続体力学向け DSL」を用いれ2 は宮城県での石巻市の位置を示している．

ば，数学記述（例えば，行列，テンソルの式を tex 形式で記述）から各アクセラレータタイプ向けに最適化されたコードを生成するので，ユーザーが連続体物理モデルを自由にカスタマイズすることができる．この「連続体力学向け DSL」の機能によって，スパコンの知識が無い人でも，スパコンを利用して，新しい連続体物理モデルを解くことができるようになるため，従来と比較にならないくらい高精度で大規模な連続体力

図 1．東アジアでの石巻市の位置学系シミュレーションを，短い開発期間で行うことがで（©2013 Google）
2

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

第 1 の解析は，震源で発生する波源から沿岸部までの津波伝播計算である．解析領域は，667km×525km で，要素数（三角形数）17,970,848 の約 250m のメッシュを用いてモデルを作成した．海底座標は，JTOPO30 のグリッド水深データ（1km の解像度）を用い，時間ステップ 0.25 秒で 4,000 秒間の計算をした．初期値には，東日本大震災時の津波の初期水位を用いた．
計算モデルのサイズから，第 1 の解析は MPI 並列を
図 2．宮城県での石巻市の位置行わず，OpenMP 並列のみ行うことにした．本解析を（©2013 Google）

Intel core i7 960 (3.2GHz) 4core を搭載した PC で解析した結果，7.7 時間かかった．図 5 は，1,000 秒後の

図 3 は石巻市から提供して貰った 5m 間隔の標高デ水位である．

ータを元に作成した石巻市の地形モデルである．テクスチャーには，NTT 空間情報から購入した航空写真を用いている．

図 3．石巻市の地形モデル（©NTT 空間情報（航空写真））
図 4 は石巻市の市街地モデルである．Google が提供する建物 3D から，東日本大震災前の町並みを再図 5．1,000 秒後の水位現したポリゴンモデルである．テクスチャーには，NTT第 2 段階の津波遡上解析の入力条件に，第１段階空間情報から購入した航空写真を用いている．の粘性浅水波方程式による解析結果を用いた．第 2段階の津波遡上解析は，解析領域 4.5km×3km のモデル（解析モデル a）と，解析領域 10.5km×10km（解析モデル b）のモデルを作成し，それぞれ，800 秒間と2000 秒間の解析を行った．
図 6 では，粒子直径 1m の石巻市に津波が遡上する解析を行った結果である．実時間 800 秒の解析に，
図 4．石巻市の市街地モデル

東京大学情報基盤研究開発センターの FX10 の 144（©NTT 空間情報（航空写真））

計算ノードを用いて約 7 日間かかり，最大 2.6 億粒子の解析となった．
3

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

州大学情報基盤研究開発センターの CX400 の 32 計
算ノードを用いて約 1 ヶ月かかり，最大 4 億粒子の解析となった．本解析では，固定されていた水色とピンク色の 2 個のタンクタンクが，60 秒後に解放され流されていく解析である．

図 6．解析モデル a（400 秒後）
図 7 では，粒子直径 2m の石巻市に津波が遡上する解析を行った結果である．実時間 2000 秒の解析に，筑波大学計算科学研究センターの T2K-Tsukuba の32 計算ノードを用いて約 13 日間かかり，最大 1.3 億

図 8．解析モデル c（90 秒後）粒子の解析となった．
図 9 では，粒子直径 0.3m の石巻市街地に津波が侵入する解析を行った．実時間 400 秒の解析に，九州大学情報基盤研究開発センターの CX400 の 32 計算ノードを用いて約 1 ヶ月かかり，最大 4 億粒子の解析となった．本解析では，固定されていた 431 個の地上構造物が，120 秒後に解放され流されていく解析である．図 9 では 431 個の地上構造物が上に向って流されている様子が見て取れる．

図 7．解析モデル b（800 秒後）
第 3 段階の解析は，遡上した津波が浮遊物を運搬し，構造物に衝突する解析である．第 3 段階の解析で
図 9．解析モデル d（240 秒後）は，2 種類の解析モデルを作成した．一つ目の解析では，400m×550m の市街地をモデル化して，2 個の浮遊物(10m 級タンク)を流す市街地浸水解析を 200

第 2,3 段階の解析モデルは次の通りである．

秒間行った（解析モデル c）．二つ目の解析では， （解析モデル a）第 2 段階の 800 秒間の津波遡上解析(4.5km×3km)

660m×810m の市街地をモデル化して，431 個の浮遊物(地上構造物)を流す市街地浸水解析を 400 秒間行 （解析モデル b）第 2 段階の 2,000 秒間の津波遡上解った．現段階では，浮遊物は剛体であり，構造解析は析(10.5km×10km)

行っていない．

（解析モデル c）第 3 段階の 2 個の浮遊物(10m 級タン

図 8 では，粒子直径 0.2m の石巻市街地に津波がク)ある 200 秒間の市街地浸水解析(400m×550m)

侵入する解析を行った．実時間 200 秒の解析に，九 （解析モデル d）第 3 段階の 432 個の浮遊物(地上構造4学際大規模情報基盤共同利用・共同研究拠点 平成 25 年度共同研究 最終報告書 2014 年 5 月物)ある 400 秒間の市街地浸水解析(660m×810m)【中間報告書後の成果】

各解析の計算時間は，表 1 にまとめる．

図 10 では，解析モデル c を用いて，地上構造物を本研究課題で，第 2,3 段階の解析で最も計算量を線形弾性体として，圧力による応力解析を行った．流必要とする，並列 MPS 法コードのピーク性能値を示

体・剛体粒子から地上構造物への力の渡し方は，流す．表 2 は MPS 法コードのピーク性能値である．九大体・剛体粒子の圧力勾配から求めた面荷重を外力とCX400 はプロファイラーが無いので東大 FX10 で測定して与える片方向連成方法を用いた．構造解析には，した推定値である．SPARC 系と Intel 系向けにチュー

ADVENTURE_Solid を用いた．本解析では，名古屋ニングされた連続体向け DSL を用いてチューニング

大学 CX400 を用いた．図 10-a,b,c のように時間が経した結果，全体で東大 FX10 では７％，九大 CX400過する（浸水領域が広がるに）につれて，圧力によっでは 12％でているので，連続体計算としては十分なて応力分布が発生している地上構造物が広がってい結果であると言える．

ることが見て取れる．地上構造物 1 つに注目した所，水流の変化が激しい，角の部分に大きな応力が発生

表 1．各解析モデルの計算時間
解析
対象abccd

計算機
東大 FX10－144 ノード
筑波大 T2K－32 ノード
東大 FX10－600 ノード
九大 CX400－32 ノード
九大 CX400－
32 ノードしていることが見て取れる．

粒子直粒子

計算時

径数

間1m

2.6 億

7 日間2m

1.3 億

13 日

0.5m

4000
万6 時間
図 10-a．地上構造物の応力解析（30 秒後）

0.2m

4億

1 ヶ月

0.3m

4億

1 ヶ月

表 2．MPS 法コードのピーク性能値
東大 FX10

九大 CX400

Flops/PeakFlops/Peak7.426212.8124

11.139619.0215

Collision (例外処理)

7.422412.5120

圧力計算

7.2425

12.6185

7.055512.1003

Total
重力･粘性による力の計算 (空間 2 次精度)

圧力勾配による力の計
算 (空間 1 次精度)

図 10-b．地上構造物の応力解析（75 秒後）

図 10-c．地上構造物の応力解析（100 秒後）
5学際大規模情報基盤共同利用・共同研究拠点 平成 25 年度共同研究 最終報告書 2014 年 5 月

非ダイレクトメモリアクセスから，粒子データの図 11 は，JST-CREST プロジェクト研究課題「ポソートによる直接参照のダイレクトメモリアクセ

ストペタスケールシミュレーションのための階層スにアルゴリズムが変更されている．

分割型数値解法ライブラリ開発」で開発されたマIntel 系 CPU,MIC のコンパイラオプションは-fast

ルチプラットホーム大規模領域分割型並列可視化である．CUDA 向けチューニングでは，CUDA のライブラリ VSCG（Versatile Scientiﬁc Computing

一般的な最適化をおこなってある（constant メモGraphics library）を用いて生成した 15,000×20,000

リやキャッシュの活用など）．SPARC64 向けチュピクセルの画像である．VSCG がこれまでの非構ーニングでは，novrec，swp，simd の 3 つの最適化

造格子向けに持っていた機能を，粒子系シミュレ指示子を最内ループに設定して，コンパイラオプータに対応できるよう拡張した．これまでは，解ションに-Kfast,openmp,parallel –Kocl を設定して，析モデル c で生成された 8TB のバイナリデータをソフトウェアパイプラインと SIMD ベクトル化が九州・東京間でインターネットを用いたデータ転有効にしてある．今回は，CPU 向けと GPU 向け

送に 1 ヶ月かかっていた．VSCG ライブラリを用

チューニングは進んでいるが，MIC 向けチューニ

いれば，10 万×10 万ピクセルを超える高解像度画ングは，まだ行っていない前提での測定結果の比像をスパコン上で作成することができる．また，較である．

VSCG は依存ライブラリが数学関数のみであるこ

空間微分離散化モデルには，標準 MPS の簡易微とから，将来を含めた多様な計算機環境で動作す分モデルと，空間高次精度モデル(ラプラシアン 2ることが可能である．

次多項式精度，勾配 1 次多項式精度)を用いた．空間高次精度モデルは，近傍粒子の相対位置ベクトルの外積の重ね合わせをとるモーメントマトリッ
クスの逆行列を求めるので，演算量が非常に多いのが特徴である．
表 3 から，標準 MPS モデルと高次精度モデルの間で，CPU と MIC では似た傾向あることが分かる．
一方で，GPU は，標準 MPS モデルでは非常に高速であるにも関わらず，高次精度モデルでは，非常に低速になっている．PTX でスレッド当りのレジスタ数を確認すると，標準 MPS モデルでは 50以下である一方で，高次精度モデルでは 200 以上もレジスタを必要とする．このため，高次精度モ

図 11．VSCG による可視化結果

デルを GPU で解くと，レジスタ数が不足して，グローバルメモリへスワップが発生し演算性能が低

表 3 は，影響半径＝(平均粒子間距離×4.1)として， 下していると推測できる．CPU や MIC では，標準132,651 粒子のダムブレイク解析の 1 タイムステッ

MPS モデルと高次精度モデルの間で演算量の増加

プ当りの計算時間を，Intel Core i7-3820，Intel Xeon

ほどの計算時間に差は発生しておらず，メモリバ

E5-2697 v2 x2(名古屋大学 CX400)，SPARC64TM

ンド幅ネックの演算になっていることが分かる．

Ixfx(東京大学 FX10)，NIDIA Tesla K20m(九州大学

表 3 の結果から，標準 MPS モデルを用いる場合は，CX400)，Intel Xeon Phi 3120P(名古屋大学 CX400)

GPU を用いると最も高速に解けるが，高次精度モを用いて比較した結果である．

デルを用いたいときは，Xeon や MIC を用いる方

いずれのコードも，リンクリストの間接参照のが高速に解けることが分かる．
6

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

表 3．MPS 法コードのプロセッサ間の比較

プロセッサ

理論性能

標準 MPS モデル

高次精度モデル(GFLOPS)

Time(s)Speed upTime(s)Speed upCore i7-3820

115.2

1.358

1.000

1.839

1.000

Xeon E5-2697 v2 x2

518.4

0.310

4.3780.414

4.446

Tesla K20m

1174.7

0.158

8.612

1.333

1.379

Xeon Phi 3120P

1003.2

0.444

3.059

0.614

2.995

5.

(2) 当初計画の達成状況について

研究成果リスト

今年度の開発進捗は 90％とする．課題申請書に(1) 学術論文（投稿中のものは「投稿中」と明記）記述した研究計画で達成できなかった項目は下記[1] 荻野正雄, “領域分割法の対角スケーリング前の 2 項目である．

処理に関する考察”, Transaction of JSCES, Vol.

[1] 【10 月～3 月】3. 粒子法を「階層分割型数値

2013,20130013,2013(DOI:10.11421/jsces.2013.20130013).

計算フレームワーク」を用いて MPS 半陰解法

[2] Hiroshi Kanayama, Masao Ogino, Shin-ichiroを開発し，MPS 陽解法とのトレードオフを検討する．
【未達成箇所：達成率 50％】MPS 半Sugimoto

and

陰解法を開発は行ったが，比較検討までにはMagnetostatic Domain Decomposition Analysis

至らなかった．

BasedontheSeigo

Terada,MINRES

Large-ScaleMethod,

IEEE

Transactions on Magnetics, vol. 49, No. 5,[2] 【10 月～3 月】6. HDDM 構造解析：1000 万要
素のモデルで物理時間 15 分間の解析を行う．

pp.1565-1568,

May2013.

【未達成箇所：達成率 90％】構造解析は 100

10.1109/TMAG.2013.2238612)(DOI:万要素程度の解析に留まった．この項目は主[3] Kohei Murotani, Shin-ichiro Sugimoto, Hiroshiに連成解析のためのインターフェイス部分のKawai and Shinobu Yoshimura, Hierarchical開発が目的であった．更に，今回は粒子サイ

Domain Decomposition with Parallel Meshズとの関係で 100 万要素程度でも十分な要素Refinement for Billions-of-DOF Scale Finite数となったため，達成率を 90％とした．

ElementAnalyses, International Journal ofComputational Methods, published via internet,

4.

21

今後の展望

June2013.

(DOI:

10.1142/S0219876213500618).

今年度は，九大 CX400 の TeslaK20 と名大 CX400に Xeon Phi の利用のおかげで，アクセラレータ向け[4] Hiroshi KANAYAMA and Hiroshi DAN, A開発が大きく進んだ．今後も引き続き，名大 CX400 を

Tsunami Simulation of Hakata Bay Using the用いて，Xeon Phi 向けチューニングを進めて行っていViscous Shallow-water Equations, Japan Journal

く予定である．

of Industrial and Applied Mathematics, Volume

本研究課題では，地上構造物は弾性体として扱っ

30, Issue 3, pp 605-624, 2013/11. (DOI:たが，浮遊物は剛体のままであった．今後は，浮遊物

10.1007/s13160-013-0111-7)も弾性体として解析を行えるようにするのが次の目的[5] A.M.M. Mukaddes, Masao Ogino, and Ryuji

である．更には，浮遊物と地上構造物の弾性衝突を

Shioya, “Performance Evaluation of Domain

解けるようにするのが本研究の最終目的の一つとな

Decomposition Method with Sparse Storage

る．

Schemes in Modern Supercomputer”, Special7

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

IJCM Issue based on ICCM2012 (in press).
[6] N.

MITSUME,S.

decomposition explicit MPS method (InvitedYOSHIMURA,K.

Lecture)”,5thAsia

PacificCongresson

MUROTANI and T. YAMADA: MPS-FEM

Computational Mechanics & 4th International

PARTITIONED COUPLING APPROACH FOR

Symposium

FLUID-STRUCTURE INTERACTION WITH

(APCOM & ISCM 2013), Singapore, December

FREE SURFACE FLOW, International Journal of11-14, 2013.

Computational Methods, (in press).

[2] Masao

on

Ogino,

Computational“AMechanics

Balancing

Domain[7] Kohei Murotani, Seiichi Koshizuka, Tasuku

Decomposition Method Combined with Diagonal

Tamai, Kazuya Shibata, Naoto Mitsume, Shinobu

Scaling Preconditioner for Multi-materials (KeyYoshimura, Satoshi Tanaka, Kyoko Hasegawa,note)”,

Eiichi

Computational Mechanics & 4th InternationalNagaiandToshimitsu

Fujisawa:

5thAsia

Congress

Symposium

Objects(APCOM & ISCM 2013), Singapore, DecemberaHierarchical

DomainComputationalon

Large-scale Tsunami Analysis with FloatingUsingon

Pacific

Mechanics

Decomposition Explicit MPS Method, Journal of

11-14, 2013.

Advanced Simulation in Science and Engineering(4) 国内会議発表(JASSE), 「投稿中」.

[1] 室谷浩平，越塚誠一，金山寛，安重晃，藤澤智[8] KentaroTanaka,

SatoshiTanaka,Kyoko

光，玉井佑, 柴田和也, 三目直登，吉村忍：津Hasegawa, Kohei Murotani, Seiichi Koshizuka:

波波源から浅水波方程式を用いて生成した入力

Translucent Visual Analysis of Large Scale 3D

境界条件による MPS 法市街地遡上解析，日本

Point Data Generated by Particle Fluid, Journal of

機械学会 第 26 回計算力学講演会（CMD2013），

Advanced Simulation in Science and Engineering

佐賀大学，2013 年 11 月 2 日--4 日．

(JASSE), 「投稿中」.

[2] 室谷浩平，越塚誠一，玉井佑，山田祥徳，柴田(2) 国際会議プロシーディングス

和也，三目直登，吉村忍，田中覚，長谷川 恭子，藤澤 智光：粒子系シミュレーション並列化(3) 国際会議発表[1] Kohei

フレームワークの開発と大規模津波解析，第 18

MUROTANI,Seiichi

KOSHIZUKA,

回計算工学講演会，東京大学生産技術研究所,Kazuya SHIBATA, Tasuku TAMAI, Naoto

2013 年 6 月 19 日--21 日.

MITSUME, Shinobu YOSHIMURA, Satoshi
TANAKA, Kyoko HASEGAWA, and Toshimitsu(5) その他（特許，プレス発表，著書等）

FUJISAWA,realistic[1] 越塚誠一，柴田和也，室谷浩平：粒子法入門，

tsunami analysis

with floating objects on

丸善出版，256 ページ，2014 年 6 月下旬初版

Ishinomaki

using出版．

“Large-scalecityparallelhierarchicaldomain8