From 340d74fb708e345f5aa52dc85dbfa9929cbe37f6 Mon Sep 17 00:00:00 2001 From: bruxillensis Date: Wed, 29 Jul 2026 18:37:58 -0400 Subject: [PATCH] fix: spatially varying plane sources and PMC walls under mirror symmetry (#425) Four independent defects made a symmetry-reduced simulation wrong for any source with spatial variation along a symmetry axis. All were invisible to the existing symmetry tests, which only use UniformPlaneSource with normalize_by_energy=False and direction="+" -- a combination blind to every one of them. 1. Profile center re-centered on the reduced quadrant. TFSFPlaneSource._get_center derived the transverse center from the object's clipped extent, so a Gaussian spot moved to the middle of the kept half instead of staying on the symmetry plane. The reduction discarded the pre-clip extent, so the source could not know it had been cut; reduce_resolved_slices now reports the per-axis low-side clip and place_objects records it as SimulationObject._symmetry_clip_low. 2. Normalization summed over the clipped plane. normalize_by_energy and the Gaussian profile's own sum normalization both divided by a quarter-plane sum, inflating the injected amplitude by 2**(n_axes/2) -- exactly 2x for two symmetric axes. 3. Profile mirrored about `center` along the vertical axis. v_basis = cross(wave_vector, u_basis) has a sign that depends on the propagation axis and direction, so the sampled profile was reflected for e.g. direction="-" along z (and, oppositely, direction="+" along y). Predates the symmetry feature by ~14 months and is invisible for any profile symmetric about its own center; it also blocked fix 1, since a correctly off-center profile got flipped out of the array. 4. The PMC symmetry wall over-constrained the field. Along a symmetry axis the six components split by Yee offset: tangential E and normal H sit at offset 0, normal E and tangential H at offset 1/2. The mirror plane lands on the reduced domain's low cell edge, so the tangential-H node a magnetic mirror must zero is at local -1 -- outside the array, where zero-padding already supplies it. Inserting a PerfectMagneticConductor instead zeroed tangential H at local 0, a full cell inside the domain. PEC still needs its wall (tangential E has no node on the plane). make_symmetry_walls now builds walls only for electric (-1) planes. Validation. A symmetry-reduced run is compared cell-for-cell against the half of an equivalent full-domain run that it is supposed to represent (no unfolding involved). The figure quoted is the maximum absolute deviation of E over that half, divided by the peak |E| of the full-domain run -- so 0.0019 means the worst cell in the reduced domain is off by 0.19% of the full-domain field peak. Source is a Gaussian beam except where a dipole group is named. Values are before -> after this change: PEC on x/y/z (single axis) 0.0019 -> 0.0019 (unchanged) PMC on x/y/z (single axis) 0.2768 -> 0.0019 PEC+PMC pairs (4 orientations) 0.2753 -> 0.0019 PMC on x/y (mirror dipole pair) 0.6819 -> 0.0002 PMC+PMC (mirror dipole quad) 1.8677 -> 0.0003 PEC-only axes are untouched by this commit; they are listed for contrast, since they were already correct and are what made the PMC error stand out. The residual ~0.2% on the beam cases is dominated by PML reflection and beam truncation, not by the symmetry plane -- the dipole cases, which have neither, land at 0.02-0.03%. The PMC fix converges under grid refinement -- 0.00142, 0.00053, 0.00032 at 200, 100 and 50 nm spacing -- so it is a consistent discretization rather than a cancellation that happens to work at one resolution. With all four applied the injected reduced profile equals the kept half of the full-domain profile to float32 roundoff (~1e-6 relative). Known remaining limitation, not addressed here: unfold_fields applies one array flip to all three components, but the correct reconstruction map differs by Yee staggering class. Measured over the PML-free interior, three components per axis reconstruct to ~0.1% of the full-domain peak and the other three to 2-12%, and which class is exact flips between axes. The reduced simulation itself is correct; only the unfolded reconstruction is affected. --- src/fdtdx/fdtd/initialization.py | 30 ++-- src/fdtdx/fdtd/symmetry.py | 38 ++-- src/fdtdx/objects/object.py | 9 + .../objects/sources/linear_polarization.py | 11 +- src/fdtdx/objects/sources/tfsf.py | 28 ++- .../fdtd/test_symmetry_plane_source.py | 147 ++++++++++++++++ .../fdtd/test_symmetry_reduction.py | 12 +- .../test_symmetry_gaussian_source.py | 166 ++++++++++++++++++ 8 files changed, 400 insertions(+), 41 deletions(-) create mode 100644 tests/integration/fdtd/test_symmetry_plane_source.py create mode 100644 tests/simulation/physics/boundaries/test_symmetry_gaussian_source.py diff --git a/src/fdtdx/fdtd/initialization.py b/src/fdtdx/fdtd/initialization.py index 041d7344..70d45647 100644 --- a/src/fdtdx/fdtd/initialization.py +++ b/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. diff --git a/src/fdtdx/fdtd/symmetry.py b/src/fdtdx/fdtd/symmetry.py index f7101586..9f76f45d 100644 --- a/src/fdtdx/fdtd/symmetry.py +++ b/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,13 +145,17 @@ 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] @@ -159,6 +163,7 @@ def reduce_resolved_slices( 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, diff --git a/src/fdtdx/objects/object.py b/src/fdtdx/objects/object.py index e99c16f6..e4291a56 100644 --- a/src/fdtdx/objects/object.py +++ b/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: diff --git a/src/fdtdx/objects/sources/linear_polarization.py b/src/fdtdx/objects/sources/linear_polarization.py index 1e159841..22582b47 100644 --- a/src/fdtdx/objects/sources/linear_polarization.py +++ b/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 diff --git a/src/fdtdx/objects/sources/tfsf.py b/src/fdtdx/objects/sources/tfsf.py index efefcb50..dccc4e68 100644 --- a/src/fdtdx/objects/sources/tfsf.py +++ b/src/fdtdx/objects/sources/tfsf.py @@ -527,6 +527,11 @@ 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. @@ -534,22 +539,29 @@ def _get_center(self, key: jax.Array) -> jax.Array: # shape(2,) 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( diff --git a/tests/integration/fdtd/test_symmetry_plane_source.py b/tests/integration/fdtd/test_symmetry_plane_source.py new file mode 100644 index 00000000..30ba5093 --- /dev/null +++ b/tests/integration/fdtd/test_symmetry_plane_source.py @@ -0,0 +1,147 @@ +"""Regression tests for spatially varying plane sources under mirror symmetry (issue #425). + +Three defects were fixed: + +1. ``TFSFPlaneSource._get_center`` computed the transverse profile center from the + *clipped* extent, so a Gaussian spot re-centered on the reduced quadrant instead of + staying on the symmetry plane. +2. ``normalize_by_energy`` (and the Gaussian profile's own sum normalization) summed over + the clipped plane, inflating the injected amplitude by ``2**(n_axes/2)``. +3. ``v_basis = cross(wave_vector, u_basis)`` mirrored the sampled profile about ``center`` + along the vertical axis for some axis/direction combinations. Predates symmetry; only + observable for a profile that is not symmetric about its own center. +""" + +import jax +import jax.numpy as jnp +import numpy as np +import pytest + +import fdtdx + +_L = 2e-6 +_H = 3e-6 +_SPACING = 200e-9 +_WAVELENGTH = 1.5e-6 + + +def _build(symmetry, *, direction="-", normalize=False, source_cls=None, **source_kwargs): + config = fdtdx.SimulationConfig( + time=20e-15, + dtype=jnp.float32, + courant_factor=0.99, + grid=fdtdx.UniformGrid(spacing=_SPACING), + symmetry=symmetry, + ) + object_list, constraints = [], [] + volume = fdtdx.SimulationVolume( + partial_real_shape=(2 * _L, 2 * _L, _H), + material=fdtdx.Material(permittivity=1.0, permeability=1.0), + ) + object_list.append(volume) + + cls = source_cls or fdtdx.GaussianPlaneSource + src = cls( + name="Source", + partial_grid_shape=(None, None, 1), + fixed_E_polarization_vector=(1, 0, 0), + wave_character=fdtdx.WaveCharacter(wavelength=_WAVELENGTH), + direction=direction, + normalize_by_energy=normalize, + **source_kwargs, + ) + object_list.append(src) + own = 1 if direction == "-" else -1 + margin = -0.6e-6 if direction == "-" else 0.6e-6 + constraints.extend( + [ + src.place_relative_to( + volume, axes=(2,), own_positions=(own,), other_positions=(own,), margins=(margin,) + ), + src.place_at_center(volume, axes=(0, 1)), + src.size_relative_to(volume, axes=(0, 1), proportions=(1.0, 1.0)), + ] + ) + + bound_cfg = fdtdx.BoundaryConfig.from_uniform_bound(boundary_type="pml", thickness=5) + bound_dict, c_list = fdtdx.boundary_objects_from_config(bound_cfg, volume) + constraints.extend(c_list) + object_list.extend(bound_dict.values()) + + key = jax.random.PRNGKey(0) + objs, arrays, params, config, _ = fdtdx.place_objects( + object_list=object_list, config=config, constraints=constraints, key=key + ) + arrays, objs, _ = fdtdx.apply_params(arrays, objs, params, key=key) + return objs, arrays, config + + +def _injected_profile(objs): + """The transverse |E| profile the source injects, squeezed to 2D.""" + src = next(s for s in objs.sources if s.name == "Source") + return np.abs(np.asarray(src._E)[0, :, :, 0]), src + + +class _OffCenterSource(fdtdx.UniformPlaneSource): + """Single bright cell at a deliberately off-center transverse index.""" + + def _get_amplitude_raw(self, center): + del center + idx = [0, 0, 0] + idx[self.horizontal_axis] = 2 + idx[self.vertical_axis] = 3 + return jnp.zeros(self.grid_shape).at[tuple(idx)].set(1.0) + + +class TestProfileOrientation: + """Fix 3: the sampled profile must not be mirrored about ``center``.""" + + @pytest.mark.parametrize("direction", ["+", "-"]) + def test_off_center_profile_is_not_mirrored(self, direction): + objs, _, _ = _build((0, 0, 0), direction=direction, source_cls=_OffCenterSource, amplitude=1.0) + prof, _ = _injected_profile(objs) + peak = np.unravel_index(np.argmax(prof), prof.shape) + assert peak == (2, 3), f"profile mirrored for direction={direction!r}: peak at {peak}" + assert np.count_nonzero(prof > 1e-12) == 1 + + +class TestSymmetryProfileCentering: + """Fixes 1 and 2: the reduced profile must equal the kept half of the full-domain profile.""" + + @pytest.mark.parametrize("symmetry", [(-1, 0, 0), (0, 1, 0), (-1, 1, 0)]) + @pytest.mark.parametrize("normalize", [False, True]) + def test_reduced_profile_matches_full_half(self, symmetry, normalize): + full_objs, _, _ = _build((0, 0, 0), normalize=normalize, radius=_L, std=1 / 3) + red_objs, _, _ = _build(symmetry, normalize=normalize, radius=_L, std=1 / 3) + full, _ = _injected_profile(full_objs) + red, _ = _injected_profile(red_objs) + + half = full + if symmetry[0] != 0: + half = half[half.shape[0] // 2 :, :] + if symmetry[1] != 0: + half = half[:, half.shape[1] // 2 :] + + assert red.shape == half.shape + np.testing.assert_allclose(red, half, rtol=0, atol=1e-5 * full.max()) + + @pytest.mark.parametrize("symmetry", [(-1, 0, 0), (0, 1, 0), (-1, 1, 0)]) + def test_profile_peaks_on_the_symmetry_plane(self, symmetry): + objs, _, _ = _build(symmetry, radius=_L, std=1 / 3) + prof, _ = _injected_profile(objs) + peak = np.unravel_index(np.argmax(prof), prof.shape) + for a, arr_axis in ((0, 0), (1, 1)): + if symmetry[a] != 0: + assert peak[arr_axis] == 0, f"peak {peak} is off the {a}-symmetry plane" + + def test_clip_low_recorded(self): + objs, _, _ = _build((-1, 1, 0), radius=_L, std=1 / 3) + _, src = _injected_profile(objs) + assert src._symmetry_clip_low == (10, 10, 0) + assert src.unreduced_grid_shape == (20, 20, 1) + + def test_clip_low_zero_without_symmetry(self): + objs, _, _ = _build((0, 0, 0), radius=_L, std=1 / 3) + _, src = _injected_profile(objs) + assert src._symmetry_clip_low == (0, 0, 0) + assert src.unreduced_grid_shape == src.grid_shape diff --git a/tests/integration/fdtd/test_symmetry_reduction.py b/tests/integration/fdtd/test_symmetry_reduction.py index 6f46e217..b2410b22 100644 --- a/tests/integration/fdtd/test_symmetry_reduction.py +++ b/tests/integration/fdtd/test_symmetry_reduction.py @@ -105,14 +105,16 @@ def test_symmetry_reduces_volume_and_inserts_walls(): assert oc.volume.grid_shape == (20, 10, 10) assert arrays.fields.E.shape == (3, 20, 10, 10) - # PEC wall on y, PMC wall on z, each a single cell on the symmetry plane (min edge). + # PEC wall on y, a single cell on the symmetry plane (min edge). pec = oc.pec_objects - pmc = oc.pmc_objects - assert len(pec) == 1 and len(pmc) == 1 + assert len(pec) == 1 assert pec[0].axis == 1 and pec[0].direction == "-" assert pec[0]._grid_slice_tuple[1] == (0, 1) - assert pmc[0].axis == 2 and pmc[0].direction == "-" - assert pmc[0]._grid_slice_tuple[2] == (0, 1) + + # No wall on the PMC z plane: the tangential-H node a magnetic mirror zeroes lies at local + # -1, outside the reduced array, where zero-padding already supplies it. A PMC object would + # zero tangential H a full cell inside the domain instead. + assert len(oc.pmc_objects) == 0 def test_symmetry_clips_and_drops_objects(): diff --git a/tests/simulation/physics/boundaries/test_symmetry_gaussian_source.py b/tests/simulation/physics/boundaries/test_symmetry_gaussian_source.py new file mode 100644 index 00000000..cfe7656f --- /dev/null +++ b/tests/simulation/physics/boundaries/test_symmetry_gaussian_source.py @@ -0,0 +1,166 @@ +"""A symmetry-reduced Gaussian beam, unfolded, must reproduce the full-domain run (issue #425). + +The existing symmetry simulation tests all use ``UniformPlaneSource`` with +``normalize_by_energy=False`` and ``direction="+"`` — a combination that cannot observe a +misplaced profile center, an inflated normalization, or a mirrored profile. This exercises +a spatially varying source instead. + +A magnetic (PMC) symmetry plane deliberately gets no wall object: see +``test_reduced_simulation_matches_kept_half``. +""" + +import jax +import jax.numpy as jnp +import numpy as np +import pytest + +import fdtdx +import fdtdx.fdtd.symmetry + +_L = 2e-6 +_H = 3e-6 +_SPACING = 200e-9 +_WAVELENGTH = 1.5e-6 +_TIME = 30e-15 + + +def _run(symmetry, *, polarization, direction="-", normalize=True): + config = fdtdx.SimulationConfig( + time=_TIME, + dtype=jnp.float32, + courant_factor=0.99, + grid=fdtdx.UniformGrid(spacing=_SPACING), + symmetry=symmetry, + ) + object_list, constraints = [], [] + volume = fdtdx.SimulationVolume( + partial_real_shape=(2 * _L, 2 * _L, _H), + material=fdtdx.Material(permittivity=1.0, permeability=1.0), + ) + object_list.append(volume) + + src = fdtdx.GaussianPlaneSource( + name="Source", + partial_grid_shape=(None, None, 1), + fixed_E_polarization_vector=polarization, + wave_character=fdtdx.WaveCharacter(wavelength=_WAVELENGTH), + direction=direction, + radius=_L, + std=1 / 3, + normalize_by_energy=normalize, + ) + object_list.append(src) + own = 1 if direction == "-" else -1 + margin = -0.6e-6 if direction == "-" else 0.6e-6 + constraints.extend( + [ + src.place_relative_to(volume, axes=(2,), own_positions=(own,), other_positions=(own,), margins=(margin,)), + src.place_at_center(volume, axes=(0, 1)), + src.size_relative_to(volume, axes=(0, 1), proportions=(1.0, 1.0)), + ] + ) + + bound_cfg = fdtdx.BoundaryConfig.from_uniform_bound(boundary_type="pml", thickness=5) + bound_dict, c_list = fdtdx.boundary_objects_from_config(bound_cfg, volume) + constraints.extend(c_list) + object_list.extend(bound_dict.values()) + + key = jax.random.PRNGKey(0) + objs, arrays, params, config, _ = fdtdx.place_objects( + object_list=object_list, config=config, constraints=constraints, key=key + ) + arrays, objs, _ = fdtdx.apply_params(arrays, objs, params, key=key) + _, arrays = fdtdx.run_fdtd(arrays=arrays, objects=objs, config=config, key=key) + return arrays + + +def _interior_error(reference, candidate): + """Max |diff| over the PML-free interior, relative to the reference peak.""" + nx, ny, nz = reference.shape[1:] + sl = (slice(8, nx - 8), slice(8, ny - 8), slice(6, nz - 6)) + a, b = reference[0][sl], candidate[0][sl] + return float(np.abs(a - b).max() / np.abs(a).max()) + + +@pytest.mark.parametrize("direction", ["-", "+"]) +@pytest.mark.parametrize("normalize", [True, False]) +def test_pec_symmetry_gaussian_matches_full_domain(direction, normalize): + """An x-polarized Gaussian beam with a PEC x-symmetry plane, unfolded, matches full domain.""" + pol = (1, 0, 0) + symmetry = (-1, 0, 0) + full = np.asarray(_run((0, 0, 0), polarization=pol, direction=direction, normalize=normalize).fields.E) + reduced = _run(symmetry, polarization=pol, direction=direction, normalize=normalize).fields.E + + # Half the domain in x, so the amplitude must not pick up a normalization factor. + assert reduced.shape[1] == full.shape[1] // 2 + unfolded = np.asarray(fdtdx.unfold_fields(reduced, symmetry, "E")) + assert unfolded.shape == full.shape + assert _interior_error(full, unfolded) < 0.02 + + +def test_pec_symmetry_preserves_amplitude_scale(): + """normalize_by_energy must not inflate the reduced-domain amplitude (fix 2).""" + pol = (1, 0, 0) + symmetry = (-1, 0, 0) + full = np.asarray(_run((0, 0, 0), polarization=pol).fields.E) + unfolded = np.asarray(fdtdx.unfold_fields(_run(symmetry, polarization=pol).fields.E, symmetry, "E")) + nx, ny, nz = full.shape[1:] + sl = (slice(8, nx - 8), slice(8, ny - 8), slice(6, nz - 6)) + a, b = full[0][sl], unfolded[0][sl] + scale = float(np.sum(a * b) / np.sum(b * b)) + assert abs(scale - 1.0) < 0.05, f"amplitude off by {scale:.4f}x" + + +def _kept_half(field, symmetry): + """The half of a full-domain array that the reduction keeps.""" + out = field + for a in range(3): + if symmetry[a] != 0: + half = out.shape[1 + a] // 2 + out = np.take(out, np.arange(half, out.shape[1 + a]), axis=1 + a) + return out + + +@pytest.mark.parametrize( + "symmetry,polarization", + [ + ((0, 1, 0), (1, 0, 0)), # PMC on y (x-polarized beam) + ((1, 0, 0), (0, 1, 0)), # PMC on x (y-polarized beam) + ((-1, 1, 0), (1, 0, 0)), # PEC x + PMC y + ((1, -1, 0), (0, 1, 0)), # PMC x + PEC y + ], +) +def test_reduced_simulation_matches_kept_half(symmetry, polarization): + """The reduced run must reproduce the half of the full domain it represents. + + A PMC symmetry plane gets no explicit wall: the tangential-H node a magnetic mirror + zeroes sits outside the reduced array, where zero-padding already supplies it. Inserting + a PerfectMagneticConductor instead zeroed tangential H a full cell inside the domain, + which showed up as ~30% error concentrated at the plane. + """ + full = np.asarray(_run((0, 0, 0), polarization=polarization).fields.E) + reduced = np.asarray(_run(symmetry, polarization=polarization).fields.E) + kept = _kept_half(full, symmetry) + assert reduced.shape == kept.shape + err = np.abs(kept - reduced).max() / np.abs(full).max() + assert err < 0.02, f"reduced run deviates from the kept half by {err:.4f}" + + +def test_no_pmc_wall_is_created_for_a_magnetic_symmetry_plane(): + from fdtdx.objects.boundaries.pmc import PerfectMagneticConductor + + config = fdtdx.SimulationConfig( + time=_TIME, + dtype=jnp.float32, + courant_factor=0.99, + grid=fdtdx.UniformGrid(spacing=_SPACING), + symmetry=(-1, 1, 0), + ) + walls = fdtdx.fdtd.symmetry.make_symmetry_walls( + config=config, + reduced_volume_shape=(8, 8, 8), + key=jax.random.PRNGKey(0), + existing_names=set(), + ) + assert not any(isinstance(w, PerfectMagneticConductor) for w in walls) + assert len(walls) == 1 # only the PEC x plane