Tidy3D first walkthrough#

Run this notebook in your browser using Binder.

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.

[1]:
# standard python imports
import numpy as np
import matplotlib.pyplot as plt
import h5py

# tidy3d imports
import tidy3d as td
from tidy3d import web

[09:48:11] WARNING  This version of Tidy3D was pip installed from the 'tidy3d-beta' repository on   __init__.py:102
                    PyPI. Future releases will be uploaded to the 'tidy3d' repository. From now on,                
                    please use 'pip install tidy3d' instead.                                                       
           INFO     Using client version: 1.8.3                                                     __init__.py:120

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.

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

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

# # pad the z direction (out of plane) with PML but leave it off x and y for periodic boundary conditions.
# boundary_spec=td.BoundarySpec.pml(x=False, y=False, z=True)

# apply a PML in all directions
boundary_spec = td.BoundarySpec.all_sides(boundary=td.PML())

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.

[3]:
# Total time to run in seconds
run_time = 2 / fwidth

Structures and materials#

Next, we initialize the simulated structure. The structure consists of two Structure objects. Each object consists of a Geometry and a Medium to define the spatial extent and material properties, respectively. 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.

Note: For best results, structures that intersect with the PML or simulation edges should extend extend all the way through. In many such cases, an “infinite” size td.inf can be used to define the size along that dimension.

[4]:
# Lossless dielectric specified directly using relative permittivity
material1 = td.Medium(permittivity=6.0)

# Lossy dielectric defined from the real and imaginary part of the refractive index
material2 = td.Medium.from_nk(n=1.5, k=0.0, freq=freq0)
# material2 = td.Medium(permittivity=2.)


# Rectangular slab, extending infinitely in x and y with medium `material1`
box = td.Structure(
    geometry=td.Box(center=[0, 0, 0], size=[td.inf, td.inf, 1]), medium=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.Structure(
    geometry=td.PolySlab(
        vertices=(2 * np.array(equi_tri_verts)).tolist(),
        # vertices=equi_tri_verts,
        slab_bounds=(0.5, 1.0),
        axis=2,
    ),
    medium=material2,
)

Sources#

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.

[5]:
psource = td.PlaneWave(
    center=(0, 0, 1.5),
    direction="-",
    size=(td.inf, td.inf, 0),
    source_time=td.GaussianPulse(freq0=freq0, fwidth=fwidth),
    pol_angle=np.pi / 2,
)

Monitors#

Finally, we can also add some monitors that will record the fields that we request during the simulation run.

The two monitor types for measuring fields are FieldMonitor and FieldTimeMonitor, which record the frequency-domain and time-domain fields, respectively.

FieldMonitor objects operate by running a discrete Fourier transform of the fields at a given set of frequencies to perform the calculation “in-place” with the time stepping. FieldMonitor objects are useful for investigating the steady-state field distribution in 2D or even 3D regions of the simulation.

FieldMonitor objects are best used to monitor the time dependence of the fields at a single point, but they can also be used to create “animations” of the field pattern evolution. Because spatially large FieldMonitor objects can lead to a very large amount of data that needs to be stored, an optional start and stop time can be supplied, as well as an interval specifying the amount of time steps between each measurement (default of 1).

[6]:
# measure time domain fields at center location, measure every 5 time steps
time_mnt = td.FieldTimeMonitor(
    center=[0, 0, 0], size=[0, 0, 0], interval=5, name="field_time"
)

# measure the steady state fields at central frequency in the xy plane and the xz plane.
freq_mnt1 = td.FieldMonitor(
    center=[0, 0, -1], size=[20, 20, 0], freqs=[freq0], name="field1"
)
freq_mnt2 = td.FieldMonitor(
    center=[0, 0, 0], size=[20, 0, 20], freqs=[freq0], name="field2"
)

Simulation#

Now we can initialize the Simulation with all the elements defined above. A nonuniform simulation grid is generated automatically based on a given minimum number of cells per wavelength in each material (10 by default), using the frequencies defined in the source.

[7]:
# Initialize simulation
sim = td.Simulation(
    size=sim_size,
    grid_spec=td.GridSpec.auto(min_steps_per_wvl=20),
    structures=[box, poly],
    sources=[psource],
    monitors=[time_mnt, freq_mnt1, freq_mnt2],
    run_time=run_time,
    boundary_spec=boundary_spec,
)

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

[8]:
for m in sim.monitors:
    m.help()

╭───────────── <class 'tidy3d.components.monitor.FieldTimeMonitor'> ─────────────╮
 :class:`Monitor` that records electromagnetic fields in the time domain.       
                                                                                
 ╭────────────────────────────────────────────────────────────────────────────╮ 
  FieldTimeMonitor(                                                           
  type='FieldTimeMonitor',                                                
  center=(0.0, 0.0, 0.0),                                                 
  size=(0.0, 0.0, 0.0),                                                   
  name='field_time',                                                      
  start=0.0,                                                              
  stop=None,                                                              
  interval=5,                                                             
  fields=('Ex', 'Ey', 'Ez', 'Hx', 'Hy', 'Hz'),                            
  interval_space=(1, 1, 1),                                               
  colocate=False                                                          
  )                                                                           
 ╰────────────────────────────────────────────────────────────────────────────╯ 
                                                                                
   bounding_box = Box(type='Box', center=(0.0, 0.0, 0.0), size=(0.0, 0.0, 0.0)) 
         bounds = ((0.0, 0.0, 0.0), (0.0, 0.0, 0.0))                            
         center = (0.0, 0.0, 0.0)                                               
       colocate = False                                                         
         fields = ('Ex', 'Ey', 'Ez', 'Hx', 'Hy', 'Hz')                          
       geometry = Box(type='Box', center=(0.0, 0.0, 0.0), size=(0.0, 0.0, 0.0)) 
       interval = 5                                                             
 interval_space = (1, 1, 1)                                                     
           name = 'field_time'                                                  
    plot_params = PlotParams(                                                   
                      alpha=0.4,                                                
                      edgecolor='orange',                                       
                      facecolor='orange',                                       
                      fill=True,                                                
                      hatch=None,                                               
                      linewidth=3.0,                                            
                      type='PlotParams'                                         
                  )                                                             
           size = (0.0, 0.0, 0.0)                                               
          start = 0.0                                                           
           stop = None                                                          
           type = 'FieldTimeMonitor'                                            
      zero_dims = [0, 1, 2]                                                     
╰────────────────────────────────────────────────────────────────────────────────╯
╭─────────────────────── <class 'tidy3d.components.monitor.FieldMonitor'> ───────────────────────╮
 :class:`Monitor` that records electromagnetic fields in the frequency domain.                  
                                                                                                
 ╭────────────────────────────────────────────────────────────────────────────────────────────╮ 
  FieldMonitor(                                                                               
  type='FieldMonitor',                                                                    
  center=(0.0, 0.0, -1.0),                                                                
  size=(20.0, 20.0, 0.0),                                                                 
  name='field1',                                                                          
  freqs=(200000000000000.0,),                                                             
  apodization=ApodizationSpec(                                                            
  │   │   start=None,                                                                         
  │   │   end=None,                                                                           
  │   │   width=None,                                                                         
  │   │   type='ApodizationSpec'                                                              
  ),                                                                                      
  fields=('Ex', 'Ey', 'Ez', 'Hx', 'Hy', 'Hz'),                                            
  interval_space=(1, 1, 1),                                                               
  colocate=False                                                                          
  )                                                                                           
 ╰────────────────────────────────────────────────────────────────────────────────────────────╯ 
                                                                                                
    apodization = ApodizationSpec(start=None, end=None, width=None, type='ApodizationSpec')     
   bounding_box = Box(type='Box', center=(0.0, 0.0, -1.0), size=(20.0, 20.0, 0.0))              
         bounds = ((-10.0, -10.0, -1.0), (10.0, 10.0, -1.0))                                    
         center = (0.0, 0.0, -1.0)                                                              
       colocate = False                                                                         
         fields = ('Ex', 'Ey', 'Ez', 'Hx', 'Hy', 'Hz')                                          
          freqs = (200000000000000.0,)                                                          
       geometry = Box(type='Box', center=(0.0, 0.0, -1.0), size=(20.0, 20.0, 0.0))              
 interval_space = (1, 1, 1)                                                                     
           name = 'field1'                                                                      
    plot_params = PlotParams(                                                                   
                      alpha=0.4,                                                                
                      edgecolor='orange',                                                       
                      facecolor='orange',                                                       
                      fill=True,                                                                
                      hatch=None,                                                               
                      linewidth=3.0,                                                            
                      type='PlotParams'                                                         
                  )                                                                             
           size = (20.0, 20.0, 0.0)                                                             
           type = 'FieldMonitor'                                                                
      zero_dims = [2]                                                                           
╰────────────────────────────────────────────────────────────────────────────────────────────────╯
╭─────────────────────── <class 'tidy3d.components.monitor.FieldMonitor'> ───────────────────────╮
 :class:`Monitor` that records electromagnetic fields in the frequency domain.                  
                                                                                                
 ╭────────────────────────────────────────────────────────────────────────────────────────────╮ 
  FieldMonitor(                                                                               
  type='FieldMonitor',                                                                    
  center=(0.0, 0.0, 0.0),                                                                 
  size=(20.0, 0.0, 20.0),                                                                 
  name='field2',                                                                          
  freqs=(200000000000000.0,),                                                             
  apodization=ApodizationSpec(                                                            
  │   │   start=None,                                                                         
  │   │   end=None,                                                                           
  │   │   width=None,                                                                         
  │   │   type='ApodizationSpec'                                                              
  ),                                                                                      
  fields=('Ex', 'Ey', 'Ez', 'Hx', 'Hy', 'Hz'),                                            
  interval_space=(1, 1, 1),                                                               
  colocate=False                                                                          
  )                                                                                           
 ╰────────────────────────────────────────────────────────────────────────────────────────────╯ 
                                                                                                
    apodization = ApodizationSpec(start=None, end=None, width=None, type='ApodizationSpec')     
   bounding_box = Box(type='Box', center=(0.0, 0.0, 0.0), size=(20.0, 0.0, 20.0))               
         bounds = ((-10.0, 0.0, -10.0), (10.0, 0.0, 10.0))                                      
         center = (0.0, 0.0, 0.0)                                                               
       colocate = False                                                                         
         fields = ('Ex', 'Ey', 'Ez', 'Hx', 'Hy', 'Hz')                                          
          freqs = (200000000000000.0,)                                                          
       geometry = Box(type='Box', center=(0.0, 0.0, 0.0), size=(20.0, 0.0, 20.0))               
 interval_space = (1, 1, 1)                                                                     
           name = 'field2'                                                                      
    plot_params = PlotParams(                                                                   
                      alpha=0.4,                                                                
                      edgecolor='orange',                                                       
                      facecolor='orange',                                                       
                      fill=True,                                                                
                      hatch=None,                                                               
                      linewidth=3.0,                                                            
                      type='PlotParams'                                                         
                  )                                                                             
           size = (20.0, 0.0, 20.0)                                                             
           type = 'FieldMonitor'                                                                
      zero_dims = [1]                                                                           
╰────────────────────────────────────────────────────────────────────────────────────────────────╯

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.

First, let’s take a look at the source time dependence.

[9]:
# Visualize source
psource.source_time.plot(np.linspace(0, run_time, 1001))
plt.show()

../_images/notebooks_Simulation_17_0.png

And now let’s visualize the simulation.

For this, 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 grey.

[10]:
fig, ax = plt.subplots(1, 3, figsize=(13, 4))
sim.plot_eps(z=0.75, freq=freq0, ax=ax[0])
sim.plot_eps(y=0.01, freq=freq0, ax=ax[1])
sim.plot_eps(x=0, freq=freq0, ax=ax[2])

           INFO     Auto meshing using wavelength 1.4990 defined from sources.                     grid_spec.py:510
[10]:
<AxesSubplot: title={'center': 'cross section at x=0.00'}, xlabel='y', ylabel='z'>
../_images/notebooks_Simulation_19_2.png

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

[11]:
fig, ax = plt.subplots(1, 3, figsize=(12, 3))
sim.plot(z=0.75, ax=ax[0])
sim.plot(y=0.01, ax=ax[1])
sim.plot(x=0, ax=ax[2])

[11]:
<AxesSubplot: title={'center': 'cross section at x=0.00'}, xlabel='y', ylabel='z'>
../_images/notebooks_Simulation_21_1.png

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.

[12]:
task_id = web.upload(sim, task_name="Simulation")

[09:48:12] INFO     Using Tidy3D credentials from stored file.                                           auth.py:70
[09:48:16] INFO     Authentication successful.                                                           auth.py:30
[09:48:17] INFO     Created task 'Simulation' with task_id 'd88ab538-1d3f-4bcc-a447-9bfc97de43fc'.    webapi.py:120

The task is still in draft status and will not run until we call the start function. Before that, we may want to check the estimated cost of the task.

[13]:
print("Max flex unit cost: ", web.estimate_cost(task_id))

Max flex unit cost:  0.025

We can now start the task, and if we want to, continously monitor its status, and wait until the run is successful. The monitor function will keep running until either a 'success' or 'error' status is returned.

[14]:
web.start(task_id)
web.monitor(task_id)

[09:48:23] INFO     Maximum FlexUnit cost: 0.025                                                      webapi.py:252
           INFO     status = queued                                                                   webapi.py:261
[09:48:26] INFO     status = preprocess                                                               webapi.py:273
[09:48:33] INFO     starting up solver                                                                webapi.py:277
[09:48:43] INFO     running solver                                                                    webapi.py:283
[09:48:47] INFO     early shutoff detected, exiting.                                                  webapi.py:294
           INFO     status = postprocess                                                              webapi.py:300
[09:48:49] INFO     status = success                                                                  webapi.py:306
           INFO     Billed FlexUnit cost: 0.025                                                       webapi.py:310

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.to_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 data/ folder.

[15]:
sim_data = web.load(task_id, path="data/sim_data.hdf5")

# Show the output of the log file
print(sim_data.log)

           INFO     downloading file "output/monitor_data.hdf5" to "data/sim_data.hdf5"               webapi.py:592
[09:48:51] INFO     loading SimulationData from data/sim_data.hdf5                                    webapi.py:414
Simulation domain Nx, Ny, Nz: [156, 157, 104]
Applied symmetries: (0, 0, 0)
Number of computational grid points: 2.6629e+06.
Using subpixel averaging: True
Number of time steps: 3.4930e+03
Automatic shutoff factor: 1.00e-05
Time step (s): 5.7275e-17


Compute source modes time (s):     0.2348
Compute monitor modes time (s):    0.0178
Rest of setup time (s):            5.1093

Running solver for 3493 time steps...
- Time step    139 / time 7.96e-15s (  4 % done), field decay: 1.00e+00
- Time step    279 / time 1.60e-14s (  8 % done), field decay: 1.00e+00
- Time step    419 / time 2.40e-14s ( 12 % done), field decay: 1.00e+00
- Time step    558 / time 3.20e-14s ( 16 % done), field decay: 1.00e+00
- Time step    698 / time 4.00e-14s ( 20 % done), field decay: 1.00e+00
- Time step    838 / time 4.80e-14s ( 24 % done), field decay: 1.00e+00
- Time step    978 / time 5.60e-14s ( 28 % done), field decay: 1.00e+00
- Time step   1117 / time 6.40e-14s ( 32 % done), field decay: 1.00e+00
- Time step   1257 / time 7.20e-14s ( 36 % done), field decay: 1.00e+00
- Time step   1389 / time 7.96e-14s ( 39 % done), field decay: 1.00e+00
- Time step   1397 / time 8.00e-14s ( 40 % done), field decay: 1.00e+00
- Time step   1536 / time 8.80e-14s ( 44 % done), field decay: 1.00e+00
- Time step   1676 / time 9.60e-14s ( 48 % done), field decay: 8.05e-01
- Time step   1816 / time 1.04e-13s ( 52 % done), field decay: 3.46e-01
- Time step   1956 / time 1.12e-13s ( 56 % done), field decay: 1.51e-01
- Time step   2095 / time 1.20e-13s ( 60 % done), field decay: 5.29e-02
- Time step   2235 / time 1.28e-13s ( 64 % done), field decay: 1.21e-02
- Time step   2375 / time 1.36e-13s ( 68 % done), field decay: 2.31e-03
- Time step   2514 / time 1.44e-13s ( 72 % done), field decay: 4.89e-04
- Time step   2654 / time 1.52e-13s ( 76 % done), field decay: 1.11e-04
- Time step   2794 / time 1.60e-13s ( 80 % done), field decay: 2.39e-05
- Time step   2934 / time 1.68e-13s ( 84 % done), field decay: 4.75e-06
Field decay smaller than shutoff factor, exiting solver.

Solver time (s):                   1.1452

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

[16]:
fig, ax = plt.subplots(1, 2, figsize=(10, 4))
sim_data.plot_field("field1", "Ey", z=-1.0, ax=ax[0], val="real")
sim_data.plot_field("field2", "Ey", ax=ax[1], val="real")

           INFO     Auto meshing using wavelength 1.4990 defined from sources.                     grid_spec.py:510
[16]:
<AxesSubplot: title={'center': 'cross section at y=0.00'}, xlabel='x', ylabel='z'>
../_images/notebooks_Simulation_32_2.png

Monitor data#

The raw field data can be accessed through indexing by monitor name directly.

For plenty of discussion on accessing and manipulating data, refer to the data visualization tutorial.

[17]:
mon1_data = sim_data["field1"]
mon1_data.Ex

[17]:
<xarray.ScalarFieldDataArray (x: 158, y: 159, z: 1, f: 1)>
array([[[[ 0.0000000e+00+0.0000000e+00j]],

        [[ 0.0000000e+00+0.0000000e+00j]],

        [[ 9.8186092e-07+4.0431851e-07j]],

        ...,

        [[ 9.9795454e-07+1.7111465e-07j]],

        [[ 1.2159646e-07+6.0561064e-08j]],

        [[ 0.0000000e+00+0.0000000e+00j]]],


       [[[ 0.0000000e+00+0.0000000e+00j]],

        [[ 0.0000000e+00+0.0000000e+00j]],

        [[ 9.8186092e-07+4.0431851e-07j]],
...
        [[-1.1717159e-07+2.6605029e-08j]],

        [[ 8.6237844e-09-7.4280364e-09j]],

        [[ 0.0000000e+00+0.0000000e+00j]]],


       [[[ 0.0000000e+00+0.0000000e+00j]],

        [[ 0.0000000e+00+0.0000000e+00j]],

        [[-8.3378872e-08+3.6394130e-08j]],

        ...,

        [[-1.1717159e-07+2.6605029e-08j]],

        [[ 8.6237844e-09-7.4280364e-09j]],

        [[ 0.0000000e+00+0.0000000e+00j]]]], dtype=complex64)
Coordinates:
  * f        (f) float64 2e+14
  * x        (x) float64 -2.379 -2.348 -2.318 -2.288 ... 2.288 2.318 2.348 2.379
  * y        (y) float64 -2.39 -2.36 -2.33 -2.3 ... 2.266 2.295 2.325 2.354
  * z        (z) float64 -1.0
Attributes:
    long_name:  field value
[18]:
ax = mon1_data.Ez.real.plot()

../_images/notebooks_Simulation_35_0.png

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

[19]:
time_data = sim_data["field_time"]
fig, ax = plt.subplots(1)
time_data.Ey.plot()
ax.set_ylabel("$E_y(t)$ [V/m]")
plt.show()

../_images/notebooks_Simulation_37_0.png

Permittivity data#

We can also query the relative permittivity in the simulation within a volume parameterized by a td.Box. The method Simulation.epsilon(box, coord_key) returns the permittivity within the specified volume.

The coord_key specifies at what locations in the yee cell to evaluate the permittivity at (eg. 'centers', 'Ey', 'Hz', etc.).

[20]:
volume = td.Box(center=(0, 0, 0.75), size=(5, 5, 0))

# at Yee cell centers
eps_centers = sim.epsilon(box=volume, coord_key="centers")

# at Ex locations in the yee cell
eps_Ex = sim.epsilon(box=volume, coord_key="Ex")

Return an xarray DataArray containing the complex-valued permittivity values at the Yee cell centers and the “Ex” within the box.

We can then plot or post-process this data as we wish.

[21]:
f, (ax1, ax2) = plt.subplots(1, 2, tight_layout=True, figsize=(10, 4))

eps_centers.real.plot(x="x", y="y", cmap="Greys", ax=ax1)
eps_Ex.real.plot(x="x", y="y", cmap="Greys", ax=ax2)
ax1.set_title("epsilon_r at centers")
ax2.set_title("epsilon_r at Ex locations")

plt.show()

../_images/notebooks_Simulation_41_0.png
[ ]: