Available Propagators

Field solvers

class struphy.propagators.poisson_solve.PoissonSolve(rho: FEECVariable | ParticlesToGrid | Callable | list | None = None, rho_coeffs: float | list | None = None)

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'] = 'M0', stab_average: str | None = None, diffusion_mat: Literal['M1', 'M1perp', 'M1para', 'M1gyro'] = 'M1', x0: StencilVector | None = None, solver: Literal['pcg', 'cg'] = 'pcg', precond: Literal['MultiGrid', 'StiffnessPreconditioner', 'MassMatrixPreconditioner', 'PolarMassMatrixPreconditioner', None] | None = None, precond_params: dict | None = None, multigrid: MultiGridOptions | None = None, solver_params: SolverParameters | None = None, enforce_compatibility: bool = True, filter_params: dict[PICVariable | SPHVariable, FilterParameters] | 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="M0") –

  Stabilization matrix multiplied by stab_eps.

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

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

  • "Id": identity operator.

• stab_average (str, default=None) – If not None, the stabilization becomes stab_mat @ (I - A), where A is the stab_mat-weighted average over the given logical directions, e.g. "eta3" or "eta2 eta3", see ImplicitDiffusion.

• 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.

• 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.OptsDiffusionPrecond, default=None) – Name of the preconditioner, see ImplicitDiffusion ("MultiGrid" for geometric multigrid, "StiffnessPreconditioner" for the Kronecker approximation of the Laplacian, "MassMatrixPreconditioner" for the Kronecker approximation of M0). Requires solver="pcg".

• precond_params (dict, default=None) – Keyword arguments passed to the constructor of the mass-matrix or stiffness preconditioner, see ImplicitDiffusion.

• multigrid (MultiGridOptions, default=None) – Options of the multigrid preconditioner (if precond="MultiGrid"). For periodic or Neumann boundary conditions with stab_eps = 0, use MultiGridOptions(nullspace="constants").

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

• enforce_compatibility (bool, default=True) –

  If True and the Poisson operator is singular (no Dirichlet boundary conditions and no polar splines, e.g. fully periodic), the constant (“net charge”) component of the right-hand side is removed before the solve,

  \[\mathbf b \leftarrow \mathbf b - \frac{\mathbf 1^\top \mathbf b}{\mathbf 1^\top \mathbb M^0 \mathbf 1}\, \mathbb M^0 \mathbf 1\,,\]

  such that the discrete compatibility condition \(\mathbf 1^\top \mathbf b = 0\) holds (\(\mathbf 1\) is the coefficient vector of the constant function). This is equivalent to a uniform neutralizing background; for particle sources, the net charge is Monte-Carlo noise. Without it, the constant mode of phi becomes (net charge) / stab_eps, and grad phi loses digits to cancellation. stab_eps then only fixes the constant of phi. Set to False if stab_eps * stab_mat is a physical term (e.g. adiabatic electrons with stab_mat="M0ad") and not just a stabilization.

• filter_params (dict[PICVariable | SPHVariable, FilterParameters], default=None) – If not None, specifies a filter to the accumulation of a specific variable.

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.gyrokinetic_poisson_solve.GyrokineticPoissonSolve(rho: FEECVariable | ParticlesToGrid | Callable | list | None = None, rho_coeffs: float | list | None = None, epsilon: float = 1.0, Z: int = 1, diagnostic: FEECVariable | None = None)

Bases: ImplicitDiffusion

Weak discretization of the Poisson equation, with adabatic response of electrons.

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

\[\frac{1}{Z\epsilon^2} \int_\Omega \frac{n_0}{T_0} \psi\, \phi\,\textrm d \mathbf x + \int_\Omega \frac{n_0}{|B_0|²} \nabla \psi^\top \,\mathbb D \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 \(\mathbb D\) can be an anisotropic (gyrokinetic) diffusion operator (or the identity), \(\epsilon \in \mathbb R\) is the gyrokinetic ratio defined in units, and Z the charge number of ions. Boundary terms from integration by parts are assumed to vanish.

The equation is discretized as

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

where \(\mathbb M^1_\textnormal{gyro}\) is the gyrokinetic mass matrix and \(\mathbb S\) is a stabilization matrix.

Parameters:

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

• stab_eps (float) – Stabilization parameter multiplied on stab_mat (default=0.0).

• stab_mat (str) – Name of the stabilizing matrix.

• rho (FEECVariable or Callable or ParticlesToGrid or list) – (List of) right-hand side source term(s) of the Poisson problem (optional, can be set with a setter later). See ImplicitDiffusion for the accepted entries; particle sources are passed as ParticlesToGrid.

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

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

class Options(stab_mat: Literal['M0', 'M0ad', 'Id', 'M0ad_withT'] = 'M0ad_withT', diffusion_mat: Literal['M1', 'M1perp', 'M1gyro'] = 'M1gyro', which_geometry: Literal['cylindrical', 'toroidal'] = 'cylindrical', x0: StencilVector | None = None, solver: Literal['pcg', 'cg'] = 'pcg', precond: Literal['MultiGrid', 'StiffnessPreconditioner', 'MassMatrixPreconditioner', 'PolarMassMatrixPreconditioner', None] | None = None, precond_params: dict | None = None, multigrid: MultiGridOptions | None = None, solver_params: SolverParameters | None = None, filter_params: dict[PICVariable | SPHVariable, FilterParameters] | None = None, enforce_compatibility: bool = False)

Bases: object

Configuration options for GyrokineticPoissonSolve.

Parameters:

• epsilon (float, default=1.0) – Gyrokinetic parameter

• Z (float, default=1.0) – Charge number of ions

• 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.

• which_geometry ({"cylindrical", "toroidal"}, default="cylindrical") – Geometry of the problem, determines the meaning of <phi> in the stabilization term: the average over eta3 (cylindrical) or over eta2 and eta3 (toroidal), i.e. stab_average of ImplicitDiffusion is "eta3" or "eta2 eta3". The average is weighted with stab_mat (for equilibrium density and temperature constant on the averaged surfaces, the flux-surface average with the volume element), which keeps the operator symmetric (WeightedAverageProjection).

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

• rho (FEECVariable or Callable or ParticlesToGrid 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.

  • ParticlesToGrid, describing a particle-to-grid (charge/current) deposition.

  • a list containing any mix of the entries above.

• 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.OptsDiffusionPrecond, default=None) – Name of the preconditioner, see ImplicitDiffusion. Requires solver="pcg". "StiffnessPreconditioner" approximates the operator with diffusion_mat, stab_mat and the average; "MultiGrid" is the geometric multigrid V-cycle, with the average re-built on every coarse level.

• precond_params (dict, default=None) – Keyword arguments passed to the constructor of the mass-matrix or stiffness preconditioner, see ImplicitDiffusion.

• multigrid (MultiGridOptions, default=None) – Options of the multigrid preconditioner (if precond="MultiGrid"). If None, defaults to MultiGridOptions().

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

• filter_params (dict[PICVariable | SPHVariable, FilterParameters], default=None) – If not None, specifies a filter to the accumulation of a specific variable.

Notes

GyrokineticPoissonSolve reuses ImplicitDiffusion internals by enforcing sigma_1 = 1 / (Z * epsilon**2), sigma_2 = 0.0, sigma_3 = 1 / epsilon and divide_by_dt = False in allocate.

property stab_average: str

Directions averaged in the stabilization, see which_geometry.

class struphy.propagators.curl_curl_solve.CurlCurlSolve(j: FEECVariable | tuple[Callable, Callable, Callable] | tuple[AccumulatorVector, Particles] | list | None = None, j_coeffs: float | list | None = None)

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 > 0\) is a scalar (the operator is symmetric positive definite). 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', x0: StencilVector | None = None, solver: Literal['pcg', 'cg'] = 'pcg', precond: Literal['MultiGrid', 'HiptmairXu', 'StiffnessPreconditioner', 'MassMatrixPreconditioner', 'PolarMassMatrixPreconditioner', None] | None = None, precond_params: dict | None = None, multigrid: MultiGridOptions | None = None, solver_params: SolverParameters | None = None)

Bases: object

Configuration options for CurlCurlSolve.

Parameters:

• sigma (float, default=1.0) – Coefficient of the stabilization/mass term on the left-hand side, curl.T @ diffusion_mat @ curl + sigma * stab_mat. Must be positive for the operator to be positive definite (values below 1e-14 are set to 1e-14).

• 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.

• 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.OptsCurlCurlPrecond, default=None) –

  Name of the preconditioner (requires solver="pcg"):

  • None: no preconditioner.

  • "MassMatrixPreconditioner": Kronecker approximation of the inverse of the 1-form mass matrix M1 (see struphy.feec.preconditioner).

  • "PolarMassMatrixPreconditioner": 2D x 1D Kronecker approximation of the inverse of M1 on polar splines.

  • "StiffnessPreconditioner": Kronecker approximation of the inverse of curl.T @ diffusion_mat @ curl + sigma * stab_mat (StiffnessPreconditioner with derivative="curl", mid=diffusion_mat and stab=stab_mat): block Jacobi with fast diagonalization, plus a correction on the gradients (kernel_correction, default True). With polar splines, this is a PolarStiffnessPreconditioner.

  • "HiptmairXu": auxiliary space preconditioner of Hiptmair and Xu (HiptmairXuPreconditioner with mid=diffusion_mat and stab=stab_mat): smoother plus corrections on the gradients and on vector-valued \(H^1\) splines.

  • "MultiGrid": geometric multigrid V-cycle (MultiGridPreconditioner) with the hybrid smoother of Hiptmair (see multigrid).

• precond_params (dict, default=None) – Keyword arguments passed to the constructor of the preconditioner, e.g. {"kernel_correction": False} for "StiffnessPreconditioner" or {"smoother": "jacobi"} for "HiptmairXu".

• multigrid (MultiGridOptions, default=None) – Options of the multigrid preconditioner (if precond="MultiGrid"). If None, defaults to MultiGridOptions(hiptmair=True, smoother_precond="jacobi"): without the Hiptmair smoother, multigrid is not robust for the curl-curl operator.

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

property sources: list[Vector | 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', 'PolarMassMatrixPreconditioner', 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', 'PolarMassMatrixPreconditioner', None] = 'MassMatrixPreconditioner', solver_params: SolverParameters | None = None, solver_M1: Literal['pcg', 'cg'] = 'pcg', precond_M1: Literal['MassMatrixPreconditioner', 'PolarMassMatrixPreconditioner', 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(b_field: FEECVariable | None = None)

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(u_space: Literal['Hcurl', 'Hdiv', 'H1vec'] = 'Hdiv', solver: Literal['pbicgstab', 'bicgstab', 'gmres'] = 'pbicgstab', precond: Literal['MassMatrixPreconditioner', 'PolarMassMatrixPreconditioner', None] = 'MassMatrixPreconditioner', solver_params: SolverParameters | None = None)

Bases: OptionsBase

Configuration options for Magnetosonic.

Parameters:

• 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', 'PolarMassMatrixPreconditioner', 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(*args, **kwargs)

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.

Note

Not available: this propagator has not been ported to the current Propagator, Accumulator and solver APIs, and the kernels hybrid_fA_Arelated and hybrid_curlA are missing.

class struphy.propagators.hall.Hall(epsilon_from: Species | None = None)

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', 'PolarMassMatrixPreconditioner', None] = 'MassMatrixPreconditioner', solver_params: SolverParameters | 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().

Nonlinear MHD solvers

class struphy.propagators.variational_density_evolve.VariationalDensityEvolve(s: FEECVariable | None = None)

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', 'PolarMassMatrixPreconditioner', None] = 'MassMatrixPreconditioner', solver_params: SolverParameters | None = None, nonlin_solver: NonlinearSolverParameters | 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 the mass-matrix solves (the non-symmetric Jacobian is always solved with “pbicgstab”).

• 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_entropy_evolve.VariationalEntropyEvolve(rho: FEECVariable | None = None)

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', 'PolarMassMatrixPreconditioner', None] = 'MassMatrixPreconditioner', solver_params: SolverParameters | None = None, nonlin_solver: NonlinearSolverParameters | 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.

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', 'PolarMassMatrixPreconditioner', 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', 'PolarMassMatrixPreconditioner', 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 (the non-symmetric Jacobian is always solved with “gmres”).

• 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(div_u: FEECVariable | None = None, u2: FEECVariable | None = None, pt3: FEECVariable | None = None, bt2: FEECVariable | None = None)

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', 'PolarMassMatrixPreconditioner', None] = 'MassMatrixPreconditioner', solver_params: SolverParameters | None = None, nonlin_solver: NonlinearSolverParameters | 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.

class struphy.propagators.variational_qb_evolve.VariationalQBEvolve(div_u: FEECVariable | None = None, u2: FEECVariable | None = None, qt3: FEECVariable | None = None, bt2: FEECVariable | None = None)

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', 'PolarMassMatrixPreconditioner', None] = 'MassMatrixPreconditioner', solver_params: SolverParameters | None = None, nonlin_solver: NonlinearSolverParameters | 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.

class struphy.propagators.variational_resistivity.VariationalResistivity(rho: FEECVariable | None = None, pt3: FEECVariable | None = None)

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', 'PolarMassMatrixPreconditioner', None] = 'MassMatrixPreconditioner', solver_params: SolverParameters | None = None, nonlin_solver: NonlinearSolverParameters | None = None, linearize_current: bool = False, eta: float = 0.0, eta_a: float = 0.0, fast: bool = False)

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="MassMatrixPreconditioner") – 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.

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

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

• fast (bool, default=False) – If True and model="linear_q", skip the thermodynamic (energy balance) update.

class struphy.propagators.variational_viscosity.VariationalViscosity(rho: FEECVariable | None = None, pt3: FEECVariable | None = None)

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', 'PolarMassMatrixPreconditioner', None] = 'MassMatrixPreconditioner', solver_params: SolverParameters | None = None, nonlin_solver: NonlinearSolverParameters | None = None, mu: float = 0.0, mu_a: float = 0.0, alpha: float = 0.0, fast: bool = False)

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="MassMatrixPreconditioner") – Preconditioner used in linear solves.

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

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

• 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.

• fast (bool, default=False) – If True and model="linear_q", skip the thermodynamic (energy balance) update.

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', 'PolarMassMatrixPreconditioner', 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', 'PolarMassMatrixPreconditioner', 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(rho: FEECVariable | ParticlesToGrid | Callable | list | None = None, rho_coeffs: float | list | None = None, diagnostic: FEECVariable | None = None)

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

Parameters:

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

  Source term(s) \(\rho_i\) on the right-hand side. Several entries can be combined by passing a list; each is weighted by the corresponding entry of rho_coeffs. Accepted entries are:

  • None: zero source.

  • FEECVariable in H1: its coefficients are mass-multiplied by \(M^0\) on the fly.

  • Callable \(\rho_i(\mathbf x)\), projected to H1 via L2Projector.

  • ParticlesToGrid, describing a

    particle-to-grid deposition. An AccumulatorVector is built from it at allocate() and re-accumulated at every call. Its accum_space selects what is deposited:

    • "H1": the kernel deposits a scalar density directly, used as-is.

    • "Hcurl": the kernel deposits a 1-form (vector) quantity \(\mathbf w\), e.g. \(n\mathbf u\); the propagator then applies the weak divergence \(-\mathbb G^\top\) to it before adding it to the right-hand side. This is used in the Chorin-projection scheme for example.

  • a list containing any mix of the entries above.

• 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, it has one entry per entry of rho (the neutralising background of a full-f accumulator gets the coefficient of its accumulator), or one entry per collected source. If None, all coefficients default to 1.0.

• diagnostic (FEECVariable, default=None) – If not None, updates at each call to the propagator, takes the value of the right-hand side. Otherwise does not provide diagnostic.

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', stab_average: str | None = None, diffusion_mat: Literal['M1', 'M1perp'] = 'M1', x0: StencilVector | None = None, solver: Literal['pcg', 'cg'] = 'pcg', precond: Literal['MultiGrid', 'StiffnessPreconditioner', 'MassMatrixPreconditioner', 'PolarMassMatrixPreconditioner', None] | None = None, precond_params: dict | None = None, multigrid: MultiGridOptions | None = None, solver_params: SolverParameters | None = None, enforce_compatibility: bool = False, filter_params: dict[PICVariable | SPHVariable, FilterParameters] | 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.

• stab_average (str, default=None) – If not None, the stabilization becomes stab_mat @ (I - A), where A is the average over the given logical directions, e.g. "eta3" or "eta2 eta3", weighted with stab_mat (WeightedAverageProjection; e.g. the flux-surface average with the volume element for "M0"), such that the operator stays symmetric. With precond="MultiGrid", the average is re-built on every coarse level.

• 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.

• 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.OptsDiffusionPrecond, default=None) –

  Name of the preconditioner:

  • None: no preconditioner.

  • "MultiGrid": geometric multigrid V-cycle (MultiGridPreconditioner, requires solver="pcg").

  • "MassMatrixPreconditioner": Kronecker approximation of the inverse of the 0-form mass matrix M0 (see struphy.feec.preconditioner).

  • "PolarMassMatrixPreconditioner": 2D x 1D Kronecker approximation of the inverse of M0 on polar splines.

  • "StiffnessPreconditioner": Kronecker approximation of the inverse of grad.T @ diffusion_mat @ grad + sigma * stab_mat (I - A) (StiffnessPreconditioner with mid=diffusion_mat, stab=stab_mat and stab_average), with sigma = sigma_1 (divided by dt if divide_by_dt=True). With polar splines, this is a PolarStiffnessPreconditioner.

• precond_params (dict, default=None) – Keyword arguments passed to the constructor of the mass-matrix or stiffness preconditioner, e.g. {"weights": "unit"} for "StiffnessPreconditioner" or {"dim_reduce": None} for "MassMatrixPreconditioner".

• multigrid (MultiGridOptions, default=None) – Options of the multigrid preconditioner (if precond="MultiGrid"). If None, defaults to MultiGridOptions(). Set nullspace="constants" for (nearly) singular problems, e.g. a periodic Poisson problem with tiny sigma_1.

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

• enforce_compatibility (bool, default=False) –

  If True and the diffusion operator is singular (no Dirichlet boundary conditions and no polar splines, i.e. the constant functions are in its kernel), the right-hand side is made compatible before the solve by removing its constant (“net charge”) component,

  \[\mathbf b \leftarrow \mathbf b - \frac{\mathbf 1^\top \mathbf b}{\mathbf 1^\top \mathbb M^0 \mathbf 1}\, \mathbb M^0 \mathbf 1\,,\]

  where \(\mathbf 1\) (all ones) is the coefficient vector of the constant function (partition of unity). This is equivalent to adding a uniform neutralizing background. It is meant for (nearly) singular problems where sigma_1 is only a small stabilization; then phi stays of order one instead of growing like (net charge) / sigma_1, which would destroy the accuracy of grad phi. For stab_mat="M0" and sigma_2 = 0 it changes phi only by a constant, for any sigma_1.

• filter_params (dict[PICVariable | SPHVariable, FilterParameters], default=None) – If not None, specifies a filter to the accumulation of a specific variable. Keyed by the pic_variable of the corresponding ParticlesToGrid source in rho.

Note

The accumulation kernel for a particle source is not configured here: it is carried by the ParticlesToGrid instance passed as (part of) rho.

property operator_is_singular: bool

True if the diffusion operator grad.T @ diffusion_mat @ grad is known to be singular with the constant functions (coefficient vector of all ones) in its kernel: no Dirichlet boundary conditions and tensor-product splines. Polar splines are conservatively treated as regular.

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("forward_euler") (Euler-Maruyama). The noise is additive and constant within a step, so multi-stage tableaux give statistically the same result at a higher cost.

• 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(u_tilde: FEECVariable | None = None)

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_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_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_v_in_force_field.PushVinForceField(force_field: FEECVariable | tuple[Callable] | None = None, potential: FEECVariable | Callable | None = None)

Bases: Propagator

Push the velocities according to

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

where \(\varepsilon \in \mathbb R\) is a constant species parameter. 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 F}^1(\boldsymbol \eta_p) \,,\]

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

class Variables

Bases: object

Container for variables advanced by PushVinForceField.

Variables:

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

class Options

Bases: OptionsBase

Configuration options for PushVinForceField.

Only particle variables are updated.

class struphy.propagators.push_vxb.PushVxB(b2_var: FEECVariable | None = None)

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

Bases: OptionsBase

Configuration options for PushVxB.

Parameters:

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

Only particle variables are updated.

class struphy.propagators.push_guiding_center_bx_estar.PushGuidingCenterBxEstar(phi: FEECVariable | None = None, b_tilde: FEECVariable | None = None)

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_dg_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(evaluate_e_field: bool = False, 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:

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

• 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(phi: FEECVariable | None = None, b_tilde: FEECVariable | None = None)

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(evaluate_e_field: bool = False, 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:

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

• 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 1N \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} \sum_p w_p \mathbb{\Lambda}^1 \cdot \left( DF^{-1} \mathbf{v}_p \right) \,,\end{split}\]

with the weights \(w_p = f_{1,p} / (N s_{0,p})\).

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} \mathbb{\Lambda}^1 \cdot \left( DF^{-1} \mathbf{v}_p \right) \, , \\[2mm] \mathbb{W} &= \frac 1N \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', 'PolarMassMatrixPreconditioner', 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', 'PolarMassMatrixPreconditioner', 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', 'PolarMassMatrixPreconditioner', 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.

Gyrokinetic-Maxwell coupling

Propagators of the gauge-invariant gyrokinetic Maxwell model of Remmerswaal, Hatzky and Sonnendrücker (J. Plasma Phys. 91, E127, 2025), see GyrokineticMaxwell.

The field unknowns are \(\mathbf E \in H(\textnormal{curl})\) and \(\mathbf B \in H(\textnormal{div})\) (the model is gauge invariant, hence formulated in the fields rather than in the potentials). All propagators share the dielectric matrix dielectric_matrix() (vacuum permittivity plus polarization), the magnetizing matrix magnetizing_matrix() (pressure anisotropy of the background) and, for a background with parallel flow, the flow operator flow_operator(). The background moments of each species (SpeciesBackground) are taken from its kinetic background.

Finite Larmor radius effects (GyroAverage) replace the point evaluations of the fields at the gyrocentres by gyro-averages (ring) and disc averages, with the couplings to the electric field applied matrix-free (GyroParticleOperator), since the gyro-averaged particle matrices are not stencil matrices.

The propagators work with tensor-product FE coefficients (no polar splines).

class struphy.propagators.gyrokinetic_maxwell.SpeciesBackground(var: PICVariable, kinetic: bool = True)

Bases: object

The background moments \(n_{0,s}\), \(u_{0,\parallel,s}\), \(v_{\textnormal{th},\parallel,s}\), \(v_{\textnormal{th},\perp,s}\) of one gyrocentre species on the logical cube, from its kinetic background (a GyroMaxwellian2D or a sum of them), or from the fluid equilibrium density (kinetic=False: \(n_0\), no flow, isotropic).

The normalized species pressures are \(p_{\parallel,s} = \alpha_s^2 n_s v_{\textnormal{th},\parallel,s}^2\) and \(p_{\perp,s} = \alpha_s^2 n_s v_{\textnormal{th},\perp,s}^2\) (light units, Eq. (5.56) of the paper).

Parameters:

• var (PICVariable) – The species’ variable ("Particles5D"), with its kinetic background(s) added.

• kinetic (bool) – Whether to take the moments from the kinetic background (default) or the density from the fluid equilibrium.

n(*etas, equil=None)

Background density \(n_{0,s}\) (units of the density unit).

u_par(*etas)

Background parallel flow \(u_{0,\parallel,s}\) (units of c).

vth_par2(*etas)

Squared parallel thermal speed.

vth_perp2(*etas)

Squared perpendicular thermal speed.

class struphy.propagators.gyrokinetic_maxwell.BackgroundWeights(backgrounds: tuple[SpeciesBackground, ...], equil)

Bases: object

The coefficients of the field operators from the species backgrounds (Eq. (2.2a) of the paper, light units):

\[\textnormal{polarization} = \sum_s \alpha_s^2 \frac{n_{0,s}}{B_0^2}\,, \qquad \textnormal{flow} = \sum_s \alpha_s^2 \frac{n_{0,s} u_{0,\parallel,s}}{B_0^2}\,, \qquad \textnormal{anisotropy} = \frac{p_{0,\parallel} - p_{0,\perp}}{B_0^2} = \sum_s \alpha_s^2 \frac{n_{0,s} (v_{\textnormal{th},\parallel,s}^2 - v_{\textnormal{th},\perp,s}^2)}{B_0^2}\,,\]

as callables of the logical coordinates (for create_weighted_mass()).

Parameters:

• backgrounds (tuple[SpeciesBackground]) – The species backgrounds.

• equil (FluidEquilibriumWithB) – The fluid equilibrium (for \(B_0\) and, with kinetic=False backgrounds, \(n_0\)).

property has_flow: bool

Whether any kinetic species background has a parallel flow (u_para != 0 on a sample grid).

property has_anisotropy: bool

Whether any kinetic species background is anisotropic (vth_para != vth_perp on a sample grid).

skin(e1, e2, e3)

The continuum limit \(\sum_s \alpha_s^2 n_{0,s}/\varepsilon_s^2\) of the skin term of the implicit parallel couplings.

struphy.propagators.gyrokinetic_maxwell.parallel_flow_from_equilibrium(equil, alpha: float, epsilon: float, density: Callable, fraction: float = 1.0)

The parallel flow \(u_{0,\parallel,s}\) a species needs to carry a fraction of the equilibrium current, Eq. (5.62) of the paper with isotropic background pressure:

\[\frac{\alpha_s^2}{\varepsilon_s} n_{0,s}\, u_{0,\parallel,s} = f_s\, \mathbf b_0 \cdot (\nabla \times \mathbf B_0)\,,\]

such that, with \(\sum_s f_s = 1\), the gyrocentre currents of the backgrounds balance \(\nabla \times \mathbf B_0\) in Ampère’s law (background_current=True of the model).

Parameters:

• equil (FluidEquilibriumWithB) – The equilibrium (provides \(\nabla \times \mathbf B_0\) and \(\mathbf b_0\)).

• alpha (float) – The species parameters.

• epsilon (float) – The species parameters.

• density (Callable) – The background density \(n_{0,s}\) of the species on the logical cube.

• fraction (float) – The fraction \(f_s\) of the parallel current carried by the species.

Returns:

u_par(e1, e2, e3), to be used as u_para=(u_par, None) of a GyroMaxwellian2D.

Return type:

Callable

struphy.propagators.gyrokinetic_maxwell.dielectric_matrix(mass_ops: WeightedMassOperators, polarization) → WeightedMassOperator

Dielectric matrix of the gyrokinetic Maxwell model,

\[\mathbb A = \int \vec{\Lambda}^1 \left[ \left(1 + \sum_s \alpha_s^2 \frac{n_{0,s}}{B_0^2}\right) G^{-1} - \sum_s \alpha_s^2 \frac{n_{0,s}}{B_0^2} \hat{\mathbf b}_0^v (\hat{\mathbf b}_0^v)^\top \right] \vec{\Lambda}^{1\top} \sqrt g\, \textnormal d \boldsymbol \eta\,,\]

i.e. the weak form of \(\mathbf I + \sum_s \alpha_s^2 \frac{n_{0,s}}{B_0^2} \Pi_\perp\) (vacuum permittivity plus the polarization of the second-order Hamiltonian), assembled as one WeightedMassOperator.

Parameters:

• mass_ops (WeightedMassOperators) – The weighted mass operators of the simulation (provide domain and equilibrium).

• polarization (float | Callable) – Either the constant \(\sum_s \alpha_s^2\) (the density is then the fluid equilibrium \(n_0\) for every species), or the callable \(\sum_s \alpha_s^2 n_{0,s}/B_0^2\) of the logical coordinates (BackgroundWeights.polarization()).

Return type:

The assembled operator; it is cached in mass_ops and shared between the propagators.

struphy.propagators.gyrokinetic_maxwell.polarization_matrix(mass_ops: WeightedMassOperators, polarization) → WeightedMassOperator

Polarization matrix \(\mathbb P = \mathbb A - \mathbb M^1\), the weak form of \(\sum_s \alpha_s^2 \frac{n_{0,s}}{B_0^2} \Pi_\perp\) on 1-forms (the coefficient \(\mathcal C(1) \Pi_\perp\) of Eq. (5.47) of the paper), assembled as one WeightedMassOperator; same polarization argument as dielectric_matrix().

struphy.propagators.gyrokinetic_maxwell.skin_preconditioner_matrix(mass_ops: WeightedMassOperators, polarization, skin, scale: float, vacuum: bool)

A weighted mass matrix approximating the operator of the implicit parallel coupling, \(\mathbb A + c\, \mathbb S\) (or \(\mathbb P_L + c\, \mathbb S\) of the Darwin model), by the weight \([v + s(\boldsymbol \eta)] G^{-1} - s(\boldsymbol \eta)\, \hat{\mathbf b}_0^v (\hat{\mathbf b}_0^v)^\top + c\, \kappa(\boldsymbol \eta)\, \hat{\mathbf b}_0^v (\hat{\mathbf b}_0^v)^\top\) with \(\kappa = \sum_s \frac{\alpha_s^2}{\varepsilon_s^2} n_{0,s}\) (the continuum limit of the skin term), for the MassMatrixPreconditioner.

Parameters:

• mass_ops (WeightedMassOperators)

• polarization (float | Callable) – As in dielectric_matrix().

• skin (Callable) – \(\kappa(\boldsymbol \eta)\) on the logical cube.

• scale (float) – The constant \(c\) (e.g. \(\Delta t^2/4\)).

• vacuum (bool) – Whether the vacuum permittivity is included (Maxwell) or not (Darwin).

struphy.propagators.gyrokinetic_maxwell.magnetizing_matrix(mass_ops: WeightedMassOperators, anisotropy=None) → WeightedMassOperator

Magnetizing matrix of the gyrokinetic Maxwell model: the weak form of \(\mathbf H = \mathbf B - \frac{p_{0,\parallel} - p_{0,\perp}}{B_0^2} \mathbf B_\perp\) on 2-forms,

\[\mathbb M^2_H = \int \vec{\Lambda}^2 \left[(1 - c)\, G + c\, \hat{\mathbf b}^1_0 (\hat{\mathbf b}^1_0)^\top\right] \vec{\Lambda}^{2\top} \frac{1}{\sqrt g}\, \textnormal d \boldsymbol \eta\,, \qquad c = \frac{p_{0,\parallel} - p_{0,\perp}}{B_0^2}\,,\]

which is \(\mathbb M^2\) for an isotropic background (anisotropy=None).

struphy.propagators.gyrokinetic_maxwell.flow_operator(mass_ops: WeightedMassOperators, flow) → WeightedMassOperator

The weighted mass operator \(\mathbb R: V_h^1 \to V_h^2\) of the background parallel flow,

\[\int \mathbf F \cdot \left(s\, \mathbf b_0 \times \mathbf E\right) \textnormal d \mathbf x = \mathbf f^\top \mathbb R\, \mathbf e\,, \qquad s = \sum_s \alpha_s^2 \frac{n_{0,s} u_{0,\parallel,s}}{B_0^2}\,,\]

for a 2-form test function \(\mathbf F\). The flow terms of Ampère’s law (2.2a), the magnetisation \(\nabla \times (s\, \mathbf b_0 \times \mathbf E)\) and the polarization current \(s\, \mathbf b_0 \times (\nabla \times \mathbf E)\), are the anti-symmetric operator \(\mathbb K = \mathbb C^\top \mathbb R - \mathbb R^\top \mathbb C\).

class struphy.propagators.gyrokinetic_maxwell.GyroAverage(n_gyro: int = 0, n_disc: int = 2)

Bases: object

Discretization of the gyro-average (ring) and the disc average of the fields at the gyrocentres.

Parameters:

• n_gyro (int) – Number of points on the gyro-ring; 0 is the zero-Larmor-radius limit (point evaluation).

• n_disc (int) – Number of Gauss-Legendre points in the radial direction of the disc average; 0 approximates the disc average by the gyro-average.

property flr: bool

Whether finite-Larmor-radius averages are used.

class struphy.propagators.gyrokinetic_maxwell.GyroParticleOperator(particles, derham, args_domain, kernel, args: tuple, scale: float = 1.0, symmetric=True)

Bases: LinearOperator

Matrix-free linear operator \(V^1_h \to V^1_h\) of the gyro-averaged coupling of the markers of one species to a 1-form field,

\[\mathbf e \mapsto c \sum_p w_p \langle \vec\Lambda^1 \rangle_p\, \mathbb K_p\, \langle \vec\Lambda^{1\top} \rangle_p \mathbf e\,,\]

with a per-marker \(3 \times 3\) matrix \(\mathbb K_p\) defined by the kernel (e.g. \(\mathbf c_p \mathbf c_p^\top\) for the parallel coupling, the E x B rotation for the drift coupling). Each product evaluates the field at the gyro-points of every marker and deposits there again, so the operator is applied inside an iterative solver (the gyro-averaged particle matrices are not stencil matrices).

Parameters:

• particles (Particles) – The markers.

• derham (Derham) – The discrete de Rham complex (tensor-product spaces).

• args_domain (DomainArguments) – Mapping arguments.

• kernel (Kernel) – A kernel with the signature (args_markers, args_derham, args_domain, vec1, vec2, vec3, e1, e2, e3, *args).

• args (tuple) – The remaining kernel arguments.

• scale (float) – The constant \(c\).

• symmetric (bool) – Whether the operator is symmetric (else it is anti-symmetric), for transpose().

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.

transpose(conjugate=False)

Transpose the LinearOperator .

If conjugate is True, return the Hermitian transpose.

tosparse()

Assemble the global matrix of the linear operator as a scipy.sparse matrix.

Default implementation calling the generic LinearOperator.toarray with is_sparse=True; see there for cost and parallel behavior. Subclasses with an explicit matrix representation should override it.

Parameters:

format (str) – One of ‘csr’ (default), ‘csc’, ‘bsr’, ‘lil’, ‘dok’, ‘coo’ or ‘dia’.

Returns:

The global matrix of shape self.shape, identical on all ranks.

Return type:

scipy.sparse matrix

toarray(**kwargs)

Assemble the global matrix of the linear operator column by column.

Column j is computed as self.dot(e_j), where e_j is the j-th canonical basis vector of the domain, in the global numbering of Vector.toarray() (C-ordered within each StencilVector, blocks concatenated). Since only dot is used, this works for any linear operator, including matrix-free ones, whose domain is a StencilVectorSpace or a (possibly nested) BlockVectorSpace thereof.

The cost is one call to dot per global degree of freedom of the domain, hence this default is meant for testing and small problems. Subclasses with an explicit matrix representation should override it.

In parallel, all ranks call dot collectively, each rank fills the rows it owns, and every rank receives the full matrix.

Parameters:

• out (numpy.ndarray, optional) – Dense array of shape self.shape into which the result is written in place. Must be None if is_sparse is True.

• is_sparse (bool) – If True, return a scipy.sparse matrix, otherwise a dense numpy.ndarray.

• format (str) – Sparse format, only used if is_sparse is True: one of ‘csr’ (default), ‘csc’, ‘bsr’, ‘lil’, ‘dok’, ‘coo’ or ‘dia’.

Returns:

The global matrix of shape self.shape, identical on all ranks.

Return type:

numpy.ndarray or scipy.sparse matrix

dot(v, out=None)

Apply the LinearOperator self to the Vector v.

The result is written to the Vector out, if provided.

Parameters:

• v (Vector) – The vector to which the linear operator (self) is applied. It must belong to the domain of self.

• out (Vector) – The vector in which the result of the operation is stored. It must belong to the codomain of self. If out is None, a new vector is created and returned.

Returns:

The result of the operation. If out is None, a new vector is returned. Otherwise, the result is stored in out and out is returned.

Return type:

Vector

class struphy.propagators.gyrokinetic_maxwell.MaxwellGyroMagnetisation(species: tuple[PICVariable, ...] = (), background: BackgroundWeights | None = None, polarization: float | None = None, background_current: bool = False, gyro: GyroAverage | None = None)

Bases: Propagator, _PerturbedFieldMixin

FEEC discretization of Maxwell’s equations of the gyrokinetic Maxwell model (Eq. (2.2) of the paper): find \(\mathbf E \in H(\textnormal{curl})\) and \(\mathbf B \in H(\textnormal{div})\) such that

\[\begin{split}\int_{\Omega} \left(\mathbf I + \sum_s \alpha_s^2 \frac{n_{0,s}}{B_0^2} \Pi_\perp \right)\frac{\partial \mathbf E}{\partial t} \cdot \mathbf F \, \textrm d \mathbf x - \int_{\Omega} \left(\mathbf B - \frac{p_{0,\parallel} - p_{0,\perp}}{B_0^2} \mathbf B_\perp + \sum_s \alpha_s^2 \frac{n_{0,s} u_{0,\parallel,s}}{B_0^2} \mathbf b_0 \times \mathbf E + \sum_s \alpha_s^2 \int f_s\, \mu\, \mathbf b_0 \,\textrm d v_\parallel \textrm d \mu \right) \cdot \nabla \times \mathbf F \,\textrm d \mathbf x - \int_\Omega \sum_s \alpha_s^2 \frac{n_{0,s} u_{0,\parallel,s}}{B_0^2} \mathbf b_0 \times (\nabla \times \mathbf E) \cdot \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}\]

where \(\mathbf B = \mathbf B_0 + \tilde{\mathbf B}\) (with background_current, else \(\tilde{\mathbf B}\) only, which assumes that the gyrocentre currents of the backgrounds balance \(\nabla \times \mathbf B_0\)). The curl term with \(\mu \mathbf b_0\) is the curl of the intrinsic magnetisation \(-\mu \mathbf b_0\) of the gyrocentres. In matrix form (dielectric_matrix() \(\mathbb A\), magnetizing_matrix() \(\mathbb M^2_H\), flow_operator() \(\mathbb K = \mathbb C^\top \mathbb R - \mathbb R^\top \mathbb C\)):

\[\mathbb A \dot{\mathbf e} = \mathbb C^\top \mathbb M^2_H (\mathbf b + \mathbf b_0) + \mathbb K \mathbf e + \mathbb C^\top (\mathcal P^b)^\top \sum_s \alpha_s^2 \mathbf m_s\,, \qquad \dot{\mathbf b} = - \mathbb C \mathbf e\,, \qquad \mathbf m_s = \sum_p w_p \mu_p \langle\langle \vec \Lambda^0 \rangle\rangle_p\,,\]

where \(\mathcal P^b\) is the basis projection PB (\(\mathbf b \mapsto \Pi^0(\mathbf b_0 \cdot \mathbf B)\)) and \(\langle\langle \cdot \rangle\rangle\) the disc average (a point evaluation in the ZLR limit), such that the magnetisation source is the exact transpose of the map \(\mathbf b \mapsto \sum_p w_p \mu_p \langle\langle B_\parallel \rangle\rangle(\boldsymbol \eta_p)\) of the particle energy: the system is Hamiltonian and the implicit mid-point rule conserves the energy exactly (for fixed markers); \(\mathbb K\) is anti-symmetric and exchanges no energy.

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

\[\left(\mathbb A + \frac{\Delta t^2}{4} \mathbb C^\top \mathbb M^2_H \mathbb C - \frac{\Delta t}{2} \mathbb K\right) \mathbf e^{n+1} = \left(\mathbb A - \frac{\Delta t^2}{4} \mathbb C^\top \mathbb M^2_H \mathbb C + \frac{\Delta t}{2} \mathbb K\right) \mathbf e^{n} + \Delta t\, \mathbb C^\top \mathbb M^2_H \mathbf b^n + \Delta t\, \mathbf s\,, \qquad \mathbf b^{n+1} = \mathbf b^n - \frac{\Delta t}{2} \mathbb C (\mathbf e^n + \mathbf e^{n+1})\,,\]

solved with pcg (pbicgstab with the flow operator).

class Variables

Bases: object

Container for the variables advanced by MaxwellGyroMagnetisation.

Variables:

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

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

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

Bases: OptionsBase

Configuration options for MaxwellGyroMagnetisation.

Parameters:

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

• solver_flow (LiteralOptions.OptsGenSolver, default="pbicgstab") – General iterative solver (with background flow, the system is not symmetric).

• precond (LiteralOptions.OptsMassPrecond, default="MassMatrixPreconditioner") – Preconditioner of the dielectric matrix.

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

class struphy.propagators.gyrokinetic_maxwell.GyroEparCoupling(b_tilde: FEECVariable | None = None, polarization=0.0, gyro: GyroAverage | None = None)

Bases: Propagator, _PerturbedFieldMixin

Parallel coupling of one gyrocentre species to the electric field: find \(\mathbf E \in H(\textnormal{curl})\) and \(f\) such that

\[\begin{split}\int_{\Omega} \left(\mathbf I + \sum_s \alpha_s^2 \frac{n_{0,s}}{B_0^2} \Pi_\perp \right)\frac{\partial \mathbf{E}}{\partial t} \cdot \mathbf{F}\,\textrm d \mathbf x &= - \frac{\alpha^2}{\varepsilon} \int_{\Omega} \int f\, v_\parallel \frac{\mathbf B^*}{B^*_\parallel} \cdot \langle \mathbf{F} \rangle \,\textrm d v_\parallel \textrm d \mu \,\textrm d \mathbf x \qquad \forall \, \mathbf{F} \in H(\textnormal{curl}) \\[2mm] \frac{\partial f}{\partial t} + \frac{1}{\varepsilon}\, \frac{\mathbf B^*}{B^*_\parallel} \cdot \langle \mathbf{E} \rangle\, \frac{\partial f}{\partial v_\parallel} &= 0\,,\end{split}\]

where \(\mathbf B^* = \mathbf B_0 + \langle \tilde{\mathbf B} \rangle + \varepsilon v_\parallel \nabla \times \mathbf b_0\) (the FLR-corrected field \(\mathbf B_1^\star = \langle \mathbf B_1 \rangle\) of the paper, Eq. (7.24)) and \(B^*_\parallel = \mathbf b_0 \cdot \mathbf B^*\); \(\langle \cdot \rangle\) is the gyro-average (a point evaluation in the ZLR limit). For each marker, in logical coordinates,

\[\mathbb A \dot{\mathbf e} = -\frac{\alpha^2}{\varepsilon} \sum_p w_p v_{\parallel,p}\, \langle \vec \Lambda^1 \rangle_p \cdot \mathbf c_p\,, \qquad \dot v_{\parallel,p} = \frac{1}{\varepsilon}\, \mathbf c_p \cdot \langle \vec \Lambda^{1\top} \rangle_p\, \mathbf e\,, \qquad \mathbf c_p = \frac{\hat{\mathbf B}^{*2}}{\sqrt g\, \hat B^*_\parallel}(\boldsymbol \eta_p)\,,\]

with \(\mathbf c_p\) frozen at \(v_{\parallel,p}^n\), such that the system is linear in \((\mathbf e, v_\parallel)\).

Time discretization: Crank-Nicolson (implicit mid-point), system size reduction via the Schur complement

\[\left(\mathbb A + \frac{\Delta t^2 \alpha^2}{4 \varepsilon^2} \mathbb S\right) \mathbf e^{n+1} = \left(\mathbb A - \frac{\Delta t^2 \alpha^2}{4 \varepsilon^2} \mathbb S\right) \mathbf e^{n} - \Delta t \frac{\alpha^2}{\varepsilon} \mathbf j\,, \qquad \mathbb S = \sum_p w_p \langle \vec \Lambda^1 \rangle_p \mathbf c_p \mathbf c_p^\top \langle \vec \Lambda^{1\top} \rangle_p\,, \qquad \mathbf j = \sum_p w_p v^n_{\parallel,p} \mathbf c_p \langle \vec \Lambda^1 \rangle_p\,,\]

with the stencil matrix \(\mathbb S\) and SchurSolver in the ZLR limit (same structure as VlasovAmpereCoupling), and the matrix-free GyroParticleOperator with finite Larmor radius. The accumulation and push kernels evaluate the identical \(\mathbf c_p\), which gives exact energy conservation.

class Variables

Bases: object

Container for the variables advanced by GyroEparCoupling.

Variables:

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

• var (PICVariable) – Gyrocentre particle variable in "Particles5D" space.

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

Bases: OptionsBase

Configuration options for GyroEparCoupling.

Parameters:

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

• precond (LiteralOptions.OptsMassPrecond, default="MassMatrixPreconditioner") – Preconditioner of the dielectric matrix.

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

class struphy.propagators.gyrokinetic_maxwell.PushGuidingCenterParallelFullB(b_tilde: FEECVariable)

Bases: PushGuidingCenterParallel

Same as PushGuidingCenterParallel (without electric potential), but with the magnetic perturbation in the full \(\mathbf B^*\):

\[\mathbf B^* = \mathbf B_0 + \tilde{\mathbf B} + \varepsilon v_\parallel \nabla \times \mathbf b_0\,,\qquad B^*_\parallel = |B_0| + \Pi^0(\mathbf b_0 \cdot \tilde{\mathbf B}) + \varepsilon v_\parallel\, \mathbf b_0 \cdot \nabla \times \mathbf b_0\,,\]

which adds the magnetic flutter (streaming along the perturbed field lines) to the parallel motion (the parent class uses \(\tilde{\mathbf B}\) only in \(B^*_\parallel\) and in \(\mathbf E^* = -\varepsilon \mu \nabla B_\parallel\)). The vector b2 of the projected equilibrium seen by the kernels is replaced by a persistent vector holding \(\mathbf B_0 + \tilde{\mathbf B}\) (tensor-product coefficients), updated before every push. Zero-Larmor-radius limit; see GyroParallelPush for finite Larmor radius.

property projected_equil

Fluid equilibrium projected on 3d Derham sequence with commuting projectors.

class struphy.propagators.gyrokinetic_maxwell.GyroParallelPush(b_tilde: FEECVariable | None = None, gyro: GyroAverage | None = None)

Bases: Propagator, _PerturbedFieldMixin

Parallel gyrocentre motion with finite Larmor radius: 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} &= -\mu_p \frac{\mathbf B^*}{B^*_\parallel} \cdot \nabla \left(B_0 + \langle\langle \tilde B_\parallel \rangle\rangle\right) (\mathbf X_p(t)) \,, \end{aligned} \right.\end{split}\]

with \(\mathbf B^* = \mathbf B_0 + \langle \tilde{\mathbf B} \rangle + \varepsilon v_\parallel \nabla \times \mathbf b_0\), \(B^*_\parallel = |B_0| + \langle \tilde B_\parallel \rangle + \varepsilon v_\parallel \mathbf b_0 \cdot \nabla \times \mathbf b_0\) (Eqs. (5.23)-(5.24) of the paper with \(\mathbf B_1^\star = \langle \mathbf B_1 \rangle\) and the disc-averaged parallel perturbation in the mirror force), by the first-order discrete-gradient (Gonzalez) scheme of PushGuidingCenterParallel, which conserves the energy \(v_\parallel^2/2 + \mu (B_0 + \langle\langle \tilde B_\parallel \rangle\rangle)\) of each marker exactly.

class Variables

Bases: object

Container for the variables advanced by GyroParallelPush.

Variables:

ions (PICVariable) – Gyrocentre particle variable in "Particles5D" space.

class Options(maxiter: int = 20, tol: float = 1e-07, mpi_sort: Literal['each', 'last', None] = 'each')

Bases: OptionsBase

Configuration options for GyroParallelPush.

Parameters:

• maxiter (int, default=20) – Maximum number of fixed-point iterations of the discrete-gradient scheme.

• tol (float, default=1e-7) – Convergence tolerance of the iterations.

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

class struphy.propagators.gyrokinetic_maxwell.GyroDriftCoupling(species: tuple[PICVariable, ...] = (), b_tilde: FEECVariable | None = None, polarization=0.0, gyro: GyroAverage | None = None)

Bases: Propagator, _PerturbedFieldMixin

Perpendicular (drift) coupling of all gyrocentre species to the electric field: find \(\mathbf E \in H(\textnormal{curl})\) and \(f_s\) such that

\[\begin{split}\int_{\Omega} \left(\mathbf I + \sum_s \alpha_s^2 \frac{n_{0,s}}{B_0^2} \Pi_\perp \right)\frac{\partial \mathbf{E}}{\partial t} \cdot \mathbf{F}\,\textrm d \mathbf x &= - \sum_s \frac{\alpha_s^2}{\varepsilon_s} \int_{\Omega} \int f_s\, \frac{\mathbf E^* \times \mathbf b_0}{B^*_\parallel} \cdot \langle \mathbf{F} \rangle \,\textrm d v_\parallel \textrm d \mu \,\textrm d \mathbf x \qquad \forall \, \mathbf{F} \in H(\textnormal{curl}) \\[2mm] \frac{\partial f_s}{\partial t} + \frac{\mathbf E^* \times \mathbf b_0}{B^*_\parallel} \cdot \nabla f_s &= 0\,, \qquad \mathbf E^* = \langle \mathbf E \rangle - \varepsilon_s \mu \nabla \left(B_0 + \langle\langle \tilde B_\parallel \rangle\rangle\right)\,,\end{split}\]

(the gyro- and disc averages are point evaluations in the ZLR limit). With the markers at \(\boldsymbol \eta^n\), the \(E \times B\) current is an anti-symmetric operator \(\mathbb Q = \sum_s \frac{\alpha_s^2}{\varepsilon_s} \sum_p w_p \langle \vec \Lambda^1 \rangle_p \frac{\epsilon_{\mu\nu\lambda} \hat b_{0,\lambda}}{\sqrt g \hat B^*_\parallel} \langle \vec \Lambda^{1\top} \rangle_p\) (hence exchanges no energy, \(\mathbf E \cdot (\mathbf E \times \mathbf b_0) = 0\)) and the \(\nabla B\) current a vector \(\mathbf j = \sum_s \frac{\alpha_s^2}{\varepsilon_s} \sum_p w_p \dot{\boldsymbol \eta}_{\nabla B,p} \langle \vec \Lambda^1 \rangle_p\):

\[\left(\mathbb A + \frac{\Delta t}{2} \mathbb Q \right) \mathbf e^{n+1} = \left(\mathbb A - \frac{\Delta t}{2} \mathbb Q \right) \mathbf e^{n} - \Delta t\, \mathbf j\,,\]

with a stencil matrix \(\mathbb Q\) in the ZLR limit and the matrix-free GyroParticleOperator with finite Larmor radius. The markers are then pushed with an explicit Runge-Kutta scheme in the mid-point field \((\mathbf e^n + \mathbf e^{n+1})/2\). The energy exchange of the \(\nabla B\) current with the particle energy \(\mu B_\parallel\) is consistent to second order.

All species are handled in one propagator: a single species’ \(E \times B\) current, balanced only by the polarization, gives a spurious stiff mode near the cyclotron frequency.

class Variables

Bases: object

Container for the variables advanced by GyroDriftCoupling.

Variables:

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

class Options(solver: Literal['pbicgstab', 'bicgstab', 'gmres'] = 'pbicgstab', precond: Literal['MassMatrixPreconditioner', 'PolarMassMatrixPreconditioner', None] = 'MassMatrixPreconditioner', solver_params: SolverParameters | None = None, butcher: ButcherTableau | None = None, mpi_sort: Literal['each', 'last', None] = 'last')

Bases: OptionsBase

Configuration options for GyroDriftCoupling.

Parameters:

• solver (LiteralOptions.OptsGenSolver, default="pbicgstab") – General iterative solver for the (non-symmetric) linear system.

• precond (LiteralOptions.OptsMassPrecond, default="MassMatrixPreconditioner") – Preconditioner of the dielectric matrix.

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

• butcher (ButcherTableau, default=None) – Butcher tableau of the explicit marker push. If None, defaults to ButcherTableau(algo="rk2").

• mpi_sort (LiteralOptions.OptsMPIsort, default="last") – MPI sorting policy for particle exchange (the markers are sorted after the full step).

Propagators of the quasi-neutral gyrokinetic Darwin model of Remmerswaal, Hatzky and Sonnendrücker (J. Plasma Phys. 91, E127, 2025, Sec. 6), see GyrokineticDarwin.

The unknowns are the potentials \(\phi_1 \in H^1\) and \(\mathbf A_1 \in H(\textnormal{curl})\) in the perpendicular Coulomb gauge (Eq. (6.12)), with the Lagrange multiplier \(\lambda \in H^1\) of the gauge condition. The field operators (polarization, magnetizing and flow operators, the gyro-averaged particle couplings) are those of struphy.propagators.gyrokinetic_maxwell.

class struphy.propagators.gyrokinetic_darwin.DarwinLongitudinalOperator(derham, mass_ops, P, directions: tuple[int, ...], stab_eps: float, tol: float, maxiter: int, perp_polarization: bool = True)

Bases: LinearOperator

The longitudinal polarization operator of the gyrokinetic Darwin model on 1-forms,

\[\mathbb P_L = \mathbb P \mathbb G \left(\mathbb G^\top \mathbb P \mathbb G + \sigma \mathbb M^0 \Pi_{\textnormal{avg}}\right)^{-1} \mathbb G^\top \mathbb P\,,\]

the weak form of the Darwin displacement \(\mathbf D^{\textnormal{Dar}} = \Pi_{L,\perp}(\mathcal C(1) \mathbf E_\perp)\) (Eqs. (6.6), (6.9a) of the paper): the \(\mathbb P\)-orthogonal projection of a 1-form on gradients, applied matrix-free with an inner (preconditioned) conjugate-gradient solve of the perpendicular Laplacian. The average projection \(\Pi_{\textnormal{avg}}\) over the perpendicular directions regularizes the kernel of the perpendicular Laplacian (the \(k_\perp = 0\) modes, which carry no longitudinal displacement). The operator is symmetric and positive semi-definite.

For the perpendicular fields that are constant in the perpendicular directions (\(k_\perp = 0\), including the uniform field) the longitudinal projection is degenerate: they are both curl- and divergence-free in the perpendicular plane, Eq. (6.6) assigns them no longitudinal part and the Darwin model gives them no inertia (a \(k_\perp = 0\) shear Alfvén wave would have no polarization current, and the net E x B current of a non-neutral marker ensemble no equation at all). With perp_polarization=True the full (transversal) polarization of these modes is retained, \(\mathbb P_L + \Pi_\perp \mathbb P \Pi_\perp\) with \(\Pi_\perp\) the mean of the coefficients of each component over the averaged directions (a symmetric projection, so the operator stays symmetric positive semi-definite and the energy \(\frac12 \mathbf e^\top \mathbb P_L \mathbf e\) includes these modes); their shear Alfvén frequency is then the \(k_\perp \to 0\) limit \(\omega = k_\parallel / \sqrt{\sum_s \alpha_s^2 n_{0,s}/B_0^2}\).

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 laplacian

The regularized perpendicular Laplacian \(\mathbb G^\top \mathbb P \mathbb G + \sigma \mathbb M^0 \Pi_{\textnormal{avg}}\).

property inverse_laplacian

Its inverse (an iterative solver).

transpose(conjugate=False)

Transpose the LinearOperator .

If conjugate is True, return the Hermitian transpose.

tosparse()

Assemble the global matrix of the linear operator as a scipy.sparse matrix.

Default implementation calling the generic LinearOperator.toarray with is_sparse=True; see there for cost and parallel behavior. Subclasses with an explicit matrix representation should override it.

Parameters:

format (str) – One of ‘csr’ (default), ‘csc’, ‘bsr’, ‘lil’, ‘dok’, ‘coo’ or ‘dia’.

Returns:

The global matrix of shape self.shape, identical on all ranks.

Return type:

scipy.sparse matrix

potential(v, out=None)

The potential \(\phi = (\mathbb G^\top \mathbb P \mathbb G + \ldots)^{-1} \mathbb G^\top \mathbb P \mathbf v\) of a 1-form v.

solve_laplacian(rhs, out=None)

Solve the regularized perpendicular Laplacian with the tolerance relative to the norm of rhs.

remove_average(phi)

Remove (in place) the part of the 0-form phi that is constant in the averaged directions, i.e. \(\phi \leftarrow (\mathbb I - \Pi_{\textnormal{avg}}) \phi\): the kernel of the perpendicular Laplacian, whose gradient carries no longitudinal displacement. Returns the norm of the removed part.

perpendicular_average(v, out=None)

The projection \(\Pi_\perp \mathbf v\) on the modes constant in the averaged directions: each component is replaced by the mean of its coefficients over these directions (an orthogonal projection).

longitudinal_part(v, out=None)

The longitudinal part \(\Pi_L \mathbf v = \mathbb G \phi(\mathbf v)\) of a 1-form (the gradient of its potential, \(\mathbb P \Pi_L = \mathbb P_L\) without the retained \(k_\perp = 0\) polarization).

longitudinal_transpose(w, out=None)

The transpose \(\Pi_L^\top \mathbf w = \mathbb P \mathbb G (\ldots)^{-1} \mathbb G^\top \mathbf w\) on a dual vector (a current).

dot(v, out=None)

Apply the LinearOperator self to the Vector v.

The result is written to the Vector out, if provided.

Parameters:

• v (Vector) – The vector to which the linear operator (self) is applied. It must belong to the domain of self.

• out (Vector) – The vector in which the result of the operation is stored. It must belong to the codomain of self. If out is None, a new vector is created and returned.

Returns:

The result of the operation. If out is None, a new vector is returned. Otherwise, the result is stored in out and out is returned.

Return type:

Vector

class struphy.propagators.gyrokinetic_darwin.DarwinFieldCoupling(species: tuple[PICVariable, ...] = (), background: BackgroundWeights | None = None, polarization=0.0, background_current: bool = False, gyro: GyroAverage | None = None, with_drifts: bool = True)

Bases: Propagator, _PerturbedFieldMixin

The quasi-neutral gyrokinetic Darwin model (Sec. 6 of the paper) in the gauge-invariant form, coupled to the gyrocentres of all species: find \(\mathbf E \in H(\textnormal{curl})\), \(\mathbf B \in H(\textnormal{div})\) and \(f_s\) such that

\[\begin{split}\frac{\partial \mathbf D^{\textnormal{Dar}}}{\partial t} &= \nabla \times \left(\mathbf B - \frac{p_{0,\parallel} - p_{0,\perp}}{B_0^2} \mathbf B_\perp + \sum_s \alpha_s^2 \int f_s\, \mu\, \mathbf b_0 \right) - \sum_s \frac{\alpha_s^2}{\varepsilon_s} \int f_s\, \dot{\mathbf R}\,, \qquad \mathbf D^{\textnormal{Dar}} = \Pi_{L,\perp}\left(\sum_s \alpha_s^2 \frac{n_{0,s}}{B_0^2} \mathbf E_\perp\right)\,, \\ \frac{\partial \mathbf B}{\partial t} &= - \nabla \times \mathbf E\,,\end{split}\]

with \(\mathbf B = \mathbf B_0 + \tilde{\mathbf B}\) (background_current) and the gyro-averaged currents of the gyrocentre motion (parallel streaming along \(\mathbf B^*\), E x B and grad-B drifts) of GyroEparCoupling and GyroDriftCoupling. The quasi-neutrality \(\nabla \cdot \mathbf D^{\textnormal{Dar}} = \rho_f\) is the divergence of the first equation, and the transversal part of the first equation (no displacement current) is the constraint of the differential-algebraic system: the fields are the gauge-invariant \(\mathbf E = -\nabla \phi - \partial_t \mathbf A\) and \(\tilde{\mathbf B} = \nabla \times \mathbf A\) of Eqs. (6.15)-(6.16), without the potentials.

In matrix form, with the longitudinal polarization operator DarwinLongitudinalOperator \(\mathbb P_L\) in place of the dielectric matrix \(\mathbb A\) of MaxwellGyroMagnetisation:

\[\mathbb P_L \dot{\mathbf e} = \mathbb C^\top \mathbb M^2_H (\mathbf b + \mathbf b_0) + \mathbb C^\top (\mathcal P^b)^\top \sum_s \alpha_s^2 \mathbf m_s - \sum_s \frac{\alpha_s^2}{\varepsilon_s} \mathbf j_s(v_\parallel) - \sum_s \mathbb Q_s \mathbf e - \mathbf j_{\nabla B}\,, \qquad \dot{\mathbf b} = -\mathbb C \mathbf e\,, \qquad \dot v_{\parallel,p} = \frac{1}{\varepsilon_s} \mathbf c_p \cdot \langle \vec \Lambda^{1\top} \rangle_p \mathbf e\,,\]

discretized with the implicit mid-point rule in \((\mathbf e, \mathbf b, v_\parallel)\), all couplings included (\(\mathbf c_p\) frozen at \(v_\parallel^n\), the markers at \(\boldsymbol \eta^n\)):

\[\left(\mathbb P_L + \frac{\Delta t^2}{4} (\mathbb K + \mathbb S) + \frac{\Delta t}{2} \mathbb Q\right) \mathbf e^{n+1} = \left(\mathbb P_L - \frac{\Delta t^2}{4} (\mathbb K + \mathbb S) - \frac{\Delta t}{2} \mathbb Q\right) \mathbf e^{n} + \Delta t \left(\mathbb C^\top \mathbb M^2_H (\mathbf b^n + \mathbf b_0) + \mathbf s_m - \mathbf j^n_\parallel - \mathbf j_{\nabla B}\right)\,,\]

with \(\mathbb K = \mathbb C^\top \mathbb M^2_H \mathbb C\), \(\mathbb S = \sum_s \frac{\alpha_s^2}{\varepsilon_s^2} \mathbb S_s\), \(\mathbb Q = \sum_s \mathbb Q_s\) (stencil matrices in the ZLR limit, matrix-free with finite Larmor radius) and the magnetisation source \(\mathbf s_m\); then \(\mathbf b^{n+1} = \mathbf b^n - \Delta t\, \mathbb C \bar{\mathbf e}\), the parallel velocities are pushed with \(\bar{\mathbf e} = (\mathbf e^n + \mathbf e^{n+1})/2\), the drifts with an explicit Runge-Kutta scheme in \(\mathbf b^n\) and the longitudinal field \(\Pi_L \mathbf e^n = -\mathbb G \boldsymbol \phi^n\) before the field solve (the couplings are accumulated at the new positions). The E x B operator and the grad-B current act on the longitudinal field as well, \(\mathbb Q \to \Pi_L^\top \mathbb Q \Pi_L\) and \(\mathbf j_{\nabla B} \to \Pi_L^\top \mathbf j_{\nabla B}\): the transversal part of the Darwin field is the multiplier of the quasi-static transversal Ampère law, not a dynamical field, and does not drive the drifts (as in the \(A_\parallel\) gyrokinetic models, where the E x B drift is electrostatic). The operator is positive definite (longitudinal part from \(\mathbb P_L\), transversal part from \(\mathbb K + \mathbb S\)), the system is solved with pbicgstab (pcg without drifts) preconditioned by the weighted mass matrix of the polarization plus the continuum limit of the skin term. The mid-point rule conserves the Darwin energy \(\frac12 \mathbf e^\top \mathbb P_L \mathbf e + \frac12 \mathbf b^\top \mathbb M^2_H \mathbf b + \sum_s \alpha_s^2 \sum_p w_p (v_\parallel^2/2 + \mu B_\parallel)\) exactly for the parallel and E x B couplings (the grad-B current exchange is second order).

class Variables

Bases: object

Container for the variables advanced by DarwinFieldCoupling.

Variables:

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

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

class Options(solver: str | None = None, precond: Literal['MassMatrixPreconditioner', 'PolarMassMatrixPreconditioner', None] = 'MassMatrixPreconditioner', solver_params: SolverParameters | None = None, inner_tol: float = 1e-12, inner_maxiter: int = 1000, kernel_average: str = 'eta1 eta2', stab_eps: float = 1.0, perp_polarization: bool = True, picard_maxiter: int = 20, picard_tol: float = 1e-10, butcher: ButcherTableau | None = None, mpi_sort: Literal['each', 'last', None] = 'last')

Bases: OptionsBase

Configuration options for DarwinFieldCoupling.

Parameters:

• solver (str, default=None) – Iterative solver: a general one (LiteralOptions.OptsGenSolver) solves the full system including the anti-symmetric E x B operator, a symmetric one (LiteralOptions.OptsSymmSolver) solves the symmetric part with Picard iterations for the E x B operator (picard_maxiter, picard_tol; they converge only for small time steps). None: "pbicgstab" with drifts, "pcg" without.

• precond (LiteralOptions.OptsMassPrecond, default="MassMatrixPreconditioner") – Preconditioner, built from the weighted mass matrix of the polarization plus the continuum limit of the skin term (skin_preconditioner_matrix()).

• solver_params (SolverParameters, default=None) – Iterative-solver controls. If None, defaults to SolverParameters(). The tolerance is relative to the norm of the right-hand side (the system has a kernel, the uniform perpendicular field, and its scale is set by the skin term).

• inner_tol (float, default=1e-12) – Tolerance of the inner solves of the perpendicular Laplacian (relative to the right-hand side).

• inner_maxiter (int, default=1000) – Maximum iterations of the inner solves.

• kernel_average (str, default="eta1 eta2") – The directions of the average projection that regularizes the kernel of the perpendicular Laplacian (the directions perpendicular to the background field; "" for the global mean only).

• stab_eps (float, default=1.0) – The regularization \(\sigma\) of that kernel.

• perp_polarization (bool, default=True) – Whether the (transversal) polarization of the fields constant in the averaged directions (\(k_\perp = 0\), where the longitudinal projection is degenerate and the Darwin model has no inertia) is retained, see DarwinLongitudinalOperator.

• picard_maxiter (int, default=20) – Maximum number of Picard iterations for the E x B operator (symmetric solvers only).

• picard_tol (float, default=1e-10) – Relative tolerance of the Picard iterations.

• butcher (ButcherTableau, default=None) – Butcher tableau of the explicit drift push, default rk2.

• mpi_sort (LiteralOptions.OptsMPIsort, default="last") – MPI sorting policy of the drift push.

property longitudinal_operator: DarwinLongitudinalOperator

The longitudinal polarization operator (available after allocate()).

MHD-kinetic coupling

Particle and FEEC variables are updated.

class struphy.propagators.current_coupling_5d_curlb.CurrentCoupling5DCurlb(b_tilde: FEECVariable | None = None)

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(ep_scale: float = 1.0, u_space: Literal['Hcurl', 'Hdiv', 'H1vec'] = 'Hdiv', solver: Literal['pcg', 'cg'] = 'pcg', precond: Literal['MassMatrixPreconditioner', 'PolarMassMatrixPreconditioner', None] = 'MassMatrixPreconditioner', solver_params: SolverParameters | None = None, filter_params: FilterParameters | None = None)

Bases: OptionsBase

Configuration options for CurrentCoupling5DCurlb.

Parameters:

• 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. Only "Hdiv" and "H1vec" are supported.

• 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(energetic_ions: PICVariable | None = None, b_tilde: FEECVariable | None = None)

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(ep_scale: float = 1.0, u_space: Literal['Hcurl', 'Hdiv', 'H1vec'] = 'Hdiv', solver: Literal['pbicgstab', 'bicgstab', 'gmres'] = 'pbicgstab', precond: Literal['MassMatrixPreconditioner', 'PolarMassMatrixPreconditioner', None] = 'MassMatrixPreconditioner', solver_params: SolverParameters | None = None, filter_params: FilterParameters | None = None)

Bases: OptionsBase

Configuration options for CurrentCoupling5DDensity.

Parameters:

• 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.OptsGenSolver, default="pbicgstab") – General 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(b_tilde: FEECVariable | None = None)

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(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', 'PolarMassMatrixPreconditioner', 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:

• 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. "Hdiv" and "H1vec" are supported for algo="explicit", only "Hdiv" for algo="discrete_gradient".

• 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(b_tilde: FEECVariable | None = None)

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(u_space: Literal['Hcurl', 'Hdiv', 'H1vec'] = 'Hdiv', solver: Literal['pcg', 'cg'] = 'pcg', precond: Literal['MassMatrixPreconditioner', 'PolarMassMatrixPreconditioner', 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:

• 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(energetic_ions: PICVariable | None = None, b_tilde: FEECVariable | None = None)

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(u_space: Literal['Hcurl', 'Hdiv', 'H1vec'] = 'Hdiv', solver: Literal['pcg', 'cg'] = 'pcg', precond: Literal['MassMatrixPreconditioner', 'PolarMassMatrixPreconditioner', 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:

• 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(energetic_ions: PICVariable | None = None)

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(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', 'PolarMassMatrixPreconditioner', None] = 'MassMatrixPreconditioner', solver_params: SolverParameters | None = None, filter_params: FilterParameters | None = None, butcher: ButcherTableau | None = None, nonlinear: bool = True)

Bases: OptionsBase

Configuration options for ShearAlfvenCurrentCoupling5D.

Parameters:

• 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="MassMatrixPreconditioner") – 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', 'PolarMassMatrixPreconditioner', 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: LinearOperator

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.

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_v_in_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_v_in_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(phi: FEECVariable | None = None)

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(c_fun: Literal['const'] = 'const', coupling: float = 1.0, kappa: float = 1.0, nu: float = 0.01, butcher: ButcherTableau | None = None, solver: Literal['pcg', 'cg'] = 'pcg', precond: Literal['MassMatrixPreconditioner', 'PolarMassMatrixPreconditioner', None] = 'MassMatrixPreconditioner', solver_params: SolverParameters | None = None)

Bases: OptionsBase

Configuration options for HasegawaWakataniStep.

Parameters:

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

• coupling (float, default=1.0) – Constant value of \(C\) when c_fun="const".

• 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 any FEEC space ("H1", "Hcurl", "Hdiv", "L2" or "H1vec").

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.

update_source(t: float)

Set the source coefficients to \(c_0\, h(t)\) at time t.

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.