From 007ebce5022493f6b25130101c39995a24c7458e Mon Sep 17 00:00:00 2001 From: Stefan Possanner Date: Fri, 23 Oct 2026 08:13:40 +0200 Subject: [PATCH 1/6] Document the discretization of MaxwellWeakAmpere and VlasovAmpereCoupling Semi-discrete systems in struphy operator notation (M1, M2, C, Lambda^1, W), Crank-Nicolson update, SchurSolver blocks A, B, C, resulting Schur system, accumulated matrix/vector and conserved quantities. Co-Authored-By: Claude Opus 5.5 --- .../propagators/maxwell_weak_ampere.py | 45 ++++++++++++++-- .../propagators/vlasov_ampere_coupling.py | 53 ++++++++++++++++++- 2 files changed, 92 insertions(+), 6 deletions(-) diff --git a/src/struphy/propagators/maxwell_weak_ampere.py b/src/struphy/propagators/maxwell_weak_ampere.py index 9c5d45486..777753235 100644 --- a/src/struphy/propagators/maxwell_weak_ampere.py +++ b/src/struphy/propagators/maxwell_weak_ampere.py @@ -30,12 +30,49 @@ class MaxwellWeakAmpere(Propagator): \\[2mm] \frac{\partial \mathbf B}{\partial t} + \nabla\times\mathbf E &= 0\,. - Time discretization: + Ampère's law holds weakly, Faraday's law strongly. With the FE coefficients + :math:`\mathbf e \in \mathbb R^{N_1}` of :math:`\mathbf E \in V^1_h` and :math:`\mathbf b \in \mathbb R^{N_2}` + of :math:`\mathbf B \in V^2_h`, the semi-discrete system reads - - implicit: Crank-Nicolson (implicit mid-point) - - explicit: explicit RK methods from ButcherTableau + .. math:: + + \mathbb M_1 \dot{\mathbf e} = \mathbb C^\top \mathbb M_2 \mathbf b\,, \qquad + \dot{\mathbf b} = - \mathbb C \mathbf e\,, + + with the mass matrices :math:`\mathbb M_1`, :math:`\mathbb M_2` + (:class:`~struphy.feec.mass.WeightedMassOperators`) and the discrete curl :math:`\mathbb C` + (``derham.curl``). It conserves the energy + :math:`\frac 12 \mathbf e^\top \mathbb M_1 \mathbf e + \frac 12 \mathbf b^\top \mathbb M_2 \mathbf b` + and the discrete divergence :math:`\mathbb D \mathbf b` (since :math:`\mathbb D \mathbb C = 0`). + + :ref:`time_discret`: + + - implicit: Crank-Nicolson (implicit mid-point), + + .. math:: + + \mathbb M_1 \frac{\mathbf e^{n+1} - \mathbf e^n}{\Delta t} = \mathbb C^\top \mathbb M_2 \frac{\mathbf b^{n+1} + \mathbf b^n}{2}\,, + \qquad + \frac{\mathbf b^{n+1} - \mathbf b^n}{\Delta t} = - \mathbb C \frac{\mathbf e^{n+1} + \mathbf e^n}{2}\,, + + which conserves both invariants exactly (up to the solver tolerance). The magnetic field is eliminated by + :class:`~struphy.linear_algebra.schur_solver.SchurSolver` with :math:`x = \mathbf e`, :math:`y = \mathbf b`, + :math:`A = \mathbb M_1`, :math:`B = -\frac 12 \mathbb C^\top \mathbb M_2` and :math:`C = \frac 12 \mathbb C`, + which gives the curl-curl system + + .. math:: + + \left(\mathbb M_1 + \frac{\Delta t^2}{4} \mathbb C^\top \mathbb M_2 \mathbb C\right) \mathbf e^{n+1} + = \left(\mathbb M_1 - \frac{\Delta t^2}{4} \mathbb C^\top \mathbb M_2 \mathbb C\right) \mathbf e^n + + \Delta t \, \mathbb C^\top \mathbb M_2 \mathbf b^n\,, + + solved with ``solver`` and ``precond`` (built from :math:`\mathbb M_1`). Then + :math:`\mathbf b^{n+1} = \mathbf b^n - \frac{\Delta t}{2} \mathbb C (\mathbf e^{n+1} + \mathbf e^n)`. - System size reduction via :class:`~struphy.linear_algebra.schur_solver.SchurSolver`. + - explicit: Runge-Kutta method from :class:`~struphy.ode.utils.ButcherTableau` for the system + :math:`\dot{\mathbf e} = \mathbb M_1^{-1} \mathbb C^\top \mathbb M_2 \mathbf b`, + :math:`\dot{\mathbf b} = - \mathbb C \mathbf e`, where :math:`\mathbb M_1^{-1}` is applied with ``solver`` and + ``precond`` in each stage. The time step is limited by the CFL condition of light waves. """ class Variables: diff --git a/src/struphy/propagators/vlasov_ampere_coupling.py b/src/struphy/propagators/vlasov_ampere_coupling.py index 02bf5ed52..00be6cc72 100644 --- a/src/struphy/propagators/vlasov_ampere_coupling.py +++ b/src/struphy/propagators/vlasov_ampere_coupling.py @@ -31,9 +31,58 @@ class VlasovAmpereCoupling(Propagator): \\[2mm] \frac{\partial f}{\partial t} + \frac{1}{\varepsilon}\, \mathbf{E} \cdot \frac{\partial f}{\partial \mathbf{v}} &= 0 . - Time discretization: Crank-Nicolson (implicit mid-point). + With :math:`f` approximated by :math:`N` markers with positions :math:`\boldsymbol \eta_p` (logical, fixed in this + step), velocities :math:`\mathbf v_p` and weights :math:`w_p` (including the factor :math:`1/N`), and the FE + coefficients :math:`\mathbf e \in \mathbb R^{N_1}` of :math:`\mathbf E \in V^1_h`, the semi-discrete system reads - System size reduction via :class:`~struphy.linear_algebra.schur_solver.SchurSolver`. + .. math:: + + \mathbb M_1 \dot{\mathbf e} = - \frac{\alpha^2}{\varepsilon} \, (\mathbf \Lambda^1)^\top \, \overline{DF^{-1}} \, \mathbb W \, \mathbf V \,, + \qquad + \dot{\mathbf V} = \frac{1}{\varepsilon} \, \overline{DF^{-\top}} \, \mathbf \Lambda^1 \, \mathbf e\,, + + where :math:`\mathbf V = (\mathbf v_p)_p \in \mathbb R^{3N}`, :math:`\mathbb W = \textnormal{diag}(w_p)` + (acting on each velocity component), :math:`\mathbf \Lambda^1 \in \mathbb R^{3N \times N_1}` holds the + :math:`V^1_h` basis functions evaluated at the markers, :math:`(\mathbf \Lambda^1 \mathbf e)_p = \hat{\mathbf E}^1(\boldsymbol \eta_p)`, + and :math:`\overline{DF^{-1}}` is block-diagonal with the inverse Jacobian matrices :math:`DF^{-1}(\boldsymbol \eta_p)`. + It conserves the energy + :math:`\frac 12 \mathbf e^\top \mathbb M_1 \mathbf e + \frac{\alpha^2}{2} \mathbf V^\top \mathbb W \mathbf V`. + + :ref:`time_discret`: Crank-Nicolson (implicit mid-point), + + .. math:: + + \mathbb M_1 \frac{\mathbf e^{n+1} - \mathbf e^n}{\Delta t} = - \frac{\alpha^2}{\varepsilon} \, (\mathbf \Lambda^1)^\top \, \overline{DF^{-1}} \, \mathbb W \, \frac{\mathbf V^{n+1} + \mathbf V^n}{2}\,, + \qquad + \frac{\mathbf V^{n+1} - \mathbf V^n}{\Delta t} = \frac{1}{\varepsilon} \, \overline{DF^{-\top}} \, \mathbf \Lambda^1 \, \frac{\mathbf e^{n+1} + \mathbf e^n}{2}\,, + + which conserves the energy exactly (up to the solver tolerance; for constant weights, i.e. without control + variate). The velocities are eliminated by :class:`~struphy.linear_algebra.schur_solver.SchurSolver` with + :math:`x = \mathbf e`, :math:`y = \mathbf V`, :math:`A = \mathbb M_1`, + :math:`B = \frac{\alpha^2}{2 \varepsilon} (\mathbf \Lambda^1)^\top \overline{DF^{-1}} \mathbb W` and + :math:`C = - \frac{1}{2 \varepsilon} \overline{DF^{-\top}} \mathbf \Lambda^1`, which gives + + .. math:: + + \left(\mathbb M_1 + \frac{\Delta t^2 \alpha^2}{4 \varepsilon^2} \mathbb A\right) \mathbf e^{n+1} + = \left(\mathbb M_1 - \frac{\Delta t^2 \alpha^2}{4 \varepsilon^2} \mathbb A\right) \mathbf e^n + - \Delta t \frac{\alpha^2}{\varepsilon} \, \mathbf j\,, + + with the accumulated matrix and vector (:class:`~struphy.pic.accumulation.particles_to_grid.Accumulator`, + kernel :func:`~struphy.pic.accumulation.kernels.vlasov_maxwell.vlasov_maxwell`) + + .. math:: + + \mathbb A = (\mathbf \Lambda^1)^\top \, \overline{DF^{-1}} \, \mathbb W \, \overline{DF^{-\top}} \, \mathbf \Lambda^1 + = \sum_p w_p \, \mathbf \Lambda^1_p{}^\top G^{-1}(\boldsymbol \eta_p) \, \mathbf \Lambda^1_p\,, + \qquad + \mathbf j = (\mathbf \Lambda^1)^\top \, \overline{DF^{-1}} \, \mathbb W \, \mathbf V^n\,. + + :math:`\mathbb A` is a particle approximation of the :math:`n`-weighted mass matrix :math:`\mathbb M_1`, + so the Schur operator is a mass matrix with weight :math:`1 + \Delta t^2 \alpha^2 n / (4 \varepsilon^2)`, + solved with ``solver`` and ``precond`` (built from :math:`\mathbb M_1`). The velocities are then updated by + :func:`~struphy.pic.pushing.kernels.push_v_with_efield.push_v_with_efield` with the mid-point field + :math:`\frac 12 (\mathbf e^{n+1} + \mathbf e^n)`. """ class Variables: From 693a4b4c842476ec095a81f0d28480089814071c Mon Sep 17 00:00:00 2001 From: Stefan Possanner Date: Fri, 23 Oct 2026 08:39:43 +0200 Subject: [PATCH 2/6] MaxwellWeakAmpere: StiffnessPreconditioner as default precond [preview docs] The Schur operator M1 + dt^2/4 C^T M2 C is preconditioned by StiffnessPreconditioner(curl, mid=dt^2/4 M2, stab=M1), built per time step in prepare_time_step. Iterations hardly depend on dt (hollow cylinder 16x16x4, p=3: 39 instead of 411 PCG iterations at dt = 100 h); for dt <~ h the mass preconditioner stays slightly better. algo="explicit" falls back to the mass preconditioner, solver="cg" drops the preconditioner with a warning. New test file test_maxwell_weak_ampere.py. Co-Authored-By: Claude Opus 5.5 --- src/struphy/io/options.py | 3 + .../propagators/maxwell_weak_ampere.py | 61 ++++++-- .../tests/test_maxwell_weak_ampere.py | 140 ++++++++++++++++++ 3 files changed, 193 insertions(+), 11 deletions(-) create mode 100644 src/struphy/propagators/tests/test_maxwell_weak_ampere.py diff --git a/src/struphy/io/options.py b/src/struphy/io/options.py index fbde739ce..1c3c19948 100644 --- a/src/struphy/io/options.py +++ b/src/struphy/io/options.py @@ -80,6 +80,9 @@ class LiteralOptions: OptsShearAlfvenPrecond = Literal[ "MassMatrixPreconditioner", "PolarMassMatrixPreconditioner", "MultiGrid", "StiffnessPreconditioner", None ] + OptsMaxwellPrecond = Literal[ + "StiffnessPreconditioner", "MassMatrixPreconditioner", "PolarMassMatrixPreconditioner", None + ] OptsShearAlfvenCurlCurlPrecond = Literal[ "HiptmairXu", "MassMatrixPreconditioner", "PolarMassMatrixPreconditioner", None ] diff --git a/src/struphy/propagators/maxwell_weak_ampere.py b/src/struphy/propagators/maxwell_weak_ampere.py index 777753235..36e1bd237 100644 --- a/src/struphy/propagators/maxwell_weak_ampere.py +++ b/src/struphy/propagators/maxwell_weak_ampere.py @@ -66,13 +66,21 @@ class MaxwellWeakAmpere(Propagator): = \left(\mathbb M_1 - \frac{\Delta t^2}{4} \mathbb C^\top \mathbb M_2 \mathbb C\right) \mathbf e^n + \Delta t \, \mathbb C^\top \mathbb M_2 \mathbf b^n\,, - solved with ``solver`` and ``precond`` (built from :math:`\mathbb M_1`). Then + solved with ``solver`` and ``precond``. By default (``precond="StiffnessPreconditioner"``) the Schur operator + is preconditioned by :class:`~struphy.feec.preconditioner.StiffnessPreconditioner` for + :math:`\mathbb C^\top \left(\frac{\Delta t^2}{4} \mathbb M_2\right) \mathbb C + \mathbb M_1`, built for the + current time step (:meth:`prepare_time_step`). Its number of iterations hardly depends on :math:`\Delta t`, + whereas with the mass matrix preconditioner of :math:`\mathbb M_1` it grows with + :math:`\Delta t / h` (16x16x4 cells, p=3, hollow cylinder, :math:`\Delta t = 10 h`: 35 instead of 257 PCG + iterations, :math:`\Delta t = 100 h`: 39 instead of 411). For :math:`\Delta t \lesssim h` the mass matrix + preconditioner needs fewer iterations. Then :math:`\mathbf b^{n+1} = \mathbf b^n - \frac{\Delta t}{2} \mathbb C (\mathbf e^{n+1} + \mathbf e^n)`. - explicit: Runge-Kutta method from :class:`~struphy.ode.utils.ButcherTableau` for the system :math:`\dot{\mathbf e} = \mathbb M_1^{-1} \mathbb C^\top \mathbb M_2 \mathbf b`, :math:`\dot{\mathbf b} = - \mathbb C \mathbf e`, where :math:`\mathbb M_1^{-1}` is applied with ``solver`` and - ``precond`` in each stage. The time step is limited by the CFL condition of light waves. + ``precond`` in each stage (``"StiffnessPreconditioner"`` is replaced by the mass matrix preconditioner of + :math:`\mathbb M_1` here). The time step is limited by the CFL condition of light waves. """ class Variables: @@ -115,7 +123,7 @@ def __init__(self): @dataclass(repr=False) class Options(OptionsBase): - """Configuration options for :class:`MaxwellWeakAmpere`. + r"""Configuration options for :class:`MaxwellWeakAmpere`. Parameters ---------- @@ -129,9 +137,21 @@ class Options(OptionsBase): Name of the symmetric iterative solver passed to :func:`psydac.linalg.solvers.inverse`. - precond : LiteralOptions.OptsMassPrecond, default="MassMatrixPreconditioner" - Name of the preconditioner configuration. - Currently used only for implicit time stepping. + precond : LiteralOptions.OptsMaxwellPrecond, default="StiffnessPreconditioner" + Preconditioner of the linear solves. + + - ``"StiffnessPreconditioner"``: Kronecker (fast diagonalization) approximation of the inverse of the + Schur operator :math:`\mathbb M_1 + \frac{\Delta t^2}{4} \mathbb C^\top \mathbb M_2 \mathbb C` + (:class:`~struphy.feec.preconditioner.StiffnessPreconditioner` with ``derivative="curl"``; with polar + splines a :class:`~struphy.feec.preconditioner.PolarStiffnessPreconditioner`), rebuilt when the time + step changes. For ``algo="explicit"``, the mass matrix preconditioner of :math:`\mathbb M_1` is used + instead. (With polar splines, the implicit step currently fails in the + :class:`~struphy.linear_algebra.schur_solver.SchurSolver` with any preconditioner.) + - ``"MassMatrixPreconditioner"``, ``"PolarMassMatrixPreconditioner"``: preconditioner of + :math:`\mathbb M_1`, better for small time steps :math:`\Delta t \lesssim h`, e.g. in CFL-limited + PIC runs. + - ``None``: no preconditioner (set automatically, with a warning, for ``solver="cg"``, which takes no + preconditioner). solver_params : SolverParameters, default=None Iterative-solver controls (for example ``tol``, ``maxiter``, @@ -155,7 +175,7 @@ class Options(OptionsBase): # propagator options algo: OptsAlgo = "implicit" solver: LiteralOptions.OptsSymmSolver = "pcg" - precond: LiteralOptions.OptsMassPrecond = "MassMatrixPreconditioner" + precond: LiteralOptions.OptsMaxwellPrecond = "StiffnessPreconditioner" solver_params: SolverParameters = None butcher: ButcherTableau = None @@ -163,7 +183,10 @@ def __post_init__(self): # checks check_option(self.algo, self.OptsAlgo) check_option(self.solver, LiteralOptions.OptsSymmSolver) - check_option(self.precond, LiteralOptions.OptsMassPrecond) + check_option(self.precond, LiteralOptions.OptsMaxwellPrecond) + if self.solver == "cg" and self.precond is not None: + logger.warning(f'solver="cg" takes no preconditioner; ignoring precond="{self.precond}".') + self.precond = None # defaults if self.solver_params is None: @@ -191,11 +214,15 @@ def allocate(self): M2 = self.mass_ops.M2 curl = self.derham.curl - # Preconditioner for M1 + ... - if self.options.precond is None: + # Preconditioner (the one of the Schur operator depends on dt and is built in prepare_time_step) + self._precond = self.options.precond + self._schur_precond_dt = None + if self._precond == "StiffnessPreconditioner" and self.options.algo == "explicit": + self._precond = "PolarMassMatrixPreconditioner" if self.derham.polar_splines else "MassMatrixPreconditioner" + if self._precond in (None, "StiffnessPreconditioner"): pc = None else: - pc_class = getattr(preconditioner, self.options.precond) + pc_class = getattr(preconditioner, self._precond) pc = pc_class(M1) if self.options.algo == "implicit": @@ -257,6 +284,15 @@ def f2(t, y1, y2, out: BlockVector = out2): self._e_tmp2 = self.variables.e.spline.vector.space.zeros() self._b_tmp1 = self.variables.b.spline.vector.space.zeros() + def prepare_time_step(self, dt: float): + """Stiffness preconditioner of the Schur operator M1 + dt^2/4 C^T M2 C for the time step ``dt``: built in the + setup, rebuilt when dt changes.""" + if self._precond == "StiffnessPreconditioner" and dt != self._schur_precond_dt: + self._schur_solver.precond = preconditioner.StiffnessPreconditioner( + self.mass_ops, "curl", sigma=1.0, mid=(dt**2 / 4) * self.mass_ops.M2, stab="M1" + ) + self._schur_precond_dt = dt + @profile def __call__(self, dt): # current FE coeffs @@ -264,6 +300,9 @@ def __call__(self, dt): bn = self.variables.b.spline.vector if self.options.algo == "implicit": + # preconditioner for the current time step (no-op unless dt changed since the setup) + self.prepare_time_step(dt) + # solve for new e coeffs self._B.dot(bn, out=self._byn) diff --git a/src/struphy/propagators/tests/test_maxwell_weak_ampere.py b/src/struphy/propagators/tests/test_maxwell_weak_ampere.py new file mode 100644 index 000000000..4b73b95b7 --- /dev/null +++ b/src/struphy/propagators/tests/test_maxwell_weak_ampere.py @@ -0,0 +1,140 @@ +import cunumpy as xp +import pytest +from maybempi import MPI + +from struphy import domains +from struphy.feec.mass import WeightedMassOperators +from struphy.feec.psydac_derham import Derham +from struphy.io.options import DerhamOptions +from struphy.linear_algebra.solver import SolverParameters +from struphy.models.variables import FEECVariable +from struphy.propagators.base import Propagator +from struphy.propagators.maxwell_weak_ampere import MaxwellWeakAmpere +from struphy.topology.grids import TensorProductGrid + +comm = MPI.COMM_WORLD + + +def _setup(mapping): + polar = mapping == "IGAPolarCylinder" + if mapping == "Cuboid": + dom = domains.Cuboid() + elif mapping == "HollowCylinder": + dom = domains.HollowCylinder(a1=0.1, a2=1.0, Lz=3.0) + else: + dom = domains.IGAPolarCylinder(num_elements=(12, 12), degree=(2, 2), a=1.0, Lz=3.0) + derham = Derham( + grid=TensorProductGrid(num_elements=(12, 12, 4)), + options=DerhamOptions( + degree=(2, 2, 1), + bcs=(("free", "dirichlet") if polar else ("dirichlet", "dirichlet"), None, None), + polar_splines=polar, + ), + comm=comm, + domain=dom, + ) + Propagator.derham = derham + Propagator.domain = dom + Propagator.mass_ops = WeightedMassOperators(derham=derham, domain=dom) + return derham, dom + + +def _propagator(derham, dom, seed=1, **options): + e = FEECVariable(space="Hcurl") + e.allocate(derham=derham, domain=dom) + b = FEECVariable(space="Hdiv") + b.allocate(derham=derham, domain=dom) + # smooth logical components (projected, so that they are also in the polar spaces) + funcs = [ + lambda e1, e2, e3, k=k: ( + xp.sin(xp.pi * e1) * xp.cos(2 * xp.pi * (e2 + k * seed / 7)) * xp.sin(2 * xp.pi * e3 + k) + ) + for k in range(3) + ] + e.spline.vector = derham.boundary_ops["1"].dot(derham.P1(funcs)) + b.spline.vector = derham.boundary_ops["2"].dot(derham.P2(funcs[::-1])) + prop = MaxwellWeakAmpere() + prop.variables.e = e + prop.variables.b = b + prop.options = prop.Options(**options) + prop.allocate() + return prop + + +def _energy(prop): + e = prop.variables.e.spline.vector + b = prop.variables.b.spline.vector + return 0.5 * e.inner(prop.mass_ops.M1.dot(e)) + 0.5 * b.inner(prop.mass_ops.M2.dot(b)) + + +@pytest.mark.parametrize("mapping", ["Cuboid", "HollowCylinder"]) +def test_stiffness_preconditioner(mapping): + """Implicit step with the stiffness and the mass matrix preconditioner: same solution, exact energy conservation, + and fewer iterations with the stiffness preconditioner for dt >> h.""" + derham, dom = _setup(mapping) + dt = 1.0 + solver_params = SolverParameters(tol=1e-12, maxiter=3000) + + niter = {} + solution = {} + mass_pc = "PolarMassMatrixPreconditioner" if derham.polar_splines else "MassMatrixPreconditioner" + for precond in (mass_pc, "StiffnessPreconditioner"): + prop = _propagator(derham, dom, precond=precond, solver_params=solver_params) + prop.prepare_time_step(dt) + e = prop.variables.e.spline.vector + b = prop.variables.b.spline.vector + _, info = prop._schur_solver(e, prop._B.dot(b), dt) + assert info["success"] + niter[precond] = info["niter"] + + energy_0 = _energy(prop) + prop(dt) + assert abs(_energy(prop) - energy_0) < 1e-8 * energy_0 + solution[precond] = prop.variables.e.spline.vector.copy() + + diff = solution["StiffnessPreconditioner"] - solution[mass_pc] + assert diff.inner(diff) ** 0.5 < 1e-8 * solution[mass_pc].inner(solution[mass_pc]) ** 0.5 + assert niter["StiffnessPreconditioner"] < niter[mass_pc] + + +def test_stiffness_preconditioner_dt_change(): + """The stiffness preconditioner is rebuilt only when dt changes.""" + derham, dom = _setup("Cuboid") + prop = _propagator(derham, dom) + prop(0.1) + pc = prop._schur_solver.precond + prop(0.1) + assert prop._schur_solver.precond is pc + prop(0.2) + assert prop._schur_solver.precond is not pc + assert prop._schur_precond_dt == 0.2 + + +def test_polar_stiffness_preconditioner(): + """With polar splines, the default preconditioner is a PolarStiffnessPreconditioner. (No solve: the implicit step + with polar splines currently fails in SchurSolver with any preconditioner.)""" + from struphy.feec.preconditioner import PolarStiffnessPreconditioner + + derham, dom = _setup("IGAPolarCylinder") + prop = _propagator(derham, dom) + prop.prepare_time_step(1.0) + assert isinstance(prop._schur_solver.precond, PolarStiffnessPreconditioner) + + +def test_explicit_default_precond(): + """For algo="explicit", the default stiffness preconditioner falls back to the mass matrix preconditioner.""" + derham, dom = _setup("Cuboid") + prop = _propagator(derham, dom, algo="explicit") + assert prop._precond == "MassMatrixPreconditioner" + energy_0 = _energy(prop) + prop(0.01) + assert abs(_energy(prop) - energy_0) < 1e-3 * energy_0 + + +@pytest.mark.mpi_skip +def test_precond_options(): + """Invalid preconditioners are rejected; solver="cg" takes no preconditioner.""" + with pytest.raises(AssertionError): + MaxwellWeakAmpere.Options(precond="HiptmairXu") + assert MaxwellWeakAmpere.Options().precond == "StiffnessPreconditioner" + assert MaxwellWeakAmpere.Options(solver="cg").precond is None From f98a69526f4656c645bf108f191ad70a3c09fe48 Mon Sep 17 00:00:00 2001 From: Stefan Possanner Date: Fri, 23 Oct 2026 09:30:20 +0200 Subject: [PATCH 3/6] Docs: working links in MaxwellWeakAmpere/VlasovAmpereCoupling, document preconditioners and solvers [preview docs] - remove the iteration-count sentence from the MaxwellWeakAmpere docstring - full-path cross-references (Derham.curl, SchurSolver, prepare_time_step, kernel functions in their *_kernels modules); feectools inverse as literal - dev-feec.rst: new sections Preconditioners (struphy.feec.preconditioner, struphy.polar.kron_2dx1d, multigrid preconditioner) and Linear solvers (schur_solver, solver) - dev-pic.rst: vlasov_maxwell accumulation kernel, new Pusher kernels section with push_v_with_efield Co-Authored-By: Claude Opus 5.5 --- doc/sections/subsections/dev-feec.rst | 37 +++++++++++++++++++ doc/sections/subsections/dev-pic.rst | 8 ++++ .../propagators/maxwell_weak_ampere.py | 14 +++---- .../propagators/vlasov_ampere_coupling.py | 10 ++--- 4 files changed, 55 insertions(+), 14 deletions(-) diff --git a/doc/sections/subsections/dev-feec.rst b/doc/sections/subsections/dev-feec.rst index 98ae440b1..6cebc7079 100644 --- a/doc/sections/subsections/dev-feec.rst +++ b/doc/sections/subsections/dev-feec.rst @@ -46,3 +46,40 @@ Kronecker approximations of operators :special-members: :show-inheritance: :exclude-members: __init__ + + +Preconditioners +^^^^^^^^^^^^^^^ + +.. automodule:: struphy.feec.preconditioner + :members: + :show-inheritance: + :exclude-members: __init__ + +Polar 2D x 1D Kronecker solvers used by the polar preconditioners: + +.. automodule:: struphy.polar.kron_2dx1d + :members: + :show-inheritance: + :exclude-members: __init__ + +Geometric multigrid preconditioner: + +.. automodule:: struphy.linear_algebra.multigrid.preconditioner + :members: + :show-inheritance: + :exclude-members: __init__ + + +Linear solvers +^^^^^^^^^^^^^^ + +.. automodule:: struphy.linear_algebra.schur_solver + :members: + :show-inheritance: + :exclude-members: __init__ + +.. automodule:: struphy.linear_algebra.solver + :members: + :show-inheritance: + :exclude-members: __init__ diff --git a/doc/sections/subsections/dev-pic.rst b/doc/sections/subsections/dev-pic.rst index 9192a6686..993614ce8 100644 --- a/doc/sections/subsections/dev-pic.rst +++ b/doc/sections/subsections/dev-pic.rst @@ -55,6 +55,8 @@ Accumulation kernels .. automodule:: struphy.pic.accumulation.kernels +.. autofunction:: struphy.pic.accumulation.kernels.vlasov_maxwell.vlasov_maxwell_kernels.vlasov_maxwell + Pusher class ^^^^^^^^^^^^ @@ -101,3 +103,9 @@ are zero. Initialization setups must use ``alpha=0`` (the default). :special-members: :show-inheritance: :exclude-members: __init__ + + +Pusher kernels +^^^^^^^^^^^^^^ + +.. autofunction:: struphy.pic.pushing.kernels.push_v_with_efield.push_v_with_efield_kernels.push_v_with_efield diff --git a/src/struphy/propagators/maxwell_weak_ampere.py b/src/struphy/propagators/maxwell_weak_ampere.py index 36e1bd237..88caeee23 100644 --- a/src/struphy/propagators/maxwell_weak_ampere.py +++ b/src/struphy/propagators/maxwell_weak_ampere.py @@ -41,7 +41,7 @@ class MaxwellWeakAmpere(Propagator): with the mass matrices :math:`\mathbb M_1`, :math:`\mathbb M_2` (:class:`~struphy.feec.mass.WeightedMassOperators`) and the discrete curl :math:`\mathbb C` - (``derham.curl``). It conserves the energy + (:attr:`~struphy.feec.psydac_derham.Derham.curl`). It conserves the energy :math:`\frac 12 \mathbf e^\top \mathbb M_1 \mathbf e + \frac 12 \mathbf b^\top \mathbb M_2 \mathbf b` and the discrete divergence :math:`\mathbb D \mathbf b` (since :math:`\mathbb D \mathbb C = 0`). @@ -69,11 +69,7 @@ class MaxwellWeakAmpere(Propagator): solved with ``solver`` and ``precond``. By default (``precond="StiffnessPreconditioner"``) the Schur operator is preconditioned by :class:`~struphy.feec.preconditioner.StiffnessPreconditioner` for :math:`\mathbb C^\top \left(\frac{\Delta t^2}{4} \mathbb M_2\right) \mathbb C + \mathbb M_1`, built for the - current time step (:meth:`prepare_time_step`). Its number of iterations hardly depends on :math:`\Delta t`, - whereas with the mass matrix preconditioner of :math:`\mathbb M_1` it grows with - :math:`\Delta t / h` (16x16x4 cells, p=3, hollow cylinder, :math:`\Delta t = 10 h`: 35 instead of 257 PCG - iterations, :math:`\Delta t = 100 h`: 39 instead of 411). For :math:`\Delta t \lesssim h` the mass matrix - preconditioner needs fewer iterations. Then + current time step (:meth:`~struphy.propagators.maxwell_weak_ampere.MaxwellWeakAmpere.prepare_time_step`). Then :math:`\mathbf b^{n+1} = \mathbf b^n - \frac{\Delta t}{2} \mathbb C (\mathbf e^{n+1} + \mathbf e^n)`. - explicit: Runge-Kutta method from :class:`~struphy.ode.utils.ButcherTableau` for the system @@ -123,7 +119,7 @@ def __init__(self): @dataclass(repr=False) class Options(OptionsBase): - r"""Configuration options for :class:`MaxwellWeakAmpere`. + r"""Configuration options for :class:`~struphy.propagators.maxwell_weak_ampere.MaxwellWeakAmpere`. Parameters ---------- @@ -135,7 +131,7 @@ class Options(OptionsBase): solver : LiteralOptions.OptsSymmSolver, default="pcg" Name of the symmetric iterative solver passed to - :func:`psydac.linalg.solvers.inverse`. + ``feectools.linalg.solvers.inverse``. precond : LiteralOptions.OptsMaxwellPrecond, default="StiffnessPreconditioner" Preconditioner of the linear solves. @@ -165,7 +161,7 @@ class Options(OptionsBase): Notes ----- - System size reduction is performed via :class:`SchurSolver` for + System size reduction is performed via :class:`~struphy.linear_algebra.schur_solver.SchurSolver` for implicit time stepping to eliminate the magnetic field and reduce the system to the electric field equation. """ diff --git a/src/struphy/propagators/vlasov_ampere_coupling.py b/src/struphy/propagators/vlasov_ampere_coupling.py index 00be6cc72..65c9cc282 100644 --- a/src/struphy/propagators/vlasov_ampere_coupling.py +++ b/src/struphy/propagators/vlasov_ampere_coupling.py @@ -69,7 +69,7 @@ class VlasovAmpereCoupling(Propagator): - \Delta t \frac{\alpha^2}{\varepsilon} \, \mathbf j\,, with the accumulated matrix and vector (:class:`~struphy.pic.accumulation.particles_to_grid.Accumulator`, - kernel :func:`~struphy.pic.accumulation.kernels.vlasov_maxwell.vlasov_maxwell`) + kernel :func:`~struphy.pic.accumulation.kernels.vlasov_maxwell.vlasov_maxwell_kernels.vlasov_maxwell`) .. math:: @@ -81,12 +81,12 @@ class VlasovAmpereCoupling(Propagator): :math:`\mathbb A` is a particle approximation of the :math:`n`-weighted mass matrix :math:`\mathbb M_1`, so the Schur operator is a mass matrix with weight :math:`1 + \Delta t^2 \alpha^2 n / (4 \varepsilon^2)`, solved with ``solver`` and ``precond`` (built from :math:`\mathbb M_1`). The velocities are then updated by - :func:`~struphy.pic.pushing.kernels.push_v_with_efield.push_v_with_efield` with the mid-point field + :func:`~struphy.pic.pushing.kernels.push_v_with_efield.push_v_with_efield_kernels.push_v_with_efield` with the mid-point field :math:`\frac 12 (\mathbf e^{n+1} + \mathbf e^n)`. """ class Variables: - """Container for variables advanced by :class:`VlasovAmpereCoupling`. + """Container for variables advanced by :class:`~struphy.propagators.vlasov_ampere_coupling.VlasovAmpereCoupling`. Attributes ---------- @@ -125,12 +125,12 @@ def __init__(self): @dataclass(repr=False) class Options(OptionsBase): - """Configuration options for :class:`VlasovAmpereCoupling`. + """Configuration options for :class:`~struphy.propagators.vlasov_ampere_coupling.VlasovAmpereCoupling`. Parameters ---------- solver : LiteralOptions.OptsSymmSolver, default="pcg" - Symmetric iterative solver used by :class:`SchurSolver`. + Symmetric iterative solver used by :class:`~struphy.linear_algebra.schur_solver.SchurSolver`. precond : LiteralOptions.OptsMassPrecond, default="MassMatrixPreconditioner" Preconditioner for electric-field mass matrix block. solver_params : SolverParameters, default=None From 3edcd69eb628bffed1e171e9c2e73851e8d187e2 Mon Sep 17 00:00:00 2001 From: Stefan Possanner Date: Fri, 23 Oct 2026 09:35:35 +0200 Subject: [PATCH 4/6] Fix the SchurSolver docstrings (broken math, parameters) [preview docs] The raw docstrings of SchurSolver, SchurSolverFull and SchurSolverFull3 used \\text and \\right (a LaTeX line break in a raw string), which broke the math; rewritten with pmatrix. Markdown links replaced by plain names, and the SchurSolver Parameters section now matches the signature (precond, solver_params). Co-Authored-By: Claude Opus 5.5 --- src/struphy/linear_algebra/schur_solver.py | 95 +++++++++------------- 1 file changed, 39 insertions(+), 56 deletions(-) diff --git a/src/struphy/linear_algebra/schur_solver.py b/src/struphy/linear_algebra/schur_solver.py index 586421cf1..177afe8ea 100644 --- a/src/struphy/linear_algebra/schur_solver.py +++ b/src/struphy/linear_algebra/schur_solver.py @@ -12,42 +12,38 @@ class SchurSolver: .. math:: - \left( \matrix{ - A & \Delta t B \cr - \Delta t C & \\text{Id} - } \\right) - \left( \matrix{ - x^{n+1} \cr y^{n+1} - } \\right) + \begin{pmatrix} A & \Delta t\, B \\ \Delta t\, C & \mathbb I \end{pmatrix} + \begin{pmatrix} x^{n+1} \\ y^{n+1} \end{pmatrix} = - \left( \matrix{ - A & - \Delta t B \cr - - \Delta t C & \\text{Id} - } \\right) - \left( \matrix{ - x^n \cr y^n - } \\right) + \begin{pmatrix} A & - \Delta t\, B \\ - \Delta t\, C & \mathbb I \end{pmatrix} + \begin{pmatrix} x^n \\ y^n \end{pmatrix} - using the Schur complement :math:`S = A - \Delta t^2 BC`, where Id is the identity matrix - and :math:`(x^n, y^n)` is given. The solution is given by + using the Schur complement :math:`S = A - \Delta t^2 BC`, where :math:`\mathbb I` is the identity matrix + and :math:`(x^n, y^n)` is given. The solution is .. math:: - x^{n+1} = S^{-1} \left[ (A + \Delta t^2 BC) \, x^n - 2 \Delta t B \, y^n \\right] \,. + x^{n+1} = S^{-1} \left[ (A + \Delta t^2 BC) \, x^n - 2 \Delta t\, B \, y^n \right] \,. + + The caller passes :math:`B y^n` and computes :math:`y^{n+1} = y^n - \Delta t\, C (x^{n+1} + x^n)` itself. Parameters ---------- A : LinearOperator - Upper left block from [[A B], [C Id]]. + Upper left block :math:`A`. BC : LinearOperator - Product from [[A B], [C Id]]. + Product :math:`BC` of the off-diagonal blocks (same domain and codomain as :math:`A`). solver_name : str - See [feectools.linalg.solvers](https://github.com/pyccel/psydac/blob/535717c6f5ea328aacbbbbcc2d582a92b31c9377/psydac/linalg/solvers.py#L47) for possible names. + Name of the iterative solver, passed to ``feectools.linalg.solvers.inverse`` (e.g. ``"pcg"``, ``"cg"``, + ``"pbicgstab"``, ``"gmres"``). - **solver_params : - Must correspond to the chosen solver. + precond : LinearOperator, optional + Preconditioner of the Schur complement :math:`S` (can be replaced later with :attr:`precond`). + + solver_params : SolverParameters, optional + Solver controls (``tol``, ``maxiter``, ...); defaults to ``SolverParameters()``. """ def __init__( @@ -186,37 +182,31 @@ class SchurSolverFull: .. math:: - \left( \matrix{ - A & B \cr - C & \\text{Id} - } \\right) - \left( \matrix{ - x \cr y - } \\right) + \begin{pmatrix} A & B \\ C & \mathbb I \end{pmatrix} + \begin{pmatrix} x \\ y \end{pmatrix} = - \left( \matrix{ - b_x \cr b_y - } \\right) + \begin{pmatrix} b_x \\ b_y \end{pmatrix} - using the Schur complement :math:`S = A - BC`, where Id is the identity matrix - and :math:`(b_x, b_y)^T` is given. The solution is given by + using the Schur complement :math:`S = A - BC`, where :math:`\mathbb I` is the identity matrix + and :math:`(b_x, b_y)^\top` is given. The solution is .. math:: x &= S^{-1} \, (b_x - B b_y ) \,, - + \\ y &= b_y - C x \,. Parameters ---------- M : BlockLinearOperator - Matrix [[A B], [C Id]]. + Block matrix :math:`[[A, B], [C, \mathbb I]]`. solver_name : str - See [feectools.linalg.solvers](https://github.com/pyccel/psydac/blob/535717c6f5ea328aacbbbbcc2d582a92b31c9377/psydac/linalg/solvers.py#L47) for possible names. + Name of the iterative solver, passed to ``feectools.linalg.solvers.inverse`` (e.g. ``"pcg"``, ``"cg"``, + ``"pbicgstab"``, ``"gmres"``). **solver_params : - Must correspond to the chosen solver. + Keyword arguments of the chosen solver. """ def __init__(self, M, solver_name, **solver_params): @@ -295,40 +285,33 @@ class SchurSolverFull3: .. math:: - \left( \matrix{ - A & B & D \cr - C & \\text{Id} & 0 \cr - E & 0 & \\text{Id} \cr - } \\right) - \left( \matrix{ - x \cr y \cr z - } \\right) + \begin{pmatrix} A & B & D \\ C & \mathbb I & 0 \\ E & 0 & \mathbb I \end{pmatrix} + \begin{pmatrix} x \\ y \\ z \end{pmatrix} = - \left( \matrix{ - b_x \cr b_y \cr b_z - } \\right) + \begin{pmatrix} b_x \\ b_y \\ b_z \end{pmatrix} - using the Schur complement :math:`S = A - BC - DE`, where Id is the identity matrix - and :math:`(b_x, b_y, b_z)^T` is given. The solution is given by + using the Schur complement :math:`S = A - BC - DE`, where :math:`\mathbb I` is the identity matrix + and :math:`(b_x, b_y, b_z)^\top` is given. The solution is .. math:: x &= S^{-1} \, (b_x - B b_y - D b_z) \,, - + \\ y &= b_y - C x \,, - + \\ z &= b_z - E x \,. Parameters ---------- M : BlockLinearOperator - Matrix [[A B D], [C Id 0], [E 0 Id]]. + Block matrix :math:`[[A, B, D], [C, \mathbb I, 0], [E, 0, \mathbb I]]`. solver_name : str - See [feectools.linalg.solvers](https://github.com/pyccel/psydac/blob/535717c6f5ea328aacbbbbcc2d582a92b31c9377/psydac/linalg/solvers.py#L47) for possible names. + Name of the iterative solver, passed to ``feectools.linalg.solvers.inverse`` (e.g. ``"pcg"``, ``"cg"``, + ``"pbicgstab"``, ``"gmres"``). **solver_params : - Must correspond to the chosen solver. + Keyword arguments of the chosen solver. """ def __init__(self, M, solver_name, **solver_params): From b8d04fa90b18063742229e845124b2bf79c86960 Mon Sep 17 00:00:00 2001 From: Stefan Possanner Date: Fri, 23 Oct 2026 09:36:49 +0200 Subject: [PATCH 5/6] StiffnessPreconditioner docstring: references for fast diagonalization [preview docs] Lynch, Rice and Thomas (1964) for the method, Sangalli and Tani (2016) for preconditioning spline stiffness matrices; short description of the generalized eigendecomposition used by KroneckerSumSolver. Co-Authored-By: Claude Opus 5.5 --- src/struphy/feec/preconditioner.py | 16 ++++++++++++++-- 1 file changed, 14 insertions(+), 2 deletions(-) diff --git a/src/struphy/feec/preconditioner.py b/src/struphy/feec/preconditioner.py index cb75576eb..7dd3153f6 100644 --- a/src/struphy/feec/preconditioner.py +++ b/src/struphy/feec/preconditioner.py @@ -624,8 +624,12 @@ class StiffnessPreconditioner(KroneckerPreconditioner): with the 1d stiffness matrices :math:`S_d = \mathbb d_d^T M^D_d \mathbb d_d` in the directions :math:`N_c` in which component c is differentiated (its B-spline directions). These blocks - are inverted exactly by fast diagonalization - (:class:`~feectools.linalg.kron.KroneckerSumSolver`). For the gradient this is the exact + are inverted exactly by fast diagonalization (Lynch, Rice and Thomas 1964; for spline stiffness matrices + see Sangalli and Tani 2016, in the References below; implemented in + :class:`~feectools.linalg.kron.KroneckerSumSolver`): with the generalized eigendecompositions + :math:`S_d U_d = M_d U_d \Lambda_d`, :math:`U_d^\top M_d U_d = I`, the inverse of a Kronecker sum is + :math:`(U_1 \otimes U_2 \otimes U_3) \, D^{-1} \, (U_1 \otimes U_2 \otimes U_3)^\top` with a diagonal + :math:`D`, at the cost of a few 1d matrix products per direction. For the gradient this is the exact inverse on the logical cube. For curl and div, the off-diagonal blocks are neglected (block Jacobi, :math:`P_{BJ}`). @@ -734,6 +738,14 @@ class StiffnessPreconditioner(KroneckerPreconditioner): stab_average : str | None Gradient only: the directions averaged in the stabilization :math:`\mathbb S (\mathbb I - \mathbb A)`, e.g. ``"eta3"``, ``"eta2 eta3"`` or ``"eta1 eta2 eta3"`` (global mean); None (default) for no averaging. + + References + ---------- + R. E. Lynch, J. R. Rice and D. H. Thomas, Direct solution of partial difference equations by tensor product + methods, Numer. Math. 6, 185-199 (1964), `doi:10.1007/BF01386067 `_. + + G. Sangalli and M. Tani, Isogeometric preconditioners based on fast solvers for the Sylvester equation, + SIAM J. Sci. Comput. 38(6), A3644-A3671 (2016), `doi:10.1137/16M1062788 `_. """ _FORMS = {"grad": 0, "curl": 1, "div": 2} From 8dd6dc0fb71400c7506016ee5fd39838281d10ae Mon Sep 17 00:00:00 2001 From: Stefan Possanner Date: Fri, 23 Oct 2026 09:41:38 +0200 Subject: [PATCH 6/6] Review fixes: SchurSolver property docstrings, define D in fast diagonalization [preview docs] Co-Authored-By: Claude Opus 5.5 --- src/struphy/feec/preconditioner.py | 7 ++++--- src/struphy/linear_algebra/schur_solver.py | 12 ++++++------ 2 files changed, 10 insertions(+), 9 deletions(-) diff --git a/src/struphy/feec/preconditioner.py b/src/struphy/feec/preconditioner.py index 7dd3153f6..c3c9da0c7 100644 --- a/src/struphy/feec/preconditioner.py +++ b/src/struphy/feec/preconditioner.py @@ -627,9 +627,10 @@ class StiffnessPreconditioner(KroneckerPreconditioner): are inverted exactly by fast diagonalization (Lynch, Rice and Thomas 1964; for spline stiffness matrices see Sangalli and Tani 2016, in the References below; implemented in :class:`~feectools.linalg.kron.KroneckerSumSolver`): with the generalized eigendecompositions - :math:`S_d U_d = M_d U_d \Lambda_d`, :math:`U_d^\top M_d U_d = I`, the inverse of a Kronecker sum is - :math:`(U_1 \otimes U_2 \otimes U_3) \, D^{-1} \, (U_1 \otimes U_2 \otimes U_3)^\top` with a diagonal - :math:`D`, at the cost of a few 1d matrix products per direction. For the gradient this is the exact + :math:`S_d U_d = M_d U_d \Lambda_d`, :math:`U_d^T M_d U_d = I`, and :math:`U = U_1 \otimes U_2 \otimes U_3`, + the block is :math:`A_{cc} = U^{-T} D \, U^{-1}` with the diagonal matrix :math:`D = U^T A_{cc} U` (sums of + Kronecker products of the :math:`\Lambda_d` and identities, plus :math:`\sigma`), hence + :math:`A_{cc}^{-1} = U D^{-1} U^T`, at the cost of a few 1d matrix products per direction. For the gradient this is the exact inverse on the logical cube. For curl and div, the off-diagonal blocks are neglected (block Jacobi, :math:`P_{BJ}`). diff --git a/src/struphy/linear_algebra/schur_solver.py b/src/struphy/linear_algebra/schur_solver.py index 177afe8ea..a46b98fe8 100644 --- a/src/struphy/linear_algebra/schur_solver.py +++ b/src/struphy/linear_algebra/schur_solver.py @@ -87,12 +87,12 @@ def __init__( @property def A(self): - """Upper left block from [[A B], [C Id]].""" + """Upper left block :math:`A`.""" return self._A @property def BC(self): - """Product from [[A B], [C Id]].""" + """Product :math:`BC` of the off-diagonal blocks.""" return self._BC @property @@ -110,12 +110,12 @@ def precond(self, pc): @A.setter def A(self, a): - """Upper left block from [[A B], [C Id]].""" + """Upper left block :math:`A`.""" self._A = a @BC.setter def BC(self, bc): - """Product from [[A B], [C Id]].""" + """Product :math:`BC` of the off-diagonal blocks.""" self._BC = bc @profile @@ -206,7 +206,7 @@ class SchurSolverFull: ``"pbicgstab"``, ``"gmres"``). **solver_params : - Keyword arguments of the chosen solver. + Keyword arguments of the chosen solver; must contain ``pc`` (the preconditioner, None for none). """ def __init__(self, M, solver_name, **solver_params): @@ -311,7 +311,7 @@ class SchurSolverFull3: ``"pbicgstab"``, ``"gmres"``). **solver_params : - Keyword arguments of the chosen solver. + Keyword arguments of the chosen solver; must contain ``pc`` (the preconditioner, None for none). """ def __init__(self, M, solver_name, **solver_params):