三次元弾性解析コード (/3) プログラムの概要 22 年夏学期 中島研吾 科学技術計算 Ⅰ(82-27) コンピュータ科学特別講義 Ⅰ(8-2)
FEM3D-Part 2 対象とする問題 弾性体 Z ヤング率 E ポアソン比 ν U Z Z Z Y Y 直方体 一辺長さの立方体 ( 六面体 ) 要素 各方向に Y Z 個 境界条件 対称条件 U @ U Y @Y U Z @Z 強制変位 U Z @ZZ ma move
FEM3D-Part 3 構成式 : 応力 ~ひずみ関係 ヤング率 E 応力とひずみは比例 比例定数をヤング率 Eとする ( 各物質に固有の値 ) σ σ Eε ε E ポアソン比 ν 方向に荷重をかけると 横方向 (YZ) (YZ) にも変形 縮み割合をポアソン比 ν とする 各物質に固有の値 ε 金属では.3 程度 水 :.5 ゴム : ほぼ.5 非圧縮 νε σ ν E ε νε ε σ
FEM3D-Part YZν.3 とすると Z U Z Z ε ε. ε νε.3 Z u u ε. 3 Y Y Y Y
FEM3D-Part 5 有限要素法の処理 支配方程式 ガラーキン法 : 弱形式 要素単位の積分 要素マトリクス生成 全体マトリクス生成 境界条件適用 連立一次方程式
FEM3D-Part 6 初期化 有限要素法の処理 : プログラム 制御変数読み込み 座標読み込み 要素生成 (: 節点数 ICELO: 要素数 ) 配列初期化 ( 全体マトリクス 要素マトリクス ) 要素 全体マトリクスマッピング (IndeItem) マトリクス生成 要素単位の処理 (do cel ICELO) 要素マトリクス計算 全体マトリクスへの重ね合わせ 境界条件の処理 連立一次方程式 共役勾配法 (CG) 応力計算
FEM3D-Part 7 三次元要素の定式化 三次元弾性力学方程式 ガラーキン法 要素マトリクス生成 宿題 プログラムの実行 データ構造 プログラムの構成 計算効率について Paraew による可視化
FEM3D-Part 8 二次元への拡張 : 三角形要素 任意の形状を扱うことができる 特に一次要素は精度が悪く 一部の問題を除いてあまり使用されない
FEM3D-Part 9 二次元への拡張 : 四角形要素 一次元要素と同じ形状関数を 軸に適用することによって 四角形要素の定式化は可能である 三角形と比較して特に低次要素の精度はよい しかしながら 各辺が座標軸に平行な長方形でなければならない 差分法と変わらない 3 このような形状を扱うことができない 2
FEM3D-Part アイソパラメトリック要素 (/3) 各要素を 自然座標系 (η) の正方形要素 [± ±] に変換する 3 2 η 3 - - 2 各要素の全体座標系 (global coordnate)() における座標成分を 自然座標系における形状関数 [] ( 従属変数の内挿に使うのと同じ []) を使用して変換する場合 このような要素をアイソパラメトリック要素 (soparametrc element) という
FEM3D-Part アイソパラメトリック要素 (2/3) 3 η 3 2 - - 2 各節点の座標 :( ) ( 2 2 ) ( 3 3 ) ( ) 各節点における 方向変位 :u u 2 u 3 u u ( η) u η) ( ( η)
FEM3D-Part 2 アイソパラメトリック要素 (3/3) η - - 高次の補間関数を使えば 曲線 曲面も扱うことが可能となる そういう意味で 自然座標系 と呼んでいる
FEM3D-Part 3 2D 自然座標系の形状関数 (/2) 自然座標系における正方形上の内挿多項式は下式で与えられる : u α α 2 α3η α η 各節点での条件より : η 3 - - 2 u u2 u3 u α α 2 u u2 u3 u α3 α u u u2 u u2 u 3 3 u u
FEM3D-Part 2D 自然座標系の形状関数 (2/2) 元の式に代入して u について整理すると以下のようになる : η 3 u u u u u 2 2 3 3 u 形状関数 は以下のようになる : - 2 - ( η) ( )( η) 2( η) ( )( η) 3 ( η) ( )( η) ( η) ( )( η) 双一次 (b-lnear) 要素とも呼ばれる 各節点における の値を計算してみよ Y 方向変位 v についても同様 v v 2v2 3v3 v
FEM3D-Part 5 三次元への拡張 四面体要素 : 二次元における三角形要素 任意の形状を扱うことができる 特に一次要素は精度が悪く 一部の問題を除いてあまり使用されない 高次の四面体要素は広く使用されている 本講義では低次六面体要素 ( アイソパラメトリック要素 ) を使用する tr-lnear 変位法 自由度 : 変位 各節点上で 3 成分 (uvw) が定義される
FEM3D-Part 6 3D 自然座標系の形状関数 ( η ζ ) 8 2 ( η ζ ) η ζ 8 3 ( η ζ ) η ζ 8 ( η ζ ) η ζ 8 8 u ( η ζ ) u 8 v ( η ζ ) v 8 ( )( η)( ζ ) w ( η ζ ) w ( )( )( ) ( )( )( ) ( )( )( ) 5 ( η ζ ) 8 6 ( η ζ ) η ζ 8 7 ( η ζ ) η ζ 8 8 ( η ζ ) η ζ 8 ( ) 8 7 ( ) 5 6 ( ) ( η ζ ) ( ) ( ) 3 2 ( )( η)( ζ ) ( )( )( ) ( )( )( ) ( )( )( ) ( ) ( ) ( )
FEM3D-Part 7 三次元要素の定式化 三次元弾性力学方程式 ガラーキン法 要素マトリクス生成 宿題 プログラムの実行 データ構造 プログラムの構成 計算効率について Paraew による可視化
FEM3D-Part 8 弾性力学の支配方程式 つりあい式 ひずみ ~ 変位関係式 ひずみ ~ 応力関係式
FEM3D-Part 9 三次元のつりあい式な成応力の独立な成分は 6 つ τ τ { τ σ τ τ τ σ σ τ τ τ τ σ τ τ τ τ τ τ σ Y τ σ τ Z σ τ τ Z
FEM3D-Part 2 垂直ひずみ~ 変位の関係 PQ P Q ε u d u d ( u ) d d u d R P u / P d R v / Q Q ε ε ε u u v w
FEM3D-Part 2 せん断ひずみ~ 変位の関係 R R γ u v d u / P v / Q γ γ v w w u P d Q
FEM3D-Part 22 応力 ひずみ関係ず関係 σ σ ν ν ν ν ε ε ( ) E τ σ σ ν ν ν ν ν γ ε ε 2 ( ) ( ) ( ) E τ τ ν ν γ γ 2 2 2 ( ) τ ν γ 2
FEM3D-Part 23 ひずみ 応力関係ず関係 ν ν ν ε ε ν ν ν ν ν ν σ σ ( )( ) ( ) ( ) E γ ε ν ν ν ν τ σ 2 2 2 2 ( ) ( ) γ γ ν ν τ τ 2 2 2 2 2 [ ] D { [ ]{ ε σ D 非圧縮性材料 (ν~.5) の場合 特別な扱い必要
FEM3D-Part 2 - 方向のつりあい式に注目 σ τ τ σ τ τ ガラーキン法 [ ] { σ τ τ d [ ] { [ ] { [ ] { σ d τ d τ d [ ] { d
FEM3D-Part 25 - 方向のつりあい式 (/3) [ ] { σ [ ] { d σ d [ ] { σ Q [ ] { ] σ d [ ] { σ S n S ds n ds [ ] { ] [ ] { [ ] σ d σ d { σ 同様にして以下の弱形式が得られる ( 表面積分省略 ): 一次元ガウスの公式 d 部分積分の公式 [ ] { [ ] { [ ] { σ d τ d τ d [ ] { d [ ] { [ ] { [ ] σ d τ d { τ d [ ] { d
FEM3D-Part 26 - 方向のつりあい式 (2/3) [ ] { [ ] { [ ] d τ d { τ d [ ] { d (*) σ σ E ν ε νε ν 2ν νε ( )( ) ( ) [ ] E ν u νv w ν 2ν ν ( )( ) ( ) [ ] τ τ E 2 γ E [ ] u v ( ν ) 2( ν ) E E 2 γ [ ] u w ( ν ) 2( ν )
FEM3D-Part 27 要素内の変位分布形状関数要素内の変位分布 形状関数 [ ] [ ] [ ] { { { W w v U u [ ] [ ] [ ] { { { W w v U u [ ] [ ] [ ] { { { U u U u U u [ ] [ ] [ ] [ ] [ ] { { { { { { { v v U u U u U u [ ] [ ] { { W w W w E ( )( ) ( )[ ] [ ] [ ] ( ) { { { 2 E E W U E ν ν ν ν ν σ ( ) ( ) [ ] [ ] ( ) { { 2 2 U E E ν γ ν τ ( ) ( ) [ ] [ ] ( ) { { 2 2 W U E E ν γ ν τ
FEM3D-Part 28 - 方向のつりあい式 (3/3) 方向のつりあい式 (3/3) ( ) ( )( ) ( )( ) ( ) ν ν 2 2 2 E b E a E D とすると [ ] [ ] [ ] { { { W a a U D σ ( )( ) ( )( ) ( ) ν ν ν ν ν 2 2 2 [ ] [ ] [ ] [ ] [ ] [ ] [ ] { { { { { { { W b U b b U b τ τ [ ] [ ] { { W b U b τ () [ ] { [ ] { [ ] { [ ] { * d d d d τ τ σ [ ] [ ] [ ] [ ] [ ] [ ] ( ) { { U d b D [ ] [ ] [ ] [ ] ( ) { { [ ] [ ] [ ] [ ] { { W d b a d b a [ ] { d
FEM3D-Part 29 Y- 方向のつりあい式 [ ][ ] [ ] [ ] [ ] [ ] ( ) { { d b D [ ] [ ] [ ] [ ] ( ) { { [ ] [ ] [ ] [ ] { { [ ] { d Y W d b a U d b a [ ] { d Y Z- 方向のつりあい式 [ ][ ] [ ] [ ] [ ] [ ] ( ) { { d b U d b W d b D [ ] [ ] [ ] [ ] ( ) { { [ ] [ ] [ ] [ ] { { [ ] { d Z d b a U d b a [ ] { d Z
FEM3D-Part 3 3 つのつりあい式 変位の 3 成分を未知数 33 つの方程式 3つの方程式はカップルしている ( 独立では無い ) 形状関数ベクトル []:8 行列 [] [] [ ] [ ] 等 :8 8 行列 σ τ τ τ τ σ τ τ σ Y Z
FEM3D-Part 3 要素マトリクス :2 2 行列 各節点上の (uvw) 成分が物理的にも強くカップルしているので3 自由度をまとめて扱う :8 8 行列 j ( ) 8 7 ( ) ( ) 5 6 ) ( ) 3 ( ) 2 ( η ζ ) ( ) ( ) j j 2 j3 j2 j 22 j 32 j3 j 23 j 33 ( j K8)
FEM3D-Part 32 要素マトリクス :2 2 行列 各節点上の (uvw) 成分を独立に扱うことも可能であるが 8 8とする利点については来週以降も説明する j ( ) 8 7 ( ) ( ) 5 6 ) ( ) 3 ( ) 2 ( η ζ ) ( ) ( ) j j 2 j3 j2 j 22 j 32 j3 j 23 j 33 ( j K8)
FEM3D-Part 33 要素マトリクス :-j j 成分 方向 (/3) j j2 j 3 8 7 j 2 j3 j 22 j 32 j 23 j 33 ( j K 8 ) { D[ ][ ] b( [ ] [ ] [ ] [ ]) d{ U { a [ ] [ ] b ( [ ] [ ] ) d { { a[ ] [ ] b[ ] [ ] d{ W j ( ) ( ) 3 ( ) ( ) 5 6 ( ) 2 ( η ζ ) ( ) { D ( ) j b j j { a j b j j 2 { a j b j j3 d d d ( ) ( )
FEM3D-Part 3 要素マトリクス :-j j 成分 Y 方向 (2/3) j j 2 j3 j2 j 22 j 32 j3 j 23 j 33 ( j K 8 ) { D[ ][ ] b( [ ] [ ] [ ] [ ]) d{ { a [ ] [ ] b ( [ ] [ ] ) d { U { a[ ] [ ] b[ ] [ ] d{ W { a j b j j 2 j 22 ( ) 8 ( ) 7 3 ( ) ( ) 5 6 ( ) 2 ( η ζ ) ( ) d { D ( ) j b j j { a j b j j 23 d d ( ) ( )
FEM3D-Part 35 要素マトリクス :-j j 成分 Z 方向 (3/3) j j 2 j3 j2 j 22 j 32 j3 j 23 j 33 ( j K 8 ) { D[ ][ ] b( [ ] [ ] [ ] [ ]) d{ W { a [ ] [ ] b ( [ ] [ ] ) d { U { a[ ] [ ] b[ ] [ ] d{ { a j b j j 3 { a j b j j 32 j33 ( ) 8 ( ) 7 3 ( ) ( ) 5 6 ( ) 2 ( η ζ ) ( ) d d { D ( ) j b j j d ( ) ( )
FEM3D-Part 36 あとは積分を求めれば良いあ積分を求良 自然座標系 (ηζ) 上で定義 ガウス積分公式が ηζ 使える ( 三次元 ) しかし 微分が M L d d d f I ) ( ζ η ζ η [ ] j j M j L f W W W ) ( ζ η ) ( ζ LM: ηζ 方向の積分点数積分点座標値 j W W W ) ( j ζ η : 積分点での重み係数 : 積分点の座標値
FEM3D-Part 37 自然座標系における偏微分 (/) 自然座標系における偏微分 (/) 偏微分の公式より以下のようになる : ζ η ) ( η η η η ζ η ) ( ζ ζ ζ ζ ζ η η η η η ) ( ζ ζ ζ ζ は定義より簡単に求められるが ζ η は定義より簡単に求められるが を実際の計算で使用する
FEM3D-Part 38 自然座標系における偏微分 (2/) 自然座標系における偏微分 (2/) マトリックス表示すると : [ ] J [ ] ζ ζ ζ η η η ζ η ζ ζ ζ ζ [ ] J : ヤコビのマトリクス [ ] J : ヤコビのマトリクス (Jacob matr Jacoban)
FEM3D-Part 39 自然座標系における偏微分 (3/) 自然座標系における偏微分 (3/) の定義より簡単に求められる J J 8 8 2 8 8 J 8 8 3 J J 8 8 22 8 8 2 η η η η η η J 8 8 23 η η η J J 8 8 32 8 8 3 ζ ζ ζ ζ ζ ζ J 8 8 33 ζ ζ ζ
FEM3D-Part 自然座標系における偏微分 (/) 自然座標系における偏微分 (/) 従って下記のように偏微分を計算できる ヤコビアン (3 3 行列 ) の逆行列を求める [ ] η η η η η J ζ ζ ζ ζ ζ
FEM3D-Part 要素単位での積分 j j 2 j3 j2 j 22 j 32 j3 j 23 j 33 ( j K8) ( ) 8 7 ( ) 3 ( ) ( ) 5 6 ( ) 2 ( η ζ ) ( ) ( ) ( ) { ( ) j D j b j j d j j j D b d
FEM3D-Part 2 自然座標系での積分 j j 2 j3 j2 j 22 j 32 j3 j 23 j 33 ( j K8) ( ) 8 7 ( ) 3 ( ) ( ) 5 6 ( ) 2 ( η ζ ) ( ) ( ) ( ) j j j D b d j j j D b ddd D j j j b dt det J d dηdζζ
FEM3D-Part 3 ガウスの積分公式 j j 2 j3 j2 j 22 j 32 j3 j 23 j 33 ( j K8) ( ) 8 7 ( ) 3 ( ) ( ) 5 6 ( ) 2 ( η ζ ) ( ) ( ) ( ) D j b j j det J ddηdζ I L M j f ( η ζ ) ddηdζ [ W ] W j W f ( η j ζ )
FEM3D-Part 残りの手順 ここまでで 要素ごとの積分が可能となる あとは : 全体マトリクスへの重ね合わせ 境界条件処理 連立一次方程式を解く 詳細は来週以降の講義でプログラムの内容を解説 詳細は来週以降の講義でプログラムの内容を解説しながら説明する
FEM3D-Part 5 要素 全体マトリクス重ね合わせ要素全体トリク重ね合わ 5 6 3 (2) (2) ) ( ) ( ) ( ) ( (2) (2) (2) { ]{ [ f φ 2 5 6 2 (2) (2) 3 (2) 2 (2) (2) (2) 3 (2) 2 (2) (2) (2) (2) (2) (2) 3 (2) 33 (2) 32 (2) 3 (2) 2 (2) 23 (2) 22 (2) 2 (2) (2) 3 (2) 2 (2) f f f φ φ φ 3 3 2 () () () () () () () () () { ]{ [ f f φ φ (2) (2) (2) (2) 3 (2) 2 (2) f φ 2 2 2 () () 3 () 2 () () 3 () 2 () () 3 () 2 () () 3 () 33 () 32 () 3 () 2 () 23 () 22 () 2 3 2 f f f f φ φ φ φ 2 3 2 f φ Φ Φ { ]{ [ B D F K Φ Φ Φ Φ 3 2 3 2 3 2 B B B B D D D D Φ Φ Φ 6 5 6 5 6 5 B B B D D D
FEM3D-Part 6 三次元要素の定式化 三次元弾性力学方程式 ガラーキン法 要素マトリクス生成 宿題 プログラムの実行 データ構造 プログラムの構成 計算効率について Paraew による可視化
FEM3D-Part 7 宿題 ガウスの積分公式を使用して以下の四角形の面積を求めよ ( プログラムを作って 計算機で計算してください ) 3 2 :(..) 2:(. 2.) 3:(3. 5.) :(2..) I d det J ddζ ζ
FEM3D-Part 8 ヒント (/2) 座標値によってヤコビアン ( ヤコビの行列 ) を計算 ガウスの積分公式 (n2) に代入する I m n f j η ) d d η [ W ] W j f ( j ) ( η mplct REAL*8 (A-HO-Z) real*8 W(2) real*8 POI(2) W().d W(2).d POI() -.577352692d POI(2).577352692d SUM.d do jp 2 do p 2 FC F(POI(p)POI(jp)) SUM SUM W (p)*w (jp)*fc enddo enddo
FEM3D-Part 9 ヒント (2/2) [ ] η η η η J J det η η η η η η η η ( )( ) ( )( ) η η η η ) ( ) ( 2 ( )( ) ( )( ) η η η η ) ( ) ( 3
FEM3D-Part 5 三次元要素の定式化 三次元弾性力学方程式 ガラーキン法 要素マトリクス生成 宿題 プログラムの実行 データ構造 プログラムの構成 計算効率について Paraew による可視化
FEM3D-Part 5 対象とする問題 弾性体 Z ヤング率 E ポアソン比 ν U Z Z Z Y Y 直方体 一辺長さの立方体 ( 六面体 ) 要素 各方向に Y Z 個 境界条件 対称条件 U @ U Y @Y U Z @Z 強制変位 U Z @ZZ ma move
FEM3D-Part 52 YZν.3 とすると Z U Z Z ε ε. ε νε.3 Z u u ε. 3 Y Y Y Y
FEM3D-Part 53 ファイルコピー インストール (/2) 三次元弾性解析コード >$ cd <$fem> >$ cp /home3/sengon/documents/class/fem/fem3d.tar. >$ tar vf fem3d.tar >$ cd fem3d >$ ls c f run FEM インストール (C) >$ cd <$fem>/fem3d/c >$ mae >$ ls../run/sol sol FEM インストール (FORRA) >$ cd <$fem>/fem3d/f >$ mae >$ ls../run/sol sol
FEM3D-Part 5 ファイルコピー インストール (2/2) メッシュジェネレータインストール >$ cd <$fem>/fem3d/run >$ g95 O3 mgcube.f o mgcube
FEM3D-Part 55 計算の流れメッシュ生成 計算 ファイル名称固定 mgcube メッシュジェネレータ cube. メッシュファイル sol FEM 本体 IPU.DA 制御データ test.np 可視化用出力変位 3 成分 σ
FEM3D-Part 56 メッシュ生成 Z U Z >$ cd <$fem>/fem3d/run / >$./mgcube Y Z 各辺長さを訊いてくる このように 入れてみる Z >$ ls cube. 生成を確認 Y Y cube.
FEM3D-Part 57 制御ファイル :IPU.DA IPU.DA cube. fname MEHOD PRECOD terprema( 不使用 ) 2 IER..3 ELAS POISSO fname: メッシュファイル名 MEHOD: 反復解法 ( に固定 ) PRECOD: 前処理手法 :Bloc-LU-GS:Bloc 対角スケーリング terprema: ( 不使用 ) IER: ELAS: 反復回数上限 ヤング率 POISSO: ポアソン比 あとで.999などの場合を試してみよ
FEM3D-Part 58 >$ cd <$fem>/fem3d/run >$./sol ( 反復法の収束履歴 ) 33 2.28867E-8 3.32592E-8 35 7.383E-9 実行 ### DISPLACEME at (maymazma)) 33-3.E- -3.2E-.E >$ ls test.np 生成を確認 test.np Z U Z この点での変位が表示されている 33 番節点 ( 3 ) Z Y Y
FEM3D-Part 59 三次元要素の定式化 三次元弾性力学方程式 ガラーキン法 要素マトリクス生成 宿題 プログラムの実行 データ構造 プログラムの構成 計算効率について Paraew による可視化
FEM3D-Part 6 メッシュファイル構成 :cube. 番号は から始まっている 節点データ 節点数 節点番号 座標 要素データ 要素数 要素タイプ 要素番号 材料番号 コネクティビティ 節点グループデータ グループ数 グループ内節点数 グループ名 グループ内節点
FEM3D-Part 6 メッシュ生成コード :mgcube.f(/5) mplct REAL*8 (A-HO-Z) real(nd8) dmenson(::) allocatable :: Y real(nd8) dmenson(::) allocatable :: Y character(len8) :: GRIDFILE HHH nteger dmenson(::) allocatable :: IW!C!C -------!C II.!C -------!C wrte (**) 'YZ' read (**) YZ P YP Y ZP Z D.d IODO P*YP*ZP ICELO *Y *Z IBODO P*YP allocate (IW(IODO)) IW 方向節点数 Y 方向節点数 Z 方向節点数 総節点数総要素数 Y 平面の節点数
FEM3D-Part 62 メッシュ生成コード :mgcube.f(2/5)!c cou b do ZP do j YP cou cou (-)*IBODO (j-)*p IW(coub) enddo enddo cou b 2 do ZP j do P cou cou (-)*IBODO (j-)*p IW(coub) enddo enddo cou b 3 do j YP do P cou cou (-)*IBODO (j-)*p IW(coub) enddo enddo cou b ZP do j YP do P cou cou (-)*IBODO (j-)*p IW(coub) enddo enddo mn の節点を IW(b) に格納 (byp*zp) Z Y U Z Z Y
FEM3D-Part 63 メッシュ生成コード :mgcube.f(2/5)!c cou b do ZP do j YP cou cou (-)*IBODO (j-)*p IW(coub) enddo enddo cou b 2 do ZP j do P cou cou (-)*IBODO (j-)*p IW(coub) enddo enddo cou b 3 do j YP do P cou cou (-)*IBODO (j-)*p IW(coub) enddo enddo cou b ZP do j YP do P cou cou (-)*IBODO (j-)*p IW(coub) enddo enddo YYmnの節点をIW(b2) に格納 (bp*zp) ZP) Z U Z Z Y Y
FEM3D-Part 6 メッシュ生成コード :mgcube.f(2/5)!c cou b do ZP do j YP cou cou (-)*IBODO (j-)*p IW(coub) enddo enddo Z U Z cou b 2 do ZP Z j do P cou cou (-)*IBODO (j-)*p IW(coub) Y enddo enddo cou Y b 3 do j YP do P cou cou (-)*IBODO (j-)*p IW(coub) enddo enddo cou b ZP do j YP do P cou cou (-)*IBODO (j-)*p IW(coub) enddo enddo ZZmn の節点を IW(b3) に格納 (bp*yp) ZZmaの節点をIW(b) に格納 (bp*yp) YP)
FEM3D-Part 65 メッシュ生成コード :mgcube.f(3/5)!c!c -------------!C GeoFEM data!c -------------!C wrte (**) 'GeoFEM grdfle name?' GRIDFILE 'cube.' open (2 fle GRIDFILE status'unnown'form'formatted') wrte(2'()') IODO cou do ZP do j YP do P dfloat(-)*d YY dfloat(j-)*d ZZ dfloat(-)*d cou cou wrte (2'(3(pe6.6))') cou YY ZZ enddo enddo enddo wrte(2'()') ICELO IELMYPL 36 wrte(2'()') (IELMYPL ICELO) 節点数節点番号 座標 要素数要素タイプ ( 使わないデータ )
FEM3D-Part 66 cube.: 節点データ 要素数 要素タイプ (YZ) 25 5*5*5.E.E.E 2.E.E.E 3 2.E.E.E 3.E.E.E 5.E.E.E 6.E.E.E 7.E.E.E 8 2.E.E.E 9 3.E.E.E ( 途中省略 ) 2.EE.E.E 22.E.E.E 23 2.E.E.E 2 3.E.E.E 25.E.E.E 6 ** 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 move
FEM3D-Part 67 メッシュ生成コード :mgcube.f(/5) cou mat do Z do j Y do cou cou n (-)*IBODO (j-)*p n2 n n3 n2 P n n3 - n5 n IBODO n6 n2 IBODO n7 n3 IBODO n8 n IBODO wrte (2'()') cou mat n n2 n3 n & & n5 n6 n7 n8 enddo enddo enddo mat: 材料番号 () 8 7 5 6 3 2
FEM3D-Part 68 cube.: 要素データ 2 7 6 26 27 32 3 2 2 3 8 7 27 28 33 32 3 3 9 8 28 29 3 33 5 9 29 3 35 3 5 6 7 2 3 32 37 36 6 7 8 3 2 32 33 38 37 7 8 9 3 33 3 39 38 8 9 5 3 35 39 9 2 7 6 36 37 2 2 3 8 7 37 38 3 2 3 9 8 38 39 3 2 5 2 9 39 5 3 6 7 22 2 2 7 6 ( 途中省略 ) 2 62 63 68 67 87 88 93 92 3 63 6 69 68 88 89 9 93 6 65 7 69 89 9 95 9 5 66 67 72 7 9 92 97 96 6 67 68 73 72 92 93 98 97 7 68 69 7 73 93 9 99 98 8 69 7 75 7 9 95 99 9 76 77 82 8 2 7 6 5 77 78 83 82 2 3 8 7 5 78 79 8 83 3 9 8 52 79 8 85 8 5 9 53 8 82 87 86 6 7 2 5 82 83 88 87 7 8 3 2 55 83 8 89 88 8 9 3 56 8 85 9 89 9 5 57 86 87 92 9 2 7 6 58 87 88 93 92 2 3 8 7 59 88 89 9 93 3 9 8 6 89 9 95 9 5 2 9 6 9 92 97 96 6 7 22 2 62 92 93 98 97 7 8 23 22 63 93 9 99 98 8 9 2 23 6 9 95 99 9 2 25 2
FEM3D-Part 69 メッシュ生成コード :mgcube.f(5/5) IGO IB YP*ZP IB2 P*ZP IB IB3 P*YP IB2 IB P*YP IB3 wrte (2'()') IGO wrte (2'()') IB IB2 IB3 IB HHH 'mn' wrte (2'(a8)') HHH wrte (2'()') (IW() YP*ZP) HHH 'Ymn' wrte (2'(a8)') HHH wrte (2'()') ( ) (IW(2) ( P*ZP) HHH 'Zmn' wrte (2'(a8)') HHH wrte (2'()') (IW(3) P*YP) HHH 'Zma' wrte (2'(a8)') HHH wrte (2'()') (IW() P*YP) IGO IB グループ総数 (mnymnzmnzma) Zma) 累積数 ( 以下略 )!C stop end deallocate (IW) close (2)
FEM3D-Part 7 cube.: 節点グループデータ mn Ymn Zmn Zma 25 5 75 6 6 2 26 3 36 6 5 56 6 66 7 76 8 86 9 96 6 6 2 2 3 5 26 27 28 29 3 5 52 53 5 55 76 77 78 79 8 2 3 5 2 3 5 6 7 8 9 2 3 5 6 7 8 9 2 2 22 23 2 25 2 3 5 6 7 8 9 2 3 5 6 7 8 9 2 2 22 23 2 25 ( 以下使用せず )
FEM3D-Part 7 メッシュ生成 実は技術的には大きな課題 複雑形状 大規模メッシュ 並列化が難しい 市販のメッシュ生成アプリケーション FEMAP CAD データとのインタフェース move
FEM3D-Part 72 三次元要素の定式化 三次元弾性力学方程式 ガラーキン法 要素マトリクス生成 宿題 プログラムの実行 データ構造 プログラムの構成 計算効率について Paraew による可視化
FEM3D-Part 73 有限要素法の処理 : プログラム 初期化 制御変数読み込み 座標読み込み 要素生成 (: 節点数 E: 要素数 ) 配列初期化 ( 全体マトリクス 要素マトリクス ) 要素 全体マトリクスマッピング (IndeItem) マトリクス生成 要素単位の処理 (do cel E) 要素マトリクス計算 全体マトリクスへの重ね合わせ 境界条件の処理 連立一次方程式 共役勾配法 (CG) 応力計算
FEM3D-Part 7 test メインプログラム nput_cntl 制御データ入力 nput_grd メッシュファイル入力 fnd_node 節点探索 mat_con 行列コネクティビティ生成 msor ソート mat_con 行列コネクティビティ生成 mat_ass_manass man 係数行列生成 jacob ヤコビアン計算 mat_ass_bc 境界条件処理 三次元弾性解析コードfem3d の構成 solve33 線形ソルバー制御 recover_stress 応力計算 output_ucd 可視化処理 cg_3 CG 法計算 jacob ヤコビアン計算
FEM3D-Part 75 全体処理 #nclude <stdo.h> #nclude <stdlb.h> FILE* fp_log; #defne GLOBAL_ALUE_DEFIE #nclude "pfem_utl.h" etern vod IPU_CL(); etern vod IPU_GRID(); etern vod MA_CO(); etern vod MA_CO(); etern vod MA_ASS_MAI(); etern vod MA_ASS_BC(); ASS etern vod SOLE33(); etern vod RECOER_SRESS(); etern vod OUPU_UCD(); nt man() { /** Logfle for debug **/ f( (fp_logfopen("log.log""w")) ULL){ fprntf(stdout"nput fle cannot be opened! n"); et(); IPU_CL(); IPU_GRID(); MA_CO(); MA_CO(); MA_ASS_MAI(); MA_ASS_BC() ; SOLE33(); RECOER_SRESS(); OUPU_UCD() ;
FEM3D-Part 76 Global 変数表 :pfem_utl.h/f(/3) 変数名種別サイズ I/O 内容 fname C [8] I メッシュファイル名 P I I 節点数 ICELO I I 要素数 ODGRPtot I I 節点グループ数 YZ R [][3] I 節点座標 ICELOD I [ICELO][8] I 要素コネクティビティ ODGRP_IDE I [ODGRPtot] I 各節点グループに含まれる節点数 ( 累積 ) ODGRP_IEM I [ODGRP_IDE[ODG RPO]] I 節点グループに含まれる節点 ODGRP_AME C8 [ODGRP_IDE[ODG RPO]] I 節点グループ名 L U I O 各節点非対角成分数 ( 上三角 下三角 ) PL PU I O 非対角成分総数 ( 上三角 下三角 ) D R [9*] O 全体行列 : 対角ブロック B R [3*] O 右辺ベクトル 未知数ベクトル
FEM3D-Part 77 Global 変数表 :pfem_utl.h/f(2/3) 変数名種別サイズ I/O 内容 ALUG R [9*] O 全体行列 : 対角ブロックの完全 LU 分解 ALAU R [9*PL][9*PU] O 全体行列 : 上 下三角成分 ndel ndeu I [] O 全体行列 : 非零非対角ブロック数 teml temu I [PL] [PU] O 全体行列 : 上 下三角ブロック ( 列番号 ) IL IU I [] O 各節点の上 下三角ブロック数 IALIAUIAU I [][L][][U][][U] O 各節点の上 下三角ブロック ( 列番号 ) IWK I [][2] O ワーク用配列 MEHOD I I 反復解法 ( に固定 ) PRECOD I I 前処理手法 (: ブロック SSOR: ブロック対角スケーリング ) IER IERactual I I 反復回数の上限 実際の反復回数 RESID R I 打ち切り誤差 (.e-8 に設定 ) SIGMA_DIAG R I LU 分解時の対角成分係数 (. に設定 ) pfemiarra I [] O 諸定数 ( 整数 ) pfemrarra R [] O 諸定数 ( 実数 )
FEM3D-Part 78 Global 変数表 :pfem_utl.h/f(3/3) 変数名種別サイズ I/O 内容 O8th R I.25 PQ PE P R [2][2][8] O 各ガウス積分点における η ζ ( ~ 8) POS WEI R [2] O 各ガウス積分点の座標 重み係数 COL COL2 I [] O ソート用ワーク配列 SHAPE R [2][2][2][8] O 各ガウス積分点における形状関数 (~8) P PY PZ R [2][2][2][8] O 各ガウス積分点における ( ~ 8 ) DEJ R [2][2][2] O 各ガウス積分点におけるヤコビアン行列式 ELAS POISSO R I ヤング率 ポアソン比 SIGMA_ AU_ R [][3] O 節点における垂直 せん断応力成分
FEM3D-Part 79 制御ファイル入力 :IPU_CL #nclude <stdo.h> #nclude <stdlb.h> #nclude "pfem_utl.h" /** **/ vod IPU_CL() { FILE *fp; f( (fpfopen("ipu.da""r")) ULL){ fprntf(stdout"nput fle cannot be opened! n"); et(); fscanf(fp"%s"fname); f(f " f fscanf(fp"%d %d"&mehod&precod); fscanf(fp"%d"&terprema); fscanf(fp"%d"&ier); fscanf(fp "%lf %lf" &ELAS &POISSO); fclose(fp); f( ( terprema < ) ){ terprema ; f( ( terprema > ) ){ terprema ; SIGMA_DIAG.; SIGMA.; RESID.e-8; SE ; pfemrarra[] RESID; pfemrarra[] SIGMA_DIAG; pfemrarra[2] SIGMA; IPU.DA cube. fname MEHOD PRECOD 2 terprema( 不使用 ) IER..3 ELAS POISSO pfemiarra[] IER; pfemiarra[] MEHOD; pfemiarra[2] PRECOD; pfemiarra[3] SE; pfemiarra[] terprema;
FEM3D-Part 8 メッシュ入力 :IPU_GRID(/2) #nclude <stdo.h> #nclude <stdlb.h> #nclude "pfem_utl.h" #nclude "allocate.h" vod IPU_GRID() { FILE *fp; nt jnncelse; nt YPEIMA; /** **/ f( (fpfopen(fname"r")) ULL){ fprntf(stdout"nput f( fle cannot be opened! n"); et(); ODE fscanf(fp"%d"&); P; YZ(KREAL**)allocate_matr(seof(KREAL)3); for(;<;){ for(j;j<3;j){ YZ[][j].; YZ[][3] for(;<;){ for( ;<;){ fscanf(fp"%d %lf %lf %lf"&&yz[][]&yz[][]&yz[][2]);
FEM3D-Part 8 allocate deallocate 関数 #nclude <stdo.h> #nclude <stdlb.h> vod* allocate_vector(nt sent m) { vod d*; *a; f ( ( a(vod * )malloc( m * se ) ) ULL ) { fprntf(stdout"error:memor does not enough! n vector n"); et(); return a; allocate を FORRA 並みに簡単にやるための関数 vod deallocate_vector(vod *a) { free( a ); vod** allocate_matr(nt sent mnt n) { vod **aa; nt ; f((aa(vod ( (vod ** )malloc( m * seof(vod*) )) ) ULL ) { fprntf(stdout"error:memor does not enough! aa n matr n"); et(); f ( ( aa[](vod * )malloc( m * n * se ) ) ULL ) { fprntf(stdout"error:memor does not enough! n matr n"); et(); for(;<m;) aa[](char*)aa[-]se*n; return aa; vod deallocate matr(vod **aa) vod deallocate_matr(vod **aa) { free( aa );
FEM3D-Part 82 メッシュ入力 :IPU_GRID(2/2) /** **/ ELEME fscanf(fp"%d"&icelo); ICELOD(KI**)allocate_matr(seof(KI)ICELO8); t( ICELO for(;<icelo;) fscanf(fp"%d"&ype); ICELOD[][j] の中身としては から始まる通し節点番号がそのまま読み込まれている 要素番号は から番号付け /** **/ for(cel;cel<icelo;cel){ fscanf(fp"%d %d %d %d %d %d %d %d %d %d"&&ima &ICELOD[cel][]&ICELOD[cel][]&ICELOD[cel][2]&ICELOD[cel][3] l][] l][2] l][3] &ICELOD[cel][]&ICELOD[cel][5]&ICELOD[cel][6]&ICELOD[cel][7]); ODE grp. nfo. fscanf(fp"%d"&odgrptot); ODGRP_IDE(KI* )allocate_vector(seof(ki)odgrptot); ODGRP_AME (CHAR8*)allocate_vector(seof(CHAR8)ODGRPtot); for(;<odgrptot;) ODGRP_IDE[]; for(;<odgrptot;) fscanf(fp"%d"&odgrp_ide[]); nnodgrp_ide[odgrptot]; ODGRP_IEM(KI*)allocate_vector(seof(KI)nn); for(;<odgrptot;){ S ODGRP_IDE[]; E ODGRP_IDE[]; fscanf(fp"%s"odgrp_ame[].name); nn E - S; f( nn! ){ for(s;<e;) fscanf(fp"%d"&odgrp_iem[]); fclose(fp);
FEM3D-Part 83 メッシュ入力 :IPU_GRID(2/2) /** **/ ELEME fscanf(fp"%d"&icelo); ICELOD(KI**)allocate_matr(seof(KI)ICELO8); t( ICELO for(;<icelo;) fscanf(fp"%d"&ype); /** **/ for(cel;cel<icelo;cel){ fscanf(fp"%d %d %d %d %d %d %d %d %d %d"&&ima &ICELOD[cel][]&ICELOD[cel][]&ICELOD[cel][2]&ICELOD[cel][3] l][] l][2] l][3] &ICELOD[cel][]&ICELOD[cel][5]&ICELOD[cel][6]&ICELOD[cel][7]); ODE grp. nfo. fscanf(fp"%d"&odgrptot); ODGRP_IDE(KI* )allocate_vector(seof(ki)odgrptot); ODGRP_AME (CHAR8*)allocate_vector(seof(CHAR8)ODGRPtot); for(;<odgrptot;) ODGRP_IDE[]; for(;<odgrptot;) fscanf(fp"%d"&odgrp_ide[]); nnodgrp_ide[odgrptot]; ODGRP_IEM(KI*)allocate_vector(seof(KI)nn); for(;<odgrptot;){ S ODGRP_IDE[]; E ODGRP_IDE[]; fscanf(fp"%s"odgrp_ame[].name); nn E - S; f( nn! ){ for(s;<e;) fscanf(fp"%d"&odgrp_iem[]); 節点グループの中身も から始まる通し節点番号がそのまま読み込まれている fclose(fp);
FEM3D-Part 8 fem3d: いくつかの特徴 非対角成分 上三角 下三角を分けて記憶 U ブロックとして記憶 ベクトル : 節点 3 成分 L 行列 : 各ブロック 9 成分 行列の各成分ではなく 節点上の3 変数に基づくブロックとして処理する
有限要素法で得られるマトリクス FEM3D-Part 85 疎行列 Φ Φ Φ 3 2 3 2 F F F D D D 疎行列 が多い A( j) のように正方行列の Φ Φ Φ Φ 7 6 5 3 7 6 5 3 F F F F D D D D A(j) のように正方行列の全成分を記憶することは疎行列では非効率的 Φ Φ Φ Φ Φ 2 9 8 2 9 8 F F F F F D D D D D 行列では非効率的 密 行列向け Φ Φ Φ Φ Φ 6 5 3 2 6 5 3 2 F F F F F D D D D D 有限要素法 : 非零非対角成分の数は高々 数百 規模 例えば未知数が 8 個あるとすると記憶容量 ( ワード数 ) は例えば未知数が 個あるとすると記憶容量 ( ワド数 ) は 正方行列 :O( 6 ) 非零非対角成分数 :O( ) 非零成分のみ記憶するのが効率的
FEM3D-Part 86 d.fd.c におけるマトリクス関連変数 変数名型サイズ内容 I - 未知数総数 PLU I - 連立一次方程式係数マトリクス非対角成分総数 Dag(:) R 連立一次方程式係数マトリクス対角成分 U(:) R 連立一次方程式未知数ベクトル Rhs(:) R 連立一次方程式右辺ベクトル Inde(:) I : Item(:) I PLU AMat(:) R PLU 係数マトリクス非対角成分要素番号用一次元圧縮配列 ( 非対角成分数 ) 係数マトリクス非対角成分要素番号用一次元圧縮配列 ( 非対角成分要素 ( 列 ) 番号 ) 係数マトリクス非対角成分要素番号用一次元圧縮配列 ( 非対角成分 ) 非零非対角成分のみを格納する Compressed Row Storage 法を使用している
行列ベクトル積への適用 ( 非零 ) 非対角成分のみを格納疎行列向け方法 FEM3D-Part 87 ( 非零 ) 非対角成分のみを格納 疎行列向け方法 Compressed Row Storage (CRS) 対角成分 ( 実数 ) Dag () 対角成分 ( 実数 ) Inde() 非対角成分数に関する一次元配列 ( 通し番号 ) ( 整数 ) () 非対角成分の要素 ( 列 ) 番号 Item() 非対角成分の要素 ( 列 ) 番号 ( 整数 nde()) AMat() 非対角成分 ( 実数 d ()) ( 実数 nde()) {Y [A]{ Φ Φ 2 2 F F D D do Y() Dag()*() d I d ( ) I d () Φ Φ Φ Φ Φ Φ 8 7 6 5 3 8 7 6 5 3 F F F F F F D D D D D D do Inde(-) Inde() Y() Y() Amat()*(Item()) enddo Φ Φ Φ Φ Φ Φ 3 2 9 3 2 9 F F F F F F D D D D D D enddo Φ Φ 6 5 6 5 F F D D
FEM3D-Part 88 行列ベクトル積への適用 ( 非零 ) 非対角成分のみを格納 疎行列向け方法 Compressed Row Storage (CRS) {Q[A]{P for(;<;){ W[Q][] Dag[] * W[P][]; for(inde[];<inde[];){ W[Q][] AMat[]*W[P][Item[]]; []]
FEM3D-Part 89 行列ベクトル積 : 密行列 とても簡単 a a2... a a a a a a 2 22 2 2......... 2 2 M M a a a a 2 a a2... a a {Y [A]{ do j Y(j).d do Y(j) Y(j) A(j)*() enddo enddo
FEM3D-Part 9 Compressed Row Storage (CRS) 2 3 5 6 7 8. 2. 3.2 2.3 3.6 2.5 3.7 9. 3 5.7.5 3.. 9.8 2.5 2.7 5 3. 9.5..5.3 6 6.5 2. 9.5 7 6. 2.5. 23. 3. 8 9.5.3 9.6 3. 5. 3
FEM3D-Part 9 Compressed Row Storage (CRS):C 2 3 5 6 7 2 3 5 6 7. 2 2. 32 3.2 8.3 3.6 2.5 3.7 9. 対角成分 3 5 7 Dag[]. 5.7.5 3. Dag[] 3.6 2 6 Dag[2] 5.7. 98 9.8 25 2.5 27 2.7 Dag[3] 9.8 3 5 Dag[].5 Dag[5] 2. 3. 9.5..5.3 Dag[6] 23. 2 6 Dag[7] 5.3 6.5 2. 9.5 2 5 6 6 6. 25 2.5. 23. 3. 2 5 6 7 9.5.3 9.6 3. 5.3 2 3 5 7
FEM3D-Part 92 Compressed Row Storage (CRS):C 2 3 5 6 7 2 3 5 6 7. 2 2. 32 3.2 3.6.3 2.5 3.7 9. 3 5 7 5.7.5 3. 2 6 98 9.8. 25 2.5 27 2.7 3 5.5 3. 9.5..3 2 6 2. 6.5 9.5 5 2 6 23. 6 6. 25 2.5. 3. 6 2 5 7 5.3 9.5.3 9.6 3. 7 2 3 5
FEM3D-Part 93 Compressed Row Storage (CRS):C 非対角成分数 Inde[]. 3.6 2 2..3 32 3.2 2.5 3.7 9. 3 5 7 2 Inde[] 2 Inde[2] 6 2 5.7.5 3. 2 6 2 Inde[3] 8 3 98 9.8. 25 2.5 27 2.7 3 5 3 Inde[].5 3. 9.5..3 2 6 Inde[5] 5 5 2. 6.5 9.5 5 2 6 2 Inde[6] 7 6 23. 6 6. 25 2.5. 3. 6 2 5 7 Inde[7] 2 7 5.3 9.5.3 9.6 3. Inde[8] 25 PLU 25 7 2 3 5 (Inde[]) Inde[]~Inde[]- 番目が 行目の非対角成分
FEM3D-Part 9 Compressed Row Storage (CRS):C 2 3 5 6 7 非対角成分数 Inde[]. 2 2. 32 3.2 2 Inde[] 2 3.6.3 2.5 3.7 9. 3 5 7 Inde[2] 6 5.7.5 3. 2 6 2 Inde[3] 8 98 9.8. 25 2.5 27 2.7 3 5 3 Inde[].5 3. 9.5..3 2 6 Inde[5] 5 2. 6.5 9.5 5 2 6 2 Inde[6] 7 23. 6 6. 25 2.5. 3. 6 2 5 7 Inde[7] 2 5.3 9.5.3 9.6 3. 7 2 3 5 Inde[8] 25 Inde[]~Inde[]- 番目が 行目の非対角成分 PLU 25 (Inde[])
FEM3D-Part 95 Compressed Row Storage (CRS):C 2 3 5 6 7. 2 2. 32 3.2 3.6.3 2.5 3.7 9. 3 5 7 5.7.5 3. 2 6 98 9.8. 25 2.5 27 2.7 3 5.5 3. 9.5..3 2 6 2. 6.5 9.5 5 2 6 23. 6 6. 25 2.5. 3. 6 2 5 7 5.3 9.5.3 9.6 3. 7 2 3 5 例 : Item[ 6] AMat[ 6].5 Item[8] 2 AMat[8] 2.5
FEM3D-Part 96 Compressed Row Storage (CRS):C 2 3 5 6 7. 2 2. 32 3.2 3.6.3 2.5 3.7 9. Dag [] 対角成分 ( 実数 []) Inde[] 非対角成分数に関する一次元配列 22 333 5 755 ( 通し番号 )( 整数 []) 5.7.5 3. Item[] 非対角成分の要素 ( 列 ) 番号 2 6 67 ( 整数 [Inde[]]) 98 9.8. 25 2.5 27 2.7 AMat[] 非対角成分 3 8 9 5 ( 実数 [Inde[]]).5 3. 9.5..3 2236 3 {Y[A]{ 2. 6.5 9.5 for(;<;){ 5 25 66 Y[] Dag[] * []; 23. 6 6. 25 2.5. 3. for(inde[];<inde[];){ 6 7 28 5972 Y[] AMat[]*[Item[]]; 5.3 9.5.3 9.6 3. 7 2 222 32352 2
FEM3D-Part 97 ブロックとして記憶 (/3) 記憶容量が減る Inde Item に関する記憶容量を削減できる j j j j
FEM3D-Part 98 ブロックとして記憶 (2/3) 計算効率 間接参照 ( メモリに負担 ) と計算の比が大きくなる ベクトル スカラー共に効く :2 倍以上の性能 連続領域 キャッシュに載る ループあたりの計算量増加 do 3* do Y() D()*() (3*-2) do nde(-) nde() 2 (3*-) tem() 3 (3*) enddo enddo Y() Y() AMA()*() Y(3*-2) D(9*-8)*D(9*-7)*2D(9*-6)*3 Y(3*-) D(9*-5)*D(9*-)*2D(9*-3)*3 Y(3*I ) D(9*-2)*D(9*-)*2D(9*I )*3 do nde(-) nde() enddo enddo tem() (3*-2) 2 (3*-) 3 (3*) Y(3*-2) Y(3*-2)AMA(9*-8)*AMA(9*-7)*2 & AMA(9*-6)*3 Y(3*-) Y(3*-)AMA(9*-5)*AMA(9*-)*2 & AMA(9*-3)*3 Y(3*I ) Y(3*I )AMA(9*-2)*AMA(9*-)*2 & AMA(9* )*3
FEM3D-Part 99 三次元要素の定式化 三次元弾性力学方程式 ガラーキン法 要素マトリクス生成 宿題 プログラムの実行 データ構造 プログラムの構成 計算効率について Paraew による可視化
FEM3D-Part マイクロプロセッサの動向 CPU 性能 メモリバンド幅のギャップ http://www.streambench.org/
FEM3D-Part スカラープロセッサ CPU- キャッシュ - メモリの階層構造 FAS CPU Regster Cache 小容量 (MB): 一時置き場高価大きい ( 億以上のトランジスタ ) SLOW Man Memor 大容量 (GB) 廉価
FEM3D-Part 2 ベクトルプロセッサ ベクトルレジスタと高速メモリ ector Processor 単純構造のDOループの並列処理 単純 大規模な演算に適している ector Regster do A() B() C() enddo er FAS Man Memor
FEM3D-Part 3 典型的な挙動 : CG 法 ( 疎行列 FEM) 3..E2 2.5 % of pea GFL LOPS 2..5..5 8 % of pea GFL LOPS.E.E..E.E5.E6.E7 IBM-SP3: DOF: Problem Se 問題サイズが小さい場合はキャッシュの影響のため性能が良い.E-.E.E5.E6.E7 DOF: Problem Se Earth Smulator: 大規模な問題ほどベクトル長が長くなり 性能が高い
FEM3D-Part プロセッサに応じたチューニング メモリ参照の最適化 に尽きる j A(j)
FEM3D-Part 5 プロセッサに応じたチューニング ( 続き ) ベクトルプロセッサ ループ長を大きくとる スカラープロセッサ キャッシュを有効利用 細切れにしたデータを扱う PC クラスタなどでは キャッシュサイズを大きくできない一方で メモリレイテンシ ( 立ち上がりオーバーヘッド ) の減少 メモリバンド幅の増加の傾向 しかしマルチコア化で帳消し 共通事項 メモリアクセスの連続性 局所性 計算順序の変更によって計算結果が変わる可能性について注意すること
FEM3D-Part 6 Scalar スカラープロセッサの代表的なチューニング ループアンローリング ループのオーバーヘッドの削減 ロード ストアの削減 ブロック化 キャッシュミスの削減
FEM3D-Part 7 Scalar ループアンローリング ロード ストアの削減 (/) ループ処理に対する演算の割合が増す do j do A() A() B()*C(j) enddo enddo do j - 2 do A() A() B()*C(j) A() A() B()*C(j) enddo enddo 2K 東大での計算時間.2338E- 3.85938E- 2.6788E- do j -3 do A() A() B()*C(j) A() A() B()*C(j) A() A() B()*C(j2) A() A() B()*C(j3) enddo enddo
FEM3D-Part 8 Scalar ループアンローリング ロード ストアの削減 (2/) ロード : メモリ キャッシュ レジスタ ストア : ロードの逆 ロード ストアが少ないほど効率良い FAS CPU Regster Cache SLOW Man Memor
FEM3D-Part 9 Scalar ループアンローリング ロード ストアの削減 (3/) do j do A() A() B()*C(j) ストア ロード ロードロード enddo enddo A()B()C(j) に対して 各ループでロード ストアが () ()C( j) に対して 各ルプでドストアが発生する
FEM3D-Part Scalar ループアンローリング ロード ストアの削減 (/) do j -3 do A() A() B()*C(j) ロードロードロード A() A() B()*C(j) A() A() B()*C(j2) A() A() B()*C(j3) ストア enddo enddo 同一ループ内で複数回表れている場合には 最初にロード 最後にストアが実施され その間レジスターに保持される 計算順序に注意
FEM3D-Part Scalar ループ入替によるメモリ参照最適化 (/2) YPE-A do do j A(j) A(j) B(j) enddo enddo YPE-B do j do A(j) A(j) B(j) enddo enddo j A(j) FORRA では A(j) のアドレスは A() ) A(2) A(3) A() A(2) A(22) A() A(2) A() のように並んでいる C は逆 :A[][] A[][] A[][2] A[-][] A[-][] A[-][-] この順番にアクセスしないと効率悪い ( ベクトル スカラーに共通 )
FEM3D-Part 2 Scalar ループ入替によるメモリ参照最適化 (2/2) YPE-A do do j A(j) A(j) B(j) enddo enddo YPE-B do j do enddo enddo A(j) A(j) B(j) YPE-A for (j; j<; j){ for (; <; ){ A[][j] A[][j] B[][j]; YPE-B for (; <; ){ for (j; j<; j){ A[][j] A[][j] B[][j]; 2K 東大での計算時間 (FORRA) ### ### 52 A 3.25E-2 B 3.9625E-3 ### ### 2 A 2.3375E- B 7.825E-3 ### ### 536 A 3.76563E- B.5625E-2 2 ### ### 28 A 9.296875E- B 3.25E-2 ### ### 256 A 9.6875E- B.6875E-2 ### ### 372 A 2.523E B 7.2875E-2 ### ### 358 A.92875E B 9.765625E-2 ### ### 96 A 3.8688E B.25E-
FEM3D-Part 3 Scalar ブロック化によるキャッシュミス削減 (/7) do do j A(j) A(j) B(j) enddo enddo for (j;j<;j){ for (;<;){ A[j][] A[j][] B[][j]; 計算の処理の都合でこのような順番で計算せざるを得ない場合
FEM3D-Part Scalar キャッシュの有効利用 キャッシュ 実際は6bte~28bteの細かいキャッシュラインに分かれている キャッシュライン単位でメモリへのリクエストが行われる FAS SLOW CPU Regster Cache Man Memor 小容量 (MB) 高価大きい ( 億以上のトランジスタ ) 大容量 (GB) 廉価 LB(ranslate Looasde Buffer) アドレス変換バッファ 仮想アドレスから実アドレスへの変換機能 LB 用のキャッシュ 通常 28 8 bte 程度 : リンク時に可変 キャッシュ が有効に使われていれば LB ミスも起こりにくい
FEM3D-Part 5 Scalar ブロック化によるキャッシュミス削減 (2/7) AB でメモリアクセスパターンが相反 特に B は効率が悪い j j A[j][] B[][j]
FEM3D-Part 6 Scalar ブロック化によるキャッシュミス削減 (3/7) 例えば キャッシュラインサイズを ワードとすると配列の値は以下のようにキャッシュに転送される j A[j][]
FEM3D-Part 7 Scalar ブロック化によるキャッシュミス削減 (/7) したがって A[][][ ] をアクセスしたら A[][][ ] A[][] ] A[][2] ] A[][3] が A[][9] をアクセスしたらA[][9] A[2][] A[2][] A[2][2] がそれぞれキャッシュ上にあるということになる 2 j 2 2 22 3 A[j][] 9
FEM3D-Part 8 Scalar ブロック化によるキャッシュミス削減 (5/7) したがって 以下のようなブロック型のパターンでアクセスすれば キャッシュを有効利用できる j j A[j][] B[][j]
FEM3D-Part 9 Scalar ブロック化によるキャッシュミス削減 (6/7) で囲んだ部分はキャッシュに載っている j j A[j][] B[][j]
FEM3D-Part 2 Scalar ブロック化によるキャッシュミス削減 (7/7) 2 2 ブロック for (j;j<;j){ for (;<;){ A[j][] A[j][] B[][j]; for (j;j<-;j2){ for (;<-;2){ A[j ][ ] A[j ][ ] B[ ][j ]; A[j ][] A[j ][] B[ ][j]; A[j][ ] A[j][ ] B[][j ]; A[j][] A[j][] B[][j]; 2K 東大での実行時間 (FORRA) ### ### 2 BASIC 2.73375E-2 22.5625E-2 ### ### 536 BASIC 6.25E-2 22 3.55625E-2 ### ### 372 BASIC 2.57825E- 22.8375E- ### ### 358 BASIC 3.7938E- 22 2.325E- ### ### 96 BASIC 8.375E- 22.375E-
FEM3D-Part 2 Scalar チューニング : まとめ スカラープロセッサ 密行列 疎行列の場合はもっと難しい ( 研究途上の課題 ) しかし 基本的な考え方は変わらない メモリアクセスの効率化
FEM3D-Part 22 ブロックとして記憶 (3/3) 計算の安定化 対角成分で割るのではなく 対角ブロックの完全 LU 分解を求めて解く 特に悪条件問題で有効 j j j j
FEM3D-Part 23 三次元要素の定式化 三次元弾性力学方程式 ガラーキン法 要素マトリクス生成 宿題 プログラムの実行 データ構造 プログラムの構成 計算効率について Paraew による可視化
FEM3D-Part 2 Paraew ファイルを開く 図の表示 イメージファイルの保存
25 UCD Format (/3) Unstructured Cell Data 要素の種類キーワード 点 線 三角形 四角形 四面体 角錐 三角柱 六面体 二次要素 線 2 三角形 2 四角形 2 四面体 2 pt lne tr quad tet pr prsm he lne2 tr2 quad2 tet2 角錐 2 pr2 三角柱 2 六面体 2 prsm2 he2
26 UCD Format (2/3) Orgnall for AS mcroas Etenson of the UCD fle s np here are two tpes of formats. Onl old tpe can be read b Paraew.
27 UCD Format (3/3): Old Format ( 全節点数 ) ( 全要素数 ) ( 各節点のデータ数 ) ( 各要素のデータ数 ) ( モデルのデータ数 ) ( 節点番号 ) ( 座標 ) (Y 座標 ) (Z 座標 ) ( 節点番号 2) ( 座標 ) (Y 座標 ) (Z 座標 ) ( 要素番号 ) ( 材料番号 ) ( 要素の種類 ) ( 要素を構成する節点のつながり ) ( 要素番号 2) ( 材料番号 ) ( 要素の種類 ) ( 要素を構成する節点のつながり ) ( 要素のデータ成分数 ) ( 成分 の構成数 ) ( 成分 2の構成数 ) ( 各成分の構成数 ) ( 要素データ成分 のラベル )( 単位 ) ( 要素データ成分 2のラベル )( 単位 ) ( 各要素データ成分のラベル )( 単位 ) ( 要素番号 ) ( 要素データ ) ( 要素データ 2) ( 要素番号 2) ( 要素データ ) ( 要素データ 2) ( 節点のデータ成分数 ) ( 成分 の構成数 ) ( 成分 2 の構成数 ) ( 各成分の構成数 ) ( 節点データ成分 のラベル )( 単位 ) ( 節点データ成分 2 のラベル )( 単位 ) ( 各節点データ成分のラベル )( 単位 ) ( 節点番号 ) ( 節点データ) ( 節点データ2) ( 節点番号 2) ( 節点データ ) ( 節点データ 2)