局所細分化メッシュに基づく並列有限要素法における前処理付き反復法
Preconditioned Iterative Methods for Parallel Finite-Element
Applications with Adaptive Mesh Refinement

中島 研吾（1）

兵藤 守（2）（1）東京大学大学院 理学系研究科 地球惑星科学専攻（2）地球シミュレータセンター 固体地球シミュレーション研究グループ
日本応用数理学会「行列・固有値問題の解法とその応用」研究部会 第2回研究会
2006年11月16日，国立情報学研究所

2

本発表の概要• 背景– 地震発生サイクルシミュレーション– 局所細分化，負荷分散• 問題設定• 計算結果• まとめ– 今後の予定等

JSIAM06-NOV

3

地震発生サイクルシミュレーション• プレート境界における準静的応力蓄積過程• もともとのアプローチ– 非線形接触問題をNewton-Raphson法によって解く– ALM法（Augmented Lagrangean,拡大ラグランジェ法）による拘束条件：ペナルティ数– 領域分割による並列有限要素法• GeoFEM，地球シミュレータ• 前処理付き反復法 [Nakajima, 2003]，[Nakajima, 2006]

JSIAM06-NOV

4現状• 2003年以降，全てをFEMで解こうという考え方から方針転換– 接触モデリングの困難さ– 境界条件• 伝統的な境界積分法，境界要素法を使用する– 単位すべり応答関数（ある場所ξ，t=0における単位すべり量が時間t-s, 場所xにおける応力変化を与える）による寄与の積分（重ね合わせ）– 不均質場における「単位すべり応答関数」を計算する必要がある• 有限要素法の使用

JSIAM06-NOV

5FEMによる「単位すべり応答関数」計算• GeoFEM• 分割節点法（Split Node Method）による，食い違い変位（Dislocation）モデル– 断層をはさんだ変位の食い違い量を仮定し，仮想変位に相当する節点力を外力項として与える– 特殊な要素は不要，定式化も簡単，観測値の利用

JSIAM06-NOV

6FEMによる「単位すべり応答関数」計算（続き）• 不均質粘弾性モデル– Maxwell粘弾性体εσdε 1 dσ 1=+ σdt E dt η• ちなみに弾性体は・・・

tt

1ε= σE
JSIAM06-NOV

7

局所細分化の必要性• 1km×1kmの応答関数を求めるためには，数十m幅のメッシュが必要• 数百km×数百kmの領域を扱うため，全領域をこの大きさのメッシュで切ることは不可能– 局所細分化：Adaptive Mesh Refinement （AMR）

JSIAM06-NOV

8六面体要素のAMR• 実は六面体要素のAMRは難しい• 八分木（二次元では四分木）は細分化レベルの境界でHanging Nodeを生じてしまう– 四面体，三角形ではこのようなことはないが，断層のモデリングには六面体要素が望ましく，六面体要素が一般的に用いられている

JSIAM06-NOV

Refinement Patterns of Tetrahedron
3 edges ona same
surface aremarked

1 edgemarked2 children(Binary)
2 edges ona samesurface aremarkedmore than
2 edges on
different
surfaces aremarkedmore than
4 edges aremarkedJSIAM06-NOV

3rd edge is
automaticallymarked4 children(Quadtree)

6 edges areautomaticallymarked8 children(Octree)910

対処法：27分木 [Wada, 2000]• 各辺の分割単位を2ではなく，3とする。
– 細分化レベル境界でHanging Nodeを生じない– 二次元：9分木– 三次元：27分木JSIAM06-NOV

1127分木 における分割パターン[Wada, 2000]

JSIAM06-NOV

1227分木 による細分化例
初期レベル

JSIAM06-NOV

1327分木 による細分化例
細分化レベル1

JSIAM06-NOV

1427分木 による細分化例
細分化レベル2

JSIAM06-NOV

1527分木 による細分化例
細分化レベル3

JSIAM06-NOV

1627分木 による細分化例
細分化レベル3（拡大図）

JSIAM06-NOV

1727分木• 膨大な細分化メッシュ• 8000要素 ⇒ 1つの要素を3レベルの分割しただけで67,000要素以上– となりあうメッシュの分割レベルの差は1以下とする場合:精度• これまであまり使われなかった主な理由

JSIAM06-NOV

18「並列」27分木• 要素数，節点数が100万を越えるような初期メッシュを局所的に細分化する必要がある– 細分化レベルを増やしていくと，1プロセッサで扱うことは困難• 色々な細分化パターンを実施する必要がある– 「すべり応答関数」を算出するのが目的であるため，多数のメッシュパターンを「高速に」生成する必要がある• 並列計算– 初期メッシュ，領域分割– 「並列」細分化– 動的負荷分散（Dynamic LoadBalancing）＋データマイグレー
ション

マントル
四国

フィリピン海プレートzyx伊豆太平洋プレート
マントル

JSIAM06-NOV

Supersonic Flow around a Sphere
M=1.40 Uniform Flow, Ideal Gas, Re=106
4 PE cases (tetrahedron only)

PE0
PE1PE2PE3
JSIAM06-NOV

Initial Grid

1-Lev. Adapted

2-Lev. Adapted

before DRAMA
137
137
136
136
-

before
793
696
668
448

before
3834
2769
2703
1390

DRAMA
652
650
652
651

DRAMA
2527
2526
2522
2524
19

20

「並列」27分木の開発• 並列細分化– Double Numberingの導入 [Nakajima & Fingberg, 2001]– 局所的な細分化情報の通信（未実装）– 四面体と比較すると膨大な組み合わせを考慮する必要がある• ParMETISによる動的負荷分散– http://glaros.dtc.umn.edu/gkhome/metis/parmetis/– マルチレベル的な手法，節点数がバランスするような手法• データマイグレーション– 節点座標，要素コネクティビティ，物性値，グループ情報，通信テー
ブル等の領域間の移動（実装済み，最適化はしていない）– 領域間オーバーラップ深さの自動変更機能
JSIAM06-NOV

21

METIS
http://www-users.cs.umn.edu/~karypis/metis/• マルチレベルグラフ理論に基づいた方法

JSIAM06-NOV

22

METIS
http://www-users.cs.umn.edu/~karypis/metis/• マルチレベルグラフ理論に基づいた方法– 特に通信（edge-cut）が少ない分割を提供する– 安定，高速– フリーウェア，他のプログラムに組み込むことも容易• 色々な種類がある– k-METIS
通信量（edge-cut）最小– p-METIS
領域間バランス最適化– ParMETIS
並列版– 領域分割だけでなく，オーダリング，データマイニングなど色々な分野に使用されている• 接触，衝突問題における並列接触面探索
JSIAM06-NOV

23• 背景– 地震発生サイクルシミュレーション– 局所細分化，負荷分散• 問題設定• 計算結果• まとめ– 今後の予定等

JSIAM06-NOV

24

懸 念
• 細分化によって，要素の幅が大きく分布するため，反復法の収束が悪化するのではないか・・・• ParMETISを使用して，ブラックボックス的に領域分割を実施すると，並列時の収束が悪化するのではないか・・・

JSIAM06-NOV

25

例 題 （1/3）• 三次元弾性問題• 初期メッシュ– 1領域 203要素，64領域（約150万自由度）– z方向• 上10層：弾性体（リソスフェア，地殻）• 下10層：粘弾性体（アセノスフェア，マントル）– 剛性率の低い弾性体（ヤング率：小さい，ポアソン比：大きい）としてモデル化– ｔ→∞で成立する近似 [Fukahata, 2006]

リソスフェア（地殻）弾性体
アセノスフェア（マントル）
粘弾性体
JSIAM06-NOV

26例 題 （2/3）• 細分化メッシュ– 4領域に1領域の割合で，1要素を3段階に細分化– 70,000メッシュ近い領域が出てくる– 全体として，約150万要素（節点），450万自由度（正確には1,480,960節点，4,442,880自由度）– ParMETISによる動的負荷分散＋データマイグレーション

リソスフェア（地殻）弾性体
アセノスフェア（マントル）
粘弾性体
JSIAM06-NOV

27

例 題 （3/3）• 境界条件– 右図参照– 現実的なもの（食い違い変位）ではないzUniform Distributed Force inz-dirrection @ z=Zmin
Uy=0 @ y=Ymin• 物性値

Ux=0 @ x=Xmin– リソスフェア• ヤング率 1.00，ポアソン比：0.24

Nz-1
Ny-1

– アセノスフェア• ヤング率 0.01，ポアソン比：0.47
を基準

Uz=0 @ z=Zmin

yNx-1xJSIAM06-NOV

28

つりあい方程式∂σ x ∂τ xy ∂τ zx+++X =0∂x∂y∂z∂σ y ∂τ xy ∂τ yz+++Y = 0∂y∂x∂z∂σ z ∂τ zx ∂τ yz+++Z =0
∂z
∂x
∂y

σ x , σ y , σ z : 垂直応力
τ xy ,τ yz ,τ zx : せん断応力
X ,Y , Z :物体力zxyJSIAM06-NOV

29

変位－ひずみ関係式∂u∂v∂w
εx = , εy = , εz =
∂x
∂y∂z
∂v ∂w
∂w ∂u
∂u ∂v
γ xy =, γ zx =+ , γ yz =++∂z ∂y∂x ∂z∂y ∂xJSIAM06-NOV

30

応力－ひずみ関係式（弾性）σ x = λe + 2με x ( x j = x, y, z )jjτ x x = μγ x x ( xi , x j = x, y, z )iji je = εx +εy +εzνEEλ=, μ=(1 +ν )(1 − 2ν )
2 (1 +ν )E: ヤング率
ν: ポアソン比（0~0.5）
JSIAM06-NOV

31

ポアソン比 ννΔL

ΔL

ポアソン比＝0.50：
伸びた分だけ縮む ⇒非圧縮：解きにくい
特殊な要素を使うことがある（ν=0.4999･･･とすると非物理的挙動）
JSIAM06-NOV

32

応力－ひずみ関係式（弾性）{σ} = [D]{ε}000 ⎤⎡1 ν /(1 −ν ) ν /(1 −ν )⎢⎥1/(1−)000νν⎢⎥⎢
1000 ⎥⎢⎥
1 − 2ν
νE (1 −ν ) ⎢00 ⎥[D] =2 (1 −ν )⎥(1 +ν )(1 − 2ν ) ⎢1 − 2ν
⎢
⎥
(.)
0Sym⎢
⎥
2 (1 −ν )
⎢⎥
1−2ν
⎢
⎥
2 (1 −ν ) ⎥⎦ポアソン比＝0.50： ⎢⎣伸びた分だけ縮む ⇒非圧縮：解きにくい
特殊な要素を使うことがある（ν=0.4999･･･とすると非物理的挙動）
JSIAM06-NOV

33

計算手法，環境等• 前処理付き反復法– GPBiCG– Block SGS （Symmetric Gauss-Seidel）• 並列，ポアソン比大の場合BILU(0)，BILU(1)，BILUT等では解けない場合あり• 並列計算
– 一層のオーバーラップ（次ページ参照）– ブロックヤコビ型局所前処理• オーダリング– ベクトルプロセッサ，多段階並列化– 初期– Multicolor，Cuthilｌ-McKee• PCクラスタ– Opteron 275（2.2GHz）×64コア，PGIコンパイラ＋MPI

JSIAM06-NOV

34

Local Data Structure
Node-based Partitioning
internal nodes - elements - external nodes
PE#12122

PE#023

24

4

PE#1

517121118

13

7619

２

PE#3

３15

9４

PE#27

PE#0
2

113

1078

911101251413451010

1238911121091112963１620

148151225
1166848465512

PE#37712PE#23
JSIAM06-NOV35オーダリング手法

CM（Cuthill-McKee）：- 要素間結合グラフをもとにオーダリング- 色間の依存性も考慮

MC（Multicolor，4色）：- 要素間結合グラフをもとにオーダリング- 並列性，独立性重視

JSIAM06-NOV

36

色数と計算効率（ICCG）3D弾性問題
1283節点，「地球シミュレータ」1ノード
500

150

PDJDS/MC
125sec.

Iterations

450

400350PDJDS/MC
10075

PDJDS/CM-RCM

PDJDS/CM-RCM

30050
110100100010000

Colors

10100100010000

Colors• オーダリング手法

3025

GFLOPS

1– PDJDS/MC：マルチカラー– PDJDS/CM-RCM

20151050110100

Colors100010000• 色数が増えると反復回数は減る• FLOPS値は低下，したがって，色数が増えても計算時間が減るとは限らない。
JSIAM06-NOV

37

反復回数と色数の関係
Incompatible Nodesとは?
Doi, S. (NEC) et al.
影響の伝わり方212223242516171819201112131415678910１2345

他から影響を受けない点が「Incompatible Node」
少ない方が収束が早い
JSIAM06-NOV

38

RCMの場合1519222425

1014182123691317203581216

１24711

JSIAM06-NOV

39Red-Blackの場合並列性は高いがincompatible node多い
ILU系前処理，Gauss-Seidelで反復回数増加
11241225

13219221023619720816417518１142153

JSIAM06-NOV

40

4色の場合
並列性は弱まるがincompatible nodeは減少
ILU系前処理，Gauss-Seidelで反復回数減少
6131925723512182416224111791521310１814202

JSIAM06-NOV

41• 背景– 地震発生サイクルシミュレーション– 局所細分化，負荷分散• 問題設定• 計算結果• まとめ– 今後の予定等

JSIAM06-NOV

予備計算：1領域• メッシュ1：なだらかな細分化– 203から初めて1要素を3段階に分割– 約68,600節点，67,000要素– ヤング率：1.00，ポアソン比：0.24• メッシュ2：均等– メッシュ数がほぼ「メッシュ1」と同じ均等な403メッシュ– 約69,000節点，64,000要素– ヤング率：1.00，ポアソン比：0.24，0.47，0.495• メッシュ3：急激な細分化– メッシュ2と同じ規模だが，要素長さを急激に変化• 隣接した部分で：100：3（ほぼ3段階の細分化に相当）• x，y，z方向に5層ずつ– ヤング率：1.00，ポアソン比：0.24
JSIAM06-NOV42予備計算：反復回数400●■▲ メッシュ1（ν=0.24）○□△ メッシュ2（ν=0.24）●■▲ メッシュ2（ν=0.47）●■▲ メッシュ3（ν=0.24）

ITERATIONS

300200○ Multicolor オーダリング□ Initial オーダリング△ Cuthill-McKee

1000
1.E+01

1.E+02

1.E+03

1.E+04COLOR#

メッシュ1と2（均質）の差はほとんど無い
JSIAM06-NOV

メッシュ1

メッシュ3
43

予備計算：反復回数
ポアソン比の影響
700○□△ メッシュ2（ν=0.24）●■▲ メッシュ2（ν=0.47）●■▲ メッシュ2（ν=0.495）

600

ITERATIONS

500400300○ Multicolor オーダリング□ Initial オーダリング△ Cuthill-McKee

200
1000
1.E+01

1.E+02

1.E+03

1.E+04COLOR#

ポアソン比が0.50に近づくと急激に収束悪化
色数増加による収束の向上顕著：incompatible nodes
ただ，解けても物理的にはあまり意味がない
JSIAM06-NOV44

64領域の計算：負荷分散の効果：節点数
ParMETISの係数の影響，アセノスフェアのν=0.47△▲ 負荷分散前の最大，最小○● 負荷分散後の最大，最小

8.E+04NODE#

6.E+04係数が1.0に近いほど厳密な負荷分散を適用する

4.E+04

本ケースの場合，節点数がバランスするように設定してある

2.E+040.E+00
1.00

1.10

1.20

1.30

1.40

1.50

1.60

FACTOR
JSIAM06-NOV45

64領域の計算：負荷分散の効果：要素数
ParMETISの係数の影響，アセノスフェアのν=0.47△▲ 負荷分散前の最大，最小○● 負荷分散後の最大，最小8.E+04

Element#

6.E+04係数が1.0に近いほど厳密な負荷分散を適用する

4.E+04

2.E+040.E+00
1.00

1.10

1.20

1.30

1.40

1.50

1.60

FACTOR
JSIAM06-NOV

46

64領域の計算：負荷分散の効果：
非ゼロ成分数，edge-cut
ParMETISの係数の影響，アセノスフェアのν=0.47
非ゼロ成分数△▲ 負荷分散前の最大，最小○● 負荷分散後の最大，最小

Off-Diagonal#, edge-cut#

2.00E+06

1.50E+06

edge-cut▲ 負荷分散前● 負荷分散後

1.00E+06

5.00E+05

0.00E+00
1.00

1.10

1.20

1.30

1.40

1.50

負荷バランスが悪いほど，edge-cut
1.60 数は少ない：並列計算の収束は良い

FACTOR
JSIAM06-NOV47

64領域の計算：負荷分散の効果：
非ゼロ成分数，edge-cut
ParMETISの係数の影響，アセノスフェアのν=0.47
非ゼロ成分数△▲ 負荷分散前の最大，最小○● 負荷分散後の最大，最小

Off-Diagonal#, edge-cut#

2.00E+06

1.50E+06

edge-cut▲ 負荷分散前● 負荷分散後

1.00E+06

5.00E+05

0.00E+00
1.00

1.10

1.20

1.30

FACTOR
JSIAM06-NOV

1.40

1.50

負荷バランスが悪いほど，非ゼロ
1.60 成分の最大値が増加⇒全体として計算時間がかかる
48

EDGE-CUTとは ?
z 辺の両端の節点（または要素）が異なった領域に属している場合，「EDGE-CUTが生じている。」という。
z EDGE-CUTが少ないほど，通信は少ない。

EDGE-CUT無し

EDGE-CUT有り

JSIAM06-NOV

49

64領域の計算：負荷分散の効果：
反復回数，計算時間
ParMETISの係数の影響，アセノスフェアのν=0.47
反復回数▲ 負荷分散前● 負荷分散後500

ITERATIONS, sec.

400

計算時間（sec.）△ 負荷分散前○ 負荷分散後300

200

10001.00

1.10

1.20

1.30

FACTOR
JSIAM06-NOV

1.40

1.50

edge-cut数が少ないほど並列計算
1.60 の収束は良い：計算時間は負荷分散している状態でほぼ同じ
50

64領域の計算：負荷分散の効果：反復回数
係数=1.20，アセノスフェアポアソン比と色数の影響
ポアソン比：0.470● Multicolor■ 初期オーダリング▲ Cuthill-McKee

800

ITERATIONS

600

ポアソン比：0.495○ Multicolor□ 初期オーダリング△ Cuthill-McKee

400

200

0
1.E+01

1.E+02

1.E+03

COLOR#JSIAM06-NOV

1.E+04Multicolorにおける色数と反復回数の関係は一様ではない
51

64領域の計算：負荷分散の効果：計算時間係数=1.20，アセノスフェアポアソン比と色数の影響
ポアソン比：0.470● Multicolor■ 初期オーダリング▲ Cuthill-McKee

150100sec.

ポアソン比：0.495○ Multicolor□ 初期オーダリング△ Cuthill-McKee

50

01.E+01

1.E+02

1.E+03

1.E+04COLOR#JSIAM06-NOV

52

反復回数が減っているのに計算時間が増加している ?
150

800

100sec.

ITERATIONS

600

400

50
200

01.E+01

01.E+02

1.E+03

COLOR#JSIAM06-NOV

1.E+04

1.E+01

1.E+02

1.E+03

1.E+04COLOR#
53

1反復あたりの計算時間を測定してみるとポアソン比：0.470● Multicolor■ 初期オーダリング▲ Cuthill-McKee

0.25

sec./iteration

0.20

ポアソン比：0.495○ Multicolor□ 初期オーダリング△ Cuthill-McKee

0.150.100.05

色数が少ないときの性能が悪い

0.00
1.E+01

1.E+02

1.E+03

1.E+04COLOR#JSIAM06-NOV54

1領域の計算について測定
前処理部分（全体の50%以上）1反復あたり計算時間
メッシュ1：細分化メッシュ● Multicolor■ 初期オーダリング▲ Cuthill-McKee

0.30

MINV / iteration (sec.)

0.25

0.20
0.15

メッシュ2：均質メッシュ○ Multicolor□ 初期オーダリング△ Cuthill-McKee

0.10

メッシュ1がより顕著
0.05

0.00
1.E+01

1.E+02

1.E+03

1.E+04COLOR#JSIAM06-NOV55

1領域の計算について測定
上三角成分または下三角成分が無い行の総数
メッシュ1：細分化メッシュ● Multicolor■ 初期オーダリング▲ Cuthill-McKee

2.0E+04

Node #

1.5E+04

1.0E+04

メッシュ2：均質メッシュ○ Multicolor□ 初期オーダリング△ Cuthill-McKee

5.0E+03

0.0E+00
1.E+01

1.E+02

1.E+03

1.E+04COLOR#JSIAM06-NOV

56

1領域の計算について測定
上三角成分または下三角成分が無い行の総数
上三角成分，または下三角成分のいずれかしか無い行では，最も内側のループの計算量が増加しメモリに負担がかかる。
メッシュ1は非対角成分の数がもともと多いため，この傾向が顕著になると思われる。
⇒ESでは別の挙動 ?

2.0E+04

Node #

1.5E+04

1.0E+04

5.0E+03

0.0E+00
1.E+01

1.E+02

1.E+03

COLOR#JSIAM06-NOV

1.E+04

do i= 1, N
do j= 1, INL(i)Y(i)= Y(i) + AL(j,i)*X(IAL(j,i))
enddo
do j= 1, INU(i)Y(i)= Y(i) + AU(j,i)*X(IAU(j,i))
enddo
enddo
57

58

• 背景– 地震発生サイクルシミュレーション– 局所細分化，負荷分散• 問題設定• 計算結果• まとめ– 今後の予定等

JSIAM06-NOV

59

まとめ
• 並列FEMによる「すべり応答関数」計算– AMR– 並列27分木– ParMETISによる負荷分散＋データマイグレーション• 大規模（でもないが）弾性問題– ポアソン比増加による収束性能悪化はあるが顕著では無い– 局所Block SGSで十分解ける• メッシュ幅変化がなだらかなため• 負荷分散と収束性のトレードオフ• 色数と「性能」の関係
JSIAM06-NOV

60

今後の予定
• 並列27分木の完全実装• 「地球シミュレータ」への移植，実計算• ポアソン比の扱い• 「高級な」前処理の適用– 不均質な場合– 領域間オーバーラップ

JSIAM06-NOV