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

(1)

六面体メッシュの

適合

並

細

適合

並

細

適合型並列局所細分化と

負荷分散

Conforming parallel adaptive mesh refinement for

hexahedral elements

中島研吾

東京大学情報基盤センター 2008年並列／分散／協調処理に関する『佐賀』サマー・ワークショップ（SWoPP佐賀2008） 2008年8月5日～7日

(2)

SWoPP2008 2

概要

• AMRとData Migration

• 断層における応力蓄積シミュレーション

• 単位すべり応答関数計算用フレームワーク

– 六面体メッシュと局所細分化 – 負荷分散

• 評価

評価

– PCクラスタ – T2Kオープンスパコン（東大）T2Kオプンスン（東大）

• まとめ

(3)

SWoPP2008 SWoPP2008 33

適応格子法（

適応格子法（AMR

AMR）（

）（1/2

1/2）

）

適応格子法（

適応格子法（AMR

AMR）（

）（1/2

1/2）

）

• 差分法有限要素法などで空間を離散化する場合

差分法，有限要素法などで空間を離散化する場合，

精度良い解を求めるためには，細かいメッシュが必

要である

要である。

• 細かいメッシュを切ると，それだけ，計算量，必要記

憶容量が増加する。

• 細かいメッシュを必要とするのは解析領域全体で

• 細かいメッシュを必要とするのは，解析領域全体で

はなく，ごく一部である。

衝撃波渦剥離

– 衝撃波，渦，剥離

– 応力集中，亀裂

– 従属変数が急激に変化する領域

(4)

適応格子法（

適応格子法（AMR

AMR）（

）（2/2

2/2）

）

適応格子法（

適応格子法（AMR

AMR）（

）（2/2

2/2）

）

• しかしながら多くの場合これらの現象が生じる領

しかしながら，多くの場合，これらの現象が生じる領

域は，前もって分かっていない。

比較的粗いメシから始めて計算結果に応じて

• 比較的粗いメッシュから始めて，計算結果に応じて，

局所的な細分化（refinement）を自動的に実施し，

精度良い解を効率的に求める手法が，適応格子法

（Adaptive Mesh Refinement，AMR）である。

（ dapt e

es

e

e e t，

）である。

– 非定常問題の場合，細かいメッシュが必要でなくなった

領域のメッシュを粗くする（coarsening）場合もありうる

領域のメッシュを粗くする（coarsening）場合もありうる。

(5)

適応格子法の応用例（

適応格子法の応用例（1/4

1/4）

）

適応格子法の応用例（

適応格子法の応用例（1/4

1/4）

）

構造力学，有限要素法

Mises Stress

Initial Mesh Mises Stress Adapted Mesh

(6)

適応格子法の応用例（

適応格子法の応用例（2/4

2/4）

）

適応格子法の応用例（

適応格子法の応用例（2/4

2/4）

）

流体力学，有限要素法

2-Lev. Adapted 1-Lev. Adapted Initial Grid

(7)

適応格子法の応用例（

適応格子法の応用例（3/4

3/4）

）

適応格子法の応用例（

適応格子法の応用例（3/4

3/4）

）

流体力学，ブロック格子

• Rayleigh-Taylor Instability

• Block-structured adaptive

grid

• 宇宙物理学の世界ではよく使

われているようである。

(8)

適応格子法の応用例（

適応格子法の応用例（4/4

4/4）

）

適応格子法の応用例（

適応格子法の応用例（4/4

4/4）

）

可視化，ボリュームレンダリング

L Chen (Tsinghua U ) H Matsui (UC Berkeley)

• FEMによる地球外核のMHD解析結果

非構造格子の計算結果をブロク格子にマピングしてレンダ

L.Chen (Tsinghua U.), H.Matsui (UC Berkeley)

• 非構造格子の計算結果をブロック格子にマッピングしてレンダ

(9)

適応格子法の応用例（

適応格子法の応用例（4/4

4/4）

）

適応格子法の応用例（

適応格子法の応用例（4/4

4/4）

）

可視化，ボリュームレンダリング（続き）

かなり複雑な

状も適

能ある（

L.Chen (Tsinghua U.)

• かなり複雑な形状にも適用可能である（Mobile Pentium

(10)

適応格子法の現状・問題点

• 並列計算

局所的な細分化によって負荷バランスが均等でなくな

– 局所的な細分化によって，負荷バランスが均等でなくな

る可能性がある ⇒ 動的負荷分散（Dynamic Load

Balancing：DLB）

• 形状再現性

– 例えば，曲面上のメッシュを細分化する場合，もとの形

状を忠実に再現する必要がある。

– CADとの連動。

– 近似的な曲面表現手法の開発

近似的な曲面表現手法の開発。

(11)

Supersonic Flow around a Sphere

M 1 40 Uniform Flow Ideal Gas Re 10

66

M=1.40 Uniform Flow, Ideal Gas, Re=10

66

4 PE cases (tetrahedron only)

2-Lev. Adapted 1-Lev. Adapted Initial Grid before DRAMA 3834 2527 2769 2526 before DRAMA 793 652 696 650 before DRAMA PE0 137 PE1 137 -SWoPP2008 11 2703 2522 1390 2524 668 652 448 651 PE2 136 PE3 136

(12)

-SWoPP2008 SWoPP2008 1212

負荷分散のためのプロセス

AMR

• AMR

• 領域再分割（Repartitioing）

p

g

– ParMETIS（U.Minnessota），JOSTLE（U.Greenwich）

• デタ再構成デタマイグレション（Data

• データ再構成，データマイグレーション（Data

Migration）

– Zoltan（Sandia）

(13)

SWoPP2008 13

Procedures for Parallel AMR

• Start from initial distributed mesh

(14)

SWoPP2008 14

Procedures for Parallel AMR

• Start from initial distributed mesh

(15)

SWoPP2008 15

Procedures for Parallel AMR

• Start from initial distributed mesh

(16)

SWoPP2008 16

Procedures for Parallel AMR

• Parallel AMR

– zones to be refined are specified – based on templates

Load Imbalance !!

(17)

SWoPP2008 17

Procedures for Parallel AMR

• Repartitioning for Load Balancing

(18)

SWoPP2008 18

Procedures for Parallel AMR

• Data Migration

(19)

データマイグレーション

デタイグ

ション

デタイグ

ション

• 領域再分割（Repartitioning）で決まる

領

p

g

決

のは，「節点」の所属先（宛先）

自由度が存在する

– 自由度が存在する

• 以下については別途移動する

– 節点の座標，境界条件

– 接続する要素，その上の節点

接続する要素，その上の節点

– 要素物性等

領域間通信テブル

– 領域間通信テーブル

• アプリケーションに依存：標準化困難?

• 結構時間がかかる

(20)

SWoPP2008 20

Local Data Structure

Node-based Partitioning

internal nodes - elements - external nodes

21 22 23 24 25 6 7 15 PE#0 4 5PE#16 12 21 22 23 24 25 PE#0 PE#1 21 22 23 24 25 16 17 18 19 20 4 5 14 13 PE#0 3 11 1 2 16 17 18 19 20 16 17 18 19 20 16 20 11 12 13 14 15 1 2 3 10 7 8 9 10 10 12 11 16 20 11 12 13 14 15 16 20 11 12 13 14 155 6 7 8 9 10 8 9 11 12 10 9 11 12 6 9 5 5 6 7 8 9 10 5 6 7 8 9 10 １２３４ 5 5 6 8 4 3 4 8 6 １２３４ 5 １２３４ 5 １２３４ 5 7 1 2 3 PE#2 1 2 7 PE#3 １２３４ 5 PE#2 PE#3 １２３４ 5

(21)

SWoPP2008 21

Adaptive Mesh Refinement (AMR) and

Load Balancing (LB)

Some slides from:

N k ji K Fi b J Ok d H (2001) P ll l 3D Ad i Nakajima, K., Fingberg, J. Okuda, H. (2001), Parallel 3D Adaptive

Compressible Navier-Stokes Solver in GeoFEM with Dynamic Load-Balancing by DRAMA Library Lecture Notes in Computer Science Balancing by DRAMA Library, Lecture Notes in Computer Science 2110, 183-193.

*DRAMA Library: Repartitioning Library developed by NEC Europe

DR

AM

A

DR

AM

A

DR

A

DR

(22)

SWoPP2008 22

Supersonic Flow around a Sphere

M=1 40 Uniform Flow Ideal Gas Re=10

6

M=1.40 Uniform Flow, Ideal Gas, Re=10

6

4 PE cases (tetrahedron only)

2-Lev. Adapted 1-Lev. Adapted Initial Grid before DRAMA 3834 2527 2769 2526 before DRAMA 793 652 696 650 before DRAMA PE0 137 PE1 137 -2703 2522 1390 2524 668 652 448 651 PE2 136 PE3 136

(23)

Refinement Patterns of Tetrahedron

NO “H

i

” N d

NO “H

i

” N d

NO “Hanging” Nodes

3 edges on a same 1 edge marked a same surface are marked 2 edges on a same 2 children (Binary) a same surface are marked 4 children 3rd edge is _{4 children} (Quadtree) more than 2 edges on different surfaces are marked 3 d edge s automatically marked marked more than SWoPP2008 23 8 children (Octree) more than 4 edges are

marked _{automatically}6 edges are marked

(24)

Directional Refinement of Prisms

SWoPP2008 24

Initial Grid

Initial Grid Edge CutEdge Cut Embedded grid obtained after quadtreeEmbedded grid obtained after quadtreeand binary divisionsand binary divisions of the triangular faces

(25)

SWoPP2008 25

Results (Hybrid) : LAMP Cluster

8 PE cases with 1 level adapted meshes

8 PE cases with 1-level adapted meshes

Initial Mesh : 16,050 nodes, 76,800 meshes

Final Mesh : 29,852 nodes, 117,694 meshes

P

P--JOSTLE : Vicinity of the Sphere Wall SurfaceJOSTLE : Vicinity of the Sphere Wall Surface

(26)

SWoPP2008 26

Results (Tet.Only) : LAMP Cluster

8 PE cases with 2-level adapted meshes

Final Mesh : 10,240 nodes, 69,462 tet.

P

P--JOSTLEJOSTLE PP--MMEETTIISS k k--wayway NO NO Repartition Repartition RCB RCB bucket bucket

(27)

SWoPP2008 27

Results (Tet.Only) : LAMP Cluster

16 PE cases with 1 level adapted meshes

16 PE cases with 1-level adapted meshes

Initial Mesh : 16,050 nodes, 92,160 tet.

Final Mesh : 47,074 nodes, 306,236 tet.

Elaplsed Time for 1,000 steps

Repartitiong

Repartitiong Edge #Edge # Repartitiong

Repartitiong Methods

Methods (min./max.)(min./max.)Edge #Edge # Total EdgecutTotal Edgecut Time (sec.)Time (sec.)

NO Repartition NO Repartition _{10,576/48,495}_{10,576/48,495} _39,888_39,888 _1,683_1,683 NO Repartition NO Repartition P P--JOSTLEJOSTLE 10,576/48,495 10,576/48,495 21,089/22,233 21,089/22,233 39,888 39,888 25,085 25,085 1,683 1,683 874 874 Para M Para MEETTIISS 21,201/22,63021,201/22,630 26,27426,274 880 880 RCB simple RCB simple _{22,520/23,090}_{22,520/23,090} _41,980_41,980 ₈₉₉₈₉₉ RCB bucket RCB bucket _{21,231/23,269}_{21,231/23,269} _37,192_37,192 ₉₂₆₉₂₆

(28)

SWoPP2008 28 Results Results p A D epHYBRID CNTL. A D A P R A M P T A Grid Grid

(29)

SWoPP2008 29 Results Results p A D epHYBRID CNTL. A D A P R A M P T A Grid Grid

AM

A

AM

A

DR

AM

DR

AM

(30)

SWoPP2008 30

GeoFEM: FY.1998-2002

http://geofem.tokyo.rist.or.jp/

• 文部科学省「科学技術振興調整費総合研究」

– 「高精度の地球変動予測のための並列ソフトウェア開

発に関する研究」の一部

• 固体地球シミュレーション用並列有限要素法プ

ラットフォーム

ラットフォ

ム

– 並列I/O，並列線形ソルバー，並列可視化をサポート

固体地球地球内部⇒コアマントル対流地殻変動

– 固体地球：地球内部⇒コア・マントル対流，地殻変動

• HPCと自然科学の緊密な協力

(31)

心残り

DRAMA（Repartitioning）のあとのData Migration

• DRAMA（Repartitioning）のあとのData Migration

(32)

本研究の目的

以下の機能を有する並列メッシ生成フレムワク

• 以下の機能を有する並列メッシュ生成フレームワーク

– 並列細分化（AMR）

– Repartitioning

– Data Migration

Data Migration

• Data Migrationを「自作」

ブ

プ

– 標準化ライブラリのためのプロトタイプ

– とにかく自分で実装してみる

(33)

SWoPP2008 33

• AMRとData Migration

• 断層における応力蓄積シミュレーション

• 単位すべり応答関数計算用フレームワーク

• 評価

評価

• まとめ

(34)

34

地震発生サイクルシミュレーション

• プレート境界における準静的応力蓄積過程

もも

プ

• もともとのアプローチ

– 非線形接触問題をNewton-Raphson法によって解く – ALM法（Augmented Lagrangean,拡大ラグランジェ法）による拘束条件：ペナルティ数 – 領域分割による並列有限要素法 • GeoFEM，地球シミュレータ前処理付き反復法 [N k ji 2003] [N k ji 2006] • 前処理付き反復法 [Nakajima, 2003]，[Nakajima, 2006] SWoPP2008

(35)

35

現状

• 2003年以降全てをFEMで解こうという考え方から方針転換

現状

• 2003年以降，全てをFEMで解こうという考え方から方針転換

– 接触モデリングの困難さ境界条件 – 境界条件

伝統的な境

積分法境

素法を使する

• 伝統的な境界積分法，境界要素法を使用する

– 単位すべり応答関数（ある場所，t=0における単位すべり量が時間場所おる応力変を与るる寄与積分重ね合わ t-s, 場所xにおける応力変化を与える）による寄与の積分（重ね合わせ）不均質場における「単位すべり応答関数を計算する必要がある – 不均質場における「単位すべり応答関数」を計算する必要がある • 有限要素法の使用 • 〔Hori et al 2004〕 SWoPP2008 • 〔Hori et al 2004〕

(36)

SWoPP2008 36

BEM/BIEM Approach (1/3)

pp

(

)

[Hashimoto & Matsu’ura, 2000]

Elastic lithosphere divided into 2 parts

i fi it l l ti l i t f i fi it l l ti l i t f

infinitely long vertical interface infinitely long vertical interface

interaction between the 2 parts : fault slip interaction between the 2 parts : fault slip

Fault Slip

p

“

w

”

v

vplpl steady plate motion ratesteady plate motion rate

u u perturbationperturbation Plate A

)

t

,

(

u

v

)

t

,

(

w

x



_pl



x

Plate B

(37)

SWoPP2008 37

BEM/BIEM Approach (2/3)

pp

(

)

[Hashimoto & Matsu’ura, 2000]

Shear Stress (



) due to Fault Slip(w)

H

H Unit “Slip Response Function”, Unit “Slip Response Function”, analyticallypp pp analyticallyyy yy calc

calcuulatedlated BEFORE simulation)BEFORE simulation) 1st term : steady plate motion

1st term : steady plate motion

2 d t li t b ti

2nd term: slip perturbation 2nd term: slip perturbation





 

u

(

ξ

,

)

_H

₍

ξ

₀

₎

_d

ξd

)

(

)

(

t steady plate

motion slip perturbation





















₍

_,

_t

₎

₍

_,

_t

₎

 

u

(

ξ

,

)

_H

₍

_,

_t

_;

ξ

_,

₀

₎

_d

ξd

0 s 0

x

Constitutive Relation between

“

w

”

and

“

u

”

Aochi and Matsu

‘

ura (1999) : New Modelura (1999) : New Model



x



x

,

t

)

f

w

(

,

t

);

(





Aochi and Matsu

(38)

SWoPP2008 38

BEM/BIEM Approach (3/3)

pp

(

)

[Hashimoto & Matsu’ura, 2000]

Spatial Discretization

p

slip

“

uu

”

is defined by combination of is defined by combination of

“

MM

”

Cubic Spline Cubic Spline Basis. Basis.

)

(

)

t

(

)

t

(

M





a

(

t

)

(

)

t

,

(

u

_m 1 m m

x









Solve Equations

Construct Nonlinear Coupled Equation Construct Nonlinear Coupled Equationpp qq Linealized by Levenberg

Linealized by Levenberg--Marquardt Method using Marquardt Method using

“

NN

”

points points N M N M N>M N>M

solve linearized equation for

“

aamm(t)(t)

”

at each iterationat each iteration

requires M*M dense matrix inversion requires M*M dense matrix inversion

(39)

SWoPP2008 39

Unit “Slip Response Function”: H

p





















₍

_,

_t

₎

₍

_,

_t

₎

 

u

(

ξ

,

)

_H

₍

_,

_t

_;

ξ

_,

₀

₎

_d

ξd

t 0

x

steady plate

motion







slip perturbation

 

0 s

• “H” is a function of time and location

          ( ,t) ( ,t) _{ } u(ξ, )H( ,t ;ξ,0)dξd t 0 s 0 x x x

– Unit slip rate at (X’=x, T’=0) gives stress perturbation “H” at (X=, T=t-)

– Values of “H” for every combination of (X,T;X’,T’) have to be calculated BEFORE computation.

– “H” can be analytically defined for simple geometries but it’s difficult for heterogeneous/complicated geometries.

(40)

SWoPP2008 40

Current Status (After FY.2003) (2/2)

Ch

f St t

Change of Strategy





















₍

_,

_t

₎

₍

_,

_t

₎

 

u

(

ξ

,

)

_H

₍

_,

_t

_;

ξ

_,

₀

₎

_d

ξd

t 0

x

steady plate

motion







slip perturbation

 

0 s

• Calculate “H” by FEM !!

          ( ,t) ( ,t) _{ } u(ξ, )H( ,t ;ξ,0)dξd t 0 s 0 x x x

(41)

41

FEMによる「単位すべり応答関数」計算

• GeoFEM

FEMによる「単位すべり応答関数」計算

• GeoFEM

• 分割節点法（Split Node Method）による，食い違い変位

（Di l

ti

）モデル

（Dislocation）モデル

– 断層をはさんだ変位の食い違い量を仮定し，仮想変位に相当する節点力を外力項として与える点力を外力項として与える – 特殊な要素は不要，定式化も簡単，観測値の利用 SWoPP2008

(42)

42

FEMによる「単位すべり応答関数」計算

FEMによる単位す

り応答関数」計算

（続き）

• 不均質粘弾性モデル

– Maxwell粘弾性体Maxwell粘弾性体





1 d

 









1

1 



dt

d

E

dt

d

• ちなみに弾性体は・・・

t t





E

1 

SWoPP2008

E

(43)

43

局所細分化の必要性

• 1km×1kmの応答関数を求めるためには数十m幅のメッ

局所細分化の必要性

• 1km×1kmの応答関数を求めるためには，数十m幅のメッ

シュが必要

数百k ×数百k の領域を扱うため全領域をこの大きさの

• 数百km×数百kmの領域を扱うため，全領域をこの大きさの

メッシュで切ることは不可能

局所細分化（）

– 局所細分化：Adaptive Mesh Refinement （AMR）

(44)

SWoPP2008 44

• AMRとData Migration

• 断層における応力蓄積シミュレーション

• 単位すべり応答関数計算用フレームワーク

• 評価

評価

• まとめ

(45)

45

単位すべり応答関数計算の現状

• 「粗い」初期全体メッシュ（O(10

5

~10

6

)要素）

単位す

り応答関数計算

現状

粗」初期

_{体ッシ（ (}

_)要素）

• 局所細分化

• 領域分割

• 六面体メッシュ（Tri-Linear）

– コミュニティーの要請 – 精度，アプリケーション

• 必要とする単位すべり応答関数のセットの数にもよるが，

必要とする単位すり応答関数のセットの数にもよるが，

数千～数万種類の計算が必要

– 高速 SWoPP2008 高速 – 局所細分化された分散メッシュファイルは不要

(46)

46

単位すべり応答関数計算用フレームワーク

• プロセス

単位す

り応答関数計算用

ワク

伊豆マントル伊豆マントル – 「粗い」分散メッシュ（数百m～1km） – 局所細分化四国伊四国伊 – 負荷分散 – 単位すべり応答関数フィリピン海プレート太平洋プレートマントルフィリピン海プレート太平洋プレートマントル x y z x y z

• 必要とする単位すべり応答関数のセットの数にもよるが，数千

必要とする単位すり応答関数のセットの数にもよるが，数千

～数万種類の計算が必要

– 高速高速 – 局所細分化された分散メッシュファイルは不要

• 六面体メッシュ（Tri-Linear）

SWoPP2008

六面体メッシュ（Tri-Linear）

(47)

SWoPP2008 47

単位すべり応答関数計算用フレームワーク

単位す

り応答関数計算用

ワク

Initial Distributed Mesh Files

AMR

Repartitioning (ParMETIS) Location of

Refinement

Data Migration & Management L ll R fi d Di t ib t d

Locally Refined Distributed Mesh Files (on memory)

Parallel FEM Parallel FEM (GeoFEM) Unit Response Unit Response Function

(48)

48

六面体要素のAMR

実は六面体要素のAMRは難しい

六面体要素のAMR

• 実は六面体要素のAMRは難しい

• 八分木（二次元では四分木）は細分化レベルの境界で

Hanging Nodeを生じてしまう

– 四面体，三角形ではこのようなことはないが，断層のモデリングには六面体要素が望ましく，六面体要素が一般的に用いられている SWoPP2008

(49)

Refinement Patterns of Tetrahedron

3 edges on a same 1 edge marked a same surface are marked 2 edges on a same 2 children (Binary) a same surface are marked 4 children 3rd edge is _{4 children} (Quadtree) more than 2 edges on different surfaces are marked 3 d edge s automatically marked marked more than SWoPP2008 49 8 children (Octree) more than 4 edges are

marked _{automatically}6 edges are marked

(50)

50

対処法：27分木 [Schneiders et al 1996]

対処法：27分木 [Schneiders et al, 1996]

• Quadtree/Octree with Conforming Refinement

• 各辺の分割単位を2ではなく，3とする。

各辺の分割単位を2ではなく，3とする。

– 細分化レベル境界でHanging Nodeを生じない – 二次元：9分木二次元：9分木

– 三次元：27分木

(51)

51

(a) (a)

Octree/Quadtree: Hanging Nodes to be

refined refined

(b)

Octree/Quadtree

with Conforming Refinement: NO Hanging Nodes

SWoPP2008

(52)

52

27分木における分割パターン

27分木における分割パタ

ン

[Schneiders et al, 1996]

A B C D SWoPP2008

(53)

53

27分木による細分化例

初期レベル

(54)

54

27分木による細分化例

細分化レベル1

(55)

55

27分木による細分化例

細分化レベル2

(56)

56

27分木による細分化例

細分化レベル3

(57)

57

27分木による細分化例

細分化レベル3（拡大図）

(58)

SWoPP2008 58

(59)

59 SWoPP2008 速度：背景色：位置：モード：透明度：

(60)

60

単位すべり応答関数計算用フレームワーク

• 並列細分化（AMR）

単位すべり応答関数計算用フレムワク

– Double Numberingの導入 [Nakajima & Fingberg, 2001] – 局所的な細分化情報の通信局所的な細分化情報通信 – 四面体と比較すると膨大な組み合わせを考慮する必要がある

• ParMETISによる領域再分割（Repartitioning）

ParMETISによる領域再分割（Repartitioning）

– http://glaros.dtc.umn.edu/gkhome/metis/parmetis/ マルチレベル的な手法節点数がバランスするような手法 – マルチレベル的な手法，節点数がバランスするような手法

– ParMETIS_V3_PartKway

• データマイグレーション（Data Migration）

– 節点座標，要素コネクティビティ，物性値，グループ情報，通信テー SWoPP2008 ブル等の領域間の移動 – 領域間オーバーラップ深さの自動変更機能（未実装）

(61)

SWoPP2008 61

細分化におけるルール（1/3）

• Double Numbering：並列計算向け

ロカル番号 – ID(i,1)=ローカル番号 – ID(i,2)=所属領域番号節点要素面辺 – 節点，要素，面，辺 – ローカル番号＋グローバル番号より便利

• 隣接要素

– 2レベル以上異ならない⇒自然に満たされる

(62)

SWoPP2008 62

細分化におけるルール（2/3）

• 連続して異方性のある分割を行わない

C 次は必ず分割を適用 – A,B,Cの次は必ずDの分割を適用 A B C D A B C D

(63)

SWoPP2008 63

細分化におけるルール（3/3）

• 節点を重複して生成するのを避けるため，辺，面，要

素に所属する節点を決めておく

(64)

64

M

E

T

I

S

http://www-users.cs.umn.edu/~karypis/metis/

チベ

グラ理論に基づた方法

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

SWoPP2008

(65)

65

M

E

T

I

S

http://www-users.cs.umn.edu/~karypis/metis/

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

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

– 特に通信（edge-cut）が少ない分割を提供する安定高速 – 安定，高速 – フリーウェア，他のプログラムに組み込むことも容易

• 色々な種類がある

– k-METIS 通信量（edge-cut）最小 – p-METIS 領域間バランス最適化 – ParMETIS 並列版 – 領域分割だけでなく，オーダリング，データマイニングなど色々な分野に使用されている SWoPP2008 • 接触，衝突問題における並列接触面探索

• ParMETIS：並列版

(66)

66

データマイグレーション（Data Migration）

• データマイグレーション（Data Migration）

デタマイグレション（Data Migration）

• デ

タマイグレション（Data Migration）

– 節点座標，要素コネクティビティ，物性値，グループ情報，通信テーブル等の領域間の移動ブル等の領域間の移動 – 領域間オーバーラップ深さの自動変更機能（未実装）

• 極力，All-to-All型の通信を避ける

• 隣接領域との通信のみに限定

SWoPP2008

(67)

SWoPP2008 67

• AMRとData Migration

• 断層における応力蓄積シミュレーション

• 単位すべり応答関数計算用フレームワーク

• 評価

評価

• まとめ

(68)

SWoPP2008 68

性能評価（ParMETIS+Data Migration）

性能評価（

_g

）

Initial Distributed Mesh Files

AMR

Repartitioning (ParMETIS) Location of

Refinement

Data Migration & Management L ll R fi d Di t ib t d

Locally Refined Distributed Mesh Files (on memory)

Parallel FEM Parallel FEM (GeoFEM) Unit Response Unit Response Function

(69)

SWoPP2008 69

例題

• 三次元弾性問題

初期メシ

例題

• 初期メッシュ

– 203 節点／コア，～512コア（4M節点）体体 – 立方体メッシュ（六面体）

• 細分化

– 4領域に1領域の割合で，1～2個の要素を細分化 • 16領域であればそのうち4領域が細分化 – 2レベル～4レベルケース細分化レベル最小節点数最大節点数平均節点数 1 2 8,000 17,290 10,324 2 3 68,560 23,140 3 3 , 247 648 67 192 3 3 247,648 67,192 4 4 812,320 209,080

(70)

SWoPP2008 70

性能評価手法

• Weak Scaling

16～512領域 – 16～512領域

• GeoFEM-type Parallel Distributed Data Structure

N d b d 1 L O l i – Node-based, 1-Layer Overlapping

安

• 評価の目安

(71)

SWoPP2008 71

(72)

SWoPP2008 72

ハードウェア等

• PCクラスタ（VT社製）

ハドウェア等

• PCクラスタ（VT社製）

– Dual-core Opteron 275 (2.2GHz) Cluster up to 64 cores Infiniband Gigabit Ethernet

– Infiniband, Gigabit-Ethernet

• T2Kオープンスパコン（東大）Hitachi HA8000

– Quad-core Opteron (2.3GHz) , up to 512 cores (32 nodes)

• F90+C+ MPI

– Flat MPI

– PGI（PCクラスタ） – Hitachi（HA8000）

(73)

SWoPP2008 73

PCクラスタ（Infiniband, Pathscale MPI）

ケース2（3レベル），初期実装

■ New Pointers ■ Data Migration 1.50 ■ Data Migration ■ Double Numbering ■ ParMETIS 1.00 s ec. 0.50 s 0.00 4 8 16 32 48 64 #

• Data MigrationはParMETISと比べて計算時間は少ないが

領域数増加とともに増大

core#

領域数増加とともに増大

(74)

SWoPP2008 74

デタマイグレシンの各プロセスと特徴

データマイグレーションの各プロセスと特徴

A：領域間通信，B：非スケーラブル処理

領域

通信，

非

名称機能 A B 名称機能 A B ■ New Local ID 新局所番号生成 ○ △ ■ Comm. Table 領域間通信テーブルサイズ決定 × △ サイズ決定 ■ Data Transfer 領域間データ移動 ◎ △ ■ Data Removal 不要データ消去 × △

(75)

SWoPP2008 75

非スケラブル処理

非スケーラブル処理

do ip= 1, PETOT do i= 1, N (…) enddo enddo enddo

(76)

SWoPP2008 76

Data Migrationの中での内訳

0.60 ■ New Local ID 0 40 0.50 ■ New Local ID

■ Comm. Tab. Size

■ Data Transfer 0.30 0.40 sec. ■ a a a s e ■ Data Removal 0.10 0.20 0.00 0.10 4 8 16 32 48 64 4 8 16 32 48 64 core#

(77)

SWoPP2008 77

Data Migration：改良点

名称変更点名称変更点 New Local ID 無し Comm Table 探索アルゴリズム効率化 Comm. Table 探索アルゴリズム効率化 Data Transfer 送受信バッファの初期化省略探索アルゴリズム効率化 Data Removal 探索アルゴリズム効率化送受信バッファの初期化省略

(78)

SWoPP2008 78

改良：64ノードで半分以下に

0 60 0 60

改良：64ノドで半分以下に

0.50 0.60 0.50 0.60 Original Improved 0.30 0.40 sec. 0.30 0.40 sec. 0.10 0.20 0.10 0.20 0.00 4 8 16 32 48 64 core# 0.00 4 8 16 32 48 64 core# ■ New Local ID

■ _{Comm. Tab. Size}

■ Data Transfer

(79)

SWoPP2008 79

Comm. Tab. Size

■ original 0.40 ■ original ■ improved 0.30 0.20 sec. 0 00 0.10 0.00 4 8 16 32 48 64 core# 名称変更点名称変更点 New Local ID 無し Comm. Table 探索アルゴリズム効率化 D T f 送受信バフの初期化省略 Data Transfer 送受信バッファの初期化省略 Data Removal 探索アルゴリズム効率化_{送受信バッファの初期化省略}

(80)

SWoPP2008 80

Data Transfer

0.15

Data Transfer

■ original 0.10 ■ original ■ improved 0.05 sec. 0 00 0.00 4 8 16 32 48 64 core# 名称変更点名称変更点 New Local ID 無し Comm. Table 探索アルゴリズム効率化 D T f 送受信バフの初期化省略 Data Transfer 送受信バッファの初期化省略 Data Removal 探索アルゴリズム効率化_{送受信バッファの初期化省略}

(81)

SWoPP2008 81

Data Removal

0.08

不要データ消去

■ original 0.06 ■ original ■ improved 0.04 sec. 0 00 0.02 0.00 4 8 16 32 48 64 core# 名称変更点名称変更点 New Local ID 無し Comm. Table 探索アルゴリズム効率化 D T f 送受信バフの初期化省略 Data Transfer 送受信バッファの初期化省略 Data Removal 探索アルゴリズム効率化_{送受信バッファの初期化省略}

(82)

SWoPP2008 82

PCクラスタネトワ

ク MPIの影響

PCクラスタ：ネットワーク，MPIの影響

■ IB+PathScale MPI ■ IB+MPICH ■ GB+MPICH ■ GB+MPICH 1.50 1.50 1.00 1.00

ParMETIS Data Migration

0.50 sec. 0.50 sec. 0.00 0.00 4 8 16 32 48 64 core# 0.00 4 8 16 32 48 64 core#

(83)

SWoPP2008 83

ケース1（T2K）

ケス1（T2K）

0.80 0.25 0.60 0.20 0.40 sec. 0.10 0.15 sec. 0 00 0.20 0 00 0.05 0.00 16 32 64 128 256 384 512 core# 0.00 16 32 64 128 256 384 512 core# ■ New Pointers ■ Data Migration ■ D bl N b i ■ New Local ID

■ Comm. Tab. Size ■ Double Numbering

■ ParMETIS

■ Data Transfer

(84)

SWoPP2008 84

ケース2（T2K）

0.40 1.50

@64cores ParMETIS：0.54sec.⇒0.57，DM：0.23⇒0.13

0.30 1.00 0.20 sec. 0.50 sec. 0 00 0.10 0 00 0.50 0.00 16 32 64 128 256 384 512 core# 0.00 16 32 64 128 256 384 512 core# ■ New Pointers ■ Data Migration ■ D bl N b i ■ New Local ID

■ ParMETIS

■ Data Transfer

(85)

SWoPP2008 85

ケース4（T2K）

4.00 5.00

ケス4（T2K）

3.00 4.00 2.00 sec. 2.00 3.00 sec. 0 00 1.00 0 00 1.00 0.00 16 32 64 128 256 384 512 core# 0.00 16 32 64 128 256 384 512 core# ■ New Pointers ■ Data Migration ■ D bl N b i ■ New Local ID

■ ParMETIS

■ Data Transfer

(86)

SWoPP2008 86

ケース4（T2K）

5.00

ケス4（T2K）

15.0 4.00 10.0 2.00 3.00 sec. 5.0 sec. 0 00 1.00 0 0 5.0 0.00 16 32 64 128 256 384 512 core# 0.0 16 32 64 128 256 384 512 core# ■ New Pointers ■ Data Migration ■ D bl N b i ■ New Local ID

■ ParMETIS

■ Data Transfer

(87)

SWoPP2008 87

デタマイグレシンの各プロセスと特徴

データマイグレーションの各プロセスと特徴

コア数増大による計算時間増大

数増

算

増

名称機能 1 2 3 4 名称機能 1 2 3 4 ■ New Local ID 新局所番号生成 ○ ○ ○ ○ ■ Comm. Table 領域間通信テーブルサイズ決定サイズ決定 ■ Data Transfer 領域間データ移動 ○ △ ■ Data Removal 不要データ消去 ○ ○ ○ ○

(88)

SWoPP2008 88

Example of “non-scalable” processes

Example of non scalable processes

in “NEW Local ID”

do ip= 1, PETOT

if (ip-1 eq my rank) then if (ip-1.eq.my_rank) then do i= 1, NODEtotLorg nn= NODE_ID_NEW(i,2) t( +1) t( +1) + 1 gcount(nn+1)= gcount(nn+1) + 1 NODE_ID_NEW(i,1)= gcount(nn+1) enddo endif

call MPI_BCAST ( gcount, PETOT, MPI_INTEGER, ip-1, &

MPI_COMM_WORLD, ierr )_ _

(89)

SWoPP2008 89

MAX Number of “destinations” for

data transfer

30 ● Case-1 ○ Case-2 ▲ Case-3 20 A TI O N △ Case-4 10 # : D E S T IN A 10 MA X # 0 10 100 1000 core#

(90)

SWoPP2008 90

• AMRとData Migration

• 断層における応力蓄積シミュレーション

• 単位すべり応答関数計算用フレームワーク

• 評価

評価

• まとめ

(91)

91

問題点

• 非スケーラブルプロセス

問題点

● Case-1 ○

非ケラル

– New Local ID – Data Removal ○ Case-2 ▲ Case-3 △ Case-4 Data Removal – 当面はこれでも問題無いが 1.30 1.40 tio

• ParMETIS

T2Kが遅い（64coresで 1.20 a la nce R a t – T2Kが遅い（64coresで 0.53sec⇒0.57sec.） – ロードインバランス 1.00 1.10 Load I m b a – ロドインバランス • ケース3，4はあまり現実的ではないが 0.90 10 100 1000 SWoPP2008 core#

(92)

92

今後の展望

• 改良

今後

展望

改良

– 非スケーラブル処理の部分

• 標準化

多様なデタ構造を受容きるようなもある必要がある – 多様なデータ構造を受容できるようなものである必要がある – かつ高速

• METIS，ParMETIS

– 他の関数，パラメータ – いつまでもこれに頼っていて良いのか？ SWoPP2008 までもれ頼てて良の

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

六面体メッシュの

六面体メッシュの

適合

並

細

適合

並

細

適合型並列局所細分化と

適合型並列局所細分化と

負荷分散

負荷分散

負荷分散

負荷分散

Conforming parallel adaptive mesh refinement for

hexahedral elements

hexahedral elements

中島 研吾

中島 研吾

概 要

概 要

• AMRとData Migration

• AMRとData Migration

• 断層における応力蓄積シミュレーション

• 単位すべり応答関数計算用フレームワーク

• 評価

評価

• まとめ

適応格子法（

適応格子法（AMR

AMR）（

）（1/2

1/2）

）

適応格子法（

適応格子法（AMR

AMR）（

）（1/2

1/2）

）

•

差分法 有限要素法などで空間を離散化する場合

差分法，有限要素法などで空間を離散化する場合，

精度良い解を求めるためには，細かいメッシュが必

要である

要である。

•

細かいメッシュを切ると，それだけ，計算量，必要記

憶容量が増加する。

•

細かいメッシュを必要とするのは 解析領域全体で

•

細かいメッシュを必要とするのは，解析領域全体で

はなく，ごく一部である。

衝撃波 渦 剥離

– 衝撃波，渦，剥離

– 応力集中，亀裂

– 従属変数が急激に変化する領域

適応格子法（

適応格子法（AMR

AMR）（

）（2/2

2/2）

）

適応格子法（

適応格子法（AMR

AMR）（

）（2/2

2/2）

）

•

しかしながら 多くの場合 これらの現象が生じる領

しかしながら，多くの場合，これらの現象が生じる領

域は，前もって分かっていない。

比較的粗いメ シ から始めて 計算結果に応じて

•

比較的粗いメッシュから始めて，計算結果に応じて，

局所的な細分化（refinement）を自動的に実施し，

精度良い解を効率的に求める手法が，適応格子法

中島研吾

中島研吾

概要

概要

差分法有限要素法などで空間を離散化する場合

細かいメッシュを必要とするのは解析領域全体で

衝撃波渦剥離

しかしながら多くの場合これらの現象が生じる領

比較的粗いメシから始めて計算結果に応じて

非構造格子の計算結果をブロク格子にマピングしてレンダ