From ce9b7db2f131f2e3186d24f931694981442fd2ab Mon Sep 17 00:00:00 2001 From: dmarek Date: Thu, 28 Mar 2024 13:37:15 -0400 Subject: [PATCH] Improved documentation. Revised api of path integrals and impedance calculator to accept monitor data directly and error if the data is the wrong type. Added more testing coverage. Extended utils.py to create random data for time and mode solver monitors. Plus fixed some ruff errors. --- docs/api/plugins/index.rst | 2 + docs/api/plugins/microwave.rst | 13 ++ tests/test_plugins/test_microwave.py | 131 +++++++++++----- tests/utils.py | 73 ++++++++- tidy3d/plugins/microwave/__init__.py | 11 +- .../plugins/microwave/impedance_calculator.py | 42 ++++-- tidy3d/plugins/microwave/path_integrals.py | 140 +++++++++++------- .../smatrix/component_modelers/terminal.py | 13 +- 8 files changed, 306 insertions(+), 119 deletions(-) create mode 100644 docs/api/plugins/microwave.rst diff --git a/docs/api/plugins/index.rst b/docs/api/plugins/index.rst index e25a596631..a70209857d 100644 --- a/docs/api/plugins/index.rst +++ b/docs/api/plugins/index.rst @@ -12,6 +12,7 @@ Plugins adjoint waveguide design + microwave .. include:: /api/plugins/mode_solver.rst @@ -22,3 +23,4 @@ Plugins .. include:: /api/plugins/adjoint.rst .. include:: /api/plugins/waveguide.rst .. include:: /api/plugins/design.rst +.. include:: /api/plugins/microwave.rst \ No newline at end of file diff --git a/docs/api/plugins/microwave.rst b/docs/api/plugins/microwave.rst new file mode 100644 index 0000000000..dc963f691f --- /dev/null +++ b/docs/api/plugins/microwave.rst @@ -0,0 +1,13 @@ +.. currentmodule:: tidy3d + +Microwave +---------------------------- + +.. autosummary:: + :toctree: ../_autosummary/ + :template: module.rst + + tidy3d.plugins.microwave.AxisAlignedPathIntegral + tidy3d.plugins.microwave.VoltageIntegralAxisAligned + tidy3d.plugins.microwave.CurrentIntegralAxisAligned + tidy3d.plugins.microwave.ImpedanceCalculator \ No newline at end of file diff --git a/tests/test_plugins/test_microwave.py b/tests/test_plugins/test_microwave.py index 9ddf8b6174..bad8ff8941 100644 --- a/tests/test_plugins/test_microwave.py +++ b/tests/test_plugins/test_microwave.py @@ -4,7 +4,11 @@ import tidy3d as td from tidy3d import FieldData from tidy3d.constants import ETA_0 -from tidy3d.plugins.microwave import VoltageIntegralAA, CurrentIntegralAA, ImpedanceCalculator +from tidy3d.plugins.microwave import ( + VoltageIntegralAxisAligned, + CurrentIntegralAxisAligned, + ImpedanceCalculator, +) import pydantic.v1 as pydantic from tidy3d.exceptions import DataError from ..utils import run_emulated @@ -17,7 +21,7 @@ FSTOP = 1.5e9 F0 = (FSTART + FSTOP) / 2 FWIDTH = FSTOP - FSTART -FS = np.linspace(FSTART, FSTOP, 5) +FS = np.linspace(FSTART, FSTOP, 3) FIELD_MONITOR = td.FieldMonitor( size=MON_SIZE, fields=FIELDS, name="strip_field", freqs=FS, colocate=False ) @@ -33,8 +37,13 @@ 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" + td.FieldTimeMonitor(center=(0, 0, 0), size=(1, 1, 0), colocate=False, name="field_time"), + td.ModeSolverMonitor( + center=(0, 0, 0), + size=(1, 1, 0), + freqs=FS, + mode_spec=td.ModeSpec(num_modes=2), + name="mode", ), ], sources=[ @@ -44,9 +53,11 @@ source_time=td.GaussianPulse(freq0=F0, fwidth=FWIDTH), ) ], - run_time=2e-12, + run_time=5e-16, ) +SIM_Z_DATA = run_emulated(SIM_Z) + """ Generate the data arrays for testing path integral computations """ @@ -105,13 +116,13 @@ def test_voltage_integral_axes(axis): size = [0, 0, 0] size[axis] = length center = [0, 0, 0] - voltage_integral = VoltageIntegralAA( + voltage_integral = VoltageIntegralAxisAligned( center=center, size=size, + sign="+", ) - sim = SIM_Z - sim_data = run_emulated(sim) - _ = voltage_integral.compute_voltage(sim_data["field"].field_components) + + _ = voltage_integral.compute_voltage(SIM_Z_DATA["field"]) @pytest.mark.parametrize("axis", [0, 1, 2]) @@ -120,13 +131,12 @@ def test_current_integral_axes(axis): size = [length, length, length] size[axis] = 0.0 center = [0, 0, 0] - current_integral = CurrentIntegralAA( + current_integral = CurrentIntegralAxisAligned( center=center, size=size, + sign="+", ) - sim = SIM_Z - sim_data = run_emulated(sim) - _ = current_integral.compute_current(sim_data["field"].field_components) + _ = current_integral.compute_current(SIM_Z_DATA["field"]) def test_voltage_integral_toggles(): @@ -134,16 +144,14 @@ def test_voltage_integral_toggles(): size = [0, 0, 0] size[0] = length center = [0, 0, 0] - voltage_integral = VoltageIntegralAA( + voltage_integral = VoltageIntegralAxisAligned( 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) + _ = voltage_integral.compute_voltage(SIM_Z_DATA["field"]) def test_current_integral_toggles(): @@ -151,16 +159,14 @@ def test_current_integral_toggles(): size = [length, length, length] size[0] = 0.0 center = [0, 0, 0] - current_integral = CurrentIntegralAA( + current_integral = CurrentIntegralAxisAligned( 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) + _ = current_integral.compute_current(SIM_Z_DATA["field"]) def test_voltage_missing_fields(): @@ -168,14 +174,14 @@ def test_voltage_missing_fields(): size = [0, 0, 0] size[1] = length center = [0, 0, 0] - voltage_integral = VoltageIntegralAA( + voltage_integral = VoltageIntegralAxisAligned( center=center, size=size, + sign="+", ) - sim = SIM_Z - sim_data = run_emulated(sim) + with pytest.raises(DataError): - _ = voltage_integral.compute_voltage(sim_data["ExHx"].field_components) + _ = voltage_integral.compute_voltage(SIM_Z_DATA["ExHx"]) def test_current_missing_fields(): @@ -183,14 +189,42 @@ def test_current_missing_fields(): size = [length, length, length] size[0] = 0.0 center = [0, 0, 0] - current_integral = CurrentIntegralAA( + current_integral = CurrentIntegralAxisAligned( center=center, size=size, + sign="+", ) - sim = SIM_Z - sim_data = run_emulated(sim) + with pytest.raises(DataError): - _ = current_integral.compute_current(sim_data["ExHx"].field_components) + _ = current_integral.compute_current(SIM_Z_DATA["ExHx"]) + + +def test_time_monitor_voltage_integral(): + length = 0.5 + size = [0, 0, 0] + size[1] = length + center = [0, 0, 0] + voltage_integral = VoltageIntegralAxisAligned( + center=center, + size=size, + sign="+", + ) + + voltage_integral.compute_voltage(SIM_Z_DATA["field_time"]) + + +def test_mode_solver_monitor_voltage_integral(): + length = 0.5 + size = [0, 0, 0] + size[1] = length + center = [0, 0, 0] + voltage_integral = VoltageIntegralAxisAligned( + center=center, + size=size, + sign="+", + ) + + voltage_integral.compute_voltage(SIM_Z_DATA["mode"]) def test_tiny_voltage_path(): @@ -198,10 +232,11 @@ def test_tiny_voltage_path(): 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) + voltage_integral = VoltageIntegralAxisAligned( + center=center, size=size, sign="+", extrapolate_to_endpoints=True + ) + + _ = voltage_integral.compute_voltage(SIM_Z_DATA["field"]) def test_impedance_calculator(): @@ -209,16 +244,42 @@ def test_impedance_calculator(): _ = ImpedanceCalculator(voltage_integral=None, current_integral=None) +def test_impedance_calculator_on_time_data(): + # Setup path integrals + length = 0.5 + size = [0, length, 0] + size[1] = length + center = [0, 0, 0] + voltage_integral = VoltageIntegralAxisAligned( + center=center, size=size, sign="+", extrapolate_to_endpoints=True + ) + + size = [length, length, 0] + current_integral = CurrentIntegralAxisAligned(center=center, size=size, sign="+") + + # Compute impedance using the tool + Z_calc = ImpedanceCalculator( + voltage_integral=voltage_integral, current_integral=current_integral + ) + _ = Z_calc.compute_impedance(SIM_Z_DATA["field_time"]) + Z_calc = ImpedanceCalculator(voltage_integral=voltage_integral, current_integral=None) + _ = Z_calc.compute_impedance(SIM_Z_DATA["field_time"]) + Z_calc = ImpedanceCalculator(voltage_integral=None, current_integral=current_integral) + _ = Z_calc.compute_impedance(SIM_Z_DATA["field_time"]) + + 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) + voltage_integral = VoltageIntegralAxisAligned( + center=center, size=size, sign="+", extrapolate_to_endpoints=True + ) size = [STRIP_WIDTH * 1.25, STRIP_HEIGHT / 2, 0] center = [0, 0, 0] - current_integral = CurrentIntegralAA(center=center, size=size) + current_integral = CurrentIntegralAxisAligned(center=center, size=size, sign="+") def impedance_of_stripline(width, height): # Assuming no fringing fields, is the same as a parallel plate diff --git a/tests/utils.py b/tests/utils.py index 1fe9932f5e..311cf4f120 100644 --- a/tests/utils.py +++ b/tests/utils.py @@ -10,6 +10,7 @@ from tidy3d.log import _get_level_int from tidy3d.web import BatchData from tidy3d.components.base import Tidy3dBaseModel +from tidy3d import ModeIndexDataArray """ utilities shared between all tests """ np.random.seed(4) @@ -78,7 +79,7 @@ def cartesian_to_unstructured( xyz = [array.x, array.y, array.z] lens = [len(coord) for coord in xyz] - num_len_zero = sum(l == 1 for l in lens) + num_len_zero = sum(length == 1 for length in lens) if num_len_zero == 1: normal_axis = lens.index(1) @@ -245,7 +246,7 @@ def cartesian_to_unstructured( def make_spatial_data( size, bounds, - lims=[0, 1], + lims=(0, 1), seed_data=None, unstructured=False, perturbation=0.1, @@ -875,6 +876,72 @@ def make_field_data(monitor: td.FieldMonitor) -> td.FieldData: **field_cmps, ) + def make_field_time_data(monitor: td.FieldTimeMonitor) -> td.FieldTimeData: + """make a random FieldTimeData from a FieldTimeMonitor.""" + field_cmps = {} + coords = {} + grid = simulation.discretize_monitor(monitor) + tmesh = simulation.tmesh + for field_name in monitor.fields: + spatial_coords_dict = grid[field_name].dict() + + for axis, dim in enumerate("xyz"): + if monitor.size[axis] == 0: + coords[dim] = [monitor.center[axis]] + else: + coords[dim] = np.array(spatial_coords_dict[dim]) + + (idx_begin, idx_end) = monitor.time_inds(tmesh) + tcoords = tmesh[idx_begin:idx_end] + coords["t"] = tcoords + field_cmps[field_name] = make_data( + coords=coords, data_array_type=td.ScalarFieldTimeDataArray, is_complex=False + ) + + return td.FieldTimeData( + monitor=monitor, + symmetry=(0, 0, 0), + symmetry_center=simulation.center, + grid_expanded=grid, + **field_cmps, + ) + + def make_mode_solver_data(monitor: td.ModeSolverMonitor) -> td.ModeSolverData: + """make a random ModeSolverData from a ModeSolverMonitor.""" + field_cmps = {} + coords = {} + grid = simulation.discretize_monitor(monitor) + index_coords = {} + index_coords["f"] = list(monitor.freqs) + index_coords["mode_index"] = np.arange(monitor.mode_spec.num_modes) + index_data_shape = (len(index_coords["f"]), len(index_coords["mode_index"])) + index_data = ModeIndexDataArray( + (1 + 1j) * np.random.random(index_data_shape), coords=index_coords + ) + for field_name in ["Ex", "Ey", "Ez", "Hx", "Hy", "Hz"]: + spatial_coords_dict = grid[field_name].dict() + + for axis, dim in enumerate("xyz"): + if monitor.size[axis] == 0: + coords[dim] = [monitor.center[axis]] + else: + coords[dim] = np.array(spatial_coords_dict[dim]) + + coords["f"] = list(monitor.freqs) + coords["mode_index"] = index_coords["mode_index"] + field_cmps[field_name] = make_data( + coords=coords, data_array_type=td.ScalarModeFieldDataArray, is_complex=True + ) + + return td.ModeSolverData( + monitor=monitor, + symmetry=(0, 0, 0), + symmetry_center=simulation.center, + grid_expanded=grid, + n_complex=index_data, + **field_cmps, + ) + def make_eps_data(monitor: td.PermittivityMonitor) -> td.PermittivityData: """make a random PermittivityData from a PermittivityMonitor.""" field_mnt = td.FieldMonitor(**monitor.dict(exclude={"type", "fields"})) @@ -920,6 +987,8 @@ def make_mode_data(monitor: td.ModeMonitor) -> td.ModeData: MONITOR_MAKER_MAP = { td.FieldMonitor: make_field_data, + td.FieldTimeMonitor: make_field_time_data, + td.ModeSolverMonitor: make_mode_solver_data, td.ModeMonitor: make_mode_data, td.PermittivityMonitor: make_eps_data, td.DiffractionMonitor: make_diff_data, diff --git a/tidy3d/plugins/microwave/__init__.py b/tidy3d/plugins/microwave/__init__.py index a7a6808ad2..5dbb3a1ae6 100644 --- a/tidy3d/plugins/microwave/__init__.py +++ b/tidy3d/plugins/microwave/__init__.py @@ -1,10 +1,15 @@ """ Imports from microwave plugin. """ -from .path_integrals import VoltageIntegralAA, CurrentIntegralAA +from .path_integrals import ( + AxisAlignedPathIntegral, + VoltageIntegralAxisAligned, + CurrentIntegralAxisAligned, +) from .impedance_calculator import ImpedanceCalculator __all__ = [ - "VoltageIntegralAA", - "CurrentIntegralAA", + "AxisAlignedPathIntegral", + "VoltageIntegralAxisAligned", + "CurrentIntegralAxisAligned", "ImpedanceCalculator", ] diff --git a/tidy3d/plugins/microwave/impedance_calculator.py b/tidy3d/plugins/microwave/impedance_calculator.py index 309db82a71..1a03da3269 100644 --- a/tidy3d/plugins/microwave/impedance_calculator.py +++ b/tidy3d/plugins/microwave/impedance_calculator.py @@ -6,40 +6,48 @@ import pydantic.v1 as pd import numpy as np from typing import Optional -from ...components.data.monitor_data import FieldData, FieldTimeData, ModeSolverData +from ...components.data.monitor_data import FieldTimeData from ...components.base import Tidy3dBaseModel -from ...exceptions import ValidationError +from ...exceptions import DataError, ValidationError -from .path_integrals import VoltageIntegralAA, CurrentIntegralAA +from .path_integrals import VoltageIntegralAxisAligned, CurrentIntegralAxisAligned +from .path_integrals import MonitorDataTypes, IntegralResultTypes class ImpedanceCalculator(Tidy3dBaseModel): """Tool for computing the characteristic impedance of a transmission line.""" - voltage_integral: Optional[VoltageIntegralAA] = pd.Field( + voltage_integral: Optional[VoltageIntegralAxisAligned] = pd.Field( ..., title="Voltage Integral", - description="Integral for computing voltage.", + description="Definition of path integral for computing voltage.", ) - current_integral: Optional[CurrentIntegralAA] = pd.Field( + current_integral: Optional[CurrentIntegralAxisAligned] = pd.Field( ..., - title="", - description="Integral for computing current.", + title="Current Integral", + description="Definition of contour integral for computing current.", ) - def compute_impedance(self, em_field: FieldData | ModeSolverData | FieldTimeData): + def compute_impedance(self, em_field: MonitorDataTypes) -> IntegralResultTypes: + """Compute impedance for the supplied ``em_field`` using ``voltage_integral`` and + ``current_integral``. If only a single integral has been defined, impedance is + computed using the total flux in ``em_field``.""" + if not isinstance(em_field, MonitorDataTypes): + raise DataError("'em_field' type not supported by impedance calculator.") + # 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) + voltage = self.voltage_integral.compute_voltage(em_field) if self.current_integral: - current = self.current_integral.compute_current(em_field.field_components) + current = self.current_integral.compute_current(em_field) # 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)) + # 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,11 +61,13 @@ def compute_impedance(self, em_field: FieldData | ModeSolverData | FieldTimeData else: current = np.conj(2 * flux / voltage) - return voltage / current + impedance = voltage / current + return impedance @pd.validator("current_integral", always=True) def check_voltage_or_current(cls, val, values): - """Ensure that 'voltage_integral' and/or 'current_integral' were provided.""" + """Raise validation error if both ``voltage_integral`` and ``current_integral`` + were not provided.""" if not values.get("voltage_integral") and not val: raise ValidationError( "Atleast one of 'voltage_integral' or 'current_integral' must be provided." diff --git a/tidy3d/plugins/microwave/path_integrals.py b/tidy3d/plugins/microwave/path_integrals.py index af75306115..b4f6a2462b 100644 --- a/tidy3d/plugins/microwave/path_integrals.py +++ b/tidy3d/plugins/microwave/path_integrals.py @@ -4,28 +4,39 @@ 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 AbstractSpatialDataArray -from ...components.base import cached_property +from abc import ABC, abstractmethod +from typing import List, Tuple, Union + +from ...components.data.monitor_data import FieldData, FieldTimeData, ModeSolverData +from ...components.data.data_array import ( + ScalarFieldDataArray, + ScalarFieldTimeDataArray, + ScalarModeFieldDataArray, +) +from ...components.data.data_array import FreqDataArray, TimeDataArray, FreqModeDataArray +from ...components.base import cached_property, Tidy3dBaseModel from ...components.types import Axis, Direction from ...components.geometry.base import Box from ...components.validators import assert_line, assert_plane -from ...exceptions import DataError +from ...exceptions import Tidy3dError, DataError +MonitorDataTypes = Union[FieldData, FieldTimeData, ModeSolverData] +EMScalarFieldType = Union[ScalarFieldDataArray, ScalarFieldTimeDataArray, ScalarModeFieldDataArray] +IntegralResultTypes = Union[FreqDataArray, FreqModeDataArray, TimeDataArray] -class AbstractAxesRH: + +class AbstractAxesRH(Tidy3dBaseModel, ABC): """Represents an axis-aligned right-handed coordinate system with one axis preferred.""" @cached_property + @abstractmethod def main_axis(self) -> Axis: - """Subclasses should implement this method.""" - raise NotImplementedError() + """Get the preferred axis.""" @cached_property def remaining_axes(self) -> Tuple[Axis, Axis]: - """Axes in plane, ordered to maintain a right-handed coordinate system""" - axes = [0, 1, 2] + """Get in-plane axes, ordered to maintain a right-handed coordinate system.""" + axes: List[Axis] = [0, 1, 2] axes.pop(self.main_axis) if self.main_axis == 1: return (axes[1], axes[0]) @@ -34,29 +45,30 @@ def remaining_axes(self) -> Tuple[Axis, Axis]: class AxisAlignedPathIntegral(AbstractAxesRH, Box): - """Class for defining the simplest type of path integral which is aligned with Cartesian axes.""" + """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", + title="Extrapolate to Endpoints", 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. ", + "the field is likely discontinuous. When this field is ``True``, fields that are outside and on the bounds " + "of the integral are ignored. Should be enabled when computing voltage between two conductors.", ) snap_path_to_grid: bool = pd.Field( False, - title="Snap path to grid", + 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.", + "a field. When this field is ``True``, the integration path will be snapped to the grid.", ) - def compute_integral(self, scalar_field: AbstractSpatialDataArray): - """Computes the defined integral given the input `scalar_field`.""" + def compute_integral(self, scalar_field: EMScalarFieldType) -> IntegralResultTypes: + """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") + raise DataError("Scalar field does not cover the integration domain.") coord = "xyz"[self.main_axis] scalar_field = self._get_field_along_path(scalar_field) @@ -87,12 +99,17 @@ def compute_integral(self, scalar_field: AbstractSpatialDataArray): scalar_field = scalar_field.interp( coords_interp, method=method, kwargs={"fill_value": "extrapolate"} ) - return scalar_field.integrate(coord=coord) + result = scalar_field.integrate(coord=coord) + if isinstance(scalar_field, ScalarFieldDataArray): + return FreqDataArray(data=result.data, coords=result.coords) + elif isinstance(scalar_field, ScalarFieldTimeDataArray): + return TimeDataArray(data=result.data, coords=result.coords) + else: + assert isinstance(scalar_field, ScalarModeFieldDataArray) + return FreqModeDataArray(data=result.data, coords=result.coords) - def _get_field_along_path( - self, scalar_field: AbstractSpatialDataArray - ) -> AbstractSpatialDataArray: - """Returns a selection of the input `scalar_field` ready for integration.""" + def _get_field_along_path(self, scalar_field: EMScalarFieldType) -> EMScalarFieldType: + """Returns a selection of the input ``scalar_field`` ready for integration.""" axis1 = self.remaining_axes[0] axis2 = self.remaining_axes[1] coord1 = "xyz"[axis1] @@ -126,68 +143,74 @@ def _get_field_along_path( scalar_field = scalar_field.interp( coord2dict, method="linear", kwargs={"bounds_error": True} ) - + # Remove unneeded coordinates + scalar_field = scalar_field.reset_coords(drop=True) return scalar_field @cached_property 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 + for index, value in enumerate(self.size): + if value != 0: + return index + raise Tidy3dError("Failed to identify axis.") -class VoltageIntegralAA(AxisAlignedPathIntegral): +class VoltageIntegralAxisAligned(AxisAlignedPathIntegral): """Class for computing the voltage between two points defined by an axis-aligned line.""" sign: Direction = pd.Field( - "+", - title="Direction of path integral.", + ..., + 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): + def compute_voltage(self, em_field: MonitorDataTypes) -> IntegralResultTypes: """Compute voltage along path defined by a line.""" + if not isinstance(em_field, MonitorDataTypes): + raise DataError("'em_field' type not supported.") 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: - raise DataError(f"field_name '{field_name}' not found") - e_field = em_field[f"E{e_component}"] - integral = self.compute_integral(e_field) + if field_name not in em_field.field_components: + raise DataError(f"'field_name' '{field_name}' not found.") + e_field = em_field.field_components[field_name] + # V = -integral(E) + voltage = self.compute_integral(e_field) - voltage = -integral - if self.sign == "-": + if self.sign == "+": voltage *= -1 # Return data array of voltage while keeping coordinates of frequency|time|mode index return voltage -class CurrentIntegralAA(AbstractAxesRH, Box): - """Class for computing conduction current via Ampere's Circuital Law on an axis-aligned loop.""" +class CurrentIntegralAxisAligned(AbstractAxesRH, Box): + """Class for computing conduction current via Ampère's circuital law on an axis-aligned loop.""" _plane_validator = assert_plane() sign: Direction = pd.Field( - "+", - title="Direction of contour integral", + ..., + 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.", + 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="Snap contour to grid", - description="This parameter is passed to `AxisAlignedPathIntegral` objects when computing the contour integral.", + 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): + def compute_current(self, em_field: MonitorDataTypes) -> IntegralResultTypes: """Compute current flowing in loop defined by the outer edge of a rectangle.""" - + if not isinstance(em_field, MonitorDataTypes): + raise DataError("'em_field' type not supported.") ax1 = self.remaining_axes[0] ax2 = self.remaining_axes[1] h_component = "xyz"[ax1] @@ -195,12 +218,12 @@ def compute_current(self, em_field: FieldDataset | ModeSolverDataset | FieldTime 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}"] + if h_field_name not in em_field.field_components: + raise DataError(f"'field_name' '{h_field_name}' not found.") + if v_field_name not in em_field.field_components: + raise DataError(f"'field_name' '{v_field_name}' not found.") + h_horizontal = em_field.field_components[h_field_name] + h_vertical = em_field.field_components[v_field_name] # Decompose contour into path integrals (bottom, right, top, left) = self._to_path_integrals(h_horizontal, h_vertical) @@ -214,17 +237,20 @@ def compute_current(self, em_field: FieldDataset | ModeSolverDataset | FieldTime if self.sign == "-": current *= -1 - # Return data array of current while keeping coordinates of frequency|time|mode index - current = current.drop_vars((h_component, v_component)) + return current @cached_property 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 + for index, value in enumerate(self.size): + if value == 0: + return index + raise Tidy3dError("Failed to identify axis.") def _to_path_integrals(self, h_horizontal, h_vertical) -> Tuple[AxisAlignedPathIntegral, ...]: + """Returns four ``AxisAlignedPathIntegral`` instances, which represent a contour + integral around the surface defined by ``self.size``.""" ax1 = self.remaining_axes[0] ax2 = self.remaining_axes[1] diff --git a/tidy3d/plugins/smatrix/component_modelers/terminal.py b/tidy3d/plugins/smatrix/component_modelers/terminal.py index e2faec7118..6d5aa88277 100644 --- a/tidy3d/plugins/smatrix/component_modelers/terminal.py +++ b/tidy3d/plugins/smatrix/component_modelers/terminal.py @@ -23,7 +23,7 @@ from .base import AbstractComponentModeler, FWIDTH_FRAC from ..ports.lumped import LumpedPortDataArray, LumpedPort -from ...microwave import VoltageIntegralAA, CurrentIntegralAA +from ...microwave import VoltageIntegralAxisAligned, CurrentIntegralAxisAligned class TerminalComponentModeler(AbstractComponentModeler): @@ -215,14 +215,14 @@ def port_voltage(port: LumpedPort, sim_data: SimulationData) -> xr.DataArray: size = list(port.size) size[port.current_axis] = 0 - voltage_integral = VoltageIntegralAA( + voltage_integral = VoltageIntegralAxisAligned( center=port.center, size=size, extrapolate_to_endpoints=True, snap_path_to_grid=True, sign="+", ) - voltage = voltage_integral.compute_voltage(field_data.field_components) + voltage = voltage_integral.compute_voltage(field_data) # Return data array of voltage with coordinates of frequency return voltage @@ -242,7 +242,8 @@ def port_current(port: LumpedPort, sim_data: SimulationData) -> xr.DataArray: # Get h field tangent to resistive sheet h_component = "xyz"[port.current_axis] inject_component = "xyz"[port.injection_axis] - field_data = sim_data[f"{port.name}_H{h_component}"].field_components + monitor_data = sim_data[f"{port.name}_H{h_component}"] + field_data = monitor_data.field_components h_field = field_data[f"H{h_component}"] # Coordinates as numpy array for h_field along curren and injection axis h_coords_along_current = h_field.coords[h_component].values @@ -290,14 +291,14 @@ def port_current(port: LumpedPort, sim_data: SimulationData) -> xr.DataArray: size[port.injection_axis] = dcap # H field is continuous at integral bounds, so extrapolation is turned off - I_integral = CurrentIntegralAA( + I_integral = CurrentIntegralAxisAligned( center=center, size=size, sign="+", extrapolate_to_endpoints=False, snap_contour_to_grid=True, ) - return I_integral.compute_current(field_data) + return I_integral.compute_current(monitor_data) def port_ab(port: LumpedPort, sim_data: SimulationData): """Helper to compute the port incident and reflected power waves."""