Waveguide Crossing¶
This example uses of PhotonForge to optimize a waveguide crossing based on [1]. A Monte Carlo parameter sweep is used to explore the initial design space and provide a good optimization starting point. The local optimization uses gradients computed automatically through Tidy3D’s adjoint solver, available for parametric components since PhotonForge 1.4.4 (see the Autograd Optimization guide).
References
Sujith Chandran, Marcus Dahlem, Yusheng Bian, Paulo Moreira, Ajey P. Jacob, Michal Rakowski, Andy Stricker, Karen Nummy, Colleen Meagher, Bo Peng, Abu Thomas, Shuren Hu, Jan Petykiewicz, Zoey Sowinski, Won Suk Lee, Rod Augur, Dave Riggs, Ted Letavic, Anthony Yu, Ken Giewont, John Pellerin, and Jaime Viegas, “Beam shaping for ultra-compact waveguide crossings on monolithic silicon photonics platform,” Opt. Lett. 45, 6230-6233 (2020), doi: 10.1364/OL.402446
This example has a similar Tidy3D version: Waveguide crossing based on cosine tapers.
[1]:
import autograd as ag
import autograd.numpy as np # Important to import autograd's numpy module for tracing derivatives
import photonforge as pf
import tidy3d as td
from matplotlib import pyplot as plt
# We raise the Tidy3D logging level to reduce output noise during the
# optimization loop. Warnings can flag real problems, so only suppress them
# once you have confirmed they are not relevant to your case.
td.config.logging.level = "ERROR"
For this demonstration, the basic_technology is enough to show the main concepts behind it. Other technologies can be equally used to optimize the crossing for other PDKs.
[2]:
tech = pf.basic_technology()
pf.config.default_technology = tech
Geometry Parametrization¶
PhotonForge provides a flexible crossing parametric component that can be directly used to re-create the cosine-tapered geometry proposed in the main reference.
We create a custom parametric component function with the parameters we’re interested in optimizing as variables. The decorator is what makes the component differentiable later on: when the function is called with parameters traced by Autograd, it automatically returns a differentiable version of the component (the function body itself always receives plain values, so the expression string for added_width works as usual).
Two details support the gradient computation:
The total arm span is kept fixed by compensating
arm_lengthwithextra_length, so the port positions do not depend on the optimization parameters (gradients are computed for the structures in the simulation, not for the port monitors).The Tidy3D model is set to remove verbosity and increase the default run time because, depending on the exact parameters, the geometry can become resonant.
[3]:
total_arm = 3.5
@pf.parametric_component
def sine_crossing(*, amplitude, arm_length):
model = pf.Tidy3DModel(verbose=False, run_time=3e-12)
crossing = pf.parametric.crossing(
port_spec="Strip",
arm_length=arm_length,
extra_length=total_arm - arm_length,
added_width=f"{amplitude} * sin(pi * u)",
model=model,
)
return crossing
sine_crossing(arm_length=3, amplitude=1)
[3]:
Monte Carlo Exploration¶
The initial exploration of the parameter space is done with the help of the Monte Carlo tools in PhotonForge. We can calculate a range of S parameters for variations of the crossing by setting its parameters as random variables within a fixed value range. The ranges are quite wide because we know almost nothing a priori about the design, therefore the number of samples is also considerably high to allow proper exploration of the whole space.
Note that because we are only interested in the transmission from P0 to P2, we explicitly set the ports and modes we want to run as sources for the Tidy3D simulations using the inputs argument.
[4]:
freqs = pf.C_0 / np.linspace(1.45, 1.65, 21)
variables, results = pf.monte_carlo.s_matrix(
sine_crossing,
freqs,
("amplitude", "component", {"value_range": (0.1, 1.5)}),
("arm_length", "component", {"value_range": (0.5, 3)}),
random_samples=30,
component_kwargs={"amplitude": 1, "arm_length": 3},
model_kwargs={"inputs": ["P0@0"]},
random_seed=0,
)
The results include the full S matrix for each crossing variation. To help us make sense of all this data, we can define a figure of merit for the design and plot it versus the crossing parameters. In this case, a reasonable definition is to use the minimal transmission value within the frequency range we’re interested in. Let’s first create a result table with only the parameters and our figure of merit, instead of the whole S matrix objects:
[5]:
names = [v.name for v in variables]
table = np.array([(*v, np.abs(s["P0@0", "P2@0"]).min() ** 2) for *v, s in results])
table
[5]:
array([[1.89293245, 1.30670413, 0.05249639],
[1.78530457, 0.95793781, 0.33610302],
[1.57431442, 0.15939063, 0.73723963],
[1.39931042, 1.10239893, 0.15817887],
[1.74379972, 1.32907849, 0.0217574 ],
[0.66835153, 0.56466724, 0.43499857],
[2.07250101, 1.05543504, 0.49590571],
[2.45512328, 0.64501438, 0.84326877],
[2.55332675, 0.68507742, 0.87322951],
[2.12682069, 0.21543469, 0.71566396],
[0.98118614, 0.76343902, 0.35722034],
[2.66823248, 0.38698483, 0.71412424],
[1.46807783, 1.44878793, 0.21811395],
[1.94662566, 0.60877627, 0.93363979],
[2.92168546, 0.83047987, 0.85789326],
[2.35086521, 0.36836535, 0.71803131],
[0.86364256, 0.28845985, 0.70117002],
[2.24903115, 1.48379688, 0.02436091],
[2.26344246, 1.15381968, 0.53983602],
[2.87219256, 1.0157638 , 0.59663377],
[0.59253126, 0.90813608, 0.43886585],
[1.26638186, 0.50425125, 0.67516803],
[1.00809177, 0.73810412, 0.34653502],
[0.58137162, 1.24378247, 0.49406468],
[2.80040379, 1.17426609, 0.55669658],
[1.10690058, 0.44344074, 0.6420322 ],
[1.17656064, 0.88775554, 0.28933405],
[0.75534563, 0.12637548, 0.68833767],
[2.60544912, 1.4027804 , 0.36152831],
[1.66009647, 0.24387858, 0.74642182]])
[6]:
plt.scatter(table[:, 0], table[:, 1], c=table[:, 2])
plt.xlabel(names[0])
plt.ylabel(names[1])
_ = plt.colorbar()
Local Optimization¶
Around the best result from the exploration stage, we switch to gradient-based optimization: because the parameters are traced by Autograd, the s_matrix call returns S parameters that gradients can flow through, computed by PhotonForge via Tidy3D’s adjoint solver. Each gradient evaluation costs one forward and one adjoint simulation, independently of the number of parameters — see the Autograd Optimization guide for details.
While the exploration stage ranks designs by their minimal transmission over the band, for the local refinement we use the transmission at the central frequency only: the crossing response is smooth over this range, single-frequency gradients require a single adjoint simulation, and the result files stay small because only one frequency of volumetric field data is involved. The full band is verified after the optimization. We also restrict the simulation inputs to the only excited port, as
before.
[7]:
freq0 = pf.C_0 / 1.55
def objective(params):
# Map by name: the variable order in the Monte Carlo results table is not
# guaranteed to match the order they were declared in
kwargs = dict(zip(names, params))
c = sine_crossing(**kwargs)
s_matrix = c.s_matrix([freq0], show_progress=False, model_kwargs={"inputs": ["P0@0"]})
transmission = np.abs(s_matrix["P0@0", "P2@0"]) ** 2
return transmission[0]
For the actual optimization run, we use a few iterations of the Adam algorithm, clipping the parameters to their exploration ranges after each update.
Note that, near convergence, the objective fluctuates by a small amount due to re-meshing noise from sub-grid geometry changes, so there is no benefit in running many iterations with very small updates: if needed, the optimization can always be restarted from the last best parameters with a finer simulation grid.
[8]:
best_variation = np.argmax(table[:, 2])
params = np.array([table[best_variation, 0], table[best_variation, 1]])
bounds = tuple(np.array([v.value_spec["value_range"] for v in variables]).T)
learning_rate = 0.05
value_and_grad = ag.value_and_grad(objective)
m = np.zeros_like(params)
v = np.zeros_like(params)
history = []
for i in range(1, 9):
value, gradient = value_and_grad(params)
history.append((params, value))
print(f"Iteration {i}: transmission = {value:.3f} at {np.round(params, 3)}")
# Adam update (gradient ascent)
m = 0.9 * m + 0.1 * gradient
v = 0.999 * v + 0.001 * gradient**2
step = (m / (1 - 0.9**i)) / (np.sqrt(v / (1 - 0.999**i)) + 1e-8)
params = np.clip(params + learning_rate * step, *bounds)
best_params, best_value = max(history, key=lambda h: h[1])
print(f"Best: {best_value:.3f} at {np.round(best_params, 3)}")
Iteration 1: transmission = 0.943 at [1.947 0.609]
Iteration 2: transmission = 0.942 at [1.997 0.559]
Iteration 3: transmission = 0.940 at [1.978 0.534]
Iteration 4: transmission = 0.948 at [1.945 0.559]
Iteration 5: transmission = 0.948 at [1.913 0.58 ]
Iteration 6: transmission = 0.948 at [1.894 0.58 ]
Iteration 7: transmission = 0.951 at [1.887 0.567]
Iteration 8: transmission = 0.951 at [1.885 0.55 ]
Best: 0.951 at [1.887 0.567]
Inspecting the Optimized Component¶
We can recreate the best parameters found in the optimization history for our technology and reset its model to remove the helper symmetries and include a field monitor if we want to fully inspect the component. Of course, with our parameters being layout dimensions, we should not use an arbitrary number of decimals, so we round them appropriately: in this case, 2 decimals (10 nm grid) should suffice.
[9]:
kwargs = dict(zip(names, np.round(best_params, decimals=2)))
crossing = sine_crossing(**kwargs)
tidy3d_model = pf.Tidy3DModel(
monitors=[
td.FieldMonitor(
name="field",
center=(0, 0, 0.11),
size=(td.inf, td.inf, 0),
freqs=[freqs.mean()],
)
]
)
crossing.add_model(tidy3d_model, "Tidy3D")
s_matrix = crossing.s_matrix(freqs)
fig, ax = pf.plot_s_matrix(s_matrix)
_ = fig.suptitle(
f"Amplitude: {kwargs['amplitude']} μm, arm length: {kwargs['arm_length']} μm"
)
14:39:08 Eastern Daylight Time Loading simulation from local cache. View cached task using web UI at 'https://tidy3d.simulation.cloud/workbench?taskId =fdve-97b4903f-f702-410d-9e9a-78fb0c2c5689'.
Loading simulation from local cache. View cached task using web UI at 'https://tidy3d.simulation.cloud/workbench?taskId =fdve-a3789ac2-2bf6-4efe-8f00-616fe97b7334'.
Loading simulation from local cache. View cached task using web UI at 'https://tidy3d.simulation.cloud/workbench?taskId =fdve-a8fe831d-5e58-4b15-b6da-733391279bda'.
14:39:09 Eastern Daylight Time Loading simulation from local cache. View cached task using web UI at 'https://tidy3d.simulation.cloud/workbench?taskId =fdve-c3e35e5d-a07a-4f38-bdf0-8ea0e397fdf5'.
In order to plot the fields, we can use the batch_data function from the Tidy3D model to get the BatchData that includes the fields from our extra monitor:
[10]:
sim_data = tidy3d_model.batch_data(crossing, freqs)
_ = sim_data["P0@0"].plot_field("field", "Hz", robust=False)
Progress: 100%
Loading simulation from local cache. View cached task using web UI at 'https://tidy3d.simulation.cloud/workbench?taskId =fdve-97b4903f-f702-410d-9e9a-78fb0c2c5689'.
Loading simulation from local cache. View cached task using web UI at 'https://tidy3d.simulation.cloud/workbench?taskId =fdve-a3789ac2-2bf6-4efe-8f00-616fe97b7334'.
Loading simulation from local cache. View cached task using web UI at 'https://tidy3d.simulation.cloud/workbench?taskId =fdve-a8fe831d-5e58-4b15-b6da-733391279bda'.
Loading simulation from local cache. View cached task using web UI at 'https://tidy3d.simulation.cloud/workbench?taskId =fdve-c3e35e5d-a07a-4f38-bdf0-8ea0e397fdf5'.
