Inverse design optimization of a bilayer grating coupler#
This example demonstrates the co-optimization of two separate design regions using tidy3d inverse design. Instead of the adjoint plugin, we use the native automatic differentiation support built into tidy3d versions 2.7 and later. In this framework, regular tidy3d components and web.run() functions are differentiable without any change of syntax.
We optimize a grating coupler for coupling efficiency with respect to
the pattern partially etched into a Si coupling layer.
the pattern etched into a SiN layer suspended above the coupling region.
Our example will roughly be following this paper.
[1]:
import autograd
import autograd.numpy as np
import matplotlib.pylab as plt
import tidy3d as td
import tidy3d.web as web
# define some length scales in tidy3d default of microns
um = 1.0
nm = 1e-3
Set up simulation#
First, we define everything “static” (not optimizing) in the simulation, including parameters and tidy3d components.
Spectral Parameters#
[2]:
wavelength = 1310 * nm
freq = td.C_0 / wavelength
fwidth = freq / 10
run_time = 50 / fwidth
Geometric Parameters#
[3]:
# device parameters
thick_substrate = 2.0 * um
thick_Si = 220 * nm
thick_SiN = 600 * nm
thick_etch = 160 * nm
space_SiSiN = 100 * nm
width_waveguide = 500 * nm
# size of design region
size_design_x = 4 * um
size_design_y = 4 * um
# spacings between things (technically infinite, but finite for sim)
length_waveguide = 1 * wavelength
thick_Si_substrate = 550 * nm
space_cladding_top = 1 * wavelength
space_design_edge = 1 * wavelength
# compute sim size
size_sim_x = length_waveguide + size_design_x + space_design_edge
size_sim_y = space_design_edge + size_design_y + space_design_edge
size_sim_z = (
thick_Si_substrate + thick_substrate + thick_Si + space_SiSiN + thick_SiN + space_cladding_top
)
size_sim = (size_sim_x, size_sim_y, size_sim_z)
# references to various locations
sim_top_z = size_sim_z / 2.0
sim_bot_z = -size_sim_z / 2.0
center_SiN_z = sim_top_z - space_cladding_top - thick_SiN / 2.0
center_Si_z = sim_top_z - space_cladding_top - thick_SiN - space_SiSiN - thick_Si / 2.0
center_Si_etch_z = sim_top_z - space_cladding_top - thick_SiN - space_SiSiN - thick_etch / 2.0
center_source_z = sim_top_z - space_cladding_top / 2.0
center_Si_substrate_z = sim_bot_z + thick_Si_substrate / 2
center_waveguide_x = -size_sim_x / 2 + length_waveguide / 2.0
center_design_x = size_sim_x / 2.0 - space_design_edge - size_design_x / 2.0
Material Parameters#
[4]:
# dispersionless materials for now. note: can also do `td.Medium.from_nk(n=n, k=k, freq=freq)`
n_SiO2 = 1.44
n_SiN = 2.0
n_Si = 3.5
eps_Si = n_Si**2
eps_SiN = n_SiN**2
eps_SiO2 = n_SiO2**2
medium_Si = td.Medium(permittivity=eps_Si)
medium_SiN = td.Medium(permittivity=eps_SiN)
medium_SiO2 = td.Medium(permittivity=eps_SiO2)
Source Parameters#
[5]:
spot_size = np.sqrt(size_design_x**2 + size_design_y**2)
tilt_angle_deg = 0.0
Static Components#
[6]:
waveguide = td.Structure(
geometry=td.Box(
center=(-size_sim_x + length_waveguide + size_design_x / 2.0, 0, center_Si_z),
size=(size_sim_x, width_waveguide, thick_Si),
),
medium=medium_Si,
)
design_region = td.Structure(
geometry=td.Box(
center=(center_design_x, 0, center_Si_z),
size=(size_design_x, size_design_y, thick_Si),
),
medium=medium_Si,
)
slab_SiN = td.Structure(
geometry=td.Box(
center=(center_design_x, 0, center_SiN_z),
size=(size_design_x, size_design_y, thick_SiN),
),
medium=medium_SiN,
)
substrate_Si = td.Structure(
geometry=td.Box(
center=(0, 0, -size_sim_z + thick_Si_substrate),
size=(td.inf, td.inf, size_sim_z),
),
medium=medium_Si,
)
gaussian_beam = td.GaussianBeam(
center=(center_design_x, 0, center_source_z),
size=(size_design_x, size_design_y, 0),
source_time=td.GaussianPulse(freq0=freq, fwidth=fwidth),
pol_angle=np.pi / 2,
angle_theta=tilt_angle_deg * np.pi / 180.0,
direction="-",
num_freqs=7,
waist_radius=spot_size / 2,
)
mode_monitor = td.ModeMonitor(
size=(0, width_waveguide * 6, thick_Si * 6),
center=(center_waveguide_x, 0, center_Si_z),
name="mode",
freqs=[freq],
mode_spec=td.ModeSpec(),
)
field_monitor_xy = td.FieldMonitor(
center=(0, 0, center_Si_z),
size=(td.inf, td.inf, 0),
freqs=[freq],
name="field_xy",
)
field_monitor_xz = td.FieldMonitor(
center=(0, 0, 0),
size=(td.inf, 0, td.inf),
freqs=[freq],
name="field_xz",
)
sim_no_etch = td.Simulation(
size=size_sim,
run_time=run_time,
medium=medium_SiO2,
structures=[waveguide, design_region, slab_SiN, substrate_Si],
grid_spec=td.GridSpec.auto(wavelength=wavelength),
sources=[gaussian_beam],
monitors=[mode_monitor, field_monitor_xy, field_monitor_xz],
)
[7]:
# shift the plot positions a little so the field monitors don't overlap plots
shift_plot = 0.01
f, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(10, 4), tight_layout=True)
sim_no_etch.plot(x=0.0 + shift_plot, ax=ax1)
sim_no_etch.plot(y=0.0 + shift_plot, ax=ax2)
sim_no_etch.plot(z=center_Si_z + shift_plot, ax=ax3)
plt.show()
Inspect waveguide modes#
We will run the mode solver to determine the proper mode_index corresponding to the fundamental TE mode.
[8]:
from tidy3d.plugins.mode import ModeSolver
from tidy3d.plugins.mode.web import run as run_mode_solver
num_modes = 3
mode_spec = td.ModeSpec(num_modes=num_modes)
mode_solver = ModeSolver(
simulation=sim_no_etch,
plane=mode_monitor.geometry,
mode_spec=mode_spec,
freqs=[freq],
)
modes = run_mode_solver(mode_solver, reduce_simulation=True)
21:17:22 CEST Mode solver created with task_id='fdve-d206abf6-79aa-4848-8390-c86411cce662', solver_id='mo-9d65902a-3151-40cf-9a9c-eca9f176e23b'.
21:17:26 CEST Mode solver status: queued
21:17:28 CEST Mode solver status: running
21:17:34 CEST Mode solver status: success
[9]:
fig, axs = plt.subplots(num_modes, 3, figsize=(10, 6), tight_layout=True)
for mode_index in range(num_modes):
vmax = 1.1 * max(
abs(modes.field_components[key].sel(mode_index=mode_index)).max()
for key in ("Ex", "Ey", "Ez")
)
for field_name, ax in zip(("Ex", "Ey", "Ez"), axs[mode_index]):
field = modes.field_components[field_name].sel(mode_index=mode_index)
field.real.squeeze().plot.imshow(x="y", y="z", label="Real", ax=ax)
ax.set_title(f"index={mode_index}, {field_name}")
ax.set_aspect("equal")
print("Effective index of computed modes: ", np.array(modes.n_eff))
# mode_index=0 corresponds to our mode of interest, so we'll keep our mode spec with default num_modes=1.
Effective index of computed modes: [[2.70037825 2.18933488 1.81996429]]
Include Etch Regions#
Next, we’ll define the functions describing our two etched design regions.
Starting from raw optimization parameters, we will use filtering and projection techniques to obtain an “etch density” pattern. We then rescale this to the etched and un-etched permittivity values. Finally, we package the etch pattern as a Structure containing a CustomMedium.
[10]:
from tidy3d.plugins.autograd import make_filter_and_project, rescale
# resolution of the design region pixels
pixel_size_Si = 15 * nm
pixel_size_SiN = 15 * nm
# radius of the circular filter (um) (higher = larger features)
radius_Si = 100 * nm
radius_SiN = 100 * nm
# projection strengths (higher = more binarized)
beta_Si = 5
beta_SiN = 5
beta_penalty = 5
nx_Si = int(size_design_x // pixel_size_Si)
ny_Si = int(size_design_y // pixel_size_Si)
nx_SiN = int(size_design_x // pixel_size_SiN)
ny_SiN = int(size_design_y // pixel_size_SiN)
etch_Si_geometry = geometry = td.Box(
size=(size_design_x, size_design_y, thick_etch),
center=(center_design_x, 0, center_Si_etch_z),
)
def get_density(params: np.ndarray, radius: float, beta: float, pixel_size: float) -> np.ndarray:
"""Generic function to get the etch density as a function of design parameters, using filter and projection."""
filter_project = make_filter_and_project(radius, pixel_size, beta=beta)
return filter_project(params)
def get_density_Si(params: np.ndarray) -> td.Structure:
"""Get the density of the Si etch as a function of its design parameters."""
return get_density(params=params, radius=radius_Si, beta=beta_Si, pixel_size=pixel_size_Si)
def get_density_SiN(params: np.ndarray) -> td.Structure:
"""Get the density of the SiN etch as a function of its design parameters."""
return get_density(params=params, radius=radius_SiN, beta=beta_SiN, pixel_size=pixel_size_SiN)
def get_permittivity(density: np.ndarray, eps_max: float) -> np.ndarray:
"""Get the permittivity array between SiO2 and the max, given an array of material densities between 0 and 1."""
return rescale(density, eps_SiO2, eps_max)
def make_etch_structure(density: np.ndarray, eps_max: float, geometry: td.Box) -> td.Structure:
"""Make a `td.Structure` containing a `td.CustomMedium` corresponding to this density array, given a geometry."""
permittivity_data = get_permittivity(density, eps_max=eps_max)
rmin, rmax = geometry.bounds
coords = {}
for key, pt_min, pt_max, num_pts in zip("xyz", rmin, rmax, density.shape):
coord_edges = np.linspace(pt_min, pt_max, num_pts + 1)
coord_centers = (coord_edges[1:] + coord_edges[:-1]) / 2.0
coords[key] = coord_centers
custom_medium = td.CustomMedium(
permittivity=td.ScalarFieldDataArray(
permittivity_data,
coords=coords,
)
)
return td.Structure(geometry=geometry, medium=custom_medium)
def get_etch_Si(params: np.ndarray) -> td.Structure:
"""Get the etch region for Si, using the Si parameters."""
density = get_density_Si(params)
return make_etch_structure(density=density, eps_max=eps_Si, geometry=etch_Si_geometry)
def get_etch_SiN(params: np.ndarray) -> td.Structure:
"""Get the etch region for SiN, using the SiN parameters."""
density = get_density_SiN(params)
return make_etch_structure(density=density, eps_max=eps_SiN, geometry=slab_SiN.geometry)
Optimization Parameterization#
Next we will write a few other functions to help define our Simulation itself from the design parameters.
[11]:
def make_sim_with_etch(params_Si: np.ndarray, params_SiN: np.ndarray) -> td.Simulation:
"""Make a simulation as a function of the design parameters for Si and SiN etch regions."""
etch_Si = get_etch_Si(params_Si)
etch_SiN = get_etch_SiN(params_SiN)
# add uniform mesh override structures to simulation
design_region_mesh_Si = td.MeshOverrideStructure(
geometry=design_region.geometry,
dl=[pixel_size_Si] * 3,
enforce=True,
)
design_region_mesh_SiN = td.MeshOverrideStructure(
geometry=slab_SiN.geometry,
dl=[pixel_size_SiN] * 3,
enforce=True,
)
grid_spec = sim_no_etch.grid_spec.updated_copy(
override_structures=list(sim_no_etch.grid_spec.override_structures)
+ [design_region_mesh_Si, design_region_mesh_SiN]
)
return sim_no_etch.updated_copy(
structures=list(sim_no_etch.structures) + [etch_Si, etch_SiN],
grid_spec=grid_spec,
)
Let’s create some initial parameters to test our parameterization with.
[12]:
def make_symmetric_y(arr: np.ndarray) -> np.ndarray:
"""make an array symmetric in y."""
return (arr + np.fliplr(arr)) / 2.0
params0_Si = make_symmetric_y(np.random.random((nx_Si, ny_Si, 1)))
params0_SiN = make_symmetric_y(np.random.random((nx_SiN, ny_SiN, 1)))
sim_etch_random = make_sim_with_etch(params0_Si, params0_SiN)
f, axes = f, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(10, 4), tight_layout=True)
sim_etch_random.plot_eps(x=0.0 + shift_plot, monitor_alpha=0.0, ax=ax1)
sim_etch_random.plot_eps(y=0.0 + shift_plot, monitor_alpha=0.0, ax=ax2)
sim_etch_random.plot_eps(z=center_Si_z + shift_plot, monitor_alpha=0.0, ax=ax3)
for ax in axes:
ax.set_aspect("equal")
plt.show()
[13]:
sim_data = web.run(sim_etch_random, task_name="check fields")
21:17:38 CEST Created task 'check fields' with task_id 'fdve-39de21d3-987d-4d72-9b26-1805ec0527b0' and task_type 'FDTD'.
View task using web UI at 'https://tidy3d.simulation.cloud/workbench?taskId=fdve-39de21d3-98 7d-4d72-9b26-1805ec0527b0'.
Task folder: 'default'.
21:17:44 CEST Maximum FlexCredit cost: 0.477. Minimum cost depends on task execution details. Use 'web.real_cost(task_id)' to get the billed FlexCredit cost after a simulation run.
21:17:45 CEST status = queued
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 cloud.
21:18:01 CEST status = preprocess
21:18:07 CEST starting up solver
running solver
21:18:24 CEST early shutoff detected at 20%, exiting.
status = postprocess
21:18:27 CEST status = success
21:18:29 CEST View simulation result at 'https://tidy3d.simulation.cloud/workbench?taskId=fdve-39de21d3-98 7d-4d72-9b26-1805ec0527b0'.
21:18:32 CEST loading simulation from simulation_data.hdf5
[14]:
f, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 4), tight_layout=True)
sim_data.plot_field("field_xy", field_name="E", val="abs^2", ax=ax1)
sim_data.plot_field("field_xz", field_name="E", val="abs^2", ax=ax2)
plt.show()
Define Objective Function#
The next step involves defining our objective function, which we seek to maximize with respect to the parameters.
We’ll compute the coupling efficiency for the grating coupler, and also include a penalty for each of the design regions depending on how well they satisfy minimum feature size criteria.
[15]:
from tidy3d.plugins.autograd import make_erosion_dilation_penalty
penalty_fn_Si = make_erosion_dilation_penalty(radius_Si, pixel_size_Si, beta=beta_penalty)
penalty_fn_SiN = make_erosion_dilation_penalty(radius_SiN, pixel_size_SiN, beta=beta_penalty)
def penalty_Si(params: np.ndarray) -> float:
"""Define the erosion dilation invariance penalty for Si density parameters."""
density = get_density_Si(params)
return penalty_fn_Si(density)
def penalty_SiN(params: np.ndarray) -> float:
"""Define the erosion dilation invariance penalty for SiN density parameters."""
density = get_density_SiN(params)
return penalty_fn_SiN(density)
def coupling_efficiency(sim_data: td.SimulationData) -> float:
"""Coupling efficiency into the waveguide as a function of the simulation output data."""
output_amps = sim_data["mode"].amps
amp = output_amps.sel(direction="-", f=freq, mode_index=0).values
return np.sum(np.abs(amp) ** 2)
[16]:
def objective(
params_Si: np.ndarray,
params_SiN,
verbose: bool = False,
include_field_mnts: bool = False,
) -> float:
"""Combined objective function over the full set of parameters."""
# coupling efficiency calculation, through a differentiable simulation
sim = make_sim_with_etch(params_Si, params_SiN)
if not include_field_mnts:
sim = sim.updated_copy(monitors=[mode_monitor])
data = web.run(sim, task_name="coupler", verbose=verbose)
efficiency = coupling_efficiency(data)
# fabrication penalty calculation for both design regions
penalty_val_Si = penalty_Si(params_Si)
penalty_val_SiN = penalty_SiN(params_SiN)
# total objective (to maximize)
coupling_objective = 2 * efficiency
penalty_objective = 0.5 * (penalty_val_Si + penalty_val_SiN)
return coupling_objective - penalty_objective
Let’s test this out and inspect the gradients for each design region as a function of space.
[17]:
val_grad_fn = autograd.value_and_grad(objective, argnum=(0, 1))
[18]:
val, grads = val_grad_fn(params0_Si, params0_SiN, verbose=True)
21:18:33 CEST Created task 'coupler' with task_id 'fdve-f8b0e55c-da99-4b0e-9712-304cdc7c448a' and task_type 'FDTD'.
View task using web UI at 'https://tidy3d.simulation.cloud/workbench?taskId=fdve-f8b0e55c-da 99-4b0e-9712-304cdc7c448a'.
Task folder: 'default'.
21:18:36 CEST Maximum FlexCredit cost: 0.486. Minimum cost depends on task execution details. Use 'web.real_cost(task_id)' to get the billed FlexCredit cost after a simulation run.
21:18:37 CEST status = queued
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 cloud.
21:19:02 CEST status = preprocess
21:19:09 CEST starting up solver
running solver
21:19:21 CEST early shutoff detected at 20%, exiting.
status = postprocess
21:19:33 CEST status = success
21:19:35 CEST View simulation result at 'https://tidy3d.simulation.cloud/workbench?taskId=fdve-f8b0e55c-da 99-4b0e-9712-304cdc7c448a'.
21:19:45 CEST loading simulation from simulation_data.hdf5
21:19:48 CEST Started working on Batch containing 1 tasks.
21:19:50 CEST Maximum FlexCredit cost: 0.486 for the whole batch.
Use 'Batch.real_cost()' to get the billed FlexCredit cost after the Batch has completed.
21:20:44 CEST Batch complete.
[19]:
grad_Si, grad_SiN = grads
print(f"objective function value = {val}")
f, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 4), tight_layout=True)
vmag1 = np.max(abs(grad_Si))
im1 = ax1.imshow(np.flipud(np.squeeze(grad_Si)).T, cmap="PiYG", vmax=vmag1, vmin=-vmag1)
plt.colorbar(im1, ax=ax1)
ax1.set_title("gradient w.r.t. Si parameters")
vmag2 = np.max(abs(grad_SiN))
im2 = ax2.imshow(np.flipud(np.squeeze(grad_SiN)).T, cmap="PiYG", vmax=vmag2, vmin=-vmag2)
plt.colorbar(im2, ax=ax2)
ax2.set_title("gradient w.r.t. SiN parameters")
plt.show()
objective function value = -0.8576508342522808
For reference: The green regions represent where we would need to locally increase parameters to increase the objective function and the purple regions represent where we would need to locally decrease parameters to increase our objective function.
Optimization#
Let’s use our objective function gradient calculator to perform optimization. Like in the adjoint notebooks, we use optax to do Adam optimization in a for loop, feeding it our gradient computed at each iteration.
[20]:
import optax
# hyperparameters
num_steps = 30
learning_rate = 0.25
def flatten_arrays(arr_Si: np.ndarray, arr_SiN: np.ndarray) -> np.ndarray:
"""Put arr into a 1D array for optax."""
return np.concatenate([arr_Si.flatten(), arr_SiN.flatten()])
def unflatten_array(arr: np.ndarray) -> tuple[np.ndarray, np.ndarray]:
"""Unflatten flattened array into Si and SiN arrays."""
arr_Si = arr[: params0_Si.size].reshape(params0_Si.shape)
arr_SiN = arr[params0_Si.size :].reshape(params0_SiN.shape)
return arr_Si, arr_SiN
def plot_density(density: np.ndarray, ax):
"""Plot the density of the device."""
arr = np.flipud(1 - density.squeeze().T)
ax.imshow(arr, cmap="gray", vmin=0, vmax=1, interpolation="none")
return ax
params = flatten_arrays(params0_Si, params0_SiN)
# initialize adam optimizer with starting parameters (all combined)
optimizer = optax.adam(learning_rate=learning_rate)
opt_state = optimizer.init(params)
# store history
objective_values_history = []
params_history = [(params0_Si, params0_SiN)]
# optimization loop
for i in range(num_steps):
# unpack parameters
params_Si, params_SiN = unflatten_array(np.array(params))
# plot the densities, to monitor
f, (ax1, ax2) = plt.subplots(1, 2, figsize=(4, 2))
plot_density(get_density_Si(params_Si), ax=ax1)
plot_density(get_density_SiN(params_SiN), ax=ax2)
ax1.set_title("Si etch density")
ax2.set_title("SiN etch density")
ax1.axis("off")
ax2.axis("off")
plt.show()
val, grads = val_grad_fn(params_Si, params_SiN, verbose=False)
grad_Si, grad_SiN = grads
gradient = flatten_arrays(grad_Si, grad_SiN)
# outputs
print(f"step = {i + 1}")
print(f"\tobjective = {val:.4e}")
print(f"\tgrad_norm = {np.linalg.norm(gradient):.4e}")
# compute and apply updates to the optimizer based on gradient (-1 sign to maximize obj_fn)
updates, opt_state = optimizer.update(-gradient, opt_state, params)
params[:] = optax.apply_updates(params, updates)
# clip the parameters between 0 and 1
np.clip(params, 0.0, 1.0, out=params)
# save history
objective_values_history.append(val)
params_history.append((params_Si, params_SiN))
step = 1
objective = -8.5765e-01
grad_norm = 2.5621e-03
step = 2
objective = -1.1963e-01
grad_norm = 2.1515e-02
step = 3
objective = 1.0083e-01
grad_norm = 3.1901e-02
step = 4
objective = 2.6022e-01
grad_norm = 3.6418e-02
step = 5
objective = 3.9154e-01
grad_norm = 3.0264e-02
step = 6
objective = 4.8566e-01
grad_norm = 2.6367e-02
step = 7
objective = 5.5079e-01
grad_norm = 1.9033e-02
step = 8
objective = 5.9703e-01
grad_norm = 2.2765e-02
step = 9
objective = 6.5384e-01
grad_norm = 1.5830e-02
step = 10
objective = 6.7951e-01
grad_norm = 2.3458e-02
step = 11
objective = 7.2627e-01
grad_norm = 1.4495e-02
step = 12
objective = 7.5166e-01
grad_norm = 1.0774e-02
step = 13
objective = 7.6559e-01
grad_norm = 1.3285e-02
step = 14
objective = 7.8622e-01
grad_norm = 1.1526e-02
step = 15
objective = 8.0981e-01
grad_norm = 9.1037e-03
step = 16
objective = 8.3001e-01
grad_norm = 7.8232e-03
step = 17
objective = 8.4326e-01
grad_norm = 8.2218e-03
step = 18
objective = 8.5366e-01
grad_norm = 7.5081e-03
step = 19
objective = 8.6417e-01
grad_norm = 7.3905e-03
step = 20
objective = 8.7278e-01
grad_norm = 8.7542e-03
step = 21
objective = 8.8494e-01
grad_norm = 6.6962e-03
step = 22
objective = 8.9623e-01
grad_norm = 5.4868e-03
step = 23
objective = 9.0515e-01
grad_norm = 7.0182e-03
step = 24
objective = 9.1622e-01
grad_norm = 6.7522e-03
step = 25
objective = 9.2721e-01
grad_norm = 5.6352e-03
step = 26
objective = 9.3764e-01
grad_norm = 5.2267e-03
step = 27
objective = 9.4586e-01
grad_norm = 5.7851e-03
step = 28
objective = 9.5139e-01
grad_norm = 9.7382e-03
step = 29
objective = 9.4952e-01
grad_norm = 2.2997e-02
step = 30
objective = 9.2232e-01
grad_norm = 3.8045e-02
Analysis#
Let’s compute the objective for the last set of parameters and plot our results.
[27]:
params_final = params_history[-3]
objective_value_final = objective(*params_final)
objective_values_history.append(objective_value_final)
[28]:
plt.plot(objective_values_history)
plt.xlabel("iterations")
plt.ylabel("objective function")
plt.show()
[29]:
sim_final = make_sim_with_etch(*params_final)
sim_data_final = web.run(sim_final, task_name="coupler final")
07:44:34 CEST Created task 'coupler final' with task_id 'fdve-9c2fe59a-9e45-4d0f-a6e0-7e3f729eb82b' and task_type 'FDTD'.
View task using web UI at 'https://tidy3d.simulation.cloud/workbench?taskId=fdve-9c2fe59a-9e 45-4d0f-a6e0-7e3f729eb82b'.
Task folder: 'default'.
07:44:37 CEST Maximum FlexCredit cost: 0.477. Minimum cost depends on task execution details. Use 'web.real_cost(task_id)' to get the billed FlexCredit cost after a simulation run.
07:44:38 CEST status = queued
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 cloud.
07:44:50 CEST status = preprocess
07:44:57 CEST starting up solver
running solver
07:45:12 CEST early shutoff detected at 20%, exiting.
status = postprocess
07:45:15 CEST status = success
07:45:17 CEST View simulation result at 'https://tidy3d.simulation.cloud/workbench?taskId=fdve-9c2fe59a-9e 45-4d0f-a6e0-7e3f729eb82b'.
07:45:20 CEST loading simulation from simulation_data.hdf5
[30]:
f, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, figsize=(12, 8), tight_layout=False)
params_Si, params_SiN = params_final
density_Si = get_density_Si(params_Si)
density_SiN = get_density_SiN(params_SiN)
ax1.imshow(np.flipud(1 - density_Si.squeeze().T), cmap="grey")
ax2.imshow(np.flipud(1 - density_SiN.squeeze().T), cmap="grey")
ax1.set_title("Si etch pattern")
ax2.set_title("SiN etch pattern")
sim_data_final.plot_field("field_xy", field_name="E", val="abs^2", ax=ax3)
sim_data_final.plot_field("field_xz", field_name="E", val="abs^2", ax=ax4)
ax3.set_title("intensity coupler (xy) plane")
ax4.set_title("intensity side (xz) plane")
plt.show()
[31]:
efficiency_final = coupling_efficiency(sim_data_final)
print(f"final coupling efficiency = {100 * efficiency_final:.3f}%")
penalty_val_Si = penalty_Si(params_Si)
penalty_val_SiN = penalty_SiN(params_SiN)
print(f"final penalty (Si) = {penalty_val_Si:.3f}")
print(f"final penalty (SiN) = {penalty_val_SiN:.3f}")
final coupling efficiency = 51.292%
final penalty (Si) = 0.127
final penalty (SiN) = 0.022
[26]:
# optionally save the coupler history if we want to load it later.
# import pickle
# with open("coupler.pkl", "wb") as file:
# pickle.dump(params_history, file)