From 1bcd02e8fe37b46942ed38d338191ee7315acac7 Mon Sep 17 00:00:00 2001 From: George-Guryev-flxcmp Date: Thu, 4 Sep 2025 11:00:43 -0400 Subject: [PATCH 1/2] feat: implement structure extrusion when waveport is defined on a boundary --- CHANGELOG.md | 2 + .../test_terminal_component_modeler.py | 58 ++++- .../smatrix/component_modelers/terminal.py | 239 +++++++++++++++++- tidy3d/plugins/smatrix/ports/wave.py | 4 + 4 files changed, 301 insertions(+), 2 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 56d6dc97ac..b9a4096393 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -8,6 +8,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ## [Unreleased] ### Added +- Add automatic structure extrusion for waveports defined on boundaries, controlled by the `extrude_structures` field in `WavePort`. +- The extrusion method, implemented in `TerminalComponentModeler`, ensures that mode sources, absorbers, and PEC frames are fully contained within the extruded structures; extrusion occurs only when `extrude_structures` is set to `True`. - 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/tests/test_plugins/smatrix/test_terminal_component_modeler.py b/tests/test_plugins/smatrix/test_terminal_component_modeler.py index 420b06a4a9..bb2e60bd13 100644 --- a/tests/test_plugins/smatrix/test_terminal_component_modeler.py +++ b/tests/test_plugins/smatrix/test_terminal_component_modeler.py @@ -33,7 +33,10 @@ 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, +) mm = 1e3 @@ -1351,3 +1354,56 @@ 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) + 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 + sims = list(tcm.sim_dict.values()) + + # loop over simulations + for sim in sims: + # get injection axis that would be used to extrude structure + inj_axis = sim.sources[0].injection_axis + + # 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[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)) diff --git a/tidy3d/plugins/smatrix/component_modelers/terminal.py b/tidy3d/plugins/smatrix/component_modelers/terminal.py index 257f0dae6c..a5a2d46fbc 100644 --- a/tidy3d/plugins/smatrix/component_modelers/terminal.py +++ b/tidy3d/plugins/smatrix/component_modelers/terminal.py @@ -7,6 +7,7 @@ import numpy as np import pydantic.v1 as pd +from tidy3d import Box, ClipOperation, GeometryGroup, GridSpec, PolySlab, Structure from tidy3d.components.base import cached_property from tidy3d.components.boundary import BroadbandModeABCSpec from tidy3d.components.geometry.utils_2d import snap_coordinate_to_grid @@ -260,7 +261,9 @@ 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) - sim_dict[task_name] = sim_with_src + + # extrude structures if necessary and update simulation + sim_dict[task_name] = self._extrude_port_structures(sim_with_src) # Check final simulations for grid size at ports for _, sim in sim_dict.items(): @@ -472,5 +475,239 @@ def get_radiation_monitor_by_name(self, monitor_name: str) -> DirectivityMonitor return monitor raise Tidy3dKeyError(f"No radiation monitor named '{monitor_name}'.") + def get_antenna_metrics_data( + self, + port_amplitudes: Optional[dict[str, complex]] = None, + monitor_name: Optional[str] = None, + ) -> AntennaMetricsData: + """Calculate antenna parameters using superposition of fields from multiple port excitations. + + The method computes the radiated far fields and port excitation power wave amplitudes + for a superposition of port excitations, which can be used to analyze antenna radiation + characteristics. + + Parameters + ---------- + port_amplitudes : dict[str, complex] = None + Dictionary mapping port names to their desired excitation amplitudes, ``a``. For each port, + :math:`\\frac{1}{2}|a|^2` represents the incident power from that port into the system. + If ``None``, uses only the first port without any scaling of the raw simulation data. + When ``None`` is passed as a port amplitude, the raw simulation data is used for that port. + Note that in this method ``a`` represents the incident wave amplitude + using the power wave definition in [2]. + monitor_name : str = None + Name of the :class:`.DirectivityMonitor` to use for calculating far fields. + If None, uses the first monitor in `radiation_monitors`. + + Returns + ------- + :class:`.AntennaMetricsData` + Container with antenna parameters including directivity, gain, and radiation efficiency, + computed from the superposition of fields from all excited ports. + """ + # Use the first port as default if none specified + if port_amplitudes is None: + port_amplitudes = {self.ports[0].name: None} + + # Check port names, and create map from port to amplitude + port_dict = {} + for key in port_amplitudes.keys(): + port, _ = self.network_dict[key] + port_dict[port] = port_amplitudes[key] + # Get the radiation monitor, use first as default + # if none specified + if monitor_name is None: + rad_mon = self.radiation_monitors[0] + else: + rad_mon = self.get_radiation_monitor_by_name(monitor_name) + + # Create data arrays for holding the superposition of all port power wave amplitudes + f = list(rad_mon.freqs) + coords = {"f": f, "port": list(self.matrix_indices_monitor)} + a_sum = PortDataArray( + np.zeros((len(f), len(self.matrix_indices_monitor)), dtype=complex), coords=coords + ) + b_sum = a_sum.copy() + # Retrieve associated simulation data + combined_directivity_data = None + for port, amplitude in port_dict.items(): + if amplitude == 0.0: + continue + sim_data_port = self.batch_data[self._task_name(port=port)] + radiation_data = sim_data_port[rad_mon.name] + + a, b = self.compute_wave_amplitudes_at_each_port( + self.port_reference_impedances, sim_data_port, s_param_def="power" + ) + # Select a possible subset of frequencies + a = a.sel(f=f) + b = b.sel(f=f) + a_raw = a.sel(port=self.network_index(port)) + + if amplitude is None: + # No scaling performed when amplitude is None + scaled_directivity_data = sim_data_port[rad_mon.name] + scale_factor = 1.0 + else: + scaled_directivity_data = self._monitor_data_at_port_amplitude( + port, sim_data_port, radiation_data, amplitude + ) + scale_factor = amplitude / a_raw + a = scale_factor * a + b = scale_factor * b + + # Combine the possibly scaled directivity data and the power wave amplitudes + if combined_directivity_data is None: + combined_directivity_data = scaled_directivity_data + else: + combined_directivity_data = combined_directivity_data + scaled_directivity_data + a_sum += a + b_sum += b + + # Compute and add power measures to results + power_incident = np.real(0.5 * a_sum * np.conj(a_sum)).sum(dim="port") + power_reflected = np.real(0.5 * b_sum * np.conj(b_sum)).sum(dim="port") + return AntennaMetricsData.from_directivity_data( + combined_directivity_data, power_incident, power_reflected + ) + + 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 coordinated of the simulation grid + coords = sim.grid.boundaries.to_list + + mode_sources = [] + + # 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: + # update center here (example) + inj_axis = port.injection_axis + + port_center = list(port.center) + + idx = np.abs(port_center[inj_axis] - coords[inj_axis]).argmin() + port_center[inj_axis] = coords[inj_axis][idx] + + port = port.updated_copy(center=tuple(port_center)) + + mode_src_pos = port.center[port.injection_axis] + self._shift_value_signed(port) + + # then convert to source + mode_sources.append(port.to_source(self._source_time, snap_center=mode_src_pos)) + + # clip indices to a valid range + def _clip(i, lo, hi): + return int(max(lo, min(hi, i))) + + new_structures = [] + + # loop over individual mode sources associated with waveports + for mode in mode_sources: + direction = mode.direction + inj_axis = mode.injection_axis + span_indx = np.array(sim.grid.discretize_inds(mode)) + + target_val = mode.center[inj_axis] + + bnd_coords = coords[inj_axis] + + offset = mode.frame.length + sim.internal_absorbers[0].grid_shift + 1 + + # get total number of boundaries along injection direction + n_axis = len(bnd_coords) - 1 + + # define indicies of cells to be used for extrusion + if direction == "+": + idx = np.searchsorted(bnd_coords, target_val, side="left") - 1 + left = _clip(idx - 1, 0, n_axis) + right = _clip(idx + offset, 0, n_axis) + else: + idx = np.searchsorted(bnd_coords, target_val, side="right") + left = _clip(idx - offset, 0, n_axis) + right = _clip(idx + 1, 0, n_axis) + + # get indices for extrusion box boundaries + span_indx[inj_axis][0] = left + span_indx[inj_axis][1] = right + + # get bounding box bounds + box_bounds = [[c[beg], c[end]] for c, (beg, end) in zip(coords, span_indx)] + + # construct extrusion bounding box from bounds + box = Box.from_bounds(*np.transpose(box_bounds)) + + # get bounding box faces orthogonal to extrusion direction + slices = box.surfaces(box.size, box.center) + slice_plane_left = slices[2 * inj_axis] + slice_plane_right = slices[2 * inj_axis + 1] + + # 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 + left_geom = slice_plane_left.intersections_with(structure.geometry) + right_geom = slice_plane_right.intersections_with(structure.geometry) + shapely_geom = left_geom or right_geom or [] + + new_geoms = [] + + # loop over identified geometries and extrude them + for polygon in shapely_geom: + # construct outer shell of an extruded geometry first + exterior_vertices = np.array(polygon.exterior.coords) + outer_shell = PolySlab( + axis=inj_axis, slab_bounds=box_bounds[inj_axis], vertices=exterior_vertices + ) + + # construct innner shells that represent holes + hole_polyslabs = [ + PolySlab( + axis=inj_axis, + slab_bounds=box_bounds[inj_axis], + 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 + + new_geoms.append(extruded_slab_new) + + new_geoms.append(structure.geometry) + + new_struct = Structure( + geometry=GeometryGroup(geometries=new_geoms), medium=structure.medium + ) + new_structures.append(new_struct) + + # return simulation with added extruded structures + return sim.updated_copy(grid_spec=GridSpec.from_grid(sim.grid), structures=new_structures) + TerminalComponentModeler.update_forward_refs() + diff --git a/tidy3d/plugins/smatrix/ports/wave.py b/tidy3d/plugins/smatrix/ports/wave.py index 291cb5421c..ea31fec5e5 100644 --- a/tidy3d/plugins/smatrix/ports/wave.py +++ b/tidy3d/plugins/smatrix/ports/wave.py @@ -98,6 +98,10 @@ class WavePort(AbstractTerminalPort, Box): "If ``ABCBoundary`` or ``ModeABCBoundary``, a mode absorber is placed in the port with the specified boundary conditions.", ) + extrude_structures: bool = pd.Field( + False, title="Extrusion flag", description="Extrude structures attached to wave port." + ) + def _mode_voltage_coefficients(self, mode_data: ModeData) -> FreqModeDataArray: """Calculates scaling coefficients to convert mode amplitudes to the total port voltage. From 156ae248bc8b65bf1c72939335e97bc973704a7f Mon Sep 17 00:00:00 2001 From: George-Guryev-flxcmp Date: Thu, 4 Sep 2025 12:48:34 -0400 Subject: [PATCH 2/2] --amend --- CHANGELOG.md | 3 +- schemas/TerminalComponentModeler.json | 4 + .../smatrix/terminal_component_modeler_def.py | 144 +++++++++ .../test_terminal_component_modeler.py | 111 +++++-- tidy3d/components/geometry/utils.py | 3 +- tidy3d/components/simulation.py | 55 ++-- .../smatrix/component_modelers/terminal.py | 291 ++++++------------ tidy3d/plugins/smatrix/ports/wave.py | 20 +- 8 files changed, 372 insertions(+), 259 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index b9a4096393..78b11a32ac 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -8,8 +8,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ## [Unreleased] ### Added -- Add automatic structure extrusion for waveports defined on boundaries, controlled by the `extrude_structures` field in `WavePort`. -- The extrusion method, implemented in `TerminalComponentModeler`, ensures that mode sources, absorbers, and PEC frames are fully contained within the extruded structures; extrusion occurs only when `extrude_structures` is set to `True`. +- 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 bb2e60bd13..9a918ab847 100644 --- a/tests/test_plugins/smatrix/test_terminal_component_modeler.py +++ b/tests/test_plugins/smatrix/test_terminal_component_modeler.py @@ -36,6 +36,7 @@ from .terminal_component_modeler_def import ( make_coaxial_component_modeler, make_component_modeler, + make_differential_stripline_modeler, ) mm = 1e3 @@ -1370,40 +1371,102 @@ def test_wave_port_extrusion_coaxial(): 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 - sims = list(tcm.sim_dict.values()) + sim = tcm.base_sim - # loop over simulations - for sim in sims: - # get injection axis that would be used to extrude structure - inj_axis = sim.sources[0].injection_axis + # 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 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[0].slab_bounds), - np.max(sim.structures[2].geometry.geometries[0].slab_bounds), - ] + # 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) - pec_bnds = [] + # 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)] - # 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 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)] + # 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)) + # 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 a5a2d46fbc..f0941bd8f3 100644 --- a/tidy3d/plugins/smatrix/component_modelers/terminal.py +++ b/tidy3d/plugins/smatrix/component_modelers/terminal.py @@ -7,9 +7,10 @@ import numpy as np import pydantic.v1 as pd -from tidy3d import Box, ClipOperation, GeometryGroup, GridSpec, PolySlab, Structure +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 @@ -17,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 ( @@ -261,9 +262,8 @@ 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) - - # extrude structures if necessary and update simulation - sim_dict[task_name] = self._extrude_port_structures(sim_with_src) + # update simulation + sim_dict[task_name] = sim_with_src # Check final simulations for grid size at ports for _, sim in sim_dict.items(): @@ -347,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.""" @@ -366,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 @@ -475,102 +482,6 @@ def get_radiation_monitor_by_name(self, monitor_name: str) -> DirectivityMonitor return monitor raise Tidy3dKeyError(f"No radiation monitor named '{monitor_name}'.") - def get_antenna_metrics_data( - self, - port_amplitudes: Optional[dict[str, complex]] = None, - monitor_name: Optional[str] = None, - ) -> AntennaMetricsData: - """Calculate antenna parameters using superposition of fields from multiple port excitations. - - The method computes the radiated far fields and port excitation power wave amplitudes - for a superposition of port excitations, which can be used to analyze antenna radiation - characteristics. - - Parameters - ---------- - port_amplitudes : dict[str, complex] = None - Dictionary mapping port names to their desired excitation amplitudes, ``a``. For each port, - :math:`\\frac{1}{2}|a|^2` represents the incident power from that port into the system. - If ``None``, uses only the first port without any scaling of the raw simulation data. - When ``None`` is passed as a port amplitude, the raw simulation data is used for that port. - Note that in this method ``a`` represents the incident wave amplitude - using the power wave definition in [2]. - monitor_name : str = None - Name of the :class:`.DirectivityMonitor` to use for calculating far fields. - If None, uses the first monitor in `radiation_monitors`. - - Returns - ------- - :class:`.AntennaMetricsData` - Container with antenna parameters including directivity, gain, and radiation efficiency, - computed from the superposition of fields from all excited ports. - """ - # Use the first port as default if none specified - if port_amplitudes is None: - port_amplitudes = {self.ports[0].name: None} - - # Check port names, and create map from port to amplitude - port_dict = {} - for key in port_amplitudes.keys(): - port, _ = self.network_dict[key] - port_dict[port] = port_amplitudes[key] - # Get the radiation monitor, use first as default - # if none specified - if monitor_name is None: - rad_mon = self.radiation_monitors[0] - else: - rad_mon = self.get_radiation_monitor_by_name(monitor_name) - - # Create data arrays for holding the superposition of all port power wave amplitudes - f = list(rad_mon.freqs) - coords = {"f": f, "port": list(self.matrix_indices_monitor)} - a_sum = PortDataArray( - np.zeros((len(f), len(self.matrix_indices_monitor)), dtype=complex), coords=coords - ) - b_sum = a_sum.copy() - # Retrieve associated simulation data - combined_directivity_data = None - for port, amplitude in port_dict.items(): - if amplitude == 0.0: - continue - sim_data_port = self.batch_data[self._task_name(port=port)] - radiation_data = sim_data_port[rad_mon.name] - - a, b = self.compute_wave_amplitudes_at_each_port( - self.port_reference_impedances, sim_data_port, s_param_def="power" - ) - # Select a possible subset of frequencies - a = a.sel(f=f) - b = b.sel(f=f) - a_raw = a.sel(port=self.network_index(port)) - - if amplitude is None: - # No scaling performed when amplitude is None - scaled_directivity_data = sim_data_port[rad_mon.name] - scale_factor = 1.0 - else: - scaled_directivity_data = self._monitor_data_at_port_amplitude( - port, sim_data_port, radiation_data, amplitude - ) - scale_factor = amplitude / a_raw - a = scale_factor * a - b = scale_factor * b - - # Combine the possibly scaled directivity data and the power wave amplitudes - if combined_directivity_data is None: - combined_directivity_data = scaled_directivity_data - else: - combined_directivity_data = combined_directivity_data + scaled_directivity_data - a_sum += a - b_sum += b - - # Compute and add power measures to results - power_incident = np.real(0.5 * a_sum * np.conj(a_sum)).sum(dim="port") - power_reflected = np.real(0.5 * b_sum * np.conj(b_sum)).sum(dim="port") - return AntennaMetricsData.from_directivity_data( - combined_directivity_data, power_incident, power_reflected - ) - def _extrude_port_structures(self, sim: Simulation) -> Simulation: """ Extrude structures intersecting a port plane when a wave port lies on a structure boundary. @@ -591,123 +502,97 @@ def _extrude_port_structures(self, sim: Simulation) -> Simulation: Updated simulation with extruded structures added to ``simulation.structures``. """ - # get coordinated of the simulation grid - coords = sim.grid.boundaries.to_list - - mode_sources = [] - # 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: - # update center here (example) - inj_axis = port.injection_axis - - port_center = list(port.center) - - idx = np.abs(port_center[inj_axis] - coords[inj_axis]).argmin() - port_center[inj_axis] = coords[inj_axis][idx] - - port = port.updated_copy(center=tuple(port_center)) + # 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, + ) - mode_src_pos = port.center[port.injection_axis] + self._shift_value_signed(port) + # 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) - # then convert to source - mode_sources.append(port.to_source(self._source_time, snap_center=mode_src_pos)) + # 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)] - # clip indices to a valid range - def _clip(i, lo, hi): - return int(max(lo, min(hi, i))) + # get extrusion extent along injection axis + extrude_to = cutting_plane.center[inj_axis] - new_structures = [] + # 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) - # loop over individual mode sources associated with waveports - for mode in mode_sources: - direction = mode.direction - inj_axis = mode.injection_axis - span_indx = np.array(sim.grid.discretize_inds(mode)) + # define extrusion bounds + extrusion_bounds = [center[inj_axis], extrude_to][::sign] - target_val = mode.center[inj_axis] + new_structures = [] - bnd_coords = coords[inj_axis] + # 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 - offset = mode.frame.length + sim.internal_absorbers[0].grid_shift + 1 + shapely_geom = cutting_plane.intersections_with(structure.geometry) - # get total number of boundaries along injection direction - n_axis = len(bnd_coords) - 1 + polygon_list = [] + for geom in shapely_geom: + polygon_list = polygon_list + ClipOperation.to_polygon_list(geom) - # define indicies of cells to be used for extrusion - if direction == "+": - idx = np.searchsorted(bnd_coords, target_val, side="left") - 1 - left = _clip(idx - 1, 0, n_axis) - right = _clip(idx + offset, 0, n_axis) - else: - idx = np.searchsorted(bnd_coords, target_val, side="right") - left = _clip(idx - offset, 0, n_axis) - right = _clip(idx + 1, 0, n_axis) - - # get indices for extrusion box boundaries - span_indx[inj_axis][0] = left - span_indx[inj_axis][1] = right - - # get bounding box bounds - box_bounds = [[c[beg], c[end]] for c, (beg, end) in zip(coords, span_indx)] - - # construct extrusion bounding box from bounds - box = Box.from_bounds(*np.transpose(box_bounds)) - - # get bounding box faces orthogonal to extrusion direction - slices = box.surfaces(box.size, box.center) - slice_plane_left = slices[2 * inj_axis] - slice_plane_right = slices[2 * inj_axis + 1] - - # 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 - left_geom = slice_plane_left.intersections_with(structure.geometry) - right_geom = slice_plane_right.intersections_with(structure.geometry) - shapely_geom = left_geom or right_geom or [] - - new_geoms = [] - - # loop over identified geometries and extrude them - for polygon in shapely_geom: - # construct outer shell of an extruded geometry first - exterior_vertices = np.array(polygon.exterior.coords) - outer_shell = PolySlab( - axis=inj_axis, slab_bounds=box_bounds[inj_axis], vertices=exterior_vertices - ) - - # construct innner shells that represent holes - hole_polyslabs = [ - PolySlab( - axis=inj_axis, - slab_bounds=box_bounds[inj_axis], - 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 + 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 ) - else: - extruded_slab_new = outer_shell - new_geoms.append(extruded_slab_new) - - new_geoms.append(structure.geometry) + # 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) - new_struct = Structure( - geometry=GeometryGroup(geometries=new_geoms), medium=structure.medium + # update structures in simulation while keeping the same grid + sim = sim.updated_copy( + grid_spec=GridSpec.from_grid(sim.grid), structures=new_structures ) - new_structures.append(new_struct) - # return simulation with added extruded structures - return sim.updated_copy(grid_spec=GridSpec.from_grid(sim.grid), structures=new_structures) - + return sim -TerminalComponentModeler.update_forward_refs() +TerminalComponentModeler.update_forward_refs() diff --git a/tidy3d/plugins/smatrix/ports/wave.py b/tidy3d/plugins/smatrix/ports/wave.py index ea31fec5e5..5077f6ba01 100644 --- a/tidy3d/plugins/smatrix/ports/wave.py +++ b/tidy3d/plugins/smatrix/ports/wave.py @@ -87,19 +87,21 @@ 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="Extrusion flag", description="Extrude structures attached to wave port." + 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: @@ -324,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