fix(mode): fixed bug where mode solver was not respecting location of PEC boundary conditions from the Simulation

Author: dmarek-flex

CHANGELOG.md

@@ -97,6 +97,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.

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,328 @@ 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 grid 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]
+
+    # Get solver grid boundaries (may extend beyond simulation)
+    solver_grid = ms._solver_grid
+    solver_x = solver_grid.boundaries.x
+    solver_z = solver_grid.boundaries.z
+
+    # Check if the solver grid extends beyond simulation bounds
+    grid_extends_x_min = solver_x[0] < sim_x_min - fp_eps
+    grid_extends_x_max = solver_x[-1] > sim_x_max + fp_eps
+    grid_extends_z_min = solver_z[0] < sim_z_min - fp_eps
+    grid_extends_z_max = solver_z[-1] > sim_z_max + fp_eps
+
+    # 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
+        if grid_extends_x_min:
+            x_outside_min = field_coords_x[field_coords_x < sim_x_min - 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))}"
+                )
+
+        if grid_extends_x_max:
+            x_outside_max = field_coords_x[field_coords_x > sim_x_max + 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
+        if grid_extends_z_min:
+            z_outside_min = field_coords_z[field_coords_z < sim_z_min - 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))}"
+                )
+
+        if grid_extends_z_max:
+            z_outside_max = field_coords_z[field_coords_z > sim_z_max + 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, atol=fp_eps)]
+        x_at_max = field_coords_x[np.isclose(field_coords_x, sim_x_max, 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, atol=fp_eps)]
+        z_at_max = field_coords_z[np.isclose(field_coords_z, sim_z_max, 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, atol=fp_eps)]
+    x_at_max = field_coords_x[np.isclose(field_coords_x, sim_x_max, 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, atol=fp_eps)]
+    z_at_max = field_coords_z[np.isclose(field_coords_z, sim_z_max, 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