Tidy3D first walkthrough

See on github, run on colab, or just follow along with the output below.

Our first tutorial focuses on illustrating the basic setup, run, and analysis of a Tidy3D simulation. In this example, we will simulate a plane wave impinging on dielectric slab with a triangular pillar made of a lossy dielectric sitting on top. First, we import everything needed.

# get the most recent version of tidy3d
!pip install -q --upgrade tidy3d

# make sure notebook plots inline
%matplotlib inline

# standard python imports
import numpy as np
import matplotlib.pyplot as plt
import h5py

# tidy3D import
import tidy3d as td
from tidy3d import web

First, we initialize some general simulation parameters. We note that the PML layers extend beyond the simulation domain, making the total simulation size larger - as opposed to some solvers in which the PML is covering part of the user-defined simulation domain.

# Resolution in x, y, z (points per micron)
resolution = 20

# Simulation domain size (in micron)
sim_size = [4, 4, 4]

# Central frequency and bandwidth of pulsed excitation, in Hz
fcen = 2e14
fwidth = 1e13

# Number of PML layers to use along each of the three directions.
pml_layers = [12, 12, 12]

The run time of a simulation depends a lot on whether there are any long-lived resonances. In our example here, there is no strong resonance. Thus, we do not need to run the simulation much longer than after the sources have decayed. We thus set the run time based on the source bandwidth.

# Total time to run in seconds
run_time = 2/fwidth

Structures and materials

Next, we initialize the simulated structure. The structure consists of two geometric objects. Each object is made of a material. Note that the size of any object (structure, source, or monitor) can extend beyond the simulation domain, and is truncated at the edges of that domain. For best results with PML layers, structures must extend all the way through the PMLs. In some such cases, an “infinite” size td.inf can be used to define the size.

# Lossless dielectric
material1 = td.Medium(epsilon=6.)
# Lossy dielectric defined from the real and imaginary part of the refractive index
material2 = td.Medium(n=1.5, k=0.1, freq=fcen)

# Rectangular slab
box = td.Box(center=[0, 0, 0], size=[td.inf, td.inf, 1], material=material1)
# Triangle in the xy-plane with a finite extent in z
equi_tri_verts = [[-1/2, -1/4],
                  [1/2, -1/4],
                  [0, np.sqrt(3)/2 - 1/4]]
poly = td.PolySlab(


Next, we define a source injecting a normal-incidence plane-wave from above. The time dependence of the source is a Gaussian pulse. A source can be added to multiple simulations. After we add the source to a specific simulation, such that the total run time is known, we can use in-built plotting tools to visualize its time- and frequency-dependence, which we will show below.

psource = td.PlaneWave(
    source_time = td.GaussianPulse(
Tidy3D PlaneWave:
name            = None
injection_axis  = -z
position        = 1.50
source_time     = GaussianPulse(
    frequency  = 2.00e+14,
    fwidth     = 1.00e+13,
    offset     = 5.00)
polarization    = y
amplitude       = 1.00e+00


Finally, we can also add some monitors that will record the fields that we request during the simulation run. The two main types are time monitors that record the time-domain fields, and frequency monitors that record a running discrete Fourier transform of the fields at a given set of frequencies. Time monitors are best used to monitor the time dependence of the fields at a single point. Spatially large time monitors can lead to a very large amount of data that needs to be stored. Frequency monitors on the other hand are great for investigating the steady-state field distribution in 2D or even 3D regions of the simulation.

time_mnt = td.TimeMonitor(center=[0, 0, 0], size=[0, 0, 0])
freq_mnt1 = td.FreqMonitor(center=[0, 0, -1], size=[20, 20, 0], freqs=[fcen])
freq_mnt2 = td.FreqMonitor(center=[0, 0, 0], size=[20, 0, 20], freqs=[fcen])


Now we can initialize the simulation with all the elements defined above.

# Initialize simulation
sim = td.Simulation(size=sim_size,
                    structures=[box, poly],
                    monitors=[time_mnt, freq_mnt1, freq_mnt2],
Initializing simulation...
Mesh step (micron): [5.00e-02, 5.00e-02, 5.00e-02].
Simulation domain in number of grid points: [104, 104, 104].
Total number of computational grid points: 1.12e+06.
Total number of time steps: 2308.
Estimated data size (GB) of monitor monitor: 0.0001.
Estimated data size (GB) of monitor monitor_1: 0.0005.
Estimated data size (GB) of monitor monitor_2: 0.0005.

We can check the simulation monitors just to make sure everything looks right.

[Tidy3D TimeMonitor: {
name     = None
center   = [0.0000, 0.0000, 0.0000]
size     = [0.0000, 0.0000, 0.0000]
t_start  = 0.00e+00,
t_stop   = None
t_step   = None
Store: ['e', 'h']
, Tidy3D FreqMonitor: {
name     = None
center   = [0.0000, 0.0000, -1.0000]
size     = [20.0000, 20.0000, 0.0000]
freqs    = [2.00e+14]
Store: ['e', 'h']
, Tidy3D FreqMonitor: {
name     = None
center   = [0.0000, 0.0000, 0.0000]
size     = [20.0000, 0.0000, 20.0000]
freqs    = [2.00e+14]
Store: ['e', 'h']

Visualization functions

We can now use the some in-built plotting functions to make sure that we have set up the simulation as we desire.

# Visualize source

To visualize the structures in the simulation, we will plot three cross sections at z=0.75, y=0, and x=0, respectively. The relative permittivity of objects is plotted in greyscale. By default, sources are overlayed in green, monitors in yellow, and PML boundaries in orange. The transparency of the overlay can be controlled using source_alpha, monitor_alpha, and pml_alpha. As an example, in the middle plot (y=0), we turn off the monitor overlay since the monitor covers the whole cross-section, while in the right-most plot, we make the source more visible.

# Set the same color scale in all plots
clim = (1, 6)

fig, ax = plt.subplots(1, 3, figsize=(13, 4))
sim.viz_eps_2D(normal='z', position=0.75, ax=ax[0], clim=clim);
sim.viz_eps_2D(normal='y', ax=ax[1], clim=clim, monitor_alpha=0.);
sim.viz_eps_2D(normal='x', ax=ax[2], clim=clim, source_alpha=0.9);
<AxesSubplot:title={'center':'yz-plane at x=0.00e+00um'}, xlabel='y (um)', ylabel='z (um)'>

Alternatively, we can also plot the structures with a fake color based on the material they are made of.

fig, ax = plt.subplots(1, 3, figsize=(12, 3))
sim.viz_mat_2D(normal='z', position=0.75, ax=ax[0]);
sim.viz_mat_2D(normal='y', ax=ax[1], monitor_alpha=0);
sim.viz_mat_2D(normal='x', ax=ax[2], source_alpha=0.9, legend=True);
<AxesSubplot:title={'center':'yz-plane at x=0.00e+00um'}, xlabel='y (um)', ylabel='z (um)'>

Running through the web API

Now that the simulation is constructed, we can run it using the web API of Tidy3D. First, we submit the project. Note that we can give it a custom name.

project = web.new_project(sim.export(), task_name='Simulation test')
Uploading the json file...

We can monitor the status of all our projects by listing them in chronological order of submission.

# Check status (of the last 5 submitted projects)
Project name: Simulation test-16b87539b3ba6838-16b87539ed2b41a0
Task ID     : 28313f5f-ec38-4bd4-bc5b-ed3104425903
Submit time : 2021:11:17:22:15:56
Status      : queued
Project name: Task_1636051009-16b46bba1ae99b38
Task ID     : d92cd1ac-6d33-4164-9f86-038885b312be
Submit time : 2021:11:04:18:36:51
Status      : success
Project name: Task_1636049328-16b46a32d86db088
Task ID     : b55fedd8-0407-412f-bb77-184894f19d86
Submit time : 2021:11:04:18:08:51
Status      : error
Project name: Task_1636048673-16b4699a6af6d8d8
Task ID     : 419c4139-4e53-42cc-99e2-698881d60d07
Submit time : 2021:11:04:17:57:56
Status      : deleted
Project name: Task_1635993296-16b4373cf414f7e8
Task ID     : 1de9f277-5c03-40fc-8b93-1f4d1310f6eb
Submit time : 2021:11:04:02:34:59
Status      : deleted

Or, we can continously monitor the status of the current project, and wait until the run is successful. The monitor_project() function will keep running until either a 'success' or 'error' status is returned.

Project 'Simulation test-16b87539b3ba6838-16b87539ed2b41a0' status: success...

Loading and analyzing data

After a successful run, we can download the results and load them into our simulation model. We use the download_results function from our web API, which downloads a single hdf5 file containing all the monitor data, a log file, and a json file defining the original simulation (same as what you’ll get if you run sim.export_json() on the current object). Optionally, you can provide a folder in which to store the files. In the example below, the results are stored in the out folder.

web.download_results(project['taskId'], target_folder='output/')
# Show the output of the log file
with open("output/tidy3d.log") as f:
Simulation domain Nx, Ny, Nz: [104, 104, 104]
Applied symmetries: [0, 0, 0]
Number of computational grid points: 1.1249e+06.
Using subpixel averaging: True
Number of time steps: 2308
Automatic shutoff factor: 1.00e-05
Time step (s): 8.6662e-17

Compute source modes time (s):     0.0972
Compute monitor modes time (s):    0.1224

Rest of setup time (s):            0.0605

Starting solver...
- Time step     92 / time 7.97e-15s (  4 % done), field decay: 1.00e+00
- Time step    184 / time 1.59e-14s (  8 % done), field decay: 1.00e+00
- Time step    276 / time 2.39e-14s ( 12 % done), field decay: 1.00e+00
- Time step    369 / time 3.20e-14s ( 16 % done), field decay: 1.00e+00
- Time step    461 / time 4.00e-14s ( 20 % done), field decay: 1.00e+00
- Time step    553 / time 4.79e-14s ( 24 % done), field decay: 1.00e+00
- Time step    646 / time 5.60e-14s ( 28 % done), field decay: 1.00e+00
- Time step    738 / time 6.40e-14s ( 32 % done), field decay: 1.00e+00
- Time step    830 / time 7.19e-14s ( 36 % done), field decay: 1.00e+00
- Time step    918 / time 7.96e-14s ( 39 % done), field decay: 1.00e+00
- Time step    923 / time 8.00e-14s ( 40 % done), field decay: 1.00e+00
- Time step   1015 / time 8.80e-14s ( 44 % done), field decay: 1.00e+00
- Time step   1107 / time 9.59e-14s ( 48 % done), field decay: 8.93e-01
- Time step   1200 / time 1.04e-13s ( 52 % done), field decay: 3.56e-01
- Time step   1292 / time 1.12e-13s ( 56 % done), field decay: 1.34e-01
- Time step   1384 / time 1.20e-13s ( 60 % done), field decay: 4.49e-02
- Time step   1477 / time 1.28e-13s ( 64 % done), field decay: 9.47e-03
- Time step   1569 / time 1.36e-13s ( 68 % done), field decay: 1.50e-03
- Time step   1661 / time 1.44e-13s ( 72 % done), field decay: 2.60e-04
- Time step   1754 / time 1.52e-13s ( 76 % done), field decay: 5.43e-05
- Time step   1846 / time 1.60e-13s ( 80 % done), field decay: 1.15e-05
- Time step   1938 / time 1.68e-13s ( 84 % done), field decay: 2.11e-06
Field decay smaller than shutoff factor, exiting solver.

Solver time (s):                   3.9134
Post-processing time (s):          0.0232

Load a simulation from JSON file

In case we are returning to a project for which we do not have an active simulation object, we can re-initialize the simulation from the freshly downloaded json file. Below, we load the file from the project we just ran, and load the ouput data downloaded from the server.

NB: when running the notebook from the beginning, it is not needed to call import_json(), as the current Simulation is stored in memory. Here we do it to illustrate how this function can be used, essentially re-initializing the same Simulation.

sim = td.Simulation.import_json("output/simulation.json")
Initializing simulation...
Mesh step (micron): [5.00e-02, 5.00e-02, 5.00e-02].
Simulation domain in number of grid points: [104, 104, 104].
Total number of computational grid points: 1.12e+06.
Total number of time steps: 2308.
Estimated data size (GB) of monitor monitor: 0.0001.
Estimated data size (GB) of monitor monitor_1: 0.0005.
Estimated data size (GB) of monitor monitor_2: 0.0005.
Applying source normalization to all frequency monitors using source index 0.

Visualization functions

Finally, we can now use the in-built visualization tools to examine the results. Below, we plot the y-component of the field recorded by the two frequency monitors (this is the dominant component since the source is y-polarized).

monitors = sim.monitors

fig, ax = plt.subplots(1, 2, figsize=(10, 4))
sim.viz_field_2D(monitors[1], ax=ax[0], cbar=True, comp='y', val='re');
sim.viz_field_2D(monitors[2], ax=ax[1], cbar=True, comp='y', val='re');
<matplotlib.collections.QuadMesh at 0x7f132868b2e8>

Monitor data in Simulation object and in hdf5 file

The raw data is also accessible as numpy arrays and can be queried using sim.data(monitor).

mdata = sim.data(monitors[1])
print("Shape of E-field array stored in monitor 1: ", mdata['E'].shape)
print("Size of corresponding mesh in x           : ", mdata['xmesh'].size)
print("Size of corresponding mesh in y           : ", mdata['ymesh'].size)
print("Size of corresponding mesh in z           : ", mdata['zmesh'].size)
print("Number of frequencies                     : ", mdata['freqs'].size)
Shape of E-field array stored in monitor 1:  (3, 104, 104, 1, 1)
Size of corresponding mesh in x           :  104
Size of corresponding mesh in y           :  104
Size of corresponding mesh in z           :  1
Number of frequencies                     :  1

We can use this raw data for example to also plot the time-domain fields recorded in the TimeMonitor, which look largely like a delayed version of the source input, indicating that no resonant features were excited.

fig, ax = plt.subplots(1)
tdata = sim.data(monitors[0])
tmesh = tdata['tmesh']
ax.plot(tmesh, tdata['E'][1, 0, 0, 0, :])
ax.set_xlim(tmesh[0], tmesh[-1])
ax.set_xlabel("Time [s]")
ax.set_ylabel("Ey [a. u.]")

Finally, the raw data is also accessible through the hdf5 file. Still, note that the preferred method to access the data is through sim.data(). If you do access through the hdf5 file, the spatial mesh is not directly available. Instead, the 'indspan' dataset is an array of shape (3, 2) that gives the (beginning, end) index of the points of the simulation grid at which the fields were computed.

# Open the data file and print stored groups and datasets
data_file = h5py.File("output/monitor_data.hdf5", "r")
print("Data groups stored in file:   ", list(data_file.keys()))
print("Datasets in 'monitor_2': ", list(data_file['monitor_2'].keys()))

# Read the E-field of the second FreqMonitor
E_hdf5 = np.array(data_file['monitor_2']['E'])
indspan = np.array(data_file['monitor_2']['indspan'])
x_hdf5 = sim.grid.mesh[0][indspan[0, 0]:indspan[0, 1]]
y_hdf5 = sim.grid.mesh[1][indspan[1, 0]:indspan[1, 1]]
z_hdf5 = sim.grid.mesh[2][indspan[2, 0]:indspan[2, 1]]

# Close the file
Data groups stored in file:    ['diverged', 'monitor', 'monitor_1', 'monitor_2']
Datasets in 'monitor_2':  ['E', 'H', 'freqs', 'indspan', 'tmesh']

Permittivity data

We can also query the relative permittivity in the simulation within a volume of a given center and size. The method Simulation.epsilon returns the permittivity at all Yee grid locations inside the specified volume. Instead of manually setting a volume, we can also input a monitor, in which case the fields in that monitor volume are returned.

eps, mesh = sim.epsilon(freq_mnt2)

The permittivity data is returned at the same locations that the E field at which the E field is stored in the monitor.

print("Shape of permittivity array (x, y, z)              : ", eps.shape)
print("Shape of E-field (polarization, x, y, z, frequency): ", E_hdf5.shape)
Shape of permittivity array (x, y, z)              :  (104, 1, 104)
Shape of E-field (polarization, x, y, z, frequency):  (3, 104, 1, 104, 1)

We can use this data to approximately reproduce the figure we get from the built-in visualization. Note that in the in-build function, permittivity of 1 (vacuum) is automatically made fully transparent.

fig, ax = plt.subplots(1, 2, figsize=(10, 4))
# Built-in Tidy3D viz
sim.viz_field_2D(sim.monitors[2], ax=ax[0], cbar=True, comp='y', val='re', eps_alpha=0.3);

# Manually from data loaded from the hdf5 file
ax[1].imshow(np.real(E_hdf5[1, :, 0, :, 0]).T, cmap='RdBu', origin='lower',
                 extent=[x_hdf5[0], x_hdf5[-1], z_hdf5[0], z_hdf5[-1]])
# Overlay the permittivity data
ax[1].imshow(np.real(eps[:, 0, :]).T, alpha=0.3,
                 extent=[mesh[0][0], mesh[0][-1], mesh[2][0], mesh[2][-1]],
                 origin='lower', cmap='Greys')