[tests] add some more tests around numerical self-consistency

meanas/test/test_eigensolvers_numerics.py

@@ -0,0 +1,47 @@
+import numpy
+from numpy.linalg import norm
+from scipy import sparse
+import scipy.sparse.linalg as spalg
+
+from ..eigensolvers import rayleigh_quotient_iteration, signed_eigensolve
+
+
+def test_rayleigh_quotient_iteration_with_linear_operator() -> None:
+    operator = sparse.diags([5.0, 3.0, 1.0, -2.0]).tocsr()
+    linear_operator = spalg.LinearOperator(
+        shape=operator.shape,
+        dtype=complex,
+        matvec=lambda vv: operator @ vv,
+        )
+
+    def dense_solver(
+            shifted_operator: spalg.LinearOperator,
+            rhs: numpy.ndarray,
+            ) -> numpy.ndarray:
+        basis = numpy.eye(operator.shape[0], dtype=complex)
+        columns = [shifted_operator.matvec(basis[:, ii]) for ii in range(operator.shape[0])]
+        dense_matrix = numpy.column_stack(columns)
+        return numpy.linalg.lstsq(dense_matrix, rhs, rcond=None)[0]
+
+    guess = numpy.array([0.0, 1.0, 1e-6, 0.0], dtype=complex)
+    eigval, eigvec = rayleigh_quotient_iteration(
+        linear_operator,
+        guess,
+        iterations=8,
+        solver=dense_solver,
+        )
+
+    residual = norm(operator @ eigvec - eigval * eigvec)
+    assert abs(eigval - 3.0) < 1e-12
+    assert residual < 1e-12
+
+
+def test_signed_eigensolve_negative_returns_largest_negative_mode() -> None:
+    operator = sparse.diags([5.0, 3.0, 1.0, -2.0]).tocsr()
+
+    eigvals, eigvecs = signed_eigensolve(operator, how_many=1, negative=True)
+
+    assert eigvals.shape == (1,)
+    assert eigvecs.shape == (4, 1)
+    assert abs(eigvals[0] + 2.0) < 1e-12
+    assert abs(eigvecs[3, 0]) > 0.99

meanas/test/test_waveguide_2d_numerics.py

@@ -0,0 +1,108 @@
+import numpy
+from numpy.linalg import norm
+
+from ..fdmath import vec
+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)]
+
+
+def build_asymmetric_epsilon() -> numpy.ndarray:
+    epsilon = numpy.ones((3, *GRID_SHAPE), dtype=float)
+    epsilon[:, 2, 1] = 2.0
+    return vec(epsilon)
+
+
+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)
+
+    e_xys, wavenumbers = waveguide_2d.solve_modes(
+        [0, 1],
+        omega=OMEGA,
+        dxes=DXES_2D,
+        epsilon=epsilon,
+        )
+
+    assert numpy.all(numpy.diff(numpy.real(wavenumbers)) <= 0)
+
+    for e_xy, wavenumber in zip(e_xys, wavenumbers, strict=True):
+        residual = norm(operator_e @ e_xy - (wavenumber ** 2) * e_xy) / norm(e_xy)
+        assert residual < 1e-6
+
+
+def test_waveguide_2d_normalized_fields_are_consistent() -> None:
+    epsilon = build_asymmetric_epsilon()
+    operator_h = waveguide_2d.operator_h(OMEGA, DXES_2D, epsilon)
+
+    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,
+        )
+    h_xy = numpy.concatenate(numpy.split(h_field, 3)[:2])
+
+    overlap = waveguide_2d.inner_product(e_field, h_field, DXES_2D, conj_h=True)
+    h_residual = norm(operator_h @ h_xy - (wavenumber ** 2) * h_xy) / norm(h_xy)
+
+    assert abs(overlap.real - 1.0) < 1e-10
+    assert abs(overlap.imag) < 1e-10
+    assert waveguide_2d.e_err(e_field, wavenumber, OMEGA, DXES_2D, epsilon) < 1e-6
+    assert waveguide_2d.h_err(h_field, wavenumber, OMEGA, DXES_2D, epsilon) < 1e-6
+    assert h_residual < 1e-6
+
+
+def test_waveguide_2d_sensitivity_matches_finite_difference() -> None:
+    epsilon = build_asymmetric_epsilon()
+    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,
+        )
+    sensitivity = waveguide_2d.sensitivity(
+        e_field,
+        h_field,
+        wavenumber=wavenumber,
+        omega=OMEGA,
+        dxes=DXES_2D,
+        epsilon=epsilon,
+        )
+
+    target_index = int(numpy.argmax(numpy.abs(sensitivity)))
+    delta = 1e-4
+    epsilon_perturbed = epsilon.copy()
+    epsilon_perturbed[target_index] += delta
+
+    _, perturbed_wavenumber = waveguide_2d.solve_mode(
+        0,
+        omega=OMEGA,
+        dxes=DXES_2D,
+        epsilon=epsilon_perturbed,
+        )
+    finite_difference = (perturbed_wavenumber - wavenumber) / delta
+
+    numpy.testing.assert_allclose(
+        sensitivity[target_index],
+        finite_difference,
+        rtol=0.1,
+        atol=1e-6,
+        )

meanas/test/test_waveguide_mode_helpers.py

@@ -0,0 +1,103 @@
+import numpy
+from numpy.linalg import norm
+
+from ..fdmath import vec
+from ..fdfd import waveguide_3d, waveguide_cyl
+
+
+OMEGA = 1 / 1500
+
+
+def test_waveguide_3d_solve_mode_and_expand_e_are_phase_consistent() -> None:
+    epsilon = numpy.ones((3, 5, 5, 1), dtype=float)
+    dxes = [[numpy.ones(5), numpy.ones(5), numpy.ones(1)] for _ in range(2)]
+    axis = 0
+    polarity = 1
+    slices = (slice(0, 1), slice(None), slice(None))
+
+    result = waveguide_3d.solve_mode(
+        0,
+        omega=OMEGA,
+        dxes=dxes,
+        axis=axis,
+        polarity=polarity,
+        slices=slices,
+        epsilon=epsilon,
+        )
+    expanded = waveguide_3d.expand_e(
+        E=result['E'],
+        wavenumber=result['wavenumber'],
+        dxes=dxes,
+        axis=axis,
+        polarity=polarity,
+        slices=slices,
+        )
+
+    dx_prop = 0.5 * numpy.array([dx[2][slices[2]] for dx in dxes]).sum()
+    expected_wavenumber = 2 / dx_prop * numpy.arcsin(result['wavenumber_2d'] * dx_prop / 2)
+    solved_slice = (slice(None), *slices)
+
+    assert result['E'].shape == epsilon.shape
+    assert result['H'].shape == epsilon.shape
+    assert numpy.isfinite(result['E']).all()
+    assert numpy.isfinite(result['H']).all()
+    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(
+        numpy.abs(result['E']).argmax(),
+        result['E'].shape,
+        )
+    values = expanded[component, :, y_index, z_index]
+    ratios = values[1:] / values[:-1]
+    expected_ratio = numpy.exp(-1j * result['wavenumber'])
+
+    numpy.testing.assert_allclose(ratios, expected_ratio, rtol=1e-6, atol=1e-9)
+
+
+def test_waveguide_cyl_solved_modes_are_ordered_and_low_residual() -> None:
+    shape = (5, 5)
+    dxes = [[numpy.ones(shape[0]), numpy.ones(shape[1])] for _ in range(2)]
+    epsilon = vec(numpy.ones((3, *shape), dtype=float))
+    rmin = 10.0
+
+    e_xys, angular_wavenumbers = waveguide_cyl.solve_modes(
+        [0, 1],
+        omega=OMEGA,
+        dxes=dxes,
+        epsilon=epsilon,
+        rmin=rmin,
+        )
+    operator = waveguide_cyl.cylindrical_operator(OMEGA, dxes, epsilon, rmin=rmin)
+
+    assert numpy.all(numpy.diff(numpy.real(angular_wavenumbers)) <= 0)
+
+    for e_xy, angular_wavenumber in zip(e_xys, angular_wavenumbers, strict=True):
+        eigenvalue = (angular_wavenumber / rmin) ** 2
+        residual = norm(operator @ e_xy - eigenvalue * e_xy) / norm(e_xy)
+        assert residual < 1e-6
+
+
+def test_waveguide_cyl_linear_wavenumbers_are_finite_and_ordered() -> None:
+    shape = (5, 5)
+    dxes = [[numpy.ones(shape[0]), numpy.ones(shape[1])] for _ in range(2)]
+    epsilon = vec(numpy.ones((3, *shape), dtype=float))
+
+    e_xys, angular_wavenumbers = waveguide_cyl.solve_modes(
+        [0, 1],
+        omega=OMEGA,
+        dxes=dxes,
+        epsilon=epsilon,
+        rmin=10.0,
+        )
+    linear_wavenumbers = waveguide_cyl.linear_wavenumbers(
+        e_xys,
+        angular_wavenumbers,
+        epsilon=epsilon,
+        dxes=dxes,
+        rmin=10.0,
+        )
+
+    assert numpy.isfinite(linear_wavenumbers).all()
+    assert numpy.all(numpy.real(linear_wavenumbers) > 0)
+    assert numpy.all(numpy.diff(numpy.real(linear_wavenumbers)) <= 0)