Microsoft PowerPoint - KN-2006NOV16.ppt

Size: px

Start display at page:

Download "Microsoft PowerPoint - KN-2006NOV16.ppt"

こうしろうごみぶち
5 years ago
Views:

1 局所細分化メッシュに基づく並列有限要素法における前処理付き反復法 Preconditioned Iterative Methods for Parallel Finite-Element Applications with Adaptive Mesh Refinement 中島研吾 (1) 兵藤守 (2) (1) 東京大学大学院理学系研究科地球惑星科学専攻 (2) 地球シミュレータセンター固体地球シミュレーション研究グループ日本応用数理学会行列固有値問題の解法とその応用研究部会第 2 回研究会 2006 年 11 月 16 日, 国立情報学研究所

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

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

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

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

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

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

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

9 Refinement Patterns of Tetrahedron 1 edge marked 3 edges on a same surface are marked 2 edges on a same surface are marked more than 2 edges on different surfaces are marked 2 children (Binary) 3rd edge is automatically marked 4 children (Quadtree) more than 4 edges are marked 6 edges are automatically marked 8 children (Octree) 9

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

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

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

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

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

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

16 16 27 分木による細分化例細分化レベル 3( 拡大図 )

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

18 18 並列 27 分木要素数, 節点数が 100 万を越えるような初期メッシュを局所的に細分化する必要がある細分化レベルを増やしていくと,1 プロセッサで扱うことは困難色々な細分化パターンを実施する必要があるすべり応答関数を算出するのが目的であるため, 多数のメッシュパターンを高速に生成する必要がある並列計算初期メッシュ, 領域分割並列細分化動的負荷分散 (Dynamic Load Balancing)+データマイグレーションマントル四国フィリピン海プレート z y x 伊豆太平洋プレートマントル

19 Supersonic Flow around a Sphere M=1.40 Uniform Flow, Ideal Gas, Re= PE cases (tetrahedron only) Initial Grid 1-Lev. Adapted 2-Lev. Adapted before DRAMA PE PE PE PE before DRAMA before DRAMA

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

21 METIS 21 マルチレベルグラフ理論に基づいた方法

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

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

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

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

26 26 例題 (2/3) 細分化メッシュ 4 領域に 1 領域の割合で,1 要素を 3 段階に細分化 70,000メッシュ近い領域が出てくる全体として, 約 150 万要素 ( 節点 ),450 万自由度 ( 正確には 1,480,960 節点,4,442,880 自由度 ) ParMETISによる動的負荷分散 +データマイグレーションリソスフェア ( 地殻 ) 弾性体アセノスフェア ( マントル ) 粘弾性体

27 27 境界条件右図参照現実的なもの ( 食い違い変位 ) ではない例題 (3/3) z Uniform Distributed Force in z=z min 物性値リソスフェアヤング率 1.00, ポアソン比 :0.24 アセノスフェアヤング率 0.01, ポアソン比 :0.47 を基準 U y y=y min N z -1 N y -1 U z z=z min x U x x=x min y N x -1

28 28 つりあい方程式 = = = Z y x z Y z x y X z y x yz zx z yz xy y zx xy x τ τ σ τ τ σ τ τ σ :,, :,, zx yz xy z y x τ τ τ σ σ σ 垂直応力せん断応力 :,, Z Y X 物体力 x z y

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

30 30 応力 - ひずみ関係式 ( 弾性 ) ) 2(1, ) 2 )(1 (1 ),,, ( ),, ( 2 ν μ ν ν ν λ ε ε ε μγ τ με λ σ + = + = + + = = = = + = E E e z y x x x z y x x e z y x j i x x x x j x x j i j i j j E: ヤング率 ν: ポアソン比 (0~0.5)

31 31 ポアソン比 ν νδl ΔL ポアソン比 =0.50: 伸びた分だけ縮む非圧縮 : 解きにくい特殊な要素を使うことがある (ν= とすると非物理的挙動)

32 32 応力 - ひずみ関係式 ( 弾性 ) {} [ ]{} [ ] + = = ) 2( ) 2(1 2 1.) ( 0 0 ) 2( ) /( ) /(1 ) /(1 1 ) 2 )(1 (1 ) (1 ν ν ν ν ν ν ν ν ν ν ν ν ν ν ν ν Sym E D ε D σ ポアソン比 =0.50: 伸びた分だけ縮む非圧縮 : 解きにくい特殊な要素を使うことがある (ν= とすると非物理的挙動 )

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

34 Local Data Structure Node-based Partitioning internal nodes - elements - external nodes PE#1 PE#0 PE# PE# PE#3 PE#2 1 2 PE# PE#2

35 35 オーダリング手法 CM(Cuthill-McKee): - 要素間結合グラフをもとにオーダリング - 色間の依存性も考慮 MC(Multicolor,4 色 ): - 要素間結合グラフをもとにオーダリング - 並列性, 独立性重視

36 36 ICCG 3D PDJDS/MC 125 Iterations 400 sec. 100 PDJDS/MC GFLOPS PDJDS/CM-RCM Colors Colors Colors PDJDS/CM-RCM オーダリング手法 PDJDS/MC: マルチカラー PDJDS/CM-RCM 色数が増えると反復回数は減る FLOPS 値は低下, したがって, 色数が増えても計算時間が減るとは限らない

37 反復回数と色数の関係 37 Incompatible Nodes とは? Doi, S. (NEC) et al 影響の伝わり方他から影響を受けない点が Incompatible Node 少ない方が収束が早い

38 38 RCM の場合

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

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

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

42 メッシュ 1: なだらかな細分化予備計算 :1 領域 20 3 から初めて1 要素を3 段階に分割約 68,600 節点,67,000 要素ヤング率 :1.00, ポアソン比 :0.24 メッシュ 2: 均等メッシュ数がほぼメッシュ1 と同じ均等な40 3 メッシュ約 69,000 節点,64,000 要素ヤング率 :1.00, ポアソン比 :0.24,0.47,0.495 メッシュ 3: 急激な細分化メッシュ 2 と同じ規模だが, 要素長さを急激に変化隣接した部分で :100:3( ほぼ3 段階の細分化に相当 ) x,y,z 方向に5 層ずつヤング率 :1.00, ポアソン比 :

予備計算 : 反復回数 ITERATIONS 400 300 200 100 0 1.E+01 1.E+02 1.E+03 1.E+04 COLOR# メッシュ 1(ν=0.24) メッシュ 2(ν=0.

43 予備計算 : 反復回数 ITERATIONS E+01 1.E+02 1.E+03 1.E+04 COLOR# メッシュ 1(ν=0.24) メッシュ 2(ν=0.24) メッシュ 2(ν=0.47) メッシュ 3(ν=0.24) Multicolor オーダリング Initial オーダリング Cuthill-McKee メッシュ 1 と 2( 均質 ) の差はほとんど無いメッシュ 1 メッシュ 3 43

44 予備計算 : 反復回数ポアソン比の影響 ITERATIONS E+01 1.E+02 1.E+03 1.E+04 COLOR# メッシュ 2(ν=0.24) メッシュ 2(ν=0.47) メッシュ 2(ν=0.495) Multicolor オーダリング Initial オーダリング Cuthill-McKee ポアソン比が0.50に近づくと急激に収束悪化色数増加による収束の向上顕著 :incompatible nodes ただ, 解けても物理的にはあまり意味がない 44

45 64 領域の計算 : 負荷分散の効果 : 節点数 ParMETIS の係数の影響, アセノスフェアの ν=0.47 NODE# 8.E+04 6.E+04 4.E+04 2.E+04 負荷分散前の最大, 最小負荷分散後の最大, 最小係数が 1.0 に近いほど厳密な負荷分散を適用する本ケースの場合, 節点数がバランスするように設定してある 0.E FACTOR 45

46 64 領域の計算 : 負荷分散の効果 : 要素数 ParMETIS の係数の影響, アセノスフェアの ν= E+04 負荷分散前の最大, 最小負荷分散後の最大, 最小 Element# 6.E+04 4.E+04 係数が 1.0 に近いほど厳密な負荷分散を適用する 2.E+04 0.E FACTOR 46

47 64 領域の計算 : 負荷分散の効果 : 非ゼロ成分数,edge-cut ParMETIS の係数の影響, アセノスフェアの ν=0.47 Off-Diagonal#, edge-cut# 2.00E E E E E FACTOR 非ゼロ成分数負荷分散前の最大, 最小負荷分散後の最大, 最小 edge-cut 負荷分散前負荷分散後負荷バランスが悪いほど,edge-cut 数は少ない : 並列計算の収束は良い 47

48 64 領域の計算 : 負荷分散の効果 : 非ゼロ成分数,edge-cut ParMETIS の係数の影響, アセノスフェアの ν=0.47 Off-Diagonal#, edge-cut# 2.00E E E E E FACTOR 非ゼロ成分数負荷分散前の最大, 最小負荷分散後の最大, 最小 edge-cut 負荷分散前負荷分散後負荷バランスが悪いほど, 非ゼロ成分の最大値が増加全体として計算時間がかかる 48

49 EDGE-CUT とは? 辺の両端の節点 ( または要素 ) が異なった領域に属している場合, EDGE-CUTが生じているという EDGE-CUTが少ないほど, 通信は少ない EDGE-CUT 無し EDGE-CUT 有り 49

50 64 領域の計算 : 負荷分散の効果 : 反復回数, 計算時間 ParMETIS の係数の影響, アセノスフェアの ν=0.47 ITERATIONS, sec 反復回数負荷分散前負荷分散後計算時間 (sec.) 負荷分散前負荷分散後 FACTOR edge-cut 数が少ないほど並列計算の収束は良い : 計算時間は負荷分散している状態でほぼ同じ 50

51 64 領域の計算 : 負荷分散の効果 : 反復回数係数 =1.20, アセノスフェアポアソン比と色数の影響 ITERATIONS ポアソン比 :0.470 Multicolor 初期オーダリング Cuthill-McKee ポアソン比 :0.495 Multicolor 初期オーダリング Cuthill-McKee 0 1.E+01 1.E+02 1.E+03 1.E+04 COLOR# Multicolor における色数と反復回数の関係は一様ではない 51

52 64 領域の計算 : 負荷分散の効果 : 計算時間係数 =1.20, アセノスフェアポアソン比と色数の影響 sec ポアソン比 :0.470 Multicolor 初期オーダリング Cuthill-McKee ポアソン比 :0.495 Multicolor 初期オーダリング Cuthill-McKee 0 1.E+01 1.E+02 1.E+03 1.E+04 COLOR# 52

53 反復回数が減っているのに計算時間が増加している? ITERATIONS 400 sec E+01 1.E+02 1.E+03 1.E+04 1.E+01 1.E+02 1.E+03 1.E+04 COLOR# COLOR# 53

54 1 反復あたりの計算時間を測定してみるとポアソン比 :0.470 Multicolor 初期オーダリング Cuthill-McKee sec./iteration ポアソン比 :0.495 Multicolor 初期オーダリング Cuthill-McKee E+01 1.E+02 1.E+03 1.E+04 COLOR# 色数が少ないときの性能が悪い 54

55 1 領域の計算について測定前処理部分 ( 全体の 50% 以上 )1 反復あたり計算時間 MINV / iteration (sec.) メッシュ 1 がより顕著 1.E+01 1.E+02 1.E+03 1.E+04 COLOR# メッシュ 1: 細分化メッシュ Multicolor 初期オーダリング Cuthill-McKee メッシュ 2: 均質メッシュ Multicolor 初期オーダリング Cuthill-McKee 55

1 領域の計算について測定上三角成分または下三角成分が無い行の総数 2.0E+04 1.

Cuthill-McKee Node # 1.0E+04 5.0E+03 0.0E+00 1.E+01 1.

56 1 領域の計算について測定上三角成分または下三角成分が無い行の総数 2.0E E+04 メッシュ 1: 細分化メッシュ Multicolor 初期オーダリング Cuthill-McKee Node # 1.0E E E+00 1.E+01 1.E+02 1.E+03 1.E+04 COLOR# メッシュ 2: 均質メッシュ Multicolor 初期オーダリング Cuthill-McKee 56

57 1 領域の計算について測定上三角成分または下三角成分が無い行の総数 Node # 2.0E E E+04 上三角成分, または下三角成分のいずれかしか無い行では, 最も内側のループの計算量が増加しメモリに負担がかかるメッシュ 1 は非対角成分の数がもともと多いため, この傾向が顕著になると思われる ES では別の挙動? 5.0E E+00 1.E+01 1.E+02 1.E+03 1.E+04 COLOR# 57

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

59 59 まとめ FEM AMR 27 ParMETIS Block SGS

60 60 今後の予定 27

Microsoft PowerPoint - AMR.ppt [互換モード]

Microsoft PowerPoint - AMR.ppt [互換モード] 六面体メッシュの適合型並列局所細分化と負荷分散 Conforming parallel adaptive mesh refinement for hexahedral elements 中島研吾東京大学情報基盤センター 2008 年並列 / 分散 / 協調処理に関する佐賀サマーワークショップ (SWoPP 佐賀 2008) 2008 年 8 月 5 日 ~7 日 SWoPP2008 2