並列有限要素法による三次元定常熱伝導解析プログラム（2/2）Fortran編
中島 研吾
東京大学情報基盤センター

pFEM3D-2

2

対象とする問題：三次元定常熱伝導∂  ∂T  ∂  ∂T  ∂  ∂T  ɺ +  λλ +  λ + Q ( x, y , z ) = 0∂x  ∂x  ∂y  ∂y  ∂z  ∂z 

• 定常熱伝導＋発熱• 一様な熱伝導率 λ
• 直方体

Z

T=0@Z=zmax– 一辺長さ1の立方体（六面体）要素– 各方向にNX・NY・NZ個
NZ• 境界条件– T=0@Z=zmax

Y
NX
NY
X

• 体積当たり発熱量は位置（メッシュの中心の座標 xc,yc）に依存– Qɺ ( x, y, z ) = QVOL xC + yC

pFEM3D-2

3

有限要素法の処理• 支配方程式• ガラーキン法：弱形式• 要素単位の積分– 要素マトリクス生成• 全体マトリクス生成• 境界条件適用• 連立一次方程式

pFEM3D-2

4

並列有限要素法の処理：プログラム• 初期化– 制御変数読み込み– 座標読み込み⇒要素生成（N:節点数，NE：要素数）– 配列初期化（全体マトリクス，要素マトリクス）– 要素⇒全体マトリクスマッピング（Index，Item）• マトリクス生成– 要素単位の処理（do icel= 1, NE）• 要素マトリクス計算• 全体マトリクスへの重ね合わせ– 境界条件の処理• 連立一次方程式– 共役勾配法（CG）

pFEM3D-2

5

並列有限要素法の手順（並列計算実行）
pfem3d/run/
INPUT.DAT

pfem3d/mesh/<HEADER>.*

pfem3d/run/sol

局所分散メッシュファイル

pfem3d/run/test.inp
ParaView 出力：名称固定

pFEM3D-2

6制御ファイル：INPUT.DAT../mesh/aaa
2000
1.0 1.0
1.0e-08• HEADER：• ITER：• COND：• QVOL：• RESID：

HEADER
ITER
COND, QVOLRESID
局所分散ファイルヘッダ名<HEADER>.my_rank
反復回数上限熱伝導率
体積当たり発熱量係数
反復法の収束判定値∂  ∂T  ∂  ∂T  ∂  ∂T  ɺ +  λλ +  λ + Q ( x, y , z ) = 0∂x  ∂x  ∂y  ∂y  ∂z  ∂z Qɺ ( x, y, z ) = QVOL xC + yC

pFEM3D-2

7

pFEM/pfem3d/run/t01x16x2.sh
16-cores/socket, 32-cores/node#PJM -N "Flatx16x2“#PJM -L rscgrp=tutorial#PJM –L node=1#PJM –mpi proc=32#PJM –L elapse=00:15:00#PJM –g gt00
#PJM -j#PJM -e err#PJM -o t01x16x2_0001.lst

Job Name
Name of “QUEUE”
node#
Total # of MPI Processes
TimeGroup Name (Wallet)
Standard Error
Standard Output

mpiexec.hydra -n ${PJM_MPI_PROC} ./sol
mpiexec.hydra -n ${PJM_MPI_PROC} numactl -l ./sol#PJM -L node=1/--mpi proc= 32#PJM -L node=2/--mpi proc= 64#PJM -L node=4/--mpi proc=128#PJM -L node=8/--mpi proc=256

Memory
96 GB

DDR4
DDR4
DDR4
DDR4
DDR4
DDR4

2933 MHz×6ch
140.8 GB/sec

Intel® Xeon®
Platinum 8280
(Cascade Lake, CLX)
2.7GHz, 28-Cores
2.419 TFLOPS
Soc. #0: 0th-27th cores

UPI
UPI
UPI

1-nodes, 32-processes
2-nodes, 64-processes
4-nodes, 128-processes
8-nodes, 256-processes
Intel® Xeon®
Platinum 8280
(Cascade Lake, CLX)
2.7GHz, 28-Cores
2.419 TFLOPS
Soc. #1: 28th-55th cores

Ultra Path Interconnect
10.4 GT/sec ×3= 124.8 GB/sec

DDR4
DDR4
DDR4
DDR4
DDR4
DDR4

Memory96 GB2933 MHz×6ch
140.8 GB/sec

pFEM3D-2

8

pFEM/pfem3d/run/p20.sh, k20.sh
20-cores/socket, 40-cores/node#PJM -N "Flatx20x2“#PJM -L rscgrp=tutorial#PJM -L node=1#PJM --mpi proc=40#PJM -L elapse=00:15:00#PJM -g gt00
#PJM -j#PJM -e err#PJM -o p01x20x2_0001.lst

Job Name
Name of “QUEUE”
node#
Total # of MPI Processes
TimeGroup Name (Wallet)
Standard Error
Standard Output

mpiexec.hydra -n ${PJM_MPI_PROC} ./sol
mpiexec.hydra -n ${PJM_MPI_PROC} numactl -l ./sol#PJM -L node=1/--mpi proc= 40#PJM -L node=2/--mpi proc= 80#PJM -L node=4/--mpi proc=160#PJM -L node=8/--mpi proc=320

Memory
96 GB

DDR4
DDR4
DDR4
DDR4
DDR4
DDR4

2933 MHz×6ch
140.8 GB/sec

Intel® Xeon®
Platinum 8280
(Cascade Lake, CLX)
2.7GHz, 28-Cores
2.419 TFLOPS
Soc. #0: 0th-27th cores

UPI
UPI
UPI

1-nodes, 40-processes
2-nodes, 80-processes
4-nodes, 160-processes
8-nodes, 320-processes
Intel® Xeon®
Platinum 8280
(Cascade Lake, CLX)
2.7GHz, 28-Cores
2.419 TFLOPS
Soc. #1: 28th-55th cores

Ultra Path Interconnect
10.4 GT/sec ×3= 124.8 GB/sec

DDR4
DDR4
DDR4
DDR4
DDR4
DDR4

Memory96 GB2933 MHz×6ch
140.8 GB/sec

pFEM3D-2

9

pFEM/pfem3d/run/p24.sh, k24.sh
24-cores/socket, 48-cores/node#PJM -N "Flatx24x2“#PJM -L rscgrp=tutorial#PJM -L node=1#PJM --mpi proc=48#PJM -L elapse=00:15:00#PJM -g gt00
#PJM -j#PJM -e err#PJM -o p01x24x2_0001.lst

Job Name
Name of “QUEUE”
node#
Total # of MPI Processes
TimeGroup Name (Wallet)
Standard Error
Standard Output

mpiexec.hydra -n ${PJM_MPI_PROC} ./sol
mpiexec.hydra -n ${PJM_MPI_PROC} numactl -l ./sol#PJM -L node=1/--mpi proc= 48#PJM -L node=2/--mpi proc= 96#PJM -L node=4/--mpi proc=192#PJM -L node=8/--mpi proc=384

Memory
96 GB

DDR4
DDR4
DDR4
DDR4
DDR4
DDR4

2933 MHz×6ch
140.8 GB/sec

Intel® Xeon®
Platinum 8280
(Cascade Lake, CLX)
2.7GHz, 28-Cores
2.419 TFLOPS
Soc. #0: 0th-27th cores

UPI
UPI
UPI

1-nodes, 48-processes
2-nodes, 96-processes
4-nodes, 192-processes
8-nodes, 384-processes
Intel® Xeon®
Platinum 8280
(Cascade Lake, CLX)
2.7GHz, 28-Cores
2.419 TFLOPS
Soc. #1: 28th-55th cores

Ultra Path Interconnect
10.4 GT/sec ×3= 124.8 GB/sec

DDR4
DDR4
DDR4
DDR4
DDR4
DDR4

Memory96 GB2933 MHz×6ch
140.8 GB/sec

pFEM3D-2

10

pFEM/pfem3d/run/p28.sh, k28.sh
28-cores/socket, 56-cores/node#PJM -N "Flatx24x2“#PJM -L rscgrp=tutorial#PJM -L node=1#PJM --mpi proc=56#PJM -L elapse=00:15:00#PJM -g gt00
#PJM -j#PJM -e err#PJM -o p01x28x2_0001.lst

Job Name
Name of “QUEUE”
node#
Total # of MPI Processes
TimeGroup Name (Wallet)
Standard Error
Standard Output

mpiexec.hydra -n ${PJM_MPI_PROC} ./sol
mpiexec.hydra -n ${PJM_MPI_PROC} numactl -l ./sol#PJM -L node=1/--mpi proc= 56#PJM -L node=2/--mpi proc=112#PJM -L node=4/--mpi proc=224#PJM -L node=8/--mpi proc=448

Memory
96 GB

DDR4
DDR4
DDR4
DDR4
DDR4
DDR4

2933 MHz×6ch
140.8 GB/sec

Intel® Xeon®
Platinum 8280
(Cascade Lake, CLX)
2.7GHz, 28-Cores
2.419 TFLOPS
Soc. #0: 0th-27th cores

UPI
UPI
UPI

1-nodes, 56-processes
2-nodes, 112-processes
4-nodes, 224-processes
8-nodes, 448-processes
Intel® Xeon®
Platinum 8280
(Cascade Lake, CLX)
2.7GHz, 28-Cores
2.419 TFLOPS
Soc. #1: 28th-55th cores

Ultra Path Interconnect
10.4 GT/sec ×3= 124.8 GB/sec

DDR4
DDR4
DDR4
DDR4
DDR4
DDR4

Memory96 GB2933 MHz×6ch
140.8 GB/sec

pFEM3D-2

11

test1

input_cntl

メインプログラム制御データ入力

input_griddefine_file_name

メッシュファイル入力

局所ファイル名

mat_con0mSORT

行列コネクティビティ生成

ソート

mat_con1

find_node

行列コネクティビティ生成

節点探索

mat_ass_main

jacobi係数行列生成ヤコビアン計算

mat_ass_bc
境界条件処理三次元並列熱伝
導解析コードheat3Dp

solve11cg

線形ソルバー制御

CG法計算

output_ucd
可視化処理

pFEM3D-2

12

全体処理
program heat3Dp
use solver11
use pfem_util
implicit REAL*8(A-H,O-Z)
call PFEM_INIT
call INPUT_CNTL
call INPUT_GRID
call MAT_CON0
call MAT_CON1
call MAT_ASS_MAIN
call MAT_ASS_BC
call SOLVE11
call OUTPUT_UCD
call PFEM_FINALIZE
end program heat3Dp

pFEM3D-2

13

Global変数表：pfem_util.f（1/4）
変数名

種別

サイズ

内 容

I/O
I

メッシュファイル名I

I

節点数（N：内点，NP：内点＋外点）

ICELTOTII

要素数NODGRPtot

I

I

節点グループ数

XYZ

R(NP,3)I

節点座標

ICELNODI(ICELTOT,8)I

要素コネクティビティ

NODGRP_INDEX

I(0:NODGRPtot)

I

各節点グループに含まれる節点数（累積）NODGRP_ITEMI(NODGRP_INDEX(NODG
RPTOT)

I

節点グループに含まれる節点

NODGRP_NAME

C80(NODGRPTOT)I

節点グループ名

NLU

IO

各節点非対角成分数

NPLU

I

O

非対角成分総数

D

R

(NP)

O

全体行列：対角ブロック

B,X

R

(NP)

O

右辺ベクトル，未知数ベクトル

fname

C

N,NP

(80)

pFEM3D-2

14

Global変数表：pfem_util.f（2/4）
変数名

種別

サイズ

I/O

内 容

AMAT

R

(NPLU)

O

全体行列：非零非対角成分

index

I

(0:NP)

O

全体行列：非零非対角成分数

item

I

(NPLU)

O

全体行列：非零非対角成分（列番号）

INLU

I

(NP)

O

各節点の非零非対角成分数

IALU

I

(NP,NLU)

O

各節点の非零非対角成分数（列番号）

IWKX

I

(NP,2)

O

ワーク用配列ITER, ITERactualI

I

反復回数の上限，実際の反復回数

RESID

R

I

打ち切り誤差（1.e-8に設定）

pFEM3D-2

15

Global変数表：pfem_util.f（3/4）
変数名

種別

O8th

R

PNQ, PNE, PNT

R

POS, WEI

サイズ

内 容

I/O
I

=0.125(2,2,8)O

各ガウス積分点における

R

(2,2)

O

各ガウス積分点の座標，重み係数

NCOL1, NCOL2

I

(100)

O

ソート用ワーク配列

SHAPE

R

(2,2,2,8)O

各ガウス積分点における形状関数 Ni (i=1~8)

PNX, PNY, PNZ

R(2,2,2,8)O

各ガウス積分点における

DETJ

R

(2,2,2)O

各ガウス積分点におけるヤコビアン行列式

COND, QVOL

R

I

熱伝導率，体積当たり発熱量係数∂N i ∂N i ∂N i
(i = 1 ~ 8),,∂ξ ∂η ∂ζ∂N i ∂N i ∂N i
(i = 1 ~ 8),,∂x ∂y ∂zpFEM3D-2

16

Global変数表：pfem_util.f（4/4）
変数名

種別

サイズ

内 容

I/O

PETOT

I

O

領域数（MPIプロセス数）

my_rank

I

O

MPIプロセス番号

errno

I

O

エラーフラグ

NEIBPETOT

I

I

隣接領域数

NEIBPEI(NEIBPETOT)I

隣接領域番号

IMPORT_INDEXEXPORT_INEDXI(0:NEIBPETOT)

I

送信，受信テーブルのサイズ（一次元圧縮配列）

IMPORT_ITEM

I(Npimport)

I

受信テーブル（外点）（NPimport=IMPORT_INDEX(NEIBPETOT))）EXPORT_ITEM

I(Npexport)

I

送信テーブル（境界点）（NPexport=EXPORT_INDEX(NEIBPETOT))）ICELTOT_INT

I

I

領域所属要素数

intELEM_list

II

領域所属要素のリスト:(ICELTOT_INT)

可視化に使用pFEM3D-2

17

開始，終了：MPI_Init/Finalize
subroutine PFEM_INIT
use pfem_util
implicit REAL*8 (A-H,O-Z)
call MPI_INIT(ierr)
call MPI_COMM_SIZE (MPI_COMM_WORLD, PETOT, ierr )call MPI_COMM_RANK (MPI_COMM_WORLD, my_rank, ierr )pfemRarray= 0.d0
pfemIarray= 0
returnend
subroutine PFEM_FINALIZE
use pfem_util
implicit REAL*8 (A-H,O-Z)
call MPI_FINALIZE (errno)
if (my_rank.eq.0) stop ' * normal termination'
returnend

pFEM3D-2

18制御ファイル入力：INPUT_CNTL
subroutine INPUT_CNTL
use pfem_util
implicit REAL*8 (A-H,O-Z)
if (my_rank.eq.0) thenopen (11,file= 'INPUT.DAT', status='unknown')
read (11,'(a80)') HEADER
read (11,*) ITER
read (11,*) COND, QVOL
read (11,*) RESID
close (11)endifcall MPI_BCAST (HEADER, 80, MPI_CHARACTER, 0, MPI_COMM_WORLD,ierr)call MPI_BCAST (ITER , 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)call MPI_BCAST (COND , 1, MPI_DOUBLE_PRECISION, 0,&
MPI_COMM_WORLD, ierr)call MPI_BCAST (QVOL , 1, MPI_DOUBLE_PRECISION, 0,&
MPI_COMM_WORLD, ierr)call MPI_BCAST (RESID , 1, MPI_DOUBLE_PRECISION, 0,&
MPI_COMM_WORLD, ierr)
pfemRarray(1)= RESID
pfemIarray(1)= ITER
returnend

pFEM3D-2

19

メッシュ入力：INPUT_GRID（1/3）
subroutine INPUT_GRID
use pfem_util
implicit REAL*8 (A-H,O-Z)
call define_file_name (HEADER, fname, my_rank)
open (11, file= fname, status= 'unknown', form= 'formatted')!C!C-- NEIB-PE
read (11,'(10i10)') kkk
read (11,'(10i10)') NEIBPETOT
allocate (NEIBPE(NEIBPETOT))
read (11,'(10i10)') (NEIBPE(i), i= 1, NEIBPETOT)
do i= 1, NEIBPETOT
if (NEIBPE(i).gt.PETOT-1) thencall ERROR_EXIT (202, my_rank)
endif
enddo

pFEM3D-2

20

分散メッシュファイル名：
DEFINE_FILE_NAME
HEADER+ランク番号，実際は106個まで可能
subroutine DEFINE_FILE_NAME (HEADERo, filename, my_rank)
character (len=80) ::character (len=80) ::character (len= 1) ::character (len= 2) ::character (len= 3) ::integer:: LENGTH, ID

HEADERo, filename
HEADER
SUBindex1
SUBindex2
SUBindex3

HEADER= adjustL (HEADERo)LENGTH= len_trim(HEADER)
if (my_rank.le.9) thenID= 1write(SUBindex1 ,'(i1.1)') my_rank
else if (my_rank.le.99) thenID= 2write(SUBindex2 ,'(i2.2)') my_rank
else if (my_rank.le.999) thenID= 3write(SUBindex3 ,'(i3.3)') my_rank
endif
if (ID.eq.1) filename= HEADER(1:LENGTH)//'.'//SUBindex1
if (ID.eq.2) filename= HEADER(1:LENGTH)//'.'//SUBindex2
if (ID.eq.3) filename= HEADER(1:LENGTH)//'.'//SUBindex3
end subroutine define_file_namepFEM3D-2

21

allocate, deallocate関数：C言語#include <stdio.h>#include <stdlib.h>
void* allocate_vector(int size,int m){void *a;if ( ( a=(void * )malloc( m * size ) ) == NULL ) {fprintf(stdout,"Error:Memory does not enough! in vector ¥n");exit(1);}
return a;}
allocateをFORTRAN並みに
void deallocate_vector(void *a){free( a );}簡単にやるための関数

void** allocate_matrix(int size,int m,int n){void **aa;
int i;
if ( ( aa=(void ** )malloc( m * sizeof(void*) ) ) == NULL ) {fprintf(stdout,"Error:Memory does not enough! aa in matrix ¥n");exit(1);}
if ( ( aa[0]=(void * )malloc( m * n * size ) ) == NULL ) {fprintf(stdout,"Error:Memory does not enough! in matrix ¥n");exit(1);}for(i=1;i<m;i++) aa[i]=(char*)aa[i-1]+size*n;
return aa;}
void deallocate_matrix(void **aa){
free( aa );}pFEM3D-2

22

メッシュ入力：INPUT_GRID（2/3）!C!C-- NODE
read (11,'(10i10)') NP, N
allocate (XYZ(NP,3), NODE_ID(NP,2))
XYZ= 0.d0
do i= 1, NP
read (11,*) NODE_ID(i,1), NODE_ID(i,2), (XYZ(i,kk),kk=1,3)
enddo!C!C-- ELEMENT
read (11,*)ICELTOT, ICELTOT_INT

allocate (ICELNOD(ICELTOT,8), intELEM_list(ICELTOT))
allocate (ELEM_ID(ICELTOT,2))
read (11,'(10i10)') (NTYPE, i= 1, ICELTOT)
do icel= 1, ICELTOT
read (11,'(i10,2i5,8i10)') (ELEM_ID(icel,jj),jj=1,2),&
IMAT, (ICELNOD(icel,k), k= 1, 8)
enddo
read (11,'(10i10)') (intELEM_list(ic0), ic0= 1, ICELTOT_INT)&pFEM3D-2

23

メッシュ入力：INPUT_GRID（3/3）!C-- COMMUNICATION tableallocate (IMPORT_INDEX(0:NEIBPETOT))
allocate (EXPORT_INDEX(0:NEIBPETOT))IMPORT_INDEX= 0
EXPORT_INDEX= 0
if (PETOT.ne.1) thenread (11,'(10i10)') (IMPORT_INDEX(i), i= 1, NEIBPETOT)
nn= IMPORT_INDEX(NEIBPETOT)
allocate (IMPORT_ITEM(nn))
do i= 1, nn
read (11,*) IMPORT_ITEM(i)enddoread (11,'(10i10)') (EXPORT_INDEX(i), i= 1, NEIBPETOT)
nn= EXPORT_INDEX(NEIBPETOT)
allocate (EXPORT_ITEM(nn))
do i= 1, nn
read (11,*) EXPORT_ITEM(i)
enddo
endif
!C-- NODE grp. info.
read (11,'(10i10)') NODGRPtot
allocate (NODGRP_INDEX(0:NODGRPtot),NODGRP_NAME(NODGRPtot))
NODGRP_INDEX= 0
read (11,'(10i10)') (NODGRP_INDEX(i), i= 1, NODGRPtot)
nn= NODGRP_INDEX(NODGRPtot)
allocate (NODGRP_ITEM(nn))
do k= 1, NODGRPtot
iS= NODGRP_INDEX(k-1) + 1
iE= NODGRP_INDEX(k )read (11,'(a80)') NODGRP_NAME(k)nn= iE - iS + 1
if (nn.ne.0) thenread (11,'(10i10)') (NODGRP_ITEM(kk),kk=iS, iE)
endif
enddo

pFEM3D-2

24

並列有限要素法の処理：プログラム• 初期化– 制御変数読み込み– 座標読み込み⇒要素生成（N:節点数，NE：要素数）– 配列初期化（全体マトリクス，要素マトリクス）– 要素⇒全体マトリクスマッピング（Index，Item）• マトリクス生成– 要素単位の処理（do icel= 1, NE）• 要素マトリクス計算• 全体マトリクスへの重ね合わせ– 境界条件の処理• 連立一次方程式– 共役勾配法（CG）

pFEM3D-2

25

test1

input_cntl

メインプログラム制御データ入力

1プロセッサの場合（heat3D）とほとんど変わらない

input_griddefine_file_name

メッシュファイル入力

局所ファイル名

mat_con0mSORT

行列コネクティビティ生成

ソート

mat_con1

find_node

行列コネクティビティ生成

節点探索

mat_ass_main

jacobi係数行列生成ヤコビアン計算

mat_ass_bc
境界条件処理三次元並列熱伝
導解析コードheat3Dp

solve11cg

線形ソルバー制御

CG法計算

output_ucd
可視化処理

pFEM3D-2

26

マトリクス生成まで• 一次元のときは，index，itemに関連した情報を簡単に作ることができた– 非ゼロ非対角成分の数は2– 番号が自分に対して：+1と-1• 三次元の場合はもっと複雑– 非ゼロ非対角ブロックの数は7～26（現在の形状）– 実際はもっと複雑– 前以て，非ゼロ非対角ブロックの数はわからない

pFEM3D-2

27

マトリクス生成まで• 一次元のときは，index，itemに関連した情報を簡単に作ることができた– 非ゼロ非対角成分の数は2– 番号が自分に対して：+1と-1• 三次元の場合はもっと複雑– 非ゼロ非対角ブロックの数は7～26（現在の形状）– 実際はもっと複雑– 前以て，非ゼロ非対角ブロックの数はわからない• INLU(N)，IALU(N,NLU) を使って非ゼロ非対角成分数を予備的に勘定する

pFEM3D-2

28

全体処理
program heat3Dp
use solver11
use pfem_util
implicit REAL*8(A-H,O-Z)
call PFEM_INIT
call INPUT_CNTL
call INPUT_GRID
call MAT_CON0
call MAT_CON1
call MAT_ASS_MAIN
call MAT_ASS_BC
call SOLVE11
call OUTPUT_UCD
call PFEM_FINALIZE
end program heat3Dp

MAT_CON0: INLU, IALU生成
MAT_CON1: index, item生成

pFEM3D-2

29

MAT_CON0：全体構成(− 1,+1,+1) 8

do icel= 1, ICELTOT
8節点相互の関係から，
INL，INU，IAL，IAUを生成（FIND_NODE）(− 1,−1,+1) 5

7 (+ 1,+1,+1)
6 (+ 1,−1,+1)

enddo(− 1,+1,−1)
4
1

13

14

7
9

8
10

4
5

11

6

12

6
7

2
2

16

9

5

1
1

15

8

3
3

4

(ξ ,η , ζ ) = (− 1,−1,−1)

2

3 (+ 1,+1,−1)(+ 1,−1,−1)pFEM3D-2

30

行列コネクティビティ生成：
MAT_CON0（1/4）!C!C***!C*** MAT_CON0!C***!C
subroutine MAT_CON0
use pfem_util
implicit REAL*8 (A-H,O-Z)
NLU= 26
allocate (INLU(NP), IALU(NP,NLU))INLU= 0IALU= 0

NLU：
各節点における
非ゼロ非対角成分の最大数（接続する節点数）
今の問題の場合はわかっているので，
このようにできる不明の場合の実装：⇒レポート課題

pFEM3D-2

31

行列コネクティビティ生成：
MAT_CON0（1/4）!C!C***!C*** MAT_CON0!C***!C
subroutine MAT_CON0
use pfem_util
implicit REAL*8 (A-H,O-Z)
NLU= 26
allocate (INLU(NP), IALU(NP,NLU))INLU= 0IALU= 0

変数名

サイズ

内

容

INLU

(NP)

各節点の非零非対角成
分数

IALU

(NP,NLU)

各節点の非零非対角成分（列番号）

pFEM3D-2

32

行列コネクティビティ生成：
MAT_CON0（2/4）
do icel= 1, ICELTOT
in1= ICELNOD(icel,1)
in2= ICELNOD(icel,2)
in3= ICELNOD(icel,3)
in4= ICELNOD(icel,4)
in5= ICELNOD(icel,5)
in6= ICELNOD(icel,6)
in7= ICELNOD(icel,7)
in8= ICELNOD(icel,8)call FIND_TS_NODE (in1,in2)
call FIND_TS_NODE (in1,in3)
call FIND_TS_NODE (in1,in4)
call FIND_TS_NODE (in1,in5)
call FIND_TS_NODE (in1,in6)
call FIND_TS_NODE (in1,in7)
call FIND_TS_NODE (in1,in8)
call FIND_TS_NODE (in2,in1)
call FIND_TS_NODE (in2,in3)
call FIND_TS_NODE (in2,in4)
call FIND_TS_NODE (in2,in5)
call FIND_TS_NODE (in2,in6)
call FIND_TS_NODE (in2,in7)
call FIND_TS_NODE (in2,in8)
call FIND_TS_NODE (in3,in1)
call FIND_TS_NODE (in3,in2)
call FIND_TS_NODE (in3,in4)
call FIND_TS_NODE (in3,in5)
call FIND_TS_NODE (in3,in6)
call FIND_TS_NODE (in3,in7)
call FIND_TS_NODE (in3,in8)(− 1,+1,+1) 8(− 1,−1,+1) 5

7 (+ 1,+1,+1)
6 (+ 1,−1,+1)(− 1,+1,−1)
4
1(ξ ,η , ζ ) = (− 1,−1,−1)

2

3 (+ 1,+1,−1)(+ 1,−1,−1)

FEM3D

33

節点探索：FIND_TS_NODEINLU,IAU探索：一次元ではこの部分は手動!C!C***!C*** FIND_TS_NODE!C***!Csubroutine FIND_TS_NODE (ip1,ip2)
do kk= 1, INLU(ip1)
if (ip2.eq.IALU(ip1,kk)) return
enddo
icou= INLU(ip1) + 1
IALU(ip1,icou)= ip2
INLU(ip1)= icou
returnend subroutine FIND_TS_NODE

変数名

サイズ

内

容

INLU

(N)

各節点の非零非対角成
分数

IALU

(N,NLU)

各節点の非零非対角成分（列番号）

FEM3D34

節点探索：FIND_TS_NODE
一次元ではこの部分は手動，要素#2!C!C***!C*** FIND_TS_NODE!C***!Csubroutine FIND_TS_NODE (ip1,ip2)
do kk= 1, INLU(ip1)
if (ip2.eq.IALU(ip1,kk)) return
enddo
icou= INLU(ip1) + 1
IALU(ip1,icou)= ip2
INLU(ip1)= icou
returnend subroutine FIND_TS_NODE

ip2がIALU(ip1,kk)に含まれているかチェック

13

14

7
9

8
10

4
5

ip1: No.6 node
ip2: No.7 node

11

6

12

6
7

2
2

16

9

5

1
1

15

8

3
3

4

FEM3D

35

節点探索：FIND_TS_NODE
一次元ではこの部分は手動，要素#2!C!C***!C*** FIND_TS_NODE!C***!Csubroutine FIND_TS_NODE (ip1,ip2)
do kk= 1, INLU(ip1)
if (ip2.eq.IALU(ip1,kk)) return
enddo
icou= INLU(ip1) + 1
IALU(ip1,icou)= ip2
INLU(ip1)= icou
returnend subroutine FIND_TS_NODEIALUに含まれていない
場合は，INLUに1を加えてIALUに格納
13

14

7
9

8
10

4
5

ip1: No.6 node
ip2: No.7 node

11

6

12

6
7

2
2

16

9

5

1
1

15

8

3
3

4

FEM3D36

節点探索：FIND_TS_NODE
一次元ではこの部分は手動，要素#5!C!C***!C*** FIND_TS_NODE!C***!Csubroutine FIND_TS_NODE (ip1,ip2)
do kk= 1, INLU(ip1)
if (ip2.eq.IALU(ip1,kk)) return
enddo
icou= INLU(ip1) + 1
IALU(ip1,icou)= ip2
INLU(ip1)= icou
returnend subroutine FIND_TS_NODE

既にIALUに含まれている場合は，次のペアへ

13

14

7
9

8
10

4
5

ip1: No.6 node
ip2: No.7 node

11

6

12

6
7

2
2

16

9

5

1
1

15

8

3
3

4

FEM3D37

行列コネクティビティ生成：
MAT_CON0（3/4）
call FIND_TS_NODE (in4,in1)
call FIND_TS_NODE (in4,in2)
call FIND_TS_NODE (in4,in3)
call FIND_TS_NODE (in4,in5)
call FIND_TS_NODE (in4,in6)
call FIND_TS_NODE (in4,in7)
call FIND_TS_NODE (in4,in8)
call FIND_TS_NODE (in5,in1)
call FIND_TS_NODE (in5,in2)
call FIND_TS_NODE (in5,in3)
call FIND_TS_NODE (in5,in4)
call FIND_TS_NODE (in5,in6)
call FIND_TS_NODE (in5,in7)
call FIND_TS_NODE (in5,in8)
call FIND_TS_NODE (in6,in1)
call FIND_TS_NODE (in6,in2)
call FIND_TS_NODE (in6,in3)
call FIND_TS_NODE (in6,in4)
call FIND_TS_NODE (in6,in5)
call FIND_TS_NODE (in6,in7)
call FIND_TS_NODE (in6,in8)
call FIND_TS_NODE (in7,in1)
call FIND_TS_NODE (in7,in2)
call FIND_TS_NODE (in7,in3)
call FIND_TS_NODE (in7,in4)
call FIND_TS_NODE (in7,in5)
call FIND_TS_NODE (in7,in6)
call FIND_TS_NODE (in7,in8)(− 1,+1,+1) 8(− 1,−1,+1) 5

7 (+ 1,+1,+1)
6 (+ 1,−1,+1)(− 1,+1,−1)
4
1(ξ ,η , ζ ) = (− 1,−1,−1)

2

3 (+ 1,+1,−1)(+ 1,−1,−1)

FEM3D38

行列コネクティビティ生成：
MAT_CON0（4/4）
call FIND_TS_NODE (in8,in1)
call FIND_TS_NODE (in8,in2)
call FIND_TS_NODE (in8,in3)
call FIND_TS_NODE (in8,in4)
call FIND_TS_NODE (in8,in5)
call FIND_TS_NODE (in8,in6)
call FIND_TS_NODE (in8,in7)
enddo

各節点において，
IALU(i,k)が小さい番号から大きい番号に並ぶようにソート（単純なバブルソート）
せいぜい100程度のものをソートする

do in= 1, N
NN= INLU(in)do k= 1, NN
NCOL1(k)= IALU(in,k)enddo
call mSORT (NCOL1, NCOL2, NN)
do k= NN, 1, -1
IALU(in,NN-k+1)= NCOL1(NCOL2(k))
enddo
enddo

FEM3D

39

CRS形式への変換：MAT_CON1!C!C***!C*** MAT_CON1!C***!C
subroutine MAT_CON1
use pfem_util
implicit REAL*8 (A-H,O-Z)
allocate (index(0:NP))
index= 0Ci

do i= 1, NP
index(i)= index(i-1) + INLU(i)enddo

index[i + 1] =  INLU[k ]k =0

index[0] = 0

NPLU= index(NP)allocate (item(NPLU))FORTRANido i= 1, NP
do k= 1, INLU(i)kk = k + index(i-1)
item(kk)=
IALU(i,k)
enddo
enddo
deallocate (INLU, IALU)
end subroutine MAT_CON1

index (i ) =  INLU(k )k =1

index (0) = 0

FEM3D

40

CRS形式への変換：MAT_CON1!C!C***!C*** MAT_CON1!C***!C
subroutine MAT_CON1
use pfem_util
implicit REAL*8 (A-H,O-Z)
allocate (index(0:NP))
index= 0
do i= 1, NP
index(i)= index(i-1) + INLU(i)enddo
NPLU= index(NP)allocate (item(NPLU))
do i= 1, NP
do k= 1, INLU(i)kk = k + index(i-1)
item(kk)=
IALU(i,k)
enddo
enddo
deallocate (INLU, IALU)
end subroutine MAT_CON1

NPLU=index(NP)itemのサイズ
非ゼロ非対角成分総数

FEM3D

41

CRS形式への変換：MAT_CON1!C!C***!C*** MAT_CON1!C***!C
subroutine MAT_CON1
use pfem_util
implicit REAL*8 (A-H,O-Z)
allocate (index(0:NP))
index= 0
do i= 1, NP
index(i)= index(i-1) + INLU(i)enddo
NPLU= index(NP)allocate (item(NPLU))
do i= 1, NP
do k= 1, INLU(i)kk = k + index(i-1)
item(kk)=
IALU(i,k)
enddo
enddo
deallocate (INLU, IALU)
end subroutine MAT_CON1

itemに1から始まる
節点番号を記憶FEM3D

42

CRS形式への変換：MAT_CON1!C!C***!C*** MAT_CON1!C***!C
subroutine MAT_CON1
use pfem_util
implicit REAL*8 (A-H,O-Z)
allocate (index(0:NP))
index= 0
do i= 1, NP
index(i)= index(i-1) + INLU(i)enddo
NPLU= index(NP)allocate (item(NPLU))
do i= 1, NP
do k= 1, INLU(i)kk = k + index(i-1)
item(kk)=
IALU(i,k)
enddo
enddo
deallocate (INLU, IALU)
end subroutine MAT_CON1

これらはもはや不要

pFEM3D-2

43

全体処理
program heat3Dp
use solver11
use pfem_util
implicit REAL*8(A-H,O-Z)
call PFEM_INIT
call INPUT_CNTL
call INPUT_GRID
call MAT_CON0
call MAT_CON1
call MAT_ASS_MAIN
call MAT_ASS_BC
call SOLVE11
call OUTPUT_UCD
call PFEM_FINALIZE
end program heat3Dp

pFEM3D-2

44

MAT_ASS_MAIN：全体構成
do kpn= 1, 2
do jpn= 1, 2
do ipn= 1, 2

ガウス積分点番号（ζ方向）
ガウス積分点番号（η方向）
ガウス積分点番号（ξ方向）
ガウス積分点（8個）における形状関数，
およびその「自然座標系」における微分の算出

enddo
enddo
enddo
do icel= 1, ICELTOT

要素ループ
8節点の座標から，ガウス積分点における，形状関数の「全体座標系」における微分，およびヤコビアンを算出（JACOBI）

do ie= 1, 8
do je= 1, 8

局所節点番号局所節点番号
全体節点番号：ip, jp
Aip,jpのitemLUにおけるアドレス：kk

do kpn= 1, 2
do jpn= 1, 2
do ipn= 1, 2

ガウス積分点番号（ζ方向）ie
ガウス積分点番号（η方向）
ガウス積分点番号（ξ方向）
要素積分⇒要素行列成分計算，全体行列への足しこみ

enddo
enddo
enddo
enddo
enddo
enddo

je

pFEM3D-245

係数行列：MAT_ASS_MAIN（1/6）!C!C***!C*** MAT_ASS_MAIN!C***!Csubroutine MAT_ASS_MAIN
use pfem_util
implicit REAL*8 (A-H,O-Z)
integer(kind=kint), dimension(8) :: nodLOCAL

allocate (AMAT(NPLU))
allocate (B(NP), D(NP), X(NP))
AMAT= 0.d0
B= 0.d0
X= 0.d0
D= 0.d0
WEI(1)= +1.0000000000D+00
WEI(2)= +1.0000000000D+00
POS(1)= -0.5773502692D+00
POS(2)= +0.5773502692D+00

係数行列（非零非対角成分）
右辺ベクトル
未知数ベクトル
係数行列（対角成分）

pFEM3D-246

係数行列：MAT_ASS_MAIN（1/6）!C!C***!C*** MAT_ASS_MAIN!C***!Csubroutine MAT_ASS_MAIN
use pfem_util
implicit REAL*8 (A-H,O-Z)
integer(kind=kint), dimension(8) :: nodLOCAL

allocate (AMAT(NPLU))
allocate (B(NP), D(NP), X(NP))
AMAT= 0.d0
B= 0.d0
X= 0.d0
D= 0.d0
WEI(1)= +1.0000000000D+00
WEI(2)= +1.0000000000D+00
POS(1)= -0.5773502692D+00
POS(2)= +0.5773502692D+00

POS： 積分点座標
WEI： 重み係数

係数行列（非零非対角成分）
右辺ベクトル
未知数ベクトル
係数行列（対角成分）

pFEM3D-247

係数行列：MAT_ASS_MAIN（2/6）!C!C-- INIT.
!C
PNQ
!C
PNE
!C
PNT
!C

- 1st-order derivative of shape function by QSI- 1st-order derivative of shape function by ETA- 1st-order derivative of shape function by ZET

do kp= 1, 2
do jp= 1, 2
do ip= 1, 2
QP1= 1.d0 + POS(ip)QM1= 1.d0 - POS(ip)EP1= 1.d0 + POS(jp)EM1= 1.d0 - POS(jp)TP1= 1.d0 + POS(kp)TM1= 1.d0 - POS(kp)
SHAPE(ip,jp,kp,1)= O8th * QM1 * EM1 * TM1
SHAPE(ip,jp,kp,2)= O8th * QP1 * EM1 * TM1
SHAPE(ip,jp,kp,3)= O8th * QP1 * EP1 * TM1
SHAPE(ip,jp,kp,4)= O8th * QM1 * EP1 * TM1
SHAPE(ip,jp,kp,5)= O8th * QM1 * EM1 * TP1
SHAPE(ip,jp,kp,6)= O8th * QP1 * EM1 * TP1
SHAPE(ip,jp,kp,7)= O8th * QP1 * EP1 * TP1
SHAPE(ip,jp,kp,8)= O8th * QP1 * EP1 * TP1

pFEM3D-248

係数行列：MAT_ASS_MAIN（2/6）!C!C-- INIT.
!C
PNQ
!C
PNE
!C
PNT
!C

- 1st-order derivative of shape function by QSI- 1st-order derivative of shape function by ETA- 1st-order derivative of shape function by ZET

do kp= 1, 2
do jp= 1, 2
do ip= 1, 2
QP1= 1.d0 + POS(ip)QM1= 1.d0 - POS(ip)EP1= 1.d0 + POS(jp)EM1= 1.d0 - POS(jp)TP1= 1.d0 + POS(kp)TM1= 1.d0 - POS(kp)
SHAPE(ip,jp,kp,1)= O8th * QM1 * EM1 * TM1
SHAPE(ip,jp,kp,2)= O8th * QP1 * EM1 * TM1
SHAPE(ip,jp,kp,3)= O8th * QP1 * EP1 * TM1
SHAPE(ip,jp,kp,4)= O8th * QM1 * EP1 * TM1
SHAPE(ip,jp,kp,5)= O8th * QM1 * EM1 * TP1
SHAPE(ip,jp,kp,6)= O8th * QP1 * EM1 * TP1
SHAPE(ip,jp,kp,7)= O8th * QP1 * EP1 * TP1
SHAPE(ip,jp,kp,8)= O8th * QP1 * EP1 * TP1

QP1(i ) = (1 + ξ i ), QM1(i ) = (1 − ξ i )EP1( j) = (1 + η j ), EM1( j) = (1 − ηi )TP1(k ) = (1 + ζ k ), TM1(k ) = (1 − ζ k )pFEM3D-249

係数行列：MAT_ASS_MAIN（2/6）!C!C-- INIT.
!C
PNQ
!C
PNE
!C
PNT
!C

- 1st-order derivative of shape function by QSI- 1st-order derivative of shape function by ETA- 1st-order derivative of shape function by ZET

do kp= 1, 2
do jp= 1, 2
do ip= 1, 2
QP1= 1.d0 + POS(ip)QM1= 1.d0 - POS(ip)EP1= 1.d0 + POS(jp)EM1= 1.d0 - POS(jp)TP1= 1.d0 + POS(kp)TM1= 1.d0 - POS(kp)
SHAPE(ip,jp,kp,1)= O8th * QM1 * EM1 * TM1
SHAPE(ip,jp,kp,2)= O8th * QP1 * EM1 * TM1
SHAPE(ip,jp,kp,3)= O8th * QP1 * EP1 * TM1
SHAPE(ip,jp,kp,4)= O8th * QM1 * EP1 * TM1
SHAPE(ip,jp,kp,5)= O8th * QM1 * EM1 * TP1
SHAPE(ip,jp,kp,6)= O8th * QP1 * EM1 * TP1
SHAPE(ip,jp,kp,7)= O8th * QP1 * EP1 * TP1
SHAPE(ip,jp,kp,8)= O8th * QP1 * EP1 * TP1(− 1,+1,+1) 8(− 1,−1,+1) 5

7 (+ 1,+1,+1)
6 (+ 1,−1,+1)(− 1,+1,−1)
4
1(ξ ,η , ζ ) = (− 1,−1,−1)

2

3 (+ 1,+1,−1)(+ 1,−1,−1)pFEM3D-250

係数行列：MAT_ASS_MAIN（2/6）!C!C-- INIT.
!C
PNQ
!C
PNE
!C
PNT
!C

- 1st-order derivative of shape function by QSI- 1st-order derivative of shape function by ETA- 1st-order derivative of shape function by ZET

do kp= 1, 2
do jp= 1, 2
do ip= 1, 2
QP1= 1.d0 + POS(ip)QM1= 1.d0 - POS(ip)EP1= 1.d0 + POS(jp)EM1= 1.d0 - POS(jp)TP1= 1.d0 + POS(kp)TM1= 1.d0 - POS(kp)
SHAPE(ip,jp,kp,1)= O8th * QM1 * EM1 * TM1
SHAPE(ip,jp,kp,2)= O8th * QP1 * EM1 * TM1
SHAPE(ip,jp,kp,3)= O8th * QP1 * EP1 * TM1
SHAPE(ip,jp,kp,4)= O8th * QM1 * EP1 * TM1
SHAPE(ip,jp,kp,5)= O8th * QM1 * EM1 * TP1
SHAPE(ip,jp,kp,6)= O8th * QP1 * EM1 * TP1
SHAPE(ip,jp,kp,7)= O8th * QP1 * EP1 * TP1
SHAPE(ip,jp,kp,8)= O8th * QP1 * EP1 * TP1

1(1 − ξ )(1 − η )(1 − ζ )81
N 2 (ξ ,η , ζ ) = (1 + ξ )(1 − η )(1 − ζ )81
N 3 (ξ ,η , ζ ) = (1 + ξ )(1 + η )(1 − ζ )81
N 4 (ξ ,η , ζ ) = (1 − ξ )(1 + η )(1 − ζ )8
N1 (ξ ,η , ζ ) =1(1 − ξ )(1 − η )(1 + ζ )81
N 6 (ξ ,η , ζ ) = (1 + ξ )(1 − η )(1 + ζ )81
N 7 (ξ ,η , ζ ) = (1 + ξ )(1 + η )(1 + ζ )81
N 8 (ξ ,η , ζ ) = (1 − ξ )(1 + η )(1 + ζ )8
N 5 (ξ ,η , ζ ) =pFEM3D-251

係数行列：MAT_ASS_MAIN（3/6）
PNQ(jp,kp,1)= - O8th * EM1 * TM1PNQ(jp,kp,2)= + O8th * EM1 * TM1PNQ(jp,kp,3)= + O8th * EP1 * TM1PNQ(jp,kp,4)= - O8th * EP1 * TM1PNQ(jp,kp,5)= - O8th * EM1 * TP1PNQ(jp,kp,6)= + O8th * EM1 * TP1PNQ(jp,kp,7)= + O8th * EP1 * TP1PNQ(jp,kp,8)= - O8th * EP1 * TP1
PNE(ip,kp,1)= - O8th * QM1 * TM1
PNE(ip,kp,2)= - O8th * QP1 * TM1
PNE(ip,kp,3)= + O8th * QP1 * TM1
PNE(ip,kp,4)= + O8th * QM1 * TM1
PNE(ip,kp,5)= - O8th * QM1 * TP1
PNE(ip,kp,6)= - O8th * QP1 * TP1
PNE(ip,kp,7)= + O8th * QP1 * TP1
PNE(ip,kp,8)= + O8th * QM1 * TP1PNT(ip,jp,1)= - O8th * QM1 * EM1PNT(ip,jp,2)= - O8th * QP1 * EM1PNT(ip,jp,3)= - O8th * QP1 * EP1PNT(ip,jp,4)= - O8th * QM1 * EP1PNT(ip,jp,5)= + O8th * QM1 * EM1PNT(ip,jp,6)= + O8th * QP1 * EM1PNT(ip,jp,7)= + O8th * QP1 * EP1PNT(ip,jp,8)= + O8th * QM1 * EP1
enddo
enddo
enddo
do icel= 1, ICELTOTCOND0= COND
in1= ICELNOD(icel,1)
in2= ICELNOD(icel,2)
in3= ICELNOD(icel,3)
in4= ICELNOD(icel,4)
in5= ICELNOD(icel,5)
in6= ICELNOD(icel,6)
in7= ICELNOD(icel,7)
in8= ICELNOD(icel,8)∂N l(ξ = ξ i ,η = η j , ζ = ζ k )∂ξ
∂N
PNE (i, k ) = l (ξ = ξ i ,η = η j , ζ = ζ k )∂η
∂N
PNT (i, j ) = l (ξ = ξ i ,η = η j , ζ = ζ k )∂ζ
PNQ( j , k ) =∂N11(ξ i ,η j , ζ k ) = − (1 − η j )(1 − ζ k )
∂ξ
8
∂N 2
1
(ξ i ,η j , ζ k ) = + (1 − η j )(1 − ζ k )
∂ξ
8
∂N 3
1
(ξ i ,η j , ζ k ) = + (1 + η j )(1 − ζ k )
∂ξ
8
∂N 3
1
(ξ i ,η j , ζ k ) = − (1 + η j )(1 − ζ k )∂ξ
8(ξ i ,η j , ζ k ) における形状関数の一階微分

pFEM3D-252係数行列：MAT_ASS_MAIN（3/6）
PNQ(jp,kp,1)= - O8th * EM1 * TM1PNQ(jp,kp,2)= + O8th * EM1 * TM1PNQ(jp,kp,3)= + O8th * EP1 * TM1PNQ(jp,kp,4)= - O8th * EP1 * TM1PNQ(jp,kp,5)= - O8th * EM1 * TP1PNQ(jp,kp,6)= + O8th * EM1 * TP1PNQ(jp,kp,7)= + O8th * EP1 * TP1PNQ(jp,kp,8)= - O8th * EP1 * TP1
PNE(ip,kp,1)= - O8th * QM1 * TM1
PNE(ip,kp,2)= - O8th * QP1 * TM1
PNE(ip,kp,3)= + O8th * QP1 * TM1
PNE(ip,kp,4)= + O8th * QM1 * TM1
PNE(ip,kp,5)= - O8th * QM1 * TP1
PNE(ip,kp,6)= - O8th * QP1 * TP1
PNE(ip,kp,7)= + O8th * QP1 * TP1
PNE(ip,kp,8)= + O8th * QM1 * TP1PNT(ip,jp,1)= - O8th * QM1 * EM1PNT(ip,jp,2)= - O8th * QP1 * EM1PNT(ip,jp,3)= - O8th * QP1 * EP1PNT(ip,jp,4)= - O8th * QM1 * EP1PNT(ip,jp,5)= + O8th * QM1 * EM1PNT(ip,jp,6)= + O8th * QP1 * EM1PNT(ip,jp,7)= + O8th * QP1 * EP1PNT(ip,jp,8)= + O8th * QM1 * EP1
enddo
enddo
enddo
do icel= 1, ICELTOTCOND0= COND
in1= ICELNOD(icel,1)
in2= ICELNOD(icel,2)
in3= ICELNOD(icel,3)
in4= ICELNOD(icel,4)
in5= ICELNOD(icel,5)
in6= ICELNOD(icel,6)
in7= ICELNOD(icel,7)
in8= ICELNOD(icel,8)(− 1,+1,+1) 8(− 1,−1,+1) 5

7 (+ 1,+1,+1)
6 (+ 1,−1,+1)(− 1,+1,−1)
4
1(ξ ,η , ζ ) = (− 1,−1,−1)

2

3 (+ 1,+1,−1)(+ 1,−1,−1)pFEM3D-253

係数行列：MAT_ASS_MAIN（4/6）
nodLOCAL(1)= in1
nodLOCAL(2)= in2
nodLOCAL(3)= in3
nodLOCAL(4)= in4
nodLOCAL(5)= in5
nodLOCAL(6)= in6
nodLOCAL(7)= in7
nodLOCAL(8)= in8

&

&
&
&

8節点の節点番号

X1= XYZ(in1,1)
X2= XYZ(in2,1)
X3= XYZ(in3,1)
X4= XYZ(in4,1)
X5= XYZ(in5,1)
X6= XYZ(in6,1)
X7= XYZ(in7,1)
X8= XYZ(in8,1)
Y1= XYZ(in1,2)
Y2= XYZ(in2,2)
Y3= XYZ(in3,2)
Y4= XYZ(in4,2)
Y5= XYZ(in5,2)
Y6= XYZ(in6,2)
Y7= XYZ(in7,2)
Y8= XYZ(in8,2)
QVC= O8th * (X1+X2+X3+X4+X5+X6+X7+X8+
Y1+Y2+Y3+Y4+Y5+Y6+Y7+Y8)
Z1= XYZ(in1,3)
Z2= XYZ(in2,3)
Z3= XYZ(in3,3)
Z4= XYZ(in4,3)
Z5= XYZ(in5,3)
Z6= XYZ(in6,3)
Z7= XYZ(in7,3)
Z8= XYZ(in8,3)call JACOBI (DETJ, PNQ, PNE, PNT, PNX, PNY, PNZ,
X1, X2, X3, X4, X5, X6, X7, X8,
Y1, Y2, Y3, Y4, Y5, Y6, Y7, Y8,
Z1, Z2, Z3, Z4, Z5, Z6, Z7, Z8 )(− 1,+1,+1) 8(− 1,−1,+1) 5

7 (+ 1,+1,+1)
6 (+ 1,−1,+1)(− 1,+1,−1)
4
1(ξ ,η , ζ ) = (− 1,−1,−1)

&
&
&

2

3 (+ 1,+1,−1)(+ 1,−1,−1)pFEM3D-254

係数行列：MAT_ASS_MAIN（4/6）
nodLOCAL(1)= in1
nodLOCAL(2)= in2
nodLOCAL(3)= in3
nodLOCAL(4)= in4
nodLOCAL(5)= in5
nodLOCAL(6)= in6
nodLOCAL(7)= in7
nodLOCAL(8)= in8

&

&
&
&

X1= XYZ(in1,1)
X2= XYZ(in2,1)
X3= XYZ(in3,1)
X4= XYZ(in4,1)
8節点のX座標
X5= XYZ(in5,1)
X6= XYZ(in6,1)
X7= XYZ(in7,1)
X8= XYZ(in8,1)
Y1= XYZ(in1,2)
Y2= XYZ(in2,2)
Y3= XYZ(in3,2)
Y4= XYZ(in4,2)
8節点のY座標
Y5= XYZ(in5,2)
Y6= XYZ(in6,2)
Y7= XYZ(in7,2)
Y8= XYZ(in8,2)
QVC= O8th * (X1+X2+X3+X4+X5+X6+X7+X8+
Y1+Y2+Y3+Y4+Y5+Y6+Y7+Y8)
Z1= XYZ(in1,3)
Z2= XYZ(in2,3)
Z3= XYZ(in3,3)
Z4= XYZ(in4,3)
8節点のZ座標
Z5= XYZ(in5,3)
Z6= XYZ(in6,3)
Z7= XYZ(in7,3)
Z8= XYZ(in8,3)call JACOBI (DETJ, PNQ, PNE, PNT, PNX, PNY, PNZ,
X1, X2, X3, X4, X5, X6, X7, X8,
Y1, Y2, Y3, Y4, Y5, Y6, Y7, Y8,
Z1, Z2, Z3, Z4, Z5, Z6, Z7, Z8 )(− 1,+1,+1) 8(− 1,−1,+1) 5

7 (+ 1,+1,+1)
6 (+ 1,−1,+1)(− 1,+1,−1)
4
1(ξ ,η , ζ ) = (− 1,−1,−1)

&
&
&

2

3 (+ 1,+1,−1)(+ 1,−1,−1)pFEM3D-255

係数行列：MAT_ASS_MAIN（4/6）
nodLOCAL(1)= in1
nodLOCAL(2)= in2
nodLOCAL(3)= in3
nodLOCAL(4)= in4
nodLOCAL(5)= in5
nodLOCAL(6)= in6
nodLOCAL(7)= in7
nodLOCAL(8)= in8
X1= XYZ(in1,1)
X2= XYZ(in2,1)
X3= XYZ(in3,1)
X4= XYZ(in4,1)
8節点のX座標
X5= XYZ(in5,1)
X6= XYZ(in6,1)
X7= XYZ(in7,1)
X8= XYZ(in8,1)
Y1= XYZ(in1,2)
Y2= XYZ(in2,2)
Y3= XYZ(in3,2)
Y4= XYZ(in4,2)
8節点のY座標
Y5= XYZ(in5,2)
Y6= XYZ(in6,2)
Y7= XYZ(in7,2)
Y8= XYZ(in8,2)
QVC= O8th * (X1+X2+X3+X4+X5+X6+X7+X8+&
Y1+Y2+Y3+Y4+Y5+Y6+Y7+Y8)
Z1= XYZ(in1,3)
Z2= XYZ(in2,3)
Z3= XYZ(in3,3)
Z4= XYZ(in4,3)
Z5= XYZ(in5,3)
Z6= XYZ(in6,3)
Z7= XYZ(in7,3)
Z8= XYZ(in8,3)&&&(− 1,+1,+1) 8

7 (+ 1,+1,+1)(− 1,−1,+1) 5

6 (+ 1,−1,+1)(− 1,+1,−1)
4
1(ξ ,η , ζ ) = (− 1,−1,−1)

2

3 (+ 1,+1,−1)(+ 1,−1,−1)∂  ∂T  ∂  ∂T  ∂  ∂T  ɺ +  λλ +  λ + Q ( x, y , z ) = 0∂x  ∂x  ∂y  ∂y  ∂z  ∂z Qɺ (x, y, z ) = QVOL xC + yC

call JACOBI (DETJ, PNQ, PNE, PNT, PNX, PNY, PNZ,&
X1, X2, X3, X4, X5, X6, X7, X8, 体積当たり発熱量は位置（メッシュの中心&
Y1, Y2, Y3, Y4, Y5, Y6, Y7, Y8,&
Z1, Z2, Z3, Z4, Z5, Z6, Z7, Z8 ) の座標 xc,yc）に依存

pFEM3D-256

係数行列：MAT_ASS_MAIN（4/6）
nodLOCAL(1)= in1
nodLOCAL(2)= in2
nodLOCAL(3)= in3
nodLOCAL(4)= in4
nodLOCAL(5)= in5
nodLOCAL(6)= in6
nodLOCAL(7)= in7
nodLOCAL(8)= in8
X1= XYZ(in1,1)
X2= XYZ(in2,1)
X3= XYZ(in3,1)
X4= XYZ(in4,1)
X5= XYZ(in5,1)
X6= XYZ(in6,1)
X7= XYZ(in7,1)
X8= XYZ(in8,1)
Y1= XYZ(in1,2)
Y2= XYZ(in2,2)
Y3= XYZ(in3,2)
Y4= XYZ(in4,2)
Y5= XYZ(in5,2)
Y6= XYZ(in6,2)
Y7= XYZ(in7,2)
Y8= XYZ(in8,2)
QVC= O8th * (X1+X2+X3+X4+X5+X6+X7+X8+&
Y1+Y2+Y3+Y4+Y5+Y6+Y7+Y8)
Z1= XYZ(in1,3)
Z2= XYZ(in2,3)
Z3= XYZ(in3,3)
Z4= XYZ(in4,3)
Z5= XYZ(in5,3)
Z6= XYZ(in6,3)
Z7= XYZ(in7,3)
Z8= XYZ(in8,3)&&&(− 1,+1,+1) 8

7 (+ 1,+1,+1)(− 1,−1,+1) 5

6 (+ 1,−1,+1)(− 1,+1,−1)
4
1(ξ ,η , ζ ) = (− 1,−1,−1)

2

3 (+ 1,+1,−1)(+ 1,−1,−1)∂  ∂T  ∂  ∂T  ∂  ∂T  ɺ +  λλ +  λ + Q ( x, y , z ) = 0∂x  ∂x  ∂y  ∂y  ∂z  ∂z Qɺ (x, y, z ) = QVOL xC + yC

call JACOBI (DETJ, PNQ, PNE, PNT, PNX, PNY, PNZ,&
QVC = xC + y&C
X1, X2, X3, X4, X5, X6, X7, X8,
Y1, Y2, Y3, Y4, Y5, Y6, Y7, Y8,&
Z1, Z2, Z3, Z4, Z5, Z6, Z7, Z8 )

pFEM3D-257

係数行列：MAT_ASS_MAIN（4/6）
nodLOCAL(1)= in1
nodLOCAL(2)= in2
nodLOCAL(3)= in3
nodLOCAL(4)= in4
nodLOCAL(5)= in5
nodLOCAL(6)= in6
nodLOCAL(7)= in7
nodLOCAL(8)= in8

&

&
&
&

X1= XYZ(in1,1)
X2= XYZ(in2,1)
X3= XYZ(in3,1)
X4= XYZ(in4,1)
X5= XYZ(in5,1)
X6= XYZ(in6,1)
X7= XYZ(in7,1)
X8= XYZ(in8,1)
Y1= XYZ(in1,2)
Y2= XYZ(in2,2)
Y3= XYZ(in3,2)
Y4= XYZ(in4,2)
Y5= XYZ(in5,2)
Y6= XYZ(in6,2)
Y7= XYZ(in7,2)
Y8= XYZ(in8,2)
QVC= O8th * (X1+X2+X3+X4+X5+X6+X7+X8+
Y1+Y2+Y3+Y4+Y5+Y6+Y7+Y8)
Z1= XYZ(in1,3)
Z2= XYZ(in2,3)
Z3= XYZ(in3,3)
Z4= XYZ(in4,3)
Z5= XYZ(in5,3)
Z6= XYZ(in6,3)
Z7= XYZ(in7,3)
Z8= XYZ(in8,3)call JACOBI (DETJ, PNQ, PNE, PNT, PNX, PNY, PNZ,
X1, X2, X3, X4, X5, X6, X7, X8,
Y1, Y2, Y3, Y4, Y5, Y6, Y7, Y8,
Z1, Z2, Z3, Z4, Z5, Z6, Z7, Z8 )&&&

pFEM3D-258

係数行列：MAT_ASS_MAIN（5/6）!C!C== CONSTRUCT the GLOBAL MATRIX
do ie= 1, 8
ip = nodLOCAL(ie)do je= 1, 8
jp = nodLOCAL(je)

全体行列の非対角成分

Aip , jp

kk= 0
if (jp.ne.ip) theniiS= index(ip-1) + 1
iiE= index(ip )do k= iiS, iiE
if ( item(k).eq.jp ) then
kk= k
exit
endif
enddo
endif

(− 1,+1,+1) 8(− 1,−1,+1) 5

6 (+ 1,−1,+1)
4(ξ ,η , ζ ) = (− 1,−1,−1)

ip= nodLOCAL(ie)jp= nodLOCAL(je)
7 (+ 1,+1,+1)(− 1,+1,−1)1kk：itemにおけるアドレス2

3 (+ 1,+1,−1)(+ 1,−1,−1)pFEM3D-2

59

要素マトリクス：8×8行列
j[k ] (i, j = 1…8)ij

i(− 1,+1,+1) 8(− 1,−1,+1) 5

7 (+ 1,+1,+1)
6 (+ 1,−1,+1)(− 1,+1,−1)
4
1(ξ ,η , ζ ) = (− 1,−1,−1)

2

3 (+ 1,+1,−1)(+ 1,−1,−1)pFEM3D-260

係数行列：MAT_ASS_MAIN（5/6）!C!C== CONSTRUCT the GLOBAL MATRIX
do ie= 1, 8
ip = nodLOCAL(ie)do je= 1, 8
jp = nodLOCAL(je)要素マトリクス（ie~je）
全体マトリクス（ip~jp）の関係

kk= 0
if (jp.ne.ip) theniiS= index(ip-1) + 1
iiE= index(ip )do k= iiS, iiE
if ( item(k).eq.jp ) then
kk= k
exit
endif
enddo
endif

(− 1,+1,+1) 8(− 1,−1,+1) 5

7 (+ 1,+1,+1)
6 (+ 1,−1,+1)(− 1,+1,−1)
4
1(ξ ,η , ζ ) = (− 1,−1,−1)kk：itemにおけるアドレス2

3 (+ 1,+1,−1)(+ 1,−1,−1)pFEM3D-261

係数行列：MAT_ASS_MAIN（6/6）
QV0 = 0.d0
COEFij= 0.d0
do kpn= 1, 2
do jpn= 1, 2
do ipn= 1, 2
coef= dabs(DETJ(ipn,jpn,kpn))*WEI(ipn)*WEI(jpn)*WEI(kpn)
PNXi= PNX(ipn,jpn,kpn,ie)
PNYi= PNY(ipn,jpn,kpn,ie)
PNZi= PNZ(ipn,jpn,kpn,ie)
PNXj= PNX(ipn,jpn,kpn,je)
PNYj= PNY(ipn,jpn,kpn,je)
PNZj= PNZ(ipn,jpn,kpn,je)
COEFij= COEFij + coef * COND0 *(PNXi*PNXj+PNYi*PNYj+PNZi*PNZj)&

SHi= SHAPE(ipn,jpn,kpn,ie)
QV0= QV0 + SHi * QVOL * coef
enddo
enddo
enddo
if (jp.eq.ip) thenD(ip)= D(ip) + COEFij
B(ip)= B(ip) + QV0*QVC
else
AMAT(kk)= AMAT(kk) + COEFij
endif
enddo
+1 +1 +1
enddo
 ∂N i ∂N j
∂N i ∂N j∂N i ∂N j 
enddo
−
+λ
+λ
λ
 det J
returnend  

−1 −1 −1

∂x ∂x

∂y ∂y

∂z ∂z dξdηdζ

pFEM3D-262

係数行列：MAT_ASS_MAIN（6/6）
QV0 = 0.d0
COEFij= 0.d0
do kpn= 1, 2
do jpn= 1, 2
do ipn= 1, 2
coef= dabs(DETJ(ipn,jpn,kpn))*WEI(ipn)*WEI(jpn)*WEI(kpn)
PNXi= PNX(ipn,jpn,kpn,ie)
PNYi= PNY(ipn,jpn,kpn,ie)
PNZi= PNZ(ipn,jpn,kpn,ie)
PNXj= PNX(ipn,jpn,kpn,je)
PNYj= PNY(ipn,jpn,kpn,je)
PNZj= PNZ(ipn,jpn,kpn,je)
COEFij= COEFij + coef * COND0 *(PNXi*PNXj+PNYi*PNYj+PNZi*PNZj)&

SHi= SHAPE(ipn,jpn,kpn,ie)
QV0= QV0 + SHi * QVOL * coef
enddo
enddo
enddo

+1 +1 +1

I=

  f (ξ ,η , ζ ) dξdηdζ

−1 −1 −1  [L
M
N
if (jp.eq.ip) thenD(ip)= D(ip) + COEFij=
Wi ⋅ W j ⋅Wk ⋅ f (ξ i ,η j , ζ k )B(ip)= B(ip) + QV0*QVC
i =1
j =1 k =1
else
AMAT(kk)= AMAT(kk) + COEFij
endif
enddo
+1 +1 +1
enddo
 ∂N i ∂N j
∂N i ∂N j∂N i ∂N j 
enddo
−
+λ
+λ
λ
 det J dξdηdζreturnend  

−1 −1 −1

∂x ∂x

∂y ∂y

∂z ∂z ]

FEM3D63

MAT_ASS_MAIN (6/6)
do kpn= 1, 2
do jpn= 1, 2
do ipn= 1, 2
coef= dabs(DETJ(ipn,jpn,kpn))*WEI(ipn)*WEI(jpn)*WEI(kpn)
PNXi= PNX(ipn,jpn,kpn,ie)
PNYi= PNY(ipn,jpn,kpn,ie)
PNZi= PNZ(ipn,jpn,kpn,ie)
PNXj= PNX(ipn,jpn,kpn,je)
PNYj= PNY(ipn,jpn,kpn,je)
PNZj= PNZ(ipn,jpn,kpn,je)

coef = Wi ⋅ W j ⋅ Wk ⋅ det J (ξ i ,η j , ζ k )+1 +1 +1

  f (ξ ,η , ζ ) dξdηdζ

I=

−1 −1 −1
L

M

N

i =1

j =1

k =1

[

=    Wi ⋅ W j ⋅ Wk ⋅ f (ξ i ,η j , ζ k )COEFij= COEFij + coef * COND0 * (PNXi*PNXj+PNYi*PNYj+PNZi*PNZj)
enddo
enddo
enddo
 ∂N i ∂N j
∂N i ∂N j∂N i ∂N j 
+
+
λ
λ
λ
 det J dξdηdζ
−1 −1 −1  ∂x ∂x∂y ∂y∂z ∂z +1 +1 +1]

FEM3D64

MAT_ASS_MAIN (6/6)
do kpn= 1, 2
do jpn= 1, 2
do ipn= 1, 2
coef= dabs(DETJ(ipn,jpn,kpn))*WEI(ipn)*WEI(jpn)*WEI(kpn)
PNXi= PNX(ipn,jpn,kpn,ie)
PNYi= PNY(ipn,jpn,kpn,ie)
PNZi= PNZ(ipn,jpn,kpn,ie)
PNXj= PNX(ipn,jpn,kpn,je)
PNYj= PNY(ipn,jpn,kpn,je)
PNZj= PNZ(ipn,jpn,kpn,je)

coef = Wi ⋅ W j ⋅ Wk ⋅ det J (ξ i ,η j , ζ k )+1 +1 +1

  f (ξ ,η , ζ ) dξdηdζ

I=

−1 −1 −1
L

M

N

i =1

j =1

k =1

[

=    Wi ⋅ W j ⋅ Wk ⋅ f (ξ i ,η j , ζ k )COEFij= COEFij + coef * COND0 * (PNXi*PNXj+PNYi*PNYj+PNZi*PNZj)
enddo
enddo
enddo
 ∂N i ∂N j
∂N i ∂N j∂N i ∂N j 
+
+
λ
λ
λ
 det J dξdηdζ
−1 −1 −1  ∂x ∂x∂y ∂y∂z ∂z +1 +1 +1]

pFEM3D-265

係数行列：MAT_ASS_MAIN（6/6）
QV0= 0.d0
COEFij= 0.d0
do kpn= 1, 2
do jpn= 1, 2
do ipn= 1, 2
coef= dabs(DETJ(ipn,jpn,kpn))*WEI(ipn)*WEI(jpn)*WEI(kpn)
PNXi= PNX(ipn,jpn,kpn,ie)
PNYi= PNY(ipn,jpn,kpn,ie)
PNZi= PNZ(ipn,jpn,kpn,ie)
PNXj= PNX(ipn,jpn,kpn,je)
PNYj= PNY(ipn,jpn,kpn,je)
PNZj= PNZ(ipn,jpn,kpn,je)[k ] (i, j = 1…8)

COEFij= COEFij + coef * COND0 *(PNXi*PNXj+PNYi*PNYj+PNZi*PNZj)

&

ij

j

SHi= SHAPE(ipn,jpn,kpn,ie)
QV0= QV0 + SHi * QVOL * coef
enddo
enddo
enddo
if (jp.eq.ip) thenD(ip)= D(ip) + COEFij
B(ip)= B(ip) + QV0*QVC
else
AMAT(kk)= AMAT(kk) + COEFij
endif
enddo
enddo
enddo
returnend

ipFEM3D-266

係数行列：MAT_ASS_MAIN（6/6）
QV0 = 0.d0
COEFij= 0.d0
do kpn= 1, 2
do jpn= 1, 2
do ipn= 1, 2
coef= dabs(DETJ(ipn,jpn,kpn))*WEI(ipn)*WEI(jpn)*WEI(kpn)
PNXi= PNX(ipn,jpn,kpn,ie)
PNYi= PNY(ipn,jpn,kpn,ie)
PNZi= PNZ(ipn,jpn,kpn,ie)
PNXj= PNX(ipn,jpn,kpn,je)
PNYj= PNY(ipn,jpn,kpn,je)
PNZj= PNZ(ipn,jpn,kpn,je)
COEFij= COEFij + coef * COND0 *(PNXi*PNXj+PNYi*PNYj+PNZi*PNZj)&

SHi= SHAPE(ipn,jpn,kpn,ie)
QV0= QV0 + SHi * QVOL * coef
enddo
enddo
enddo
if (jp.eq.ip) thenD(ip)= D(ip) + COEFij
B(ip)= B(ip) + QV0*QVC
else
AMAT(kk)= AMAT(kk) + COEFij
endif
enddo
enddo
enddo
returnend[k ](e) {φ }(e) = { f }(e)[ f ](e) =  Qɺ [N ]T dVVQɺ ( x, y, z ) = QVOL xC + yC

QVC = xC + yC

QV 0 =  QVOL[N ] dVTV

[ f ](e) = QV 0 ⋅ QVC

pFEM3D-2

67

MAT_ASS_BC：全体構成
do i= 1, NP

節点ループ（ディリクレ）境界条件を設定する節点をマーク（IWKX）

enddo
節点ループ
do i= 1, NP
if (IWKX(i,1).eq.1) then

マークされた節点だったら対応する右辺ベクトル（B）の成分，対角成分（D）の成分の修正（行・列）Z

do k= index(i-1)+1, index(i)対応する非零非対角成分（AMAT）の成分の修正（行）

enddo
endif
enddo

T=0@Z=zmax

do i= 1, NP 節点ループ
do k= index (i-1)+1, index (i)if (IWKX(item (k),1).eq.1) then

NZ

対応する非零非対角成分の節点がマークされていたら対応する右辺ベクトル，非零非対角成分（AMAT）の成分の修正（列）

endif
enddo
enddo

Y
NX
NY
X

pFEM3D-2

68

境界条件：MAT_ASS_BC（1/2）
subroutine MAT_ASS_BC
use pfem_util
implicit REAL*8 (A-H,O-Z)
allocate (IWKX(NP,2))
IWKX= 0!C
!C== Z=Zmax
do in= 1, NP
IWKX(in,1)= 0
enddo

節点グループ名が「Zmax」である
節点inにおいて：
IWKX(in,1)= 1とする
ib0= -1
do ib0= 1, NODGRPtot
if (NODGRP_NAME(ib0).eq.'Zmax') exit
enddo
do ib= NODGRP_INDEX(ib0-1)+1, NODGRP_INDEX(ib0)in= NODGRP_ITEM(ib)
IWKX(in,1)= 1
enddopFEM3D-2

69

境界条件：MAT_ASS_BC（2/2）
do in= 1, NP
if (IWKX(in,1).eq.1) thenB(in)= 0.d0
D(in)= 1.d0
iS= index(in-1) + 1
iE= index(in )do k= iS, iE
AMAT(k)= 0.d0
enddo
endif
enddo

!C==

do in= 1, NP
iS= index(in-1) + 1
iE= index(in )do k= iS, iE
if (IWKX(item(k),1).eq.1) then
AMAT(k)= 0.d0
endif
enddo
enddo
returnend

pFEM3D-270

境界条件：MAT_ASS_BC（2/2）
do in= 1, NP
if (IWKX(in,1).eq.1) thenB(in)= 0.d0
D(in)= 1.d0
iS= index(in-1) + 1
iE= index(in )do k= iS, iE
AMAT(k)= 0.d0
enddo
endif
enddo

!C==

ゼロクリア !!

do in= 1, NP
iS= index(in-1) + 1
iE= index(in )do k= iS, iE
if (IWKX(item(k),1).eq.1) then
AMAT(k)= 0.d0
endif
enddo
enddo
returnend

IWKX(in,1)=1となる節点に対して対角成分=1，右辺=0，
非零対角成分=0

ここでやっていることも1CPUの時と全く変わらないpFEM3D-2

71

境界条件：MAT_ASS_BC（2/2）
do in= 1, NP
if (IWKX(in,1).eq.1) thenB(in)= 0.d0
D(in)= 1.d0
iS= index(in-1) + 1
iE= index(in )do k= iS, iE
AMAT(k)= 0.d0
enddo
endif
enddo

!C==

IWKX(in,1)=1となる節点に対して対角成分=1，右辺=0，
非零対角成分=0
消去，
ゼロクリア

do in= 1, NP
iS= index(in-1) + 1
iE= index(in )do k= iS, iE
if (IWKX(item(k),1).eq.1) then
AMAT(k)= 0.d0
endif
enddo
enddo
returnend

ここでやっていることも1CPUの時と全く変わらないpFEM3D-2

72

並列有限要素法の処理：プログラム• 初期化– 制御変数読み込み– 座標読み込み⇒要素生成（N:節点数，NE：要素数）– 配列初期化（全体マトリクス，要素マトリクス）– 要素⇒全体マトリクスマッピング（Index，Item）• マトリクス生成– 要素単位の処理（do icel= 1, NE）• 要素マトリクス計算• 全体マトリクスへの重ね合わせ– 境界条件の処理• 連立一次方程式– 共役勾配法（CG）

pFEM3D-2

73

test1

input_cntl

メインプログラム制御データ入力

input_griddefine_file_name

メッシュファイル入力

局所ファイル名

mat_con0mSORT

行列コネクティビティ生成

ソート

mat_con1

find_node

行列コネクティビティ生成

節点探索

mat_ass_main

jacobi係数行列生成ヤコビアン計算

mat_ass_bc
境界条件処理三次元並列熱伝
導解析コードheat3Dp

solve11cg

線形ソルバー制御

CG法計算

output_ucd
可視化処理

pFEM3D-2

74

全体処理
program heat3Dp
use solver11
use pfem_util
implicit REAL*8(A-H,O-Z)
call PFEM_INIT
call INPUT_CNTL
call INPUT_GRID
call MAT_CON0
call MAT_CON1
call MAT_ASS_MAIN
call MAT_ASS_BC
call SOLVE11
call OUTPUT_UCD
call PFEM_FINALIZE
end program heat3Dp

pFEM3D-2

75

SOLVE11
module SOLVER11
contains
subroutine SOLVE11
use pfem_util
use solver_CG
implicit REAL*8 (A-H,O-Z)
integer :: ERROR, ICFLAG
character(len=char_length) :: BUF
data ICFLAG/0/!C!C +------------+!C | PARAMETERs |!C +------------+
!C===
ITER
= pfemIarray(1)CG法の最大反復回数
RESID= pfemRarray(1)CG法の収束判定値!C===!C!C +------------------+!C | ITERATIVE solver |!C +------------------+!C===
call CG&( N, NP, NPLU, D, AMAT, index, item, B, X, RESID,&
ITER, ERROR, my_rank,&
NEIBPETOT, NEIBPE, IMPORT_INDEX, IMPORT_ITEM,&EXPORT_INDEX, EXPORT_ITEM)
ITERactual= ITER!C===
end subroutine SOLVE11
end module SOLVER11

&
&
&
&

pFEM3D-2

76

前処理付き共役勾配法
Preconditioned Conjugate Gradient Method （CG）
Compute r(0)= b-[A]x(0)for i= 1, 2, …
solve [M]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)r(i)= r(i-1) - αiq(i)check convergence |r|end

前処理：対角スケーリング

pFEM3D-2

77

対角スケーリング，点ヤコビ前処理• 前処理行列として，もとの行列の対角成分のみを取り出した行列を前処理行列 [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)という場合に逆行列を簡単に求めることができる。

pFEM3D-2

78

CG法（1/6）
subroutine CG(N, NP, NPLU, D, AMAT, index, item, B, X, RESID,
ITER, ERROR, my_rank,
NEIBPETOT, NEIBPE, IMPORT_INDEX, IMPORT_ITEM,EXPORT_INDEX, EXPORT_ITEM)

&
&
&
&

use

&
&
&
&

solver_SR

implicit REAL*8(A-H,O-Z)
include 'precision.inc'include 'mpif.h'
integer(kind=kint ), intent(in):: N, NP, NPLU, my_rank
integer(kind=kint ), intent(in):: NEIBPETOT
integer(kind=kint ), intent(inout):: ITER, ERROR
real (kind=kreal), intent(inout):: RESID
real(kind=kreal), dimension(NP) , intent(inout):: B, X, Dreal(kind=kreal), dimension(NPLU), intent(inout):: AMAT
integer(kind=kint ), dimension(0:NP),intent(in) :: index
integer(kind=kint ), dimension(NPLU),intent(in) :: item
integer(kind=kint ), pointer :: NEIBPE(:)
integer(kind=kint ), pointer :: IMPORT_INDEX(:), IMPORT_ITEM(:)
integer(kind=kint ), pointer :: EXPORT_INDEX(:), EXPORT_ITEM(:)
real(kind=kreal), dimension(:),
real(kind=kreal), dimension(:,:),

allocatable:: WS, WR
allocatable:: WWinteger(kind=kint), parameter :: R= 1
integer(kind=kint), parameter :: Z= 2
integer(kind=kint), parameter :: Q= 2
integer(kind=kint), parameter :: P= 3
integer(kind=kint), parameter :: DD= 4
integer(kind=kint ) :: MAXIT
real (kind=kreal) :: TOL, W, SS

送信バッファ，受信バッファ

pFEM3D-2

79

CG法（2/6）
COMMtime= 0.d0
COMPtime= 0.d0
ERROR= 0
allocate (WW(NP,4),WR(NP),WS(NP))
MAXIT
TOL

= ITER= RESID

X = 0.d0
WS= 0.d0
WR= 0.d0!C!C +-----------------------+!C | {r0}= {b} - [A]{xini} |!C +-----------------------+!C===
call SOLVER_SEND_RECV( NP, NEIBPETOT, NEIBPE, IMPORT_INDEX, IMPORT_ITEM,EXPORT_INDEX, EXPORT_ITEM, WS, WR, X , my_rank)&&

do j= 1, NWW(j,DD)= 1.d0/D(j)WVAL= B(j) - D(j)*X(j)do k= index(j-1)+1, index(j)i= item(k)WVAL= WVAL - AMAT(k)*X(i)enddo
WW(j,R)= WVAL
enddo

Compute r(0)= b-[A]x(0)for i= 1, 2, …
solve [M]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)
endif
&
&
q(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

BNRM20= 0.d0
do i= 1, N
BNRM20= BNRM20 + B(i)**2
enddo
call MPI_Allreduce (BNRM20, BNRM2, 1, MPI_DOUBLE_PRECISION,&MPI_SUM, MPI_COMM_WORLD, ierr)&

pFEM3D-2

80

SOLVER_SEND_RECV（1/2）
subroutine&&&

SOLVER_SEND_RECV&( N, NEIBPETOT, NEIBPE, IMPORT_INDEX, IMPORT_ITEM,&EXPORT_INDEX, EXPORT_ITEM,&
WS, WR, X, my_rank)

implicit REAL*8 (A-H,O-Z)
include 'mpif.h'include 'precision.inc'integer(kind=kint ), intent(in) :: N
integer(kind=kint ), intent(in) :: NEIBPETOT
integer(kind=kint ), pointer :: NEIBPE(:)
integer(kind=kint ), pointer :: IMPORT_INDEX(:)
integer(kind=kint ), pointer :: IMPORT_ITEM (:)
integer(kind=kint ), pointer :: EXPORT_INDEX(:)
integer(kind=kint ), pointer :: EXPORT_ITEM (:)
real(kind=kreal), dimension(N), intent(inout):: WS
real(kind=kreal), dimension(N), intent(inout):: WR
real(kind=kreal), dimension(N), intent(inout):: X
integer, intent(in) :: my_rank
integer(kind=kint ), dimension(:,:), save, allocatable :: sta1, sta2, req1, req2
integer(kind=kint ), save :: NFLAG
data NFLAG/0/
if (NFLAG.eq.0) thenallocate (sta1(MPI_STATUS_SIZE,NEIBPETOT), sta2(MPI_STATUS_SIZE,NEIBPETOT))
allocate (req1(NEIBPETOT), req2(NEIBPETOT))NFLAG= 1
endif
do neib= 1, NEIBPETOT
istart= EXPORT_INDEX(neib-1)
inum = EXPORT_INDEX(neib ) - istart
do k= istart+1, istart+inum
ii = EXPORT_ITEM(k)WS(k)= X(ii)enddo
call MPI_Isend (WS(istart+1), inum, MPI_DOUBLE_PRECISION,
NEIBPE(neib), 0, MPI_COMM_WORLD, req1(neib),
ierr)
enddo

&
&

&
&

pFEM3D-2

81

SOLVER_SEND_RECV（2/2）
do neib= 1, NEIBPETOT
istart= IMPORT_INDEX(neib-1)
inum = IMPORT_INDEX(neib ) - istartcall MPI_Irecv (WR(istart+1), inum, MPI_DOUBLE_PRECISION,&
NEIBPE(neib), 0, MPI_COMM_WORLD, req2(neib),
&
ierr)
enddo
call MPI_Waitall (NEIBPETOT, req2, sta2, ierr)do neib= 1, NEIBPETOT
istart= IMPORT_INDEX(neib-1)
inum = IMPORT_INDEX(neib ) - istart
do k= istart+1, istart+inum
ii = IMPORT_ITEM(k)X(ii)= WR(k)
enddo
enddo
call MPI_Waitall (NEIBPETOT, req1, sta1, ierr)
end subroutine solver_send_recv
end module
solver_SR&&

pFEM3D-2

82

CG法（3/6）
do iter= 1, MAXIT

Compute r(0)= b-[A]x(0)for i= 1, 2, …
solve [M]z(i-1)= r(i-1)ρi-1= r(i-1) z(i-1)
if i=1!C!C +---------------+
p(1)= z(0)!C | {RHO}= {r}{z} |!C +---------------+else
RHO0= 0.d0
βi-1= ρi-1/ρi-2
do i= 1, N
RHO0= RHO0 + WW(i,R)*WW(i,Z)p(i)= z(i-1) + βi-1 p(i-1)
enddo
endif
call MPI_Allreduce (RHO0, RHO, 1, MPI_DOUBLE_PRECISION,
q(i)= [A]p(i)&MPI_SUM, MPI_COMM_WORLD, ierr)
αi = ρi-1/p(i)q(i)!C!C +-----------------------------+
x(i)= x(i-1) + αip(i)!C | {p} = {z} if
ITER=1|!C | BETA= RHO / RHO1 otherwise |r(i)= r(i-1) - αiq(i)!C +-----------------------------+
if ( ITER.eq.1 ) thencheck convergence |r|do i= 1, N
WW(i,P)= WW(i,Z)
end!C!C +----------------+!C | {z}= [Minv]{r} |!C +----------------+
do i= 1, N
WW(i,Z)= WW(i,R) * WW(i,DD)
enddo

enddo
else
BETA= RHO / RHO1
do i= 1, N
WW(i,P)= WW(i,Z) + BETA*WW(i,P)
enddo
endif

pFEM3D-2

83

CG法（4/6）!C!C +-------------+!C | {q}= [A]{p} |!C +-------------+
call SOLVER_SEND_RECV&( NP, NEIBPETOT, NEIBPE, IMPORT_INDEX, IMPORT_ITEM,&EXPORT_INDEX, EXPORT_ITEM, WS, WR, WW(1,P),&
my_rank)
do j= 1, N
WVAL= D(j)*WW(j,P)
do k= index(j-1)+1, index(j)i= item(k)X1= WW(3*i-2,P)
X2= WW(3*i-1,P)
X3= WW(3*i ,P)
WVAL= WVAL + AMAT(k)*WW(i,P)
enddo
WW(j,Q)= WVAL
enddo
!C
!C +---------------------+!C | ALPHA= RHO / {p}{q} |!C +---------------------+
C10= 0.d0
do i= 1, N
C10= C10 + WW(i,P)*WW(i,Q)
enddo
call MPI_Allreduce (C10, C1, 1, MPI_DOUBLE_PRECISION,&MPI_SUM, MPI_COMM_WORLD, ierr)
ALPHA= RHO / C1

Compute r(0)= b-[A]x(0)for i=& 1, 2, …
&
solve[M]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)r(i)= r(i-1) - αiq(i)check convergence |r|end

pFEM3D-2

84

CG法（5/6）!C!C +----------------------+!C | {x}= {x} + ALPHA*{p} |!C | {r}= {r} - ALPHA*{q} |!C +----------------------+
do i= 1, N
X(i) = X (i) + ALPHA * WW(i,P)
WW(i,R)= WW(i,R) - ALPHA * WW(i,Q)
enddo
DNRM20= 0.d0
do i= 1, N
DNRM20= DNRM20 + WW(i,R)**2
enddo
call MPI_Allreduce (DNRM20, DNRM2, 1,&
MPI_DOUBLE_PRECISION,&MPI_SUM, MPI_COMM_WORLD, ierr)
RESID= dsqrt(DNRM2/BNRM2)
if ( RESID.le.TOL ) exit
if ( ITER .eq.MAXIT ) ERROR= -300
RHO1 = RHO
enddo

!C===

30 continue
call SOLVER_SEND_RECV( NP, NEIBPETOT, NEIBPE, IMPORT_INDEX, IMPORT_ITEM,EXPORT_INDEX, EXPORT_ITEM, WS, WR, X , my_rank)&&

deallocate (WW,WR,WS)
end subroutineCG
end module
solver_CG

Compute r(0)= b-[A]x(0)for i= 1, 2, …
solve [M]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)r(i)= r(i-1) - αiq(i)checkconvergence |r|&end &pFEM3D-2

85

CG法（6/6）!C!C +----------------------+!C | {x}= {x} + ALPHA*{p} |!C | {r}= {r} - ALPHA*{q} |!C +----------------------+
do i= 1, N
X(i) = X (i)+ ALPHA * WW(i,P)
WW(i,R)= WW(i,R) - ALPHA * WW(i,Q)
enddo
DNRM20= 0.d0
do i= 1, N
DNRM20= DNRM20 + WW(i,R)**2
enddo
call MPI_Allreduce (DNRM20, DNRM2, 1,&
MPI_DOUBLE_PRECISION,&MPI_SUM, MPI_COMM_WORLD, ierr)
RESID= dsqrt(DNRM2/BNRM2)
if ( RESID.le.TOL) exitif ( ITER .eq.MAXIT ) ERROR= -300
RHO1 = RHO
!C===

enddo

30 continue
call SOLVER_SEND_RECV( NP, NEIBPETOT, NEIBPE, IMPORT_INDEX, IMPORT_ITEM,EXPORT_INDEX, EXPORT_ITEM, WS, WR, X , my_rank)&&

deallocate (WW,WR,WS)
end subroutineCG
end module
solver_CG&&

外点に最新のX（温度）代入

pFEM3D-2

86

OUTPUT_UCD• 各領域からintELEM_listに所属する要素の情報を集める– intELEM_list： 各節点の属する領域番号のうち最も若い番号の領域に所属する・・・と見なす• MPI_Allgathervを使って以下の情報を一箇所に集める– 節点：節点座標，温度– 要素：要素コネクティビティ（要素を構成する節点）• 節点の情報は一部重複• 問題規模が大きくなると困難– あまり賢いやり方ではない– 並列可視化pFEM3D-2

87

Example: Strong Scaling: Fortran• 128×128×128 nodes, 2,097,152 DOF• 32～896 cores, Flat MPI• Linear Solver
Memory
96 GB

DDR4
DDR4
DDR4
DDR4
DDR4
DDR4

2933 MHz×6ch
140.8 GB/sec

Intel® Xeon®
Platinum 8280

(Cascade Lake, CLX)
2.7GHz, 28-Cores
2.419 TFLOPS
Soc. #0: 0th-27th cores

UPI
UPI
UPI

Intel® Xeon®
Platinum 8280
(Cascade Lake, CLX)
2.7GHz, 28-Cores
2.419 TFLOPS
Soc. #1: 28th-55th cores

Ultra Path Interconnect
10.4 GT/sec ×3= 124.8 GB/sec

DDR4
DDR4
DDR4
DDR4
DDR4
DDR4

Memory96 GB2933 MHz×6ch
140.8 GB/sec

900
800
700
600

Speed-Up◆ 1x16 proc’s for each Socket
1x16x2 (=32) for each Node
○ 1x20x2 (=40), pMETIS● 1x24x2 (=48), pMETIS● 1x28x2 (=56), pMETIS△ 1x20x2 (=40), kMETIS▲ 1x24x2 (=48), kMETIS▲ 1x28x2 (=56), kMETIS

1x16
1x20-pMETIS
1x24-pMETIS
1x28-pMETIS
1x20-kMETIS
1x24-kMETIS
1x28-kMETIS

500
400
300
200

Performance of ◆ without
NUMA at 32-cores= 32.0

100
0
0

100

200

300

400

500

Core #

600

700

800

900

pFEM3D-2

88

Example: Strong Scaling: Fortran• Best Case at 16 nodes (with/without NUMA Ctrl)– Difference is not so large• 48-cores/node●▲ are the best
• No significant difference between kMeTis & pMeTis
900
800
700
600

Speed-Up◆ 1x16 proc’s for each Socket
1x16x2 (=32) for each Node
○ 1x20x2 (=40), pMETIS● 1x24x2 (=48), pMETIS● 1x28x2 (=56), pMETIS△ 1x20x2 (=40), kMETIS▲ 1x24x2 (=48), kMETIS▲ 1x28x2 (=56), kMETIS

1x16
1x20-pMETIS
1x24-pMETIS
1x28-pMETIS
1x20-kMETIS
1x24-kMETIS
1x28-kMETIS

500
400
300
200

Performance of ◆ without
NUMA at 32-cores= 32.0

100
0
0

100

200

300

400

500

Core #

600

700

800

900

pFEM3D-2

89

Example: Strong Scaling: Fortran• 1x16: Parallel Mesh Generator
128 128 128
4
4
2
pcube

128 128 128
4
4
4
pcube

128 128 128
8
4
4
pcube

128 128 128
8
8
4
pcube

node=
1
procs= 32

node=2
procs= 64

node=
4
procs=128

node=8
procs=256

900
800
700
600

Speed-Up

Node
#

1x16
1x20-pMETIS
1x24-pMETIS
1x28-pMETIS
1x20-kMETIS
1x24-kMETIS
1x28-kMETIS

500
400
300
200
100
0
0

100

200

300

400

500

Core #

600

700

800

900

p20
k20

p24
k24

p28
k28

1

40

48

56

2

80

96

112

4

160

192

224

8

320

384

448

pFEM3D-2

90

課題（1/2）• 自分で問題設定を行い，「sol」の挙動を分析してみよ• 例– Strong Scaling• 問題サイズを固定，PE数を変化させて時間（全体，各部分）を測定。

– Weak Scaling• PEあたりの問題サイズを固定，1反復あたりの計算時間を求める。

– 考慮すべき項目• 問題サイズ• 領域分割手法（RCB，K-METIS，P-METIS，1D～3D）の影響。

– メッシュ生成，領域分割におけるFX10の性能低い• 1283くらいが限界（領域分割に15分以上かかる）• 「*.inp」の出力に時間がかかる場合がある。
– OUTPUT_UCDの呼び出しのコメントアウト– src，partpFEM3D-2

91

１D～3D分割
どの方法が良いか考えて見よ

1D型

2D型

3D型

pFEM3D-2

92

１D～3D分割
通信量の総和（各辺4N，8領域とする）

1D型

2D型

3D型

16 N2 x 7 = 112 N2

16 N2 x 4 = 64 N2

16 N2 x 3 = 48 N2

pFEM3D-2

93

１D～3D分割
mesh.inp

1D型

2D型

3D型

64 64 64
8 1 1
pcube

64 64 64
4 2 1
pcube

64 64 64
2 2 2
pcube

pFEM3D-2

94

課題（2/2）• 領域間通信（solver_SR）の性能改善ができないかどうか考えてみよ。
– Recv. Bufferへのコピーを効率的に実施できないか?
do neib= 1, NEIBPETOT
istart= IMPORT_INDEX(neib-1)
inum = IMPORT_INDEX(neib ) - istartcall MPI_Irecv (WR(istart+1), inum, MPI_DOUBLE_PRECISION,&
NEIBPE(neib), 0, MPI_COMM_WORLD, req2(neib),
&
ierr)
enddo
call MPI_Waitall (NEIBPETOT, req2, sta2, ierr)do neib= 1, NEIBPETOT
istart= IMPORT_INDEX(neib-1)
inum = IMPORT_INDEX(neib ) - istart
do k= istart+1, istart+inum
ii = IMPORT_ITEM(k)X(ii)= WR(k)
enddo
enddo

&
&

95

pFEM3D-1

SEND/RECV (Original)!C!C-- INIT.
allocate (sta1(MPI_STATUS_SIZE,NEIBPETOT), sta2(MPI_STATUS_SIZE,NEIBPETOT))
allocate (req1(NEIBPETOT), req2(NEIBPETOT))!C!C-- SEND
do neib= 1, NEIBPETOT
istart= STACK_EXPORT(neib-1)
inum = STACK_EXPORT(neib ) - istart
do k= istart+1, istart+inum
WS(k)= X(NOD_EXPORT(k))
enddocall MPI_ISEND (WS(istart+1), inum, MPI_DOUBLE_PRECISION,&&
NEIBPE(neib), 0, MPI_COMM_WORLD, req1(neib), ierr)enddo!C!C-- RECEIVE
do neib= 1, NEIBPETOT
istart= STACK_IMPORT(neib-1)
inum = STACK_IMPORT(neib ) - istartcall MPI_IRECV (WR(istart+1), inum, MPI_DOUBLE_PRECISION,&&
NEIBPE(neib), 0, MPI_COMM_WORLD, req2(neib), ierr)enddo
call MPI_WAITALL (NEIBPETOT, req2, sta2, ierr)do neib= 1, NEIBPETOT
istart= STACK_IMPORT(neib-1)
inum = STACK_IMPORT(neib ) - istart
do k= istart+1, istart+inum
X(NOD_IMPORT(k))= WR(k)
enddo
enddo
call MPI_WAITALL (NEIBPETOT, req1, sta1, ierr)

96

pFEM3D-1

If numbering of external nodes is
continuous in each neighboringprocess ...
84

81

85

82

83

86

88

87

96

57

58

59

60

61

62

63

64

73

95

49

50

51

52

53

54

55

56

74

94

41

42

43

44

45

46

47

48

80

93

33

34

35

36

37

38

39

40

79

92

25

26

27

28

29

30

31

32

78

91

17

18

19

20

21

22

23

24

77

90

9

10

11

12

13

14

15

16

76

89

1

2

3

4

5

6

7

8

75

65

66

67

68

69

70

71

72

Parallel FEM 3D-2
pFEM3D-1

97

SEND/RECV (NEW:1)!C!C-- INIT.
allocate (sta1(MPI_STATUS_SIZE,2*NEIBPETOT))
allocate (req1(2*NEIBPETOT))!C!C-- SEND
do neib= 1, NEIBPETOT
istart= STACK_EXPORT(neib-1)
inum = STACK_EXPORT(neib ) - istart
do k= istart+1, istart+inum
WS(k)= X(NOD_EXPORT(k))
enddo
enddo
do neib= 1, NEIBPETOT
istart= STACK_EXPORT(neib-1)
inum = STACK_EXPORT(neib ) - istartcall MPI_ISEND (WS(istart+1), inum, MPI_DOUBLE_PRECISION,&
NEIBPE(neib), 0, MPI_COMM_WORLD, req1(neib), ierr)
enddo

&

!C
!C-- RECEIVE
do neib= 1, NEIBPETOT
inum = STACK_IMPORT(neib) - STACK_IMPORT(neib-1)
istart= NOD_IMPORT(STACK_IMPORT(neib-1)+1)&

call MPI_IRECV (X(istart), inum, MPI_DOUBLE_PRECISION,&
NEIBPE(neib), 0, MPI_COMM_WORLD, req1(NEIBPETOT+neib), ierr)enddo
call MPI_WAITALL (2*NEIBPETOT, req1, sta1, ierr)

pFEM3D-1

98

SEND/RECV (NEW:2), N0: int. node #!C!C-- INIT.
allocate (sta1(MPI_STATUS_SIZE,2*NEIBPETOT))
allocate (req1(2*NEIBPETOT))!C!C-- SEND
do neib= 1, NEIBPETOT
istart= STACK_EXPORT(neib-1)
inum = STACK_EXPORT(neib ) - istart
do k= istart+1, istart+inum
WS(k)= X(NOD_EXPORT(k))
enddo
enddo

N0: 内点総数

do neib= 1, NEIBPETOT
istart= STACK_EXPORT(neib-1)
inum = STACK_EXPORT(neib ) - istartcall MPI_ISEND (WS(istart+1), inum, MPI_DOUBLE_PRECISION,&
NEIBPE(neib), 0, MPI_COMM_WORLD, req1(neib), ierr)
enddo

&

!C
!C-- RECEIVE
do neib= 1, NEIBPETOT
inum = STACK_IMPORT(neib) - STACK_IMPORT(neib-1)
istart= STACK_IMPORT(neib-1) + N0 + 1&

call MPI_IRECV (X(istart), inum, MPI_DOUBLE_PRECISION,&
NEIBPE(neib), 0, MPI_COMM_WORLD, req1(NEIBPETOT+neib), ierr)enddo
call MPI_WAITALL (2*NEIBPETOT, req1, sta1, ierr)