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

偏微分方程式の差分計算

（移流方程式）

今日の内容



差分法



_{1次関数の差分}



共有メモリの利用



_{2次元移流方程式}



_{gnuplotによる結果の表示}



ダブルバッファリング

数値計算



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



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



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

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



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



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



支配方程式



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



流体の支配方程式

0













i i

x

u

t





j ij i j i j i

x

x

p

x

u

u

t

u





























i i i ij i i i i i

x

q

x

u

x

p

u

x

E

u

t

E



































（質量保存式）

（エネルギ保存式）

（運動量保存式）

移流方程式



流体中の物質の移動を表す方程式



流れている水の中に落ちたインクの移動等



時刻

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



時間進行に伴い，

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



時間積分しながら

fの分布を求める

0

)

,

(

)

,

(













x

t

x

f

u

t

t

x

f

0

)

,

,

(

)

,

,

(

)

,

,

(



















y

t

y

x

f

v

x

t

y

x

f

u

t

t

y

x

f

_{u : x方向速度}

v : y方向速度

差分法



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



微分の定義



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 ₂ 3 3 ₃

!

3

!

2

)

(

)

(

dx

f

d

Δx

dx

f

d

Δx

dx

df

Δx

x

f

Δx

x

f































2 2 ₂ 3 3 ₃

_

!

3

!

2

1

)

(

)

(

dx

f

d

Δx

dx

f

d

Δx

Δx

Δx

x

f

Δx

x

f

dx

df

差分法（理論的なお話）



空間打ち切り誤差



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



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































2 2 ₂ 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 3 ₃

_

!

3

!

2

dx

f

d

Δx

dx

f

d

Δx

差分法（理論的なお話）



誤差の主要項



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





xは小さい→



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



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

O(



x)





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

)

(

!

3

!

2

3 3 2 2 2

Δx

O

dx

f

d

Δx

dx

f

d

Δx

_

























_

差分法（理論的なお話）



差分の取り方



テイラー展開で整理

Δx

Δx

x

f

Δx

x

f

Δx

Δx

x

f

Δx

x

f

dx

df

Δx

2

)

(

)

(

2

)

(

)

(

lim

0











































Δx

Δx

x

f

x

f

Δx

x

f

Δx

x

f

(

)

(

)

(

)

(

)

2

1





























3 3 ₃

_

!

3

2

2

1

2

)

(

)

(

dx

f

d

Δx

Δx

Δx

Δx

x

f

Δx

x

f

dx

df

差分法の概念図

f

(x)

x

x

=0 ･･･ x−



x

x x+



x



x

Δx

x

f

Δx

x

f

(



)



(

)

Δx

Δx

x

f

x

f

(

)



(



)

Δx

Δx

x

f

Δx

x

f

2

)

(

)

(







中心差分

前進差分

後退差分

差分法の概念図

f

(x)

x

i

=0 ･･･ i−1 i

i

+1

x

=0 ･･･ (i−1)



x

i



x

(i+1)



x



x

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

を配列f[]で保持

差分法の概念図

f[i]

i

i

=0 ･･･ i−1 i

i

+1

dx

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

を配列f[]で保持

f[i]

f[i‐1]

f[i+1]

中心差分

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

/(2*dx)

差分法の実装



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

dfdx[i] f[i]

＋

＋

＋

＋

＋

＋

Δx

f

f

Δx

Δx

x

f

Δx

x

f

dx

df

_{i i} x

2

2

)

(

)

(







_

_1



_1



Δx 2 1

差分法の実装



計算領域内部



_{dfdx[i]=(f[i+1]‐f[i‐1])/(2.0*dx);}

dfdx[i] f[i]

＋

＋

＋

＋

Δx 2 1 ×−1

差分法の実装



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



_{dfdx[0  ]=(‐3*f[0  ]+4*f[0+1]‐f[0+2])/(2*dx);}



_{dfdx[N‐1]=( 3*f[N‐1]‐4*f[N‐2]+f[N‐3])/(2*dx);}



_{2階微分値が一定と仮定して関数を補外した事に相当}

dfdx[i] f[i]

＋

Δx 2 1 ×−3 ×4 ×−1

差分法の実装



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



_{dfdx[0  ]=(‐3*f[0  ]+4*f[0+1]‐f[0+2])/(2*dx);}



_{dfdx[N‐1]=( 3*f[N‐1]‐4*f[N‐2]+f[N‐3])/(2*dx);}



_{2階微分値が一定と仮定して関数を補外した事に相当}

dfdx[i] f[i]

＋

Δx 2 1 ×−4 ×3

差分法の実装



境界条件(境界で微分値が0と規定)



_{dfdx[0  ]=0.0;}



_{dfdx[N‐1]=0.0;}



そのまま値を代入すればよい

0.0

0.0

dfdx[i] f[i]

＋

Δx 2 1 ×−4 ×3

＋

Δx 2 1 ×−3 ×4 ×−1

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

CPUプログラム

differentiate.c

関数の離散化



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

f

(x)

x



x

0

2



L

_x

0からL

_x

の間に設けられた点の数

N

_x

=L

_x

/



x

+ 1

実行結果

   f, df /dx

GPUへの移植



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



_{dfdx[i]=(f[i+1]‐f[i‐1])/(2.0*dx);}



境界を計算するスレッド



_{dfdx[0  ]=(‐3*f[0  ]+4*f[0+1]‐f[0+2])/(2*dx);}



_{dfdx[N‐1]=( 3*f[N‐1]‐4*f[N‐2]+f[N‐3])/(2*dx);}

dfdx[i] f[i]

＋

＋

＋

＋

＋

＋

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

GPUプログラム

differentiate.cu

2階微分からの変更箇所

int main(void){ double *host_f,*host_dfdx; double *f,*dfdx; host_f =(double *)malloc(Nbytes); host_dfdx=(double *)malloc(Nbytes); cudaMalloc((void **)&f,Nbytes); cudaMalloc((void **)&dfdx,Nbytes); init(host_f); cudaMemcpy(f, host_f, Nbytes, cudaMemcpyHostToDevice); differentiate<<<NB, NT>>>(f,dfdx); cudaMemcpy(host_dfdx, dfdx, Nbytes, cudaMemcpyDeviceToHost); free(host_f); free(host_dfdx); cudaFree(f); cudaFree(dfdx); return 0; }

GPUプログラム

differentiate.cu

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



差分法



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

i, i+1を参照



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



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

1

1

2

2

1

1

dfdx[i] f[i]

＋

＋

＋

＋

Δx 2 1

参照される回数

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



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



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

を伴うif分岐は高負荷



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

2

2

3

3

2

2

dfdx[i] f[i]

＋

＋

＋

＋

＋

＋

Δx 2 1

参照される回数

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



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



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



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

dfdx[i] f[i] 共有メモリ

blockIdx.x=0

blockIdx.x=1

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

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



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



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



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

dfdx[i] f[i] 共有メモリ

blockIdx.x=0

blockIdx.x=1

何らかの方法で

境界条件を反映

何らかの方法で

境界条件を反映

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



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



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



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

dfdx[i] f[i] 共有メモリ ＋ ＋ ＋ ＋ ＋ ＋

blockIdx.x=0

blockIdx.x=1

全スレッドが同じ式

で中心差分を計算

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

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

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

differentiate_shared.cu

共有メモリの宣言と代入

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

__shared__ double sf[1+NT+1]; //右と左の袖領域を追加して宣言 int tx = threadIdx.x+1; sf[tx] = f[i]; __syncthreads(); f[i] sf[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 f[0] f[1] f[2]       f[3] f[4] f[5]

袖領域の処理

if(blockIdx.x> 0       && threadIdx.x==0       ) sf[tx‐1] = f[i‐1]; if(blockIdx.x< gridDim.x‐1 && threadIdx.x==blockDim.x‐1) sf[tx+1] = f[i+1]; f[i] sf[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 f[0] f[1] f[2]         f[3] f[4] f[5]

袖領域の処理

if(blockIdx.x> 0       && threadIdx.x==0       ) sf[tx‐1] = f[i‐1]; if(blockIdx.x< gridDim.x‐1 && threadIdx.x==blockDim.x‐1) sf[tx+1] = f[i+1]; f[i] sf[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 f[0] f[1] f[2]        f[2] f[3] f[4] f[5]

袖領域の処理

if(blockIdx.x> 0       && threadIdx.x==0       ) sf[tx‐1] = f[i‐1]; if(blockIdx.x< gridDim.x‐1 && threadIdx.x==blockDim.x‐1) sf[tx+1] = f[i+1]; f[i] sf[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 f[0] f[1] f[2] f[2] f[3] f[4] f[5] threadIdx.x= 0    1    2

袖領域の処理

if(blockIdx.x> 0       && threadIdx.x==0       ) sf[tx‐1] = f[i‐1]; if(blockIdx.x< gridDim.x‐1 && threadIdx.x==blockDim.x‐1) sf[tx+1] = f[i+1]; f[i] sf[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 f[0] f[1] f[2] f[3] f[2] f[3] f[4] f[5] threadIdx.x=2

境界条件の処理（差分式で計算）

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

3

f

3

f

f

f

_







Δx

f

f

Δx

f

f

f

dx

df

i

2

2

4

3

_{0 1 2 1 1} 0  













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

f

₋₁

を決定

f[‐1]

境界条件の処理（差分式で計算）

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

3

3

_



_



_



_{N N N} N

f

f

f

f

Δx

f

f

Δx

f

f

f

dx

df

_{N N N N N} N i

2

2

4

3

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











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

f

_N

を決定

f[‐1] f[6]

境界条件の処理（勾配が

0）

if(blockIdx.x==0       && threadIdx.x==0       )sf[tx‐1] = sf[tx+1]; if(blockIdx.x==gridDim.x‐1 && threadIdx.x==blockDim.x‐1)sf[tx+1] = sf[tx‐1]; sf[NT+2] tx=   0    1    2    3   4       0    1    2    3    4 blockIdx.x=0 blockIdx.x=1 f[0] f[1] f[2] f[3] f[2] f[3] f[4] f[5] 1 1

f

f

_



0

2

1 1 0







 

Δx

f

f

dx

df

i

境界での中心差分式が0となるようにf

₋₁

を決定

f[‐1]

境界条件の処理（勾配が

0）

if(blockIdx.x==0       && threadIdx.x==0       )sf[tx‐1] = sf[tx+1]; if(blockIdx.x==gridDim.x‐1 && threadIdx.x==blockDim.x‐1)sf[tx+1] = sf[tx‐1]; sf[NT+2] tx=   0    1    2    3   4       0    1    2    3    4 blockIdx.x=0 blockIdx.x=1 f[0] f[1] f[2] f[3] f[2] f[3] f[4] f[5] 2 



_N N

f

f

0

2

2 1







  

Δx

f

f

dx

df

_{N N} N i

境界での中心差分式が0となるようにf

_N

を決定

f[‐1] f[6]

中心差分の計算

__syncthreads(); dfdx[i] = ( sf[tx+1]‐sf[tx‐1])/(2.0*dx); sf[NT+2] tx=   0    1    2    3   4       0    1    2    3    4 blockIdx.x=0 blockIdx.x=1 f[0] f[1] f[2] f[3] f[2] f[3] f[4] f[5] dfdx[i] ＋ ＋ ＋ ＋ ＋ ＋

_{全スレッドが同じ式}

で中心差分を計算

Δx 2 1 f[‐1] f[6]

2次元への拡張



x方向1階偏微分



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



y方向１階偏微分



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

Δx

y

Δx

x

f

y

Δx

x

f

x

f

2

)

,

(

)

,

(













Δy

Δy

y

x

f

Δy

y

x

f

y

f

2

)

,

(

)

,

(













時間積分



時間微分項の離散化



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



右辺の

t+



tの項を移行

Δt

t

y

x

f

Δt

t

y

x

f

t

f

_

(

,

,



)



(

,

,

)



































y

f

v

x

f

u

t

f

t

f

Δt

t

y

x

f

Δt

t

y

x

f











)

(

,

,

)

,

,

(

移流方程式を代入





























y

t

y

x

f

v

x

t

y

x

f

u

Δt

t

y

x

f

Δt

t

y

x

f

(

,

,

)

(

,

,

)

(

,

,

)

(

,

,

)

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



簡略化した表現



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

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 ,

)

,

,

(

x

y

t



Δt



f

_in_j

f

離散化された移流方程式



連続系



離散系





t秒後の値

0

)

,

,

(

)

,

,

(

)

,

,

(

)

,

,

(

)

,

,

(



















y

t

y

x

f

t

y

x

v

x

t

y

x

f

t

y

x

u

t

t

y

x

f

0

2

2

1 , 1 , , , 1 , 1 , , 1 , 



_







_







_

Δy

f

f

v

Δx

f

f

u

Δt

f

f

_{n i}n_{j i}n_j j i n j i n j i n j i n j i n j i

















_









    

Δy

f

f

v

Δx

f

f

u

Δ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

2

1 , 1 , , , 1 , 1 , , 1 ,

移流方程式

x

y

i

i

+1

i

−1

j

j

+1

j

−1



x



y

x

y



x



y

t

t+



t

i

i

+1

i

−1

j

j

+1

j

−1

















_









    

Δy

f

f

v

Δx

f

f

u

Δ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

2

1 , 1 , , , 1 , 1 , , 1 ,

移流方程式の安定性



プログラムを正しく作成しても正常な計算結果が得

られない場合がある



安定条件



二つの条件を満たすことが必要（結果が正しいかは別）



_{CFL条件(Courant‐Friedrichs‐Lewy条件)}

1



Δx

Δt

u



1

Δy

Δt

v

















_









    

Δy

f

f

v

Δx

f

f

u

Δ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

2

1 , 1 , , , 1 , 1 , , 1 ,

計算手順

1.

計算条件の決定



計算領域の大きさ

L

_x

, L

_y 

計算領域の分割数（離散点の個数）

N

_x

, N

_y 

離散点同士の間隔（格子間隔）



x,



y



計算時間間隔



t

2.

初期値の決定



fの初期分布の決定



速度

u, vの分布の決定

3.

差分値の計算



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



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





t秒後のfを計算

CPUプログラム



計算条件の決定



計算領域の大きさ

L

_x

, L

_y



計算領域の分割数

（離散点の個数）

_N

_x

_{, N}

_y



離散点同士の間隔

（格子間隔）



_x,



_y



計算時間間隔



t

（CFL条件から逆算）

#include<stdio.h> #include<stdlib.h> #include<math.h> #define Lx 2.0 #define Ly 2.0 #define Nx 64 #define Ny Nx #define dx (Lx/(Nx‐1)) #define dy (Ly/(Ny‐1)) #define CFL 0.01 #define CONV (2.0*M_PI) #define Lt 1.0 #define dt (CFL*dx/CONV) #define Nt ((int)(Lt/(dt)))

CPUプログラム



初期条件



_{Rotating Cone問題}



境界条件

























2

1

0

2

1

1

2

cos

2

1

r

r

r

f



2

2

x

v

y

u





2

2











2 2

_

_

₀

_.

₅



x

y

r

x

y

r

f

1

0

0



f

正確ではないが計算を安定させるために便宜上

f

=0とする

CPUプログラム



初期条件



_{Rotating Cone問題}



境界条件

2

2

x

y

0



f





2 2

_

_

₀

_.

₅



x

y

r

























2

1

0

2

1

1

2

cos

2

1

r

r

r

f



正確ではないが計算を安定させるために便宜上

f

=0とする

x

v

y

u





2

2







#include<stdio.h> #include<stdlib.h> #include<math.h> #define Lx 2.0 #define Ly Lx #define Nx 128 #define Ny Nx #define dx (Lx/(Nx‐1)) #define dy dx #define CFL 0.01 #define CONV (2.0*M_PI) #define Lt 1.0 #define dt (CFL*dx/CONV) #define Nt ((int)(Lt/(dt))) void init(double *f,double *d_f_dx, double *d_f_dy, double *u, double *v); void differentiate(double *f,double *d_f_dx, double *d_f_dy); void integrate(double *f,double *d_f_dx, double *d_f_dy, double *u, double *v,double *fnew); void update(double *f,double *fnew); int main(void){ double *f, *d_f_dx, *d_f_dy,*fnew; double *u,*v; int n; f = (double *)malloc(Nx*Ny*sizeof(double)); d_f_dx = (double *)malloc(Nx*Ny*sizeof(double)); d_f_dy = (double *)malloc(Nx*Ny*sizeof(double)); u = (double *)malloc(Nx*Ny*sizeof(double)); v = (double *)malloc(Nx*Ny*sizeof(double)); fnew = (double *)malloc(Nx*Ny*sizeof(double)); init(f,d_f_dx,d_f_dy,u,v); for(n=1;n<=Nt;n++){ printf("%d/%d¥n",n,Nt); differentiate(f,d_f_dx,d_f_dy); integrate(f,d_f_dx,d_f_dy, u, v,fnew); update(f,fnew); } return 0; }

CPUプログラム

rotating_cone.c

void init(double *f,double *d_f_dx,double *d_f_dy,  double *u, double *v){ int i,j; double x,y,r; for(i=0;i<Nx;i++){ for(j=0;j<Ny;j++){ x = (double)i*dx ‐ 1.0; y = (double)j*dy ‐ 1.0; r = sqrt(x*x + (y‐0.5)*(y‐0.5)); f[i+Nx*j] = 0.0; if(r<0.5) f[i+Nx*j] = 0.5*(cos(2.0*M_PI*r)+1.0); d_f_dx[i+Nx*j] = 0.0; d_f_dy[i+Nx*j] = 0.0; u[i+Nx*j] = ‐(2.0*M_PI*y); v[i+Nx*j] =  (2.0*M_PI*x); } } } void differentiate(double *f,double *d_f_dx, double *d_f_dy){ int i,j; for(i=1;i<Nx‐1;i++){ for(j=1;j<Ny‐1;j++){ d_f_dx[i+Nx*j] = (f[(i+1)+Nx*j]‐f[(i‐1)+Nx*j])/(2.0*dx); } } for(i=1;i<Nx‐1;i++){ for(j=1;j<Ny‐1;j++){ d_f_dy[i+Nx*j] = (f[i+Nx*(j+1)]‐f[i+Nx*(j‐1)])/(2.0*dy); } } }

CPUプログラム

rotating_cone.c

void integrate(double *f,double *d_f_dx, double *d_f_dy, double *u, double *v,double *fnew){ int i,j; for(i=0;i<Nx;i++){ for(j=0;j<Ny;j++){ fnew[i+Nx*j] = f[i+Nx*j] ‐ dt*(  u[i+Nx*j]*d_f_dx[i+Nx*j] + v[i+Nx*j]*d_f_dy[i+Nx*j]); } } } void update(double *f,double *fnew){ int i,j; for(i=0;i<Nx;i++){ for(j=0;j<Ny;j++){ f[i+Nx*j] = fnew[i+Nx*j]; } } }

CPUプログラム

rotating_cone.c

CPUプログラム



差分計算



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

void differentiate(double *f, double *d_f_dx, double *d_f_dy){ int i,j; for(i=1;i<Nx‐1;i++){ for(j=1;j<Ny‐1;j++){ d_f_dx[i+Nx*j] = (f[(i+1)+Nx*j]‐f[(i‐1)+Nx*j])/(2.0*dx); } } for(i=1;i<Nx‐1;i++){ for(j=1;j<Ny‐1;j++){ d_f_dy[i+Nx*j] = (f[i+Nx*(j+1)]‐f[i+Nx*(j‐1)])/(2.0*dy); } } }

差分計算のメモリ参照



ある1点i,jのx方向差分を計算するためにx(i)方

向に2点のfを参照

f[]

d_f_dx[]

差分計算のメモリ参照



ある1点i,jのx方向差分を計算するためにx(i)方

向に2点のfを参照

f[]

d_f_dx[]

差分計算のメモリ参照



ある1点i,jのx方向差分を計算するためにx(i)方

向に2点のfを参照

f[]

d_f_dx[]

差分計算のメモリ参照



ある1点i,jのy方向差分を計算するためにy(j)方

向に2点のfを参照

j−

1

jj

+1

f[]

d_f_dy[]

差分計算のメモリ参照



ある1点i,jのy方向差分を計算するためにy(j)方

向に2点のfを参照

f[]

d_f_dy[]

j−

1

jj

+1

差分計算のメモリ参照



ある1点i,jのy方向差分を計算するためにy(j)方

向に2点のfを参照

f[]

d_f_dy[]

j−

1

jj

+1

差分計算のメモリ参照



ある1点i,jのy方向差分を計算するためにy(j)方

向に2点のfを参照

f[]

d_f_dy[]

j−

1

jj

+1

CPUプログラム



境界条件



fのx, y方向差分



いずれの時刻においても

0



境界で

f=0を満たすためには差分

が

0でなければならない

f,d_f_dx,d_f_dy





























y

f

v

x

f

u

Δt

f

f

n 1 n

CPUプログラム



fの積分

















_









    

Δy

f

f

v

Δx

f

f

u

Δ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

2

1 , 1 , , , 1 , 1 , , 1 ,

void integrate(double *f,double *d_f_dx,double *d_f_dy,

double *u, double *v,double *f_new){ int i,j,ij; for(i=0;i<Nx;i++){ for(j=0;j<Ny;j++){ ij = i+Nx*j; f_new[ij] = f[ij] ‐ dt*(u[ij]*d_f_dx[ij] + v[ij]*d_f_dy[ij]); } } }

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(i=0;i<Nx;i++){ for(j=0;j<Ny;j++){ ij =  i+Nx*j; f[ij] = f_new[ij]; } } }

#include<stdio.h> #include<stdlib.h> #include<math.h> #define Lx 2.0 #define Ly Lx #define Nx 64 #define Ny Nx #define dx (Lx/(Nx‐1)) #define dy dx #define Nbytes (Nx*Ny*sizeof(double)) #define CFL 0.01 #define CONV (2.0*M_PI) #define Lt 1.0 #define dt (CFL*dx/CONV) #define Nt ((int)(Lt/(dt))) #include "conv1.cu" int main(void){ double *dev_f,*dev_fnew; double *dev_d_f_dx, *dev_d_f_dy; double *dev_u,*dev_v; dim3 Thread, Block; int n; cudaMalloc( (void **)&dev_f , Nbytes ); cudaMalloc( (void **)&dev_d_f_dx, Nbytes ); cudaMalloc( (void **)&dev_d_f_dy, Nbytes ); cudaMalloc( (void **)&dev_u , Nbytes ); cudaMalloc( (void **)&dev_v , Nbytes ); cudaMalloc( (void **)&dev_fnew , Nbytes ); Thread = dim3(THREADX,THREADY,1); Block  = dim3(BLOCKX ,BLOCKY, 1); init<<<Block, Thread>>>(dev_f, dev_d_f_dx,dev_d_f_dy,dev_u,dev_v);

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

rotating_cone.cu

for(n=1;n<=Nt;n++){ printf("%d/%d¥n",n,Nt); differentiate<<<Block,Thread>>> (dev_f,dev_d_f_dx,dev_d_f_dy); integrate<<<Block,Thread>>> (dev_f,dev_d_f_dx,dev_d_f_dy, dev_u, dev_v, dev_fnew); update<<<Block,Thread>>> (dev_f,dev_fnew); } cudaFree(dev_f); cudaFree(dev_d_f_dx); cudaFree(dev_d_f_dy); cudaFree(dev_u); cudaFree(dev_v); cudaFree(dev_fnew); return 0; }

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

rotating_cone.cu

#define THREADX 16 #define THREADY 16 #define BLOCKX  (Nx/THREADX) #define BLOCKY  (Ny/THREADY) __global__ void init(double *f,double *d_f_dx, double *d_f_dy, double *u, double *v){ int i,j,ij; double x,y,r; i = blockIdx.x*blockDim.x + threadIdx.x; j = blockIdx.y*blockDim.y + threadIdx.y; ij = i+Nx*j; x = (double)i*dx ‐ 1.0; y = (double)j*dy ‐ 1.0; r = sqrt(x*x + (y‐0.5)*(y‐0.5)); f[ij] = 0.0; if(r<0.5) f[ij] = 0.5*(cos(2.0*M_PI*r)+1.0); d_f_dx[ij] = 0.0; d_f_dy[ij] = 0.0; u[ij] = ‐(2.0*M_PI*y); v[ij] =  (2.0*M_PI*x); } __global__ void differentiate(double *f, double *d_f_dx,double *d_f_dy){ 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); d_f_dx[ij] = (f[ip1]‐f[im1])/(2.0*dx); d_f_dy[ij] = (f[jp1]‐f[jm1])/(2.0*dy); } }

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

conv1.cu

__global__ void integrate(double *f,double *d_f_dx, double *d_f_dy, double *u, double *v,double *fnew){ int i,j,ij; i = blockIdx.x*blockDim.x + threadIdx.x; j = blockIdx.y*blockDim.y + threadIdx.y; ij = i+Nx*j; fnew[ij] = f[ij] ‐ dt* (u[ij]*d_f_dx[ij] + v[ij]*d_f_dy[ij]); } __global__ void update(double *f, double *fnew){ int i,j,ij; i = blockIdx.x*blockDim.x + threadIdx.x; j = blockIdx.y*blockDim.y + threadIdx.y; ij =  i+Nx*j; f[ij] = fnew[ij]; }

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

conv1.cu

gnuplotによる結果の表示



_{2次元，3次元データをプロットするアプリケーション}



コマンドラインで命令を実行してグラフを描画



関数の描画，ファイルから読み込んだデータの表示が可能



_{tesla??では正しく動作しないため，grouseで実行}

すること

gnuplotによる結果の表示



表示可能なファイルフォーマット



データはスペース区切り



_{3次元表示では列（または行）の区切りとして空白行が必要}

‐1.000000 ‐1.000000 0.000000 ‐1.000000 ‐0.968254 0.000000 ... ‐1.000000 0.968254 0.000000 ‐1.000000 1.000000 0.000000 ‐0.968254 ‐1.000000 0.000000 ‐0.968254 ‐0.968254 0.000000 ... ‐0.968254 0.968254 0.000000 ‐0.968254 1.000000 0.000000 ‐0.936508 ‐1.000000 0.000000 ‐0.936508 ‐0.968254 0.000000

gnuplotによる結果の表示



スクリプトの実行



_{gnuplotのコマンドを記述したファイルを実行}



_{load 'スクリプトファイル名'}



_anim.gpl



結果を単色の等値面で3次元表示



_{anim_color.gpl}



結果をカラーの等値面で3次元表示



_{anim_2d.gpl}



結果をカラーの等高線で2次元表示

gnuplotによる結果の表示



スクリプトanim.gplの内容

set xrange [‐1:1] x軸の表示範囲を‐1~1に固定 set yrange [‐1:1] y軸の表示範囲を‐1~1に固定 set zrange [‐0.2:1.2] z軸の表示範囲を‐0.2~1.2に固定 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 [‐1:1] x軸の表示範囲を‐1~1に固定 set yrange [‐1:1] y軸の表示範囲を‐1~1に固定 set zrange [‐0.2:1.2] z軸の表示範囲を‐0.2~1.2に固定 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 カラー表示する 以下でファイルを読み込み，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 [‐1:1] x軸の表示範囲を‐1~1に固定 set yrange [‐1:1] y軸の表示範囲を‐1~1に固定 set zrange [‐0.2:1.2] z軸の表示範囲を‐0.2~1.2に固定 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による結果の表示

ダブルバッファリングによる値の更新



ダブルバッファリング



ゲーム等で画面のチラつきを抑える技法とは異なる



値の更新を，二つの配列を用いて簡単に実行する方法



移流方程式（拡散方程式）の時間積分



時刻nのデータからn+1のデータを計算



時刻n+1のデータを時刻nのデータに上書き



メモリコピーが発生



ポインタを利用



アドレスを格納する配列を宣言し，配列添字を切替

int a,b;

int *ptr[2];

int curr, next;

curr = 0; next = 1;

ptr[curr]=&a;

prt[next]=&b;

アドレスを格納する配列を用いる方法



時刻n，時刻n+1の変数を用意せず，ポインタ変数を

取り扱う配列を利用



時刻n用とするか時刻n+1用途するかを配列添字で区別



_{int型変数curr(currentの略), nextを利用}



配列添字にcurrを用いると時刻n用



配列添字にnextを用いると時刻n+1用

a

1

2

b

ptr[curr]

aのア ドレス bのア ドレス

ptr[next]

結合 ポインタ変数が矢印先端 の変数のアドレスを保持 しており，間接参照演算 子*を利用してその変数 の値にアクセスできる 値 アド レス 変数 ポインタ 変数

‐‐‐

ptr[curr]はaのアドレスを参照

ptr[next]はbのアドレスを参照

‐‐‐

curr=next;   

//currは1

next=1‐curr; 

//nextは1‐1=0

ptr[curr]はbのアドレスを参照

ptr[next]はaのアドレスを参照

‐‐‐

curr=next;   

//currは0

next=1‐curr; 

//nextは1‐0=1

アドレスを格納する配列を用いる方法



以下の2式で配列添字curr, nextの値を切替



_{curr = next;}



_{next = 1‐curr;}

a

1

2

b

ptr[next

(=0)

]

aのア ドレス bのア ドレス

ptr[curr

(=1)

]

アドレスを格納する配列を用いる方法



無駄がない



時刻n+1の値をnとすることを表現（curr=next）

double *f[2];

int curr=0,next=1;

f[curr

(=0)

] = (double *)malloc(...);

f[next

(=1)

] = (double *)malloc(...);

differentiate(f[curr

(=0)

]);

integrate(f[curr

(=0)

],f[next

(=1)

]);

//添字を変更

curr = next;

next = 1‐curr;

//時刻nと時刻n+1を交換して関数を実行

differentiate(f[curr

(=1)

]);

integrate(f[curr

(=1)

],f[next

(=0)

]);

アドレスを格納する配列を用いる方法



関数integrateの実行



実行後，時刻nの値は不要

f[0]

時刻n

f[1]

時刻n+1

integrate(f[curr

(=0)

],...,f[next

(=1)

]);

アドレスを格納する配列を用いる方法



関数integrateの実行



実行後，時刻nの値は不要

f[0]

時刻n+1

f[1]

時刻n

integrate(f[curr

(=1)

],...,f[next

(=0)

]);

#include<stdio.h> #include<stdlib.h> #include<math.h> #define Lx 2.0 #define Ly Lx #define Nx 64 #define Ny Nx #define dx (Lx/(Nx‐1)) #define dy dx #define Nbytes (Nx*Ny*sizeof(double)) #define CFL 0.01 #define CONV (2.0*M_PI) #define Lt 1.0 #define dt (CFL*dx/CONV) #define Nt ((int)(Lt/(dt))) #include "conv1.cu" int main(void){ double *dev_d_f_dx, *dev_d_f_dy; double *dev_u,*dev_v; dim3 Thread, Block; int n; double *dev_f[2]; int curr=0,next=1; cudaMalloc( (void **)&dev_d_f_dx, Nbytes ); cudaMalloc( (void **)&dev_d_f_dy, Nbytes ); cudaMalloc( (void **)&dev_u , Nbytes ); cudaMalloc( (void **)&dev_v , Nbytes ); cudaMalloc( (void **)&dev_f[curr], Nbytes ); cudaMalloc( (void **)&dev_f[next], Nbytes ); Thread = dim3(THREADX,THREADY,1); Block  = dim3(BLOCKX ,BLOCKY, 1); init<<<Block, Thread>>>(dev_f[curr], dev_d_f_dx,dev_d_f_dy,dev_u,dev_v);

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

rotating_cone2.cu

for(n=1;n<=Nt;n++){ printf("%d/%d¥n",n,Nt); differentiate<<<Block,Thread>>> (dev_f[curr],dev_d_f_dx,dev_d_f_dy); integrate<<<Block,Thread>>> (dev_f[curr],dev_d_f_dx,dev_d_f_dy,  dev_u, dev_v, dev_f[next]); curr = next; next = 1‐curr; } cudaFree(dev_f[curr]); cudaFree(dev_f[next]); cudaFree(dev_d_f_dx); cudaFree(dev_d_f_dy); cudaFree(dev_u); cudaFree(dev_v); return 0; }

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

rotating_cone2.cu

レポート課題

3（提出期限は1学期末）



_{GPUで2次元移流拡散方程式を計算するプログラム}

を作成せよ



移流拡散方程式



必ずcudaSetDevice命令を利用してGPUを選択せよ



初期条件，境界条件，計算終了時刻はRotating Coneと

同じとする



カーネルは独自に最適化せよ



配列のサイズも自身で決定せよ













































2 2 2 2

y

f

x

f

y

f

v

x

f

u

t

f

_