From 2797fe6109c002351a590766256d7a4ce17d069c Mon Sep 17 00:00:00 2001 From: dmarek Date: Fri, 8 Mar 2024 09:24:17 -0500 Subject: [PATCH] adding microwave plugin with path integrals for computing voltage/current/impedance reused path integrals in smatrix plugin --- CHANGELOG.md | 1 + docs/api/plugins/index.rst | 2 + docs/api/plugins/microwave.rst | 13 + tests/test_plugins/test_microwave.py | 305 ++++++++++++++++ tests/utils.py | 73 +++- tidy3d/components/validators.py | 13 + tidy3d/plugins/microwave/__init__.py | 15 + .../plugins/microwave/impedance_calculator.py | 75 ++++ tidy3d/plugins/microwave/path_integrals.py | 325 ++++++++++++++++++ .../smatrix/component_modelers/terminal.py | 154 ++++----- tidy3d/plugins/smatrix/ports/lumped.py | 5 +- 11 files changed, 889 insertions(+), 92 deletions(-) create mode 100644 docs/api/plugins/microwave.rst create mode 100644 tests/test_plugins/test_microwave.py create mode 100644 tidy3d/plugins/microwave/__init__.py create mode 100644 tidy3d/plugins/microwave/impedance_calculator.py create mode 100644 tidy3d/plugins/microwave/path_integrals.py diff --git a/CHANGELOG.md b/CHANGELOG.md index c1615a256..973e13fc2 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -15,6 +15,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ## [2.6.1] - 2024-03-07 ### Added +- Introduced the `microwave` plugin which includes `ImpedanceCalculator` for computing the characteristic impedance of transmission lines. - `Simulation` now accepts `LumpedElementType`, which currently only supports the `LumpedResistor` type. `LumpedPort` together with `LumpedResistor` make up the new `TerminalComponentModeler` in the `smatrix` plugin. - Uniaxial medium Lithium niobate to material library. - Added support for conformal mesh methods near PEC structures that can be specified through the field `pec_conformal_mesh_spec` in the `Simulation` class. diff --git a/docs/api/plugins/index.rst b/docs/api/plugins/index.rst index e25a59663..a70209857 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 000000000..dc963f691 --- /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 new file mode 100644 index 000000000..bad8ff894 --- /dev/null +++ b/tests/test_plugins/test_microwave.py @@ -0,0 +1,305 @@ +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 ( + VoltageIntegralAxisAligned, + CurrentIntegralAxisAligned, + 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, 3) +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.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=[ + td.PointDipole( + center=(0, 0, 0), + polarization="Ex", + source_time=td.GaussianPulse(freq0=F0, fwidth=FWIDTH), + ) + ], + run_time=5e-16, +) + +SIM_Z_DATA = run_emulated(SIM_Z) + +""" 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 = VoltageIntegralAxisAligned( + center=center, + size=size, + sign="+", + ) + + _ = voltage_integral.compute_voltage(SIM_Z_DATA["field"]) + + +@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 = CurrentIntegralAxisAligned( + center=center, + size=size, + sign="+", + ) + _ = current_integral.compute_current(SIM_Z_DATA["field"]) + + +def test_voltage_integral_toggles(): + length = 0.5 + size = [0, 0, 0] + size[0] = length + center = [0, 0, 0] + voltage_integral = VoltageIntegralAxisAligned( + center=center, + size=size, + extrapolate_to_endpoints=True, + snap_path_to_grid=True, + sign="-", + ) + _ = voltage_integral.compute_voltage(SIM_Z_DATA["field"]) + + +def test_current_integral_toggles(): + length = 0.5 + size = [length, length, length] + size[0] = 0.0 + center = [0, 0, 0] + current_integral = CurrentIntegralAxisAligned( + center=center, + size=size, + extrapolate_to_endpoints=True, + snap_contour_to_grid=True, + sign="-", + ) + _ = current_integral.compute_current(SIM_Z_DATA["field"]) + + +def test_voltage_missing_fields(): + length = 0.5 + size = [0, 0, 0] + size[1] = length + center = [0, 0, 0] + voltage_integral = VoltageIntegralAxisAligned( + center=center, + size=size, + sign="+", + ) + + with pytest.raises(DataError): + _ = voltage_integral.compute_voltage(SIM_Z_DATA["ExHx"]) + + +def test_current_missing_fields(): + length = 0.5 + size = [length, length, length] + size[0] = 0.0 + center = [0, 0, 0] + current_integral = CurrentIntegralAxisAligned( + center=center, + size=size, + sign="+", + ) + + with pytest.raises(DataError): + _ = 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(): + length = 0.02 + size = [0, 0, 0] + size[1] = length + center = [0, 0, 0] + voltage_integral = VoltageIntegralAxisAligned( + center=center, size=size, sign="+", extrapolate_to_endpoints=True + ) + + _ = voltage_integral.compute_voltage(SIM_Z_DATA["field"]) + + +def test_impedance_calculator(): + with pytest.raises(pydantic.ValidationError): + _ = 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 = 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 = CurrentIntegralAxisAligned(center=center, size=size, sign="+") + + 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)) diff --git a/tests/utils.py b/tests/utils.py index 1fe9932f5..311cf4f12 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/components/validators.py b/tidy3d/components/validators.py index a2d9b7162..4b163f23c 100644 --- a/tidy3d/components/validators.py +++ b/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""" diff --git a/tidy3d/plugins/microwave/__init__.py b/tidy3d/plugins/microwave/__init__.py new file mode 100644 index 000000000..5dbb3a1ae --- /dev/null +++ b/tidy3d/plugins/microwave/__init__.py @@ -0,0 +1,15 @@ +""" Imports from microwave plugin. """ + +from .path_integrals import ( + AxisAlignedPathIntegral, + VoltageIntegralAxisAligned, + CurrentIntegralAxisAligned, +) +from .impedance_calculator import ImpedanceCalculator + +__all__ = [ + "AxisAlignedPathIntegral", + "VoltageIntegralAxisAligned", + "CurrentIntegralAxisAligned", + "ImpedanceCalculator", +] diff --git a/tidy3d/plugins/microwave/impedance_calculator.py b/tidy3d/plugins/microwave/impedance_calculator.py new file mode 100644 index 000000000..cc7464f7f --- /dev/null +++ b/tidy3d/plugins/microwave/impedance_calculator.py @@ -0,0 +1,75 @@ +"""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 ...exceptions import DataError, ValidationError + +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[VoltageIntegralAxisAligned] = pd.Field( + ..., + title="Voltage Integral", + description="Definition of path integral for computing voltage.", + ) + + current_integral: Optional[CurrentIntegralAxisAligned] = pd.Field( + ..., + title="Current Integral", + description="Definition of contour integral for computing current.", + ) + + 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, (FieldData, FieldTimeData, ModeSolverData)): + 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) + if self.current_integral: + 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)) + 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) + + impedance = voltage / current + return impedance + + @pd.validator("current_integral", always=True) + def check_voltage_or_current(cls, val, values): + """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." + ) + return val diff --git a/tidy3d/plugins/microwave/path_integrals.py b/tidy3d/plugins/microwave/path_integrals.py new file mode 100644 index 000000000..7044a0d5d --- /dev/null +++ b/tidy3d/plugins/microwave/path_integrals.py @@ -0,0 +1,325 @@ +"""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 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 Tidy3dError, DataError + +MonitorDataTypes = Union[FieldData, FieldTimeData, ModeSolverData] +EMScalarFieldType = Union[ScalarFieldDataArray, ScalarFieldTimeDataArray, ScalarModeFieldDataArray] +IntegralResultTypes = Union[FreqDataArray, FreqModeDataArray, TimeDataArray] + + +class AbstractAxesRH(Tidy3dBaseModel, ABC): + """Represents an axis-aligned right-handed coordinate system with one axis preferred.""" + + @cached_property + @abstractmethod + def main_axis(self) -> Axis: + """Get the preferred axis.""" + + @cached_property + def remaining_axes(self) -> Tuple[Axis, Axis]: + """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]) + 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 a material interface, " + "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", + description="It might be desireable to integrate exactly along the Yee grid associated with " + "a field. When this field is ``True``, the integration path will be snapped to the grid.", + ) + + 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.") + coord = "xyz"[self.main_axis] + + scalar_field = self._get_field_along_path(scalar_field) + # Get the boundaries + 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 + 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, 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=method, kwargs={"fill_value": "extrapolate"} + ) + 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: 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] + 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} + ) + # 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.""" + for index, value in enumerate(self.size): + if value != 0: + return index + raise Tidy3dError("Failed to identify axis.") + + +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", + description="Positive indicates V=Vb-Va where position b has a larger coordinate along the axis of integration.", + ) + + def compute_voltage(self, em_field: MonitorDataTypes) -> IntegralResultTypes: + """Compute voltage along path defined by a line.""" + if not isinstance(em_field, (FieldData, FieldTimeData, ModeSolverData)): + 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.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) + + if self.sign == "+": + voltage *= -1 + # Return data array of voltage while keeping coordinates of frequency|time|mode index + return voltage + + +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", + 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="Snap Contour to Grid", + description="This parameter is passed to ``AxisAlignedPathIntegral`` objects when computing the contour integral.", + ) + + 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, (FieldData, FieldTimeData, ModeSolverData)): + raise DataError("'em_field' type not supported.") + 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.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) + + 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 current + + @cached_property + def main_axis(self) -> Axis: + """Axis normal to loop""" + 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] + + 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=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=self.extrapolate_to_endpoints, + 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=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=self.extrapolate_to_endpoints, + snap_path_to_grid=self.snap_contour_to_grid, + ) + + return (bottom, right, top, left) diff --git a/tidy3d/plugins/smatrix/component_modelers/terminal.py b/tidy3d/plugins/smatrix/component_modelers/terminal.py index 326466621..6d5aa8827 100644 --- a/tidy3d/plugins/smatrix/component_modelers/terminal.py +++ b/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 VoltageIntegralAxisAligned, CurrentIntegralAxisAligned + class TerminalComponentModeler(AbstractComponentModeler): """Tool for modeling two-terminal multiport devices and computing port parameters @@ -194,40 +196,33 @@ 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.""" 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 = VoltageIntegralAxisAligned( + 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) # Return data array of voltage with coordinates of frequency return voltage @@ -246,83 +241,64 @@ 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] + 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 + 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) + (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 + # 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 = CurrentIntegralAxisAligned( + 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(monitor_data) def port_ab(port: LumpedPort, sim_data: SimulationData): """Helper to compute the port incident and reflected power waves.""" diff --git a/tidy3d/plugins/smatrix/ports/lumped.py b/tidy3d/plugins/smatrix/ports/lumped.py index 961e478fd..9e6089318 100644 --- a/tidy3d/plugins/smatrix/ports/lumped.py +++ b/tidy3d/plugins/smatrix/ports/lumped.py @@ -1,10 +1,12 @@ """Class and custom data array for representing a scattering matrix port based on lumped circuit elements.""" + import pydantic.v1 as pd import numpy as np from typing import Optional from ....constants import OHM from ....components.geometry.base import Box +from ....components.geometry.utils_2d import increment_float from ....components.types import Complex, FreqArray, Axis from ....components.base import cached_property from ....components.lumped_element import LumpedResistor @@ -159,8 +161,9 @@ def to_current_monitor(self, freqs: FreqArray, snap_center: float) -> FieldMonit h_cap_component = "xyz"[self.injection_axis] # Size of current monitor needs to encompass the current carrying 2D sheet # Needs to have a nonzero thickness so a closed loop of gridpoints around the 2D sheet can be formed + dl = 2 * (increment_float(center[self.injection_axis], 1.0) - center[self.injection_axis]) current_mon_size = list(self.size) - current_mon_size[self.injection_axis] = 1e-10 + current_mon_size[self.injection_axis] = dl current_mon_size[self.voltage_axis] = 0.0 # Create a current monitor return FieldMonitor(