reused path integrals in smatrix plugin

Author: dmarek-flex

tidy3d/plugins/microwave/path_integrals.py

@@ -6,28 +6,44 @@
 import numpy as np
 from typing import Tuple
 from ...components.data.dataset import FieldDataset, FieldTimeDataset, ModeSolverDataset
-from ...components.data.data_array import (
-    ScalarFieldDataArray,
-    ScalarModeFieldDataArray,
-    ScalarFieldTimeDataArray,
-)
+from ...components.data.data_array import AbstractSpatialDataArray
 from ...components.base import cached_property
 from ...components.types import Axis, Direction
 from ...components.geometry.base import Box
 from ...components.validators import assert_line, assert_plane
 from ...exceptions import DataError
 
 
-class AxisAlignedPathIntegral(Box):
+class AbstractAxesRH:
+    """Represents an axis-aligned right-handed coordinate system with one axis preferred."""
+
+    @cached_property
+    def main_axis(self) -> Axis:
+        """Subclasses should implement this method."""
+        raise NotImplementedError()
+
+    @cached_property
+    def remaining_axes(self) -> Tuple[Axis, Axis]:
+        """Axes in plane, ordered to maintain a right-handed coordinate system"""
+        axes = [0, 1, 2]
+        axes.pop(self.main_axis)
+        if self.main_axis == 1:
+            return (axes[1], axes[0])
+        else:
+            return (axes[0], axes[1])
+
+
+class AxisAlignedPathIntegral(AbstractAxesRH, Box):
     """Class for defining the simplest type of path integral which is aligned with Cartesian axes."""
 
     _line_validator = assert_line()
 
     extrapolate_to_endpoints: bool = pd.Field(
         False,
         title="Extrapolate to endpoints",
-        description="If the endpoints of the path integral terminate at or near an interface, the field is likely discontinuous. "
-        "This option ignores fields outside the bounds of the integral.",
+        description="If the endpoints of the path integral terminate at or near a material interface, "
+        "the field is likely discontinuous. This option ignores fields outside and on the bounds "
+        "of the integral. Should be turned on when computing voltage between two conductors. ",
     )
 
     snap_path_to_grid: bool = pd.Field(
@@ -37,15 +53,16 @@ class AxisAlignedPathIntegral(Box):
         "a field. If enabled, the integration path will be snapped to the grid.",
     )
 
-    def compute_integral(self, scalar_field):
+    def compute_integral(self, scalar_field: AbstractSpatialDataArray):
+        """Computes the defined integral given the input `scalar_field`."""
         if not scalar_field.does_cover(self.bounds):
             raise DataError("scalar field does not cover the integration domain")
-        coord = "xyz"[self.integration_axis]
+        coord = "xyz"[self.main_axis]
 
         scalar_field = self._get_field_along_path(scalar_field)
         # Get the boundaries
-        min_bound = self.bounds[0][self.integration_axis]
-        max_bound = self.bounds[1][self.integration_axis]
+        min_bound = self.bounds[0][self.main_axis]
+        max_bound = self.bounds[1][self.main_axis]
 
         if self.extrapolate_to_endpoints:
             # Remove field outside the boundaries
@@ -62,16 +79,20 @@ def compute_integral(self, scalar_field):
         coords_interp = np.concatenate((coords_interp, coordinates))
         coords_interp = np.concatenate((coords_interp, [max_bound]))
         coords_interp = {coord: coords_interp}
-        # Use extrapolation for the 2 additional endpoints
+
+        # Use extrapolation for the 2 additional endpoints, unless there is only a single sample point
+        method = "linear"
+        if len(coordinates) == 1 and self.extrapolate_to_endpoints:
+            method = "nearest"
         scalar_field = scalar_field.interp(
-            coords_interp, method="linear", kwargs={"fill_value": "extrapolate"}
+            coords_interp, method=method, kwargs={"fill_value": "extrapolate"}
         )
         return scalar_field.integrate(coord=coord)
 
     def _get_field_along_path(
-        self,
-        scalar_field: ScalarFieldDataArray | ScalarModeFieldDataArray | ScalarFieldTimeDataArray,
-    ) -> ScalarFieldDataArray | ScalarModeFieldDataArray | ScalarFieldTimeDataArray:
+        self, scalar_field: AbstractSpatialDataArray
+    ) -> AbstractSpatialDataArray:
+        """Returns a selection of the input `scalar_field` ready for integration."""
         axis1 = self.remaining_axes[0]
         axis2 = self.remaining_axes[1]
         coord1 = "xyz"[axis1]
@@ -109,18 +130,11 @@ def _get_field_along_path(
         return scalar_field
 
     @cached_property
-    def integration_axis(self) -> Axis:
-        """Axis for performing integration"""
+    def main_axis(self) -> Axis:
+        """Axis for performing integration."""
         val = next((index for index, value in enumerate(self.size) if value != 0), None)
         return val
 
-    @cached_property
-    def remaining_axes(self) -> Tuple[Axis, Axis]:
-        """Axes perpendicular to the voltage axis"""
-        axes = [0, 1, 2]
-        axes.pop(self.integration_axis)
-        return (axes[0], axes[1])
-
 
 class VoltageIntegralAA(AxisAlignedPathIntegral):
     """Class for computing the voltage between two points defined by an axis-aligned line."""
@@ -133,7 +147,7 @@ class VoltageIntegralAA(AxisAlignedPathIntegral):
 
     def compute_voltage(self, em_field: FieldDataset | ModeSolverDataset | FieldTimeDataset):
         """Compute voltage along path defined by a line."""
-        e_component = "xyz"[self.integration_axis]
+        e_component = "xyz"[self.main_axis]
         field_name = f"E{e_component}"
         # Validate that the field is present
         if field_name not in em_field:
@@ -148,22 +162,27 @@ def compute_voltage(self, em_field: FieldDataset | ModeSolverDataset | FieldTime
         return voltage
 
 
-class CurrentIntegralAA(Box):
+class CurrentIntegralAA(AbstractAxesRH, Box):
     """Class for computing conduction current via Ampere's Circuital Law on an axis-aligned loop."""
 
     _plane_validator = assert_plane()
 
     sign: Direction = pd.Field(
         "+",
-        title="Direction of contour integral.",
-        description="Positive indicates V=Vb-Va where position b has a larger coordinate along the axis of integration.",
+        title="Direction of contour integral",
+        description="Positive indicates current flowing in the positive normal axis direction.",
+    )
+
+    extrapolate_to_endpoints: bool = pd.Field(
+        False,
+        title="Extrapolate to endpoints",
+        description="This parameter is passed to `AxisAlignedPathIntegral` objects when computing the contour integral.",
     )
 
     snap_contour_to_grid: bool = pd.Field(
         False,
-        title="",
-        description="It might be desireable to integrate exactly along the Yee grid associated with "
-        "the fields. If enabled, the integration path will be snapped to the grid.",
+        title="Snap contour to grid",
+        description="This parameter is passed to `AxisAlignedPathIntegral` objects when computing the contour integral.",
     )
 
     def compute_current(self, em_field: FieldDataset | ModeSolverDataset | FieldTimeDataset):
@@ -200,21 +219,11 @@ def compute_current(self, em_field: FieldDataset | ModeSolverDataset | FieldTime
         return current
 
     @cached_property
-    def normal_axis(self) -> Axis:
+    def main_axis(self) -> Axis:
         """Axis normal to loop"""
         val = next((index for index, value in enumerate(self.size) if value == 0), None)
         return val
 
-    @cached_property
-    def remaining_axes(self) -> Tuple[Axis, Axis]:
-        """Axes in integration plane, ordered to maintain a right-handed coordinate system"""
-        axes = [0, 1, 2]
-        axes.pop(self.normal_axis)
-        if self.normal_axis == 1:
-            return (axes[1], axes[0])
-        else:
-            return (axes[0], axes[1])
-
     def _to_path_integrals(self, h_horizontal, h_vertical) -> Tuple[AxisAlignedPathIntegral, ...]:
         ax1 = self.remaining_axes[0]
         ax2 = self.remaining_axes[1]
@@ -255,14 +264,14 @@ def _to_path_integrals(self, h_horizontal, h_vertical) -> Tuple[AxisAlignedPathI
         bottom = AxisAlignedPathIntegral(
             center=path_center,
             size=path_size,
-            extrapolate_to_endpoints=False,
+            extrapolate_to_endpoints=self.extrapolate_to_endpoints,
             snap_path_to_grid=self.snap_contour_to_grid,
         )
         path_center[ax2] = top_bound
         top = AxisAlignedPathIntegral(
             center=path_center,
             size=path_size,
-            extrapolate_to_endpoints=False,
+            extrapolate_to_endpoints=self.extrapolate_to_endpoints,
             snap_path_to_grid=self.snap_contour_to_grid,
         )
 
@@ -276,14 +285,14 @@ def _to_path_integrals(self, h_horizontal, h_vertical) -> Tuple[AxisAlignedPathI
         left = AxisAlignedPathIntegral(
             center=path_center,
             size=path_size,
-            extrapolate_to_endpoints=False,
+            extrapolate_to_endpoints=self.extrapolate_to_endpoints,
             snap_path_to_grid=self.snap_contour_to_grid,
         )
         path_center[ax1] = right_bound
         right = AxisAlignedPathIntegral(
             center=path_center,
             size=path_size,
-            extrapolate_to_endpoints=False,
+            extrapolate_to_endpoints=self.extrapolate_to_endpoints,
             snap_path_to_grid=self.snap_contour_to_grid,
         )
 

tidy3d/plugins/smatrix/component_modelers/terminal.py

@@ -23,6 +23,8 @@
 from .base import AbstractComponentModeler, FWIDTH_FRAC
 from ..ports.lumped import LumpedPortDataArray, LumpedPort
 
+from ...microwave import VoltageIntegralAA, CurrentIntegralAA
+
 
 class TerminalComponentModeler(AbstractComponentModeler):
     """Tool for modeling two-terminal multiport devices and computing port parameters
@@ -206,28 +208,17 @@ def port_voltage(port: LumpedPort, sim_data: SimulationData) -> xr.DataArray:
             """Helper to compute voltage across the port."""
             e_component = "xyz"[port.voltage_axis]
             field_data = sim_data[f"{port.name}_E{e_component}"]
-            e_field = field_data.field_components[f"E{e_component}"]
-            # Get the boundaries of the port along the voltage axis
-            min_port_bound = port.bounds[0][port.voltage_axis]
-            max_port_bound = port.bounds[1][port.voltage_axis]
-            # Remove E field outside the port region
-            e_field = e_field.sel({e_component: slice(min_port_bound, max_port_bound)})
-            # Ignore values on the port boundary, which are likely considered within the conductor
-            e_field = e_field.drop_sel(
-                {e_component: (min_port_bound, max_port_bound)}, errors="ignore"
-            )
-            e_coords = [e_field.x, e_field.y, e_field.z]
-            # Integration is along the original coordinates plus two additional
-            # endpoints corresponding to the precise bounds of the port
-            e_coords_interp = np.array([min_port_bound])
-            e_coords_interp = np.concatenate((e_coords_interp, e_coords[port.voltage_axis].values))
-            e_coords_interp = np.concatenate((e_coords_interp, [max_port_bound]))
-            e_coords_interp = {e_component: e_coords_interp}
-            # Use extrapolation for the 2 additional endpoints
-            e_field = e_field.interp(
-                **e_coords_interp, method="linear", kwargs={"fill_value": "extrapolate"}
+
+            size = list(port.size)
+            size[port.current_axis] = 0
+            voltage_integral = VoltageIntegralAA(
+                center=port.center,
+                size=size,
+                extrapolate_to_endpoints=True,
+                snap_path_to_grid=True,
+                sign="+",
             )
-            voltage = -e_field.integrate(coord=e_component).squeeze(drop=True)
+            voltage = voltage_integral.compute_voltage(field_data.field_components)
             # Return data array of voltage with coordinates of frequency
             return voltage
 
@@ -246,83 +237,59 @@ def port_current(port: LumpedPort, sim_data: SimulationData) -> xr.DataArray:
 
             # Get h field tangent to resistive sheet
             h_component = "xyz"[port.current_axis]
-            orth_component = "xyz"[port.injection_axis]
-            field_data = sim_data[f"{port.name}_H{h_component}"]
-            h_field = field_data.field_components[f"H{h_component}"]
-            h_coords = [h_field.x, h_field.y, h_field.z]
+            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
+            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.
-            h_cap_field = field_data.field_components[f"H{orth_component}"]
-            h_cap_coords = [h_cap_field.x, h_cap_field.y, h_cap_field.z]
+            h_cap_field = field_data[f"H{inject_component}"]
+            # Coordinates of h_cap field as numpy arrays
+            h_cap_coords_along_current = h_cap_field.coords[h_component].values
+            h_cap_coords_along_injection = h_cap_field.coords[inject_component].values
 
             # Use the coordinates of h_cap since it lies on the same grid that the
             # lumped resistor is snapped to
             orth_index = np.argmin(
-                np.abs(h_cap_coords[port.injection_axis].values - port.center[port.injection_axis])
+                np.abs(h_cap_coords_along_injection - port.center[port.injection_axis])
             )
+            inject_center = h_cap_coords_along_injection[orth_index]
             # Some sanity checks, tangent H field coordinates should be directly above
             # and below the coordinates of the resistive sheet
             assert orth_index > 0
-            assert (
-                h_cap_coords[port.injection_axis].values[orth_index]
-                < h_coords[port.injection_axis].values[orth_index]
-            )
-            assert (
-                h_coords[port.injection_axis].values[orth_index - 1]
-                < h_cap_coords[port.injection_axis].values[orth_index]
-            )
-
-            # Extract field just below and just above sheet
-            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
+            assert inject_center < h_coords_along_injection[orth_index]
+            assert h_coords_along_injection[orth_index - 1] < inject_center
+            # Distance between the h1_field and h2_field, a single cell size
+            dcap = h_coords_along_injection[orth_index] - h_coords_along_injection[orth_index - 1]
 
+            # Next find the size in the current_axis direction
             # Find exact bounds of port taking into consideration the Yee grid
-            np_coords = h_coords[port.current_axis].values
+            # 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]
-            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, 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
-            dcap = (
-                h_coords[port.injection_axis].values[orth_index]
-                - h_coords[port.injection_axis].values[orth_index - 1]
-            )
-
-            h_min_bound = hcap_minus.coords[h_component].values
-            h_max_bound = hcap_plus.coords[h_component].values
-            h_coords_interp = {
-                h_component: np.linspace(
-                    h_min_bound,
-                    h_max_bound,
-                    len(h_coords[port.current_axis] + 2),
-                )
-            }
-            # Integration that corresponds to the tangent H field
-            h_field = h_field.interp(**h_coords_interp)
-            current = h_field.integrate(coord=h_component).squeeze(drop=True)
-
-            # Integration that corresponds with the contribution to current from cap contours
-            hcap_current = (
-                ((hcap_plus - hcap_minus) * dcap).squeeze(drop=True).reset_coords(drop=True)
+            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]
+
+            # Setup axis aligned contour integral, which is defined by a plane
+            # The path integral is snapped to the grid, so center and size will
+            # be slightly modified when compared to the original port.
+            center = list(port.center)
+            center[port.injection_axis] = inject_center
+            center[port.current_axis] = (h_max_bound + h_min_bound) / 2
+            size = [0, 0, 0]
+            size[port.current_axis] = h_max_bound - h_min_bound
+            size[port.injection_axis] = dcap
+
+            # H field is continuous at integral bounds, so extrapolation is turned off
+            I_integral = CurrentIntegralAA(
+                center=center,
+                size=size,
+                sign="+",
+                extrapolate_to_endpoints=False,
+                snap_contour_to_grid=True,
             )
-            # Add the contribution from the hcap integral
-            current = current + hcap_current
-            # Make sure we compute current flowing from plus to minus voltage
-            if port.current_axis != (port.voltage_axis + 1) % 3:
-                current *= -1
-            # Return data array of current with coordinates of frequency
-            return current
+            return I_integral.compute_current(field_data)
 
         def port_ab(port: LumpedPort, sim_data: SimulationData):
             """Helper to compute the port incident and reflected power waves."""