---
title: 第10回  配信講義 計算科学技術特論B (2024)
tags: 
author: [R-CCS　計算科学研究推進室](https://www.docswell.com/user/2300203199)
site: [Docswell](https://www.docswell.com/)
thumbnail: https://bcdn.docswell.com/page/DEYGLL88EM.jpg?width=480
description: 第10回 6月20日　大規模地震シミュレーション１ 本講義では科学技術計算の一例として、地震分野における大規模有限要素法シミュレーション手法について説明する。ここでは、有限要素法の基礎、大規模線形連立方程式の解法の一種である共役勾配法、及び、共役勾配法の並列計算方法について説明する。
published: June 14, 24
canonical: https://www.docswell.com/s/2300203199/KVVGD2-2024-06-14-130229
---
# Page. 1

![Page Image](https://bcdn.docswell.com/page/DEYGLL88EM.jpg)

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


# Page. 2

![Page Image](https://bcdn.docswell.com/page/VJNG44Q978.jpg)

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


# Page. 3

![Page Image](https://bcdn.docswell.com/page/YE9WQQD3E3.jpg)

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


# Page. 4

![Page Image](https://bcdn.docswell.com/page/GE8NGG3L7D.jpg)

地震シミュレーション概要
• 断層～地殻～地盤～構造物
• 領域サイズ: 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


# Page. 5

![Page Image](https://bcdn.docswell.com/page/LELPGGVQ7R.jpg)

解析例
• ターゲット：断層～地殻～地盤～構造物～社会活動
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


# Page. 6

![Page Image](https://bcdn.docswell.com/page/4JM4QQNK7W.jpg)

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


# Page. 7

![Page Image](https://bcdn.docswell.com/page/PJRW88D679.jpg)

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


# Page. 8

![Page Image](https://bcdn.docswell.com/page/PEXP881D7X.jpg)

有限要素法の基礎
16


# Page. 9

![Page Image](https://bcdn.docswell.com/page/3EKDKK2DJD.jpg)

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


# Page. 10

![Page Image](https://bcdn.docswell.com/page/L73RZZ4PJ5.jpg)

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


# Page. 11

![Page Image](https://bcdn.docswell.com/page/87DYRRQ37G.jpg)

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


# Page. 12

![Page Image](https://bcdn.docswell.com/page/VJP6WWLPE8.jpg)

有限要素法（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
𝑥𝑥
𝑥𝑥


# Page. 13

![Page Image](https://bcdn.docswell.com/page/2EV388QV7Q.jpg)

有限要素法（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
𝑥𝑥
𝑥𝑥


# Page. 14

![Page Image](https://bcdn.docswell.com/page/57G355W1JL.jpg)

有限要素法（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
𝑥𝑥
𝑥𝑥


# Page. 15

![Page Image](https://bcdn.docswell.com/page/4EQ5ZZ3NJP.jpg)

有限要素法（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
𝑥𝑥
𝑥𝑥


# Page. 16

![Page Image](https://bcdn.docswell.com/page/KJ49331371.jpg)

有限要素法（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
𝑥𝑥
𝑥𝑥


# Page. 17

![Page Image](https://bcdn.docswell.com/page/LE1K11GZJG.jpg)

有限要素法（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
𝑥𝑥
𝑥𝑥


# Page. 18

![Page Image](https://bcdn.docswell.com/page/GEWD88K6J2.jpg)

有限要素法（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
𝑥𝑥
𝑥𝑥


# Page. 19

![Page Image](https://bcdn.docswell.com/page/47ZG88ZRE3.jpg)

有限要素法（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


# Page. 20

![Page Image](https://bcdn.docswell.com/page/YJ6YPPZ17V.jpg)

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


# Page. 21

![Page Image](https://bcdn.docswell.com/page/GJ5GKKWLJ4.jpg)

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


# Page. 22

![Page Image](https://bcdn.docswell.com/page/LE3RZZ4D75.jpg)

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


# Page. 23

![Page Image](https://bcdn.docswell.com/page/8EDYRRQ2EG.jpg)

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


# Page. 24

![Page Image](https://bcdn.docswell.com/page/V7P6WWLLJ8.jpg)

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


# Page. 25

![Page Image](https://bcdn.docswell.com/page/2JV388Q6EQ.jpg)

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


# Page. 26

![Page Image](https://bcdn.docswell.com/page/5EG355W27L.jpg)

共役勾配法の並列化
34


# Page. 27

![Page Image](https://bcdn.docswell.com/page/4JQ5ZZ397P.jpg)

共役勾配法の並列化
• ベクトル内積と行列ベクトル積を並列計算する
• 良く使われる方法：要素ベースの分割（領域境界で節点は共有す
る）
• グラフ分割ソフトウェアを利用する(e.g., METIS)
分割前のメッシュ
分割されたメッシュ
内積：
for(i=0;i&lt;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


# Page. 28

![Page Image](https://bcdn.docswell.com/page/K74933NRE1.jpg)

共役勾配法の前処理
36


# Page. 29

![Page Image](https://bcdn.docswell.com/page/LJ1K115V7G.jpg)

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


# Page. 30

![Page Image](https://bcdn.docswell.com/page/GJWD885372.jpg)

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


# Page. 31

![Page Image](https://bcdn.docswell.com/page/4EZG88NGJ3.jpg)

共役勾配法における前処理(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


# Page. 32

![Page Image](https://bcdn.docswell.com/page/Y76YPP8QEV.jpg)

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


# Page. 33

![Page Image](https://bcdn.docswell.com/page/G75GKK9274.jpg)

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


# Page. 34

![Page Image](https://bcdn.docswell.com/page/9J2YWW21JR.jpg)

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


# Page. 35

![Page Image](https://bcdn.docswell.com/page/DEYGLLY4EM.jpg)

ソルバー例１：地震時の地盤増幅解析＠京コンピュータ
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 &amp; 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


# Page. 36

![Page Image](https://bcdn.docswell.com/page/VJNG44P678.jpg)

ソルバー例１： Performance comparison on full K computer
• Compare with PCGE &amp; 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


# Page. 37

![Page Image](https://bcdn.docswell.com/page/YE9WQQ3RE3.jpg)

ソルバー例１： 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 &amp; scalable solver algorithm + fast Element-by-Element method
• Enables very good scalability &amp; peak-performance &amp; convergency &amp; 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


# Page. 38

![Page Image](https://bcdn.docswell.com/page/GE8NGGMG7D.jpg)

ソルバー例２： 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


# Page. 39

![Page Image](https://bcdn.docswell.com/page/LELPGG3X7R.jpg)

ソルバー例２：地殻変動解析＠京コンピュータ
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)


# Page. 40

![Page Image](https://bcdn.docswell.com/page/4JM4QQD37W.jpg)

ソルバー例２： 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


# Page. 41

![Page Image](https://bcdn.docswell.com/page/PJRW884R79.jpg)

ソルバー例２： 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


# Page. 42

![Page Image](https://bcdn.docswell.com/page/PEXP882Y7X.jpg)

ソルバー例２： 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


# Page. 43

![Page Image](https://bcdn.docswell.com/page/3EKDKKM4JD.jpg)

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


# Page. 44

![Page Image](https://bcdn.docswell.com/page/L73RZZYDJ5.jpg)

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


