Mur吸収境界条件の準備作業をリスト3-3-1に示します。
この関数は引数をmode=0/1として2回呼ばれます。
最初のmode=0で作業配列の大きさを求め、配列をmallocします。
配列の大きさはNyNzのオーダーです。
次のmode=1で係数の作業配列を代入します。
なお、3,4,7行目はMPIによる領域分割も考慮しています。
領域分割をしないときはiMin=0, iMax=Nxです。
本コードは冗長な部分もありますが、準備作業は計算時間に影響しないので、
並列計算も考慮して読みやすさを優先しています。
リスト3-3-1 Mur吸収境界条件の準備作業(Hx成分、setupMurHx.c)
1 void setupMurHx(int mode)
2 {
3 int64_t num = 0;
4 for (int i = iMin - 0; i <= iMax; i++) {
5 for (int j = jMin - 1; j <= jMax; j++) {
6 for (int k = kMin - 1; k <= kMax; k++) {
7 if ((((j < 0) || (j >= Ny)) && (k >= 0) && (k < Nz)) ||
8 (((k < 0) || (k >= Nz)) && (j >= 0) && (j < Ny))) {
9 if (mode == 1) {
10 fMurHx[num].i = i;
11 fMurHx[num].j = j;
12 fMurHx[num].k = k;
13 unsigned char m = 0;
14 double d = 0;
15 int i1 = 0, j1 = 0, k1 = 0;
16 if (j < 0) {
17 m = IEZ(i, j + 1, k );
18 d = Yn[1] - Yn[0];
19 i1 = i;
20 j1 = j + 1;
21 k1 = k;
22 }
23 else if (j >= Ny) {
24 m = IEZ(i, j, k );
25 d = Yn[Ny] - Yn[Ny - 1];
26 i1 = i;
27 j1 = j - 1;
28 k1 = k;
29 }
30 else if (k < 0) {
31 m = IEY(i, j, k + 1);
32 d = Zn[1] - Zn[0];
33 i1 = i;
34 j1 = j;
35 k1 = k + 1;
36 }
37 else if (k >= Nz) {
38 m = IEY(i, j, k );
39 d = Zn[Nz] - Zn[Nz - 1];
40 i1 = i;
41 j1 = j;
42 k1 = k - 1;
43 }
44 fMurHx[num].g = factorMur(d, m);
45 fMurHx[num].i1 = i1;
46 fMurHx[num].j1 = j1;
47 fMurHx[num].k1 = k1;
48 }
49 num++;
50 }
51 }
52 }
53 }
54
55 // array size
56 if (mode == 0) {
57 numMurHx = num;
58 }
59 }
リスト3-3-2 Mur吸収境界条件の係数計算部(setup.c)
1 float factorMur(double d, int m)
2 {
3 if (m != PEC) {
4 double vdt = (C * Dt) / sqrt(Material[m].epsr * Material[m].amur);
5 return (float)((vdt - d) / (vdt + d));
6 }
7 else {
8 return -1;
9 }
10 }
反復計算の中で呼ばれるMur吸収境界条件をリスト3-3-3に示します。
12行目で境界の半セル内側の磁界を次回の計算のために保存しています。
リスト3-3-3 Mur吸収境界条件(Hx成分、murHx.c)
1 void murHx(void)
2 {
3 for (int64_t n = 0; n < numMurHx; n++) {
4 int i = fMurHx[n].i;
5 int j = fMurHx[n].j;
6 int k = fMurHx[n].k;
7 int i1 = fMurHx[n].i1;
8 int j1 = fMurHx[n].j1;
9 int k1 = fMurHx[n].k1;
10 HX(i, j, k) = fMurHx[n].f
11 + fMurHx[n].g * (HX(i1, j1, k1) - HX(i, j, k));
12 fMurHx[n].f = HX(i1, j1, k1);
13 }
14 }