[docs] describe corrected FDFD solver and modal conventions

docs/api/fdfd.md

@@ -14,6 +14,47 @@
 
 ::: meanas.fdfd.farfield
 
+## Bloch magnetic materials
+
+Bloch state vectors store the transverse Fourier coefficients of magnetic flux
+`B = mu * H`. They equal `H` when permeability is omitted. Pass the same `mu`
+to `bloch.hmn_2_hxyz(...)` and `bloch.hmn_2_exyz(...)` when converting solved
+vectors to fields. The eigensolver uses the projected `mu**-1` energy metric
+for finite positive real permeability, including anisotropic and spatially
+varying profiles. Its preconditioner remains finite at the Gamma point.
+
+## Far-field normalization
+
+The near-to-far transform includes the cell area `dx * dy`, so normalized
+far-field amplitudes describe a Fourier integral and converge under spatial
+refinement. This changes amplitudes previously returned as unweighted FFT sums.
+The inverse applies the matching Fourier quadrature and returns spacings for
+the requested output shape. Both transforms center their origins at index
+`size // 2`; padding and cropping preserve that sample. The inverse does not
+modify input spectra.
+
+Transforms use the library's `exp(+i omega t)` time convention: outgoing plane
+waves vary as `exp(-i k dot r)`, so the angular projection uses a positive
+Fourier phase. The inverse uses a negative phase. Far-field vector pairs are
+ordered as polar and azimuthal components; near-field pairs are Cartesian.
+
+## Solver and port calibration
+
+The default QMR wrapper raises `RuntimeError` on breakdown or failure to
+converge. PEC rows impose zero electric fields with unit diagonal entries;
+`solvers.generic` zeroes the corresponding source and initial-guess entries.
+H-only and EH operators impose analogous constraints for PMC locations.
+
+FDFD port calibration uses the complete incoming-amplitude matrix from all
+excitations, solving `S @ incoming_amplitudes = raw_outgoing_amplitudes`.
+Diagnostics include these matrices and the incoming/outgoing flux Gram
+matrices. Guided powers retain cross terms between nonorthogonal modes.
+
+Lossless modes at or above the propagation-grid band edge require a finer
+mesh for power normalization. Lorentz-normalized evanescent modes use the
+decaying spatial branch. Expansion on graded grids uses each electric
+component's own Yee intervals.
+
 ## SCPML conventions
 
 `stretch_with_scpml` now uses `-1` for the low-index face and `+1` for the

meanas/fdfd/bloch.py

@@ -1,97 +1,44 @@
 """
-Bloch eigenmode solver/operators
+Bloch eigenmode solvers and operators in a spatial Fourier basis.
 
-This module contains functions for generating and solving the
+The two stored coefficients at each wavevector represent transverse magnetic
- 3D Bloch eigenproblem. The approach is to transform the problem
+flux B = mu H. They coincide with H when mu is omitted. Using B ensures
- into the (spatial) fourier domain, transforming the equation
+k dot B = 0 even in anisotropic or spatially varying magnetic materials; H
+itself need not be transverse. The historical h_mn function names are retained.
+Pass the same mu to hmn_2_hxyz and hmn_2_exyz when reconstructing these modes.
 
-    1/mu * curl(1/eps * curl(H_eigenmode)) = (w/c)^2 H_eigenmode
+For each Fourier wavevector, generate_kmn constructs unit vectors m and n with
+k cross m = |k| n. The Cartesian Fourier field is B_k = b_m m + b_n n.
+The Maxwell operator applies, in order:
 
- into
+    H_k = fftn(ifftn(B_k) / mu)
+    E_k = fftn(ifftn(k cross H_k) / epsilon)
+    A(B_k) = k cross E_k
 
-    conv(1/mu_k, ik x conv(1/eps_k, ik x H_k)) = (w/c)^2 H_k
+Material division is performed in real space and curls in Fourier space.
+The explicit factors of i are omitted, giving the historical negative squared
+frequency eigenvalues. Wavevectors and frequencies use cycles per length:
+frequency = sqrt(-real(eigenvalue)). Reconstructed electric fields retain the
+historical overall frequency scaling and are not power-normalized.
 
- where:
+For finite positive real permeability, the flux operator is self-adjoint in
+the projected magnetic energy metric M = P mu^-1 P, where P projects onto the
+transverse Fourier subspace. The block variational solve uses that metric in
+its Rayleigh quotients, gradients, line searches, and final Ritz extraction.
+The inverse-curl preconditioner regularizes zero wavevectors so static modes
+remain accessible at Gamma. Nonzero-frequency modes have no static Fourier
+flux component at Gamma; the two stored zero-wavevector basis vectors choose
+two representatives of the static sector.
 
-  - the `_k` subscript denotes a 3D fourier transformed field
+Typical usage:
-  - each component of `H_k` corresponds to a plane wave with wavevector `k`
-  - `x` is the cross product
-  - `conv()` denotes convolution
 
- Since `k` and `H` are orthogonal for each plane wave, we can use each
+    values, vectors = eigsolve(5, k0, reciprocal_lattice, epsilon, mu=mu)
- `k` to create an orthogonal basis (k, m, n), with `k x m = n`, and
+    frequency = numpy.sqrt(-numpy.real(values[0]))
- `|m| = |n| = 1`. The cross products are then simplified as follows:
+    electric = hmn_2_exyz(k0, reciprocal_lattice, epsilon, mu=mu)(vectors[0])
-
+    magnetic = hmn_2_hxyz(k0, reciprocal_lattice, epsilon, mu=mu)(vectors[0])
-  - `h` is shorthand for `H_k`
-  - `(...)_xyz` denotes the `(x, y, z)` basis
-  - `(...)_kmn` denotes the `(k, m, n)` basis
-  - `hm` is the component of `h` in the `m` direction, etc.
-
-  We know
-
-    k @ h = kx hx + ky hy + kz hz = 0 = hk
-    h = hk + hm + hn = hm + hn
-    k = kk + km + kn = kk = |k|
-
-  We can write
-
-    k x h = (ky hz - kz hy,
-             kz hx - kx hz,
-             kx hy - ky hx)_xyz
-          = ((k x h) @ k, (k x h) @ m, (k x h) @ n)_kmn
-          = (0, (m x k) @ h, (n x k) @ h)_kmn         # triple product ordering
-          = (0, kk (-n @ h), kk (m @ h))_kmn          # (m x k) = -|k| n, etc.
-          = |k| (0, -h @ n, h @ m)_kmn
-
-  which gives us a straightforward way to perform the cross product
-  while simultaneously transforming into the `_kmn` basis.
-  We can also write
-
-    k x h = (km hn - kn hm,
-             kn hk - kk hn,
-             kk hm - km hk)_kmn
-          = (0, -kk hn, kk hm)_kmn
-          = (-kk hn)(mx, my, mz)_xyz + (kk hm)(nx, ny, nz)_xyz
-          = |k| (hm * (nx, ny, nz)_xyz
-               - hn * (mx, my, mz)_xyz)
-
-  which gives us a way to perform the cross product while simultaneously
-  trasnforming back into the `_xyz` basis.
-
- We can also simplify `conv(X_k, Y_k)` as `fftn(X * ifftn(Y_k))`.
-
- Using these results and storing `H_k` as `h = (hm, hn)`, we have
-
-    e_xyz = fftn(1/eps * ifftn(|k| (hm * n - hn * m)))
-    b_mn = |k| (-e_xyz @ n, e_xyz @ m)
-    h_mn = fftn(1/mu * ifftn(b_m * m + b_n * n))
-
- which forms the operator from the left side of the equation.
-
- We can then use a preconditioned block Rayleigh iteration algorithm, as in
-  SG Johnson and JD Joannopoulos, Block-iterative frequency-domain methods
-  for Maxwell's equations in a planewave basis, Optics Express 8, 3, 173-190 (2001)
- (similar to that used in MPB) to find the eigenvectors for this operator.
-
- ===
-
- Typically you will want to do something like
-
-    recip_lattice = numpy.diag(1/numpy.array(epsilon[0].shape * dx))
-    n, v = bloch.eigsolve(5, k0, recip_lattice, epsilon)
-    f = numpy.sqrt(-numpy.real(n[0]))
-    n_eff = norm(recip_lattice @ k0) / f
-
-    v2e = bloch.hmn_2_exyz(k0, recip_lattice, epsilon)
-    e_field = v2e(v[0])
-
-    k, f = find_k(frequency=1/1550,
-                  tolerance=(1/1550 - 1/1551),
-                  direction=[1, 0, 0],
-                  G_matrix=recip_lattice,
-                  epsilon=epsilon,
-                  band=0)
 
+The block iteration follows the approach of S. G. Johnson and J. D.
+Joannopoulos, Optics Express 8, 173-190 (2001).
 """
 
 from typing import Any, cast
@@ -174,7 +121,7 @@ def generate_kmn(
     m = numpy.broadcast_to([0, 1, 0], tuple(shape[:3]) + (3,)).astype(float)
     n = numpy.broadcast_to([0, 0, 1], tuple(shape[:3]) + (3,)).astype(float)
 
-    xy_non0 = numpy.any(k_xyz[:, :, :, 0:1] != 0, axis=3)
+    xy_non0 = numpy.any(k_xyz[:, :, :, :2] != 0, axis=3)
     if numpy.any(xy_non0):
         u = numpy.cross(k_xyz[xy_non0], [0, 0, 1])
         m[xy_non0, :] = u / norm(u, axis=1)[:, None]
@@ -195,13 +142,20 @@ def maxwell_operator(
         mu: fdfield | None = None
         ) -> Callable[[NDArray[numpy.complex128]], NDArray[numpy.complex128]]:
     """
-    Generate the Maxwell operator
+    Generate the Maxwell operator on transverse magnetic-flux coefficients.
 
-        conv(1/mu_k, ik x conv(1/eps_k, ik x ___))
+        k x conv(1/eps_k, k x conv(1/mu_k, ___))
 
-    which is the spatial-frequency-space representation of
+    The corresponding real-space curl-curl problem is
 
-        1/mu * curl(1/eps * curl(___))
+        curl(1/eps * curl(1/mu * ___))
+
+    with the overall factor ``-(2*pi)**2`` omitted for the module's
+    cycles-per-length frequency convention.
+
+    The state represents B = mu H, whose Fourier components are transverse.
+    It equals H when mu is omitted. Eigenvalues retain the historical negative
+    squared-frequency convention (the two explicit factors of i are omitted).
 
     The operator is a function that acts on a vector h_mn of size `2 * epsilon[0].size`
 
@@ -230,23 +184,28 @@ def maxwell_operator(
         """
         Maxwell operator for Bloch eigenmode simulation.
 
-        h is complex 2-field in (m, n) basis, vectorized
+        h contains the two transverse magnetic-flux coefficients.
 
         Args:
-            h: Raveled h_mn; size `2 * epsilon[0].size`.
+            h: Raveled b_mn; size `2 * epsilon[0].size`. Not modified.
-                Altered in-place.
 
         Returns:
-            Raveled conv(1/mu_k, ik x conv(1/eps_k, ik x h_mn)), returned
+            Raveled transverse flux image under the Maxwell operator.
-            and overwritten in-place of `h`.
         """
         hin_m, hin_n = (hi.reshape(shape) for hi in numpy.split(h, 2))
 
         #{d,e,h}_xyz fields are complex 3-fields in (1/x, 1/y, 1/z) basis
 
-        # cross product and transform into xyz basis
+        # Convert transverse B into H before taking either curl. H itself need
-        d_xyz = (n * hin_m
+        # not be transverse for anisotropic or spatially varying permeability.
-               - m * hin_n) * k_mag
+        if mu is None:
+            h_m, h_n = hin_m, hin_n
+        else:
+            b_xyz = m * hin_m + n * hin_n
+            h_xyz = fftn(ifftn(b_xyz, axes=range(3)) / mu, axes=range(3))
+            h_m = numpy.sum(h_xyz * m, axis=3, keepdims=True)
+            h_n = numpy.sum(h_xyz * n, axis=3, keepdims=True)
+        d_xyz = (n * h_m - m * h_n) * k_mag
 
         # divide by epsilon
         temp = ifftn(d_xyz, axes=range(3))   # reuses d_xyz if using pyfftw
@@ -257,23 +216,7 @@ def maxwell_operator(
         b_m = numpy.sum(e_xyz * n, axis=3, keepdims=True) * -k_mag
         b_n = numpy.sum(e_xyz * m, axis=3, keepdims=True) * +k_mag
 
-        if mu is None:
+        return _assemble_hmn_vector(b_m, b_n)
-            h_m, h_n = b_m, b_n
-        else:
-            # transform from mn to xyz
-            b_xyz = (m * b_m
-                   + n * b_n)    # noqa
-
-            # divide by mu
-            temp = ifftn(b_xyz, axes=range(3))
-            temp /= mu
-            h_xyz = fftn(temp, axes=range(3))
-
-            # transform back to mn
-            h_m = numpy.sum(h_xyz * m, axis=3)
-            h_n = numpy.sum(h_xyz * n, axis=3)
-
-        return _assemble_hmn_vector(h_m, h_n)
 
     return operator
 
@@ -282,12 +225,12 @@ def hmn_2_exyz(
         k0: ArrayLike,
         G_matrix: ArrayLike,
         epsilon: fdfield,
+        mu: fdfield | None = None,
         ) -> Callable[[NDArray[numpy.complex128]], cfdfield_t]:
     """
-    Generate an operator which converts a vectorized spatial-frequency-space
+    Convert transverse Fourier flux coefficients into an unnormalized E field.
-     `h_mn` into an E-field distribution, i.e.
 
-        ifft(conv(1/eps_k, ik x h_mn))
+        ifft(conv(1/eps_k, k x conv(1/mu_k, b_mn)))
 
     The operator is a function that acts on a vector `h_mn` of size `2 * epsilon[0].size`.
 
@@ -298,17 +241,26 @@ def hmn_2_exyz(
         G_matrix: 3x3 matrix, with reciprocal lattice vectors as columns.
         epsilon: Dielectric constant distribution for the simulation.
                  All fields are sampled at cell centers (i.e., NOT Yee-gridded)
+        mu: Permeability used by the solve. The input coefficients describe B;
+            divide by mu to recover H before computing the electric field.
 
     Returns:
         Function for converting `h_mn` into `E_xyz`
     """
     shape = epsilon[0].shape + (1,)
     epsilon = numpy.moveaxis(epsilon, 0, -1)
+    if mu is not None:
+        mu = numpy.moveaxis(mu, 0, -1)
 
     k_mag, m, n = generate_kmn(k0, G_matrix, shape)
 
     def operator(h: NDArray[numpy.complex128]) -> cfdfield_t:
         hin_m, hin_n = (hi.reshape(shape) for hi in numpy.split(h, 2))
+        if mu is not None:
+            b_xyz = m * hin_m + n * hin_n
+            h_xyz = fftn(ifftn(b_xyz, axes=range(3)) / mu, axes=range(3))
+            hin_m = numpy.sum(h_xyz * m, axis=3, keepdims=True)
+            hin_n = numpy.sum(h_xyz * n, axis=3, keepdims=True)
         d_xyz = (n * hin_m
                - m * hin_n) * k_mag
 
@@ -323,12 +275,12 @@ def hmn_2_hxyz(
         k0: ArrayLike,
         G_matrix: ArrayLike,
         epsilon: fdfield,
+        mu: fdfield | None = None,
         ) -> Callable[[NDArray[numpy.complex128]], cfdfield_t]:
     """
-    Generate an operator which converts a vectorized spatial-frequency-space
+    Convert transverse Fourier flux coefficients into a physical H field.
-     `h_mn` into an H-field distribution, i.e.
 
-        ifft(h_mn)
+        ifft(b_mn) / mu
 
     The operator is a function that acts on a vector `h_mn` of size `2 * epsilon[0].size`.
 
@@ -339,6 +291,8 @@ def hmn_2_hxyz(
         G_matrix: 3x3 matrix, with reciprocal lattice vectors as columns.
         epsilon: Dielectric constant distribution for the simulation.
                  Only `epsilon[0].shape` is used.
+        mu: Permeability used by the solve. Input coefficients describe B;
+            the returned physical H is B / mu. Default is unit permeability.
 
     Returns:
         Function for converting `h_mn` into `H_xyz`
@@ -350,7 +304,10 @@ def hmn_2_hxyz(
         hin_m, hin_n = (hi.reshape(shape) for hi in numpy.split(h, 2))
         h_xyz = (m * hin_m
                + n * hin_n)
-        return cfdfield_t(numpy.array([ifftn(hi) for hi in numpy.moveaxis(h_xyz, 3, 0)]))
+        field = numpy.array([ifftn(hi) for hi in numpy.moveaxis(h_xyz, 3, 0)])
+        if mu is not None:
+            field /= mu
+        return cfdfield_t(field)
 
     return operator
 
@@ -362,11 +319,16 @@ def inverse_maxwell_operator_approx(
         mu: fdfield | None = None,
         ) -> Callable[[NDArray[numpy.complex128]], NDArray[numpy.complex128]]:
     """
-    Generate an approximate inverse of the Maxwell operator,
+    Generate a positive approximate inverse of minus the flux Maxwell operator,
 
-        ik x conv(eps_k, ik x conv(mu_k, ___))
+        project_mn(conv(mu_k, cross_inv(conv(eps_k, cross_inv(___)))))
 
-     which can be used to improve the speed of ARPACK in shift-invert mode.
+    At k = 0, inverse crosses use the smallest nonzero grid wavevector scale
+    (or one on a single-point grid). This finite regularization keeps the
+    preconditioner positive on static directions, which must remain accessible
+    to the eigensolver at Gamma.
+
+    It is used for the preconditioned block variational solve.
 
     See the `meanas.fdfd.bloch` docstring for more information.
 
@@ -387,6 +349,9 @@ def inverse_maxwell_operator_approx(
         mu = numpy.moveaxis(mu, 0, -1)
 
     k_mag, m, n = generate_kmn(k0, G_matrix, shape)
+    nonzero = k_mag[k_mag != 0]
+    reference = float(nonzero.min()) if nonzero.size else 1.0
+    k_inverse = numpy.divide(1, k_mag, out=numpy.full_like(k_mag, 1 / reference), where=k_mag != 0)
 
     def operator(h: NDArray[numpy.complex128]) -> NDArray[numpy.complex128]:
         """
@@ -398,31 +363,14 @@ def inverse_maxwell_operator_approx(
             h: Raveled h_mn; size `2 * epsilon[0].size`.
 
         Returns:
-            Raveled ik x conv(eps_k, ik x conv(mu_k, h_mn))
+            Raveled, regularized approximate inverse flux image.
         """
         hin_m, hin_n = (hi.reshape(shape) for hi in numpy.split(h, 2))
 
         #{d,e,h}_xyz fields are complex 3-fields in (1/x, 1/y, 1/z) basis
 
-        if mu is None:
-            b_m, b_n = hin_m, hin_n
-        else:
-            # transform from mn to xyz
-            h_xyz = (m * hin_m
-                   + n * hin_n)  # noqa
-
-            # multiply by mu
-            temp = ifftn(h_xyz, axes=range(3))
-            temp *= mu
-            b_xyz = fftn(temp, axes=range(3))
-
-            # transform back to mn
-            b_m = numpy.sum(b_xyz * m, axis=3, keepdims=True)
-            b_n = numpy.sum(b_xyz * n, axis=3, keepdims=True)
-
         # cross product and transform into xyz basis
-        e_xyz = (n * b_m
+        e_xyz = (n * hin_m - m * hin_n) * k_inverse
-               - m * b_n) / k_mag
 
         # multiply by epsilon
         temp = ifftn(e_xyz, axes=range(3))
@@ -430,14 +378,41 @@ def inverse_maxwell_operator_approx(
         d_xyz = fftn(temp, axes=range(3))
 
         # cross product and transform into mn basis   crossinv_t2c
-        h_m = numpy.sum(d_xyz * n, axis=3, keepdims=True) / +k_mag
+        h_m = numpy.sum(d_xyz * n, axis=3, keepdims=True) * k_inverse
-        h_n = numpy.sum(d_xyz * m, axis=3, keepdims=True) / -k_mag
+        h_n = numpy.sum(d_xyz * m, axis=3, keepdims=True) * -k_inverse
+        if mu is not None:
+            h_xyz = m * h_m + n * h_n
+            b_xyz = fftn(ifftn(h_xyz, axes=range(3)) * mu, axes=range(3))
+            h_m = numpy.sum(b_xyz * m, axis=3, keepdims=True)
+            h_n = numpy.sum(b_xyz * n, axis=3, keepdims=True)
 
         return _assemble_hmn_vector(h_m, h_n)
 
     return operator
 
 
+def _transverse_material_operator(
+        k0: ArrayLike,
+        G_matrix: ArrayLike,
+        material: fdfield,
+        *,
+        inverse: bool,
+        ) -> Callable[[NDArray], NDArray]:
+    """Apply P material^(+/-1) P on the transverse Fourier subspace."""
+    shape = material[0].shape + (1,)
+    _k_mag, m, n = generate_kmn(k0, G_matrix, shape)
+    values = numpy.moveaxis(material, 0, -1)
+
+    def apply(h: NDArray) -> NDArray:
+        hm, hn = (part.reshape(shape) for part in numpy.split(h, 2))
+        field = ifftn(m * hm + n * hn, axes=range(3))
+        field = field / values if inverse else field * values
+        transformed = fftn(field, axes=range(3))
+        return _assemble_hmn_vector(numpy.sum(transformed * m, axis=3), numpy.sum(transformed * n, axis=3))
+
+    return apply
+
+
 def find_k(
         frequency: float,
         tolerance: float,
@@ -528,18 +503,20 @@ def eigsolve(
         G_matrix: 3x3 matrix, with reciprocal lattice vectors as columns.
         epsilon: Dielectric constant distribution for the simulation.
                  All fields are sampled at cell centers (i.e., NOT Yee-gridded)
-        mu: Magnetic permability distribution for the simulation.
+        mu: Finite positive real diagonal permeability. Default is one.
-            Default `None` (1 everywhere).
+        tolerance: Relative projected-gradient tolerance in the magnetic energy
-        tolerance: Solver stops when fractional change in the objective
+            metric. Flat objective traces trigger a conjugate-direction restart.
-                   `trace(Z.H @ A @ Z @ inv(Z Z.H))` is smaller than the tolerance
+        max_iters: Maximum iteration count; returns Ritz pairs from the final
-        max_iters: TODO
+            subspace if the limit is reached.
-        reset_iters: TODO
+        reset_iters: Periodic conjugate-direction restart interval.
-        callback: TODO
+        callback: Optional no-argument callback after each subspace update.
-        y0: TODO, initial guess
+        y0: Initial flux coefficients with shape (num_modes, 2 * grid_size).
+            Copied before use.
 
     Returns:
         `(eigenvalues, eigenvectors)` where `eigenvalues[i]` corresponds to the
-        vector `eigenvectors[i, :]`
+        transverse flux vector `eigenvectors[i, :]`. Pass mu to the field
+        converters when it was supplied to this solve.
     """
     k0 = numpy.asarray(k0)
 
@@ -556,6 +533,24 @@ def eigsolve(
     scipy_op = spalg.LinearOperator(dtype=complex, shape=(h_size, h_size), matvec=mop)
     scipy_iop = spalg.LinearOperator(dtype=complex, shape=(h_size, h_size), matvec=imop)
 
+    metric = inverse_metric = None
+    if mu is not None:
+        if numpy.iscomplexobj(mu) or not numpy.isfinite(mu).all() or numpy.any(mu <= 0):
+            raise ValueError('Bloch variational solving requires finite positive real permeability')
+        metric = spalg.LinearOperator(
+            shape=(h_size, h_size), dtype=complex,
+            matvec=_transverse_material_operator(k0, G_matrix, mu, inverse=True),
+        )
+        inverse_metric = spalg.LinearOperator(
+            shape=(h_size, h_size), dtype=complex,
+            matvec=_transverse_material_operator(k0, G_matrix, mu, inverse=False),
+        )
+
+    def mass(vectors: NDArray) -> NDArray:
+        # P mu^-1 P makes the flux operator self-adjoint in the magnetic
+        # energy inner product, even though it is not Euclidean-Hermitian.
+        return vectors if metric is None else metric @ vectors
+
     y_shape = (h_size, num_modes)
 
     prev_E = 0.0
@@ -569,11 +564,11 @@ def eigsolve(
     if y0 is None:
         Z = rng.random(y_shape) + 1j * rng.random(y_shape)
     else:
-        Z = numpy.asarray(y0).T
+        Z = numpy.array(y0, dtype=complex, copy=True).T
 
     while True:
         Z *= num_modes / norm(Z)
-        ZtZ = Z.conj().T @ Z
+        ZtZ = Z.conj().T @ mass(Z)
         try:
             U = numpy.linalg.inv(ZtZ)
         except numpy.linalg.LinAlgError:
@@ -592,27 +587,28 @@ def eigsolve(
 
     for i in range(max_iters):
         Zt = numpy.conj(Z.T, out=Zt)
-        ZtZ = Zt @ Z
+        ZtZ = Zt @ mass(Z)
         U = numpy.linalg.inv(ZtZ)
         AZ = scipy_op @ Z.copy()
-        ZtAZ = Zt @ AZ
+        ZtAZ = Zt @ mass(AZ)
         ZtAZU = ZtAZ @ U
         E_signed = real(trace(ZtAZU))
         sgn = numpy.sign(E_signed)
         E = numpy.abs(E_signed)
 
         # G = AZU projected onto the space orthonormal to Z via (1 - ZUZt)
-        G = (AZ @ U - Z @ U @ ZtAZU) * sgn
+        G = mass(AZ @ U - Z @ U @ ZtAZU) * sgn
 
         if i > 0 and abs(E - prev_E) < tolerance * 0.5 * (E + prev_E + 1e-7):
-            logger.info(
+            # A flat line search alone does not imply a stationary subspace.
-                'Optimization succeded: '
+            # Restart with steepest descent if conjugacy has been lost.
-                f'[change in trace] {abs(E - prev_E)} - 5e-8 '
+            prev_traceGtKG = 0
-                f'< {tolerance} [tolerance] * {(E + prev_E) / 2} [value of trace]'
-                )
-            break
 
-        KG = scipy_iop @ G          # Preconditioned steepest descent direction
+        if norm(G) <= tolerance * max(norm(mass(AZ) @ U), 1e-7):
+            break
+        # Approximate the inverse energy Hessian by a symmetric sandwich of
+        # P mu P around the inverse-curl / epsilon preconditioner.
+        KG = scipy_iop @ (G if inverse_metric is None else inverse_metric @ G)
         traceGtKG = _rtrace_AtB(G, KG)
 
         if prev_traceGtKG == 0 or i % reset_iters == 0:
@@ -642,11 +638,11 @@ def eigsolve(
         #   Qi = inv(Q) = U'
 
         AD = scipy_op @ D.copy()
-        DtD = D.conj().T @ D
+        DtD = D.conj().T @ mass(D)
-        DtAD = D.conj().T @ AD
+        DtAD = D.conj().T @ mass(AD)
 
-        symZtD = _symmetrize(Zt @ D)
+        symZtD = _symmetrize(Zt @ mass(D))
-        symZtAD = _symmetrize(Zt @ AD)
+        symZtAD = _symmetrize(Zt @ mass(AD))
 
         Qi_memo: list[float | None] = [None, None]
 
@@ -737,8 +733,10 @@ def eigsolve(
         result = scipy.optimize.minimize_scalar(
             trace_func,
             method='bounded',
-            bounds=(0, pi),
+            # Center the period on zero so the line search can resolve small
-            options={'xatol': tolerance},
+            # steps without the relative angular floor near theta = pi.
+            bounds=(-pi / 2, pi / 2),
+            options={'xatol': min(tolerance, 1e-12)},
             )
         new_E = result.fun
         theta = result.x
@@ -757,9 +755,9 @@ def eigsolve(
     #
     # Recover eigenvectors from Z
     #
-    U = numpy.linalg.inv(ZtZ)
+    U = numpy.linalg.inv(Z.conj().T @ mass(Z))
     Y = Z @ scipy.linalg.sqrtm(U).astype(numpy.complex128)
-    W = Y.conj().T @ (scipy_op @ Y.copy())
+    W = Y.conj().T @ mass(scipy_op @ Y.copy())
 
     eigvals, W_eigvecs = numpy.linalg.eig(W)
     eigvecs = Y @ W_eigvecs
@@ -769,10 +767,10 @@ def eigsolve(
         n = eigvals[i]
         v /= norm(v)
         Av = numpy.asarray(scipy_op @ v.copy()).reshape(-1)
-        eigness = norm(Av - (v.conj() @ Av) * v)
+        eigness = norm(Av - n * v)
-        f = numpy.sqrt(-numpy.real(n))
+        f = numpy.sqrt(max(0.0, -numpy.real(n)))
-        df = numpy.sqrt(-numpy.real(n) + eigness)
+        df = numpy.sqrt(max(0.0, -numpy.real(n) + eigness))
-        neff_err = kmag * (1 / df - 1 / f)
+        neff_err = kmag * (1 / df - 1 / f) if f > 0 and df > 0 else numpy.inf
         logger.info(f'eigness {i}: {eigness}\n neff_err: {neff_err}')
 
     order = numpy.argsort(numpy.abs(eigvals))

meanas/fdfd/farfield.py

@@ -1,8 +1,9 @@
 """
 Functions for performing near-to-farfield transformation (and the reverse).
 """
-from typing import Any, cast
+from typing import Any
 from collections.abc import Sequence
+import operator
 import numpy
 from numpy.fft import fft2, fftshift, fftfreq, ifft2, ifftshift
 from numpy import pi
@@ -13,6 +14,34 @@ type farfield_slice = NDArray[complexfloating]
 type transverse_slice_pair = Sequence[farfield_slice]
 
 
+def _padded_shape(shape: tuple[int, ...], padded_size: Sequence[int] | int | None) -> tuple[int, int]:
+    if len(shape) != 2 or any(n == 0 for n in shape):
+        raise ValueError('fields must be nonempty two-dimensional arrays')
+    if padded_size is None:
+        sizes = tuple(1 << (n - 1).bit_length() for n in shape)
+    else:
+        raw_sizes = (padded_size, padded_size) if numpy.isscalar(padded_size) else padded_size
+        sizes = tuple(operator.index(n) for n in raw_sizes)
+    if len(sizes) != 2 or any(n <= 0 for n in sizes):
+        raise ValueError('padded_size must contain two positive integer sizes')
+    return sizes[0], sizes[1]
+
+
+def _resize_centered(field: farfield_slice, shape: tuple[int, int]) -> farfield_slice:
+    """Pad or crop while preserving the sample at index size // 2 as the origin."""
+    if field.shape == shape:
+        return field
+    result = numpy.zeros(shape, dtype=field.dtype)
+    source, target = [], []
+    for old, new in zip(field.shape, shape, strict=True):
+        start_old, start_new = max(0, old // 2 - new // 2), max(0, new // 2 - old // 2)
+        count = min(old - start_old, new - start_new)
+        source.append(slice(start_old, start_old + count))
+        target.append(slice(start_new, start_new + count))
+    result[tuple(target)] = field[tuple(source)]
+    return result
+
+
 def near_to_farfield(
         E_near: transverse_slice_pair,
         H_near: transverse_slice_pair,
@@ -24,7 +53,11 @@ def near_to_farfield(
     Compute the farfield, i.e. the distribution of the fields after propagation
       through several wavelengths of uniform medium.
 
-    The input fields should be complex phasors.
+    The input fields should be complex phasors. The spatial origin is the
+    sample at ``(nx // 2, ny // 2)``. Padding/cropping preserves that origin.
+    Fourier integrals include the cell area ``dx * dy``.
+    Phasors use ``exp(+i omega t)``; an outgoing plane wave has spatial
+    dependence ``exp(-i k dot r)`` and is projected with a positive Fourier phase.
 
     Args:
         E_near: List of 2 ndarrays containing the 2D phasor field slices for the transverse
@@ -40,9 +73,9 @@ def near_to_farfield(
     Returns:
         Dict with keys
 
-        -   `E_far`: Normalized E-field farfield; multiply by
+        -   `E`: Normalized polar and azimuthal E-field components; multiply by
                 (i k exp(-i k r) / (4 pi r)) to get the actual field value.
-        -   `H_far`: Normalized H-field farfield; multiply by
+        -   `H`: Normalized polar and azimuthal H-field components; multiply by
                 (i k exp(-i k r) / (4 pi r)) to get the actual field value.
         -   `kx`, `ky`: Wavevector values corresponding to the x- and y- axes in E_far and H_far,
                 normalized to wavelength (dimensionless).
@@ -62,14 +95,10 @@ def near_to_farfield(
     if not all(s == f.shape for f in [*E_near, *H_near]):
         raise Exception('All fields must be the same shape!')
 
-    if padded_size is None:
+    padded_shape = _padded_shape(s, padded_size)
-        padded_size = (2**numpy.ceil(numpy.log2(s))).astype(int)
-    if not hasattr(padded_size, '__len__'):
-        padded_size = (padded_size, padded_size)            # type: ignore  # checked if sequence
-    padded_shape = cast('Sequence[int]', padded_size)
 
-    En_fft = [fftshift(fft2(fftshift(Eni), s=padded_shape)) for Eni in E_near]
+    En_fft = [fftshift(ifft2(ifftshift(_resize_centered(Eni, padded_shape)), norm='forward')) * dx * dy for Eni in E_near]
-    Hn_fft = [fftshift(fft2(fftshift(Hni), s=padded_shape)) for Hni in H_near]
+    Hn_fft = [fftshift(ifft2(ifftshift(_resize_centered(Hni, padded_shape)), norm='forward')) * dx * dy for Hni in H_near]
 
     # Propagation vectors kx, ky
     k  = 2 * pi
@@ -112,8 +141,8 @@ def near_to_farfield(
     outputs = {
         'E': E_far,
         'H': H_far,
-        'dkx': float(kxx[1] - kxx[0]),
+        'dkx': float(2 * pi / (padded_shape[0] * dx)),
-        'dky': float(kyy[1] - kyy[0]),
+        'dky': float(2 * pi / (padded_shape[1] * dy)),
         'kx': kx,
         'ky': ky,
         'theta': theta,
@@ -131,18 +160,17 @@ def far_to_nearfield(
         padded_size: Sequence[int] | int | None = None
         ) -> dict[str, Any]:
     """
-    Compute the farfield, i.e. the distribution of the fields after propagation
+    Reconstruct the propagating near-field spectrum from a normalized farfield.
-      through several wavelengths of uniform medium.
 
-    The input fields should be complex phasors.
+    Inputs are not modified. Spectral padding/cropping preserves zero frequency
+    at index ``size // 2``. The output spacings follow its resized shape, and
+    inverse Fourier quadrature restores the cell-area normalization.
 
     Args:
-        E_far: List of 2 ndarrays containing the 2D phasor field slices for the transverse
+        E_far: List of 2 ndarrays containing polar and azimuthal E components.
-                E fields (e.g. [Ex, Ey] for calculating the nearfield toward the z-direction).
                 Fields should be normalized so that
                 E_far = E_far_actual / (i k exp(-i k r) / (4 pi r))
-        H_far: List of 2 ndarrays containing the 2D phasor field slices for the transverse
+        H_far: List of 2 ndarrays containing polar and azimuthal H components.
-                H fields (e.g. [Hx, hy] for calculating the nearfield toward the z-direction).
                 Fields should be normalized so that
                 H_far = H_far_actual / (i k exp(-i k r) / (4 pi r))
         dkx: kx discretization, in units of wavelength.
@@ -168,11 +196,7 @@ def far_to_nearfield(
     if not all(s == f.shape for f in [*E_far, *H_far]):
         raise Exception('All fields must be the same shape!')
 
-    if padded_size is None:
+    padded_shape = _padded_shape(s, padded_size)
-        padded_size = (2 ** numpy.ceil(numpy.log2(s))).astype(int)
-    if not hasattr(padded_size, '__len__'):
-        padded_size = (padded_size, padded_size)            # type: ignore  # checked if sequence
-    padded_shape = cast('Sequence[int]', padded_size)
 
     k = 2 * pi
     kxs = dkx * fftshift(fftfreq(s[0], d=1 / s[0]))
@@ -197,9 +221,7 @@ def far_to_nearfield(
     invalid_ind = kxy2 >= k * k
     theta[invalid_ind] = 0
     phi[invalid_ind] = 0
-    for i in range(2):
+    E_far = [numpy.where(invalid_ind, 0, field) for field in E_far]
-        E_far[i][invalid_ind] = 0
-        H_far[i][invalid_ind] = 0
 
     # Normalized vector potentials N, L
     L = [0.5 * E_far[1],
@@ -237,11 +259,11 @@ def far_to_nearfield(
             ),
         ]
 
-    E_near = [ifftshift(ifft2(ifftshift(Ei), s=padded_shape)) for Ei in En_fft]
+    dx = 2 * pi / (padded_shape[0] * dkx)
-    H_near = [ifftshift(ifft2(ifftshift(Hi), s=padded_shape)) for Hi in Hn_fft]
+    dy = 2 * pi / (padded_shape[1] * dky)
 
-    dx = 2 * pi / (s[0] * dkx)
+    E_near = [fftshift(fft2(ifftshift(_resize_centered(Ei, padded_shape)), norm='forward')) / (dx * dy) for Ei in En_fft]
-    dy = 2 * pi / (s[1] * dky)
+    H_near = [fftshift(fft2(ifftshift(_resize_centered(Hi, padded_shape)), norm='forward')) / (dx * dy) for Hi in Hn_fft]
 
     outputs = {
         'E': E_near,

meanas/fdfd/functional.py

@@ -226,6 +226,8 @@ def poynting_e_cross_h(dxes: dx_lists_t) -> Callable[[cfdfield, cfdfield], cfdfi
         For time-average power, call it as `f(E, H.conj())` and take `Re(...) / 2`.
     """
     def exh(e: cfdfield, h: cfdfield) -> cfdfield_t:
+        dtype = numpy.result_type(e, h, *dxes[0], *dxes[1])
+        e, h = numpy.asarray(e, dtype=dtype), numpy.asarray(h, dtype=dtype)
         s = numpy.empty_like(e)
         ex = e[0] * dxes[0][0][:, None, None]
         ey = e[1] * dxes[0][1][None, :, None]

meanas/fdfd/operators.py

@@ -57,6 +57,8 @@ def e_full(
         (del x (1/mu * del x) - omega**2 * epsilon) E = -i * omega * J
 
     To make this matrix symmetric, use the preconditioners from `e_full_preconditioners()`.
+    PEC rows impose E = 0 with unit diagonal entries; the corresponding
+    right-hand-side entries must be zero (as done by solvers.generic).
 
     Args:
         omega: Angular frequency of the simulation
@@ -93,6 +95,8 @@ def e_full(
         m_div = sparse.diags_array(1 / mu)
 
     op = pe @ (ch @ pm @ m_div @ ce - omega**2 * e) @ pe
+    if pec is not None:
+        op += sparse.diags_array(numpy.where(pec, 1, 0))
     return op
 
 
@@ -137,6 +141,9 @@ def h_full(
 
         del x (1/epsilon * del x) - omega**2 * mu
 
+    PMC rows impose H = 0 with unit diagonal entries. Their right-hand-side
+    entries must be zero.
+
      for use with the H-field, with wave equation
      $$ (\nabla \times (\frac{1}{\epsilon} \nabla \times) - \omega^2 \mu) E = \imath \omega M $$
 
@@ -177,6 +184,8 @@ def h_full(
         m = sparse.diags_array(mu)
 
     A = pm @ (ce @ pe @ e_div @ ch - omega**2 * m) @ pm
+    if pmc is not None:
+        A += sparse.diags_array(numpy.where(pmc, 1, 0))
     return A
 
 
@@ -199,6 +208,9 @@ def eh_full(
         [[-i * omega * epsilon,  del x         ],
          [del x,                 i * omega * mu]]
 
+    PEC electric rows and PMC magnetic rows impose zero fields with unit
+    diagonal entries. The corresponding source entries must be zero.
+
     for use with a field vector of the form `cat(vec(E), vec(H))`:
     $$  \begin{bmatrix}
         -\imath \omega \epsilon  &  \nabla \times      \\
@@ -247,8 +259,10 @@ def eh_full(
     A1 = pe @ curl_back(dxes[1]) @ pm
     A2 = pm @ curl_forward(dxes[0]) @ pe
 
-    A = sparse.block_array([[-iwe, A1],
+    electric = -iwe if pec is None else -iwe + sparse.diags_array(numpy.where(pec, 1, 0))
-                            [A2,  iwm]])
+    magnetic = iwm if pmc is None else iwm + sparse.diags_array(numpy.where(pmc, 1, 0))
+    A = sparse.block_array([[electric, A1],
+                            [A2, magnetic]])
     return A
 
 

meanas/fdfd/solvers.py

@@ -32,7 +32,7 @@ def _scipy_qmr(
         kwargs: Passed as **kwargs to the wrapped function
 
     Returns:
-        Guess for solution (returned even if didn't converge)
+        Converged solution. Raises RuntimeError on breakdown or nonconvergence.
     """
 
     #
@@ -61,7 +61,10 @@ def _scipy_qmr(
     #
     # Run the actual solve
     #
-    x, _ = scipy.sparse.linalg.qmr(A, b, **kwargs)
+    x, info = scipy.sparse.linalg.qmr(A, b, **kwargs)
+    if info != 0:
+        residual = norm(A @ x - b)
+        raise RuntimeError(f'QMR failed to converge ({info=}, residual={residual:.6g})')
     return x
 
 
@@ -80,7 +83,7 @@ def generic(
         E_guess: vcfdfield | None = None,
         ) -> vcfdfield_t:
     """
-    Conjugate gradient FDFD solver using CSR sparse matrices.
+    Preconditioned FDFD solver using CSR sparse matrices.
 
     All ndarray arguments should be 1D arrays, as returned by `meanas.fdmath.vectorization.vec()`.
 
@@ -92,7 +95,8 @@ def generic(
         epsilon: Dielectric constant distribution (at E-field locations)
         mu: Magnetic permeability distribution (at H-field locations)
         pec: Perfect electric conductor distribution
-             (at E-field locations; non-zero value indicates PEC is present)
+             (at E-field locations; non-zero value indicates PEC is present).
+             Sources and initial guesses are zeroed at PEC locations.
         pmc: Perfect magnetic conductor distribution
              (at H-field locations; non-zero value indicates PMC is present)
         adjoint: If true, solves the adjoint problem.
@@ -101,8 +105,8 @@ def generic(
                       `b`: `ArrayLike`;
                       `x`: `ArrayLike`;
                 Default is a wrapped version of `scipy.sparse.linalg.qmr()`
-                 which doesn't return convergence info and logs the residual
+                 which raises on nonconvergence and logs the residual every
-                 every 100 iterations.
+                 100 iterations.
         matrix_solver_opts: Passed as kwargs to `matrix_solver(...)`
         E_guess: Guess at the solution E-field. `matrix_solver` must accept an
                 `x0` argument with the same purpose.
@@ -111,10 +115,11 @@ def generic(
         E-field which solves the system.
     """
 
-    if matrix_solver_opts is None:
+    matrix_solver_opts = {} if matrix_solver_opts is None else dict(matrix_solver_opts)
-        matrix_solver_opts = dict()
 
     b0 = -1j * omega * J
+    if pec is not None:
+        b0 = numpy.where(pec, 0, b0)
     A0 = operators.e_full(omega, dxes, epsilon=epsilon, mu=mu, pec=pec, pmc=pmc)
 
     Pl, Pr = operators.e_full_preconditioners(dxes)
@@ -127,17 +132,23 @@ def generic(
         b = Pl @ b0
 
     if E_guess is not None:
+        if pec is not None:
+            E_guess = numpy.where(pec, 0, E_guess)
         if adjoint:
             x0 = Pr.T.conjugate() @ E_guess
         else:
             x0 = Pl @ E_guess
         matrix_solver_opts['x0'] = x0
 
-    x = matrix_solver(A.tocsr(), b, **matrix_solver_opts)
+    x = numpy.asarray(matrix_solver(A.tocsr(), b, **matrix_solver_opts))
+    if x.shape != b.shape or not numpy.isfinite(x).all():
+        raise RuntimeError('matrix_solver must return a finite solution vector matching the source shape')
 
     if adjoint:
         x0 = Pl.T.conjugate() @ x
     else:
         x0 = Pr @ x
 
+    if pec is not None:
+        x0 = numpy.where(pec, 0, x0)
     return vcfdfield_t(x0)

meanas/fdfd/waveguide_2d.py

@@ -801,74 +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}
-    $$
-
-    However, $H_{yx}^\star$ is actually 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)
@@ -883,16 +849,17 @@ def sensitivity(
     mu_x, mu_y, _mu_z = numpy.split(mu, 3)
     mu_yx = sparse.diags_array(numpy.hstack((mu_y, mu_x)))
 
-    da_exxhyy = vec(dxes[1][0][:, None] * dxes[0][1][None, :])
+    da_exxhyy = vec(dxes[0][0][:, None] * dxes[1][1][None, :])
-    da_eyyhxx = vec(dxes[1][1][None, :] * dxes[0][0][:, None])
+    da_eyyhxx = vec(dxes[1][0][:, None] * dxes[0][1][None, :])
-    ev_xy = numpy.concatenate(numpy.split(e_norm, 3)[:2]) * numpy.concatenate([da_exxhyy, da_eyyhxx])
+    ev_xy = numpy.concatenate(numpy.split(e_norm, 3)[:2])
-    hx, hy, hz = numpy.split(h_norm, 3)
+    hx, hy, _hz = numpy.split(h_norm, 3)
-    hv_yx_conj = numpy.conj(numpy.concatenate([hy, -hx]))
+    # The Yee-weighted H projection is the left eigenvector on a graded grid.
+    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,15 +157,21 @@ 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)
     e = unvec(ve, shape)
     h = unvec(vh, shape)
 
-    # Adjust for propagation direction
+    # Reversing beta changes Ez and transverse H, while Hz is unchanged.
-    h *= polarity           # type: ignore          # mypy issue with numpy
+    e[2] *= polarity
+    h[:2] *= polarity
 
     # Apply phase shift to H-field
     h[:2] *= numpy.exp(-1j * polarity * 0.5 * wavenumber * dx_prop)
@@ -371,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)
@@ -391,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)
 
 
@@ -549,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,
@@ -636,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')
@@ -709,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] = []
@@ -746,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,
@@ -780,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/fdmath/functional.py

@@ -71,6 +71,12 @@ def deriv_back(
     return tuple(_derivative(axis, False, None if dx_h is None else dx_h[axis]) for axis in range(3))
 
 
+def _promote_field(field: NDArray, widths: Sequence[NDArray] | None) -> NDArray:
+    """Use one calculation dtype for all components of a curl."""
+    dtype = field.dtype if widths is None else numpy.result_type(field, *widths)
+    return numpy.asarray(field, dtype=dtype)
+
+
 def curl_forward(
         dx_e: Sequence[NDArray[floating | complexfloating]] | None = None,
         ) -> Callable[[NDArray[floating | complexfloating]], NDArray[floating | complexfloating]]:
@@ -89,8 +95,8 @@ def curl_forward(
     Dx, Dy, Dz = deriv_forward(dx_e)
 
     def ce_fun(e: NDArray[floating | complexfloating]) -> NDArray[floating | complexfloating]:
-        dtype = e.dtype if dx_e is None else numpy.result_type(e, *dx_e)
+        e = _promote_field(e, dx_e)
-        output = numpy.empty_like(e, dtype=dtype)
+        output = numpy.empty_like(e)
         output[0] = Dy(e[2])
         output[1] = Dz(e[0])
         output[2] = Dx(e[1])
@@ -120,8 +126,8 @@ def curl_back(
     Dx, Dy, Dz = deriv_back(dx_h)
 
     def ch_fun(h: NDArray[floating | complexfloating]) -> NDArray[floating | complexfloating]:
-        dtype = h.dtype if dx_h is None else numpy.result_type(h, *dx_h)
+        h = _promote_field(h, dx_h)
-        output = numpy.empty_like(h, dtype=dtype)
+        output = numpy.empty_like(h)
         output[0] = Dy(h[2])
         output[1] = Dz(h[0])
         output[2] = Dx(h[1])
@@ -139,6 +145,7 @@ def curl_forward_parts(
     Dx, Dy, Dz = deriv_forward(dx_e)
 
     def mkparts_fwd(e: fdfield) -> tuple[tuple[fdfield, fdfield], ...]:
+        e = _promote_field(e, dx_e)
         return ((-Dz(e[1]),  Dy(e[2])),
                 ( Dz(e[0]), -Dx(e[2])),
                 (-Dy(e[0]),  Dx(e[1])))
@@ -152,6 +159,7 @@ def curl_back_parts(
     Dx, Dy, Dz = deriv_back(dx_h)
 
     def mkparts_back(h: fdfield) -> tuple[tuple[fdfield, fdfield], ...]:
+        h = _promote_field(h, dx_h)
         return ((-Dz(h[1]),  Dy(h[2])),
                 ( Dz(h[0]), -Dx(h[2])),
                 (-Dy(h[0]),  Dx(h[1])))

meanas/fdmath/operators.py

@@ -120,7 +120,7 @@ def deriv_forward(
     def deriv(axis: int) -> sparse.sparray:
         return shift_circ(axis, shape, 1) - sparse.eye_array(n)
 
-    Ds = [sparse.diags_array(+1 / dx.ravel(order='C')) @ deriv(a)
+    Ds = [sparse.diags_array(numpy.reciprocal(dx.ravel(order='C'), dtype=numpy.result_type(dx, numpy.float64))) @ deriv(a)
           for a, dx in enumerate(dx_e_expanded)]
 
     return Ds
@@ -147,7 +147,7 @@ def deriv_back(
     def deriv(axis: int) -> sparse.sparray:
         return shift_circ(axis, shape, -1) - sparse.eye_array(n)
 
-    Ds = [sparse.diags_array(-1 / dx.ravel(order='C')) @ deriv(a)
+    Ds = [sparse.diags_array(-numpy.reciprocal(dx.ravel(order='C'), dtype=numpy.result_type(dx, numpy.float64))) @ deriv(a)
           for a, dx in enumerate(dx_h_expanded)]
 
     return Ds

meanas/test/test_bloch_foundations.py

@@ -1,10 +1,61 @@
 import numpy
+import pytest
 from numpy.linalg import norm
 
 from ..fdfd import bloch
 from ._bloch_case import EPSILON, G_MATRIX, H_MN, K0_AXIAL, K0_GENERAL, MU, SHAPE, ZERO_H_MN
 from .utils import assert_close
 
+
+@pytest.mark.parametrize('k0', [[0.2, 0, 0], [0, 0.2, 0], [0, 0, 0.2], [0, 0.2, 0.1], [0, 0, 0]])
+def test_kmn_basis_is_transverse_on_axes_and_at_gamma(k0) -> None:
+    k0 = numpy.asarray(k0)
+    shape = (2, 3, 2)
+    _magnitude, m, n = bloch.generate_kmn(k0, numpy.eye(3), shape)
+    reciprocal_indices = numpy.moveaxis(numpy.array(numpy.meshgrid(
+        *(numpy.fft.fftfreq(size, 1 / size) for size in shape), indexing='ij',
+    )), 0, -1)
+    k = k0 - reciprocal_indices
+    assert_close(numpy.sum(k * m, axis=-1), 0, atol=1e-14)
+    assert_close(numpy.sum(k * n, axis=-1), 0, atol=1e-14)
+    assert_close(numpy.sum(m * n, axis=-1), 0, atol=1e-14)
+    assert_close(norm(m, axis=-1), 1)
+    assert_close(norm(n, axis=-1), 1)
+
+
+@pytest.mark.parametrize('mu', [None, MU])
+def test_inverse_preconditioner_is_finite_at_gamma(mu) -> None:
+    operator = bloch.inverse_maxwell_operator_approx(numpy.zeros(3), G_MATRIX, EPSILON, mu)
+    with numpy.errstate(divide='raise', invalid='raise', over='raise'):
+        result = operator(H_MN)
+        zero = operator(ZERO_H_MN)
+    assert numpy.isfinite(result).all()
+    assert_close(zero, 0)
+
+
+def test_anisotropic_magnetic_bloch_modes_obey_full_fourier_maxwell_equations() -> None:
+    shape = EPSILON.shape[1:]
+    size = 2 * numpy.prod(shape)
+    operator = bloch.maxwell_operator(K0_GENERAL, G_MATRIX, EPSILON, MU)
+    matrix = numpy.column_stack([operator(column).ravel() for column in numpy.eye(size, dtype=complex)])
+    values, vectors = numpy.linalg.eig(matrix)
+    reciprocal_indices = numpy.moveaxis(numpy.array(numpy.meshgrid(
+        *(numpy.fft.fftfreq(n, 1 / n) for n in shape), indexing='ij',
+    )), 0, -1)
+    k = K0_GENERAL - reciprocal_indices
+    epsilon, mu = numpy.moveaxis(EPSILON, 0, -1), numpy.moveaxis(MU, 0, -1)
+    h_converter = bloch.hmn_2_hxyz(K0_GENERAL, G_MATRIX, EPSILON, MU)
+    e_converter = bloch.hmn_2_exyz(K0_GENERAL, G_MATRIX, EPSILON, MU)
+    for value, vector in zip(values, vectors.T, strict=True):
+        h = numpy.moveaxis(h_converter(vector), 0, -1)
+        b = mu * h
+        h_fourier = numpy.fft.fftn(h, axes=(0, 1, 2))
+        electric = numpy.fft.ifftn(numpy.cross(k, h_fourier), axes=(0, 1, 2)) / epsilon
+        curl_electric = numpy.fft.ifftn(numpy.cross(k, numpy.fft.fftn(electric, axes=(0, 1, 2))), axes=(0, 1, 2))
+        assert_close(curl_electric, value * b, atol=1e-13, rtol=1e-11)
+        assert_close(numpy.sum(k * numpy.fft.fftn(b, axes=(0, 1, 2)), axis=-1), 0, atol=1e-13)
+        assert_close(numpy.moveaxis(e_converter(vector), 0, -1), electric, atol=1e-13)
+
 def test_generate_kmn_general_case_returns_orthonormal_basis() -> None:
     k_mag, m_vecs, n_vecs = bloch.generate_kmn(K0_GENERAL, G_MATRIX, SHAPE)
 

meanas/test/test_bloch_interactions.py

@@ -4,7 +4,7 @@ from numpy.testing import assert_allclose
 from types import SimpleNamespace
 
 from ..fdfd import bloch
-from ._bloch_case import EPSILON, G_MATRIX, H_SIZE, K0_X, Y0, Y0_TWO_MODE, build_overlap_fixture
+from ._bloch_case import EPSILON, G_MATRIX, H_SIZE, K0_X, K0_GENERAL, MU, Y0, Y0_TWO_MODE, build_overlap_fixture
 from .utils import assert_close
 
 
@@ -81,6 +81,42 @@ def test_eigsolve_returns_finite_modes_with_small_residual() -> None:
     assert callback_count > 0
 
 
+@pytest.mark.parametrize('count', [1, 2])
+def test_iteration_limit_returns_consistent_ritz_pairs_without_mutating_guess(count) -> None:
+    rng = numpy.random.default_rng(11)
+    guess = rng.normal(size=(count, H_SIZE)) + 1j * rng.normal(size=(count, H_SIZE))
+    original = guess.copy()
+    values, vectors = bloch.eigsolve(count, K0_X, G_MATRIX, EPSILON, y0=guess, max_iters=1)
+    operator = bloch.maxwell_operator(K0_X, G_MATRIX, EPSILON)
+    for value, vector in zip(values, vectors, strict=True):
+        rayleigh = numpy.vdot(vector, operator(vector).ravel()) / numpy.vdot(vector, vector)
+        assert_close(value, rayleigh, atol=1e-13)
+    assert_allclose(vectors.conj() @ vectors.T, numpy.eye(count), atol=1e-12)
+    numpy.testing.assert_array_equal(guess, original)
+
+
+@pytest.mark.parametrize('mu', [None, MU])
+def test_eigsolve_at_gamma_returns_a_static_mode_with_small_residual(mu) -> None:
+    with numpy.errstate(divide='raise', invalid='raise', over='raise'):
+        values, vectors = bloch.eigsolve(1, numpy.zeros(3), G_MATRIX, EPSILON, mu, y0=Y0, max_iters=200, tolerance=1e-9)
+    residual = bloch.maxwell_operator(numpy.zeros(3), G_MATRIX, EPSILON, mu)(vectors[0]).ravel() - values[0] * vectors[0]
+    assert abs(values[0]) < 1e-12
+    assert numpy.linalg.norm(residual) < 1e-9
+
+
+def test_anisotropic_magnetic_eigensolve_matches_dense_reference() -> None:
+    operator = bloch.maxwell_operator(K0_GENERAL, G_MATRIX, EPSILON, MU)
+    matrix = numpy.column_stack([operator(column).ravel() for column in numpy.eye(H_SIZE, dtype=complex)])
+    reference = numpy.linalg.eigvals(matrix)
+    reference = reference[numpy.argsort(numpy.abs(reference))[:3]]
+    rng = numpy.random.default_rng(9)
+    guess = rng.normal(size=(3, H_SIZE)) + 1j * rng.normal(size=(3, H_SIZE))
+    values, vectors = bloch.eigsolve(3, K0_GENERAL, G_MATRIX, EPSILON, MU, y0=guess, max_iters=300, tolerance=1e-10)
+    assert_allclose(values, reference, atol=1e-11, rtol=1e-9)
+    for value, vector in zip(values, vectors, strict=True):
+        assert numpy.linalg.norm(matrix @ vector - value * vector) < 1e-9
+
+
 def test_eigsolve_without_initial_guess_returns_finite_modes() -> None:
     eigvals, eigvecs = bloch.eigsolve(
         1,

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_fdfd.py

@@ -3,6 +3,7 @@ import dataclasses
 import pytest       # type: ignore
 import numpy
 from numpy.typing import NDArray
+from scipy.sparse.linalg import spsolve
 #from numpy.testing import assert_allclose, assert_array_equal
 
 from .. import fdfd
@@ -145,7 +146,9 @@ def sim(
         omega = omega,
         dxes = dxes,
         epsilon = eps_vec,
-        matrix_solver_opts = dict(atol=1e-15, rtol=1e-11),
+        # These small, nearly static systems are ill-conditioned. Test the
+        # Maxwell residual and flux identities independently of QMR convergence.
+        matrix_solver = spsolve,
         )
     e = unvec(e_vec, shape[1:])
 

meanas/test/test_fdfd_algebra_helpers.py

@@ -34,7 +34,7 @@ def _dense_e_full(mu: numpy.ndarray | None) -> numpy.ndarray:
         curl_term = curl_term / mu
     curl_term = pm * curl_term
     curl_term = ch(curl_term)
-    return pe * (curl_term - OMEGA**2 * EPSILON * masked_e)
+    return pe * (curl_term - OMEGA**2 * EPSILON * masked_e) + numpy.where(PEC, E_FIELD, 0)
 
 
 def _dense_h_full(mu: numpy.ndarray | None) -> numpy.ndarray:
@@ -48,7 +48,7 @@ def _dense_h_full(mu: numpy.ndarray | None) -> numpy.ndarray:
     curl_term = ch(masked_h)
     curl_term = pe * (curl_term / EPSILON)
     curl_term = ce(curl_term)
-    return pm * (curl_term - OMEGA**2 * magnetic * masked_h)
+    return pm * (curl_term - OMEGA**2 * magnetic * masked_h) + numpy.where(PMC, H_FIELD, 0)
 
 
 def test_h_full_matches_dense_reference_with_and_without_mu() -> None:
@@ -95,6 +95,8 @@ def test_eh_full_matches_manual_block_operator_with_masks() -> None:
 
     dense_e = pe * ch(pm * H_FIELD) - pe * (1j * OMEGA * EPSILON * (pe * E_FIELD))
     dense_h = pm * ce(pe * E_FIELD) + pm * (1j * OMEGA * MU * (pm * H_FIELD))
+    dense_e += numpy.where(PEC, E_FIELD, 0)
+    dense_h += numpy.where(PMC, H_FIELD, 0)
 
     assert_fields_close(matrix_e, dense_e, atol=1e-10, rtol=1e-10)
     assert_fields_close(matrix_h, dense_h, atol=1e-10, rtol=1e-10)

meanas/test/test_fdfd_farfield.py

@@ -98,3 +98,76 @@ def test_farfield_roundtrip_supports_rectangular_arrays() -> None:
     assert restored['dy'] == pytest.approx(0.3)
     assert numpy.isfinite(restored['E'][0]).all()
     assert numpy.isfinite(restored['H'][0]).all()
+
+
+@pytest.mark.parametrize('shape', [(4, 4), (4, 5), (5, 4), (5, 7)])
+def test_farfield_padding_preserves_aperture_origin(shape) -> None:
+    rng = numpy.random.default_rng(7)
+    fields = [rng.normal(size=shape) + 1j * rng.normal(size=shape) for _ in range(4)]
+    target = (8, 10)
+    padding = [(new // 2 - old // 2, new - old - (new // 2 - old // 2)) for old, new in zip(shape, target, strict=True)]
+    padded = [numpy.pad(field, padding) for field in fields]
+    automatic = farfield.near_to_farfield(fields[:2], fields[2:], 0.5, 0.7, padded_size=target)
+    manual = farfield.near_to_farfield(padded[:2], padded[2:], 0.5, 0.7, padded_size=target)
+    for key in ('E', 'H'):
+        numpy.testing.assert_allclose(automatic[key], manual[key], atol=1e-13)
+
+
+@pytest.mark.parametrize('size', [4, 8])
+def test_farfield_amplitude_is_independent_of_aperture_sampling(size) -> None:
+    # A constant 2 by 2 wavelength aperture has integral 4 on either mesh.
+    ones = numpy.ones((size, size), dtype=complex)
+    zeros = numpy.zeros_like(ones)
+    result = farfield.near_to_farfield([ones, zeros], [zeros, ones], 2 / size, 2 / size, padded_size=size)
+    assert result['E'][0][size // 2, size // 2] == pytest.approx(8)
+
+
+@pytest.mark.parametrize('shape', [(4, 4), (16, 12), (7, 9)])
+def test_inverse_resizing_matches_plane_wave_quadrature(shape) -> None:
+    spectrum_shape = (8, 6)
+    dkx, dky = 0.4, 0.3
+    amplitude = 1.7 + 0.2j
+    e = [numpy.zeros(spectrum_shape, dtype=complex) for _ in range(2)]
+    e[0][spectrum_shape[0] // 2 + 1, spectrum_shape[1] // 2] = 2 * amplitude
+    h = [-e[1], e[0]]
+    result = farfield.far_to_nearfield(e, h, dkx, dky, padded_size=shape)
+    dx, dy = 2 * numpy.pi / (shape[0] * dkx), 2 * numpy.pi / (shape[1] * dky)
+    assert result['dx'] == pytest.approx(dx)
+    assert result['dy'] == pytest.approx(dy)
+    x = (numpy.arange(shape[0]) - shape[0] // 2) * dx
+    expected = numpy.broadcast_to(amplitude * dkx * dky / (2 * numpy.pi) ** 2 * numpy.exp(-1j * dkx * x[:, None]), shape)
+    numpy.testing.assert_allclose(result['E'][0], expected, atol=1e-14)
+    numpy.testing.assert_allclose(result['H'][1], expected / numpy.sqrt(1 - (dkx / (2 * numpy.pi)) ** 2), atol=1e-14)
+
+
+def test_inverse_accepts_readonly_spectra_without_modifying_them() -> None:
+    rng = numpy.random.default_rng(9)
+    fields = [rng.normal(size=(8, 8)).astype(complex) for _ in range(4)]
+    originals = [field.copy() for field in fields]
+    for field in fields:
+        field.flags.writeable = False
+    result = farfield.far_to_nearfield(fields[:2], fields[2:], dkx=3, dky=3)
+    assert numpy.isfinite(result['E']).all()
+    for field, original in zip(fields, originals, strict=True):
+        numpy.testing.assert_array_equal(field, original)
+
+
+@pytest.mark.parametrize('direction', [-1, 1])
+def test_farfield_direction_and_roundtrip_match_outgoing_phasors(direction) -> None:
+    shape = (16, 8)
+    dx = dy = 0.25
+    k = 2 * numpy.pi
+    kx = direction * 2 * numpy.pi / (shape[0] * dx)
+    kz = numpy.sqrt(k * k - kx * kx)
+    x = (numpy.arange(shape[0]) - shape[0] // 2) * dx
+    phase = numpy.broadcast_to(numpy.exp(-1j * kx * x[:, None]), shape)
+    zeros = numpy.zeros_like(phase)
+    electric, magnetic = [kz / k * phase, zeros], [zeros, phase]
+    result = farfield.near_to_farfield(electric, magnetic, dx, dy, padded_size=shape)
+    intensity = numpy.abs(result['E'][0]) ** 2 + numpy.abs(result['E'][1]) ** 2
+    peak = numpy.unravel_index(numpy.argmax(intensity), shape)
+    assert result['kx'][peak] == pytest.approx(kx)
+    assert result['phi'][peak] == pytest.approx(numpy.arccos(kz / k))
+    restored = farfield.far_to_nearfield(result['E'], result['H'], result['dkx'], result['dky'], padded_size=shape)
+    numpy.testing.assert_allclose(restored['E'], electric, atol=1e-13)
+    numpy.testing.assert_allclose(restored['H'], magnetic, atol=1e-13)

meanas/test/test_fdfd_functional.py

@@ -1,4 +1,5 @@
 import numpy
+import pytest
 
 from ..fdmath import unvec, vec
 from ..fdfd import functional, operators
@@ -120,3 +121,17 @@ def test_poynting_e_cross_h_matches_sparse_operator() -> None:
     functional_result = functional.poynting_e_cross_h(DXES)(E_FIELD, H_FIELD)
 
     assert_fields_match(functional_result, matrix_result)
+
+
+@pytest.mark.parametrize('wide_input', ['electric', 'magnetic', 'widths'])
+def test_poynting_uses_promoted_precision_before_products(wide_input) -> None:
+    shape = (2, 2, 2)
+    e = numpy.zeros((3, *shape), dtype=numpy.complex128 if wide_input == 'electric' else numpy.complex64)
+    h = numpy.zeros((3, *shape), dtype=numpy.complex128 if wide_input == 'magnetic' else numpy.complex64)
+    e[1], h[2] = 1e20, 1e20j
+    widths = [[numpy.ones(n, dtype=numpy.float64 if wide_input == 'widths' else numpy.float32) for n in shape] for _ in range(2)]
+    with numpy.errstate(over='raise', invalid='raise'):
+        result = functional.poynting_e_cross_h(widths)(e, h)
+    expected = unvec(operators.poynting_e_cross(vec(e.astype(complex)), widths) @ vec(h.astype(complex)), shape)
+    assert result.dtype == numpy.complex128
+    assert_fields_close(result, expected, atol=0, rtol=1e-14)

meanas/test/test_fdfd_solvers.py

@@ -1,8 +1,10 @@
 from typing import cast
 
 import numpy
+import pytest
+from scipy.sparse.linalg import spsolve
 
-from ..fdfd import solvers
+from ..fdfd import solvers, operators
 from ._solver_cases import solver_plumbing_case
 from .utils import assert_close
 
@@ -126,3 +128,41 @@ def test_generic_without_guess_does_not_inject_x0(monkeypatch) -> None:
 
     assert 'x0' not in cast(dict[str, object], captured['kwargs'])
     assert_close(result, case.pr @ numpy.array([1.0, -1.0]))
+
+
+def test_default_solver_reports_breakdown_in_nonsingular_maxwell_system() -> None:
+    dxes = [[numpy.ones(2), numpy.ones(1), numpy.ones(1)] for _ in range(2)]
+    epsilon = numpy.full(6, 2.0)
+    current = numpy.zeros(6, dtype=complex)
+    current[2] = 1
+    with pytest.raises(RuntimeError, match='QMR failed to converge.*info=-14'):
+        solvers.generic(1.0, dxes, current, epsilon)
+    direct = solvers.generic(1.0, dxes, current, epsilon, matrix_solver=spsolve)
+    assert_close(operators.e_full(1.0, dxes, epsilon) @ direct, -1j * current)
+
+
+def test_qmr_iteration_limit_is_not_silently_accepted(monkeypatch) -> None:
+    monkeypatch.setattr(solvers.scipy.sparse.linalg, 'qmr', lambda *args, **kwargs: (numpy.zeros(2), 1))
+    with pytest.raises(RuntimeError, match='info=1'):
+        solvers._scipy_qmr(solver_plumbing_case().a0, numpy.ones(2))
+
+
+@pytest.mark.parametrize('adjoint', [False, True])
+@pytest.mark.parametrize('direct', [False, True])
+def test_pec_constraints_remove_sources_and_initial_fields(adjoint, direct) -> None:
+    dxes = [[numpy.ones(2), numpy.ones(1), numpy.ones(1)] for _ in range(2)]
+    epsilon = numpy.full(6, 2.0)
+    current = numpy.array([7, 0, 1, 0, 0, 0], dtype=complex)
+    pec = numpy.array([1, 0, 0, 0, 0, 0])
+    guess = numpy.ones(6, dtype=complex)
+    options = {'rtol': 1e-10}
+    kwargs = {'matrix_solver': spsolve} if direct else {'E_guess': guess, 'matrix_solver_opts': options}
+    field = solvers.generic(0.7, dxes, current, epsilon, pec=pec, adjoint=adjoint, **kwargs)
+    assert field[0] == 0
+    assert numpy.isfinite(field).all()
+    matrix = operators.e_full(0.7, dxes, epsilon, pec=pec)
+    if adjoint:
+        matrix = matrix.conj().T
+    assert_close(matrix @ field, -0.7j * numpy.where(pec, 0, current), atol=1e-10)
+    assert 'x0' not in options
+    assert guess[0] == 1

meanas/test/test_fdmath_functional.py

@@ -87,6 +87,24 @@ def test_curls_preserve_small_residual_between_large_derivative_terms(name) -> N
     assert_fields_close(result, expected, atol=1e-7, rtol=1e-12)
 
 
+@pytest.mark.parametrize('name', ['curl_forward', 'curl_back'])
+@pytest.mark.parametrize('complex_field', [False, True])
+def test_curls_use_common_precision_for_mixed_width_dtypes(name, complex_field) -> None:
+    shape = (1, 2, 2)
+    field = numpy.zeros((3, *shape), dtype=numpy.complex64 if complex_field else numpy.float32)
+    field[2, :, 1, :] = 1e8
+    field[1, :, :, 1] = 1e8
+    widths = [numpy.ones(1), numpy.full(2, 0.7, dtype=numpy.float32), numpy.full(2, float(numpy.float32(0.7)))]
+    wide_widths = [dx.astype(numpy.float64) for dx in widths]
+    expected = unvec(getattr(fd_operators, name)(wide_widths) @ vec(field), shape)
+    assert expected[0, 0, 0, 0] == 0
+    result = getattr(fd_functional, name)(widths)(field)
+    assert_fields_close(result, expected, atol=1e-8, rtol=1e-12)
+    parts = getattr(fd_functional, name + '_parts')(widths)(field)
+    for component, (left, right) in enumerate(parts):
+        assert_close(left + right, expected[component], atol=1e-8, rtol=1e-12)
+
+
 @pytest.mark.parametrize('name', ['deriv_forward', 'deriv_back'])
 @pytest.mark.parametrize('axis', [0, 1, 2])
 @pytest.mark.parametrize('spaced', [False, True])

meanas/test/test_fdmath_operators.py

@@ -91,6 +91,22 @@ def test_vec_cross_matches_pointwise_cross_product() -> None:
     assert_close(matrix_result, expected)
 
 
+@pytest.mark.parametrize('name', ['deriv_forward', 'deriv_back'])
+@pytest.mark.parametrize('complex_widths', [False, True])
+def test_sparse_derivatives_widen_reciprocals_before_overflow(name, complex_widths) -> None:
+    dtype = numpy.complex64 if complex_widths else numpy.float32
+    width = numpy.array(1e-39 * (1 + 1j) if complex_widths else 1e-39, dtype=dtype)
+    widths = [numpy.full(2, width, dtype=dtype), numpy.ones(1), numpy.ones(1)]
+    wide_widths = [dx.astype(numpy.complex128 if complex_widths else numpy.float64) for dx in widths]
+    field = numpy.array([0, 1e-39])
+    with numpy.errstate(over='raise', invalid='raise'):
+        matrix = getattr(operators, name)(widths)[0]
+        result = matrix @ field
+    assert numpy.isfinite(matrix.data).all()
+    expected = getattr(operators, name)(wide_widths)[0] @ field
+    assert_close(result, expected, atol=1e-12, rtol=1e-12)
+
+
 def test_avg_forward_matches_half_sum_with_forward_neighbor() -> None:
     matrix_result = _apply_scalar_matrix(operators.avg_forward(1, SHAPE))
     expected = 0.5 * (SCALAR_FIELD + numpy.roll(SCALAR_FIELD, -1, axis=1))

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,9 +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
@@ -120,6 +122,38 @@ def test_waveguide_2d_sensitivity_matches_finite_difference() -> None:
     assert 0.4 < ratio < 1.8
 
 
+@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
+    epsilon = vec(1 + rng.uniform(0, 3, (3, *shape)))
+    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])
+    e, h = waveguide_2d.normalized_fields_e(vectors[:, selected], beta, omega, dxes, epsilon)
+    sensitivity = waveguide_2d.sensitivity(e, h, beta, omega, dxes, epsilon)
+    delta = 1e-5
+    finite_difference = numpy.empty(epsilon.size, dtype=complex)
+    for index in range(epsilon.size):
+        roots = []
+        for direction in (-1, 1):
+            perturbed = epsilon.copy()
+            perturbed[index] += direction * delta
+            roots_squared = linalg.eigvals(waveguide_2d.operator_e(omega, dxes, perturbed).toarray())
+            candidates = numpy.sqrt(roots_squared)
+            roots.append(candidates[numpy.argmin(numpy.abs(candidates - beta))])
+        finite_difference[index] = (roots[1] - roots[0]) / (2 * delta)
+    assert_allclose(sensitivity, finite_difference, atol=1e-8, rtol=1e-5)
+
+
 def test_waveguide_2d_normalized_fields_h_are_finite_and_unit_normalized_with_mu() -> None:
     epsilon = build_asymmetric_epsilon()
     mu = build_mu_profile()

meanas/test/test_waveguide_fdtd_fdfd.py

@@ -36,7 +36,7 @@ SCATTERING_SOURCE_SLICES = (slice(4, 5), slice(None), slice(None))
 SCATTERING_REFLECT_SLICES = (slice(10, 11), slice(None), slice(None))
 SCATTERING_TRANSMIT_SLICES = (slice(29, 30), slice(None), slice(None))
 SCATTERING_STEP_X = 18
-SCATTERING_WARMUP_PERIODS = 10
+SCATTERING_WARMUP_PERIODS = 30
 SCATTERING_ACCUMULATION_PERIODS = 10
 
 
@@ -221,6 +221,21 @@ def _build_cpml_params() -> list[list[dict[str, numpy.ndarray | float]]]:
     ]
 
 
+def _build_cpml_frequency_dxes(base_dxes: list[list[numpy.ndarray]]) -> list[list[numpy.ndarray]]:
+    """Match the discrete CPML recurrence, including its Yee sample offsets."""
+    widths = [[dx.astype(complex) for dx in family] for family in base_dxes]
+    for axis, faces in enumerate(_build_cpml_params()):
+        for face in faces:
+            # psi[n] = p0 psi[n-1] + p1 p2 D[n]. Its phasor response adds
+            # p1 / (1 - p0 exp(-i omega dt)) to the scaled derivative.
+            # Forward electric derivatives use the H-update CPML parameters.
+            for family, parameter in enumerate(('param_h', 'param_e')):
+                p0, p1, p2 = face[parameter]
+                response = p2 * (1 + p1 / (1 - p0 * numpy.exp(-1j * OMEGA * DT)))
+                widths[family][axis][face['region'][axis]] /= response.ravel()
+    return widths
+
+
 def _build_complex_pulse_waveform(total_steps: int) -> tuple[numpy.ndarray, complex]:
     source_phasor, _delay = gaussian_packet(wl=WAVELENGTH, dwl=PULSE_DWL, dt=DT, turn_on=1e-5)
     aa, cc, ss = source_phasor(numpy.arange(total_steps) + 0.5)
@@ -450,7 +465,10 @@ def _run_straight_waveguide_case(variant: str) -> WaveguideCalibrationResult:
 def _run_width_step_scattering_case() -> WaveguideScatteringResult:
     epsilon = _build_scattering_epsilon()
     base_dxes = _build_uniform_dxes(SCATTERING_SHAPE)
-    stretched_dxes = _build_stretched_dxes(base_dxes)
+    # A correctly signed backward mode resolves the weak reflected wave;
+    # incident leakage no longer hides PML differences or startup transients.
+    stretched_dxes = _build_cpml_frequency_dxes(base_dxes)
+    discrete_frequency = 2 * numpy.sin(OMEGA * DT / 2) / DT
 
     source_mode = waveguide_3d.solve_mode(
         0,
@@ -539,14 +557,14 @@ def _run_width_step_scattering_case() -> WaveguideScatteringResult:
     e_fdfd = unvec(
         fdfd.solvers.generic(
             J=vec(j_ph),
-            omega=OMEGA,
+            omega=discrete_frequency,
             dxes=stretched_dxes,
             epsilon=vec(epsilon),
             matrix_solver_opts={'atol': 1e-10, 'rtol': 1e-7},
         ),
         SCATTERING_SHAPE[1:],
     )
-    h_fdfd = functional.e2h(OMEGA, stretched_dxes)(e_fdfd)
+    h_fdfd = functional.e2h(discrete_frequency, stretched_dxes)(e_fdfd)
 
     reflected_td = complex(vec(e_ph) @ vec(reflected_overlap).conj())
     reflected_fd = complex(vec(e_fdfd) @ vec(reflected_overlap).conj())

meanas/test/test_waveguide_mode_helpers.py

@@ -2,17 +2,49 @@ import contextlib
 import io
 import numpy
 from numpy.linalg import norm
+from scipy import linalg
 import pytest
 import warnings
 
 from ..fdmath import vec, unvec
 from ..fdmath.operators import deriv_forward, deriv_back
+from ..fdmath.functional import curl_forward, curl_back
 from ..fdfd import waveguide_2d, waveguide_3d, waveguide_cyl
 
 
 OMEGA = 1 / 1500
 
 
+@pytest.mark.parametrize('axis', [0, 1, 2])
+@pytest.mark.parametrize('polarity', [-1, 1])
+def test_embedded_hybrid_mode_obeys_both_maxwell_curls(axis, polarity) -> None:
+    omega = 0.7
+    # Cyclic rotation preserves the orientation of the vector cross product.
+    permutation = numpy.roll(numpy.arange(3), axis + 1)
+    rng = numpy.random.default_rng(4)
+    canonical = numpy.repeat((1 + rng.uniform(0, 3, (3, 4, 3)))[..., None], 9, axis=3)
+    epsilon = canonical[permutation].transpose(0, *(permutation + 1))
+    shape = epsilon.shape[1:]
+    dxes = [[numpy.full(n, 0.4) for n in shape] for _ in range(2)]
+    cross_dxes = [[numpy.full(4, 0.4), numpy.full(3, 0.4)] for _ in range(2)]
+    values, vectors = linalg.eig(waveguide_2d.operator_e(omega, cross_dxes, vec(canonical[:, :, :, 0])).toarray())
+    selected = numpy.argmax(values.real)
+    slices = [slice(None)] * 3
+    slices[axis] = slice(4, 5)
+    mode = waveguide_3d.mode_from_2d(
+        vectors[:, selected], numpy.sqrt(values[selected]), omega=omega, dxes=dxes,
+        axis=axis, polarity=polarity, slices=slices, epsilon=epsilon,
+    )
+    e, h = [waveguide_3d.expand_e(mode[field], mode['wavenumber'], dxes, axis, polarity, slices) for field in ('E', 'H')]
+    assert norm(e[axis]) > 0.01 * norm(e)
+    assert norm(h[axis]) > 0.01 * norm(h)
+    region = [slice(None)] * 4
+    region[axis + 1] = slice(1, -1)
+    ce, ch = curl_forward(dxes[0])(e), curl_back(dxes[1])(h)
+    assert relative_residual(ce[tuple(region)], (-1j * omega * h)[tuple(region)]) < 1e-11
+    assert relative_residual(ch[tuple(region)], (1j * omega * epsilon * e)[tuple(region)]) < 1e-11
+
+
 def build_waveguide_3d_mode(
         *,
         slice_start: int,
@@ -121,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)
 
@@ -390,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
@@ -437,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: