pec frame

Author: dbochkov-flexcompute

tidy3d/__init__.py

@@ -421,6 +421,7 @@ def set_logging_level(level: str) -> None:
     "PMC",
     "PML",
     "TFSF",
+    "ABCBoundary",
     "Absorber",
     "AbsorberParams",
     "AbstractFieldProjectionData",
@@ -593,6 +594,7 @@ def set_logging_level(level: str) -> None:
     "Medium2D",
     "MediumMediumInterface",
     "MeshOverrideStructure",
+    "ModeABCBoundary",
     "ModeAmpsDataArray",
     "ModeData",
     "ModeIndexDataArray",
@@ -716,6 +718,4 @@ def set_logging_level(level: str) -> None:
     "set_logging_console",
     "set_logging_file",
     "wavelengths",
-    "ABCBoundary",
-    "ModeABCBoundary",
 ]

tidy3d/components/boundary.py

@@ -3,7 +3,7 @@
 from __future__ import annotations
 
 from abc import ABC
-from typing import Union
+from typing import Optional, Union
 
 import numpy as np
 import pydantic.v1 as pd
@@ -46,6 +46,8 @@ def _warn_num_layers(cls, val):
 
     return _warn_num_layers
 
+DEFAULT_MODE_SPEC_MODE_ABC = ModeSpec()
+
 
 class BoundaryEdge(ABC, Tidy3dBaseModel):
     """Electromagnetic boundary condition at a domain edge."""
@@ -119,7 +121,7 @@ class ModeABCBoundary(AbstractABCBoundary):
     """One-way wave equation absorbing boundary conditions for absorbing a waveguide mode."""
 
     mode_spec: ModeSpec = pd.Field(
-        ModeSpec(),
+        DEFAULT_MODE_SPEC_MODE_ABC,
         title="Mode Specification",
         description="Parameters that determine the modes computed by the mode solver.",
     )
@@ -894,7 +896,7 @@ def abc(
     def mode_abc(
         cls,
         plane: Box,
-        mode_spec: ModeSpec = ModeSpec(),
+        mode_spec: ModeSpec = DEFAULT_MODE_SPEC_MODE_ABC,
         mode_index: pd.NonNegativeInt = 0,
         frequency: Optional[pd.PositiveFloat] = None,
     ):

tidy3d/components/simulation.py

@@ -49,7 +49,7 @@
 from .data.unstructured.tetrahedral import TetrahedralGridDataset
 from .data.unstructured.triangular import TriangularGridDataset
 from .data.utils import CustomSpatialDataType
-from .geometry.base import Box, Geometry
+from .geometry.base import Box, ClipOperation, Geometry, GeometryGroup
 from .geometry.mesh import TriangleMesh
 from .geometry.utils import flatten_groups, traverse_geometries
 from .geometry.utils_2d import get_bounds, get_thickened_geom, snap_coordinate_to_grid, subdivide
@@ -67,6 +67,7 @@
     Medium2D,
     MediumType,
     MediumType3D,
+    PECMedium,
 )
 from .monitor import (
     AbstractFieldProjectionMonitor,
@@ -175,6 +176,10 @@
 # RF frequency warning
 RF_FREQ_WARNING = 300e9
 
+# length and thickness of optional PEC frames around mode sources (in cells)
+MODE_PEC_FRAME_LENGTH = 2
+MODE_PEC_FRAME_THICKNESS = 1e-3
+
 
 def validate_boundaries_for_zero_dims(warn_on_change: bool = True):
     """Error if absorbing boundaries, bloch boundaries, unmatching pec/pmc, or symmetry is used along a zero dimension."""
@@ -849,7 +854,7 @@ def pml_thicknesses(self) -> list[tuple[float, float]]:
 
         Returns
         -------
-        List[Tuple[float, float]]
+        list[Tuple[float, float]]
             List containing the absorber thickness (micron) in - and + boundaries.
         """
         num_layers = self.num_pml_layers
@@ -867,7 +872,7 @@ def internal_override_structures(self) -> list[MeshOverrideStructure]:
 
         Returns
         -------
-        List[MeshOverrideSructure]
+        list[MeshOverrideSructure]
             List of override structures.
         """
         wavelength = self.grid_spec.get_wavelength(self.sources)
@@ -884,7 +889,7 @@ def internal_snapping_points(self) -> list[CoordinateOptional]:
 
         Returns
         -------
-        List[CoordinateOptional]
+        list[CoordinateOptional]
             List of snapping points coordinates.
         """
         return self.grid_spec.internal_snapping_points(
@@ -1223,7 +1228,7 @@ def _grid_and_snapping_lines(self) -> tuple[Grid, list[CoordinateOptional]]:
 
         Returns
         -------
-        Tuple[:class:`.Grid`, List[CoordinateOptional]]
+        Tuple[:class:`.Grid`, list[CoordinateOptional]]
             :class:`.Grid` storing the spatial locations relevant to the simulation
             the list of snapping points generated during iterative gap meshing.
         """
@@ -1281,7 +1286,7 @@ def _gap_meshing_snapping_lines(self) -> list[CoordinateOptional]:
 
         Returns
         -------
-        List[CoordinateOptional]
+        list[CoordinateOptional]
             List of snapping lines resolving thin gaps and strips.
         """
 
@@ -1345,7 +1350,7 @@ def num_pml_layers(self) -> list[tuple[float, float]]:
 
         Returns
         -------
-        List[Tuple[float, float]]
+        list[Tuple[float, float]]
             List containing the number of absorber layers in - and + boundaries.
         """
         num_layers = [[0, 0], [0, 0], [0, 0]]
@@ -3285,13 +3290,13 @@ def _projection_monitors_homogeneous(cls, val, values):
     @classmethod
     def _get_mediums_on_abc(
         cls, boundary_spec, medium, center, size, structures
-    ) -> Tuple[
-        List[MediumType3D],
-        List[MediumType3D],
-        List[MediumType3D],
-        List[MediumType3D],
-        List[MediumType3D],
-        List[MediumType3D],
+    ) -> tuple[
+        list[MediumType3D],
+        list[MediumType3D],
+        list[MediumType3D],
+        list[MediumType3D],
+        list[MediumType3D],
+        list[MediumType3D],
     ]:
         """For each ABC boundary that needs an automatic medium detection (permittivity=None)
         determine mediums it crosses.
@@ -3310,7 +3315,7 @@ def _get_mediums_on_abc(
             center=structure_bg.geometry.center, size=structure_bg.geometry.size
         )
 
-        total_structures = [structure_bg] + list(structures)
+        total_structures = [structure_bg, *list(structures)]
 
         mediums = []
         for boundary, surface in zip(np.ravel(boundary_spec.to_list), surfaces):
@@ -4544,7 +4549,7 @@ def mediums(self) -> set[MediumType]:
 
         Returns
         -------
-        List[:class:`.AbstractMedium`]
+        set[:class:`.AbstractMedium`]
             Set of distinct mediums in the simulation.
         """
         log.warning(
@@ -4606,12 +4611,12 @@ def intersecting_media(
         -------
         test_object : :class:`.Box`
             Object for which intersecting media are to be detected.
-        structures : List[:class:`.AbstractMedium`]
+        structures : tuple[:class:`.AbstractMedium`]
             List of structures whose media will be tested.
 
         Returns
         -------
-        List[:class:`.AbstractMedium`]
+        tuple[:class:`.AbstractMedium`]
             Set of distinct mediums that intersect with the given planar object.
         """
 
@@ -4633,12 +4638,12 @@ def intersecting_structures(
         -------
         test_object : :class:`.Box`
             Object for which intersecting media are to be detected.
-        structures : List[:class:`.AbstractMedium`]
+        structures : tuple[:class:`.AbstractMedium`]
             List of structures whose media will be tested.
 
         Returns
         -------
-        List[:class:`.Structure`]
+        tuple[:class:`.Structure`]
             Set of distinct structures that intersect with the given surface, or with the surfaces
             of the given volume.
         """
@@ -5451,3 +5456,101 @@ def from_scene(cls, scene: Scene, **kwargs) -> Simulation:
         )
 
     _boundaries_for_zero_dims = validate_boundaries_for_zero_dims()
+
+    def _make_pec_frame(self, mode_source) -> Structure:
+        """Make a pec frame around a mode source."""
+
+        coords = self.grid.boundaries.to_list
+        axis = mode_source.injection_axis
+        direction = mode_source.direction
+        length = mode_source.pec_frame
+
+        span_inds = np.array(self.grid.discretize_inds(mode_source))
+        if direction == "+":
+            span_inds[axis][1] += length - 1
+        else:
+            span_inds[axis][0] -= length - 1
+
+        box_bounds = [
+            [
+                c[beg],
+                c[end],
+            ]
+            for c, (beg, end) in zip(coords, span_inds)
+        ]
+
+        prev_cell = span_inds[axis][0] - 1
+        if prev_cell >= 0:
+            box_bounds[axis][0] = (1 - MODE_PEC_FRAME_THICKNESS) * box_bounds[axis][0] + MODE_PEC_FRAME_THICKNESS * coords[axis][prev_cell]
+
+        next_cell = span_inds[axis][1] + 1
+        if next_cell <= len(coords[axis]) - 1 :
+            box_bounds[axis][1] = (1 - MODE_PEC_FRAME_THICKNESS) * box_bounds[axis][1] + MODE_PEC_FRAME_THICKNESS * coords[axis][next_cell]
+
+        box = Box.from_bounds(*np.transpose(box_bounds))
+        
+        surfaces = Box.surfaces(box.size, box.center)
+        del surfaces[2 * axis: 2 * axis + 2]
+
+        structure = Structure(
+            geometry=GeometryGroup(
+                geometries=surfaces,
+            ),
+            medium=PECMedium(),
+        )
+
+
+        # bounds_outer = [
+        #     [
+        #         (1 - MODE_PEC_FRAME_THICKNESS) * c[beg]
+        #         + MODE_PEC_FRAME_THICKNESS * c[max(0, beg - 1)],
+        #         (1 - MODE_PEC_FRAME_THICKNESS) * c[end]
+        #         + MODE_PEC_FRAME_THICKNESS * c[min(len(c) - 1, end + 1)],
+        #     ]
+        #     for c, (beg, end) in zip(coords, span_inds)
+        # ]
+        # bounds_inner = [
+        #     [
+        #         (1 - MODE_PEC_FRAME_THICKNESS) * c[beg] + MODE_PEC_FRAME_THICKNESS * c[beg + 1],
+        #         (1 - MODE_PEC_FRAME_THICKNESS) * c[end] + MODE_PEC_FRAME_THICKNESS * c[end - 1],
+        #     ]
+        #     for c, (beg, end) in zip(coords, span_inds)
+        # ]
+        # bounds_inner[axis] = [-inf, inf]
+        # structure = Structure(
+        #     geometry=ClipOperation(
+        #         geometry_a=Box.from_bounds(*np.transpose(bounds_outer)),
+        #         geometry_b=Box.from_bounds(*np.transpose(bounds_inner)),
+        #         operation="difference",
+        #     ),
+        #     medium=PECMedium(),
+        # )
+        return structure
+
+    @cached_property
+    def with_mode_source_pec_frames(self) -> Simulation:
+        """Return an instance with added pec frames around mode sources."""
+
+        pec_frames = [
+            self._make_pec_frame(src)
+            for src in self.sources
+            if isinstance(src, ModeSource) and src.pec_frame > 0
+        ]
+
+        if len(pec_frames) == 0:
+            return self
+
+        return self.updated_copy(
+            grid_spec=GridSpec.from_grid(self.grid), structures=list(self.structures) + pec_frames
+        )
+
+    def _validate_with_mode_source_pec_frames(self):
+        """Validate that after adding pec frames simulation setup is still valid."""
+
+        try:
+            _ = self.with_mode_source_pec_frames
+        except Exception:
+            log.error(
+                "Simulation fails after requested mode source PEC frames are added. "
+                "Please inspec '.with_mode_source_pec_frames'."
+            )

tidy3d/components/source/field.py

@@ -403,6 +403,12 @@ class ModeSource(DirectionalSource, PlanarSource, BroadbandSource):
         "``num_modes`` in the solver will be set to ``mode_index + 1``.",
     )
 
+    pec_frame: pydantic.NonNegativeInt = pydantic.Field(
+        0,
+        title="PEC Frame.",
+        description="Add a thin pec frame around the source during FDTD run.",
+    )
+
     @cached_property
     def angle_theta(self):
         """Polar angle of propagation."""

tidy3d/plugins/smatrix/ports/wave.py

@@ -82,6 +82,12 @@ class WavePort(AbstractTerminalPort, Box):
         description="Use conjugated or non-conjugated dot product for mode decomposition.",
     )
 
+    pec_frame: pd.NonNegativeInt = pd.Field(
+        0,
+        title="PEC Frame.",
+        description="Add a thin pec frame around the source during FDTD run.",
+    )
+
     def _mode_voltage_coefficients(self, mode_data: ModeData) -> FreqModeDataArray:
         """Calculates scaling coefficients to convert mode amplitudes
         to the total port voltage.
@@ -139,6 +145,7 @@ def to_source(
             mode_index=self.mode_index,
             direction=self.direction,
             name=self.name,
+            pec_frame=self.pec_frame,
         )
 
     def to_monitors(