質問をすることでしか得られない、回答やアドバイスがある。

15分調べてもわからないことは、質問しよう!

新規登録して質問してみよう
ただいま回答率
85.50%
C

C言語は、1972年にAT&Tベル研究所の、デニス・リッチーが主体となって作成したプログラミング言語です。 B言語の後継言語として開発されたことからC言語と命名。そのため、表記法などはB言語やALGOLに近いとされています。 Cの拡張版であるC++言語とともに、現在世界中でもっとも普及されているプログラミング言語です。

CUDA

CUDAは並列計算プラットフォームであり、Nvidia GPU(Graphics Processing Units)向けのプログラミングモデルです。CUDAは様々なプログラミング言語、ライブラリ、APIを通してNvidiaにインターフェイスを提供します。

Q&A

解決済

1回答

597閲覧

cudaを用いたFDTD計算でのMurの1次吸収境界条件について

run919

総合スコア10

C

C言語は、1972年にAT&Tベル研究所の、デニス・リッチーが主体となって作成したプログラミング言語です。 B言語の後継言語として開発されたことからC言語と命名。そのため、表記法などはB言語やALGOLに近いとされています。 Cの拡張版であるC++言語とともに、現在世界中でもっとも普及されているプログラミング言語です。

CUDA

CUDAは並列計算プラットフォームであり、Nvidia GPU(Graphics Processing Units)向けのプログラミングモデルです。CUDAは様々なプログラミング言語、ライブラリ、APIを通してNvidiaにインターフェイスを提供します。

0グッド

1クリップ

投稿2018/10/13 10:09

編集2018/10/13 10:43

前提・実現したいこと

NVIDIA社のcudaをvisualsudio2015用いて3次元FDTD法を練習しています。今回、サイズ41×41×41の立方体中のx=20,y=20の位置における0<z<41のEzに電圧を常時印加した場合に、z=20におけるxy平面のEzの値を見ているのですが、普通にCでこのプログラムを作った場合はきれいに最後まで円形の画像が出力され続けるのですが、cudaでは出力画像でx=0とy=0付近においてmurの1次吸収境界条件の効きが悪く、それが原因で発散しているような状況となっています。(吸収境界条件を入れずに完全導体で囲んだ場合はきれいに反射された画像が続く)
正直自分ではもうどこが間違っているのかわからないので、質問させていただきました。おそらく吸収境界条件が間違っていると考えているのですが。
動画ではなく、画像ですがいくつか掲載させていただきます。(プログラムのdt(電界計算+磁界計算時間)は2e-12s)

イメージ説明
___________________t=0.8e-10s
イメージ説明
___________________t=3.0e-10s
イメージ説明
___________________t=3.2e-10s
イメージ説明
___________________t=3.4e-10s
イメージ説明
___________________t=3.6e-10s
![イメージ説明
___________________t=4.0e-10s
イメージ説明
___________________t=2.0e-9s

該当のソースコード

/*--------------------------------------------メイン計算部分------------------------------------------------------------*/ while (t[0]<2e-9) { culcEx << <grid, block >> > (Ex_d, oldEx_d, oldHy_d, oldHz_d, A_d, d_d); culcEy << <grid, block >> > (Ey_d, oldEy_d, oldHx_d, oldHz_d, A_d, d_d); culcEz << <grid, block >> > (Ez_d, oldEz_d, oldHx_d, oldHy_d, A_d, d_d); cudaDeviceSynchronize(); vin[0] = 1.0*sin(2.0e10*M_PI*t[0]); cudaMemcpy(vin_d, vin, sizeof(double), cudaMemcpyHostToDevice); Init << <grid, block >> > ( Ez_d, t_d, vin_d); cudaDeviceSynchronize(); Init2 << <grid, block >> > (Ex_d, Ey_d, Ez_d); cudaDeviceSynchronize(); MurEx << <grid, block >> > (Ex_d, oldEx_d, A_d); MurEy << <grid, block >> > (Ey_d, oldEy_d, A_d); MurEz << <grid, block >> > (Ez_d, oldEz_d, A_d); cudaDeviceSynchronize(); culct << <2,2 >> > (t_d); t[0] += (dt / 2); cudaDeviceSynchronize(); culcHx << <grid, block >> > (Ey_d, Ez_d, Hx_d, oldHx_d, A_d, d_d); culcHy << <grid, block >> > (Ex_d, Ez_d, Hy_d, oldHy_d, A_d, d_d); culcHz << <grid, block >> > (Ex_d, Ey_d, Hz_d, oldHz_d, A_d, d_d); cudaDeviceSynchronize(); change << <grid, block >> > (Ex_d, Ey_d, Ez_d, oldEx_d, oldEy_d, oldEz_d, Hx_d, Hy_d, Hz_d, oldHx_d, oldHy_d, oldHz_d); cudaDeviceSynchronize(); culct << <2, 2 >> > (t_d); t[0] += (dt / 2); cudaDeviceSynchronize(); output << <grid, block >> > (out_d, Ez_d); cudaDeviceSynchronize(); } /*---------------------------------------------------------------------------------------------------------------------*/ /*------------------吸収境界条件部分----------------------------*/ __global__ void MurEx(double *Ex, double *oldEx, double *A){ int i = blockIdx.x*blockDim.x + threadIdx.x; int j = blockIdx.y*blockDim.y + threadIdx.y; int k = blockIdx.z*blockDim.z + threadIdx.z; if (i < nx && j < ny && k < nz) { Ex[(k * ny * nx + j * nx + 0)] = oldEx[(k * ny * nx + j * nx + 1)] + A[3] * (Ex[(k * ny * nx + j * nx + 1)] - oldEx[(k * ny * nx + j * nx + 0)]); Ex[(k * ny * nx + j * nx + nx - 1)] = oldEx[(k * ny * nx + j * nx + nx - 2)] + A[3] * (Ex[(k * ny * nx + j * nx + nx - 2)] - oldEx[(k * ny * nx + j * nx + nx - 1)]); } __syncthreads(); } __global__ void MurEy( double *Ey, double *oldEy, double *A){ int i = blockIdx.x*blockDim.x + threadIdx.x; int j = blockIdx.y*blockDim.y + threadIdx.y; int k = blockIdx.z*blockDim.z + threadIdx.z; if (i < nx && j < ny && k < nz) { Ey[(k * ny * nx + 0 * nx + i)] = oldEy[(k * ny * nx + 1 * nx + i)] + A[3] * (Ey[(k * ny * nx + 1 * nx + i)] - oldEy[(k * ny * nx + 0 * nx + i)]); Ey[(k * ny * nx + (ny - 1) * nx + i)] = oldEy[(k * ny * nx + (ny - 2) * nx + i)] + A[3] * (Ey[(k * ny * nx + (ny - 2) * nx + i)] - oldEy[(k * ny * nx + (ny - 1) * nx + i)]); } __syncthreads(); } __global__ void MurEz( double *Ez, double *oldEz, double *A){ int i = blockIdx.x*blockDim.x + threadIdx.x; int j = blockIdx.y*blockDim.y + threadIdx.y; int k = blockIdx.z*blockDim.z + threadIdx.z; if (i < nx && j < ny && k < nz) { Ez[(k * ny * nx + j * nx + 0)] = oldEz[(k * ny * nx + j * nx + 1)] + A[3] * (Ez[(k * ny * nx + j * nx + 1)] - oldEz[(k * ny * nx + j * nx + 0)]); Ez[(k * ny * nx + j * nx + nx - 1)] = oldEz[(k * ny * nx + j * nx + nx - 2)] + A[3] * (Ez[(k * ny * nx + j * nx + nx - 2)] - oldEz[(k * ny * nx + j * nx + nx - 1)]); Ez[(k * ny * nx + 0 * nx + i)] = oldEz[(k * ny * nx + 1 * nx + i)] + A[3] * (Ez[(k * ny * nx + 1 * nx + i)] - oldEz[(k * ny * nx + 0 * nx + i)]); Ez[(k * ny * nx + (ny - 1) * nx + i)] = oldEz[(k * ny * nx + (ny - 2) * nx + i)] + A[3] * (Ez[(k * ny * nx + (ny - 2) * nx + i)] - oldEz[(k * ny * nx + (ny - 1) * nx + i)]); } __syncthreads(); } /*------------------------------------------------------------------------------*/ //その他補足(プログラムより抜粋) ♯define nx 41 ♯define ny 41 ♯define nz 41 dim3 grid; grid.x = 3; grid.y = 3; grid.z = 24; dim3 block; block.x = 16; block.y = 16; block.z = 2; /*--------------------------電磁界計算用係数設定------------------------*/ A[0] = ((1 - (sig*dt) / (2 * ep)) / (1 + (sig*dt) / (2 * ep))); A[1] = (dt / ep) / (1 + (sig*dt / (2 * ep))); A[2] = dt / mu; A[3] = (c*dt - dx) / (c*dt + dx); /*----------------------------------------------- 0:電界計算,1:電界計算,2:磁界計算,3:Mur吸収境界条件 -------------------------------------------------*/ d[0] = 2.0e-12; d[1] = 1.5e-3; d[2] = 1.5e-3; d[3] = 1.5e-3; /*----------------------------------------------- 0:dt,1:dx,2:dy,3:dz -------------------------------------------------*/ /*----------------------------電界計算----------------------*/ __global__ void culcEx(double *Ex, double *oldEx, double *oldHy, double *oldHz, double *A, double *d) { int i = blockIdx.x*blockDim.x + threadIdx.x; int j = blockIdx.y*blockDim.y + threadIdx.y+1; int k = blockIdx.z*blockDim.z + threadIdx.z+1; if ((0 < j) && (0 < k) && (i < nx) && (j < ny) && (k < nz)) { Ex[(k * ny * nx + j * nx + i)] = A[0] * oldEx[(k * ny * nx + j * nx + i)] + A[1] * (((oldHz[(k * ny * nx + j * nx + i)] - oldHz[(k * ny * nx + (j - 1) * nx + i)]) / d[2]) - ((oldHy[(k * ny * nx + j * nx + i)] - oldHy[((k - 1) * ny * nx + j * nx + i)]) / d[3])); //4 } __syncthreads(); } __global__ void culcEy(double *Ey, double *oldEy, double *oldHx, double *oldHz, double *A, double *d) { int i = blockIdx.x*blockDim.x + threadIdx.x+1; int j = blockIdx.y*blockDim.y + threadIdx.y; int k = blockIdx.z*blockDim.z + threadIdx.z+1; if ((0<i) && (0<k) && (i < nx) && (j < ny) && (k < nz)) { Ey[(k * ny * nx + j * nx + i)] = A[0] * oldEy[(k * ny * nx + j * nx + i)] + A[1] * (((oldHx[(k * ny * nx + j * nx + i)] - oldHx[((k - 1) * ny * nx + j * nx + i)]) / d[3]) - ((oldHz[(k * ny * nx + j * nx + i)] - oldHz[(k * ny * nx + j * nx + (i - 1))]) / d[3])); //4 } __syncthreads(); } __global__ void culcEz(double *Ez, double *oldEz, double *oldHx, double *oldHy, double *A, double *d) { int i = blockIdx.x*blockDim.x + threadIdx.x+1; int j = blockIdx.y*blockDim.y + threadIdx.y+1; int k = blockIdx.z*blockDim.z + threadIdx.z; if ((0<i) && (0<j) && (i < nx) && (j < ny) && (k < nz)) { Ez[(k * ny * nx + j * nx + i)] = A[0] * oldEz[(k * ny * nx + j * nx + i)] + A[1] * (((oldHy[(k * ny * nx + j * nx + i)] - oldHy[(k * ny * nx + j * nx + (i - 1))]) / d[1]) - ((oldHx[(k * ny * nx + j * nx + i)] - oldHx[(k * ny * nx + (j - 1) * nx + i)]) / d[2])); //4 } __syncthreads(); } /*----------------------------磁界計算----------------------*/ __global__ void culcHx(double *Ey, double *Ez, double *Hx, double *oldHx, double *A, double *d) { int i = blockIdx.x*blockDim.x + threadIdx.x; int j = blockIdx.y*blockDim.y + threadIdx.y; int k = blockIdx.z*blockDim.z + threadIdx.z; if ((0<i)&&(i<nx) && (j<ny - 1) && (k<nz - 1)) { Hx[k * ny * nx + j * nx + i] = oldHx[k * ny * nx + j * nx + i] - (A[2] * (((Ez[k * ny * nx + (j + 1) * nx + i] - Ez[k * ny * nx + j * nx + i]) / d[2]) - ((Ey[(k + 1) * ny * nx + j * nx + i] - Ey[k * ny * nx + j * nx + i]) / d[3]))); //5 } __syncthreads(); } __global__ void culcHy(double *Ex, double *Ez, double *Hy, double *oldHy, double *A, double *d) { int i = blockIdx.x*blockDim.x + threadIdx.x; int j = blockIdx.y*blockDim.y + threadIdx.y; int k = blockIdx.z*blockDim.z + threadIdx.z; if ((0<j) && (i<nx - 1) && (j<ny) && (k<nz - 1)) { Hy[k * ny * nx + j * nx + i] = oldHy[k * ny * nx + j * nx + i]- (A[2] * (((Ex[(k + 1) * ny * nx + j * nx + i] - Ex[k * ny * nx + j * nx + i]) / d[3]) - ((Ez[k * ny * nx + j * nx + (i + 1)] - Ez[k * ny * nx + j * nx + i]) / d[1]))); //6 } __syncthreads(); } __global__ void culcHz(double *Ex, double *Ey, double *Hz, double *oldHz, double *A, double *d) { int i = blockIdx.x*blockDim.x + threadIdx.x; int j = blockIdx.y*blockDim.y + threadIdx.y; int k = blockIdx.z*blockDim.z + threadIdx.z; if ((0<j) && (i<nx - 1) && j<(ny - 1) && k<nz) { Hz[k * ny * nx + j * nx + i] = oldHz[k * ny * nx + j * nx + i]- (A[2] * (((Ey[k * ny * nx + j * nx + (i + 1)] - Ey[k * ny * nx + j * nx + i]) / d[1]) - ((Ex[k * ny * nx + (j + 1) * nx + i] - Ex[k * ny * nx + j * nx + i]) / d[2]))); //6 } __syncthreads(); }

試したこと

使用している関数のif条件を様々なパターンで試しました。
0から始まるところを1にそろえてみたり、部分的に1にしてみたり、
i = blockIdx.x*blockDim.x + threadIdx.x+1;のようにしてみたり
メイン関数中のMurE関数がくるところを変更してCPUで計算させてもみましたがダメでした。

補足情報(FW/ツールのバージョンなど)

cuda runtime 9.2
windows10
visualstudio2015 update3

気になる質問をクリップする

クリップした質問は、後からいつでもMYページで確認できます。

またクリップした質問に回答があった際、通知やメールを受け取ることができます。

バッドをするには、ログインかつ

こちらの条件を満たす必要があります。

ruei

2018/10/13 10:16

①電場分布の時刻変化の動画を載せられませんか?②コードは「```」で囲ってください
ruei

2018/10/13 10:18

発生している問題・エラーメッセージも、サンプルのままになっています。
guest

回答1

0

自己解決

Murの吸収境界条件において、Ex,Eyの計算が逆になっていた。

投稿2018/10/17 04:58

run919

総合スコア10

バッドをするには、ログインかつ

こちらの条件を満たす必要があります。

あなたの回答

tips

太字

斜体

打ち消し線

見出し

引用テキストの挿入

コードの挿入

リンクの挿入

リストの挿入

番号リストの挿入

表の挿入

水平線の挿入

プレビュー

15分調べてもわからないことは
teratailで質問しよう!

ただいまの回答率
85.50%

質問をまとめることで
思考を整理して素早く解決

テンプレート機能で
簡単に質問をまとめる

質問する

関連した質問