From a08e9ce71b1e95b5eb4e8f1e40306ac3503db607 Mon Sep 17 00:00:00 2001 From: Forgejo Actions Date: Sun, 18 Oct 2026 23:21:29 -0700 Subject: [PATCH] [eme] characterize orientation bias under modal truncation Compare mixed and Galerkin testing against exact polarization/Fresnel interfaces and independent driven FDFD solves. Show complete-basis recovery, omitted-channel power, and the limits of longitudinal refinement, orientation averaging, and Euclidean least squares. Keep the default formulation unchanged and document its approximation limits. Validation: 1689 tests passed, 584 skipped; Ruff and documentation build passed. --- docs/api/waveguides.md | 62 ++++++++ meanas/fdfd/eme.py | 8 ++ meanas/test/test_eme_orientation.py | 210 ++++++++++++++++++++++++++++ 3 files changed, 280 insertions(+) create mode 100644 meanas/test/test_eme_orientation.py diff --git a/docs/api/waveguides.md b/docs/api/waveguides.md index 464baba..ba2d1de 100644 --- a/docs/api/waveguides.md +++ b/docs/api/waveguides.md @@ -89,3 +89,65 @@ devices. Singular interface and cascade systems also raise contextual errors rather than silently substituting pseudoinverses. ::: meanas.fdfd.eme + +## Orientation bias under modal truncation + +Mixed testing is a projected approximation: H continuity is tested on the left +and E continuity on the right. Swapping the ports changes the testing spaces +when the retained modal subspaces are incomplete. This effect survives checks +of reciprocity and unitarity. + +A reproducible Maxwell counterexample is a uniform isotropic medium with two +degenerate linear polarizations. Keep one polarization on each port, rotated +so their overlap is `a=0.8`. There is no material interface: true reflection is +zero and transmission into the retained receiver is `a`. The other polarization +carries the remaining transmitted power, `1-a**2=0.36`. + +| Calculation | R_left | R_right | T_left_to_right | T_right_to_left | +| --- | ---: | ---: | ---: | ---: | +| Exact full basis, then restrict ports | 0 | 0 | 0.8 | 0.8 | +| Truncated mixed | 0.219512 | -0.219512 | 0.975610 | 0.975610 | +| Swap ports, mixed, then reorder | -0.219512 | 0.219512 | 0.975610 | 0.975610 | +| Truncated left Galerkin | 0 | 0 | 1.25 | 0.8 | +| Euclidean E/H least squares | 0 | 0 | 0.8 | 0.8 | + +An independent driven 3D FDFD solve with PML agrees with the exact reference +from either incidence direction: transmission is 0.8 within `1e-12`, with +reflection below `2e-8` in the test fixture. Keeping both polarizations makes +mixed testing match the exact complete scattering matrix and restores +port-relabeling covariance. + +For the single retained polarization, mixed testing yields +`r=(1-a**2)/(1+a**2)` and `t=2*a/(1+a**2)`. The spurious reflection is quadratic +in a small angle between adjacent retained polarizations. However, longitudinal +refinement alone is insufficient: rotating a single retained vector through a +uniform degenerate medium makes the sampled chain approach artificial following +of that vector (`|T|` approaches 1), while the physical endpoint overlap remains +0.8. Keeping the complete degenerate block recovers the physical result at every +sampling density tested. + +Two tempting corrections are insufficient. Averaging the two port orientations +cancels reflection in this example but leaves transmission at 0.975610. Plain +Euclidean E/H least squares gives the correct uniform-medium result, but is not +a general replacement. A second fixture uses exact degenerate modes of diagonal +positive epsilon/mu tensors, with electric scaling `(1,2)` on the left and +`(1,3)` on the right and a 30-degree polarization rotation on both ports. After +single-mode truncation, least squares has unequal directional transmissions +(approximately 0.850839 and 1.092163) and a maximum singular value of 1.104947. +The complete Fresnel solution is reciprocal and lossless. + +The mixed default remains unchanged. Retain complete coupled/degenerate blocks +when the geometry can mix them, check modal-basis convergence as well as section +sampling, and compare against independent full-field results. Neither orientation +averaging nor gain/reciprocity postprocessing recovers omitted-channel physics. +Orientation discrepancy is a useful diagnostic, not a bound on transmission +error. These examples use real power-orthonormal propagating modes; for general +complex modes, the reciprocity constraint involves the modal orthogonality +matrix, as derived by [Svendsen, Haakestad, and Skaar](https://arxiv.org/abs/1301.2458). + +Reproduce the analytic, full-basis, candidate-formulation, taper-refinement, and +driven-FDFD checks with: + +```sh +OPENBLAS_NUM_THREADS=1 python -m pytest meanas/test/test_eme_orientation.py -q +``` diff --git a/meanas/fdfd/eme.py b/meanas/fdfd/eme.py index 87b6225..7c4980b 100644 --- a/meanas/fdfd/eme.py +++ b/meanas/fdfd/eme.py @@ -512,6 +512,12 @@ def get_s( use the primal fields. Wavenumbers label the modes and are not used in the interface continuity equations. Forward modes imply backward (E, -H). + Mixed testing assigns different test families to the two ports. Incomplete + modal subspaces can therefore give different approximations after a port + exchange, even when each result is reciprocal and unitary. Establish basis + convergence, especially within coupled or degenerate blocks; do not infer + physical accuracy from these matrix constraints alone. + ``force_nogain`` caps Euclidean singular values; ``force_reciprocal`` symmetrizes S. These optional numerical projections are not general power or reciprocity tests for open, evanescent, or non-power-normalized bases. @@ -734,6 +740,8 @@ def get_taper_s( Each interval propagates half its length in either adjacent basis. Accuracy is controlled by section sampling and basis/PML convergence, not ODE tolerances. Standard outgoing/incoming [left, right] ordering is used. + Longitudinal refinement alone cannot repair omitted coupled/degenerate + modes or remove the orientation bias of incomplete interface projections. """ prepared = prepare_taper(sections, dxes, formulation=formulation, measure=measure) return evaluate_taper_s(prepared, force_nogain=force_nogain, force_reciprocal=force_reciprocal) diff --git a/meanas/test/test_eme_orientation.py b/meanas/test/test_eme_orientation.py new file mode 100644 index 0000000..d3d7eaa --- /dev/null +++ b/meanas/test/test_eme_orientation.py @@ -0,0 +1,210 @@ +"""Characterize truncation bias against exact reciprocal Maxwell interfaces. + +These tests intentionally document mixed/Galerkin projection limitations; they +are not passivity-enforcement tests. Both complete modal bases and an independent +driven FDFD solve establish the reference for the truncated polarization case. +""" +import numpy as np +import pytest +from numpy.testing import assert_allclose +from scipy.sparse.linalg import spsolve + +from ..fdmath import vec, unvec +from ..fdfd import eme, scpml, solvers, waveguide_2d, waveguide_3d + + +DXES = [[np.ones(1), np.ones(1)] for _ in range(2)] + + +def _rotation(angle): + return np.array([[np.cos(angle), -np.sin(angle)], [np.sin(angle), np.cos(angle)]]) + + +def _polarizations(angle, scales=(1., 1.)): + """Two power-orthonormal signed-beta=1 modes on a one-cell transverse grid. + + E_xy = sqrt(2) diag(scales) Q, and J^-1 H_xy = sqrt(2) diag(1/scales) Q. + They solve Maxwell with diagonal epsilon=(1/sx^2,1/sy^2,1) and + mu=(sy^2,sx^2,1) at omega=1, including nonparallel E/H traces. + """ + rotation = _rotation(angle) + electric = np.diag(scales) @ rotation + magnetic = np.diag(1 / np.asarray(scales)) @ rotation + return [ + (np.sqrt(2) * np.array([electric[0, i], electric[1, i], 0], complex), + np.sqrt(2) * np.array([-magnetic[1, i], magnetic[0, i], 0], complex)) + for i in range(2) + ] + + +def _scattering(left, right, formulation='mixed'): + return eme.get_s(left, np.ones(len(left)), right, np.ones(len(right)), DXES, formulation=formulation) + + +def _reorder_ports(matrix): + n = matrix.shape[0] // 2 + order = [*range(n, 2*n), *range(n)] + return matrix[np.ix_(order, order)] + + +def _fresnel_reference(left_angle, right_angle, left_scales=(1., 1.), right_scales=(1., 1.)): + """Independent exact two-polarization interface in the supplied port bases.""" + yl, yr = 1 / np.asarray(left_scales)**2, 1 / np.asarray(right_scales)**2 + reflection = np.diag((yl - yr) / (yl + yr)) + transmission = np.diag(2 * np.sqrt(yl * yr) / (yl + yr)) + ql, qr = _rotation(left_angle), _rotation(right_angle) + return np.block([[ql.T @ reflection @ ql, ql.T @ transmission @ qr], + [qr.T @ transmission @ ql, -qr.T @ reflection @ qr]]) + + +def _trace_least_squares(left, right): + """Research comparator: unweighted Euclidean tangential E/H least squares. + + This is deliberately not a public formulation: correcting the isotropic + polarization example does not establish reciprocity or passivity generally. + """ + el = np.column_stack([e[:2] for e, _h in left]) + hl = np.column_stack([h[:2] for _e, h in left]) + er = np.column_stack([e[:2] for e, _h in right]) + hr = np.column_stack([h[:2] for _e, h in right]) + outgoing = np.block([[el, -er], [-hl, -hr]]) + incoming = np.block([[-el, er], [-hl, -hr]]) + return np.linalg.lstsq(outgoing, incoming, rcond=None)[0] + + +@pytest.mark.parametrize('overlap', [.5, .8, .95]) +def test_truncated_mixed_polarizations_show_orientation_bias_despite_unitarity(overlap): + angle = np.arccos(overlap) + left, right = _polarizations(0)[:1], _polarizations(angle)[:1] + result = _scattering(left, right) + reversed_result = _reorder_ports(_scattering(right, left)) + reflection = (1-overlap**2) / (1+overlap**2) + transmission = 2*overlap / (1+overlap**2) + assert_allclose(result, [[reflection, transmission], [transmission, -reflection]], atol=1e-14) + assert_allclose(reversed_result, [[-reflection, transmission], [transmission, reflection]], atol=1e-14) + assert_allclose(result.T @ result, np.eye(2), atol=1e-14) + assert_allclose(result, result.T, atol=1e-14) + assert_allclose(np.linalg.norm(result-reversed_result), 2*np.sqrt(2)*reflection, atol=1e-14) + + +@pytest.mark.parametrize(('left_scales', 'right_scales'), [((1., 1.), (1., 1.)), ((1., 2.), (1., 3.))]) +def test_complete_maxwell_bases_recover_fresnel_and_port_relabeling(left_scales, right_scales): + left_angle, right_angle = np.pi/6, np.arccos(.8) + left, right = _polarizations(left_angle, left_scales), _polarizations(right_angle, right_scales) + # Establish that the synthetic fields are actual degenerate Maxwell modes, + # not arbitrary profiles or mixed vectors with different eigenvalues. + for modes, scales in ((left, left_scales), (right, right_scales)): + sx, sy = scales + epsilon, mu = np.array([1/sx**2, 1/sy**2, 1]), np.array([sy**2, sx**2, 1]) + operator = waveguide_2d.operator_e(1, DXES, epsilon, mu) + for e, h in modes: + assert_allclose(operator @ e[:2], e[:2], atol=1e-14) + assert_allclose(waveguide_2d.exy2h(1, 1, DXES, epsilon, mu) @ e[:2], h, atol=1e-14) + gram = [[waveguide_2d.lorentz_overlap(a, b, DXES) for b in modes] for a in modes] + assert_allclose(gram, np.eye(2), atol=1e-14) + result = _scattering(left, right) + expected = _fresnel_reference(left_angle, right_angle, left_scales, right_scales) + assert_allclose(result, expected, atol=1e-14) + assert_allclose(_reorder_ports(_scattering(right, left)), expected, atol=1e-14) + assert_allclose(result.T @ result, np.eye(4), atol=1e-14) + assert_allclose(result, result.T, atol=1e-14) + + +def test_missing_polarization_accounts_for_retained_submatrix_power_deficit(): + overlap = .8 + full = _scattering(_polarizations(0), _polarizations(np.arccos(overlap))) + retained = full[np.ix_([0, 2], [0, 2])] + assert_allclose(retained, [[0, overlap], [overlap, 0]], atol=1e-14) + assert_allclose(abs(full[3, 0])**2, 1-overlap**2, atol=1e-14) + assert_allclose(np.sum(abs(full[:, 0])**2), 1, atol=1e-14) + assert_allclose(np.sum(abs(retained[:, 0])**2), overlap**2, atol=1e-14) + + +def test_posthoc_orientation_averaging_does_not_recover_physical_transmission(): + overlap = .8 + left, right = _polarizations(0)[:1], _polarizations(np.arccos(overlap))[:1] + averaged = .5 * (_scattering(left, right) + _reorder_ports(_scattering(right, left))) + assert_allclose(averaged[0, 0], 0, atol=1e-14) + assert_allclose(averaged[1, 1], 0, atol=1e-14) + assert abs(averaged[1, 0] - overlap) > .17 + + +def test_galerkin_testing_does_not_resolve_truncated_polarization_error(): + overlap = .8 + left, right = _polarizations(0)[:1], _polarizations(np.arccos(overlap))[:1] + result = _scattering(left, right, 'galerkin') + assert_allclose(result, [[0, overlap], [1/overlap, 0]], atol=1e-14) + assert not np.allclose(result, result.T) + + +def test_trace_least_squares_recovers_uniform_reference_but_is_not_a_general_fix(): + overlap = .8 + left, right = _polarizations(0)[:1], _polarizations(np.arccos(overlap))[:1] + assert_allclose(_trace_least_squares(left, right), [[0, overlap], [overlap, 0]], atol=1e-14) + angle = np.pi/6 + left, right = _polarizations(angle, (1., 2.)), _polarizations(angle, (1., 3.)) + reference = _fresnel_reference(angle, angle, (1., 2.), (1., 3.)) + assert_allclose(_trace_least_squares(left, right), reference, atol=1e-14) + truncated = _trace_least_squares(left[:1], right[:1]) + # All traces come from the exact reciprocal, lossless Maxwell bases above. + # Euclidean residual minimization nevertheless violates both constraints + # after truncation. Port-relabeling covariance alone is not sufficient. + assert np.linalg.norm(truncated-truncated.T) > .3 + assert np.linalg.svd(truncated, compute_uv=False).max() > 1.10 + assert_allclose(_reorder_ports(_trace_least_squares(right[:1], left[:1])), truncated, atol=1e-14) + + +@pytest.mark.parametrize('side', [0, 1]) +def test_driven_fdfd_confirms_uniform_polarization_reference(side): + omega, thickness = np.pi/8, 24 + size = 2*thickness + 64 + physical = [[np.ones(1), np.ones(1), np.ones(size)] for _ in range(2)] + epsilon = np.ones((3, 1, 1, size)) + solve_widths = physical + for polarity in (-1, 1): + solve_widths = scpml.stretch_with_scpml(solve_widths, 2, polarity, omega, thickness=thickness) + overlap = .8 + profiles = [np.array([1., 0.]), np.array([overlap, np.sqrt(1-overlap**2)])] + polarity = 1 if side == 0 else -1 + source_index = thickness+4 if side == 0 else thickness+52 + slices = (slice(None), slice(None), slice(source_index, source_index+1)) + mode = waveguide_3d.mode_from_2d(profiles[side], omega, omega=omega, dxes=physical, + axis=2, polarity=polarity, slices=slices, epsilon=epsilon) + current = waveguide_3d.compute_source(mode['E'], mode['wavenumber'], omega, physical, 2, polarity, slices, epsilon) + field = unvec(solvers.generic(omega, solve_widths, vec(current), vec(epsilon), matrix_solver=spsolve), epsilon.shape[1:]) + left_index, right_index = thickness+12, thickness+40 + beta = 2*np.arcsin(omega/2) + phase = np.exp(-1j*beta) + + def waves(index): + e0, e1 = field[:2, 0, 0, index], field[:2, 0, 0, index+1] + forward = (e1 - e0/phase) / (phase - 1/phase) + return forward, e0-forward + + fl, bl = waves(left_index) + fr, br = waves(right_index) + incoming = np.dot(profiles[side], fl if side == 0 else br) + reflected = np.dot(profiles[side], bl if side == 0 else fr) / incoming + transmitted = np.dot(profiles[1-side], fr if side == 0 else bl) / incoming + transmitted *= np.exp(1j*beta*(right_index-left_index)) + assert abs(reflected) < 5e-8 + assert_allclose(transmitted, overlap, atol=1e-12) + + +def test_longitudinal_refinement_cannot_repair_an_omitted_degenerate_polarization(): + angle = np.arccos(.8) + reflection_magnitudes = [] + for count in (3, 9, 33): + complete = [eme.ModeSection(z, _polarizations(angle*z), [1., 1.]) for z in np.linspace(0, 1, count)] + truncated = [eme.ModeSection(section.z, section.modes[:1], [1.]) for section in complete] + full = eme.get_taper_s(complete, DXES) + expected = _fresnel_reference(0, angle).astype(complex) * np.exp(-1j) + assert_allclose(full, expected, atol=1e-13) + small = eme.get_taper_s(truncated, DXES) + reflection_magnitudes.append(abs(small[0, 0])) + # There is no material taper: the basis is merely rotated within a + # degenerate subspace. Restricting the exact result gives |T|=.8, while + # the truncated chain approaches artificial adiabatic following, |T|=1. + assert abs(small[1, 0]) > .99 + assert reflection_magnitudes[-1] < reflection_magnitudes[0] / 10 + assert_allclose(abs(expected[2, 0]), .8, atol=1e-14)