[tests] test more 2D waveguide results

meanas/test/test_waveguide_2d_numerics.py

@@ -1,5 +1,6 @@
 import numpy
 from numpy.linalg import norm
+from numpy.testing import assert_allclose
 
 from ..fdmath import vec
 from ..fdfd import waveguide_2d
@@ -8,6 +9,10 @@ from ..fdfd import waveguide_2d
 OMEGA = 1 / 1500
 GRID_SHAPE = (5, 5)
 DXES_2D = [[numpy.ones(GRID_SHAPE[0]), numpy.ones(GRID_SHAPE[1])] for _ in range(2)]
+DXES_2D_NONUNIFORM = [[
+    numpy.array([1.0, 1.2, 0.9, 1.1, 1.3]),
+    numpy.array([0.8, 1.1, 1.0, 1.2, 0.9]),
+] for _ in range(2)]
 
 
 def build_asymmetric_epsilon() -> numpy.ndarray:
@@ -16,6 +21,10 @@ def build_asymmetric_epsilon() -> numpy.ndarray:
     return vec(epsilon)
 
 
+def build_mu_profile() -> numpy.ndarray:
+    return numpy.linspace(1.5, 2.2, 3 * GRID_SHAPE[0] * GRID_SHAPE[1])
+
+
 def test_waveguide_2d_solved_modes_are_ordered_and_low_residual() -> None:
     epsilon = build_asymmetric_epsilon()
     operator_e = waveguide_2d.operator_e(OMEGA, DXES_2D, epsilon)
@@ -109,3 +118,167 @@ def test_waveguide_2d_sensitivity_matches_finite_difference() -> None:
     assert sensitivity[target_index].real > 0
     assert finite_difference.real > 0
     assert 0.4 < ratio < 1.8
+
+
+def test_waveguide_2d_normalized_fields_h_are_finite_and_unit_normalized_with_mu() -> None:
+    epsilon = build_asymmetric_epsilon()
+    mu = build_mu_profile()
+
+    e_xy, wavenumber = waveguide_2d.solve_mode(
+        0,
+        omega=OMEGA,
+        dxes=DXES_2D,
+        epsilon=epsilon,
+        )
+    _e_ref, h_ref = waveguide_2d.normalized_fields_e(
+        e_xy,
+        wavenumber=wavenumber,
+        omega=OMEGA,
+        dxes=DXES_2D,
+        epsilon=epsilon,
+        mu=mu,
+        )
+    h_xy = numpy.concatenate(numpy.split(h_ref, 3)[:2])
+
+    e_field, h_field = waveguide_2d.normalized_fields_h(
+        h_xy,
+        wavenumber=wavenumber,
+        omega=OMEGA,
+        dxes=DXES_2D,
+        epsilon=epsilon,
+        mu=mu,
+        )
+    overlap = waveguide_2d.inner_product(e_field, h_field, DXES_2D, conj_h=True)
+
+    assert e_field.shape == (3 * GRID_SHAPE[0] * GRID_SHAPE[1],)
+    assert h_field.shape == (3 * GRID_SHAPE[0] * GRID_SHAPE[1],)
+    assert numpy.isfinite(e_field).all()
+    assert numpy.isfinite(h_field).all()
+    assert abs(overlap.real - 1.0) < 1e-10
+    assert abs(overlap.imag) < 1e-10
+
+
+def test_waveguide_2d_helper_operators_with_mu_return_finite_outputs() -> None:
+    epsilon = build_asymmetric_epsilon()
+    mu = build_mu_profile()
+
+    e_xy, wavenumber = waveguide_2d.solve_mode(
+        0,
+        omega=OMEGA,
+        dxes=DXES_2D,
+        epsilon=epsilon,
+        )
+    _e_ref, h_ref = waveguide_2d.normalized_fields_e(
+        e_xy,
+        wavenumber=wavenumber,
+        omega=OMEGA,
+        dxes=DXES_2D,
+        epsilon=epsilon,
+        mu=mu,
+        )
+    h_xy = numpy.concatenate(numpy.split(h_ref, 3)[:2])
+    n_pts = GRID_SHAPE[0] * GRID_SHAPE[1]
+
+    operators = [
+        ('exy2h', waveguide_2d.exy2h(wavenumber, OMEGA, DXES_2D, epsilon, mu), e_xy),
+        ('hxy2e', waveguide_2d.hxy2e(wavenumber, OMEGA, DXES_2D, epsilon, mu), h_xy),
+        ('hxy2h', waveguide_2d.hxy2h(wavenumber, DXES_2D, mu), h_xy),
+    ]
+
+    for _name, operator, vector in operators:
+        result = operator @ vector
+        assert operator.shape == (3 * n_pts, 2 * n_pts)
+        assert numpy.isfinite(operator.data).all()
+        assert result.shape == (3 * n_pts,)
+        assert numpy.isfinite(result).all()
+
+
+def test_waveguide_2d_error_helpers_with_mu_return_finite_values() -> None:
+    epsilon = build_asymmetric_epsilon()
+    mu = build_mu_profile()
+
+    e_xy, wavenumber = waveguide_2d.solve_mode(
+        0,
+        omega=OMEGA,
+        dxes=DXES_2D,
+        epsilon=epsilon,
+        )
+    e_field, h_field = waveguide_2d.normalized_fields_e(
+        e_xy,
+        wavenumber=wavenumber,
+        omega=OMEGA,
+        dxes=DXES_2D,
+        epsilon=epsilon,
+        mu=mu,
+        )
+
+    h_error = waveguide_2d.h_err(h_field, wavenumber, OMEGA, DXES_2D, epsilon, mu)
+    e_error = waveguide_2d.e_err(e_field, wavenumber, OMEGA, DXES_2D, epsilon, mu)
+
+    assert numpy.isfinite(h_error)
+    assert numpy.isfinite(e_error)
+    assert 0 < h_error < 0.1
+    assert 0 < e_error < 2.0
+
+
+def test_waveguide_2d_inner_product_phase_and_conjugation_branches() -> None:
+    epsilon = build_asymmetric_epsilon()
+    mu = build_mu_profile()
+
+    e_xy, wavenumber = waveguide_2d.solve_mode(
+        0,
+        omega=OMEGA,
+        dxes=DXES_2D,
+        epsilon=epsilon,
+        )
+    e_field, h_field = waveguide_2d.normalized_fields_e(
+        e_xy,
+        wavenumber=wavenumber,
+        omega=OMEGA,
+        dxes=DXES_2D,
+        epsilon=epsilon,
+        mu=mu,
+        )
+
+    overlap_no_conj = waveguide_2d.inner_product(e_field, h_field, DXES_2D, conj_h=False)
+    overlap_conj = waveguide_2d.inner_product(e_field, h_field, DXES_2D, conj_h=True)
+    overlap_phase = waveguide_2d.inner_product(e_field, h_field, DXES_2D, conj_h=True, prop_phase=0.3)
+
+    assert numpy.isfinite(overlap_no_conj)
+    assert numpy.isfinite(overlap_phase)
+    assert abs(overlap_no_conj.real - 1.0) < 1e-10
+    assert abs(overlap_no_conj.imag) < 1e-10
+    assert_allclose(overlap_phase / overlap_conj, numpy.exp(-0.15j), atol=1e-12, rtol=1e-12)
+
+
+def test_waveguide_2d_inner_product_trapezoid_branch_on_nonuniform_grid() -> None:
+    epsilon = build_asymmetric_epsilon()
+
+    e_xy, wavenumber = waveguide_2d.solve_mode(
+        0,
+        omega=OMEGA,
+        dxes=DXES_2D_NONUNIFORM,
+        epsilon=epsilon,
+        )
+    e_field, h_field = waveguide_2d.normalized_fields_e(
+        e_xy,
+        wavenumber=wavenumber,
+        omega=OMEGA,
+        dxes=DXES_2D_NONUNIFORM,
+        epsilon=epsilon,
+        )
+
+    overlap_rect = waveguide_2d.inner_product(e_field, h_field, DXES_2D_NONUNIFORM, conj_h=True)
+    overlap_trap = waveguide_2d.inner_product(
+        e_field,
+        h_field,
+        DXES_2D_NONUNIFORM,
+        conj_h=True,
+        trapezoid=True,
+        )
+
+    assert numpy.isfinite(overlap_rect)
+    assert numpy.isfinite(overlap_trap)
+    assert abs(overlap_rect.imag) < 1e-10
+    assert abs(overlap_trap.imag) < 1e-10
+    assert abs(overlap_rect - overlap_trap) > 0.1