TIDY3D
LEARNING CENTER

Lectures

Tutorials

Applications

Introduction to perfectly matched layer (PML)

By Weiliang Jin, Zongfu Yu and Shanhui Fan

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 in FDTD simulations.

- 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.

Share On:

Visit public projects:

Download .ipynb View Presentation Slides

Additional information: This was updated in Apr 13, 2022

FDTD_Python_Tutorial_06

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.

Setup¶

First we import Tidy3d and the other packages needed. If it is not installed, please perform the installation and account authentication by following the detailed instructions here.

In [1]:

                                        # 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.

In [2]:

                                        # 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¶

In [3]:

                                        # 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.

In [4]:

                                        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¶

In [5]:

                                        sim_data = web.run(sim, task_name='lecture06_pml_reflection', path=f'data/data_pml_reflection.hdf5')

                                      

[11:43:02] Created task 'lecture06_pml_reflection' with task_id                 webapi.py:139
           'fdve-37831c8b-1d55-47b9-8500-928ccb8a2f50v1'.

[11:43:07] status = queued                                                      webapi.py:269

[11:43:12] status = preprocess                                                  webapi.py:263

[11:43:18] Maximum FlexCredit cost: 0.025. Use 'web.real_cost(task_id)' to get  webapi.py:286
           the billed FlexCredit cost after a simulation run.

           starting up solver                                                   webapi.py:290

[11:43:19] running solver                                                       webapi.py:300

[11:43:32] early shutoff detected, exiting.                                     webapi.py:313

           status = postprocess                                                 webapi.py:330

[11:43:37] status = success                                                     webapi.py:337

[11:43:40] loading SimulationData from data/data_pml_reflection.hdf5            webapi.py:509

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).

In [6]:

                                        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.

In [7]:

                                        # 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])

                                      

In [8]:

                                        # 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.

In [9]:

                                        # 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¶

In [10]:

                                        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¶

In [11]:

                                        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',
)

                                      

In [12]:

                                        # 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¶

In [13]:

                                        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')

[11:43:40] Created task 'lecture06_FP_default_PML' with task_id                 webapi.py:139
           'fdve-c9b48ff3-8ce7-40d6-9288-110878939cb6v1'.

[11:43:42] status = queued                                                      webapi.py:269

[11:43:46] status = preprocess                                                  webapi.py:263

[11:43:53] Maximum FlexCredit cost: 0.025. Use 'web.real_cost(task_id)' to get  webapi.py:286
           the billed FlexCredit cost after a simulation run.

           starting up solver                                                   webapi.py:290

           running solver                                                       webapi.py:300

[11:44:01] early shutoff detected, exiting.                                     webapi.py:313

           status = postprocess                                                 webapi.py:330

[11:44:06] status = success                                                     webapi.py:337

[11:44:07] loading SimulationData from data/data_FP_default.hdf5                webapi.py:509

[11:44:07] Created task 'lecture06_FP_2PMLlayer' with task_id                   webapi.py:139
           'fdve-ade4dac9-e2b1-4aa2-a64a-393e022a788bv1'.

[11:44:08] status = queued                                                      webapi.py:269

[11:44:13] status = preprocess                                                  webapi.py:263

[11:44:18] Maximum FlexCredit cost: 0.025. Use 'web.real_cost(task_id)' to get  webapi.py:286
           the billed FlexCredit cost after a simulation run.

           starting up solver                                                   webapi.py:290

           running solver                                                       webapi.py:300

[11:44:34] early shutoff detected, exiting.                                     webapi.py:313

[11:44:35] status = postprocess                                                 webapi.py:330

[11:44:38] status = success                                                     webapi.py:337

[11:44:39] loading SimulationData from data/data_FP_2layer.hdf5                 webapi.py:509

Postprocess and Plot¶

Now, we compute the transmitted flux and plot the transmission spectrum and compare them to the analytics.

In [14]:

                                        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.

In [ ]:

Download Intro to PML

TIDY3D LEARNING CENTER

Introduction to perfectly matched layer (PML)

Tutorial 6: Truncate computational domain with PML¶

Setup¶

Creat Simulation¶

Run the simulations¶

Post Run Analysis¶

Impact of parasitic PML reflection to Fabry-Perot cavity simulation¶

Analytical Method¶

Create Simulation¶

Run the simulations¶

Postprocess and Plot¶

TIDY3D
LEARNING CENTER