From 9ac24892d6e1a789361c1a6b1785fba7567541cf Mon Sep 17 00:00:00 2001 From: Jan Petykiewicz Date: Fri, 17 Apr 2026 22:30:42 -0700 Subject: [PATCH] [tests] test more 2D waveguide results --- meanas/test/test_waveguide_2d_numerics.py | 173 ++++++++++++++++++++++ 1 file changed, 173 insertions(+) diff --git a/meanas/test/test_waveguide_2d_numerics.py b/meanas/test/test_waveguide_2d_numerics.py index 0bc5bf2..0005584 100644 --- a/meanas/test/test_waveguide_2d_numerics.py +++ b/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