PML吸収境界条件の準備作業をリスト3-4-1に示します。
この関数は引数をmode=0/1として2回呼ばれます。
最初のmode=0で作業配列の大きさを求め、配列をmallocします。
次のmode=1で係数の作業配列を代入します。
なお、8-10行目はMPIによる領域分割も考慮しています。
19-24行目で境界での物性値番号を取得しています。
本コードは冗長な部分もありますが、準備作業は計算時間に影響しないので読みやすさを優先しています。
リスト3-4-1 PML吸収境界条件の準備作業(Ex成分、setupPmlEx.c)
1 void setupPmlEx(int mode)
2 {
3 int lx = cPML.l;
4 int ly = cPML.l;
5 int lz = cPML.l;
6
7 int64_t num = 0;
8 for (int i = iMin - lx + 0; i < iMax + lx; i++) {
9 for (int j = jMin - ly + 1; j < jMax + ly; j++) {
10 for (int k = kMin - lz + 1; k < kMax + lz; k++) {
11 if ((i < 0) || (Nx <= i) ||
12 (j < 0) || (Ny < j) ||
13 (k < 0) || (Nz < k)) {
14 if (mode == 1) {
15 fPmlEx[num].i = i;
16 fPmlEx[num].j = j;
17 fPmlEx[num].k = k;
18 int m = 0;
19 if (i < 0) m = IEX(0, j, k );
20 else if (Nx <= i) m = IEX(Nx - 1, j, k );
21 else if (j < 0) m = IEX(i, 0, k );
22 else if (Ny < j) m = IEX(i, Ny, k );
23 else if (k < 0) m = IEX(i, j, 0 );
24 else if (Nz < k) m = IEX(i, j, Nz );
25 assert((m >= 0) && (m <= UCHAR_MAX));
26 fPmlEx[num].m = (unsigned char)m;
27 }
28 num++;
29 }
30 }
31 }
32 }
33
34 // array size
35 if (mode == 0) {
36 numPmlEx = num;
37 }
38 }
PML吸収境界条件による使用メモリー増加分は Nx, Ny, Nz >> L のとき次式で評価されます。
{(Nx+2L)(Ny+2L)(Nz+2L) - (Nx Ny Nz)} * 6 * (4+4+4+4+4+4+1) ～ 300 (NxNy + NyNz + NzNx) L バイト
例えばNx=Ny=Nz=100, L=5層のとき使用メモリー増加分は45MBになります。
これは吸収境界条件を除く使用メモリーが30MBであることからPMLにより使用メモリーが2.5倍になることを意味します。
ただし、この比はセル数が大きくなると小さくなります。
反復計算の中で呼ばれるPML吸収境界条件をリスト3-4-2に示します。
リスト3-4-2 PML吸収境界条件(Ex成分、pmlEx.c)
1 void pmlEx(void)
2 {
3 int ly = cPML.l;
4 int lz = cPML.l;
5
6 for (int64_t n = 0; n < numPmlEx; n++) {
7 int i = fPmlEx[n].i;
8 int j = fPmlEx[n].j;
9 int k = fPmlEx[n].k;
10 unsigned char m = fPmlEx[n].m;
11
12 float dhz = HZ(i, j, k) - HZ(i, j - 1, k);
13 float ry = (m != PEC) ? RYn[MIN(MAX(j, 0), Ny )] * rPmlE[m] : 0;
14 Exy[n] = (Exy[n] + (ry * dhz)) * gPmlYn[j + ly];
15
16 float dhy = HY(i, j, k) - HY(i, j, k - 1);
17 float rz = (m != PEC) ? RZn[MIN(MAX(k, 0), Nz )] * rPmlE[m] : 0;
18 Exz[n] = (Exz[n] - (rz * dhy)) * gPmlZn[k + lz];
19
20 EX(i, j, k) = Exy[n] + Exz[n];
21 }
22 }