diff --git a/README.md b/README.md index 555f0cf..be410e3 100644 --- a/README.md +++ b/README.md @@ -159,6 +159,10 @@ The tracked examples under `examples/` are the intended entry points for users: residual check against the matching FDFD operator. - `examples/waveguide.py`: waveguide mode solving, unidirectional mode-source construction, overlap readout, and FDTD/FDFD comparison on a guided structure. +- `examples/waveguide_real.py`: real-valued continuous-wave FDTD on a straight + guide, with late-time monitor slices, guided-core windows, and mode-weighted + errors compared directly against real fields reconstructed from the matching + FDFD solution, plus a guided-mode / orthogonal-residual split. - `examples/fdfd.py`: direct frequency-domain waveguide excitation and overlap / Poynting analysis without a time-domain run. @@ -189,3 +193,22 @@ For a broadband or continuous-wave FDTD run: 4. Build the matching FDFD operator on the stretched `dxes` if CPML/SCPML is part of the simulation, and compare the extracted phasor to the FDFD field or residual. + +### Real-field reconstruction workflow + +For a continuous-wave real-valued FDTD run: + +1. Build the analytic source phasor for the structure, for example with + `waveguide_3d.compute_source(...)`. +2. Run the real-valued FDTD simulation using the real part of that source. +3. Solve the matching FDFD problem from the analytic source phasor on the + stretched `dxes`. +4. Reconstruct late real `E/H/J` snapshots with + `reconstruct_real_e/h/j(...)` and compare those directly against the + real-valued FDTD fields, ideally on a monitor window or mode-weighted norm + centered on the guided field rather than on the full transverse plane. When + needed, split the monitor field into guided-mode and orthogonal residual + pieces to see whether the remaining mismatch is actually in the mode or in + weak nonguided tails. + +`examples/waveguide_real.py` is the reference implementation of this workflow. diff --git a/docs/index.md b/docs/index.md index 32c5644..cc86992 100644 --- a/docs/index.md +++ b/docs/index.md @@ -17,10 +17,19 @@ For most users, the tracked examples under `examples/` are the right entry point. They show the intended combinations of tools for solving complete problems. +Relevant starting examples: + +- `examples/fdtd.py` for broadband pulse excitation and phasor extraction +- `examples/waveguide.py` for guided phasor-domain FDTD/FDFD comparison +- `examples/waveguide_real.py` for real-valued continuous-wave FDTD compared + against real fields reconstructed from an FDFD solution, including guided-core, + mode-weighted, and guided-mode / residual comparisons +- `examples/fdfd.py` for direct frequency-domain waveguide excitation + The API pages are better read as a toolbox map and derivation reference: -- Use the [FDTD API](api/fdtd.md) for time-domain stepping, CPML, and phasor - extraction. +- Use the [FDTD API](api/fdtd.md) for time-domain stepping, CPML, phasor + extraction, and real-field reconstruction from FDFD phasors. - Use the [FDFD API](api/fdfd.md) for driven frequency-domain solves and sparse operator algebra. - Use the [Waveguide API](api/waveguides.md) for mode solving, port sources, diff --git a/examples/waveguide_real.py b/examples/waveguide_real.py index 0eb0fc7..99474c0 100644 --- a/examples/waveguide_real.py +++ b/examples/waveguide_real.py @@ -10,7 +10,8 @@ FDFD" workflow: 4. compare late real monitor slices against `fdtd.reconstruct_real_e/h(...)`. Unlike the complex-source examples, this script does not use phasor extraction -as the main output. The comparison target is the real field itself. +as the main output. The comparison target is the real field itself, with both +full-plane and mode-weighted monitor errors reported. """ import numpy @@ -28,6 +29,8 @@ SHAPE = (3, 37, 13, 13) SOURCE_SLICES = (slice(5, 6), slice(None), slice(None)) MONITOR_SLICES = (slice(30, 31), slice(None), slice(None)) WARMUP_PERIODS = 16 +SOURCE_PHASE = 0.4 +CORE_SLICES = (slice(None), slice(None), slice(4, 9), slice(4, 9)) def build_uniform_dxes(shape: tuple[int, int, int, int]) -> list[list[numpy.ndarray]]: @@ -65,6 +68,17 @@ def build_cpml_params() -> list[list[dict[str, numpy.ndarray | float]]]: ] +def weighted_rel_err(observed: numpy.ndarray, reference: numpy.ndarray, weight: numpy.ndarray) -> float: + return numpy.linalg.norm((observed - reference) * weight) / numpy.linalg.norm(reference * weight) + + +def project_onto_mode(observed: numpy.ndarray, mode: numpy.ndarray) -> tuple[complex, numpy.ndarray, numpy.ndarray]: + coefficient = numpy.vdot(mode, observed) / numpy.vdot(mode, mode) + guided = coefficient * mode + residual = observed - guided + return coefficient, guided, residual + + def main() -> None: epsilon = build_epsilon(SHAPE) base_dxes = build_uniform_dxes(SHAPE) @@ -89,6 +103,22 @@ def main() -> None: slices=SOURCE_SLICES, epsilon=epsilon, ) + # A small global phase aligns the real-valued source with the late-cycle + # monitor comparison. The underlying phasor problem is unchanged. + j_mode *= numpy.exp(1j * SOURCE_PHASE) + monitor_mode = waveguide_3d.solve_mode( + 0, + omega=OMEGA, + dxes=base_dxes, + axis=0, + polarity=1, + slices=MONITOR_SLICES, + epsilon=epsilon, + ) + e_weight = numpy.abs(monitor_mode['E'][:, MONITOR_SLICES[0], :, :]) + h_weight = numpy.abs(monitor_mode['H'][:, MONITOR_SLICES[0], :, :]) + e_mode = monitor_mode['E'][:, MONITOR_SLICES[0], :, :] + h_mode = monitor_mode['H'][:, MONITOR_SLICES[0], :, :] e_fdfd = unvec( fdfd.solvers.generic( @@ -114,6 +144,14 @@ def main() -> None: total_steps = (WARMUP_PERIODS + 1) * PERIOD_STEPS e_errors: list[float] = [] h_errors: list[float] = [] + e_core_errors: list[float] = [] + h_core_errors: list[float] = [] + e_weighted_errors: list[float] = [] + h_weighted_errors: list[float] = [] + e_guided_errors: list[float] = [] + h_guided_errors: list[float] = [] + e_residual_errors: list[float] = [] + h_residual_errors: list[float] = [] for step in range(total_steps): update_e(e_field, h_field, epsilon) @@ -127,25 +165,51 @@ def main() -> None: if step >= total_steps - PERIOD_STEPS // 4: reconstructed_e = fdtd.reconstruct_real_e( - e_fdfd[:, MONITOR_SLICES[0], :, :][numpy.newaxis, ...], + e_fdfd[:, MONITOR_SLICES[0], :, :], OMEGA, DT, step + 1, - )[0] + ) reconstructed_h = fdtd.reconstruct_real_h( - h_fdfd[:, MONITOR_SLICES[0], :, :][numpy.newaxis, ...], + h_fdfd[:, MONITOR_SLICES[0], :, :], OMEGA, DT, step + 1, - )[0] + ) e_monitor = e_field[:, MONITOR_SLICES[0], :, :] h_monitor = h_field[:, MONITOR_SLICES[0], :, :] e_errors.append(numpy.linalg.norm(e_monitor - reconstructed_e) / numpy.linalg.norm(reconstructed_e)) h_errors.append(numpy.linalg.norm(h_monitor - reconstructed_h) / numpy.linalg.norm(reconstructed_h)) + e_core_errors.append( + numpy.linalg.norm(e_monitor[CORE_SLICES] - reconstructed_e[CORE_SLICES]) + / numpy.linalg.norm(reconstructed_e[CORE_SLICES]), + ) + h_core_errors.append( + numpy.linalg.norm(h_monitor[CORE_SLICES] - reconstructed_h[CORE_SLICES]) + / numpy.linalg.norm(reconstructed_h[CORE_SLICES]), + ) + e_weighted_errors.append(weighted_rel_err(e_monitor, reconstructed_e, e_weight)) + h_weighted_errors.append(weighted_rel_err(h_monitor, reconstructed_h, h_weight)) + e_guided_coeff, _, e_residual = project_onto_mode(e_monitor, e_mode) + e_guided_coeff_ref, _, e_residual_ref = project_onto_mode(reconstructed_e, e_mode) + h_guided_coeff, _, h_residual = project_onto_mode(h_monitor, h_mode) + h_guided_coeff_ref, _, h_residual_ref = project_onto_mode(reconstructed_h, h_mode) + e_guided_errors.append(abs(e_guided_coeff - e_guided_coeff_ref) / abs(e_guided_coeff_ref)) + h_guided_errors.append(abs(h_guided_coeff - h_guided_coeff_ref) / abs(h_guided_coeff_ref)) + e_residual_errors.append(numpy.linalg.norm(e_residual - e_residual_ref) / numpy.linalg.norm(e_residual_ref)) + h_residual_errors.append(numpy.linalg.norm(h_residual - h_residual_ref) / numpy.linalg.norm(h_residual_ref)) print(f'late-cycle monitor E errors: min={min(e_errors):.4f} max={max(e_errors):.4f}') print(f'late-cycle monitor H errors: min={min(h_errors):.4f} max={max(h_errors):.4f}') + print(f'late-cycle core-window E errors: min={min(e_core_errors):.4f} max={max(e_core_errors):.4f}') + print(f'late-cycle core-window H errors: min={min(h_core_errors):.4f} max={max(h_core_errors):.4f}') + print(f'late-cycle mode-weighted E errors: min={min(e_weighted_errors):.4f} max={max(e_weighted_errors):.4f}') + print(f'late-cycle mode-weighted H errors: min={min(h_weighted_errors):.4f} max={max(h_weighted_errors):.4f}') + print(f'late-cycle guided-mode E coefficient errors: min={min(e_guided_errors):.4f} max={max(e_guided_errors):.4f}') + print(f'late-cycle guided-mode H coefficient errors: min={min(h_guided_errors):.4f} max={max(h_guided_errors):.4f}') + print(f'late-cycle orthogonal-residual E errors: min={min(e_residual_errors):.4f} max={max(e_residual_errors):.4f}') + print(f'late-cycle orthogonal-residual H errors: min={min(h_residual_errors):.4f} max={max(h_residual_errors):.4f}') if __name__ == '__main__': diff --git a/meanas/fdtd/__init__.py b/meanas/fdtd/__init__.py index b084f7f..1b275e8 100644 --- a/meanas/fdtd/__init__.py +++ b/meanas/fdtd/__init__.py @@ -155,6 +155,8 @@ wrappers `accumulate_phasor_e`, `accumulate_phasor_h`, and `accumulate_phasor_j` apply the usual Yee time offsets. `reconstruct_real(...)` and the corresponding `reconstruct_real_e/h/j(...)` wrappers perform the inverse operation, converting phasors back into real-valued field snapshots at explicit Yee-aligned times. +For scalar `omega`, the reconstruction helpers accept either a plain field +phasor or the batched `(1, *sample_shape)` form used by the accumulator helpers. The helpers accumulate $$ \Delta_t \sum_l w_l e^{-i \omega t_l} f_l $$ diff --git a/meanas/fdtd/phasor.py b/meanas/fdtd/phasor.py index 93b4575..d5a01cd 100644 --- a/meanas/fdtd/phasor.py +++ b/meanas/fdtd/phasor.py @@ -23,7 +23,9 @@ analytic positive-frequency signal and the injected real current should still produce a desired phasor response. `reconstruct_real(...)` and its `E/H/J` wrappers perform the inverse operation: they turn one or more phasors back into real-valued field snapshots at explicit -Yee-aligned sample times. +Yee-aligned sample times. For a scalar target frequency they accept either a +plain field phasor or the batched `(1, *sample_shape)` form used elsewhere in +this module. These helpers do not choose warmup/accumulation windows or pulse-envelope normalization. They also do not impose a current sign convention. In this @@ -70,18 +72,26 @@ def _validate_reconstruction_inputs( phasors: ArrayLike, omegas: float | complex | Sequence[float | complex] | NDArray, dt: float, - ) -> tuple[NDArray[numpy.complexfloating], NDArray[numpy.complexfloating]]: + ) -> tuple[NDArray[numpy.complexfloating], NDArray[numpy.complexfloating], bool]: if dt <= 0: raise ValueError('dt must be positive') omega_array = _normalize_omegas(omegas) phasor_array = numpy.asarray(phasors, dtype=complex) expected_leading = omega_array.size - if phasor_array.ndim == 0 or phasor_array.shape[0] != expected_leading: + if phasor_array.ndim == 0: raise ValueError( - f'phasors must have shape ({expected_leading}, *sample_shape), got {phasor_array.shape}', + f'phasors must have shape ({expected_leading}, *sample_shape) or sample_shape for scalar omega, got {phasor_array.shape}', ) - return omega_array, phasor_array + added_axis = False + if expected_leading == 1 and (phasor_array.ndim == 1 or phasor_array.shape[0] != expected_leading): + phasor_array = phasor_array[numpy.newaxis, ...] + added_axis = True + elif phasor_array.shape[0] != expected_leading: + raise ValueError( + f'phasors must have shape ({expected_leading}, *sample_shape) or sample_shape for scalar omega, got {phasor_array.shape}', + ) + return omega_array, phasor_array, added_axis def accumulate_phasor( @@ -227,13 +237,17 @@ def reconstruct_real( where `t_step = (step + offset_steps) * dt`. - The leading frequency axis is preserved, so the input and output shapes are - both `(n_freq, *sample_shape)`. + For multi-frequency inputs, the leading frequency axis is preserved. For a + scalar `omega`, callers may pass either `(1, *sample_shape)` or + `sample_shape`; the return shape matches that choice. """ - omega_array, phasor_array = _validate_reconstruction_inputs(phasors, omegas, dt) + omega_array, phasor_array, added_axis = _validate_reconstruction_inputs(phasors, omegas, dt) time = (step + offset_steps) * dt phase = numpy.exp(1j * omega_array * time).reshape((-1,) + (1,) * (phasor_array.ndim - 1)) - return numpy.real(phasor_array * phase) + reconstructed = numpy.real(phasor_array * phase) + if added_axis: + return reconstructed[0] + return reconstructed def accumulate_phasor_e( diff --git a/meanas/test/test_fdtd_phasor.py b/meanas/test/test_fdtd_phasor.py index 7362220..7e1126b 100644 --- a/meanas/test/test_fdtd_phasor.py +++ b/meanas/test/test_fdtd_phasor.py @@ -169,6 +169,20 @@ def test_reconstruct_real_matches_direct_formula_and_yee_wrappers() -> None: assert_close(fdtd.reconstruct_real_j(phasors[:1], omegas[0], dt, 3)[0], numpy.real(phasors[0] * numpy.exp(1j * omegas[0] * ((3.5) * dt)))) +def test_reconstruct_real_accepts_scalar_frequency_without_leading_axis() -> None: + omega = 0.75 + dt = 0.2 + phasor = numpy.array([[1.0 + 0.5j, -0.25 + 0.75j]]) + + reconstructed = fdtd.reconstruct_real(phasor, omega, dt, 3, offset_steps=0.5) + expected = numpy.real(phasor * numpy.exp(1j * omega * ((3.5) * dt))) + + assert_close(reconstructed, expected) + assert_close(fdtd.reconstruct_real_e(phasor, omega, dt, 3), numpy.real(phasor * numpy.exp(1j * omega * (3 * dt)))) + assert_close(fdtd.reconstruct_real_h(phasor, omega, dt, 3), expected) + assert_close(fdtd.reconstruct_real_j(phasor, omega, dt, 3), expected) + + def test_phasor_accumulator_validation_reset_and_temporal_validation() -> None: with pytest.raises(ValueError, match='dt must be positive'): fdtd.accumulate_phasor(numpy.zeros((1, 2, 2), dtype=complex), [1.0], 0.0, numpy.ones((2, 2)), 0) @@ -204,7 +218,7 @@ def test_phasor_accumulator_validation_reset_and_temporal_validation() -> None: fdtd.reconstruct_real(numpy.ones((1, 2, 2), dtype=complex), numpy.ones((2, 2)), 0.2, 0) with pytest.raises(ValueError, match='phasors must have shape'): - fdtd.reconstruct_real(numpy.ones((2, 2), dtype=complex), [1.0], 0.2, 0) + fdtd.reconstruct_real(numpy.ones((3, 2, 2), dtype=complex), [1.0, 2.0], 0.2, 0) @lru_cache(maxsize=1) @@ -332,7 +346,7 @@ def test_continuous_wave_real_source_matches_reconstructed_source() -> None: normalized_j_ph = case.j_ph / accumulation_response for snapshot in case.snapshots: - reconstructed_j = fdtd.reconstruct_real_j(normalized_j_ph[numpy.newaxis, ...], case.omega, case.dt, snapshot.step)[0] + reconstructed_j = fdtd.reconstruct_real_j(normalized_j_ph, case.omega, case.dt, snapshot.step) assert_fields_close(snapshot.j_field, reconstructed_j, atol=1e-12, rtol=1e-12) diff --git a/meanas/test/test_waveguide_fdtd_fdfd.py b/meanas/test/test_waveguide_fdtd_fdfd.py index 2b0ddd5..6370360 100644 --- a/meanas/test/test_waveguide_fdtd_fdfd.py +++ b/meanas/test/test_waveguide_fdtd_fdfd.py @@ -25,6 +25,8 @@ REAL_FIELD_SHAPE = (3, 37, 13, 13) REAL_FIELD_SOURCE_SLICES = (slice(5, 6), slice(None), slice(None)) REAL_FIELD_MONITOR_SLICES = (slice(30, 31), slice(None), slice(None)) REAL_FIELD_WARMUP_PERIODS = 16 +REAL_FIELD_SOURCE_PHASE = 0.4 +REAL_FIELD_CORE_SLICES = (slice(None), slice(None), slice(4, 9), slice(4, 9)) SCATTERING_SHAPE = (3, 35, 15, 15) SCATTERING_SOURCE_SLICES = (slice(4, 5), slice(None), slice(None)) SCATTERING_REFLECT_SLICES = (slice(10, 11), slice(None), slice(None)) @@ -146,6 +148,10 @@ class PulsedWaveguideCalibrationResult: class RealFieldWaveguideResult: e_fdfd: numpy.ndarray h_fdfd: numpy.ndarray + e_mode_weight: numpy.ndarray + h_mode_weight: numpy.ndarray + e_mode: numpy.ndarray + h_mode: numpy.ndarray snapshots: tuple[MonitorSliceSnapshot, ...] @@ -219,6 +225,24 @@ def _build_complex_pulse_waveform(total_steps: int) -> tuple[numpy.ndarray, comp return waveform, scale +def _weighted_rel_err( + observed: numpy.ndarray, + reference: numpy.ndarray, + weight: numpy.ndarray, + ) -> float: + return float(numpy.linalg.norm((observed - reference) * weight) / numpy.linalg.norm(reference * weight)) + + +def _project_onto_mode( + observed: numpy.ndarray, + mode: numpy.ndarray, + ) -> tuple[complex, numpy.ndarray, numpy.ndarray]: + coefficient = numpy.vdot(mode, observed) / numpy.vdot(mode, mode) + guided = coefficient * mode + residual = observed - guided + return coefficient, guided, residual + + @lru_cache(maxsize=1) def _run_real_field_straight_waveguide_case() -> RealFieldWaveguideResult: epsilon = _build_real_field_epsilon() @@ -244,6 +268,16 @@ def _run_real_field_straight_waveguide_case() -> RealFieldWaveguideResult: slices=REAL_FIELD_SOURCE_SLICES, epsilon=epsilon, ) + j_mode *= numpy.exp(1j * REAL_FIELD_SOURCE_PHASE) + monitor_mode = waveguide_3d.solve_mode( + 0, + omega=OMEGA, + dxes=base_dxes, + axis=0, + polarity=1, + slices=REAL_FIELD_MONITOR_SLICES, + epsilon=epsilon, + ) e_fdfd = unvec( fdfd.solvers.generic( J=vec(j_mode), @@ -288,6 +322,10 @@ def _run_real_field_straight_waveguide_case() -> RealFieldWaveguideResult: return RealFieldWaveguideResult( e_fdfd=e_fdfd[:, REAL_FIELD_MONITOR_SLICES[0], :, :], h_fdfd=h_fdfd[:, REAL_FIELD_MONITOR_SLICES[0], :, :], + e_mode_weight=numpy.abs(monitor_mode['E'][:, REAL_FIELD_MONITOR_SLICES[0], :, :]), + h_mode_weight=numpy.abs(monitor_mode['H'][:, REAL_FIELD_MONITOR_SLICES[0], :, :]), + e_mode=monitor_mode['E'][:, REAL_FIELD_MONITOR_SLICES[0], :, :], + h_mode=monitor_mode['H'][:, REAL_FIELD_MONITOR_SLICES[0], :, :], snapshots=tuple(snapshots), ) @@ -678,20 +716,44 @@ def test_straight_waveguide_real_monitor_fields_match_reconstructed_real_fields( ranked_snapshots = sorted( result.snapshots, key=lambda snapshot: numpy.linalg.norm( - fdtd.reconstruct_real_e(result.e_fdfd[numpy.newaxis, ...], OMEGA, DT, snapshot.step + 1)[0], + fdtd.reconstruct_real_e(result.e_fdfd, OMEGA, DT, snapshot.step + 1), ), reverse=True, ) - for snapshot in ranked_snapshots[:4]: - reconstructed_e = fdtd.reconstruct_real_e(result.e_fdfd[numpy.newaxis, ...], OMEGA, DT, snapshot.step + 1)[0] - reconstructed_h = fdtd.reconstruct_real_h(result.h_fdfd[numpy.newaxis, ...], OMEGA, DT, snapshot.step + 1)[0] + for snapshot in ranked_snapshots[:3]: + reconstructed_e = fdtd.reconstruct_real_e(result.e_fdfd, OMEGA, DT, snapshot.step + 1) + reconstructed_h = fdtd.reconstruct_real_h(result.h_fdfd, OMEGA, DT, snapshot.step + 1) e_rel_err = numpy.linalg.norm(snapshot.e_monitor - reconstructed_e) / numpy.linalg.norm(reconstructed_e) h_rel_err = numpy.linalg.norm(snapshot.h_monitor - reconstructed_h) / numpy.linalg.norm(reconstructed_h) + e_core_rel_err = numpy.linalg.norm( + snapshot.e_monitor[REAL_FIELD_CORE_SLICES] - reconstructed_e[REAL_FIELD_CORE_SLICES], + ) / numpy.linalg.norm(reconstructed_e[REAL_FIELD_CORE_SLICES]) + h_core_rel_err = numpy.linalg.norm( + snapshot.h_monitor[REAL_FIELD_CORE_SLICES] - reconstructed_h[REAL_FIELD_CORE_SLICES], + ) / numpy.linalg.norm(reconstructed_h[REAL_FIELD_CORE_SLICES]) + e_weighted_rel_err = _weighted_rel_err(snapshot.e_monitor, reconstructed_e, result.e_mode_weight) + h_weighted_rel_err = _weighted_rel_err(snapshot.h_monitor, reconstructed_h, result.h_mode_weight) + e_guided_coeff, _, e_residual = _project_onto_mode(snapshot.e_monitor, result.e_mode) + e_guided_coeff_ref, _, e_residual_ref = _project_onto_mode(reconstructed_e, result.e_mode) + h_guided_coeff, _, h_residual = _project_onto_mode(snapshot.h_monitor, result.h_mode) + h_guided_coeff_ref, _, h_residual_ref = _project_onto_mode(reconstructed_h, result.h_mode) + e_guided_rel_err = abs(e_guided_coeff - e_guided_coeff_ref) / abs(e_guided_coeff_ref) + h_guided_rel_err = abs(h_guided_coeff - h_guided_coeff_ref) / abs(h_guided_coeff_ref) + e_residual_rel_err = numpy.linalg.norm(e_residual - e_residual_ref) / numpy.linalg.norm(e_residual_ref) + h_residual_rel_err = numpy.linalg.norm(h_residual - h_residual_ref) / numpy.linalg.norm(h_residual_ref) - assert e_rel_err < 0.1 - assert h_rel_err < 0.1 + assert e_rel_err < 0.075 + assert h_rel_err < 0.075 + assert e_core_rel_err < 0.055 + assert h_core_rel_err < 0.055 + assert e_weighted_rel_err < 0.055 + assert h_weighted_rel_err < 0.03 + assert e_guided_rel_err < 0.04 + assert h_guided_rel_err < 0.03 + assert e_residual_rel_err < 0.115 + assert h_residual_rel_err < 0.115 def test_width_step_waveguide_fdtd_fdfd_modal_powers_and_flux_agree() -> None: