[waveguide_2d] improve test coverage

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 < h_error < 1e-6
-    assert 0 < e_error < 2.0
+    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,

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)
+    epsilon = numpy.ones((3, *shape), dtype=float)
-    dxes = [[numpy.ones(5), numpy.ones(5), numpy.ones(1)] for _ in range(2)]
+    if dxes is None:
-    slices = (slice(slice_start, slice_start + 1), slice(None), slice(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:
+def relative_residual(lhs: numpy.ndarray, rhs: numpy.ndarray) -> float:
-    epsilon, dxes, slices, result = build_waveguide_3d_mode(slice_start=0, polarity=1)
+    return float(norm(lhs - rhs) / (norm(lhs) + norm(rhs)))
-    axis = 0
+
+
+@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