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

Transcription

1 数値計算の基礎 一次元熱伝導方程式の差分法 による解法 年夏季集中講義 中島研吾 並列計算プログラミング 先端計算機演習 66-49

2 概要 一次元熱伝導方程式と差分法 連立一次方程式の解法 反復法について 共役勾配法 CG 法 前処理について 共役勾配法を使用した一次元熱伝導解析プログラムについて

3 d d 一次元熱伝導方程式 / 支配方程式 : 簡単のため熱伝導率 BF d BF m m 一様体積発熱 BF 断熱 実際は以下のような離散化をしているので注意が必要 断熱となっているのはこの面, しかし温度は計算されない

4 4 BF 一次元熱伝導方程式 / 解析解 BF m 断熱となっているのはこの面, しかし温度は計算されない XX X m X.d, メッシュ数 5, とすると,X m 49.5, の点の X 座標は 49. となる BF.d とすると での温度は :

5 5 念のため 差分について 差分法 :Fte Dfferece ethod マクロな微分 微分係数を数値的に近似する手法 以下のような一次元系を考える -

6 6 直感的 というか安易な定義直感的というか安易な定義 と の中点 における微分係数 と の中点 における微分係数 - d d d / となると微分係数の定義そのもの における二階微分係数 / / d d d d d d

7 7 厳密な定義 :Tyor 展開 / 厳密な定義 :Tyor 展開 / - K!! K!!

8 8 厳密な定義 :Tyor 展開 / 厳密な定義 :Tyor 展開 / - 前進差分 K!! 打ち切り誤差 K!! 打ち切り誤差が のオーダー 次精度 一次精度 後退差分 K!! 打ち切り誤差が K!! 打ち切り誤差が のオーダー 一次精度

9 9 厳密な定義 :Tyor 展開 / 厳密な定義 :Tyor 展開 / - 中央差分, 中心差分 K!! K!!!! 打ち切り誤差が K! 打ち切り誤差が のオーダー 二次精度

10 安易な定義 : 実は二次精度だった安易な定義 : 実は二次精度だった - / / K / / / /! /! / / K / / / /! /! / / K! / 打ち切り誤差が のオーダー / /! オダ 二次精度 二点間の中点で二次精度, それ以外の点では一次精度 ということもできる が均一でない場合も同様のことが起こる

11 一次元熱伝導方程式 / 要素単位の線形方程式 差分法による離散化 d d d d d 差分法による離散化各要素における線形方程式は以下のような形になる / / d d d d 各要素における線形方程式は以下のような形になる N BF BF d d N BF N BF N BF A A A R D,, A A A R D

12 連立一次方程式 :N88 の場合未知数 8, 方程式の数 8 : A A BF D R : A A A BF : A 4 : A D R 5 : A 6 : A A 4 A 5 4 A 6 5 A D D D D A R 4 4 A 5 5 A 6 6 A R R R 4 BF 4 5 BF4 5 6 BF5 6 7 BF6 7 : A 7 6 A 7 7 A 7 8 BF7 D R 8 : A 8 7 A 8 8 BF8 D

13 連立一次方程式 :N88 の場合自分とその周囲のみに非ゼロ成分 各行 つ : 疎行列 AD AR A AD AR A AD AR A4 AD4 AR4 A5 AD5 AR5 A6 AD6 AR6 A7 AD7 AR7 このような 疎な 行列を係数とする連立一次方程式を解く [A]{}{b} {b} A A 8 A8 AD8 8 A A BF, A D D, A R R

14 4 係数行列は疎行列 Sprse p tr が多い 非ゼロ成分のみを記憶する方法 Compressed Row Storge, CRS 後述 密行列と比較して取り扱いが困難

15 5 一次元熱伝導方程式と差分法 連立一次方程式の解法 反復法について 共役勾配法 CG 法 前処理について 共役勾配法を使用した一次元熱伝導解析プログラムについて

16 6 科学技術計算における大規模線形 方程式の解法 多くの科学技術計算は, 最終的に大規模線形方程式 Ab を解くことに帰着される : 線形ソルバー mportt, epesve アプリケーションに応じて様々な手法が提案されている 疎行列 sprse, 密行列 dese 直接法 drect, 反復法 tertve 本講義 演習では主に, 疎行列, 反復法について扱う 密行列 dese グローバルな相互作用あり :BE, スペクトル法,O,D 気液 疎行列 sprse ローカルな相互作用 :FE,FD,D 固, 高速多重極展開付 BE

17 密行列 遠隔領域も含め, 多数領域との相互作用あり境界要素法, スペクトル法,D 法 7

18 疎行列 近接領域のみとの相互作用差分法, 有限要素法 8

19 9 線形ソルバー にどう向き合うか? これについてはあとで

20 直接法 Drect ethod Gss の消去法, 完全 U 分解 逆行列 A - を直接求める 利点 安定, 幅広いアプリケーションに適用可能 Prt Pvotg 疎行列, 密行列いずれにも適用可能 欠点 反復法よりもメモリ, 計算時間を必要とする 密行列の場合,ON の計算量 大規模な計算向けではない ON の記憶容量,ON の計算量

21 U 分解法 : 完全 U 分解法 直接法の一種 逆行列を直接求める手法 逆行列 に相当するものを保存しておけるので, 右辺が変わったときに計算時間を節約できる 逆行列を求める際にF- もとの行列ではであったところに値が入る が生じる

22 U 分解による連立一次方程式の解法 A が 行列のとき A を次式のように表すことを あるいはそのような と U そのものを A の U 分解という あるいは そのような と U そのものを A の U 分解という. O O O U A t t f t A :ower trgr prt of mtr A U:Upper trgr prt of mtr A

23 連立一次方程式の行列表現連立次方程式の行列表現元の連立次方程式の般形 元の連立一次方程式の一般形 b b b b 行列表現行列表現 b b b A b O b A A b

24 U 分解による解法 U 分解による解法 U 分解 エル ユーブンカイ A が 行列のとき A を次式のように表すことを A が 行列のとき A を次式のように表すことを あるいは そのような と U そのものを A の U 分解という. O O O U A t t f t A 4 :ower trgr prt of mtr A U:Upper trgr prt of mtr A

25 U 分解を用いた Ab の解法 A U となる A の U 分解 と U を求める y b の解 yを求める 前進代入 U の解 を求める 後退代入 y この が A b の解となる Q A U y b 5

26 yb の解法 y b の解法 b y b b y y b y O b y b y y b y y b y y b y y y y b y y b y 6 芋づる式に解が求まる : 前進代入

27 Uy の解法 U y の解法 y U y y y O y / / y,, y y,, / y y / y j j 7 芋づる式に解が求まる : 後退代入

28 U 分解の求め方 U 分解の求め方 O O O 8

29 j { の 行目 } {U の j 列目 } j { の 行目 } {U の j 列目 } j j j j j O O j O j j j jj j O O O O 9 O O

30 U 分解の求め方 U 分解の求め方 O O O j :[] の第 行と [U] の第 j 列の内積 [U] の第 行のみ残る j :[] の第 行と,[U] の第 j 列の内積 [U] の第 行のみ残る,,,,,,

31 U 分解の求め方 O :[] の第 行と,[U] の第 列の内積,,, O O,,,,,,,,,

32 U 分解の求め方 O O O j :[] の第 行と,[U] の第 j 列の内積,,,,,,,,,,,,,,,,,

33 U 分解の求め方 O 4 :[] の第 行と,[U] の第 列の内積,,, O O,,,,,,,,,,,,,, 4,, 4,,

34 数値例 A 第 行,,, 4 4 第 列 /, / 4 4 / 第 行 6, 第 列, 4 4

35 数値例 続き A 第 行 , 第 列 第 4 行 第 4 列 行 列 行 列 の順に求める式を作っていく. 5

36 数値例 続き 数値例 続き 結局 A U U 6

37 7 定常 sttory 法 反復法 Itertve ethod 反復計算中, 解ベクトル以外の変数は変化せず SOR,Gss-Sede,Jcob など 概して遅い 非定常 osttory 法 拘束, 最適化条件が加わる Kryov 部分空間 sbspce への写像を基底として使用するため, Kryov 部分空間法とも呼ばれる CGCojgte Grdet: 共役勾配法 BCGSTABB-Cojgte Cojgte Grdet Stbzed GRESGeerzed m Resd

38 8 利点 反復法 Itertve ethod 続き 直接法と比較して, メモリ使用量, 計算量が少ない 並列計算には適している 欠点 収束性が, アプリケーション, 境界条件の影響を受けやすい 前処理 precodtog が重要

39 9 一次元熱伝導方程式と差分法 連立一次方程式の解法 反復法について 共役勾配法 CG 法 前処理について 共役勾配法を使用した一次元熱伝導解析プログラムについて

40 4 代表的な反復法 : 共役勾配法 Cojgte Grdet 法, 略して CG 法 最も代表的な 非定常 反復法 対称正定値行列 Symmetrc Postve Defte:SPD 任意のベクトル {} に対して {} T [A]{}> 全対角成分 >, 全固有値 >, 全部分行列式 >と同値 ガラーキン法 熱伝導, 弾性, ねじり : 本コードの場合も SPD アルゴリズム 最急降下法 Steepest Descet ethod の変種 - α p : 反復解,p : 探索ベクトル,α : 定数 厳密解をyとするとき {-y} T [A]{-y} を最小とするような {} を求める 詳細は参考文献参照 例えば : 森正武 数値解析 第 版 共立出版

41 4 共役勾配法のアルゴリズム Compte r b-[a] for,, z - r - ρ - r - z - f p z ese 行列ベクトル積 ベクトル内積 ベクトル定数倍の加減 ed β - ρ - /ρ - p z - β - edf q [A]p α ρ - /p q - α p p - r r - - α q chec covergece r : ベクトル α : スカラー

42 4 共役勾配法のアルゴリズム Compte r b-[a] for,, z - r - ρ - r - z - f p z ese β - ρ - /ρ - p z - β - edf q [A]p α ρ - /p q - α p r r - - α q chec covergece r ed p - 行列ベクトル積 ベクトル内積 ベクトル定数倍の加減 : ベクトル α : スカラー

43 4 共役勾配法のアルゴリズム Compte r b-[a] for,, z - r - ρ r - z - - f p z ese β - ρ - /ρ - p z - β - edf q [A]p α ρ - /p q - α p r r - - α q chec covergece r ed p - 行列ベクトル積 ベクトル内積 ベクトル定数倍の加減 : ベクトル α : スカラー

44 44 共役勾配法のアルゴリズム Compte r b-[a] for,, z - r - ρ - r - z - f p z ese β - ρ - /ρ - p z - β - edf q [A]p α ρ - /p q - α p r r - - α q chec covergece r ed p - 行列ベクトル積 ベクトル内積 ベクトル定数倍の加減 : ベクトル α : スカラー

45 45 共役勾配法のアルゴリズム /5 yを厳密解 Ayb y とするとき, 下式を最小にするを求める : T y [ A] y T y [ A] y, A y, A, Ay y, Ay, A, Ay y, Ay, A, b y, b 定数 従って, 下記 f を最小にする を求めればよい : f,, A b f, h f h, A b h Ah 任意のベクトル h

46 46 b A f,, Ah h b A h f h f,, 任意のベクトル h Ah h b A h f h f,, b h h A h h f,, b h b Ah h A h,,,, b h b Ah h Ah A h A,,,,,, Ah h b A h f Ah h b h A h b A,,,,, Ah h b A h f,,

47 47 共役勾配法のアルゴリズム /5 共役勾配法のアルゴリズム /5 CG 法は任意の から始めて,f の最小値を逐次探索する p α f 今, 番目の近似値 と探索方向 p が決まったとすると : p α f を最小にするためには :,, f A b p Ap p p f α α α f を最小にするためには :,, f p p p p f,, A r p A A b p p f α α,, Ap p Ap p α A b r は第 近似に対する残差 A b r は第 近似に対する残差

48 48 共役勾配法のアルゴリズム /5 残差 r も以下の式によって計算できる : r r α Ap 探索方向を以下の漸化式によって求める : p r β p, r p 本当のところは下記のように 回目に厳密解 y が求まれば良いのであるが, 解がわかっていない場合は困難 y α p

49 49 共役勾配法のアルゴリズム 4/5 共役勾配法のアルゴリズム 4/5 ところで, 下式のような都合の良い直交関係がある :, 下式ような都合良直交関係ある, y Ap [ ],,, Ap A b p p A b p A b p A Ay p y Ap α α [ ],,,,, Ap p r p Ap r p Ap A b p p A b p α α α α,, Ap p r p Qα 従って以下が成立する A A A α, p p,,, Ap p p Ap y Ap α

50 5 共役勾配法のアルゴリズム 5/5 p, Ap r β,,, p Ap r Ap β p Ap r, Ap β p, Ap p, Ap 勾配ベクトル p が A に関して共役である 勾配ベクトル p, 残差ベクトル r についても以下の関係が成立 証明略 : 帰納法 j j p, Ap j, r, r j p, r r, r N 次元空間で互いに直交で一次独立な残差ベクトル r は N 個しか存在しない, 従って共役勾配法は未知数が N 個のときに N 回以内に収束する 実際は丸め誤差の影響がある

51 5 α,β 実際は α,β はもうちょっと簡単な形に変形できる : α Q p, b A p, r r, r p, Ap p, Ap p, Ap p, r r, r r, Ap r, r p, Ap r, r r, r r r, r r, Ap β Q α α

52 5 共役勾配法のアルゴリズム Compte r b-[a] for,, z - r - ρ - r - z - f p z ese 行列ベクトル積 ベクトル内積 ベクトル定数倍の加減 ed β - ρ - /ρ - p z - β - edf q [A]p α ρ - /p q - α p p - r r - - α q chec covergece r : ベクトル α : スカラー

53 5 前処理 precodtog とは? 反復法の収束は係数行列の固有値分布に依存 固有値分布が少なく, かつに近いほど収束が早い 単位行列 条件数 codto mber 対称正定 最大最小固有値の比 条件数が に近いほど収束しやすい もとの係数行列 [A] に良く似た前処理行列 [] を適用することによって固有値分布を改善する 前処理行列 [] によって元の方程式 [A]{}{b} {b} を [A ]{ }{b } へと変換する ここで [A ][] - [A], {b }[] - {b} である [A ][] - [A] が単位行列に近ければ良い, ということになる 前処理 は密行列, 疎行列ともに使用するが, 普通は疎行列前処理 は密行列, 疎行列ともに使用するが, 普通は疎行列を対象にすることが多い

54 前処理付共役勾配法 Precodtoed Cojgte Grdet ethod PCG 54 Compte r b-[a] for,, sove []z - r - ρ - r - z - f p z ese β - ρ - /ρ - p z - β - edf q [A]p α ρ - /p q - α p r r - - α q chec covergece r ed z - 実際にやるべき計算は : { z} [ ] { r } 近似逆行列 の計算が必要 : [ ] [ ] A, [ ] [ A] 究極の前処理 : 本当の逆行列 [ ] [ A], [ ] [ A ] 対角スケーリング : 簡単 弱い [ ] [ ] D, [ ] [ D]

55 55 IU, IC 最もよく使用されている前処理 疎行列用 不完全 U 分解 Icompete U Fctorzto 不完全コレスキー分解 Icompete Choesy Fctorzto 対称行列 不完全な直接法 もとの行列が疎でも, 逆行列は疎とは限らない f- もとの行列と同じ非ゼロパターン f- 無し を持っているのが IU,IC

56 56 IU, IC 最もよく使用されている前処理 疎行列用 Icompete U Fctorzto Icompete Choesy Fctorzto 対称行列 IU : eep o-zero ptter of the org coeffcet mtr do, do, - f, NoZeroA the : / edf do j, f,j NoZeroA the j : j - * j edf eddo eddo eddo 不完全な直接法 もとの行列が疎でも, 逆行列は疎とは限らない f- もとの行列と同じ非ゼロパターン f- 無し を持っているのがIU,IC

57 U 分解法 : 完全 U 分解法 直接法の一種 逆行列を直接求める手法 逆行列 に相当するものを保存しておけるので, 右辺が変わったときに計算時間を節約できる 逆行列を求める際にF- もとの行列ではであったところに値が入る が生じる 不完全 U 分解法 F- の発生を制限して, 前処理に使う手法 不完全な逆行列, 少し弱い直接法 57

58 58 実例 : 差分法による熱伝導等 5 点差分

59 59 実例 : 差分法による熱伝導等 5 点差分

60 6 実例 : 係数マトリクス X

61 6 実例 : 解

62 6 ECCS8ログイン, ディレクトリ作成 ディレクトリ作成 >$ cd Docmets この下に作る >$ mdr EPSsmmer 好きな名前でよい >$ cd EPSsmmer このディレクトリを本講義では <$E-EPS> と呼ぶ基本的にファイル類はこのディレクトリにコピー, 解凍する この下に課題番号に応じて S, S, S-ref などのディレクトリを作る <$S> <$E-EPS>/S <$S> <$E-EPS>/S EPS>/S TK <$T-EPS>

63 6 ファイルコピー 一次元熱伝導方程式 >$ cd <$E-EPS> >$ cp /home/sego/docmets/css/epssmmer/tro.tr. >$ tr vf tro.tr >$ cd tro

64 64 完全 U 分解, 不完全 U 分解 >$ cd <$E-EPS>/tro.f.f 完全 U 分解によるプログラム FORTRAN 不完全 U 分解 f- 無し によるプログラム FORTRAN.f をコンパイル, リンクしたもの.f をコンパイル, リンクしたもの C ユーザーも FORTRAN でお願いします

65 65 完全 U 分解したマトリクス./ とタイプ もとのマトリクス U 分解したマトリクス [][U] 同時に表示 [] 対角成分 省略 f-が生じている もともとだった成分が非ゼロになっている

66 66 不完全 U 分解したマトリクス f- 無し./ とタイプ 不完全 U 分解したマトリクス f- 無し [][U] 同時に表示 [] 対角成分 省略 完全 U 分解したマトリクス [][U] 同時に表示 [] 対角成分 省略 f-が生じている もともと だった成分が非ゼロになっている

67 67 解の比較 : ちょっと違う 不完全 U 分解 完全 U 分解

68 68 IU, IC 前処理 F- を全く考慮しない 不完全な 分解 記憶容量, 計算量削減 これを解くと 不完全な 解が得られるが, 本来の解とそれほどずれているわけではない 問題に依存する

69 69 大規模線形ソルバーの動向 反復法がより広く使用されるようになりつつある コアを超えるような大規模並列システムでは直接法は並列性能が並出ない : 逆にそれより小さければ直接法でもOKということになる 密行列も反復法で解くような試みがなされている 密行列を使わないで済ませられるようなアルゴリズムの開発 高速多重極展開 Fst tpoe 遠方からの効果をクラスタリング, あるいは無視 密行列 メモリースケーラブルではない 前処理付き反復法 precodtoed tertve sovers 安定した前処理の必要性 安定した前処理は概して 並列化 が困難 詳しくはCEで

70 7 一次元熱伝導方程式と差分法 連立一次方程式の解法 反復法について 共役勾配法 CG 法 前処理について 共役勾配法を使用した一次元熱伝導解析プログラムについて

71 7 d d 再び一次元熱伝導方程式 / 支配方程式 : 熱伝導率 一様 BF d BF m m 一様体積発熱 BF 断熱 以下のような離散化 要素中心で従属変数を定義 をしているので注意が必要 断熱となっているのはこの面, しかし温度は計算されない

72 7 BF 一次元熱伝導方程式 / 解析解 BF m 断熱となっているのはこの面, しかし温度は計算されない XX X m T X.d, メッシュ数 5, とすると,X m 49.5, の点の X 座標は 49. となる BF.d とすると での温度は :

73 7 一次元熱伝導方程式 / 連立一次方程式 差分法による離散化 / / d d d d d 各要素における線形方程式は以下のような形になる / / d 各要素における線形方程式は以下のような形になる これを共役勾配法 Cojgte Grdet 法 で解く N BF BF d d N BF N BF N BF A A A R D N BF,, A A A R D

74 74 係数行列の格納形式 各要素における線形方程式 A A N A D A R BF N, A D, A R AD AR A AD AR A AD AR A4 AD4 AR4 A5 AD5 AR5 A6 AD6 AR6 7 A7 AD7 AR7 より一般的な形で格納すると INDEX DIAG PHI 対角成分 INDEX [ AAT PHI ITE ] RHS,, K, N 非対角成分 8 A8 AD8

75 75 N8 の場合 : 連立一次方程式自分とその周囲のみに非ゼロ成分 : 疎行列 A D A R BF A A BF D R A A D D A R R A A D A R BF A A A BF BF A D R A A 4 BF D R 4 5 A 4 A D D 4 A R R 4 A 5 A D 5 A R BF4 BF5 A 4 A 4 4 A 4 5 BF4 D R A 5 4 A 5 5 A 5 6 BF5 D R 6 A 6 A D D 6 A R R 6 6 BF6 A 6 5 A 6 6 A 6 7 BF6 D R 7 A 7 A D 7 A R 7 7 BF7 A 7 6 A 7 7 A 7 8 BF7 D R 8 D 8 BF8 A 8 7 A D 8 8 BF8 A 8 A D 8 INDEX DIAG PHI 対角成分 [ AAT PHI ITE ] RHS,, K, N INDEX 非対角成分

76 76 未知数 :N 個とすると 密行列 疎行列 密行列 :N 個の成分を記憶する必要あり 疎行列 : 非ゼロの部分だけ記憶すれば, 一次元問題の場合 N 個の成分で済む N 6 としたときの必要記憶容量 密行列 TB 疎行列 後掲の INDEX,ITE も含む B

77 77 Compressed Row Storge CRS 一般的な形

78 78 係数行列の格納形式 /8 非ゼロ成分のみを格納, 疎行列向け方法 Compressed Row Storge CRS DIAG PHI INDEX INDEX [ AAT PHI ITE ] RHS,, K, N DIAG 対角成分 実数,,N INDEX 非対角成分数に関する一次元配列 整数,,N ITE 非対角成分の要素 列 番号 整数,, INDEXN AAT 非対角成分 実数,,, INDEXN A D A R A A D A R A A D A R A 4 A D 4 A R 4 A 5 A D 5 A R 5 A 6 A D 6 A R 6 A 7 A D 7 A R 7 A 8 A D 8

79 79 Compressed Row Storge CRS NODE_tot 8 対角成分 D. D.6 D 5.7 D4 9.8 D5.5 D6.4 D7. D8 5.

80 8 Compressed Row Storge CRS

81 Compressed Row Storge CRS 非対角 de 成分数 4.4. de de de de de de de de8 5 NPU den de-~de 番目が 行目の非対角成分

82 Compressed Row Storge CRS 非対角成分数 de 4.4., 5, de , 4,44 6,55 8,66 4 de ,7 7,8 de ,9 5, 6, de ,,,4 7,5 4 de ,6 7,7 de ,8,9 6, 8, 4 de de8 5 NPU 5,, 4,4 6,5 den de-~de 番目が 行目の非対角成分

83 8 Compressed Row Storge CRS , 5, , 4,44 6,55 8, ,7 7, ,9 5, 6, ,,,4 7, ,6 7, ,8,9 6, 8, ,, 4,4 6,5 例 : tem 7 5, AAT 7.5 tem9, AAT9.5

84 84 Compressed Row Storge CRS. 4.4., 5, D 対角成分 実数,,N de 非対角成分数に関する一次元配列, 4,44 6,55 8,66 通し番号 整数,,N tem 非対角成分の要素 列 番号 5,7 7,8 整数,, den AAT 非対角成分 4,9 5, 6, 実数,, den ,,,4 7,5 {Y} [A]{X} do, N 6,6 7,7 Y D*X do de-, de 7,8,9 6, 8, Y Y AAT*Xtem 8 8,, 4,4 6,5 eddo eddo

85 85 行列ベクトル積 : 密行列 とても簡単..., N, N y N,, N y N N N N,,, N, N N yn N, N,... N, N N, N N yn {Y} [A]{X} do j, N Yj.d do, N Yj Yj A,j*X eddo eddo

86 86 行列ベクトル積への適用非ゼロ成分のみを格納, 疎行列向け方法 DIAG 対角成分 実数,,N DIAG PHI INDEX INDEX 非対角成分数に関する一次元配列 AAT PHI ITE RHS,, K, 整数,,N INDEX ITE 非対角成分の要素 列 番号 整数,, INDEXN AAT 非対角成分 A D A R 実数,, INDEXN A A D A R {Y} [A]{X} [ ] N 4 A A D A R A 4 A D 4 A R 4 do, N Y D*X 5 A 5 A D 5 A R 5 do INDEX-, INDEX 6 A 6 A D 6 A R 6 Y Y AAT*XITE 7 A eddo 7 A D 7 A R 7 eddo 8 A 8 A D 8

87 87 DIAG PHI INDEX 係数行列の格納形式 /8 [ AAT PHI ITE ] 対角成分 :DIAG INDEX RHS, DIAG, K, N DIAG A D DIAG A D DIAG A D DIAG8 A D 8 対角成分 実数,,N A D A R A A D A R A A D A R A 4 A D 4 A R 4 A 5 A D 5 A R 5 A 6 A D 6 A R 6 A 7 A D 7 A R 7 8 A 8 A D 8

88 88 係数行列の格納形式 /8 INDEX, AAT の関係 A D A R AD A A D A R AD A A D A R 4 AD 5 4 A 4 A D 4 A R AD4 7 5 A 5 A D 5 A R AD A 6 A D 6 A R 6 A 7 A D 7 A R AD6 AD7 8 A 8 A D AD8 INDEX, INDEX5, 5 ITE4, AAT4 A ITE54, AAT5 A R

89 89 係数行列の格納形式 /8 INDEX, AAT の関係 A D A R AD A A D A R AD A A D A R 4 AD 5 4 A 4 A D 4 A R AD4 7 5 A 5 A D 5 A R AD A 6 A D 6 A R 6 A 7 A D 7 A R AD6 AD7 8 A 8 A D AD8 INDEX, INDEX5, 5 ITE4, AAT4 A ITE54, AAT5 A R

90 9 係数行列の格納形式 4/8 非対角成分 : DIAG PHI INDEX INDEX RHS, [ AAT PHI ITE ], K, N INDEX 非対角成分数に関する一次元配列 整数,,N ITE 非対角成分の要素 列 番号 整数,, INDEXN AAT 非対角成分 実数,, INDEXN INDEX INDEX ITE, AAT A R AD 4 AD 5 6 AD4 7 8 AD5 9 AD6 AD7 AD 4 AD8

91 9 係数行列の格納形式 5/8 非対角成分 : DIAG PHI INDEX INDEX RHS, [ AAT PHI ITE ], K, N INDEX 非対角成分数に関する一次元配列 整数,,N ITE 非対角成分の要素 列 番号 整数,, INDEXN AAT 非対角成分 実数,, INDEXN INDEX INDEX ITE, AAT A ITE, AAT A R AD AD AD 5 6 AD4 7 8 AD5 9 AD6 AD7 4 AD8

92 9 係数行列の格納形式 6/8 非対角成分 : DIAG PHI INDEX INDEX RHS, [ AAT PHI ITE ], K, N INDEX 非対角成分数に関する一次元配列 整数,,N ITE 非対角成分の要素 列 番号 整数,, INDEXN AAT 非対角成分 実数,, INDEXN INDEX INDEX5 ITE4, AAT4 A ITE54, AAT5 A R AD AD 4 AD 5 6 AD4 7 8 AD5 9 AD6 AD7 4 AD8

93 9 係数行列の格納形式 7/8 非対角成分 :7 DIAG PHI INDEX INDEX RHS, [ AAT PHI ITE ], K, N INDEX 非対角成分数に関する一次元配列 整数,,N ITE 非対角成分の要素 列 番号 整数,, INDEXN AAT 非対角成分 実数,, INDEXN INDEX6 INDEX7 ITE6, AAT A 7 ITE8, AAT A R AD AD 4 AD 5 6 AD4 7 8 AD5 9 AD6 AD7 4 AD8

94 94 係数行列の格納形式 8/8 非対角成分 :8 DIAG PHI INDEX INDEX RHS, [ AAT PHI ITE ], K, N INDEX 非対角成分数に関する一次元配列 整数,,N ITE 非対角成分の要素 列 番号 整数,, INDEXN AAT 非対角成分 実数,, INDEXN INDEX7 INDEX84 ITE47, AAT4 A AD AD 4 AD AD AD5 9 6 AD6 7 AD7 8 4 AD8

95 95 プログラムのコンパイル, 実行法 プログラムのありか $ cd <$E-EPS>/tro $ s het_cg.f $ cd <$E-EPS>/tro $ s het_cg.c コンパイル, 実行法 $ cd <$E-EPS>/tro $ cd <$E-9S>/tro $ g95 het_cg.f o cg $ cc het_cg.c o cg $./cg $./cg pt.dt N メッシュ数.d.d d, BF メッシュ幅, 体積発熱量 5 ITERm 最大反復回数.e-7 EPS 打切誤差

96 96 CG 法による一次元熱伝導方程式 計算プログラム 前処理付き共役勾配法 Precodtoed Cojgte Grdet ethod 対称正定行列用 対角スケーリング Dgo Scg [A] の対角成分のみからなる行列を前処理行列 [] とする 簡単 な問題にしか適用できない 点ヤコビ Pot Jcob 前処理とも呼ばれる Compte r b-[a] for,, sove []z - r - ρ - r - z - f p z ese β - ρ - /ρ - p z - β - p edf q [A]p α ρ - /p q - α p r r - - α q chec covergece r ed p -

97 97 対角スケーリング, 点ヤコビ前処理 前処理行列として, もとの行列の対角成分のみを取り出した行列を前処理行列 [] とする 対角スケーリング, 点ヤコビ pot-jcob 前処理 D... D DN... D N [ ] sove []z - r - という場合に逆行列を簡単に求めることができる

98 98 CG 法による一次元熱伝導方程式 プログラム概要 /7!C!C D Posso Eqto Sover by!c CG Cojgte Grdet ethod!c!c d/ddphi/d BF!C PHI@!C progrm CG_po mpct REA*8 A-H,O-Z teger :: N, ITERm teger :: R, Z, P, Q, DD pt.dt N メッシュ数 d, BF メッシュ幅, 体積発熱量 ITERm 最大反復回数 EPS 打切誤差 red8 :: d, OEGA, RESID, dphi, dphim, BF, EPS red8, dmeso:, octbe :: PHI, RHS red8, dmeso:, octbe :: DIAG, AAT red8, dmeso:,:, octbe :: W teger, dmeso:, octbe :: INDEX, ITE!C!C-- INIT. ope, fe'pt.dt', stts'ow' red,* N red,* dx, BF red,* ITERm red,* EPS cose

99 99 CG 法による一次元熱伝導方程式 プログラム概要 / octe PHIN, DIAGN, AAT*N-, RHSN octe INDEX:N, ITE*N-, WN,4 PHI.d AAT.d/dX DIAG -.d/dx RHS -BF * dx 4 d d d V BF V 非対角成分の数は ただし,,N N のときは したがって, 非対角成分の総数は *N m A D A R A A D A R A A D A R A 4 A D 4 A R 4 A 5 A D 5 A R 5 A 6 A D 6 A R 6 A 7 A D 7 A R 7 A 8 A D 8

100 CG 法による一次元熱伝導方程式 プログラム概要 / octe PHIN, DIAGN, AAT*N-, RHSN octe INDEX:N, ITE*N-, WN,4 PHI.d AAT.d/dX DIAG -.d/dx RHS -BF * dx 4 V BF V BF 非対角成分の数は ただし,,N N のときは したがって, 非対角成分の総数は *N A D A R A A D A R A A D A R A 4 A D 4 A R 4 A 5 A D 5 A R 5 A 6 A D 6 A R 6 A 7 A D 7 A R 7 A 8 A D 8

101 CG 法による一次元熱伝導方程式 プログラム概要 / octe PHIN, DIAGN, AAT*N-, RHSN octe INDEX:N, ITE*N-, WN,4 PHI.d AAT.d/dX DIAG -.d/dx RHS -BF * dx 4 BF A A D A R 非対角成分の数は ただし,,N N のときは したがって, 非対角成分の総数は *N A D A R A A D A R A A D A R A 4 A D 4 A R 4 BF RHS A 5 A D 5 A R 5 A 6 A D 6 A R 6 A 7 A D 7 A R 7 A 8 A D 8

102 CG 法による一次元熱伝導方程式 プログラム概要 /7 INDEX!C!C-- CONNECTIVITY INDEX INDEX INDEXN A D A R A A D A R do, N INDEX INDEX INDEX- 4 eddo 5 do, N 非対角成分の数は 6 js INDEX- f.eq. the ただし,,Nのときは 7 ITEjS 8 AATjS.d DIAG.d このルールに従って ルに従って INDEX RHS.d ese f 生成 & &.eq.n the ITEjS - DIAG -.d/dx ese ITEjS - ITEjS f -.eq. the AATjS.d edf edf eddo A A D A R A 4 A D 4 A R 4 A 5 A D 5 A R 5 A 6 A D 6 A R 6 A 7 A D 7 A R 7 A 8 A D 8

103 CG 法による一次元熱伝導方程式 プログラム概要 /7 INDEX !C!C-- CONNECTIVITY AD INDEX INDEX INDEXN do, N INDEX INDEX INDEX- eddo 4 5 AD AD 4 AD 5 6 do, N 7 js INDEX- f.eq. the 8 ITEjS AATjS.d DIAG.d INDEX RHS.d INDEX ese f ITE, AAT AR & &.eq.n the INDEX ITE, AAT A ITEjS - ITE, AAT AR DIAG -.d/dx INDEX5 ese ITE4, AAT4 A ITEjS - ITE54, AAT5 AR ITEjS f -.eq. the edf edf eddo AATjS.d 6 AD4 7 8 AD5 9 AD6 AD7 4 AD8 INDEXN- *N- ITE*N-4 N-, AAT*N-4 AN- ITE*N- N, AAT*N- ARN- INDEXN *N- ITE*N- N-, AATN*N- AN

104 CG 法による一次元熱伝導方程式 4 プログラム概要 /7 INDEX!C!C-- CONNECTIVITY INDEX INDEX INDEXN A D A R A A D A R do, N INDEX INDEX INDEX- eddo 4 A A D A R 5 do, N 6 js INDEX- f.eq. the 7 ITEjS 8 AATjS.d DIAG.d 境界条件 : 後述 RHS.d ese f & &.eq.n the ITEjS - DIAG -.d/dx 境界条件 N ese ITEjS - ITEjS それ以外 f -.eq. the edf edf eddo AATjS.d A 4 A D 4 A R 4 A 5 A D 5 A R 5 A 6 A D 6 A R 6 A 7 A D 7 A R 7 A 8 A D 8

105 CG 法による一次元熱伝導方程式 5 プログラム概要 /7 INDEX !C A D A R!C-- CONNECTIVITY D R INDEX INDEX INDEXN A A D A R A A D A R do, N INDEX INDEX INDEX- eddo 4 5 A 4 A D 4 A R 4 6 do, N 7 js INDEX- f.eq. the 8 ITEjS AATjS.d DIAG.d 境界条件 : 後述 RHS.d ese f & &.eq.n the ITEjS - DIAG -.d/dx 境界条件 N ese ITEjS - ITEjS それ以外 f -.eq. the edf edf eddo AATjS.d 固定境界条件が指定されている点を非対角成分として持っている場合非対角成分をゼロクリアする A 5 A D 5 A R 5 A 6 A D 6 A R 6 A 7 A D 7 A R 7 A 8 A D 8

106 6 境界条件の処理 :N N js INDEX- ITEjS - AATjS.d/d DIAG -.d/dx デフォルト値 d N m d N N N BF N N BF BF N N N DIAG AATjS RHS N N 境界面で断熱条件が成立するためには, N N を満たすような仮想的な要素があると都合が良い

107 7 境界条件の処理 : /4 ITEjS AATjS.d DIAG.d RHS.d DIAG RHS これでは係数行列が対称とならないため共役勾配法が使えない

108 8 境界条件の処理 :/4 ITEjS AATjS.d DIAG.d RHS.d DIAG RHS これでは係数行列が対称とならないため共役勾配法が使えない

109 9 における方程式 境界条件の処理 :/4 RHS ここで は既知の値

110 における方程式 境界条件の処理 :4/4 RHS ここで は既知の値 RHS このケースの場合は, であるため, 右辺 RHS の修正は不要である

111 固定境界条件の処理 : 消去 の点でPHIの値に固定されている場合 do, N f.eq. the DIAG.d RHS PHI do j INDEX-, INDEX AATj.d eddo edf edf ゼロクリア do, N f.e. the do j INDEX-, INDEX jj ITEj f jj.eq. the RHS RHS - AATj*PHI AATj.d edf eddo eddo eddo

112 固定境界条件の処理 : 消去 の点でPHIの値に固定されている場合 do, N f.eq. the DIAG.d RHS PHI do j INDEX-, INDEX AATj.d eddo edf edf do, N f.e. the do j INDEX-, INDEX jj ITEj f jj.eq. the RHS RHS - AATj*PHI AATj.d edf eddo eddo eddo ゼロクリア 今回の問題では PHIXjj に相当するものがだったのでこの右辺の処理は不要であった

113 固定境界条件の処理 : 具体例 点 で固定境界条件を適用しているものとする ~ , とすると 8 DIAG AAT AAT 4 DIAG 4 AAT AAT RHS AAT 47 7 RHS 4

114 4 固定境界条件の処理 : 具体例 点 で固定境界条件を適用しているものとする ~ , とすると 8 DIAG AAT AAT 4 DIAG 4 AAT AAT RHS AAT RHS 4 4 DIAG AA66 RHS AAT ~ 4 ~ 4

115 5 CG 法による一次元熱伝導方程式!C!C !C CG tertos!c !C R Z Q P DD 4!C!C-- {r} {b} - [A]{} プログラム概要 4/7 do, N 対角成分の逆数 前処理用 W,DD.D / DIAG その都度, 除算をすると効率が eddo 悪いため, 予め配列に格納 do, N W,R DIAG*PHI do j INDEX-, INDEX eddo eddo W,R W,R AATj*PHIITEj BNR.D do, N BNR BNR RHS ** W,R RHS - W,R eddo!c********************************************************************

116 6 CG 法による一次元熱伝導方程式!C!C !C CG tertos!c !C R Z Q P DD 4 do, N W,DD.D / DIAG eddo プログラム概要 4/7 Compte r b-[a] for,, sove []z - r - ρ - r - z - f p z ese β - ρ - /ρ - p z - β - edf q [A]p α ρ - /p q - α p W, W,R {r} W, W,Z {z} W, W,Q {q} α p W, W,P {p} W,4 W,DD /DIAG ed p - r r - - α q chec covergece r

117 7 CG 法による一次元熱伝導方程式!C!C !C CG tertos!c !C R Z Q P DD 4 do, N W,DD.D / DIAG eddo!c!c-- {r} {b} - [A]{} プログラム概要 4/7 do, N W,R DIAG*PHI do j INDEX-, INDEX eddo eddo W,R W,R AATj*PHIITEj BNR.D do, N BNR BNR RHS ** W,R RHS - W,R eddo Compte r b-[a] for,, sove []z - r - ρ - r - z - f p z ese β - ρ - /ρ - p z - β - edf q [A]p α ρ - /p q - α p r r - - α q chec covergece r ed BNR: 収束判定用 p -

118 8 CG 法による一次元熱伝導方程式 プログラム概要 5/7 do ter, ITERm!C!C-- {z} [v]{r} do, N W,Z W,DD * W,R eddo!c!c-- RHO {r}{z} RHO.d do, N RHO RHO W,R*W,Z eddo!c!c-- {p} {z} f ITER!C BETA RHO / RHO otherwse f ter.eq. the do, N W,P W,Z eddo ese BETA RHO / RHO do, N W,P W,Z BETA*W,P eddo edf Compte r b-[a] for,, sove []z - r - ρ - r - z - f p z ese β - ρ - /ρ - p z - β - edf q [A]p α ρ - /p q - α p r r - - α q chec covergece r ed p -

119 9 CG 法による一次元熱伝導方程式 プログラム概要 5/7 do ter, ITERm!C!C-- {z} [v]{r} do, N W,Z W,DD * W,R eddo!c!c-- RHO {r}{z} RHO.d do, N RHO RHO W,R*W,Z eddo!c!c-- {p} {z} f ITER!C BETA RHO / RHO otherwse f ter.eq. the do, N W,P W,Z eddo ese BETA RHO / RHO do, N W,P W,Z BETA*W,P eddo edf Compte r b-[a] for,, sove []z - r - ρ - r - z - f p z ese β - ρ - /ρ - p z - β - edf q [A]p α ρ - /p q - α p r r - - α q chec covergece r ed p -

120 CG 法による一次元熱伝導方程式 プログラム概要 5/7 do ter, ITERm!C!C-- {z} [v]{r} do, N W,Z W,DD * W,R eddo!c!c-- RHO {r}{z} RHO.d do, N RHO RHO W,R*W,Z eddo!c!c-- {p} {z} f ITER!C BETA RHO / RHO otherwse f ter.eq. the do, N W,P W,Z eddo ese BETA RHO / RHO do, N W,P W,Z BETA*W,P eddo edf Compte r b-[a] for,, sove []z - r - ρ - r - z - f p z ese β - ρ - /ρ - p z - β - edf q [A]p α ρ - /p q - α p r r - - α q chec covergece r ed p -

121 CG 法による一次元熱伝導方程式 プログラム概要 6/7!C!C-- {q} [A]{p} do, N W,Q DIAG*W,P do j INDEX-, INDEX eddo eddo W,Q W,Q AATj*WITEj WITEj,P!C!C-- APHA RHO / {p}{q} } C.d do, N C C W,P*W,Q eddo APHA RHO / C!C!C-- {} {} APHA*{p}!C {r} {r} - APHA*{q} do, N PHI PHI APHA * W,P W,R W,R - APHA * W,Q eddo Compte r b-[a] for,, sove []z - r - ρ - r - z - f p z ese β - ρ - /ρ - p z - β - edf q [A]p α ρ - /p q - α p r r - - α q chec covergece r ed p -

122 CG 法による一次元熱伝導方程式 プログラム概要 6/7!C!C-- {q} [A]{p} do, N W,Q DIAG*W,P do j INDEX-, INDEX eddo eddo W,Q W,Q AATj*WITEj WITEj,P!C!C-- APHA RHO / {p}{q} } C.d do, N C C W,P*W,Q eddo APHA RHO / C!C!C-- {} {} APHA*{p}!C {r} {r} - APHA*{q} do, N PHI PHI APHA * W,P W,R W,R - APHA * W,Q eddo Compte r b-[a] for,, sove []z - r - ρ - r - z - f p z ese β - ρ - /ρ - p z - β - edf q [A]p α ρ - /p q - α p r r - - α q chec covergece r ed p -

123 CG 法による一次元熱伝導方程式 プログラム概要 6/7!C!C-- {q} [A]{p} do, N W,Q DIAG*W,P do j INDEX-, INDEX eddo eddo W,Q W,Q AATj*WITEj WITEj,P!C!C-- APHA RHO / {p}{q} } C.d do, N C C W,P*W,Q eddo APHA RHO / C!C!C-- {} {} APHA*{p}!C {r} {r} - APHA*{q} do, N PHI PHI APHA * W,P W,R W,R - APHA * W,Q eddo Compte r b-[a] for,, sove []z - r - ρ - r - z - f p z ese β - ρ - /ρ - p z - β - edf q [A]p α ρ - /p q - α p r r - - α q chec covergece r ed p -

124 4 CG 法による一次元熱伝導方程式 DNR. do, N DNR DNR W,R** eddo RESID dsqrtdnr/bnr wrte *, ter, RESID formt 5, pe6.6 プログラム概要 7/7 f RESID.e.EPS goto 9 RHO RHO Compte r b-[a] for,, sove []z - r - ρ - r - z - f p z ese β - ρ - /ρ - p z - β - edf q [A]p α ρ - /p q eddo!c******************************************************************** IER 9 cote!c!c!c-- OUTPUT do, N wrte *,'8, pe6.6', PHI eddo ed progrm CG_po ed p - - α p r r - - α q chec covergece r

125 5 CG 法による一次元熱伝導方程式 DNR. do, N DNR DNR W,R** eddo RESID dsqrtdnr/bnr wrte r b-[a] *, ter, RESID formt 5, pe6.6 DNR r プログラム概要 7/7 f BNR b RESID.e.EPS goto 9 RHO RHO BNR: 無次元化に使用 Compte r b-[a] for,, sove []z - r - ρ - r - z - f p z ese β - ρ - /ρ - p z - β - edf q [A]p α ρ - /p q eddo RESID r / b!c******************************************************************** IER 9 cote!c!c!c-- OUTPUT do, N wrte *,'8, pe6.6', PHI eddo ed progrm CG_po ed p - - α p r r - - α q chec covergece r

126 6 反復法の収束判定 解を得られた という判定 解の推定値 適切な条件のもとで, のプロセスを繰り返すことによって, は正しい解に収束していく [A]{}{b} という方程式を解いているので, 残差ノルム b-a ~A となれば収束したとみなすことができる 通常は b で無次元化した 残差ノルム が予め設定した値 ε より小さくなった場合に収束したと みなす ε の値は要求される精度によって異なる b A b < ε b A b j j j, b j b j