Models

Full-orbit Vlasov

class struphy.models.VlasovAmpereOneSpecies(base_units: BaseUnits = BaseUnits(x=1.0, B=1.0, n=1.0, kBT=None), charge_number: int = 1, mass_number: float = 1.0, alpha: float | None = None, epsilon: float | None = None, with_B0: bool = True)

Bases: StruphyModel

Vlasov-Ampère system for a single kinetic species in an electric field.

Parameters:

• base_units (BaseUnits) – Base units for normalization (default: BaseUnits())

• charge_number (int) – Charge number (in units of the positive elementary charge) of the ion species (default: 1)

• mass_number (float) – Mass number (in units or Proton mass) of the ion species (default: 1.0)

• alpha (float, optional) – Dimensionless parameter: plasma frequency / cyclotron frequency. If None, computed from units and charge/mass numbers.

• epsilon (float, optional) – Normalized cyclotron period: 1 / (cyclotron frequency × time unit). If None, computed from units and charge/mass numbers.

• with_B0 (bool) – Whether to include the effect of a background magnetic field B0 (default: True)

classmethod doc_pde()

PDEs solved by model:

Vlasov equation:

\[\frac{\partial f}{\partial t} + \mathbf{v} \cdot \nabla f + \frac{1}{\varepsilon} \left( \mathbf{E} + \mathbf{v} \times \mathbf{B}_0 \right) \cdot \frac{\partial f}{\partial \mathbf{v}} = 0\]

Ampère’s law:

\[-\frac{\partial \mathbf{E}}{\partial t} = \frac{\alpha^2}{\varepsilon} \int_{\mathbb{R}^3} \mathbf{v} f \, \mathrm{d}^3 \mathbf{v}\]

Initial Poisson equation: At \(t=0\), solve weakly for the electric potential \(\phi\):

\[\begin{split}\int_{\Omega} \nabla \psi^{\top} \cdot \nabla \phi \, \mathrm{d} \mathbf{x} = \frac{\alpha^2}{\varepsilon} \int_{\Omega} \int_{\mathbb{R}^3} \psi \, (f - f_0) \, \mathrm{d}^3 \mathbf{v} \, \mathrm{d} \mathbf{x} \qquad \forall \ \psi \in H^1 \\[2mm] \mathbf{E}(t=0) = -\nabla \phi(t=0)\end{split}\]

classmethod doc_normalization()

Velocity and field normalizations:

\[\hat{v} = c, \qquad \hat{E} = \hat{B} \hat{v}, \qquad \hat{\phi} = \hat{E} \hat{x}\]

Dimensionless parameters:

\[\alpha = \frac{\hat{\omega}_\mathrm{p}}{\hat{\omega}_\mathrm{c}}, \qquad \varepsilon = \frac{1}{\hat{\omega}_\mathrm{c} \hat{t}}\]

where

\[\hat{\omega}_\mathrm{p} = \sqrt{\frac{\hat{n} (Ze)^2}{\epsilon_0 (A m_\mathrm{H})}}, \qquad \hat{\omega}_\mathrm{c} = \frac{(Ze) \hat{B}}{(A m_\mathrm{H})}\]

For electrons: \(Z = -1\), \(A = 1/1836\).

classmethod doc_scalar_quantities()

The following scalars are tracked during simulation:

• Electric field energy: \(E_E = \frac{1}{2} \int |\mathbf{E}|^2 \, \mathrm{d}V\)

• Kinetic energy: \(E_f = \frac{\alpha^2}{2N} \sum_p w_p |\mathbf{v}_p|^2\)

• Total energy: \(E_\mathrm{tot} = E_E + E_f\)

classmethod doc_discretization()

Time integration is performed by the following propagators (in sequence):

1. PushEta

2. PushVxB (if with_B0 is True)

3. VlasovAmpereCoupling

classmethod doc_long_description()

This model couples the kinetic Vlasov equation for the particle distribution function with Ampère’s law for the electric field evolution. It includes the effect of a static background magnetic field \(\mathbf{B}_0\) and solves the initial Poisson equation to satisfy Gauss’s law at \(t=0\). The model uses a particle-in-cell (PIC) method with finite element exterior calculus (FEEC) for the electromagnetic fields.

All field variables are perturbations around a reference equilibrium distribution. The model enables studies of kinetic instabilities, particle-wave interactions, and nonlinear plasma physics without assuming a fluid approximation.

classmethod doc_examples()

Create and initialize a Vlasov-Ampère model:

from struphy.models import VlasovAmpereOneSpecies

model = VlasovAmpereOneSpecies()

# Access fields
model.em_fields.e_field
model.em_fields.phi
model.kinetic_ions.var

# Access tracked scalar quantities
model.scalars["electric_energy"]
model.scalars["kinetic_energy"]
model.scalars["total_energy"]

classmethod doc_use_cases()

This model is appropriate for:

• studying kinetic instabilities in collisionless plasmas

• particle-wave interactions and nonlinear plasma physics

• initial value problems with prescribed equilibrium distributions

• benchmark problems for kinetic-field coupling schemes

• verification of PIC methods in simplified geometries

classmethod doc_cannot_be_used_for()

This model is not suitable for:

• multi-species systems (requires extension)

• collision operator effects

• electromagnetic wave propagation (no magnetic field evolution)

• drift-reduced or gyrokinetic approximations (full 6D Vlasov equation)

class struphy.models.LinearVlasovAmpereOneSpecies(base_units: BaseUnits = BaseUnits(x=1.0, B=1.0, n=1.0, kBT=None), charge_number: int = 1, mass_number: float = 1.0, alpha: float | None = None, epsilon: float | None = None, with_B0: bool = True, with_E0: bool = True)

Bases: StruphyModel

Linearized Vlasov-Ampère equations for one kinetic species around a Maxwellian background.

Parameters:

• base_units (BaseUnits) – Base units for normalization (default: BaseUnits())

• charge_number (int) – Charge number (in units of the positive elementary charge) of the species (default: 1)

• mass_number (float) – Mass number (in units of Proton mass) of the species (default: 1.0)

• alpha (float, optional) – Dimensionless parameter: plasma frequency / cyclotron frequency. If None, computed from units and charge/mass numbers.

• epsilon (float, optional) – Normalized cyclotron period: 1 / (cyclotron frequency × time unit). If None, computed from units and charge/mass numbers.

• with_B0 (bool) – Whether to include the effect of a background magnetic field B0 (default: True)

• with_E0 (bool) – Whether to include the effect of a background electric field E0 (default: True)

classmethod doc_pde()

PDEs solved by model:

Linearized Ampère’s law:

\[\frac{\partial \tilde{\mathbf{E}}}{\partial t} = -\frac{\alpha^2}{\varepsilon} \int_{\mathbb{R}^3} \mathbf{v} \tilde{f} \, \textrm{d}^3 \mathbf{v}\]

Linearized Vlasov equation:

\[\frac{\partial \tilde{f}}{\partial t} + \mathbf{v} \cdot \nabla \tilde{f} + \frac{1}{\varepsilon} \left( \mathbf{E}_0 + \mathbf{v} \times \mathbf{B}_0 \right) \cdot \frac{\partial \tilde{f}}{\partial \mathbf{v}} = \frac{1}{v_{\text{th}}^2 \varepsilon} \tilde{\mathbf{E}} \cdot \mathbf{v} f_0\]

where \(Z=-1\) and \(A=1/1836\) for electrons. The background distribution function \(f_0\) is a uniform Maxwellian

\[f_0 = \frac{n_0(\mathbf{x})}{\left( \sqrt{2 \pi} v_{\text{th}} \right)^3} \exp \left( - \frac{|\mathbf{v}|^2}{2 v_{\text{th}}^2} \right)\]

and the background electric field has to verify the following compatibility condition between with background density

\[\nabla_{\mathbf{x}} \ln (n_0(\mathbf{x})) = \frac{1}{v_{\text{th}}^2 \varepsilon} \mathbf{E}_0\]

At initial time the weak Poisson equation is solved once to weakly satisfy Gauss’ law,

\[\begin{split}\int_\Omega \nabla \psi^\top \cdot \nabla \phi \, \textrm{d} \mathbf{x} &= \frac{\alpha^2}{\varepsilon} \int_\Omega \int_{\mathbb{R}^3} \psi \, \tilde{f} \, \text{d}^3 \mathbf{v} \, \textrm{d} \mathbf{x} \qquad \forall \ \psi \in H^1 \\[2mm] \tilde{\mathbf{E}}(t=0) &= -\nabla \phi(t=0)\end{split}\]

Moreover, it is assumed that

\[\int_{\mathbb{R}^3} \mathbf{v} f_0 \, \text{d}^3 \mathbf{v} = 0\]

classmethod doc_normalization()

The light speed defines the velocity scale:

\[\hat v = c,\qquad \hat E = \hat B \hat v,\qquad \hat\phi = \hat E \hat x.\]

The species parameters are \(\alpha=\hat\Omega_p/\hat\Omega_c\) and \(\varepsilon=1/(\hat\Omega_c\hat t)\).

classmethod doc_scalar_quantities()

The following scalars are tracked during simulation:

• Electric field energy: en_E

• Perturbation particle energy: en_w

• Total energy: en_tot

classmethod doc_discretization()

Time integration is performed by the following propagators (in sequence):

1. PushEta

2. PushVinForceField (if with_E0 is True)

3. EfieldWeightsCoupling

4. PushVxB (if with_B0 is True)

classmethod doc_long_description()

LinearVlasovAmpereOneSpecies is the linearized delta-f counterpart of VlasovAmpereOneSpecies. It is intended for weakly perturbed kinetic problems around a prescribed Maxwellian equilibrium and is especially useful for linear instability or damping studies.

classmethod doc_examples()

Create and initialize the linear Vlasov-Ampère model:

from struphy.models import LinearVlasovAmpereOneSpecies

model = LinearVlasovAmpereOneSpecies()
model.em_fields.e_field
model.em_fields.phi
model.kinetic_ions.var

classmethod doc_use_cases()

This model is appropriate for:

• linear delta-f kinetic instabilities

• weakly perturbed Landau-like damping studies with electrostatic fields

• verification of linearized field-particle coupling

classmethod doc_cannot_be_used_for()

This model is not suitable for:

• strongly nonlinear departures from the chosen equilibrium

• multi-species kinetic coupling

• fully electromagnetic magnetic-field evolution

• equilibria that are not compatible with the built-in Maxwellian assumptions

class struphy.models.VlasovMaxwellOneSpecies(base_units: BaseUnits = BaseUnits(x=1.0, B=1.0, n=1.0, kBT=None), charge_number: int = 1, mass_number: float = 1.0, alpha: float | None = None, epsilon: float | None = None, measure_gauss_law: bool = False)

Bases: StruphyModel

Vlasov-Maxwell equations for one kinetic species.

Parameters:

• base_units (BaseUnits) – Base units for normalization (default: BaseUnits())

• charge_number (int) – Charge number (in units of the positive elementary charge) of the species (default: 1)

• mass_number (float) – Mass number (in units of Proton mass) of the species (default: 1.0)

• alpha (float, optional) – Dimensionless parameter: plasma frequency / cyclotron frequency. If None, computed from units and charge/mass numbers.

• epsilon (float, optional) – Normalized cyclotron period: 1 / (cyclotron frequency × time unit). If None, computed from units and charge/mass numbers.

• measure_gauss_law (bool) – Whether to track the Gauss-law error as a scalar quantity (default: False)

classmethod doc_pde()

PDEs solved by model:

Vlasov equation:

\[\frac{\partial f}{\partial t} + \mathbf{v} \cdot \nabla f + \frac{1}{\varepsilon} \left( \mathbf{E} + \mathbf{v} \times \left( \mathbf{B} + \mathbf{B}_0 \right) \right) \cdot \frac{\partial f}{\partial \mathbf{v}} = 0\]

Ampère’s law:

\[-\frac{\partial \mathbf{E}}{\partial t} + \nabla \times \mathbf{B} = \frac{\alpha^2}{\varepsilon} \int_{\mathbb{R}^3} \mathbf{v} f \, \text{d}^3 \mathbf{v}\]

Faraday’s law:

\[\frac{\partial \mathbf{B}}{\partial t} + \nabla \times \mathbf{E} = 0\]

where \(Z=-1\) and \(A=1/1836\) for electrons.

At initial time the weak Poisson equation is solved once to weakly satisfy Gauss’ law,

\[\begin{split}\int_{\Omega} \nabla \psi^{\top} \cdot \nabla \phi \, \textrm{d} \mathbf{x} = \frac{\alpha^2}{\varepsilon} \int_{\Omega} \int_{\mathbb{R}^3} \psi \, (f - f_0) \, \text{d}^3 \mathbf{v} \, \textrm{d} \mathbf{x} \qquad \forall \ \psi \in H^1 \\[2mm] \mathbf{E}(t=0) = -\nabla \phi(t=0)\end{split}\]

Moreover, it is assumed that

\[\nabla \times \mathbf{B}_0 = \frac{\alpha^2}{\varepsilon} \int_{\mathbb{R}^3} \mathbf{v} f_0 \, \text{d}^3 \mathbf{v}\]

where \(\mathbf{B}_0\) is the static equilibrium magnetic field.

classmethod doc_normalization()

The model uses the light speed as reference velocity:

\[\hat v = c,\qquad \hat E = \hat B \hat v,\qquad \hat\phi = \hat E \hat x.\]

The species parameters are \(\alpha=\hat\Omega_p/\hat\Omega_c\) and \(\varepsilon=1/(\hat\Omega_c\hat t)\).

classmethod doc_scalar_quantities()

The following scalars are tracked during simulation:

• Electric field energy: en_E

• Magnetic field energy: en_B

• Particle kinetic energy: en_f

• Total energy: en_tot

• Optional Gauss-law diagnostic: gauss_error

classmethod doc_discretization()

Time integration is performed by the following propagators (in sequence):

1. MaxwellWeakAmpere

2. PushEta

3. PushVxB

4. VlasovAmpereCoupling

classmethod doc_long_description()

VlasovMaxwellOneSpecies is the fully electromagnetic one-species PIC model in Struphy. It evolves particles and fields self-consistently and supports an optional control-variate formulation for the field coupling.

classmethod doc_examples()

Create and initialize a Vlasov-Maxwell model:

from struphy.models import VlasovMaxwellOneSpecies

model = VlasovMaxwellOneSpecies()
model.em_fields.e_field
model.em_fields.b_field
model.kinetic_ions.var

classmethod doc_use_cases()

This model is appropriate for:

• self-consistent electromagnetic kinetic simulations

• one-species PIC benchmarks

• wave-particle interaction studies with evolving magnetic fields

• verification of the full Vlasov-Maxwell splitting

classmethod doc_cannot_be_used_for()

This model is not suitable for:

• multi-species plasma dynamics without extension

• collisional kinetic closures

• reduced electrostatic-only models where magnetic evolution is unnecessary

• linearized delta-f studies that should use the dedicated linear models

class struphy.models.VlasovMaxwellQuasiNeutral(base_units: BaseUnits = BaseUnits(x=1.0, B=1.0, n=1.0, kBT=None), charge_number: int = 1, mass_number: float = 1.0, alpha: float | None = None, epsilon: float | None = None, lam: float | None = None, kappa: float | None = None, measure_gauss_law: bool = False)

Bases: VlasovMaxwellOneSpecies

Vlasov-Maxwell equations for one kinetic species in the cyclotron normalization, with an implicit, energy conserving field-particle coupling that also holds in the quasi-neutral limit (no displacement current).

Same equations as VlasovMaxwellOneSpecies, normalized with the cyclotron frequency (velocity scale "cyclotron"), so that the coefficients \(\lambda^2\) of the displacement current and \(\kappa^2\) of the plasma current appear in Ampère’s law. Ampère’s law, Faraday’s law and the electric force on the particles are solved together in one Crank-Nicolson step (VlasovMaxwellCoupling); lam=0.0 is the quasi-neutral limit.

Parameters:

• base_units (BaseUnits) – Base units for normalization (default: BaseUnits())

• charge_number (int) – Charge number (in units of the positive elementary charge) of the species (default: 1)

• mass_number (float) – Mass number (in units of Proton mass) of the species (default: 1.0)

• alpha (float, optional) – Dimensionless parameter: plasma frequency / cyclotron frequency. If None, computed from units and charge/mass numbers. It does not enter the equations of this model (which use lam and kappa = alpha lam); when overriding lam and kappa, alpha is not adapted.

• epsilon (float, optional) – Normalized cyclotron period: 1 / (cyclotron frequency × time unit). If None, computed from units and charge/mass numbers (equal to 1 in this model’s normalization).

• lam (float, optional) – Cyclotron frequency × length unit / speed of light (coefficient lam^2 of the displacement current; 0 is the quasi-neutral limit). If None, computed from units and charge/mass numbers.

• kappa (float, optional) – Plasma frequency × length unit / speed of light (coefficient kappa^2 of the plasma current). If None, computed from units and charge/mass numbers.

• measure_gauss_law (bool) – Whether to track the Gauss-law error as a scalar quantity (default: False; not available for lam = 0)

classmethod doc_pde()

PDEs solved by model:

Vlasov equation:

\[\frac{\partial f}{\partial t} + \mathbf{v} \cdot \nabla f + \left( \mathbf{E} + \mathbf{v} \times \left( \mathbf{B} + \mathbf{B}_0 \right) \right) \cdot \frac{\partial f}{\partial \mathbf{v}} = 0\]

Ampère’s law:

\[-\lambda^2 \frac{\partial \mathbf{E}}{\partial t} + \nabla \times \mathbf{B} = \kappa^2 \int_{\mathbb{R}^3} \mathbf{v} f \, \text{d}^3 \mathbf{v}\]

Faraday’s law:

\[\frac{\partial \mathbf{B}}{\partial t} + \nabla \times \mathbf{E} = 0\]

For \(\lambda = 0\) (quasi-neutral limit) Ampère’s law becomes \(\nabla \times \mathbf B = \kappa^2 \int \mathbf v f \, \textrm d^3 \mathbf v\), which implies \(\nabla \cdot \mathbf J = 0\); the electric field is then determined by this constraint, not by its own evolution. The coefficients of the coupling terms can be changed by the options lam, c1, c2, c3 of the propagator VlasovMaxwellCoupling.

At initial time the weak Poisson equation (Gauss’ law \(\lambda^2 \nabla \cdot \mathbf E = \kappa^2 \rho\)) is solved once for \(\lambda > 0\),

\[\begin{split}\lambda^2 \int_{\Omega} \nabla \psi^{\top} \cdot \nabla \phi \, \textrm{d} \mathbf{x} = \kappa^2 \int_{\Omega} \int_{\mathbb{R}^3} \psi \, (f - f_0) \, \text{d}^3 \mathbf{v} \, \textrm{d} \mathbf{x} \qquad \forall \ \psi \in H^1 \\[2mm] \mathbf{E}(t=0) = -\nabla \phi(t=0)\end{split}\]

(for \(\lambda = 0\) the initial electric field does not enter the time stepping). Moreover, it is assumed that \(\nabla \times \mathbf{B}_0 = \kappa^2 \int_{\mathbb{R}^3} \mathbf{v} f_0 \, \text{d}^3 \mathbf{v}\), where \(\mathbf{B}_0\) is the static equilibrium magnetic field.

classmethod doc_normalization()

The model uses the cyclotron velocity of the species as reference velocity:

\[\hat v = \hat\Omega_c \hat x,\qquad \hat t = \frac{1}{\hat\Omega_c},\qquad \hat E = \hat B \hat v,\qquad \hat\phi = \hat E \hat x.\]

The species parameters (see Equation parameters) are

\[\lambda = \frac{\hat\Omega_c \hat x}{c}\,,\qquad \kappa = \frac{\hat\Omega_p \hat x}{c} = \alpha \lambda\,,\]

where \(c/\hat\Omega_c\) and \(c/\hat\Omega_p\) are the inertial lengths (\(\lambda\) carries the sign of the charge, like \(\hat\Omega_c\); only \(\lambda^2\) enters). The Vlasov equation and Faraday’s law then have no coefficients, and Ampère’s law \(\frac{1}{c^2} \partial_t \mathbf E = \nabla \times \mathbf B - \mu_0 \mathbf J\) becomes \(\lambda^2 \partial_t \mathbf E = \nabla \times \mathbf B - \kappa^2 \mathbf J\). The quasi-neutral limit is \(\lambda \to 0\) at fixed \(\kappa\). With \(\lambda = \kappa = 1\) the equations are those of VlasovMaxwellOneSpecies with \(\alpha = 1\), i.e. lengths in units of \(c/\hat\Omega_p\) and \(\hat\Omega_c = \hat\Omega_p\).

classmethod doc_scalar_quantities()

The following scalars are tracked during simulation (the parts of the energy conserved by the coupling step, with the coefficients \(\lambda^2, c_1, c_2, c_3\) of VlasovMaxwellCoupling):

• Electric field energy \(\frac{\lambda^2}{2} \int |\mathbf E|^2 \textrm d \mathbf x\): en_E

• Magnetic field energy \(\frac{c_1}{2} \int |\mathbf B|^2 \textrm d \mathbf x\): en_B

• Particle kinetic energy \(\frac{c_2}{2 c_3} \sum_p w_p |\mathbf v_p|^2\): en_f (\(\frac{\kappa^2}{2} \sum_p w_p |\mathbf v_p|^2\) for the default coefficients)

• Total energy: en_tot (conserved by the coupling step up to the solver tolerance)

• Optional Gauss-law diagnostic: gauss_error (NaN for \(\lambda = 0\))

classmethod doc_discretization()

Time integration is performed by the following propagators (in sequence):

1. VlasovMaxwellCoupling

2. PushEta

3. PushVxB

classmethod doc_long_description()

VlasovMaxwellQuasiNeutral solves the Vlasov-Maxwell system in the cyclotron normalization, with one implicit step for Maxwell’s equations and the electric force, which conserves the energy exactly and allows time steps larger than the inverse plasma frequency and the light-wave CFL limit. The coefficient \(\lambda^2\) of the displacement current can be set to zero (lam=0), which gives the quasi-neutral (Darwin-type) limit, in which the electric field is the Lagrange multiplier of a divergence-free current.

classmethod doc_examples()

Create a quasi-neutral Vlasov-Maxwell model (lengths in units of the skin depth):

from struphy.models import VlasovMaxwellQuasiNeutral

model = VlasovMaxwellQuasiNeutral(lam=0.0, kappa=1.0)
model.em_fields.b_field
model.kinetic_ions.var

classmethod doc_use_cases()

This model is appropriate for:

• electromagnetic kinetic simulations with time steps above the plasma period

• the quasi-neutral limit of the Vlasov-Maxwell system (lam=0.0)

• comparison with VlasovMaxwellOneSpecies

classmethod doc_cannot_be_used_for()

This model is not suitable for:

• the quasi-neutral limit with control-variate (delta-f) weights or with vacuum regions

• multi-species plasma dynamics without extension

• collisional kinetic closures

class struphy.models.LinearVlasovMaxwellOneSpecies(base_units: BaseUnits = BaseUnits(x=1.0, B=1.0, n=1.0, kBT=None), charge_number: int = 1, mass_number: float = 1.0, alpha: float | None = None, epsilon: float | None = None, with_B0: bool = True, with_E0: bool = True)

Bases: LinearVlasovAmpereOneSpecies

Linearized Vlasov-Maxwell equations for one kinetic species around a Maxwellian background.

Parameters:

• base_units (BaseUnits) – Base units for normalization (default: BaseUnits())

• charge_number (int) – Charge number (in units of the positive elementary charge) of the species (default: 1)

• mass_number (float) – Mass number (in units of Proton mass) of the species (default: 1.0)

• alpha (float, optional) – Dimensionless parameter: plasma frequency / cyclotron frequency. If None, computed from units and charge/mass numbers.

• epsilon (float, optional) – Normalized cyclotron period: 1 / (cyclotron frequency × time unit). If None, computed from units and charge/mass numbers.

• with_B0 (bool) – Whether to include the effect of a background magnetic field B0 (default: True)

• with_E0 (bool) – Whether to include the effect of a background electric field E0 (default: True)

classmethod doc_pde()

PDEs solved by model:

Linearized Ampère’s law:

\[\frac{\partial \tilde{\mathbf{E}}}{\partial t} = \nabla \times \tilde{\mathbf{B}} - \frac{\alpha^2}{\varepsilon} \int_{\mathbb{R}^3} \mathbf{v} \tilde{f} \, \textrm{d}^3 \mathbf{v}\]

Linearized Faraday’s law:

\[\frac{\partial \tilde{\mathbf{B}}}{\partial t} = -\nabla \times \tilde{\mathbf{E}}\]

Linearized Vlasov equation:

\[\frac{\partial \tilde{f}}{\partial t} + \mathbf{v} \cdot \nabla \tilde{f} + \frac{1}{\varepsilon} \left( \mathbf{E}_0 + \mathbf{v} \times \mathbf{B}_0 \right) \cdot \frac{\partial \tilde{f}}{\partial \mathbf{v}} = \frac{1}{v_{\text{th}}^2 \varepsilon} \tilde{\mathbf{E}} \cdot \mathbf{v} f_0\]

where \(Z=-1\) and \(A=1/1836\) for electrons. The background distribution function \(f_0\) is a uniform Maxwellian

\[f_0 = \frac{n_0(\mathbf{x})}{\left( \sqrt{2 \pi} v_{\text{th}} \right)^3} \exp \left( - \frac{|\mathbf{v}|^2}{2 v_{\text{th}}^2} \right)\]

and the background electric field has to verify the following compatibility condition between with background density

\[\nabla_{\mathbf{x}} \ln (n_0(\mathbf{x})) = \frac{1}{v_{\text{th}}^2 \varepsilon} \mathbf{E}_0\]

At initial time the weak Poisson equation is solved once to weakly satisfy Gauss’ law,

\[\begin{split}\int_\Omega \nabla \psi^\top \cdot \nabla \phi \, \textrm{d} \mathbf{x} &= \frac{\alpha^2}{\varepsilon} \int_\Omega \int_{\mathbb{R}^3} \psi \, \tilde{f} \, \text{d}^3 \mathbf{v} \, \textrm{d} \mathbf{x} \qquad \forall \ \psi \in H^1 \\[2mm] \tilde{\mathbf{E}}(t=0) &= -\nabla \phi(t=0)\end{split}\]

Moreover, it is assumed that

\[\int_{\mathbb{R}^3} \mathbf{v} f_0 \, \text{d}^3 \mathbf{v} = 0\]

classmethod doc_normalization()

The normalization matches the linear Vlasov-Ampère model:

\[\hat v = c,\qquad \hat E = \hat B \hat v,\qquad \hat\phi = \hat E \hat x.\]

The model retains the same \(\alpha\) and \(\varepsilon\) parameters while adding magnetic perturbation dynamics.

classmethod doc_scalar_quantities()

The following scalars are tracked during simulation:

• Electric field energy: en_E

• Magnetic field energy: en_B

• Perturbation particle energy: en_w

• Total energy: en_tot

classmethod doc_discretization()

Time integration is performed by the following propagators (in sequence):

1. PushEta

2. PushVinForceField (if with_E0 is True)

3. EfieldWeightsCoupling

4. PushVxB (if with_B0 is True)

5. MaxwellWeakAmpere

classmethod doc_long_description()

LinearVlasovMaxwellOneSpecies is the electromagnetic extension of the linear delta-f Vlasov model. It is intended for weakly perturbed kinetic electromagnetic waves and instabilities around a prescribed equilibrium.

classmethod doc_examples()

Create and initialize the linear Vlasov-Maxwell model:

from struphy.models import LinearVlasovMaxwellOneSpecies

model = LinearVlasovMaxwellOneSpecies()
model.em_fields.e_field
model.em_fields.b_field
model.kinetic_ions.var

classmethod doc_use_cases()

This model is appropriate for:

• linear electromagnetic kinetic wave studies

• delta-f verification of linear Vlasov-Maxwell coupling

• weakly nonlinear regimes where equilibrium perturbations stay small

classmethod doc_cannot_be_used_for()

This model is not suitable for:

• strongly nonlinear fully kinetic plasma dynamics

• multi-species electromagnetic coupling

• problems without a well-defined linearization background

• collisional kinetic closures

class struphy.models.VlasovPoissonAdiabatic(base_units: BaseUnits = BaseUnits(x=1.0, B=1.0, n=1.0, kBT=1.0), charge_number: int = 1, mass_number: float = 1.0, alpha: float | None = None, epsilon: float | None = None, with_B0: bool = True)

Bases: StruphyModel

Electrostatic full-orbit (6D) model for a single ion species with adiabatic electrons.

The full-orbit counterpart of DriftKineticElectrostaticAdiabatic: the ions are not reduced to guiding centers, so the polarization density is not modelled but comes out of the ion orbits. The electrons enter through the adiabatic (Boltzmann) response in the Poisson equation, optionally with the same flux-surface average as in the drift-kinetic model.

Parameters:

• base_units (BaseUnits) – Base units for normalization (default: BaseUnits(kBT=1.0)).

• charge_number (int) – Charge number (in units of the positive elementary charge) of the ion species (default: 1).

• mass_number (float) – Mass number (in units of the proton mass) of the ion species (default: 1.0).

• alpha (float, optional) – Ratio of plasma to cyclotron frequency \(\hat \omega_p / \hat \Omega_c\). If None, computed from base_units and the charge/mass numbers.

• epsilon (float, optional) – Normalized inverse cyclotron frequency \(1 / (\hat{\Omega}_\mathrm{c} \hat{t})\). If None, computed from base_units and the charge/mass numbers.

• with_B0 (bool) – Whether to include the Lorentz force of the background magnetic field (default: True).

classmethod doc_pde()

PDEs solved by model:

Vlasov equation:

\[\frac{\partial f}{\partial t} + \mathbf{v} \cdot \nabla f + \frac{1}{\varepsilon} \left( -\nabla \phi + \mathbf{v} \times \mathbf{B}_0 \right) \cdot \frac{\partial f}{\partial \mathbf{v}} = 0\]

Poisson equation with adiabatic electrons:

\[-\frac{\varepsilon}{\alpha^2} \Delta \phi + \frac{1}{Z \varepsilon} \frac{n_0}{T_0} \left( \phi - \langle \phi \rangle \right) = \int_{\mathbb{R}^3} f \, \mathrm{d}^3 \mathbf{v} - n_0\]

where \(f(\mathbf{x}, \mathbf{v}, t)\) is the ion distribution function, \(n_0\) the equilibrium ion density, \(T_0\) the electron temperature and \(Z\) the ion charge number. The vacuum term \(-\varepsilon \alpha^{-2} \Delta \phi\) is dropped in the quasi-neutral limit \(\alpha \to \infty\) (option quasi_neutral). The flux-surface average \(\langle \phi \rangle\) (zonal correction, as in DriftKineticElectrostaticAdiabatic) is optional; by default \(\langle \phi \rangle = 0\).

classmethod doc_normalization()

The reference speed is the ion thermal speed and the electrostatic fields are scaled accordingly:

\[\hat v = \sqrt{k_B \hat T / m_i},\qquad \hat E = \hat v \hat B,\qquad \hat\phi = \hat E \hat x,\]

with temperatures in units of \(k_B \hat T\). The dimensionless parameters are

\[\varepsilon = \frac{1}{\hat\Omega_c\hat t} = \frac{\hat \rho}{\hat x}, \qquad \alpha = \frac{\hat{\omega}_\mathrm{p}}{\hat{\Omega}_\mathrm{c}},\]

such that \(\varepsilon / \alpha\) is the normalized Debye length.

classmethod doc_scalar_quantities()

The following scalars are tracked during simulation (their sum is conserved by the model equations):

• Field energy: en_phi \(= \frac{1}{2\alpha^2} \int |\nabla \phi|^2 \, \mathrm d \mathbf x + \frac{1}{2 Z \varepsilon^2} \int \frac{n_0}{T_0} \phi (\phi - \langle \phi \rangle) \, \mathrm d \mathbf x\) (the first term only if quasi_neutral=False)

• Kinetic energy: en_particles \(= \frac{1}{2} \sum_p w_p |\mathbf{v}_p|^2\) (of \(f - f_0\) with the control variate)

• Total energy: en_tot

classmethod doc_discretization()

Time integration is performed by the following propagators (in sequence):

1. PushEta

2. AdiabaticPoissonSolve

3. PushVinForceField

4. PushVxB (if with_B0 is True)

classmethod doc_long_description()

This model is the full-orbit (6D) counterpart of the electrostatic drift-kinetic model DriftKineticElectrostaticAdiabatic. The ions follow full Lorentz-force orbits in the self-consistent electrostatic field and a static background magnetic field. Electrons are not evolved; they enter through a linearized Boltzmann response in the Poisson equation, optionally without response to the flux-surface-averaged potential (zonal correction), as in the drift-kinetic model. The potential is recomputed from the marker positions at each time step.

In the strong-field limit \(\varepsilon \to 0\) with drift ordering, the linear dispersion relation of this model reduces to that of the drift-kinetic model (up to finite-Larmor-radius corrections of order \(k_\perp^2 \rho_i^2\)), which makes it a reference for verifying guiding-center models. The time step must resolve the cyclotron motion, \(\Delta t \lesssim \varepsilon\). In the quasi-neutral limit it must also resolve the ion-acoustic frequency of the smallest cells, \(\Delta t\, c_s \pi / h \lesssim 2\) (there is no Debye cutoff).

classmethod doc_examples()

Create and initialize the full-orbit adiabatic-electron model:

from struphy.models import VlasovPoissonAdiabatic

model = VlasovPoissonAdiabatic()
model.em_fields.phi
model.kinetic_ions.var

# quasi-neutral (default) or with the vacuum term, optionally with the zonal correction
model.propagators.poisson.options = model.propagators.poisson.Options(
    quasi_neutral=False, which_geometry="cylindrical"
)

classmethod doc_use_cases()

This model is appropriate for:

• electrostatic ion dynamics with adiabatic electrons, including cyclotron and finite-Larmor-radius effects

• ion-acoustic waves, ion Bernstein waves, drift waves and slab ITG modes in the full-orbit description

• verification of drift-kinetic and gyrokinetic models in the strong-field limit

classmethod doc_cannot_be_used_for()

This model is not suitable for:

• electromagnetic dynamics with evolving magnetic perturbations

• electron kinetic effects beyond the adiabatic closure

• long simulations in strongly magnetized plasmas, where resolving the cyclotron motion is too expensive

• multi-species kinetic coupling without extending the model

class struphy.models.VlasovPoissonAdiabaticPauli(base_units: BaseUnits = BaseUnits(x=1.0, B=1.0, n=1.0, kBT=1.0), charge_number: int = 1, mass_number: float = 1.0, alpha: float | None = None, epsilon: float | None = None, with_poisson: bool = True)

Bases: VlasovPoissonAdiabatic

Guiding-center ion dynamics with adiabatic electrons, computed from classical Pauli particles.

The ions are full-orbit (6D) markers with a magnetic moment (classical Pauli particles, Xiao and Qin 2021), loaded on the slow manifold without gyration. They follow guiding-center orbits with the full-orbit pushers of VlasovPoissonAdiabatic, with time steps larger than the cyclotron period. With the control variate (delta-f), the polarization density must be added in the Poisson equation (AdiabaticPoissonSolve.Options(polarization=True)); without the Poisson solve (with_poisson=False) the model computes guiding-center orbits in the equilibrium magnetic field.

Parameters:

• base_units (BaseUnits) – Base units for normalization (default: BaseUnits(kBT=1.0)).

• charge_number (int) – Charge number (in units of the positive elementary charge) of the ion species (default: 1).

• mass_number (float) – Mass number (in units of the proton mass) of the ion species (default: 1.0).

• alpha (float, optional) – Ratio of plasma to cyclotron frequency \(\hat \omega_p / \hat \Omega_c\). If None, computed from base_units and the charge/mass numbers.

• epsilon (float, optional) – Normalized inverse cyclotron frequency \(1 / (\hat{\Omega}_\mathrm{c} \hat{t})\). If None, computed from base_units and the charge/mass numbers.

• with_poisson (bool) – Whether to solve the Poisson equation for the self-consistent potential (default: True); if False, the markers move in the equilibrium magnetic field only (orbit studies).

classmethod doc_pde()

PDEs solved by model:

Vlasov equation for classical Pauli particles, with the magnetic moment \(\mu \geq 0\) as a parameter of the distribution function \(f(\mathbf{x}, \mathbf{v}, \mu, t)\):

\[\frac{\partial f}{\partial t} + \mathbf{v} \cdot \nabla f + \left[ \frac{1}{\varepsilon} \left( \mathbf{E}^\dagger + \mathbf{v} \times \mathbf{B}_0 \right) \right] \cdot \frac{\partial f}{\partial \mathbf{v}} = 0\,, \qquad \mathbf{E}^\dagger = -\nabla \left( \phi + \varepsilon \mu |\mathbf{B}_0| \right)\,,\]

with \(f\) concentrated on the slow manifold \(\mathbf v \times \mathbf b_0 \approx 0\). Poisson equation with adiabatic electrons:

\[-\frac{\varepsilon}{\alpha^2} \Delta \phi + \frac{1}{Z \varepsilon} \frac{n_0}{T_0} \left( \phi - \langle \phi \rangle \right) = \int_0^\infty \int_{\mathbb{R}^3} f \, \mathrm{d}^3 \mathbf{v}\, \mathrm{d} \mu - n_0\,,\]

where \(n_0\) is the equilibrium ion density, \(T_0\) the electron temperature and \(Z\) the ion charge number. On the slow manifold the marker positions are guiding-center positions, so the charge density is the guiding-center density (no gyro-average). With full-f weights the polarization density comes out of the marker orbits (compression by the polarization drift). With the control variate it is removed from the deposited density (see the long description) and is added in the quasi-neutrality equation, as for guiding centers (option polarization of the Poisson solve):

\[-\varepsilon \nabla \cdot \left( \frac{n_0}{|\mathbf B_0|^2} \nabla_\perp \phi \right) + \frac{1}{Z \varepsilon} \frac{n_0}{T_0} \left( \phi - \langle \phi \rangle \right) = \int_0^\infty \int_{\mathbb{R}^3} (f - f_0) \, \mathrm{d}^3 \mathbf{v}\, \mathrm{d} \mu\,.\]

As in VlasovPoissonAdiabatic, the vacuum term is dropped in the quasi-neutral limit (option quasi_neutral) and the flux-surface average \(\langle \phi \rangle\) is optional. With with_poisson=False the Poisson equation is not solved and \(\phi = 0\).

classmethod doc_normalization()

As in VlasovPoissonAdiabatic: the reference speed is the ion thermal speed,

\[\hat v = \sqrt{k_B \hat T / m_i},\qquad \hat E = \hat v \hat B,\qquad \hat\phi = \hat E \hat x,\qquad \hat \mu = \frac{m_i \hat v^2}{\hat B}\,,\]

with the dimensionless parameters \(\varepsilon = 1/(\hat\Omega_c\hat t)\) and \(\alpha = \hat{\omega}_\mathrm{p}/\hat{\Omega}_\mathrm{c}\). The normalized magnetic moment is \(\mu = v_\perp^2 / (2 |\mathbf B_0|)\), as in the guiding-center models.

classmethod doc_scalar_quantities()

The following scalars are tracked during simulation (the sum of the energies is conserved by the model equations):

• Field energy (only if with_poisson): en_phi \(= \frac{1}{2\alpha^2} \int |\nabla \phi|^2 \, \mathrm d \mathbf x + \frac{1}{2 Z \varepsilon^2} \int \frac{n_0}{T_0} \phi (\phi - \langle \phi \rangle) \, \mathrm d \mathbf x\) (the first term only if quasi_neutral=False), plus the polarization energy \(\frac{1}{2} \int \frac{n_0}{|\mathbf B_0|^2} |\nabla_\perp \phi|^2 \, \mathrm d \mathbf x\) with polarization=True

• Kinetic energy: en_particles \(= \frac{1}{2} \sum_p w_p |\mathbf{v}_p|^2\)

• Magnetic-moment energy: en_mu \(= \sum_p w_p \mu_p |\mathbf B_0(\mathbf x_p)|\)

• Total energy: en_tot = en_phi + en_particles + en_mu; with the control variate it is conserved only approximately (the E x B kinetic energy of the background is represented by the polarization energy, that of the delta-f weights remains in en_particles)

• Distance from the slow manifold: mu_prime_max \(= \max_p |\mathbf v_p \times \mathbf b_0|^2 / (2 |\mathbf B_0|)\), Eq. (7) of Xiao and Qin; it must stay small compared to the magnetic moments.

classmethod doc_discretization()

Time integration is performed by the following propagators (in sequence):

1. PushEta

2. AdiabaticPoissonSolve

3. PushVinForceField (with magnetic_moment=True)

4. PushVxB

Without the Poisson solve (with_poisson=False), propagator 2 is left out.

classmethod doc_long_description()

Guiding-center dynamics computed with the slow manifold of the classical Pauli particle (J. Xiao and H. Qin, “Slow manifolds of classical Pauli particle enable structure-preserving geometric algorithms for guiding center dynamics”, Comput. Phys. Commun. 265, 107981 (2021)). A classical Pauli particle is a charged particle with an intrinsic magnetic moment \(\mu\), with the Lagrangian \(L = \frac{1}{2}|\dot{\mathbf x}|^2 + \dot{\mathbf x} \cdot \mathbf A - \phi - \mu B\). When its perpendicular velocity vanishes initially, its effective magnetic moment \(\mu' = |\dot{\mathbf x} \times \mathbf b|^2 / (2B)\) is an adiabatic invariant and stays close to zero: the particle moves on a slow manifold, along which its position follows the guiding-center dynamics of a particle with magnetic moment \(\mu\). Unlike the guiding-center Lagrangian, the Lagrangian of the classical Pauli particle is regular, so the full-orbit Boris scheme can be used with time steps larger than the cyclotron period.

The model is the counterpart of VlasovPoissonAdiabatic (full orbits, with finite-Larmor-radius effects) and DriftKineticElectrostaticAdiabatic (guiding centers with the polarization density in the Poisson equation): the ions have no gyration, hence no finite-Larmor-radius effects, and the polarization drift comes out of the orbits. The magnetic moment of each marker is taken from its drawn perpendicular velocity (PauliParticles6D), such that the markers represent the gyro-average of the initial condition, given in the guiding-center position.

Polarization and delta-f. The polarization drift on the slow manifold is compressible, so the phase-space density of the markers is not constant along their orbits. With full-f weights this compression gives the polarization density. The control variate subtracts the background \(f_0\) along the markers, which removes the compression part \(n_0 \nabla \cdot \boldsymbol \xi\) of the polarization density (\(\boldsymbol \xi\) the polarization displacement) from the deposited charge; the advection part \(\boldsymbol \xi \cdot \nabla n_0\) remains in the weights. Without correction, a slab drift wave then has the frequency without polarization and becomes unstable at grid scale. With the guiding-center polarization term in the Poisson equation (polarization=True, the symmetric operator \(-\nabla \cdot (n_0 |\mathbf B_0|^{-2} \nabla_\perp \phi)\) of the drift-kinetic model, which also contains the advection part), the delta-f slab drift wave has the drift-kinetic frequency, also with time steps of 8 and 16 cyclotron periods. A weak instability remains, which needs the density gradient and is independent of the time step: in the slab example \(\mu'\) stays below \(5 \cdot 10^{-3}\) up to t = 25 and grows by a factor of about 5 per 5 time units and the field energy grows from t = 25 on (about one wave period). A likely cause is that the compression part of the polarization is now instantaneous (field equation), while its advection part follows the marker orbits, including their cyclotron response. Removing only the advection part from the polarization term does not change it. Full-f simulations of small perturbations are dominated by the noise, since the perpendicular ion response is that of cold ions.

classmethod doc_examples()

Create and initialize the model:

from struphy.models import VlasovPoissonAdiabaticPauli

model = VlasovPoissonAdiabaticPauli(epsilon=0.1)
model.em_fields.phi
model.kinetic_ions.var

# with the control variate (delta-f): polarization density in the quasi-neutrality equation
model.propagators.poisson.options = model.propagators.poisson.Options(polarization=True)

# guiding-center orbits in the equilibrium magnetic field, without the Poisson solve
model = VlasovPoissonAdiabaticPauli(epsilon=0.01, with_poisson=False)

# Boris rotation (Cayley transform) for time steps larger than the cyclotron period
model.propagators.push_vxb.options = model.propagators.push_vxb.Options(algo="implicit")

classmethod doc_use_cases()

This model is appropriate for:

• electrostatic guiding-center ion dynamics with adiabatic electrons (delta-f with polarization=True), with time steps larger than the cyclotron period (drift waves, slab ITG modes, …)

• guiding-center orbits (banana orbits, grad-B and curvature drifts) in general magnetic fields, with a full-orbit pusher and time steps larger than the cyclotron period (with_poisson=False)

• verification of guiding-center pushers

classmethod doc_cannot_be_used_for()

This model is not suitable for:

• small perturbations with full-f weights, which are dominated by the noise (use the control variate)

• finite-Larmor-radius effects (use VlasovPoissonAdiabatic)

• magnetic fields without scale separation, where the magnetic moment is not an adiabatic invariant

• time-dependent magnetic fields and electromagnetic dynamics

Gyrokinetic

class struphy.models.DriftKineticElectrostaticAdiabatic(base_units: BaseUnits = BaseUnits(x=1.0, B=1.0, n=1.0, kBT=1.0), charge_number: int = 1, mass_number: float = 1.0, epsilon: float | None = None, alpha: float | None = None, use_diagnostic_poisson: bool = False)

Bases: StruphyModel

Electrostatic drift-kinetic model for a single ion species with adiabatic electrons.

Evolves the guiding-center distribution \(f(\mathbf{X}, v_\parallel, \mu, t)\) under both the \(\mathbf{E}^* \times \mathbf{B}\) drift and parallel streaming along the background magnetic field. Electrons are not evolved kinetically; instead they contribute through an adiabatic response term in the quasi-neutrality (Poisson) equation. The optional control-variate method can be activated for the field solve.

Parameters:

• base_units (BaseUnits) – Base units for normalization (default: BaseUnits(kBT=1.0)).

• charge_number (int) – Charge number (in units of the positive elementary charge) of the ion species (default: 1).

• mass_number (float) – Mass number (in units of the proton mass) of the ion species (default: 1.0).

• epsilon (float, optional) – Normalized inverse cyclotron frequency: \(1 / (\hat{\Omega}_\mathrm{c} \hat{t})\). If None, computed from base_units and the charge/mass numbers.

classmethod doc_pde()

PDEs solved by model:

Drift-kinetic equation:

\[\frac{\partial f}{\partial t} + \left[ v_\parallel \frac{\mathbf{B}^*}{B^*_\parallel} + \frac{\mathbf{E}^* \times \mathbf{b}_0}{B^*_\parallel} \right] \cdot \frac{\partial f}{\partial \mathbf{X}} + \left[ \frac{1}{\varepsilon} \frac{\mathbf{B}^*}{B^*_\parallel} \cdot \mathbf{E}^* \right] \cdot \frac{\partial f}{\partial v_\parallel} = 0\]

Poisson equation:

\[-\nabla_\perp \cdot \left( \frac{n_0}{|B_0|^2} \nabla_\perp \phi \right) + \frac{1}{\varepsilon} n_0 \left( 1 + \frac{1}{Z \varepsilon} \frac{1}{T_0} \phi \right) = \frac{1}{\varepsilon} \int f B^*_\parallel \, \textnormal{d} v_\parallel \textnormal{d} \mu\]

where \(f(\mathbf{X}, v_\parallel, \mu, t)\) is the guiding center distribution and

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

The control variate method can be activated in the Poisson equation; if enabled, the following Poisson equation is solved:

\[-\nabla_\perp \cdot \left( \frac{n_0}{|B_0|^2} \nabla_\perp \phi \right) + \frac{1}{Z \varepsilon^2} \frac{n_0}{T_0} \phi = \frac{1}{\varepsilon} \int (f - f_0) B^*_\parallel \, \textnormal{d} v_\parallel \textnormal{d} \mu\]

classmethod doc_normalization()

The reference speed is the ion thermal speed and the electrostatic fields are scaled accordingly:

\[\hat v = \hat v_i,\qquad \hat E = \hat v_i \hat B,\qquad \hat\phi = \hat E \hat x.\]

The small parameter is \(\varepsilon = 1/(\hat\Omega_c\hat t)\).

classmethod doc_scalar_quantities()

The following scalars are tracked during simulation:

• Field energy: en_phi

• Guiding-center particle energy: en_particles

• Total energy: en_tot

classmethod doc_discretization()

Time integration is performed by the following propagators (in sequence):

1. PushGuidingCenterBxEstar

2. GyrokineticPoissonSolve

3. PushGuidingCenterParallel

classmethod doc_long_description()

This model is an electrostatic drift-kinetic reduction for strongly magnetized ions in a fixed magnetic equilibrium. Electrons are not evolved kinetically; instead they enter through the adiabatic response in the quasi-neutrality solve. The implementation supports control variates for the field solve.

classmethod doc_examples()

Create and initialize the drift-kinetic adiabatic-electron model:

from struphy.models import DriftKineticElectrostaticAdiabatic

model = DriftKineticElectrostaticAdiabatic()
model.em_fields.phi
model.kinetic_ions.var

classmethod doc_use_cases()

This model is appropriate for:

• electrostatic drift-kinetic ion turbulence studies

• strongly magnetized plasmas with adiabatic electrons

• guiding-center PIC verification in realistic magnetic geometry

• low-frequency regimes where full gyrophase resolution is unnecessary

classmethod doc_cannot_be_used_for()

This model is not suitable for:

• fully electromagnetic dynamics with evolving magnetic perturbations

• electron kinetic effects beyond the adiabatic closure

• problems that require resolving full cyclotron motion

• multi-species kinetic coupling without extending the model

class struphy.models.DriftKineticMaxwell(base_units: BaseUnits = BaseUnits(x=1.0, B=1.0, n=1.0, kBT=1.0), charge_number: int = 1, mass_number: float = 1.0, epsilon: float | None = None, lam: float | None = None, kappa: float | None = None)

Bases: StruphyModel

Electromagnetic drift-kinetic equation for one guiding-center species coupled to Maxwell’s equations, in the thermal normalization, with an implicit, energy conserving field-particle coupling that also holds in the quasi-neutral limit.

The electric field enters the guiding-center equations through the gauge-invariant \(\mathbf E^* = \mathbf E - \varepsilon \mu \nabla |B|\). Ampère’s law contains the guiding-center current and the (linearized) polarization current; the displacement current has the coefficient \(\lambda^2 \varepsilon^2 = \hat v_\textrm{th}^2/c^2\), option lam of the coupling propagator DriftKineticMaxwellCoupling; lam=0.0 is the quasi-neutral limit.

Parameters:

• base_units (BaseUnits) – Base units for normalization (default: BaseUnits(kBT=1.0); the thermal energy defines the velocity unit).

• charge_number (int) – Charge number (in units of the positive elementary charge) of the species (default: 1)

• mass_number (float) – Mass number (in units of Proton mass) of the species (default: 1.0)

• epsilon (float, optional) – Normalized cyclotron period: 1 / (cyclotron frequency × time unit). If None, computed from units and charge/mass numbers.

• lam (float, optional) – Cyclotron frequency × length unit / speed of light. If None, computed from units and charge/mass numbers.

• kappa (float, optional) – Plasma frequency × length unit / speed of light. If None, computed from units and charge/mass numbers.

classmethod doc_pde()

PDEs solved by model:

Drift-kinetic equation for \(f(\mathbf X, v_\parallel, \mu, t)\):

\[\frac{\partial f}{\partial t} + \frac{v_\parallel \mathbf B^* + \mathbf E^* \times \mathbf b_0}{B^*_\parallel} \cdot \nabla f + \frac{1}{\varepsilon} \frac{\mathbf B^* \cdot \mathbf E^*}{B^*_\parallel} \frac{\partial f}{\partial v_\parallel} = 0\,,\]

with

\[\mathbf E^* = \mathbf E - \varepsilon \mu \nabla |\mathbf B|\,,\qquad \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^*\,,\]

where \(|\mathbf B| = \mathbf b_0 \cdot (\mathbf B_0 + \tilde{\mathbf B})\): the total magnetic field in the first term of \(\mathbf B^*\) (parallel streaming along the perturbed field lines), the curvature term with the equilibrium unit vector \(\mathbf b_0\).

Ampère’s law, with the displacement current and the (linearized) polarization current of the guiding centers:

\[\varepsilon^2 \lambda^2 \frac{\partial \mathbf E}{\partial t} + \varepsilon^2 \kappa^2 \frac{n_0}{|\mathbf B_0|^2} \frac{\partial \mathbf E_\perp}{\partial t} = \nabla \times \tilde{\mathbf B} - \varepsilon \kappa^2 \iint f B^*_\parallel \left( v_\parallel \frac{\mathbf B^*}{B^*_\parallel} + \frac{\mathbf E \times \mathbf b_0}{B^*_\parallel} \right) \textrm d v_\parallel \textrm d \mu\,,\]

Faraday’s law:

\[\frac{\partial \tilde{\mathbf B}}{\partial t} + \nabla \times \mathbf E = 0\,.\]

\(\mathbf B_0\) is a static equilibrium field, \(n_0\) the equilibrium density, and \(\mathbf E_\perp = \mathbf E - (\mathbf b_0 \cdot \mathbf E) \mathbf b_0\). The current in Ampère’s law is the guiding-center current without the magnetic (grad-\(B\)) drift, whose work on \(\mathbf E\) is therefore missing where \(\nabla |B| \neq 0\). The coefficients \(\varepsilon\lambda\) (option lam, 0 is the quasi-neutral limit), \((\varepsilon\kappa)^2\) (option chi) and \(\varepsilon\kappa^2\) (option c2) can be changed in DriftKineticMaxwellCoupling. Without the polarization current the quasi-neutral limit would be singular for flute modes (\(k_\parallel = 0\)).

classmethod doc_normalization()

The model uses the thermal velocity of the species as reference velocity:

\[\hat v = \hat v_\textrm{th} = \sqrt{\frac{k_\textrm{B} \hat T}{m}}\,,\qquad \hat t = \frac{\hat x}{\hat v}\,, \qquad \hat E = \hat v \hat B\,,\qquad \hat \mu = \frac{m \hat v^2}{\hat B}\,.\]

The species parameters (see Equation parameters) are

\[\varepsilon = \frac{1}{\hat\Omega_c \hat t} = \frac{\hat v}{\hat\Omega_c \hat x}\,,\qquad \lambda = \frac{\hat\Omega_c \hat x}{c}\,,\qquad \kappa = \frac{\hat\Omega_p \hat x}{c}\,.\]

The drift-kinetic equation and Faraday’s law then have only the coefficient \(\varepsilon\) (\(\frac{m \hat v}{q \hat B \hat x} = \varepsilon\) in \(\mathbf B^*\) and \(\mathbf E^*\), \(\frac{q \hat B \hat t}{m} = \frac{1}{\varepsilon}\) in the parallel force). Ampère’s law \(\frac{1}{c^2} \partial_t \mathbf E + \mu_0 \frac{n m}{|B|^2} \partial_t \mathbf E_\perp = \nabla \times \mathbf B - \mu_0 q \mathbf J\), multiplied by \(\hat x / \hat B\), has the coefficients

\[\frac{\hat v^2}{c^2} = \varepsilon^2 \lambda^2\,,\qquad \frac{\mu_0 \hat n m \hat v^2}{\hat B^2} = \frac{\hat v^2}{\hat v_A^2} = \varepsilon^2 \kappa^2\,,\qquad \frac{\mu_0 q \hat n \hat v \hat x}{\hat B} = \varepsilon \kappa^2\,.\]

Their ratio \(\kappa^2/\lambda^2 = \hat\Omega_p^2/\hat\Omega_c^2 = c^2/\hat v_A^2\) compares the polarization with the vacuum dielectric, and \(1/\kappa\) is the inertial length \(c/\hat\Omega_p\) in units of \(\hat x\). For example, deuterium with \(\hat x = 1\) m, \(\hat B = 5\) T, \(\hat n = 10^{20}\) m^-3, \(k_\textrm{B} \hat T = 10\) keV gives \(\varepsilon = 2.9 \cdot 10^{-3}\), \(\lambda = 0.80\), \(\kappa = 31\), i.e. \(\varepsilon^2 \lambda^2 = 5.3 \cdot 10^{-6}\), \(\varepsilon^2 \kappa^2 = 8.0 \cdot 10^{-3}\) and \(\varepsilon \kappa^2 = 2.8\).

classmethod doc_scalar_quantities()

The following scalars are tracked during simulation (the parts of the energy conserved by the coupling step, with the coefficients of DriftKineticMaxwellCoupling, here the defaults):

• Electric field energy (displacement and polarization) \(\frac{\varepsilon^2 \lambda^2}{2} \| \mathbf E \|^2 + \frac{\varepsilon^2 \kappa^2}{2} \| \mathbf E_\perp \|^2_{n_0/|B_0|^2}\): en_E

• Magnetic field energy \(\frac{1}{2} \| \tilde{\mathbf B} \|^2\): en_B

• Particle energy \(\varepsilon^2 \kappa^2 \sum_p w_p (v_{\parallel,p}^2/2 + \mu_p |B_0|)\): en_f

• Total energy: en_tot (conserved by the coupling step up to the solver tolerance)

classmethod doc_discretization()

Time integration is performed by the following propagators (in sequence):

1. DriftKineticMaxwellCoupling

2. PushGuidingCenterParallel

3. PushGuidingCenterBxEstar

classmethod doc_long_description()

DriftKineticMaxwell couples one drift-kinetic (guiding-center) species to Maxwell’s equations in the gauge-invariant \(v_\parallel\) formulation (the electric field, not a potential, enters \(\mathbf E^*\)), normalized with the thermal velocity of the species. Maxwell’s equations and the parallel acceleration are solved together in one implicit, energy conserving step, which allows time steps above the plasma period and has no cancellation problem. Since the drift-kinetic ordering removes the polarization drift, its (linearized) polarization current is added to Ampère’s law; with it, the quasi-neutral limit without displacement current is well-posed and gives the inertial Alfvén wave. The guiding centers then stream along the field and drift with \(\mathbf E^* \times \mathbf b_0\). The total energy is conserved exactly by the splitting if \(\nabla |B_0| = 0\) and either \(\mu = 0\) or \(\mathbf b_0 \cdot \tilde{\mathbf B} = 0\).

classmethod doc_examples()

Create a drift-kinetic Maxwell model in the quasi-neutral limit:

from struphy.models import DriftKineticMaxwell

model = DriftKineticMaxwell()
prop = model.propagators.coupling
prop.options = prop.Options(lam=0.0)

classmethod doc_use_cases()

This model is appropriate for:

• electromagnetic drift-kinetic simulations with time steps above the plasma period

• shear (inertial) Alfvén waves with kinetic parallel dynamics

• the quasi-neutral limit (lam=0.0) with full-f weights

classmethod doc_cannot_be_used_for()

This model is not suitable for:

• polar splines (not supported yet by the coupling propagator)

• the quasi-neutral limit with control-variate (delta-f) weights or with vacuum regions

• finite Larmor radius effects beyond the linearized long-wavelength polarization

• exact energy conservation in inhomogeneous \(|B_0|\) (the grad-\(B\) drift current is missing)

class struphy.models.ToyDrift(base_units: BaseUnits = BaseUnits(x=1.0, B=1.0, n=1.0, kBT=1.0), charge_number: int = 1, mass_number: float = 1.0, epsilon: float | None = None, alpha: float | None = None)

Bases: StruphyModel

Electrostatic drift toy model for a single ion species in a given background magnetic field.

Parameters:

• base_units (BaseUnits) – Base units for normalization (default: BaseUnits(kBT=1.0)).

• charge_number (int) – Charge number (in units of the positive elementary charge) of the ion species (default: 1).

• mass_number (float) – Mass number (in units of the proton mass) of the ion species (default: 1.0).

• epsilon (float, optional) – Normalized cyclotron period: \(1 / (\hat{\omega}_\mathrm{c} \hat{t})\). If None, computed from base_units and the charge/mass numbers.

classmethod doc_pde()

PDEs solved by model:

Drift equation:

\[\frac{\partial f}{\partial t} + \frac{\mathbf{E} \times \mathbf{b}_0}{B^{*}_{\parallel}} \cdot \frac{\partial f}{\partial \mathbf{X}} = 0\]

Poisson equation:

\[-\nabla \cdot \nabla \phi = \int f B^*_\parallel \, \textnormal{d} v_\parallel \textnormal{d} \mu\]

where \(f(\mathbf{X}, v_\parallel, \mu, t)\) is the guiding center distribution and

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

The control variate method can be activated in the Poisson equation; if enabled, the following Poisson equation is solved:

\[-\nabla \cdot \nabla \phi = \int (f - f_0) B^*_\parallel \, \textnormal{d} v_\parallel \textnormal{d} \mu\]

classmethod doc_normalization()

The reference speed is the ion thermal velocity:

\[\hat v = \hat v_i,\qquad \hat E = \hat v_i \hat B,\qquad \hat\phi = \hat E \hat x.\]

classmethod doc_scalar_quantities()

The following scalars are tracked during simulation:

• Electrostatic field energy: en_phi

• Particle kinetic energy: en_particles

• Total energy: en_tot

classmethod doc_discretization()

Time integration is performed by the following propagators (in sequence):

1. PoissonSolve

2. PushGuidingCenterBxEstar

classmethod doc_long_description()

ToyDrift is a stripped-down guiding-center model used to isolate the electrostatic drift part of the dynamics. It is intended for algorithm prototyping and reduced verification problems rather than production drift-kinetic studies.

classmethod doc_examples()

Create and initialize the toy drift model:

from struphy.models import ToyDrift

model = ToyDrift()
model.em_fields.phi
model.kinetic_ions.var

classmethod doc_use_cases()

This model is appropriate for:

• reduced electrostatic guiding-center benchmarks

• testing the field solve plus \(\mathbf E\times\mathbf B\) pusher

• algorithm prototyping before moving to the full drift-kinetic model

classmethod doc_cannot_be_used_for()

This model is not suitable for:

• full drift-kinetic dynamics with parallel streaming

• electromagnetic perturbations

• high-fidelity turbulence studies

• full-orbit kinetic physics

MHD

class struphy.models.LinearMHD(base_units: BaseUnits = BaseUnits(x=1.0, B=1.0, n=1.0, kBT=None), mass_number: float = 1.0)

Bases: StruphyModel

Linear ideal MHD with zero-flow equilibrium for magnetohydrodynamic wave propagation.

Parameters:

• base_units (BaseUnits) – Base units for normalization (default: BaseUnits())

• mass_number (float) – Mass number (in units of Proton mass) of the ion species (default: 1.0)

classmethod doc_pde()

PDEs solved by model:

Continuity (mass conservation):

\[\frac{\partial \tilde{\rho}}{\partial t} + \nabla \cdot (\rho_0 \tilde{\mathbf{U}}) = 0\]

Momentum (Lorentz force):

\[\rho_0 \frac{\partial \tilde{\mathbf{U}}}{\partial t} + \nabla \tilde{p} = (\nabla \times \tilde{\mathbf{B}}) \times \mathbf{B}_0 + (\nabla \times \mathbf{B}_0) \times \tilde{\mathbf{B}}\]

Energy (adiabatic process):

\[\frac{\partial \tilde{p}}{\partial t} + \nabla \cdot (p_0 \tilde{\mathbf{U}}) + \frac{2}{3} p_0 \nabla \cdot \tilde{\mathbf{U}} = 0\]

Induction (Faraday’s law):

\[\frac{\partial \tilde{\mathbf{B}}}{\partial t} - \nabla \times (\tilde{\mathbf{U}} \times \mathbf{B}_0) = 0.\]

classmethod doc_normalization()

All velocities are normalized by the Alfvén velocity:

\[\hat{U} = \hat{v}_\mathrm{A} = \frac{\hat{B}_0}{\sqrt{\mu_0 \rho_0}}\]

The model therefore uses the Alfvén velocity as its characteristic speed scale.

classmethod doc_scalar_quantities()

The following scalars are tracked during simulation:

• Kinetic energy (perturbation): \(E_U = \frac{1}{2} \int \rho_0 |\tilde{\mathbf{U}}|^2 \, \mathrm{d}V\)

• Magnetic energy (perturbation): \(E_B = \frac{1}{2} \int \frac{|\tilde{\mathbf{B}}|^2}{\mu_0} \, \mathrm{d}V\)

• Mean pressure perturbation: \(E_p = \int \frac{\tilde{p}}{\gamma - 1} \, \mathrm{d}V\) (en_p). This is a linear diagnostic, not an energy.

• Thermal perturbation energy: \(E_{\mathrm{thermal}} = \frac{1}{2}\int \frac{\tilde{p}^2}{\gamma p_0} \, \mathrm{d}V\) (en_thermal), with \(\gamma = 5/3\)

• Total perturbed energy: \(E_{\mathrm{tot}} = E_U + E_B + E_{\mathrm{thermal}}\) (en_tot)

• Equilibrium magnetic energy: \(E_{B0} = \frac{1}{2} \int \frac{|\mathbf{B}_0|^2}{\mu_0} \, \mathrm{d}V\)

• Equilibrium internal energy: \(E_{p0} = \int \frac{p_0}{\gamma - 1} \, \mathrm{d}V\)

• Total magnetic energy: \(E_{B,\mathrm{tot}} = \frac{1}{2} \int \frac{|\mathbf{B}_0 + \tilde{\mathbf{B}}|^2}{\mu_0} \, \mathrm{d}V\)

classmethod doc_discretization()

Time integration is performed by the following propagators (in sequence):

1. ShearAlfvenPropagator

2. Magnetosonic

classmethod doc_long_description()

This model simulates small-amplitude perturbations in a magnetized plasma with a static equilibrium magnetic field \(\mathbf{B}_0\) and zero background flow. The model solves the linearized ideal magnetohydrodynamic equations, which couple fluid dynamics (density, velocity, pressure) with magnetic field evolution.

The linear MHD system supports three wave families:

• Alfvén waves: incompressible shear waves propagating along \(\mathbf{B}_0\)

• Fast magnetosonic waves: compressible waves with phase velocity above the Alfvén speed

• Slow magnetosonic waves: compressible waves with phase velocity below the Alfvén speed

All evolved quantities are perturbations around a stationary equilibrium:

• \(\tilde{\rho}\) - density perturbation

• \(\tilde{\mathbf{U}}\) - velocity perturbation

• \(\tilde{p}\) - pressure perturbation

• \(\tilde{\mathbf{B}}\) - magnetic field perturbation

The corresponding equilibrium quantities are:

• \(\rho_0\) - background density

• \(p_0\) - background pressure

• \(\mathbf{B}_0\) - background magnetic field

classmethod doc_examples()

Create and initialize a linear MHD model:

from struphy.models import LinearMHD

model = LinearMHD()

# Access fields
model.em_fields.b_field
model.mhd.density
model.mhd.velocity
model.mhd.pressure

# Access tracked scalar quantities
model.scalars["en_U"]
model.scalars["en_B"]
model.scalars["en_p"]

classmethod doc_use_cases()

This model is appropriate for studying linear wave propagation in a static magnetized plasma, including Alfvén and magnetosonic dynamics around a prescribed equilibrium.

Typical use cases include:

• verification of linear MHD wave dispersion and mode structure

• perturbative studies around static equilibria

• benchmark problems for ideal-MHD field-fluid coupling

classmethod doc_cannot_be_used_for()

This model is not suitable for:

• nonlinear MHD dynamics with finite-amplitude perturbations

• equilibria with non-zero background flow

• dissipative effects such as resistivity or viscosity

• kinetic or particle-field effects beyond the fluid MHD description

class struphy.models.LinearMHDCurlCurl(base_units: BaseUnits = BaseUnits(x=1.0, B=1.0, n=1.0, kBT=None), mass_number: float = 1.0)

Bases: LinearMHD

Linear ideal MHD with zero-flow equilibrium, with the magnetic field in \(H(\textnormal{curl})\).

Same equations as LinearMHD, but with a weak induction equation, such that the implicit shear-Alfvén step leads to a curl-curl system for the magnetic field, which is preconditioned by the Hiptmair-Xu preconditioner (as in ShearAlfvenCurlCurl).

Parameters:

• base_units (BaseUnits) – Base units for normalization (default: BaseUnits())

• mass_number (float) – Mass number (in units of Proton mass) of the ion species (default: 1.0)

classmethod doc_pde()

PDEs solved by model:

Continuity (mass conservation):

\[\frac{\partial \tilde{\rho}}{\partial t} + \nabla \cdot (\rho_0 \tilde{\mathbf{U}}) = 0\]

Momentum (Lorentz force):

\[\rho_0 \frac{\partial \tilde{\mathbf{U}}}{\partial t} + \nabla \tilde{p} = (\nabla \times \tilde{\mathbf{B}}) \times \mathbf{B}_0 + (\nabla \times \mathbf{B}_0) \times \tilde{\mathbf{B}}\]

Energy (adiabatic process):

\[\frac{\partial \tilde{p}}{\partial t} + \nabla \cdot (p_0 \tilde{\mathbf{U}}) + \frac{2}{3} p_0 \nabla \cdot \tilde{\mathbf{U}} = 0\]

Induction (Faraday’s law, solved weakly in \(H(\textnormal{curl})\)):

\[\frac{\partial \tilde{\mathbf{B}}}{\partial t} - \nabla \times (\tilde{\mathbf{U}} \times \mathbf{B}_0) = 0.\]

classmethod doc_normalization()

All velocities are normalized by the Alfvén velocity:

\[\hat{U} = \hat{v}_\mathrm{A} = \frac{\hat{B}_0}{\sqrt{\mu_0 \rho_0}}\]

The model therefore uses the Alfvén velocity as its characteristic speed scale.

classmethod doc_scalar_quantities()

The following scalars are tracked during simulation:

• Kinetic energy (perturbation): \(E_U = \frac{1}{2} \int \rho_0 |\tilde{\mathbf{U}}|^2 \, \mathrm{d}V\)

• Magnetic energy (perturbation, magnetic field as a 1-form): \(E_B = \frac{1}{2} \int \frac{|\tilde{\mathbf{B}}|^2}{\mu_0} \, \mathrm{d}V\)

• Mean pressure perturbation: \(E_p = \int \frac{\tilde{p}}{\gamma - 1} \, \mathrm{d}V\) (en_p). This is a linear diagnostic, not an energy.

• Thermal perturbation energy: \(E_{\mathrm{thermal}} = \frac{1}{2}\int \frac{\tilde{p}^2}{\gamma p_0} \, \mathrm{d}V\) (en_thermal), with \(\gamma = 5/3\)

• Total perturbed energy: \(E_{\mathrm{tot}} = E_U + E_B + E_{\mathrm{thermal}}\) (en_tot). The shear-Alfvén step damps it numerically (the damping vanishes with the resolution along the background field; exactly conserving shear-Alfvén step in LinearMHDCurlCurlwithProj).

• Equilibrium magnetic energy: \(E_{B0} = \frac{1}{2} \int \frac{|\mathbf{B}_0|^2}{\mu_0} \, \mathrm{d}V\) (with \(\mathbf B_0\) projected as a 1-form)

• Equilibrium internal energy: \(E_{p0} = \int \frac{p_0}{\gamma - 1} \, \mathrm{d}V\)

• Total magnetic energy: \(E_{B,\mathrm{tot}} = \frac{1}{2} \int \frac{|\mathbf{B}_0 + \tilde{\mathbf{B}}|^2}{\mu_0} \, \mathrm{d}V\)

classmethod doc_discretization()

Time integration is performed by the following propagators (in sequence):

1. ShearAlfvenCurlCurlPropagator

2. Magnetosonic

classmethod doc_long_description()

LinearMHDCurlCurl solves the same equations as LinearMHD, with the perturbed magnetic field in \(H(\textnormal{curl})\) instead of \(H(\textnormal{div})\). The splitting is the same: the shear-Alfvén step of ShearAlfvenCurlCurl (weak induction equation, field-aligned curl-curl system for the magnetic field, Hiptmair-Xu preconditioner by default), followed by the magnetosonic step, in which the force \((\nabla \times \mathbf B_0) \times \tilde{\mathbf B}\) of the 1-form \(\tilde{\mathbf B}\) enters through the mixed mass matrix M21J.

The shear-Alfvén step is not exactly energy conserving: it damps the waves numerically, with a damping that vanishes quickly with the resolution along the background field (see LinearMHDCurlCurlwithProj for the exactly energy conserving version).

Boundary conditions: with essential boundary conditions in the de Rham sequence, the tangential components of the magnetic field (a 1-form) vanish on that boundary, which is unphysical at a wall with \(\mathbf B_0 \cdot \mathbf n = 0\). Free boundary conditions are therefore recommended: they impose \(\tilde{\mathbf U} \cdot \mathbf n = 0\) weakly in the shear-Alfvén step and \(\tilde p = 0\) weakly in the magnetosonic step (natural boundary condition of the pressure term).

All evolved quantities are perturbations around a stationary equilibrium:

• \(\tilde{\rho}\) - density perturbation

• \(\tilde{\mathbf{U}}\) - velocity perturbation (2-form)

• \(\tilde{p}\) - pressure perturbation

• \(\tilde{\mathbf{B}}\) - magnetic field perturbation (1-form)

classmethod doc_examples()

Create and initialize a linear MHD curl-curl model:

from struphy.models import LinearMHDCurlCurl

model = LinearMHDCurlCurl()

# Access fields
model.em_fields.b_field
model.mhd.density
model.mhd.velocity
model.mhd.pressure

# Access tracked scalar quantities
model.scalars["en_U"]
model.scalars["en_B"]
model.scalars["en_tot"]

classmethod doc_use_cases()

This model is appropriate for:

• linear MHD wave propagation with large time steps (fast curl-curl solve of the shear-Alfvén step)

• comparison with LinearMHD (magnetic field in \(H(\textnormal{div})\))

• testing curl-curl (Hiptmair-Xu) preconditioning in linear MHD

classmethod doc_cannot_be_used_for()

This model is not suitable for:

• nonlinear MHD dynamics with finite-amplitude perturbations

• equilibria with non-zero background flow

• dissipative effects such as resistivity or viscosity

• kinetic or particle-field effects beyond the fluid MHD description

• polar splines with the Hiptmair-Xu preconditioner

class struphy.models.LinearMHDCurlCurlwithProj(base_units: BaseUnits = BaseUnits(x=1.0, B=1.0, n=1.0, kBT=None), mass_number: float = 1.0)

Bases: LinearMHDCurlCurl

Linear ideal MHD with zero-flow equilibrium, with the magnetic field in \(H(\textnormal{curl})\) and the exactly energy conserving shear-Alfvén step.

Same as LinearMHDCurlCurl, but the curl-curl term of the implicit shear-Alfvén step contains the \(\rho_0\)-weighted \(L^2\) projection of the Lorentz force onto \(H(\textnormal{div})\) (as in ShearAlfvenCurlCurlwithProj), at the price of an inner mass matrix solve in every application of the Schur operator (and no effective curl-curl preconditioner).

Parameters:

• base_units (BaseUnits) – Base units for normalization (default: BaseUnits())

• mass_number (float) – Mass number (in units of Proton mass) of the ion species (default: 1.0)

classmethod doc_scalar_quantities()

The following scalars are tracked during simulation:

• Kinetic energy (perturbation): \(E_U = \frac{1}{2} \int \rho_0 |\tilde{\mathbf{U}}|^2 \, \mathrm{d}V\)

• Magnetic energy (perturbation, magnetic field as a 1-form): \(E_B = \frac{1}{2} \int \frac{|\tilde{\mathbf{B}}|^2}{\mu_0} \, \mathrm{d}V\)

• Mean pressure perturbation: \(E_p = \int \frac{\tilde{p}}{\gamma - 1} \, \mathrm{d}V\) (en_p). This is a linear diagnostic, not an energy.

• Thermal perturbation energy: \(E_{\mathrm{thermal}} = \frac{1}{2}\int \frac{\tilde{p}^2}{\gamma p_0} \, \mathrm{d}V\) (en_thermal), with \(\gamma = 5/3\)

• Total perturbed energy: \(E_{\mathrm{tot}} = E_U + E_B + E_{\mathrm{thermal}}\) (en_tot), conserved by the shear-Alfvén step up to the solver tolerances

• Equilibrium magnetic energy: \(E_{B0} = \frac{1}{2} \int \frac{|\mathbf{B}_0|^2}{\mu_0} \, \mathrm{d}V\) (with \(\mathbf B_0\) projected as a 1-form)

• Equilibrium internal energy: \(E_{p0} = \int \frac{p_0}{\gamma - 1} \, \mathrm{d}V\)

• Total magnetic energy: \(E_{B,\mathrm{tot}} = \frac{1}{2} \int \frac{|\mathbf{B}_0 + \tilde{\mathbf{B}}|^2}{\mu_0} \, \mathrm{d}V\)

classmethod doc_discretization()

Time integration is performed by the following propagators (in sequence):

1. ShearAlfvenCurlCurlwithProjPropagator

2. Magnetosonic

classmethod doc_long_description()

LinearMHDCurlCurlwithProj solves the same equations as LinearMHDCurlCurl (magnetic field in \(H(\textnormal{curl})\), weak induction equation, momentum equation projected onto \(H(\textnormal{div})\)), with the projection also in the curl-curl term of the implicit shear-Alfvén step, such that this step conserves the perturbed energy exactly. The Schur complement then contains the inverse of the density-weighted mass matrix of 2-forms. The magnetosonic step is that of LinearMHDCurlCurl. Free boundary conditions are recommended for bounded domains (see LinearMHDCurlCurl).

classmethod doc_examples()

Create and initialize the energy conserving linear MHD curl-curl model:

from struphy.models import LinearMHDCurlCurlwithProj

model = LinearMHDCurlCurlwithProj()
model.em_fields.b_field
model.mhd.velocity
model.scalars["en_tot"]

classmethod doc_use_cases()

This model is appropriate for:

• linear MHD wave propagation with exact energy conservation of the shear-Alfvén step

• comparison with LinearMHDCurlCurl and LinearMHD

classmethod doc_cannot_be_used_for()

This model is not suitable for:

• nonlinear MHD dynamics with finite-amplitude perturbations

• equilibria with non-zero background flow

• dissipative effects such as resistivity or viscosity

• kinetic or particle-field effects beyond the fluid MHD description

• large problems where the inner mass matrix solves of the shear-Alfvén step are too expensive (use LinearMHDCurlCurl)

classmethod doc_pde()

PDEs solved by model:

Continuity (mass conservation):

\[\frac{\partial \tilde{\rho}}{\partial t} + \nabla \cdot (\rho_0 \tilde{\mathbf{U}}) = 0\]

Momentum (Lorentz force):

\[\rho_0 \frac{\partial \tilde{\mathbf{U}}}{\partial t} + \nabla \tilde{p} = (\nabla \times \tilde{\mathbf{B}}) \times \mathbf{B}_0 + (\nabla \times \mathbf{B}_0) \times \tilde{\mathbf{B}}\]

Energy (adiabatic process):

\[\frac{\partial \tilde{p}}{\partial t} + \nabla \cdot (p_0 \tilde{\mathbf{U}}) + \frac{2}{3} p_0 \nabla \cdot \tilde{\mathbf{U}} = 0\]

Induction (Faraday’s law, solved weakly in \(H(\textnormal{curl})\)):

\[\frac{\partial \tilde{\mathbf{B}}}{\partial t} - \nabla \times (\tilde{\mathbf{U}} \times \mathbf{B}_0) = 0.\]

classmethod doc_normalization()

All velocities are normalized by the Alfvén velocity:

\[\hat{U} = \hat{v}_\mathrm{A} = \frac{\hat{B}_0}{\sqrt{\mu_0 \rho_0}}\]

The model therefore uses the Alfvén velocity as its characteristic speed scale.

class struphy.models.LinearExtendedMHDuniform(base_units: BaseUnits = BaseUnits(x=1.0, B=1.0, n=1.0, kBT=None), charge_number: int = 1, mass_number: float = 1.0, epsilon: float | None = None)

Bases: StruphyModel

Linear extended MHD with zero-flow equilibrium, for uniform background conditions only.

Parameters:

• base_units (BaseUnits) – Base units for normalization (default: BaseUnits())

• charge_number (int) – Charge number (in units of the positive elementary charge) of the ion species (default: 1)

• mass_number (float) – Mass number (in units of Proton mass) of the ion species (default: 1.0)

• epsilon (float, optional) – Normalized cyclotron period: 1 / (cyclotron frequency × time unit). If None, computed from units and charge/mass numbers.

classmethod doc_pde()

PDEs solved by model:

Continuity:

\[\frac{\partial \tilde{\rho}}{\partial t} + \nabla \cdot (\rho_0 \tilde{\mathbf{U}}) = 0\]

Momentum:

\[\rho_0 \frac{\partial \tilde{\mathbf{U}}}{\partial t} + \nabla \tilde{p} = (\nabla \times \tilde{\mathbf{B}}) \times \mathbf{B}_0\]

Pressure:

\[\frac{\partial \tilde{p}}{\partial t} + \frac{5}{3} p_0 \nabla \cdot \tilde{\mathbf{U}} = 0\]

Induction:

\[\frac{\partial \tilde{\mathbf{B}}}{\partial t} - \nabla \times \left( \tilde{\mathbf{U}} \times \mathbf{B}_0 - \frac{1}{\varepsilon} \frac{\nabla \times \tilde{\mathbf{B}}}{\rho_0} \times \mathbf{B}_0 \right) = 0\]

classmethod doc_normalization()

All velocities are normalized by the Alfvén speed,

\[\hat U = \hat v_A,\]

and Hall effects enter through \(\varepsilon = 1/(\hat\Omega_c\hat t)\).

classmethod doc_scalar_quantities()

The following scalars are tracked during simulation:

• Flow kinetic energy: en_U

• Mean pressure perturbation (linear diagnostic): en_p

• Thermal perturbation energy: en_thermal

• Magnetic perturbation energy: en_B

• Background pressure energy: en_p_eq

• Background magnetic energy: en_B_eq

• Total magnetic energy: en_B_tot

• Total perturbation energy (kinetic + magnetic + thermal): en_tot

• Magnetic helicity-like invariant: helicity

classmethod doc_discretization()

Time integration is performed by the following propagators (in sequence):

1. ShearAlfvenB1

2. Hall

3. MagnetosonicUniform

classmethod doc_long_description()

LinearExtendedMHDuniform extends the linear MHD model with Hall dynamics for a uniform equilibrium. It is targeted at wave propagation and mode studies where the equilibrium is simple but finite Hall corrections are important.

classmethod doc_examples()

Create and initialize a linear extended-MHD model:

from struphy.models import LinearExtendedMHDuniform

model = LinearExtendedMHDuniform()
model.em_fields.b_field
model.mhd.density
model.mhd.velocity
model.mhd.pressure

classmethod doc_use_cases()

This model is appropriate for:

• linear Hall-MHD wave studies

• uniform-equilibrium benchmarks for extended MHD

• comparison against ideal-MHD dispersion in the Hall-corrected regime

classmethod doc_cannot_be_used_for()

This model is not suitable for:

• nonlinear extended-MHD turbulence

• non-uniform equilibrium configurations

• kinetic ion/electron effects beyond the Hall correction

• dissipation-dominated problems with explicit viscosity or resistivity

class struphy.models.ShearAlfven(base_units: BaseUnits = BaseUnits(x=1.0, B=1.0, n=1.0, kBT=None), mass_number: float = 1.0)

Bases: StruphyModel

Shear-Alfvén part of LinearMHD with zero-flow equilibrium.

Parameters:

• base_units (BaseUnits) – Base units for normalization (default: BaseUnits())

• mass_number (float) – Mass number (in units of Proton mass) of the ion species (default: 1.0)

classmethod doc_pde()

PDEs solved by model:

Momentum:

\[\rho_0 \frac{\partial \tilde{\mathbf{U}}}{\partial t} = (\nabla \times \tilde{\mathbf{B}}) \times \mathbf{B}_0\]

Induction:

\[\frac{\partial \tilde{\mathbf{B}}}{\partial t} - \nabla \times (\tilde{\mathbf{U}} \times \mathbf{B}_0) = 0\]

classmethod doc_normalization()

All velocities are normalized by the Alfvén speed,

\[\hat U = \hat v_A.\]

classmethod doc_scalar_quantities()

The following scalars are tracked during simulation:

• Perturbed kinetic energy: en_U

• Perturbed magnetic energy: en_B

• Background magnetic energy: en_B_eq

• Total magnetic energy including background: en_B_tot

• Perturbation total energy: en_tot

• Total energy including background contribution: en_tot2

classmethod doc_discretization()

Time integration is performed by the following propagators (in sequence):

1. ShearAlfvenPropagator

classmethod doc_long_description()

ShearAlfven is the incompressible magnetic-tension subsystem extracted from linear MHD. It is useful when one wants to isolate Alfvénic propagation from the full linear MHD model.

classmethod doc_examples()

Create and initialize a shear-Alfvén model:

from struphy.models import ShearAlfven

model = ShearAlfven()
model.em_fields.b_field
model.mhd.velocity

classmethod doc_use_cases()

This model is appropriate for:

• isolated shear-Alfvén wave propagation

• verification of the ShearAlfven propagator

• reduced linear-MHD studies without compressibility

classmethod doc_cannot_be_used_for()

This model is not suitable for:

• magnetosonic or compressible wave physics

• nonlinear MHD dynamics

• kinetic particle coupling

• resistive, viscous, or Hall-MHD effects

class struphy.models.ShearAlfvenCurlCurl(base_units: BaseUnits = BaseUnits(x=1.0, B=1.0, n=1.0, kBT=None), mass_number: float = 1.0)

Bases: StruphyModel

Shear-Alfvén part of linear MHD with zero-flow equilibrium, with the magnetic field in \(H(\textnormal{curl})\).

Same equations as ShearAlfven, but with a weak induction equation, such that the implicit time step leads to a curl-curl system for the magnetic field, which is preconditioned by the Hiptmair-Xu preconditioner.

Parameters:

• base_units (BaseUnits) – Base units for normalization (default: BaseUnits())

• mass_number (float) – Mass number (in units of Proton mass) of the ion species (default: 1.0)

classmethod doc_pde()

PDEs solved by model:

Momentum:

\[\rho_0 \frac{\partial \tilde{\mathbf{U}}}{\partial t} = (\nabla \times \tilde{\mathbf{B}}) \times \mathbf{B}_0\]

Induction:

\[\frac{\partial \tilde{\mathbf{B}}}{\partial t} - \nabla \times (\tilde{\mathbf{U}} \times \mathbf{B}_0) = 0\]

classmethod doc_normalization()

All velocities are normalized by the Alfvén speed,

\[\hat U = \hat v_A.\]

classmethod doc_scalar_quantities()

The following scalars are tracked during simulation:

• Perturbed kinetic energy: en_U

• Perturbed magnetic energy: en_B

• Background magnetic energy: en_B_eq

• Total magnetic energy including background: en_B_tot

• Perturbation total energy: en_tot (not increased by the implicit time step: numerical damping that vanishes with the resolution along the background field; exactly conserved in ShearAlfvenCurlCurlwithProj)

• Total energy including background contribution: en_tot2

classmethod doc_discretization()

Time integration is performed by the following propagators (in sequence):

1. ShearAlfvenCurlCurlPropagator

classmethod doc_long_description()

ShearAlfvenCurlCurl solves the same equations as ShearAlfven, with the perturbed magnetic field in \(H(\textnormal{curl})\) instead of \(H(\textnormal{div})\). The momentum equation is the \(\rho_0\)-weighted \(L^2\) projection onto \(H(\textnormal{div})\), and the induction equation is weak. In the implicit time step, the velocity update is inserted into the induction equation before the projection, which gives a field-aligned curl-curl system for the magnetic field, preconditioned by the Hiptmair-Xu auxiliary space preconditioner. This step is not exactly energy conserving: it damps the waves numerically, with a damping that vanishes quickly with the resolution along the background field (see ShearAlfvenCurlCurlwithProj for the exactly energy conserving version). With essential boundary conditions in the de Rham sequence, the tangential components of the magnetic field (a 1-form) vanish on that boundary (stronger than \(\mathbf B \cdot \mathbf n = 0\) of ShearAlfven).

classmethod doc_examples()

Create and initialize a shear-Alfvén curl-curl model:

from struphy.models import ShearAlfvenCurlCurl

model = ShearAlfvenCurlCurl()
model.em_fields.b_field
model.mhd.velocity

classmethod doc_use_cases()

This model is appropriate for:

• isolated shear-Alfvén wave propagation with large time steps

• testing curl-curl (Hiptmair-Xu) preconditioning of the shear-Alfvén step

• reduced linear-MHD studies without compressibility

classmethod doc_cannot_be_used_for()

This model is not suitable for:

• magnetosonic or compressible wave physics

• nonlinear MHD dynamics

• kinetic particle coupling

• resistive, viscous, or Hall-MHD effects

• polar splines with the Hiptmair-Xu preconditioner

class struphy.models.ShearAlfvenCurlCurlwithProj(base_units: BaseUnits = BaseUnits(x=1.0, B=1.0, n=1.0, kBT=None), mass_number: float = 1.0)

Bases: ShearAlfvenCurlCurl

Shear-Alfvén part of linear MHD with zero-flow equilibrium, with the magnetic field in \(H(\textnormal{curl})\) and the exactly energy conserving time step.

Same as ShearAlfvenCurlCurl, but the curl-curl term of the implicit time step contains the \(\rho_0\)-weighted \(L^2\) projection of the Lorentz force onto \(H(\textnormal{div})\), which conserves the perturbed energy exactly (up to the solver tolerances), at the price of an inner mass matrix solve in every application of the Schur operator (and no effective curl-curl preconditioner).

Parameters:

• base_units (BaseUnits) – Base units for normalization (default: BaseUnits())

• mass_number (float) – Mass number (in units of Proton mass) of the ion species (default: 1.0)

classmethod doc_scalar_quantities()

The following scalars are tracked during simulation:

• Perturbed kinetic energy: en_U

• Perturbed magnetic energy: en_B

• Background magnetic energy: en_B_eq

• Total magnetic energy including background: en_B_tot

• Perturbation total energy: en_tot (conserved by the implicit time step up to the solver tolerances)

• Total energy including background contribution: en_tot2

classmethod doc_discretization()

Time integration is performed by the following propagators (in sequence):

1. ShearAlfvenCurlCurlwithProjPropagator

classmethod doc_long_description()

ShearAlfvenCurlCurlwithProj solves the same equations as ShearAlfvenCurlCurl (magnetic field in \(H(\textnormal{curl})\), weak induction equation, momentum equation projected onto \(H(\textnormal{div})\)), with the projection also in the curl-curl term of the implicit time step, such that the perturbed energy is conserved exactly. The Schur complement then contains the inverse of the density-weighted mass matrix of 2-forms.

classmethod doc_examples()

Create and initialize the energy conserving shear-Alfvén curl-curl model:

from struphy.models import ShearAlfvenCurlCurlwithProj

model = ShearAlfvenCurlCurlwithProj()
model.em_fields.b_field
model.mhd.velocity

classmethod doc_use_cases()

This model is appropriate for:

• isolated shear-Alfvén wave propagation with exact energy conservation

• comparison with ShearAlfvenCurlCurl

classmethod doc_pde()

PDEs solved by model:

Momentum:

\[\rho_0 \frac{\partial \tilde{\mathbf{U}}}{\partial t} = (\nabla \times \tilde{\mathbf{B}}) \times \mathbf{B}_0\]

Induction:

\[\frac{\partial \tilde{\mathbf{B}}}{\partial t} - \nabla \times (\tilde{\mathbf{U}} \times \mathbf{B}_0) = 0\]

classmethod doc_normalization()

All velocities are normalized by the Alfvén speed,

\[\hat U = \hat v_A.\]

classmethod doc_cannot_be_used_for()

This model is not suitable for:

• magnetosonic or compressible wave physics

• nonlinear MHD dynamics

• kinetic particle coupling

• resistive, viscous, or Hall-MHD effects

• polar splines with the Hiptmair-Xu preconditioner

class struphy.models.VariationalMHD(base_units: BaseUnits = BaseUnits(x=1.0, B=1.0, n=1.0, kBT=None), mass_number: float = 1.0, with_viscosity: bool = True, with_resistivity: bool = True, divdiv_alpha: float = 0.0)

Bases: StruphyModel

Full (non-linear) visco-resistive MHD equations discretized with a variational method.

Parameters:

• base_units (BaseUnits) – Base units for normalization (default: BaseUnits())

• mass_number (float) – Mass number (in units of Proton mass) of the fluid species (default: 1.0)

• with_viscosity (bool) – Whether to include viscous dissipation (default: True)

• with_resistivity (bool) – Whether to include resistive dissipation (default: True)

• divdiv_alpha (float) – Coefficient alpha >= 0 of the div-div regularization alpha int rho (div u)^2 dx of the kinetic energy (Hamiltonian regularization, penalizes compressions; default: 0.0), see kinetic_metric()

classmethod doc_pde()

PDEs solved by model:

Continuity:

\[\partial_t \rho + \nabla \cdot (\rho \mathbf{u}) = 0\]

Momentum:

\[\partial_t (\rho \mathbf{u}) + \nabla \cdot (\rho \mathbf{u} \otimes \mathbf{u}) + \rho \nabla \frac{(\rho \mathcal{U}(\rho, s))}{\partial \rho} + s \nabla \frac{(\rho \mathcal{U}(\rho, s))}{\partial s} + \mathbf{B} \times \nabla \times \mathbf{B} - \nabla \cdot \left( (\mu + \mu_a(\mathbf{x})) \nabla \mathbf{u} \right) = 0\]

Entropy:

\[\partial_t s + \nabla \cdot (s \mathbf{u}) = \frac{1}{T} \left( (\mu + \mu_a(\mathbf{x})) |\nabla \mathbf{u}|^2 + (\eta + \eta_a(\mathbf{x})) |\nabla \times \mathbf{B}|^2 \right)\]

Induction:

\[\partial_t \mathbf{B} + \nabla \times (\mathbf{B} \times \mathbf{u}) + \nabla \times (\eta + \eta_a(\mathbf{x})) \nabla \times \mathbf{B} = 0\]

where the internal energy per unit mass is \(\mathcal U(\rho) = \rho^{\gamma-1} \exp(s / \rho)\), and \(\mu_a(\mathbf{x})\) and \(\eta_a(\mathbf{x})\) are artificial viscosity and resistivity coefficients.

With 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\) (Hamiltonian regularization of compressions). The equations follow from the same variational principle with the momentum \(\rho \mathbf u - 2 \alpha \nabla (\rho \nabla \cdot \mathbf u)\) instead of \(\rho \mathbf u\), and the density derivative \(\alpha (\nabla \cdot \mathbf u)^2\) of the kinetic energy enters the force like \(|\mathbf u|^2/2\).

classmethod doc_normalization()

The model uses Alfvén-speed normalization together with entropy-based thermodynamic units inherited from the chosen internal-energy functional.

classmethod doc_scalar_quantities()

The following scalars are tracked during simulation:

• Kinetic energy: en_U

• Thermodynamic energy: en_thermo

• Magnetic energy: en_mag

• Div-div part of the kinetic energy (only for divdiv_alpha > 0): en_divdiv

• Total energy: en_tot

• Total density / entropy: dens_tot, entr_tot

• Divergence diagnostic: tot_div_B

classmethod doc_discretization()

Time integration is performed by the following propagators (in sequence):

1. VariationalDensityEvolve

2. VariationalMomentumAdvection

3. VariationalEntropyEvolve

4. VariationalMagFieldEvolve

5. VariationalViscosity (if with_viscosity is True)

6. VariationalResistivity (if with_resistivity is True)

classmethod doc_long_description()

VariationalMHD is the full nonlinear dissipative MHD model in the entropy formulation. It is the most complete variational MHD model in this family when both viscosity and resistivity are required.

classmethod doc_examples()

Create and initialize the full visco-resistive MHD model:

from struphy.models import VariationalMHD

model = VariationalMHD()
model.mhd.density
model.mhd.velocity
model.mhd.entropy
model.em_fields.b_field

classmethod doc_use_cases()

This model is appropriate for:

• nonlinear dissipative MHD simulations

• entropy-based variational MHD benchmarks

• studies where both viscosity and resistivity matter

classmethod doc_cannot_be_used_for()

This model is not suitable for:

• ideal nondissipative MHD if strict conservation is required

• kinetic or hybrid particle effects

• reduced linear perturbation studies better served by the linear models

class struphy.models.VariationalMHD_with_p(base_units: BaseUnits = BaseUnits(x=1.0, B=1.0, n=1.0, kBT=None), mass_number: float = 1.0, with_viscosity: bool = True, with_resistivity: bool = True)

Bases: StruphyModel

Full (non-linear) visco-resistive MHD equations, with the pressure variable, discretized with a variational method.

Parameters:

• base_units (BaseUnits) – Base units for normalization (default: BaseUnits())

• mass_number (float) – Mass number (in units of Proton mass) of the fluid species (default: 1.0)

• with_viscosity (bool) – Whether to include viscous dissipation (default: True)

• with_resistivity (bool) – Whether to include resistive dissipation (default: True)

classmethod doc_pde()

PDEs solved by model:

Continuity:

\[\partial_t \rho + \nabla \cdot (\rho \mathbf{u}) = 0\]

Momentum:

\[\partial_t (\rho \mathbf{u}) + \nabla \cdot (\rho \mathbf{u} \otimes \mathbf{u}) + \frac{1}{\gamma - 1} \nabla p + \mathbf{B} \times \nabla \times \mathbf{B} - \nabla \cdot \left( (\mu + \mu_a(\mathbf{x})) \nabla \mathbf{u} \right) = 0\]

Pressure:

\[\partial_t p + u \cdot \nabla p + \gamma p \nabla \cdot u = \frac{1}{\gamma - 1} \left( (\mu + \mu_a(\mathbf{x})) |\nabla \mathbf{u}|^2 + (\eta + \eta_a(\mathbf{x})) |\nabla \times \mathbf{B}|^2 \right)\]

Induction:

\[\partial_t \mathbf{B} + \nabla \times (\mathbf{B} \times \mathbf{u}) + \nabla \times (\eta + \eta_a(\mathbf{x})) \nabla \times \mathbf{B} = 0\]

Here \(\mu_a(\mathbf{x})\) and \(\eta_a(\mathbf{x})\) are artificial viscosity and resistivity coefficients.

classmethod doc_normalization()

The characteristic flow speed is the Alfvén speed,

\[\hat u = \hat v_A.\]

classmethod doc_scalar_quantities()

The following scalars are tracked during simulation:

• Kinetic energy: en_U

• Thermodynamic energy: en_thermo

• Magnetic energy: en_mag

• Total energy: en_tot

• Total density: dens_tot

• Divergence diagnostic: tot_div_B

classmethod doc_discretization()

Time integration is performed by the following propagators (in sequence):

1. VariationalDensityEvolve

2. VariationalMomentumAdvection

3. VariationalPBEvolve

4. VariationalViscosity (if with_viscosity is True)

5. VariationalResistivity (if with_resistivity is True)

classmethod doc_long_description()

This is the pressure-based full nonlinear dissipative MHD model. It is the natural counterpart to the entropy- and q-based full MHD variants when pressure is preferred as primitive thermodynamic variable.

classmethod doc_examples()

Create and initialize the pressure-based visco-resistive MHD model:

from struphy.models import VariationalMHD_with_p

model = VariationalMHD_with_p()
model.mhd.pressure
model.em_fields.b_field

classmethod doc_use_cases()

This model is appropriate for:

• full nonlinear dissipative MHD in p-formulation

• comparison against entropy- and q-based full MHD models

• variational resistive/viscous MHD benchmarks

classmethod doc_cannot_be_used_for()

This model is not suitable for:

• ideal MHD without dissipation

• kinetic plasma physics

• reduced linear perturbation studies

• thermodynamic formulations that require entropy or q as primitive variables

class struphy.models.VariationalMHD_with_q(base_units: BaseUnits = BaseUnits(x=1.0, B=1.0, n=1.0, kBT=None), mass_number: float = 1.0, with_viscosity: bool = True, with_resistivity: bool = True)

Bases: StruphyModel

Full (non-linear) visco-resistive MHD equations, with the q variable (square root of the pressure), discretized with a variational method.

Parameters:

• base_units (BaseUnits) – Base units for normalization (default: BaseUnits())

• mass_number (float) – Mass number (in units of Proton mass) of the fluid species (default: 1.0)

• with_viscosity (bool) – Whether to include viscous dissipation (default: True)

• with_resistivity (bool) – Whether to include resistive dissipation (default: True)

classmethod doc_pde()

PDEs solved by model:

Continuity:

\[\partial_t \rho + \nabla \cdot (\rho \mathbf{u}) = 0\]

Momentum:

\[\partial_t (\rho \mathbf{u}) + \nabla \cdot (\rho \mathbf{u} \otimes \mathbf{u}) + \frac{2 q}{\gamma - 1} \nabla q + \mathbf{B} \times \nabla \times \mathbf{B} - \nabla \cdot \left( (\mu + \mu_a(\mathbf{x})) \nabla \mathbf{u} \right) = 0\]

Energy-like variable:

\[\partial_t q + \nabla \cdot (q \mathbf{u}) + \left( \frac{\gamma}{2} - 1 \right) q \nabla \cdot \mathbf{u} = \frac{2 q}{\gamma - 1} \left( (\mu + \mu_a(\mathbf{x})) |\nabla \mathbf{u}|^2 + (\eta + \eta_a(\mathbf{x})) |\nabla \times \mathbf{B}|^2 \right)\]

Induction:

\[\partial_t \mathbf{B} + \nabla \times (\mathbf{B} \times \mathbf{u}) + \nabla \times (\eta + \eta_a(\mathbf{x})) \nabla \times \mathbf{B} = 0\]

Here \(\mu_a(\mathbf{x})\) and \(\eta_a(\mathbf{x})\) are artificial viscosity and resistivity coefficients.

classmethod doc_normalization()

The flow variables are normalized with the Alfvén speed:

\[\hat u = \hat v_A.\]

classmethod doc_scalar_quantities()

The following scalars are tracked during simulation:

• Kinetic energy: en_U

• Thermodynamic q-energy: en_thermo

• Magnetic energy: en_mag

• Total energy: en_tot

• Total density: dens_tot

• Divergence diagnostic: tot_div_B

classmethod doc_discretization()

Time integration is performed by the following propagators (in sequence):

1. VariationalDensityEvolve

2. VariationalMomentumAdvection

3. VariationalQBEvolve

4. VariationalViscosity (if with_viscosity is True)

5. VariationalResistivity (if with_resistivity is True)

classmethod doc_long_description()

This model is the full nonlinear q-based dissipative MHD formulation. It is useful when the square-root pressure variable is preferred in the variational discretization while still retaining full nonlinear magnetic dynamics.

classmethod doc_examples()

Create and initialize the q-based visco-resistive MHD model:

from struphy.models import VariationalMHD_with_q

model = VariationalMHD_with_q()
model.mhd.sqrt_p
model.em_fields.b_field

classmethod doc_use_cases()

This model is appropriate for:

• full nonlinear dissipative MHD in q-form

• comparison against the p- and entropy-based full models

• variational MHD studies where q is numerically advantageous

classmethod doc_cannot_be_used_for()

This model is not suitable for:

• ideal dissipation-free MHD

• kinetic or hybrid particle physics

• reduced linear perturbation problems

• pressure or entropy primitive-variable studies

class struphy.models.VariationalPerturbativeMHD(base_units: BaseUnits = BaseUnits(x=1.0, B=1.0, n=1.0, kBT=None), mass_number: float = 1.0, linear: bool = False, with_viscosity: bool = True, with_resistivity: bool = True)

Bases: StruphyModel

Perturbative (delta-f) or linear visco-resistive MHD equations, with the pressure variable, discretized with a variational method.

Parameters:

• base_units (BaseUnits) – Base units for normalization (default: BaseUnits())

• mass_number (float) – Mass number (in units of Proton mass) of the fluid species (default: 1.0)

• linear (bool) – Whether to solve the linearized equations, i.e. to drop the terms that are nonlinear in the perturbations (default: False)

• with_viscosity (bool) – Whether to include viscous dissipation (default: True)

• with_resistivity (bool) – Whether to include resistive dissipation (default: True)

classmethod doc_pde()

PDEs solved by model:

The unknowns are the perturbations \((\tilde\rho, \tilde p, \tilde{\mathbf B})\) with respect to a zero-flow MHD equilibrium \((\rho_0, p_0, \mathbf B_0)\), and the velocity \(\mathbf u\).

Continuity:

\[\partial_t \tilde{\rho} + \nabla \cdot \left( (\tilde{\rho} + \rho_0) \tilde{\mathbf{u}} \right) = 0\]

Momentum:

\[\partial_t \left( (\tilde{\rho} + \rho_0) \tilde{\mathbf{u}} \right) + \nabla \cdot \left( (\tilde{\rho} + \rho_0) \tilde{\mathbf{u}} \otimes \tilde{\mathbf{u}} \right) + \nabla \tilde{p} + \mathbf{B}_0 \times \nabla \times \tilde{\mathbf{B}} + \tilde{\mathbf{B}} \times \nabla \times \mathbf{B}_0 + \tilde{\mathbf{B}} \times \nabla \times \tilde{\mathbf{B}} - \nabla \cdot \left( (\mu + \mu_a(\mathbf{x})) \nabla \tilde{\mathbf{u}} \right) = 0\]

Pressure:

\[\partial_t \tilde{p} + \tilde{\mathbf{u}} \cdot \nabla (\tilde{p} + p_0) + \gamma (\tilde{p} + p_0) \nabla \cdot \tilde{\mathbf{u}} = (\gamma - 1) \left( (\mu + \mu_a(\mathbf{x})) |\nabla \tilde{\mathbf{u}}|^2 + (\eta + \eta_a(\mathbf{x})) |\nabla \times \tilde{\mathbf{B}}|^2 \right)\]

Induction:

\[\partial_t \tilde{\mathbf{B}} + \nabla \times \left( (\tilde{\mathbf{B}} + \mathbf{B}_0) \times \tilde{\mathbf{u}} \right) + \nabla \times (\eta + \eta_a(\mathbf{x})) \nabla \times \tilde{\mathbf{B}} = 0\]

Here \(\mu_a(\mathbf{x})\) and \(\eta_a(\mathbf{x})\) are artificial viscosity and resistivity coefficients.

Linear model (linear=True): the terms that are nonlinear in the perturbations are dropped,

\[\begin{split}&\partial_t \tilde{\rho} + \nabla \cdot (\rho_0 \tilde{\mathbf{u}}) = 0\,, \\[2mm] &\partial_t (\rho_0 \tilde{\mathbf{u}}) + \nabla \tilde{p} + \mathbf{B}_0 \times \nabla \times \tilde{\mathbf{B}} + \tilde{\mathbf{B}} \times \nabla \times \mathbf{B}_0 - \nabla \cdot \left( (\mu + \mu_a(\mathbf{x})) \nabla \tilde{\mathbf{u}} \right) = 0\,, \\[2mm] &\partial_t \tilde{p} + \tilde{\mathbf{u}} \cdot \nabla p_0 + \gamma p_0 \nabla \cdot \tilde{\mathbf{u}} = 0\,, \\[2mm] &\partial_t \tilde{\mathbf{B}} + \nabla \times (\mathbf{B}_0 \times \tilde{\mathbf{u}}) + \nabla \times (\eta + \eta_a(\mathbf{x})) \nabla \times \tilde{\mathbf{B}} = 0\,.\end{split}\]

In the linear model the viscous and resistive heating is not fed back into the pressure (the dissipative propagators with model="linear_p" update only \(\tilde{\mathbf{u}}\) and \(\tilde{\mathbf{B}}\)).

classmethod doc_normalization()

The velocity scale is Alfvénic:

\[\hat u = \hat v_A.\]

classmethod doc_scalar_quantities()

The following scalars are tracked during simulation:

• Kinetic energy: en_U

• Thermal energies: en_thermo, en_thermo_l1

• Magnetic energies: en_mag_1, en_mag_2, en_mag_l1

• Total energies: en_tot, en_tot_l1

classmethod doc_discretization()

Time integration is performed by the following propagators (in sequence):

1. VariationalDensityEvolve

2. VariationalMomentumAdvection (if linear is False)

3. VariationalPBEvolve

4. VariationalViscosity (if with_viscosity is True)

5. VariationalResistivity (if with_resistivity is True)

The option model of the variational propagators is "perturbative" ("full_p" for the dissipative ones), or "linear" ("linear_p") if linear is True.

classmethod doc_long_description()

VariationalPerturbativeMHD is the dissipative perturbative MHD model using pressure as thermodynamic variable. It is intended for departures from a background equilibrium that are not strictly linear but still naturally formulated as perturbations; with linear=True it solves linear MHD around the equilibrium.

classmethod doc_examples()

Create and initialize the visco-resistive perturbative MHD model and its linear version:

from struphy.models import VariationalPerturbativeMHD

model = VariationalPerturbativeMHD()
linear_model = VariationalPerturbativeMHD(linear=True)
model.mhd.pressure
model.em_fields.b_field

classmethod doc_use_cases()

This model is appropriate for:

• perturbative nonlinear MHD around a prescribed equilibrium

• linear MHD waves and instabilities (linear=True)

• comparing linear and perturbative variational formulations

classmethod doc_cannot_be_used_for()

This model is not suitable for:

• full-f MHD evolution without a reference equilibrium

• kinetic or particle-resolved plasma effects

• entropy- or q-based thermodynamic formulations

class struphy.models.VariationalPerturbativeMHD_with_q(base_units: BaseUnits = BaseUnits(x=1.0, B=1.0, n=1.0, kBT=None), mass_number: float = 1.0, linear: bool = False, with_viscosity: bool = True, with_resistivity: bool = True)

Bases: StruphyModel

Perturbative (delta-f) or linear visco-resistive MHD equations, with the q variable (square root of the pressure), discretized with a variational method.

Parameters:

• base_units (BaseUnits) – Base units for normalization (default: BaseUnits())

• mass_number (float) – Mass number (in units of Proton mass) of the fluid species (default: 1.0)

• linear (bool) – Whether to solve the linearized equations, i.e. to drop the terms that are nonlinear in the perturbations (default: False)

• with_viscosity (bool) – Whether to include viscous dissipation (default: True)

• with_resistivity (bool) – Whether to include resistive dissipation (default: True)

classmethod doc_pde()

PDEs solved by model:

The unknowns are the perturbations \((\tilde\rho, \tilde q, \tilde{\mathbf B})\) with respect to a zero-flow MHD equilibrium \((\rho_0, q_0, \mathbf B_0)\) with \(q_0 = \sqrt{p_0}\), and the velocity \(\mathbf u\). The total quantities are \(\rho = \rho_0 + \tilde\rho\), \(q = q_0 + \tilde q = \sqrt{p}\) and \(\mathbf B = \mathbf B_0 + \tilde{\mathbf B}\).

Continuity:

\[\partial_t \tilde{\rho} + \nabla \cdot \left( (\rho_0 + \tilde{\rho}) \mathbf{u} \right) = 0\]

Momentum:

\[\partial_t (\rho \mathbf{u}) + \nabla \cdot (\rho\, \mathbf{u} \otimes \mathbf{u}) + \nabla \left( 2 q_0 \tilde{q} + \tilde{q}^2 \right) = (\nabla \times \tilde{\mathbf{B}}) \times \mathbf{B}_0 + (\nabla \times \mathbf{B}_0) \times \tilde{\mathbf{B}} + (\nabla \times \tilde{\mathbf{B}}) \times \tilde{\mathbf{B}} + \nabla \cdot \left( (\mu + \mu_a(\mathbf{x})) \nabla \mathbf{u} \right)\]

Square root of the pressure:

\[\partial_t \tilde{q} + \mathbf{u} \cdot \nabla (q_0 + \tilde{q}) + \frac{\gamma}{2} (q_0 + \tilde{q}) \nabla \cdot \mathbf{u} = 0\]

Induction:

\[\partial_t \tilde{\mathbf{B}} - \nabla \times \left( \mathbf{u} \times (\mathbf{B}_0 + \tilde{\mathbf{B}}) \right) + \nabla \times \left( (\eta + \eta_a(\mathbf{x})) \nabla \times \tilde{\mathbf{B}} \right) = 0\]

Here \(\mu_a(\mathbf{x})\) and \(\eta_a(\mathbf{x})\) are artificial viscosity and resistivity coefficients. Since \(\nabla(2 q_0 \tilde q + \tilde q^2) = \nabla \tilde p\) with \(\tilde p = p - p_0\), the ideal part (\(\mu = \eta = 0\)) is the perturbative ideal MHD system written in \(q\) instead of \(p\).

Linear model (linear=True): the terms that are nonlinear in the perturbations are dropped, i.e. \(\rho \to \rho_0\) in the continuity and momentum equations (no \(\mathbf u \cdot \nabla \mathbf u\)), the pressure force is \(\nabla (2 q_0 \tilde q) = 2 q_0 \nabla \tilde q + 2 \tilde q \nabla q_0\), the Lorentz force has no \((\nabla \times \tilde{\mathbf B}) \times \tilde{\mathbf B}\), and \(\tilde q\), \(\tilde{\mathbf B}\) are advected with \(q_0\) and \(\mathbf B_0\) only:

\[\partial_t \tilde{q} + \mathbf{u} \cdot \nabla q_0 + \frac{\gamma}{2} q_0 \nabla \cdot \mathbf{u} = 0\,, \qquad \partial_t \tilde{\mathbf{B}} - \nabla \times (\mathbf{u} \times \mathbf{B}_0) + \nabla \times \left( (\eta + \eta_a(\mathbf{x})) \nabla \times \tilde{\mathbf{B}} \right) = 0\,.\]

The equations derive from the Lagrangian (with \(\rho_0\) instead of \(\rho\) in the linear model)

\[L = \int_\Omega \left( \rho \frac{|\mathbf{u}|^2}{2} - \frac{\tilde{q}^2 + 2 q_0 \hat{q}}{\gamma - 1} - \frac{|\tilde{\mathbf{B}}|^2}{2} - \mathbf{B}_0 \cdot \hat{\mathbf{B}} \right) \textrm d \mathbf{x}\,,\]

with the auxiliary variables \(\hat q\) (diagnostics.qt3) and \(\hat{\mathbf B}\) (diagnostics.bt2) advected by the perturbations only (in both models),

\[\partial_t \hat{q} + \mathbf{u} \cdot \nabla \tilde{q} + \frac{\gamma}{2} \tilde{q} \nabla \cdot \mathbf{u} = \frac{\gamma - 1}{2 q_0} \left( (\mu + \mu_a(\mathbf{x})) |\nabla \mathbf{u}|^2 + (\eta + \eta_a(\mathbf{x})) |\nabla \times \tilde{\mathbf{B}}|^2 \right)\,, \qquad \partial_t \hat{\mathbf{B}} - \nabla \times (\mathbf{u} \times \tilde{\mathbf{B}}) = 0\,.\]

They do not feed back on the other unknowns, but they make the energy

\[H = \int_\Omega \left( \rho \frac{|\mathbf{u}|^2}{2} + \frac{\tilde{q}^2 + 2 q_0 \hat{q}}{\gamma - 1} + \frac{|\tilde{\mathbf{B}}|^2}{2} + \mathbf{B}_0 \cdot \hat{\mathbf{B}} \right) \textrm d \mathbf{x}\]

a conserved quantity. Note that the heat produced by viscosity and resistivity goes into \(\hat q\), not into \(\tilde q\) (in the linear model, the heating is skipped with the option fast=True of the dissipative propagators).

classmethod doc_normalization()

The model uses Alfvén-speed normalization:

\[\hat u = \hat v_A.\]

classmethod doc_scalar_quantities()

The following scalars are tracked during simulation:

• Kinetic energy: en_U \(= \int \rho |\mathbf u|^2/2\) (\(\rho_0\) in the linear model)

• Magnetic energies: en_mag_1 \(= \int |\tilde{\mathbf B}|^2/2\), en_mag_2 \(= \int \mathbf B_0 \cdot \hat{\mathbf B}\)

• Thermodynamic energies: en_thermo_1 \(= \int \tilde q^2/(\gamma - 1)\), en_thermo_2 \(= \int 2 q_0 \hat q/(\gamma - 1)\)

• Total energy: en_tot (conserved)

classmethod doc_discretization()

Time integration is performed by the following propagators (in sequence):

1. VariationalDensityEvolve

2. VariationalMomentumAdvection (if linear is False)

3. VariationalQBEvolve

4. VariationalViscosity (if with_viscosity is True)

5. VariationalResistivity (if with_resistivity is True)

The option model of the variational propagators is "perturbative_q", or "linear_q" if linear is True.

classmethod doc_long_description()

MHD written for the perturbations with respect to a zero-flow equilibrium, in the variable \(q = \sqrt{p}\), with optional viscosity and resistivity. By default all nonlinear terms are kept, so the model is equivalent to full MHD, but the equilibrium force balance is subtracted analytically: a zero perturbation stays exactly zero, also when the discrete equilibrium is not in exact force balance. With linear=True the model solves linear MHD around the equilibrium. Using \(q\) makes the Lagrangian quadratic in the perturbed variables, which allows the variational discretization with the auxiliary variables \(\hat q\) and \(\hat{\mathbf B}\). Set with_viscosity=False and with_resistivity=False for ideal MHD.

classmethod doc_examples()

Create the ideal perturbative q-MHD model and its linear version:

from struphy.models import VariationalPerturbativeMHD_with_q

model = VariationalPerturbativeMHD_with_q(with_viscosity=False, with_resistivity=False)
linear_model = VariationalPerturbativeMHD_with_q(linear=True, with_viscosity=False, with_resistivity=False)
model.mhd.sqrt_p
model.em_fields.b_field

classmethod doc_use_cases()

This model is appropriate for:

• nonlinear MHD dynamics around a given (possibly non-uniform) zero-flow equilibrium

• linear MHD waves and instabilities (linear=True)

• instabilities and waves where the discrete equilibrium force must not drive spurious flows

• comparison with VariationalMHD_with_q (full-f), with VariationalPerturbativeMHD (pressure) and with LinearMHD

classmethod doc_cannot_be_used_for()

This model is not suitable for:

• equilibria with flow

• dissipative runs where the viscous and resistive heating must raise the pressure (it goes into \(\hat q\), see doc_pde)

• entropy-based thermodynamic evolution

• kinetic or hybrid particle coupling

Multi-fluid

class struphy.models.ColdPlasma(base_units: BaseUnits = BaseUnits(x=1.0, B=1.0, n=1.0, kBT=None), charge_number: int = -1, mass_number: float = 0.0005446623093681918, alpha: float | None = None, epsilon: float | None = None)

Bases: StruphyModel

Cold plasma model: electron-fluid current coupled with Maxwell’s equations via a cold-plasma Ohm’s law.

Parameters:

• base_units (BaseUnits) – Base units for normalization (default: BaseUnits())

• charge_number (int) – Charge number (in units of the positive elementary charge) of the electron species (default: -1)

• mass_number (float) – Mass number (in units of Proton mass) of the electron species (default: 1/1836)

• alpha (float, optional) – Dimensionless parameter: plasma frequency / cyclotron frequency. If None, computed from units and charge/mass numbers.

• epsilon (float, optional) – Normalized cyclotron period: 1 / (cyclotron frequency × time unit). If None, computed from units and charge/mass numbers.

classmethod doc_pde()

PDEs solved by model:

Cold-plasma current:

\[\frac{1}{n_0} \frac{\partial \mathbf{j}}{\partial t} = \frac{1}{\varepsilon} \mathbf{E} + \frac{1}{\varepsilon n_0} \mathbf{j} \times \mathbf{B}_0\]

Faraday’s law:

\[\frac{\partial \mathbf{B}}{\partial t} + \nabla \times \mathbf{E} = 0\]

Ampère’s law:

\[-\frac{\partial \mathbf{E}}{\partial t} + \nabla \times \mathbf{B} = \frac{\alpha^2}{\varepsilon} \mathbf{j}\]

where \((n_0, \mathbf{B}_0)\) denotes an inhomogeneous background.

classmethod doc_normalization()

Velocities are normalized with the speed of light and the fields satisfy

\[\hat v = c,\qquad \hat E = c \hat B.\]

The dimensionless plasma parameters are the cold-species \(\alpha = \hat\Omega_\mathrm{p}/\hat\Omega_\mathrm{c}\) and \(\varepsilon = 1/(\hat\Omega_\mathrm{c}\hat t)\).

classmethod doc_scalar_quantities()

The following scalars are tracked during simulation:

• Electric field energy: electric_energy

• Magnetic field energy: magnetic_energy

• Cold-current energy: kinetic_energy

• Total energy: total_energy

classmethod doc_discretization()

Time integration is performed by the following propagators (in sequence):

1. MaxwellWeakAmpere

2. OhmCold

3. JxBCold

classmethod doc_long_description()

This is a fluid electromagnetic model for a cold plasma response. The electron fluid is represented only through its current, so thermal pressure and kinetic velocity-space effects are omitted. It is useful as a reduced model between vacuum Maxwell and fully kinetic Vlasov-Maxwell dynamics.

classmethod doc_examples()

Create and initialize a cold-plasma model:

from struphy.models import ColdPlasma

model = ColdPlasma()
model.em_fields.e_field
model.em_fields.b_field
model.electrons.current

classmethod doc_use_cases()

This model is appropriate for:

• cold-plasma wave propagation studies

• electromagnetic benchmarks with a fluid current response

• regimes where thermal pressure can be neglected

• algorithm verification for Maxwell plus current coupling

classmethod doc_cannot_be_used_for()

This model is not suitable for:

• finite-temperature or pressure-driven plasma dynamics

• kinetic resonances or velocity-space instabilities

• multi-species hybrid or fully kinetic problems

• collisional closures beyond the built-in cold-plasma approximation

class struphy.models.TwoFluidQuasiNeutralToy(base_units: BaseUnits = BaseUnits(x=1.0, B=1.0, n=1.0, kBT=1.0), ion_charge_number: int = 1, ion_mass_number: float = 1.0, ion_epsilon: float | None = None, electron_charge_number: int = 1, electron_mass_number: float = 1.0, electron_epsilon: float | None = None)

Bases: StruphyModel

Linearized, quasi-neutral two-fluid model with zero electron inertia.

Parameters:

• base_units (BaseUnits) – Base units for normalization (default: BaseUnits(kBT=1.0))

• ion_charge_number (int) – Charge number (in units of the positive elementary charge) of the ion species (default: 1)

• ion_mass_number (float) – Mass number (in units of Proton mass) of the ion species (default: 1.0)

• ion_epsilon (float, optional) – Normalized cyclotron period of the ion species. If None, computed from units and charge/mass numbers.

• electron_charge_number (int) – Charge number (in units of the positive elementary charge) of the electron species (default: 1)

• electron_mass_number (float) – Mass number (in units of Proton mass) of the electron species (default: 1.0)

• electron_epsilon (float, optional) – Normalized cyclotron period of the electron species. If None, computed from units and charge/mass numbers.

classmethod doc_pde()

PDEs solved by model:

Ion momentum:

\[\frac{\partial \mathbf{u}}{\partial t} = -\nabla \phi + \frac{\mathbf{u} \times \mathbf{B}_0}{\varepsilon} + \nu \Delta \mathbf{u} + \mathbf{f}\]

Electron momentum:

\[0 = \nabla \phi - \frac{\mathbf{u}_e \times \mathbf{B}_0}{\varepsilon} + \nu_e \Delta \mathbf{u}_e + \mathbf{f}_e\]

Quasi-neutrality constraint:

\[\nabla \cdot (\mathbf{u} - \mathbf{u}_e) = 0\]

where \(\mathbf{B}_0\) is a static magnetic field and \(\mathbf{f}, \mathbf{f}_e\) are given forcing terms.

classmethod doc_normalization()

Thermal-speed scaling is used:

\[\hat u = \hat v_\mathrm{th},\qquad e\hat\phi = m \hat v_\mathrm{th}^2.\]

classmethod doc_scalar_quantities()

The following scalars are tracked during simulation:

• No default scalar diagnostics are defined by this model.

classmethod doc_discretization()

Time integration is performed by the following propagators (in sequence):

1. TwoFluidQuasiNeutralFull

classmethod doc_long_description()

TwoFluidQuasiNeutralToy is a reduced linear two-fluid benchmark with zero electron inertia. It is meant for studying the quasi-neutral solve and the coupled ion/electron velocity response in a simplified setting.

References

[1] Juan Vicente Gutiérrez-Santacreu, Omar Maj, Marco Restelli: Finite element discretization of a Stokes-like model arising in plasma physics, Journal of Computational Physics 2018.

classmethod doc_examples()

Create and initialize the quasi-neutral toy model:

from struphy.models import TwoFluidQuasiNeutralToy

model = TwoFluidQuasiNeutralToy()
model.em_fields.phi
model.ions.u
model.electrons.u

classmethod doc_use_cases()

This model is appropriate for:

• linear quasi-neutral two-fluid benchmarks

• Stokes-like plasma model verification

• testing coupled velocity-potential FEEC solvers

classmethod doc_cannot_be_used_for()

This model is not suitable for:

• nonlinear two-fluid dynamics

• finite electron inertia effects

• kinetic phase-space phenomena

• self-consistent electromagnetic wave propagation

class struphy.models.PressureLessSPH(base_units: BaseUnits = BaseUnits(x=1.0, B=1.0, n=1.0, kBT=None), charge_number: int = 1, mass_number: float = 1.0, epsilon: float | None = None)

Bases: StruphyModel

Particle discretization of pressureless Euler flow with external forcing.

Parameters:

• base_units (BaseUnits, optional) – Reference units used to derive model normalization constants.

• charge_number (int, optional) – Species charge number (default is 1).

• mass_number (float, optional) – Species mass number (default is 1.0).

• epsilon (float, optional) – Scaling parameter for force field (default is None).

classmethod doc_pde()

PDEs solved by model:

Continuity:

\[\partial_t \rho + \nabla \cdot (\rho \mathbf{u}) = 0\]

Momentum:

\[\partial_t (\rho \mathbf{u}) + \nabla \cdot (\rho \mathbf{u} \otimes \mathbf{u}) = \mathbf{F}\]

where \(\mathbf{F}\) is an external force.

classmethod doc_normalization()

Velocity and field normalizations:

\[\hat{u} = 1\,\textrm{m/s}, \qquad \hat{F} = \frac{m\hat{n}\hat{u}}{\hat{t}}\,.\]

classmethod doc_scalar_quantities()

The following scalars are tracked during simulation:

• SPH kinetic energy: en_kin

classmethod doc_discretization()

Time integration is performed by the following propagators (in sequence):

1. PushEta

2. PushVinForceField

classmethod doc_long_description()

PressureLessSPH is a meshfree particle model for a pressureless fluid. It is primarily useful as a simple SPH benchmark or as a reduced particle transport model in a prescribed force field or potential.

classmethod doc_examples()

Create and initialize a pressureless SPH model:

from struphy.models import PressureLessSPH

model = PressureLessSPH()
model.cold_fluid.var

classmethod doc_use_cases()

This model is appropriate for:

• meshfree pressureless-flow benchmarks

• straight-line particle transport with external forcing

• testing SPH particle loading and diagnostics

classmethod doc_cannot_be_used_for()

This model is not suitable for:

• compressible fluids with pressure forces

• FEEC-based grid discretizations

• electromagnetic plasma dynamics

• viscous or thermal closures

class struphy.models.ViscousEulerSPH(base_units: BaseUnits = BaseUnits(x=1.0, B=1.0, n=1.0, kBT=1.0), charge_number: int = 1, mass_number: float = 1.0, with_B0: bool = True, with_p: bool = True, with_viscosity: bool = True)

Bases: StruphyModel

Euler equations with viscosity discretized with smoothed particle hydrodynamics (SPH).

Parameters:

• base_units (BaseUnits) – Base units for normalization (default: BaseUnits(kBT=1.0))

• charge_number (int) – Charge number (in units of the positive elementary charge) of the fluid species (default: 1)

• mass_number (float) – Mass number (in units of Proton mass) of the fluid species (default: 1.0)

• with_B0 (bool) – Whether to include the effect of a background magnetic field B0 (default: True)

• with_p (bool) – Whether to include pressure forces (default: True)

• with_viscosity (bool) – Whether to include viscous dissipation (default: True)

classmethod doc_pde()

PDEs solved by model:

Continuity:

\[\partial_t \rho + \nabla \cdot (\rho \mathbf{u}) = 0\]

Momentum:

\[\rho (\partial_t \mathbf{u} + \mathbf{u} \cdot \nabla \mathbf{u}) = -\nabla \left( \rho^2 \frac{\partial \mathcal{U}(\rho, S)}{\partial \rho} \right) - \nabla \cdot \boldsymbol{\pi}\]

Entropy:

\[\partial_t S + \mathbf{u} \cdot \nabla S = 0\]

where \(S\) denotes the entropy per unit mass and \(\boldsymbol{\pi}\) is the viscous stress tensor.

The viscous stress tensor for a Newtonian fluid is given by:

\[\boldsymbol{\sigma} = -\mu \left( \nabla \mathbf{u} + (\nabla \mathbf{u})^T - \frac{2}{3} (\nabla \cdot \mathbf{u}) \mathbf{I} \right)\]

where \(\mu\) is the dynamic (shear) viscosity and \(\mathbf{I}\) is the identity tensor.

The internal energy per unit mass can be defined in two ways:

\[\mathrm{isothermal:} \qquad \mathcal{U}(\rho, S) = \kappa(S) \log \rho\]

\[\mathrm{polytropic:} \qquad \mathcal{U}(\rho, S) = \kappa(S) \frac{\rho^{\gamma - 1}}{\gamma - 1}\]

classmethod doc_normalization()

The characteristic speed is thermal:

\[\hat u = \hat v_\mathrm{th}.\]

classmethod doc_scalar_quantities()

The following scalars are tracked during simulation:

• SPH kinetic energy: en_kin

classmethod doc_discretization()

Time integration is performed by the following propagators (in sequence):

1. PushEta

2. PushVxB (if with_B0 is True)

3. PushVinSPHpressure (if with_p is True)

4. PushVinViscousPotential (if with_viscosity is True)

classmethod doc_long_description()

ViscousEulerSPH is the particle-based SPH fluid model with optional pressure and viscosity contributions. It is intended for meshfree fluid experiments rather than FEEC-based field simulations.

classmethod doc_examples()

Create and initialize the viscous Euler SPH model:

from struphy.models import ViscousEulerSPH

model = ViscousEulerSPH()
model.euler_fluid.var

classmethod doc_use_cases()

This model is appropriate for:

• SPH verification with pressure and viscosity

• meshfree compressible-fluid experiments

• testing particle-based viscosity and pressure pushers

classmethod doc_cannot_be_used_for()

This model is not suitable for:

• FEEC grid-based fluid or MHD formulations

• entropy-resolved thermodynamic evolution

• kinetic plasma physics

• studies that require exact field-based conservation structures

class struphy.models.IncompressibleNavierStokesSPH(base_units: BaseUnits = BaseUnits(x=1.0, B=1.0, n=1.0, kBT=1.0), charge_number: int = 1, mass_number: float = 1.0, with_B0: bool = True, with_viscosity: bool = True)

Bases: StruphyModel

Incompressible Navier-Stokes equations discretized with smoothed particle hydrodynamics (SPH).

Parameters:

• base_units (BaseUnits) – Base units for normalization (default: BaseUnits(kBT=1.0))

• charge_number (int) – Charge number (in units of the positive elementary charge) of the fluid species (default: 1)

• mass_number (float) – Mass number (in units of Proton mass) of the fluid species (default: 1.0)

• with_B0 (bool) – Whether to include the effect of a background magnetic field B0 (default: True)

• with_viscosity (bool) – Whether to include viscous dissipation (default: True)

classmethod doc_pde()

PDEs solved by model:

Incompressible Navier-Stokes:

\[\begin{split}\rho (\partial_t \mathbf{u} + \mathbf{u} \cdot \nabla \mathbf{u}) = -\nabla p - \nabla \cdot \boldsymbol{\sigma} \\[2mm] \nabla \cdot \mathbf{u} = 0\,.\end{split}\]

where \(\boldsymbol{\sigma}\) is the viscous stress tensor.

The viscous stress tensor for a Newtonian fluid is given by:

\[\boldsymbol{\sigma} = -\mu \left( \nabla \mathbf{u} + (\nabla \mathbf{u})^T - \frac{2}{3} (\nabla \cdot \mathbf{u}) \mathbf{I} \right)\]

where \(\mu\) is the dynamic (shear) viscosity and \(\mathbf{I}\) is the identity tensor.

classmethod doc_normalization()

The characteristic speed is thermal:

\[\hat u = \hat v_\mathrm{th}.\]

classmethod doc_scalar_quantities()

The following scalars are tracked during simulation:

• SPH kinetic energy: en_kin

classmethod doc_discretization()

Time integration is performed by the following propagators (in sequence):

1. PushEta

2. PushVxB (if with_B0 is True)

3. PushVinViscousPotential (if with_viscosity is True)

4. PoissonSolve

5. PushVinForceField

classmethod doc_long_description()

IncompressibleNavierStokesSPH is the particle-based SPH fluid model…

classmethod doc_examples()

Create and initialize the viscous Euler SPH model:

from struphy.models import ViscousEulerSPH

model = ViscousEulerSPH()
model.euler_fluid.var

classmethod doc_use_cases()

This model is appropriate for:

• SPH verification with pressure and viscosity

• meshfree compressible-fluid experiments

• testing particle-based viscosity and pressure pushers

classmethod doc_cannot_be_used_for()

This model is not suitable for:

• FEEC grid-based fluid or MHD formulations

• entropy-resolved thermodynamic evolution

• kinetic plasma physics

• studies that require exact field-based conservation structures

class struphy.models.HasegawaWakatani(base_units: BaseUnits = BaseUnits(x=1.0, B=1.0, n=1.0, kBT=None), mass_number: float = 1.0)

Bases: StruphyModel

Hasegawa-Wakatani equations in 2D for drift-wave turbulence.

Parameters:

• base_units (BaseUnits) – Base units for normalization (default: BaseUnits())

• mass_number (float) – Mass number (in units of Proton mass) of the species (default: 1.0)

classmethod doc_pde()

PDEs solved by model:

Density equation:

\[\frac{\partial n}{\partial t} = C (\phi - n) - [\phi, n] - \kappa \, \partial_y \phi + \nu \, \nabla^{2N} n\]

Vorticity equation:

\[\frac{\partial \omega}{\partial t} = C (\phi - n) - [\phi, \omega] + \nu \, \nabla^{2N} \omega\]

Potential equation:

\[\Delta \phi = \omega\]

where \([\phi, n] = \partial_x \phi \partial_y n - \partial_y \phi \partial_x n\), \(C = C(x, y)\) and \(\kappa\) and \(\nu\) are constants (at the moment only \(N=1\) is available).

classmethod doc_normalization()

The electrostatic potential is scaled by the thermal-speed unit,

\[\hat u = \hat v_\mathrm{th},\qquad \hat\phi = \hat u \hat x.\]

classmethod doc_scalar_quantities()

The following scalars are tracked during simulation:

• No default scalar diagnostics are defined by this model.

classmethod doc_discretization()

Time integration is performed by the following propagators (in sequence):

1. PoissonSolve

2. HasegawaWakataniStep

classmethod doc_long_description()

This is a reduced 2D drift-wave turbulence model. The electrostatic potential is recovered from the vorticity at each step, and the coupled density-vorticity dynamics capture resistive drift-wave behavior without resolving full kinetic physics.

classmethod doc_examples()

Create and initialize a Hasegawa-Wakatani model:

from struphy.models import HasegawaWakatani

model = HasegawaWakatani()
model.em_fields.phi
model.plasma.density
model.plasma.vorticity

classmethod doc_use_cases()

This model is appropriate for:

• 2D resistive drift-wave turbulence studies

• reduced electrostatic edge-plasma benchmarks

• testing Poisson-coupled advection-diffusion solvers

classmethod doc_cannot_be_used_for()

This model is not suitable for:

• three-dimensional electromagnetic turbulence

• full kinetic Landau or cyclotron physics

• multi-species warm-fluid closures beyond the reduced HW system

• self-consistent magnetic-field evolution

class struphy.models.ViscousFluid(base_units: BaseUnits = BaseUnits(x=1.0, B=1.0, n=1.0, kBT=None), mass_number: float = 1.0, with_viscosity: bool = True)

Bases: StruphyModel

Full (non-linear) viscous Navier-Stokes equations discretized with a variational method.

Parameters:

• base_units (BaseUnits) – Base units for normalization (default: BaseUnits())

• mass_number (float) – Mass number (in units of Proton mass) of the fluid species (default: 1.0)

• with_viscosity (bool) – Whether to include viscous dissipation (default: True)

classmethod doc_pde()

PDEs solved by model:

Continuity:

\[\partial_t \rho + \nabla \cdot (\rho \mathbf{u}) = 0\]

Momentum:

\[\partial_t (\rho \mathbf{u}) + \nabla \cdot (\rho \mathbf{u} \otimes \mathbf{u}) + \rho \nabla \frac{(\rho \mathcal{U}(\rho, s))}{\partial \rho} + s \nabla \frac{(\rho \mathcal{U}(\rho, s))}{\partial s} - \nabla \cdot \left( (\mu + \mu_a(\mathbf{x})) \nabla \mathbf{u} \right) = 0\]

Entropy:

\[\partial_t s + \nabla \cdot (s \mathbf{u}) = \frac{1}{T} \left( (\mu + \mu_a(\mathbf{x})) |\nabla \mathbf{u}|^2 \right)\]

where the internal energy per unit mass is \(\mathcal U(\rho) = \rho^{\gamma-1} \exp(s / \rho)\). and \(\mu_a(\mathbf{x})\) is an artificial viscosity coefficient.

classmethod doc_normalization()

The model uses Alfvén-speed scaling for the velocity and entropy-based thermodynamic units for the internal-energy closure.

classmethod doc_scalar_quantities()

The following scalars are tracked during simulation:

• Kinetic energy: en_U

• Thermodynamic energy: en_thermo

• Total energy: en_tot

• Total density / entropy: dens_tot, entr_tot

classmethod doc_discretization()

Time integration is performed by the following propagators (in sequence):

1. VariationalDensityEvolve

2. VariationalMomentumAdvection

3. VariationalEntropyEvolve

4. VariationalViscosity (if with_viscosity is True)

classmethod doc_long_description()

ViscousFluid is the non-magnetic viscous member of the variational fluid family. It is suited for compressible hydrodynamics with entropy transport and viscous heating but without magnetic effects.

classmethod doc_examples()

Create and initialize a viscous-fluid model:

from struphy.models import ViscousFluid

model = ViscousFluid()
model.fluid.density
model.fluid.velocity
model.fluid.entropy

classmethod doc_use_cases()

This model is appropriate for:

• nonlinear viscous compressible hydrodynamics

• entropy-based Navier-Stokes benchmarks

• testing variational viscous closures without magnetism

classmethod doc_cannot_be_used_for()

This model is not suitable for:

• magnetic or MHD dynamics

• inviscid strictly conservative benchmarks

• kinetic particle effects

class struphy.models.VariationalBarotropicFluid(base_units: BaseUnits = BaseUnits(x=1.0, B=1.0, n=1.0, kBT=None), mass_number: float = 1.0)

Bases: StruphyModel

Barotropic fluid equations discretized with a variational method.

Parameters:

• base_units (BaseUnits) – Base units for normalization (default: BaseUnits())

• mass_number (float) – Mass number (in units of Proton mass) of the fluid species (default: 1.0)

classmethod doc_pde()

PDEs solved by model:

Continuity:

\[\partial_t \rho + \nabla \cdot (\rho \mathbf{u}) = 0\]

Momentum:

\[\partial_t (\rho \mathbf{u}) + \nabla \cdot (\rho \mathbf{u} \otimes \mathbf{u}) + \rho \nabla \frac{(\rho \mathcal{U}(\rho))}{\partial \rho} = 0\]

where the internal energy per unit mass is \(\mathcal U(\rho) = \rho/2\).

classmethod doc_normalization()

The reference flow speed is the Alfvén speed and the internal energy is scaled with density:

\[\hat u = \hat v_A,\qquad \hat{\mathcal U}=\hat\rho/2.\]

classmethod doc_scalar_quantities()

The following scalars are tracked during simulation:

• Kinetic energy: kinetic_energy

• Thermodynamic energy: thermo_energy

• Total energy: total_energy

classmethod doc_discretization()

Time integration is performed by the following propagators (in sequence):

1. VariationalDensityEvolve

2. VariationalMomentumAdvection

classmethod doc_long_description()

This is the simplest variational compressible-fluid model in Struphy. It keeps only density and velocity and uses a barotropic closure rather than a separate entropy or pressure evolution equation.

classmethod doc_examples()

Create and initialize a variational barotropic-fluid model:

from struphy.models import VariationalBarotropicFluid

model = VariationalBarotropicFluid()
model.fluid.density
model.fluid.velocity

classmethod doc_use_cases()

This model is appropriate for:

• barotropic compressible-flow benchmarks

• testing variational density and momentum propagators

• reduced fluid studies without entropy evolution

classmethod doc_cannot_be_used_for()

This model is not suitable for:

• non-barotropic thermodynamics

• magnetic-field dynamics

• viscous or resistive effects

• kinetic plasma phenomena

class struphy.models.VariationalCompressibleFluid(base_units: BaseUnits = BaseUnits(x=1.0, B=1.0, n=1.0, kBT=None), mass_number: float = 1.0, divdiv_alpha: float = 0.0)

Bases: StruphyModel

Fully compressible fluid equations discretized with a variational method.

Parameters:

• base_units (BaseUnits) – Base units for normalization (default: BaseUnits())

• mass_number (float) – Mass number (in units of Proton mass) of the fluid species (default: 1.0)

• divdiv_alpha (float) – Coefficient alpha >= 0 of the div-div regularization alpha int rho (div u)^2 dx of the kinetic energy (Hamiltonian regularization, penalizes compressions; default: 0.0), see kinetic_metric()

classmethod doc_pde()

PDEs solved by model:

Continuity:

\[\partial_t \rho + \nabla \cdot (\rho \mathbf{u}) = 0\]

Momentum:

\[\partial_t (\rho \mathbf{u}) + \nabla \cdot (\rho \mathbf{u} \otimes \mathbf{u}) + \rho \nabla \frac{(\rho \mathcal{U}(\rho, s))}{\partial \rho} + s \nabla \frac{(\rho \mathcal{U}(\rho, s))}{\partial s} = 0\]

Entropy:

\[\partial_t s + \nabla \cdot (s \mathbf{u}) = 0\]

where the internal energy per unit mass is \(\mathcal U(\rho) = \rho^{\gamma-1} \exp(s / \rho)\).

With 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\) (Hamiltonian regularization of compressions). The equations follow from the same variational principle with the momentum \(\rho \mathbf u - 2 \alpha \nabla (\rho \nabla \cdot \mathbf u)\) instead of \(\rho \mathbf u\), and the density derivative \(\alpha (\nabla \cdot \mathbf u)^2\) of the kinetic energy enters the force like \(|\mathbf u|^2/2\).

classmethod doc_normalization()

The model uses Alfvén-speed scaling for the flow variables together with separate units for internal energy and entropy:

\[\hat u = \hat v_A,\qquad \hat{\mathcal U}=K,\qquad \hat s=\hat\rho C_v.\]

classmethod doc_scalar_quantities()

The following scalars are tracked during simulation:

• Kinetic energy: en_U

• Thermodynamic energy: en_thermo

• Div-div part of the kinetic energy (only for divdiv_alpha > 0): en_divdiv

• Total energy: en_tot

classmethod doc_discretization()

Time integration is performed by the following propagators (in sequence):

1. VariationalDensityEvolve

2. VariationalMomentumAdvection

3. VariationalEntropyEvolve

classmethod doc_long_description()

VariationalCompressibleFluid is the entropy-based compressible fluid model in the variational family. It is the natural non-magnetic counterpart to the full variational MHD models.

classmethod doc_examples()

Create and initialize a variational compressible-fluid model:

from struphy.models import VariationalCompressibleFluid

model = VariationalCompressibleFluid()
model.fluid.density
model.fluid.velocity
model.fluid.entropy

classmethod doc_use_cases()

This model is appropriate for:

• compressible fluid benchmarks with entropy transport

• testing the variational fluid propagator stack

• non-magnetic hydrodynamics with conservative thermodynamics

classmethod doc_cannot_be_used_for()

This model is not suitable for:

• magnetic-field evolution or MHD coupling

• pressureless or barotropic-only reductions

• kinetic or particle-based transport physics

class struphy.models.VariationalPressurelessFluid(base_units: BaseUnits = BaseUnits(x=1.0, B=1.0, n=1.0, kBT=None), mass_number: float = 1.0)

Bases: StruphyModel

Pressure-less fluid equations discretized with a variational method.

Parameters:

• base_units (BaseUnits) – Base units for normalization (default: BaseUnits())

• mass_number (float) – Mass number (in units of Proton mass) of the fluid species (default: 1.0)

classmethod doc_pde()

PDEs solved by model:

Continuity:

\[\partial_t \rho + \nabla \cdot (\rho \mathbf{u}) = 0\]

Momentum:

\[\partial_t (\rho \mathbf{u}) + \nabla \cdot (\rho \mathbf{u} \otimes \mathbf{u}) = 0\]

classmethod doc_normalization()

The flow speed is normalized with the Alfvén speed:

\[\hat u = \hat v_A.\]

classmethod doc_scalar_quantities()

The following scalars are tracked during simulation:

• Kinetic energy: kinetic_energy

classmethod doc_discretization()

Time integration is performed by the following propagators (in sequence):

1. VariationalDensityEvolve

2. VariationalMomentumAdvection

classmethod doc_long_description()

This is the pressureless limit of the variational fluid hierarchy. It is intended as a reduced benchmark and as a simple transport model with conservative density and momentum updates.

classmethod doc_examples()

Create and initialize a pressureless variational-fluid model:

from struphy.models import VariationalPressurelessFluid

model = VariationalPressurelessFluid()
model.fluid.density
model.fluid.velocity

classmethod doc_use_cases()

This model is appropriate for:

• pressureless compressible benchmarks

• testing the minimal variational fluid update chain

• reduced transport problems without thermodynamics

classmethod doc_cannot_be_used_for()

This model is not suitable for:

• pressure- or entropy-driven flow

• magnetic-field coupling

• viscous/resistive dissipation

• kinetic particle physics

Fluid-kinetic hybrids

class struphy.models.ColdPlasmaVlasov(base_units: BaseUnits = BaseUnits(x=1.0, B=1.0, n=1.0, kBT=None), thermal_charge_number: int = -1, thermal_mass_number: float = 0.0005446623093681918, hot_charge_number: int = -1, hot_mass_number: float = 0.0005446623093681918, thermal_alpha: float | None = None, thermal_epsilon: float | None = None, hot_epsilon: float | None = None)

Bases: StruphyModel

Hybrid cold-plasma model: a cold electron fluid coupled with a hot kinetic species and Maxwell’s equations.

Parameters:

• base_units (BaseUnits) – Base units for normalization (default: BaseUnits())

• thermal_charge_number (int) – Charge number of the cold (fluid) electron species (default: -1)

• thermal_mass_number (float) – Mass number of the cold (fluid) electron species (default: 1/1836)

• hot_charge_number (int) – Charge number of the hot (kinetic) electron species (default: -1)

• hot_mass_number (float) – Mass number of the hot (kinetic) electron species (default: 1/1836)

• thermal_alpha (float, optional) – Dimensionless parameter: cold-species plasma frequency / cyclotron frequency. If None, computed from units and charge/mass numbers.

• thermal_epsilon (float, optional) – Normalized cyclotron period of the cold species. If None, computed from units and charge/mass numbers.

• hot_epsilon (float, optional) – Normalized cyclotron period of the hot species. If None, computed from units and charge/mass numbers.

classmethod doc_pde()

PDEs solved by model:

Hot Vlasov species:

\[\frac{\partial f}{\partial t} + \mathbf{v} \cdot \nabla f + \frac{1}{\varepsilon_\textnormal{h}} \Big[ \mathbf{E} + \mathbf{v} \times \left( \mathbf{B} + \mathbf{B}_0 \right) \Big] \cdot \frac{\partial f}{\partial \mathbf{v}} = 0\]

Cold-plasma current:

\[\frac{1}{n_0} \frac{\partial \mathbf{j}_\textnormal{c}}{\partial t} = \frac{1}{\varepsilon_\textnormal{c}} \mathbf{E} + \frac{1}{\varepsilon_\textnormal{c} n_0} \mathbf{j}_\textnormal{c} \times \mathbf{B}_0\]

Faraday’s law:

\[\frac{\partial \mathbf{B}}{\partial t} + \nabla \times \mathbf{E} = 0\]

Ampère’s law:

\[-\frac{\partial \mathbf{E}}{\partial t} + \nabla \times \mathbf{B} = \frac{\alpha^2}{\varepsilon_\textnormal{h}} \left( \mathbf{j}_\textnormal{c} + \int_{\mathbb{R}^3} \mathbf{v} f \, \text{d}^3 \mathbf{v} \right)\]

where \((n_0, \mathbf{B}_0)\) denotes an inhomogeneous background.

At initial time the Poisson equation is solved once to weakly satisfy the Gauss law:

\[\nabla \cdot \mathbf{E} = \nu \frac{\alpha^2}{\varepsilon_\textnormal{h}} \int_{\mathbb{R}^3} f \, \text{d}^3 \mathbf{v}\]

classmethod doc_normalization()

Velocities are normalized with the speed of light and the kinetic background is scaled as

\[\hat v = c,\qquad \hat E = c \hat B,\qquad \hat f = \hat n / c^3.\]

The model distinguishes cold and hot cyclotron scales through \(\varepsilon_\mathrm{c}\) and \(\varepsilon_\mathrm{h}\).

classmethod doc_scalar_quantities()

The following scalars are tracked during simulation:

• Electric field energy: en_E

• Magnetic field energy: en_B

• Cold-current energy: en_J

• Hot-particle kinetic energy: en_f

• Total energy: en_tot

classmethod doc_discretization()

Time integration is performed by the following propagators (in sequence):

1. MaxwellWeakAmpere

2. OhmCold

3. JxBCold

4. PushEta

5. PushVxB

6. VlasovAmpereCoupling

classmethod doc_long_description()

ColdPlasmaVlasov is an electromagnetic hybrid model that keeps a cold fluid response for the thermal electrons and a kinetic PIC description for a hot species. It is designed for problems where the energetic population matters kinetically but the bulk response can still be treated as cold.

classmethod doc_examples()

Create and initialize the cold-plasma plus Vlasov model:

from struphy.models import ColdPlasmaVlasov

model = ColdPlasmaVlasov()
model.em_fields.e_field
model.em_fields.b_field
model.thermal_elec.current
model.hot_elec.var

classmethod doc_use_cases()

This model is appropriate for:

• hybrid electromagnetic problems with one hot kinetic species

• energetic-particle interaction with a cold background plasma

• reduced-cost alternatives to fully kinetic multi-population models

• benchmarks of fluid-kinetic current coupling

classmethod doc_cannot_be_used_for()

This model is not suitable for:

• fully warm-fluid or pressure-anisotropic background dynamics

• all-species kinetic simulations

• collision operators or detailed dissipative closures

• electrostatic-only reductions where magnetic evolution is irrelevant

class struphy.models.LinearMHDVlasovCC(base_units: BaseUnits = BaseUnits(x=1.0, B=1.0, n=1.0, kBT=None), mhd_mass_number: float = 1.0, hot_charge_number: int = 1, hot_mass_number: float = 1.0, hot_epsilon: float | None = None)

Bases: StruphyModel

Hybrid linear MHD coupled with energetic ions (6D Vlasov) via the current-coupling scheme.

Parameters:

• base_units (BaseUnits) – Base units for normalization (default: BaseUnits())

• mhd_mass_number (float) – Mass number (in units of Proton mass) of the MHD bulk species (default: 1.0)

• hot_charge_number (int) – Charge number (in units of the positive elementary charge) of the energetic ion species (default: 1)

• hot_mass_number (float) – Mass number (in units of Proton mass) of the energetic ion species (default: 1.0)

• hot_epsilon (float, optional) – Normalized cyclotron period of the energetic ion species. If None, computed from units and charge/mass numbers.

classmethod doc_pde()

PDEs solved by model:

MHD continuity:

\[\frac{\partial \tilde{\rho}}{\partial t} + \nabla \cdot (\rho_0 \tilde{\mathbf{U}}) = 0\]

MHD momentum:

\[\rho_0 \frac{\partial \tilde{\mathbf{U}}}{\partial t} + \nabla \tilde{p} = (\nabla \times \tilde{\mathbf{B}}) \times \mathbf{B}_0 + \mathbf{J}_0 \times \tilde{\mathbf{B}} \color{blue} + \frac{A_\textnormal{h}}{A_\textnormal{b}} \frac{1}{\varepsilon} \left( n_\textnormal{h} \tilde{\mathbf{U}} - n_\textnormal{h} \mathbf{u}_\textnormal{h} \right) \times (\mathbf{B}_0 + \tilde{\mathbf{B}}) \color{black}\]

MHD pressure:

\[\frac{\partial \tilde{p}}{\partial t} + (\gamma - 1) \nabla \cdot (p_0 \tilde{\mathbf{U}}) + p_0 \nabla \cdot \tilde{\mathbf{U}} = 0\]

MHD induction:

\[\frac{\partial \tilde{\mathbf{B}}}{\partial t} = \nabla \times (\tilde{\mathbf{U}} \times \mathbf{B}_0), \qquad \nabla \cdot \tilde{\mathbf{B}} = 0\]

Energetic-particle Vlasov equation:

\[\frac{\partial f_\textnormal{h}}{\partial t} + \mathbf{v} \cdot \nabla f_\textnormal{h} + \frac{1}{\varepsilon} \left[ \color{blue} (\mathbf{B}_0 + \tilde{\mathbf{B}}) \times \tilde{\mathbf{U}} \color{black} + \mathbf{v} \times (\mathbf{B}_0 + \tilde{\mathbf{B}}) \right] \cdot \frac{\partial f_\textnormal{h}}{\partial \mathbf{v}} = 0\]

Energetic-particle moments:

\[n_\textnormal{h} = \int_{\mathbb{R}^3} f_\textnormal{h} \, \textnormal{d}^3 \mathbf{v}, \qquad n_\textnormal{h} \mathbf{u}_\textnormal{h} = \int_{\mathbb{R}^3} f_\textnormal{h} \mathbf{v} \, \textnormal{d}^3 \mathbf{v}\]

where \(\mathbf{J}_0 = \nabla\times\mathbf{B}_0\).

classmethod doc_normalization()

Both fluid and particle velocities are normalized with the Alfvén speed,

\[\hat U = \hat v = \hat v_A,\qquad \hat f_h = \hat n / \hat v_A^3.\]

The hot-species cyclotron parameter is \(\varepsilon = 1/(\hat\Omega_{c,\mathrm{hot}}\hat t)\).

classmethod doc_scalar_quantities()

The following scalars are tracked during simulation:

• MHD kinetic energy: en_U

• Thermal pressure energy: en_p

• Magnetic energy: en_B

• Energetic-particle energy: en_f

• Total energy: en_tot

• Lost particles: n_lost_particles

classmethod doc_discretization()

Time integration is performed by the following propagators (in sequence):

1. CurrentCoupling6DDensity

2. ShearAlfvenPropagator

3. CurrentCoupling6DCurrent

4. PushEta

5. PushVxB

6. Magnetosonic

classmethod doc_long_description()

LinearMHDVlasovCC is a current-coupling hybrid model for small perturbations around an MHD equilibrium. It is intended for energetic particle effects on linear MHD waves and instabilities without the cost of a fully kinetic background plasma.

classmethod doc_examples()

Create and initialize the linear MHD-Vlasov current-coupling model:

from struphy.models import LinearMHDVlasovCC

model = LinearMHDVlasovCC()
model.em_fields.b_field
model.mhd.velocity
model.energetic_ions.var

classmethod doc_use_cases()

This model is appropriate for:

• linear energetic-particle coupling to MHD modes

• current-coupling hybrid verification

• reduced-cost studies of fast-ion effects on wave propagation

classmethod doc_cannot_be_used_for()

This model is not suitable for:

• strongly nonlinear dynamics far from equilibrium

• full multi-species kinetic plasma evolution

• pressure-coupling closures

• collisional or dissipative MHD effects not present in the equations

class struphy.models.LinearMHDVlasovPC(base_units: BaseUnits = BaseUnits(x=1.0, B=1.0, n=1.0, kBT=None), mhd_mass_number: float = 1.0, hot_charge_number: int = 1, hot_mass_number: float = 1.0, hot_epsilon: float | None = None, turn_off: tuple[str, ...] = (None,))

Bases: StruphyModel

Hybrid linear MHD coupled with energetic ions (6D Vlasov) via the pressure-coupling scheme.

Parameters:

• base_units (BaseUnits) – Base units for normalization (default: BaseUnits())

• mhd_mass_number (float) – Mass number (in units of Proton mass) of the MHD bulk species (default: 1.0)

• hot_charge_number (int) – Charge number (in units of the positive elementary charge) of the energetic ion species (default: 1)

• hot_mass_number (float) – Mass number (in units of Proton mass) of the energetic ion species (default: 1.0)

• hot_epsilon (float, optional) – Normalized cyclotron period of the energetic ion species. If None, computed from units and charge/mass numbers.

• turn_off (tuple[str, ...]) – Names of coupling terms to turn off (default: (None,))

classmethod doc_pde()

PDEs solved by model:

MHD continuity:

\[\frac{\partial \tilde{\rho}}{\partial t} + \nabla \cdot (\rho_0 \tilde{\mathbf{U}}) = 0\]

MHD momentum:

\[\rho_0 \frac{\partial \tilde{\mathbf{U}}}{\partial t} + \nabla \tilde{p} + \frac{A_\textnormal{h}}{A_\textnormal{b}} \nabla \cdot \tilde{\mathbb{P}}_{\textnormal{h},\perp} = (\nabla \times \tilde{\mathbf{B}}) \times \mathbf{B}_0 + \mathbf{J}_0 \times \tilde{\mathbf{B}}\]

\[\mathbf{J}_0 = \nabla \times \mathbf{B}_0\]

MHD pressure:

\[\frac{\partial \tilde{p}}{\partial t} + \nabla \cdot (p_0 \tilde{\mathbf{U}}) + \frac{2}{3} p_0 \nabla \cdot \tilde{\mathbf{U}} = 0\]

MHD induction:

\[\frac{\partial \tilde{\mathbf{B}}}{\partial t} - \nabla \times (\tilde{\mathbf{U}} \times \mathbf{B}_0) = 0\]

Energetic-particle Vlasov equation:

\[\frac{\partial f_\textnormal{h}}{\partial t} + (\mathbf{v} + \tilde{\mathbf{U}}_\perp) \cdot \nabla f_\textnormal{h} + \left[ \frac{1}{\epsilon} \mathbf{v} \times (\mathbf{B}_0 + \tilde{\mathbf{B}}) - \nabla \tilde{\mathbf{U}}_\perp \cdot \mathbf{v} \right] \cdot \frac{\partial f_\textnormal{h}}{\partial \mathbf{v}} = 0\]

Perpendicular pressure tensor:

\[\tilde{\mathbb{P}}_{\textnormal{h},\perp} = \int \mathbf{v}_\perp \mathbf{v}_\perp^\top f_\textnormal{h} \, d \mathbf{v}\]

classmethod doc_normalization()

Fluid and hot-particle velocities are normalized with the bulk Alfvén speed. The kinetic pressure tensor is scaled consistently with \(A_h m_H \hat n \hat v_A^2\), and the hot cyclotron parameter is \(\varepsilon\).

classmethod doc_scalar_quantities()

The following scalars are tracked during simulation:

• MHD kinetic energy: en_U

• Thermal pressure energy: en_p

• Magnetic energy: en_B

• Energetic-particle kinetic energy: en_f

• Total energy: en_tot

• Lost particles: n_lost_particles

classmethod doc_discretization()

Time integration is performed by the following propagators (in sequence):

1. PushEtaPC

2. PushVxB

3. PressureCoupling6D

4. ShearAlfvenPropagator

5. Magnetosonic

classmethod doc_long_description()

LinearMHDVlasovPC is the pressure-coupling counterpart to the current-coupling hybrid model. It is targeted at linear problems where energetic-particle pressure anisotropy is the relevant feedback channel on the bulk MHD dynamics.

classmethod doc_examples()

Create and initialize the linear MHD-Vlasov pressure-coupling model:

from struphy.models import LinearMHDVlasovPC

model = LinearMHDVlasovPC()
model.em_fields.b_field
model.mhd.velocity
model.energetic_ions.var

classmethod doc_use_cases()

This model is appropriate for:

• linear pressure-coupling hybrid studies

• energetic-particle pressure feedback on MHD modes

• verification of PushEtaPC and pressure-coupling operators

classmethod doc_cannot_be_used_for()

This model is not suitable for:

• current-coupling closure studies

• nonlinear hybrid turbulence

• dissipative/resistive MHD

• fully kinetic treatment of the bulk plasma

class struphy.models.LinearMHDDriftkineticCC(base_units: BaseUnits = BaseUnits(x=1.0, B=1.0, n=1.0, kBT=None), mhd_mass_number: float = 1.0, hot_charge_number: int = 1, hot_mass_number: float = 1.0, hot_epsilon: float | None = None, turn_off: tuple[str, ...] = (None,))

Bases: StruphyModel

Hybrid linear ideal MHD coupled with energetic ions (5D drift-kinetic) via the current-coupling scheme.

Parameters:

• base_units (BaseUnits) – Base units for normalization (default: BaseUnits())

• mhd_mass_number (float) – Mass number (in units of Proton mass) of the MHD bulk species (default: 1.0)

• hot_charge_number (int) – Charge number (in units of the positive elementary charge) of the energetic ion species (default: 1)

• hot_mass_number (float) – Mass number (in units of Proton mass) of the energetic ion species (default: 1.0)

• hot_epsilon (float, optional) – Normalized cyclotron period of the energetic ion species. If None, computed from units and charge/mass numbers.

• turn_off (tuple[str, ...]) – Names of coupling terms to turn off (default: (None,))

classmethod doc_pde()

PDEs solved by model:

MHD continuity:

\[\frac{\partial \tilde{\rho}}{\partial t} + \nabla \cdot (\rho_0 \tilde{\mathbf{U}}) = 0\]

MHD momentum:

\[\rho_0 \frac{\partial \tilde{\mathbf{U}}}{\partial t} + \nabla \tilde{p} = (\nabla \times \tilde{\mathbf{B}}) \times \mathbf{B} + (\nabla \times \mathbf{B}_0) \times \tilde{\mathbf{B}} + \frac{A_\textnormal{h}}{A_\textnormal{b}} \left[ \frac{1}{\epsilon} n_\textnormal{gc} \tilde{\mathbf{U}} - \frac{1}{\epsilon} \mathbf{J}_\textnormal{gc} - \nabla \times \mathbf{M}_\textnormal{gc} \right] \times \mathbf{B}\]

MHD pressure:

\[\frac{\partial \tilde{p}}{\partial t} + \nabla \cdot (p_0 \tilde{\mathbf{U}}) + \frac{2}{3} p_0 \nabla \cdot \tilde{\mathbf{U}} = 0\]

MHD induction:

\[\frac{\partial \tilde{\mathbf{B}}}{\partial t} - \nabla \times (\tilde{\mathbf{U}} \times \mathbf{B}) = 0\]

Energetic-particle drift-kinetic equation:

\[\frac{\partial f_\textnormal{h}}{\partial t} + \frac{1}{B_\parallel^*} \left( v_\parallel \mathbf{B}^* - \mathbf{b}_0 \times \mathbf{E}^* \right) \cdot \nabla f_\textnormal{h} + \frac{1}{\epsilon} \frac{1}{B_\parallel^*} (\mathbf{B}^* \cdot \mathbf{E}^*) \frac{\partial f_\textnormal{h}}{\partial v_\parallel} = 0\]

Energetic-particle moments:

\[n_\textnormal{gc} = \int f_\textnormal{h} B_\parallel^* \, \textnormal{d} v_\parallel \textnormal{d} \mu\]

\[\mathbf{J}_\textnormal{gc} = \int \frac{f_\textnormal{h}}{B_\parallel^*} \left( v_\parallel \mathbf{B}^* - \mathbf{b}_0 \times \mathbf{E}^* \right) \, \textnormal{d} v_\parallel \textnormal{d} \mu\]

\[\mathbf{M}_\textnormal{gc} = -\int f_\textnormal{h} B_\parallel^* \mu \mathbf{b}_0 \, \textnormal{d} v_\parallel \textnormal{d} \mu\]

where

\[\begin{split}B^*_\parallel = \mathbf{b}_0 \cdot \mathbf{B}^* \\[2mm] \mathbf{B}^* &= \mathbf{B} + \epsilon v_\parallel \nabla \times \mathbf{b}_0 \\[2mm] \mathbf{E}^* &= -\tilde{\mathbf{U}} \times \mathbf{B} - \epsilon \mu \nabla (\mathbf{b}_0 \cdot \mathbf{B})\end{split}\]

classmethod doc_normalization()

The bulk and energetic-particle flow scales are normalized with the bulk Alfvén speed, while the magnetic moment carries its own unit:

\[\hat U = \hat v = \hat v_{A,\mathrm{bulk}},\qquad \hat\mu = A_h m_H \hat v_h^2 / \hat B.\]

classmethod doc_scalar_quantities()

The following scalars are tracked during simulation:

• MHD kinetic energy: en_U

• Thermal pressure energy: en_p

• Magnetic energy: en_B

• Parallel energetic-particle energy: en_fv

• Magnetic-moment energetic-particle energy: en_fB

• Total energy: en_tot

• Lost particles: n_lost_particles

classmethod doc_discretization()

Time integration is performed by the following propagators (in sequence):

1. PushGuidingCenterBxEstar

2. PushGuidingCenterParallel

3. ShearAlfvenCurrentCoupling5D

4. Magnetosonic

5. CurrentCoupling5DDensity

6. CurrentCoupling5DGradB

7. CurrentCoupling5DCurlb

The default parameter file of this model (struphy params LinearMHDDriftkineticCC) sets the energy-conserving discrete gradient methods of the reference discretization for the guiding-center pushers: algo="discrete_gradient_2nd_order" (implicit mid-point) for 1 and algo="discrete_gradient_1st_order_newton" (Itoh-Abe with Newton iterations) for 2. The explicit Runge-Kutta methods (algo="explicit") are cheaper but not energy conserving.

classmethod doc_long_description()

LinearMHDDriftkineticCC is the reduced-kinetic hybrid current-coupling model for energetic ions. It is useful when gyrophase averaging is acceptable but energetic-particle feedback on linear MHD still has to be retained.

classmethod doc_examples()

Create and initialize the linear MHD plus drift-kinetic CC model:

from struphy.models import LinearMHDDriftkineticCC

model = LinearMHDDriftkineticCC()
model.em_fields.b_field
model.mhd.velocity
model.energetic_ions.var

classmethod doc_use_cases()

This model is appropriate for:

• linear energetic-ion effects with guiding-center reduction

• current-coupling hybrid mode studies in strong magnetic fields

• verification of 5D hybrid coupling operators

classmethod doc_cannot_be_used_for()

This model is not suitable for:

• full 6D energetic-particle dynamics

• nonlinear hybrid turbulence

• resistive or viscous MHD closures

• regimes where gyrophase resolution is required

Particle diffusion

class struphy.models.DeterministicParticleDiffusion(base_units: BaseUnits = BaseUnits(x=1.0, B=1.0, n=1.0, kBT=None))

Bases: StruphyModel

Diffusion equation discretized with a deterministic particle method.

Parameters:

base_units (BaseUnits) – Base units for normalization (default: BaseUnits())

classmethod doc_pde()

PDEs solved by model:

Find \(u : \mathbb{R} \times \Omega \to \mathbb{R}^+\) such that

\[\frac{\partial u}{\partial t} + \nabla \cdot \left( \mathbf{F}(u) u \right) = 0, \qquad \mathbf{F}(u) = -\mathbb{D} \frac{\nabla u}{u}\]

where \(\mathbb{D} : \Omega \to \mathbb{R}^{3 \times 3}\) is a positive diffusion matrix. At the moment only matrices of the form \(D * Id\) are implemented, where \(D > 0\) is a positive diffusion coefficient.

classmethod doc_normalization()

The diffusion coefficient defines the normalization,

\[\hat D = \hat x^2 / \hat t.\]

No separate plasma velocity scale is used in this model.

classmethod doc_scalar_quantities()

The following scalars are tracked during simulation:

• No default scalar diagnostics are defined by this model.

classmethod doc_discretization()

Time integration is performed by the following propagators (in sequence):

1. PushDeterministicDiffusion

classmethod doc_long_description()

This is a particle discretization of diffusion where particles follow a deterministic drift derived from the current estimate of the density. It is intended mainly for diffusion-method development and verification rather than plasma dynamics.

classmethod doc_examples()

Create and initialize a deterministic diffusion model:

from struphy.models import DeterministicParticleDiffusion

model = DeterministicParticleDiffusion()
model.hydrogen.var

classmethod doc_use_cases()

This model is appropriate for:

• deterministic particle discretizations of diffusion

• transport benchmarks with positive scalar densities

• numerical comparison against stochastic diffusion methods

classmethod doc_cannot_be_used_for()

This model is not suitable for:

• electromagnetic plasma dynamics

• kinetic Vlasov problems in phase space

• nonlinear fluid systems with pressure or momentum evolution

• diffusion tensors outside the currently supported simplified forms

class struphy.models.RandomParticleDiffusion(base_units: BaseUnits = BaseUnits(x=1.0, B=1.0, n=1.0, kBT=None))

Bases: StruphyModel

Diffusion equation discretized with a random particle method, via a Wiener process.

Parameters:

base_units (BaseUnits) – Base units for normalization (default: BaseUnits())

classmethod doc_pde()

PDEs solved by model:

Find \(u : \mathbb{R} \times \Omega \to \mathbb{R}^+\) such that

\[\frac{\partial u}{\partial t} - D \, \Delta u = 0\]

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

classmethod doc_normalization()

The natural scaling is set by the diffusion coefficient:

\[\hat D = \hat x^2 / \hat t.\]

classmethod doc_scalar_quantities()

The following scalars are tracked during simulation:

• No default scalar diagnostics are defined by this model.

classmethod doc_discretization()

Time integration is performed by the following propagators (in sequence):

1. PushRandomDiffusion

classmethod doc_long_description()

This model is the stochastic counterpart of the deterministic particle diffusion model. It is intended for diffusion-method development and comparisons between random-walk and deterministic transport strategies.

classmethod doc_examples()

Create and initialize a random diffusion model:

from struphy.models import RandomParticleDiffusion

model = RandomParticleDiffusion()
model.hydrogen.var

classmethod doc_use_cases()

This model is appropriate for:

• stochastic particle diffusion benchmarks

• Monte-Carlo transport verification

• comparison against deterministic diffusion solvers

classmethod doc_cannot_be_used_for()

This model is not suitable for:

• electromagnetic or fluid plasma dynamics

• deterministic advection-dominated transport

• anisotropic plasma kinetics in phase space

Pure fields

class struphy.models.Poisson(base_units: BaseUnits = BaseUnits(x=1.0, B=1.0, n=1.0, kBT=None), with_t_dep_source=False)

Bases: StruphyModel

Weak discretization of Poisson’s equation with a diffusion matrix, stabilization and an optional time-dependent right-hand side.

Parameters:

• base_units (BaseUnits) – Base units for normalization (default: BaseUnits())

• with_t_dep_source (bool) – Whether the right-hand side source term is time-dependent (default: False)

classmethod doc_pde()

PDEs solved by model:

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

\[-\nabla \cdot D_0(\mathbf{x}) \nabla \phi + n_0(\mathbf{x}) \phi = \rho(t, \mathbf{x})\]

where \(n_0, \rho(t) : \Omega \to \mathbb{R}\) are real-valued functions, \(\rho(t)\) is parametrized by time \(t\), and \(D_0 : \Omega \to \mathbb{R}^{3 \times 3}\) is a positive matrix. Boundary terms from integration by parts are assumed to vanish.

classmethod doc_normalization()

The coefficient scaling is

\[\hat D = \hat n / \hat x^2,\qquad \hat\rho = \hat n.\]

No dedicated velocity normalization is used.

classmethod doc_scalar_quantities()

The following scalars are tracked during simulation:

• No default scalar diagnostics are defined by this model.

classmethod doc_discretization()

Time integration is performed by the following propagators (in sequence):

1. TimeDependentSource (if with_t_dep_source is True)

2. PoissonSolve

classmethod doc_long_description()

Poisson is the standalone elliptic field-solve model used for weak diffusion/Poisson problems. It is also the building block reused by other models for initial electrostatic solves.

classmethod doc_examples()

Create and initialize a Poisson model:

from struphy.models import Poisson

model = Poisson()
model.em_fields.phi
model.em_fields.source

classmethod doc_use_cases()

This model is appropriate for:

• elliptic benchmark problems

• electrostatic field solves with prescribed source terms

• testing weak Poisson discretizations and boundary handling

classmethod doc_cannot_be_used_for()

This model is not suitable for:

• hyperbolic time-dependent wave propagation

• self-consistent kinetic plasma evolution on its own

• magnetic-field dynamics or full Maxwell coupling

class struphy.models.CurlCurl(base_units: BaseUnits = BaseUnits(x=1.0, B=1.0, n=1.0, kBT=None), with_t_dep_source=False)

Bases: StruphyModel

Weak discretization of the curl-curl problem with a mass (stabilization) term and an optional time-dependent right-hand side.

Parameters:

• base_units (BaseUnits) – Base units for normalization (default: BaseUnits())

• with_t_dep_source (bool) – Whether the right-hand side source term is time-dependent (default: False)

classmethod doc_pde()

PDEs solved by model:

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

\[\nabla \times \nabla \times \mathbf E + \sigma\, \mathbf E = \mathbf J(t, \mathbf{x})\]

where \(\sigma > 0\) is a scalar and \(\mathbf J(t) : \Omega \to \mathbb{R}^3\) is a source parametrized by time \(t\). Boundary terms from integration by parts are assumed to vanish.

classmethod doc_normalization()

The coefficient scaling is

\[\hat \sigma = 1 / \hat x^2,\qquad \hat J = \hat E / \hat x^2.\]

No dedicated velocity normalization is used.

classmethod doc_scalar_quantities()

The following scalars are tracked during simulation:

• No default scalar diagnostics are defined by this model.

classmethod doc_discretization()

Time integration is performed by the following propagators (in sequence):

1. TimeDependentSource (if with_t_dep_source is True)

2. CurlCurlSolve

classmethod doc_long_description()

CurlCurl is the standalone elliptic field-solve model for the stabilized curl-curl problem in \(H(\textnormal{curl})\), e.g. from implicit time stepping of Maxwell’s equations or from magnetostatics. It is used to test and profile solvers and preconditioners for this operator, whose large kernel (the gradients) makes standard preconditioners ineffective.

classmethod doc_examples()

Create and initialize a CurlCurl model:

from struphy.models import CurlCurl

model = CurlCurl()
model.em_fields.e_field
model.em_fields.source

classmethod doc_use_cases()

This model is appropriate for:

• elliptic benchmark problems in \(H(\textnormal{curl})\)

• testing and profiling curl-curl solvers and preconditioners

• field solves with prescribed current sources

classmethod doc_cannot_be_used_for()

This model is not suitable for:

• hyperbolic time-dependent wave propagation

• self-consistent kinetic plasma evolution on its own

• the singular case \(\sigma = 0\) (it is regularized with \(\sigma = 10^{-14}\))

class struphy.models.GyrokineticPoisson(base_units: BaseUnits = BaseUnits(x=1.0, B=1.0, n=1.0, kBT=None), epsilon: float = 1.0, Z: int = 1, with_t_dep_source=False)

Bases: StruphyModel

Weak discretization of the gyrokinetic Poisson equation with adiabatic electrons and an optional time-dependent right-hand side.

Parameters:

• base_units (BaseUnits) – Base units for normalization (default: BaseUnits())

• epsilon (float) – Gyrokinetic parameter (default: 1.0)

• Z (int) – Charge number of the ions (default: 1)

• with_t_dep_source (bool) – Whether the right-hand side source term is time-dependent (default: False)

classmethod doc_pde()

PDEs solved by model:

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

\[\frac{1}{Z\epsilon^2} \frac{n_0}{T_0} \left( \phi - \langle \phi \rangle \right) - \nabla \cdot \left( \frac{n_0}{|B_0|^2} \nabla_\perp \phi \right) = \frac{1}{\epsilon} \rho(t, \mathbf{x})\,,\]

where \(n_0, T_0, |B_0| : \Omega \to \mathbb{R}\) are the equilibrium density, temperature and magnetic field strength, \(\nabla_\perp = (\mathbb{1} - \mathbf b_0 \mathbf b_0^\top) \nabla\) is the gradient perpendicular to the unit vector \(\mathbf b_0 = \mathbf B_0 / |B_0|\), and \(\langle \phi \rangle\) is the \(n_0/T_0\)-weighted average over the toroidal (cylindrical geometry) or the poloidal and toroidal (toroidal geometry) angle. \(\epsilon\) is the gyrokinetic parameter, \(Z\) the charge number of the ions, and the right-hand side \(\rho(t)\) is parametrized by time \(t\). Boundary terms from integration by parts are assumed to vanish.

classmethod doc_normalization()

The coefficient scaling is

\[\hat \rho = \hat n.\]

No dedicated velocity normalization is used.

classmethod doc_scalar_quantities()

The following scalars are tracked during simulation:

• No default scalar diagnostics are defined by this model.

classmethod doc_discretization()

Time integration is performed by the following propagators (in sequence):

1. TimeDependentSource (if with_t_dep_source is True)

2. GyrokineticPoissonSolve

classmethod doc_long_description()

GyrokineticPoisson is the standalone elliptic field-solve model for the gyrokinetic Poisson equation with adiabatic electrons, the field equation of DriftKineticElectrostaticAdiabatic. It solves the field equation with a prescribed (possibly time-dependent) source, e.g. for verification and profiling of the solver in a given MHD equilibrium.

classmethod doc_examples()

Create and initialize a GyrokineticPoisson model:

from struphy.models import GyrokineticPoisson

model = GyrokineticPoisson(epsilon=1.0, Z=1)
model.em_fields.phi
model.em_fields.source

classmethod doc_use_cases()

This model is appropriate for:

• verification and profiling of the gyrokinetic Poisson solver

• electrostatic potentials of a prescribed gyrokinetic charge density in an MHD equilibrium

classmethod doc_cannot_be_used_for()

This model is not suitable for:

• self-consistent kinetic plasma evolution on its own

• electromagnetic or magnetic-field dynamics

class struphy.models.Maxwell(base_units: BaseUnits = BaseUnits(x=1.0, B=1.0, n=1.0, kBT=None))

Bases: StruphyModel

Maxwell’s equations in vacuum for electromagnetic field evolution.

Parameters:

base_units (BaseUnits) – Base units for normalization (default: BaseUnits())

classmethod doc_pde()

PDEs solved by model:

Ampère’s law (no current):

\[\frac{\partial \mathbf{E}}{\partial t} - \nabla \times \mathbf{B} = 0\]

Faraday’s law:

\[\frac{\partial \mathbf{B}}{\partial t} + \nabla \times \mathbf{E} = 0\]

classmethod doc_normalization()

Velocity and fields are normalized as:

\[\hat v = c\,,\qquad \hat E = c \hat B\]

where \(c\) is the speed of light.

classmethod doc_scalar_quantities()

The following scalars are tracked during simulation:

• Electric energy: \(E_E = \frac{1}{2} \int |\mathbf E|^2 \, dV\)

• Magnetic energy: \(E_B = \frac{1}{2} \int |\mathbf B|^2 \, dV\)

• Total energy: \(E_{total} = E_E + E_B\)

classmethod doc_discretization()

Time integration is performed by the following propagators (in sequence):

1. MaxwellWeakAmpere

classmethod doc_long_description()

This model simulates the propagation of electromagnetic waves in vacuum using Maxwell’s equations without sources. It uses a finite element exterior calculus (FEEC) formulation with the electric field in H(curl) and the magnetic field in H(div) spaces.

classmethod doc_examples()

Create and initialize a Maxwell model:

from struphy.models import Maxwell

model = Maxwell()
# Fields are accessible via:
model.em_fields.e_field
model.em_fields.b_field

classmethod doc_use_cases()

Propagation of electromagnetic waves in vacuum.

classmethod doc_cannot_be_used_for()

Plasma dynamics, plasma-field interactions, or any scenario involving charged particles. This model does not include any particle species or coupling to matter.

Particle tracing

class struphy.models.Vlasov(base_units: BaseUnits = BaseUnits(x=1.0, B=1.0, n=1.0, kBT=None), charge_number: int = 1, mass_number: float = 1.0)

Bases: StruphyModel

Vlasov equation for a single species in a static background magnetic field.

Parameters:

• base_units (BaseUnits) – Base units for normalization (default: BaseUnits())

• charge_number (int) – Charge number (in units of the positive elementary charge) of the species (default: 1)

• mass_number (float) – Mass number (in units of Proton mass) of the species (default: 1.0)

classmethod doc_pde()

PDEs solved by model:

Vlasov equation:

\[\frac{\partial f}{\partial t} + \mathbf{v} \cdot \nabla f + \left( \mathbf{v} \times \mathbf{B}_0 \right) \cdot \frac{\partial f}{\partial \mathbf{v}} = 0\]

classmethod doc_normalization()

The characteristic speed is the cyclotron scale

\[\hat v = \hat\Omega_c \hat x.\]

classmethod doc_scalar_quantities()

The following scalars are tracked during simulation:

• Particle kinetic energy: kinetic_energy

classmethod doc_discretization()

Time integration is performed by the following propagators (in sequence):

1. PushVxB

2. PushEta

classmethod doc_long_description()

Vlasov is the simplest kinetic test-particle model in the 6D hierarchy. It evolves particles in a static magnetic background without electric or magnetic self-consistency.

classmethod doc_examples()

Create and initialize a Vlasov test-particle model:

from struphy.models import Vlasov

model = Vlasov()
model.kinetic_ions.var

classmethod doc_use_cases()

This model is appropriate for:

• test-particle motion in prescribed magnetic fields

• verification of the Boris-like VxB and PushEta splitting

• reduced kinetic transport studies without field feedback

classmethod doc_cannot_be_used_for()

This model is not suitable for:

• self-consistent electrostatic or electromagnetic coupling

• collisional kinetic dynamics

• guiding-center reduction studies

• fluid or MHD-scale closures

class struphy.models.GuidingCenter(base_units: BaseUnits = BaseUnits(x=1.0, B=1.0, n=1.0, kBT=None), charge_number: int = 1, mass_number: float = 1.0, epsilon: float | None = None)

Bases: StruphyModel

Guiding-center equation for a single species in a static background magnetic field.

Parameters:

• base_units (BaseUnits) – Base units for normalization (default: BaseUnits())

• charge_number (int) – Charge number (in units of the positive elementary charge) of the species (default: 1)

• mass_number (float) – Mass number (in units of Proton mass) of the species (default: 1.0)

• epsilon (float, optional) – Normalized cyclotron period: 1 / (cyclotron frequency × time unit). If None, computed from units and charge/mass numbers.

classmethod doc_pde()

PDEs solved by model:

Guiding-center Vlasov equation:

\[\frac{\partial f}{\partial t} + \left[ v_\parallel \frac{\mathbf{B}^*}{B^*_\parallel} + \frac{\mathbf{E}^* \times \mathbf{b}_0}{B^*_\parallel} \right] \cdot \frac{\partial f}{\partial \mathbf{X}} + \left[ \frac{1}{\epsilon} \frac{\mathbf{B}^*}{B^*_\parallel} \cdot \mathbf{E}^* \right] \cdot \frac{\partial f}{\partial v_\parallel} = 0\]

where \(f(\mathbf{X}, v_\parallel, \mu, t)\) is the guiding center distribution and

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

classmethod doc_normalization()

The reference speed is Alfvénic and the key small parameter is the normalized cyclotron time:

\[\hat v = \hat v_A,\qquad \varepsilon = 1/(\hat\Omega_c \hat t).\]

classmethod doc_scalar_quantities()

The following scalars are tracked during simulation:

• Parallel kinetic energy: en_fv

• Magnetic-moment energy contribution: en_fB

• Total particle energy: en_tot

• Lost markers: n_lost_particles

classmethod doc_discretization()

Time integration is performed by the following propagators (in sequence):

1. PushGuidingCenterBxEstar

2. PushGuidingCenterParallel

classmethod doc_long_description()

GuidingCenter is the reduced kinetic model used when fast gyromotion is averaged out and only guiding-center motion needs to be resolved. It evolves particles in a fixed magnetic background without any self-consistent field solve.

classmethod doc_examples()

Create and initialize a guiding-center model:

from struphy.models import GuidingCenter

model = GuidingCenter()
model.kinetic_ions.var

classmethod doc_use_cases()

This model is appropriate for:

• test-particle guiding-center orbit calculations

• strongly magnetized plasmas with scale separation from gyromotion

• verification of 5D guiding-center pushers and diagnostics

classmethod doc_cannot_be_used_for()

This model is not suitable for:

• self-consistent electric or magnetic field evolution

• full-orbit particle dynamics with resolved gyrophase

• collisional transport or source terms not present in the equation

• fluid closures or MHD force balance