【Python】FDTD 法で z 方向を進む 1 次元の電磁波解析 #Numpy
真空中の \(z\) 方向へ進む電磁波を FDTD (Finite Difference Time Domain) 法で解析する。
前回，FDTD 法で解析するコードを Python の標準ライブラリのみで記載した。今回は，Numpy を使って，時間刻み以外では for を使わずに，FDTD 法で解析を行った。
解析コード
\(z\) 方向は 0 ～ 2 000 mm の解析領域を設ける。\(z\) = 1 000 mm の位置に，広がり 100 mm のガウシアンパルスの電界を波源として与える。
このパルスが真空中を伝搬していく様を FDTD 法でシミュレーションする。
ライブラリの読み込み
解析に使用する Python のライブラリをインポートする。今回は，標準ライブラリ math に加えて，Numpy も使用する。
import math
import numpy as npパラメータの設定
解析に用いるパラメータを設定する。
ep=8.85418782e-12#真空の誘電率
mu=4*math.pi*1e-7#真空の透磁率
cc=1/math.sqrt(ep*mu)#光速
dd=2e-3#メッシュサイズ 2 mm
NZ=1000#z 方向の格子数
CFL=0.5#CFL
dt=CFL*dd/cc#時間刻み
NT=1500#時間ステップ
ee=dt/ep/dd
hh=dt/mu/dd
cz=(cc*dt-dd)/(cc*dt+dd)
z0=NZ/2*dd#ガウシアンパルスの中心
sig=100e-3#ガウシアンパルスの広がり初期条件の設定 ～ z 軸，電界，磁界の配列～
\(z\) 軸，\(z\) 方向へ進む電磁界として \(E_x\)，\(E_y\)，\(H_x\)，\(H_y\) の配列を用意する。ここで，電界 \(E_x\) にガウシアンパルスを与えている。
ここから，Numpy を活用することで，シンプルな記述で配列が準備できる。
zz=np.arange(0,NZ*dd,dd)
Ex=np.exp(-(zz - z0)**2 /sig/sig)
Ey=np.zeros(NZ)
Hx=np.zeros(NZ-1)
Hy=np.zeros(NZ-1)FDTD 法
FDTD 法では，電界の更新，磁界の更新を順に行う。また，解析領域の境界には，Mur の 1 次吸収境界条件を適用した。
for t in range(0,NT+1):
# Update Ex,Ey
oldEx0=Ex[1]
oldEy0=Ey[1]
oldExN=Ex[NZ-2]
oldEyN=Ey[NZ-2]
Ex[1:NZ-1]=Ex[1:NZ-1]-ee*(Hy[1:NZ-1]-Hy[0:NZ-2])
Ey[1:NZ-1]=Ey[1:NZ-1]+ee*(Hx[1:NZ-1]-Hx[0:NZ-2])
# Mur's first order absorbing boundary condition
Ex[0]=oldEx0+cz*(Ex[1]-Ex[0])
Ey[0]=oldEy0+cz*(Ey[1]-Ey[0])
Ex[NZ-1]=oldExN+cz*(Ex[NZ-2]-Ex[NZ-1])
Ey[NZ-1]=oldEyN+cz*(Ey[NZ-2]-Ey[NZ-1])
# Update Hx,Hy
oldHx0=Hx[1]
oldHy0=Hy[1]
oldHxN=Hx[NZ-3]
oldHyN=Hy[NZ-3]
Hx[1:NZ-2]=Hx[1:NZ-2]+hh*(Ey[2:NZ-1]-Ey[1:NZ-2])
Hy[1:NZ-2]=Hy[1:NZ-2]-hh*(Ex[2:NZ-1]-Ex[1:NZ-2])
# Mur's first order absorbing boundary condition
Hx[0]=oldHx0+cz*(Hx[1]-Hx[0])
Hy[0]=oldHy0+cz*(Hy[1]-Hy[0])
Hx[NZ-2]=oldHxN+cz*(Hx[NZ-3]-Hx[NZ-2])
Hy[NZ-2]=oldHyN+cz*(Hy[NZ-3]-Hy[NZ-2])解析結果
真空中の \(z\) 方向へ進む電磁波を FDTD 法で解析した結果を示す。電界 \(E_x\)，磁界 \(E_y\) の伝搬をアニメーションで示している。
解析時間の比較
Python の time モジュールを用いて，解析時間の計測を行う。そして，Numpy を使った場合，Numpy を使わない場合の解析時間を比較する。
import time
start=time.time()
# パラメータの設定
# z 軸，電界，磁界の配列
# FDTD 法
end=time.time()
time_diff=end-start
print(time_diff)解析時間の実績
解析時間の実績を次表に示す。Numpy を使用した場合の平均解析時間 0.048 秒に対し，使用しない場合は 1.363 秒である。Numpy 使用有無で約 28.7 倍の差がある。
FDTD 法のような配列計算を多く実施する場合，Numpy を使うのが有効だ。
| 回数 | Numpy 使用 | Numpy 不使用 |
|---|---|---|
| 1 回目 | 0.0473024845123291 | 1.40471911430359 |
| 2 回目 | 0.0498981475830078 | 1.38675308227539 |
| 3 回目 | 0.0468921661376953 | 1.37017369270325 |
| 4 回目 | 0.0468785762786865 | 1.3435435295105 |
| 5 回目 | 0.0468978881835938 | 1.31075477600098 |
| 平均 | 0.0475738525390625 | 1.36318883895874 |
（参考）解析に使用したPC，Python
今回，解析に使用した PC のスペック，Python のバージョンを以下に示す。
- プロセッサ : AMD Ryzen 7 2700 Eight-Core Processor 3.20 GHz
- 実装 RAM : 16.0 GB
- OS : Windows 11 Home 22H2
- Python 3.11.1
- numpy 1.24.1
まとめ
FDTD 法で \(z\) 方向を進む 1 次元の電磁波の伝搬を解析した。今回，Numpy を使うことで，少ない行数で FDTD 法を実装できた。また，解析時間は約 1/30 にすることができた。
