MPIによる並列アプリケーション
開発法入門（I）
2011年5月19日・20日
2011年5月19日
20日
中島研吾
東京大学情報基盤センター
T2Kオープンスパコン（東大）並列プログラミング講習会

T2K-FVM-01

1

本講習の目的• 1CPU用シリアルアプリケーションの並列化（Flat用シリアルア リケ ション 並列化（
MPI））
– 並列分散データ構造の重要性• 並列アプリケ
並列アプリケーションを使用したシミュレーションの手順
ションを使用したシミュレ ションの手順を一通り体験する
– データ生成（初期データ，領域分割）
データ生成（初期データ 領域分割）– 計算– ポスト処理• 対象– 有限体積法（または直接差分法）による熱伝導解析コード– 共役勾配法（CG法）による連立一次方程式求解
共役勾 法（法）
る連立 次方程式求解

T2K-FVM-01

2

中島研吾の紹介• 所属– 東京大学情報基盤センター– スーパーコンピューティング研究部門
ス パンピ
テ ング研究部門• 大学院情報理工学系研究科数理情報学専攻（兼担）• 専門分野– 計算力学，（並列）数値線形代数– 並列計算プログラミングモデル，（並列）適応格子
並列計算プ グ
グ デ（並列）適応格• これまで関わったプロジェクト– GeoFEM： 地球シミュレータ関連，並列有限要素法– HPC-MW： 戦略的基盤ソフトウェア– ppOpen-HPC：自動チューニング機構を有するアプリケーショ
ン開発・実行環境（FY.2011～）• お気軽に質問ください– nakajima(at)cc.u-tokyo.ac.jp

T2K-FVM-01

3• 5月19日（木）：1000~1200

スケジュ ル – T2Kオープンスパコン（東大）へのログイン
スケジュール• 5月19日（木）：1300~1730（予定）– 並列アプリケーション開発入門（I）
並列アプリケ ション開発入門（I）• 並列プログラミングの学び方• 有限体積法• 1CPU用計算プログラム
計算プ グ ム eps_fvm– 並列アプリケーション開発入門（II）• MPI超入門• 並列分散メッシュデータ• 5月14日（金）：0930~1200– 並列アプリケーション開発入門（II）（続き）• 並列分散メッシュデータ（続き）• 領域分割• 5月14日（金）：1300～1700– 並列アプリケーション開発入門（III）
並列アプリケ ション開発入門（III）• 並列FVMコードの開発– 実習

T2K-FVM-01

4

謝辞• サイバネットシステム株式会社
サイ ネッ シ テ 株式会社

T2K-FVM-01

5

準備（Windows）• T2Kへのログイン• WinSCP等インストール
MicroAVSインスト ル• MicroAVSインストール

T2K-FVM-01

6

注意等• 講義室内では飲食厳禁です。今後のこともありますので
ご理解，ご協力のほどをお願いいたします。
– 自動販売機は工学部12号館前，武田先端知ビル内。

• アンケ
アンケートにご協力ください
トにご協力ください。

T2K-FVM-01

7• 背景：並列プログラミングの学び方• 有限体積法の概要• １CPU用プログラム eps_fvm– メッシュ生成，形状データ
生成
状デ タ– 計算実行例– プログラムの内容

8

T2K-FVM-01

背 景• 「T2Kオープンスパコン」
T2Kオ プンスパコン」 ，，「次世代スーパーコンピュー
次世代ス パ コンピュタ」 等の開発を背景に，大規模並列シミュレーションへの期待は，産学において 層高まっている。
の期待は，産学において一層高まっている。
• 並列計算機を使いこなすためには，「並列プログラミン
並列計算機を使いこなすためには 「並列プログラミン
グ」の習得が必須である。
– 並列計算機を使いこなす，アプリーション分野の研究者，技
並列計算機を使なす
プリ シ 分野 研究者 技
術者の育成

9

T2K-FVM-01

4S型人材育成戦略「4つのS S t
「4つのS」：System，Stage,
St
Status,St t
Style（1/2）
St l （1/2）• System– SMASH– 科学技術計算の真髄• Stage：4つの段階– 並列プログラミングへの道– ③が最も重要，かつ教育困難③が最も重要 かつ教育困難• 現状はアルゴリズム中心（連続体
力学）

Science

A S① 計算機リテラシー
プログラミング言語

Modeling
d li

A S H

② 科学技術計算基礎論・演習
科学技術計算基礎論 演習（数値解析・プログラミングの基礎）

Algorithm

ここに厚い壁

Software

S M A S H

③ アプリケーションの実用的プログラ
ミング（有限要素法等）

Hardware

S M A S H

④ 並列プログラミング

10

T2K-FVM-01

並列プログラミングへの道
A S① 計算機リテラシー
グラ ング言語
プログラミング言語

A S H② 科学技術計算基礎論・演習（数値解析 プログラミングの基礎）（数値解析・プログラミングの基礎）
ここに厚い壁S M A S H

③ アプリケーションの実用的プログラ
ミング（有限要素法等）

S M A S H

④ 並列プログラミング• ③が重要であるが，プログラミング教育はほとんど実施
されされていないなT2K-FVM-01③③（実用プログラミング）を学ぶには?
• プログラミング能力をつけるためには，徹底して実アプリ
ケーションコードのソースを「読む」能力をつけることが重要（特に学部4年～大学院初級）– 英語，漢文の音読のごとく• 並列計算で大事なのはMPI等の文法ではなく，並列デー
タ構造等の設計。対象アプリケ ション，アルゴリズムにタ構造等の設計。対象アプリケーション，アルゴリズムに対する深い理解が必要⇒③が重要• これまでの経験で効果は確認済

11

12

T2K-FVM-01

科学技術アプリケーション• 局所的手法– 差分法，有限要素法，有限体積法– 疎行列• 大域的手法– 境界要素法，スペクトル法等
– 密行列

13

T2K-FVM-01

本講習の方針• まずは1CPU用のシリアルアプリケーションコードを徹底理解– アルゴリズムのやさしい有限体積法を選択– 熱伝導方程式• アプリケーションに適した並列計算のためのデータ構造
アプリケ ションに適した並列計算のためのデ タ構造（＝局所分散データ構造）を考える• 慣れているFORTRANで解説していますが，C言語版の慣れているFORTRANで解説していますが C言語版のアプリケーションも準備，ほとんど各行解説するので，余り支障は無いでしょう• MPIそのものについてはあまり詳しく説明しません。MPI
入門向け講習会は年2 3回程度実施されています
入門向け講習会は年2～3回程度実施されています。

T2K-FVM-01

14• 背景：並列プログラミングの学び方• 有限体積法の概要• １CPU用プログラム eps_fvm– メッシュ生成，形状データ
生成
状デ タ– 計算実行例– プログラムの内容

T2K-FVM-01

15

対象とするアプリケーションの概要• 支配方程式：三次元定常熱伝導方程式– 物性の温度依存性無し，等方性∂ ⎛ ∂T ⎞ ∂ ⎛ ∂T ⎞ ∂ ⎛ ∂T ⎞⎟⎟ + ⎜ λ⎜λ⎟ + ⎜⎜ λ⎟+Q = 0∂x ⎝ ∂x ⎠ ∂y ⎝ ∂y ⎠ ∂z ⎝ ∂z ⎠• 有限体積法（Finite Volume Method
Method，FVM）による空間
FVM）による空間
離散化– 任意形状の要素，要素中心で変数を定義。
任意形状の要素 要素中心で変数を定義– 直接差分法（Direct Finite Difference Method）とも呼ばれる。

• 境界条件– ディリクレ（温度固定），ノイマン（境界熱流束），体積発熱• 反復法による連立一次方程式解法
復法 よる連立 次方程式解法– 共役勾配法（CG）＋対角スケーリング（点ヤコビ）前処理

T2K-FVM-01

16

一次元熱伝導方程式：中央差分
次元熱伝導方程式 中央差分
熱伝導率λ一定d ⎛ dT ⎞⎜λ⎟ +Q = 0dx ⎝ dx ⎠i

Τi-1

Τi
Δx

Τi+1
Δx

⎛ dT ⎞⎛ dT ⎞
Ti +1 − Ti
T −T
− ⎜λ
⎟
⎜λ
⎟
λ
− λ i i −1d ⎛ dT ⎞ ⎝ dx ⎠i +1/ 2 ⎝ dx ⎠i −1/ 2
Δx
Δx
=
⎜λ
⎟ =
dx ⎝ dx ⎠i
Δx
Δx⎛ Ti +1 − 2Ti + Ti −1 ⎞
= λ⎜
⎟
2
Δx
⎝
⎠

T2K-FVM-0117

有限体積法による空間離散化
熱流束に関するつりあい式
隣接要素との熱伝導

∑d
k

a

Sik
ik

λi

+

(
Tk − Ti ) + ∑ Sid q&id + Vi Q& i = 0
d
ki

λk

d

要素境界面 体積発熱
通過熱流束

b
dai
dbi

Sib
dia

Sia

dib

i

dic

Sic

Sid
qid

dci

c

λ ：熱伝導率
Vi ：要素体積
S ：表面面積
表面面積
dij ：要素中心から表面までの距離
q ：表面フラックス
Q ：体積発熱
体積発熱

T2K-FVM-01

18

一次元差分法との比較（1/3）
次元差分法との比較（1/3）

Δh

λi

λa
dai
a

dia

λb
dib

i

Sia

一辺の長さΔhの正方形メッシュ
接触面積
接触面積：
Sik= Δh
要素体積：
Vi = Δh2
接触面までの距離： dij=Δh/2
各メッシュの熱伝導率： λi

dbi
Sib

b
厚さ：1

この面を通過する熱量：Qsib
Tb − Ti
Tb − Ti
Qsib = −λ⋅ Sib = −⋅ Sib
d ib d bi
d ib + d bi+

Tb − Ti
Qsib = −⋅S
d ib d bi ib
+

λi

λb

λ

λ

熱伝導率が一様の場合
フーリエの法則
面を通過する熱量（熱伝導率）×（温度勾配）＝ー（熱伝導率）×（温度勾配）

熱伝導率が要素ごとに異なる場合 調和平均
異なる場合：調和平均

T2K-FVM-01

19

一次元差分法との比較（2/3）
次元差分法との比較（2/3）

Δh

λi

λa
dai
a

Sia

dia

λb
dib

dbi

i

b

Sib

厚さ：1

∑d
k

両辺をViで割る：
この部分に注目すると

1
Vi

Sik
ik

λi

∑
k

+

一辺の長さΔhの正方形メッシュ
接触面積
接触面積：
Sik= Δh
要素体積：
Vi = Δh2
接触面までの距離： dij=Δh/2
各メッシュの熱伝導率： λi(
Tk − Ti ) + ∑ Sid q&id + Vi Q& i = 0
d
ki

d

λk

1
(
Tk − Ti ) +d ik d ki
Vi
+
Sik

λi

λk

∑
d

Sid q&id + Q& i = 0

T2K-FVM-01

20

一次元差分法との比較（3/3）
次元差分法との比較（3/3）

Δh

λi

λa
dai
a

dia

Sia

熱伝導率一様として

1
Vi

∑
k

λb
dib

i

dbi
Sib

b
厚さ：1

一辺の長さΔhの正方形メッシュ
接触面積
接触面積：
Sik= Δh
要素体積：
Vi = Δh2
接触面までの距離： dij=Δh/2
各メッシュの熱伝導率： λi1Δh(T − T ) =(T − T )d ik d ki k i (Δh )2 k = a ,b Δh / 2 Δh / 2 k i
+
+

∑

Sik

λi

λk

λ

λ

1
Δh
1
Δh
1
(
(
=
Tk − Ti ) =Tk − Ti ) =λ (Tk − Ti )222(Δh ) k =a ,b Δh + Δh(Δh ) k =a ,b Δh(Δh ) k =a ,b
2λ 2λ
λ
⎡ 1
⎤
⎡ Ta − 2Ti + Tb ⎤1
(Ta − Ti ) +(Tb − Ti )⎥ = λ ⋅ ⎢= λ ⋅⎢
⎥
2
2
2
(
)
(
)
(
)
Δ
Δ
h
Δ
h
h
⎦
⎣
⎣
⎦

∑

∑

∑

T2K-FVM-01

21

有限体積法
• セル中心型有限体積法– 「直接差分法」– ガウスグリーン型有限体積法• 任意形状の要素を扱うことが可能– 均等メッシュでない場合は空間について一次精度• 各要素の体積，要素間の接続関係（接触面積，重心（本当は 外心」）から接触面 の距離）などを予め指定する
当は「外心」）から接触面への距離）などを予め指定する必要があるため，メッシュ作成は多少面倒• 境界面，ディリクレ境界条件（温度固定）の扱いが問題
境界面 ディリクレ境界条件（温度固定）の扱いが問題

T2K-FVM-01

22

ガウスの定理：Gauss’s Theorem⎛ ∂U ∂V ∂W ⎞⎟⎟ dV = (Unx + Vn y + Wnz ) dS⎜⎜++∂x ∂y ∂z ⎠
V ⎝
S

∫

∫

n

• 三次元デカルト座標（x,y,z）
三次元デカルト座標（）• 滑らかな閉曲面Sによって囲まれたV• V内で定義される，3つの連続関数– U(x,y,z)( y )，V(x,y,z)( y )，W(x,y,z)( y )• 曲面S上で外向きに引いた法線ベクトルn– nx，ny，nz：方向余弦

SV

T2K-FVM-0123ベクトル表記すると• ガウスの定理
ガ
定
T
∇
⋅
w
dV
=
w
∫
∫ n dS

VS• グリーンの定理∫

VvΔu dV =

T
(
)
(
v
∇
u
n
dS
−
∇
∫
∫ v )(∇u ) dVTSVT2K-FVM-0124

有限体積法による空間離散化
各要素の境界面を通過する熱量に着目
隣接要素との熱伝導

∑d
k

a

Sik
ik

λi

+

(
Tk − Ti ) + ∑ Sid q&id + Vi Q& i = 0
d
ki

λk

d

要素境界面 体積発熱
通過熱流束

b
dai
dbi

Sib
dia

Sia

dib

i

dic

Sic

Sid
qid

dci

c

λ ：熱伝導率
Vi ：要素体積
S ：表面面積
表面面積
dij ：要素中心から表面までの距離
q ：表面フラックス
Q ：体積発熱
体積発熱

T2K-FVM-01

25

ディリクレ境界条件（境界温度固定）

単純に要素中心で指定
単純に要素中心指定（精度悪化）

境界面からの等価なフラックスを加算（面積，重心との距離必要）「薄い」要素を設定「薄要素を設定（～10-20）

T2K-FVM-01

26

ディリクレ境界条件（境界温度固定）

単純に要素中心で指定
単純に要素中心指定（精度悪化）

SB

λi

TB

Ti

dB

境界面からの等価なフラックスを加算（面積，重心との距離必要）
現在はこの方式を使用

SB(TB − Ti )q=d B / λi

「薄い」要素を設定「薄要素を設定（～10-20）

T2K-FVM-0127

有限体積法による空間離散化
熱流束に関するつりあい式：ディリクレ境界条件考慮
隣接要素との熱伝導

∑d
k

Sik
ik

λi

+

d ki

温度固定境界(Tk − Ti ) + ∑
e

λk

∑

Sie
(
TiBe − Ti ) +Sid q&id + Vi Q& i = 0
d ie
d

λi

要素境界面 体積発熱
通過熱流束

a
b
dai
dbi

Sib
Sia

dia

die

TiBe

Sie
Sid

dib

i

dic
dci

Sic
qid

c

λ ：熱伝導率
Vi ：要素体積
S ：表面面積
表面面積
dij ：要素中心から表面までの距離
q ：表面フラックス
Q ：体積発熱
体積発熱
TiB ：境界温度

T2K-FVM-0128

有限体積法による空間離散化
隣接している要素の影響のみ受ける
隣接要素との熱伝導

∑d
k

Sik
ik

λi

+

d ki

温度固定境界(Tk − Ti ) + ∑
e

λk

∑

Sie
(
TiBe − Ti ) +Sid q&id + Vi Q& i = 0
d ie
d

λi

要素境界面 体積発熱
通過熱流束

a
b
dai
dbi

Sib
Sia

dia

die

TiBe

Sie
Sid

dib

i

dic
dci

Sic
qid

c

λ ：熱伝導率
Vi ：要素体積
S ：表面面積
表面面積
dij ：要素中心から表面までの距離
q ：表面フラックス
Q ：体積発熱
体積発熱
TiB ：境界温度

T2K-FVM-0129係数行列は疎行列（SparsepMatrix）• 0が多い多• 通常「数値解析」などで扱うのは「密行列」• 非ゼロ成分のみを記憶する方法– Compressed Row Storage, CRS（後述）• 密行列と比較して取り扱いが困難

T2K-FVM-01

30

現行コ ドにおける制限
現行コードにおける制限• 要素形状は六面体のみ。
– 各要素は各辺長さ1の立方体– 全体形状，境界条件等はメッシュジェネレータで規定• 隣接要素数は要素各面に接する6面までとする。
– プログラムを変更すれば増やすことは可能
プ グラムを変更すれば増やす とは可能

• 境界条件– ディリクレ境界条件– ノイマン境界条件（面フラックス）– 体積発熱• 要素番号は1から始まって連続していなければならない。
要素番号は1から始ま て連続していなければならない

T2K-FVM-01

31

現在対象としている形状，境界条件
プログラムでは，後掲のデータ形式に従えば，任意の形状を扱うことができる
X=Xmax
ディリクレ境界条件
T=0
NZ

一様発熱
（体積フラ クス 1 0）
（体積フラックス=1.0）

X=Xminノイマン境界条件-λ(dT/dX)=1.

NY

z
y

NX
ｘ

T2K-FVM-01

32• 背景：並列プログラミングの学び方• 有限体積法の概要• １CPU用プログラム eps_fvm– メッシュ生成，形状データ
生成
状デ タ– 計算実行例– プログラムの内容

T2K-FVM-01

33

T2Kで利用できるコンパイラ• 日立コンパイラ（デフォルト）
日立 ンパイラ（デフ ルト）• Intel• PGI• GNU• 本講義では「日立コンパイラ」を使用する– mpif90 -Oss –noparallel s1-3.f– mpicc -Os –noparallel s1-3
s1-3.cc

T2K-FVM-01

34

ファイルシステム（1/2）• /home– 各コースの制限値まで利用可能（クラス全体で4TB）– バックアップなし– 全ノードで共有• /short– 容量無制限– 最終アクセスから5日で消える– 全ノードで共有
全
ド 共有• /tmp– ノードあたり
ドあ
140GB– ノードローカル– バッチジョブ終了時に消えるバ チジ ブ終了時に消える

T2K-FVM-01

35

ファイルシステム（2/2）• 分散ファイルシステム
分散フ イルシステム HSFS
– 全ノードで共有するファイルシステムにはHSFSが使われている– 16台のファイルサーバに負荷分散– ユーザーには単一のディレクトリツリーが見える
ユ ザ には単 のディレクトリツリ が見える
– ファイルストライプモード
ファイルストライプモ ド• ファイル単位で処理するサーバを割り当てる
同じファイルに多数のノ ドからアクセスすると破綻する• 同じファイルに多数のノードからアクセスすると破綻する• 現在はLustreも使えます：速いです– /lustre/username

T2K-FVM-01

36

必要ファイルインストール：T2K
Lustreを使いたければ下記へ移動$> cd /$
/lustre/t003xx
ust e/t003（自分のUID））（自分のU
FORTRAN$> cp
p /home/t00000/fvm-f.tar .
$> tar xvf fvm-f.tar
C$> cp /home/t00000/fvm-c.tar .
$> tar xvf fvm-c.tar
それぞれ直下に fvm-f，fvm-c というディレクトリができます。
それぞれ直
と うデ
ク が きますこれらを今後<$FVM>と呼びます。

T2K-FVM-01

37.bash_profile_p$> cd*****************************************$> cp /home/t00000/.bash_profile .
または ~/.bash_profileに以下の2行を挿入~/ bash profileに以下の2行を挿入
PATH="$PATH":/home/t00000/bin
export PATH******************************************$> source .bash_p$
profile$> printenv
PATHに以下が入っていることを確認
home/t00000/bin
次回からはログイン時に自動的に起動

T2K-FVM-01

38

メッシュ生成からやってみよう$> cd <$FVM>/run$> cat fvmmg.ctrl
g
2 2 3$> eps_fvm_mg$> ls fvm_entire*fvm entire*
fvm_entire_mesh.dat
fvm_entire_mesh.inp
fvm_entire_mesh.inp
p_g
geo

• 入力
– fvmmg.ctrl：各座標軸方向の要素数
fvmmg ctrl：各座標軸方向の要素数 （NX,NY,NZ）（NX NY NZ）#MGCTRL• この場合2×2×3=12個の要素が生成される• 出力– fvm_entire_mesh.dat: 初期全体メッシュデータ– fvm_entire_mesh.inp: AVS表示用ファイル#S-GRID#S-GRID-UCD• メッシュ形状をmicroAVSで表示するためのファイル– fvm_entire_mesh.inp-geo:AVS表示用形状ファイル #S-GRID-GEO

T2K-FVM-01

39

メッシュ生成：専用ツ ル使用
メッシュ生成：専用ツール使用#MGCTRL<$FVM>/run/fvmmg.ctrl

各座標軸方向の要素数（名称固定）

メッシュジェネレータ
eps fvm mg
eps_fvm_mg#S-GRID<$FVM>/<$FVM>/run//fvm_entire_mesh.dat#S-GRID-UCD<$FVM>/<$FVM>/run//fvm_entire_mesh.inp

fvm_entire_mesh.inp-geo

初期全体メッシュデータ（名称固定）

初期全体メッシュデータ
UCDファイル（AVS用）（名称固定）

初期全体メッシュデータ
UCDファイル形状部（AVS用）（名称固定）#S-GRID-GEO<$FVM>/run/T2K-FVM-01

初期全体メッシュデータ（例）（1/6）
http://nkl cc u tokyo ac jp/tutorial/part tutorial
http://nkl.cc.u-tokyo.ac.jp/tutorial/part_tutorial

12
1
2
3
4
5
6
7
8
9
10
11
12
20
1
1
1
2
2
3
3
4
5
5
5
6
6
7
7
8
9
9
10
11
6
2
4
6
8
10
12
6
1
3
5
7
9
11
4
5
6
7
8

1.000000E+00
1.000000E+00
1.000000E+00
1.000000E+00
1.000000E+00
1.000000E+00
1.000000E+00
1.000000E+00
1.000000E+00
1.000000E+00
1.000000E+00
1.000000E+00
2
3
5
4
6
4
7
8
6
7
9
8
10
8
11
12
10
11
12
12

1.000000E+00
1.000000E+00
1.000000E+00
1.000000E+00
1.000000E+00
1.000000E+00
1.000000E+00
1.000000E+00
1.000000E+00
1.000000E+00
1.000000E+00
1.000000E+00

11.000000E+00
000000E+00
1.000000E+00
1.000000E+00
1.000000E+00
1.000000E+00
1.000000E+00
1.000000E+00
1.000000E+00
1.000000E+00
1.000000E+00
1.000000E+00
1.000000E+00
1.000000E+00
1.000000E+00
1.000000E+00
1.000000E+00
1.000000E+00
1.000000E+00
1.000000E+00
1.000000E+00

55.000000E-01
000000E 01
5.000000E-01
5.000000E-01
5.000000E-01
5.000000E-01
5.000000E-01
5.000000E-01
5.000000E-01
5.000000E-01
5.000000E 01
5.000000E-01
5.000000E-01
5.000000E-01
5.000000E-01
5.000000E-01
5.000000E-01
5.000000E-01
5.000000E-01
5.000000E-01
5.000000E-01
5.000000E-01

1.000000E+00
1.000000E+00
1.000000E+00
1.000000E+00
11.000000E+00
000000E 00
1.000000E+00

5.000000E-01
5.000000E-01
5.000000E-01
5.000000E-01
55.000000E-01
000000E 01
5.000000E-01

1.000000E+00
1.000000E+00
1.000000E+00
1.000000E+00
1.000000E+00
11.000000E+00
000000E+00

1.000000E+00
1.000000E+00
1.000000E+00
1.000000E+00
1.000000E+00
11.000000E+00
000000E+00

1.000000E+00
1.000000E+00
1.000000E+00
1.000000E+00

5.000000E-01
1.500000E+00
5.000000E-01
1.500000E+00
5.000000E-01
1.500000E+00
5.000000E-01
1.500000E+00
5.000000E-01
1.500000E+00
5.000000E-01
1.500000E+00

5.000000E-01
5.000000E-01
1.500000E+00
1.500000E+00
5.000000E-01
5.000000E-01
1.500000E+00
1.500000E+00
5.000000E-01
5.000000E-01
1.500000E+00
1.500000E+00

55.000000E-01
000000E 01
5.000000E-01
5.000000E-01
5.000000E-01
5.000000E-01
5.000000E-01
5.000000E-01
5.000000E-01
5.000000E-01
5.000000E 01
5.000000E-01
5.000000E-01
5.000000E-01
5.000000E-01
5.000000E-01
5.000000E-01
5.000000E-01
5.000000E-01
5.000000E-01
5.000000E-01
5.000000E-01

0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
00.000000E+00
000000E 00
0.000000E+00

5.000000E-01
5.000000E-01
5.000000E-01
5.000000E-01
1.500000E+00
1.500000E+00
1.500000E+00
1.500000E+00
2.500000E+00
2.500000E+00
2.500000E+00
2.500000E+00

各要素情報

要素間コネクティビティ面に隣接する要素

ディリクレ境界条件
温度固定
ノイマン境界条件
境界熱流束
体積発熱

40

T2K-FVM-01

41

変数表
変数名

型

配列サイズ

NODE tot
NODE_tot

I

-

総要素数

NODE_VOL(:)

R

NODE_tot

要素体積

NODE_COND(:)

R

NODE_tot

要素熱伝導率NODE_XYZ(:)

R

3*NODE_tot

要素重心座標（3次元）

CONN_totI-

コネクティビティ総数

CONN_node(:)

I

2*CONN_tot

コネクティビティ構成要素

CONN_COEF(:)

R

CONN_tot

コネクティビティ係数

FIX_NODE_tot

I-

ディリクレ境界条件適用要素数

FIX_NODE_ID(:)

I

FIX_NODE_tot

ディリクレ境界条件適用要素番号

FIX_NODE_COEF(:)

R

FIX_NODE_tot

ディリクレ境界条件係数

FIX_NODE_VAL(:)

R

FIX_NODE_tot

ディリクレ境界条件値

SURF_NODE_tot

I-

ノイマン境界条件適用要素数

SURF_NODE_ID(:)

I

SURF_NODE_tot

ノイマン境界条件適用要素番号

SURF_NODE_FLUX(:)

R

SURF_NODE_tot

ノイマン境界条件フラックス

BODY_NODE_tot

I-

体積発熱境界条件適用要素数

BODY_NODE_ID(:)

I

BODY_NODE_tot

体積発熱境界条件適用要素番号

BODY NODE FLUX( )BODY_NODE_FLUX(:)

R

BODY NODEBODY_NODE_tot

体積発熱境界条件
体積発熱境界条件フラックスク

内容

T2K-FVM-01

42

初期全体メッシュデータ（例）（2/6）
各要素
12 要素数:NODE_tot
11.000000E+00
21.000000E+00
31.000000E+00
41.000000E+00
51.000000E+00
61.000000E+00
71.000000E+00
81.000000E+00
91.000000E+00
10
11.000000E+00
000000E+00
111.000000E+00
121.000000E+00
要素番号 要素体積
NODE_VOL(i)

9
5

z
y

1
ｘ

11
7
3

1.000000E+00
1.000000E+00
1.000000E+00
1.000000E+00
1.000000E+00
1.000000E+00
1.000000E+00
1.000000E+00
1.000000E+00
11.000000E+00
000000E+00
1.000000E+00
1.000000E+00
要素熱伝導率NODE_COND(i)5.000000E-01
5.000000E-01
5.000000E-01
1.500000E+00
5.000000E-01
5.000000E-01
5.000000E-01
1.500000E+00
5.000000E-01
1.500000E+00
1.500000E+00
5.000000E-01
5.000000E-01
5.000000E-01
1.500000E+00
1.500000E+00
5.000000E-01
1.500000E+00
5.000000E-01
1.500000E+00
1.500000E+00
1.500000E+00
1.500000E+00
1.500000E+00
5.000000E-01
5.000000E-01
2.500000E+00
11.500000E+00
500000E+00
5.000000E-01
5 000000E-01
2.500000E+00
2 500000E+00
5.000000E-01
1.500000E+00
2.500000E+00
1.500000E+00
1.500000E+00
2.500000E+00
要素中心X座標
要素中心Y座標
要素中心Z座標
NODE_XYZ(3*i-2) NODE_XYZ(3*i-1) NODE_XYZ(3*i)

1210read (IUNIT, '(10i10)') NODE_tot
do i=i 1, NODE_tot
NODE tot
read (IUNIT,'(i10,5e16.6)')ii, NODE_VOL(i), NODE_COND(i),(NODE_XYZ(3*i-3+k), k=1, 3)
enddo

6
2

8
4

&
&

T2K-FVM-01

43

初期全体メッシュデータ（例）（3/6）
S
要素間コネクティビティ：ある面に接続する2要素
20
1
1
1
2
2
3
3
4
5
5
5
6
6
7
7
8
9
9
10
11
E1

コネクティビティ総数： CONN_tot
21.000000E+00
31.000000E+00
51.000000E+00
41.000000E+00
61.000000E+00
411.000000E+00
000000E+00
71.000000E+00
81.000000E+00
61.000000E+00
71.000000E+00
91.000000E+00
81.000000E+00
101.000000E+00
81.000000E+00
111.000000E+00
121.000000E+00
101.000000E+00
111.000000E+00
121.000000E+00
121.000000E+00
1 000000E+00
E2
S:要素境界面積

E1= CONN_NODE(2*ic-1)
E2= CONN_NODE(2*ic)

5.000000E
5.000000E-01
01
5.000000E-01
5.000000E-01
5.000000E-01
5.000000E-01
55.000000E-01
000000E 01
5.000000E-01
5.000000E-01
5.000000E-01
5.000000E-01
5.000000E-01
5.000000E-01
5.000000E-01
5.000000E-01
5.000000E
5.000000E-01
01
5.000000E-01
5.000000E-01
5.000000E-01
5.000000E-01
55.000000E-01
000000E 01
d1：E1重心～
境界面距離5.000000E
5.000000E-01
01
5.000000E-01
5.000000E-01
5.000000E-01
5.000000E-01
5.000000E-01
5 000000E 01
5.000000E-01
5.000000E-01
5.000000E-01
5.000000E-01
5.000000E-01
5.000000E-01
5.000000E-01
5.000000E-01
5.000000E
5.000000E-01
01
5.000000E-01
5.000000E-01
5.000000E-01
5.000000E-01
5.000000E-01
5 000000E 01
d2：E2重心～
境界面距離

d1

d2

E1

E2

i

ik

ki

k

係数行列の非対角成分

read (IUNIT,'(10i10)') CONN_tot
do i=i 1, CONN_tot
CONN tot
read (IUNIT,'( 2i10, 3e16.6)')(CONN_NODE(2*i-2+k), k= 1, 2),
AREA, d1, d2
enddo

&
&

T2K-FVM-01

44

初期全体メッシュデータ（例）（3/6）
要素間コネクティビティ：5番要素に着目
1
5
5
5
E1

5
6
7
9
E2

1.0E+00
5.0E-01
1.0E+00
5.0E-01
1.0E+00
5.0E-011
1.0E+000E+00
55.0E-01
0E-01
S:要素境界面積 d1：E1重心～
境界面距離

5.0E-01
5.0E-01
5.0E-0155.0E-01
0E-01
d2：E2重心～
境界面距離

各要素は一辺の長さ1の立方体：• 要素境界面積は1×1＝1• 各要素重心から境界面までの距離は0.50

9
5
1

11

12
10

9

7

8

5

4

1

3

6
2

7

6

d1

d2

E1

E2

i

ik

ki

k

係数行列の非対角成分

T2K-FVM-01

45

初期全体メッシュデータ（例）（4/6）
ディリクレ境界条件（表面温度固定）
6 ディリクレ境界条件を与える要素数： FIX_NODE_tot
21.000000E+00
5.000000E-01
0.000000E+00
41.000000E+00
5.000000E-01
0.000000E+00
61.000000E+00
5.000000E-01
0.000000E+00
81.000000E+00
5.000000E-01
0.000000E+00
101.000000E+00
5.000000E-01
0.000000E+00
121.000000E+00
5.000000E-01
0.000000E+00
要素番号 境界面面積S
境界面と要素重心距離d 境界値TB
FIX NODE ID(ib)FIX_NODE_ID(ib)FIX
FIX_NODE_VAL(ib)
NODE VAL(ib)951read (IUNIT,'(10i10)') FIX_NODE_tot
do i= 1, FIX_NODE_tot
read (IUNIT, '(i10,(i10, 3e16.6)')
3e16.6) )FIX_NODE_ID(i), S, d,
FIX_NODE_VAL(i)
enddo

11

12
10

7

8

3

6
2

4

S

&
&

d
T=TB

T2K-FVM-01

46

1つの「要素」に複数のディリクレ
境界条件を与えることも実は可能
7 ディリクレ境界条件を与える要素数： FIX_NODE_tot
21.000000E+00
5.000000E-01
0.000000E+00
41.000000E+00
5.000000E-01
0.000000E+00
611.000000E+00
000000E+00
55.000000E
000000E-0101
00.000000E+00
000000E+00
81.000000E+00
5.000000E-01
0.000000E+00
101.000000E+00
5.000000E-01
0.000000E+00
121.000000E+00
5.000000E-01
0.000000E+00
101.000000E+00
5.000000E-01
1.000000E+01
要素番号 境界面面積S
境界面と要素重心距離d 境界値TB
FIX_NODE_ID(ib)FIX NODE VAL(ib)FIX_NODE_VAL(ib)

9
5
1

11

12
10

7

8

3

6
2

4

TB=10

ある表面の温度を規定するだけであって，その要素の温度自体（要素中心）を規定するわけではない。「要素数」というよりは「面数」という方がより正確がより正確。

10

TB=0

T2K-FVM-0147

有限体積法による空間離散化
熱流束に関するつりあい式
隣接要素との熱伝導

∑d
k

Sik
ik

λi

+

d ki

温度固定境界(Tk − Ti ) + ∑
e

λk

∑

Sie
(
TiBe − Ti ) +Sid q&id + Vi Q& i = 0
d ie
d

λi

要素境界面 体積発熱
通過熱流束

a
b
dai
dbi

Sib
Sia

dia

die

TiBe

Sie
Sid

dib

i

dic
dci

Sic
qid

c

λ ：熱伝導率
Vi ：要素体積
S ：表面面積
表面面積
dij ：要素中心から表面までの距離
q ：表面フラックス
Q ：体積発熱
体積発熱
TiB ：境界温度

T2K-FVM-01

48

各要素に適用できる境界条件の数に
制限は無い
q2

TBｂ

i
TBa

q1

T2K-FVM-01

49

初期全体メッシュデータ（例）（5/6）
ノイマン境界条件（表面熱フラックス：W/m2）
6 ノイマン境界条件を与える要素数： SURF_NODE_tot
11.000000E+00
1.000000E+00
31.000000E+00
1.000000E+00
51.000000E+00
1.000000E+00
71.000000E+00
1.000000E+00
91.000000E+00
1.000000E+00
111.000000E+00
1.000000E+00
要素番号 境界表面積S
表面フラックスQSSURF NODE ID(ib)
SURF_NODE_ID(ib)read (IUNIT,'(10i10)') SURF_NODE_tot
do i= 1, SURF_NODE_tot(IUNIT,, '(i10,(, 3e16.6)')) )read (
SURF_NODE_ID(i), S, QS
enddo

9
5
1
S
&

QS

1つの要素に複数のノイマン境界条件を適用することは可能（複数の面に異なった熱流束を
与える）。
与える）

11
7
3

12
10
6
2

8
4

T2K-FVM-01

50

初期全体メッシュデータ（例）（6/6）
体積発熱（W/m3）
4 体積発熱を与える要素数
51.000000E+00
61.000000E+00
71.000000E+00
81.000000E+00
要素番号 体積フラックスQV

read (IUNIT,'(10i10)') BODY_NODE_tot__
do i= 1, BODY_NODE_tot
read (IUNIT, '(i10, 3e16.6)') icel, QV
enddo

9
5
1

11

12
10

7

8

3

6
2

4

QV

同じ要素に対して，複数回体積発熱を指定した場合には， 番最後に指定した値が適用さた場合には，一番最後に指定した値が適用される。

T2K-FVM-01

51

MicroAVSによる表示#S-GRID-UCD:<$FVM>/run/fvm_entire_mesh.inp__• AVS，microAVSのUCDファイルフォーマット
（Unstructured Cell Data）を使用
Data）を使用。
– http://kgt.cybernet.co.jp/

T2K-FVM-01

52

MicroAVSの使用• 手順は以下– http://nkl.cc.u-tokyo.ac.jp/10e/mavs/

T2K-FVM-01

53

UCDフォ マットについて（1/4）
UCDフォーマットについて（1/4）
MicroAVSで扱うことのできる要素形状
要素の種類

キ ワ ド
キーワード

点

pt

線

line

三角形

tri

四角形

quad

四面体

tet

角錐

pyr

三角柱

prism

六面体

hex

二次要素
線2

line2

三角形2

tri2

四角形2

quad2

四面体2

tet2

角錐2

pyr2
2

三角柱2

prism2六面体2

hex2

T2K-FVM-01

54

UCDフォ マットについて（2/4）
UCDフォーマットについて（2/4）
書式の概要• ファイルの拡張子–デ
データファイルの拡張子はタファイルの拡張子は “.inp
inp”

• 書式
– アスキーファイル
アスキ ファイル
– 時系列（複数ステップ）に対応したものが標準フォーマット
– MicroAVS Ver
Ver.6.0（現在はVer.13.0）まで使用していた旧
6 0（現在はVer 13 0）まで使用していた旧フォーマット（単一ステップデータのための書式）も読み込むことは可能• 実は一部これを使用している

T2K-FVM-01

55

UCDフォ マットについて（4/3）
UCDフォーマットについて（4/3）
書式の概要(要素のデ タ成分数) (成分1の構成数) (成分2の構成数) ・・・(各成分の構成数)(要素のデータ成分数)(各成分の構成数)(要素データ成分1のラベル)，(単位)(要素データ成分2のラベル)，(単位)(コメント行)(ステップ数)(データの繰り返しタイプ)(ステップ番号1) (コメント)(全節点数) (全要素数)(節点番号1) (X座標) (Y座標) (Z座標)(節点番号2) (X座標) (Y座標) (Z座標)

・
・
・

・
・
・

(各要素データ成分のラベル)，(単位)(要素番号1) (要素データ1) (要素データ2) ・・・・・(要素番号2) (要素データ1) (要素データ2) ・・・・・(要素番号1) (材料番号) (要素の種類) (要素を構成する節点のつながり)(要素番号2)) (材料番号) (要素の種類)(要素番号(要素 種類) (要素を構成す(要素を構成する節点のつながり)
節点
り)

・
・
・

・
・
・
(各節点のデータ数) (各要素のデータ数)(節点のデータ成分数)(節点のデ
タ成分数) (成分1の構成数) (成分2の構成数) ・・・(各成分の構成数)(節点データ成分1のラベル)，(単位)(節点データ成分2のラベル)，(単位)・・・(各節点データ成分のラベル)，(単位)(節点番号1) (節点データ1) (節点データ2) ・・・・・(節点番号2) (節点データ1) (節点データ2) ・・・・・・・・(ステップ番号2) (コメント) (全節点数) (全要素数)
・・・

T2K-FVM-01

56

UCDフォ マットについて（4/4）
UCDフォーマットについて（4/4）
旧フォーマット(全節点数) (全要素数) (各節点のデータ数) (各要素のデータ数) (モデルのデータ
数)(節点番号1) (X座標) (Y座標) (Z座標)(節点番号2) (X座標) (Y座標) (Z座標)・・・(要素のデータ成分数)
デ(成分1の構成数) (成分2の構成数) ・・・(各成分の構成数)
各(要素データ成分1のラベル)，(単位)(要素データ成分2のラベル)，(単位)・・・(要素番号1) (材料番号) (要素の種類) (要素を構成する節点のつながり)(要素番号2) (材料番号) (要素の種類) (要素を構成する節点のつながり)(各要素データ成分のラベル)，(単位)(要素番号1) (要素データ1) (要素データ2) ・・・・・(要素番号2) (要素データ1) (要素データ2) ・・・・・

・
・
・

・
・
・

(節点のデータ成分数) (成分1の構成数) (成分2の構成数) ・・・(各成分の構成数)(節点データ成分1のラベル)，(単位)(節点データ成分2のラベル)，(単位)・・・(各節点データ成分のラベル)，(単位)(節点番号1) (節点データ1) (節点データ2) ・・・・・(節点番号2) (節点データ1) (節点データ2) ・・・・・・・・

T2K-FVM-01

57

UCDデータの例
デ タ 例
#
1
data
step1
36 12
1 0.0 0.0
2 1.0 0.0
3 2.0 0.0
4 0.0 1.0
5 1.0
1 0 1.0
1 0
6 2.0 1.0…
31 0.0 1.0
32 1.0 1.0
33 2.0 1.0
34 0.0 2.0
35 1.0 2.0
36 2.0 2.0
1 1 hex1
2 1 hex2
3 1 hex4…
10 1 hex 20
11 1 hex 22
12 1 hex 23
0 1
1 1
Mesh,
1 1.002 1.003 1.00
4 1
1.00
00
…
10 1.00
11 1.00
12 1.00

新フォーマット
0.0
0.0
0.0
0.0
0
0.0
0
0.0
3.0
3.0
3.0
3.0
3.0
3.0
2
3
5

5
6
8

4
5
7

10
11
13

11
12
14

14
15
17

13
14
16

21
23
24

24
26
27

23
25
26

29
31
32

30
32
33

33
35
36

32
34
35

36 12 0 1
1 0.0 0.0
2 1.0 0.0
3 2.0 0.0
4 0.0 1.0
5 1.0 1.0
6 2.0 1.0…
31 0.0 1.0
32 1
1.0
0 1
1.0
0
33 2.0 1.0
34 0.0 2.0
35 1.0 2.0
36 2.0 2.0
1 1 hexe1
2 1 hex2
3 1 hex4…
10 1 hex 20
11 1 hex 22
12 1 hhex 23
1 1
Mesh,
1 1.002 1.003 1.004 1.00…
10 1.00
11 1.00
12 1.00

0
0.0
0.0
0.0
0.0
0.0
0.0

旧フォーマット

3.0
3.0
3 0
3.0
3.0
3.0
3.0
2
3
5

5
6
8

4
5
7

10
11
13

11
12
14

14
15
17

13
14
16

21
23
24

24
26
27

23
25
26

29
31
32

30
32
33

33
35
36

32
34
35

T2K-FVM-01

58• 背景：並列プログラミングの学び方• 有限体積法の概要• １CPU用プログラム eps_fvm– メッシュ生成，形状データ
生成
状デ タ– 計算実行例– プログラムの内容

T2K-FVM-01

59

「eps
eps_fvm」計算の実行
fvm」計算の実行$> cd <$FVM>/serial$> make…
$> cd ../run$> ls eps_fvm
eps_fvmf

ロードモジュールの生成確認$> ./eps_fvm• 実行形式 eps_fvmは<$FVM>/run に生成される。
• 入力– <$FVM>/run/fvmmg.ctrl#MGCTRL– <$FVM>/run/fvm_entire_mesh.dat#S-GRID– ファイルを<$FVM>/runに予め生成しておく必要がある。

• 出力– <$FVM>/run/fvm_entire_mesh_results.inp（名称固定）• MicroAVS用結果出力#S-GRID-R-UCD

T2K-FVM-01

60

シミュレーションコード：eps
シミュレ
ションコ ド：eps_fvm
fvm#MGCTRL<$FVM>/run/fvmmg.ctrl#S-GRID<$FVM>/run/$FVM //
fvm_entire_mesh.dat

各座標軸方向の要素数（名称固定）

初期全体メッシュデータ（名称固定）

eps_fvm#S-GRID-R-UCD<$FVM>/run/fvm_entire_mesh_results.inp

結果ファイルAVS表示用（名称固定）

T2K-FVM-01

61

計算結果

• 「fvm_entire_mesh_results.inp」をMicroAVSによって処理することができる。

T2K-FVM-01

62• 背景：並列プログラミングの学び方• 有限体積法の概要• １CPU用プログラム eps_fvm– メッシュ生成，形状データ
生成
状デ タ– 計算実行例– プログラムの内容

T2K-FVM-01

63

プログラムの概要$> cd <$FVM>/serial$> catt eps_fvm.f
ff

メインプログラム

program eps_fvm
use eps_fvm_all
implicit REAL*8 (A-H,O-Z)
call eps
eps_fvm_input_grid
fvm input gridcall poi_gen
call eps_fvm_solve
call output_ucd
end program eps_fvm

T2K-FVM-01

64

「eps
p _fvm」処理：メイン（1/2）
」処（ ）eps_fvm_all
program eps_fvm
fuse eps_fvm_all
implicit REAL*8 (A-H,O-Z)
call eps_fvm_input_grid
call poi_gen
call eps_fvm_solver
call
ll output_ucd
t td

変数の内容について記載したモ
ジュールブロック（コモンブロックのうなも ）。
ようなもの）。
詳細はチュートリアル参照

end program eps_fvm

http://nkl.cc.u-tokyo.ac.jp/tutorial/eps_fvm_tutorial//////
http://nkl.cc.u-tokyo.ac.jp/tutorial/eps_fvm_tutorial.tar

module eps_fvm_alluse eps_fvm_util
use eps_fvm_pcg
fuse appl_cntl
end module eps_fvm_all

T2K-FVM-01

65

eps fvm の処理
eps_fvm#MGCTRL<$FVM>/run/fvmmg.ctrl#S-GRID<$FVM>/run/$FVM //
fvm_entire_mesh.dat

各座標軸方向の要素数（名称固定）

初期全体メッシュデータ（名称固定）

eps_fvm#S-GRID-R-UCD#SGRID R UCD<$FVM>/run/fvm_entire_mesh_results.inp

結果ファイルAVS表示用（名称固定）

T2K-FVM-01

66

「eps
p _fvm」処理：メイン（2/2）
」処
（ ）
その他
program eps_fvm
fuse eps_fvm_all
implicit REAL*8 (A-H,O-Z)
call eps_fvm_input_grid
call poi_gen
call eps_fvm_solver
call
ll output_ucd
t td
end program eps_fvm

メッシュ読み込み(#S-GRID)
マトリクス生成
線形ソルバー
AVS用結果ファイル書き出し(S GRID R UCD)
AVS用結果ファイル書き出し(S-GRID-R-UCD)

T2K-FVM-01

67

program eps_fvm
fuse eps_fvm_all
implicitp
REAL*8 (A-H,O-Z)
call eps_fvm_input_grid
call poi_gen
call
ll eps_fvm_solverflcall output_ucd
end program eps_fvm

メッシュ読み込み(#S-GRID)
マトリクス生成
線形ソルバ
線形ソルバー
AVS用結果ファイル書き出し(S-GRID-R-UCD)

T2K-FVM-01

68

メッシュ関連変数（eps fvm util）
メッシュ関連変数（eps_fvm_util）
変数名

型

配列サイズ

NODE tot
NODE_tot

I

-

要素数

NODE_VOL(:)

R

NODE_tot

要素体積

NODE_COND(:)

R

NODE_tot

要素熱伝導率NODE_XYZ(:)

R

3*NODE_tot

要素重心座標（3次元）

CONN_totI-

コネクティビティ総数

CONN_node(:)

I

2*CONN_tot

コネクティビティ構成要素

CONN_COEF(:)

R

CONN_tot

コネクティビティ係数

FIX_NODE_tot

I-

ディリクレ境界条件適用要素数

FIX_NODE_ID(:)

I

FIX_NODE_tot

ディリクレ境界条件適用要素番号

FIX_NODE_COEF(:)

R

FIX_NODE_tot

ディリクレ境界条件係数

FIX_NODE_VAL(:)

R

FIX_NODE_tot

ディリクレ境界条件値

SURF_NODE_tot

I-

ノイマン境界条件適用要素数

SURF_NODE_ID(:)

I

SURF_NODE_tot

ノイマン境界条件適用要素番号

SURF_NODE_FLUX(:)

R

SURF_NODE_tot

ノイマン境界条件フラックス

BODY_NODE_tot

I-

体積発熱境界条件適用要素数

BODY_NODE_ID(:)

I

BODY_NODE_tot

体積発熱境界条件適用要素番号

BODY NODE FLUX( )BODY_NODE_FLUX(:)

R

BODY NODEBODY_NODE_tot

体積発熱境界条件
体積発熱境界条件フラックスク

内容

T2K-FVM-01

69「eps fvm」処理：メッシュ読み込み（1/5）
「eps_fvm」処理：メッシュ読み込み（1/5）
subroutine eps_fvm_input_grid
use hpcmw_eps_fvm_all
implicit REALREAL*88 (A-H,O-Z)(A H,O Z)
character(len=NAME_LEN) :: member
character(len=80) :: LINE!C!C +------------+!C | MESH INPUT |!C +------------+!C===
open (22, file='fvmmg.ctrl', status='unknown')
read (22,*) NX, NY, NZ
close (22)各座標軸方向の要素数 AVS出力
要素数，AVS出力
用に必要

IUNIT= 11
open (IUNIT,file= ‘fvm_entire_mesh.dat', status='unknown')!C!C!C-NODEread (IUNIT, '(10i10)') NODE_tot

各要素情報の読み込み

if (NODE_tot.ne.NX*NY*NZ) thenwrite (*,'(a)') "incosistent test.grid and test.mesh !!!"
call hpcmw_eps_fvm_abort
endif
allocate (NODE_VOL(NODE_tot), NODE_COND(NODE_tot),&
NODE_XYZ(3*NODE_tot))&

do i= 1,
, NODE_tot
read (IUNIT,'(i10,5e16.6)') ii, NODE_VOL(i), NODE_COND(i),&(NODE_XYZ(3*i-3+k), k=1, 3)
enddo

&

NODE_tot：
総要素数

T2K-FVM-01

70「eps fvm」処理：メッシュ読み込み（1/5）
「eps_fvm」処理：メッシュ読み込み（1/5）!C!C-- NODE
read (IUNIT, '(10i10)')(10i10) ) NODENODE_tottot
if (NODE_tot.ne.NX*NY*NZ) thenwrite (*,'(a)') "incosistent test.grid and test.mesh !!!"
call hpcmw_eps_fvm_abort
endif
allocate (NODE_VOL(NODE_tot), NODE_COND(NODE_tot),
NODE_XYZ(3*NODE_tot))&

do i= 1, NODE_tot
read (IUNIT,'(i10,5e16.6)') ii, NODE_VOL(i), NODE_COND(i),&(NODE(NODE_XYZ(3*i-3+k),
XYZ(3*i 3+k) k=1,
k=1 3)
enddo&&

12 要素数:NODE_tot
11.000000E+00
21.000000E+00
31.000000E+00
41.000000E+00
51.000000E+00
61.000000E+00
71.000000E+00
81.000000E+00
91.000000E+00
101.000000E+00
11
11.000000E+00
000000E+00
121.000000E+00
要素番号 要素体積
NODE_VOL(i)1.000000E+00
1.000000E+00
1.000000E+00
1.000000E+00
1.000000E+00
1.000000E+00
1.000000E+00
1.000000E+00
1.000000E+00
1.000000E+00
11.000000E+00
000000E+00
1.000000E+00
要素熱伝導率NODE_COND(i)5.000000E-01
5.000000E-01
5.000000E-01
1.500000E+00
5.000000E-01
5.000000E-01
5.000000E-01
1.500000E+00
5.000000E-01
1.500000E+00
1.500000E+00
5.000000E-01
5.000000E-01
5.000000E-01
1.500000E+00
1.500000E+00
5.000000E-01
1.500000E+00
5.000000E-01
1.500000E+00
1.500000E+00
1.500000E+00
1.500000E+00
1.500000E+00
5.000000E-01
5.000000E-01
2.500000E+00
1.500000E+00
5.000000E-01
2.500000E+00
5.000000E-01
5 000000E 01
1.500000E+00
1 500000E+00
22.500000E+00
500000E+00
1.500000E+00
1.500000E+00
2.500000E+00
要素中心X座標
要素中心Y座標
要素中心Z座標
NODE_XYZ(3*i-2) NODE_XYZ(3*i-1) NODE_XYZ(3*i)

T2K-FVM-01

71「eps_fvm」処理：メッシュ読み込み（2/5）
コネクティビティ情報!C!C-- CONNECTION
read (IUNIT,'(10i10)')((10 10) ) CONN_tot
allocate (CONN_NODE(2*CONN_tot), CONN_COEF(CONN_tot))
do i= 1, CONN_tot
read (IUNIT,'( 2i10, 3e16.6)') in1, in2, AREA, D1, D2
CONN_NODE(2*i-1)= in1
CONN_NODE(2*i )= in2
C1 = NODE_COND(in1)
C2 = NODE_COND(in2)
CONN_COEF(i)= AREA / ( D1/C1 + D2/C2 )
enddo
20
1
1
1
2
…
6
7
7
8
9
9
10
11
E1

コネクティビティ総数： CONN_tot
21.000000E+00
31.000000E+00
51.000000E+00
411.000000E+00
000000E+00
10
8
11
12
10
11
12
12
E2

1.000000E+00
1.000000E+00
1.000000E+00
11.000000E+00
000000E+00
1.000000E+00
1.000000E+00
1.000000E+00
1.000000E+00
S:要素境界面積

E1= CONN_NODE(2*ic-1)
E2= CONN_NODE(2*ic)

5.000000E-01
5.000000E-01
5.000000E-01
55.000000E
000000E-01
01

5.000000E-01
5.000000E-01
5.000000E-01
5.000000E
5 000000E-01
01

5.000000E-01
5.000000E-01
5.000000E-01
5.000000E-01
5 000000E 01
5.000000E-01
5.000000E-01
5.000000E-01
5.000000E-01
d1：E1重心～
境界面距離5.000000E-01
5.000000E-01
5.000000E-01
5.000000E-01
5 000000E 01
5.000000E-01
5.000000E-01
5.000000E-01
5.000000E-01
d2：E2重心～
境界面距離

CONN_tot：
総コネクティビティ数

T2K-FVM-0172

有限体積法による空間離散化
熱流束に関するつりあい式
隣接要素との熱伝導

∑d
k

Sik
ik

λi

+

d ki

温度固定境界(Tk − Ti ) + ∑
e

λk

∑

Sie
(
TiBe − Ti ) +Sid q&id + Vi Q& i = 0
d ie
d

λi

要素境界面 体積発熱
通過熱流束

a
b
dai
dbi

Sib
Sia

dia

die

TiBe

Sie
Sid

dib

i

dic
dci

Sic
qid

c

λ ：熱伝導率
Vi ：要素体積
S ：表面面積
表面面積
dij ：要素中心から表面までの距離
q ：表面フラックス
Q ：体積発熱
体積発熱
TiB ：境界温度

T2K-FVM-01

73「eps_fvm」処理：メッシュ読み込み（2/5）
コネクティビティ情報!C!C-- CONNECTION!Cread (IUNIT,'(10i10)') CONN_tot
allocate (CONN_NODE(2*CONN_tot), CONN_COEF(CONN_tot))
do i= 1, CONN_tot
read (IUNIT,'( 2i10, 3e16.6)') in1, in2, AREA, D1, D2
CONN_NODE(2*i-1)= in1
CONN NODE(2*i )= in2
CONN_NODE(2*iC1 = NODE_COND(in1)
C2 = NODE_COND(in2)
CONN_COEF(i)= AREA / ( D1/C1 + D2/C2 )enddo

CONN_tot：
総コネクティビティ数（これを2倍すると非対角成
分総数）
非対角成分の計算を実施
している

AREA
D1

D2

in1

in2
λ=C1

∑d
k

S ik
ik

λi

+

i

ik

ki

k

λ=C2

S
(
Tk − Ti ) + ∑ ie (TiBe − Ti ) + ∑ Sid q&id + Vi Q& i = 0
d
d
ki

λk

CONN_COEF

e

ie

λi

d

係数行列の非対角成分

T2K-FVM-01

74「eps_fvm」処理：メッシュ読み込み（3/5）
ディリクレ境界条件!C!C-- DIRICHLET!C
read (IUNIT,'(10i10)') FIX_NODE_tot
allocate (FIX_NODE_ID(FIX_NODE_tot), FIX_NODE_COEF(FIX_NODE_tot))
allocate (FIX_NODE_VAL(FIX_NODE_tot))
do i= 1, FIX_NODE_tot
read (IUNIT,(IUNIT '(i10,'(i10 3e16
3e16.6)')6)')&
FIX_NODE_ID(i), AREA, DIST, FIX_NODE_VAL(i): Tk
icel= FIX_NODE_ID(i)COND= NODE_COND(icel)
FIX_NODE_COEF(i)= AREA / (DIST/COND)
enddo
6 ディリクレ境界条件を与える要素数： FIX_NODE_tot
21.000000E+00
5.000000E-01
0.000000E+00
41.000000E+00
5.000000E-01
0.000000E+00
61.000000E+00
1.000000E 00
5.000000E
5.000000E-01
01
0.000000E
0.000000E+0000
81.000000E+00
5.000000E-01
0.000000E+00
101.000000E+00
5.000000E-01
0.000000E+00
121.000000E+00
5.000000E-01
0.000000E+00
要素番号 境界面面積S
境界面と要素重心距離d 境界値TB
FIX NODE ID(ib)FIX_NODE_ID(ib)FIX
FIX_NODE_VAL(ib)
NODE VAL(ib)

&

∑d
k

ik

λi

+

S
(
Tk − Ti ) + ∑ ie (TiBe − Ti ) + ∑ Sid q&id + Vi Q& i = 0
d
d
ki

λk

e

ie

λi

d

TB DISTTi

FIX_NODE_COEF

∑
e

Sik

AREA

Sie
(
TiBe − Ti )
d ie

λi

FIX_NODE_VAL

T2K-FVM-01

75「eps_fvm」処理：メッシュ読み込み（4/5）
ノイマン境界条件!C!C-- NEUMANN!C
read (IUNIT,'(10i10)') SURF_NODE_tot
allocate& (SURF_NODE_ID (SURF_NODE_tot), SURF_NODE_FLUX(SURF_NODE_tot))&

do i= 1, SURF_NODE_tot__
read (IUNIT, '(i10, 3e16.6)') SURF_NODE_ID(i), AREA, FLUXSURF_NODE_FLUX(i)= AREA*FLUX
enddo
6 ノイマン境界条件を与える要素数： SURF_NODE_tot
11.000000E+00
1.000000E+00
31.000000E+00
1.000000E+00
51.000000E+00
1.000000E+00
71.000000E+00
1.000000E+00
91.000000E+00
1.000000E+00
111.000000E+00
1.000000E+00
要素番号 境界表面積S
表面フラックスQSSURF_NODE_ID(ib)

∑d
k

Sik
ik

λi

+

S
(
Tk − Ti ) + ∑ ie (TiBe − Ti ) + ∑ Sid q&id + Vi Q& i = 0
d
d
ki

λk

e

ie

λi

d

SURF_NODE_FLUX

AREA
FLUX

T2K-FVM-01

76「eps_fvm」処理：メッシュ読み込み（5/5）
体積発熱!C!C-- BODY FLUX
read (IUNIT,'(10i10)') BODY_NODE_tot
allocate (BODY_NODE_FLUX(NODE_tot))
d i
do
i= 1
1, BODY
BODY_NODE_tot
NODE t tread (IUNIT, '(i10, 3e16.6)') icel, FLUXBODY_NODE_FLUX(icel)= FLUX * NODE_VOL(icel)
enddo

∑d
k

Sik
ik

λi

+

4 体積発熱を与える要素数
51.000000E+00
61.000000E+00
71.000000E+00
81.000000E+00
要素番号 体積フラックスQVS(Tk − Ti ) + ∑ ie (TiBe − Ti ) + ∑ Sid q&id + Vi Q& i = 0
d
d
ki

λk

e

ie

λi

d

BODY_NODE_FLUX

FLUX

T2K-FVM-01

77

program eps_fvm
fuse eps_fvm_all
implicitp
REAL*8 (A-H,O-Z)
call eps_fvm_input_grid
call poi_gen
call
ll eps_fvm_solverflcall output_ucd
end program eps_fvm

メッシュ読み込み(#S-GRID)
マトリクス生成
線形ソルバ
線形ソルバー
AVS用結果ファイル書き出し(S-GRID-R-UCD)

T2K-FVM-01

78

マトリクス関連変数表（eps fvm pcg）
マトリクス関連変数表（eps_fvm_pcg）
変数名

型

サイズ

内容
連立一次方程式係数マトリクス非対角成分総数

NPLU

I

-

D(:)

R

NODE_tot

連立一次方程式係数マトリクス対角成分

PHI(:)

R

NODE_tot

連立一次方程式未知数ベクトル

BFORCE(:)

R

NODE_tot

連立一次方程式右辺ベクトル

index(:)

I

0:NODE tot
0:NODE_tot

係数マトリクス非対角成分要素番号用一次元圧縮配列（非対角成分数）
係数マトリクス非対角成分要素番号用
次元圧縮配列（非対角成分数）

item(:)I

NPLU

係数マトリクス非対角成分要素番号用一次元圧縮配列（非対角成分要素（列）番号）

AMAT(:)R

NPLU

係数マトリクス非対角成分要素番号用一次元圧縮配列（非対角成分）
係数マトリクス非対角成分要素番号用
次元圧縮配列（非対角成分）

非零非対角成分のみを格納する
Compressed Row Storage法を使用している。

T2K-FVM-01

79

行列ベクトル積への適用（非零）非対角成分のみを格納，疎行列向け方法（非零）非対角成分のみを格納
疎行列向け方法
Compressed Row Storage (CRS)
D(i) 対角成分（実数対角成分（実数，i=1,N））
index(i) 非対角成分数に関する一次元配列（通し番号）（整数，i=0,N）iitem(k)(k) 非対角成分の要素（列）番号（整数，k=1, index(N)）
AMAT(k) 非対角成分（実数 k 1 i（実数，k=1,
index(N)）
d (N)）{Y}= [A]{X}
do i= 1, N
Y(i)= D(i)*X(i)d k=dok iindex(i-1)+1,
d (i 1) 1 index(i)i d (i)Y(i)= Y(i) + AMAT(k)*X(item(k))
enddo
enddo

T2K-FVM-0180

行列ベクトル積：密行列⇒とても簡単
⎡ a11
⎢ a
⎢ 21
⎢ ...
⎢
⎢a N −1,1⎢⎣ a N ,1

a12
a22

...

a1, N −1
a2, N −1...
a N −1, 2
aN ,2

a N −1, N −1...

a N , N −1

a1, N ⎧⎤ x1 ⎫ ⎧ y1 ⎫a2, N ⎪⎥⎥ x2 ⎪ ⎪ y2 ⎪
⎪⎪
⎪⎪ ⎪⎪
⎪⎪
... ⎨⎥ M ⎬ = ⎨ M ⎬⎥⎪a N −1, N ⎥ x N −1 ⎪ ⎪ y N −1 ⎪⎪⎪ ⎪⎪a N , N ⎪⎩⎥⎦ x N ⎭⎪ ⎩⎪ y N ⎭⎪{Y}= [A]{X}
do j= 1, N
Y(j)= 0.d0
do i= 1,
1 N
Y(j)= Y(j) + A(i,j)*X(i)
enddo
enddo
dd

T2K-FVM-01

81

Compressedp
Row Storageg (CRS)()
12

1 ⎡1.1

2.4

2

3.6

⎢4.3
⎢
3 ⎢ 0
⎢
4
⎢0
5 ⎢ 3.1
⎢
6 ⎢ 0
7 ⎢ 0
⎢
8 ⎢
⎣0

0
4.1
9.5
0
6.4
9.5

3

4

5

6

7

8

0 ⎤
0
2.5 0
3.7
0
9.1 ⎥⎥
5.7
0 1.5
0
3.1
0 ⎥
⎥
0 9.8 2.5 2.7
0
0 ⎥10.4 0 11.5 0
4.3
0 ⎥
⎥
6.5 0
0 12.4 9.5
0 ⎥
2.5 0
0
1.4 23.1 13.1⎥⎥
1.3 9.6 0
3.1
0 51.3⎥⎦
0

0

3.2

0

0

T2K-FVM-01

82

Compressedp
Row Storageg (CRS)
(
)
1
2

1

2

1.1
1
1
①
4.3
①

2.4
2
4
②
3.6
②

4
3.1
①

4.1
4
1
②
9.5
②

6
7
8

6.4
6
4
②
9.5
②

4

5

6

7

8
NODE tot= 8
NODE_tot=

3.2
3
2
⑤
2.5
④
5.7
③

3

5

3

3.7
⑥

9.1
⑧

1.5
3.1
⑤
⑦
9.8
9
8 2.5
25 2
2.7
7
④
⑤
⑥
10.4
11.5
4.3
③
⑤
⑦
6.5
12.4 9.5
③
⑥
⑦
25
2.5
1 4 23.1
1.4
23 1 13.1
13 1
③
⑥
⑦
⑧
1.3 9.6
3.1
51.3
③
④
⑥
⑧

対角成分
D(1)= 1
1.1
1
D(2)= 3.6
D(3)= 5.7
D(4)= 9.8
D(4)
D(5)= 11.5
D(6)= 12.4
D(7)=( ) 23.1
D(8)= 51.3

T2K-FVM-01

83

Compressedp
Row Storageg (CRS)
(
)
1
2
3
4
5
6
7
8

1

2

1.1
1
1
①
3.6 4.3
②
①
5.7
③
98
9.8
④
11.5 3.1
⑤
①
12.4
⑥
23 1
23.1
⑦
51.3
⑧

2.4
2
4
②

3

4

5

6

7

8

3.2
3
2
⑤
2.5
④

4.1
4
1
②
9.5 10.4②③
6.5③
6 4 2.5
6.4
25②③
9.5 1.3 9.6
②
③
④

3.7
⑥
1.5
⑤
25
2.5
⑤

9.1
⑧
3.1
⑦

2.7
2
7
⑥
4.3
⑦
9.5
⑦
1.4
1
4
⑥
3.1
⑥

13.1
13
1
⑧

T2K-FVM-01

84

Compressedp
Row Storageg (CRS)
(
)
非対角
成分数
1
2
3
4
5
6
7
8

1.1
1
1
①
3.6
②
5.7
③
98
9.8
④
11.5
⑤
12.4
⑥
23 1
23.1
⑦
51.3
⑧

2.4
2
4
②
4.3
①
1.5
⑤
41
4.1
②
3.1
①
6.5
③
64
6.4
②
9.5
②

3.2
3
2
⑤
2.5 3.7 9.1
④
⑥
⑧
3.1
⑦
2 5 2.7
2.5
27⑤⑥
9.5 10.4 4.3
②
③
⑦
9.5
⑦
2 5 1.4
2.5
1 4 1313.1
1③⑥⑧
1.3 9.6 3.1③④⑥

index(0)=02index(1)=24index(2)=62index(3)=83index(4)= 11

4index(5)=( ) 152index(6)= 17

4index(7)= 21

4index(8)=( ) 25index(i-1)+1~index(i)番目がi行目の非対角成分

NPLU= 25(=index(N))

T2K-FVM-01

85

Compressedp
Row Storageg (CRS)
(
)
非対角
成分数
1
2
3
4
5
6
7
8

1.1
1
1 2.4
2 4 3.2
32① ②,1 ⑤,2
3.6 4.3 2.5 3.7 9.1② ①,3① 3 ④,4④ 4 ⑥,5⑥ 5 ⑧,6⑧6
5.7 1.5 3.1③ ⑤,7 ⑦,8
9 8 4.1
9.8
4 1 2.52 5 2.7
27④ ②,9 ⑤,10 ⑥,11
11.5 3.1 9.5 10.4 4.3⑤ ①,12① 12 ②,13② 13 ③,14③ 14 ⑦,15⑦ 15
12.4 6.5 9.5⑥ ③,16 ⑦,17
23 1 6.4
23.1
6 4 2.5
2 5 1.4
1 4 13.1
13 1⑦ ②,18 ③,19 ⑥,20 ⑧,21
51.3 9.5 1.3 9.6 3.1② 22 ③,23③ 23 ④,24④ 24 ⑥,25⑥ 25⑧ ②,22index(0)=02index(1)=24index(2)=62index(3)=83index(4)= 11

4index(5)=( ) 152index(6)= 17

4index(7)= 21

4index(8)=( ) 25index(i-1)+1~index(i)番目がi行目の非対角成分

NPLU= 25(=index(N))

T2K-FVM-01

86

Compressedp
Row Storageg (CRS)
(
)
1
2
3
4
5
6
7
8

1.1
1
1 2.4
2 4 3.2
32① ②,1 ⑤,2
3.6 4.3 2.5 3.7 9.1② ①,3① 3 ④,4④ 4 ⑥,5⑥ 5 ⑧,6⑧6
5.7 1.5 3.1③ ⑤,7 ⑦,8
9 8 4.1
9.8
4 1 2.52 5 2.7
27④ ②,9 ⑤,10 ⑥,11
11.5 3.1 9.5 10.4 4.3⑤ ①,12① 12 ②,13② 13 ③,14③ 14 ⑦,15⑦ 15
12.4 6.5 9.5⑥ ③,16 ⑦,17
23 1 6.4
23.1
6 4 2.5
2 5 1.4
1 4 13.1
13 1⑦ ②,18 ③,19 ⑥,20 ⑧,21
51.3 9.5 1.3 9.6 3.1② 22 ③,23③ 23 ④,24④ 24 ⑥,25⑥ 25⑧ ②,22例：
item( 7)7)= 5, AMAT( 7)=7) 1.5item(19)= 3, AMAT(19)= 2.5

T2K-FVM-01

87

Compressedp
Row Storageg (CRS)
(
)
1
2
3
4
5
6
7
8

1.1
1
1 2.4
2 4 3.2
32① ②,1 ⑤,2
3.6 4.3 2.5 3.7 9.1② ①,3① 3 ④,4④ 4 ⑥,5⑥ 5 ⑧,6⑧6
5.7 1.5 3.1③ ⑤,7 ⑦,8
9 8 4.1
9.8
4 1 2.52 5 2.7
27④ ②,9 ⑤,10 ⑥,11
11.5 3.1 9.5 10.4 4.3⑤ ①,12① 12 ②,13② 13 ③,14③ 14 ⑦,15⑦ 15
12.4 6.5 9.5⑥ ③,16 ⑦,17
23 1 6.4
23.1
6 4 2.5
2 5 1.4
1 4 13.1
13 1⑦ ②,18 ③,19 ⑥,20 ⑧,21
51.3 9.5 1.3 9.6 3.1② 22 ③,23③ 23 ④,24④ 24 ⑥,25⑥ 25⑧ ②,22D(i) 対角成分（実数，i=1,N）index(i) 非対角成分数に関する一次元配列（通し番号）（整数 i 0 N）
（通し番号）（整数，i=0,N）
item(k) 非対角成分の要素（列）番号（整数，k=1, index(N)）
AMAT(k) 非対角成分（実数，k=1, index(N)）{Y}={} [A]{X}[ ]{ }
do i= 1, N
Y(i)= D(i)*X(i)d kk= iindex(i-1)+1,do
d (i 1) 1 iindex(i)d (i)Y(i)= Y(i) + AMAT(k)*X(item(k))
enddo
enddo

T2K-FVM-01「eps fvm」処理：マトリクス生成（1/5）
「eps_fvm」処理：マトリクス生成（1/5）!C!C***!C*** POI_GEN
!C
!C***
!C
!C
generate COEF. MATRIX for POISSON equations!C
subroutine POI_GEN
use eps
eps_fvm_all
fvm all
implicit REAL*8 (A-H,O-Z)
integer, pointer :: IWKX(:,:)!C!C +-------+!C | INIT. |!C ++-------+
+
!C===
!C
!C-- MATRIX
nn = NODE_tot
allocate (BFORCE(nn), D(nn), PHI(nn))
allocate (index(0:nn))
BFORCE= 0.d0
PHI= 0.d0
D= 0.d0
index= 088

T2K-FVM-01

89「eps fvm」処理：マトリクス生成（2/5）
「eps_fvm」処理：マトリクス生成（2/5）!C!C-- ETC.
allocate (IWKX(NODE_tot,6))
IWKX= 0do ic= 1, CONN_tot
in1= CONN_NODE(2*ic-1)
in2= CONN_NODE(2*ic )ik1= index(in1) + 1
IWKX (in1,ik1)= ic
index(in1)= ik1
ik2= index(in2) + 1
IWKX (in2,ik2)= ic
i d (i 2
index(in2))= ik2enddo
do i= 1, nn
index(i)= index(i-1) + index(i)enddo
NPLU= index(nn)allocate (item(NPLU), AMAT(NPLU))（非零）非対角成分数算出（非零）非対角成分抽出
index(in)：
各要素の非対角成分数（通し番号）IWKX(in,1-6)：
各要素の非対角成分のコネクティビティID（CONN_NODE）
非対角要素数を増やした場合は
IWKXのサイズを増やせばよい（現在は隣接要素数（現在は隣接要素数が6までに限ま限定されている）T2K-FVM-01

90「eps fvm」処理：マトリクス生成（2/5）
「eps_fvm」処理：マトリクス生成（2/5）!C!C-- ETC.
allocate (IWKX(NODE_tot,6))
IWKX= 0do ic= 1, CONN_tot
in1= CONN_NODE(2*ic-1)
in2= CONN_NODE(2*ic )ik1= index(in1) + 1
IWKX (in1,ik1)= ic
index(in1)= ik1
ik2= index(in2) + 1
IWKX (in2,ik2)= ic
i d (i 2
index(in2))= ik2enddo
do i= 1, nn
index(i)= index(i-1) + index(i)enddoindex(in)：
各要素の非対角成分数（通し番号）

NPLU= index(nn)allocate (item(NPLU), AMAT(NPLU))

NPLU：（非零）非対角成分の総数

T2K-FVM-01

91「eps fvm」処理：マトリクス生成（3/5）
「eps_fvm」処理：マトリクス生成（3/5）
do i= 1, NODE_tot
do j= 1, index(i)-index(i-1)
k= index(i-1) + j
ic = IWKX(i,j)
in1= CONN_NODE(2*ic-1)
in2= CONN_NODE(2*ic )if (in1.eq.i) thenitem(k)= in2
else
item(k)= in1
endif
enddo
enddo

index(in)：
各要素の非対角成分数
各要素
非対角成分数（通し番号）IWKX(in,1-6)：
各要素の非対角成分のコネクティビティID（CONN_NODE）

!C
!C===
!C
!C +----------------------------+!C | INTERIOR NODEs + BODY FLUX |!C +----------------------------+!C===
do icel= 1, NODE_tot
BFORCE(icel)= BFORCE(icel) + BODY_NODE_FLUX(icel)
enddo
do i= 1, NODE_tot
do j= index(i-1)+1, index(i)icon= IWKX(i,j-index(i-1))
AMAT(j)= -CONN_COEF(icon)
D(i)= D(i) + CONN_COEF(icon)
enddo
enddo
deallocate (IWKX)!C===

item(k)：非対角成分（列番号）T2K-FVM-01

92

有限体積法による空間離散化
熱流束に関するつりあい式
隣接要素との熱伝導

∑d
k

Sik
ik

λi

+

d ki

温度固定境界(Tk − Ti ) + ∑
e

λk

∑

Sie
(
TiBe − Ti ) +Sid q&id + Vi Q& i = 0
d ie
d

λi

要素境界面 体積発熱
通過熱流束

a
b
dai
dbi

Sib
Sia

dia

die

TiBe

Sie
Sid

dib

i

dic
dci

Sic
qid

c

λ ：熱伝導率
Vi ：要素体積
S ：表面面積
表面面積
dij ：要素中心から表面までの距離
q ：表面フラックス
Q ：体積発熱
体積発熱
TiB ：境界温度

T2K-FVM-01

93

全体マトリクスの生成
要素iに関する釣り合い

∑d
k

Sik
ik

λi

+

d ki

(Tk − Ti ) + ∑λk

隣接要素との熱伝導

e

∑

Sie
(
TiBe − Ti ) +Sid q&id + Vi Q& i = 0
d ie
d

λi

温度固定境界

要素境界面 体積
通過熱流束 発熱

T2K-FVM-01

94

全体マトリクスの生成
要素iに関する釣り合い
−∑
k

Sik
d ik

λi

+

Sie
(
(
Tk − Ti ) − ∑TiBe − Ti ) − ∑ Sid q&id − Vi Q& i = 0
d ki
d ie
e
d

λk

隣接要素との熱伝導λi

温度固定境界

要素境界面 体積
通過熱流束 発熱

T2K-FVM-01

95

全体マトリクスの生成
要素iに関する釣り合い
−∑
k

−

d ik

λi

∑d
k

Sik
+

Sie
(
(
Tk − Ti ) − ∑TiBe − Ti ) − ∑ Sid q&id − Vi Q& i = 0
d ki
d ie
e
d

λk

Sik
ik

λi

+

d ki

λk

Tk +

λi

∑d
k

Sik
ik

λi

+

d ki

λk

Ti −

∑
e

∑

∑

Sie
Sie
TiBe +Ti =Sid q&id + Vi Q& i
d iei
d iei
e
d

λi

定数項
定数項：右辺へ移項
移項λi

T2K-FVM-01

96

全体マトリクスの生成
要素iに関する釣り合い
−∑
k

−

d ik

λi

∑d
k

Sik
+

Sie
(
(
Tk − Ti ) − ∑TiBe − Ti ) − ∑ Sid q&id − Vi Q& i = 0
d ki
d ie
e
d

λk

Sik
ik

λi

+

d ki

λi

Tk +

∑d
k

λk

Sik
ik

λi

+

d ki

Ti −

λk

∑
e

∑

∑

Sie
Sie
TiBe +Ti =Sid q&id + Vi Q& i
d iei
d iei
e
d

λi

λi

⎡
⎤
⎡
⎤
⎢
⎢
⎥
Sik
Sie ⎥
Sik
Sie
&
⎢
⎥Ti − ⎢+
Tk ⎥ =
Sid q&id + Vi Qi +TiBe
d ie ⎥d ie⎢ k d ik + d ki
⎢ k d ik + d ki ⎥ d
e
e
⎢⎣
⎢⎣
λi λk
λi ⎥⎦
λi λk ⎥⎦λi∑∑

D（対角成分）∑AMAT（非対角成分）∑∑

BFORCE（右辺）T2K-FVM-01

97「eps_fvm」処理：マトリクス生成（3/5）
体積発熱
do i= 1, NODE_tot_
do j= 1, index(i)-index(i-1)
k= index(i-1) + j
ic = IWKX(i,j)
in1= CONN_NODE(2*ic-1)
in2= CONN_NODE(2*ic )if (in1.eq.i) thenitem(k)= in2
else
item(k)= in1
endif
e
d
enddo
enddo!C===

⎡
⎤
⎡
⎤
⎢
⎢
⎥
Sik
Sik
Sie ⎥
⎢
⎥Ti − ⎢
+
Tk ⎥
d
d
d
d
d
ie ⎥
⎢ k ik + ki
⎢ k ik + ki ⎥
e
⎢⎣
⎢⎣
λi λk
λi ⎥⎦
λi λk ⎥⎦
S ie
&
&
=
TiBe
Sid qid + Vi Qi +die
BFORCE（右辺）
（右辺）
d
e

∑

∑

∑

!C
BODY_NODE_FLUX!C +----------------------------+!C | INTERIOR NODEs + BODY FLUX |!C +----------------------------+!C===
do icel= 1, NODE_tot
BFORCE(icel)= BFORCE(icel) + BODY_NODE_FLUX(icel)
enddo
do i= 1, NODE_tot
do j= index(i-1)+1, index(i)icon= IWKX(i,j-index(i-1))
AMAT(j)= -CONN_COEF(icon)
D(i)= D(i) + CONN_COEF(icon)
enddo
enddo
deallocate (IWKX)
!C===

∑

∑

λi

T2K-FVM-01

98

「eps
p _fvm」処理：マトリクス生成（3/5）
内部熱伝導
do i= 1,
, NODE_tot
do j= 1, index(i)-index(i-1)
k= index(i-1) + j
ic = IWKX(i,j)
in1= CONN_NODE(2*ic-1)
in2= CONN_NODE(2in2CONN NODE(2*ic
ic )if (in1.eq.i) thenitem(k)= in2
else
item(k)= in1
endif
enddo
enddo!C===D（対角成分）

⎡
⎤
⎡
⎤
⎢
⎢
⎥
S ik
Sie ⎥
Sik
⎢
⎥Ti − ⎢+
Tk ⎥d ie ⎥⎢ k d ik + d ki
⎢ k d ik + d ki ⎥
e
⎢⎣
λi λk
λi ⎥⎦
λi λk ⎥⎦
⎣⎢

∑

∑

=

!C===

∑

CONN_COEF∑

S q& + Vi Q& i +

!C
id id
!C +----------------------------+d!C | INTERIOR NODEs + BODY FLUX |!C +----------------------------+!C===
do icel= 1, NODE_tot
BFORCE(icel)= BFORCE(icel) + BODY_NODE_FLUX(icel)
enddo
do i= 1, NODE_tot
do j= index(i-1)+1, index(i)icon= IWKX(i,j-index(i-1))
AMAT(j)= -CONN_COEF(icon)
D(i)= D(i) + CONN_COEF(icon)
enddo
enddo
deallocate (IWKX)AMAT（非対角成分）

CONN_COEF

∑
e

S ie
TiBe
d ie

λi

T2K-FVM-01

99「eps_fvm」処理：マトリクス生成（4/5）
ディリクレ境界条件（表面温度固定）!C!C +-----------+!C | DIRICHLET |!C +-----------+!C===
do i= 1, FIX_NODE_tot
icel= FIX_NODE_ID(i)D(icel)=() D(icel)() + FIX_NODE_COEF(i)( )BFORCE(icel)= BFORCE(icel) + FIX_NODE_COEF(i)*FIX_NODE_VAL(i)
enddo
!C===

⎡
⎤
⎡
⎤
D（対角成分）
⎢
⎢
⎥
Sik
Sie ⎥
Sik
⎢
⎥Ti − ⎢+
Tk ⎥d ie ⎥⎢ k d ik + d ki
⎢ k d ik + d ki ⎥
e
⎢⎣
⎢⎣
λi λk
λi ⎥⎦
λi λk ⎥⎦∑∑∑FIX_NODE_COEF=∑d

Sid q&id + Vi Q& i +

∑
e

S ie
TiBe BFORCE（右辺）
d ie

λi

FIX_NODE_COEF * FIX_NODE_VAL

TB
TiT2K-FVM-01

100「eps_fvm」処理：マトリクス生成（5/5）
ノイマン境界条件（表面熱流束）!C!C +--------------+!C | SURFACE FLUX |!C +--------------+!C===
do i= 1, SURF_NODE_tot
icel= SURF_NODE_ID(i)( )BFORCE(icel)= BFORCE(icel) + SURF_NODE_FLUX(i)
enddo
!C===
return
end

⎡
⎤
⎡
⎤
⎢
⎢
⎥
S ik
Sie ⎥
Sik
⎢
⎥Ti − ⎢
+
Tk ⎥
d
d
d
d
d
ie ⎥
⎢ k ik + ki
⎢ k ik + ki ⎥
e
⎢⎣
⎢⎣
λi λk
λi ⎥⎦
λi λk ⎥⎦
S ie=
Sid q&id + Vi Q& i +TiBe BFORCE（右辺）
d ie
d
e

∑

∑

∑

SURF_NODE_FLUX

∑

∑

λi

T2K-FVM-01

101

program eps_fvm
fuse eps_fvm_all
implicitp
REAL*8 (A-H,O-Z)
call eps_fvm_input_grid
call poi_gen
call eps_fvm_solver
eps fvm solver
call output_ucd
end program eps_fvm

メッシュ読み込み(#S-GRID)
マトリクス生成
線形ソルバ
線形ソルバー
AVS用結果ファイル書き出し(S-GRID-R-UCD)

T2K-FVM-01

102「eps fvm」処理：ソルバー「eps_fvm」処理：ソルバ
subroutine eps_fvm_solver
use eps_fvm_all__
implicit REAL*8 (A-H,O-Z)
EPS = 1.d-8
ITR = NODE_tot
call eps
p _fvm_solver_CG&( NODE_tot, NPLU, D, BFORCE, PHI, EPS,&
ITR, IER, index, item, AMAT, COMMtime)
ISET= ISET + 1
open (11, file=‘fvmmg.ctrl', status='unknown')
read (11,(11,*)) NX, NY, NZ
close (11)iS= NX*NY*NZ/2 + NX*NY/2
do i= iS+1, iS+NX
write (*,'(i8,3(1pe16.6))') i, PHI(i)enddo
end subroutine eps_fvm_solver

&
&

103

科学技術計算における大規模線形方程式の解法• 多くの科学技術計算は，最終的に大規模線形方程式
Ax=bを解くことに帰着される。
– important, expensive• アプリケーションに応じて様々な手法が提案されている
アプリケ ションに応じて様々な手法が提案されている– 疎行列（sparse），密行列（dense）– 直接法（direct），反復法（iterative）
直接法（direct） 反復法（iterative）• 密行列（dense）– グローバルな相互作用：BEM，スペクトル法，MO，MD（気液）
グロ バルな相互作用：BEM スペクトル法 MO MD（気液）• 疎行列（sparse）– ローカルな相互作用：FEM，FDM，MD（固），高速多重極展開付
カ な相互作用（固） 高速多重極展開付
BEM

104

直接法（Direct Method）• Gaussの消去法，完全LU分解– 逆行列A-1を直接求める• 利点– 安定，幅広いアプリケ
安定 幅広いアプリケーションに適用可能
ションに適用可能• Partial Pivoting– 疎行列，密行列いずれにも適用可能• 欠点– 反復法よりもメモリ，計算時間を必要とする• 密行列の場合，O（N3 ）の計算量– 大規模な計算向けではない• O（N2 ）の記憶容量，O（N3 ）の計算量

105

反復法（Iterative Method）• 定常（stationary）法– 反復計算中，解ベクトル以外の変数は変化せず– SOR，Gauss
SOR Gauss-Seidel
Seidel，Jacobiなど
Jacobiなど– 概して遅い• 非定常（nonstationary）法– 拘束，最適化条件が加わる
拘束 最適化条件が加わる– Krylov部分空間（subspace）への写像を基底として使用するため Krylov部分空間法とも呼ばれる
め，Krylov部分空間法とも呼ばれる– CG（Conjugate Gradient：共役勾配法）– BiCGSTAB（Bi-Conjugate
BiCGSTAB（Bi Conjugate Gradient Stabilized）– GMRES（Generalized Minimal Residual）

106

反復法（Iterative Method）（続き）• 利点– 直接法と比較して，メモリ使用量，計算量が少ない。
– 並列計算には適している。

• 欠点– 収束性が，アプリケーション，境界条件の影響を受けやすい。
– 前処理（preconditioning）が重要。

T2K-FVM-01

107

代表的な反復法：共役勾配法• Conjugate Gradient法，略して「CG」法– 最も代表的な「非定常」反復法• 対称正定値行列（Symmetric Positive Definite：SPD）
向け– 任意のベクトル{x}に対して{x}T[A]{x}>0– 全対角成分>０，全固有値>0，全部分行列式>0と同値
全対角成分>０ 全固有値>0 全部分行列式>0と同値– 熱伝導，弾性，ねじり・・・本コードの場合もSPD• アルゴリズム– 最急降下法（Steepest Descent Method）の変種– x(i)= x(i-1) + αip(i)• x(i)：反復解，p(i)：探索ベクトル，αi：定数）– 厳密解をyとするとき {x-y}T[A]{x-y}を最小とするような{x}を求める。
– 詳細は例えば：森正武「数値解析（第2版）」（共立出版）

T2K-FVM-01108

共役勾配法のアルゴリズム
Compute r(0)= bb-[A]x[A]x(0)for i= 1, 2, …
z(i-1)= r(i-1)ρii-11= r(i-1) z(i-1)
if i=1p(1)= z(0)else
βi-1= ρi-1/ρi-2p(i)= z(i-1) + βi-1 p(i-1)
endifq(i)= [A]p(i)αi = ρi-1/p(i)q(i)x(i)= x(i-1) + αip(i)(i 1) - α q(i)r(i)= r(i-1)
icheck convergence |r|end• 行列ベクトル積• ベクトル内積• ベクトル定数倍の加減x(i) ：ベクトル
αi ：スカラー

T2K-FVM-01109

共役勾配法のアルゴリズム
Compute r(0)= bb-[A]x[A]x(0)for i= 1, 2, …
z(i-1)= r(i-1)ρii-11= r(i-1) z(i-1)
if i=1p(1)= z(0)else
βi-1= ρi-1/ρi-2p(i)= z(i-1) + βi-1 p(i-1)
endifq(i)= [A]p(i)αi = ρi-1/p(i)q(i)x(i)= x(i-1) + αip(i)(i 1) - α q(i)r(i)= r(i-1)
icheck convergence |r|end• 行列ベクトル積• ベクトル内積• ベクトル定数倍の加減x(i) ：ベクトル
αi ：スカラー

T2K-FVM-01110

共役勾配法のアルゴリズム
Compute r(0)= bb-[A]x[A]x(0)for i= 1, 2, …
z(i-1)= r(i-1)ρi-1= r(i-1) z(i-1)
if i=1p(1)= z(0)else
βi-1= ρi-1/ρi-2p(i)= z(i-1) + βi-1 p(i-1)
endifq(i)= [A]p(i)αi = ρi-1/p(i)q(i)x(i)= x(i-1) + αip(i)(i 1) - α q(i)r(i)= r(i-1)
icheck convergence |r|end• 行列ベクトル積• ベクトル内積• ベクトル定数倍の加減x(i) ：ベクトル
αi ：スカラー

T2K-FVM-01111

共役勾配法のアルゴリズム
Compute r(0)= bb-[A]x[A]x(0)for i= 1, 2, …
z(i-1)= r(i-1)ρii-11= r(i-1) z(i-1)
if i=1p(1)= z(0)else
βi-1= ρi-1/ρi-2p(i)= z(i-1) + βi-1 p(i-1)
endifq(i)= [A]p(i)αi = ρi-1/p(i)q(i)x(i)= x(i-1) + αip(i)(i 1) - α q(i)r(i)= r(i-1)
icheck convergence |r|end• 行列ベクトル積• ベクトル内積• ベクトル定数倍の加減（DAXPY）

x(i) ：ベクトル
αi ：スカラー

T2K-FVM-01

112

前処理（preconditioning）とは?
• 反復法の収束は係数行列の固有値分布に依存– 固有値分布が少なく，かつ1に近いほど収束が早い（単位行列）– 条件数（condition number）（対称正定）＝最大最小固有値比• 条件数が1に近いほど収束しやすい• もとの係数行列[A]に良く似た前処理行列[M]を適用す
る とることによって固有値分布を改善する。
固有値分布を改善する– 前処理行列[M]によって元の方程式[A]{x}={b}を-1
1[A]，[A’]{x’}={b’}へと変換する。ここで[A’]=[M]と変換する ここで{b’}=[M]-1{b} である。
– [A’]=[M][A’] [M]-11[A]が単位行列に近ければ良いことになる。
[A]が単位行列に近ければ良いことになる[A’]= [A][M]-1のように右から乗ずることもある• 「前処理」は密行列，疎行列ともに使用するが，普通は疎「前処理」は密行列 疎行列ともに使用するが 普通は疎行列を対象にすることが多い。

T2K-FVM-01

113

前処理付き共役勾配法のアルゴリズム
Compute r(0)= b-[A]x(0)for i= 1, 2, …
solve [M]z(i-1)= r(i-1)(i 1) z(i-1)(i 1)
ρi-1= r(i-1)
if i=1p(1)= z(0)xelse
βi-1= ρi-1/ρi-2p(i)= z(i-1) + βi-1 p(i-1)
endifq(i)= [A]p(i)αi = ρi-1/p(i)q(i)x(i)= x(i-1) + αip(i)r(i)= r(i-1) - αiq(i)check convergence |r|end

T2K-FVM-01

114

前処理付共役勾配法
Preconditioned Conjugate Gradient Method （PCG）
Compute r(0)= bb-[A]x[A]x(0)for i= 1, 2, …
solve [M]z(i-1)= r(i-1)ρii-11= r(i-1) z(i-1)
if i=1p(1)= z(0)else
βi-1= ρi-1/ρi-2p(i)= z(i-1) + βi-1 z(i-1)
endifq(i)= [A]p(i)αi = ρi-1/p(i)q(i)x(i)= x(i-1) + αip(i)(i 1) - α q(i)r(i)= r(i-1)
icheck convergence |r|end

実際にやるべき計算は実際にやるべき計算は：{z} = [M ]−1{r}「近似逆行列」の計算が必要：[M ] ≈ [A] , [M ] ≈ [A]−1−1

究極の前処理：本当の逆行列[M ]−1 = [A]−1 , [M ] = [A]対角スケーリング：簡単＝弱い[M ] = [D] , [M ] = [D]−1−1

T2K-FVM-01

115

ILU(0), IC(0)• 最もよく使用されている前処理（疎行列用）
最もよく使用される前処理（疎行列用）– 不完全LU分解• Incomplete LU Factorization– 不完全コレスキー分解• Incomplete Cholesky Factorization（対称行列）• 不完全な直接法– もとの行列が疎でも，逆行列は疎とは限らない。
– fill-in– もとの行列と同じ非ゼロパターン（fill-in無し）を持っているのがILU（0），IC（0）T2K-FVM-01

116

対角スケーリング
対角スケ
リング，点ヤコビ前処理
点ヤコビ前処理• 前処理行列として，もとの行列の対角成分のみを取り出した行列を前処理行列 [M] とする。
– 対角スケーリング，点ヤコビ（point-Jacobi）前処理
⎡ D1
⎢0
⎢
[M ] = ⎢ ...
⎢
⎢0
⎢⎣ 0

0

...

D2

0
0

...
0
0

DN −1
...

0

0 ⎤
0 ⎥⎥
... ⎥
⎥
0 ⎥
DN ⎥⎦

• solve [M]z(i-1)= r(i-1)という場合に逆行列を簡単に求めることができる。
• 簡単な問題では収束する。

T2K-FVM-01

117「eps fvm」処理：ソルバー（1/7）
「eps_fvm」処理：ソルバ
（1/7）
!C
!C***
!C*** CG!C***!C
subroutine eps_fvm_solver_CG&( N, NPLU, D, B, X, EPS, ITR, IER,&
index, item, COEF)
use eps_fvm_util
implicit REAL*8 (A-H,O-Z)real(kind=kreal), dimension(N) :: Dreal(kind=kreal), dimension(N) :: Breal(kind=kreal),
real(kind
kreal), dimension(N) :: X
integer, dimension(0:N) :: index
integer, dimension(NPLU):: item
real(kind=kreal), dimension(NPLU):: COEF
real(kind=kreal)
real(kind
kreal) :: EPS
integer :: ITR, IER
integer :: P, Q, R, Z, DD
real(kind=kreal), dimension(:,:), allocatable, save :: W&&

T2K-FVM-01

118「eps fvm」処理：ソルバー（2/7）
「eps_fvm」処理：ソルバ（2/7）!C!C +-------+!C | INIT. |!C +-------+!C===if (.not.allocated(W)) thenallocate (W(N,4))
endif
X= 0.d0
W= 0.d0
R = 1
Z = 2
Q = 2
P = 3
DD= 4
do i= 1, N
W(i,DD)= 1.0D0 / D(i)enddo
IER = 0Tcomm= 0.d0

W(i,1)= W(i,R) ⇒ {r}W(i 2)= W(iW(i,2)=
W(i,Z)Z) ⇒ {z}W(i,2)= W(i,Q) ⇒ {q}W(i,3)= W(i,P) ⇒ {p}!C===

W(i,4)= W(i,DD) ⇒ 1/DIAG

Compute r(0)= b-[A]x(0)
for i
i= 1,
1 2
2, …
solve [M]z(i-1)= r(i-1)ρi-1= r(i-1) z(i-1)
if i
i=1
1
p(1)= z(0)else
βii-11= ρii-11/ρ
ρii-22p(i)= z(i-1) + βi-1 p(i-1)
endifq(i)= [A]p(i)αi = ρi-1/p(i)q(i)x(i)= x(i-1) + αip(i)r(i)= r(i-1) - αiq(i)check convergence |r|end

T2K-FVM-01

119「eps fvm」処理：ソルバー（2/7）
「eps_fvm」処理：ソルバ（2/7）!C!C +-------+!C | INIT. |!C +-------+!C===if (.not.allocated(W)) thenallocate (W(NP,4))
endif
X= 0.d0
W= 0.d0
R = 1
Z = 2
Q = 2
P = 3
DD= 4
do i= 1, N
W(i,DD)= 1.0D0 / D(i)enddo
IER = 0Tcomm= 0.d0!C===

対角成分の逆数（前処理用）
その都度，除算をすると効率が悪いため，予め配列に格納

Compute r(0)= b-[A]x(0)
for i
i= 1,
1 2
2, …
solve [M]z(i-1)= r(i-1)ρi-1= r(i-1) z(i-1)
if i
i=1
1
p(1)= z(0)else
βii-11= ρii-11/ρ
ρii-22p(i)= z(i-1) + βi-1 p(i-1)
endifq(i)= [A]p(i)αi = ρi-1/p(i)q(i)x(i)= x(i-1) + αip(i)r(i)= r(i-1) - αiq(i)check convergence |r|end

T2K-FVM-01

120「eps fvm」処理：ソルバー（3/7）
「eps_fvm」処理：ソルバ（3/7）!C!C +-----------------------+!C | {r0}= {b} - [A]{xini} |!C +-----------------------+!C===
do i= 1, N
W(i,R) = D(i)*X(i)do jj= index(i-1)+1,() , index(i)( )W(i,R) = W(i,R) + COEF(j) * X(item(j))
enddo
enddo
BNRM2= 0.0D0
do i=
i 1, NBNRM2= BNRM2 + B(i) **2
W(i,R)= B(i) - W(i,R)
enddo
!C===

Compute r(0)= b-[A]x(0)
for i
i= 1,
1 2
2, …
solve [M]z(i-1)= r(i-1)ρi-1= r(i-1) z(i-1)
if i
i=1
1
p(1)= z(0)else
βii-11= ρii-11/ρ
ρii-22p(i)= z(i-1) + βi-1 p(i-1)
endifq(i)= [A]p(i)αi = ρi-1/p(i)q(i)x(i)= x(i-1) + αip(i)r(i)= r(i-1) - αiq(i)check convergence |r|end

T2K-FVM-01

121「eps fvm」処理：ソルバー（3/7）
「eps_fvm」処理：ソルバ（3/7）!C!C +-----------------------+!C | {r0}= {b} - [A]{xini} |!C +-----------------------+!C===
do i= 1, N
W(i,R) = D(i)*X(i)do jj= index(i-1)+1,() , index(i)( )W(i,R) = W(i,R) + COEF(j) * X(item(j))
enddo
enddo
BNRM2= 0.0D0
do i=
i 1, NBNRM2= BNRM2 + B(i) **2
W(i,R)= B(i) - W(i,R)
enddo
!C===

BNRM2=|b|2
あとで収束判定に使用

Compute r(0)= b-[A]x(0)
for i
i= 1,
1 2
2, …
solve [M]z(i-1)= r(i-1)ρi-1= r(i-1) z(i-1)
if i
i=1
1
p(1)= z(0)else
βii-11= ρii-11/ρ
ρii-22p(i)= z(i-1) + βi-1 p(i-1)
endifq(i)= [A]p(i)αi = ρi-1/p(i)q(i)x(i)= x(i-1) + αip(i)r(i)= r(i-1) - αiq(i)check convergence |r|end

T2K-FVM-01

122「eps fvm」処理：ソルバー（4/7）
「eps_fvm」処理：ソルバ（4/7）!C!C*********************************
do L= 1, ITR!C!C +----------------+!C | {z}= [Minv]{r} |!C +----------------+!C===
do i= 1,, N
W(i,Z)= W(i,DD) * W(i,R)
enddo
!C===
!C
!C +-------------+++!C | RHO= {r}{z} |!C +-------------+!C===
RHO= 0.d0
do i= 1, N
RHO= RHO + W(iRHOW(i,R)R)*W(i
W(i,Z)
Z)
enddo
!C===

ITERATION

Compute r(0)= b-[A]x(0)
for i
i= 1,
1 2
2, …
solve [M]z(i-1)= r(i-1)ρi-1= r(i-1) z(i-1)
if i
i=1
1
p(1)= z(0)else
βii-11= ρii-11/ρ
ρii-22p(i)= z(i-1) + βi-1 p(i-1)
endifq(i)= [A]p(i)αi = ρi-1/p(i)q(i)x(i)= x(i-1) + αip(i)r(i)= r(i-1) - αiq(i)check convergence |r|end

T2K-FVM-01

123「eps fvm」処理：ソルバー（5/7）
「eps_fvm」処理：ソルバ（5/7）!C!C +-----------------------------+!C | {p} = {z} if
ITER=1|!C | BETA= RHO / RHO1 otherwise |!C +-----------------------------+!C===
if ( L.eq.1 ) thendo i= 1, N
W(i,P)=( , ) W(i,Z)
( , )
enddo
else
BETA= RHO / RHO1
do i= 1, N
W(i,P)= W(i,Z) + BETA*W(i,P)
enddo
endif
!C===
!C
!C +-------------+!C | {q}={q} [A]{p} |!C +-------------+!C===
do i= 1, N
W(i,Q) = D(i) * W(i,P)
do j= index(i-1)+1, index(i)W(i Q) = W(i,Q)
W(i,Q)W(i Q) + COEF(j) * W(item(j)W(item(j),P)
P)
enddo
enddo
!C===

Compute r(0)= b-[A]x(0)
for i
i= 1,
1 2
2, …
solve [M]z(i-1)= r(i-1)ρi-1= r(i-1) z(i-1)
if i
i=1
1
p(1)= z(0)else
βii-11= ρii-11/ρii-22
p(i)= z(i-1) + βi-1 p(i-1)
endifq(i)= [A]p(i)αi = ρi-1/p(i)q(i)x(i)= x(i-1) + αip(i)r(i)= r(i-1) - αiq(i)check convergence |r|end

T2K-FVM-01

124「eps fvm」処理：ソルバー（6/7）
「eps_fvm」処理：ソルバ（6/7）!C!C +---------------------+!C | ALPHA= RHO / {p}{q} |!C +---------------------+!C===
C1= 0.d0
do i= 1, N
C1= C1 + W(i,P)*W(i,Q)
enddo
ALPHA= RHO / C1
!C===
!C +----------------------+!C | {x}= {x} + ALPHA*{p} |!C | {r}={r} {r} - ALPHA
ALPHA*{q}{q} |!C +----------------------+!C===
do i= 1, N
X(i) = X(i)+ ALPHA * W(i,P)
W(i,R)= W(i,R) - ALPHA * W(i,Q)
enddo

Compute r(0)= b-[A]x(0)
for i
i= 1,
1 2
2, …
solve [M]z(i-1)= r(i-1)ρi-1= r(i-1) z(i-1)
if i
i=1
1
p(1)= z(0)else
βii-11= ρii-11/ρ
ρii-22p(i)= z(i-1) + βi-1 p(i-1)
endifq(i)= [A]p(i)αi = ρi-1/p(i)q(i)x(i)= x(i-1) + αip(i)r(i)= r(i-1) - αiq(i)check convergence |r|end

T2K-FVM-01

125「eps fvm」処理：ソルバー（7/7）
「eps_fvm」処理：ソルバ（7/7）
DNRM2 = 0.0
do i= 1, N
DNRM2= DNRM2 + W(i,R)**2
enddo
RESID= dsqrt(DNRM2/BNRM2)
1000

if (my_rank.eq.0) write (*, 1000) L, RESID
format ((i5,, 1pe16.6)
p)
if ( RESID.le.EPS) goto 900
RHO1 = RHO
enddo
IER = 1900 continue
ITR= LEPS RESID
EPS=

r= b-[A]x
DNRM2=|r|
2 | |2
BNRM2=|b|2
RESID |r|/|b|
RESID=

returnend subroutine hpcmw_eps_fvm_solver_CG

Compute r(0)= b-[A]x(0)
for i
i= 1,
1 2
2, …
solve [M]z(i-1)= r(i-1)ρi-1= r(i-1) z(i-1)
if i
i=1
1
p(1)= z(0)else
βii-11= ρii-11/ρ
ρii-22p(i)= z(i-1) + βi-1 p(i-1)
endifq(i)= [A]p(i)αi = ρi-1/p(i)q(i)x(i)= x(i-1) + αip(i)r(i)= r(i-1) - αiq(i)check convergence |r|end

T2K-FVM-01

126

いろいろなメッシュでやってみよう• MicroAVSの使い方に慣れる。
c o Sの使い方に慣れる。