adding microwave plugin with path integrals

added utility for computing impedance from EM fields

reused path integrals in smatrix plugin

CHANGELOG.md

@@ -12,6 +12,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 ## [2.6.1] - 2024-03-07
 
 ### Added
+- Introduced `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.

tests/test_plugins/test_microwave.py

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

tidy3d/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,65 @@
+"""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 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 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)
+
+        return voltage / current
+
+    @pd.validator("current_integral", always=True)
+    def check_voltage_or_current(cls, val, values):
+        """Ensure that 'voltage_integral' and/or 'current_integral' were provided."""
+        if not values.get("voltage_integral") and not val:
+            raise ValidationError(
+                "Atleast one of 'voltage_integral' or 'current_integral' must be provided."
+            )
+        return val