[waveguide] correct modal reconstruction and calibrate full port bases

meanas/fdfd/waveguide_2d.py

@@ -801,77 +801,40 @@ def sensitivity(
         wavenumber: complex,
         omega: complex,
         dxes: dx_lists2_t,
-        epsilon: vfdslice,
+        epsilon: vfdslice | vcfdslice,
-        mu: vfdslice | None = None,
+        mu: vfdslice | vcfdslice | None = None,
         ) -> vcfdslice_t:
     r"""
-      Given a waveguide structure (`dxes`, `epsilon`, `mu`) and mode fields
+    Calculate the analytic sensitivity of an isolated mode's beta to epsilon.
-    (`e_norm`, `h_norm`, `wavenumber`, `omega`), calculates the sensitivity of the wavenumber
-    $\beta$ to changes in the dielectric structure $\epsilon$.
 
-    The output is a vector of the same size as `vec(epsilon)`, with each element specifying the
+    For the reciprocal diagonal waveguide problem ``A_E E_xy = beta**2 E_xy``,
-    sensitivity of `wavenumber` to changes in the corresponding element in `vec(epsilon)`, i.e.
+    the Yee-weighted Lorentz left eigenvector is the unconjugated row
 
-    $$sens_{i} = \frac{\partial\beta}{\partial\epsilon_i}$$
+    $$ L^T = [H_y \Delta x_e \Delta y_h,\; -H_x \Delta x_h \Delta y_e]. $$
 
-    An adjoint approach is used to calculate the sensitivity; the derivation is provided here:
+    Differentiating the eigenproblem and projecting with this left vector gives
 
-    Starting with the eigenvalue equation
+    $$ \frac{\partial\beta}{\partial\epsilon_i}
-
+       = \frac{L^T (\partial_{\epsilon_i} A_E) E_{xy}}
-    $$\beta^2 E_{xy} = A_E E_{xy}$$
+              {2\beta L^T E_{xy}}. $$
-
-    where $A_E$ is the waveguide operator from `operator_e()`, and $E_{xy} = \begin{bmatrix} E_x \\
-                                                                                             E_y \end{bmatrix}$,
-    we can differentiate with respect to one of the $\epsilon$ elements (i.e. at one Yee grid point), $\epsilon_i$:
-
-    $$
-    (2 \beta) \partial_{\epsilon_i}(\beta) E_{xy} + \beta^2 \partial_{\epsilon_i} E_{xy}
-        = \partial_{\epsilon_i}(A_E) E_{xy} + A_E \partial_{\epsilon_i} E_{xy}
-    $$
-
-    We then multiply by $H_{yx}^\star = \begin{bmatrix}H_y^\star \\ -H_x^\star \end{bmatrix}$ from the left:
-
-    $$
-    (2 \beta) \partial_{\epsilon_i}(\beta) H_{yx}^\star E_{xy} + \beta^2 H_{yx}^\star \partial_{\epsilon_i} E_{xy}
-        = H_{yx}^\star \partial_{\epsilon_i}(A_E) E_{xy} + H_{yx}^\star A_E \partial_{\epsilon_i} E_{xy}
-    $$
-
-    On a graded Yee grid, the left projection includes the quadrature weights:
-    $H_y^\star \Delta x_e \Delta y_h$ and $-H_x^\star \Delta x_h \Delta y_e$.
-    These weights multiply the left vector, before it acts on any derivative operators.
-    With these weights understood, $H_{yx}^\star$ is a left-eigenvector of $A_E$. This can be verified by inspecting
-    the form of `operator_h` ($A_H$) and comparing its conjugate transpose to `operator_e` ($A_E$). Also, note
-    $H_{yx}^\star \cdot E_{xy} = H^\star \times E$ recalls the mode orthogonality relation. See doi:10.5194/ars-9-85-201
-    for a similar approach. Therefore,
-
-    $$
-    H_{yx}^\star A_E \partial_{\epsilon_i} E_{xy} = \beta^2 H_{yx}^\star \partial_{\epsilon_i} E_{xy}
-    $$
-
-    and we can simplify to
-
-    $$
-    \partial_{\epsilon_i}(\beta)
-        = \frac{1}{2 \beta} \frac{H_{yx}^\star \partial_{\epsilon_i}(A_E) E_{xy} }{H_{yx}^\star E_{xy}}
-    $$
-
-    This expression can be quickly calculated for all $i$ by writing out the various terms of
-    $\partial_{\epsilon_i} A_E$ and recognizing that the vector-matrix-vector products (i.e. scalars)
-    $sens_i = \vec{v}_{left} \partial_{\epsilon_i} (\epsilon_{xyz}) \vec{v}_{right}$, indexed by $i$, can be expressed as
-    elementwise multiplications $\vec{sens} = \vec{v}_{left} \star \vec{v}_{right}$
 
+    This pairing is bilinear: neither H nor the coordinate weights are
+    conjugated. It therefore applies to complex PML stretches and reciprocal
+    lossy materials as well as lossless grids. Quadrature weights multiply the
+    left projection before derivative operators act. Overall field scaling
+    cancels between numerator and denominator.
 
     Args:
-        e_norm: Normalized, vectorized E_xyz field for the mode. E.g. as returned by `normalized_fields_e`.
+        e_norm: Vectorized E_xyz eigenfield for the mode.
-        h_norm: Normalized, vectorized H_xyz field for the mode. E.g. as returned by `normalized_fields_e`.
+        h_norm: Matching vectorized H_xyz eigenfield at the same phase reference.
-        wavenumber: Propagation constant for the mode. The z-axis is assumed to be continuous (i.e. without numerical dispersion).
+        wavenumber: Nonzero beta of the continuous-z eigenproblem.
-        omega: The angular frequency of the system.
+        omega: Angular frequency of the system.
-        dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.fdmath.types` (2D)
+        dxes: Two-dimensional Yee widths, including any complex PML stretches.
-        epsilon: Vectorized dielectric constant grid
+        epsilon: Vectorized diagonal permittivity.
-        mu: Vectorized magnetic permeability grid (default 1 everywhere)
+        mu: Optional vectorized diagonal permeability; defaults to one.
 
     Returns:
-        Sparse matrix representation of the operator.
+        Vector ``d beta / d epsilon_i`` in the same component order as epsilon.
     """
     if mu is None:
         mu = numpy.ones_like(epsilon)
@@ -891,12 +854,12 @@ def sensitivity(
     ev_xy = numpy.concatenate(numpy.split(e_norm, 3)[:2])
     hx, hy, _hz = numpy.split(h_norm, 3)
     # The Yee-weighted H projection is the left eigenvector on a graded grid.
-    hv_yx_conj = numpy.conj(numpy.concatenate([hy, -hx])) * numpy.concatenate([da_exxhyy, da_eyyhxx])
+    hv_yx = numpy.concatenate([hy, -hx]) * numpy.concatenate([da_exxhyy, da_eyyhxx])
 
-    sens_xy1 = (hv_yx_conj @ (omega * omega * mu_yx)) * ev_xy
+    sens_xy1 = (hv_yx @ (omega * omega * mu_yx)) * ev_xy
-    sens_xy2 = (hv_yx_conj @ sparse.vstack((Dfx, Dfy)) @ eps_z_inv @ sparse.hstack((Dbx, Dby))) * ev_xy
+    sens_xy2 = (hv_yx @ sparse.vstack((Dfx, Dfy)) @ eps_z_inv @ sparse.hstack((Dbx, Dby))) * ev_xy
-    sens_z   = (hv_yx_conj @ sparse.vstack((Dfx, Dfy)) @ (-eps_z_inv * eps_z_inv)) * (sparse.hstack((Dbx, Dby)) @ eps_xy @ ev_xy)
+    sens_z   = (hv_yx @ sparse.vstack((Dfx, Dfy)) @ (-eps_z_inv * eps_z_inv)) * (sparse.hstack((Dbx, Dby)) @ eps_xy @ ev_xy)
-    norm = hv_yx_conj @ ev_xy
+    norm = hv_yx @ ev_xy
 
     sens_tot = numpy.concatenate([sens_xy1 + sens_xy2, sens_z]) / (2 * wavenumber * norm)
     return vcfdslice_t(sens_tot)

meanas/fdfd/waveguide_3d.py

@@ -143,6 +143,9 @@ def mode_from_2d(
 
     ``polarity`` is +1 toward increasing indices and -1 toward decreasing
     indices. ``normalization`` follows waveguide_2d.normalized_fields_e.
+    Lossless modes at or beyond the longitudinal grid band edge cannot be
+    power-normalized. Lorentz normalization can represent these modes
+    beyond that edge and selects the decaying branch for real beta.
     """
     if polarity not in (-1, 1):
         raise ValueError('polarity must be -1 or +1')
@@ -154,7 +157,12 @@ def mode_from_2d(
         raise ValueError('the 2D eigenvector and wavenumber must be finite and match the cross-section')
     wavenumber_2d = complex(wavenumber_2d)
     # Correct wavenumber to account for numerical dispersion.
-    wavenumber = 2 / dx_prop * numpy.arcsin(wavenumber_2d * dx_prop / 2)
+    argument = wavenumber_2d * dx_prop / 2
+    if numpy.isreal(argument) and abs(argument.real) >= 1:
+        if normalization == 'power':
+            raise ValueError('mode is at or beyond the longitudinal grid band edge; refine the propagation grid')
+        argument = complex(argument.real, -0.0)
+    wavenumber = 2 / dx_prop * numpy.arcsin(argument)
 
     shape = [d.size for d in args_2d['dxes'][0]]
     ve, vh = waveguide_2d.normalized_fields_e(e_xy, wavenumber=wavenumber_2d, prop_phase=dx_prop * wavenumber, normalization=normalization, **args_2d)
@@ -372,12 +380,18 @@ def expand_e(
         $$
 
         to each copied slice.
+
+        Transverse electric samples use forward intervals dx_e[i] between
+        indices i and i+1. The longitudinal electric component uses backward
+        intervals dx_h[i] between indices i-1 and i on its staggered grid.
     """
     slices = tuple(slices)
 
     # Determine phase factors for parallel slices
     a_shape = numpy.roll([1, -1, 1, 1], axis)
-    a_E = numpy.real(dxes[0][axis]).cumsum()
+    # dx_e[i] spans transverse electric samples i -> i+1. The longitudinal
+    # component lies on the complementary grid, with backward intervals dx_h.
+    a_E = numpy.concatenate(([0.0], numpy.real(dxes[0][axis][:-1]).cumsum()))
     r_E = a_E - a_E[slices[axis]]
     iphi = polarity * -1j * wavenumber
     phase_E = numpy.exp(iphi * r_E).reshape(a_shape)
@@ -392,6 +406,9 @@ def expand_e(
     slices_in = (slice(None), *slices)
 
     E_expanded[slices_exp] = phase_E * numpy.array(E)[slices_in]
+    a_long = numpy.real(dxes[1][axis]).cumsum()
+    phase_long = numpy.exp(iphi * (a_long - a_long[slices[axis]])).reshape(a_shape[1:])
+    E_expanded[(axis, *slices_exp_l)] = phase_long * E[(axis, *slices)]
     return cfdfield_t(E_expanded)
 
 
@@ -550,10 +567,41 @@ def _decompose_port_field(
     matrix, overlaps, mode_count = decomposition
     field_vec = vec(field)
     measured = numpy.asarray([field_vec @ vec(overlap).conj() for overlap in overlaps])
-    coeffs, _resid, _rank, _sing = numpy.linalg.lstsq(matrix, measured, rcond=None)
+    coeffs, _resid, rank, _sing = numpy.linalg.lstsq(matrix, measured, rcond=None)
+    if rank != matrix.shape[1]:
+        raise numpy.linalg.LinAlgError('monitor modal basis is rank deficient')
     return coeffs[:mode_count], coeffs[mode_count:]
 
 
+def _port_power_gram(
+        modes: Sequence[Waveguide3DMode],
+        dxes: dx_lists_t,
+        axis: int,
+        polarity: int,
+        slices: Sequence[slice],
+        ) -> NDArray[numpy.complex128]:
+    """Hermitian flux metric, including interference between nonorthogonal modes."""
+    u, v = (axis + 1) % 3, (axis + 2) % 3
+
+    def width(family: int, component: int) -> NDArray:
+        values = numpy.real(dxes[family][component][slices[component]])
+        shape = [1, 1, 1]
+        shape[component] = values.size
+        return values.reshape(shape)
+
+    area_uv = width(0, u) * width(1, v)
+    area_vu = width(1, u) * width(0, v)
+    dx = float(numpy.real(dxes[0][axis][slices[axis]]).item())
+    gram = numpy.empty((len(modes), len(modes)), dtype=complex)
+    for row, test in enumerate(modes):
+        h = test['H'][(slice(None), *slices)].conj()
+        for column, mode in enumerate(modes):
+            e = mode['E'][(slice(None), *slices)]
+            phase = numpy.exp(-1j * polarity * mode['wavenumber'] * dx)
+            gram[row, column] = 0.5 * polarity * phase * numpy.sum(e[u] * h[v] * area_uv - e[v] * h[u] * area_vu)
+    return 0.5 * (gram + gram.conj().T)
+
+
 def _copy_cross_section(
         target: NDArray,
         cross_section: NDArray,
@@ -637,8 +685,11 @@ def get_s_fdfd(
     Extract a two-port guided-mode scattering matrix from local 3D FDFD solves.
 
     The returned S-matrix is ordered as `[[R_L, T_R_to_L], [T_L_to_R, R_R]]`, with
-    reference planes at the supplied monitor slices. Each column is normalized
+    reference planes at the supplied monitor slices. Calibration solves
-    by the incident modal amplitude measured on the corresponding input monitor.
+    ``S @ incoming_amplitudes = raw_outgoing_amplitudes`` using every measured
+    incoming coefficient on both ports, including mixed or degenerate modes.
+    Diagnostic guided powers use the full input/output flux Gram matrices,
+    retaining interference between nonorthogonal modal basis vectors.
     """
     if incident_amplitude_atol <= 0:
         raise ValueError('incident_amplitude_atol must be positive')
@@ -710,8 +761,8 @@ def get_s_fdfd(
         slices=right_monitor,
     )
 
-    s = numpy.zeros((n_total, n_total), dtype=complex)
+    raw_outgoing = numpy.zeros((n_total, n_total), dtype=complex)
-    raw_outgoing = numpy.zeros_like(s)
+    incoming_amplitudes = numpy.zeros_like(raw_outgoing)
     incident_amplitudes = numpy.zeros((n_total,), dtype=complex)
     electric_fields: list[cfdfield_t] = []
     magnetic_fields: list[cfdfield_t] = []
@@ -747,26 +798,39 @@ def get_s_fdfd(
         e_field = unvec(solvers.generic(**solve_kwargs), shape)
         left_in, left_out = _decompose_port_field(left_decomposition, e_field)
         right_in, right_out = _decompose_port_field(right_decomposition, e_field)
-        if column < n_left:
+        incoming_column = numpy.concatenate((left_in, right_in))
-            incident = left_in[column]
+        if numpy.linalg.norm(incoming_column) <= incident_amplitude_atol:
-        else:
-            incident = right_in[column - n_left]
-        if abs(incident) <= incident_amplitude_atol:
             raise ValueError(f'incident modal amplitude for column {column} is zero')
 
         raw_column = numpy.concatenate((left_out, right_out))
 
-        incident_amplitudes[column] = incident
+        incident_amplitudes[column] = incoming_column[column]
+        incoming_amplitudes[:, column] = incoming_column
         raw_outgoing[:, column] = raw_column
-        s[:, column] = raw_column / incident
         if return_diagnostics:
             electric_fields.append(cfdfield_t(e_field))
             magnetic_fields.append(cfdfield_t(functional.e2h(omega, solve_dxes, mu=mu)(e_field)))
 
+    if not numpy.isfinite(incoming_amplitudes).all() or not numpy.isfinite(raw_outgoing).all():
+        raise ValueError('modal calibration requires finite incoming and outgoing amplitudes')
+    try:
+        s = numpy.linalg.solve(incoming_amplitudes.T, raw_outgoing.T).T
+    except numpy.linalg.LinAlgError as error:
+        raise numpy.linalg.LinAlgError('incoming excitation matrix is singular') from error
+
     if not return_diagnostics:
         return s
 
-    guided_powers = numpy.sum(numpy.abs(s) ** 2, axis=0).real
+    zero = numpy.zeros((n_left, n_right), dtype=complex)
+    input_power_gram = numpy.block([
+        [_port_power_gram(left_in_modes, dxes, axis, +1, left_monitor), zero],
+        [zero.T, _port_power_gram(right_in_modes, dxes, axis, -1, right_monitor)],
+    ])
+    output_power_gram = numpy.block([
+        [_port_power_gram(left_out_modes, dxes, axis, -1, left_monitor), zero],
+        [zero.T, _port_power_gram(right_out_modes, dxes, axis, +1, right_monitor)],
+    ])
+    guided_powers = numpy.diag(s.conj().T @ output_power_gram @ s).real / numpy.diag(input_power_gram).real
     diagnostics: dict[str, Any] = {
         'axis': axis,
         'left_mode_numbers': left_modes_i,
@@ -781,6 +845,9 @@ def get_s_fdfd(
         'electric_fields': tuple(electric_fields),
         'magnetic_fields': tuple(magnetic_fields),
         'incident_amplitudes': incident_amplitudes,
+        'incoming_amplitudes': incoming_amplitudes,
+        'incoming_power_gram': input_power_gram,
+        'outgoing_power_gram': output_power_gram,
         'raw_outgoing_amplitudes': raw_outgoing,
         'guided_powers': guided_powers,
         'unresolved_powers': 1.0 - guided_powers,

meanas/fdfd/waveguide_cyl.py

@@ -86,7 +86,7 @@ large-radius-compatible identities
 
 $$
 E_z =
-\frac{1}{\imath \beta} T_a \epsilon_{zz}^{-1}
+\frac{1}{\imath \beta} \epsilon_{zz}^{-1}
 \begin{bmatrix}
 \hat{\partial}_r T_b \epsilon_{rr} &
 \hat{\partial}_y T_a \epsilon_{yy}
@@ -304,6 +304,11 @@ def linear_wavenumbers(
     Calculate linear wavenumbers (1/distance) based on angular wavenumbers (1/rad)
       and the mode's energy distribution.
 
+    The effective radius is weighted by transverse electric energy, using the
+    physical Yee cell areas for each component. Radial E uses dx_e * dy_h;
+    transverse E uses dx_h * dy_e. Real coordinate widths define these areas
+    and radii when complex PML stretching is present.
+
     Args:
         e_xys: Vectorized mode fields with shape (num_modes, 2 * x *y)
         angular_wavenumbers: Wavenumbers assuming fields have theta-dependence of
@@ -321,14 +326,19 @@ def linear_wavenumbers(
 
     shape2d = (len(dxes[0][0]), len(dxes[0][1]))
     epsilon2d = unvec(epsilon, shape2d)[:2]
-    ra = rmin + numpy.cumsum(dxes[0][0])
+    physical_widths = [[numpy.real(dx) for dx in family] for family in dxes]
-    rb = rmin + dxes[0][0] / 2.0 + numpy.concatenate((
+    ra = rmin + numpy.cumsum(physical_widths[0][0])
-        numpy.zeros(1, dtype=dxes[1][0].dtype),
+    rb = rmin + physical_widths[0][0] / 2.0 + numpy.concatenate((
-        numpy.cumsum(dxes[1][0][:-1]),
+        numpy.zeros(1),
+        numpy.cumsum(physical_widths[1][0][:-1]),
         ))
+    areas = numpy.stack((
+        numpy.outer(physical_widths[0][0], physical_widths[1][1]),
+        numpy.outer(physical_widths[1][0], physical_widths[0][1]),
+    ))
     for ii in range(angular_wavenumbers.size):
         efield = unvec(e_xys[ii], shape2d, 2)
-        energy = numpy.real((efield * efield.conj()) * epsilon2d)
+        energy = numpy.real((efield * efield.conj()) * epsilon2d) * areas
         er_energy_vs_r = energy[0].sum(axis=1)
         ey_energy_vs_r = energy[1].sum(axis=1)
         energy_vs_r = er_energy_vs_r + ey_energy_vs_r
@@ -405,7 +415,7 @@ def exy2e(
 
     n_pts = dxes[0][0].size * dxes[0][1].size
     exy2ez = (
-        Ta @ epsilon_z_inv
+        epsilon_z_inv
         @ sparse.hstack((Dbx @ Tb @ epsilon_x,
                          Dby @ Ta @ epsilon_y))
         / (1j * wavenumber)

meanas/test/test_eme_numerics.py

@@ -2,6 +2,7 @@ from typing import cast
 
 import numpy
 import pytest
+from scipy.sparse.linalg import spsolve
 
 from ..fdmath import vec
 from ..fdfd import eme, waveguide_2d, waveguide_3d, waveguide_cyl
@@ -547,6 +548,50 @@ def test_get_s_fdfd_from_2d_ports_identical_guides_are_nearly_transmissive() ->
     assert abs(abs(ss[0, 1]) - 1.0) < 0.03
 
 
+def test_multimode_fdfd_calibrates_the_full_incoming_matrix(monkeypatch) -> None:
+    # Exactly degenerate polarizations have different source/monitor bases.
+    # A uniform guide must still have zero interface reflection and conserve flux.
+    original_solve = waveguide_2d.solve_mode
+    seed = 0
+
+    def seeded_solve(*args, **kwargs):
+        nonlocal seed
+        seed += 1
+        return original_solve(*args, **kwargs, rng=seed)
+
+    decompositions = []
+    original_decomposition = waveguide_3d._port_decomposition
+
+    def capture_decomposition(*args, **kwargs):
+        result = original_decomposition(*args, **kwargs)
+        decompositions.append(result)
+        return result
+
+    monkeypatch.setattr(waveguide_2d, 'solve_mode', seeded_solve)
+    monkeypatch.setattr(waveguide_3d, '_port_decomposition', capture_decomposition)
+    epsilon = numpy.ones((3, 2, 2))
+    widths = [[numpy.full(2, 0.5), numpy.full(2, 0.5)] for _ in range(2)]
+    ss, diagnostics = waveguide_3d.get_s_fdfd_from_2d_ports(
+        epsilon, epsilon, omega=0.5, dxes_2d=widths, mode_numbers=[0, 1],
+        pml_cells=4, pml_buffer_cells=1, source_monitor_cells=2, monitor_interface_cells=2,
+        source_scales=[2, 0.3, 1j, 0.7j], matrix_solver=spsolve, return_diagnostics=True,
+    )
+    incoming, outgoing = [], []
+    for field in diagnostics['electric_fields']:
+        li, lo = waveguide_3d._decompose_port_field(decompositions[0], field)
+        ri, ro = waveguide_3d._decompose_port_field(decompositions[1], field)
+        incoming.append(numpy.concatenate((li, ri)))
+        outgoing.append(numpy.concatenate((lo, ro)))
+    incoming, outgoing = numpy.column_stack(incoming), numpy.column_stack(outgoing)
+    assert numpy.linalg.norm(incoming - numpy.diag(numpy.diag(incoming))) > 0.1
+    assert_close(ss @ incoming, outgoing, atol=1e-11, rtol=1e-11)
+    assert_close(diagnostics['incoming_amplitudes'], incoming, atol=1e-12)
+    assert_close(ss[:2, :2], 0, atol=1e-11)
+    assert_close(ss[2:, 2:], 0, atol=1e-11)
+    assert_close(ss.conj().T @ diagnostics['outgoing_power_gram'] @ ss, diagnostics['incoming_power_gram'], atol=1e-11)
+    assert_close(diagnostics['guided_powers'], 1, atol=1e-11)
+
+
 def test_get_s_fdfd_from_2d_ports_mismatch_returns_finite_diagnostics() -> None:
     left = _fdfd_cross_section(2.25)
     right = _fdfd_cross_section(1.5)

meanas/test/test_waveguide_2d_modes.py

@@ -1261,6 +1261,7 @@ def test_richer_complex_basis_reduces_driven_fdfd_monitor_residual() -> None:
         dxes_2d=dxes_2d,
         mode_numbers=[0],
         axis=2,
+        prop_dx=0.15,  # Resolve both port modes below the longitudinal band edge.
         pml_cells=4,
         pml_buffer_cells=2,
         source_monitor_cells=4,

meanas/test/test_waveguide_2d_numerics.py

@@ -1,10 +1,11 @@
 import numpy
+import pytest
 from numpy.linalg import norm
 from numpy.testing import assert_allclose
 from scipy import linalg
 
 from ..fdmath import vec
-from ..fdfd import waveguide_2d
+from ..fdfd import waveguide_2d, scpml
 
 
 OMEGA = 1 / 1500
@@ -121,7 +122,8 @@ def test_waveguide_2d_sensitivity_matches_finite_difference() -> None:
     assert 0.4 < ratio < 1.8
 
 
-def test_sensitivity_matches_central_differences_on_graded_yee_mesh() -> None:
+@pytest.mark.parametrize('absorption', ['none', 'pml', 'material'])
+def test_sensitivity_matches_central_differences_on_graded_yee_mesh(absorption) -> None:
     rng = numpy.random.default_rng(4)
     shape = (4, 3)
     omega = 0.7
@@ -129,6 +131,10 @@ def test_sensitivity_matches_central_differences_on_graded_yee_mesh() -> None:
     magnetic_widths = [rng.uniform(0.2, 1.2, n) for n in shape]
     electric_widths = [0.5 * (dx + numpy.roll(dx, -1)) for dx in magnetic_widths]
     dxes = [electric_widths, magnetic_widths]
+    if absorption == 'pml':
+        dxes = scpml.stretch_with_scpml(dxes, axis=0, polarity=-1, omega=omega, thickness=1, max_imag_ratio=0.3)
+    elif absorption == 'material':
+        epsilon = epsilon * (1 - 0.15j * rng.uniform(0.5, 1.5, epsilon.size))
     values, vectors = linalg.eig(waveguide_2d.operator_e(omega, dxes, epsilon).toarray())
     selected = numpy.argmax(values.real)
     beta = numpy.sqrt(values[selected])

meanas/test/test_waveguide_mode_helpers.py

@@ -153,8 +153,8 @@ def test_waveguide_3d_solve_mode_and_expand_e_are_phase_consistent(
     line_selection[axis] = slice(None)
     values = expanded[component][tuple(line_selection)]
     ratios = values[1:] / values[:-1]
-    phase_positions = numpy.real(dxes[0][axis]).cumsum()
+    intervals = numpy.real(dxes[1][axis][1:] if component == axis else dxes[0][axis][:-1])
-    expected_ratios = numpy.exp(-1j * result['wavenumber'] * numpy.diff(phase_positions))
+    expected_ratios = numpy.exp(-1j * result['wavenumber'] * intervals)
 
     numpy.testing.assert_allclose(ratios, expected_ratios, rtol=1e-6, atol=1e-9)
 
@@ -422,7 +422,7 @@ def test_waveguide_cyl_reconstruction_is_finite_radius_consistent(use_mu: bool)
     eps_r, eps_y, eps_z = numpy.split(epsilon, 3)
     mu_r, mu_y, mu_z = numpy.split(mu_full, 3)
 
-    ez_rhs = ta * (
+    ez_rhs = (
         (Dbx @ (tb * eps_r * er))
         + (Dby @ (ta * eps_y * ey))
         ) / eps_z
@@ -469,8 +469,79 @@ def test_waveguide_cyl_linear_wavenumbers_use_component_radii() -> None:
         rmin=rmin,
         )
 
-    numpy.testing.assert_allclose(er_linear[0], angular_wavenumber[0] / rb.mean())
+    numpy.testing.assert_allclose(er_linear[0], angular_wavenumber[0] / numpy.average(rb, weights=dxes[0][0]))
-    numpy.testing.assert_allclose(ey_linear[0], angular_wavenumber[0] / ra.mean())
+    numpy.testing.assert_allclose(ey_linear[0], angular_wavenumber[0] / numpy.average(ra, weights=dxes[1][0]))
+
+
+@pytest.mark.parametrize('use_mu', [False, True])
+def test_cylindrical_modes_obey_the_full_ampere_curl(use_mu) -> None:
+    shape, omega, rmin = (4, 3), 0.7, 2.0
+    dxes = [[numpy.full(n, 0.4) for n in shape] for _ in range(2)]
+    rng = numpy.random.default_rng(4)
+    epsilon = vec(1 + rng.uniform(0, 3, (3, *shape)))
+    mu = vec(1 + rng.uniform(0, 2, (3, *shape))) if use_mu else numpy.ones_like(epsilon)
+    matrix = waveguide_cyl.cylindrical_operator(omega, dxes, epsilon, rmin, mu).toarray()
+    values, vectors = linalg.eig(matrix)
+    selected = numpy.argmax(values.real)
+    beta = numpy.sqrt(values[selected])
+    electric = waveguide_cyl.exy2e(beta * rmin, dxes, rmin, epsilon) @ vectors[:, selected]
+    magnetic = waveguide_cyl.e2h(beta * rmin, omega, dxes, rmin, mu) @ electric
+    ta, tb = (metric.diagonal() for metric in waveguide_cyl.dxes2T(dxes, rmin))
+    dbx, dby = deriv_back(dxes[1])
+    hr, hy, hz = numpy.split(magnetic, 3)
+    # Physical cylindrical Ampere curl, independently of electric reconstruction.
+    curl_h = numpy.concatenate((
+        1j * beta * hy / tb + dby @ hz,
+        -1j * beta * hr / ta - (dbx @ (tb * hz)) / ta,
+        -dby @ hr + dbx @ hy,
+    ))
+    numpy.testing.assert_allclose(curl_h, 1j * omega * epsilon * electric, atol=1e-11, rtol=1e-10)
+
+
+def test_cylindrical_mean_radius_is_invariant_under_constant_density_refinement() -> None:
+    for widths in ([1.0, 3.0], [1.0, 1.0, 1.0, 1.0]):
+        n = len(widths)
+        dxes = [[numpy.array(widths), numpy.ones(1)] for _ in range(2)]
+        electric_xy = numpy.concatenate((numpy.ones(n), numpy.zeros(n)))[None, :].astype(complex)
+        beta = waveguide_cyl.linear_wavenumbers(electric_xy, [12.0], numpy.ones(3 * n), dxes, rmin=10.0)
+        numpy.testing.assert_allclose(beta, 1.0)
+
+
+@pytest.mark.parametrize('axis', [0, 1, 2])
+@pytest.mark.parametrize('polarity', [-1, 1])
+def test_graded_expansion_uses_each_components_actual_yee_intervals(axis, polarity) -> None:
+    shape = [1, 1, 1]
+    shape[axis] = 4
+    dxes = [[numpy.ones(n) for n in shape] for _ in range(2)]
+    dxes[0][axis] = numpy.array([1.0, 2.0, 3.0, 4.0])
+    dxes[1][axis] = numpy.array([0.6, 0.8, 1.0, 1.2])
+    source = [slice(None)] * 3
+    source[axis] = slice(1, 2)
+    field = numpy.zeros((3, *shape), dtype=complex)
+    field[(slice(None), *source)] = 1 + 0.5j
+    expanded = waveguide_3d.expand_e(field, 0.2, dxes, axis, polarity, source)
+    for component in range(3):
+        line = expanded[component].ravel()
+        intervals = dxes[1][axis][1:] if component == axis else dxes[0][axis][:-1]
+        numpy.testing.assert_allclose(line[1:] / line[:-1], numpy.exp(-1j * polarity * 0.2 * intervals))
+    numpy.testing.assert_array_equal(expanded[(slice(None), *source)], field[(slice(None), *source)])
+
+
+@pytest.mark.parametrize('beta', [2.0, 3.0])
+def test_grid_band_edge_rejects_power_normalization_and_selects_decay(beta) -> None:
+    shape = (1, 1, 5)
+    dxes = [[numpy.ones(n) for n in shape] for _ in range(2)]
+    kwargs = dict(omega=1.0, dxes=dxes, axis=2, polarity=1,
+                  slices=(slice(None), slice(None), slice(0, 1)), epsilon=numpy.full((3, *shape), beta ** 2))
+    with pytest.raises(ValueError, match='grid band edge'):
+        waveguide_3d.mode_from_2d(numpy.array([1, 0], dtype=complex), beta, **kwargs)
+    mode = waveguide_3d.mode_from_2d(numpy.array([1, 0], dtype=complex), beta, normalization='lorentz', **kwargs)
+    assert mode['wavenumber'].imag <= 0
+    expanded = waveguide_3d.expand_e(mode['E'], mode['wavenumber'], dxes, 2, 1, kwargs['slices'])
+    ratio = abs(expanded[0, 0, 0, 1] / expanded[0, 0, 0, 0])
+    assert ratio <= 1 + 1e-14
+    if beta > 2:
+        assert ratio < 1
 
 
 def test_waveguide_cyl_exy2e_and_exy2h_return_finite_full_fields() -> None: