4 次精度 FDTD 法の並列計算による大規模電波伝搬解析
園田12

潤 1 ，小林敬生 2 ，佐藤源之 3

仙台電波工業高等専門学校

電子工学科

Korea Institute of Geoscience and Mineral Resources
3
東北大学 東北アジア研究センター

はじめに1

近年の 3G/4G 携帯電話等の無線通信の高周波化や、UWB (Ultra Wide Band) の民間利用開放などの無線通信システムの高度化に伴い、通信システム最適設計のための電波伝搬解析が重要視されており、FDTD(Finite-Difference Time-Domain) 法 [1][2]が時間領域の解析手法として広く用いられている。FDTD 法は Maxwell の方程式を時間と空間で差分化し、時間毎に解析領域内の電磁界を計算する方法であり、電磁界の過渡応答解析が可能なことが特長である。しかしながら、屋内外における高周波電波伝搬問題等の大規模解析を行う場合には、多くのメモリを必要とし計算時間がかかる問題がある。この計算時間の問題について、FDTD 法は領域分割型の手法であり並列計算に適していることから、様々な並列計算に関する研究が行われている [3]-[6]。
また、従来の FDTD 法では空間を 2 次精度で差分化するため、大規模な領域の問題を解析するためにセルサイズを大きくした場合には、位相誤差が生じ伝搬距離が長くなるにつれ誤差が蓄積し増大する問題がある。この問題に対して、空間をより高次の4 次精度で差分化する FDTD(2,4) 法等が提案されている [7]-[10] が、FDTD(2,4) 法を用いた大規模伝搬解析に関する研究は、計算量や使用メモリが大きくなること等の理由からこれまでほとんど行われていない。これに対して我々は、PC クラスタを用いたFDTD(2,4) 法の並列計算により室内電波伝搬解析を行っている [11] が、FDTD(2,4)
法の並列計算では従来の FDTD に比べ通信量が 3 倍になるため、高効率な並列計算が行えていなかった。
そこで本稿では、PC クラスタ等の分散メモリ型並列計算機において FDTD(2,4) 法の並列計算を行う場合に問題となる、通信量増加による計算効率の低下を解決するために、共有メモリ型のスーパーコンピュータ SX-7 を用いた FDTD(2,4) 法による大規模電波伝搬の高速・高効率な解析について報告する。

2FDTD 法の大規模電波伝搬解析における数値分散誤差

2.1

FDTD 法による大規模電波伝搬解析の例

従来の FDTD 法を用いて大規模領域の解析を行う場合、数値分散誤差を軽減するためにはセルサイズを小さくする必要があり使用メモリが増大する。この例として、

1

A
Tr
B

yz

C
Human body

PEC

Concrete W=2.7m=48 λx
L = 30.0m = 530 λ
図 1: FDTD 法による室内 UWB パルス伝搬の解析モデル
表 1: 室内 UWB 伝搬解析における解析パラメータ

analysis region

L = 30.0 m, W = 2.7 m

concretehuman bodyεrc = 4.0 − j0.2elliptical cylinder (1.0m × 0.3 m) , εrh = 30.0 − j1.5

UWB pulse
space increment ∆stime increment ∆ttotal cells
required memoryτ = 0.16 ns (fmax =5.3 GHz, λ0 =0.06 m)∆s = 0.001 m (∆s = λ0 /57 = λg /10)
∆t = 0.02 ns
30000 × 2700 × 63.89 Gbyte

図 1 に示す室内廊下 (L = 30.0 m × W = 2.7 m) における UWB パルスの電波伝搬解析を示す。表 1 に解析パラメータを示す。図 1 の解析モデルにおいて、壁を厚さ0.1m のコンクリートとし、コンクリートの比誘電率 εrc を εrc = 4.0 − j0.2 とした。
扉と壁中の鉄筋を完全導体（PEC: perfectly electrically conductor）、室内は自由空間とし、最右端に Mur の吸収境界条件を用いた。また、人体のモデルとして、波源Tr から 5m の位置に比誘電率 εrh = 30.0 − j1.5 の楕円柱 (1.0m × 0.3m) A、B を 2体、A、B から 10m 離れた位置に C を設定した。波源 Tr を扉から 1m の位置に設定し、半値幅 0.16 ns のガウシアンパルスを電流源とした。ここで、ガウシアンパルスの最大周波数 fmax を-10 dB の周波数で定義すれば、fmax = 5.29 GHz となり、解析領域を波長で表せば、530λ × 48λ となる。FDTD 法におけるセルサイズ ∆s は、最大周波数 fmax = 5.29 GHz と比誘電率 εr = 30.0 の人体モデルにおいて λg /10 となるように ∆s = 1 mm とし、時間ステップ ∆t は CFL 安定条件から ∆t = 0.02 nsとした。これより、図 1 の解析モデルの解析には、30000 × 2700 × 6 個の配列を使用し、総メモリ約 4 Gbyte を要することになり、PC 1 台では計算が困難になる。

図 2 に UWB パルス送信後 20 ns、60 ns 後の室内廊下における電界分布 |Ez | を示す。図 2(a)、(b) は人体が存在する場合、(c) は人体が存在しない場合である。図 2(b)、(c) より、人体が存在する場合には電界分布 |Ez | が変化することが確認できる。
2

(a)

(b)

(c)
図 2: FDTD 法による室内廊下における UWB パルスの伝搬解析
布）(a) 20 ns 後

2.2(b) 60 ns 後(c) 60 ns 後（人体なし）（|Ez | の空間分

FDTD 法の数値分散誤差

ここで、FDTD 法を用いた大規模解析における数値分散誤差を議論するために、分散関係式を求め長距離伝搬による数値分散誤差を計算する。
従来の FDTD 法では、マクスウェルの方程式を式 (1) のような空間および時間の 2次精度中心差分で差分方程式に変換する。ここで、∆x は x 方向の空間ステップ（セルサイズ）である。

F n (i + 21 , j, k) − F n (i − 12 , j, k)∂F (i, j, k, t)=+ O(∆x2 )
∂x
∆x(1)

従来の FDTD 法では空間を式 (1) の 2 次精度で中心差分するため、波の伝搬速度3が方向によって異なり、位相誤差や分散誤差と呼ばれる誤差が生じる [2]。分散誤差を低減する方法として、テーラー展開によって得られる式 (1) の中心差分式を、より高次の項まで考慮する高次 FDTD 法が提案されている [7]。例えば、空間 4 次精度のFDTD(2,4) 法では式 (2) の中心差分式を用いる。
∂F (i, j, k, t)
9 F n (i + 12 , j, k) − F n (i − 12 , j, k)=∂x8∆x
1 F n (i + 32 , j, k) − F n (i − 32 , j, k)−+ O(∆x4 )
(2)
24
∆x
FDTD 法および FDTD(2,4) 法の分散関係式は、それぞれ 式 (3)、(4) のように求めることができる。
Ã!#2 "
Ã
!#2µ·
¶¸2 "
ω∆tk̃x ∆xk̃y ∆y1
11
sin
sin+
sin
=c∆t
2∆x2
∆y2"
Ã
!#21k̃z ∆z+sin(3)∆z2·1sinc∆tµω∆t2¶¸2"

ÃÃ!
Ã! !#2k̃x ∆x
1
1 9
3k̃x ∆x
=sin
−sin
∆x 8
2242"
ÃÃ!
Ã! !#2
11 9k̃y ∆y
3k̃y ∆y−+
sin
sin
∆y 8
2242"
ÃÃ!
Ã! !#2
1 9k̃z ∆z
13k̃z ∆z
+
sin
−sin
∆z 8
2242(4)

ここで、k̃ は数値的波数であり、θ を z 軸とのなす角、φ を x 軸とのなす角とすると、k̃x = k̃ sin θ cos φ、k̃y = k̃ sin θ sin φ、k̃z = k̃ cos θ で表される。
図 3 に式 (3)、(4) の分散関係式から得られる数値的波数 k̃ と、式 (5) より求めた

FDTD 法および FDTD(2,4) 法の数値分散誤差特性を示す。
¶µλ−1eφ = 2 πλ̃(5)図 3 は、セルサイズ（1 波長あたりのセルの分割数）による分散誤差特性であり、セルサイズを ∆x = ∆y = ∆z 、クーラン数 v ∆t/∆s を 0.3 、伝搬方向を θ = 90◦ 、φ =0◦ とした結果である。図 3 より、同じセルサイズでは FDTD(2,4) 法は従来の FDTD法に比べおおよそ一桁精度がよいことが分かる。また、図 3 より、FDTD(2,4) 法でセルサイズ ∆s を λ/10 とした場合、従来の FDTD 法を用いて FDTD(2,4) 法と同程度の精度で解析を行うには、セルサイズ ∆s を λ/45 にする必要があり、総セル数はFDTD(2,4) 法の約 90 倍必要になる。このことから、FDTD(2,4) 法は従来の FDTD法に比べ大規模な電波伝搬解析を低コストで解析できることが分かる。

4

absolute dispersion error [deg./λ]100

conventional FDTD
FDTD(2,4)

1010.10.010.001

020406080100cells/λ
図 3: FDTD および FDTD(2,4) 法の数値分散誤差

2.3

長距離伝搬における数値分散誤差の実際

図 4 に図 3 の結果から計算した伝搬距離の増加による数値分散誤差の蓄積を示す。
図 4 より、一般に知られているようにセルサイズを 波長 λ の 1/10 とすると、従来のFDTD 法では 1 波長伝搬あたり 6◦ の数値分散誤差が生じ、30λ の伝搬により 167◦の数値分散誤差を生じ位相がほぼ反転する。一方、FDTD(2,4) 法では 1 波長伝搬あたり 0.3◦ の数値分散誤差が生じ、30λ の伝搬においても 8.3◦ の数値分散誤差が生じることがわかる。
図 5 に実際の FDTD 法の計算において生じる数値分散誤差を示す。図 5 は、FDTD法および FDTD(2,4) 法で計算した微小ダイポールアンテナから放射された f = 1 GHzの正弦波の伝搬であり、セルサイズ ∆s = λ/10 とした場合の 30λ の観測点における電界 Ez の時間応答波形である。ここで、解析解は文献 [12] の微小ダイポールアンテナの放射界（時間領域）を用いた。図 4 において、伝搬距離に比例した数値分散誤差が生じ、30λ の伝搬により従来の FDTD 法では位相がほぼ反転することを示したが、図 5 より実際の計算においてもこの数値分散誤差が確認できる。

3

空間 4 次精度 FDTD 法の並列計算による大規模電波伝搬解析

3.1

FDTD(2,4) 法の並列計算

FDTD(2,4) 法を用いて大規模な領域の電波伝搬解析を行う場合には、より多くのメモリを必要とし、また、高速に計算するためには、大容量のメモリおよび複数のCPU による並列計算が必要になる。ここでは、複数台の計算機をネットワーク接続
510000

conventional FDTD, ∆s=λ/10
conventional FDTD, ∆s=λ/20
FDTD(2,4), ∆s=λ/10
FDTD(2,4), ∆s=λ/20

absolute dispersion error [deg.]

1000

1001010.1

λ10 λ100 λpropagation distance1000 λ図 4: FDTD および (2,4) 法における伝搬距離による数値分散誤差の蓄積

3.0

conventional FDTD
FDTD(2,4)
exact

electric field Ey [V/m]

2.0

1.0

0.0

-1.0-2.0-3.0
32.0

32.5

33.0

33.5

34.0
time [ns]34.5

35.035.536.0

図 5: 伝搬距離 30λ で生じる数値分散誤差の実際 (∆s = λ/10)

6analysis region

p-1

pp+1

NN

N/nz

yx

Hy
Hz

EyEz

HyHz

Ey
Ez

FDTD

EyEz
Hy x2
Hz x2

Ey x2
Ez x2
HyHz

Ey
Ez
Hy x2
Hz x2

Ey x2
Ez x2
HyHz

FDTD(2,4)

図 6: FDTD 法の並列計算における計算機間の電磁界の送受信
表 2: FDTD 法および FDTD(2,4) 法の並列計算の計算コスト

conventional FDTD
FDTD(2,4)

calculation [flop/(step · cell)]
E
H
total

communication [/step]
E
H
total

75
102

2N 26N 2214896
150

2N 26N 24N 212N 2する分散メモリ型並列計算環境である PC クラスタを用いた FDTD(2,4) 法の並列計算について示す。

PC クラスタのような分散メモリ型並列計算機を用いて並列計算を行う場合には、各 PC 間で相互にデータ送受信が必要になる。ここで、図 6 に示すような N × N × Nの解析領域を n 台の PC で並列計算することを考え、並列計算の領域分割法としてyz 平面による 1 次元分割を考える。FDTD 法の並列計算では、電界 Ey 、Ez および電界 Hy 、Hz 計 4N 2 個のデータ送受信が必要であるが、FDTD(2,4) 法の並列計算
では、図 6 に示すように PC 間で計 12N 2 個の電磁界の送受信が必要になる。以上のように、FDTD(2,4) 法の並列計算では、従来の FDTD 法に比べデータ送受信量が3 倍増加する。
表 2 に FDTD 法および FDTD(2,4) 法の並列計算における計算コストを示す。計算量については、時間 1 ステップ 1 セルあたりの浮動小数点演算数 flop/(step · cell)で考えると、式 (1)、(2) から導出される差分方程式より FDTD(2,4) 法は FDTD 法に比べ約 1.6 倍増加する。通信量については、時間 1 ステップあたりに送受信するセル数を考えれば、FDTD(2,4) 法は FDTD 法に比べ約 3 倍増加する。表 2 に示すよう

7

に FDTD(2,4) 法の並列計算では、従来の FDTD 法に比べ計算量や通信量がともに増加するが、同程度の精度を得るのに必要なセル数は前述したしたように約 1/90 に削減できるので、全体としての計算量は約 1/60 程度に減少できる。しかしながら、通信量増加が並列計算の効率に大きな影響を及ぼすので、高速なデータ送受信が必要である。

3.2

PC クラスタにおける並列計算特性

表 3 に示す市販の計算機 (CPU: Pentium 4 3.0 GHz 、メモリ: 1.5 Gbyte) 16 台を用いた PC クラスタによる FDTD(2,4) 法の並列計算特性を示す。計算機間のネットワークとして、100 Base-TX (NIC: RTL8139) を用いる。OS として、Linux (kernel2.4.20-8) を用いる。PC クラスタのような分散メモリ型並列計算機で計算機間の通信を行うライブラリとして MPI (Message Passing Interface) の実装である mpich-1.2.6を用いる。
図 7 に PC クラスタにおける従来の FDTD 法および FDTD(2,4) 法の並列計算の計算時間を示す。図 7 は、周囲を完全導体で囲まれた解析領域 160 × 160 × 160 セルの自由空間における電波伝搬を、FDTD 法および FDTD(2,4) 法を用いて計算した計算時間である。セルサイズ ∆s および時間ステップ ∆t は、それぞれ ∆s = 0.01 m、∆t = 0.01 ns、計算時間は 1000 ステップ (10 ns) とした。図 7 より、FDTD 法では計算機 8 台以上、FDTD(2,4) 法では 4 台以上で使用計算機台数を増加すると計算時間も増加する現象が見られる。これは、全計算量に占める通信量の比が大きくなるためであり、また、FDTD(2,4) 法は FDTD 法に比べ 3 倍の通信量となるので、台数を増加させるごとに通信量の影響が大きくなるためである。これは通信比を少なくすれば改善でき、解析領域が大きいほど並列計算の効果が得られる。この通信量増加による計算時間の増大は、ギガビットネットワーク等の高速・高帯域のネットワーク環境を使用することにより改善される。
表 3: 計算環境

CPU, cachememoryPentium4 3.0GHz,
1.5 GB

number of PCs
total memory

1624 GB

NIC

RealTek RTL8139 (100Mbps)

OS
message passing

Linux (kernel 2.4.20-8)
MPI (mpich-1.2.6)

compiler

mpicc (gcc version 3.2.2)

83500conventional FDTD
FDTD(2,4)

computational time [s]

3000
2500
2000
1500
1000
500
00246810121416

number of PCs

図 7: PC クラスタにおける FDTD 法および FDTD(2,4) 法の並列計算の計算時間

3.3

スーパーコンピュータ SX-7 との比較

図 7 に示したように、通信量の多い FDTD(2,4) 法の場合、PC クラスタ等の汎用ネットワークによる分散メモリ型並列計算機では高い計算効率が得られないことから、東北大学情報シナジーセンターのスーパーコンピュータ SX-7 を用いて FDTD(2,4)
法の計算を行う。表 4 にそれぞれの計算機の仕様を示す。スーパーコンピュータ SX-7
は、総 CPU 数が 240、メモリ 1920 Gbyte の高性能計算機である。本稿では 16 台のCPU を用い、この CPU 数におけるジョブクラスで使用可能なメモリは 128 Gbyte である。本稿で使用した PC クラスタは CPU 数が 16、メモリ 24 Gbyte 、100 Base-TXのネットワークを持つが、通信量の多いアプリケーションでは通信性能がボトルネックになる。
スーパーコンピュータと PC クラスタを用いた FDTD(2,4) 法の並列計算により、周囲を完全導体で囲まれた解析領域 160 × 160 × 160 セルの電波伝搬解析の時間 5000
ステップに要する計算時間を図 8 に示す。図 8 において、スーパーコンピュータは 16
表 4: 使用するスーパーコンピュータのジョブクラスと PC クラスタの仕様
super computer SX-7 PC cluster (Pentium4 3.0 GHz × 16)

number of CPU

1616

main memory
parallelize

128 Gbyte
auto (sxcc -Pauto )24 Gbyte
Message Passing (MPI, mpich-1.2.6)

9

100000

computational time [s]

10000

1000

10010

PC cluster

SX-7

conventional FDTD

PC cluster

SX-7

FDTD(2,4)

図 8: スーパーコンピュータと PC クラスタによる FDTD 並列計算の計算時間の比較

台の CPU を使用した場合の計算時間である。図 8 より、FDTD 法や FDTD(2,4) 法ともスーパーコンピュータは、PC クラスタに比べ高速に計算できることが分かる。
特に、通信量の多い FDTD(2,4) 法の計算においては、PC クラスタでは、スーパーコンピュータに比べ約 340 倍の計算時間を要している。また、スーパーコンピュータを用いた FDTD(2,4) 法の計算時間を見ると、FDTD 法の 1.53 倍程度の増加であり、表 2 に示した 1 ステップあたりの計算量の増加量 1.56 倍と同程度の値である。これはデータ通信が非常に高速であるためであり、通信量増加の影響をほとんど受けていないことが分かる。

4

まとめ
本稿では、FDTD(2,4) 法を用いて大規模電波伝搬解析を低コストで行うことを目的に、スーパーコンピュータ SX-7 による高速・高効率な計算について示した。大規
模な電波伝搬を FDTD 法で解析する場合には、数値分散により誤差が蓄積するため精度よく解析するためには、従来の手法ではセルサイズを小さくする必要があり計算コストが増大する問題があった。これに対し、空間 4 次精度の FDTD(2,4) 法では、同程度の計算精度を得るのに従来の FDTD 法に比べ、使用メモリはおおよそ 1/90、計算量は 1/50 程度に削減できる。FDTD(2,4) 法の計算を PC クラスタ等の分散メモリ型並列計算機で行う場合には、通信量増加のため高効率な計算は困難であるのに対し、スーパーコンピュータでは通信量増加の影響をほとんど受けることなく高速・高効率な計算が行えることを示した。

10謝辞
本研究は、東北大学情報シナジーセンターの大規模科学計算システムを利用して行われた。ここに感謝の意を表す。

参考文献[1] K. S. Yee, “Numerical Solution of Initial Boundary Value Problems Involving
Maxwell’s Equations in Isotropic Media,” IEEE Trans. Antennas and Propagat., 14, 2, pp. 302–307 (1966-5)[2] A. Taflove, Computational Electrodynamics:

The Finite-Difference Time-Domain Method Artech House, (1995)[3] V. Varadarajan and R. Mittra, “Finite-Difference Time-Domain (FDTD) Analysis Using Distributed Computing,” IEEE Trans. Microwave and Guided Wave
Lett., 4, 5, pp.144–145 (1994-5)[4] C. Guiffaut and K.Mahdjoubi, “A Parallel FDTD Algorithm Using the MPI
Library,” IEEE Antennas and Propagation Magazine, 43, 2, pp.94–103 (20014)[5] 打矢 匡, 柏 達也, “並列型スーパーコンピュータを用いた FDTD 並列計算,” 信学論 C, Vol. J84-C, pp.1122–1125 (2001-11)[6] 園田 潤, ヘテロ PC クラスタにおける FDTD 並列計算の自動最適負荷分散, 信
学論 B, J87-B, No. 5, pp. 760–764, (2004-5)[7] P. G. Petropoulos, “Phase Error Control for FD-TD Methods of Second andFourth Order Accuracy,” IEEE Trans. Antennas and Propagat., 42, 6, pp.
859-862 (1994-6)[8] M. F. Hadi and M. Piket-May, “A Modified FDTD (2, 4) Scheme for Modeling Electrically Large Structures with High-Phase Accuracy,” IEEE Trans.
Antennas and Propagat., 45, 2, pp. 254–264 (1997-2)[9] S. V. Georgakopoulos, R. A. Renaut, C. A. Balanis and C. R. Birtcher, “AHybrid Fourth-Order FDTD Utilizing a Second-Order FDTD Subgrid,” IEEE
Trans. Microwave and Wireless Components Lett., 11, 11, pp. 462–464 (200111)[10] K. L. Shlager and J. B. Schneider, “Comparison of the Dispersion Propertiesof Several Low-Dispersion Finite-Diffe rence Time-Domain Algorithms,” IEEE
Trans. Antennas and Propagat., 51, 3, pp. 642–653 (2003-3)11[11] 園田 潤, 佐藤源之, “高次 FDTD 法の並列計算による大規模電波伝搬解析と光電界センサを用いた室内 UWB パルス伝搬実験,” 電気学会電磁界理論研究会 EMT

05-25, pp.1–6, (2005-7)[12] Glenn S. Smith, An Introduction to Classical Electromagnetic Radiation. Cambridge University Press, 199712