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

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

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.

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_length with extra_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.

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)

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.

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,
)

Starting sample 1 of 30…
Starting sample 2 of 30…
Starting sample 3 of 30…
Starting sample 4 of 30…
Starting sample 5 of 30…
Starting sample 6 of 30…
Starting sample 7 of 30…
Starting sample 8 of 30…
Starting sample 9 of 30…
Starting sample 10 of 30…
Starting sample 11 of 30…
Starting sample 12 of 30…
Starting sample 13 of 30…
Starting sample 14 of 30…
Starting sample 15 of 30…
Starting sample 16 of 30…
Starting sample 17 of 30…
Starting sample 18 of 30…
Starting sample 19 of 30…
Starting sample 20 of 30…
Starting sample 21 of 30…
Starting sample 22 of 30…
Starting sample 23 of 30…
Starting sample 24 of 30…
Starting sample 25 of 30…
Starting sample 26 of 30…
Starting sample 27 of 30…
Starting sample 28 of 30…
Starting sample 29 of 30…
Starting sample 30 of 30…
Sample 1 done.
Sample 2 done.
Sample 3 done.
Sample 4 done.
Sample 5 done.
Sample 6 done.
Sample 7 done.
Sample 8 done.
Sample 9 done.
Sample 10 done.
Sample 12 done.
Sample 13 done.
Sample 11 done.
Sample 15 done.
Sample 14 done.
Sample 16 done.
Sample 17 done.
Sample 18 done.
Sample 22 done.
Sample 19 done.
Sample 20 done.
Sample 21 done.
Sample 23 done.
Sample 24 done.
Sample 25 done.
Sample 28 done.
Sample 26 done.
Sample 27 done.
Sample 29 done.
Sample 30 done.
All samples done!

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:

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

array([[1.89293245, 1.30670413, 0.05251185],
       [1.78530457, 0.95793781, 0.33611872],
       [1.57431442, 0.15939063, 0.73721869],
       [1.39931042, 1.10239893, 0.15817226],
       [1.74379972, 1.32907849, 0.0217593 ],
       [0.66835153, 0.56466724, 0.43496782],
       [2.07250101, 1.05543504, 0.49591299],
       [2.45512328, 0.64501438, 0.84323516],
       [2.55332675, 0.68507742, 0.87319082],
       [2.12682069, 0.21543469, 0.71566228],
       [0.98118614, 0.76343902, 0.35722436],
       [2.66823248, 0.38698483, 0.71411322],
       [1.46807783, 1.44878793, 0.21810369],
       [1.94662566, 0.60877627, 0.93366074],
       [2.92168546, 0.83047987, 0.85787784],
       [2.35086521, 0.36836535, 0.71801924],
       [0.86364256, 0.28845985, 0.70114407],
       [2.24903115, 1.48379688, 0.02434806],
       [2.26344246, 1.15381968, 0.5398749 ],
       [2.87219256, 1.0157638 , 0.59674977],
       [0.59253126, 0.90813608, 0.43885717],
       [1.26638186, 0.50425125, 0.67519111],
       [1.00809177, 0.73810412, 0.34646837],
       [0.58137162, 1.24378247, 0.49403132],
       [2.80040379, 1.17426609, 0.55671789],
       [1.10690058, 0.44344074, 0.64204262],
       [1.17656064, 0.88775554, 0.28930014],
       [0.75534563, 0.12637548, 0.68830613],
       [2.60544912, 1.4027804 , 0.36149046],
       [1.66009647, 0.24387858, 0.74640803]])

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.

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.

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.939 at [1.978 0.534]
Iteration 4: transmission = 0.948 at [1.945 0.559]
Iteration 5: transmission = 0.948 at [1.913 0.582]
Iteration 6: transmission = 0.946 at [1.893 0.586]
Iteration 7: transmission = 0.949 at [1.889 0.575]
Iteration 8: transmission = 0.951 at [1.892 0.557]
Best: 0.951 at [1.892 0.557]

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.

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"
)

Uploading task 'P0@0…'
Uploading task 'P1@0…'
Uploading task 'P2@0…'
Uploading task 'P3@0…'
Starting task 'P0@0': https://tidy3d.simulation.cloud/workbench?taskId=fdve-fda5265f-4c31-4265-96c9-3017b3b71c12
Starting task 'P1@0': https://tidy3d.simulation.cloud/workbench?taskId=fdve-c0a7bc91-e303-49f5-a1b0-98000f0ba423
Starting task 'P2@0': https://tidy3d.simulation.cloud/workbench?taskId=fdve-3f59ccea-1445-4d9d-b970-46125170ca6b
Downloading data from 'P0@0'…
Downloading data from 'P1@0'…
Downloading data from 'P2@0'…
Starting task 'P3@0': https://tidy3d.simulation.cloud/workbench?taskId=fdve-9e5fd6b8-2578-438b-8610-94bf5f43c68f
Downloading data from 'P3@0'…
Progress: 100%

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:

sim_data = tidy3d_model.batch_data(crossing, freqs)
_ = sim_data["P0@0"].plot_field("field", "Hz", robust=False)

Progress: 100%

07:34:27 -03 Loading simulation from local cache. View cached task using web UI
             at
             'https://tidy3d.simulation.cloud/workbench?taskId=fdve-fda5265f-4c3
             1-4265-96c9-3017b3b71c12'.

             Loading simulation from local cache. View cached task using web UI
             at
             'https://tidy3d.simulation.cloud/workbench?taskId=fdve-c0a7bc91-e30
             3-49f5-a1b0-98000f0ba423'.

             Loading simulation from local cache. View cached task using web UI
             at
             'https://tidy3d.simulation.cloud/workbench?taskId=fdve-3f59ccea-144
             5-4d9d-b970-46125170ca6b'.

             Loading simulation from local cache. View cached task using web UI
             at
             'https://tidy3d.simulation.cloud/workbench?taskId=fdve-9e5fd6b8-257
             8-438b-8610-94bf5f43c68f'.