偏微分方程式の差分計算
今回の内容
差分法
1階微分，2階微分に対する差分法
2次元拡散方程式
gnuplotによる結果の表示
if分岐の書き方による実行時間の変化
高速化に利用できるいくつかのテクニック
前回授業
ビットマップを使った画像処理
配列の1要素が物理的な配置に対応
配列の1要素に物理的なデータが定義
数値計算（差分法）
計算機を利用して数学・物理学的問題の解を計算
微積分を計算機で扱える形に変換
差分法では微分を差分に置き換え
処理自体は単純
精度を上げるために計算量が増加
シミュレーションの歴史と進歩
数値計算（差分法）
物理空間に観測点を設置
観測点での物理量の変化を計算
偏微分方程式（拡散方程式）
物質の拡散を表す方程式
水の中に落ちたインクの拡散，金属中の熱伝導等
時刻
t=0におけるfの分布（初期値）が既知
時間進行に伴い，
fがどのように変化するかを計算
時間積分しながら
fの分布を求める
2
2
(
,
)
)
,
(
x
t
x
f
t
t
x
f
2
2
2
2
(
,
,
)
(
,
,
)
)
,
,
(
y
t
y
x
f
x
t
y
x
f
t
t
y
x
f
差分法
計算機で微分を計算する方法の一つ
微分の定義
xの関数fについて，
xだけ離れた2点間の傾きを計算し，
2点の間隔を
無限小に近づけたときの極限
差分近似
関数をある間隔でサンプリング
その間隔
x
が
fの変化に対して十分小さい
と仮定
Δx
x
f
Δx
x
f
dx
df
Δx
)
(
)
(
lim
0
Δx
x
f
Δx
x
f
dx
df
(
)
(
)
差分法（理論的なお話）
差分の誤差
limを排除したことでどの程度の誤差が入るのか
関数
fのテイラー展開を利用
Δx
x
f
Δx
x
f
Δx
x
f
Δx
x
f
dx
df
Δx
)
(
)
(
)
(
)
(
lim
0
2
2
2
3
3
3
!
3
!
2
)
(
)
(
dx
f
d
Δx
dx
f
d
Δx
dx
df
Δx
x
f
Δx
x
f
2
2
2
3
3
3
!
3
!
2
1
)
(
)
(
dx
f
d
Δx
dx
f
d
Δx
Δx
Δx
x
f
Δx
x
f
dx
df
差分法（理論的なお話）
空間打ち切り誤差
定義とテイラー展開の比較
テイラー展開を有限項で打ち切ったことによる誤差
2
2
2
3
3
3
!
3
!
2
1
)
(
)
(
dx
f
d
Δx
dx
f
d
Δx
Δx
Δx
x
f
Δx
x
f
dx
df
Δx
x
f
Δx
x
f
Δx
x
f
Δx
x
f
dx
df
Δx
)
(
)
(
)
(
)
(
lim
0
誤差
2
2
2
3
3
!
3
!
2
dx
f
d
Δx
dx
f
d
Δx
2階微分の離散化
テイラー展開を応用
x方向にx離れた2点で展開
2式を足すと1階微分が消滅
2
2
2
3
3
3
!
3
!
2
)
(
)
(
x
f
Δx
x
f
Δx
x
f
Δx
x
f
Δx
x
f
2
2
2
3
3
3
!
3
!
2
)
(
)
(
x
f
Δx
x
f
Δx
x
f
Δx
x
f
Δx
x
f
2
2
2
!
2
2
)
(
2
)
(
)
(
x
f
Δx
x
f
Δx
x
f
Δx
x
f
2階微分の離散化
2階微分の式に整理
2
2
2
(
)
2
(
)
(
)
Δx
Δx
x
f
x
f
Δx
x
f
x
f
x
f(x)
f(x−
x)
f(x)
x
x=0 ･･･ x−
x
x x+
x
f(x+
x)
2階微分の離散化
2階微分の式に整理
2
2
2
(
)
2
(
)
(
)
Δx
Δx
x
f
x
f
Δx
x
f
x
f
dx
サンプリングされた関数値
を配列f[]で保持
f[i]
f[i‐1]
f[i+1]
2階の中心差分
(f[i+1]‐2*f[i]+f[i‐1])
/(dx*dx)
f[i]
i
2階微分の離散化
サンプリングされた関数値fの保持
x方向長さ（計算領域の長さ）を等間隔xに分割
間隔
xの両端に離散点（観測点）を置き，そこでfの値を定義
原点から順に離散点に番号を付与し，iで表現
離散点iにおけるfの値を配列f[i]に保存
配列f[]
i=0 ･･･
i−1 i
i+1
f[i]
i
離散点
配列（メモリ）
差分法の実装
計算領域内部
d2fdx2[i]=(f[i‐1]‐2*f[i]+f[i+1])/dxdx;
境界条件(関数値が無いため処理を変更)
d2fdx2[0 ]=(2*f[0 ]‐5*f[1 ]+4*f[2 ]‐f[3 ])/dxdx;
d2fdx2[N‐1]=(2*f[N‐1]‐5*f[N‐2]+4*f[N‐3]‐f[N‐4])/dxdx;
d2fdx2[i]
f[i]
＋
＋
＋
＋
＋
2 1 Δx＋
差分法の実装
計算領域内部
d2fdx2[i]=(f[i‐1]‐2*f[i]+f[i+1])/dxdx;
境界条件(関数値が無いため処理を変更)
d2fdx2[0 ]=(2*f[0 ]‐5*f[1 ]+4*f[2 ]‐f[3 ])/dxdx;
d2fdx2[N‐1]=(2*f[N‐1]‐5*f[N‐2]+4*f[N‐3]‐f[N‐4])/dxdx;
d2fdx2[i]
f[i]
＋
2 1 Δx差分法の実装
計算領域内部
d2fdx2[i]=(f[i‐1]‐2*f[i]+f[i+1])/dxdx;
境界条件(関数値が無いため処理を変更)
d2fdx2[0 ]=(2*f[0 ]‐5*f[1 ]+4*f[2 ]‐f[3 ])/dxdx;
d2fdx2[N‐1]=(2*f[N‐1]‐5*f[N‐2]+4*f[N‐3]‐f[N‐4])/dxdx;
d2fdx2[i]
f[i]
2 1 Δx＋
#include<stdlib.h>
#include<math.h>/*‐lmオプションが必要*/
#define Lx
(2.0*M_PI)
#define Nx
(256)
#define dx
(Lx/(Nx‐1))
#define Nbytes
(Nx*sizeof(double))
#define dxdx
(dx*dx)
void init(double *f){
int i;
for(i=0; i<Nx; i++){
f[i] = sin(i*dx);
}
}
void differentiate(double *f,
double *d2fdx2){
int i;
d2fdx2[0] = ( 2.0*f[0]
‐5.0*f[1]
+4.0*f[2]
‐
f[3])/dxdx;
for(i=1; i<Nx‐1; i++)
d2fdx2[i] = ( f[i+1]
‐2.0*f[i
]
+ f[i‐1])/dxdx;
d2fdx2[Nx‐1]=(‐
f[Nx‐4]
+4.0*f[Nx‐3]
‐5.0*f[Nx‐2]
+2.0*f[Nx‐1])/dxdx;
}
int main(void){
double *f,*d2fdx2;
f = (double *)malloc(Nbytes);
d2fdx2 = (double *)malloc(Nbytes);
init(f);
differentiate(f,d2fdx2);
return 0;
}
CPUプログラム
differentiate.c
関数の離散化
計算領域の長さと離散点の数，離散点の間隔の関係
f(x)
x
x
0
2
L
x0からL
xの間に設けられた点の数
N
x=L
x/
x + 1
実行結果
f,
d2f
/dx2
GPUへの移植
計算領域内部を計算するスレッド
d2fdx2[i]=(f[i‐1]‐2*f[i]+f[i+1])/dxdx;
境界を計算するスレッド
d2fdx2[0 ]=(2*f[0 ]‐5*f[1 ]+4*f[2 ]‐f[3 ])/dxdx;
d2fdx2[N‐1]=(2*f[N‐1]‐5*f[N‐2]+4*f[N‐3]‐f[N‐4])/dxdx;
d2fdx2[i]
f[i]
＋
＋
＋
＋
＋
＋
2 1 Δx＋
＋
#include<stdio.h>
#include<stdlib.h>
#include<math.h>/*‐lmオプションが必要*/
#define Lx
(2.0*M_PI)
#define Nx
(256)
#define dx
(Lx/(Nx‐1))
#define Nbytes
(Nx*sizeof(double))
#define dxdx
(dx*dx)
#define NT (128)
#define NB (Nx/NT)
void init(double *f){
int i;
for(i=0; i<Nx; i++){
f[i] = sin(i*dx);
}
}
__global__
void differentiate
(double *f, double *d2fdx2){
int i = blockIdx.x*blockDim.x
+ threadIdx.x;
if(i==0)
d2fdx2[i] = ( 2.0*f[i
]
‐5.0*f[i+1]
+4.0*f[i+2]
‐
f[i+3])/dxdx;
if(0<i && i<Nx‐1)
d2fdx2[i] = ( f[i+1]
‐2.0*f[i
]
+ f[i‐1])/dxdx;
if(i==Nx‐1)
d2fdx2[i]=(‐
f[i‐3]
+4.0*f[i‐2]
‐5.0*f[i‐1]
+2.0*f[i
])/dxdx;
}
GPUプログラム
differentiate.cu
int main(void){
double *host_f,*host_d2fdx2;
double *f,*d2fdx2;
host_f =(double *)malloc(Nbytes);
cudaMalloc((void **)&f,Nbytes);
cudaMalloc((void **)&d2fdx2,Nbytes);
init(host_f);
cudaMemcpy(f, host_f, Nbytes, cudaMemcpyHostToDevice);
differentiate<<<NB, NT>>>(f, d2fdx2);
//host_d2fdx2=(double *)malloc(Nbytes);
//cudaMemcpy(host_d2fdx2, d2fdx2, Nbytes, cudaMemcpyDeviceToHost);
//for(int i=0; i<Nx; i++)printf("%f,%f,%f¥n",i*dx,host_f[i],host_d2fdx2[i]);
free(host_f);
free(host_d2fdx2);
cudaFree(f);
cudaFree(d2fdx2);
return 0;
}
GPUプログラム
differentiate.cu
2次元への拡張
拡散方程式
1次元
2次元
2
2
(
,
)
)
,
(
x
t
x
f
t
t
x
f
2
2
2
2
(
,
,
)
(
,
,
)
)
,
,
(
y
t
y
x
f
x
t
y
x
f
t
t
y
x
f
differentiate.cuで計算
どのように2次元
に拡張するか
どのように離散化
するか
2次元への拡張
x方向2階偏微分
y方向を固定してx方向に偏微分
y方向2階偏微分
x方向を固定してy方向に偏微分
2
2
2
(
,
)
2
(
,
)
(
,
)
Δx
y
Δx
x
f
y
x
f
y
Δx
x
f
x
f
2
2
2
(
,
)
2
(
,
)
(
,
)
Δy
Δy
y
x
f
y
x
f
Δy
y
x
f
y
f
時間積分
時間微分項の離散化
時間微分項を前進差分で離散化
右辺の
t+tの項を移行
Δt
t
y
x
f
Δt
t
y
x
f
t
f
(
,
,
)
(
,
,
)
f
t
f
2
t
f
Δt
t
y
x
f
Δt
t
y
x
f
)
(
,
,
)
,
,
(
拡散方程式を代入
)
,
,
(
)
,
,
(
)
,
,
(
x
y
t
Δt
f
x
y
t
Δt
2
f
x
y
t
f
f(x, y, t)の2階微分を計算できればf(x, y, t+
t)が求められる
離散化された方程式の記述
簡略化した表現
配列との対応をとるため下付き添字
i, jを利用
時間は，上付き添字
nを利用
j
i
f
y
x
f
(
,
)
,
j
i
f
y
Δx
x
f
(
,
)
1
,
1
,
)
,
(
x
y
Δy
f
i
j
f
n
j
i
f
t
y
x
f
(
,
,
)
,
1
,
)
,
,
(
n
j
i
f
Δt
t
y
x
f
離散化された拡散方程式
連続系
離散系
t秒後の値
2
2
2
2
(
,
,
)
(
,
,
)
)
,
,
(
y
t
y
x
f
x
t
y
x
f
t
t
y
x
f
2
1
,
,
1
,
2
,
1
,
,
1
,
1
,
2
2
Δy
f
f
f
Δx
f
f
f
Δt
f
f
n
j
i
n
j
i
n
j
i
n
j
i
n
j
i
n
j
i
n
j
i
n
j
i
2
1
,
,
1
,
2
,
1
,
,
1
,
1
,
2
2
Δy
f
f
f
Δx
f
f
f
Δt
f
f
n
j
i
n
j
i
n
j
i
n
j
i
n
j
i
n
j
i
n
j
i
n
j
i
離散化された関数値の保持
各方向を等間隔
x,
yに分割し，離散点を配置
離散点で関数値を定義し，配列に保存
配列f[][]
i
離散点
配列（メモリ）
j=0
･･･
j
−
1
j
j+1
拡散方程式（熱伝導方程式）
x
y
i
i+1
i−1
j
j+1
j−1
x
y
x
y
x
y
t
t+
t
2
1
,
,
1
,
2
,
1
,
,
1
,
1
,
2
2
Δy
f
f
f
Δx
f
f
f
Δt
f
f
n
j
i
n
j
i
n
j
i
n
j
i
n
j
i
n
j
i
n
j
i
n
j
i
i
i+1
i−1
j
j+1
j−1
拡散方程式の安定性
プログラムを正しく作成しても正常な計算結果が得られ
ない場合がある
安定条件
拡散の強さを表す係数
（拡散係数）を使った形
二つの条件を満たすことが必要条件（結果が正しいかは別）
5
.
0
2
Δx
Δt
v
2
0
.
5
Δy
Δt
v
2
1
,
,
1
,
2
,
1
,
,
1
,
1
,
2
2
Δy
f
f
f
Δx
f
f
f
Δt
f
f
n
j
i
n
j
i
n
j
i
n
j
i
n
j
i
n
j
i
n
j
i
n
j
i
計算手順
1.
計算条件の決定
計算領域の大きさ
L
x
, L
y
計算領域の分割数（離散点の個数）
N
x
, N
y
離散点同士の間隔（格子間隔）
x, y
計算時間間隔
t
2.
初期値の決定
fの初期分布の決定
3.
差分値の計算
fの分布からx, y方向の2階微分値を計算
境界条件に基づいて境界の値を決定
t秒後のfを計算
CPUプログラム
計算条件の決定
計算領域の大きさ
L
x
, L
y
計算領域の分割数
（離散点の個数）
N
x
, N
y
離散点同士の間隔
（格子間隔）
x, y
計算時間間隔
t
#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#define Lx (1.0)
#define Ly (1.0)
#define Nx 128
#define Ny 128
#define dx (Lx/(Nx‐1))
#define dy (Ly/(Ny‐1))
#define dt 0.0001
CPUプログラム
初期条件
境界条件
)
2
.
0
(
1
)
2
.
0
(
0
r
r
f
1
1
x
y
r
0
n
f
n
n
n: 外向き法線ベクトル
j i j i j i j if
f
Δx
f
f
, 1 , 1 , 1 , 10
2
1 , 1 , 1 , 1 ,0
2
j i j i j i j if
f
Δy
f
f
2 2y
x
r
#include<stdlib.h> #include<math.h> #define Lx (1.0) #define Ly (1.0) #define Nx 128 #define Ny 128 #define dx (Lx/(Nx‐1)) #define dy (Ly/(Ny‐1)) #define dt 0.0001 #define endT (1.0) #define Nt (int)(endT/dt) #define DIFF (0.01) #define dxdx (dx*dx) #define dydy (dy*dy) #define Nbytes (Nx*Ny*sizeof(double)) void init(double *, double *, double *); void laplacian(double *, double *); void integrate(double *, double *, double *); void update(double *, double *); int main(void){ double *f,*f_new,*f_lap,x,y; int i,j,n; f = (double *)malloc(Nx*Ny*sizeof(double)); f_new = (double *)malloc(Nx*Ny*sizeof(double)); f_lap = (double *)malloc(Nx*Ny*sizeof(double)); init(f, f_lap, f_new); for(n=0;n<Nt;n++){ laplacian(f,f_lap); integrate(f,f_lap,f_new); update(f,f_new); } return 0; }
CPUプログラム
diffusion.c
void init(double *f, double *f_lap, double *f_new){ int i,j,ij; double x,y; for(j=0;j<Ny;j++){ for(i=0;i<Nx;i++){ f [i+Nx*j] = 0.0; f_lap[i+Nx*j] = 0.0; f_new[i+Nx*j] = 0.0; x=i*dx‐Lx/2.0; y=j*dy‐Ly/2.0; if(sqrt(x*x+y*y)<=0.2) f[i+Nx*j] = 1.0; } } } void laplacian(double *f, double *f_lap){ int i,j,ij,ip1j,im1j,ijp1,ijm1; for(j=0;j<Ny;j++){ for(i=0;i<Nx;i++){ ij = i+Nx*j; im1j=i‐1+Nx* j; if(i‐1< 0 )im1j=i+1+Nx*j; ip1j=i+1+Nx* j; if(i+1>=Nx)ip1j=i‐1+Nx*j; ijm1=i +Nx*(j‐1);if(j‐1< 0 )ijm1=i+Nx*(j+1); ijp1=i +Nx*(j+1);if(j+1>=Ny)ijp1=i+Nx*(j‐1); f_lap[ij] = (f[ip1j]‐2.0*f[ij]+f[im1j])/dxdx +(f[ijp1]‐2.0*f[ij]+f[ijm1])/dydy; } } } void integrate (double *f, double *f_lap, double *f_new){ int i,j,ij; for(j=0;j<Ny;j++){ for(i=0;i<Nx;i++){ ij = i+Nx*j; f_new[ij] = f[ij] + dt*DIFF*f_lap[ij]; } } } void update(double *f, double *f_new){ int i,j,ij; for(j=0;j<Ny;j++){ for(i=0;i<Nx;i++){ ij = i+Nx*j; f[ij] = f_new[ij]; } } }
CPUプログラム
diffusion.c
CPUプログラム
差分計算
x方向，y方向偏微分を個別に計算して加算
void laplacian(double *f, double *f_lap){
int i,j,ij,ip1j,im1j,ijp1,ijm1;
for(j=0;j<Ny;j++){
for(i=0;i<Nx;i++){
ij
=i
+Nx* j;
im1j=i‐1+Nx* j; if(i‐1< 0 )im1j=i+1+Nx*j;
ip1j=i+1+Nx* j; if(i+1>=Nx)ip1j=i‐1+Nx*j;
ijm1=i
+Nx*(j‐1);if(j‐1< 0 )ijm1=i+Nx*(j+1);
ijp1=i
+Nx*(j+1);if(j+1>=Ny)ijp1=i+Nx*(j‐1);
f_lap[ij] = (f[ip1j]‐2.0*f[ij]+f[im1j])/dxdx
+(f[ijp1]‐2.0*f[ij]+f[ijm1])/dydy;
}
}
}
ラプラシアン計算のメモリ参照
ある1点i,jのラプラシアンを計算するために，周囲5点
のfを参照
f[]
f_lap[]
i−1 i i+1
j−
1
jj
+1
ラプラシアン計算のメモリ参照
ある1点i,jのラプラシアンを計算するために，周囲5点
のfを参照
i−1 i i+1
j−
1
jj
+1
f[]
f_lap[]
ラプラシアン計算のメモリ参照
ある1点i,jのラプラシアンを計算するために，周囲5点
のfを参照
i−1 i i+1
j−
1
jj
+1
f[]
f_lap[]
ラプラシアン計算のメモリ参照
ある1点i,jのラプラシアンを計算するために，周囲5点
のfを参照
i−1 i i+1
j−
1
jj
+1
f[]
f_lap[]
ラプラシアン計算のメモリ参照
ある1点i,jのラプラシアンを計算するために，周囲5点
のfを参照
全てのfを参照し，領域内部のf_lapを計算
CPUプログラム
境界条件
法線方向の1階微分が0
境界での2階の差分式
f_lap[]
0
n
f
n: 外向き法線ベクトル
j i j i j i j if
f
Δx
f
f
, 1 , 1 , 1 , 10
2
1 , 1 , 1 , 1 ,0
2
j i j i j i j if
f
Δy
f
f
2 , , 1 2 , , 1 2 , 1 , , 12
2
or
2
2
2
Δx
f
f
Δx
f
f
Δx
f
f
f
i j i j i j i j
i j i j
i j
2 , 1 , 2 , 1 , 2 1 , , 1 ,2
2
or
2
2
2
Δy
f
f
Δy
f
f
Δy
f
f
f
i j i j i j i j
i j i j
i j
CPUプログラム
境界条件
参照する点が境界からはみ出して
いる場合は添字を修正
2階微分を計算する式は変更なし
右端ではi+1をi‐1に修正
左端ではi‐1をi+1に修正
上端ではj+1をj‐1に修正
下端ではj‐1をj+1に修正
f_lap[]
if(i‐1< 0 )im1j=i+1+Nx*j;
if(i+1>=Nx)ip1j=i‐1+Nx*j;
if(j‐1< 0 )ijm1=i
+Nx*(j+1);
if(j+1>=Ny)ijp1=i
+Nx*(j‐1);
CPUプログラム
fの積分
void integrate(double *f, double *f_lap, double *f_new){
int i,j,ij;
for(j=0;j<Ny;j++){
for(i=0;i<Nx;i++){
ij = i+Nx*j;
f_new[ij] = f[ij] + dt*DIFF*f_lap[ij];
}
}
}
2
1
,
,
1
,
2
,
1
,
,
1
,
1
,
2
2
Δy
f
f
f
Δx
f
f
f
Δt
f
f
n
j
i
n
j
i
n
j
i
n
j
i
n
j
i
n
j
i
n
j
i
n
j
i
CPUプログラム
fの更新
f
n
から
f
n+1
を計算
f
n+1
から
f
n+2
を計算
今の時刻から次の時刻を求める
求められた次の時刻を今の時刻と見なし，次の時刻を求める
void update(double *f, double *f_new){
int i,j,ij;
for(j=0;j<Ny;j++){
for(i=0;i<Nx;i++){
ij = i+Nx*j;
f[ij] = f_new[ij];
}
}
}
同じアルゴリズム
#include<stdio.h> #include<stdlib.h> #include<math.h> #define Lx (1.0) #define Ly (1.0) #define Nx 2048 #define Ny Nx #define dx (Lx/(Nx‐1)) #define dy (Ly/(Ny‐1)) #define dt 0.000001 #define endT (1.0) #define Nt (int)(endT/dt) #define DIFF (0.01) #define dxdx (dx*dx) #define dydy (dy*dy) #define Nbytes (Nx*Ny*sizeof(double)) #include "dif1.cu" //#include "dif2.cu" int main(void){ int n; double *dev_f,*dev_f_new,*dev_f_lap; dim3 Thread, Block; if(DIFF*dt/dxdx > 0.5){ printf("configuration error¥n"); exit(1); } cudaMalloc( (void**)&dev_f , Nbytes ); cudaMalloc( (void**)&dev_f_new, Nbytes ); cudaMalloc( (void**)&dev_f_lap, Nbytes ); Thread = dim3(THREADX,THREADY,1); Block = dim3(BLOCKX ,BLOCKY, 1); init<<<Block, Thread>>> (dev_f, dev_f_lap, dev_f_new); for(n=1;n<=Nt;n++){ laplacian<<<Block, Thread>>>(dev_f,dev_f_lap); integrate<<<Block, Thread>>> (dev_f,dev_f_lap,dev_f_new); update<<<Block, Thread>>>(dev_f,dev_f_new); } return 0; }
GPUプログラム（CPU処理用共通部分）
diffusion.cu
#define THREADX 16 #define THREADY 16 #define BLOCKX (Nx/THREADX) #define BLOCKY (Ny/THREADY) __global__ void laplacian(double *f, double *f_lap){ int i,j,ij,im1j,ip1j,ijm1,ijp1; i = blockIdx.x*blockDim.x + threadIdx.x; j = blockIdx.y*blockDim.y + threadIdx.y; ij =i +Nx* j; im1j=i‐1+Nx* j; if(i‐1< 0 )im1j=i+1+Nx*j; ip1j=i+1+Nx* j; if(i+1>=Nx)ip1j=i‐1+Nx*j; ijm1=i +Nx*(j‐1);if(j‐1< 0 )ijm1=i+Nx*(j+1); ijp1=i +Nx*(j+1);if(j+1>=Ny)ijp1=i+Nx*(j‐1); f_lap[ij] = (f[ip1j]‐2.0*f[ij]+f[im1j])/dxdx +(f[ijp1]‐2.0*f[ij]+f[ijm1])/dydy; } __global__ void integrate (double *f, double *f_lap, double *f_new){ int i,j,ij; i = blockIdx.x*blockDim.x + threadIdx.x; j = blockIdx.y*blockDim.y + threadIdx.y; ij = i+Nx*j; f_new[ij] = f[ij] + dt*DIFF*f_lap[ij]; } __global__ void update(double *f, double *f_new){ int i,j,ij; i = blockIdx.x*blockDim.x + threadIdx.x; j = blockIdx.y*blockDim.y + threadIdx.y; ij = i+Nx*j; f[ij] = f_new[ij]; } __global__ void init (double *f, double *f_lap, double *f_new){ int i,j,ij; double x,y; i = blockIdx.x*blockDim.x + threadIdx.x; j = blockIdx.y*blockDim.y + threadIdx.y; ij = i+Nx*j; f [ij] = 0.0; f_lap[ij] = 0.0; f_new[ij] = 0.0; x=i*dx‐Lx/2.0; y=j*dy‐Ly/2.0; if(sqrt(x*x+y*y)<=0.2) f[ij] = 1.0; }
GPUプログラム（1スレッドが1点を計算）
dif1.cu
gnuplotによる結果の表示
2次元，3次元データをプロットするアプリケーション
コマンドラインで命令を実行してグラフを描画
関数の描画，ファイルから読み込んだデータの表示が可能
tesla??では正しく動作しないため，grouseで実行する
こと
gnuplotによる結果の表示
表示可能なファイルフォーマット
データはスペース区切り
3次元表示では列（または行）の区切りとして空白行が必要
‐0.500000 ‐0.500000 0.000000e+00
‐0.499022 ‐0.500000 0.000000e+00
...
0.499022 ‐0.500000 0.000000e+00
0.500000 ‐0.500000 0.000000e+00
<‐改行
‐0.500000 ‐0.499022 0.000000e+00
‐0.499022 ‐0.499022 0.000000e+00
...
0.499022 ‐0.499022 0.000000e+00
0.500000 ‐0.499022 0.000000e+00
<‐改行
‐0.500000 ‐0.498045 0.000000e+00
‐0.499022 ‐0.498045 0.000000e+00
gnuplotによる結果の表示
スクリプトの実行
gnuplotのコマンドを記述したファイルを実行
load 'スクリプトファイル名'
anim.gpl
結果を単色の等値面で3次元表示
anim_color.gpl
結果をカラーの等値面で3次元表示
anim_2d.gpl
結果をカラーの等高線で2次元表示
gnuplotによる結果の表示
スクリプトanim.gplの内容
set xrange [‐0.5:0.5]
x軸の表示範囲を‐0.5~0.5に固定
set yrange [‐0.5:0.5]
y軸の表示範囲を‐0.5~0.5に固定
set zrange [0:1]
z軸の表示範囲を0~1に固定
set ticslevel 0
xy平面とz軸の最小値表示位置の差
set hidden3d
隠線消去をする
set view 45,30
視点をx方向45°,y方向30°に設定
set size 1,0.75
グラフの縦横比を1:0.75
unset contour
2次元等高線は表示しない
set surface
3次元で等値面を表示
unset pm3d
カラー表示はしない
以下でファイルを読み込み，3次元表示
splot 'f0.00.txt' using 1:2:3 with line title "t=0.00s"
...
splot 'f1.00.txt' using 1:2:3 with line title "t=1.00s"
gnuplotによる結果の表示
gnuplotによる結果の表示
スクリプトanim_color.gplの内容
set xrange [‐0.5:0.5]
x軸の表示範囲を‐0.5~0.5に固定
set yrange [‐0.5:0.5]
y軸の表示範囲を‐0.5~0.5に固定
set zrange [0:1]
z軸の表示範囲を0~1に固定
set ticslevel 0
xy平面とz軸の最小値表示位置の差
set hidden3d
隠線消去をする
set view 45,30
視点をx方向45°,y方向30°に設定
set size 1,0.75
グラフの縦横比を1:0.75
unset contour
2次元等高線は表示しない
set surface
3次元で等値面を表示
set pm3d
カラー表示する
set cbrange[0:1]
カラーバーの表示範囲を0~1に固定
以下でファイルを読み込み，3次元表示
splot 'f0.00.txt' using 1:2:3 with line title "t=0.00s"
...
splot 'f1.00.txt' using 1:2:3 with line title "t=1.00s"
gnuplotによる結果の表示
gnuplotによる結果の表示
スクリプトanim_2d.gplの内容
set xrange [‐0.5:0.5]
x軸の表示範囲を‐0.5~0.5に固定
set yrange [‐0.5:0.5]
y軸の表示範囲を‐0.5~0.5に固定
set zrange [0:1]
z軸の表示範囲を0~1に固定
set view 0,0
真上から見下ろす
set size 1,1
グラフの縦横比を1:1
set contour
2次元等高線を表示
unset surface
3次元で等値面を表示しない
unset pm3d
カラー表示しない
set cntrparam levels incremental 0,0.1,1
等高線を0から1まで0.1刻みに設定
以下でファイルを読み込み，3次元表示
splot 'f0.00.txt' using 1:2:3 with line title "t=0.00s"
...
splot 'f1.00.txt' using 1:2:3 with line title "t=1.00s"
gnuplotによる結果の表示
様々な境界条件の計算
共有メモリを利用してキャッシュを模擬
共有メモリに付加的な領域を追加
袖領域
様々な境界条件の計算
共有メモリを利用してキャッシュを模擬
共有メモリに付加的な領域を追加
様々な境界条件の計算
共有メモリを利用してキャッシュを模擬
共有メモリに付加的な領域を追加
f[]
sf[][]
f_lap[]
付加的な領域（袖領域）の取り扱い
データがグローバルメモリに存在する場合は，グローバ
ルメモリから読み込み
f[]
sf[][]
f_lap[]
付加的な領域（袖領域）の取り扱い
グローバルメモリに無い場合は境界条件から決定
f[]
sf[][]
境界条件
1階微分が0
0
2
1 , 1 ,
Δy
f
f
y
f
i j i j 1 , 1 ,j i j if
f
0
2
, 1 , 1
Δx
f
f
x
u
i j i j j i j if
f
1,
1,境界条件
2階微分が0
0
2
2 , 1 , , 1 2 2
Δx
f
f
f
x
f
i j i j i j j i j i j if
f
f
1,
2
,
1, 1 , , 1 ,j2
i j i j if
f
f
0
2
2 1 , , 1 , 2 2
Δy
f
f
f
y
f
i j i j i j境界条件
2階微分を片側差分で計算
3 , 2 , 1 , , 1 ,j
4
i j
6
i j
4
i j
i j if
f
f
f
f
2 3 , 2 , 1 , , 1 , 2 22
5
4
Δy
f
f
f
f
y
f
i j i j i j i j j i
2 , 3 , 2 , 1 , , 1 2 22
5
4
Δx
f
f
f
f
x
f
i j i j i j i j j i
j i j i j i j i j if
f
f
f
f
1,
4
,
6
1,
4
2,
3,#define THREADX 16 #define THREADY 16 #define BLOCKX (Nx/THREADX) #define BLOCKY (Ny/THREADY) __global__ void laplacian(double *f, double *f_lap){ int i,j,ij; int tx,ty; __shared__ float sf[1+THREADX+1][1+THREADY+1]; i = blockIdx.x*blockDim.x + threadIdx.x; j = blockIdx.y*blockDim.y + threadIdx.y; tx = threadIdx.x + 1; ty = threadIdx.y + 1; ij = i+Nx*j; sf[tx][ty] = f[ij]; __syncthreads(); if(blockIdx.x == 0 && threadIdx.x == 0) sf[tx‐1][ty] = sf[tx+1][ty]; if(blockIdx.x != 0 && threadIdx.x == 0) sf[tx‐1][ty] = f[i‐1+Nx*j]; if(blockIdx.x == gridDim.x‐1 && threadIdx.x == blockDim.x‐1) sf[tx+1][ty] = sf[tx‐1][ty]; if(blockIdx.x != gridDim.x‐1 && threadIdx.x == blockDim.x‐1) sf[tx+1][ty] = f[i+1+Nx*j]; if(blockIdx.y == 0 && threadIdx.y == 0) sf[tx][ty‐1] = sf[tx][ty+1]; if(blockIdx.y != 0 && threadIdx.y == 0) sf[tx][ty‐1] = f[i+Nx*(j‐1)]; if(blockIdx.y == gridDim.y‐1 && threadIdx.y == blockDim.y‐1) sf[tx][ty+1] = sf[tx][ty‐1]; if(blockIdx.y != gridDim.y‐1 && threadIdx.y == blockDim.y‐1) sf[tx][ty+1] = f[i+Nx*(j+1)]; __syncthreads(); f_lap[ij] = (sf[tx‐1][ty ]‐2.0*sf[tx][ty]+sf[tx+1][ty ])/dxdx +(sf[tx ][ty‐1]‐2.0*sf[tx][ty]+sf[tx ][ty+1])/dydy; }
GPUプログラム（ラプラシアンkernel）
dif2.cu
GPUプログラム（ラプラシアンkernel）
ブロック内のスレッド数＋袖領域分の共有メモ
リを確保
袖領域があるために添字の対応が変化
添字の対応を考えないと，必要なデータを袖領域に置い
てしまう
__syncthreads()を呼んでスレッドを同期
共有メモリにデータが正しく書き込まれた事を保証
__shared__ float sf[
1+
THREADX
+1
][
1+
THREADY
+1
];
i = blockIdx.x*blockDim.x + threadIdx.x;
j = blockIdx.y*blockDim.y + threadIdx.y;
tx = threadIdx.x;
ty = threadIdx.y;
ij = i+Ny*j;
sf[tx][ty] = f[ij];
__syncthreads();
f[]
sf[][]
GPUプログラム（ラプラシアンkernel）
ブロック内のスレッド数＋袖領域分の共有メモ
リを確保
袖領域があるために添字の対応が変化
添字の対応を考えないと，必要なデータを袖領域に置い
てしまう
__syncthreads()を呼んでスレッドを同期
共有メモリにデータが正しく書き込まれた事を保証
__shared__ float sf[
1+
THREADX
+1
][
1+
THREADY
+1
];
i = blockIdx.x*blockDim.x + threadIdx.x;
j = blockIdx.y*blockDim.y + threadIdx.y;
tx = threadIdx.x
+ 1
;
ty = threadIdx.y;
ij = i+Ny*j;
sf[tx][ty] = f[ij];
__syncthreads();
f[]
sf[][]
GPUプログラム（ラプラシアンkernel）
ブロック内のスレッド数＋袖領域分の共有メモ
リを確保
袖領域があるために添字の対応が変化
添字の対応を考えないと，必要なデータを袖領域に置い
てしまう
__syncthreads()を呼んでスレッドを同期
共有メモリにデータが正しく書き込まれた事を保証
__shared__ float sf[
1+
THREADX
+1
][
1+
THREADY
+1
];
i = blockIdx.x*blockDim.x + threadIdx.x;
j = blockIdx.y*blockDim.y + threadIdx.y;
tx = threadIdx.x
+ 1
;
ty = threadIdx.y
+ 1
;
ij = i+Ny*j;
sf[tx][ty] = f[ij];
__syncthreads();
f[]
sf[][]
+1
GPUプログラム（ラプラシアンkernel）
袖領域の設定
if(blockIdx.x == 0 && threadIdx.x == 0 ) sf[tx‐1][ty] = sf[tx+1][ty]; if(blockIdx.x != 0 && threadIdx.x == 0 ) sf[tx‐1][ty] = f[i‐1+Nx*j];
if(blockIdx.x == gridDim.x‐1 && threadIdx.x == blockDim.x‐1) sf[tx+1][ty] = sf[tx‐1][ty]; if(blockIdx.x != gridDim.x‐1 && threadIdx.x == blockDim.x‐1) sf[tx+1][ty] = f[i+1+Nx*j];
if(blockIdx.y == 0 && threadIdx.y == 0 ) sf[tx][ty‐1] = sf[tx][ty+1]; if(blockIdx.y != 0 && threadIdx.y == 0 ) sf[tx][ty‐1] = f[i+Nx*(j‐1)];
if(blockIdx.y == gridDim.y‐1 && threadIdx.y == blockDim.y‐1) sf[tx][ty+1] = sf[tx][ty‐1]; if(blockIdx.y != gridDim.y‐1 && threadIdx.y == blockDim.y‐1) sf[tx][ty+1] = f[i+Nx*(j+1)];
__syncthreads();
グローバルメモリにデータがある箇所はグロー
バルメモリから読み込み
グローバルメモリにデータがない箇所は2階微
分が0になるように袖領域の値を決定
ブロックが境界に接しているか否かで処理を切替
GPUプログラム（ラプラシアンkernel）
共有メモリのデータを利用してラプラシアンを計算
if分岐を排除
f_lap[ij] = (sf[tx‐1][ty ]‐2.0*sf[tx][ty]+sf[tx+1][ty ])/dxdx
+(sf[tx ][ty‐1]‐2.0*sf[tx][ty]+sf[tx ][ty+1])/dydy;
f_lap[]
sf[][]
GPUプログラムの評価
共有メモリを使用したdif2.cu
dif1.cu（共有メモリ不使用）と比較して遅い
Fermi世代以降のGPUはキャッシュを搭載
共有メモリを使っても速くならない
共有メモリへ明示的にデータを移動
余分な処理により負荷が増加
せっかくなのでif文の書き方について評価
GPUプログラムの評価
if分岐の書き方による実行速度の変化
ブロックとスレッドの条件を同時に記述
Warp内のスレッドが分岐すると実行速度が著しく低下
GPUはブロック単位で分岐，Warp単位で分岐しても実行速度
の低下は抑えられる
上の書き方は最適ではない？
ブロック単位の分岐とスレッド単位の分岐を同時に記述
if(blockIdx.x == 0 && threadIdx.x == 0 ) sf[tx‐1][ty] = sf[tx+1][ty]; if(blockIdx.x != 0 && threadIdx.x == 0 ) sf[tx‐1][ty] = f[i‐1+Nx*j];
if(blockIdx.x == gridDim.x‐1 && threadIdx.x == blockDim.x‐1) sf[tx+1][ty] = sf[tx‐1][ty]; if(blockIdx.x != gridDim.x‐1 && threadIdx.x == blockDim.x‐1) sf[tx+1][ty] = f[i+1+Nx*j];
if(blockIdx.y == 0 && threadIdx.y == 0 ) sf[tx][ty‐1] = sf[tx][ty+1]; if(blockIdx.y != 0 && threadIdx.y == 0 ) sf[tx][ty‐1] = f[i+Nx*(j‐1)];
if(blockIdx.y == gridDim.y‐1 && threadIdx.y == blockDim.y‐1) sf[tx][ty+1] = sf[tx][ty‐1]; if(blockIdx.y != gridDim.y‐1 && threadIdx.y == blockDim.y‐1) sf[tx][ty+1] = f[i+Nx*(j+1)];
GPUプログラムの評価
if分岐の書き方による実行速度の変化
ブロックとスレッドの条件を分け，ifを入れ子に
ブロック単位で分岐した後，スレッド単位で分岐
ブロック単位での分岐による速度低下が抑えられそうな気がする
ブロック単位の分岐を記述する8個のifは2個一組（同じ色の
行）で排他的に分岐
blockIdx.x==0が成立すると，blockIdx!=0は絶対に成立しない
if(blockIdx.x == 0 )if(threadIdx.x == 0 ) sf[tx‐1][ty] = sf[tx+1][ty]; if(blockIdx.x != 0 )if(threadIdx.x == 0 ) sf[tx‐1][ty] = f[i‐1+Nx*j];
if(blockIdx.x == gridDim.x‐1)if(threadIdx.x == blockDim.x‐1) sf[tx+1][ty] = sf[tx‐1][ty]; if(blockIdx.x != gridDim.x‐1)if(threadIdx.x == blockDim.x‐1) sf[tx+1][ty] = f[i+1+Nx*j];
if(blockIdx.y == 0 )if(threadIdx.y == 0 ) sf[tx][ty‐1] = sf[tx][ty+1]; if(blockIdx.y != 0 )if(threadIdx.y == 0 ) sf[tx][ty‐1] = f[i+Nx*(j‐1)];
if(blockIdx.y == gridDim.y‐1)if(threadIdx.y == blockDim.y‐1) sf[tx][ty+1] = sf[tx][ty‐1]; if(blockIdx.y != gridDim.y‐1)if(threadIdx.y == blockDim.y‐1) sf[tx][ty+1] = f[i+Nx*(j+1)];
GPUプログラムの評価
if分岐の書き方による実行速度の変化
2個一組のifの片方をelseに書き換え
2個一組のifにおける条件判断は1回限り
2回ifの条件判断を行うよりは速いような気がする
if(blockIdx.x == 0 ){if(threadIdx.x == 0 ) sf[tx‐1][ty] = sf[tx+1][ty];} else {if(threadIdx.x == 0 ) sf[tx‐1][ty] = f[i‐1+Nx*j];}
if(blockIdx.x == gridDim.x‐1){if(threadIdx.x == blockDim.x‐1) sf[tx+1][ty] = sf[tx‐1][ty];} else {if(threadIdx.x == blockDim.x‐1) sf[tx+1][ty] = f[i+1+Nx*j];}
if(blockIdx.y == 0 ){if(threadIdx.y == 0 ) sf[tx][ty‐1] = sf[tx][ty+1];} else {if(threadIdx.y == 0 ) sf[tx][ty‐1] = f[i+Nx*(j‐1)];}
if(blockIdx.y == gridDim.y‐1){if(threadIdx.y == blockDim.y‐1) sf[tx][ty+1] = sf[tx][ty‐1];} else {if(threadIdx.y == blockDim.y‐1) sf[tx][ty+1] = f[i+Nx*(j+1)];}
