"""Broadband FDTD/EME transmission comparison for an X-cut TFLN Euler bend. The model uses a literature-typical TFLN rib stack and compares an explicit curved FDTD bend against the Flexcompute anisotropic-bend EME construction. The final EME setup is N=10 longitudinal bend cells and M=16 modes per cell. All dimensions are in micrometers and frequencies are in Hz. The analysis reuses matching X-cut local result files when present. Commands: python tfln_euler_bend_broadband_final.py build python tfln_euler_bend_broadband_final.py estimate python tfln_euler_bend_broadband_final.py run python tfln_euler_bend_broadband_final.py analyze """ from __future__ import annotations import argparse from concurrent.futures import ThreadPoolExecutor from dataclasses import dataclass import json from pathlib import Path import gdstk import matplotlib matplotlib.use("Agg") import matplotlib.pyplot as plt import numpy as np import tidy3d as td from scipy.integrate import cumulative_trapezoid from scipy.interpolate import CubicSpline # ----------------------------------------------------------------------------- # Device and numerical parameters. # ----------------------------------------------------------------------------- WAVELENGTHS_UM = np.linspace(1.52, 1.58, 6) FREQUENCIES_HZ = np.sort(td.C_0 / WAVELENGTHS_UM) GRID_WAVELENGTH_UM = float(WAVELENGTHS_UM.min()) CENTER_WAVELENGTH_UM = 1.55 CENTER_FREQUENCY_HZ = td.C_0 / CENTER_WAVELENGTH_UM # Typical thin-film lithium-niobate rib geometry. LN_THICKNESS_UM = 0.60 ETCH_DEPTH_UM = 0.30 SLAB_THICKNESS_UM = LN_THICKNESS_UM - ETCH_DEPTH_UM RIDGE_TOP_WIDTH_UM = 0.90 SIDEWALL_FROM_HORIZONTAL_DEG = 70.0 SIDEWALL_ANGLE_RAD = np.deg2rad(90.0 - SIDEWALL_FROM_HORIZONTAL_DEG) # 90-degree symmetric Euler bend. MIN_BEND_RADIUS_UM = 10.0 BEND_ANGLE_RAD = np.pi / 2 INPUT_LEAD_UM = 12.0 OUTPUT_LEAD_UM = 12.0 # Final EME resolution after mode-count, bend-cell, and mode-plane checks. EME_NUM_CELLS = 10 EME_NUM_MODES = 16 # Mode planes and mesh. MODE_WINDOW_Y_UM = 10.0 MODE_WINDOW_Z_UM = 8.0 # This window is large enough for the guided and relevant radiation content, # while avoiding unstable high-order radiation modes seen with a 12 x 9 um # EME plane. It also matches the radius-sweep setup. EME_PLANE_Y_UM = 9.0 EME_PLANE_Z_UM = 6.0 TARGET_NEFF = 2.00 STEPS_PER_WAVELENGTH = 40 PATH_TOLERANCE_UM = 1e-4 PATH_SAMPLES = 4001 TASK_PREFIX = "tfln_euler_xcut_r10_tefund_final_broadband_1520_1580_m16" OUTPUT_DIR = Path("results") MANIFEST_PATH = OUTPUT_DIR / f"{TASK_PREFIX}_tasks.json" @dataclass(frozen=True) class EulerCenterline: s: np.ndarray x: np.ndarray y: np.ndarray curvature: np.ndarray @property def length(self) -> float: return float(self.s[-1]) @property def x_end(self) -> float: return float(self.x[-1]) @property def y_end(self) -> float: return float(self.y[-1]) def build_euler_centerline() -> EulerCenterline: """Construct a symmetric 0-to-maximum-to-0 curvature Euler turn.""" half_length = BEND_ANGLE_RAD * MIN_BEND_RADIUS_UM s = np.linspace(0.0, 2.0 * half_length, PATH_SAMPLES) curvature = np.where( s <= half_length, s / (MIN_BEND_RADIUS_UM * half_length), (2.0 * half_length - s) / (MIN_BEND_RADIUS_UM * half_length), ) tangent_angle = cumulative_trapezoid(curvature, s, initial=0.0) x = cumulative_trapezoid(np.cos(tangent_angle), s, initial=0.0) y = cumulative_trapezoid(np.sin(tangent_angle), s, initial=0.0) return EulerCenterline(s=s, x=x, y=y, curvature=curvature) def ln_medium(): """X-cut LiNbO3 with its optic axis along global y in the bend plane.""" return td.material_library["LiNbO3"]["Zelmon1997"](1) def sio2_medium(): try: return td.material_library["SiO2"]["Palik_NoLoss"] except KeyError: return td.material_library["SiO2"]["Palik_Lossless"] def mode_sort_spec(*, track_freq: str | None = None): """Select the highest-neff quasi-TE mode at every frequency.""" return td.ModeSortSpec( filter_key="TE_fraction", filter_reference=0.5, filter_order="over", sort_key="n_eff", sort_order="descending", track_freq=track_freq, ) def mode_spec( *, bend_radius: float | None = None, for_eme: bool = False, ) -> td.ModeSpec: cls = td.EMEModeSpec if for_eme else td.ModeSpec kwargs = dict( num_modes=EME_NUM_MODES if for_eme else 2, target_neff=TARGET_NEFF, num_pml=(12, 12), bend_radius=bend_radius, bend_axis=1 if bend_radius is not None else None, # Track the full EME basis from the central frequency so broadband # high-order radiation-mode crossings do not create spectral kinks. sort_spec=mode_sort_spec(track_freq="central" if for_eme else None), precision="double", ) if for_eme: # Keep the crystal tensor fixed in the laboratory frame as the local # propagation direction rotates, matching both FDTD and the official # AnisotropicBendsEME notebook. kwargs["bend_medium_frame"] = "global" kwargs["increasing_mode_tolerance"] = 1e-3 # Solve all six requested frequencies explicitly; mode tracking keeps # the broadband basis consistent without spectral interpolation. kwargs["interp_spec"] = None return cls(**kwargs) def stack_structures(ridge_geometry: td.Geometry) -> list[td.Structure]: medium = ln_medium() slab = td.Structure( geometry=td.Box( center=(0, 0, -SLAB_THICKNESS_UM / 2), size=(td.inf, td.inf, SLAB_THICKNESS_UM), ), medium=medium, ) ridge = td.Structure(geometry=ridge_geometry, medium=medium) return [slab, ridge] def ridge_polygon(centerline: EulerCenterline) -> np.ndarray: """Create one smooth ridge polygon along the Euler path and leads.""" x_spline = CubicSpline(centerline.s, centerline.x) y_spline = CubicSpline(centerline.s, centerline.y) def curve(u: float) -> tuple[float, float]: s = float(np.clip(u, 0.0, 1.0) * centerline.length) return float(x_spline(s)), float(y_spline(s)) def gradient(u: float) -> tuple[float, float]: s = float(np.clip(u, 0.0, 1.0) * centerline.length) return ( float(x_spline(s, 1) * centerline.length), float(y_spline(s, 1) * centerline.length), ) cell = gdstk.Cell("tfln_final_euler_ridge") path = gdstk.RobustPath( (-INPUT_LEAD_UM, 0.0), RIDGE_TOP_WIDTH_UM, tolerance=PATH_TOLERANCE_UM, max_evals=100000, ends="flush", ) path.segment((0.0, 0.0)) path.parametric( curve, path_gradient=gradient, width=RIDGE_TOP_WIDTH_UM, relative=False, ) path.segment((centerline.x_end, centerline.y_end + OUTPUT_LEAD_UM)) cell.add(path) polygons = cell.get_polygons() if len(polygons) != 1: raise RuntimeError(f"Expected one ridge polygon, got {len(polygons)}.") return np.asarray(polygons[0].points, dtype=float) def ridge_polyslab(vertices: np.ndarray) -> td.PolySlab: return td.PolySlab( axis=2, slab_bounds=(0.0, LN_THICKNESS_UM), vertices=vertices, sidewall_angle=SIDEWALL_ANGLE_RAD, reference_plane="top", ) def build_fdtd(centerline: EulerCenterline) -> td.Simulation: polygon = ridge_polygon(centerline) structures = stack_structures(ridge_polyslab(polygon)) input_x = -3.0 output_y = centerline.y_end + 3.0 input_plane = td.Box( center=(input_x, 0.0, LN_THICKNESS_UM / 2), size=(0.0, MODE_WINDOW_Y_UM, MODE_WINDOW_Z_UM), ) output_plane = td.Box( center=(centerline.x_end, output_y, LN_THICKNESS_UM / 2), size=(MODE_WINDOW_Y_UM, 0.0, MODE_WINDOW_Z_UM), ) source = td.ModeSource( center=input_plane.center, size=input_plane.size, source_time=td.GaussianPulse( freq0=CENTER_FREQUENCY_HZ, fwidth=0.1 * CENTER_FREQUENCY_HZ, ), mode_spec=mode_spec(), direction="+", mode_index=0, name="source", ) input_x_monitor = input_x + 1.0 mode_in = td.ModeMonitor( center=(input_x_monitor, 0.0, LN_THICKNESS_UM / 2), size=input_plane.size, freqs=FREQUENCIES_HZ, mode_spec=mode_spec(), store_fields_direction="+", name="mode_in", ) mode_out = td.ModeMonitor( center=output_plane.center, size=output_plane.size, freqs=FREQUENCIES_HZ, mode_spec=mode_spec(), store_fields_direction="+", name="mode_out", ) x_min, x_max = -6.0, centerline.x_end + 4.0 y_min, y_max = -4.0, centerline.y_end + 6.0 z_min, z_max = -5.0, 4.0 center = ( (x_min + x_max) / 2, (y_min + y_max) / 2, (z_min + z_max) / 2, ) return td.Simulation( center=center, size=(x_max - x_min, y_max - y_min, z_max - z_min), medium=sio2_medium(), structures=structures, sources=[source], monitors=[mode_in, mode_out], boundary_spec=td.BoundarySpec( x=td.Boundary.absorber(num_layers=80), y=td.Boundary.absorber(num_layers=80), z=td.Boundary.pml(num_layers=16), ), grid_spec=td.GridSpec.auto( wavelength=GRID_WAVELENGTH_UM, min_steps_per_wvl=STEPS_PER_WAVELENGTH, ), run_time=td.RunTimeSpec(quality_factor=1), symmetry=(0, 0, 0), ) def build_eme(centerline: EulerCenterline) -> td.EMESimulation: large = 1e6 ridge_vertices = np.asarray( [ [-large, -RIDGE_TOP_WIDTH_UM / 2], [large, -RIDGE_TOP_WIDTH_UM / 2], [large, RIDGE_TOP_WIDTH_UM / 2], [-large, RIDGE_TOP_WIDTH_UM / 2], ], dtype=float, ) structures = stack_structures(ridge_polyslab(ridge_vertices)) straight_mode = mode_spec(for_eme=True) bend_lengths = np.full(EME_NUM_CELLS, centerline.length / EME_NUM_CELLS) bend_s = (np.arange(EME_NUM_CELLS) + 0.5) * bend_lengths[0] radii = np.interp(bend_s, centerline.s, centerline.curvature) # Negative sign follows the Tidy3D convention for this counter-clockwise # x-to-y turn with the curvature center on the positive transverse side. radii = -1.0 / np.maximum(radii, 1e-12) bend_modes = [ mode_spec( bend_radius=float(radius), for_eme=True, ) for radius in radii ] mode_specs = [straight_mode, *bend_modes, straight_mode] lengths = np.concatenate(([INPUT_LEAD_UM], bend_lengths, [OUTPUT_LEAD_UM])) eme_grid = td.EMEExplicitGrid( boundaries=np.cumsum(lengths)[:-1], mode_specs=mode_specs, ) total_length = float(lengths.sum()) return td.EMESimulation( center=(total_length / 2, 0.0, 0.0), size=(total_length, EME_PLANE_Y_UM, EME_PLANE_Z_UM), medium=sio2_medium(), structures=structures, axis=0, freqs=FREQUENCIES_HZ, eme_grid_spec=eme_grid, grid_spec=td.GridSpec.auto( wavelength=GRID_WAVELENGTH_UM, min_steps_per_wvl=STEPS_PER_WAVELENGTH, ), store_port_modes=False, sweep_spec=None, constraint="passive", ) def build_models() -> tuple[td.Simulation, td.EMESimulation]: centerline = build_euler_centerline() return build_fdtd(centerline), build_eme(centerline) def upload_and_estimate( manifest_path: Path = MANIFEST_PATH, ) -> dict[str, object]: """Upload both models and save exact cost estimates without starting them.""" from tidy3d import web if manifest_path.exists(): manifest = json.loads(manifest_path.read_text()) print(f"Existing task manifest: {manifest_path}") for kind, item in manifest["tasks"].items(): print( f" {kind.upper()}: {item['estimated_cost_flexcredits']:.3f} " f"FlexCredits ยท {item['task_id']}" ) return manifest fdtd_sim, eme_sim = build_models() simulations = { "fdtd": fdtd_sim, "eme": eme_sim, } tasks = {} for kind, simulation in simulations.items(): task_id = web.upload( simulation, task_name=f"{TASK_PREFIX}_{kind}", verbose=True, ) cost = float(web.estimate_cost(task_id, verbose=True)) tasks[kind] = { "task_id": task_id, "task_name": f"{TASK_PREFIX}_{kind}", "estimated_cost_flexcredits": cost, "result_path": str(OUTPUT_DIR / f"{TASK_PREFIX}_{kind}.hdf5"), } aggregate = float( sum(item["estimated_cost_flexcredits"] for item in tasks.values()) ) manifest = { "task_prefix": TASK_PREFIX, "wavelength_um": WAVELENGTHS_UM.tolist(), "crystal_cut": "X-cut", "optic_axis_global": "y", "settings": { "steps_per_wavelength": STEPS_PER_WAVELENGTH, "effective_radius_um": MIN_BEND_RADIUS_UM, "eme_bend_cells": EME_NUM_CELLS, "eme_modes": EME_NUM_MODES, "eme_plane_um": [EME_PLANE_Y_UM, EME_PLANE_Z_UM], "eme_bend_medium_frame": "global", "eme_mode_tracking": "central", }, "tasks": tasks, "aggregate_estimated_cost_flexcredits": aggregate, } OUTPUT_DIR.mkdir(parents=True, exist_ok=True) manifest_path.write_text(json.dumps(manifest, indent=2)) print(f"Aggregate estimate: {aggregate:.3f} FlexCredits") print(f"Task manifest: {manifest_path}") return manifest def start_tasks(manifest_path: Path = MANIFEST_PATH) -> dict[str, str]: """Start the estimated tasks after explicit cost approval.""" from tidy3d import web manifest = json.loads(manifest_path.read_text()) def start_one(item: tuple[str, dict[str, object]]) -> tuple[str, str]: kind, task = item task_id = str(task["task_id"]) status = str(web.get_info(task_id).status).lower() if status == "draft": web.start(task_id) elif status not in {"queued", "pre", "running", "post", "success"}: raise RuntimeError( f"Cannot start {kind.upper()} task from status '{status}'." ) return kind, f"https://tidy3d.simulation.cloud/workbench?taskId={task_id}" with ThreadPoolExecutor(max_workers=2) as executor: return dict(executor.map(start_one, manifest["tasks"].items())) def _fdtd_transmission(path: Path) -> tuple[np.ndarray, np.ndarray]: import h5py with h5py.File(path, "r") as file: data = file["data"] amp_groups = [ data[name]["amps"] for name in sorted(data.keys(), key=int) if "amps" in data[name] ] if len(amp_groups) != 2: raise RuntimeError("Expected two FDTD modal monitor groups.") input_amps = amp_groups[0]["__xarray_dataarray_variable__"][0, :, 0] output_amps = amp_groups[1]["__xarray_dataarray_variable__"][0, :, 0] freqs = np.asarray(amp_groups[0]["f"][:], dtype=float) wavelengths = td.C_0 / freqs order = np.argsort(wavelengths) return wavelengths[order], ( np.abs(output_amps) ** 2 / np.abs(input_amps) ** 2 )[order] def _eme_transmission(path: Path) -> tuple[np.ndarray, np.ndarray]: import h5py with h5py.File(path, "r") as file: s21 = np.asarray(file["smatrix/S21/__xarray_dataarray_variable__"][:]) freqs = np.asarray(file["smatrix/S21/f"][:], dtype=float) if s21.ndim != 4 or s21.shape[1] != 1: raise RuntimeError("Unexpected standalone EME S21 shape.") wavelengths = td.C_0 / freqs order = np.argsort(wavelengths) return wavelengths[order], np.abs(s21[:, 0, 0, 0])[order] ** 2 def _loss_db(transmission: np.ndarray) -> np.ndarray: return -10.0 * np.log10(transmission) def analyze( output_dir: Path = OUTPUT_DIR, manifest_path: Path = MANIFEST_PATH, ) -> dict[str, object]: """Load completed results and plot broadband linear transmission.""" from tidy3d import web output_dir.mkdir(parents=True, exist_ok=True) paths = { kind: output_dir / f"{TASK_PREFIX}_{kind}.hdf5" for kind in ("fdtd", "eme") } manifest = ( json.loads(manifest_path.read_text()) if manifest_path.exists() else None ) for kind, path in paths.items(): if not path.exists(): if manifest is None: raise FileNotFoundError( f"Missing {path} and no task manifest is available." ) task_id = str(manifest["tasks"][kind]["task_id"]) web.download(task_id, path=str(path), verbose=True) fdtd_wavelengths, fdtd_t = _fdtd_transmission(paths["fdtd"]) eme_wavelengths, eme_t = _eme_transmission(paths["eme"]) if not np.allclose(fdtd_wavelengths, eme_wavelengths, rtol=1e-9, atol=1e-12): raise RuntimeError("FDTD and EME frequency grids do not match.") wavelengths = fdtd_wavelengths wavelengths_nm = 1e3 * wavelengths fdtd_loss = _loss_db(fdtd_t) eme_loss = _loss_db(eme_t) error_percent = np.abs(eme_t - fdtd_t) / fdtd_t * 100.0 plt.rcParams.update( { "font.family": "DejaVu Sans", "font.size": 14, "axes.labelsize": 16, "legend.fontsize": 14, "xtick.labelsize": 14, "ytick.labelsize": 14, } ) fig, ax = plt.subplots(figsize=(8.4, 5.4), constrained_layout=True) fig.patch.set_facecolor("white") blue = "#164E8C" red = "#C9364A" ax.set_facecolor("#FBFCFE") ax.grid(True, color="#B8C2CC", alpha=0.35, linewidth=0.8) ax.spines["top"].set_visible(False) ax.spines["right"].set_visible(False) ax.axhline(1.0, color="#64748B", linewidth=1.0, alpha=0.75) eme_line, = ax.plot( wavelengths_nm, eme_t, linewidth=2.7, marker="s", markersize=6.5, color=blue, zorder=2, label="EME", ) fdtd_line, = ax.plot( wavelengths_nm, fdtd_t, linestyle=(0, (5, 3)), linewidth=2.5, marker="o", markersize=7, color=red, zorder=4, label="FDTD", ) ax.set_xlim(float(wavelengths_nm.min() - 3), float(wavelengths_nm.max() + 3)) ax.set_ylim(0.5, 1.0) ax.set_xlabel("Wavelength (nm)") ax.set_ylabel("Fundamental-mode transmission, T") ax.set_xticks(wavelengths_nm) ax.legend( handles=[fdtd_line, eme_line], loc="lower left", frameon=True, framealpha=0.95, ) plot_path = output_dir / f"{TASK_PREFIX}_transmission.png" fig.savefig(plot_path, dpi=220, bbox_inches="tight") plt.close(fig) results = { "wavelength_um": wavelengths.tolist(), "wavelength_nm": wavelengths_nm.tolist(), "effective_radius_um": MIN_BEND_RADIUS_UM, "crystal_cut": "X-cut", "optic_axis_global": "y", "fdtd": { "transmission": fdtd_t.tolist(), "bend_loss_db": fdtd_loss.tolist(), }, "eme": { "num_cells": EME_NUM_CELLS, "num_modes": EME_NUM_MODES, "plane_um": [EME_PLANE_Y_UM, EME_PLANE_Z_UM], "bend_medium_frame": "global", "mode_tracking": "central", "transmission": eme_t.tolist(), "bend_loss_db": eme_loss.tolist(), }, "relative_transmission_error_percent": error_percent.tolist(), "max_relative_transmission_error_percent": float(error_percent.max()), "mean_relative_transmission_error_percent": float(error_percent.mean()), } json_path = output_dir / f"{TASK_PREFIX}_transmission.json" with json_path.open("w") as file: json.dump(results, file, indent=2) print( f"FDTD/EME comparison: max error={error_percent.max():.3f}%, " f"mean error={error_percent.mean():.3f}%" ) print(f"Plot: {plot_path}") print(f"Data: {json_path}") return results def print_build_summary( fdtd_model: td.Simulation, eme_model: td.EMESimulation, ) -> None: """Print the key physical and numerical settings.""" print("Models built locally.") print(" crystal: X-cut LiNbO3, optic axis along global y") print(f" wavelengths: {(1e3 * WAVELENGTHS_UM).tolist()} nm") print(f" effective radius: {MIN_BEND_RADIUS_UM:g} um") print(f" mesh: wavelength/{STEPS_PER_WAVELENGTH}") print(f" FDTD domain: {tuple(round(value, 3) for value in fdtd_model.size)} um") print(f" EME: N={EME_NUM_CELLS}, M={EME_NUM_MODES}") print(f" EME plane: {EME_PLANE_Y_UM:g} x {EME_PLANE_Z_UM:g} um") print(" EME broadband mode tracking: central frequency") print(f" EME cells: {len(eme_model.eme_grid.mode_specs)} total") def main() -> None: parser = argparse.ArgumentParser(description=__doc__) parser.add_argument( "command", nargs="?", choices=("build", "estimate", "run", "analyze"), default="build", ) args = parser.parse_args() if args.command == "build": fdtd_model, eme_model = build_models() print_build_summary(fdtd_model, eme_model) elif args.command == "estimate": upload_and_estimate() elif args.command == "run": urls = start_tasks() for kind, url in urls.items(): print(f"{kind.upper()}: {url}") else: analyze() if __name__ == "__main__": main()