diff --git a/CHANGELOG.md b/CHANGELOG.md index 3b82677c60..5df9fcc582 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -45,6 +45,12 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 - Changed the default value of `fill_value` in `UnstructuredGridDataset.interp()` from `0` to `"extrapolate"`. This means points outside the mesh will now use nearest-neighbor extrapolation instead of being filled with zeros. ### 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. - Fixed `CustomMedium` gradient calculation when field coordinates exactly align with boundaries. - Fixed adjoint simulation `grid_spec` to align exactly with forward simulation for correct `FieldData` adjoint source power. diff --git a/tests/test_components/test_microwave.py b/tests/test_components/test_microwave.py index f43aab6d73..b7ca47d8c3 100644 --- a/tests/test_components/test_microwave.py +++ b/tests/test_components/test_microwave.py @@ -918,26 +918,61 @@ def test_mode_plane_analyzer_mode_bounds(mode_size, symmetry): colocate=True, freqs=[freq0], ) - mode_solver_boundaries = mms._solver_grid.boundaries.to_list mode_plane_analyzer = ModePlaneAnalyzer( center=mode_center, size=mode_size, field_data_colocated=False, ) - mode_plane_limits = mode_plane_analyzer._get_mode_limits(sim.grid, sim.symmetry) - - for dim in (0, 1): - solver_dim_boundaries = mode_solver_boundaries[dim] - # TODO: Need the second check because the mode solver erroneously adds - # an extra grid cell even when touching the simulation boundary - assert ( - solver_dim_boundaries[0] == mode_plane_limits[0][dim] - or mode_plane_limits[0][dim] == sim.bounds[0][dim] - ) - assert ( - solver_dim_boundaries[-1] == mode_plane_limits[1][dim] - or mode_plane_limits[1][dim] == sim.bounds[1][dim] - ) + # Verify that ModePlaneAnalyzer PEC bounds match ModeSolver PEC bounds + mode_symmetry_3d = mode_plane_analyzer._get_mode_symmetry(sim.bounding_box, sim.symmetry) + analyzer_pec = mode_plane_analyzer._get_pec_boundary_positions(sim.grid, mode_symmetry_3d) + solver_pec = mms._pec_boundary_positions + _, tangential_axes = td.Box.pop_axis([0, 1, 2], mms.normal_axis) + for i, ax in enumerate(tangential_axes): + # Max side should always match + assert np.isclose(analyzer_pec[1][ax], solver_pec[1][i]) + # Min side matches only when no symmetry (with symmetry, analyzer + # moves min to center for half-domain conductor filtering) + if mode_symmetry_3d[ax] == 0: + assert np.isclose(analyzer_pec[0][ax], solver_pec[0][i]) + + +def test_mode_plane_analyzer_ground_plane_filtered(): + """Test that a ground plane extending beyond the mode plane is correctly filtered + as shorted to PEC boundary. Regression test for SCRF-2777.""" + # Microstrip: narrow signal trace + wide ground plane, with PEC boundaries + signal_trace = td.Structure( + geometry=td.Box(center=(0, 0, 1.1 * mm), size=(0.4 * mm, td.inf, 35)), + medium=td.PEC, + ) + ground_plane = td.Structure( + geometry=td.Box(center=(0, 0, 0), size=(1.5 * mm, td.inf, 35)), + medium=td.PEC, + ) + substrate = td.Structure( + geometry=td.Box(center=(0, 0, 0.55 * mm), size=(2 * mm, td.inf, 1.1 * mm)), + medium=td.Medium(permittivity=4.4), + ) + sim = td.Simulation( + center=(0, 0, 1 * mm), + size=(2 * mm, 2 * mm, 3 * mm), + grid_spec=td.GridSpec.uniform(dl=0.05 * mm), + structures=[substrate, signal_trace, ground_plane], + run_time=1e-12, + boundary_spec=td.BoundarySpec( + x=td.Boundary.pec(), y=td.Boundary.pec(), z=td.Boundary.pec() + ), + ) + # Mode plane is narrower than the ground plane in x + mode_plane = td.Box(center=(0, 0, 1 * mm), size=(1.2 * mm, 0, 2.5 * mm)) + analyzer = ModePlaneAnalyzer( + center=mode_plane.center, size=mode_plane.size, field_data_colocated=False + ) + bounding_boxes, conductors = analyzer.get_conductor_bounding_boxes( + sim.structures, sim.grid, sim.symmetry, sim.bounding_box + ) + # Only the signal trace should remain; ground plane is shorted to PEC + assert len(conductors) == 1 def test_impedance_spec_validation(): diff --git a/tests/test_plugins/test_mode_solver.py b/tests/test_plugins/test_mode_solver.py index b30fb2b37c..2f902d4452 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, ValidationError from tidy3d.plugins.mode import ModeSolver from tidy3d.plugins.mode.mode_solver import MODE_MONITOR_NAME @@ -1823,3 +1824,402 @@ def test_mode_data_fill_fraction_box_requires_intersection(): data, _ = make_fill_fraction_mode_data() with pytest.raises(ValidationError): data.fill_fraction(td.Box(center=(0.0, 2.0, 0.0), size=(1.0, 1.0, 1.0))) + + +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._pec_boundary_positions + assert ctx.num_records == expected_warnings + + +def test_mode_solver_pec_boundary_interior_plane(): + """Test that PEC snaps to the mode plane edges for an interior plane (smaller than simulation).""" + freq0 = td.C_0 / 1.55 + dl = 0.05 + + # Large simulation with PEC boundaries + simulation = td.Simulation( + size=(4.0, 0.0, 4.0), + grid_spec=td.GridSpec.uniform(dl=dl), + run_time=1e-14, + boundary_spec=td.BoundarySpec( + x=td.Boundary.pec(), + y=td.Boundary.periodic(), + z=td.Boundary.pec(), + ), + sources=[ + td.PointDipole( + center=(0, 0, 0), + source_time=td.GaussianPulse(freq0=freq0, fwidth=freq0 / 10), + polarization="Ex", + ) + ], + ) + + # Interior plane much smaller than the simulation + plane = td.Box(center=(0, 0, 0), size=(2.0, 0, 1.5)) + + ms = ModeSolver( + simulation=simulation, + plane=plane, + mode_spec=td.ModeSpec(num_modes=1, precision="double"), + freqs=[freq0], + direction="+", + colocate=False, + ) + + pec_bounds = ms._pec_boundary_positions + sim_bounds = simulation.bounds + + # PEC should snap to plane edges, NOT simulation edges + # x-axis: plane spans [-1.0, 1.0], sim spans [-2.0, 2.0] + assert pec_bounds[0][0] <= plane.bounds[0][0] # min_x snapped <= plane_min_x + assert pec_bounds[1][0] >= plane.bounds[1][0] # max_x snapped >= plane_max_x + assert pec_bounds[0][0] > sim_bounds[0][0] # NOT at sim edge + assert pec_bounds[1][0] < sim_bounds[1][0] # NOT at sim edge + + # z-axis: plane spans [-0.75, 0.75], sim spans [-2.0, 2.0] + assert pec_bounds[0][1] <= plane.bounds[0][2] + assert pec_bounds[1][1] >= plane.bounds[1][2] + assert pec_bounds[0][1] > sim_bounds[0][2] + assert pec_bounds[1][1] < sim_bounds[1][2] + + # PEC should be close to plane edges (within one grid cell) + assert abs(pec_bounds[0][0] - plane.bounds[0][0]) <= dl + assert abs(pec_bounds[1][0] - plane.bounds[1][0]) <= dl + assert abs(pec_bounds[0][1] - plane.bounds[0][2]) <= dl + assert abs(pec_bounds[1][1] - plane.bounds[1][2]) <= dl + + +def test_mode_solver_pec_boundary_at_sim_edge(): + """Test that PEC at simulation edge gives same behavior as before (plane spans full sim).""" + freq0 = td.C_0 / 1.55 + + sim_size = (2.0, 0.0, 1.5) + 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(), + z=td.Boundary.pec(), + ), + sources=[ + td.PointDipole( + center=(0, 0, 0), + source_time=td.GaussianPulse(freq0=freq0, fwidth=freq0 / 10), + polarization="Ex", + ) + ], + ) + + # Plane spans full simulation + plane = td.Box(center=(0, 0, 0), size=(sim_size[0], 0, sim_size[2])) + + ms = ModeSolver( + simulation=simulation, + plane=plane, + mode_spec=td.ModeSpec(num_modes=1), + freqs=[freq0], + direction="+", + ) + + pec_bounds = ms._pec_boundary_positions + sim_grid = simulation.grid.boundaries.to_list + + # PEC should be at (or very close to) simulation grid boundaries + assert np.isclose(pec_bounds[0][0], sim_grid[0][0], rtol=fp_eps) + assert np.isclose(pec_bounds[1][0], sim_grid[0][-1], rtol=fp_eps) + assert np.isclose(pec_bounds[0][1], sim_grid[2][0], rtol=fp_eps) + assert np.isclose(pec_bounds[1][1], sim_grid[2][-1], rtol=fp_eps) + + +def test_mode_solver_pec_boundary_1d_mode_solve(): + """Test that degenerate tangential axes (2D simulation) use Off behavior and mode solve completes. + + In a 2D simulation (zero size along z), the z-axis has num_cells <= 1. When z is tangential + to the mode plane (propagation along y), the PEC snapping should use SnapBehavior.Off for that + degenerate axis to avoid issues with snap_box_to_grid's force-expansion logic. + """ + freq0 = td.C_0 / 1.55 + + # 2D simulation: zero size in z → num_cells <= 1 along z + # Normal axis is y (propagation), tangential axes are x and z. + # z is degenerate (1 cell). + simulation = td.Simulation( + size=(2.0, 3.0, 0.0), + grid_spec=td.GridSpec.uniform(dl=0.05), + run_time=1e-14, + boundary_spec=td.BoundarySpec( + x=td.Boundary.pec(), + y=td.Boundary.pec(), + z=td.Boundary.periodic(), + ), + sources=[ + td.PointDipole( + center=(0, 0, 0), + source_time=td.GaussianPulse(freq0=freq0, fwidth=freq0 / 10), + polarization="Ex", + ) + ], + ) + + # Mode plane normal to y, tangential to x (non-degenerate) and z (degenerate). + # Plane z-size is nonzero (required for a valid planar Box), but the sim has 0 z-size + # so num_cells[z] <= 1 and SnapBehavior.Off will be used for z. + plane = td.Box(center=(0, 0, 0), size=(2.0, 0, 1.0)) + + ms = ModeSolver( + simulation=simulation, + plane=plane, + mode_spec=td.ModeSpec(num_modes=1), + freqs=[freq0], + direction="+", + ) + + # Should not raise — degenerate z-axis uses SnapBehavior.Off + pec_bounds = ms._pec_boundary_positions + assert pec_bounds is not None + + # The mode solve should complete without error + data = ms.data + assert data is not None diff --git a/tidy3d/components/microwave/path_integrals/mode_plane_analyzer.py b/tidy3d/components/microwave/path_integrals/mode_plane_analyzer.py index c5ab6c7e7c..7a070fb5ff 100644 --- a/tidy3d/components/microwave/path_integrals/mode_plane_analyzer.py +++ b/tidy3d/components/microwave/path_integrals/mode_plane_analyzer.py @@ -89,6 +89,34 @@ def _get_mode_limits( return (tuple(min_b_2d_list), max_b) + def _get_pec_boundary_positions( + self, sim_grid: Grid, mode_symmetry: tuple[Symmetry, Symmetry, Symmetry] + ) -> Bound: + """PEC boundary positions snapped to nearest grid boundary on or outside the plane. + + Uses ``SnapBehavior.Expand`` (no margin) so the PEC is placed at the grid + boundary closest to the mode plane edge. This matches the positions used by + the mode solver for PEC enforcement during eigensolve. + + For symmetry axes, the min bound is moved to the symmetry center so the + PEC line is placed where the half-domain boundary is. The filter method + ``_filter_conductors_touching_sim_bounds`` then removes that line for + PMC symmetry (symmetry=1) but keeps it for PEC-like symmetry (symmetry=-1). + """ + behavior = [SnapBehavior.Off] * 3 + location = [SnapLocation.Boundary] * 3 + for ax in range(3): + if ax != self._normal_axis and sim_grid.num_cells[ax] > 1: + behavior[ax] = SnapBehavior.Expand + snap_spec = SnappingSpec(location=location, behavior=behavior) + pec_box = snap_box_to_grid(sim_grid, self.geometry, snap_spec=snap_spec) + min_b, max_b = pec_box.bounds + min_b_list = list(min_b) + for dim in range(3): + if mode_symmetry[dim] != 0: + min_b_list[dim] = self.center[dim] + return (tuple(min_b_list), max_b) + def _get_isolated_conductors_as_shapely( self, plane: Box, @@ -125,7 +153,7 @@ def is_conductor(med: Medium) -> bool: def _filter_conductors_touching_sim_bounds( self, - mode_limits: Bound, + pec_bounds: Bound, mode_symmetry_3d: tuple[Symmetry, Symmetry, Symmetry], conductor_polygons: list[Shapely], ) -> list[Shapely]: @@ -135,8 +163,8 @@ def _filter_conductors_touching_sim_bounds( Parameters ---------- - mode_limits : Bound - The locations of the boundary conditions. + pec_bounds : Bound + PEC boundary positions snapped to the grid. mode_symmetry_3d : tuple[Symmetry, Symmetry, Symmetry] Symmetry settings for the mode solver plane. conductor_polygons : list[Shapely] @@ -148,7 +176,7 @@ def _filter_conductors_touching_sim_bounds( list[Shapely] The filtered list of shapely geometries, where structures "shorted" to PEC boundaries have been removed. """ - min_b_3d, max_b_3d = mode_limits[0], mode_limits[1] + min_b_3d, max_b_3d = pec_bounds[0], pec_bounds[1] _, mode_symmetry = Geometry.pop_axis(mode_symmetry_3d, self._normal_axis) _, min_b = Geometry.pop_axis(min_b_3d, self._normal_axis) _, max_b = Geometry.pop_axis(max_b_3d, self._normal_axis) @@ -219,8 +247,12 @@ def bounding_box_from_shapely(geom: Shapely) -> Box: intersection_plane, structures ) + # Use PEC-snapped bounds (tighter than mode_limits) for filtering conductors + # that are "shorted" to PEC boundaries. This matches the PEC positions used + # by the mode solver during eigensolve. + pec_bounds = self._get_pec_boundary_positions(grid, mode_symmetry_3d) filtered_conductor_shapely = self._filter_conductors_touching_sim_bounds( - (min_b_3d, max_b_3d), mode_symmetry_3d, isolated_conductor_shapely + pec_bounds, mode_symmetry_3d, isolated_conductor_shapely ) if len(filtered_conductor_shapely) < 1: @@ -243,7 +275,6 @@ def bounding_box_from_shapely(geom: Shapely) -> Box: box_snapped = snap_box_to_grid(grid, box, snap_spec) bounding_boxes.append(box_snapped) - # TODO Improve these checks once FXC-4112-PEC-boundary-position-not-respected-by-ModeSolver is merged for bounding_box in bounding_boxes: if self._check_box_intersects_with_conductors(isolated_conductor_shapely, bounding_box): raise SetupError( diff --git a/tidy3d/components/mode/mode_solver.py b/tidy3d/components/mode/mode_solver.py index 08a161a6b4..1a3866d708 100644 --- a/tidy3d/components/mode/mode_solver.py +++ b/tidy3d/components/mode/mode_solver.py @@ -16,7 +16,17 @@ Tidy3dBaseModel, cached_property, ) -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 ( ModeIndexDataArray, ScalarModeFieldCylindricalDataArray, @@ -81,6 +91,7 @@ Ax, Axis, Axis2D, + Bound2D, EpsSpecType, PlotScale, Symmetry, @@ -1201,6 +1212,108 @@ def _reduced_simulation_copy_with_fallback(self) -> ModeSolver: ) return self + @cached_property + def _pec_boundary_positions(self) -> Bound2D: + """PEC boundary positions snapped to nearest grid boundary on or outside the plane. + + For each tangential axis, the PEC boundary is placed at the grid boundary + closest to (but not inside) the mode plane bounds. This ensures correct PEC + placement for interior wave ports where the mode plane is smaller than the + simulation domain. + + Also logs warnings if the PEC position coincides with a simulation boundary + that uses an unsupported condition (PMC, Periodic, Bloch). + + Returns + ------- + Bound2D + ((min_0, min_1), (max_0, max_1)) positions along tangential axes. + """ + from tidy3d.components.geometry.utils import ( + SnapBehavior, + SnapLocation, + SnappingSpec, + snap_box_to_grid, + ) + + _, tangential_axes = Box.pop_axis([0, 1, 2], self.normal_axis) + + # Build snap spec: Expand on tangential axes, Off on normal/degenerate axes + behavior = [SnapBehavior.Off, SnapBehavior.Off, SnapBehavior.Off] + location = [SnapLocation.Boundary, SnapLocation.Boundary, SnapLocation.Boundary] + for ax in tangential_axes: + # Skip degenerate axes (1D mode solves in 2D simulations) — + # solver grid is fixed to all available cells, PEC truncation is meaningless + if self.simulation.grid.num_cells[ax] > 1: + behavior[ax] = SnapBehavior.Expand + + snap_spec = SnappingSpec( + location=tuple(location), + behavior=tuple(behavior), + ) + snapped = snap_box_to_grid(self.simulation.grid, self.plane, snap_spec) + + # Extract tangential-axis bounds as Bound2D + pos_min = [snapped.bounds[0][ax] for ax in tangential_axes] + pos_max = [snapped.bounds[1][ax] for ax in tangential_axes] + + # Warn if PEC position coincides with a simulation boundary using unsupported BC + 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 + + for axis in tangential_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] + + # Check if solver grid extends beyond (below/above) simulation boundary + intersects_min = solver_min < sim_min and not np.isclose( + solver_min, sim_min, rtol=fp_eps, atol=fp_eps + ) + 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." + ) + + # 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." + ) + + 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 @@ -1217,7 +1330,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 @@ -1291,7 +1405,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 @@ -1371,10 +1487,16 @@ def _colocate_data(self, mode_solver_data: ModeSolverData) -> ModeSolverData: colocate_coords = self._get_colocation_coordinates() + min_bound, max_bound = self._pec_boundary_positions + _, plane_dims = self.plane.pop_axis("xyz", self.normal_axis) + slice_dict = {plane_dims[axis]: slice(min_bound[axis], max_bound[axis]) for axis in [0, 1]} # Colocate input data to new coordinates data_dict_colocated = {} for key, field in mode_solver_data.symmetry_expanded.field_components.items(): - data_dict_colocated[key] = field.interp(**colocate_coords).astype(field.dtype) + field_sliced = field.sel(slice_dict) + data_dict_colocated[key] = field_sliced.interp( + **colocate_coords, kwargs={"fill_value": "extrapolate"} + ).astype(field.dtype) # Update data mode_solver_monitor = self.to_mode_solver_monitor( @@ -1601,14 +1723,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) @@ -1626,7 +1753,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 = [] @@ -1634,7 +1763,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) @@ -1692,6 +1824,7 @@ def _solve_single_freq( direction=self.direction, precision=self._precision, plane_center=self.plane_center_tangential(self.plane), + sim_pec_bound=self._pec_boundary_positions, ) fields = self._postprocess_solver_fields( @@ -1756,6 +1889,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._pec_boundary_positions, ) fields = self._postprocess_solver_fields( diff --git a/tidy3d/components/mode/solver.py b/tidy3d/components/mode/solver.py index ee6ac0eae7..e223bc490e 100644 --- a/tidy3d/components/mode/solver.py +++ b/tidy3d/components/mode/solver.py @@ -20,7 +20,7 @@ from scipy import sparse as sp - from tidy3d.components.types import EpsSpecType, ModeSolverType + from tidy3d.components.types import Bound2D, EpsSpecType, ModeSolverType # Consider vec to be complex if norm(vec.imag)/norm(vec) > TOL_COMPLEX TOL_COMPLEX = 1e-10 @@ -67,6 +67,7 @@ def compute_modes( direction: Literal["+", "-"] = "+", solver_basis_fields: Optional[ArrayComplex] = None, plane_center: Optional[tuple[float, float]] = None, + sim_pec_bound: Optional[Bound2D] = None, ) -> tuple[ArrayComplex, ArrayComplex, EpsSpecType]: """ Solve for the modes of a waveguide cross-section. @@ -106,6 +107,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 ------- @@ -124,6 +128,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 @@ -301,6 +366,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) @@ -1092,6 +1170,39 @@ def eigs_to_effective_index( raise RuntimeError(f"Unidentified 'mode_solver_type={mode_solver_type}'.") + @staticmethod + def _truncate_medium_data( + mat_data: Union[np.ndarray, tuple[np.ndarray, ...]], + cell_min: list[int], + cell_max: list[int], + ) -> Union[np.ndarray, tuple[np.ndarray, ...]]: + """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, np.ndarray): + # 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: Union[ArrayComplex, Sequence[ArrayComplex]], diff --git a/tidy3d/components/types/__init__.py b/tidy3d/components/types/__init__.py index b68e245b88..f661aba700 100644 --- a/tidy3d/components/types/__init__.py +++ b/tidy3d/components/types/__init__.py @@ -19,6 +19,7 @@ Axis, Axis2D, Bound, + Bound2D, BoxSurface, ClipOperationType, ColormapType, @@ -82,6 +83,7 @@ "Axis", "Axis2D", "Bound", + "Bound2D", "BoxSurface", "ClipOperationType", "ColormapType", diff --git a/tidy3d/components/types/base.py b/tidy3d/components/types/base.py index c27c490f00..6e4f70d5f9 100644 --- a/tidy3d/components/types/base.py +++ b/tidy3d/components/types/base.py @@ -245,6 +245,7 @@ def _parse_complex(v: Any) -> complex: CoordinateOptional = tuple[Optional[float], Optional[float], Optional[float]] Coordinate2D = tuple[float, float] Bound = tuple[Coordinate, Coordinate] +Bound2D = tuple[Coordinate2D, Coordinate2D] GridSize = Union[PositiveFloat, tuple[PositiveFloat, ...]] Axis = Annotated[Literal[0, 1, 2], BeforeValidator(_dtype2python)] Axis2D = Annotated[Literal[0, 1], BeforeValidator(_dtype2python)]