Available Propagators

Field solvers

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

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 adiabatic response of electrons.

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

\[\frac{1}{Z\epsilon^2} \int_\Omega \frac{n_0}{T_0} \psi\, (\phi - \langle \phi \rangle)\,\textrm d \mathbf x + \int_\Omega \frac{n_0}{|B_0|^2} \nabla \psi^\top \,\mathbb D \nabla \phi \,\textrm d \mathbf x = \frac{1}{\epsilon} \sum_i c_i \int_\Omega \psi\, \rho_i(\mathbf x)\,\textrm d \mathbf x \qquad \forall \ \psi \in H^1\,,\]

where \(c_i\) are the rho_coeffs, \(\mathbb D\) can be an anisotropic (gyrokinetic) diffusion operator (or the identity; the weight \(n_0/|B_0|^2\) is that of diffusion_mat="M1gyro", "M1" and "M1perp" are unweighted), \(\epsilon \in \mathbb R\) is the gyrokinetic ratio defined in units, Z the charge number of ions and \(\langle \phi \rangle\) the (weighted) average over \(\eta_3\) or over \(\eta_2, \eta_3\) (which_geometry, or zero for which_geometry=None): the electrons do not respond to the flux-surface averaged potential. For a full-f particle source (no control variate) the neutralizing background \(-n_0\) is added to \(\rho_i\) automatically (see ImplicitDiffusion). Boundary terms from integration by parts are assumed to vanish.

The equation is discretized as

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

where \(\mathbb M^0_{\textnormal{ad},T}\) is the mass matrix weighted with \(n_0/T_0\) (stab_mat, default "M0ad_withT"), \(\mathbb A\) the weighted average projection and \(\mathbb M^1_\textnormal{gyro}\) the gyrokinetic mass matrix (diffusion_mat).

Parameters:

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

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

• epsilon (float) – Gyrokinetic parameter.

• Z (int) – Charge number of the ions.

• diagnostic (FEECVariable, optional) – If given, the L2 projection of the right-hand side is stored in it.

class Options(stab_mat: Literal['M0', 'M0ad', 'Id', 'M0ad_withT'] = 'M0ad_withT', diffusion_mat: Literal['M1', 'M1perp', 'M1gyro'] = 'M1gyro', which_geometry: Literal['cylindrical', 'toroidal', None] = '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:

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

  Mass matrix of the adiabatic term, multiplied by 1 / (Z * epsilon**2).

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

  • "M0ad": adiabatic-electron weighted 0-form mass operator (weight \(n_0\)).

  • "Id": identity operator.

  • "M0ad_withT": adiabatic-electron weighted 0-form mass operator (weight \(n_0/T_0\)).

• which_geometry ({"cylindrical", "toroidal", None}, 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". With None there is no average, \(\langle \phi \rangle = 0\): the electrons respond adiabatically to all modes (e.g. for a comparison with the full-orbit model in a slab, see AdiabaticPoissonSolve). 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"}, default="M1gyro") – Diffusion matrix.

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

Directions averaged in the stabilization, see which_geometry.

class struphy.propagators.adiabatic_poisson_solve.AdiabaticPoissonSolve(rho: FEECVariable | ParticlesToGrid | Callable | list | None = None, rho_coeffs: float | list | None = None, epsilon: float = 1.0, alpha: float = 1.0, Z: int = 1, diagnostic: FEECVariable | None = None)

Bases: ImplicitDiffusion

Weak discretization of the Poisson equation for full-orbit ions with adiabatic electrons.

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

\[\frac{1}{Z\varepsilon^2} \int_\Omega \frac{n_0}{T_0}\, \psi\, (\phi - \langle \phi \rangle)\,\textrm d \mathbf x + \frac{1}{\alpha^2} \int_\Omega \nabla \psi \cdot \nabla \phi \,\textrm d \mathbf x = \frac{1}{\varepsilon} \int_\Omega \psi\, (n_i - n_0)\,\textrm d \mathbf x \qquad \forall \ \psi \in H^1\,,\]

where \(n_i = \int f \,\textrm d^3 \mathbf v\) is the ion density, \(n_0\) and \(T_0\) are the equilibrium ion density and electron temperature, \(Z\) is the ion charge number, \(\varepsilon = 1/(\hat \Omega_c \hat t)\) and \(\alpha = \hat \omega_p / \hat \Omega_c\). The first term is the linearized Boltzmann response of the electrons. The optional flux-surface average \(\langle \phi \rangle\) (which_geometry) is the zonal correction of GyrokineticPoissonSolve: electrons do not respond to the averaged potential. The second term is the vacuum (Debye) term; it is dropped in the quasi-neutral limit \(\alpha \to \infty\) (quasi_neutral=True), where the equation becomes a weighted \(L^2\) projection. Boundary terms from integration by parts are assumed to vanish.

With polarization=True (quasi-neutral only), the polarization density of guiding centers is added as in GyrokineticPoissonSolve,

\[\frac{1}{Z\varepsilon^2} \int_\Omega \frac{n_0}{T_0}\, \psi\, (\phi - \langle \phi \rangle)\,\textrm d \mathbf x + \int_\Omega \frac{n_0}{|\mathbf B_0|^2} \nabla_\perp \psi \cdot \nabla_\perp \phi \,\textrm d \mathbf x = \frac{1}{\varepsilon} \int_\Omega \psi\, (n_i - n_0)\,\textrm d \mathbf x\,,\]

discretized with the gyrokinetic mass matrix \(\mathbb M^1_\textnormal{gyro}\) in place of \(\alpha^{-2} \mathbb M^1\). This is needed when \(n_i\) is the density of the marker positions and these are guiding-center positions whose density does not contain the polarization density (e.g. classical Pauli particles with the control variate, see VlasovPoissonAdiabaticPauli).

In the quasi-neutral limit the ion-acoustic branch \(\omega = k c_s\) has no Debye cutoff, so an explicit particle push needs \(\Delta t\, c_s k_\textnormal{max} \lesssim 2\), with \(k_\textnormal{max} \approx \pi / h\) for the smallest cell size \(h\); refining the grid requires a smaller time step.

The equation is discretized as (multiplied by \(\alpha^2\) if quasi_neutral=False)

\[\left( \frac{1}{Z\varepsilon^2}\,\mathbb M^0_{\textnormal{ad},T} (\mathbb I - \mathbb A) + \frac{1}{\alpha^2} \mathbb G^\top \mathbb M^1 \mathbb G \right)\, \boldsymbol\phi = \frac{1}{\varepsilon} (\Lambda^0, n_i - n_0 )_{L^2}\,,\]

where \(\mathbb M^0_{\textnormal{ad},T}\) is the mass matrix weighted with \(n_0/T_0\) (stab_mat="M0ad_withT") and \(\mathbb A\) is the weighted average projection (\(\mathbb A = 0\) if which_geometry=None). The particle density is deposited with the marker weights; the background \(n_0\) is subtracted either by the control variate or, for full-f markers, by an explicit background source (see ImplicitDiffusion).

class Options(quasi_neutral: bool = True, polarization: bool = False, which_geometry: Literal['cylindrical', 'toroidal', None] | None = None, stab_mat: Literal['M0', 'M0ad', 'M0ad_withT', 'Id'] = 'M0ad_withT', 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: OptionsBase

Configuration options for AdiabaticPoissonSolve.

Parameters:

• quasi_neutral (bool, default=True) – If True, drop the vacuum term \(\alpha^{-2}\Delta \phi\) (limit \(\alpha \to \infty\)); the system is then a weighted mass-matrix solve (unless polarization=True).

• polarization (bool, default=False) – If True, add the polarization term \(-\nabla \cdot (n_0 |\mathbf B_0|^{-2} \nabla_\perp \phi)\) (diffusion matrix "M1gyro"). Requires quasi_neutral=True. Only for markers at guiding-center positions whose deposited density lacks the polarization density (classical Pauli particles with the control variate); for full-orbit ions or full-f Pauli particles the polarization density is counted twice.

• which_geometry ({None, "cylindrical", "toroidal"}, default=None) – Zonal correction \(\phi - \langle \phi \rangle\) in the adiabatic electron response: no average (None), average over eta3 ("cylindrical") or over eta2 and eta3 ("toroidal"), as in GyrokineticPoissonSolve. Requires quasi_neutral=False or polarization=True: without a diffusion term the averaged potential is undetermined.

• stab_mat ({"M0", "M0ad", "M0ad_withT", "Id"}, default="M0ad_withT") – Mass matrix of the adiabatic term ("M0ad_withT" has the weight \(n_0/T_0\)).

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

• solver (LiteralOptions.OptsSymmSolver, default="pcg") – Name of the symmetric iterative solver.

• precond (LiteralOptions.OptsDiffusionPrecond, default=None) – Name of the preconditioner, see ImplicitDiffusion. "MultiGrid" and "StiffnessPreconditioner" need a diffusion term (quasi_neutral=False or polarization=True).

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

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

• solver_params (SolverParameters, default=None) – Iterative-solver controls. 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

AdiabaticPoissonSolve reuses ImplicitDiffusion with sigma_2 = 0, divide_by_dt = False, sigma_1 = 1 / (Z * epsilon**2), sigma_3 = 1 / epsilon and diffusion_mat=None (quasi-neutral), or sigma_1 = alpha**2 / (Z * epsilon**2), sigma_3 = alpha**2 / epsilon and diffusion_mat="M1", or, with polarization=True, sigma_1 = 1 / (Z * epsilon**2), sigma_3 = 1 / epsilon and diffusion_mat="M1gyro".

property stab_average: str | None

Directions averaged in the adiabatic term, see which_geometry.

property diffusion_mat: str | None

none in the quasi-neutral limit, "M1gyro" with the polarization term.

Type:

Diffusion matrix of ImplicitDiffusion

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', 'MultiGrid', 'StiffnessPreconditioner', None] = 'StiffnessPreconditioner', exb_correction: bool = True, multigrid: MultiGridOptions | None = None, 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.OptsShearAlfvenPrecond, default="StiffnessPreconditioner") – Preconditioner of the linear solve: of the mass-matrix block ("MassMatrixPreconditioner", "PolarMassMatrixPreconditioner"), or of the whole Schur operator \(\mathbb M^{2,n} + \frac{\Delta t^2}{4} \mathbb T^\top \mathbb C^\top \mathbb M^2 \mathbb C \mathbb T\): a geometric multigrid V-cycle ("MultiGrid", MultiGridPreconditioner) or the Kronecker (fast diagonalization) preconditioner ("StiffnessPreconditioner", the default, ShearAlfvenSchurPreconditioner: block Jacobi in the velocity components with the auxiliary space correction of the \(\mathbf E \times \mathbf B\) flows if exb_correction, see there). Both need the implicit algorithm with solver="pcg". The Kronecker preconditioner also needs u_space="Hdiv", tensor-product splines (no polar splines) and the global commuting projectors; otherwise "StiffnessPreconditioner" falls back (with a warning) to "MassMatrixPreconditioner" ("PolarMassMatrixPreconditioner" with polar splines). The Schur operator depends on the time step, so these preconditioners are built in the setup (prepare_time_step()) and updated if the time step changes. The multigrid preconditioner is experimental: on the TAE example (HollowTorus, 8x32x4) it needs more PCG iterations than the mass matrix preconditioner and is much slower (the smoother does not treat the large kernel of the curl term on 2-forms). The Kronecker preconditioner (with exb_correction) needs about 11 instead of 31 PCG iterations per step there (dt = 0.1), with larger gains for larger time steps.

• exb_correction (bool, default=True) – For precond="StiffnessPreconditioner": add the auxiliary space correction \(\mathbb N P_\phi \mathbb N^\top\) (AuxiliarySpacePreconditioner) of the \(\mathbf E \times \mathbf B\) flows \(\mathbf B_0 \times \nabla \phi\), \(\mathbb N = \hat \Pi^2[\mathbf B_0 \times \,\cdot\,] \, \mathbb G\) with \(\phi \in V^0_h\). The block Jacobi preconditioner in the velocity components does not reduce the incompressible shear Alfven modes with small parallel wave number (their velocity components cancel in the compression), which this correction resolves; \(P_\phi\) is the Kronecker (fast diagonalization) preconditioner of \(\mathbb N^\top S \mathbb N\) on \(V^0_h\).

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

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

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

prepare_time_step(dt: float)

Preconditioner of the Schur operator for the time step dt (multigrid or Kronecker): built in the setup, updated or rebuilt when dt changes.

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.shear_alfven_curl_curl_propagator.ShearAlfvenCurlCurlPropagator

Bases: Propagator

FEEC discretization of the following equations: find \(\tilde{\mathbf U} \in H(\textnormal{div})\) and \(\tilde{\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 \left[(\nabla \times \tilde{\mathbf B}) \times \mathbf B_0\right] \cdot \mathbf V \,\textnormal d \mathbf x \qquad \forall \,\mathbf V \in H(\textnormal{div})\,, \\[2mm] &\int_\Omega \frac{\partial \tilde{\mathbf{B}}}{\partial t} \cdot \mathbf C\,\textnormal d \mathbf x = - \int_\Omega \left[(\nabla \times \mathbf C) \times \mathbf B_0\right] \cdot \tilde{\mathbf{U}} \,\textrm d \mathbf x \qquad \forall \, \mathbf C \in H(\textrm{curl})\,.\end{split}\]

The momentum equation holds strongly up to the \(\rho_0\)-weighted \(L^2\) projection onto \(H(\textnormal{div})\), the induction equation is weak. In discrete form,

\[\mathbb M^\rho_2 \dot{\mathbf u} = - \mathbb M^{2,B} \mathbb C \mathbf b\,, \qquad \mathbb M_1 \dot{\mathbf b} = - \mathbb C^\top \mathbb M^{2,B} \mathbf u\,,\]

where \(\mathbb M^{2,B}\) (antisymmetric) is M2B. This is the discrete system of ShearAlfvenB1.

Time discretization:

• implicit: Crank-Nicolson (implicit mid-point) without the projection in the curl-curl term: the update of the velocity, \(\tilde{\mathbf U}^{n+1} = \tilde{\mathbf U}^n + \frac{\Delta t}{2 \rho_0} (\nabla \times (\tilde{\mathbf B}^{n+1} + \tilde{\mathbf B}^n)) \times \mathbf B_0\), is inserted into the weak induction equation before the projection, which gives the curl-curl system

  \[\left(\mathbb M_1 + \frac{\Delta t^2}{4} \mathbb C^\top \mathbb M^{2,\textnormal{gyro}} \mathbb C\right) \mathbf b^{n+1} = \left(\mathbb M_1 - \frac{\Delta t^2}{4} \mathbb C^\top \mathbb M^{2,\textnormal{gyro}} \mathbb C\right) \mathbf b^n - \Delta t \, \mathbb C^\top \mathbb M^{2,B} \mathbf u^n\,,\]

  with M2gyro, the mass matrix of \(\int \frac{1}{\rho_0} \left[(\nabla \times \mathbf C) \times \mathbf B_0\right] \cdot \left[(\nabla \times \tilde{\mathbf B}) \times \mathbf B_0\right] \textnormal d \mathbf x\). The velocity is then updated with the projection, \(\mathbf u^{n+1} = \mathbf u^n - \frac{\Delta t}{2} (\mathbb M^\rho_2)^{-1} \mathbb M^{2,B} \mathbb C (\mathbf b^{n+1} + \mathbf b^n)\). The perturbed energy \(\frac 12 \mathbf u^\top \mathbb M^\rho_2 \mathbf u + \frac 12 \mathbf b^\top \mathbb M_1 \mathbf b\) is therefore not exactly conserved (the projection error of the Lorentz force enters); see ShearAlfvenCurlCurlwithProjPropagator for the exactly energy conserving version. The system is solved by SchurSolver with \(x = \mathbf b\), \(y = \mathbf u\). With precond="HiptmairXu", it is preconditioned by HiptmairXuPreconditioner for the curl-curl operator (scaled by \(\sigma = 4/\Delta t^2\); TAE example 8x32x4: 37 instead of 706 PCG iterations per step for \(\Delta t = 2\)), optionally (aux_correction) plus the field-aligned auxiliary space corrections of AuxiliarySpacePreconditioner: perpendicular perturbations \(\mathbf B_0 \times \nabla \psi\) and parallel perturbations \(\phi \mathbf B_0\) with \(\psi, \phi \in V^0_h\), whose operators are approximated by the Kronecker StiffnessPreconditioner.

• explicit: explicit RK methods from ButcherTableau (for the semi-discrete system above).

class Variables

Bases: object

Container for variables advanced by ShearAlfvenCurlCurlPropagator.

Variables:

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

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

class Options(algo: Literal['implicit', 'explicit'] = 'implicit', solver: Literal['pcg', 'cg'] = 'pcg', precond: Literal['HiptmairXu', 'MassMatrixPreconditioner', 'PolarMassMatrixPreconditioner', None] = 'HiptmairXu', aux_correction: bool = False, hx_aux_weight: float = 1.0, solver_params: SolverParameters | None = None, solver_M2n: Literal['pcg', 'cg'] = 'pcg', precond_M2n: Literal['MassMatrixPreconditioner', 'PolarMassMatrixPreconditioner', None] = 'MassMatrixPreconditioner', solver_params_M2n: SolverParameters | None = None, butcher: ButcherTableau | None = None)

Bases: OptionsBase

Configuration options for ShearAlfvenCurlCurlPropagator.

Parameters:

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

• solver (LiteralOptions.OptsSymmSolver, default="pcg") – Solver of the Schur system in b (implicit) or of M1 (explicit).

• precond (LiteralOptions.OptsShearAlfvenCurlCurlPrecond, default="HiptmairXu") – Preconditioner of the Schur system: Hiptmair-Xu for the curl-curl operator with M2gyro in the middle ("HiptmairXu", implicit algorithm with solver="pcg" only, no polar splines; it depends on the time step, so it is built in the setup (prepare_time_step()) and rebuilt when the time step changes), or a preconditioner of the mass matrix M1 ("MassMatrixPreconditioner", "PolarMassMatrixPreconditioner"). With algo="explicit", "HiptmairXu" is replaced by "MassMatrixPreconditioner" (for M1).

• aux_correction (bool, default=False) – With precond="HiptmairXu": add the field-aligned auxiliary space corrections of the perturbations \(\mathbf B_0 \times \nabla \psi\) and \(\phi \mathbf B_0\) (\(\psi, \phi \in V^0_h\)), which the Hiptmair-Xu preconditioner of the anisotropic curl-curl operator reduces only slowly. They lower the condition number (TAE geometry, 4x12x3, dt = 1: 34 to 19) but hardly the PCG iterations (tilted slab 16x16x8, dt = 2: 148 to 125; TAE example 8x32x4: unchanged), at about twice the cost per iteration.

• hx_aux_weight (float, default=1.0) – Factor of the vector auxiliary space correction of the Hiptmair-Xu preconditioner (aux_weight of HiptmairXuPreconditioner, 0 switches it off; M2gyro vanishes along the background field, so the isotropic vector Laplacian of the correction does not match the curl-curl term).

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

• solver_M2n (LiteralOptions.OptsSymmSolver, default="pcg") – Solver of the inner M2n inversion.

• precond_M2n (LiteralOptions.OptsMassPrecond, default="MassMatrixPreconditioner") – Preconditioner of the inner M2n inversion.

• solver_params_M2n (SolverParameters, default=None) – Solver controls of the inner M2n inversion; defaults to SolverParameters(tol=1e-14) (the outer Krylov solver needs an accurate inner inverse).

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

prepare_time_step(dt: float)

Hiptmair-Xu preconditioner of S ~ dt^2/4 (C^T M2gyro C + 4/dt^2 M1) for the time step dt: built in the setup, rebuilt when dt changes (sigma = 4/dt^2 with the unscaled M2gyro in the middle, as the vector auxiliary space correction of Hiptmair-Xu assumes a unit-scale curl-curl term).

class struphy.propagators.shear_alfven_curl_curl_propagator.ShearAlfvenCurlCurlwithProjPropagator

Bases: ShearAlfvenCurlCurlPropagator

FEEC discretization of the equations of ShearAlfvenCurlCurlPropagator (same semi-discrete system, \(\tilde{\mathbf U} \in H(\textnormal{div})\), \(\tilde{\mathbf B} \in H(\textnormal{curl})\)), with the projection in the curl-curl term, which conserves the perturbed energy \(\frac 12 \mathbf u^\top \mathbb M^\rho_2 \mathbf u + \frac 12 \mathbf b^\top \mathbb M_1 \mathbf b\) exactly (up to the solver tolerances).

Time discretization:

• implicit: Crank-Nicolson (implicit mid-point). System size reduction via SchurSolver with \(x = \mathbf b\), \(y = \mathbf u\); the Schur complement is

  \[S = \mathbb M_1 + \frac{\Delta t^2}{4} \mathbb C^\top (\mathbb M^{2,B})^\top (\mathbb M^\rho_2)^{-1} \mathbb M^{2,B} \mathbb C\,,\]

  where \((\mathbb M^\rho_2)^{-1}\) is an inner iterative solve. The curl-curl operator with M2gyro (the continuous counterpart) is not spectrally equivalent to \(S\) (it overestimates \(S\) on the fields whose Lorentz force is removed by the projection), so the Hiptmair-Xu preconditioner (precond="HiptmairXu") does not pay off here; the default is the mass matrix preconditioner of \(\mathbb M_1\).

• explicit: explicit RK methods from ButcherTableau.

class Options(algo: Literal['implicit', 'explicit'] = 'implicit', solver: Literal['pcg', 'cg'] = 'pcg', precond: Literal['HiptmairXu', 'MassMatrixPreconditioner', 'PolarMassMatrixPreconditioner', None] = 'MassMatrixPreconditioner', aux_correction: bool = False, hx_aux_weight: float = 1.0, solver_params: SolverParameters | None = None, solver_M2n: Literal['pcg', 'cg'] = 'pcg', precond_M2n: Literal['MassMatrixPreconditioner', 'PolarMassMatrixPreconditioner', None] = 'MassMatrixPreconditioner', solver_params_M2n: SolverParameters | None = None, butcher: ButcherTableau | None = None)

Bases: Options

Configuration options for ShearAlfvenCurlCurlwithProjPropagator: as for ShearAlfvenCurlCurlPropagator, with the default precond="MassMatrixPreconditioner".

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}\]

The magnetic field perturbation \(\tilde{\mathbf B}\) is not updated (it enters only the Lorentz force) and is either a 2-form (\(H(\textnormal{div})\), default, with M2J) or a 1-form (\(H(\textnormal{curl})\), with the mixed mass matrix M21J; only for \(\tilde{\mathbf U} \in H(\textnormal{div})\)).

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, divdiv_alpha: float = 0.0)

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\} .\]

With a div-div regularization (divdiv_alpha \(=\alpha > 0\)), the kinetic energy is \(\int \rho |\mathbf u|^2/2 + \alpha \rho (\nabla \cdot \mathbf u)^2 \,\textrm d \mathbf x\): the mass matrix \(\mathbb M^v[\hat\rho]\) is replaced by the metric \(\mathbb A[\hat\rho]\) of kinetic_metric() and \(\alpha \, d^{n+1} d^n / g\) is added to the argument of \(\hat l^3\), with \(d^k = (\mathbb D \mathbb U^v \mathbf u^k)^\top \vec{\boldsymbol \Lambda}^3\). Since \(\mathbb A\) is linear in \(\hat\rho\), the scheme conserves the regularized energy.

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', 'perturbative', 'linear_q', 'perturbative_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", "perturbative", "linear_q", "perturbative_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. The solves are done for normalized right-hand sides, so tol is relative.

• nonlin_solver (NonlinearSolverParameters, default=None) – Nonlinear iteration controls. The iteration stops when the norm of the residual is below tol times the norm of the state at \(t^n\) if that is smaller than one (absolute otherwise), or when it stagnates (with a warning).

class struphy.propagators.variational_entropy_evolve.VariationalEntropyEvolve(rho: FEECVariable | None = None, divdiv_alpha: float = 0.0)

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 \.\]

With a div-div regularization (divdiv_alpha \(> 0\)), the mass matrix \(\mathbb M^v[\hat\rho]\) is replaced by the metric \(\mathbb A[\hat\rho]\) of the regularized kinetic energy, see kinetic_metric().

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. The solves are done for normalized right-hand sides, so tol is relative.

• nonlin_solver (NonlinearSolverParameters, default=None) – Nonlinear iteration controls. The iteration stops when the norm of the residual is below tol times the norm of the state at \(t^n\) if that is smaller than one (absolute otherwise), or when it stagnates (with a warning).

class struphy.propagators.variational_mag_field_evolve.VariationalMagFieldEvolve(divdiv_alpha: float = 0.0)

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 \,.\]

With a div-div regularization (divdiv_alpha \(> 0\)), the mass matrix \(\mathbb M^v[\hat\rho]\) is replaced by the metric \(\mathbb A[\hat\rho]\) of the regularized kinetic energy, see kinetic_metric().

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. The solves are done for normalized right-hand sides, so tol is relative.

• nonlin_solver (NonlinearSolverParameters, default=None) – Nonlinear iteration controls. The iteration stops when the norm of the residual is below tol times the norm of the state at \(t^n\) if that is smaller than one (absolute otherwise), or when it stagnates (with a warning).

class struphy.propagators.variational_momentum_advection.VariationalMomentumAdvection(divdiv_alpha: float = 0.0)

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 \,.\]

With a div-div regularization (divdiv_alpha \(> 0\)), the mass matrix \(\mathbb M^v[\hat\rho]\) is replaced by the metric \(\mathbb A[\hat\rho]\) of the regularized kinetic energy, see kinetic_metric().

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}} \cdot \vec{\boldsymbol \Lambda}^v])^\top \mathbb M^2 \mathbf B^{n+\frac{1}{2}} - (\mathbb D \hat{\Pi}^{2}[\hat{p_h^{n}} \cdot \vec{\boldsymbol \Lambda}^v])^\top \hat{l}^3(\frac{g}{\gamma-1}) = 0 ~ , \\[2mm] &\frac{\mathbf b^{n+1}- \mathbf b^n}{\Delta t} + \mathbb C \hat{\Pi}^{1}[\hat{\mathbf B_h^{n}} \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}} \cdot \vec{\boldsymbol \Lambda}^v]] + (\gamma - 1)\hat{\Pi}^{3}[\hat{p_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} = (\mathbf{b}^{n})^\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} = (\boldsymbol p^{n})^\top \vec{\boldsymbol \Lambda}^3 \in V_h^3 \,.\]

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

The weights of the transport operators are taken at \(t^n\) (in all equations, which conserves the energy); the step is therefore only first-order accurate in time.

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', 'perturbative'] = '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", "perturbative"}, 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. The Newton systems are solved for a normalized right-hand side, so tol is relative to the norm of the residual.

• nonlin_solver (NonlinearSolverParameters, default=None) – Nonlinear iteration controls. The iteration stops when the norm of the residual is below tol times the norm of the state at \(t^n\) if that is smaller than one (absolute otherwise), or when it stagnates (with a warning).

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\), \(q \in L^2\) and \(\mathbf B \in H(\textnormal{div})\) such that, with the transport operators

\[\mathcal T_{\bar{\mathbf B}}(\mathbf v) = \nabla \times (\bar{\mathbf B} \times \mathbf v)\,, \qquad \mathcal S_{\bar q}(\mathbf v) = \nabla \cdot (\bar q \mathbf v) + \left(\frac \gamma 2 - 1\right) \bar q\, \nabla \cdot \mathbf v\,,\]

the following holds, depending on the option model:

• "full_q" (full MHD in \(q = \sqrt p\)):

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

• "linear_q" (\(\mathbf B\) and \(q\) are the perturbations of the equilibrium \(\mathbf B_0\), \(q_0\)):

  \[\begin{split}&\int_\Omega \partial_t (\rho_0 \mathbf u) \cdot \mathbf v\,\textrm d \mathbf x - \int_\Omega \left(\mathbf B \cdot \mathcal T_{\mathbf B_0}(\mathbf v) + \mathbf B_0 \cdot \mathcal T_{\mathbf B}(\mathbf v)\right)\textrm d \mathbf x - \int_\Omega \frac{2}{\gamma -1} \left(q\, \mathcal S_{q_0}(\mathbf v) + q_0\, \mathcal S_q(\mathbf v)\right)\textrm d \mathbf x = 0\,, \\[2mm] &\partial_t \mathbf B + \mathcal T_{\mathbf B_0}(\mathbf u) = 0 \,, \qquad \partial_t q + \mathcal S_{q_0}(\mathbf u) = 0 \,.\end{split}\]

• "perturbative_q" (\(\mathbf B\) and \(q\) are the perturbations, all nonlinear terms are kept):

  \[\begin{split}&\int_\Omega \partial_t (\rho \mathbf u) \cdot \mathbf v\,\textrm d \mathbf x - \int_\Omega \left(\mathbf B \cdot \mathcal T_{\mathbf B_0 + \mathbf B}(\mathbf v) + \mathbf B_0 \cdot \mathcal T_{\mathbf B}(\mathbf v)\right)\textrm d \mathbf x - \int_\Omega \frac{2}{\gamma -1} \left(q\, \mathcal S_{q_0 + q}(\mathbf v) + q_0\, \mathcal S_q(\mathbf v)\right)\textrm d \mathbf x = 0\,, \\[2mm] &\partial_t \mathbf B + \mathcal T_{\mathbf B_0 + \mathbf B}(\mathbf u) = 0 \,, \qquad \partial_t q + \mathcal S_{q_0 + q}(\mathbf u) = 0 \,.\end{split}\]

The second-order terms come from the auxiliary variables \(\hat{\mathbf B}\) (bt2) and \(\hat q\) (qt3), which, if given, are advanced with \(\partial_t \hat{\mathbf B} + \mathcal T_{\mathbf B}(\mathbf u) = 0\) and \(\partial_t \hat q + \mathcal S_q(\mathbf u) = 0\). They enter the conserved energy of the linear and perturbative models as \(\int \mathbf B_0 \cdot \hat{\mathbf B} + 2 q_0 \hat q / (\gamma - 1)\). For "full_q" the mass density \(\rho\) is the variable of VariationalDensityEvolve, for "perturbative_q" it is \(\rho_0\) plus that variable.

The scheme below is written for "full_q"; for the other models the operators and linear forms are replaced by the sums above. Below, hats denote pullbacks to the logical domain (not the auxiliary variables). The weights of the transport operators are taken at the same time in the momentum, induction and pressure equations, which makes the step energy conserving: at \(t^n\) as written below (first order in time), or at the midpoint \(t^{n+1/2}\) with the option implicit_weights (second order in time, the default for "linear_q" and "perturbative_q").

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{2 \hat{q}}{\gamma -1} \left(\nabla \cdot (\hat{q} \hat{\mathbf{v}}) + (\gamma/2 - 1) \hat{q} \nabla \cdot \hat{\mathbf{v}}\right) \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{q} \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}} \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} = (\mathbf{b}^{n})^\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} = (\boldsymbol q^{n})^\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', 'perturbative_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, implicit_weights: bool | None = None)

Bases: OptionsBase

Configuration options for VariationalQBEvolve.

Parameters:

• model ({"full_q", "linear_q", "perturbative_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. The Newton systems are solved for a normalized right-hand side, so tol is relative to the norm of the residual.

• nonlin_solver (NonlinearSolverParameters, default=None) – Nonlinear/Picard iteration controls. The iteration stops when the norm of the residual is below tol times the norm of the state at \(t^n\) if that is smaller than one (absolute otherwise), or when it stagnates (with a warning).

• implicit_weights (bool, default=None) – If True, the weights of the transport operators (\(\mathbf B\) and \(q\)) are taken at the midpoint \(t^{n+1/2}\) and updated in every Picard iteration, which makes the step second-order accurate in time. If False, they are taken at \(t^n\); the step is then only first-order accurate (e.g. for sound waves). Both conserve the energy. None means True for "linear_q" and "perturbative_q", and False for "full_q". With True, the Jacobian of the Picard iteration does not contain the dependence of the operators on their weights, so a step needs several iterations (typically 2 to 10, more for large time steps or amplitudes) instead of one, each with a Schur solve; False is the cheaper first-order choice.

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', 'perturbative_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", "perturbative_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. The solves are done for normalized right-hand sides, so tol is relative.

• nonlin_solver (NonlinearSolverParameters, default=None) – Nonlinear iteration controls. The iteration stops when the norm of the residual is below tol times the norm of the state at \(t^n\) if that is smaller than one (absolute otherwise), or when it stagnates (with a warning).

• 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, divdiv_alpha: float = 0.0)

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 ~.\]

With a div-div regularization (divdiv_alpha \(> 0\)), the mass matrix \(\mathbb M^v[\hat\rho]\) is replaced by the metric \(\mathbb A[\hat\rho]\) of the regularized kinetic energy, see kinetic_metric().

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', 'perturbative_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", "perturbative_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. The solves are done for normalized right-hand sides, so tol is relative.

• nonlin_solver (NonlinearSolverParameters, default=None) – Nonlinear iteration controls. The iteration stops when the norm of the residual is below tol times the norm of the state at \(t^n\) if that is smaller than one (absolute otherwise), or when it stagnates (with a warning).

• 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, delta_f_source: bool = False)

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 "H1" accumulator gets the coefficient of its accumulator), or one entry per collected source (no background source is collected with delta_f_source=True, the control variate or delta-f markers). 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.

• delta_f_source (bool, default=False) – Deposit particle sources with the delta-f weights \((f - f_0)/s\) also for full-f markers (without the control variate): the weights are updated before and the full-f weights restored after the deposit, and no background source is added. The Monte-Carlo noise is then that of \(f - f_0\) instead of \(f\), which is much smaller for small perturbations; meant for the initial Poisson solve of full-f models.

Note

The deposited particle density is \(f - f_0\) with the control variate (the weights are updated before each deposit), with delta-f markers (DeltaFParticles6D) and with delta_f_source=True. Otherwise (full-f), for "H1" accumulators the neutralizing background \(-n_0\) is added as a separate source with the coefficient of the particle source.

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', 'M0ad_withT', 'Id'] = 'M0', stab_average: str | None = None, diffusion_mat: Literal['M1', 'M1perp', None] = '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", "M0ad_withT", "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 (weight \(n_0\)).

  • "M0ad_withT": adiabatic-electron weighted 0-form mass operator (weight \(n_0/T_0\)).

  • "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", None} 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. With None there is no diffusion term and the system is sigma_1 * stab_mat only (a weighted L2 projection); then precond cannot be "MultiGrid" or "StiffnessPreconditioner".

• 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, magnetic_moment: bool = False)

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

With magnetic_moment=True, each marker carries a magnetic moment \(\mu_p\) (classical Pauli particle, see PauliParticles6D) and the force is augmented by the magnetic-moment force of the equilibrium field,

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

i.e. \(\mathbf F = -\nabla(\phi + \varepsilon \mu_p |\mathbf B_0|)\), where \(\nabla |\mathbf B_0|\) is the projected 1-form gradB1 of the equilibrium (the same as in the guiding-center pushers).

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, e_field: 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. Instead of the potential, an electric field \(\mathbf E \in H(\textnormal{curl})\) can be passed (e_field, explicit algorithm); then \(\mathbf E^* = \mathbf E - \varepsilon \mu_p \nabla |\mathbf B|\) (gauge-invariant electromagnetic formulation). 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, e_field: 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. Instead of the potential, an electric field \(\mathbf E \in H(\textnormal{curl})\) can be passed (e_field, explicit algorithm); then \(\mathbf E^* = \mathbf E - \varepsilon \mu_p \nabla |\mathbf B|\) (gauge-invariant electromagnetic formulation). 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, matrix_free: bool = False)

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

• matrix_free (bool, default=False) – Apply the particle matrix matrix-free, with its rank-one filling stored per marker (MarkerOperator), instead of assembling it: one spline evaluation and one deposition per marker and solver iteration instead of the products of all pairs of basis functions per marker. Needs 4 auxiliary marker columns. Ignored on the CuPy backend.

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}\]

With \(f\) approximated by \(N\) markers with positions \(\boldsymbol \eta_p\) (logical, fixed in this step), velocities \(\mathbf v_p\) and weights \(w_p\) (including the factor \(1/N\)), and the FE coefficients \(\mathbf e \in \mathbb R^{N_1}\) of \(\mathbf E \in V^1_h\), the semi-discrete system reads

\[\mathbb M_1 \dot{\mathbf e} = - \frac{\alpha^2}{\varepsilon} \, (\mathbf \Lambda^1)^\top \, \overline{DF^{-1}} \, \mathbb W \, \mathbf V \,, \qquad \dot{\mathbf V} = \frac{1}{\varepsilon} \, \overline{DF^{-\top}} \, \mathbf \Lambda^1 \, \mathbf e\,,\]

where \(\mathbf V = (\mathbf v_p)_p \in \mathbb R^{3N}\), \(\mathbb W = \textnormal{diag}(w_p)\) (acting on each velocity component), \(\mathbf \Lambda^1 \in \mathbb R^{3N \times N_1}\) holds the \(V^1_h\) basis functions evaluated at the markers, \((\mathbf \Lambda^1 \mathbf e)_p = \hat{\mathbf E}^1(\boldsymbol \eta_p)\), and \(\overline{DF^{-1}}\) is block-diagonal with the inverse Jacobian matrices \(DF^{-1}(\boldsymbol \eta_p)\). It conserves the energy \(\frac 12 \mathbf e^\top \mathbb M_1 \mathbf e + \frac{\alpha^2}{2} \mathbf V^\top \mathbb W \mathbf V\).

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

\[\mathbb M_1 \frac{\mathbf e^{n+1} - \mathbf e^n}{\Delta t} = - \frac{\alpha^2}{\varepsilon} \, (\mathbf \Lambda^1)^\top \, \overline{DF^{-1}} \, \mathbb W \, \frac{\mathbf V^{n+1} + \mathbf V^n}{2}\,, \qquad \frac{\mathbf V^{n+1} - \mathbf V^n}{\Delta t} = \frac{1}{\varepsilon} \, \overline{DF^{-\top}} \, \mathbf \Lambda^1 \, \frac{\mathbf e^{n+1} + \mathbf e^n}{2}\,,\]

which conserves the energy exactly (up to the solver tolerance; for constant weights, i.e. without control variate). The velocities are eliminated by SchurSolver with \(x = \mathbf e\), \(y = \mathbf V\), \(A = \mathbb M_1\), \(B = \frac{\alpha^2}{2 \varepsilon} (\mathbf \Lambda^1)^\top \overline{DF^{-1}} \mathbb W\) and \(C = - \frac{1}{2 \varepsilon} \overline{DF^{-\top}} \mathbf \Lambda^1\), which gives

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

with the accumulated matrix and vector (Accumulator, kernel vlasov_maxwell())

\[\mathbb A = (\mathbf \Lambda^1)^\top \, \overline{DF^{-1}} \, \mathbb W \, \overline{DF^{-\top}} \, \mathbf \Lambda^1 = \sum_p w_p \, \mathbf \Lambda^1_p{}^\top G^{-1}(\boldsymbol \eta_p) \, \mathbf \Lambda^1_p\,, \qquad \mathbf j = (\mathbf \Lambda^1)^\top \, \overline{DF^{-1}} \, \mathbb W \, \mathbf V^n\,.\]

\(\mathbb A\) is a particle approximation of the \(n\)-weighted mass matrix \(\mathbb M_1\), so the Schur operator is a mass matrix with weight \(1 + \Delta t^2 \alpha^2 n / (4 \varepsilon^2)\), solved with solver and precond (built from \(\mathbb M_1\)). The velocities are then updated by push_v_with_efield() with the mid-point field \(\frac 12 (\mathbf e^{n+1} + \mathbf e^n)\).

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, matrix_free: bool = False)

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.

• matrix_free (bool, default=False) – Apply the particle matrix matrix-free, with the filling \(w_p G^{-1}(oldsymbol \eta_p)\) stored per marker (MarkerOperator), instead of assembling it: one spline evaluation and one deposition per marker and solver iteration instead of the products of all pairs of basis functions per marker. Needs 6 auxiliary marker columns. Ignored on the CuPy backend.

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}\]

Ampère’s law holds weakly, Faraday’s law strongly. With the FE coefficients \(\mathbf e \in \mathbb R^{N_1}\) of \(\mathbf E \in V^1_h\) and \(\mathbf b \in \mathbb R^{N_2}\) of \(\mathbf B \in V^2_h\), the semi-discrete system reads

\[\mathbb M_1 \dot{\mathbf e} = \mathbb C^\top \mathbb M_2 \mathbf b\,, \qquad \dot{\mathbf b} = - \mathbb C \mathbf e\,,\]

with the mass matrices \(\mathbb M_1\), \(\mathbb M_2\) (WeightedMassOperators) and the discrete curl \(\mathbb C\) (curl). It conserves the energy \(\frac 12 \mathbf e^\top \mathbb M_1 \mathbf e + \frac 12 \mathbf b^\top \mathbb M_2 \mathbf b\) and the discrete divergence \(\mathbb D \mathbf b\) (since \(\mathbb D \mathbb C = 0\)).

Time discretization:

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

  \[\mathbb M_1 \frac{\mathbf e^{n+1} - \mathbf e^n}{\Delta t} = \mathbb C^\top \mathbb M_2 \frac{\mathbf b^{n+1} + \mathbf b^n}{2}\,, \qquad \frac{\mathbf b^{n+1} - \mathbf b^n}{\Delta t} = - \mathbb C \frac{\mathbf e^{n+1} + \mathbf e^n}{2}\,,\]

  which conserves both invariants exactly (up to the solver tolerance). The magnetic field is eliminated by SchurSolver with \(x = \mathbf e\), \(y = \mathbf b\), \(A = \mathbb M_1\), \(B = -\frac 12 \mathbb C^\top \mathbb M_2\) and \(C = \frac 12 \mathbb C\), which gives the curl-curl system

  \[\left(\mathbb M_1 + \frac{\Delta t^2}{4} \mathbb C^\top \mathbb M_2 \mathbb C\right) \mathbf e^{n+1} = \left(\mathbb M_1 - \frac{\Delta t^2}{4} \mathbb C^\top \mathbb M_2 \mathbb C\right) \mathbf e^n + \Delta t \, \mathbb C^\top \mathbb M_2 \mathbf b^n\,,\]

  solved with solver and precond. By default (precond="StiffnessPreconditioner") the Schur operator is preconditioned by StiffnessPreconditioner for \(\mathbb C^\top \left(\frac{\Delta t^2}{4} \mathbb M_2\right) \mathbb C + \mathbb M_1\), built for the current time step (prepare_time_step()). Then \(\mathbf b^{n+1} = \mathbf b^n - \frac{\Delta t}{2} \mathbb C (\mathbf e^{n+1} + \mathbf e^n)\).

• explicit: Runge-Kutta method from ButcherTableau for the system \(\dot{\mathbf e} = \mathbb M_1^{-1} \mathbb C^\top \mathbb M_2 \mathbf b\), \(\dot{\mathbf b} = - \mathbb C \mathbf e\), where \(\mathbb M_1^{-1}\) is applied with solver and precond in each stage ("StiffnessPreconditioner" is replaced by the mass matrix preconditioner of \(\mathbb M_1\) here). The time step is limited by the CFL condition of light waves.

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['StiffnessPreconditioner', 'MassMatrixPreconditioner', 'PolarMassMatrixPreconditioner', None] = 'StiffnessPreconditioner', 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 feectools.linalg.solvers.inverse.

• precond (LiteralOptions.OptsMaxwellPrecond, default="StiffnessPreconditioner") –

  Preconditioner of the linear solves.

  • "StiffnessPreconditioner": Kronecker (fast diagonalization) approximation of the inverse of the Schur operator \(\mathbb M_1 + \frac{\Delta t^2}{4} \mathbb C^\top \mathbb M_2 \mathbb C\) (StiffnessPreconditioner with derivative="curl"; with polar splines a PolarStiffnessPreconditioner), rebuilt when the time step changes. For algo="explicit", the mass matrix preconditioner of \(\mathbb M_1\) is used instead. (With polar splines, the implicit step currently fails in the SchurSolver with any preconditioner.)

  • "MassMatrixPreconditioner", "PolarMassMatrixPreconditioner": preconditioner of \(\mathbb M_1\), better for small time steps \(\Delta t \lesssim h\), e.g. in CFL-limited PIC runs.

  • None: no preconditioner (set automatically, with a warning, for solver="cg", which takes no preconditioner).

• solver_params (SolverParameters, default=None) – Iterative-solver controls (for example tol, maxiter, 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.

prepare_time_step(dt: float)

Stiffness preconditioner of the Schur operator M1 + dt^2/4 C^T M2 C for the time step dt: built in the setup, rebuilt when dt changes.

Particle and FEEC variables are updated.

class struphy.propagators.vlasov_maxwell_coupling.VlasovMaxwellCoupling

Bases: Propagator

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

\[\begin{split}\lambda^2 \int_{\Omega} \frac{\partial \mathbf E}{\partial t} \cdot \mathbf F \,\textrm d \mathbf x &= c_1 \int_{\Omega} \mathbf B \cdot \nabla \times \mathbf F \,\textrm d \mathbf x - c_2 \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 \mathbf B}{\partial t} + \nabla\times\mathbf E &= 0\,, \\[2mm] \frac{\partial f}{\partial t} + c_3\, \mathbf{E} \cdot \frac{\partial f}{\partial \mathbf{v}} &= 0\,,\end{split}\]

with \(\lambda \geq 0\), \(c_1 > 0\) and \(c_2 c_3 > 0\). It solves the steps of MaxwellWeakAmpere and VlasovAmpereCoupling in one implicit step (for \(\lambda = 1\), \(c_1 = 1\), \(c_2 = \alpha^2/\varepsilon\), \(c_3 = 1/\varepsilon\), the light velocity scale). The defaults are the coefficients of the species’ equation parameters, \(\lambda = \varepsilon \lambda_\textrm{s}\), \(c_2 = \kappa_\textrm{s}^2 \varepsilon\), \(c_3 = 1/\varepsilon\), i.e. \(\lambda = \lambda_\textrm{s}\), \(c_2 = \kappa_\textrm{s}^2\), \(c_3 = 1\) for the cyclotron velocity scale. For \(\lambda = 0\) the displacement current vanishes, \(c_1 \nabla \times \mathbf B = c_2 \mathbf J\) (quasi-neutral limit).

With the FE coefficients \(\mathbf e \in \mathbb R^{N_1}\) of \(\mathbf E \in V^1_h\) and \(\mathbf b \in \mathbb R^{N_2}\) of \(\mathbf B \in V^2_h\), and \(f\) approximated by \(N\) markers with positions \(\boldsymbol \eta_p\) (logical, fixed in this step), velocities \(\mathbf V = (\mathbf v_p)_p \in \mathbb R^{3N}\) and weights \(\mathbb W = \textnormal{diag}(w_p)\), the semi-discrete system reads

\[\lambda^2\, \mathbb M_1 \dot{\mathbf e} = c_1\, \mathbb C^\top \mathbb M_2 \mathbf b - c_2\, (\mathbf \Lambda^1)^\top \, \overline{DF^{-1}} \, \mathbb W \, \mathbf V \,, \qquad \dot{\mathbf b} = - \mathbb C \mathbf e\,, \qquad \dot{\mathbf V} = c_3\, \overline{DF^{-\top}} \, \mathbf \Lambda^1 \, \mathbf e\,,\]

with the notation of VlasovAmpereCoupling. It conserves the energy (energy())

\[H = \frac{\lambda^2}{2}\, \mathbf e^\top \mathbb M_1 \mathbf e + \frac{c_1}{2}\, \mathbf b^\top \mathbb M_2 \mathbf b + \frac{c_2}{2 c_3}\, \mathbf V^\top \mathbb W \mathbf V\,.\]

Time discretization: Crank-Nicolson (implicit mid-point), with \(\mathbf e^\textnormal{mid} = \frac 12 (\mathbf e^{n+1} + \mathbf e^n)\) etc.,

\[\lambda^2\, \mathbb M_1 \frac{\mathbf e^{n+1} - \mathbf e^n}{\Delta t} = c_1\, \mathbb C^\top \mathbb M_2 \mathbf b^\textnormal{mid} - c_2\, (\mathbf \Lambda^1)^\top \, \overline{DF^{-1}} \, \mathbb W \, \mathbf V^\textnormal{mid}\,, \qquad \frac{\mathbf b^{n+1} - \mathbf b^n}{\Delta t} = - \mathbb C \mathbf e^\textnormal{mid}\,, \qquad \frac{\mathbf V^{n+1} - \mathbf V^n}{\Delta t} = c_3\, \overline{DF^{-\top}} \, \mathbf \Lambda^1 \, \mathbf e^\textnormal{mid}\,,\]

which conserves \(H\) exactly (up to the solver tolerance; for constant weights, i.e. without control variate). Eliminating \(\mathbf b^\textnormal{mid}\) and \(\mathbf V^\textnormal{mid}\) gives the Schur system for the mid-point field

\[\left(\lambda^2\, \mathbb M_1 + \frac{\Delta t^2}{4}\, \mathbb K\right) \mathbf e^\textnormal{mid} = \lambda^2\, \mathbb M_1 \mathbf e^n + \frac{\Delta t}{2} \left( c_1\, \mathbb C^\top \mathbb M_2 \mathbf b^n - c_2\, \mathbf j \right)\,, \qquad \mathbb K = c_1\, \mathbb C^\top \mathbb M_2 \mathbb C + c_2 c_3\, \mathbb A\,,\]

with the accumulated matrix \(\mathbb A\) and vector \(\mathbf j\) of VlasovAmpereCoupling (kernel vlasov_maxwell()). It is solved with solver and precond. Then

\[\mathbf b^{n+1} = \mathbf b^n - \Delta t\, \mathbb C \mathbf e^\textnormal{mid}\,, \qquad \mathbf e^{n+1} = 2 \mathbf e^\textnormal{mid} - \mathbf e^n \quad (\lambda > 0)\,,\]

and the velocities are updated by push_v_with_efield() with \(\mathbf e^\textnormal{mid}\).

Quasi-neutral limit \(\lambda = 0\). The curl-curl operator \(\mathbb C^\top \mathbb M_2 \mathbb C\) is singular on the discrete gradients (and harmonic fields); this kernel is absorbed by the plasma term \(c_2 c_3 \mathbb A\), a particle approximation of the density-weighted mass matrix \(\mathbb M_1\). The Schur operator is therefore positive definite for positive weights (full-f, no control variate) and plasma everywhere in the domain. The scheme then imposes \(c_1\, \mathbb C^\top \mathbb M_2 \mathbf b^\textnormal{mid} = c_2\, (\mathbf \Lambda^1)^\top \overline{DF^{-1}} \mathbb W \mathbf V^\textnormal{mid}\), hence a weakly divergence-free mid-point current (multiply by \(\mathbb G^\top\)); the gradient part of \(\mathbf E\) is its Lagrange multiplier. \(\mathbf e^n\) drops out of the system and only \(\mathbf e^\textnormal{mid}\) is determined: the propagator stores \(\mathbf e^{n+1} = \mathbf e^\textnormal{mid}\) (the Crank-Nicolson update \(2 \mathbf e^\textnormal{mid} - \mathbf e^n\) would carry any error of \(\mathbf e^n\) as an undamped even-odd oscillation). The energy \(H\) (without the electric part) is still conserved exactly.

The constraint holds at the mid-point only: its residual \(\mathbf r^n = c_1 \mathbb C^\top \mathbb M_2 \mathbf b^n - c_2 \mathbf j(\mathbf V^n)\) (linear in \(\mathbf V\) at fixed positions) satisfies \(\mathbf r^{n+1} = - \mathbf r^n\). Initial data that violate the constraint (e.g. by marker noise in \(\mathbf J\)) would thus oscillate with period 2 without damping. With initial_projection (default), the setup replaces \((\mathbf b^0, \mathbf V^0)\) by the mid-point state of one step, \(\mathbf b^0 - \frac{\Delta t}{2} \mathbb C \mathbf e^\textnormal{mid}\) and \(\mathbf V^0 + \frac{\Delta t}{2} c_3 \overline{DF^{-\top}} \mathbf \Lambda^1 \mathbf e^\textnormal{mid}\) (an implicit Euler step of \(\Delta t / 2\)), which satisfies the constraint exactly and does not depend on \(\Delta t\): a projection onto the constraint. The particle pushers of the splitting change \(\mathbf r\) again by \(O(\Delta t)\) in every step.

Preconditioner. With \(\mathbb A \approx \bar n\, \mathbb M_1\), the Schur operator is approximately the stabilized curl-curl operator \(\mathbb C^\top \left(\frac{\Delta t^2 c_1}{4} \mathbb M_2\right) \mathbb C + \left(\lambda^2 + \frac{\Delta t^2 c_2 c_3 \bar n}{4}\right) \mathbb M_1\), preconditioned (default) by StiffnessPreconditioner, built for the current time step (prepare_time_step()). Its stabilization stays positive for \(\lambda = 0\).

class Variables

Bases: object

Container for variables advanced by VlasovMaxwellCoupling.

Variables:

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

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

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

class Options(lam: float | None = None, c1: float = 1.0, c2: float | None = None, c3: float | None = None, solver: Literal['pcg', 'cg'] = 'pcg', precond: Literal['StiffnessPreconditioner', 'MassMatrixPreconditioner', 'PolarMassMatrixPreconditioner', None] = 'StiffnessPreconditioner', solver_params: SolverParameters | None = None, n_ref: float | None = None, matrix_free: bool = False, initial_projection: bool = True)

Bases: OptionsBase

Configuration options for VlasovMaxwellCoupling.

Parameters:

• lam (float, default=None) – \(\lambda \geq 0\), with the coefficient \(\lambda^2\) of the displacement current \(\partial_t \mathbf E\); lam=0.0 is the quasi-neutral limit (needs full-f weights, no control variate). If None, \(\varepsilon \lambda_\textrm{s}\) of the species (lam and epsilon of its equation parameters).

• c1 (float, default=1.0) – Coefficient \(c_1 > 0\) of \(\nabla \times \mathbf B\) in Ampère’s law.

• c2 (float, default=None) – Coefficient of the current density in Ampère’s law. If None, \(\kappa_\textrm{s}^2 \varepsilon\) of the species (kappa and epsilon of its equation parameters).

• c3 (float, default=None) – Coefficient of the electric force in the Vlasov equation (\(c_2 c_3 > 0\)). If None, \(1/\varepsilon\) of the species.

• solver (LiteralOptions.OptsSymmSolver, default="pcg") – Symmetric iterative solver of the Schur system, passed to feectools.linalg.solvers.inverse.

• precond (LiteralOptions.OptsMaxwellPrecond, default="StiffnessPreconditioner") –

  Preconditioner of the Schur system.

  • "StiffnessPreconditioner": StiffnessPreconditioner with derivative="curl" for the stabilized curl-curl approximation of the Schur operator (see the class docstring), rebuilt when the time step changes.

  • "MassMatrixPreconditioner", "PolarMassMatrixPreconditioner": preconditioner of \(\mathbb M_1\), sufficient for \(\Delta t \lesssim h\) (polar splines are not supported by this propagator yet).

  • None: no preconditioner (set automatically, with a warning, for solver="cg").

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

• n_ref (float, default=None) – Reference density \(\bar n\) in the stiffness preconditioner. If None, the Rayleigh quotient \(\mathbf 1^\top \mathbb A \mathbf 1 / \mathbf 1^\top \mathbb M_1 \mathbf 1\) of the accumulated matrix at the setup (a density average).

• matrix_free (bool, default=False) – Apply the particle matrix \(\mathbb A\) matrix-free, with the filling \(w_p G^{-1}(\boldsymbol \eta_p)\) stored per marker (MarkerOperator), instead of assembling it: one spline evaluation and one deposition per marker and solver iteration instead of the products of all pairs of basis functions per marker. Needs 6 auxiliary marker columns. Ignored on the CuPy backend.

• initial_projection (bool, default=True) – Only for lam=0: in the setup (first prepare_time_step()), replace \((\mathbf b, \mathbf V)\) by the mid-point state of one step, which satisfies \(c_1 \nabla \times \mathbf B = c_2 \mathbf J\) exactly (see the class docstring). Without it, initial data that violate the constraint lead to an undamped period-2 oscillation of \(\mathbf b\) and \(\mathbf V\). Not applied to a restored (restart) state.

prepare_time_step(dt: float)

Schur operator \(\lambda^2 \mathbb M_1 + \frac{\Delta t^2}{4} \mathbb K\), its preconditioner and solver for the time step dt: built in the setup, rebuilt when dt changes. The particle matrix enters by reference (assembled, or matrix-free reading the fillings stored in the markers), so it is updated by the accumulation in every step. For lam=0 (and initial_projection), the first call also projects the state onto the constraint \(c_1 \mathbb C^\top \mathbb M_2 \mathbf b = c_2 \mathbf j\).

energies() → tuple[float, float, float]

The three parts \(\frac{\lambda^2}{2} \mathbf e^\top \mathbb M_1 \mathbf e\), \(\frac{c_1}{2} \mathbf b^\top \mathbb M_2 \mathbf b\) and \(\frac{c_2}{2 c_3} \mathbf V^\top \mathbb W \mathbf V\) of the energy energy() (summed over all MPI ranks; a collective operation).

energy() → float

The energy \(H = \frac{\lambda^2}{2} \mathbf e^\top \mathbb M_1 \mathbf e + \frac{c_1}{2} \mathbf b^\top \mathbb M_2 \mathbf b + \frac{c_2}{2 c_3} \mathbf V^\top \mathbb W \mathbf V\) conserved by the propagator (summed over all MPI ranks; a collective operation).

Particle and FEEC variables are updated.

class struphy.propagators.drift_kinetic_maxwell_coupling.DriftKineticMaxwellCoupling

Bases: Propagator

PIC-FEEC discretization of the following equations: find \(\mathbf E \in H(\textnormal{curl})\), \(\tilde{\mathbf B} \in H(\textnormal{div})\) and the guiding-center distribution \(f(\mathbf X, v_\parallel, \mu)\) (positions fixed in this step) such that

\[\begin{split}\lambda^2 \int_{\Omega} \frac{\partial \mathbf E}{\partial t} \cdot \mathbf F \,\textrm d \mathbf x + \chi \int_{\Omega} \frac{n_0}{|B_0|^2} \frac{\partial \mathbf E_\perp}{\partial t} \cdot \mathbf F_\perp \,\textrm d \mathbf x &= c_1 \int_{\Omega} \tilde{\mathbf B} \cdot \nabla \times \mathbf F \,\textrm d \mathbf x - c_2 \int_{\Omega} \iint f^\textnormal{vol} \left( v_\parallel \frac{\mathbf B^*}{B^*_\parallel} + \frac{\mathbf E \times \mathbf b_0}{B^*_\parallel} \right) \cdot \mathbf F \,\textrm d v_\parallel \textrm d \mu \,\textrm d \mathbf x \qquad \forall \, \mathbf F \in H(\textnormal{curl})\,, \\[2mm] \frac{\partial \tilde{\mathbf B}}{\partial t} + \nabla\times\mathbf E &= 0\,, \\[2mm] \frac{\partial f}{\partial t} + c_3\, \frac{\mathbf B^*}{B^*_\parallel} \cdot \mathbf E\, \frac{\partial f}{\partial v_\parallel} &= 0\,,\end{split}\]

with \(\mathbf E_\perp = \mathbf E - (\mathbf b_0 \cdot \mathbf E) \mathbf b_0\), the phase-space density \(f^\textnormal{vol} = B^*_\parallel f\) (the marker weights) and, as in PushGuidingCenterParallel,

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

with the total magnetic field in the first term (parallel streaming along the perturbed field lines) and the curvature term of the equilibrium.

The current is the guiding-center current \(\iint f^\textnormal{vol} \dot{\mathbf X}\) without the magnetic drifts (which move the positions): parallel streaming, including the curvature drift in \(\mathbf B^*\), and (option exb_current) the \(E \times B\) drift. The coefficients are \(\lambda \geq 0\), \(\chi \geq 0\) with \(\lambda^2 + \chi > 0\), \(c_1 > 0\) and \(c_2 c_3 > 0\). The defaults are the coefficients of the species’ equation parameters, \(\lambda = \varepsilon \lambda_\textrm{s}\), \(\chi = (\varepsilon \kappa_\textrm{s})^2\), \(c_1 = 1\), \(c_2 = \kappa_\textrm{s}^2 \varepsilon\), \(c_3 = 1/\varepsilon\), i.e. for the thermal velocity scale of DriftKineticMaxwell \(\lambda = \hat v_\textrm{th}/c\), \(\chi = \hat v_\textrm{th}^2/\hat v_A^2\). The \(\chi\) term is the (linearized, long-wavelength) polarization current of the guiding centers, which the drift-kinetic ordering removes from \(\dot{\mathbf X}\): without it, the perpendicular electrostatic field is undetermined for \(\lambda = 0\) (flute modes), and with only the displacement current the shear Alfvén wave travels at the speed of light.

With the FE coefficients \(\mathbf e \in \mathbb R^{N_1}\) of \(\mathbf E \in V^1_h\) and \(\mathbf b \in \mathbb R^{N_2}\) of \(\tilde{\mathbf B} \in V^2_h\), and \(N\) markers with logical positions \(\boldsymbol \eta_p\) (fixed), parallel velocities \(v_{\parallel,p}\) and weights \(w_p\), the semi-discrete system reads

\[\mathbb M_E \dot{\mathbf e} = c_1\, \mathbb C^\top \mathbb M_2 \mathbf b - c_2 \left( \mathbf j(\mathbf V) + \mathbb S\, \mathbf e \right) \,, \qquad \dot{\mathbf b} = - \mathbb C \mathbf e\,, \qquad \dot v_{\parallel,p} = c_3\, \mathbf a_p \cdot \mathbf \Lambda^1_p \mathbf e\,,\]

with \(\mathbb M_E = \lambda^2\, \mathbb M_1 + \chi\, \mathbb M_1^\perp\) (\(\mathbb M_1^\perp\) = M1gyro), \(\mathbf a_p = \hat{\mathbf B}^{*2} / \hat B^{*3}_\parallel (\boldsymbol \eta_p)\), \(\mathbf \Lambda^1_p\) the basis of \(V^1_h\) at \(\boldsymbol \eta_p\), and

\[\mathbf j(\mathbf V) = \sum_p w_p v_{\parallel,p}\, (\mathbf \Lambda^1_p)^\top \mathbf a_p\,, \qquad \mathbb A_\parallel = \sum_p w_p\, (\mathbf \Lambda^1_p)^\top \mathbf a_p \mathbf a_p^\top \mathbf \Lambda^1_p\,, \qquad \mathbf f^\top \mathbb S\, \mathbf e = \sum_p \frac{w_p}{\hat B^{*3}_\parallel(\boldsymbol \eta_p)}\, (\mathbf \Lambda^1_p \mathbf f) \cdot \big( (\mathbf \Lambda^1_p \mathbf e) \times \hat{\mathbf b}^1_0 (\boldsymbol \eta_p) \big)\,.\]

\(\mathbb A_\parallel\) is symmetric positive semi-definite (the parallel response), \(\mathbb S\) is antisymmetric (the \(E \times B\) current does no work). The system conserves the energy

\[H = \frac{1}{2}\, \mathbf e^\top \mathbb M_E \mathbf e + \frac{c_1}{2}\, \mathbf b^\top \mathbb M_2 \mathbf b + \frac{c_2}{c_3} \sum_p w_p \left( \frac{v_{\parallel,p}^2}{2} + \mu_p |B_0|(\boldsymbol \eta_p) \right) \qquad (:meth:`~struphy.propagators.drift_kinetic_maxwell_coupling.DriftKineticMaxwellCoupling.energy`)\,.\]

Time discretization: Crank-Nicolson (implicit mid-point), with \(\mathbf a_p\) and \(\mathbb S\) evaluated at the beginning of the step, which conserves \(H\) exactly (up to the solver tolerance; for constant weights). Eliminating \(\mathbf b^\textnormal{mid}\) and \(\mathbf V^\textnormal{mid}\) gives the Schur system for the mid-point field

\[\left( \mathbb M_E + \frac{\Delta t^2}{4} \left( c_1\, \mathbb C^\top \mathbb M_2 \mathbb C + c_2 c_3\, \mathbb A_\parallel \right) + \frac{\Delta t}{2}\, c_2\, \mathbb S \right) \mathbf e^\textnormal{mid} = \mathbb M_E \mathbf e^n + \frac{\Delta t}{2} \left( c_1\, \mathbb C^\top \mathbb M_2 \mathbf b^n - c_2\, \mathbf j(\mathbf V^n) \right)\,,\]

accumulated by dk_maxwell() (or matrix-free) and solved with solver and precond (non-symmetric with the \(E \times B\) current). Then

\[\mathbf b^{n+1} = \mathbf b^n - \Delta t\, \mathbb C \mathbf e^\textnormal{mid}\,, \qquad v_{\parallel,p}^{n+1} = v_{\parallel,p}^n + \Delta t\, c_3\, \mathbf a_p \cdot \mathbf \Lambda^1_p \mathbf e^\textnormal{mid}\,, \qquad \mathbf e^{n+1} = 2 \mathbf e^\textnormal{mid} - \mathbf e^n\,,\]

with the velocity update by push_vpar_with_efield().

Quasi-neutral limit \(\lambda = 0\). The parallel part of \(\mathbf e\) is then determined by \(\mathbb A_\parallel\) alone (which needs positive weights, i.e. full-f, and markers everywhere), the perpendicular part by \(\chi\, \mathbb M_1^\perp\) and the curl-curl operator. The perpendicular field still carries the energy \(\frac 12 \mathbf e^\top \mathbb M_E \mathbf e\), so the Crank-Nicolson update is kept (storing \(\mathbf e^\textnormal{mid}\) would dissipate it). The parallel part of the stored \(\mathbf e^{n+1}\) is not determined by the scheme and may oscillate from step to step; it does not feed back, since the steps see \(\mathbf e^n\) only through \(\mathbb M_E \mathbf e^n\) and the \(E \times B\) push only through \(\mathbf E_\perp\). With initial_projection (default), the setup replaces \((\mathbf b^0, \mathbf V^0)\) by the mid-point state of one step, which removes the undamped period-2 oscillation of the constraint residual.

Preconditioner. With \(\mathbb A_\parallel \approx \bar n\, \mathbb M_1^\parallel\) (M1para), the symmetric part of the Schur operator is approximately \(\mathbb C^\top \left(\frac{\Delta t^2 c_1}{4} \mathbb M_2\right) \mathbb C + \mathbb M_E + \frac{\Delta t^2 c_2 c_3 \bar n}{4} \mathbb M_1^\parallel\), preconditioned (default) by StiffnessPreconditioner, built for the current time step (prepare_time_step()).

class Variables

Bases: object

Container for variables advanced by DriftKineticMaxwellCoupling.

Variables:

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

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

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

class Options(lam: float | None = None, chi: float | None = None, c1: float = 1.0, c2: float | None = None, c3: float | None = None, exb_current: bool = True, solver: Literal['pbicgstab', 'bicgstab', 'gmres', 'pcg', 'cg'] = 'pbicgstab', precond: Literal['StiffnessPreconditioner', 'MassMatrixPreconditioner', 'PolarMassMatrixPreconditioner', None] = 'StiffnessPreconditioner', solver_params: SolverParameters | None = None, n_ref: float | None = None, matrix_free: bool = False, initial_projection: bool = True)

Bases: OptionsBase

Configuration options for DriftKineticMaxwellCoupling.

Parameters:

• lam (float, default=None) – The displacement current has the coefficient \(\lambda^2\) (lam \(= \lambda \geq 0\)); lam=0.0 is the quasi-neutral limit (needs chi > 0 and full-f weights, no control variate). If None, \(\varepsilon \lambda_\textrm{s}\) of the species (lam and epsilon of its equation parameters).

• chi (float, default=None) – Coefficient \(\chi \geq 0\) of the polarization current \(\frac{n_0}{|B_0|^2} \partial_t \mathbf E_\perp\). If None, \((\varepsilon \kappa_\textrm{s})^2\) of the species (the polarization of the guiding centers themselves); a larger value mimics the polarization of a heavier species. lam**2 + chi must be positive. The density \(n_0\) is the one of the MHD equilibrium (M1gyro), not of the kinetic background: choose them consistently.

• c1 (float, default=1.0) – Coefficient \(c_1 > 0\) of \(\nabla \times \tilde{\mathbf B}\) in Ampère’s law.

• c2 (float, default=None) – Coefficient of the current density in Ampère’s law. If None, \(\kappa_\textrm{s}^2 \varepsilon\) of the species (kappa and epsilon of its equation parameters).

• c3 (float, default=None) – Coefficient of the parallel electric force (\(c_2 c_3 > 0\)). If None, \(1/\varepsilon\) of the species.

• exb_current (bool, default=True) – Include the \(E \times B\) current (antisymmetric operator \(\mathbb S\)). The Schur system is then non-symmetric and needs a general solver. Switch it off if the \(E \times B\) current is compensated (e.g. by other species).

• solver (OptsSolver, default="pbicgstab") – Iterative solver of the Schur system, passed to feectools.linalg.solvers.inverse. The symmetric solvers "pcg" and "cg" need exb_current=False.

• precond (LiteralOptions.OptsMaxwellPrecond, default="StiffnessPreconditioner") –

  Preconditioner of the Schur system (only for "pcg" and "pbicgstab"; set to None otherwise).

  • "StiffnessPreconditioner": StiffnessPreconditioner with derivative="curl" for the stabilized curl-curl approximation of the Schur operator (see the class docstring), rebuilt when the time step changes.

  • "MassMatrixPreconditioner", "PolarMassMatrixPreconditioner": preconditioner of \(\mathbb M_1\).

  • None: no preconditioner.

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

• n_ref (float, default=None) – Reference density \(\bar n\) in the stiffness preconditioner. If None, the Rayleigh quotient \(\mathbf 1^\top \mathbb A_\parallel \mathbf 1 / \mathbf 1^\top \mathbb M_1^\parallel \mathbf 1\) at the setup (a density average).

• matrix_free (bool, default=False) – Apply the particle matrices \(\mathbb A_\parallel\) and \(\mathbb S\) matrix-free, with their fillings stored per marker (MarkerOperator), instead of assembling them. Needs 7 auxiliary marker columns.

• initial_projection (bool, default=True) – Only for lam=0: in the setup (first prepare_time_step()), replace \((\mathbf b, \mathbf V)\) by the mid-point state of one step, which satisfies the constraint of the quasi-neutral limit exactly. Not applied to a restored (restart) state.

prepare_time_step(dt: float)

Schur operator, its preconditioner and solver for the time step dt: built in the setup, rebuilt when dt changes. The particle matrices enter by reference (assembled, or matrix-free reading the fillings stored in the markers), so they are updated by the accumulation in every step. For lam=0 (and initial_projection), the first call also projects the state onto the constraint.

energies() → tuple[float, float, float]

The three parts \(\frac{1}{2} \mathbf e^\top \mathbb M_E \mathbf e\) (displacement and polarization), \(\frac{c_1}{2} \mathbf b^\top \mathbb M_2 \mathbf b\) and \(\frac{c_2}{c_3} \sum_p w_p (v_{\parallel,p}^2/2 + \mu_p |B_0|(\boldsymbol \eta_p))\) of the energy energy() (summed over all MPI ranks; a collective operation).

energy() → float

The energy \(H = \frac{1}{2} \mathbf e^\top \mathbb M_E \mathbf e + \frac{c_1}{2} \mathbf b^\top \mathbb M_2 \mathbf b + \frac{c_2}{c_3} \sum_p w_p (v_{\parallel,p}^2/2 + \mu_p |B_0|(\boldsymbol \eta_p))\) conserved by the propagator (summed over all MPI ranks; a collective operation).

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 | 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, matrix_free: bool = True)

Bases: OptionsBase

Configuration options for CurrentCoupling5DCurlb.

Parameters:

• ep_scale (float, default=None) – Factor of the energetic-particle contributions (accumulation): the mass ratio \(A_\textnormal{h}/A_\textnormal{b}\) of the energetic ions and the MHD bulk if None.

• 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) – Filter of the accumulated (particle-to-grid) vectors; defaults to FilterParameters() (no filter). The three-point (binomial) filter FilterParameters(use_filter="three_point") reduces the noise, but the propagator then no longer conserves the energy exactly (the filter is applied to the accumulated vectors only, not to the fields seen by the particles).

• matrix_free (bool, default=True) – Apply the accumulated matrix \(\sum_p c_p (\vec \Lambda_p \cdot \mathbf e_p) (\vec \Lambda_p \cdot \mathbf e_p)^\top\) (rank one per marker) in each iteration of the solver from data stored per marker (MarkerOperator), instead of assembling it: one spline evaluation and one deposition per marker and iteration instead of the 6 blocks of size \(n_\textnormal{basis}^2\) per marker (much cheaper for cubic splines). Same result up to the solver tolerance.

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 | None = None, 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, matrix_free: bool = True)

Bases: OptionsBase

Configuration options for CurrentCoupling5DDensity.

Parameters:

• ep_scale (float, default=None) – Factor of the energetic-particle contributions (accumulation): the mass ratio \(A_\textnormal{h}/A_\textnormal{b}\) of the energetic ions and the MHD bulk if None.

• 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) – Filter of the accumulated (particle-to-grid) vectors; defaults to FilterParameters() (no filter). The three-point (binomial) filter FilterParameters(use_filter="three_point") reduces the noise, but the propagator then no longer conserves the energy exactly (the filter is applied to the accumulated vectors only, not to the fields seen by the particles).

• matrix_free (bool, default=True) – Apply the particle matrix \(\sum_p \vec \Lambda_p F_p \vec \Lambda_p^\top\) (antisymmetric 3x3 \(F_p\)) in each iteration of the solver from data stored per marker (MarkerOperator) instead of assembling it: one spline evaluation and one deposition per marker and iteration instead of 3 blocks of size \(n_\textnormal{basis}^2\) per marker. Same result up to the solver tolerance.

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 | None = None, 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=None) – Factor of the energetic-particle contributions (accumulation): the mass ratio \(A_\textnormal{h}/A_\textnormal{b}\) of the energetic ions and the MHD bulk if None.

• 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) – Filter of the accumulated (particle-to-grid) vectors; defaults to FilterParameters() (no filter). The three-point (binomial) filter FilterParameters(use_filter="three_point") reduces the noise, but the propagator then no longer conserves the energy exactly (the filter is applied to the accumulated vectors only, not to the fields seen by the particles).

• 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), matrix_free: bool = False)

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.

• matrix_free (bool, default=False) – Apply the particle matrix matrix-free, with its symmetric filling stored per marker (MarkerOperator), instead of assembling it: one spline evaluation and one deposition per marker and solver iteration instead of the products of all pairs of basis functions per marker. Needs 6 auxiliary marker columns. Ignored on the CuPy backend.

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['pbicgstab', 'bicgstab', 'gmres'] = 'pbicgstab', 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), matrix_free: bool = False)

Bases: OptionsBase

Configuration options for CurrentCoupling6DDensity.

Parameters:

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

• solver (LiteralOptions.OptsGenSolver, default="pbicgstab") – Iterative solver for the linear system, which is not symmetric (the density matrix is antisymmetric).

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

• matrix_free (bool, default=False) – Apply the density matrix matrix-free, with its antisymmetric filling stored per marker (MarkerOperator), instead of assembling it: one spline evaluation and one deposition per marker and solver iteration instead of the products of all pairs of basis functions per marker. Needs 3 auxiliary marker columns. Ignored on the CuPy backend.

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^1 \left[\frac{\mathbf B^2}{\sqrt{g}} \times \vec \Lambda^2\right]\) with the total field \(\mathbf B = \mathbf B_0 + \tilde{\mathbf B}\) (updated in each step) if nonlinear, else the equilibrium field \(\mathbf B_0\), 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 | None = None, u_space: Literal['Hcurl', 'Hdiv', 'H1vec'] = 'Hdiv', algo: Literal['implicit', 'explicit'] = 'implicit', solver: Literal['pcg', 'cg'] = 'pcg', precond: Literal['MassMatrixPreconditioner', 'PolarMassMatrixPreconditioner', 'MultiGrid', 'StiffnessPreconditioner', None] = 'StiffnessPreconditioner', exb_correction: bool = True, 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=None) – Factor of the energetic-particle contributions (accumulation): the mass ratio \(A_\textnormal{h}/A_\textnormal{b}\) of the energetic ions and the MHD bulk if None.

• 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.OptsShearAlfvenPrecond, default="StiffnessPreconditioner") – Preconditioner of the linear solve: of the mass-matrix block ("MassMatrixPreconditioner", "PolarMassMatrixPreconditioner") or the Kronecker (fast diagonalization) preconditioner of the whole Schur operator ("StiffnessPreconditioner", the default, ShearAlfvenSchurPreconditioner, built from the equilibrium part \(\mathbb T\) of the basis projection operator; with nonlinear=True it approximates the Schur operator with \(\mathbb T + \mathbb T^B\), since \(|\tilde{\mathbf B}| \ll |\mathbf B_0|\)). The Kronecker preconditioner needs algo="implicit", solver="pcg", u_space="Hdiv", tensor-product splines (no polar splines) and the global commuting projectors; otherwise it falls back (with a warning) to "MassMatrixPreconditioner" ("PolarMassMatrixPreconditioner" with polar splines). "MultiGrid" is not available here.

• exb_correction (bool, default=True) – For precond="StiffnessPreconditioner": add the auxiliary space correction of the \(\mathbf E \times \mathbf B\) flows, see ShearAlfvenPropagator.

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

• filter_params (FilterParameters, default=None) – Filter of the accumulated (particle-to-grid) vectors; defaults to FilterParameters() (no filter). The three-point (binomial) filter FilterParameters(use_filter="three_point") reduces the noise, but the propagator then no longer conserves the energy exactly (the filter is applied to the accumulated vectors only, not to the fields seen by the particles).

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

• nonlinear (bool, default=True) – If True, the basis projection operator \(\mathcal T^2\) uses the total magnetic field \(\mathbf B_0 + \tilde{\mathbf B}\) (one operator, re-weighted in each step), else the equilibrium field \(\mathbf B_0\) only.

prepare_time_step(dt: float)

Kronecker preconditioner of the Schur operator for the time step dt: built in the setup, rebuilt when dt changes.

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, matrix_free: 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.

• matrix_free (bool, default=True) – Apply the 3x3 block particle matrix matrix-free, with its filling \(w_p G^{-1}(oldsymbol \eta_p)\) stored per marker (MarkerOperator with n_blocks=3), instead of assembling its six symmetric blocks: per marker and solver iteration, one evaluation of each block and one deposition into each block instead of the products of all pairs of basis functions per marker. Needs 6 auxiliary marker columns. Ignored on the CuPy backend.

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 | MarkerOperator) – 3x3 nested list of the accumulated pressure terms (six of them distinct), or the matrix-free MarkerOperator of the whole 3x3 block matrix (n_blocks=3).

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.