Parameter scan

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

diagram

In this notebook, we will show an example of using tidy3d to evaluate device performance over a set of many design parameters.

This example will also provide a walkthrough of Tidy3D’s Job and Batch features for managing both individual simulations and sets of simulations.

For demonstration, we look at the splitting ratio of a directional coupler as we vary the coupling length between two waveguides.

[1]:
# get the most recent version of tidy3d
!pip install -q --upgrade tidy3d

# gdspy is also needed for gds file manipulation
!pip install -q gdspy

# make sure notebook plots inline
%matplotlib inline

# standard python imports
import numpy as np
import matplotlib.pyplot as plt
import os
import tqdm
import gdspy

# tidy3D imports
import tidy3d as td
from tidy3d import web

# set tidy3d to only print warning information to reduce verbosity
td.logging_level('warning')

Setup

First we set up some global parameters

[2]:
# wavelength / frequency
lambda0 = 1.550                     # all length scales in microns
freq0 = td.constants.C_0 / lambda0
fwidth = freq0 / 10

# Permittivity of waveguide and substrate
wg_n = 3.48
sub_n = 1.45
mat_wg = td.Medium(n=wg_n)
mat_sub = td.Medium(n=sub_n)

# Waveguide dimensions

# Waveguide height
wg_height = 0.22
# Waveguide width
wg_width = 0.45
# Waveguide separation in the beginning/end
wg_spacing_in = 8
# Total device length along propagation direction
device_length = 100
# Length of the bend region
bend_length = 16
# space between waveguide and PML
pml_spacing = 1
# Mesh step in all directions
mesh_step = 0.040

Define waveguide bends and coupler

Here is where we define our directional coupler shape programmatically in terms of the geometric parameters

[3]:
def bend_pts(bend_length, width, npts=10):
    """ Set of points describing a tanh bend from (0, 0) to (length, width)"""
    x = np.linspace(0, bend_length, npts)
    y = width*(1 + np.tanh(6*(x/bend_length - 0.5)))/2
    return np.stack((x, y), axis=1)

def arm_pts(length, width, coup_length, bend_length, npts_bend=30):
    """ Set of points defining one arm of an integrated coupler """
    ### Make the right half of the coupler arm first
    # Make bend and offset by coup_length/2
    bend = bend_pts(bend_length, width, npts_bend)
    bend[:, 0] += coup_length / 2
    # Add starting point as (0, 0)
    right_half = np.concatenate(([[0, 0]], bend))
    # Add an extra point to make sure waveguide is straight past the bend
    right_half = np.concatenate((right_half, [[right_half[-1, 0] + 0.1, width]]))
    # Add end point as (length/2, width)
    right_half = np.concatenate((right_half, [[length/2, width]]))

    # Make the left half by reflecting and omitting the (0, 0) point
    left_half = np.copy(right_half)[1:, :]
    left_half[:, 0] = -left_half[::-1, 0]
    left_half[:, 1] = left_half[::-1, 1]

    return np.concatenate((left_half, right_half), axis=0)

def make_coupler(
    length,
    wg_spacing_in,
    wg_width,
    wg_spacing_coup,
    coup_length,
    bend_length,
    npts_bend=30):
    """ Make an integrated coupler using the gdspy FlexPath object. """

    # Compute one arm of the coupler
    arm_width = (wg_spacing_in - wg_width - wg_spacing_coup)/2
    arm = arm_pts(length, arm_width, coup_length, bend_length, npts_bend)
    # Reflect and offset bottom arm
    coup_bot = np.copy(arm)
    coup_bot[:, 1] = -coup_bot[::-1, 1] - wg_width/2 - wg_spacing_coup/2
    # Offset top arm
    coup_top = np.copy(arm)
    coup_top[:, 1] += wg_width/2 + wg_spacing_coup/2

    # Create waveguides as GDS paths
    path_bot = gdspy.FlexPath(coup_bot, wg_width, layer=1, datatype=0)
    path_top = gdspy.FlexPath(coup_top, wg_width, layer=1, datatype=1)

    return [path_bot, path_top]

Create Simulation and Submit Job

The following function creates a tidy3d simulation object for a set of design parameters.

Note that the simulation has not been run yet, just created.

[4]:
def make_sim(coup_length, wg_spacing_coup):
    """ gets the parameters from the scan,
        creates a simulation,
        exports it to the server to run,
        returns a taskID handle to use to get the results later
    """

    gdspy.current_library = gdspy.GdsLibrary()
    lib = gdspy.GdsLibrary()

    # Geometry must be placed in GDS cells to import into Tidy3D
    coup_cell = lib.new_cell('Coupler')

    substrate = gdspy.Rectangle(
        (-device_length/2, -wg_spacing_in/2-10),
        (device_length/2, wg_spacing_in/2+10),
        layer=0)
    coup_cell.add(substrate)

    # Add the coupler to a gdspy cell
    gds_coup = make_coupler(
        device_length,
        wg_spacing_in,
        wg_width,
        wg_spacing_coup,
        coup_length,
        bend_length)
    coup_cell.add(gds_coup)

    # Substrate
    oxide = td.GdsSlab(
        material=mat_sub,
        gds_cell=coup_cell,
        gds_layer=0,
        z_min=-10,
        z_max=0)

    # Waveguides (import all datatypes if gds_dtype not specified)
    coupler = td.GdsSlab(
        material=mat_wg,
        gds_cell=coup_cell,
        gds_layer=1,
        z_cent=wg_height/2,
        z_size=wg_height,)

    # Simulation size along propagation direction
    sim_length = 2 + 2*bend_length + coup_length

    # Spacing between waveguides and PML
    sim_size = [
        sim_length,
        wg_spacing_in + wg_width + 2*pml_spacing,
        wg_height + 2*pml_spacing]


    # source
    src_pos = -sim_length/2 + 0.5
    msource = td.ModeSource(
        center=[src_pos , wg_spacing_in / 2 , wg_height / 2],
        size=[0, 3, 2],
        source_time = td.GaussianPulse(
            frequency=freq0,
            fwidth=fwidth),
        direction='forward')

    domain_monitor = td.FreqMonitor(
        center = [0,0,wg_height/2],
        size = [sim_size[0], sim_size[1], 0],
        freqs = [freq0]
    )

    mon_in = td.ModeMonitor(
        center=[(src_pos + 0.5), wg_spacing_in / 2, wg_height / 2],
        size=[0, 3, 2],
        freqs=[freq0],
        Nmodes=1,
        name='in')
    mon_ref = td.ModeMonitor(
        center=[(src_pos + 0.5), -wg_spacing_in / 2, wg_height / 2],
        size=[0, 3, 2],
        freqs=[freq0],
        Nmodes=1,
        name='refect')
    mon_top = td.ModeMonitor(
        center=[-(src_pos + 0.5), wg_spacing_in / 2, wg_height / 2],
        size=[0, 3, 2],
        freqs=[freq0],
        Nmodes=1,
        name='top')
    mon_bot = td.ModeMonitor(
        center=[-(src_pos + 0.5), -wg_spacing_in / 2, wg_height / 2],
        size=[0, 3, 2],
        freqs=[freq0],
        Nmodes=1,
        name='bottom')
    wg_monitors = [mon_in, mon_ref, mon_top, mon_bot]

    # initialize the simulation
    sim = td.Simulation(
        size=sim_size,
        mesh_step=mesh_step,
        structures=[oxide, coupler],
        sources=[msource],
        monitors=[domain_monitor] + wg_monitors,
        run_time=20/fwidth,
        pml_layers=[12, 12, 12])

    # set the modes
    sim.compute_modes(msource, Nmodes=2)
    sim.set_mode(msource, mode_ind=0)

    return sim

Inspect Simulation

Let’s create and inspect a single simulation to make sure it was defined correctly before doing the full scan.

[5]:
# Length of the coupling region
coup_length = 10

# Waveguide separation in the coupling region
wg_spacing_coup = 0.10

sim = make_sim(coup_length, wg_spacing_coup)
[6]:
# visualize geometry
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 4))
sim.viz_mat_2D(normal='z', position=wg_height/2, ax=ax1);
sim.viz_mat_2D(normal='x', ax=ax2, source_alpha=1);
ax2.set_xlim([-3, 3])
plt.show()
../_images/examples_ParameterScan_10_0.png

Create and Submit Job

The Job object provides an interface for managing simulations.

job = Job(simulation) will create a job and upload the simulation to our server to run.

Then, one may call various methods of job to monitor progress, download results, and get information.

For more information, refer to the API reference.

[7]:
# create job, upload sim to server to begin running
job = web.Job(sim, task_name='CouplerVerify')

# monitor progress of the job
job.monitor()

# download the results and load them into a simulation
sim = job.load_results()
Project 'CouplerVerify-16b90ecdcbef1d80' status: success...

Postprocessing

The following function takes a completed simulation (with data loaded into it) and computes the quantities of interest.

For this case, we measure both the total transmission in the right ports and also the ratio of power between the top and bottom ports.

[8]:
def measure_transmission(sim, verbose=False):

    # download results if the job has finished

    if verbose:
        with open("output/tidy3d.log") as f:
             print(f.read())

    def get_power(monitor):
        # gets the mode amplitudes and returns the total forward and backward power in monitor.
        f, b = sim.data(monitor)['mode_amps']
        F, B = np.abs(f)**2, np.abs(b)**2
        return F, B

    # get the momitors and the normalization (input) power
    _, incident, reflect, top, bot = sim.monitors
    norm, _ = get_power(incident)
    norm = np.squeeze(norm)

    # computes S matrix
    S = np.zeros((2, 2))
    S[0,0] = get_power(incident)[-1]
    S[1,0] = get_power(reflect)[-1]
    S[0,1] = get_power(top)[0]
    S[1,1] = get_power(bot)[0]
    S = S/norm

    # computes quantities of interest
    split_ratio = S[0, 1] / (S[0, 1] + S[1, 1])
    efficiency = (S[0, 1] + S[1, 1]) / norm

    if verbose:
        print(f'split ratio of {(split_ratio * 100):.2f}% to top port')
        print(f'efficiency of {(efficiency * 100):.2f}% to transmission port')

    return split_ratio, efficiency
[9]:
# monitor and test out the measure_transmission function the results of the single run
split_ratio, efficiency = measure_transmission(sim, verbose=True)
Simulation domain Nx, Ny, Nz: [190, 158, 94]
Applied symmetries: [0, 0, 0]
Number of computational grid points: 2.8219e+06.
Using subpixel averaging: True
Number of time steps: 19092
Automatic shutoff factor: 1.00e-05
Time step (s): 5.4164e-17

Compute source modes time (s):     2.1846
Compute monitor modes time (s):    6.8366

Rest of setup time (s):            0.0635

Starting solver...
- Time step    760 / time 4.12e-14s (  3 % done), field decay: 1.00e+00
- Time step    763 / time 4.13e-14s (  4 % done), field decay: 1.00e+00
- Time step   1527 / time 8.27e-14s (  8 % done), field decay: 1.00e+00
- Time step   2291 / time 1.24e-13s ( 12 % done), field decay: 9.37e-03
- Time step   3054 / time 1.65e-13s ( 16 % done), field decay: 2.39e-05
- Time step   3818 / time 2.07e-13s ( 20 % done), field decay: 7.85e-07
Field decay smaller than shutoff factor, exiting solver.

Solver time (s):                   4.3262
Post-processing time (s):          0.1982

split ratio of 26.02% to top port
efficiency of 99.26% to transmission port
[10]:
fig, ax = plt.subplots(1, 1, figsize=(16, 3))
sim.viz_field_2D(sim.monitors[0], ax=ax, cbar=True, comp='y', val='re', eps_alpha=0.2)
plt.show()
../_images/examples_ParameterScan_16_0.png
[11]:
fig, axes = plt.subplots(1, 4, figsize=(16, 3))
for m, ax in zip(sim.monitors[1:], axes):
    im = sim.viz_field_2D(m, ax=ax, cbar=True, comp='y', val='abs', eps_alpha=0.2);
../_images/examples_ParameterScan_17_0.png

1D Parameter Scan

Now we will scan through the coupling length parameter to see the effect on splitting ratio.

To do this, we will create a list of simulations corresponding to each parameter combination.

We will use this list to create a Batch object, which has similar functionality to Job but allows one to manage a set of jobs.

First, we create arrays to store the input and output values.

[12]:
# create variables to store parameters, simulation information, results
Nl = 11

ls = np.linspace(5, 12, Nl)
split_ratios = np.zeros(Nl)
efficiencies = np.zeros(Nl)

Create Batch

We now create our list of simulations and use them to initialize a Batch.

For more information, refer to the API reference.

[13]:
# submit all jobs
sims = [make_sim(l, wg_spacing_coup) for l in ls]
batch = web.Batch(sims)

Monitor Batch

Here we can perform real-time monitoring of how many of the jobs in the batch have completed.

[14]:
batch.monitor()
Percentage of jobs completed: : 100%|██████████| 11/11 [07:41<00:00, 41.93s/it]

Download and Compute Results

When the jobs are all done, we can download results to create another list of simulations where the data has been loaded into them.

[15]:
# get results from all jobs
sims_loaded = batch.load_results()

loading results for job (1/11)


loading results for job (2/11)


loading results for job (3/11)


loading results for job (4/11)


loading results for job (5/11)


loading results for job (6/11)


loading results for job (7/11)


loading results for job (8/11)


loading results for job (9/11)


loading results for job (10/11)


loading results for job (11/11)

Plot Reults

Finally, we can compute the output quantities and load them into the arrays we created initally.

Then we may plot the results.

[16]:
for i, sim in enumerate(sims_loaded):
    split_ratio, efficiency = measure_transmission(sim)
    split_ratios[i] = split_ratio
    efficiencies[i] = efficiency
[17]:
plt.plot(ls, 100*split_ratios, 'k')
plt.plot(ls, 50 * np.ones_like(ls), 'gray', linestyle='--')
plt.xlabel('coupling length (um)')
plt.ylabel('splitting ratio (%)')
plt.ylim(0, 100)
plt.show()
../_images/examples_ParameterScan_28_0.png
[18]:
plt.plot(ls, efficiencies)
plt.xlabel('coupling length (um)')
plt.ylabel('$P_{out} / P_{in}$')
plt.ylim((.9, 1))
plt.title('transmission into right ports')
plt.show()
../_images/examples_ParameterScan_29_0.png
[19]:
fig, axes = plt.subplots(Nl, 1, figsize=(6, 30))
for (sim, ax) in zip(sims_loaded, axes):
    sim.viz_field_2D(sim.monitors[0], ax=ax, cbar=True, comp='y', val='abs', eps_alpha=0.2)
plt.show()
../_images/examples_ParameterScan_30_0.png

Final Remarks

Batches provide some other convenient functionality for managing large numbers of jobs.

For example, one can save the batch information to file and load the batch at a later time, if needing to disconnect from the service while the jobs are running.

For more reference, please refer to our documentation.