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

reorg of ModeData, the impedance is now calculated by the ModeSolver class avoiding the need for a duplicate terminal spec

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

@@ -1579,11 +1579,19 @@ 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'.",
     )
 
+    Z0: Optional[FreqModeDataArray] = pd.Field(
+        None,
+        title="Characteristic Impedance",
+        description="Optional quantity calculated for transmission lines."
+        "The characteristic impedance is only calculated when a :class:``TerminalSpec`` "
+        "is provided to the :class:``ModeSpec`` associated with this data.",
+    )
+
     @pd.validator("eps_spec", always=True)
     @skip_if_fields_missing(["monitor"])
     def eps_spec_match_mode_spec(cls, val, values):
@@ -2032,6 +2040,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.Z0 is not None:
+            info["Re(Z0)"] = self.Z0.real
+            info["Im(Z0)"] = self.Z0.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)