Models

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

Bases: StruphyModel

Maxwell’s equations in vacuum for electromagnetic field evolution.

Parameters:

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

classmethod doc_pde()

PDEs solved by model:

Ampère’s law (no current):

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

Faraday’s law:

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

classmethod doc_normalization()

Velocity and fields are normalized as:

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

where \(c\) is the speed of light.

classmethod doc_scalar_quantities()

The following scalars are tracked during simulation:

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

• Magnetic energy: \(E_B = \frac{1}{2} \int |\mathbf B|^2 \, dV\)

• Total energy: \(E_{total} = E_E + E_B\)

classmethod doc_discretization()

Propagators:

1. Maxwell

classmethod doc_long_description()

This model simulates the propagation of electromagnetic waves in vacuum using Maxwell’s equations without sources. It uses a finite element exterior calculus (FEEC) formulation with the electric field in H(curl) and the magnetic field in H(div) spaces.

classmethod doc_examples()

Create and initialize a Maxwell model:

from struphy.models import Maxwell

model = Maxwell()
# Fields are accessible via:
model.em_fields.e_field
model.em_fields.b_field

classmethod doc_use_cases()

Propagation of electromagnetic waves in vacuum.

classmethod doc_cannot_be_used_for()

Plasma dynamics, plasma-field interactions, or any scenario involving charged particles. This model does not include any particle species or coupling to matter.

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

Bases: StruphyModel

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

Parameters:

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

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

classmethod doc_pde()

PDEs solved by model:

Continuity (mass conservation):

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

Momentum (Lorentz force):

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

Energy (adiabatic process):

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

Induction (Faraday’s law):

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

classmethod doc_normalization()

All velocities are normalized by the Alfvén velocity:

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

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

classmethod doc_scalar_quantities()

The following scalars are tracked during simulation:

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

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

• Internal energy (perturbation): \(E_p = \int \frac{\tilde{p}}{\gamma - 1} \, \mathrm{d}V\) with \(\gamma = 5/3\)

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

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

Propagators:

1. ShearAlfvenPropagator

2. Magnetosonic

classmethod doc_long_description()

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

The linear MHD system supports three wave families:

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

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

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

All evolved quantities are perturbations around a stationary equilibrium:

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

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

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

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

The corresponding equilibrium quantities are:

• \(\rho_0\) - background density

• \(p_0\) - background pressure

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

classmethod doc_examples()

Create and initialize a linear MHD model:

from struphy.models import LinearMHD

model = LinearMHD()

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

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

classmethod doc_use_cases()

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

Typical use cases include:

• verification of linear MHD wave dispersion and mode structure

• perturbative studies around static equilibria

• benchmark problems for ideal-MHD field-fluid coupling

classmethod doc_cannot_be_used_for()

This model is not suitable for:

• nonlinear MHD dynamics with finite-amplitude perturbations

• equilibria with non-zero background flow

• dissipative effects such as resistivity or viscosity

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

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

Bases: StruphyModel

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

Parameters:

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

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

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

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

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

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

classmethod doc_pde()

PDEs solved by model:

Vlasov equation:

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

Ampère’s law:

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

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

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

classmethod doc_normalization()

Velocity and field normalizations:

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

Dimensionless parameters:

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

where

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

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

classmethod doc_scalar_quantities()

The following scalars are tracked during simulation:

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

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

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

classmethod doc_discretization()

classmethod doc_long_description()

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

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

classmethod doc_examples()

Create and initialize a Vlasov-Ampère model:

from struphy.models import VlasovAmpereOneSpecies

model = VlasovAmpereOneSpecies()

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

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

classmethod doc_use_cases()

This model is appropriate for:

• studying kinetic instabilities in collisionless plasmas

• particle-wave interactions and nonlinear plasma physics

• initial value problems with prescribed equilibrium distributions

• benchmark problems for kinetic-field coupling schemes

• verification of PIC methods in simplified geometries

classmethod doc_cannot_be_used_for()

This model is not suitable for:

• multi-species systems (requires extension)

• collision operator effects

• electromagnetic wave propagation (no magnetic field evolution)

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

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

Bases: StruphyModel

Cold plasma model.

Normalization:

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

Equations:

\[\begin{split}\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\,, \\[2mm] &\frac{\partial \mathbf B}{\partial t} + \nabla\times\mathbf E = 0\,, \\[2mm] -&\frac{\partial \mathbf E}{\partial t} + \nabla\times\mathbf B = \frac{\alpha^2}{\varepsilon} \mathbf j \,,\end{split}\]

where \((n_0,\mathbf B_0)\) denotes a (inhomogeneous) background and

\[\alpha = \frac{\hat \Omega_\textnormal{p}}{\hat \Omega_\textnormal{c}}\,, \qquad \varepsilon = \frac{1}{\hat \Omega_\textnormal{c} \hat t}\,.\]

Propagators (called in sequence):

1. Maxwell

2. OhmCold

3. JxBCold

Model info:

classmethod doc_pde()

PDEs solved by model:

Cold-plasma current:

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

Faraday’s law:

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

Ampère’s law:

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

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

classmethod doc_normalization()

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

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

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

classmethod doc_scalar_quantities()

The following scalars are tracked during simulation:

• Electric field energy: electric_energy

• Magnetic field energy: magnetic_energy

• Cold-current energy: kinetic_energy

• Total energy: total_energy

classmethod doc_discretization()

classmethod doc_long_description()

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

classmethod doc_examples()

Create and initialize a cold-plasma model:

from struphy.models import ColdPlasma

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

classmethod doc_use_cases()

This model is appropriate for:

• cold-plasma wave propagation studies

• electromagnetic benchmarks with a fluid current response

• regimes where thermal pressure can be neglected

• algorithm verification for Maxwell plus current coupling

classmethod doc_cannot_be_used_for()

This model is not suitable for:

• finite-temperature or pressure-driven plasma dynamics

• kinetic resonances or velocity-space instabilities

• multi-species hybrid or fully kinetic problems

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

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

Bases: StruphyModel

Cold plasma hybrid model.

Normalization:

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

Equations:

\[\begin{split}&\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 \,, \\[2mm] \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\,, \\[2mm] &\frac{\partial \mathbf B}{\partial t} + \nabla\times\mathbf E = 0\,, \\[2mm] -&\frac{\partial \mathbf E}{\partial t} + \nabla\times\mathbf B = \frac{\alpha^2}{\varepsilon_\textnormal{h}} \left( \mathbf j_\textnormal{c} + \int_{\mathbb{R}^3} \mathbf{v} f \, \text{d}^3 \mathbf{v} \right) \,,\end{split}\]

where \((n_0,\mathbf B_0)\) denotes a (inhomogeneous) background and

\[\alpha = \frac{\hat \Omega_\textnormal{p,cold}}{\hat \Omega_\textnormal{c,cold}}\,, \qquad \varepsilon_\textnormal{c} = \frac{1}{\hat \Omega_\textnormal{c,cold} \hat t}\,, \qquad \varepsilon_\textnormal{h} = \frac{1}{\hat \Omega_\textnormal{c,hot} \hat t} \,.\]

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

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

Propagators (called in sequence):

1. Maxwell

2. OhmCold

3. JxBCold

4. PushVxB

5. PushEta

6. VlasovAmpereCoupling

classmethod doc_pde()

PDEs solved by model:

Hot Vlasov species:

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

Cold-plasma current:

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

Faraday’s law:

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

Ampère’s law:

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

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

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

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

classmethod doc_normalization()

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

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

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

classmethod doc_scalar_quantities()

The following scalars are tracked during simulation:

• Electric field energy: en_E

• Magnetic field energy: en_B

• Cold-current energy: en_J

• Hot-particle kinetic energy: en_f

• Total energy: en_tot

classmethod doc_discretization()

classmethod doc_long_description()

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

classmethod doc_examples()

Create and initialize the cold-plasma plus Vlasov model:

from struphy.models import ColdPlasmaVlasov

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

classmethod doc_use_cases()

This model is appropriate for:

• hybrid electromagnetic problems with one hot kinetic species

• energetic-particle interaction with a cold background plasma

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

• benchmarks of fluid-kinetic current coupling

classmethod doc_cannot_be_used_for()

This model is not suitable for:

• fully warm-fluid or pressure-anisotropic background dynamics

• all-species kinetic simulations

• collision operators or detailed dissipative closures

• electrostatic-only reductions where magnetic evolution is irrelevant

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

Bases: StruphyModel

Diffusion equation discretized with a deterministic particle method; the solution is \(L^2\)-projected onto \(V^0 \subset H^1\) to compute the flux.

Normalization:

\[\hat D := \frac{\hat x^2}{\hat t } \,.\]

Equations: 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.

Propagators (called in sequence):

1. PushDeterministicDiffusion

Model info:

classmethod doc_pde()

PDEs solved by model:

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

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

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

classmethod doc_normalization()

The diffusion coefficient defines the normalization,

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

No separate plasma velocity scale is used in this model.

classmethod doc_scalar_quantities()

The following scalars are tracked during simulation:

• No default scalar diagnostics are defined by this model.

classmethod doc_discretization()

classmethod doc_long_description()

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

classmethod doc_examples()

Create and initialize a deterministic diffusion model:

from struphy.models import DeterministicParticleDiffusion

model = DeterministicParticleDiffusion()
model.hydrogen.var

classmethod doc_use_cases()

This model is appropriate for:

• deterministic particle discretizations of diffusion

• transport benchmarks with positive scalar densities

• numerical comparison against stochastic diffusion methods

classmethod doc_cannot_be_used_for()

This model is not suitable for:

• electromagnetic plasma dynamics

• kinetic Vlasov problems in phase space

• nonlinear fluid systems with pressure or momentum evolution

• diffusion tensors outside the currently supported simplified forms

class struphy.models.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)

Bases: StruphyModel

Drift-kinetic equation for one ion species in static background magnetic field, coupled to quasi-neutrality equation with adiabatic electrons.

Normalization:

\[\hat v = \hat v_\textrm{i} = \sqrt{\frac{k_B \hat T_\textrm{i}}{m_\textrm{i}}}\,,\qquad \hat E = \hat v_\textrm{i}\hat B\,,\qquad \hat \phi = \hat E \hat x \,.\]

Equations:

\[\begin{split}&\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\,. \\[2mm] - &\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 \,.\end{split}\]

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

and with the normalization parameters

\[\varepsilon := \frac{1}{\hat \Omega_\textrm{c} \hat t}\,,\qquad \hat \Omega_\textrm{c} = \frac{q_\textrm{i} \hat B}{m_\textrm{i}} \,.\]

Notes

• The Control variate method in the Poisson equation is optional; in case it is enabled via the parameter file, the following Poisson equation is solved:

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

\[\int \frac{n_0}{|B_0|^2} \nabla_\perp \psi \cdot \nabla_\perp \phi\,\textrm d \mathbf x + \frac{1}{Z\varepsilon^2} \int \frac{n_0}{T_{0}} \psi \phi \,\textrm d \mathbf x = \frac 1 \varepsilon \int \int \psi \, (f - f_0) B^*_\parallel \,\textrm d \mathbf x\,\textnormal d v_\parallel \textnormal d \mu \qquad \forall \ \psi \in H^1\,.\]

Propagators (called in sequence):

1. ImplicitDiffusion

2. PushGuidingCenterBxEstar

3. PushGuidingCenterParallel

Model info:

classmethod doc_pde()

PDEs solved by model:

Drift-kinetic equation:

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

Poisson equation:

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

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

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

Notes

• The control_var in the Poisson equation is optional; in case it is enabled via the parameter file, the following Poisson equation is solved:

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

\[\int \frac{n_0}{|B_0|^2} \nabla_\perp \psi \cdot \nabla_\perp \phi \, \textrm{d} \mathbf{x} + \frac{1}{Z \varepsilon^2} \int \frac{n_0}{T_0} \psi \phi \, \textrm{d} \mathbf{x} = \frac{1}{\varepsilon} \int \int \psi \, (f - f_0) B^*_\parallel \, \textrm{d} \mathbf{x} \, \textnormal{d} v_\parallel \textnormal{d} \mu \qquad \forall \ \psi \in H^1\]

classmethod doc_normalization()

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

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

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

classmethod doc_scalar_quantities()

The following scalars are tracked during simulation:

• Field energy: en_phi

• Guiding-center particle energy: en_particles

• Total energy: en_tot

classmethod doc_discretization()

classmethod doc_long_description()

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

classmethod doc_examples()

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

from struphy.models import DriftKineticElectrostaticAdiabatic

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

classmethod doc_use_cases()

This model is appropriate for:

• electrostatic drift-kinetic ion turbulence studies

• strongly magnetized plasmas with adiabatic electrons

• guiding-center PIC verification in realistic magnetic geometry

• low-frequency regimes where full gyrophase resolution is unnecessary

classmethod doc_cannot_be_used_for()

This model is not suitable for:

• fully electromagnetic dynamics with evolving magnetic perturbations

• electron kinetic effects beyond the adiabatic closure

• problems that require resolving full cyclotron motion

• multi-species kinetic coupling without extending the model

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

Bases: StruphyModel

Guiding-center equation in static background magnetic field.

Normalization:

\[\hat v = \hat v_\textnormal{A} \,.\]

Equations:

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

Moreover,

\[\epsilon = \frac{1 }{ \hat \Omega_{\textnormal{c}} \hat t}\,,\qquad \textnormal{with} \qquad\hat \Omega_{\textnormal{c}} = \frac{Ze \hat B}{A m_\textnormal{H}}\,.\]

Propagators (called in sequence):

1. PushGuidingCenterBxEstar

2. PushGuidingCenterParallel

classmethod doc_pde()

PDEs solved by model:

Guiding-center Vlasov equation:

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

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

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

classmethod doc_normalization()

The reference speed is Alfvénic and the key small parameter is the normalized cyclotron time:

\[\hat v = \hat v_A,\qquad \varepsilon = 1/(\hat\Omega_c \hat t).\]

classmethod doc_scalar_quantities()

The following scalars are tracked during simulation:

• Parallel kinetic energy: en_fv

• Magnetic-moment energy contribution: en_fB

• Total particle energy: en_tot

• Lost markers: n_lost_particles

classmethod doc_discretization()

classmethod doc_long_description()

GuidingCenter is the reduced kinetic model used when fast gyromotion is averaged out and only guiding-center motion needs to be resolved. It evolves particles in a fixed magnetic background without any self-consistent field solve.

classmethod doc_examples()

Create and initialize a guiding-center model:

from struphy.models import GuidingCenter

model = GuidingCenter()
model.kinetic_ions.var

classmethod doc_use_cases()

This model is appropriate for:

• test-particle guiding-center orbit calculations

• strongly magnetized plasmas with scale separation from gyromotion

• verification of 5D guiding-center pushers and diagnostics

classmethod doc_cannot_be_used_for()

This model is not suitable for:

• self-consistent electric or magnetic field evolution

• full-orbit particle dynamics with resolved gyrophase

• collisional transport or source terms not present in the equation

• fluid closures or MHD force balance

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

Bases: StruphyModel

Hasegawa-Wakatani equations in 2D.

Normalization:

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

Equations:

\[\begin{split}&\frac{\partial n}{\partial t} = C (\phi - n) - [\phi, n] - \kappa\, \partial_y \phi + \nu\, \nabla^{2N} n\,, \\[2mm] &\frac{\partial \omega}{\partial t} = C (\phi - n) - [\phi, \omega] + \nu\, \nabla^{2N} \omega \,, \\[3mm] &\Delta \phi = \omega\,,\end{split}\]

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

Propagators (called in sequence):

1. PoissonFieldSolve

2. HasegawaWakataniStep

Model info:

classmethod doc_pde()

PDEs solved by model:

Density equation:

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

Vorticity equation:

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

Potential equation:

\[\Delta \phi = \omega\]

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

classmethod doc_normalization()

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

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

classmethod doc_scalar_quantities()

The following scalars are tracked during simulation:

• No default scalar diagnostics are defined by this model.

classmethod doc_discretization()

classmethod doc_long_description()

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

classmethod doc_examples()

Create and initialize a Hasegawa-Wakatani model:

from struphy.models import HasegawaWakatani

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

classmethod doc_use_cases()

This model is appropriate for:

• 2D resistive drift-wave turbulence studies

• reduced electrostatic edge-plasma benchmarks

• testing Poisson-coupled advection-diffusion solvers

classmethod doc_cannot_be_used_for()

This model is not suitable for:

• three-dimensional electromagnetic turbulence

• full kinetic Landau or cyclotron physics

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

• self-consistent magnetic-field evolution

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

Bases: StruphyModel

Linear extended MHD with zero-flow equilibrium (\(\mathbf U_0 = 0\)). For uniform background conditions only.

Normalization:

\[\hat U = \hat v_\textnormal{A} \,.\]

Equations:

\[\begin{split}&\frac{\partial \tilde \rho}{\partial t}+\nabla\cdot(\rho_0 \tilde{\mathbf{U}})=0\,, \\[2mm] \rho_0&\frac{\partial \tilde{\mathbf{U}}}{\partial t} + \nabla \tilde p =(\nabla\times \tilde{\mathbf{B}})\times\mathbf{B}_0 \,, \\[2mm] &\frac{\partial \tilde p}{\partial t} + \frac{5}{3}\,p_{0}\nabla\cdot \tilde{\mathbf{U}}=0\,, \\[2mm] &\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\,.\end{split}\]

where

\[\varepsilon = \frac{1}{\hat \Omega_{\textnormal{c}} \hat t}\,,\qquad \textnormal{with} \qquad\hat \Omega_{\textnormal{c}} = \frac{Ze \hat B}{A m_\textnormal{H}}\,.\]

Propagators (called in sequence):

1. ShearAlfvenB1

2. Hall

3. MagnetosonicUniform

Model info:

classmethod doc_pde()

PDEs solved by model:

Continuity:

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

Momentum:

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

Pressure:

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

Induction:

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

classmethod doc_normalization()

All velocities are normalized by the Alfvén speed,

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

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

classmethod doc_scalar_quantities()

The following scalars are tracked during simulation:

• Flow kinetic energy: en_U

• Thermal pressure energy: en_p

• 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: en_tot

• Magnetic helicity-like invariant: helicity

classmethod doc_discretization()

classmethod doc_long_description()

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

classmethod doc_examples()

Create and initialize a linear extended-MHD model:

from struphy.models import LinearExtendedMHDuniform

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

classmethod doc_use_cases()

This model is appropriate for:

• linear Hall-MHD wave studies

• uniform-equilibrium benchmarks for extended MHD

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

classmethod doc_cannot_be_used_for()

This model is not suitable for:

• nonlinear extended-MHD turbulence

• non-uniform equilibrium configurations

• kinetic ion/electron effects beyond the Hall correction

• dissipation-dominated problems with explicit viscosity or resistivity

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

Bases: StruphyModel

Hybrid linear ideal MHD + energetic ions (5D Driftkinetic) with current coupling scheme.

Normalization:

\[\hat U = \hat v =: \hat v_\textnormal{A, bulk} \,, \qquad \hat f_\textnormal{h} = \frac{\hat n}{\hat v_\textnormal{h} \hat \mu \hat B} \,,\qquad \hat \mu = \frac{A_\textnormal{h} m_\textnormal{H} \hat v_\textnormal{h}^2}{\hat B} \,.\]

Equations:

\[\begin{split}\begin{align} \textnormal{MHD} &\left\{ \begin{aligned} &\frac{\partial \tilde{\rho}}{\partial t}+\nabla\cdot(\rho_{0} \tilde{\mathbf{U}})=0\,, \\ \rho_{0} &\frac{\partial \tilde{\mathbf{U}}}{\partial t} - \tilde p\, \nabla = (\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} \,, \\ &\frac{\partial \tilde p}{\partial t} + \nabla\cdot(p_0 \tilde{\mathbf{U}}) + \frac{2}{3}\,p_0\nabla\cdot \tilde{\mathbf{U}}=0\,, \\ &\frac{\partial \tilde{\mathbf{B}}}{\partial t} - \nabla\times(\tilde{\mathbf{U}} \times \mathbf{B}) = 0\,, \end{aligned} \right. \\[2mm] \textnormal{EPs}\,\, &\left\{\,\, \begin{aligned} \quad &\frac{\partial f_\textnormal{h}}{\partial t} + \frac{1}{B_\parallel^*}(v_\parallel \mathbf{B}^* - \mathbf{b}_0 \times \mathbf{E}^*)\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\,, \\ & n_\textnormal{gc} = \int f_\textnormal{h} B_\parallel^* \,\textnormal dv_\parallel \textnormal d\mu \,, \\ & \mathbf{J}_\textnormal{gc} = \int \frac{f_\textnormal{h}}{B_\parallel^*}(v_\parallel \mathbf{B}^* - \mathbf{b}_0 \times \mathbf{E}^*) \,\textnormal dv_\parallel \textnormal d\mu \,, \\ & \mathbf{M}_\textnormal{gc} = - \int f_\textnormal{h} B_\parallel^* \mu \mathbf{b}_0 \,\textnormal dv_\parallel \textnormal d\mu \,, \end{aligned} \right. \end{align}\end{split}\]

where

\[\begin{split}\begin{align} 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{align}\end{split}\]

with the normalization parameter

\[\epsilon = \frac{1}{\hat \Omega_\textnormal{c,hot} \hat t} \,, \qquad \hat \Omega_\textnormal{c,hot} = \frac{Z_\textnormal{h} e \hat B}{A_\textnormal{h} m_\textnormal{H}} \,.\]

Propagators (called in sequence):

1. PushGuidingCenterBxEstar

2. PushGuidingCenterParallel

3. CurrentCoupling5DGradB

4. CurrentCoupling5DCurlb

5. CurrentCoupling5DDensity

6. ShearAlfvenCurrentCoupling5D

7. Magnetosonic

Model info:

classmethod doc_pde()

PDEs solved by model:

MHD continuity:

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

MHD momentum:

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

MHD pressure:

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

MHD induction:

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

Energetic-particle drift-kinetic equation:

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

Energetic-particle moments:

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

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

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

where

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

classmethod doc_normalization()

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

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

classmethod doc_scalar_quantities()

The following scalars are tracked during simulation:

• MHD kinetic energy: en_U

• Thermal pressure energy: en_p

• Magnetic energy: en_B

• Parallel energetic-particle energy: en_fv

• Magnetic-moment energetic-particle energy: en_fB

• Total energy: en_tot

• Lost particles: n_lost_particles

classmethod doc_discretization()

classmethod doc_long_description()

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

classmethod doc_examples()

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

from struphy.models import LinearMHDDriftkineticCC

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

classmethod doc_use_cases()

This model is appropriate for:

• linear energetic-ion effects with guiding-center reduction

• current-coupling hybrid mode studies in strong magnetic fields

• verification of 5D hybrid coupling operators

classmethod doc_cannot_be_used_for()

This model is not suitable for:

• full 6D energetic-particle dynamics

• nonlinear hybrid turbulence

• resistive or viscous MHD closures

• regimes where gyrophase resolution is required

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

Bases: StruphyModel

Hybrid linear MHD + energetic ions (6D Vlasov) with current coupling scheme.

Normalization:

\[\hat U = \hat v = \hat v_\textnormal{A} \,, \qquad \hat f_\textnormal{h} = \frac{\hat n}{\hat v_\textnormal{A}^3} \,.\]

Equations:

\[\begin{split}\begin{align} \textnormal{MHD}\,\, &\left\{\,\, \begin{aligned} &\frac{\partial \tilde{\rho}}{\partial t}+\nabla\cdot(\rho_0 \tilde{\mathbf{U}})=0\,, \\[2mm] \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}\,, \\[2mm] &\frac{\partial \tilde p}{\partial t} + (\gamma-1)\nabla\cdot(p_0 \tilde{\mathbf{U}}) + p_0\nabla\cdot \tilde{\mathbf{U}}=0\,, \\[2mm] &\frac{\partial \tilde{\mathbf{B}}}{\partial t} = \nabla\times(\tilde{\mathbf{U}} \times \mathbf{B}_0)\,,\qquad \nabla\cdot\tilde{\mathbf{B}}=0\,, \end{aligned} \right. \\[2mm] \textnormal{EPs}\,\, &\left\{\,\, \begin{aligned} &\quad\,\,\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\,, \\[2mm] &\quad\,\,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\,, \end{aligned} \right. \end{align}\end{split}\]

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

\[\varepsilon = \frac{1}{\hat \Omega_{\textnormal{c,hot}} \hat t}\,,\qquad \textnormal{with} \qquad\hat \Omega_{\textnormal{c,hot}} = \frac{Z_\textnormal{h}e \hat B}{A_\textnormal{h} m_\textnormal{H}}\,.\]

Propagators (called in sequence):

1. CurrentCoupling6DDensity

2. ShearAlfvenPropagator

3. CurrentCoupling6DCurrent

4. PushEta

5. PushVxB

6. Magnetosonic

Model info:

classmethod doc_pde()

PDEs solved by model:

MHD continuity:

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

MHD momentum:

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

MHD pressure:

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

MHD induction:

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

Energetic-particle Vlasov equation:

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

Energetic-particle moments:

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

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

classmethod doc_normalization()

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

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

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

classmethod doc_scalar_quantities()

The following scalars are tracked during simulation:

• MHD kinetic energy: en_U

• Thermal pressure energy: en_p

• Magnetic energy: en_B

• Energetic-particle energy: en_f

• Total energy: en_tot

• Lost particles: n_lost_particles

classmethod doc_discretization()

classmethod doc_long_description()

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

classmethod doc_examples()

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

from struphy.models import LinearMHDVlasovCC

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

classmethod doc_use_cases()

This model is appropriate for:

• linear energetic-particle coupling to MHD modes

• current-coupling hybrid verification

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

classmethod doc_cannot_be_used_for()

This model is not suitable for:

• strongly nonlinear dynamics far from equilibrium

• full multi-species kinetic plasma evolution

• pressure-coupling closures

• collisional or dissipative MHD effects not present in the equations

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

Bases: StruphyModel

Hybrid linear MHD + energetic ions (6D Vlasov) with pressure coupling scheme.

Normalization:

\[\hat U = \hat v =: \hat v_\textnormal{A, bulk} \,, \qquad \hat f_\textnormal{h} = \frac{\hat n}{\hat v_\textnormal{A}^3} \,,\qquad \hat{\mathbb{P}}_\textnormal{h} = A_\textnormal{h}m_\textnormal{H}\hat n \hat v_\textnormal{A}^2\,,\]

Implemented equations:

\[\begin{split}\begin{align} \textnormal{MHD} &\left\{ \begin{aligned} &\frac{\partial \tilde{\rho}}{\partial t}+\nabla\cdot(\rho_0 \tilde{\mathbf{U}})=0\,, \\ \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}} \,, \qquad \mathbf{J}_0 = \nabla\times\mathbf{B}_0\,, \\ &\frac{\partial \tilde p}{\partial t} + \nabla\cdot(p_0 \tilde{\mathbf{U}}) + \frac{2}{3}\,p_0\nabla\cdot \tilde{\mathbf{U}}=0\,, \\ &\frac{\partial \tilde{\mathbf{B}}}{\partial t} - \nabla\times(\tilde{\mathbf{U}} \times \mathbf{B}_0) = 0\,, \end{aligned} \right. \\[2mm] \textnormal{EPs}\,\, &\left\{\,\, \begin{aligned} &\quad\,\,\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\,, \\ &\quad\,\,\tilde{\mathbb{P}}_{\textnormal{h},\perp} = \int \mathbf{v}_\perp\mathbf{v}^\top_\perp f_\textnormal{h} d\mathbf{v} \,, \end{aligned} \right. \end{align}\end{split}\]

where

\[\epsilon = \frac{\hat \omega}{2 \pi \, \hat \Omega_{\textnormal{c,hot}}} \,,\qquad \textnormal{with} \qquad\hat \Omega_{\textnormal{c,hot}} = \frac{Z_\textnormal{h}e \hat B}{A_\textnormal{h} m_\textnormal{H}}\,.\]

Propagators (called in sequence):

1. PushEtaPC

2. PushVxB

3. PressureCoupling6D

4. ShearAlfvenPropagator

5. Magnetosonic

Model info:

classmethod doc_pde()

PDEs solved by model:

MHD continuity:

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

MHD momentum:

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

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

MHD pressure:

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

MHD induction:

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

Energetic-particle Vlasov equation:

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

Perpendicular pressure tensor:

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

classmethod doc_normalization()

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

classmethod doc_scalar_quantities()

The following scalars are tracked during simulation:

• MHD kinetic energy: en_U

• Thermal pressure energy: en_p

• Magnetic energy: en_B

• Energetic-particle kinetic energy: en_f

• Total energy: en_tot

• Lost particles: n_lost_particles

classmethod doc_discretization()

classmethod doc_long_description()

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

classmethod doc_examples()

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

from struphy.models import LinearMHDVlasovPC

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

classmethod doc_use_cases()

This model is appropriate for:

• linear pressure-coupling hybrid studies

• energetic-particle pressure feedback on MHD modes

• verification of PushEtaPC and pressure-coupling operators

classmethod doc_cannot_be_used_for()

This model is not suitable for:

• current-coupling closure studies

• nonlinear hybrid turbulence

• dissipative/resistive MHD

• fully kinetic treatment of the bulk plasma

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

Bases: StruphyModel

Linearized Vlasov-Ampère equations for one species.

Normalization:

\[\begin{align} \hat v = c \,, \qquad \hat E = \hat B \hat v\,,\qquad \hat \phi = \hat E \hat x \,. \end{align}\]

Equations:

\[\begin{split}\begin{align} & \frac{\partial \tilde{\mathbf E}}{\partial t} = - \frac{\alpha^2}{\varepsilon} \int_{\mathbb R^3} \mathbf{v} \tilde f\, \textrm d^3 \mathbf v \,, \\[2mm] & \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 \,, \end{align}\end{split}\]

with the normalization parameter

\[\alpha = \frac{\hat \Omega_\textnormal{p}}{\hat \Omega_\textnormal{c}}\,,\qquad \varepsilon = \frac{1}{\hat \Omega_\textnormal{c} \hat t} \,,\qquad \textnormal{with} \qquad \hat\Omega_\textnormal{p} = \sqrt{\frac{\hat n (Ze)^2}{\epsilon_0 (A m_\textnormal{H})}} \,,\qquad \hat \Omega_{\textnormal{c}} = \frac{(Ze) \hat B}{(A m_\textnormal{H})}\,,\]

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}\begin{align} \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{align}\end{split}\]

Moreover, it is assumed that

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

Propagators (called in sequence):

1. PushEta

2. PushVinEfield

3. EfieldWeightsCoupling

4. PushVxB

Model info:

classmethod doc_pde()

PDEs solved by model:

Linearized Ampère’s law:

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

Linearized Vlasov equation:

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

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

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

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

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

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

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

Moreover, it is assumed that

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

classmethod doc_normalization()

The light speed defines the velocity scale:

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

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

classmethod doc_scalar_quantities()

The following scalars are tracked during simulation:

• Electric field energy: en_E

• Perturbation particle energy: en_w

• Total energy: en_tot

classmethod doc_discretization()

classmethod doc_long_description()

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

classmethod doc_examples()

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

from struphy.models import LinearVlasovAmpereOneSpecies

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

classmethod doc_use_cases()

This model is appropriate for:

• linear delta-f kinetic instabilities

• weakly perturbed Landau-like damping studies with electrostatic fields

• verification of linearized field-particle coupling

classmethod doc_cannot_be_used_for()

This model is not suitable for:

• strongly nonlinear departures from the chosen equilibrium

• multi-species kinetic coupling

• fully electromagnetic magnetic-field evolution

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

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

Bases: LinearVlasovAmpereOneSpecies

Linearized Vlasov-Ampère equations for one species.

Normalization:

\[\begin{align} \hat v = c \,, \qquad \hat E = \hat B \hat v\,,\qquad \hat \phi = \hat E \hat x \,. \end{align}\]

Equations:

\[\begin{split}\begin{align} & \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 \,, \\[2mm] & \frac{\partial \tilde{\mathbf B}}{\partial t} = - \nabla \times \tilde{\mathbf E} \,, \\[2mm] & \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 \,, \end{align}\end{split}\]

with the normalization parameter

\[\alpha = \frac{\hat \Omega_\textnormal{p}}{\hat \Omega_\textnormal{c}}\,,\qquad \varepsilon = \frac{1}{\hat \Omega_\textnormal{c} \hat t} \,,\qquad \textnormal{with} \qquad \hat\Omega_\textnormal{p} = \sqrt{\frac{\hat n (Ze)^2}{\epsilon_0 (A m_\textnormal{H})}} \,,\qquad \hat \Omega_{\textnormal{c}} = \frac{(Ze) \hat B}{(A m_\textnormal{H})}\,,\]

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}\begin{align} \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{align}\end{split}\]

Moreover, it is assumed that

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

Propagators (called in sequence):

1. PushEta

2. PushVinEfield

3. EfieldWeightsCoupling

4. PushVxB

5. Maxwell

Model info:

classmethod doc_pde()

PDEs solved by model:

Linearized Ampère’s law:

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

Linearized Faraday’s law:

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

Linearized Vlasov equation:

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

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

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

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

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

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

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

Moreover, it is assumed that

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

classmethod doc_normalization()

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

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

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

classmethod doc_scalar_quantities()

The following scalars are tracked during simulation:

• Electric field energy: en_E

• Magnetic field energy: en_B

• Perturbation particle energy: en_w

• Total energy: en_tot

classmethod doc_discretization()

classmethod doc_long_description()

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

classmethod doc_examples()

Create and initialize the linear Vlasov-Maxwell model:

from struphy.models import LinearVlasovMaxwellOneSpecies

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

classmethod doc_use_cases()

This model is appropriate for:

• linear electromagnetic kinetic wave studies

• delta-f verification of linear Vlasov-Maxwell coupling

• weakly nonlinear regimes where equilibrium perturbations stay small

classmethod doc_cannot_be_used_for()

This model is not suitable for:

• strongly nonlinear fully kinetic plasma dynamics

• multi-species electromagnetic coupling

• problems without a well-defined linearization background

• collisional kinetic closures

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

Bases: StruphyModel

Weak discretization of Poisson’s equation with diffusion matrix, stabilization and time-depedent right-hand side.

Normalization:

\[\hat D = \frac{\hat n}{\hat x^2}\,,\qquad \hat \rho = \hat n \,.\]

Equations: 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)\) parametrized with 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.

Propagators (called in sequence):

1. TimeDependentSource

2. ImplicitDiffusion

Model info:

classmethod doc_pde()

PDEs solved by model:

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

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

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

classmethod doc_normalization()

The coefficient scaling is

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

No dedicated velocity normalization is used.

classmethod doc_scalar_quantities()

The following scalars are tracked during simulation:

• No default scalar diagnostics are defined by this model.

classmethod doc_discretization()

classmethod doc_long_description()

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

classmethod doc_examples()

Create and initialize a Poisson model:

from struphy.models import Poisson

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

classmethod doc_use_cases()

This model is appropriate for:

• elliptic benchmark problems

• electrostatic field solves with prescribed source terms

• testing weak Poisson discretizations and boundary handling

classmethod doc_cannot_be_used_for()

This model is not suitable for:

• hyperbolic time-dependent wave propagation

• self-consistent kinetic plasma evolution on its own

• magnetic-field dynamics or full Maxwell coupling

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

Bases: StruphyModel

Pressureless fluid discretized with smoothed particle hydrodynamics

Equations:

\[\begin{split}&\partial_t \rho + \nabla \cdot ( \rho \mathbf u ) = 0 \,, \\[4mm] &\partial_t (\rho \mathbf u) + \nabla \cdot (\rho \mathbf u \otimes \mathbf u) = - \nabla \phi_0 \,,\end{split}\]

where \(\phi_0\) is a static external potential.

Propagators (called in sequence):

1. PushEta

This is discretized by particles going in straight lines.

classmethod doc_pde()

PDEs solved by model:

Continuity:

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

Momentum:

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

where \(\phi_0\) is a static external potential.

classmethod doc_normalization()

No special field normalization is introduced beyond the particle units inherited from the simulation setup. This model does not define a separate wave or plasma velocity scale.

classmethod doc_scalar_quantities()

The following scalars are tracked during simulation:

• SPH kinetic energy: en_kin

classmethod doc_discretization()

classmethod doc_long_description()

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

classmethod doc_examples()

Create and initialize a pressureless SPH model:

from struphy.models import PressureLessSPH

model = PressureLessSPH()
model.cold_fluid.var

classmethod doc_use_cases()

This model is appropriate for:

• meshfree pressureless-flow benchmarks

• straight-line particle transport with external forcing

• testing SPH particle loading and diagnostics

classmethod doc_cannot_be_used_for()

This model is not suitable for:

• compressible fluids with pressure forces

• FEEC-based grid discretizations

• electromagnetic plasma dynamics

• viscous or thermal closures

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

Bases: StruphyModel

Diffusion equation discretized with a (random) particle method; the diffusion is computed through a Wiener process.

Normalization:

\[\hat D := \frac{\hat x^2}{\hat t } \,.\]

Equations: 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.

Propagators (called in sequence):

1. PushRandomDiffusion

Model info:

classmethod doc_pde()

PDEs solved by model:

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

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

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

classmethod doc_normalization()

The natural scaling is set by the diffusion coefficient:

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

classmethod doc_scalar_quantities()

The following scalars are tracked during simulation:

• No default scalar diagnostics are defined by this model.

classmethod doc_discretization()

classmethod doc_long_description()

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

classmethod doc_examples()

Create and initialize a random diffusion model:

from struphy.models import RandomParticleDiffusion

model = RandomParticleDiffusion()
model.hydrogen.var

classmethod doc_use_cases()

This model is appropriate for:

• stochastic particle diffusion benchmarks

• Monte-Carlo transport verification

• comparison against deterministic diffusion solvers

classmethod doc_cannot_be_used_for()

This model is not suitable for:

• electromagnetic or fluid plasma dynamics

• deterministic advection-dominated transport

• anisotropic plasma kinetics in phase space

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

Bases: StruphyModel

ShearAlfven propagator from LinearMHD with zero-flow equilibrium (\(\mathbf U_0 = 0\)).

Normalization:

\[\hat U = \hat v_\textnormal{A} \,.\]

Equations:

\[ \begin{align}\begin{aligned}\rho_0&\frac{\partial \tilde{\mathbf{U}}}{\partial t} = (\nabla\times \tilde{\mathbf{B}})\times\mathbf{B}_0\,,\\&\frac{\partial \tilde{\mathbf{B}}}{\partial t} - \nabla\times(\tilde{\mathbf{U}} \times \mathbf{B}_0) = 0\,.\end{aligned}\end{align} \]

Propagators (called in sequence):

1. ShearAlfvenPropagator

Model info:

classmethod doc_pde()

PDEs solved by model:

Momentum:

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

Induction:

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

classmethod doc_normalization()

All velocities are normalized by the Alfvén speed,

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

classmethod doc_scalar_quantities()

The following scalars are tracked during simulation:

• Perturbed kinetic energy: en_U

• Perturbed magnetic energy: en_B

• Background magnetic energy: en_B_eq

• Total magnetic energy including background: en_B_tot

• Perturbation total energy: en_tot

• Total energy including background contribution: en_tot2

classmethod doc_discretization()

classmethod doc_long_description()

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

classmethod doc_examples()

Create and initialize a shear-Alfvén model:

from struphy.models import ShearAlfven

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

classmethod doc_use_cases()

This model is appropriate for:

• isolated shear-Alfvén wave propagation

• verification of the ShearAlfven propagator

• reduced linear-MHD studies without compressibility

classmethod doc_cannot_be_used_for()

This model is not suitable for:

• magnetosonic or compressible wave physics

• nonlinear MHD dynamics

• kinetic particle coupling

• resistive, viscous, or Hall-MHD effects

class struphy.models.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)

Bases: StruphyModel

Drift-kinetic equation for one ion species in static background magnetic field.

Normalization:

\[\hat v = \hat v_\textrm{i} = \sqrt{\frac{k_B \hat T_\textrm{i}}{m_\textrm{i}}}\,,\qquad \hat E = \hat v_\textrm{i}\hat B\,,\qquad \hat \phi = \hat E \hat x \,.\]

Equations:

\[\begin{split}&\frac{\partial f}{\partial t} + \frac{\mathbf{E} \times \mathbf{b}_0}{B^*_\parallel} \cdot \frac{\partial f}{\partial \mathbf{X}} = 0\,. \\[2mm] - &\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 \,.\end{split}\]

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

and with the normalization parameters

\[\varepsilon := \frac{1}{\hat \Omega_\textrm{c} \hat t}\,,\qquad \hat \Omega_\textrm{c} = \frac{q_\textrm{i} \hat B}{m_\textrm{i}} \,.\]

Notes

• The Control variate method in the Poisson equation is optional; in case it is enabled via the parameter file, the following Poisson equation is solved:

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

\[\int \frac{n_0}{|B_0|^2} \nabla_\perp \psi \cdot \nabla_\perp \phi\,\textrm d \mathbf x + \frac{1}{Z\varepsilon^2} \int \frac{n_0}{T_{0}} \psi \phi \,\textrm d \mathbf x = \frac 1 \varepsilon \int \int \psi \, (f - f_0) B^*_\parallel \,\textrm d \mathbf x\,\textnormal d v_\parallel \textnormal d \mu \qquad \forall \ \psi \in H^1\,.\]

Propagators (called in sequence):

1. ImplicitDiffusion

2. PushGuidingCenterBxEstar

Model info:

classmethod doc_pde()

PDEs solved by model:

Drift equation:

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

Poisson equation:

\[-\nabla_\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, \qquad \mathbf{B}^* = \mathbf{B}_0 + \varepsilon v_\parallel \nabla \times \mathbf{b}_0, \qquad B^*_\parallel = \mathbf{B}^* \cdot \mathbf{b}_0\]

Notes

• The control_var in the Poisson equation is optional; in case it is enabled via the parameter file, the following Poisson equation is solved:

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

\[\int \frac{n_0}{|B_0|^2} \nabla_\perp \psi \cdot \nabla_\perp \phi \, \textrm{d} \mathbf{x} + \frac{1}{Z \varepsilon^2} \int \frac{n_0}{T_0} \psi \phi \, \textrm{d} \mathbf{x} = \frac{1}{\varepsilon} \int \int \psi \, (f - f_0) B^*_\parallel \, \textrm{d} \mathbf{x} \, \textnormal{d} v_\parallel \textnormal{d} \mu \qquad \forall \ \psi \in H^1\]

classmethod doc_normalization()

The reference speed is the ion thermal velocity:

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

classmethod doc_scalar_quantities()

The following scalars are tracked during simulation:

• Electrostatic field energy: en_phi

• Particle kinetic energy: en_particles

• Total energy: en_tot

classmethod doc_discretization()

classmethod doc_long_description()

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

classmethod doc_examples()

Create and initialize the toy drift model:

from struphy.models import ToyDrift

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

classmethod doc_use_cases()

This model is appropriate for:

• reduced electrostatic guiding-center benchmarks

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

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

classmethod doc_cannot_be_used_for()

This model is not suitable for:

• full drift-kinetic dynamics with parallel streaming

• electromagnetic perturbations

• high-fidelity turbulence studies

• full-orbit kinetic physics

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

Bases: StruphyModel

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

Normalization:

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

Equations:

\[\begin{split}\frac{\partial \mathbf u}{\partial t} &= - \nabla \phi + \frac{\mathbf u \times \mathbf B_0}{\varepsilon} + \nu \Delta \mathbf u + \mathbf f\,, \\[2mm] 0 &= \nabla \phi - \frac{\mathbf u_e \times \mathbf B_0}{\varepsilon} + \nu_e \Delta \mathbf u_e + \mathbf f_e \,, \\[3mm] \nabla & \cdot (\mathbf u - \mathbf u_e) = 0\,,\end{split}\]

where \(\mathbf B_0\) is a static magnetic field and \(\mathbf f, \mathbf f_e\) are given forcing terms, and with the normalization parameter

\[\varepsilon = \frac{1}{\hat \Omega_\textnormal{c} \hat t} \,,\qquad \textnormal{with} \,,\qquad \hat \Omega_{\textnormal{c}} = \frac{(Ze) \hat B}{(A m_\textnormal{H})}\,,\]

Propagators (called in sequence):

1. TwoFluidQuasiNeutralFull

Model info:

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

PDEs solved by model:

Ion momentum:

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

Electron momentum:

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

Quasi-neutrality constraint:

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

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

classmethod doc_normalization()

Thermal-speed scaling is used:

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

classmethod doc_scalar_quantities()

The following scalars are tracked during simulation:

• No default scalar diagnostics are defined by this model.

classmethod doc_discretization()

classmethod doc_long_description()

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

classmethod doc_examples()

Create and initialize the quasi-neutral toy model:

from struphy.models import TwoFluidQuasiNeutralToy

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

classmethod doc_use_cases()

This model is appropriate for:

• linear quasi-neutral two-fluid benchmarks

• Stokes-like plasma model verification

• testing coupled velocity-potential FEEC solvers

classmethod doc_cannot_be_used_for()

This model is not suitable for:

• nonlinear two-fluid dynamics

• finite electron inertia effects

• kinetic phase-space phenomena

• self-consistent electromagnetic wave propagation

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

Bases: StruphyModel

Barotropic fluid equations discretized with a variational method.

Normalization:

\[\hat u = \hat v_\textnormal{A} \qquad \hat{\mathcal U} = \frac{\hat \rho}{2} \,.\]

Equations:

\[\begin{split}&\partial_t \rho + \nabla \cdot ( \rho \mathbf u ) = 0 \,, \\[4mm] &\partial_t (\rho \mathbf u) + \nabla \cdot (\rho \mathbf u \otimes \mathbf u) + \rho \nabla \frac{(\rho \mathcal U (\rho))}{\partial \rho} = 0 \,.\end{split}\]

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

Propagators (called in sequence):

1. VariationalDensityEvolve

2. VariationalMomentumAdvection

Model info:

classmethod doc_pde()

PDEs solved by model:

Continuity:

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

Momentum:

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

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

classmethod doc_normalization()

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

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

classmethod doc_scalar_quantities()

The following scalars are tracked during simulation:

• Kinetic energy: kinetic_energy

• Thermodynamic energy: thermo_energy

• Total energy: total_energy

classmethod doc_discretization()

classmethod doc_long_description()

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

classmethod doc_examples()

Create and initialize a variational barotropic-fluid model:

from struphy.models import VariationalBarotropicFluid

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

classmethod doc_use_cases()

This model is appropriate for:

• barotropic compressible-flow benchmarks

• testing variational density and momentum propagators

• reduced fluid studies without entropy evolution

classmethod doc_cannot_be_used_for()

This model is not suitable for:

• non-barotropic thermodynamics

• magnetic-field dynamics

• viscous or resistive effects

• kinetic plasma phenomena

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

Bases: StruphyModel

Fully compressible fluid equations discretized with a variational method.

Normalization:

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

Equations:

\[\begin{split}&\partial_t \rho + \nabla \cdot ( \rho \mathbf u ) = 0 \,, \\[4mm] &\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 \,, \\[4mm] &\partial_t s + \nabla \cdot ( s \mathbf u ) = 0 \,,\end{split}\]

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

Propagators (called in sequence):

1. VariationalDensityEvolve

2. VariationalMomentumAdvection

3. VariationalEntropyEvolve

Model info:

classmethod doc_pde()

PDEs solved by model:

Continuity:

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

Momentum:

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

Entropy:

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

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

classmethod doc_normalization()

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

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

classmethod doc_scalar_quantities()

The following scalars are tracked during simulation:

• Kinetic energy: en_U

• Thermodynamic energy: en_thermo

• Total energy: en_tot

classmethod doc_discretization()

classmethod doc_long_description()

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

classmethod doc_examples()

Create and initialize a variational compressible-fluid model:

from struphy.models import VariationalCompressibleFluid

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

classmethod doc_use_cases()

This model is appropriate for:

• compressible fluid benchmarks with entropy transport

• testing the variational fluid propagator stack

• non-magnetic hydrodynamics with conservative thermodynamics

classmethod doc_cannot_be_used_for()

This model is not suitable for:

• magnetic-field evolution or MHD coupling

• pressureless or barotropic-only reductions

• kinetic or particle-based transport physics

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

Bases: StruphyModel

Pressure-less fluid equations discretized with a variational method.

Normalization:

\[\hat u = \hat v_\textnormal{A} \,.\]

Equations:

\[\begin{split}&\partial_t \rho + \nabla \cdot ( \rho \mathbf u ) = 0 \,, \\[4mm] &\partial_t (\rho \mathbf u) + \nabla \cdot (\rho \mathbf u \otimes \mathbf u) = 0 \,.\end{split}\]

Propagators (called in sequence):

1. VariationalDensityEvolve

2. VariationalMomentumAdvection

Model info:

classmethod doc_pde()

PDEs solved by model:

Continuity:

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

Momentum:

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

classmethod doc_normalization()

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

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

classmethod doc_scalar_quantities()

The following scalars are tracked during simulation:

• Kinetic energy: kinetic_energy

classmethod doc_discretization()

classmethod doc_long_description()

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

classmethod doc_examples()

Create and initialize a pressureless variational-fluid model:

from struphy.models import VariationalPressurelessFluid

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

classmethod doc_use_cases()

This model is appropriate for:

• pressureless compressible benchmarks

• testing the minimal variational fluid update chain

• reduced transport problems without thermodynamics

classmethod doc_cannot_be_used_for()

This model is not suitable for:

• pressure- or entropy-driven flow

• magnetic-field coupling

• viscous/resistive dissipation

• kinetic particle physics

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

Bases: StruphyModel

\(\delta f\) visco-resistive MHD equations discretized with a variational method.

Normalization:

\[\hat u = \hat v_\textnormal{A}\,.\]

Equations:

\[\begin{split}&\partial_t \tilde{\rho} + \nabla \cdot ( (\tilde{\rho}+\rho_0) \tilde{\mathbf u} ) = 0 \,, \\[4mm] &\partial_t ((\tilde{\rho}+\rho_0) \tilde{\mathbf u}) + \nabla \cdot ((\tilde{\rho}+\rho_0) \tilde{\mathbf u} \otimes \tilde{\mathbf u}) + \frac{1}{\gamma -1} \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 \,, \\[4mm] &\partial_t \tilde{p} + \tilde{\mathbf u} \cdot \nabla (\tilde{p} + p_0) + \gamma (\tilde{p} + p_0) \nabla \cdot \tilde{\mathbf u} = \frac{1}{(\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) \,, \\[4mm] &\partial_t \tilde{\mathbf B} + \nabla \times ( (\tilde{\mathbf B} + \mathbf B_0) \times \tilde{\mathbf u} ) + \nabla \times (\eta + \eta_a(\mathbf x)) \nabla \times \tilde{\mathbf B} = 0 \,,\end{split}\]

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

Propagators (called in sequence):

1. VariationalDensityEvolve

2. VariationalMomentumAdvection

3. VariationalPBEvolve

4. VariationalViscosity

5. VariationalResistivity

Model info:

classmethod doc_pde()

PDEs solved by model:

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) + \frac{1}{\gamma - 1} \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}} = \frac{1}{\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.

classmethod doc_normalization()

The velocity scale is Alfvénic:

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

classmethod doc_scalar_quantities()

The following scalars are tracked during simulation:

• Kinetic energy: en_U

• Thermal energies: en_thermo, en_thermo_l1

• Magnetic energies: en_mag_1, en_mag_2, en_mag_l1

• Total energies: en_tot, en_tot_l1

classmethod doc_discretization()

classmethod doc_long_description()

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

classmethod doc_examples()

Create and initialize the visco-resistive delta f MHD model:

from struphy.models import ViscoResistiveDeltafMHD

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

classmethod doc_use_cases()

This model is appropriate for:

• perturbative nonlinear MHD around a prescribed equilibrium

• dissipative MHD verification with background fields

• comparing linear and delta f variational formulations

classmethod doc_cannot_be_used_for()

This model is not suitable for:

• fully full-f MHD evolution far from the chosen equilibrium

• kinetic or particle-resolved plasma effects

• entropy- or q-based thermodynamic formulations

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

Bases: StruphyModel

Linear visco-resistive MHD equations discretized with a variational method.

Normalization:

\[\hat u = \hat v_\textnormal{A}\,.\]

Equations:

\[\begin{split}&\partial_t \tilde{\rho} + \nabla \cdot ( \rho_0 \tilde{\mathbf u} ) = 0 \,, \\[4mm] &\partial_t (\rho_0 \tilde{\mathbf u}) + \frac{2 q_0}{\gamma -1} \nabla \tilde{q} + \frac{2 \tilde{q}}{\gamma -1} \nabla q_0 + \frac{2 \tilde{q}}{\gamma -1} \nabla \tilde{q} + \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 \,, \\[4mm] &\partial_t \tilde{q} + \cdot(\nabla (q_0 + \tilde{q}) \mathbf u) + (\gamma/2 -1) (q_0 + \tilde{q}) \nabla \cdot u = 0 \,, \\[4mm] &\partial_t \tilde{\mathbf B} + \nabla \times ( \mathbf (B_0 + \tilde{\mathbf B}) \times \tilde{\mathbf u} ) + \nabla \times (\eta + \eta_a(\mathbf x)) \nabla \times \tilde{\mathbf B} = 0 \,,\end{split}\]

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

Propagators (called in sequence):

1. VariationalDensityEvolve

2. VariationalMomentumAdvection

3. VariationalQBEvolve

4. VariationalViscosity

5. VariationalResistivity

Model info:

classmethod doc_pde()

PDEs solved by model:

Continuity:

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

Momentum:

\[\partial_t (\rho_0 \tilde{\mathbf{u}}) + \frac{2 q_0}{\gamma - 1} \nabla \tilde{q} + \frac{2 \tilde{q}}{\gamma - 1} \nabla q_0 + \frac{2 \tilde{q}}{\gamma - 1} \nabla \tilde{q} + \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\]

Energy-like variable:

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

Induction:

\[\partial_t \tilde{\mathbf{B}} + \nabla \times \left( (\mathbf{B}_0 + \tilde{\mathbf{B}}) \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.

classmethod doc_normalization()

The model uses Alfvén-speed normalization:

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

classmethod doc_scalar_quantities()

The following scalars are tracked during simulation:

• Kinetic energy: en_U

• Magnetic energies: en_mag_1, en_mag_2

• Thermodynamic q-energies: en_thermo_1, en_thermo_2

• Total energy: en_tot

classmethod doc_discretization()

classmethod doc_long_description()

This model is the delta f dissipative MHD formulation in the q-variable representation. It is aimed at perturbative nonlinear MHD studies where the q-formulation is preferred numerically.

classmethod doc_examples()

Create and initialize the visco-resistive delta f q-MHD model:

from struphy.models import ViscoResistiveDeltafMHD_with_q

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

classmethod doc_use_cases()

This model is appropriate for:

• dissipative delta f MHD with q-thermodynamics

• comparison against the p-based delta f model

• nonlinear perturbative benchmarks with viscosity and resistivity

classmethod doc_cannot_be_used_for()

This model is not suitable for:

• full-f MHD far from equilibrium

• entropy-based thermodynamic evolution

• kinetic or hybrid particle coupling

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

Bases: StruphyModel

Linear visco-resistive MHD equations discretized with a variational method.

Normalization:

\[\hat u = \hat v_\textnormal{A}\,.\]

Equations:

\[\begin{split}&\partial_t \tilde{\rho} + \nabla \cdot ( \rho_0 \tilde{\mathbf u} ) = 0 \,, \\[4mm] &\partial_t (\rho_0 \tilde{\mathbf u}) + \frac{1}{\gamma -1} \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 \,, \\[4mm] &\partial_t \tilde{p} + \tilde{\mathbf u} \cdot \nabla p_0 + \gamma p_0 \nabla \cdot \tilde{\mathbf u} = \frac{1}{(\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) \,, \\[4mm] &\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}\]

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

Propagators (called in sequence):

1. VariationalDensityEvolve

2. VariationalPBEvolve

3. VariationalViscosity

4. VariationalResistivity

Model info:

classmethod doc_pde()

PDEs solved by model:

Continuity:

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

Momentum:

\[\partial_t (\rho_0 \tilde{\mathbf{u}}) + \frac{1}{\gamma - 1} \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\]

Pressure:

\[\partial_t \tilde{p} + \tilde{\mathbf{u}} \cdot \nabla p_0 + \gamma p_0 \nabla \cdot \tilde{\mathbf{u}} = \frac{1}{\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 (\mathbf{B}_0 \times \tilde{\mathbf{u}}) + \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.

classmethod doc_normalization()

The characteristic velocity is the Alfvén speed,

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

classmethod doc_scalar_quantities()

The following scalars are tracked during simulation:

• Kinetic energy: en_U

• Thermal perturbation energy: en_thermo

• Magnetic energies: en_mag_1, en_mag_2

• Total energy diagnostics: en_tot, en_tot_l1

• Auxiliary linearized diagnostics: en_thermo_l1, en_mag_l1

classmethod doc_discretization()

classmethod doc_long_description()

This is the linear dissipative MHD model with pressure as primitive thermodynamic variable. It is intended for small-amplitude perturbations around an equilibrium while retaining explicit viscosity and resistivity.

classmethod doc_examples()

Create and initialize a linear visco-resistive MHD model:

from struphy.models import ViscoResistiveLinearMHD

model = ViscoResistiveLinearMHD()
model.em_fields.b_field
model.mhd.velocity
model.mhd.pressure

classmethod doc_use_cases()

This model is appropriate for:

• linear dissipative MHD wave studies

• equilibrium perturbation benchmarks with viscosity and resistivity

• verification of linear variational MHD operators

classmethod doc_cannot_be_used_for()

This model is not suitable for:

• strongly nonlinear MHD dynamics

• fully kinetic plasma effects

• entropy- or q-based thermodynamic formulations

• ideal MHD benchmarks where dissipation must be absent by construction

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

Bases: StruphyModel

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

Normalization:

\[\hat u = \hat v_\textnormal{A}\,.\]

Equations:

\[\begin{split}&\partial_t \tilde{\rho} + \nabla \cdot ( \rho_0 \tilde{\mathbf u} ) = 0 \,, \\[4mm] &\partial_t (\rho_0 \tilde{\mathbf u}) + \frac{2 q_0}{\gamma -1} \nabla \tilde{q} + \frac{2 \tilde{q}}{\gamma -1} \nabla q_0 + \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 \,, \\[4mm] &\partial_t \tilde{q} + \cdot(\nabla q_0 \mathbf u) + (\gamma/2 -1) q_0 \nabla \cdot u = 0 \,, \\[4mm] &\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}\]

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

Propagators (called in sequence):

1. VariationalDensityEvolve

2. VariationalQBEvolve

3. VariationalViscosity

4. VariationalResistivity

Model info:

classmethod doc_pde()

PDEs solved by model:

Continuity:

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

Momentum:

\[\partial_t (\rho_0 \tilde{\mathbf{u}}) + \frac{2 q_0}{\gamma - 1} \nabla \tilde{q} + \frac{2 \tilde{q}}{\gamma - 1} \nabla q_0 + \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\]

Energy-like variable:

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

Induction:

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

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

classmethod doc_normalization()

The flow normalization is Alfvénic:

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

classmethod doc_scalar_quantities()

The following scalars are tracked during simulation:

• Kinetic energy: en_U

• Magnetic energies: en_mag_1, en_mag_2

• Thermodynamic q-energies: en_thermo_1, en_thermo_2

• Total energy: en_tot

classmethod doc_discretization()

classmethod doc_long_description()

This variant of the linear dissipative MHD model uses the square root of the pressure as primary thermodynamic variable, which is convenient for the corresponding variational discretization.

classmethod doc_examples()

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

from struphy.models import ViscoResistiveLinearMHD_with_q

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

classmethod doc_use_cases()

This model is appropriate for:

• linear dissipative MHD in the q-formulation

• comparing p- and q-based variational discretizations

• verification of linear q/B propagators

classmethod doc_cannot_be_used_for()

This model is not suitable for:

• nonlinear MHD dynamics

• entropy-based thermodynamics

• kinetic plasma coupling

• ideal nondissipative studies

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

Bases: StruphyModel

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

Normalization:

\[\hat u = \hat v_\textnormal{A}\,, \qquad \hat{\mathcal U} = \frac{\hat{\mathbf B}^2}{\hat \rho \mu_0 (\gamma-1)} \,,\qquad \hat s = \hat \rho\ \textrm{ln}\left(\frac{\hat{\mathbf B}^2}{\mu_0 (\gamma -1) \hat{\rho}}\right) \,.\]

Equations:

\[\begin{split}&\partial_t \rho + \nabla \cdot ( \rho \mathbf u ) = 0 \,, \\[4mm] &\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 \,, \\[4mm] &\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) \,, \\[4mm] &\partial_t \mathbf B + \nabla \times ( \mathbf B \times \mathbf u ) + \nabla \times (\eta + \eta_a(\mathbf x)) \nabla \times \mathbf B = 0 \,,\end{split}\]

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.

Propagators (called in sequence):

1. VariationalDensityEvolve

2. VariationalMomentumAdvection

3. VariationalEntropyEvolve

4. VariationalMagFieldEvolve

5. VariationalViscosity

6. VariationalResistivity

Model info:

classmethod doc_pde()

PDEs solved by model:

Continuity:

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

Momentum:

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

Entropy:

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

Induction:

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

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

classmethod doc_normalization()

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

classmethod doc_scalar_quantities()

The following scalars are tracked during simulation:

• Kinetic energy: en_U

• Thermodynamic energy: en_thermo

• Magnetic energy: en_mag

• Total energy: en_tot

• Total density / entropy: dens_tot, entr_tot

• Divergence diagnostic: tot_div_B

classmethod doc_discretization()

classmethod doc_long_description()

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

classmethod doc_examples()

Create and initialize the full visco-resistive MHD model:

from struphy.models import ViscoResistiveMHD

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

classmethod doc_use_cases()

This model is appropriate for:

• nonlinear dissipative MHD simulations

• entropy-based variational MHD benchmarks

• studies where both viscosity and resistivity matter

classmethod doc_cannot_be_used_for()

This model is not suitable for:

• ideal nondissipative MHD if strict conservation is required

• kinetic or hybrid particle effects

• reduced linear perturbation studies better served by the linear models

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

Bases: StruphyModel

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

Normalization:

\[\hat u = \hat v_\textnormal{A}\,.\]

Equations:

\[\begin{split}&\partial_t \rho + \nabla \cdot ( \rho \mathbf u ) = 0 \,, \\[4mm] &\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 \,, \\[4mm] &\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) \,, \\[4mm] &\partial_t \mathbf B + \nabla \times ( \mathbf B \times \mathbf u ) + \nabla \times (\eta + \eta_a(\mathbf x)) \nabla \times \mathbf B = 0 \,,\end{split}\]

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

Propagators (called in sequence):

1. VariationalDensityEvolve

2. VariationalMomentumAdvection

3. VariationalPBEvolve

4. VariationalViscosity

5. VariationalResistivity

Model info:

classmethod doc_pde()

PDEs solved by model:

Continuity:

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

Momentum:

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

Pressure:

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

Induction:

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

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

classmethod doc_normalization()

The characteristic flow speed is the Alfvén speed,

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

classmethod doc_scalar_quantities()

The following scalars are tracked during simulation:

• Kinetic energy: en_U

• Thermodynamic energy: en_thermo

• Magnetic energy: en_mag

• Total energy: en_tot

• Total density: dens_tot

• Divergence diagnostic: tot_div_B

classmethod doc_discretization()

classmethod doc_long_description()

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

classmethod doc_examples()

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

from struphy.models import ViscoResistiveMHD_with_p

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

classmethod doc_use_cases()

This model is appropriate for:

• full nonlinear dissipative MHD in p-formulation

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

• variational resistive/viscous MHD benchmarks

classmethod doc_cannot_be_used_for()

This model is not suitable for:

• ideal MHD without dissipation

• kinetic plasma physics

• reduced linear perturbation studies

• thermodynamic formulations that require entropy or q as primitive variables

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

Bases: StruphyModel

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

Normalization:

\[\hat u = \hat v_\textnormal{A}\,.\]

Equations:

\[\begin{split}&\partial_t \rho + \nabla \cdot ( \rho \mathbf u ) = 0 \,, \\[4mm] &\partial_t (\rho \mathbf u) + \nabla \cdot (\rho \mathbf u \otimes \mathbf u) + \frac{2q}{\gamma -1} \nabla q + \mathbf B \times \nabla \times \mathbf B - \nabla \cdot \left((\mu+\mu_a(\mathbf x)) \nabla \mathbf u \right) = 0 \,, \\[4mm] &\partial_t q + \cdot(\nabla q \mathbf u) + (\gamma/2 -1) q \nabla \cdot 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) \,, \\[4mm] &\partial_t \mathbf B + \nabla \times ( \mathbf B \times \mathbf u ) + \nabla \times (\eta + \eta_a(\mathbf x)) \nabla \times \mathbf B = 0 \,,\end{split}\]

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

Propagators (called in sequence):

1. VariationalDensityEvolve

1. VariationalDensityEvolve

2. VariationalMomentumAdvection

3. VariationalQBEvolve

4. VariationalViscosity

5. VariationalResistivity

Model info:

classmethod doc_pde()

PDEs solved by model:

Continuity:

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

Momentum:

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

Energy-like variable:

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

Induction:

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

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

classmethod doc_normalization()

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

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

classmethod doc_scalar_quantities()

The following scalars are tracked during simulation:

• Kinetic energy: en_U

• Thermodynamic q-energy: en_thermo

• Magnetic energy: en_mag

• Total energy: en_tot

• Total density: dens_tot

• Divergence diagnostic: tot_div_B

classmethod doc_discretization()

classmethod doc_long_description()

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

classmethod doc_examples()

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

from struphy.models import ViscoResistiveMHD_with_q

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

classmethod doc_use_cases()

This model is appropriate for:

• full nonlinear dissipative MHD in q-form

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

• variational MHD studies where q is numerically advantageous

classmethod doc_cannot_be_used_for()

This model is not suitable for:

• ideal dissipation-free MHD

• kinetic or hybrid particle physics

• reduced linear perturbation problems

• pressure or entropy primitive-variable studies

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

Bases: StruphyModel

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

Normalization:

\[\hat u = \hat v_\textnormal{th} \,.\]

Equations:

\[\begin{split}\begin{align} \partial_t \rho + \nabla \cdot (\rho \mathbf u) &= 0\,, \\[2mm] \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}\,, \\[2mm] \partial_t S + \mathbf u \cdot \nabla S &= 0\,, \end{align}\end{split}\]

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

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

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

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

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

\[ \begin{align}\begin{aligned}\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}\,.\end{aligned}\end{align} \]

Propagators (called in sequence):

1. PushEta

2. PushVxB

3. PushVinSPHpressure

4. PushVinViscousPotential

classmethod doc_pde()

PDEs solved by model:

Continuity:

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

Momentum:

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

Entropy:

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

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

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

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

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

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

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

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

classmethod doc_normalization()

The characteristic speed is thermal:

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

classmethod doc_scalar_quantities()

The following scalars are tracked during simulation:

• SPH kinetic energy: en_kin

classmethod doc_discretization()

classmethod doc_long_description()

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

classmethod doc_examples()

Create and initialize the viscous Euler SPH model:

from struphy.models import ViscousEulerSPH

model = ViscousEulerSPH()
model.euler_fluid.var

classmethod doc_use_cases()

This model is appropriate for:

• SPH verification with pressure and viscosity

• meshfree compressible-fluid experiments

• testing particle-based viscosity and pressure pushers

classmethod doc_cannot_be_used_for()

This model is not suitable for:

• FEEC grid-based fluid or MHD formulations

• entropy-resolved thermodynamic evolution

• kinetic plasma physics

• studies that require exact field-based conservation structures

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

Bases: StruphyModel

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

Normalization:

\[\hat u = \hat v_\textnormal{A}\,, \qquad \hat{\mathcal U} = \frac{\hat{\mathbf B}^2}{\hat \rho \mu_0 (\gamma-1)} \,,\qquad \hat s = \hat \rho\ \textrm{ln}\left(\frac{\hat{\mathbf B}^2}{\mu_0 (\gamma -1) \hat{\rho}}\right) \,.\]

Equations:

\[\begin{split}&\partial_t \rho + \nabla \cdot ( \rho \mathbf u ) = 0 \,, \\[4mm] &\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 \,, \\[4mm] &\partial_t s + \nabla \cdot ( s \mathbf u ) = \frac{1}{T}\left((\mu+\mu_a(\mathbf x)) |\nabla \mathbf u|^2 \right) \,,\end{split}\]

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.

Propagators (called in sequence):

1. VariationalDensityEvolve

2. VariationalMomentumAdvection

3. VariationalEntropyEvolve

4. VariationalViscosity

Model info:

classmethod doc_pde()

PDEs solved by model:

Continuity:

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

Momentum:

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

Entropy:

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

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

classmethod doc_normalization()

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

classmethod doc_scalar_quantities()

The following scalars are tracked during simulation:

• Kinetic energy: en_U

• Thermodynamic energy: en_thermo

• Total energy: en_tot

• Total density / entropy: dens_tot, entr_tot

classmethod doc_discretization()

classmethod doc_long_description()

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

classmethod doc_examples()

Create and initialize a viscous-fluid model:

from struphy.models import ViscousFluid

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

classmethod doc_use_cases()

This model is appropriate for:

• nonlinear viscous compressible hydrodynamics

• entropy-based Navier-Stokes benchmarks

• testing variational viscous closures without magnetism

classmethod doc_cannot_be_used_for()

This model is not suitable for:

• magnetic or MHD dynamics

• inviscid strictly conservative benchmarks

• kinetic particle effects

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

Bases: StruphyModel

Vlasov equation in static background magnetic field.

Normalization:

\[\hat v = \hat \Omega_\textnormal{c} \hat x\,.\]

Equations:

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

Propagators (called in sequence):

1. PushVxB

2. PushEta

classmethod doc_pde()

PDEs solved by model:

Vlasov equation:

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

classmethod doc_normalization()

The characteristic speed is the cyclotron scale

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

classmethod doc_scalar_quantities()

The following scalars are tracked during simulation:

• Particle kinetic energy: kinetic_energy

classmethod doc_discretization()

classmethod doc_long_description()

Vlasov is the simplest kinetic test-particle model in the 6D hierarchy. It evolves particles in a static magnetic background without electric or magnetic self-consistency.

classmethod doc_examples()

Create and initialize a Vlasov test-particle model:

from struphy.models import Vlasov

model = Vlasov()
model.kinetic_ions.var

classmethod doc_use_cases()

This model is appropriate for:

• test-particle motion in prescribed magnetic fields

• verification of the Boris-like VxB and PushEta splitting

• reduced kinetic transport studies without field feedback

classmethod doc_cannot_be_used_for()

This model is not suitable for:

• self-consistent electrostatic or electromagnetic coupling

• collisional kinetic dynamics

• guiding-center reduction studies

• fluid or MHD-scale closures

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

Bases: StruphyModel

Vlasov-Maxwell equations for one species.

Normalization:

\[\begin{align} \hat v = c \,, \qquad \hat E = \hat B \hat v\,,\qquad \hat \phi = \hat E \hat x \,. \end{align}\]

Equations:

\[\begin{split}&\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 \,, \\[2mm] -&\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}\,, \\[2mm] &\frac{\partial \mathbf{B}}{\partial t} + \nabla \times \mathbf{E} = 0 \,,\end{split}\]

with the normalization parameters

\[\alpha = \frac{\hat \Omega_\textnormal{p}}{\hat \Omega_\textnormal{c}}\,,\qquad \varepsilon = \frac{1}{\hat \Omega_\textnormal{c} \hat t} \,,\qquad \textnormal{with} \qquad \hat\Omega_\textnormal{p} = \sqrt{\frac{\hat n (Ze)^2}{\epsilon_0 (A m_\textnormal{H})}} \,,\qquad \hat \Omega_{\textnormal{c}} = \frac{(Ze) \hat B}{(A m_\textnormal{H})}\,,\]

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}\begin{align} \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{align}\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 equilibirum magnetic field.

Notes

• The Control variate method for Ampère’s law is optional; in case it is enabled via the parameter file, the following system is solved:

Find \((\mathbf E, \tilde{\mathbf B}, f) \in H(\textnormal{curl}) \times H(\textnormal{div}) \times C^\infty\) such that

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

where \(\tilde{\mathbf B} = \mathbf B - \mathbf B_0\) denotes the magnetic perturbation.

Propagators (called in sequence):

1. Maxwell

2. PushEta

3. PushVxB

4. VlasovAmpereCoupling

Model info:

classmethod doc_pde()

PDEs solved by model:

Vlasov equation:

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

Ampère’s law:

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

Faraday’s law:

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

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

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

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

Moreover, it is assumed that

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

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

classmethod doc_normalization()

The model uses the light speed as reference velocity:

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

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

classmethod doc_scalar_quantities()

The following scalars are tracked during simulation:

• Electric field energy: en_E

• Magnetic field energy: en_B

• Particle kinetic energy: en_f

• Total energy: en_tot

• Optional Gauss-law diagnostic: gauss_error

classmethod doc_discretization()

classmethod doc_long_description()

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

classmethod doc_examples()

Create and initialize a Vlasov-Maxwell model:

from struphy.models import VlasovMaxwellOneSpecies

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

classmethod doc_use_cases()

This model is appropriate for:

• self-consistent electromagnetic kinetic simulations

• one-species PIC benchmarks

• wave-particle interaction studies with evolving magnetic fields

• verification of the full Vlasov-Maxwell splitting

classmethod doc_cannot_be_used_for()

This model is not suitable for:

• multi-species plasma dynamics without extension

• collisional kinetic closures

• reduced electrostatic-only models where magnetic evolution is unnecessary

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