Digital 1×2 splitter: shallow-hole inverse design

This notebook reconstructs Tidy3D Autograd24DigitalSplitter.ipynb with FDTDX. It uses the same 20×10 independent shallow-hole permittivities, uniform ε=6.5 start, five wavelengths, electric y symmetry, physical upper-port FOM, and two published 25-update stages. If the binary device has not converged after those stages, the same update continues within the benchmark’s device-independent 3× iteration cap. It does not use optimized Tidy3D pixels or a handcrafted warm start.

The retained optimization history comes from the benchmark artifact. This notebook then performs a fresh five-frequency FDTDX field solve of the fully digitized radius-sweep winner, so the fabricated geometry, Yee material, fields, and spectrum are inspectable.

Imports and accelerator

import os
os.environ.setdefault("XLA_PYTHON_CLIENT_PREALLOCATE", "false")
import json
from pathlib import Path
from time import perf_counter
import warnings
import matplotlib.pyplot as plt
from matplotlib.patches import Circle, Rectangle
import numpy as np
import pandas as pd
import jax
import yaml
from IPython import get_ipython

from benchmarks.invdes_digital_splitter_runner import (
    DESIGN_LENGTH_UM, DL_UM, DL_Z_UM, EPS_AIR, EPS_SILICON,
    HOLE_DEPTH_UM, HOLE_RADIUS_UM, PARAMETER_SHAPE, PIXEL_UM,
    SILICON_HEIGHT_UM, WAVELENGTHS_UM, initial_parameters,
    optimize_digital_splitter, simulate_digital_splitter,
)

get_ipython().run_line_magic("matplotlib", "inline")
warnings.filterwarnings("ignore", category=FutureWarning, module=r"tidy3d\.components\.mode")
plt.rcParams.update({"figure.figsize": (10, 4.8), "figure.dpi": 120})
meta = yaml.safe_load(Path("benchmarks/cases/invdes_power_splitter/case.yaml").read_text())
print(f"JAX devices: {jax.devices()}")

JAX devices: [CudaDevice(id=0)]

1. Published physical and numerical parameters

pd.Series({
    "wavelengths_um": tuple(WAVELENGTHS_UM),
    "design_region_um": (DESIGN_LENGTH_UM, DESIGN_LENGTH_UM),
    "independent_holes": PARAMETER_SHAPE,
    "pitch_um": PIXEL_UM,
    "hole_radius_um": HOLE_RADIUS_UM,
    "hole_depth_um": HOLE_DEPTH_UM,
    "silicon_height_um": SILICON_HEIGHT_UM,
    "FDTD_dx_dy_um": DL_UM,
    "FDTD_dz_um": DL_Z_UM,
    "symmetry": (0, -1, 0),
    "published_updates": 50,
    "automatic_iteration_cap": meta["invdes"]["max_iterations"],
}, name="value").to_frame()

	value
wavelengths_um	(1.53, 1.54, 1.55, 1.56, 1.57)
design_region_um	(2.6, 2.6)
independent_holes	(20, 10)
pitch_um	0.13
hole_radius_um	0.045
hole_depth_um	0.14
silicon_height_um	0.22
FDTD_dx_dy_um	0.045
FDTD_dz_um	0.05
symmetry	(0, -1, 0)
published_updates	50
automatic_iteration_cap	150

2. Exact hole geometry and honest uniform start

parameters0 = initial_parameters()
x_centers = -DESIGN_LENGTH_UM / 2 + (np.arange(PARAMETER_SHAPE[0]) + 0.5) * PIXEL_UM
y_centers = (np.arange(PARAMETER_SHAPE[1]) + 0.5) * PIXEL_UM

fig, axes = plt.subplots(1, 2, figsize=(12, 5))
ax = axes[0]
ax.add_patch(Rectangle((-DESIGN_LENGTH_UM / 2, -DESIGN_LENGTH_UM / 2),
                       DESIGN_LENGTH_UM, DESIGN_LENGTH_UM,
                       color="#55c2b5", alpha=.45, label="220 nm silicon slab"))
for x in x_centers:
    for y in y_centers:
        for signed_y in (-y, y):
            ax.add_patch(Circle((x, signed_y), HOLE_RADIUS_UM,
                                facecolor="#8221a8", edgecolor="none", alpha=.72))
ax.set(xlim=(-1.5, 1.5), ylim=(-1.5, 1.5), aspect="equal",
       xlabel="x (µm)", ylabel="y (µm)", title="200 independent holes + y mirror")
ax.legend(loc="upper right")

image = axes[1].imshow(parameters0.T, origin="lower", cmap="viridis",
                       extent=(-1.3, 1.3, 0, 1.3), vmin=EPS_AIR, vmax=EPS_SILICON,
                       aspect="auto")
axes[1].set(xlabel="x (µm)", ylabel="retained y (µm)",
            title="Published uniform ε=6.5 start")
fig.colorbar(image, ax=axes[1], label="hole permittivity")
plt.tight_layout(); plt.show()

3. Optimization and symmetry-normalized objective

optimize_digital_splitter executes the full differentiable 3D FDTD loop. Set RUN_FULL_OPTIMIZATION=True to repeat the cold-start adjoints. The retained notebook defaults to the audited benchmark result to avoid paying that GPU cost twice; the best design receives a fresh forward solve below.

The physical upper-port power is half the reduced output/source ratio because the source plane is clipped by y symmetry and the upper-output plane is not.

RUN_FULL_OPTIMIZATION = False
if RUN_FULL_OPTIMIZATION:
    result = optimize_digital_splitter(meta)
else:
    progress = json.loads(Path("progress.json").read_text())
    result = progress["cases"]["invdes_power_splitter"]
assert result["status"] in {"pass", "fail"}, result
print("Scientific result:", result["status"].upper())
display(pd.DataFrame(result["metrics"]))

Scientific result: FAIL

	name	compare	got	ref	pass	error	atol	rtol	max_got	max_ref
0	best_fom_parity	fom_ge	0.432533	0.448495	False	0.015962	0.0	0.0	NaN	NaN
1	optimizer_improvement	fom_ge	0.324736	0.000100	True	-0.324636	0.0	0.0	NaN	NaN
2	passivity_total_transmission	upper_bound	NaN	NaN	True	0.000000	0.0	0.0	0.924030	1.05
3	source_normalization_drift	upper_bound	NaN	NaN	True	0.000000	0.0	0.0	0.018771	0.05
4	final_binary_fraction	fom_ge	1.000000	1.000000	True	0.000000	0.0	0.0	NaN	NaN

4. Convergence, passivity, and source loading

extras = result["extras"]
history = np.asarray(extras["fom_history"])
source_drift = np.max(np.abs(np.asarray(extras["source_power_ratios"]) - 1), axis=1)
total_peak = np.max(np.asarray(extras["total_transmission_spectra"]), axis=1)
golden = np.load("benchmarks/goldens/invdes_power_splitter.npz", allow_pickle=True)

fig, axes = plt.subplots(1, 2, figsize=(12, 4.5))
axes[0].plot(history, "o-", ms=3, label="FDTDX physical upper port")
axes[0].plot(np.asarray(golden["fom_history"]), "--", lw=1.4, label="Tidy3D notebook")
if "checkpoint_iterations" in extras:
    axes[0].scatter(extras["checkpoint_iterations"], extras["checkpoint_binary_foms"],
                    marker="s", color="#087d83", zorder=4,
                    label="fresh fully binary checkpoints")
axes[0].axhline(extras["threshold"], color="#d29b00", ls=":", label="95% threshold")
axes[0].set(xlabel="update", ylabel="FOM", title="Same one-output objective")
axes[0].legend(); axes[0].grid(alpha=.25)

axes[1].plot(total_peak, label="maximum total transmission")
axes[1].plot(source_drift, label="source-loading change")
axes[1].axhline(1.05, color="red", ls=":", label="passivity limit")
axes[1].set(xlabel="update", ylabel="fraction", title="Independent numerical audits")
axes[1].legend(); axes[1].grid(alpha=.25)
plt.tight_layout(); plt.show()

5. Continuous controls and fully binary fabricated design

artifact = np.load(extras["best_design_artifact"], allow_pickle=True)
best = np.asarray(
    artifact["selected_checkpoint_hole_permittivities"]
    if "selected_checkpoint_hole_permittivities" in artifact
    else artifact["hole_permittivities"]
)
binary_epsilon = np.asarray(artifact["binary_hole_permittivities"])
binary_radii = np.asarray(artifact["binary_hole_radii_um"])
fig, axes = plt.subplots(1, 2, figsize=(12, 5))
image = axes[0].imshow(best.T, origin="lower", cmap="viridis", extent=(-1.3, 1.3, 0, 1.3),
                       vmin=EPS_AIR, vmax=EPS_SILICON, aspect="auto")
axes[0].set(xlabel="x (µm)", ylabel="retained y (µm)", title="Controls at selected binary checkpoint")
fig.colorbar(image, ax=axes[0], label="hole permittivity")

ax = axes[1]
ax.add_patch(Rectangle((-1.3, -1.3), 2.6, 2.6, color="#55c2b5", alpha=.35))
for i, x in enumerate(x_centers):
    for j, y in enumerate(y_centers):
        for signed_y in (-y, y):
            radius = binary_radii[i, j]
            if radius > 0:
                ax.add_patch(Circle((x, signed_y), radius,
                                    facecolor="white", edgecolor="#18364a", lw=.25))
ax.set(xlim=(-1.4, 1.4), ylim=(-1.4, 1.4), aspect="equal",
       xlabel="x (µm)", ylabel="y (µm)", title="Accepted binary Si/air radii")
plt.tight_layout(); plt.show()
print(f"Binary material fraction: {extras['final_binary_fraction']:.6f}")
print(f"Selected reduced radius: {extras['best_reduced_radius_um'] * 1e3:.1f} nm")

Binary material fraction: 1.000000
Selected reduced radius: 30.0 nm

6. Fresh five-frequency FDTDX validation solve

started = perf_counter()
validation = simulate_digital_splitter(
    binary_epsilon, hole_radii_um=binary_radii, record_fields=True
)
print(f"Fresh forward solve: {perf_counter() - started:.1f} s")
print(f"Reduced grid: {validation['grid']}; steps: {validation['steps']:,}")
pd.DataFrame({
    "wavelength_um": validation["wavelengths_um"],
    "upper_port_power": validation["upper_port_power"],
    "total_two_port_transmission": validation["total_transmission"],
})

Fresh forward solve: 30.3 s
Reduced grid: (124, 62, 67); steps: 12,161

	wavelength_um	upper_port_power	total_two_port_transmission
0	1.53	0.417988	0.835975
1	1.54	0.443194	0.886387
2	1.55	0.456906	0.913811
3	1.56	0.447226	0.894452
4	1.57	0.397354	0.794708

7. Yee-grid material and fields at every design wavelength

extent = validation["extent_um"]
epsilon = validation["permittivity_xy"]
fields = validation["field_intensity_xy"]
x = np.linspace(extent[0], extent[1], epsilon.shape[0])
y = np.linspace(extent[2], extent[3], epsilon.shape[1])

fig, ax = plt.subplots(figsize=(9, 6))
image = ax.imshow(epsilon.T, origin="lower", extent=extent, cmap="viridis", aspect="equal")
ax.set(xlabel="x (µm)", ylabel="y (µm)", title="Fully binary design on the actual Yee grid")
fig.colorbar(image, ax=ax, label="relative permittivity"); plt.show()

fig, axes = plt.subplots(2, 3, figsize=(15, 9), sharex=True, sharey=True)
for ax, wavelength, intensity in zip(axes.flat, WAVELENGTHS_UM, fields, strict=False):
    normalized = intensity / (intensity.max() + 1e-30)
    image = ax.imshow(normalized.T, origin="lower", extent=extent, cmap="magma",
                      vmin=0, vmax=1, aspect="equal")
    ax.contour(x, y, epsilon.T, levels=[3.0], colors="cyan", linewidths=.35)
    ax.set_title(f"λ = {wavelength:.2f} µm")
axes.flat[-1].axis("off")
for ax in axes[-1, :2]: ax.set_xlabel("x (µm)")
for ax in axes[:, 0]: ax.set_ylabel("y (µm)")
fig.colorbar(image, ax=axes, label="normalized |E|²", shrink=.72)
fig.suptitle("Fresh binary-design phasor fields—not terminal-time snapshots")
plt.show()

8. Spectral result and insertion loss

wavelengths = np.asarray(validation["wavelengths_um"])
upper = np.asarray(validation["upper_port_power"])
total = np.asarray(validation["total_transmission"])
fig, axes = plt.subplots(1, 2, figsize=(12, 4.5))
axes[0].plot(wavelengths, upper, "o-", label="upper output")
axes[0].plot(wavelengths, total, "s-", label="two-output total")
axes[0].axhline(.5, color="gray", ls=":", label="one-port upper bound")
axes[0].set(xlabel="wavelength (µm)", ylabel="modal power", title="Fresh best-design spectrum")
axes[0].legend(); axes[0].grid(alpha=.25)
axes[1].plot(wavelengths, 10 * np.log10(np.maximum(total, 1e-12)), "o-", color="#8221a8")
axes[1].set(xlabel="wavelength (µm)", ylabel="insertion loss (dB)",
            title="Same 2× convention as the Tidy3D notebook")
axes[1].grid(alpha=.25)
plt.tight_layout(); plt.show()

fig, ax = plt.subplots(figsize=(6.5, 4.2))
ax.plot(np.asarray(extras["digitized_radius_sweep_um"]) * 1e3,
        extras["digitized_fom_sweep"], "o-", color="#087d83")
ax.axhline(extras["ref_fom"], color="#8221a8", ls="--", label="pinned Tidy3D")
ax.axhline(extras["threshold"], color="#d29b00", ls=":", label="95% gate")
ax.set(xlabel="reduced-hole radius (nm)", ylabel="fresh binary FOM",
       title="Published fabrication-radius sweep")
ax.legend(); ax.grid(alpha=.25); plt.show()

Reproduce and refine

Run uv run fdtdx-bench run --case invdes_power_splitter to repeat the optimizer. The 45 nm target mesh is GPU-memory adapted and automatically aligned to the substrate, etch-bottom, and slab-top interfaces. Production sign-off still needs refinement toward the notebook’s ~22.4 nm silicon mesh and wider process-corner validation after digitization.

print("uv run fdtdx-bench run --case invdes_power_splitter")

uv run fdtdx-bench run --case invdes_power_splitter