"""`ged.sim` — verify a recipe's assumed group index against a REAL mode solve. GED's photonic recipes take the waveguide group index ``n_g`` as a GIVEN literature constant (the ring recipe is ``R = c / (2*pi*n_g*FSR)`` with ``n_g = 4.2``). This module closes that loop honestly: it runs a real electromagnetic **mode solve** of the waveguide cross-section and checks whether the simulated ``n_g`` reproduces the recipe's assumed value within a physical tolerance. A faithful sim confirms the anchor; a recipe whose assumed ``n_g`` no longer matches the physics is a DETECTED violation (REJECTED), never a crash. Pure-core / thin-edge, per the repo convention: * ``verify_group_index(simulated_n_g, recipe, tolerance)`` — PURE: reuses ``solve._correctness_verdict``. No network, no clock, no disk. The mutmut target. * ``ModeRunner`` — the injectable EDGE (a ``Callable[[WaveguideSpec], ModeResult]``), exactly the ``serve.py`` ``pack_materializer`` pattern: inject a fake in tests (CI stays offline, free, deterministic), bind a real backend in production. * ``femwell_runner`` — a REAL, LOCAL FEM mode solve (the free path; no API key, no cloud cost). ``tidy3d_runner`` — the Tidy3D cloud adapter (needs a Tidy3D API key + cloud credits). Both lazy-import their heavy backend (the ``[sim]`` extra), so ``ged.sim``'s pure core imports with the stdlib alone. A real solve costs real compute and (for tidy3d) money, so it is OPT-IN: nothing here runs a backend unless a runner is explicitly invoked. The default everywhere is no runner -> the ``n_g`` claim stays UNVERIFIED (planned, not available) — the same honesty boundary as ``validate --strict`` and the ``code-patching-arena`` gate. """ from __future__ import annotations from dataclasses import dataclass from typing import Callable, Mapping from .solve import _correctness_verdict # A simulated n_g may legitimately differ from a rounded literature constant by a # couple percent (mesh, material model, geometry idealization). This is a PHYSICAL # tolerance — deliberately looser than the 1e-3 used for pure re-derivation identity # and the 1e-2 anchor check. A constant-index model (no material dispersion) lands # ~3% off and correctly REJECTS; a dispersion-correct solve lands <1% and ACCEPTS. GROUP_INDEX_TOLERANCE = 0.02 # Verdict when a recipe carries no n_g to check (bragg / coupler / phase-shifter). NO_GROUP_INDEX = "UNVERIFIED" @dataclass(frozen=True) class WaveguideSpec: """The cross-section a mode solve needs: core size (um), the wavelength (um), the core/clad refractive indices at that wavelength, and their linear material dispersion ``dn/dlambda`` (per um). Pure data — a runner turns it into a real simulation. The dispersion slopes matter: the group index ``n_g = n_eff - lambda*dn_eff/dlambda`` has BOTH a waveguide-dispersion term (geometry) and a material-dispersion term (``dn/dlambda`` of the core). A solve that holds the index constant captures only the first and lands a few percent low — for silicon at 1550 nm that is the difference between ~4.05 (REJECTED) and ~4.20 (ACCEPTED vs the recipe's 4.2).""" core_width_um: float core_thickness_um: float wavelength_um: float core_index: float clad_index: float core_index_slope_per_um: float = 0.0 # dn_core/dlambda (material dispersion) clad_index_slope_per_um: float = 0.0 # dn_clad/dlambda def core_index_at(self, wavelength_um: float) -> float: return self.core_index + self.core_index_slope_per_um * (wavelength_um - self.wavelength_um) def clad_index_at(self, wavelength_um: float) -> float: return self.clad_index + self.clad_index_slope_per_um * (wavelength_um - self.wavelength_um) @dataclass(frozen=True) class ModeResult: """What a mode solve returns for the fundamental mode: the effective and group indices. A runner (real or fake) produces this; the pure verifier consumes it.""" n_eff: float n_group: float backend: str # "femwell" | "tidy3d" | "fake" — provenance of the number # The injectable edge: spec -> result. Production binds femwell_runner / tidy3d_runner; # tests bind a fake returning a captured float. Never call a real backend in pure logic. ModeRunner = Callable[[WaveguideSpec], ModeResult] def verify_group_index( simulated_n_g: float, recipe: Mapping[str, object], tolerance: float = GROUP_INDEX_TOLERANCE, ) -> str: """PURE: does a simulated group index reproduce the recipe's assumed ``n_g``? ACCEPTED / REJECTED within ``tolerance`` (relative), reusing the same numeric oracle every other GED layer uses. UNVERIFIED when the recipe carries no ``n_g`` to check (bragg / coupler / phase-shifter store ``n_eff`` / ``kappa`` / ``delta_n`` instead) — a recipe without the quantity is honestly unverified, never a KeyError and never a vacuous pass.""" params = recipe.get("default_params") if not isinstance(params, Mapping): return NO_GROUP_INDEX assumed = params.get("n_g") if not isinstance(assumed, (int, float)): return NO_GROUP_INDEX return _correctness_verdict(simulated_n_g, float(assumed), tolerance, "relative") def recipe_has_group_index(recipe: Mapping[str, object]) -> bool: """PURE: whether ``recipe`` declares an ``n_g`` a mode solve could verify.""" params = recipe.get("default_params") return isinstance(params, Mapping) and isinstance(params.get("n_g"), (int, float)) def verify_recipe_group_index( recipe: Mapping[str, object], runner: ModeRunner, spec: WaveguideSpec, tolerance: float = GROUP_INDEX_TOLERANCE, ) -> tuple[str, ModeResult]: """EDGE: run ``runner`` on ``spec`` (a real or fake mode solve), then verify the simulated group index against the recipe — returns ``(verdict, ModeResult)``. The runner is the only impure thing; the verdict is the pure oracle above. The caller supplies the runner, so CI injects a fake and production injects a real backend.""" result = runner(spec) verdict = verify_group_index(result.n_group, recipe, tolerance) return verdict, result # --- real backends (lazy-imported; the [sim] extra) ------------------------- # These are thin adapters at the very edge: each imports its heavy library INSIDE # the function so `import ged.sim` never pulls femwell/tidy3d, and each returns the # shared ModeResult so the pure verifier is backend-agnostic. They are NOT exercised # in CI (a real solve blows the 2 GiB / 30 s test limits) — an opt-in script or the # live goal probe drives them. def femwell_runner(spec: WaveguideSpec, *, mesh_resolution_um: float = 0.015) -> ModeResult: """REAL local FEM mode solve via femwell (free; no API key, no cloud cost). Solves the waveguide eigenproblem on a meshed cross-section, then computes the group index by central difference of n_eff over wavelength (``n_g = n_eff - lambda * dn_eff/dlambda``). Material dispersion enters through ``spec.core_index_at`` / ``spec.clad_index_at``, which shift the indices across the +/- wavelength step — so the result includes BOTH waveguide and material dispersion. For the dispersion-correct Si strip this lands within ~0.5% of the cloud-FDTD reference (and ~0.02% of the recipe's 4.2); a zero-slope spec omits material dispersion and lands ~3% low. Heavy import is local so the module's pure core stays stdlib-only.""" from collections import OrderedDict import numpy as np from femwell.maxwell.waveguide import compute_modes from femwell.mesh import mesh_from_OrderedDict from shapely.geometry import box from skfem import Basis, ElementTriP0 from skfem.io.meshio import from_meshio def n_eff_at(wavelength: float) -> float: w, t = spec.core_width_um, spec.core_thickness_um polygons = OrderedDict( core=box(-w / 2, 0, w / 2, t), clad=box(-2.0, -2.0, 2.0, t + 2.0), ) resolutions = {"core": {"resolution": mesh_resolution_um, "distance": 0.6}} mesh = from_meshio( mesh_from_OrderedDict(polygons, resolutions, default_resolution_max=0.25) ) basis = Basis(mesh, ElementTriP0()) eps = basis.zeros() eps[basis.get_dofs(elements="core")] = spec.core_index_at(wavelength) ** 2 eps[basis.get_dofs(elements="clad")] = spec.clad_index_at(wavelength) ** 2 modes = compute_modes(basis, eps, wavelength=wavelength, num_modes=1, order=2) return float(np.real(modes[0].n_eff)) wl = spec.wavelength_um step = 0.01 n0 = n_eff_at(wl) dn_dwl = (n_eff_at(wl + step) - n_eff_at(wl - step)) / (2 * step) n_group = n0 - wl * dn_dwl return ModeResult(n_eff=n0, n_group=n_group, backend="femwell") def tidy3d_runner( spec: WaveguideSpec, *, core_material: str = "cSi", core_variant: str = "Li1993_293K", clad_material: str = "SiO2", clad_variant: str = "Horiba", ) -> ModeResult: """REAL cloud mode solve via the Tidy3D cloud FDTD/mode solver. Needs a Tidy3D API key and costs a few cloud credits (cents) per run — the group index is only returned by the SERVER run (the local solver omits it), so this submits to the cloud. Uses Tidy3D's DISPERSIVE material library (cSi / SiO2) rather than the spec's scalar indices, so material dispersion — the term a constant-index model omits — is included by construction (the documented minimal path). Returns the same ModeResult as femwell, so the pure verifier and the goal criterion are backend-agnostic. Heavy/networked import lives here at the edge.""" import tidy3d as td from tidy3d.plugins import waveguide strip = waveguide.RectangularDielectric( wavelength=spec.wavelength_um, core_width=spec.core_width_um, core_thickness=spec.core_thickness_um, core_medium=td.material_library[core_material][core_variant], clad_medium=td.material_library[clad_material][clad_variant], mode_spec=td.ModeSpec(num_modes=1, group_index_step=True), ) # .n_eff / .n_group are (num_modes,)-shaped; take the fundamental mode. n_eff = float(strip.n_eff.values.flatten()[0].real) n_group = float(strip.n_group.values.flatten()[0].real) return ModeResult(n_eff=n_eff, n_group=n_group, backend="tidy3d") # The verified reference cross-section: a 500x220 nm SOI strip at 1550 nm with # literature silicon dispersion, for which femwell returns n_g ~= 4.20 (recipe 4.2) # and tidy3d ~= 4.18. Bundled so the live probe / CLI has a concrete spec to run. SOI_STRIP_1550 = WaveguideSpec( core_width_um=0.5, core_thickness_um=0.22, wavelength_um=1.55, core_index=3.476, # crystalline Si @ 1550 nm (Li 1993) clad_index=1.444, # SiO2 @ 1550 nm core_index_slope_per_um=-0.086, # dn_Si/dlambda @ 1550 nm (Li 1993 Sellmeier) clad_index_slope_per_um=-0.010, # dn_SiO2/dlambda @ 1550 nm ) __all__ = [ "GROUP_INDEX_TOLERANCE", "NO_GROUP_INDEX", "WaveguideSpec", "ModeResult", "ModeRunner", "SOI_STRIP_1550", "verify_group_index", "recipe_has_group_index", "verify_recipe_group_index", "femwell_runner", "tidy3d_runner", ]