From 4f75d6384282d439017a3bddc79359301ca0ba65 Mon Sep 17 00:00:00 2001 From: Forgejo Actions Date: Wed, 28 Oct 2026 18:30:42 -0700 Subject: [PATCH] [tests] verify hybrid mode signs against matched CPML scattering references --- meanas/test/test_waveguide_fdtd_fdfd.py | 26 ++++++++++++++++++---- meanas/test/test_waveguide_mode_helpers.py | 5 +++-- 2 files changed, 25 insertions(+), 6 deletions(-) diff --git a/meanas/test/test_waveguide_fdtd_fdfd.py b/meanas/test/test_waveguide_fdtd_fdfd.py index 0488197..c1d52b8 100644 --- a/meanas/test/test_waveguide_fdtd_fdfd.py +++ b/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()) diff --git a/meanas/test/test_waveguide_mode_helpers.py b/meanas/test/test_waveguide_mode_helpers.py index 0531a59..e5a45a6 100644 --- a/meanas/test/test_waveguide_mode_helpers.py +++ b/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)) - canonical[:, 1:3, :, :] = 4 + rng = numpy.random.default_rng(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)