Setting up boundary conditions
Contents
Setting up boundary conditions#
This notebook will give a tutorial on setting up boundary conditions in Tidy3D.
If you are new to the finite-difference time-domain (FDTD) method, we highly recommend going through our FDTD101 tutorials. For simulation examples, please visit our examples page. FDTD simulations can diverge due to various reasons. If you run into any simulation divergence issues, please follow the steps outlined in our troubleshooting guide to resolve it.
[1]:
# standard python imports
import numpy as np
import matplotlib.pyplot as plt
# tidy3d imports
import tidy3d as td
import tidy3d.web as web
Define Simulation Parameters#
First, we’ll define some basc simulation parameters, the size of the domain, and the discretization resolution.
[2]:
# Define material properties
medium = td.Medium(permittivity=2)
wavelength = 1
f0 = td.C_0 / wavelength / np.sqrt(medium.permittivity)
# Set the domain size in x, y, and z
domain_size = 12 * wavelength
# create the geometry
geometry = []
# construct simulation size array
sim_size = (domain_size, domain_size, domain_size)
# Bandwidth in Hz
fwidth = f0 / 40.0
# Gaussian source offset; the source peak is at time t = offset/fwidth
offset = 4.0
# time dependence of sources
source_time = td.GaussianPulse(freq0=f0, fwidth=fwidth, offset=offset)
# Simulation run time past the source decay (around t=2*offset/fwidth)
run_time = 50 / fwidth
Create sources and monitors#
To study the effect of the various boundary conditions, we’ll define a point dipole source and a series of frequency- and time-domain monitors in the volume of the simulation domain and at its edges.
[3]:
# create a point dipole source
dipole = td.PointDipole(
center=(0, 0, 0),
source_time=source_time,
polarization="Ex",
)
# these monitors will be used to plot fields on planes through the middle of the domain in the frequency domain
monitor_xz_freq = td.FieldMonitor(
center=(0, 0, 0), size=(domain_size, 0, domain_size), freqs=[f0], name="xz_freq"
)
monitor_yz_freq = td.FieldMonitor(
center=(0, 0, 0), size=(0, domain_size, domain_size), freqs=[f0], name="yz_freq"
)
monitor_xy_freq = td.FieldMonitor(
center=(0, 0, 0), size=(domain_size, domain_size, 0), freqs=[f0], name="xy_freq"
)
# this monitor will be used to plot fields on a plane through the middle of the domain in the time domain
monitor_xz_time = td.FieldTimeMonitor(
center=(0, 0, 0), size=(domain_size, 0, domain_size), interval=50, name="xz_time"
)
[10:18:30] WARNING: Default value for the field monitor monitor.py:261 'colocate' setting has changed to 'True' in Tidy3D 2.4.0. All field components will be colocated to the grid boundaries. Set to 'False' to get the raw fields on the Yee grid instead.
WARNING: Default value for the field monitor monitor.py:261 'colocate' setting has changed to 'True' in Tidy3D 2.4.0. All field components will be colocated to the grid boundaries. Set to 'False' to get the raw fields on the Yee grid instead.
WARNING: Default value for the field monitor monitor.py:261 'colocate' setting has changed to 'True' in Tidy3D 2.4.0. All field components will be colocated to the grid boundaries. Set to 'False' to get the raw fields on the Yee grid instead.
WARNING: Default value for the field monitor monitor.py:261 'colocate' setting has changed to 'True' in Tidy3D 2.4.0. All field components will be colocated to the grid boundaries. Set to 'False' to get the raw fields on the Yee grid instead.
Boundary specifications#
A BoundarySpec object defines the boundary conditions applied on each of the 6 domain edges, and is provided as an input to the simulation. In the following sections, we’ll explore several different features within BoundarySpec and different ways of defining it.
A BoundarySpec consists of three Boundary objects, each defining the boundaries on the plus and minus side of each dimension. A number of convenience functions are available to quickly define various types of boundaries, which will be demonstrated below.
Example 1: Default PML boundaries along some dimensions#
In most cases, one just wants to specify whether there are absorbing PML layers along any of the x, y, z dimensions. By default, Tidy3D simulations have PML boundaries on all sides.
[4]:
# initialize the simulation object with default boundaries (PML on all sides)
sim = td.Simulation(
size=sim_size,
sources=[dipole],
monitors=[monitor_xz_time],
run_time=run_time,
)
# Visualize the geometry
fig, ax1 = plt.subplots(figsize=(4, 4))
sim.plot(y=0, ax=ax1)
plt.show()
If we want to explicitly set the boundaries, we can use the all_sides() method. This can be used to set any type of boundary condition on all sides of the simulation.
[5]:
# define a basic boundary spec setting PML in all directions
bspec_pml = td.BoundarySpec.all_sides(boundary=td.PML())
# initialize the simulation object with the above boundary spec and source
sim = td.Simulation(
size=sim_size,
sources=[dipole],
monitors=[monitor_xz_time],
run_time=run_time,
boundary_spec=bspec_pml,
)
# Visualize the geometry
fig, ax1 = plt.subplots(figsize=(4, 4))
sim.plot(y=0, ax=ax1)
plt.show()
We can run this simulation and visualize the results. We can observe the effect of the PML by looking at the fields in the time domain as they impinge on the boundaries. The figure shows that the fields are absorbed by the PML layers as expected, and no reflections are observed.
[6]:
sim_data = web.run(
sim, task_name="bc_example1", path="data/bc_example1.hdf5", verbose=True
)
fig, ax = plt.subplots(tight_layout=True, figsize=(5, 4))
sim_data.plot_field(
field_monitor_name="xz_time", field_name="Ex", y=0, val="abs", t=2e-13, ax=ax
)
plt.show()
View task using web UI at webapi.py:190 'https://tidy3d.simulation.cloud/workbench?taskId=fdve- 3077d877-17e7-4314-a8aa-78ec653ff3dfv1'.
We can also set PML on specified sides only by calling the BoundarySpec.pml() method, e.g. BoundarySpec.pml(x=False, y=False, z=False). This will put PML along dimensions where the argument is True, and set periodic boundaries along the other dimensions.
Let’s try putting pml on only the x edge.
[7]:
# define a basic boundary spec setting PML in along x only
bspec_pml = td.BoundarySpec.pml(x=True)
sim = td.Simulation(
size=sim_size,
sources=[dipole],
monitors=[monitor_xz_time],
run_time=run_time,
boundary_spec=bspec_pml,
)
ax = sim.plot(y=0)
plt.show()
Example 2: different boundaries on different edges#
Here, we’ll place a PML along the y and z directions. In the x direction, we’ll place a PML on the left side (x-minus), and a PEC on the right (x-plus). To specify individual boundary conditions along different dimensions, instead of BoundarySpec, the class Boundary is used, which defines the plus and minus boundaries along a single dimension.
We’ll test this set of boundaries by placing an angled Gaussian beam near the lower edge of the domain and observing the field bounce off the PEC in x-plus.
[8]:
# this defines a 1D boundary with a PML on both the plus and minus sides, which can be reused for both y and z directions
boundary_yz = td.Boundary.pml(num_layers=15)
# this defines a 1D boundary with a PML on the minus side and a PEC on the plus side
boundary_x = td.Boundary(minus=td.PML(), plus=td.PECBoundary())
# now just set these in the boundary spec along the appropriate dimensions
bspec_pml_pec = td.BoundarySpec(x=boundary_x, y=boundary_yz, z=boundary_yz)
# create the Gaussian beam source
buffer_source = (
domain_size / 10
) # distance between the source and the bottom of the domain
gaussian_beam = td.GaussianBeam(
center=(0, 0, -domain_size / 2 + buffer_source),
size=(td.inf, td.inf, 0),
source_time=source_time,
direction="+",
pol_angle=0,
angle_theta=np.pi / 4.0,
angle_phi=np.pi / 8.0,
waist_radius=wavelength * 2,
waist_distance=-wavelength * 4,
)
# initialize the simulation object with the above boundary spec and source
sim = td.Simulation(
size=sim_size,
sources=[gaussian_beam],
monitors=[monitor_xz_time],
run_time=run_time,
boundary_spec=bspec_pml_pec,
)
We can verify that the PML is placed only on the left hand side and not on the right.
[9]:
# Visualize the geometry
fig, ax1 = plt.subplots(figsize=(4, 4))
sim.plot(y=0, ax=ax1)
plt.show()
Run Simulation#
[10]:
sim_data = web.run(
sim, task_name="bc_example2", path="data/bc_example2.hdf5", verbose=True
)
View task using web UI at webapi.py:190 'https://tidy3d.simulation.cloud/workbench?taskId=fdve- fb8fa337-e05a-46d6-b609-cb37f6fd2470v1'.
Visualize results#
We can observe the effect of the PEC on x-plus by looking at the fields in the time domain as they bounce off the PEC boundary, as shown in the figure. Furthermore, the z-component of the electric field, which is tangential to the PEC boundary, goes to 0 at the boundary as expected.
[11]:
fig, (ax1, ax2) = plt.subplots(1, 2, tight_layout=True, figsize=(10, 4))
sim_data.plot_field(
field_monitor_name="xz_time", field_name="Ex", y=0, val="abs", t=2e-13, ax=ax1
)
sim_data.plot_field(
field_monitor_name="xz_time", field_name="Ez", y=0, val="abs", t=2e-13, ax=ax2
)
plt.show()
Example 3: different boundaries on different edges#
Next, let’s consider an even more general setup where all 6 boundaries of the domain are individually specified, and different types of boundaries are mixed in the plus and minus sides along each dimension.
We’ll test this set of boundaries by placing the point dipole at the center of the domain again, this time studying the fields at the edges of the domain and checking that the various tangential field components satisfy the PEC and PMC boundary conditions as expected from theory.
Note that when a periodic boundary is applied across from a PEC (PMC), the periodic boundary is just replaced by the PEC (PMC). This is an equivalent setup but a warning is issued to avoid confusion.
[12]:
# define the boundary spec
bspec_general = td.BoundarySpec(
x=td.Boundary(minus=td.Periodic(), plus=td.PECBoundary()),
y=td.Boundary(minus=td.PECBoundary(), plus=td.PMCBoundary()),
z=td.Boundary(minus=td.Absorber(), plus=td.PML()),
)
# initialize the simulation object with the above boundary spec and the dipole source from earlier
sim = td.Simulation(
size=sim_size,
sources=[dipole],
monitors=[monitor_xz_freq, monitor_yz_freq, monitor_xy_freq],
run_time=run_time,
boundary_spec=bspec_general,
)
# Visualize the geometry
fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(11, 4))
sim.plot(x=0, ax=ax1)
sim.plot(y=0, ax=ax2)
sim.plot(z=0, ax=ax3)
plt.show()
[10:19:24] WARNING: A periodic boundary condition was specified boundary.py:400 on the opposite side of a perfect electric or magnetic conductor boundary. This periodic boundary condition will be replaced by the perfect electric or magnetic conductor across from it.
Run Simulation#
[13]:
sim_data = web.run(
sim, task_name="bc_example3", path="data/bc_example3.hdf5", verbose=True
)
View task using web UI at webapi.py:190 'https://tidy3d.simulation.cloud/workbench?taskId=fdve- 8c9d7a27-d931-481c-8c6b-25fb98b492f7v1'.
Visualize results#
In x-plus, we have a PEC, so Ey and Ez should be nearly zero on the x-plus boundary.
In x-minus we have a periodic condition across from the PEC in x-plus, which just means that there’s also a PEC in x-minus. So again Ey and Ez should be nearly zero on the x-minus boundary.
In y-minus, we have another PEC, so Ex and Ez should be nearly zero on the y-minus boundary.
In y-plus, we have a PMC, so Hx and Hz should be nearly zero on the y-plus boundary.
Each of these cases is verified below.
[14]:
fig, (ax1, ax2) = plt.subplots(1, 2, tight_layout=True, figsize=(10, 4))
sim_data.plot_field(
field_monitor_name="xz_freq", field_name="Ey", y=0, val="abs", f=f0, ax=ax1
)
sim_data.plot_field(
field_monitor_name="xz_freq", field_name="Ez", y=0, val="abs", f=f0, ax=ax2
)
fig.suptitle("Notice that both Ey and Ez go to zero on both sides", fontsize=16)
fig, (ax1, ax2) = plt.subplots(1, 2, tight_layout=True, figsize=(10, 4))
sim_data.plot_field(
field_monitor_name="yz_freq", field_name="Ex", x=0, val="abs", f=f0, ax=ax1
)
sim_data.plot_field(
field_monitor_name="yz_freq", field_name="Ez", x=0, val="abs", f=f0, ax=ax2
)
fig.suptitle(
"Notice that both Ex and Ez go to zero on the left (PEC) side", fontsize=16
)
fig, (ax1, ax2) = plt.subplots(1, 2, tight_layout=True, figsize=(10, 4))
sim_data.plot_field(
field_monitor_name="yz_freq", field_name="Hx", x=0, val="abs", f=f0, ax=ax1
)
sim_data.plot_field(
field_monitor_name="yz_freq", field_name="Hz", x=0, val="abs", f=f0, ax=ax2
)
fig.suptitle(
"Notice that both Hx and Hz go to zero on the right (PMC) side", fontsize=16
)
plt.show()
Example 4: Bloch boundary conditions along x and y to allow injecting an angled plane wave source#
Finally, we’ll use the Bloch boundary condition in combination with a plane wave source to demonstrate the injection of an angled plane wave. The Bloch boundaries are used along the x and y dimensions, while PMLs are used along z. The Bloch vectors are automatically computed based on the angles defined in the plane wave source. A dielectric background medium with a relative permittivity of 2 is used.
We’ll test this set of boundaries by placing an angled plane wave near the lower edge of the domain and observing the plane wave in the frequency domain.
[15]:
# First, define the plane wave source, since it is needed to define the Bloch boundary in this case.
# Note that in general, the Bloch boundary can also be defined by just providing a bandstructure-normalized Bloch vector.
buffer_source = (
domain_size / 10
) # distance between the source and the bottom of the domain
plane_wave = td.PlaneWave(
center=(0, 0, -domain_size / 2 + buffer_source),
size=(td.inf, td.inf, 0),
source_time=source_time,
direction="+",
pol_angle=0,
angle_theta=np.pi / 3.0,
angle_phi=np.pi / 6.0,
)
# create the Bloch boundaries
bloch_x = td.Boundary.bloch_from_source(
source=plane_wave, domain_size=sim_size[0], axis=0, medium=medium
)
bloch_y = td.Boundary.bloch_from_source(
source=plane_wave, domain_size=sim_size[1], axis=1, medium=medium
)
bspec_bloch = td.BoundarySpec(x=bloch_x, y=bloch_y, z=td.Boundary.pml())
# initialize the simulation object with the above boundary spec and source
sim = td.Simulation(
size=sim_size,
sources=[plane_wave],
monitors=[monitor_xz_freq, monitor_yz_freq, monitor_xy_freq],
run_time=run_time,
boundary_spec=bspec_bloch,
medium=medium,
)
# Visualize the geometry
fig, ax1 = plt.subplots(figsize=(4, 4))
sim.plot(y=0, ax=ax1)
plt.show()
Run Simulation#
[16]:
sim_data = web.run(sim, task_name="bc_example4", path="data/bc_example4.hdf5")
View task using web UI at webapi.py:190 'https://tidy3d.simulation.cloud/workbench?taskId=fdve- 33eec131-5063-4f3d-9dc9-e94062ff1e5av1'.
Visualize results#
In shown in the plots, the plane wave is successfully injected at the specified angles without reflections.
[17]:
fig, (ax1, ax2, ax3) = plt.subplots(1, 3, tight_layout=True, figsize=(13, 4))
sim_data.plot_field(
field_monitor_name="xz_freq", field_name="Ey", y=0, val="real", f=f0, ax=ax1
)
sim_data.plot_field(
field_monitor_name="xz_freq", field_name="Ey", y=0, val="real", f=f0, ax=ax2
)
sim_data.plot_field(
field_monitor_name="xz_freq", field_name="Ez", y=0, val="real", f=f0, ax=ax3
)
fig, (ax1, ax2, ax3) = plt.subplots(1, 3, tight_layout=True, figsize=(13, 4))
sim_data.plot_field(
field_monitor_name="yz_freq", field_name="Ex", x=0, val="real", f=f0, ax=ax1
)
sim_data.plot_field(
field_monitor_name="yz_freq", field_name="Ey", x=0, val="real", f=f0, ax=ax2
)
sim_data.plot_field(
field_monitor_name="yz_freq", field_name="Ez", x=0, val="real", f=f0, ax=ax3
)
fig, (ax1, ax2, ax3) = plt.subplots(1, 3, tight_layout=True, figsize=(13, 4))
sim_data.plot_field(
field_monitor_name="xy_freq", field_name="Ex", z=0, val="real", f=f0, ax=ax1
)
sim_data.plot_field(
field_monitor_name="xy_freq", field_name="Ey", z=0, val="real", f=f0, ax=ax2
)
sim_data.plot_field(
field_monitor_name="xy_freq", field_name="Ez", z=0, val="real", f=f0, ax=ax3
)
fig.suptitle("Observe the phi=pi/6 angle w.r.t. the x axis", fontsize=16)
plt.show()
Summary#
The fastest way to define boundary conditions with periodic and absorbing boundaries on some sides is through
td.BoundarySpec.pml(y=True)(i.e., for PML along y only), or throughtd.BoundarySpec.all_sides(boundary=PML())to set the same boundary on all sides (in this example, PML on all sides).To further fine tune the global boundary specifications, use BoundarySpec(x, y, z), where x, y, z are Boundary objects.
Boundary objects have some convenience methods, like
.pml()for defining the boundaries along an edge with more parameters.Boundary objects can be further specified by supplying boundary conditions to the
plusandminusfields.
[ ]: