講義１：離散化と並列化

201４年３月５日
神戸大学大学院システム情報学研究科計算科学専攻
横川三津夫

講義の内容計算機シミュレーション
非定常問題の時間微分項の離散化
定常問題
線型方程式の反復解法
実行時間の評価

2014/3/5

RIKEN AICS HPC Spring School 2014：離散化とその周辺

1

計算機シミュレーション
現象
現象のモデル化
数理モデル

偏微分方程式など

離散化，アルゴリズム

離散化（計算格子も選択）

計算

プログラム，コンピュータ•有限差分法，有限要素法，スペクトル法など•足し算，掛け算だけの式で記述

結果の検証
2014/3/5

RIKEN AICS HPC Spring School 2014：離散化とその周辺

2

熱伝導問題（非定常）20cm• 長さ20ｃｍのアルミニウム棒の左端を温めたとき，1分後の右端の温度は何度になっているのか？ 棒の温度がどんな風に変化するかを知りたい．
▫ 最初の棒の温度は20℃（初期状態）▫ 左端を加熱し，温度100℃にする．
▫ 熱伝導度 9.7x10-5m2/s

長さ方向の温度変化だけを見れば良い． ⇒ 1次元問題

2014/3/5

RIKEN AICS HPC Spring School 2014：離散化とその周辺

3モデル化（1次元非定常熱伝導問題）

20cm
座標 x02468101214161820• 長さ20ｃｍのアルミニウム棒の左端を温めたとき，1分後の右端の温度は何度になっているのか？ 棒の温度がどんな風に変化するかを知りたい．
▫ 最初の棒の温度は20℃（初期状態）▫ 左端を加熱し，温度100℃にする．
▫ 熱伝導度 9.7x10-5m2/s

長さ方向の温度変化だけを見れば良い． ⇒ 1次元問題
1次元の数理モデル（偏微分方程式）で表現
1次元非定常熱伝導方程式

2014/3/5uua 2tx

RIKEN AICS HPC Spring School 2014：離散化とその周辺

24離散化：連続な量を有限個の点で表現することu

番号 jj =01234x  2t0x  j x u

56u
u
a 2
txut0j7x89105  x  10u tj0  t  u tj0

2t0  Δtj

定義

 u  r (ut0j

離散化t

時間微分：前進オイラー法
空間微分：有限差分法

t0
j 1

 2u  u
t0
jt0
j 1)

ut0  Δtj

au tj01  2u tj0  u tj01(x) 2

att0
t0t0
u (u2uuj 1jj 1 )
2(x)
t0
j=r

2014/3/5

RIKEN AICS HPC Spring School 2014：離散化とその周辺

5

時刻ｔ0 の温度から⊿ｔ 経過した後の温度を求める式

ut0  Δtj

 u  r (ut0
j12u

時刻 t0 + ⊿ｔ0
2014/3/5

 2u  ut0
jt0
j 1

乗算： 2回加減算： 3回)

j =3に注目

時刻 t00t0
j 1

1

2

3t0
2

u

4t0
3

uut0  Δt
3

3

4

5

6

7

8

9105678910t0
4

RIKEN AICS HPC Spring School 2014：離散化とその周辺

6

プログラミング言語FORTRAN77によるプログラム
real*8格子分割数 20
格子間隔 ⊿x=1cm
時間間隔 ⊿ｔ=0.5秒
熱伝導度 a=9.7ｘ10‐5

初期温度を設定する．

n = 20dx = 0.2/real(n)dt = 0.5
a = 9.7e‐5
r = a*dt/(dx*dx)
do j = 0, 20
u_old(j) = 20.0
enddo

左端を100度にする．

u_old(0)120ｘ0.5秒=60秒まで繰り返し

do

各格子点の値を求める．

結果を書き出す．

u_new(0:20), u_old(0:20)

u tj0  Δt  u tj0  r (u tj01  2u tj0  u tj01 )= 100.0

it = 1, 120
do j = 1, 19
u_new(j) = (1.0–2.0*r)*u_old(j) + r*(u_old(j‐1)+u_old(j+1))
enndo
do j = 1, 19
u_old(j) = u_new(j)enddou_old[20] = u_old[19]write(6,*) dt*it, (u_old(j), j=0,20)
enddo
stopend

2014/3/5

RIKEN AICS HPC Spring School 2014：離散化とその周辺

7

0秒から60秒までの温度変化（10秒間隔）
100

t=0秒
t=10秒
t=20秒

温度（℃）

80

t=30秒

49℃

60

t=40秒
t=50秒
t=60秒

30℃

40

200
0.0

アルミ棒の左端
2014/3/5

5.0

10.0

15.0

位置（cm）RIKEN AICS HPC Spring School 2014：離散化とその周辺

20.0

アルミ棒の右端
8

非定常問題（時間発展方程式）に対する時間微分項の離散化：
初期条件

は空間の微分演算子を含む項，境界条件

境界条件が与えられる．

,2∆
∆
初期条件

,"
2014/3/5

RIKEN AICS HPC Spring School 2014：離散化とその周辺

9時間発展（非定常）方程式に対する時間微分項の離散化，，は時間刻み幅（時間ステップ）

陽解法
過去の情報だけで，次の時刻を求める方法

時刻 の値だけを用いて，時刻の値を求める方法

次の時刻の値を求めるために，空間の微分式が成立するとして求める方法

・後退オイラー法など∆∆

・ルンゲ・クッタ法など
・アダムス-バッシュフォース法

多段法

2014/3/5の値

・前進オイラー法

単段法

複数の過去の時刻 ，⋯ の値を用いて，時刻の値を求める方法

陰解法，・リープ-フロッグ法など

・アダムス-モールトン法など

2∆
RIKEN AICS HPC Spring School 2014：離散化とその周辺

10

定常問題（熱伝導方程式を例として）

熱伝導方程式,の定常状態の解（ ⟶ ∞の時の解）はどうなるか？

＆境界条件【境界条件】• Dirichlet（ディリクレ）問題： 対象となる領域の境界上で解の値が指定される場合• Neumann（ノイマン）問題： 境界上で法線方向の微分が指定されている
場合
2014/3/5

RIKEN AICS HPC Spring School 2014：離散化とその周辺

112次元定常熱伝導問題• 2次元正方形領域 [0,1]×[0,1] での熱伝導問題• 境界条件

2014/3/50,0011,001,0001,1sin01

RIKEN AICS HPC Spring School 2014：離散化とその周辺

12差分式のテイラー級数展開からの導出（2変数）,,,,,,

⋯

2!
,,

⋯2!

これらの式から．．．
,,,,2,,,
2014/3/5,,,,,2,,RIKEN AICS HPC Spring School 2014：離散化とその周辺

132次元定常熱伝導問題の離散化（その１）1 等分：

方向を1/1 ， 方向を1 等分：

1/1,

境界条件として既知の値,内部の格子点（未知数）,内部の格子点での関係式を作る．
,,,,,2014/3/5,,,,,RIKEN AICS HPC Spring School 2014：離散化とその周辺

14

離散化
偏微分方程式の中の偏導関数をその差分近似で置き換えることにより，偏微分方程式が連立一次方程式になる．

を，点,,,とすると2014/3/5で離散化,,,,,,,,RIKEN AICS HPC Spring School 2014：離散化とその周辺,,15

, を辞書式順序で並べたベクトルを,,4,とすると，点 , に関する離散式は,,,0,,⋮,,⋮,,,,,

⋮

⋮

,,

⋮,

⋮
⋮

,0⋯

1

⋯ 1離れている,⋮

4 1⋯ 1 0離れている⋮⋮⋮,0,,,,,,⋮,

2014/3/50,

⋮⋮

,

⋮⋯

,⋮,

⋮
⋮

,,,

RIKEN AICS HPC Spring School 2014：離散化とその周辺

16

行列表現
4 -1-1 4 -1-1 4 -1-1-1,-1,⋮-1 4 -1-1 4-1,,,-1,⋮-1-14 -1-1 4 -1-1 4 -1-1-1-1,,

⋮
⋮
-1

⋮,,-1

-1 4 -1-1 4-1-1-1-1-1-1-1,

⋮
⋮

4 -1
-1 4 -1-1 4 -1-1
2014/3/5⋮,,-1 4 -1-1 4RIKEN AICS HPC Spring School 2014：離散化とその周辺⋮,,17

1次元問題に対する行列,次連立一次方程式が導かれる

行列形式：

2014/3/5

RIKEN AICS HPC Spring School 2014：離散化とその周辺

18線形方程式とその解法＝:

正方行列，正則かつ＝， ＝: 次元ベクトル次線形方程式＝ ，またはを解く．

直接法： 行列操作により厳密解を求める方法
ガウス消去法
LU分解法（対称行列用 コレスキー法）など
2014/3/5

RIKEN AICS HPC Spring School 2014：離散化とその周辺

19

連立一次方程式の反復解法

反復法
を真の解 ∗ の近似値とし，反復ベクトルの系列→→→⋯→により，真の解に近づける方法
反復ベクトルの系列の作り方によりいろいろな解法がある．
定常反復法
ヤコビ法（Jacobi法）
ガウス・ザイデル法（Gauss-Seidel法）
逐次過大緩和法（Successive Over-Relaxation：SOR）法
交互方向法（Alternating-Direction Implicit Iterative Method: ADI法）など

クリロフ（Krylov）空間法
共役勾配法（Conjugate Gradient Method，CG法）→ 講義4

双共役勾配法（Bi-CG法）
一般化最小残差法（GMRES法）など
2014/3/5

RIKEN AICS HPC Spring School 2014：離散化とその周辺

20

行列分離

００００
行列分離：

＝: 対角行列: 狭義下三角行列: 狭義上三角行列

2014/3/5

RIKEN AICS HPC Spring School 2014：離散化とその周辺

21

(1) ヤコビ法（Jacobi Method）

ヤコビ行列

で，反復ベクトル系列を生成する．

/,

1

2014/3/5,0

RIKEN AICS HPC Spring School 2014：離散化とその周辺

22

(2) ガウス・ザイデル法（Gauss-Seidel Method）

ガウス・ザイデル行列/2014/3/5

1

RIKEN AICS HPC Spring School 2014：離散化とその周辺,023

(3) SOR法（Successive Over-Relaxation Method）

1,0

これらを整理すると，

11SOR行列2014/3/5

ω

RIKEN AICS HPC Spring School 2014：離散化とその周辺

24熱伝導問題への反復法の適用（ヤコビ法）

ヤコビ法
,,,,,

⁄4,,,

⁄4

→ 演習１

ガウス・ザイデル法,
2014/3/5,

RIKEN AICS HPC Spring School 2014：離散化とその周辺

25

ヤコビ法

,から,を計算するので，どの点から計算しても同じ結果になる．
領域分割によって並列化が可能．

2014/3/5

RIKEN AICS HPC Spring School 2014：離散化とその周辺

26

ヤコビ法

プロセス1

プロセス0

2014/3/5

RIKEN AICS HPC Spring School 2014：離散化とその周辺

27

ガウス・ザイデル法,の計算には，更新された,，,が必要．このままでは並列化出来ない．
2014/3/5

RIKEN AICS HPC Spring School 2014：離散化とその周辺

28

ガウス・ザイデル法

プロセス1プロセス0

計算の順番に着目して，並列化可能なプログラムを見つける．
マルチカラー法2014/3/5

RIKEN AICS HPC Spring School 2014：離散化とその周辺

29

ヤコビ法のMPI並列化の方針（その１）• 領域分割によりMPI並列化（プロセス並列化）を行う．
▫ 各プロセスは，部分の計算を担当する．
境界は条件として与えられているので，計算しないことに注意

2014/3/5

RIKEN AICS HPC Spring School 2014：離散化とその周辺30ヤコビ法のMPI並列化の方針（その２）4 -1-1 4 -1-1 4 -1-1-1,-1,⋮-1 4 -1-1 4-1,,,-1,⋮-1-14 -1-1 4 -1-1 4 -1-1-1-1,,

⋮
⋮
-1

⋮,,

-1-1 4 -1-1 4-1-1-1-1-1-1-1,⋮⋮

4 -1-1 4 -1
-1 4 -1-1
2014/3/5⋮,,-1 4 -1-1 4

RIKEN AICS HPC Spring School 2014：離散化とその周辺⋮,,31

ヤコビ法のMPI並列化の方針（その３）• 各プロセスでの計算には，上下の部分の値,,4,,,が必要であることに注意．

0▫ スレッド並列では，の部分に同時にアクセスしないようにする．
▫ プロセス並列では，上下部分の値を転送する必要がある．

2014/3/5

RIKEN AICS HPC Spring School 2014：離散化とその周辺

32

実行時間の評価
並列実行数とともに，実行時間が減少することが期待される．ただし，並列化により，計算結果が異ならないよう注意する．
並列数 n に対し，計算時間は 1/n になることが期待
実行時間を計測して確認．

スケーラビリティ
弱スケーリング（weak scaling）
並列実行単位（スレッド，またはプロセス）あたりの問題サイズを一定に保ったまま（したがって問題のサイズは大きくなる），並列実行数（スレッド数，プロセス数）を増加させた時に，実行時間がどのように変化するかをみる性能評価指標
強スケーリング（strong scaling）
全体の問題サイズを固定したまま，並列実行数（スレッド数，プロセス数）を増加させた時に，実行時間がどのように変化するかをみる性能評価
指標
2014/3/5

RIKEN AICS HPC Spring School 2014：離散化とその周辺

33

弱スケーリング
並列実行単位（スレッド，またはプロセス）あたりの問題サイズを一定に保ったまま，並列実行数（スレッド数，プロセス数）を増加させた時に，実行時間がどのように変化するかをみる性能評価指標

1プロセス
2014/3/5

2プロセス

4プロセス

RIKEN AICS HPC Spring School 2014：離散化とその周辺

8プロセス
34

強スケーリング全体の問題サイズを固定したまま，並列実行数（スレッド数，プロセス数）を増加させた時に，実行時間がどのように変化するかをみる性能評価指標

2プロセス

2014/3/5

4プロセス

RIKEN AICS HPC Spring School 2014：離散化とその周辺

8プロセス

35