diff --git a/docs/api/plugins/smatrix.rst b/docs/api/plugins/smatrix.rst index 4cabba2080..39114121d4 100644 --- a/docs/api/plugins/smatrix.rst +++ b/docs/api/plugins/smatrix.rst @@ -11,6 +11,7 @@ Scattering Matrix Calculator tidy3d.plugins.smatrix.Port tidy3d.plugins.smatrix.ModalPortDataArray tidy3d.plugins.smatrix.TerminalComponentModeler + tidy3d.plugins.smatrix.TerminalPortDataArray tidy3d.plugins.smatrix.LumpedPort - tidy3d.plugins.smatrix.LumpedPortDataArray - tidy3d.plugins.smatrix.CoaxialLumpedPort \ No newline at end of file + tidy3d.plugins.smatrix.CoaxialLumpedPort + tidy3d.plugins.smatrix.WavePort \ No newline at end of file diff --git a/tests/test_plugins/terminal_component_modeler_def.py b/tests/test_plugins/terminal_component_modeler_def.py index de762b8104..697688eb08 100644 --- a/tests/test_plugins/terminal_component_modeler_def.py +++ b/tests/test_plugins/terminal_component_modeler_def.py @@ -1,9 +1,13 @@ +from typing import Union + import numpy as np import tidy3d as td +import tidy3d.plugins.microwave as microwave from tidy3d.plugins.smatrix import ( CoaxialLumpedPort, LumpedPort, TerminalComponentModeler, + WavePort, ) # Microstrip dimensions @@ -25,8 +29,8 @@ freq_stop = 10e9 # Coaxial dimensions -Rinner = 0.2768 * 1e-3 -Router = 1.0 * 1e-3 +Rinner = 0.2768 * 1e3 +Router = 1.0 * 1e3 def make_simulation(planar_pec: bool, length: float = None, auto_grid: bool = True): @@ -176,7 +180,13 @@ def make_coaxial_simulation(length: float = None, auto_grid: bool = True): # Spatial grid specification if auto_grid: - grid_spec = td.GridSpec.auto(min_steps_per_wvl=10, wavelength=td.C_0 / freq_stop) + mesh_override = td.MeshOverrideStructure( + geometry=td.Box(center=(0, 0, 0), size=(2 * Router, 2 * Router, coax_length)), + dl=(50, 50, 10e3), + ) + grid_spec = td.GridSpec.auto( + min_steps_per_wvl=10, wavelength=td.C_0 / freq_stop, override_structures=[mesh_override] + ) else: grid_spec = td.GridSpec.uniform(wavelength0 / 11) @@ -246,6 +256,10 @@ def make_coaxial_component_modeler( length: float = None, port_refinement: bool = True, auto_grid: bool = True, + port_types: tuple[Union[CoaxialLumpedPort, WavePort], Union[CoaxialLumpedPort, WavePort]] = ( + CoaxialLumpedPort, + CoaxialLumpedPort, + ), **kwargs, ): if length: @@ -255,24 +269,55 @@ def make_coaxial_component_modeler( sim = make_coaxial_simulation(length=coax_length, auto_grid=auto_grid) - center_src1 = [0, 0, -coax_length / 2] - - port_cells = None - if port_refinement: - port_cells = 21 + def make_port(center, direction, type, name) -> Union[CoaxialLumpedPort, WavePort]: + if type is CoaxialLumpedPort: + port_cells = None + if port_refinement: + port_cells = 21 + port = CoaxialLumpedPort( + center=center, + outer_diameter=2 * Router, + inner_diameter=2 * Rinner, + normal_axis=2, + direction=direction, + name=name, + num_grid_cells=port_cells, + impedance=reference_impedance, + ) + else: + mean_radius = (Router + Rinner) / 2 + voltage_center = list(center) + voltage_center[0] += mean_radius + voltage_size = [Router - Rinner, 0, 0] + + port = WavePort( + center=center, + size=[2 * Router, 2 * Router, 0], + direction=direction, + name=name, + mode_spec=td.ModeSpec(num_modes=1), + mode_index=0, + voltage_integral=microwave.VoltageIntegralAxisAligned( + center=voltage_center, + size=voltage_size, + extrapolate_to_endpoints=True, + snap_path_to_grid=True, + sign="+", + ), + current_integral=microwave.CustomCurrentIntegral2D.from_circular_path( + center=center, + radius=mean_radius, + num_points=41, + normal_axis=2, + clockwise=False, + ), + ) + return port - port_1 = CoaxialLumpedPort( - center=center_src1, - outer_diameter=2 * Router, - inner_diameter=2 * Rinner, - normal_axis=2, - direction="+", - name="coax_port_1", - num_grid_cells=port_cells, - impedance=reference_impedance, - ) + center_src1 = [0, 0, -coax_length / 2] + port_1 = make_port(center_src1, direction="+", type=port_types[0], name="coax_port_1") center_src2 = [0, 0, coax_length / 2] - port_2 = port_1.updated_copy(name="coax_port_2", center=center_src2, direction="-") + port_2 = make_port(center_src2, direction="-", type=port_types[1], name="coax_port_2") ports = [port_1, port_2] freqs = np.linspace(freq_start, freq_stop, 100) diff --git a/tests/test_plugins/test_microwave.py b/tests/test_plugins/test_microwave.py index ca0302ce17..a978853e91 100644 --- a/tests/test_plugins/test_microwave.py +++ b/tests/test_plugins/test_microwave.py @@ -1,5 +1,6 @@ """Test the microwave plugin.""" +import matplotlib.pyplot as plt import numpy as np import pydantic.v1 as pydantic import pytest @@ -33,6 +34,8 @@ ) STRIP_WIDTH = 1.5 STRIP_HEIGHT = 0.5 +COAX_R1 = 0.04 +COAX_R2 = 0.5 SIM_Z = td.Simulation( size=(2, 1, 1), @@ -103,6 +106,51 @@ def make_stripline_scalar_field_data_array(grid_key: str): return td.ScalarFieldDataArray(values, coords=dict(x=XS, y=YS, z=ZS, f=FS)) +def make_coaxial_field_data_array(grid_key: str): + """Populate FIELD_MONITOR with a coaxial transmission line mode.""" + + # Get a normalized electric field that represents the electric field within a coaxial cable transmission line. + def compute_coax_radial_electric(rin, rout, x, y, is_x): + # Radial distance + r = np.sqrt((x) ** 2 + (y) ** 2) + # Remove division by 0 + r_valid = np.where(r == 0.0, 1, r) + # Compute current density so that the total current + # is 1 flowing through surfaces of constant r + # Extra r is for changing to Cartesian coordinates + denominator = 2 * np.pi * r_valid**2 + if is_x: + Exy = np.where(r <= rin, 0, (x / denominator)) + Exy = np.where(r >= rout, 0, Exy) + else: + Exy = np.where(r <= rin, 0, (y / denominator)) + Exy = np.where(r >= rout, 0, Exy) + return Exy + + 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))) + + XGRID = np.broadcast_to(XGRID, values.shape) + YGRID = np.broadcast_to(YGRID, values.shape) + + is_x = grid_key[1] == "x" + if grid_key[0] == "E": + field = compute_coax_radial_electric(COAX_R1, COAX_R2, XGRID, YGRID, is_x) + else: + field = compute_coax_radial_electric(COAX_R1, COAX_R2, XGRID, YGRID, not is_x) + # H field is perpendicular and oriented for positive z propagation + # We want to compute Hx which is -Ey/eta0, Hy which is Ex/eta0 + if is_x: + field /= -ETA_0 + else: + field /= ETA_0 + + return td.ScalarFieldDataArray(field, coords=dict(x=XS, y=YS, z=ZS, f=FS)) + + def make_field_data(): return FieldData( monitor=FIELD_MONITOR, @@ -116,6 +164,19 @@ def make_field_data(): ) +def make_coax_field_data(): + return FieldData( + monitor=FIELD_MONITOR, + Ex=make_coaxial_field_data_array("Ex"), + Ey=make_coaxial_field_data_array("Ey"), + Hx=make_coaxial_field_data_array("Hx"), + Hy=make_coaxial_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): """Check VoltageIntegralAxisAligned runs.""" @@ -210,6 +271,15 @@ def test_current_missing_fields(): with pytest.raises(DataError): _ = current_integral.compute_current(SIM_Z_DATA["ExHx"]) + current_integral = CurrentIntegralAxisAligned( + center=center, + size=[length, length, 0], + sign="+", + ) + + with pytest.raises(DataError): + _ = current_integral.compute_current(SIM_Z_DATA["ExHx"]) + def test_time_monitor_voltage_integral(): """Check VoltageIntegralAxisAligned runs on time domain data.""" @@ -417,6 +487,20 @@ def test_fields_missing_custom_current_integral(): current_integral.compute_current(SIM_Z_DATA["ExHx"]) +def test_wrong_monitor_data(): + """Check that the function arguments to integrals are correctly validated.""" + length = 0.5 + size = [0, 0, 0] + size[1] = length + # Make box + vertices = [(0.2, -0.2), (0.2, 0.2), (-0.2, 0.2), (-0.2, -0.2), (0.2, -0.2)] + current_integral = CustomCurrentIntegral2D(axis=2, position=0, vertices=vertices) + with pytest.raises(DataError): + current_integral.compute_current(SIM_Z.sources[0].source_time) + with pytest.raises(DataError): + current_integral.compute_current(em_field=SIM_Z.sources[0].source_time) + + def test_time_monitor_custom_current_integral(): length = 0.5 size = [0, 0, 0] @@ -435,3 +519,109 @@ def test_mode_solver_custom_current_integral(): vertices = [(0.2, -0.2), (0.2, 0.2), (-0.2, 0.2), (-0.2, -0.2), (0.2, -0.2)] current_integral = CustomCurrentIntegral2D(axis=2, position=0, vertices=vertices) current_integral.compute_current(SIM_Z_DATA["mode"]) + + +def test_custom_current_integral_normal_y(): + current_integral = CustomCurrentIntegral2D.from_circular_path( + center=(0, 0, 0), radius=0.4, num_points=31, normal_axis=1, clockwise=False + ) + current_integral.compute_current(SIM_Z_DATA["field"]) + + +def test_custom_path_integral_accuracy(): + """Test the accuracy of the custom path integral.""" + field_data = make_coax_field_data() + + current_integral = CustomCurrentIntegral2D.from_circular_path( + center=(0, 0, 0), radius=0.4, num_points=31, normal_axis=2, clockwise=False + ) + current = current_integral.compute_current(field_data) + assert np.allclose(current.values, 1.0 / ETA_0, rtol=0.01) + + current_integral = CustomCurrentIntegral2D.from_circular_path( + center=(0, 0, 0), radius=0.4, num_points=31, normal_axis=2, clockwise=True + ) + current = current_integral.compute_current(field_data) + assert np.allclose(current.values, -1.0 / ETA_0, rtol=0.01) + + +def test_impedance_accuracy_on_coaxial(): + """Test the accuracy of the ImpedanceCalculator.""" + field_data = make_coax_field_data() + # Setup path integrals + mean_radius = (COAX_R2 + COAX_R1) * 0.5 + size = [COAX_R2 - COAX_R1, 0, 0] + center = [mean_radius, 0, 0] + + voltage_integral = VoltageIntegralAxisAligned( + center=center, size=size, sign="-", extrapolate_to_endpoints=True, snap_path_to_grid=True + ) + current_integral = CustomCurrentIntegral2D.from_circular_path( + center=(0, 0, 0), radius=mean_radius, num_points=31, normal_axis=2, clockwise=False + ) + + def impedance_of_coaxial_cable(r1, r2, wave_impedance=td.ETA_0): + return np.log(r2 / r1) * wave_impedance / 2 / np.pi + + Z_analytic = impedance_of_coaxial_cable(COAX_R1, COAX_R2) + + # Compute impedance using the tool + Z_calculator = ImpedanceCalculator( + voltage_integral=voltage_integral, current_integral=current_integral + ) + Z_calc = Z_calculator.compute_impedance(field_data) + assert np.allclose(Z_calc, Z_analytic, rtol=0.04) + + +def test_path_integral_plotting(): + """Test that all types of path integrals correctly plot themselves.""" + + mean_radius = (COAX_R2 + COAX_R1) * 0.5 + size = [COAX_R2 - COAX_R1, 0, 0] + center = [mean_radius, 0, 0] + + voltage_integral = VoltageIntegralAxisAligned( + center=center, size=size, sign="-", extrapolate_to_endpoints=True, snap_path_to_grid=True + ) + + current_integral = CustomCurrentIntegral2D.from_circular_path( + center=(0, 0, 0), radius=0.4, num_points=31, normal_axis=2, clockwise=False + ) + + ax = voltage_integral.plot(z=0) + current_integral.plot(z=0, ax=ax) + plt.close() + + # Test off center plotting + ax = voltage_integral.plot(z=2) + current_integral.plot(z=2, ax=ax) + plt.close() + + # Plot + voltage_integral = CustomVoltageIntegral2D( + axis=1, position=0, vertices=[(-1, -1), (0, 0), (1, 1)] + ) + + current_integral = CurrentIntegralAxisAligned( + center=(0, 0, 0), + size=(2, 0, 1), + sign="-", + extrapolate_to_endpoints=False, + snap_contour_to_grid=False, + ) + + ax = voltage_integral.plot(y=0) + current_integral.plot(y=0, ax=ax) + plt.close() + + # Test off center plotting + ax = voltage_integral.plot(y=2) + current_integral.plot(y=2, ax=ax) + plt.close() + + +def test_creation_from_terminal_positions(): + """Test creating an VoltageIntegralAxisAligned using terminal positions.""" + _ = VoltageIntegralAxisAligned.from_terminal_positions( + plus_terminal=2, minus_terminal=1, y=2.2, z=1 + ) diff --git a/tests/test_plugins/test_terminal_component_modeler.py b/tests/test_plugins/test_terminal_component_modeler.py index b23935f5b4..d855f5a322 100644 --- a/tests/test_plugins/test_terminal_component_modeler.py +++ b/tests/test_plugins/test_terminal_component_modeler.py @@ -4,12 +4,15 @@ import pytest import tidy3d as td from tidy3d.exceptions import SetupError, Tidy3dKeyError +from tidy3d.plugins.microwave import CustomCurrentIntegral2D, VoltageIntegralAxisAligned from tidy3d.plugins.smatrix import ( AbstractComponentModeler, CoaxialLumpedPort, LumpedPort, - LumpedPortDataArray, + PortDataArray, TerminalComponentModeler, + TerminalPortDataArray, + WavePort, ) from ..utils import run_emulated @@ -21,9 +24,16 @@ def run_component_modeler(monkeypatch, modeler: TerminalComponentModeler): batch_data = {task_name: run_emulated(sim) for task_name, sim in sim_dict.items()} monkeypatch.setattr(TerminalComponentModeler, "_run_sims", lambda self, path_dir: batch_data) # for the random data, the power wave matrix might be singular, leading to an error - # during inversion, so monkeypatch the inv method so that it operates on a dummy - # identity matrix - monkeypatch.setattr(AbstractComponentModeler, "inv", lambda matrix: np.eye(len(modeler.ports))) + # during inversion, so monkeypatch the inv method so that it operates on a random matrix + monkeypatch.setattr( + AbstractComponentModeler, "inv", lambda matrix: np.random.rand(*matrix.shape) + 1e-4 + ) + # Need to do something similar when computing F + monkeypatch.setattr( + TerminalComponentModeler, + "_compute_F", + lambda matrix: 1.0 / (2.0 * np.sqrt(np.abs(matrix) + 1e-4)), + ) s_matrix = modeler.run(path_dir=modeler.path_dir) return s_matrix @@ -105,10 +115,10 @@ def test_s_to_z_component_modeler(): S22 = ((Z11 + Z0) * (Z22 - Z0) - Z12 * Z21) / deltaZ port_names = ["lumped_port_1", "lumped_port_2"] - freqs = [1e8] + freqs = [1e8, 2e8, 3e8] values = np.array( - [[[S11, S12], [S21, S22]]], + 3 * [[[S11, S12], [S21, S22]]], dtype=complex, ) # Put coords in opposite order to check reordering @@ -118,8 +128,22 @@ def test_s_to_z_component_modeler(): port_in=port_names, ) - s_matrix = LumpedPortDataArray(data=values, coords=coords) - z_matrix = AbstractComponentModeler.s_to_z(s_matrix, reference=Z0) + s_matrix = TerminalPortDataArray(data=values, coords=coords) + z_matrix = TerminalComponentModeler.s_to_z(s_matrix, reference=Z0) + z_matrix_at_f = z_matrix.sel(f=1e8) + assert np.isclose(z_matrix_at_f[0, 0], Z11) + assert np.isclose(z_matrix_at_f[0, 1], Z12) + assert np.isclose(z_matrix_at_f[1, 0], Z21) + assert np.isclose(z_matrix_at_f[1, 1], Z22) + + # test version with different port reference impedances + values = np.full((len(freqs), len(port_names)), Z0) + coords = dict( + f=np.array(freqs), + port=port_names, + ) + z_port_matrix = PortDataArray(data=values, coords=coords) + z_matrix = TerminalComponentModeler.s_to_z(s_matrix, reference=z_port_matrix) z_matrix_at_f = z_matrix.sel(f=1e8) assert np.isclose(z_matrix_at_f[0, 0], Z11) assert np.isclose(z_matrix_at_f[0, 1], Z12) @@ -139,9 +163,9 @@ def test_ab_to_s_component_modeler(): a_values = np.eye(2, 2) a_values = np.reshape(a_values, (1, 2, 2)) b_values = (1 + 1j) * np.random.random((1, 2, 2)) - a_matrix = LumpedPortDataArray(data=a_values, coords=coords) - b_matrix = LumpedPortDataArray(data=b_values, coords=coords) - S_matrix = AbstractComponentModeler.ab_to_s(a_matrix, b_matrix) + a_matrix = TerminalPortDataArray(data=a_values, coords=coords) + b_matrix = TerminalPortDataArray(data=b_values, coords=coords) + S_matrix = TerminalComponentModeler.ab_to_s(a_matrix, b_matrix) assert np.isclose(S_matrix, b_matrix).all() @@ -221,3 +245,139 @@ def test_validate_coaxial_port_diameters(): num_grid_cells=None, impedance=50, ) + + +def test_make_coaxial_component_modeler_with_wave_ports(tmp_path): + """Checks that the terminal component modeler is created successfully with wave ports.""" + _ = make_coaxial_component_modeler( + path_dir=str(tmp_path), port_types=(WavePort, WavePort), auto_grid=True + ) + + +def test_run_coaxial_component_modeler_with_wave_ports(monkeypatch, tmp_path): + """Checks that the terminal component modeler runs with wave ports.""" + modeler = make_coaxial_component_modeler( + path_dir=str(tmp_path), port_types=(WavePort, WavePort), auto_grid=True + ) + s_matrix = run_component_modeler(monkeypatch, modeler) + + shape_one_port = (len(modeler.freqs), len(modeler.ports)) + shape_both_ports = (len(modeler.freqs),) + for port_in in modeler.ports: + for port_out in modeler.ports: + coords_in = dict(port_in=port_in.name) + coords_out = dict(port_out=port_out.name) + + assert np.all( + s_matrix.sel(**coords_in).values.shape == shape_one_port + ), "source index not present in S matrix" + assert np.all( + s_matrix.sel(**coords_in).sel(**coords_out).values.shape == shape_both_ports + ), "monitor index not present in S matrix" + + +def test_run_mixed_component_modeler_with_wave_ports(monkeypatch, tmp_path): + """Checks the terminal component modeler will allow mixed ports.""" + modeler = make_coaxial_component_modeler( + path_dir=str(tmp_path), port_types=(CoaxialLumpedPort, WavePort), auto_grid=True + ) + s_matrix = run_component_modeler(monkeypatch, modeler) + + shape_one_port = (len(modeler.freqs), len(modeler.ports)) + shape_both_ports = (len(modeler.freqs),) + for port_in in modeler.ports: + for port_out in modeler.ports: + coords_in = dict(port_in=port_in.name) + coords_out = dict(port_out=port_out.name) + + assert np.all( + s_matrix.sel(**coords_in).values.shape == shape_one_port + ), "source index not present in S matrix" + assert np.all( + s_matrix.sel(**coords_in).sel(**coords_out).values.shape == shape_both_ports + ), "monitor index not present in S matrix" + + +def test_wave_port_path_integral_validation(): + """Checks that wave port will ensure path integrals are within the bounds of the port.""" + size_port = [2, 2, 0] + center_port = [0, 0, -10] + + voltage_path = VoltageIntegralAxisAligned( + center=(0.5, 0, -10), + size=(1.0, 0, 0), + extrapolate_to_endpoints=True, + snap_path_to_grid=True, + sign="+", + ) + + custom_current_path = CustomCurrentIntegral2D.from_circular_path( + center=center_port, radius=0.5, num_points=21, normal_axis=2, clockwise=False + ) + + mode_spec = td.ModeSpec(num_modes=1, target_neff=1.8) + + _ = WavePort( + center=center_port, + size=size_port, + name="wave_port_1", + mode_spec=mode_spec, + direction="+", + voltage_integral=voltage_path, + current_integral=None, + ) + + _ = WavePort( + center=center_port, + size=size_port, + name="wave_port_1", + mode_spec=mode_spec, + direction="+", + voltage_integral=None, + current_integral=custom_current_path, + ) + + with pytest.raises(pydantic.ValidationError): + _ = WavePort( + center=center_port, + size=size_port, + name="wave_port_1", + mode_spec=mode_spec, + direction="+", + voltage_integral=None, + current_integral=None, + ) + + voltage_path = voltage_path.updated_copy(size=(4, 0, 0)) + with pytest.raises(pydantic.ValidationError): + _ = WavePort( + center=center_port, + size=size_port, + name="wave_port_1", + mode_spec=mode_spec, + direction="+", + voltage_integral=voltage_path, + current_integral=None, + ) + + custom_current_path = CustomCurrentIntegral2D.from_circular_path( + center=center_port, radius=3, num_points=21, normal_axis=2, clockwise=False + ) + with pytest.raises(pydantic.ValidationError): + _ = WavePort( + center=center_port, + size=size_port, + name="wave_port_1", + mode_spec=mode_spec, + direction="+", + voltage_integral=None, + current_integral=custom_current_path, + ) + + +def test_wave_port_to_mode_solver(tmp_path): + """Checks that wave port can be converted to a mode solver.""" + modeler = make_coaxial_component_modeler( + path_dir=str(tmp_path), port_types=(WavePort, WavePort), auto_grid=True + ) + _ = modeler.ports[0].to_mode_solver(modeler.simulation, freqs=[1e9, 2e9, 3e9]) diff --git a/tidy3d/components/data/monitor_data.py b/tidy3d/components/data/monitor_data.py index 5835762a10..771dff874b 100644 --- a/tidy3d/components/data/monitor_data.py +++ b/tidy3d/components/data/monitor_data.py @@ -1389,7 +1389,11 @@ def _isel(self, **isel_kwargs): newly created data.""" update_dict = dict(self._grid_correction_dict, **self.field_components) - update_dict = {key: field.isel(**isel_kwargs) for key, field in update_dict.items()} + update_dict = { + key: field.isel(**isel_kwargs) + for key, field in update_dict.items() + if isinstance(field, DataArray) + } return self._updated(update=update_dict) def _assign_coords(self, **assign_coords_kwargs): diff --git a/tidy3d/components/geometry/base.py b/tidy3d/components/geometry/base.py index ba858a615f..1043ee0cd2 100644 --- a/tidy3d/components/geometry/base.py +++ b/tidy3d/components/geometry/base.py @@ -684,6 +684,29 @@ def parse_xyz_kwargs(**xyz) -> Tuple[Axis, float]: axis = "xyz".index(axis_label) return axis, position + @staticmethod + def parse_two_xyz_kwargs(**xyz) -> List[Tuple[Axis, float]]: + """Turns x,y,z kwargs into indices of axes and the position along each axis. + + Parameters + ---------- + x : float = None + Position in x direction, only two of x,y,z can be specified to define line. + y : float = None + Position in y direction, only two of x,y,z can be specified to define line. + z : float = None + Position in z direction, only two of x,y,z can be specified to define line. + + Returns + ------- + [(int, float), (int, float)] + Index into xyz axis (0,1,2) and position along that axis. + """ + xyz_filtered = {k: v for k, v in xyz.items() if v is not None} + assert len(xyz_filtered) == 2, "exactly two kwarg in [x,y,z] must be specified." + xyz_list = list(xyz_filtered.items()) + return [("xyz".index(axis_label), position) for axis_label, position in xyz_list] + @staticmethod def rotate_points(points: ArrayFloat3D, axis: Coordinate, angle: float) -> ArrayFloat3D: """Rotate a set of points in 3D. diff --git a/tidy3d/components/viz.py b/tidy3d/components/viz.py index b26db20f8f..7631c0f9fa 100644 --- a/tidy3d/components/viz.py +++ b/tidy3d/components/viz.py @@ -73,18 +73,15 @@ def _plot(*args, **kwargs) -> Ax: """ plot parameters """ -class PlotParams(Tidy3dBaseModel): - """Stores plotting parameters / specifications for a given model.""" +class AbstractPlotParams(Tidy3dBaseModel): + """Abstract class for storing plotting parameters. + Corresponds with select properties of ``matplotlib.artist.Artist``. + """ alpha: Any = pd.Field(1.0, title="Opacity") - edgecolor: Any = pd.Field(None, title="Edge Color", alias="ec") - facecolor: Any = pd.Field(None, title="Face Color", alias="fc") - fill: bool = pd.Field(True, title="Is Filled") - hatch: str = pd.Field(None, title="Hatch Style") zorder: float = pd.Field(None, title="Display Order") - linewidth: pd.NonNegativeFloat = pd.Field(1, title="Line Width", alias="lw") - def include_kwargs(self, **kwargs) -> PlotParams: + def include_kwargs(self, **kwargs) -> AbstractPlotParams: """Update the plot params with supplied kwargs.""" update_dict = { key: value @@ -101,6 +98,32 @@ def to_kwargs(self) -> dict: return kwarg_dict +class PathPlotParams(AbstractPlotParams): + """Stores plotting parameters / specifications for a path. + Corresponds with select properties of ``matplotlib.lines.Line2D``. + """ + + color: Any = pd.Field(None, title="Color", alias="c") + linewidth: pd.NonNegativeFloat = pd.Field(2, title="Line Width", alias="lw") + linestyle: str = pd.Field("--", title="Line Style", alias="ls") + marker: Any = pd.Field("o", title="Marker Style") + markeredgecolor: Any = pd.Field(None, title="Marker Edge Color", alias="mec") + markerfacecolor: Any = pd.Field(None, title="Marker Face Color", alias="mfc") + markersize: pd.NonNegativeFloat = pd.Field(10, title="Marker Size", alias="ms") + + +class PlotParams(AbstractPlotParams): + """Stores plotting parameters / specifications for a given model. + Corresponds with select properties of ``matplotlib.patches.Patch``. + """ + + edgecolor: Any = pd.Field(None, title="Edge Color", alias="ec") + facecolor: Any = pd.Field(None, title="Face Color", alias="fc") + fill: bool = pd.Field(True, title="Is Filled") + hatch: str = pd.Field(None, title="Hatch Style") + linewidth: pd.NonNegativeFloat = pd.Field(1, title="Line Width", alias="lw") + + # defaults for different tidy3d objects plot_params_geometry = PlotParams() plot_params_structure = PlotParams() diff --git a/tidy3d/constants.py b/tidy3d/constants.py index 9867738e8b..7c2da4213d 100644 --- a/tidy3d/constants.py +++ b/tidy3d/constants.py @@ -164,7 +164,12 @@ OHM = "ohm" """ -SI unit of resistance.. +SI unit of resistance. +""" + +AMP = "A" +""" +SI unit of electric current. """ THERMAL_CONDUCTIVITY = "W/(um*K)" diff --git a/tidy3d/plugins/microwave/__init__.py b/tidy3d/plugins/microwave/__init__.py index 4c69532254..c9d308e1a7 100644 --- a/tidy3d/plugins/microwave/__init__.py +++ b/tidy3d/plugins/microwave/__init__.py @@ -6,7 +6,7 @@ CustomPathIntegral2D, CustomVoltageIntegral2D, ) -from .impedance_calculator import ImpedanceCalculator +from .impedance_calculator import CurrentIntegralTypes, ImpedanceCalculator, VoltageIntegralTypes from .path_integrals import ( AxisAlignedPathIntegral, CurrentIntegralAxisAligned, @@ -21,6 +21,8 @@ "CurrentIntegralAxisAligned", "CustomVoltageIntegral2D", "CustomCurrentIntegral2D", + "VoltageIntegralTypes", + "CurrentIntegralTypes", "ImpedanceCalculator", "models", "rf_material_library", diff --git a/tidy3d/plugins/microwave/custom_path_integrals.py b/tidy3d/plugins/microwave/custom_path_integrals.py index 14946b7f94..1fca7a80a4 100644 --- a/tidy3d/plugins/microwave/custom_path_integrals.py +++ b/tidy3d/plugins/microwave/custom_path_integrals.py @@ -11,10 +11,26 @@ from ...components.base import cached_property from ...components.data.data_array import FreqDataArray, FreqModeDataArray, TimeDataArray from ...components.data.monitor_data import FieldData, FieldTimeData, ModeSolverData -from ...components.types import ArrayFloat2D, Axis +from ...components.geometry.base import Geometry +from ...components.types import ArrayFloat2D, Ax, Axis, Bound, Coordinate +from ...components.viz import add_ax_if_none from ...constants import MICROMETER, fp_eps from ...exceptions import DataError, SetupError -from .path_integrals import AbstractAxesRH, IntegralResultTypes, MonitorDataTypes +from .path_integrals import ( + AbstractAxesRH, + CurrentIntegralAxisAligned, + IntegralResultTypes, + MonitorDataTypes, + VoltageIntegralAxisAligned, + _check_em_field_supported, +) +from .viz import ( + ARROW_CURRENT, + plot_params_current_path, + plot_params_voltage_minus, + plot_params_voltage_path, + plot_params_voltage_plus, +) FieldParameter = Literal["E", "H"] @@ -72,8 +88,6 @@ def compute_integral( :class:`.IntegralResultTypes` Result of integral over remaining dimensions (frequency, time, mode indices). """ - if not isinstance(em_field, (FieldData, FieldTimeData, ModeSolverData)): - raise DataError("'em_field' type not supported.") (dim1, dim2, dim3) = self.local_dims @@ -138,8 +152,61 @@ def _compute_dl_component(coord_array: xr.DataArray, closed_contour=False) -> np dl[0] = dl[-1] = grad_end[1] return dl + @classmethod + def from_circular_path( + cls, center: Coordinate, radius: float, num_points: int, normal_axis: Axis, clockwise: bool + ) -> CustomPathIntegral2D: + """Creates a ``CustomPathIntegral2D`` from a circular path given a desired number of points + along the perimeter. + + Parameters + ---------- + center : Coordinate + The center of the circle. + radius : float + The radius of the circle. + num_points : int + THe number of equidistant points to use along the perimeter of the circle. + normal_axis : Axis + The axis normal to the defined circle. + clockwise : bool + When ``True``, the points will be ordered clockwise with respect to the positive + direction of the ``normal_axis``. + + Returns + ------- + :class:`.CustomPathIntegral2D` + A path integral defined on a circular path. + """ + + def generate_circle_coordinates(radius: float, num_points: int, clockwise: bool): + """Helper for generating x,y vertices around a circle in the local coordinate frame.""" + sign = 1.0 + if clockwise: + sign = -1.0 + angles = np.linspace(0, sign * 2 * np.pi, num_points, endpoint=True) + xt = radius * np.cos(angles) + yt = radius * np.sin(angles) + return (xt, yt) + + # Get transverse axes + normal_center, trans_center = Geometry.pop_axis(center, normal_axis) + + # These x,y coordinates in the local coordinate frame + if normal_axis == 1: + # Handle special case when y is the axis that is popped + clockwise = not clockwise + xt, yt = generate_circle_coordinates(radius, num_points, clockwise) + xt += trans_center[0] + yt += trans_center[1] + circle_vertices = np.column_stack((xt, yt)) + # Close the contour exactly + circle_vertices[-1, :] = circle_vertices[0, :] + return cls(axis=normal_axis, position=normal_center, vertices=circle_vertices) + @cached_property def is_closed_contour(self) -> bool: + """Returns ``true`` when the first vertex equals the last vertex.""" return np.isclose( self.vertices[0, :], self.vertices[-1, :], @@ -163,6 +230,15 @@ def _correct_shape(cls, val): ) return val + @cached_property + def bounds(self) -> Bound: + """Helper to get the geometric bounding box of the path integral.""" + path_min = np.amin(self.vertices, axis=0) + path_max = np.amax(self.vertices, axis=0) + min_bound = Geometry.unpop_axis(self.position, path_min, self.axis) + max_bound = Geometry.unpop_axis(self.position, path_max, self.axis) + return (min_bound, max_bound) + class CustomVoltageIntegral2D(CustomPathIntegral2D): """Class for computing the voltage between two points defined by a custom path. @@ -189,9 +265,60 @@ def compute_voltage(self, em_field: MonitorDataTypes) -> IntegralResultTypes: :class:`.IntegralResultTypes` Result of voltage computation over remaining dimensions (frequency, time, mode indices). """ + _check_em_field_supported(em_field=em_field) voltage = -1.0 * self.compute_integral(field="E", em_field=em_field) + voltage = VoltageIntegralAxisAligned._set_data_array_attributes(voltage) return voltage + @add_ax_if_none + def plot( + self, + x: float = None, + y: float = None, + z: float = None, + ax: Ax = None, + **path_kwargs, + ) -> Ax: + """Plot path integral at single (x,y,z) coordinate. + + Parameters + ---------- + x : float = None + Position of plane in x direction, only one of x,y,z can be specified to define plane. + y : float = None + Position of plane in y direction, only one of x,y,z can be specified to define plane. + z : float = None + Position of plane in z direction, only one of x,y,z can be specified to define plane. + ax : matplotlib.axes._subplots.Axes = None + Matplotlib axes to plot on, if not specified, one is created. + **path_kwargs + Optional keyword arguments passed to the matplotlib plotting of the line. + For details on accepted values, refer to + `Matplotlib's documentation `_. + + Returns + ------- + matplotlib.axes._subplots.Axes + The supplied or created matplotlib axes. + """ + axis, position = Geometry.parse_xyz_kwargs(x=x, y=y, z=z) + if axis != self.main_axis or not np.isclose(position, self.position, rtol=fp_eps): + return ax + + plot_params = plot_params_voltage_path.include_kwargs(**path_kwargs) + plot_kwargs = plot_params.to_kwargs() + xs = self.vertices[:, 0] + ys = self.vertices[:, 1] + ax.plot(xs, ys, markevery=[0, -1], **plot_kwargs) + + # Plot special end points + end_kwargs = plot_params_voltage_plus.include_kwargs(**path_kwargs).to_kwargs() + start_kwargs = plot_params_voltage_minus.include_kwargs(**path_kwargs).to_kwargs() + ax.plot(xs[0], ys[0], **start_kwargs) + ax.plot(xs[-1], ys[-1], **end_kwargs) + + return ax + class CustomCurrentIntegral2D(CustomPathIntegral2D): """Class for computing conduction current via Ampère's circuital law on a custom path. @@ -211,5 +338,57 @@ def compute_current(self, em_field: MonitorDataTypes) -> IntegralResultTypes: :class:`.IntegralResultTypes` Result of current computation over remaining dimensions (frequency, time, mode indices). """ + _check_em_field_supported(em_field=em_field) current = self.compute_integral(field="H", em_field=em_field) + current = CurrentIntegralAxisAligned._set_data_array_attributes(current) return current + + @add_ax_if_none + def plot( + self, + x: float = None, + y: float = None, + z: float = None, + ax: Ax = None, + **path_kwargs, + ) -> Ax: + """Plot path integral at single (x,y,z) coordinate. + + Parameters + ---------- + x : float = None + Position of plane in x direction, only one of x,y,z can be specified to define plane. + y : float = None + Position of plane in y direction, only one of x,y,z can be specified to define plane. + z : float = None + Position of plane in z direction, only one of x,y,z can be specified to define plane. + ax : matplotlib.axes._subplots.Axes = None + Matplotlib axes to plot on, if not specified, one is created. + **path_kwargs + Optional keyword arguments passed to the matplotlib plotting of the line. + For details on accepted values, refer to + `Matplotlib's documentation `_. + + Returns + ------- + matplotlib.axes._subplots.Axes + The supplied or created matplotlib axes. + """ + axis, position = Geometry.parse_xyz_kwargs(x=x, y=y, z=z) + if axis != self.main_axis or not np.isclose(position, self.position, rtol=fp_eps): + return ax + + plot_params = plot_params_current_path.include_kwargs(**path_kwargs) + plot_kwargs = plot_params.to_kwargs() + xs = self.vertices[:, 0] + ys = self.vertices[:, 1] + ax.plot(xs, ys, **plot_kwargs) + + # Add arrow at start of contour + ax.annotate( + "", + xytext=(xs[0], ys[0]), + xy=(xs[1], ys[1]), + arrowprops=ARROW_CURRENT, + ) + return ax diff --git a/tidy3d/plugins/microwave/impedance_calculator.py b/tidy3d/plugins/microwave/impedance_calculator.py index 430108a2dc..0d38874872 100644 --- a/tidy3d/plugins/microwave/impedance_calculator.py +++ b/tidy3d/plugins/microwave/impedance_calculator.py @@ -8,14 +8,16 @@ import pydantic.v1 as pd from ...components.base import Tidy3dBaseModel -from ...components.data.monitor_data import FieldData, FieldTimeData, ModeSolverData -from ...exceptions import DataError, ValidationError +from ...components.data.monitor_data import FieldTimeData +from ...constants import OHM +from ...exceptions import ValidationError from .custom_path_integrals import CustomCurrentIntegral2D, CustomVoltageIntegral2D from .path_integrals import ( CurrentIntegralAxisAligned, IntegralResultTypes, MonitorDataTypes, VoltageIntegralAxisAligned, + _check_em_field_supported, ) VoltageIntegralTypes = Union[VoltageIntegralAxisAligned, CustomVoltageIntegral2D] @@ -53,9 +55,7 @@ def compute_impedance(self, em_field: MonitorDataTypes) -> IntegralResultTypes: :class:`.IntegralResultTypes` Result of impedance computation over remaining dimensions (frequency, time, mode indices). """ - - if not isinstance(em_field, (FieldData, FieldTimeData, ModeSolverData)): - raise DataError("'em_field' type not supported by impedance calculator.") + _check_em_field_supported(em_field=em_field) # If both voltage and current integrals have been defined then impedance is computed directly if self.voltage_integral: @@ -63,11 +63,12 @@ def compute_impedance(self, em_field: MonitorDataTypes) -> IntegralResultTypes: 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 only one of the integrals has been provided, then the computation falls back to using + # total power (flux) with Ohm's law to compute the missing quantity. 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): @@ -82,14 +83,21 @@ def compute_impedance(self, em_field: MonitorDataTypes) -> IntegralResultTypes: current = np.conj(2 * flux / voltage) impedance = voltage / current + impedance = ImpedanceCalculator._set_data_array_attributes(impedance) 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.""" + are not provided.""" if not values.get("voltage_integral") and not val: raise ValidationError( - "Atleast one of 'voltage_integral' or 'current_integral' must be provided." + "At least one of 'voltage_integral' or 'current_integral' must be provided." ) return val + + @staticmethod + def _set_data_array_attributes(data_array: IntegralResultTypes) -> IntegralResultTypes: + """Helper to set additional metadata for ``IntegralResultTypes``.""" + data_array.name = "Z0" + return data_array.assign_attrs(units=OHM, long_name="characteristic impedance") diff --git a/tidy3d/plugins/microwave/path_integrals.py b/tidy3d/plugins/microwave/path_integrals.py index bfd05acc92..76e769ed7c 100644 --- a/tidy3d/plugins/microwave/path_integrals.py +++ b/tidy3d/plugins/microwave/path_integrals.py @@ -3,10 +3,11 @@ from __future__ import annotations from abc import ABC, abstractmethod -from typing import Union +from typing import Any, Union import numpy as np import pydantic.v1 as pd +import shapely as shapely from ...components.base import Tidy3dBaseModel, cached_property from ...components.data.data_array import ( @@ -18,16 +19,34 @@ TimeDataArray, ) from ...components.data.monitor_data import FieldData, FieldTimeData, ModeSolverData -from ...components.geometry.base import Box -from ...components.types import Axis, Direction +from ...components.geometry.base import Box, Geometry +from ...components.types import Ax, Axis, Coordinate2D, Direction from ...components.validators import assert_line, assert_plane +from ...components.viz import add_ax_if_none +from ...constants import AMP, VOLT, fp_eps from ...exceptions import DataError, Tidy3dError +from .viz import ( + ARROW_CURRENT, + plot_params_current_path, + plot_params_voltage_minus, + plot_params_voltage_path, + plot_params_voltage_plus, +) MonitorDataTypes = Union[FieldData, FieldTimeData, ModeSolverData] EMScalarFieldType = Union[ScalarFieldDataArray, ScalarFieldTimeDataArray, ScalarModeFieldDataArray] IntegralResultTypes = Union[FreqDataArray, FreqModeDataArray, TimeDataArray] +def _check_em_field_supported(em_field: Any): + """Function for validating correct data arrays.""" + if not isinstance(em_field, (FieldData, FieldTimeData, ModeSolverData)): + raise DataError( + "'em_field' type not supported. Supported types are " + "'FieldData', 'FieldTimeData', 'ModeSolverData'." + ) + + class AbstractAxesRH(Tidy3dBaseModel, ABC): """Represents an axis-aligned right-handed coordinate system with one axis preferred. Typically `main_axis` would refer to the normal axis of a plane. @@ -171,6 +190,17 @@ def main_axis(self) -> Axis: return index raise Tidy3dError("Failed to identify axis.") + def _vertices_2D(self, axis: Axis) -> tuple[Coordinate2D, Coordinate2D]: + """Returns the two vertices of this path in the plane defined by ``axis``.""" + min = self.bounds[0] + max = self.bounds[1] + _, min = Box.pop_axis(min, axis) + _, max = Box.pop_axis(max, axis) + + u = [min[0], max[0]] + v = [min[1], max[1]] + return (u, v) + class VoltageIntegralAxisAligned(AxisAlignedPathIntegral): """Class for computing the voltage between two points defined by an axis-aligned line.""" @@ -183,22 +213,138 @@ class VoltageIntegralAxisAligned(AxisAlignedPathIntegral): 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.") + _check_em_field_supported(em_field=em_field) 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 + + voltage = VoltageIntegralAxisAligned._set_data_array_attributes(voltage) # Return data array of voltage while keeping coordinates of frequency|time|mode index return voltage + @staticmethod + def _set_data_array_attributes(data_array: IntegralResultTypes) -> IntegralResultTypes: + """Add explanatory attributes to the data array.""" + data_array.name = "V" + return data_array.assign_attrs(units=VOLT, long_name="voltage") + + @staticmethod + def from_terminal_positions( + plus_terminal: float, + minus_terminal: float, + x: float = None, + y: float = None, + z: float = None, + extrapolate_to_endpoints: bool = True, + snap_path_to_grid: bool = True, + ) -> VoltageIntegralAxisAligned: + """Helper to create a :class:`VoltageIntegralAxisAligned` from two coordinates that + define a line and two positions indicating the endpoints of the path integral. + + Parameters + ---------- + plus_terminal : float + Position along the voltage axis of the positive terminal. + minus_terminal : float + Position along the voltage axis of the negative terminal. + x : float = None + Position in x direction, only two of x,y,z can be specified to define line. + y : float = None + Position in y direction, only two of x,y,z can be specified to define line. + z : float = None + Position in z direction, only two of x,y,z can be specified to define line. + extrapolate_to_endpoints: bool = True + Passed directly to :class:`VoltageIntegralAxisAligned` + snap_path_to_grid: bool = True + Passed directly to :class:`VoltageIntegralAxisAligned` + + Returns + ------- + VoltageIntegralAxisAligned + The created path integral for computing voltage between the two terminals. + """ + axis_positions = Geometry.parse_two_xyz_kwargs(x=x, y=y, z=z) + # Calculate center and size of the future box + midpoint = (plus_terminal + minus_terminal) / 2 + length = np.abs(plus_terminal - minus_terminal) + center = [midpoint, midpoint, midpoint] + size = [length, length, length] + for axis, position in axis_positions: + size[axis] = 0 + center[axis] = position + + direction = "+" + if plus_terminal < minus_terminal: + direction = "-" + + return VoltageIntegralAxisAligned( + center=center, + size=size, + extrapolate_to_endpoints=extrapolate_to_endpoints, + snap_path_to_grid=snap_path_to_grid, + sign=direction, + ) + + @add_ax_if_none + def plot( + self, + x: float = None, + y: float = None, + z: float = None, + ax: Ax = None, + **path_kwargs, + ) -> Ax: + """Plot path integral at single (x,y,z) coordinate. + + Parameters + ---------- + x : float = None + Position of plane in x direction, only one of x,y,z can be specified to define plane. + y : float = None + Position of plane in y direction, only one of x,y,z can be specified to define plane. + z : float = None + Position of plane in z direction, only one of x,y,z can be specified to define plane. + ax : matplotlib.axes._subplots.Axes = None + Matplotlib axes to plot on, if not specified, one is created. + **path_kwargs + Optional keyword arguments passed to the matplotlib plotting of the line. + For details on accepted values, refer to + `Matplotlib's documentation `_. + + Returns + ------- + matplotlib.axes._subplots.Axes + The supplied or created matplotlib axes. + """ + axis, position = self.parse_xyz_kwargs(x=x, y=y, z=z) + if axis == self.main_axis or not np.isclose(position, self.center[axis], rtol=fp_eps): + return ax + + (xs, ys) = self._vertices_2D(axis) + # Plot the path + plot_params = plot_params_voltage_path.include_kwargs(**path_kwargs) + plot_kwargs = plot_params.to_kwargs() + ax.plot(xs, ys, markevery=[0, -1], **plot_kwargs) + + # Plot special end points + end_kwargs = plot_params_voltage_plus.include_kwargs(**path_kwargs).to_kwargs() + start_kwargs = plot_params_voltage_minus.include_kwargs(**path_kwargs).to_kwargs() + + if self.sign == "-": + start_kwargs, end_kwargs = end_kwargs, start_kwargs + + ax.plot(xs[0], ys[0], **start_kwargs) + ax.plot(xs[1], ys[1], **end_kwargs) + return ax + class CurrentIntegralAxisAligned(AbstractAxesRH, Box): """Class for computing conduction current via Ampère's circuital law on an axis-aligned loop.""" @@ -214,19 +360,18 @@ class CurrentIntegralAxisAligned(AbstractAxesRH, Box): extrapolate_to_endpoints: bool = pd.Field( False, title="Extrapolate to Endpoints", - description="This parameter is passed to ``AxisAlignedPathIntegral`` objects when computing the contour integral.", + description="This parameter is passed to :class:`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.", + description="This parameter is passed to :class:`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.") + _check_em_field_supported(em_field=em_field) ax1 = self.remaining_axes[0] ax2 = self.remaining_axes[1] h_component = "xyz"[ax1] @@ -253,7 +398,7 @@ def compute_current(self, em_field: MonitorDataTypes) -> IntegralResultTypes: if self.sign == "-": current *= -1 - + current = CurrentIntegralAxisAligned._set_data_array_attributes(current) return current @cached_property @@ -264,13 +409,17 @@ def main_axis(self) -> Axis: return index raise Tidy3dError("Failed to identify axis.") - def _to_path_integrals(self, h_horizontal, h_vertical) -> tuple[AxisAlignedPathIntegral, ...]: + def _to_path_integrals( + self, h_horizontal=None, h_vertical=None + ) -> 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: + horizontal_passed = h_horizontal is not None + vertical_passed = h_vertical is not None + if self.snap_contour_to_grid and horizontal_passed and vertical_passed: (coord1, coord2) = self.remaining_dims # Locations where horizontal paths will be snapped @@ -339,3 +488,87 @@ def _to_path_integrals(self, h_horizontal, h_vertical) -> tuple[AxisAlignedPathI ) return (bottom, right, top, left) + + @staticmethod + def _set_data_array_attributes(data_array: IntegralResultTypes) -> IntegralResultTypes: + """Add explanatory attributes to the data array.""" + data_array.name = "I" + return data_array.assign_attrs(units=AMP, long_name="current") + + @add_ax_if_none + def plot( + self, + x: float = None, + y: float = None, + z: float = None, + ax: Ax = None, + **path_kwargs, + ) -> Ax: + """Plot path integral at single (x,y,z) coordinate. + + Parameters + ---------- + x : float = None + Position of plane in x direction, only one of x,y,z can be specified to define plane. + y : float = None + Position of plane in y direction, only one of x,y,z can be specified to define plane. + z : float = None + Position of plane in z direction, only one of x,y,z can be specified to define plane. + ax : matplotlib.axes._subplots.Axes = None + Matplotlib axes to plot on, if not specified, one is created. + **path_kwargs + Optional keyword arguments passed to the matplotlib plotting of the line. + For details on accepted values, refer to + `Matplotlib's documentation `_. + + Returns + ------- + matplotlib.axes._subplots.Axes + The supplied or created matplotlib axes. + """ + axis, position = self.parse_xyz_kwargs(x=x, y=y, z=z) + if axis != self.main_axis or not np.isclose(position, self.center[axis], rtol=fp_eps): + return ax + + plot_params = plot_params_current_path.include_kwargs(**path_kwargs) + plot_kwargs = plot_params.to_kwargs() + path_integrals = self._to_path_integrals() + # Plot the path + for path in path_integrals: + (xs, ys) = path._vertices_2D(axis) + ax.plot(xs, ys, **plot_kwargs) + + (ax1, ax2) = self.remaining_axes + + # Add arrow to bottom path, unless right path is longer + arrow_path = path_integrals[0] + if self.size[ax2] > self.size[ax1]: + arrow_path = path_integrals[1] + + (xs, ys) = arrow_path._vertices_2D(axis) + X = (xs[0] + xs[1]) / 2 + Y = (ys[0] + ys[1]) / 2 + center = np.array([X, Y]) + dx = xs[1] - xs[0] + dy = ys[1] - ys[0] + direction = np.array([dx, dy]) + segment_length = np.linalg.norm(direction) + unit_dir = direction / segment_length + + # Change direction of arrow depending on sign of current definition + if self.sign == "-": + unit_dir *= -1.0 + # Change direction of arrow when the "y" axis is dropped, + # since the plotted coordinate system will be left-handed (x, z) + if self.main_axis == 1: + unit_dir *= -1.0 + + start = center - unit_dir * segment_length + end = center + ax.annotate( + "", + xytext=(start[0], start[1]), + xy=(end[0], end[1]), + arrowprops=ARROW_CURRENT, + ) + return ax diff --git a/tidy3d/plugins/microwave/viz.py b/tidy3d/plugins/microwave/viz.py new file mode 100644 index 0000000000..04c61b812e --- /dev/null +++ b/tidy3d/plugins/microwave/viz.py @@ -0,0 +1,54 @@ +"""Utilities for plotting microwave components""" + +from numpy import inf + +from ...components.viz import PathPlotParams + +""" Constants """ +VOLTAGE_COLOR = "red" +CURRENT_COLOR = "blue" +PATH_LINEWIDTH = 2 +ARROW_CURRENT = dict( + arrowstyle="-|>", + mutation_scale=32, + linestyle="", + lw=PATH_LINEWIDTH, + color=CURRENT_COLOR, +) + +plot_params_voltage_path = PathPlotParams( + alpha=1.0, + zorder=inf, + color=VOLTAGE_COLOR, + linestyle="--", + linewidth=PATH_LINEWIDTH, + marker="o", + markersize=10, + markeredgecolor=VOLTAGE_COLOR, + markerfacecolor="white", +) + +plot_params_voltage_plus = PathPlotParams( + alpha=1.0, + zorder=inf, + color=VOLTAGE_COLOR, + marker="+", + markersize=6, +) + +plot_params_voltage_minus = PathPlotParams( + alpha=1.0, + zorder=inf, + color=VOLTAGE_COLOR, + marker="_", + markersize=6, +) + +plot_params_current_path = PathPlotParams( + alpha=1.0, + zorder=inf, + color=CURRENT_COLOR, + linestyle="--", + linewidth=PATH_LINEWIDTH, + marker="", +) diff --git a/tidy3d/plugins/smatrix/__init__.py b/tidy3d/plugins/smatrix/__init__.py index 75d02b8fd1..e88ac42c54 100644 --- a/tidy3d/plugins/smatrix/__init__.py +++ b/tidy3d/plugins/smatrix/__init__.py @@ -1,10 +1,15 @@ """Imports from scattering matrix plugin.""" from .component_modelers.modal import AbstractComponentModeler, ComponentModeler, ModalPortDataArray -from .component_modelers.terminal import LumpedPortDataArray, TerminalComponentModeler +from .component_modelers.terminal import ( + PortDataArray, + TerminalComponentModeler, + TerminalPortDataArray, +) from .ports.coaxial_lumped import CoaxialLumpedPort from .ports.modal import Port from .ports.rectangular_lumped import LumpedPort +from .ports.wave import WavePort __all__ = [ "AbstractComponentModeler", @@ -14,5 +19,7 @@ "TerminalComponentModeler", "CoaxialLumpedPort", "LumpedPort", - "LumpedPortDataArray", + "WavePort", + "TerminalPortDataArray", + "PortDataArray", ] diff --git a/tidy3d/plugins/smatrix/component_modelers/base.py b/tidy3d/plugins/smatrix/component_modelers/base.py index 295ab4b050..e15134418a 100644 --- a/tidy3d/plugins/smatrix/component_modelers/base.py +++ b/tidy3d/plugins/smatrix/component_modelers/base.py @@ -13,18 +13,23 @@ from ....components.base import Tidy3dBaseModel, cached_property from ....components.data.data_array import DataArray from ....components.simulation import Simulation -from ....components.types import Complex, FreqArray +from ....components.types import FreqArray from ....constants import HERTZ from ....exceptions import SetupError, Tidy3dKeyError from ....log import log from ....web.api.container import Batch, BatchData +from ..ports.coaxial_lumped import CoaxialLumpedPort from ..ports.modal import Port from ..ports.rectangular_lumped import LumpedPort +from ..ports.wave import WavePort # fwidth of gaussian pulse in units of central frequency FWIDTH_FRAC = 1.0 / 10 DEFAULT_DATA_DIR = "." +LumpedPortType = Union[LumpedPort, CoaxialLumpedPort] +TerminalPortType = Union[LumpedPortType, WavePort] + class AbstractComponentModeler(ABC, Tidy3dBaseModel): """Tool for modeling devices and computing port parameters.""" @@ -35,7 +40,7 @@ class AbstractComponentModeler(ABC, Tidy3dBaseModel): description="Simulation describing the device without any sources present.", ) - ports: Tuple[Union[Port, LumpedPort], ...] = pd.Field( + ports: Tuple[Union[Port, TerminalPortType], ...] = pd.Field( (), title="Ports", description="Collection of ports describing the scattering matrix elements. " @@ -214,43 +219,46 @@ def load(self, path_dir: str = DEFAULT_DATA_DIR) -> DataArray: batch_data = BatchData.load(path_dir=path_dir) return self._construct_smatrix(batch_data=batch_data) - @staticmethod - def s_to_z(s_matrix: DataArray, reference: Complex) -> DataArray: - """Get the impedance matrix given the scattering matrix and a reference impedance.""" - - # move the input and output port dimensions to the end, for ease of matrix operations - dims = list(s_matrix.dims) - dims.append(dims.pop(dims.index("port_out"))) - dims.append(dims.pop(dims.index("port_in"))) - z_matrix = s_matrix.copy(deep=True).transpose(*dims) - s_vals = z_matrix.values - - eye = np.eye(len(s_matrix.port_out.values), len(s_matrix.port_in.values)) - z_vals = np.matmul(AbstractComponentModeler.inv(eye - s_vals), (eye + s_vals)) * reference - - z_matrix.data = z_vals - return z_matrix - - @staticmethod - def ab_to_s(a_matrix: DataArray, b_matrix: DataArray) -> DataArray: - """Get the scattering matrix given the power wave matrices.""" - - # move the input and output port dimensions to the end, for ease of matrix operations - assert a_matrix.dims == b_matrix.dims - dims = list(a_matrix.dims) - dims.append(dims.pop(dims.index("port_out"))) - dims.append(dims.pop(dims.index("port_in"))) - - s_matrix = a_matrix.copy(deep=True).transpose(*dims) - a_vals = s_matrix.copy(deep=True).transpose(*dims).values - b_vals = b_matrix.copy(deep=True).transpose(*dims).values - - s_vals = np.matmul(b_vals, AbstractComponentModeler.inv(a_vals)) - - s_matrix.data = s_vals - return s_matrix - @staticmethod def inv(matrix: DataArray): """Helper to invert a port matrix.""" return np.linalg.inv(matrix) + + def _shift_value_signed(self, port: Union[Port, WavePort]) -> float: + """How far (signed) to shift the source from the monitor.""" + + # get the grid boundaries and sizes along port normal from the simulation + normal_axis = port.size.index(0.0) + grid = self.simulation.grid + grid_boundaries = grid.boundaries.to_list[normal_axis] + grid_centers = grid.centers.to_list[normal_axis] + + # get the index of the grid cell where the port lies + port_position = port.center[normal_axis] + port_pos_gt_grid_bounds = np.argwhere(port_position > grid_boundaries) + + # no port index can be determined + if len(port_pos_gt_grid_bounds) == 0: + raise SetupError(f"Port position '{port_position}' outside of simulation bounds.") + port_index = port_pos_gt_grid_bounds[-1] + + # shift the port to the left + if port.direction == "+": + shifted_index = port_index - 2 + if shifted_index < 0: + raise SetupError( + f"Port {port.name} normal is too close to boundary " + f"on -{'xyz'[normal_axis]} side." + ) + + # shift the port to the right + else: + shifted_index = port_index + 2 + if shifted_index >= len(grid_centers): + raise SetupError( + f"Port {port.name} normal is too close to boundary " + f"on +{'xyz'[normal_axis]} side." + ) + + new_pos = grid_centers[shifted_index] + return new_pos - port_position diff --git a/tidy3d/plugins/smatrix/component_modelers/modal.py b/tidy3d/plugins/smatrix/component_modelers/modal.py index 6a03f4e348..18e6c3ee27 100644 --- a/tidy3d/plugins/smatrix/component_modelers/modal.py +++ b/tidy3d/plugins/smatrix/component_modelers/modal.py @@ -194,45 +194,6 @@ def to_source(self, port: Port, mode_index: int) -> List[ModeSource]: name=port.name, ) - def _shift_value_signed(self, port: Port) -> float: - """How far (signed) to shift the source from the monitor.""" - - # get the grid boundaries and sizes along port normal from the simulation - normal_axis = port.size.index(0.0) - grid = self.simulation.grid - grid_boundaries = grid.boundaries.to_list[normal_axis] - grid_centers = grid.centers.to_list[normal_axis] - - # get the index of the grid cell where the port lies - port_position = port.center[normal_axis] - port_pos_gt_grid_bounds = np.argwhere(port_position > grid_boundaries) - - # no port index can be determined - if len(port_pos_gt_grid_bounds) == 0: - raise SetupError(f"Port position '{port_position}' outside of simulation bounds.") - port_index = port_pos_gt_grid_bounds[-1] - - # shift the port to the left - if port.direction == "+": - shifted_index = port_index - 2 - if shifted_index < 0: - raise SetupError( - f"Port {port.name} normal is too close to boundary " - f"on -{'xyz'[normal_axis]} side." - ) - - # shift the port to the right - else: - shifted_index = port_index + 2 - if shifted_index >= len(grid_centers): - raise SetupError( - f"Port {port.name} normal is too close to boundary " - f"on +{'xyz'[normal_axis]} side." - ) - - new_pos = grid_centers[shifted_index] - return new_pos - port_position - def shift_port(self, port: Port) -> Port: """Generate a new port shifted by the shift amount in normal direction.""" diff --git a/tidy3d/plugins/smatrix/component_modelers/terminal.py b/tidy3d/plugins/smatrix/component_modelers/terminal.py index 8b074ec9f8..074c9bdacc 100644 --- a/tidy3d/plugins/smatrix/component_modelers/terminal.py +++ b/tidy3d/plugins/smatrix/component_modelers/terminal.py @@ -8,30 +8,46 @@ import pydantic.v1 as pd from ....components.base import cached_property +from ....components.data.data_array import DataArray from ....components.data.sim_data import SimulationData from ....components.geometry.utils_2d import snap_coordinate_to_grid from ....components.simulation import Simulation from ....components.source import GaussianPulse from ....components.types import Ax from ....components.viz import add_ax_if_none, equal_aspect -from ....constants import C_0, fp_eps +from ....constants import C_0, OHM from ....exceptions import ValidationError from ....web.api.container import BatchData -from ..ports.base_lumped import LumpedPortDataArray -from ..ports.coaxial_lumped import CoaxialLumpedPort -from ..ports.rectangular_lumped import LumpedPort -from .base import FWIDTH_FRAC, AbstractComponentModeler +from ..ports.base_lumped import AbstractLumpedPort +from ..ports.base_terminal import AbstractTerminalPort, TerminalPortDataArray +from ..ports.wave import WavePort +from .base import FWIDTH_FRAC, AbstractComponentModeler, TerminalPortType + + +class PortDataArray(DataArray): + """Array of values over dimensions of frequency and port name. + + Example + ------- + >>> f = [2e9, 3e9, 4e9] + >>> ports = ["port1", "port2"] + >>> coords = dict(f=f, port=ports) + >>> fd = PortDataArray((1+1j) * np.random.random((3, 2)), coords=coords) + """ + + __slots__ = () + _dims = ("f", "port") class TerminalComponentModeler(AbstractComponentModeler): """Tool for modeling two-terminal multiport devices and computing port parameters - with lumped ports.""" + with lumped and wave ports.""" - ports: Tuple[Union[LumpedPort, CoaxialLumpedPort], ...] = pd.Field( + ports: Tuple[TerminalPortType, ...] = pd.Field( (), - title="Lumped ports", - description="Collection of lumped ports associated with the network. " - "For each port, one simulation will be run with a lumped port source.", + title="Terminal Ports", + description="Collection of lumped and wave ports associated with the network. " + "For each port, one simulation will be run with a source that is associated with the port.", ) @equal_aspect @@ -43,9 +59,7 @@ def plot_sim( plot_sources = [] for port_source in self.ports: - source_0 = port_source.to_source( - self._source_time, snap_center=None, grid=self.simulation.grid - ) + source_0 = port_source.to_source(self._source_time) plot_sources.append(source_0) sim_plot = self.simulation.copy(update=dict(sources=plot_sources)) return sim_plot.plot(x=x, y=y, z=z, ax=ax, **kwargs) @@ -59,9 +73,7 @@ def plot_sim_eps( plot_sources = [] for port_source in self.ports: - source_0 = port_source.to_source( - self._source_time, snap_center=None, grid=self.simulation.grid - ) + source_0 = port_source.to_source(self._source_time) plot_sources.append(source_0) sim_plot = self.simulation.copy(update=dict(sources=plot_sources)) return sim_plot.plot_eps(x=x, y=y, z=z, ax=ax, **kwargs) @@ -72,14 +84,7 @@ def sim_dict(self) -> Dict[str, Simulation]: sim_dict = {} - port_voltage_monitors = [ - port.to_voltage_monitor(self.freqs, snap_center=None) for port in self.ports - ] - port_current_monitors = [ - port.to_current_monitor(self.freqs, snap_center=None) for port in self.ports - ] - lumped_resistors = [port.to_load(snap_center=None) for port in self.ports] - + lumped_resistors = [port.to_load() for port in self._lumped_ports] # Create a mesh override for each port in case refinement is needed. # The port is a flat surface, but when computing the port current, # we'll eventually integrate the magnetic field just above and below @@ -92,12 +97,10 @@ def sim_dict(self) -> Dict[str, Simulation]: # override regions are defined, one above and one below the analytical # port region. mesh_overrides = [] - for port, lumped_resistor in zip(self.ports, lumped_resistors): + for port, lumped_resistor in zip(self._lumped_ports, lumped_resistors): if port.num_grid_cells: mesh_overrides.extend(lumped_resistor.to_mesh_overrides()) - new_mnts = list(self.simulation.monitors) + port_voltage_monitors + port_current_monitors - # also, use the highest frequency in the simulation to define the grid, rather than the # source's central frequency, to ensure an accurate solution over the entire range grid_spec = self.simulation.grid_spec.copy( @@ -108,67 +111,63 @@ def sim_dict(self) -> Dict[str, Simulation]: } ) - # Checking if snapping is required needs the simulation to be created, because the - # elements may impact the final grid discretization - snap_and_recreate = False - snap_centers = [] - for port in self.ports: - update_dict = dict( - monitors=new_mnts, - lumped_elements=lumped_resistors, - grid_spec=grid_spec, - ) - - sim_copy = self.simulation.copy(update=update_dict) - port_source = port.to_source(self._source_time, snap_center=None, grid=sim_copy.grid) - sim_copy = sim_copy.updated_copy(sources=[port_source]) - task_name = self._task_name(port=port) - sim_dict[task_name] = sim_copy - - # Check if snapping to grid is needed + # Make an initial simulation with new grid_spec to determine where LumpedPorts are snapped + sim_wo_source = self.simulation.updated_copy(grid_spec=grid_spec) + snap_centers = dict() + for port in self._lumped_ports: port_center_on_axis = port.center[port.injection_axis] new_port_center = snap_coordinate_to_grid( - sim_copy.grid, port_center_on_axis, port.injection_axis + sim_wo_source.grid, port_center_on_axis, port.injection_axis ) - snap_centers.append(new_port_center) - if not np.isclose(port_center_on_axis, new_port_center, fp_eps, fp_eps): - snap_and_recreate = True - - # Check if snapping was needed and if it was recreate the simulations - if snap_and_recreate: - sim_dict.clear() - port_voltage_monitors = [ - port.to_voltage_monitor(self.freqs, snap_center=center) - for port, center in zip(self.ports, snap_centers) - ] - port_current_monitors = [ - port.to_current_monitor(self.freqs, snap_center=center) - for port, center in zip(self.ports, snap_centers) - ] - lumped_resistors = [ - port.to_load(snap_center=center) for port, center in zip(self.ports, snap_centers) - ] - new_mnts = ( - list(self.simulation.monitors) + port_voltage_monitors + port_current_monitors + snap_centers[port.name] = new_port_center + + # Create monitors and snap to the center positions + field_monitors = [ + mon + for port in self.ports + for mon in port.to_field_monitors(self.freqs, snap_center=snap_centers.get(port.name)) + ] + + new_mnts = list(self.simulation.monitors) + field_monitors + + lumped_resistors = [ + port.to_load(snap_center=snap_centers[port.name]) for port in self._lumped_ports + ] + + update_dict = dict( + monitors=new_mnts, + lumped_elements=lumped_resistors, + ) + + # This is the new default simulation will all shared components added + sim_wo_source = sim_wo_source.copy(update=update_dict) + + # Next, simulations are generated that include the source corresponding with the excitation port + for port in self._lumped_ports: + port_source = port.to_source( + self._source_time, snap_center=snap_centers[port.name], grid=sim_wo_source.grid ) - for port, center in zip(self.ports, snap_centers): - update_dict = dict( - monitors=new_mnts, - lumped_elements=lumped_resistors, - grid_spec=grid_spec, - ) - - sim_copy = self.simulation.copy(update=update_dict) - port_source = port.to_source( - self._source_time, snap_center=center, grid=sim_copy.grid - ) - sim_copy = sim_copy.updated_copy(sources=[port_source]) - task_name = self._task_name(port=port) - sim_dict[task_name] = sim_copy - - # Check final simulations for grid size at ports + task_name = self._task_name(port=port) + sim_dict[task_name] = sim_wo_source.updated_copy(sources=[port_source]) + + # Check final simulations for grid size at lumped ports for _, sim in sim_dict.items(): - TerminalComponentModeler._check_grid_size_at_ports(sim, self.ports) + TerminalComponentModeler._check_grid_size_at_ports(sim, self._lumped_ports) + + # Now, create simulations with wave port sources and mode solver monitors for computing port modes + for wave_port in self._wave_ports: + mode_monitor = wave_port.to_mode_solver_monitor(freqs=self.freqs) + # Source is placed just before the field monitor of the port + mode_src_pos = wave_port.center[wave_port.injection_axis] + self._shift_value_signed( + wave_port + ) + port_source = wave_port.to_source(self._source_time, snap_center=mode_src_pos) + + new_mnts_for_wave = new_mnts + [mode_monitor] + update_dict = dict(monitors=new_mnts_for_wave, sources=[port_source]) + + task_name = self._task_name(port=wave_port) + sim_dict[task_name] = sim_wo_source.copy(update=update_dict) return sim_dict @cached_property @@ -181,51 +180,65 @@ def _source_time(self): freq0=freq0, fwidth=fwidth, remove_dc_component=self.remove_dc_component ) - def _construct_smatrix(self, batch_data: BatchData) -> LumpedPortDataArray: + def _construct_smatrix(self, batch_data: BatchData) -> TerminalPortDataArray: """Post process ``BatchData`` to generate scattering matrix.""" port_names = [port.name for port in self.ports] values = np.zeros( - (len(port_names), len(port_names), len(self.freqs)), + (len(self.freqs), len(port_names), len(port_names)), dtype=complex, ) coords = dict( + f=np.array(self.freqs), port_out=port_names, port_in=port_names, - f=np.array(self.freqs), ) - a_matrix = LumpedPortDataArray(values, coords=coords) + a_matrix = TerminalPortDataArray(values, coords=coords) b_matrix = a_matrix.copy(deep=True) + V_matrix = a_matrix.copy(deep=True) + I_matrix = a_matrix.copy(deep=True) - def port_ab(port: Union[LumpedPort, CoaxialLumpedPort], sim_data: SimulationData): - """Helper to compute the port incident and reflected power waves.""" - voltage = port.compute_voltage(sim_data) - current = port.compute_current(sim_data) + def port_VI(port_out: AbstractTerminalPort, sim_data: SimulationData): + """Helper to compute the port voltages and currents.""" + voltage = port_out.compute_voltage(sim_data) + current = port_out.compute_current(sim_data) + return voltage, current - power_a = (voltage + port.impedance * current) / 2 / np.sqrt(np.real(port.impedance)) - power_b = (voltage - port.impedance * current) / 2 / np.sqrt(np.real(port.impedance)) - return power_a, power_b + # Tabulate the reference impedances at each port and frequency + port_impedances = self.port_reference_impedances(batch_data=batch_data) # loop through source ports for port_in in self.ports: sim_data = batch_data[self._task_name(port=port_in)] - for port_out in self.ports: - a_out, b_out = port_ab(port_out, sim_data) - a_matrix.loc[ + V_out, I_out = port_VI(port_out, sim_data) + V_matrix.loc[ dict( port_in=port_in.name, port_out=port_out.name, ) - ] = a_out + ] = V_out - b_matrix.loc[ + I_matrix.loc[ dict( port_in=port_in.name, port_out=port_out.name, ) - ] = b_out + ] = I_out + + # Reshape arrays so that broadcasting can be used to make F and Z act as diagonal matrices at each frequency + # Ensure data arrays have the correct data layout + V_numpy = V_matrix.transpose(*TerminalPortDataArray._dims).values + I_numpy = I_matrix.transpose(*TerminalPortDataArray._dims).values + Z_numpy = port_impedances.transpose(*PortDataArray._dims).values.reshape( + (len(self.freqs), len(port_names), 1) + ) + F_numpy = TerminalComponentModeler._compute_F(Z_numpy) + + # Equation 4.67 - Pozar - Microwave Engineering 4ed + a_matrix.values = F_numpy * (V_numpy + Z_numpy * I_numpy) + b_matrix.values = F_numpy * (V_numpy - np.conj(Z_numpy) * I_numpy) s_matrix = self.ab_to_s(a_matrix, b_matrix) @@ -242,10 +255,104 @@ def _validate_3d_simulation(cls, val): return val @staticmethod - def _check_grid_size_at_ports( - simulation: Simulation, ports: list[Union[LumpedPort, CoaxialLumpedPort]] - ): + def _check_grid_size_at_ports(simulation: Simulation, ports: list[Union[AbstractLumpedPort]]): """Raises :class:`.SetupError` if the grid is too coarse at port locations""" yee_grid = simulation.grid.yee for port in ports: port._check_grid_size(yee_grid) + + @staticmethod + def ab_to_s( + a_matrix: TerminalPortDataArray, b_matrix: TerminalPortDataArray + ) -> TerminalPortDataArray: + """Get the scattering matrix given the power wave matrices.""" + # Ensure dimensions are ordered properly + a_matrix = a_matrix.transpose(*TerminalPortDataArray._dims) + b_matrix = b_matrix.transpose(*TerminalPortDataArray._dims) + + s_matrix = a_matrix.copy(deep=True) + a_vals = s_matrix.copy(deep=True).values + b_vals = b_matrix.copy(deep=True).values + + s_vals = np.matmul(b_vals, AbstractComponentModeler.inv(a_vals)) + + s_matrix.data = s_vals + return s_matrix + + @staticmethod + def s_to_z( + s_matrix: TerminalPortDataArray, reference: Union[complex, PortDataArray] + ) -> DataArray: + """Get the impedance matrix given the scattering matrix and a reference impedance.""" + + # Ensure dimensions are ordered properly + z_matrix = s_matrix.transpose(*TerminalPortDataArray._dims).copy(deep=True) + s_vals = z_matrix.values + eye = np.eye(len(s_matrix.port_out.values), len(s_matrix.port_in.values)) + if isinstance(reference, PortDataArray): + # From Equation 4.68 - Pozar - Microwave Engineering 4ed + # Ensure that Zport, F, and Finv act as diagonal matrices when multiplying by left or right + shape_left = (len(s_matrix.f), len(s_matrix.port_out), 1) + shape_right = (len(s_matrix.f), 1, len(s_matrix.port_in)) + Zport = reference.values.reshape(shape_right) + F = TerminalComponentModeler._compute_F(Zport).reshape(shape_right) + Finv = (1.0 / F).reshape(shape_left) + FinvSF = Finv * s_vals * F + RHS = eye * np.conj(Zport) + FinvSF * Zport + LHS = eye - FinvSF + z_vals = np.matmul(AbstractComponentModeler.inv(LHS), RHS) + else: + # Simpler case when all port impedances are the same + z_vals = ( + np.matmul(AbstractComponentModeler.inv(eye - s_vals), (eye + s_vals)) * reference + ) + + z_matrix.data = z_vals + return z_matrix + + def port_reference_impedances(self, batch_data: BatchData) -> PortDataArray: + """Tabulates the reference impedance of each port at each frequency.""" + port_names = [port.name for port in self.ports] + + values = np.zeros( + (len(self.freqs), len(port_names)), + dtype=complex, + ) + coords = dict(f=np.array(self.freqs), port=port_names) + port_impedances = PortDataArray(values, coords=coords) + for port in self.ports: + if isinstance(port, WavePort): + # Mode solver data for each wave port is stored in its associated SimulationData + sim_data_port = batch_data[self._task_name(port=port)] + # WavePorts have a port impedance calculated from its associated modal field distribution + # and is frequency dependent. + impedances = port.compute_port_impedance(sim_data_port).values + port_impedances.loc[dict(port=port.name)] = impedances.squeeze() + else: + # LumpedPorts have a constant reference impedance + port_impedances.loc[dict(port=port.name)] = np.full(len(self.freqs), port.impedance) + + port_impedances = TerminalComponentModeler._set_port_data_array_attributes(port_impedances) + return port_impedances + + @staticmethod + def _compute_F(Z_numpy: np.array): + """Helper to convert port impedance matrix to F, which is used for + computing generalized scattering parameters.""" + return 1.0 / (2.0 * np.sqrt(np.real(Z_numpy))) + + @cached_property + def _lumped_ports(self) -> list[AbstractLumpedPort]: + """A list of all lumped ports in the ``TerminalComponentModeler``""" + return [port for port in self.ports if isinstance(port, AbstractLumpedPort)] + + @cached_property + def _wave_ports(self) -> list[WavePort]: + """A list of all wave ports in the ``TerminalComponentModeler``""" + return [port for port in self.ports if isinstance(port, WavePort)] + + @staticmethod + def _set_port_data_array_attributes(data_array: PortDataArray) -> PortDataArray: + """Helper to set additional metadata for ``PortDataArray``.""" + data_array.name = "Z0" + return data_array.assign_attrs(units=OHM, long_name="characteristic impedance") diff --git a/tidy3d/plugins/smatrix/ports/base_lumped.py b/tidy3d/plugins/smatrix/ports/base_lumped.py index a353fcc8bf..60574f0615 100644 --- a/tidy3d/plugins/smatrix/ports/base_lumped.py +++ b/tidy3d/plugins/smatrix/ports/base_lumped.py @@ -1,55 +1,24 @@ """Class and custom data array for representing a scattering matrix port based on lumped circuit elements.""" -from abc import ABC, abstractmethod +from abc import abstractmethod from typing import Optional import pydantic.v1 as pd -from ....components.base import Tidy3dBaseModel, cached_property -from ....components.data.data_array import DataArray, FreqDataArray -from ....components.data.sim_data import SimulationData -from ....components.grid.grid import Grid, YeeGrid +from ....components.base import cached_property +from ....components.grid.grid import YeeGrid from ....components.lumped_element import AbstractLumpedResistor from ....components.monitor import FieldMonitor -from ....components.source import GaussianPulse, UniformCurrentSource from ....components.types import Complex, FreqArray from ....constants import OHM +from .base_terminal import AbstractTerminalPort DEFAULT_PORT_NUM_CELLS = 3 DEFAULT_REFERENCE_IMPEDANCE = 50 -class LumpedPortDataArray(DataArray): - """Port parameter matrix elements for lumped ports. - - Example - ------- - >>> import numpy as np - >>> ports_in = ['port1', 'port2'] - >>> ports_out = ['port1', 'port2'] - >>> f = [2e14] - >>> coords = dict( - ... port_in=ports_in, - ... port_out=ports_out, - ... f=f - ... ) - >>> fd = LumpedPortDataArray((1 + 1j) * np.random.random((2, 2, 1)), coords=coords) - """ - - __slots__ = () - _dims = ("port_out", "port_in", "f") - _data_attrs = {"long_name": "lumped port matrix element"} - - -class AbstractLumpedPort(Tidy3dBaseModel, ABC): - """Class representing a single lumped port""" - - name: str = pd.Field( - ..., - title="Name", - description="Unique name for the port.", - min_length=1, - ) +class AbstractLumpedPort(AbstractTerminalPort): + """Class representing a single lumped port.""" impedance: Complex = pd.Field( DEFAULT_REFERENCE_IMPEDANCE, @@ -74,37 +43,25 @@ def _voltage_monitor_name(self) -> str: def _current_monitor_name(self) -> str: return f"{self.name}_current" - @cached_property - @abstractmethod - def injection_axis(self): - """Injection axis of the port.""" - - @abstractmethod - def to_source( - self, source_time: GaussianPulse, snap_center: float, grid=Grid - ) -> UniformCurrentSource: - """Create a current source from the lumped port.""" - @abstractmethod - def to_load(self, snap_center: float) -> AbstractLumpedResistor: + def to_load(self, snap_center: float = None) -> AbstractLumpedResistor: """Create a load resistor from the lumped port.""" @abstractmethod - def to_voltage_monitor(self, freqs: FreqArray, snap_center: float) -> FieldMonitor: + def to_voltage_monitor(self, freqs: FreqArray, snap_center: float = None) -> FieldMonitor: """Field monitor to compute port voltage.""" @abstractmethod - def to_current_monitor(self, freqs: FreqArray, snap_center: float) -> FieldMonitor: + def to_current_monitor(self, freqs: FreqArray, snap_center: float = None) -> FieldMonitor: """Field monitor to compute port current.""" - @abstractmethod - def compute_voltage(self, sim_data: SimulationData) -> FreqDataArray: - """Helper to compute voltage across the port.""" - - @abstractmethod - def compute_current(self, sim_data: SimulationData) -> FreqDataArray: - """Helper to compute current flowing through the port.""" + def to_field_monitors(self, freqs: FreqArray, snap_center: float = None) -> list[FieldMonitor]: + """Field monitors to compute port voltage and current.""" + return [ + self.to_voltage_monitor(freqs, snap_center), + self.to_current_monitor(freqs, snap_center), + ] @abstractmethod - def _check_grid_size(yee_grid: YeeGrid): - """Raises :class:``SetupError`` if the grid is too coarse at port locations""" + def _check_grid_size(self, yee_grid: YeeGrid): + """Raises :class:`SetupError` if the grid is too coarse at port locations.""" diff --git a/tidy3d/plugins/smatrix/ports/base_terminal.py b/tidy3d/plugins/smatrix/ports/base_terminal.py new file mode 100644 index 0000000000..41d5d0ba2d --- /dev/null +++ b/tidy3d/plugins/smatrix/ports/base_terminal.py @@ -0,0 +1,72 @@ +"""Class and custom data array for representing a scattering-matrix port, which is defined by a pair of terminals.""" + +from abc import ABC, abstractmethod + +import pydantic.v1 as pd + +from ....components.base import Tidy3dBaseModel, cached_property +from ....components.data.data_array import DataArray, FreqDataArray +from ....components.data.sim_data import SimulationData +from ....components.grid.grid import Grid +from ....components.monitor import FieldMonitor +from ....components.source import GaussianPulse, Source +from ....components.types import FreqArray + + +class TerminalPortDataArray(DataArray): + """Port parameter matrix elements for terminal-based ports. + + Example + ------- + >>> import numpy as np + >>> ports_in = ['port1', 'port2'] + >>> ports_out = ['port1', 'port2'] + >>> f = [2e14] + >>> coords = dict( + ... f=f, + ... port_out=ports_out, + ... port_in=ports_in, + ... ) + >>> fd = TerminalPortDataArray((1 + 1j) * np.random.random((1, 2, 2)), coords=coords) + """ + + __slots__ = () + _dims = ("f", "port_out", "port_in") + _data_attrs = {"long_name": "terminal-based port matrix element"} + + +class AbstractTerminalPort(Tidy3dBaseModel, ABC): + """Class representing a single terminal-based port. All terminal ports must provide methods + for computing voltage and current. These quantities represent the voltage between the + terminals, and the current flowing from one terminal into the other. + """ + + name: str = pd.Field( + ..., + title="Name", + description="Unique name for the port.", + min_length=1, + ) + + @cached_property + @abstractmethod + def injection_axis(self): + """Injection axis of the port.""" + + @abstractmethod + def to_source( + self, source_time: GaussianPulse, snap_center: float = None, grid: Grid = None + ) -> Source: + """Create a current source from a terminal-based port.""" + + @abstractmethod + def to_field_monitors(self, freqs: FreqArray, snap_center: float = None) -> list[FieldMonitor]: + """Field monitors to compute port voltage and current.""" + + @abstractmethod + def compute_voltage(self, sim_data: SimulationData) -> FreqDataArray: + """Helper to compute voltage across the port.""" + + @abstractmethod + def compute_current(self, sim_data: SimulationData) -> FreqDataArray: + """Helper to compute current flowing into the port.""" diff --git a/tidy3d/plugins/smatrix/ports/coaxial_lumped.py b/tidy3d/plugins/smatrix/ports/coaxial_lumped.py index c4c01a430c..9a4982852d 100644 --- a/tidy3d/plugins/smatrix/ports/coaxial_lumped.py +++ b/tidy3d/plugins/smatrix/ports/coaxial_lumped.py @@ -1,4 +1,4 @@ -"""Lumped port specialization with an annuluar geometry for exciting coaxial ports.""" +"""Lumped port specialization with an annular geometry for exciting coaxial ports.""" import numpy as np import pydantic.v1 as pd @@ -21,9 +21,11 @@ from ...microwave.path_integrals import AbstractAxesRH from .base_lumped import AbstractLumpedPort +DEFAULT_COAX_SOURCE_NUM_POINTS = 11 + class CoaxialLumpedPort(AbstractLumpedPort, AbstractAxesRH): - """Class representing a single coaxial lumped port + """Class representing a single coaxial lumped port. Example ------- @@ -79,7 +81,7 @@ def main_axis(self): @cached_property def injection_axis(self): - """Required for inheriting from AbstractLumpedPort.""" + """Required for inheriting from AbstractTerminalPort.""" return self.normal_axis @pd.validator("center", always=True) @@ -103,7 +105,7 @@ def _ensure_inner_diameter_is_smaller(cls, val, values): return val def to_source( - self, source_time: GaussianPulse, snap_center: float, grid: Grid + self, source_time: GaussianPulse, snap_center: float = None, grid: Grid = None ) -> CustomCurrentSource: """Create a current source from the lumped port.""" # Discretized source amps are manually zeroed out later if they @@ -122,9 +124,14 @@ def to_source( size = [self.outer_diameter] * 3 size[self.injection_axis] = 0 bounding_box = Box(center=self.center, size=size) - inds = grid.discretize_inds(box=bounding_box) - num1 = inds[trans_axes[0]][1] - inds[trans_axes[0]][0] - num2 = inds[trans_axes[1]][1] - inds[trans_axes[1]][0] + + num1 = DEFAULT_COAX_SOURCE_NUM_POINTS + num2 = DEFAULT_COAX_SOURCE_NUM_POINTS + + if grid: + inds = grid.discretize_inds(box=bounding_box) + num1 = inds[trans_axes[0]][1] - inds[trans_axes[0]][0] + num2 = inds[trans_axes[1]][1] - inds[trans_axes[1]][0] # Get a normalized current density that is flowing radially from inner circle to outer circle # Total current is normalized to 1 @@ -187,10 +194,10 @@ def compute_coax_current(rin, rout, x, y): current_dataset=dataset_E, ) - def to_load(self, snap_center: float) -> CoaxialLumpedResistor: + def to_load(self, snap_center: float = None) -> CoaxialLumpedResistor: """Create a load resistor from the lumped port.""" # 2D materials are currently snapped to the grid, so snapping here is not needed. - # It is done here so plots of the simulation will more accurately portray the setup + # Snapping is done here so plots of the simulation will more accurately portray the setup. center = list(self.center) if snap_center: center[self.injection_axis] = snap_center @@ -204,7 +211,7 @@ def to_load(self, snap_center: float) -> CoaxialLumpedResistor: name=f"{self.name}_resistor", ) - def to_voltage_monitor(self, freqs: FreqArray, snap_center: float) -> FieldMonitor: + def to_voltage_monitor(self, freqs: FreqArray, snap_center: float = None) -> FieldMonitor: """Field monitor to compute port voltage.""" center = list(self.center) if snap_center: @@ -226,7 +233,7 @@ def to_voltage_monitor(self, freqs: FreqArray, snap_center: float) -> FieldMonit colocate=False, ) - def to_current_monitor(self, freqs: FreqArray, snap_center: float) -> FieldMonitor: + def to_current_monitor(self, freqs: FreqArray, snap_center: float = None) -> FieldMonitor: """Field monitor to compute port current.""" center = list(self.center) if snap_center: @@ -285,15 +292,7 @@ def compute_current(self, sim_data: SimulationData) -> FreqDataArray: # Loops around inner conductive circle conductor field_data = sim_data[self._current_monitor_name] - # Helper for generating x,y vertices around a circle in the local coordinate frame - def generate_circle_coordinates(radius, num_points): - angles = np.linspace(0, 2 * np.pi, num_points, endpoint=True) - xt = radius * np.cos(angles) - yt = radius * np.sin(angles) - return (xt, yt) - # Get transverse axes - trans_axes = self.remaining_axes (coord1, coord2, coord3) = self.local_dims # Just need a rough estimate for the number of cells @@ -309,31 +308,22 @@ def generate_circle_coordinates(radius, num_points): # One extra point is used to close the loop. num_path_coords = round(np.pi * num_coords / 4) * 4 + 1 - # These x,y coordinates are relate to the local coordinate frame - xt, yt = generate_circle_coordinates( - (self.outer_diameter + self.inner_diameter) / 4, num_path_coords - ) - xt += self.center[trans_axes[0]] - yt += self.center[trans_axes[1]] - - circle_vertices = np.column_stack((xt, yt)) - # Close the contour exactly - circle_vertices[-1, :] = circle_vertices[0, :] - # Get the coordinates normal to port and select positions just on either side of the port normal_coords = field_coords[coord3].values upper_bound = np.searchsorted(normal_coords, self.center[self.injection_axis]) lower_bound = upper_bound - 1 + center = list(self.center) # We need to choose which side of the port to place the path integral, # which depends on which way the inner conductor is extending if self.direction == "+": - path_pos = normal_coords[upper_bound] + center[self.injection_axis] = normal_coords[upper_bound] else: - path_pos = normal_coords[lower_bound] + center[self.injection_axis] = normal_coords[lower_bound] + radius = (self.outer_diameter + self.inner_diameter) / 4 # Setup the path integral and integrate the H field - path_integral = CustomCurrentIntegral2D( - axis=self.injection_axis, position=path_pos, vertices=circle_vertices + path_integral = CustomCurrentIntegral2D.from_circular_path( + self.center, radius, num_path_coords, self.injection_axis, False ) current = path_integral.compute_current(field_data) @@ -346,6 +336,7 @@ def generate_circle_coordinates(radius, num_points): @cached_property def _voltage_axis(self) -> Axis: + """Helper to return the chosen voltage axis. We arbitrarily choose the first in-plane axis.""" return self.remaining_axes[0] @cached_property diff --git a/tidy3d/plugins/smatrix/ports/rectangular_lumped.py b/tidy3d/plugins/smatrix/ports/rectangular_lumped.py index 9594a4b8a8..fd7db9f220 100644 --- a/tidy3d/plugins/smatrix/ports/rectangular_lumped.py +++ b/tidy3d/plugins/smatrix/ports/rectangular_lumped.py @@ -21,7 +21,7 @@ class LumpedPort(AbstractLumpedPort, Box): - """Class representing a single rectangular lumped port + """Class representing a single rectangular lumped port. Example ------- @@ -61,7 +61,7 @@ def current_axis(self) -> Axis: return 3 - self.injection_axis - self.voltage_axis def to_source( - self, source_time: GaussianPulse, snap_center: float, grid: Grid + self, source_time: GaussianPulse, snap_center: float = None, grid: Grid = None ) -> UniformCurrentSource: """Create a current source from the lumped port.""" # Discretized source amps are manually zeroed out later if they @@ -80,7 +80,7 @@ def to_source( confine_to_bounds=True, ) - def to_load(self, snap_center: float) -> LumpedResistor: + def to_load(self, snap_center: float = None) -> LumpedResistor: """Create a load resistor from the lumped port.""" # 2D materials are currently snapped to the grid, so snapping here is not needed. # It is done here so plots of the simulation will more accurately portray the setup @@ -96,7 +96,7 @@ def to_load(self, snap_center: float) -> LumpedResistor: voltage_axis=self.voltage_axis, ) - def to_voltage_monitor(self, freqs: FreqArray, snap_center: float) -> FieldMonitor: + def to_voltage_monitor(self, freqs: FreqArray, snap_center: float = None) -> FieldMonitor: """Field monitor to compute port voltage.""" center = list(self.center) if snap_center: @@ -117,7 +117,7 @@ def to_voltage_monitor(self, freqs: FreqArray, snap_center: float) -> FieldMonit colocate=False, ) - def to_current_monitor(self, freqs: FreqArray, snap_center: float) -> FieldMonitor: + def to_current_monitor(self, freqs: FreqArray, snap_center: float = None) -> FieldMonitor: """Field monitor to compute port current.""" center = list(self.center) if snap_center: @@ -249,7 +249,7 @@ def _select_within_bounds(coords: np.array, min, max): return (min_idx, max_idx - 1) def _check_grid_size(self, yee_grid: YeeGrid): - """Raises :class:``SetupError`` if the grid is too coarse at port locations.""" + """Raises :class:`SetupError` if the grid is too coarse at port locations""" e_component = "xyz"[self.voltage_axis] e_yee_grid = yee_grid.grid_dict[f"E{e_component}"] coords = e_yee_grid.to_dict[e_component] diff --git a/tidy3d/plugins/smatrix/ports/wave.py b/tidy3d/plugins/smatrix/ports/wave.py new file mode 100644 index 0000000000..6ee439de03 --- /dev/null +++ b/tidy3d/plugins/smatrix/ports/wave.py @@ -0,0 +1,187 @@ +"""Class and custom data array for representing a scattering matrix wave port.""" + +from typing import Optional, Union + +import numpy as np +import pydantic.v1 as pd + +from ....components.base import cached_property +from ....components.data.data_array import FreqDataArray, FreqModeDataArray +from ....components.data.monitor_data import ModeSolverData +from ....components.data.sim_data import SimulationData +from ....components.geometry.base import Box +from ....components.monitor import FieldMonitor, ModeSolverMonitor +from ....components.simulation import Simulation +from ....components.source import GaussianPulse, ModeSource, ModeSpec +from ....components.types import Bound, Direction, FreqArray +from ....exceptions import ValidationError +from ...microwave import CurrentIntegralTypes, ImpedanceCalculator, VoltageIntegralTypes +from ...mode import ModeSolver +from .base_terminal import AbstractTerminalPort + + +class WavePort(AbstractTerminalPort, Box): + """Class representing a single wave port""" + + direction: Direction = pd.Field( + ..., + title="Direction", + description="'+' or '-', defining which direction is considered 'input'.", + ) + + mode_spec: ModeSpec = pd.Field( + ModeSpec(), + title="Mode Specification", + description="Parameters to feed to mode solver which determine modes measured by monitor.", + ) + + mode_index: pd.NonNegativeInt = pd.Field( + 0, + title="Mode Index", + description="Index into the collection of modes returned by mode solver. " + " Specifies which mode to inject using this source. " + "If larger than ``mode_spec.num_modes``, " + "``num_modes`` in the solver will be set to ``mode_index + 1``.", + ) + + voltage_integral: Optional[VoltageIntegralTypes] = pd.Field( + None, + title="Voltage Integral", + description="Definition of voltage integral used to compute voltage and the characteristic impedance.", + ) + + current_integral: Optional[CurrentIntegralTypes] = pd.Field( + None, + title="Current Integral", + description="Definition of current integral used to compute current and the characteristic impedance.", + ) + + @cached_property + def injection_axis(self): + """Injection axis of the port.""" + return self.size.index(0.0) + + @cached_property + def _field_monitor_name(self) -> str: + """Return the name of the :class:`.FieldMonitor` associated with this port.""" + return f"{self.name}_field" + + @cached_property + def _mode_monitor_name(self) -> str: + """Return the name of the :class:`.ModeMonitor` associated with this port.""" + return f"{self.name}_mode" + + def to_source(self, source_time: GaussianPulse, snap_center: float = None) -> ModeSource: + """Create a mode source from the wave port.""" + center = list(self.center) + if snap_center: + center[self.injection_axis] = snap_center + return ModeSource( + center=center, + size=self.size, + source_time=source_time, + mode_spec=self.mode_spec, + mode_index=self.mode_index, + direction=self.direction, + name=self.name, + ) + + def to_field_monitors(self, freqs: FreqArray, snap_center: float = None) -> list[FieldMonitor]: + """Field monitor to compute port voltage and current.""" + center = list(self.center) + if snap_center: + center[self.injection_axis] = snap_center + field_mon = FieldMonitor( + center=center, + size=self.size, + freqs=freqs, + name=self._field_monitor_name, + colocate=False, + ) + return [field_mon] + + def to_mode_solver_monitor(self, freqs: FreqArray) -> ModeSolverMonitor: + """Mode solver monitor to compute modes that will be used to + compute characteristic impedances.""" + mode_mon = ModeSolverMonitor( + center=self.center, + size=self.size, + freqs=freqs, + name=self._mode_monitor_name, + colocate=False, + mode_spec=self.mode_spec, + direction=self.direction, + ) + return mode_mon + + def to_mode_solver(self, simulation: Simulation, freqs: FreqArray) -> ModeSolver: + """Helper to create a :class:`.ModeSolver` instance.""" + mode_solver = ModeSolver( + simulation=simulation, + plane=self, + mode_spec=self.mode_spec, + freqs=freqs, + direction=self.direction, + colocate=False, + ) + return mode_solver + + def compute_voltage(self, sim_data: SimulationData) -> FreqDataArray: + """Helper to compute voltage across the port.""" + field_monitor = sim_data[self._field_monitor_name] + return self.voltage_integral.compute_voltage(field_monitor) + + def compute_current(self, sim_data: SimulationData) -> FreqDataArray: + """Helper to compute current flowing through the port.""" + field_monitor = sim_data[self._field_monitor_name] + return self.current_integral.compute_current(field_monitor) + + def compute_port_impedance( + self, sim_mode_data: Union[SimulationData, ModeSolverData] + ) -> FreqModeDataArray: + """Helper to compute impedance of port. The port impedance is computed from the + transmission line mode, which should be TEM or at least quasi-TEM.""" + impedance_calc = ImpedanceCalculator( + voltage_integral=self.voltage_integral, current_integral=self.current_integral + ) + if isinstance(sim_mode_data, SimulationData): + mode_solver_data = sim_mode_data[self._mode_monitor_name] + else: + mode_solver_data = sim_mode_data + + # Filter out unwanted modes to reduce impedance computation effort + mode_solver_data = mode_solver_data._isel(mode_index=[self.mode_index]) + impedance_array = impedance_calc.compute_impedance(mode_solver_data) + return impedance_array + + @staticmethod + def _within_port_bounds(path_bounds: Bound, port_bounds: Bound) -> bool: + """Helper to check if one bounding box is completely within the other.""" + path_min = np.array(path_bounds[0]) + path_max = np.array(path_bounds[1]) + bound_min = np.array(port_bounds[0]) + bound_max = np.array(port_bounds[1]) + return (bound_min <= path_min).all() and (bound_max >= path_max).all() + + @pd.validator("voltage_integral", "current_integral") + def _validate_path_integrals_within_port(cls, val, values): + """Raise ``ValidationError`` when the supplied path integrals are not within the port bounds.""" + center = values["center"] + size = values["size"] + box = Box(center=center, size=size) + if val and not WavePort._within_port_bounds(val.bounds, box.bounds): + raise ValidationError( + f"'{cls.__name__}' must be setup with all path integrals defined within the bounds " + f"of the port. Path bounds are '{val.bounds}', but port bounds are '{box.bounds}'." + ) + return val + + @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( + "At least one of 'voltage_integral' or 'current_integral' must be provided." + ) + return val