MPIによる並列アプリケーション
開発法入門（１）

中島研吾、伊東 聰、大島聡史
東京大学情報基盤センター

T2K-FVM-01

1

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

2

T2K-FVM-01

MPI応用編の方針• まずは1CPU用のシリアルアプリケーションコードを徹底理解– アルゴリズムのやさしい有限体積法を選択– 熱伝導方程式• アプリケーションに適した並列計算のためのデータ構造（＝局所分散データ構造）を考える• 慣れているFORTRANで解説していますが，C言語版のアプリケーションも準備，ほとんど各行解説するので，余り支障は無いでしょう• MPIそのものについても適宜説明しますT2K-FVM-01

3• 有限体積法の概要• １CPU用プログラム eps_fvm– メッシュ生成，形状データ– 計算実行例– プログラムの内容

T2K-FVM-01

4

対象とするアプリケーションの概要• 支配方程式：三次元定常熱伝導方程式– 物性の温度依存性無し，等方性  T    T    T       Q  0
x  x  y  y  z  z 

• 有限体積法（Finite Volume Method，FVM）による空間
離散化– 任意形状の要素，要素中心で変数を定義。
– 直接差分法（Direct Finite Difference Method）とも呼ばれる。

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

T2K-FVM-01

5

一次元熱伝導方程式：中央差分
熱伝導率一定d  dT  Q  0dx  dx i

Ti-1

Ti

x

Ti+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 T  2Ti  Ti 1 
   i 1

2
x



T2K-FVM-01

6

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

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

7

一次元差分法との比較（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
Qsib  

Tb  Ti
T T Sib   b i  Sib
d ib d bi
d ib  d biTb  Ti
Qsib   Sib
d ib d bi


i

b





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

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

T2K-FVM-01

8

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

h

i

a
dai
a

Sia

dia

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

b
dib

dbi

i

b

Sib

厚さ：1

d
k

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

1
Vi

Sik
ik

i



d
k


Tk  Ti    Sid qid  Vi Q i  0
d
ki

k

Sik
ik

i

d



1

Tk  Ti    Sid qid  Q i  0
d
ki

k

Vi

d

T2K-FVM-01

9

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

h

i

a
dai
a

dia

Sia

熱伝導率一様として

1
Vi


k

b
dib

i

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

dbi
Sib

b
厚さ：1

1
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




T

T

T

T

 Tk  Ti 
k
i
k
i
2
2
2
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

10

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

T2K-FVM-01

11

ガウスの定理：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)，V(x,y,z)，W(x,y,z)• 曲面S上で外向きに引いた法線ベクトルn– nx，ny，nz：方向余弦

SV

T2K-FVM-0112ベクトル表記すると• ガウスの定理
T


w
dV

w

 n dS
VS• グリーンの定理

Vvu dV 

T



v

u
n
dS



 vu  dVTSVT2K-FVM-0113

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

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

14

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

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

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

T2K-FVM-01

15

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

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

SB

i

TB

Ti

dB

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

SBTB  Ti qd B / i「薄い」要素を設定（～10-20）

T2K-FVM-0116

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

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-0117

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

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-0118係数行列は疎行列（Sparse Matrix）• 0が多い• 通常「数値解析」などで扱うのは「密行列」• 非ゼロ成分のみを記憶する方法– Compressed Row Storage, CRS（後述）• 密行列と比較して取り扱いが困難

19

疎行列：0が多い• A(i,j)のように正方行列の全成分を記憶することは疎行列では非効率的– 「密」行列向け

D
X




X

X

















  1   F1     F 
D X
X X X
 2   2    3   F3 
X D X
X X X
   
X D
X X
   4   F4    5   F5 
X
D X
X X
   
X X
X D X
X X X
   6   F6     F 
X X X
X D X
X X X
 7   7 
X X
X D
X X
   8   F8      F 
X X
D X
X X
 9   9 
 10   F10 
X X X
X D X
X X X
   
X X X
X D X
X X X  11   F11 
X X
X D
X X  12   F12    
X X
D X
 13   F13     F 
X X X
X D X
  14   14 
X X X
X D X  15   F15    X X
X D  16   F16 X

X

X• 有限体積法：非零非対角成分の数は高々「数百」規模– 例えば未知数が108個あるとすると記憶容量（ワード数）は• 正方行列：O(1016)• 非零非対角成分数：O(1010)• 非零成分のみ記憶するのが効率的
1D-Part1

T2K-FVM-01

20

現行コードにおける制限• 要素形状は六面体のみ。
– 各要素は各辺長さ1の立方体– 全体形状，境界条件等はメッシュジェネレータで規定• 隣接要素数は要素各面に接する6面までとする。
– プログラムを変更すれば増やすことは可能• 境界条件– ディリクレ境界条件– ノイマン境界条件（面フラックス）– 体積発熱• 要素番号は1から始まって連続していなければならない。

T2K-FVM-01

21

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

NZ

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

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

NY

z

y

NX
ｘ

T2K-FVM-01

22• 有限体積法の概要• １CPU用プログラム eps_fvm– メッシュ生成，形状データ– 計算実行例– プログラムの内容

T2K-FVM-01

23

T2Kで利用できるコンパイラ• 日立コンパイラ（デフォルト）• Intel
• PGI
• GNU

• 本講義では「日立コンパイラ」を使用する– mpif90 -Oss –noparallel s1-3.f– mpicc -Os –noparallel s1-3.c

T2K-FVM-01

24

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

T2K-FVM-01

25

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

T2K-FVM-01

26

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

T2K-FVM-01

27.bash_profile$> cd*****************************************$> cp /home/t00007/.bash_profile .
または ~/.bash_profileに以下の2行を挿入
PATH="$PATH":/home/t00007/bin
export PATH******************************************$> source .bash_profile$> printenv

PATHに以下が入っていることを確認
home/t00007/bin
次回からはログイン時に自動的に起動

T2K-FVM-01

28

メッシュ生成からやってみよう$> cd <$FVM>/run$> cat fvmmg.ctrl
2 2 3$> eps_fvm_mg$> ls fvm_entire*
fvm_entire_mesh.dat
fvm_entire_mesh.inp
fvm_entire_mesh.inp_geo• 入力– fvmmg.ctrl：各座標軸方向の要素数 （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

29

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

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

メッシュジェネレータ
eps_fvm_mg#S-GRID<$FVM>/run/fvm_entire_mesh.dat#S-GRID-UCD<$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

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

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
1.000000E+00

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
1.000000E+00
1.000000E+00

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
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

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

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
0.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

各要素情報

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

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

30

T2K-FVM-01

31

変数表
変数名

型

配列サイズ

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(:)

R

BODY_NODE_tot

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

内容T2K-FVM-01

32

初期全体メッシュデータ（例）（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
101.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
1.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
1.500000E+00
5.000000E-01
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= 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

6
2

8

4

&
&

T2K-FVM-01

33

初期全体メッシュデータ（例）（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
41.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
E2
S:要素境界面積

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

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
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
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
d2：E2重心～
境界面距離

d1

d2

E1

E2

i

ik

ki

k

係数行列の非対角成分

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

&
&

T2K-FVM-01

34

初期全体メッシュデータ（例）（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-01
1.0E+00
5.0E-01
S:要素境界面積 d1：E1重心～
境界面距離

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

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

9
5
1

11

12
10

7

8

5

4

1

3

6
2

d1

d2

E1

E2

i

ik

ki

k

9
7

6

係数行列の非対角成分

T2K-FVM-01

35

初期全体メッシュデータ（例）（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_VAL(ib)read (IUNIT,'(10i10)') FIX_NODE_tot
do i= 1, FIX_NODE_tot
read (IUNIT, '(i10, 3e16.6)')
FIX_NODE_ID(i), S, d,
FIX_NODE_VAL(i)
enddo

9
5
1

11

12
10

7

8

3

6
2

4

S

&
&

d
T=TB

T2K-FVM-01

36

1つの「要素」に複数のディリクレ
境界条件を与えることも実は可能
7 ディリクレ境界条件を与える要素数： 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
101.000000E+00
5.000000E-01
1.000000E+01
要素番号 境界面面積S
境界面と要素重心距離d 境界値TB
FIX_NODE_ID(ib)FIX_NODE_VAL(ib)

9
5
1

11

12
10

7

8

3

6
2

4

TB=10

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

10

TB=0

T2K-FVM-0137

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

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

38

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

TBｂ

i
TBa

q1

T2K-FVM-01

39

初期全体メッシュデータ（例）（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)read (IUNIT,'(10i10)') SURF_NODE_tot
do i= 1, SURF_NODE_tot
read (IUNIT, '(i10, 3e16.6)')
SURF_NODE_ID(i), S, QS
enddo

9
5
1
S
&

QS

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

11

12
10

7

8

3

6
2

4

T2K-FVM-01

40

初期全体メッシュデータ（例）（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

41

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

T2K-FVM-01

42

MicroAVSの使用• 手順は以下– http://nkl.cc.u-tokyo.ac.jp/10e/mavs/• MicroAVSについて– 本講習会では，10分バージョンを（株）サイバネットシステム様のご厚意で利用できます– 別途，評価版を一定期間利用することも可能です– 利用規定の遵守をお願いいたします• AVS/Expressに関しては本センターで利用可能です– 入力ファイル形式はUCDでMicroAVSと同じものが利用可能

T2K-FVM-01

43

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

キーワード

点

pt

線

line

三角形

tri

四角形

quad

四面体

tet

角錐

pyr

三角柱

prism

六面体

hex

二次要素
線2

line2

三角形2

tri2

四角形2

quad2

四面体2

tet2

角錐2

pyr2

三角柱2

prism2六面体2

hex2

T2K-FVM-01

44

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

T2K-FVM-01

45

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

・
・
・

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

・
・
・

(ステップ番号2) (コメント) (全節点数) (全要素数)
・・・

T2K-FVM-01

46

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

47

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
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.004 1.00…
10 1.00
11 1.00
12 1.00

新フォーマット
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.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
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
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

48• 有限体積法の概要• １CPU用プログラム eps_fvm– メッシュ生成，形状データ– 計算実行例– プログラムの内容

T2K-FVM-01

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

ロードモジュールの生成確認$> ./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

50

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

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

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

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

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

T2K-FVM-01

51

計算結果

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

T2K-FVM-01

52• 有限体積法の概要• １CPU用プログラム eps_fvm– メッシュ生成，形状データ– 計算実行例– プログラムの内容

T2K-FVM-01

53

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

メインプログラム

program eps_fvm
use eps_fvm_all
implicit REAL*8 (A-H,O-Z)
call eps_fvm_input_grid
call poi_gen
call eps_fvm_solve

call output_ucd
end program eps_fvm

T2K-FVM-01

54「eps_fvm」処理：メイン（1/2）
eps_fvm_all
program eps_fvm
use eps_fvm_all
implicit REAL*8 (A-H,O-Z)

call eps_fvm_input_grid
call poi_gen
call eps_fvm_solver
call output_ucd

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

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
use appl_cntl
end module eps_fvm_all

T2K-FVM-01

55「eps_fvm」処理：メイン（2/2）
その他
program eps_fvm
use eps_fvm_all
implicit REAL*8 (A-H,O-Z)

call eps_fvm_input_grid
call poi_gen
call eps_fvm_solver
call output_ucd
end program eps_fvm

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

T2K-FVM-01

56

program eps_fvm
use eps_fvm_all
implicit REAL*8 (A-H,O-Z)

call eps_fvm_input_grid
call poi_gen
call eps_fvm_solver
call output_ucd
end program eps_fvm

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

T2K-FVM-01

57

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

型

配列サイズ

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(:)

R

BODY_NODE_tot

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

内容T2K-FVM-01

58「eps_fvm」処理：メッシュ読み込み（1/5）
subroutine eps_fvm_input_grid
use hpcmw_eps_fvm_all
implicit REAL*8 (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出力
用に必要

IUNIT= 11
open (IUNIT,file= ‘fvm_entire_mesh.dat', status='unknown')!C!C-- NODE
read (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

59「eps_fvm」処理：メッシュ読み込み（1/5）!C!C-- NODE
read (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

&

&

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
111.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
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
1.500000E+00
5.000000E-01
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)

T2K-FVM-01

60「eps_fvm」処理：メッシュ読み込み（2/5）
コネクティビティ情報!C!C-- CONNECTION
read (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
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
41.000000E+00
10
8
11
12
10
11
12
12
E2

1.000000E+00
1.000000E+00
1.000000E+00
1.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
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
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
d2：E2重心～
境界面距離

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

T2K-FVM-0161

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

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

62「eps_fvm」処理：メッシュ読み込み（2/5）
コネクティビティ情報!C!C-- CONNECTION
read (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
C1 = 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

Sik
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

63「eps_fvm」処理：メッシュ読み込み（3/5）
ディリクレ境界条件!C!C-- DIRICHLET
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, '(i10, 3e16.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
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_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

64「eps_fvm」処理：メッシュ読み込み（4/5）
ノイマン境界条件!C!C-- NEUMANN
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

65「eps_fvm」処理：メッシュ読み込み（5/5）
体積発熱!C!C-- BODY FLUX
read (IUNIT,'(10i10)') BODY_NODE_tot
allocate (BODY_NODE_FLUX(NODE_tot))
do i= 1, BODY_NODE_tot
read (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

66

program eps_fvm
use eps_fvm_all
implicit REAL*8 (A-H,O-Z)

call eps_fvm_input_grid
call poi_gen
call eps_fvm_solver
call output_ucd
end program eps_fvm

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

T2K-FVM-01

67

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

型

サイズ

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

NPLU

I

-

D(:)

R

NODE_tot

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

PHI(:)

R

NODE_tot

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

BFORCE(:)

R

NODE_tot

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

index(:)

I

0:NODE_tot

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

item(:)I

NPLU

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

AMAT(:)R

NPLU

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

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

T2K-FVM-01

68

行列ベクトル積への適用（非零）非対角成分のみを格納，疎行列向け方法
Compressed Row Storage (CRS)
D(i) 対角成分（実数，i=1,N）index(i) 非対角成分数に関する一次元配列（通し番号）（整数，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)do k= index(i-1)+1, index(i)Y(i)= Y(i) + AMAT(k)*X(item(k))
enddo
enddo

T2K-FVM-01

69

行列ベクトル積：密行列⇒とても簡単
 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 

 

...       
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, N
Y(j)= Y(j) + A(i,j)*X(i)
enddo
enddo

T2K-FVM-01

70

Compressed Row Storage (CRS)
1

2

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
03.1
051.3
0

0

3.2

0

0

T2K-FVM-01

71

Compressed Row Storage (CRS)
1
2

1

2

1.1
①
4.3
①

2.4
②
3.6
②

4

3.1
①

4.1
②
9.5
②

6
7
8

6.4
②
9.5
②

4

5

6

7

8
NODE_tot= 8

3.2
⑤
2.5
④

5.7
③

3

5

3

3.7
⑥

9.1
⑧

1.5
3.1
⑤
⑦
9.8 2.5 2.7
④
⑤
⑥
10.4
11.5
4.3
③
⑤
⑦
6.5
12.4 9.5③⑥⑦
2.5
1.4 23.1 13.1③⑥⑦⑧
1.3 9.6
3.1
51.3
③
④
⑥
⑧

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

T2K-FVM-01

72

Compressed Row Storage (CRS)
1
2

3
4
5
6
7
8

1

2

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

2.4
②

3

4

5

6

7

8

3.2
⑤
2.5
④

4.1
②
9.5 10.4②③
6.5③
6.4 2.5②③
9.5 1.3 9.6
②
③
④

3.7
⑥

1.5
⑤
2.5
⑤

9.1
⑧

3.1
⑦
2.7
⑥

4.3
⑦
9.5
⑦
1.4
⑥
3.1
⑥

13.1
⑧

T2K-FVM-01

73

Compressed Row Storage (CRS)
非対角
成分数
1
2

3
4
5
6
7
8

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

2.4
②
4.3
①
1.5
⑤
4.1
②
3.1
①
6.5
③
6.4
②
9.5
②

3.2
⑤
2.5 3.7 9.1
④
⑥
⑧
3.1
⑦
2.5 2.7⑤⑥
9.5 10.4 4.3
②
③
⑦
9.5
⑦
2.5 1.4 13.1③⑥⑧
1.3 9.6 3.1③④⑥

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

4index(5)= 15

2index(6)= 17

4index(7)= 21

4index(8)= 25

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

NPLU= 25(=index(N))

T2K-FVM-01

74

Compressed Row Storage (CRS)
非対角
成分数
1
2

3
4
5
6
7
8

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

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

4index(5)= 15

2index(6)= 17

4index(7)= 21

4index(8)= 25

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

NPLU= 25(=index(N))

T2K-FVM-01

75

Compressed Row Storage (CRS)
1
2

3
4
5
6
7
8

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

例：
item( 7)= 5, AMAT( 7)= 1.5
item(19)= 3, AMAT(19)= 2.5

T2K-FVM-01

76

Compressed Row Storage (CRS)
1
2

3
4
5
6
7
8

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

D(i) 対角成分（実数，i=1,N）index(i) 非対角成分数に関する一次元配列（通し番号）（整数，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)do k= index(i-1)+1, index(i)Y(i)= Y(i) + AMAT(k)*X(item(k))
enddo
enddo

T2K-FVM-01

77

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

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

78

全体マトリクスの生成
要素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

79

全体マトリクスの生成
要素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

80

全体マトリクスの生成
要素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 ie
d ie
e
d

i

定数項：右辺へ移項

i

T2K-FVM-01

81

全体マトリクスの生成
要素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 ie
d ie
e
d

i

i








Sik
Sie 
Sik
Sie

Ti  
Tk  
Sid qid  Vi Q i 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「eps_fvm」処理：マトリクス生成（1/5）!C!C***!C*** POI_GEN
!C***
!C
!C
generate COEF. MATRIX for POISSON equations!C
subroutine POI_GEN
use eps_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= 082

T2K-FVM-01

83「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
index(in2)= ik2enddo
do i= 1, NODE_tot
index(i)= index(i-1) + index(i)enddo
NPLU= index(NODE_tot)
allocate (item(NPLU), AMAT(NPLU))（非零）非対角成分数算出（非零）非対角成分抽出
index(in)：
各要素の非対角成分数（通し番号）IWKX(in,1-6)：
各要素の非対角成分のコネクティビティID（CONN_NODE）
非対角要素数を増やした場合は
IWKXのサイズを増やせばよい（現在は隣接要素数が6までに限定されている）T2K-FVM-01

84「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
index(in2)= ik2enddo
do i= 1, NODE_tot
index(i)= index(i-1) + index(i)enddoindex(in)：
各要素の非対角成分数（通し番号）

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

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

T2K-FVM-01

85「eps_fvm」処理：マトリクス生成（3/5）
do i= 1, NODE_tot
do j= 1, index(i)-index(i-1)
k= index(i-1) + j

index(in)：
各要素の非対角成分数（通し番号）

ic = IWKX(i,j)
in1= CONN_NODE(2*ic-1)
in2= CONN_NODE(2*ic )!C===

if (in1.eq.i) thenitem(k)= in2
else
item(k)= in1
endif
enddo
enddo

ここまでがCRS形式
行列の準備

!C
!C +----------------------------+!C | INTERIOR NODEs + BODY FLUX |!C +----------------------------+!C===
do icel= 1, NODE_tot
BFORCE(icel)= BFORCE(icel) + BODY_NODE_FLUX(icel)
enddo

!C===

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)IWKX(in,1-6)：
各要素の非対角成分のコネクティビティID（CONN_NODE）
item(k)：非対角成分（列番号）T2K-FVM-01

86「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!C===








Sik
Sie 
Sik

Ti  Tk d ie 
 k d ik  d ki
 k d ik  d ki 
e


i k
i 
i k Sie

Sid qid  Vi Q i 
TiBe
d
ie
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

!C===

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)





i

T2K-FVM-01

87「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!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 







!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

Sie
T
d ie iBei

T2K-FVM-01

88「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

Sie
TiBe BFORCE（右辺）
d ie

i

FIX_NODE_COEF * FIX_NODE_VAL

TB
TiT2K-FVM-01

89「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








Sik
Sie 
Sik

Ti  

Tk 
d
d
d
d
d
ie 
 k ik  ki
 k ik  ki 
e


i k
i 
i k Sie

Sid qid  Vi Q i 
TiBe BFORCE（右辺）
d ie
d
e







SURF_NODE_FLUX





i

T2K-FVM-01

90

program eps_fvm
use eps_fvm_all
implicit REAL*8 (A-H,O-Z)

call eps_fvm_input_grid
call poi_gen
call eps_fvm_solver
call output_ucd
end program eps_fvm

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

T2K-FVM-01

91「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_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,*) 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

&
&

92

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

93

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

94

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

95

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

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

T2K-FVM-01

96

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

T2K-FVM-0197

共役勾配法のアルゴリズム
Compute r(0)= b-[A]x(0)for i= 1, 2, …
z(i-1)= r(i-1)
ri-1= r(i-1) z(i-1)
if i=1p(1)= z(0)elsebi-1= ri-1/ri-2p(i)= z(i-1) + bi-1 p(i-1)
endifq(i)= [A]p(i)ai = ri-1/p(i)q(i)x(i)= x(i-1) + aip(i)r(i)= r(i-1) - aiq(i)check convergence |r|end• 行列ベクトル積• ベクトル内積• ベクトル定数倍の加減x(i) ：ベクトル
ai ：スカラー

T2K-FVM-0198

共役勾配法のアルゴリズム
Compute r(0)= b-[A]x(0)for i= 1, 2, …
z(i-1)= r(i-1)
ri-1= r(i-1) z(i-1)
if i=1p(1)= z(0)elsebi-1= ri-1/ri-2p(i)= z(i-1) + bi-1 p(i-1)
endifq(i)= [A]p(i)ai = ri-1/p(i)q(i)x(i)= x(i-1) + aip(i)r(i)= r(i-1) - aiq(i)check convergence |r|end• 行列ベクトル積• ベクトル内積• ベクトル定数倍の加減x(i) ：ベクトル
ai ：スカラー

T2K-FVM-0199

共役勾配法のアルゴリズム
Compute r(0)= b-[A]x(0)for i= 1, 2, …
z(i-1)= r(i-1)
ri-1= r(i-1) z(i-1)
if i=1p(1)= z(0)elsebi-1= ri-1/ri-2p(i)= z(i-1) + bi-1 p(i-1)
endifq(i)= [A]p(i)ai = ri-1/p(i)q(i)x(i)= x(i-1) + aip(i)r(i)= r(i-1) - aiq(i)check convergence |r|end• 行列ベクトル積• ベクトル内積• ベクトル定数倍の加減x(i) ：ベクトル
ai ：スカラー

T2K-FVM-01100

共役勾配法のアルゴリズム
Compute r(0)= b-[A]x(0)for i= 1, 2, …
z(i-1)= r(i-1)
ri-1= r(i-1) z(i-1)
if i=1p(1)= z(0)elsebi-1= ri-1/ri-2p(i)= z(i-1) + b i-1 p(i-1)
endifq(i)= [A]p(i)ai = ri-1/p(i)q(i)x(i)= x(i-1) + aip(i)r(i)= r(i-1) - aiq(i)check convergence |r|end• 行列ベクトル積• ベクトル内積• ベクトル定数倍の加減（DAXPY）

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

T2K-FVM-01

101

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

T2K-FVM-01

102

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

T2K-FVM-01

103

前処理付共役勾配法
Preconditioned Conjugate Gradient Method （PCG）
Compute r(0)= b-[A]x(0)for i= 1, 2, …
solve [M]z(i-1)= r(i-1)
ri-1= r(i-1) z(i-1)
if i=1p(1)= z(0)elsebi-1= ri-1/ri-2p(i)= z(i-1) + bi-1 z(i-1)
endifq(i)= [A]p(i)ai = ri-1/p(i)q(i)x(i)= x(i-1) + aip(i)r(i)= r(i-1) - aiq(i)check convergence |r|end

実際にやるべき計算は：z  M 1r
「近似逆行列」の計算が必要：M   A , M   A1

1

究極の前処理：本当の逆行列M 1  A1 , M   A

対角スケーリング：簡単＝弱いM 1  D1 , M   D

T2K-FVM-01

104

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

105

対角スケーリング，点ヤコビ前処理• 前処理行列として，もとの行列の対角成分のみを取り出した行列を前処理行列 [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

106「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), dimension(N) :: X
integer, dimension(0:N) :: index
integer, dimension(NPLU):: item
real(kind=kreal), dimension(NPLU):: COEF&&

係数行列の配列名が変更
本体：
AMAT
ソルバー：COEF

real(kind=kreal) :: EPS

integer :: ITR, IER
integer :: P, Q, R, Z, DD
real(kind=kreal), dimension(:,:), allocatable, save :: W

T2K-FVM-01

107「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(i,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= 1, 2, …
solve [M]z(i-1)= r(i-1)
ri-1= r(i-1) z(i-1)
if i=1p(1)= z(0)elsebi-1= ri-1/ri-2p(i)= z(i-1) + bi-1 p(i-1)
endifq(i)= [A]p(i)ai = ri-1/p(i)q(i)x(i)= x(i-1) + aip(i)r(i)= r(i-1) - aiq(i)check convergence |r|end

T2K-FVM-01

108「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

!C===

IER = 0Tcomm= 0.d0

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

Compute r(0)= b-[A]x(0)for i= 1, 2, …
solve [M]z(i-1)= r(i-1)
ri-1= r(i-1) z(i-1)
if i=1p(1)= z(0)elsebi-1= ri-1/ri-2p(i)= z(i-1) + bi-1 p(i-1)
endifq(i)= [A]p(i)ai = ri-1/p(i)q(i)x(i)= x(i-1) + aip(i)r(i)= r(i-1) - aiq(i)check convergence |r|end

T2K-FVM-01

109「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 j= index(i-1)+1, index(i)W(i,R) = W(i,R) + COEF(j) * X(item(j))
enddo
enddo

!C===

BNRM2= 0.0D0
do i= 1, N
BNRM2= BNRM2 + B(i) **2
W(i,R)= B(i) - W(i,R)
enddo

Compute r(0)= b-[A]x(0)for i= 1, 2, …
solve [M]z(i-1)= r(i-1)
ri-1= r(i-1) z(i-1)
if i=1p(1)= z(0)elsebi-1= ri-1/ri-2p(i)= z(i-1) + bi-1 p(i-1)
endifq(i)= [A]p(i)ai = ri-1/p(i)q(i)x(i)= x(i-1) + aip(i)r(i)= r(i-1) - aiq(i)check convergence |r|end

T2K-FVM-01

110「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 j= index(i-1)+1, index(i)W(i,R) = W(i,R) + COEF(j) * X(item(j))
enddo
enddo

!C===

BNRM2= 0.0D0
do i= 1, N
BNRM2= BNRM2 + B(i) **2
W(i,R)= B(i) - W(i,R)
enddo

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

Compute r(0)= b-[A]x(0)for i= 1, 2, …
solve [M]z(i-1)= r(i-1)
ri-1= r(i-1) z(i-1)
if i=1p(1)= z(0)elsebi-1= ri-1/ri-2p(i)= z(i-1) + bi-1 p(i-1)
endifq(i)= [A]p(i)ai = ri-1/p(i)q(i)x(i)= x(i-1) + aip(i)r(i)= r(i-1) - aiq(i)check convergence |r|end

T2K-FVM-01

111「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(i,R)*W(i,Z)
enddo
!C===

ITERATION

Compute r(0)= b-[A]x(0)for i= 1, 2, …
solve [M]z(i-1)= r(i-1)
ri-1= r(i-1) z(i-1)
if i=1p(1)= z(0)elsebi-1= ri-1/ri-2p(i)= z(i-1) + bi-1 p(i-1)
endifq(i)= [A]p(i)ai = ri-1/p(i)q(i)x(i)= x(i-1) + aip(i)r(i)= r(i-1) - aiq(i)check convergence |r|end

T2K-FVM-01

112「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}= [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) + COEF(j) * W(item(j),P)
enddo
enddo
!C===

Compute r(0)= b-[A]x(0)for i= 1, 2, …
solve [M]z(i-1)= r(i-1)
ri-1= r(i-1) z(i-1)
if i=1p(1)= z(0)else
b i-1= ri-1/ri-2p(i)= z(i-1) + b i-1 p(i-1)
endifq(i)= [A]p(i)ai = ri-1/p(i)q(i)x(i)= x(i-1) + aip(i)r(i)= r(i-1) - aiq(i)check convergence |r|end

T2K-FVM-01

113「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} - ALPHA*{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= 1, 2, …
solve [M]z(i-1)= r(i-1)
ri-1= r(i-1) z(i-1)
if i=1p(1)= z(0)elsebi-1= ri-1/ri-2p(i)= z(i-1) + bi-1 p(i-1)
endifq(i)= [A]p(i)ai = ri-1/p(i)q(i)x(i)= x(i-1) + aip(i)r(i)= r(i-1) - aiq(i)check convergence |r|end

T2K-FVM-01

114「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)
if ( RESID.le.EPS) goto 900
RHO1 = RHO

enddo
IER = 1900 continue
ITR= L
EPS= RESID

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

RESID= |r|/|b|returnend subroutine hpcmw_eps_fvm_solver_CG

Compute r(0)= b-[A]x(0)for i= 1, 2, …
solve [M]z(i-1)= r(i-1)
ri-1= r(i-1) z(i-1)
if i=1p(1)= z(0)elsebi-1= ri-1/ri-2p(i)= z(i-1) + bi-1 p(i-1)
endifq(i)= [A]p(i)ai = ri-1/p(i)q(i)x(i)= x(i-1) + aip(i)r(i)= r(i-1) - aiq(i)check convergence |r|end

T2K-FVM-01

115

いろいろなメッシュでやってみよう