[fdfd] report solver failures and enforce conductor constraints

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],
-                            [A2,  iwm]])
+    electric = -iwe if pec is None else -iwe + sparse.diags_array(numpy.where(pec, 1, 0))
+    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
-                 every 100 iterations.
+                 which raises on nonconvergence and logs the residual every
+                 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 = dict()
+    matrix_solver_opts = {} if matrix_solver_opts is None else dict(matrix_solver_opts)
 
     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/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_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