diff --git a/CHANGELOG.md b/CHANGELOG.md index f094397115..b5f58ed17f 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -10,6 +10,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### Added - Added `symmetrize_mirror`, `symmetrize_rotation`, `symmetrize_diagonal` functions to the autograd plugin. They can be used for enforcing symmetries in topology optimization. - Klayout plugin automatically finds klayout installation path at common locations. +- Added `EMESimulationData.coeffs` to store coefficients from the EME solver, including mode overlaps, interface S matrices, and effective propagation indices. ### Changed - Removed validator that would warn if `PerturbationMedium` values could become numerically unstable, since an error will anyway be raised if this actually happens when the medium is converted using actual perturbation data. diff --git a/schemas/EMESimulation.json b/schemas/EMESimulation.json index 42e7a1ac19..df704ebeeb 100644 --- a/schemas/EMESimulation.json +++ b/schemas/EMESimulation.json @@ -12879,6 +12879,10 @@ }, "type": "array" }, + "store_coeffs": { + "default": true, + "type": "boolean" + }, "store_port_modes": { "default": true, "type": "boolean" diff --git a/tests/sims/full_fdtd.h5 b/tests/sims/full_fdtd.h5 index 4f383a0e73..c53806fe81 100644 Binary files a/tests/sims/full_fdtd.h5 and b/tests/sims/full_fdtd.h5 differ diff --git a/tests/test_components/test_eme.py b/tests/test_components/test_eme.py index 058ae384ac..3910b3a3e4 100644 --- a/tests/test_components/test_eme.py +++ b/tests/test_components/test_eme.py @@ -311,16 +311,16 @@ def test_eme_simulation(): # test warning for not providing wavelength in autogrid grid_spec = td.GridSpec.auto(min_steps_per_wvl=20) + sim = sim.updated_copy(grid_spec=grid_spec) with AssertLogLevel("INFO", contains_str="wavelength"): - sim = sim.updated_copy(grid_spec=grid_spec) + _ = sim.updated_copy(monitors=[]) # multiple freqs are ok, but not for autogrid _ = sim.updated_copy( grid_spec=td.GridSpec.uniform(dl=0.2), freqs=[10000000000.0, *list(sim.freqs)] ) with AssertLogLevel("INFO", contains_str="wavelength"): _ = sim.updated_copy( - freqs=[*list(sim.freqs), 10000000000.0], - grid_spec=grid_spec, + freqs=[*list(sim.freqs), 10000000000.0], grid_spec=grid_spec, monitors=[] ) # test port offsets @@ -558,12 +558,14 @@ def test_eme_simulation(): constraint="passive", eme_grid_spec=td.EMEUniformGrid(num_cells=1, mode_spec=td.EMEModeSpec(num_modes=40)), grid_spec=sim.grid_spec.updated_copy(wavelength=1), + monitors=[], ) sim_good.validate_pre_upload() sim_good = sim.updated_copy( constraint=None, eme_grid_spec=td.EMEUniformGrid(num_cells=1, mode_spec=td.EMEModeSpec(num_modes=60)), grid_spec=sim.grid_spec.updated_copy(wavelength=1), + monitors=[], ) sim_good.validate_pre_upload() # warn about num modes with constraint @@ -731,6 +733,47 @@ def _get_eme_smatrix_data_array(num_modes_in=2, num_modes_out=3, num_freqs=2, nu return smatrix_entry +def _get_eme_interface_smatrix_data_array( + num_modes_in=2, num_modes_out=3, num_freqs=2, num_sweep=0 +): + if num_modes_in != 0: + mode_index_in = np.arange(num_modes_in) + else: + mode_index_in = [0] + if num_modes_out != 0: + mode_index_out = np.arange(num_modes_out) + else: + mode_index_out = [0] + if num_sweep != 0: + sweep_index = np.arange(num_sweep) + else: + sweep_index = [0] + eme_cell_index = np.arange(3) + + f = td.C_0 * np.linspace(1, 2, num_freqs) + + data = (1 + 1j) * np.random.random( + (len(f), len(sweep_index), len(eme_cell_index), len(mode_index_out), len(mode_index_in)) + ) + coords = { + "f": f, + "sweep_index": sweep_index, + "eme_cell_index": eme_cell_index, + "mode_index_out": mode_index_out, + "mode_index_in": mode_index_in, + } + smatrix_entry = td.EMEInterfaceSMatrixDataArray(data, coords=coords) + + if num_modes_in == 0: + smatrix_entry = smatrix_entry.drop_vars("mode_index_in") + if num_modes_out == 0: + smatrix_entry = smatrix_entry.drop_vars("mode_index_out") + if num_sweep == 0: + smatrix_entry = smatrix_entry.drop_vars("sweep_index") + + return smatrix_entry + + def _get_eme_smatrix_dataset(num_modes_1=3, num_modes_2=4, num_sweep=0): S11 = _get_eme_smatrix_data_array( num_modes_in=num_modes_1, num_modes_out=num_modes_1, num_sweep=num_sweep @@ -747,6 +790,35 @@ def _get_eme_smatrix_dataset(num_modes_1=3, num_modes_2=4, num_sweep=0): return td.EMESMatrixDataset(S11=S11, S12=S12, S21=S21, S22=S22) +def _get_eme_interface_smatrix_dataset(num_modes_1=3, num_modes_2=4, num_sweep=0): + S11 = _get_eme_interface_smatrix_data_array( + num_modes_in=num_modes_1, num_modes_out=num_modes_1, num_sweep=num_sweep + ) + S12 = _get_eme_interface_smatrix_data_array( + num_modes_in=num_modes_2, num_modes_out=num_modes_1, num_sweep=num_sweep + ) + S21 = _get_eme_interface_smatrix_data_array( + num_modes_in=num_modes_1, num_modes_out=num_modes_2, num_sweep=num_sweep + ) + S22 = _get_eme_interface_smatrix_data_array( + num_modes_in=num_modes_2, num_modes_out=num_modes_2, num_sweep=num_sweep + ) + return td.EMEInterfaceSMatrixDataset(S11=S11, S12=S12, S21=S21, S22=S22) + + +def _get_eme_overlaps_dataset(num_modes_1=3, num_modes_2=4, num_sweep=0): + O11 = _get_eme_interface_smatrix_data_array( + num_modes_in=num_modes_1, num_modes_out=num_modes_1, num_sweep=num_sweep + ) + O12 = _get_eme_interface_smatrix_data_array( + num_modes_in=num_modes_2, num_modes_out=num_modes_1, num_sweep=num_sweep + ) + O21 = _get_eme_interface_smatrix_data_array( + num_modes_in=num_modes_1, num_modes_out=num_modes_2, num_sweep=num_sweep + ) + return td.EMEOverlapDataset(O11=O11, O12=O12, O21=O21) + + def _get_eme_coeff_data_array(num_sweep=0): f = [2e14] mode_index_out = [0, 1] @@ -787,7 +859,26 @@ def _get_eme_coeff_data_array(num_sweep=0): def _get_eme_coeff_dataset(num_sweep=0): A = _get_eme_coeff_data_array(num_sweep=num_sweep) B = _get_eme_coeff_data_array(num_sweep=num_sweep) - return td.EMECoefficientDataset(A=A, B=B) + flux = _get_eme_flux_data_array(num_sweep=num_sweep) + n_complex = _get_eme_mode_index_data_array(num_sweep=num_sweep) + interface_smatrices = _get_eme_interface_smatrix_dataset(num_sweep=num_sweep) + overlaps = _get_eme_overlaps_dataset(num_sweep=num_sweep) + return td.EMECoefficientDataset( + A=A, + B=B, + flux=flux, + n_complex=n_complex, + interface_smatrices=interface_smatrices, + overlaps=overlaps, + ) + + +def test_eme_normalize_coeff_dataset(): + coeffs = _get_eme_coeff_dataset() + coeffs_normalized = coeffs.normalized_copy + assert coeffs_normalized.flux is None + with pytest.raises(ValidationError): + _ = coeffs_normalized.normalized_copy def test_eme_coeff_data_array(): @@ -819,6 +910,29 @@ def _get_eme_mode_index_data_array(num_sweep=0): return data +def _get_eme_flux_data_array(num_sweep=0): + f = [td.C_0, 3e14] + mode_index = np.arange(10) + eme_cell_index = np.arange(7) + if num_sweep != 0: + sweep_index = np.arange(num_sweep) + else: + sweep_index = [0] + coords = { + "f": f, + "sweep_index": sweep_index, + "eme_cell_index": eme_cell_index, + "mode_index": mode_index, + } + data = td.EMEFluxDataArray( + np.random.random((len(f), len(sweep_index), len(eme_cell_index), len(mode_index))), + coords=coords, + ) + if num_sweep == 0: + data = data.drop_vars("sweep_index") + return data + + def test_eme_mode_index_data_array(): _ = _get_eme_mode_index_data_array() @@ -1305,7 +1419,7 @@ def test_eme_periodicity(): # remove the field monitor, now it passes desired_cell_index_pairs = set([(i, i + 1) for i in range(6)] + [(5, 1)]) - with AssertLogLevel(None): + with AssertLogLevel("WARNING", contains_str="deprecated"): sim = sim.updated_copy( monitors=[m for m in sim.monitors if not isinstance(m, td.EMEFieldMonitor)] ) diff --git a/tidy3d/__init__.py b/tidy3d/__init__.py index 18f7d99047..13934377f7 100644 --- a/tidy3d/__init__.py +++ b/tidy3d/__init__.py @@ -185,6 +185,8 @@ ChargeDataArray, DiffractionDataArray, EMECoefficientDataArray, + EMEFluxDataArray, + EMEInterfaceSMatrixDataArray, EMEModeIndexDataArray, EMEScalarFieldDataArray, EMEScalarModeFieldDataArray, @@ -243,7 +245,9 @@ from .components.eme.data.dataset import ( EMECoefficientDataset, EMEFieldDataset, + EMEInterfaceSMatrixDataset, EMEModeSolverDataset, + EMEOverlapDataset, EMESMatrixDataset, ) from .components.eme.data.monitor_data import EMECoefficientData, EMEFieldData, EMEModeSolverData @@ -605,8 +609,11 @@ def set_logging_level(level: str) -> None: "EMEFieldData", "EMEFieldDataset", "EMEFieldMonitor", + "EMEFluxDataArray", "EMEFreqSweep", "EMEGrid", + "EMEInterfaceSMatrixDataArray", + "EMEInterfaceSMatrixDataset", "EMELengthSweep", "EMEModeIndexDataArray", "EMEModeSolverData", @@ -615,6 +622,7 @@ def set_logging_level(level: str) -> None: "EMEModeSpec", "EMEModeSweep", "EMEMonitor", + "EMEOverlapDataset", "EMEPeriodicitySweep", "EMESMatrixDataArray", "EMESMatrixDataset", diff --git a/tidy3d/components/data/data_array.py b/tidy3d/components/data/data_array.py index ea6385a796..c6d6d9e38a 100644 --- a/tidy3d/components/data/data_array.py +++ b/tidy3d/components/data/data_array.py @@ -1188,6 +1188,31 @@ class EMESMatrixDataArray(DataArray): _data_attrs = {"long_name": "scattering matrix element"} +class EMEInterfaceSMatrixDataArray(DataArray): + """Scattering matrix elements at a single cell interface for a fixed pair of ports, + possibly with an extra sweep index. + Example + ------- + >>> mode_index_in = [0, 1] + >>> mode_index_out = [0, 1, 2] + >>> eme_cell_index = [2, 4] + >>> f = [2e14] + >>> sweep_index = np.arange(10) + >>> coords = dict( + ... f=f, + ... sweep_index=sweep_index, + ... eme_cell_index=eme_cell_index, + ... mode_index_out=mode_index_out, + ... mode_index_in=mode_index_in, + ... ) + >>> fd = EMEInterfaceSMatrixDataArray((1 + 1j) * np.random.random((1, 10, 2, 3, 2)), coords=coords) + """ + + __slots__ = () + _dims = ("f", "sweep_index", "eme_cell_index", "mode_index_out", "mode_index_in") + _data_attrs = {"long_name": "scattering matrix element"} + + class EMEModeIndexDataArray(DataArray): """Complex-valued effective propagation index of an EME mode, also indexed by EME cell. @@ -1206,6 +1231,24 @@ class EMEModeIndexDataArray(DataArray): _data_attrs = {"long_name": "Propagation index"} +class EMEFluxDataArray(DataArray): + """Power flux of an EME mode, also indexed by EME cell. + + Example + ------- + >>> f = [2e14, 3e14] + >>> sweep_index = np.arange(2) + >>> eme_cell_index = np.arange(5) + >>> mode_index = np.arange(4) + >>> coords = dict(f=f, sweep_index=sweep_index, eme_cell_index=eme_cell_index, mode_index=mode_index) + >>> data = EMEFluxDataArray(np.random.random((2,2,5,4)), coords=coords) + """ + + __slots__ = () + _dims = ("f", "sweep_index", "eme_cell_index", "mode_index") + _data_attrs = {"units": WATT, "long_name": "flux"} + + class ChargeDataArray(DataArray): """Charge data array. @@ -1595,8 +1638,10 @@ def _make_impedance_data_array(result: DataArray) -> ImpedanceResultType: EMEScalarFieldDataArray, EMEScalarModeFieldDataArray, EMESMatrixDataArray, + EMEInterfaceSMatrixDataArray, EMECoefficientDataArray, EMEModeIndexDataArray, + EMEFluxDataArray, EMEFreqModeDataArray, ChargeDataArray, SteadyVoltageDataArray, diff --git a/tidy3d/components/eme/data/dataset.py b/tidy3d/components/eme/data/dataset.py index 33f45cef2c..3b3e7a3722 100644 --- a/tidy3d/components/eme/data/dataset.py +++ b/tidy3d/components/eme/data/dataset.py @@ -2,16 +2,23 @@ from __future__ import annotations +from typing import Optional + +import numpy as np import pydantic.v1 as pd +from tidy3d.components.base import cached_property from tidy3d.components.data.data_array import ( EMECoefficientDataArray, + EMEFluxDataArray, + EMEInterfaceSMatrixDataArray, EMEModeIndexDataArray, EMEScalarFieldDataArray, EMEScalarModeFieldDataArray, EMESMatrixDataArray, ) from tidy3d.components.data.dataset import Dataset, ElectromagneticFieldDataset +from tidy3d.exceptions import ValidationError class EMESMatrixDataset(Dataset): @@ -39,22 +46,143 @@ class EMESMatrixDataset(Dataset): ) +class EMEInterfaceSMatrixDataset(Dataset): + """Dataset storing S matrices associated with EME cell interfaces.""" + + S11: EMEInterfaceSMatrixDataArray = pd.Field( + ..., + title="S11 matrix", + description="S matrix relating output modes at port 1 to input modes at port 1.", + ) + S12: EMEInterfaceSMatrixDataArray = pd.Field( + ..., + title="S12 matrix", + description="S matrix relating output modes at port 1 to input modes at port 2.", + ) + S21: EMEInterfaceSMatrixDataArray = pd.Field( + ..., + title="S21 matrix", + description="S matrix relating output modes at port 2 to input modes at port 1.", + ) + S22: EMEInterfaceSMatrixDataArray = pd.Field( + ..., + title="S22 matrix", + description="S matrix relating output modes at port 2 to input modes at port 2.", + ) + + +class EMEOverlapDataset(Dataset): + """Dataset storing overlaps between EME modes. + + ``Oij`` is the unconjugated overlap computed using the E field of cell ``i`` + and the H field of cell ``j``. + + For consistency with ``Sij``, ``mode_index_out`` refers to the mode index + in cell ``i``, and ``mode_index_in`` refers to the mode index in cell ``j``. + """ + + O11: EMEInterfaceSMatrixDataArray = pd.Field( + ..., + title="O11 matrix", + description="Overlap integral between E field and H field in the same cell.", + ) + O12: EMEInterfaceSMatrixDataArray = pd.Field( + ..., + title="O12 matrix", + description="Overlap integral between E field on side 1 and H field on side 2.", + ) + O21: EMEInterfaceSMatrixDataArray = pd.Field( + ..., + title="O21 matrix", + description="Overlap integral between E field on side 2 and H field on side 1.", + ) + + class EMECoefficientDataset(Dataset): - """Dataset storing expansion coefficients for the modes in a cell. + """Dataset storing various coefficients related to the EME simulation. + These coefficients can be used for debugging or optimization. + + The ``A`` and ``B`` fields store the expansion coefficients for the modes in a cell. These are defined at the cell centers. + + The ``n_complex`` and ``flux`` fields respectively store the complex-valued effective + propagation index and the power flux associated with the modes used in the + EME calculation. + + The ``interface_Sij`` fields store the S matrices associated with the interfaces + between EME cells. """ - A: EMECoefficientDataArray = pd.Field( - ..., + A: Optional[EMECoefficientDataArray] = pd.Field( + None, title="A coefficient", description="Coefficient for forward mode in this cell.", ) - B: EMECoefficientDataArray = pd.Field( - ..., + + B: Optional[EMECoefficientDataArray] = pd.Field( + None, title="B coefficient", description="Coefficient for backward mode in this cell.", ) + n_complex: Optional[EMEModeIndexDataArray] = pd.Field( + None, + title="Propagation Index", + description="Complex-valued effective propagation indices associated with the EME modes.", + ) + + flux: Optional[EMEFluxDataArray] = pd.Field( + None, + title="Flux", + description="Power flux of the EME modes.", + ) + + interface_smatrices: Optional[EMEInterfaceSMatrixDataset] = pd.Field( + None, + title="Interface S Matrices", + description="S matrices associated with the interfaces between EME cells.", + ) + + overlaps: Optional[EMEOverlapDataset] = pd.Field( + None, title="Overlaps", description="Overlaps between EME modes." + ) + + @cached_property + def normalized_copy(self) -> EMECoefficientDataset: + """Return a copy of the ``EMECoefficientDataset`` where + the ``A`` and ``B`` coefficients as well as the ``interface_smatrices`` + are normalized by flux.""" + if self.flux is None: + raise ValidationError( + "The 'flux' field of the 'EMECoefficientDataset' is 'None', " + "so normalization cannot be performed." + ) + fields = {"A": self.A, "B": self.B} + flux_out = self.flux.rename(mode_index="mode_index_out") + for key, field in fields.items(): + if field is not None: + fields[key] = field * np.sqrt(flux_out) + if self.interface_smatrices is not None: + num_cells = len(self.flux.eme_cell_index) + flux1 = self.flux.isel(eme_cell_index=np.arange(num_cells - 1)) + flux2 = self.flux.isel(eme_cell_index=np.arange(1, num_cells)) + flux2 = flux2.assign_coords(eme_cell_index=np.arange(num_cells - 1)) + interface_S12 = self.interface_smatrices.S12 + flux_out = flux1.rename(mode_index="mode_index_out") + flux_in = flux2.rename(mode_index="mode_index_in") + interface_S12 = interface_S12 * np.sqrt(flux_out / flux_in) + interface_S21 = self.interface_smatrices.S21 + flux_out = flux2.rename(mode_index="mode_index_out") + flux_in = flux1.rename(mode_index="mode_index_in") + interface_S21 = interface_S21 * np.sqrt(flux_out / flux_in) + + fields["interface_smatrices"] = self.interface_smatrices.updated_copy( + S12=interface_S12, S21=interface_S21 + ) + # for safety to prevent normalizing twice + fields["flux"] = None + return self.updated_copy(**fields) + class EMEFieldDataset(ElectromagneticFieldDataset): """Dataset storing scalar components of E and H fields as a function of freq, mode_index, and port_index.""" diff --git a/tidy3d/components/eme/data/sim_data.py b/tidy3d/components/eme/data/sim_data.py index 05c1e6e739..750c088f2a 100644 --- a/tidy3d/components/eme/data/sim_data.py +++ b/tidy3d/components/eme/data/sim_data.py @@ -17,7 +17,7 @@ from tidy3d.exceptions import SetupError from tidy3d.log import log -from .dataset import EMESMatrixDataset +from .dataset import EMECoefficientDataset, EMESMatrixDataset from .monitor_data import EMEFieldData, EMEModeSolverData, EMEMonitorDataType @@ -39,6 +39,12 @@ class EMESimulationData(AbstractYeeGridSimulationData): None, title="S Matrix", description="Scattering matrix of the EME simulation." ) + coeffs: Optional[EMECoefficientDataset] = pd.Field( + None, + title="Coefficients", + description="Coefficients from the EME simulation. Useful for debugging and optimization.", + ) + port_modes_raw: Optional[EMEModeSolverData] = pd.Field( None, title="Port Modes", diff --git a/tidy3d/components/eme/simulation.py b/tidy3d/components/eme/simulation.py index 9f9c0841c7..cc15504f9a 100644 --- a/tidy3d/components/eme/simulation.py +++ b/tidy3d/components/eme/simulation.py @@ -233,6 +233,13 @@ class EMESimulation(AbstractYeeGridSimulation): "Required to find scattering matrix in basis besides the computational basis.", ) + store_coeffs: bool = pd.Field( + True, + title="Store Coefficients", + description="Whether to store the internal coefficients from the EME simulation. " + "The results are stored in 'EMESimulationData.coeffs'.", + ) + normalize: bool = pd.Field( True, title="Normalize Scattering Matrix", @@ -752,6 +759,12 @@ def _validate_sweep_spec(self) -> None: def _validate_monitor_setup(self) -> None: """Check monitor setup.""" for i, monitor in enumerate(self.monitors): + if isinstance(monitor, EMECoefficientMonitor): + log.warning( + "'EMECoefficientMonitor' is deprecated. " + "The full coefficient data is stored in " + "'EMESimulationData.coeffs'." + ) if isinstance(monitor, EMEMonitor): _ = self._monitor_eme_cell_indices(monitor=monitor) if (