diff --git a/nbs/meep_grating_coupler.ipynb b/nbs/meep_grating_coupler.ipynb new file mode 100644 index 00000000..703a0a76 --- /dev/null +++ b/nbs/meep_grating_coupler.ipynb @@ -0,0 +1,266 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "0", + "metadata": {}, + "source": [ + "# Grating Coupler Simulation with FiberSource\n", + "\n", + "[MEEP](https://meep.readthedocs.io/) is an open-source FDTD electromagnetic simulator. This notebook demonstrates using the `gsim.meep` **FiberSource** API to simulate fiber-to-chip coupling through a grating coupler.\n", + "\n", + "Unlike S-parameter simulations that use eigenmode sources at waveguide ports, grating coupler simulations launch a **Gaussian beam from above** (simulating a fiber) and measure the power coupled into the waveguide port via eigenmode decomposition.\n", + "\n", + "**Requirements:**\n", + "\n", + "- UBC PDK: `uv pip install ubcpdk`\n", + "- [GDSFactory+](https://gdsfactory.com) account for cloud simulation" + ] + }, + { + "cell_type": "markdown", + "id": "1", + "metadata": {}, + "source": [ + "### Load a grating coupler from UBC PDK" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2", + "metadata": {}, + "outputs": [], + "source": [ + "from ubcpdk import PDK, cells\n", + "\n", + "PDK.activate()\n", + "\n", + "c = cells.ebeam_gc_te1550()\n", + "c" + ] + }, + { + "cell_type": "markdown", + "id": "3", + "metadata": {}, + "source": [ + "### Configure simulation with FiberSource\n", + "\n", + "Instead of `ModeSource` (which excites a waveguide eigenmode at a port), we use `FiberSource` to launch a Gaussian beam from above the grating. Key parameters:\n", + "\n", + "- **beam_waist**: Gaussian beam waist radius (SMF-28 mode field diameter / 2 = 5.2 um)\n", + "- **angle_theta**: Fiber tilt angle from vertical (typically 8-15 degrees for grating couplers)\n", + "- **z_offset**: Distance above the chip surface to place the source plane\n", + "- **polarization**: TE or TM\n", + "\n", + "The result is a `CouplingResult` instead of `SParameterResult`, giving coupling efficiency (CE) per port." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4", + "metadata": {}, + "outputs": [], + "source": [ + "from gsim import meep\n", + "\n", + "sim = meep.Simulation()\n", + "\n", + "sim.geometry(component=c, z_crop=\"auto\")\n", + "sim.materials = {\"si\": 3.47, \"SiO2\": 1.44}\n", + "\n", + "sim.source = meep.FiberSource(\n", + " wavelength=1.55,\n", + " wavelength_span=0.07,\n", + " num_freqs=21,\n", + " beam_waist=5.2,\n", + " angle_theta=0.0,\n", + " z_offset=2.0,\n", + " position=[-25, 0],\n", + ")\n", + "\n", + "sim.monitors = [\"o1\"]\n", + "sim.domain(pml=1.0, margin=1.0, margin_z_above=3.0)\n", + "sim.solver(resolution=20, save_animation=False, verbose_interval=5.0)\n", + "sim.solver.stop_after_sources(time=50)\n", + "sim.solver.animation_plane = \"xz\"\n", + "\n", + "print(sim.validate_config())" + ] + }, + { + "cell_type": "markdown", + "id": "5", + "metadata": {}, + "source": [ + "### Preview geometry" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6", + "metadata": {}, + "outputs": [], + "source": [ + "sim.plot_2d(slices=\"xyz\")" + ] + }, + { + "cell_type": "markdown", + "id": "7", + "metadata": {}, + "source": [ + "### Run simulation on cloud" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8", + "metadata": {}, + "outputs": [], + "source": [ + "result = sim.run()" + ] + }, + { + "cell_type": "markdown", + "id": "9", + "metadata": {}, + "source": [ + "### Plot coupling efficiency" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "10", + "metadata": {}, + "outputs": [], + "source": [ + "result.plot(db=True)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "11", + "metadata": {}, + "outputs": [], + "source": [ + "print(\"Peak coupling efficiency (dB):\")\n", + "for port, ce_db in result.peak_ce.items():\n", + " print(f\" {port}: {ce_db:.2f} dB\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "12", + "metadata": {}, + "outputs": [], + "source": [ + "result.show_animation()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "fd1b0fc6", + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "from matplotlib.animation import FuncAnimation\n", + "from matplotlib.colors import LogNorm\n", + "from IPython.display import HTML\n", + "\n", + "data = result.load_animation_data()\n", + "fields_abs = np.abs(data[\"fields\"])\n", + "vmax = fields_abs.max()\n", + "floor = vmax * 1e-4 # 4 orders of magnitude below peak\n", + "\n", + "fig, ax = plt.subplots(figsize=(6, 4))\n", + "ax.imshow(data[\"eps_data\"].T, origin=\"lower\", extent=data[\"extent\"], cmap=\"binary\")\n", + "im = ax.imshow(\n", + " np.clip(fields_abs[0].T, floor, None),\n", + " origin=\"lower\",\n", + " extent=data[\"extent\"],\n", + " cmap=\"hot\",\n", + " alpha=0.8,\n", + " norm=LogNorm(vmin=floor, vmax=vmax),\n", + ")\n", + "fig.colorbar(im, ax=ax, label=\"|Ey|\")\n", + "ax.set_xlabel(\"x (µm)\")\n", + "ax.set_ylabel(\"y (µm)\")\n", + "title = ax.set_title(f\"|Ey| (log) t={data['times'][0]:.2f}\")\n", + "\n", + "\n", + "def update(i):\n", + " im.set_data(np.clip(fields_abs[i].T, floor, None))\n", + " title.set_text(f\"|Ey| (log) t={data['times'][i]:.2f}\")\n", + " return [im, title]\n", + "\n", + "\n", + "ani = FuncAnimation(fig, update, frames=len(data[\"times\"]), interval=67)\n", + "plt.close(fig)\n", + "HTML(ani.to_jshtml())" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e6037483", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c305697c", + "metadata": {}, + "outputs": [], + "source": [ + "from gsim import gcloud\n", + "\n", + "job_id = gcloud.upload(\"gc_diagnostic_v2a\", \"meep\")\n", + "gcloud.start(job_id)\n", + "result = gcloud.wait_for_results(job_id, verbose=\"full\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c6aa0f29", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "gsim", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.12.10" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/pyproject.toml b/pyproject.toml index 8a2478e1..7303e7bb 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -83,7 +83,7 @@ src = "src/gsim/__init__.py" commit-template = "Release {new-version}" [tool.codespell] -ignore-words-list = "doubleclick" +ignore-words-list = "doubleclick,TE" [tool.interrogate] docstring-style = "google" diff --git a/src/gsim/common/viz/render2d.py b/src/gsim/common/viz/render2d.py index ecf57ba5..9ceac400 100644 --- a/src/gsim/common/viz/render2d.py +++ b/src/gsim/common/viz/render2d.py @@ -15,7 +15,7 @@ import matplotlib.pyplot as plt import numpy as np -from matplotlib.patches import Polygon, Rectangle +from matplotlib.patches import Circle, FancyArrowPatch, Polygon, Rectangle from gsim.common.geometry_model import GeometryModel @@ -618,6 +618,11 @@ def _draw_overlay( elif y is not None: _draw_overlay_xz(ax, cmin, cmax, pml, overlay.ports, y) + # Draw fiber source overlay + fiber = getattr(overlay, "fiber_source", None) + if fiber is not None: + _draw_fiber_source(ax, fiber, x=x, y=y, z=z) + def _draw_overlay_xy( ax: plt.Axes, @@ -846,6 +851,154 @@ def _draw_overlay_xz( ) +_FIBER_COLOR = (0.2, 0.7, 0.2) # green +_FIBER_BEAM_ALPHA = 0.15 + + +def _draw_fiber_source( + ax: plt.Axes, + fiber: Any, + *, + x: float | None, + y: float | None, + z: float | None, +) -> None: + """Draw fiber source Gaussian beam indicator on 2D slices. + + - **XZ** (y-slice): beam cone from source plane toward chip + - **YZ** (x-slice): beam cone if beam center is near x-slice + - **XY** (z-slice): beam footprint circle at the slice z + """ + fx, fy = fiber.position + fz = fiber.z_position + w0 = fiber.beam_waist + theta = np.radians(fiber.angle_theta) + phi = np.radians(fiber.angle_phi) + sign = -1.0 if fiber.direction == "down" else 1.0 + + if y is not None: + # XZ view — draw beam if center is near the y-slice + if abs(fy - y) > w0: + return + _draw_fiber_beam_side(ax, fx, fz, w0, theta, phi, sign, axis="xz") + elif x is not None: + # YZ view — draw beam if center is near the x-slice + if abs(fx - x) > w0: + return + _draw_fiber_beam_side(ax, fy, fz, w0, theta, phi, sign, axis="yz") + elif z is not None: + # XY view — draw beam footprint circle + _draw_fiber_beam_xy(ax, fx, fy, fz, z, w0, theta, phi, sign) + + +def _draw_fiber_beam_side( + ax: plt.Axes, + h_center: float, + z_src: float, + w0: float, + theta: float, + phi: float, + sign: float, + axis: str = "xz", +) -> None: + """Draw the Gaussian beam cone on XZ or YZ side views.""" + # Beam propagation length (visual, extends to source z ± some distance) + beam_len = 4.0 * w0 + + # Direction components: theta is tilt from vertical, phi selects plane + phi_comp = np.cos(phi) if axis == "xz" else np.sin(phi) + dh = np.sin(theta) * phi_comp + dz = sign * np.cos(theta) + + # Beam center line endpoints + h0, z0 = h_center, z_src + h1 = h0 + dh * beam_len + z1 = z0 + dz * beam_len + + # Draw beam envelope (expanding cone) + # At source: width = w0, at end: width = w0 * 1.5 (divergence visual) + perp_h = -dz # perpendicular in 2D + perp_z = dh + + norm = np.sqrt(perp_h**2 + perp_z**2) + if norm > 0: + perp_h /= norm + perp_z /= norm + + cone = Polygon( + [ + (h0 + perp_h * w0, z0 + perp_z * w0), + (h0 - perp_h * w0, z0 - perp_z * w0), + (h1 - perp_h * w0 * 1.5, z1 - perp_z * w0 * 1.5), + (h1 + perp_h * w0 * 1.5, z1 + perp_z * w0 * 1.5), + ], + facecolor=(*_FIBER_COLOR, _FIBER_BEAM_ALPHA), + edgecolor=(*_FIBER_COLOR, 0.5), + linewidth=0.8, + zorder=92, + label="Fiber source", + ) + ax.add_patch(cone) + + # Central beam axis arrow + ax.add_patch( + FancyArrowPatch( + (h0, z0), + (h1, z1), + arrowstyle="->", + color=_FIBER_COLOR, + linewidth=1.5, + zorder=93, + ) + ) + + # Source plane marker (horizontal line at z_src) + ax.plot( + [h0 - w0, h0 + w0], + [z0, z0], + color=_FIBER_COLOR, + linewidth=2, + linestyle="-", + zorder=93, + ) + + +def _draw_fiber_beam_xy( + ax: plt.Axes, + fx: float, + fy: float, + fz: float, + z_slice: float, + w0: float, + theta: float, + phi: float, + sign: float, +) -> None: + """Draw the Gaussian beam footprint on an XY (z-slice) view.""" + # Offset the beam center based on propagation from source to slice + dz = z_slice - fz + if abs(dz) < 1e-6: + cx, cy = fx, fy + else: + # Lateral shift due to tilt + cx = fx + dz * np.tan(theta) * np.cos(phi) / sign + cy = fy + dz * np.tan(theta) * np.sin(phi) / sign + + ax.add_patch( + Circle( + (cx, cy), + w0, + facecolor=(*_FIBER_COLOR, _FIBER_BEAM_ALPHA), + edgecolor=_FIBER_COLOR, + linewidth=1.5, + linestyle="--", + zorder=92, + label="Fiber source", + ) + ) + ax.plot(cx, cy, "+", color=_FIBER_COLOR, markersize=8, zorder=93) + + def _add_pml_rect( ax: plt.Axes, x: float, diff --git a/src/gsim/meep/ROADMAP.md b/src/gsim/meep/ROADMAP.md new file mode 100644 index 00000000..f1318030 --- /dev/null +++ b/src/gsim/meep/ROADMAP.md @@ -0,0 +1,16 @@ +# gsim.meep Roadmap + +Features planned beyond the current S-parameter extraction workflow. + +## Planned + +- [ ] **Grating coupler simulation** — near-to-far-field transformation, fiber mode overlap, coupling efficiency vs + wavelength +- [ ] **Resonator & Q-factor extraction** — Harminv-based resonance finding, Q/FSR extraction, SAX circuit model + integration +- [ ] **Parameter sweep framework** — sweep any sim parameter, parallel cloud runs, result aggregation +- [ ] **Convergence testing** — automated resolution/domain sweeps, convergence metric, pass/fail report +- [ ] **Broadband flux monitors** — total transmission/reflection via `add_flux`, for filters/gratings/stacks +- [ ] **Dispersive materials** — Lorentz/Drude models, MEEP built-in material library (~20 fitted materials) +- [ ] **Band structure / photonic crystal** — Bloch-periodic boundaries, k-point sweeps, band diagram plotting +- [ ] **Adjoint optimization / inverse design** — topology optimization via `meep.adjoint`, design region, GDS export diff --git a/src/gsim/meep/__init__.py b/src/gsim/meep/__init__.py index 4bcd6055..13e6b0dc 100644 --- a/src/gsim/meep/__init__.py +++ b/src/gsim/meep/__init__.py @@ -21,8 +21,11 @@ from gsim.gcloud import RunResult, register_result_parser from gsim.meep.models import ( FDTD, + CouplingResult, Domain, DomainConfig, + FiberSource, + FiberSourceConfig, Geometry, Material, ModeSource, @@ -37,8 +40,13 @@ from gsim.meep.simulation import BuildResult, Simulation -def _parse_meep_result(run_result: RunResult) -> SParameterResult: - """Parse MEEP cloud results into an SParameterResult.""" +def _parse_meep_result(run_result: RunResult) -> SParameterResult | CouplingResult: + """Parse MEEP cloud results into SParameterResult or CouplingResult.""" + # Check for fiber source results first + ce_csv = run_result.files.get("coupling_efficiency.csv") + if ce_csv is not None: + return CouplingResult.from_csv(ce_csv) + # Standard S-parameter results csv_path = run_result.files.get("s_parameters.csv") if csv_path is not None: return SParameterResult.from_csv(csv_path) @@ -50,8 +58,11 @@ def _parse_meep_result(run_result: RunResult) -> SParameterResult: __all__ = [ "FDTD", "BuildResult", + "CouplingResult", "Domain", "DomainConfig", + "FiberSource", + "FiberSourceConfig", "Geometry", "Material", "ModeSource", diff --git a/src/gsim/meep/models/__init__.py b/src/gsim/meep/models/__init__.py index d7cd302d..f2735612 100644 --- a/src/gsim/meep/models/__init__.py +++ b/src/gsim/meep/models/__init__.py @@ -3,6 +3,7 @@ from gsim.meep.models.api import ( FDTD, Domain, + FiberSource, Geometry, Material, ModeSource, @@ -12,6 +13,7 @@ AccuracyConfig, DiagnosticsConfig, DomainConfig, + FiberSourceConfig, LayerStackEntry, MaterialData, PortData, @@ -22,7 +24,7 @@ SymmetryEntry, WavelengthConfig, ) -from gsim.meep.models.results import SParameterResult +from gsim.meep.models.results import CouplingResult, SParameterResult # Backward compatibility alias FDTDConfig = WavelengthConfig @@ -30,10 +32,13 @@ __all__ = [ "FDTD", "AccuracyConfig", + "CouplingResult", "DiagnosticsConfig", "Domain", "DomainConfig", "FDTDConfig", + "FiberSource", + "FiberSourceConfig", "Geometry", "LayerStackEntry", "Material", diff --git a/src/gsim/meep/models/api.py b/src/gsim/meep/models/api.py index 3b7ea65c..c1ca962d 100644 --- a/src/gsim/meep/models/api.py +++ b/src/gsim/meep/models/api.py @@ -90,6 +90,80 @@ def __call__(self, **kwargs: Any) -> ModeSource: return self +class FiberSource(BaseModel): + """Fiber-to-chip coupling via the reciprocal method. + + Simulates grating coupler coupling efficiency by launching an + eigenmode source from the waveguide port and measuring coupling + into a fiber-like mode above the grating. By reciprocity, this + gives the same S-parameters as fiber→waveguide. + + The fiber monitor parameters (angle_theta, beam_waist, position, + z_offset) define where and how the fiber mode is measured above + the chip, not a physical source. + """ + + model_config = ConfigDict(validate_assignment=True) + + port: str | None = Field( + default=None, + description="Waveguide port to excite. None = auto-select first port.", + ) + wavelength: float = Field( + default=1.55, + gt=0, + description="Center wavelength in um", + ) + wavelength_span: float = Field( + default=0.1, + ge=0, + description="Wavelength span of the measurement frequency grid in um.", + ) + num_freqs: int = Field( + default=11, + ge=1, + description="Number of frequency points", + ) + beam_waist: float = Field( + default=5.2, + gt=0, + description="Gaussian beam waist radius in um (SMF-28 MFD/2 ≈ 5.2)", + ) + angle_theta: float = Field( + default=10.0, + ge=0, + lt=90, + description="Fiber angle from vertical in degrees", + ) + angle_phi: float = Field( + default=0.0, + description="Azimuthal angle in degrees (0 = tilted in XZ plane)", + ) + polarization: Literal["TE", "TM"] = Field( + default="TE", + description="Fiber polarization", + ) + position: list[float] | None = Field( + default=None, + description="Fiber center [x, y] above grating. None = auto (bbox center).", + ) + z_offset: float = Field( + default=1.0, + gt=0, + description="Distance above core top to place fiber monitor plane (um)", + ) + direction: Literal["down", "up"] = Field( + default="down", + description="Fiber direction ('down' = fiber above chip looking down)", + ) + + def __call__(self, **kwargs: Any) -> FiberSource: + """Update fields in place. Returns self for chaining.""" + for k, v in kwargs.items(): + setattr(self, k, v) + return self + + # --------------------------------------------------------------------------- # Domain # --------------------------------------------------------------------------- @@ -331,5 +405,6 @@ def __call__(self, **kwargs: Any) -> FDTD: save_epsilon_raw: bool = Field(default=False) save_animation: bool = Field(default=False) animation_interval: float = Field(default=0.5, gt=0) + animation_plane: Literal["xy", "xz"] = Field(default="xy") preview_only: bool = Field(default=False) verbose_interval: float = Field(default=0, ge=0) diff --git a/src/gsim/meep/models/config.py b/src/gsim/meep/models/config.py index f7cc7073..6d00e127 100644 --- a/src/gsim/meep/models/config.py +++ b/src/gsim/meep/models/config.py @@ -8,9 +8,17 @@ import json from pathlib import Path -from typing import Any, Literal +from typing import Annotated, Any, Literal -from pydantic import BaseModel, ConfigDict, Field, PrivateAttr, computed_field +from pydantic import ( + BaseModel, + ConfigDict, + Discriminator, + Field, + PrivateAttr, + Tag, + computed_field, +) class SymmetryEntry(BaseModel): @@ -111,7 +119,7 @@ class StoppingConfig(BaseModel): class SourceConfig(BaseModel): - """Source excitation configuration. + """Eigenmode source excitation configuration. Controls the Gaussian source bandwidth and which port is excited. When ``bandwidth`` is ``None`` (auto), ``compute_fwidth`` returns a @@ -122,6 +130,10 @@ class SourceConfig(BaseModel): model_config = ConfigDict(validate_assignment=True) + source_type: Literal["mode"] = Field( + default="mode", + description="Source type discriminator.", + ) bandwidth: float | None = Field( default=None, description=( @@ -161,6 +173,39 @@ def compute_fwidth(self, fcen: float, monitor_df: float) -> float: return max(3 * monitor_df, 0.2 * fcen) +class FiberSourceConfig(BaseModel): + """Reciprocal fiber coupling config for grating coupler simulation. + + Sent as JSON to the cloud runner. The runner builds an + ``mp.EigenModeSource`` at the waveguide ``port`` and a fiber + mode monitor at ``z_position``. By reciprocity, the coupling + into the fiber mode equals the fiber→waveguide coupling. + """ + + model_config = ConfigDict(validate_assignment=True) + + source_type: Literal["fiber"] = Field( + default="fiber", + description="Source type discriminator.", + ) + port: str | None = Field( + default=None, + description="Waveguide port for EigenModeSource. None = auto (first port).", + ) + beam_waist: float = Field(gt=0, description="Fiber mode waist radius in um") + angle_theta: float = Field( + ge=0, lt=90, description="Fiber angle from vertical (deg)" + ) + angle_phi: float = Field(description="Azimuthal angle (deg)") + polarization: Literal["TE", "TM"] + direction: Literal["down", "up"] + position: list[float] = Field(description="Fiber center [x, y] above grating") + z_position: float = Field( + description="Absolute z-coordinate of fiber monitor plane" + ) + fwidth: float = Field(gt=0, description="Source fwidth in frequency units") + + class WavelengthConfig(BaseModel): """Wavelength and frequency settings for MEEP FDTD simulation. @@ -251,6 +296,10 @@ class DiagnosticsConfig(BaseModel): gt=0, description="MEEP time units between animation frames", ) + animation_plane: Literal["xy", "xz"] = Field( + default="xy", + description="Cross-section plane for field animation: 'xy' or 'xz'", + ) preview_only: bool = Field( description="Init sim and save geometry diagnostics, skip FDTD run", ) @@ -322,7 +371,11 @@ class SimConfig(BaseModel): ports: list[PortData] materials: dict[str, MaterialData] wavelength: WavelengthConfig = Field(serialization_alias="fdtd") - source: SourceConfig + source: Annotated[ + Annotated[SourceConfig, Tag("mode")] + | Annotated[FiberSourceConfig, Tag("fiber")], + Discriminator("source_type"), + ] stopping: StoppingConfig resolution: ResolutionConfig domain: DomainConfig diff --git a/src/gsim/meep/models/results.py b/src/gsim/meep/models/results.py index 9a8bd6b1..037ff48d 100644 --- a/src/gsim/meep/models/results.py +++ b/src/gsim/meep/models/results.py @@ -106,16 +106,12 @@ def from_csv(cls, path: str | Path) -> SParameterResult: ("geometry_yz", "meep_geometry_yz.png"), ("fields_xy", "meep_fields_xy.png"), ("animation", "meep_animation.mp4"), + ("animation_data", "meep_animation_data.npz"), ]: img_path = path.parent / filename if img_path.exists(): diagnostic_images[key] = str(img_path) - # Detect animation frame PNGs - frame_pngs = sorted(path.parent.glob("meep_frame_*.png")) - if frame_pngs: - diagnostic_images["animation_frames"] = str(path.parent) - return cls( wavelengths=wavelengths, s_params=s_params, @@ -156,16 +152,12 @@ def from_directory(cls, directory: str | Path) -> SParameterResult: ("geometry_yz", "meep_geometry_yz.png"), ("fields_xy", "meep_fields_xy.png"), ("animation", "meep_animation.mp4"), + ("animation_data", "meep_animation_data.npz"), ]: img_path = directory / filename if img_path.exists(): diagnostic_images[key] = str(img_path) - # Detect animation frame PNGs - frame_pngs = sorted(directory.glob("meep_frame_*.png")) - if frame_pngs: - diagnostic_images["animation_frames"] = str(directory) - return cls( debug_info=debug_info, diagnostic_images=diagnostic_images, @@ -196,6 +188,30 @@ def show_animation(self) -> None: display(Video(mp4_path, embed=True, mimetype="video/mp4")) + def load_animation_data(self) -> dict[str, Any] | None: + """Load raw animation frame data from the consolidated .npz. + + Returns: + Dict with keys ``fields`` (N, H, W), ``times`` (N,), + ``eps_data`` (H, W), ``extent`` (4,), ``plane`` (str), + or None if not available. + """ + npz_path = self.diagnostic_images.get("animation_data") + if npz_path is None: + logger.info("No animation data .npz available.") + return None + + import numpy as np + + data = np.load(npz_path) + return { + "fields": data["fields"], + "times": data["times"], + "eps_data": data["eps_data"], + "extent": data["extent"], + "plane": str(data["plane"]), + } + def plot( self, db: bool = True, @@ -366,3 +382,178 @@ def _is_reflection(name: str) -> bool: ), ) return fig + + +class CouplingResult(BaseModel): + """Coupling efficiency results from fiber-to-chip Gaussian beam simulation. + + Parses CSV output from the cloud runner and provides visualization. + """ + + model_config = ConfigDict(validate_assignment=True) + + wavelengths: list[float] = Field(default_factory=list) + coupling_efficiency: dict[str, list[float]] = Field( + default_factory=dict, + description="Port name -> CE values (linear, 0-1) per wavelength", + ) + debug_info: dict[str, Any] = Field( + default_factory=dict, + description="Debug diagnostics from meep_debug.json (if available)", + ) + diagnostic_images: dict[str, str] = Field( + default_factory=dict, + description="Diagnostic image paths: key -> filepath", + ) + + @classmethod + def from_csv(cls, path: str | Path) -> CouplingResult: + """Parse coupling efficiency results from CSV file. + + Expected CSV format: + wavelength,o1,o2,... + 1.5,0.25,0.01,... + + Values are linear coupling efficiency (0-1). + + Args: + path: Path to CSV file + + Returns: + CouplingResult instance + """ + path = Path(path) + wavelengths: list[float] = [] + ce: dict[str, list[float]] = {} + + with open(path) as f: + reader = csv.DictReader(f) + if reader.fieldnames is None: + return cls() + + port_names = [c for c in reader.fieldnames if c != "wavelength"] + for name in port_names: + ce[name] = [] + + for row in reader: + wavelengths.append(float(row["wavelength"])) + for name in port_names: + ce[name].append(float(row[name])) + + # Auto-load debug log if present alongside CSV + debug_info: dict[str, Any] = {} + debug_path = path.parent / "meep_debug.json" + if debug_path.exists(): + with contextlib.suppress(json.JSONDecodeError, OSError): + debug_info = json.loads(debug_path.read_text()) + + # Auto-detect diagnostic PNGs + diagnostic_images: dict[str, str] = {} + for key, filename in [ + ("geometry_xy", "meep_geometry_xy.png"), + ("geometry_xz", "meep_geometry_xz.png"), + ("geometry_yz", "meep_geometry_yz.png"), + ("fields_xy", "meep_fields_xy.png"), + ("animation", "meep_animation.mp4"), + ("animation_data", "meep_animation_data.npz"), + ]: + img_path = path.parent / filename + if img_path.exists(): + diagnostic_images[key] = str(img_path) + + return cls( + wavelengths=wavelengths, + coupling_efficiency=ce, + debug_info=debug_info, + diagnostic_images=diagnostic_images, + ) + + @property + def peak_ce(self) -> dict[str, float]: + """Peak coupling efficiency per port in dB. + + Returns: + Dict of port_name -> peak CE in dB. + """ + import math + + result: dict[str, float] = {} + for name, values in self.coupling_efficiency.items(): + if values: + peak_linear = max(values) + result[name] = ( + 10 * math.log10(peak_linear) if peak_linear > 0 else -100.0 + ) + else: + result[name] = -100.0 + return result + + def show_animation(self) -> None: + """Display field animation MP4 in Jupyter.""" + mp4_path = self.diagnostic_images.get("animation") + if mp4_path is None: + logger.info("No animation MP4 available.") + return + + from IPython.display import Video, display + + display(Video(mp4_path, embed=True, mimetype="video/mp4")) + + def load_animation_data(self) -> dict[str, Any] | None: + """Load raw animation frame data from the consolidated .npz. + + Returns: + Dict with keys ``fields`` (N, H, W), ``times`` (N,), + ``eps_data`` (H, W), ``extent`` (4,), ``plane`` (str), + or None if not available. + """ + npz_path = self.diagnostic_images.get("animation_data") + if npz_path is None: + logger.info("No animation data .npz available.") + return None + + import numpy as np + + data = np.load(npz_path) + return { + "fields": data["fields"], + "times": data["times"], + "eps_data": data["eps_data"], + "extent": data["extent"], + "plane": str(data["plane"]), + } + + def plot(self, db: bool = True, **kwargs: Any) -> Any: + """Plot coupling efficiency vs wavelength. + + Args: + db: If True, plot in dB scale (10*log10) + **kwargs: Passed to matplotlib plot() + + Returns: + matplotlib Figure + """ + import matplotlib.pyplot as plt + + fig, ax = plt.subplots() + plt.close(fig) # prevent double display in notebooks + + ylabel = "Coupling Efficiency (dB)" if db else "Coupling Efficiency" + for name, values in self.coupling_efficiency.items(): + if db: + import math + + y_vals = [10 * math.log10(v) if v > 0 else -100 for v in values] + else: + y_vals = values + + ax.plot(self.wavelengths, y_vals, ".-", label=name, **kwargs) + + ax.set_xlabel("Wavelength (um)") + ax.set_ylabel(ylabel) + ax.legend() + ax.grid(True, alpha=0.3) + ax.set_title("Coupling Efficiency") + fig.tight_layout() + + return fig diff --git a/src/gsim/meep/overlay.py b/src/gsim/meep/overlay.py index e0cf97e9..f8f9c612 100644 --- a/src/gsim/meep/overlay.py +++ b/src/gsim/meep/overlay.py @@ -55,6 +55,27 @@ class DielectricOverlay: zmax: float +@dataclass(frozen=True) +class FiberSourceOverlay: + """Fiber (Gaussian beam) source metadata for 2D overlay rendering. + + Attributes: + position: (x, y) beam center on the chip surface. + z_position: Absolute z-coordinate of the source plane. + beam_waist: Gaussian beam waist radius in um. + angle_theta: Polar angle from vertical in degrees. + angle_phi: Azimuthal angle in degrees (0 = tilted in XZ plane). + direction: "down" (into chip) or "up". + """ + + position: tuple[float, float] + z_position: float + beam_waist: float + angle_theta: float + angle_phi: float + direction: str + + @dataclass(frozen=True) class SimOverlay: """Simulation cell metadata for 2D visualization overlays. @@ -65,6 +86,7 @@ class SimOverlay: dpml: PML absorber thickness in um. ports: List of port overlays for rendering. dielectrics: List of background dielectric slabs for rendering. + fiber_source: Optional fiber source overlay. """ cell_min: tuple[float, float, float] @@ -72,6 +94,7 @@ class SimOverlay: dpml: float ports: list[PortOverlay] = field(default_factory=list) dielectrics: list[DielectricOverlay] = field(default_factory=list) + fiber_source: FiberSourceOverlay | None = None def build_sim_overlay( @@ -81,6 +104,7 @@ def build_sim_overlay( z_span: float | None = None, dielectrics: list[dict[str, Any]] | None = None, component_bbox: tuple[float, float, float, float] | None = None, + fiber_source: FiberSourceOverlay | None = None, ) -> SimOverlay: """Build a SimOverlay from geometry model, domain config, and port data. @@ -94,6 +118,7 @@ def build_sim_overlay( port extension. When provided, cell XY boundaries are computed from this bbox instead of the geometry model bbox (which may include extended waveguides). + fiber_source: Optional fiber source overlay metadata. Returns: SimOverlay with computed cell boundaries and port overlays. @@ -161,4 +186,5 @@ def build_sim_overlay( dpml=dpml, ports=ports, dielectrics=diel_overlays, + fiber_source=fiber_source, ) diff --git a/src/gsim/meep/ports.py b/src/gsim/meep/ports.py index 59897ddf..a66dcd4c 100644 --- a/src/gsim/meep/ports.py +++ b/src/gsim/meep/ports.py @@ -44,6 +44,8 @@ def extract_port_info( component: Component, layer_stack: LayerStack, source_port: str | None = None, + *, + mark_source: bool = True, ) -> list[PortData]: """Extract port information from a gdsfactory component. @@ -51,6 +53,8 @@ def extract_port_info( component: gdsfactory Component with ports layer_stack: LayerStack to determine z-coordinates source_port: Name of the source port. If None, first port is the source. + mark_source: If False, no port gets ``is_source=True`` (used when + the source is a FiberSource rather than a port eigenmode). Returns: List of PortData objects ready for JSON serialization @@ -62,7 +66,12 @@ def extract_port_info( for i, gf_port in enumerate(component.ports): normal_axis, direction = get_port_normal(gf_port.orientation) - is_source = gf_port.name == source_port if source_port is not None else i == 0 + if mark_source: + is_source = ( + gf_port.name == source_port if source_port is not None else i == 0 + ) + else: + is_source = False ports.append( PortData( diff --git a/src/gsim/meep/script.py b/src/gsim/meep/script.py index ab709cf4..882db396 100644 --- a/src/gsim/meep/script.py +++ b/src/gsim/meep/script.py @@ -491,6 +491,210 @@ def build_monitors(config, sim): return monitors +# --------------------------------------------------------------------------- +# Fiber source: Gaussian beam + incident flux +# --------------------------------------------------------------------------- + +def build_fiber_monitor(config, sim): + """Build a fiber mode monitor for reciprocal coupling measurement. + + Creates a horizontal mode monitor at the fiber position above (or + below) the grating. The runner will later call + ``get_eigenmode_coefficients`` on this monitor with an oblique + kpoint matching the fiber angle to extract the coupling coefficient. + + Returns: + meep ModeMonitor object + """ + fdtd = config["fdtd"] + fcen = fdtd["fcen"] + df = fdtd["df"] + nfreq = fdtd["num_freqs"] + src_cfg = config["source"] + position = src_cfg["position"] + z_position = src_cfg["z_position"] + domain = config["domain"] + dpml = domain["dpml"] + + sx = sim.cell_size.x - 2 * dpml + sy = sim.cell_size.y - 2 * dpml + + center = mp.Vector3(position[0], position[1], z_position) + size = mp.Vector3(sx, sy, 0) + + fiber_mon = sim.add_mode_monitor( + fcen, df, nfreq, + mp.ModeRegion(center=center, size=size), + ) + return fiber_mon + + +def _get_n_at_z(config, z): + """Find the refractive index of the dielectric region at height z. + + Walks through the dielectrics list in the config and returns + the refractive index of the first region that contains z. + Falls back to 1.0 (air) if no region contains z. + """ + materials = config.get("materials", {}) + for diel in config.get("dielectrics", []): + if diel["zmin"] <= z <= diel["zmax"]: + mat_name = diel["material"] + mat = materials.get(mat_name) + if mat and "refractive_index" in mat: + return mat["refractive_index"] + return 1.0 + + +def extract_fiber_coupling(config, sim, monitors, fiber_monitor): + """Extract fiber coupling via reciprocal eigenmode decomposition. + + An EigenModeSource at the waveguide port launches light toward the + grating. This function measures: + + 1. The incident eigenmode coefficient at the source port (for + normalization). + 2. The coupling into the fiber mode at the fiber monitor plane + using ``get_eigenmode_coefficients`` with ``direction=NO_DIRECTION`` + and an oblique kpoint matching the fiber angle (gplugins approach). + + CE = |alpha_fiber|^2 / |alpha_incident|^2 + + Returns: + (ce_params, debug_data) tuple + """ + fdtd = config["fdtd"] + fcen = fdtd["fcen"] + nfreq = fdtd["num_freqs"] + df = fdtd["df"] + freqs = np.linspace(fcen - df / 2, fcen + df / 2, nfreq) + + src_cfg = config["source"] + theta_deg = src_cfg["angle_theta"] + phi_deg = src_cfg["angle_phi"] + direction = src_cfg["direction"] + z_position = src_cfg["z_position"] + + # --- Incident coefficient at waveguide source port --- + source_port = None + for p in config["ports"]: + if p["is_source"]: + source_port = p + break + if source_port is None: + logger.error("No source port found for fiber coupling extraction") + return {}, {} + + src_kp = _port_kpoint(source_port) + src_ob = sim.get_eigenmode_coefficients( + monitors[source_port["name"]], [1], eig_parity=mp.NO_PARITY, + kpoint_func=lambda f, n, kp=src_kp: kp, + ) + src_dir = source_port["direction"] + # Incoming direction for the source port + in_idx = 0 if src_dir == "+" else 1 + incident_coeffs = src_ob.alpha[0, :, in_idx] + + # --- Fiber mode kpoint --- + theta = math.radians(theta_deg) + phi = math.radians(phi_deg) + n_clad = _get_n_at_z(config, z_position) + logger.info("Fiber monitor n_clad=%.3f at z=%.3f", n_clad, z_position) + + # k-direction for the fiber mode (unit vector) + if direction == "down": + kx = math.sin(theta) * math.cos(phi) + ky = math.sin(theta) * math.sin(phi) + kz = -math.cos(theta) + else: + kx = math.sin(theta) * math.cos(phi) + ky = math.sin(theta) * math.sin(phi) + kz = math.cos(theta) + + # Scale kpoint by fcen * n_clad (gplugins convention) + fiber_kpoint = mp.Vector3(kx, ky, kz) * (fcen * n_clad) + + # --- Eigenmode decomposition at fiber monitor --- + fiber_ob = sim.get_eigenmode_coefficients( + fiber_monitor, [1], + eig_parity=mp.NO_PARITY, + direction=mp.NO_DIRECTION, + kpoint_func=lambda f, n, kp=fiber_kpoint: kp, + ) + + # With NO_DIRECTION, pick the coefficient with larger magnitude + a0 = fiber_ob.alpha[0, :, 0] + a1 = fiber_ob.alpha[0, :, 1] + fiber_alpha = np.where(np.abs(a0) >= np.abs(a1), a0, a1) + + debug_data = { + "incident_coefficients": { + "port": source_port["name"], + "alpha_mag": [float(abs(c)) for c in incident_coeffs], + }, + "fiber_coefficients": { + "alpha_0_mag": [float(abs(a)) for a in a0], + "alpha_1_mag": [float(abs(a)) for a in a1], + "alpha_used_mag": [float(abs(a)) for a in fiber_alpha], + "n_clad": n_clad, + "kpoint": [float(fiber_kpoint.x), float(fiber_kpoint.y), + float(fiber_kpoint.z)], + }, + "eigenmode_info": {}, + } + + # --- CE = |alpha_fiber|^2 / |alpha_incident|^2 --- + ce_params = {} + ce = np.zeros(nfreq) + for fi in range(nfreq): + if abs(incident_coeffs[fi]) > 0: + ce[fi] = float( + abs(fiber_alpha[fi]) ** 2 / abs(incident_coeffs[fi]) ** 2 + ) + ce_params["fiber"] = ce + + debug_data["eigenmode_info"]["fiber"] = { + "alpha_mag": [float(abs(a)) for a in fiber_alpha], + "ce": [float(c) for c in ce], + } + + logger.info( + "Fiber CE at center freq: %.4e (%.2f dB)", + ce[nfreq // 2], + 10 * np.log10(max(ce[nfreq // 2], 1e-30)), + ) + + return ce_params, debug_data + + +def save_ce_results(config, ce_params, output_path="coupling_efficiency.csv"): + """Save coupling efficiency to CSV. Only rank 0 writes.""" + if not mp.am_master(): + return + fdtd = config["fdtd"] + fcen = fdtd["fcen"] + df = fdtd["df"] + nfreq = fdtd["num_freqs"] + + freqs = np.linspace(fcen - df / 2, fcen + df / 2, nfreq) + wavelengths = 1.0 / freqs + + port_names = sorted(ce_params.keys()) + fieldnames = ["wavelength"] + port_names + + with open(output_path, "w", newline="") as f: + writer = csv.DictWriter(f, fieldnames=fieldnames) + writer.writeheader() + + for i, wl in enumerate(wavelengths): + row = {"wavelength": f"{wl:.6f}"} + for name in port_names: + row[name] = f"{float(ce_params[name][i]):.8f}" + writer.writerow(row) + + logger.info("Coupling efficiency saved to %s", output_path) + + # --------------------------------------------------------------------------- # S-parameter extraction # --------------------------------------------------------------------------- @@ -699,6 +903,8 @@ def save_debug_log(config, s_params, debug_data, wall_seconds=0.0, "stopping_mode": stopping["mode"], }, "eigenmode_info": debug_data.get("eigenmode_info", {}), + "incident_flux": debug_data.get("incident_flux", []), + "incident_flux_raw": debug_data.get("incident_flux_raw", []), "incident_coefficients": debug_data.get("incident_coefficients", {}), "raw_coefficients": debug_data.get("raw_coefficients", {}), "power_conservation": debug_data.get("power_conservation", []), @@ -817,29 +1023,139 @@ def save_field_snapshot(sim, config, cell_center): logger.warning("Field snapshot failed: %s", e) -def save_animation_field(sim, xy_plane, frame_counter): - """Save raw 2D field data during time-stepping (no plotting). +def save_animation_field(sim, anim_plane, frame_counter, plane_id="xy"): + """Save animation frame via MEEP's native HDF5 output (MPI-safe). - Saves a compressed .npz with the Ey field array and timestamp. - Rendering with a globally fixed colorbar happens after sim.run(). + Uses ``mp.output_efield_y`` with ``mp.in_volume`` which is + properly MPI-collective without the ``get_array`` vol-mismatch + bug on cross-section planes. HDF5 files are rendered to PNG + after ``sim.run()`` completes. Returns: Incremented frame counter. """ - field_data = np.real(sim.get_array(vol=xy_plane, component=mp.Ey)) + os.makedirs("frames", exist_ok=True) + # MEEP's output functions handle MPI natively + mp.output_efield_y(sim, mp.Volume( + center=anim_plane.center, size=anim_plane.size, + )) + # MEEP writes ey-TTTTTT.TT.h5 — rename to our naming convention t = sim.meep_time() if mp.am_master(): - os.makedirs("frames", exist_ok=True) - np.savez_compressed( - f"frames/meep_field_{frame_counter:04d}.npz", - field=field_data, - time=t, - ) + # Find the most recently written h5 file + import glob + h5_files = sorted(glob.glob("ey-*.h5"), key=os.path.getmtime) + if h5_files: + latest = h5_files[-1] + dest = f"frames/field_{frame_counter:04d}.h5" + os.rename(latest, dest) return frame_counter + 1 -def render_animation_frames(eps_data, extent): - """Render saved field .npz files into PNGs with fixed global colorbar. +def render_h5_animation_frames(plane_id="xy"): + """Render HDF5 field data to PNG frames after simulation completes. + + Reads per-step .h5 files saved by ``save_animation_field`` and + renders them with a globally consistent colorbar. + """ + if not HAS_MATPLOTLIB or not mp.am_master(): + return + import glob + + try: + import h5py + except ImportError: + logger.warning("h5py not available — skipping animation frame rendering") + return + + h5_files = sorted(glob.glob("frames/field_*.h5")) + if not h5_files: + logger.warning("No HDF5 field data found for animation") + return + + # First pass: find global min/max for consistent colorbar + vmax = 0 + for path in h5_files: + with h5py.File(path, "r") as hf: + key = list(hf.keys())[0] + data = hf[key][:] + vmax = max(vmax, np.max(np.abs(data))) + + if vmax == 0: + vmax = 1.0 + + # Second pass: render frames + for i, path in enumerate(h5_files): + with h5py.File(path, "r") as hf: + key = list(hf.keys())[0] + data = hf[key][:] + + fig, ax = plt.subplots(figsize=(12, 4)) + im = ax.imshow( + data.T if data.ndim == 2 else data.squeeze().T, + origin="lower", cmap="RdBu_r", + vmin=-vmax, vmax=vmax, aspect="auto", + ) + ax.set_title(f"Ey ({plane_id.upper()}) frame {i}") + plt.colorbar(im, ax=ax, label="Ey") + fig.savefig( + f"frames/meep_frame_{i:04d}.png", + dpi=100, bbox_inches="tight", + ) + plt.close(fig) + + logger.info("Rendered %d animation frames from HDF5 data", len(h5_files)) + + +def consolidate_animation_data(eps_data, extent, plane="xy"): + """Consolidate per-frame .npz files into a single meep_animation_data.npz. + + Loads all per-frame field snapshots, bundles them with the epsilon + background and extent into one file, then removes the per-frame files. + + Returns: + Path to the consolidated file, or None if no frames found. + """ + import glob + + npz_files = sorted(glob.glob("frames/meep_field_*.npz")) + if not npz_files: + logger.warning("No field data files found to consolidate") + return None + + fields = [] + times = [] + for path in npz_files: + d = np.load(path) + fields.append(d["field"]) + times.append(float(d["time"])) + + np.savez_compressed( + "meep_animation_data.npz", + fields=np.array(fields), + times=np.array(times), + eps_data=eps_data, + extent=np.array(extent), + plane=plane, + ) + + # Clean up per-frame intermediates + for path in npz_files: + os.remove(path) + # Remove frames dir if empty + try: + os.rmdir("frames") + except OSError: + pass + + logger.info( + "Consolidated %d frames into meep_animation_data.npz", len(npz_files) + ) + return "meep_animation_data.npz" + + +def render_animation_frames(anim_data_path="meep_animation_data.npz"): + """Render frames from consolidated .npz into PNGs for MP4 compilation. Two-pass: first finds the global field maximum across all frames, then renders every frame with the same vmin/vmax so field decay is @@ -847,47 +1163,48 @@ def render_animation_frames(eps_data, extent): Call only on master rank after sim.run(). """ - import glob - if not HAS_MATPLOTLIB: logger.warning( - "matplotlib not available, .npz field files kept but not rendered" + "matplotlib not available, animation data kept in .npz" ) return - npz_files = sorted(glob.glob("frames/meep_field_*.npz")) - if not npz_files: - logger.warning("No field data files found to render") + if not os.path.exists(anim_data_path): + logger.warning("No animation data file found: %s", anim_data_path) return - # Pass 1 — global max - global_max = 0.0 - for path in npz_files: - d = np.load(path) - global_max = max(global_max, float(np.max(np.abs(d["field"])))) + data = np.load(anim_data_path) + fields = data["fields"] + times = data["times"] + eps_data = data["eps_data"] + extent = data["extent"].tolist() + plane = str(data["plane"]) + + n_frames = len(fields) + global_max = float(np.max(np.abs(fields))) if global_max == 0: global_max = 1.0 logger.info( "Rendering %d frames (Ey global max = %.4g) ...", - len(npz_files), global_max, + n_frames, global_max, ) from mpl_toolkits.axes_grid1 import make_axes_locatable - # Pass 2 — render each frame - for i, path in enumerate(npz_files): - d = np.load(path) - field = d["field"] - t = float(d["time"]) + xlabel = "x (um)" + ylabel = "z (um)" if plane == "xz" else "y (um)" + + os.makedirs("frames", exist_ok=True) + for i in range(n_frames): + field = fields[i] + t = float(times[i]) fig, ax = plt.subplots(1, 1, figsize=(5, 4)) - # Epsilon background (grayscale) ax.imshow( eps_data.T, origin="lower", extent=extent, cmap="binary", interpolation="none", ) - # Field overlay with fixed colorbar im = ax.imshow( field.T, origin="lower", extent=extent, cmap="RdBu", interpolation="spline36", alpha=0.8, @@ -896,27 +1213,25 @@ def render_animation_frames(eps_data, extent): divider = make_axes_locatable(ax) cax = divider.append_axes("right", size="4%", pad=0.06) fig.colorbar(im, cax=cax, label="Ey") - ax.set_title(f"Ey t={t:.2f}") - ax.set_xlabel("x (um)") - ax.set_ylabel("y (um)") + ax.set_title(f"Ey ({plane.upper()}) t={t:.2f}") + ax.set_xlabel(xlabel) + ax.set_ylabel(ylabel) fig.tight_layout() fig.savefig(f"frames/meep_frame_{i:04d}.png", dpi=150) plt.close(fig) - # Clean up .npz intermediates - for path in npz_files: - os.remove(path) - - logger.info("Rendered %d frames with fixed colorbar", len(npz_files)) + logger.info("Rendered %d frames with fixed colorbar", n_frames) def compile_animation_mp4(fps=15): """Stitch meep_frame_*.png into meep_animation.mp4 via ffmpeg. - Falls back gracefully if ffmpeg is not available — frame PNGs - are still kept as individual files. + On success, removes the intermediate frame PNGs (raw data is + preserved in meep_animation_data.npz). Falls back gracefully + if ffmpeg is not available. """ import glob + import shutil import subprocess frames = sorted(glob.glob("frames/meep_frame_*.png")) @@ -939,11 +1254,13 @@ def compile_animation_mp4(fps=15): capture_output=True, ) logger.info("Saved meep_animation.mp4") + # Clean up intermediate PNGs — raw data lives in .npz + shutil.rmtree("frames", ignore_errors=True) except FileNotFoundError: - logger.warning("ffmpeg not found — frame PNGs saved but MP4 not created") + logger.warning("ffmpeg not found — frame PNGs kept in frames/") except subprocess.CalledProcessError as e: logger.warning("ffmpeg failed: %s", e.stderr.decode()[:500]) - logger.info("Frame PNGs are still available in frames/") + logger.info("Frame PNGs kept in frames/") def save_epsilon_raw(sim, config, cell_center): @@ -991,12 +1308,8 @@ def main(): # Background slabs first, then patterned prisms (later objects take precedence) geometry = background_slabs + geometry - logger.info("Building sources...") - sources = build_sources(config) - - if not sources: - logger.error("No source port found in config") - sys.exit(1) + source_type = config["source"].get("source_type", "mode") + logger.info("Source type: %s", source_type) resolution = config["resolution"]["pixels_per_um"] fdtd = config["fdtd"] @@ -1025,6 +1338,15 @@ def main(): dpml = domain["dpml"] margin_xy = domain["margin_xy"] + # For fiber source, ensure the cell extends to include the fiber monitor z + if source_type == "fiber": + mon_z = config["source"]["z_position"] + mon_dir = config["source"]["direction"] + if mon_dir == "down" and mon_z > z_max: + z_max = mon_z + 0.5 # 0.5 um margin above fiber monitor + elif mon_dir == "up" and mon_z < z_min: + z_min = mon_z - 0.5 + # XY: margin_xy is gap between geometry bbox and PML cell_x = (bbox_right - bbox_left) + 2 * (margin_xy + dpml) cell_y = (bbox_top - bbox_bottom) + 2 * (margin_xy + dpml) @@ -1037,6 +1359,14 @@ def main(): (z_max + z_min) / 2, ) + # Build eigenmode sources — same for both mode and fiber (reciprocal) paths + logger.info("Building eigenmode sources...") + sources = build_sources(config) + + if not sources: + logger.error("No source found in config") + sys.exit(1) + logger.info("Cell size: %.2f x %.2f x %.2f um", cell_x, cell_y, cell_z) logger.info("PML: %.2f um, margin_xy: %.2f", dpml, margin_xy) logger.info("Resolution: %s pixels/um", resolution) @@ -1079,6 +1409,9 @@ def main(): spx_tol = accuracy["subpixel_tol"] if spx_tol != 1e-4: sim_kwargs["subpixel_tol"] = spx_tol + # Suppress meep's C++ geometry dump (vertex lists for every prism). + # All important info is already logged by this script. + mp.verbosity(0) sim = mp.Simulation(**sim_kwargs) # --- Diagnostics & preview mode --- @@ -1104,6 +1437,12 @@ def main(): logger.info("Building monitors...") monitors = build_monitors(config, sim) + # For fiber source, build the fiber mode monitor (reciprocal method) + fiber_monitor = None + if source_type == "fiber": + logger.info("Building fiber mode monitor (reciprocal)...") + fiber_monitor = build_fiber_monitor(config, sim) + stopping = config["stopping"] run_after = stopping["run_after_sources"] @@ -1133,24 +1472,50 @@ def _verbose_print(sim_obj): _frame_counter = [0] # mutable container for closure _anim_plane = None + _anim_plane_id = diagnostics.get("animation_plane", "xy") if diag_animation: z_min_anim = min(l["zmin"] for l in config["layer_stack"]) z_max_anim = max(l["zmax"] for l in config["layer_stack"]) z_core_anim = (z_min_anim + z_max_anim) / 2 - _anim_plane = mp.Volume( - center=mp.Vector3(cell_center.x, cell_center.y, z_core_anim), - size=mp.Vector3(sim.cell_size.x, sim.cell_size.y, 0), - ) + + # Snap coordinates to grid and shrink slightly to avoid MPI + # chunk boundary issues with get_array on cross-section planes. + dx = 1.0 / resolution + def _snap(val): + return round(val * resolution) / resolution + + if _anim_plane_id == "xz": + _anim_plane = mp.Volume( + center=mp.Vector3( + _snap(cell_center.x), _snap(cell_center.y), + _snap(cell_center.z), + ), + size=mp.Vector3( + sim.cell_size.x - 2 * dx, 0, + sim.cell_size.z - 2 * dx, + ), + ) + else: + _anim_plane = mp.Volume( + center=mp.Vector3( + _snap(cell_center.x), _snap(cell_center.y), + _snap(z_core_anim), + ), + size=mp.Vector3( + sim.cell_size.x - 2 * dx, + sim.cell_size.y - 2 * dx, 0, + ), + ) def _capture_frame(sim_obj): _frame_counter[0] = save_animation_field( - sim_obj, _anim_plane, _frame_counter[0] + sim_obj, _anim_plane, _frame_counter[0], _anim_plane_id ) step_funcs.append(mp.at_every(animation_interval, _capture_frame)) logger.info( - "Animation: saving field data every %s time units", - animation_interval, + "Animation (%s plane): saving field data every %s time units", + _anim_plane_id, animation_interval, ) stop_mode = stopping["mode"] @@ -1227,29 +1592,33 @@ def _capture_frame(sim_obj): if diag_fields: save_field_snapshot(sim, config, cell_center) - if diag_animation and _frame_counter[0] > 0 and _anim_plane is not None: - # get_array is collective — all ranks must call - eps_data = sim.get_array(vol=_anim_plane, component=mp.Dielectric) + if diag_animation and _frame_counter[0] > 0: + # Render HDF5 field data to PNG frames, then compile MP4 + render_h5_animation_frames(plane=_anim_plane_id) if mp.am_master(): - _ctr = _anim_plane.center - _sz = _anim_plane.size - _extent = [ - _ctr.x - _sz.x / 2, _ctr.x + _sz.x / 2, - _ctr.y - _sz.y / 2, _ctr.y + _sz.y / 2, - ] - render_animation_frames(eps_data, _extent) compile_animation_mp4() - logger.info("Extracting S-parameters...") - s_params, debug_data = extract_s_params(config, sim, monitors) + if source_type == "fiber" and fiber_monitor is not None: + logger.info("Extracting fiber coupling (reciprocal method)...") + ce_params, debug_data = extract_fiber_coupling( + config, sim, monitors, fiber_monitor + ) + debug_data["_meep_time"] = sim.meep_time() + debug_data["_timesteps"] = sim.timestep() + debug_data["_cell_size"] = [cell_x, cell_y, cell_z] + + save_ce_results(config, ce_params) + save_debug_log(config, {}, debug_data, wall_seconds=wall_seconds) + else: + logger.info("Extracting S-parameters...") + s_params, debug_data = extract_s_params(config, sim, monitors) + debug_data["_meep_time"] = sim.meep_time() + debug_data["_timesteps"] = sim.timestep() + debug_data["_cell_size"] = [cell_x, cell_y, cell_z] - # Attach simulation metadata to debug_data - debug_data["_meep_time"] = sim.meep_time() - debug_data["_timesteps"] = sim.timestep() - debug_data["_cell_size"] = [cell_x, cell_y, cell_z] + save_results(config, s_params) + save_debug_log(config, s_params, debug_data, wall_seconds=wall_seconds) - save_results(config, s_params) - save_debug_log(config, s_params, debug_data, wall_seconds=wall_seconds) logger.info("Done!") diff --git a/src/gsim/meep/simulation.py b/src/gsim/meep/simulation.py index 871cb688..b8088038 100644 --- a/src/gsim/meep/simulation.py +++ b/src/gsim/meep/simulation.py @@ -16,6 +16,7 @@ from gsim.meep.models.api import ( FDTD, Domain, + FiberSource, Geometry, Material, ModeSource, @@ -77,7 +78,7 @@ class Simulation(BaseModel): geometry: Geometry = Field(default_factory=Geometry) materials: dict[str, float | Material] = Field(default_factory=dict) - source: ModeSource = Field(default_factory=ModeSource) + source: ModeSource | FiberSource = Field(default_factory=ModeSource) monitors: list[str] = Field(default_factory=list) domain: Domain = Field(default_factory=Domain) solver: FDTD = Field(default_factory=FDTD) @@ -149,7 +150,7 @@ def validate_config(self) -> Any: ports = list(self.geometry.component.ports) if not ports: errors.append("Component has no ports.") - elif self.source.port is not None: + elif isinstance(self.source, ModeSource) and self.source.port is not None: port_names = [p.name for p in ports] if self.source.port not in port_names: errors.append( @@ -192,6 +193,18 @@ def validate_config(self) -> Any: elif s.stopping == "fixed": warnings_list.append(f"Stopping: fixed (time={s.max_time})") + # Fiber source: warn if z_offset is large relative to margin_z_above + if ( + isinstance(self.source, FiberSource) + and self.source.z_offset >= self.domain.margin_z_above + ): + warnings_list.append( + f"FiberSource z_offset ({self.source.z_offset}) >= " + f"margin_z_above ({self.domain.margin_z_above}). " + f"Consider increasing margin_z_above to ensure the " + f"source is well within the simulation cell." + ) + return ValidationResult( valid=len(errors) == 0, errors=errors, warnings=warnings_list ) @@ -313,11 +326,68 @@ def _source_config(self) -> Any: """Translate ModeSource → SourceConfig.""" from gsim.meep.models.config import SourceConfig + if not isinstance(self.source, ModeSource): + raise TypeError("_source_config() requires ModeSource") return SourceConfig( bandwidth=None, port=self.source.port, ) + def _fiber_source_config(self, stack: Any) -> Any: + """Translate FiberSource → FiberSourceConfig (reciprocal method). + + Resolves fiber monitor position and the waveguide port for + the EigenModeSource. + """ + from gsim.meep.models.config import FiberSourceConfig + + if not isinstance(self.source, FiberSource): + raise TypeError("_fiber_source_config() requires FiberSource") + + src = self.source + + # Resolve position: auto = component bbox center + if src.position is not None: + position = list(src.position) + else: + bbox = self.geometry.component.dbbox() + position = [ + (bbox.left + bbox.right) / 2.0, + (bbox.bottom + bbox.top) / 2.0, + ] + + # Compute absolute z from core layer top + offset. + from gsim.meep.ports import _find_highest_n_layer + + core_layer, _ = _find_highest_n_layer(stack) + if core_layer is not None: + z_ref_top = core_layer.zmax + z_ref_bottom = core_layer.zmin + else: + z_ref_top = max(layer.zmax for layer in stack.layers.values()) + z_ref_bottom = min(layer.zmin for layer in stack.layers.values()) + + if src.direction == "down": + z_position = z_ref_top + src.z_offset + else: + z_position = z_ref_bottom - src.z_offset + + # Compute fwidth (same as ModeSource — eigenmode source) + wl_cfg = self._wavelength_config() + fwidth = max(3 * wl_cfg.df, 0.2 * wl_cfg.fcen) + + return FiberSourceConfig( + port=src.port, + beam_waist=src.beam_waist, + angle_theta=src.angle_theta, + angle_phi=src.angle_phi, + polarization=src.polarization, + direction=src.direction, + position=position, + z_position=z_position, + fwidth=fwidth, + ) + def _stopping_config(self) -> Any: """Translate FDTD stopping fields → StoppingConfig.""" from gsim.meep.models.config import StoppingConfig @@ -376,6 +446,7 @@ def _diagnostics_config(self) -> Any: save_epsilon_raw=self.solver.save_epsilon_raw, save_animation=self.solver.save_animation, animation_interval=self.solver.animation_interval, + animation_plane=self.solver.animation_plane, preview_only=self.solver.preview_only, verbose_interval=self.solver.verbose_interval, ) @@ -436,7 +507,11 @@ def build_config(self) -> BuildResult: # Build config objects domain_cfg = self._domain_config() wl_cfg = self._wavelength_config() - source_cfg = self._source_config() + is_fiber = isinstance(self.source, FiberSource) + if is_fiber: + source_cfg = self._fiber_source_config(stack) + else: + source_cfg = self._source_config() stopping_cfg = self._stopping_config() resolution_cfg = self._resolution_config() accuracy_cfg = self._accuracy_config() @@ -487,10 +562,21 @@ def build_config(self) -> BuildResult: ) used_materials.add(diel["material"]) - # Extract port info from original component - port_infos = extract_port_info( - original_component, stack, source_port=source_cfg.port - ) + # Extract port info from original component. + # For fiber (reciprocal): the waveguide port IS the EigenModeSource. + if is_fiber: + fiber_port = source_cfg.port + port_infos = extract_port_info( + original_component, stack, source_port=fiber_port + ) + # Resolve auto-selected port name back into config + if fiber_port is None and port_infos: + resolved_port = next((p.name for p in port_infos if p.is_source), None) + source_cfg = source_cfg.model_copy(update={"port": resolved_port}) + else: + port_infos = extract_port_info( + original_component, stack, source_port=source_cfg.port + ) # Resolve materials material_data = resolve_materials( @@ -498,8 +584,11 @@ def build_config(self) -> BuildResult: ) # Compute source fwidth - fwidth = source_cfg.compute_fwidth(wl_cfg.fcen, wl_cfg.df) - source_for_config = source_cfg.model_copy(update={"fwidth": fwidth}) + if is_fiber: + source_for_config = source_cfg + else: + fwidth = source_cfg.compute_fwidth(wl_cfg.fcen, wl_cfg.df) + source_for_config = source_cfg.model_copy(update={"fwidth": fwidth}) # Translate domain.symmetries → SymmetryEntry for config symmetry_entries = [ @@ -684,8 +773,8 @@ def run( If ``False``, upload + start and return the ``job_id``. Returns: - ``SParameterResult`` when ``wait=True``, or ``job_id`` string - when ``wait=False``. + ``SParameterResult`` or ``CouplingResult`` when ``wait=True``, + or ``job_id`` string when ``wait=False``. """ self.upload(verbose=False) self.start(verbose=verbose != "quiet") @@ -709,14 +798,30 @@ def plot_2d(self, **kwargs: Any) -> Any: result = self.build_config() + # Build fiber source overlay if applicable + fiber_overlay = None + if isinstance(self.source, FiberSource): + from gsim.meep.overlay import FiberSourceOverlay + + src_cfg = result.config.source + fiber_overlay = FiberSourceOverlay( + position=(src_cfg.position[0], src_cfg.position[1]), + z_position=src_cfg.z_position, + beam_waist=src_cfg.beam_waist, + angle_theta=src_cfg.angle_theta, + angle_phi=src_cfg.angle_phi, + direction=src_cfg.direction, + ) + return plot_2d( component=result.component, stack=self.geometry.stack, domain_config=result.config.domain, - source_port=result.config.source.port, + source_port=getattr(result.config.source, "port", None), extend_ports_length=0, port_data=result.config.ports, component_bbox=result.config.component_bbox, + fiber_source=fiber_overlay, **kwargs, ) diff --git a/src/gsim/meep/viz.py b/src/gsim/meep/viz.py index 0781b9ea..f9448af8 100644 --- a/src/gsim/meep/viz.py +++ b/src/gsim/meep/viz.py @@ -165,6 +165,7 @@ def build_overlay( source_port: str | None = None, port_data: list | None = None, component_bbox: list[float] | tuple[float, ...] | None = None, + fiber_source: Any | None = None, ) -> Any: """Build a SimOverlay from config, if stack is available. @@ -179,6 +180,7 @@ def build_overlay( component_bbox: Original component bbox ``[xmin, ymin, xmax, ymax]`` from :meth:`Simulation.build_config`. When provided, cell boundaries are computed from this instead of ``component.dbbox()``. + fiber_source: Optional FiberSourceOverlay for rendering. Returns: SimOverlay or None if stack isn't configured. @@ -217,6 +219,7 @@ def build_overlay( port_data, dielectrics=dielectrics, component_bbox=orig_bbox, + fiber_source=fiber_source, ) @@ -273,6 +276,7 @@ def plot_2d( extend_ports_length: float | None = None, port_data: list | None = None, component_bbox: list[float] | tuple[float, ...] | None = None, + fiber_source: Any | None = None, ) -> plt.Axes | None: """Plot 2D cross-sections of the MEEP geometry. @@ -309,5 +313,6 @@ def plot_2d( source_port, port_data=port_data, component_bbox=component_bbox, + fiber_source=fiber_source, ) return plot_prism_slices(gm, x, y, z, ax, legend, slices, overlay=overlay) diff --git a/tests/meep/test_meep_models.py b/tests/meep/test_meep_models.py index 9d298bda..a02afb5e 100644 --- a/tests/meep/test_meep_models.py +++ b/tests/meep/test_meep_models.py @@ -9,7 +9,9 @@ from pydantic import ValidationError from gsim.meep import ( + CouplingResult, DomainConfig, + FiberSourceConfig, ResolutionConfig, SimConfig, SourceConfig, @@ -1410,3 +1412,307 @@ def test_plot_with_dielectric_overlay(self): assert "oxide" in patch_labels plt.close(fig) + + +# --------------------------------------------------------------------------- +# FiberSourceConfig tests +# --------------------------------------------------------------------------- + + +class TestFiberSourceConfig: + """Test FiberSourceConfig model.""" + + def test_creation(self): + cfg = FiberSourceConfig( + beam_waist=5.2, + angle_theta=10.0, + angle_phi=0.0, + polarization="TE", + direction="down", + position=[0.0, 0.0], + z_position=1.5, + fwidth=0.13, + ) + assert cfg.source_type == "fiber" + assert cfg.beam_waist == 5.2 + assert cfg.angle_theta == 10.0 + assert cfg.z_position == 1.5 + assert cfg.fwidth == 0.13 + + def test_model_dump(self): + cfg = FiberSourceConfig( + beam_waist=5.2, + angle_theta=10.0, + angle_phi=0.0, + polarization="TE", + direction="down", + position=[0.0, 0.0], + z_position=1.5, + fwidth=0.13, + ) + d = cfg.model_dump() + assert d["source_type"] == "fiber" + assert d["beam_waist"] == 5.2 + assert d["z_position"] == 1.5 + + def test_source_type_discriminator(self): + """SourceConfig and FiberSourceConfig have distinct source_type.""" + mode_cfg = SourceConfig() + assert mode_cfg.source_type == "mode" + + fiber_cfg = FiberSourceConfig( + beam_waist=5.2, + angle_theta=10.0, + angle_phi=0.0, + polarization="TE", + direction="down", + position=[0.0, 0.0], + z_position=1.5, + fwidth=0.13, + ) + assert fiber_cfg.source_type == "fiber" + + +# --------------------------------------------------------------------------- +# CouplingResult tests +# --------------------------------------------------------------------------- + + +class TestCouplingResult: + """Test CouplingResult model.""" + + def test_empty_result(self): + result = CouplingResult() + assert result.wavelengths == [] + assert result.coupling_efficiency == {} + + def test_from_csv(self, tmp_path): + csv_path = tmp_path / "coupling_efficiency.csv" + csv_path.write_text( + "wavelength,o1,o2\n" + "1.500000,0.25000000,0.01000000\n" + "1.550000,0.30000000,0.00500000\n" + ) + + result = CouplingResult.from_csv(csv_path) + assert len(result.wavelengths) == 2 + assert "o1" in result.coupling_efficiency + assert "o2" in result.coupling_efficiency + assert len(result.coupling_efficiency["o1"]) == 2 + assert abs(result.coupling_efficiency["o1"][0] - 0.25) < 1e-6 + assert abs(result.coupling_efficiency["o2"][1] - 0.005) < 1e-6 + + def test_peak_ce(self): + import math + + result = CouplingResult( + wavelengths=[1.5, 1.55], + coupling_efficiency={ + "o1": [0.25, 0.30], + "o2": [0.01, 0.005], + }, + ) + peaks = result.peak_ce + assert abs(peaks["o1"] - 10 * math.log10(0.30)) < 1e-6 + assert abs(peaks["o2"] - 10 * math.log10(0.01)) < 1e-6 + + def test_plot(self): + import matplotlib as mpl + + mpl.use("Agg") + + result = CouplingResult( + wavelengths=[1.5, 1.55], + coupling_efficiency={"o1": [0.25, 0.30]}, + ) + fig = result.plot(db=True) + assert fig is not None + + def test_plot_linear(self): + import matplotlib as mpl + + mpl.use("Agg") + + result = CouplingResult( + wavelengths=[1.5, 1.55], + coupling_efficiency={"o1": [0.25, 0.30]}, + ) + fig = result.plot(db=False) + assert fig is not None + + +# --------------------------------------------------------------------------- +# Script fiber source tests +# --------------------------------------------------------------------------- + + +class TestScriptFiberSource: + """Test that the runner script includes reciprocal fiber coupling support.""" + + def test_script_has_fiber_monitor_builder(self): + from gsim.meep.script import generate_meep_script + + script = generate_meep_script() + assert "build_fiber_monitor" in script + + def test_script_has_fiber_coupling_extractor(self): + from gsim.meep.script import generate_meep_script + + script = generate_meep_script() + assert "extract_fiber_coupling" in script + assert "save_ce_results" in script + + def test_script_has_n_at_z_helper(self): + from gsim.meep.script import generate_meep_script + + script = generate_meep_script() + assert "_get_n_at_z" in script + + def test_script_branches_on_source_type(self): + from gsim.meep.script import generate_meep_script + + script = generate_meep_script() + assert "source_type" in script + assert '"fiber"' in script + + def test_script_still_valid_python(self): + from gsim.meep.script import generate_meep_script + + script = generate_meep_script() + ast.parse(script) + + +# --------------------------------------------------------------------------- +# SimConfig with FiberSourceConfig tests +# --------------------------------------------------------------------------- + + +class TestSimConfigFiberSource: + """Test SimConfig accepts FiberSourceConfig via discriminated union.""" + + def test_json_roundtrip_fiber_source(self, tmp_path): + from gsim.meep.models.config import AccuracyConfig, DiagnosticsConfig + + fiber_cfg = FiberSourceConfig( + beam_waist=5.2, + angle_theta=10.0, + angle_phi=0.0, + polarization="TE", + direction="down", + position=[0.0, 0.0], + z_position=1.5, + fwidth=0.13, + ) + + cfg = SimConfig( + gds_filename="layout.gds", + verbose_interval=0, + layer_stack=[], + dielectrics=[], + ports=[], + materials={}, + wavelength=WavelengthConfig(wavelength=1.55, bandwidth=0.1, num_freqs=11), + source=fiber_cfg, + stopping=StoppingConfig( + mode="fixed", + max_time=100.0, + decay_dt=50.0, + decay_component="Ey", + threshold=0.05, + dft_min_run_time=100, + ), + resolution=ResolutionConfig(pixels_per_um=32), + domain=DomainConfig( + dpml=1.0, + margin_xy=0.5, + margin_z_above=0.5, + margin_z_below=0.5, + port_margin=0.5, + extend_ports=0.0, + source_port_offset=0.1, + distance_source_to_monitors=0.2, + ), + accuracy=AccuracyConfig( + eps_averaging=False, + subpixel_maxeval=0, + subpixel_tol=1e-4, + simplify_tol=0.0, + ), + diagnostics=DiagnosticsConfig( + save_geometry=True, + save_fields=True, + save_epsilon_raw=False, + save_animation=False, + animation_interval=0.5, + preview_only=False, + verbose_interval=0, + ), + symmetries=[], + ) + + path = tmp_path / "config.json" + cfg.to_json(path) + data = json.loads(path.read_text()) + + assert data["source"]["source_type"] == "fiber" + assert data["source"]["beam_waist"] == 5.2 + assert data["source"]["angle_theta"] == 10.0 + assert data["source"]["z_position"] == 1.5 + + def test_json_roundtrip_mode_source(self, tmp_path): + """Existing mode source still works with the discriminated union.""" + from gsim.meep.models.config import AccuracyConfig, DiagnosticsConfig + + mode_cfg = SourceConfig() + + cfg = SimConfig( + gds_filename="layout.gds", + verbose_interval=0, + layer_stack=[], + dielectrics=[], + ports=[], + materials={}, + wavelength=WavelengthConfig(wavelength=1.55, bandwidth=0.1, num_freqs=11), + source=mode_cfg, + stopping=StoppingConfig( + mode="fixed", + max_time=100.0, + decay_dt=50.0, + decay_component="Ey", + threshold=0.05, + dft_min_run_time=100, + ), + resolution=ResolutionConfig(pixels_per_um=32), + domain=DomainConfig( + dpml=1.0, + margin_xy=0.5, + margin_z_above=0.5, + margin_z_below=0.5, + port_margin=0.5, + extend_ports=0.0, + source_port_offset=0.1, + distance_source_to_monitors=0.2, + ), + accuracy=AccuracyConfig( + eps_averaging=False, + subpixel_maxeval=0, + subpixel_tol=1e-4, + simplify_tol=0.0, + ), + diagnostics=DiagnosticsConfig( + save_geometry=True, + save_fields=True, + save_epsilon_raw=False, + save_animation=False, + animation_interval=0.5, + preview_only=False, + verbose_interval=0, + ), + symmetries=[], + ) + + path = tmp_path / "config.json" + cfg.to_json(path) + data = json.loads(path.read_text()) + + assert data["source"]["source_type"] == "mode" diff --git a/tests/meep/test_simulation.py b/tests/meep/test_simulation.py index b836a357..2d584825 100644 --- a/tests/meep/test_simulation.py +++ b/tests/meep/test_simulation.py @@ -8,6 +8,7 @@ from gsim.meep import ( FDTD, Domain, + FiberSource, Geometry, Material, ModeSource, @@ -215,6 +216,7 @@ class TestFieldAssignment: def test_source_port(self): sim = Simulation() + assert isinstance(sim.source, ModeSource) sim.source.port = "o1" assert sim.source.port == "o1" @@ -306,6 +308,7 @@ def test_source_callable(self): sim = Simulation() result = sim.source(port="o1", wavelength=1.31, wavelength_span=0.05) assert result is sim.source + assert isinstance(sim.source, ModeSource) assert sim.source.port == "o1" assert sim.source.wavelength == 1.31 assert sim.source.wavelength_span == 0.05 @@ -541,3 +544,140 @@ def test_simulation_instantiation(self): assert sim.materials == {} assert sim.monitors == [] assert sim.solver.resolution == 32 + + +# --------------------------------------------------------------------------- +# FiberSource model tests +# --------------------------------------------------------------------------- + + +class TestFiberSource: + """Test FiberSource model defaults and validation.""" + + def test_defaults(self): + s = FiberSource() + assert s.wavelength == 1.55 + assert s.wavelength_span == 0.1 + assert s.num_freqs == 11 + assert s.beam_waist == 5.2 + assert s.angle_theta == 10.0 + assert s.angle_phi == 0.0 + assert s.polarization == "TE" + assert s.position is None + assert s.z_offset == 1.0 + assert s.direction == "down" + + def test_custom(self): + s = FiberSource( + wavelength=1.31, + beam_waist=4.5, + angle_theta=15.0, + polarization="TM", + z_offset=2.0, + direction="up", + ) + assert s.wavelength == 1.31 + assert s.beam_waist == 4.5 + assert s.angle_theta == 15.0 + assert s.polarization == "TM" + assert s.z_offset == 2.0 + assert s.direction == "up" + + def test_callable_api(self): + s = FiberSource() + result = s(beam_waist=4.0, angle_theta=8.0) + assert result is s + assert s.beam_waist == 4.0 + assert s.angle_theta == 8.0 + + def test_beam_waist_positive(self): + with pytest.raises(ValidationError): + FiberSource(beam_waist=0) + + def test_angle_theta_bounds(self): + with pytest.raises(ValidationError): + FiberSource(angle_theta=90.0) + with pytest.raises(ValidationError): + FiberSource(angle_theta=-1.0) + + def test_z_offset_positive(self): + with pytest.raises(ValidationError): + FiberSource(z_offset=0) + + def test_position_custom(self): + s = FiberSource(position=[1.0, 2.0]) + assert s.position == [1.0, 2.0] + + +class TestFiberSourceSimulation: + """Test FiberSource integration with Simulation.""" + + def test_source_assignment(self): + sim = Simulation() + sim.source = FiberSource(beam_waist=5.2, angle_theta=10.0) + assert isinstance(sim.source, FiberSource) + assert sim.source.beam_waist == 5.2 + + def test_source_default_is_mode(self): + sim = Simulation() + assert isinstance(sim.source, ModeSource) + + def test_validate_fiber_skips_port_check(self): + """FiberSource validation should not require a source port.""" + sim = Simulation() + sim.source = FiberSource() + # Validation still fails because no component, but not due to source port + result = sim.validate_config() + assert not result.valid + assert any("No component" in e for e in result.errors) + assert not any("Source port" in e for e in result.errors) + + def test_wavelength_config_from_fiber(self): + sim = Simulation() + sim.source = FiberSource(wavelength=1.31, wavelength_span=0.05, num_freqs=21) + wl = sim._wavelength_config() + assert wl.wavelength == 1.31 + assert wl.bandwidth == 0.05 + assert wl.num_freqs == 21 + + def test_fiber_z_position_uses_core_layer(self): + """z_position should reference core (highest n) layer, not topmost layer.""" + from gsim.common.stack.extractor import Layer, LayerStack + + # Create a stack with core at z=0..0.22 and metal at z=1.8..2.5 + stack = LayerStack( + layers={ + "core": Layer( + name="core", + gds_layer=(1, 0), + zmin=0.0, + zmax=0.22, + thickness=0.22, + material="si", + layer_type="dielectric", + ), + "metal": Layer( + name="metal", + gds_layer=(12, 0), + zmin=1.8, + zmax=2.5, + thickness=0.7, + material="Aluminum", + layer_type="conductor", + ), + } + ) + sim = Simulation() + sim.source = FiberSource(z_offset=2.0, direction="down", position=[0.0, 0.0]) + cfg = sim._fiber_source_config(stack) + # z_position should be core_top + z_offset = 0.22 + 2.0 = 2.22 + # NOT metal_top + z_offset = 2.5 + 2.0 = 4.5 + assert cfg.z_position == pytest.approx(2.22, abs=0.01) + + def test_fiber_validation_warns_z_offset(self): + """Warn when z_offset >= margin_z_above.""" + sim = Simulation() + sim.source = FiberSource(z_offset=5.0) + sim.domain(margin_z_above=3.0) + result = sim.validate_config() + assert any("z_offset" in w for w in result.warnings)