3.2 電磁界の更新

3.2.1 必要メモリー

電磁界はEx/Ey/Ez/Hx/Hy/Hzの6成分があります。
これらはメモリー削減のためとGPUを考慮して単精度(float)とします。
2.3で述べたように電磁界を更新するにはこれ以外に、電界の1成分につきc₁,c₂の2個、磁界の1成分につきd₁,d₂の2個の係数が必要になります。
ここでは必要メモリーの削減のために、 Yee格子の電磁界6成分の各点での物性値番号をメモリーに保存し、係数c₁,c₂,d₁,d₂ は物性値番号の配列として間接的にアクセスします。
物性値の数を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(void)
     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			for (int k = 0; k <= Nz; k++) {
     7				const int64_t m = iEx[n];
     8				Ex[n] = C1[m] * Ex[n]
     9				      + C2[m] * (RYn[j] * (Hz[n] - Hz[n - Nj])
    10				               - RZn[k] * (Hy[n] - Hy[n - Nk]));
    11				n++;
    12			}
    13		}
    14		}
    15	}