Vlasov Family#

Models for the Vlasov equation of one kinetic species in 6D phase space (full orbits). Each model extends its parent by new terms of the equations and inherits its species, propagators and scalars (see Anatomy of a model).

Vlasov

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)[source]#

Bases: StruphyModel

Vlasov equation for a single species in a static background magnetic field.

The models of the Vlasov family inherit from this model: they call super().__init__() and extend the hooks of _setup() (species, propagators, variable links, scalars).

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()[source]#

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()[source]#: The characteristic speed is the cyclotron scale

\[\hat v = \hat\Omega_c \hat x.\]

classmethod doc_scalar_quantities()[source]#

The following scalars are tracked during simulation:

Particle kinetic energy: en_f

classmethod doc_long_description()[source]#: 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()[source]#

Create and initialize a Vlasov test-particle model:

from struphy.models import Vlasov

model = Vlasov()
model.kinetic_ions.var

classmethod doc_use_cases()[source]#

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()[source]#

This model is not suitable for:

self-consistent electrostatic or electromagnetic coupling
collisional kinetic dynamics
guiding-center reduction studies
fluid or MHD-scale closures

classmethod doc_discretization() → str#

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

push_vxb: PushVxB
push_eta: PushEta

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)[source]#

Bases: Vlasov

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

Extends Vlasov by the self-consistent electric field (Ampère’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 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.

classmethod doc_pde()[source]#

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()[source]#

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()[source]#

The following scalars are tracked during simulation:

Electric field energy en_E: \(\frac{1}{2} \int |\mathbf{E}|^2 \, \mathrm{d}V\)
Kinetic energy en_f: \(\frac{\alpha^2}{2N} \sum_p w_p |\mathbf{v}_p|^2\)
Total energy en_tot: en_E + en_f

classmethod doc_long_description()[source]#

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()[source]#

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["en_E"]
model.scalars["en_f"]
model.scalars["en_tot"]

classmethod doc_use_cases()[source]#

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()[source]#

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)

classmethod doc_discretization() → str#

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

push_eta: PushEta
push_vxb: PushVxB
coupling_va: VlasovAmpereCoupling

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)[source]#

Bases: VlasovAmpereOneSpecies

Vlasov-Maxwell equations for one kinetic species.

Extends VlasovAmpereOneSpecies by the magnetic field (Faraday’s law and the curl of the magnetic field in Ampère’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 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()[source]#

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()[source]#

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()[source]#

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_long_description()[source]#: 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()[source]#

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()[source]#

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()[source]#

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

classmethod doc_discretization() → str#

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

maxwell: MaxwellWeakAmpere
push_eta: PushEta
push_vxb: PushVxB
coupling_va: VlasovAmpereCoupling

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)[source]#

Bases: Vlasov

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. It inherits the particle propagators of Vlasov (with the Lorentz force of the total magnetic field \(\mathbf B_0 + \mathbf B\)).

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()[source]#

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()[source]#

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()[source]#

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_long_description()[source]#: 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()[source]#

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()[source]#

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()[source]#

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

classmethod doc_discretization() → str#

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

coupling: VlasovMaxwellCoupling
push_eta: PushEta
push_vxb: PushVxB

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)[source]#

Bases: Vlasov

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. Extends Vlasov by the electrostatic force and the Poisson equation.

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.

classmethod doc_pde()[source]#

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()[source]#

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()[source]#

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_f \(= \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()[source]#

Time integration is performed by the following propagators (in sequence). The Poisson solve comes right after the drift, such that with Strang splitting both velocity kicks use the potential of the current marker positions (drift-kick-kick-drift). Since a Strang step ends with the half drift, the saved potential and en_phi belong to the marker positions half a step earlier.

push_eta: PushEta
poisson: AdiabaticPoissonSolve
push_v_E: PushVinForceField
push_vxb: PushVxB

classmethod doc_long_description()[source]#

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()[source]#

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()[source]#

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()[source]#

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)[source]#

Bases: PauliMixin, 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)). Extends VlasovPoissonAdiabatic by the magnetic-moment force; for guiding-center orbits without the Poisson equation see VlasovPauli.

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.

classmethod doc_pde()[source]#

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.

classmethod doc_normalization()[source]#

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()[source]#

The following scalars are tracked during simulation (the sum of the energies 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), 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_f \(= \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_f + 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_f)
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()[source]#

Time integration is performed by the following propagators (in sequence). With Strang splitting, a time step is a drift by dt/2, the Poisson solve, a kick by dt/2 with the modified field E† = -grad(phi + epsilon mu |B0|), the rotation by v x B0 over the full dt, a second kick by dt/2 and a drift by dt/2. With PushVxB(algo="implicit") (Crank-Nicolson, i.e. the Cayley rotation) this is the Boris scheme of the classical Pauli particle (BAP2 in Xiao and Qin 2021), which gives the E x B and grad-B drifts for any dt / epsilon. The exact rotation (algo="analytic") and Lie-Trotter splitting do not; they need dt < epsilon.

The markers are loaded with PauliParticles6D: the magnetic moment of each marker is computed from its drawn perpendicular velocity, which is then set to zero.

push_eta: PushEta
poisson: AdiabaticPoissonSolve
push_v_E: PushVinForceField
push_vxb: PushVxB

classmethod doc_long_description()[source]#

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()[source]#

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)

# 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()[source]#

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, …)
verification of guiding-center pushers and models in the electrostatic limit

For guiding-center orbits in the equilibrium magnetic field (without the Poisson equation) use VlasovPauli.

classmethod doc_cannot_be_used_for()[source]#

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

class struphy.models.VlasovPauli(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)[source]#

Bases: PauliMixin, Vlasov

Guiding-center orbits in a static magnetic field, computed from classical Pauli particles.

The particles are full-orbit (6D) markers with a magnetic moment (classical Pauli particles, Xiao and Qin 2021), loaded on the slow manifold without gyration. Extends Vlasov by the magnetic-moment (mirror) force \(-\mu \nabla |\mathbf B_0|\). The markers follow guiding-center orbits with full-orbit pushers, with time steps larger than the cyclotron period. With the self-consistent electrostatic potential and adiabatic electrons, see VlasovPoissonAdiabaticPauli.

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 species (default: 1).
mass_number (float) – Mass number (in units of the proton mass) of the 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()[source]#

PDEs solved by model:

Vlasov equation for classical Pauli particles in a static magnetic field \(\mathbf B_0\), 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} \mathbf{v} \times \mathbf{B}_0 - \mu \nabla |\mathbf{B}_0| \right) \cdot \frac{\partial f}{\partial \mathbf{v}} = 0\,,\]

with \(f\) concentrated on the slow manifold \(\mathbf v \times \mathbf b_0 \approx 0\).

Conserved quantities: the number of particles and the energy

\[\frac{\textrm d}{\textrm d t} \int_0^\infty \int_{\mathbb{R}^3} \int_\Omega \left( \frac{1}{2} |\mathbf v|^2 + \mu |\mathbf B_0| \right) f \, \textrm d \mathbf x \, \textrm d^3 \mathbf v \, \textrm d \mu = 0\,,\]

for periodic boundary conditions or particles that do not leave the domain.

classmethod doc_normalization()[source]#

The reference speed is the thermal speed,

\[\hat v = \sqrt{k_B \hat T / m}\,,\qquad \hat \mu = \frac{m \hat v^2}{\hat B}\,,\]

with the dimensionless parameter \(\varepsilon = 1/(\hat\Omega_c\hat t)\). The normalized magnetic moment is \(\mu = v_\perp^2 / (2 |\mathbf B_0|)\), as in the guiding-center models.

classmethod doc_scalar_quantities()[source]#

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

Kinetic energy: en_f \(= \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_f + en_mu
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()[source]#

Time integration is performed by the following propagators (in sequence). With Strang splitting, a time step is a drift by dt/2, a kick by dt/2 with the magnetic-moment force -mu grad |B0|, the rotation by v x B0 over the full dt, a second kick by dt/2 and a drift by dt/2. With PushVxB(algo="implicit") (Crank-Nicolson, i.e. the Cayley rotation) this is the Boris scheme of the classical Pauli particle (BAP2 in Xiao and Qin 2021), which gives the grad-B and curvature drifts for any dt / epsilon. The exact rotation (algo="analytic") and Lie-Trotter splitting do not; they need dt < epsilon.

The markers are loaded with PauliParticles6D: the magnetic moment of each marker is computed from its drawn perpendicular velocity, which is then set to zero.

push_eta: PushEta
push_v_mu: PushVinForceField
push_vxb: PushVxB

classmethod doc_long_description()[source]#: Guiding-center orbits 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 - \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.

classmethod doc_examples()[source]#

Guiding-center orbits in the equilibrium magnetic field:

from struphy.models import VlasovPauli

model = VlasovPauli(epsilon=0.01)
model.kinetic_ions.var

# 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()[source]#

This model is appropriate for:

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
verification of guiding-center pushers

classmethod doc_cannot_be_used_for()[source]#

This model is not suitable for:

self-consistent fields (use VlasovPoissonAdiabaticPauli for the electrostatic case)
finite-Larmor-radius effects (use Vlasov)
magnetic fields without scale separation, where the magnetic moment is not an adiabatic invariant
time-dependent magnetic fields

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)[source]#

Bases: Vlasov

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

Extends Vlasov, written for the perturbation \(\tilde f\) (delta-f markers), by the force of the background electric field and the coupling of the perturbed electric field to the background (linearized Ampère’s law and the source term of the linearized Vlasov equation).

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.

classmethod doc_pde()[source]#

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()[source]#

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()[source]#

The following scalars are tracked during simulation:

Electric field energy: en_E
Perturbation particle energy: en_f
Total energy: en_tot

classmethod doc_long_description()[source]#: 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()[source]#

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()[source]#

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()[source]#

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

classmethod doc_discretization() → str#

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

push_eta: PushEta
push_vinE: PushVinForceField
coupling_Eweights: EfieldWeightsCoupling
push_vxb: PushVxB

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)[source]#

Bases: LinearVlasovAmpereOneSpecies

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

Extends LinearVlasovAmpereOneSpecies by the magnetic field perturbation (Faraday’s law and the curl of the magnetic field in Ampère’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 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.

classmethod doc_pde()[source]#

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()[source]#

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()[source]#

The following scalars are tracked during simulation:

Electric field energy: en_E
Magnetic field energy: en_B
Perturbation particle energy: en_f
Total energy: en_tot

classmethod doc_long_description()[source]#: 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()[source]#

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()[source]#

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()[source]#

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

classmethod doc_discretization() → str#

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

push_eta: PushEta
push_vinE: PushVinForceField
coupling_Eweights: EfieldWeightsCoupling
push_vxb: PushVxB
maxwell: MaxwellWeakAmpere

Guiding-Center Family#

Models for the guiding-center (drift-kinetic) equation of one kinetic species in 5D phase space. Each model extends its parent by new terms of the equations and inherits its species, propagators and scalars (see Anatomy of a model).

GuidingCenter
- DriftKineticElectrostaticAdiabatic
- DriftKineticMaxwell

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)[source]#

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()[source]#

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()[source]#: 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()[source]#

The following scalars are tracked during simulation:

Parallel kinetic energy: en_f
Magnetic-moment energy: en_mu
Total particle energy: en_tot
Lost markers: n_lost_particles

classmethod doc_long_description()[source]#: 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()[source]#

Create and initialize a guiding-center model:

from struphy.models import GuidingCenter

model = GuidingCenter()
model.kinetic_ions.var

classmethod doc_use_cases()[source]#

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()[source]#

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

classmethod doc_discretization() → str#

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

push_bxe: PushGuidingCenterBxEstar
push_parallel: PushGuidingCenterParallel

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)[source]#

Bases: GuidingCenter

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. Extends GuidingCenter by the electrostatic potential (in \(\mathbf E^*\)) and the quasi-neutrality equation.

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()[source]#

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()[source]#

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()[source]#

The following scalars are tracked during simulation:

Integral of the averaged squared potential (ITG diagnostic): phi_integral
Field energy: en_phi (one half of the quadratic form of the quasi-neutrality operator, see quadratic_form())
Parallel kinetic energy: en_f
Magnetic-moment energy: en_mu
Total energy: en_tot = en_phi + en_f + en_mu
Lost markers: n_lost_particles

classmethod doc_discretization()[source]#

Time integration is performed by the following propagators (in sequence). The Poisson solve sits between the two pushers, which both move the guiding centers. With Strang splitting, bxe(dt/2), poisson, para(dt), poisson, bxe(dt/2), the parallel push and the closing bxe(dt/2) use a freshly solved potential; the two bxe(dt/2) at a step boundary act as one bxe(dt) with the potential solved before it. With Lie-Trotter splitting the parallel push uses a fresh potential and the E* x b push the one of the previous step. Since a time step ends with a push, the saved potential and en_phi belong to the marker positions before that last push.

push_bxe: PushGuidingCenterBxEstar
gc_poisson: GyrokineticPoissonSolve
push_parallel: PushGuidingCenterParallel

classmethod doc_long_description()[source]#: 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()[source]#

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()[source]#

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()[source]#

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)[source]#

Bases: GuidingCenter

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. Extends GuidingCenter by the electromagnetic fields.

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()[source]#

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()[source]#

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()[source]#

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; unlike in GuidingCenter it contains the magnetic-moment energy (there is no en_mu) and the coefficients of the coupling step
Total energy: en_tot (conserved by the coupling step up to the solver tolerance)
Lost markers: n_lost_particles

classmethod doc_discretization()[source]#

Time integration is performed by the following propagators (in sequence). Each term of the drift-kinetic equation is applied by exactly one propagator; the parallel electric force only by the coupling step (the model passes no potential and no electric field to push_parallel, and asserts this):

\(\dot v_\parallel \ni \frac{1}{\varepsilon} \mathbf b^* \cdot \mathbf E\): coupling step (push_vpar_with_efield);
\(\dot v_\parallel \ni -\mu\, \mathbf b^* \cdot \nabla |B|\) (mirror force): push_parallel;
\(\dot{\mathbf X} \ni v_\parallel \mathbf b^*\) (streaming, curvature drift): push_parallel; its current enters Ampère’s law in the coupling step;
\(\dot{\mathbf X} \ni \mathbf E \times \mathbf b_0 / B^*_\parallel\): push_bxe; its current enters Ampère’s law in the coupling step (option exb_current);
\(\dot{\mathbf X} \ni -\varepsilon \mu \nabla |B| \times \mathbf b_0 / B^*_\parallel\) (grad-\(B\) drift): push_bxe; its current is not included.

coupling: DriftKineticMaxwellCoupling
push_parallel: PushGuidingCenterParallel
push_bxe: PushGuidingCenterBxEstar

classmethod doc_long_description()[source]#: 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()[source]#

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()[source]#

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()[source]#

This model is not suitable for:

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)

LinearEuler Family#

Models for small perturbations of a static (zero-flow) fluid or plasma equilibrium. Each model extends its parent by new terms of the equations and inherits its species, propagators and scalars (see Anatomy of a model).

LinearEuler
- LinearMHD

class struphy.models.LinearEuler(base_units: BaseUnits = BaseUnits(x=1.0, B=1.0, n=1.0, kBT=1.0), mass_number: float = 1.0)[source]#

Bases: StruphyModel

Linearized compressible Euler equations around a static (zero-flow) equilibrium, for sound waves.

The models of the LinearEuler family inherit from this model: they call super().__init__() and extend the hooks of _setup() (species, propagators, variable links, scalars).

Parameters:

base_units (BaseUnits) – Base units for normalization (default: BaseUnits(kBT=1.0)).
mass_number (float) – Mass number (in units of Proton mass) of the fluid species (default: 1.0)

classmethod doc_pde()[source]#

PDEs solved by model:

Continuity (mass conservation):

\[\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} = 0\]

Energy (adiabatic process):

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

Here, \(\tilde \rho\), \(\tilde{\mathbf U}\), \(\tilde p\) are the perturbed density, velocity and pressure, \(\rho_0\), \(p_0\) the density and pressure of a static equilibrium (zero flow, \(\mathbf U_0 = 0\)), and \(\gamma = 5/3\). Force balance of the equilibrium requires \(\nabla p_0 = 0\) (the density \(\rho_0\) can vary); the term \(\tilde{\mathbf U} \cdot \nabla p_0\) of the energy equation is nevertheless kept, since it is part of the models of the family whose equilibria carry a pressure gradient (e.g. \(\nabla p_0 = \mathbf J_0 \times \mathbf B_0\) in LinearMHD). For constant \(\rho_0\), \(p_0\) the perturbations are sound waves with the speed \(c_\textnormal{s} = \sqrt{\gamma p_0 / \rho_0}\).

Conserved quantities (for \(\tilde{\mathbf U} \cdot \mathbf n = 0\) on the boundary, or periodic boundaries):

Mass perturbation:

\[\frac{\textnormal d}{\textnormal d t} \int \tilde \rho \, \textnormal d \mathbf x = 0\]

Perturbed energy, for an equilibrium with \(\nabla p_0 = 0\):

\[\frac{\textnormal d}{\textnormal d t} \int \left( \frac{\rho_0}{2} |\tilde{\mathbf U}|^2 + \frac{\tilde p^2}{2 \gamma p_0} \right) \textnormal d \mathbf x = 0\]

classmethod doc_normalization()[source]#

All velocities are normalized by the thermal velocity of the bulk species, defined by the unit of the temperature \(\hat{k_B T}\) (BaseUnits.kBT):

\[\hat{U} = \hat{v}_\mathrm{th} = \sqrt{\frac{\hat{k_B T}}{A m_\mathrm{H}}}\,, \qquad \hat{\rho} = A m_\mathrm{H} \hat{n}\,, \qquad \hat{p} = \hat{\rho}\, \hat{v}_\mathrm{th}^2 = \hat{n} \hat{k_B T}\,, \qquad \hat{t} = \frac{\hat{x}}{\hat{v}_\mathrm{th}}\,.\]

The equations then have no dimensionless parameters.

classmethod doc_scalar_quantities()[source]#

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\) (en_U)
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\)
Equilibrium internal energy: \(E_{p0} = \int \frac{p_0}{\gamma - 1} \, \mathrm{d}V\) (en_p_eq)
Total perturbed energy: \(E_{\mathrm{tot}} = E_U + E_{\mathrm{thermal}}\) (en_tot). The sound wave step conserves it up to the solver tolerances; the background pressure step changes it only for \(\nabla p_0 \neq 0\).

classmethod doc_long_description()[source]#

This model simulates small-amplitude perturbations of a compressible, inviscid fluid around a static equilibrium. It solves the linearized Euler equations for the density, velocity and pressure perturbations. Its only waves are sound waves; with a non-uniform equilibrium density they are refracted.

This is the parent model of the LinearEuler family: LinearMHD adds the magnetic field (shear-Alfvén waves and the background current), and inherits the species, propagators and scalars of this model.

The time step is split into the sound wave step, which conserves the perturbed energy \(\frac 12 \int \rho_0 |\tilde{\mathbf U}|^2 \,\textnormal d \mathbf x + \frac{1}{2\gamma} \int \tilde p^2 / p_0 \,\textnormal d \mathbf x\) exactly (Crank-Nicolson with the \(1/p_0\)-weighted pressure equation), and the exact-in-time step of the term \(- \tilde{\mathbf U} \cdot \nabla p_0\), which vanishes for an equilibrium in force balance. The two steps do not commute; they are integrated with symmetric (Strang) splitting, the sound wave step in the middle, such that the time step is second order in time (also with Time(split_algo="LieTrotter")).

The equilibrium is given by an MHD or fluid equilibrium (equil) providing \(n_0\) and \(p_0\); a magnetic field of the equilibrium is not used. With free boundary conditions the pressure term imposes \(\tilde p = 0\) weakly at the boundary; with essential boundary conditions on the velocity (Dirichlet in the de Rham sequence), \(\tilde{\mathbf U} \cdot \mathbf n = 0\) holds strongly (rigid wall).

classmethod doc_examples()[source]#

Create and initialize a linear Euler model:

from struphy import BaseUnits
from struphy.models import LinearEuler

model = LinearEuler(base_units=BaseUnits(kBT=1.0))

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

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

classmethod doc_use_cases()[source]#

This model is appropriate for studying linear sound waves in a static compressible fluid.

Typical use cases include:

verification of the sound wave dispersion relation and of eigenmodes in curvilinear geometries
sound waves in non-uniform equilibrium densities (refraction)
testing the fluid part of the linear MHD models of the family

classmethod doc_cannot_be_used_for()[source]#

This model is not suitable for:

nonlinear fluid dynamics with finite-amplitude perturbations (shocks)
equilibria with non-zero background flow
dissipative effects such as viscosity or heat conduction
magnetized plasmas (see LinearMHD)

classmethod doc_discretization() → str#

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

bg_pressure: BackgroundPressureGradient
sound: SoundWaveLinear

class struphy.models.LinearMHD(base_units: BaseUnits = BaseUnits(x=1.0, B=1.0, n=1.0, kBT=None), mass_number: float = 1.0, b_space: Literal['Hdiv', 'Hcurl'] = 'Hdiv', projected: bool = False)[source]#

Bases: LinearEuler

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

Extends LinearEuler by the magnetic field: the Lorentz force and the induction equation.

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)
b_space ({"Hdiv", "Hcurl"}) – Space of the magnetic field (default: “Hdiv”). "Hdiv": strong induction equation, implicit shear-Alfvén solve for the velocity (ShearAlfvenPropagator). "Hcurl": weak induction equation, implicit curl-curl solve for the magnetic field (ShearAlfvenCurlCurlPropagator).
projected (bool) – Only with b_space="Hcurl": the curl-curl term of the implicit shear-Alfvén step contains the projection of the Lorentz force, which makes this step exactly energy conserving at the price of inner mass matrix solves (default: False).

classmethod doc_pde()[source]#

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

Here, \(\tilde \rho\), \(\tilde{\mathbf U}\), \(\tilde p\), \(\tilde{\mathbf B}\) are the perturbed density, velocity, pressure and magnetic field, and \(\rho_0\), \(p_0\), \(\mathbf B_0\) those of a static MHD equilibrium; \(\gamma = 5/3\).

Conserved quantities (for \(\tilde{\mathbf U} \cdot \mathbf n = 0\) and \(\mathbf B_0 \cdot \mathbf n = 0\) on the boundary, or periodic boundaries):

Mass perturbation and divergence of the magnetic field:

\[\frac{\textnormal d}{\textnormal d t} \int \tilde \rho \, \textnormal d \mathbf x = 0\,, \qquad \frac{\partial}{\partial t} \nabla \cdot \tilde{\mathbf B} = 0\]

Perturbed energy, for an equilibrium with \(\nabla \times \mathbf B_0 = 0\) and \(\nabla p_0 = 0\):

\[\frac{\textnormal d}{\textnormal d t} \int \left( \frac{\rho_0}{2} |\tilde{\mathbf U}|^2 + \frac 12 |\tilde{\mathbf B}|^2 + \frac{\tilde p^2}{2 \gamma p_0} \right) \textnormal d \mathbf x = 0\]

For general equilibria and perturbations generated by a plasma displacement \(\boldsymbol \xi\) (\(\partial_t \boldsymbol \xi = \tilde{\mathbf U}\), \(\tilde{\mathbf B} = \nabla \times (\boldsymbol \xi \times \mathbf B_0)\), \(\tilde p = - \boldsymbol \xi \cdot \nabla p_0 - \gamma p_0 \nabla \cdot \boldsymbol \xi\), e.g. for a pure velocity perturbation at \(t = 0\)), the conserved energy is the sum of the kinetic energy and the potential energy \(\delta W(\boldsymbol \xi)\) (energy principle of ideal MHD).

classmethod doc_normalization()[source]#

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()[source]#

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). Its change by the shear-Alfvén step is zero (up to the solver tolerances) with b_space="Hdiv" and with b_space="Hcurl", projected=True; with b_space="Hcurl", projected=False this step damps it numerically (the damping vanishes with the resolution along the background field). Its change by the sound wave step is zero too (up to the solver tolerances); the background current and background pressure steps change it only for \(\nabla \times \mathbf B_0 \neq 0\) or \(\nabla p_0 \neq 0\).
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()[source]#

Time integration is performed by the following propagators (in sequence). The shear-Alfvén step shear_alf is ShearAlfvenPropagator for b_space="Hdiv" (listed below) and ShearAlfvenCurlCurlPropagator for b_space="Hcurl" (at the end):

shear_alf: ShearAlfvenPropagator
bg_current: BackgroundCurrent
bg_pressure: BackgroundPressureGradient
sound: SoundWaveLinear

classmethod doc_long_description()[source]#

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

Two discretizations of the magnetic field are available through the flag b_space. With "Hdiv" (default), \(\tilde{\mathbf B}\) is a 2-form and the induction equation holds strongly. With "Hcurl", \(\tilde{\mathbf B}\) is a 1-form, the induction equation is weak and the implicit shear-Alfvén step solves a field-aligned curl-curl system for the magnetic field, preconditioned by Hiptmair-Xu, which allows large time steps; the force \((\nabla \times \mathbf B_0) \times \tilde{\mathbf B}\) of the background current step then enters through the mixed mass matrix M21J. This shear-Alfvén step damps the waves numerically, with a damping that vanishes quickly with the resolution along the background field; projected=True gives the exactly energy conserving version, at the price of an inner mass matrix solve in every application of the Schur operator.

Boundary conditions with b_space="Hcurl": essential boundary conditions in the de Rham sequence set the tangential components of the magnetic field to zero 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 sound wave step (natural boundary condition of the pressure term).

The species, the sound wave step and the background pressure step are inherited from LinearEuler. The compressional part of the time step thus consists of three steps: the sound wave step, which conserves the energy \(\frac 12 \int \rho_0 |\tilde{\mathbf U}|^2 \,\textnormal d \mathbf x + \frac{1}{2\gamma} \int \tilde p^2 / p_0 \,\textnormal d \mathbf x\) exactly (Crank-Nicolson with the \(1/p_0\)-weighted pressure equation), and the steps of the two terms that change the perturbed energy for general equilibria, \((\nabla \times \mathbf B_0) \times \tilde{\mathbf B}\) in the momentum equation and \(- \tilde{\mathbf U} \cdot \nabla p_0\) in the pressure equation, which are solved exactly in time. The three steps do not commute; they form one block of the splitting with symmetric (Strang) splitting, the sound wave step in the middle, so that the compressional part is second order in time as the former unsplit magnetosonic step (with Lie-Trotter sequencing of the three steps the frequencies of modes coupled to the compressional dynamics, e.g. TAEs in a torus, are wrong for large time steps). For an equilibrium with \(\nabla p_0 = 0\) and \(\nabla \times \mathbf B_0 = 0\), the scheme is thus energy conserving up to the solver tolerances (with b_space="Hdiv" or projected=True).

classmethod doc_examples()[source]#

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"]

# magnetic field in H(curl), exactly energy conserving shear-Alfvén step
model_cc = LinearMHD(b_space="Hcurl", projected=True)

classmethod doc_use_cases()[source]#

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
large time steps with the curl-curl discretization (b_space="Hcurl")

classmethod doc_cannot_be_used_for()[source]#

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 (b_space="Hcurl")

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)[source]#

Bases: LinearMHD

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

Extends LinearMHD (magnetic field in H(div)) by energetic ions, whose density and current enter the MHD momentum equation.

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()[source]#

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} + \nabla \cdot (p_0 \tilde{\mathbf{U}}) + (\gamma - 1) 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\), \(\gamma = 5/3\), the equations of LinearMHD are extended by the energetic ions; the terms in blue couple them to the MHD equations.

classmethod doc_normalization()[source]#

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()[source]#

The following scalars are tracked during simulation:

The scalars of LinearMHD: en_U, en_p, en_thermal, en_p_eq, en_B, en_B_eq, en_B_tot
Energetic-particle kinetic energy: en_f
Total energy: en_tot = en_U + en_thermal + en_B + en_f
Lost particles: n_lost_particles

classmethod doc_long_description()[source]#: 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()[source]#

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()[source]#

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()[source]#

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

classmethod doc_discretization()#

Time integration is performed by the following propagators (in sequence). The shear-Alfvén step shear_alf is ShearAlfvenPropagator for b_space="Hdiv" (listed below) and ShearAlfvenCurlCurlPropagator for b_space="Hcurl" (at the end):

couple_dens: CurrentCoupling6DDensity
shear_alf: ShearAlfvenPropagator
couple_curr: CurrentCoupling6DCurrent
push_eta: PushEta
push_vxb: PushVxB
bg_current: BackgroundCurrent
bg_pressure: BackgroundPressureGradient
sound: SoundWaveLinear

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)[source]#

Bases: LinearMHD

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

Extends LinearMHD (magnetic field in H(div)) by energetic ions, whose perpendicular pressure tensor enters the MHD momentum equation.

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()[source]#

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()[source]#: 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()[source]#

The following scalars are tracked during simulation:

The scalars of LinearMHD: en_U, en_p, en_thermal, en_p_eq, en_B, en_B_eq, en_B_tot
Energetic-particle kinetic energy: en_f
Total energy: en_tot = en_U + en_thermal + en_B + en_f
Lost particles: n_lost_particles

classmethod doc_long_description()[source]#: 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()[source]#

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()[source]#

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()[source]#

This model is not suitable for:

current-coupling closure studies
nonlinear hybrid turbulence
dissipative/resistive MHD
fully kinetic treatment of the bulk plasma

classmethod doc_discretization()#

Time integration is performed by the following propagators (in sequence). The shear-Alfvén step shear_alf is ShearAlfvenPropagator for b_space="Hdiv" (listed below) and ShearAlfvenCurlCurlPropagator for b_space="Hcurl" (at the end):

push_eta_pc: PushEtaPC
push_vxb: PushVxB
pc6d: PressureCoupling6D
shear_alf: ShearAlfvenPropagator
bg_current: BackgroundCurrent
bg_pressure: BackgroundPressureGradient
sound: SoundWaveLinear

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)[source]#

Bases: LinearMHD

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

Extends LinearMHD (magnetic field in H(div)) by energetic ions in the guiding-center approximation, whose current and magnetization enter the MHD momentum equation. The shear-Alfvén step shear_alf of LinearMHD is replaced by ShearAlfvenCurrentCoupling5D, which includes the coupling to the guiding-center density.

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()[source]#

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()[source]#: 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()[source]#

The following scalars are tracked during simulation:

The scalars of LinearMHD: en_U, en_p, en_thermal, en_p_eq, en_B, en_B_eq, en_B_tot
Parallel energetic-particle energy: en_fv
Magnetic-moment energetic-particle energy: en_fB
Total energy: en_tot = en_U + en_thermal + en_B + en_fv + en_fB
Lost particles: n_lost_particles

classmethod doc_discretization()[source]#

Time integration is performed by the following propagators (in sequence). 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 push_bxe and algo="discrete_gradient_1st_order_newton" (Itoh-Abe with Newton iterations) for push_parallel. The explicit Runge-Kutta methods (algo="explicit") are cheaper but not energy conserving.

push_bxe: PushGuidingCenterBxEstar
push_parallel: PushGuidingCenterParallel
shear_alf: ShearAlfvenCurrentCoupling5D
bg_current: BackgroundCurrent
bg_pressure: BackgroundPressureGradient
sound: SoundWaveLinear
cc5d_density: CurrentCoupling5DDensity
cc5d_gradb: CurrentCoupling5DGradB
cc5d_curlb: CurrentCoupling5DCurlb

classmethod doc_long_description()[source]#: 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()[source]#

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()[source]#

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()[source]#

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

LinearEulerWithQ Family#

The equations of the LinearEuler family written with the square root of the pressure, \(q = \sqrt p\), as thermodynamic variable (pressure perturbation \(\tilde p = 2 q_0 \tilde q\)). The magnetic field of LinearMHDWithQ is the one of LinearMHD.

LinearEulerWithQ
- LinearMHDWithQ
  - PerturbativeMHDWithQ
- PerturbativeEulerWithQ

The perturbative models add the terms that are nonlinear in the perturbations, each in its own energy conserving step, to the steps of their linear parent.

class struphy.models.LinearEulerWithQ(base_units: BaseUnits = BaseUnits(x=1.0, B=1.0, n=1.0, kBT=1.0), mass_number: float = 1.0)[source]#

Bases: StruphyModel

Linearized compressible Euler equations around a static (zero-flow) equilibrium, for sound waves, with the square root of the pressure, \(q = \sqrt p\), as thermodynamic variable.

The models of the LinearEulerWithQ family inherit from this model: they call super().__init__() and extend the hooks of _setup() (species, propagators, variable links, scalars).

Parameters:

base_units (BaseUnits) – Base units for normalization (default: BaseUnits(kBT=1.0)).
mass_number (float) – Mass number (in units of Proton mass) of the fluid species (default: 1.0)

classmethod doc_pde()[source]#

PDEs solved by model:

Continuity (mass conservation):

\[\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 (2 q_0 \tilde q) = 0\]

Square root of the pressure (adiabatic process):

\[\frac{\partial \tilde{q}}{\partial t} + \tilde{\mathbf{U}} \cdot \nabla q_0 + \frac{\gamma}{2} q_0 \nabla \cdot \tilde{\mathbf{U}} = 0\,.\]

Here, \(\tilde \rho\), \(\tilde{\mathbf U}\), \(\tilde q\) are the perturbed density, velocity and square root of the pressure, \(\rho_0\), \(q_0 = \sqrt{p_0}\) those of a static equilibrium (zero flow, \(\mathbf U_0 = 0\)), and \(\gamma = 5/3\). The equation of \(q\) is the linearization of \(\partial_t q + \mathbf U \cdot \nabla q + \frac{\gamma}{2} q \nabla \cdot \mathbf U = 0\), which follows from the adiabatic pressure equation with \(p = q^2\); the pressure perturbation is \(\tilde p = 2 q_0 \tilde q\), and the equations are those of LinearEuler written in \(\tilde q\). Force balance of the equilibrium requires \(\nabla q_0 = 0\); the term \(\tilde{\mathbf U} \cdot \nabla q_0\) is nevertheless kept, since it is part of the models of the family whose equilibria carry a pressure gradient (e.g. LinearMHDWithQ). For constant \(\rho_0\), \(q_0\) the perturbations are sound waves with the speed \(c_\textnormal{s} = \sqrt{\gamma q_0^2 / \rho_0}\).

Conserved quantities (for \(\tilde{\mathbf U} \cdot \mathbf n = 0\) on the boundary, or periodic boundaries):

Mass perturbation:

\[\frac{\textnormal d}{\textnormal d t} \int \tilde \rho \, \textnormal d \mathbf x = 0\]

Perturbed energy, for an equilibrium with \(\nabla q_0 = 0\):

\[\frac{\textnormal d}{\textnormal d t} \int \left( \frac{\rho_0}{2} |\tilde{\mathbf U}|^2 + \frac{2}{\gamma} \tilde q^2 \right) \textnormal d \mathbf x = 0\]

classmethod doc_normalization()[source]#

All velocities are normalized by the thermal velocity of the bulk species, defined by the unit of the temperature \(\hat{k_B T}\) (BaseUnits.kBT):

\[\hat{U} = \hat{v}_\mathrm{th} = \sqrt{\frac{\hat{k_B T}}{A m_\mathrm{H}}}\,, \qquad \hat{\rho} = A m_\mathrm{H} \hat{n}\,, \qquad \hat{q} = \sqrt{\hat{p}} = \sqrt{\hat{\rho}}\, \hat{v}_\mathrm{th}\,, \qquad \hat{t} = \frac{\hat{x}}{\hat{v}_\mathrm{th}}\,.\]

The equations then have no dimensionless parameters.

classmethod doc_scalar_quantities()[source]#

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\) (en_U)
Thermal perturbation energy: \(E_{\mathrm{thermal}} = \frac{2}{\gamma}\int \tilde{q}^2 \, \mathrm{d}V\) (en_thermal), with \(\gamma = 5/3\); it equals \(\frac 12 \int \tilde p^2 / (\gamma p_0)\) for \(\tilde p = 2 q_0 \tilde q\)
Equilibrium internal energy: \(E_{p0} = \int \frac{p_0}{\gamma - 1} \, \mathrm{d}V\) (en_p_eq)
Total perturbed energy: \(E_{\mathrm{tot}} = E_U + E_{\mathrm{thermal}}\) (en_tot). The sound wave step conserves it up to the solver tolerances; the background step changes it only for \(\nabla q_0 \neq 0\).

classmethod doc_long_description()[source]#

This model simulates small-amplitude perturbations of a compressible, inviscid fluid around a static equilibrium. It solves the linearized Euler equations for the perturbations of the density, the velocity and the square root of the pressure, \(q = \sqrt p\). Its only waves are sound waves; with a non-uniform equilibrium density they are refracted. It solves the same equations as LinearEuler, written in \(\tilde q = \tilde p / (2 q_0)\); it is the parent of the family with the variable \(q\), whose nonlinear (perturbative) models keep the terms quadratic in \(\tilde q\) in the energy.

LinearMHDWithQ adds the magnetic field (shear-Alfvén waves and the background current), and inherits the species, propagators and scalars of this model.

The time step is split into the sound wave step, which conserves the perturbed energy \(\frac 12 \int \rho_0 |\tilde{\mathbf U}|^2 \,\textnormal d \mathbf x + \frac{2}{\gamma} \int \tilde q^2 \,\textnormal d \mathbf x\) exactly (Crank-Nicolson; the force is the exact transpose of the operator of the \(\tilde q\) equation, so the discrete spectrum is purely imaginary and the step does not damp or amplify the waves), and the exact-in-time step of the term \(- \tilde{\mathbf U} \cdot \nabla q_0\), which vanishes for an equilibrium in force balance. The two steps do not commute; they are integrated with symmetric (Strang) splitting, the sound wave step in the middle, such that the time step is second order in time (also with Time(split_algo="LieTrotter")).

The equilibrium is given by an MHD or fluid equilibrium (equil) providing \(n_0\) and \(p_0\); a magnetic field of the equilibrium is not used. With free boundary conditions the pressure term imposes \(\tilde q = 0\) weakly at the boundary; with essential boundary conditions on the velocity (Dirichlet in the de Rham sequence), \(\tilde{\mathbf U} \cdot \mathbf n = 0\) holds strongly (rigid wall).

classmethod doc_examples()[source]#

Create and initialize a linear Euler model in \(q\):

from struphy import BaseUnits
from struphy.models import LinearEulerWithQ

model = LinearEulerWithQ(base_units=BaseUnits(kBT=1.0))

# Access fields
model.mhd.density
model.mhd.velocity
model.mhd.sqrt_p

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

classmethod doc_use_cases()[source]#

This model is appropriate for studying linear sound waves in a static compressible fluid with the variable \(q = \sqrt p\).

Typical use cases include:

verification of the sound wave dispersion relation and of eigenmodes in curvilinear geometries
sound waves in non-uniform equilibrium densities (refraction)
the linear limit of the perturbative models in \(q\) and testing the fluid part of LinearMHDWithQ

classmethod doc_cannot_be_used_for()[source]#

This model is not suitable for:

nonlinear fluid dynamics with finite-amplitude perturbations (shocks)
equilibria with non-zero background flow
dissipative effects such as viscosity or heat conduction
magnetized plasmas (see LinearMHDWithQ)

classmethod doc_discretization() → str#

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

bg_q: BackgroundQGradient
sound: SoundWaveLinearWithQ

class struphy.models.LinearMHDWithQ(base_units: BaseUnits = BaseUnits(x=1.0, B=1.0, n=1.0, kBT=None), mass_number: float = 1.0, b_space: Literal['Hdiv', 'Hcurl'] = 'Hdiv', projected: bool = False)[source]#

Bases: LinearEulerWithQ

Linear ideal MHD with zero-flow equilibrium, with the square root of the pressure, \(q = \sqrt p\), as thermodynamic variable.

Extends LinearEulerWithQ by the magnetic field: the Lorentz force and the induction equation. The magnetic field (species em_fields, propagators shear_alf and bg_current, scalars) is the one of LinearMHD.

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)
b_space ({"Hdiv", "Hcurl"}) – Space of the magnetic field (default: “Hdiv”). "Hdiv": strong induction equation, implicit shear-Alfvén solve for the velocity (ShearAlfvenPropagator). "Hcurl": weak induction equation, implicit curl-curl solve for the magnetic field (ShearAlfvenCurlCurlPropagator).
projected (bool) – Only with b_space="Hcurl": the curl-curl term of the implicit shear-Alfvén step contains the projection of the Lorentz force, which makes this step exactly energy conserving at the price of inner mass matrix solves (default: False).

classmethod doc_pde()[source]#

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 (2 q_0 \tilde q) = (\nabla \times \tilde{\mathbf{B}}) \times \mathbf{B}_0 + (\nabla \times \mathbf{B}_0) \times \tilde{\mathbf{B}}\]

Square root of the pressure (adiabatic process):

\[\frac{\partial \tilde{q}}{\partial t} + \tilde{\mathbf{U}} \cdot \nabla q_0 + \frac{\gamma}{2} q_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.\]

Here, \(\tilde \rho\), \(\tilde{\mathbf U}\), \(\tilde q\), \(\tilde{\mathbf B}\) are the perturbed density, velocity, square root of the pressure and magnetic field, and \(\rho_0\), \(q_0 = \sqrt{p_0}\), \(\mathbf B_0\) those of a static MHD equilibrium; \(\gamma = 5/3\). The pressure perturbation is \(\tilde p = 2 q_0 \tilde q\), and the equations are those of LinearMHD written in \(\tilde q\).

Conserved quantities (for \(\tilde{\mathbf U} \cdot \mathbf n = 0\) and \(\mathbf B_0 \cdot \mathbf n = 0\) on the boundary, or periodic boundaries):

Mass perturbation and divergence of the magnetic field:

\[\frac{\textnormal d}{\textnormal d t} \int \tilde \rho \, \textnormal d \mathbf x = 0\,, \qquad \frac{\partial}{\partial t} \nabla \cdot \tilde{\mathbf B} = 0\]

Perturbed energy, for an equilibrium with \(\nabla \times \mathbf B_0 = 0\) and \(\nabla q_0 = 0\):

\[\frac{\textnormal d}{\textnormal d t} \int \left( \frac{\rho_0}{2} |\tilde{\mathbf U}|^2 + \frac 12 |\tilde{\mathbf B}|^2 + \frac{2}{\gamma} \tilde q^2 \right) \textnormal d \mathbf x = 0\]

For general equilibria and perturbations generated by a plasma displacement \(\boldsymbol \xi\), the conserved energy is the sum of the kinetic energy and the potential energy \(\delta W(\boldsymbol \xi)\) (energy principle of ideal MHD), as for LinearMHD.

classmethod doc_normalization()[source]#

All velocities are normalized by the Alfvén velocity:

\[\hat{U} = \hat{v}_\mathrm{A} = \frac{\hat{B}_0}{\sqrt{\mu_0 \rho_0}}\,, \qquad \hat{q} = \sqrt{\hat p} = \frac{\hat B_0}{\sqrt{\mu_0}}\]

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

classmethod doc_scalar_quantities()[source]#

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\) (en_U)
Magnetic energy (perturbation): \(E_B = \frac{1}{2} \int \frac{|\tilde{\mathbf{B}}|^2}{\mu_0} \, \mathrm{d}V\) (en_B)
Thermal perturbation energy: \(E_{\mathrm{thermal}} = \frac{2}{\gamma}\int \tilde{q}^2 \, \mathrm{d}V\) (en_thermal), with \(\gamma = 5/3\)
Total perturbed energy: \(E_{\mathrm{tot}} = E_U + E_B + E_{\mathrm{thermal}}\) (en_tot). Its change by the shear-Alfvén step is zero (up to the solver tolerances) with b_space="Hdiv" and with b_space="Hcurl", projected=True; with b_space="Hcurl", projected=False this step damps it numerically (the damping vanishes with the resolution along the background field). Its change by the sound wave step is zero too (up to the solver tolerances); the background current and background \(q\) steps change it only for \(\nabla \times \mathbf B_0 \neq 0\) or \(\nabla q_0 \neq 0\).
Equilibrium magnetic energy: \(E_{B0} = \frac{1}{2} \int \frac{|\mathbf{B}_0|^2}{\mu_0} \, \mathrm{d}V\) (en_B_eq)
Equilibrium internal energy: \(E_{p0} = \int \frac{p_0}{\gamma - 1} \, \mathrm{d}V\) (en_p_eq)
Total magnetic energy: \(E_{B,\mathrm{tot}} = \frac{1}{2} \int \frac{|\mathbf{B}_0 + \tilde{\mathbf{B}}|^2}{\mu_0} \, \mathrm{d}V\) (en_B_tot)

classmethod doc_discretization()[source]#

Time integration is performed by the following propagators (in sequence). The shear-Alfvén step shear_alf is ShearAlfvenPropagator for b_space="Hdiv" (listed below) and ShearAlfvenCurlCurlPropagator for b_space="Hcurl" (at the end):

classmethod doc_long_description()[source]#

This model simulates small-amplitude perturbations in a magnetized plasma with a static equilibrium magnetic field \(\mathbf{B}_0\) and zero background flow. It solves the linearized ideal MHD equations of LinearMHD with the square root of the pressure, \(q = \sqrt p\), instead of the pressure; it is the linear limit of the perturbative MHD model in \(q\), and supports the same waves (shear Alfvén, fast and slow magnetosonic).

The magnetic field and its two discretizations (flag b_space: "Hdiv" with the strong induction equation, "Hcurl" with the weak induction equation and the field-aligned curl-curl solve, exactly energy conserving with projected=True) are those of LinearMHD, shared through the utilities in struphy.models.utils. The species, the sound wave step and the background \(q\) step are inherited from LinearEulerWithQ.

The compressional part of the time step consists of three steps: the sound wave step, which conserves the energy \(\frac 12 \int \rho_0 |\tilde{\mathbf U}|^2 \,\textnormal d \mathbf x + \frac{2}{\gamma} \int \tilde q^2 \,\textnormal d \mathbf x\) exactly (Crank-Nicolson), and the steps of the two terms that change the perturbed energy for general equilibria, \((\nabla \times \mathbf B_0) \times \tilde{\mathbf B}\) in the momentum equation and \(- \tilde{\mathbf U} \cdot \nabla q_0\) in the equation of \(\tilde q\), which are solved exactly in time. The three steps form one block of the splitting with symmetric (Strang) splitting, the sound wave step in the middle, so that the compressional part is second order in time. For an equilibrium with \(\nabla p_0 = 0\) and \(\nabla \times \mathbf B_0 = 0\), the scheme is thus energy conserving up to the solver tolerances (with b_space="Hdiv" or projected=True).

Boundary conditions with b_space="Hcurl": free boundary conditions are recommended, as for LinearMHD.

classmethod doc_examples()[source]#

Create and initialize a linear MHD model in \(q\):

from struphy.models import LinearMHDWithQ

model = LinearMHDWithQ()

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

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

# magnetic field in H(curl), exactly energy conserving shear-Alfvén step
model_cc = LinearMHDWithQ(b_space="Hcurl", projected=True)

classmethod doc_use_cases()[source]#

This model is appropriate for studying linear wave propagation in a static magnetized plasma with the variable \(q = \sqrt p\), including Alfvén and magnetosonic dynamics around a prescribed equilibrium.

Typical use cases include:

verification of linear MHD wave dispersion and mode structure
the linear limit of the perturbative MHD model in \(q\)
comparison with LinearMHD (same equations in \(\tilde p\))

classmethod doc_cannot_be_used_for()[source]#

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 (b_space="Hcurl")

class struphy.models.PerturbativeEulerWithQ(base_units: BaseUnits = BaseUnits(x=1.0, B=1.0, n=1.0, kBT=1.0), mass_number: float = 1.0)[source]#

Bases: LinearEulerWithQ

Compressible Euler equations for the perturbations of a static (zero-flow) equilibrium, with all nonlinear terms, and the square root of the pressure, \(q = \sqrt p\), as thermodynamic variable.

Extends LinearEulerWithQ by the terms that are nonlinear in the perturbations: the inertia of the density perturbation, the advection of mass and momentum and the nonlinear sound wave terms.

Parameters:

base_units (BaseUnits) – Base units for normalization (default: BaseUnits(kBT=1.0)).
mass_number (float) – Mass number (in units of Proton mass) of the fluid species (default: 1.0)

classmethod doc_pde()[source]#

PDEs solved by model:

Continuity (mass conservation):

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

Momentum:

\[(\rho_0 + \tilde \rho) \left( \frac{\partial \tilde{\mathbf{U}}}{\partial t} + \tilde{\mathbf U} \cdot \nabla \tilde{\mathbf U} \right) + \nabla (2 q_0 \tilde q + \tilde q^2) = 0\]

Square root of the pressure (adiabatic process):

\[\frac{\partial \tilde{q}}{\partial t} + \tilde{\mathbf{U}} \cdot \nabla (q_0 + \tilde q) + \frac{\gamma}{2} (q_0 + \tilde q) \nabla \cdot \tilde{\mathbf{U}} = 0\,.\]

Here, \(\tilde \rho\), \(\tilde{\mathbf U}\), \(\tilde q\) are the perturbed density, velocity and square root of the pressure (not necessarily small), \(\rho_0\), \(q_0 = \sqrt{p_0}\) those of a static equilibrium (zero flow) with \(\nabla p_0 = 0\), and \(\gamma = 5/3\). These are the compressible Euler equations \(\partial_t \rho + \nabla \cdot (\rho \mathbf U) = 0\), \(\rho (\partial_t \mathbf U + \mathbf U \cdot \nabla \mathbf U) + \nabla q^2 = 0\), \(\partial_t q + \mathbf U \cdot \nabla q + \frac{\gamma}{2} q \nabla \cdot \mathbf U = 0\) for \(\rho = \rho_0 + \tilde \rho\), \(\mathbf U = \tilde{\mathbf U}\), \(q = q_0 + \tilde q\), with the equilibrium subtracted. Dropping the terms quadratic in the perturbations gives LinearEulerWithQ.

Conserved quantities (for \(\tilde{\mathbf U} \cdot \mathbf n = 0\) on the boundary, or periodic boundaries):

Mass perturbation:

\[\frac{\textnormal d}{\textnormal d t} \int \tilde \rho \, \textnormal d \mathbf x = 0\]

Energy of the perturbation (kinetic plus internal energy, the internal energy of the equilibrium subtracted), for an equilibrium with \(\nabla q_0 = 0\):

\[\frac{\textnormal d}{\textnormal d t} \int \left( \frac{\rho_0 + \tilde \rho}{2} |\tilde{\mathbf U}|^2 + \frac{2 q_0 \tilde q + \tilde q^2}{\gamma - 1} \right) \textnormal d \mathbf x = 0\]

classmethod doc_normalization()[source]#

All velocities are normalized by the thermal velocity of the bulk species, defined by the unit of the temperature \(\hat{k_B T}\) (BaseUnits.kBT):

\[\hat{U} = \hat{v}_\mathrm{th} = \sqrt{\frac{\hat{k_B T}}{A m_\mathrm{H}}}\,, \qquad \hat{\rho} = A m_\mathrm{H} \hat{n}\,, \qquad \hat{q} = \sqrt{\hat{p}} = \sqrt{\hat{\rho}}\, \hat{v}_\mathrm{th}\,, \qquad \hat{t} = \frac{\hat{x}}{\hat{v}_\mathrm{th}}\,.\]

The equations then have no dimensionless parameters.

classmethod doc_scalar_quantities()[source]#

The following scalars are tracked during simulation:

Kinetic energy: \(E_U = \frac{1}{2} \int (\rho_0 + \tilde \rho) |\tilde{\mathbf{U}}|^2 \, \mathrm{d}V\) (en_U), with the total density (projected \(\rho_0\))
Quadratic thermal energy of the linear sound wave step: \(\frac{2}{\gamma}\int \tilde{q}^2 \, \mathrm{d}V\) (en_thermal)
Quadratic thermal energy of the nonlinear sound wave step: \(\frac{1}{\gamma - 1}\int \tilde{q}^2 \, \mathrm{d}V\) (en_q2)
Internal energy of the perturbation: \(E_\mathrm{int} = \int \frac{2 q_0 \tilde q + \tilde q^2}{\gamma - 1} \, \mathrm{d}V = \int \frac{q^2 - q_0^2}{\gamma - 1} \, \mathrm{d}V\) (en_internal)
Equilibrium internal energy: \(E_{p0} = \int \frac{p_0}{\gamma - 1} \, \mathrm{d}V\) (en_p_eq)
Energy of the perturbation: \(E_{\mathrm{tot}} = E_U + E_\mathrm{int}\) (en_tot), conserved by the equations for \(\nabla q_0 = 0\). Each step of the splitting conserves its own energy exactly (the linear sound wave step \(E_U\) + en_thermal, the nonlinear one \(E_U\) + en_q2, the advection step \(E_U\)), so \(E_{\mathrm{tot}}\) changes by the splitting error only.

classmethod doc_long_description()[source]#

This model simulates finite-amplitude perturbations of a compressible, inviscid fluid around a static equilibrium (e.g. nonlinear sound waves up to the formation of shocks). It solves the full compressible Euler equations for the perturbations of the density, the velocity and the square root of the pressure, \(q = \sqrt p\), written as the linear system of LinearEulerWithQ plus the nonlinear terms.

The time step reuses the steps of the linear model and adds the steps of the nonlinear terms:

the steps of LinearEulerWithQ (background \(q\) step and linear sound wave step, one symmetric block), with the inertia of the total density \(\rho_0 + \tilde \rho^n\) instead of \(\rho_0\) and without the continuity equation;
nonlin_sound: the terms of the sound waves quadratic in \(\tilde q\), with the transport operator \(\mathcal S_{\tilde q}\) and its transpose, so that it conserves \(\frac 12 \int \rho |\tilde{\mathbf U}|^2 + \frac{1}{\gamma - 1} \int \tilde q^2\) exactly;
mom_adv: the continuity equation and the advection of the momentum, which conserves the mass and \(\frac 12 \int \rho |\tilde{\mathbf U}|^2\) exactly (discrete gradient, Picard iteration).

All steps are energy conserving, each for its own energy; the nonlinear steps have the weights of the nonlinear terms at \(t^n\) by default (Options(weights="frozen")) or at \(t^{n+1/2}\) (weights="midpoint", Picard iteration, second order in time). The inertia of every step is that of the total density at \(t^n\), which is not changed except in mom_adv. No numerical dissipation is added, so shocks are not captured; the model is meant for smooth solutions.

classmethod doc_examples()[source]#

Create and initialize a perturbative Euler model in \(q\):

from struphy import BaseUnits
from struphy.models import PerturbativeEulerWithQ

model = PerturbativeEulerWithQ(base_units=BaseUnits(kBT=1.0))

# Access fields
model.mhd.density
model.mhd.velocity
model.mhd.sqrt_p

# second order in time: midpoint weights of the nonlinear steps
props = model.propagators
props.nonlin_sound.options = props.nonlin_sound.Options(weights="midpoint")
props.mom_adv.options = props.mom_adv.Options(weights="midpoint")

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

classmethod doc_use_cases()[source]#

This model is appropriate for finite-amplitude perturbations of a static compressible fluid.

Typical use cases include:

nonlinear sound waves (steepening, before the formation of shocks)
pressureless and compressible advection of density and momentum
the nonlinear fluid part of PerturbativeMHDWithQ

classmethod doc_cannot_be_used_for()[source]#

This model is not suitable for:

shocks and other discontinuities (no numerical dissipation)
equilibria with non-zero background flow, or with a pressure gradient (no force balance without a magnetic field)
dissipative effects such as viscosity or heat conduction
magnetized plasmas (see PerturbativeMHDWithQ)

classmethod doc_discretization() → str#

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

bg_q: BackgroundQGradient
sound: SoundWaveLinearWithQ
nonlin_sound: SoundWaveNonlinearWithQ
mom_adv: MomentumAdvection

class struphy.models.PerturbativeMHDWithQ(base_units: BaseUnits = BaseUnits(x=1.0, B=1.0, n=1.0, kBT=None), mass_number: float = 1.0)[source]#

Bases: LinearMHDWithQ

Ideal MHD equations for the perturbations of a static (zero-flow) equilibrium, with all nonlinear terms, and the square root of the pressure, \(q = \sqrt p\), as thermodynamic variable.

Extends LinearMHDWithQ by the terms that are nonlinear in the perturbations: the inertia of the density perturbation, the advection of mass and momentum, the nonlinear sound wave terms and the nonlinear terms of the Lorentz force and of the induction equation. The magnetic field is in \(H(\textnormal{div})\).

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()[source]#

PDEs solved by model:

Continuity (mass conservation):

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

Momentum (Lorentz force):

\[(\rho_0 + \tilde \rho) \left( \frac{\partial \tilde{\mathbf{U}}}{\partial t} + \tilde{\mathbf U} \cdot \nabla \tilde{\mathbf U} \right) + \nabla (2 q_0 \tilde q + \tilde q^2) = (\nabla \times \tilde{\mathbf{B}}) \times (\mathbf{B}_0 + \tilde{\mathbf B}) + (\nabla \times \mathbf{B}_0) \times \tilde{\mathbf{B}}\]

Square root of the pressure (adiabatic process):

\[\frac{\partial \tilde{q}}{\partial t} + \tilde{\mathbf{U}} \cdot \nabla (q_0 + \tilde q) + \frac{\gamma}{2} (q_0 + \tilde q) \nabla \cdot \tilde{\mathbf{U}} = 0\]

Induction (Faraday’s law):

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

Here, \(\tilde \rho\), \(\tilde{\mathbf U}\), \(\tilde q\), \(\tilde{\mathbf B}\) are the perturbed density, velocity, square root of the pressure and magnetic field (not necessarily small), and \(\rho_0\), \(q_0 = \sqrt{p_0}\), \(\mathbf B_0\) those of a static MHD equilibrium, \(\nabla q_0^2 = (\nabla \times \mathbf B_0) \times \mathbf B_0\); \(\gamma = 5/3\). These are the ideal MHD equations for \(\rho = \rho_0 + \tilde \rho\), \(\mathbf U = \tilde{\mathbf U}\), \(q = q_0 + \tilde q\), \(\mathbf B = \mathbf B_0 + \tilde{\mathbf B}\), with the equilibrium force balance subtracted. Dropping the terms quadratic in the perturbations gives LinearMHDWithQ.

Conserved quantities (for \(\tilde{\mathbf U} \cdot \mathbf n = 0\) and \(\mathbf B \cdot \mathbf n = 0\) on the boundary, or periodic boundaries):

Mass perturbation and divergence of the magnetic field:

\[\frac{\textnormal d}{\textnormal d t} \int \tilde \rho \, \textnormal d \mathbf x = 0\,, \qquad \frac{\partial}{\partial t} \nabla \cdot \tilde{\mathbf B} = 0\]

Energy of the perturbation, for an equilibrium with \(\nabla q_0 = 0\) and \(\nabla \times \mathbf B_0 = 0\):

\[\frac{\textnormal d}{\textnormal d t} \int \left( \frac{\rho_0 + \tilde \rho}{2} |\tilde{\mathbf U}|^2 + \frac{2 q_0 \tilde q + \tilde q^2}{\gamma - 1} + \mathbf B_0 \cdot \tilde{\mathbf B} + \frac 12 |\tilde{\mathbf B}|^2 \right) \textnormal d \mathbf x = 0\]

classmethod doc_discretization()[source]#

Time integration is performed by the following propagators (in sequence; the magnetic field is in \(H( extnormal{div})\)):

shear_alf: ShearAlfvenPropagator
nonlin_alf: ShearAlfvenNonlinear
bg_current: BackgroundCurrent
bg_q: BackgroundQGradient
sound: SoundWaveLinearWithQ
nonlin_sound: SoundWaveNonlinearWithQ
mom_adv: MomentumAdvection

classmethod doc_normalization()[source]#

All velocities are normalized by the Alfvén velocity:

\[\hat{U} = \hat{v}_\mathrm{A} = \frac{\hat{B}_0}{\sqrt{\mu_0 \rho_0}}\,, \qquad \hat{q} = \sqrt{\hat p} = \frac{\hat B_0}{\sqrt{\mu_0}}\]

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

classmethod doc_scalar_quantities()[source]#

The following scalars are tracked during simulation:

Kinetic energy: \(E_U = \frac{1}{2} \int (\rho_0 + \tilde \rho) |\tilde{\mathbf{U}}|^2 \, \mathrm{d}V\) (en_U), with the total density (projected \(\rho_0\))
Quadratic thermal energy of the linear sound wave step: \(\frac{2}{\gamma}\int \tilde{q}^2 \, \mathrm{d}V\) (en_thermal); of the nonlinear sound wave step: \(\frac{1}{\gamma - 1}\int \tilde{q}^2 \, \mathrm{d}V\) (en_q2)
Internal energy of the perturbation: \(E_\mathrm{int} = \int \frac{2 q_0 \tilde q + \tilde q^2}{\gamma - 1} \, \mathrm{d}V\) (en_internal)
Magnetic energy (perturbation): \(E_B = \frac{1}{2} \int |\tilde{\mathbf{B}}|^2 \, \mathrm{d}V\) (en_B), equilibrium magnetic energy (en_B_eq), total magnetic energy \(\frac 12 \int |\mathbf B_0 + \tilde{\mathbf B}|^2 \, \mathrm{d}V\) (en_B_tot)
Equilibrium internal energy: \(E_{p0} = \int \frac{p_0}{\gamma - 1} \, \mathrm{d}V\) (en_p_eq)
Energy of the perturbation: \(E_{\mathrm{tot}} = E_U + E_\mathrm{int} + \frac 12 \int (|\mathbf B_0 + \tilde{\mathbf B}|^2 - |\mathbf B_0|^2) \, \mathrm{d}V\) (en_tot), conserved by the equations for \(\nabla q_0 = 0\), \(\nabla \times \mathbf B_0 = 0\). Each step of the splitting conserves its own energy exactly, so \(E_{\mathrm{tot}}\) changes by the splitting error only.

classmethod doc_long_description()[source]#

This model simulates finite-amplitude perturbations in a magnetized plasma around a static MHD equilibrium. It solves the ideal MHD equations for the perturbations of the density, the velocity, the square root of the pressure, \(q = \sqrt p\), and the magnetic field, written as the linear system of LinearMHDWithQ plus the nonlinear terms.

The time step reuses the steps of the linear model (with the inertia of the total density \(\rho_0 + \tilde \rho^n\) and without the continuity equation) and adds the steps of the nonlinear terms, each energy conserving for its own energy:

nonlin_alf: the shear-Alfvén step with the weight \(\mathbf B_0 \to \tilde{\mathbf B}\), i.e. the force \((\nabla \times \tilde{\mathbf B}) \times \tilde{\mathbf B}\) and the induction term \(\nabla \times (\tilde{\mathbf U} \times \tilde{\mathbf B})\);
nonlin_sound: the terms of the sound waves quadratic in \(\tilde q\);
mom_adv: the continuity equation and the advection of the momentum (discrete gradient, Picard iteration).

The weights of the nonlinear terms are at \(t^n\) by default (Options(weights="frozen")) or at \(t^{n+1/2}\) (weights="midpoint", Picard iteration, second order in time). The magnetic field is in \(H(\textnormal{div})\) (the shear-Alfvén step of LinearMHD with b_space="Hdiv"). No numerical dissipation is added, so the model is meant for smooth solutions.

classmethod doc_examples()[source]#

Create and initialize a perturbative MHD model in \(q\):

from struphy.models import PerturbativeMHDWithQ

model = PerturbativeMHDWithQ()

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

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

classmethod doc_use_cases()[source]#

This model is appropriate for finite-amplitude perturbations of a static magnetized plasma.

Typical use cases include:

nonlinear Alfvén and magnetosonic waves
nonlinear saturation of instabilities of a static equilibrium (in the absence of shocks)
comparison with LinearMHDWithQ in the small-amplitude limit

classmethod doc_cannot_be_used_for()[source]#

This model is not suitable for:

shocks and other discontinuities (no numerical dissipation)
equilibria with non-zero background flow
dissipative effects such as resistivity or viscosity (reconnection)
the magnetic field in \(H(\textnormal{curl})\)
kinetic or particle-field effects beyond the fluid MHD description

Maxwell Family#

Maxwell’s equations and their extensions by plasma currents. Each model extends its parent by new terms of the equations and inherits its species, propagators and scalars (see Anatomy of a model).

Maxwell
- ColdPlasma
  - ColdPlasmaVlasov

class struphy.models.Maxwell(base_units: BaseUnits = BaseUnits(x=1.0, B=1.0, n=1.0, kBT=None))[source]#

Bases: StruphyModel

Maxwell’s equations in vacuum for electromagnetic field evolution.

The models of the Maxwell family inherit from this model: they call super().__init__() and extend the hooks of _setup() (species, propagators, variable links, scalars).

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

classmethod doc_pde()[source]#

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

Conserved quantities (for \(\mathbf E \times \mathbf n = 0\) on the boundary, or periodic boundaries):

Electromagnetic energy:

\[\frac{\textnormal d}{\textnormal d t} \frac 12 \int \left( |\mathbf E|^2 + |\mathbf B|^2 \right) \textnormal d \mathbf x = 0\]

classmethod doc_normalization()[source]#

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()[source]#

The following scalars are tracked during simulation:

Electric energy: \(E_E = \frac{1}{2} \int |\mathbf E|^2 \, dV\) (en_E)
Magnetic energy: \(E_B = \frac{1}{2} \int |\mathbf B|^2 \, dV\) (en_B)
Total energy: \(E_{tot} = E_E + E_B\) (en_tot)

classmethod doc_long_description()[source]#

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.

This is the parent model of the Maxwell family: ColdPlasma adds the current of a cold electron fluid, and ColdPlasmaVlasov adds a hot kinetic species; both inherit the species, propagators and scalars of this model.

classmethod doc_examples()[source]#

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()[source]#: Propagation of electromagnetic waves in vacuum.

classmethod doc_cannot_be_used_for()[source]#: Plasma dynamics, plasma-field interactions, or any scenario involving charged particles. This model does not include any particle species or coupling to matter.

classmethod doc_discretization() → str#

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

maxwell: MaxwellWeakAmpere

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)[source]#

Bases: Maxwell

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

Extends Maxwell by the current of a cold electron fluid.

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()[source]#

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.

Conserved quantities (for \(\mathbf E \times \mathbf n = 0\) on the boundary, or periodic boundaries):

Energy:

\[\frac{\textnormal d}{\textnormal d t} \frac 12 \int \left( |\mathbf E|^2 + |\mathbf B|^2 + \frac{\alpha^2}{n_0} |\mathbf j|^2 \right) \textnormal d \mathbf x = 0\]

classmethod doc_normalization()[source]#

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()[source]#

The following scalars are tracked during simulation:

Electric field energy: \(\frac 12 \int |\mathbf E|^2 \, dV\) (en_E)
Magnetic field energy: \(\frac 12 \int |\mathbf B|^2 \, dV\) (en_B)
Cold-current energy: \(\frac{\alpha^2}{2} \int |\mathbf j|^2 / n_0 \, dV\) (en_J)
Total energy: en_tot (sum of the above)

classmethod doc_long_description()[source]#

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.

In the Maxwell family it extends Maxwell by the cold electron current; its child ColdPlasmaVlasov adds a hot kinetic species.

classmethod doc_examples()[source]#

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()[source]#

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()[source]#

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

classmethod doc_discretization() → str#

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

maxwell: MaxwellWeakAmpere
ohm: OhmCold
jxb: JxBCold

class struphy.models.ColdPlasmaVlasov(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, hot_charge_number: int = -1, hot_mass_number: float = 0.0005446623093681918, hot_epsilon: float | None = None)[source]#

Bases: ColdPlasma

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

Extends ColdPlasma by a hot kinetic electron species (Vlasov equation) and its current in Ampère’s law.

Parameters:

base_units (BaseUnits) – Base units for normalization (default: BaseUnits())
charge_number (int) – Charge number of the cold (fluid) electron species (default: -1)
mass_number (float) – Mass number of the cold (fluid) electron species (default: 1/1836)
alpha (float, optional) – Dimensionless parameter: cold-species plasma frequency / cyclotron frequency. If None, computed from units and charge/mass numbers.
epsilon (float, optional) – Normalized cyclotron period of the cold species. If None, computed from units and charge/mass numbers.
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)
hot_epsilon (float, optional) – Normalized cyclotron period of the hot species. If None, computed from units and charge/mass numbers.

classmethod doc_pde()[source]#

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_\textnormal{c}^2}{\varepsilon_\textnormal{c}} \mathbf{j}_\textnormal{c} + \frac{\alpha_\textnormal{h}^2}{\varepsilon_\textnormal{h}} \int_{\mathbb{R}^3} \mathbf{v} f \, \text{d}^3 \mathbf{v}\]

where \((n_0, \mathbf{B}_0)\) denotes an inhomogeneous background, and \(\alpha_\textnormal{c}, \varepsilon_\textnormal{c}\) and \(\alpha_\textnormal{h}, \varepsilon_\textnormal{h}\) are the parameters of the cold and the hot species.

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

\[\nabla \cdot \mathbf{E} = \frac{\alpha_\textnormal{h}^2}{\varepsilon_\textnormal{h}} \int_{\mathbb{R}^3} (f - f_0) \, \text{d}^3 \mathbf{v}\]

Conserved quantities (for \(\mathbf E \times \mathbf n = 0\) on the boundary, or periodic boundaries):

Energy:

\[\frac{\textnormal d}{\textnormal d t} \left[ \frac 12 \int \left( |\mathbf E|^2 + |\mathbf B|^2 + \frac{\alpha_\textnormal{c}^2}{n_0} |\mathbf j_\textnormal{c}|^2 \right) \textnormal d \mathbf x + \frac{\alpha_\textnormal{h}^2}{2} \int \int_{\mathbb{R}^3} |\mathbf v|^2 f \, \text{d}^3 \mathbf{v} \, \textnormal d \mathbf x \right] = 0\]

classmethod doc_normalization()[source]#

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()[source]#

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 (sum of the above)

classmethod doc_long_description()[source]#

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.

In the Maxwell family it extends ColdPlasma by the hot species and inherits its species (em_fields, electrons), propagators and scalars.

classmethod doc_examples()[source]#

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.electrons.current
model.hot_electrons.var

classmethod doc_use_cases()[source]#

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()[source]#

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

classmethod doc_discretization() → str#

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

maxwell: MaxwellWeakAmpere
ohm: OhmCold
jxb: JxBCold
push_eta: PushEta
push_vxb: PushVxB
coupling_va: VlasovAmpereCoupling

SPH Family#

Fluid models discretized with smoothed particle hydrodynamics (SPH). Each model extends its parent by new terms of the equations and inherits its species, propagators and scalars (see Anatomy of a model).

ColdEulerSPH
- EulerSPH
  - ViscousEulerSPH
- IncompressibleNavierStokesSPH

class struphy.models.ColdEulerSPH(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)[source]#

Bases: StruphyModel

Pressureless (cold) Euler equations with external forces, discretized with smoothed particle hydrodynamics (SPH).

The models of the SPH family inherit from this model: they call super().__init__() and extend the hooks of _setup() (species, propagators, variable links, scalars). Formerly named PressureLessSPH.

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 of the forces (default is None, computed from the units).

classmethod doc_pde()[source]#

PDEs solved by model:

Continuity:

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

Momentum:

\[\partial_t \mathbf{u} + (\mathbf{u} \cdot \nabla) \mathbf{u} = \frac{1}{\varepsilon} \left( \mathbf{F} + \mathbf{u} \times \mathbf{B}_0 \right)\]

where \(\mathbf{F}\) is an external force field per unit charge, normalized like an electric field (e.g. \(\mathbf F = -\nabla \phi\) with a given potential \(\phi\)), \(\mathbf{B}_0\) the magnetic field of the equilibrium (if any) and \(\varepsilon\) the species parameter epsilon.

Conserved quantities (periodic boundaries or no outflow):

Mass:

\[\frac{\textnormal d}{\textnormal d t} \int \rho \, \textnormal d \mathbf x = 0\]

Kinetic energy, for \(\mathbf F = 0\) (the Lorentz force does no work):

\[\frac{\textnormal d}{\textnormal d t} \frac 12 \int \rho |\mathbf u|^2 \, \textnormal d \mathbf x = 0\]

classmethod doc_normalization()[source]#

There is no characteristic velocity: the velocity unit \(\hat u\) follows from the base units (\(\hat t = \hat x / \hat u\)). The force field and the magnetic field are normalized like an electric and a magnetic field,

\[\hat{F} = \hat{u} \hat{B}\,, \qquad \varepsilon = \frac{1}{\hat \Omega_\textnormal{c} \hat t}\,,\]

where \(\hat \Omega_\textnormal{c} = q \hat B / m\) is the cyclotron frequency of the species.

classmethod doc_scalar_quantities()[source]#

The following scalars are tracked during simulation:

SPH kinetic energy: \(\frac 12 \int \rho |\mathbf u|^2 \, dV\) (en_kin)

classmethod doc_long_description()[source]#

ColdEulerSPH (formerly 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.

This is the parent model of the SPH family: EulerSPH adds the pressure force (and its child ViscousEulerSPH the viscous force), and IncompressibleNavierStokesSPH the incompressibility constraint and the viscous force; all inherit the species, propagators and scalars of this model.

classmethod doc_examples()[source]#

Create and initialize a cold Euler SPH model:

from struphy.models import ColdEulerSPH

model = ColdEulerSPH()
model.fluid.var

classmethod doc_use_cases()[source]#

This model is appropriate for:

meshfree pressureless-flow benchmarks
particle transport with external forcing and in a background magnetic field
testing SPH particle loading and diagnostics

classmethod doc_cannot_be_used_for()[source]#

This model is not suitable for:

fluids with pressure forces (see EulerSPH)
FEEC-based grid discretizations
self-consistent electromagnetic plasma dynamics
viscous or thermal closures

classmethod doc_discretization() → str#

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

push_eta: PushEta
push_vxb: PushVxB
push_v: PushVinForceField

class struphy.models.EulerSPH(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)[source]#

Bases: ColdEulerSPH

Euler equations (compressible, inviscid) discretized with smoothed particle hydrodynamics (SPH).

Extends ColdEulerSPH by the pressure force.

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)

classmethod doc_pde()[source]#

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) + \rho \mathbf g + \frac{\rho}{\varepsilon} \left( \mathbf{F} + \mathbf{u} \times \mathbf{B}_0 \right)\]

Entropy:

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

where \(S\) denotes the entropy per unit mass, \(\mathbf g\) a constant gravity, \(\mathbf{F}\) an external force field per unit charge, \(\mathbf{B}_0\) the magnetic field of the equilibrium (if any) and \(\varepsilon\) the species parameter epsilon (see ColdEulerSPH). 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}\]

Conserved quantities (periodic boundaries or no outflow):

Mass:

\[\frac{\textnormal d}{\textnormal d t} \int \rho \, \textnormal d \mathbf x = 0\]

Energy, for \(\mathbf F = 0\) and \(\mathbf g = 0\) (the Lorentz force does no work):

\[\frac{\textnormal d}{\textnormal d t} \int \rho \left( \frac 12 |\mathbf u|^2 + \mathcal U(\rho, S) \right) \textnormal d \mathbf x = 0\]

classmethod doc_normalization()[source]#: The characteristic speed is thermal:

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

classmethod doc_scalar_quantities()[source]#

The following scalars are tracked during simulation:

SPH kinetic energy: \(\frac 12 \int \rho |\mathbf u|^2 \, dV\) (en_kin)

classmethod doc_long_description()[source]#

EulerSPH is the particle-based SPH model of a compressible, inviscid fluid. The pressure force follows from the internal energy of the SPH density (isothermal or polytropic closure, selected with the options of push_sph_p). It is intended for meshfree fluid experiments rather than FEEC-based field simulations.

In the SPH family it extends ColdEulerSPH by the pressure force; its child ViscousEulerSPH adds the viscous force.

classmethod doc_examples()[source]#

Create and initialize the Euler SPH model:

from struphy.models import EulerSPH

model = EulerSPH()
model.fluid.var

classmethod doc_use_cases()[source]#

This model is appropriate for:

SPH verification with pressure forces: sound waves, gas expansion
meshfree compressible-fluid experiments
testing the SPH pressure pusher

classmethod doc_cannot_be_used_for()[source]#

This model is not suitable for:

viscous flows (see ViscousEulerSPH)
FEEC grid-based fluid or MHD formulations
entropy-resolved thermodynamic evolution
kinetic plasma physics

classmethod doc_discretization() → str#

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

push_eta: PushEta
push_vxb: PushVxB
push_v: PushVinForceField
push_sph_p: PushVinSPHpressure

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)[source]#

Bases: EulerSPH

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

Extends EulerSPH by the viscous force. Without pressure force (viscous diffusion of a pressureless fluid), set kappa=0.0 in the options of push_sph_p.

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)

classmethod doc_pde()[source]#

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} + \rho \mathbf g + \frac{\rho}{\varepsilon} \left( \mathbf{F} + \mathbf{u} \times \mathbf{B}_0 \right)\]

Entropy:

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

where \(S\) denotes the entropy per unit mass, \(\boldsymbol{\pi}\) is the viscous stress tensor and \(\mathbf g\), \(\mathbf F\), \(\mathbf B_0\), \(\varepsilon\) are as in EulerSPH.

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

\[\boldsymbol{\pi} = -\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}\]

Conserved quantities (periodic boundaries or no outflow):

Mass:

\[\frac{\textnormal d}{\textnormal d t} \int \rho \, \textnormal d \mathbf x = 0\]

The energy of EulerSPH is dissipated by the viscosity.

classmethod doc_normalization()[source]#: The characteristic speed is thermal:

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

classmethod doc_scalar_quantities()[source]#

The following scalars are tracked during simulation:

SPH kinetic energy: \(\frac 12 \int \rho |\mathbf u|^2 \, dV\) (en_kin)

classmethod doc_long_description()[source]#

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

In the SPH family it extends EulerSPH by the viscous force and inherits its species, propagators and scalars.

classmethod doc_examples()[source]#

Create and initialize the viscous Euler SPH model:

from struphy.models import ViscousEulerSPH

model = ViscousEulerSPH()
model.fluid.var

classmethod doc_use_cases()[source]#

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()[source]#

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

classmethod doc_discretization() → str#

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

push_eta: PushEta
push_vxb: PushVxB
push_v: PushVinForceField
push_sph_p: PushVinSPHpressure
push_viscous: PushVinViscousPotential

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)[source]#

Bases: ColdEulerSPH

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

Extends ColdEulerSPH by the viscous force and the incompressibility constraint, whose Lagrange multiplier is the pressure. Without viscosity (incompressible Euler equations), set mu=0.0 in the options of push_viscous.

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)

classmethod doc_pde()[source]#

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} + \frac{\rho}{\varepsilon} \left( \mathbf{F} + \mathbf{u} \times \mathbf{B}_0 \right) \\[2mm] \nabla \cdot \mathbf{u} = 0\,,\end{split}\]

where the pressure \(p\) is the Lagrange multiplier of the incompressibility constraint, \(\boldsymbol{\sigma}\) is the viscous stress tensor, and \(\mathbf F\), \(\mathbf B_0\), \(\varepsilon\) are as in ColdEulerSPH (here \(\varepsilon = 1\)).

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.

Conserved quantities (periodic boundaries or no outflow):

Mass (volume):

\[\frac{\textnormal d}{\textnormal d t} \int \rho \, \textnormal d \mathbf x = 0\]

Kinetic energy for \(\mu = 0\), \(\mathbf F = 0\) (the Lorentz force does no work); it is dissipated by the viscosity:

\[\frac{\textnormal d}{\textnormal d t} \frac 12 \int \rho |\mathbf u|^2 \, \textnormal d \mathbf x = \int \boldsymbol \sigma : \nabla \mathbf u \, \textnormal d \mathbf x \leq 0\]

classmethod doc_normalization()[source]#: The characteristic speed is thermal:

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

classmethod doc_scalar_quantities()[source]#

The following scalars are tracked during simulation:

SPH kinetic energy: \(\frac 12 \int \rho |\mathbf u|^2 \, dV\) (en_kin)

classmethod doc_discretization() → str[source]#

Time integration is performed by the following propagators (in sequence). The pressure is computed with Chorin’s projection method: pressure_poisson solves a Poisson equation for the pressure with the divergence of the velocity (divided by the time step) on the right-hand side, and chorin_projection pushes the velocities with the pressure force, which makes them divergence free. This needs Time(split_algo="LieTrotter").

push_eta: PushEta
push_vxb: PushVxB
push_v: PushVinForceField
push_viscous: PushVinViscousPotential
pressure_poisson: PoissonSolve
chorin_projection: PushVinForceField

classmethod doc_long_description()[source]#

IncompressibleNavierStokesSPH is the particle-based SPH model of an incompressible, viscous fluid. The pressure is not given by an equation of state but is the Lagrange multiplier of the incompressibility constraint, computed on the FEEC grid with Chorin’s projection method.

In the SPH family it extends ColdEulerSPH by the viscous force and the incompressibility constraint, and inherits its species, propagators and scalars.

classmethod doc_examples()[source]#

Create and initialize the incompressible Navier-Stokes SPH model:

from struphy.models import IncompressibleNavierStokesSPH

model = IncompressibleNavierStokesSPH()
model.fluid.var
model.lagrange_multiplier.pressure

classmethod doc_use_cases()[source]#

This model is appropriate for:

incompressible flows, e.g. the relaxation of a shear flow in a channel with no-slip walls
SPH verification with a projection method for the pressure
testing the coupling of SPH particles to FEEC Poisson solves

classmethod doc_cannot_be_used_for()[source]#

This model is not suitable for:

compressible flows (see EulerSPH)
FEEC grid-based fluid or MHD formulations
Strang splitting (the projection needs LieTrotter)
kinetic plasma physics

Other models#

Models that are not yet arranged in families of inherited models.

Gyrokinetic#

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)[source]#

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()[source]#

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()[source]#: 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()[source]#

The following scalars are tracked during simulation:

Electrostatic field energy: en_phi
Particle kinetic energy: en_particles
Total energy: en_tot

classmethod doc_discretization()[source]#

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

classmethod doc_long_description()[source]#: 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()[source]#

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()[source]#

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()[source]#

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.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)[source]#

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()[source]#

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()[source]#

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()[source]#

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()[source]#

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

classmethod doc_long_description()[source]#: 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()[source]#

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()[source]#

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()[source]#

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, b_space: Literal['Hdiv', 'Hcurl'] = 'Hdiv', projected: bool = False)[source]#

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)
b_space ({"Hdiv", "Hcurl"}) – Space of the magnetic field (default: “Hdiv”). "Hdiv": strong induction equation, implicit solve for the velocity (ShearAlfvenPropagator). "Hcurl": weak induction equation, implicit curl-curl solve for the magnetic field (ShearAlfvenCurlCurlPropagator).
projected (bool) – Only with b_space="Hcurl": the curl-curl term of the implicit step contains the projection of the Lorentz force, which makes the step exactly energy conserving at the price of inner mass matrix solves (default: False).

classmethod doc_pde()[source]#

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

Here, \(\tilde{\mathbf U}\) is the perturbed velocity, \(\tilde{\mathbf B}\) the perturbed magnetic field, and \(\rho_0\), \(\mathbf B_0\) are the mass density and magnetic field of a static MHD equilibrium.

Conserved quantities:

Perturbed energy (for \(\tilde{\mathbf U} \cdot \mathbf n = 0\) and \(\mathbf B_0 \cdot \mathbf n = 0\) on the boundary, or periodic boundaries):

\[\frac{\textnormal d}{\textnormal d t} \int \left( \frac{\rho_0}{2} |\tilde{\mathbf U}|^2 + \frac 12 |\tilde{\mathbf B}|^2 \right) \textnormal d \mathbf x = 0\]

Divergence of the magnetic field:

\[\frac{\partial}{\partial t} \nabla \cdot \tilde{\mathbf B} = 0\]

classmethod doc_normalization()[source]#: All velocities are normalized by the Alfvén speed,

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

classmethod doc_scalar_quantities()[source]#

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) with b_space="Hdiv" and with b_space="Hcurl", projected=True; with b_space="Hcurl", projected=False it is not increased (numerical damping that vanishes with the resolution along the background field).
Total energy including background contribution: en_tot2

classmethod doc_discretization()[source]#

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

ShearAlfvenPropagator (b_space="Hdiv") or ShearAlfvenCurlCurlPropagator (b_space="Hcurl")

classmethod doc_long_description()[source]#: 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. Two discretizations are available through the flag b_space. With "Hdiv" (default), the magnetic field is a 2-form, the induction equation holds strongly and the implicit step solves for the velocity. With "Hcurl", the magnetic field is a 1-form, the induction equation is weak and the implicit step solves a field-aligned curl-curl system for the magnetic field, preconditioned by the Hiptmair-Xu auxiliary space preconditioner, which allows large time steps. This step damps the waves numerically, with a damping that vanishes quickly with the resolution along the background field; projected=True gives the exactly energy conserving version, at the price of an inner mass matrix solve in every application of the Schur operator (and no effective curl-curl preconditioner). With b_space="Hcurl", essential boundary conditions in the de Rham sequence set the tangential components of the magnetic field to zero on that boundary (stronger than \(\mathbf B \cdot \mathbf n = 0\)), so free boundary conditions are recommended for bounded domains.

classmethod doc_examples()[source]#

Create and initialize a shear-Alfvén model, with the magnetic field in \(H(\textnormal{div})\) or in \(H(\textnormal{curl})\) (exactly energy conserving step):

from struphy.models import ShearAlfven

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

model_cc = ShearAlfven(b_space="Hcurl", projected=True)
model_cc.propagators.shear_alf.options = model_cc.propagators.shear_alf.Options(
    precond="MassMatrixPreconditioner"
)

classmethod doc_use_cases()[source]#

This model is appropriate for:

isolated shear-Alfvén wave propagation
large time steps with the curl-curl discretization (b_space="Hcurl")
verification of the shear-Alfvén propagators
reduced linear-MHD studies without compressibility

classmethod doc_cannot_be_used_for()[source]#

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 (b_space="Hcurl")

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)[source]#

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()[source]#

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()[source]#: The model uses Alfvén-speed normalization together with entropy-based thermodynamic units inherited from the chosen internal-energy functional.

classmethod doc_scalar_quantities()[source]#

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()[source]#

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

VariationalDensityEvolve
VariationalMomentumAdvection
VariationalEntropyEvolve
VariationalMagFieldEvolve
VariationalViscosity (if with_viscosity is True)
VariationalResistivity (if with_resistivity is True)

classmethod doc_long_description()[source]#: 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()[source]#

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()[source]#

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()[source]#

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)[source]#

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()[source]#

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()[source]#: The characteristic flow speed is the Alfvén speed,

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

classmethod doc_scalar_quantities()[source]#

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()[source]#

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

VariationalDensityEvolve
VariationalMomentumAdvection
VariationalPBEvolve
VariationalViscosity (if with_viscosity is True)
VariationalResistivity (if with_resistivity is True)

classmethod doc_long_description()[source]#: 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()[source]#

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()[source]#

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()[source]#

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)[source]#

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()[source]#

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()[source]#: The flow variables are normalized with the Alfvén speed:

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

classmethod doc_scalar_quantities()[source]#

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()[source]#

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

VariationalDensityEvolve
VariationalMomentumAdvection
VariationalQBEvolve
VariationalViscosity (if with_viscosity is True)
VariationalResistivity (if with_resistivity is True)

classmethod doc_long_description()[source]#: 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()[source]#

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()[source]#

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()[source]#

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, crank_nicolson: bool = False, u_space: str = 'H1vec')[source]#

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)
crank_nicolson (bool) – Only for linear=True: whether to advance the velocity, the magnetic field and the pressure with one linear solve per substep (variat_lin_alfven and variat_lin_sound, energy conserving, by default with midpoint weights of the second-order terms, option weights) instead of the Picard iteration of variat_pb (default: False)
u_space (str) – Space of the velocity, "H1vec" (vector field) or "Hdiv" (2-form, as in LinearMHD; for the nonlinear model, and for the linear model with one linear solve per substep, linear=True, crank_nicolson=True) (default: "H1vec")

classmethod doc_pde()[source]#

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, which is quadratic in the perturbations, is not fed back into the pressure perturbation \(\tilde p\): it goes into the second-order auxiliary pressure \(\hat p\) (diagnostic pt3), whose energy \(\int \hat p/(\gamma - 1)\) is part of the conserved total energy en_tot. So with dissipation the total energy is still conserved up to the solver tolerances.

Second-order magnetic energy with one linear solve per substep (crank_nicolson=True): the term \(\int \mathbf B_0 \cdot \hat{\mathbf B}\) is written as \(\int \mathbf J_0 \cdot \hat{\mathbf A}\), with the auxiliary vector potential \(\hat{\mathbf A}\) (diagnostics.at1), \(\partial_t \hat{\mathbf A} = - \tilde{\mathbf B} \times \mathbf u\), \(\hat{\mathbf B} = \nabla \times \hat{\mathbf A}\) (diagnostics.bt2) and the equilibrium current \(\mathbf J_0 = \nabla \times \mathbf B_0\). Both are the same in the continuum (\(\hat{\mathbf A} \times \mathbf n = 0\) at the wall); discretely the force \(\mathbf J_0 \times \tilde{\mathbf B}\) then uses the current of the dual complex, the curl of the 1-form \(\star \mathbf B_0\) without wall conditions, instead of the weak curl of \(\mathbf B_0\), which vanishes tangentially at a wall and causes a spurious wall instability (see VariationalLinearAlfvenEvolve).

Velocity in \(H(\textnormal{div})\) (linear=True, crank_nicolson=True, u_space="Hdiv"): the linear model can also be discretized with the velocity as a 2-form, as in LinearMHD; this is possible because the linear model has no advection of the momentum (no commutator), the only term that needs \(\mathbf u \in (H^1)^3\). The variational structure is the same, so the energy is conserved up to the solver tolerances; the boundary condition is \(\mathbf u \cdot \mathbf n = 0\) instead of \(\mathbf u = 0\). It is available with one linear solve per substep only (crank_nicolson=True). The viscous term is then discretized with the divergence and the curl of \(\mathbf u\), \(\int \mu (\nabla \cdot \mathbf u\, \nabla \cdot \mathbf v + \nabla \times \mathbf u \cdot \nabla \times \mathbf v)\), which equals \(\int \mu \nabla \mathbf u : \nabla \mathbf v\) up to boundary terms. The curl is the weak curl, as the current of the resistivity, so the wall condition is perfect slip, \((\nabla \times \mathbf u) \times \mathbf n = 0\), and the tangential velocity at the wall stays free.

Velocity in \(H(\textnormal{div})\) for the nonlinear model (u_space="Hdiv"): the advection of the momentum (the commutator of the Euler-Poincaré equations) is discretized with the vector fields acting as derivations on an auxiliary space that contains the contravariant velocity components, so that the bracket needs no projection (HdivBracketOperator). It needs the degree \(p \geq 2\), except in periodic directions with one element, and no polar splines. The velocity mass matrix is the 2-form mass matrix of the total density \(\rho_0 + \tilde\rho\), and the energy is conserved up to the solver tolerances, as with the velocity in \((H^1)^3\).

classmethod doc_normalization()[source]#: The velocity scale is Alfvénic:

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

classmethod doc_scalar_quantities()[source]#

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()[source]#

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

VariationalDensityEvolve
VariationalMomentumAdvection (if linear is False)
VariationalPBEvolve, or, if crank_nicolson is True (linear model only), VariationalLinearAlfvenEvolve and VariationalLinearSoundEvolve (one linear solve each, energy conserving, with midpoint weights of the second-order terms by default)
VariationalViscosity (if with_viscosity is True)
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. With u_space="Hdiv" the propagators of the linear model use the velocity mass matrix M2n and, for the equilibrium terms, the basis projection operators of ShearAlfvenPropagator and Magnetosonic; those of the nonlinear model use the 2-form mass matrix of the total density (VelocityMass) and the transport operators of variational_utilities with u_space="Hdiv".

classmethod doc_long_description()[source]#: 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, with crank_nicolson=True in addition without Picard iteration (one linear solve per substep, still energy conserving).

classmethod doc_examples()[source]#

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)
linear_cn_model = VariationalPerturbativeMHD(linear=True, crank_nicolson=True)
linear_model_hdiv = VariationalPerturbativeMHD(linear=True, crank_nicolson=True, u_space="Hdiv")
model_hdiv = VariationalPerturbativeMHD(u_space="Hdiv")
model.mhd.pressure
model.em_fields.b_field

classmethod doc_use_cases()[source]#

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
comparing the variational linear MHD with the velocity in \(H(\textnormal{div})\) (u_space="Hdiv") with LinearMHD
perturbative nonlinear MHD in toroidal geometries with the velocity in \(H(\textnormal{div})\) (u_space="Hdiv"), whose wall condition is only \(\mathbf u \cdot \mathbf n = 0\)

classmethod doc_cannot_be_used_for()[source]#

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, split_sound: bool = True, crank_nicolson: bool = False, u_space: str = 'H1vec')[source]#

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)
split_sound (bool) – Whether to advance the coupling of the velocity to the magnetic field (shear-Alfvén/fast waves, variat_alfven) and to q (sound waves, variat_sound) in two substeps, both energy conserving (default: True). With False they are advanced together in variat_qb, without the splitting error but with a larger and stiffer Newton system.
crank_nicolson (bool) – Only for linear=True: whether to advance the velocity, the magnetic field and q with one linear solve per substep (variat_lin_alfven and variat_lin_sound, energy conserving; midpoint weights of the second-order term by default in variat_lin_alfven, option weights, weights at t^n in variat_lin_sound) instead of the Picard iteration of variat_qb (default: False)
u_space (str) – Space of the velocity, "H1vec" (vector field) or, only for the ideal linear model with one linear solve per substep (linear=True, crank_nicolson=True, with_viscosity=False, with_resistivity=False), "Hdiv" (2-form, as in LinearMHD) (default: "H1vec")

classmethod doc_pde()[source]#

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

Second-order magnetic energy with one linear solve per substep (crank_nicolson=True): the term \(\int \mathbf B_0 \cdot \hat{\mathbf B}\) is written as \(\int \mathbf J_0 \cdot \hat{\mathbf A}\), with the auxiliary vector potential \(\hat{\mathbf A}\) (diagnostics.at1), \(\partial_t \hat{\mathbf A} = - \tilde{\mathbf B} \times \mathbf u\), \(\hat{\mathbf B} = \nabla \times \hat{\mathbf A}\) (diagnostics.bt2) and the equilibrium current \(\mathbf J_0 = \nabla \times \mathbf B_0\). Both are the same in the continuum (\(\hat{\mathbf A} \times \mathbf n = 0\) at the wall); discretely the force \(\mathbf J_0 \times \tilde{\mathbf B}\) then uses the current of the dual complex, the curl of the 1-form \(\star \mathbf B_0\) without wall conditions, instead of the weak curl of \(\mathbf B_0\), which vanishes tangentially at a wall and causes a spurious wall instability (see VariationalLinearAlfvenEvolve).

Velocity in \(H(\textnormal{div})\) (linear=True, crank_nicolson=True, u_space="Hdiv"): the linear model can also be discretized with the velocity as a 2-form, as in LinearMHD; this is possible because the linear model has no advection of the momentum (no commutator), the only term that needs \(\mathbf u \in (H^1)^3\). The variational structure is the same, so the energy is conserved up to the solver tolerances; the boundary condition is \(\mathbf u \cdot \mathbf n = 0\) instead of \(\mathbf u = 0\). It is available for the ideal model with one linear solve per substep only (crank_nicolson=True, with_viscosity=False, with_resistivity=False).

classmethod doc_normalization()[source]#: The model uses Alfvén-speed normalization:

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

classmethod doc_scalar_quantities()[source]#

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}\) (\(\int \mathbf J_0 \cdot \hat{\mathbf A}\) with crank_nicolson=True)
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()[source]#

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

VariationalDensityEvolve
VariationalMomentumAdvection (if linear is False)
VariationalAlfvenEvolve (shear-Alfvén/fast waves) and 4. VariationalSoundEvolve (sound waves), if split_sound is True (default); otherwise VariationalQBEvolve (both together); or, if crank_nicolson is True (linear model only), VariationalLinearAlfvenEvolve and VariationalLinearSoundEvolve (one linear solve each, energy conserving; midpoint weights of the second-order term by default in the Alfvén step, weights at \(t^n\) in the sound step)

VariationalViscosity (if with_viscosity is True)
VariationalResistivity (if with_resistivity is True)

The option model of the variational propagators is "perturbative_q", or "linear_q" if linear is True ("perturbative" or "linear" for VariationalAlfvenEvolve). With u_space="Hdiv" the propagators use the velocity mass matrix M2n and, for the equilibrium terms of the magnetic field and of the density, the basis projection operators of ShearAlfvenPropagator and Magnetosonic.

classmethod doc_long_description()[source]#: 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; with crank_nicolson=True in addition, without Picard iteration (one linear solve per substep, still energy conserving). 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()[source]#

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)
linear_cn_model = VariationalPerturbativeMHD_with_q(
    linear=True, crank_nicolson=True, with_viscosity=False, with_resistivity=False
)
linear_model_hdiv = VariationalPerturbativeMHD_with_q(
    linear=True, crank_nicolson=True, with_viscosity=False, with_resistivity=False, u_space="Hdiv"
)
model.mhd.sqrt_p
model.em_fields.b_field

classmethod doc_use_cases()[source]#

This model is appropriate for:

nonlinear MHD dynamics around a given (possibly non-uniform) zero-flow equilibrium
linear MHD waves and instabilities (linear=True), also with the velocity in \(H(\textnormal{div})\) (u_space="Hdiv")
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()[source]#

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.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)[source]#

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()[source]#

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()[source]#: 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()[source]#

The following scalars are tracked during simulation:

No default scalar diagnostics are defined by this model.

classmethod doc_discretization()[source]#

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

TwoFluidQuasiNeutralFull

classmethod doc_long_description()[source]#: 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()[source]#

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()[source]#

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()[source]#

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.HasegawaWakatani(base_units: BaseUnits = BaseUnits(x=1.0, B=1.0, n=1.0, kBT=None), mass_number: float = 1.0)[source]#

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()[source]#

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()[source]#: 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()[source]#

The following scalars are tracked during simulation:

No default scalar diagnostics are defined by this model.

classmethod doc_discretization()[source]#

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

classmethod doc_long_description()[source]#: 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()[source]#

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()[source]#

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()[source]#

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)[source]#

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()[source]#

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()[source]#: The model uses Alfvén-speed scaling for the velocity and entropy-based thermodynamic units for the internal-energy closure.

classmethod doc_scalar_quantities()[source]#

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()[source]#

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

VariationalDensityEvolve
VariationalMomentumAdvection
VariationalEntropyEvolve
VariationalViscosity (if with_viscosity is True)

classmethod doc_long_description()[source]#: 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()[source]#

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()[source]#

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()[source]#

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)[source]#

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()[source]#

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()[source]#: 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()[source]#

The following scalars are tracked during simulation:

Kinetic energy: kinetic_energy
Thermodynamic energy: thermo_energy
Total energy: total_energy

classmethod doc_discretization()[source]#

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

classmethod doc_long_description()[source]#: 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()[source]#

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()[source]#

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()[source]#

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)[source]#

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()[source]#

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()[source]#: 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()[source]#

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()[source]#

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

classmethod doc_long_description()[source]#: 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()[source]#

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()[source]#

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()[source]#

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)[source]#

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()[source]#

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()[source]#: The flow speed is normalized with the Alfvén speed:

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

classmethod doc_scalar_quantities()[source]#

The following scalars are tracked during simulation:

Kinetic energy: kinetic_energy

classmethod doc_discretization()[source]#

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

classmethod doc_long_description()[source]#: 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()[source]#

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()[source]#

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()[source]#

This model is not suitable for:

pressure- or entropy-driven flow
magnetic-field coupling
viscous/resistive dissipation
kinetic particle physics

Particle diffusion#

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

Bases: StruphyModel

Diffusion equation discretized with a deterministic particle method.

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

classmethod doc_pde()[source]#

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()[source]#

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()[source]#

The following scalars are tracked during simulation:

No default scalar diagnostics are defined by this model.

classmethod doc_discretization()[source]#

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

PushDeterministicDiffusion

classmethod doc_long_description()[source]#: 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()[source]#

Create and initialize a deterministic diffusion model:

from struphy.models import DeterministicParticleDiffusion

model = DeterministicParticleDiffusion()
model.hydrogen.var

classmethod doc_use_cases()[source]#

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()[source]#

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))[source]#

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()[source]#

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()[source]#: The natural scaling is set by the diffusion coefficient:

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

classmethod doc_scalar_quantities()[source]#

The following scalars are tracked during simulation:

No default scalar diagnostics are defined by this model.

classmethod doc_discretization()[source]#

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

PushRandomDiffusion

classmethod doc_long_description()[source]#: 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()[source]#

Create and initialize a random diffusion model:

from struphy.models import RandomParticleDiffusion

model = RandomParticleDiffusion()
model.hydrogen.var

classmethod doc_use_cases()[source]#

This model is appropriate for:

stochastic particle diffusion benchmarks
Monte-Carlo transport verification
comparison against deterministic diffusion solvers

classmethod doc_cannot_be_used_for()[source]#

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)[source]#

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()[source]#

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()[source]#

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()[source]#

The following scalars are tracked during simulation:

No default scalar diagnostics are defined by this model.

classmethod doc_discretization()[source]#

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

TimeDependentSource (if with_t_dep_source is True)
PoissonSolve

classmethod doc_long_description()[source]#: 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()[source]#

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()[source]#

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()[source]#

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)[source]#

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()[source]#

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()[source]#

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()[source]#

The following scalars are tracked during simulation:

No default scalar diagnostics are defined by this model.

classmethod doc_discretization()[source]#

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

TimeDependentSource (if with_t_dep_source is True)
CurlCurlSolve

classmethod doc_long_description()[source]#: 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()[source]#

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()[source]#

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()[source]#

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)[source]#

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()[source]#

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()[source]#

The coefficient scaling is

\[\hat \rho = \hat n.\]

No dedicated velocity normalization is used.

classmethod doc_scalar_quantities()[source]#

The following scalars are tracked during simulation:

No default scalar diagnostics are defined by this model.

classmethod doc_discretization()[source]#

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

TimeDependentSource (if with_t_dep_source is True)
GyrokineticPoissonSolve

classmethod doc_long_description()[source]#: 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()[source]#

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()[source]#

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()[source]#

This model is not suitable for:

self-consistent kinetic plasma evolution on its own
electromagnetic or magnetic-field dynamics

Table of Contents

Previous topic

Next topic

Vlasov Family#

Guiding-Center Family#

LinearEuler Family#

LinearEulerWithQ Family#

Maxwell Family#

SPH Family#

Other models#

Gyrokinetic#

MHD#

Multi-fluid#

Particle diffusion#

Pure fields#