From f3500e497b7cd5b63a2d126c2c131cab06aa1774 Mon Sep 17 00:00:00 2001 From: Forgejo Actions Date: Thu, 25 Jun 2026 18:57:17 -0700 Subject: [PATCH] [waveguide_2d] improve test coverage --- meanas/test/test_waveguide_2d_numerics.py | 8 +- meanas/test/test_waveguide_mode_helpers.py | 127 +++++++++++++++++++-- 2 files changed, 121 insertions(+), 14 deletions(-) diff --git a/meanas/test/test_waveguide_2d_numerics.py b/meanas/test/test_waveguide_2d_numerics.py index 0005584..16ac613 100644 --- a/meanas/test/test_waveguide_2d_numerics.py +++ b/meanas/test/test_waveguide_2d_numerics.py @@ -129,6 +129,7 @@ def test_waveguide_2d_normalized_fields_h_are_finite_and_unit_normalized_with_mu omega=OMEGA, dxes=DXES_2D, epsilon=epsilon, + mu=mu, ) _e_ref, h_ref = waveguide_2d.normalized_fields_e( e_xy, @@ -167,6 +168,7 @@ def test_waveguide_2d_helper_operators_with_mu_return_finite_outputs() -> None: omega=OMEGA, dxes=DXES_2D, epsilon=epsilon, + mu=mu, ) _e_ref, h_ref = waveguide_2d.normalized_fields_e( e_xy, @@ -202,6 +204,7 @@ def test_waveguide_2d_error_helpers_with_mu_return_finite_values() -> None: omega=OMEGA, dxes=DXES_2D, epsilon=epsilon, + mu=mu, ) e_field, h_field = waveguide_2d.normalized_fields_e( e_xy, @@ -217,8 +220,8 @@ def test_waveguide_2d_error_helpers_with_mu_return_finite_values() -> None: assert numpy.isfinite(h_error) assert numpy.isfinite(e_error) - assert 0 < h_error < 0.1 - assert 0 < e_error < 2.0 + assert 0 < h_error < 1e-6 + assert 0 < e_error < 1e-6 def test_waveguide_2d_inner_product_phase_and_conjugation_branches() -> None: @@ -230,6 +233,7 @@ def test_waveguide_2d_inner_product_phase_and_conjugation_branches() -> None: omega=OMEGA, dxes=DXES_2D, epsilon=epsilon, + mu=mu, ) e_field, h_field = waveguide_2d.normalized_fields_e( e_xy, diff --git a/meanas/test/test_waveguide_mode_helpers.py b/meanas/test/test_waveguide_mode_helpers.py index dc51de0..c375073 100644 --- a/meanas/test/test_waveguide_mode_helpers.py +++ b/meanas/test/test_waveguide_mode_helpers.py @@ -6,6 +6,7 @@ import pytest import warnings from ..fdmath import vec, unvec +from ..fdmath.operators import deriv_forward, deriv_back from ..fdfd import waveguide_2d, waveguide_3d, waveguide_cyl @@ -16,15 +17,21 @@ def build_waveguide_3d_mode( *, slice_start: int, polarity: int, + axis: int = 0, + shape: tuple[int, int, int] = (5, 5, 1), + dxes: list[list[numpy.ndarray]] | None = None, ) -> tuple[numpy.ndarray, list[list[numpy.ndarray]], tuple[slice, slice, slice], waveguide_3d.Waveguide3DMode]: - epsilon = numpy.ones((3, 5, 5, 1), dtype=float) - dxes = [[numpy.ones(5), numpy.ones(5), numpy.ones(1)] for _ in range(2)] - slices = (slice(slice_start, slice_start + 1), slice(None), slice(None)) + epsilon = numpy.ones((3, *shape), dtype=float) + if dxes is None: + dxes = [[numpy.ones(size) for size in shape] for _ in range(2)] + slices_l = [slice(None), slice(None), slice(None)] + slices_l[axis] = slice(slice_start, slice_start + 1) + slices = tuple(slices_l) result = waveguide_3d.solve_mode( 0, omega=OMEGA, dxes=dxes, - axis=0, + axis=axis, polarity=polarity, slices=slices, epsilon=epsilon, @@ -55,10 +62,37 @@ def build_waveguide_cyl_mu_profile() -> numpy.ndarray: return numpy.linspace(1.5, 2.2, 3 * 5 * 5) -def test_waveguide_3d_solve_mode_and_expand_e_are_phase_consistent() -> None: - epsilon, dxes, slices, result = build_waveguide_3d_mode(slice_start=0, polarity=1) - axis = 0 +def relative_residual(lhs: numpy.ndarray, rhs: numpy.ndarray) -> float: + return float(norm(lhs - rhs) / (norm(lhs) + norm(rhs))) + + +@pytest.mark.parametrize( + ('axis', 'shape'), + [ + (0, (5, 5, 1)), + (1, (5, 5, 1)), + (2, (5, 5, 5)), + ], + ) +def test_waveguide_3d_solve_mode_and_expand_e_are_phase_consistent( + axis: int, + shape: tuple[int, int, int], + ) -> None: polarity = 1 + dxes = [ + [ + numpy.linspace(0.9 + 0.1 * eh + 0.05 * aa, 1.2 + 0.1 * eh + 0.05 * aa, size) + for aa, size in enumerate(shape) + ] + for eh in range(2) + ] + epsilon, dxes, slices, result = build_waveguide_3d_mode( + slice_start=2, + polarity=polarity, + axis=axis, + shape=shape, + dxes=dxes, + ) expanded = waveguide_3d.expand_e( E=result['E'], wavenumber=result['wavenumber'], @@ -68,7 +102,7 @@ def test_waveguide_3d_solve_mode_and_expand_e_are_phase_consistent() -> None: slices=slices, ) - dx_prop = 0.5 * numpy.array([dx[2][slices[2]] for dx in dxes]).sum() + dx_prop = 0.5 * numpy.array([dx[axis][slices[axis]] for dx in dxes]).sum() expected_wavenumber = 2 / dx_prop * numpy.arcsin(result['wavenumber_2d'] * dx_prop / 2) solved_slice = (slice(None), *slices) @@ -79,15 +113,18 @@ def test_waveguide_3d_solve_mode_and_expand_e_are_phase_consistent() -> None: assert abs(result['wavenumber'] - expected_wavenumber) < 1e-12 assert numpy.allclose(expanded[solved_slice], result['E'][solved_slice]) - component, _x, y_index, z_index = numpy.unravel_index( + component, *indices = numpy.unravel_index( numpy.abs(result['E']).argmax(), result['E'].shape, ) - values = expanded[component, :, y_index, z_index] + line_selection: list[int | slice] = list(indices) + line_selection[axis] = slice(None) + values = expanded[component][tuple(line_selection)] ratios = values[1:] / values[:-1] - expected_ratio = numpy.exp(-1j * result['wavenumber']) + phase_positions = numpy.real(dxes[0][axis]).cumsum() + expected_ratios = numpy.exp(-1j * result['wavenumber'] * numpy.diff(phase_positions)) - numpy.testing.assert_allclose(ratios, expected_ratio, rtol=1e-6, atol=1e-9) + numpy.testing.assert_allclose(ratios, expected_ratios, rtol=1e-6, atol=1e-9) @pytest.mark.parametrize( @@ -117,6 +154,17 @@ def test_waveguide_3d_compute_overlap_e_uses_adjacent_window( assert nonzero[:, 1].min() == expected_range[0] assert nonzero[:, 1].max() == expected_range[1] + expanded = waveguide_3d.expand_e( + E=result['E'], + wavenumber=result['wavenumber'], + dxes=dxes, + axis=0, + polarity=polarity, + slices=slices, + ) + numpy.testing.assert_allclose(numpy.sum(overlap * expanded.conj()), 1.0, atol=1e-12) + numpy.testing.assert_allclose(vec(expanded) @ vec(overlap).conj(), 1.0, atol=1e-12) + @pytest.mark.parametrize( ('polarity', 'slice_start', 'expected_index'), @@ -301,6 +349,61 @@ def test_waveguide_cyl_reconstruction_matches_2d_limit(use_mu: bool) -> None: numpy.testing.assert_allclose(cyl_h, straight_h, rtol=1e-8, atol=1e-8) +@pytest.mark.parametrize('use_mu', [False, True]) +def test_waveguide_cyl_reconstruction_is_finite_radius_consistent(use_mu: bool) -> None: + dxes, epsilon, rmin = build_waveguide_cyl_fixture(nonuniform=True, asymmetric=True) + mu = build_waveguide_cyl_mu_profile() if use_mu else None + mu_full = numpy.ones_like(epsilon) if mu is None else mu + + e_xy, angular_wavenumber = waveguide_cyl.solve_mode( + 0, + omega=OMEGA, + dxes=dxes, + epsilon=epsilon, + rmin=rmin, + mu=mu, + ) + e_field = waveguide_cyl.exy2e( + angular_wavenumber=angular_wavenumber, + dxes=dxes, + rmin=rmin, + epsilon=epsilon, + ) @ e_xy + h_field = waveguide_cyl.exy2h( + angular_wavenumber=angular_wavenumber, + omega=OMEGA, + dxes=dxes, + rmin=rmin, + epsilon=epsilon, + mu=mu, + ) @ e_xy + + Dfx, Dfy = deriv_forward(dxes[0]) + Dbx, Dby = deriv_back(dxes[1]) + Ta, Tb = waveguide_cyl.dxes2T(dxes, rmin) + ta = Ta.diagonal() + tb = Tb.diagonal() + beta = angular_wavenumber / rmin + + er, ey, ez = numpy.split(e_field, 3) + hr, hy, hz = numpy.split(h_field, 3) + eps_r, eps_y, eps_z = numpy.split(epsilon, 3) + mu_r, mu_y, mu_z = numpy.split(mu_full, 3) + + ez_rhs = ta * ( + (Dbx @ (tb * eps_r * er)) + + (Dby @ (ta * eps_y * ey)) + ) / eps_z + hr_rhs = 1j * beta * ey / ta + Dfy @ ez + hy_rhs = -1j * beta * er / tb - (Dfx @ (ta * ez)) / tb + hz_rhs = -Dfy @ er + Dfx @ ey + + assert relative_residual(1j * beta * ez, ez_rhs) < 1e-7 + assert relative_residual(-1j * OMEGA * mu_r * hr, hr_rhs) < 1e-7 + assert relative_residual(-1j * OMEGA * mu_y * hy, hy_rhs) < 1e-7 + assert relative_residual(-1j * OMEGA * mu_z * hz, hz_rhs) < 1e-7 + + def test_waveguide_cyl_linear_wavenumbers_use_component_radii() -> None: dxes, epsilon, rmin = build_waveguide_cyl_fixture(nonuniform=True) nx = dxes[0][0].size