From 4384ea858454119e362a2e2b0d897e7f63e028b1 Mon Sep 17 00:00:00 2001 From: dmarek Date: Fri, 14 Nov 2025 11:54:28 -0500 Subject: [PATCH] fix(mode): fixed bug where mode solver was not respecting location of PEC boundary conditions from the Simulation --- CHANGELOG.md | 1 + tests/test_plugins/test_mode_solver.py | 323 +++++++++++++++++++++++++ tidy3d/components/mode/mode_solver.py | 140 ++++++++++- tidy3d/components/mode/solver.py | 110 +++++++++ 4 files changed, 567 insertions(+), 7 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index bf12238276..fca9af74dd 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -108,6 +108,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 - Removed mode solver warnings about evaluating permittivity of a `Medium2D`. - Maximum number of grid points in an EME simulation is now based solely on transverse grid points. Maximum number of EME cells is unchanged. - Fixed handling of polygons with holes (interiors) in `subdivide()` function. The function now properly converts polygons with interiors into `PolySlab` geometries using subtraction operations with `GeometryGroup`. +- Fixed mode solver grid extending beyond simulation domain boundaries, causing incorrect PEC boundary application. The solver now properly truncates the computational grid to simulation bounds and zero-pads fields outside the domain. ### Removed - Removed deprecated `use_complex_fields` parameter from `TwoPhotonAbsorption` and `KerrNonlinearity`. Parameters `beta` and `n2` are now real-valued only, as is `n0` if specified. diff --git a/tests/test_plugins/test_mode_solver.py b/tests/test_plugins/test_mode_solver.py index 97155fdd56..9cf945dcef 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 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 @@ -1503,3 +1504,325 @@ def test_sort_spec_track_freq(): assert np.allclose(modes_lowest.Ex.abs, modes_lowest_retracked.Ex.abs) assert np.all(modes_lowest.n_eff == modes_lowest_retracked.n_eff) assert np.all(modes_lowest.n_group == modes_lowest_retracked.n_group) + + +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))}" + ) + + +def test_mode_solver_boundary_warning_pmc(): + """Test that warnings are emitted when mode solver plane intersects unsupported PMC boundaries.""" + # Create a simulation with PMC boundaries on x + simulation = td.Simulation( + size=(4, 3, 3), + grid_spec=td.GridSpec.auto(wavelength=1.0), + structures=[WAVEGUIDE], + run_time=1e-12, + boundary_spec=td.BoundarySpec( + x=td.Boundary.pmc(), + y=td.Boundary.pec(), + z=td.Boundary.pec(), + ), + ) + + # Mode plane in xz plane that extends to simulation x boundaries + # The plane is larger than simulation, so it will intersect the PMC boundaries + plane = td.Box(center=(0, 0, 0), size=(6, 0, 6)) + + 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", contains_str="PMC") as ctx: + _ = ms._sim_boundary_positions + # Should warn for both x- and x+ sides + assert ctx.num_records == 2 + + +def test_mode_solver_boundary_warning_periodic(): + """Test that warnings are emitted when mode solver plane intersects unsupported periodic boundaries.""" + # Create a simulation with periodic boundaries on z + simulation = td.Simulation( + size=(4, 3, 3), + grid_spec=td.GridSpec.auto(wavelength=1.0), + structures=[WAVEGUIDE], + run_time=1e-12, + boundary_spec=td.BoundarySpec( + x=td.Boundary.pec(), + y=td.Boundary.pec(), + z=td.Boundary.periodic(), + ), + ) + + # Mode plane in xz plane that extends to simulation z boundaries + plane = td.Box(center=(0, 0, 0), size=(6, 0, 6)) + + 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", contains_str="periodic") as ctx: + _ = ms._sim_boundary_positions + # Should warn for both z- and z+ sides + assert ctx.num_records == 2 + + +def test_mode_solver_boundary_warning_bloch(): + """Test that warnings are emitted when mode solver plane intersects unsupported Bloch boundaries.""" + # Create a simulation with Bloch boundaries on x + simulation = td.Simulation( + size=(4, 3, 3), + grid_spec=td.GridSpec.auto(wavelength=1.0), + structures=[WAVEGUIDE], + run_time=1e-12, + boundary_spec=td.BoundarySpec( + x=td.Boundary.bloch(bloch_vec=0.1), + y=td.Boundary.pec(), + z=td.Boundary.pec(), + ), + ) + + # Mode plane in xz plane that extends to simulation x boundaries + plane = td.Box(center=(0, 0, 0), size=(6, 0, 6)) + + 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", contains_str="Bloch") as ctx: + _ = ms._sim_boundary_positions + # Should warn for both x- and x+ sides + assert ctx.num_records == 2 + + +def test_mode_solver_boundary_no_warning_when_not_intersecting(): + """Test that no warnings are emitted when mode solver plane doesn't intersect unsupported boundaries.""" + # Create a simulation with PMC boundaries on x + simulation = td.Simulation( + size=(4, 3, 3), + grid_spec=td.GridSpec.auto(wavelength=1.0), + structures=[WAVEGUIDE], + run_time=1e-12, + boundary_spec=td.BoundarySpec( + x=td.Boundary.pmc(), + y=td.Boundary.pec(), + z=td.Boundary.pec(), + ), + ) + + # Mode plane in xz plane that is smaller than simulation x bounds + # The plane doesn't extend to the PMC boundaries + plane = td.Box(center=(0, 0, 0), size=(2, 0, 2)) + + 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 - should NOT trigger any warning + with AssertLogLevel(None) as ctx: + _ = ms._sim_boundary_positions + assert ctx.num_records == 0 diff --git a/tidy3d/components/mode/mode_solver.py b/tidy3d/components/mode/mode_solver.py index 5acc5f4e64..b004065b67 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, @@ -1171,6 +1181,99 @@ def _reduced_simulation_copy_with_fallback(self) -> ModeSolver: ) return self + @cached_property + def _sim_boundary_positions( + self, + ) -> tuple[tuple[float, float], tuple[float, float]]: + """Get simulation boundary positions for the mode plane's tangential axes. + + This method checks the simulation's boundary conditions and returns positions + only for sides that have supported boundaries. For unsupported boundary types, + a warning is logged only if the mode solver plane intersects that boundary. + + The mode solver supports: + - PEC boundaries: Fields are truncated and zero-padded at these boundaries + - Symmetry: Handled via solver_symmetry parameter + - PML/Absorber: Converted to PEC in reduced_simulation_copy (already handled) + + Unsupported boundary types that will trigger a warning (if intersected): + - PMC: Mode solver doesn't implement PMC boundary conditions (treated as PEC) + - Periodic/Bloch: Mode solver doesn't implement periodic field wrapping + - ABC: Mode solver doesn't implement absorbing boundary conditions + + Returns + ------- + tuple[tuple[float, float], tuple[float, float]] + (pos_min, pos_max) where: + - pos_min = (x_min, y_min) positions for supported boundaries at min edges (None otherwise) + - pos_max = (x_max, y_max) positions for supported boundaries at max edges (None otherwise) + """ + 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 +1290,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 +1365,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 +1677,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 +1707,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 +1717,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) @@ -1647,6 +1763,9 @@ def _solve_single_freq( if not LOCAL_SOLVER_IMPORTED: raise ImportError(IMPORT_ERROR_MSG) + # Get simulation boundary positions for potential grid truncation + sim_pec_pos_min, sim_pec_pos_max = self._sim_boundary_positions + solver_fields, n_complex, eps_spec = compute_modes( eps_cross=self._solver_eps(freq), coords=coords, @@ -1656,6 +1775,8 @@ def _solve_single_freq( direction=self.direction, precision=self._precision, plane_center=self.plane_center_tangential(self.plane), + sim_pec_pos_min=sim_pec_pos_min, + sim_pec_pos_max=sim_pec_pos_max, ) fields = self._postprocess_solver_fields( @@ -1708,6 +1829,9 @@ def _solve_single_freq_relative( fields=basis_fields, normal_axis=self.normal_axis, plane=self.plane ) + # Get simulation boundary positions for potential grid truncation + sim_pec_pos_min, sim_pec_pos_max = self._sim_boundary_positions + solver_fields, n_complex, eps_spec = compute_modes( eps_cross=self._solver_eps(freq), coords=coords, @@ -1718,6 +1842,8 @@ def _solve_single_freq_relative( solver_basis_fields=solver_basis_fields, precision=self._precision, plane_center=self.plane_center_tangential(self.plane), + sim_pec_pos_min=sim_pec_pos_min, + sim_pec_pos_max=sim_pec_pos_max, ) fields = self._postprocess_solver_fields( diff --git a/tidy3d/components/mode/solver.py b/tidy3d/components/mode/solver.py index bb107d9d16..219095d10d 100644 --- a/tidy3d/components/mode/solver.py +++ b/tidy3d/components/mode/solver.py @@ -54,6 +54,8 @@ def compute_modes( direction="+", solver_basis_fields=None, plane_center: Optional[tuple[float, float]] = None, + sim_pec_pos_min: Optional[tuple[float, float]] = None, + sim_pec_pos_max: Optional[tuple[float, float]] = None, ) -> tuple[Numpy, Numpy, EpsSpecType]: """ Solve for the modes of a waveguide cross-section. @@ -93,6 +95,12 @@ 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_pos_min : Optional[tuple[float, float]] + Simulation PEC boundary positions (x_min, y_min). If the solver grid extends beyond a PEC + position, it will be truncated and fields outside will be zero-padded. + sim_pec_pos_max : Optional[tuple[float, float]] + Simulation PEC boundary positions (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 ------- @@ -111,6 +119,66 @@ 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_pos_min is not None: + # Find first coord index >= sim_pec_pos_min (for boundary coords) + for i in range(2): + idx = np.searchsorted(coords[i], sim_pec_pos_min[i]) + # Only trim if sim_pec_pos is actually inside the coord range + if idx > 0 and not np.isclose(coords[i][idx - 1], sim_pec_pos_min[i]): + trim_min[i] = idx + + if sim_pec_pos_max is not None: + # Find last coord index <= sim_pec_pos_max (for boundary coords) + for i in range(2): + idx = np.searchsorted(coords[i], sim_pec_pos_max[i], side="right") + # Only trim if sim_pec_pos is actually inside the coord range + if idx < len(coords[i]) and not np.isclose(coords[i][idx], sim_pec_pos_max[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 @@ -287,6 +355,19 @@ def compute_modes( H = H.reshape((3, Nx, Ny, 1, num_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) @@ -1029,6 +1110,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): """