From cde834582d4b5f40c36d7efa393e641144fc2017 Mon Sep 17 00:00:00 2001 From: Forgejo Actions Date: Wed, 28 Oct 2026 19:03:24 -0700 Subject: [PATCH] [fdfd] report solver failures and enforce conductor constraints --- meanas/fdfd/operators.py | 18 ++++++++-- meanas/fdfd/solvers.py | 29 +++++++++++----- meanas/test/test_fdfd.py | 5 ++- meanas/test/test_fdfd_algebra_helpers.py | 6 ++-- meanas/test/test_fdfd_solvers.py | 42 +++++++++++++++++++++++- 5 files changed, 85 insertions(+), 15 deletions(-) diff --git a/meanas/fdfd/operators.py b/meanas/fdfd/operators.py index 6425a29..6ff661f 100644 --- a/meanas/fdfd/operators.py +++ b/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 diff --git a/meanas/fdfd/solvers.py b/meanas/fdfd/solvers.py index e1f394c..55bef79 100644 --- a/meanas/fdfd/solvers.py +++ b/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) diff --git a/meanas/test/test_fdfd.py b/meanas/test/test_fdfd.py index 82587b4..9db4f05 100644 --- a/meanas/test/test_fdfd.py +++ b/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:]) diff --git a/meanas/test/test_fdfd_algebra_helpers.py b/meanas/test/test_fdfd_algebra_helpers.py index bf6e973..e3370dc 100644 --- a/meanas/test/test_fdfd_algebra_helpers.py +++ b/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) diff --git a/meanas/test/test_fdfd_solvers.py b/meanas/test/test_fdfd_solvers.py index de39d70..c7e60eb 100644 --- a/meanas/test/test_fdfd_solvers.py +++ b/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