diff --git a/.gitignore b/.gitignore index 0de6b17b..afd2c414 100644 --- a/.gitignore +++ b/.gitignore @@ -1,6 +1,7 @@ # Simulation outputs sim-data-*/ palace-sim-*/ +meep-sim-*/ # files *.c diff --git a/mkdocs.yml b/mkdocs.yml index 9a7b0d92..b3289192 100644 --- a/mkdocs.yml +++ b/mkdocs.yml @@ -20,6 +20,7 @@ nav: - home: index.md - examples: - Simple Example: nbs/example.md + - MEEP Y-Branch: nbs/meep_ybranch.md - Palace CPW: nbs/palace_demo_cpw.md - Palace Microstrip: nbs/palace_demo_microstrip.md - api reference: diff --git a/nbs/meep_ybranch.ipynb b/nbs/meep_ybranch.ipynb new file mode 100644 index 00000000..480f15d0 --- /dev/null +++ b/nbs/meep_ybranch.ipynb @@ -0,0 +1,155 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "0", + "metadata": {}, + "source": [ + "# Running MEEP Simulations\n", + "\n", + "[MEEP](https://meep.readthedocs.io/) is an open-source FDTD electromagnetic simulator. This notebook demonstrates using the `gsim.meep` API to run an S-parameter simulation on a photonic Y-branch.\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 pcell 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_y_1550()\n", + "c" + ] + }, + { + "cell_type": "markdown", + "id": "3", + "metadata": {}, + "source": [ + "### Configure and run simulation" + ] + }, + { + "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\n", + "sim.geometry.z_crop = \"auto\"\n", + "\n", + "sim.materials = {\"si\": 3.47, \"SiO2\": 1.44}\n", + "\n", + "sim.source.port = \"o1\"\n", + "sim.source.wavelength = 1.55\n", + "sim.source.bandwidth = 0.01\n", + "sim.source.num_freqs = 11\n", + "\n", + "sim.monitors = [\"o1\", \"o2\", \"o3\"]\n", + "\n", + "sim.domain.pml = 1.0\n", + "sim.domain.margin = 0.5\n", + "sim.domain.margin_z_above = 0.5\n", + "sim.domain.margin_z_below = 0.5\n", + "sim.domain.port_margin = 0.5\n", + "\n", + "sim.solver.resolution = 20\n", + "sim.solver.stopping = \"dft_decay\"\n", + "sim.solver.max_time = 200\n", + "sim.solver.stopping_threshold = 1e-3\n", + "sim.solver.stopping_min_time = 100\n", + "sim.solver.subpixel = False\n", + "sim.solver.simplify_tol = 0.01\n", + "sim.solver.save_geometry = True\n", + "sim.solver.save_fields = True\n", + "sim.solver.save_animation = True\n", + "sim.solver.verbose_interval = 5.0\n", + "\n", + "print(sim.validate_config())" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5", + "metadata": {}, + "outputs": [], + "source": [ + "sim.plot_2d(slices=\"xyz\")" + ] + }, + { + "cell_type": "markdown", + "id": "6", + "metadata": {}, + "source": [ + "### Run simulation on cloud" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7", + "metadata": {}, + "outputs": [], + "source": [ + "# Run on GDSFactory+ cloud\n", + "result = sim.run()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8", + "metadata": {}, + "outputs": [], + "source": [ + "result.plot(db=True)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9", + "metadata": {}, + "outputs": [], + "source": [ + "result.show_animation()" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "name": "python" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/pyproject.toml b/pyproject.toml index 23a1aa3a..d8359e0a 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -57,7 +57,8 @@ dev = [ "ty>=0.0.13" ] docs = [ - "ihp-gdsfactory>=0.1.4" + "ihp-gdsfactory>=0.1.4", + "ubcpdk @ git+https://github.com/gdsfactory/ubc.git" ] [[tool.bver.file]] @@ -72,6 +73,9 @@ src = "src/gsim/__init__.py" [tool.bver.git] commit-template = "Release {new-version}" +[tool.codespell] +ignore-words-list = "doubleclick" + [tool.mypy] mypy_path = ["src"] python_version = "3.12" @@ -92,6 +96,16 @@ disallow_untyped_defs = false module = "gsim.palace.*" warn_return_any = false +[[tool.mypy.overrides]] +check_untyped_defs = false +disable_error_code = ["no-untyped-def", "no-untyped-call", "var-annotated", "assignment"] +disallow_any_generics = false +disallow_untyped_defs = false +# The MEEP subpackage wraps cloud-only MEEP simulation. Relaxed mypy +# to match the Palace approach during early development. +module = "gsim.meep.*" +warn_return_any = false + [[tool.mypy.overrides]] ignore_errors = true module = "tests.*" @@ -198,6 +212,11 @@ select = ["ALL"] "INP001", # implicit-namespace-package "PERF401" # use-list-extend (keeping for readability) ] +"nbs/*.py" = [ + "EXE001", # shebang-present-but-not-executable + "INP001", # implicit-namespace-package + "T201" # print (natural in scripts) +] "scripts/*.py" = [ "ANN", # flake8-annotations "ARG001", # unused-function-argument @@ -220,6 +239,7 @@ select = ["ALL"] "PLC0415", # allow imports inside tests "PT011", # pytest-raises-too-broad "S101", # assert + "SLF001", # private-member-access (testing internal state) "T201" # print ] @@ -228,3 +248,23 @@ convention = "google" [tool.setuptools.packages.find] where = ["src"] + +[tool.ty.analysis] +# Optional runtime dependencies that are not always installed locally +allowed-unresolved-imports = ["open3d", "open3d.**", "meep", "meep.**"] + +[[tool.ty.overrides]] +include = ["tests/**"] + +[tool.ty.overrides.rules] +possibly-missing-attribute = "ignore" + +[tool.ty.rules] +missing-argument = "ignore" +possibly-missing-attribute = "ignore" +# Suppress rules that generate many false-positives in this codebase +# (optional deps, third-party stubs, matplotlib Axes|None, etc.) +unused-ignore-comment = "ignore" + +[tool.ty.src] +include = ["src"] diff --git a/src/gsim/__init__.py b/src/gsim/__init__.py index 5d34b757..02862c91 100644 --- a/src/gsim/__init__.py +++ b/src/gsim/__init__.py @@ -5,6 +5,7 @@ Currently includes: - palace: Palace EM simulation API + - meep: MEEP photonic FDTD simulation API """ from __future__ import annotations diff --git a/src/gsim/common/__init__.py b/src/gsim/common/__init__.py index affa0707..e1f52134 100644 --- a/src/gsim/common/__init__.py +++ b/src/gsim/common/__init__.py @@ -11,6 +11,7 @@ from __future__ import annotations from gsim.common.geometry import Geometry +from gsim.common.geometry_model import GeometryModel, Prism, extract_geometry_model from gsim.common.stack import ( MATERIALS_DB, Layer, @@ -37,13 +38,16 @@ __all__ = [ "MATERIALS_DB", "Geometry", + "GeometryModel", "Layer", "LayerStack", "MaterialProperties", + "Prism", "Stack", "StackLayer", "ValidationResult", "extract_from_pdk", + "extract_geometry_model", "extract_layer_stack", "get_material_properties", "get_stack", diff --git a/src/gsim/common/geometry_model.py b/src/gsim/common/geometry_model.py new file mode 100644 index 00000000..f15d8f19 --- /dev/null +++ b/src/gsim/common/geometry_model.py @@ -0,0 +1,419 @@ +"""Solver-agnostic geometry model for 3D visualization. + +Provides generic Prism and GeometryModel dataclasses that represent +extruded 2D polygons without any solver dependency (no meep, no tidy3d). +These are suitable for client-side visualization and geometry inspection. + +Classes: + Prism: A 2D polygon extruded along z. + GeometryModel: Collection of prisms organized by layer. + +Functions: + extract_geometry_model: Convert a LayeredComponentBase into a GeometryModel. +""" + +from __future__ import annotations + +import logging +import warnings +from dataclasses import dataclass, field +from typing import TYPE_CHECKING, Any + +import numpy as np +from scipy.spatial import Delaunay +from shapely.geometry import Point + +if TYPE_CHECKING: + from shapely import Polygon + + from gsim.common.layered_component import LayeredComponentBase + +logger = logging.getLogger(__name__) + + +# --------------------------------------------------------------------------- +# Prism +# --------------------------------------------------------------------------- + + +@dataclass(frozen=True) +class Prism: + """A 2D polygon extruded in z. No solver dependency. + + Attributes: + vertices: (N, 2) numpy array of xy coordinates defining the polygon. + z_base: Bottom z coordinate of the extrusion. + z_top: Top z coordinate of the extrusion. + layer_name: Name of the layer this prism belongs to. + material: Name of the material assigned to this prism. + sidewall_angle: Sidewall taper angle in degrees (gdsfactory convention). + original_polygon: Optional reference to the source Shapely polygon. + """ + + vertices: np.ndarray # (N, 2) xy coords + z_base: float + z_top: float + layer_name: str = "" + material: str = "" + sidewall_angle: float = 0.0 + original_polygon: Any = field(default=None, repr=False, compare=False) + + @property + def height(self) -> float: + """Extrusion height (z_top - z_base).""" + return self.z_top - self.z_base + + @property + def z_center(self) -> float: + """Center z coordinate.""" + return (self.z_base + self.z_top) / 2.0 + + +# --------------------------------------------------------------------------- +# GeometryModel +# --------------------------------------------------------------------------- + + +@dataclass +class GeometryModel: + """Complete 3D geometry: layers of prisms, ready for visualization. + + Attributes: + prisms: Mapping from layer name to list of Prism objects. + bbox: Axis-aligned 3D bounding box as ((xmin, ymin, zmin), (xmax, ymax, zmax)). + layer_bboxes: Optional per-layer bounding boxes for 2D slice logic. + layer_mesh_orders: Optional mapping of layer_name -> mesh_order for + z-ordering in 2D plots. + """ + + prisms: dict[str, list[Prism]] + bbox: tuple[tuple[float, float, float], tuple[float, float, float]] + layer_bboxes: dict[ + str, + tuple[tuple[float, float, float], tuple[float, float, float]], + ] = field(default_factory=dict) + layer_mesh_orders: dict[str, int] = field(default_factory=dict) + + @property + def all_prisms(self) -> list[Prism]: + """Flat list of all prisms across all layers.""" + return [p for layer_prisms in self.prisms.values() for p in layer_prisms] + + @property + def layer_names(self) -> list[str]: + """Ordered list of layer names that contain prisms.""" + return list(self.prisms.keys()) + + @property + def size(self) -> tuple[float, float, float]: + """(dx, dy, dz) extent of the bounding box.""" + mn, mx = self.bbox + return (mx[0] - mn[0], mx[1] - mn[1], mx[2] - mn[2]) + + def get_layer_center(self, layer_name: str) -> tuple[float, float, float]: + """Return the center of a layer's bounding box. + + Falls back to the geometry-wide center if no per-layer bbox is stored. + """ + bb = self.get_layer_bbox(layer_name) + mn, mx = bb + return ( + (mn[0] + mx[0]) / 2.0, + (mn[1] + mx[1]) / 2.0, + (mn[2] + mx[2]) / 2.0, + ) + + def get_layer_bbox( + self, + layer_name: str, + ) -> tuple[tuple[float, float, float], tuple[float, float, float]]: + """Return the bounding box for a specific layer. + + Falls back to the geometry-wide bbox if no per-layer bbox is stored. + """ + return self.layer_bboxes.get(layer_name, self.bbox) + + +# --------------------------------------------------------------------------- +# Triangulation helpers (ported from fdtd/geometry/core.py, no meep) +# --------------------------------------------------------------------------- + + +def _triangulate_polygon_with_holes( + polygon: Polygon, + z_base: float, + z_top: float, + sidewall_angle: float, + layer_name: str, + material: str, +) -> list[Prism]: + """Create triangular Prisms from a polygon with holes via Delaunay. + + Collects all boundary points (exterior + interiors), runs Delaunay + triangulation, and keeps only triangles whose centroid lies inside + the original polygon. + """ + all_points: list[tuple[float, float]] = [] + all_points.extend(list(polygon.exterior.coords[:-1])) + for interior in polygon.interiors: + all_points.extend(list(interior.coords[:-1])) + + if len(all_points) < 3: + # Degenerate case -- fall back to exterior-only prism + return _prism_from_exterior( + polygon, + z_base, + z_top, + sidewall_angle, + layer_name, + material, + ) + + # Optimised ring triangulation for single-hole, similarly-sized rings + if ( + len(polygon.interiors) == 1 + and len(all_points) < 200 + and abs( + len(list(polygon.exterior.coords)) - len(list(polygon.interiors[0].coords)), + ) + < 10 + ): + return _ring_triangulation( + polygon, + z_base, + z_top, + sidewall_angle, + layer_name, + material, + ) + + points_2d = np.array(all_points) + tri = Delaunay(points_2d) + + triangular_prisms: list[Prism] = [] + for simplex in tri.simplices: + triangle_pts = points_2d[simplex] + centroid = np.mean(triangle_pts, axis=0) + if polygon.contains(Point(centroid[0], centroid[1])): + triangular_prisms.append( + Prism( + vertices=triangle_pts.copy(), + z_base=z_base, + z_top=z_top, + layer_name=layer_name, + material=material, + sidewall_angle=sidewall_angle, + original_polygon=polygon, + ), + ) + + if not triangular_prisms: + warnings.warn( + "No valid triangles found for holed polygon; " + "falling back to exterior-only prism.", + stacklevel=2, + ) + return _prism_from_exterior( + polygon, + z_base, + z_top, + sidewall_angle, + layer_name, + material, + ) + + logger.debug( + "Created %d triangular prisms for polygon with %d holes", + len(triangular_prisms), + len(polygon.interiors), + ) + return triangular_prisms + + +def _ring_triangulation( + polygon: Polygon, + z_base: float, + z_top: float, + sidewall_angle: float, + layer_name: str, + material: str, +) -> list[Prism]: + """Efficient fan triangulation for a polygon with a single hole. + + Pairs exterior vertices with interior vertices proportionally and + creates two triangles per exterior edge, covering the ring area. + """ + exterior_coords = list(polygon.exterior.coords[:-1]) + interior_coords = list(polygon.interiors[0].coords[:-1]) + + n_outer = len(exterior_coords) + n_inner = len(interior_coords) + + triangular_prisms: list[Prism] = [] + + for i in range(n_outer): + next_i = (i + 1) % n_outer + inner_i = int(i * n_inner / n_outer) % n_inner + inner_next = int(next_i * n_inner / n_outer) % n_inner + + tri1 = np.array( + [exterior_coords[i], exterior_coords[next_i], interior_coords[inner_i]], + ) + triangular_prisms.append( + Prism( + vertices=tri1, + z_base=z_base, + z_top=z_top, + layer_name=layer_name, + material=material, + sidewall_angle=sidewall_angle, + original_polygon=polygon, + ), + ) + + tri2 = np.array( + [ + exterior_coords[next_i], + interior_coords[inner_next], + interior_coords[inner_i], + ], + ) + triangular_prisms.append( + Prism( + vertices=tri2, + z_base=z_base, + z_top=z_top, + layer_name=layer_name, + material=material, + sidewall_angle=sidewall_angle, + original_polygon=polygon, + ), + ) + + logger.debug( + "Ring triangulation: %d prisms (from %d boundary points)", + len(triangular_prisms), + n_outer + n_inner, + ) + return triangular_prisms + + +def _prism_from_exterior( + polygon: Polygon, + z_base: float, + z_top: float, + sidewall_angle: float, + layer_name: str, + material: str, +) -> list[Prism]: + """Create a single Prism from a polygon's exterior ring (ignoring holes).""" + coords = np.array(polygon.exterior.coords[:-1]) + return [ + Prism( + vertices=coords, + z_base=z_base, + z_top=z_top, + layer_name=layer_name, + material=material, + sidewall_angle=sidewall_angle, + original_polygon=polygon, + ), + ] + + +# --------------------------------------------------------------------------- +# Main extraction function +# --------------------------------------------------------------------------- + + +def extract_geometry_model( + layered_component: LayeredComponentBase, +) -> GeometryModel: + """Convert a LayeredComponentBase into a GeometryModel with generic Prisms. + + For each geometry layer (sorted by mesh_order): + 1. Retrieve the merged Shapely polygon from *layered_component.polygons*. + 2. Compute z_base / z_top from *get_layer_bbox*. + 3. Iterate sub-polygons for MultiPolygon geometries. + 4. Handle polygons with holes via Delaunay triangulation. + 5. Produce Prism objects with (N, 2) numpy vertex arrays. + + Args: + layered_component: A LayeredComponentBase (or subclass) instance + that provides polygons, geometry_layers, and get_layer_bbox. + + Returns: + A GeometryModel containing all extracted prisms and the overall + 3D bounding box. + """ + all_prisms: dict[str, list[Prism]] = {} + layer_bboxes: dict[ + str, + tuple[tuple[float, float, float], tuple[float, float, float]], + ] = {} + layer_mesh_orders: dict[str, int] = {} + + # Sort layers by mesh_order (ascending) for consistent rendering order + sorted_layers = sorted( + layered_component.geometry_layers.items(), + key=lambda item: item[1].mesh_order, + ) + + for name, layer in sorted_layers: + shape = layered_component.polygons[name] + bbox = layered_component.get_layer_bbox(name) + layer_bboxes[name] = bbox + layer_mesh_orders[name] = ( + layer.mesh_order if layer.mesh_order is not None else 0 + ) + z_base = bbox[0][2] + z_top = bbox[1][2] + + sidewall_angle = layer.sidewall_angle or 0.0 + material_name = str(layer.material) if layer.material else "" + + # Normalise to a list of individual Polygon objects + if hasattr(shape, "geoms"): + polygons: list[Polygon] = list(shape.geoms) + else: + polygons = [shape] + + layer_prisms: list[Prism] = [] + + for polygon in polygons: + if polygon.is_empty or not polygon.is_valid: + continue + + if hasattr(polygon, "interiors") and polygon.interiors: + layer_prisms.extend( + _triangulate_polygon_with_holes( + polygon, + z_base=z_base, + z_top=z_top, + sidewall_angle=sidewall_angle, + layer_name=name, + material=material_name, + ), + ) + else: + coords = np.array(polygon.exterior.coords[:-1]) + layer_prisms.append( + Prism( + vertices=coords, + z_base=z_base, + z_top=z_top, + layer_name=name, + material=material_name, + sidewall_angle=sidewall_angle, + original_polygon=polygon, + ), + ) + + all_prisms[name] = layer_prisms + + return GeometryModel( + prisms=all_prisms, + bbox=layered_component.bbox, + layer_bboxes=layer_bboxes, + layer_mesh_orders=layer_mesh_orders, + ) diff --git a/src/gsim/common/layered_component.py b/src/gsim/common/layered_component.py new file mode 100644 index 00000000..479e2c14 --- /dev/null +++ b/src/gsim/common/layered_component.py @@ -0,0 +1,300 @@ +"""Layered component base class for 3D geometry from GDS + PDK. + +Ported from gplugins.common.base_models.component.LayeredComponentBase +to eliminate the gplugins dependency. Combines a gdsfactory Component +with a LayerStack and provides: + - polygon extraction with DerivedLayer support (via `fuse_polygons`) + - z-coordinate computation per layer + - bounding box, center, size properties + - port extension and padding +""" + +from __future__ import annotations + +from functools import cached_property +from hashlib import md5 + +import gdsfactory as gf +import numpy as np +from gdsfactory.component import Component +from gdsfactory.pdk import get_layer_name +from gdsfactory.technology import LayerLevel, LayerStack +from pydantic import BaseModel, ConfigDict, NonNegativeFloat, computed_field + +from gsim.common.polygon import cleanup_component +from gsim.common.types import AnyShapelyPolygon, GFComponent + + +class LayeredComponentBase(BaseModel): + """Base class combining a gdsfactory Component with a LayerStack. + + Provides polygon extraction, z-coordinate management, bounding box, + port handling, and layer ordering — everything needed for 3D geometry + construction without any solver dependency. + """ + + model_config = ConfigDict( + frozen=True, + extra="forbid", + arbitrary_types_allowed=True, + allow_inf_nan=False, + validate_return=True, + ) + + component: GFComponent + layer_stack: LayerStack + extend_ports: NonNegativeFloat = 0.0 + port_offset: float = 0.0 + pad_xy_inner: float = 0.0 + pad_xy_outer: NonNegativeFloat = 0.0 + pad_z_inner: float = 0.0 + pad_z_outer: NonNegativeFloat = 0.0 + wafer_layer: tuple[int, int] = (999, 0) + slice_stack: tuple[int, int | None] = (0, None) + + def __hash__(self) -> int: + if not hasattr(self, "_hash"): + dump = str.encode(self.model_dump_json()) + object.__setattr__(self, "_hash", int(md5(dump).hexdigest()[:15], 16)) + return self._hash # type: ignore[return-value] + + # ------------------------------------------------------------------ + # Padding helpers + # ------------------------------------------------------------------ + + @property + def pad_xy(self) -> float: + """Total XY padding.""" + return self.pad_xy_inner + self.pad_xy_outer + + @property + def pad_z(self) -> float: + """Total Z padding.""" + return self.pad_z_inner + self.pad_z_outer + + # ------------------------------------------------------------------ + # GDS component with extended ports + wafer layer + # ------------------------------------------------------------------ + + @cached_property + def gds_component(self) -> GFComponent: + """Build GDS component with extended ports and wafer layer.""" + c = Component() + c << gf.components.extend_ports( + self.component, length=self.extend_ports + self.pad_xy + ) + (xmin, ymin), (xmax, ymax) = self._gds_bbox + delta = self.pad_xy_outer + points = np.array( + [ + [xmin - delta, ymin - delta], + [xmax + delta, ymin - delta], + [xmax + delta, ymax + delta], + [xmin - delta, ymax + delta], + ] + ) + c.add_polygon(points, layer=self.wafer_layer) + c.add_ports(self.ports) + c.copy_child_info(self.component) + return c + + @cached_property + def _gds_bbox(self) -> tuple[tuple[float, float], tuple[float, float]]: + c = gf.components.extend_ports( + self.component, length=self.extend_ports + self.pad_xy_inner + ) + unchanged = np.isclose( + np.abs(np.round(c.bbox_np() - self.component.bbox_np(), 3)), 0 + ) + bbox = ( + c.bbox_np() + unchanged * np.array([[-1, -1], [1, 1]]) * self.pad_xy_inner + ) + return tuple(map(tuple, bbox)) # type: ignore[return-value] + + # ------------------------------------------------------------------ + # Ports + # ------------------------------------------------------------------ + + @cached_property + def ports(self) -> tuple[gf.Port, ...]: + """Component ports offset by extension and padding.""" + p = tuple( + p.copy_polar( + self.extend_ports + self.pad_xy_inner - self.port_offset, + orientation=p.orientation, + ) + for p in self.component.ports + ) + for pi, po in zip(self.component.ports, p, strict=True): + po.angle = pi.angle + return p + + @cached_property + def port_names(self) -> tuple[str, ...]: + """Names of all ports.""" + return tuple(p.name for p in self.ports if p.name is not None) + + @cached_property + def port_centers(self) -> tuple[tuple[float, float, float], ...]: + """3D center coordinates of all ports.""" + return tuple(self.get_port_center(p) for p in self.ports) + + def get_port_center(self, port: gf.Port) -> tuple[float, float, float]: + """Return 3D center of a port averaged over its layers.""" + layers = self.get_port_layers(port) + return ( + *port.dcenter, + np.mean([self.get_layer_center(layer)[2] for layer in layers]), + ) + + def get_port_layers(self, port: gf.Port) -> list[str]: + """Return layer names associated with a port.""" + layer_name = get_layer_name(port.layer) + if "_intent" in layer_name: + layer_name = layer_name.replace("_intent", "") + + derived_layers = [] + for l_name, level in self.layer_stack.layers.items(): + if layer_name in str(level.layer): + derived_layers.append(l_name) + return derived_layers + + # ------------------------------------------------------------------ + # Polygon extraction (uses gsim.common.polygon) + # ------------------------------------------------------------------ + + @computed_field + @cached_property + def polygons(self) -> dict[str, AnyShapelyPolygon]: + """Cleaned-up Shapely polygons keyed by layer name.""" + return cleanup_component( + self.gds_component, + self.layer_stack, + round_tol=3, + simplify_tol=1e-3, + ) + + # ------------------------------------------------------------------ + # Layer helpers + # ------------------------------------------------------------------ + + @cached_property + def geometry_layers(self) -> dict[str, LayerLevel]: + """Non-empty layers with valid z-coordinates, sorted by height.""" + layers = { + k: v + for k, v in self.layer_stack.layers.items() + if not self.polygons[k].is_empty + and v.zmin is not None + and v.thickness is not None + } + layers = dict(sorted(layers.items(), key=lambda x: x[1].zmin + x[1].thickness)) + return dict(tuple(layers.items())[slice(*self.slice_stack)]) + + @cached_property + def bottom_layer(self) -> str: + """Name of the lowest layer in the stack.""" + return min( + self.geometry_layers.items(), + key=lambda item: min(item[1].zmin, item[1].zmin + item[1].thickness), + )[0] + + @cached_property + def top_layer(self) -> str: + """Name of the highest layer in the stack.""" + return max( + self.geometry_layers.items(), + key=lambda item: max(item[1].zmin, item[1].zmin + item[1].thickness), + )[0] + + @cached_property + def device_layers(self) -> tuple[str, ...]: + """Layer names present in the component's GDS layers.""" + return tuple( + k + for k, v in self.layer_stack.layers.items() + if v.layer in self.component.layers + ) + + # ------------------------------------------------------------------ + # Bounding box / geometry + # ------------------------------------------------------------------ + + @property + def xmin(self) -> float: + """Minimum x coordinate.""" + return self._gds_bbox[0][0] + + @property + def xmax(self) -> float: + """Maximum x coordinate.""" + return self._gds_bbox[1][0] + + @property + def ymin(self) -> float: + """Minimum y coordinate.""" + return self._gds_bbox[0][1] + + @property + def ymax(self) -> float: + """Maximum y coordinate.""" + return self._gds_bbox[1][1] + + @cached_property + def zmin(self) -> float: + """Minimum z coordinate including inner padding.""" + return ( + min( + min(layer.zmin, layer.zmin + layer.thickness) + for layer in self.geometry_layers.values() + ) + - self.pad_z_inner + ) + + @cached_property + def zmax(self) -> float: + """Maximum z coordinate including inner padding.""" + return ( + max( + max(layer.zmin, layer.zmin + layer.thickness) + for layer in self.geometry_layers.values() + ) + + self.pad_z_inner + ) + + @property + def bbox( + self, + ) -> tuple[tuple[float, float, float], tuple[float, float, float]]: + """3D bounding box as ((xmin, ymin, zmin), (xmax, ymax, zmax)).""" + return (*self._gds_bbox[0], self.zmin), (*self._gds_bbox[1], self.zmax) + + @property + def center(self) -> tuple[float, float, float]: + """3D center of the bounding box.""" + return tuple(np.mean(self.bbox, axis=0)) # type: ignore[return-value] + + @property + def size(self) -> tuple[float, float, float]: + """3D size of the bounding box (dx, dy, dz).""" + return tuple(np.squeeze(np.diff(self.bbox, axis=0))) # type: ignore[return-value] + + def get_layer_bbox( + self, layername: str + ) -> tuple[tuple[float, float, float], tuple[float, float, float]]: + """Return 3D bounding box of a single layer.""" + layer = self.layer_stack[layername] + bounds_xy = self.polygons[layername].bounds + zmin, zmax = sorted([layer.zmin, layer.zmin + layer.thickness]) + + if layername == self.bottom_layer: + zmin -= self.pad_z + if layername == self.top_layer: + zmax += self.pad_z + + return (*bounds_xy[:2], zmin), (*bounds_xy[2:], zmax) + + def get_layer_center(self, layername: str) -> tuple[float, float, float]: + """Return 3D center of a single layer.""" + bbox = self.get_layer_bbox(layername) + return tuple(np.mean(bbox, axis=0)) # type: ignore[return-value] diff --git a/src/gsim/common/polygon.py b/src/gsim/common/polygon.py new file mode 100644 index 00000000..88f6a162 --- /dev/null +++ b/src/gsim/common/polygon.py @@ -0,0 +1,114 @@ +"""Polygon extraction and processing for layered GDS components. + +Ported from gplugins.common.base_models.component to avoid the gplugins +dependency. Uses gdsfactory's DerivedLayer/LogicalLayer `.get_shapes()` +to resolve boolean layer operations (e.g., WG - DEEP_ETCH) and returns +merged Shapely polygons. +""" + +from __future__ import annotations + +from typing import TYPE_CHECKING + +import shapely +import shapely.geometry +import shapely.ops + +if TYPE_CHECKING: + from gdsfactory.technology import LayerStack + + from gsim.common.types import AnyShapelyPolygon, GFComponent + + +def round_coordinates(geom: AnyShapelyPolygon, ndigits: int = 4) -> AnyShapelyPolygon: + """Round polygon coordinates to eliminate floating point noise.""" + + def _round_coords(x: float, y: float, z: float | None = None) -> list[float]: + x = round(x, ndigits) + y = round(y, ndigits) + if z is not None: + z = round(z, ndigits) + return [c for c in (x, y, z) if c is not None] + + return shapely.ops.transform(_round_coords, geom) + + +def fuse_polygons( + component: GFComponent, + layer, + round_tol: int = 4, + simplify_tol: float = 1e-4, +) -> AnyShapelyPolygon: + """Extract and merge all polygons for a layer from a component. + + Calls ``layer.get_shapes(component)`` which handles DerivedLayer boolean + operations (KLayout Region). Converts the resulting KLayout polygons to + Shapely, including holes, then merges and simplifies. + + Args: + component: gdsfactory Component + layer: gdsfactory LayerLevel (has ``.get_shapes()`` via its + ``.layer`` attribute which is a LogicalLayer or DerivedLayer) + round_tol: decimal places for coordinate rounding + simplify_tol: tolerance for polygon simplification + + Returns: + Merged Shapely Polygon or MultiPolygon + """ + layer_region = layer.get_shapes(component) + + shapely_polygons = [] + for klayout_polygon in layer_region.each_merged(): + exterior_points = [ + (point.x / 1000, point.y / 1000) + for point in klayout_polygon.each_point_hull() + ] + interior_points = [] + for hole_index in range(klayout_polygon.holes()): + hole_points = [ + (point.x / 1000, point.y / 1000) + for point in klayout_polygon.each_point_hole(hole_index) + ] + interior_points.append(hole_points) + + shapely_polygons.append( + round_coordinates( + shapely.geometry.Polygon(shell=exterior_points, holes=interior_points), + round_tol, + ) + ) + + return shapely.ops.unary_union(shapely_polygons).simplify( + simplify_tol, preserve_topology=False + ) + + +def cleanup_component( + component: GFComponent, + layer_stack: LayerStack, + round_tol: int = 2, + simplify_tol: float = 1e-2, +) -> dict[str, AnyShapelyPolygon]: + """Extract and fuse polygons for every layer in a stack. + + Args: + component: gdsfactory Component + layer_stack: gdsfactory LayerStack + round_tol: decimal places for coordinate rounding + simplify_tol: polygon simplification tolerance + + Returns: + Dict mapping layer name to merged Shapely polygon + """ + layer_stack_dict = layer_stack.to_dict() + + return { + layername: fuse_polygons( + component, + layer["layer"], + round_tol=round_tol, + simplify_tol=simplify_tol, + ) + for layername, layer in layer_stack_dict.items() + if layer["layer"] is not None + } diff --git a/src/gsim/common/stack/extractor.py b/src/gsim/common/stack/extractor.py index 638763ca..fd8a6090 100644 --- a/src/gsim/common/stack/extractor.py +++ b/src/gsim/common/stack/extractor.py @@ -35,6 +35,7 @@ class Layer(BaseModel): thickness: float # um material: str layer_type: Literal["conductor", "via", "dielectric", "substrate"] + sidewall_angle: float = 0.0 # degrees mesh_resolution: str | float = "medium" def get_mesh_size(self, base_size: float = 1.0) -> float: @@ -58,7 +59,7 @@ def get_mesh_size(self, base_size: float = 1.0) -> float: def to_dict(self) -> dict: """Convert to dictionary for YAML output.""" - return { + d = { "gds_layer": list(self.gds_layer), "zmin": round(self.zmin, 4), "zmax": round(self.zmax, 4), @@ -67,6 +68,9 @@ def to_dict(self) -> dict: "type": self.layer_type, "mesh_resolution": self.mesh_resolution, } + if self.sidewall_angle != 0.0: + d["sidewall_angle"] = self.sidewall_angle + return d class ValidationResult(BaseModel): @@ -387,9 +391,10 @@ def extract_layer_stack( zmin = layer_level.zmin if layer_level.zmin is not None else 0.0 thickness = layer_level.thickness if layer_level.thickness is not None else 0.0 zmax = zmin + thickness - material = layer_level.material if layer_level.material else "unknown" + material = layer_level.material or "unknown" gds_layer = _get_gds_layer_tuple(layer_level) layer_type = _classify_layer_type(layer_name, material) + sidewall_angle = getattr(layer_level, "sidewall_angle", 0.0) or 0.0 if layer_type == "substrate" and not include_substrate: continue @@ -402,6 +407,7 @@ def extract_layer_stack( thickness=thickness, material=material, layer_type=layer_type, + sidewall_angle=float(sidewall_angle), ) stack.layers[layer_name] = layer diff --git a/src/gsim/common/stack/materials.py b/src/gsim/common/stack/materials.py index 7d9f5c95..e27ee6c5 100644 --- a/src/gsim/common/stack/materials.py +++ b/src/gsim/common/stack/materials.py @@ -21,6 +21,10 @@ class MaterialProperties(BaseModel): permittivity: float | None = Field(default=None, ge=1.0) # relative permittivity loss_tangent: float | None = Field(default=None, ge=0, le=1) + # Optical properties (for photonic simulation, e.g. MEEP) + refractive_index: float | None = Field(default=None, gt=0) + extinction_coeff: float | None = Field(default=None, ge=0) + def to_dict(self) -> dict[str, object]: """Convert to dictionary for YAML output.""" d: dict[str, object] = {"type": self.type} @@ -30,6 +34,10 @@ def to_dict(self) -> dict[str, object]: d["permittivity"] = self.permittivity if self.loss_tangent is not None: d["loss_tangent"] = self.loss_tangent + if self.refractive_index is not None: + d["refractive_index"] = self.refractive_index + if self.extinction_coeff is not None: + d["extinction_coeff"] = self.extinction_coeff return d @classmethod @@ -46,6 +54,24 @@ def dielectric( type="dielectric", permittivity=permittivity, loss_tangent=loss_tangent ) + @classmethod + def optical( + cls, + refractive_index: float, + extinction_coeff: float = 0.0, + ) -> MaterialProperties: + """Create a material with optical properties for photonic simulation. + + Args: + refractive_index: Refractive index (n) + extinction_coeff: Extinction coefficient (k), default 0 + """ + return cls( + type="dielectric", + refractive_index=refractive_index, + extinction_coeff=extinction_coeff, + ) + # Material properties database # Sources: @@ -82,6 +108,7 @@ def dielectric( type="dielectric", permittivity=4.1, # Matches gds2palace IHP SG13G2 loss_tangent=0.0, + refractive_index=1.44, ), "passive": MaterialProperties( type="dielectric", @@ -92,6 +119,7 @@ def dielectric( type="dielectric", permittivity=7.5, loss_tangent=0.001, + refractive_index=2.0, ), "polyimide": MaterialProperties( type="dielectric", @@ -102,22 +130,26 @@ def dielectric( type="dielectric", permittivity=1.0, loss_tangent=0.0, + refractive_index=1.0, ), "vacuum": MaterialProperties( type="dielectric", permittivity=1.0, loss_tangent=0.0, + refractive_index=1.0, ), # Semiconductors (conductivity values from gds2palace IHP SG13G2) "silicon": MaterialProperties( type="semiconductor", permittivity=11.9, conductivity=2.0, # ~50 Ω·cm substrate (matches gds2palace) + refractive_index=3.47, ), "si": MaterialProperties( type="semiconductor", permittivity=11.9, conductivity=2.0, + refractive_index=3.47, ), } diff --git a/src/gsim/common/stack/visualization.py b/src/gsim/common/stack/visualization.py index 4d323ab5..18a99731 100644 --- a/src/gsim/common/stack/visualization.py +++ b/src/gsim/common/stack/visualization.py @@ -88,7 +88,7 @@ def parse_layer_stack(layer_stack: GfLayerStack) -> list[StackLayer]: zmin = level.zmin if level.zmin is not None else 0.0 thickness = level.thickness if level.thickness is not None else 0.0 zmax = zmin + thickness - material = level.material if level.material else None + material = level.material or None gds_layer = _get_gds_layer_number(level) layer_type = _classify_layer(name) diff --git a/src/gsim/common/types.py b/src/gsim/common/types.py new file mode 100644 index 00000000..312a8ed3 --- /dev/null +++ b/src/gsim/common/types.py @@ -0,0 +1,10 @@ +"""Type aliases shared across gsim modules.""" + +from __future__ import annotations + +from typing import Any + +from shapely import MultiPolygon, Polygon + +type AnyShapelyPolygon = Polygon | MultiPolygon +type GFComponent = Any # gdsfactory Component (avoid hard import) diff --git a/src/gsim/common/viz/__init__.py b/src/gsim/common/viz/__init__.py new file mode 100644 index 00000000..e472cc85 --- /dev/null +++ b/src/gsim/common/viz/__init__.py @@ -0,0 +1,31 @@ +"""Common visualization utilities for 3D and 2D geometry rendering. + +All functions accept a ``GeometryModel`` (from ``gsim.common.geometry_model``) +and do **not** depend on any solver (meep, tidy3d, etc.). + +3D backends: + - PyVista (desktop, ``plot_prisms_3d``) + - Open3D + Plotly (Jupyter, ``plot_prisms_3d_open3d``) + - Three.js / FastAPI (browser, ``serve_threejs_visualization``) + +2D backends: + - matplotlib (``plot_prism_slices``) +""" + +from gsim.common.viz.render2d import plot_prism_slices +from gsim.common.viz.render3d import ( + create_web_export, + export_3d_mesh, + plot_prisms_3d, + plot_prisms_3d_open3d, + serve_threejs_visualization, +) + +__all__ = [ + "create_web_export", + "export_3d_mesh", + "plot_prism_slices", + "plot_prisms_3d", + "plot_prisms_3d_open3d", + "serve_threejs_visualization", +] diff --git a/src/gsim/common/viz/_colors.py b/src/gsim/common/viz/_colors.py new file mode 100644 index 00000000..d2fdffd7 --- /dev/null +++ b/src/gsim/common/viz/_colors.py @@ -0,0 +1,53 @@ +"""Color and layer utilities shared across 3D rendering backends. + +These helpers are consumed internally by the render3d_* modules. +""" + +from __future__ import annotations + + +def generate_layer_colors( + layer_names: list[str], +) -> dict[str, tuple[float, float, float]]: + """Generate distinct RGB colours for each layer using a matplotlib colormap. + + Returns: + Mapping from layer name to (r, g, b) with values in [0, 1]. + """ + import matplotlib.pyplot as plt + + cmap = plt.cm.get_cmap("tab10" if len(layer_names) <= 10 else "tab20") + colors: dict[str, tuple[float, float, float]] = {} + for i, name in enumerate(layer_names): + rgb = cmap(i / max(len(layer_names) - 1, 1))[:3] + colors[name] = rgb + return colors + + +def generate_layer_colors_with_opacity( + layer_names: list[str], + layer_opacity: dict[str, float] | None = None, +) -> tuple[dict[str, list[float]], dict[str, float]]: + """Generate RGB colors and a separate opacity dict. + + Args: + layer_names: Ordered list of layer names. + layer_opacity: Optional per-layer opacity override. + Defaults to core=1.0, everything else=0.2. + + Returns: + (colors_dict, opacity_dict) where colors has [r, g, b] in [0, 1] + and opacity has float in [0, 1]. + """ + import matplotlib.pyplot as plt + + if layer_opacity is None: + layer_opacity = {name: 1.0 if name == "core" else 0.2 for name in layer_names} + + cmap = plt.cm.get_cmap("tab10" if len(layer_names) <= 10 else "tab20") + colors: dict[str, list[float]] = {} + for i, name in enumerate(layer_names): + rgb = cmap(i / max(len(layer_names) - 1, 1))[:3] + colors[name] = list(rgb) + + return colors, layer_opacity diff --git a/src/gsim/common/viz/_mesh_helpers.py b/src/gsim/common/viz/_mesh_helpers.py new file mode 100644 index 00000000..690b579c --- /dev/null +++ b/src/gsim/common/viz/_mesh_helpers.py @@ -0,0 +1,74 @@ +"""Low-level mesh construction helpers for 3D rendering. + +These operate on generic Prism objects and produce numpy vertex/face arrays +(or PyVista/Open3D meshes) without any solver-specific code. +""" + +from __future__ import annotations + +import numpy as np + +from gsim.common.geometry_model import GeometryModel, Prism + +# --------------------------------------------------------------------------- +# Generic vertex helpers +# --------------------------------------------------------------------------- + + +def prism_base_top_vertices( + prism: Prism, +) -> tuple[np.ndarray, np.ndarray]: + """Return (base_vertices, top_vertices) as (N, 3) numpy arrays.""" + xy = prism.vertices # (N, 2) + n = len(xy) + base = np.column_stack([xy, np.full(n, prism.z_base)]) + top = np.column_stack([xy, np.full(n, prism.z_top)]) + return base, top + + +# --------------------------------------------------------------------------- +# Simulation-box corner/line helpers +# --------------------------------------------------------------------------- + + +def simulation_box_corners( + geometry_model: GeometryModel, +) -> tuple[np.ndarray, list[list[int]]]: + """Return 8 corners and 12 edge-index pairs for the simulation bbox. + + Returns: + (points, lines) where points is (8, 3) and lines is a list of + [i, j] index pairs. + """ + mn = np.array(geometry_model.bbox[0]) + mx = np.array(geometry_model.bbox[1]) + + points = np.array( + [ + mn, + [mx[0], mn[1], mn[2]], + [mx[0], mx[1], mn[2]], + [mn[0], mx[1], mn[2]], + [mn[0], mn[1], mx[2]], + [mx[0], mn[1], mx[2]], + mx, + [mn[0], mx[1], mx[2]], + ] + ) + + lines = [ + [0, 1], + [1, 2], + [2, 3], + [3, 0], # bottom + [4, 5], + [5, 6], + [6, 7], + [7, 4], # top + [0, 4], + [1, 5], + [2, 6], + [3, 7], # verticals + ] + + return points, lines diff --git a/src/gsim/common/viz/render2d.py b/src/gsim/common/viz/render2d.py new file mode 100644 index 00000000..c5a397a7 --- /dev/null +++ b/src/gsim/common/viz/render2d.py @@ -0,0 +1,947 @@ +"""2D cross-sectional plotting for GeometryModel. + +Provides matplotlib-based 2D slicing views (XY, XZ, YZ) of a generic +GeometryModel without any solver-specific imports. + +When a ``SimOverlay`` is provided, the plot also draws: +- Simulation cell boundary (dashed black) +- PML regions (semi-transparent orange) +- Port markers (source in red, monitors in blue) +""" + +from __future__ import annotations + +from typing import Any + +import matplotlib.pyplot as plt +import numpy as np +from matplotlib.patches import Polygon, Rectangle + +from gsim.common.geometry_model import GeometryModel + +# --------------------------------------------------------------------------- +# Public API +# --------------------------------------------------------------------------- + + +def plot_prism_slices( + geometry_model: GeometryModel, + x: float | str | None = None, + y: float | str | None = None, + z: float | str = "core", + ax: plt.Axes | None = None, + legend: bool = True, + slices: str = "z", + *, + overlay: Any | None = None, +) -> plt.Axes | None: + """Plot cross sections of a GeometryModel with multi-view support. + + Args: + geometry_model: GeometryModel with prisms and bbox. + x: X-coordinate (or layer name) for the slice plane. + y: Y-coordinate (or layer name) for the slice plane. + z: Z-coordinate (or layer name) for the slice plane. + ax: Axes to draw on. If ``None``, a new figure is created. + legend: Whether to show the legend. + slices: Which slice(s) to plot -- "x", "y", "z", or combinations + like "xy", "xz", "yz", "xyz". + overlay: Optional SimOverlay with sim cell / PML / port metadata. + + Returns: + ``plt.Axes`` when *ax* was provided, otherwise ``None`` + (the figure is shown directly). + """ + slices_to_plot = sorted(set(slices.lower())) + if not all(s in "xyz" for s in slices_to_plot): + raise ValueError(f"slices must only contain 'x', 'y', 'z'. Got: {slices}") + + if ax is not None: + if len(slices_to_plot) > 1: + raise ValueError("Cannot plot multiple slices when ax is provided") + slice_axis = slices_to_plot[0] + if slice_axis == "x": + x_val = x if x is not None else "core" + return _plot_single_prism_slice( + geometry_model, + x=x_val, + y=None, + z=None, + ax=ax, + legend=legend, + overlay=overlay, + ) + if slice_axis == "y": + y_val = y if y is not None else "core" + return _plot_single_prism_slice( + geometry_model, + x=None, + y=y_val, + z=None, + ax=ax, + legend=legend, + overlay=overlay, + ) + if slice_axis == "z": + return _plot_single_prism_slice( + geometry_model, + x=None, + y=None, + z=z, + ax=ax, + legend=legend, + overlay=overlay, + ) + + _plot_multi_view( + geometry_model, + slices_to_plot, + x, + y, + z, + show_legend=legend, + overlay=overlay, + ) + return None + + +# --------------------------------------------------------------------------- +# Multi-view helper +# --------------------------------------------------------------------------- + + +def _plot_multi_view( + geometry_model: GeometryModel, + slices_to_plot: list[str], + x: float | str | None, + y: float | str | None, + z: float | str | None, + show_legend: bool = True, + overlay: Any | None = None, +) -> None: + """Create multi-view plot with a shared legend panel.""" + num_plots = len(slices_to_plot) + fig = plt.figure(constrained_layout=True) + gs = fig.add_gridspec(ncols=2, nrows=num_plots, width_ratios=(3, 1)) + + axes: list[plt.Axes] = [fig.add_subplot(gs[i, 0]) for i in range(num_plots)] + + for ax_i, slice_axis in zip(axes, slices_to_plot, strict=True): + if slice_axis == "x": + x_val = x if x is not None else "core" + _plot_single_prism_slice( + geometry_model, + x=x_val, + y=None, + z=None, + ax=ax_i, + legend=False, + overlay=overlay, + ) + elif slice_axis == "y": + y_val = y if y is not None else "core" + _plot_single_prism_slice( + geometry_model, + x=None, + y=y_val, + z=None, + ax=ax_i, + legend=False, + overlay=overlay, + ) + elif slice_axis == "z": + _plot_single_prism_slice( + geometry_model, + x=None, + y=None, + z=z, + ax=ax_i, + legend=False, + overlay=overlay, + ) + + if show_legend: + all_handles: list = [] + all_labels: list[str] = [] + seen: set[str] = set() + + for ax in axes: + handles, labels = ax.get_legend_handles_labels() + for handle, label in zip(handles, labels, strict=True): + if label not in seen: + all_handles.append(handle) + all_labels.append(label) + seen.add(label) + + legend_row = num_plots // 2 + axl = fig.add_subplot(gs[legend_row, 1]) + if all_handles: + axl.legend(all_handles, all_labels, loc="center") + axl.axis("off") + + plt.show() + + +# --------------------------------------------------------------------------- +# Single-slice renderer +# --------------------------------------------------------------------------- + + +def _plot_single_prism_slice( + geometry_model: GeometryModel, + x: float | str | None = None, + y: float | str | None = None, + z: float | str | None = None, + ax: plt.Axes | None = None, + legend: bool = True, + overlay: Any | None = None, +) -> plt.Axes: + """Plot a single cross-section using generic Prisms.""" + if ax is None: + _, ax = plt.subplots() + + # Resolve layer-name shortcuts to coordinates + x, y, z = ( + geometry_model.get_layer_center(c)[i] if isinstance(c, str) else c + for i, c in enumerate((x, y, z)) + ) + + active_count = sum([x is not None, y is not None, z is not None]) + if active_count != 1: + raise ValueError("Specify exactly one of x, y, or z for the slice plane") + + # Layer colours + layer_names = geometry_model.layer_names + colors = dict( + zip( + layer_names, + plt.colormaps.get_cmap("tab10")( + np.linspace(0, 1, max(len(layer_names), 1)) + ), + strict=True, + ) + ) + + # Compute total polygon area per layer for draw-order tiebreaking + # and alpha assignment. Larger area → background → draw first / more + # transparent. Smaller area → patterned → draw last / most opaque. + def _layer_area(name: str) -> float: + if name not in geometry_model.prisms: + return float("inf") + return sum( + p.original_polygon.area + for p in geometry_model.prisms[name] + if p.original_polygon is not None + ) + + layer_areas: dict[str, float] = {n: _layer_area(n) for n in layer_names} + + # Sort layers for drawing: lower mesh_order first (background), then + # larger polygon area first (background). This ensures patterned + # layers (small area, like waveguide cores) are drawn on top of + # fill layers (large area, like cladding rectangles). + sorted_names = sorted( + layer_names, + key=lambda n: ( + geometry_model.layer_mesh_orders.get(n, 0), + -layer_areas.get(n, 0.0), + ), + ) + + # Mesh-order → zorder: lower mesh_order = background = lower zorder + mesh_orders = np.unique( + [geometry_model.layer_mesh_orders.get(n, 0) for n in sorted_names] + ) + order_map = dict(zip(mesh_orders, range(0, -len(mesh_orders), -1), strict=True)) + + # Assign a small per-layer zorder bump so that within the same + # mesh_order, layers with smaller area (patterned) are on top. + _names_by_mesh: dict[int, list[str]] = {} + for n in sorted_names: + mo = geometry_model.layer_mesh_orders.get(n, 0) + _names_by_mesh.setdefault(mo, []).append(n) + + _layer_zorder_bump: dict[str, float] = {} + for names_in_group in _names_by_mesh.values(): + for rank, n in enumerate(names_in_group): + _layer_zorder_bump[n] = rank * 0.01 + + # Alpha mapping: patterned layers (smaller area) get higher alpha, + # background layers (larger area) get lower alpha. + max_area = max(layer_areas.values()) if layer_areas else 1.0 + _layer_alpha: dict[str, float] = {} + for n in layer_names: + ratio = layer_areas.get(n, max_area) / max_area if max_area > 0 else 1.0 + # Map: ratio 1.0 (largest area) → alpha 0.35, ratio ~0 → alpha 0.90 + _layer_alpha[n] = 0.90 - 0.55 * ratio + + xmin, xmax = np.inf, -np.inf + ymin, ymax = np.inf, -np.inf + + # First pass: compute axis limits from layer bboxes + for name in sorted_names: + try: + bbox = geometry_model.get_layer_bbox(name) + if z is not None: + xmin = min(xmin, bbox[0][0]) + xmax = max(xmax, bbox[1][0]) + ymin = min(ymin, bbox[0][1]) + ymax = max(ymax, bbox[1][1]) + elif x is not None: + ymin = min(ymin, bbox[0][1]) + ymax = max(ymax, bbox[1][1]) + elif y is not None: + xmin = min(xmin, bbox[0][0]) + xmax = max(xmax, bbox[1][0]) + except Exception: # noqa: S112 + continue + + # Second pass: draw patches + for name in sorted_names: + base_color = colors.get(name, "lightgray") + mesh_order = geometry_model.layer_mesh_orders.get(name, 0) + bbox = geometry_model.get_layer_bbox(name) + layer_zorder = order_map.get(mesh_order, 0) + _layer_zorder_bump.get(name, 0) + alpha = _layer_alpha.get(name, 0.7) + + # Build RGBA face colour with per-layer alpha + if hasattr(base_color, "__len__") and len(base_color) >= 3: + fc = (*base_color[:3], alpha) + else: + fc = base_color # fallback for named colours + + if z is not None: + z_min_layer = bbox[0][2] + z_max_layer = bbox[1][2] + if not (z_min_layer <= z <= z_max_layer): + continue + + if name in geometry_model.prisms: + prisms = geometry_model.prisms[name] + for idx, prism in enumerate(prisms): + if not (prism.z_base <= z <= prism.z_top): + continue + + xy_points = prism.vertices.tolist() + patch = Polygon( + xy_points, + facecolor=fc, + edgecolor="k", + linewidth=0.5, + label=name if idx == 0 else None, + zorder=layer_zorder, + ) + ax.add_patch(patch) + else: + rect = Rectangle( + (bbox[0][0], bbox[0][1]), + bbox[1][0] - bbox[0][0], + bbox[1][1] - bbox[0][1], + facecolor=fc, + edgecolor="k", + linewidth=0.5, + label=name, + zorder=layer_zorder, + ) + ax.add_patch(rect) + + elif x is not None: + # Intersect each prism polygon with the x=const plane + # to get actual y-ranges, then draw rectangles in YZ + drawn = _draw_prism_cross_sections( + ax, + geometry_model, + name, + "x", + x, + fc, + layer_zorder, + ) + if not drawn: + # Fallback to bbox if no prism data + rect = Rectangle( + (bbox[0][1], bbox[0][2]), + bbox[1][1] - bbox[0][1], + bbox[1][2] - bbox[0][2], + facecolor=fc, + edgecolor="k", + linewidth=0.5, + label=name, + zorder=layer_zorder, + ) + ax.add_patch(rect) + + elif y is not None: + # Intersect each prism polygon with the y=const plane + # to get actual x-ranges, then draw rectangles in XZ + drawn = _draw_prism_cross_sections( + ax, + geometry_model, + name, + "y", + y, + fc, + layer_zorder, + ) + if not drawn: + # Fallback to bbox if no prism data + rect = Rectangle( + (bbox[0][0], bbox[0][2]), + bbox[1][0] - bbox[0][0], + bbox[1][2] - bbox[0][2], + facecolor=fc, + edgecolor="k", + linewidth=0.5, + label=name, + zorder=layer_zorder, + ) + ax.add_patch(rect) + + # Axis labels and simulation-box outline + size = list(geometry_model.size) + cmin = list(geometry_model.bbox[0]) + + if z is not None: + size = size[:2] + cmin = cmin[:2] + xlabel, ylabel = "x (um)", "y (um)" + ax.set_title(f"XY cross section at z={z:.2f}") + elif x is not None: + size = [size[1], size[2]] + cmin = [cmin[1], cmin[2]] + xlabel, ylabel = "y (um)", "z (um)" + ax.set_title(f"YZ cross section at x={x:.2f}") + xmin, xmax = cmin[0], cmin[0] + size[0] + ymin, ymax = cmin[1], cmin[1] + size[1] + elif y is not None: + size = [size[0], size[2]] + cmin = [cmin[0], cmin[2]] + xlabel, ylabel = "x (um)", "z (um)" + ax.set_title(f"XZ cross section at y={y:.2f}") + ymin, ymax = cmin[1], cmin[1] + size[1] + + # Draw overlay or fallback geometry bbox + if overlay is not None: + _draw_overlay(ax, overlay, x=x, y=y, z=z) + # Expand axis limits to include full sim cell + if z is not None: + xmin = min(xmin, overlay.cell_min[0]) + xmax = max(xmax, overlay.cell_max[0]) + ymin = min(ymin, overlay.cell_min[1]) + ymax = max(ymax, overlay.cell_max[1]) + elif x is not None: + xmin = min(xmin, overlay.cell_min[1]) + xmax = max(xmax, overlay.cell_max[1]) + ymin = min(ymin, overlay.cell_min[2]) + ymax = max(ymax, overlay.cell_max[2]) + elif y is not None: + xmin = min(xmin, overlay.cell_min[0]) + xmax = max(xmax, overlay.cell_max[0]) + ymin = min(ymin, overlay.cell_min[2]) + ymax = max(ymax, overlay.cell_max[2]) + else: + sim_roi = Rectangle( + tuple(cmin), # type: ignore[arg-type] + *size, + facecolor="none", + edgecolor="k", + linestyle="--", + linewidth=1, + label="Simulation", + ) + ax.add_patch(sim_roi) + + ax.set_xlabel(xlabel) + ax.set_ylabel(ylabel) + ax.set_xlim(xmin, xmax) + ax.set_ylim(ymin, ymax) + ax.set_aspect("equal") + + if legend: + ax.legend(fancybox=True, framealpha=1.0) + + return ax + + +# --------------------------------------------------------------------------- +# Prism cross-section helpers (for x/y slice views) +# --------------------------------------------------------------------------- + + +def _draw_prism_cross_sections( + ax: plt.Axes, + geometry_model: GeometryModel, + layer_name: str, + slice_axis: str, + slice_coord: float, + fc: Any, + layer_zorder: float, +) -> bool: + """Intersect prism polygons with a slice plane, draw rectangles. + + For a y-slice at y=y0, each prism polygon is intersected with the + line y=y0 to find x-segments. Each segment becomes a rectangle in + XZ space: (x_start, z_base) to (x_end, z_top). + + For an x-slice at x=x0, similarly finds y-segments → YZ rectangles. + + Returns True if at least one patch was drawn. + """ + from shapely.geometry import LineString, MultiLineString + from shapely.geometry import Polygon as ShapelyPolygon + + if layer_name not in geometry_model.prisms: + return False + + prisms = geometry_model.prisms[layer_name] + if not prisms: + return False + + drawn = False + first = True + + for prism in prisms: + poly = ShapelyPolygon(prism.vertices) + if poly.is_empty or not poly.is_valid: + continue + + # Build a long line through the polygon at the slice coordinate + bounds = poly.bounds # (minx, miny, maxx, maxy) + margin = 1.0 + if slice_axis == "y": + line = LineString( + [ + (bounds[0] - margin, slice_coord), + (bounds[2] + margin, slice_coord), + ] + ) + else: # x-slice + line = LineString( + [ + (slice_coord, bounds[1] - margin), + (slice_coord, bounds[3] + margin), + ] + ) + + intersection = poly.intersection(line) + if intersection.is_empty: + continue + + # Collect line segments from the intersection + segments = [] + if isinstance(intersection, LineString): + segments.append(intersection) + elif isinstance(intersection, MultiLineString): + segments.extend(intersection.geoms) + + for seg in segments: + coords = list(seg.coords) + if len(coords) < 2: + continue + + if slice_axis == "y": + # XZ view: horizontal axis = x, vertical = z + h_vals = [c[0] for c in coords] + h_min, h_max = min(h_vals), max(h_vals) + else: + # YZ view: horizontal axis = y, vertical = z + h_vals = [c[1] for c in coords] + h_min, h_max = min(h_vals), max(h_vals) + + if h_max - h_min < 1e-9: + continue + + rect = Rectangle( + (h_min, prism.z_base), + h_max - h_min, + prism.z_top - prism.z_base, + facecolor=fc, + edgecolor="k", + linewidth=0.5, + label=layer_name if first else None, + zorder=layer_zorder, + ) + ax.add_patch(rect) + drawn = True + first = False + + return drawn + + +# --------------------------------------------------------------------------- +# Overlay drawing +# --------------------------------------------------------------------------- + +_PML_COLOR = (1.0, 0.6, 0.0, 0.20) # semi-transparent orange +_PML_EDGE = (1.0, 0.6, 0.0, 0.50) +_SRC_COLOR = "red" +_MON_COLOR = "royalblue" + +# Dielectric background colours — keyed by lowercase material name +_DIELECTRIC_COLORS: dict[str, tuple[float, float, float, float]] = { + "sio2": (0.68, 0.85, 0.90, 0.25), # light blue + "silicon": (0.60, 0.60, 0.65, 0.30), # warm grey + "si": (0.60, 0.60, 0.65, 0.30), + "air": (1.00, 1.00, 1.00, 0.08), # near-transparent +} +_DIELECTRIC_DEFAULT_COLOR = (0.80, 0.80, 0.80, 0.20) # light grey fallback +_DIELECTRIC_ZORDER = -10 + + +def _draw_overlay( + ax: plt.Axes, + overlay: Any, + *, + x: float | None, + y: float | None, + z: float | None, +) -> None: + """Draw dielectric backgrounds, cell boundary, PML regions, and port markers.""" + cmin = overlay.cell_min + cmax = overlay.cell_max + pml = overlay.dpml + + # Draw dielectric background slabs first (lowest zorder) + dielectrics = getattr(overlay, "dielectrics", []) + if dielectrics: + if z is not None: + _draw_dielectrics_xy(ax, dielectrics, z, cmin, cmax) + elif x is not None: + _draw_dielectrics_side(ax, dielectrics, cmin[1], cmax[1], "yz") + elif y is not None: + _draw_dielectrics_side(ax, dielectrics, cmin[0], cmax[0], "xz") + + if z is not None: + _draw_overlay_xy(ax, cmin, cmax, pml, overlay.ports) + elif x is not None: + _draw_overlay_yz(ax, cmin, cmax, pml, overlay.ports, x) + elif y is not None: + _draw_overlay_xz(ax, cmin, cmax, pml, overlay.ports, y) + + +def _draw_overlay_xy( + ax: plt.Axes, + cmin: tuple[float, float, float], + cmax: tuple[float, float, float], + pml: float, + ports: list, +) -> None: + """Draw overlay elements for an XY (z-slice) view.""" + x0, y0 = cmin[0], cmin[1] + x1, y1 = cmax[0], cmax[1] + w, h = x1 - x0, y1 - y0 + + # Sim cell boundary + ax.add_patch( + Rectangle( + (x0, y0), + w, + h, + facecolor="none", + edgecolor="k", + linestyle="--", + linewidth=1, + label="Sim cell", + zorder=90, + ) + ) + + # PML rectangles (4 edges) + _add_pml_rect(ax, x0, y0, pml, h, label="PML") # left + _add_pml_rect(ax, x1 - pml, y0, pml, h) # right + _add_pml_rect(ax, x0 + pml, y0, w - 2 * pml, pml) # bottom + _add_pml_rect(ax, x0 + pml, y1 - pml, w - 2 * pml, pml) # top + + # Ports + labeled: set[str] = set() + for port in ports: + cx, cy, _cz = port.center + color = _SRC_COLOR if port.is_source else _MON_COLOR + legend_key = "Source" if port.is_source else "Monitor" + label = legend_key if legend_key not in labeled else None + labeled.add(legend_key) + hw = port.width / 2 + + if port.normal_axis == 0: # x-normal → vertical line + ax.plot( + [cx, cx], + [cy - hw, cy + hw], + color=color, + linewidth=2, + zorder=95, + label=label, + ) + dx = 0.15 if port.direction == "+" else -0.15 + ax.annotate( + "", + xy=(cx + dx, cy), + xytext=(cx, cy), + arrowprops=dict(arrowstyle="->", color=color, lw=1.5), + zorder=96, + ) + else: # y-normal → horizontal line + ax.plot( + [cx - hw, cx + hw], + [cy, cy], + color=color, + linewidth=2, + zorder=95, + label=label, + ) + dy = 0.15 if port.direction == "+" else -0.15 + ax.annotate( + "", + xy=(cx, cy + dy), + xytext=(cx, cy), + arrowprops=dict(arrowstyle="->", color=color, lw=1.5), + zorder=96, + ) + + ax.annotate( + port.name, + (cx, cy), + fontsize=7, + ha="center", + va="bottom", + color=color, + zorder=96, + xytext=(0, 4), + textcoords="offset points", + ) + + +def _draw_overlay_yz( + ax: plt.Axes, + cmin: tuple[float, float, float], + cmax: tuple[float, float, float], + pml: float, + ports: list, + x_slice: float, +) -> None: + """Draw overlay elements for a YZ (x-slice) view.""" + y0, z0 = cmin[1], cmin[2] + y1, z1 = cmax[1], cmax[2] + w, h = y1 - y0, z1 - z0 + + # Sim cell boundary + ax.add_patch( + Rectangle( + (y0, z0), + w, + h, + facecolor="none", + edgecolor="k", + linestyle="--", + linewidth=1, + label="Sim cell", + zorder=90, + ) + ) + + # PML rectangles + _add_pml_rect(ax, y0, z0, pml, h, label="PML") # left + _add_pml_rect(ax, y1 - pml, z0, pml, h) # right + _add_pml_rect(ax, y0 + pml, z0, w - 2 * pml, pml) # bottom + _add_pml_rect(ax, y0 + pml, z1 - pml, w - 2 * pml, pml) # top + + # Ports that intersect this x-slice (x-normal ports at x_slice) + labeled: set[str] = set() + for port in ports: + cx, cy, cz = port.center + if port.normal_axis == 0 and abs(cx - x_slice) < 0.01: + color = _SRC_COLOR if port.is_source else _MON_COLOR + legend_key = "Source" if port.is_source else "Monitor" + label = legend_key if legend_key not in labeled else None + labeled.add(legend_key) + hw = port.width / 2 + hz = port.z_span / 2 + ax.add_patch( + Rectangle( + (cy - hw, cz - hz), + port.width, + port.z_span, + facecolor="none", + edgecolor=color, + linewidth=1.5, + label=label, + zorder=95, + ) + ) + ax.annotate( + port.name, + (cy, cz + hz), + fontsize=7, + ha="center", + va="bottom", + color=color, + zorder=96, + ) + + +def _draw_overlay_xz( + ax: plt.Axes, + cmin: tuple[float, float, float], + cmax: tuple[float, float, float], + pml: float, + ports: list, + y_slice: float, +) -> None: + """Draw overlay elements for an XZ (y-slice) view.""" + x0, z0 = cmin[0], cmin[2] + x1, z1 = cmax[0], cmax[2] + w, h = x1 - x0, z1 - z0 + + # Sim cell boundary + ax.add_patch( + Rectangle( + (x0, z0), + w, + h, + facecolor="none", + edgecolor="k", + linestyle="--", + linewidth=1, + label="Sim cell", + zorder=90, + ) + ) + + # PML rectangles + _add_pml_rect(ax, x0, z0, pml, h, label="PML") # left + _add_pml_rect(ax, x1 - pml, z0, pml, h) # right + _add_pml_rect(ax, x0 + pml, z0, w - 2 * pml, pml) # bottom + _add_pml_rect(ax, x0 + pml, z1 - pml, w - 2 * pml, pml) # top + + # Ports that intersect this y-slice (y-normal ports at y_slice) + labeled: set[str] = set() + for port in ports: + cx, cy, cz = port.center + if port.normal_axis == 1 and abs(cy - y_slice) < 0.01: + color = _SRC_COLOR if port.is_source else _MON_COLOR + legend_key = "Source" if port.is_source else "Monitor" + label = legend_key if legend_key not in labeled else None + labeled.add(legend_key) + hw = port.width / 2 + hz = port.z_span / 2 + ax.add_patch( + Rectangle( + (cx - hw, cz - hz), + port.width, + port.z_span, + facecolor="none", + edgecolor=color, + linewidth=1.5, + label=label, + zorder=95, + ) + ) + ax.annotate( + port.name, + (cx, cz + hz), + fontsize=7, + ha="center", + va="bottom", + color=color, + zorder=96, + ) + + +def _add_pml_rect( + ax: plt.Axes, + x: float, + y: float, + w: float, + h: float, + label: str | None = None, +) -> None: + """Add a single semi-transparent PML rectangle.""" + if w <= 0 or h <= 0: + return + ax.add_patch( + Rectangle( + (x, y), + w, + h, + facecolor=_PML_COLOR, + edgecolor=_PML_EDGE, + linewidth=0.5, + label=label, + zorder=80, + ) + ) + + +# --------------------------------------------------------------------------- +# Dielectric background drawing +# --------------------------------------------------------------------------- + + +def _diel_color(material: str) -> tuple[float, float, float, float]: + """Look up the fill colour for a dielectric material name.""" + return _DIELECTRIC_COLORS.get(material.lower(), _DIELECTRIC_DEFAULT_COLOR) + + +def _draw_dielectrics_side( + ax: plt.Axes, + dielectrics: list, + h_min: float, + h_max: float, + _axis: str = "xz", +) -> None: + """Draw dielectric background bands for a side view (XZ or YZ). + + Each dielectric is a horizontal band spanning the full horizontal extent + of the simulation cell at the dielectric's z-range. + """ + labeled: set[str] = set() + for diel in dielectrics: + zmin, zmax = diel.zmin, diel.zmax + if zmax <= zmin: + continue + color = _diel_color(diel.material) + label = diel.name if diel.name not in labeled else None + labeled.add(diel.name) + ax.add_patch( + Rectangle( + (h_min, zmin), + h_max - h_min, + zmax - zmin, + facecolor=color, + edgecolor="none", + label=label, + zorder=_DIELECTRIC_ZORDER, + ) + ) + + +def _draw_dielectrics_xy( + ax: plt.Axes, + dielectrics: list, + z_slice: float, + cmin: tuple[float, float, float], + cmax: tuple[float, float, float], +) -> None: + """Draw dielectric background for an XY (z-slice) view. + + If the z-slice is within a dielectric slab, fill the entire cell with + that material's colour. If multiple slabs overlap, draw all of them + (later ones on top). + """ + labeled: set[str] = set() + for diel in dielectrics: + if not (diel.zmin <= z_slice <= diel.zmax): + continue + color = _diel_color(diel.material) + label = diel.name if diel.name not in labeled else None + labeled.add(diel.name) + ax.add_patch( + Rectangle( + (cmin[0], cmin[1]), + cmax[0] - cmin[0], + cmax[1] - cmin[1], + facecolor=color, + edgecolor="none", + label=label, + zorder=_DIELECTRIC_ZORDER, + ) + ) diff --git a/src/gsim/common/viz/render3d.py b/src/gsim/common/viz/render3d.py new file mode 100644 index 00000000..ae71fbcc --- /dev/null +++ b/src/gsim/common/viz/render3d.py @@ -0,0 +1,23 @@ +"""Unified 3D rendering API for GeometryModel. + +Re-exports the public functions from the backend-specific modules so that +callers can simply do:: + + from gsim.common.viz.render3d import plot_prisms_3d, plot_prisms_3d_open3d +""" + +from gsim.common.viz.render3d_open3d import plot_prisms_3d_open3d +from gsim.common.viz.render3d_pyvista import ( + create_web_export, + export_3d_mesh, + plot_prisms_3d, +) +from gsim.common.viz.render3d_threejs import serve_threejs_visualization + +__all__ = [ + "create_web_export", + "export_3d_mesh", + "plot_prisms_3d", + "plot_prisms_3d_open3d", + "serve_threejs_visualization", +] diff --git a/src/gsim/common/viz/render3d_open3d.py b/src/gsim/common/viz/render3d_open3d.py new file mode 100644 index 00000000..528211b3 --- /dev/null +++ b/src/gsim/common/viz/render3d_open3d.py @@ -0,0 +1,441 @@ +"""Open3D / Plotly-based 3D rendering for GeometryModel. + +Provides Jupyter-compatible interactive visualisation via Open3D meshes +displayed through Plotly. +""" + +from __future__ import annotations + +from typing import Any + +import numpy as np + +from gsim.common.geometry_model import GeometryModel, Prism +from gsim.common.viz._colors import generate_layer_colors_with_opacity +from gsim.common.viz._mesh_helpers import ( + prism_base_top_vertices, + simulation_box_corners, +) + +try: + import open3d as o3d + + OPEN3D_AVAILABLE = True +except ImportError: + OPEN3D_AVAILABLE = False + + +# --------------------------------------------------------------------------- +# Public API +# --------------------------------------------------------------------------- + + +def plot_prisms_3d_open3d( + geometry_model: GeometryModel, + *, + show_edges: bool = False, + color_by_layer: bool = True, + show_simulation_box: bool = True, + notebook: bool = True, + layer_opacity: dict[str, float] | None = None, + **kwargs: Any, +) -> None: + """Create interactive 3D visualisation using Open3D + Plotly. + + Args: + geometry_model: GeometryModel containing prisms and bbox. + show_edges: Show wireframe edges. + color_by_layer: Colour each layer differently. + show_simulation_box: Draw the simulation box. + notebook: Display inside Jupyter notebook. + layer_opacity: Per-layer opacity override (default: core=1.0, else 0.2). + **kwargs: Extra Plotly figure options. + """ + if not OPEN3D_AVAILABLE: + raise ImportError("Open3D is required. Install with: pip install open3d") + + try: + import plotly.graph_objects as go + except ImportError as err: + raise ImportError( + "Plotly is required for Open3D notebook visualization. " + "Install with: pip install plotly" + ) from err + + layer_meshes = _convert_prisms_to_open3d(geometry_model) + colors, opacity_dict = generate_layer_colors_with_opacity( + list(layer_meshes.keys()), layer_opacity + ) + + plotly_meshes: list[Any] = [] + + for layer_name, meshes in layer_meshes.items(): + layer_color = colors[layer_name] if color_by_layer else [0.7, 0.7, 0.7] + layer_opacity_val = opacity_dict.get(layer_name, 0.8) + + for i, mesh in enumerate(meshes): + if color_by_layer: + mesh.paint_uniform_color(layer_color[:3]) + + plotly_meshes.append( + _mesh_to_mesh3d( + mesh, + opacity=layer_opacity_val, + name=f"{layer_name}_{i}", + color=layer_color[:3] if color_by_layer else [0.7, 0.7, 0.7], + ) + ) + + if show_edges: + plotly_meshes.append( + _wireframe_to_scatter3d(mesh, name=f"{layer_name}_edges_{i}") + ) + + if show_simulation_box: + plotly_meshes.append(_create_simulation_box_plotly(geometry_model)) + + fig = go.Figure(data=plotly_meshes) + + # Compute axis ranges for zoom + all_x: list[float] = [] + all_y: list[float] = [] + all_z: list[float] = [] + for meshes in layer_meshes.values(): + for mesh in meshes: + verts = np.asarray(mesh.vertices) + all_x.extend(verts[:, 0]) + all_y.extend(verts[:, 1]) + all_z.extend(verts[:, 2]) + + if all_x: + cx = np.mean([min(all_x), max(all_x)]) + cy = np.mean([min(all_y), max(all_y)]) + cz = np.mean([min(all_z), max(all_z)]) + rs = max( + max(all_x) - min(all_x), + max(all_y) - min(all_y), + max(all_z) - min(all_z), + ) + else: + cx = cy = cz = 0.0 + rs = 10.0 + + fig.update_layout( + scene=dict( + xaxis_title="X (um)", + yaxis_title="Y (um)", + zaxis_title="Z (um)", + aspectmode="data", + camera=dict( + eye=dict(x=1.5, y=1.5, z=1.5), + center=dict(x=0, y=0, z=0), + projection=dict(type="perspective"), + ), + xaxis=dict(range=[cx - rs * 0.6, cx + rs * 0.6]), + yaxis=dict(range=[cy - rs * 0.6, cy + rs * 0.6]), + zaxis=dict(range=[cz - rs * 0.6, cz + rs * 0.6]), + ), + title="3D Geometry Visualization", + dragmode="orbit", + **kwargs, + ) + + config = { + "scrollZoom": True, + "doubleClick": "reset+autosize", + "modeBarButtonsToRemove": ["pan2d", "lasso2d"], + "displayModeBar": True, + "responsive": True, + } + + if notebook: + fig.show(config=config) + else: + fig.write_html("geometry_3d.html", config=config) + import webbrowser + + webbrowser.open("geometry_3d.html") + + +# --------------------------------------------------------------------------- +# Open3D mesh conversion +# --------------------------------------------------------------------------- + + +def _convert_prisms_to_open3d( + geometry_model: GeometryModel, +) -> dict[str, list[Any]]: + """Convert generic Prisms to Open3D meshes, with triangle-merge optimisation.""" + layer_meshes: dict[str, list[Any]] = {} + + for layer_name, prisms in geometry_model.prisms.items(): + if not prisms: + continue + + triangular_count = sum(1 for p in prisms if len(p.vertices) == 3) + + if triangular_count > 100: + triangular = [p for p in prisms if len(p.vertices) == 3] + non_triangular = [p for p in prisms if len(p.vertices) != 3] + + meshes: list[Any] = [] + if triangular: + merged = _merge_triangular_prisms_o3d(triangular) + if merged is not None: + meshes.append(merged) + for prism in non_triangular: + base, top = prism_base_top_vertices(prism) + meshes.append(_create_prism_mesh_o3d(base, top, prism)) + else: + meshes = [] + for prism in prisms: + base, top = prism_base_top_vertices(prism) + meshes.append(_create_prism_mesh_o3d(base, top, prism)) + + layer_meshes[layer_name] = meshes + + return layer_meshes + + +def _merge_triangular_prisms_o3d(prisms: list[Prism]) -> Any: + """Merge many triangular prisms into one Open3D mesh for performance.""" + all_vertices: list[np.ndarray] = [] + all_triangles: list[list[int]] = [] + offset = 0 + + for prism in prisms: + base, top = prism_base_top_vertices(prism) + verts = np.vstack([base, top]) + all_vertices.append(verts) + n = len(base) + + # bottom + all_triangles.append([offset, offset + 2, offset + 1]) + # top + all_triangles.append([offset + n, offset + n + 1, offset + n + 2]) + # sides (2 tris per quad) + for i in range(n): + ni = (i + 1) % n + all_triangles.append([offset + i, offset + ni, offset + ni + n]) + all_triangles.append([offset + i, offset + ni + n, offset + i + n]) + + offset += len(verts) + + if all_vertices: + combined = np.vstack(all_vertices) + mesh = o3d.geometry.TriangleMesh() + mesh.vertices = o3d.utility.Vector3dVector(combined) + mesh.triangles = o3d.utility.Vector3iVector(all_triangles) + mesh.compute_vertex_normals() + return mesh + return None + + +def _create_prism_mesh_o3d( + base_vertices: np.ndarray, + top_vertices: np.ndarray, + prism: Prism | None = None, +) -> Any: + """Create Open3D mesh from base/top 3D vertices, with hole/concave support.""" + if prism is not None and prism.original_polygon is not None: + poly = prism.original_polygon + has_holes = hasattr(poly, "interiors") and poly.interiors + is_concave = ( + not has_holes + and hasattr(poly, "convex_hull") + and poly.convex_hull.area > poly.area * 1.001 + ) + if has_holes or is_concave: + return _create_prism_mesh_with_holes_o3d(poly, base_vertices, top_vertices) + + n = len(base_vertices) + all_verts = np.vstack([base_vertices, top_vertices]) + # bottom (fan) + faces: list[list[int]] = [[0, i + 1, i] for i in range(1, n - 1)] + # top (fan) + faces.extend([n, n + i, n + i + 1] for i in range(1, n - 1)) + # sides + for i in range(n): + ni = (i + 1) % n + faces.append([i, ni, ni + n]) + faces.append([i, ni + n, i + n]) + + mesh = o3d.geometry.TriangleMesh() + mesh.vertices = o3d.utility.Vector3dVector(all_verts) + mesh.triangles = o3d.utility.Vector3iVector(faces) + mesh.compute_vertex_normals() + return mesh + + +def _create_prism_mesh_with_holes_o3d( + shapely_polygon: Any, + base_vertices: np.ndarray, + top_vertices: np.ndarray, +) -> Any: + """Create Open3D mesh from a Shapely polygon with holes via Delaunay.""" + try: + import shapely.geometry as sg + from scipy.spatial import Delaunay + except ImportError: + return _create_prism_mesh_o3d(base_vertices, top_vertices) + + all_points: list[tuple[float, float]] = [] + boundary_segments: list[list[int]] = [] + + exterior_coords = list(shapely_polygon.exterior.coords[:-1]) + start_idx = 0 + all_points.extend(exterior_coords) + boundary_segments = [ + [start_idx + i, start_idx + (i + 1) % len(exterior_coords)] + for i in range(len(exterior_coords)) + ] + + for interior in shapely_polygon.interiors: + interior_coords = list(interior.coords[:-1]) + start_idx = len(all_points) + all_points.extend(interior_coords) + boundary_segments.extend( + [start_idx + i, start_idx + (i + 1) % len(interior_coords)] + for i in range(len(interior_coords)) + ) + + if len(all_points) < 3: + return _create_prism_mesh_o3d(base_vertices, top_vertices) + + points_2d = np.array(all_points) + tri = Delaunay(points_2d) + + valid_triangles = [ + simplex + for simplex in tri.simplices + if shapely_polygon.contains(sg.Point(*np.mean(points_2d[simplex], axis=0))) + ] + + if not valid_triangles: + return _create_prism_mesh_o3d(base_vertices, top_vertices) + + z_base = base_vertices[0, 2] if len(base_vertices) > 0 else 0 + z_top = top_vertices[0, 2] if len(top_vertices) > 0 else z_base + 1 + + verts_3d: list[list[float]] = [[pt[0], pt[1], z_base] for pt in points_2d] + [ + [pt[0], pt[1], z_top] for pt in points_2d + ] + + n_pts = len(points_2d) + faces_3d: list[list[int]] = [ + [tri_idx[0], tri_idx[1], tri_idx[2]] for tri_idx in valid_triangles + ] + faces_3d.extend( + [tri_idx[0] + n_pts, tri_idx[2] + n_pts, tri_idx[1] + n_pts] + for tri_idx in valid_triangles + ) + for seg in boundary_segments: + i, j = seg + faces_3d.append([i, j, j + n_pts]) + faces_3d.append([i, j + n_pts, i + n_pts]) + + try: + mesh = o3d.geometry.TriangleMesh() + mesh.vertices = o3d.utility.Vector3dVector(verts_3d) + mesh.triangles = o3d.utility.Vector3iVector(faces_3d) + mesh.compute_vertex_normals() + except Exception: + return _create_prism_mesh_o3d(base_vertices, top_vertices) + else: + return mesh + + +# --------------------------------------------------------------------------- +# Plotly helpers +# --------------------------------------------------------------------------- + + +def _mesh_to_mesh3d( + mesh: Any, + opacity: float = 1.0, + name: str = "", + color: list[float] | None = None, +) -> Any: + """Convert Open3D mesh to Plotly Mesh3d.""" + import plotly.graph_objects as go + + verts = np.asarray(mesh.vertices) + tris = np.asarray(mesh.triangles) + + if color is not None: + c = (np.array(color) * 255).astype(int) + color_str = f"rgb({c[0]},{c[1]},{c[2]})" + elif len(mesh.vertex_colors): + c = (np.asarray(mesh.vertex_colors)[0] * 255).astype(int) + color_str = f"rgb({c[0]},{c[1]},{c[2]})" + else: + color_str = "rgb(180,180,180)" + + return go.Mesh3d( + x=verts[:, 0], + y=verts[:, 1], + z=verts[:, 2], + i=tris[:, 0], + j=tris[:, 1], + k=tris[:, 2], + color=color_str, + opacity=opacity, + name=name, + showlegend=bool(name), + ) + + +def _wireframe_to_scatter3d(mesh: Any, name: str = "") -> Any: + """Convert Open3D mesh edges to Plotly Scatter3d.""" + import plotly.graph_objects as go + + line_set = o3d.geometry.LineSet.create_from_triangle_mesh(mesh) + points = np.asarray(line_set.points) + lines = np.asarray(line_set.lines) + + x_lines: list[float | None] = [] + y_lines: list[float | None] = [] + z_lines: list[float | None] = [] + + for line in lines: + p1, p2 = points[line[0]], points[line[1]] + x_lines.extend([p1[0], p2[0], None]) + y_lines.extend([p1[1], p2[1], None]) + z_lines.extend([p1[2], p2[2], None]) + + return go.Scatter3d( + x=x_lines, + y=y_lines, + z=z_lines, + mode="lines", + line=dict(color="black", width=2), + name=name, + showlegend=bool(name), + ) + + +def _create_simulation_box_plotly(geometry_model: GeometryModel) -> Any: + """Create a Plotly Scatter3d wireframe box for the simulation bbox.""" + import plotly.graph_objects as go + + points, lines = simulation_box_corners(geometry_model) + + x_lines: list[float | None] = [] + y_lines: list[float | None] = [] + z_lines: list[float | None] = [] + + for line in lines: + p1, p2 = points[line[0]], points[line[1]] + x_lines.extend([p1[0], p2[0], None]) + y_lines.extend([p1[1], p2[1], None]) + z_lines.extend([p1[2], p2[2], None]) + + return go.Scatter3d( + x=x_lines, + y=y_lines, + z=z_lines, + mode="lines", + line=dict(color="black", width=2, dash="dot"), + name="Simulation Box", + showlegend=True, + ) diff --git a/src/gsim/common/viz/render3d_pyvista.py b/src/gsim/common/viz/render3d_pyvista.py new file mode 100644 index 00000000..f614e5d1 --- /dev/null +++ b/src/gsim/common/viz/render3d_pyvista.py @@ -0,0 +1,385 @@ +"""PyVista-based 3D rendering for GeometryModel. + +Provides desktop-quality interactive visualisation using the PyVista library. +""" + +from __future__ import annotations + +from typing import Any + +import numpy as np + +from gsim.common.geometry_model import GeometryModel, Prism +from gsim.common.viz._colors import generate_layer_colors +from gsim.common.viz._mesh_helpers import ( + prism_base_top_vertices, +) + +try: + import pyvista as pv + + PYVISTA_AVAILABLE = True +except ImportError: + PYVISTA_AVAILABLE = False + + +# --------------------------------------------------------------------------- +# Public API +# --------------------------------------------------------------------------- + + +def plot_prisms_3d( + geometry_model: GeometryModel, + *, + show_edges: bool = True, + opacity: float = 0.8, # noqa: ARG001 + color_by_layer: bool = True, + show_simulation_box: bool = True, + camera_position: str | None = "isometric", + notebook: bool = True, + theme: str = "default", + **kwargs: Any, +) -> Any | None: + """Create interactive 3D visualisation of prisms using PyVista. + + Args: + geometry_model: A GeometryModel with prisms and bbox. + show_edges: Whether to show edges of the prisms. + opacity: Base opacity (0.0-1.0). Core layer forced opaque. + color_by_layer: Colour by layer name. + show_simulation_box: Draw the simulation bounding box. + camera_position: "isometric", "xy", "xz", "yz", or custom tuple. + notebook: Whether running inside Jupyter. + theme: PyVista theme ("default", "dark", "document"). + **kwargs: Extra args forwarded to ``pv.Plotter``. + + Returns: + PyVista plotter object for further customisation. + """ + if not PYVISTA_AVAILABLE: + raise ImportError( + "PyVista is required for 3D visualization. " + "Install with: pip install pyvista" + ) + + plotter = ( + pv.Plotter(notebook=notebook, **kwargs) if notebook else pv.Plotter(**kwargs) + ) + + if theme == "dark": + pv.set_plot_theme("dark") + elif theme == "document": + pv.set_plot_theme("document") + + layer_meshes = _convert_prisms_to_meshes(geometry_model) + colors = generate_layer_colors(list(layer_meshes.keys())) + + for layer_name, meshes in layer_meshes.items(): + color = colors[layer_name] if color_by_layer else None + layer_opacity = 1.0 if layer_name == "core" else 0.2 + + for mesh in meshes: + plotter.add_mesh( + mesh, + color=color, + opacity=layer_opacity, + show_edges=show_edges, + label=layer_name, + name=f"{layer_name}_{id(mesh)}", + ) + + if show_simulation_box: + sim_box = _create_simulation_box_pv(geometry_model) + plotter.add_mesh( + sim_box, + style="wireframe", + color="black", + line_width=2, + label="Simulation Box", + ) + + _set_camera_position(plotter, camera_position) + + if color_by_layer and len(layer_meshes) > 1: + try: # noqa: SIM105 + plotter.add_legend() + except Exception: + pass + + if notebook: + try: + pv.set_jupyter_backend("trame") + return plotter.show() + except Exception: + pv.set_jupyter_backend("static") + return plotter.show() + else: + return plotter.show() + + +def export_3d_mesh( + geometry_model: GeometryModel, + filename: str, + fmt: str = "auto", +) -> None: + """Export 3D geometry to mesh file (STL, PLY, OBJ, VTK, glTF).""" + if not PYVISTA_AVAILABLE: + raise ImportError( + "PyVista is required for mesh export. Install with: pip install pyvista" + ) + + layer_meshes = _convert_prisms_to_meshes(geometry_model) + combined = pv.MultiBlock() + for layer_name, meshes in layer_meshes.items(): + for i, mesh in enumerate(meshes): + combined[f"{layer_name}_{i}"] = mesh + + if fmt == "auto": + fmt = filename.rsplit(".", 1)[-1].lower() + + if fmt == "stl": + combined.combine().save(filename) + elif fmt in ("ply", "obj", "vtk"): + combined.save(filename) + elif fmt == "gltf": + try: + combined.save(filename) + except Exception as e: + raise ValueError( + f"glTF export failed. May need additional dependencies: {e}" + ) from e + else: + raise ValueError(f"Unsupported format: {fmt}") + + +def create_web_export( + geometry_model: GeometryModel, + filename: str = "geometry_3d.html", + title: str = "3D Geometry Visualization", # noqa: ARG001 +) -> str: + """Export 3D visualisation as standalone HTML (via PyVista).""" + if not PYVISTA_AVAILABLE: + raise ImportError("PyVista is required for web export.") + + plotter = pv.Plotter(notebook=False, off_screen=True) + layer_meshes = _convert_prisms_to_meshes(geometry_model) + colors = generate_layer_colors(list(layer_meshes.keys())) + + for layer_name, meshes in layer_meshes.items(): + color = colors[layer_name] + for mesh in meshes: + plotter.add_mesh(mesh, color=color, opacity=0.8) + + plotter.export_html(filename) + return filename + + +# --------------------------------------------------------------------------- +# Internal helpers +# --------------------------------------------------------------------------- + + +def _convert_prisms_to_meshes( + geometry_model: GeometryModel, +) -> dict[str, list[Any]]: + """Convert generic Prisms to PyVista meshes, with triangle-merge optimisation.""" + layer_meshes: dict[str, list[Any]] = {} + + for layer_name, prisms in geometry_model.prisms.items(): + triangular_count = sum(1 for p in prisms if len(p.vertices) == 3) + + if triangular_count > 100: + triangular = [p for p in prisms if len(p.vertices) == 3] + non_triangular = [p for p in prisms if len(p.vertices) != 3] + + meshes: list[Any] = [] + if triangular: + merged = _merge_triangular_prisms(triangular) + if merged is not None: + meshes.append(merged) + + for prism in non_triangular: + base, top = prism_base_top_vertices(prism) + meshes.append(_create_prism_mesh(base, top, prism)) + + layer_meshes[layer_name] = meshes + else: + meshes = [] + for prism in prisms: + base, top = prism_base_top_vertices(prism) + meshes.append(_create_prism_mesh(base, top, prism)) + layer_meshes[layer_name] = meshes + + return layer_meshes + + +def _merge_triangular_prisms(prisms: list[Prism]) -> Any: + """Merge many triangular prisms into one PyVista mesh for performance.""" + all_vertices: list[np.ndarray] = [] + all_faces: list[int] = [] + offset = 0 + + for prism in prisms: + base, top = prism_base_top_vertices(prism) + verts = np.vstack([base, top]) + all_vertices.append(verts) + n = len(base) + + # bottom face + all_faces.extend([3, offset + 2, offset + 1, offset + 0]) + # top face + all_faces.extend([3, offset + n, offset + n + 1, offset + n + 2]) + # side quads + for i in range(n): + ni = (i + 1) % n + all_faces.extend( + [ + 4, + offset + i, + offset + ni, + offset + ni + n, + offset + i + n, + ] + ) + offset += len(verts) + + if all_vertices: + combined = np.vstack(all_vertices) + return pv.PolyData(combined, all_faces) + return None + + +def _create_prism_mesh( + base_vertices: np.ndarray, + top_vertices: np.ndarray, + prism: Prism | None = None, +) -> Any: + """Create a PyVista mesh from base/top 3D vertices, with hole/concave support.""" + if prism is not None and prism.original_polygon is not None: + poly = prism.original_polygon + has_holes = hasattr(poly, "interiors") and poly.interiors + is_concave = ( + not has_holes + and hasattr(poly, "convex_hull") + and poly.convex_hull.area > poly.area * 1.001 + ) + if has_holes or is_concave: + return _create_prism_mesh_with_holes(poly, base_vertices, top_vertices) + + n = len(base_vertices) + all_verts = np.vstack([base_vertices, top_vertices]) + + faces: list[int] = [] + # bottom + faces.extend([n, *list(range(n))[::-1]]) + # top + faces.extend([n, *[i + n for i in range(n)]]) + # sides + for i in range(n): + ni = (i + 1) % n + faces.extend([4, i, ni, ni + n, i + n]) + + return pv.PolyData(all_verts, faces) + + +def _create_prism_mesh_with_holes( + shapely_polygon: Any, + base_vertices: np.ndarray, + top_vertices: np.ndarray, +) -> Any: + """Create PyVista mesh from a Shapely polygon with holes using Delaunay.""" + try: + import shapely.geometry as sg + from scipy.spatial import Delaunay + except ImportError: + return _create_prism_mesh(base_vertices, top_vertices) + + all_points: list[tuple[float, float]] = [] + boundary_segments: list[list[int]] = [] + + exterior_coords = list(shapely_polygon.exterior.coords[:-1]) + start_idx = 0 + all_points.extend(exterior_coords) + boundary_segments = [ + [start_idx + i, start_idx + (i + 1) % len(exterior_coords)] + for i in range(len(exterior_coords)) + ] + + for interior in shapely_polygon.interiors: + interior_coords = list(interior.coords[:-1]) + start_idx = len(all_points) + all_points.extend(interior_coords) + boundary_segments.extend( + [start_idx + i, start_idx + (i + 1) % len(interior_coords)] + for i in range(len(interior_coords)) + ) + + if len(all_points) < 3: + return _create_prism_mesh(base_vertices, top_vertices) + + points_2d = np.array(all_points) + tri = Delaunay(points_2d) + + valid_triangles = [ + simplex + for simplex in tri.simplices + if shapely_polygon.contains(sg.Point(*np.mean(points_2d[simplex], axis=0))) + ] + + if not valid_triangles: + return _create_prism_mesh(base_vertices, top_vertices) + + z_base = base_vertices[0, 2] if len(base_vertices) > 0 else 0 + z_top = top_vertices[0, 2] if len(top_vertices) > 0 else z_base + 1 + + verts_3d: list[list[float]] = [[pt[0], pt[1], z_base] for pt in points_2d] + [ + [pt[0], pt[1], z_top] for pt in points_2d + ] + + n_pts = len(points_2d) + faces_pv: list[int] = [] + + for tri_idx in valid_triangles: + faces_pv.extend([3, tri_idx[0], tri_idx[1], tri_idx[2]]) + for tri_idx in valid_triangles: + faces_pv.extend( + [ + 3, + tri_idx[0] + n_pts, + tri_idx[2] + n_pts, + tri_idx[1] + n_pts, + ] + ) + for seg in boundary_segments: + i, j = seg + faces_pv.extend([4, i, j, j + n_pts, i + n_pts]) + + try: + return pv.PolyData(verts_3d, faces_pv) + except Exception: + return _create_prism_mesh(base_vertices, top_vertices) + + +def _create_simulation_box_pv(geometry_model: GeometryModel) -> Any: + """Create a wireframe PyVista box for the simulation bbox.""" + mn = geometry_model.bbox[0] + mx = geometry_model.bbox[1] + bounds = [mn[0], mx[0], mn[1], mx[1], mn[2], mx[2]] + return pv.Box(bounds=bounds) + + +def _set_camera_position(plotter: Any, position: str | None) -> None: + """Configure camera position for optimal viewing.""" + if position == "isometric": + plotter.camera_position = "iso" + elif position == "xy": + plotter.view_xy() + elif position == "xz": + plotter.view_xz() + elif position == "yz": + plotter.view_yz() + elif isinstance(position, (tuple, list)) and len(position) == 3: + plotter.camera_position = position + else: + plotter.camera_position = "iso" + plotter.reset_camera() diff --git a/src/gsim/common/viz/render3d_threejs.py b/src/gsim/common/viz/render3d_threejs.py new file mode 100644 index 00000000..d9daf3c0 --- /dev/null +++ b/src/gsim/common/viz/render3d_threejs.py @@ -0,0 +1,233 @@ +"""Three.js / FastAPI-based 3D rendering for GeometryModel. + +Provides browser-based interactive 3D visualisation via a local FastAPI server +that serves a Three.js HTML page. +""" + +from __future__ import annotations + +import json +import logging +from pathlib import Path +from typing import Any + +import numpy as np + +from gsim.common.geometry_model import GeometryModel +from gsim.common.viz._colors import generate_layer_colors_with_opacity +from gsim.common.viz._mesh_helpers import simulation_box_corners +from gsim.common.viz.render3d_open3d import _convert_prisms_to_open3d + +logger = logging.getLogger(__name__) + +# --------------------------------------------------------------------------- +# Public API +# --------------------------------------------------------------------------- + + +def serve_threejs_visualization( + geometry_model: GeometryModel, + *, + show_edges: bool = False, + color_by_layer: bool = True, + show_simulation_box: bool = True, + layer_opacity: dict[str, float] | None = None, + port: int = 8000, + auto_open: bool = True, + show_stats: bool = False, + **kwargs: Any, +) -> str: + """Start a FastAPI server with a Three.js 3D viewer. + + Args: + geometry_model: GeometryModel containing prisms and bbox. + show_edges: Show wireframe edges. + color_by_layer: Colour each layer differently. + show_simulation_box: Draw the simulation box. + layer_opacity: Per-layer opacity override. + port: Port to serve on (default 8000). + auto_open: Open the browser automatically. + show_stats: Show FPS counter. + **kwargs: Extra Three.js options. + + Returns: + URL of the running server. + """ + try: + import threading + import time + import webbrowser + + import uvicorn + from fastapi import FastAPI + from fastapi.responses import HTMLResponse + except ImportError as err: + raise ImportError( + "FastAPI and uvicorn required. Install with: pip install fastapi uvicorn" + ) from err + + app = FastAPI(title="3D Geometry Viewer") + + layer_meshes = _convert_prisms_to_open3d(geometry_model) + colors, opacity_dict = generate_layer_colors_with_opacity( + list(layer_meshes.keys()), layer_opacity + ) + + logger.info("Converting geometry: %d layers found", len(layer_meshes)) + for layer_name, meshes in layer_meshes.items(): + logger.info(" Layer '%s': %d meshes", layer_name, len(meshes)) + + threejs_data = _convert_to_threejs_data( + layer_meshes, colors, opacity_dict, color_by_layer + ) + + if show_simulation_box: + threejs_data["simulation_box"] = _create_simulation_box_threejs(geometry_model) + + total_vertices = 0 + total_faces = 0 + for layer in threejs_data.get("layers", []): + for mesh in layer.get("meshes", []): + total_vertices += len(mesh.get("vertices", [])) // 3 + total_faces += len(mesh.get("faces", [])) // 3 + logger.info( + "Three.js data prepared: %d vertices, %d faces", + total_vertices, + total_faces, + ) + + @app.get("/", response_class=HTMLResponse) + def get_visualization(): + """Return the Three.js viewer HTML page.""" + return _generate_threejs_html( + threejs_data, + show_edges=show_edges, + show_stats=show_stats, + **kwargs, + ) + + @app.get("/api/geometry") + def get_geometry_data(): + """Return geometry mesh data as JSON.""" + return threejs_data + + @app.get("/api/info") + def get_info(): + """Return summary info about the geometry.""" + return { + "layers": len(threejs_data.get("layers", [])), + "total_meshes": sum( + len(layer.get("meshes", [])) for layer in threejs_data.get("layers", []) + ), + "has_simulation_box": "simulation_box" in threejs_data, + "layer_names": [ + layer.get("name") for layer in threejs_data.get("layers", []) + ], + } + + server_url = f"http://localhost:{port}" + + def run_server(): + try: + uvicorn.run(app, host="127.0.0.1", port=port, log_level="info") + except Exception as e: + logger.info("Server error: %s", e) + + server_thread = threading.Thread(target=run_server, daemon=True) + server_thread.start() + time.sleep(1) + + if auto_open: + webbrowser.open(server_url) + + logger.info("Server started at: %s", server_url) + logger.info("Geometry API available at: %s/api/geometry", server_url) + logger.info("Server running in background. Keep Python session alive to view.") + + return server_url + + +# --------------------------------------------------------------------------- +# Internal helpers +# --------------------------------------------------------------------------- + + +def _convert_to_threejs_data( + layer_meshes: dict[str, list[Any]], + colors: dict[str, list[float]], + opacity_dict: dict[str, float], + color_by_layer: bool, +) -> dict[str, Any]: + """Convert Open3D meshes to Three.js-compatible JSON data.""" + threejs_data: dict[str, Any] = {"layers": []} + + for layer_name, meshes in layer_meshes.items(): + layer_color = colors[layer_name] if color_by_layer else [0.7, 0.7, 0.7] + layer_opacity = opacity_dict.get(layer_name, 0.8) + + layer_data: dict[str, Any] = { + "name": layer_name, + "color": [int(c * 255) for c in layer_color[:3]], + "opacity": layer_opacity, + "meshes": [], + } + + for i, mesh in enumerate(meshes): + vertices = np.asarray(mesh.vertices).flatten().tolist() + faces = np.asarray(mesh.triangles).flatten().tolist() + + layer_data["meshes"].append( + { + "vertices": vertices, + "faces": faces, + "id": f"{layer_name}_{i}", + } + ) + + threejs_data["layers"].append(layer_data) + + return threejs_data + + +def _create_simulation_box_threejs( + geometry_model: GeometryModel, +) -> dict[str, Any]: + """Create simulation box data for Three.js.""" + points, lines = simulation_box_corners(geometry_model) + return { + "vertices": points.flatten().tolist(), + "lines": [idx for pair in lines for idx in pair], + "color": [0, 0, 0], + } + + +def _generate_threejs_html( + threejs_data: dict[str, Any], + *, + show_edges: bool = False, + show_stats: bool = False, + **kwargs: Any, +) -> str: + """Generate HTML using the viewer.html template.""" + template_path = Path(__file__).parent / "templates" / "viewer.html" + + with open(template_path) as f: + template = f.read() + + width = kwargs.get("width", "100vw") + height = kwargs.get("height", "100vh") + background_color = kwargs.get("background_color", "#f0f0f0") + show_wireframe = str(show_edges).lower() + stats_display = "block" if show_stats else "none" + + data_json = json.dumps(threejs_data, indent=2) + + return template.format( + geometry_data=data_json, + show_wireframe=show_wireframe, + show_stats=str(show_stats).lower(), + width=width, + height=height, + background_color=background_color, + stats_display=stats_display, + ) diff --git a/src/gsim/common/viz/templates/viewer.html b/src/gsim/common/viz/templates/viewer.html new file mode 100644 index 00000000..23323878 --- /dev/null +++ b/src/gsim/common/viz/templates/viewer.html @@ -0,0 +1,482 @@ + + + + + + 3D Geometry - Three.js Viewer + + + +
+
+

🔬 3D Geometry Visualization

+

Mouse: Rotate view

+

Wheel: Zoom in/out (enhanced sensitivity)

+

Right-click: Pan camera

+

Double-click: Reset view

+
+ +
+ +
+ +
+ +
+ +
+

📊 Layers

+
+
+ +
+
FPS: --
+
Vertices: --
+
Faces: --
+
+
+ + + + + + + + + diff --git a/src/gsim/gcloud.py b/src/gsim/gcloud.py index 95e70dff..2b2cb508 100644 --- a/src/gsim/gcloud.py +++ b/src/gsim/gcloud.py @@ -110,7 +110,7 @@ def upload_simulation_dir(input_dir: str | Path, job_type: str): def run_simulation( output_dir: str | Path, - job_type: Literal["palace"] = "palace", + job_type: Literal["palace", "meep"] = "palace", verbose: bool = True, on_started: Callable | None = None, ) -> dict[str, Path]: diff --git a/src/gsim/meep/__init__.py b/src/gsim/meep/__init__.py new file mode 100644 index 00000000..dec415b4 --- /dev/null +++ b/src/gsim/meep/__init__.py @@ -0,0 +1,52 @@ +"""MEEP photonic FDTD simulation module. + +Provides a declarative API (``Simulation``) for configuring and running +MEEP FDTD simulations on the GDSFactory+ cloud. No local MEEP +installation required. + +Example:: + + from gsim import meep + + sim = meep.Simulation() + sim.geometry.component = ybranch + sim.source.port = "o1" + sim.monitors = ["o1", "o2"] + sim.solver.stopping = "dft_decay" + sim.solver.max_time = 200 + result = sim.run("./meep-sim") +""" + +from gsim.meep.models import ( + FDTD, + Domain, + DomainConfig, + Geometry, + Material, + ModeSource, + ResolutionConfig, + SimConfig, + SourceConfig, + SParameterResult, + Symmetry, + # Config models used by Simulation.write_config() + WavelengthConfig, +) +from gsim.meep.simulation import BuildResult, Simulation + +__all__ = [ + "FDTD", + "BuildResult", + "Domain", + "DomainConfig", + "Geometry", + "Material", + "ModeSource", + "ResolutionConfig", + "SParameterResult", + "SimConfig", + "Simulation", + "SourceConfig", + "Symmetry", + "WavelengthConfig", +] diff --git a/src/gsim/meep/materials.py b/src/gsim/meep/materials.py new file mode 100644 index 00000000..25fa9c4b --- /dev/null +++ b/src/gsim/meep/materials.py @@ -0,0 +1,68 @@ +"""Material resolution for MEEP simulation. + +Resolves material names from the common materials database to optical +properties needed for photonic FDTD simulation. No MEEP imports. +""" + +from __future__ import annotations + +import warnings + +from gsim.common.stack.materials import MaterialProperties, get_material_properties +from gsim.meep.models.config import MaterialData + + +def resolve_materials( + used_material_names: set[str], + overrides: dict[str, MaterialProperties] | None = None, +) -> dict[str, MaterialData]: + """Resolve material names to optical properties for MEEP. + + Only resolves materials that actually have geometry (i.e., polygons + extracted from the component). Materials present in the GDS but + missing optical data emit a warning and are excluded from the config. + + Args: + used_material_names: Material names that appear in extracted geometry + overrides: User-supplied material property overrides (from set_material()) + + Returns: + Dict mapping material name to MaterialData + """ + overrides = overrides or {} + materials: dict[str, MaterialData] = {} + + for name in sorted(used_material_names): + # Check user overrides first + if name in overrides: + props = overrides[name] + if props.refractive_index is None: + warnings.warn( + f"Material override '{name}' has no refractive_index — " + f"skipping. Use set_material('{name}', refractive_index=...)", + stacklevel=2, + ) + continue + materials[name] = MaterialData( + refractive_index=props.refractive_index, + extinction_coeff=props.extinction_coeff or 0.0, + ) + continue + + # Look up in common DB + db_props = get_material_properties(name) + if db_props is not None and db_props.refractive_index is not None: + materials[name] = MaterialData( + refractive_index=db_props.refractive_index, + extinction_coeff=db_props.extinction_coeff or 0.0, + ) + continue + + warnings.warn( + f"Material '{name}' has no optical properties (refractive_index) " + f"— layer will be omitted from simulation. " + f"Use sim.set_material('{name}', refractive_index=...) to include it.", + stacklevel=2, + ) + + return materials diff --git a/src/gsim/meep/meep-overview.md b/src/gsim/meep/meep-overview.md new file mode 100644 index 00000000..b15cbabd --- /dev/null +++ b/src/gsim/meep/meep-overview.md @@ -0,0 +1,135 @@ +# gsim.meep — Overview & Roadmap + +## Architecture + +``` +GDS + PDK ──► 3D Geometry ──► Client Viz ──► SimConfig JSON ──► Cloud Runner +(common/) (GeometryModel) (no meep) (no meep) (meep in Docker) +``` + +**Key constraint:** gsim is the client SDK — no meep dependency. MEEP runs only in Docker/cloud. + +**Upload package:** `layout.gds` + `sim_config.json` + `run_meep.py` -> Docker -> `s_parameters.csv` + `meep_debug.json` + diagnostic PNGs. + +--- + +## API: Declarative `Simulation` + +```python +from gsim import meep + +sim = meep.Simulation() +sim.geometry.component = ybranch +sim.geometry.z_crop = "auto" +sim.materials = {"si": 3.47, "SiO2": 1.44} # float shorthand +sim.source.port = "o1" +sim.source.wavelength = 1.55 +sim.source.bandwidth = 0.1 +sim.source.num_freqs = 11 +sim.monitors = ["o1", "o2"] +sim.domain.pml = 1.0 +sim.domain.margin = 0.5 +sim.solver.resolution = 32 +sim.solver.stopping = "dft_decay" +sim.solver.max_time = 200 +sim.solver.simplify_tol = 0.01 +sim.solver.save_geometry = True +sim.solver.save_fields = True +sim.plot_2d(slices="xyz") +result = sim.run("./meep-sim") +``` + +### Design principles + +- **6 typed physics objects** -- `Geometry`, `Material`, `ModeSource`, `Domain`, `FDTD` + `monitors: list[str]` -- assigned to a `Simulation` container. No ordering dependencies. +- **Attribute style or constructor style** -- `sim.source.port = "o1"` or `sim.source = ModeSource(port="o1")`. Both work via Pydantic `validate_assignment=True`. +- **Float shorthand for materials** -- `{"si": 3.47}` auto-normalizes to `Material(n=3.47)` via validator. +- **Flat stopping fields** -- `stopping` mode string + flat fields (`max_time`, `stopping_threshold`, etc.) on `FDTD`. +- **Source defines spectral window** -- `WavelengthConfig` derived from `ModeSource.wavelength/bandwidth/num_freqs`. Monitors are just port name strings. +- **JSON contract unchanged** -- `write_config()` translates new API -> existing `SimConfig` -> JSON. Runner template untouched. + +--- + +## Module Structure + +### `gsim.meep` -- Public API + +| Module | Purpose | +| ------------------- | ----------------------------------------------------------------------------------------------------------------------- | +| `__init__.py` | Exports: `Simulation`, all model classes | +| `models/api.py` | Declarative models: `Geometry`, `Material`, `ModeSource`, `Domain`, `FDTD`, `Symmetry` | +| `simulation.py` | `Simulation` container -- `write_config()`, `run()`, `validate_config()`, `plot_2d()`/`plot_3d()`, `estimate_meep_np()` | +| `viz.py` | Standalone viz helpers -- `build_geometry_model()`, `plot_2d()`, `plot_3d()` (calls `common/viz`) | +| `models/config.py` | `SimConfig` + all sub-configs (JSON serialization layer) | +| `models/results.py` | `SParameterResult` -- CSV + debug JSON + diagnostic PNGs | +| `ports.py` | `extract_port_info()` -- port center/direction/normal from gdsfactory | +| `materials.py` | `resolve_materials()` -- material names -> (n, k) via common DB | +| `script.py` | `generate_meep_script()` -- cloud runner template (string in Python) | +| `overlay.py` | `SimOverlay` + `PortOverlay` + `DielectricOverlay` -- viz metadata | + +### `gsim.common` -- Shared infrastructure + +Solver-agnostic: `LayeredComponentBase`, `GeometryModel`/`Prism`, `LayerStack`, `MaterialProperties`, `viz/` (Matplotlib 2D, PyVista/Open3D/Three.js 3D). + +### Visualization pipeline + +``` +Simulation.plot_2d() + -> viz.build_geometry_model(component, stack, domain_config) + -> LayeredComponentBase + extract_geometry_model() + -> crop_geometry_model() if stack is z-cropped + -> viz.build_overlay(gm, component, stack, domain_config) + -> overlay.build_sim_overlay() + -> common/viz/render2d.plot_prism_slices(gm, overlay=overlay) +``` + +- **Z slices** -- draw actual prism polygon outlines at the slice plane (XY view) +- **X/Y slices** -- compute Shapely line-polygon intersections for accurate cross-sections (XZ/YZ views); each intersection segment becomes a rectangle from z_base to z_top +- **Overlay** -- sim cell boundary, PML regions, port markers, dielectric backgrounds + +--- + +## Key Design Decisions + +- **GDS-file approach** -- Send raw GDS to cloud (not polygon coords in JSON). Complex geometries have 1000s of nodes; DerivedLayers need gdsfactory to resolve. +- **Port z-center = highest refractive index** -- Photonic ports center on waveguide core (not conductor layer like RF). +- **Z-crop** -- Auto-crops stack around core layer. Full UBC stack is 15um; cropped to ~1.5um (massive compute savings). +- **Port extension into PML** -- `gf.components.extend_ports()` at `write_config()` time. Original bbox stored in `SimConfig.component_bbox` for correct cell sizing. +- **Symmetries disabled for S-params** -- `add_mode_monitor` uses `use_symmetry=false` internally; `get_eigenmode_coefficients` doesn't apply `S.multiplicity()`. Source port coefficients underestimated ~2x. gplugins also never uses `mp.Mirror`. +- **`dft_min_run_time` default 100** -- Prevents false convergence before pulse traverses device. Must exceed `device_length * n_group`. +- **Stack resolved lazily** -- `_ensure_stack()` falls back to `get_stack()` (active PDK defaults) when no explicit stack kwargs are set. Same path for `write_config()` and viz. + +--- + +## TODO + +### Near-term + +- [ ] **Cloud end-to-end test** -- Full round-trip (upload -> run meep -> download results) hasn't been tested with real cloud instance. +- [ ] **Monitor z-span touches PML** -- After z-crop, layer stack z-extent exactly matches PML inner boundary. At coarse resolution, outermost monitor pixels may straddle PML. Fix: shrink monitor z-span by ~`2/resolution` inset. +- [ ] **Empty `core2` layer in config** -- UBC PDK's ebeam_y_1550 includes a `core2` entry (GDS layer 31) with zero geometry. Harmless but clutters config. Consider filtering empty layers. + +### Medium-term + +- [ ] **Custom monitors** -- Monitors are port-name strings only. Add `FieldMonitor(center, size)` and `FluxMonitor` for custom measurement locations. +- [ ] **Per-port margin/mode_index** -- `port_margin` is global in `Domain`. Could add per-port control for margin and higher-order modes. +- [ ] **Typed boundaries** -- Only PML today. Add `Periodic`, `Bloch`, `PEC`, `PMC` per-axis boundary specs. +- [ ] **Mesh refinement regions** -- Single global `resolution`. Add local refinement for thin features. +- [ ] **Dispersive materials** -- `Material(n, k)` is non-dispersive. Add Sellmeier/Lorentz/Drude support. +- [ ] **`updated_copy()` for parameter sweeps** -- Pydantic `model_copy(update={...})` makes this nearly free. + +### Deferred + +- [ ] **JSON schema v2** -- Restructure JSON groups (`"solver.wavelength"`, per-port margins, etc.) with version field for runner backward compat. Requires atomic client + runner update. +- [ ] **Port symmetries** -- gplugins-style: run fewer source ports, copy S-params between symmetric port pairs (e.g. S31=S21 for Y-branch). +- [ ] **`add_flux` + `eig_parity` path** -- Alternative to `add_mode_monitor` that works correctly with `mp.Mirror` symmetry. Needs auto-detection of correct parity from symmetry config and waveguide polarization. + +--- + +## Docker / Cloud + +- **Base:** `continuumio/miniconda3` -> conda env with `pymeep=*=mpi_mpich_*`, `gdsfactory`, `nlopt` +- **Entrypoint:** downloads input -> `mpirun -np $NP python run_meep.py` -> uploads outputs +- **`meep_np` hardcoded to 2** -- auto-compute logic exists (`total_voxels // 200k`) but `lscpu` reports host cores inside containers, not Batch-allocated vCPUs. TODO: pass Batch vCPU allocation as `MEEP_NP` env var from job submission, then restore auto-compute with proper clamping +- **Instance recommendation:** `c7a.16xlarge` (AMD EPYC Genoa, DDR5, ~460 GB/s mem BW) for typical photonics; `hpc7a.96xlarge` for large problems +- **Local testing:** `./nbs/test_meep_local.sh` (regenerate -> Docker build -> run -> results) diff --git a/src/gsim/meep/models/__init__.py b/src/gsim/meep/models/__init__.py new file mode 100644 index 00000000..d7cd302d --- /dev/null +++ b/src/gsim/meep/models/__init__.py @@ -0,0 +1,51 @@ +"""MEEP simulation models.""" + +from gsim.meep.models.api import ( + FDTD, + Domain, + Geometry, + Material, + ModeSource, + Symmetry, +) +from gsim.meep.models.config import ( + AccuracyConfig, + DiagnosticsConfig, + DomainConfig, + LayerStackEntry, + MaterialData, + PortData, + ResolutionConfig, + SimConfig, + SourceConfig, + StoppingConfig, + SymmetryEntry, + WavelengthConfig, +) +from gsim.meep.models.results import SParameterResult + +# Backward compatibility alias +FDTDConfig = WavelengthConfig + +__all__ = [ + "FDTD", + "AccuracyConfig", + "DiagnosticsConfig", + "Domain", + "DomainConfig", + "FDTDConfig", + "Geometry", + "LayerStackEntry", + "Material", + "MaterialData", + "ModeSource", + "PortData", + "ResolutionConfig", + "SParameterResult", + "SimConfig", + "SourceConfig", + "StoppingConfig", + "Symmetry", + "SymmetryEntry", + "WavelengthConfig", +] diff --git a/src/gsim/meep/models/api.py b/src/gsim/meep/models/api.py new file mode 100644 index 00000000..2cee14cc --- /dev/null +++ b/src/gsim/meep/models/api.py @@ -0,0 +1,183 @@ +"""Declarative API models for MEEP photonic simulation. + +Typed physics objects that map to the underlying SimConfig JSON contract. +These models are the user-facing API; the translation to SimConfig happens +in ``Simulation.write_config()``. +""" + +from __future__ import annotations + +from typing import Any, Literal + +from pydantic import BaseModel, ConfigDict, Field + +# --------------------------------------------------------------------------- +# Geometry +# --------------------------------------------------------------------------- + + +class Geometry(BaseModel): + """Physical layout: component + layer stack + optional z-crop.""" + + model_config = ConfigDict( + validate_assignment=True, + arbitrary_types_allowed=True, + ) + + component: Any = None + stack: Any = None + z_crop: str | None = Field( + default=None, + description='Z-crop mode: "auto" | layer_name | None (no crop)', + ) + + +# --------------------------------------------------------------------------- +# Material +# --------------------------------------------------------------------------- + + +class Material(BaseModel): + """Optical material properties.""" + + model_config = ConfigDict(validate_assignment=True) + + n: float = Field(gt=0, description="Refractive index") + k: float = Field(default=0.0, ge=0, description="Extinction coefficient") + + +# --------------------------------------------------------------------------- +# Source +# --------------------------------------------------------------------------- + + +class ModeSource(BaseModel): + """Mode source excitation and spectral measurement window.""" + + model_config = ConfigDict(validate_assignment=True) + + port: str | None = Field( + default=None, + description="Source port name. None = auto-select first port.", + ) + wavelength: float = Field( + default=1.55, + gt=0, + description="Center wavelength in um", + ) + bandwidth: float = Field( + default=0.1, + ge=0, + description="Measurement wavelength bandwidth in um", + ) + num_freqs: int = Field( + default=11, + ge=1, + description="Number of frequency points", + ) + + +# --------------------------------------------------------------------------- +# Domain +# --------------------------------------------------------------------------- + + +class Symmetry(BaseModel): + """Mirror symmetry plane.""" + + model_config = ConfigDict(validate_assignment=True) + + direction: Literal["X", "Y", "Z"] + phase: Literal[1, -1] = Field(default=1) + + +class Domain(BaseModel): + """Computational domain sizing: PML + margins + symmetries.""" + + model_config = ConfigDict(validate_assignment=True) + + pml: float = Field(default=1.0, ge=0, description="PML thickness in um") + margin: float = Field( + default=0.5, + ge=0, + description="XY margin between geometry and PML in um", + ) + margin_z_above: float = Field( + default=0.5, ge=0, description="Z margin above core in um" + ) + margin_z_below: float = Field( + default=0.5, ge=0, description="Z margin below core in um" + ) + port_margin: float = Field( + default=0.5, + ge=0, + description="Margin on each side of port width for monitors (um)", + ) + extend_ports: float = Field( + default=0.0, + ge=0, + description="Extend ports into PML (um). 0 = auto (margin + pml).", + ) + symmetries: list[Symmetry] = Field( + default_factory=list, + description="Mirror symmetry planes. Not yet used in production runs.", + ) + + +# --------------------------------------------------------------------------- +# FDTD solver +# --------------------------------------------------------------------------- + + +class FDTD(BaseModel): + """Solver numerics: resolution, stopping, subpixel, diagnostics.""" + + model_config = ConfigDict(validate_assignment=True) + + resolution: int = Field(default=32, ge=4, description="Pixels per micrometer") + + # Stopping criteria (flat fields instead of variant classes) + stopping: Literal["fixed", "decay", "dft_decay"] = Field( + default="dft_decay", + description="Stopping mode: fixed time, field decay, or DFT convergence", + ) + max_time: float = Field( + default=200.0, gt=0, description="Max run time after sources (um/c)" + ) + stopping_threshold: float = Field( + default=1e-3, gt=0, lt=1, description="Decay/convergence threshold" + ) + stopping_min_time: float = Field( + default=100.0, ge=0, description="Min run time for dft_decay mode" + ) + stopping_component: str = Field( + default="Ey", description="Field component for decay mode" + ) + stopping_dt: float = Field( + default=50.0, gt=0, description="Decay measurement window for decay mode" + ) + stopping_monitor_port: str | None = Field( + default=None, description="Port to monitor for decay mode" + ) + + subpixel: bool = Field(default=False, description="Toggle subpixel averaging") + subpixel_maxeval: int = Field( + default=0, ge=0, description="Cap on integration evaluations (0=unlimited)" + ) + subpixel_tol: float = Field( + default=1e-4, gt=0, description="Subpixel integration tolerance" + ) + simplify_tol: float = Field( + default=0.0, + ge=0, + description="Shapely simplification tolerance in um (0=off)", + ) + + # Diagnostics — output control for plots, fields, animations + save_geometry: bool = Field(default=True) + save_fields: bool = Field(default=True) + save_epsilon_raw: bool = Field(default=False) + save_animation: bool = Field(default=False) + animation_interval: float = Field(default=0.5, gt=0) + 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 new file mode 100644 index 00000000..52abcfa7 --- /dev/null +++ b/src/gsim/meep/models/config.py @@ -0,0 +1,341 @@ +"""Configuration models for MEEP photonic simulation. + +Defines serializable Pydantic models for the complete simulation config +that gets written as JSON and consumed by the cloud MEEP runner script. +""" + +from __future__ import annotations + +import json +from pathlib import Path +from typing import Any, Literal + +from pydantic import BaseModel, ConfigDict, Field, computed_field + + +class SymmetryEntry(BaseModel): + """One MEEP mirror symmetry plane.""" + + model_config = ConfigDict(validate_assignment=True) + + direction: Literal["X", "Y", "Z"] + phase: Literal[1, -1] = Field(default=1) + + +class DomainConfig(BaseModel): + """Simulation domain sizing: margins around geometry + PML thickness. + + Margins control how much material (from the layer stack) is kept around + the waveguide core. ``set_z_crop()`` uses ``margin_z_above`` / + ``margin_z_below`` to determine the crop window. Along XY the margin + is the gap between the geometry bounding-box and the PML inner edge. + + Cell size formula: + cell_x = bbox_width + 2*(margin_xy + dpml) + cell_y = bbox_height + 2*(margin_xy + dpml) + cell_z = z_extent + 2*dpml (z-margins baked into z_extent via set_z_crop) + """ + + model_config = ConfigDict(validate_assignment=True) + + dpml: float = Field(default=1.0, ge=0, description="PML thickness in um") + margin_xy: float = Field( + default=0.5, ge=0, description="XY margin between geometry and PML in um" + ) + margin_z_above: float = Field( + default=0.5, ge=0, description="Z margin above core kept by set_z_crop in um" + ) + margin_z_below: float = Field( + default=0.5, ge=0, description="Z margin below core kept by set_z_crop in um" + ) + port_margin: float = Field( + default=0.5, + ge=0, + description="Margin on each side of port waveguide width for mode monitors in um", + ) + extend_ports: float = Field( + default=0.0, + ge=0, + description="Length to extend waveguide ports into PML in um. " + "0 = auto (margin_xy + dpml).", + ) + + +class StoppingConfig(BaseModel): + """Controls when the MEEP simulation stops. + + ``fixed`` mode runs for a fixed time after sources turn off. + ``decay`` mode monitors field decay at a point and stops when the + fields have decayed by ``threshold``, with ``max_time`` as + a numeric time cap (whichever fires first). + ``dft_decay`` mode monitors convergence of all DFT monitors and + stops when they stabilize, with built-in min/max time bounds. + Best for S-parameter extraction. + """ + + model_config = ConfigDict(validate_assignment=True) + + mode: Literal["fixed", "decay", "dft_decay"] = Field(default="fixed") + max_time: float = Field(default=100.0, gt=0, serialization_alias="run_after_sources") + decay_dt: float = Field(default=50.0, gt=0) + decay_component: str = Field(default="Ey") + threshold: float = Field(default=1e-3, gt=0, lt=1, serialization_alias="decay_by") + decay_monitor_port: str | None = Field(default=None) + dft_min_run_time: float = Field( + default=100, + ge=0, + description="Minimum run time after sources for dft_decay mode. " + "Must exceed pulse transit time through the device to avoid " + "false convergence on near-zero fields at output ports.", + ) + + +class SourceConfig(BaseModel): + """Source excitation configuration. + + Controls the Gaussian source bandwidth and which port is excited. + When ``bandwidth`` is ``None`` (auto), ``compute_fwidth`` returns a + bandwidth ~3x wider than the monitor frequency span (matching + gplugins' ``dfcen=0.2`` convention) so edge frequencies receive + adequate spectral power. + """ + + model_config = ConfigDict(validate_assignment=True) + + bandwidth: float | None = Field( + default=None, + description="Source Gaussian bandwidth in wavelength um. None = auto (~3x monitor bw).", + ) + port: str | None = Field( + default=None, + description="Source port name. None = auto-select first port.", + ) + fwidth: float = Field( + default=0.0, + ge=0, + description="Computed source fwidth in frequency units. " + "Set automatically by write_config(); 0 = not yet computed.", + ) + + def compute_fwidth(self, fcen: float, monitor_df: float) -> float: + """Compute Gaussian source fwidth in frequency units. + + When auto (bandwidth=None), returns ``max(3 * monitor_df, 0.2 * fcen)`` + to ensure edge frequencies have enough spectral power. + + Args: + fcen: Center frequency (1/um). + monitor_df: Monitor frequency span (1/um). + + Returns: + Source fwidth in frequency units (1/um). + """ + if self.bandwidth is not None: + # Convert wavelength bandwidth to frequency bandwidth + wl_center = 1.0 / fcen + wl_min = wl_center - self.bandwidth / 2 + wl_max = wl_center + self.bandwidth / 2 + return 1.0 / wl_min - 1.0 / wl_max + return max(3 * monitor_df, 0.2 * fcen) + + +class WavelengthConfig(BaseModel): + """Wavelength and frequency settings for MEEP FDTD simulation. + + MEEP uses normalized units where c = 1 and lengths are in um. + Frequency f = 1/wavelength (in 1/um). + """ + + model_config = ConfigDict(validate_assignment=True) + + wavelength: float = Field(default=1.55, gt=0, description="Center wavelength in um") + bandwidth: float = Field( + default=0.1, ge=0, description="Wavelength bandwidth in um" + ) + num_freqs: int = Field(default=11, ge=1, description="Number of frequency points") + + @computed_field # type: ignore[prop-decorator] + @property + def fcen(self) -> float: + """Center frequency in MEEP units (1/um, since c=1).""" + return 1.0 / self.wavelength + + @computed_field # type: ignore[prop-decorator] + @property + def df(self) -> float: + """Frequency width in MEEP units.""" + wl_min = self.wavelength - self.bandwidth / 2 + wl_max = self.wavelength + self.bandwidth / 2 + f_max = 1.0 / wl_min + f_min = 1.0 / wl_max + return f_max - f_min + + +class ResolutionConfig(BaseModel): + """MEEP grid resolution configuration.""" + + model_config = ConfigDict(validate_assignment=True) + + pixels_per_um: int = Field( + default=32, ge=4, description="Grid pixels per micrometer" + ) + + @classmethod + def coarse(cls) -> ResolutionConfig: + """Coarse resolution (16 pixels/um) for quick tests.""" + return cls(pixels_per_um=16) + + @classmethod + def default(cls) -> ResolutionConfig: + """Default resolution (32 pixels/um).""" + return cls(pixels_per_um=32) + + @classmethod + def fine(cls) -> ResolutionConfig: + """Fine resolution (64 pixels/um) for production runs.""" + return cls(pixels_per_um=64) + + +class AccuracyConfig(BaseModel): + """Controls MEEP subpixel averaging and polygon simplification. + + ``eps_averaging`` toggles MEEP's subpixel smoothing (expensive for + complex polygons). ``subpixel_maxeval`` / ``subpixel_tol`` tune the + integration that implements it. ``simplify_tol`` applies Shapely + polygon simplification before extrusion, dramatically reducing vertex + counts on dense GDS curves. + """ + + model_config = ConfigDict(validate_assignment=True) + + eps_averaging: bool = Field(default=False, description="Toggle subpixel averaging") + subpixel_maxeval: int = Field( + default=0, ge=0, description="Cap on integration evaluations (0=unlimited)" + ) + subpixel_tol: float = Field( + default=1e-4, gt=0, description="Subpixel integration tolerance" + ) + simplify_tol: float = Field( + default=0.0, ge=0, description="Shapely simplification tolerance in um (0=off)" + ) + + +class DiagnosticsConfig(BaseModel): + """Controls diagnostic outputs from the MEEP runner.""" + + model_config = ConfigDict(validate_assignment=True) + + save_geometry: bool = Field(default=True, description="Pre-run geometry cross-section plots") + save_fields: bool = Field(default=True, description="Post-run field snapshot") + save_epsilon_raw: bool = Field( + default=False, description="Raw epsilon .npy (advanced)" + ) + save_animation: bool = Field( + default=False, description="Field animation MP4 (heavy, needs ffmpeg)" + ) + animation_interval: float = Field( + default=0.5, gt=0, + description="MEEP time units between animation frames", + ) + preview_only: bool = Field( + default=False, + description="Init sim and save geometry diagnostics, skip FDTD run", + ) + verbose_interval: float = Field( + default=0, ge=0, + description="MEEP time units between progress prints (0=off)", + ) + + +class PortData(BaseModel): + """Serializable port data for the config JSON.""" + + model_config = ConfigDict(validate_assignment=True) + + name: str + center: list[float] = Field(description="[x, y, z] center coordinates") + orientation: float = Field(description="Port orientation in degrees") + width: float = Field(gt=0) + normal_axis: int = Field(ge=0, le=1, description="0=x, 1=y") + direction: Literal["+", "-"] = Field(description="Direction along normal axis") + is_source: bool = False + + +class LayerStackEntry(BaseModel): + """One layer in the stack config sent to the cloud runner. + + The cloud runner reads these entries alongside layout.gds to + extract polygons per layer and extrude them to 3D prisms. + """ + + model_config = ConfigDict(validate_assignment=True) + + layer_name: str + gds_layer: list[int] = Field(description="[layer_number, datatype]") + zmin: float + zmax: float + material: str + sidewall_angle: float = Field(default=0.0, description="Sidewall angle in degrees") + + +class MaterialData(BaseModel): + """Optical material data for config JSON.""" + + model_config = ConfigDict(validate_assignment=True) + + refractive_index: float = Field(gt=0) + extinction_coeff: float = Field(default=0.0, ge=0) + + +class SimConfig(BaseModel): + """Complete serializable simulation config written as JSON. + + This is the top-level config that the cloud MEEP runner reads. + The geometry is NOT included here — it's in the GDS file. + The layer_stack tells the runner how to extrude each GDS layer. + """ + + model_config = ConfigDict(validate_assignment=True) + + gds_filename: str = Field( + default="layout.gds", description="GDS file with 2D layout" + ) + component_bbox: list[float] | None = Field( + default=None, + description="Original component bbox [xmin, ymin, xmax, ymax] before port extension.", + ) + layer_stack: list[LayerStackEntry] = Field(default_factory=list) + dielectrics: list[dict[str, Any]] = Field(default_factory=list) + ports: list[PortData] = Field(default_factory=list) + materials: dict[str, MaterialData] = Field(default_factory=dict) + wavelength: WavelengthConfig = Field( + default_factory=WavelengthConfig, serialization_alias="fdtd" + ) + source: SourceConfig = Field(default_factory=SourceConfig) + stopping: StoppingConfig = Field(default_factory=StoppingConfig) + resolution: ResolutionConfig = Field(default_factory=ResolutionConfig) + domain: DomainConfig = Field(default_factory=DomainConfig) + accuracy: AccuracyConfig = Field(default_factory=AccuracyConfig) + verbose_interval: float = Field( + default=0, ge=0, description="MEEP time units between progress prints (0=off)" + ) + diagnostics: DiagnosticsConfig = Field(default_factory=DiagnosticsConfig) + symmetries: list[SymmetryEntry] = Field(default_factory=list) + split_chunks_evenly: bool = Field(default=False) + meep_np: int = Field( + default=1, ge=1, description="Recommended MPI process count" + ) + + def to_json(self, path: str | Path) -> Path: + """Write config to JSON file. + + Args: + path: Output file path + + Returns: + Path to the written file + """ + path = Path(path) + path.parent.mkdir(parents=True, exist_ok=True) + path.write_text(json.dumps(self.model_dump(by_alias=True), indent=2)) + return path diff --git a/src/gsim/meep/models/results.py b/src/gsim/meep/models/results.py new file mode 100644 index 00000000..d1ae62b6 --- /dev/null +++ b/src/gsim/meep/models/results.py @@ -0,0 +1,233 @@ +"""Result models for MEEP simulation output.""" + +from __future__ import annotations + +import contextlib +import csv +import json +import logging +from pathlib import Path +from typing import Any + +from pydantic import BaseModel, ConfigDict, Field + +logger = logging.getLogger(__name__) + + +class SParameterResult(BaseModel): + """S-parameter results from MEEP simulation. + + Parses CSV output from the cloud runner and provides + visualization via matplotlib. + """ + + model_config = ConfigDict(validate_assignment=True) + + wavelengths: list[float] = Field(default_factory=list) + s_params: dict[str, list[complex]] = Field( + default_factory=dict, + description="S-param name -> list of complex values per wavelength", + ) + port_names: list[str] = Field(default_factory=list) + debug_info: dict[str, Any] = Field( + default_factory=dict, + description="Eigenmode 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) -> SParameterResult: + """Parse S-parameter results from CSV file. + + Expected CSV format: + wavelength,S11_mag,S11_phase,S21_mag,S21_phase,... + 1.5, 0.1, -30.0, 0.9, 45.0, ... + + Automatically loads ``meep_debug.json`` from the same directory + if it exists, populating the ``debug_info`` field. + + Args: + path: Path to CSV file + + Returns: + SParameterResult instance + """ + import cmath + + path = Path(path) + wavelengths: list[float] = [] + s_params: dict[str, list[complex]] = {} + port_names: set[str] = set() + + with open(path) as f: + reader = csv.DictReader(f) + if reader.fieldnames is None: + return cls() + + # Discover S-param columns + sparam_names: list[str] = [] + for col in reader.fieldnames: + if col.endswith("_mag"): + name = col.removesuffix("_mag") + sparam_names.append(name) + s_params[name] = [] + + for row in reader: + wavelengths.append(float(row["wavelength"])) + for name in sparam_names: + mag = float(row[f"{name}_mag"]) + phase_deg = float(row[f"{name}_phase"]) + phase_rad = cmath.pi * phase_deg / 180.0 + s_params[name].append(cmath.rect(mag, phase_rad)) + + # Extract port names from S-param names (e.g., "S11" -> port 1) + for name in sparam_names: + # S-param format: "Sij" where i,j are port indices + if len(name) >= 3 and name[0] == "S": + indices = name[1:] + for idx_char in indices: + port_names.add(f"port_{idx_char}") + + # 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"), + ]: + 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, + port_names=sorted(port_names), + debug_info=debug_info, + diagnostic_images=diagnostic_images, + ) + + @classmethod + def from_directory(cls, directory: str | Path) -> SParameterResult: + """Load from directory — handles preview-only with no CSV. + + If ``s_parameters.csv`` exists, delegates to ``from_csv()``. + Otherwise loads only debug info and diagnostic images (preview mode). + + Args: + directory: Path to results directory + + Returns: + SParameterResult instance + """ + directory = Path(directory) + csv_path = directory / "s_parameters.csv" + if csv_path.exists(): + return cls.from_csv(csv_path) + + # Preview-only: load debug + images only + debug_info: dict[str, Any] = {} + debug_path = directory / "meep_debug.json" + if debug_path.exists(): + with contextlib.suppress(json.JSONDecodeError, OSError): + debug_info = json.loads(debug_path.read_text()) + + 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"), + ]: + 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, + ) + + def show_diagnostics(self) -> None: + """Display diagnostic images in Jupyter.""" + from IPython.display import Image, display + + if not self.diagnostic_images: + logger.info("No diagnostic images available.") + return + + for name, img_path in sorted(self.diagnostic_images.items()): + if not img_path.endswith((".png", ".jpg", ".jpeg", ".gif")): + continue # skip video/directories; use show_animation() for MP4 + logger.info("--- %s ---", name) + display(Image(filename=img_path)) + + 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 plot(self, db: bool = True, **kwargs: Any) -> Any: + """Plot S-parameters vs wavelength. + + Args: + db: If True, plot in dB scale + **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 = "|S| (dB)" if db else "|S|" + for name, values in self.s_params.items(): + magnitudes = [abs(v) for v in values] + if db: + import math + + y_vals = [20 * math.log10(m) if m > 0 else -100 for m in magnitudes] + else: + y_vals = magnitudes + + 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("S-Parameters") + fig.tight_layout() + + return fig diff --git a/src/gsim/meep/overlay.py b/src/gsim/meep/overlay.py new file mode 100644 index 00000000..1ce51d43 --- /dev/null +++ b/src/gsim/meep/overlay.py @@ -0,0 +1,163 @@ +"""Simulation overlay metadata for 2D visualization. + +Provides dataclasses that describe the simulation cell boundaries, PML regions, +port locations, and dielectric background layers for rendering on top of +geometry cross-sections. +""" + +from __future__ import annotations + +from dataclasses import dataclass, field +from typing import TYPE_CHECKING, Any + +if TYPE_CHECKING: + from gsim.common.geometry_model import GeometryModel + from gsim.meep.models.config import DomainConfig, PortData + + +@dataclass(frozen=True) +class PortOverlay: + """Port location metadata for 2D overlay rendering. + + Attributes: + name: Port name (e.g. "o1"). + center: (x, y, z) center of the port monitor. + width: Transverse width of the port monitor in um. + normal_axis: 0 for x-normal, 1 for y-normal. + direction: "+" or "-" along the normal axis. + is_source: Whether this port is the excitation source. + z_span: Height of the port monitor in z. + """ + + name: str + center: tuple[float, float, float] + width: float + normal_axis: int + direction: str + is_source: bool + z_span: float + + +@dataclass(frozen=True) +class DielectricOverlay: + """Background dielectric slab metadata for 2D overlay rendering. + + Attributes: + name: Dielectric region name (e.g. "oxide", "substrate"). + material: Material name (e.g. "SiO2", "silicon"). + zmin: Bottom z-coordinate of the slab in um. + zmax: Top z-coordinate of the slab in um. + """ + + name: str + material: str + zmin: float + zmax: float + + +@dataclass(frozen=True) +class SimOverlay: + """Simulation cell metadata for 2D visualization overlays. + + Attributes: + cell_min: (xmin, ymin, zmin) of the full simulation cell. + cell_max: (xmax, ymax, zmax) of the full simulation cell. + dpml: PML absorber thickness in um. + ports: List of port overlays for rendering. + dielectrics: List of background dielectric slabs for rendering. + """ + + cell_min: tuple[float, float, float] + cell_max: tuple[float, float, float] + dpml: float + ports: list[PortOverlay] = field(default_factory=list) + dielectrics: list[DielectricOverlay] = field(default_factory=list) + + +def build_sim_overlay( + geometry_model: GeometryModel, + domain_config: DomainConfig, + port_data: list[PortData], + z_span: float | None = None, + dielectrics: list[dict[str, Any]] | None = None, + component_bbox: tuple[float, float, float, float] | None = None, +) -> SimOverlay: + """Build a SimOverlay from geometry model, domain config, and port data. + + Args: + geometry_model: The geometry model providing the geometry bbox. + domain_config: Domain / PML configuration. + port_data: List of PortData objects from port extraction. + z_span: Port monitor z-span. If None, computed from geometry bbox. + dielectrics: List of dielectric dicts from stack.dielectrics. + component_bbox: Original component bbox (xmin, ymin, xmax, ymax) before + port extension. When provided, cell XY boundaries are computed from + this bbox instead of the geometry model bbox (which may include + extended waveguides). + + Returns: + SimOverlay with computed cell boundaries and port overlays. + """ + gmin, gmax = geometry_model.bbox + dpml = domain_config.dpml + margin_xy = domain_config.margin_xy + + # XY: use original component bbox if available (port extension changes geometry bbox) + if component_bbox is not None: + xy_min_x, xy_min_y = component_bbox[0], component_bbox[1] + xy_max_x, xy_max_y = component_bbox[2], component_bbox[3] + else: + xy_min_x, xy_min_y = gmin[0], gmin[1] + xy_max_x, xy_max_y = gmax[0], gmax[1] + + # XY: margin_xy is gap between geometry bbox and PML + # Z: margin_z_above/below is already baked into the geometry bbox via set_z_crop(), + # so only add dpml beyond the geometry z-extent + cell_min = ( + xy_min_x - margin_xy - dpml, + xy_min_y - margin_xy - dpml, + gmin[2] - dpml, + ) + cell_max = ( + xy_max_x + margin_xy + dpml, + xy_max_y + margin_xy + dpml, + gmax[2] + dpml, + ) + + if z_span is None: + z_span = gmax[2] - gmin[2] + + port_margin = domain_config.port_margin + + ports: list[PortOverlay] = [ + PortOverlay( + name=p.name, + center=tuple(p.center), # type: ignore[arg-type] + width=p.width + 2 * port_margin, + normal_axis=p.normal_axis, + direction=p.direction, + is_source=p.is_source, + z_span=z_span, + ) + for p in port_data + ] + + diel_overlays: list[DielectricOverlay] = [] + if dielectrics: + for d in dielectrics: + diel_overlays.append( + DielectricOverlay( + name=d["name"], + material=d["material"], + zmin=d["zmin"], + zmax=d["zmax"], + ) + ) + + return SimOverlay( + cell_min=cell_min, + cell_max=cell_max, + dpml=dpml, + ports=ports, + dielectrics=diel_overlays, + ) diff --git a/src/gsim/meep/ports.py b/src/gsim/meep/ports.py new file mode 100644 index 00000000..3594b77e --- /dev/null +++ b/src/gsim/meep/ports.py @@ -0,0 +1,136 @@ +"""Port extraction for MEEP simulation. + +Extracts port information from a gdsfactory component into a +serializable format for the MEEP config JSON. +""" + +from __future__ import annotations + +from typing import TYPE_CHECKING, Any, Literal + +from gsim.meep.models.config import PortData + +if TYPE_CHECKING: + from gdsfactory.component import Component + + from gsim.common import LayerStack + + +def get_port_normal(orientation: float) -> tuple[int, Literal["+", "-"]]: + """Get port normal axis and direction from orientation angle. + + Args: + orientation: Port orientation in degrees (0, 90, 180, 270) + + Returns: + Tuple of (axis_index, direction) where axis_index is 0=x or 1=y + + Raises: + ValueError: If orientation is not a standard angle + """ + ort = round(orientation) % 360 + if ort == 0: + return 0, "-" + if ort == 90: + return 1, "-" + if ort == 180: + return 0, "+" + if ort == 270: + return 1, "+" + raise ValueError(f"Invalid port orientation: {orientation}") + + +def extract_port_info( + component: Component, + layer_stack: LayerStack, + source_port: str | None = None, +) -> list[PortData]: + """Extract port information from a gdsfactory component. + + Args: + 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. + + Returns: + List of PortData objects ready for JSON serialization + """ + ports: list[PortData] = [] + + z_center = _get_z_center(layer_stack) + + 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 + + ports.append( + PortData( + name=gf_port.name, + center=[ + float(gf_port.center[0]), + float(gf_port.center[1]), + z_center, + ], + orientation=float(gf_port.orientation), + width=float(gf_port.width), + normal_axis=normal_axis, + direction=direction, + is_source=is_source, + ) + ) + + return ports + + +def _find_highest_n_layer(layer_stack: LayerStack) -> tuple[Any, float]: + """Find the layer with the highest refractive index. + + Args: + layer_stack: LayerStack from gsim.common + + Returns: + (best_layer, best_n) -- (None, 0.0) if no optical data. + """ + from gsim.common.stack.materials import get_material_properties + + best_layer = None + best_n = 0.0 + + for layer in layer_stack.layers.values(): + props = get_material_properties(layer.material) + if ( + props is not None + and props.refractive_index is not None + and props.refractive_index > best_n + ): + best_n = props.refractive_index + best_layer = layer + + return best_layer, best_n + + +def _get_z_center(layer_stack: LayerStack) -> float: + """Get z-center for ports from the layer stack. + + For photonic simulation, uses the midpoint of the layer with the + highest refractive index (waveguide core). Falls back to the + midpoint of all layers if no optical data is available. + + Args: + layer_stack: LayerStack from gsim.common + + Returns: + z-center coordinate in um + """ + best_layer, best_n = _find_highest_n_layer(layer_stack) + + if best_layer is not None and best_n > 1.5: + return (best_layer.zmin + best_layer.zmax) / 2.0 + + # Fall back to midpoint of all layers + if not layer_stack.layers: + return 0.0 + all_zmin = min(l.zmin for l in layer_stack.layers.values()) + all_zmax = max(l.zmax for l in layer_stack.layers.values()) + return (all_zmin + all_zmax) / 2.0 diff --git a/src/gsim/meep/script.py b/src/gsim/meep/script.py new file mode 100644 index 00000000..3d08b362 --- /dev/null +++ b/src/gsim/meep/script.py @@ -0,0 +1,1156 @@ +"""MEEP runner script generation. + +Generates a self-contained Python script that runs on the cloud MEEP instance. +The script reads a JSON config + GDS file, builds the MEEP simulation, runs it, +and outputs S-parameters as CSV. +""" + +from __future__ import annotations + +_MEEP_RUNNER_TEMPLATE = '''\ +#!/usr/bin/env python3 +"""Auto-generated MEEP runner script. + +Reads simulation config from JSON + layout from GDS, builds geometry, +runs FDTD, and extracts S-parameters via mode decomposition. + +Cloud dependencies: meep, gdsfactory, numpy, scipy, shapely +""" + +import csv +import cmath +import json +import logging +import math +import sys +import time + +import meep as mp +import numpy as np + +try: + import matplotlib + matplotlib.use("Agg") + import matplotlib.pyplot as plt + HAS_MATPLOTLIB = True +except ImportError: + HAS_MATPLOTLIB = False + +logging.basicConfig(level=logging.INFO, format="%(levelname)s: %(message)s") +logger = logging.getLogger("meep_runner") + + +# --------------------------------------------------------------------------- +# Config loading +# --------------------------------------------------------------------------- + +def load_config(path): + """Load simulation config from JSON file.""" + with open(path) as f: + return json.load(f) + + +# --------------------------------------------------------------------------- +# Materials +# --------------------------------------------------------------------------- + +def build_materials(config): + """Build MEEP material objects from config.""" + materials = {} + for name, props in config["materials"].items(): + n = props["refractive_index"] + k = props.get("extinction_coeff", 0.0) + if k > 0: + materials[name] = mp.Medium(index=n, D_conductivity=2 * cmath.pi * k / n) + else: + materials[name] = mp.Medium(index=n) + return materials + + +# --------------------------------------------------------------------------- +# Geometry: GDS + layer_stack -> mp.Prism objects +# --------------------------------------------------------------------------- + +def load_gds_component(gds_filename): + """Load GDS file as a gdsfactory Component.""" + import gdsfactory as gf + gf.gpdk.PDK.activate() + return gf.import_gds(gds_filename) + + +def _klayout_poly_to_coords(poly_obj, dbu): + """Convert a KLayout polygon object to a list of (x, y) coordinate tuples. + + Args: + poly_obj: KLayout PolygonWithProperties or similar object + dbu: Database unit scaling factor (converts integer coords to microns) + + Returns: + list of (float, float) coordinate tuples + """ + coords = [] + for pt in poly_obj.each_point_hull(): + coords.append((pt.x * dbu, pt.y * dbu)) + return coords + + +def extract_layer_polygons(component, gds_layer, simplify_tol=0.0): + """Extract merged Shapely polygons for a GDS layer from a component. + + Uses gdsfactory / KLayout to get polygons for the given (layer, datatype) + tuple, then converts to Shapely with hole support. + + Handles both old gdsfactory (dict keys are (layer, datatype) tuples, + values are numpy arrays) and new gdsfactory 9.x (dict keys are integer + KLayout layer indices, values are PolygonWithProperties objects). + + Args: + component: gdsfactory Component + gds_layer: [layer_number, datatype] list + simplify_tol: Shapely simplification tolerance in um (0=off) + + Returns: + list of Shapely Polygon objects + """ + from shapely.geometry import Polygon + from shapely.ops import unary_union + + layer_tuple = tuple(gds_layer) + raw = component.get_polygons(layers=(layer_tuple,), merge=True) + + # Collect all polygon objects from the dict, regardless of key type. + # Old gdsfactory: keys are (layer, datatype) tuples, values are ndarray lists + # New gdsfactory 9.x: keys are integer layer indices, values are + # PolygonWithProperties lists + all_objects = [] + if isinstance(raw, dict): + for v in raw.values(): + if isinstance(v, list): + all_objects.extend(v) + else: + all_objects.append(v) + elif raw is not None: + all_objects = list(raw) if hasattr(raw, '__iter__') else [raw] + + dbu = getattr(getattr(component, 'kcl', None), 'dbu', 0.001) + + shapely_polygons = [] + for obj in all_objects: + try: + # KLayout PolygonWithProperties — convert via each_point_hull + if hasattr(obj, 'each_point_hull'): + coords = _klayout_poly_to_coords(obj, dbu) + # Numpy array — legacy path + elif hasattr(obj, 'shape') and len(obj) >= 3: + coords = [(float(p[0]), float(p[1])) for p in obj] + else: + continue + + if len(coords) < 3: + continue + poly = Polygon(coords) + if poly.is_valid and not poly.is_empty: + shapely_polygons.append(poly) + except Exception: + continue + + if not shapely_polygons: + return [] + + merged = unary_union(shapely_polygons) + if simplify_tol > 0: + merged = merged.simplify(simplify_tol, preserve_topology=True) + if hasattr(merged, "geoms"): + return list(merged.geoms) + return [merged] if not merged.is_empty else [] + + +def triangulate_polygon_with_holes(polygon): + """Triangulate a polygon with holes using Delaunay triangulation. + + Returns list of triangles as coordinate lists, each triangle being + [[x1,y1], [x2,y2], [x3,y3]]. + """ + from scipy.spatial import Delaunay + from shapely.geometry import Point + + all_points = list(polygon.exterior.coords[:-1]) + for interior in polygon.interiors: + all_points.extend(list(interior.coords[:-1])) + + if len(all_points) < 3: + return [list(polygon.exterior.coords[:-1])] + + points_2d = np.array(all_points) + tri = Delaunay(points_2d) + + triangles = [] + for simplex in tri.simplices: + triangle_pts = points_2d[simplex] + centroid = np.mean(triangle_pts, axis=0) + if polygon.contains(Point(centroid)): + triangles.append(triangle_pts.tolist()) + + return triangles if triangles else [list(polygon.exterior.coords[:-1])] + + +def build_background_slabs(config, materials): + """Build background mp.Block slabs from dielectric entries. + + These infinite-XY slabs fill the simulation cell at each z-range with + the correct cladding/substrate material. They must come FIRST in the + geometry list so that patterned prisms (added later) take precedence. + """ + slabs = [] + for diel in sorted(config.get("dielectrics", []), key=lambda d: d["zmin"]): + mat = materials.get(diel["material"]) + if mat is None: + continue + zmin = diel["zmin"] + zmax = diel["zmax"] + thickness = zmax - zmin + if thickness <= 0: + continue + block = mp.Block( + size=mp.Vector3(mp.inf, mp.inf, thickness), + center=mp.Vector3(0, 0, (zmin + zmax) / 2), + material=mat, + ) + slabs.append(block) + return slabs + + +def build_symmetries(config): + """Build MEEP symmetry objects from config.""" + direction_map = {"X": mp.X, "Y": mp.Y, "Z": mp.Z} + symmetries = [] + for sym in config.get("symmetries", []): + symmetries.append( + mp.Mirror(direction_map[sym["direction"]], phase=sym.get("phase", 1)) + ) + return symmetries + + +# Map string component names to MEEP field components +_COMPONENT_MAP = { + "Ex": mp.Ex, "Ey": mp.Ey, "Ez": mp.Ez, + "Hx": mp.Hx, "Hy": mp.Hy, "Hz": mp.Hz, +} + + +def resolve_decay_monitor_point(config): + """Return center mp.Vector3 for decay monitoring. + + Uses the port named by stopping.decay_monitor_port if set, + otherwise picks the first non-source port. Falls back to + the first port if all are sources. + """ + stopping = config.get("stopping", {}) + target_name = stopping.get("decay_monitor_port") + ports = config.get("ports", []) + + if target_name: + for p in ports: + if p["name"] == target_name: + return mp.Vector3(*p["center"]) + + # Fallback: first non-source port, or first port + for p in ports: + if not p.get("is_source", False): + return mp.Vector3(*p["center"]) + if ports: + return mp.Vector3(*ports[0]["center"]) + + return mp.Vector3(0, 0, 0) + + +def build_geometry(config, materials): + """Build MEEP geometry from GDS file + layer stack config. + + For each layer in layer_stack: + 1. Extract polygons from the GDS for that layer's gds_layer + 2. Extrude to 3D as mp.Prism with correct z-range and material + 3. Handle polygon holes via Delaunay triangulation + """ + gds_filename = config.get("gds_filename", "layout.gds") + component = load_gds_component(gds_filename) + + accuracy = config.get("accuracy", {}) + simplify_tol = accuracy.get("simplify_tol", 0.0) + + geometry = [] + total_vertices = 0 + + for layer_entry in config["layer_stack"]: + material_name = layer_entry["material"] + mat = materials.get(material_name, mp.Medium()) + zmin = layer_entry["zmin"] + zmax = layer_entry["zmax"] + height = zmax - zmin + gds_layer = layer_entry["gds_layer"] + sidewall_angle_deg = layer_entry.get("sidewall_angle", 0.0) + sw_rad = math.radians(sidewall_angle_deg) if sidewall_angle_deg else 0 + + if height <= 0: + continue + + polygons = extract_layer_polygons( + component, gds_layer, simplify_tol=simplify_tol, + ) + + for polygon in polygons: + if polygon.is_empty or not polygon.is_valid: + continue + + if hasattr(polygon, "interiors") and polygon.interiors: + # Polygon has holes — triangulate + triangles = triangulate_polygon_with_holes(polygon) + for tri_coords in triangles: + vertices = [mp.Vector3(p[0], p[1], zmin) for p in tri_coords] + total_vertices += len(vertices) + prism = mp.Prism( + vertices=vertices, + height=height, + material=mat, + sidewall_angle=sw_rad, + ) + geometry.append(prism) + else: + # Simple polygon — direct extrusion + coords = list(polygon.exterior.coords[:-1]) + vertices = [mp.Vector3(p[0], p[1], zmin) for p in coords] + total_vertices += len(vertices) + prism = mp.Prism( + vertices=vertices, + height=height, + material=mat, + sidewall_angle=sw_rad, + ) + geometry.append(prism) + + logger.info("Total vertices across all prisms: %d", total_vertices) + return geometry, component + + +# --------------------------------------------------------------------------- +# Sources and monitors +# --------------------------------------------------------------------------- + +def get_port_z_span(config): + """Get z-span for ports from layer stack.""" + zmin = min(l["zmin"] for l in config["layer_stack"]) + zmax = max(l["zmax"] for l in config["layer_stack"]) + return zmax - zmin + + +def build_sources(config): + """Build MEEP source from config port data.""" + fdtd = config["fdtd"] + fcen = fdtd["fcen"] + df = fdtd["df"] + fwidth = config["source"]["fwidth"] + z_span = get_port_z_span(config) + port_margin = config.get("domain", {}).get("port_margin", 0.5) + + sources = [] + for port in config["ports"]: + if not port["is_source"]: + continue + + center = mp.Vector3(*port["center"]) + normal_axis = port["normal_axis"] + direction = port["direction"] + + size = [0, 0, 0] + transverse_axis = 1 - normal_axis + size[transverse_axis] = port["width"] + 2 * port_margin + size[2] = z_span + + # Propagation axis + prop_axis = mp.X if normal_axis == 0 else mp.Y + + # eig_kpoint points INTO the device (direction of incoming mode) + if normal_axis == 0: + kpoint = mp.Vector3(x=1) if direction == "+" else mp.Vector3(x=-1) + else: + kpoint = mp.Vector3(y=1) if direction == "+" else mp.Vector3(y=-1) + + eig_src = mp.EigenModeSource( + src=mp.GaussianSource(frequency=fcen, fwidth=fwidth), + center=center, + size=mp.Vector3(*size), + eig_band=1, + direction=prop_axis, + eig_kpoint=kpoint, + eig_match_freq=True, + ) + sources.append(eig_src) + + return sources + + +def build_monitors(config, sim): + """Build mode monitors at all ports and return flux regions.""" + fdtd = config["fdtd"] + fcen = fdtd["fcen"] + df = fdtd["df"] + nfreq = fdtd["num_freqs"] + z_span = get_port_z_span(config) + port_margin = config.get("domain", {}).get("port_margin", 0.5) + + monitors = {} + for port in config["ports"]: + center = mp.Vector3(*port["center"]) + normal_axis = port["normal_axis"] + + size = [0, 0, 0] + transverse_axis = 1 - normal_axis + size[transverse_axis] = port["width"] + 2 * port_margin + size[2] = z_span + + flux = sim.add_mode_monitor( + fcen, df, nfreq, + mp.ModeRegion(center=center, size=mp.Vector3(*size)), + ) + monitors[port["name"]] = flux + + return monitors + + +# --------------------------------------------------------------------------- +# S-parameter extraction +# --------------------------------------------------------------------------- + +def _port_kpoint(port): + """Return positive-axis kpoint vector for eigenmode decomposition. + + Always points along +normal_axis so that: + alpha[band, freq, 0] = forward (+axis) coefficient + alpha[band, freq, 1] = backward (-axis) coefficient + """ + if port["normal_axis"] == 0: + return mp.Vector3(x=1) + return mp.Vector3(y=1) + + +def extract_s_params(config, sim, monitors): + """Extract S-parameters via mode decomposition. + + Uses MEEP eigenmode coefficients with explicit kpoint_func to anchor + the forward/backward convention: + alpha[band, freq, 0] = forward (+normal_axis) coefficient + alpha[band, freq, 1] = backward (-normal_axis) coefficient + + Port "direction" field = direction of incoming mode along normal_axis: + "+" -> incoming goes +normal, outgoing (reflected/transmitted) goes -normal + "-" -> incoming goes -normal, outgoing (reflected/transmitted) goes +normal + + Returns: + (s_params, debug_data) tuple where debug_data contains eigenmode + diagnostics for post-run analysis. + """ + port_names = [p["name"] for p in config["ports"]] + ports = {p["name"]: p for p in config["ports"]} + source_port = None + for p in config["ports"]: + if p["is_source"]: + source_port = p["name"] + break + + if not source_port: + return {}, {} + + fdtd = config["fdtd"] + fcen = fdtd["fcen"] + nfreq = fdtd["num_freqs"] + df = fdtd["df"] + freqs = np.linspace(fcen - df / 2, fcen + df / 2, nfreq) + + debug_data = { + "eigenmode_info": {}, + "raw_coefficients": {}, + "incident_coefficients": {}, + } + + def _incoming_idx(direction): + """Alpha index for the incoming (incident) mode at a port.""" + return 0 if direction == "+" else 1 + + def _outgoing_idx(direction): + """Alpha index for the outgoing (reflected/transmitted) mode at a port.""" + return 1 if direction == "+" else 0 + + def _collect_eigenmode_debug(port_name, ob, freqs): + """Collect eigenmode diagnostics from coefficients object.""" + info = {"band": 1} + try: + kdom = ob.kdom + kdom_list = [[float(kdom[i].x), float(kdom[i].y), float(kdom[i].z)] + for i in range(len(kdom))] + info["kdom"] = kdom_list + info["n_eff"] = [float(np.linalg.norm(kdom_list[i])) / float(freqs[i]) + if float(freqs[i]) > 0 else 0.0 + for i in range(min(len(kdom_list), len(freqs)))] + except Exception: + info["kdom"] = [] + info["n_eff"] = [] + try: + cg = ob.cg + info["group_velocity"] = [float(cg[i]) for i in range(len(cg))] + except Exception: + info["group_velocity"] = [] + return info + + # Get incident coefficient at source port for normalization + src_dir = ports[source_port]["direction"] + src_kp = _port_kpoint(ports[source_port]) + src_ob = sim.get_eigenmode_coefficients( + monitors[source_port], [1], eig_parity=mp.NO_PARITY, + kpoint_func=lambda f, n, kp=src_kp: kp, + ) + incident_coeffs = src_ob.alpha[0, :, _incoming_idx(src_dir)] + + debug_data["eigenmode_info"][source_port] = _collect_eigenmode_debug( + source_port, src_ob, freqs) + debug_data["incident_coefficients"] = { + "port": source_port, + "magnitudes": [float(abs(c)) for c in incident_coeffs], + } + debug_data["raw_coefficients"][source_port] = { + "forward_mag": [float(abs(src_ob.alpha[0, i, 0])) for i in range(len(freqs))], + "backward_mag": [float(abs(src_ob.alpha[0, i, 1])) for i in range(len(freqs))], + } + + s_params = {} + + for i, port_i in enumerate(port_names): + for j, port_j in enumerate(port_names): + if port_j != source_port: + continue + + s_name = f"S{i+1}{j+1}" + + port_kp = _port_kpoint(ports[port_i]) + ob = sim.get_eigenmode_coefficients( + monitors[port_i], [1], eig_parity=mp.NO_PARITY, + kpoint_func=lambda f, n, kp=port_kp: kp, + ) + + # Collect debug info (skip source port — already collected) + if port_i != source_port: + debug_data["eigenmode_info"][port_i] = _collect_eigenmode_debug( + port_i, ob, freqs) + nf = len(freqs) + debug_data["raw_coefficients"][port_i] = { + "forward_mag": [ + float(abs(ob.alpha[0, k, 0])) for k in range(nf) + ], + "backward_mag": [ + float(abs(ob.alpha[0, k, 1])) for k in range(nf) + ], + } + + # Outgoing direction: reflected at source port, transmitted at output ports + port_dir = ports[port_i]["direction"] + alpha = ob.alpha[0, :, _outgoing_idx(port_dir)] + + s_params[s_name] = alpha / incident_coeffs + + # Power conservation: sum |Sij|^2 per frequency + power_conservation = np.zeros(len(freqs)) + for s_name, s_vals in s_params.items(): + power_conservation += np.abs(s_vals) ** 2 + debug_data["power_conservation"] = [float(p) for p in power_conservation] + + return s_params, debug_data + + +def save_results(config, s_params, output_path="s_parameters.csv"): + """Save S-parameters to CSV file. 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 + + fieldnames = ["wavelength"] + for name in sorted(s_params.keys()): + fieldnames.extend([f"{name}_mag", f"{name}_phase"]) + + 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 sorted(s_params.keys()): + val = s_params[name][i] + row[f"{name}_mag"] = f"{abs(val):.6f}" + row[f"{name}_phase"] = f"{cmath.phase(val) * 180 / cmath.pi:.4f}" + writer.writerow(row) + + logger.info("S-parameters saved to %s", output_path) + + +def save_debug_log(config, s_params, debug_data, wall_seconds=0.0, + output_path="meep_debug.json"): + """Save eigenmode diagnostics as JSON for post-run analysis. Only rank 0 writes.""" + if not mp.am_master(): + return + fdtd = config["fdtd"] + domain = config.get("domain", {}) + resolution = config.get("resolution", {}).get("pixels_per_um", 0) + stopping = config.get("stopping", {}) + + meep_time = 0.0 + timesteps = 0 + try: + # These are set by caller if available + meep_time = debug_data.get("_meep_time", 0.0) + timesteps = debug_data.get("_timesteps", 0) + except Exception: + pass + + log = { + "metadata": { + "resolution": resolution, + "cell_size": debug_data.get("_cell_size", []), + "meep_time": meep_time, + "timesteps": timesteps, + "wall_seconds": wall_seconds, + "stopping_mode": stopping.get("mode", "fixed"), + }, + "eigenmode_info": debug_data.get("eigenmode_info", {}), + "incident_coefficients": debug_data.get("incident_coefficients", {}), + "raw_coefficients": debug_data.get("raw_coefficients", {}), + "power_conservation": debug_data.get("power_conservation", []), + } + + with open(output_path, "w") as f: + json.dump(log, f, indent=2) + logger.info("Debug log saved to %s", output_path) + + +# --------------------------------------------------------------------------- +# Diagnostics +# --------------------------------------------------------------------------- + +def save_geometry_diagnostics(sim, config, cell_center): + """Save geometry cross-section plots showing epsilon, sources, monitors, PML.""" + if not HAS_MATPLOTLIB: + logger.warning("matplotlib not available, skipping geometry diagnostics") + return + if not mp.am_master(): + # plot2D is collective — all ranks call it, only master saves + pass + + domain = config.get("domain", {}) + dpml = domain.get("dpml", 1.0) + z_min = min(l["zmin"] for l in config["layer_stack"]) + z_max = max(l["zmax"] for l in config["layer_stack"]) + z_core = (z_min + z_max) / 2 + + # XY cross-section at z=core center + try: + fig, ax = plt.subplots(1, 1, figsize=(10, 8)) + xy_plane = mp.Volume( + center=mp.Vector3(cell_center.x, cell_center.y, z_core), + size=mp.Vector3(sim.cell_size.x, sim.cell_size.y, 0), + ) + sim.plot2D(ax=ax, output_plane=xy_plane) + ax.set_title(f"XY cross-section at z={z_core:.3f} um (core center)") + ax.set_xlabel("x (um)") + ax.set_ylabel("y (um)") + fig.tight_layout() + if mp.am_master(): + fig.savefig("meep_geometry_xy.png", dpi=150) + logger.info("Saved meep_geometry_xy.png") + plt.close(fig) + except Exception as e: + logger.warning("XY geometry plot failed: %s", e) + + # For 3D sims: XZ and YZ cross-sections + if sim.cell_size.z > 0: + # XZ at y=center + try: + fig, ax = plt.subplots(1, 1, figsize=(10, 6)) + xz_plane = mp.Volume( + center=mp.Vector3(cell_center.x, cell_center.y, cell_center.z), + size=mp.Vector3(sim.cell_size.x, 0, sim.cell_size.z), + ) + sim.plot2D(ax=ax, output_plane=xz_plane) + ax.set_title(f"XZ cross-section at y={cell_center.y:.3f} um") + ax.set_xlabel("x (um)") + ax.set_ylabel("z (um)") + fig.tight_layout() + if mp.am_master(): + fig.savefig("meep_geometry_xz.png", dpi=150) + logger.info("Saved meep_geometry_xz.png") + plt.close(fig) + except Exception as e: + logger.warning("XZ geometry plot failed: %s", e) + + # YZ at x=center + try: + fig, ax = plt.subplots(1, 1, figsize=(8, 6)) + yz_plane = mp.Volume( + center=mp.Vector3(cell_center.x, cell_center.y, cell_center.z), + size=mp.Vector3(0, sim.cell_size.y, sim.cell_size.z), + ) + sim.plot2D(ax=ax, output_plane=yz_plane) + ax.set_title(f"YZ cross-section at x={cell_center.x:.3f} um") + ax.set_xlabel("y (um)") + ax.set_ylabel("z (um)") + fig.tight_layout() + if mp.am_master(): + fig.savefig("meep_geometry_yz.png", dpi=150) + logger.info("Saved meep_geometry_yz.png") + plt.close(fig) + except Exception as e: + logger.warning("YZ geometry plot failed: %s", e) + + +def save_field_snapshot(sim, config, cell_center): + """Save post-run field snapshot (Ey overlaid on epsilon).""" + if not HAS_MATPLOTLIB: + logger.warning("matplotlib not available, skipping field snapshot") + return + + z_min = min(l["zmin"] for l in config["layer_stack"]) + z_max = max(l["zmax"] for l in config["layer_stack"]) + z_core = (z_min + z_max) / 2 + + try: + fig, ax = plt.subplots(1, 1, figsize=(10, 8)) + xy_plane = mp.Volume( + center=mp.Vector3(cell_center.x, cell_center.y, z_core), + size=mp.Vector3(sim.cell_size.x, sim.cell_size.y, 0), + ) + sim.plot2D(ax=ax, output_plane=xy_plane, fields=mp.Ey) + ax.set_title(f"Ey field at z={z_core:.3f} um (post-run)") + ax.set_xlabel("x (um)") + ax.set_ylabel("y (um)") + fig.tight_layout() + if mp.am_master(): + fig.savefig("meep_fields_xy.png", dpi=150) + logger.info("Saved meep_fields_xy.png") + plt.close(fig) + except Exception as e: + 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). + + Saves a compressed .npz with the Ey field array and timestamp. + Rendering with a globally fixed colorbar happens after sim.run(). + + Returns: + Incremented frame counter. + """ + field_data = np.real(sim.get_array(vol=xy_plane, component=mp.Ey)) + t = sim.meep_time() + if mp.am_master(): + np.savez_compressed( + f"meep_field_{frame_counter:04d}.npz", + field=field_data, + time=t, + ) + return frame_counter + 1 + + +def render_animation_frames(eps_data, extent): + """Render saved field .npz files into PNGs with fixed global colorbar. + + Two-pass: first finds the global field maximum across all frames, + then renders every frame with the same vmin/vmax so field decay is + clearly visible. + + Call only on master rank after sim.run(). + """ + import glob + import os + + if not HAS_MATPLOTLIB: + logger.warning( + "matplotlib not available, .npz field files kept but not rendered" + ) + return + + npz_files = sorted(glob.glob("meep_field_*.npz")) + if not npz_files: + logger.warning("No field data files found to render") + 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"])))) + if global_max == 0: + global_max = 1.0 + + logger.info( + "Rendering %d frames (Ey global max = %.4g) ...", + len(npz_files), 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"]) + + 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, + vmin=-global_max, vmax=global_max, + ) + 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)") + fig.tight_layout() + fig.savefig(f"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)) + + +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. + """ + import glob + import subprocess + + frames = sorted(glob.glob("meep_frame_*.png")) + if not frames: + logger.warning("No animation frames found to compile") + return + + logger.info("Compiling %d frames into meep_animation.mp4 ...", len(frames)) + try: + subprocess.run( + [ + "ffmpeg", "-y", + "-framerate", str(fps), + "-i", "meep_frame_%04d.png", + "-c:v", "libx264", + "-pix_fmt", "yuv420p", + "meep_animation.mp4", + ], + check=True, + capture_output=True, + ) + logger.info("Saved meep_animation.mp4") + except FileNotFoundError: + logger.warning("ffmpeg not found — frame PNGs saved but MP4 not created") + except subprocess.CalledProcessError as e: + logger.warning("ffmpeg failed: %s", e.stderr.decode()[:500]) + logger.info("Frame PNGs are still available as meep_frame_*.png") + + +def save_epsilon_raw(sim, config, cell_center): + """Save raw epsilon array as .npy for XY slice at core center.""" + z_min = min(l["zmin"] for l in config["layer_stack"]) + z_max = max(l["zmax"] for l in config["layer_stack"]) + z_core = (z_min + z_max) / 2 + + try: + xy_plane = mp.Volume( + center=mp.Vector3(cell_center.x, cell_center.y, z_core), + size=mp.Vector3(sim.cell_size.x, sim.cell_size.y, 0), + ) + eps_data = sim.get_array(vol=xy_plane, component=mp.Dielectric) + if mp.am_master(): + np.save("meep_epsilon_xy.npy", eps_data) + logger.info("Saved meep_epsilon_xy.npy (shape=%s)", eps_data.shape) + except Exception as e: + logger.warning("Epsilon raw save failed: %s", e) + + +# --------------------------------------------------------------------------- +# Main +# --------------------------------------------------------------------------- + +def main(): + """Run MEEP simulation.""" + config_path = "%%CONFIG_FILENAME%%" + config = load_config(config_path) + + logger.info("Building materials...") + materials = build_materials(config) + + logger.info("Building background dielectric slabs...") + background_slabs = build_background_slabs(config, materials) + logger.info("Created %d background slabs", len(background_slabs)) + + logger.info("Building geometry from GDS + layer stack...") + geometry, component = build_geometry(config, materials) + logger.info( + "Created %d prisms from %d layers", + len(geometry), len(config['layer_stack']), + ) + + # 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) + + resolution = config["resolution"]["pixels_per_um"] + fdtd = config["fdtd"] + + # Compute simulation cell from component bounds + layer z-range + # Use original component bbox if available (port extension changes GDS bbox) + component_bbox = config.get("component_bbox") + if component_bbox is not None: + bbox_left, bbox_bottom, bbox_right, bbox_top = component_bbox + else: + bbox = component.dbbox() + bbox_left, bbox_right = bbox.left, bbox.right + bbox_bottom, bbox_top = bbox.bottom, bbox.top + + z_min = min(l["zmin"] for l in config["layer_stack"]) + z_max = max(l["zmax"] for l in config["layer_stack"]) + + domain = config.get("domain", {}) + dpml = domain.get("dpml", 1.0) + margin_xy = domain.get("margin_xy", 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) + # Z: margin_z_above/below is already baked into layer_stack via set_z_crop(), + # so only add dpml beyond the stack extent + cell_z = (z_max - z_min) + 2 * dpml + cell_center = mp.Vector3( + (bbox_right + bbox_left) / 2, + (bbox_top + bbox_bottom) / 2, + (z_max + z_min) / 2, + ) + + 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) + + accuracy = config.get("accuracy", {}) + diagnostics = config.get("diagnostics", {}) + preview_only = diagnostics.get("preview_only", False) + + # In preview mode, skip expensive subpixel averaging + eps_avg = False if preview_only else accuracy.get("eps_averaging", True) + + # Symmetries are NOT used for S-parameter extraction runs. + # MEEP's get_eigenmode_coefficients with add_mode_monitor produces + # incorrect normalization when the source monitor straddles a mirror + # symmetry plane (coefficients underestimated by ~2x). This is + # consistent with gplugins which also never uses mp.Mirror for + # S-parameter extraction. Symmetries are only applied in preview-only + # mode (geometry validation, no FDTD/S-params). + cfg_symmetries = config.get("symmetries", []) + if cfg_symmetries and not preview_only: + logger.info("Symmetries present in config but IGNORED for S-parameter " + "extraction (causes incorrect eigenmode normalization). " + "Symmetries are only used in preview-only mode.") + use_symmetries = build_symmetries(config) if preview_only else [] + + sim_kwargs = dict( + cell_size=mp.Vector3(cell_x, cell_y, cell_z), + geometry_center=cell_center, + geometry=geometry, + sources=sources, + resolution=resolution, + boundary_layers=[mp.PML(dpml)], + symmetries=use_symmetries, + split_chunks_evenly=config.get("split_chunks_evenly", False), + eps_averaging=eps_avg, + ) + spx_maxeval = accuracy.get("subpixel_maxeval", 0) + if spx_maxeval > 0: + sim_kwargs["subpixel_maxeval"] = spx_maxeval + spx_tol = accuracy.get("subpixel_tol", 1e-4) + if spx_tol != 1e-4: + sim_kwargs["subpixel_tol"] = spx_tol + sim = mp.Simulation(**sim_kwargs) + + # --- Diagnostics & preview mode --- + diag_geometry = diagnostics.get("save_geometry", True) + diag_fields = diagnostics.get("save_fields", True) + diag_epsilon = diagnostics.get("save_epsilon_raw", False) + + if diag_geometry or diag_epsilon or preview_only: + logger.info("Initializing simulation for diagnostics...") + sim.init_sim() + if diag_geometry or preview_only: + save_geometry_diagnostics(sim, config, cell_center) + if diag_epsilon: + save_epsilon_raw(sim, config, cell_center) + + if preview_only: + logger.info("MEEP_PREVIEW_ONLY=1 — skipping simulation run.") + save_debug_log(config, {}, {"_meep_time": 0, "_timesteps": 0, + "_cell_size": [cell_x, cell_y, cell_z]}) + logger.info("Preview complete.") + sys.exit(0) + + logger.info("Building monitors...") + monitors = build_monitors(config, sim) + + stopping = config.get("stopping", {}) + run_after = stopping.get("run_after_sources", 100) + + # Build verbose step functions + step_funcs = [] + verbose_interval = config.get("verbose_interval", 0) + if verbose_interval > 0: + _wall_start = time.time() + def _verbose_print(sim_obj): + elapsed = time.time() - _wall_start + logger.info("t=%.2f | wall=%.1fs", sim_obj.meep_time(), elapsed) + step_funcs.append(mp.at_every(verbose_interval, _verbose_print)) + + # Animation field capture step function (raw data, no plotting yet) + diag_animation = diagnostics.get("save_animation", False) + animation_interval = diagnostics.get("animation_interval", 0.5) + _frame_counter = [0] # mutable container for closure + _anim_plane = None + + 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), + ) + + def _capture_frame(sim_obj): + _frame_counter[0] = save_animation_field( + sim_obj, _anim_plane, _frame_counter[0] + ) + + step_funcs.append(mp.at_every(animation_interval, _capture_frame)) + logger.info( + "Animation: saving field data every %s time units", + animation_interval, + ) + + stop_mode = stopping.get("mode", "fixed") + + wall_start = time.time() + + if stop_mode == "dft_decay": + decay_by = stopping.get("decay_by", 1e-3) + min_time = stopping.get("dft_min_run_time", 100) + logger.info( + "Running simulation (dft_decay mode: tol=%s, " + "min=%.1f, max=%.1f)...", + decay_by, min_time, run_after, + ) + sim.run( + *step_funcs, + until_after_sources=mp.stop_when_dft_decayed( + tol=decay_by, + minimum_run_time=min_time, + maximum_run_time=run_after, + ), + ) + elif stop_mode == "decay": + dt = stopping.get("decay_dt", 50.0) + comp_name = stopping.get("decay_component", "Ey") + comp = _COMPONENT_MAP.get(comp_name, mp.Ey) + decay_by = stopping.get("decay_by", 1e-3) + monitor_pt = resolve_decay_monitor_point(config) + logger.info("Running simulation (decay mode: component=%s, dt=%s, " + "decay_by=%s, cap=%.1f)...", comp_name, dt, decay_by, run_after) + + # Decay condition + numeric time cap (list = OR logic, first wins) + decay_fn = mp.stop_when_fields_decayed(dt, comp, monitor_pt, decay_by) + sim.run(*step_funcs, until_after_sources=[decay_fn, run_after]) + else: + logger.info("Running simulation (until_after_sources=%.1f)...", run_after) + sim.run(*step_funcs, until_after_sources=run_after) + + wall_seconds = time.time() - wall_start + + 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 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) + + # 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) + logger.info("Done!") + + +if __name__ == "__main__": + main() +''' + + +def generate_meep_script(config_filename: str = "sim_config.json") -> str: + """Generate the Python script that runs on the cloud MEEP instance. + + The script is self-contained and reads a JSON config + GDS file to: + 1. Load GDS layout via gdsfactory + 2. Extract polygons per layer, handling holes via triangulation + 3. Create mp.Prism objects with correct materials and sidewall angles + 4. Create EigenModeSource at the source port + 5. Create mode monitors at all ports + 6. Build mp.Simulation with PML boundaries + 7. Run and extract S-parameters + 8. Save results as CSV + + Args: + config_filename: Name of the JSON config file (default: sim_config.json) + + Returns: + Python script as a string + """ + return _MEEP_RUNNER_TEMPLATE.replace("%%CONFIG_FILENAME%%", config_filename) diff --git a/src/gsim/meep/simulation.py b/src/gsim/meep/simulation.py new file mode 100644 index 00000000..ea938ae8 --- /dev/null +++ b/src/gsim/meep/simulation.py @@ -0,0 +1,699 @@ +"""Declarative Simulation container for MEEP photonic FDTD. + +Translates the user-facing declarative API objects into the existing +``SimConfig`` JSON contract consumed by the cloud runner. +""" + +from __future__ import annotations + +import logging +import math +from dataclasses import dataclass +from pathlib import Path +from typing import Any + +from pydantic import BaseModel, ConfigDict, Field, PrivateAttr, field_validator + +from gsim.meep.models.api import ( + FDTD, + Domain, + Geometry, + Material, + ModeSource, +) + +logger = logging.getLogger(__name__) + + +# --------------------------------------------------------------------------- +# BuildResult +# --------------------------------------------------------------------------- + + +@dataclass +class BuildResult: + """Result of :meth:`Simulation.build_config` — single source of truth. + + Attributes: + config: Full serializable SimConfig. + component: Extended component (what meep actually simulates). + original_component: Original component before port extension. + """ + + config: Any # SimConfig + component: Any # gdsfactory Component (extended) + original_component: Any # gdsfactory Component (original) + + +# --------------------------------------------------------------------------- +# Utility +# --------------------------------------------------------------------------- + + +def estimate_meep_np( + cell_x: float, + cell_y: float, + cell_z: float, + pixels_per_um: int, + *, + max_cores: int = 24, + min_voxels_per_proc: int = 200_000, +) -> int: + """Estimate optimal MPI process count from problem size. + + Args: + cell_x: Cell size along X in um. + cell_y: Cell size along Y in um. + cell_z: Cell size along Z in um. + pixels_per_um: Grid resolution (pixels per um). + max_cores: Maximum physical cores available (default 24 for c6i.12xlarge). + min_voxels_per_proc: Minimum voxels per process to keep MPI overhead < ~5%. + + Returns: + Recommended number of MPI processes (1 to max_cores). + """ + nx = math.ceil(cell_x * pixels_per_um) + ny = math.ceil(cell_y * pixels_per_um) + nz = max(1, math.ceil(cell_z * pixels_per_um)) + total_voxels = nx * ny * nz + np_est = total_voxels // min_voxels_per_proc + return max(1, min(np_est, max_cores)) + + +class Simulation(BaseModel): + """Declarative MEEP FDTD simulation container. + + Assigns typed physics objects, then calls ``write_config()`` to + produce the JSON + GDS + runner consumed by the cloud engine. + + Example:: + + from gsim import meep + + sim = meep.Simulation() + sim.geometry.component = ybranch + sim.geometry.stack = stack + sim.materials = {"si": 3.47, "sio2": 1.44} + sim.source.port = "o1" + sim.monitors = ["o1", "o2"] + sim.solver.stopping = "dft_decay" + sim.solver.max_time = 200 + result = sim.run("./meep-sim") + """ + + model_config = ConfigDict( + validate_assignment=True, + arbitrary_types_allowed=True, + ) + + geometry: Geometry = Field(default_factory=Geometry) + materials: dict[str, float | Material] = Field(default_factory=dict) + source: ModeSource = Field(default_factory=ModeSource) + monitors: list[str] = Field(default_factory=list) + domain: Domain = Field(default_factory=Domain) + solver: FDTD = Field(default_factory=FDTD) + + # Private: kwargs captured from geometry.stack when it's a string/path + _stack_kwargs: dict[str, Any] = PrivateAttr(default_factory=dict) + + # ------------------------------------------------------------------------- + # Validators + # ------------------------------------------------------------------------- + + @field_validator("materials", mode="before") + @classmethod + def _normalize_materials( + cls, + v: dict[str, float | Material | dict], + ) -> dict[str, float | Material]: + """Accept float shorthand: ``{"si": 3.47}`` → ``Material(n=3.47)``.""" + out: dict[str, float | Material] = {} + for name, val in v.items(): + if isinstance(val, (int, float)): + out[name] = Material(n=float(val)) + elif isinstance(val, dict): + out[name] = Material(**val) + else: + out[name] = val + return out + + # ------------------------------------------------------------------------- + # Resolved materials helper + # ------------------------------------------------------------------------- + + def _resolved_materials(self) -> dict[str, Material]: + """Return materials dict with all values normalized to Material.""" + out: dict[str, Material] = {} + for name, val in self.materials.items(): + if isinstance(val, (int, float)): + out[name] = Material(n=float(val)) + else: + out[name] = val + return out + + # ------------------------------------------------------------------------- + # Validation + # ------------------------------------------------------------------------- + + def validate_config(self) -> Any: + """Validate the simulation configuration. + + Returns: + ValidationResult with errors/warnings. + """ + from gsim.common import ValidationResult + + errors: list[str] = [] + warnings_list: list[str] = [] + + if self.geometry.component is None: + errors.append("No component set. Assign sim.geometry.component first.") + + if self.geometry.component is not None: + ports = list(self.geometry.component.ports) + if not ports: + errors.append("Component has no ports.") + elif self.source.port is not None: + port_names = [p.name for p in ports] + if self.source.port not in port_names: + errors.append( + f"Source port '{self.source.port}' not found. " + f"Available: {port_names}" + ) + + # Validate monitor port names + if ports and self.monitors: + port_names = [p.name for p in ports] + errors.extend( + f"Monitor port '{m}' not found. Available: {port_names}" + for m in self.monitors + if m not in port_names + ) + + if self.geometry.stack is None: + warnings_list.append( + "No stack configured. Will use active PDK with defaults." + ) + + return ValidationResult( + valid=len(errors) == 0, errors=errors, warnings=warnings_list + ) + + # ------------------------------------------------------------------------- + # Internal: stack resolution + # ------------------------------------------------------------------------- + + def _ensure_stack(self) -> None: + """Lazily resolve the layer stack if not yet built.""" + if self.geometry.stack is not None: + return + + from gsim.common.stack import get_stack + + if self._stack_kwargs: + yaml_path = self._stack_kwargs.pop("yaml_path", None) + self.geometry.stack = get_stack(yaml_path=yaml_path, **self._stack_kwargs) + self._stack_kwargs["yaml_path"] = yaml_path + else: + # Fall back to active PDK defaults + self.geometry.stack = get_stack() + + # ------------------------------------------------------------------------- + # Internal: z-crop + # ------------------------------------------------------------------------- + + def _apply_z_crop(self) -> None: + """Apply z-crop to the stack if geometry.z_crop is set. + + Only applies once per stack — after cropping, sets z_crop to None + to prevent double-cropping on subsequent calls. + """ + if self.geometry.z_crop is None: + return + + from gsim.common.stack.extractor import Layer, LayerStack + from gsim.meep.ports import _find_highest_n_layer + + stack = self.geometry.stack + if stack is None: + raise ValueError("No stack configured for z-crop.") + + # Find reference layer + ref: Layer | None = None + if self.geometry.z_crop == "auto": + ref, best_n = _find_highest_n_layer(stack) + if ref is None or best_n <= 1.5: + raise ValueError( + "Could not auto-detect core layer (no layer with n > 1.5). " + "Set geometry.z_crop to an explicit layer name." + ) + else: + layer_name = self.geometry.z_crop + if layer_name not in stack.layers: + raise ValueError( + f"Layer '{layer_name}' not found. " + f"Available: {list(stack.layers.keys())}" + ) + ref = stack.layers[layer_name] + + z_lo = ref.zmin - self.domain.margin_z_below + z_hi = ref.zmax + self.domain.margin_z_above + + # Filter and clip layers + cropped: dict[str, Layer] = {} + for name, layer in stack.layers.items(): + if layer.zmax <= z_lo or layer.zmin >= z_hi: + continue + new_zmin = max(layer.zmin, z_lo) + new_zmax = min(layer.zmax, z_hi) + cropped[name] = layer.model_copy( + update={ + "zmin": new_zmin, + "zmax": new_zmax, + "thickness": new_zmax - new_zmin, + } + ) + + # Crop dielectrics + cropped_dielectrics = [] + for diel in stack.dielectrics: + if diel["zmax"] <= z_lo or diel["zmin"] >= z_hi: + continue + cropped_dielectrics.append( + { + **diel, + "zmin": max(diel["zmin"], z_lo), + "zmax": min(diel["zmax"], z_hi), + } + ) + + self.geometry.stack = LayerStack( + pdk_name=stack.pdk_name, + units=stack.units, + layers=cropped, + materials=stack.materials, + dielectrics=cropped_dielectrics, + simulation=stack.simulation, + ) + # Mark as applied by clearing the z_crop setting + self.geometry.z_crop = None + + # ------------------------------------------------------------------------- + # Internal: translate to config objects + # ------------------------------------------------------------------------- + + def _wavelength_config(self) -> Any: + """Derive WavelengthConfig from source.""" + from gsim.meep.models.config import WavelengthConfig + + return WavelengthConfig( + wavelength=self.source.wavelength, + bandwidth=self.source.bandwidth, + num_freqs=self.source.num_freqs, + ) + + def _source_config(self) -> Any: + """Translate ModeSource → SourceConfig.""" + from gsim.meep.models.config import SourceConfig + + return SourceConfig( + bandwidth=None, + port=self.source.port, + ) + + def _stopping_config(self) -> Any: + """Translate FDTD stopping fields → StoppingConfig.""" + from gsim.meep.models.config import StoppingConfig + + s = self.solver + return StoppingConfig( + mode=s.stopping, + max_time=s.max_time, + threshold=s.stopping_threshold, + dft_min_run_time=s.stopping_min_time, + decay_component=s.stopping_component, + decay_dt=s.stopping_dt, + decay_monitor_port=s.stopping_monitor_port, + ) + + def _domain_config(self) -> Any: + """Translate Domain → DomainConfig.""" + from gsim.meep.models.config import DomainConfig + + return DomainConfig( + dpml=self.domain.pml, + margin_xy=self.domain.margin, + margin_z_above=self.domain.margin_z_above, + margin_z_below=self.domain.margin_z_below, + port_margin=self.domain.port_margin, + extend_ports=self.domain.extend_ports, + ) + + def _resolution_config(self) -> Any: + """Translate FDTD.resolution → ResolutionConfig.""" + from gsim.meep.models.config import ResolutionConfig + + return ResolutionConfig(pixels_per_um=self.solver.resolution) + + def _accuracy_config(self) -> Any: + """Translate FDTD accuracy fields → AccuracyConfig.""" + from gsim.meep.models.config import AccuracyConfig + + return AccuracyConfig( + eps_averaging=self.solver.subpixel, + subpixel_maxeval=self.solver.subpixel_maxeval, + subpixel_tol=self.solver.subpixel_tol, + simplify_tol=self.solver.simplify_tol, + ) + + def _diagnostics_config(self) -> Any: + """Translate FDTD diagnostic fields → DiagnosticsConfig.""" + from gsim.meep.models.config import DiagnosticsConfig + + return DiagnosticsConfig( + save_geometry=self.solver.save_geometry, + save_fields=self.solver.save_fields, + save_epsilon_raw=self.solver.save_epsilon_raw, + save_animation=self.solver.save_animation, + animation_interval=self.solver.animation_interval, + preview_only=self.solver.preview_only, + verbose_interval=self.solver.verbose_interval, + ) + + def _material_overrides(self) -> dict[str, Any]: + """Convert materials dict to MaterialProperties overrides.""" + from gsim.common.stack.materials import MaterialProperties + + overrides: dict[str, MaterialProperties] = {} + for name, val in self._resolved_materials().items(): + overrides[name] = MaterialProperties( + type="dielectric", + refractive_index=val.n, + extinction_coeff=val.k, + ) + return overrides + + # ------------------------------------------------------------------------- + # build_config — single source of truth + # ------------------------------------------------------------------------- + + def build_config(self) -> BuildResult: + """Build the complete simulation config (single source of truth). + + All computation — validation, stack resolution, z-crop, port + extension, material resolution, MPI estimation — happens here. + Both :meth:`write_config` and the viz methods consume this output. + + Returns: + BuildResult with SimConfig, extended component, and original. + + Raises: + ValueError: If config is invalid. + """ + from gsim.meep.materials import resolve_materials + from gsim.meep.models.config import LayerStackEntry, SimConfig, SymmetryEntry + from gsim.meep.ports import extract_port_info + + validation = self.validate_config() + if not validation.valid: + raise ValueError("Invalid configuration:\n" + "\n".join(validation.errors)) + + # Resolve stack + self._ensure_stack() + if self.geometry.stack is None: + raise ValueError("Stack resolution failed.") + if self.geometry.component is None: + raise ValueError("No geometry set.") + + # Apply z-crop if requested + self._apply_z_crop() + + import gdsfactory as gf + + original_component = self.geometry.component.copy() + stack = self.geometry.stack + + # Build config objects + domain_cfg = self._domain_config() + wl_cfg = self._wavelength_config() + source_cfg = self._source_config() + stopping_cfg = self._stopping_config() + resolution_cfg = self._resolution_config() + accuracy_cfg = self._accuracy_config() + diagnostics_cfg = self._diagnostics_config() + + # Compute port extension length + extend_length = domain_cfg.extend_ports + if extend_length == 0.0: + extend_length = domain_cfg.margin_xy + domain_cfg.dpml + + # Extend waveguide ports into PML region + original_bbox: list[float] | None = None + if extend_length > 0: + bbox = original_component.dbbox() + original_bbox = [bbox.left, bbox.bottom, bbox.right, bbox.top] + component = gf.components.extend_ports( + original_component, length=extend_length + ) + else: + component = original_component + + # Build layer stack entries + layer_stack_entries = [] + used_materials: set[str] = set() + for layer_name, layer in stack.layers.items(): + layer_stack_entries.append( + LayerStackEntry( + layer_name=layer_name, + gds_layer=list(layer.gds_layer), + zmin=layer.zmin, + zmax=layer.zmax, + material=layer.material, + sidewall_angle=layer.sidewall_angle, + ) + ) + used_materials.add(layer.material) + + # Build dielectric entries + dielectric_entries = [] + for diel in stack.dielectrics: + dielectric_entries.append( + { + "name": diel["name"], + "zmin": diel["zmin"], + "zmax": diel["zmax"], + "material": diel["material"], + } + ) + used_materials.add(diel["material"]) + + # Extract port info from original component + port_infos = extract_port_info( + original_component, stack, source_port=source_cfg.port + ) + + # Resolve materials + material_data = resolve_materials( + used_materials, overrides=self._material_overrides() + ) + + # Compute source fwidth + 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 = [ + SymmetryEntry(direction=s.direction, phase=s.phase) + for s in self.domain.symmetries + ] + if symmetry_entries: + import warnings + + warnings.warn( + "Symmetries are not yet used in production S-parameter runs " + "(only applied in preview-only mode).", + stacklevel=2, + ) + + # Build SimConfig + sim_config = SimConfig( + gds_filename="layout.gds", + component_bbox=original_bbox, + layer_stack=layer_stack_entries, + dielectrics=dielectric_entries, + ports=port_infos, + materials=material_data, + wavelength=wl_cfg, + source=source_for_config, + stopping=stopping_cfg, + resolution=resolution_cfg, + domain=domain_cfg, + accuracy=accuracy_cfg, + diagnostics=diagnostics_cfg, + verbose_interval=diagnostics_cfg.verbose_interval, + symmetries=symmetry_entries, + ) + + # Estimate MPI process count + dpml = domain_cfg.dpml + margin_xy = domain_cfg.margin_xy + if original_bbox is not None: + bbox_w = original_bbox[2] - original_bbox[0] + bbox_h = original_bbox[3] - original_bbox[1] + else: + bbox = component.dbbox() + bbox_w = bbox.right - bbox.left + bbox_h = bbox.top - bbox.bottom + z_min = min(e.zmin for e in layer_stack_entries) + z_max = max(e.zmax for e in layer_stack_entries) + cell_x = bbox_w + 2 * (margin_xy + dpml) + cell_y = bbox_h + 2 * (margin_xy + dpml) + cell_z = (z_max - z_min) + 2 * dpml + + _meep_np_est = estimate_meep_np( + cell_x, cell_y, cell_z, resolution_cfg.pixels_per_um + ) + # TODO: use _meep_np_est once Batch vCPU allocation is passed to the + # container (lscpu reports host cores, not container vCPUs → OOM). + meep_np = 2 + sim_config.meep_np = meep_np + logger.info( + "meep_np=%d (estimated %d, cell %.1f x %.1f x %.1f um, res %d)", + meep_np, + _meep_np_est, + cell_x, + cell_y, + cell_z, + resolution_cfg.pixels_per_um, + ) + + return BuildResult( + config=sim_config, + component=component, + original_component=original_component, + ) + + # ------------------------------------------------------------------------- + # write_config + # ------------------------------------------------------------------------- + + def write_config(self, output_dir: str | Path) -> Path: + """Serialize simulation config to output directory. + + Thin wrapper around :meth:`build_config` — writes GDS, JSON, and + the runner script. + + Args: + output_dir: Directory to write layout.gds, sim_config.json, run_meep.py. + + Returns: + Path to the output directory. + + Raises: + ValueError: If config is invalid. + """ + from gsim.meep.script import generate_meep_script + + output_dir = Path(output_dir) + output_dir.mkdir(parents=True, exist_ok=True) + + result = self.build_config() + + # Write extended component GDS + result.component.write_gds(output_dir / "layout.gds") + + # Write JSON config + result.config.to_json(output_dir / "sim_config.json") + + # Write runner script + script_path = output_dir / "run_meep.py" + script_content = generate_meep_script(config_filename="sim_config.json") + script_path.write_text(script_content) + + logger.info("Config written to %s", output_dir) + return output_dir + + # ------------------------------------------------------------------------- + # run + # ------------------------------------------------------------------------- + + def run(self, output_dir: str | Path | None = None, *, verbose: bool = True) -> Any: + """Run MEEP simulation on the cloud. + + Args: + output_dir: Directory to write config and download results. + If None, a temporary directory is created automatically. + verbose: Print progress info. + + Returns: + SParameterResult with parsed S-parameters. + """ + import tempfile + + from gsim.gcloud import run_simulation + from gsim.meep.models.results import SParameterResult + + if output_dir is None: + output_dir = Path(tempfile.mkdtemp(prefix="meep_")) + + output_dir = self.write_config(output_dir) + + results = run_simulation( + output_dir, + job_type="meep", + verbose=verbose, + ) + + csv_path = results.get("s_parameters.csv") + if csv_path is not None: + return SParameterResult.from_csv(csv_path) + + logger.warning("No s_parameters.csv found in results") + return SParameterResult() + + # ------------------------------------------------------------------------- + # Visualization + # ------------------------------------------------------------------------- + + def plot_2d(self, **kwargs: Any) -> Any: + """Plot 2D cross-sections of the geometry. + + Uses :meth:`build_config` so the plot shows exactly what meep + processes — including extended ports and PML boundaries. + + Accepts the same keyword arguments as :func:`gsim.meep.viz.plot_2d`. + """ + from gsim.meep.viz import plot_2d + + result = self.build_config() + + return plot_2d( + component=result.component, + stack=self.geometry.stack, + domain_config=result.config.domain, + source_port=result.config.source.port, + extend_ports_length=0, + port_data=result.config.ports, + component_bbox=result.config.component_bbox, + **kwargs, + ) + + def plot_3d(self, **kwargs: Any) -> Any: + """Plot 3D visualization of the geometry. + + Uses :meth:`build_config` so the plot shows exactly what meep + processes — including extended ports. + + Accepts the same keyword arguments as :func:`gsim.meep.viz.plot_3d`. + """ + from gsim.meep.viz import plot_3d + + result = self.build_config() + + return plot_3d( + component=result.component, + stack=self.geometry.stack, + domain_config=result.config.domain, + extend_ports_length=0, + **kwargs, + ) diff --git a/src/gsim/meep/viz.py b/src/gsim/meep/viz.py new file mode 100644 index 00000000..6e6115da --- /dev/null +++ b/src/gsim/meep/viz.py @@ -0,0 +1,304 @@ +"""Visualization helpers for MEEP simulations. + +Standalone functions that build geometry models and render 2D/3D plots. +Called directly by ``Simulation.plot_2d()`` / ``plot_3d()`` — no legacy +``MeepSim`` dependency. +""" + +from __future__ import annotations + +from typing import TYPE_CHECKING, Any + +import matplotlib.pyplot as plt + +if TYPE_CHECKING: + from gdsfactory.component import Component + + from gsim.common import LayerStack + from gsim.common.geometry_model import GeometryModel + from gsim.meep.models.config import DomainConfig + + +# --------------------------------------------------------------------------- +# Geometry model construction +# --------------------------------------------------------------------------- + + +def build_geometry_model( + component: Component, + stack: LayerStack | None, + domain_config: DomainConfig, + extend_ports_length: float | None = None, +) -> GeometryModel: + """Build a GeometryModel from a component + stack for visualization. + + Uses the gdsfactory LayerStack from the active PDK (not the gsim + LayerStack), because ``LayeredComponentBase`` needs gdsfactory's + ``LayerStack`` for polygon extraction via ``DerivedLayer.get_shapes()``. + + When ``domain_config.extend_ports`` is configured, the waveguide + ports are extended into the PML region so the visualization matches + what the runner will simulate. + + Args: + component: gdsfactory Component to visualize. + stack: gsim LayerStack (used for z-crop clipping). May be None. + domain_config: Domain configuration for port extension length. + extend_ports_length: Override port extension length. Pass ``0`` + when the component has already been extended by + :meth:`Simulation.build_config`. ``None`` (default) computes + the length from *domain_config* as before. + + Returns: + GeometryModel ready for visualization. + + Raises: + ValueError: If no active PDK with a layer_stack. + """ + import gdsfactory as gf + + from gsim.common.geometry_model import extract_geometry_model + from gsim.common.layered_component import LayeredComponentBase + + pdk = gf.get_active_pdk() + gf_layer_stack = pdk.layer_stack + if gf_layer_stack is None: + raise ValueError( + "Active PDK has no layer_stack. Activate a PDK with a layer stack first." + ) + + # Compute port extension length + if extend_ports_length is not None: + extend_length = extend_ports_length + else: + # Default: same logic as build_config + extend_length = domain_config.extend_ports + if extend_length == 0.0: + extend_length = domain_config.margin_xy + domain_config.dpml + + lc = LayeredComponentBase( + component=component, + layer_stack=gf_layer_stack, + extend_ports=extend_length, + ) + gm = extract_geometry_model(lc) + + # If the stack has been z-cropped, clip the geometry model to match + if stack is not None: + gm = crop_geometry_model(gm, stack) + + return gm + + +def crop_geometry_model( + gm: GeometryModel, + stack: LayerStack, +) -> GeometryModel: + """Clip a GeometryModel's prisms and bbox to the stack z-range. + + Args: + gm: Source geometry model. + stack: LayerStack whose z-extent defines the crop window. + + Returns: + New GeometryModel clipped to the stack z-range. + """ + from gsim.common.geometry_model import GeometryModel, Prism + + z_lo = min(layer.zmin for layer in stack.layers.values()) + z_hi = max(layer.zmax for layer in stack.layers.values()) + + cropped_prisms: dict[str, list[Prism]] = {} + for layer_name, prism_list in gm.prisms.items(): + clipped = [] + for p in prism_list: + if p.z_top <= z_lo or p.z_base >= z_hi: + continue + clipped.append( + Prism( + vertices=p.vertices, + z_base=max(p.z_base, z_lo), + z_top=min(p.z_top, z_hi), + layer_name=p.layer_name, + material=p.material, + sidewall_angle=p.sidewall_angle, + original_polygon=p.original_polygon, + ) + ) + if clipped: + cropped_prisms[layer_name] = clipped + + # Recompute bbox + old_min, old_max = gm.bbox + new_bbox = ( + (old_min[0], old_min[1], z_lo), + (old_max[0], old_max[1], z_hi), + ) + + return GeometryModel( + prisms=cropped_prisms, + bbox=new_bbox, + layer_bboxes=gm.layer_bboxes, + layer_mesh_orders=gm.layer_mesh_orders, + ) + + +# --------------------------------------------------------------------------- +# Overlay construction +# --------------------------------------------------------------------------- + + +def build_overlay( + geometry_model: GeometryModel, + component: Component, + stack: LayerStack | None, + domain_config: DomainConfig, + source_port: str | None = None, + port_data: list | None = None, + component_bbox: list[float] | tuple[float, ...] | None = None, +) -> Any: + """Build a SimOverlay from config, if stack is available. + + Args: + geometry_model: The geometry model (for bbox). + component: Component (for bbox and port extraction fallback). + stack: gsim LayerStack (needed for dielectrics). May be None. + domain_config: Domain configuration. + source_port: Source port name (or None for auto). + port_data: Pre-computed port data from :meth:`Simulation.build_config`. + When provided, skips port extraction (avoids duplicate work). + 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()``. + + Returns: + SimOverlay or None if stack isn't configured. + """ + from gsim.meep.overlay import build_sim_overlay + + if stack is None: + return None + + # Use pre-computed port data or extract from component + if port_data is None: + try: + from gsim.meep.ports import extract_port_info + + comp_copy = component.copy() + port_data = extract_port_info(comp_copy, stack, source_port=source_port) + except Exception: + port_data = [] + + # Use pre-computed component bbox or extract from component + if component_bbox is not None: + orig_bbox = ( + component_bbox[0], + component_bbox[1], + component_bbox[2], + component_bbox[3], + ) + else: + bbox = component.dbbox() + orig_bbox = (bbox.left, bbox.bottom, bbox.right, bbox.top) + + dielectrics = stack.dielectrics if stack else [] + return build_sim_overlay( + geometry_model, + domain_config, + port_data, + dielectrics=dielectrics, + component_bbox=orig_bbox, + ) + + +# --------------------------------------------------------------------------- +# Public plot functions +# --------------------------------------------------------------------------- + + +def plot_3d( + component: Component, + stack: LayerStack | None, + domain_config: DomainConfig, + backend: str = "open3d", + extend_ports_length: float | None = None, + **kwargs: Any, +) -> Any: + """Create interactive 3D visualization of MEEP geometry. + + Args: + component: gdsfactory Component to visualize. + stack: gsim LayerStack (may be None). + domain_config: Domain configuration. + backend: "open3d" (Jupyter/VS Code) or "pyvista" (desktop). + extend_ports_length: Override port extension length (pass 0 if + the component is already extended). + **kwargs: Extra args forwarded to the backend renderer. + + Returns: + Renderer-specific widget or plotter object. + """ + from gsim.common.viz import plot_prisms_3d, plot_prisms_3d_open3d + + gm = build_geometry_model( + component, stack, domain_config, extend_ports_length=extend_ports_length + ) + if backend == "pyvista": + return plot_prisms_3d(gm, **kwargs) + if backend == "open3d": + return plot_prisms_3d_open3d(gm, **kwargs) + raise ValueError(f"Unsupported backend: {backend}. Use 'open3d' or 'pyvista'") + + +def plot_2d( + component: Component, + stack: LayerStack | None, + domain_config: DomainConfig, + source_port: str | None = None, + x: float | str | None = None, + y: float | str | None = None, + z: float | str = "core", + ax: plt.Axes | None = None, + legend: bool = True, + slices: str = "z", + extend_ports_length: float | None = None, + port_data: list | None = None, + component_bbox: list[float] | tuple[float, ...] | None = None, +) -> plt.Axes | None: + """Plot 2D cross-sections of the MEEP geometry. + + Args: + component: gdsfactory Component to visualize. + stack: gsim LayerStack (may be None). + domain_config: Domain configuration. + source_port: Source port name (for overlay). + x: X-coordinate or layer name for slice plane. + y: Y-coordinate or layer name for slice plane. + z: Z-coordinate or layer name for slice plane. + ax: Axes to draw on. If None, a new figure is created. + legend: Whether to show the legend. + slices: Slice direction(s) — "x", "y", "z", or combinations. + extend_ports_length: Override port extension length (pass 0 if + the component is already extended). + port_data: Pre-computed port data (skips re-extraction). + component_bbox: Original component bbox ``[xmin, ymin, xmax, ymax]`` + (for correct cell boundary computation with extended ports). + + Returns: + ``plt.Axes`` when *ax* was provided, otherwise ``None``. + """ + from gsim.common.viz import plot_prism_slices + + gm = build_geometry_model( + component, stack, domain_config, extend_ports_length=extend_ports_length + ) + overlay = build_overlay( + gm, + component, + stack, + domain_config, + source_port, + port_data=port_data, + component_bbox=component_bbox, + ) + return plot_prism_slices(gm, x, y, z, ax, legend, slices, overlay=overlay) diff --git a/src/gsim/palace/models/ports.py b/src/gsim/palace/models/ports.py index a1ba047f..3b2c09ca 100644 --- a/src/gsim/palace/models/ports.py +++ b/src/gsim/palace/models/ports.py @@ -85,7 +85,7 @@ class CPWPortConfig(BaseModel): @property def effective_name(self) -> str: """Get the effective port name.""" - return self.name if self.name else f"cpw_{self.lower}" + return self.name or f"cpw_{self.lower}" class TerminalConfig(BaseModel): diff --git a/tests/meep/__init__.py b/tests/meep/__init__.py new file mode 100644 index 00000000..e69de29b diff --git a/tests/meep/test_meep_models.py b/tests/meep/test_meep_models.py new file mode 100644 index 00000000..a3ea3e98 --- /dev/null +++ b/tests/meep/test_meep_models.py @@ -0,0 +1,1135 @@ +"""Tests for MEEP config models, ports, materials, script generation, and overlays.""" + +from __future__ import annotations + +import ast +import json + +import pytest +from pydantic import ValidationError + +from gsim.meep import ( + DomainConfig, + ResolutionConfig, + SimConfig, + SourceConfig, + SParameterResult, + WavelengthConfig, +) +from gsim.meep.models.config import ( + LayerStackEntry, + MaterialData, + PortData, + StoppingConfig, + SymmetryEntry, +) + +# --------------------------------------------------------------------------- +# Config model tests +# --------------------------------------------------------------------------- + + +class TestWavelengthConfig: + """Test WavelengthConfig frequency/wavelength conversion.""" + + def test_defaults(self): + cfg = WavelengthConfig() + assert cfg.wavelength == 1.55 + assert cfg.bandwidth == 0.1 + assert cfg.num_freqs == 11 + + def test_fcen(self): + cfg = WavelengthConfig(wavelength=1.55) + assert abs(cfg.fcen - 1.0 / 1.55) < 1e-10 + + def test_df(self): + cfg = WavelengthConfig(wavelength=1.55, bandwidth=0.1) + wl_min = 1.55 - 0.05 + wl_max = 1.55 + 0.05 + expected_df = 1.0 / wl_min - 1.0 / wl_max + assert abs(cfg.df - expected_df) < 1e-10 + + def test_model_dump(self): + cfg = WavelengthConfig() + d = cfg.model_dump() + assert "wavelength" in d + assert "fcen" in d + assert "df" in d + assert "num_freqs" in d + # stopping and run_after_sources are no longer in fdtd dict + assert "stopping" not in d + assert "run_after_sources" not in d + + +class TestResolutionConfig: + """Test ResolutionConfig presets.""" + + def test_default(self): + cfg = ResolutionConfig.default() + assert cfg.pixels_per_um == 32 + + def test_coarse(self): + cfg = ResolutionConfig.coarse() + assert cfg.pixels_per_um == 16 + + def test_fine(self): + cfg = ResolutionConfig.fine() + assert cfg.pixels_per_um == 64 + + def test_custom(self): + cfg = ResolutionConfig(pixels_per_um=48) + assert cfg.pixels_per_um == 48 + + def test_model_dump(self): + cfg = ResolutionConfig() + d = cfg.model_dump() + assert d["pixels_per_um"] == 32 + + +class TestSimConfig: + """Test SimConfig serialization.""" + + def test_to_json(self, tmp_path): + wl_cfg = WavelengthConfig() + fwidth = SourceConfig().compute_fwidth(wl_cfg.fcen, wl_cfg.df) + source_cfg = SourceConfig(fwidth=fwidth) + stopping_cfg = StoppingConfig(mode="dft_decay", max_time=200.0) + cfg = SimConfig( + gds_filename="layout.gds", + layer_stack=[ # ty: ignore[invalid-argument-type] + { + "layer_name": "core", + "gds_layer": [1, 0], + "zmin": 0.0, + "zmax": 0.22, + "material": "si", + "sidewall_angle": 0.0, + } + ], + ports=[ # ty: ignore[invalid-argument-type] + { + "name": "o1", + "center": [0, 0, 0.11], + "orientation": 0, + "width": 0.5, + "normal_axis": 0, + "direction": "-", + "is_source": True, + } + ], + materials={"si": {"refractive_index": 3.47, "extinction_coeff": 0.0}}, # ty: ignore[invalid-argument-type] + wavelength=wl_cfg, + source=source_cfg, + stopping=stopping_cfg, + resolution=ResolutionConfig(), + ) + path = tmp_path / "config.json" + cfg.to_json(path) + assert path.exists() + + data = json.loads(path.read_text()) + assert "layer_stack" in data + assert "ports" in data + assert "materials" in data + assert data["gds_filename"] == "layout.gds" + assert data["materials"]["si"]["refractive_index"] == 3.47 + assert data["layer_stack"][0]["layer_name"] == "core" + # JSON keys use serialization aliases (fdtd, run_after_sources, decay_by) + assert "source" in data + assert data["source"]["fwidth"] > data["fdtd"]["df"] + assert "stopping" in data + assert data["stopping"]["mode"] == "dft_decay" + assert data["stopping"]["run_after_sources"] == 200.0 + + +class TestPortData: + """Test PortData model.""" + + def test_creation(self): + p = PortData( + name="o1", + center=[0.0, 0.0, 0.11], + orientation=0.0, + width=0.5, + normal_axis=0, + direction="-", + is_source=True, + ) + assert p.name == "o1" + assert p.is_source + + def test_model_dump(self): + p = PortData( + name="o1", + center=[0.0, 0.0, 0.11], + orientation=180.0, + width=0.5, + normal_axis=0, + direction="+", + ) + d = p.model_dump() + assert d["name"] == "o1" + assert d["direction"] == "+" + + +class TestLayerStackEntry: + """Test LayerStackEntry model.""" + + def test_creation(self): + entry = LayerStackEntry( + layer_name="core", + gds_layer=[1, 0], + zmin=0.0, + zmax=0.22, + material="si", + ) + assert entry.layer_name == "core" + assert entry.gds_layer == [1, 0] + assert entry.sidewall_angle == 0.0 + + def test_with_sidewall_angle(self): + entry = LayerStackEntry( + layer_name="core", + gds_layer=[1, 0], + zmin=0.0, + zmax=0.22, + material="si", + sidewall_angle=10.0, + ) + assert entry.sidewall_angle == 10.0 + + def test_model_dump(self): + entry = LayerStackEntry( + layer_name="clad", + gds_layer=[2, 0], + zmin=-0.5, + zmax=0.5, + material="SiO2", + ) + d = entry.model_dump() + assert d["layer_name"] == "clad" + assert d["gds_layer"] == [2, 0] + + +class TestMaterialData: + """Test MaterialData model.""" + + def test_creation(self): + m = MaterialData(refractive_index=3.47) + assert m.refractive_index == 3.47 + assert m.extinction_coeff == 0.0 + + def test_with_extinction(self): + m = MaterialData(refractive_index=3.47, extinction_coeff=0.01) + assert m.extinction_coeff == 0.01 + + +# --------------------------------------------------------------------------- +# Port extraction tests +# --------------------------------------------------------------------------- + + +class TestPortExtraction: + """Test port info extraction.""" + + def test_get_port_normal(self): + from gsim.meep.ports import get_port_normal + + assert get_port_normal(0) == (0, "-") + assert get_port_normal(90) == (1, "-") + assert get_port_normal(180) == (0, "+") + assert get_port_normal(270) == (1, "+") + + def test_get_port_normal_invalid(self): + from gsim.meep.ports import get_port_normal + + with pytest.raises(ValueError, match="Invalid port orientation"): + get_port_normal(45) + + +# --------------------------------------------------------------------------- +# Port z-center tests +# --------------------------------------------------------------------------- + + +class TestPortZCenter: + """Test photonic port z-center logic.""" + + def test_z_center_uses_highest_n(self): + """Port z-center should use layer with highest refractive index.""" + from gsim.common.stack.extractor import Layer, LayerStack + + stack = LayerStack( + layers={ + "clad": Layer( + name="clad", + gds_layer=(2, 0), + zmin=-1.0, + zmax=1.0, + thickness=2.0, + material="SiO2", + layer_type="dielectric", + ), + "core": Layer( + name="core", + gds_layer=(1, 0), + zmin=0.0, + zmax=0.22, + thickness=0.22, + material="si", + layer_type="dielectric", + ), + } + ) + + from gsim.meep.ports import _get_z_center + + z = _get_z_center(stack) + # Si (n=3.47) has highest index, so z-center should be midpoint of core + assert abs(z - 0.11) < 1e-6 + + def test_z_center_fallback(self): + """Falls back to midpoint of all layers when no optical data.""" + from gsim.common.stack.extractor import Layer, LayerStack + + stack = LayerStack( + layers={ + "unknown_layer": Layer( + name="unknown_layer", + gds_layer=(99, 0), + zmin=0.0, + zmax=1.0, + thickness=1.0, + material="unknown_material", + layer_type="dielectric", + ), + } + ) + + from gsim.meep.ports import _get_z_center + + z = _get_z_center(stack) + assert abs(z - 0.5) < 1e-6 + + +# --------------------------------------------------------------------------- +# Materials resolution tests +# --------------------------------------------------------------------------- + + +class TestMaterialsResolution: + """Test material resolution from common DB.""" + + def test_resolve_known_material(self): + from gsim.common.stack.materials import get_material_properties + + props = get_material_properties("silicon") + assert props is not None + assert props.refractive_index == 3.47 + + def test_resolve_sio2(self): + from gsim.common.stack.materials import get_material_properties + + props = get_material_properties("SiO2") + assert props is not None + assert props.refractive_index == 1.44 + + def test_optical_classmethod(self): + from gsim.common.stack.materials import MaterialProperties + + mat = MaterialProperties.optical(refractive_index=2.5) + assert mat.refractive_index == 2.5 + assert mat.extinction_coeff == 0.0 + assert mat.type == "dielectric" + + +# --------------------------------------------------------------------------- +# Script generation tests +# --------------------------------------------------------------------------- + + +class TestScriptGeneration: + """Test MEEP runner script generation.""" + + def test_generates_valid_python(self): + from gsim.meep.script import generate_meep_script + + script = generate_meep_script() + ast.parse(script) + + def test_contains_config_filename(self): + from gsim.meep.script import generate_meep_script + + script = generate_meep_script("my_config.json") + assert "my_config.json" in script + + def test_has_main_function(self): + from gsim.meep.script import generate_meep_script + + script = generate_meep_script() + assert "def main():" in script + assert 'if __name__ == "__main__":' in script + + def test_has_gds_loading(self): + from gsim.meep.script import generate_meep_script + + script = generate_meep_script() + assert "import gdsfactory" in script + assert "import_gds" in script + + def test_has_triangulation(self): + from gsim.meep.script import generate_meep_script + + script = generate_meep_script() + assert "triangulate_polygon_with_holes" in script + assert "Delaunay" in script + + +# --------------------------------------------------------------------------- +# SParameterResult tests +# --------------------------------------------------------------------------- + + +class TestSParameterResult: + """Test S-parameter result parsing.""" + + def test_from_csv(self, tmp_path): + csv_path = tmp_path / "s_parameters.csv" + csv_path.write_text( + "wavelength,S11_mag,S11_phase,S21_mag,S21_phase\n" + "1.500000,0.100000,-30.0000,0.900000,45.0000\n" + "1.550000,0.150000,-25.0000,0.850000,50.0000\n" + ) + + result = SParameterResult.from_csv(csv_path) + assert len(result.wavelengths) == 2 + assert "S11" in result.s_params + assert "S21" in result.s_params + assert len(result.s_params["S11"]) == 2 + assert abs(abs(result.s_params["S11"][0]) - 0.1) < 1e-6 + + def test_empty_result(self): + result = SParameterResult() + assert result.wavelengths == [] + assert result.s_params == {} + + +# --------------------------------------------------------------------------- +# Layer sidewall_angle test +# --------------------------------------------------------------------------- + + +class TestLayerSidewallAngle: + """Test sidewall_angle field on Layer model.""" + + def test_default_zero(self): + from gsim.common.stack.extractor import Layer + + layer = Layer( + name="test", + gds_layer=(1, 0), + zmin=0.0, + zmax=0.22, + thickness=0.22, + material="si", + layer_type="dielectric", + ) + assert layer.sidewall_angle == 0.0 + + def test_custom_angle(self): + from gsim.common.stack.extractor import Layer + + layer = Layer( + name="test", + gds_layer=(1, 0), + zmin=0.0, + zmax=0.22, + thickness=0.22, + material="si", + layer_type="dielectric", + sidewall_angle=10.0, + ) + assert layer.sidewall_angle == 10.0 + + def test_to_dict_includes_nonzero_angle(self): + from gsim.common.stack.extractor import Layer + + layer = Layer( + name="test", + gds_layer=(1, 0), + zmin=0.0, + zmax=0.22, + thickness=0.22, + material="si", + layer_type="dielectric", + sidewall_angle=5.0, + ) + d = layer.to_dict() + assert d["sidewall_angle"] == 5.0 + + def test_to_dict_omits_zero_angle(self): + from gsim.common.stack.extractor import Layer + + layer = Layer( + name="test", + gds_layer=(1, 0), + zmin=0.0, + zmax=0.22, + thickness=0.22, + material="si", + layer_type="dielectric", + ) + d = layer.to_dict() + assert "sidewall_angle" not in d + + +# --------------------------------------------------------------------------- +# Import test +# --------------------------------------------------------------------------- + + +class TestImportWithoutMeep: + """Verify module imports work without meep installed.""" + + def test_import_all_public_api(self): + from gsim.meep import ( + DomainConfig, + ResolutionConfig, + SimConfig, + Simulation, + SourceConfig, + SParameterResult, + WavelengthConfig, + ) + + assert WavelengthConfig is not None + assert DomainConfig is not None + assert ResolutionConfig is not None + assert SParameterResult is not None + assert SimConfig is not None + assert SourceConfig is not None + assert Simulation is not None + + +# --------------------------------------------------------------------------- +# DomainConfig tests +# --------------------------------------------------------------------------- + + +class TestDomainConfig: + """Test DomainConfig model.""" + + def test_defaults(self): + cfg = DomainConfig() + assert cfg.dpml == 1.0 + assert cfg.margin_xy == 0.5 + assert cfg.margin_z_above == 0.5 + assert cfg.margin_z_below == 0.5 + assert cfg.port_margin == 0.5 + assert cfg.extend_ports == 0.0 + + def test_custom(self): + cfg = DomainConfig( + dpml=0.5, margin_xy=0.2, margin_z_above=0.3, margin_z_below=0.4 + ) + assert cfg.dpml == 0.5 + assert cfg.margin_xy == 0.2 + assert cfg.margin_z_above == 0.3 + assert cfg.margin_z_below == 0.4 + + def test_extend_ports_custom(self): + cfg = DomainConfig(extend_ports=2.5) + assert cfg.extend_ports == 2.5 + + def test_extend_ports_serialization(self): + cfg = DomainConfig(extend_ports=3.0) + d = cfg.model_dump() + assert d["extend_ports"] == 3.0 + + def test_model_dump(self): + cfg = DomainConfig( + dpml=2.0, margin_xy=0.5, margin_z_above=1.0, margin_z_below=1.5 + ) + d = cfg.model_dump() + assert d["dpml"] == 2.0 + assert d["margin_xy"] == 0.5 + assert d["margin_z_above"] == 1.0 + assert d["margin_z_below"] == 1.5 + + +class TestSimConfigComponentBbox: + """Test SimConfig.component_bbox field.""" + + def test_default_none(self): + cfg = SimConfig() + assert cfg.component_bbox is None + + def test_with_bbox(self): + cfg = SimConfig(component_bbox=[-5.0, -2.0, 5.0, 2.0]) + assert cfg.component_bbox == [-5.0, -2.0, 5.0, 2.0] + + def test_json_roundtrip(self, tmp_path): + cfg = SimConfig(component_bbox=[-1.0, -0.5, 1.0, 0.5]) + path = tmp_path / "config.json" + cfg.to_json(path) + data = json.loads(path.read_text()) + assert data["component_bbox"] == [-1.0, -0.5, 1.0, 0.5] + + def test_json_roundtrip_none(self, tmp_path): + cfg = SimConfig() + path = tmp_path / "config.json" + cfg.to_json(path) + data = json.loads(path.read_text()) + assert data["component_bbox"] is None + + +class TestDomainInSimConfig: + """Test domain serialization in SimConfig JSON.""" + + def test_domain_in_sim_config_json(self, tmp_path): + cfg = SimConfig( + gds_filename="layout.gds", + layer_stack=[], + ports=[], + materials={}, + wavelength=WavelengthConfig(), + resolution=ResolutionConfig(), + domain=DomainConfig(dpml=0.5, margin_xy=0.2), + ) + path = tmp_path / "config.json" + cfg.to_json(path) + + data = json.loads(path.read_text()) + assert "domain" in data + assert data["domain"]["dpml"] == 0.5 + assert data["domain"]["margin_xy"] == 0.2 + + +# --------------------------------------------------------------------------- +# SymmetryEntry tests +# --------------------------------------------------------------------------- + + +class TestSymmetryEntry: + """Test SymmetryEntry model.""" + + def test_creation(self): + s = SymmetryEntry(direction="X", phase=-1) + assert s.direction == "X" + assert s.phase == -1 + + def test_default_phase(self): + s = SymmetryEntry(direction="Y") + assert s.phase == 1 + + def test_model_dump(self): + s = SymmetryEntry(direction="Z", phase=-1) + d = s.model_dump() + assert d == {"direction": "Z", "phase": -1} + + def test_invalid_direction(self): + with pytest.raises(ValidationError): + SymmetryEntry(direction="W") # ty: ignore[invalid-argument-type] + + def test_invalid_phase(self): + with pytest.raises(ValidationError): + SymmetryEntry(direction="X", phase=2) # ty: ignore[invalid-argument-type] + + +# --------------------------------------------------------------------------- +# StoppingConfig tests +# --------------------------------------------------------------------------- + + +class TestStoppingConfig: + """Test StoppingConfig model.""" + + def test_default_is_fixed(self): + cfg = StoppingConfig() + assert cfg.mode == "fixed" + assert cfg.max_time == 100.0 + assert cfg.decay_dt == 50.0 + assert cfg.decay_component == "Ey" + assert cfg.threshold == 1e-3 + assert cfg.decay_monitor_port is None + + def test_decay_mode(self): + cfg = StoppingConfig(mode="decay", threshold=1e-4, decay_dt=25.0) + assert cfg.mode == "decay" + assert cfg.threshold == 1e-4 + assert cfg.decay_dt == 25.0 + + def test_invalid_mode(self): + with pytest.raises(ValidationError): + StoppingConfig(mode="invalid") # ty: ignore[invalid-argument-type] + + def test_threshold_bounds(self): + with pytest.raises(ValidationError): + StoppingConfig(threshold=0) + with pytest.raises(ValidationError): + StoppingConfig(threshold=1.0) + + +class TestWavelengthConfigStopping: + """Test that WavelengthConfig no longer embeds StoppingConfig.""" + + def test_model_dump_excludes_stopping(self): + cfg = WavelengthConfig() + d = cfg.model_dump() + assert "stopping" not in d + assert "run_after_sources" not in d + + +# --------------------------------------------------------------------------- +# SourceConfig tests +# --------------------------------------------------------------------------- + + +class TestSourceConfig: + """Test SourceConfig model.""" + + def test_defaults(self): + cfg = SourceConfig() + assert cfg.bandwidth is None + assert cfg.port is None + + def test_auto_fwidth(self): + """Auto fwidth should be max(3*monitor_df, 0.2*fcen).""" + cfg = SourceConfig() + fcen = 1.0 / 1.55 + monitor_df = 0.042 # typical for 0.1um bandwidth at 1550nm + fwidth = cfg.compute_fwidth(fcen, monitor_df) + expected = max(3 * monitor_df, 0.2 * fcen) + assert abs(fwidth - expected) < 1e-10 + assert fwidth > monitor_df # always wider than monitor + + def test_explicit_bandwidth(self): + """Explicit wavelength bandwidth converted to frequency.""" + cfg = SourceConfig(bandwidth=0.3) + fcen = 1.0 / 1.55 + fwidth = cfg.compute_fwidth(fcen, 0.042) + # Should convert 0.3um wavelength bw to freq bw + wl_min = 1.55 - 0.15 + wl_max = 1.55 + 0.15 + expected = 1.0 / wl_min - 1.0 / wl_max + assert abs(fwidth - expected) < 1e-10 + + def test_model_dump(self): + cfg = SourceConfig(port="o1") + fcen = 1.0 / 1.55 + fwidth = cfg.compute_fwidth(fcen, 0.042) + cfg_with_fwidth = cfg.model_copy(update={"fwidth": fwidth}) + d = cfg_with_fwidth.model_dump() + assert "fwidth" in d + assert "port" in d + assert d["port"] == "o1" + assert d["fwidth"] > 0 + + def test_auto_fwidth_wider_than_monitor(self): + """Auto source fwidth should always be wider than monitor df.""" + cfg = SourceConfig() + fcen = 1.0 / 1.55 + wl = WavelengthConfig(wavelength=1.55, bandwidth=0.1) + fwidth = cfg.compute_fwidth(fcen, wl.df) + assert fwidth > wl.df + + +# --------------------------------------------------------------------------- +# Script symmetry tests +# --------------------------------------------------------------------------- + + +class TestScriptSymmetry: + """Test that the runner script includes symmetry support.""" + + def test_script_has_build_symmetries(self): + from gsim.meep.script import generate_meep_script + + script = generate_meep_script() + assert "build_symmetries" in script + + def test_script_has_mp_mirror(self): + from gsim.meep.script import generate_meep_script + + script = generate_meep_script() + assert "mp.Mirror" in script + + def test_script_has_split_chunks_evenly(self): + from gsim.meep.script import generate_meep_script + + script = generate_meep_script() + assert "split_chunks_evenly" in script + + +# --------------------------------------------------------------------------- +# Script decay tests +# --------------------------------------------------------------------------- + + +class TestScriptDecay: + """Test that the runner script includes decay stopping logic.""" + + def test_script_has_stop_when_fields_decayed(self): + from gsim.meep.script import generate_meep_script + + script = generate_meep_script() + assert "stop_when_fields_decayed" in script + + def test_script_has_component_map(self): + from gsim.meep.script import generate_meep_script + + script = generate_meep_script() + assert "_COMPONENT_MAP" in script + assert "mp.Ez" in script + assert "mp.Ey" in script + + def test_script_valid_python(self): + from gsim.meep.script import generate_meep_script + + script = generate_meep_script() + ast.parse(script) + + +# --------------------------------------------------------------------------- +# Overlay tests +# --------------------------------------------------------------------------- + + +class TestOverlay: + """Test SimOverlay + PortOverlay dataclasses and build_sim_overlay.""" + + def test_port_overlay_creation(self): + from gsim.meep.overlay import PortOverlay + + p = PortOverlay( + name="o1", + center=(0.0, 0.0, 0.11), + width=0.5, + normal_axis=0, + direction="-", + is_source=True, + z_span=0.22, + ) + assert p.name == "o1" + assert p.is_source + assert p.z_span == 0.22 + + def test_sim_overlay_creation(self): + from gsim.meep.overlay import SimOverlay + + ov = SimOverlay( + cell_min=(-5.0, -2.0, -1.0), + cell_max=(5.0, 2.0, 1.0), + dpml=1.0, + ports=[], + ) + assert ov.cell_min == (-5.0, -2.0, -1.0) + assert ov.dpml == 1.0 + + def test_dielectric_overlay_creation(self): + from gsim.meep.overlay import DielectricOverlay + + d = DielectricOverlay( + name="oxide", + material="SiO2", + zmin=-1.0, + zmax=0.0, + ) + assert d.name == "oxide" + assert d.material == "SiO2" + assert d.zmin == -1.0 + + def test_build_sim_overlay(self): + import numpy as np + + from gsim.common.geometry_model import GeometryModel, Prism + from gsim.meep.overlay import build_sim_overlay + + gm = GeometryModel( + prisms={ + "core": [ + Prism( + vertices=np.array([[-2, -1], [2, -1], [2, 1], [-2, 1]]), + z_base=0.0, + z_top=0.22, + layer_name="core", + ) + ] + }, + bbox=((-2.0, -1.0, 0.0), (2.0, 1.0, 0.22)), + ) + + domain_cfg = DomainConfig( + dpml=1.0, margin_xy=0.5, margin_z_above=0.0, margin_z_below=0.0 + ) + + port_data = [ + PortData( + name="o1", + center=[-2.0, 0.0, 0.11], + orientation=0.0, + width=0.5, + normal_axis=0, + direction="-", + is_source=True, + ), + PortData( + name="o2", + center=[2.0, 0.0, 0.11], + orientation=180.0, + width=0.5, + normal_axis=0, + direction="+", + is_source=False, + ), + ] + + overlay = build_sim_overlay(gm, domain_cfg, port_data) + + # cell_min = geo_min - (margin_xy + dpml) for xy, - dpml for z + assert overlay.cell_min[0] == pytest.approx(-3.5) # -2 - (0.5 + 1.0) + assert overlay.cell_min[1] == pytest.approx(-2.5) # -1 - (0.5 + 1.0) + assert overlay.cell_min[2] == pytest.approx(-1.0) # 0 - 1.0 + assert overlay.cell_max[0] == pytest.approx(3.5) + assert overlay.cell_max[1] == pytest.approx(2.5) + assert overlay.cell_max[2] == pytest.approx(1.22) # 0.22 + 1.0 + + assert len(overlay.ports) == 2 + assert overlay.ports[0].is_source + assert not overlay.ports[1].is_source + assert overlay.dpml == 1.0 + + def test_build_sim_overlay_with_dielectrics(self): + import numpy as np + + from gsim.common.geometry_model import GeometryModel, Prism + from gsim.meep.overlay import build_sim_overlay + + gm = GeometryModel( + prisms={ + "core": [ + Prism( + vertices=np.array([[-2, -1], [2, -1], [2, 1], [-2, 1]]), + z_base=0.0, + z_top=0.22, + layer_name="core", + ) + ] + }, + bbox=((-2.0, -1.0, -1.0), (2.0, 1.0, 1.22)), + ) + + domain_cfg = DomainConfig(dpml=1.0, margin_xy=0.5) + dielectrics = [ + {"name": "substrate", "material": "silicon", "zmin": -1.0, "zmax": 0.0}, + {"name": "oxide", "material": "SiO2", "zmin": 0.0, "zmax": 1.22}, + ] + + overlay = build_sim_overlay(gm, domain_cfg, [], dielectrics=dielectrics) + + assert len(overlay.dielectrics) == 2 + assert overlay.dielectrics[0].name == "substrate" + assert overlay.dielectrics[0].material == "silicon" + assert overlay.dielectrics[1].name == "oxide" + assert overlay.dielectrics[1].material == "SiO2" + + +# --------------------------------------------------------------------------- +# Script domain config tests +# --------------------------------------------------------------------------- + + +class TestScriptDomainConfig: + """Test that the runner script reads domain config.""" + + def test_script_reads_domain(self): + from gsim.meep.script import generate_meep_script + + script = generate_meep_script() + assert "domain" in script + assert "dpml" in script + assert "margin_xy" in script + + def test_script_uses_dpml_from_config(self): + """Verify the script no longer hardcodes padding = 1.0.""" + from gsim.meep.script import generate_meep_script + + script = generate_meep_script() + assert "padding = 1.0" not in script + + def test_script_has_background_slabs(self): + """Verify the script has build_background_slabs function.""" + from gsim.meep.script import generate_meep_script + + script = generate_meep_script() + assert "build_background_slabs" in script + assert "mp.Block" in script + assert "mp.inf" in script + + def test_script_handles_component_bbox(self): + """Verify the runner script uses component_bbox for cell sizing.""" + from gsim.meep.script import generate_meep_script + + script = generate_meep_script() + assert "component_bbox" in script + assert "bbox_left" in script + assert "bbox_right" in script + assert "bbox_top" in script + assert "bbox_bottom" in script + + +# --------------------------------------------------------------------------- +# Render2d overlay tests +# --------------------------------------------------------------------------- + + +class TestRender2dOverlay: + """Test overlay drawing in render2d.""" + + def test_plot_prism_slices_with_overlay(self): + """Verify plot_prism_slices accepts overlay kwarg without error.""" + import matplotlib as mpl + + mpl.use("Agg") + import matplotlib.pyplot as plt + import numpy as np + + from gsim.common.geometry_model import GeometryModel, Prism + from gsim.common.viz.render2d import plot_prism_slices + from gsim.meep.overlay import PortOverlay, SimOverlay + + gm = GeometryModel( + prisms={ + "core": [ + Prism( + vertices=np.array([[-2, -1], [2, -1], [2, 1], [-2, 1]]), + z_base=0.0, + z_top=0.22, + layer_name="core", + ) + ] + }, + bbox=((-2.0, -1.0, 0.0), (2.0, 1.0, 0.22)), + ) + + overlay = SimOverlay( + cell_min=(-3.0, -2.0, -1.0), + cell_max=(3.0, 2.0, 1.22), + dpml=1.0, + ports=[ + PortOverlay( + name="o1", + center=(-2.0, 0.0, 0.11), + width=0.5, + normal_axis=0, + direction="-", + is_source=True, + z_span=0.22, + ), + ], + ) + + fig, ax = plt.subplots() + result = plot_prism_slices(gm, z=0.11, ax=ax, overlay=overlay) + assert result is ax + + # Check that sim cell boundary was drawn (check patches) + patch_labels = [p.get_label() for p in ax.patches] + assert "Sim cell" in patch_labels + + plt.close(fig) + + def test_plot_without_overlay_fallback(self): + """Without overlay, should still draw geometry bbox.""" + import matplotlib as mpl + + mpl.use("Agg") + import matplotlib.pyplot as plt + import numpy as np + + from gsim.common.geometry_model import GeometryModel, Prism + from gsim.common.viz.render2d import plot_prism_slices + + gm = GeometryModel( + prisms={ + "core": [ + Prism( + vertices=np.array([[-2, -1], [2, -1], [2, 1], [-2, 1]]), + z_base=0.0, + z_top=0.22, + layer_name="core", + ) + ] + }, + bbox=((-2.0, -1.0, 0.0), (2.0, 1.0, 0.22)), + ) + + fig, ax = plt.subplots() + result = plot_prism_slices(gm, z=0.11, ax=ax) + assert result is ax + + patch_labels = [p.get_label() for p in ax.patches] + assert "Simulation" in patch_labels + + plt.close(fig) + + def test_plot_with_dielectric_overlay(self): + """Verify dielectric overlays render without error.""" + import matplotlib as mpl + + mpl.use("Agg") + import matplotlib.pyplot as plt + import numpy as np + + from gsim.common.geometry_model import GeometryModel, Prism + from gsim.common.viz.render2d import plot_prism_slices + from gsim.meep.overlay import DielectricOverlay, SimOverlay + + gm = GeometryModel( + prisms={ + "core": [ + Prism( + vertices=np.array([[-2, -1], [2, -1], [2, 1], [-2, 1]]), + z_base=0.0, + z_top=0.22, + layer_name="core", + ) + ] + }, + bbox=((-2.0, -1.0, -1.0), (2.0, 1.0, 1.22)), + ) + + overlay = SimOverlay( + cell_min=(-3.0, -2.0, -2.0), + cell_max=(3.0, 2.0, 2.22), + dpml=1.0, + ports=[], + dielectrics=[ + DielectricOverlay( + name="substrate", material="silicon", zmin=-1.0, zmax=0.0 + ), + DielectricOverlay(name="oxide", material="SiO2", zmin=0.0, zmax=1.22), + ], + ) + + # Test XZ view (side view — should draw horizontal bands) + fig, ax = plt.subplots() + result = plot_prism_slices(gm, y=0.0, ax=ax, slices="y", overlay=overlay) + assert result is ax + + patch_labels = [p.get_label() for p in ax.patches] + assert "substrate" in patch_labels + assert "oxide" in patch_labels + + plt.close(fig) + + # Test XY view (top view — should draw background fill for active z) + fig, ax = plt.subplots() + result = plot_prism_slices(gm, z=0.11, ax=ax, slices="z", overlay=overlay) + assert result is ax + + patch_labels = [p.get_label() for p in ax.patches] + # z=0.11 is within oxide (0.0 to 1.22) but not substrate (-1.0 to 0.0) + assert "oxide" in patch_labels + + plt.close(fig) diff --git a/tests/meep/test_simulation.py b/tests/meep/test_simulation.py new file mode 100644 index 00000000..b84009cc --- /dev/null +++ b/tests/meep/test_simulation.py @@ -0,0 +1,401 @@ +"""Tests for the declarative Simulation API.""" + +from __future__ import annotations + +import pytest +from pydantic import ValidationError + +from gsim.meep import ( + FDTD, + Domain, + Geometry, + Material, + ModeSource, + Simulation, + Symmetry, +) + +# --------------------------------------------------------------------------- +# Model construction & defaults +# --------------------------------------------------------------------------- + + +class TestMaterial: + def test_basic(self): + m = Material(n=3.47) + assert m.n == 3.47 + assert m.k == 0.0 + + def test_with_extinction(self): + m = Material(n=3.47, k=0.01) + assert m.k == 0.01 + + def test_n_must_be_positive(self): + with pytest.raises(ValidationError): + Material(n=0) + + def test_k_must_be_non_negative(self): + with pytest.raises(ValidationError): + Material(n=1.5, k=-0.1) + + +class TestModeSource: + def test_defaults(self): + s = ModeSource() + assert s.port is None + assert s.wavelength == 1.55 + assert s.bandwidth == 0.1 + assert s.num_freqs == 11 + + def test_custom(self): + s = ModeSource(port="o1", wavelength=1.31, bandwidth=0.05, num_freqs=21) + assert s.port == "o1" + assert s.wavelength == 1.31 + assert s.bandwidth == 0.05 + assert s.num_freqs == 21 + + +class TestDomain: + def test_defaults(self): + d = Domain() + assert d.pml == 1.0 + assert d.margin == 0.5 + assert d.margin_z_above == 0.5 + assert d.margin_z_below == 0.5 + assert d.port_margin == 0.5 + assert d.extend_ports == 0.0 + assert d.symmetries == [] + + def test_custom(self): + d = Domain(pml=0.5, margin=0.2, port_margin=0.3) + assert d.pml == 0.5 + assert d.margin == 0.2 + assert d.port_margin == 0.3 + + def test_symmetries(self): + d = Domain(symmetries=[Symmetry(direction="Y", phase=-1)]) + assert len(d.symmetries) == 1 + assert d.symmetries[0].direction == "Y" + assert d.symmetries[0].phase == -1 + + +class TestStoppingFields: + def test_defaults(self): + f = FDTD() + assert f.stopping == "dft_decay" + assert f.max_time == 200.0 + assert f.stopping_threshold == 1e-3 + assert f.stopping_min_time == 100.0 + assert f.stopping_component == "Ey" + assert f.stopping_dt == 50.0 + assert f.stopping_monitor_port is None + + def test_fixed_mode(self): + f = FDTD(stopping="fixed", max_time=100) + assert f.stopping == "fixed" + assert f.max_time == 100 + + def test_decay_mode(self): + f = FDTD(stopping="decay", stopping_threshold=1e-4, stopping_dt=25) + assert f.stopping == "decay" + assert f.stopping_threshold == 1e-4 + assert f.stopping_dt == 25 + + def test_dft_decay_custom(self): + f = FDTD(max_time=300, stopping_threshold=1e-4, stopping_min_time=80) + assert f.stopping == "dft_decay" + assert f.max_time == 300 + assert f.stopping_threshold == 1e-4 + assert f.stopping_min_time == 80 + + def test_invalid_mode(self): + with pytest.raises(ValidationError): + FDTD(stopping="invalid") # ty: ignore[invalid-argument-type] + + +class TestFDTD: + def test_defaults(self): + f = FDTD() + assert f.resolution == 32 + assert f.stopping == "dft_decay" + assert f.max_time == 200.0 + assert f.subpixel is False + assert f.simplify_tol == 0.0 + + def test_with_custom_stopping(self): + f = FDTD(stopping="decay", stopping_threshold=1e-4) + assert f.stopping == "decay" + assert f.stopping_threshold == 1e-4 + + def test_resolution_custom(self): + f = FDTD(resolution=64) + assert f.resolution == 64 + + +class TestFDTDDiagnostics: + def test_defaults(self): + f = FDTD() + assert f.save_geometry is True + assert f.save_fields is True + assert f.save_animation is False + assert f.preview_only is False + assert f.verbose_interval == 0 + + def test_custom(self): + f = FDTD(save_animation=True, verbose_interval=5.0, preview_only=True) + assert f.save_animation is True + assert f.verbose_interval == 5.0 + assert f.preview_only is True + + +# --------------------------------------------------------------------------- +# Material normalization +# --------------------------------------------------------------------------- + + +class TestMaterialNormalization: + def test_float_shorthand(self): + sim = Simulation(materials={"si": 3.47, "sio2": 1.44}) + assert isinstance(sim.materials["si"], Material) + assert sim.materials["si"].n == 3.47 + assert isinstance(sim.materials["sio2"], Material) + assert sim.materials["sio2"].n == 1.44 + + def test_material_object(self): + sim = Simulation(materials={"si": Material(n=3.47, k=0.01)}) + assert sim.materials["si"].n == 3.47 + assert sim.materials["si"].k == 0.01 + + def test_dict_shorthand(self): + sim = Simulation(materials={"si": {"n": 3.47, "k": 0.01}}) # ty: ignore[invalid-argument-type] + assert isinstance(sim.materials["si"], Material) + assert sim.materials["si"].n == 3.47 + + def test_assignment_normalization(self): + sim = Simulation() + sim.materials = {"si": 3.47} + assert isinstance(sim.materials["si"], Material) + assert sim.materials["si"].n == 3.47 + + +# --------------------------------------------------------------------------- +# Field-by-field assignment +# --------------------------------------------------------------------------- + + +class TestFieldAssignment: + def test_source_port(self): + sim = Simulation() + sim.source.port = "o1" + assert sim.source.port == "o1" + + def test_source_wavelength(self): + sim = Simulation() + sim.source.wavelength = 1.31 + assert sim.source.wavelength == 1.31 + + def test_domain_pml(self): + sim = Simulation() + sim.domain.pml = 0.5 + assert sim.domain.pml == 0.5 + + def test_solver_resolution(self): + sim = Simulation() + sim.solver.resolution = 64 + assert sim.solver.resolution == 64 + + def test_solver_stopping_replace(self): + sim = Simulation() + sim.solver.stopping = "decay" + sim.solver.stopping_threshold = 1e-4 + assert sim.solver.stopping == "decay" + assert sim.solver.stopping_threshold == 1e-4 + + def test_geometry_component(self): + sim = Simulation() + sim.geometry.component = "placeholder" + assert sim.geometry.component == "placeholder" + + +class TestWholeObjectAssignment: + def test_source(self): + sim = Simulation() + sim.source = ModeSource(port="o2", wavelength=1.31) + assert sim.source.port == "o2" + assert sim.source.wavelength == 1.31 + + def test_domain(self): + sim = Simulation() + sim.domain = Domain(pml=0.5, margin=0.2) + assert sim.domain.pml == 0.5 + assert sim.domain.margin == 0.2 + + def test_solver(self): + sim = Simulation() + sim.solver = FDTD(resolution=64, stopping="dft_decay") + assert sim.solver.resolution == 64 + assert sim.solver.stopping == "dft_decay" + + +# --------------------------------------------------------------------------- +# Validation +# --------------------------------------------------------------------------- + + +class TestValidation: + def test_missing_component(self): + sim = Simulation() + result = sim.validate_config() + assert not result.valid + assert any("No component" in e for e in result.errors) + + def test_no_stack_warning(self): + sim = Simulation() + result = sim.validate_config() + assert any("No stack" in w for w in result.warnings) + + def test_string_monitors(self): + sim = Simulation() + sim.monitors = ["o1", "o2", "o3"] + assert sim.monitors == ["o1", "o2", "o3"] + + def test_write_config_requires_output_dir(self): + sim = Simulation() + with pytest.raises(TypeError): + sim.write_config() + + +# --------------------------------------------------------------------------- +# Wavelength derivation +# --------------------------------------------------------------------------- + + +class TestWavelengthDerivation: + def test_from_source_defaults(self): + sim = Simulation() + wl = sim._wavelength_config() + assert wl.wavelength == 1.55 + assert wl.bandwidth == 0.1 + assert wl.num_freqs == 11 + + def test_from_source_custom(self): + sim = Simulation() + sim.source = ModeSource(wavelength=1.31, bandwidth=0.05, num_freqs=21) + wl = sim._wavelength_config() + assert wl.wavelength == 1.31 + assert wl.bandwidth == 0.05 + assert wl.num_freqs == 21 + + +# --------------------------------------------------------------------------- +# Config translation +# --------------------------------------------------------------------------- + + +class TestConfigTranslation: + def test_stopping_fixed(self): + sim = Simulation() + sim.solver.stopping = "fixed" + sim.solver.max_time = 200 + cfg = sim._stopping_config() + assert cfg.mode == "fixed" + assert cfg.max_time == 200 + + def test_stopping_decay(self): + sim = Simulation() + sim.solver.stopping = "decay" + sim.solver.stopping_threshold = 1e-4 + sim.solver.stopping_dt = 25 + sim.solver.stopping_monitor_port = "o2" + cfg = sim._stopping_config() + assert cfg.mode == "decay" + assert cfg.threshold == 1e-4 + assert cfg.decay_dt == 25 + assert cfg.decay_monitor_port == "o2" + + def test_stopping_dft_decay(self): + sim = Simulation() + sim.solver.stopping = "dft_decay" + sim.solver.max_time = 200 + sim.solver.stopping_threshold = 1e-4 + sim.solver.stopping_min_time = 80 + cfg = sim._stopping_config() + assert cfg.mode == "dft_decay" + assert cfg.max_time == 200 + assert cfg.threshold == 1e-4 + assert cfg.dft_min_run_time == 80 + + def test_domain_translation(self): + sim = Simulation() + sim.domain = Domain(pml=0.5, margin=0.3, margin_z_above=1.0, margin_z_below=0.8) + cfg = sim._domain_config() + assert cfg.dpml == 0.5 + assert cfg.margin_xy == 0.3 + assert cfg.margin_z_above == 1.0 + assert cfg.margin_z_below == 0.8 + + def test_resolution_translation(self): + sim = Simulation() + sim.solver.resolution = 64 + cfg = sim._resolution_config() + assert cfg.pixels_per_um == 64 + + def test_accuracy_translation(self): + sim = Simulation() + sim.solver.subpixel = True + sim.solver.simplify_tol = 0.01 + cfg = sim._accuracy_config() + assert cfg.eps_averaging is True + assert cfg.simplify_tol == 0.01 + + def test_source_translation(self): + sim = Simulation() + sim.source = ModeSource(port="o1", bandwidth=0.05) + cfg = sim._source_config() + assert cfg.port == "o1" + assert cfg.bandwidth is None # always auto fwidth + + def test_diagnostics_translation(self): + sim = Simulation() + sim.solver.save_animation = True + sim.solver.preview_only = True + cfg = sim._diagnostics_config() + assert cfg.save_animation is True + assert cfg.preview_only is True + + +# --------------------------------------------------------------------------- +# Import tests +# --------------------------------------------------------------------------- + + +class TestImports: + def test_import_all_new_api(self): + from gsim.meep import ( + FDTD, + Domain, + Material, + ModeSource, + Simulation, + Symmetry, + ) + + assert all( + cls is not None + for cls in [ + Domain, + FDTD, + Geometry, + Material, + ModeSource, + Simulation, + Symmetry, + ] + ) + + def test_simulation_instantiation(self): + sim = Simulation() + assert sim.geometry.component is None + assert sim.materials == {} + assert sim.monitors == [] + assert sim.solver.resolution == 32