From 87010c8a97f91a0218fade8a04a717301e49002d Mon Sep 17 00:00:00 2001 From: dmarek Date: Mon, 24 Jun 2024 10:00:43 -0400 Subject: [PATCH] fixed bug when coaxial ports are snapped to grid cell boundaries. added additional tests to both rectangular and coaxial type lumped ports --- CHANGELOG.md | 2 +- .../terminal_component_modeler_def.py | 52 ++++++-------- .../test_terminal_component_modeler.py | 67 +++++++++++++++---- tidy3d/plugins/smatrix/ports/base_lumped.py | 14 +++- .../plugins/smatrix/ports/coaxial_lumped.py | 36 +++++----- 5 files changed, 110 insertions(+), 61 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index bc267e81f..f1aac8974 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -16,7 +16,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### Fixed - Bug where boundary layers would be plotted too small in 2D simulations. - Bug when plotting transformed geometries. -- Bug when placing path integrals in the `CoaxialLumpedPort`. +- Bug when snapping `CoaxialLumpedPort` to grid cell boundaries. ## [2.7.0] - 2024-06-17 diff --git a/tests/test_plugins/terminal_component_modeler_def.py b/tests/test_plugins/terminal_component_modeler_def.py index de762b810..2928dc632 100644 --- a/tests/test_plugins/terminal_component_modeler_def.py +++ b/tests/test_plugins/terminal_component_modeler_def.py @@ -7,12 +7,12 @@ ) # Microstrip dimensions -unit = 1e6 -default_strip_length = 75e-3 * unit -strip_width = 3e-3 * unit -gap = 1e-3 * unit +mm = 1e3 +default_strip_length = 75 * mm +strip_width = 3 * mm +gap = 1 * mm gnd_width = strip_width * 8 -metal_thickness = 0.2e-3 * unit +metal_thickness = 0.2 * mm # Microstrip materials pec = td.PECMedium() @@ -25,11 +25,11 @@ freq_stop = 10e9 # Coaxial dimensions -Rinner = 0.2768 * 1e-3 -Router = 1.0 * 1e-3 +Rinner = 0.2768 * mm +Router = 1.0 * mm -def make_simulation(planar_pec: bool, length: float = None, auto_grid: bool = True): +def make_simulation(planar_pec: bool, length: float = None, grid_spec: td.GridSpec = None): if length: strip_length = length else: @@ -49,10 +49,8 @@ def make_simulation(planar_pec: bool, length: float = None, auto_grid: bool = Tr run_time = 60 / fwidth # Spatial grid specification - if auto_grid: + if not grid_spec: grid_spec = td.GridSpec.auto(min_steps_per_wvl=10, wavelength=td.C_0 / freq_stop) - else: - grid_spec = td.GridSpec.uniform(wavelength0 / 11) # Make structures strip = td.Structure( @@ -109,7 +107,7 @@ def make_component_modeler( reference_impedance: complex = 50, length: float = None, port_refinement: bool = True, - auto_grid: bool = True, + grid_spec: td.GridSpec = None, **kwargs, ): if length: @@ -117,7 +115,7 @@ def make_component_modeler( else: strip_length = default_strip_length - sim = make_simulation(planar_pec, length=length, auto_grid=auto_grid) + sim = make_simulation(planar_pec, length=length, grid_spec=grid_spec) if planar_pec: height = 0 @@ -162,11 +160,9 @@ def make_component_modeler( return modeler -def make_coaxial_simulation(length: float = None, auto_grid: bool = True): - if length: - coax_length = length - else: - coax_length = default_strip_length +def make_coaxial_simulation(length: float = None, grid_spec: td.GridSpec = None): + if not length: + length = default_strip_length # wavelength / frequency freq0 = (freq_start + freq_stop) / 2 @@ -175,10 +171,8 @@ def make_coaxial_simulation(length: float = None, auto_grid: bool = True): run_time = 60 / fwidth # Spatial grid specification - if auto_grid: + if not grid_spec: grid_spec = td.GridSpec.auto(min_steps_per_wvl=10, wavelength=td.C_0 / freq_stop) - else: - grid_spec = td.GridSpec.uniform(wavelength0 / 11) # Make structures inner_conductor = td.Cylinder( @@ -223,7 +217,7 @@ def make_coaxial_simulation(length: float = None, auto_grid: bool = True): size_sim = [ 2 * Router, 2 * Router, - coax_length + 0.5 * wavelength0, + length + 0.5 * wavelength0, ] sim = td.Simulation( @@ -245,17 +239,15 @@ def make_coaxial_component_modeler( reference_impedance: complex = 50, length: float = None, port_refinement: bool = True, - auto_grid: bool = True, + grid_spec: td.GridSpec = None, **kwargs, ): - if length: - coax_length = length - else: - coax_length = default_strip_length + if not length: + length = default_strip_length - sim = make_coaxial_simulation(length=coax_length, auto_grid=auto_grid) + sim = make_coaxial_simulation(length=length, grid_spec=grid_spec) - center_src1 = [0, 0, -coax_length / 2] + center_src1 = [0, 0, -length / 2] port_cells = None if port_refinement: @@ -271,7 +263,7 @@ def make_coaxial_component_modeler( num_grid_cells=port_cells, impedance=reference_impedance, ) - center_src2 = [0, 0, coax_length / 2] + center_src2 = [0, 0, length / 2] port_2 = port_1.updated_copy(name="coax_port_2", center=center_src2, direction="-") ports = [port_1, port_2] freqs = np.linspace(freq_start, freq_stop, 100) diff --git a/tests/test_plugins/test_terminal_component_modeler.py b/tests/test_plugins/test_terminal_component_modeler.py index 05ba24c4c..167a1f54e 100644 --- a/tests/test_plugins/test_terminal_component_modeler.py +++ b/tests/test_plugins/test_terminal_component_modeler.py @@ -11,6 +11,7 @@ LumpedPortDataArray, TerminalComponentModeler, ) +from tidy3d.plugins.smatrix.ports.base_lumped import AbstractLumpedPort from ..utils import run_emulated from .terminal_component_modeler_def import make_coaxial_component_modeler, make_component_modeler @@ -28,6 +29,31 @@ def run_component_modeler(monkeypatch, modeler: TerminalComponentModeler): return s_matrix +def check_lumped_port_components_snapped_correctly(modeler: TerminalComponentModeler): + """Given an instance of a ``TerminalComponentModeler``, check that all simulation components + have been snapped exactly to the position of the load resistor. + """ + sim_dict = modeler.sim_dict + num_ports = len(modeler.ports) + # Check to make sure all components are exactly aligned along the normal axis + for src_port, src_idx, src_sim in zip(modeler.ports, range(num_ports), sim_dict.values()): + assert isinstance(src_port, AbstractLumpedPort) + monitor_dict = {monitor.name: monitor for monitor in src_sim.monitors} + normal_axis = src_port.injection_axis + center_load = src_sim.lumped_elements[src_idx].center[normal_axis] + assert len(src_sim.sources) == 1 + center_source = src_sim.sources[0].center[normal_axis] + assert center_load == center_source + for port, idx in zip(modeler.ports, range(num_ports)): + assert isinstance(port, AbstractLumpedPort) + normal_axis = port.injection_axis + center_load = src_sim.lumped_elements[idx].center[normal_axis] + center_voltage_monitor = monitor_dict[port._voltage_monitor_name].center[normal_axis] + center_current_monitor = monitor_dict[port._current_monitor_name].center[normal_axis] + assert center_load == center_voltage_monitor + assert center_load == center_current_monitor + + def test_validate_no_sources(tmp_path): modeler = make_component_modeler(planar_pec=True, path_dir=str(tmp_path)) source = td.PointDipole( @@ -145,13 +171,17 @@ def test_ab_to_s_component_modeler(): assert np.isclose(S_matrix, b_matrix).all() -def test_port_snapping(monkeypatch, tmp_path): +def test_port_snapping(tmp_path): + """Make sure that the snapping behavior of the load resistor is mirrored + by all other components in the modeler simulations with rectangular ports. + """ + y_z_grid = td.UniformGrid(dl=0.1 * 1e3) + x_grid = td.UniformGrid(dl=11 * 1e3) + grid_spec = td.GridSpec(grid_x=x_grid, grid_y=y_z_grid, grid_z=y_z_grid) modeler = make_component_modeler( - planar_pec=True, path_dir=str(tmp_path), port_refinement=False, auto_grid=False + planar_pec=True, path_dir=str(tmp_path), port_refinement=False, grid_spec=grid_spec ) - # Without port refinement the grid is much too coarse for these port sizes - with pytest.raises(SetupError): - _ = run_component_modeler(monkeypatch, modeler) + check_lumped_port_components_snapped_correctly(modeler=modeler) def test_coarse_grid_at_port(monkeypatch, tmp_path): @@ -187,9 +217,7 @@ def test_run_coaxial_component_modeler(monkeypatch, tmp_path): def test_coarse_grid_at_coaxial_port(monkeypatch, tmp_path): - modeler = make_coaxial_component_modeler( - path_dir=str(tmp_path), port_refinement=False, auto_grid=False - ) + modeler = make_coaxial_component_modeler(path_dir=str(tmp_path), port_refinement=False) # Without port refinement the grid is much too coarse for these port sizes with pytest.raises(SetupError): _ = run_component_modeler(monkeypatch, modeler) @@ -226,7 +254,8 @@ def test_validate_coaxial_port_diameters(): @pytest.mark.parametrize("direction", ["+", "-"]) def test_current_integral_positioning_coaxial_port(direction): """Make sure the positioning of the current integral used by the CoaxialLumpedPort is correct, - when the coordinates and port position do not exactly match. + when the coordinates and port position do not exactly match. This requires that the port is + snapped correctly to cell boundaries. """ # Test coordinates from a failing case normal_coords = np.array( @@ -237,8 +266,9 @@ def test_current_integral_positioning_coaxial_port(direction): -14009.999999999978, ] ) - - normal_port_position = -14030 + # The port center should be snapped to cell boundaries which is the midpoint of + # adjacent transverse magnetic field locations + normal_port_position = (normal_coords[2] + normal_coords[3]) / 2 path_pos = CoaxialLumpedPort._determine_current_integral_pos( normal_port_position, normal_coords, direction ) @@ -246,4 +276,17 @@ def test_current_integral_positioning_coaxial_port(direction): if direction == "+": assert path_pos == normal_coords[3] else: - assert path_pos == normal_coords[1] + assert path_pos == normal_coords[2] + + +def test_coaxial_port_snapping(tmp_path): + """Make sure that the snapping behavior of the load resistor is mirrored + by all other components in the modeler simulations with coaxial ports. + """ + x_y_grid = td.UniformGrid(dl=0.1 * 1e3) + z_grid = td.UniformGrid(dl=11 * 1e3) + grid_spec = td.GridSpec(grid_x=x_y_grid, grid_y=x_y_grid, grid_z=z_grid) + modeler = make_coaxial_component_modeler( + path_dir=str(tmp_path), port_refinement=False, grid_spec=grid_spec + ) + check_lumped_port_components_snapped_correctly(modeler=modeler) diff --git a/tidy3d/plugins/smatrix/ports/base_lumped.py b/tidy3d/plugins/smatrix/ports/base_lumped.py index a353fcc8b..e921d5488 100644 --- a/tidy3d/plugins/smatrix/ports/base_lumped.py +++ b/tidy3d/plugins/smatrix/ports/base_lumped.py @@ -8,11 +8,12 @@ from ....components.base import Tidy3dBaseModel, cached_property from ....components.data.data_array import DataArray, FreqDataArray from ....components.data.sim_data import SimulationData +from ....components.geometry.utils_2d import snap_coordinate_to_grid from ....components.grid.grid import Grid, YeeGrid from ....components.lumped_element import AbstractLumpedResistor from ....components.monitor import FieldMonitor from ....components.source import GaussianPulse, UniformCurrentSource -from ....components.types import Complex, FreqArray +from ....components.types import Complex, Coordinate, FreqArray from ....constants import OHM DEFAULT_PORT_NUM_CELLS = 3 @@ -74,6 +75,17 @@ def _voltage_monitor_name(self) -> str: def _current_monitor_name(self) -> str: return f"{self.name}_current" + def snapped_center(self, grid: Grid) -> Coordinate: + """Get the exact center of this port after snapping along the injection axis. + Ports are snapped to the nearest Yee cell boundary to match the exact position + of the ``AbstractLumpedResistor". + """ + center = list(self.center) + normal_axis = self.injection_axis + normal_port_center = center[normal_axis] + center[normal_axis] = snap_coordinate_to_grid(grid, normal_port_center, normal_axis) + return tuple(center) + @cached_property @abstractmethod def injection_axis(self): diff --git a/tidy3d/plugins/smatrix/ports/coaxial_lumped.py b/tidy3d/plugins/smatrix/ports/coaxial_lumped.py index f892e9ea7..fdbe1d31d 100644 --- a/tidy3d/plugins/smatrix/ports/coaxial_lumped.py +++ b/tidy3d/plugins/smatrix/ports/coaxial_lumped.py @@ -178,7 +178,7 @@ def compute_coax_current(rin, rout, x, y): dataset_E = FieldDataset(**kwargs) return CustomCurrentSource( - center=self.center, + center=center, size=(self.outer_diameter, self.outer_diameter, 0), source_time=source_time, name=self.name, @@ -218,7 +218,7 @@ def to_voltage_monitor(self, freqs: FreqArray, snap_center: float) -> FieldMonit # Create a voltage monitor return FieldMonitor( - center=self._voltage_path_center, + center=self._voltage_path_center(center), size=self._voltage_path_size, freqs=freqs, fields=[E1, E2], @@ -261,11 +261,11 @@ def compute_voltage(self, sim_data: SimulationData) -> FreqDataArray: We arbitrarily choose the positive first in-plane axis as the location for the path. Any of the four possible choices should give the same result. """ - + exact_port_center = self.snapped_center(sim_data.simulation.grid) field_data = sim_data[self._voltage_monitor_name] voltage_integral = VoltageIntegralAxisAligned( - center=self._voltage_path_center, + center=self._voltage_path_center(exact_port_center), size=self._voltage_path_size, extrapolate_to_endpoints=True, snap_path_to_grid=True, @@ -281,7 +281,7 @@ def compute_current(self, sim_data: SimulationData) -> FreqDataArray: The contour is a closed loop around the inner conductor. It is positioned at the midpoint between inner and outer radius of the annulus. """ - + exact_port_center = self.snapped_center(sim_data.simulation.grid) # Loops around inner conductive circle conductor field_data = sim_data[self._current_monitor_name] @@ -313,8 +313,8 @@ def generate_circle_coordinates(radius, num_points): xt, yt = generate_circle_coordinates( (self.outer_diameter + self.inner_diameter) / 4, num_path_coords ) - xt += self.center[trans_axes[0]] - yt += self.center[trans_axes[1]] + xt += exact_port_center[trans_axes[0]] + yt += exact_port_center[trans_axes[1]] circle_vertices = np.column_stack((xt, yt)) # Close the contour exactly @@ -322,7 +322,10 @@ def generate_circle_coordinates(radius, num_points): # Get the coordinates normal to port and select positions just on either side of the port normal_coords = field_coords[coord3].values - normal_port_position = self.center[self.injection_axis] + normal_port_position = exact_port_center[self.injection_axis] + # The exact center position of the port should coincide with a Yee cell boundary, so we + # want to select magnetic field positions a half-step on either side, + # depending on the direction. path_pos = CoaxialLumpedPort._determine_current_integral_pos( normal_port_position, normal_coords, self.direction ) @@ -341,31 +344,30 @@ def generate_circle_coordinates(radius, num_points): @staticmethod def _determine_current_integral_pos( - port_center: float, normal_coords: np.array, direction: Direction + snapped_center: float, normal_coords: np.array, direction: Direction ) -> float: """Helper to locate where the current integral should be placed in relation to the normal axis of the port. """ - idx = np.abs(normal_coords - port_center).argmin() + upper_bound = np.searchsorted(normal_coords, snapped_center) + lower_bound = upper_bound - 1 # We need to choose which side of the port to place the path integral, - # which depends on which way the inner conductor is extending if direction == "+": - return normal_coords[idx + 1] + return normal_coords[upper_bound] else: - return normal_coords[idx - 1] + return normal_coords[lower_bound] @cached_property def _voltage_axis(self) -> Axis: return self.remaining_axes[0] - @cached_property - def _voltage_path_center(self) -> Coordinate: + def _voltage_path_center(self, port_center: Coordinate) -> Coordinate: """We arbitrarily choose the positive first in-plane axis as the location for the path. Any of the four possible choices should give the same result. """ - center = list(self.center) + center = list(port_center) center[self._voltage_axis] += (self.outer_diameter + self.inner_diameter) / 4 - return center + return tuple(center) @cached_property def _voltage_path_size(self) -> Size: