diff --git a/CHANGELOG.md b/CHANGELOG.md index 9f84be6c38..d2974c928a 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -12,6 +12,12 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### Changed ### Fixed +- Fix to `outer_dot` when frequencies stored in the data were not in increasing order. Previously, the result would be provided with re-sorted frequencies, which would not match the order of the original data. +- Fixed bug where an extra spatial coordinate could appear in `complex_flux` and `ImpedanceCalculator` results. +- Fixed normal for `Box` shape gradient computation to always point outward from boundary which is needed for correct PEC handling. +- Fixed `Box` gradients within `GeometryGroup` where the group intersection boundaries were forwarded. +- Fixed `Box` gradients to use automatic permittivity detection for inside/outside permittivity. +- Fixed bug where `ModeSolver` and `Simulation` did not place PEC boundaries at the same location. The solver now truncates the computational grid to simulation bounds and zero-pads fields outside the domain. ## [2.10.0] - 2025-12-18 diff --git a/tests/test_plugins/test_mode_solver.py b/tests/test_plugins/test_mode_solver.py index cdc412513d..c5031a6950 100644 --- a/tests/test_plugins/test_mode_solver.py +++ b/tests/test_plugins/test_mode_solver.py @@ -15,6 +15,7 @@ from tidy3d.components.mode.derivatives import create_sfactor_b, create_sfactor_f from tidy3d.components.mode.solver import TOL_DEGENERATE_CANDIDATE, EigSolver, compute_modes from tidy3d.components.mode_spec import MODE_DATA_KEYS +from tidy3d.constants import fp_eps from tidy3d.exceptions import DataError, SetupError from tidy3d.plugins.mode import ModeSolver from tidy3d.plugins.mode.mode_solver import MODE_MONITOR_NAME @@ -1561,3 +1562,247 @@ def test_degenerate_mode_processing(): msg = f"Found {len(indices)} off-diagonal values > {threshold}:\n" msg += "\n".join(f" |S[{i},{j}]| = {np.abs(S[i, j]):.4e}" for i, j in indices) assert not np.any(problem_mask), msg + + +def test_mode_solver_pec_boundary_truncation(): + """Test that fields are correctly zero outside simulation bounds and satisfy PEC boundary conditions. + + This test creates a rectangular waveguide where the waveguide walls are modeled using the + PEC boundary conditions of the simulation domain. The mode solver grid may extend beyond + the simulation bounds due to _discretize_inds_monitor adding extra cells for interpolation. + This test verifies that: + 1. Fields outside the simulation boundaries are exactly 0 + 2. Electric fields tangential to the boundary are 0 at the boundary + 3. Magnetic fields normal to the boundary are 0 at the boundary + """ + # Create a simulation with PEC boundaries (default) + # The simulation size defines the waveguide cross-section + sim_size = (2.0, 0.0, 1.5) # Waveguide in y-direction, cross-section in x-z plane + freq0 = td.C_0 / 1.55 + + # Simple simulation with just air - the waveguide walls are the PEC boundaries + simulation = td.Simulation( + size=sim_size, + grid_spec=td.GridSpec.uniform(dl=0.05), + run_time=1e-14, + boundary_spec=td.BoundarySpec( + x=td.Boundary.pec(), + y=td.Boundary.periodic(), # Propagation direction + z=td.Boundary.pec(), + ), + sources=[ + td.PointDipole( + center=(0, 0, 0), + source_time=td.GaussianPulse(freq0=freq0, fwidth=freq0 / 10), + polarization="Ex", + ) + ], + ) + + # Mode plane perpendicular to propagation direction + plane = td.Box(center=(0, 0, 0), size=(sim_size[0], 0, sim_size[2])) + + mode_spec = td.ModeSpec( + num_modes=2, + precision="double", + ) + + ms = ModeSolver( + simulation=simulation, + plane=plane, + mode_spec=mode_spec, + freqs=[freq0], + direction="+", + colocate=False, + ) + + # Get the mode solver data + data = ms.data + + # Get simulation boundaries + sim_bounds = simulation.bounds + sim_x_min, sim_x_max = sim_bounds[0][0], sim_bounds[1][0] + sim_z_min, sim_z_max = sim_bounds[0][2], sim_bounds[1][2] + + # Test 1: Fields outside simulation boundaries should be exactly 0 (zero-padded) + for field_name in ("Ex", "Ey", "Ez", "Hx", "Hy", "Hz"): + field = data.field_components[field_name] + field_coords_x = field.coords["x"].values + field_coords_z = field.coords["z"].values + + # Check fields at x positions outside simulation + x_outside_min = field_coords_x[ + (field_coords_x < sim_x_min) + & ~np.isclose(field_coords_x, sim_x_min, rtol=fp_eps, atol=fp_eps) + ] + if len(x_outside_min) > 0: + field_outside = field.sel(x=x_outside_min) + assert np.all(field_outside.values == 0.0), ( + f"{field_name} should be exactly 0 outside x_min boundary, got {np.max(np.abs(field_outside.values))}" + ) + + x_outside_max = field_coords_x[ + (field_coords_x > sim_x_max) + & ~np.isclose(field_coords_x, sim_x_max, rtol=fp_eps, atol=fp_eps) + ] + if len(x_outside_max) > 0: + field_outside = field.sel(x=x_outside_max) + assert np.all(field_outside.values == 0.0), ( + f"{field_name} should be exactly 0 outside x_max boundary, got {np.max(np.abs(field_outside.values))}" + ) + + # Check fields at z positions outside simulation + z_outside_min = field_coords_z[ + (field_coords_z < sim_z_min) + & ~np.isclose(field_coords_z, sim_z_min, rtol=fp_eps, atol=fp_eps) + ] + if len(z_outside_min) > 0: + field_outside = field.sel(z=z_outside_min) + assert np.all(field_outside.values == 0.0), ( + f"{field_name} should be exactly 0 outside z_min boundary, got {np.max(np.abs(field_outside.values))}" + ) + + z_outside_max = field_coords_z[ + (field_coords_z > sim_z_max) + & ~np.isclose(field_coords_z, sim_z_max, rtol=fp_eps, atol=fp_eps) + ] + if len(z_outside_max) > 0: + field_outside = field.sel(z=z_outside_max) + assert np.all(field_outside.values == 0.0), ( + f"{field_name} should be exactly 0 outside z_max boundary, got {np.max(np.abs(field_outside.values))}" + ) + + # Test 2: Tangential E-fields at PEC boundaries should be exactly 0 + # At x boundaries: Ey and Ez are tangential + # At z boundaries: Ex and Ey are tangential + + # Get field coordinates exactly at boundaries + for field_name in ("Ey", "Ez"): + field = data.field_components[field_name] + field_coords_x = field.coords["x"].values + # Find coordinates exactly at boundaries + x_at_min = field_coords_x[np.isclose(field_coords_x, sim_x_min, rtol=fp_eps, atol=fp_eps)] + x_at_max = field_coords_x[np.isclose(field_coords_x, sim_x_max, rtol=fp_eps, atol=fp_eps)] + + if len(x_at_min) > 0: + field_at_boundary = field.sel(x=x_at_min[0]) + assert np.all(field_at_boundary.values == 0.0), ( + f"{field_name} (tangential) should be exactly 0 at x_min PEC boundary, got {np.max(np.abs(field_at_boundary.values))}" + ) + + if len(x_at_max) > 0: + field_at_boundary = field.sel(x=x_at_max[0]) + assert np.all(field_at_boundary.values == 0.0), ( + f"{field_name} (tangential) should be exactly 0 at x_max PEC boundary, got {np.max(np.abs(field_at_boundary.values))}" + ) + + for field_name in ("Ex", "Ey"): + field = data.field_components[field_name] + field_coords_z = field.coords["z"].values + # Find coordinates exactly at boundaries + z_at_min = field_coords_z[np.isclose(field_coords_z, sim_z_min, rtol=fp_eps, atol=fp_eps)] + z_at_max = field_coords_z[np.isclose(field_coords_z, sim_z_max, rtol=fp_eps, atol=fp_eps)] + + if len(z_at_min) > 0: + field_at_boundary = field.sel(z=z_at_min[0]) + assert np.all(field_at_boundary.values == 0.0), ( + f"{field_name} (tangential) should be exactly 0 at z_min PEC boundary, got {np.max(np.abs(field_at_boundary.values))}" + ) + + if len(z_at_max) > 0: + field_at_boundary = field.sel(z=z_at_max[0]) + assert np.all(field_at_boundary.values == 0.0), ( + f"{field_name} (tangential) should be exactly 0 at z_max PEC boundary, got {np.max(np.abs(field_at_boundary.values))}" + ) + + # Test 3: Normal H-fields at PEC boundaries should be exactly 0 + # At x boundaries: Hx is normal + # At z boundaries: Hz is normal + field = data.field_components["Hx"] + field_coords_x = field.coords["x"].values + x_at_min = field_coords_x[np.isclose(field_coords_x, sim_x_min, rtol=fp_eps, atol=fp_eps)] + x_at_max = field_coords_x[np.isclose(field_coords_x, sim_x_max, rtol=fp_eps, atol=fp_eps)] + + if len(x_at_min) > 0: + field_at_boundary = field.sel(x=x_at_min[0]) + assert np.all(field_at_boundary.values == 0.0), ( + f"Hx (normal) should be exactly 0 at x_min PEC boundary, got {np.max(np.abs(field_at_boundary.values))}" + ) + + if len(x_at_max) > 0: + field_at_boundary = field.sel(x=x_at_max[0]) + assert np.all(field_at_boundary.values == 0.0), ( + f"Hx (normal) should be exactly 0 at x_max PEC boundary, got {np.max(np.abs(field_at_boundary.values))}" + ) + + field = data.field_components["Hz"] + field_coords_z = field.coords["z"].values + z_at_min = field_coords_z[np.isclose(field_coords_z, sim_z_min, rtol=fp_eps, atol=fp_eps)] + z_at_max = field_coords_z[np.isclose(field_coords_z, sim_z_max, rtol=fp_eps, atol=fp_eps)] + + if len(z_at_min) > 0: + field_at_boundary = field.sel(z=z_at_min[0]) + assert np.all(field_at_boundary.values == 0.0), ( + f"Hz (normal) should be exactly 0 at z_min PEC boundary, got {np.max(np.abs(field_at_boundary.values))}" + ) + + if len(z_at_max) > 0: + field_at_boundary = field.sel(z=z_at_max[0]) + assert np.all(field_at_boundary.values == 0.0), ( + f"Hz (normal) should be exactly 0 at z_max PEC boundary, got {np.max(np.abs(field_at_boundary.values))}" + ) + + +@pytest.mark.parametrize( + "boundary_type,boundary_axis,warning_str,plane_size,expected_warnings", + [ + ("pmc", "x", "PMC", (6, 0, 6), 2), + ("periodic", "z", "periodic", (6, 0, 6), 2), + ("bloch", "x", "Bloch", (6, 0, 6), 2), + ("pmc", "x", None, (2, 0, 2), 0), # No warning when plane doesn't intersect boundary + ], + ids=["pmc", "periodic", "bloch", "no_intersection"], +) +def test_mode_solver_boundary_warning( + boundary_type, boundary_axis, warning_str, plane_size, expected_warnings +): + """Test that warnings are emitted when mode solver plane intersects unsupported boundaries.""" + # Set up boundary spec based on test parameters + if boundary_type == "pmc": + test_boundary = td.Boundary.pmc() + elif boundary_type == "periodic": + test_boundary = td.Boundary.periodic() + elif boundary_type == "bloch": + test_boundary = td.Boundary.bloch(bloch_vec=0.1) + + # Apply the test boundary to the specified axis, PEC to others + boundary_spec = td.BoundarySpec( + x=test_boundary if boundary_axis == "x" else td.Boundary.pec(), + y=td.Boundary.pec(), + z=test_boundary if boundary_axis == "z" else td.Boundary.pec(), + ) + + simulation = td.Simulation( + size=(4, 3, 3), + grid_spec=td.GridSpec.auto(wavelength=1.0), + structures=[WAVEGUIDE], + run_time=1e-12, + boundary_spec=boundary_spec, + ) + + plane = td.Box(center=(0, 0, 0), size=plane_size) + + mode_spec = td.ModeSpec(num_modes=1) + ms = ModeSolver( + simulation=simulation, + plane=plane, + mode_spec=mode_spec, + freqs=[td.C_0 / 1.0], + ) + + # Access the cached property to trigger the warning + with AssertLogLevel( + "WARNING" if expected_warnings > 0 else None, contains_str=warning_str + ) as ctx: + _ = ms._sim_boundary_positions + assert ctx.num_records == expected_warnings diff --git a/tidy3d/components/mode/mode_solver.py b/tidy3d/components/mode/mode_solver.py index 5acc5f4e64..12a571b82f 100644 --- a/tidy3d/components/mode/mode_solver.py +++ b/tidy3d/components/mode/mode_solver.py @@ -17,7 +17,17 @@ cached_property, skip_if_fields_missing, ) -from tidy3d.components.boundary import PML, Absorber, Boundary, BoundarySpec, PECBoundary, StablePML +from tidy3d.components.boundary import ( + PML, + Absorber, + BlochBoundary, + Boundary, + BoundarySpec, + PECBoundary, + Periodic, + PMCBoundary, + StablePML, +) from tidy3d.components.data.data_array import ( FreqModeDataArray, ModeIndexDataArray, @@ -65,6 +75,7 @@ Ax, Axis, Axis2D, + Bound2D, Direction, EMField, EpsSpecType, @@ -1171,6 +1182,88 @@ def _reduced_simulation_copy_with_fallback(self) -> ModeSolver: ) return self + @cached_property + def _sim_boundary_positions(self) -> Bound2D: + """Get simulation boundary positions for the mode plane's tangential axes. + + Returns the simulation boundary positions and logs warnings if the mode solver + grid extends beyond boundaries with unsupported conditions (PMC, Periodic, Bloch). + PEC boundaries are fully supported. PML/Absorber/ABC boundaries are silently + ignored since including them would introduce too many warnings in existing simulations. + + .. TODO Consolidate all boundary conditions in the underlying simulation and the mode solver. + + Returns + ------- + Bound2D + ((min_0, min_1), (max_0, max_1)) positions along tangential axes. + """ + trans_axes = [0, 1, 2] + trans_axes.remove(self.normal_axis) + axis_names = ["x", "y", "z"] + + sim_grid_list = self.simulation.grid.boundaries.to_list + solver_grid_list = self._solver_grid.boundaries.to_list + bspec = self.simulation.boundary_spec + + pos_min = [None, None] + pos_max = [None, None] + + for i, axis in enumerate(trans_axes): + axis_name = axis_names[axis] + boundary = bspec[axis_name] + + sim_min = sim_grid_list[axis][0] + sim_max = sim_grid_list[axis][-1] + solver_min = solver_grid_list[axis][0] + solver_max = solver_grid_list[axis][-1] + pos_min[i] = sim_min + pos_max[i] = sim_max + # Check if mode solver plane intersects simulation boundaries + # Min side: intersects if solver grid extends beyond (below) sim boundary + intersects_min = solver_min < sim_min and not np.isclose( + solver_min, sim_min, rtol=fp_eps, atol=fp_eps + ) + # Max side: intersects if solver grid extends beyond (above) sim boundary + intersects_max = solver_max > sim_max and not np.isclose( + solver_max, sim_max, rtol=fp_eps, atol=fp_eps + ) + + # Check minus side + bc_minus = boundary.minus + if isinstance(bc_minus, PMCBoundary): + if intersects_min: + log.warning( + f"Mode solver plane intersects simulation '{axis_name}-' boundary with " + f"unsupported PMC boundary condition. The mode solver will apply PEC instead." + ) + elif isinstance(bc_minus, (Periodic, BlochBoundary)): + if intersects_min: + log.warning( + f"Mode solver plane intersects simulation '{axis_name}-' boundary with " + f"unsupported periodic/Bloch boundary condition. " + f"Fields will not wrap around periodically." + ) + # ABC boundaries: no truncation, no warning (fields can extend) + + # Check plus side + bc_plus = boundary.plus + if isinstance(bc_plus, PMCBoundary): + if intersects_max: + log.warning( + f"Mode solver plane intersects simulation '{axis_name}+' boundary with " + f"unsupported PMC boundary condition. The mode solver will apply PEC instead." + ) + elif isinstance(bc_plus, (Periodic, BlochBoundary)): + if intersects_max: + log.warning( + f"Mode solver plane intersects simulation '{axis_name}+' boundary with " + f"unsupported periodic/Bloch boundary condition. " + f"Fields will not wrap around periodically." + ) + # ABC boundaries: no truncation, no warning (fields can extend) + return (tuple(pos_min), tuple(pos_max)) + def _data_on_yee_grid(self) -> ModeSolverData: """Solve for all modes, and construct data with fields on the Yee grid.""" solver = self._reduced_simulation_copy_with_fallback @@ -1187,7 +1280,8 @@ def _data_on_yee_grid(self) -> ModeSolverData: # Compute and store the modes at all frequencies n_complex, fields, eps_spec = solver._solve_all_freqs( - coords=_solver_coords, symmetry=solver.solver_symmetry + coords=_solver_coords, + symmetry=solver.solver_symmetry, ) # start a dictionary storing the data arrays for the ModeSolverData @@ -1261,7 +1355,9 @@ def _data_on_yee_grid_relative(self, basis: ModeSolverData) -> ModeSolverData: # Compute and store the modes at all frequencies n_complex, fields, eps_spec = self._solve_all_freqs_relative( - coords=_solver_coords, symmetry=self.solver_symmetry, basis_fields=basis_fields + coords=_solver_coords, + symmetry=self.solver_symmetry, + basis_fields=basis_fields, ) # start a dictionary storing the data arrays for the ModeSolverData @@ -1571,14 +1667,19 @@ def _solve_all_freqs( """Call the mode solver at all requested frequencies.""" if tidy3d_extras["use_local_subpixel"]: subpixel_ms = tidy3d_extras["mod"].SubpixelModeSolver.from_mode_solver(self) - return subpixel_ms._solve_all_freqs(coords=coords, symmetry=symmetry) + return subpixel_ms._solve_all_freqs( + coords=coords, + symmetry=symmetry, + ) fields = [] n_complex = [] eps_spec = [] for freq in self.freqs: n_freq, fields_freq, eps_spec_freq = self._solve_single_freq( - freq=freq, coords=coords, symmetry=symmetry + freq=freq, + coords=coords, + symmetry=symmetry, ) fields.append(fields_freq) n_complex.append(n_freq) @@ -1596,7 +1697,9 @@ def _solve_all_freqs_relative( if tidy3d_extras["use_local_subpixel"]: subpixel_ms = tidy3d_extras["mod"].SubpixelModeSolver.from_mode_solver(self) return subpixel_ms._solve_all_freqs_relative( - coords=coords, symmetry=symmetry, basis_fields=basis_fields + coords=coords, + symmetry=symmetry, + basis_fields=basis_fields, ) fields = [] @@ -1604,7 +1707,10 @@ def _solve_all_freqs_relative( eps_spec = [] for freq, basis_fields_freq in zip(self.freqs, basis_fields): n_freq, fields_freq, eps_spec_freq = self._solve_single_freq_relative( - freq=freq, coords=coords, symmetry=symmetry, basis_fields=basis_fields_freq + freq=freq, + coords=coords, + symmetry=symmetry, + basis_fields=basis_fields_freq, ) fields.append(fields_freq) n_complex.append(n_freq) @@ -1656,6 +1762,7 @@ def _solve_single_freq( direction=self.direction, precision=self._precision, plane_center=self.plane_center_tangential(self.plane), + sim_pec_bound=self._sim_boundary_positions, ) fields = self._postprocess_solver_fields( @@ -1718,6 +1825,7 @@ def _solve_single_freq_relative( solver_basis_fields=solver_basis_fields, precision=self._precision, plane_center=self.plane_center_tangential(self.plane), + sim_pec_bound=self._sim_boundary_positions, ) fields = self._postprocess_solver_fields( diff --git a/tidy3d/components/mode/solver.py b/tidy3d/components/mode/solver.py index e533c3d084..ab10df1574 100644 --- a/tidy3d/components/mode/solver.py +++ b/tidy3d/components/mode/solver.py @@ -7,7 +7,7 @@ import numpy as np from tidy3d.components.base import Tidy3dBaseModel -from tidy3d.components.types import EpsSpecType, ModeSolverType, Numpy +from tidy3d.components.types import Bound2D, EpsSpecType, ModeSolverType, Numpy from tidy3d.constants import C_0, ETA_0, fp_eps, pec_val from .derivatives import create_d_matrices as d_mats @@ -58,6 +58,7 @@ def compute_modes( direction="+", solver_basis_fields=None, plane_center: Optional[tuple[float, float]] = None, + sim_pec_bound: Optional[Bound2D] = None, ) -> tuple[Numpy, Numpy, EpsSpecType]: """ Solve for the modes of a waveguide cross-section. @@ -97,6 +98,9 @@ def compute_modes( The center of the mode plane along the tangential axes of the global simulation. Used in case of bend modes to offset the coordinates correctly w.r.t. the bend radius, which is assumed to refer to the distance from the bend center to the mode plane center. + sim_pec_bound : Optional[Bound2D] + Simulation PEC boundary positions ((x_min, y_min), (x_max, y_max)). If the solver grid extends + beyond a PEC position, it will be truncated and fields outside will be zero-padded. Returns ------- @@ -115,6 +119,67 @@ def compute_modes( k0 = omega / C_0 enable_incidence_matrices = False # Experimental feature, always off for now + # Determine truncation indices based on simulation PEC boundary positions + # Store original shapes for later zero-padding + original_Nx = len(coords[0]) - 1 + original_Ny = len(coords[1]) - 1 + + # Find indices where coords fall within simulation PEC boundaries + trim_min = [0, 0] + trim_max = [len(coords[0]), len(coords[1])] + + if sim_pec_bound is not None: + # Find first coord index >= min_bound (for boundary coords) + min_bound = sim_pec_bound[0] + for i in range(2): + idx = np.searchsorted(coords[i], min_bound[i], side="left") + # Only trim if min_bound is actually inside the coord range + if idx > 0 and not np.isclose(coords[i][idx - 1], min_bound[i]): + trim_min[i] = idx + + # Find last coord index <= max_bound (for boundary coords) + max_bound = sim_pec_bound[1] + for i in range(2): + idx = np.searchsorted(coords[i], max_bound[i], side="right") + # Only trim if max_bound is actually inside the coord range + if idx < len(coords[i]) and not np.isclose(coords[i][idx], max_bound[i]): + trim_max[i] = idx + + # Check if truncation is needed + needs_truncation = trim_min != [0, 0] or trim_max != [len(coords[0]), len(coords[1])] + + if needs_truncation: + # Truncate coords + coords = [ + coords[0][trim_min[0] : trim_max[0]], + coords[1][trim_min[1] : trim_max[1]], + ] + + # Calculate cell indices for eps truncation (cells are between boundaries) + cell_min = [trim_min[0], trim_min[1]] + cell_max = [trim_max[0] - 1, trim_max[1] - 1] + + # Truncate eps_cross + eps_cross = cls._truncate_medium_data(eps_cross, cell_min, cell_max) + + # Truncate mu_cross if provided + if mu_cross is not None: + mu_cross = cls._truncate_medium_data(mu_cross, cell_min, cell_max) + + # Truncate split_curl_scaling if provided + if split_curl_scaling is not None: + split_curl_scaling = tuple( + s[cell_min[0] : cell_max[0], cell_min[1] : cell_max[1]] + for s in split_curl_scaling + ) + + # Truncate solver_basis_fields if provided + # Shape is (6, Nx, Ny, 1, num_modes) for E and H fields + if solver_basis_fields is not None: + solver_basis_fields = solver_basis_fields[ + :, cell_min[0] : cell_max[0], cell_min[1] : cell_max[1], :, : + ] + eps_formated = cls.format_medium_data(eps_cross) eps_xx, eps_xy, eps_xz, eps_yx, eps_yy, eps_yz, eps_zx, eps_zy, eps_zz = eps_formated @@ -292,6 +357,19 @@ def compute_modes( fields = np.stack((E, H), axis=0) + # Zero-pad fields back to original size if truncation was applied + if needs_truncation: + # fields shape: (2, 3, Nx, Ny, 1, num_modes) for E and H stacked + pad_width = ( + (0, 0), # E/H axis + (0, 0), # field component axis (Ex, Ey, Ez) + (trim_min[0], original_Nx - trim_max[0] + 1), # x axis padding + (trim_min[1], original_Ny - trim_max[1] + 1), # y axis padding + (0, 0), # normal axis (always 1) + (0, 0), # modes axis + ) + fields = np.pad(fields, pad_width, mode="constant", constant_values=0.0) + neff = neff * np.linalg.norm(kp_to_k) keff = keff * np.linalg.norm(kp_to_k) @@ -1066,6 +1144,35 @@ def eigs_to_effective_index(cls, eig_list: Numpy, mode_solver_type: ModeSolverTy raise RuntimeError(f"Unidentified 'mode_solver_type={mode_solver_type}'.") + @staticmethod + def _truncate_medium_data(mat_data, cell_min, cell_max): + """Truncate medium data (eps or mu) to the specified cell range. + + Parameters + ---------- + mat_data : array_like or tuple of array_like + Either a single 2D array or nine 2D arrays (for tensorial data). + cell_min : list[int] + [x_min, y_min] cell indices for truncation start. + cell_max : list[int] + [x_max, y_max] cell indices for truncation end. + + Returns + ------- + Truncated medium data in the same format as input. + """ + x_slice = slice(cell_min[0], cell_max[0]) + y_slice = slice(cell_min[1], cell_max[1]) + + if isinstance(mat_data, Numpy): + # Tensorial format: shape (9, Nx, Ny) + return mat_data[:, x_slice, y_slice] + if len(mat_data) == 9: + # Nine separate 2D arrays + return tuple(arr[x_slice, y_slice] for arr in mat_data) + + raise ValueError("Wrong input to mode solver permittivity/permeability truncation!") + @staticmethod def format_medium_data(mat_data): """ diff --git a/tidy3d/components/types/__init__.py b/tidy3d/components/types/__init__.py index 0d7c70e94c..027516c71c 100644 --- a/tidy3d/components/types/__init__.py +++ b/tidy3d/components/types/__init__.py @@ -19,6 +19,7 @@ Axis, Axis2D, Bound, + Bound2D, BoxSurface, ClipOperationType, ColormapType, @@ -88,6 +89,7 @@ "Axis", "Axis2D", "Bound", + "Bound2D", "BoxSurface", "ClipOperationType", "ColormapType", diff --git a/tidy3d/components/types/base.py b/tidy3d/components/types/base.py index b49e5d7086..dd9b1105bd 100644 --- a/tidy3d/components/types/base.py +++ b/tidy3d/components/types/base.py @@ -192,6 +192,7 @@ def __modify_schema__(cls, field_schema) -> None: CoordinateOptional = tuple[Optional[float], Optional[float], Optional[float]] Coordinate2D = tuple[float, float] Bound = tuple[Coordinate, Coordinate] +Bound2D = tuple[Coordinate2D, Coordinate2D] GridSize = Union[pydantic.PositiveFloat, tuple[pydantic.PositiveFloat, ...]] Axis = Literal[0, 1, 2] Axis2D = Literal[0, 1]