MaxwellWeakAmpere: StiffnessPreconditioner as default; document the discretization of MaxwellWeakAmpere and VlasovAmpereCoupling

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__

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

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 <https://doi.org/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 <https://doi.org/10.1137/16M1062788>`_.
     """
 
     _FORMS = {"grad": 0, "curl": 1, "div": 2}

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
     ]

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):