FDTD法を用いた電磁界解析にOpenMPを導入

質問

FDTD法を用いた電磁界解析にOpenMPを導入するにはどうしたらよいのでしょうか？
ちなみに３次元のプログラムを作成しており、電磁界解析のエンジン部分は下記の様になっています。

for( kt=0; kt<=it ; kt++)
{ /*　※※※forの括弧開始※※※　*/

/*********** 1ステップ前の電磁界[jt-1]の値の定式 ***********/
for( jx=0; jx<=ix; jx++ )
for( jy=0; jy<=iy; jy++ )
for( jz=0; jz<=iz; jz++ )
{
exg[jx][jy][jz] = ex[jx][jy][jz];
eyg[jx][jy][jz] = ey[jx][jy][jz];
ezg[jx][jy][jz] = ez[jx][jy][jz];

}

/*********** 磁界(hx)の計算 ***********/
for( jx=0; jx<=ix; jx++ )
for( jy=0; jy<=iy-1; jy++ )
for( jz=0; jz<=iz-1; jz++ )
{
double dey = (ey[jx][jy][jz+1]-ey[jx][jy][jz])/dz;
double dez = (ez[jx][jy+1][jz]-ez[jx][jy][jz])/dy;
hx[jx][jy][jz] = hx[jx][jy][jz]+dt/amuy[jx][jy][jz]*(dey-dez);

}

/*********** 磁界(hy)の計算 ***********/
for( jx=0; jx<=ix-1; jx++ )
for( jy=0; jy<=iy; jy++ )
for( jz=0; jz<=iz-1; jz++ )
{
double dez = (ez[jx+1][jy][jz]-ez[jx][jy][jz])/dx;
double dex = (ex[jx][jy][jz+1]-ex[jx][jy][jz])/dz;
hy[jx][jy][jz] = hy[jx][jy][jz]+dt/amuz[jx][jy][jz]*(dez-dex);

}

/*********** 磁界(hz)の計算 ***********/
for( jx=0; jx<=ix-1; jx++ )
for( jy=0; jy<=iy-1; jy++ )
for( jz=0; jz<=iz; jz++ )
{
double dex = (ex[jx][jy+1][jz]-ex[jx][jy][jz])/dy;
double dey = (ey[jx+1][jy][jz]-ey[jx][jy][jz])/dx;
hz[jx][jy][jz] = hz[jx][jy][jz]+dt/amuz[jx][jy][jz]*(dex-dey);
}

/*********** 励振 ***********/
for( jx=1; jx<=ix-1; jx++ )
for( jz=1; jz<=iz-2; jz++ )
{
ez[jx][80][jz] = sin( 2.0*pi*(kt%100)/100 )*sin(pi*float(jx)/float(ix));
}
/*********** 電界(Ez)の計算 ***********/
for( jx=1; jx<=ix-1; jx++ )
for( jy=1; jy<=iy-1; jy++ )
for( jz=0; jz<=iz-1; jz++ )
{
double term1 = (1.0-sigmax[jx][jy][jz]*dt/(2.0*epsirx[jx][jy][jz]))
/(1.0+sigmax[jx][jy][jz]*dt/(2.0*epsirx[jx][jy][jz]));
double term2 = 1.0/(1.0+sigmax[jx][jy][jz]*dt/(2.0*epsirx[jx][jy][jz]));
double dhy = (hy[jx][jy][jz]-hy[jx-1][jy][jz])/dx;
double dhx = (hx[jx][jy][jz]-hx[jx][jy-1][jz])/dy;
ez[jx][jy][jz] = term1*ez[jx][jy][jz]+dt/epsirx[jx][jy][jz]*term2*(dhy-dhx);
}
・
・
・
この後電界の計算が同様に続きます。

どなたかこのプログラムにOpenMPを導入するにはどうすればよいか教えていただきたいです。

回答No.4

ki073

>三重の場合はどのようにすればよいのでしょうか？
変数をカンマで区切って並べるだけです
#pragma omp parallel for private(lwjy, lwjz)
のような感じです

回答No.3

ki073

No.1のお礼欄のソースコードですが、
例えば***領域II***の場合、ループの中に
( cos(lwux)*cos((lwuz/lwtz)*(lwkz0-lwjz)*lwdz-lwfhi)/exp((lwwx/lwtx)*((lwkx0-lwjx)*lwdx-lwtx)) )*lwaz;
とありますが、
cos(lwux)は二重ループの外に
/exp((lwwx/lwtx)*(lwkx0-lwjx)*lwdx-lwtx)*lwaz
は１つ外のループに移せますので、それでかなり計算量が減らせます。三角関数やexpはかなり時間がかかるので、結構効果が期待できます。

コンパイラにもよるのですが、gcc(4.4.7)だとループの外に出す最適化はコンパイラでは行われていないようで、書き換え効果はかなり出そうです。
私が普段使っているpgiのコンパイラだと、最適化の過程である程度は自動的にループの外に出してくれるのでわざわざ出すことも無さそうですが。
ソースコードが極端に見ずらくならない範囲で、ソースコードレベルで最適化されたら良いと思います。

openMP化をする前に、このような最適化を行うのが良いと思います。

openMP(4coreとして)と最適化の効果で、10倍以上は速くなるような気がします。

回答No.2

ki073

No.1です。
二重、三重のループになっていますが、内側のループの変数はprivate宣言が必要なはずです。
No.1にリンクしているOpenMP入門(1)のp 20あたりを参考に

それと、コンパイラのオプション(-fopenmp など)と環境変数の設定(OMP_NUM_THREADS)も忘れていませんか？
スレッド生成のオーバーヘッドがありますので、場合によっては効果がないことがあります。

openMPはprivate指定を忘れると正常な結果がでませんので、openMPを使わない状態の結果との比較して、結果が正しいか確認が必要です。

確認ですが、
使用しているOS, CPUの種類(マルチCPUの場合は全体のcore数も)、使用コンパイラ、
それとlwixの値はどの程度でしょうか？

お礼

ありごとうございます。
二重ループの場合、
#pragma omp parallel for private(lwjy)
を書き込めば良いのはわかったんですが、三重の場合はどのようにすればよいのでしょうか？

またOSはWindows 7、CPUはIntelのコア2、コンパイラはMicrosoft Visual C++ 2010 Expressです。

回答No.1

ki073

まずOpenMP入門(1)と(2)を見てください。
http://www2.cc.kyushu-u.ac.jp/scp/system/library/OpenMP/OpenMP.html

要するに、ループの順序が入れ替わっても問題ないところを並列化すればよいのです。
一番外の
for( kt=0; kt<=it ; kt++)
は順序を入れ替えると駄目ですが、その中のループは入れ替え可能です。
出来るだけ外側の方がスレッドの生成が減らせますので、
for( jx=0; jx<=ix; jx++ )
のループを並列化させます。

多分速くすることが目的でしょうから、もう少し改良できるところがあります。
電界(Ez)の計算の中の、
sigmax[jx][jy][jz]*dt/(2.0*epsirx[jx][jy][jz])
は何度かでてきますので、1つに
それと、励振の
ez[jx][80][jz] = sin( 2.0*pi*(kt%100)/100 )*sin(pi*float(jx)/float(ix));
は値が変わらないので、ループの外側であらかじめ計算しておいて、この位置で代入すれば計算量が減らせます。
これだけの範囲では分からないのですが、ループの外側に出せる計算があればそれにより計算量を減らせます。

細かいことですが、
/dx, /dyなどの定数の割り算がありますが、割り算はかけ算などに比べて計算時間がかなり必要なので、逆数をあらかじめ計算しておき、その値を掛けることで多少は速くなります。ただしこの変換はコンパイラの最適化で自動的に行われる可能性もあります。

お礼

解答ありがとうございます。
とてもためになりました。
高速化のために改良できるように頑張ってみます。
ちなみに実際に適用したいプログラムのソースコードはこんな感じです。

for( lwkt=0; lwkt<=lwit ; lwkt++)
{ /*　※※※forの括弧開始※※※　*/

/*********** 1ステップ前の電磁界[jt-1]の値の定式 *******/
#pragma omp for
for( lwjx=0; lwjx<=lwix; lwjx++ )
for( lwjy=0; lwjy<=lwiy; lwjy++ )
for( lwjz=0; lwjz<=lwiz; lwjz++ )
{
lwexg[lwjx][lwjy][lwjz] = lwex[lwjx][lwjy][lwjz];
lweyg[lwjx][lwjy][lwjz] = lwey[lwjx][lwjy][lwjz];
lwezg[lwjx][lwjy][lwjz] = lwez[lwjx][lwjy][lwjz];
}

/*********** 光波磁界の計算 ***********/

#pragma omp for
for( lwjx=0; lwjx<=lwix; lwjx++ )
for( lwjy=0; lwjy<=lwiy-1; lwjy++ )
for( lwjz=0; lwjz<=lwiz-1; lwjz++ )
{
long double lwdey = (lwey[lwjx][lwjy][lwjz+1]-lwey[lwjx][lwjy][lwjz])/lwdz;
long double lwdez = (lwez[lwjx][lwjy+1][lwjz]-lwez[lwjx][lwjy][lwjz])/lwdy;
lwhx[lwjx][lwjy][lwjz] = lwhx[lwjx][lwjy][lwjz]
+ lwdt/lwamux[lwjx][lwjy][lwjz] * (lwdey-lwdez);
}
#pragma omp for
for( lwjx=0; lwjx<=lwix-1; lwjx++ )
for( lwjy=0; lwjy<=lwiy; lwjy++ )
for( lwjz=0; lwjz<=lwiz-1; lwjz++ )
{
long double lwdez = (lwez[lwjx+1][lwjy][lwjz]-lwez[lwjx][lwjy][lwjz])/lwdx;
long double lwdex = (lwex[lwjx][lwjy][lwjz+1]-lwex[lwjx][lwjy][lwjz])/lwdz;
lwhy[lwjx][lwjy][lwjz] = lwhy[lwjx][lwjy][lwjz]
+ lwdt/lwamuy[lwjx][lwjy][lwjz]*(lwdez-lwdex);
}
#pragma omp for
for( lwjx=0; lwjx<=lwix-1; lwjx++ )
for( lwjy=0; lwjy<=lwiy-1; lwjy++ )
for( lwjz=0; lwjz<=lwiz; lwjz++ )
{
long double lwdex = (lwex[lwjx][lwjy+1][lwjz]-lwex[lwjx][lwjy][lwjz])/lwdy;
long double lwdey = (lwey[lwjx+1][lwjy][lwjz]-lwey[lwjx][lwjy][lwjz])/lwdx;
lwhz[lwjx][lwjy][lwjz] = lwhz[lwjx][lwjy][lwjz]
+ lwdt/lwamuz[lwjx][lwjy][lwjz]*(lwdex-lwdey);
}

/*********** 光波励振 ***********/
long double lwaz = sin( 0.04*pi*(lwkt%50));

//******* 領域II *******
#pragma omp for
for( lwjx=0; lwjx<lwkx1; lwjx++ )
for( lwjz=lwkz1; lwjz<=lwkz2; lwjz++ )
{
lwez[lwjx][lwky0][lwjz]
= ( cos(lwux)*cos((lwuz/lwtz)*(lwkz0-lwjz)*lwdz-lwfhi)
/exp((lwwx/lwtx)*((lwkx0-lwjx)*lwdx-lwtx)) )*lwaz;
}
//******* 領域III *******
#pragma omp for
for( lwjx=lwkx2+1; lwjx<=lwix; lwjx++ )
for( lwjz=lwkz1; lwjz<=lwkz2; lwjz++ )
{
lwez[lwjx][lwky0][lwjz]
= ( cos(lwux)*cos((lwuz/lwtz)*(lwkz0-lwjz)*lwdz-lwfhi)
/exp((lwwx/lwtx)*((lwjx-lwkx0)*lwdx-lwtx)) )*lwaz;
}
//******* 領域IV *******
#pragma omp for
for( lwjx=lwkx1; lwjx<=lwkx2; lwjx++ )
for( lwjz=0; lwjz<lwkz1; lwjz++ )
{
lwez[lwjx][lwky0][lwjz]
= ( (cos((lwux/lwtx)*(lwkx0-lwjx)*lwdx)*cos(lwuz-lwfhi))
/(exp((lww2z/lwtz)*((lwkz0-lwjz)*lwdz-lwtz))) )*lwaz;
}
//******* 領域V *******
#pragma omp for
for( lwjx=lwkx1; lwjx<=lwkx2; lwjx++ )
for( lwjz=lwkz2+1; lwjz<=lwiz; lwjz++ )
{
lwez[lwjx][lwky0][lwjz]
= ( (cos((lwux/lwtx)*(lwkx0-lwjx)*lwdx)*cos((lwuz/lwtz)*(lwkz0-lwkz2)*lwdz-lwfhi))
/(exp((lww0z/lwtz)*((lwjz-lwkz0)*lwdz-lwtz))) )*lwaz;
}
//******* 領域I *******
#pragma omp for
for( lwjx=lwkx1; lwjx<=lwkx2; lwjx++ )
for( lwjz=lwkz1; lwjz<=lwkz2; lwjz++ )
{
lwez[lwjx][lwky0][lwjz]
= cos((lwux/lwtx)*(lwkx0-lwjx)*lwdx)*cos((lwuz/lwtz)*(lwkz0-lwjz)*lwdz-lwfhi)*lwaz;
}

//******* 領域II-V *******
#pragma omp for
for( lwjx=0; lwjx<lwkx1; lwjx++ )
for( lwjz=lwkz2+1; lwjz<=lwiz; lwjz++ )
{
lwez[lwjx][lwky0][lwjz]
= ( (cos(lwux)*cos((lwuz/lwtz)*(lwkz0-lwkz2)*lwdz-lwfhi))
/(exp((lwwx/lwtx)*((lwkx0-lwjx)*lwdx-lwtx))*exp((lww0z/lwtz)*((lwjz-lwkz0)*lwdz-lwtz))) )*lwaz;
}
//******* 領域II-IV *******
#pragma omp for
for( lwjx=0; lwjx<lwkx1; lwjx++ )
for( lwjz=0; lwjz<lwkz1; lwjz++ )
{
lwez[lwjx][lwky0][lwjz]
= ( (cos(lwux)*cos(lwuz-lwfhi))
/(exp((lwwx/lwtx)*((lwkx0-lwjx)*lwdx-lwtx))*exp((lww2z/lwtz)*((lwkz0-lwjz)*lwdz-lwtz))) )*lwaz;
}
・
・
・
この後電界の計算が同じように書かれています。
このようにfor文の前に#pragma omp forを導入したのですがあまり短縮はできませんでした。