（様式 20）成果報告書

TSUBAME 共同利用 令和元年度

学術利用成果報告書

利用課題名 LRnLA アルゴリズムを用いた物理シミュレーション
英文：Simulation of Physical Processes with LRnLA Algorithms

善甫 康成
Yasunari Zempo
法政大学 情報科学部
Computer and Information Sciences
http://cis.k.hosei.ac.jp/
地震解析のための DiamondSeism コードは、毎秒最大 14 億セル更新のパフォーマンスで、層状媒体の地震弾性波をシミュレートできる性能を持っている。 この性能は LRnLA 並列アルゴリズムを使用使うことで達成されたものである。 昨年報告したシミュレーションを継続し、更に多くの地震情報を取得している。 表面波、屈折波、二次波など総てのタイプの弾性波を考慮に入れた 3D の有限差分モデルにより、地震応答の正確な予測ができるようになった。3D
クロスステンシルを持つ新しいモデルを用いて LRnLA コードの最適化のテストを実行した。
この新しいアルゴリズムは、今後様々な物理シミュレーションコードの開発で使用できると期待できる。
The code DiamondSeism simulates seismic elastic waves in a layered media with aperformance up to 1.4 billion cell updates per second. This performance is achieved with theuse of parallel LRnLA algorithms. We continue the simulation batch that has been reported
in the previous year, and obtain new seismograms in the extended domain. Taking accountof all types of elastic waves: surface wave, refracted waves, secondary waves, finite difference
modelling in 3D media gives accurate predictions of the seismic response. We have
performed tests for optimization of LRnLA codes on a model with a 3D cross stencil. The newalgorithm will be used in further physical simulation code development.

Keywords: LRnLA Algorithm, Hooke’s law, synthetic seismogram, finite difference, crossstencil

Depth (km)4 次近似の数値スキームを使用すれば、空間ステッ
背景と目的
プは一番短い波長の 8 分の 1 程度で十分であり、必永久凍土(permafrost)がある地域での地震波ト要な精度が得られるモデルリングが可能である。不モグラフィーを用いる際に注意しなければならないことの一つに不凍土層(Talik)の存在の有無がある。 凍土層の P 波の伝播速度は𝑉𝑉𝑝𝑝 ~2400 m/s, S 波 の伝不凍土層があると地震計の記録に干渉によるノイズ
播速度は𝑉𝑉𝑠𝑠 ~1550 m/s であるので、長さが~10 kmが発生することが知られている。3 次元全波数値シ
程 度 の シ ミ ュ レー シ ョ ンを 行 う と 、 デ ータ 量 はミュレーションを使用すると、推測した位置にある50GB 以上にる。また幾つかの震源位置での地震解
不凍土層を含む地震層モデル
Layerから理論的な地震波が得られ0.0
Talik1る。その記録を解析すること
2-1.0
3により不凍土層に起因した干4-2.0
渉の検出が可能になる。更に
5地震波トモグラフィーの解析-3.06について補正を行い、地震波-4.0を正確にとらえることが可能-5.0である。
𝑥𝑥020406080
100120140160
180 km
計算コストの観点から、不𝑧
凍土層は低速度層であるの
図 1．地震層モデル。 Layer1 の 141200m-144925m が Talik 層でで、解析の解像度は空間ステある。地震源は 134925m から 151200m の空間にあるとした。
ップによって決まり、それがモデルの全計算量を 決める。

（様式 20）成果報告書

𝜌 (g/cm−3)020406080100120140160180

図 2 密度の空間依存性

4.23.8km(a)3.4𝑉𝑉𝑝𝑝 ( × 103 m/s)

内容
地震波の伝播の数値モデルはフックの法則とニュートン運動の法則からなる[2]。モデルは、応力テンソルと変位速度ベクトルの時間発展をあらわす。
使用されている Finite-Difference 数値スキームは空間で 4 次近似、時間で二次近似である[3]。境界条件として PML を使用し、地震波源のため TF/SF 法を使用している。なお本モデルは以前、DTmaxwell コード[4]を基本に TSUBAME2.5 とTSUBAME3.0 上で開発してきた解析を実施してきた実績のあるものである。
図 1 に地震層モデルを示す。地表から地球内部に向かって主要な構造が 6 層ある.これらの層を用
い図 2 に示す層のパラメータ依存性を考慮しシミュレーションを行った。密度の空間依存性、P 波および S 波の伝播速度との空間依存性を、それぞれ
図 2、図 3 に示す。本プロジェクトで用いたモデルは 3 次元データを用いた 3 次元解析であるが、1断面をとり 2 次元の表示としている。伝播速度が低い Talik 層は Layer1(赤線)の中の 141200 m −144925 mに示してある。
本コードの特徴は、LRnLA アルゴリズムの使用していることである。この数値計算スキームでは局所依存性を考慮し、メモリバンド幅への負担をできるだけ最小化・緩和する。また LRnLA アルゴリズ
ムは、全メモリ階層と全並列法を考慮し temporal
blocking[5] の 考 え 方 を 更 に 進 め た 形 に な って いる。 GPU の場合、計算のデータを比較的に大きいレジスタファイル内に閉込め CUDA-スレッドの並列性を使う[4]。

2.62.42.22.01.8
1.61.41.21.03.0

図 4. Roofline 解析による性能評価

2.62.2020406080100120140160180km(b)2.42.2𝑉𝑉𝑠𝑠 ( × 103 m/s)

析の結果を得るには、必要 な計算を 1000 ケース以上行う必要がある。
このためには効率的なアルゴリズムを持つ並列コードと高性能計算機の利用が不可欠となる。本プロ
ジェクトでは、GPU を使用して地震波の伝播現象をシミュレーションする高性能コード
DiamondSeism を many-GPU の 環 境 で あ る
TSUBAME3.0 に移植し不凍土層を含む地震層モデル計算を実施した[1]。
加えて並列効率とメモリ効率を更に向上させるため、新しい LRnLA アルゴリズムと実装法を開発し、その試験を実施した。

2.01.81.6
1.4020406080100120140160180図 3 P 波(𝑎𝑎)および S 波(𝑎𝑎)の伝播速度との空間依存性km

本プロジェクトでは、スレッドの並列に加え、新しい試みとして CUDA-warp 並列と shared メモリアクセス数を減じて実装効率増加させるという実装法を導入した。
図５に RoofLine[6] 解析による性能評価の結果を示す。３次元クロスステンシルを用いたコードのパフォーマンスの見積と NVidia P100GPU での実施結果である。黒い矢印は LRnLA アルゴリズムを用いない従来型コードの最大パフォーマンスの見積である。これは全データ(23.3GB)をデバイスメモリ(HBM2)に保存し、1 回全データのアプデートのため、数値毎にロードして、再度 1 回セーブすることが必要であるということから算出した理想的な場合の見積結果である。
LRnLA アルゴリズムを用いた場合の上限を赤い矢印で示す。従来型のコードとは様子が大きく異なる。その理由は、デバイスメモリからのロードは一番時間が掛かるが、一旦ロードすると、既にロードされているデータを用い多数の計算が行われるため、operational intensity が非常に高く保たれるからである。単純に LRnLA アルゴリズムを用いた結果を図 4 中に桃色の矢印で示す(一番右)。3 次元のデータを DiamondTorre アルゴリズムにて計算した
時 の デ ー タ 量 から 算 出 した 推 測 値 で あ る。 ま た
DiamondTorre 間のデータ交換のため、LLC cacheを使用可能と仮定した結果を、二番目の’Block’という紫色の矢印で示している。また CUDA-blockの中の計算では CUDA-スレッドの並列化を使用するので、更に DiamondTorre 間のデータ交換のた（様式 20）成果報告書
め、LLC cache を使用可能と仮定した結果が、右から二番目の’Block’という紫色の線で示した結果である．かなり処理能力が向上していることがわかる。
また CUDA-block の中の計算では CUDA-スレッドの並列化を使用する。Warp-間のデータ交換には
shared メモリを使用する。これを利用した場合の上限値を見積ったのが緑の矢印で示した値である，もちろん CUDA-warp 中のデータ効果は shuffle ベ
クトル演算を用いるため、更に性能が向上することが期待できる（赤の矢印）。その結果、パフォーマンスの上限値は 490 cell update per second であると見積ることができる。実際のテストでは、オーバーヘッドのためと思われるが、260 cell update persecond であった。

参考文献[1] T. Levchenko, V. Rok, V. Levchenko, A.Perepelkina,
Y. Zempo, “Computer modelling specifics of thegeologicalstructure
withcontrasting
inhomogeneities under the permafrost conditions”
GEO Eurasia-2019, Materials of the International
geologic and geophysics conference and exhibition(Feb. 4-7, 2019, Moscow), pp. 814-817
https://www.gece.moscow/https://drive.google.com/file/d/1Z7o5H8ZQ6EYqgp
atcbBEGlaKiEMFiarh/view?usp=sharing[2] Aki, Keiiti, and Paul G. Richards, “Quantitativeseismology”, 2002.
[3] Levander Alan R. Fourth-order finite-difference P-W
seismograms // Geophysics. 53(11), p1425-1436.,
1988[4] A. Zakirov, V. Levchenko, A. Perepelkina, Y. Zempo,“High performance FDTD algorithm for GPGPU
supercomputers”, J. Phys.: Conf. Ser. 759 012100,
2016[5] Endo, Toshio. "Applying Recursive Tempora
Blocking for Stencil Computations to Deeper
Memory Hierarchy." 2018 IEEE 7th Non-Volatile
Memory Systems and Applications Symposium(NVMSA). IEEE, 2018.
[6] Williams, Samuel, Andrew Waterman, and DavidPatterson.
"Roofline:
aninsightful
visualperformance model for multicore architectures."Communications of the ACM 52.4 (2009): 65-76.

𝑥𝑥 (km)結果および考察TSUBAME3.0 の ノ ー ド 内 の 4 つ の NVIDIA
Tesla P100 を用い計算を行った。
震源の振動数は 30 Hz である。モデルの 134925
m から 151200 m までの空間内で 50 m ごとに 325個の地震源があるとした。地震波センサーは地表に25 m 毎に 650 個設置し、シミュレーション領域は地震計を設置した線に沿って長さ 10 km であり、深 さ は 5 km で あ る 。 空 間 の メ ッ シ ュ 数 は1500×600×256 個である。地震計の線の垂線の方向
で 128 メッシュ内には PML 層も含んでいるが、球形の波面を持つ地震波の伝播距離にからすると減衰には十分であると考えられる。
計算データは 13 GB となるので、TSUBAME のノード内では、 4 つの独立なシミュレーションを実実施することが可能である。時間ステップは~0.667 秒であるため、5 秒の伝播を図るためには15000 ステップ必要である。現在の パ フ ォ ー マン ス は 3.1 × 109 cell update per second であり一つの計算は 12 分程度で済む。
計算結果の例として、図 4 では理論的に求めた地震計の記録が表している。これは、𝑡𝑡 − 𝑥𝑥のグラフで変位の速度の𝑉𝑉𝑧𝑧 である。z-軸は下に向いている。
不凍土層による反射が良く見てとれる。

まとめと今後の課題
我々が、これまでに開発した LRnLA アルゴリズ
ムを使ったコードを地震波の伝播計算に展開した。
また LRnLS に関して GPU のハードに合わせた新しいアルゴリズムの開発も行った。その特徴は、3次元データにおける計算の局所化とレジスター内でのデータ交換である。CUDA を用い実施法の単純化がさらに必要ではありさらに改良を続ける予定である。完成すれば、電磁場、地震波のシミュレーシ
ョンで使用可能になり、現在使用している地震波解析コード DiamondSeism より高いパフォーマンスが期待できる。

𝑡𝑡 (s)

図 6. 理論的に求めた地震波。震源付近での変位の速度の大きさを表している。