diff --git a/README.md b/README.md index 74837b64..ed4edd5c 100644 --- a/README.md +++ b/README.md @@ -14,6 +14,7 @@ with minimal boilerplate. | Module | Solver | Method | Use Case | | ------------- | ------------------------------------------- | ------ | ------------------------------------------------------ | +| `gsim.fdtd` | GDSFactory FDTD | FDTD | PDK-native mesh and runtime configuration generation | | `gsim.palace` | [Palace](https://awslabs.github.io/palace/) | FEM | RF/microwave, impedance extraction, driven simulations | | `gsim.meep` | [Meep](https://meep.readthedocs.io/) | FDTD | Photonic components, S-parameters, mode propagation | diff --git a/docs/api/fdtd.md b/docs/api/fdtd.md new file mode 100644 index 00000000..9923f582 --- /dev/null +++ b/docs/api/fdtd.md @@ -0,0 +1,75 @@ +# GDSFactory FDTD API + +`gsim.fdtd` generates the coarse tetrahedral mesh and validated `config.json` +consumed by GDSFactory FDTD. The backend voxelizes this mesh onto its own Yee +grid, so the Gmsh mesh does not need to resolve the electromagnetic fields. + +## PDK-native workflow + +Pass the PDK module when it exposes project-level `MATERIAL_CARDS`; otherwise, +pass a PDK object or use the active PDK. Material names are resolved exactly, +using the project's cards first and gsim's built-in cards as fallbacks. + +```python +import gpdk + +from gsim import fdtd + +simulation = fdtd.Simulation(pdk=gpdk) +simulation.geometry("mmi1x2") +artifacts = simulation.write("fdtd_output") + +print(artifacts.mesh_path) # fdtd_output/mesh.msh +print(artifacts.config_path) # fdtd_output/config.json +``` + +The generated mesh is ASCII Gmsh MSH 2.2 with linear tetrahedra for material +regions and linear triangles for `port_` groups. Geometry and wavelength +values in the artifacts are in nanometers. PML extrusion is left to GDSFactory +FDTD. + +## Initial geometry limits + +The backend supports disconnected polygons, polygon holes, axis-aligned guided +ports, and vertical or constant-angle sidewalls. Tapered layers use a small +number of midpoint-sampled prisms selected from the Yee-cell size, keeping the +lateral approximation error below one quarter cell while leaving field +resolution to the backend voxelizer. + +Vertical `vertical_te` and `vertical_tm` ports are free-space apertures rather +than material-owned eigenmode ports. By default a vertical port becomes a +plane/fiber monitor while the first guided port is excited. Select the vertical +port explicitly to generate a Gaussian-beam source: + +```python +simulation = fdtd.Simulation(pdk=gpdk, default_port="o2") +simulation.geometry("grating_coupler_elliptical") +simulation.write("fdtd_output/grating") +``` + +The aperture defaults to a square using the port width, top-facing `+z`, with a +beam waist equal to half the aperture width. Override these policies with +`vertical_port_axis`, `vertical_port_aperture_width_um`, and +`vertical_port_waist_radius_um`. + +The initial implementation rejects unsupported `bias`/`z_to_bias` profiles and +lossy material snapshots because config schema version 1 accepts only real +scalar refractive indices. + +## Reference + +::: gsim.fdtd.Simulation + options: + show_source: false + +::: gsim.fdtd.SimulationArtifacts + options: + show_source: false + +::: gsim.fdtd.MeshManifest + options: + show_source: false + +::: gsim.fdtd.FDTDConfig + options: + show_source: false diff --git a/docs/index.md b/docs/index.md index 31d97a65..25988ac7 100644 --- a/docs/index.md +++ b/docs/index.md @@ -12,6 +12,7 @@ with minimal boilerplate. | Module | Solver | Method | Use Case | | ------------- | ------------------------------------------- | ------ | ------------------------------------------------------ | +| `gsim.fdtd` | GDSFactory FDTD | FDTD | PDK-native mesh and runtime configuration generation | | `gsim.palace` | [Palace](https://awslabs.github.io/palace/) | FEM | RF/microwave, impedance extraction, driven simulations | | `gsim.meep` | [Meep](https://meep.readthedocs.io/) | FDTD | Photonic components, S-parameters, mode propagation | @@ -62,5 +63,5 @@ result = sim.run() ## API Reference -See the API docs for full details: [Palace](api/palace.md), [Meep](api/meep.md), [Common](api/common.md), -[Cloud](api/cloud.md). +See the API docs for full details: [GDSFactory FDTD](api/fdtd.md), [Palace](api/palace.md), [Meep](api/meep.md), +[Common](api/common.md), [Cloud](api/cloud.md). diff --git a/docs/zensical.toml b/docs/zensical.toml index b1c4a71c..5337891f 100644 --- a/docs/zensical.toml +++ b/docs/zensical.toml @@ -29,6 +29,7 @@ nav = [ ] }, { "Changelog" = "CHANGELOG.md" }, { "API Reference" = [ + { "GDSFactory FDTD" = "api/fdtd.md" }, { "Palace" = "api/palace.md" }, { "Meep" = "api/meep.md" }, { "Common" = "api/common.md" }, diff --git a/src/gsim/__init__.py b/src/gsim/__init__.py index 100085c8..1d312d49 100644 --- a/src/gsim/__init__.py +++ b/src/gsim/__init__.py @@ -6,16 +6,19 @@ Currently includes: - palace: Palace EM simulation API - meep: MEEP photonic FDTD simulation API + - fdtd: PDK-native GDSFactory FDTD artifact generation """ from __future__ import annotations +from gsim import fdtd as fdtd from gsim.gcloud import get_status, wait_for_results __version__ = "0.3.0" __all__ = [ "__version__", + "fdtd", "get_status", "wait_for_results", ] diff --git a/src/gsim/common/materials/__init__.py b/src/gsim/common/materials/__init__.py index f1a54c3d..a4be6767 100644 --- a/src/gsim/common/materials/__init__.py +++ b/src/gsim/common/materials/__init__.py @@ -8,6 +8,7 @@ ) from gsim.common.materials.si_li_293k import SI_LI_293K from gsim.common.materials.si_salzberg import SI_SALZBERG +from gsim.common.materials.sin_luke import SIN_LUKE from gsim.common.materials.sio2_malitson import SIO2_MALITSON from gsim.common.materials.snapshots import ( MaterialModelError, @@ -20,6 +21,7 @@ __all__ = [ "GSIM_MATERIAL_CARDS", + "SIN_LUKE", "SIO2_MALITSON", "SI_LI_293K", "SI_SALZBERG", diff --git a/src/gsim/common/materials/_helpers.py b/src/gsim/common/materials/_helpers.py index 8b1df180..7a158dbc 100644 --- a/src/gsim/common/materials/_helpers.py +++ b/src/gsim/common/materials/_helpers.py @@ -30,7 +30,7 @@ def wavelength_validity(minimum_um: float, maximum_um: float) -> Validity: def material_card( name: str, permittivity: Index | Sellmeier, - temperature_ref: float, + temperature_ref: float | None, ) -> MaterialCard: """Build a compact optical material card.""" provenance = Provenance( diff --git a/src/gsim/common/materials/registry.py b/src/gsim/common/materials/registry.py index 157965a4..48ced866 100644 --- a/src/gsim/common/materials/registry.py +++ b/src/gsim/common/materials/registry.py @@ -8,6 +8,7 @@ from gsim.common.materials.si_li_293k import SI_LI_293K from gsim.common.materials.si_salzberg import SI_SALZBERG +from gsim.common.materials.sin_luke import SIN_LUKE from gsim.common.materials.sio2_malitson import SIO2_MALITSON MaterialSource = Literal["project", "gsim"] @@ -16,6 +17,8 @@ "Si": SI_SALZBERG.model_copy(update={"name": "Si"}), "Si-Salzberg": SI_SALZBERG, "Si-Li-293K": SI_LI_293K, + "SiN": SIN_LUKE.model_copy(update={"name": "SiN"}), + "SiN-Luke": SIN_LUKE, "SiO2": SIO2_MALITSON.model_copy(update={"name": "SiO2"}), "SiO2-Malitson": SIO2_MALITSON, } diff --git a/src/gsim/common/materials/sin_luke.py b/src/gsim/common/materials/sin_luke.py new file mode 100644 index 00000000..5b72ae3c --- /dev/null +++ b/src/gsim/common/materials/sin_luke.py @@ -0,0 +1,22 @@ +"""Luke et al. silicon-nitride model.""" + +from pdk_schema import Sellmeier, SellmeierTerm + +from gsim.common.materials._helpers import material_card, wavelength_validity + +SIN_LUKE = material_card( + name="SiN-Luke", + temperature_ref=None, + permittivity=Sellmeier( + validity=wavelength_validity(0.310, 5.504), + variation=None, + conductivity=None, + terms=( + SellmeierTerm(b=3.0249, c_um=0.1353406), + SellmeierTerm(b=40314.0, c_um=1239.842), + ), + offset=0.0, + ), +) + +__all__ = ["SIN_LUKE"] diff --git a/src/gsim/common/pdk/models.py b/src/gsim/common/pdk/models.py index e23a972e..b6872158 100644 --- a/src/gsim/common/pdk/models.py +++ b/src/gsim/common/pdk/models.py @@ -61,8 +61,14 @@ class ResolvedPort: width: float orientation: float normal: tuple[int, int, int] - layer_key: str - material: str + port_type: str + layer_key: str | None + material: str | None + + @property + def is_vertical(self) -> bool: + """Return whether this is a free-space vertical optical port.""" + return self.port_type.startswith("vertical_") @dataclass(frozen=True) diff --git a/src/gsim/common/pdk/resolve.py b/src/gsim/common/pdk/resolve.py index cebc1464..221493f5 100644 --- a/src/gsim/common/pdk/resolve.py +++ b/src/gsim/common/pdk/resolve.py @@ -190,6 +190,17 @@ def _axis_aligned_orientation( return aligned, (round(cos(angle)), round(sin(angle)), 0) +def _port_orientation_and_normal( + port: Any, +) -> tuple[float, tuple[int, int, int]]: + """Resolve guided-port normals while preserving vertical-port semantics.""" + port_type = str(getattr(port, "port_type", "")) + if port_type.startswith("vertical_"): + orientation = 0.0 if port.orientation is None else float(port.orientation) + return orientation % 360.0, (0, 0, 1) + return _axis_aligned_orientation(port.name, port.orientation) + + def _port_layer(port: Any) -> tuple[int, int]: """Return a concrete GDS tuple for a component port.""" if isinstance(port.layer, int): @@ -227,15 +238,17 @@ def _resolved_ports( """Map every component port to one resolved physical layer.""" resolved: dict[str, ResolvedPort] = {} for port in component.ports: - orientation, normal = _axis_aligned_orientation(port.name, port.orientation) + orientation, normal = _port_orientation_and_normal(port) + port_type = str(getattr(port, "port_type", "")) candidates = _port_candidates(port, layers) - if not candidates: + is_vertical = port_type.startswith("vertical_") + if not candidates and not is_vertical: raise UnsupportedPortError( f"Port {port.name!r} on layer {_port_layer(port)} does not map to " "resolved LayerStack geometry." ) - layer = candidates[0] - z_lower, z_upper = layer.z_bounds + layer = candidates[0] if candidates else None + z_lower, z_upper = layer.z_bounds if layer is not None else (0.0, 0.0) resolved[port.name] = ResolvedPort( name=port.name, center=( @@ -246,8 +259,9 @@ def _resolved_ports( width=float(port.width), orientation=orientation, normal=normal, - layer_key=layer.key, - material=layer.material, + port_type=port_type, + layer_key=layer.key if layer is not None else None, + material=layer.material if layer is not None else None, ) return resolved @@ -297,7 +311,7 @@ def resolve_passive_pcell( raise LayerResolutionError( f"Could not evaluate derived layers: {error}" ) from error - layers = _resolved_layers(derived_component, layer_stack) + layers = _resolved_layers(resolved_component, layer_stack) project_cards = _resolve_project_cards(pdk_object, pdk_or_module) materials = {} for material_name in dict.fromkeys(layer.material for layer in layers.values()): diff --git a/src/gsim/common/polygon.py b/src/gsim/common/polygon.py index 93c754db..f6e37d94 100644 --- a/src/gsim/common/polygon.py +++ b/src/gsim/common/polygon.py @@ -56,25 +56,7 @@ def fuse_polygons( Merged Shapely Polygon or MultiPolygon """ source_layer = getattr(layer, "layer", layer) - derived_layer = getattr(layer, "derived_layer", None) - - # Simulation export materializes derived layers onto their GDS targets. - # Prefer those polygons when present: re-evaluating the source Boolean - # expression loses LayerLevel.background semantics (notably full-height - # grating teeth). Fall back to evaluating the source expression for the - # ordinary, unmaterialized visualization path. - layer_region = None - if derived_layer is not None: - target = tuple(derived_layer.layer) - layer_index = component.kcl.layer(*target) - target_region = component.kdb_cell.begin_shapes_rec(layer_index) - if not target_region.at_end(): - from kfactory import kdb - - layer_region = kdb.Region(target_region) - - if layer_region is None: - layer_region = source_layer.get_shapes(component) + layer_region = source_layer.get_shapes(component) shapely_polygons = [] for klayout_polygon in layer_region.each_merged(): diff --git a/src/gsim/fdtd/__init__.py b/src/gsim/fdtd/__init__.py new file mode 100644 index 00000000..45664f13 --- /dev/null +++ b/src/gsim/fdtd/__init__.py @@ -0,0 +1,21 @@ +"""Passive artifact generation for GDSFactory FDTD.""" + +from gsim.fdtd.config import FDTDConfig +from gsim.fdtd.models import ( + FDTDArtifactError, + FDTDConfigError, + FDTDGeometryError, + MeshManifest, + SimulationArtifacts, +) +from gsim.fdtd.simulation import Simulation + +__all__ = [ + "FDTDArtifactError", + "FDTDConfig", + "FDTDConfigError", + "FDTDGeometryError", + "MeshManifest", + "Simulation", + "SimulationArtifacts", +] diff --git a/src/gsim/fdtd/config.py b/src/gsim/fdtd/config.py new file mode 100644 index 00000000..464b2f2a --- /dev/null +++ b/src/gsim/fdtd/config.py @@ -0,0 +1,311 @@ +"""Validated GDSFactory FDTD schema-version-1 configuration models.""" + +from __future__ import annotations + +from collections.abc import Mapping +from typing import Literal + +from pydantic import BaseModel, ConfigDict, Field, model_validator + +from gsim.common.materials import MaterialSnapshot +from gsim.fdtd.models import FDTDConfigError, MeshManifest + + +class _StrictModel(BaseModel): + """Base model that rejects fields GDSFactory FDTD does not understand.""" + + model_config = ConfigDict(extra="forbid", allow_inf_nan=False) + + +class MaterialConfig(_StrictModel): + """Scalar real optical material supported by GDSFactory FDTD schema v1.""" + + refractive_index: float = Field(gt=0) + + +class RegionConfig(_StrictModel): + """Material assignment for a three-dimensional physical group.""" + + phys_group: int = Field(gt=0) + material: str = Field(min_length=1) + priority: int = Field(ge=0) + + +class PortConfig(_StrictModel): + """Layer assignment and outward normal for a port physical group.""" + + phys_group: int = Field(gt=0) + layer: str = Field(min_length=1) + normal: tuple[int, int, int] + + @model_validator(mode="after") + def validate_axis_aligned_normal(self) -> PortConfig: + """Require exactly one signed unit-axis component.""" + if sum(component != 0 for component in self.normal) != 1 or any( + component not in {-1, 0, 1} for component in self.normal + ): + raise ValueError("port normal must be one signed Cartesian unit axis") + return self + + +class GeometryConfig(_StrictModel): + """All mesh physical groups consumed by GDSFactory FDTD.""" + + volumes: dict[str, RegionConfig] = Field(min_length=1) + layers: dict[str, RegionConfig] = Field(min_length=1) + ports: dict[str, PortConfig] = Field(default_factory=dict) + + +Vector3 = tuple[float, float, float] +SignedAxis = Literal["+x", "-x", "+y", "-y", "+z", "-z"] + + +class GaussianBeamConfig(_StrictModel): + """Free-space Gaussian beam injected through an axis-aligned aperture.""" + + region_min: Vector3 + region_max: Vector3 + aperture_normal: SignedAxis + propagation_direction: Vector3 + e_polarization: Vector3 + focal_point: Vector3 + waist_radius: float = Field(gt=0) + refractive_index: float = Field(gt=0) + + +class FiberModeConfig(_StrictModel): + """Analytic Gaussian fiber profile used by a plane monitor.""" + + propagation_direction: Vector3 + e_polarization: Vector3 + focal_point: Vector3 + waist_radius: float = Field(gt=0) + refractive_index: float = Field(gt=0) + + +class PlaneMonitorConfig(_StrictModel): + """Flux and optional fiber-overlap monitor on an arbitrary plane.""" + + name: str = Field(min_length=1) + region_min: Vector3 + region_max: Vector3 + normal: SignedAxis + flux: bool = True + wavelengths: list[float] | None = None + fiber_mode: FiberModeConfig | None = None + + @model_validator(mode="after") + def validate_plane(self) -> PlaneMonitorConfig: + """Require ordered bounds and zero thickness along the normal axis.""" + if any( + lower > upper + for lower, upper in zip(self.region_min, self.region_max, strict=True) + ): + raise ValueError("monitor region_min cannot exceed region_max") + normal_axis = {"x": 0, "y": 1, "z": 2}[self.normal[-1]] + if self.region_min[normal_axis] != self.region_max[normal_axis]: + raise ValueError("monitor region must be planar along its normal axis") + if self.wavelengths is not None and any( + wavelength <= 0 for wavelength in self.wavelengths + ): + raise ValueError("monitor wavelengths must be positive") + return self + + +class ExcitationConfig(_StrictModel): + """Initial eigenmode pulse configuration.""" + + type: Literal["eigenmode", "gaussian_beam"] = "eigenmode" + waveform: Literal["pulse", "continuous_wave"] = "pulse" + center_wavelength: float = Field(gt=0) + wavelength_halfspan: float = Field(ge=0) + num_wavelengths: int = Field(ge=1) + amplitude: float = 1.0 + default_port: str | None = Field(default=None, min_length=1) + gaussian_beam: GaussianBeamConfig | None = None + + @model_validator(mode="after") + def validate_wavelength_span(self) -> ExcitationConfig: + """Keep the wavelength sweep positive.""" + if self.wavelength_halfspan >= self.center_wavelength: + raise ValueError("wavelength_halfspan must be smaller than the center") + if self.waveform == "continuous_wave" and self.num_wavelengths != 1: + raise ValueError("continuous_wave requires num_wavelengths=1") + if self.amplitude == 0: + raise ValueError("excitation amplitude cannot be zero") + if self.type == "eigenmode": + if self.default_port is None: + raise ValueError("eigenmode excitation requires default_port") + if self.gaussian_beam is not None: + raise ValueError("eigenmode excitation cannot include gaussian_beam") + else: + if self.gaussian_beam is None: + raise ValueError("gaussian_beam excitation requires its settings") + if self.default_port is not None: + raise ValueError("gaussian_beam excitation cannot use default_port") + return self + + +class GridConfig(_StrictModel): + """Yee-grid and PML settings.""" + + nanometers_per_cell: float = Field(gt=0) + pml_cells: int = Field(ge=0) + + +class RunConfig(_StrictModel): + """GDSFactory FDTD termination controls.""" + + max_timesteps: int | None = Field(default=None, gt=0) + energy_decay_fraction: float = Field(gt=0, lt=1) + max_wall_seconds: float = Field(ge=0) + + +class FDTDConfig(_StrictModel): + """Complete GDSFactory FDTD runtime configuration.""" + + schema_version: Literal[1] = 1 + mesh_file: Literal["mesh.msh"] = "mesh.msh" + length_scale_meters: float = Field(default=1e-9, ge=1e-9, le=1e-9) + background_refractive_index: float = Field(gt=0) + materials: dict[str, MaterialConfig] = Field(min_length=1) + geometry: GeometryConfig + excitation: ExcitationConfig + monitors: list[PlaneMonitorConfig] = Field(default_factory=list) + grid: GridConfig + run: RunConfig + + @model_validator(mode="after") + def validate_references(self) -> FDTDConfig: + """Require all material, layer, and port references to exist.""" + material_names = set(self.materials) + for group_name, region in { + **self.geometry.volumes, + **self.geometry.layers, + }.items(): + if region.material not in material_names: + raise ValueError( + f"geometry group {group_name!r} references unknown material " + f"{region.material!r}" + ) + layer_names = set(self.geometry.layers) + for port_name, port in self.geometry.ports.items(): + if port.layer not in layer_names: + raise ValueError( + f"port {port_name!r} references unknown layer {port.layer!r}" + ) + if ( + self.excitation.type == "eigenmode" + and self.excitation.default_port not in self.geometry.ports + ): + raise ValueError( + f"default_port {self.excitation.default_port!r} is not declared" + ) + monitor_names = [monitor.name for monitor in self.monitors] + if len(monitor_names) != len(set(monitor_names)): + raise ValueError("monitor names must be unique") + return self + + +def _material_config(snapshot: MaterialSnapshot) -> MaterialConfig: + """Convert one lossless scalar snapshot to the GDSFactory FDTD schema.""" + if snapshot.extinction_coefficient != 0: + raise FDTDConfigError( + f"Material {snapshot.material_name!r} has extinction coefficient " + f"{snapshot.extinction_coefficient}; GDSFactory FDTD schema v1 " + "supports only " + "lossless real refractive indices." + ) + return MaterialConfig(refractive_index=snapshot.refractive_index) + + +def build_fdtd_config( + manifest: MeshManifest, + material_snapshots: Mapping[str, MaterialSnapshot], + *, + background_material: str, + center_wavelength_nm: float, + wavelength_halfspan_nm: float, + num_wavelengths: int, + default_port: str | None, + nanometers_per_cell: float, + pml_cells: int, + max_timesteps: int | None, + energy_decay_fraction: float, + max_wall_seconds: float, + gaussian_beam: GaussianBeamConfig | None = None, + monitors: list[PlaneMonitorConfig] | None = None, +) -> FDTDConfig: + """Build and cross-validate a GDSFactory FDTD config from a mesh manifest.""" + if background_material not in material_snapshots: + raise FDTDConfigError( + f"Background material {background_material!r} has no snapshot." + ) + materials = { + name: _material_config(snapshot) + for name, snapshot in material_snapshots.items() + } + return FDTDConfig( + background_refractive_index=materials[background_material].refractive_index, + materials=materials, + geometry=GeometryConfig( + volumes={ + name: RegionConfig( + phys_group=group.physical_tag, + material=group.material, + priority=group.priority, + ) + for name, group in manifest.volumes.items() + }, + layers={ + name: RegionConfig( + phys_group=group.physical_tag, + material=group.material, + priority=group.priority, + ) + for name, group in manifest.layers.items() + }, + ports={ + name: PortConfig( + phys_group=group.physical_tag, + layer=group.layer, + normal=group.normal, + ) + for name, group in manifest.ports.items() + }, + ), + excitation=ExcitationConfig( + type="gaussian_beam" if gaussian_beam is not None else "eigenmode", + center_wavelength=center_wavelength_nm, + wavelength_halfspan=wavelength_halfspan_nm, + num_wavelengths=num_wavelengths, + default_port=default_port, + gaussian_beam=gaussian_beam, + ), + monitors=monitors or [], + grid=GridConfig( + nanometers_per_cell=nanometers_per_cell, + pml_cells=pml_cells, + ), + run=RunConfig( + max_timesteps=max_timesteps, + energy_decay_fraction=energy_decay_fraction, + max_wall_seconds=max_wall_seconds, + ), + ) + + +__all__ = [ + "ExcitationConfig", + "FDTDConfig", + "FiberModeConfig", + "GaussianBeamConfig", + "GeometryConfig", + "GridConfig", + "MaterialConfig", + "PlaneMonitorConfig", + "PortConfig", + "RegionConfig", + "RunConfig", + "build_fdtd_config", +] diff --git a/src/gsim/fdtd/mesh.py b/src/gsim/fdtd/mesh.py new file mode 100644 index 00000000..ef7d2ebf --- /dev/null +++ b/src/gsim/fdtd/mesh.py @@ -0,0 +1,321 @@ +"""Coarse Gmsh artifact generation for GDSFactory FDTD voxelization.""" + +from __future__ import annotations + +from collections.abc import Mapping +from pathlib import Path +from typing import Any + +import gmsh + +from gsim.common.pdk import ResolvedLayer, ResolvedPassivePcell, ResolvedPort +from gsim.fdtd.mesh_geometry import ( + GEOMETRY_TOLERANCE_NM, + UM_TO_NM, + add_layer_volumes, +) +from gsim.fdtd.mesh_validation import validate_mesh +from gsim.fdtd.models import ( + FDTDGeometryError, + MeshGroup, + MeshManifest, + PortMeshGroup, +) + +_GMSH_OPTIONS_CHANGED = ( + "General.Terminal", + "Mesh.Binary", + "Mesh.ElementOrder", + "Mesh.MeshSizeExtendFromBoundary", + "Mesh.MeshSizeFromCurvature", + "Mesh.MeshSizeMax", + "Mesh.MeshSizeMin", + "Mesh.MshFileVersion", + "Mesh.SaveAll", +) + + +def _snapshot_gmsh_options() -> dict[str, float]: + """Capture options that FDTD meshing temporarily changes.""" + return { + option_name: gmsh.option.getNumber(option_name) + for option_name in _GMSH_OPTIONS_CHANGED + } + + +def _restore_gmsh_options(option_values: Mapping[str, float]) -> None: + """Restore options owned by an existing caller Gmsh session.""" + for option_name, value in option_values.items(): + gmsh.option.setNumber(option_name, value) + + +def _priority_by_mesh_order( + layers: Mapping[str, ResolvedLayer], +) -> dict[str, int]: + """Invert lower-wins PDK mesh order into higher-wins GDSFactory FDTD priority.""" + unique_orders = sorted({layer.mesh_order for layer in layers.values()}) + order_priority = { + mesh_order: len(unique_orders) - index + for index, mesh_order in enumerate(unique_orders) + } + return {name: order_priority[layer.mesh_order] for name, layer in layers.items()} + + +def background_bounds_nm( + resolved: ResolvedPassivePcell, + background_material: str, + padding_um: float, +) -> tuple[float, float, float, float, float, float]: + """Build a port-aligned background box from PDK and component bounds.""" + lower, upper = resolved.bounds + port_axes = { + next(index for index, value in enumerate(port.normal) if value) + for port in resolved.ports.values() + if not port.is_vertical + } + x_padding = 0.0 if 0 in port_axes else padding_um + y_padding = 0.0 if 1 in port_axes else padding_um + + background_z_bounds = [] + for level in resolved.layer_stack.layers.values(): + if level.material != background_material or level.thickness == 0: + continue + level_zmax = float(level.zmin + level.thickness) + background_z_bounds.append( + (min(float(level.zmin), level_zmax), max(float(level.zmin), level_zmax)) + ) + if background_z_bounds: + z_lower = min(lower[2], *(bounds[0] for bounds in background_z_bounds)) + z_upper = max(upper[2], *(bounds[1] for bounds in background_z_bounds)) + else: + z_lower = lower[2] - padding_um + z_upper = upper[2] + padding_um + + return ( + (lower[0] - x_padding) * UM_TO_NM, + (lower[1] - y_padding) * UM_TO_NM, + z_lower * UM_TO_NM, + (upper[0] + x_padding) * UM_TO_NM, + (upper[1] + y_padding) * UM_TO_NM, + z_upper * UM_TO_NM, + ) + + +def _add_physical_group(dimension: int, tags: list[int], name: str) -> int: + """Create a named Gmsh physical group and return its actual tag.""" + if not tags: + raise FDTDGeometryError(f"Physical group {name!r} has no entities.") + physical_tag = gmsh.model.addPhysicalGroup(dimension, tags) + gmsh.model.setPhysicalName(dimension, physical_tag, name) + return physical_tag + + +def _port_surface_tags( + port: Any, + volume_tags: list[int], + claimed_surfaces: set[int], +) -> list[int]: + """Find the owning layer boundary face at an axis-aligned port plane.""" + normal_axis = next(index for index, value in enumerate(port.normal) if value) + if normal_axis not in {0, 1}: + raise FDTDGeometryError( + f"Guided port {port.name!r} must have an in-plane normal." + ) + target_nm = port.center[normal_axis] * UM_TO_NM + transverse_axis = 1 - normal_axis + transverse_center_nm = port.center[transverse_axis] * UM_TO_NM + candidates: list[int] = [] + boundary_bounds: list[tuple[int, tuple[float, ...]]] = [] + for volume_tag in volume_tags: + for dimension, surface_tag in gmsh.model.getBoundary( + [(3, volume_tag)], + combined=False, + oriented=False, + recursive=False, + ): + if dimension != 2 or surface_tag in claimed_surfaces: + continue + bounds = gmsh.model.getBoundingBox(2, surface_tag) + boundary_bounds.append((surface_tag, bounds)) + if ( + abs(bounds[normal_axis] - target_nm) > GEOMETRY_TOLERANCE_NM + or abs(bounds[normal_axis + 3] - target_nm) > GEOMETRY_TOLERANCE_NM + ): + continue + if ( + bounds[transverse_axis] - GEOMETRY_TOLERANCE_NM + <= transverse_center_nm + <= bounds[transverse_axis + 3] + GEOMETRY_TOLERANCE_NM + ): + candidates.append(surface_tag) + if not candidates: + nearest_bounds = sorted( + boundary_bounds, + key=lambda item: min( + abs(item[1][normal_axis] - target_nm), + abs(item[1][normal_axis + 3] - target_nm), + ), + )[:3] + raise FDTDGeometryError( + f"Port {port.name!r} does not coincide with a boundary face of " + f"layer {port.layer_key!r}; nearest boundary bounds are {nearest_bounds}." + ) + candidates = sorted(set(candidates)) + claimed_surfaces.update(candidates) + return candidates + + +def _validate_port_on_background_face( + port: Any, + background_bounds: tuple[float, float, float, float, float, float], +) -> None: + """Require each port plane to lie on the material-union AABB face.""" + axis = next(index for index, value in enumerate(port.normal) if value) + side = 0 if port.normal[axis] < 0 else 3 + background_face = background_bounds[axis + side] + port_coordinate = port.center[axis] * UM_TO_NM + if abs(background_face - port_coordinate) > GEOMETRY_TOLERANCE_NM: + raise FDTDGeometryError( + f"Port {port.name!r} is not on the background domain face required " + "for unambiguous GDSFactory FDTD port extrusion." + ) + + +def _guided_layer_key(port: ResolvedPort) -> str: + """Return the physical owner required by a guided port.""" + if port.layer_key is None: + raise FDTDGeometryError( + f"Guided port {port.name!r} has no owning physical layer." + ) + return port.layer_key + + +def generate_mesh( + resolved: ResolvedPassivePcell, + mesh_path: Path, + *, + background_material: str, + background_padding_um: float, + mesh_size_nm: float, + nanometers_per_cell: float, +) -> MeshManifest: + """Generate and validate a coarse GDSFactory FDTD tetrahedral mesh.""" + if background_padding_um <= 0: + raise FDTDGeometryError("background_padding_um must be positive.") + if mesh_size_nm <= 0: + raise FDTDGeometryError("mesh_size_nm must be positive.") + if nanometers_per_cell <= 0: + raise FDTDGeometryError("nanometers_per_cell must be positive.") + if "background" in resolved.layers: + raise FDTDGeometryError( + "Layer name 'background' is reserved by GDSFactory FDTD." + ) + + initialized_here = not bool(gmsh.isInitialized()) + caller_option_values = {} if initialized_here else _snapshot_gmsh_options() + if initialized_here: + gmsh.initialize() + else: + gmsh.clear() + try: + gmsh.option.setNumber("General.Terminal", 0) + gmsh.model.add("gsim_fdtd") + kernel = gmsh.model.occ + background_bounds = background_bounds_nm( + resolved, + background_material, + background_padding_um, + ) + guided_ports = { + name: port for name, port in resolved.ports.items() if not port.is_vertical + } + for port in guided_ports.values(): + _validate_port_on_background_face(port, background_bounds) + + xmin, ymin, zmin, xmax, ymax, zmax = background_bounds + background_tag = kernel.addBox( + xmin, + ymin, + zmin, + xmax - xmin, + ymax - ymin, + zmax - zmin, + ) + layer_volume_tags = { + name: add_layer_volumes( + kernel, + layer, + [port for port in guided_ports.values() if port.layer_key == name], + nanometers_per_cell=nanometers_per_cell, + ) + for name, layer in resolved.layers.items() + } + kernel.synchronize() + + background_physical_tag = _add_physical_group(3, [background_tag], "background") + priorities = _priority_by_mesh_order(resolved.layers) + layer_groups = { + name: MeshGroup( + name=name, + physical_tag=_add_physical_group(3, tags, name), + material=resolved.layers[name].material, + priority=priorities[name], + ) + for name, tags in layer_volume_tags.items() + } + claimed_surfaces: set[int] = set() + port_groups = {} + for name, port in guided_ports.items(): + layer_key = _guided_layer_key(port) + surface_tags = _port_surface_tags( + port, + layer_volume_tags[layer_key], + claimed_surfaces, + ) + physical_name = f"port_{name}" + port_groups[name] = PortMeshGroup( + name=name, + physical_name=physical_name, + physical_tag=_add_physical_group(2, surface_tags, physical_name), + layer=layer_key, + normal=port.normal, + ) + manifest = MeshManifest( + volumes={ + "background": MeshGroup( + name="background", + physical_tag=background_physical_tag, + material=background_material, + priority=0, + ) + }, + layers=layer_groups, + ports=port_groups, + ) + + gmsh.option.setNumber("Mesh.MshFileVersion", 2.2) + gmsh.option.setNumber("Mesh.Binary", 0) + gmsh.option.setNumber("Mesh.ElementOrder", 1) + gmsh.option.setNumber("Mesh.SaveAll", 0) + gmsh.option.setNumber("Mesh.MeshSizeMin", mesh_size_nm) + gmsh.option.setNumber("Mesh.MeshSizeMax", mesh_size_nm) + gmsh.option.setNumber("Mesh.MeshSizeFromCurvature", 0) + gmsh.option.setNumber("Mesh.MeshSizeExtendFromBoundary", 0) + gmsh.model.mesh.generate(3) + gmsh.write(str(mesh_path)) + except FDTDGeometryError: + raise + except Exception as error: + raise FDTDGeometryError(f"Gmsh mesh generation failed: {error}") from error + finally: + if initialized_here: + gmsh.finalize() + else: + gmsh.clear() + _restore_gmsh_options(caller_option_values) + + validate_mesh(mesh_path, manifest) + return manifest + + +__all__ = ["background_bounds_nm", "generate_mesh", "validate_mesh"] diff --git a/src/gsim/fdtd/mesh_geometry.py b/src/gsim/fdtd/mesh_geometry.py new file mode 100644 index 00000000..f75316eb --- /dev/null +++ b/src/gsim/fdtd/mesh_geometry.py @@ -0,0 +1,274 @@ +"""Robust OCC geometry construction for coarse GDSFactory FDTD meshes.""" + +from __future__ import annotations + +from collections import defaultdict +from collections.abc import Iterable +from math import ceil, radians, tan +from typing import Any + +from shapely.geometry import GeometryCollection, MultiPolygon, Polygon, box +from shapely.geometry.base import BaseGeometry +from shapely.ops import unary_union +from shapely.validation import explain_validity + +from gsim.common.pdk import ResolvedLayer, ResolvedPort +from gsim.fdtd.models import FDTDGeometryError +from gsim.palace.mesh.gmsh_utils import extrude_polygon + +UM_TO_NM = 1000.0 +GEOMETRY_TOLERANCE_NM = 1e-3 +DEFAULT_MAX_SIDEWALL_ERROR_FRACTION = 0.25 + + +def _polygon_sort_key(polygon: Polygon) -> tuple[float, ...]: + """Return a stable ordering key for disconnected polygon members.""" + centroid = polygon.centroid + return ( + round(centroid.x, 12), + round(centroid.y, 12), + round(polygon.area, 12), + *(round(value, 12) for value in polygon.bounds), + ) + + +def iter_polygons(geometry: BaseGeometry, *, layer_key: str) -> list[Polygon]: + """Return every valid polygonal member while retaining interior rings.""" + if isinstance(geometry, Polygon): + polygons = [geometry] + elif isinstance(geometry, MultiPolygon): + polygons = list(geometry.geoms) + elif isinstance(geometry, GeometryCollection): + polygons = [part for part in geometry.geoms if isinstance(part, Polygon)] + else: + polygons = [] + if not polygons: + raise FDTDGeometryError( + f"Layer {layer_key!r} has unsupported or empty " + f"{type(geometry).__name__} geometry." + ) + for polygon in polygons: + if polygon.is_empty or not polygon.is_valid or polygon.area <= 0: + raise FDTDGeometryError( + f"Layer {layer_key!r} has invalid polygon geometry: " + f"{explain_validity(polygon)}." + ) + return sorted(polygons, key=_polygon_sort_key) + + +def _sidewall_offset_um(layer: ResolvedLayer, normalized_z: float) -> float: + """Evaluate the PDK's linear lateral offset at normalized z.""" + return ( + (layer.width_to_z - normalized_z) + * abs(layer.thickness) + * tan(radians(layer.sidewall_angle)) + ) + + +def sidewall_slice_count( + layer: ResolvedLayer, + nanometers_per_cell: float, + *, + max_error_fraction: float = DEFAULT_MAX_SIDEWALL_ERROR_FRACTION, +) -> int: + """Choose midpoint slices that bound lateral error relative to a Yee cell.""" + if nanometers_per_cell <= 0: + raise ValueError("nanometers_per_cell must be positive.") + if not 0 < max_error_fraction <= 1: + raise ValueError("max_error_fraction must be in the interval (0, 1].") + total_displacement_nm = ( + abs(layer.thickness) * abs(tan(radians(layer.sidewall_angle))) * UM_TO_NM + ) + if total_displacement_nm == 0: + return 1 + maximum_error_nm = max_error_fraction * nanometers_per_cell + return max(1, ceil(total_displacement_nm / (2 * maximum_error_nm))) + + +def _port_groups( + ports: Iterable[ResolvedPort], +) -> dict[tuple[int, int, float], list[ResolvedPort]]: + """Group ports that share one axis-aligned end plane.""" + groups: dict[tuple[int, int, float], list[ResolvedPort]] = defaultdict(list) + for port in ports: + normal_axis = next(index for index, value in enumerate(port.normal) if value) + if normal_axis not in {0, 1}: + raise FDTDGeometryError( + f"Guided port {port.name!r} is not in the component plane." + ) + groups[ + (normal_axis, port.normal[normal_axis], port.center[normal_axis]) + ].append(port) + return groups + + +def _clip_to_port_plane( + geometry: BaseGeometry, + *, + normal_axis: int, + normal_sign: int, + target: float, + maximum_port_width: float, +) -> BaseGeometry: + """Clip geometry at one outward-facing port plane.""" + xmin, ymin, xmax, ymax = geometry.bounds + margin = max(xmax - xmin, ymax - ymin, maximum_port_width, 1.0) + if normal_axis == 0 and normal_sign < 0: + clip = box(target, ymin - margin, xmax + margin, ymax + margin) + elif normal_axis == 0: + clip = box(xmin - margin, ymin - margin, target, ymax + margin) + elif normal_sign < 0: + clip = box(xmin - margin, target, xmax + margin, ymax + margin) + else: + clip = box(xmin - margin, ymin - margin, xmax + margin, target) + return geometry.intersection(clip) + + +def _condition_profile_at_ports( + geometry: BaseGeometry, + ports: list[ResolvedPort], + offset_um: float, +) -> BaseGeometry: + """Create all same-plane port stubs together, then clip once per plane.""" + conditioned = geometry + extension_epsilon_um = GEOMETRY_TOLERANCE_NM / UM_TO_NM + for (normal_axis, normal_sign, target), grouped_ports in _port_groups( + ports + ).items(): + transverse_axis = 1 - normal_axis + extensions = [] + for port in grouped_ports: + half_width = port.width / 2 + offset_um + if half_width <= 0: + raise FDTDGeometryError( + f"Layer {port.layer_key!r} sidewall closes port {port.name!r}." + ) + transverse_lower = port.center[transverse_axis] - half_width + transverse_upper = port.center[transverse_axis] + half_width + inward = target - normal_sign * (abs(offset_um) + extension_epsilon_um) + if normal_axis == 0: + extension = box( + min(target, inward), + transverse_lower, + max(target, inward), + transverse_upper, + ) + else: + extension = box( + transverse_lower, + min(target, inward), + transverse_upper, + max(target, inward), + ) + extensions.append(extension) + conditioned = unary_union([conditioned, *extensions]) + conditioned = _clip_to_port_plane( + conditioned, + normal_axis=normal_axis, + normal_sign=normal_sign, + target=target, + maximum_port_width=max(port.width for port in grouped_ports), + ) + return conditioned + + +def _validate_layer(layer: ResolvedLayer) -> None: + """Reject fabrication profiles the initial mesh writer cannot evaluate.""" + if layer.bias not in (None, 0, 0.0) or layer.z_to_bias is not None: + raise FDTDGeometryError( + f"Layer {layer.key!r} uses bias or z_to_bias, which is not supported " + "by the initial GDSFactory FDTD mesh writer." + ) + if not 0 <= layer.width_to_z <= 1: + raise FDTDGeometryError( + f"Layer {layer.key!r} width_to_z must be between 0 and 1." + ) + if abs(layer.sidewall_angle) >= 80: + raise FDTDGeometryError( + f"Layer {layer.key!r} sidewall angle is too steep to mesh safely." + ) + + +def _scaled_ring( + coordinates: Iterable[tuple[float, ...]], +) -> list[tuple[float, float]]: + """Convert one closed Shapely ring from micrometers to nanometers.""" + points = [ + (float(point[0]) * UM_TO_NM, float(point[1]) * UM_TO_NM) + for point in coordinates + ] + if len(points) >= 2 and points[0] == points[-1]: + points.pop() + return points + + +def _add_polygon_prism( + kernel: Any, + polygon: Polygon, + *, + z_lower_um: float, + z_upper_um: float, + layer_key: str, +) -> int: + """Extrude one polygon while retaining every interior ring as a hole.""" + exterior = _scaled_ring(polygon.exterior.coords) + holes = [] + for interior in polygon.interiors: + points = _scaled_ring(interior.coords) + holes.append(([point[0] for point in points], [point[1] for point in points])) + volume_tag = extrude_polygon( + kernel, + [point[0] for point in exterior], + [point[1] for point in exterior], + z_lower_um * UM_TO_NM, + (z_upper_um - z_lower_um) * UM_TO_NM, + holes=holes, + ) + if volume_tag is None: + raise FDTDGeometryError(f"Could not extrude layer {layer_key!r}.") + return volume_tag + + +def add_layer_volumes( + kernel: Any, + layer: ResolvedLayer, + ports: list[ResolvedPort], + *, + nanometers_per_cell: float, +) -> list[int]: + """Create all unfused midpoint-slice prisms for one logical layer.""" + _validate_layer(layer) + slice_count = sidewall_slice_count(layer, nanometers_per_cell) + z_lower_um, z_upper_um = layer.z_bounds + volume_tags = [] + for slice_index in range(slice_count): + lower_fraction = slice_index / slice_count + upper_fraction = (slice_index + 1) / slice_count + midpoint_fraction = (lower_fraction + upper_fraction) / 2 + offset_um = _sidewall_offset_um(layer, midpoint_fraction) + profile = layer.geometry.buffer(offset_um, join_style=2) + profile = _condition_profile_at_ports(profile, ports, offset_um) + slice_z_lower = z_lower_um + lower_fraction * (z_upper_um - z_lower_um) + slice_z_upper = z_lower_um + upper_fraction * (z_upper_um - z_lower_um) + volume_tags.extend( + _add_polygon_prism( + kernel, + polygon, + z_lower_um=slice_z_lower, + z_upper_um=slice_z_upper, + layer_key=layer.key, + ) + for polygon in iter_polygons(profile, layer_key=layer.key) + ) + if not volume_tags: + raise FDTDGeometryError(f"Layer {layer.key!r} produced no volumes.") + return volume_tags + + +__all__ = [ + "GEOMETRY_TOLERANCE_NM", + "UM_TO_NM", + "add_layer_volumes", + "iter_polygons", + "sidewall_slice_count", +] diff --git a/src/gsim/fdtd/mesh_validation.py b/src/gsim/fdtd/mesh_validation.py new file mode 100644 index 00000000..61aa2506 --- /dev/null +++ b/src/gsim/fdtd/mesh_validation.py @@ -0,0 +1,73 @@ +"""Strict validation for GDSFactory FDTD-compatible Gmsh artifacts.""" + +from __future__ import annotations + +from pathlib import Path + +import meshio + +from gsim.fdtd.models import FDTDGeometryError, MeshManifest + + +def _manifest_names(manifest: MeshManifest) -> set[str]: + """Return all physical names expected in a mesh.""" + return { + *manifest.volumes, + *manifest.layers, + *(port.physical_name for port in manifest.ports.values()), + } + + +def validate_mesh(mesh_path: Path, manifest: MeshManifest) -> None: + """Preflight the strict MSH 2.2 groups and elements GDSFactory FDTD reads.""" + lines = mesh_path.read_text(encoding="utf8").splitlines() + try: + mesh_format_index = lines.index("$MeshFormat") + nodes_index = lines.index("$Nodes") + except ValueError as error: + raise FDTDGeometryError("Mesh is missing MSH 2.2 sections.") from error + if lines[mesh_format_index + 1].strip() != "2.2 0 8": + raise FDTDGeometryError("GDSFactory FDTD requires ASCII Gmsh MSH 2.2 output.") + node_count = int(lines[nodes_index + 1]) + node_ids = [ + int(lines[nodes_index + 2 + index].split()[0]) for index in range(node_count) + ] + if node_ids != list(range(1, node_count + 1)): + raise FDTDGeometryError( + "GDSFactory FDTD requires sequential one-based node IDs." + ) + + mesh = meshio.read(mesh_path) + actual_names = set(mesh.field_data) + expected_names = _manifest_names(manifest) + if actual_names != expected_names: + raise FDTDGeometryError( + f"Mesh physical names {sorted(actual_names)} do not match manifest " + f"{sorted(expected_names)}." + ) + physical_data = mesh.cell_data_dict.get("gmsh:physical", {}) + for group in [*manifest.volumes.values(), *manifest.layers.values()]: + field_tag, dimension = mesh.field_data[group.name] + tetra_tags = physical_data.get("tetra", []) + if ( + dimension != 3 + or field_tag != group.physical_tag + or not any(tag == field_tag for tag in tetra_tags) + ): + raise FDTDGeometryError( + f"Volume group {group.name!r} has no linear tetrahedra." + ) + for port in manifest.ports.values(): + field_tag, dimension = mesh.field_data[port.physical_name] + triangle_tags = physical_data.get("triangle", []) + if ( + dimension != 2 + or field_tag != port.physical_tag + or not any(tag == field_tag for tag in triangle_tags) + ): + raise FDTDGeometryError( + f"Port group {port.physical_name!r} has no linear triangles." + ) + + +__all__ = ["validate_mesh"] diff --git a/src/gsim/fdtd/models.py b/src/gsim/fdtd/models.py new file mode 100644 index 00000000..ceb3a36f --- /dev/null +++ b/src/gsim/fdtd/models.py @@ -0,0 +1,68 @@ +"""Shared models and errors for GDSFactory FDTD artifact generation.""" + +from __future__ import annotations + +from dataclasses import dataclass +from pathlib import Path + + +class FDTDArtifactError(ValueError): + """Base error for invalid or unsupported GDSFactory FDTD artifacts.""" + + +class FDTDGeometryError(FDTDArtifactError): + """Raised when geometry cannot produce a valid GDSFactory FDTD mesh.""" + + +class FDTDConfigError(FDTDArtifactError): + """Raised when settings cannot produce a valid GDSFactory FDTD config.""" + + +@dataclass(frozen=True) +class MeshGroup: + """One named three-dimensional Gmsh physical group.""" + + name: str + physical_tag: int + material: str + priority: int + + +@dataclass(frozen=True) +class PortMeshGroup: + """One named two-dimensional Gmsh port physical group.""" + + name: str + physical_name: str + physical_tag: int + layer: str + normal: tuple[int, int, int] + + +@dataclass(frozen=True) +class MeshManifest: + """Authoritative physical-group mapping emitted with a Gmsh mesh.""" + + volumes: dict[str, MeshGroup] + layers: dict[str, MeshGroup] + ports: dict[str, PortMeshGroup] + + +@dataclass(frozen=True) +class SimulationArtifacts: + """Paths and metadata produced by :meth:`Simulation.write`.""" + + mesh_path: Path + config_path: Path + manifest: MeshManifest + + +__all__ = [ + "FDTDArtifactError", + "FDTDConfigError", + "FDTDGeometryError", + "MeshGroup", + "MeshManifest", + "PortMeshGroup", + "SimulationArtifacts", +] diff --git a/src/gsim/fdtd/simulation.py b/src/gsim/fdtd/simulation.py new file mode 100644 index 00000000..f02867e4 --- /dev/null +++ b/src/gsim/fdtd/simulation.py @@ -0,0 +1,328 @@ +"""Public Simulation workflow for GDSFactory FDTD artifact generation.""" + +from __future__ import annotations + +from collections.abc import Mapping +from dataclasses import dataclass, field +from math import cos, isfinite, radians, sin +from pathlib import Path +from typing import Any, Literal + +from pydantic import ValidationError + +from gsim.common.materials import ( + MaterialResolutionError, + MaterialSnapshot, + get_project_material_cards, + resolve_material_snapshot, +) +from gsim.common.pdk import ResolvedPassivePcell, resolve_passive_pcell +from gsim.fdtd.config import ( + FDTDConfig, + FiberModeConfig, + GaussianBeamConfig, + PlaneMonitorConfig, + build_fdtd_config, +) +from gsim.fdtd.mesh import background_bounds_nm, generate_mesh +from gsim.fdtd.models import ( + FDTDConfigError, + FDTDGeometryError, + MeshManifest, + SimulationArtifacts, +) + + +@dataclass +class Simulation: + """Generate coarse Gmsh and config artifacts for GDSFactory FDTD runs.""" + + pdk: Any | None = None + wavelength_um: float = 1.55 + background_material: str = "SiO2" + nanometers_per_cell: float = 31.25 + pml_cells: int = 32 + wavelength_halfspan_um: float = 0.05 + num_wavelengths: int = 11 + default_port: str | None = None + background_padding_um: float = 1.0 + mesh_size_nm: float = 500.0 + vertical_port_axis: Literal["+z", "-z"] = "+z" + vertical_port_aperture_width_um: float | None = None + vertical_port_waist_radius_um: float | None = None + max_timesteps: int | None = None + energy_decay_fraction: float = 1e-6 + max_wall_seconds: float = 3600.0 + _resolved: ResolvedPassivePcell | None = field( + default=None, + init=False, + repr=False, + ) + + def __post_init__(self) -> None: + """Validate constructor settings that do not depend on geometry.""" + positive_values = { + "wavelength_um": self.wavelength_um, + "nanometers_per_cell": self.nanometers_per_cell, + "background_padding_um": self.background_padding_um, + "mesh_size_nm": self.mesh_size_nm, + } + for name, value in positive_values.items(): + if not isfinite(value) or value <= 0: + raise ValueError(f"{name} must be finite and positive.") + if not self.background_material: + raise ValueError("background_material cannot be empty.") + if self.pml_cells < 0: + raise ValueError("pml_cells cannot be negative.") + if not 0 <= self.wavelength_halfspan_um < self.wavelength_um: + raise ValueError( + "wavelength_halfspan_um must be nonnegative and smaller than " + "wavelength_um." + ) + if self.num_wavelengths < 1: + raise ValueError("num_wavelengths must be at least 1.") + if self.max_timesteps is not None and self.max_timesteps <= 0: + raise ValueError("max_timesteps must be positive when provided.") + if not 0 < self.energy_decay_fraction < 1: + raise ValueError("energy_decay_fraction must be between 0 and 1.") + if self.max_wall_seconds < 0: + raise ValueError("max_wall_seconds cannot be negative.") + if self.vertical_port_axis not in {"+z", "-z"}: + raise ValueError("vertical_port_axis must be '+z' or '-z'.") + for name, value in { + "vertical_port_aperture_width_um": self.vertical_port_aperture_width_um, + "vertical_port_waist_radius_um": self.vertical_port_waist_radius_um, + }.items(): + if value is not None and (not isfinite(value) or value <= 0): + raise ValueError(f"{name} must be finite and positive when provided.") + + @property + def resolved(self) -> ResolvedPassivePcell: + """Return stored canonical geometry or fail before geometry setup.""" + if self._resolved is None: + raise FDTDGeometryError( + "No geometry is configured. Call Simulation.geometry(...) first." + ) + return self._resolved + + def geometry( + self, + component: Any, + *, + settings: Mapping[str, Any] | None = None, + ) -> ResolvedPassivePcell: + """Resolve and store a component through the canonical PDK boundary.""" + self._resolved = resolve_passive_pcell( + component, + pdk=self.pdk, + settings=settings, + wavelength_um=self.wavelength_um, + ) + return self._resolved + + def _material_snapshots(self) -> dict[str, MaterialSnapshot]: + """Add a strict project-first background snapshot to layer snapshots.""" + snapshots = dict(self.resolved.materials) + if self.background_material in snapshots: + return snapshots + try: + project_cards = get_project_material_cards(self.pdk) + snapshots[self.background_material] = resolve_material_snapshot( + self.background_material, + self.wavelength_um, + project_cards, + ) + except MaterialResolutionError as error: + raise FDTDConfigError( + f"Could not resolve background material " + f"{self.background_material!r}: {error}" + ) from error + return snapshots + + def _config( + self, + manifest: MeshManifest, + material_snapshots: Mapping[str, MaterialSnapshot], + ) -> FDTDConfig: + """Build the GDSFactory FDTD schema after mesh group tags are known.""" + if not self.resolved.ports: + raise FDTDConfigError( + "GDSFactory FDTD simulations require at least one port." + ) + selected_port_name = self._selected_port_name() + selected_port = self.resolved.ports[selected_port_name] + vertical_configs = self._vertical_port_configs(material_snapshots) + gaussian_beam = ( + vertical_configs[selected_port_name][0] + if selected_port.is_vertical + else None + ) + default_port = None if selected_port.is_vertical else selected_port_name + monitors = [monitor for _, monitor in vertical_configs.values()] + try: + return build_fdtd_config( + manifest, + material_snapshots, + background_material=self.background_material, + center_wavelength_nm=self.wavelength_um * 1000, + wavelength_halfspan_nm=self.wavelength_halfspan_um * 1000, + num_wavelengths=self.num_wavelengths, + default_port=default_port, + nanometers_per_cell=self.nanometers_per_cell, + pml_cells=self.pml_cells, + max_timesteps=self.max_timesteps, + energy_decay_fraction=self.energy_decay_fraction, + max_wall_seconds=self.max_wall_seconds, + gaussian_beam=gaussian_beam, + monitors=monitors, + ) + except ValidationError as error: + raise FDTDConfigError( + f"Invalid GDSFactory FDTD configuration: {error}" + ) from error + + def _selected_port_name(self) -> str: + """Choose an explicit port or prefer the first guided port by default.""" + if self.default_port is not None: + if self.default_port not in self.resolved.ports: + raise FDTDConfigError( + f"default_port {self.default_port!r} is not present on the " + "resolved component." + ) + return self.default_port + for name, port in self.resolved.ports.items(): + if not port.is_vertical: + return name + return next(iter(self.resolved.ports)) + + @staticmethod + def _vertical_polarization(port: Any) -> tuple[float, float, float]: + """Map vertical TE/TM and its in-plane orientation to an E direction.""" + angle = radians(port.orientation) + if port.port_type == "vertical_te": + vector = (-sin(angle), cos(angle), 0.0) + elif port.port_type == "vertical_tm": + vector = (cos(angle), sin(angle), 0.0) + else: + raise FDTDConfigError(f"Unsupported vertical port type {port.port_type!r}.") + return ( + round(vector[0], 12), + round(vector[1], 12), + round(vector[2], 12), + ) + + def _vertical_port_configs( + self, + material_snapshots: Mapping[str, MaterialSnapshot], + ) -> dict[str, tuple[GaussianBeamConfig, PlaneMonitorConfig]]: + """Translate vertical markers into free-space source/monitor settings.""" + vertical_ports = { + name: port for name, port in self.resolved.ports.items() if port.is_vertical + } + if not vertical_ports: + return {} + domain_bounds = background_bounds_nm( + self.resolved, + self.background_material, + self.background_padding_um, + ) + background_index = material_snapshots[self.background_material].refractive_index + outward_sign = 1 if self.vertical_port_axis == "+z" else -1 + inward_axis = "-z" if outward_sign > 0 else "+z" + propagation_inward = (0.0, 0.0, float(-outward_sign)) + propagation_outward = (0.0, 0.0, float(outward_sign)) + aperture_z_nm = domain_bounds[5] if outward_sign > 0 else domain_bounds[2] + device_z_nm = ( + self.resolved.bounds[1][2] + if outward_sign > 0 + else self.resolved.bounds[0][2] + ) * 1000 + configs = {} + for name, port in vertical_ports.items(): + aperture_width_um = self.vertical_port_aperture_width_um or port.width + waist_radius_um = ( + self.vertical_port_waist_radius_um or aperture_width_um / 2 + ) + half_width_nm = aperture_width_um * 500 + center_x_nm = port.center[0] * 1000 + center_y_nm = port.center[1] * 1000 + region_min = ( + center_x_nm - half_width_nm, + center_y_nm - half_width_nm, + aperture_z_nm, + ) + region_max = ( + center_x_nm + half_width_nm, + center_y_nm + half_width_nm, + aperture_z_nm, + ) + if ( + region_min[0] < domain_bounds[0] + or region_min[1] < domain_bounds[1] + or region_max[0] > domain_bounds[3] + or region_max[1] > domain_bounds[4] + ): + raise FDTDConfigError( + f"Vertical port {name!r} aperture exceeds the background " + "domain; increase background_padding_um or reduce " + "vertical_port_aperture_width_um." + ) + focal_point = (center_x_nm, center_y_nm, device_z_nm) + polarization = self._vertical_polarization(port) + common_fiber = FiberModeConfig( + propagation_direction=propagation_outward, + e_polarization=polarization, + focal_point=focal_point, + waist_radius=waist_radius_um * 1000, + refractive_index=background_index, + ) + configs[name] = ( + GaussianBeamConfig( + region_min=region_min, + region_max=region_max, + aperture_normal=inward_axis, + propagation_direction=propagation_inward, + e_polarization=polarization, + focal_point=focal_point, + waist_radius=waist_radius_um * 1000, + refractive_index=background_index, + ), + PlaneMonitorConfig( + name=name, + region_min=region_min, + region_max=region_max, + normal=self.vertical_port_axis, + fiber_mode=common_fiber, + ), + ) + return configs + + def write(self, output_dir: str | Path) -> SimulationArtifacts: + """Write ``mesh.msh`` and ``config.json`` into an output directory.""" + resolved = self.resolved + directory = Path(output_dir) + directory.mkdir(parents=True, exist_ok=True) + mesh_path = directory / "mesh.msh" + config_path = directory / "config.json" + material_snapshots = self._material_snapshots() + manifest = generate_mesh( + resolved, + mesh_path, + background_material=self.background_material, + background_padding_um=self.background_padding_um, + mesh_size_nm=self.mesh_size_nm, + nanometers_per_cell=self.nanometers_per_cell, + ) + config = self._config(manifest, material_snapshots) + config_path.write_text( + config.model_dump_json(indent=2, exclude_none=True) + "\n", + encoding="utf8", + ) + return SimulationArtifacts( + mesh_path=mesh_path, + config_path=config_path, + manifest=manifest, + ) + + +__all__ = ["Simulation"] diff --git a/tests/common/materials/test_snapshots.py b/tests/common/materials/test_snapshots.py index 30b8583a..7e790231 100644 --- a/tests/common/materials/test_snapshots.py +++ b/tests/common/materials/test_snapshots.py @@ -17,6 +17,7 @@ [ ("Si-Salzberg", 3.477723756), ("Si-Li-293K", 3.4757), + ("SiN-Luke", 1.996279731714), ("SiO2-Malitson", 1.444023622), ], ) @@ -37,6 +38,13 @@ def test_tabulated_material_interpolates() -> None: assert snapshot.refractive_index == pytest.approx(3.4778) +def test_sin_fallback_alias_uses_luke_card() -> None: + snapshot = resolve_material_snapshot("SiN", 1.55, {}) + + assert snapshot.refractive_index == pytest.approx(1.996279731714) + assert snapshot.source == "gsim" + + def test_project_card_overrides_fallback_and_missing_card_uses_fallback() -> None: project_si = GSIM_MATERIAL_CARDS["Si-Li-293K"].model_copy(update={"name": "Si"}) diff --git a/tests/common/pdk/test_resolve.py b/tests/common/pdk/test_resolve.py index a8d5bd57..04478fa2 100644 --- a/tests/common/pdk/test_resolve.py +++ b/tests/common/pdk/test_resolve.py @@ -138,6 +138,18 @@ def test_preserves_derived_geometry_and_authoritative_layer_fields() -> None: assert (10, 0) in result.derived_component.layers +def test_preexisting_derived_target_does_not_replace_source_geometry() -> None: + component = demo_component() + component.add_polygon( + [(-10, -10), (10, -10), (10, 10), (-10, 10)], + layer=(10, 0), + ) + + result = resolve_passive_pcell(component, pdk=make_test_pdk()) + + assert result.layers["core_key"].geometry.area == pytest.approx(0.875) + + def test_accepts_callable_and_instantiated_components() -> None: pdk = make_test_pdk() @@ -221,6 +233,26 @@ def test_rejects_non_axis_aligned_ports() -> None: ) +def test_preserves_vertical_port_type_and_in_plane_polarization_angle() -> None: + component = demo_component() + component.add_port( + name="fiber", + center=(1.0, 0.1), + width=10, + orientation=45, + layer=(1, 0), + port_type="vertical_te", + ) + + result = resolve_passive_pcell(component, pdk=make_test_pdk()) + + fiber = result.ports["fiber"] + assert fiber.is_vertical + assert fiber.port_type == "vertical_te" + assert fiber.orientation == pytest.approx(45) + assert fiber.normal == (0, 0, 1) + + def test_rejects_settings_for_component_instance() -> None: with pytest.raises(ComponentResolutionError, match="Settings cannot"): resolve_passive_pcell( diff --git a/tests/fdtd/__init__.py b/tests/fdtd/__init__.py new file mode 100644 index 00000000..03733887 --- /dev/null +++ b/tests/fdtd/__init__.py @@ -0,0 +1 @@ +"""Tests for passive GDSFactory FDTD artifact generation.""" diff --git a/tests/fdtd/conftest.py b/tests/fdtd/conftest.py new file mode 100644 index 00000000..739b0a17 --- /dev/null +++ b/tests/fdtd/conftest.py @@ -0,0 +1,95 @@ +"""Fixtures for PDK-native GDSFactory FDTD artifact tests.""" + +from __future__ import annotations + +from types import SimpleNamespace + +import gdsfactory as gf +import pytest +from gdsfactory.technology import LayerLevel, LayerStack, LogicalLayer + +from gsim.common.materials import GSIM_MATERIAL_CARDS + + +def straight_component(length: float = 2.0) -> gf.Component: + """Return a minimal two-port waveguide on the test core layer.""" + component = gf.Component() + component.add_polygon( + [(0, -0.25), (length, -0.25), (length, 0.25), (0, 0.25)], + layer=(1, 0), + ) + component.add_port( + name="o1", + center=(0, 0), + width=0.5, + orientation=180, + layer=(1, 0), + ) + component.add_port( + name="o2", + center=(length, 0), + width=0.5, + orientation=0, + layer=(1, 0), + ) + return component + + +def vertical_coupler_component(length: float = 2.0) -> gf.Component: + """Return a waveguide with one guided and one vertical optical port.""" + component = gf.Component() + component.add_polygon( + [(0, -0.25), (length, -0.25), (length, 0.25), (0, 0.25)], + layer=(1, 0), + ) + component.add_port( + name="o1", + center=(0, 0), + width=0.5, + orientation=180, + layer=(1, 0), + ) + component.add_port( + name="o2", + center=(length / 2, 0), + width=1, + orientation=0, + layer=(1, 0), + port_type="vertical_te", + ) + return component + + +@pytest.fixture +def fdtd_pdk_module() -> SimpleNamespace: + """Return a PDK module with project Si and fallback-only SiO2.""" + layer_stack = LayerStack( + layers={ + "core": LayerLevel( + layer=LogicalLayer(layer=(1, 0)), + thickness=0.22, + zmin=0, + sidewall_angle=10, + width_to_z=0.5, + mesh_order=2, + material="Si", + ), + "buried_oxide": LayerLevel( + layer=LogicalLayer(layer=(99, 0)), + thickness=2, + zmin=-1, + mesh_order=4, + material="SiO2", + ), + } + ) + pdk = gf.Pdk( + name="fdtd_test_pdk", + cells={ + "straight": straight_component, + "vertical_coupler": vertical_coupler_component, + }, + layer_stack=layer_stack, + ) + project_si = GSIM_MATERIAL_CARDS["Si-Li-293K"].model_copy(update={"name": "Si"}) + return SimpleNamespace(PDK=pdk, MATERIAL_CARDS={"Si": project_si}) diff --git a/tests/fdtd/test_config.py b/tests/fdtd/test_config.py new file mode 100644 index 00000000..e9673dd1 --- /dev/null +++ b/tests/fdtd/test_config.py @@ -0,0 +1,165 @@ +"""Tests for the strict GDSFactory FDTD configuration boundary.""" + +from __future__ import annotations + +from dataclasses import replace + +import pytest +from pydantic import ValidationError + +from gsim.common.materials import resolve_material_snapshot +from gsim.fdtd.config import ( + FDTDConfig, + FiberModeConfig, + GaussianBeamConfig, + PlaneMonitorConfig, + build_fdtd_config, +) +from gsim.fdtd.models import ( + FDTDConfigError, + MeshGroup, + MeshManifest, + PortMeshGroup, +) + + +def _manifest() -> MeshManifest: + return MeshManifest( + volumes={ + "background": MeshGroup( + name="background", + physical_tag=11, + material="SiO2", + priority=0, + ) + }, + layers={ + "core": MeshGroup( + name="core", + physical_tag=17, + material="Si", + priority=1, + ) + }, + ports={ + "o1": PortMeshGroup( + name="o1", + physical_name="port_o1", + physical_tag=23, + layer="core", + normal=(-1, 0, 0), + ) + }, + ) + + +def _config() -> FDTDConfig: + snapshots = { + name: resolve_material_snapshot(name, 1.55, {}) for name in ("Si", "SiO2") + } + return build_fdtd_config( + _manifest(), + snapshots, + background_material="SiO2", + center_wavelength_nm=1550, + wavelength_halfspan_nm=50, + num_wavelengths=3, + default_port="o1", + nanometers_per_cell=31.25, + pml_cells=16, + max_timesteps=None, + energy_decay_fraction=1e-6, + max_wall_seconds=3600, + ) + + +def test_config_uses_manifest_tags_and_rejects_extra_fields() -> None: + config = _config() + + assert config.length_scale_meters == 1e-9 + assert config.geometry.volumes["background"].phys_group == 11 + assert config.geometry.layers["core"].phys_group == 17 + assert config.geometry.ports["o1"].phys_group == 23 + + document = config.model_dump() + document["unsupported"] = True + with pytest.raises(ValidationError, match="Extra inputs are not permitted"): + FDTDConfig.model_validate(document) + + document = config.model_dump() + document["excitation"]["default_port"] = "missing" + with pytest.raises(ValidationError, match="is not declared"): + FDTDConfig.model_validate(document) + + +def test_config_rejects_lossy_scalar_materials() -> None: + silicon = resolve_material_snapshot("Si", 1.55, {}) + silica = resolve_material_snapshot("SiO2", 1.55, {}) + lossy_silicon = replace(silicon, extinction_coefficient=0.01) + + with pytest.raises(FDTDConfigError, match="lossless real"): + build_fdtd_config( + _manifest(), + {"Si": lossy_silicon, "SiO2": silica}, + background_material="SiO2", + center_wavelength_nm=1550, + wavelength_halfspan_nm=50, + num_wavelengths=3, + default_port="o1", + nanometers_per_cell=31.25, + pml_cells=16, + max_timesteps=None, + energy_decay_fraction=1e-6, + max_wall_seconds=3600, + ) + + +def test_config_supports_gaussian_beam_and_fiber_monitor() -> None: + beam = GaussianBeamConfig( + region_min=(0, 0, 3000), + region_max=(10000, 10000, 3000), + aperture_normal="-z", + propagation_direction=(0, 0, -1), + e_polarization=(0, 1, 0), + focal_point=(5000, 5000, 220), + waist_radius=5000, + refractive_index=1.444, + ) + monitor = PlaneMonitorConfig( + name="fiber", + region_min=(0, 0, 3000), + region_max=(10000, 10000, 3000), + normal="+z", + fiber_mode=FiberModeConfig( + propagation_direction=(0, 0, 1), + e_polarization=(0, 1, 0), + focal_point=(5000, 5000, 220), + waist_radius=5000, + refractive_index=1.444, + ), + ) + snapshots = { + name: resolve_material_snapshot(name, 1.55, {}) for name in ("Si", "SiO2") + } + + config = build_fdtd_config( + _manifest(), + snapshots, + background_material="SiO2", + center_wavelength_nm=1550, + wavelength_halfspan_nm=50, + num_wavelengths=3, + default_port=None, + nanometers_per_cell=31.25, + pml_cells=16, + max_timesteps=None, + energy_decay_fraction=1e-6, + max_wall_seconds=3600, + gaussian_beam=beam, + monitors=[monitor], + ) + + assert config.excitation.type == "gaussian_beam" + assert config.excitation.gaussian_beam == beam + assert config.excitation.default_port is None + assert config.monitors == [monitor] diff --git a/tests/fdtd/test_mesh.py b/tests/fdtd/test_mesh.py new file mode 100644 index 00000000..7b15c162 --- /dev/null +++ b/tests/fdtd/test_mesh.py @@ -0,0 +1,73 @@ +"""Tests for solver-specific mesh semantics.""" + +from __future__ import annotations + +from types import SimpleNamespace + +import pytest +from shapely.geometry import MultiPolygon, Polygon + +from gsim.common.pdk import ResolvedLayer +from gsim.fdtd.mesh import _priority_by_mesh_order +from gsim.fdtd.mesh_geometry import iter_polygons, sidewall_slice_count + + +def test_lower_pdk_mesh_order_becomes_higher_fdtd_priority() -> None: + layers = { + "core": SimpleNamespace(mesh_order=1), + "same_order": SimpleNamespace(mesh_order=1), + "slab": SimpleNamespace(mesh_order=4), + "cladding": SimpleNamespace(mesh_order=7), + } + + priorities = _priority_by_mesh_order(layers) + + assert priorities == { + "core": 3, + "same_order": 3, + "slab": 2, + "cladding": 1, + } + + +def _resolved_layer(geometry: Polygon | MultiPolygon) -> ResolvedLayer: + return ResolvedLayer( + key="core", + declared_name="core", + layer=(1, 0), + derived_layer=None, + geometry=geometry, + material="Si", + zmin=0, + thickness=0.22, + zmax=0.22, + sidewall_angle=10, + width_to_z=0.5, + bias=None, + z_to_bias=None, + mesh_order=1, + ) + + +def test_sidewall_slices_bound_error_to_quarter_cell() -> None: + layer = _resolved_layer(Polygon([(0, 0), (1, 0), (1, 1), (0, 1)])) + + slice_count = sidewall_slice_count(layer, nanometers_per_cell=31.25) + + assert slice_count == 3 + total_displacement_nm = 38.7919357558623 + assert total_displacement_nm / (2 * slice_count) < 31.25 / 4 + + +def test_disconnected_polygons_retain_holes() -> None: + ring = Polygon( + [(0, 0), (4, 0), (4, 4), (0, 4)], + holes=[[(1, 1), (3, 1), (3, 3), (1, 3)]], + ) + bus = Polygon([(5, 0), (7, 0), (7, 1), (5, 1)]) + + polygons = iter_polygons(MultiPolygon([ring, bus]), layer_key="core") + + assert len(polygons) == 2 + assert sum(len(polygon.interiors) for polygon in polygons) == 1 + assert sum(polygon.area for polygon in polygons) == pytest.approx(14) diff --git a/tests/fdtd/test_simulation.py b/tests/fdtd/test_simulation.py new file mode 100644 index 00000000..654bfd26 --- /dev/null +++ b/tests/fdtd/test_simulation.py @@ -0,0 +1,170 @@ +"""Behavioral tests for the public GDSFactory FDTD artifact workflow.""" + +from __future__ import annotations + +import json + +import meshio +import numpy as np +import pytest + +from gsim import fdtd +from gsim.fdtd.models import FDTDGeometryError + + +def _physical_group_points(mesh: meshio.Mesh, name: str, cell_type: str) -> np.ndarray: + """Return mesh points used by one named physical group.""" + physical_tag = mesh.field_data[name][0] + cells = mesh.cells_dict[cell_type] + tags = mesh.cell_data_dict["gmsh:physical"][cell_type] + point_indices = np.unique(cells[tags == physical_tag].ravel()) + return mesh.points[point_indices] + + +def test_write_uses_project_material_then_fallback_and_valid_mesh( + tmp_path, + fdtd_pdk_module, +) -> None: + simulation = fdtd.Simulation( + pdk=fdtd_pdk_module, + mesh_size_nm=750, + background_padding_um=0.25, + pml_cells=16, + num_wavelengths=3, + ) + resolved = simulation.geometry("straight", settings={"length": 2.0}) + + assert resolved.materials["Si"].source == "project" + assert resolved.materials["Si"].refractive_index == pytest.approx(3.4757) + + artifacts = simulation.write(tmp_path) + document = json.loads(artifacts.config_path.read_text(encoding="utf8")) + mesh = meshio.read(artifacts.mesh_path) + + assert document["schema_version"] == 1 + assert document["mesh_file"] == "mesh.msh" + assert document["length_scale_meters"] == 1e-9 + assert document["materials"]["Si"]["refractive_index"] == pytest.approx(3.4757) + assert document["materials"]["SiO2"]["refractive_index"] == pytest.approx( + 1.4440236217 + ) + assert document["geometry"]["ports"]["o1"]["normal"] == [-1, 0, 0] + assert document["geometry"]["ports"]["o2"]["normal"] == [1, 0, 0] + + expected_names = {"background", "core", "port_o1", "port_o2"} + assert set(mesh.field_data) == expected_names + assert {block.type for block in mesh.cells} == {"triangle", "tetra"} + for group in artifacts.manifest.volumes.values(): + assert mesh.field_data[group.name].tolist() == [group.physical_tag, 3] + for group in artifacts.manifest.layers.values(): + assert mesh.field_data[group.name].tolist() == [group.physical_tag, 3] + for port in artifacts.manifest.ports.values(): + assert mesh.field_data[port.physical_name].tolist() == [ + port.physical_tag, + 2, + ] + + port_points = _physical_group_points(mesh, "port_o1", "triangle") + assert port_points[:, 0] == pytest.approx(0) + assert port_points[:, 1].min() == pytest.approx(-262.930645, abs=1e-3) + assert port_points[:, 1].max() == pytest.approx(262.930645, abs=1e-3) + assert port_points[:, 2].min() == pytest.approx(0) + assert port_points[:, 2].max() == pytest.approx(220) + + mesh_header = artifacts.mesh_path.read_text(encoding="utf8").splitlines() + assert mesh_header[mesh_header.index("$MeshFormat") + 1] == "2.2 0 8" + + +def test_write_requires_geometry(tmp_path) -> None: + with pytest.raises(FDTDGeometryError, match=r"Call Simulation\.geometry"): + fdtd.Simulation().write(tmp_path) + + +def test_vertical_port_becomes_fiber_monitor_not_material_port( + tmp_path, + fdtd_pdk_module, +) -> None: + simulation = fdtd.Simulation( + pdk=fdtd_pdk_module, + mesh_size_nm=750, + background_padding_um=0.25, + ) + simulation.geometry("vertical_coupler") + + artifacts = simulation.write(tmp_path) + document = json.loads(artifacts.config_path.read_text(encoding="utf8")) + mesh = meshio.read(artifacts.mesh_path) + + assert document["excitation"]["type"] == "eigenmode" + assert document["excitation"]["default_port"] == "o1" + assert set(document["geometry"]["ports"]) == {"o1"} + assert "port_o2" not in mesh.field_data + assert document["monitors"][0]["name"] == "o2" + assert document["monitors"][0]["normal"] == "+z" + assert document["monitors"][0]["fiber_mode"]["e_polarization"] == [ + -0.0, + 1.0, + 0.0, + ] + + +def test_vertical_default_port_uses_gaussian_beam( + tmp_path, + fdtd_pdk_module, +) -> None: + simulation = fdtd.Simulation( + pdk=fdtd_pdk_module, + default_port="o2", + mesh_size_nm=750, + background_padding_um=0.25, + ) + simulation.geometry("vertical_coupler") + + artifacts = simulation.write(tmp_path) + document = json.loads(artifacts.config_path.read_text(encoding="utf8")) + + assert document["excitation"]["type"] == "gaussian_beam" + assert "default_port" not in document["excitation"] + assert document["excitation"]["gaussian_beam"]["aperture_normal"] == "-z" + assert document["excitation"]["gaussian_beam"]["propagation_direction"] == [ + 0.0, + 0.0, + -1.0, + ] + + +def test_write_restores_options_from_existing_gmsh_session( + tmp_path, + fdtd_pdk_module, +) -> None: + import gmsh + + original_options = { + "Mesh.Binary": 1.0, + "Mesh.ElementOrder": 2.0, + "Mesh.MeshSizeExtendFromBoundary": 1.0, + "Mesh.MeshSizeFromCurvature": 1.0, + "Mesh.MeshSizeMax": 34.5, + "Mesh.MeshSizeMin": 12.5, + "Mesh.MshFileVersion": 4.1, + "Mesh.SaveAll": 1.0, + } + gmsh.initialize() + try: + for option_name, value in original_options.items(): + gmsh.option.setNumber(option_name, value) + + simulation = fdtd.Simulation( + pdk=fdtd_pdk_module, + mesh_size_nm=750, + background_padding_um=0.25, + ) + simulation.geometry("straight", settings={"length": 2.0}) + simulation.write(tmp_path) + + assert gmsh.isInitialized() + for option_name, expected in original_options.items(): + assert gmsh.option.getNumber(option_name) == pytest.approx(expected) + finally: + gmsh.clear() + gmsh.finalize() diff --git a/tests/palace/test_patterned_dielectrics.py b/tests/palace/test_patterned_dielectrics.py index af41d239..7d748ab0 100644 --- a/tests/palace/test_patterned_dielectrics.py +++ b/tests/palace/test_patterned_dielectrics.py @@ -4,12 +4,26 @@ from types import SimpleNamespace +import pytest from gdsfactory.technology import LayerLevel from gsim.common.stack.extractor import Layer, LayerStack, extract_layer_stack from gsim.palace.mesh.geometry import build_entities +@pytest.fixture +def gmsh_session(): + """Provide an isolated Gmsh session and always release it.""" + import gmsh + + gmsh.initialize() + try: + yield gmsh + finally: + gmsh.clear() + gmsh.finalize() + + def _fake_gf_stack(): """Build a minimal gdsfactory-like stack for extractor tests.""" return SimpleNamespace( @@ -99,7 +113,7 @@ def test_build_entities_prioritizes_patterned_dielectrics_over_background_boxes( assert orders["air"] == 4 -def test_add_patterned_dielectrics_skips_covered_dielectric_layers(): +def test_add_patterned_dielectrics_skips_covered_dielectric_layers(gmsh_session): """Layers covered by bulk dielectric boxes must not be re-extruded.""" import gdsfactory as gf import klayout.db as kdb @@ -163,9 +177,7 @@ def test_add_patterned_dielectrics_skips_covered_dielectric_layers(): geometry = extract_geometry(c, stack) - import gmsh - - gmsh.initialize() + gmsh = gmsh_session gmsh.option.setNumber("General.Verbosity", 0) gmsh.model.add("test") kernel = gmsh.model.occ