Performing parallel / batch processing of simulations#


Note: the cost of running the entire notebook is larger than 1 FlexCredit.

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.

Note: as of version 1.8, the tidy3d.web.run_async function handles the same functionality as the batch, with a simpler syntax. As such, it could be a good alternative for parameter scan depending on how your script is set up.

For demonstration, we look at the splitting ratio of a directional coupler as we vary the coupling length between two waveguides. The sidewall of the waveguides is slanted, deviating from the vertical direction by sidewall_angle.

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

# tidy3D imports
import tidy3d as td
from tidy3d import web


First we set up some global parameters

# 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(permittivity=wg_n**2)
mat_sub = td.Medium(permittivity=sub_n**2)

# Waveguide dimensions

# Waveguide height
wg_height = 0.22
# Waveguide width
wg_width = 0.45
# Waveguide separation in the beginning/end
wg_spacing_in = 8
# Reference plane where the cross section of the device is defined
reference_plane = "bottom"
# Angle of the sidewall deviating from the vertical ones, positive values for the base larger than the top
sidewall_angle = np.pi / 6
# 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
# resolution control: minimum number of grid cells per wavelength in each material
grid_cells_per_wvl = 16

Define waveguide bends and coupler#

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

def tanh_interp(max_arg):
    """Interpolator for tanh with adjustable extension"""
    scale = 1 / np.tanh(max_arg)
    return lambda u: 0.5 * (1 + scale * np.tanh(max_arg * (u * 2 - 1)))

def make_coupler(
    """Make an integrated coupler using the gdstk RobustPath object."""
    # bend interpolator
    interp = tanh_interp(3)
    delta = wg_width + wg_spacing_coup - wg_spacing_in
    offset = lambda u: wg_spacing_in + interp(u) * delta

    coup = gdstk.RobustPath(
        (-0.5 * length, 0),
        (wg_width, wg_width),
        datatype=[0, 1],
    coup.segment((-0.5 * coup_length - bend_length, 0))
        (-0.5 * coup_length, 0),
        offset=[lambda u: -0.5 * offset(u), lambda u: 0.5 * offset(u)],
    coup.segment((0.5 * coup_length, 0))
        (0.5 * coup_length + bend_length, 0),
        offset=[lambda u: -0.5 * offset(1 - u), lambda u: 0.5 * offset(1 - u)],
    coup.segment((0.5 * length, 0))
    return coup

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.

def make_sim(coup_length, wg_spacing_coup, domain_field=False):
    """Make a simulation with a given length of the coupling region and
    distance between the waveguides in that region. If ``domain_field``
    is True, a 2D in-plane field monitor will be added.

    # Geometry must be placed in GDS cells to import into Tidy3D
    coup_cell = gdstk.Cell("Coupler")

    substrate = gdstk.rectangle(
        (-device_length / 2, -wg_spacing_in / 2 - 10),
        (device_length / 2, wg_spacing_in / 2 + 10),

    # Add the coupler to a gdstk cell
    gds_coup = make_coupler(

    # Substrate
    (oxide_geo,) = td.PolySlab.from_gds(
        slab_bounds=(-10, 0),

    oxide = td.Structure(geometry=oxide_geo, medium=mat_sub)

    # Waveguides (import all datatypes if gds_dtype not specified)
    coupler1_geo, coupler2_geo = td.PolySlab.from_gds(
        slab_bounds=(0, wg_height),

    coupler1 = td.Structure(geometry=coupler1_geo, medium=mat_wg)

    coupler2 = td.Structure(geometry=coupler2_geo, medium=mat_wg)

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

    # Spacing between waveguides and PML
    sim_size = [
        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(freq0=freq0, fwidth=fwidth),

    mon_in = td.ModeMonitor(
        center=[(src_pos + 0.5), wg_spacing_in / 2, wg_height / 2],
        size=[0, 3, 2],
    mon_ref_bot = td.ModeMonitor(
        center=[(src_pos + 0.5), -wg_spacing_in / 2, wg_height / 2],
        size=[0, 3, 2],
    mon_top = td.ModeMonitor(
        center=[-(src_pos + 0.5), wg_spacing_in / 2, wg_height / 2],
        size=[0, 3, 2],
    mon_bot = td.ModeMonitor(
        center=[-(src_pos + 0.5), -wg_spacing_in / 2, wg_height / 2],
        size=[0, 3, 2],
    monitors = [mon_in, mon_ref_bot, mon_top, mon_bot]

    if domain_field == True:
        domain_monitor = td.FieldMonitor(
            center=[0, 0, wg_height / 2],
            size=[td.inf, td.inf, 0],

    # initialize the simulation
    sim = td.Simulation(
        structures=[oxide, coupler1, coupler2],
        run_time=50 / fwidth,

    return sim

Inspect Simulation#

Let’s create and inspect a single simulation to make sure it was defined correctly before doing the full scan. The sidewalls of the waveguides deviate from the vertical direction by 30 degree. We also add an in-plane field monitor to have a look at the fields evolution in this one simulation. We will not use such a monitor in the batch to avoid storing unnecesarrily large amounts of data.

# 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, domain_field=True)

14:13:28 EDT WARNING: Default value for the field monitor 'colocate' setting has
             changed to 'True' in Tidy3D 2.4.0. All field components will be    
             colocated to the grid boundaries. Set to 'False' to get the raw    
             fields on the Yee grid instead.                                    
# visualize geometry
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 4))
sim.plot(z=wg_height / 2 + 0.01, ax=ax1)
sim.plot(x=0.1, ax=ax2)
ax2.set_xlim([-3, 3])


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.

# create job, upload sim to server to begin running
job = web.Job(simulation=sim, task_name="CouplerVerify", verbose=True)

# download the results and load them into a simulation
sim_data ="data/sim_data.hdf5")

14:13:29 EDT Created task 'CouplerVerify' with task_id
             'fdve-2a4191a7-5424-47c1-a71c-da54cf33977fv1' and task_type 'FDTD'.
14:13:31 EDT status = queued
14:13:45 EDT status = preprocess
14:13:52 EDT Maximum FlexCredit cost: 0.570. Use 'web.real_cost(task_id)' to get
             the billed FlexCredit cost after a simulation run.
             starting up solver
             running solver
             To cancel the simulation, use 'web.abort(task_id)' or
             'web.delete(task_id)' or abort/delete the task in the web UI.
             Terminating the Python script will not stop the job running on the
14:17:16 EDT early shutoff detected, exiting.
             status = postprocess
14:17:26 EDT status = success
14:17:32 EDT loading simulation from data/sim_data.hdf5


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.

def measure_transmission(sim_data):
    """Constructs a "row" of the scattering matrix when sourced from top left port"""

    input_amp = sim_data["in"].amps.sel(direction="+")

    amps = np.zeros(4, dtype=complex)
    directions = ("-", "-", "+", "+")
    for i, (monitor, direction) in enumerate(
        zip(sim_data.simulation.monitors[:4], directions)
        amp = sim_data[].amps.sel(direction=direction)
        amp_normalized = amp / input_amp
        amps[i] = np.squeeze(amp_normalized.values)

    return amps

# monitor and test out the measure_transmission function the results of the single run
amps_arms = measure_transmission(sim_data)
print("mode amplitudes in each port:\n")
for amp, monitor in zip(amps_arms, sim_data.simulation.monitors[:-1]):
    print(f'\tmonitor     = "{}"')
    print(f"\tamplitude^2 = {abs(amp)**2:.2f}")
    print(f"\tphase       = {(np.angle(amp)):.2f} (rad)\n")

mode amplitudes in each port:

        monitor     = "in"
        amplitude^2 = 0.00
        phase       = -2.01 (rad)

        monitor     = "refect_bottom"
        amplitude^2 = 0.00
        phase       = 0.76 (rad)

        monitor     = "top"
        amplitude^2 = 0.95
        phase       = -0.37 (rad)

        monitor     = "bottom"
        amplitude^2 = 0.04
        phase       = 1.21 (rad)

fig, ax = plt.subplots(1, 1, figsize=(16, 3))
sim_data.plot_field("field", "Ey", z=wg_height / 2, freq=freq0, ax=ax)

             WARNING: 'freq' supplied to 'plot_field', frequency selection key  
             renamed to 'f' and 'freq' will error in future release, please     
             update your local script to use 'f=value'.                         

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.

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

# submit all jobs
sims = {f"l={l:.2f}": make_sim(l, wg_spacing_coup) for l in ls}
batch = web.Batch(simulations=sims, verbose=True)

14:17:34 EDT Created task 'l=5.00' with task_id
             'fdve-13b89b71-9f45-4099-a898-2606cb2c7e28v1' and task_type 'FDTD'.
14:17:35 EDT Created task 'l=5.70' with task_id
             'fdve-2bd47601-4724-442f-9365-234e9f6a13a3v1' and task_type 'FDTD'.
             Created task 'l=6.40' with task_id
             'fdve-df571a70-1377-43dd-82cb-04d633a24df2v1' and task_type 'FDTD'.
14:17:36 EDT Created task 'l=7.10' with task_id
             'fdve-914f9472-87e9-4e7b-91a4-ca52eb0e8693v1' and task_type 'FDTD'.
14:17:37 EDT Created task 'l=7.80' with task_id
             'fdve-f39701aa-1da0-4d70-8241-4dd418309ac7v1' and task_type 'FDTD'.
14:17:38 EDT Created task 'l=8.50' with task_id
             'fdve-77765c7e-852d-4d13-8490-34acda934185v1' and task_type 'FDTD'.
14:17:39 EDT Created task 'l=9.20' with task_id
             'fdve-88b3ce05-af0b-4757-aa6d-d79efaf9f6c7v1' and task_type 'FDTD'.
             Created task 'l=9.90' with task_id
             'fdve-b5c19a5f-994a-4992-a1bd-d1baee3c6522v1' and task_type 'FDTD'.
14:17:40 EDT Created task 'l=10.60' with task_id
             'fdve-72d051b4-93ca-4d8e-92d1-f1cbcd635e93v1' and task_type 'FDTD'.
14:17:41 EDT Created task 'l=11.30' with task_id
             'fdve-52d604a9-d9be-40a1-a3c4-4ddceb9c61abv1' and task_type 'FDTD'.
14:17:42 EDT Created task 'l=12.00' with task_id
             'fdve-346473bf-5a80-4db9-91dd-3e8f97bd62ebv1' and task_type 'FDTD'.

Monitor Batch#

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

batch_results ="data")

14:17:46 EDT Started working on Batch.
14:18:58 EDT Maximum FlexCredit cost: 6.043 for the whole batch. Use
             'Batch.real_cost()' to get the billed FlexCredit cost after the
             Batch has completed.
14:25:37 EDT Batch complete.

Load and visualize Results#

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

Then we may plot the results.

amps_batch = []
for task_name, sim_data in batch_results.items():
    amps_arms_i = np.array(measure_transmission(sim_data))
amps_batch = np.stack(amps_batch, axis=1)
print(amps_batch.shape)  # (4, Nl)

14:25:39 EDT loading simulation from
14:25:40 EDT loading simulation from
14:25:41 EDT loading simulation from
14:25:42 EDT loading simulation from
14:25:43 EDT loading simulation from
14:25:44 EDT loading simulation from
14:25:45 EDT loading simulation from
14:25:46 EDT loading simulation from
14:25:47 EDT loading simulation from
14:25:48 EDT loading simulation from
14:25:49 EDT loading simulation from
(4, 11)
[[-6.88294062e-03+1.61016672e-04j  2.80545910e-03-3.25850838e-03j
  -6.05623712e-03-1.35837074e-03j -2.54948119e-05+3.26426641e-03j
  -2.64603665e-03-1.04155616e-02j -3.92398865e-03+3.56062722e-03j
   7.26028148e-04-4.00447556e-03j -4.50129367e-03-8.15524753e-04j
   6.45059889e-04-2.07259622e-03j -4.68488718e-03-1.09918831e-03j
 [ 7.60738204e-03-1.41929535e-03j -3.56046381e-03+3.17859946e-03j
   2.74558993e-03+2.55031515e-03j  1.24389155e-03+9.01951524e-04j
   1.17309161e-03+8.00086485e-03j  2.06995208e-03-2.68645085e-03j
  -1.63114672e-03+8.07276798e-03j  4.17566549e-03-1.16747389e-03j
  -2.49207660e-03+2.07816053e-03j  2.59925766e-03+7.49679312e-03j
 [ 4.47001365e-01+4.22382004e-01j  6.76803732e-01-2.61152916e-01j
   2.86673130e-02-8.19183765e-01j -8.13009030e-01-3.75187854e-01j
  -7.02815644e-01+6.40098058e-01j  3.28939822e-01+9.27131908e-01j
   9.91166611e-01+6.25283978e-02j  4.60598499e-01-8.64562351e-01j
  -5.90526169e-01-7.34374594e-01j -8.50714066e-01+2.34752546e-01j
 [ 5.36546422e-01-5.62458676e-01j -2.40372525e-01-6.31713375e-01j
  -5.57351973e-01-2.21623999e-02j -1.80638765e-01+3.86319694e-01j
   1.90658062e-01+2.11407195e-01j  1.28452325e-01-4.50088044e-02j
  -9.56415432e-04+1.63354725e-02j  1.47599188e-01+7.88078674e-02j
   2.45300488e-01-1.97485389e-01j -1.21799454e-01-4.38410802e-01j
powers = abs(amps_batch) ** 2
power_top = powers[2]
power_bot = powers[3]
power_out = power_top + power_bot

plt.plot(ls, 100 * power_top, label="Top")
plt.plot(ls, 100 * power_bot, label="Bottom")
plt.plot(ls, 100 * power_out, label="Top + Bottom")
plt.xlabel("Coupling length (µm)")
plt.ylabel("Power ratio (%)")
plt.ylim(0, 100)


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.

# save batch metadata

# load batch metadata into a new batch
loaded_batch = web.Batch.from_file("data/batch_data.json")

For more reference, refer to our documentation.

Parameter Scan using Signac#

For users who might need more structure to their parameter scans, the open source tool “signac” is a fantastic resource. Here we will reproduce some of the previous parameter scan using signac to give an introduction to how to apply it to Tidy3D projects. You can see detailed tutorials and examples at their documentation page.

After importing the package, we need to define a project, which also tells signac to store our parameter scan data in a local directory of our choice.

import signac

# make the project and give it a path to save the data into
project = signac.init_project(path="data/splitter")

With siganc (and more generally in parameter sweeps), it is very useful to have a single function to define the inputs and outputs of your parameter scan. In our case, we have a single input (coupling length l) and two outputs: coupling efficiency (eff) and the split ratio (ratio). We write a function to link these inputs and outputs through a Tidy3D simulation.

def compute(l: float):
    """compute efficiency and split ratio as a function of the coupling length"""
    sim = make_sim(l, wg_spacing_coup)
    task_name = f"SWEEP_l={l:.3f}"
    sim_data =, task_name=task_name, verbose=False)
    amps_arms = np.array(measure_transmission(sim_data))
    powers = np.abs(amps_arms)**2
    efficiency = np.sum(powers)
    ratio_0 = powers[2] / efficiency
    return efficiency, ratio_0

The signac project contains a set of job objects, which each define a specific data point of our parameter scan. As such, we define another function that wraps compute but instead simply accepts a job instance, which makes our lives easier for managing parameter sweeps in signac.

This function basically will open a job, grab its inputs (l value in our case), compute the output quantities, and save them to both the job.document record and also the .txt files storing the parameter scan results.

def compute_transmission(job):
    l = job.statepoint()["l"]
    print(f"Computing transmission of job: [l={l:.1f}]", job)
    eff, ratio = compute(**job.statepoint())
    job.document["eff"] = float(eff)
    job.document["ratio"] = float(ratio)
    with open(job.fn("eff.txt"), "w") as file:
        file.write(str(float(eff)) + "\n")
    with open(job.fn("ratio.txt"), "w") as file:
        file.write(str(float(ratio)) + "\n")

With this, we can start our parameter scan, we simply loop through our desired l values, construct a job for each, and pass it to our compute_transmission function.

Note: If intermediate simulation results are needed, it should be possible to modify compute to return also the SimulationData associated with our task, which we could then store in the job.document or elsewhere.

We will compute only 3 values between our ranges of 5 and 12 to save time.

Note, the job instances each have their own unique ID, similar to Tidy3D task_id, this is what is being printed below with each job being computed.

for l in np.linspace(5, 12, 3):
    statepoint = {"l": float(l)}
    job = project.open_job(statepoint)
Computing transmission of job: [l=5.0] abb6bcf88ee5474914d016b3c1fcc7b3
Computing transmission of job: [l=8.5] c261dda788d2af5c1645cb17feb7854a
Computing transmission of job: [l=12.0] bc0f33313f39cc49b60566dcb7f169f9

Let’s take a look at our project and what data we’ve computed and stored so far.

for job in project:
    print(job.statepoint(), job.document)
{'l': 8.5} {'eff': 0.9863404768883945, 'ratio': 0.981177516614684}
{'l': 5.0} {'eff': 0.9825661842137496, 'ratio': 0.38492766600414147}
{'l': 12.0} {'eff': 0.9848524347794334, 'ratio': 0.65295559679756}

signac also provides the ability to create a “linked view”, which gives us a human-readable filesystem for looking at our jobs.


# let's print out some info stored in the directory we just created
!ls data/splitter/view/l/
!ls data/splitter/view/l/5.0/job/
!cat data/splitter/view/l/5.0/job/eff.txt
12.0 5.0  8.5
eff.txt                  signac_job_document.json
ratio.txt                signac_statepoint.json

We can also initialize many data points to compute later, and feed them to the compute_transmission function with a very basic “caching” implemented by checking if our inputs already exist in the record.

def init_statepoints(n):
    for l in np.linspace(5, 12, n):
        sp = {"l": float(l)}
        job = project.open_job(sp)
        print("initialize", job)

# make 5 points between 5 and 12, note, 3 have already been computed
initialize abb6bcf88ee5474914d016b3c1fcc7b3
initialize 07222394eb1a1b59f9e869ef871a3c29
initialize c261dda788d2af5c1645cb17feb7854a
initialize bdbbc3c4500db5773fa3cd350b6195ed
initialize bc0f33313f39cc49b60566dcb7f169f9

After initializing our statepoints (input values), they are stored in our project and we can loop through them and compute any that don’t have records.

for job in project:
    if "eff" not in job.document or "ratio" not in job.document:
        print(f" - skipping job: {job} ")
 - skipping job: c261dda788d2af5c1645cb17feb7854a
 - skipping job: abb6bcf88ee5474914d016b3c1fcc7b3
Computing transmission of job: [l=6.8] 07222394eb1a1b59f9e869ef871a3c29
Computing transmission of job: [l=10.2] bdbbc3c4500db5773fa3cd350b6195ed
 - skipping job: bc0f33313f39cc49b60566dcb7f169f9

While we used td.web.Batch in our original example, signac lets us also leverage parallel processing tools in python to perform something similar.

Let’s initialize 9 total statepoints now (5 have already been computed) and feed them to a ThreadPool. We notice that the jobs will be computed in parallel depending on how many threads are available on your machine.


from multiprocessing.pool import ThreadPool

# make a convenience function to just call compute_transmission only for uncomputed jobs
def compute_transmission_cached(job):
    if "eff" not in job.document or "ratio" not in job.document:

with ThreadPool() as pool:, list(project))
initialize abb6bcf88ee5474914d016b3c1fcc7b3
initialize dca27bcde6b7bfede0d701809749e789
initialize 07222394eb1a1b59f9e869ef871a3c29
initialize b22f487baf2db695e6b15444d42e31ff
initialize c261dda788d2af5c1645cb17feb7854a
initialize c7a1813a4fd1fd835898df062c3262c9
initialize bdbbc3c4500db5773fa3cd350b6195ed
initialize 9af2b1986858046540a9da68dbf4f425
initialize bc0f33313f39cc49b60566dcb7f169f9
Computing transmission of job: [l=7.6] b22f487baf2db695e6b15444d42e31ff
Computing transmission of job: [l=11.1] 9af2b1986858046540a9da68dbf4f425
Computing transmission of job: [l=5.9] dca27bcde6b7bfede0d701809749e789
Computing transmission of job: [l=9.4] c7a1813a4fd1fd835898df062c3262c9
for job in project:
    print(job.statepoint(), job.document)
{'l': 8.5} {'eff': 0.9863404768883945, 'ratio': 0.981177516614684}
{'l': 7.625} {'eff': 0.9848879706230304, 'ratio': 0.8947588124079515}
{'l': 5.0} {'eff': 0.9825661842137496, 'ratio': 0.38492766600414147}
{'l': 6.75} {'eff': 0.9824394427673719, 'ratio': 0.7526833195705787}
{'l': 11.125} {'eff': 0.9864516067050044, 'ratio': 0.8200631095893199}
{'l': 5.875} {'eff': 0.9825398650603947, 'ratio': 0.5741899484506372}
{'l': 10.25} {'eff': 0.9870634987508813, 'ratio': 0.9402364889170417}
{'l': 9.375} {'eff': 0.9873992720742533, 'ratio': 0.9971436153934053}
{'l': 12.0} {'eff': 0.9848524347794334, 'ratio': 0.65295559679756}

Finally, we can consolidate and plot our results.

ls =  np.array([job.statepoint()["l"] for job in project])
effs = np.array([job.document["eff"] for job in project])
ratios = np.array([job.document["ratio"] for job in project])
plt.scatter(ls, 100 * ratios, label="split ratio")
plt.scatter(ls, 100 * effs, label="efficiency")
plt.xlabel("Coupling length (µm)")
plt.ylabel("value (%)")
plt.ylim(0, 100)


For more information, we highly recommend you visit signac documentation page, which includes explanations about all of the many other features not covered here, which could assist you in your parameter scans using Tidy3D.

[ ]: