bug fix: spatially varying plane sources and PMC walls under mirror symme…

src/fdtdx/fdtd/initialization.py

@@ -162,8 +162,9 @@ def place_objects(
     # non-symmetric path is unchanged.
     dropped_names: set[str] = set()
     reduced_volume_shape = None
+    symmetry_clip_low: dict[str, tuple[int, int, int]] = {}
     if config.has_symmetry:
-        resolved_slices, dropped_names, reduced_volume_shape = reduce_resolved_slices(
+        resolved_slices, dropped_names, reduced_volume_shape, symmetry_clip_low = reduce_resolved_slices(
             resolved_slices=resolved_slices,
             object_map=object_map,
             config=config,
@@ -201,25 +202,26 @@ def place_objects(
         obj = object_map[name]
         assert key is not None
         key, subkey = jax.random.split(key)
-        placed_objects.append(
-            obj.place_on_grid(
-                grid_slice_tuple=slice_tuple,
-                config=config,
-                key=subkey,
-            )
+        placed = obj.place_on_grid(
+            grid_slice_tuple=slice_tuple,
+            config=config,
+            key=subkey,
         )
+        if name in symmetry_clip_low:
+            placed = placed.aset("_symmetry_clip_low", symmetry_clip_low[name])
+        placed_objects.append(placed)
 
     # Step 7: Place volume first (index 0)
     assert key is not None
     key, subkey = jax.random.split(key)
-    placed_objects.insert(
-        0,
-        volume_obj.place_on_grid(
-            grid_slice_tuple=resolved_slices[volume_obj.name],
-            config=config,
-            key=subkey,
-        ),
+    placed_volume = volume_obj.place_on_grid(
+        grid_slice_tuple=resolved_slices[volume_obj.name],
+        config=config,
+        key=subkey,
     )
+    if volume_obj.name in symmetry_clip_low:
+        placed_volume = placed_volume.aset("_symmetry_clip_low", symmetry_clip_low[volume_obj.name])
+    placed_objects.insert(0, placed_volume)
 
     # Step 8: Insert the PEC/PMC symmetry walls and forward the per-axis condition to mode
     # sources/detectors, then warn that the simulation now runs on the reduced domain.

src/fdtdx/fdtd/symmetry.py

@@ -30,7 +30,6 @@
 from fdtdx.objects.boundaries.bloch import BlochBoundary
 from fdtdx.objects.boundaries.boundary import BaseBoundary
 from fdtdx.objects.boundaries.pec import PerfectElectricConductor
-from fdtdx.objects.boundaries.pmc import PerfectMagneticConductor
 from fdtdx.objects.detectors.diffractive import DiffractiveDetector
 from fdtdx.objects.detectors.energy import EnergyDetector
 from fdtdx.objects.detectors.field import FieldDetector
@@ -105,7 +104,7 @@ def reduce_resolved_slices(
     object_map: dict[str, SimulationObject],
     config: SimulationConfig,
     volume_name: str,
-) -> tuple[dict[str, SliceTuple3D], set[str], tuple[int, int, int]]:
+) -> tuple[dict[str, SliceTuple3D], set[str], tuple[int, int, int], dict[str, tuple[int, int, int]]]:
     """Clip full-domain grid slices onto the symmetric (kept upper) half along each axis.
 
     Args:
@@ -115,9 +114,10 @@ def reduce_resolved_slices(
         volume_name (str): Name of the simulation volume object.
 
     Returns:
-        tuple: ``(new_slices, dropped_names, reduced_volume_shape)`` where ``new_slices`` excludes
-        the dropped objects, ``dropped_names`` lists objects entirely in the discarded half, and
-        ``reduced_volume_shape`` is the reduced volume grid shape ``(Nx', Ny', Nz')``.
+        tuple: ``(new_slices, dropped_names, reduced_volume_shape, clip_low)`` where ``new_slices``
+        excludes the dropped objects, ``dropped_names`` lists objects entirely in the discarded half,
+        ``reduced_volume_shape`` is the reduced volume grid shape ``(Nx', Ny', Nz')``, and
+        ``clip_low`` maps object name → cells removed from its low side per axis.
     """
     symmetry = config.symmetry
     vol_slice = resolved_slices[volume_name]
@@ -145,20 +145,25 @@ def reduce_resolved_slices(
 
     new_slices: dict[str, SliceTuple3D] = {}
     dropped: set[str] = set()
+    clip_low: dict[str, tuple[int, int, int]] = {}
     for name, sl in resolved_slices.items():
         if name == volume_name:
             new_slices[name] = (new_vol[0], new_vol[1], new_vol[2])
+            vol_cut = [mid_abs[a] - vol_slice[a][0] if symmetry[a] != 0 else 0 for a in range(3)]
+            clip_low[name] = (vol_cut[0], vol_cut[1], vol_cut[2])
             continue
 
         obj = object_map[name]
         clipped: list[tuple[int, int]] = []
+        cut: list[int] = []
         drop = False
         for a in range(3):
             s0, s1 = sl[a]
             if symmetry[a] != 0:
                 m = mid_abs[a]
                 ns0 = max(s0, m) - m
                 ns1 = min(s1, vol_slice[a][1]) - m
+                cut.append(max(0, min(m, s1) - s0))
                 if ns1 <= ns0:
                     drop = True
                 elif s0 < m < s1 and (s0 + s1) != 2 * m:
@@ -170,6 +175,7 @@ def reduce_resolved_slices(
                     )
                 clipped.append((ns0, ns1))
             else:
+                cut.append(0)
                 clipped.append((s0, s1))
 
         if drop:
@@ -189,8 +195,9 @@ def reduce_resolved_slices(
             continue
 
         new_slices[name] = (clipped[0], clipped[1], clipped[2])
+        clip_low[name] = (cut[0], cut[1], cut[2])
 
-    return new_slices, dropped, reduced_volume_shape
+    return new_slices, dropped, reduced_volume_shape, clip_low
 
 
 def make_symmetry_walls(
@@ -199,7 +206,16 @@ def make_symmetry_walls(
     key: jax.Array,
     existing_names: set[str],
 ) -> list[SimulationObject]:
-    """Create the PEC/PMC wall objects sitting on each symmetry plane (reduced min edge).
+    """Create the PEC wall objects sitting on each electric symmetry plane (reduced min edge).
+
+    Only ``symmetry[a] == -1`` (PEC / electric mirror) gets an explicit wall. The symmetry plane
+    lies on the reduced domain's low cell *edge*, i.e. at local coordinate ``-0.5``. Along the
+    axis, tangential ``E`` has Yee offset 0, so its first node sits half a cell inside the domain
+    and the electric mirror has to be imposed explicitly. Tangential ``H`` has offset ``1/2``, so
+    the node that a magnetic mirror forces to zero sits at local ``-1`` — outside the array, where
+    zero-padding already supplies it. A ``PerfectMagneticConductor`` object would instead zero
+    tangential ``H`` at local ``0``, a full cell inside the domain, over-constraining the field
+    (measured ~30% error near the plane versus ~0.1% with the natural termination).
 
     Args:
         config (SimulationConfig): Simulation config; ``config.symmetry`` selects the axes/types.
@@ -208,12 +224,13 @@ def make_symmetry_walls(
         existing_names (set[str]): Names already in use, so wall names can be made unique.
 
     Returns:
-        list[SimulationObject]: PEC/PMC boundary objects, already placed on the reduced grid.
+        list[SimulationObject]: PEC boundary objects, already placed on the reduced grid. PMC
+        (``+1``) axes contribute nothing — their mirror is the zero-padded array edge.
     """
     walls: list[SimulationObject] = []
     used = set(existing_names)
     for a in range(3):
-        if config.symmetry[a] == 0:
+        if config.symmetry[a] != -1:
             continue
         name = f"_sym_wall_{_AXIS_NAMES[a]}"
         suffix = 0
@@ -227,8 +244,7 @@ def make_symmetry_walls(
             1 if a == 1 else None,
             1 if a == 2 else None,
         )
-        wall_cls = PerfectElectricConductor if config.symmetry[a] == -1 else PerfectMagneticConductor
-        wall = wall_cls(
+        wall = PerfectElectricConductor(
             name=name,
             axis=a,
             partial_grid_shape=partial_grid_shape,

src/fdtdx/objects/object.py

@@ -245,6 +245,15 @@ class SimulationObject(TreeClass, ABC):
     )
     _config: SimulationConfig = private_field()
 
+    #: Cells removed from this object's low side per axis by mirror-symmetry reduction.
+    _symmetry_clip_low: tuple[int, int, int] = frozen_private_field(default=(0, 0, 0))
+
+    @property
+    def unreduced_grid_shape(self) -> GridShape3D:
+        """This object's grid shape before symmetry reduction clipped it."""
+        shape, clip = self.grid_shape, self._symmetry_clip_low
+        return (shape[0] + clip[0], shape[1] + clip[1], shape[2] + clip[2])
+
     @property
     def grid_slice_tuple(self) -> SliceTuple3D:
         if self._grid_slice_tuple == INVALID_SLICE_TUPLE_3D:

src/fdtdx/objects/sources/linear_polarization.py

@@ -203,6 +203,9 @@ def apply(
         u_basis = h_axis - jnp.dot(h_axis, wave_vector) * wave_vector
         u_basis = u_basis / jnp.linalg.norm(u_basis)
         v_basis = jnp.cross(wave_vector, u_basis)
+        # Orient along +vertical_axis; the raw cross product sign flips the sampled
+        # profile about ``center`` for some axis/direction combinations.
+        v_basis = v_basis * jnp.where(v_basis[self.vertical_axis] >= 0, 1.0, -1.0).astype(v_basis.dtype)
 
         # projection
         def project(point):
@@ -241,7 +244,8 @@ def project(point):
                 inv_permittivity=inv_permittivities,
                 inv_permeability=inv_permeabilities,
             )
-            total_energy_root = jnp.sqrt(energy.sum())
+            # Under symmetry the plane is clipped; normalize by the full-domain energy.
+            total_energy_root = jnp.sqrt(energy.sum() * self._symmetry_fold_factor())
             E = E / total_energy_root
             H = H / total_energy_root
 
@@ -366,7 +370,7 @@ def _get_amplitude_raw(
             h_widths = horizontal_edges[1:] - horizontal_edges[:-1]
             v_widths = vertical_edges[1:] - vertical_edges[:-1]
             cell_areas = h_widths[:, None] * v_widths[None, :]
-            profile_2d = profile_2d / (profile_2d * cell_areas).sum()
+            profile_2d = profile_2d / ((profile_2d * cell_areas).sum() * self._symmetry_fold_factor())
             return jnp.expand_dims(profile_2d, axis=self.propagation_axis)
 
         grid_radius = self.radius / self._config.uniform_spacing()
@@ -378,7 +382,8 @@ def _get_amplitude_raw(
             radii=(grid_radius, grid_radius),
             std=self.std,
         )
-        return profile
+        # ``_gauss_profile`` normalizes over its own cells; rescale to the full-domain sum.
+        return profile / self._symmetry_fold_factor()
 
 
 @autoinit

src/fdtdx/objects/sources/tfsf.py

@@ -527,29 +527,41 @@ def _get_azimuth_elevation(
         )
         return azimuth_radians, elevation_radians
 
+    def _symmetry_fold_factor(self) -> int:
+        """Factor relating a transverse sum over the reduced plane to the full-domain sum."""
+        clip = self._symmetry_clip_low
+        return 2 ** sum(1 for axis in (self.horizontal_axis, self.vertical_axis) if clip[axis] > 0)
+
     def _get_center(self, key: jax.Array) -> jax.Array:  # shape(2,)
         """Calculate the randomized source center.
 
         Uniform-grid sources use the legacy index-space center.  On a
         non-uniform grid, tilted projections and Gaussian profiles are sampled in
         physical coordinates, so the returned center is measured in metres from
         the source-slice lower edge along the transverse axes.
+
+        Under mirror symmetry the source is clipped to the kept half, so the center is
+        the *full-domain* center expressed in reduced-local coordinates: it lands on the
+        symmetry plane at the slice's lower edge rather than mid-slice.
         """
+        clip = self._symmetry_clip_low
         if self._config.has_nonuniform_grid:
             grid = self._config.resolved_grid
             assert grid is not None
-            local_edges = []
+            local_centers = []
             for axis in (self.horizontal_axis, self.vertical_axis):
                 lower, upper = self.grid_slice_tuple[axis]
                 edges = grid.edges(axis)[lower : upper + 1]
-                local_edges.append(edges - edges[0])
-            center_horizontal = 0.5 * local_edges[0][-1]
-            center_vertical = 0.5 * local_edges[1][-1]
+                edges = edges - edges[0]
+                # A clipped axis has the symmetry plane at its lower edge (local 0.0).
+                local_centers.append(jnp.asarray(0.0) if clip[axis] > 0 else 0.5 * edges[-1])
+            center_horizontal, center_vertical = local_centers
         else:
-            horizontal_size = self.grid_shape[self.horizontal_axis]
-            vertical_size = self.grid_shape[self.vertical_axis]
-            center_horizontal = (horizontal_size - 1) / 2
-            center_vertical = (vertical_size - 1) / 2
+            index_centers = []
+            for axis in (self.horizontal_axis, self.vertical_axis):
+                full_size = self.grid_shape[axis] + clip[axis]
+                index_centers.append((full_size - 1) / 2 - clip[axis])
+            center_horizontal, center_vertical = index_centers
 
         key, subkey = jax.random.split(key)
         horizontal_offset = jax.random.uniform(