Python 3D FDTD模拟器

Python 3D FDTD模拟器

翻译自：https://github.com/flaport/fdtd
未获得作者翻译授权，只是为方便自己查看。

一个用Python编写的三维电磁FDTD模拟器。FDTD模拟器有一个可选的PyTorch后端，支持GPU上的FDTD模拟。

安装

fdtd库可以用pip安装：

pip install fdtd

可以通过克隆存储库来安装开发版本：

git clone http://github.com/flaport/fdtd

并与pip连接

pip install -e fdtd

可以安装开发依赖项：

pip install -e fdtd[dev]

依赖项

• python 3.6+
• numpy
• scipy
• matplotlib
• tqdm
• pytorch (optional)

贡献

欢迎所有改进或添加(例如新的对象、源或检测器)。请提出拉拔请求 😊。

文档

在这里阅读文档： https://fdtd.readthedocs.org

导入

fdtd库简单地导入如下：

import fdtd

设置后端

fdtd库允许选择后端。numpy后端是默认的，但也有几个额外的PyTorch后端:

• "numpy" (默认为float64数组)
• "torch" (默认为float64数组)
• "torch.float32"
• "torch.float64"
• "torch.cuda" (默认为float64数组)
• "torch.cuda.float32"
• "torch.cuda.float64"

例如，这是如何选择torch后端：

fdtd.set_backend("torch")

一般来说，numpy后端首选用于具有float64精度的标准CPU计算。一般来说，在FDTD模拟中，float64的精度总是优于float32，然而，float32可能会提供显著的性能提升。

cuda后端只适用于带有GPU的计算机。

FDTD网格

FDTD网格定义了仿真区域。

# signature
fdtd.Grid(
    shape: Tuple[Number, Number, Number],
    grid_spacing: float = 155e-9,
    permittivity: float = 1.0,
    permeability: float = 1.0,
    courant_number: float = None,
)

网格是由它的shape定义的，它只是一个Number类型(整数或浮点数)的3D元组。如果形状以浮动形式给出，则表示网格的宽度、高度和长度(以米为单位)。如果形状以整数形式给出，则用grid_spacing表示网格的宽度、高度和长度。在内部，这些数字将被转换为三个整数：grid.Nx, grid.Ny 和grid.Nz。

可以给出一个grid_spacing。出于稳定性考虑，建议选择比网格中最小波长至少小10倍的网格间距。这意味着对于包含波长为1550nm的光源和折射率为3.1的材料的栅格，推荐的grid_spacing间距是50pm

对于有一下形式的permittivity 和permeability 的浮点数或数组：

• (grid.Nx, grid.Ny, grid.Nz)
• or (grid.Nx, grid.Ny, grid.Nz, 1)
• or (grid.Nx, grid.Ny, grid.Nz, 3)

在最后一种情况下，这种形状暗示了每个长轴(所谓的单轴或双轴材料)的不同介电常数的可能性。在内部，这些变量将(出于性能原因)转换为它们的反向 grid.inverse_permittivity阵列和grid.inverse_permeability阵列的形状(grid.Nx, grid.Ny, grid.Nz, 3)。在制作网格之后，可以改变这些数组。

最后，网格的courant_number决定了仿真的time_step与grid_spacing之间的关系。如果没有给出，则选择为 Courant-Friedrichs-Lewy Condition所允许的最大值：
1 为 1D 模拟, 1/√2 为2D 模拟 ， 1/√3 为 3D模拟 (根据网格的形状计算维数)。出于稳定性考虑，建议不要更改此值。

grid = fdtd.Grid(
    shape = (25e-6, 15e-6, 1), # 25um x 15um x 1 (grid_spacing) --> 2D FDTD
)
print(grid)

Grid(shape=(161,97,1), grid_spacing=1.55e-07, courant_number=0.70)

向网格中添加对象

另一种局部改变栅极 permittivity或permeability的方法是向栅极添加一个Object。

# signature
fdtd.Object(
    permittivity: Tensorlike,
    name: str = None
)

一个对象用修改的更新方程定义网格的一部分，允许引入例如吸收材料或双轴材料，轴之间的混合通过Pockels coefficients或更多。在这种情况下，我们将使一个物体的permittivity不同于它所在的网格。

就像网格一样，Object期望permittivity为float或以下可能形状的数组

• (obj.Nx, obj.Ny, obj.Nz)
• or (obj.Nx, obj.Ny, obj.Nz, 1)
• or (obj.Nx, obj.Ny, obj.Nz, 3)

注意， obj.Nx, obj.Ny and obj.Nz 没有给出给对象构造函数。它们是从它在网格中的位置派生出来的:

grid[11:32, 30:84, 0] = fdtd.Object(permittivity=1.7**2, name="object")

这里发生了几件事。首先，对象在网格中被赋予空间[11:32,30:84,0]。因为它被赋予了这个空间，对象的Nx ， Ny和Nz会被自动设置。此外，通过为对象提供一个名称，该名称将在网格中可用:

print(grid.object)

    Object(name='object')
        @ x=11:32, y=30:84, z=0:1

第二个对象可以添加到网格中:

grid[13e-6:18e-6, 5e-6:8e-6, 0] = fdtd.Object(permittivity=1.5**2)

这里选择了一个带有浮点数的切片。在对象注册期间，这些浮点数将被整数Nx、Ny和Nz替换。由于该对象没有接收到名称，因此该对象不能作为网格的属性使用。然而，它仍然可以通过 grid.objects 列表:

print(grid.objects)

[Object(name='object'), Object(name=None)]

这个列表存储所有对象(例如类型为fdtd.Object的对象)，按照它们被添加到网格的顺序

向网格中添加一个源

类似于将对象添加到网格中，一个fdtd.LineSource也可以添加:

# signature
fdtd.LineSource(
    period: Number = 15, # timesteps or seconds
    power: float = 1.0,
    phase_shift: float = 0.0,
    name: str = None,
)

就像fdtd.Object一样,一个fdtd.LineSource 的大小由它在网格中的位置来定义:

grid[7.5e-6:8.0e-6, 11.8e-6:13.0e-6, 0] = fdtd.LineSource(
    period = 1550e-9 / (3e8), name="source"
)

然而，需要注意的是，在这种情况下，一个LineSource被添加到网格中，也就是说，源跨越了由切片定义的立方体的对角线。在内部，这些片将被转换为列表，以确保以下行为:

print(grid.source)

    LineSource(period=14, power=1.0, phase_shift=0.0, name='source')
        @ x=[48, ... , 51], y=[76, ... , 83], z=[0, ... , 0]

请注意，我们也可以首先提供列表来索引网格。这个特性对于创建任意形状的LineSource非常有用。

给网格添加一个检测器

# signature
fdtd.LineDetector(
    name=None
)

向网格中添加检测器的工作原理与添加源相同

grid[12e-6, :, 0] = fdtd.LineDetector(name="detector")

print(grid.detector)

    LineDetector(name='detector')
        @ x=[77, ... , 77], y=[0, ... , 96], z=[0, ... , 0]

添加网格边界

# signature
fdtd.PML(
    a: float = 1e-8, # stability factor
    name: str = None
)

虽然，有一个对象，源和探测器，以进行FDTD模拟在原则上是足够的，人们还需要定义一个网格边界，以防止场被反射。其中一个可以添加到网格的边界是 Perfectly Matched Layer or PML完美匹配层。这些基本上是吸收边界。

# x boundaries
grid[0:10, :, :] = fdtd.PML(name="pml_xlow")
grid[-10:, :, :] = fdtd.PML(name="pml_xhigh")

# y boundaries
grid[:, 0:10, :] = fdtd.PML(name="pml_ylow")
grid[:, -10:, :] = fdtd.PML(name="pml_yhigh")

网格的总结

可以通过打印出网格来显示网格的简单摘要：

print(grid)

Grid(shape=(161,97,1), grid_spacing=1.55e-07, courant_number=0.70)

sources:
    LineSource(period=14, power=1.0, phase_shift=0.0, name='source')
        @ x=[48, ... , 51], y=[76, ... , 83], z=[0, ... , 0]

detectors:
    LineDetector(name='detector')
        @ x=[77, ... , 77], y=[0, ... , 96], z=[0, ... , 0]

boundaries:
    PML(name='pml_xlow')
        @ x=0:10, y=:, z=:
    PML(name='pml_xhigh')
        @ x=-10:, y=:, z=:
    PML(name='pml_ylow')
        @ x=:, y=0:10, z=:
    PML(name='pml_yhigh')
        @ x=:, y=-10:, z=:

objects:
    Object(name='object')
        @ x=11:32, y=30:84, z=0:1
    Object(name=None)
        @ x=84:116, y=32:52, z=0:1

运行一个仿真

运行一个模拟就像使用grid.run方法一样简单。

grid.run(
    total_time: Number,
    progress_bar: bool = True
)

与网格中的长度一样，模拟的total_time可以指定为整数(time_steps的数量)或浮点数(以秒为单位)。

grid.run(total_time=100)

网格可视化

Le让我们把网格形象化。这可以通过grid.visualize来实现方法:

# signature
grid.visualize(
    grid,
    x=None,
    y=None,
    z=None,
    cmap="Blues",
    pbcolor="C3",
    pmlcolor=(0, 0, 0, 0.1),
    objcolor=(1, 0, 0, 0.1),
    srccolor="C0",
    detcolor="C2",
    show=True,
)

默认情况下，此方法将可视化网格中的所有对象，以及在某个x、y或z平面上当前time_step的功率。通过设置show=False，可以禁用matplotlib映像的即时可视化。

grid.visualize(z=0)

原理

对麦克斯韦方程组的FDTD离散化的尽可能快的解释。

更新方程

电磁时域有限差分(FDTD)求解器可以求解随时间变化的麦克斯韦方程组

    curl(H) = ε*ε0*dE/dt
    curl(E) = -µ*µ0*dH/dt

这两个方程分别称为安培定律 和法拉第定律。

其中ε和µ分别为相对介电常数张量和渗透率张量。ε0和µ0分别为真空介电常数和磁导率，其平方根分别被E和H吸收；表示为 E := √ε0*E 和 H := õ0*H。

这样，麦克斯韦方程可以写成更新方程:

    E  += c*dt*inv(ε)*curl(H)
    H  -= c*dt*inv(µ)*curl(E)

然后电场和磁场可以在一个交错的yee坐标网格上离散，在3D中看起来像这样:

根据Yee离散算法，网格上固有两种类型的场:整数网格位置上的E型场和半整数网格位置上的H型场。

这些交错坐标的美妙之处在于，它们可以用一种非常自然的方式来表示电场和磁场的旋度：H型场的旋度将是E型场，反之亦然。

这样，E的旋度就可以写成

curl(E)[m,n,p]=(dEz/dy - dEy/dz, dEx/dz - dEz/dx, dEy/dx - dEx/dy)[m,n,p]
              =( ((Ez[m,n+1,p]-Ez[m,n,p])/dy - (Ey[m,n,p+1]-Ey[m,n,p])/dz),
               ((Ex[m,n,p+1]-Ex[m,n,p])/dz - (Ez[m+1,n,p]-Ez[m,n,p])/dx),
               ((Ey[m+1,n,p]-Ey[m,n,p])/dx - (Ex[m,n+1,p]-Ex[m,n,p])/dy) )
              =(1/du)*( ((Ez[m,n+1,p]-Ez[m,n,p]) - (Ey[m,n,p+1]-Ey[m,n,p])), [assume dx=dy=dz=du]
               ((Ex[m,n,p+1]-Ex[m,n,p]) - (Ez[m+1,n,p]-Ez[m,n,p])),
               ((Ey[m+1,n,p]-Ey[m,n,p]) - (Ex[m,n+1,p]-Ex[m,n,p])) )

这可以用数组切片来有效地表示(注意，因子(1/du)被省略了)

def curl_E(E):
    curl_E = np.zeros(E.shape)
    curl_E[:,:-1,:,0] += E[:,1:,:,2] - E[:,:-1,:,2]
    curl_E[:,:,:-1,0] -= E[:,:,1:,1] - E[:,:,:-1,1]

    curl_E[:,:,:-1,1] += E[:,:,1:,0] - E[:,:,:-1,0]
    curl_E[:-1,:,:,1] -= E[1:,:,:,2] - E[:-1,:,:,2]

    curl_E[:-1,:,:,2] += E[1:,:,:,1] - E[:-1,:,:,1]
    curl_E[:,:-1,:,2] -= E[:,1:,:,0] - E[:,:-1,:,0]
    return curl_E

H的旋度可以用类似的方法得到(再次注意，省略了因子(1/du)):

def curl_H(H):
    curl_H = np.zeros(H.shape)

    curl_H[:,1:,:,0] += H[:,1:,:,2] - H[:,:-1,:,2]
    curl_H[:,:,1:,0] -= H[:,:,1:,1] - H[:,:,:-1,1]

    curl_H[:,:,1:,1] += H[:,:,1:,0] - H[:,:,:-1,0]
    curl_H[1:,:,:,1] -= H[1:,:,:,2] - H[:-1,:,:,2]

    curl_H[1:,:,:,2] += H[1:,:,:,1] - H[:-1,:,:,1]
    curl_H[:,1:,:,2] -= H[:,1:,:,0] - H[:,:-1,:,0]
    return curl_H

更新方程现在可以重写为

    E  += (c*dt/du)*inv(ε)*curl_H
    H  -= (c*dt/du)*inv(µ)*curl_E

这个数(c*dt/du)是一个无量纲参数，称为Courant数sc。出于稳定性的原因，Courant数应该总是小于1/√D, D是模拟的维数。这可以直观地理解为信息在网格中的传播速度总是低于光速的条件。在这里描述的FDTD方法中，信息只能传播到邻近的网格单元(通过旋度的应用)。因此，在一个D维立方体(2D的正方形，3D的立方体)的对角线上移动需要D个时间步长，那么Courant条件就会自动遵循这个事实，即对角线的长度是1/√D。

这将产生FDTD算法的最终更新方程：

    E  += sc*inv(ε)*curl_H
    H  -= sc*inv(µ)*curl_E

这是它的实现方式：

class Grid:
    # ... [initialization]

    def step(self):
        self.update_E()
        self.update_H()

    def update_E(self):
        self.E += self.courant_number * self.inverse_permittivity * curl_H(self.H)

    def update_H(self):
        self.H -= self.courant_number * self.inverse_permeability * curl_E(self.E)

源

安培定律可以更新为包含电流密度:

    curl(H) = J + ε*ε0*dE/dt

再次进行常规替换sc:= c*dt/du, E:=√ε0*E, H:=õ0*H，可将更新方程修改为包含电流密度：

    E += sc*inv(ε)*curl_H - dt*inv(ε)*J/√ε0

最后做一个替换Es:= -dt*inv(ε)*J/√ε0使我们可以用一种非常简洁的方式来表示:

    E += sc*inv(ε)*curl_H + Es

我们定义Es为电场源项。

通常定义一个磁场源项 Hs也是有用的，如果它存在，它将从磁电流密度导出。同理，法拉第的更新方程可以改写为

    H  -= sc*inv(µ)*curl_E + Hs

class Source:
    # ... [initialization]
    def update_E(self):
        # electric source function here

    def update_H(self):
        # magnetic source function here

class Grid:
    # ... [initialization]
    def update_E(self):
        # ... [electric field update equation]
        for source in self.sources:
            source.update_E()

    def update_H(self):
        # ... [magnetic field update equation]
        for source in self.sources:
            source.update_H()

损耗介质

当材料具有导电性 σe时，导电电流将保证介质是有损耗的。安培定律与传导电流的关系

    curl(H) = σe*E + ε*ε0*dE/dt

用通常的替换法，这就变成了:

    E(t+dt) - E(t) = sc*inv(ε)*curl_H(t+dt/2) - dt*inv(ε)*σe*E(t+dt/2)/ε0

这个更新方程取决于半整数时间步长(磁场时间步长)上的电场。我们需要做一个替换，将电场插值到这个时间步长:

    (1 + 0.5*dt*inv(ε)*σ/√ε0)*E(t+dt) = sc*inv(ε)*curl_H(t+dt/2) + (1 - 0.5*dt*inv(ε)*σe/ε0)*E(t)

代入σ:= inv(ε)*σe/ε0，得到新的更新方程:

    f = 0.5*dt*σ
    E *= inv(1 + f) * (1 - f)
    E += inv(1 + f)*sc*inv(ε)*curl_H

注意，介电常数张量ε越复杂，这个算法就会消耗越多的时间。因此，有时通过引入(非物理)磁导率将吸收转移到磁畴是正确的决定，因为磁导张量µ通常等于1。

将σ:= inv(µ)*σm/µ0代入，得到磁场更新方程:

    f = 0.5*dt*σ
    H *= inv(1 + f) * (1 - f)
    H += inv(1 + f)*sc*inv(µ)*curl_E

能量密度和坡印亭矢量

电磁能量密度可由此式给出：

    e = (1/2)*ε*ε0*E**2 + (1/2)*µ*µ0*H**2

通过上面的替换，这就变成了模拟单元:

    e = (1/2)*ε*E**2 + (1/2)*µ*H**2

坡印亭矢量由此式给出：

    P = E×H

模拟单元变成

    P = c*E×H

源Es引入的能量可以通过跟踪能量密度的变化得到

    de = ε*Es·E + (1/2)*ε*Es**2

这也可以从Poyntings能量守恒定律中得到:

    de/dt = -grad(S) - J·E

第一项描述了体积内能量的再分配第二项描述了电流密度引入的能量。

注意:虽然它是非物理的，但也可以引入一个磁源。这一来源将引入下列能量:

    de = ε*Hs·H + (1/2)*µ*Hs**2

由于µ张量通常等于1，使用磁源项通常更有效。

类似地，我们也可以用下列方法来跟踪由于电导率而吸收的能量:

    f = 0.5*dt*σ
    Enoabs = E + sc*inv(ε)*curl_H
    E *= inv(1 + f) * (1 - f)
    E += inv(1 + f)*sc*inv(ε)*curl_H
    dE = Enoabs - E
    e_abs += ε*E*dE + 0.5*ε*dE**2

或者如果我们想通过磁，磁导电性来记录吸收的能量：

    f = 0.5*dt*inv(µ)*σ
    Hnoabs = E + sc*inv(µ)*curl_E
    H *= inv(1 + f) * (1 - f)
    H += inv(1 + f)*sc*inv(µ)*curl_E
    dH = Hnoabs - H
    e_abs += µ*H*dH + 0.5*µ*dH**2

能量密度中的电项和磁项通常具有相同的大小。因此，引入导磁率σm所吸收的能量与引入导磁率σe所吸收的能量相同，如果:

    inv(µ)*σm/µ0 = inv(ε)*σe/ε0

边界条件

周期性边界条件

假设我们想要沿x方向的周期性边界条件，那么我们必须确保Xlow和Xhigh的场是相同的。这必须在执行更新方程式后执行:

注意电场E依赖于curl_H，这意味着E的第一个指标不会通过更新方程更新。这些指标需要通过周期边界条件来设定。具体来说，需要将E[0]设为E[-1]。对于磁场，情况正好相反:H依赖于curl_E，这意味着它的最后指标不会被设置。这需要通过边界条件来完成:H[-1]需要被设为H[0]:

class PeriodicBoundaryX:
    # ... [initialization]
    def update_E(self):
        self.grid.E[0, :, :, :] = self.grid.E[-1, :, :, :]

    def update_H(self):
        self.grid.H[-1, :, :, :] = self.grid.H[0, :, :, :]

class Grid:
    # ... [initialization]
    def update_E(self):
        # ... [electric field update equation]
        # ... [electric field source update equations]
        for boundary in self.boundaries:
            boundary.update_E()

    def update_H(self):
        # ... [magnetic field update equation]
        # ... [magnetic field source update equations]
        for boundary in self.boundaries:
            boundary.update_H()

完全匹配层

完美匹配层(PML)是在FDTD网格中引入吸收边界条件的最新技术。PML是电网中的阻抗匹配吸收区。结果表明，为了保持阻抗匹配条件，PML只能在一个方向上吸收。这就是为什么PML实际上是一种非物理材料。

考虑Ez分量的安培定律，其中常用的替换E:=√ε0*E, H:=õ0*H和σ:= inv(ε)*σ E /ε0已经引入:

    ε*dEz/dt + ε*σ*Ez = c*dHy/dx - c*dHx/dy

这就变成了频域:

    iω*ε*Ez + ε*σ*Ez = c*dHy/dx - c*dHx/dy

们可以把这个方程分解成x波和y波:

    iω*ε*Ezx + ε*σx*Ezx = iω*ε*(1 + σx/iω)*Ezx = c*dHy/dx
    iω*ε*Ezy + ε*σy*Ezy = iω*ε*(1 + σy/iω)*Ezy = -c*dHx/dy

我们可以这样定义s算子

    Su = 1 + σu/iω          with u in {x, y, z}

一般来说，我们倾向于在Su中加入稳定因子au和比例因子ku:

    Su = ku + σu/(iω+au)    with u in {x, y, z}

将Ez的两个方程相加，然后除以各自的S算符得到

    iω*ε*Ez = (c/Sx)*dHy/dx - (c/Sy)*dHx/dy

把它转换回时域

    ε*dEz/dt = c*sx[*]dHy/dx - c*sx[*]dHx/dy

其中sx为(1/ sx)的傅里叶反变换，[*]为卷积。
su的表达式可以被证明[经过一些推导]如下所示:

    su = (1/ku)*δ(t) + Cu(t)    with u in {x, y, z}

其中δ(t)为狄拉克函数，C(t)为指数衰减函数：

    Cu(t) = -(σu/ku**2)*exp(-(au+σu/ku)*t)     for all t > 0 and u in {x, y, z}

代入得到：

    dEz/dt = (c/kx)*inv(ε)*dHy/dx - (c/ky)*inv(ε)*dHx/dy + c*inv(ε)*Cx[*]dHy/dx - c*inv(ε)*Cx[*]dHx/dy
           = (c/kx)*inv(ε)*dHy/dx - (c/ky)*inv(ε)*dHx/dy + c*inv(ε)*Фez/du      with du=dx=dy=dz

这可以写成一个更新等式：

    Ez += (1/kx)*sc*inv(ε)*dHy - (1/ky)*sc*inv(ε)*dHx + sc*inv(ε)*Фez

我们把Фeu定义为：

    Фeu = Ψeuv - Ψezw           with u, v, w in {x, y, z}

取Hw在v方向上的导数来更新分量Eu的卷积Ψeuv：

    Ψeuv = dv*Cv[*]dHw/dv     with u, v, w in {x, y, z}

这可以重写[在一些推导之后]作为更新方程本身：

     Ψeuv = bv*Ψeuv + cv*dv*(dHw/dv)
          = bv*Ψeuv + cv*dHw            with u, v, w in {x, y, z}

其中常量bu和cu被推导为：

    bu = exp(-(au + σu/ku)*dt)              with u in {x, y, z}
    cu = σu*(bu - 1)/(σu*ku + au*ku**2)     with u in {x, y, z}

电场的最终PML算法现在变成：

1. 使用Ψ部分的更新方程更新Фe=[Фex， Фey， Фez]。
2. 用正常的方法更新电场。
3. 将Фe加到电场中。

或者作为python代码：

class PML(Boundary):
    # ... [initialization]
    def update_phi_E(self): # update convolution
        self.psi_Ex *= self.bE
        self.psi_Ey *= self.bE
        self.psi_Ez *= self.bE

        c = self.cE
        Hx = self.grid.H[self.locx]
        Hy = self.grid.H[self.locy]
        Hz = self.grid.H[self.locz]

        self.psi_Ex[:, 1:, :, 1] += (Hz[:, 1:, :] - Hz[:, :-1, :]) * c[:, 1:, :, 1]
        self.psi_Ex[:, :, 1:, 2] += (Hy[:, :, 1:] - Hy[:, :, :-1]) * c[:, :, 1:, 2]

        self.psi_Ey[:, :, 1:, 2] += (Hx[:, :, 1:] - Hx[:, :, :-1]) * c[:, :, 1:, 2]
        self.psi_Ey[1:, :, :, 0] += (Hz[1:, :, :] - Hz[:-1, :, :]) * c[1:, :, :, 0]

        self.psi_Ez[1:, :, :, 0] += (Hy[1:, :, :] - Hy[:-1, :, :]) * c[1:, :, :, 0]
        self.psi_Ez[:, 1:, :, 1] += (Hx[:, 1:, :] - Hx[:, :-1, :]) * c[:, 1:, :, 1]

        self.phi_E[..., 0] = self.psi_Ex[..., 1] - self.psi_Ex[..., 2]
        self.phi_E[..., 1] = self.psi_Ey[..., 2] - self.psi_Ey[..., 0]
        self.phi_E[..., 2] = self.psi_Ez[..., 0] - self.psi_Ez[..., 1]

    def update_E(self): # update PML located at self.loc
        self.grid.E[self.loc] += (
            self.grid.courant_number
            * self.grid.inverse_permittivity[self.loc]
            * self.phi_E
        )

class Grid:
    # ... [initialization]
    def update_E(self):
        for boundary in self.boundaries:
            boundary.update_phi_E()
        # ... [electric field update equation]
        # ... [electric field source update equations]
        for boundary in self.boundaries:
            boundary.update_E()

磁场也是一样的。

These update equations for the PML were based on
Schneider, Chap. 11.

License

© Floris laporte - MIT License