diff --git a/.gitignore b/.gitignore index 556d1b05e..7da5e4365 100644 --- a/.gitignore +++ b/.gitignore @@ -129,3 +129,4 @@ htmlcov/ # Specific file and folder exclusions .idea +.vscode diff --git a/CHANGELOG.md b/CHANGELOG.md index 16964dc3a..1933e5229 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -6,6 +6,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ## [Unreleased] ### Added +- `Simulation` now accepts `LumpedElementType`, which currently only supports the `LumpedResistor` type. `LumpedPort` together with `LumpedResistor` make up the new `TerminalComponentModeler` in the `smatrix` plugin. - Uniaxial medium Lithium niobate to material library. - Added support for conformal mesh methods near PEC structures that can be specified through the field `pec_conformal_mesh_spec` in the `Simulation` class. diff --git a/docs/api/abstract_models.rst b/docs/api/abstract_models.rst index ab3bcdbc3..72f26a8b1 100644 --- a/docs/api/abstract_models.rst +++ b/docs/api/abstract_models.rst @@ -40,6 +40,7 @@ These are some classes that are used to organize the tidy3d components, but aren tidy3d.components.monitor.AbstractFluxMonitor tidy3d.components.monitor.PlanarMonitor tidy3d.components.monitor.AbstractFieldProjectionMonitor + tidy3d.components.lumped_element.LumpedElement tidy3d.components.grid.grid_spec.GridSpec1d tidy3d.components.data.data_array.DataArray tidy3d.components.data.monitor_data.MonitorData diff --git a/docs/api/index.rst b/docs/api/index.rst index 816cecb31..662ba2ecf 100644 --- a/docs/api/index.rst +++ b/docs/api/index.rst @@ -14,6 +14,7 @@ API |:computer:| sources monitors mode + lumped_elements discretization output_data scene @@ -34,6 +35,7 @@ API |:computer:| .. include:: /api/sources.rst .. include:: /api/monitors.rst .. include:: /api/mode.rst +.. include:: /api/lumped_elements.rst .. include:: /api/discretization.rst .. include:: /api/output_data.rst .. include:: /api/scene.rst diff --git a/docs/api/lumped_elements.rst b/docs/api/lumped_elements.rst new file mode 100644 index 000000000..6bc5c7410 --- /dev/null +++ b/docs/api/lumped_elements.rst @@ -0,0 +1,12 @@ +.. currentmodule:: tidy3d + +Lumped elements +=============== + +Passive elements +---------------- +.. autosummary:: + :toctree: _autosummary/ + :template: module.rst + + tidy3d.LumpedResistor diff --git a/docs/faq b/docs/faq index c3ebbbe99..b9430b7d0 160000 --- a/docs/faq +++ b/docs/faq @@ -1 +1 @@ -Subproject commit c3ebbbe99d7d4cc5255d58131927c5cc1c18ff11 +Subproject commit b9430b7d0244d3a3b699eb0d5752ab1783a6d51a diff --git a/tests/sims/simulation_sample.json b/tests/sims/simulation_sample.json index 3df0c16a5..5dfb578f5 100644 --- a/tests/sims/simulation_sample.json +++ b/tests/sims/simulation_sample.json @@ -1015,14 +1015,14 @@ }, "transform": [ [ - 0.9659258262890683, - -0.25881904510252074, + 0.9659258262890684, + -0.2588190451025208, 0.0, 0.0 ], [ - 0.25881904510252074, - 0.9659258262890683, + 0.2588190451025208, + 0.9659258262890684, 0.0, 0.0 ], @@ -1084,6 +1084,7 @@ "remove_dc_component": true }, "interpolate": true, + "confine_to_bounds": false, "polarization": "Hx" }, { @@ -1109,6 +1110,7 @@ "remove_dc_component": true }, "interpolate": true, + "confine_to_bounds": false, "polarization": "Ex" }, { @@ -1302,6 +1304,7 @@ "remove_dc_component": true }, "interpolate": true, + "confine_to_bounds": false, "current_dataset": { "type": "FieldDataset", "Ex": "ScalarFieldDataArray", @@ -1366,6 +1369,7 @@ } }, "interpolate": true, + "confine_to_bounds": false, "polarization": "Hx" } ], @@ -1510,7 +1514,7 @@ 1, 1 ], - "colocate": true, + "colocate": 1, "freqs": [ 200000000000000.0, 250000000000000.0 @@ -1542,7 +1546,7 @@ 1, 1 ], - "colocate": true, + "colocate": 1, "start": 0.0, "stop": null, "interval": 1, @@ -1644,7 +1648,7 @@ 1, 1 ], - "colocate": true, + "colocate": 1, "freqs": [ 250000000000000.0, 300000000000000.0 @@ -1794,7 +1798,7 @@ 1, 1 ], - "colocate": true, + "colocate": 1, "freqs": [ 250000000000000.0, 300000000000000.0 @@ -1851,7 +1855,7 @@ 1, 1 ], - "colocate": true, + "colocate": 1, "freqs": [ 250000000000000.0, 300000000000000.0 @@ -1905,7 +1909,7 @@ 1, 1 ], - "colocate": true, + "colocate": 1, "freqs": [ 250000000000000.0, 300000000000000.0 @@ -2225,6 +2229,7 @@ "version": "2.6.0", "courant": 0.8, "normalize_index": 0, + "lumped_elements": [], "shutoff": 0.0001, "subpixel": false, "run_time": 1e-12 diff --git a/tests/test_components/test_medium.py b/tests/test_components/test_medium.py index 0916ea286..4eba082b1 100644 --- a/tests/test_components/test_medium.py +++ b/tests/test_components/test_medium.py @@ -739,3 +739,44 @@ def test_nonlinear_medium(log_capture): medium2d = td.Medium2D(ss=medium, tt=medium) with pytest.raises(ValidationError): medium2d = td.Medium2D(ss=modulated, tt=modulated) + + +def test_lumped_resistor(): + resistor = td.LumpedResistor( + resistance=50.0, + center=[0, 0, 0], + size=[2, 0, 3], + voltage_axis=0, + name="R", + ) + _ = resistor.sheet_conductance + normal_axis = resistor.normal_axis + assert normal_axis == 1 + + # error if voltage axis is not in plane with the resistor + with pytest.raises(pydantic.ValidationError): + _ = td.LumpedResistor( + resistance=50.0, + center=[0, 0, 0], + size=[2, 0, 3], + voltage_axis=1, + name="R", + ) + + # error if not planar + with pytest.raises(pydantic.ValidationError): + _ = td.LumpedResistor( + resistance=50.0, + center=[0, 0, 0], + size=[0, 0, 3], + voltage_axis=2, + name="R", + ) + with pytest.raises(pydantic.ValidationError): + _ = td.LumpedResistor( + resistance=50.0, + center=[0, 0, 0], + size=[2, 1, 3], + voltage_axis=2, + name="R", + ) diff --git a/tests/test_components/test_simulation.py b/tests/test_components/test_simulation.py index e048bd1cc..8b91f1ed2 100644 --- a/tests/test_components/test_simulation.py +++ b/tests/test_components/test_simulation.py @@ -2623,3 +2623,82 @@ def test_advanced_material_intersection(): ) # it's ok if these are both present as long as they don't intersect sim = sim.updated_copy(structures=[struct1, struct2]) + + +def test_num_lumped_elements(): + """Make sure we error if too many lumped elements supplied.""" + + resistor = td.LumpedResistor( + size=(0, 1, 2), center=(0, 0, 0), name="R1", voltage_axis=2, resistance=75 + ) + grid_spec = td.GridSpec.auto(wavelength=1.0) + + _ = td.Simulation( + size=(5, 5, 5), + grid_spec=grid_spec, + structures=[], + lumped_elements=[resistor] * MAX_NUM_MEDIUMS, + run_time=1e-12, + ) + with pytest.raises(pydantic.ValidationError): + _ = td.Simulation( + size=(5, 5, 5), + grid_spec=grid_spec, + structures=[], + lumped_elements=[resistor] * (MAX_NUM_MEDIUMS + 1), + run_time=1e-12, + ) + + +def test_validate_lumped_elements(): + resistor = td.LumpedResistor( + size=(0, 1, 2), center=(0, 0, 0), name="R1", voltage_axis=2, resistance=75 + ) + + _ = td.Simulation( + size=(1, 2, 3), + run_time=1e-12, + grid_spec=td.GridSpec.uniform(dl=0.1), + lumped_elements=[resistor], + ) + # error for 1D/2D simulation with lumped elements + with pytest.raises(pydantic.ValidationError): + td.Simulation( + size=(1, 0, 3), + run_time=1e-12, + grid_spec=td.GridSpec.uniform(dl=0.1), + lumped_elements=[resistor], + ) + + with pytest.raises(pydantic.ValidationError): + td.Simulation( + size=(1, 0, 0), + run_time=1e-12, + grid_spec=td.GridSpec.uniform(dl=0.1), + lumped_elements=[resistor], + ) + + +def test_suggested_mesh_overrides(): + resistor = td.LumpedResistor( + size=(0, 1, 2), center=(0, 0, 0), name="R1", voltage_axis=2, resistance=75 + ) + sim = td.Simulation( + size=(1, 2, 3), + run_time=1e-12, + grid_spec=td.GridSpec.uniform(dl=0.1), + lumped_elements=[resistor], + ) + + suggested_mesh_overrides = sim.suggest_mesh_overrides() + assert len(suggested_mesh_overrides) == 2 + grid_spec = sim.grid_spec.copy( + update={ + "override_structures": list(sim.grid_spec.override_structures) + + suggested_mesh_overrides, + } + ) + + _ = sim.updated_copy( + grid_spec=grid_spec, + ) diff --git a/tests/test_plugins/terminal_component_modeler_def.py b/tests/test_plugins/terminal_component_modeler_def.py new file mode 100644 index 000000000..f0c2bf7b5 --- /dev/null +++ b/tests/test_plugins/terminal_component_modeler_def.py @@ -0,0 +1,148 @@ +import numpy as np + +import tidy3d as td +from tidy3d.plugins.smatrix import ( + LumpedPort, + TerminalComponentModeler, +) + +# Microstrip dimensions +unit = 1e6 +default_strip_length = 75e-3 * unit +strip_width = 3e-3 * unit +gap = 1e-3 * unit +gnd_width = strip_width * 8 +metal_thickness = 0.2e-3 * unit + +# Microstrip materials +pec = td.PECMedium() +pec_cond = td.Medium(conductivity=1e10) +pec2d = td.Medium2D(ss=pec_cond, tt=pec_cond) +diel = td.Medium(permittivity=4.4) + +# Frequency setup +freq_start = 1e8 +freq_stop = 10e9 + + +def make_simulation(planar_pec: bool, length: float = None): + if length: + strip_length = length + else: + strip_length = default_strip_length + + if planar_pec: + height = 0 + metal = pec2d + else: + height = metal_thickness + metal = pec + + # wavelength / frequency + freq0 = (freq_start + freq_stop) / 2 + fwidth = freq_stop - freq_start + wavelength0 = td.C_0 / freq0 + run_time = 60 / fwidth + + # Spatial grid specification + grid_spec = td.GridSpec.auto(min_steps_per_wvl=10, wavelength=td.C_0 / freq_stop) + + # Make structures + strip = td.Structure( + geometry=td.Box( + center=[0, 0, height + gap + height / 2], + size=[strip_length, strip_width, height], + ), + medium=metal, + ) + + ground = td.Structure( + geometry=td.Box( + center=[0, 0, height / 2], + size=[strip_length, gnd_width, height], + ), + medium=metal, + ) + + substrate = td.Structure( + geometry=td.Box( + center=[0, 0, height + gap / 2], + size=[strip_length, gnd_width, gap], + ), + medium=diel, + ) + + structures = [substrate, strip, ground] + + # Make simulation + center_sim = [0, 0, height + gap / 2 + gap * 2] + size_sim = [ + strip_length + 0.5 * wavelength0, + gnd_width + 0.5 * wavelength0, + 2 * height + gap + 0.5 * wavelength0, + ] + + sim = td.Simulation( + center=center_sim, + size=size_sim, + grid_spec=grid_spec, + structures=structures, + sources=[], + monitors=[], + run_time=run_time, + boundary_spec=td.BoundarySpec.all_sides(boundary=td.PML()), + shutoff=1e-4, + ) + + return sim + + +def make_component_modeler( + planar_pec: bool, reference_impedance: complex = 50, length: float = None, **kwargs +): + if length: + strip_length = length + else: + strip_length = default_strip_length + + sim = make_simulation(planar_pec, length=length) + + if planar_pec: + height = 0 + else: + height = metal_thickness + + center_src1 = [-strip_length / 2, 0, height + gap / 2] + size_src1 = [0, strip_width, gap] + + center_src2 = [strip_length / 2, 0, height + gap / 2] + size_src2 = [0, strip_width, gap] + + port_cells = np.ceil(gap / (metal_thickness / 1)) + + port_1 = LumpedPort( + center=center_src1, + size=size_src1, + voltage_axis=2, + name="lumped_port_1", + num_grid_cells=port_cells, + impedance=reference_impedance, + ) + + port_2 = LumpedPort( + center=center_src2, + size=size_src2, + voltage_axis=2, + name="lumped_port_2", + num_grid_cells=port_cells, + impedance=reference_impedance, + ) + + ports = [port_1, port_2] + freqs = np.linspace(freq_start, freq_stop, 100) + + modeler = TerminalComponentModeler( + simulation=sim, ports=ports, freqs=freqs, remove_dc_component=False, verbose=True, **kwargs + ) + + return modeler diff --git a/tests/test_plugins/test_component_modeler.py b/tests/test_plugins/test_component_modeler.py index d05066d31..16e21c482 100644 --- a/tests/test_plugins/test_component_modeler.py +++ b/tests/test_plugins/test_component_modeler.py @@ -6,7 +6,10 @@ import tidy3d as td from tidy3d.web.api.container import Batch -from tidy3d.plugins.smatrix.smatrix import Port, ComponentModeler +from tidy3d.plugins.smatrix import ( + Port, + ComponentModeler, +) from tidy3d.exceptions import SetupError, Tidy3dKeyError from ..utils import run_emulated @@ -59,7 +62,9 @@ def make_coupler( # bend interpolator interp = tanh_interp(3) delta = wg_width + wg_spacing_coup - wg_spacing_in - offset = lambda u: wg_spacing_in + interp(u) * delta + + def offset(u): + return wg_spacing_in + interp(u) * delta coup = gdstk.RobustPath( (-0.5 * length, 0), @@ -213,7 +218,7 @@ def test_validate_no_sources(tmp_path): def test_element_mappings_none(tmp_path): modeler = make_component_modeler(path_dir=str(tmp_path)) modeler = modeler.updated_copy(ports=[], element_mappings=()) - modeler.matrix_indices_run_sim + _ = modeler.matrix_indices_run_sim def test_no_port(tmp_path): @@ -308,11 +313,11 @@ def test_component_modeler_run_only(monkeypatch): def _test_mappings(element_mappings, s_matrix): """Makes sure the mappings are reflected in a given S matrix.""" - for (i, j), (k, l), mult_by in element_mappings: + for (i, j), (k, L), mult_by in element_mappings: (port_out_from, mode_index_out_from) = i (port_in_from, mode_index_in_from) = j (port_out_to, mode_index_out_to) = k - (port_in_to, mode_index_in_to) = l + (port_in_to, mode_index_in_to) = L coords_from = dict( port_in=port_in_from, @@ -376,3 +381,4 @@ def test_mapping_exclusion(monkeypatch, tmp_path): def test_batch_filename(tmp_path): modeler = make_component_modeler(path_dir=str(tmp_path)) path = modeler._batch_path + assert path diff --git a/tests/test_plugins/test_terminal_component_modeler.py b/tests/test_plugins/test_terminal_component_modeler.py new file mode 100644 index 000000000..765a4cd55 --- /dev/null +++ b/tests/test_plugins/test_terminal_component_modeler.py @@ -0,0 +1,140 @@ +import pytest +import numpy as np +import pydantic.v1 as pydantic +import matplotlib.pyplot as plt + +import tidy3d as td +from tidy3d.plugins.smatrix import ( + AbstractComponentModeler, + LumpedPortDataArray, + TerminalComponentModeler, +) +from tidy3d.exceptions import Tidy3dKeyError +from ..utils import run_emulated +from .terminal_component_modeler_def import make_component_modeler + + +def run_component_modeler(monkeypatch, modeler: TerminalComponentModeler): + sim_dict = modeler.sim_dict + 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))) + s_matrix = modeler.run(path_dir=modeler.path_dir) + return s_matrix + + +def test_validate_no_sources(tmp_path): + modeler = make_component_modeler(planar_pec=True, path_dir=str(tmp_path)) + source = td.PointDipole( + source_time=td.GaussianPulse(freq0=2e14, fwidth=1e14), polarization="Ex" + ) + sim_w_source = modeler.simulation.copy(update=dict(sources=(source,))) + with pytest.raises(pydantic.ValidationError): + _ = modeler.copy(update=dict(simulation=sim_w_source)) + + +def test_no_port(tmp_path): + modeler = make_component_modeler(planar_pec=True, path_dir=str(tmp_path)) + _ = modeler.ports + with pytest.raises(Tidy3dKeyError): + modeler.get_port_by_name(port_name="NOT_A_PORT") + + +def test_plot_sim(tmp_path): + modeler = make_component_modeler(planar_pec=False, path_dir=str(tmp_path)) + modeler.plot_sim(z=0) + plt.close() + + +def test_plot_sim_eps(tmp_path): + modeler = make_component_modeler(planar_pec=False, path_dir=str(tmp_path)) + modeler.plot_sim_eps(z=0) + plt.close() + + +def test_make_component_modeler(tmp_path): + _ = make_component_modeler(planar_pec=False, path_dir=str(tmp_path)) + + +def test_run(monkeypatch, tmp_path): + modeler = make_component_modeler(planar_pec=True, path_dir=str(tmp_path)) + monkeypatch.setattr(TerminalComponentModeler, "run", lambda self, path_dir: None) + modeler.run(path_dir=str(tmp_path)) + + +def test_run_component_modeler(monkeypatch, tmp_path): + modeler = make_component_modeler(planar_pec=True, path_dir=str(tmp_path)) + s_matrix = run_component_modeler(monkeypatch, modeler) + + 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) != 0), "source index not present in S matrix" + assert np.all( + s_matrix.sel(**coords_in).sel(**coords_out) != 0 + ), "monitor index not present in S matrix" + + +def test_s_to_z_component_modeler(): + # Test case is 2 port T network with reference impedance of 50 Ohm + A = 20 + 30j + B = 50 - 15j + C = 60 + + Z11 = A + C + Z21 = C + Z12 = C + Z22 = B + C + + Z0 = 50.0 + # Manual creation of S parameters Pozar Table 4.2 + deltaZ = (Z11 + Z0) * (Z22 + Z0) - Z12 * Z21 + S11 = ((Z11 - Z0) * (Z22 + Z0) - Z12 * Z21) / deltaZ + S12 = (2 * Z12 * Z0) / deltaZ + S21 = (2 * Z21 * Z0) / deltaZ + S22 = ((Z11 + Z0) * (Z22 - Z0) - Z12 * Z21) / deltaZ + + port_names = ["lumped_port_1", "lumped_port_2"] + freqs = [1e8] + + values = np.array( + [[[S11, S12], [S21, S22]]], + dtype=complex, + ) + # Put coords in opposite order to check reordering + coords = dict( + f=np.array(freqs), + port_out=port_names, + port_in=port_names, + ) + + s_matrix = LumpedPortDataArray(data=values, coords=coords) + z_matrix = AbstractComponentModeler.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) + + +def test_ab_to_s_component_modeler(): + coords = dict( + f=np.array([1e8]), + port_out=["lumped_port_1", "lumped_port_2"], + port_in=["lumped_port_1", "lumped_port_2"], + ) + # Common case is reference impedance matched to loads, which means ideally + # the a matrix would be an identity matrix, and as a result the s matrix will be + # given directly by the b_matrix + 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) + assert np.isclose(S_matrix, b_matrix).all() diff --git a/tidy3d/__init__.py b/tidy3d/__init__.py index 7e388179e..59c951e88 100644 --- a/tidy3d/__init__.py +++ b/tidy3d/__init__.py @@ -52,6 +52,9 @@ from .components.monitor import FieldProjectionKSpaceMonitor, FieldProjectionSurface from .components.monitor import DiffractionMonitor +# lumped elements +from .components.lumped_element import LumpedResistor + # simulation from .components.simulation import Simulation @@ -300,6 +303,7 @@ def set_logging_level(level: str) -> None: "config", "__version__", "Updater", + "LumpedResistor", "Scene", "StructureStructureInterface", "StructureBoundary", diff --git a/tidy3d/components/geometry/utils_2d.py b/tidy3d/components/geometry/utils_2d.py index 9da52ab49..afff6c123 100644 --- a/tidy3d/components/geometry/utils_2d.py +++ b/tidy3d/components/geometry/utils_2d.py @@ -4,11 +4,12 @@ from typing import Tuple, List from ..types import Axis -from ...constants import fp_eps, inf +from ...constants import inf from ...exceptions import ValidationError from ..geometry.base import Geometry, Box, ClipOperation from ..geometry.primitives import Cylinder from ..geometry.polyslab import PolySlab +from ..grid.grid import Grid from ..scene import Scene from ..structure import Structure @@ -17,6 +18,31 @@ DIST_NEIGHBOR_REL_2D_MED = 1e-5 +def increment_float(val: np.float32, sign) -> np.float32: + """Applies a small positive or negative shift to a 32bit float using numpy.nextafter,""" + """but additionally handles some corner cases.""" + # Infinity is left unchanged + if val == inf or val == -inf: + return val + + if sign >= 0: + sign = 1 + else: + sign = -1 + # Numpy seems to skip over the increment from -0.0 and +0.0 + # which is different from c++ + val_inc = np.nextafter(val, sign * inf, dtype=np.float32) + + return np.float32(val_inc) + + +def snap_coordinate_to_grid(grid: Grid, center: float, axis: Axis) -> float: + """2D materials are snapped to grid along their normal axis""" + new_centers = grid.boundaries.to_list[axis] + new_center = new_centers[np.argmin(abs(new_centers - center))] + return new_center + + def get_bounds(geom: Geometry, axis: Axis) -> Tuple[float, float]: """Get the bounds of a geometry in the axis direction.""" return (geom.bounds[0][axis], geom.bounds[1][axis]) @@ -74,10 +100,8 @@ def get_neighbors( bounds = [list(i) for i in geom_shifted.bounds] _, tan_dirs = Geometry.pop_axis([0, 1, 2], axis=axis) for dim in tan_dirs: - if bounds[0][dim] != -inf: - bounds[0][dim] += fp_eps * max(np.abs(bounds[0][dim]), 1.0) - if bounds[1][dim] != inf: - bounds[1][dim] -= fp_eps * max(np.abs(bounds[1][dim]), 1.0) + bounds[0][dim] = increment_float(bounds[0][dim], 1.0) + bounds[1][dim] = increment_float(bounds[1][dim], -1.0) structures_side = Scene.intersecting_structures(Box.from_bounds(*bounds), structures) diff --git a/tidy3d/components/lumped_element.py b/tidy3d/components/lumped_element.py new file mode 100644 index 000000000..634495d5b --- /dev/null +++ b/tidy3d/components/lumped_element.py @@ -0,0 +1,130 @@ +"""Defines lumped elements that should be included in the simulation.""" +from __future__ import annotations +from abc import ABC +from typing import Union, List, Optional + +import pydantic.v1 as pydantic + +from .base import cached_property, skip_if_fields_missing + +from .geometry.base import Box +from ..components.structure import MeshOverrideStructure +from .types import Axis +from .viz import PlotParams, plot_params_lumped_element +from ..constants import OHM +from ..components.validators import validate_name_str, assert_plane +from ..exceptions import ValidationError + +DEFAULT_LUMPED_ELEMENT_NUM_CELLS = 3 + + +class LumpedElement(Box, ABC): + """Base class describing lumped elements.""" + + name: str = pydantic.Field( + ..., + title="Name", + description="Unique name for the lumped element.", + min_length=1, + ) + + voltage_axis: Axis = pydantic.Field( + ..., + title="Voltage Drop Axis", + description="Specifies the axis along which the component is oriented and along which the " + "associated voltage drop will occur. Must be in the plane of the element.", + ) + + num_grid_cells: Optional[pydantic.PositiveInt] = pydantic.Field( + DEFAULT_LUMPED_ELEMENT_NUM_CELLS, + title="Lumped element grid cells", + description="Number of mesh grid cells associated with the lumped element along each direction. " + "Used in generating the suggested list of :class:`MeshOverrideStructure` objects." + "A value of ``None`` will turn off mesh refinement suggestions.", + ) + + _name_validator = validate_name_str() + _plane_validator = assert_plane() + + @cached_property + def normal_axis(self): + """Normal axis of the lumped element.""" + return self.size.index(0.0) + + @cached_property + def plot_params(self) -> PlotParams: + """Default parameters for plotting a :class:`LumpedElement` object.""" + return plot_params_lumped_element + + def to_mesh_overrides(self) -> List[MeshOverrideStructure]: + """Creates a suggested :class:`MeshOverrideStructure` list that could be added to the :class:`Simulation`""" + # Create a mesh override when refinement is needed. + # An important use case is when a LumpedResistor is used with a LumpedPort + # The port is a flat surface, but when computing the port current, + # we'll eventually integrate the magnetic field just above and below + # this surface, so the mesh override needs to ensure that the mesh + # is fine enough not only in plane, but also in the normal direction. + # So in the normal direction, we'll make sure there are at least + # 2 cell layers above and below whose size is the same as the in-plane + # cell size in the override region. Also, to ensure that the port itself + # is aligned with a grid boundary in the normal direction, two separate + # override regions are defined, one above and one below the analytical + # port region. + mesh_overrides = [] + if self.num_grid_cells: + override_dl = [self.size[self.voltage_axis] / self.num_grid_cells] * 3 + override_dl[self.normal_axis] /= 2 + override_size = list(self.size) + override_size[override_size.index(0)] = 2 * override_dl[self.normal_axis] + override_center_above = list(self.center) + override_center_above[self.normal_axis] += override_dl[self.normal_axis] + override_center_below = list(self.center) + override_center_below[self.normal_axis] -= override_dl[self.normal_axis] + mesh_overrides.append( + MeshOverrideStructure( + geometry=Box(center=override_center_below, size=override_size), + dl=override_dl, + ) + ) + mesh_overrides.append( + MeshOverrideStructure( + geometry=Box(center=override_center_above, size=override_size), + dl=override_dl, + ) + ) + return mesh_overrides + + @pydantic.validator("voltage_axis", always=True) + @skip_if_fields_missing(["name", "size"]) + def _voltage_axis_in_plane(cls, val, values): + """Ensure voltage drop axis is in the plane of the lumped element.""" + name = values.get("name") + size = values.get("size") + if size.count(0.0) == 1 and size.index(0.0) == val: + # if not planar, then a separate validator should be triggered, not this one + raise ValidationError( + f"'voltage_axis' must be in the plane of lumped element '{name}'." + ) + return val + + +class LumpedResistor(LumpedElement): + """Class representing a lumped resistor. Lumped resistors are appended to the list of structures in the simulation + as :class:`Medium2D` with the appropriate conductivity given their size and voltage axis.""" + + resistance: pydantic.PositiveFloat = pydantic.Field( + ..., + title="Resistance", + description="Resistance value in ohms.", + unit=OHM, + ) + + @cached_property + def sheet_conductance(self): + """Effective sheet conductance.""" + lateral_axis = 3 - self.voltage_axis - self.normal_axis + return self.size[self.voltage_axis] / self.size[lateral_axis] / self.resistance + + +# lumped elements allowed in Simulation.lumped_elements +LumpedElementType = Union[LumpedResistor,] diff --git a/tidy3d/components/simulation.py b/tidy3d/components/simulation.py index 82f81adb0..f9e77337c 100644 --- a/tidy3d/components/simulation.py +++ b/tidy3d/components/simulation.py @@ -16,30 +16,31 @@ from .geometry.base import Geometry, Box from .geometry.mesh import TriangleMesh from .geometry.utils import flatten_groups, traverse_geometries -from .geometry.utils_2d import get_bounds, set_bounds, get_thickened_geom, subdivide +from .geometry.utils_2d import get_bounds, increment_float, set_bounds, get_thickened_geom +from .geometry.utils_2d import subdivide, snap_coordinate_to_grid from .types import Ax, FreqBound, Axis, annotate_type, InterpMethod, Symmetry from .types import Literal, TYPE_TAG_STR from .grid.grid import Coords1D, Grid, Coords from .grid.grid_spec import GridSpec, UniformGrid, AutoGrid, CustomGrid from .grid.grid_spec import ConformalMeshSpecType, StaircasingConformalMeshSpec from .medium import MediumType, AbstractMedium -from .medium import AbstractCustomMedium, Medium2D +from .medium import AbstractCustomMedium, Medium, Medium2D, MediumType3D from .medium import AnisotropicMedium, FullyAnisotropicMedium, AbstractPerturbationMedium from .boundary import BoundarySpec, BlochBoundary, PECBoundary, PMCBoundary, Periodic from .boundary import PML, StablePML, Absorber, AbsorberSpec -from .structure import Structure +from .structure import Structure, MeshOverrideStructure from .source import SourceType, PlaneWave, GaussianBeam, AstigmaticGaussianBeam, CustomFieldSource from .source import CustomCurrentSource, CustomSourceTime, ContinuousWave from .source import TFSF, Source, ModeSource -from .medium import Medium, MediumType3D from .monitor import MonitorType, Monitor, FreqMonitor, SurfaceIntegrationMonitor from .monitor import AbstractModeMonitor, FieldMonitor, TimeMonitor from .monitor import PermittivityMonitor, DiffractionMonitor, AbstractFieldProjectionMonitor from .monitor import FieldProjectionAngleMonitor, FieldProjectionKSpaceMonitor +from .lumped_element import LumpedElementType, LumpedResistor from .data.dataset import Dataset from .data.data_array import SpatialDataArray from .viz import add_ax_if_none, equal_aspect -from .scene import Scene +from .scene import Scene, MAX_NUM_MEDIUMS from .viz import PlotParams from .viz import plot_params_pml, plot_params_override_structures @@ -316,6 +317,42 @@ class Simulation(AbstractSimulation): * `Numerical dispersion in FDTD `_ """ + lumped_elements: Tuple[LumpedElementType, ...] = pydantic.Field( + (), + title="Lumped Elements", + description="Tuple of lumped elements in the simulation. " + "Note: only :class:`tidy3d.LumpedResistor` is supported currently.", + ) + """ + Tuple of lumped elements in the simulation. + + Example + ------- + Simple application reference: + + .. code-block:: python + + Simulation( + ... + lumped_elements=[ + LumpedResistor( + size=(0, 3, 1), + center=(0, 0, 0), + voltage_axis=2, + resistance=50, + name="resistor_1", + ) + ], + ... + ) + + See Also + -------- + + `Index <../lumped_elements.html>`_: + Available lumped element types. + """ + grid_spec: GridSpec = pydantic.Field( GridSpec(), title="Grid Specification", @@ -1572,6 +1609,35 @@ def _check_normalize_index(cls, val, values): return val + @pydantic.validator("lumped_elements", always=True) + @skip_if_fields_missing(["structures"]) + def _validate_num_lumped_elements(cls, val, values): + """Error if too many lumped elements present.""" + + if val is None: + return val + structures = values.get("structures") + mediums = {structure.medium for structure in structures} + total_num_mediums = len(val) + len(mediums) + if total_num_mediums > MAX_NUM_MEDIUMS: + raise ValidationError( + f"Tidy3D only supports {MAX_NUM_MEDIUMS} distinct lumped elements and structures." + f"{total_num_mediums} were supplied." + ) + + return val + + @pydantic.validator("lumped_elements") + @skip_if_fields_missing(["size"]) + def _check_3d_simulation_with_lumped_elements(cls, val, values): + """Error if Simulation contained lumped elements and is not a 3D simulation""" + size = values.get("size") + if val and size.count(0.0) > 0: + raise ValidationError( + f"'{cls.__name__}' must be a 3D simulation when a 'LumpedElement' is present." + ) + return val + """ Post-init validators """ def _post_init_validators(self) -> None: @@ -2360,6 +2426,7 @@ def plot( ax: Ax = None, source_alpha: float = None, monitor_alpha: float = None, + lumped_element_alpha: float = None, hlim: Tuple[float, float] = None, vlim: Tuple[float, float] = None, **patch_kwargs, @@ -2378,6 +2445,8 @@ def plot( Opacity of the sources. If ``None``, uses Tidy3d default. monitor_alpha : float = None Opacity of the monitors. If ``None``, uses Tidy3d default. + lumped_element_alpha : float = None + Opacity of the lumped elements. If ``None``, uses Tidy3d default. ax : matplotlib.axes._subplots.Axes = None Matplotlib axes to plot on, if not specified, one is created. hlim : Tuple[float, float] = None @@ -2404,6 +2473,9 @@ def plot( ax = self.plot_structures(ax=ax, x=x, y=y, z=z, hlim=hlim, vlim=vlim) ax = self.plot_sources(ax=ax, x=x, y=y, z=z, hlim=hlim, vlim=vlim, alpha=source_alpha) ax = self.plot_monitors(ax=ax, x=x, y=y, z=z, hlim=hlim, vlim=vlim, alpha=monitor_alpha) + ax = self.plot_lumped_elements( + ax=ax, x=x, y=y, z=z, hlim=hlim, vlim=vlim, alpha=lumped_element_alpha + ) ax = self.plot_symmetries(ax=ax, x=x, y=y, z=z, hlim=hlim, vlim=vlim) ax = self.plot_pml(ax=ax, x=x, y=y, z=z, hlim=hlim, vlim=vlim) ax = Scene._set_plot_bounds( @@ -2423,6 +2495,7 @@ def plot_eps( alpha: float = None, source_alpha: float = None, monitor_alpha: float = None, + lumped_element_alpha: float = None, hlim: Tuple[float, float] = None, vlim: Tuple[float, float] = None, ax: Ax = None, @@ -2448,6 +2521,8 @@ def plot_eps( Opacity of the sources. If ``None``, uses Tidy3d default. monitor_alpha : float = None Opacity of the monitors. If ``None``, uses Tidy3d default. + lumped_element_alpha : float = None + Opacity of the lumped elements. If ``None``, uses Tidy3d default. ax : matplotlib.axes._subplots.Axes = None Matplotlib axes to plot on, if not specified, one is created. hlim : Tuple[float, float] = None @@ -2484,6 +2559,9 @@ def plot_eps( ) ax = self.plot_sources(ax=ax, x=x, y=y, z=z, hlim=hlim, vlim=vlim, alpha=source_alpha) ax = self.plot_monitors(ax=ax, x=x, y=y, z=z, hlim=hlim, vlim=vlim, alpha=monitor_alpha) + ax = self.plot_lumped_elements( + ax=ax, x=x, y=y, z=z, hlim=hlim, vlim=vlim, alpha=lumped_element_alpha + ) ax = self.plot_symmetries(ax=ax, x=x, y=y, z=z, hlim=hlim, vlim=vlim) ax = self.plot_pml(ax=ax, x=x, y=y, z=z, hlim=hlim, vlim=vlim) ax = Scene._set_plot_bounds( @@ -2894,6 +2972,51 @@ def set_plot_params(boundary_edge, lim, side, thickness): return ax + @equal_aspect + @add_ax_if_none + def plot_lumped_elements( + self, + x: float = None, + y: float = None, + z: float = None, + hlim: Tuple[float, float] = None, + vlim: Tuple[float, float] = None, + alpha: float = None, + ax: Ax = None, + ) -> Ax: + """Plot each of simulation's lumped elements on a plane defined by one + nonzero x,y,z coordinate. + + Parameters + ---------- + x : float = None + position of plane in x direction, only one of x, y, z must be specified to define plane. + y : float = None + position of plane in y direction, only one of x, y, z must be specified to define plane. + z : float = None + position of plane in z direction, only one of x, y, z must be specified to define plane. + hlim : Tuple[float, float] = None + The x range if plotting on xy or xz planes, y range if plotting on yz plane. + vlim : Tuple[float, float] = None + The z range if plotting on xz or yz planes, y plane if plotting on xy plane. + alpha : float = None + Opacity of the lumped element, If ``None`` uses Tidy3d default. + ax : matplotlib.axes._subplots.Axes = None + Matplotlib axes to plot on, if not specified, one is created. + + Returns + ------- + matplotlib.axes._subplots.Axes + The supplied or created matplotlib axes. + """ + bounds = self.bounds + for element in self.lumped_elements: + ax = element.plot(x=x, y=y, z=z, alpha=alpha, ax=ax, sim_bounds=bounds) + ax = Scene._set_plot_bounds( + bounds=self.simulation_bounds, ax=ax, x=x, y=y, z=z, hlim=hlim, vlim=vlim + ) + return ax + @cached_property def frequency_range(self) -> FreqBound: """Range of frequencies spanning all sources' frequency dependence. @@ -3412,7 +3535,10 @@ def _volumetric_structures_grid(self, grid: Grid) -> Tuple[Structure]: """Generate a tuple of structures wherein any 2D materials are converted to 3D volumetric equivalents, using ``grid`` as the simulation grid.""" - if not any(isinstance(medium, Medium2D) for medium in self.scene.mediums): + if ( + not any(isinstance(medium, Medium2D) for medium in self.scene.mediums) + and not self.lumped_elements + ): return self.structures def get_dls(geom: Geometry, axis: Axis, num_dls: int) -> List[float]: @@ -3438,13 +3564,38 @@ def get_dls(geom: Geometry, axis: Axis, num_dls: int) -> List[float]: def snap_to_grid(geom: Geometry, axis: Axis) -> Geometry: """Snap a 2D material to the Yee grid.""" - new_centers = self._discretize_grid( - Box.from_bounds(*geom.bounds), grid=grid - ).boundaries.to_list[axis] - new_center = new_centers[np.argmin(abs(new_centers - get_bounds(geom, axis)[0]))] - return set_bounds(geom, (new_center, new_center), axis) + center = get_bounds(geom, axis)[0] + assert get_bounds(geom, axis)[0] == get_bounds(geom, axis)[1] + snapped_center = snap_coordinate_to_grid(self.grid, center, axis) + return set_bounds(geom, (snapped_center, snapped_center), axis) + + lumped_structures = [] + for lumped_element in self.lumped_elements: + _, tan_dirs = self.pop_axis([0, 1, 2], axis=lumped_element.normal_axis) + + if isinstance(lumped_element, LumpedResistor): + conductivity = lumped_element.sheet_conductance + + if tan_dirs[0] == lumped_element.voltage_axis: + medium_dict = { + "ss": Medium(conductivity=conductivity), + "tt": self.medium, + } + else: + medium_dict = { + "tt": Medium(conductivity=conductivity), + "ss": self.medium, + } + lumped_structures.append( + Structure( + geometry=Box(size=lumped_element.size, center=lumped_element.center), + medium=Medium2D(**medium_dict), + ) + ) # Begin volumetric structures grid + all_structures = list(self.structures) + lumped_structures + # For 1D and 2D simulations, a nonzero size is needed for the polygon operations in subdivide placeholder_size = tuple(i if i > 0 else inf for i in self.geometry.size) simulation_placeholder_geometry = self.geometry.updated_copy( @@ -3456,7 +3607,7 @@ def snap_to_grid(geom: Geometry, axis: Axis) -> Geometry: ) background_structures = [simulation_background] new_structures = [] - for structure in self.structures: + for structure in all_structures: if not isinstance(structure.medium, Medium2D): # found a 3D material; keep it background_structures.append(structure) @@ -3488,8 +3639,9 @@ def snap_to_grid(geom: Geometry, axis: Axis) -> Geometry: temp_structure = structure.updated_copy(geometry=temp_geometry, medium=new_medium) if structure.medium.is_pec: - pec_delta = fp_eps * max(np.abs(snapped_center), 1.0) - new_bounds = (snapped_center - pec_delta, snapped_center + pec_delta) + pec_plus = increment_float(snapped_center, 1.0) + pec_minus = increment_float(snapped_center, -1.0) + new_bounds = (pec_minus, pec_plus) new_geometry = set_bounds(snapped_geometry, bounds=new_bounds, axis=axis) new_structure = structure.updated_copy(geometry=new_geometry, medium=new_medium) @@ -3852,3 +4004,15 @@ def subsection( ) return new_sim + + def suggest_mesh_overrides(self, **kwargs) -> List[MeshOverrideStructure]: + """Generate a :class:.`MeshOverrideStructure` `List` which is automatically generated + from structures in the simulation. + """ + mesh_overrides = [] + + # For now we can suggest MeshOverrideStructures for lumped elements. + for lumped_element in self.lumped_elements: + mesh_overrides.extend(lumped_element.to_mesh_overrides()) + + return mesh_overrides diff --git a/tidy3d/components/source.py b/tidy3d/components/source.py index 804f6917f..ed3748b1d 100644 --- a/tidy3d/components/source.py +++ b/tidy3d/components/source.py @@ -468,6 +468,15 @@ class ReverseInterpolatedSource(Source): "placement at the specified location using linear interpolation.", ) + confine_to_bounds: bool = pydantic.Field( + False, + title="Confine to Analytical Bounds", + description="If ``True``, any source amplitudes which, after discretization, fall beyond " + "the bounding box of the source are zeroed out, but only along directions where " + "the source has a non-zero extent. The bounding box is inclusive. Should be set ```True`` " + "when the current source is being used to excite a current in a conductive material.", + ) + class UniformCurrentSource(CurrentSource, ReverseInterpolatedSource): """Source in a rectangular volume with uniform time dependence. diff --git a/tidy3d/components/viz.py b/tidy3d/components/viz.py index 1c04bdc57..7da9aa40d 100644 --- a/tidy3d/components/viz.py +++ b/tidy3d/components/viz.py @@ -116,6 +116,9 @@ def to_kwargs(self) -> dict: ) plot_params_fluid = PlotParams(facecolor="white", edgecolor="lightsteelblue", lw=0.4, hatch="xx") plot_params_grid = PlotParams(edgecolor="black", lw=0.2) +plot_params_lumped_element = PlotParams( + alpha=0.4, facecolor="mediumblue", edgecolor="mediumblue", lw=3 +) # stores color of simulation.structures for given index in simulation.medium_map MEDIUM_CMAP = [ diff --git a/tidy3d/constants.py b/tidy3d/constants.py index b0a8a7d41..ffbebf6e1 100644 --- a/tidy3d/constants.py +++ b/tidy3d/constants.py @@ -162,6 +162,11 @@ Picosecond per (nanometer kilometer). """ +OHM = "ohm" +""" +SI unit of resistance.. +""" + THERMAL_CONDUCTIVITY = "W/(um*K)" """ Watts per (micrometer Kelvin). diff --git a/tidy3d/plugins/smatrix/__init__.py b/tidy3d/plugins/smatrix/__init__.py index 1d88704b0..622614d47 100644 --- a/tidy3d/plugins/smatrix/__init__.py +++ b/tidy3d/plugins/smatrix/__init__.py @@ -1,5 +1,17 @@ """ Imports from scattering matrix plugin. """ -from .smatrix import ComponentModeler, Port, SMatrixDataArray +from .ports.modal import Port +from .ports.lumped import LumpedPort +from .component_modelers.modal import AbstractComponentModeler +from .component_modelers.modal import ComponentModeler, ModalPortDataArray +from .component_modelers.terminal import TerminalComponentModeler, LumpedPortDataArray -__all__ = ["ComponentModeler", "Port", "SMatrixDataArray"] +__all__ = [ + "AbstractComponentModeler", + "ComponentModeler", + "Port", + "ModalPortDataArray", + "TerminalComponentModeler", + "LumpedPort", + "LumpedPortDataArray", +] diff --git a/tidy3d/plugins/smatrix/component_modelers/base.py b/tidy3d/plugins/smatrix/component_modelers/base.py new file mode 100644 index 000000000..cc24bbf56 --- /dev/null +++ b/tidy3d/plugins/smatrix/component_modelers/base.py @@ -0,0 +1,228 @@ +"""Base class for generating an S matrix automatically from tidy3d simulations and port definitions.""" +from __future__ import annotations +from typing import Tuple, Dict, Union +import os +from abc import ABC, abstractmethod + +import pydantic.v1 as pd +import numpy as np + +from ....constants import HERTZ +from ....components.simulation import Simulation +from ....components.data.data_array import DataArray +from ....components.types import Complex, FreqArray +from ....components.base import Tidy3dBaseModel, cached_property +from ....exceptions import SetupError, Tidy3dKeyError +from ....log import log +from ....web.api.container import BatchData, Batch + +from ..ports.modal import Port +from ..ports.lumped import LumpedPort + +# fwidth of gaussian pulse in units of central frequency +FWIDTH_FRAC = 1.0 / 10 +DEFAULT_DATA_DIR = "." + + +class AbstractComponentModeler(ABC, Tidy3dBaseModel): + """Tool for modeling devices and computing port parameters.""" + + simulation: Simulation = pd.Field( + ..., + title="Simulation", + description="Simulation describing the device without any sources present.", + ) + + ports: Tuple[Union[Port, LumpedPort], ...] = pd.Field( + (), + title="Ports", + description="Collection of ports describing the scattering matrix elements. " + "For each input mode, one simulation will be run with a modal source.", + ) + + freqs: FreqArray = pd.Field( + ..., + title="Frequencies", + description="Array or list of frequencies at which to compute port parameters.", + units=HERTZ, + ) + + remove_dc_component: bool = pd.Field( + True, + title="Remove DC Component", + description="Whether to remove the DC component in the Gaussian pulse spectrum. " + "If ``True``, the Gaussian pulse is modified at low frequencies to zero out the " + "DC component, which is usually desirable so that the fields will decay. However, " + "for broadband simulations, it may be better to have non-vanishing source power " + "near zero frequency. Setting this to ``False`` results in an unmodified Gaussian " + "pulse spectrum which can have a nonzero DC component.", + ) + + folder_name: str = pd.Field( + "default", + title="Folder Name", + description="Name of the folder for the tasks on web.", + ) + + verbose: bool = pd.Field( + False, + title="Verbosity", + description="Whether the :class:`.ComponentModeler` should print status and progressbars.", + ) + + callback_url: str = pd.Field( + None, + title="Callback URL", + description="Http PUT url to receive simulation finish event. " + "The body content is a json file with fields " + "``{'id', 'status', 'name', 'workUnit', 'solverVersion'}``.", + ) + + path_dir: str = pd.Field( + DEFAULT_DATA_DIR, + title="Directory Path", + description="Base directory where data and batch will be downloaded.", + ) + + solver_version: str = pd.Field( + None, + title="Solver Version", + description_str="Custom solver version to use. " + "If not supplied, uses default for the current front end version.", + ) + + @pd.validator("simulation", always=True) + def _sim_has_no_sources(cls, val): + """Make sure simulation has no sources as they interfere with tool.""" + if len(val.sources) > 0: + raise SetupError("'AbstractComponentModeler.simulation' must not have any sources.") + return val + + @staticmethod + def _task_name(port: Port, mode_index: int = None) -> str: + """The name of a task, determined by the port of the source and mode index, if given.""" + if mode_index is not None: + return f"smatrix_{port.name}_{mode_index}" + return f"smatrix_{port.name}" + + @cached_property + def sim_dict(self) -> Dict[str, Simulation]: + """Generate all the :class:`Simulation` objects for the S matrix calculation.""" + + @cached_property + def batch(self) -> Batch: + """Batch associated with this component modeler.""" + + # first try loading the batch from file, if it exists + batch_path = self._batch_path + if os.path.exists(batch_path): + return Batch.from_file(fname=batch_path) + + return Batch( + simulations=self.sim_dict, + folder_name=self.folder_name, + callback_url=self.callback_url, + verbose=self.verbose, + solver_version=self.solver_version, + ) + + @cached_property + def batch_path(self) -> str: + """Path to the batch saved to file.""" + + return self.batch._batch_path(path_dir=DEFAULT_DATA_DIR) + + def get_path_dir(self, path_dir: str) -> None: + """Check whether the supplied 'path_dir' matches the internal field value.""" + + if path_dir not in (DEFAULT_DATA_DIR, self.path_dir): + log.warning( + f"'ComponentModeler' method was supplied a 'path_dir' of '{path_dir}' " + f"when its internal 'path_dir' field was set to '{self.path_dir}'. " + "The passed value will be deprecated in later versions. " + "Please set the internal 'path_dir' field to the desired value and " + "remove the 'path_dir' from the method argument. " + f"Using supplied '{path_dir}'." + ) + return path_dir + + return self.path_dir + + @cached_property + def _batch_path(self) -> str: + """Where we store the batch for this ComponentModeler instance after the run.""" + return os.path.join(self.path_dir, "batch" + str(hash(self)) + ".json") + + def _run_sims(self, path_dir: str = DEFAULT_DATA_DIR) -> BatchData: + """Run :class:`Simulations` for each port and return the batch after saving.""" + batch = self.batch + + batch_data = batch.run(path_dir=path_dir) + batch.to_file(self._batch_path) + return batch_data + + def get_port_by_name(self, port_name: str) -> Port: + """Get the port from the name.""" + ports = [port for port in self.ports if port.name == port_name] + if len(ports) == 0: + raise Tidy3dKeyError(f'Port "{port_name}" not found.') + return ports[0] + + @abstractmethod + def _construct_smatrix(self, batch_data: BatchData) -> DataArray: + """Post process :class:`BatchData` to generate scattering matrix.""" + + def run(self, path_dir: str = DEFAULT_DATA_DIR) -> DataArray: + """Solves for the scattering matrix of the system.""" + path_dir = self.get_path_dir(path_dir) + + batch_data = self._run_sims(path_dir=path_dir) + return self._construct_smatrix(batch_data=batch_data) + + def load(self, path_dir: str = DEFAULT_DATA_DIR) -> DataArray: + """Load a scattering matrix from saved :class:`BatchData` object.""" + path_dir = self.get_path_dir(path_dir) + + 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) diff --git a/tidy3d/plugins/smatrix/smatrix.py b/tidy3d/plugins/smatrix/component_modelers/modal.py similarity index 69% rename from tidy3d/plugins/smatrix/smatrix.py rename to tidy3d/plugins/smatrix/component_modelers/modal.py index 498106244..f77196241 100644 --- a/tidy3d/plugins/smatrix/smatrix.py +++ b/tidy3d/plugins/smatrix/component_modelers/modal.py @@ -1,83 +1,32 @@ -"""Tools for generating an S matrix automatically from tidy3d simulation and port definitions.""" +"""Tool for generating an S matrix automatically from a Tidy3d simulation and modal port definitions.""" +# TODO: The names "ComponentModeler" and "Port" should be changed to "ModalComponentModeler" and +# "ModalPort" to explicitly differentiate these from "TerminalComponentModeler" and "LumpedPort". from __future__ import annotations from typing import List, Tuple, Optional, Dict -import os import pydantic.v1 as pd import numpy as np -from ...constants import HERTZ -from ...components.simulation import Simulation -from ...components.geometry.base import Box -from ...components.mode import ModeSpec -from ...components.monitor import ModeMonitor -from ...components.source import ModeSource, GaussianPulse -from ...components.data.sim_data import SimulationData -from ...components.data.data_array import DataArray -from ...components.types import Direction, Ax, Complex, FreqArray -from ...components.viz import add_ax_if_none, equal_aspect -from ...components.base import Tidy3dBaseModel, cached_property -from ...exceptions import SetupError, Tidy3dKeyError -from ...log import log -from ...web.api.container import BatchData, Batch - -# fwidth of gaussian pulse in units of central frequency -FWIDTH_FRAC = 1.0 / 10 -DEFAULT_DATA_DIR = "." - - -class Port(Box): - """Specifies a port in the scattering matrix.""" - - direction: Direction = pd.Field( - ..., - title="Direction", - description="'+' or '-', defining which direction is considered 'input'.", - ) - mode_spec: ModeSpec = pd.Field( - ModeSpec(), - title="Mode Specification", - description="Specifies how the mode solver will solve for the modes of the port.", - ) - name: str = pd.Field( - ..., - title="Name", - description="Unique name for the port.", - min_length=1, - ) +from ....components.simulation import Simulation +from ....components.monitor import ModeMonitor +from ....components.source import ModeSource, GaussianPulse +from ....components.data.sim_data import SimulationData +from ....components.types import Ax, Complex +from ....components.viz import add_ax_if_none, equal_aspect +from ....components.base import cached_property +from ....exceptions import SetupError +from ....web.api.container import BatchData + +from .base import AbstractComponentModeler, FWIDTH_FRAC +from ..ports.modal import ModalPortDataArray, Port MatrixIndex = Tuple[str, pd.NonNegativeInt] # the 'i' in S_ij Element = Tuple[MatrixIndex, MatrixIndex] # the 'ij' in S_ij -class SMatrixDataArray(DataArray): - """Scattering matrix elements. - - Example - ------- - >>> port_in = ['port1', 'port2'] - >>> port_out = ['port1', 'port2'] - >>> mode_index_in = [0, 1] - >>> mode_index_out = [0, 1] - >>> f = [2e14] - >>> coords = dict( - ... port_in=ports_in, - ... port_out=ports_out, - ... mode_index_in=mode_index_in, - ... mode_index_out=mode_index_out, - ... f=f - ... ) - >>> fd = SMatrixDataArray((1 + 1j) * np.random.random((2, 2, 2, 2, 1)), coords=coords) - """ - - __slots__ = () - _dims = ("port_out", "mode_index_out", "port_in", "mode_index_in", "f") - _data_attrs = {"long_name": "scattering matrix element"} - - -class ComponentModeler(Tidy3dBaseModel): +class ComponentModeler(AbstractComponentModeler): """ Tool for modeling devices and computing scattering matrix elements. @@ -90,29 +39,13 @@ class ComponentModeler(Tidy3dBaseModel): * `Computing the scattering matrix of a device <../../notebooks/SMatrix.html>`_ """ - simulation: Simulation = pd.Field( - ..., - title="Simulation", - description="Simulation describing the device without any sources present.", - ) ports: Tuple[Port, ...] = pd.Field( (), title="Ports", description="Collection of ports describing the scattering matrix elements. " "For each input mode, one simulation will be run with a modal source.", ) - freqs: FreqArray = pd.Field( - ..., - title="Frequencies", - description="Array or list of frequencies at which to evaluate the scattering matrix.", - units=HERTZ, - ) - folder_name: str = pd.Field( - "default", - title="Folder Name", - description="Name of the folder for the tasks on web.", - ) element_mappings: Tuple[Tuple[Element, Element, Complex], ...] = pd.Field( (), title="Element Mappings", @@ -125,6 +58,7 @@ class ComponentModeler(Tidy3dBaseModel): " ``element_mappings``, the simulation corresponding to this column " "is skipped automatically.", ) + run_only: Optional[Tuple[MatrixIndex, ...]] = pd.Field( None, title="Run Only", @@ -142,6 +76,7 @@ class ComponentModeler(Tidy3dBaseModel): title="Verbosity", description="Whether the :class:`.ComponentModeler` should print status and progressbars.", ) + callback_url: str = pd.Field( None, title="Callback URL", @@ -149,11 +84,6 @@ class ComponentModeler(Tidy3dBaseModel): "The body content is a json file with fields " "``{'id', 'status', 'name', 'workUnit', 'solverVersion'}``.", ) - path_dir: str = pd.Field( - DEFAULT_DATA_DIR, - title="Directory Path", - description="Base directory where data and batch will be downloaded.", - ) @pd.validator("simulation", always=True) def _sim_has_no_sources(cls, val): @@ -223,12 +153,22 @@ def matrix_indices_run_sim(self) -> Tuple[MatrixIndex, ...]: return source_indices_needed - def get_port_by_name(self, port_name: str) -> Port: - """Get the port from the name.""" - ports = [port for port in self.ports if port.name == port_name] - if len(ports) == 0: - raise Tidy3dKeyError(f'Port "{port_name}" not found.') - return ports[0] + @cached_property + def port_names(self) -> Tuple[List[str], List[str]]: + """List of port names for inputs and outputs, respectively.""" + + def get_port_names(matrix_elements: Tuple[str, int]) -> List[str]: + """Get the port names from a list of (port name, mode index).""" + port_names = [] + for port_name, _ in matrix_elements: + if port_name not in port_names: + port_names.append(port_name) + return port_names + + port_names_in = get_port_names(self.matrix_indices_source) + port_names_out = get_port_names(self.matrix_indices_monitor) + + return port_names_out, port_names_in def to_monitor(self, port: Port) -> ModeMonitor: """Creates a mode monitor from a given port.""" @@ -303,16 +243,9 @@ def shift_port(self, port: Port) -> Port: port_shifted = port.copy(update=dict(center=center_shifted)) return port_shifted - @staticmethod - def _task_name(port: Port, mode_index: int) -> str: - """The name of a task, determined by the port of the source and mode index.""" - return f"smatrix_{port.name}_{mode_index}" - @equal_aspect @add_ax_if_none - def plot_sim( - self, x: float = None, y: float = None, z: float = None, ax: Ax = None, **kwargs - ) -> Ax: + def plot_sim(self, x: float = None, y: float = None, z: float = None, ax: Ax = None) -> Ax: """Plot a :class:`Simulation` with all sources added for each port, for troubleshooting.""" plot_sources = [] @@ -320,7 +253,7 @@ def plot_sim( mode_source_0 = self.to_source(port=port_source, mode_index=0) plot_sources.append(mode_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) + return sim_plot.plot(x=x, y=y, z=z, ax=ax) @equal_aspect @add_ax_if_none @@ -336,57 +269,6 @@ def plot_sim_eps( sim_plot = self.simulation.copy(update=dict(sources=plot_sources)) return sim_plot.plot_eps(x=x, y=y, z=z, ax=ax, **kwargs) - @cached_property - def batch(self) -> Batch: - """Batch associated with this component modeler.""" - - # first try loading the batch from file, if it exists - batch_path = self._batch_path - if os.path.exists(batch_path): - return Batch.from_file(fname=batch_path) - - return Batch( - simulations=self.sim_dict, - folder_name=self.folder_name, - callback_url=self.callback_url, - verbose=self.verbose, - ) - - @cached_property - def batch_path(self) -> str: - """Path to the batch saved to file.""" - - return self.batch._batch_path(path_dir=DEFAULT_DATA_DIR) - - def get_path_dir(self, path_dir: str) -> None: - """Check whether the supplied 'path_dir' matches the internal field value.""" - - if path_dir not in (DEFAULT_DATA_DIR, self.path_dir): - log.warning( - f"'ComponentModeler' method was supplied a 'path_dir' of '{path_dir}' " - f"when its internal 'path_dir' field was set to '{self.path_dir}'. " - "The passed value will be deprecated in later versions. " - "Please set the internal 'path_dir' field to the desired value and " - "remove the 'path_dir' from the method argument. " - f"Using supplied '{path_dir}'." - ) - return path_dir - - return self.path_dir - - @cached_property - def _batch_path(self) -> str: - """Where we store the batch for this ComponentModeler instance after the run.""" - hash_str = self._hash_self() - return os.path.join(self.path_dir, "batch" + hash_str + ".json") - - def _run_sims(self, path_dir: str = DEFAULT_DATA_DIR) -> BatchData: - """Run :class:`Simulations` for each port and return the batch after saving.""" - batch = self.batch - batch_data = batch.run(path_dir=path_dir) - batch.to_file(self._batch_path) - return batch_data - def _normalization_factor(self, port_source: Port, sim_data: SimulationData) -> complex: """Compute the normalization amplitude based on the measured input mode amplitude.""" @@ -414,24 +296,7 @@ def get_max_mode_indices(matrix_elements: Tuple[str, int]) -> int: return max_mode_index_out, max_mode_index_in - @cached_property - def port_names(self) -> Tuple[List[str], List[str]]: - """List of port names for inputs and outputs, respectively.""" - - def get_port_names(matrix_elements: Tuple[str, int]) -> List[str]: - """Get the port names from a list of (port name, mode index).""" - port_names = [] - for port_name, _ in matrix_elements: - if port_name not in port_names: - port_names.append(port_name) - return port_names - - port_names_in = get_port_names(self.matrix_indices_source) - port_names_out = get_port_names(self.matrix_indices_monitor) - - return port_names_out, port_names_in - - def _construct_smatrix(self, batch_data: BatchData) -> SMatrixDataArray: + def _construct_smatrix(self, batch_data: BatchData) -> ModalPortDataArray: """Post process `BatchData` to generate scattering matrix.""" max_mode_index_out, max_mode_index_in = self.max_mode_index @@ -450,7 +315,7 @@ def _construct_smatrix(self, batch_data: BatchData) -> SMatrixDataArray: mode_index_in=range(num_modes_in), f=np.array(self.freqs), ) - s_matrix = SMatrixDataArray(values, coords=coords) + s_matrix = ModalPortDataArray(values, coords=coords) # loop through source ports for col_index in self.matrix_indices_run_sim: @@ -502,17 +367,3 @@ def _construct_smatrix(self, batch_data: BatchData) -> SMatrixDataArray: s_matrix.loc[coords_to] = mult_by * s_matrix.loc[coords_from].values return s_matrix - - def run(self, path_dir: str = DEFAULT_DATA_DIR) -> SMatrixDataArray: - """Solves for the scattering matrix of the system.""" - path_dir = self.get_path_dir(path_dir) - - batch_data = self._run_sims(path_dir=path_dir) - return self._construct_smatrix(batch_data=batch_data) - - def load(self, path_dir: str = DEFAULT_DATA_DIR) -> SMatrixDataArray: - """Load a scattering matrix from saved :class:`BatchData` object.""" - path_dir = self.get_path_dir(path_dir) - - batch_data = BatchData.load(path_dir=path_dir) - return self._construct_smatrix(batch_data=batch_data) diff --git a/tidy3d/plugins/smatrix/component_modelers/terminal.py b/tidy3d/plugins/smatrix/component_modelers/terminal.py new file mode 100644 index 000000000..eae6302be --- /dev/null +++ b/tidy3d/plugins/smatrix/component_modelers/terminal.py @@ -0,0 +1,335 @@ +"""Tool for generating an S matrix automatically from a Tidy3d simulation and lumped port definitions.""" +from __future__ import annotations + +from typing import Tuple, Dict + +import pydantic.v1 as pd +import numpy as np +import xarray as xr + +from ....constants import C_0, fp_eps +from ....components.simulation import Simulation +from ....components.geometry.utils_2d import snap_coordinate_to_grid +from ....components.data.sim_data import SimulationData +from ....components.source import GaussianPulse +from ....components.types import Ax +from ....components.viz import add_ax_if_none, equal_aspect +from ....components.base import cached_property +from ....exceptions import ValidationError +from ....web.api.container import BatchData + +from .base import AbstractComponentModeler, FWIDTH_FRAC +from ..ports.lumped import LumpedPortDataArray, LumpedPort + + +class TerminalComponentModeler(AbstractComponentModeler): + """Tool for modeling two-terminal multiport devices and computing port parameters + with lumped ports.""" + + ports: Tuple[LumpedPort, ...] = 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.", + ) + + @equal_aspect + @add_ax_if_none + def plot_sim( + self, x: float = None, y: float = None, z: float = None, ax: Ax = None, **kwargs + ) -> Ax: + """Plot a :class:`Simulation` with all sources added for each port, for troubleshooting.""" + + plot_sources = [] + for port_source in self.ports: + source_0 = port_source.to_source(self._source_time, snap_center=None) + 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) + + @equal_aspect + @add_ax_if_none + def plot_sim_eps( + self, x: float = None, y: float = None, z: float = None, ax: Ax = None, **kwargs + ) -> Ax: + """Plot permittivity of the :class:`Simulation` with all sources added for each port.""" + + plot_sources = [] + for port_source in self.ports: + source_0 = port_source.to_source(self._source_time, snap_center=None) + 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) + + @cached_property + def sim_dict(self) -> Dict[str, Simulation]: + """Generate all the :class:`Simulation` objects for the port parameter calculation.""" + + 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] + + # 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 + # this surface, so the mesh override needs to ensure that the mesh + # is fine enough not only in plane, but also in the normal direction. + # So in the normal direction, we'll make sure there are at least + # 2 cell layers above and below whose size is the same as the in-plane + # cell size in the override region. Also, to ensure that the port itself + # is aligned with a grid boundary in the normal direction, two separate + # 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): + 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( + update={ + "wavelength": C_0 / np.max(self.freqs), + "override_structures": list(self.simulation.grid_spec.override_structures) + + mesh_overrides, + } + ) + + # 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: + port_source = port.to_source(self._source_time, snap_center=None) + update_dict = dict( + sources=[port_source], + monitors=new_mnts, + lumped_elements=lumped_resistors, + grid_spec=grid_spec, + ) + + sim_copy = self.simulation.copy(update=update_dict) + task_name = self._task_name(port=port) + sim_dict[task_name] = sim_copy + + # Check if snapping to grid is needed + 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 + ) + 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 + ) + for port, center in zip(self.ports, snap_centers): + port_source = self.to_source(self._source_time, port=port, snap_center=center) + update_dict = dict( + sources=[port_source], + monitors=new_mnts, + lumped_elements=lumped_resistors, + grid_spec=grid_spec, + ) + + sim_copy = self.simulation.copy(update=update_dict) + task_name = self._task_name(port=port) + sim_dict[task_name] = sim_copy + + return sim_dict + + @cached_property + def _source_time(self): + """Helper to create a time domain pulse for the frequeny range of interest.""" + freq0 = np.mean(self.freqs) + fdiff = max(self.freqs) - min(self.freqs) + fwidth = max(fdiff, freq0 * FWIDTH_FRAC) + return GaussianPulse( + freq0=freq0, fwidth=fwidth, remove_dc_component=self.remove_dc_component + ) + + def _construct_smatrix(self, batch_data: BatchData) -> LumpedPortDataArray: + """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)), + dtype=complex, + ) + coords = dict( + port_out=port_names, + port_in=port_names, + f=np.array(self.freqs), + ) + a_matrix = LumpedPortDataArray(values, coords=coords) + b_matrix = a_matrix.copy(deep=True) + + def port_voltage(port: LumpedPort, sim_data: SimulationData) -> xr.DataArray: + """Helper to compute voltage across the port.""" + e_component = "xyz"[port.voltage_axis] + field_data = sim_data[f"{port.name}_E{e_component}"] + e_field = field_data.field_components[f"E{e_component}"] + e_coords = [e_field.x, e_field.y, e_field.z] + # interpolate E locations to coincide with port bounds along the integration path + e_coords_interp = { + e_component: np.linspace( + port.bounds[0][port.voltage_axis], + port.bounds[1][port.voltage_axis], + len(e_coords[port.voltage_axis]), + ) + } + e_field = e_field.interp(**e_coords_interp) + voltage = -e_field.integrate(coord=e_component).squeeze(drop=True) + # Return data array of voltage with coordinates of frequency + return voltage + + def port_current(port: LumpedPort, sim_data: SimulationData) -> xr.DataArray: + """Helper to compute current flowing through the port.""" + + # Diagram of contour integral, dashed line indicates location of sheet resistance + # and electric field used for voltage computation. Voltage axis is out-of-page. + # + # current_axis = -> + # injection_axis = ^ + # + # | h2_field -> | + # h_cap_minus ^ ------------------------------------------- h_cap_plus ^ + # | h1_field -> | + + # Get h field tangent to resistive sheet + h_component = "xyz"[port.current_axis] + orth_component = "xyz"[port.injection_axis] + field_data = sim_data[f"{port.name}_H{h_component}"] + h_field = field_data.field_components[f"H{h_component}"] + h_coords = [h_field.x, h_field.y, h_field.z] + # h_cap represents the very short section (single cell) of the H contour that + # is in the injection_axis direction. It is needed to fully enclose the sheet. + h_cap_field = field_data.field_components[f"H{orth_component}"] + h_cap_coords = [h_cap_field.x, h_cap_field.y, h_cap_field.z] + + # Use the coordinates of h_cap since it lies on the same grid that the + # lumped resistor is snapped to + orth_index = np.argmin( + np.abs(h_cap_coords[port.injection_axis].values - port.center[port.injection_axis]) + ) + # Some sanity checks, tangent H field coordinates should be directly above + # and below the coordinates of the resistive sheet + assert orth_index > 0 + assert ( + h_cap_coords[port.injection_axis].values[orth_index] + < h_coords[port.injection_axis].values[orth_index] + ) + assert ( + h_coords[port.injection_axis].values[orth_index - 1] + < h_cap_coords[port.injection_axis].values[orth_index] + ) + + # Extract field just below and just above sheet + h1_field = h_field.isel({orth_component: orth_index - 1}) + h2_field = h_field.isel({orth_component: orth_index}) + h_field = h1_field - h2_field + # Extract cap field which is coincident with sheet + h_cap = h_cap_field.isel({orth_component: orth_index}) + + # Need to make sure to use the nearest coordinate that is + # at least greater than the port bounds + hcap_minus = h_cap.sel({h_component: slice(-np.inf, port.bounds[0][port.current_axis])}) + hcap_plus = h_cap.sel({h_component: slice(port.bounds[1][port.current_axis], np.inf)}) + hcap_minus = hcap_minus.isel({h_component: -1}) + hcap_plus = hcap_plus.isel({h_component: 0}) + # Length of integration along the h_cap contour is a single cell width + dcap = ( + h_coords[port.injection_axis].values[orth_index] + - h_coords[port.injection_axis].values[orth_index - 1] + ) + + h_min_bound = hcap_minus.coords[h_component].values + h_max_bound = hcap_plus.coords[h_component].values + h_coords_interp = { + h_component: np.linspace( + h_min_bound, + h_max_bound, + len(h_coords[port.current_axis] + 2), + ) + } + # Integration that corresponds to the tangent H field + h_field = h_field.interp(**h_coords_interp) + current = h_field.integrate(coord=h_component).squeeze(drop=True) + + # Integration that corresponds with the contribution to current from cap contours + hcap_current = ( + ((hcap_plus - hcap_minus) * dcap).squeeze(drop=True).reset_coords(drop=True) + ) + # Add the contribution from the hcap integral + current = current + hcap_current + # Make sure we compute current flowing from plus to minus voltage + if port.current_axis != (port.voltage_axis + 1) % 3: + current *= -1 + # Return data array of current with coordinates of frequency + return current + + def port_ab(port: LumpedPort, sim_data: SimulationData): + """Helper to compute the port incident and reflected power waves.""" + voltage = port_voltage(port, sim_data) + current = port_current(port, sim_data) + 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 + + # 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[ + dict( + port_in=port_in.name, + port_out=port_out.name, + ) + ] = a_out + + b_matrix.loc[ + dict( + port_in=port_in.name, + port_out=port_out.name, + ) + ] = b_out + + s_matrix = self.ab_to_s(a_matrix, b_matrix) + + return s_matrix + + @pd.validator("simulation") + def _validate_3d_simulation(cls, val): + """Error if Simulation is not a 3D simulation""" + + if val.size.count(0.0) > 0: + raise ValidationError( + f"'{cls.__name__}' must be setup with a 3D simulation with all sizes greater than 0." + ) + return val diff --git a/tidy3d/plugins/smatrix/ports/lumped.py b/tidy3d/plugins/smatrix/ports/lumped.py new file mode 100644 index 000000000..961e478fd --- /dev/null +++ b/tidy3d/plugins/smatrix/ports/lumped.py @@ -0,0 +1,173 @@ +"""Class and custom data array for representing a scattering matrix port based on lumped circuit elements.""" +import pydantic.v1 as pd +import numpy as np +from typing import Optional + +from ....constants import OHM +from ....components.geometry.base import Box +from ....components.types import Complex, FreqArray, Axis +from ....components.base import cached_property +from ....components.lumped_element import LumpedResistor +from ....components.monitor import FieldMonitor +from ....components.source import UniformCurrentSource, GaussianPulse +from ....components.validators import assert_plane +from ....components.data.data_array import DataArray +from ....exceptions import ValidationError + +DEFAULT_PORT_NUM_CELLS = 3 + + +class LumpedPortDataArray(DataArray): + """Port parameter matrix elements for lumped ports. + + Example + ------- + >>> port_in = ['port1', 'port2'] + >>> port_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 LumpedPort(Box): + """Class representing a single lumped port""" + + name: str = pd.Field( + ..., + title="Name", + description="Unique name for the port.", + min_length=1, + ) + + voltage_axis: Axis = pd.Field( + ..., + title="Voltage Integration Axis", + description="Specifies the axis along which the E-field line integral is performed when " + "computing the port voltage. The integration axis must lie in the plane of the port.", + ) + + impedance: Complex = pd.Field( + 50, + title="Reference impedance", + description="Reference port impedance for scattering parameter computation.", + units=OHM, + ) + + num_grid_cells: Optional[pd.PositiveInt] = pd.Field( + DEFAULT_PORT_NUM_CELLS, + title="Port grid cells", + description="Number of mesh grid cells associated with the port along each direction, " + "which are added through automatic mesh refinement. " + "A value of ``None`` will turn off automatic mesh refinement.", + ) + + _plane_validator = assert_plane() + + @cached_property + def injection_axis(self): + """Injection axis of the port.""" + return self.size.index(0.0) + + @pd.validator("voltage_axis", always=True) + def _voltage_axis_in_plane(cls, val, values): + """Ensure voltage integration axis is in the port's plane.""" + size = values.get("size") + if val == size.index(0.0): + raise ValidationError("'voltage_axis' must lie in the port's plane.") + return val + + @cached_property + def current_axis(self) -> Axis: + """Integration axis for computing the port current via the magnetic field.""" + return 3 - self.injection_axis - self.voltage_axis + + def to_source( + self, + source_time: GaussianPulse, + snap_center: float, + ) -> UniformCurrentSource: + """Create a current source from the lumped port.""" + # Discretized source amps are manually zeroed out later if they + # fall on Yee grid locations outside the analytical source region. + component = "xyz"[self.voltage_axis] + center = list(self.center) + if snap_center: + center[self.injection_axis] = snap_center + return UniformCurrentSource( + center=center, + size=self.size, + source_time=source_time, + polarization=f"E{component}", + name=self.name, + interpolate=True, + confine_to_bounds=True, + ) + + def to_load(self, snap_center: float) -> 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 + center = list(self.center) + if snap_center: + center[self.injection_axis] = snap_center + return LumpedResistor( + center=center, + size=self.size, + num_grid_cells=self.num_grid_cells, + resistance=np.real(self.impedance), + name=f"{self.name}_resistor", + voltage_axis=self.voltage_axis, + ) + + def to_voltage_monitor(self, freqs: FreqArray, snap_center: float) -> FieldMonitor: + """Field monitor to compute port voltage.""" + center = list(self.center) + if snap_center: + center[self.injection_axis] = snap_center + + e_component = "xyz"[self.voltage_axis] + # Size of voltage monitor can essentially be 1D from ground to signal conductor + voltage_mon_size = list(self.size) + voltage_mon_size[self.injection_axis] = 0.0 + voltage_mon_size[self.current_axis] = 0.0 + # Create a voltage monitor + return FieldMonitor( + center=center, + size=voltage_mon_size, + freqs=freqs, + fields=[f"E{e_component}"], + name=f"{self.name}_E{e_component}", + colocate=False, + ) + + def to_current_monitor(self, freqs: FreqArray, snap_center: float) -> FieldMonitor: + """Field monitor to compute port current.""" + center = list(self.center) + if snap_center: + center[self.injection_axis] = snap_center + + h_component = "xyz"[self.current_axis] + h_cap_component = "xyz"[self.injection_axis] + # Size of current monitor needs to encompass the current carrying 2D sheet + # Needs to have a nonzero thickness so a closed loop of gridpoints around the 2D sheet can be formed + current_mon_size = list(self.size) + current_mon_size[self.injection_axis] = 1e-10 + current_mon_size[self.voltage_axis] = 0.0 + # Create a current monitor + return FieldMonitor( + center=center, + size=current_mon_size, + freqs=freqs, + fields=[f"H{h_component}", f"H{h_cap_component}"], + name=f"{self.name}_H{h_component}", + colocate=False, + ) diff --git a/tidy3d/plugins/smatrix/ports/modal.py b/tidy3d/plugins/smatrix/ports/modal.py new file mode 100644 index 000000000..1e51427ea --- /dev/null +++ b/tidy3d/plugins/smatrix/ports/modal.py @@ -0,0 +1,53 @@ +"""Class and custom data array for representing a scattering matrix port based on waveguide modes.""" +import pydantic.v1 as pd + +from ....components.geometry.base import Box +from ....components.mode import ModeSpec +from ....components.types import Direction +from ....components.data.data_array import DataArray + + +class ModalPortDataArray(DataArray): + """Port parameter matrix elements for modal ports. + + Example + ------- + >>> port_in = ['port1', 'port2'] + >>> port_out = ['port1', 'port2'] + >>> mode_index_in = [0, 1] + >>> mode_index_out = [0, 1] + >>> f = [2e14] + >>> coords = dict( + ... port_in=ports_in, + ... port_out=ports_out, + ... mode_index_in=mode_index_in, + ... mode_index_out=mode_index_out, + ... f=f + ... ) + >>> fd = ModalPortDataArray((1 + 1j) * np.random.random((2, 2, 2, 2, 1)), coords=coords) + """ + + __slots__ = () + _dims = ("port_out", "mode_index_out", "port_in", "mode_index_in", "f") + _data_attrs = {"long_name": "modal port matrix element"} + + +class Port(Box): + """Specifies a port in the scattering matrix.""" + + direction: Direction = pd.Field( + ..., + title="Direction", + description="'+' or '-', defining which direction is considered 'input'.", + ) + mode_spec: ModeSpec = pd.Field( + ModeSpec(), + title="Mode Specification", + description="Specifies how the mode solver will solve for the modes of the port.", + ) + name: str = pd.Field( + ..., + title="Name", + description="Unique name for the port.", + min_length=1, + ) diff --git a/tidy3d/web/api/container.py b/tidy3d/web/api/container.py index f0d087754..56a06bf1f 100644 --- a/tidy3d/web/api/container.py +++ b/tidy3d/web/api/container.py @@ -307,7 +307,7 @@ def estimate_cost(self, verbose: bool = True) -> float: Cost is calculated assuming the simulation runs for the full ``run_time``. If early shut-off is triggered, the cost is adjusted proportionately. """ - return web.estimate_cost(self.task_id, verbose=verbose) + return web.estimate_cost(self.task_id, verbose=verbose, solver_version=self.solver_version) class BatchData(Tidy3dBaseModel): diff --git a/tidy3d/web/api/webapi.py b/tidy3d/web/api/webapi.py index fcd1e4f00..aae3d73e9 100644 --- a/tidy3d/web/api/webapi.py +++ b/tidy3d/web/api/webapi.py @@ -777,7 +777,7 @@ def get_tasks( @wait_for_connection -def estimate_cost(task_id: str, verbose: bool = True) -> float: +def estimate_cost(task_id: str, verbose: bool = True, solver_version: str = None) -> float: """Compute the maximum FlexCredit charge for a given task. Parameters @@ -786,6 +786,8 @@ def estimate_cost(task_id: str, verbose: bool = True) -> float: Unique identifier of task on server. Returned by :meth:`upload`. verbose : bool = True Whether to log the cost and helpful messages. + solver_version : str = None + Target solver version. Returns ------- @@ -825,7 +827,7 @@ def estimate_cost(task_id: str, verbose: bool = True) -> float: if not task: raise ValueError("Task not found.") - task.estimate_cost() + task.estimate_cost(solver_version=solver_version) task_info = get_info(task_id) status = task_info.metadataStatus