VlasovMethods

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

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}$.

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.

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.

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

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

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.

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.

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! writes those values into a buffer and returns the index of the first, before wrapping; basis_index is what wraps it, and is the identity where the basis is not periodic.

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.

$(\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.