adding microwave plugin

voltage computation completed

finished up current and voltage computations with extra options for snapping to grid

added the impedance calculator tool, and slight reorganization

Author: dmarek-flex

tidy3d/components/validators.py

@@ -46,6 +46,19 @@
 MIN_FREQUENCY = 1e5
 
 
+def assert_line():
+    """makes sure a field's ``size`` attribute has exactly 2 zeros"""
+
+    @pydantic.validator("size", allow_reuse=True, always=True)
+    def is_line(cls, val):
+        """Raise validation error if not 1 dimensional."""
+        if val.count(0.0) != 2:
+            raise ValidationError(f"'{cls.__name__}' object must be a line, given size={val}")
+        return val
+
+    return is_line
+
+
 def assert_plane():
     """makes sure a field's ``size`` attribute has exactly 1 zero"""
 

tidy3d/plugins/microwave/__init__.py

@@ -0,0 +1,10 @@
+""" Imports from microwave plugin. """
+
+from .path_integrals import VoltageIntegralAA, CurrentIntegralAA
+from .impedance_calculator import ImpedanceCalculator
+
+__all__ = [
+    "VoltageIntegralAA",
+    "CurrentIntegralAA",
+    "ImpedanceCalculator",
+]

tidy3d/plugins/microwave/impedance_calculator.py

@@ -0,0 +1,63 @@
+"""Class for computing characteristic impedance of transmission lines."""
+
+from __future__ import annotations
+
+
+import pydantic.v1 as pd
+import numpy as np
+from typing import Optional
+from ...components.data.monitor_data import FieldData, FieldTimeData, ModeSolverData
+
+from ...components.base import Tidy3dBaseModel
+from ...components.validators import ValidationError
+
+from .path_integrals import VoltageIntegralAA, CurrentIntegralAA
+
+
+class ImpedanceCalculator(Tidy3dBaseModel):
+    """Tool for computing the characteristic impedance of a transmission line."""
+
+    voltage_integral: Optional[VoltageIntegralAA] = pd.Field(
+        ...,
+        title="Voltage Integral",
+        description="Integral for computing voltage.",
+    )
+
+    current_integral: Optional[CurrentIntegralAA] = pd.Field(
+        ...,
+        title="",
+        description="Integral for computing current.",
+    )
+
+    def compute_impedance(self, em_field: FieldData | ModeSolverData | FieldTimeData):
+        # If both voltage and current integrals have been defined then impedance is computed directly
+        if self.voltage_integral:
+            voltage = self.voltage_integral.compute_voltage(em_field.field_components)
+        if self.current_integral:
+            current = self.current_integral.compute_current(em_field.field_components)
+
+        # 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 not self.voltage_integral:
+            flux = em_field.flux
+            if isinstance(em_field, FieldTimeData):
+                voltage = flux / current
+            else:
+                voltage = 2 * flux / np.conj(current)
+        if not self.current_integral:
+            flux = em_field.flux
+            if isinstance(em_field, FieldTimeData):
+                current = flux / voltage
+            else:
+                current = np.conj(2 * flux / voltage)
+
+        return voltage / current
+
+    @pd.validator("voltage_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:
+            raise ValidationError(
+                "Atleast one of 'voltage_integral' or 'current_integral' must be provided."
+            )
+        return val

tidy3d/plugins/microwave/path_integrals.py

@@ -0,0 +1,290 @@
+"""Helper classes for performing path integrals with fields on the Yee grid"""
+
+from __future__ import annotations
+
+import pydantic.v1 as pd
+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.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 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.",
+    )
+
+    snap_path_to_grid: bool = pd.Field(
+        False,
+        title="Snap path to grid",
+        description="It might be desireable to integrate exactly along the Yee grid associated with "
+        "a field. If enabled, the integration path will be snapped to the grid.",
+    )
+
+    def compute_integral(self, 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]
+
+        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]
+
+        if self.extrapolate_to_endpoints:
+            # Remove field outside the boundaries
+            scalar_field = scalar_field.sel({coord: slice(min_bound, max_bound)})
+            # Ignore values on the boundary (sel is inclusive)
+            scalar_field = scalar_field.drop_sel({coord: (min_bound, max_bound)}, errors="ignore")
+            coordinates = scalar_field.coords[coord].values
+        else:
+            coordinates = scalar_field.coords[coord].sel({coord: slice(min_bound, max_bound)})
+
+        # Integration is along the original coordinates plus ensure that
+        # endpoints corresponding to the precise bounds of the port are included
+        coords_interp = np.array([min_bound])
+        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
+        scalar_field = scalar_field.interp(
+            coords_interp, method="linear", kwargs={"fill_value": "extrapolate"}
+        )
+        return scalar_field.integrate(coord=coord)
+
+    def _get_field_along_path(
+        self,
+        scalar_field: ScalarFieldDataArray | ScalarModeFieldDataArray | ScalarFieldTimeDataArray,
+    ) -> ScalarFieldDataArray | ScalarModeFieldDataArray | ScalarFieldTimeDataArray:
+        axis1 = self.remaining_axes[0]
+        axis2 = self.remaining_axes[1]
+        coord1 = "xyz"[axis1]
+        coord2 = "xyz"[axis2]
+
+        if self.snap_path_to_grid:
+            # Coordinates that are not integrated
+            remaining_coords = {
+                coord1: self.center[axis1],
+                coord2: self.center[axis2],
+            }
+            # Select field nearest to center of integration line
+            scalar_field = scalar_field.sel(
+                remaining_coords,
+                method="nearest",
+                drop=False,
+            )
+        else:
+            # Try to interpolate unless there is only a single coordinate
+            coord1dict = {coord1: self.center[axis1]}
+            if scalar_field.sizes[coord1] == 1:
+                scalar_field = scalar_field.sel(coord1dict, method="nearest")
+            else:
+                scalar_field = scalar_field.interp(
+                    coord1dict, method="linear", kwargs={"bounds_error": True}
+                )
+            coord2dict = {coord2: self.center[axis2]}
+            if scalar_field.sizes[coord2] == 1:
+                scalar_field = scalar_field.sel(coord2dict, method="nearest")
+            else:
+                scalar_field = scalar_field.interp(
+                    coord2dict, method="linear", kwargs={"bounds_error": True}
+                )
+
+        return scalar_field
+
+    @cached_property
+    def integration_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."""
+
+    sign: Direction = pd.Field(
+        "+",
+        title="Direction of path integral.",
+        description="Positive indicates V=Vb-Va where position b has a larger coordinate along the axis of integration.",
+    )
+
+    def compute_voltage(self, em_field: FieldDataset | ModeSolverDataset | FieldTimeDataset):
+        """Compute voltage along path defined by a line."""
+        e_component = "xyz"[self.integration_axis]
+        field_name = f"E{e_component}"
+        # Validate that the field is present
+        if field_name not in em_field:
+            raise DataError(f"field_name '{field_name}' not found")
+        e_field = em_field[f"E{e_component}"]
+        integral = self.compute_integral(e_field)
+
+        voltage = -integral
+        if self.sign == "-":
+            voltage *= -1
+        # Return data array of voltage while keeping coordinates of frequency|time|mode index
+        return voltage
+
+
+class CurrentIntegralAA(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.",
+    )
+
+    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.",
+    )
+
+    def compute_current(self, em_field: FieldDataset | ModeSolverDataset | FieldTimeDataset):
+        """Compute current flowing in loop defined by the outer edge of a rectangle."""
+
+        ax1 = self.remaining_axes[0]
+        ax2 = self.remaining_axes[1]
+        h_component = "xyz"[ax1]
+        v_component = "xyz"[ax2]
+        h_field_name = f"H{h_component}"
+        v_field_name = f"H{v_component}"
+        # Validate that fields are present
+        if h_field_name not in em_field:
+            raise DataError(f"field_name '{h_field_name}' not found")
+        if v_field_name not in em_field:
+            raise DataError(f"field_name '{v_field_name}' not found")
+        h_horizontal = em_field[f"H{h_component}"]
+        h_vertical = em_field[f"H{v_component}"]
+
+        # Decompose contour into path integrals
+        (bottom, right, top, left) = self._to_path_integrals(h_horizontal, h_vertical)
+
+        current = 0
+        # Compute and add contributions from each part of the contour
+        current += bottom.compute_integral(h_horizontal)
+        current += right.compute_integral(h_vertical)
+        current -= top.compute_integral(h_horizontal)
+        current -= left.compute_integral(h_vertical)
+
+        if self.sign == "-":
+            current *= -1
+        # Return data array of current while keeping coordinates of frequency|time|mode index
+        current = current.drop((h_component, v_component))
+        return current
+
+    @cached_property
+    def normal_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]
+
+        if self.snap_contour_to_grid:
+            coord1 = "xyz"[ax1]
+            coord2 = "xyz"[ax2]
+            # Locations where horizontal paths will be snapped
+            v_bounds = [
+                self.center[ax2] - self.size[ax2] / 2,
+                self.center[ax2] + self.size[ax2] / 2,
+            ]
+            h_snaps = h_horizontal.sel({coord2: v_bounds}, method="nearest").coords[coord2].values
+            # Locations where vertical paths will be snapped
+            h_bounds = [
+                self.center[ax1] - self.size[ax1] / 2,
+                self.center[ax1] + self.size[ax1] / 2,
+            ]
+            v_snaps = h_vertical.sel({coord1: h_bounds}, method="nearest").coords[coord1].values
+
+            bottom_bound = h_snaps[0]
+            top_bound = h_snaps[1]
+            left_bound = v_snaps[0]
+            right_bound = v_snaps[1]
+        else:
+            bottom_bound = self.bounds[0][ax2]
+            top_bound = self.bounds[1][ax2]
+            left_bound = self.bounds[0][ax1]
+            right_bound = self.bounds[1][ax1]
+
+        # Horizontal paths
+        path_size = list(self.size)
+        path_size[ax1] = right_bound - left_bound
+        path_size[ax2] = 0
+        path_center = list(self.center)
+        path_center[ax2] = bottom_bound
+
+        bottom = AxisAlignedPathIntegral(
+            center=path_center,
+            size=path_size,
+            extrapolate_to_endpoints=False,
+            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,
+            snap_path_to_grid=self.snap_contour_to_grid,
+        )
+
+        # Vertical paths
+        path_size = list(self.size)
+        path_size[ax1] = 0
+        path_size[ax2] = top_bound - bottom_bound
+        path_center = list(self.center)
+
+        path_center[ax1] = left_bound
+        left = AxisAlignedPathIntegral(
+            center=path_center,
+            size=path_size,
+            extrapolate_to_endpoints=False,
+            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,
+            snap_path_to_grid=self.snap_contour_to_grid,
+        )
+
+        return (bottom, right, top, left)