add unit tests for path integrals

tests/test_plugins/test_microwave.py

@@ -0,0 +1,244 @@
+import pytest
+import numpy as np
+
+import tidy3d as td
+from tidy3d import FieldData
+from tidy3d.constants import ETA_0
+from tidy3d.plugins.microwave import VoltageIntegralAA, CurrentIntegralAA, ImpedanceCalculator
+import pydantic.v1 as pydantic
+from tidy3d.exceptions import DataError
+from ..utils import run_emulated
+
+
+# Using similar code as "test_data/test_data_arrays.py"
+MON_SIZE = (2, 1, 0)
+FIELDS = ("Ex", "Ey", "Hx", "Hy")
+FSTART = 0.5e9
+FSTOP = 1.5e9
+F0 = (FSTART + FSTOP) / 2
+FWIDTH = FSTOP - FSTART
+FS = np.linspace(FSTART, FSTOP, 5)
+FIELD_MONITOR = td.FieldMonitor(
+    size=MON_SIZE, fields=FIELDS, name="strip_field", freqs=FS, colocate=False
+)
+STRIP_WIDTH = 1.5
+STRIP_HEIGHT = 0.5
+
+SIM_Z = td.Simulation(
+    size=(2, 1, 1),
+    grid_spec=td.GridSpec.uniform(dl=0.04),
+    monitors=[
+        FIELD_MONITOR,
+        td.FieldMonitor(center=(0, 0, 0), size=(1, 1, 1), freqs=FS, name="field"),
+        td.FieldMonitor(
+            center=(0, 0, 0), size=(1, 1, 1), freqs=FS, fields=["Ex", "Hx"], name="ExHx"
+        ),
+        td.ModeMonitor(
+            center=(0, 0, 0), size=(1, 1, 0), freqs=FS, mode_spec=td.ModeSpec(), name="mode"
+        ),
+    ],
+    sources=[
+        td.PointDipole(
+            center=(0, 0, 0),
+            polarization="Ex",
+            source_time=td.GaussianPulse(freq0=F0, fwidth=FWIDTH),
+        )
+    ],
+    run_time=2e-12,
+)
+
+""" Generate the data arrays for testing path integral computations """
+
+
+def get_xyz(
+    monitor: td.components.monitor.MonitorType, grid_key: str
+) -> tuple[list[float], list[float], list[float]]:
+    grid = SIM_Z.discretize_monitor(monitor)
+    x, y, z = grid[grid_key].to_list
+    return x, y, z
+
+
+def make_stripline_scalar_field_data_array(grid_key: str):
+    """Populate FIELD_MONITOR with a idealized stripline mode, where fringing fields are assumed 0."""
+    XS, YS, ZS = get_xyz(FIELD_MONITOR, grid_key)
+    XGRID, YGRID = np.meshgrid(XS, YS, indexing="ij")
+    XGRID = XGRID.reshape((len(XS), len(YS), 1, 1))
+    YGRID = YGRID.reshape((len(XS), len(YS), 1, 1))
+    values = np.zeros((len(XS), len(YS), len(ZS), len(FS)))
+    ones = np.ones((len(XS), len(YS), len(ZS), len(FS)))
+    XGRID = np.broadcast_to(XGRID, values.shape)
+    YGRID = np.broadcast_to(YGRID, values.shape)
+
+    # Numpy masks for quickly determining location
+    above_in_strip = np.logical_and(YGRID >= 0, YGRID <= STRIP_HEIGHT / 2)
+    below_in_strip = np.logical_and(YGRID < 0, YGRID >= -STRIP_HEIGHT / 2)
+    within_strip_width = np.logical_and(XGRID >= -STRIP_WIDTH / 2, XGRID < STRIP_WIDTH / 2)
+    above_and_within = np.logical_and(above_in_strip, within_strip_width)
+    below_and_within = np.logical_and(below_in_strip, within_strip_width)
+    # E field is perpendicular to strip surface and magnetic field is parallel
+    if grid_key == "Ey":
+        values = np.where(above_and_within, ones, values)
+        values = np.where(below_and_within, -ones, values)
+    elif grid_key == "Hx":
+        values = np.where(above_and_within, -ones / ETA_0, values)
+        values = np.where(below_and_within, ones / ETA_0, values)
+
+    return td.ScalarFieldDataArray(values, coords=dict(x=XS, y=YS, z=ZS, f=FS))
+
+
+def make_field_data():
+    return FieldData(
+        monitor=FIELD_MONITOR,
+        Ex=make_stripline_scalar_field_data_array("Ex"),
+        Ey=make_stripline_scalar_field_data_array("Ey"),
+        Hx=make_stripline_scalar_field_data_array("Hx"),
+        Hy=make_stripline_scalar_field_data_array("Hy"),
+        symmetry=SIM_Z.symmetry,
+        symmetry_center=SIM_Z.center,
+        grid_expanded=SIM_Z.discretize_monitor(FIELD_MONITOR),
+    )
+
+
+@pytest.mark.parametrize("axis", [0, 1, 2])
+def test_voltage_integral_axes(axis):
+    length = 0.5
+    size = [0, 0, 0]
+    size[axis] = length
+    center = [0, 0, 0]
+    voltage_integral = VoltageIntegralAA(
+        center=center,
+        size=size,
+    )
+    sim = SIM_Z
+    sim_data = run_emulated(sim)
+    _ = voltage_integral.compute_voltage(sim_data["field"].field_components)
+
+
+@pytest.mark.parametrize("axis", [0, 1, 2])
+def test_current_integral_axes(axis):
+    length = 0.5
+    size = [length, length, length]
+    size[axis] = 0.0
+    center = [0, 0, 0]
+    current_integral = CurrentIntegralAA(
+        center=center,
+        size=size,
+    )
+    sim = SIM_Z
+    sim_data = run_emulated(sim)
+    _ = current_integral.compute_current(sim_data["field"].field_components)
+
+
+def test_voltage_integral_toggles():
+    length = 0.5
+    size = [0, 0, 0]
+    size[0] = length
+    center = [0, 0, 0]
+    voltage_integral = VoltageIntegralAA(
+        center=center,
+        size=size,
+        extrapolate_to_endpoints=True,
+        snap_path_to_grid=True,
+        sign="-",
+    )
+    sim = SIM_Z
+    sim_data = run_emulated(sim)
+    _ = voltage_integral.compute_voltage(sim_data["field"].field_components)
+
+
+def test_current_integral_toggles():
+    length = 0.5
+    size = [length, length, length]
+    size[0] = 0.0
+    center = [0, 0, 0]
+    current_integral = CurrentIntegralAA(
+        center=center,
+        size=size,
+        extrapolate_to_endpoints=True,
+        snap_contour_to_grid=True,
+        sign="-",
+    )
+    sim = SIM_Z
+    sim_data = run_emulated(sim)
+    _ = current_integral.compute_current(sim_data["field"].field_components)
+
+
+def test_voltage_missing_fields():
+    length = 0.5
+    size = [0, 0, 0]
+    size[1] = length
+    center = [0, 0, 0]
+    voltage_integral = VoltageIntegralAA(
+        center=center,
+        size=size,
+    )
+    sim = SIM_Z
+    sim_data = run_emulated(sim)
+    with pytest.raises(DataError):
+        _ = voltage_integral.compute_voltage(sim_data["ExHx"].field_components)
+
+
+def test_current_missing_fields():
+    length = 0.5
+    size = [length, length, length]
+    size[0] = 0.0
+    center = [0, 0, 0]
+    current_integral = CurrentIntegralAA(
+        center=center,
+        size=size,
+    )
+    sim = SIM_Z
+    sim_data = run_emulated(sim)
+    with pytest.raises(DataError):
+        _ = current_integral.compute_current(sim_data["ExHx"].field_components)
+
+
+def test_tiny_voltage_path():
+    length = 0.02
+    size = [0, 0, 0]
+    size[1] = length
+    center = [0, 0, 0]
+    voltage_integral = VoltageIntegralAA(center=center, size=size, extrapolate_to_endpoints=True)
+    sim = SIM_Z
+    sim_data = run_emulated(sim)
+    _ = voltage_integral.compute_voltage(sim_data["field"].field_components)
+
+
+def test_impedance_calculator():
+    with pytest.raises(pydantic.ValidationError):
+        _ = ImpedanceCalculator(voltage_integral=None, current_integral=None)
+
+
+def test_impedance_accuracy():
+    field_data = make_field_data()
+    # Setup path integrals
+    size = [0, STRIP_HEIGHT / 2, 0]
+    center = [0, -STRIP_HEIGHT / 4, 0]
+    voltage_integral = VoltageIntegralAA(center=center, size=size, extrapolate_to_endpoints=True)
+
+    size = [STRIP_WIDTH * 1.25, STRIP_HEIGHT / 2, 0]
+    center = [0, 0, 0]
+    current_integral = CurrentIntegralAA(center=center, size=size)
+
+    def impedance_of_stripline(width, height):
+        # Assuming no fringing fields, is the same as a parallel plate
+        # with half the height and carrying twice the current
+        Z0_parallel_plate = 0.5 * height / width * td.ETA_0
+        return Z0_parallel_plate / 2
+
+    analytic_impedance = impedance_of_stripline(STRIP_WIDTH, STRIP_HEIGHT)
+
+    # Compute impedance using the tool
+    Z_calc = ImpedanceCalculator(
+        voltage_integral=voltage_integral, current_integral=current_integral
+    )
+    Z1 = Z_calc.compute_impedance(field_data)
+    Z_calc = ImpedanceCalculator(voltage_integral=voltage_integral, current_integral=None)
+    Z2 = Z_calc.compute_impedance(field_data)
+    Z_calc = ImpedanceCalculator(voltage_integral=None, current_integral=current_integral)
+    Z3 = Z_calc.compute_impedance(field_data)
+
+    # Computation that uses the flux is less accurate, due to staircasing the field
+    assert np.all(np.isclose(Z1, analytic_impedance, rtol=0.02))
+    assert np.all(np.isclose(Z2, analytic_impedance, atol=3.5))
+    assert np.all(np.isclose(Z3, analytic_impedance, atol=3.5))

tidy3d/plugins/microwave/impedance_calculator.py

@@ -9,7 +9,7 @@
 from ...components.data.monitor_data import FieldData, FieldTimeData, ModeSolverData
 
 from ...components.base import Tidy3dBaseModel
-from ...components.validators import ValidationError
+from ...exceptions import ValidationError
 
 from .path_integrals import VoltageIntegralAA, CurrentIntegralAA
 
@@ -38,6 +38,8 @@ def compute_impedance(self, em_field: FieldData | ModeSolverData | FieldTimeData
 
         # If only one of the integrals has been provided then fall back to using total power (flux)
         # with Ohm's law. The input field should cover an area large enough to render the flux computation accurate.
+        # If the input field is a time signal, then it is real and flux corresponds to the instantaneous power.
+        # Otherwise the input field is in frequency domain, where flux indicates the time-averaged power 0.5*Re(V*conj(I))
         if not self.voltage_integral:
             flux = em_field.flux
             if isinstance(em_field, FieldTimeData):
@@ -53,10 +55,10 @@ def compute_impedance(self, em_field: FieldData | ModeSolverData | FieldTimeData
 
         return voltage / current
 
-    @pd.validator("voltage_integral", always=True)
+    @pd.validator("current_integral", always=True)
     def check_voltage_or_current(cls, val, values):
         """Ensure that 'voltage_integral' and/or 'current_integral' were provided."""
-        if not values.get("current_integral") and not val:
+        if not values.get("voltage_integral") and not val:
             raise ValidationError(
                 "Atleast one of 'voltage_integral' or 'current_integral' must be provided."
             )

tidy3d/plugins/smatrix/component_modelers/terminal.py

@@ -196,13 +196,17 @@ def _construct_smatrix(self, batch_data: BatchData) -> LumpedPortDataArray:
         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
+            """Helper to return indices of coordinates within min and max bounds,"""
+            """including a tolerance. xarray does not have this functionality yet. """
+            min_idx = np.searchsorted(coords, min, "left")
+            # If a coordinate is close enough, it is considered included
+            if min_idx > 0 and np.isclose(coords[min_idx - 1], min, rtol=fp_eps, atol=fp_eps):
+                min_idx -= 1
+            max_idx = np.searchsorted(coords, max, "left")
+            if max_idx < len(coords) and np.isclose(coords[max_idx], max, rtol=fp_eps, atol=fp_eps):
+                max_idx += 1
+
+            return (min_idx, max_idx - 1)
 
         def port_voltage(port: LumpedPort, sim_data: SimulationData) -> xr.DataArray:
             """Helper to compute voltage across the port."""
@@ -240,7 +244,8 @@ def port_current(port: LumpedPort, sim_data: SimulationData) -> xr.DataArray:
             inject_component = "xyz"[port.injection_axis]
             field_data = sim_data[f"{port.name}_H{h_component}"].field_components
             h_field = field_data[f"H{h_component}"]
-            # Coordinates as numpy array for h_field along injection axis
+            # Coordinates as numpy array for h_field along curren and injection axis
+            h_coords_along_current = h_field.coords[h_component].values
             h_coords_along_injection = h_field.coords[inject_component].values
             # h_cap represents the very short section (single cell) of the H contour that
             # is in the injection_axis direction. It is needed to fully enclose the sheet.
@@ -268,8 +273,11 @@ def port_current(port: LumpedPort, sim_data: SimulationData) -> xr.DataArray:
             # Select bounds carefully and allow for h_cap very close to the port bounds
             port_min = port.bounds[0][port.current_axis]
             port_max = port.bounds[1][port.current_axis]
-            h_min_bound = select_within_bounds(h_cap_coords_along_current, -np.inf, port_min)[-1]
-            h_max_bound = select_within_bounds(h_cap_coords_along_current, port_max, np.inf)[0]
+
+            (idx_min, idx_max) = select_within_bounds(h_coords_along_current, port_min, port_max)
+            # Use these indices to select the exact positions of the h_cap field
+            h_min_bound = h_cap_coords_along_current[idx_min - 1]
+            h_max_bound = h_cap_coords_along_current[idx_max]
 
             # Setup axis aligned contour integral, which is defined by a plane
             # The path integral is snapped to the grid, so center and size will