Available Propagators

Field solvers

class struphy.propagators.poisson_solve.PoissonSolve

Bases: ImplicitDiffusion

Weak discretization of the (stabilized) Poisson equation: find \(\phi \in H^1\) such that

\[\epsilon \int_\Omega \psi\, \phi\,\textrm d \mathbf x + \int_\Omega \nabla \psi^\top \, \nabla \phi \,\textrm d \mathbf x = \sum_i \int_\Omega \psi\, \rho_i(\mathbf x)\,\textrm d \mathbf x \qquad \forall \ \psi \in H^1\,,\]

where \(\epsilon \in \mathbb R\) is a stabilization parameter. Boundary terms from integration by parts are assumed to vanish. The equation is discretized as

\[\left( \epsilon\,\mathbb S + \mathbb G^\top \mathbb M^1 \mathbb G \right)\, \boldsymbol{\phi} = \sum_i(\Lambda^0, \rho_i )_{L^2}\,,\]

where \(\mathbb M^1\) is the \(H(\textnormal{curl})\)-mass matrix and \(\mathbb S\) is a stabilization matrix.

class Options(stab_eps: float = 0.0, stab_mat: Literal['M0', 'M0ad', 'Id'] = 'Id', diffusion_mat: Literal['M1', 'M1perp', 'M1para', 'M1gyro'] = 'M1', rho: FEECVariable | Callable | tuple[AccumulatorVector, Particles] | list | None = None, rho_coeffs: float | list | None = None, x0: StencilVector | None = None, solver: Literal['pcg', 'cg'] = 'pcg', precond: Literal['MassMatrixPreconditioner', 'MassMatrixDiagonalPreconditioner', None] = 'MassMatrixPreconditioner', solver_params: SolverParameters | None = None)

Bases: OptionsBase

Configuration options for Poisson.

Parameters:

• stab_eps (float, default=0.0) – Stabilization weight used for the mass-like term in the Poisson operator. Internally mapped to sigma_1 = stab_eps in the parent ImplicitDiffusion formulation.

• stab_mat ({"M0", "M0ad", "Id"}, default="Id") –

  Stabilization matrix multiplied by stab_eps.

  • "M0": standard weighted 0-form mass operator.

  • "M0ad": adiabatic-electron weighted 0-form mass operator.

  • "Id": identity operator.

• diffusion_mat ({"M1", "M1perp", "M1gyro"}, defaults="M1") –

  Diffusion matrix.

  • "M1": standard weighted 1-form mass operator.

  • "M1perp": weighted 1-form mass operator perpendicular to magnetic field.

  • "M1para": weighted 1-form mass operator parallele to magnetic field.

  • "M1gyro": weighted 1-form mass operator used in gyrokinetic model.

• rho (FEECVariable or Callable or tuple or list, default=None) –

  Right-hand side source term(s) of the Poisson problem. Accepted entries are:

  • None: zero source.

  • FEECVariable in H1.

  • Callable to be projected to H1 via L2Projector.

  • AccumulatorVector.

  • a list containing any mix of the entries above.

  The tuple form is accepted by typing for compatibility with other propagator interfaces that pair particle data with accumulators.

• rho_coeffs (float or list, default=None) – Multiplicative coefficient(s) applied to rho. If None, coefficients default to 1.0 for all sources.

• x0 (StencilVector, default=None) – Initial guess for the iterative linear solver.

• solver (LiteralOptions.OptsSymmSolver, default="pcg") – Name of the symmetric iterative solver passed to psydac.linalg.solvers.inverse().

• precond (LiteralOptions.OptsMassPrecond, default="MassMatrixPreconditioner") – Name of the preconditioner configuration. Currently this class inherits the same behavior as ImplicitDiffusion, where pc=None is used internally.

• solver_params (SolverParameters, default=None) – Iterative-solver controls (for example tol, maxiter, verbose, info, recycle). If None, defaults to SolverParameters().

Notes

Poisson.Options reuses ImplicitDiffusion internals by enforcing sigma_2 = 0.0, sigma_3 = 1.0, divide_by_dt = False and diffusion_mat = "M1" in __post_init__.

class struphy.propagators.curl_curl_solve.CurlCurlSolve

Bases: Propagator

Weak discretization of the curl-curl problem.

Find \(\mathbf E \in H(\textnormal{curl})\) such that

\[\int_\Omega \nabla \times \mathbf F \cdot \nabla \times \mathbf E\,\textrm d \mathbf x - \sigma \int_\Omega \mathbf F \cdot \mathbf E\,\textrm d \mathbf x = \sum_i \int_\Omega \mathbf F \cdot \mathbf J _i\,\textrm d \mathbf x \qquad \forall \,\mathbf F \in H(\textnormal{curl})\,,\]

where \(\mathbf J _i:\Omega \to \mathbb R^3\) are real-valued and \(\sigma \in \mathbb R / \{0\}\) is a scalar. Boundary terms from integration by parts are assumed to vanish. The equation is discretized as

\[\left( \mathbb C^\top \cdot \mathbb M^2 \cdot \mathbb C - \sigma \mathbb M^1 \right)\, \boldsymbol e^{n+1} =\sum_i \boldsymbol j_i\,,\]

where \(\mathbb M^1\) and \(\mathbb M^2\) are WeightedMassOperators and \(\boldsymbol j_i\) is the vector of coefficients of the projection of \(\mathbf J_i\) into the discrete space \(V^1_h\subset H(\textnormal{curl})\).

class Variables

Bases: object

Variables:

e (FEECVariable) – Vector-valued solution in "Hcurl" space.

class Options(sigma: float = 1.0, stab_mat: Literal['M1', 'Id'] = 'M1', diffusion_mat: Literal['M2'] = 'M2', j: FEECVariable | tuple[Callable, Callable, Callable] | tuple[AccumulatorVector, Particles] | list | None = None, j_coeffs: float | list | None = None, x0: StencilVector | None = None, solver: Literal['pcg', 'cg'] = 'pcg', precond: Literal['MassMatrixPreconditioner', 'MassMatrixDiagonalPreconditioner', None] = 'MassMatrixPreconditioner', solver_params: SolverParameters | None = None)

Bases: object

Configuration options for CurlCurlSolve.

Parameters:

• sigma (float, default=1.0) – Coefficient multiplying the stabilization/mass contribution on the left-hand side.

• stab_mat ({"M1", "Id"}, default="M1") –

  Stabilization operator used in the term weighted by sigma.

  • "M1": standard weighted 1-form mass operator.

  • "Id": identity operator.

• diffusion_mat ("M2" or WeightedMassOperator, default="M2") – Diffusion metric in the bilinear form curl.T @ diffusion_mat @ curl. You can pass the name of a pre-built operator in mass_ops or a custom WeightedMassOperator compatible with the codomain of curl.

• j (FEECVariable or list of Callables or tuple or list, default=None) –

  Source term(s) on the right-hand side. Accepted entries are:

  • None: zero source.

  • FEECVariable in "Hcurl".

  • list of Callable``s to be projected to ``"Hcurl" via L2Projector.

  • AccumulatorVector.

  • a list containing any mix of the entries above.

  The tuple form is accepted by typing for compatibility with other propagator interfaces that pair particle data with accumulators.

• j_coeffs (float or list, default=None) – Multiplicative coefficient(s) for j sources. If a scalar is provided, it is applied to a single source. If a sequence is provided, its length must match the number of collected sources. If None, all coefficients default to 1.0.

• x0 (StencilVector, default=None) – Initial guess for the iterative linear solver.

• solver (LiteralOptions.OptsSymmSolver, default="pcg") – Name of the symmetric iterative solver passed to psydac.linalg.solvers.inverse().

• precond (LiteralOptions.OptsMassPrecond, default="MassMatrixPreconditioner") – Name of the preconditioner configuration. Currently this class sets pc=None internally, so this option is reserved for compatibility and future extensions.

• solver_params (SolverParameters, default=None) – Iterative-solver controls (for example tol, maxiter, verbose, info, recycle). If None, defaults to SolverParameters().

property sources: list[StencilVector | FEECVariable | AccumulatorVector]

Right-hand side of the equation (sources).

property coeffs: list[float]

Same length as self.sources. Coefficients multiplied with sources before solve (default is 1.0).

property x0

feectools.linalg.stencil.StencilVector or struphy.polar.basic.PolarVector. First guess of the iterative solver.

class struphy.propagators.adiabatic_phi.AdiabaticPhi(phi: StencilVector, *, A_mat: WeightedMassOperator = 'M0', rho: StencilVector | tuple | None = None, sigma_1: float = 1.0, x0: StencilVector | None = None, **params)

Bases: Propagator

Electrostatic potential for adiabatic electrons, computed from

\[n_e = n_{e0}\,\exp \left( \frac{e \phi}{k_B T_e} \right) \approx n_{e0} \left( 1 + \frac{e \phi}{k_B T_{e0}} \right)\,,\]

where \(n_{e0}\) and \(T_{e0}\) denote electron equilibrium density and temperature, respectively.

This is solved in weak form: find \(\phi \in H^1\) such that

\[\int_\Omega \psi\, \frac{n_{e0}(\mathbf x)}{T_{e0}(\mathbf x)}\phi \,\textrm d \mathbf x = \int_\Omega \psi\, (n_{e}(\mathbf x) - n_{e0}(\mathbf x))\,\textrm d \mathbf x \qquad \forall \ \psi \in H^1\,.\]

The equation is discretized as

\[\sigma_1 \mathbb M^0_{n/T} \boldsymbol \phi = (\Lambda^0, n_{e} - n_{e0} )_{L^2}\,,\]

where \(M^0_{n/T}\) is a WeightedMassOperator and \(\sigma_1\) is a normalization parameter.

Parameters:

• phi (StencilVector) – FE coefficients of the solution as a discrete 0-form.

• A_mat (WeightedMassOperator) – The matrix to invert.

• rho (StencilVector or tuple) – Right-hand side FE coefficients of a 0-form (optional, can be set with a setter later). Can be either a) StencilVector or b) 2-tuple. In case b) the first tuple entry must be AccumulatorVector, and the second entry must be Particles.

• sigma_1 (float) – Normalization parameter.

• x0 (StencilVector) – Initial guess for the iterative solver (optional, can be set with a setter later).

• **params (dict) – Parameters for the iterative solver (see __init__ for details).

property rho

Right-hand side FE coefficients of a 0-form. Can be either a) StencilVector or b) 2-tuple. In the latter case, the first tuple entry must be AccumulatorVector, and the second entry must be Particles.

property x0

feectools.linalg.stencil.StencilVector or struphy.polar.basic.PolarVector. First guess of the iterative solver.

Linear MHD solvers

class struphy.propagators.shear_alfven_propagator.ShearAlfvenPropagator

Bases: Propagator

FEEC discretization of the following equations: find \(\tilde{\mathbf{U}} \in \{H(\textnormal{curl}), H(\textnormal{div}), (H^1)^3\}\) and \(\tilde{\mathbf{B}} \in H(\textnormal{div})\) such that

\[\begin{split}\int_{\Omega} \rho_0\frac{\partial \tilde{\mathbf{U}}}{\partial t} \cdot \mathbf V\,\textnormal d \mathbf x &= \int_{\Omega} \tilde{\mathbf{B}} \cdot \nabla \times (\mathbf{B}_0 \times \mathbf V) \,\textnormal d \mathbf x \qquad \forall \,\mathbf V \in \{H(\textnormal{curl}), H(\textnormal{div}), (H^1)^3\} \\[2mm] \frac{\partial \tilde{\mathbf{B}}}{\partial t} &= \nabla\times(\tilde{\mathbf{U}} \times \mathbf{B}_0)\end{split}\]

Time discretization:

• implicit: Crank-Nicolson (implicit mid-point)

• explicit: explicit RK methods from ButcherTableau

System size reduction via SchurSolver.

class Variables

Bases: object

Container for variables advanced by ShearAlfvenPropagator.

Variables:

• u (FEECVariable) – Velocity variable in one of "Hcurl", "Hdiv", or "H1vec".

• b (FEECVariable) – Magnetic-field variable in "Hdiv" space.

class Options(u_space: Literal['Hcurl', 'Hdiv', 'H1vec'] = 'Hdiv', algo: Literal['implicit', 'explicit'] = 'implicit', solver: Literal['pcg', 'cg'] = 'pcg', precond: Literal['MassMatrixPreconditioner', 'MassMatrixDiagonalPreconditioner', None] = 'MassMatrixPreconditioner', solver_params: SolverParameters | None = None, butcher: ButcherTableau | None = None)

Bases: OptionsBase

Configuration options for ShearAlfvenPropagator.

Parameters:

• u_space (LiteralOptions.OptsVecSpace, default="Hdiv") – FEEC space used for the velocity variable.

• algo ({"implicit", "explicit"}, default="implicit") – Time stepping algorithm.

• solver (LiteralOptions.OptsSymmSolver, default="pcg") – Symmetric iterative solver used by implicit or explicit operators.

• precond (LiteralOptions.OptsMassPrecond, default="MassMatrixPreconditioner") – Preconditioner for the mass-matrix block.

• solver_params (SolverParameters, default=None) – Solver controls; defaults to SolverParameters().

• butcher (ButcherTableau, default=None) – Explicit Runge-Kutta tableau when algo="explicit".

class struphy.propagators.shear_alfven_b1.ShearAlfvenB1

Bases: Propagator

FEEC discretization of the following equations: find \(\mathbf U \in \{H(\textnormal{curl}), H(\textnormal{div}), (H^1)^3\}\) and \(\mathbf B \in H(\textnormal{curl})\) such that

\[\begin{split}&\int_\Omega \rho_0\frac{\partial \tilde{\mathbf{U}}}{\partial t} \cdot \mathbf V\,\textnormal d \mathbf x =\int_\Omega \tilde{\mathbf B } \cdot \nabla \times (\mathbf{B}_0 \times \mathbf V) \,\textnormal d \mathbf x \qquad \forall \,\mathbf V \in \{H(\textnormal{curl}), H(\textnormal{div}), (H^1)^3\}\,, \\[2mm] &\int_\Omega \frac{\partial \tilde{\mathbf{B}}}{\partial t} \cdot \mathbf C\,\textnormal d \mathbf x - \int_\Omega \mathbf C \cdot\nabla\times \left( \tilde{\mathbf{U}} \times \mathbf{B}_0 \right) \textrm d \mathbf x = 0 \qquad \forall \, \mathbf C \in H(\textrm{curl})\,.\end{split}\]

Time discretization: Crank-Nicolson (implicit mid-point). System size reduction via SchurSolver:

\[\begin{split}\begin{bmatrix} \mathbf u^{n+1} - \mathbf u^n \\ \mathbf b^{n+1} - \mathbf b^n \end{bmatrix} = \frac{\Delta t}{2} \begin{bmatrix} 0 & (\mathbb M^\rho_2)^{-1} \mathcal {T^2}^\top \mathbb C \mathbb M_1^{-1}\\ - \mathbb M_1^{-1} \mathbb C^\top \mathcal {T^2} (\mathbb M^\rho_2)^{-1} & 0 \end{bmatrix} \begin{bmatrix} {\mathbb M^\rho_2}(\mathbf u^{n+1} + \mathbf u^n) \\ \mathbb M_1(\mathbf b^{n+1} + \mathbf b^n) \end{bmatrix} ,\end{split}\]

where \(\mathbb M^\rho_2\) is a weighted mass matrix in 2-space, the weight being \(\rho_0\), the MHD equilibirum density.

class Variables

Bases: object

Container for variables advanced by ShearAlfvenB1.

Variables:

• u (FEECVariable) – Velocity variable in "Hdiv" space.

• b (FEECVariable) – Magnetic-field variable in "Hcurl" space.

class Options(solver: Literal['pcg', 'cg'] = 'pcg', precond: Literal['MassMatrixPreconditioner', 'MassMatrixDiagonalPreconditioner', None] = 'MassMatrixPreconditioner', solver_params: SolverParameters | None = None, solver_M1: Literal['pcg', 'cg'] = 'pcg', precond_M1: Literal['MassMatrixPreconditioner', 'MassMatrixDiagonalPreconditioner', None] = 'MassMatrixPreconditioner', solver_params_M1: SolverParameters | None = None)

Bases: OptionsBase

Configuration options for ShearAlfvenB1.

Parameters:

• solver (LiteralOptions.OptsSymmSolver, default="pcg") – Solver used for the Schur-reduced system.

• precond (LiteralOptions.OptsMassPrecond, default="MassMatrixPreconditioner") – Preconditioner for M2n block.

• solver_params (SolverParameters, default=None) – Solver controls for the Schur solve.

• solver_M1 (LiteralOptions.OptsSymmSolver, default="pcg") – Solver used for inverting M1.

• precond_M1 (LiteralOptions.OptsMassPrecond, default="MassMatrixPreconditioner") – Preconditioner for M1 inversion.

• solver_params_M1 (SolverParameters, default=None) – Solver controls for M1 inversion.

class struphy.propagators.magnetosonic.Magnetosonic

Bases: Propagator

FEEC discretization of the following equations: find \(\tilde{\rho} \in L^2, \tilde{\mathbf U} \in \{H(\textnormal{curl}), H(\textnormal{div}), (H^1)^3\}, \tilde p \in L^2\) such that

\[\begin{split}\frac{\partial \tilde{\rho}}{\partial t} + \nabla\cdot(\rho_0 \tilde{\mathbf{U}}) &= 0 \\[2mm] \int_{\Omega} \rho_0\frac{\partial \tilde{\mathbf{U}}}{\partial t} \cdot \mathbf V\,\textrm{d}\mathbf x - \int_{\Omega} \tilde p\,\nabla\cdot\mathbf V\,\textrm{d}\mathbf x &= \int_{\Omega} (\nabla\times\mathbf{B}_0)\times\tilde{\mathbf{B}}\cdot\mathbf V\,\textrm{d}\mathbf x \qquad \forall\,\mathbf V \in \{H(\textnormal{curl}), H(\textnormal{div}), (H^1)^3\} \\[2mm] \frac{\partial \tilde p}{\partial t} + \nabla\cdot(p_0 \tilde{\mathbf{U}}) + \frac{2}{3}\,p_0\nabla\cdot\tilde{\mathbf{U}} &= 0\end{split}\]

Time discretization: Crank-Nicolson (implicit mid-point).

System size reduction via SchurSolver:

class Variables

Bases: object

Container for variables advanced by Magnetosonic.

Variables:

• n (FEECVariable) – Density perturbation variable in "L2" space.

• u (FEECVariable) – Velocity perturbation variable in one of "Hcurl", "Hdiv", or "H1vec".

• p (FEECVariable) – Pressure perturbation variable in "L2" space.

class Options(b_field: FEECVariable | None = None, u_space: Literal['Hcurl', 'Hdiv', 'H1vec'] = 'Hdiv', solver: Literal['pbicgstab', 'bicgstab', 'gmres'] = 'pbicgstab', precond: Literal['MassMatrixPreconditioner', 'MassMatrixDiagonalPreconditioner', None] = 'MassMatrixPreconditioner', solver_params: SolverParameters | None = None)

Bases: OptionsBase

Configuration options for Magnetosonic.

Parameters:

• b_field (FEECVariable, default=None) – Background magnetic-field variable in "Hdiv" space used in the Lorentz-force coupling. If None, a default FEECVariable(space="Hdiv") is created.

• u_space (LiteralOptions.OptsVecSpace, default="Hdiv") – FEEC space used for the velocity variable u.

• solver (LiteralOptions.OptsGenSolver, default="pbicgstab") – Iterative solver used by SchurSolver.

• precond (LiteralOptions.OptsMassPrecond, default="MassMatrixPreconditioner") – Preconditioner applied to the mass matrix block associated with u_space.

• solver_params (SolverParameters, default=None) – Iterative-solver controls. If None, defaults to SolverParameters().

class struphy.propagators.magnetosonic_uniform.MagnetosonicUniform

Bases: Propagator

FEEC discretization of the following equations: find \(\tilde \rho \in L^2, \tilde{\mathbf U} \in \{H(\textnormal{curl}), H(\textnormal{div}), (H^1)^3\}, \tilde p \in L^2\) such that

\[ \begin{align}\begin{aligned}&\frac{\partial \tilde \rho}{\partial t}+\nabla\cdot(\rho_0 \tilde{\mathbf{U}})=0\,,\\\int \rho_0&\frac{\partial \tilde{\mathbf{U}}}{\partial t} \cdot \mathbf V\,\textrm d \mathbf x - \int \tilde p\, \nabla \cdot \mathbf V \,\textrm d \mathbf x = 0 \qquad \forall \ \mathbf V \in \{H(\textnormal{curl}), H(\textnormal{div}), (H^1)^3\}\,,\\&\frac{\partial \tilde p}{\partial t} + \frac{5}{3}\,p_0\nabla\cdot \tilde{\mathbf{U}}=0\,.\end{aligned}\end{align} \]

Time discretization: Crank-Nicolson (implicit mid-point). System size reduction via SchurSolver:

\[\begin{split}\begin{bmatrix} \mathbf u^{n+1} - \mathbf u^n \\ \mathbf p^{n+1}_i - \mathbf p^n_i \end{bmatrix} = \frac{\Delta t}{2} \begin{bmatrix} 0 & (\mathbb M^\rho_2)^{-1} \mathbb D^\top \mathbb M_3 \\ - \gamma \mathcal K^3 \mathbb D & 0 \end{bmatrix} \begin{bmatrix} (\mathbf u^{n+1} + \mathbf u^n) \\ (\mathbf p^{n+1}_i + \mathbf p^n_i) \end{bmatrix} ,\end{split}\]

where \(\mathbb M^\rho_2\) is a weighted mass matrix in 2-space, the weight being the MHD equilibirum density \(\rho_0\). Furthermore, \(\mathcal K^3\) is the basis projection operator given by :

\[\mathcal{K}^3_{ijk,mno} := \hat{\Pi}^3_{ijk} \left[ \frac{\hat{p}^3_{\text{eq}}}{\sqrt{g}}\Lambda^3_{mno} \right] \,.\]

The solution of the above system is based on the Schur complement.

Decoupled density update:

\[\boldsymbol{\rho}^{n+1} = \boldsymbol{\rho}^n - \frac{\Delta t}{2} \mathcal Q \mathbb D (\mathbf u^{n+1} + \mathbf u^n) \,.\]

Parameters:

• n (feectools.linalg.stencil.StencilVector) – FE coefficients of a discrete 3-form.

• u (feectools.linalg.block.BlockVector) – FE coefficients of MHD velocity 2-form.

• p (feectools.linalg.stencil.StencilVector) –

  FE coefficients of a discrete 3-form.

  **paramsdict

  Solver- and/or other parameters for this splitting step.

class Variables

Bases: object

Container for variables advanced by MagnetosonicUniform.

Variables:

• n (FEECVariable) – Density perturbation variable in "L2" space.

• u (FEECVariable) – Velocity perturbation variable in one of "Hcurl", "Hdiv", or "H1vec".

• p (FEECVariable) – Pressure perturbation variable in "L2" space.

class Options(solver: Literal['pbicgstab', 'bicgstab', 'gmres'] = 'pbicgstab', precond: Literal['MassMatrixPreconditioner', 'MassMatrixDiagonalPreconditioner', None] = 'MassMatrixPreconditioner', solver_params: SolverParameters | None = None)

Bases: OptionsBase

Configuration options for MagnetosonicUniform.

Parameters:

• solver (LiteralOptions.OptsGenSolver, default="pbicgstab") – Iterative solver used by SchurSolver.

• precond (LiteralOptions.OptsMassPrecond, default="MassMatrixPreconditioner") – Preconditioner applied to the M2n block.

• solver_params (SolverParameters, default=None) – Iterative-solver controls. If None, defaults to SolverParameters().

class struphy.propagators.faraday_extended.FaradayExtended(a, **params)

Bases: Propagator

Equations: Faraday’s law

\[\begin{split}\begin{align*} & \frac{\partial {\mathbf A}}{\partial t} = - \frac{\nabla \times (\nabla \times {\mathbf A} + {\mathbf B}_0) }{n} \times (\nabla \times {\mathbf A} + {\mathbf B}_0) - \frac{\int ({\mathbf A} - {\mathbf p}f \mathrm{d}{\mathbf p})}{n} \times (\nabla \times {\mathbf A} + {\mathbf B}_0), \\ & n = \int f \mathrm{d}{\mathbf p}. \end{align*}\end{split}\]

Mid-point rule:

\[\begin{split}\begin{align*} & \left[ \mathbb{M}_1 - \frac{\Delta t}{2} \mathbb{F}(\hat{n}^0_h, {\mathbf a}^{n+\frac{1}{2}}) \mathbb{M}_1^{-1} (\mathbb{P}_1^\top \mathbb{W} \mathbb{P}_1 + \mathbb{C}^\top \mathbb{M}_2 \mathbb{C} ) \right] {\mathbf a}^{n+1} \\ & = \mathbb{M}_1 {\mathbf a}^n + \frac{\Delta t}{2} \mathbb{F}(\hat{n}^0_h, {\mathbf a}^{n+\frac{1}{2}}) \mathbb{M}_1^{-1} (\mathbb{P}_1^\top \mathbb{W} \mathbb{P}_1 + \mathbb{C}^\top \mathbb{M}_2 \mathbb{C} ) {\mathbf a}^{n+1} \\ & - \Delta t \mathbb{F}(\hat{n}^0_h, {\mathbf a}^{n+\frac{1}{2}}) \mathbb{M}_1^{-1} \mathbb{P}_1^\top \mathbb{W} {\mathbf P}^n\\ & + \Delta t \mathbb{F}(\hat{n}^0_h, {\mathbf a}^{n+\frac{1}{2}}) \mathbb{M}_1^{-1} \mathbb{C}^\top \mathbb{M}_2 {\mathbf b}_0\\ & \mathbb{F}_{ij} = - \int \frac{1}{\hat{n}^0_h \sqrt{g}} G (\nabla \times {\mathbf A} + {\mathbf B}_0) \cdot (\Lambda^1_i \times \Lambda^1_j) \mathrm{d}{\boldsymbol \eta}. \end{align*}\end{split}\]

Parameters:

• a (feectools.linalg.block.BlockVector) – FE coefficients of vector potential.

• **params (dict) – Solver- and/or other parameters for this splitting step.

class struphy.propagators.hall.Hall

Bases: Propagator

FEEC discretization of the following equations: find \(\mathbf B \in H(\textnormal{curl})\) such that

\[\int_\Omega \frac{\partial \tilde{\mathbf{B}}}{\partial t} \cdot \mathbf C\,\textnormal d \mathbf x + \frac{1}{\varepsilon} \int_\Omega \nabla\times \mathbf C \cdot \left( \frac{\nabla\times \tilde{\mathbf{B}}}{\rho_0}\times \mathbf{B}_0 \right) \textrm d \mathbf x = 0 \qquad \forall \, \mathbf C \in H(\textrm{curl})\]

Time discretization: Crank-Nicolson (implicit mid-point):

\[\mathbf b^{n+1} - \mathbf b^n = \frac{\Delta t}{2} \mathbb M_1^{-1} \mathbb C^\top \mathbb M^{\mathcal{T},\rho}_2 \mathbb C (\mathbf b^{n+1} + \mathbf b^n) ,\]

where \(\mathbb M^{\mathcal{T},\rho}_2\) is a weighted mass matrix in 2-space, the weight being \(\frac{\mathcal{T}}{\rho_0}\), the MHD equilibirum density \(\rho_0\) as a 0-form, and rotation matrix \(\mathcal{T} \vec v = \vec B^2_{\textnormal{eq}} \times \vec v\,,\). The solution of the above system is based on the Pre-conditioned Biconjugate Gradient Stabilized algortihm (PBiConjugateGradientStab).

class Variables

Bases: object

Container for variables advanced by Hall.

Variables:

b (FEECVariable) – Magnetic-field variable in "Hcurl" space.

class Options(solver: Literal['pbicgstab', 'bicgstab', 'gmres'] = 'pbicgstab', precond: Literal['MassMatrixPreconditioner', 'MassMatrixDiagonalPreconditioner', None] = 'MassMatrixPreconditioner', solver_params: SolverParameters | None = None, epsilon_from: Species | None = None)

Bases: OptionsBase

Configuration options for Hall.

Parameters:

• solver (LiteralOptions.OptsGenSolver, default="pbicgstab") – General iterative solver used for the linear Hall update.

• precond (LiteralOptions.OptsMassPrecond, default="MassMatrixPreconditioner") – Preconditioner applied to the M1 block.

• solver_params (SolverParameters, default=None) – Iterative-solver controls. If None, defaults to SolverParameters().

• epsilon_from (Species, default=None) – Species object from which epsilon is taken. If None, epsilon defaults to 1.0.

Nonlinear MHD solvers

class struphy.propagators.variational_density_evolve.VariationalDensityEvolve

Bases: Propagator

FEEC discretization of the following equations: find \(\rho \in L^2\) and \(\mathbf u \in (H^1)^3\) such that

\[\begin{split}&\partial_t \rho + \nabla \cdot ( \tilde{\rho} \mathbf u ) = 0 \,, \\[4mm] &\int_\Omega \partial_t (\rho \mathbf u) \cdot \mathbf v\,\textrm d \mathbf x + \int_\Omega \left(\frac{|\mathbf u|^2}{2} - \frac{\partial(\rho \mathcal U(\rho))}{\partial \rho}\right) \nabla \cdot (\tilde{\rho} \mathbf v) \,\textrm d \mathbf x = 0 \qquad \forall \, \mathbf v \in (H^1)^3\,.\end{split}\]

Where \(\tilde{\rho}\) is either \(\rho\) for full-f models, \(\rho_0\) for linear models or \(\rho_0+\rho\) for \(\delta f\) models.

In the case of linear model, the second equation is not updated.

On the logical domain:

\[\begin{split}\begin{align} &\partial_t \hat{\rho}^3 + \nabla \cdot ( \hat{\rho}^3 \hat{\mathbf{u}} ) = 0 ~ , \\[4mm] &\int_{\hat{\Omega}} \partial_t ( \hat{\rho}^3 \hat{\mathbf{u}}) \cdot G \hat{\mathbf{v}} \, \textrm d \boldsymbol \eta + \int_{\hat{\Omega}} \left( \frac{| DF \hat{\mathbf{u}} |^2}{2} - \frac{\partial (\hat{\rho}^3 \mathcal U)}{\partial \hat{\rho}^3} \right) \nabla \cdot (\hat{\rho}^3 \hat{\mathbf{v}}) \, \textrm d \boldsymbol \eta = 0 ~ , \\[2mm] \end{align}\end{split}\]

where \(\mathcal U\) depends on the chosen model. It is discretized as

\[\begin{split}\begin{align} &\frac{\mathbb M^v[\hat{\rho}_h^{n+1}] \mathbf u^{n+1}- \mathbb M^v[\hat{\rho}_h^n] \mathbf u^n}{\Delta t} + (\mathbb D \hat{\Pi}^{2}[\hat{\rho_h^{n}} \vec{\boldsymbol \Lambda}^v])^\top \hat{l}^3\left(\frac{DF \hat{\mathbf{u}}_h^{n+1} \cdot DF \hat{\mathbf{u}}_h^{n}}{2} - \frac{\hat{\rho}_h^{n+1}\mathcal U(\hat{\rho}_h^{n+1})-\hat{\rho}_h^{n}\mathcal U(\hat{\rho}_h^{n})}{\hat{\rho}_h^{n+1}-\hat{\rho}_h^n} \right) = 0 ~ , \\[2mm] &\frac{\boldsymbol \rho^{n+1}- \boldsymbol \rho^n}{\Delta t} + \mathbb D \hat{\Pi}^{2}[\hat{\rho_h^{n}} \vec{\boldsymbol \Lambda}^v] \mathbf u^{n+1/2} = 0 ~ , \\[2mm] \end{align}\end{split}\]

where \(\hat{l}^3(f)\) denotes the vector representing the linear form \(v_h \mapsto \int_{\hat{\Omega}} f(\boldsymbol \eta) v_h(\boldsymbol \eta) d \boldsymbol \eta\), that is the vector with components

\[\hat{l}^3(f)_{ijk}=\int_{\hat{\Omega}} f \Lambda^3_{ijk} \textrm d \boldsymbol \eta\]

and the weights in the the BasisProjectionOperator and the WeightedMassOperator are given by

\[\hat{\mathbf{u}}_h^{k} = (\mathbf{u}^{k})^\top \vec{\boldsymbol \Lambda}^v \in (V_h^0)^3 \, \text{for k in} \{n, n+1/2, n+1\}, \qquad \hat{\rho}_h^{k} = (\rho^{k})^\top \vec{\boldsymbol \Lambda}^3 \in V_h^3 \, \text{for k in} \{n, n+1/2, n+1\} .\]

class Variables

Bases: object

Container for variables advanced by VariationalDensityEvolve.

Variables:

• rho (FEECVariable) – Density variable in "L2" space.

• u (FEECVariable) – Velocity variable in "H1vec" space.

class Options(model: Literal['pressureless', 'barotropic', 'full', 'full_p', 'full_q', 'linear', 'deltaf', 'linear_q', 'deltaf_q'] = 'barotropic', gamma: float = 1.6666666666666667, solver: Literal['pcg', 'cg'] = 'pcg', precond: Literal['MassMatrixPreconditioner', 'MassMatrixDiagonalPreconditioner', None] = 'MassMatrixPreconditioner', solver_params: SolverParameters | None = None, nonlin_solver: NonlinearSolverParameters | None = None, s: FEECVariable | None = None)

Bases: OptionsBase

Configuration options for VariationalDensityEvolve.

Parameters:

• model ({"pressureless", "barotropic", "full", "full_p", "full_q", "linear", "deltaf", "linear_q", "deltaf_q"}, default="barotropic") – Density/thermodynamic model variant.

• gamma (float, default=5/3) – Adiabatic index.

• solver (LiteralOptions.OptsSymmSolver, default="pcg") – Linear solver for implicit substeps.

• precond (LiteralOptions.OptsMassPrecond, default="MassMatrixPreconditioner") – Preconditioner used in linear solves.

• solver_params (SolverParameters, default=None) – Linear-solver controls.

• nonlin_solver (NonlinearSolverParameters, default=None) – Nonlinear iteration controls.

• s (FEECVariable, default=None) – Entropy-like variable required by model="full".

class struphy.propagators.variational_entropy_evolve.VariationalEntropyEvolve

Bases: Propagator

FEEC discretization of the following equations: find \(\mathbf u \in (H^1)^3\) and \(s \in L^2\) such that

\[\begin{split}&\int_\Omega \partial_t (\rho \mathbf u) \cdot \mathbf v\,\textrm d \mathbf x - \int_\Omega \frac{\partial(\rho \mathcal U)}{\partial s} \nabla \cdot (s \mathbf v) \,\textrm d \mathbf x = 0 \qquad \forall \, \mathbf v \in (H^1)^3\,, \\[4mm] &\partial_t s + \nabla \cdot ( s \mathbf u ) = 0 \,.\end{split}\]

On the logical domain:

\[\begin{split}\begin{align} &\int_{\hat{\Omega}} \partial_t ( \hat{\rho}^3 \hat{\mathbf{u}}) \cdot G \hat{\mathbf{v}} \, \textrm d \boldsymbol \eta - \int_{\hat{\Omega}} \left(\frac{\partial \hat{\rho}^3 \mathcal U}{\partial \hat{s}} \right) \nabla \cdot (\hat{s} \hat{\mathbf{v}}) \, \textrm d \boldsymbol \eta = 0 ~ , \\[2mm] &\partial_t \hat{s} + \nabla \cdot ( \hat{s} \hat{\mathbf{u}} ) = 0 ~ , \end{align}\end{split}\]

where \(\mathcal U\) depends on the chosen model. It is discretized as

\[\begin{split}\begin{align} &\mathbb M^v[\hat{\rho}_h^{n}] \frac{ \mathbf u^{n+1}-\mathbf u^n}{\Delta t} - (\mathbb D \hat{\Pi}^{2}[\hat{s_h^{n}} \vec{\boldsymbol \Lambda}^v])^\top \hat{l}^3\left( \frac{\hat{\rho}_h^{n}\mathcal U(\hat{\rho}_h^{n},\hat{s}_h^{n+1})-\hat{\rho}_h^{n}\mathcal U(\hat{\rho}_h^{n},\hat{s}_h^{n})}{\hat{s}_h^{n+1}-\hat{s}_h^n} \right) = 0 ~ , \\[2mm] &\frac{\mathbf s^{n+1}- \mathbf s^n}{\Delta t} + \mathbb D \hat{\Pi}^{2}[\hat{s_h^{n}} \vec{\boldsymbol \Lambda}^v] \mathbf u^{n+1/2} = 0 ~ , \\[2mm] \end{align}\end{split}\]

where \(\hat{l}^3(f)\) denotes the vector representing the linear form \(v_h \mapsto \int_{\hat{\Omega}} f(\boldsymbol \eta) v_h(\boldsymbol \eta) d \boldsymbol \eta\), that is the vector with components

\[\hat{l}^3(f)_{ijk}=\int_{\hat{\Omega}} f \Lambda^3_{ijk} \textrm d \boldsymbol \eta\]

and the weights in the the BasisProjectionOperator and the WeightedMassOperator are given by

\[\hat{\mathbf{u}}_h^{k} = (\mathbf{u}^{k})^\top \vec{\boldsymbol \Lambda}^v \in (V_h^0)^3 \, \text{for k in} \{n, n+1/2, n+1\}, \qquad \hat{s}_h^{k} = (s^{k})^\top \vec{\boldsymbol \Lambda}^3 \in V_h^3 \, \text{for k in} \{n, n+1/2, n+1\} \qquad \hat{\rho}_h^{n} = (\rho^{n})^\top \vec{\boldsymbol \Lambda}^3 \in V_h^3 \.\]

class Variables

Bases: object

Container for variables advanced by VariationalEntropyEvolve.

Variables:

• s (FEECVariable) – Entropy variable in "L2" space.

• u (FEECVariable) – Velocity variable in "H1vec" space.

class Options(model: Literal['full'] = 'full', gamma: float = 1.6666666666666667, solver: Literal['pcg', 'cg'] = 'pcg', precond: Literal['MassMatrixPreconditioner', 'MassMatrixDiagonalPreconditioner', None] = 'MassMatrixPreconditioner', solver_params: SolverParameters | None = None, nonlin_solver: NonlinearSolverParameters | None = None, rho: FEECVariable | None = None)

Bases: OptionsBase

Configuration options for VariationalEntropyEvolve.

Parameters:

• model ({"full"}, default="full") – Entropy-evolution model selection.

• gamma (float, default=5/3) – Adiabatic index.

• solver (LiteralOptions.OptsSymmSolver, default="pcg") – Linear solver for implicit substeps.

• precond (LiteralOptions.OptsMassPrecond, default="MassMatrixPreconditioner") – Preconditioner used in linear solves.

• solver_params (SolverParameters, default=None) – Linear-solver controls.

• nonlin_solver (NonlinearSolverParameters, default=None) – Nonlinear iteration controls.

• rho (FEECVariable, default=None) – Density variable used in variational forms.

class struphy.propagators.variational_mag_field_evolve.VariationalMagFieldEvolve

Bases: Propagator

FEEC discretization of the following equations: find \(\mathbf u \in (H^1)^3\) and \(\mathbf B \in H(\textnormal{div})\) such that

\[\begin{split}&\int_\Omega \partial_t (\rho \mathbf u) \cdot \mathbf v\,\textrm d \mathbf x - \int_\Omega \mathbf B \cdot \nabla \times (\mathbf \tilde{B} \times \mathbf v) \,\textrm d \mathbf x = 0 \qquad \forall \, \mathbf v \in (H^1)^3\,, \\[4mm] &\partial_t \mathbf B + \nabla \cdot ( \mathbf \tilde{B} \times \mathbf u ) = 0 \,.\end{split}\]

Where \(\tilde{\mathbf B}\) is either \(\mathbf B\) for full-f models, \(\mathbf B_0\) for linear models or \(\mathbf B_0+\mathbf B\) for \(\delta f\) models.

On the logical domain:

\[\begin{split}\begin{align} &\int_{\hat{\Omega}} \partial_t ( \hat{\rho}^3 \hat{\mathbf{u}}) \cdot G \hat{\mathbf{v}} \, \textrm d \boldsymbol \eta - \int_{\hat{\Omega}} \hat{\mathbf{B}}^2 \cdot G \,\nabla \times (\hat{\mathbf{B}}^2 \times \hat{\mathbf{v}}) \,\frac{1}{\sqrt g}\, \textrm d \boldsymbol \eta = 0 ~ , \\[2mm] &\partial_t \hat{\mathbf{B}}^2 + \nabla \times (\hat{\mathbf{B}}^2 \times \hat{\mathbf{u}}) = 0 ~ . \end{align}\end{split}\]

It is discretized as

\[\begin{split}\begin{align} &\mathbb M^v[\hat{\rho}_h^{n}] \frac{ \mathbf u^{n+1}-\mathbf u^n}{\Delta t} - (\mathbb C \hat{\Pi}^{1}[B_h^{n+1}} \cdot \vec{\boldsymbol \Lambda}^v])^\top \mathbb M^2 B^{n+\frac{1}{2}} \big) = 0 ~ , \\[2mm] &\frac{\mathbf b^{n+1}- \mathbf b^n}{\Delta t} + \mathbb C \hat{\Pi}^{1}[\hat{B_h^{n}} \cdot \vec{\boldsymbol \Lambda}^v]] \mathbf u^{n+1/2} = 0 ~ , \end{align}\end{split}\]

where weights in the the BasisProjectionOperator and the WeightedMassOperator are given by

\[\hat{\mathbf{B}}_h^{n+1/2} = (\mathbf{b}^{n+1/2})^\top \vec{\boldsymbol \Lambda}^2 \in V_h^2 \, \qquad \hat{\rho}_h^{n} = (\boldsymbol \rho^{n})^\top \vec{\boldsymbol \Lambda}^3 \in V_h^3 \,.\]

class Variables

Bases: object

Container for variables advanced by VariationalMagFieldEvolve.

Variables:

• u (FEECVariable) – Velocity variable in "H1vec" space.

• b (FEECVariable) – Magnetic-field variable in "Hdiv" space.

class Options(model: Literal['full', 'full_p', 'linear'] = 'full', solver: Literal['pcg', 'cg'] = 'pcg', precond: Literal['MassMatrixPreconditioner', 'MassMatrixDiagonalPreconditioner', None] = 'MassMatrixPreconditioner', solver_params: SolverParameters | None = None, nonlin_solver: NonlinearSolverParameters | None = None)

Bases: OptionsBase

Configuration options for VariationalMagFieldEvolve.

Parameters:

• model ({"full", "full_p", "linear"}, default="full") – Magnetic-field evolution model variant.

• solver (LiteralOptions.OptsSymmSolver, default="pcg") – Linear solver for implicit substeps.

• precond (LiteralOptions.OptsMassPrecond, default="MassMatrixPreconditioner") – Preconditioner used in linear solves.

• solver_params (SolverParameters, default=None) – Linear-solver controls.

• nonlin_solver (NonlinearSolverParameters, default=None) – Nonlinear iteration controls.

class struphy.propagators.variational_momentum_advection.VariationalMomentumAdvection

Bases: Propagator

FEEC discretization of the following equations: find \(\mathbf u \in (H^1)^3\) such that

\[\int_{\Omega} \partial_t ( \rho \mathbf{u}) \cdot \mathbf{v} \,\textrm d \mathbf x - \int_{\Omega} \rho \mathbf{u} \cdot [\mathbf{u}, \mathbf{v}] \, \textrm d \mathbf x = 0 \,.\]

On the logical domain:

\[\int_{\hat{\Omega}} \partial_t ( \hat{\rho}^3 \hat{\mathbf{u}}) \cdot G \hat{\mathbf{v}} \,\textrm d \boldsymbol \eta - \int_{\hat{\Omega}} \hat{\rho}^3 \hat{\mathbf{u}} \cdot G [\hat{\mathbf{u}}, \hat{\mathbf{v}}] \, \textrm d \boldsymbol \eta = 0 \,,\]

which is discretized as

\[\mathbb M^v[\hat{\rho}_h^n] \frac{\mathbf u^{n+1}- \mathbf u^n}{\Delta t} - \left(\sum_{\mu} (\hat{\Pi}^{0}[\hat{\mathbf u}_h^{n+1/2} \cdot \vec{\boldsymbol \Lambda}^1] \mathbb G P_{\mu} - \hat{\Pi}^0[\hat{\mathbf A}^1_{\mu,h} \cdot \vec{\boldsymbol \Lambda}^v])^\top P_i \right) \mathbb M^v[\hat{\rho}_h^n] \mathbf u^{n} = 0 ~ .\]

where \(P_\mu\) stand for the CoordinateProjector and the weights in the the two BasisProjectionOperator and the WeightedMassOperator are given by

\[\hat{\mathbf{u}}_h^{n+1/2} = (\mathbf{u}^{n+1/2})^\top \vec{\boldsymbol \Lambda}^v \in (V_h^0)^3 \,, \qquad \hat{\mathbf A}^1_{\mu,h} = \nabla P_\mu((\mathbf u^{n+1/2})^\top \vec{\boldsymbol \Lambda}^v)] \in V_h^1\,, \qquad \hat{\rho}_h^{n} = (\rho^{n})^\top \vec{\boldsymbol \Lambda}^3 \in V_h^3 \,.\]

class Variables

Bases: object

Container for variables advanced by VariationalMomentumAdvection.

Variables:

u (FEECVariable) – Velocity variable in "H1vec" space.

class Options(solver: Literal['pcg', 'cg'] = 'pcg', precond: Literal['MassMatrixPreconditioner', 'MassMatrixDiagonalPreconditioner', None] = 'MassMatrixPreconditioner', solver_params: SolverParameters | None = None, nonlin_solver: NonlinearSolverParameters | None = None)

Bases: OptionsBase

Configuration options for VariationalMomentumAdvection.

Parameters:

• solver (LiteralOptions.OptsSymmSolver, default="pcg") – Linear solver for mass-matrix related solves.

• precond (LiteralOptions.OptsMassPrecond, default="MassMatrixPreconditioner") – Preconditioner used in linear solves.

• solver_params (SolverParameters, default=None) – Linear-solver controls.

• nonlin_solver (NonlinearSolverParameters, default=None) – Nonlinear iteration controls (Picard/Newton).

class struphy.propagators.variational_pb_evolve.VariationalPBEvolve

Bases: Propagator

FEEC discretization of the following equations: find \(\mathbf u \in (H^1)^3\), \(p \in L^2\) and \(\mathbf B \in H(\textnormal{div})\) such that

\[\begin{split}&\int_\Omega \partial_t (\rho \mathbf u) \cdot \mathbf v\,\textrm d \mathbf x - \int_\Omega \mathbf B \cdot \nabla \times (\tilde{\mathbf B} \times \mathbf v) - \int_\Omega \frac{1}{\gamma -1} (\nabla \cdot (\tilde{p} \mathbf v))\,\textrm d \mathbf x = 0 \qquad \forall \, \mathbf v \in (H^1)^3\,, \\[4mm] &\partial_t \mathbf B + \nabla \cdot ( \tilde{\mathbf B} \times \mathbf u ) = 0 \,, \\[4mm] &\partial_t p + \nabla \cdot(\tilde{p} \mathbf u) + (\gamma - 1) \tilde{p} \nabla \cdot u = 0 \,.\end{split}\]

Where \(\tilde{\mathbf B}\) (resp. \(\tilde{p}\)) is either \(\mathbf B\) (resp. \(p\)) for full-f models, \(\mathbf B_0\) (resp. \(p_0\)) for linear models or \(\mathbf B_0+\mathbf B\) (resp. \(p_0+p\)) for \(\delta f\) models.

On the logical domain:

\[\begin{split}\begin{align} &\int_{\hat{\Omega}} \partial_t ( \hat{\rho}^3 \hat{\mathbf{u}}) \cdot G \hat{\mathbf{v}} \, \textrm d \boldsymbol \eta - \int_{\hat{\Omega}} \hat{\mathbf{B}}^2 \cdot G \,\nabla \times (\hat{\mathbf{B}}^2 \times \hat{\mathbf{v}}) \,\frac{1}{\sqrt g}\, - \frac{g}{\gamma -1} \nabla \cdot (\hat{p} \hat{v}) \textrm d \boldsymbol \eta = 0 ~ , \\[2mm] &\partial_t \hat{\mathbf{B}}^2 + \nabla \times (\hat{\mathbf{B}}^2 \times \hat{\mathbf{u}}) = 0 ~ , \\[2mm] &\partial_t \hat{p} + \nabla \cdot (\hat{p} \hat{u}) + (\gamma - 1 ) \hat{p} \nabla \cdot G \hat{u} = 0 \, \end{align}\end{split}\]

It is discretized as

\[\begin{split}\begin{align} &\mathbb M^v[\hat{\rho}_h^{n}] \frac{ \mathbf u^{n+1}-\mathbf u^n}{\Delta t} - (\mathbb C \hat{\Pi}^{1}[\hat{\mathbf B_h^{n+\frac{1}{2}}} \cdot \vec{\boldsymbol \Lambda}^v])^\top \mathbb M^2 \mathbf B^{n+\frac{1}{2}} = 0 ~ , - (\mathbb D \hat{\Pi}^{2}[\hat{p_h^{n+\frac{1}{2}}} \cdot \vec{\boldsymbol \Lambda}^v])^\top \hat{l}^3(\frac{g}{\gamma-1}) \\[2mm] &\frac{\mathbf b^{n+1}- \mathbf b^n}{\Delta t} + \mathbb C \hat{\Pi}^{1}[\hat{\mathbf B_h^{n+\frac{1}{2}}} \cdot \vec{\boldsymbol \Lambda}^v]] \mathbf u^{n+1/2} = 0 ~ , \\[2mm] &\frac{\mathbf p^{n+1}- \mathbf p^n}{\Delta t} + \big(\mathbb D \hat{\Pi}^{2}[\hat{p_h^{n+\frac{1}{2}}} \cdot \vec{\boldsymbol \Lambda}^v]] + (\gamma - 1)\hat{\Pi}^{3}[\hat{p_h^{n+\frac{1}{2}}} \cdot \vec{\boldsymbol \Lambda}^3] \mathbb D \mathcal{U}^v \big) \mathbf u^{n+1/2}= 0 ~ , \\[2mm] \end{align}\end{split}\]

with

\[\hat{l}^3(f)_{ijk}=\int_{\hat{\Omega}} f \Lambda^3_{ijk} \textrm d \boldsymbol \eta\]

where weights in the the BasisProjectionOperator and the WeightedMassOperator are given by

\[\hat{\mathbf{B}}_h^{n+1/2} = (\mathbf{b}^{n+\frac{1}{2}})^\top \vec{\boldsymbol \Lambda}^2 \in V_h^2 \, \qquad \hat{\rho}_h^{n} = (\boldsymbol \rho^{n})^\top \vec{\boldsymbol \Lambda}^3 \in V_h^3 \, \qquad \hat{p}_h^{n+1/2} = (\boldsymbol p^{n+1/2})^\top \vec{\boldsymbol \Lambda}^3 \in V_h^3 \,.\]

and \(\mathcal{U}^v\) is BasisProjectionOperators.

class Variables

Bases: object

Container for variables advanced by VariationalPBEvolve.

Variables:

• p (FEECVariable) – Pressure variable in "L2" space.

• u (FEECVariable) – Velocity variable in "H1vec" space.

• b (FEECVariable) – Magnetic-field variable in "Hdiv" space.

class Options(model: Literal['full_p', 'linear', 'deltaf'] = 'full_p', gamma: float = 1.6666666666666667, solver: Literal['pcg', 'cg'] = 'pcg', precond: Literal['MassMatrixPreconditioner', 'MassMatrixDiagonalPreconditioner', None] = 'MassMatrixPreconditioner', solver_params: SolverParameters | None = None, nonlin_solver: NonlinearSolverParameters | None = None, div_u: FEECVariable | None = None, u2: FEECVariable | None = None, pt3: FEECVariable | None = None, bt2: FEECVariable | None = None)

Bases: OptionsBase

Configuration options for VariationalPBEvolve.

Parameters:

• model ({"full_p", "linear", "deltaf"}, default="full_p") – Pressure-magnetic coupled model variant.

• gamma (float, default=5/3) – Adiabatic index.

• solver (LiteralOptions.OptsSymmSolver, default="pcg") – Linear solver for implicit substeps.

• precond (LiteralOptions.OptsMassPrecond, default="MassMatrixPreconditioner") – Preconditioner used in linear solves.

• solver_params (SolverParameters, default=None) – Linear-solver controls.

• nonlin_solver (NonlinearSolverParameters, default=None) – Nonlinear iteration controls.

• div_u (FEECVariable, default=None) – Optional external divergence-of-velocity field.

• u2 (FEECVariable, default=None) – Optional external velocity in 2-form representation.

• pt3 (FEECVariable, default=None) – Optional pressure background field.

• bt2 (FEECVariable, default=None) – Optional magnetic background field.

class struphy.propagators.variational_qb_evolve.VariationalQBEvolve

Bases: Propagator

FEEC discretization of the following equations: find \(\mathbf u \in (H^1)^3\), \(p \in L^2\) and \(\mathbf B \in H(\textnormal{div})\) such that

\[\begin{split}&\int_\Omega \partial_t (\rho \mathbf u) \cdot \mathbf v\,\textrm d \mathbf x - \int_\Omega \mathbf B \cdot \nabla \times (\tilde{\mathbf B} \times \mathbf v) - \int_\Omega \frac{2 q}{\gamma -1} (\nabla \cdot (\tilde{q} \mathbf v))\,\textrm d \mathbf x = 0 \qquad \forall \, \mathbf v \in (H^1)^3\,, \\[4mm] &\partial_t \mathbf B + \nabla \cdot ( \tilde{\mathbf B} \times \mathbf u ) = 0 \,, \\[4mm] &\partial_t q + \nabla \cdot(\tilde{q} \mathbf u) + (\gamma/2 - 1) \tilde{q} \nabla \cdot u = 0 \,.\end{split}\]

Where \(\tilde{\mathbf B}\) (resp. \(\tilde{q}\)) is either \(\mathbf B\) (resp. \(q\)) for full-f models, \(\mathbf B_0\) (resp. \(q_0\)) for linear models or \(\mathbf B_0+\mathbf B\) (resp. \(q_0+q\)) for \(\delta f\) models.

On the logical domain:

\[\begin{split}\begin{align} &\int_{\hat{\Omega}} \partial_t ( \hat{\rho}^3 \hat{\mathbf{u}}) \cdot G \hat{\mathbf{v}} \, \textrm d \boldsymbol \eta - \int_{\hat{\Omega}} \hat{\mathbf{B}}^2 \cdot G \,\nabla \times (\hat{\mathbf{B}}^2 \times \hat{\mathbf{v}}) \,\frac{1}{\sqrt g}\, - \frac{q}{\gamma -1} \nabla \cdot (\hat{q} \hat{v}) \textrm d \boldsymbol \eta = 0 ~ , \\[2mm] &\partial_t \hat{\mathbf{B}}^2 + \nabla \times (\hat{\mathbf{B}}^2 \times \hat{\mathbf{u}}) = 0 ~ , \\[2mm] &\partial_t \hat{q} + \nabla \cdot (\hat{q} \hat{u}) + (\gamma/2 - 1 ) \hat{p} \nabla \cdot G \hat{u} = 0 \, \end{align}\end{split}\]

It is discretized as

\[\begin{split}\begin{align} &\mathbb M^v[\hat{\rho}_h^{n}] \frac{ \mathbf u^{n+1}-\mathbf u^n}{\Delta t} - (\mathbb C \hat{\Pi}^{1}[\hat{\mathbf B_h^{n}} \cdot \vec{\boldsymbol \Lambda}^v])^\top \mathbb M^2 \mathbf B^{n+\frac{1}{2}}- \Big(\big(\mathbb D \hat{\Pi}^{2}[\hat{q_h^n} \cdot \vec{\boldsymbol \Lambda}^v]] + (\gamma/2 - 1)\hat{\Pi}^{3}[\hat{q_h^n} \cdot \vec{\boldsymbol \Lambda}^3] \mathbb D \mathcal{U}^v \big) \mathbf v \Big)^\top \hat{l}^3(\frac{2q^{n+\frac{1}{2}}}{\gamma-1}) = 0 ~ , \\[2mm] &\frac{\mathbf b^{n+1}- \mathbf b^n}{\Delta t} + \mathbb C \hat{\Pi}^{1}[\hat{\mathbf B_h^{n+\frac{1}{2}}} \cdot \vec{\boldsymbol \Lambda}^v]] \mathbf u^{n+1/2} = 0 ~ , \\[2mm] &\frac{\mathbf q^{n+1}- \mathbf q^n}{\Delta t} + \big(\mathbb D \hat{\Pi}^{2}[\hat{q_h^n} \cdot \vec{\boldsymbol \Lambda}^v]] + (\gamma/2 - 1)\hat{\Pi}^{3}[\hat{q_h^n} \cdot \vec{\boldsymbol \Lambda}^3] \mathbb D \mathcal{U}^v \big) \mathbf u^{n+1/2}= 0 ~ , \\[2mm] \end{align}\end{split}\]

with

\[\hat{l}^3(f)_{ijk}=\int_{\hat{\Omega}} f \Lambda^3_{ijk} \textrm d \boldsymbol \eta\]

where weights in the the BasisProjectionOperator and the WeightedMassOperator are given by

\[\hat{\mathbf{B}}_h^{n+1/2} = (\mathbf{b}^{n+\frac{1}{2}})^\top \vec{\boldsymbol \Lambda}^2 \in V_h^2 \, \qquad \hat{\rho}_h^{n} = (\boldsymbol \rho^{n})^\top \vec{\boldsymbol \Lambda}^3 \in V_h^3 \, \qquad \hat{q}_h^{n+1/2} = (\boldsymbol q^{n+1/2})^\top \vec{\boldsymbol \Lambda}^3 \in V_h^3 \,.\]

and \(\mathcal{U}^v\) is BasisProjectionOperators.

class Variables

Bases: object

Container for variables advanced by VariationalQBEvolve.

Variables:

• q (FEECVariable) – Thermal-pressure-like scalar in "L2" space.

• u (FEECVariable) – Velocity variable in "H1vec" space.

• b (FEECVariable) – Magnetic-field variable in "Hdiv" space.

class Options(model: Literal['full_q', 'linear_q', 'deltaf_q'] = 'full_q', gamma: float = 1.6666666666666667, solver: Literal['pcg', 'cg'] = 'pcg', precond: Literal['MassMatrixPreconditioner', 'MassMatrixDiagonalPreconditioner', None] = 'MassMatrixPreconditioner', solver_params: SolverParameters | None = None, nonlin_solver: NonlinearSolverParameters | None = None, div_u: FEECVariable | None = None, u2: FEECVariable | None = None, qt3: FEECVariable | None = None, bt2: FEECVariable | None = None)

Bases: OptionsBase

Configuration options for VariationalQBEvolve.

Parameters:

• model ({"full_q", "linear_q", "deltaf_q"}, default="full_q") – Coupled q-B evolution model variant.

• gamma (float, default=5/3) – Adiabatic index.

• solver (LiteralOptions.OptsSymmSolver, default="pcg") – Linear solver for implicit substeps.

• precond (LiteralOptions.OptsMassPrecond, default="MassMatrixPreconditioner") – Preconditioner used in linear solves.

• solver_params (SolverParameters, default=None) – Linear-solver controls.

• nonlin_solver (NonlinearSolverParameters, default=None) – Nonlinear/Picard iteration controls.

• div_u (FEECVariable, default=None) – Optional external divergence-of-velocity field.

• u2 (FEECVariable, default=None) – Optional external velocity in 2-form representation.

• qt3 (FEECVariable, default=None) – Optional q background field.

• bt2 (FEECVariable, default=None) – Optional magnetic background field.

class struphy.propagators.variational_resistivity.VariationalResistivity

Bases: Propagator

FEEC discretization of the following equations: find \(s \in L^2\) and \(\mathbf B \in H(\textrm{div})\) such that

\[\begin{split}&\int_\Omega \partial_t \mathbf B \cdot \mathbf v \,\textrm d \mathbf x + \int_\Omega (\eta + \eta_a(\mathbf x)) \nabla \times \mathbf B \cdot \nabla \times \mathbf v \,\textrm d \mathbf x = 0 \qquad \forall \, \mathbf v \in H(\textrm{div}) \,, \\[4mm] &\int_\Omega \frac{\partial \mathcal U}{\partial s} \partial_t s \, q\,\textrm d \mathbf x - \int_\Omega (\eta + \eta_a(\mathbf x)) |\nabla \times \mathbf B|^2 \, q \,\textrm d \mathbf x = 0 \qquad \forall \, q \in L^2\, \text{if using } s, \\[4mm] &\int_\Omega \frac{1}{\gamma - 1} \partial_t p \, q\,\textrm d \mathbf x - \int_\Omega (\eta + \eta_a(\mathbf x)) |\nabla \times \mathbf B|^2 \, q \,\textrm d \mathbf x = 0 \qquad \forall \, q \in L^2\, \text{if using } p.\end{split}\]

With \(\eta_a(\mathbf x) = \eta_a |\nabla \times \mathbf B(\mathbf x)|\)

On the logical domain:

\[\begin{split}\begin{align} &\partial_t \hat{\boldsymbol B} - \eta \Delta \hat{\boldsymbol B} = 0 ~ , \\[2mm] &\int_{\hat{\Omega}} \partial_t (\hat{\rho} \mathcal U(\hat{\rho}, \hat{s})) \hat{w} \,\frac{1}{\sqrt g}\, \textrm d \boldsymbol \eta - \int_{\hat{\Omega}} (\eta + \eta_a(\mathbf \eta)) |DF^{-T}\tilde{\nabla} \times \hat{\boldsymbol B}|^2 \hat{w} \, \textrm d \boldsymbol \eta = 0 ~ , \text{if using } s, \\[2mm] &\int_{\hat{\Omega}} \partial_t (\frac{1}{\gamma -1} \hat{p} ) \hat{w} \,\frac{1}{\sqrt g}\, \textrm d \boldsymbol \eta - \int_{\hat{\Omega}} (\eta + \eta_a(\mathbf \eta)) |DF^{-T}\tilde{\nabla} \times \hat{\boldsymbol B}|^2 \hat{w} \, \textrm d \boldsymbol \eta = 0 ~, \text{if using } p. \end{align}\end{split}\]

It is discretized as

\[\begin{split}\begin{align} &\frac{\mathbf B^{n+1}-\mathbf B^n}{\Delta t} + \, \mathbb C \mathbb M_1^{-1} (\eta M_1 + \eta_a M_1[|\nabla \times \mathbf B|]) M_1^{-1} \mathbb C^T \mathbb M_2 \mathbf B^{n+1} = 0 ~ , \\[2mm] &\frac{P^{3}(\rho e(s^{n+1})- P^{3}(\rho e(s^{n}))}{\Delta t} - P^3((\eta + \eta_a(\mathbf x)) DF^{-T} \tilde{\mathbb C} \frac{ \mathbf B^{n+1}+\mathbf B^n}{2} \cdot DF^{-T} \tilde{\mathbb C} \mathbf B^{n+1}) = 0 ~ , \text{if using } s, \\[2mm] &\frac{1}{\gamma -1}\frac{p^{n+1}- p^{n}}{\Delta t} - P^3((\eta + \eta_a(\mathbf x)) DF^{-T} \tilde{\mathbb C} \frac{ \mathbf B^{n+1}+\mathbf B^n}{2} \cdot DF^{-T} \tilde{\mathbb C} \mathbf B^{n+1}) = 0 ~ , \text{if using } p. \end{align}\end{split}\]

where $P^3$ denotes the $L^2$ projection in the last space of the de Rham sequence.

class Variables

Bases: object

Container for variables advanced by VariationalResistivity.

Variables:

• s (FEECVariable) – Thermodynamic scalar variable in "L2" space.

• b (FEECVariable) – Magnetic-field variable in "Hdiv" space.

class Options(model: Literal['full', 'full_p', 'full_q', 'linear_p', 'linear_q', 'deltaf_q'] = 'full', gamma: float = 1.6666666666666667, solver: Literal['pcg', 'cg'] = 'pcg', precond: Literal['MassMatrixPreconditioner', 'MassMatrixDiagonalPreconditioner', None] = 'MassMatrixDiagonalPreconditioner', solver_params: SolverParameters | None = None, nonlin_solver: NonlinearSolverParameters | None = None, linearize_current: bool = False, rho: FEECVariable | None = None, pt3: FEECVariable | None = None, eta: float = 0.0, eta_a: float = 0.0)

Bases: OptionsBase

Configuration options for VariationalResistivity.

Parameters:

• model ({"full", "full_p", "full_q", "linear_p", "linear_q", "deltaf_q"}, default="full") – Thermodynamic model variant.

• gamma (float, default=5/3) – Adiabatic index.

• solver (LiteralOptions.OptsSymmSolver, default="pcg") – Linear solver for implicit subproblems.

• precond (LiteralOptions.OptsMassPrecond, default="MassMatrixDiagonalPreconditioner") – Preconditioner used in linear solves.

• solver_params (SolverParameters, default=None) – Linear-solver controls.

• nonlin_solver (NonlinearSolverParameters, default=None) – Nonlinear iteration controls.

• linearize_current (bool, default=False) – If True, linearize current terms around background fields.

• rho (FEECVariable, default=None) – Density variable used by variational forms.

• pt3 (FEECVariable, default=None) – Optional pressure-like background field.

• eta (float, default=0.0) – Physical resistivity coefficient.

• eta_a (float, default=0.0) – Artificial-resistivity coefficient.

class struphy.propagators.variational_viscosity.VariationalViscosity

Bases: Propagator

FEEC discretization of the following equations: find \(s \in L^2\) and \(\mathbf u \in (H^1)^3\) such that

\[\begin{split}&\int_\Omega \partial_t (\rho \mathbf u) \cdot \mathbf v\,\textrm d \mathbf x + \int_\Omega (\mu + \mu_a(\mathbf x)) \nabla \mathbf u : \nabla \mathbf v \,\textrm d \mathbf x = 0 \qquad \forall \, \mathbf v \in (H^1)^3 \,, \\[4mm] &\int_\Omega \frac{\partial \mathcal U}{\partial s} \partial_t s \, q \,\textrm d \mathbf x - \mu \int_\Omega |\nabla \mathbf u|^2 \, q \,\textrm d \mathbf x = 0 \qquad \forall \, q \in L^2\,\text{if using } s, \\[4mm] &\int_\Omega \frac{1}{\gamma - 1} \partial_t p \, q\,\textrm d \mathbf x - \mu \int_\Omega |\nabla \mathbf u|^2 \, q \,\textrm d \mathbf x = 0 \qquad \forall \, q \in L^2\, \text{if using } p.\end{split}\]

With \(\mu_a(\mathbf x) = \mu_a |\nabla \mathbf u(\mathbf x)|\)

On the logical domain:

\[\begin{split}\begin{align} &\int_{\hat{\Omega}} \partial_t ( \hat{\rho}^3 \hat{\mathbf{u}}) \cdot G \hat{\mathbf{v}} \, \textrm d \boldsymbol \eta + \mu \int_{\hat{\Omega}} \nabla (DF \hat{\mathbf{u}}) : \nabla (DF \hat{\mathbf{v}}) \,\frac{1}{\sqrt g}\, \textrm d \boldsymbol \eta = 0 ~ , \\[2mm] &\int_{\hat{\Omega}} \partial_t (\hat{\rho} \hat{e}(\hat{\rho}, \hat{s})) \hat{w} \,\frac{1}{\sqrt g}\, \textrm d \boldsymbol \eta - \int_{\hat{\Omega}} (\mu + \mu_a(\boldsymbol \eta)) \nabla (DF \hat{\mathbf{u}}) : \nabla (DF \hat{\mathbf{u}}) \hat{w} \, \textrm d \boldsymbol \eta = 0 ~ , \text{if using } s, \\[2mm] &\int_{\hat{\Omega}} \partial_t (\frac{1}{\gamma -1} \hat{p} ) \hat{w} \,\frac{1}{\sqrt g}\, \textrm d \boldsymbol \eta - \int_{\hat{\Omega}} (\mu + \mu_a(\boldsymbol \eta)) \nabla (DF \hat{\mathbf{u}}) : \nabla (DF \hat{\mathbf{u}}) \hat{w} \, \textrm d \boldsymbol \eta = 0 ~, \text{if using } p. \end{align}\end{split}\]

It is discretized as

\[\begin{split}\begin{align} &\mathbb M^v[\hat{\rho}_h^{n}] \frac{ \mathbf u^{n+1}-\mathbf u^n}{\Delta t} + \sum_\nu (\mathbb G \mathcal{X}^v_\nu)^T (\mu \mathbb M_0 + \mu_a \mathbb M_0[|\nabla u|] \mathbb G \mathcal{X}^v_\nu \mathbf u^{n+1} = 0 ~ , \\[2mm] &\frac{P^{3}(\hat{\rho}_h^{n}\mathcal U(\hat{\rho}_h^{n},\hat{s}_h^{n}))- P^{3}(\hat{\rho}_h^{n}\mathcal U(\hat{\rho}_h^{n},\hat{s}_h^{n+1}))}{\Delta t} - \mu P^3(\sum_\nu DF \mathcal{X}^v_\nu \frac{ \mathbf u^{n+1}+\mathbf u^n}{2} \cdot DF \mathcal{X}^v_\nu \mathbf u^{n+1}) = 0 ~ , \text{if using } s, \\[2mm] &\frac{1}{\gamma -1}\frac{p^{n+1}- p^{n}}{\Delta t} - \mu P^3(\sum_\nu DF \mathcal{X}^v_\nu \frac{ \mathbf u^{n+1}+\mathbf u^n}{2} \cdot DF \mathcal{X}^v_\nu \mathbf u^{n+1}) = 0 ~ , \text{if using } p. \end{align}\end{split}\]

where $P^3$ denotes the $L^2$ projection in the last space of the de Rham sequence and the weights in \(\mathbb M_0[|\nabla u|]\) are given by

\[P^0(g \sqrt{\sum_\nu |(\mathbb G \mathcal{X}^v_\nu \mathbb u)^\top \vec{\boldsymbol \Lambda}^0 |^2]})^\top \vec{\boldsymbol \Lambda}^0 ~.\]

class Variables

Bases: object

Container for variables advanced by VariationalViscosity.

Variables:

• s (FEECVariable) – Thermodynamic scalar variable in "L2" space.

• u (FEECVariable) – Velocity variable in "H1vec" space.

class Options(model: Literal['full', 'full_p', 'full_q', 'linear_p', 'linear_q', 'deltaf_q'] = 'full', gamma: float = 1.6666666666666667, solver: Literal['pcg', 'cg'] = 'pcg', precond: Literal['MassMatrixPreconditioner', 'MassMatrixDiagonalPreconditioner', None] = 'MassMatrixDiagonalPreconditioner', solver_params: SolverParameters | None = None, nonlin_solver: NonlinearSolverParameters | None = None, rho: FEECVariable | None = None, pt3: FEECVariable | None = None, mu: float = 0.0, mu_a: float = 0.0, alpha: float = 0.0)

Bases: OptionsBase

Configuration options for VariationalViscosity.

Parameters:

• model ({"full", "full_p", "full_q", "linear_p", "linear_q", "deltaf_q"}, default="full") – Thermodynamic model variant.

• gamma (float, default=5/3) – Adiabatic index.

• solver (LiteralOptions.OptsSymmSolver, default="pcg") – Linear solver for implicit subproblems.

• precond (LiteralOptions.OptsMassPrecond, default="MassMatrixDiagonalPreconditioner") – Preconditioner used in linear solves.

• solver_params (SolverParameters, default=None) – Linear-solver controls.

• nonlin_solver (NonlinearSolverParameters, default=None) – Nonlinear iteration controls.

• rho (FEECVariable, default=None) – Density variable used by variational forms.

• pt3 (FEECVariable, default=None) – Optional equilibrium/background pressure-like field.

• mu (float, default=0.0) – Physical viscosity coefficient.

• mu_a (float, default=0.0) – Artificial-viscosity coefficient.

• alpha (float, default=0.0) – Optional linear damping/regularization parameter.

Cold plasma

class struphy.propagators.ohm_cold.OhmCold

Bases: Propagator

FEEC discretization of the following equations: find \(\mathbf j \in H(\textnormal{curl})\) and \(\mathbf E \in H(\textnormal{curl})\) such that

\[\begin{split}\int_\Omega \frac{1}{n_0} \frac{\partial \mathbf j}{\partial t} \cdot \mathbf F \,\textrm d \mathbf x &= \frac{1}{\varepsilon} \int_\Omega \mathbf E \cdot \mathbf F \,\textrm d \mathbf x \qquad \forall \,\mathbf F \in H(\textnormal{curl})\,, \\[2mm] -\frac{\partial \mathbf E}{\partial t} &= \frac{\alpha^2}{\varepsilon} \mathbf j \,,\end{split}\]

Time discretization: Crank-Nicolson (implicit mid-point). System size reduction via SchurSolver, such that

\[\begin{split}\begin{bmatrix} \mathbf j^{n+1} - \mathbf j^n \\ \mathbf e^{n+1} - \mathbf e^n \end{bmatrix} = \frac{\Delta t}{2} \begin{bmatrix} 0 & \frac{1}{\varepsilon} \mathbb M_{1/n_0}^{-1} \\ - \frac{1}{\varepsilon} \mathbb M_{1/n_0}^{-1} & 0 \end{bmatrix} \begin{bmatrix} \alpha^2 \mathbb M_{1/n_0} (\mathbf j^{n+1} + \mathbf j^{n}) \\ \mathbb M_1 (\mathbf e^{n+1} + \mathbf e^{n}) \end{bmatrix} \,.\end{split}\]

class Variables

Bases: object

Container for variables advanced by OhmCold.

Variables:

• j (FEECVariable) – Current variable in "Hcurl" space.

• e (FEECVariable) – Electric-field variable in "Hcurl" space.

class Options(solver: Literal['pcg', 'cg'] = 'pcg', precond: Literal['MassMatrixPreconditioner', 'MassMatrixDiagonalPreconditioner', None] = 'MassMatrixPreconditioner', solver_params: SolverParameters | None = None)

Bases: OptionsBase

Configuration options for OhmCold.

Parameters:

• solver (LiteralOptions.OptsSymmSolver, default="pcg") – Symmetric iterative solver used by SchurSolver.

• precond (LiteralOptions.OptsMassPrecond, default="MassMatrixPreconditioner") – Preconditioner applied to the weighted mass matrix M1ninv.

• solver_params (SolverParameters, default=None) – Iterative-solver controls. If None, defaults to SolverParameters().

class struphy.propagators.jxb_cold.JxBCold

Bases: Propagator

FEEC discretization of the following equations: find \(\mathbf j \in H(\textnormal{curl})\) such that

\[\int_\Omega \frac{1}{n_0} \frac{\partial \mathbf j}{\partial t} \cdot \mathbf F \,\textrm d \mathbf x = \frac{1}{\varepsilon} \int_\Omega \frac{1}{n_0} (\mathbf j \times \mathbf B_0) \cdot \mathbf F \,\textrm d \mathbf x \qquad \forall \,\mathbf F \in H(\textnormal{curl})\,,\]

Time discretization: Crank-Nicolson (implicit mid-point), such that

\[\mathbb M_{1/n_0} \left( \mathbf j^{n+1} - \mathbf j^n \right) = \frac{\Delta t}{2} \frac{1}{\varepsilon} \mathbb M_{B_0/n_0} \left( \mathbf j^{n+1} - \mathbf j^n \right)\,.\]

class Variables

Bases: object

Container for variables advanced by JxBCold.

Variables:

j (FEECVariable) – Current variable in "Hcurl" space.

class Options(solver: Literal['pcg', 'cg'] = 'pcg', precond: Literal['MassMatrixPreconditioner', 'MassMatrixDiagonalPreconditioner', None] = 'MassMatrixPreconditioner', solver_params: SolverParameters | None = None)

Bases: OptionsBase

Configuration options for JxBCold.

Parameters:

• solver (LiteralOptions.OptsSymmSolver, default="pcg") – Symmetric iterative solver used for the linear current update.

• precond (LiteralOptions.OptsMassPrecond, default="MassMatrixPreconditioner") – Preconditioner applied to the weighted mass matrix M1ninv.

• solver_params (SolverParameters, default=None) – Iterative-solver controls. If None, defaults to SolverParameters().

Diffusion

class struphy.propagators.implicit_diffusion.ImplicitDiffusion

Bases: Propagator

Weak, implicit discretization of the diffusion (or heat) equation (can be used as a Poisson solver too).

Find \(\phi \in H^1\) such that

\[\int_\Omega \psi\, n_0(\mathbf x)\frac{\partial \phi}{\partial t}\,\textrm d \mathbf x + \int_\Omega \nabla \psi^\top D_0(\mathbf x) \nabla \phi \,\textrm d \mathbf x = \sum_i \int_\Omega \psi\, \rho_i(\mathbf x)\,\textrm d \mathbf x \qquad \forall \ \psi \in H^1\,,\]

where \(n_0, \rho_i:\Omega \to \mathbb R\) are real-valued functions and \(D_0:\Omega \to \mathbb R^{3\times 3}\) is a positive diffusion matrix. Boundary terms from integration by parts are assumed to vanish. The equation is discretized as

\[\left( \frac{\sigma_1}{\Delta t} \mathbb M^0_{n_0} + \mathbb G^\top \mathbb M^1_{D_0} \mathbb G \right)\, \boldsymbol\phi^{n+1} = \frac{\sigma_2}{\Delta t} \mathbb M^0_{n_0} \boldsymbol\phi^{n} + \frac{\sigma_3}{\Delta t} \sum_i(\Lambda^0, \rho_i )_{L^2}\,,\]

where \(M^0_{n_0}\) and \(M^1_{D_0}\) are WeightedMassOperators and \(\sigma_1, \sigma_2, \sigma_3 \in \mathbb R\) are artificial parameters that can be tuned to change the model (see Notes).

Notes

• \(\sigma_1=\sigma_2=0\) and \(\sigma_3 = \Delta t\): Poisson solver with a given charge density \(\sum_i\rho_i\).

• \(\sigma_2=0\) and \(\sigma_1 = \sigma_3 = \Delta t\) : Poisson with adiabatic electrons.

• \(\sigma_1=\sigma_2=1\) and \(\sigma_3 = 0\): Implicit heat equation.

class Variables

Bases: object

Variables:

phi (FEECVariable) – Scalar solution field in H^1 space.

class Options(sigma_1: float = 1.0, sigma_2: float = 0.0, sigma_3: float = 1.0, divide_by_dt: bool = False, stab_mat: Literal['M0', 'M0ad', 'Id'] = 'M0', diffusion_mat: Literal['M1', 'M1perp'] = 'M1', rho: FEECVariable | Callable | tuple[AccumulatorVector, Particles] | list | None = None, rho_coeffs: float | list | None = None, x0: StencilVector | None = None, solver: Literal['pcg', 'cg'] = 'pcg', precond: Literal['MassMatrixPreconditioner', 'MassMatrixDiagonalPreconditioner', None] = 'MassMatrixPreconditioner', solver_params: SolverParameters | None = None)

Bases: OptionsBase

Configuration options for ImplicitDiffusion.

Parameters:

• sigma_1 (float, default=1.0) –

  Coefficient multiplying the stabilization/mass contribution on the left-hand side.

  With divide_by_dt=True it is interpreted as sigma_1 / dt in the assembled linear system.

• sigma_2 (float, default=0.0) –

  Coefficient multiplying the previous solution contribution stab_mat * phi^n on the right-hand side.

  With divide_by_dt=True it is interpreted as sigma_2 / dt.

• sigma_3 (float, default=1.0) –

  Coefficient multiplying the source terms rho on the right-hand side.

  With divide_by_dt=True it is interpreted as sigma_3 / dt.

• divide_by_dt (bool, default=False) – If True, divide sigma_1, sigma_2 and sigma_3 by the time step passed to __call__(dt) before assembling the system.

• stab_mat ({"M0", "M0ad", "Id"}, default="M0") –

  Stabilization operator used in the term weighted by sigma_1 and sigma_2.

  • "M0": standard weighted 0-form mass operator.

  • "M0ad": adiabatic-electron weighted 0-form mass operator.

  • "Id": identity operator.

• diffusion_mat ({"M1", "M1perp"} or WeightedMassOperator, default="M1") – Diffusion metric in the bilinear form grad.T @ diffusion_mat @ grad. You can pass the name of a pre-built operator in mass_ops or a custom WeightedMassOperator compatible with the codomain of grad.

• rho (FEECVariable or Callable or tuple or list, default=None) –

  Source term(s) on the right-hand side. Accepted entries are:

  • None: zero source.

  • FEECVariable in H1.

  • Callable to be projected to H1 via L2Projector.

  • AccumulatorVector.

  • a list containing any mix of the entries above.

  The tuple form is accepted by typing for compatibility with other propagator interfaces that pair particle data with accumulators.

• rho_coeffs (float or list, default=None) – Multiplicative coefficient(s) for rho sources. If a scalar is provided, it is applied to a single source. If a sequence is provided, its length must match the number of collected sources. If None, all coefficients default to 1.0.

• x0 (StencilVector, default=None) – Initial guess for the iterative linear solver.

• solver (LiteralOptions.OptsSymmSolver, default="pcg") – Name of the symmetric iterative solver passed to psydac.linalg.solvers.inverse().

• precond (LiteralOptions.OptsMassPrecond, default="MassMatrixPreconditioner") – Name of the preconditioner configuration. Currently this class sets pc=None internally, so this option is reserved for compatibility and future extensions.

• solver_params (SolverParameters, default=None) – Iterative-solver controls (for example tol, maxiter, verbose, info, recycle). If None, defaults to SolverParameters().

property sources: list[StencilVector | FEECVariable | AccumulatorVector]

Right-hand side of the equation (sources).

property coeffs: list[float]

Same length as self.sources. Coefficients multiplied with sources before solve (default is 1.0).

property x0

feectools.linalg.stencil.StencilVector or struphy.polar.basic.PolarVector. First guess of the iterative solver.

Only particle variables are updated.

class struphy.propagators.push_deterministic_diffusion.PushDeterministicDiffusion

Bases: Propagator

For each marker \(p\), solves

\[\frac{\textnormal d \mathbf x_p(t)}{\textnormal d t} = - D \, \frac{\nabla u}{ u}\mathbf (\mathbf x_p(t))\,,\]

in logical space given by \(\mathbf x = F(\boldsymbol \eta)\):

\[\frac{\textnormal d \boldsymbol \eta_p(t)}{\textnormal d t} = - G\, D \, \frac{\nabla \Pi^0_{L^2}u_h}{\Pi^0_{L^2} u_h}\mathbf (\boldsymbol \eta_p(t))\,, \qquad [\Pi^0_{L^2, ijk} u_h](\boldsymbol \eta_p) = \frac 1N \sum_{p} w_p \boldsymbol \Lambda^0_{ijk}(\boldsymbol \eta_p)\,,\]

where \(D>0\) is a positive diffusion coefficient.

Available algorithms:

• Explicit from ButcherTableau

class Variables

Bases: object

Container for variables advanced by PushDeterministicDiffusion.

Variables:

var (PICVariable) – Particle variable in "Particles3D" space whose marker positions are advanced.

class Options(butcher: ButcherTableau | None = None, bc_type: tuple = ('periodic', 'periodic', 'periodic'), diff_coeff: float = 1.0)

Bases: OptionsBase

Configuration options for PushDeterministicDiffusion.

Parameters:

• butcher (ButcherTableau, default=None) – Butcher tableau used for explicit Runge-Kutta integration. If None, defaults to ButcherTableau().

• bc_type (tuple, default=("periodic", "periodic", "periodic")) – Boundary-condition types per logical coordinate.

• diff_coeff (float, default=1.0) – Positive deterministic diffusion coefficient \(D\).

Only particle variables are updated.

class struphy.propagators.push_random_diffusion.PushRandomDiffusion

Bases: Propagator

For each marker \(p\), solves

\[\textnormal d \mathbf x_p(t) = \sqrt{2 D} \, \textnormal d \mathbf B_{t}\,,\]

where \(D>0\) is a positive diffusion coefficient and \(\textnormal d \mathbf B_{t}\) is a Wiener process,

\[\mathbf B_{t + \Delta t} - \mathbf B_{t} = \sqrt{\Delta t} \,\mathcal N(0;1)\,,\]

with \(\mathcal N(0;1)\) denoting the standard normal distribution with mean zero and variance one.

Available algorithms:

• forward_euler (1st order)

class Variables

Bases: object

Container for variables advanced by PushRandomDiffusion.

Variables:

var (PICVariable) – Particle variable in "Particles3D" space whose marker positions are advanced.

class Options(butcher: ButcherTableau | None = None, bc_type: tuple = ('periodic', 'periodic', 'periodic'), diff_coeff: float = 1.0)

Bases: OptionsBase

Configuration options for PushRandomDiffusion.

Parameters:

• butcher (ButcherTableau, default=None) – Butcher tableau used for explicit integration. If None, defaults to ButcherTableau().

• bc_type (tuple, default=("periodic", "periodic", "periodic")) – Boundary-condition types per logical coordinate.

• diff_coeff (float, default=1.0) – Positive diffusion coefficient \(D\) used in the stochastic increment.

Particle pushing

Only particle variables are updated.

class struphy.propagators.push_eta.PushEta

Bases: Propagator

For each marker \(p\), solves

\[\frac{\textnormal d \mathbf{x}_p(t)}{\textnormal d t} = \mathbf{v}_p\,,\]

in logical space given by \(\mathbf{x} = F(\boldsymbol{\eta})\):

\[\frac{\textnormal d \boldsymbol{\eta}_p(t)}{\textnormal d t} = DF^{-1}(\boldsymbol{\eta}_p(t)) \,\mathbf{v}_p\,,\]

for constant \(\mathbf{v}_p\). Available algorithms:

• Explicit RK from ButcherTableau

class Variables

Bases: object

Container for variables advanced by PushEta.

Variables:

var (PICVariable or SPHVariable) – Particle variable whose marker positions are advanced.

class Options(butcher: ButcherTableau | None = None)

Bases: OptionsBase

Configuration options for PushEta.

Parameters:

butcher (ButcherTableau, default=None) – Butcher tableau used for explicit Runge-Kutta marker pushing. If None, defaults to ButcherTableau().

Only particle variables are updated.

class struphy.propagators.push_eta_pc.PushEtaPC

Bases: Propagator

For each marker \(p\), solves

\[\frac{\textnormal d \mathbf x_p(t)}{\textnormal d t} = \mathbf v_p + \mathbf U (\mathbf x_p(t))\,,\]

for constant \(\mathbf v_p\) and \(\mathbf U\) in logical space given by \(\mathbf x = F(\boldsymbol \eta)\):

\[\frac{\textnormal d \boldsymbol \eta_p(t)}{\textnormal d t} = DF^{-1}(\boldsymbol \eta_p(t)) \,\mathbf v_p + \textnormal{vec}(\hat{\mathbf U}) \,,\]

where

\[\textnormal{vec}( \hat{\mathbf U}^{1}) = G^{-1}\hat{\mathbf U}^{1}\,,\qquad \textnormal{vec}( \hat{\mathbf U}^{2}) = \frac{\hat{\mathbf U}^{2}}{\sqrt g}\,, \qquad \textnormal{vec}( \hat{\mathbf U}) = \hat{\mathbf U}\,.\]

Available algorithms:

• rk4 (4th order, default)

• forward_euler (1st order)

• heun2 (2nd order)

• rk2 (2nd order)

• heun3 (3rd order)

class Variables

Bases: object

Container for variables advanced by PushEtaPC.

Variables:

var (PICVariable or SPHVariable) – Particle variable whose marker positions are advanced.

class Options(butcher: ButcherTableau | None = None, use_perp_model: bool = True, u_tilde: FEECVariable | None = None, u_space: Literal['Hcurl', 'Hdiv', 'H1vec'] = 'Hdiv')

Bases: OptionsBase

Configuration options for PushEtaPC.

Parameters:

• butcher (ButcherTableau, default=None) – Butcher tableau used for explicit Runge-Kutta marker pushing. If None, defaults to ButcherTableau().

• use_perp_model (bool, default=True) – Flag forwarded to the particle kernel to select the perpendicular model formulation.

• u_tilde (FEECVariable, default=None) – Flow field variable used in the advection term.

• u_space (LiteralOptions.OptsVecSpace, default="Hdiv") – FEEC space used to interpret u_tilde in the pusher kernel.

Only particle variables are updated.

class struphy.propagators.push_vin_efield.PushVinEfield

Bases: Propagator

Push the velocities according to

\[\frac{\text{d} \mathbf{v}_p}{\text{d} t} = \frac{1}{\varepsilon}\mathbf{E}(\mathbf{x}_p) \,,\]

where \(\varepsilon \in \mathbb R\) is a constant. In logical coordinates, given by \(\mathbf x = F(\boldsymbol \eta)\):

\[\frac{\text{d} \mathbf{v}_p}{\text{d} t} = \frac{1}{\varepsilon}DF^{-\top}\hat{\mathbf E}^1(\boldsymbol \eta_p) \,,\]

which is solved analytically. \(\mathbf E\) can optionally be defined through a potential, \(\mathbf E = - \nabla \phi\).

class Variables

Bases: object

Container for variables advanced by PushVinEfield.

Variables:

var (PICVariable or SPHVariable) – Particle variable whose velocities are advanced.

class Options(e_field: FEECVariable | tuple[Callable] | None = None, phi: FEECVariable | Callable | None = None)

Bases: OptionsBase

Configuration options for PushVinEfield.

Parameters:

• e_field (FEECVariable or tuple[Callable], default=None) – Electric field used in velocity pushing. Accepted forms are an Hcurl FEEC variable or a tuple of callables to be projected.

• phi (FEECVariable or Callable, default=None) – Optional electrostatic potential from which the electric field is built as -grad(phi). If provided, it overrides e_field.

Only particle variables are updated.

class struphy.propagators.push_vxb.PushVxB

Bases: Propagator

For each marker \(p\), solves

\[\frac{\textnormal d \mathbf{v}_p(t)}{\textnormal d t} = \frac{1}{\varepsilon} \, \mathbf{v}_p(t) \times (\mathbf{B} + \mathbf{B}_{\text{add}}) \,,\]

where \(\varepsilon = 1/(\hat{\Omega}_c \hat t)\) is a constant scaling factor, and for rotation vector \(\mathbf{B}\) and optional, additional fixed rotation vector \(\mathbf{B}_{\text{add}}\), both given as a 2-form:

\[\mathbf{B} = \frac{DF\, \hat{\mathbf{B}}^2}{\sqrt{g}}\,.\]

Available algorithms:

• analytic

• implicit.

class Variables

Bases: object

Container for variables advanced by PushVxB.

Variables:

ions (PICVariable or SPHVariable) – Particle variable whose velocities are advanced.

class Options(algo: Literal['analytic', 'implicit'] = 'analytic', b2_var: FEECVariable | None = None)

Bases: OptionsBase

Configuration options for PushVxB.

Parameters:

• algo ({"analytic", "implicit"}, default="analytic") – Time stepping algorithm used for the velocity-rotation update.

• b2_var (FEECVariable, default=None) – Optional additional magnetic-field 2-form added to the projected equilibrium field before pushing.

Only particle variables are updated.

class struphy.propagators.push_guiding_center_bx_estar.PushGuidingCenterBxEstar

Bases: Propagator

For each marker \(p\), solves

\[\frac{\textnormal d \mathbf{X}_p(t)}{\textnormal d t} = \frac{\mathbf{E}^* \times \mathbf{b}_0}{B_{\parallel}^*} (\mathbf{X}_p(t)) \,,\]

where

\[\mathbf E^* = -\nabla \phi - \varepsilon \mu_p \nabla |\mathbf B|\,,\qquad \mathbf B^* = \mathbf B + \varepsilon v_\parallel \nabla \times \mathbf{b}_0\,,\qquad B^*_\parallel = \mathbf B^* \cdot \mathbf{b}_0\,,\]

Here, \(\mathbf B = \mathbf{B}_0 + \tilde{\mathbf B}\) can be the full magnetic field (equilibrium + perturbation). The electric potential \(\phi\) and/or the magnetic perturbation \(\tilde{\mathbf B}\) can be ignored by passing None. In logical space this is given by \(\mathbf X = F(\boldsymbol \eta)\):

\[\frac{\textnormal d \boldsymbol{\eta}_p(t)}{\textnormal d t} = \frac{\hat{\mathbf E}^{*1} \times \hat{\mathbf b}^1_0}{\sqrt{g}\,\hat{B}_{\parallel}^{*}} (\boldsymbol \eta_p(t)) \,.\]

Available algorithms:

• Explicit from ButcherTableau

• push_gc_bxEstar_discrete_gradient_1st_order()

• push_gc_bxEstar_discrete_gradient_1st_order_newton()

• push_gc_bxEstar_discrete_gradient_2nd_order()

class Variables

Bases: object

Container for variables advanced by PushGuidingCenterBxEstar.

Variables:

ions (PICVariable) – Guiding-center particle variable in "Particles5D" space.

class Options(phi: FEECVariable | None = None, evaluate_e_field: bool = False, b_tilde: FEECVariable | None = None, algo: Literal['discrete_gradient_2nd_order', 'discrete_gradient_1st_order', 'discrete_gradient_1st_order_newton', 'explicit'] = 'discrete_gradient_1st_order', butcher: ButcherTableau | None = None, maxiter: int = 20, tol: float = 1e-07, mpi_sort: Literal['each', 'last', None] = 'each')

Bases: OptionsBase

Configuration options for PushGuidingCenterBxEstar.

Parameters:

• phi (FEECVariable, default=None) – Electrostatic potential variable in "H1" space. If None, defaults to FEECVariable(space="H1").

• evaluate_e_field (bool, default=False) – If True, evaluate and include electric-field contributions in drift-kinetic kernels.

• b_tilde (FEECVariable, default=None) – Optional magnetic perturbation variable added to the equilibrium magnetic field.

• algo ({"discrete_gradient_2nd_order", "discrete_gradient_1st_order", "discrete_gradient_1st_order_newton", "explicit"}, default="discrete_gradient_1st_order") – Guiding-center pushing algorithm.

• butcher (ButcherTableau, default=None) – Butcher tableau used in explicit mode. If None and algo="explicit", defaults to ButcherTableau().

• maxiter (int, default=20) – Maximum number of fixed-point or Newton iterations in discrete-gradient modes.

• tol (float, default=1e-7) – Convergence tolerance for iterative discrete-gradient updates.

• mpi_sort (LiteralOptions.OptsMPIsort, default="each") – MPI sorting policy for particle exchange.

Only particle variables are updated.

class struphy.propagators.push_guiding_center_parallel.PushGuidingCenterParallel

Bases: Propagator

For each marker \(p\), solves

\[\begin{split}\left\{ \begin{aligned} \frac{\textnormal d \mathbf X_p(t)}{\textnormal d t} &= v_{\parallel,p}(t) \frac{\mathbf B^*}{B^*_\parallel}(\mathbf X_p(t)) \,, \\ \frac{\textnormal d v_{\parallel,p}(t)}{\textnormal d t} &= \frac{1}{\varepsilon} \frac{\mathbf B^*}{B^*_\parallel} \cdot \mathbf E^* (\mathbf X_p(t)) \,, \end{aligned} \right.\end{split}\]

where

\[\mathbf E^* = -\nabla \phi - \varepsilon \mu_p \nabla |\mathbf B|\,,\qquad \mathbf B^* = \mathbf B + \varepsilon v_\parallel \nabla \times \mathbf b_0\,,\qquad B^*_\parallel = \mathbf B^* \cdot \mathbf b_0\,,\]

where \(\mathbf B = \mathbf B_0 + \tilde{\mathbf B}\) can be the full magnetic field (equilibrium + perturbation). The electric potential phi and/or the magnetic perturbation b_tilde can be ignored by passing None. In logical space this is given by \(\mathbf X = F(\boldsymbol \eta)\):

\[\begin{split}\left\{ \begin{aligned} \frac{\textnormal d \boldsymbol \eta_p(t)}{\textnormal d t} &= v_{\parallel,p}(t) \frac{\hat{\mathbf B}^{*2}}{\hat B^{*3}_\parallel}(\boldsymbol \eta_p(t)) \,, \\ \frac{\textnormal d v_{\parallel,p}(t)}{\textnormal d t} &= \frac{1}{\varepsilon} \frac{\hat{\mathbf B}^{*2}}{\hat B^{*3}_\parallel} \cdot \hat{\mathbf E}^{*1} (\boldsymbol \eta_p(t)) \,. \end{aligned} \right.\end{split}\]

Available algorithms:

• Explicit from ButcherTableau

• push_gc_Bstar_discrete_gradient_1st_order()

• push_gc_Bstar_discrete_gradient_1st_order_newton()

• push_gc_Bstar_discrete_gradient_2nd_order()

class Variables

Bases: object

Container for variables advanced by PushGuidingCenterParallel.

Variables:

ions (PICVariable) – Guiding-center particle variable in "Particles5D" space.

class Options(phi: FEECVariable | None = None, evaluate_e_field: bool = False, b_tilde: FEECVariable | None = None, algo: Literal['discrete_gradient_2nd_order', 'discrete_gradient_1st_order', 'discrete_gradient_1st_order_newton', 'explicit'] = 'discrete_gradient_1st_order', butcher: ButcherTableau | None = None, maxiter: int = 20, tol: float = 1e-07, mpi_sort: Literal['each', 'last', None] = 'each')

Bases: OptionsBase

Configuration options for PushGuidingCenterParallel.

Parameters:

• phi (FEECVariable, default=None) – Electrostatic potential variable in "H1" space. If None, defaults to FEECVariable(space="H1").

• evaluate_e_field (bool, default=False) – If True, evaluate and include electric-field contributions in drift-kinetic kernels.

• b_tilde (FEECVariable, default=None) – Optional magnetic perturbation variable added to the equilibrium magnetic field.

• algo ({"discrete_gradient_2nd_order", "discrete_gradient_1st_order", "discrete_gradient_1st_order_newton", "explicit"}, default="discrete_gradient_1st_order") – Guiding-center pushing algorithm.

• butcher (ButcherTableau, default=None) – Butcher tableau used in explicit mode. If None and algo="explicit", defaults to ButcherTableau().

• maxiter (int, default=20) – Maximum number of fixed-point or Newton iterations in discrete-gradient modes.

• tol (float, default=1e-7) – Convergence tolerance for iterative discrete-gradient updates.

• mpi_sort (LiteralOptions.OptsMPIsort, default="each") – MPI sorting policy for particle exchange.

Weights updates (delta-f schemes)

Particle and FEEC variables are updated.

class struphy.propagators.efield_weights_coupling.EfieldWeightsCoupling

Bases: Propagator

Solves the following substep

\[\begin{split}\begin{align} & \frac{\partial \mathbf{E}}{\partial t} = - \frac{\alpha^2}{\varepsilon} \int \mathbf{v} f_1 \, \text{d} \mathbf{v} \, , \\[2mm] & \frac{\partial f_1}{\partial t} = \frac{1}{v_{\text{th}}^2 \varepsilon} \, \mathbf{E} \cdot \mathbf{v} f_0 \, , \end{align}\end{split}\]

which after discretization and in curvilinear coordinates reads

\[\begin{split}\frac{\text{d}}{\text{d} t} w_p &= \frac{f_{0,p}}{s_{0, p}} \frac{1}{v_{\text{th}}^2 \varepsilon} \left[ DF^{-T} (\mathbb{\Lambda}^1)^T \mathbf{e} \right] \cdot \mathbf{v}_p \, , \\[2mm] \frac{\text{d}}{\text{d} t} \mathbb{M}_1 \mathbf{e} &= - \frac{\alpha^2}{\varepsilon} \frac 1N \sum_p w_p \mathbb{\Lambda}^1 \cdot \left( DF^{-1} \mathbf{v}_p \right) \,.\end{split}\]

This is solved using the Crank-Nicolson method

\[\begin{split}\begin{bmatrix} \mathbb{M}_1 \left( \mathbf{e}^{n+1} - \mathbf{e}^n \right) \mathbf{W}^{n+1} - \mathbf{W}^n \end{bmatrix} = \frac{\Delta t}{2} \begin{bmatrix} 0 & - \mathbb{E} \\ \mathbb{W} & 0 \end{bmatrix} \begin{bmatrix} \mathbf{e}^{n+1} + \mathbf{e}^n \\ \mathbf{V}^{n+1} + \mathbf{V}^n \end{bmatrix} \, ,\end{split}\]

where

\[\begin{split}\mathbb{E} &= \frac{\alpha^2}{\varepsilon} \frac 1N \mathbb{\Lambda}^1 \cdot \left( DF^{-1} \mathbf{v}_p \right) \, , \\[2mm] \mathbb{W} &= \frac{f_{0,p}}{s_{0,p}} \frac{1}{v_\text{th}^2 \varepsilon} \left( DF^{-1} \mathbf{v}_p \right) \cdot \left(\mathbb{\Lambda}^1\right)^T \, ,\end{split}\]

based on the SchurSolver.

The accumulation matrix and vector assembled in Accumulator are

\[BC = \mathbb{E} \mathbb{W} \, , \qquad By_n = \mathbb{E} \mathbf{W} \,.\]

class Variables

Bases: object

Container for variables advanced by EfieldWeightsCoupling.

Variables:

• e (FEECVariable) – Electric-field variable in "Hcurl" space.

• ions (PICVariable) – Particle variable in "Particles6D" or "DeltaFParticles6D" space.

class Options(alpha: float = 1.0, kappa: float = 1.0, solver: Literal['pcg', 'cg'] = 'pcg', precond: Literal['MassMatrixPreconditioner', 'MassMatrixDiagonalPreconditioner', None] = 'MassMatrixPreconditioner', solver_params: SolverParameters | None = None)

Bases: OptionsBase

Configuration options for EfieldWeightsCoupling.

Parameters:

• alpha (float, default=1.0) – Coupling scaling factor used in the electric-field equation.

• kappa (float, default=1.0) – Coupling scaling factor used in the particle-weight update.

• solver (LiteralOptions.OptsSymmSolver, default="pcg") – Symmetric iterative solver used by SchurSolver.

• precond (LiteralOptions.OptsMassPrecond, default="MassMatrixPreconditioner") – Preconditioner applied to the electric-field mass matrix block.

• solver_params (SolverParameters, default=None) – Iterative-solver controls. If None, defaults to SolverParameters().

Vlasov-Maxwell coupling

Particle and FEEC variables are updated.

class struphy.propagators.vlasov_ampere_coupling.VlasovAmpereCoupling

Bases: Propagator

PIC-FEEC discretization of the following equations: find \(\mathbf{E} \in H(\textnormal{curl})\) and \(f\) such that

\[\begin{split}-\int_{\Omega} \frac{\partial \mathbf{E}}{\partial t} \cdot \mathbf{F}\,\textrm d \mathbf x &= \frac{\alpha^2}{\varepsilon} \int_{\Omega} \int_{\mathbb{R}^3} f \mathbf{v} \cdot \mathbf{F} \, \text{d}^3 \mathbf{v} \,\textrm d \mathbf x \qquad \forall \, \mathbf{F} \in H(\textnormal{curl}) \\[2mm] \frac{\partial f}{\partial t} + \frac{1}{\varepsilon}\, \mathbf{E} \cdot \frac{\partial f}{\partial \mathbf{v}} &= 0 .\end{split}\]

Time discretization: Crank-Nicolson (implicit mid-point).

System size reduction via SchurSolver.

class Variables

Bases: object

Container for variables advanced by VlasovAmpereCoupling.

Variables:

• e (FEECVariable) – Electric-field variable in "Hcurl" space.

• ions (PICVariable) – Particle variable in "Particles6D" space.

class Options(solver: Literal['pcg', 'cg'] = 'pcg', precond: Literal['MassMatrixPreconditioner', 'MassMatrixDiagonalPreconditioner', None] = 'MassMatrixPreconditioner', solver_params: SolverParameters | None = None)

Bases: OptionsBase

Configuration options for VlasovAmpereCoupling.

Parameters:

• solver (LiteralOptions.OptsSymmSolver, default="pcg") – Symmetric iterative solver used by SchurSolver.

• precond (LiteralOptions.OptsMassPrecond, default="MassMatrixPreconditioner") – Preconditioner for electric-field mass matrix block.

• solver_params (SolverParameters, default=None) – Solver controls.

class struphy.propagators.maxwell_weak_ampere.MaxwellWeakAmpere

Bases: Propagator

FEEC discretization of the following equations: find \(\mathbf E \in H(\textnormal{curl})\) and \(\mathbf B \in H(\textnormal{div})\) such that

\[\begin{split}\int_{\Omega} \frac{\partial \mathbf E}{\partial t} \cdot \mathbf F \, \textrm d \mathbf x - \int_{\Omega} \mathbf B \cdot \nabla \times \mathbf F \,\textrm d \mathbf x &= 0\,, \qquad \forall \, \mathbf F \in H(\textnormal{curl}) \\[2mm] \frac{\partial \mathbf B}{\partial t} + \nabla\times\mathbf E &= 0\,.\end{split}\]

Time discretization:

• implicit: Crank-Nicolson (implicit mid-point)

• explicit: explicit RK methods from ButcherTableau

System size reduction via SchurSolver.

class Variables

Bases: object

Container for Maxwell equation variables.

Variables:

• e (FEECVariable) – Electric field in H(curl) space.

• b (FEECVariable) – Magnetic field in H(div) space.

class Options(algo: Literal['implicit', 'explicit'] = 'implicit', solver: Literal['pcg', 'cg'] = 'pcg', precond: Literal['MassMatrixPreconditioner', 'MassMatrixDiagonalPreconditioner', None] = 'MassMatrixPreconditioner', solver_params: SolverParameters | None = None, butcher: ButcherTableau | None = None)

Bases: OptionsBase

Configuration options for MaxwellWeakAmpere.

Parameters:

• algo ({"implicit", "explicit"}, default="implicit") –

  Time stepping scheme to use for the Maxwell equations.

  • "implicit": Crank-Nicolson (implicit mid-point) scheme.

  • "explicit": explicit Runge-Kutta methods from ButcherTableau.

• solver (LiteralOptions.OptsSymmSolver, default="pcg") – Name of the symmetric iterative solver passed to psydac.linalg.solvers.inverse().

• precond (LiteralOptions.OptsMassPrecond, default="MassMatrixPreconditioner") – Name of the preconditioner configuration. Currently used only for implicit time stepping.

• solver_params (SolverParameters, default=None) – Iterative-solver controls (for example tol, maxiter, verbose, info, recycle). If None, defaults to SolverParameters().

• butcher (ButcherTableau, default=None) – Butcher tableau for explicit Runge-Kutta methods. Only used when algo="explicit". If None, defaults to ButcherTableau().

Notes

System size reduction is performed via SchurSolver for implicit time stepping to eliminate the magnetic field and reduce the system to the electric field equation.

MHD-kinetic coupling

Particle and FEEC variables are updated.

class struphy.propagators.current_coupling_5d_curlb.CurrentCoupling5DCurlb

Bases: Propagator

FEEC discretization of the following equations: find \(\mathbf U \in \{H(\textnormal{curl}), H(\textnormal{div}), (H^1)^3\}\) and \(\mathbf B \in H(\textnormal{div})\) such that

\[\begin{split}\left\{ \begin{aligned} \int \rho_0 &\frac{\partial \tilde{\mathbf U}}{\partial t} \cdot \mathbf V\, \textnormal{d} \mathbf x = - \frac{A_\textnormal{h}}{A_b} \iint \frac{f^\text{vol}}{B^*_\parallel} v_\parallel^2 (\nabla \times \mathbf b_0) \times \mathbf B \cdot \mathbf V \, \textnormal{d} \mathbf x \textnormal{d} v_\parallel \textnormal{d} \mu \quad \forall \,\mathbf V \in \{H(\textnormal{curl}), H(\textnormal{div}), (H^1)^3\}\,, \\ &\frac{\partial f}{\partial t} = - \frac{1}{B^*_\parallel} v_\parallel (\nabla \times \mathbf b_0) \cdot (\mathbf B \times \tilde{\mathbf U}) \nabla_{v_\parallel}f \,. \end{aligned} \right.\end{split}\]

Time discretization: Crank-Nicolson (implicit mid-point). System size reduction via SchurSolver:

\[\begin{split}\begin{bmatrix} \mathbf u^{n+1} - \mathbf u^n \\ V_\parallel^{n+1} - V_\parallel^n \end{bmatrix} = \frac{\Delta t}{2} \begin{bmatrix} 0 & - (\mathbb{M}^{2,n})^{-1} \left\{ \mathbb{L}^2 \frac{1}{\bar{\sqrt{g}}} \right\}\cdot_\text{vector} \left\{\bar{b}^{\nabla \times}_0 (\bar{B}^\times_f)^\top \bar{V}_\parallel \frac{1}{\bar{\sqrt{g}}}\right\} \frac{1}{\bar B^{*0}_\parallel}) \\ \frac{1}{\bar B^{*0}_\parallel} \left\{\bar{b}^{\nabla \times}_0 (\bar{B}^\times_f)^\top \bar{V}_\parallel \frac{1}{\bar{\sqrt{g}}}\right\}\, \cdot_\text{vector} \left\{\frac{1}{\bar{\sqrt{g}}}(\mathbb{L}²)^\top\right\} (\mathbb{M}^{2,n})^{-1} & 0 \end{bmatrix} \begin{bmatrix} (\mathbb{M}^{2,n})^{-1} (\mathbf u^{n+1} + \mathbf u^n) \\ \frac{A_\textnormal{h}}{A_b} W (V_\parallel^{n+1} + V_\parallel^n) \end{bmatrix} \,.\end{split}\]

For the detail explanation of the notations, see 2022_DriftKineticCurrentCoupling.

class Variables

Bases: object

Container for variables advanced by CurrentCoupling5DCurlb.

Variables:

• u (FEECVariable) – Fluid-like FEEC variable in one of "Hcurl", "Hdiv", or "H1vec".

• energetic_ions (PICVariable) – Energetic-ion particle variable in "Particles5D" space.

class Options(b_tilde: FEECVariable | None = None, ep_scale: float = 1.0, u_space: Literal['Hcurl', 'Hdiv', 'H1vec'] = 'Hdiv', solver: Literal['pcg', 'cg'] = 'pcg', precond: Literal['MassMatrixPreconditioner', 'MassMatrixDiagonalPreconditioner', None] = 'MassMatrixPreconditioner', solver_params: SolverParameters | None = None, filter_params: FilterParameters | None = None)

Bases: OptionsBase

Configuration options for CurrentCoupling5DCurlb.

Parameters:

• b_tilde (FEECVariable, default=None) – Perturbed magnetic field variable used to build the total magnetic field in the coupling term.

• ep_scale (float, default=1.0) – Scaling factor applied to energetic-particle contributions in the accumulation kernel.

• u_space (LiteralOptions.OptsVecSpace, default="Hdiv") – FEEC space used for the unknown u variable.

• solver (LiteralOptions.OptsSymmSolver, default="pcg") – Symmetric iterative solver used by SchurSolver.

• precond (LiteralOptions.OptsMassPrecond, default="MassMatrixPreconditioner") – Preconditioner used for the FEEC mass matrix block.

• solver_params (SolverParameters, default=None) – Iterative-solver controls. If None, defaults to SolverParameters().

• filter_params (FilterParameters, default=None) – Particle-to-grid filtering parameters used by the accumulator. If None, defaults to FilterParameters().

class struphy.propagators.current_coupling_5d_density.CurrentCoupling5DDensity

Bases: Propagator

FEEC discretization of the following equations: find \(\mathbf U \in \{H(\textnormal{curl}), H(\textnormal{div}), (H^1)^3\}\) and \(\mathbf B \in H(\textnormal{div})\) such that

\[\int \rho_0 \frac{\partial \tilde{\mathbf U}}{\partial t} \cdot \mathbf V \, \textnormal{d} \mathbf{x} = \frac{A_\textnormal{h}}{A_b} \frac{1}{\epsilon} \iiint f^\text{vol} \left(1 - \frac{B_\parallel}{B^*_\parallel}\right) \tilde{\mathbf U} \times \mathbf B_f \cdot \mathbf V \, \textnormal{d} \mathbf{x} \textnormal{d} v_\parallel \textnormal{d} \mu \quad \forall \,\mathbf V \in \{H(\textnormal{curl}), H(\textnormal{div}), (H^1)^3\}\,,\]

FE coefficients update:

\[\mathbf u^{n+1} - \mathbf u^n = -\frac{A_\textnormal{h}}{A_b} \frac{1}{\epsilon} \mathbb{L}²{\mathbb{B}}^\times_f \mathbb{N}(1/g) \mathbb{W} \mathbb{N}\left(1- \frac{\hat B^0_\parallel}{\hat B^{*0} _\parallel}\right) (\mathbb{L}²)^\top \frac{\Delta t}{2} \cdot (\mathbf u^{n+1} + \mathbf u^n) \,.\]

For the detail explanation of the notations, see 2022_DriftKineticCurrentCoupling.

class Variables

Bases: object

Container for variables advanced by CurrentCoupling5DDensity.

Variables:

u (FEECVariable) – Fluid-like FEEC variable in one of "Hcurl", "Hdiv", or "H1vec".

class Options(energetic_ions: PICVariable | None = None, b_tilde: FEECVariable | None = None, ep_scale: float = 1.0, u_space: Literal['Hcurl', 'Hdiv', 'H1vec'] = 'Hdiv', solver: Literal['pcg', 'cg'] = 'pcg', precond: Literal['MassMatrixPreconditioner', 'MassMatrixDiagonalPreconditioner', None] = 'MassMatrixPreconditioner', solver_params: SolverParameters | None = None, filter_params: FilterParameters | None = None)

Bases: OptionsBase

Configuration options for CurrentCoupling5DDensity.

Parameters:

• energetic_ions (PICVariable, default=None) – Energetic-ion particle variable providing marker data in "Particles5D" space.

• b_tilde (FEECVariable, default=None) – Perturbed magnetic field variable added to the equilibrium field.

• ep_scale (float, default=1.0) – Scaling factor applied in the energetic-particle accumulation kernel.

• u_space (LiteralOptions.OptsVecSpace, default="Hdiv") – FEEC space used for the unknown u variable.

• solver (LiteralOptions.OptsSymmSolver, default="pcg") – Symmetric iterative solver used for the linear system.

• precond (LiteralOptions.OptsMassPrecond, default="MassMatrixPreconditioner") – Preconditioner applied to the FEEC mass matrix.

• solver_params (SolverParameters, default=None) – Iterative-solver controls. If None, defaults to SolverParameters().

• filter_params (FilterParameters, default=None) – Particle-to-grid filtering parameters used by the accumulator. If None, defaults to FilterParameters().

Particle and FEEC variables are updated.

class struphy.propagators.current_coupling_5d_gradb.CurrentCoupling5DGradB

Bases: Propagator

FEEC discretization of the following equations: find \(\mathbf U \in \{H(\textnormal{curl}), H(\textnormal{div}), (H^1)^3\}\) and \(\mathbf B \in H(\textnormal{div})\) such that

\[\begin{split}\left\{ \begin{aligned} \int \rho_0 &\frac{\partial \tilde{\mathbf U}}{\partial t} \cdot \mathbf V \, \textnormal{d} \mathbf x = - \frac{A_\textnormal{h}}{A_b} \iint \mu \frac{f^\text{vol}}{B^*_\parallel} (\mathbf b_0 \times \nabla B_\parallel) \times \mathbf B \cdot \mathbf V \,\textnormal{d} \mathbf x \textnormal{d} v_\parallel \textnormal{d} \mu \quad \forall \,\mathbf V \in \{H(\textnormal{curl}), H(\textnormal{div}), (H^1)^3\}\,, \\ &\frac{\partial f}{\partial t} = \frac{1}{B^*_\parallel} \left[ \mathbf b_0 \times (\tilde{\mathbf U} \times \mathbf B) \right] \cdot \nabla f\,. \end{aligned} \right.\end{split}\]

Time discretization: Explicit (‘rk4’, ‘forward_euler’, ‘heun2’, ‘rk2’, ‘heun3’).

\[\begin{split}\begin{bmatrix} \dot{\mathbf u}\\ \dot{\mathbf H} \end{bmatrix} = \begin{bmatrix} 0 & (\mathbb{M}^{2,n})^{-1} \mathbb{L}² \frac{1}{\bar{\sqrt{g}}} \frac{1}{\bar B^{*0}_\parallel}\bar{B}^\times_f \bar{G}^{-1} \bar{b}^\times_0 \bar{G}^{-1} \\ -\bar{G}^{-1} \bar{b}^\times_0 \bar{G}^{-1} \bar{B}^\times_f \frac{1}{\bar B^{*0}_\parallel} \frac{1}{\bar{\sqrt{g}}} (\mathbb{L}²)^\top (\mathbb{M}^{2,n})^{-1} & 0 \end{bmatrix} \begin{bmatrix} \mathbb M^{2,n} \mathbf u \\ \frac{A_\textnormal{h}}{A_b} \bar M \bar W \overline{\nabla B}_\parallel \end{bmatrix} \,.\end{split}\]

For the detail explanation of the notations, see 2022_DriftKineticCurrentCoupling.

class Variables

Bases: object

Container for variables advanced by CurrentCoupling5DGradB.

Variables:

• u (FEECVariable) – Fluid-like FEEC variable in one of "Hcurl", "Hdiv", or "H1vec".

• energetic_ions (PICVariable) – Energetic-ion particle variable in "Particles5D" space.

class Options(b_tilde: FEECVariable | None = None, ep_scale: float = 1.0, algo: Literal['discrete_gradient', 'explicit'] = 'explicit', butcher: ButcherTableau | None = None, u_space: Literal['Hcurl', 'Hdiv', 'H1vec'] = 'Hdiv', solver: Literal['pcg', 'cg'] = 'pcg', precond: Literal['MassMatrixPreconditioner', 'MassMatrixDiagonalPreconditioner', None] = 'MassMatrixPreconditioner', solver_params: SolverParameters | None = None, filter_params: FilterParameters | None = None, dg_solver_params: DiscreteGradientSolverParameters | None = None)

Bases: OptionsBase

Configuration options for CurrentCoupling5DGradB.

Parameters:

• b_tilde (FEECVariable, default=None) – Perturbed magnetic field variable added to the equilibrium field.

• ep_scale (float, default=1.0) – Scaling factor applied in particle accumulation.

• algo ({"discrete_gradient", "explicit"}, default="explicit") –

  Time integration algorithm.

  • "explicit": explicit Runge-Kutta update with butcher.

  • "discrete_gradient": nonlinear discrete-gradient update.

• butcher (ButcherTableau, default=None) – Butcher tableau used in explicit mode. If None and algo="explicit", defaults to ButcherTableau().

• u_space (LiteralOptions.OptsVecSpace, default="Hdiv") – FEEC space used for the unknown u variable.

• solver (LiteralOptions.OptsSymmSolver, default="pcg") – Symmetric iterative solver used for mass-matrix inversions.

• precond (LiteralOptions.OptsMassPrecond, default="MassMatrixPreconditioner") – Preconditioner used with the mass matrix.

• solver_params (SolverParameters, default=None) – Iterative-solver controls. If None, defaults to SolverParameters().

• filter_params (FilterParameters, default=None) – Particle-to-grid filtering parameters used by the accumulator. If None, defaults to FilterParameters().

• dg_solver_params (DiscreteGradientSolverParameters, default=None) – Nonlinear solver controls for discrete-gradient mode. If None and algo="discrete_gradient", defaults to DiscreteGradientSolverParameters().

Particle and FEEC variables are updated.

class struphy.propagators.current_coupling_6d_current.CurrentCoupling6DCurrent

Bases: Propagator

FEEC discretization of the following equations: find \(\tilde{\mathbf{U}} \in \{H(\textnormal{curl}), H(\textnormal{div}), (H^1)^3\}\) and \(f\) such that

\[\begin{split}\int_\Omega \rho_0 &\frac{\partial \tilde{\mathbf{U}}}{\partial t} \cdot \mathbf V \,\textrm d \mathbf x = - \frac{A_\textnormal{h}}{A_\textnormal{b}} \frac{1}{\varepsilon} \int_\Omega n_\textnormal{h}\mathbf{u}_\textnormal{h} \times(\mathbf{B}_0+\tilde{\mathbf{B}}) \cdot \mathbf V \,\textrm d \mathbf x \qquad \forall \, \mathbf V \in \{H(\textnormal{curl}), H(\textnormal{div}), (H^1)^3\}\,, \\[2mm] &\frac{\partial f_\textnormal{h}}{\partial t} + \frac{1}{\varepsilon} \Big[(\mathbf{B}_0+\tilde{\mathbf{B}})\times\tilde{\mathbf{U}} \Big] \cdot \frac{\partial f_\textnormal{h}}{\partial \mathbf{v}} =0\,, \\[2mm] &n_\textnormal{h}\mathbf{u}_\textnormal{h}=\int_{\mathbb{R}^3}f_\textnormal{h}\mathbf{v}\,\textnormal{d}^3 \mathbf v\,.\end{split}\]

Time discretization: Crank-Nicolson (implicit mid-point). System size reduction via SchurSolver.

class Variables

Bases: object

Container for variables advanced by CurrentCoupling6DCurrent.

Variables:

• ions (PICVariable) – Energetic-ion particle variable in "Particles6D" space.

• u (FEECVariable) – Fluid-like FEEC variable in one of "Hcurl", "Hdiv", or "H1vec".

class Options(b_tilde: FEECVariable | None = None, u_space: Literal['Hcurl', 'Hdiv', 'H1vec'] = 'Hdiv', solver: Literal['pcg', 'cg'] = 'pcg', precond: Literal['MassMatrixPreconditioner', 'MassMatrixDiagonalPreconditioner', None] = 'MassMatrixPreconditioner', solver_params: SolverParameters | None = None, filter_params: FilterParameters | None = None, boundary_cut: tuple = (0.0, 0.0, 0.0))

Bases: OptionsBase

Configuration options for CurrentCoupling6DCurrent.

Parameters:

• b_tilde (FEECVariable, default=None) – Perturbed magnetic field variable added to the equilibrium field.

• u_space (LiteralOptions.OptsVecSpace, default="Hdiv") – FEEC space used for the unknown u variable.

• solver (LiteralOptions.OptsSymmSolver, default="pcg") – Symmetric iterative solver used by SchurSolver.

• precond (LiteralOptions.OptsMassPrecond, default="MassMatrixPreconditioner") – Preconditioner applied to the FEEC mass matrix block.

• solver_params (SolverParameters, default=None) – Iterative-solver controls. If None, defaults to SolverParameters().

• filter_params (FilterParameters, default=None) – Particle-to-grid filtering parameters used by the accumulator.

• boundary_cut (tuple, default=(0.0, 0.0, 0.0)) – Boundary clipping parameters forwarded to accumulation/pushing kernels.

class struphy.propagators.current_coupling_6d_density.CurrentCoupling6DDensity

Bases: Propagator

FEEC discretization of the following equations: find \(\tilde{\mathbf{U}} \in \{H(\textnormal{curl}), H(\textnormal{div}), (H^1)^3\}\) such that

\[\begin{split}&\int_\Omega \rho_0 \frac{\partial \tilde{\mathbf{U}}}{\partial t} \cdot \mathbf V \,\textrm d \mathbf x = \frac{A_\textnormal{h}}{A_\textnormal{b}} \frac{1}{\varepsilon} \int_\Omega n_\textnormal{h}\tilde{\mathbf{U}} \times(\mathbf{B}_0+\tilde{\mathbf{B}}) \cdot \mathbf V \,\textrm d \mathbf x \qquad \forall \, \mathbf V \in \{H(\textnormal{curl}), H(\textnormal{div}), (H^1)^3\}\,, \\[2mm] &n_\textnormal{h}=\int_{\mathbb{R}^3}f_\textnormal{h}\,\textnormal{d}^3 \mathbf v\,.\end{split}\]

Time discretization: Crank-Nicolson (implicit mid-point).

class Variables

Bases: object

Container for variables advanced by CurrentCoupling6DDensity.

Variables:

u (FEECVariable) – Fluid-like FEEC variable in one of "Hcurl", "Hdiv", or "H1vec".

class Options(energetic_ions: PICVariable | None = None, b_tilde: FEECVariable | None = None, u_space: Literal['Hcurl', 'Hdiv', 'H1vec'] = 'Hdiv', solver: Literal['pcg', 'cg'] = 'pcg', precond: Literal['MassMatrixPreconditioner', 'MassMatrixDiagonalPreconditioner', None] = 'MassMatrixPreconditioner', solver_params: SolverParameters | None = None, filter_params: FilterParameters | None = None, boundary_cut: tuple = (0.0, 0.0, 0.0))

Bases: OptionsBase

Configuration options for CurrentCoupling6DDensity.

Parameters:

• energetic_ions (PICVariable, default=None) – Energetic-ion particle variable providing marker data in "Particles6D" space.

• b_tilde (FEECVariable, default=None) – Perturbed magnetic field variable added to the equilibrium field.

• u_space (LiteralOptions.OptsVecSpace, default="Hdiv") – FEEC space used for the unknown u variable.

• solver (LiteralOptions.OptsSymmSolver, default="pcg") – Symmetric iterative solver used for the linear system.

• precond (LiteralOptions.OptsMassPrecond, default="MassMatrixPreconditioner") – Preconditioner applied to the FEEC mass matrix.

• solver_params (SolverParameters, default=None) – Iterative-solver controls. If None, defaults to SolverParameters().

• filter_params (FilterParameters, default=None) – Particle-to-grid filtering parameters used by the accumulator.

• boundary_cut (tuple, default=(0.0, 0.0, 0.0)) – Boundary clipping parameters forwarded to the accumulation kernel.

class struphy.propagators.shear_alfven_current_coupling_5d.ShearAlfvenCurrentCoupling5D

Bases: Propagator

FEEC discretization of the following equations: find \(\mathbf U \in \{H(\textnormal{curl}), H(\textnormal{div}), (H^1)^3\}\) and \(\mathbf B \in H(\textnormal{div})\) such that

\[\begin{split}\left\{ \begin{aligned} \int \rho_0 &\frac{\partial \tilde{\mathbf U}}{\partial t} \cdot \mathbf V \, \textnormal{d} \mathbf{x} = \int \left(\tilde{\mathbf B} - \frac{A_\textnormal{h}}{A_b} \iint f^\text{vol} \mu \mathbf{b}_0\textnormal{d} \mu \textnormal{d} v_\parallel \right) \cdot \nabla \times (\tilde{\mathbf B} \times \mathbf V) \, \textnormal{d} \mathbf{x} \quad \forall \, \mathbf V \in \{H(\textnormal{curl}), H(\textnormal{div}), (H^1)^3\}\,, \,, \\ &\frac{\partial \tilde{\mathbf B}}{\partial t} = - \nabla \times (\tilde{\mathbf B} \times \tilde{\mathbf U}) \,. \end{aligned} \right.\end{split}\]

Time discretization: Crank-Nicolson (implicit mid-point). System size reduction via SchurSolver:

\[\begin{split}\begin{bmatrix} \mathbf u^{n+1} - \mathbf u^n \\ \mathbf b^{n+1} - \mathbf b^n \end{bmatrix} = \frac{\Delta t}{2} \,. \begin{bmatrix} 0 & (\mathbb M^{2,n})^{-1} \mathcal {T^2}^\top \mathbb C^\top \\ - \mathbb C \mathcal {T^2} (\mathbb M^{2,n})^{-1} & 0 \end{bmatrix} \begin{bmatrix} {\mathbb M^{2,n}}(\mathbf u^{n+1} + \mathbf u^n) \\ \mathbb M_2(\mathbf b^{n+1} + \mathbf b^n) + \sum_k^{N_p} \omega_k \mu_k \hat{\mathbf b}¹_0 (\boldsymbol \eta_k) \cdot \left(\frac{1}{\sqrt{g(\boldsymbol \eta_k)}} \vec \Lambda² (\boldsymbol \eta_k) \right) \end{bmatrix} \,,\end{split}\]

where \(\mathcal{T}^2 = \hat \Pi \left[\frac{\tilde{\mathbf B}^2}{\sqrt{g} \times \vec \Lambda^2\right]\) and \(\mathbb M^{2,n}\) is a WeightedMassOperators being weighted with \(\rho_\text{eq}\), the MHD equilibirum density.

class Variables

Bases: object

Container for variables advanced by ShearAlfvenCurrentCoupling5D.

Variables:

• u (FEECVariable) – Velocity variable in one of "Hcurl", "Hdiv", or "H1vec".

• b (FEECVariable) – Magnetic-field variable in "Hdiv" space.

class Options(energetic_ions: PICVariable | None = None, ep_scale: float = 1.0, u_space: Literal['Hcurl', 'Hdiv', 'H1vec'] = 'Hdiv', algo: Literal['implicit', 'explicit'] = 'implicit', solver: Literal['pcg', 'cg'] = 'pcg', precond: Literal['MassMatrixPreconditioner', 'MassMatrixDiagonalPreconditioner', None] = 'MassMatrixDiagonalPreconditioner', solver_params: SolverParameters | None = None, filter_params: FilterParameters | None = None, butcher: ButcherTableau | None = None, nonlinear: bool = True)

Bases: OptionsBase

Configuration options for ShearAlfvenCurrentCoupling5D.

Parameters:

• energetic_ions (PICVariable, default=None) – Energetic-ion particle variable in "Particles5D" space.

• ep_scale (float, default=1.0) – Scaling factor for particle accumulation terms.

• u_space (LiteralOptions.OptsVecSpace, default="Hdiv") – FEEC space used for velocity variable u.

• algo ({"implicit", "explicit"}, default="implicit") – Time stepping algorithm.

• solver (LiteralOptions.OptsSymmSolver, default="pcg") – Symmetric iterative solver.

• precond (LiteralOptions.OptsMassPrecond, default="MassMatrixDiagonalPreconditioner") – Preconditioner for mass-matrix block.

• solver_params (SolverParameters, default=None) – Solver controls.

• filter_params (FilterParameters, default=None) – Particle-to-grid filtering parameters.

• butcher (ButcherTableau, default=None) – Explicit Runge-Kutta tableau for algo="explicit".

• nonlinear (bool, default=True) – If True, include nonlinear TB operator updates.

Particle and FEEC variables are updated.

class struphy.propagators.pressure_coupling_6d.PressureCoupling6D

Bases: Propagator

FEEC discretization of the following equations: find \(\tilde{\mathbf{U}} \in \{H(\textnormal{curl}), H(\textnormal{div}), (H^1)^3\}\) and \(f\) such that

\[\begin{split}\int_\Omega \rho_0 &\frac{\partial \tilde{\mathbf{U}}}{\partial t} \cdot \mathbf V \,\textrm d \mathbf x = - \frac{A_\textnormal{h}}{A_\textnormal{b}} \nabla \cdot \tilde{\mathbb{P}}_{\textnormal{h},\perp} \cdot \mathbf V \,\textrm d \mathbf x \qquad \forall \, \mathbf V \in \{H(\textnormal{curl}), H(\textnormal{div}), (H^1)^3\}\,, \\[2mm] &\frac{\partial f_\textnormal{h}}{\partial t} - \left(\nabla \tilde{\mathbf U}_\perp \cdot \mathbf{v} \right) \cdot \frac{\partial f_\textnormal{h}}{\partial \mathbf{v}} =0\,, \\[2mm] &\tilde{\mathbb{P}}_{\textnormal{h},\perp} = \int_{\mathbb{R}^3}f_\textnormal{h}\mathbf{v}_\perp \mathbf{v}_\perp^\top \,\textnormal{d}^3 \mathbf v\,.\end{split}\]

Time discretization: Crank-Nicolson (implicit mid-point). System size reduction via SchurSolver:

\[\begin{split}\begin{bmatrix} u^{n+1} - u^n \\ V^{n+1} - V^n \end{bmatrix} = \frac{\Delta t}{2} \begin{bmatrix} 0 & (\mathbb M^n)^{-1} V^\top (\bar {\mathcal X})^\top \mathbb G^\top (\bar {\mathbf \Lambda}^1)^\top \bar {DF}^{-1} \\ - {DF}^{-\top} \bar {\mathbf \Lambda}^1 \mathbb G \bar {\mathcal X} V (\mathbb M^n)^{-1} & 0 \end{bmatrix} \begin{bmatrix} {\mathbb M^n}(u^{n+1} + u^n) \\ \bar W (V^{n+1} + V^{n} \end{bmatrix} \,.\end{split}\]

class Variables

Bases: object

Container for variables advanced by PressureCoupling6D.

Variables:

• u (FEECVariable) – Fluid-like FEEC variable in one of "Hcurl", "Hdiv", or "H1vec".

• energetic_ions (PICVariable) – Energetic-ion particle variable in "Particles6D" space.

class Options(ep_scale: float = 1.0, u_space: Literal['Hcurl', 'Hdiv', 'H1vec'] = 'Hdiv', solver: Literal['pcg', 'cg'] = 'pcg', precond: Literal['MassMatrixPreconditioner', 'MassMatrixDiagonalPreconditioner', None] = 'MassMatrixPreconditioner', solver_params: SolverParameters | None = None, filter_params: FilterParameters | None = None, use_perp_model: bool = True)

Bases: OptionsBase

Configuration options for PressureCoupling6D.

Parameters:

• ep_scale (float, default=1.0) – Scaling factor applied to energetic-particle pressure coupling terms.

• u_space (LiteralOptions.OptsVecSpace, default="Hdiv") – FEEC space used for the unknown u variable.

• solver (LiteralOptions.OptsSymmSolver, default="pcg") – Symmetric iterative solver used by SchurSolver.

• precond (LiteralOptions.OptsMassPrecond, default="MassMatrixPreconditioner") – Preconditioner applied to the mass matrix block for u.

• solver_params (SolverParameters, default=None) – Iterative-solver controls. If None, defaults to SolverParameters().

• filter_params (FilterParameters, default=None) – Particle-to-grid filtering parameters used by the accumulator. If None, defaults to FilterParameters().

• use_perp_model (bool, default=True) – If True, use the perpendicular pressure-coupling model; otherwise use the full-model kernels.

class GT_MAT_G(derham, MAT, transposed=False)

Bases: LinOpWithTransp

Class for defining LinearOperator corresponding to \(G^\top (\text{MAT}) G \in \mathbb{R}^{3N^0 \times 3N^0}\) where \(\text{MAT} = V^\top (\bar {\mathbf \Lambda}^1)^\top \bar{DF}^{-1} \bar{W} \bar{DF}^{-\top} \bar{\mathbf \Lambda}^1 V \in \mathbb{R}^{3N^1 \times 3N^1}\).

Parameters:

• derham (struphy.feec.psydac_derham.Derham) – Discrete de Rham sequence on the logical unit cube.

• MAT (List of StencilMatrices) – List with six of accumulated pressure terms

property domain

The domain of the linear operator - an element of Vectorspace

property codomain

The codomain of the linear operator - an element of Vectorspace

property dtype

The data type of the coefficients of the linear operator, upon convertion to matrix.

property tosparse

Convert to a sparse matrix in any of the formats supported by scipy.sparse.

property toarray

Convert to Numpy 2D array.

transpose()

Transpose the LinearOperator .

If conjugate is True, return the Hermitian transpose.

dot(v, out=None)

dot product between GT_MAT_G and v.

Parameters:

v (StencilVector or BlockVector) – Input FE coefficients from V.coeff_space.

Return type:

A StencilVector or BlockVector from W.coeff_space.

SPH

Only particle variables are updated.

class struphy.propagators.push_vin_sph_pressure.PushVinSPHpressure

Bases: Propagator

For each marker \(p\), solves

\[\frac{\textnormal d \mathbf v_p(t)}{\textnormal d t} = \mathbf g - \sum_{i=1}^N w_i \left( \frac{\partial \mathcal U}{\partial \rho}(\boldsymbol \eta_p) + \frac{\partial \mathcal U}{\partial \rho}(\boldsymbol \eta_i) \right) DF^{-\top}\nabla W_h(\boldsymbol \eta_p - \boldsymbol \eta_i) \,,\]

where \(\mathbf g\) is a constant acceleration and the second term corresponds to the pressure gradient. Here, \(\mathcal U(\rho)\) denotes the internal energy per unit mass as a function of the mass density \(\rho\) and \(DF^{-\top}\) denotes the inverse transpose Jacobian arising in the pull back of the gradient of the smoothing kernel \(W_h\) chosen from sph_smoothing_kernels. Two choices of the internal energy are implemented:

• Isothermal closure: \(\mathcal U(\rho) = \kappa \, \ln(\rho)\), where \(\kappa\) is constant.

• Polytropic closure: \(\mathcal U(\rho) = \kappa \, \rho^{\gamma - 1} / (\gamma - 1)\), where \(\kappa\) is the polytropic constant and \(\gamma = C_p / C_v\) is the polytropic index.

Time stepping is forward Euler.

class Variables

Bases: object

Container for variables advanced by PushVinSPHpressure.

Variables:

fluid (SPHVariable) – SPH particle variable in "ParticlesSPH" space.

class Options(kernel_type: Literal['trigonometric_1d', 'gaussian_1d', 'linear_1d', 'trigonometric_2d', 'gaussian_2d', 'linear_2d', 'trigonometric_3d', 'gaussian_3d', 'linear_isotropic_3d', 'linear_3d'] = 'gaussian_2d', kernel_width: tuple | None = None, algo: Literal['forward_euler'] = 'forward_euler', gravity: tuple = (0.0, 0.0, 0.0), kappa: float = 1.0, thermodynamics: Literal['isothermal', 'polytropic'] = 'isothermal')

Bases: OptionsBase

Configuration options for PushVinSPHpressure.

Parameters:

• kernel_type (LiteralOptions.OptsKernel, default="gaussian_2d") – Smoothing kernel used for density and pressure-force evaluation.

• kernel_width (tuple, default=None) – Kernel widths per logical direction. If None, defaults to (1 / n_i) based on sorting boxes.

• algo ({"forward_euler"}, default="forward_euler") – Time stepping algorithm for velocity pushing.

• gravity (tuple, default=(0.0, 0.0, 0.0)) – Constant gravity vector added in the SPH pressure push.

• kappa (float, default=1.0) – Coefficient in the internal energy function.

• thermodynamics ({"isothermal", "polytropic"}, default="isothermal") – Thermodynamic closure selecting the SPH pressure kernel.

Only particle variables are updated.

class struphy.propagators.push_vin_viscous_potential.PushVinViscousPotential

Bases: Propagator

For each marker \(p\), solves

\[\frac{\textnormal d \mathbf v_p(t)}{\textnormal d t} = \kappa_p \sum_{i=1}^N w_i \left( \frac{1}{\rho^{N,h}(\boldsymbol \eta_p)} + \frac{1}{\rho^{N,h}(\boldsymbol \eta_i)} \right) DF^{-\top}\nabla W_h(\boldsymbol \eta_p - \boldsymbol \eta_i) \,,\]

where \(DF^{-\top}\) denotes the inverse transpose Jacobian, and with the smoothed density

\[\rho^{N,h}(\boldsymbol \eta) = \frac 1N \sum_{j=1}^N w_j \, W_h(\boldsymbol \eta - \boldsymbol \eta_j)\,,\]

where \(W_h(\boldsymbol \eta)\) is a smoothing kernel from sph_smoothing_kernels. Time stepping:

• Explicit from ButcherTableau

class Variables

Bases: object

Container for variables advanced by PushVinViscousPotential.

Variables:

fluid (SPHVariable) – SPH particle variable in "ParticlesSPH" space.

class Options(kernel_type: Literal['trigonometric_1d', 'gaussian_1d', 'linear_1d', 'trigonometric_2d', 'gaussian_2d', 'linear_2d', 'trigonometric_3d', 'gaussian_3d', 'linear_isotropic_3d', 'linear_3d'] = 'gaussian_2d', kernel_width: tuple | None = None, algo: Literal['forward_euler'] = 'forward_euler', mu: float = 1.0)

Bases: OptionsBase

Configuration options for PushVinViscousPotential.

Parameters:

• kernel_type (LiteralOptions.OptsKernel, default="gaussian_2d") – Smoothing kernel used for SPH evaluations.

• kernel_width (tuple, default=None) – Kernel widths per logical direction. If None, defaults to (1 / n_i) based on sorting boxes.

• algo ({"forward_euler"}, default="forward_euler") – Time stepping algorithm for the viscous potential push.

• mu (float, default=1.0) – Dynamic viscosity coefficient used by the viscosity tensor kernel. Must be non-negative.

Miscellaneous

class struphy.propagators.hasegawa_wakatani_step.HasegawaWakataniStep

Bases: Propagator

FEEC discretization of the following equations: find \((n, \omega) \in H^1 \times H^1\) such that

\[\begin{split}&\int_\Omega\frac{\partial n}{\partial t} m \,\textrm d \mathbf x = \int_\Omega C(x, y)(\phi - n) \, m \,\textrm d \mathbf x - \int_\Omega \phi [n, m] \,\textrm d \mathbf x - \kappa \int_\Omega \partial_y \phi \,m \,\textrm d \mathbf x - \nu \int_\Omega \nabla n \cdot \nabla m \,\textrm d \mathbf x \qquad \forall m \in H^1\,, \\[2mm] &\int_\Omega\frac{\partial \omega}{\partial t} \psi \,\textrm d \mathbf x = \int_\Omega C(x, y)(\phi - n) \, \psi \,\textrm d \mathbf x - \int_\Omega \phi [\omega, \psi] \,\textrm d \mathbf x - \nu \int_\Omega \nabla \omega \cdot \nabla \psi \,\textrm d \mathbf x \qquad \forall \psi \in H^1\,,\end{split}\]

where \(\phi \in H^1\) is a given stream function, \(C = C(x, y)\), \(\kappa\) and \(\nu\) are constants and \([a, b] = \partial_x a \partial_y b - \partial_y a \partial_x b\).

Time discretization: explicit Runge-Kutta, see ODEsolverFEEC.

Parameters:

• n0 (StencilVector) – The density.

• omega0 (StencilVector) – The stream function.

• phi (SplineFuncion) – The potential.

• c_fun (str) – Defines the function c(x,y) in front of (phi - n).

• kappa (float) – Equation parameters.

• nu (float) – Equation parameters.

• algo (str) – See ButcherTableau for available algorithms.

• M0_solver (dict) – Solver parameters for M0 inversion.

class Variables

Bases: object

Container for variables advanced by HasegawaWakataniStep.

Variables:

• n (FEECVariable) – Density variable in "H1" space.

• omega (FEECVariable) – Vorticity variable in "H1" space.

class Options(phi: FEECVariable | None = None, c_fun: Literal['const'] = 'const', kappa: float = 1.0, nu: float = 0.01, butcher: ButcherTableau | None = None, solver: Literal['pcg', 'cg'] = 'pcg', precond: Literal['MassMatrixPreconditioner', 'MassMatrixDiagonalPreconditioner', None] = 'MassMatrixPreconditioner', solver_params: SolverParameters | None = None)

Bases: OptionsBase

Configuration options for HasegawaWakataniStep.

Parameters:

• phi (FEECVariable, default=None) – Stream-function variable in "H1" space. If None, a default FEECVariable(space="H1") is allocated.

• c_fun ({"const"}, default="const") – Choice of coupling profile \(C(x,y)\) used in the model.

• kappa (float, default=1.0) – Constant multiplying the background-gradient drift term.

• nu (float, default=0.01) – Diffusion coefficient in density and vorticity equations.

• butcher (ButcherTableau, default=None) – Butcher tableau for explicit Runge-Kutta integration. If None, defaults to ButcherTableau().

• solver (LiteralOptions.OptsSymmSolver, default="pcg") – Symmetric iterative solver used for M0 inversions.

• precond (LiteralOptions.OptsMassPrecond, default="MassMatrixPreconditioner") – Preconditioner used with M0 inversions.

• solver_params (SolverParameters, default=None) – Iterative-solver controls. If None, defaults to SolverParameters().

class struphy.propagators.time_dependent_source.TimeDependentSource

Bases: Propagator

Propagates a source term \(S(t) \in V_h^n\) of the form

\[S(t) = \sum_{ijk} c_{ijk} \Lambda^n_{ijk} * h(\omega t)\,,\]

where \(h(\omega t)\) is one of the functions in Notes.

Notes

• \(h(\omega t) = \cos(\omega t)\) (default)

• \(h(\omega t) = \sin(\omega t)\)

class Variables

Bases: object

Container for variables advanced by TimeDependentSource.

Variables:

source (FEECVariable) – Source coefficient field in "H1" space.

class Options(omega: float = 6.283185307179586, hfun: Literal['cos', 'sin'] = 'cos')

Bases: OptionsBase

Configuration options for TimeDependentSource.

Parameters:

• omega (float, default=2*pi) – Angular frequency of the time-dependent modulation.

• hfun ({"cos", "sin"}, default="cos") – Time modulation function applied to initial coefficients.

class struphy.propagators.two_fluid_quasi_neutral_full.TwoFluidQuasiNeutralFull

Bases: Propagator

FEEC discretization of the following equations: find \(\mathbf u \in H(\textnormal{div})\), \(\mathbf u_e \in H(\textnormal{div})\) and \(\mathbf \phi \in L^2\) such that

\[\begin{split}\int_{\Omega} \partial_t \mathbf{u}\cdot \mathbf{v} \, \textrm d\mathbf{x} &= \int_{\Omega} \phi \nabla \! \cdot \! \mathbf{v} \, \textrm d\mathbf{x} + \int_{\Omega} \mathbf{u}\! \times \! \mathbf{B}_0 \cdot \mathbf{v} \, \textrm d\mathbf{x} + \nu \int_{\Omega} \nabla \mathbf{u}\! : \! \nabla \mathbf{v} \, \textrm d\mathbf{x} + \int_{\Omega} f \mathbf{v} \, \textrm d\mathbf{x} \qquad \forall \, \mathbf{v} \in H(\textrm{div}) \,. \\[2mm] 0 &= - \int_{\Omega} \phi \nabla \! \cdot \! \mathbf{v_e} \, \textrm d\mathbf{x} - \int_{\Omega} \mathbf{u_e} \! \times \! \mathbf{B}_0 \cdot \mathbf{v_e} \, \textrm d\mathbf{x} + \nu_e \int_{\Omega} \nabla \mathbf{u_e} \!: \! \nabla \mathbf{v_e} \, \textrm d\mathbf{x} + \int_{\Omega} f_e \mathbf{v_e} \, \textrm d\mathbf{x} \qquad \forall \ \mathbf{v_e} \in H(\textrm{div}) \,. \\[2mm] 0 &= \int_{\Omega} \psi \nabla \cdot (\mathbf{u}-\mathbf{u_e}) \, \textrm d\mathbf{x} \qquad \forall \, \psi \in L^2 \,.\end{split}\]

Time discretization: fully implicit.

class Variables

Bases: object

Container for variables advanced by TwoFluidQuasiNeutralFull.

Variables:

• u (FEECVariable or None) – Ion velocity variable in "Hdiv" space.

• ue (FEECVariable or None) – Electron velocity variable in "Hdiv" space.

• phi (FEECVariable or None) – Electrostatic potential variable in "L2" space.

class Options(nu: float = 1.0, nu_e: float = 1.0, eps_norm: float = 0.001, boundary_data_u: dict[tuple[int, int], Callable] | None = None, boundary_data_ue: dict[tuple[int, int], Callable] | None = None, source_u: Callable | None = None, source_ue: Callable | None = None, stab_sigma: float | None = None, solver: Literal['pbicgstab', 'bicgstab', 'gmres'] = 'gmres', solver_params: SolverParameters | None = None)

Bases: OptionsBase

Configuration options for TwoFluidQuasiNeutralFull.

Parameters:

• nu (float, default=1.0) – Ion viscosity coefficient.

• nu_e (float, default=1.0) – Electron viscosity coefficient.

• eps_norm (float, default=1e-3) – Normalization/scaling parameter in Lorentz coupling terms.

• boundary_data_u (dict[tuple[int, int], Callable] or None, default=None) – Inhomogeneous Dirichlet data for ion velocity faces.

• boundary_data_ue (dict[tuple[int, int], Callable] or None, default=None) – Inhomogeneous Dirichlet data for electron velocity faces.

• source_u (Callable or None, default=None) – Source term for ion momentum equation.

• source_ue (Callable or None, default=None) – Source term for electron momentum equation.

• stab_sigma (float or None, default=None) – Optional stabilization coefficient for electron block.

• solver (LiteralOptions.OptsGenSolver, default="gmres") – Linear/saddle-point solver used for the global system.

• solver_params (SolverParameters or None, default=None) – Solver controls.