FXC-2087 Implement automatic structure extrusion for boundary waveports

CHANGELOG.md

@@ -10,6 +10,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ### Added
 - Added S-parameter de-embedding to `TerminalComponentModelerData`, enabling recalculation with shifted reference planes.
+- 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 support for `tidy3d-extras`, an optional plugin that enables more accurate local mode solving via subpixel averaging.
 - Added support for `symlog` and `log` scale plotting in `Scene.plot_eps()` and `Scene.plot_structures_property()` methods. The `symlog` scale provides linear behavior near zero and logarithmic behavior elsewhere, while 'log' is a base 10 logarithmic scale.
 - Added `LowFrequencySmoothingSpec` and `ModelerLowFrequencySmoothingSpec` for automatic smoothing of mode monitor data at low frequencies where DFT sampling is insufficient.

schemas/TerminalComponentModeler.json

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

tests/test_plugins/smatrix/terminal_component_modeler_def.py

@@ -5,6 +5,7 @@
 import numpy as np
 
 import tidy3d as td
+import tidy3d.plugins.microwave as mw
 from tidy3d.plugins.smatrix import (
     CoaxialLumpedPort,
     LumpedPort,
@@ -331,3 +332,145 @@ 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, L)), 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 = mw.AxisAlignedCurrentIntegral(
+        center=((se + w) / 2, 0, -waveport_z / 2), size=(2 * w, 3 * t, 0), sign="+"
+    )
+    voltage_integral = mw.AxisAlignedVoltageIntegral(
+        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

@@ -29,7 +29,11 @@
 from tidy3d.plugins.smatrix.utils import s_to_z, 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
 
@@ -1595,3 +1599,182 @@ def test_S_parameter_deembedding(monkeypatch, tmp_path):
     S_dmb_shortcut = modeler_data_LP.smatrix_deembedded(port_shifts=port_shifts_LP)
     assert not np.allclose(S_dmb.data.values, s_matrix_LP.data.values)
     assert np.allclose(S_dmb_shortcut.data.values, S_dmb.data.values)
+
+
+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[-4].geometry.geometries[0].slab_bounds),
+        np.max(sim.structures[-2].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))
+
+    # generate a new TCM simulation to test edge case when wave port plane does not intersect any structures
+    tcm = make_coaxial_component_modeler(
+        length=100000, port_types=(WavePort, WavePort), use_current=False
+    )
+
+    # 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)
+    port_2 = port_2.updated_copy(center=(0, 0, 50000), extrude_structures=True)
+
+    # move wave port plane so that is does not intersect any structures
+    port_1_center_new = (638.4, 0.0, -51000)
+
+    # update voltage integral
+    voltage_int = port_1.voltage_integral.updated_copy(center=port_1_center_new)
+    # update WavePort
+    port_1 = port_1.updated_copy(center=port_1_center_new, voltage_integral=voltage_int)
+
+    # update component modeler
+    tcm = tcm.updated_copy(ports=[port_1, port_2])
+
+    # make sure that
+    with pytest.raises(SetupError):
+        sim = tcm.base_sim
+
+
+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[-6].geometry.geometries[0].slab_bounds),
+        np.max(sim.structures[-1].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))
+
+    # test scenario when wave port extrusion is requested, but port plane does not intersect any structures
+    mil = 25.4
+    port_1_center_new = (0, 0, -2010 * mil)
+
+    # re-assemble a new component modeler
+    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)
+    port_2 = port_2.updated_copy(extrude_structures=True)
+
+    # update current and voltage integrals
+    current_int = port_1.current_integral.updated_copy(center=port_1_center_new)
+    voltage_int = port_1.voltage_integral.updated_copy(center=port_1_center_new)
+
+    # update WavePort
+    port_1 = port_1.updated_copy(
+        center=port_1_center_new, current_integral=current_int, voltage_integral=voltage_int
+    )
+
+    # update component modeler
+    tcm = tcm.updated_copy(ports=[port_1, port_2])
+
+    # make sure that the error is triggered
+    with pytest.raises(SetupError):
+        sim = tcm.base_sim
+
+    # update component modeler
+    tcm = tcm.updated_copy(ports=[port_2, port_1])
+
+    # make sure that the error is triggered even when ports are reshuffled
+    with pytest.raises(SetupError):
+        sim = tcm.base_sim

tidy3d/components/geometry/utils.py

@@ -735,10 +735,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(