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 58c3158d7..5273cc284 100644 --- a/doc/sections/subsections/dev-pic.rst +++ b/doc/sections/subsections/dev-pic.rst @@ -69,6 +69,8 @@ Accumulation kernels .. automodule:: struphy.pic.accumulation.kernels +.. autofunction:: struphy.pic.accumulation.kernels.vlasov_maxwell.vlasov_maxwell_kernels.vlasov_maxwell + Pusher class ^^^^^^^^^^^^ @@ -115,3 +117,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/feec/preconditioner.py b/src/struphy/feec/preconditioner.py index cb75576eb..c3c9da0c7 100644 --- a/src/struphy/feec/preconditioner.py +++ b/src/struphy/feec/preconditioner.py @@ -624,8 +624,13 @@ 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^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}`). @@ -734,6 +739,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} 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/linear_algebra/schur_solver.py b/src/struphy/linear_algebra/schur_solver.py index 586421cf1..a46b98fe8 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__( @@ -91,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 @@ -114,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 @@ -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; must contain ``pc`` (the preconditioner, None for none). """ 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; must contain ``pc`` (the preconditioner, None for none). """ def __init__(self, M, solver_name, **solver_params): diff --git a/src/struphy/propagators/maxwell_weak_ampere.py b/src/struphy/propagators/maxwell_weak_ampere.py index 9c5d45486..88caeee23 100644 --- a/src/struphy/propagators/maxwell_weak_ampere.py +++ b/src/struphy/propagators/maxwell_weak_ampere.py @@ -30,12 +30,53 @@ 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` + (: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`). + + :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 - System size reduction via :class:`~struphy.linear_algebra.schur_solver.SchurSolver`. + .. 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``. 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:`~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 + :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 (``"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: @@ -78,7 +119,7 @@ def __init__(self): @dataclass(repr=False) class Options(OptionsBase): - """Configuration options for :class:`MaxwellWeakAmpere`. + r"""Configuration options for :class:`~struphy.propagators.maxwell_weak_ampere.MaxwellWeakAmpere`. Parameters ---------- @@ -90,11 +131,23 @@ class Options(OptionsBase): solver : LiteralOptions.OptsSymmSolver, default="pcg" 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. + ``feectools.linalg.solvers.inverse``. + + 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``, @@ -108,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. """ @@ -118,7 +171,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 @@ -126,7 +179,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: @@ -154,11 +210,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": @@ -220,6 +280,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 @@ -227,6 +296,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 diff --git a/src/struphy/propagators/vlasov_ampere_coupling.py b/src/struphy/propagators/vlasov_ampere_coupling.py index f77cfade4..b3ce9cb54 100644 --- a/src/struphy/propagators/vlasov_ampere_coupling.py +++ b/src/struphy/propagators/vlasov_ampere_coupling.py @@ -33,13 +33,62 @@ 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_kernels.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_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 ---------- @@ -78,12 +127,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