3.2 電磁界の更新

3.2.1 必要メモリー

電磁界はEx/Ey/Ez/Hx/Hy/Hzの6成分があります。
これらはメモリー削減のためとGPUを考慮して単精度(float)とします。
2.3で述べたように電磁界を更新するにはこれ以外に、 電界の1成分につきc1,c2の2個、 磁界の1成分につきd1,d2の2個の係数が必要になります。
ここでは必要メモリーの削減のために、 Yee格子の電磁界6成分の各点での物性値番号をメモリーに保存し、 係数c1,c2,d1,d2 は物性値番号の配列として間接的にアクセスします。
物性値の数を256種類までに限定すると電磁界1成分に必要なメモリーは1バイトで済みます。
以上によりセルあたりの必要メモリーは6*(4+1)=30バイトとなります。
従って必要メモリーの合計は30NxNyNzバイトになります。
なお、物性値番号については、真空(物体が存在しない所)を0、 完全導体(PEC)を1と決め、誘電体と磁性体を2以上255以下とします。

3.2.2 電磁界変数の表現

電磁界変数は1次元配列で表現します。 配列の並びは外側からi,j,kとし、3重ループも外側からi,j,kとします。 配列のアドレスには次式のマクロを使用します。

#define NA(i,j,k) ((i)*Ni+(j)*Nj+(k)*Nk+N0)
ここで、Ni,Nj,Nk,N0は予め計算しておきます。
物性値番号も同様に表現します。

3.2.3 ソースコードの説明

電界のX成分を更新する関数はリスト3-2-1の通りです。
ここでアドレスを計算するマクロNAにおいてNk=1であることからkが1増えるとNAも1増えることを利用しています。 このようにすると一番内側のkに関するループが自動ベクトル化されやすくなります。
10,11行目のC1,C2は物性値で決まる係数です。
座標に関する係数(RYn,RZn等)は予め計算しておきます。
リスト3-2-1は非常に簡単なコードですが、 高速化のためのチューニング作業を考えると簡単であるほど望ましいです。

リスト3-2-1 Ex更新プログラム(updateEx.c)
     1	void updateEx()
     2	{
     3		for (int i = 0; i <  Nx; i++) {
     4		for (int j = 0; j <= Ny; j++) {
     5			int64_t n = NA(i, j, 0);
     6			int64_t n1 = n - Nj;
     7			int64_t n2 = n - Nk;
     8			for (int k = 0; k <= Nz; k++) {
     9				int m = iEx[n];
    10				Ex[n] = C1[m] * Ex[n]
    11				      + C2[m] * (RYn[j] * (Hz[n] - Hz[n1])
    12				               - RZn[k] * (Hy[n] - Hy[n2]));
    13				n++;
    14				n1++;
    15				n2++;
    16			}
    17		}
    18		}
    19	}