add a TerminalSpec for RF specific mode information

adding support for colocated fields

fixing bug for more complicated shapes

adding tests, fixing symmetric conditions, fixing for 2D structures

fix simulation validator for terminal spec

refactor by splitting the path integral specification away from the integral computation, now the path specification handles pretty much everything, except for the final integral computation and results preparation

rename auto patch spec to generator type name, fixed tests

add test for sim validation and improved doc strings

fix regression

fix python 3.9 tests

rebase on ruff changes

Author: dmarek-flex

tests/test_components/test_microwave.py

@@ -3645,3 +3645,54 @@ def test_create_sim_multiphysics_with_incompatibilities():
                 ),
             ],
         )
+
+
+def test_validate_terminal_spec_generation():
+    """Test that auto generation of path specs is correctly validated for currently unsupported structures."""
+    freq0 = 10e9
+    mm = 1e3
+    run_time_spec = td.RunTimeSpec(quality_factor=3.0)
+    size = (10 * mm, 10 * mm, 10 * mm)
+    size_mon = (0, 8 * mm, 8 * mm)
+
+    # Currently limited to generation of axis aligned boxes around conductors,
+    # so the path may intersect other nearby conductors, like in this coaxial cable
+    coaxial = td.Structure(
+        geometry=td.GeometryGroup(
+            geometries=(
+                td.ClipOperation(
+                    operation="difference",
+                    geometry_a=td.Cylinder(
+                        axis=0, radius=2.5 * mm, center=(0, 0, 0), length=td.inf
+                    ),
+                    geometry_b=td.Cylinder(
+                        axis=0, radius=1.3 * mm, center=(0, 0, 0), length=td.inf
+                    ),
+                ),
+                td.Cylinder(axis=0, radius=1 * mm, center=(0, 0, 0), length=td.inf),
+            )
+        ),
+        medium=td.PEC,
+    )
+    mode_spec = td.ModeSpec(num_modes=2, target_neff=1.8, terminal_spec=td.TerminalSpec())
+
+    mode_mon = td.ModeMonitor(
+        center=(0, 0, 0),
+        size=size_mon,
+        freqs=[freq0],
+        name="mode_1",
+        colocate=False,
+        mode_spec=mode_spec,
+    )
+    sim = td.Simulation(
+        run_time=run_time_spec,
+        size=size,
+        sources=[],
+        structures=[coaxial],
+        grid_spec=td.GridSpec.uniform(dl=0.1 * mm),
+        monitors=[mode_mon],
+    )
+
+    # check that validation error is caught
+    with pytest.raises(SetupError):
+        sim._validate_terminal_specs()

tests/test_components/test_simulation.py

@@ -19,6 +19,13 @@
 
 from ..utils import get_spatial_coords_dict, run_emulated
 
+MAKE_PLOTS = False
+if MAKE_PLOTS:
+    # Interative plotting for debugging
+    from matplotlib import use
+
+    use("TkAgg")
+
 # Using similar code as "test_data/test_data_arrays.py"
 MON_SIZE = (2, 1, 0)
 FIELDS = ("Ex", "Ey", "Hx", "Hy")
@@ -785,13 +792,11 @@ def test_lobe_measurements(apply_cyclic_extension, include_endpoint):
 @pytest.mark.parametrize("min_value", [0.0, 1.0])
 def test_lobe_plots(min_value):
     """Run the lobe measurer on some test data and plot the results."""
-    # Interative plotting for debugging
-    # from matplotlib import use
-    # use("TkAgg")
     theta = np.linspace(0, 2 * np.pi, 301)
     Urad = np.cos(theta) ** 2 * np.cos(3 * theta) ** 2 + min_value
     lobe_measurer = mw.LobeMeasurer(angle=theta, radiation_pattern=Urad)
     _, ax = plt.subplots(1, 1, subplot_kw={"projection": "polar"})
     ax.plot(theta, Urad, "k")
     lobe_measurer.plot(0, ax)
-    plt.show()
+    if MAKE_PLOTS:
+        plt.show()

tests/test_plugins/test_microwave.py

@@ -17,6 +17,16 @@
 from tidy3d.components.microwave.data.monitor_data import (
     AntennaMetricsData,
 )
+from tidy3d.components.microwave.path_spec import (
+    CompositeCurrentIntegralSpec,
+    CurrentIntegralAxisAlignedSpec,
+    CustomCurrentIntegral2DSpec,
+    CustomVoltageIntegral2DSpec,
+    VoltageIntegralAxisAlignedSpec,
+)
+from tidy3d.components.microwave.terminal_spec import (
+    TerminalSpec,
+)
 from tidy3d.components.spice.analysis.dc import (
     ChargeToleranceSpec,
     IsothermalSteadyChargeDCAnalysis,
@@ -439,6 +449,7 @@ def set_logging_level(level: str) -> None:
     "ChargeToleranceSpec",
     "ClipOperation",
     "CoaxialLumpedResistor",
+    "CompositeCurrentIntegralSpec",
     "ConstantDoping",
     "ConstantMobilityModel",
     "ContinuousWave",
@@ -449,8 +460,10 @@ def set_logging_level(level: str) -> None:
     "Coords1D",
     "CornerFinderSpec",
     "CurrentBC",
+    "CurrentIntegralAxisAlignedSpec",
     "CustomAnisotropicMedium",
     "CustomChargePerturbation",
+    "CustomCurrentIntegral2DSpec",
     "CustomCurrentSource",
     "CustomDebye",
     "CustomDrude",
@@ -463,6 +476,7 @@ def set_logging_level(level: str) -> None:
     "CustomPoleResidue",
     "CustomSellmeier",
     "CustomSourceTime",
+    "CustomVoltageIntegral2DSpec",
     "Cylinder",
     "DCCurrentSource",
     "DCVoltageSource",
@@ -668,6 +682,7 @@ def set_logging_level(level: str) -> None:
     "TemperatureBC",
     "TemperatureData",
     "TemperatureMonitor",
+    "TerminalSpec",
     "TetrahedralGridDataset",
     "Transformed",
     "TriangleMesh",
@@ -682,6 +697,7 @@ def set_logging_level(level: str) -> None:
     "Updater",
     "VisualizationSpec",
     "VoltageBC",
+    "VoltageIntegralAxisAlignedSpec",
     "VoltageSourceType",
     "VolumetricAveraging",
     "YeeGrid",

tidy3d/__init__.py

@@ -18,6 +18,7 @@
 from tidy3d.components.base_sim.data.monitor_data import AbstractMonitorData
 from tidy3d.components.grid.grid import Coords, Grid
 from tidy3d.components.medium import Medium, MediumType
+from tidy3d.components.microwave.terminal_spec import TerminalSpec
 from tidy3d.components.monitor import (
     AuxFieldTimeMonitor,
     DiffractionMonitor,
@@ -1579,11 +1580,18 @@ class ModeData(ModeSolverDataset, ElectromagneticFieldData):
 
     eps_spec: list[EpsSpecType] = pd.Field(
         None,
-        title="Permettivity Specification",
+        title="Permittivity Specification",
         description="Characterization of the permittivity profile on the plane where modes are "
         "computed. Possible values are 'diagonal', 'tensorial_real', 'tensorial_complex'.",
     )
 
+    terminal_spec: Optional[TerminalSpec] = pd.Field(
+        None,
+        title="Terminal Specification",
+        description="Specifies the method for defining microwave terminals, which is required "
+        "for defining the characteristic impedance of a transmission line mode.",
+    )
+
     @pd.validator("eps_spec", always=True)
     @skip_if_fields_missing(["monitor"])
     def eps_spec_match_mode_spec(cls, val, values):
@@ -2003,6 +2011,28 @@ def pol_fraction_waveguide(self) -> xr.Dataset:
 
         return xr.Dataset(data_vars={"te": te_frac, "tm": tm_frac})
 
+    @cached_property
+    def line_impedance(self) -> FreqModeDataArray:
+        from tidy3d.components.microwave.path_integral_factory import (
+            make_current_integral,
+            make_voltage_integral,
+        )
+        from tidy3d.plugins.microwave.impedance_calculator import ImpedanceCalculator
+
+        if self.terminal_spec is None:
+            return None
+        v_spec = self.terminal_spec.voltage_spec
+        i_spec = self.terminal_spec.current_spec
+        if v_spec is None and i_spec is None:
+            return None
+        v_integral = make_voltage_integral(v_spec)
+        i_integral = make_current_integral(i_spec)
+
+        impedance_calc = ImpedanceCalculator(
+            voltage_integral=v_integral, current_integral=i_integral
+        )
+        return impedance_calc.compute_impedance(self)
+
     @property
     def modes_info(self) -> xr.Dataset:
         """Dataset collecting various properties of the stored modes."""
@@ -2032,6 +2062,10 @@ def modes_info(self) -> xr.Dataset:
             info["wg TE fraction"] = self.pol_fraction_waveguide["te"]
             info["wg TM fraction"] = self.pol_fraction_waveguide["tm"]
 
+        if self.line_impedance is not None:
+            info["Re(Z0)"] = self.line_impedance.real
+            info["Im(Z0)"] = self.line_impedance.imag
+
         return xr.Dataset(data_vars=info)
 
     def to_dataframe(self) -> DataFrame:

tidy3d/components/data/monitor_data.py

@@ -2,12 +2,21 @@
 
 from __future__ import annotations
 
+from collections.abc import Iterable
 from enum import Enum
 from math import isclose
 from typing import Any, Optional, Union
 
 import numpy as np
-import pydantic
+import pydantic.v1 as pydantic
+from shapely.geometry import (
+    GeometryCollection,
+    MultiLineString,
+    MultiPoint,
+    MultiPolygon,
+    Polygon,
+)
+from shapely.geometry.base import BaseGeometry
 
 from tidy3d.components.base import Tidy3dBaseModel
 from tidy3d.components.geometry.base import Box
@@ -38,6 +47,46 @@
 ]
 
 
+def flatten_shapely_geometries(
+    geoms: Union[Shapely, Iterable[Shapely]], keep_types: tuple[type, ...] = (Polygon,)
+) -> list[Shapely]:
+    """
+    Flatten nested geometries into a flat list, while only keeping the specified types.
+
+    Recursively extracts and returns non-empty geometries of the given types from input geometries,
+    expanding any GeometryCollections or Multi* types.
+
+    Parameters
+    ----------
+    geoms : Union[Shapely, Iterable[Shapely]]
+        Input geometries to flatten.
+
+    keep_types : tuple[type, ...]
+        Geometry types to keep (e.g., (Polygon, LineString)). Default is
+        (Polygon).
+
+    Returns
+    -------
+    list[Shapely]
+        Flat list of non-empty geometries matching the specified types.
+    """
+    # Handle single Shapely object by wrapping it in a list
+    if isinstance(geoms, Shapely):
+        geoms = [geoms]
+
+    flat = []
+    for geom in geoms:
+        if geom.is_empty:
+            continue
+        if isinstance(geom, keep_types):
+            flat.append(geom)
+        elif isinstance(geom, (MultiPolygon, MultiLineString, MultiPoint, GeometryCollection)):
+            flat.extend(flatten_shapely_geometries(geom.geoms, keep_types))
+        elif isinstance(geom, BaseGeometry) and hasattr(geom, "geoms"):
+            flat.extend(flatten_shapely_geometries(geom.geoms, keep_types))
+    return flat
+
+
 def merging_geometries_on_plane(
     geometries: list[GeometryType],
     plane: Box,
@@ -373,6 +422,15 @@ class SnappingSpec(Tidy3dBaseModel):
         description="Describes how snapping positions will be chosen.",
     )
 
+    margin: Optional[
+        tuple[pydantic.NonNegativeInt, pydantic.NonNegativeInt, pydantic.NonNegativeInt]
+    ] = pydantic.Field(
+        (0, 0, 0),
+        title="Margin",
+        description="Number of additional grid points to consider when expanding or contracting "
+        "during snapping. Only applies when ``SnapBehavior`` is ``Expand`` or ``Contract``.",
+    )
+
 
 def get_closest_value(test: float, coords: np.ArrayLike, upper_bound_idx: int) -> float:
     """Helper to choose the closest value in an array to a given test value,
@@ -404,6 +462,8 @@ def get_lower_bound(
         using the index of the upper bound. If the test value is close to the upper
         bound, it assumes they are equal, and in that case the upper bound is returned.
         """
+        upper_bound_idx = min(upper_bound_idx, len(coords))
+        upper_bound_idx = max(upper_bound_idx, 0)
         if upper_bound_idx == len(coords):
             return coords[upper_bound_idx - 1]
         if upper_bound_idx == 0 or isclose(coords[upper_bound_idx], test, rel_tol=rel_tol):
@@ -417,14 +477,20 @@ def get_upper_bound(
         using the index of the upper bound. If the test value is close to the lower
         bound, it assumes they are equal, and in that case the lower bound is returned.
         """
+        upper_bound_idx = min(upper_bound_idx, len(coords))
+        upper_bound_idx = max(upper_bound_idx, 0)
         if upper_bound_idx == len(coords):
             return coords[upper_bound_idx - 1]
         if upper_bound_idx > 0 and isclose(coords[upper_bound_idx - 1], test, rel_tol=rel_tol):
             return coords[upper_bound_idx - 1]
         return coords[upper_bound_idx]
 
     def find_snapping_locations(
-        interval_min: float, interval_max: float, coords: np.ndarray, snap_type: SnapBehavior
+        interval_min: float,
+        interval_max: float,
+        coords: np.ndarray,
+        snap_type: SnapBehavior,
+        snap_margin: pydantic.NonNegativeInt,
     ) -> tuple[float, float]:
         """Helper that snaps a supplied interval [interval_min, interval_max] to a
         sorted array representing coordinate values.
@@ -436,9 +502,15 @@ def find_snapping_locations(
             min_snap = get_closest_value(interval_min, coords, min_upper_bound_idx)
             max_snap = get_closest_value(interval_max, coords, max_upper_bound_idx)
         elif snap_type == SnapBehavior.Expand:
+            min_upper_bound_idx -= snap_margin
+            max_upper_bound_idx += snap_margin
             min_snap = get_lower_bound(interval_min, coords, min_upper_bound_idx, rel_tol=rtol)
             max_snap = get_upper_bound(interval_max, coords, max_upper_bound_idx, rel_tol=rtol)
         else:  # SnapType.Contract
+            min_upper_bound_idx += snap_margin
+            max_upper_bound_idx -= snap_margin
+            if max_upper_bound_idx < min_upper_bound_idx:
+                raise SetupError("The supplied 'snap_buffer' is too large for this contraction.")
             min_snap = get_upper_bound(interval_min, coords, min_upper_bound_idx, rel_tol=rtol)
             max_snap = get_lower_bound(interval_max, coords, max_upper_bound_idx, rel_tol=rtol)
         return (min_snap, max_snap)
@@ -450,6 +522,7 @@ def find_snapping_locations(
     for axis in range(3):
         snap_location = snap_spec.location[axis]
         snap_type = snap_spec.behavior[axis]
+        snap_margin = snap_spec.margin[axis]
         if snap_type == SnapBehavior.Off:
             continue
         if snap_location == SnapLocation.Boundary:
@@ -460,7 +533,9 @@ def find_snapping_locations(
         box_min = min_b[axis]
         box_max = max_b[axis]
 
-        (new_min, new_max) = find_snapping_locations(box_min, box_max, snap_coords, snap_type)
+        (new_min, new_max) = find_snapping_locations(
+            box_min, box_max, snap_coords, snap_type, snap_margin
+        )
         min_b[axis] = new_min
         max_b[axis] = new_max
     return Box.from_bounds(min_b, max_b)