diff --git a/CHANGELOG.md b/CHANGELOG.md index 56d6dc97ac..78b11a32ac 100644 --- a/CHANGELOG.md +++ b/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. diff --git a/schemas/TerminalComponentModeler.json b/schemas/TerminalComponentModeler.json index c8bb8bd320..70b3db86e3 100644 --- a/schemas/TerminalComponentModeler.json +++ b/schemas/TerminalComponentModeler.json @@ -16129,6 +16129,10 @@ ], "type": "string" }, + "extrude_structures": { + "default": false, + "type": "boolean" + }, "frame": { "allOf": [ { diff --git a/tests/test_plugins/smatrix/terminal_component_modeler_def.py b/tests/test_plugins/smatrix/terminal_component_modeler_def.py index 609699c637..c1945527ce 100644 --- a/tests/test_plugins/smatrix/terminal_component_modeler_def.py +++ b/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 diff --git a/tests/test_plugins/smatrix/test_terminal_component_modeler.py b/tests/test_plugins/smatrix/test_terminal_component_modeler.py index 420b06a4a9..9a918ab847 100644 --- a/tests/test_plugins/smatrix/test_terminal_component_modeler.py +++ b/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)) diff --git a/tidy3d/components/geometry/utils.py b/tidy3d/components/geometry/utils.py index ad875936cf..2b80eb1536 100644 --- a/tidy3d/components/geometry/utils.py +++ b/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( diff --git a/tidy3d/components/simulation.py b/tidy3d/components/simulation.py index 60d8411733..5207870598 100644 --- a/tidy3d/components/simulation.py +++ b/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) diff --git a/tidy3d/plugins/smatrix/component_modelers/terminal.py b/tidy3d/plugins/smatrix/component_modelers/terminal.py index 257f0dae6c..f0941bd8f3 100644 --- a/tidy3d/plugins/smatrix/component_modelers/terminal.py +++ b/tidy3d/plugins/smatrix/component_modelers/terminal.py @@ -7,8 +7,10 @@ import numpy as np import pydantic.v1 as pd +from tidy3d import ClipOperation, GeometryGroup, GridSpec, PolySlab from tidy3d.components.base import cached_property from tidy3d.components.boundary import BroadbandModeABCSpec +from tidy3d.components.geometry.utils import _shift_object from tidy3d.components.geometry.utils_2d import snap_coordinate_to_grid from tidy3d.components.index import SimulationMap from tidy3d.components.monitor import DirectivityMonitor @@ -16,7 +18,7 @@ from tidy3d.components.source.time import GaussianPulse from tidy3d.components.types import Ax, Complex from tidy3d.components.viz import add_ax_if_none, equal_aspect -from tidy3d.constants import C_0, OHM +from tidy3d.constants import C_0, OHM, fp_eps from tidy3d.exceptions import SetupError, Tidy3dKeyError, ValidationError from tidy3d.log import log from tidy3d.plugins.smatrix.component_modelers.base import ( @@ -260,6 +262,7 @@ def sim_dict(self) -> SimulationMap: # Now, create simulations with wave port sources and mode solver monitors for computing port modes for network_index in self.matrix_indices_run_sim: task_name, sim_with_src = self._add_source_to_sim(network_index) + # update simulation sim_dict[task_name] = sim_with_src # Check final simulations for grid size at ports @@ -344,8 +347,14 @@ def base_sim(self) -> Simulation: "internal_absorbers": new_absorbers, } - # This is the new default simulation will all shared components added - return sim_wo_source.copy(update=update_dict) + # update base simulation with updated set of shared components + sim_wo_source = sim_wo_source.copy(update=update_dict) + + # extrude port structures + sim_wo_source = self._extrude_port_structures(sim=sim_wo_source) + + # This is the new default simulation with all shared components added + return sim_wo_source def _add_source_to_sim(self, source_index: NetworkIndex) -> tuple[str, Simulation]: """Adds the source corresponding to the ``source_index`` to the base simulation.""" @@ -363,6 +372,7 @@ def _add_source_to_sim(self, source_index: NetworkIndex) -> tuple[str, Simulatio self._source_time, snap_center=new_port_center, grid=self.base_sim.grid ) task_name = self.get_task_name(port=port, mode_index=mode_index) + return (task_name, self.base_sim.updated_copy(sources=[port_source])) @cached_property @@ -472,5 +482,117 @@ def get_radiation_monitor_by_name(self, monitor_name: str) -> DirectivityMonitor return monitor raise Tidy3dKeyError(f"No radiation monitor named '{monitor_name}'.") + def _extrude_port_structures(self, sim: Simulation) -> Simulation: + """ + Extrude structures intersecting a port plane when a wave port lies on a structure boundary. + + This method checks wave ports with ``extrude_structures==True`` and automatically extends the boundary structures + to PEC plates associated with internal absorbers in the direction opposite to the mode source. + This ensures that mode sources and internal absorbers are fully contained within the extrusion. + + Parameters + ---------- + sim : Simulation + Simulation object containing mode sources, internal absorbers, and monitors, + after mesh overrides and snapping points are applied. + + Returns + ------- + Simulation + Updated simulation with extruded structures added to ``simulation.structures``. + """ + + # get all mode sources from TerminalComponentModeler that correspond to ports with ``extrude_structures`` flag set to ``True``. + for port in self.ports: + if isinstance(port, WavePort) and port.extrude_structures: + # compute snap_center and shift the internal absorber associated with the current port + snap_center = port.center[port.injection_axis] + self._shift_value_signed(port) + absorber = port.to_absorber(snap_center=snap_center) + shifted_absorber = _shift_object( + obj=absorber, + grid=sim.grid, + bounds=sim.bounds, + direction=absorber.direction, + shift=absorber.grid_shift, + ) + + # get the PEC box with its face surfaces + (box, inj_axis, direction) = sim._pec_frame_box(shifted_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 + extrude_to = cutting_plane.center[inj_axis] + + # move cutting plane beyond the waveport plane along the `ModeSource` injection direction. + center = list(cutting_plane.center) + center[inj_axis] = port.center[inj_axis] - sign * fp_eps * box.size[inj_axis] + cutting_plane = cutting_plane.updated_copy(center=center) + + # define extrusion bounds + extrusion_bounds = [center[inj_axis], extrude_to][::sign] + + new_structures = [] + + # loop over structures and extrude those that intersect a waveport plane + for structure in sim.structures: + # get geometries that intersect the plane on which the waveport is defined + + shapely_geom = cutting_plane.intersections_with(structure.geometry) + + polygon_list = [] + for geom in shapely_geom: + polygon_list = polygon_list + ClipOperation.to_polygon_list(geom) + + new_geoms = [] + # loop over identified geometries and extrude them + for polygon in polygon_list: + # construct outer shell of an extruded geometry first + exterior_vertices = np.array(polygon.exterior.coords) + outer_shell = PolySlab( + axis=inj_axis, slab_bounds=extrusion_bounds, vertices=exterior_vertices + ) + + # construct innner shells that represent holes + hole_polyslabs = [ + PolySlab( + axis=inj_axis, + slab_bounds=extrusion_bounds, + vertices=np.array(hole.coords), + ) + for hole in polygon.interiors + ] + + # construct final geometry by removing inner holes from outer shell + if hole_polyslabs: + holes = GeometryGroup(geometries=hole_polyslabs) + extruded_slab_new = ClipOperation( + operation="difference", geometry_a=outer_shell, geometry_b=holes + ) + else: + extruded_slab_new = outer_shell + + # append extruded geometry + new_geoms.append(extruded_slab_new) + + # add the original geometry + new_geoms.append(structure.geometry) + + # update structure and add it to the list + new_struct = structure.updated_copy( + geometry=GeometryGroup(geometries=new_geoms) + ) + new_structures.append(new_struct) + + # update structures in simulation while keeping the same grid + sim = sim.updated_copy( + grid_spec=GridSpec.from_grid(sim.grid), structures=new_structures + ) + + return sim + TerminalComponentModeler.update_forward_refs() diff --git a/tidy3d/plugins/smatrix/ports/wave.py b/tidy3d/plugins/smatrix/ports/wave.py index 291cb5421c..5077f6ba01 100644 --- a/tidy3d/plugins/smatrix/ports/wave.py +++ b/tidy3d/plugins/smatrix/ports/wave.py @@ -87,17 +87,23 @@ class WavePort(AbstractTerminalPort, Box): frame: Optional[PECFrame] = pd.Field( DEFAULT_WAVE_PORT_FRAME, - title="Source Frame.", + title="Source Frame", description="Add a thin frame around the source during FDTD run for an improved injection.", ) absorber: Union[bool, ABCBoundary, ModeABCBoundary] = pd.Field( True, - title="Absorber.", + title="Absorber", description="Place a mode absorber in the port. If ``True``, an automatically generated mode absorber is placed in the port. " "If ``ABCBoundary`` or ``ModeABCBoundary``, a mode absorber is placed in the port with the specified boundary conditions.", ) + extrude_structures: bool = pd.Field( + False, + title="Extrude Structures", + description="Extrudes structures that intersect the wave port plane by a few grid cells when ``True``, improving mode injection accuracy.", + ) + def _mode_voltage_coefficients(self, mode_data: ModeData) -> FreqModeDataArray: """Calculates scaling coefficients to convert mode amplitudes to the total port voltage. @@ -320,3 +326,15 @@ def _validate_current_integral_sign(cls, val, values): f"'current_integral' sign must match the '{name}' direction '{direction}'." ) return val + + @pd.root_validator(pre=False) + def _check_absorber_if_extruding_structures(cls, values): + """Raise validation error when ``extrude_structures`` is set to ``True`` + while ``absorber`` is set to ``False``.""" + + if values.get("extrude_structures") and not values.get("absorber"): + raise ValidationError( + "Structure extrusion for a waveport requires an internal absorber. Set `absorber=True` to enable it." + ) + + return values