Implement automatic structure extrusion for boundary waveports

CHANGELOG.md

@@ -8,6 +8,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 ## [Unreleased]
 
 ### Added
+- Added optional automatic extrusion of structures intersecting with a `WavePort` via the new `extrude_structures` field, ensuring mode sources, absorbers, and PEC frames are fully contained.
 - Added rectangular and radial taper support to `RectangularAntennaArrayCalculator` for phased array amplitude weighting; refactored array factor calculation for improved clarity and performance. 
 - Selective simulation capabilities to `TerminalComponentModeler` via `run_only` and `element_mappings` fields, allowing users to run fewer simulations and extract only needed scattering matrix elements.
 - Added KLayout plugin, with DRC functionality for running design rule checks in `plugins.klayout.drc`. Supports running DRC on GDS files as well as `Geometry`, `Structure`, and `Simulation` objects.

schemas/TerminalComponentModeler.json

@@ -16129,6 +16129,10 @@
           ],
           "type": "string"
         },
+        "extrude_structures": {
+          "default": false,
+          "type": "boolean"
+        },
         "frame": {
           "allOf": [
             {

tests/test_plugins/smatrix/terminal_component_modeler_def.py

@@ -332,3 +332,147 @@ def make_port(center, direction, type, name) -> Union[CoaxialLumpedPort, WavePor
     )
 
     return modeler
+
+
+def make_differential_stripline_modeler():
+    # Frequency range (Hz)
+    f_min, f_max = (1e9, 70e9)
+
+    # Frequency sample points
+    freqs = np.linspace(f_min, f_max, 101)
+
+    # Geometry
+    mil = 25.4  # conversion to mils to microns (default unit)
+    w = 3.2 * mil  # Signal strip width
+    t = 0.7 * mil  # Conductor thickness
+    h = 10.7 * mil  # Substrate thickness
+    se = 7 * mil  # gap between edge-coupled pair
+    L = 4000 * mil  # Line length
+    len_inf = 1e6  # Effective infinity
+
+    left_end = -L / 2
+    right_end = len_inf
+
+    len_z = right_end - left_end
+    cent_z = (left_end + right_end) / 2
+    waveport_z = L
+
+    # Material properties
+    eps = 4.4  # Relative permittivity, substrate
+
+    # define media
+    med_sub = td.Medium(permittivity=eps)
+    med_metal = td.PEC
+
+    left_strip_geometry = td.Box(center=(-(se + w) / 2, 0, 0), size=(w, t, L))
+    right_strip_geometry = td.Box(center=((se + w) / 2, 0, 0), size=(w, t, L))
+
+    # Substrate
+    str_sub = td.Structure(
+        geometry=td.Box(center=(0, 0, 0), size=(len_inf, h, len_inf)), medium=med_sub
+    )
+
+    # disjoint signal strips
+    str_signal_strips = td.Structure(
+        geometry=td.GeometryGroup(geometries=[left_strip_geometry, right_strip_geometry]),
+        medium=med_metal,
+    )
+
+    # Top ground plane
+    str_gnd_top = td.Structure(
+        geometry=td.Box(center=(0, h / 2 + t / 2, 0), size=(len_inf, t, L)), medium=med_metal
+    )
+
+    # Bottom ground plane
+    str_gnd_bot = td.Structure(
+        geometry=td.Box(center=(0, -h / 2 - t / 2, 0), size=(len_inf, t, L)), medium=med_metal
+    )
+
+    # Create a LayerRefinementSpec from signal trace structures
+    lr_spec = td.LayerRefinementSpec.from_structures(
+        structures=[str_signal_strips],
+        axis=1,  # Layer normal is in y-direction
+        min_steps_along_axis=10,  # Min 10 grid cells along normal direction
+        refinement_inside_sim_only=False,  # Metal structures extend outside sim domain. Set 'False' to snap to corners outside sim.
+        bounds_snapping="bounds",  # snap grid to metal boundaries
+        corner_refinement=td.GridRefinement(
+            dl=t / 10, num_cells=2
+        ),  # snap to corners and apply added refinement
+    )
+
+    # Layer refinement for top and bottom ground planes
+    lr_spec2 = lr_spec.updated_copy(center=(0, h / 2 + t / 2, cent_z), size=(len_inf, t, len_z))
+    lr_spec3 = lr_spec.updated_copy(center=(0, -h / 2 - t / 2, cent_z), size=(len_inf, t, len_z))
+
+    # Define overall grid specification
+    grid_spec = td.GridSpec.auto(
+        wavelength=td.C_0 / f_max,
+        min_steps_per_wvl=30,
+        layer_refinement_specs=[lr_spec, lr_spec2, lr_spec3],
+    )
+
+    # boundary specs
+    boundary_spec = td.BoundarySpec(
+        x=td.Boundary.pml(),
+        y=td.Boundary.pec(),
+        z=td.Boundary.pml(),
+    )
+
+    # Define port specification
+    wave_port_mode_spec = td.ModeSpec(num_modes=1, target_neff=np.sqrt(eps))
+
+    # Define current and voltage integrals
+    current_integral = microwave.CurrentIntegralAxisAligned(
+        center=((se + w) / 2, 0, -waveport_z / 2), size=(2 * w, 3 * t, 0), sign="+"
+    )
+    voltage_integral = microwave.VoltageIntegralAxisAligned(
+        center=(0, 0, -waveport_z / 2),
+        size=(se, 0, 0),
+        extrapolate_to_endpoints=True,
+        snap_path_to_grid=True,
+        sign="+",
+    )
+
+    # Define wave ports
+    WP1 = WavePort(
+        center=(0, 0, -waveport_z / 2),
+        size=(len_inf, len_inf, 0),
+        mode_spec=wave_port_mode_spec,
+        direction="+",
+        name="WP1",
+        mode_index=0,
+        current_integral=current_integral,
+        voltage_integral=voltage_integral,
+    )
+    WP2 = WP1.updated_copy(
+        name="WP2",
+        center=(0, 0, waveport_z / 2),
+        direction="-",
+        current_integral=current_integral.updated_copy(
+            center=((se + w) / 2, 0, waveport_z / 2), sign="-"
+        ),
+        voltage_integral=voltage_integral.updated_copy(center=(0, 0, waveport_z / 2)),
+    )
+
+    # define fimulation
+    sim = td.Simulation(
+        size=(50 * mil, h + 2 * t, 1.05 * L),
+        center=(0, 0, 0),
+        grid_spec=grid_spec,
+        boundary_spec=boundary_spec,
+        structures=[str_sub, str_signal_strips, str_gnd_top, str_gnd_bot],
+        monitors=[],
+        run_time=2e-9,  # simulation run time in seconds
+        shutoff=1e-7,  # lower shutoff threshold for more accurate low frequency
+        plot_length_units="mm",
+        symmetry=(-1, 0, 0),  # odd symmetry in x-direction
+    )
+
+    # set up component modeler
+    tcm = TerminalComponentModeler(
+        simulation=sim,  # simulation, previously defined
+        ports=[WP1, WP2],  # wave ports, previously defined
+        freqs=freqs,  # S-parameter frequency points
+    )
+
+    return tcm

tests/test_plugins/smatrix/test_terminal_component_modeler.py

@@ -33,7 +33,11 @@
 from tidy3d.plugins.smatrix.utils import validate_square_matrix
 
 from ...utils import run_emulated
-from .terminal_component_modeler_def import make_coaxial_component_modeler, make_component_modeler
+from .terminal_component_modeler_def import (
+    make_coaxial_component_modeler,
+    make_component_modeler,
+    make_differential_stripline_modeler,
+)
 
 mm = 1e3
 
@@ -1351,3 +1355,118 @@ def test_wave_port_to_absorber(tmp_path):
     sim = list(modeler.sim_dict.values())[0]
     absorber = sim.internal_absorbers[0]
     assert absorber.boundary_spec == custom_boundary_spec
+
+
+def test_wave_port_extrusion_coaxial():
+    """Test extrusion of structures wave port absorber."""
+
+    # define a terminal component modeler
+    tcm = make_coaxial_component_modeler(
+        length=100000,
+        port_types=(WavePort, WavePort),
+    )
+
+    # update ports and set flag to extrude structures
+    ports = tcm.ports
+    port_1 = ports[0]
+    port_2 = ports[1]
+    port_1 = port_1.updated_copy(center=(0, 0, -50000), extrude_structures=True)
+
+    # test that structure extrusion requires an internal absorber (should raise ValidationError)
+    with pytest.raises(pd.ValidationError):
+        _ = port_2.updated_copy(center=(0, 0, 50000), extrude_structures=True, absorber=False)
+
+    # define a valid waveport
+    port_2 = port_2.updated_copy(center=(0, 0, 50000), extrude_structures=True)
+
+    # update component modeler
+    tcm = tcm.updated_copy(ports=[port_1, port_2])
+
+    # generate simulations from component modeler
+    sim = tcm.base_sim
+
+    # get injection axis that would be used to extrude structure
+    inj_axis = sim.internal_absorbers[0].size.index(0.0)
+
+    # get grid boundaries
+    bnd_coords = sim.grid.boundaries.to_list[inj_axis]
+
+    # get size of structures along injection axis directions
+    str_bnds = [
+        np.min(sim.structures[0].geometry.geometries[1].geometries[0].slab_bounds),
+        np.max(sim.structures[0].geometry.geometries[0].slab_bounds),
+    ]
+
+    pec_bnds = []
+
+    # infer placement of PEC plates beyond internal absorber
+    for absorber in sim.internal_absorbers:
+        absorber_cntr = absorber.center[inj_axis]
+        right_ind = np.searchsorted(bnd_coords, absorber_cntr, side="right")
+        left_ind = np.searchsorted(bnd_coords, absorber_cntr, side="left") - 1
+        pec_bnds.append(bnd_coords[right_ind + 1])
+        pec_bnds.append(bnd_coords[left_ind - 1])
+
+    # get range of coordinates along injection axis for PEC plates
+    pec_bnds = [np.min(pec_bnds), np.max(pec_bnds)]
+
+    # ensure that structures were extruded up to PEC plates
+    assert all(np.isclose(str_bnd, pec_bnd) for str_bnd, pec_bnd in zip(str_bnds, pec_bnds))
+
+
+def test_wave_port_extrusion_differential_stripline():
+    """Test extrusion of structures wave port absorber for differential stripline."""
+
+    tcm = make_differential_stripline_modeler()
+
+    # update ports and set flag to extrude structures
+    ports = tcm.ports
+    port_1 = ports[0]
+    port_2 = ports[1]
+    port_1 = port_1.updated_copy(extrude_structures=True)
+
+    # test that structure extrusion requires an internal absorber (should raise ValidationError)
+    with pytest.raises(pd.ValidationError):
+        _ = port_2.updated_copy(extrude_structures=True, absorber=False)
+
+    # define a valid waveport
+    port_2 = port_2.updated_copy(extrude_structures=True)
+
+    # update component modeler
+    tcm = tcm.updated_copy(ports=[port_1, port_2])
+
+    # generate simulations from component modeler
+    sim = tcm.base_sim
+
+    # get injection axis that would be used to extrude structure
+    inj_axis = sim.internal_absorbers[0].size.index(0.0)
+
+    # get grid boundaries
+    bnd_coords = sim.grid.boundaries.to_list[inj_axis]
+
+    # get size of structures along injection axis directions
+    str_bnds = [
+        np.min(sim.structures[0].geometry.geometries[1].geometries[0].slab_bounds),
+        np.max(sim.structures[0].geometry.geometries[0].slab_bounds),
+    ]
+
+    pec_bnds = []
+
+    # infer placement of PEC plates beyond internal absorber
+    for absorber in sim._shifted_internal_absorbers:
+        # get the PEC box with its face surfaces
+        (box, inj_axis, direction) = sim._pec_frame_box(absorber)
+        surfaces = box.surfaces(box.size, box.center)
+
+        # get extrusion coordinates and a cutting plane for inference of intersecting structures.
+        sign = 1 if direction == "+" else -1
+        cutting_plane = surfaces[2 * inj_axis + (1 if direction == "+" else 0)]
+
+        # get extrusion extent along injection axis
+        pec_bnds.append(cutting_plane.center[inj_axis])
+
+    # get range of coordinates along injection axis for PEC plates
+    pec_bnds = [np.min(pec_bnds), np.max(pec_bnds)]
+
+    # ensure that structures were extruded up to PEC plates
+    assert all(np.isclose(str_bnd, pec_bnd) for str_bnd, pec_bnd in zip(str_bnds, pec_bnds))

tidy3d/components/geometry/utils.py

@@ -522,10 +522,11 @@ def _shift_value_signed(
             f"{name} position '{obj_position}' is outside of simulation bounds '({grid_boundaries[0]}, {grid_boundaries[-1]})' along dimension '{'xyz'[normal_axis]}'."
         )
     obj_index = obj_pos_gt_grid_bounds[-1]
-
     # shift the obj to the left
     signed_shift = shift if direction == "+" else -shift
     if signed_shift < 0:
+        if np.isclose(obj_position, grid_boundaries[obj_index + 1]):
+            obj_index += 1
         shifted_index = obj_index + signed_shift
         if shifted_index < 0 or grid_centers[shifted_index] <= bounds[0][normal_axis]:
             raise SetupError(

tidy3d/components/simulation.py

@@ -5638,33 +5638,9 @@ def _make_pec_frame(self, obj: Union[ModeSource, InternalAbsorber]) -> Structure
         the frame is added around the injection plane. For internal absorbers, a backing pec
         plate is also added on the non-absorbing side.
         """
-        span_inds = np.array(self.grid.discretize_inds(obj))
-
-        coords = self.grid.boundaries.to_list
-        direction = obj.direction
-        if isinstance(obj, ModeSource):
-            axis = obj.injection_axis
-            length = obj.frame.length
-        else:
-            axis = obj.size.index(0.0)
-            length = 1
-
-        if direction == "+":
-            span_inds[axis][1] += length - 1
-            span_inds[axis][0] -= 1
-        else:
-            span_inds[axis][1] += 1
-            span_inds[axis][0] -= length - 1
-
-        box_bounds = [
-            [
-                c[beg],
-                c[end],
-            ]
-            for c, (beg, end) in zip(coords, span_inds)
-        ]
 
-        box = Box.from_bounds(*np.transpose(box_bounds))
+        # get pec frame bounding box, object's axis and direction
+        (box, axis, direction) = self._pec_frame_box(obj)
 
         surfaces = Box.surfaces(box.size, box.center)
         if isinstance(obj, ModeSource):
@@ -5723,3 +5699,30 @@ def _validate_finalized(self):
                 "Simulation fails after requested mode source PEC frames are added. "
                 "Please inspect '._finalized'."
             )
+
+    def _pec_frame_box(self, obj: Union[ModeSource, InternalAbsorber]) -> tuple[Box, int, str]:
+        """Return pec bounding box, frame axis and object's direction"""
+
+        span_inds = np.array(self.grid.discretize_inds(obj))
+
+        coords = self.grid.boundaries.to_list
+        direction = obj.direction
+        if isinstance(obj, ModeSource):
+            axis = obj.injection_axis
+            length = obj.frame.length
+            if direction == "+":
+                span_inds[axis][1] += length - 1
+            else:
+                span_inds[axis][0] -= length - 1
+        else:
+            axis = obj.size.index(0.0)
+
+        box_bounds = [
+            [
+                c[beg],
+                c[end],
+            ]
+            for c, (beg, end) in zip(coords, span_inds)
+        ]
+
+        return (Box.from_bounds(*np.transpose(box_bounds)), axis, direction)