fixed bug when coaxial ports are snapped to grid cell boundaries. add…

…ed additional tests to both rectangular and coaxial type lumped ports

CHANGELOG.md

@@ -21,8 +21,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 - Bug where boundary layers would be plotted too small in 2D simulations.
 - Significant speedup for field projection computations.
 - Fix numerical precision issue in `FieldProjectionCartesianMonitor`.
-- 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
 

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,15 +107,15 @@ 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:
         strip_length = length
     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)

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,13 +266,27 @@ 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
     )
 
     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)

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):

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,16 +313,19 @@ 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
         circle_vertices[-1, :] = circle_vertices[0, :]
 
         # 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: