diff --git a/meanas/test/test_eigensolvers_numerics.py b/meanas/test/test_eigensolvers_numerics.py new file mode 100644 index 0000000..8d0c5c9 --- /dev/null +++ b/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 diff --git a/meanas/test/test_waveguide_2d_numerics.py b/meanas/test/test_waveguide_2d_numerics.py new file mode 100644 index 0000000..5667edc --- /dev/null +++ b/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, + ) diff --git a/meanas/test/test_waveguide_mode_helpers.py b/meanas/test/test_waveguide_mode_helpers.py new file mode 100644 index 0000000..2bf77f2 --- /dev/null +++ b/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)