電磁界シミュレーターOpenFDTD利用講習会
主催：公益財団法人計算科学振興財団
共催：一般財団法人高度情報科学技術研究機構

株式会社EEM 大賀明夫
2022年3月2日

目次
第一部
1. OpenFDTDとは
2. FDTD法
3. OpenFDTDの高速化技術
4. OpenFDTDの計算例
第二部
5. OpenFDTDの使用法
6. FOCUSスパコンの使用法

2

1. OpenFDTDとは● 電磁界シミュレーターOpenFDTD[1]・電磁界シミュレーターをたくさんの人に使ってもらうた
めに公開した（社会的インフラ）・フリーソフト：無料で機能制限なし
・オープンソース（使用言語：C）
・動作環境：Windows（簡易GUIと実行プログラム付）またはLinux（ソースコードからコンパイルする）
・計算手法はFDTD法（時間領域差分法）、汎用的な電磁界シミュレーター
・各種の高速化技術を取り入れ大規模問題に使用することを想定している●OSSの利点・ライセンスの制限がない（いつでもどこでも何人でも使える）・自宅のPCにグラボを刺せばミニスパコン環境が実現、自宅でシミュレーション作業
・スパコン環境で自分でコンパイルして実行できる（商用ソフトでは難しい、将来は商用ソフト
もスパコンで課金利用？）
・必要な機能は自分で改変可能

・計算方法、使用方法の詳細についてはホームページで公開している
・2014年5月初版リリース
3

OpenFDTDの用途
Maxwell方程式で記述できる現象はすべて解析対象になる
地球規模からナノスケールまで（106m～10-9m）15桁に及ぶ●電波応用
1. 各種アンテナの設計
2. 携帯端末、携帯基地局アンテナ●光応用
1. 金銀ナノ粒子の近接場光学
2. フォトニック結晶光導波路

3. 放送用送信アンテナ

3. SPR（表面プラズモン共鳴）

4. ワイヤレス給電（電力伝送）
5. 屋内の電波伝搬解析（wifi、ローカル5G）

4. TERS（先端増強ラマン散乱顕微鏡）
5. SNOM（走査型近接場光学顕微鏡）

6. 交通機関のwifiサービス（新幹線、航空機、電車）
7. ミリ波、テラヘルツ波アンテナ（5G/6G）
8. レーダー断面積（軍事用、車載レーダー）

6. ナノセンサー（化学物質等）

9. 車載アンテナ（車車間通信、自動運転）
10. メタマテリアル
11. 人体の電波被爆評価
4

電磁界シミュレーターの歴史
FDTD法

モーメント法1966 K.S.Yee[2], FDTD法を最初に発表(2次元)1965 K.K.Mei[9], 数値計算法を最初に発表
1968 R.F.Harrington[10], モーメント法教科書

1975 A.Taflove[3], 3次元解析
1981 G.Mur[4], 吸収境界条件(開放領域)
1993 K.Kunz and R.J.Luebbers[5], FDTD法教科書(XFDTDへ)

1981 G.J.Burke and A.J.Poggio[11], NEC公開(OSS)
1990～ 商用ソフトが多数販売され今日に至る

1994 J.-P.Berenger[6], PML(高精度吸収境界条件)
1995 A.Taflove[7], FDTD法教科書
1998 宇野[8], FDTD法教科書
2000～ 商用ソフトが多数販売され今日に至る問題のサイズ（対数）モーメント法FDTD法
log(計算時間)

スパコン

モーメント法が有利

PC

FDTD法が有利

log(大きさ/波長)

2000

2010

2020

5

2. FDTD法
計算領域●FDTD法とは何か
・Finite Difference Time Domain Method（時間領域差分法）
・Maxwell方程式の微分形式を時間領域と周波数領域で離散化して電磁界の時間変化を求める手法
・現在の電磁界シミュレーションでは一番広く使われている手法
・有限の大きさの計算領域を考え多数のセルに分割する・得られた時間波形をフーリエ変換して周波数特性を求める・磁界から電流を計算しこれから入力インピーダンスを求める
・計算領域の境界面の電界と磁界から遠方の放射パターンが計算される
・完全導体、誘電体、磁性体を扱うことができる
・電磁界のほとんどの用途に対応できる・並列計算に適した計算手法であるために高速な計算機環境を生かすことができる

完全導体(線)

誘電体

給電点 ～

吸収境界条件完全導体(面)FDTD法の計算イメージ

6

Maxwell方程式とYee格子
Maxwell方程式∂E+σ e E∂t∂H−∇×E= μ+σ m H∂t∇ ×H=ϵ

■未知数
E:電界
H:磁界
■既知数
ε:誘電率
μ:透磁率
σe:導電率σm:導磁率●Maxwell方程式の特長
・2つの方程式にすべての情報が含まれている
・補助方程式、実験式、経験式など一切なし
・これを解けばすべての電磁気の現象が説明できるHx
Ez
Z

Ey
Y (i,j,k)
X

単位セル

Hy
Hz

Ex

●Yee格子・計算領域を多数のセルに分割する
・電界は稜線の中心の接線方向・磁界は面の中心の法線方向
・Maxwell方程式と相性がよい

7

FDTD法の離散化
En

Maxwell方程式を時間領域と空間領域で離散化する
時間的に電界Eと磁界Hを交互に計算する(蛙跳び法、leap frog法)
電界と磁界のX成分は以下のようになる
c1,c2,d1,d2は電気定数から計算される既知の場所の関数

En+1

E
時間t
H

nはタイムステップ

Hn-1/2

Y成分、Z成分についてはX→Y→Zと巡回する

Hn+1/2

時間幅については上限（安定性条件）がある（陽解法）

(

) (
)
1
1
1
1
1
1
1
1
+ c {H ( i+ , j+ , k )−H ( i+ , j− , k )}−c {H (i + , j , k + )− H ( i+ , j , k − )}
2
2
2
2
2
2
2
2
1
1
1
1
H ( i , j+ , k + )=d H
i
,
j+
,
k
+
( 2 2)
22
1
1
1
1−d {E ( i , j +1 , k + )−E (i , j , k + )}+ d {E ( i , j + , k +1)− E (i , j+ , k )}
2
2
2
2
1
1
E nx +1 i + , j , k =c 1 E nx i + , j , k
2
2

n+

2y

1
2

n+

z

2z

1
2

n−

x

n
z

n+

z

n+

2y

1
2

1

n
z

1
2

n+

y

1
2

y

1
2

x

2z

n
y

n
y

8

波源モデル・電磁界を発生させるためには何らかの波源が必要である
・波源として二つのモデルを考える(2)平面波入射モデル（散乱の計算）(1)給電点モデル（アンテナの計算）

Ei

F(t)

1.0

0.0

k
給電点 ～

パルス平面波

-1.0
0.0

0.5

1.0

t/tw

1.5

2.0

計算対象

給電点にパルス電圧を印加する（微分ガウスパルス）

9吸収境界条件

計算領域

計算対象

吸収境界条件：無限に広がる電磁波を有限の領域で正しく計算するには計算領域の境界に吸収境界条件を設定することが必要
以下の2通りがある(1)Mur一次 : 1層
(2)PML : 5～10層（メモリーと計算時間が少し増えるが精度が向上する）

10

計算手順

パルスのピーク

収束

収束状況：
・計算領域内の平均電界と平均磁界が十分小さくなったとき収束したとみなす（ピークの1/1000が目安）
・収束が不十分なときは計算精度が落ちるので常に十分収束していることを確認すること

11

計算出力(1) 近傍電磁界
近傍電磁界：計算領域内の電磁界（周波数の関数）
電磁界の時間波形をFourier変換する∞∫ E ( r , t )e− j ω t dt

E (r , ω)= 0∞∫ V in (t )e− j ω t dt
0∞

∫ H (r , t) e− j ω t dt

H ( r , ω)= 0∞∫ V in (t ) e− j ω t dt
012

計算出力(2) 入力インピーダンス
入力インピーダンスと反射係数
Vin : 給電点に加える電圧（既知数）
Iin : 給電点を流れる電流
Z0 : 給電線の特性インピーダンス（既知数、同軸線路では50Ω）

Z in (ω)=

V in (ω)
I in (ω)

Hy（入力インピーダンス）

Iz

Z in (ω)−Z 0
Γ (ω)=（反射係数）
Z in (ω)+Z 0

(

)

{ (

Hy

Hx

Z
Y

X

電流は磁界から計算される
Iz i , j ,k+

Hx

) (

)}

{ (

) (

1
1
1
1
1
1
1
1
1
=Δ y j H y i + , j , k + −H y i− , j , k +−Δ x i H x i , j + , k + − H x i , j − , k +
2
2
2
2
2
2
2
2
2

)}
13

計算出力(3) 遠方界の計算法

計算領域

E(r,θ,φ)遠方界：計算領域表面の電界と磁界の接線成分から積分で計算される

E(r)

− jkr

− jk eE {θ , ϕ } (r , θ , ϕ)=F { θ ,ϕ } (θ , ϕ)
4πrF θ (θ , ϕ)=η N θ (θ , ϕ)+ L ϕ (θ , ϕ)F ϕ (θ , ϕ)= η N ϕ (θ , ϕ)− Lθ (θ , ϕ)

dS

H(r)

N (θ , ϕ)=∫ J (r )exp ( j k r⋅^r ) dSSZ

L (θ , ϕ)=∫ M (r )exp ( j k r⋅^r ) dSr^

S

J ( r)= n×^ H (r )^ ×E (r )M (r )=− nϕ^μ0η= ϵ ≈120 π [Ω]
0

√

θ

r

θ^

k =ω √ ϵ 0 μ 0
φ

Y

X

極座標

14

計算出力(4) 遠方界特性
アンテナ：利得

平面波入射：散乱断面積21 4π r2
2
G (θ , ϕ)=|E θ (r , θ , ϕ)| +|E ϕ ( r , θ , ϕ)| }{P in 2 η
k2={|F θ (θ , ϕ)|2 +|F ϕ (θ , ϕ)|2 }8 π η P in

1
P in = Re(V in I *in )2

σ (θ , ϕ )=4 π r 2 {|Eθ ( r ,θ , ϕ )|2 +|Eϕ (r , θ , ϕ )|2 }2k0=|Fθ (θ , ϕ )|2 +|Fϕ (θ , ϕ )|2 }{4π

σ F=σ ( π −θ i , π +ϕ i )σ B=σ (θ i , ϕ i )
2π

σ T=

（前方散乱断面積）（後方散乱断面積：RCS）

π

1∫ d ϕ ∫ sin θ d θ σ (θ , ϕ )4π 00（全散乱断面積）15図形出力例

スミスチャート近傍界電界分布図

入力インピーダンス

遠方界パターン図

反射係数とVSWR

遠方界全方向パターン図

16

3. OpenFDTDの高速化技術

FXの仕様[12]2.2GHz×512/64×2(FMA)=70.4GFLOPS
70.4GFLOPS×48=3.3792TFLOPS

17

A64FX

18

OpenFDTDの高速化技術
OpenFDTDは大規模問題に適用することを想定して各種の高速化技術を取り入れている
OpenMP・マルチコアによるマルチスレッド
・共有メモリー環境

MPI
・複数ノード間の通信・分散メモリー環境
CPU

FX

CPU

CUDA
・グラフィックボードによる
GPGPU
GPU

併用可能

GPU

FX

ベクトル化
・SIMDベクトル演算(FXでは512ビッ
ト)
FX
19

領域境界

FX(A64FX)での高速化作業

通信方向(1)ベクトル化
・各コアは512ビットのSIMDベクトル演算機能を持っている
・コンパイラー(fcc)の自動ベクトル化機能を利用
・一番内側のループをベクトル化する・OpenFDTDは単精度計算であるから理想的には512/32=16倍速くなる

Y
プロセスP
X

(2)OpenMP並列化 [13]

Ey

Hy

Ez

Hz

プロセスP+1
Hx

・A64FXの48コア内で並列計算する（共有メモリー）・一番外側のループを並列化する(3)MPI並列化 [14]・計算領域をプロセス数で分割する（分散メモリー）・計算の主要部のfor文の前にOpenMPディレクティブ
を記入（反復計算ごとに通信が必要）

・コンパイラー(fcc)が並列コードを出力する・プロセス間の通信にはMPI関数を記述する
・コンパイラー(mpifcc)が並列化コードを出力する
20

OpenMPによる並列計算
1000■OpenMPのスレッド数を変えたときの計算時間（FX700、clangモード、ベンチマーク200）

計算時間 [ 秒 ]OpenMP
linear・OpenMPのスレッド数を増やすと計算時間が短縮される
・10スレッドを超えると速度向上が頭打ちになる

100

10
1

10

100

スレッド数

21

MPIとOpenMPの併用
OpenMPスレッド数

計算時間■MPIのプロセス数とOpenMPのスレッド数を

1

48407.2秒2

24

854.3秒（FX700、clangモード、ベンチマーク200）
・4×12のとき最も速く、他の組み合わせは非

3

16

18.6秒

4

12

7.9秒

常に遅い（CPUのブロック図参考）
・1×48が前ページの48スレッドの10倍遅いの
はMPIの実装に問題あり?

6

8

53.0秒

●以下では4×12を採用する [15]86137.0秒

124

188.6秒242

585.5秒48

1

1133.0秒

MPIプロセス数

変えたときの計算時間（積は48で一定）

22

複数ノードの計算時間
1000■ノード数を1～6と変えたときの計算時間

ベンチマーク 300
ベンチマーク 400
ベンチマーク 500
ベンチマーク 600（FX700、clangモード、各ノードでは4×12）●ノード数が増えると計算時間が短縮される

計算時間 [ 秒 ]

100

10

1
1

ノード数

10

23

4. OpenFDTDの計算例
4周波数共用アンテナ [16]・スマホの小型化、電波干渉低減

802.11bLTE

LTE-advanced
802.11aアンテナ形状反射係数周波数特性(1～6GHz)

全方向放射パターン(2.45GHz)
・指向性は低い
24

車載アンテナアンテナ(垂直モノポール、右端)

後

前

放射パターン(6GHz)

電界分布(6GHz、中心面)
左

右

25

新幹線

電界分布(2.45GHz, 床上80cm)

アクセスポイント(AP)
26

通勤電車

送信アンテナ

電界分布（2.45GHz、床上55cm）

27

旅客機

AP

AP

電界分布(1.5GHz)
28

オフィス
送信アンテナ

電界分布（5GHz、床上0.9m）

電界強度の頻度と累積確率分布

29

数値人体モデル

電界分布

頭部断面

Ei
k
平面波入射
(100MHz)

Brooksモデル

情報通信研究機構(NICT)

1mmセル

数値人体モデル(2mmセル、58組織)
http://emc.nict.go.jp/bio/data
人体組織誘電率DB(10Hz～100GHz)
http://niremf.ifac.cnr.it/tissprop/htmlclie/htmlclie.php
30

ロケット

ロケットの噴煙

送信アンテナ・固体燃料ロケットの噴煙はプラズマ状態なので地上との通信に支障が出る
・噴煙のプラズマパラメーターを元に電界分
布とアンテナの放射パターンを計算する
・実測値との比較から燃焼モデルについての知見が得られる・2次元モデルの計算：軸対称なので2次元モ

90m

噴煙（プラズマ）

デルでよい予測ができる

60m
電界分布[dB]

31

レーダー断面積（RCS）
RCS

散乱断面積パターン（Sバンド）●レーダー断面積（Radar Cross Section：RCS）の計算

Ei

・実機と実周波数（Xバンド）の計算を行うには多くのメモ
リーと計算時間が必要になるk（メモリー6TB：富岳200ノードで計算時間1日）
・さらに全方向のRCS（=後方散乱断面積）を計算するには入射方向を変えて繰り返し計算する必要がある
・電波吸収体の計算も可能（ステルス性）
戦闘機の電界分布32

メタマテリアル●メタマテリアルの歴史
V. G. Veselago [17], 1968. メタマテリアルを提唱
J. B. Pendry [18], 1996. 最初に実現

透磁率μ
ε<0
μ>0
・光領域の金属
・プラズマ
ε<0
μ<0
n<0
・メタマテリアル

D.R.Smith [19], 2000. 負の屈折率を実現

ε>0
μ>0
n>0
・通常の媒質
ε>0μ<0・超電導・特殊な磁性体●メタマテリアルの応用 [16][20][21][24]・超薄型アンテナ
・高機能アンテナ
誘電率ε・高性能電波吸収体
・メタサーフェス（電波伝搬制御）
・メタレンズ（負の屈折率）
・メタホログラム・バイオセンサー・クローキング（偽視）33メタマテリアル線路
L

メタマテリアル線路[20]Cstub inductor
interdigital capacitor

通常の線路（右手系）

via
C

マイクロストリップ線路
吸い込み口

湧き出し口

L

メタマテリアル線路（左手系）
進行方向

左手系区間

Ezの位相分布34

メタマテリアルアンテナ
メタループアンテナ[16]

C

誘電体基板
GND

2.25GHz (fN2)

2.55GHz (fN1)

L

3.10GHz (fT)

反射係数の周波数特性(2～4GHz)

3.40GHz (fH1)

3.80GHz (fH2) 35

メタサーフェス

入射方向●レフレクトアレー[21]・5Gで使用されるミリ波は波長が短いので回折(回り込み)が弱い
・市街地では電波伝搬を改善するためにビルの壁に反射方向を制反射方向御するシートを貼る

Snellの法則
誘電体基板
GND

スーパーセル36

金・銀1粒子σ T=( )2 π 5 ϵ r −1
3 ϵ r +2

2

D6
λ4

レイリー散乱公式(D<<λ)
1粒子に平面波が入射する計算モデル
1E-16
Au( 計算値 )Au( レイリー公式 )Ag( 計算値 )Ag( レイリー公式 )TCS[m*m]

1E-17

1E-18

1E-19

金と銀の複素比誘電率の波長特性[22]
1E-20
300

350

400

450

500

wavelength[nm]

550

600

650

700

37
金・銀1粒子の全散乱断面積の波長特性（D=20nm）

金・銀2粒子
金や銀のナノ粒子を2個近づけて並べ、これに平面波を入射するとギャップ部（～1nm）に非常に大きな光が発生する（最大100万倍）

金2粒子の電界分布

銀2粒子の電界分布

38

フォトニック結晶光導波路

光導波路

フォトニック結晶(n=3.4)

フォトニック結晶光導波路(2Dモデル)[23]通過帯域反射係数S11・透過係数S21の周波数特性

光応用で以下のケースのときは関連シミュレーターOpenTHFD（周波数領域差分法）推奨
・半無限の誘電体があるとき ・周囲媒質が空気以外のとき ・周期構造39

5. OpenFDTDの使用法
OpenFDTDはWindowsではGUI、LinuxではCUIで操作する [25]-[28]操作法の二つのモード●計算とポスト処理の分割処理●計算とポスト処理の一括処理

・計算が終わった後、ポスト処理の設定を変えて繰り返し行える

・計算の前にポスト処理の設定を行う
・半分のメモリーと計算時間で実行できる※ポスト処理：計算結果を図形出力すること※ev2/ev3ファイル：2次元/3次元図形データ（独自仕様）

40

入力データの作り方(1) テキストファイル編集 (CUI)・エディタを用いて右のようなテキストファイルを編集する
・複雑なデータを入力することが難しい(2) GUI (Graphical User Interface)・Windows環境ではGUIでデータを作成し計算する
・入力データをLinux環境に転送し計算することができる
・データ入力が容易であるが、データ量が増えると作業量も増える(3) データ作成ライブラリー使用

OpenFDTD 2 7
xmesh = -0.075 30 0.075
ymesh = -0.075 30 0.075
zmesh = -0.075 10 -0.025 11 0.025 10 0.075
geometry = 1 1 0 0 0 0 -0.025 0.025
feed = Z 0 0 0 1 0 50
frequency1 = 2e9 3e9 10
frequency2 = 3e9 3e9 0solver = 1000 100 1e-3
end
入力データ例(ダイポールアンテナ)

・OpenFDTDに同梱されている「データ作成ライブラリー」を使用してデータを作成するプログラムを作る
・プログラミングが必要であるが、大量のデータやパラメー
ターを変えながら繰り返し計算することに向いている41

GUI（Windows環境）
コマンドボタン

データ入力部（7個のタブ）●計算部(1) 全般(2) メッシュ(3) 物性値(4) 物体形状(5) 給電点・観測点●ポスト処理(6) 図形出力制御(1)(7) 図形出力制御(2)42(2) メッシュ・X/Y/Z方向のメッシュを作成する
・区間座標と区間分割数を交互に入力する・計算領域の大きさとメッシュ分割数は計
算精度と計算時間に影響を与える
・計算誤差はセルサイズに比例する
・セルサイズを半分にすると8倍のメモ
リーと16倍の計算時間が必要になる

43

(3) 物性値・集中定数●物性値
・物性値のテーブルを作成する
・比誘電率、導電率、比透磁率、導磁率
・空気とPEC(完全導体)は予め用意されているので入力不要
・物体形状で物性値を選択する●集中定数
・R/L/Cを入力する

44

(4) 物体形状
・物体形状を多数のユニットに分解して入力する（1ユニット=1ページ）
・|<,<,>,>|でユニット番号を変える
・[追加]/[挿入]/[削除]でユニット操作
・[物性値]を選択する
・[形状]と[向き]を選択する
・[座標]を入力する
・[名前]はデータの管理に使用する（オプション）
・右の図で入力した形状を確認する
・複数のユニットが重複する場所の物性
値はユニット番号が大きい方が優先される
45

(5) 波源・観測点
・アンテナでは[給電点]を入力する
・散乱問題では[平面波入射]を入力する
・[観測点]はSパラメーターの計算に
使用する

46

計算
・(1)～(5)のタブを入力した後、[計算]をクリックすると計算が行われる
・標準出力に右のように出力される・計算条件、計算結果の要約を確認する<<< OpenFDTD (CPU+OpenMP) Ver.2.7.2 >>>Thread = 8, Process = 1
Sun Dec 19 10:12:10 2021
Title = dipole antenna
Source = feed
Cells = 30 x 30 x 31 = 27900
No. of Materials = 3
No. of Geometries = 1
No. of Feeds= 1
No. of Points= 0
No. of Freq.s (1) = 3
No. of Freq.s (2) = 1
CPU Memory size
= 3 [MB]
Output filesize= 2 [MB]
ABC = Mur-1st
Dt[sec] = 9.3089e-12, Tw[sec] = 5.0800e-10, Tw/Dt = 54.572
Iterations = 1000, Convergence = 1.000e-03=== iteration start ===
step
<E>
<H>
0 0.000000 0.000000
100 4.957039 3.894715
200 0.958861 0.874313
300 0.188227 0.169974
400 0.040359 0.034269
500 0.006610 0.006112
600 0.001507 0.001271--- converged --=== input impedance ===feed #1 (Z0[ohm] = 50.00)
frequency[Hz] Rin[ohm]Xin[ohm]Gin[mS]
Bin[mS]
2.00000e+09
34.646-104.563
2.855
8.618
2.50000e+09
70.9490.435
14.094-0.0863.00000e+09
134.375
84.3575.338-3.351=== output files ===ofd.log, geom3d.htm, ofd.out=== cpu time [sec] ===part-1 :0.093
part-2 :0.009-------------------total :0.102=== normal end ===Sun Dec 19 10:12:10 2021

Ref[dB]-2.096-15.227-4.606

VSWR
8.328
1.419
3.860

47

(6) 図形出力制御(1)・計算が終了した後、図形出力制御を行い、[図形出力]→[図形出力表示]を行う・図形出力には2Dと3Dがある●図形出力制御(1)・時間特性
・周波数特性・遠方界面上・遠方界全方向48

(7) 図形出力制御(2)●図形出力制御(2)・近傍界線上・近傍界面上49

計算設定●処理モード・[分割]/[一括]を選択する（通常は[分割]）●CPU/GPU・CPU：スレッド数を指定する・GPU：既定値でよい●MPI
・CPU：PCクラスタで使用
・GPU：マルチGPUとPCクラスタで使用必要メモリー = 78 Nx Ny Nz / Np バイト （プロセス当たり、Nx, Ny, Nz : セル数、Np：MPIプロセス数）
計算時間 ∝ M Nx Ny Nz / (Np Nt) （M：タイムステップ数、Nt：スレッド数、比例定数は計算機に固有の値）
502D図形出力・2D図形出力ビューワーev2d.exe
が2D図形データev.ev2を読み込んで図形表示する
・上のボタンの意味はマウスを移動して確認する
・|<.<,>,>|：ページ操作
・++：全ページ合成
・>>,+,-：動画・クリップボードにコピー
・図形ファイルに保存

51

3D図形出力・3D図形出力ビューワーev3d.exe
が3D図形データev.ev3を読み込んで図形表示
する
・上のボタンの意味はマウスを移動して確認する
・|<.<,>,>|：ページ操作
・X,Y,Z：視点移動
・o,+,-：拡大・縮小
・クリップボードにコピー
・図形ファイルに保存
・マウス操作：移動、回転、拡大、縮小

52

データ作成ライブラリー#include "ofd_datalib.h"
int main(void){// initialize
ofd_init();// meshofd_xsection(2, -75e-3, +75e-3);
ofd_xdivision(1, 30);
ofd_ysection(2, -75e-3, +75e-3);
ofd_ydivision(1, 30);
ofd_zsection(4, -75e-3, -25e-3, +25e-3, +75e-3);
ofd_zdivision(3, 10, 11, 10);・右のようなプログラムを作成して実行するとOpenFDTDの入力データが出力される
・以下のようなときに使用する(1)データの数が多いとき// material
ofd_material(2.0, 0.0, 1.0, 0.0, "");

(2)形状を変数でパラメーター化したいとき（パラメーターを変更すると形状を一括で変形する// geometry
ofd_geometry(1, 1, 0e-3, 0e-3, 0e-3, 0e-3, -25e-3, +25e-3);

ことができる）// feedofd_feed('Z', 0e-3, 0e-3, 0e-3, 1, 0, 50);(3)入力データを変えながら繰り返し計算するバッ
チ処理を行う（system関数を使用する）※バッチ処理は多数の入力データを作成した後、ス// frequency
ofd_frequency1(2e9, 3e9, 10);
ofd_frequency2(3e9, 3e9, 0);

クリプトファイルで実行することもできる// solverofd_solver(1000, 100, 1e-3);// output
ofd_outdata("sample1.ofd");
return 0;}

53

6. FOCUSスパコンの使用法
ファイル構成：
OpenFDTD/
ofd : 実行プログラム (OpenMP対応)
ofd_mpi : 実行プログラム (MPI+OpenMP対応)
ofd_cuda : 実行プログラム (CUDA対応)
ofd_cuda_mpi : 実行プログラム (CUDA+MPI対応)
include/ : ヘッダーファイルフォルダ
src/ : ソースコートフォルダ (OpenMP対応)
mpi/ : ソースコートフォルダ (MPI＋OpenMP対応)
cuda/ : ソースコートフォルダ (CUDA対応)
cuda_mpi/ : ソースコートフォルダ (CUDA+MPI対応)
data/ : 入力データ用フォルダ/sample/ : サンプルデータ用フォルダ/benchmark/ : ベンチマークデータ用フォルダ
datalib/ : データ作成ライブラリ用ソースコードフォルダ
arf/ : 最適設計ツール用ソースコードフォルダ○実行プログラム（赤字）はビルドするとできるファイル（配布ファイルには含まれていない）○ソースコード（青字）はビルドするのに必要（作成済実行プログラムのみを使用するときは不要）

54

OpenFDTDのビルド方法$ cd src$ make -f Makefile_gcc clean ※1$ make -f Makefile_gcc
OpenMP版(ofd)が作成される$ cd ../mpi$ make -f Makefile_gcc clean$ make -f Makefile_gcc
MPI+OpenMP版(ofd_mpi)が作成される$ cd ../cuda$ make -f Makefile_linux clean$ make -f Makefile_linux
CUDA版(ofd_cuda)が作成される$ cd ../cuda_mpi$ make -f Makefile_linux clean$ make -f Makefile_linux
CUDA+MPI版(ofd_cuda_mpi)が作成される$ cd ..
※1 予め”mv Makefile_gcc Makefile”を行えば引数”-f Makefile_gcc”は不要※2 A64FXで富士通コンパイラーを使うときはMakefile_gccはMakefile_fccに置き換える※3 NEC SX Aurora TSUBASAでNECコンパイラーを使うときはMakefile_gccはMakefile_nccに置き換える
55

FOCUSとPCで共同作業する方法3通り(1) FOCUSで計算のみを行い、PCでポスト処理を行う
・PCのGUIで入力データを作成する
・入力データをFOCUSに転送する（入力データはKBのオーダー）・FOCUSで引数”-solver”をつけて計算する（ofd.outが出力される）・できたofd.outをPCに転送する・PCでポスト処理を（繰り返し）行う※一番操作性がよいがファイルサイズがGBを超える可能性があるofd.outを転送する必要がある(2) FOCUSで計算を行い、その後ポスト処理も行い、最後にPCで図形表示する
・PCのGUIで入力データを作成する
・入力データをFOCUSに転送する（入力データはKBのオーダー）・FOCUSで引数”-solver”をつけて計算する（ofd.outが出力される）・PCのGUIでポスト処理の設定を行い入力データをFOCUSに転送する
・FOCUSで引数”-post -ev”をつけてポスト処理を行う（引数”-ev”を省略するとev2d.htm,ev3d.htmが出力される）・できた図形データev.ev2,ev.ev3をPCに転送する（ファイルサイズはMBのオーダー）
・PCでビューワーev2d.exe,ev3d.exeを用いて図形表示する、必要なら3行上に戻る※ofd.outを転送する必要がないが手順が増える（問題のサイズが大きいとき使用する）

56

(3) FOCUSで計算とポスト処理を連続して行い、最後にPCで図形表示する
・PCのGUIで入力データを作成する
・入力データをFOCUSに転送する（入力データはKBのオーダー）・FOCUSで引数”-solver”も”-post”もつけずに計算する（ofd.outは出力されない）
・できた図形データev.ev2,ev.ev3をPCに転送する（ファイルサイズはMBのオーダー）
・PCでビューワーev2d.exe,ev3d.exeを用いて図形表示する※(1)(2)と比べると約半分のメモリーと計算時間ですむが、ポスト処理の設定を変えると再度計算から始める必要がある（ポスト処理の設定が確定しているときはこの方法でもよい）●参考・ofd.outはポスト処理に必要な計算結果をすべて含むバイナリファイル
・計算時間>>ポスト処理時間であるから計算結果をいったんofd.outに保存しておくとポスト処理を設定を変えながら繰り返し効率よく行うことができる
・ofd.outを名前を変えて保存しておき、後で名前を戻してポスト処理のみを行うことができる、ただしofd.outにはポスト処理設定は含まれていないのでポスト処理設定を含む入力データが同時に必要
・ofd.outはファイルサイズが大きいのでディスク容量の残量に注意
57

実行プログラムのファイル名
・ofd : Intel CPU 1ノード

（F/H/Zシステム）

・ofd_mpi : Intel CPU 複数ノード（F/H/Zシステム）※

・ofd_fcc : 富士通スパコン MPI非使用 （Xシステム）
・ofd_fcc_mpi : 富士通スパコン MPI使用 （Xシステム）※
・ofd_ncc : NECスパコン 1ノード （V/Wシステム）
・ofd_ncc_mpi : NECスパコン 複数ノード （V/Wシステム）※
・ofd_cuda : NVIDIA GPU 1ノード （Fシステム）
・ofd_cuda_mpi : NVIDIA GPU 複数ノード （Fシステム）※※MPI対応プログラムは1ノードでも計算可、計算時間（Xシステム以外）と総使用メモリーは非MPIと同じ

58

初回の作業
ログインサーバーにログイン後$ ssh ff （フロントエンドサーバーにログイン、パスワード不要）$ mkdir OpenFDTD（OpenFDTD用のディレクトリの作成）$ cd OpenFDTD$ cp -r /home1/gleg/share/OpenFDTD/data .
（サンプルデータのコピー）$ cp /home1/gleg/share/OpenFDTD/*.sh .
（計算実行スクリプトファイルのコピー※1）$ ln -s /home1/gleg/share/OpenFDTD/ofd . （実行プログラムのシンボリックリンク※2）$ ln -s /home1/gleg/share/OpenFDTD/ofd_mpi .
$ ln -s /home1/gleg/share/OpenFDTD/ofd_fcc .
$ ln -s /home1/gleg/share/OpenFDTD/ofd_fcc_mpi .
※1 自分で作るなら不要です※2 実ファイルをコピーしても構いません（”ln -s”をcpに変える）

59

計算実行スクリプト例※mpiexec_focusは一時的なパッチ、正しくはmpiexec●Xシステム1ノードで計算するとき（ofd_fcc.sh, MPI使用）#!/bin/bash#SBATCH -p x001h_lec

赤字は要編集（パーティション）#SBATCH -t 3
（計算時間上限、単位：分）#SBATCH -J ofd_fcc#SBATCH -o stdout.%J
FJ_PATH=/opt/FJSVstclanga/v1.1.0
export PATH=${FJ_PATH}/bin:${PATH}
export LD_LIBRARY_PATH=${FJ_PATH}/lib64:${LD_LIBRARY_PATH}
mpiexec_focus -n 4 ./ofd_fcc_mpi -n 12 -ev data/benchmark/benchmark100.ofd ※●Xシステム複数ノードで計算するとき（ofd_fcc_mpi.sh, MPI使用）#!/bin/bash#SBATCH -p x001h_lec#SBATCH -N 2

（パーティション）（ノード数）#SBATCH --ntasks-per-node=4（ノード当たり4プロセスが一番速い）#SBATCH -t 3
（計算時間上限、単位：分）#SBATCH -J fcc_mpi#SBATCH -o stdout.%J
FJ_PATH=/opt/FJSVstclanga/v1.1.0
export PATH=${FJ_PATH}/bin:${PATH}
export LD_LIBRARY_PATH=${FJ_PATH}/lib64:${LD_LIBRARY_PATH}
srun hostname > hostfile
mpiexec_focus -hostfile hostfile ./ofd_fcc_mpi -n 12 -ev data/benchmark/benchmark100.ofd ※60

ジョブの実行フロントエンドサーバーにログイン後$ cd OpenFDTD
スクリプトファイル（Xシステムのときはofd_fcc_mpi.sh）をエディタで編集する$ sinfo -s （稼働状況を確認する、NODESの”I”が空きノード数、下記は講習会用）
PARTITIONAVAIL TIMELIMITNODES(A/I/O/T) NODELIST
a001h_lecup
1:00:00
21/31/0/52 a[067-078,080,082-084,086-089,091095,098-101,103-117,120-127]
f001h_lecup
1:00:00
4/0/0/4 f[253-256]
f001h_p100_lec
up
1:00:00
0/0/0/0g001h_lecup
1:00:00
0/4/0/4 g[001-004]
h001h_lecup
1:00:00
2/6/0/8 h[001-005,030,034-035]
x001h_lecup
1:00:00
0/4/0/4 x[003-006]
g006m*
up
6:00
0/4/0/4 g[001-004]$ sbatch ofd_fcc_mpi.sh （ジョブ実行）$ squeue （ジョブの確認、STが現在の状態、R/PD/CG:実行中/開始待ち/終了処理中）
JOBID PARTITIONNAMEUSER ST
TIME NODES NODELIST(REASON)
1627223g006m
ofd uetk0001 R
0:03
1 g001$ scancel ジョブ番号（ジョブを取り消す、ジョブ番号はsqueueのJOBIDで確認）$ thismonth（今月の使用料を確認する）
61

出力ファイルの説明
・ofd.out : 計算結果のバイナリファイル、計算とポスト処理を分割したときのみ出力される
・ofd.log : 計算の標準出力、計算結果の確認に使用する、stdout.xxxxxxxと同じもの
・*.log : その他のログファイル、通常不要、詳しくはOpenFDTDのページの5.6参考
・ev.ev2 : 2次元図形データ、PCに転送してev2d.exeで図形表示する ※・ev.ev3 : 3次元図形データ、PCに転送してev3d.exeで図形表示する※

・geom.ev3 : 物体形状の3次元データであるがPCで確認済みなので不要※引数”-ev”をつけないときは代わりにev2d.htmとev3d.htmが出力される、HTMLファイルなのでPCに転送してブラウザで開く

62

練習問題

アンテナH=25
給電点

PEC

反射係数(2～4GHz)
誘電体εr=2

200
50

100

200

100

200
単位:mm

遠方界パターン(3GHz,X面)

63

文献[1] OpenFDTD, http://www.e-em.co.jp/OpenFDTD/[2] K.S.Yee, “Numerical solution of initial boundary value problem involving Maxwell’s equations in isotropic media”,
IEEE Trans. Antennas Propagatation, AP-14, no.3, pp.302-307, 1966.
[3] A.Taflove and M.E.Brodwin, “Numerical solution of steady-state electromagnetic scattering problems using the timedependent Maxwell’s equations,” IEEE Trans. Microwave Theory and Techniques, vol.23, pp.623-630, 1975.
[4] G.Mur, “Absorbing boundary conditions for the finite-difference approximation of the time domain electromagneticfield equations”, IEEE Trans. Electromag. Compat., EMC-23, no.11, pp.377-382, 1981.
[5] K.Kunz and R.J.Luebbers, Finite Difference Time Domain Method for Electromagnetics, CRC Press, 1993.
[6] J.-P.Berenger, "A perfectly matched layer for the absorption of electromagnetic waves", J. Computational Physics,
vol.114, pp.185-200, 1994.
[7] A.Taflove, Computational Electrodynamics, The Finite-Difference Time-Domain Method, Artech House, 1995.
[8] 宇野亨, “FDTD法による電磁界およびアンテナ解析,” コロナ社, 1998.
[9] K.K.Mei, "On the Integral Equations of Thin-Wire Antennas", IEEE Trans. Antennas & Propagation, AP-13, no.3,
pp.374-378, 1965.
[10] R.F.Harrington, Field computation by moment methods, New York Macmillan 1968.
[11] G.J.Burke and A.J.Poggio, Numerical Electromagnetic Code (NEC) - Method of Moments, Part - I : Theory,
Lawrence Livermore National Laboratory, 1981.
[12] 富士通, https://www.fujitsu.com/jp/products/computing/servers/supercomputer/[13] 菅原清文, “C/C++プログラマーのための OpenMP並列プログラミング 第2版,” カットシステム, 2012.
[14] P.パチェコ(秋葉博訳), “MPI並列プログラミング,” 培風館, 2001.

64

[15] 渡辺健介他, “スーパーコンピュータ「富岳」向けのアプリケーション開発環境,” 富士通テクニカルレビュー, 2020年
10月.
[16] H.Nakano, Low-Profile Natural and Metamaterial Antennas, IEEE Press, Wiley, 2016.
[17] V.G.Veselago, “The Electrodynamics of Media with Simultaneously Negative Values of ε and μ,” Soviet Physics
Uspekhi, vol.10, no.4, pp.509-514, 1968.
[18] J.B.Pendry, A.J.Holden, W.J.Stewart,and I.Youngs, “Extremely low frequency plasmons in metallic mesostructures,”Phys. Rev. Lett., vol.76, no.25, pp.4773-4776, 1996.
[19] D.R.Smith, W.J.Padilla, D.C.Vier, S.C.Nemat-Nasser, and S. Schultz, “Composite medium with simultaneously
negative permeability and permittivity,” Phys. Rev. Lett., vol.84, no.18, pp.4184-4187, May 2000.
[20] C.Caloz and T.Itoh, Electromagnetic Metamaterials: Transmission Line Theory and Microwave Applications, IEEE
Press, Wiley, 2006.
[21] 宇野亨, 道下尚文, “メタマテリアルアンテナの基礎,” コロナ社, 2021.
[22] P. B. Johnson and R. W. Christy, "Optical Constants of the Noble Metals," Physical Review B, vol.6, no.12, 1972.
[23] 電気学会編, ”計算電磁気学,” 培風館, 2003.
[24] 高原淳一, “メタサーフェス-新しい平面光学素子の原理と産業化への展望-,” 電子情報通信学会誌, vol.105, no.1,
pp.39-46, 2022年1月.
[25] 大賀明夫, “FDTD法の並列化技術とオープンソース化,” 電子情報通信学会技術研究報告, AP2014-41, pp.7-12, 2014.
[26] 大賀明夫, “OpenFDTDとOpenMOMで始めるはじめてのアンテナ・シミュレーション, “ RFワールド, No.39, pp.782, CQ出版, 2017年8月.
[27] 大賀明夫, “HF/VHF帯線状アンテナの特性シミュレーション,” RFワールド, No.42, pp.61-81, CQ出版, 2018年5月.
[28] 大賀明夫, “機能無制限&GPU対応! 3D電磁界シミュレータOpenFDTD,” トランジスタ技術, 2020年1月号, pp.73-87,
65
CQ出版, 2020.