第10回 配信講義 計算科学技術特論B (2024)

計算科学技術特論B 大規模地震シミュレーション1 2024/6/20 東京大学 地震研究所 藤田 航平 1

内容・スケジュール • 6/20 1300- • 地震シミュレーションの概要 • 有限要素法の基礎 • 大規模連立方程式の求解方法(共役勾配法) • 共役勾配法における並列化 • 共役勾配法における前処理 • 大規模地震シミュレーションの例 • 6/27 1300- • 有限要素法における大規模線連立程式求解方法の高速化(SIMD・マル チコア等) 4

地震シミュレーション概要 5

地震シミュレーション概要 • 断層～地殻～地盤～構造物 • 領域サイズ: 100 km • 空間分解能: 0.1 m • 従来は構造が比較的単純な地殻中の波動伝播では有限差分法が使われてきた • より複雑構造をもつ都市などを高精度で求解するため、近年では有限要素法 が広まりつつある • 都市規模の問題は数千億自由度規模の超大規模問題になるうえに、有限要素 法ではランダムアクセスが主体となり計算効率が下がる傾向にあるため、高 性能計算技術が必須 Example code: do i=1,nx do j=1,ny a(i,j)=b(i,j)+b(i+1,j)+…. enddo enddo 構造格子(差分法など) Example code: do i=1,n do j=ptr(i)+1,ptr(i+1) a(i)=a(i)+b(index(j)) enddo enddo 非構造格子(有限要素法など) 6

解析例 • ターゲット：断層～地殻～地盤～構造物～社会活動 c) Resident evacuation a) Earthquake wave propagation 0 km -7 km Ikebukuro Ueno Two million agents evacuating to nearest safe site Earthquake Post earthquake Shinjuku Tokyo station Shinbashi 京コンピュータ全系を使った 地震シミュレーション例 Shibuya b) City response simulation T. Ichimura et al., Implicit Nonlinear Wave Simulation with 1.08T DOF and 0.270T Unstructured Finite Elements to Enhance Comprehensive Earthquake Simulation, SC15 Gordon Bell Prize Finalist 7

断層-都市-社会応答シミュレーション • 断層-都市-社会応答からなる 統合地震シミュレーションを 実施 • 想定首都直下地震に対する東 京中心部の応答を計算 10 km 8

断層-都市-社会応答シミュレーションでできるように なること • 物理シミュレーションの積み重ねとして，断層から社会応答までの地 震シミュレーションを実現 • 統計・経験式ベースの被害想定に比べて，科学的合理性が向上 • 耐震補強の優先順位付け・地震保険料の合理的設定など，地震対策の高度化 につながると期待 断層-都市シミュレーション 都市シミュレーション 都市-社会シミュレーション 15

有限要素法の基礎 16

有限要素法（1次元） • ターゲット問題(支配方程式): 𝑓𝑓 𝑢𝑢(𝑥𝑥) = 0 • e.g.: 𝑓𝑓 𝑢𝑢 𝑥𝑥 𝑑𝑑 2 𝑢𝑢(𝑥𝑥) =1− for 0 < 𝑥𝑥 < 1 with boundary conditions 𝑢𝑢 0 𝑑𝑑𝑥𝑥 2 = 0, 𝑢𝑢 1 = 1/2 • 有限要素法では未知関数𝑢𝑢 𝑥𝑥 をオーバーラップの無い区間（要素）で分割する • 𝑢𝑢 𝑥𝑥 = ∑𝑖𝑖 𝑢𝑢𝑖𝑖 𝜑𝜑𝑖𝑖 (𝑥𝑥) • ここで𝑢𝑢𝑖𝑖 は定数(未知) 、𝜑𝜑𝑖𝑖 (𝑥𝑥)は形状関数(既知) • 𝑢𝑢𝑖𝑖 が決まれば𝑢𝑢 𝑥𝑥 が求まる→どうやって𝑢𝑢𝑖𝑖 を求めるか？ 節点 𝑦𝑦 𝑢𝑢 𝑥𝑥 要素 𝑥𝑥 𝑦𝑦 1 0 節点𝑖𝑖 𝜑𝜑𝑖𝑖 (𝑥𝑥) 𝑥𝑥 17

有限要素法（1次元） • 重み付き残差法 • 支配方程式(𝑓𝑓 𝑢𝑢(𝑥𝑥) = 0)をそのまま解かずに、任意の重み関数𝑤𝑤 𝑥𝑥 を 使って、以下の式に変換して解く ∫ 𝑤𝑤(𝑥𝑥)𝑓𝑓 𝑢𝑢(𝑥𝑥) 𝑑𝑑𝑑𝑑 = 0 • 𝑓𝑓 𝑢𝑢 𝑥𝑥 1 ∫0 𝑤𝑤 𝑥𝑥 𝑑𝑑 2 𝑢𝑢(𝑥𝑥) =1− の場合(0 < 𝑥𝑥 < 1)では、 𝑑𝑑𝑥𝑥 2 1 𝑑𝑑 2 𝑢𝑢 𝑥𝑥 𝑑𝑑𝑑𝑑 𝑥𝑥 𝑑𝑑𝑢𝑢 𝑥𝑥 1− 𝑑𝑑𝑑𝑑 = ∫0 𝑤𝑤 𝑥𝑥 + 𝑑𝑑𝑥𝑥 2 𝑑𝑑𝑑𝑑 𝑑𝑑𝑥𝑥 𝑑𝑑𝑑𝑑 − 𝑤𝑤 𝑥𝑥 𝑑𝑑𝑑𝑑 𝑥𝑥 1 =0 𝑑𝑑𝑑𝑑 0 • 微分の階数が減り(この場合は2→1)、結果として、低次の形状関数を 使えるようになる（弱形式） 18

有限要素法（1次元） • さらに、重み付き残差法の重み関数𝑤𝑤に、𝑢𝑢の離散化に使った形状関数を 用いる(i.e., 𝑤𝑤 𝑥𝑥 = 𝜑𝜑𝑖𝑖 (𝑥𝑥))(ガラーキン法) 1 1 𝑑𝑑 2 𝑢𝑢 𝑥𝑥 𝑑𝑑𝑑𝑑 𝑥𝑥 𝑑𝑑𝑑𝑑 𝑥𝑥 • ∫0 𝑤𝑤 𝑥𝑥 1 − 𝑑𝑑𝑑𝑑 = 𝑤𝑤 𝑥𝑥 + ∫0 𝑑𝑑𝑥𝑥 2 𝑑𝑑𝑑𝑑 𝑑𝑑𝑑𝑑 𝑑𝑑𝜑𝜑𝑗𝑗 𝑥𝑥 1 𝑑𝑑𝜑𝜑𝑖𝑖 𝑥𝑥 ∑𝑗𝑗 𝑢𝑢𝑗𝑗 = ∫0 + 𝜑𝜑𝑖𝑖 (𝑥𝑥) 𝑑𝑑𝑑𝑑 𝑑𝑑𝑑𝑑 𝑑𝑑𝑑𝑑 1 𝑑𝑑𝜑𝜑𝑖𝑖 𝑥𝑥 𝑑𝑑𝜑𝜑𝑗𝑗 𝑥𝑥 1 = ∑𝑗𝑗 ∫0 𝑑𝑑𝑑𝑑 𝑢𝑢𝑗𝑗 − ∫0 −𝜑𝜑𝑖𝑖 𝑥𝑥 𝑑𝑑𝑑𝑑 𝑑𝑑𝑑𝑑 𝑑𝑑𝑑𝑑 𝑓𝑓𝑖𝑖 𝑎𝑎𝑖𝑖𝑖𝑖 = ∑𝑗𝑗 𝑎𝑎𝑖𝑖𝑖𝑖 𝑢𝑢𝑗𝑗 − 𝑓𝑓𝑖𝑖 = 0 − 𝑤𝑤 𝑥𝑥 𝑑𝑑𝑑𝑑 𝑥𝑥 1 𝑑𝑑𝑑𝑑 𝑑𝑑𝑑𝑑 0 境界でのみ値を持つ→ 今回は𝑢𝑢が境界で固定(ディリ クレ境界条件)のため省略 • マトリクス形で書くと 𝐀𝐀𝐀𝐀 = 𝐟𝐟 • 有限要素法の形状関数を使うことで 𝐀𝐀は疎行列となる 19

有限要素法（1次元）例題 1 𝑑𝑑𝜑𝜑𝑖𝑖 𝑥𝑥 𝑑𝑑𝜑𝜑𝑗𝑗 𝑥𝑥 • 𝑎𝑎𝑖𝑖𝑖𝑖 = ∫0 𝑑𝑑𝑑𝑑 𝑑𝑑𝑑𝑑 1 𝑢𝑢1 𝑦𝑦 𝑢𝑢3 𝑢𝑢4 𝑢𝑢 𝑥𝑥 ③ ④: 節点番号 𝑥𝑥 0 1/3 3/4 1 ① ② 𝑦𝑦 𝜑𝜑1 (𝑥𝑥) 0 1/3 3/4 1 3要素（4節点）の場合 𝑦𝑦 𝑢𝑢2 𝑑𝑑𝑑𝑑 −3 𝑑𝑑𝜑𝜑1 𝑥𝑥 𝑑𝑑𝑑𝑑 0 1/3 3/4 1 1 𝑥𝑥 𝑦𝑦 𝜑𝜑2 (𝑥𝑥) 0 1/3 3/4 1 𝑥𝑥 3 𝑦𝑦 𝑑𝑑𝜑𝜑2 𝑥𝑥 𝑑𝑑𝑑𝑑 12 − 5 0 1/3 3/4 1 𝑥𝑥 𝑦𝑦 1 𝜑𝜑3 (𝑥𝑥) 0 1/3 3/4 1 12 𝑥𝑥 5 −4 𝑦𝑦 𝑑𝑑𝜑𝜑3 𝑥𝑥 𝑑𝑑𝑑𝑑 0 1/3 3/4 1 𝑥𝑥 𝑥𝑥 𝑦𝑦 1 𝜑𝜑4 (𝑥𝑥) 0 1/3 3/4 1 4 𝑦𝑦 𝑑𝑑𝜑𝜑4 𝑥𝑥 𝑑𝑑𝑑𝑑 20 0 1/3 3/4 1 𝑥𝑥 𝑥𝑥

有限要素法（1次元）例題 1 𝑑𝑑𝜑𝜑𝑖𝑖 𝑥𝑥 𝑑𝑑𝜑𝜑𝑗𝑗 𝑥𝑥 𝑑𝑑𝑑𝑑 𝑑𝑑𝑑𝑑 𝑑𝑑𝑑𝑑 1/3 𝑎𝑎11 = ∫0 (−3)(−3) 𝑑𝑑𝑑𝑑 • 𝑎𝑎𝑖𝑖𝑖𝑖 = ∫0 • • =3 1 0 1/3 3/4 1 3要素（4節点）の場合 𝑦𝑦 𝑢𝑢2 𝑢𝑢1 𝑦𝑦 𝑢𝑢3 𝑢𝑢4 𝑢𝑢 𝑥𝑥 ③ ④: 節点番号 𝑥𝑥 0 1/3 3/4 1 ① ② 𝑦𝑦 𝜑𝜑1 (𝑥𝑥) −3 𝑑𝑑𝜑𝜑1 𝑥𝑥 𝑑𝑑𝑑𝑑 0 1/3 3/4 1 1 𝑥𝑥 𝑦𝑦 𝜑𝜑2 (𝑥𝑥) 0 1/3 3/4 1 𝑥𝑥 3 𝑦𝑦 𝑑𝑑𝜑𝜑2 𝑥𝑥 𝑑𝑑𝑑𝑑 12 − 5 0 1/3 3/4 1 𝑥𝑥 𝑦𝑦 1 𝜑𝜑3 (𝑥𝑥) 0 1/3 3/4 1 12 𝑥𝑥 5 −4 𝑦𝑦 𝑑𝑑𝜑𝜑3 𝑥𝑥 𝑑𝑑𝑑𝑑 0 1/3 3/4 1 𝑥𝑥 𝑥𝑥 𝑦𝑦 1 𝜑𝜑4 (𝑥𝑥) 0 1/3 3/4 1 4 𝑦𝑦 𝑑𝑑𝜑𝜑4 𝑥𝑥 𝑑𝑑𝑑𝑑 21 0 1/3 3/4 1 𝑥𝑥 𝑥𝑥

有限要素法（1次元）例題 1 𝑑𝑑𝜑𝜑𝑖𝑖 𝑥𝑥 𝑑𝑑𝜑𝜑𝑗𝑗 𝑥𝑥 𝑑𝑑𝑑𝑑 𝑑𝑑𝑑𝑑 𝑑𝑑𝑑𝑑 1/3 𝑎𝑎12 = ∫0 (−3)(3) 𝑑𝑑𝑑𝑑 • 𝑎𝑎𝑖𝑖𝑖𝑖 = ∫0 • • = −3 1 𝑢𝑢1 𝑦𝑦 𝑢𝑢3 𝑢𝑢4 𝑢𝑢 𝑥𝑥 ③ ④: 節点番号 𝑥𝑥 0 1/3 3/4 1 ① ② 1 0 1/3 3/4 1 3要素（4節点）の場合 𝑦𝑦 𝑢𝑢2 𝑦𝑦 𝜑𝜑1 (𝑥𝑥) −3 𝑑𝑑𝜑𝜑1 𝑥𝑥 𝑑𝑑𝑑𝑑 𝑥𝑥 0 1/3 3/4 1 𝑥𝑥 𝑦𝑦 𝜑𝜑2 (𝑥𝑥) 0 1/3 3/4 1 3 𝑦𝑦 𝑑𝑑𝜑𝜑2 𝑥𝑥 𝑑𝑑𝑑𝑑 12 − 5 0 1/3 3/4 1 𝑥𝑥 𝑦𝑦 1 𝜑𝜑3 (𝑥𝑥) 0 1/3 3/4 1 12 𝑥𝑥 5 −4 𝑦𝑦 𝑑𝑑𝜑𝜑3 𝑥𝑥 𝑑𝑑𝑑𝑑 0 1/3 3/4 1 𝑥𝑥 𝑥𝑥 𝑦𝑦 1 𝜑𝜑4 (𝑥𝑥) 0 1/3 3/4 1 4 𝑦𝑦 𝑑𝑑𝜑𝜑4 𝑥𝑥 𝑑𝑑𝑑𝑑 22 0 1/3 3/4 1 𝑥𝑥 𝑥𝑥

有限要素法（1次元）例題 1 𝑑𝑑𝜑𝜑𝑖𝑖 𝑥𝑥 𝑑𝑑𝜑𝜑𝑗𝑗 𝑥𝑥 𝑑𝑑𝑑𝑑 𝑑𝑑𝑑𝑑 • 𝑎𝑎𝑖𝑖𝑖𝑖 = ∫0 • 𝑎𝑎13 = 0 𝑑𝑑𝑑𝑑 1 𝑢𝑢1 𝑦𝑦 𝑢𝑢3 𝑢𝑢4 𝑢𝑢 𝑥𝑥 ③ ④: 節点番号 𝑥𝑥 0 1/3 3/4 1 ① ② 1 0 1/3 3/4 1 3要素（4節点）の場合 𝑦𝑦 𝑢𝑢2 𝑦𝑦 𝜑𝜑1 (𝑥𝑥) −3 𝑑𝑑𝜑𝜑1 𝑥𝑥 𝑑𝑑𝑑𝑑 𝑥𝑥 0 1/3 3/4 1 𝑥𝑥 𝑦𝑦 𝜑𝜑2 (𝑥𝑥) 0 1/3 3/4 1 3 𝑦𝑦 𝑑𝑑𝜑𝜑2 𝑥𝑥 𝑑𝑑𝑑𝑑 12 − 5 0 1/3 3/4 1 𝑥𝑥 𝑦𝑦 1 𝜑𝜑3 (𝑥𝑥) 0 1/3 3/4 1 12 𝑥𝑥 5 −4 𝑦𝑦 𝑑𝑑𝜑𝜑3 𝑥𝑥 𝑑𝑑𝑑𝑑 0 1/3 3/4 1 𝑥𝑥 𝑥𝑥 𝑦𝑦 1 𝜑𝜑4 (𝑥𝑥) 0 1/3 3/4 1 4 𝑦𝑦 𝑑𝑑𝜑𝜑4 𝑥𝑥 𝑑𝑑𝑑𝑑 23 0 1/3 3/4 1 𝑥𝑥 𝑥𝑥

有限要素法（1次元）例題 1 𝑑𝑑𝜑𝜑𝑖𝑖 𝑥𝑥 𝑑𝑑𝜑𝜑𝑗𝑗 𝑥𝑥 𝑑𝑑𝑑𝑑 𝑑𝑑𝑑𝑑 • 𝑎𝑎𝑖𝑖𝑖𝑖 = ∫0 • 𝑎𝑎14 = 0 𝑑𝑑𝑑𝑑 1 𝑢𝑢1 𝑦𝑦 𝑢𝑢3 𝑢𝑢4 𝑢𝑢 𝑥𝑥 ③ ④: 節点番号 𝑥𝑥 0 1/3 3/4 1 ① ② 1 0 1/3 3/4 1 3要素（4節点）の場合 𝑦𝑦 𝑢𝑢2 𝑦𝑦 𝜑𝜑1 (𝑥𝑥) −3 𝑑𝑑𝜑𝜑1 𝑥𝑥 𝑑𝑑𝑑𝑑 𝑥𝑥 0 1/3 3/4 1 𝑥𝑥 𝑦𝑦 𝜑𝜑2 (𝑥𝑥) 0 1/3 3/4 1 3 𝑦𝑦 𝑑𝑑𝜑𝜑2 𝑥𝑥 𝑑𝑑𝑑𝑑 12 − 5 0 1/3 3/4 1 𝑥𝑥 𝑦𝑦 1 𝜑𝜑3 (𝑥𝑥) 0 1/3 3/4 1 12 𝑥𝑥 5 −4 𝑦𝑦 𝑑𝑑𝜑𝜑3 𝑥𝑥 𝑑𝑑𝑑𝑑 0 1/3 3/4 1 𝑥𝑥 𝑥𝑥 𝑦𝑦 1 𝜑𝜑4 (𝑥𝑥) 0 1/3 3/4 1 4 𝑦𝑦 𝑑𝑑𝜑𝜑4 𝑥𝑥 𝑑𝑑𝑑𝑑 24 0 1/3 3/4 1 𝑥𝑥 𝑥𝑥

有限要素法（1次元）例題 1 𝑑𝑑𝜑𝜑𝑖𝑖 𝑥𝑥 𝑑𝑑𝜑𝜑𝑗𝑗 𝑥𝑥 𝑑𝑑𝑑𝑑 𝑑𝑑𝑑𝑑 𝑑𝑑𝑑𝑑 1/3 𝑎𝑎22 = ∫0 (3)(3) 𝑑𝑑𝑑𝑑 + 3/4 12 12 ∫1/3 (− 5 )(− 5 ) 𝑑𝑑𝑑𝑑 12 27 =3+ = 5 5 • 𝑎𝑎𝑖𝑖𝑖𝑖 = ∫0 • • 1 𝑢𝑢1 𝑦𝑦 𝑢𝑢3 𝑢𝑢4 𝑢𝑢 𝑥𝑥 ③ ④: 節点番号 𝑥𝑥 0 1/3 3/4 1 ① ② 1 0 1/3 3/4 1 3要素（4節点）の場合 𝑦𝑦 𝑢𝑢2 𝑦𝑦 𝜑𝜑1 (𝑥𝑥) −3 𝑑𝑑𝜑𝜑1 𝑥𝑥 𝑑𝑑𝑑𝑑 𝑥𝑥 0 1/3 3/4 1 𝑥𝑥 𝑦𝑦 𝜑𝜑2 (𝑥𝑥) 0 1/3 3/4 1 3 𝑦𝑦 𝑑𝑑𝜑𝜑2 𝑥𝑥 𝑑𝑑𝑑𝑑 12 − 5 0 1/3 3/4 1 𝑥𝑥 𝑦𝑦 1 𝜑𝜑3 (𝑥𝑥) 0 1/3 3/4 1 12 𝑥𝑥 5 −4 𝑦𝑦 𝑑𝑑𝜑𝜑3 𝑥𝑥 𝑑𝑑𝑑𝑑 0 1/3 3/4 1 𝑥𝑥 𝑥𝑥 𝑦𝑦 1 𝜑𝜑4 (𝑥𝑥) 0 1/3 3/4 1 4 𝑦𝑦 𝑑𝑑𝜑𝜑4 𝑥𝑥 𝑑𝑑𝑑𝑑 25 0 1/3 3/4 1 𝑥𝑥 𝑥𝑥

有限要素法（1次元）例題 1 𝑑𝑑𝜑𝜑𝑖𝑖 𝑥𝑥 𝑑𝑑𝜑𝜑𝑗𝑗 𝑥𝑥 𝑑𝑑𝑑𝑑 𝑑𝑑𝑑𝑑 𝑑𝑑𝑑𝑑 3/4 12 12 𝑎𝑎23 = ∫1/3 (− )( ) 𝑑𝑑𝑑𝑑 5 5 • 𝑎𝑎𝑖𝑖𝑖𝑖 = ∫0 • • =− 12 5 1 𝑢𝑢1 𝑢𝑢3 𝑢𝑢4 0 1/3 3/4 1 𝑢𝑢 𝑥𝑥 𝑦𝑦 𝑥𝑥 1 0 1/3 3/4 1 3要素（4節点）の場合 𝑦𝑦 𝑢𝑢2 𝑦𝑦 𝜑𝜑1 (𝑥𝑥) −3 𝑑𝑑𝜑𝜑1 𝑥𝑥 𝑑𝑑𝑑𝑑 𝑥𝑥 0 1/3 3/4 1 𝑥𝑥 𝑦𝑦 𝜑𝜑2 (𝑥𝑥) 0 1/3 3/4 1 3 𝑦𝑦 𝑑𝑑𝜑𝜑2 𝑥𝑥 𝑑𝑑𝑑𝑑 12 − 5 0 1/3 3/4 1 𝑥𝑥 𝑦𝑦 1 𝜑𝜑3 (𝑥𝑥) 0 1/3 3/4 1 12 𝑥𝑥 5 −4 𝑦𝑦 𝑑𝑑𝜑𝜑3 𝑥𝑥 𝑑𝑑𝑑𝑑 0 1/3 3/4 1 𝑥𝑥 𝑥𝑥 𝑦𝑦 1 𝜑𝜑4 (𝑥𝑥) 0 1/3 3/4 1 4 𝑦𝑦 𝑑𝑑𝜑𝜑4 𝑥𝑥 𝑑𝑑𝑑𝑑 26 0 1/3 3/4 1 𝑥𝑥 𝑥𝑥

有限要素法（1次元）例題 3 −3 27 −3 0 0 12 − 0 −1/6 −3/8 5 5 , 𝑓𝑓 = • A= 12 32 −1/3 0 − −4 5 5 −1/8 0 0 −4 4 1 • あとはこれを境界条件 𝑢𝑢1 = 0, 𝑢𝑢4 = 2 を満たすように解けばよい • この場合は、中央の二式を使って、 27 12 1 3 𝑢𝑢2 − 𝑢𝑢3 + 0 × = − −3 × 0 + 5 5 2 8 12 32 1 1 𝑢𝑢 + 𝑢𝑢 − 4 × = − 0×0+ 5 2 5 3 2 3 1 9 • これを解いて、𝑢𝑢2 = , 𝑢𝑢3 = 18 32 ― 解析解 ・数値解 27

有限要素法のポイント • 有限要素法で分割することで、 行列𝐀𝐀はスパース(疎)となる • 行列の行数をNとしたとき、 • →行列の非零成分の数はO(N)になる • →行列ベクトル積のコストはO(N)になる • 並列計算機上で大規模有限要素解析を実施する際には上記の特 性に加えて、対象計算システムの特徴を踏まえてアルゴリズ ム・実装を設計 • 単ノード特性：コア数・SIMD幅等、メモリバンド幅、メモリ容量など • システム特性：計算ノード数・大域通信・隣接通信のバンド幅・レイ テンシなど 28

大規模連立方程式の求解方法 29

ソルバー（連立方程式の解法） • 直接法 • ガウスの消去法・LU分解・コレスキー分解等 • メモリ使用量がO(N2)になり大規模問題の解析には不向き • 反復法 • ヤコビ法、共役勾配法等 • 疎行列の場合、メモリ使用量がO(N)となることが多い • ここでは地震関連問題を有限要素法で離散化した際に導かれる 正定値対称行列の大規模連立方程式の求解に広く使われる共役 勾配法を説明 30

共役勾配法の概念説明 （直接法で𝐀𝐀𝐮𝐮 = 𝐟𝐟を解く場合） • 非零ベクトル𝐩𝐩, 𝐪𝐪が𝐩𝐩T 𝐀𝐀𝐪𝐪 = 0となるとき、 𝐩𝐩, 𝐪𝐪は正定値対称行列𝐀𝐀に対し て共役であるという • 共役勾配法では、互いに共役な非零ベクトルの組 𝐩𝐩1 , 𝐩𝐩2 , 𝐩𝐩3 , … , 𝐩𝐩𝑁𝑁 (i.e., 𝐩𝐩T𝑖𝑖 𝐀𝐀𝐩𝐩𝑗𝑗 = 0 𝑖𝑖 ≠ 𝑗𝑗 )を基底として解を𝐮𝐮 = ∑𝑗𝑗 𝛼𝛼𝑗𝑗 𝐩𝐩𝑗𝑗 とあらわし、 𝐩𝐩T𝑖𝑖 𝐀𝐀𝐮𝐮 = 𝐩𝐩T𝑖𝑖 𝐀𝐀 ∑𝑗𝑗 𝛼𝛼𝑗𝑗 𝐩𝐩𝑗𝑗 = 𝛼𝛼𝑖𝑖 𝐩𝐩T𝑖𝑖 𝐀𝐀𝐩𝐩𝑖𝑖 = 𝐩𝐩T𝑖𝑖 𝐟𝐟 より係数𝛼𝛼𝑖𝑖 は以下のように求められる 𝐩𝐩T 𝑖𝑖 𝐟𝐟 𝛼𝛼𝑖𝑖 = T 𝐩𝐩𝑖𝑖 𝐀𝐀 𝐩𝐩𝑖𝑖 • ポイント • 行列ベクトル積とベクトル内積を多数回実行することで求解可能 • 𝐀𝐀自体の変更は不要(=疎行列性を崩さない) • 効率的に共役ベクトルの組を求めるよう改良したものが反復法バージョン の共役勾配法（次ページで説明） 31

共役勾配法（反復法）で𝐀𝐀𝐮𝐮 = 𝐟𝐟を解く • 現在（𝑖𝑖反復目）での解𝐮𝐮𝑖𝑖 に対し残差を以下のように定義する 𝐫𝐫𝑖𝑖+1 = 𝐟𝐟 − 𝐀𝐀𝐮𝐮𝑖𝑖 • ここで、探索ベクトルが𝑖𝑖反復目までの探索ベクトルの組()と共 役となるよう、𝐫𝐫𝑖𝑖+1 に最も近い方向で共役性を満たすようにと る 𝐩𝐩𝑖𝑖+𝟏𝟏 = 𝐫𝐫𝑖𝑖+1 − 𝐩𝐩T 𝑖𝑖 𝐀𝐀 𝐫𝐫𝑖𝑖+1 𝐩𝐩T 𝑖𝑖 𝐀𝐀 𝐩𝐩𝑖𝑖 𝐩𝐩𝑖𝑖 • ポイント：次の共役ベクトル𝐩𝐩𝑖𝑖+1 を生成するために、これまで の探索ベクトル 𝐩𝐩1 , 𝐩𝐩2 , 𝐩𝐩3 , … , 𝐩𝐩𝑖𝑖 の数𝑖𝑖に比例した操作が不要 32

共役勾配法（反復法）アルゴリズム • 正定値対称行列𝐀𝐀に対して𝐀𝐀𝐮𝐮 = 𝐟𝐟を解く 𝐫𝐫0 ⇐ 𝐟𝐟 − 𝐀𝐀𝐮𝐮0 𝐩𝐩0 ⇐ 𝐫𝐫0 do k=0,1,2,… 𝛼𝛼𝑘𝑘 ⇐ 𝐫𝐫𝑘𝑘T 𝐩𝐩𝒌𝒌 𝐩𝐩T 𝑘𝑘 𝐀𝐀𝐩𝐩𝒌𝒌 𝐮𝐮𝒌𝒌+𝟏𝟏 ⇐ 𝐮𝐮𝒌𝒌 + 𝛼𝛼𝑘𝑘 𝐩𝐩𝑘𝑘 𝐫𝐫𝒌𝒌+𝟏𝟏 ⇐ 𝐫𝐫𝒌𝒌 − 𝛼𝛼𝑘𝑘 𝐀𝐀𝐩𝐩𝑘𝑘 if(|𝐫𝐫𝒌𝒌+𝟏𝟏 |/|𝐛𝐛| < 𝜀𝜀) break 𝐫𝐫 T 𝐫𝐫 𝛽𝛽𝑘𝑘 ⇐ 𝑘𝑘+1T 𝒌𝒌+𝟏𝟏 𝐫𝐫𝑘𝑘 𝐫𝐫𝒌𝒌 𝐩𝐩𝒌𝒌+𝟏𝟏 ⇐ 𝐫𝐫𝒌𝒌+1 + 𝛽𝛽𝑘𝑘 𝐩𝐩𝑘𝑘 enddo • ここで |*|はノルム、𝐮𝐮0 は初期解、𝜀𝜀は閾値、求解結果は𝐮𝐮𝒌𝒌+𝟏𝟏 • ポイント • 最大でもN反復で収束する（実際には残差が小さくなる方向に近い方向で探索ベクトルを構築する ため、 大幅に少ない反復数で収束することが多い） • 計算は行列ベクトル積とベクトル内積のみ：これらを高速化すればよい 33

共役勾配法の並列化 34

共役勾配法の並列化 • ベクトル内積と行列ベクトル積を並列計算する • 良く使われる方法：要素ベースの分割（領域境界で節点は共有す る） • グラフ分割ソフトウェアを利用する(e.g., METIS) 分割前のメッシュ 分割されたメッシュ 内積： for(i=0;i<localnodes;++i){ localresult += a[i]*b[i]*dupli[i]; } MPI_ALLREDUCE(localresult,globalresult,MPI_SUM) // 領域内ではdupli[i]=1, 境界では共有領域数をmとし たとき、dupli[i]=1/m // 右の図の場合、 // 領域内部の点（黒）はdupli[i]=1 // 共有節点（赤）はdupli[i]=1/2 ※行列ベクトル積の並列化については次回 35

共役勾配法の前処理 36

共役勾配法における前処理 • 行列𝐀𝐀の性質を改善し、共役勾配法の反復回数を減らす • 最大固有値/最小固有値をなるべく小さくする（正定値対称行列では固 有値はすべて正） • 𝐌𝐌 ≈ 𝐀𝐀−1 のようなMを使い、 𝐀𝐀𝐮𝐮 = 𝐟𝐟を解く代わりに 𝐌𝐌𝐌𝐌 𝐮𝐮 = 𝐌𝐌𝐌𝐌 を解 く • 𝐌𝐌=𝐀𝐀−1 の場合は共役勾配法は一反復で終了し、逆に𝐌𝐌 = 𝐈𝐈(単位行列)の 場合は前処理が無い場合と同じとなる(反復数削減効果なし) • 計算コスト・メモリ使用量に対して、反復数削減効果の高い𝐌𝐌を設定 することとなる 37

共役勾配法における前処理(ヤコビ法) • ヤコビ法：行列𝐀𝐀−1 の近似として、𝐀𝐀の対角成分の逆数を用いる 方法 • 計算コストが低く、メモリ使用量は少ない • 性質がそこまで悪くない問題だと収束する • 多次元問題で未知数が節点毎に複数定義されている場合、これ らをブロックとして、ブロック対角行列の逆行列を、全体の 𝐀𝐀−1 に使う方法 • ヤコビ法より計算コスト・メモリ使用量は増加するものの、近似性能 は改善 38

共役勾配法における前処理(ILU) • Lower-Upper decomposition (LU分解) • 直接法。複数の入力fに対して、𝐮𝐮 = 𝐀𝐀−1 𝐟𝐟を高効率で計算できる • 疎行列Aに適用した場合、分解後のL, Uのfill-inパターン（非零成分位置）が 𝐀𝐀のfill-inパターンと異なる（成分が増える）ことがある→コスト高・メモリ 使用量が多くなる • 大局的な依存関係がある→効率的に並列計算するのは難しい • LU分解の近似としてのIncomplete LU decomposition (ILU) • ILU(0)：fill-inなし(𝐀𝐀で零となる成分は強制的に値を0とする) • ILU(1)：fill-inパターンは𝐀𝐀2 と同じ(𝐀𝐀2 で零となる成分は強制的に値を0とす る) • など • LU分解による求解よりも精度は劣るが、メモリ使用量・計算コストは抑えら れる • ただし、LUと同じく依存関係がありそのままでは並列計算に適さない→依存 関係をなくすための近似と組み合わせることが多い 39

共役勾配法における前処理(マルチグリッ ド法) • マルチグリッド法 • N x M (M<N)の行列Pを用いて、 𝐑𝐑𝐑𝐑𝐑𝐑 𝐮𝐮c = 𝐑𝐑𝐑𝐑 を𝐮𝐮c について（粗 く）解き、これを𝐮𝐮 ≈ 𝐏𝐏𝐮𝐮c として用いる • 上記を再帰的に行うことで、多段の疎化が可能になる • 共役勾配法が苦手とする長周波成分を効率的に求めることができる • 𝐑𝐑, 𝐏𝐏の設定方法により、幾何マルチグリッド、代数マルチグ リッドに分類される • 幾何マルチグリッド(Geometric Multigrid):𝐑𝐑, 𝐏𝐏に幾何拘束を用いる場 合(e.g., 構造格子グリッドにおいて、節点を一個飛ばしで省いたコース グリッドを使う等) • 代数マルチグリッド(Algebraic Multigrid): 幾何的情報以外の情報も 使って𝐑𝐑, 𝐏𝐏を設定した場合 40

大規模並列計算用の前処理 • ヤコビ系 • 大規模問題までスケーリングが容易 • ただし反復回数は多い • ILU系 • マトリクスを格納する必要あり • 並列計算のためには工夫をする必要あり • マルチグリッド系 • 幾何マルチグリッドは低コストだが適用できる問題が限られている • 代数マルチグリッドは一般問題にも適用可能だが、マトリクスを格納する必 要あり • このようなアルゴリズム特性と、システムの特性・対象問題の特性 を踏まえて前処理を選択することとなる 41

大規模解析用有限要素法の例 42

ソルバー例１：地震時の地盤増幅解析＠京コンピュータ Solving preconditioning matrix CG loop Inner coarse loop Solving preconditioning matrix Second ordered tetrahedron Computations of outer loop Equation to be solved (double precision) Outer loop Solve system roughly using CG solver (single precision) Inner fine loop Use as initial solution Linear tetrahedron Solve system roughly using CG solver (single precision) Use for preconditioner of outer loop Second ordered tetrahedron • Solve preconditioning matrix roughly to reduce number of CG loops • Use geometric multi-grid method to reduce cost of preconditioner • Use single precision in preconditioner to reduce computation & communication T. Ichimura et al., Implicit Nonlinear Wave Simulation with 1.08T DOF and 0.270T Unstructured Finite Elements to Enhance Comprehensive Earthquake Simulation, SC15 Gordon Bell Prize Finalist 43

ソルバー例１： Performance comparison on full K computer • Compare with PCGE & SC14 solver • PCGE: CG + Element-by-Element + simple preconditioner (3x3 block diagonal Jacobi preconditioner) • SC14 : SC14 Gordon Bell Prize Finalist • SC15: developed solver • FLOP (arithmetic) count and elapsed time reduced FLOP count (x1017 FLOP) 8 6.9 1000 763 800 6 4 3.0 2 0 Elapsed time (seconds) PCGE SC14 600 1.8 SC15 328 400 200 0 89.2 PCGE SC14 SC15 [1] T. Ichimura, et al., Physics-based urban earthquake simulation enhanced by 10.7 BlnDOF x 30 K time-step unstructured FE non-linear seismic wave simulation, SC14 Gordon Bell Prize Finalist [2] T. Ichimura, et al., Implicit Nonlinear Wave Simulation with 1.08T DOF and 0.270T Unstructured Finite Elements to Enhance Comprehensive Earthquake Simulation, SC15 Gordon Bell Prize Finalist 44

ソルバー例１： Scalability • Weak-scaling: 96.6% efficiency from 9,216 cores to 663,552 cores • 18.6% of peak performance attained at full K computer (=1.97 PFLOPS) • Strong-scaling: 76% efficiency from 9,216 cores to 294,912 cores • Fast & scalable solver algorithm + fast Element-by-Element method • Enables very good scalability & peak-performance & convergency & time-tosolution on K computer 2 1 Normalized elapsed time 9216 294912 663552 # of cores 1/2 1/4 1/8 1/16 1/32 Weak-scaling Strong-scaling Ideal speedup 45

ソルバー例２： Earth’s crust deformation problem • Compute static elastic response under faulting • Input: fault slip distribution • Output: deformation at surface 0 • Converged solutions of stress and strain required for evaluating plate boundary constitutive relations • Resolution of O(10 m) required at plate boundaries, leading to O(trillion) degreesof-freedom problem • Fast and scalable solver with small memory footprint required Earth’s crust deformation problem 46

ソルバー例２：地殻変動解析＠京コンピュータ Solving preconditioning matrix Inner loop level 2 CG loop Second ordered tetrahedron Computations of outer loop Use u1←P12 u2 as initial solution Solve system roughly using CG solver (single precision) Inner loop level 0 Use u0←P01 u1 as initial solution Linear tetrahedron Solve system roughly using CG solver (single precision) Use for preconditioner of outer loop Second ordered tetrahedron K. Fujita et al., Fast and Scalable Low-Order Implicit Unstructured Finite-Element Solver for Earth’s Crust Deformation Problem, 47 PASC2017 Geometric coarsening Equation to be solved (double precision) Outer loop Inner loop level 1 Coarsened equation (P12TA1P12)u2=P12Tf1 Algebraic coarsening Solving preconditioning matrix Solve system roughly using CG solver (single precision)

ソルバー例２： Size-up efficiency on K computer • Compare performance with PCGE (CG solver with 3x3 block diagonal preconditioner) and SC14 geometric multi-grid solver 1800 9,216 18,432 36,864 73,728 147,456 294,912 # of compute cores 1400 146.5 s (20.6%) 319.7 s (20.4%) 1629 s (17.6%) 400 148.7 s (20.1%) 308.5 s (20.5%) 1573 s (17.9%) 600 151.3 s (20.1%) 311.8 s (21.3%) 1680 s (18.3%) 800 144.5 s (20.9%) 292.7 s (21.4%) 1572 s (18.1%) 1000 140.0 s (21.2%) 277.5 s (21.6%) 1546 s (18.1%) 1200 138.9 s (21.3%) 272.0 s (21.9%) 1375 s (18.5%) elapsed time of solver [s] 1600 Floating-point arithmetic efficiency to peak • 94.8% scalability from 9216 to 294912 cores Elapsed time 200 0 48

ソルバー例２： Speed-up efficiency on K computer • Good speed-up efficiency when compared with previous solvers elapsed time of solver [s] 2048 512 128 1375 s (18.5%) 272 s (21.9%) 364 s (17.5%) 178 s (21.0%) 72.2 s (20.3%) 9,216 229 s (16.6%) 73.7 s (20.3%) 139 s (21.3%) 32 8 800 s (17.7%) 18,432 38.2 s (19.4%) 36,864 111 s (16.1%) 46.3 s (19.6%) 21.8 s (17.6%) 73,728 Elapsed time Floating-point arithmetic efficiency to peak 27.2 s (17.1%) 13.2 s (14.6%) 147,456 # of compute cores 49

ソルバー例２： Performance for practical problem • Apply to East Japan model with eight crust layers with full 3D crust geometry and surface topography, and input fault slip outer inner level 0 inner level 1 inner level 2 750 500 250 0 2.20 x faster Developed GAMERA 4966 s (17.8%) 37373 ---- ---- ---- 1000 555.8 s (18.4%) 10 225 19808 ---- Overview of 381,383,242,611 degreesof-freedom unstructured finite-element model with mesh size 125 m 19.7 x faster 252.2 s (13.2%) 16 228 4350 17266 0 elapsed time of solver [s] 6000 5000 4000 # of 3000 iterations FP efficienc y Elapse d time PCGE Performance on 294,912 compute cores of K computer 2000 1000 0 50

Comparison of stress response • Stress converged at plate boundaries Cross section in b) 0 2.25x107 1.0x10-4 (mesh size 125 m) (mesh size 125 m) Structured Unstructured Crust layer interface 7.5x10-1 Pa (in log scale) 500m a) Overview of unstructured finite-element model (mesh size 250 m) (mesh size 500 m) Unstructured Unstructured b) Stress response at cross section y = 300 km 51

まとめ • 地震シミュレーションにおける有限要素法の定式化・求解方法 を概説 • 大規模有限要素解析はものづくりなど幅広い他分野で使われており、 これらの分野にも活用可能 • 大規模高速解析を実現することで、数値解の収束確認や多数 ケースの解析など、解析の信頼性向上につながると期待 • 次回(6/27)は今回説明した大規模有限要素法の高効率な実装に 関して説明する 70