Creating FDTD animations

Animations are often created from FDTD simulations to provide a more intuitive understanding of the physical phenomena being modeled. These animations can visualize the evolution of the field distribution over time, showing wave propagation, interactions, and other dynamic effects that static images cannot adequately depict. This can help researchers and engineers better understand the results of their simulations, identify anomalies, and develop or refine their models. Additionally, these animations can serve as effective educational and communication tools, making complex concepts more accessible to non-experts.

This tutorial demonstrates how to create animations in a waveguide taper simulation and display them directly in a Jupyter Notebook. We will show how to use Tidy3D convenience features to animate time-domain and frequency-domain fields.

We also provide a conprehensive list of other tutorials such as how to define boundary conditions, how to defining spatially-varying sources and structures, how to model dispersive materials, and how to visualize simulation setups and results.

If you are new to the finite-difference time-domain (FDTD) method, we highly recommend going through our FDTD101 tutorials.

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from IPython.display import HTML

import tidy3d as td
import tidy3d.web as web

Simulation Setup

For this demo, we simply simulate a linear waveguide taper. The simulation setup is directly borrowed from the waveguide size converter case study notebook.

lda0 = 1.55  # central wavelength
freq0 = td.C_0 / lda0  # central frequency
ldas = np.linspace(1.5, 1.6, 101)  # wavelength range
freqs = td.C_0 / ldas  # frequency range
fwidth = 0.5 * (np.max(freqs) - np.min(freqs))  # width of the frequency distribution

n_si = 3.48  # silicon refractive index
si = td.Medium(permittivity=n_si**2)

n_sio2 = 1.44  # silicon oxide refractive index
sio2 = td.Medium(permittivity=n_sio2**2)

w_in = 5  # input waveguide width
w_out = 0.45  # output waveguide width
t_wg = 0.11  # waveguide thickness
L_t = 10
inf_eff = 1e3  # effective infinity of the model

# define the substrate structure
sub = td.Structure(
    geometry=td.Box.from_bounds(
        rmin=(-inf_eff, -inf_eff, -inf_eff), rmax=(inf_eff, inf_eff, 0)
    ),
    medium=sio2,
)

# vertices of the taper
vertices = [
    [-inf_eff, w_in / 2],
    [0, w_in / 2],
    [L_t, w_out / 2],
    [inf_eff, w_out / 2],
    [inf_eff, -w_out / 2],
    [L_t, -w_out / 2],
    [0, -w_in / 2],
    [-inf_eff, -w_in / 2],
]

# construct the taper structure using a PolySlab
linear_taper = td.Structure(
    geometry=td.PolySlab(vertices=vertices, axis=2, slab_bounds=(0, t_wg)),
    medium=si,
)

To create a time-domain animation, we need to capture the frames at different time instances of the simulation. This can be done by using a FieldTimeMonitor. Usually, an FDTD simulation contains a large number of time steps and grid points. Recording the field at every time step and grid point will result in a large dataset. For the purpose of making time-domain animations, this is usually unnecessary. In this particular case, we will record one frame every 100 time steps, which is realized by setting interval=100 in the FieldTimeMonitor. Similarly, we set interval_space=(3,3,1) such that the field is recorded every 3 grid points. To further reduce dataset size, we will only be recording the z-component of the magnetic field.

run_time = 5e-13  # run time of the simulation

# define a field time monitor to record time-domain animation frames
time_monitor = td.FieldTimeMonitor(
    center=(0, 0, t_wg/2),
    size=(td.inf, td.inf, 0),
    fields=['Hz'],
    start=0,
    stop=run_time,
    interval=100,
    interval_space=(3,3,1),
    name='time_movie_monitor'
)

To create a frequency-domain animation, we will include a FieldMonitor in the simulation to obtain the steady-state electromagnetic fields at the end of the simulation. Afterward, we will change the field phase using a convenience function to create the frequency-domain animation.

frequency_monitor = td.FieldMonitor(
    center=(0, 0, t_wg/2),
    size=(td.inf, td.inf, 0),
    freqs=[freq0],
    fields=['Hz'],
    name='frequency_monitor'
)

Now, we can define the simulation source and create the FDTD simulation object.

# define a mode source
mode_source = td.ModeSource(
    center=(-lda0 / 2, 0, t_wg / 2),
    size=(0, 1.2 * w_in, 6 * t_wg),
    source_time=td.GaussianPulse(freq0=freq0, fwidth=fwidth),
    direction="+",
    mode_spec=td.ModeSpec(num_modes=1, target_neff=n_si),
    mode_index=0,
)

# define simulation domain size
Lx = L_t + 2 * lda0
Ly = w_in + 2 * lda0
Lz = t_wg + 1.5 * lda0
sim_size = (Lx, Ly, Lz)

# define simulation
sim = td.Simulation(
    center=(L_t / 2, 0, t_wg),
    size=sim_size,
    grid_spec=td.GridSpec.auto(min_steps_per_wvl=15, wavelength=lda0),
    structures=[linear_taper, sub],
    sources=[mode_source],
    monitors=[time_monitor, frequency_monitor],
    run_time=run_time,
    boundary_spec=td.BoundarySpec.all_sides(
        boundary=td.PML()
    ),  # pml is used in all boundaries
    symmetry=(0, -1, 0),
)  # a pec symmetry plane at y=0 can be used to reduce the grid points of the simulation

Visualize the simulation.

sim.plot(z=t_wg/2)
plt.show()

Submit the simulation.

sim_data = web.run(simulation=sim, task_name="taper", path="data/simulation_data.hdf5")

15:08:34 UTC Created task 'taper' with task_id
             'fdve-d41972d9-2882-4449-a7e2-0017c0637842' and task_type 'FDTD'.

             View task using web UI at
             'https://tidy3d.simulation.cloud/workbench?taskId=fdve-d41972d9-288
             2-4449-a7e2-0017c0637842'.

             status = success

15:08:35 UTC loading simulation from data/simulation_data.hdf5

Creating a Time-Domain Animation

Before creating the time-domain animation, we can inspect the field distribution at a few time instances.

fig, ax = plt.subplots(1,2, figsize=(10,5), tight_layout=True)

sim_data['time_movie_monitor'].Hz.sel(t=1e-13, method='nearest').plot(x='x', y='y', ax=ax[0])

sim_data['time_movie_monitor'].Hz.sel(t=2e-13, method='nearest').plot(x='x', y='y', ax=ax[1])
plt.show()

To create an animation, we need to define a function that plots each frame. Then we use matplotlib’s FuncAnimation to create the animation. For this particular animation, we use 50 frames. For detail on FuncAnimation, please reference the documentation.

After the animation is created, we can display it directly in the notebook environment by using HTML. Here we see the incident pulse from the wider waveguide is squeezed into the narrower waveguide through the taper. Some field is reflected back or leaked out of the waveguide, leading to some insertion loss, which can be minimized by using a longer taper.

t_end = sim_data['time_movie_monitor'].Hz.coords['t'][-1]  # end time of the animation
frames = 50  # number of frames

fig, ax = plt.subplots()

def animate(i):
    t = t_end*i/frames # time at each frame
    sim_data['time_movie_monitor'].Hz.sel(t=t, method='nearest').plot(x='x', y='y', ax=ax, vmin=-0.1, vmax=0.1, add_colorbar=False, cmap='seismic')

# create animation
ani = animation.FuncAnimation(fig, animate, frames=frames);
plt.close()

# display the animation
HTML(ani.to_jshtml())

Once Loop Reflect

Creating a Frequency-Domain Animation

Using the phase parameter, the plot_field function allows us to apply a phase shift to frequency-domain monitors. Below, we will apply a \(\pi\) phase shift to the Hz taper field component.

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12,5), tight_layout=True)

sim_data.plot_field("frequency_monitor", "Hz", "real", f=freq0, ax=ax1)
ax1.set_xlim(6, 11)
ax1.set_ylim(-2.5, 2.5)

sim_data.plot_field("frequency_monitor", "Hz", "real", f=freq0, ax=ax2, phase=np.pi)
ax2.set_xlim(6, 11)
ax2.set_ylim(-2.5, 2.5)

plt.show()

The apply_phase() is another convenient function for modifying the Fourier-transformed field phase. We will use it to change the field phase between \(0\) and \(-2 \pi\).

import matplotlib.colors as mcolors
from mpl_toolkits.axes_grid1 import make_axes_locatable

frames = 36
phi_range = np.linspace(0, -2*np.pi, frames)
vmax = np.amax(np.abs(sim_data['frequency_monitor'].Hz.values))*0.7
fig, ax1 = plt.subplots(1, 1, figsize=(6, 4), tight_layout=True)

def animate_ii(i):
    sim_data_p = sim_data['frequency_monitor']
    sim_data_p = sim_data_p.apply_phase(phi_range[i])
    Hz = sim_data_p.Hz
    x = Hz["x"]
    y = Hz["y"]
    Hz = np.squeeze(np.real(np.flipud(np.rot90(Hz[:,:,0]))))
    divnorm = mcolors.TwoSlopeNorm(vmin=-vmax, vcenter=0, vmax=vmax)
    image = ax1.pcolormesh(x, y, Hz, cmap='bwr', shading='gouraud', norm=divnorm)
    sim_data.simulation.plot_structures_eps(z=0, alpha=0.25, cbar=False,
                            hlim=[np.amin(x), np.amax(x)], vlim=[np.amin(y), np.amax(y)], ax=ax1)

    divider = make_axes_locatable(ax1)
    cax = divider.new_horizontal(size = "2%", pad = 0.05)
    fig.add_axes(cax)
    fig.colorbar(image, cax = cax, orientation = "vertical", label='$H_z$')

    ax1.set_xlim([np.amin(x), np.amax(x)])
    ax1.set_ylim([np.amin(y), np.amax(y)])
    ax1.set_xlabel("x ($\\mu$m)")
    ax1.set_ylabel("y ($\\mu$m)")
    ax1.set_aspect("equal")

# create animation
ani = animation.FuncAnimation(fig, animate_ii, frames=frames);
plt.close()

# display the animation
HTML(ani.to_jshtml())

Once Loop Reflect