Struct SpectralV 

pub struct SpectralV {
    pub coefficients: Vec<f64>,
    pub spatial_shape: [usize; 3],
    pub n_modes: usize,
    pub velocity_scale: f64,
    pub domain: Domain,
    pub adaptive_rescale: bool,
    pub hypercollision_nu: f64,
    pub hypercollision_order: usize,
    /* private fields */
}

Spectral-in-velocity representation using orthonormal Hermite function basis.

The distribution function is represented as: f(x,v,t) = sum_{n1,n2,n3} a_{n1,n2,n3}(x,t) * psi_{n1}(v1/sigma) * psi_{n2}(v2/sigma) * psi_{n3}(v3/sigma)

where psi_n(u) = H_n(u) * exp(-u^2/2) / sqrt(2^n * n! * sqrt(pi)) are the orthonormal Hermite functions (eigenfunctions of the Fourier transform and harmonic oscillator).

Layout

Coefficients are stored as a flat array. For each spatial cell (ix, iy, iz), there are n_modes^3 coefficients. The flat index is: spatial_flat * n_modes^3 + m0 * n_modes^2 + m1 * n_modes + m2

where spatial_flat = ix * ny * nz + iy * nz + iz.

Fields

coefficients: Vec<f64>

Hermite expansion coefficients: length = n_spatial * n_modes^3.

spatial_shape: [usize; 3]

Spatial grid dimensions [nx, ny, nz].

n_modes: usize

Number of Hermite modes per velocity dimension.

velocity_scale: f64

Velocity scaling parameter sigma (controls Hermite basis width). Chosen so that the Gaussian weight exp(-v^2/(2*sigma^2)) covers the velocity domain.

domain: Domain

Domain specification.

adaptive_rescale: bool

Whether to enable adaptive velocity rescaling each step.

hypercollision_nu: f64

Hypercollision coefficient (0 = disabled). Damps high modes to suppress recurrence: a_n = exp(-nu * dt * n^{2order}).

hypercollision_order: usize

Order of the hypercollision operator (typically 2 or 3).

Implementations

impl SpectralV

pub fn new(domain: Domain, n_modes: usize) -> Self

Create a new SpectralV with zero coefficients.

The velocity scale sigma is set to Lv / 3 so that 3-sigma covers the velocity half-extent, placing most of the Gaussian weight inside the domain.

pub fn rescale_velocity(&mut self, new_alpha: f64)

Rescale the velocity basis to a new alpha = sigma.

Transforms Hermite coefficients via the connection coefficient matrix that relates Hermite functions at different scales. This preserves the distribution function while optimizing the basis for the current state.

The optimal scale is alpha_opt = sqrt(2 * <v²> / (2*N_modes + 1)).

pub fn optimal_velocity_scale(&self) -> f64

Compute optimal velocity scale from current second velocity moment.

alpha_opt = sqrt(2 * <v²> / (2*N_modes + 1))

pub fn apply_hypercollision(&mut self, dt: f64)

Apply hypercollisional damping to suppress filamentation recurrence.

In mode space: a_n = exp(-nu * dt * (n1^{2p} + n2^{2p} + n3^{2p})) where p = hypercollision_order.

This is diagonal in mode space → O(N_spatial * N_modes³).

pub fn with_positivity_limiter(self, enabled: bool) -> Self

Enable or disable the positivity limiter via builder pattern.

When enabled, after each advect_x and advect_v call the distribution function is evaluated on a velocity quadrature grid, negative values are clipped to zero, and the Hermite coefficients are re-projected from the clipped values (using the orthonormality of Hermite functions so the Gram matrix is identity).

pub fn positivity_violations(&self) -> u64

Return the cumulative number of quadrature points where f was found negative and clipped to zero by the positivity limiter.

pub fn enforce_positivity(&mut self)

Enforce positivity of f by evaluating on a velocity quadrature grid, clipping negatives, and re-projecting onto the Hermite basis.

For each spatial cell:

1. Evaluate f at n_modes^3 quadrature points in velocity space
2. Clip any negative values to zero
3. Re-compute Hermite coefficients as inner products (Gram matrix is identity for orthonormal Hermite functions)
4. Increment violation counter by the number of clipped points

Uses rayon for parallelism over spatial cells.

pub fn hermite_raw(n: usize, x: f64) -> f64

Evaluate the raw (unnormalized) physicist’s Hermite polynomial H_n(x).

Uses the three-term recurrence: H_0(x) = 1 H_1(x) = 2x H_n(x) = 2x * H_{n-1}(x) - 2(n-1) * H_{n-2}(x)

pub fn from_snapshot( snap: &PhaseSpaceSnapshot, n_modes: usize, domain: &Domain, ) -> Self

Project a 6D snapshot onto the Hermite basis.

For each spatial cell (ix, iy, iz), the coefficients are: a_{n1,n2,n3}(x) = integral f(x,v) * psi_{n1}(v1/sigma) * psi_{n2}(v2/sigma) * psi_{n3}(v3/sigma) dv^3 / sigma^3

This integral is approximated by quadrature over the velocity grid in the snapshot.

Trait Implementations

impl PhaseSpaceRepr for SpectralV

fn compute_density(&self) -> DensityField

Compute density rho(x) = integral f(x,v) dv^3.

In the Hermite representation: rho(x) = a_{0,0,0}(x) * sigma^3 * (integral psi_0(u) du)^3

Only the zeroth mode contributes because all higher modes integrate to zero (they are orthogonal to the constant = psi_0 * norm_0, and integral H_n * w du = 0 for n > 0).

fn advect_x(&mut self, _displacement: &DisplacementField, dt: f64)

Spatial drift: advect each coefficient grid a_{n1,n2,n3}(x) by displacement v*dt.

In the Hermite representation, the spatial advection operates on each coefficient grid independently: the Hermite modes do not couple during spatial transport. Each a_{n1,n2,n3}(x) grid is shifted using semi-Lagrangian interpolation.

For a given mode, the “velocity” associated with the drift is the mean velocity of that Hermite function. However, since the Vlasov spatial drift is v * df/dx, and v is not a single value but a distribution, the correct approach in SpectralV is to use the recurrence relation: d/dt a_n = -sigma * (sqrt((n+1)/2) * da_{n+1}/dx + sqrt(n/2) * da_{n-1}/dx)

For implementation simplicity and because this avoids CFL constraints, we use the mode-coupling approach with finite differences for the spatial gradient, but applied dimension-by-dimension.

Actually, the cleaner approach: since displacement already encodes v*dt per spatial cell, we can treat each mode’s coefficient grid as being advected by the mode-weighted displacement. For zeroth mode (Gaussian), the effective velocity is zero; for first mode it is proportional to sigma, etc.

Following the standard SpectralV literature (Schumer & Holloway 1998), we implement the x-advection as mode coupling: The spatial derivative couples modes: v * psi_n(v/sigma) = sigma * (sqrt((n+1)/2) * psi_{n+1} + sqrt(n/2) * psi_{n-1})

For simplicity and correctness with the existing semi-Lagrangian infrastructure, we shift each coefficient grid by a uniform displacement derived from the displacement field (which encodes v*dt averaged over the domain). This is exact for uniform advection; for non-uniform v fields, we use the displacement field values directly.

fn advect_v(&mut self, acceleration: &AccelerationField, dt: f64)

Velocity kick: acceleration couples adjacent Hermite modes.

The recurrence relation for Hermite functions under velocity translation gives: da_n/dt = (g/sigma) * [sqrt(n/2) * a_{n-1} - sqrt((n+1)/2) * a_{n+1}]

This is applied per velocity dimension, per spatial cell. The coupling is tridiagonal in mode space and exact for constant acceleration within a cell. We use forward Euler time integration of the mode coupling.

fn moment(&self, position: &[f64; 3], order: usize) -> Tensor

Compute velocity moment at a given spatial position.

• Order 0: density rho = a_{0,0,0} * normalization
• Order 1: mean velocity <v_i> = a_{e_i} * sigma * sqrt(2) * norm / rho where e_i is the unit vector with 1 in dimension i (mode (1,0,0), (0,1,0), or (0,0,1))
• Order 2: velocity dispersion tensor (reconstructed from coefficients)

fn total_mass(&self) -> f64

Total mass M = integral f dx^3 dv^3.

Since only the zeroth Hermite mode contributes to the velocity integral, M = sum_x a_{0,0,0}(x) * norm * dx^3.

fn casimir_c2(&self) -> f64

Casimir invariant C_2 = integral f^2 dx^3 dv^3.

By Parseval’s theorem for the Hermite expansion: C_2 = sum_x sum_{n1,n2,n3} |a_{n1,n2,n3}(x)|^2 * dx^3 * sigma^3

The sigma^3 factor comes from the change of variable v -> u = v/sigma and the fact that the psi_n are orthonormal in u-space.

fn entropy(&self) -> f64

Entropy S = -integral f * ln(f) dx^3 dv^3.

No closed form exists in Hermite space; we reconstruct f on a velocity quadrature grid and integrate numerically.

fn stream_count(&self) -> StreamCountField

Stream count: number of peaks in the marginal f(v_x | x) at each spatial cell.

Reconstructs f on a 1D velocity grid (marginalised over v_y, v_z) and counts peaks.

fn velocity_distribution(&self, position: &[f64; 3]) -> Vec<f64>

Extract the local velocity distribution f(v | x) at a given spatial position.

Returns f on the velocity grid defined by the domain’s velocity resolution.

fn total_kinetic_energy(&self) -> Option<f64>

Total kinetic energy T = (1/2) integral f * v^2 dx^3 dv^3.

Using the Hermite recurrence for v^2: v_d^2 = sigma^2 * (1/2 + u_d^2) u^2 * psi_n(u) = sqrt(n(n-1)/4) * psi_{n-2} + (n + 1/2) * psi_n + sqrt((n+1)(n+2)/4) * psi_{n+2}

So integral v_d^2 * psi_{m_d}(u_d) * psi_0(u_other) du^3 is nonzero only for specific mode indices. For simplicity, we compute via quadrature.

fn to_snapshot(&self, time: f64) -> Option<PhaseSpaceSnapshot>

Reconstruct the full 6D phase-space snapshot from Hermite coefficients.

fn set_progress(&mut self, p: Arc<StepProgress>)

Attach shared progress state for intra-phase cell-level reporting. Implementations can use this to report progress via StepProgress::set_intra_progress.

fn as_any(&self) -> &dyn Any

Downcast to concrete type for implementation-specific queries (e.g. HT rank data).

fn as_any_mut(&mut self) -> &mut dyn Any

Mutable downcast for in-place modification (e.g. BUG integrator leaf updates). Returns self as &mut dyn Any. Only HtTensor overrides this; the default returns self but concrete downcasts will yield None.

fn load_snapshot( &mut self, snap: PhaseSpaceSnapshot, ) -> Result<(), CausticError>

Replace the current state with data from a dense 6D snapshot.

fn can_materialize(&self) -> bool

Whether full 6D materialization fits in available memory. Default: true. Compressed representations (HT, TT) should override to check shape vs a memory threshold.

fn memory_bytes(&self) -> usize

Approximate memory usage of this representation in bytes. Default returns 0; implementations should override for accurate tracking.