关于FDTD
FDTD是计算电磁学领域用的最多的算法之一（其他的包括MoM，以及大名鼎鼎的FEM），FDTD也可用于求解一些量子力学中的问题，其核心思想就是对电磁场中的电场、磁场分别在时间域和空间域进行差分，然后通过时间和空间上的迭代进行时域求解，这里以一维的FDTD为例，写出一维情况下的Maxwell 方程组：
\begin{aligned} - \frac { \partial H _ { y } } { \partial z } & = \varepsilon \frac { \partial E _ { x } } { \partial t } + \sigma E _ { x } \\ \frac { \partial E _ { x } } { \partial z } & = - \mu \frac { \partial H _ { y } } { \partial t } - \sigma _ { m } H _ { y } \end{aligned}
对其分别在空间和时间上进行差分
E _ { x } ^ { n- 1 } ( k ) = \mathrm { CA } ( m ) \cdot E _ { x } ^ { n } ( k ) - \mathrm { CB } ( m ) \cdot \frac { H _ { y } ^ { n + 1 / 2 } \left( k + \frac { 1 } { 2 } \right) - H _ { y } ^ { n + 1 / 2 } \left( k - \frac { 1 } { 2 } \right) } { \Delta z }
H _ { y } ^ { n+ 1 / 2 } \left( k + \frac { 1 } { 2 } \right) = \mathrm { CP } ( m ) \cdot H _ { y } ^ { n- 1 / 2 } \left( k + \frac { 1 } { 2 } \right) - \mathrm { CQ } ( m ) \cdot \frac { E _ { x }^n ( k + 1 ) - E _ { x } ^n( k ) } { \Delta z }
其中
\begin{aligned} C A ( m ) & = \frac { \frac { \varepsilon ( m ) } { \Delta t } - \frac { \sigma ( m ) } { 2 } } { \frac { \varepsilon ( m ) } { \Delta t } + \frac { \sigma ( m ) } { 2 } } = \frac { 1 - \frac { \sigma ( m ) \Delta t } { 2 \epsilon ( m ) } } { 1 + \frac { \sigma ( m ) \Delta t } { 2 \varepsilon ( m ) } } \\ C B ( m ) & = \frac { 1 } { \frac { \varepsilon ( m ) } { \Delta t } + \frac { \sigma ( m ) } { 2 } } = \frac { \frac { \Delta t } { \varepsilon ( m ) } } { 1 + \frac { \sigma ( m ) \Delta t } { 2 \varepsilon ( m ) } } \\ C P ( m ) & = \frac { \frac { \mu ( m ) } { \Delta t } - \frac { \sigma _ { m } ( m ) } { 2 } } { \frac { \mu ( m ) } { \Delta t } + \frac { \sigma _ { m } ( m ) } { 2 } } = \frac { \frac { 1 - \sigma _ { m } ( m ) \Delta t } { 2 \mu ( m ) } } { 1 + \frac { \sigma _ { m } ( m ) \Delta t } { 2 \mu ( m ) } } \end{aligned}
\mathrm { CQ } ( m ) = \frac { 1 } { \frac { \mu ( m ) } { \Delta t } + \frac { \sigma _ { m } ( m ) } { 2 } } = \frac { \frac { \Delta t } { \mu ( m ) } } { 1 + \frac { \sigma _ { m } ( m ) \Delta t } { 2 \mu ( m ) } }
其中 k 和 n 分别为空间步和时间步的索引，详细过程可以参考《电磁波时域有限差分方法》，葛德彪著。
利用Julia编写FDTD
之前有写过CPP版本的（其实是祖传代码），现在想把代码移植到julia，至少画图方便。话不多说，直接把我写好的代码放出来：
using Plots
using Revise
using SharedArrays
using LaTeXStrings
using Distributed
gr()
const sigma = 0.05
const ddx = 1e-10
const c = 3e8
const dt = ddx / (2 * c)
const epsz = 8.85419e-12
# eaf = dt * sigma / (2 * epsz * epsilon)
const eaf = 0 # 不考虑损耗
# 下面三个函数都是源
function Gaussian(n)
spread = 12.0
t0 = 40
return exp(-0.5 * ((t0 - n) / spread) ^ 2)
end
function Modulate_Gaussian(n, freq_in)
spread = 100000
t0 = 0.8 * spread
return -cos(freq_in * dt * n) * exp(-4 * pi * ((t0 - n) / spread)^2)
end
function Sin(n, freq_in)
freq_imag = freq_in * 0.1
#return sin(freq_in * dt *n) * exp(-freq_imag * dt * n)
return sin(freq_in * dt *n)
end
# 迭代电场
function calcEx(ex, hy, pluse, KE, D, ex_low_m1, ex_low_m2, ex_high_m1, ex_high_m2)
i = 1 # i用于计算AB结构一维光子晶体
for k = 2:KE
if k % D == 0
i = i + 1
end
if k == 800
ex[k] = pluse
elseif i % 2 == 0 && k < i * D && k > (i - 1) * D
ex[k] += 0.5 * (hy[k - 1] - hy[k])
else
ex[k] = ca * ex[k] + cb * (hy[k - 1] - hy[k])
end
end
#吸收边界条件
ex[1] = ex_low_m2;
ex_low_m2 = ex_low_m1;
ex_low_m1 = ex[2];
ex[KE] = ex_high_m2;
ex_high_m2 = ex_high_m1;
ex_high_m1 = ex[KE - 1];
return ex, hy, ex_low_m1, ex_low_m2, ex_high_m1, ex_high_m2
end
# 迭代磁场
function calcHy(ex, hy, KE)
for k = 1:KE-1
hy[k] += 0.5 * (ex[k] - ex[k + 1])
end
return ex, hy
end
function Time()
ex = SharedArray{Float32}(KE)
hy = SharedArray{Float32}(KE)
ex_low_m1 = 0
ex_low_m2 = 0
ex_high_m1 = 0
ex_high_m2 = 0
for n = 1:TimeSteps
pluse = Sin(n, freq_in)
ex, hy, ex_low_m1, ex_low_m2, ex_high_m1, ex_high_m2 = calcEx(ex, hy, pluse, KE, D, ex_low_m1, ex_low_m2, ex_high_m1, ex_high_m2)
ex, hy = calcHy(ex, hy, KE)
end
return ex, hy
end
epsilon = 4
ca = (1.0 - eaf) / (1.0 + eaf)
cb = 0.5 / (epsilon * (1 + eaf))
D = 1000
TimeSteps = 500000
freq_in = pi * 1e15
layer = 21 # 光子晶体层数
KE = layer * D
r_i = 100
t_i = KE - 100
@time ex, hy = Time()
x = 0.1:0.1:KE/10
plot(x, ex, xlabel = L"x(nm)", ylabel = L"E")
由于本人主要方向是光子晶体，所以这里做了一下小改动以方便计算AB结构光子晶体。
并行化
FDTD本身迭代是有多个循环的，而且既要迭代电场也要迭代磁场，两个迭代函数还一直需要互相传递数据，所以是很耗费计算资源的，能够并行化就比较好了（之前本人做过FDTD的CUDA并行，有兴趣的可以一起讨论），但是我现在在julia下试过threads和distributed效果都不太好，速度会比原来的还慢（最主要还是我菜），threads的时候这个i必须global才能计算，但是我的需求是每次循环i都初始化为1，这时候global就不太合适了。之前用传输矩阵法计算光子晶体我都是用Julia写的，而且并行效率很高，FDTD就显得麻烦一点。
其他
编程方面我实在是小菜鸡，发这个贴子只是想抛砖引玉下（没想到会是第一个计算物理的帖子）。在Python下有很多计算物理的库，我了解到的光学下的就有基于傅里叶光学的lightpipes，基于FDTD的Meep，以及用于计算光子晶体能带的MPB，Julia下光学电磁学相关的还没有了解到有，虽然这些包都可以通过pycall拿来直接用，但是也希望Julia的生态能越做越大。大家有研究光学、电磁学或者光子晶体的欢迎一起讨论鸭。