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

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

# 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=(.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'                                                                             
╰─────────────────────────────────────────────────────────────────────────────────────────────────────────────────╯
╭─────────────────────────────── <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,), fields=('Ex', 'Ey', 'Ez', 'Hx', 'Hy', 'Hz'), interval_space=(1, 1, 1),           
  colocate=False)                                                                                              
 ╰─────────────────────────────────────────────────────────────────────────────────────────────────────────────╯ 
                                                                                                                 
   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'                                                                                 
╰─────────────────────────────────────────────────────────────────────────────────────────────────────────────────╯
╭─────────────────────────────── <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,), 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=(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'                                                                                 
╰─────────────────────────────────────────────────────────────────────────────────────────────────────────────────╯

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]);
[09:04:50] INFO     Auto meshing using wavelength 1.4990 defined from sources.                     grid_spec.py:473
../_images/notebooks_Simulation_19_1.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]);
../_images/notebooks_Simulation_21_0.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')
 simulation.json ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 0.0%0.0/9.2 kB?-:--:--

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

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)
🚶  Finishing 'Simulation'...

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)
 monitor_data.hdf5 ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 97.6%4.2/4.3 MB562.7 kB/s0:00:01
[09:05:44] INFO     loading SimulationData from data/sim_data.hdf5                                    webapi.py:397
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.8430e+03
Automatic shutoff factor: 1.00e-05
Time step (s): 5.2068e-17


Compute source modes time (s):     0.0779
Compute monitor modes time (s):    0.0015
Rest of setup time (s):            3.4994

Running solver for 3843 time steps...
- Time step    153 / time 7.97e-15s (  4 % done), field decay: 1.00e+00
- Time step    307 / time 1.60e-14s (  8 % done), field decay: 1.00e+00
- Time step    461 / time 2.40e-14s ( 12 % done), field decay: 1.00e+00
- Time step    614 / time 3.20e-14s ( 16 % done), field decay: 1.00e+00
- Time step    768 / time 4.00e-14s ( 20 % done), field decay: 1.00e+00
- Time step    922 / time 4.80e-14s ( 24 % done), field decay: 1.00e+00
- Time step   1076 / time 5.60e-14s ( 28 % done), field decay: 1.00e+00
- Time step   1229 / time 6.40e-14s ( 32 % done), field decay: 1.00e+00
- Time step   1383 / time 7.20e-14s ( 36 % done), field decay: 1.00e+00
- Time step   1528 / time 7.96e-14s ( 39 % done), field decay: 1.00e+00
- Time step   1537 / time 8.00e-14s ( 40 % done), field decay: 1.00e+00
- Time step   1690 / time 8.80e-14s ( 44 % done), field decay: 1.00e+00
- Time step   1844 / time 9.60e-14s ( 48 % done), field decay: 7.18e-01
- Time step   1998 / time 1.04e-13s ( 52 % done), field decay: 2.58e-01
- Time step   2152 / time 1.12e-13s ( 56 % done), field decay: 1.12e-01
- Time step   2305 / time 1.20e-13s ( 60 % done), field decay: 4.16e-02
- Time step   2459 / time 1.28e-13s ( 64 % done), field decay: 9.96e-03
- Time step   2613 / time 1.36e-13s ( 68 % done), field decay: 2.08e-03
- Time step   2766 / time 1.44e-13s ( 72 % done), field decay: 4.91e-04
- Time step   2920 / time 1.52e-13s ( 76 % done), field decay: 1.17e-04
- Time step   3074 / time 1.60e-13s ( 80 % done), field decay: 2.51e-05
- Time step   3228 / time 1.68e-13s ( 84 % done), field decay: 4.92e-06
Field decay smaller than shutoff factor, exiting solver.

Solver time (s):                   1.2578
Post-processing time (s):          0.2348


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

[17]:
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')
../_images/notebooks_Simulation_32_0.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.

[18]:
mon1_data = sim_data['field1']
mon1_data.Ex
[18]:
<xarray.ScalarFieldDataArray (x: 158, y: 159, z: 1, f: 1)>
array([[[[ 0.00000000e+00+0.00000000e+00j]],

        [[ 0.00000000e+00+0.00000000e+00j]],

        [[ 4.08210345e-07+7.69336070e-08j]],

        ...,

        [[ 3.08751853e-07-7.09521903e-08j]],

        [[ 6.91019995e-09+4.49927794e-08j]],

        [[ 0.00000000e+00+0.00000000e+00j]]],


       [[[ 0.00000000e+00+0.00000000e+00j]],

        [[ 0.00000000e+00+0.00000000e+00j]],

        [[ 4.08210345e-07+7.69336070e-08j]],
...
        [[-3.39188149e-08+1.06002510e-07j]],

        [[-2.24265685e-08-3.78733208e-08j]],

        [[ 0.00000000e+00+0.00000000e+00j]]],


       [[[ 0.00000000e+00+0.00000000e+00j]],

        [[ 0.00000000e+00+0.00000000e+00j]],

        [[ 1.51480892e-08+1.63971981e-08j]],

        ...,

        [[-3.39188149e-08+1.06002510e-07j]],

        [[-2.24265685e-08-3.78733208e-08j]],

        [[ 0.00000000e+00+0.00000000e+00j]]]])
Coordinates:
  * 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
  * f        (f) float64 2e+14
Attributes:
    long_name:  field value
[21]:
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.

[24]:
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.).

[25]:
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.

[26]:
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
[ ]: