diff --git a/README.md b/README.md index 74837b64..7d8fe5d6 100644 --- a/README.md +++ b/README.md @@ -14,6 +14,7 @@ with minimal boilerplate. | Module | Solver | Method | Use Case | | ------------- | ------------------------------------------- | ------ | ------------------------------------------------------ | +| `gsim.fdtd` | ZapFDTD | 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..571cb343 --- /dev/null +++ b/docs/api/fdtd.md @@ -0,0 +1,54 @@ +# FDTD API + +`gsim.fdtd` generates the coarse tetrahedral mesh and validated `config.json` +consumed by ZapFDTD. 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 ZapFDTD. + +## Initial geometry limits + +The first backend supports axis-aligned optical ports, one connected polygon per +material-bearing layer, and vertical or constant-angle sidewalls. It rejects +ambiguous geometry, unsupported `bias`/`z_to_bias` profiles, and lossy material +snapshots because ZapFDTD 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.ZapConfig + options: + show_source: false diff --git a/docs/index.md b/docs/index.md index 31d97a65..9c6d366b 100644 --- a/docs/index.md +++ b/docs/index.md @@ -12,6 +12,7 @@ with minimal boilerplate. | Module | Solver | Method | Use Case | | ------------- | ------------------------------------------- | ------ | ------------------------------------------------------ | +| `gsim.fdtd` | ZapFDTD | 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: [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..94f19136 100644 --- a/docs/zensical.toml +++ b/docs/zensical.toml @@ -29,6 +29,7 @@ nav = [ ] }, { "Changelog" = "CHANGELOG.md" }, { "API Reference" = [ + { "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..da65e6a8 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 ZapFDTD 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/fdtd/__init__.py b/src/gsim/fdtd/__init__.py new file mode 100644 index 00000000..037cdc3f --- /dev/null +++ b/src/gsim/fdtd/__init__.py @@ -0,0 +1,21 @@ +"""Passive FDTD artifact generation for ZapFDTD.""" + +from gsim.fdtd.config import ZapConfig +from gsim.fdtd.models import ( + FDTDArtifactError, + FDTDConfigError, + FDTDGeometryError, + MeshManifest, + SimulationArtifacts, +) +from gsim.fdtd.simulation import Simulation + +__all__ = [ + "FDTDArtifactError", + "FDTDConfigError", + "FDTDGeometryError", + "MeshManifest", + "Simulation", + "SimulationArtifacts", + "ZapConfig", +] diff --git a/src/gsim/fdtd/config.py b/src/gsim/fdtd/config.py new file mode 100644 index 00000000..55284de8 --- /dev/null +++ b/src/gsim/fdtd/config.py @@ -0,0 +1,225 @@ +"""Validated ZapFDTD 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 ZapFDTD does not understand.""" + + model_config = ConfigDict(extra="forbid", allow_inf_nan=False) + + +class MaterialConfig(_StrictModel): + """Scalar real optical material supported by Zap 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 ZapFDTD.""" + + volumes: dict[str, RegionConfig] = Field(min_length=1) + layers: dict[str, RegionConfig] = Field(min_length=1) + ports: dict[str, PortConfig] = Field(min_length=1) + + +class ExcitationConfig(_StrictModel): + """Initial eigenmode pulse configuration.""" + + type: Literal["eigenmode"] = "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) + default_port: str = Field(min_length=1) + + @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") + 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): + """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 ZapConfig(_StrictModel): + """Complete ZapFDTD 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 + grid: GridConfig + run: RunConfig + + @model_validator(mode="after") + def validate_references(self) -> ZapConfig: + """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.default_port not in self.geometry.ports: + raise ValueError( + f"default_port {self.excitation.default_port!r} is not declared" + ) + return self + + +def _material_config(snapshot: MaterialSnapshot) -> MaterialConfig: + """Convert one lossless scalar snapshot to the Zap material schema.""" + if snapshot.extinction_coefficient != 0: + raise FDTDConfigError( + f"Material {snapshot.material_name!r} has extinction coefficient " + f"{snapshot.extinction_coefficient}; Zap schema v1 supports only " + "lossless real refractive indices." + ) + return MaterialConfig(refractive_index=snapshot.refractive_index) + + +def build_zap_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, + nanometers_per_cell: float, + pml_cells: int, + max_timesteps: int | None, + energy_decay_fraction: float, + max_wall_seconds: float, +) -> ZapConfig: + """Build and cross-validate a Zap 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 ZapConfig( + 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( + center_wavelength=center_wavelength_nm, + wavelength_halfspan=wavelength_halfspan_nm, + num_wavelengths=num_wavelengths, + default_port=default_port, + ), + 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", + "GeometryConfig", + "GridConfig", + "MaterialConfig", + "PortConfig", + "RegionConfig", + "RunConfig", + "ZapConfig", + "build_zap_config", +] diff --git a/src/gsim/fdtd/mesh.py b/src/gsim/fdtd/mesh.py new file mode 100644 index 00000000..3a783dd3 --- /dev/null +++ b/src/gsim/fdtd/mesh.py @@ -0,0 +1,483 @@ +"""Coarse Gmsh artifact generation for ZapFDTD voxelization.""" + +from __future__ import annotations + +from collections.abc import Iterable, Mapping +from math import radians, tan +from pathlib import Path +from typing import Any + +import gmsh +from shapely.geometry import MultiPolygon, Polygon, box +from shapely.geometry.polygon import orient +from shapely.validation import explain_validity + +from gsim.common.pdk import ResolvedLayer, ResolvedPassivePcell, ResolvedPort +from gsim.fdtd.mesh_validation import validate_mesh +from gsim.fdtd.models import ( + FDTDGeometryError, + MeshGroup, + MeshManifest, + PortMeshGroup, +) + +_UM_TO_NM = 1000.0 +_GEOMETRY_TOLERANCE_NM = 1e-3 + + +def _single_polygon(layer: ResolvedLayer) -> Polygon: + """Return one valid connected polygon for an initial FDTD layer.""" + geometry = layer.geometry + if isinstance(geometry, Polygon): + polygon = geometry + elif isinstance(geometry, MultiPolygon) and len(geometry.geoms) == 1: + polygon = geometry.geoms[0] + else: + count = len(geometry.geoms) if hasattr(geometry, "geoms") else 0 + raise FDTDGeometryError( + f"Layer {layer.key!r} has {count} disconnected polygons; the initial " + "FDTD backend requires one connected solid per layer." + ) + 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 polygon + + +def _scaled_ring(coordinates: Iterable[tuple[float, ...]]) -> list[tuple[float, float]]: + """Convert one 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_wire(kernel: Any, polygon: Polygon, z_nm: float) -> int: + """Create one closed OCC wire for lofting.""" + points = _scaled_ring(orient(polygon, sign=1.0).exterior.coords) + if len(points) < 3: + raise FDTDGeometryError("A loft profile has fewer than three vertices.") + first_index = min( + range(len(points)), + key=lambda index: (points[index][0], points[index][1]), + ) + points = points[first_index:] + points[:first_index] + point_tags = [kernel.addPoint(x, y, z_nm) for x, y in points] + line_tags = [ + kernel.addLine(point_tags[index], point_tags[(index + 1) % len(point_tags)]) + for index in range(len(point_tags)) + ] + return kernel.addWire(line_tags, checkClosed=True) + + +def _profile_offset_um(layer: ResolvedLayer, normalized_z: float) -> float: + """Return the PDK sidewall offset at one normalized z position.""" + return ( + (layer.width_to_z - normalized_z) + * abs(layer.thickness) + * tan(radians(layer.sidewall_angle)) + ) + + +def _condition_profile_at_ports( + profile: Polygon, + ports: list[ResolvedPort], + offset_um: float, +) -> Polygon: + """Extend and clip a profile so every port is a full planar end face.""" + 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"Port {port.name!r} is not in the component plane." + ) + transverse_axis = 1 - normal_axis + half_width = port.width / 2 + offset_um + if half_width <= 0: + raise FDTDGeometryError( + f"Layer {port.layer_key!r} sidewall offset closes port {port.name!r}." + ) + + target = port.center[normal_axis] + transverse_lower = port.center[transverse_axis] - half_width + transverse_upper = port.center[transverse_axis] + half_width + xmin, ymin, xmax, ymax = profile.bounds + epsilon = _GEOMETRY_TOLERANCE_NM / _UM_TO_NM + if normal_axis == 0: + edge = xmin if port.normal[0] < 0 else xmax + inner = max(target + epsilon, edge + epsilon) + if port.normal[0] > 0: + inner = min(target - epsilon, edge - epsilon) + extension = box( + min(target, inner), + transverse_lower, + max(target, inner), + transverse_upper, + ) + else: + edge = ymin if port.normal[1] < 0 else ymax + inner = max(target + epsilon, edge + epsilon) + if port.normal[1] > 0: + inner = min(target - epsilon, edge - epsilon) + extension = box( + transverse_lower, + min(target, inner), + transverse_upper, + max(target, inner), + ) + profile = profile.union(extension) + + xmin, ymin, xmax, ymax = profile.bounds + margin = max(xmax - xmin, ymax - ymin, port.width, 1.0) + if normal_axis == 0 and port.normal[0] < 0: + clip = box(target, ymin - margin, xmax + margin, ymax + margin) + elif normal_axis == 0: + clip = box(xmin - margin, ymin - margin, target, ymax + margin) + elif port.normal[1] < 0: + clip = box(xmin - margin, target, xmax + margin, ymax + margin) + else: + clip = box(xmin - margin, ymin - margin, xmax + margin, target) + profile = profile.intersection(clip) + + profile = profile.simplify( + 10 * _GEOMETRY_TOLERANCE_NM / _UM_TO_NM, + preserve_topology=True, + ) + if not isinstance(profile, Polygon) or profile.is_empty or not profile.is_valid: + raise FDTDGeometryError("Port conditioning produced invalid layer geometry.") + return profile + + +def _offset_profile( + polygon: Polygon, + layer: ResolvedLayer, + ports: list[ResolvedPort], + normalized_z: float, +) -> Polygon: + """Apply the PDK sidewall offset and preserve full port end faces.""" + offset_um = _profile_offset_um(layer, normalized_z) + profile = polygon.buffer(offset_um, join_style=2) + if not isinstance(profile, Polygon) or profile.is_empty or not profile.is_valid: + raise FDTDGeometryError( + f"Layer {layer.key!r} sidewall profile becomes empty or disconnected " + f"at normalized z={normalized_z:g}." + ) + return _condition_profile_at_ports(profile, ports, offset_um) + + +def _add_layer_volume( + kernel: Any, + layer: ResolvedLayer, + ports: list[ResolvedPort], +) -> list[int]: + """Create one vertical or sidewall-tapered OCC layer volume.""" + 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 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." + ) + + polygon = _single_polygon(layer) + z_lower_um, z_upper_um = layer.z_bounds + z_lower_nm = z_lower_um * _UM_TO_NM + z_upper_nm = z_upper_um * _UM_TO_NM + if layer.sidewall_angle == 0: + from gsim.palace.mesh.gmsh_utils import extrude_polygon + + polygon = _condition_profile_at_ports(polygon, ports, 0.0) + exterior = _scaled_ring(polygon.exterior.coords) + hole_coordinates = [] + for interior in polygon.interiors: + points = _scaled_ring(interior.coords) + hole_coordinates.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_nm, + z_upper_nm - z_lower_nm, + holes=hole_coordinates, + ) + if volume_tag is None: + raise FDTDGeometryError(f"Could not extrude layer {layer.key!r}.") + return [volume_tag] + + if polygon.interiors: + raise FDTDGeometryError( + f"Layer {layer.key!r} combines holes and sidewalls, which is not " + "supported by the initial FDTD mesh writer." + ) + lower_profile = _offset_profile(polygon, layer, ports, 0.0) + upper_profile = _offset_profile(polygon, layer, ports, 1.0) + lower_wire = _add_wire(kernel, lower_profile, z_lower_nm) + upper_wire = _add_wire(kernel, upper_profile, z_upper_nm) + dimtags = kernel.addThruSections( + [lower_wire, upper_wire], + makeSolid=True, + makeRuled=True, + ) + volume_tags = [tag for dimension, tag in dimtags if dimension == 3] + if not volume_tags: + raise FDTDGeometryError(f"Could not loft layer {layer.key!r}.") + return volume_tags + + +def _priority_by_mesh_order( + layers: Mapping[str, ResolvedLayer], +) -> dict[str, int]: + """Invert lower-wins PDK mesh order into higher-wins Zap 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() + } + 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) + target_nm = port.center[normal_axis] * _UM_TO_NM + center_nm = tuple(coordinate * _UM_TO_NM for coordinate in port.center) + 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 + transverse_axes = [axis for axis in range(3) if axis != normal_axis] + if all( + bounds[axis] - _GEOMETRY_TOLERANCE_NM + <= center_nm[axis] + <= bounds[axis + 3] + _GEOMETRY_TOLERANCE_NM + for axis in transverse_axes + ): + 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}." + ) + 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 ZapFDTD port extrusion." + ) + + +def generate_mesh( + resolved: ResolvedPassivePcell, + mesh_path: Path, + *, + background_material: str, + background_padding_um: float, + mesh_size_nm: float, +) -> MeshManifest: + """Generate and validate a coarse Zap-compatible 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 "background" in resolved.layers: + raise FDTDGeometryError("Layer name 'background' is reserved by FDTD.") + + initialized_here = not bool(gmsh.isInitialized()) + 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, + ) + for port in resolved.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_volume( + kernel, + layer, + [port for port in resolved.ports.values() if port.layer_key == name], + ) + 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 resolved.ports.items(): + surface_tags = _port_surface_tags( + port, + layer_volume_tags[port.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=port.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() + + validate_mesh(mesh_path, manifest) + return manifest + + +__all__ = ["generate_mesh", "validate_mesh"] diff --git a/src/gsim/fdtd/mesh_validation.py b/src/gsim/fdtd/mesh_validation.py new file mode 100644 index 00000000..6fc66a05 --- /dev/null +++ b/src/gsim/fdtd/mesh_validation.py @@ -0,0 +1,71 @@ +"""Strict validation for ZapFDTD-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 ZapFDTD 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("ZapFDTD 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("ZapFDTD 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..feaa8982 --- /dev/null +++ b/src/gsim/fdtd/models.py @@ -0,0 +1,68 @@ +"""Shared models and errors for FDTD artifact generation.""" + +from __future__ import annotations + +from dataclasses import dataclass +from pathlib import Path + + +class FDTDArtifactError(ValueError): + """Base error for invalid or unsupported FDTD artifacts.""" + + +class FDTDGeometryError(FDTDArtifactError): + """Raised when resolved geometry cannot produce a valid Zap mesh.""" + + +class FDTDConfigError(FDTDArtifactError): + """Raised when resolved settings cannot produce a valid Zap 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..cbe5d5cd --- /dev/null +++ b/src/gsim/fdtd/simulation.py @@ -0,0 +1,176 @@ +"""Public Simulation workflow for ZapFDTD artifact generation.""" + +from __future__ import annotations + +from collections.abc import Mapping +from dataclasses import dataclass, field +from math import isfinite +from pathlib import Path +from typing import Any + +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 ZapConfig, build_zap_config +from gsim.fdtd.mesh import generate_mesh +from gsim.fdtd.models import ( + FDTDConfigError, + FDTDGeometryError, + MeshManifest, + SimulationArtifacts, +) + + +@dataclass +class Simulation: + """Generate coarse Gmsh and config artifacts for passive ZapFDTD 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 + 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.") + + @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], + ) -> ZapConfig: + """Build the validated Zap schema after mesh group tags are known.""" + if not self.resolved.ports: + raise FDTDConfigError("Eigenmode FDTD requires at least one port.") + default_port = self.default_port or next(iter(self.resolved.ports)) + try: + return build_zap_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, + ) + except ValidationError as error: + raise FDTDConfigError(f"Invalid ZapFDTD configuration: {error}") from error + + 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, + ) + config = self._config(manifest, material_snapshots) + config_path.write_text(config.model_dump_json(indent=2) + "\n", encoding="utf8") + return SimulationArtifacts( + mesh_path=mesh_path, + config_path=config_path, + manifest=manifest, + ) + + +__all__ = ["Simulation"] diff --git a/tests/fdtd/__init__.py b/tests/fdtd/__init__.py new file mode 100644 index 00000000..25a2a04b --- /dev/null +++ b/tests/fdtd/__init__.py @@ -0,0 +1 @@ +"""Tests for passive ZapFDTD artifact generation.""" diff --git a/tests/fdtd/conftest.py b/tests/fdtd/conftest.py new file mode 100644 index 00000000..2897acd0 --- /dev/null +++ b/tests/fdtd/conftest.py @@ -0,0 +1,67 @@ +"""Fixtures for PDK-native 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 + + +@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}, + 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..f1776de0 --- /dev/null +++ b/tests/fdtd/test_config.py @@ -0,0 +1,108 @@ +"""Tests for the strict ZapFDTD 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 ZapConfig, build_zap_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() -> ZapConfig: + snapshots = { + name: resolve_material_snapshot(name, 1.55, {}) for name in ("Si", "SiO2") + } + return build_zap_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"): + ZapConfig.model_validate(document) + + document = config.model_dump() + document["excitation"]["default_port"] = "missing" + with pytest.raises(ValidationError, match="is not declared"): + ZapConfig.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_zap_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, + ) diff --git a/tests/fdtd/test_mesh.py b/tests/fdtd/test_mesh.py new file mode 100644 index 00000000..328e26cf --- /dev/null +++ b/tests/fdtd/test_mesh.py @@ -0,0 +1,25 @@ +"""Tests for solver-specific mesh semantics.""" + +from __future__ import annotations + +from types import SimpleNamespace + +from gsim.fdtd.mesh import _priority_by_mesh_order + + +def test_lower_pdk_mesh_order_becomes_higher_zap_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, + } diff --git a/tests/fdtd/test_simulation.py b/tests/fdtd/test_simulation.py new file mode 100644 index 00000000..c85ed554 --- /dev/null +++ b/tests/fdtd/test_simulation.py @@ -0,0 +1,80 @@ +"""Behavioral tests for the public 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(-269.39597, abs=1e-3) + assert port_points[:, 1].max() == pytest.approx(269.39597, 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)