FDTD 101: Introduction to perfectly matched layer (PML)
Lecture 6: Introduction to perfectly matched layer (PML)
Perfectly matched layer (PML) is commonly used to truncate unbounded computational region, since an ideal PML can completely absorb the incoming waves from all angles of the incidence without any reflection. In this lecture, we explain the basic idea behind PML, and show how to characterize the performance of PML.
- Explain the concept of PML through coordinate transformation.
- Show that the performance of PML should be characterized by the amplitude rather than the intensity of parasitic reflection, and present a simple approach to compute parasitic reflection at many frequencies.
- Illustrate typical signature in the spectrum when the parasitic reflection at the boundary is significant.
For Python setup used in this lecture, please see related FDTD Python Tutorial below.
Tutorial 6: Truncate computational domain with PML¶
In this notebook, we show how to measure the parasitic reflection from the perfectly matched layer (PML) boundary in order to benchmark the performance of PML. As a simple illustration, we consider a 1D system where a time monitor is placed between a broadband planewave source and PML bounday. By visualizing the evolution of the field at the time monitor, one can clearly distinguish between the incident pulse passing through the time monitor and the reflected pulse. A subsequent Fourier analysis can tell the reflection amplitude at the PML boundary over a broad range of frequencies.
Finally, taking Fabry-Parot cavity as an example, we illustrate the impact of reflection from PML.
# standard python imports
import numpy as np
import matplotlib.pyplot as plt
import tidy3d as td
from tidy3d import web
Next, let us define some key parameters.
# The setup of the 1D simulation:8
# | | | | | | | |
# |`PML`| dPML |`source`| distance |`monitor`| distance |`PML`|
# | | | | | | | |
# Number of PML layers
Npml = 2 # the default number of layers is 12
# sufficiently large distance between the source and the time monitor
distance = 100 #um
# distance between the source and PML
dPML = 1 #um
# Wavelength and frequency range of the source
freq0 = 300e12 #Hz, central frequency
fwidth = freq0/5 #frequency pulse width
lambda0 = td.C_0 / freq0
# resolution
resolution = 20
dl = lambda0/resolution
# runtime
t_stop = 5e-12
# simulation domain size
sim_size = Lx, Ly, Lz = (dl, dl, 2*distance+dPML)
Creat Simulation¶
# Uniform grid
grid_spec = td.GridSpec.uniform(dl=dl)
# Planewave source
src_z = dPML - Lz/2
source = td.PlaneWave(
center=(0,0,src_z),
size=(td.inf, td.inf, 0),
source_time = td.GaussianPulse(
freq0=freq0,
fwidth=fwidth,
),
direction='+',
name='Gaussian pulse',
)
# Time monitor at a single point
tmnt_pt = td.FieldTimeMonitor(
center=(0, 0, src_z + distance),
name='time',
size=(0, 0, 0),
)
# PML boundary
pml = td.Boundary.pml(num_layers=Npml)
bspec = td.BoundarySpec(x=td.Boundary.periodic(), y = td.Boundary.periodic(), z=pml)
Now it is time to define the simulation object.
sim = td.Simulation(
center = (0, 0, 0),
size = sim_size,
grid_spec = grid_spec,
structures = [],
sources = [source],
monitors = [tmnt_pt],
run_time = t_stop,
boundary_spec = bspec,
shutoff = 0,
)
Run the simulations¶
sim_data = web.run(sim, task_name='lecture06_pml_reflection', path=f'data/data_pml_reflection.hdf5')
Post Run Analysis¶
We visualize the time evolution of the field at the monitor plane. Two pulses are clearly visible in the figure. The amplitude of the pulse appearing first is much larger than the one appearing later. To identify the pulses, we compute analytically the time it takes for the center of the pulse to reach the monitor (red dashed line) and for it to get reflected and reach the monitor again (blue dashed line).
plt.rcParams.update({'font.size': 16})
fig, ax = plt.subplots(1,figsize=(7,4), tight_layout=True)
speed = td.C_0
time_reach_monitor = source.source_time.offset/fwidth/2/np.pi + distance/speed
time_reflection = time_reach_monitor + 2*distance/speed
ax.vlines(time_reach_monitor*1e12,-600,600,colors='r',linestyle='dashed')
ax.vlines(time_reflection*1e12,-600,600,colors='b',linestyle='dashed')
time_data = sim_data['time']
Ex = time_data.Ex
Ex_values = np.squeeze(Ex.values)
ax.plot(Ex.t*1e12,Ex_values)
ax.set_xlabel('time (ps)')
ax.set_ylabel('$E_x$')
ax.set_xlim(0,1.5)
ax.set_ylim(-600,600)
plt.show()
Now to compute the reflection amplitude across the frequency range spanned by the Gaussian source, let's perform Fourier analysis.
# First, let's select two time intervals of equal length. The first interval
# completely contains the original incident pulse, and the second interval
# completely contains the reflected pulse. We take the former to
# be [0,time_reach_PML), and the latter to be [time_reach_PML,2*time_reach_PML)
time_reach_PML = source.source_time.offset/fwidth/2/np.pi + 2* distance/speed
inds = np.argwhere(np.array(Ex.t)>=time_reach_PML)
Nt = int(inds[0])
# Now let's perform FFT in each interval, and visualize the spectrum
dt = sim.dt
df = 1/Nt/dt
fmesh = np.arange(0, Nt)*df
# original pulse
spectrum_original = 1/np.sqrt(2*np.pi)*np.fft.ifft(np.fft.ifftshift(Ex_values[:Nt]))/df
#reflected pulse
spectrum_reflection = 1/np.sqrt(2*np.pi)*np.fft.ifft(np.fft.ifftshift(Ex_values[Nt:Nt*2]))/df
# visualization of each spectrum
fig, ax = plt.subplots(1,2,figsize=(8,4), tight_layout=True)
ax[0].plot(fmesh/1e12, np.abs(spectrum_original),label='Original')
ax[0].plot(fmesh/1e12, np.abs(spectrum_reflection),label='Reflected')
ax[0].set_xlim((freq0-fwidth)/1e12,(freq0+fwidth)/1e12)
ax[0].set_ylim(0,8e-13)
ax[0].set_xlabel('Frequency (THz)')
ax[0].set_title('Spectrum')
ax[0].legend()
# visualization of reflection spectrum
ax[1].plot(fmesh/1e12, np.abs(spectrum_reflection/spectrum_original))
ax[1].set_xlim((freq0-fwidth)/1e12,(freq0+fwidth)/1e12)
ax[1].set_yscale('log')
ax[1].set_ylim(1e-6,1)
ax[1].set_title('Reflection amplitude')
ax[1].set_xlabel('Frequency (THz)')
plt.show()
Impact of parasitic PML reflection to Fabry-Perot cavity simulation¶
Finally, let's examine the impact of parasitic reflection from PML boundary to a simple Fabry-Perot cavity transmission simulation illustrated in our previous Lecture 2. Most setups are similar, so please checkout the link here for more detailed explanations.
# Wavelength and frequency range
freq_range = (200e12, 400e12)
lambda_range = (td.constants.C_0/freq_range[1], td.constants.C_0/freq_range[0])
freq0 = np.sum(freq_range)/2
# frequencies and wavelengths of monitor
Nfreq = 301
monitor_freqs = np.linspace(freq_range[0], freq_range[1], Nfreq)
monitor_lambdas = td.constants.C_0 / monitor_freqs
# central frequency, frequency pulse width and total running time
lambda0 = td.C_0 / freq0
bandwidth = 0.4 # bandwidth of source in units of delta frequency. 0.38 for broadband
freqw = bandwidth * (freq_range[1] - freq_range[0])
t_stop = 500 / freq0
# Thickness and refractive index of slab
t_slab = 0.5
n_slab = 3.5
mat_slab = td.Medium(permittivity=n_slab**2, name='silicon')
# Grid resolution (cells per um)
dl = 6e-3
# space between slabs and sources and PML
spacing = 16 * lambda_range[-1]
# simulation size
sim_size = Lx, Ly, Lz = (dl, dl, 2*spacing + t_slab)
Analytical Method¶
def slab_analytical(d, n, wvl):
""" computes transmision as a function of slab thickness (d), refractive index (n), and wavelength (wvl). """
rho = (n-1)/(n+1)
t = ((1+rho)*(1-rho)*np.exp(-2j*np.pi*n*d/wvl)) / (1 - rho**2*np.exp(-4j*np.pi*n*d/wvl))
return np.abs(t)**2
transmission_analytical = slab_analytical(t_slab, n_slab, monitor_lambdas)
fig, ax = plt.subplots(1,figsize=(5,4), tight_layout=True)
ax.plot(monitor_freqs / 1e12, transmission_analytical, 'k', label='analytical')
ax.set_xlabel('frequency (THz)')
ax.set_ylabel('Transmitted')
ax.legend()
plt.show()
Create Simulation¶
slab = td.Structure(
geometry=td.Box(
center=(0, 0, -Lz/2 + spacing + t_slab/2),
size=(td.inf, td.inf, t_slab),
),
medium=mat_slab,
name='slab',
)
# Here we define the planewave source, placed just in advance (towards negative z) of the slab
source = td.PlaneWave(
source_time = td.GaussianPulse(
freq0=freq0,
fwidth=freqw
),
size=(td.inf, td.inf, 0),
center=(0, 0, -Lz/2+spacing/2),
direction='+',
pol_angle=0,
name='planewave',
)
# We are interested in measuring the transmitted flux, so we set it to be an oversized plane.
monitor = td.FluxMonitor(
center = (0, 0, Lz/2 - spacing/2),
size = (td.inf, td.inf, 0),
freqs = monitor_freqs,
name='flux',
)
# simulation with default 12 PML layers
pml = td.Boundary.pml()
bspec = td.BoundarySpec(x=td.Boundary.periodic(), y = td.Boundary.periodic(), z=pml)
sim_default = td.Simulation(
center = (0, 0, 0),
size = sim_size,
grid_spec = td.GridSpec.uniform(dl=dl),
structures = [slab],
sources = [source],
monitors = [monitor],
run_time = t_stop,
boundary_spec = bspec,
)
# simulation with only 2 PML layers
pml_2layer = td.Boundary.pml(num_layers=2)
bspec_2layer = td.BoundarySpec(x=td.Boundary.periodic(), y = td.Boundary.periodic(), z=pml_2layer)
sim_2layer = sim_default.copy(update={'boundary_spec':bspec_2layer})
Run the simulations¶
sim_FP_data = web.run(sim_default, task_name='lecture06_FP_default_PML', path=f'data/data_FP_default.hdf5')
sim_FP_2layer_data = web.run(sim_2layer, task_name='lecture06_FP_2PMLlayer', path=f'data/data_FP_2layer.hdf5')
Postprocess and Plot¶
Now, we compute the transmitted flux and plot the transmission spectrum and compare them to the analytics.
fig, ax = plt.subplots(1,figsize=(6,4), tight_layout=True)
transmission_default = sim_FP_data['flux'].flux
transmission_2layer = sim_FP_2layer_data['flux'].flux
ax.plot(monitor_freqs/1e12, transmission_analytical, 'k', label='analytical')
ax.plot(monitor_freqs/1e12, transmission_default, 'b--', label='default')
ax.plot(monitor_freqs/1e12, transmission_2layer, 'r--', label='2layer')
ax.set_xlabel('frequency (THz)')
ax.set_xlim([200, 400])
ax.set_ylim(0,1.8)
ax.set_ylabel('Transmitted')
plt.show()
As shown in the plot, the transmission result with the default number of PML layers agrees well with the analytics; however, when there are 2 PML layers, the reflection at the PML boundary results in lots of ripples in the transmission spectrum. Therefore, suppressing parasitic reflection from the boundary with more PML layers is extremely important for obtaining high-quality simulation result.
