Fabrication-aware inverse design¶
This notebook extends the Adjoint Optimization of a Wavelength Division Multiplexer by incorporating fabrication awareness into the design process. Specifically, we integrate a differentiable computer vision model into the topology optimization loop to predict and optimize the manufacturability of photonic device designs. This approach, known as fabrication-aware inverse design (FAID), results in devices that are inherently more robust to fabrication variations.
Unlike optical proximity correction (OPC), which adjusts designs to compensate for lithography-induced distortions after the design phase, FAID directly optimizes for manufacturable designs by using fabrication predictions within the design loop. While OPC can enhance design fidelity, it may not ensure manufacturability for complex features. FAID focuses on creating features that are inherently realizable by the fabrication process, leading to designs that are more practical to manufacture.
For a deeper look at the FAID concept, refer to this paper by Lukas Chrostowski's group at the University of British Columbia.
import autograd as ag
import autograd.numpy as anp
import matplotlib.pylab as plt
import numpy as np
import tidy3d as td
import tidy3d.web as web
import prefab as pf
np.random.seed(2)
Predicting Device Performance¶
Using the prediction model, we can simulate the expected experimental performance of device designs such as the wavelength division multiplexer (WDM) in the Adjoint Optimization of a Wavelength Division Multiplexer. The WDM was optimized with a relatively large feature size and an erosion/dilation penalty to enhance robustness against fabrication variations. This optimization allows it to tolerate predicted fabrication variations quite well; however, some performance degradation is still observed. As illustrated in the figure below, the fabrication-predicted transmission spectra get distorted, insertion loss is increased, and further crosstalk is introduced. This performance degradation can be mitigated with FAID.
¶
¶
Sim Setup¶
Then we set up our basic simulation for the FAID optimization.
We have an input waveguide connected to a square design region, which has n=4 output waveguides.
The square design region is a custom medium with a pixelated permittivity grid that we wish to optimize such that input light of different wavelengths gets directed to different output ports.
As this is a SOI device, we typically define the design region and waveguides as Silicon sitting on a SiO2 substrate. For this demo, we make a 2D simulation, but it can be easily made 3D by changing the Lz parameter, adding dimension to the structures, and adding a substrate.
# material information
n_si = 3.49
n_sio2 = 1.45 # not used in 2D
n_air = 1
# channel wavelengths
wvls_design = np.array([1.270, 1.290, 1.310, 1.330])
freqs_design = td.C_0 / wvls_design
num_freqs_design = len(freqs_design)
freq_max = np.max(freqs_design)
freq_min = np.min(freqs_design)
keys = [str(i) for i in range(num_freqs_design)]
df_design = abs(np.mean(np.diff(freqs_design)))
# forward source
freq0 = np.mean(freqs_design)
wvl0 = td.C_0 / freq0
fwidth = freq_max - freq_min
run_time = 120 / fwidth
# we average the metrics over the channels with some frequency width
channel_fwidth = df_design / 2.0
channel_bounds = [
(f - channel_fwidth / 2, f + channel_fwidth / 2) for f in freqs_design
]
num_freqs_channel = 5
channel_freqs = []
for fmin, fmax in channel_bounds:
sub_freqs = np.linspace(fmin, fmax, num_freqs_channel)
channel_freqs += sub_freqs.tolist()
# size of design region
lx = 4.5
ly = 4.5
ly_single = ly / num_freqs_design
lz = td.inf
# size of waveguides
wg_width = 0.3
wg_length = 1.5
wg_spacing = 0.8
# spacing between design region and PML in y
buffer = 1.5
# size of simulation
Lx = lx + wg_length * 2
Ly = ly + buffer * 2
Lz = 0.0
# fabrication constraints (feature size and projection strength)
radius = 0.100
beta0 = 2
# resolution information
min_steps_per_wvl = 18
dl_design_region = 0.015
Base Simulation¶
First, we'll define the simulation without any design region using the "base" components that don't change over the optimization.
# define the waveguide ports
wg_in = td.Structure(
geometry=td.Box(
center=(-Lx / 2, 0, 0),
size=(wg_length * 2, wg_width, lz),
),
medium=td.Medium(permittivity=n_si**2),
)
centers_y = np.linspace(
-ly / 2.0 + ly_single / 2.0, +ly / 2.0 - ly_single / 2.0, num_freqs_design
)
mode_size = (0, 0.9 * ly_single, td.inf)
wgs_out = []
for center_y in centers_y:
wg_out = td.Structure(
geometry=td.Box(
center=(+Lx / 2, center_y, 0),
size=(wg_length * 2, wg_width, lz),
),
medium=td.Medium(permittivity=n_si**2),
)
wgs_out.append(wg_out)
# measure the mode amplitudes at each of the output ports
mnts_mode = []
for key, center_y in zip(keys, centers_y):
mnt_mode = td.ModeMonitor(
center=(Lx / 2 - wg_length / 2, center_y, 0),
size=mode_size,
freqs=channel_freqs,
mode_spec=td.ModeSpec(),
name=f"mode_{key}",
)
mnts_mode.append(mnt_mode)
# measures the flux at each of the output ports
mnts_flux = []
for key, center_y in zip(keys, centers_y):
mnt_flux = td.FluxMonitor(
center=(Lx / 2 - wg_length / 2, center_y, 0),
size=mode_size,
freqs=channel_freqs,
name=f"flux_{key}",
)
mnts_flux.append(mnt_flux)
# and a field monitor that measures fields on the z=0 plane at the design freqs
fld_mnt = td.FieldMonitor(
center=(0, 0, 0),
size=(td.inf, td.inf, 0),
freqs=freqs_design,
name="field",
)
# inject the fundamental mode into the input waveguide
mode_src = td.ModeSource(
center=(-Lx / 2 + wg_length / 2, 0, 0),
size=mode_size,
source_time=td.GaussianPulse(
freq0=freq0,
fwidth=fwidth,
),
direction="+",
mode_index=0,
)
sim_base = td.Simulation(
size=(Lx, Ly, Lz),
grid_spec=td.GridSpec.auto(
min_steps_per_wvl=min_steps_per_wvl,
wavelength=np.min(wvls_design),
),
structures=[wg_in] + wgs_out,
sources=[mode_src],
monitors=mnts_mode + mnts_flux + [fld_mnt],
boundary_spec=td.BoundarySpec.pml(x=True, y=True, z=True if Lz else False),
run_time=run_time,
)
ax = sim_base.plot_eps(z=0.01)
ax.set_aspect("equal")
Solving modes¶
Next, we want to ensure that we are injecting and measuring the right waveguide modes at each of the ports.
We'll use tidy3d's ModeSolver to analyze the modes of our input waveguide.
from tidy3d.plugins.mode import ModeSolver
from tidy3d.plugins.mode.web import run as run_mode_solver
# we'll ask for 4 modes just to inspect
num_modes = 4
# let's define how large the mode planes are and how far they are from the PML relative to the design region
mode_size = (0, ly_single, td.inf)
# make a plane corresponding to where we wish to measure the input mode
plane_in = td.Box(
center=(-Lx / 2 + wg_length / 2.0, 0, 0),
size=mode_size,
)
mode_solver = ModeSolver(
simulation=sim_base,
plane=plane_in,
freqs=freqs_design,
mode_spec=td.ModeSpec(num_modes=num_modes),
)
Next we run the mode solver on the servers.
mode_data = run_mode_solver(
mode_solver, reduce_simulation=True, results_file="data/mode_solver.hdf5"
)
15:39:02 EST Mode solver created with task_id='fdve-bd79fa16-6098-46c5-aeb7-b41164b4bf2e', solver_id='mo-840d02cc-5271-46a1-9c82-dbd5870a36ce'.
Output()
Output()
15:39:04 EST Mode solver status: success
Output()
And visualize the results.
fig, axs = plt.subplots(num_modes, 3, figsize=(12, 22), tight_layout=True)
for mode_index in range(num_modes):
vmax = 1.1 * max(
abs(mode_data.field_components[n].sel(mode_index=mode_index)).max()
for n in ("Ex", "Ey", "Ez")
)
for field_name, ax in zip(("Ex", "Ey", "Ez"), axs[mode_index]):
for freq in freqs_design:
key = f"{td.C_0 / freq * 1000 :.0f} nm"
field = mode_data.field_components[field_name].sel(
mode_index=mode_index, f=freq
)
field.real.plot(label=f"Real ({key})", ax=ax)
field.imag.plot(ls="--", label=f"Imag ({key})", ax=ax)
ax.set_title(f"index={mode_index}, {field_name}")
ax.set_ylim(-vmax, vmax)
ax.legend()
print("Effective index of computed modes: ", np.array(mode_data.n_eff))
Effective index of computed modes: [[3.1513898 2.840793 2.006425 1.0989405] [3.1436877 2.8191965 1.970435 1.0913303] [3.1359682 2.797195 1.934321 1.0847051] [3.1282325 2.7747855 1.8980899 1.0788841]]