[tests] verify hybrid mode signs against matched CPML scattering references

meanas/test/test_waveguide_fdtd_fdfd.py

@@ -36,7 +36,7 @@ SCATTERING_SOURCE_SLICES = (slice(4, 5), slice(None), slice(None))
 SCATTERING_REFLECT_SLICES = (slice(10, 11), slice(None), slice(None))
 SCATTERING_TRANSMIT_SLICES = (slice(29, 30), slice(None), slice(None))
 SCATTERING_STEP_X = 18
-SCATTERING_WARMUP_PERIODS = 10
+SCATTERING_WARMUP_PERIODS = 30
 SCATTERING_ACCUMULATION_PERIODS = 10
 
 
@@ -221,6 +221,21 @@ def _build_cpml_params() -> list[list[dict[str, numpy.ndarray | float]]]:
     ]
 
 
+def _build_cpml_frequency_dxes(base_dxes: list[list[numpy.ndarray]]) -> list[list[numpy.ndarray]]:
+    """Match the discrete CPML recurrence, including its Yee sample offsets."""
+    widths = [[dx.astype(complex) for dx in family] for family in base_dxes]
+    for axis, faces in enumerate(_build_cpml_params()):
+        for face in faces:
+            # psi[n] = p0 psi[n-1] + p1 p2 D[n]. Its phasor response adds
+            # p1 / (1 - p0 exp(-i omega dt)) to the scaled derivative.
+            # Forward electric derivatives use the H-update CPML parameters.
+            for family, parameter in enumerate(('param_h', 'param_e')):
+                p0, p1, p2 = face[parameter]
+                response = p2 * (1 + p1 / (1 - p0 * numpy.exp(-1j * OMEGA * DT)))
+                widths[family][axis][face['region'][axis]] /= response.ravel()
+    return widths
+
+
 def _build_complex_pulse_waveform(total_steps: int) -> tuple[numpy.ndarray, complex]:
     source_phasor, _delay = gaussian_packet(wl=WAVELENGTH, dwl=PULSE_DWL, dt=DT, turn_on=1e-5)
     aa, cc, ss = source_phasor(numpy.arange(total_steps) + 0.5)
@@ -450,7 +465,10 @@ def _run_straight_waveguide_case(variant: str) -> WaveguideCalibrationResult:
 def _run_width_step_scattering_case() -> WaveguideScatteringResult:
     epsilon = _build_scattering_epsilon()
     base_dxes = _build_uniform_dxes(SCATTERING_SHAPE)
-    stretched_dxes = _build_stretched_dxes(base_dxes)
+    # A correctly signed backward mode resolves the weak reflected wave;
+    # incident leakage no longer hides PML differences or startup transients.
+    stretched_dxes = _build_cpml_frequency_dxes(base_dxes)
+    discrete_frequency = 2 * numpy.sin(OMEGA * DT / 2) / DT
 
     source_mode = waveguide_3d.solve_mode(
         0,
@@ -539,14 +557,14 @@ def _run_width_step_scattering_case() -> WaveguideScatteringResult:
     e_fdfd = unvec(
         fdfd.solvers.generic(
             J=vec(j_ph),
-            omega=OMEGA,
+            omega=discrete_frequency,
             dxes=stretched_dxes,
             epsilon=vec(epsilon),
             matrix_solver_opts={'atol': 1e-10, 'rtol': 1e-7},
         ),
         SCATTERING_SHAPE[1:],
     )
-    h_fdfd = functional.e2h(OMEGA, stretched_dxes)(e_fdfd)
+    h_fdfd = functional.e2h(discrete_frequency, stretched_dxes)(e_fdfd)
 
     reflected_td = complex(vec(e_ph) @ vec(reflected_overlap).conj())
     reflected_fd = complex(vec(e_fdfd) @ vec(reflected_overlap).conj())

meanas/test/test_waveguide_mode_helpers.py

@@ -21,8 +21,8 @@ def test_embedded_hybrid_mode_obeys_both_maxwell_curls(axis, polarity) -> None:
     omega = 0.7
     # Cyclic rotation preserves the orientation of the vector cross product.
     permutation = numpy.roll(numpy.arange(3), axis + 1)
-    canonical = numpy.ones((3, 4, 3, 9))
+    rng = numpy.random.default_rng(4)
-    canonical[:, 1:3, :, :] = 4
+    canonical = numpy.repeat((1 + rng.uniform(0, 3, (3, 4, 3)))[..., None], 9, axis=3)
     epsilon = canonical[permutation].transpose(0, *(permutation + 1))
     shape = epsilon.shape[1:]
     dxes = [[numpy.full(n, 0.4) for n in shape] for _ in range(2)]
@@ -36,6 +36,7 @@ def test_embedded_hybrid_mode_obeys_both_maxwell_curls(axis, polarity) -> None:
         axis=axis, polarity=polarity, slices=slices, epsilon=epsilon,
     )
     e, h = [waveguide_3d.expand_e(mode[field], mode['wavenumber'], dxes, axis, polarity, slices) for field in ('E', 'H')]
+    assert norm(e[axis]) > 0.01 * norm(e)
     assert norm(h[axis]) > 0.01 * norm(h)
     region = [slice(None)] * 4
     region[axis + 1] = slice(1, -1)