fix bug in port snapping and correct bounds for current integral

added grid size check to TerminalComponentModeler and improved test coverage of smatrix plugin

Author: dmarek-flex

tests/test_plugins/terminal_component_modeler_def.py

@@ -25,7 +25,7 @@
 freq_stop = 10e9
 
 
-def make_simulation(planar_pec: bool, length: float = None):
+def make_simulation(planar_pec: bool, length: float = None, auto_grid: bool = True):
     if length:
         strip_length = length
     else:
@@ -45,7 +45,10 @@ def make_simulation(planar_pec: bool, length: float = None):
     run_time = 60 / fwidth
 
     # Spatial grid specification
-    grid_spec = td.GridSpec.auto(min_steps_per_wvl=10, wavelength=td.C_0 / freq_stop)
+    if auto_grid:
+        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(
@@ -98,14 +101,19 @@ def make_simulation(planar_pec: bool, length: float = None):
 
 
 def make_component_modeler(
-    planar_pec: bool, reference_impedance: complex = 50, length: float = None, **kwargs
+    planar_pec: bool,
+    reference_impedance: complex = 50,
+    length: float = None,
+    port_refinement: bool = True,
+    auto_grid: bool = True,
+    **kwargs
 ):
     if length:
         strip_length = length
     else:
         strip_length = default_strip_length
 
-    sim = make_simulation(planar_pec, length=length)
+    sim = make_simulation(planar_pec, length=length, auto_grid=auto_grid)
 
     if planar_pec:
         height = 0
@@ -118,7 +126,9 @@ def make_component_modeler(
     center_src2 = [strip_length / 2, 0, height + gap / 2]
     size_src2 = [0, strip_width, gap]
 
-    port_cells = np.ceil(gap / (metal_thickness / 1))
+    port_cells = None
+    if port_refinement:
+        port_cells = np.ceil(gap / (metal_thickness / 1))
 
     port_1 = LumpedPort(
         center=center_src1,

tests/test_plugins/test_terminal_component_modeler.py

@@ -6,10 +6,11 @@
 import tidy3d as td
 from tidy3d.plugins.smatrix import (
     AbstractComponentModeler,
+    LumpedPort,
     LumpedPortDataArray,
     TerminalComponentModeler,
 )
-from tidy3d.exceptions import Tidy3dKeyError
+from tidy3d.exceptions import Tidy3dKeyError, SetupError
 from ..utils import run_emulated
 from .terminal_component_modeler_def import make_component_modeler
 
@@ -55,8 +56,11 @@ def test_plot_sim_eps(tmp_path):
     plt.close()
 
 
-def test_make_component_modeler(tmp_path):
-    _ = make_component_modeler(planar_pec=False, path_dir=str(tmp_path))
+@pytest.mark.parametrize("port_refinement", [False, True])
+def test_make_component_modeler(tmp_path, port_refinement):
+    _ = make_component_modeler(
+        planar_pec=False, path_dir=str(tmp_path), port_refinement=port_refinement
+    )
 
 
 def test_run(monkeypatch, tmp_path):
@@ -138,3 +142,24 @@ def test_ab_to_s_component_modeler():
     b_matrix = LumpedPortDataArray(data=b_values, coords=coords)
     S_matrix = AbstractComponentModeler.ab_to_s(a_matrix, b_matrix)
     assert np.isclose(S_matrix, b_matrix).all()
+
+
+def test_port_snapping(monkeypatch, tmp_path):
+    modeler = make_component_modeler(
+        planar_pec=True, path_dir=str(tmp_path), port_refinement=False, auto_grid=False
+    )
+    # Without port refinement the grid is much too coarse for these port sizes
+    with pytest.raises(SetupError):
+        _ = run_component_modeler(monkeypatch, modeler)
+
+
+def test_coarse_grid_at_port(monkeypatch, tmp_path):
+    modeler = make_component_modeler(planar_pec=True, 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)
+
+
+def test_validate_port_voltage_axis():
+    with pytest.raises(pydantic.ValidationError):
+        LumpedPort(center=(0, 0, 0), size=(0, 1, 2), voltage_axis=0, impedance=50)

tidy3d/plugins/smatrix/component_modelers/terminal.py

@@ -9,6 +9,7 @@
 import xarray as xr
 
 from ....constants import C_0, fp_eps
+from ....exceptions import SetupError
 from ....components.simulation import Simulation
 from ....components.geometry.utils_2d import snap_coordinate_to_grid
 from ....components.data.sim_data import SimulationData
@@ -148,7 +149,7 @@ def sim_dict(self) -> Dict[str, Simulation]:
                 list(self.simulation.monitors) + port_voltage_monitors + port_current_monitors
             )
             for port, center in zip(self.ports, snap_centers):
-                port_source = self.to_source(self._source_time, port=port, snap_center=center)
+                port_source = port.to_source(self._source_time, snap_center=center)
                 update_dict = dict(
                     sources=[port_source],
                     monitors=new_mnts,
@@ -160,6 +161,9 @@ def sim_dict(self) -> Dict[str, Simulation]:
                 task_name = self._task_name(port=port)
                 sim_dict[task_name] = sim_copy
 
+        # Check final simulations for grid size at ports
+        for _, sim in sim_dict.items():
+            TerminalComponentModeler._check_grid_size_at_ports(sim, self.ports)
         return sim_dict
 
     @cached_property
@@ -189,6 +193,15 @@ def _construct_smatrix(self, batch_data: BatchData) -> LumpedPortDataArray:
         a_matrix = LumpedPortDataArray(values, coords=coords)
         b_matrix = a_matrix.copy(deep=True)
 
+        def select_within_bounds(coords: np.array, min, max):
+            """Helper to slice coordinates within min and max bounds, including a tolerance."""
+            """xarray does not have this functionality yet. """
+            np_idx = np.logical_and(coords > min, coords < max)
+            np_idx = np.logical_or(np_idx, np.isclose(coords, min, rtol=fp_eps, atol=fp_eps))
+            np_idx = np.logical_or(np_idx, np.isclose(coords, max, rtol=fp_eps, atol=fp_eps))
+            coords = coords[np_idx]
+            return coords
+
         def port_voltage(port: LumpedPort, sim_data: SimulationData) -> xr.DataArray:
             """Helper to compute voltage across the port."""
             e_component = "xyz"[port.voltage_axis]
@@ -263,13 +276,21 @@ def port_current(port: LumpedPort, sim_data: SimulationData) -> xr.DataArray:
             h1_field = h_field.isel({orth_component: orth_index - 1})
             h2_field = h_field.isel({orth_component: orth_index})
             h_field = h1_field - h2_field
+
+            # Find exact bounds of port taking into consideration the Yee grid
+            np_coords = h_coords[port.current_axis].values
+            port_min = port.bounds[0][port.current_axis]
+            port_max = port.bounds[1][port.current_axis]
+            np_coords = select_within_bounds(np_coords, port_min, port_max)
+            coord_low = np_coords[0]
+            coord_high = np_coords[-1]
             # Extract cap field which is coincident with sheet
             h_cap = h_cap_field.isel({orth_component: orth_index})
 
             # Need to make sure to use the nearest coordinate that is
             # at least greater than the port bounds
-            hcap_minus = h_cap.sel({h_component: slice(-np.inf, port.bounds[0][port.current_axis])})
-            hcap_plus = h_cap.sel({h_component: slice(port.bounds[1][port.current_axis], np.inf)})
+            hcap_minus = h_cap.sel({h_component: slice(-np.inf, coord_low)})
+            hcap_plus = h_cap.sel({h_component: slice(coord_high, np.inf)})
             hcap_minus = hcap_minus.isel({h_component: -1})
             hcap_plus = hcap_plus.isel({h_component: 0})
             # Length of integration along the h_cap contour is a single cell width
@@ -344,3 +365,21 @@ def _validate_3d_simulation(cls, val):
                 f"'{cls.__name__}' must be setup with a 3D simulation with all sizes greater than 0."
             )
         return val
+
+    @staticmethod
+    def _check_grid_size_at_ports(simulation: Simulation, ports: list[LumpedPort]):
+        """Raises :class:`SetupError` if the grid is too coarse at port locations"""
+        yee_grid = simulation.grid.yee
+        for port in ports:
+            e_component = "xyz"[port.voltage_axis]
+            e_yee_grid = yee_grid.grid_dict[f"E{e_component}"]
+            coords = e_yee_grid.to_dict[e_component]
+            min_bound = port.bounds[0][port.voltage_axis]
+            max_bound = port.bounds[1][port.voltage_axis]
+            coords_within_port = np.any(np.logical_and(coords > min_bound, coords < max_bound))
+            if not coords_within_port:
+                raise SetupError(
+                    f"Grid is too coarse along '{e_component}' direction for the lumped port "
+                    f"at location {port.center}. Either set the port's 'num_grid_cells' to "
+                    f"a nonzero integer or modify the 'GridSpec'. "
+                )