PowerPoint プレゼンテーション

４. Numerical Simulation of Seismic Waves
2-1. Equation of Motions for Elastic Waves
2-2. FDM Simulation and Numerical Error
2-3. Boundary Conditions
2-4. Large-Scale 3D Simulation by Parallel Computers
強震動シミュレーション手法
１．有限要素法（FEM）
・FEMメッシュ
・地震工学で人気（入力モデルが確定）
・Spectral Element法への拡張（1990年後半）
２．差分法（FDM）
・規則メッシュ
・地震学で人気（入力モデルの調整・検証）
３．境界要素法（BEM）
・1980年代に流行
・高速多重極展開の導入（1990年後半）
差分法による運動方程式の計算
運動方程式（釣り合いの式）
物性定数（ラメの定数、密度）
Body Force



ux  xx  xy  xz  f x ,
x
y
z
Acceleration



uy  yz  yy  yz  f y ,
x
y
z



uz  xz  yz  zz  f z .
x
y
z
Density
Stress
応力ー歪構成方程式（フックの法則）
 pq  (exx  eyy  ezz ) pq  2 epq
歪みテンソル
u y
u
u
exx  x , e yy 
, ezz  z ,
x
y
z
  ,  ,  ( x, y, z )
 u x , u y , u z ( x, y, z )
差分法による微分計算（２次精度中央差分）
u x ( x, y, z ) u x ( x  Dx, y, z )  u x ( x  Dx, y, z )

x
2Dx
Ux
Strain
1  u uy 
1  uy u 
exy   x  , e yz    z ,
2  y x 
2  z y 
1u u 
exz   x  z .
2  z x 
2Dx
x-Dx
x+Dx
x
FDM演算子の導出：テイラー展開
2次精度FDM
Dx 
d
f ( x  Dx)  f ( x)  Dx f ( x) 
dx
2
2

d
Dx 
f ( x  Dx)  f ( x)  Dx f ( x) 
dx
2
2
d2
f ( x)  O[(Dx) 2 ] (1)
dx2
d2
f ( x)  O[(Dx) 2 ] (2)
2
dx
+1/2
-1/2
From (1)-(2)
d
f ( x  Dx)  f ( x  Dx)
f ( x) 
 O[( Dx) 2 ]
dx
2Dx
f(x- Dx) f(x)
f(x+ Dx)
4次精度 FDM


d
Dx  d 2
Dx  d 3
f ( x  Dx)  f ( x)  Dx f ( x) 
f ( x) 
f ( x)  O[(Dx) 4 ]
2
3
dx
2 dx
6 dx
2
3
2


d
Dx  d
Dx  d 3
f ( x  Dx)  f ( x)  Dx f ( x) 
f ( x) 
f ( x)  O[(Dx) 4 ]
2
3
dx
2 dx
6 dx
2
3


d
2Dx  d 2
2Dx  d 3
f ( x  2Dx)  f ( x)  2Dx f ( x) 
f ( x) 
f ( x)  O[(Dx) 4 ]
dx
2 dx2
6 dx3
2Dx 2 d 2 f ( x)  2Dx 3 d 3 f ( x)  O[(Dx) 4 ]
d
f ( x  2Dx)  f ( x)  2Dx f ( x) 
dx
2 dx2
6 dx3
2
3
(1)
(2)
-1/12
8/12
(3)
-8/12
(4)
+1/12
f(x-2 Dx) f(x- Dx)
From (1)(2)(3)(4)
d
 f ( x  2Dx)  8 f ( x  Dx)  8 f ( x  Dx)  f ( x  2Dx)
f ( x) 
 O[( Dx) 4 ]
dx
12 Dx
f(x)
f(x+ Dx) f(x+2 D x)
FDMの次数と精度（２）：数値実験
(a) 4th-order FDM
Example of
2D modeling (P-SV)
PPPP
PPP
PS
Time [s]
Time [s]
4 Grid Points/Wavelength
(b) 8th-order FDM
PP
P
x
z
(d) 128th-order FDM
Time [s]
(c) 16th-order FDM
Time [s]
z
x
FDMの次数と精度（１）：実効波数
周期関数
f ( x)  exp(ikx)
k : 波数
その微分（解析解）は
df ( x)
 ik exp( ikx )  ik f ( x)
dx
解析解 (kDx)
２次中央差分による微分
df ( x)
dx
exp[ik ( x  Dx)]  exp[ik ( x  Dx)]

2Dx
sin(kDx)
i
f ( x)
Dx
 ike f ( x)
ke : 実効波数
誤差
数値解 [sin(kDx)]
FDMの高精度化：Staggered Grid （食い違い格子）の利用
Staggered-grid FDM, Virieux (1986)、Geophysics, V108, 889-901.
(b) Staggered Grid Model
(a) Ordinary Grid Model
変位、応力、物性値すべてが同一格子点に定義
Dx
Dx
Dx
変位、応力、物性値は半格子ずれた位置に定義
Dx
Dx
Dx
Dz
Dx/2 Dx/2
Dx
Dx
Dz
＊格子サイズが見かけ上、半分になる
The Staggered-Grid FDM
(a) Ordinary FDM Grid
(b) Staggered-Grid FDM
Dx
Dx/2
d
f ( x  Dx)  f ( x  Dx)
f ( x) 
dx
2Dx
d
f ( x  1 / 2Dx)  f ( x  1 / 2Dx)
f ( x) 
dx
Dx
Exercise 6. Assume a new grid configuration where the stress () and displacement (v) components are defined halfway between
each others, and the derivatives are defined on the middle. This is called the staggered grid. Currently the staggered-grid FDM
has been used practically in the wave propagation simulation. Derive the effective wavenumber for the 2nd-order central
staggered-grid FDM and compare the accuracy with that for the ordinary 2nd-order central FDM.
スタガード格子FDMの精度：実効波数の確認
・周期関数
f ( x)  exp(ikx)
k : 波数
・その微分（解析解）
df ( x)
 ik exp( ikx )  ik f ( x)
dx
解析解 (kDx)
・スタガード格子における２次FDM微分の解
Dx
Dx
)]  exp[ik ( x  )]
df ( x)
2
2

dx
Dx
Dx
sin(k )
2 f ( x)
i
Dx
2
 ike f ( x)
exp[ik ( x 
ke : 実効波数
誤差
数値解 [sin(kDx/2)]
フーリエスペクトル法 (PSM): FFTを利用した高次FDM
(1) フーリエ変換
N 1
F (lDk )  Dx f (nDx) exp(i 2nl / N )
n 0
1 N 1
f (nDx) 
 F (lDk ) exp(i2nl / N )
NDx l 0
1.FFTをかける
(2) 波数空間での解析的な微分、フーリエ逆変換
d
f (nDx)
dx
1 N 1

 i(lDx) F (lDk ) exp(i 2nl / N )
NDx l 0
2.波数を乗じる
3.FFTでもとに戻す
FFTの利用
- 高速化
乗算の計算量
N2 
 N log2 N
- 周期境界条件 (Wraparound)
FDM計算はベクトル化率が高い
スカラー計算
PSM～16th-order FDM
－FDM計算速度が数十倍に加速する
－PSM(FFT)計算はベクトル化が難しい
ベクトル計算
PSM << 128th-orderFDM
加速
自由表面の境界条件
(a) 真空条件
(b) 対照法
自由表面 (Z=0)
におけるゼロ応力
 xz   yz   zz  0
＊自由表面は難しい。。。。
解析解
FDM
高次FDM演算子：端の処理
(b) 周期境界にする
(a) 次数を低下させる
6
4
精度が低下
Wrapround
が起きる
反射波が出る
吸収境界条件 (1): Dampingの利用
－経験的（物理的背景うすい）
－簡単。でもよく効く
吸収境界条件
Cerjan et al. (1985), Geophysics,
V50, 705-808.
S
W (i)  exp[0.015(20  i) 2 ],
(i  1.2...,20)
sP
P
夢の吸収境界条件：Perfect Matched Layer (PML)
PML
鉛直
・境界面に鉛直入射する波は透過
・境界面に平行に入射する波を減衰させる
夢のPML
平行
境界のインピーダンスコントラストが０となる＝無反射
無反射境界なし
＊現実はそんなに甘くない、永遠の課題
FDMの陽解法：時間発展
Dt  0.26
1. 変位－応力の式
Dx
V max
n
n
n

1   xp  yp  zp
u 


 f pn 

  x
y
z

速度
u np1/ 2  u np1/ 2  unp Dt 時間積分１
加速度
n
p
変位
n 1
p
u
応力
n
 pq
 u np  u np1/ 2 Dt
u np1/ 2  u np1/ 2  unp1/ 2 Dt

n 1
pq
 
n
pq
 
n 1 / 2
pq
Dt
時間積分２
変位
 u xn u yn u zn 
 u np uqn 
 pq   

 



 x y

 q p 

z




n-1
n
n-1/2
n-1
n+1/2
2. 速度－応力式
速度
u np1/ 2
応力 n

1   xp  yq  zq
 u np1/ 2  


 f pn  Dt

  x
y
z

n
n
時間積分１
速度
  u n1 / 2 u n1 / 2 u n1 / 2
y
n 1
n
 pq
  pq
   x 
 z
y
z
  x
応力
n1 / 2
n1 / 2

 u p
 q 

u
   
 Dt

 pq
 q
p 



時間積分２
＊速度－応力式のほうが
計算が安定（らしい）
非弾性減衰（Ｑ）の導入
１．Damperの利用
Graves (1996), BSSA, V86, 1091-1106.

 u y u 
 u x u y u z 


  2 
 xx  (  2 )


 z   A
 y
x
y
z 
z 





  f 0t 
A  exp 

 Q 
- 強い周波数依存性Q( f )  Q0
＊ほとんどの強震動計算はこれ
が使われている（かんたん）
f
f0
- QpとQsを分離できない（Q=Qp=Qs）
２．Memory Variableの利用
・Robertsson et al., (1994), Geophysics, V59, 1444-1456.
・Hestholm (1999), Geophys. J. Int. V139, 852-878.
 P  u x u y u z 
 S
 xx  (  2 ) 


 2
   x
y
z 

 u y u z 

  rxx

 y


z


  P  u x u y u z 
  S  u y u z 
1 



rxx   rxx  (  2 )  1


 2   1




  
y
z 
z 
 
 x
 
 y
1   QS
1 
1
1  P
1
 
1 2 
,    2 ,  S 
 
Q p Q p 
 
QS   2 
- Constant Q,Q( f )  Q0
- Qp, Qs分離できる
地震動シミュレーション入力・出力
fp  

q x, y , z
M pq

 ( x, y, z )
q
入力データ:
up 
 xp
x

 yp
y

 zp
z
 fp
 pq  (exx  eyy  ezz ) pq  2 epq
z
・各点の物性値を入力
・層境界面を入力、多層構造でモデル
化
x
断層型震源（Mxz=Mzx=M0）
z
出力データ：
○３次元各点の揺れ（速度） V=（Vx, Vy, Vz）
○地表面の揺れ（２次元）
・P波、S波分離： P=div V, S=rot V
・加速度（時間微分）、変位（積分）
x
爆発型震源（Mxx=Mzz=M0)
応力
１次元波動方程式のFDM計算
 x (1)
 x
(1)
x
速度 Vx (1) Vx (1)
スタガード格子に、以下の変数を配置
x
・応力とその空間微分
・速度＊とその空間微分
c-Solving equation of motions
・物性値（剛性率、密度）
do it = 1, ntmax
 x (2)
Vx
(2)
x
 x
( 2)
x
Vx (2)
time = it*dt
if (mod(it-1,nsnap).eq.0) then
call SNAP(Vx, nx, isnap, time, dx)
isnap = isnap + 1
end if
＊各格子点の揺れの速度、
媒質の地震波伝播速度のことではない
do i = 2, nx-2
dxSx(i) = (0.00417*Sx(i-1)-1.25*Sx(i)
:
+ 1.25*Sx(i+1)0.00417*Sx(i+2))/dx
end do
dxSx(1)=0; dxSx(nx-1)=0; dxSx(nx)=0
do i = 1, nx
Vx (i) = Vx (i) + (dxSx(i))/den(i)*dt
end do
初期化
:
parameter (nx=150, dx=1.0, dt=0.125)
parameter (ntmax=2000, nsnap=20)
real Vx(nx), Sx(nx), den(nx), rig(nx)
real dxVx(nx), dxSx(nx)
character fname*17
c-Initial Condition
call MODEL (den,rig,nx)
call INIT (Vx,Sx,nx)
call SOURCE (Vx,nx)
isnap = 0
do i = 3, nx-1
dxVx(i) = (0.00417*Vx(i-2)-1.25*Vx(i-1)
+ 1.25*Vx(i)-0.00417*Vx(i+1))/dx
end do
dxVx(1)=0; dxVx(2)=0; dxVx(nx)=0
do i = 1, nx
Sx (i) = Sx (i) + rig(i)*dxVx(i)*dt
end do
end do
 x (3)
 x (4)
Vx
(3)
x
Vx (4)
Vx (3)
Dx
応力の空間微分（４次精度FDM）
 x (i )

x
c x (i  1)  d x (i )  e x (i  1)  f x (i  2)
Dx
速度値の計算
vx
n 1/ 2
1  x
 vx 
Dt
 x
n
n
速度の空間微分（４次精度FDM）
v x (i )

x
cvx (i  2)  dvx (i  1)  evx (i )  fvx (i  1)
Dx
応力値の計算
x
n 1
x
n 1/ 2
c-write (6,*) 'Last frame number= ', isnap
end
以下、繰り返し
n 1/ 2
v
 x
x
Dt
１次元波動方程式のFDM計算
モデル格子に、物性値（Vs, )を置く
subroutine MODEL (den,rig,nx)
subroutine SOURCE(Vx,nx)
real den(nx),rig(nx)
real Vx(nx)
震源（揺れの速度の初期分布）
do i = 1, nx
ix = 8
if (i.lt.25) then
do i = 1, nx
vs = 0.0
if (i.lt.(nx/2-ix)) then
ro = 2.0
Vx (i) = 0.0
else if (i.lt.75) then
else if (i.lt.(nx/2+ix)) then
vs = 0.5
Vx (i) = sin(3.14*(i-nx/2)/real(ix))
ro = 2.0
else
else if (i.lt.125) then
Vx (i) = 0.0
vs = 0.25
end if
ro = 1.5
end do
else
return
vs = 0.0
end
揺れの分布の出力
ro = 1.5
end if
subroutine SNAP(Vx, nx, isnap, time, dx)
den (i) = ro
real Vx (nx)
rig (i) = ro*vs**2
character fname*17
end do
write (fname(1:17),'(A13,I3.3)') 'results/WAVE.', isnap
return
open (16,file=fname,
end
status='unknown',form='formatted')
初期化
write (16,*) time
subroutine INIT (Vx,Sx,nx)
do i = 1, nx
real Vx(nx), Sx(nx)
write (16,*) i*dx, Vx(i)
do i = 1, nx
end do
Vx (i) = 0.0
close (16)
Sx (i) = 0.0
return
end do
end
return
end
#!/bin/csh -f
お絵かきスクリプト
gmtset TICK_LENGTH -0.1
set RANGE = '1.0/150.0/-1.5/1.5'
set SCALE = '6/2'
set INFILE = 'results/WAVE'
set OUTFILE = 'anime/Movie'
set PSFILE = tmp.ps
@ istart = 0
@ iend = 100
@ i = $istart
スナップショットの時刻と
while ($i <= $iend)
速度分布を取り出す
if ( $i < 10 ) then
head -1 ${INFILE}.00$i >tmp_t
tail +2 ${INFILE}.00$i >tmp_w
else if ( $i < 100 ) then
head -1 ${INFILE}.0$i >tmp_t
tail +2 ${INFILE}.0$i >tmp_w
else
head -1 ${INFILE}.$i >tmp_t
tail +2 ${INFILE}.$i >tmp_w
endif
#-Time
awk '{printf "120 0.5 30 0 0 0 %3.1fs\n",$1}' < tmp_t > tmp_tt
psxy tmp_w -R${RANGE} -JX${SCALE} \
-W4/0/0/100 -P -K >! ${PSFILE}
GMTで絵を描く
pstext tmp_tt -R${RANGE} -JX${SCALE} \
-Ba50f10g25:"Distance [km]":/a0.5f0.125g1:"Height
[m]":WSne /
-O >> ${PSFILE}
eps2eps ${PSFILE} tmp.eps PSをbmp連番ファイルに変
換
if ( $i < 10 ) then
convert tmp.eps ${OUTFILE}.00$i.bmp
else if ( $i < 100 ) then
convert tmp.eps ${OUTFILE}.0$i.bmp
else
convert tmp.eps ${OUTFILE}.$i.bmp
endif
echo '[ ' $i ' ]'
@ i++
Top500(http://www.top500.org)＊LINPAC（連立一次方程式の解）ベンチマークテスト
＊演算性能（FLOPS: 1秒間あたりの小数点演算回数）
＊IBM BlueGene/Lは212,992CPUの並列計算で478TFLOPS
IBM Blue Gene/L
メキシコ自治大学
溶岩公園
地球シミュレータ
＊地球シミュレータは30位に落ちた
＊PCを1345台をLANで繋げば、TOP５００に入る
中国兵馬俑
Top500(http://www.top500.org)
・スパコンは10年で1000倍高速化
・2012年には京速（10PFLOPS)計算機が完成（理化学研究所、神戸）
京速計算機
1987
Lets Note
・私のノートPCは1987年（M2）にはTop 1であった
ベクトル計算機とスカラー計算機
・CPUクロック自体は遅い（0.5GHz）
・ベクトル演算器がすさまじく高速
ベクトル計算機
ＣＰＵ（0.5GHz）
スカラー演算器
NEC SX-9
ベクトル演算器
ベクトル演算器
ベクトル演算器
ベクトル演算器
ベクトル演算器
ベクトル演算器
ベクトル演算器
32GB/s
高速アクセス
メモリーバンク
メモリー（GB)
・メモリアクセスが高速（キャッシュメモリは不要）
・さらに、メモリーアクセスを分散させるので高速
（一カ所に集中すると遅くなる；Bank Conflict）
ベクトル演算と計算の加速
DO J = 1, NY
DO I = 1, NX
DxVx (I,J) = ( Vx (I+1,J) - Vx (I 1,J) ) / (2＊Dx)
END DO
END DO
ベクトル演算機構
DOループ処理（NX回の繰り返し）
を一纏め（256個単位）に一気に
ハードウエアで実行する
・ただし、複雑な繰り返し計算（条
件文など）はベクトル化できない
●計算プログラム
＝スカラー計算部分＋ベクトル計算部分
●ベクトル化率＝ベクトル計算部分／プログラム全
体
ベクトル化率＝50% 性能向上率＝2倍
90%
8倍
95%
20倍
99%
40倍
＊ベクトル化による性能向上率を50倍と仮定
・ベクトル化率が100%に限りなく近づいて初めて効果が現れる
FDM計算はベクトル化率が高い
－FDM計算速度が数十倍に加速する
－PSM(FFT)計算はベクトル化が難しい
スカラー計算
PSM～16th-order FDM
ベクトル計算
PSM << 128th-orderFDM
スカラー計算機
・CPUは速い
・キャッシュメモリの速度、容量が命
・外部メモリのアクセス時間は、演算の数十～数百倍遅い！
スカラー計算機
ＣＰＵ（4GHz）
スカラー演算器
高速アクセス
キャッシュメモリ
（MB）
CPU速度（クロック周波数）の変化
・1970年代 Intel 8086、 5 MHz
・1990年代 Intel 80286、 33 MHz
・2000年代 Intel Pentium4、 1.5 GHz
・2007年
3.4 GHz （そろそろ、限界）
ＣＰＵ（4GHz）
スカラー演算器
数 GB/s
低速
・メモリアクセスが非常に遅い
メモリー（GB)
スカラー演算器
高速アクセス
キャッシュメモリ
（MB）
・マルチコア（Dual 、Quad ）化
スカラーコンピュータとキャッシュメモリの効
果
演算プロセッサ
1次キャッシュ
（16KB）
2次キャッシュ(L2)
（512KB）
２次キャッ
シュ溢れ
１次キャッ
シュ溢れ
理論性能と実性能の乖離！
＊キャッシュあふれが起きるとおしまい。
＊ベンチマークテストは速いのに。。。
http:/www.hpc.co.jp
外部メモリ
（2GB)
並列 FDM 計算：２つの並列化思想
地震研EICはこれに近いタイプ
共有メモリ型コンピュータ（とても高価）
(1) Thread Parallel
自動並列コンパイラー
自動並列化、効率は高くない
(2) Message Passing
分散メモリ型、PCクラスタ（いわゆるパソコン）
領域分割
MPIの利用した明示的な通信, 並列効率が高い
Thread Parallel: Doループの自動分割処理
原プログラム
read (6,*) s
do i=1, 128
a(i) = b(i) +
c(i)
end do
write (*,*) a
４CPUに並列化されたプログラム
・自動で並列コードができる。かんたん。
% f77 –parallel prog.f
% a.out
・計算の一部を並列化、残りは全CPUで同じ処理をする
ため、並列効率は高くない。
read (6,*) a, b
並列化
do i=1, 32
do i=33, 64
do i=65, 96
do i=97, 128
a(i) = b(i) + c(i) a(i) = b(i) + c(i) a(i) = b(i) + c(i) a(i) = b(i) + c(i)
end do
end do
end do
end do
write (*,*) c
並列化の効率: アムダールの法則
a
1a
プログラム＝並列化可能な部分＋並列化できない部分
－FDM計算でも結構条件がきつい
－CPUを増やせば良いというものではない
（少数精鋭主義）
●平均加速率: An
An 
－Grid による高速科学計算は幻想（錯覚）
1
a
 (1  a )
n
●並列化効率: En
a
1a
En 
An
n
・並列化効率をEn> 50％を維持するため
のCPU数の上限
並列化率
並列化部分
並列化できない部分
a=99.2%
99.6%
99.8%
99.9%
ＣＰＵ数
n=16
32
64
128
MPIを用いた並列化
MPIで最初から並列化を意識して書いたプログラム
並列化
Input (1:32)
read (*,*) b,c
並列化
read (*,*) b, c
Input (65:92)
read (*,*) b, c
Input (93:128)
read (*,*) b,c
do i=1, 32
do i=1, 32
do i=1, 32
do i=1, 32
a(i) = b(i) + c(i)
a(i) = b(i) + c(i)
a(i) = b(i) + c(i)
a(i) = b(i) + c(i)
end do
end do
end do
end do
write (*,*) a
並列化
Input (33:64)
・処理とデータを予め完全に分割する
・独立に計算。並列効率が良い
・ＣＰＵ間のデータ交換はMPIを用いる
Output (1:32)
write (*,*) a
Output (33:64)
write (*,*) a
Output (65:92)
write (*,*) a
Output (93:128)
ＭＰＩプログラミング（通信）
・MPI（Message Passing Interface）は標準仕様
・MPIプログラムで使うコマンドは意外と少ない
・送信（Send）  受信（Receive）の対応
Program 1
Program 2
Program 3
Send (to=2, b(1))
Send (to=3, b(1))
Send (to=4, b(1))
Receive (from=０, c(32))
Receive (from=1, c(32))
Receive (from=2, c(32))
メモリ内容
b(1)
b(2)
b(3)
c(1)
c(2)
c(3)
b(1)
b(2)
b(3)
c(1)
c(2)
c(3)
b(1)
b(2)
b(3)
c(1)
c(2)
c(3)
b(32)
c(32)
b(32)
c(32)
b(32)
c(32)
ＭＰＩプログラミング（通信２）
・１つのMPIプログラムを全てのCPUで走らせる
・各CPUは番号（Rank)に応じて違う動作をする
% mpif90 prog.f
% mpirun –np 3 a.out
Program 1
id=Ｒａｎｋ（）
１
Program 2
id=Ｒａｎｋ（）
２
Program 3
id=Ｒａｎｋ（）
３
Send (to=id+1,b(1))
Send (to=id+1,b(1))
Send (to=id+1, b(1))
Receive (from=id-1, a(32))
Receive (from=id-1, a(32))
Receive (from=id-1, a(32))
ＭＰＩプログラムの例
ifrom
簡単なMPIプログラムの例
include 'mpif.h' MPIの定義ファイル
double precision STAT(MPI_STATUS_SIZE)
用意されたプロセス数
初期化
call MPI_INIT (ierr)
call MPI_COMM_SIZE (MPI_COMM_WORLD, nproc, ierr)
call MPI_COMM_RANK (MPI_COMM_WORLD, myid, ierr)
ito = myid + 1 送り先（右隣）
ifrom = myid - 1受け元（左隣）
私はそのうちの何番目（0,1,…,nproc-1)
Myid=0
ito
Myid=
nproc-1
Myid=１
★nproc個のCPUで本プログラムを同時実行
（nproc数は実行時に指定する）
★自分の相対番号（myid：0～nproc-1）を知る
★myidをもとに、左右隣の通信の相手を設定
★右隣にデータ送信、左隣からデータ受信
if (ito.gt.(nproc-1)) ito = 0 右隣がなければ、左端へ
if (ifrom.lt.0) ifrom = nproc-1左隣がなければ、右端へ
sdata = myid*1.23
icount = 1
msgtag = 0
送信データ＝自分のプロセス番号＊1.23
送信データ個数＝１
送信メッセージ番号＝０
データ送受信
call MPI_SENDRECV ( sdata, icount, MPI_REAL, ito, msgtag,
:
rdata, icount, MPI_REAL, ifrom, msgtag,
:
MPI_COMM_WORLD, stat, ierr )
自分のプロセス番号と、受信データ、送り主を表示
write (6,*) 'MyID:', myid, ' Rec. Data=', rdata, ' From=', ifrom
call MPI_FINALIZE (ierr)おしまい
end
MPIプログラムの実行方法
%mpif77 test1.f コンパイル
%mpirun -np 2 a.out 実行
MyID: 0 Rec. Data= 1.230 From= 1
MyID: 1 Rec. Data= 0.000 From= 0
% mpirun -np 4 a.out
MyID: 0 Rec. Data=
MyID: 2 Rec. Data=
MyID: 3 Rec. Data=
MyID: 1 Rec. Data=
3.690
1.230
2.460
0.000
From=
From=
From=
From=
3
1
2
0
地球シミュレータ：共有メモリ／分散メモリ型並列計算機
・Node内は共有メモリ型（SMP)
・Node間は分散メモリ型
Node (8CPU)
1 2 3 4 5 6 7 8
1 2 3 4 5 6 7 8
Shared Memory
Shared Memory
640node*640node
Crossbar Switch
1 2 3 4 5 6 7 8
1 2 3 4 5 6 7 8
Shared Memory
Shared Memory
演算速度と
通信速度のバランスが重要
CPUが高速になるほ
ど通信速度が重要
地球シミュレータによる並列ＦＤＭ（差分法）計算
並列計算
並列計算効率
ノード内スレッド並列、ノード間ＭＰＩ
MPI
領域分割
地球シミュレータ
演算効率・並列化度のテスト
ベクトル計算機
（テストモデル）
・4次精度差分法
・３次元領域分割（MPI）
・格子数：1024*512*512
・タイムステップ： 800
CPU数時間[秒]
128
91
1024
16
加速率
5.7
並列化率
最大並列化度
0.999572
2,116
*ベクトル化率＝９９．６％
*演算速度＝3.1TFLOPS (Peak性能比:47%)
スカラー計算機
黙って座れば、性能が出る
CPU数時間[秒]
128
2,520
1024
344
加速率
7.3
並列化率
最大並列化度
0.999895
9,520
*演算速度＝0.14TFLOPS (Peak性能比:5%)
IBM Blue Gene/L
普通のプログラムはこんな程度
並列化率、並列化効率の評価
●加速率：An
（測定）1台のCPUでの実効時間：T1
ープログラムサイズが大きく、
普通は測定できない
An 
1
a
 (1  a )
n
a=.999572(IBM BG/L)
a=.999895 (Earth Simulator)
（測定１） n台での実効時間：Tn
（測定２） m台での実効時間：Tm
●平均加速率: An, Am
T1
1

Tn a  (1  a )
n
T
1
Am  1 
Tm a  (1  a )
m
●並列化率: a
9,520CPU
までOK
An 
a
2,116CPUまで
OK
Tn  Tm
m 1
n 1
Tn 
Tm
m
n
・次世代スパコン（京速コンピュータ）は十万CPU以上になるのは必須
・並列化率を上げるための精進が必要

Download Report