Microsoft PowerPoint - 先端GPGPUシミュレーション工学特論(web).pptx

偏微分方程式の差分計算

（拡散方程式）

今日の内容



シミュレーションの歴史と進歩



差分法



_{1階微分，2階微分に対する差分法}



_{1次関数の差分}



_{2次元拡散方程式}



付録



共有メモリの典型的な使い方

数値計算



計算機を利用して数学・物理学的問題の解を計算



微積分を計算機で扱える形に変換



処理自体はあまり複雑ではない

シミュレーションの歴史と進歩

シミュレーションの歴史と進歩



_1985年



_{2次元でのシミュレーション}

シミュレーションの歴史と進歩



_1987年

シミュレーションの歴史と進歩



_{1988年～1990年}

シミュレーションの歴史と進歩



_1993年



空気抵抗を誤差1%で予測可能

シミュレーションの歴史と進歩



_{1992年～1993年}



車体から発生する騒音のシミュレーション



ドアミラーやピラーの形状の改良



エンジンルームの冷却



_{10mm以上の部品は全て含めて解析}

シミュレーションの歴史と進歩

現象を支配する方程式（支配方程式）



現象は微分方程式によって記述



微分と積分が計算できれば現象を明らかにできる



支配方程式



場所や時間によって方程式は変わらない



流体の支配方程式

0













i i

x

u

t





j ij i j i j i

x

x

p

x

u

u

t

u





























i ij i i i

x

q

x

u

x

p

u

x

E

u

t

E



































（質量保存式）

（エネルギ保存式）

（運動量保存式）

拡散方程式



物質の拡散を表す方程式



水の中に落ちたインクの拡散，金属中の熱伝導等



時刻

t=0におけるuの分布（初期値）が既知



時間進行に伴い，

uがどのように変化するかを計算



時間積分しながら

uの分布を求める

2 2

₍

_,

₎

)

,

(

x

t

x

u

t

t

x

u











2 2 2 2

₍

_,

_,

₎

₍

_,

_,

₎

)

,

,

(

y

t

y

x

u

x

t

y

x

u

t

t

y

x

u

















差分法



計算機で微分を計算する方法の一つ



微分の定義



xの関数uについて，



xだけ離れた2点間の傾きを計算し，

₂

点の間隔を無限小に近づけたときの極限



差分近似



関数をある間隔でサンプリング



その間隔



x

が

uの変化に対して十分小さい

と仮定

Δx

x

u

Δx

x

u

dx

du

Δx

)

(

)

(

lim

0









Δx

x

u

Δx

x

u

dx

du

(



)



(

)



差分法（理論的なお話）



差分の誤差



limを排除したことでどの程度の誤差が入るのか



関数

uのテイラー展開を利用

Δx

x

u

Δx

x

u

Δx

x

u

Δx

x

u

dx

du

Δx

)

(

)

(

)

(

)

(

lim

0





























2 2₂ 3 3₃

!

3

!

2

)

(

)

(

dx

u

d

Δx

dx

u

d

Δx

dx

du

Δx

x

u

Δx

x

u































2 2₂ 3 3₃

_

!

3

!

2

1

)

(

)

(

dx

u

d

Δx

dx

u

d

Δx

Δx

Δx

x

u

Δx

x

u

dx

du

差分法（理論的なお話）



空間打ち切り誤差



定義とテイラー展開の比較



テイラー展開を有限項で打ち切ったことによる誤差































2 2₂ 3 3₃

_

!

3

!

2

1

)

(

)

(

dx

u

d

Δx

dx

u

d

Δx

Δx

Δx

x

u

Δx

x

u

dx

du

Δx

x

u

Δx

x

u

Δx

x

u

Δx

x

u

dx

du

Δx

)

(

)

(

)

(

)

(

lim

0















誤差

























2₂ 2 3₃

_

!

3

!

2

dx

u

d

Δx

dx

u

d

Δx

差分法（理論的なお話）



誤差の主要項



空間打ち切り誤差の中で最も大きい誤差は第1項





xは小さい→



xの高次項はさらに小さく，無視できる



直感的に導いた微分の数値計算法の誤差は

O(



x)





xを1/10にすれば誤差も1/10になる

)

(

!

3

!

2

3 3 2 2 2

Δx

O

dx

u

d

Δx

dx

u

d

Δx

_

























_

差分法（理論的なお話）



差分の取り方



関数値の選び方にいくつか選択肢がある



テイラー展開で整理

Δx

x

u

Δx

x

u

Δx

x

u

Δx

x

u

dx

du

Δx

)

(

)

(

)

(

)

(

lim

0





























2 2₂ 3 3₃

!

3

!

2

)

(

)

(

dx

u

d

Δx

dx

u

d

Δx

dx

du

Δx

x

u

Δx

x

u



























2 2₂ 3 3₃

_

!

3

!

2

1

)

(

)

(

dx

u

d

Δx

dx

u

d

Δx

Δx

Δx

Δx

x

u

x

u

dx

du

差分法（理論的なお話）



差分の取り方



テイラー展開で整理

Δx

Δx

x

u

Δx

x

u

Δx

Δx

x

u

Δx

x

u

dx

du

Δx

2

)

(

)

(

2

)

(

)

(

lim

0











































Δx

Δx

x

u

x

u

Δx

x

u

Δx

x

u

(

)

(

)

(

)

(

)

2

1

























3 3₃

_

!

3

2

2

1

2

)

(

)

(

dx

u

d

Δx

Δx

Δx

Δx

x

u

Δx

x

u

dx

du

差分法の概念図

u

(x)

x

x

=0 ･･･ x−

x

x x+x

x

Δx

x

u

Δx

x

u

(



)



(

)

Δx

Δx

x

u

x

u

(

)



(



)

Δx

Δx

x

u

Δx

x

u

2

)

(

)

(







中心差分

前進差分

後退差分

差分法の概念図

u

(x)

x

i

=0 ･･･ i−1 i

i

+1

x

=0 ･･･ (i−1)

x

ix

(i+1)

x

x

サンプリングされた関数値

を配列u[]で保持

差分法の概念図

u[i]

i

i=0 ･･･ i−1       i

i+1

dx

サンプリングされた関数値

を配列u[]で保持

u[i]

u[i‐1]

u[i+1]

中心差分

(u[i+1]‐u[i‐1]) 

/(2*dx)

2階微分の離散化



テイラー展開を応用



x方向に



x離れた2点で展開



_{2式を足すと1階微分が消滅}



























2 2 ₂ 3 3 ₃

!

3

!

2

)

(

)

(

x

u

Δx

x

u

Δx

x

u

Δx

x

u

Δx

x

u



























2 2 ₂ 3 3 ₃

!

3

!

2

)

(

)

(

x

u

Δx

x

u

Δx

x

u

Δx

x

u

Δx

x

u



















2 2 ₂

!

2

2

)

(

2

)

(

)

(

x

u

Δx

x

u

Δx

x

u

Δx

x

u

2階微分の離散化



_{2階微分の式に整理}

2 2 2

₍

₎

₂

₍

₎

₍

₎

Δx

Δx

x

u

x

u

Δx

x

u

x

u

_













u[i]

i

dx

サンプリングされた関数値

を配列u[]で保持

u[i]

u[i‐1]

u[i+1]

2階の中心差分

(u[i+1]‐2*u[i]+u[i‐1]) 

/(dx*dx)

差分法の実装



_{2階微分の中心差分近似}

d2udx2[i] u[i]

＋

＋

＋

＋

＋

＋

2 1 1 2 2 2

₍

₎

₂

₍

₎

₍

₎

₂

Δx

u

u

u

Δx

Δx

x

u

x

u

Δx

x

u

dx

u

d

_{i i i} x  

















2 1 Δx

＋

＋

差分法の実装



計算領域内部



_{d2udx2[i]=(u[i‐1]‐2*u[i]+u[i+1])/dxdx;}



境界条件(関数値が無いため処理を変更)



d2udx2[0  ]=(2*u[0  ]‐5*u[1  ]+4*u[2  ]‐u[3  ])/dxdx;



d2udx2[N‐1]=(2*u[N‐1]‐5*u[N‐2]+4*u[N‐3]‐u[N‐4])/dxdx;

d2udx2[i] u[i]

＋

＋

＋

＋

＋

2 1 Δx

＋

差分法の実装



計算領域内部



_{d2udx2[i]=(u[i‐1]‐2*u[i]+u[i+1])/dxdx;}



境界条件(関数値が無いため処理を変更)



d2udx2[0  ]=(2*u[0  ]‐5*u[1  ]+4*u[2  ]‐u[3  ])/dxdx;



d2udx2[N‐1]=(2*u[N‐1]‐5*u[N‐2]+4*u[N‐3]‐u[N‐4])/dxdx;

d2udx2[i] u[i]

＋

2 1 Δx

差分法の実装



計算領域内部



_{d2udx2[i]=(u[i‐1]‐2*u[i]+u[i+1])/dxdx;}



境界条件(関数値が無いため処理を変更)



d2udx2[0  ]=(2*u[0  ]‐5*u[1  ]+4*u[2  ]‐u[3  ])/dxdx;



d2udx2[N‐1]=(2*u[N‐1]‐5*u[N‐2]+4*u[N‐3]‐u[N‐4])/dxdx;

d2udx2[i] u[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 *u){ int i; for(i=0; i<Nx; i++){ u[i] = sin(i*dx); } } void differentiate(double *u, double *d2udx2){ int i; d2udx2[0] = ( 2.0*u[0] ‐5.0*u[1] +4.0*u[2] ‐ u[3])/dxdx; for(i=1; i<Nx‐1; i++) d2udx2[i] = (     u[i+1] ‐2.0*u[i ] +    u[i‐1])/dxdx; d2udx2[Nx‐1]=(‐ u[Nx‐4] +4.0*u[Nx‐3] ‐5.0*u[Nx‐2] +2.0*u[Nx‐1])/dxdx; } int main(void){ double *u,*d2udx2; u      = (double *)malloc(Nbytes); d2udx2 = (double *)malloc(Nbytes); init(u); differentiate(u,d2udx2); return 0; }

CPUプログラム

differentiate.c

関数の離散化



計算領域の長さと離散点の数，離散点の間隔の関係

u

(x)

x

x

0

2



L

_x

0からL

_x

の間に設けられた点の数

N

_x

=L

_x

/

x

+ 1

実行結果

   u, d2u /dx2

GPUへの移植



計算領域内部を計算するスレッド



_{d2udx2[i]=(u[i‐1]‐2*u[i]+u[i+1])/dxdx;}



境界を計算するスレッド



d2udx2[0  ]=(2*u[0  ]‐5*u[1  ]+4*u[2  ]‐u[3  ])/dxdx;



d2udx2[N‐1]=(2*u[N‐1]‐5*u[N‐2]+4*u[N‐3]‐u[N‐4])/dxdx;

d2udx2[i] u[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 *u){ int i; for(i=0; i<Nx; i++){ u[i] = sin(i*dx); } } __global__ void differentiate (double *u, double *d2udx2){ int i = blockIdx.x*blockDim.x + threadIdx.x; if(i==0) d2udx2[i] = ( 2.0*u[i ] ‐5.0*u[i+1] +4.0*u[i+2] ‐ u[i+3])/dxdx; if(0<i && i<Nx‐1) d2udx2[i] = (     u[i+1] ‐2.0*u[i ] +    u[i‐1])/dxdx; if(i==Nx‐1) d2udx2[i]=(‐ u[i‐3] +4.0*u[i‐2] ‐5.0*u[i‐1] +2.0*u[i ])/dxdx; }

GPUプログラム

differentiate.cu

int main(void){ double *host_u,*host_d2udx2; double *u,*d2udx2; host_u =(double *)malloc(Nbytes); cudaMalloc((void **)&u,Nbytes); cudaMalloc((void **)&d2udx2,Nbytes); init(host_u); cudaMemcpy(u, host_u, Nbytes, cudaMemcpyHostToDevice); differentiate<<<NB, NT>>>(u, d2udx2); //host_d2udx2=(double *)malloc(Nbytes); //cudaMemcpy(host_d2udx2, d2udx2, Nbytes, cudaMemcpyDeviceToHost); //for(int i=0; i<Nx; i++)printf("%f,%f,%f¥n",i*dx,host_u[i],host_d2udx2[i]); free(host_u); free(host_d2udx2); cudaFree(u); cudaFree(d2udx2); return 0; }

GPUプログラム

differentiate.cu

2次元への拡張



拡散方程式



_1次元



_2次元

2 2

₍

_,

₎

)

,

(

x

t

x

u

t

t

x

u











2 2 2 2

)

,

,

(

)

,

,

(

)

,

,

(

y

t

y

x

u

x

t

y

x

u

t

t

y

x

u

















differentiate.cuで計算

どのように2次元

に拡張するか

どのように離散化

_するか

2次元への拡張



x方向2階偏微分



y方向を固定してx方向に偏微分



y方向2階偏微分



x方向を固定してy方向に偏微分

2 2 2

₍

_,

₎

₂

₍

_,

₎

₍

_,

₎

Δx

y

Δx

x

u

y

x

u

y

Δx

x

u

x

u

_













2 2 2

₍

_,

₎

₂

₍

_,

₎

₍

_,

₎

Δy

Δy

y

x

u

y

x

u

Δy

y

x

u

y

u

_













時間積分



時間微分項の離散化



時間微分項を前進差分で離散化



右辺の

t+



tの項を移行

t

t

y

x

u

t

t

y

x

u

t

u















(

,

,

)

(

,

,

)

u

t

u

₂









t

u

t

t

y

x

u

t

t

y

x

u















)

(

,

,

)

,

,

(

拡散方程式を代入

)

,

,

(

)

,

,

(

)

,

,

(

x

y

t

t

u

x

y

t

t

2

u

x

y

t

u













u

(x, y, t)の2階微分を計算できればu(x, y, t+

t

)が求められる

離散化された方程式の記述



簡略化した表現



配列との対応をとるため下付き添字

i, jを利用



時間は，上付き添字

nを利用

j i

u

y

x

u

(

,

)



_, j i

u

y

x

x

u

(





,

)



_1, 1 ,

)

,

(

x

y





y



u

_{i j}

u

n j i

u

t

y

x

u

(

,

,

)



_, 1 ,

)

,

,

(







n j i

u

t

t

y

x

u

離散化された拡散方程式



連続系



離散系





t秒後の値

2 2 2 2

₍

_,

_,

₎

₍

_,

_,

₎

)

,

,

(

y

t

y

x

u

x

t

y

x

u

t

t

y

x

u

















2 1 , , 1 , 2 , 1 , , 1 , 1 ,

2

2

y

u

u

u

x

u

u

u

t

u

u

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

u

u

u

x

u

u

u

t

u

u

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

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

u

u

u

x

u

u

u

t

u

u

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

₂



0

.

5





y

t

v





































     2 1 , , 1 , 2 , 1 , , 1 , 1 ,

2

2

y

u

u

u

x

u

u

u

t

u

u

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.

初期値の決定



uの初期分布の決定

3.

差分値の計算



uの分布からx, y方向の2階微分値を計算



境界条件に基づいて境界の値を決定





t秒後のuを計算

CPUプログラム



計算条件の決定



計算領域の大きさ

L

_x

, L

_y



計算領域の分割数

（離散点の個数）

N

_x

, N

_y



離散点同士の間隔

（格子間隔）



_x,



_y



計算時間間隔



t

#include<stdio.h> #include<stdlib.h> #include<math.h> #define Lx (2.0*M_PI) #define Ly (2.0*M_PI) #define Nx 512 #define Ny 512 #define dx (Lx/(Nx‐1)) #define dy (Ly/(Ny‐1)) #define dt 0.001

CPUプログラム



初期条件



境界条件

for(j=0;j<Ny;j++){ for(i=0;i<Nx;i++){ x = (double)i*dx; y = (double)j*dy; u[i*Ny+j]=sin(x)*sin(y); } }

y

x

u



sin

sin

2

2

0

sin

sin

2

2

_

_

_



u

x

y

#include<stdio.h> #include<stdlib.h> #include<math.h> #define Lx (2.0*M_PI) #define Ly (2.0*M_PI) #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 (1.0) #define dxdx (dx*dx) #define dydy (dy*dy) #define Nbytes (Nx*Ny*sizeof(double)) void laplacian(double *, double *); void integrate(double *, double *, double *); void update(double *, double *); int main(void){ double *u,*unew,*ulap,x,y; int i,j,ij,ip1,im1,jp1,jm1,n; u    = (double *)malloc(Nx*Ny*sizeof(double)); unew = (double *)malloc(Nx*Ny*sizeof(double)); ulap = (double *)malloc(Nx*Ny*sizeof(double)); for(j=0;j<Ny;j++){ for(i=0;i<Nx;i++){ x = (double)i*dx; y = (double)j*dy; u[i*Ny+j]=sin(x)*sin(y); unew[i*Ny+j]=0.0f; ulap[i*Ny+j]=0.0f; } } for(n=0;n<Nt;n++){ laplacian(u,ulap); integrate(u,ulap,unew); update(u,unew); } return 0; }

CPUプログラム

diffusion.c

void laplacian(double *u, double *ulap){ int i,j,ij,ip1,im1,jp1,jm1; for(j=1;j<Ny‐1;j++){ for(i=1;i<Nx‐1;i++){ ij =  i *Ny+j; ip1 = (i+1)*Ny+j; im1 = (i‐1)*Ny+j; jp1 =  i *Ny+j+1; jm1 =  i *Ny+j‐1; ulap[ij] = (u[ip1]‐2.0*u[ij]+u[im1])/dxdx +(u[jp1]‐2.0*u[ij]+u[jm1])/dydy; } } } void integrate (double *u, double *ulap, double *unew){ int i,j,ij; for(j=0;j<Ny;j++){ for(i=0;i<Nx;i++){ ij =  i*Ny+j; unew[ij] = u[ij] + dt*DIFF*ulap[ij]; } } } void update(double *u, double *unew){ int i,j,ij; for(j=0;j<Ny;j++){ for(i=0;i<Nx;i++){ ij =  i*Ny+j; u[ij] = unew[ij]; } } }

CPUプログラム

diffusion.c

CPUプログラム



差分計算



x方向，y方向偏微分を個別に計算して加算

void laplacian(double *u, double *ulap){ int i,j,ij,ip1,im1,jp1,jm1; for(j=1;j<Ny‐1;j++){ for(i=1;i<Nx‐1;i++){ ij =  i *Ny+j; ip1 = (i+1)*Ny+j; im1 = (i‐1)*Ny+j; jp1 =  i *Ny+j+1; jm1 =  i *Ny+j‐1; ulap[ij] = (u[ip1]‐2.0*u[ij]+u[im1])/dxdx +(u[jp1]‐2.0*u[ij]+u[jm1])/dydy; } } }

ラプラシアン計算のメモリ参照



ある1点i,jのラプラシアンを計算するために，周囲

5点のuを参照

u[]

ulap[]

i

−1 i i+1

j−

1

jj

+1

ラプラシアン計算のメモリ参照



ある1点i,jのラプラシアンを計算するために，周囲

5点のuを参照

u[]

ulap[]

i

−1 i i+1

j−

1

jj

+1

ラプラシアン計算のメモリ参照



ある1点i,jのラプラシアンを計算するために，周囲

5点のuを参照

u[]

ulap[]

i

−1 i i+1

j−

1

jj

+1

ラプラシアン計算のメモリ参照



ある1点i,jのラプラシアンを計算するために，周囲

5点のuを参照

u[]

ulap[]

i

−1 i i+1

j−

1

jj

+1

ラプラシアン計算のメモリ参照



ある1点i,jのラプラシアンを計算するために，周囲

5点のuを参照



全てのuを参照し，領域内部のulapを計算

CPUプログラム



境界条件



uのラプラシアン



境界ではどの時刻においても

0



x=0, 0≤y≤2





x=2

, 0≤y≤2



0≤x≤2, y=0



0≤x≤2, y=2



変数ulapを0で初期化すれば，計算しなくてもよい

y

x

u

2

sin

sin

2







ulap[]

CPUプログラム



uの積分

void integrate (double *u, double *ulap, double *unew){ int i,j,ij; for(j=0;j<Ny;j++){ for(i=0;i<Nx;i++){ ij =  i*Ny+j; unew[ij] = u[ij] + dt*DIFF*ulap[ij]; } } }





































     2 1 , , 1 , 2 , 1 , , 1 , 1 ,

2

2

y

u

u

u

x

u

u

u

t

u

u

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プログラム



uの更新



u

n

から

u

n+1

を計算



u

n+1

から

u

n+2

を計算



今の時刻から次の時刻を求める



求められた次の時刻を今の時刻と見なし，次の時刻を求める

void update(double *u, double *unew){ int i,j,ij; for(j=0;j<Ny;j++){ for(i=0;i<Nx;i++){ ij =  i*Ny+j; u[ij] = unew[ij]; } } }

同じアルゴリズム

#include<stdio.h> #include<stdlib.h> #include<math.h> #define Lx (2.0*M_PI) #define Ly (2.0*M_PI) #define Nx 128 #define Ny Nx #define dx (Lx/(Nx‐1)) #define dy (Ly/(Ny‐1)) #define dt 0.00001 #define endT (1.0) #define Nt (int)(endT/dt) #define DIFF (1.0) #define dxdx (dx*dx) #define dydy (dy*dy) #define Nbytes (Nx*Ny*sizeof(double)) #include "dif1.cu" //#include "dif2.cu" //#include "dif3.cu" int main(void){ int n;     double *dev_u,*dev_unew,*dev_ulap; dim3 Thread, Block; if(DIFF*dt/dxdx > 0.5){ printf("configuration error¥n"); exit(1); } cudaMalloc( (void**)&dev_u , Nbytes ); cudaMalloc( (void**)&dev_unew, Nbytes ); cudaMalloc( (void**)&dev_ulap, Nbytes ); Thread = dim3(THREADX,THREADY,1); Block  = dim3(BLOCKX ,BLOCKY, 1); init<<<Block, Thread>>> (dev_u, dev_ulap, dev_unew); for(n=0;n<Nt;n++){ laplacian<<<Block, Thread>>> (dev_u,dev_ulap); integrate<<<Block, Thread>>> (dev_u,dev_ulap,dev_unew); update<<<Block, Thread>>>(dev_u,dev_unew); } return 0; }

GPUプログラム（CPU処理用共通部分）

diff.cu

#define THREADX 1 #define THREADY 1 #define BLOCKX  1 #define BLOCKY  1 __global__ void init (double *u, double *ulap, double *unew){ int i,j,ij; double x,y; for(j=0;j<Ny;j++){ for(i=0;i<Nx;i++){ ij =  i*Ny+j; x = (double)i*dx; y = (double)j*dy; u[ij]=sin(x)*sin(y); unew[ij]=0.0; ulap[ij]=0.0; } }   } __global__ void laplacian(double *u, double *ulap){ int i,j,ij,ip1,im1,jp1,jm1; for(j=1;j<Ny‐1;j++){ for(i=1;i<Nx‐1;i++){ ij =  i *Ny+j; ip1 = (i+1)*Ny+j; im1 = (i‐1)*Ny+j; jp1 =  i *Ny+j+1; jm1 =  i *Ny+j‐1; ulap[ij] = (u[ip1]‐2.0f*u[ij]+u[im1])/dxdx +(u[jp1]‐2.0f*u[ij]+u[jm1])/dydy; } } } __global__ void integrate (double *u, double *ulap, double *unew){ int i,j,ij; for(j=0;j<Ny;j++){ for(i=0;i<Nx;i++){ ij =  i *Ny+j; unew[ij] = u[ij] + dt*DIFF*ulap[ij]; } } } __global__ void update(double *u, double *unew){ int i,j,ij; for(j=0;j<Ny;j++){ for(i=0;i<Nx;i++){ ij =  i *Ny+j; u[ij] = unew[ij]; } } }

GPUプログラム（1スレッド版）

dif1.cu

2次元ブロック分割



_{1スレッドが1点のラプラシアンを計算}

blockIdx.x=0 blockIdx.x=1

blockIdx.y=0

blockIdx.y=1

gridDim.x=2

gridDim.y=2

blockDim.x=4

blockDim.y=4

threadIdx.x=

threadIdx.y=

2次元ブロック分割



_{Nx=8, Ny=8, x,y方向スレッド数4，ブロック数2}



_{i =} 

_blockIdx.x

_*

_blockDim.x

_{+ threadIdx.x}



_j = 

_blockIdx.y

_*

_blockDim.y

_{+ threadIdx.y}

(0,0)(1,0)(2,0)(3,0)(0,0) (3,3) (3,3) (0,1)(1,1)(2,1)(3,1) (0,2)(1,2)(2,2)(3,2) (0,3)(1,3)(2,3)(3,3) (3,3) (0,0) (0,0)

block(0,0)

block(1,0)

threadIdx.x

threadIdx.y

i= 0 1

2 3

4

5 6

7

j=

0

1

2

3

4

5

6

7

2次元的な配列アクセスの優先方向(CPU)



_{2次元配列の1次元配列的表現}

for(i=0;i<Nx;i++)

for(j=0;j<Ny;j++)

out[i][j]=in[i][j];

in[],out[]

Nx

Ny

j

i

for(i=0;i<Nx;i++)

for(j=0;j<Ny;j++)

out[i*Ny+j]=

in[i*Ny+j];

2次元的な配列アクセスの優先方向(GPU)



_{CUDAで2次元的に並列化してアクセスする場合}

i = blockIdx.x*blockDim.x

+ threadIdx.x;

j = blockIdx.y*blockDim.y

+ threadIdx.y;

out[

j*Nx+i

]=in[

j*Nx+i

];

in[],out[]

Nx

Ny

j

i

for(i=0;i<Nx;i++)

for(j=0;j<Ny;j++)

out[i*Ny+j]=

in[i*Ny+j];

threadIdx.x threadIdx.y

#define THREADX 16 #define THREADY 16 #define BLOCKX  (Nx/THREADX) #define BLOCKY  (Ny/THREADY) __global__ void laplacian(double *u, double *ulap){ int i,j,ij,ip1,im1,jp1,jm1; i = blockIdx.x*blockDim.x + threadIdx.x; j = blockIdx.y*blockDim.y + threadIdx.y; if( (0<i && i<Nx‐1) && (0<j && j<Ny‐1) ){ ij = i +Nx*j; ip1 = i+1+Nx*j; im1 = i‐1+Nx*j; jp1 = i +Nx*(j+1); jm1 = i +Nx*(j‐1); ulap[ij] = (u[ip1]‐2.0*u[ij]+u[im1])/dxdx +(u[jp1]‐2.0*u[ij]+u[jm1])/dydy; } } __global__ void integrate (double *u, double *ulap, double *unew){ int i,j,ij; i = blockIdx.x*blockDim.x + threadIdx.x; j = blockIdx.y*blockDim.y + threadIdx.y; ij = i+Nx*j; unew[ij] = u[ij] + dt*ulap[ij]; } __global__ void update(double *u, double *unew){ int i,j,ij; i = blockIdx.x*blockDim.x + threadIdx.x; j = blockIdx.y*blockDim.y + threadIdx.y; ij = i+Nx*j; u[ij] = unew[ij]; } __global__ void init (double *u, double *ulap, double *unew){ 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; x = (double)i*dx; y = (double)j*dy; u[ij]=sin(x)*sin(y); unew[ij]=0.0; ulap[ij]=0.0; }

GPUプログラム（1スレッドが1点を計算）

dif2.cu

境界でラプラシアンが

0にならない場合



共有メモリを利用してキャッシュを模擬



共有メモリに付加的な領域を追加

境界でラプラシアンが

0にならない場合



共有メモリを利用してキャッシュを模擬



共有メモリに付加的な領域を追加

境界でラプラシアンが

0にならない場合



共有メモリを利用してキャッシュを模擬



共有メモリに付加的な領域を追加

u[]

ulap[]

su[][]

付加的な領域（袖領域）の取り扱い



データがグローバルメモリに存在する場合は，グ

ローバルメモリから読み込み

u[]

ulap[]

su[][]

付加的な領域（袖領域）の取り扱い



グローバルメモリに無い場合は境界条件から決定

u[]

ulap[]

su[][]

境界条件



_{2階微分が0}

0

2

2 1 1 2 2













_{ }

Δx

u

u

u

x

u

_{i i i} 1 1

2

i i_ i

u

u

u

_





1 1

2

j j_ j

u

u

u

_





0

2

2 1 1 2 2













 

Δy

u

u

u

y

u

j j j

境界条件



_{2階微分を片側差分で計算}

3 2 1 1



4

j



6

j



4

j



j j

u

u

u

u

u

_ 2 3 2 1 1 2 2

₂

₅

₄

Δy

u

u

u

u

y

u

j j j j j   













 2 3 2 1 1 2 2

₂

₅

₄

Δx

u

u

u

u

x

u

_{i i i i} i   













 3 2 1 1



4

i



6

i



4

i



i i

u

u

u

u

u

_

#define THREADX 16 #define THREADY 16 #define BLOCKX  (Nx/THREADX) #define BLOCKY  (Ny/THREADY) __global__ void laplacian(double *u, double *ulap){ int i,j,ij; int tx,ty; __shared__ float su[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; su[tx][ty] = u[ij]; __syncthreads(); if(   blockIdx.x == 0 && threadIdx.x == 0       ) su[tx‐1][ty] = 2.0*su[tx][ty] ‐ su[tx+1][ty]; if(    blockIdx.x != 0 && threadIdx.x == 0       )  su[tx‐1][ty] = u[j*Nx+i‐1]; if(   blockIdx.x == gridDim.x‐1 && threadIdx.x == blockDim.x‐1) su[tx+1][ty] = 2.0*su[tx][ty] ‐ su[tx‐1][ty]; if(   blockIdx.x != gridDim.x‐1 && threadIdx.x == blockDim.x‐1) su[tx+1][ty] = u[j*Nx+i+1]; if(   blockIdx.y == 0 && threadIdx.y == 0) su[tx][ty‐1] = 2.0*su[tx][ty] ‐ su[tx][ty+1]; if(   blockIdx.y != 0 && threadIdx.y == 0)       su[tx][ty‐1] = u[(j‐1)*Nx+i]; if(   blockIdx.y == gridDim.y‐1 && threadIdx.y == blockDim.y‐1) su[tx][ty+1] = 2.0*su[tx][ty] ‐ su[tx][ty‐1]; if(   blockIdx.y != gridDim.y‐1 && threadIdx.y == blockDim.y‐1) su[tx][ty+1] = u[(j+1)*Nx+i]; __syncthreads(); ulap[ij] =  (su[tx‐1][ty  ]‐2.0*su[tx][ty]+su[tx+1][ty  ])/dxdx +(su[tx  ][ty‐1]‐2.0*su[tx][ty]+su[tx  ][ty+1])/dydy; }

GPUプログラム（ラプラシアンkernel）

dif3.cu

GPUプログラム（ラプラシアンkernel）



ブロック内のスレッド数＋袖領域分の共有メモ

リを確保



袖領域があるために添字の対応が変化

 添字の対応を考えないと，必要なデータを袖領域に置い てしまう 

_{__syncthreads()を呼んでスレッドを同期}

共有メモリにデータが正しく書き込まれた事を保証

__shared__ float su[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+Nx*j; su[tx][ty] = u[ij]; __syncthreads();

u[]

su[][]

GPUプログラム（ラプラシアンkernel）



ブロック内のスレッド数＋袖領域分の共有メモ

リを確保



袖領域があるために添字の対応が変化

 添字の対応を考えないと，必要なデータを袖領域に置い てしまう 

_{__syncthreads()を呼んでスレッドを同期}

共有メモリにデータが正しく書き込まれた事を保証

__shared__ float su[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+Nx*j; su[tx][ty] = u[ij]; __syncthreads();

u[]

su[][]

GPUプログラム（ラプラシアンkernel）



ブロック内のスレッド数＋袖領域分の共有メモ

リを確保



袖領域があるために添字の対応が変化

 添字の対応を考えないと，必要なデータを袖領域に置い てしまう 

_{__syncthreads()を呼んでスレッドを同期}

共有メモリにデータが正しく書き込まれた事を保証

__shared__ float su[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; su[tx][ty] = u[ij]; __syncthreads();

u[]

su[][]

+1

GPUプログラム（ラプラシアンkernel）



袖領域の設定



グローバルメモリにデータがある箇所はグロー

バルメモリから読み込み



グローバルメモリにデータがない箇所は2階微

分が0になるように袖領域の値を決定

 ブロックが境界に接しているか否かで処理を切替

if(blockIdx.x == 0       && threadIdx.x == 0       ) su[tx‐1][ty] = 2.0*su[tx][ty] ‐ su[tx+1][ty]; if(blockIdx.x != 0       && threadIdx.x == 0       ) su[tx‐1][ty] = u[j*Nx+i‐1];

if(blockIdx.x == gridDim.x‐1 && threadIdx.x == blockDim.x‐1) su[tx+1][ty] = 2.0*su[tx][ty] ‐ su[tx‐1][ty]; if(blockIdx.x != gridDim.x‐1 && threadIdx.x == blockDim.x‐1) su[tx+1][ty] = u[j*Nx+i+1];

if(blockIdx.y == 0       && threadIdx.y == 0       ) su[tx][ty‐1] = 2.0*su[tx][ty] ‐ su[tx][ty+1]; if(blockIdx.y != 0       && threadIdx.y == 0       ) su[tx][ty‐1] = u[(j‐1)*Nx+i];

if(blockIdx.y == gridDim.y‐1 && threadIdx.y == blockDim.y‐1) su[tx][ty+1] = 2.0*su[tx][ty] ‐ su[tx][ty‐1]; if(blockIdx.y != gridDim.y‐1 && threadIdx.y == blockDim.y‐1) su[tx][ty+1] = u[(j+1)*Nx+i];

GPUプログラム（ラプラシアンkernel）



共有メモリのデータを利用してラプラシアンを計算



_{if分岐を排除}

ulap[]

su[][]

ulap[ij] = (su[tx‐1][ty  ]‐2.0*su[tx][ty]+su[tx+1][ty  ])/dxdx +(su[tx  ][ty‐1]‐2.0*su[tx][ty]+su[tx  ][ty+1])/dydy;

GPUプログラムの評価



共有メモリを使用したdif3.cu



_{dif2.cu（共有メモリ不使用）より速いこともあれば遅いこ}

ともある



_{Fermi世代以降のGPUはキャッシュを搭載}



共有メモリを使っても速くならない



共有メモリへ明示的にデータを移動



余分な処理により負荷が増加

__global__ void update(double *u, double *unew){

int i,j,ij;

i = blockIdx.x*blockDim.x + threadIdx.x;

j = blockIdx.y*blockDim.y + threadIdx.y;

ij = i+Nx*j;

u[ij] = unew[ij];

}

その他の処理の高速化



値の更新



_{unewのデータをuにコピーしているだけ}



_{cudaMemcpyで代用可能}

その他の処理の高速化



値の更新



_{unewのデータをuにコピーしているだけ}



_{cudaMemcpyで代用可能}

for(n=0;n<Nt;n++){

laplacian<<<Block, Thread>>>(dev_u,dev_ulap);

integrate<<<Block, Thread>>>(dev_u,dev_ulap,dev_unew);

//update<<<Block, Thread>>>(dev_u,dev_unew);

cudaMemcpy(dev_u, dev_new, Nbytes, cudaMemcpyDeviceToDevice);

}

return 0;

その他の処理の高速化



初期値の計算



三角関数がそれらしい値で求められればいい場合



有効数字4桁程度

__global__ void init(double *u, double *ulap, double *unew){

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;

x = (double)i*dx;y = (double)j*dy;

u[ij]=sin(x)*sin(y);

unew[ij]=0.0;

ulap[ij]=0.0;

}

その他の処理の高速化



初期値の計算



三角関数がそれらしい値で求められればいい場合



有効数字4桁程度

__global__ void init(double *u, double *ulap, double *unew){

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;

x = (double)i*dx;y = (double)j*dy;

u[ij]=__sinf(x)*__sinf(y);

unew[ij]=0.0;

ulap[ij]=0.0;

}

FastMath関数



三角関数など数学関数を高速に計算



精度は落ちる



手動で書き換え



_{__sinf}



_{__cosf}



_{__powfなど}



_{‐use_fast_mathオプション}



_{sin,cosを利用していても，‐use‐fast‐mathオプションを}

付けることで__sinf, __cosfに置き換えられる

付録

共有メモリの典型的な使い方

共有メモリの典型的な使い方



差分法



あるスレッドが中心差分を計算するために配列の要素i‐1, 

i, i+1を参照



配列要素は複数回参照される



_{Fermi世代以降はキャッシュが利用可能}

1

1

2

2

1

1

dudx[i] u[i]

＋

＋

＋

＋

Δx 2 1

参照される回数

共有メモリの典型的な使い方



境界で異なる処理を行うためにif分岐が必要



キャッシュを搭載していてもグローバルメモリへのアクセス

を伴うif分岐は高負荷



データの再利用とif文の排除に共有メモリを利用

2

2

3

3

2

2

dudx[i] u[i]

＋

＋

＋

＋

＋

＋

Δx 2 1

参照される回数

共有メモリによる明示的なキャッシュ



グローバルメモリから共有メモリにデータをキャッシュ



共有メモリ上で境界条件を処理



中心差分の計算からifを排除

dudx[i] u[i] 共有メモリ

blockIdx.x=0

blockIdx.x=1

計算のために必要になる余分な領域（袖領域）

共有メモリによる明示的なキャッシュ



グローバルメモリから共有メモリにデータをキャッシュ



共有メモリ上で境界条件を処理



中心差分の計算からifを排除

dudx[i] u[i] 共有メモリ

blockIdx.x=0

blockIdx.x=1

何らかの方法で

境界条件を反映

何らかの方法で

境界条件を反映

共有メモリによる明示的なキャッシュ



グローバルメモリから共有メモリにデータをキャッシュ



共有メモリ上で境界条件を処理



中心差分の計算からifを排除

dudx[i] u[i] 共有メモリ ＋ ＋ ＋ ＋ ＋ ＋

blockIdx.x=0

blockIdx.x=1

全スレッドが同じ式

で中心差分を計算

#include<stdlib.h> #include<math.h>/*‐lmオプションが必要*/ #define Lx (2.0*M_PI) #define Nx (1024*1024) #define dx (Lx/(Nx‐1)) #define Nbytes (Nx*sizeof(double)) #define NT (256) #define NB (Nx/NT) void init(double *u){ int i; for(i=0; i<Nx; i++){ u[i] = sin(i*dx); } } __global__ void differentiate (double *u, double *dudx){ int i = blockIdx.x*blockDim.x + threadIdx.x; if(i==0) dudx[i]=(‐3.0*u[i ] +4.0*u[i+1] ‐ u[i+2])/(2.0*dx); if(0<i && i<Nx‐1) dudx[i] = ( u[i+1] ‐u[i‐1])/(2.0*dx); if(i==Nx‐1) dudx[i]=(     u[i‐2] ‐4.0*u[i‐1] +3.0*u[i ])/(2.0*dx); }

GPUプログラム（単純な実装）

differentiate1d.cu

int main(void){ double *host_u,*host_dudx; double *u,*dudx; host_u =(double *)malloc(Nbytes); host_dudx=(double *)malloc(Nbytes); cudaMalloc((void **)&u,Nbytes); cudaMalloc((void **)&dudx,Nbytes); init(host_u); cudaMemcpy(u,host_u,Nbytes,  cudaMemcpyHostToDevice); differentiate<<<NB, NT>>>(u,dudx); cudaMemcpy(host_dudx, dudx, Nbytes,  cudaMemcpyDeviceToHost); free(host_u); free(host_dudx); cudaFree(u); cudaFree(dudx); return 0; }

GPUプログラム（単純な実装）

differentiate1d.cu

__global__ void differentiate(double *u, double *dudx){ int i = blockIdx.x*blockDim.x + threadIdx.x;

__shared__ double su[1+NT+1]; int tx = threadIdx.x+1; su[tx] = u[i]; __syncthreads(); if(blockIdx.x> 0       && threadIdx.x==0       ) su[tx‐1] = u[i‐1]; if(blockIdx.x< gridDim.x‐1 && threadIdx.x==blockDim.x‐1) su[tx+1] = u[i+1]; if(blockIdx.x==0       && threadIdx.x==0       ) su[tx‐1] = 3.0*su[tx]‐3.0*su[tx+1]+su[tx+2]; if(blockIdx.x==gridDim.x‐1 && threadIdx.x==blockDim.x‐1) su[tx+1] = 3.0*su[tx]‐3.0*su[tx‐1]+su[tx‐2]; __syncthreads(); dudx[i] = ( su[tx+1]‐su[tx‐1])/(2.0*dx); }

共有メモリを用いた書き換え

differentiate1d_shared.cu

共有メモリの宣言と代入

int i = blockIdx.x*blockDim.x + threadIdx.x;

__shared__ double su[1+NT+1]; //右と左の袖領域を追加して宣言 int tx = threadIdx.x+1; su[tx] = u[i]; __syncthreads(); u[i] su[NT+2] i=   0    1    2      3    4    5 threadIdx.x= 0    1    2          0    1    2 tx=   0    1    2    3   4       0    1    2    3    4 NT=3 NT=3 u[0] u[1] u[2]       u[3] u[4] u[5]

袖領域の処理

if(blockIdx.x> 0       && threadIdx.x==0       ) su[tx‐1] = u[i‐1]; if(blockIdx.x< gridDim.x‐1 && threadIdx.x==blockDim.x‐1) su[tx+1] = u[i+1]; u[i] su[NT+2] i=   0    1    2      3    4    5 0    1    2=threadIdx.x tx=   0    1    2    3   4       0    1    2    3    4 blockIdx.x=0 blockIdx.x=1 u[0] u[1] u[2]         u[3] u[4] u[5]

袖領域の処理

if(blockIdx.x> 0       && threadIdx.x==0       ) su[tx‐1] = u[i‐1]; if(blockIdx.x< gridDim.x‐1 && threadIdx.x==blockDim.x‐1) su[tx+1] = u[i+1]; u[i] su[NT+2] i=   0    1    2      3    4    5 0=threadIdx.x tx=   0    1    2    3   4       0    1    2    3    4 blockIdx.x=0 blockIdx.x=1

袖領域の処理

if(blockIdx.x> 0       && threadIdx.x==0       ) su[tx‐1] = u[i‐1]; if(blockIdx.x< gridDim.x‐1 && threadIdx.x==blockDim.x‐1) su[tx+1] = u[i+1]; u[i] su[NT+2] i=   0    1    2      3    4    5 tx=   0    1    2    3   4       0    1    2    3    4 blockIdx.x=0 blockIdx.x=1

u[0] u[1] u[2] u[2] u[3] u[4] u[5] threadIdx.x= 0    1    2

袖領域の処理

if(blockIdx.x> 0       && threadIdx.x==0       ) su[tx‐1] = u[i‐1]; if(blockIdx.x< gridDim.x‐1 && threadIdx.x==blockDim.x‐1) su[tx+1] = u[i+1]; u[i] su[NT+2] i=   0    1    2      3    4    5 tx=   0    1    2    3   4       0    1    2    3    4 blockIdx.x=0 blockIdx.x=1

u[0] u[1] u[2] u[3] u[2] u[3] u[4] u[5] threadIdx.x=2

境界条件の処理

if(blockIdx.x==0       && threadIdx.x==0       ) su[tx‐1] = 3.0*su[tx]‐3.0*su[tx+1]+su[tx+2]; if(blockIdx.x==gridDim.x‐1 && threadIdx.x==blockDim.x‐1) su[tx+1] = 3.0*su[tx]‐3.0*su[tx‐1]+su[tx‐2]; su[NT+2] tx=   0    1    2    3   4       0    1    2    3    4 blockIdx.x=0 blockIdx.x=1

u[0] u[1] u[2] u[3] u[2] u[3] u[4] u[5]

2 1 0 1

3

u

3

u

u

u

_







Δx

u

u

Δx

u

u

u

dx

du

i

2

2

4

3

_{0 1 2 1 1} 0  













境界での差分式と中心差分式が一致するように

u

₋₁

を決定

u[‐1]

境界条件の処理

if(blockIdx.x==0       && threadIdx.x==0       ) su[tx‐1] = 3.0*su[tx]‐3.0*su[tx+1]+su[tx+2]; if(blockIdx.x==gridDim.x‐1 && threadIdx.x==blockDim.x‐1) su[tx+1] = 3.0*su[tx]‐3.0*su[tx‐1]+su[tx‐2]; su[NT+2] tx=   0    1    2    3   4       0    1    2    3    4 blockIdx.x=0 blockIdx.x=1

u[0] u[1] u[2] u[3] u[2] u[3] u[4] u[5]

3 2 1

3

3

_



_



_



_{N N N} N

u

u

u

u

Δx

u

u

Δx

u

u

u

dx

du

_{N N N N N} N i

2

2

4

3

_{1 2 3 2} 1      











境界での差分式と中心差分式が一致するように

u

_N

を決定

u[‐1] u[6]