VlasovMethods

AnalyticHamiltonian(β)

The Hamiltonian of the analytic quadratic-potential model, $H(q, p) = ½|p|² + ½β²|q|²$, for a state x = [q; p] of 2N entries, and the ∇H! of the reduced forms of ReducedBasisMethods. The two halves are the two halves of x.

The instance is both:

• the gradient, called as ∇H!(g, t, x, params), which writes ∇H(x) — the first block β²q and the second p — into g and returns it. It is what ReducedBasisMethods.reduced_ode and ReducedBasisMethods.reduced_pode take as their ∇H!, and it is not RBM's hamiltonian_gradient!, which is the reduced gradient of a SymplecticDEIM;
• the Hamiltonian, called as h(t, x, params), which returns H(x).

Neither call reads t or params: the model is autonomous and takes no equation parameters. The gradient accepts any element type, which Newton's JacobianAutodiff needs — the Lobatto IIIA–IIIB(2) run of the PODE form calls the vector field with ForwardDiff.Dual states.

Example

h = AnalyticHamiltonian(6.0)
z₀ = project(b, x₀)
prob = reduced_pode(b, h, (0.0, 1000.0), 1e-4, z₀)

(ent::CollisionEntropy)()

The discrete entropy

\[S_h = \int_\Omega f_s \, \log f_s \, dv ,\]

evaluated on the Gauß-Legendre grid the spline basis already carries, in any number of velocity dimensions.

ElectricField

The electric field of a one-dimensional Vlasov–Poisson problem on particles.

A field deposits particle positions with their weights, solves for the potential, and evaluates the field $E(x) = -\partial_x \phi(x)$ at particle positions. Its interface is two functions and three call forms:

update!(f, x, w, t)     # deposit the positions `x` with the weights `w` at the time `t`
efield!(f, e, x)        # write `E` at each position of `x` into `e`
f(e, x, w, t)           # `update!` followed by `efield!`
f(e, x)                 # `efield!` alone
f(x, w, t)              # both, into an `e` that call allocates

update! and efield! are called once per time step of a march, so a concrete field holds the buffers it needs and neither function allocates. energy(f) and coefficients(f) read the field energy and the potential coefficients of f. A reduced field, such as ReducedElectricField, holds a full-order field and lifts into it.

update! and efield! are not exported, so a caller reaches them through the three call forms above, as VlasovMethods.update! and VlasovMethods.efield!, or with using VlasovMethods: update!, efield!. Neither is the PoissonSolvers.update! that a PoissonField calls on its own potential: that one solves the potential's linear system, and it is a different function that happens to carry the same name.

GeometricIntegrator(model::VlasovPoisson{1, 1, <:GridDistribution}, tspan, tstep,
                    method = Integrators.ImplicitMidpoint())

The GeometricIntegrator that steps the grid 1D1V Vlasov–Poisson model model over tspan with step tstep, in the state vec(model.distribution.values).

method is any GeometricIntegrators method for an ODEProblem. ImplicitMidpoint() keeps $H$; RK4() is the one scripts/landau_damping_1d1v.jl uses, because the implicit midpoint's dense Newton Jacobian is too large on that script's $n_v$. The scheme runs in Float32 and Float64; FFTW has no Float16 plan. The vector field is _grid_1d1v_rhs!.

GridDistribution{DT, XT, VT, FT}

The distribution function of a 1D1V phase space, stored as its values at the nodes of an nx × nv grid. It is a DistributionFunction{DT, 1, 1}, beside ParticleDistribution and SplineDistribution.

The grid

• x is periodic: $x_i = x_a + (i-1) h_x$ for $i = 1, \dots, n_x$, with $h_x = (x_b - x_a) / n_x$. The node $x_b$ is the node $x_a$ again and is not stored.
• v is bounded: $v_j = v_a + (j-1) h_v$ for $j = 1, \dots, n_v$, with $h_v = (v_b - v_a) / (n_v - 1)$. Both ends of the range are nodes.

values is the nx × nv matrix $f_{ij} = f(x_i, v_j)$. Its column-major linear order, x running fastest, is the order the grid stencils such as _apply_∫dv! use, so vec(dist.values) is their state vector.

Evaluation

dist(x, v) interpolates bilinearly between the four surrounding nodes. It is periodic in x and zero outside $[v_a, v_b]$.

Moments

velocity_moments integrates over v by the rectangle rule, $\int f(x_i, v) \, dv \approx \sum_j f_{ij} h_v$. This is the integral of the interpolant plus $\tfrac{h_v}{2} (f_{i1} + f_{i n_v})$, so it is exact up to the values at the two ends of the velocity range and first order in $h_v$ where those do not vanish.

GridDistribution(nx, nv, xdomain, vdomain)

A GridDistribution of Float64 zeros on nx × nv nodes.

GridDistribution(values::AbstractMatrix, xdomain, vdomain)

Wrap the node values values of size nx × nv on the periodic interval xdomain and the bounded interval vdomain. Each domain is anything with a first and a last.

IntegratorParameters(dt, nₜ, nₛ, nₕ, nₚ, nparam)

The parameters of a reduced sweep: the timestep dt, the number of time steps nₜ, the number of saved times nₛ — the initial state included — the number of potential degrees of freedom nₕ, the number of particles nₚ and the number of sampled equation parameters nparam. t holds nₛ times, evenly spaced from 0 to dt·nₜ.

The constructor throws an ArgumentError when nₜ < 1, when nₛ < 2, and when nₛ - 1 does not divide nₜ, so that the saves of the march would write past nₛ or drop the final state. The HDF5 constructor applies the same checks to the stored counts.

Every count is Int, so an Int32 count is a MethodError. The element type T of dt is shared by P₀, Ψ, SS and IC in reduced_integrate_vp, so mixed storage types are a MethodError. The timestep scale params.χ is converted to T by the march.

The HDF5 pair is h5save(h5, IP; path = …) and IntegratorParameters(h5, path). h5save writes into the group at path, and creates it when it does not exist.

PoissonField(potential)

The full-order electric field of the PoissonSolvers Potential potential, the field a particle VlasovPoisson model integrates.

The deposit of update! is _deposit!, the field of efield! is _electric_field, and the solve is the potential's own PoissonSolvers.update!.

potential is on a B-spline basis, as the potential of a particle model is. The constructor holds the buffer of the local_width(basis) basis values that the deposit and the field evaluation fill, and the quadrature of the basis that energy integrates $(\partial_x\phi)^2$ with, so that neither allocates on a call.

PotentialReducedTensor(tensor::GridTensor{DT, 3}, Pi, Pj, Pk)

The lazy size(Pi, 2) × size(Pj, 2) × size(Pk, 2) array that projects the first two indices of the bracket tensor tensor onto the columns of Pi and Pj, and its third index onto the columns of Pk. With R = ReducedTensor(tensor, Pi, Pj),

PotentialReducedTensor[i, j, k] = ∑_K R[i, j, K] * Pk[K, k] .

Pi, Pj and Pk are matrices with one row per grid node. The nodes are counted in the order of the grid's own linear indices, the first direction fastest, as in vec(dist.values) of a GridDistribution. The grid may have any dimension.

The constructor builds ReducedTensor(tensor, Pi, Pj) first, so the checks of that constructor run here as well: Pi and Pj need one row per node, or it throws a DimensionMismatch, their element type has to be DT, or it throws an ArgumentError, and every direction of the grid needs at least 2w + 1 nodes for the tensor's stencil width w, or it throws an ArgumentError. Pk needs one row per node, or the constructor throws a DimensionMismatch, and the element type DT, or it throws an ArgumentError. An index off the array throws a BoundsError.

The tensor's contract is ReducedTensor's and is assumed, not checked: tensor is a bracket's tensor on the grid that bracket was built on, nonzero only where its stencil says so. bracket_tensor builds such a tensor from an Arakawa. A tensor written by hand on another grid than its bracket breaks the contract silently.

getindex allocates nothing, and neither does materialising the array into a preallocated output.

julia> using GeometricBrackets: Arakawa

julia> gt = bracket_tensor(Arakawa(3, 3, 0.5, 0.5));

julia> size(PotentialReducedTensor(gt, ones(9, 2), ones(9, 1), ones(9, 4)))
(2, 1, 4)

ReducedElectricField(field, Ψₓ, Ψₚ, Ψₑ, X)
ReducedElectricField(field, Ψ, X)

The reduced electric field of the full-order field: the lift, solve and project path, with no hyper-reduction. The reduced coordinates X are lifted to the particle positions by Ψₓ, the field is deposited from them and solved, evaluated at the positions lifted by Ψₚ, and projected back by Ψₑ:

update!(f, X, w, t)     # deposit  Ψₓ X  and solve
efield!(f, e, X)        # e = Ψₑ E(Ψₚ X)

A DEIMElectricField of a field basis Ψₑ with DEIM row indices p is the same path, lifted by Ψₓ[p, :] and projected by ΨₓᵀΨₑ Ψₑ[p, :]⁻¹, which evaluates the field at the kₑ interpolation points p instead of at every particle.

For an orthonormal basis Ψ, the three-argument form is the full lift and the transpose projection, ReducedElectricField(field, Ψ, X) = ReducedElectricField(field, Ψ, Ψ, Ψ', X). At full basis size, the reduced field equals the projection of the full-order field onto the basis. See ElectricField for the interface, and DEIMElectricField for the hyper-reduced path.

ReducedIntegratorCache(IP::IntegratorParameters, k)

The work buffers of the reduced march of basis size k: the reduced positions zₓ, velocities zᵥ and accelerations zₐ, each k long, and the nₚ particle weights w. One cache serves a whole march, which is why the march takes it as an argument.

The buffers follow the element type of IP, and their length is not checked against a basis: the march's k is the column count of the basis it is handed.

ScaledField(field::ElectricField, χ::Real)

The full-order field as the rescaled state of a march sees it: efield! divides the field of field by χ², and energy divides its energy by the same χ². That is not the energy of the divided field, which is quadratic in the field and would carry χ⁴. update! and coefficients forward to field unchanged, because the deposit and the solve carry no amplitude.

It is an ElectricField, so a reduced march wraps it as ReducedElectricField(ScaledField(field, χ), Ψ, X). A march takes its scale from params.χ and not from the χ stored here: it marches a ScaledField, or a reduced field over one, at the χ of its parameters (see integrate_vp), so the stored χ serves a call outside a march.

χ must be finite and positive; the constructor throws an ArgumentError otherwise, which refuses 0, a negative and NaN alike. Each division converts χ to the element type of the value it divides, so a Float32 field with a Float64 χ keeps Float32.

Snapshots{T, PS}

The snapshots of a parameter sweep: the positions X, the velocities V, the accelerations A, the potentials Φ, and the diagnostics W, K and M, on the sweep's ParameterSpace.

X, V and A are indexed [component, particle, time, sample] and have the size (nd, nₚ, nₛ, nparam), with nparam = length(paramspace). X[:, :, t, p] is the (nd, nₚ) matrix of the positions saved at the t-th of the nₛ saved times of the sample paramspace(p), X[i, j, t, p] is component i of particle j, and V and A hold the velocities and the accelerations of that state the same way. So vec(X[:, :, t, p]) is the component-fastest order of vec(particles.x) of that state, and a snapshot matrix built from these arrays has one column per saved time and sample.

Φ has the size (nh, nₛ, nparam): the potential is a scalar on the one dimension of the Poisson problem, so it carries no component axis, and nh is the number of degrees of freedom of the sweep's Poisson solver. W, K and M have the size (nₛ, nparam), one value each per saved time and sample. W is the field energy, not a particle weight: the snapshots store no weights.

The component count nd of X, V and A and the degree count nh of Φ are the leading extents of those arrays, so they are read from the arrays rather than stored beside them.

== compares all eight fields, paramspace included.

Snapshots(T, nd, nₚ, nh, nₛ, paramspace)

The zero Snapshots{T} of nd components, nₚ particles, nh potential degrees of freedom, nₛ saved times and length(paramspace) samples, which a sweep fills as it runs. nh = 0 is the sweep of a model with no potential, such as one of VlasovWithExternalFields, and its empty Φ round-trips through HDF5 like any other.

Snapshots(particles::ParticleList, poisson::PoissonSolver, nₛ, paramspace)

The zero Snapshots{T} of a sweep of particles on the Poisson solver poisson, with nₛ saved times and length(paramspace) samples: nd = size(particles.x, 1) components, nₚ = length(particles) particles and nh = length(poisson) potential degrees of freedom. Both element types are T.

particles.v needs the nd rows of particles.x, because X, V and A share one component count and the symplectic structure of the reduced model pairs positions and velocities of one size. A list that holds another count throws an ArgumentError.

SplineDistribution(xdim, vdim, nknots, s_order, domain, length_big_cell, bc = Free())

Build a SplineDistribution on nknots equally spaced breakpoints across domain, of spline order s_order, with boundary condition bc on every velocity axis.

s_order is the order $k = p+1$, the same convention as SimpleSplines.order — so s_order = 4 is the cubic basis both collision-operator manuscripts specify, and s_order = 3 is quadratic.

length_big_cell > 0 appends one oversized cell of that width at each end of domain, outside it, giving a GeneralMesh. That is the device for keeping a particle that strays out of the resolved region inside the support of the basis; with 0 the mesh is uniform on domain exactly.

domain is anything with a first and a last — a two-tuple, a vector or a range.

bc may be a SimpleSplines.BoundaryCondition, a two-tuple for the two ends, or one of the symbols :Dirichlet, :Periodic and :nothing that the drivers in scripts/ use, where :nothing means Free(). Any other symbol is an error. The default is Free() — the unconstrained clamped basis, which is the only choice that reproduces $1$, $v$ and $v^2$ and hence the only one on which the conservative schemes conserve.

The mass matrix is always the exact Galerkin one, $\mathbb{M}_{ij} = \int \varphi_i \varphi_j$; there is no lumped-mass option.

SplineDistribution{DT, XD, VD}

The distribution function represented as a spline in velocity, $f_s(v) = \sum_I \hat{f}_I \, \varphi_I(v)$, together with everything needed to project onto it and to integrate against it.

VD is the number of velocity dimensions. For VD == 1 the basis is a one-dimensional SimpleSplines B-spline basis and the coefficients are a vector; for VD > 1 it is a TensorProductBasis and the coefficients are a VD-dimensional array. Nothing else about the two cases differs — the evaluation, the projection and the mass solve are the same calls.

Fields

• spline — the SimpleSplines.Spline, which shares coefficients. A projection writes that array, and every reference to the spline sees the new function without a rebuild.
• basis — the spline basis.
• quadrature — the assembly table: the basis tabulated at the Gauß-Legendre points, the weights, and the mass operator.
• coefficients — the degrees of freedom $\hat{f}$.

The mass matrix is never formed as a Kronecker product

For VD > 1 the mass matrix is $\mathbb{M}^{(VD)} \otimes \dots \otimes \mathbb{M}^{(1)}$, and a solve is VD one-dimensional solves applied along each axis. That is an identity, not an approximation, and it is what makes the higher-dimensional cases tractable at all: forming the product and factorising it densely would mean a 1681×1681 matrix for the 41-element cubic basis of a two-dimensional velocity space, and 68921² in three dimensions, which does not fit.

The boundary condition is not a free choice

The $L^2$ projection reproduces $\int \pi(v) f \, dv$ for a polynomial $\pi$ exactly when $\pi$ lies in the span of the basis, so the mass, momentum and energy of $f_s$ agree with those of the particles only if $1$, $v$ and $v^2$ are in the span. polynomial_reproduction(basis) reports the largest degree available: p for the clamped basis (Free()), 0 for a periodic one — $v$ is not periodic — and -1 for a Dirichlet-recombined one, where not even the constants survive.

Measured on a cubic basis over $[-8,8]$ by scripts/verify_conservation.jl: Free() reproduces all three to $10^{-14}$, Periodic() gets $v$ wrong by $1.2$ and $v^2$ by $0.12$, and Dirichlet() gets even the constant wrong by $4 \times 10^{-2}$.

UniformField(E, B)

A constant electric and magnetic field, E and B in Cartesian components, stored as SVector{3,T} with T the promoted element type of the two arguments. It is constant in x and in t, so it is the field of the analytic cases: a gyro-orbit at E = 0, and a constant acceleration along B at E ∥ B.

E and B are two vectors of length three, or two SVector{3}s; a length other than three throws the DimensionMismatch of SVector{3}. A constant field of Int coordinates is written UniformField(SVector(0, 0, 0), SVector(0, 0, 0)).

VlasovHamiltonian(model)

The energy of the Vlasov–Poisson model model, as an invariant h(t, q, params),

\[h(t, q, w) = \frac{1}{2}\sum_i w_i v_i^2 + \frac{1}{2}\int (\partial_x\phi)^2 \,\mathrm{d}x ,\]

where $\phi$ solves the Poisson problem for a fresh deposit of the positions in q with the weights params.w. The call leaves that deposit in model.potential. Each acceleration substep deposits its own state before it reads the field, so this does not change a trajectory.

The field term is the Gauss quadrature of $(\partial_x\phi)^2$ on the basis of the potential, tabulated once by the constructor, so a call allocates nothing. Its terms are non-negative, so the sum keeps their relative accuracy also for a small field.

VlasovPoisson(dist::DistributionFunction, potential::Potential)

The Vlasov–Poisson model of one species with the distribution dist, whose field is the solution of potential.

For a ParticleDistribution or a SplineDistribution, the field is the potential's solve for the deposit of the distribution, and potential is on a B-spline basis.

For a 1D1V GridDistribution, the field is $\phi = G \rho$ with $\rho_i = h_v \sum_j f_{ij}$, and the constructor assembles the $n_x \times n_x$ Green's matrix $G$ of the potential's backend once. GeometricIntegrator(model, tspan, tstep, method) steps the model. The constructor throws:

• ArgumentError unless the potential's solver is a PoissonSolverFFT or a PoissonSolverMatrixFree;
• DimensionMismatch unless the potential has one degree of freedom per x-node of dist;
• ArgumentError unless the potential's x-nodes are dist.x, and unless dist has at least three nodes in each direction.

The grid scheme is periodic in v: it is valid only while f vanishes to eps(T)·maximum(f) at both ends of the velocity range.

VlasovWithExternalFields(dist, field; charge = one(T), mass = one(T))

A VlasovModel of one particle species in an external electromagnetic field: the particles of dist, a ParticleDistribution{T,3,3}, pushed by field. charge and mass are the species' charge and mass, in units in which the equation of motion reads ẍ = (q/m)(E + ẋ × B), and both are stored as T. The constructor throws an ArgumentError unless the converted mass is positive. An empty particle list is a vacuum run and not an error: the model, its SplittingMethod and run! work over the empty state.

The field is any value that defines electric_field and magnetic_field for an SVector{3} coordinate; the model puts no bound on the type of its field parameter, so the two functions are the whole contract. UniformField is the field type this package ships, and the ElectromagneticFields extension adds the two methods for that package's FieldFunctions of a Cartesian chart — FieldFunctions(ABCEquilibrium(a, b, c)) is passed to the constructor as it is. In a curvilinear chart the components are not Cartesian, so the extension defines no method and the first kick raises its MethodError.

The method is SplittingMethod(model, tspan, tstep): a Strang splitting of the drift ẋ = v and the Lorentz force v̇ = (q/m)(E + v × B), with the rotational part of the force a Boris push. One step is a half drift, two half Boris pushes with the field evaluated at the step's start, and a half drift. GeometricIntegrators.integrate(method.integrator) and run!(method, h5file) serve it as they serve the Vlasov–Poisson model.

mass_matrix(dist::SplineDistribution)

The assembled mass matrix. For VD > 1 this forms the Kronecker product, which the representation exists to avoid; it is here for tests and diagnostics, not for solving.

mass_operator(dist::SplineDistribution)

The mass operator of the basis — a factorised or circulant one for VD == 1, a KroneckerMass for VD > 1. Solve with mass_solve!(out, mass_operator(dist), rhs).

The conservative Lenard-Bernstein right-hand side,

\[\dot{v}_\alpha = - \nu \left( \frac{f_s'(v_\alpha)}{f_s(v_\alpha)} + A_1 + A_2 v_\alpha \right) ,\]

which is eq:velocity_ode of the Lenard-Bernstein manuscript with the opposite sign.

DEIMElectricField(field, Ψₓ, Ψₑ, p, X)
DEIMElectricField(field, Ψₓ, Ψₑ, Πₑ, X)

The hyper-reduced electric field of the full-order field: a ReducedElectricField lifted by Ψₓ[p, :] and projected by ΨₓᵀΨₑ Ψₑ[p, :]⁻¹, where p are the DEIM row indices of the field basis Ψₑ, from ReducedBasisMethods.get_DEIM_indices. The second form takes the DEIM interpolation matrix Πₑ of get_DEIM_interpolation_matrix, reads the row of the nonzero entry of each column, and calls the first; it does not check that Πₑ is a selection matrix. The buffers of the field follow the element type of Ψₓ and X, and X is the reduced coordinate they are sized from. With Πₑ = I[:, p], the lift of the field evaluation is

\[E_{\text{DEIM}} = \Psi_x^T \Psi_e (\Pi_e^T \Psi_e)^{-1} E(\Pi_e^T \Psi_x X),\]

so the field is evaluated at the $k_e$ interpolation points rather than at every particle. For orthonormal position and field bases, the reduced acceleration agrees with the non-hyper-reduced one to the DEIM error bound of the paper ReducedBasisMethods.get_DEIM_interpolation_matrix cites.

FullyReducedTensor(tensor::GridTensor{DT, 3}, Pi, Pj, Px)

The lazy size(Pi, 2) × size(Pj, 2) × size(Px, 2) array that projects the third index of the bracket tensor tensor onto the columns of Px, read at the x index of the node. With R = ReducedTensor(tensor, Pi, Pj) and (n₁, n₂) the grid indices of the node K,

FullyReducedTensor[i, j, k] = ∑_K R[i, j, K] * Px[n₁, k] .

Px has one row per x node, in the order of dist.x, and not one per grid node. The result is the PotentialReducedTensor whose projection has the row Px[n₁, :] for each node.

The grid needs two dimensions, or the function throws an ArgumentError. Px needs nx rows, or it throws a DimensionMismatch, and the element type DT, or it throws an ArgumentError. PotentialReducedTensor gives the other checks, the node order, the tensor's contract and the allocations.

julia> size(FullyReducedTensor(gt, ones(9, 2), ones(9, 1), ones(3, 4)))   # gt as above
(2, 1, 4)

The plain Lenard-Bernstein right-hand side,

\[\dot{v}_\alpha = - \nu \left( \frac{f_s'(v_\alpha)}{f_s(v_\alpha)} + v_\alpha \right) ,\]

i.e. the conservative operator of eq:cons_LB with the coefficients fixed at $A_1 = 0$, $A_2 = 1$. It conserves mass but not momentum or energy — that is the whole point of the conservative variant.

MaxwellianDistribution(v)

The normalised Maxwellian of unit temperature and zero mean in d = length(v) velocity dimensions,

\[f_M(v) = (2\pi)^{-d/2} \, \exp\!\left(-\tfrac{1}{2} |v|^2\right) , \qquad \int_{\mathbb{R}^d} f_M \, dv = 1 .\]

VelocityReducedMatrix(tensor::GridTensor{DT, 3}, Pi, Pj, v²)

The lazy size(Pi, 2) × size(Pj, 2) matrix that contracts the third index of the bracket tensor tensor with v²/2. With R = ReducedTensor(tensor, Pi, Pj) and (n₁, n₂) the grid indices of the node K,

VelocityReducedMatrix[i, j] = ∑_K R[i, j, K] * v²[n₂] / 2 .

v² holds the squared velocity nodes, dist.v .^ 2 of a GridDistribution: one entry per velocity node, in the order of dist.v. The result is view(prt, :, :, 1) of the PotentialReducedTensor prt whose one-column projection carries v²[n₂] / 2 in the row of each node.

The grid needs two dimensions, or the function throws an ArgumentError. v² needs nv entries, or it throws a DimensionMismatch, and the element type DT, or it throws an ArgumentError. PotentialReducedTensor gives the other checks, the node order, the tensor's contract and the allocations.

julia> size(VelocityReducedMatrix(gt, ones(9, 2), ones(9, 1), [1.0, 4.0, 9.0]))   # gt as above
(2, 1)

_get_MC̃_cubic(V, ∫dv, ∫vdv, ∫v²dv, v, ci, li, h₁, h₂)

Assemble the cubic, four-index collision tensor MC̃[o, a, b, c] on the n₁ × n₂ phase-space grid, from the evaluation matrix V, the three moment matrices ∫dv, ∫vdv and ∫v²dv of the rectangle rule, the velocity grid v, the CartesianIndices ci and the LinearIndices li of the n₁ × n₂ grid, and the grid widths h₁ and h₂.

v holds one entry per velocity node, so the call throws a DimensionMismatch under length(v) ≠ size(ci, 2). The tensor has the element type eltype(V).

Contracted with the grid values, Q₃[o] = Σ MC̃[o, a, b, c] f[a] f[b] f[c] is the Galerkin form of ρ²·C[f]. It conserves the discrete moments ∫dv·f, ∫vdv·f and ∫v²dv·f — the mass, the momentum and the energy — while f vanishes at the ends of the v-grid to eps(T). The v-stencil wraps the bounded grid with mod1, so a non-negligible end value breaks the conservation.

_get_MC̃_quadratic(V, ∫dv, ∫vdv, ∫v²dv, v, ci, li, h₁, h₂)

Assemble the quadratic, three-index collision tensor MC̃[o, a, c] on the n₁ × n₂ phase-space grid, from the same arguments as _get_MC̃_cubic. v holds one entry per velocity node, so the call throws a DimensionMismatch under length(v) ≠ size(ci, 2), and the tensor has the element type eltype(V).

Contracted with the grid values, Q₂[o] = Σ MC̃[o, a, c] f[a] f[c] is the Galerkin form of ρ²·C[f] without the ρ̂u term. It conserves the mass ∫dv·f and the energy ∫v²dv·f, but not the momentum ∫vdv·f: it pulls the mean velocity toward v = 0, by a relative 1.3e-2 over ten RK4 steps of Δt = 10⁻³ on nx = 2, v = range(-10, 10; length = 25), from the shifted Maxwellian of the invariants testset in test/gridbased/collisions.jl. The moments are conserved only while f vanishes at the ends of the v-grid to eps(T).

_green_matrix(dist::GridDistribution, potential::Potential)

The $n_x \times n_x$ Green's matrix $G$ of potential's backend on the grid of dist, which gives the field of a density in one product:

\[\rho_i = h_v \sum_j f_{ij} , \qquad \phi = G \rho .\]

It solves the backend for the $n_x$ unit vectors and symmetrises the result. The symmetry makes $\partial H / \partial f = h_x h_v h$ exact, so the discrete Hamiltonian is quadratic in $f$ and an implicit midpoint step conserves it. A spline basis is refused: its density at the nodes is not a Galerkin load vector. PoissonSolverMatrixFree bounds its residual, not its error, so in Float32 on a fine grid G is less accurate than $eps(Float32)$. The errors are listed under VlasovPoisson.

_grid_1d1v_rhs!(ḟ, t, f, params)

Arakawa's discretisation of the 1D1V Vlasov–Poisson right-hand side, as the vector field of a GeometricEquations.ODEProblem:

\[\dot{f} = - [f, h] , \qquad h_{ij} = \frac{v_j^2}{2} + \phi_i , \qquad \phi = G \rho , \qquad \rho_i = h_v \sum_j f_{ij} ,\]

for a charge of $+1$, with G from _green_matrix. params holds the bracket arakawa, the velocity nodes v and G; the state f is vec(dist.values). The bracket is doubly periodic, so a distribution that reaches the ends of the velocity range wraps mass from $v_b$ to $v_a$ (see VlasovPoisson).

The largest absolute value of D along its first axis, with the trailing axes kept.

The root sum of squares of D along its first axis, with the trailing axes kept. Unlike norm, it does not rescale, so a difference above about sqrt(floatmax(T)) overflows to Inf.

The root mean square of D along its first axis, with the trailing axes kept: the root sum of squares divided by the square root of the length of that axis.

The (n₁·n₂) × (n₃·n₄) reshape of Y; its column (p - 1) * n₃ + t is vec(Y[:, :, t, p]).

best_and_worst(γh, γr)

The indices of the equation parameter whose reduced growth rate is nearest to, and furthest from, its full-order value: the argmin and the argmax of abs.(γr .- γh). Each is the first index of a tie.

bracket_tensor(a::Arakawa{DT})

The rank-3 GridTensor of the bracket a on a's own nx × nv grid, with stencil width 1: its entry at the nodes (I, J, K) is a(I, J, K). The grid is read from a, so the tensor cannot be on another grid than its bracket. This is the tensor that PotentialReducedTensor, VelocityReducedMatrix and FullyReducedTensor take.

bracket_tensor is a bridge until VlasovMethods requires GeometricBrackets 0.3, whose four-argument PoissonTensor(DT, nx, nv, f) builds the same tensor and checks the grid.

julia> using GeometricBrackets: Arakawa

julia> size(bracket_tensor(Arakawa(3, 4, 0.5, 0.5)))
(12, 12, 12)

check_conservation_basis(dist::SplineDistribution; moments = 2)

Throw unless the basis reproduces every polynomial of degree ≤ moments, i.e. unless the moments $\int v^k f_s \, dv$ for k ≤ moments agree with the particles'.

Not called automatically: a diagnostic run may legitimately want a Dirichlet or periodic basis, and refusing to build one would be wrong. Call it where a result depends on the two representations agreeing — before reporting an entropy, or before comparing a spline moment against a particle one.

compute_J!(J, sdist, ::Landau)

The vector

\[\mathbb{L}_k = \sum_i \mathbb{M}^{-1}_{ik} \int_\Omega \varphi_i(v) \left[ 1 + \log \Big( \sum_j f_j \varphi_j(v) \Big) \right] dv\]

of the Landau manuscript's eq:defn_Lk — the $L^2$ projection of $1 + \log f_s$ onto the spline basis.

Written as a projection rather than as one quadrature per coefficient. The integrand is sampled once on the tensor-product Gauß-Legendre grid the basis already carries and contracted against the one-dimensional tabulations, so the cost is the grid size rather than the grid size times the number of coefficients. The earlier version called a gauss_quad_2d that looped over every cell pair for each of the $M$ coefficients, and took the quadrature domain from a hard-coded knots[3:(end-1)] slice — asymmetric, independent of the spline order and of the boundary condition, and a different subrange from the one the mass matrix was assembled on, so the two halves of eq:defn_Lk were not computed over the same domain.

compute_J!(J, sdist, ::MetriplecticLenardBernstein)

The vector $\mathbb{L}_k = \sum_i \mathbb{M}^{-1}_{ik} \int \varphi_i \, (1 + \log f_s) \, dv$ of the Landau manuscript's eq:defn_Lk, which is exactly the $L^2$ projection of $1 + \log f_s$ onto the basis.

This discretisation projects $1 + \log f_s$ onto the spline space and differentiates the projection, where eq:velocity_ode of the Lenard-Bernstein manuscript uses the pointwise ratio $f_s'(v_\alpha)/f_s(v_\alpha)$. The two agree only up to the projection error of the logarithm. That is deliberate and is what makes this a genuine discrete-gradient system: J is $\partial S_h / \partial f_i$ for $S_h = \int f_s \log f_s \, dv$.

compute_K!(K1, K2, v_array, sdist, landau)

The matrix $\mathbb{K}_{k,\alpha} = w_\alpha \, \nabla \varphi_k(v_\alpha)$ of the Landau manuscript's appendix, split into its two velocity components.

Every index is bounded per component before being flattened. Testing only the flat index does not bound the components: with i = 0, j = 3, M = 10 the flat index (j-1)M + i = 20 passes a 1 ≤ k ≤ M² test and decodes to (10, 2), so an out-of-support contribution is aliased onto an unrelated basis function instead of being discarded. Wrapping goes through basis_index, which is periodic-aware, so the gradient is now evaluated with the same wrap as the spline it differentiates — the two disagreed before.

compute_L!(L, sdist, landau; chunk = 256)

The symmetric matrix

\[\mathbb{L}_{ij} = \frac{1}{2} \int_\Omega \! \int_\Omega \big( \nabla \varphi_i(v) - \nabla \varphi_i(v') \big) \cdot f_s(v) \, U(v,v') \, f_s(v') \cdot \big( \nabla \varphi_j(v) - \nabla \varphi_j(v') \big) \, dv \, dv'\]

of the Landau manuscript's eq:discrete-landau-matrix.

Method

Expanding the two differences gives four terms. With $U$ symmetric in its arguments and in its indices, relabelling $v \leftrightarrow v'$ shows the first and fourth are equal and the second and third are equal, so the factor of one half cancels and

\[\mathbb{L} = \sum_{c,d} D^c \operatorname{diag}(s \odot A^{cd}) (D^d)^{\mathsf T} - \sum_{c,d} \big( D^c \operatorname{diag}(s) \big) \, U^{cd} \, \big( D^d \operatorname{diag}(s) \big)^{\mathsf T} , \qquad A^{cd}_a = \sum_b s_b \, U^{cd}_{ab} ,\]

with $s_a = w_a f_s(v_a)$ and $D^c$ the sparse gradient tables above. This is an identity, not an approximation: the same quadrature sum, rearranged.

Both terms are sparse-times-dense products. The first is a single pass over the grid; the second is the genuinely coupled one, and its $Q \times Q$ kernel is built and consumed in row blocks of chunk so that nothing of that size is ever held.

$U^{cd}$ is symmetric in $(c,d)$, so each block is swept once and its three independent components serve all four contractions. The whole assembly therefore costs $O(Q^2)$ kernel evaluations, two sweeps of the grid pairs in total: one for $A^{cd}$ and one for the blocks.

compute_coefficients(distribution, particle_dist, vp)

The coefficients $A_1$, $A_2$ of the conservative Lenard-Bernstein drift, from the requirement that the operator conserve momentum and energy:

\[\begin{pmatrix} n_h & n_h u_h \\ n_h u_h & n_h \varepsilon_h \end{pmatrix} \begin{pmatrix} A_1 \\ A_2 \end{pmatrix} = - \sum_\alpha w_\alpha \begin{pmatrix} 1 \\ v_\alpha \end{pmatrix} \frac{f_s'(v_\alpha)}{f_s(v_\alpha)} , \qquad \begin{aligned} n_h &= \textstyle\sum_\alpha w_\alpha , \\ n_h u_h &= \textstyle\sum_\alpha w_\alpha v_\alpha , \\ n_h \varepsilon_h &= \textstyle\sum_\alpha w_\alpha v_\alpha^2 , \end{aligned}\]

solved by Cramer's rule. This is eq:coefficient_lin_system of the Lenard-Bernstein manuscript.

compute_coefficients_rclb(distribution, particle_dist, vp)

The coefficients of the rescaled parametrisation, in which the coefficient of $v$ is normalised to one:

\[\dot{v}_\alpha = - \nu \left( A_1 \frac{f_s'(v_\alpha)}{f_s(v_\alpha)} + A_2 + v_\alpha \right) .\]

This is the conservative operator with time rescaled by $\varepsilon_h - u_h^2$, so that $A_1^{\mathrm{resc}} = 1 / A_2^{\mathrm{cons}} = \varepsilon_h - u_h^2$ and $A_2^{\mathrm{resc}} = A_1^{\mathrm{cons}} / A_2^{\mathrm{cons}} = -u_h$. It is not the fixed-coefficient cumulant-scaling experiment of the manuscript appendix, which sets the coefficients by hand rather than solving for them.

The particle weights appear in every sum, for the reason given at compute_coefficients.

compute_dS!(dS, J, v, sdist, ::MetriplecticLenardBernstein, pdist)

$\partial S_h / \partial v_\alpha = w_\alpha \sum_k \mathbb{L}_k \, \varphi_k'(v_\alpha)$, the Landau manuscript's eq:entropy_derivative up to its overall sign.

The weight $w_\alpha$ is read from pdist, which is why the particle distribution is passed in: it equals $1/N$ for every initialiser in examples/ but not for the importance-sampling one, so 1/length(dS) is not a substitute.

compute_entropy(f, mlb::MetriplecticLenardBernstein)

$S_h = \int f_s \log f_s \, dv$, the quantity the manuscript's entropy figures report.

Integrated on the Gauß-Legendre grid of the basis rather than adaptively, so that it is the same quadrature the discretisation itself uses. Returns the value alone; the earlier version returned a quadgk error estimate alongside it and wrote 0.5*log(f^2), so it reported $\int f \log|f|$ rather than the entropy.

Returns particles drawn from g for the bump-on-tail case.

input: nb. pf particles Nₚ, x-marginal function gₓ, parameters μ output: Particles struct

Returns particles drawn from g for the bump-on-tail case.

input: nb. pf particles Nₚ, x-marginal function gₓ, parameters μ output: Particles struct

Returns particles drawn from g for the bump-on-tail case.

input: nb. pf particles Nₚ, x-marginal function gₓ, parameters μ output: Particles struct

efield!(f::PoissonField, e, x)

Write the field $E(x_i) = -\partial_x\phi(x_i)$ at each position of x into e, and return e.

electric_field(field, t, x)
magnetic_field(field, t, x)

The electric and magnetic field of field at the time t and the position x, each an SVector{3,T} in Cartesian components, with T the element type of x.

The two functions are the interface between a VlasovWithExternalFields model and the field its particles move in. The model calls them once per particle and per kick, with x the particle's position and t the time the kick starts at, and a field type defines one method of each. There is no abstract field type and no fallback method, so a field that defines neither raises its MethodError at the first kick inside integrate.

UniformField is the one field type in this package. The ElectromagneticFields extension adds methods for that package's FieldFunctions of a Cartesian chart; see VlasovWithExternalFields.

energy(f::PoissonField)

The field energy $\frac{1}{2}\int (\partial_x\phi)^2 \,\mathrm{d}x$, by the quadrature of the potential's basis that the constructor tabulates. It is the diagnostic a snapshot march saves at every state, so it is type-stable and allocates nothing.

energy_drift(H)

The relative drift of an energy history H: (H .- H[1]) ./ H[1], the energy at each saved time against the initial one.

The energies of a conservative model are constant, so the drift is its discretisation error and H[1] is its reference. NaN and Inf in H propagate into the result.

Throws

• ArgumentError when H is empty, and when H[1] is zero, because a relative drift needs a non-zero reference.

field_basis(algorithm::Union{PODSVD, PODEVD}, s::Snapshots; parameters = NamedTuple(),
            kwargs...)

The POD ReducedBasis of the snapshot matrix of the accelerations s.A, whose nd * nₚ rows match the position block of particle_basis. algorithm is the tag PODSVD or PODEVD of ReducedBasisMethods, which VlasovMethods does not export. A PSD or complex tag would split the rows of A into positions and velocities, so it throws a MethodError.

parameters, s.paramspace and kwargs, k or tolerance, go to ReducedBasis unchanged.

field_energy_error(Wh, Wr)

The relative error of a reduced field energy against a full-order one, abs.(Wh - Wr) ./ Wh, per saved time and parameter.

A field energy is ½ϕᵀKϕ ≥ 0 by construction, so the division is by a non-negative number and no guard is needed: a Wh of 0 gives Inf, or NaN where Wr is also 0, and both propagate.

Wh and Wr must have the same size; a difference raises a DimensionMismatch naming both.

full_leapfrog_step!(x, v, a, w, efield, Δt, t)

One step Δt of the full-order leapfrog: half a drift, the deposit of the half-drifted positions and the kick of the full-order field at the time t, and the second half drift. x and v are the particle positions and velocities, w their weights, and a the acceleration buffer the field writes. The symmetric composition is second order. It is inferred and allocates nothing.

gradient_tabulations(sdist)

The sparse tables $D^c_{k,a} = \partial_c \varphi_k(v_a)$ of the two velocity components of the basis gradient at every point of the tensor-product quadrature grid, as an M × Q pair.

Built as Kronecker products of the one-dimensional tabulations the quadrature already holds: $\partial_1 (\varphi_i \otimes \varphi_j) = \varphi_i' \otimes \varphi_j$, so D¹ = kron(Φ₂, Φ₁') and D² = kron(Φ₂', Φ₁) — with the row and column orderings that LinearIndices and the grid enumeration already use. Nothing is evaluated that the quadrature had not already tabulated.

growth_rate(t, W; window = 2, maxima = 2:5)

The growth rate $γ$ of the field energy history W at the times t: half the least-squares slope of $\log W$ over the selected local maxima, because $W ∝ |E|² ∝ e^{2γt}$. The times are physical time, so this is the physical growth rate.

A maximum is a point no lower than the window points on each side, whose first differing sample on its right is lower. A flat maximum counts once, at its first index, and a rising shoulder, a plateau that ends in a higher sample, is no maximum. maxima selects which of them the fit uses, by their order along t, which drops the first few maxima; the fit runs in eltype(W), so it takes no dependency on Statistics.

A matrix W gives one γ per column, as a Vector of the element type of W.

Throws

• ArgumentError when window < 1; when maxima names an index twice, naming it, because one maximum may not enter the fit twice; when it selects fewer than two maxima, or reaches past the maxima found, each naming the number found; when it selects an index below 1; and when W is not positive at a selected maximum, where $\log W$ is not defined.
• DimensionMismatch when t and W differ in length, naming both.

integrate_vp(P₀, params, efield, SS, IP, p; save = true)

The full-order leapfrog march of an ElectricField, the same march as reduced_integrate_vp without a basis: the state of P₀ advanced IP.nₜ steps of IP.dt by full_leapfrog_step!, with one kick of the full field of efield per step. A reduced model marches the ScaledField of the same field, so the two marches of one comparison take the same field and the same χ, which is params.χ. With nsave = nₜ/(nₛ - 1), the state of every nsave-th step it, the initial state included, is written into the save div(it, nsave) + 1 of the sample p of the snapshots SS by save_full_solution.

The march holds its own buffers, so efield is the only state it shares with the caller: P₀ is copied and never written, and the march needs no cache argument.

params is the caller's, and it must hold χ. The χ scale is a timestep scale and nothing else: the step is Δt = IP.dt·T(params.χ), converted to the march's element type before the product, and the times handed to the field are unscaled: the kick of step it is at its midpoint (it - ½)·IP.dt, and the save of step it at it·IP.dt. So a time-independent field marched with (dt, χ = 2) equals the same field marched with (2dt, χ = 1) state for state.

params.χ is the one source of the scale. When efield is a ScaledField, the march rebuilds it over the same field with params.χ, so the χ it was built with is not used, and the field and its energy are divided by the χ² of the step. Any other field is marched as it is.

It returns nothing; the states are in SS, and the integration parameters hold the times.

Throws

• ArgumentError when the snapshots carry another number of components than one, and when params.χ is not finite and > 0, which refuses 0, a negative, Inf and NaN alike, as ScaledField does. The counts of IP are checked by IntegratorParameters itself.
• DimensionMismatch when the snapshots do not hold IP.nₛ saved times, and when the particle list holds another number of particles than the snapshots do.
• The NamedTuple's own error when params holds no χ. There is no default.

Cost

The loop's two calls — full_leapfrog_step! and save_full_solution — are inferred and allocate nothing. The buffers of the state are allocated once, before the loop.

kernel(v_α, v_β, ::Landau)

The Landau collision kernel

\[U(u) = \frac{1}{|u|} \left( \mathbb{I} - \frac{u \otimes u}{|u|^2} \right) , \qquad u = v_\alpha - v_\beta ,\]

the orthogonal projector onto the complement of $u$, scaled by $1/|u|$.

This is the standard Landau tensor and it is not what the manuscript's eq:landau_kernel as typeset says: that carries a prefactor $1/|u|^3$ in front of the bracketed projector, i.e. an extra $|u|^{-2}$. The intended reading is $|u|^{-3} ( |u|^2 \mathbb{I} - u \otimes u )$, which is this. The properties the conservation proofs rest on hold for the expression implemented here: $U$ is symmetric, $U(v_\alpha, v_\beta) = U(v_\beta, v_\alpha)$, and $U u = 0$ exactly.

leapfrog_step!(IC::ReducedIntegratorCache, efield, Δt, t)

One step Δt of the reduced leapfrog: half a drift, the kick of the reduced field at the half-drifted position and the time t, and the second half drift. The symmetric composition is second order. It is inferred and allocates nothing.

lift_solution(b::ReducedBasis, sol) -> (q, p)

The full-order positions and velocities of the GeometricIntegrators solution sol of the reduced PODE form in the basis b: the two blocks of the lifted trajectory, each of N rows and one column per stored time.

The solution's own halves sol.q and sol.p are the reduced coordinates, and ReducedBasisMethods.lift(b, [q; p]) is [Φq; Φp] for a cotangent lift and [Φa - Ψb; Ψa + Φb] for a complex SVD. The two blocks are views of that one lifted array. A POD basis has no PODE solution and no q and p halves to split, so it is a MethodError.

magnetic_field(field, t, x)

The magnetic field of field at the time t and the position x. See electric_field for the interface the two functions share.

particle_basis(algorithm, s::Snapshots; parameters = NamedTuple(), kwargs...)

The ReducedBasis of the stacked [Q; P] snapshot matrix of s: the nd * nₚ rows of the positions X above the nd * nₚ rows of the velocities V, one column per saved time and sample. algorithm is one of the tags PODSVD, PODEVD, CotangentLiftSVD, CotangentLiftEVD and ComplexSVD of ReducedBasisMethods, which VlasovMethods does not export.

parameters, s.paramspace and kwargs, k or tolerance, go to ReducedBasis unchanged.

plot_distribution(dist) -> Figure
plot_distribution(dist; xdomain, vdomain, nx = 100, nv = 100) -> Figure

A phase-space portrait of dist on one axis, x horizontally and v vertically.

Defined in a package extension: load CairoMakie to make it available.

using VlasovMethods, CairoMakie
fig = plot_distribution(dist)

The one-argument method takes a GridDistribution and draws dist.values as a single heatmap on the grid's own nodes, untransposed: the columns of the matrix are the v nodes, so the figure and the matrix have the same orientation.

The method with the two keywords takes a ParticleDistribution of one position and one velocity dimension and draws the weighted particle histogram on nx × nv cells of xdomain and vdomain. Cell (i, j) holds

\[F_{ij} = \frac{1}{h_x h_v} \sum_k w_k ,\]

summed over the particles that fall in it, with $h_x = (x_b - x_a) / n_x$ and $h_v = (v_b - v_a) / n_v$, so the histogram integrates to the total weight of the particles it keeps. A cell is half-open, [x_a + (i-1) h_x, x_a + i h_x) and likewise in v.

xdomain is periodic and x wraps into it. vdomain is bounded, and a particle outside it — v < v_a or v >= v_b — is dropped from the histogram rather than clamped into an end cell. A particle with a non-finite x, which the wrap cannot place, is dropped as well. A distribution with no particles throws an ArgumentError.

There are no other keywords, and the two domains are required. Use the returned Figure to set a title, a size or an axis label.

plot_particles(dist; xdomain, vdomain) -> Figure
plot_particles(SS, ts, p; xdomain, vdomain) -> Figure

The particles of a ParticleDistribution of one position and one velocity dimension, as a single scatter of (x, v).

Defined in a package extension: load CairoMakie to make it available.

using VlasovMethods, CairoMakie
fig = plot_particles(dist; xdomain = (0.0, 2π), vdomain = (-4.0, 4.0))

xdomain is periodic and x wraps into it. Every particle's velocity is plotted as it is, and the axis limits are set to xdomain horizontally and vdomain vertically, so a particle at v outside vdomain is left at that v and lies outside the axes. A distribution with no particles throws an ArgumentError.

The second method takes a Snapshots, and draws its saved state: the positions and the velocities of the one component, SS.X[1, :, ts, p] against SS.V[1, :, ts, p], as the scatter above, with x wrapped into xdomain. It plots the state that a march saved, so a frame of ts = saved_indices(nₛ)[1] shows the final state of the sweep. ts and p must index the saved times and the parameters, else an ArgumentError names the range it left.

There are no other keywords, and the two domains are required. Use the returned Figure to set a title, a size or an axis label.

plot_series(t, series; labels = nothing, xlabel = "t", ylabel = "") -> Figure

One line per element of series against the abscissae t, on one axis.

Defined in a package extension: load CairoMakie to make it available.

using VlasovMethods, CairoMakie
fig = plot_series(ts, (x2, xmax); labels = ["x2", "xmax"], ylabel = "error")
fig = plot_series(ts, (W_hfm, W_rom); labels = ["HFM", "ROM"], ylabel = "field energy")
fig = plot_series(μ, (γ_hfm, γ_rom); labels = ["HFM", "ROM"], xlabel = "parameter", ylabel = "γ")

series is any collection of vectors of length(t) values: the entries of the NamedTuple of trajectory_errors for one parameter, or a column of each entry for (nₚ, nₛ, nparam) input; the potential_error of each parameter; the field energies W of a full-order and a reduced march, the energy H of a conservative one, or its drift ΔH/H, against the saved times; or the growth_rate values of a sweep against the equation parameter. t holds the abscissa of each value, so one figure holds the quantities that share it.

labels names the lines. It is nothing, which draws no legend, or a collection with one label per element of series, in the same order, which is drawn as a legend of the axis. xlabel and ylabel label the axes. The axis is linear, so set a log scale on the axis of the returned Figure to read a decaying error, or values that span orders of magnitude.

Throws

• ArgumentError when series holds no line, or a line holds no value.
• DimensionMismatch when a series is not as long as t, and when labels does not hold one label per series.

potential_error(Φh, Φr)

The error of a reduced potential against a full-order potential: the root sum of squares of Φh - Φr over the coefficient axis, per saved time and parameter. It is the absolute $L^2$-type error of the potential, with no divisor.

Φh and Φr hold the coefficients of one potential per (ts, p), as (nh, nₛ, nparam) arrays, and the result is (nₛ, nparam) of the element type of the input. NaN and Inf propagate.

The arrays must have the same size; a difference raises a DimensionMismatch naming both.

project_Maxwellian(sdist::SplineDistribution)

Project the normalised Maxwellian onto the basis of sdist.

The projection goes through the Gauß-Legendre quadrature the basis already carries, in any number of velocity dimensions. project_function calls f with a scalar when VD == 1 and with a VD-tuple otherwise, and SVector accepts both, so one method covers every dimension.

project_function(f, sdist::SplineDistribution)

The $L^2$ projection of the function f onto the basis of sdist, written into its coefficients.

f is called with a scalar for a one-dimensional velocity space and with a VD-tuple otherwise. This is eq:numerical-example-1-coefficients of the Lenard-Bernstein manuscript, $\hat{f} = \mathbb{M}^{-1} \int f \varphi_i \, dv$, and it works in any number of velocity dimensions rather than only in two.

projection!(potential::PoissonSolvers.Potential, distribution::ParticleDistribution)

Deposit the particle weights of distribution onto the basis of potential, writing the charge density into the potential's right-hand side buffer.

Each particle contributes to the local_width(basis) basis functions that do not vanish at its position. evaluate_all! reduces a periodic argument into the domain and writes those values into a buffer, returning the index of the first; basis_index is what wraps it, and is the identity where the basis is not periodic. The particle state keeps the unwrapped position, so the reduction belongs to the evaluation.

A recombined basis has no single width: a function near an end spans the union of two parent supports, so local_width is the largest block and a cell with fewer nonzero functions leaves the tail of the buffer zero. Those padding entries carry indices past nbasis(basis), which is why the zeros are skipped rather than added.

projection(velocities, dist::ParticleDistribution, final_dist::SplineDistribution)

Project the particles onto the basis of final_dist, writing the coefficients into it and returning its Spline.

velocities is a vector for a one-dimensional velocity space and a VD × N matrix otherwise; it is passed separately from dist because the implicit solvers evaluate the right-hand side at a stage value rather than at the particles' stored positions.

reduced_integrate_vp(P₀, Ψ, params, efield, SS, IP, IC, p; save = true)

The reduced leapfrog march of the analytic quadratic model: the reduced state of P₀ — the two coordinates Ψ'·vec(P₀.x) and Ψ'·vec(P₀.v) — advanced IP.nₜ steps of IP.dt by leapfrog_step!, with one kick of the reduced acceleration of efield per step. With nsave = nₜ/(nₛ - 1), the state of every nsave-th step it, the initial state included, is written into the save div(it, nsave) + 1 of the sample p of the snapshots SS by save_solution. Ψ is the position basis, and the same basis lifts the velocities: the cotangent lift's Φ.

params is the caller's, and it must hold χ. The χ scale is a timestep scale and nothing else: the step is Δt = IP.dt·T(params.χ), converted to the march's element type before the product, and the times handed to the field are unscaled: the kick of step it is at its midpoint (it - ½)·IP.dt, and the save of step it at it·IP.dt. So a time-independent field marched with (dt, χ = 2) equals the same field marched with (2dt, χ = 1) state for state.

params.χ is the one source of the scale. When efield is a ScaledField or a ReducedElectricField over one, the march rebuilds it over the same field with params.χ, so the χ it was built with is not used. Any other field is marched as it is.

It returns nothing; the states are in SS, and the integration parameters hold the times.

Throws

• ArgumentError when the snapshots carry another number of components than one; when the basis has a zero column count or more columns than rows; and when params.χ is not > 0, which refuses 0, a negative and NaN alike. The counts of IP are checked by IntegratorParameters itself.
• DimensionMismatch when the snapshots do not hold IP.nₛ saved times.
• The NamedTuple's own error when params holds no χ. There is no default.

Cost

The loop's two calls — leapfrog_step! and save_solution — are inferred and allocate nothing.

$(\varepsilon_h - u_h^2) \, \mathbb{L} \, \partial S_h / \partial v$, i.e. the bracket of rhs_downstairs_factor! with its denominator cleared.

Useful where the common factor is supplied elsewhere; it is not $\dot{v}$, which is what rhs_downstairs_factor! returns.

The metriplectic Lenard-Bernstein right-hand side, $\dot{v} = \mathbb{L} \, \partial S_h / \partial v$ with the bracket

\[\mathbb{L}_{\alpha\beta} = - \frac{n_h}{w_\alpha} \, \delta_{\alpha\beta} + \frac{(\varepsilon_h - u_h v_\alpha) + (v_\alpha - u_h) v_\beta}{\varepsilon_h - u_h^2} .\]

$\mathbb{L}$ is symmetric and satisfies $\sum_\alpha w_\alpha \mathbb{L}_{\alpha\beta} = 0$ and $\sum_\alpha w_\alpha v_\alpha \mathbb{L}_{\alpha\beta} = 0$, so momentum and energy are exact Casimirs of the bracket whatever dS is — a stronger structure than the manuscript's, and the reason these runs conserve.

The weight is per particle, read from pdist for each $\alpha$: the two degeneracies above are what fail if a single weight is used for every particle, so with non-uniform weights the conservation the bracket is built for would be lost.

save_full_solution(SS, x, v, w, efield, p, t, ts)

Write the state of a full-order march at the time t into the save ts of the sample p of the snapshots SS: the positions x and the velocities v, the potential's coefficients, and the three diagnostics. The accelerations SS.A are not written.

The kinetic energy K = ½∑ⱼwⱼvⱼ² and the scalar momentum M = ∑ⱼwⱼvⱼ are those of a single component, nd == 1, which the march requires. It is inferred and allocates nothing.

save_solution(SS, IC, Ψ, efield, p, t, ts)

Write the state of the cache IC at the time t into the save ts of the sample p of the snapshots SS: the positions and the velocities lifted by Ψ, the potential's coefficients, and the three diagnostics. The accelerations SS.A are not written.

The kinetic energy K = ½∑ⱼwⱼvⱼ² and the scalar momentum M = ∑ⱼwⱼvⱼ are those of one component, nd == 1, which the march requires. It is inferred and allocates nothing.

saved_indices(nₛ)

The three saved states of a sweep of nₛ saved times that a plot of one trend uses: the final one nₛ, the one nearest three quarters of the way, div(3(nₛ - 1), 4) + 1, and the one nearest halfway, div(nₛ - 1, 2) + 1. The remainders take the floor, so the index never reaches past nₛ. For nₛ ≤ 3 the last two indices are the same state: saved_indices(3) == (3, 2, 2).

Throws

• ArgumentError when nₛ < 2: a sweep saves the initial state and at least one more, and a single saved time has no trend to plot.

sodeensemble(model, members, tspan, tstep)

The training data of model as a GeometricEquations.SODEEnsemble, one sample per ParticleDistribution in members.

The initial condition of a member m is (q = copy(m.particles.z),), and its parameters are its own weights, (w = copy(m.particles.w),), because an importance-sampled draw gives each member different weights. Both are copies, so a later change to m does not reach the problem. The invariant h is the VlasovHamiltonian of model.

GeometricIntegrators.integrate(ens, Strang()) solves every member into one EnsembleSolution, which h5save and h5load round-trip through the HDF5 extension of GeometricSolutions. The invariant is evaluated along a stored solution once per saved state, the initial state included, with

compute_invariant(sol.t, sol[i][:q], parameters(ens)[i], invariants(equation(ens)).h)

integrate itself evaluates no invariant.

All members share the potential and the work buffer of model. Solve them one after another, as integrate does: members solved concurrently read each other's fields, with no error.

An empty member list, and a member whose weights and state have different particle counts, throw an ArgumentError. The model and the members have one element type T; a member of another element type raises a MethodError.

trajectory_errors(Xh, Vh, Xr, Vr)

The errors of a full-order state (Xh, Vh) against a reduced state (Xr, Vr): the NamedTuple (x2, xmax, v2, vmax) of the position and the velocity errors. x2 is the root mean square of the position difference over the particle axis, the root sum of squares divided by √nₚ, so it does not change when the particles are duplicated; xmax is the largest absolute position difference, the $L^\infty$ error. v2 and vmax are the same for the velocities.

Each argument holds one saved state per (ts, p), as a (nₚ, nₛ, nparam) array or as a (nₚ, nₛ) matrix, and each result has the trailing shape, (nₛ, nparam) or (nₛ,), of the element type of its argument. NaN and Inf in the input propagate into the result.

The two arrays of one pair must have the same size; a difference raises a DimensionMismatch naming both.

update!(f::PoissonField, x, w, t)

Deposit the positions x with the weights w into the right-hand side of f.potential and re-solve it in place. The time t is accepted for the ElectricField interface and ignored: the field has no explicit time dependence.

velocity_moments(dist::GridDistribution)

The velocity moments of dist at each x-node, as the named tuple (density, momentum, energy) of three vectors of length nx:

\[n_i = \sum_j f_{ij} h_v , \qquad p_i = \sum_j v_j f_{ij} h_v , \qquad \varepsilon_i = \frac{1}{2} \sum_j v_j^2 f_{ij} h_v .\]

energy is the kinetic energy density, with the factor $\tfrac{1}{2}$. The sums are the rectangle rule of _apply_∫dv! and its moment stencils.