[bloch] solve transverse magnetic flux with a consistent energy metric

docs/api/fdfd.md

@@ -14,6 +14,25 @@
 
 ::: 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.
+
 ## 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/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,