caustic

A 6D Vlasov–Poisson solver framework for collisionless gravitational dynamics.

caustic is a modular, general-purpose library for solving the Vlasov–Poisson equations in full 6D phase space (3 spatial + 3 velocity dimensions). It targets astrophysical problems — dark matter halo formation, galaxy dynamics, tidal streams, stellar system stability — that are traditionally handled by N-body methods but suffer from artificial collisionality and loss of fine-grained phase-space structure.

The library provides a pluggable architecture where the phase-space representation, Poisson solver, time integrator, and initial condition generator can be swapped independently.

Why not N-body?

N-body simulations sample the distribution function with discrete particles. This introduces noise and artificial two-body relaxation that destroys exactly the structures a collisionless solver should preserve: caustic surfaces, thin phase-space streams, and the true velocity distribution at any point. caustic solves the governing equation directly — no particles, no sampling noise, no artificial collisionality.

Quick start

Add to your Cargo.toml:

[dependencies]
caustic = "0.0.9"

Minimal example: Plummer sphere equilibrium

use caustic::prelude::*;
use caustic::{
    FftPoisson, PlummerIC, SemiLagrangian, StrangSplitting,
    SpatialBoundType, VelocityBoundType, sample_on_grid,
};

fn main() -> Result<(), Box<dyn std::error::Error>> {
    let domain = Domain::builder()
        .spatial_extent(20.0)      // [-20, 20]^3 in natural units
        .velocity_extent(3.0)      // [-3, 3]^3
        .spatial_resolution(32)    // 32^3 spatial grid
        .velocity_resolution(32)   // 32^3 velocity grid
        .t_final(50.0)
        .spatial_bc(SpatialBoundType::Periodic)
        .velocity_bc(VelocityBoundType::Open)
        .build()?;

    // Set up a Plummer sphere: mass=1, scale_radius=1, G=1
    let ic = PlummerIC::new(1.0, 1.0, 1.0);
    let snap = sample_on_grid(&ic, &domain);

    let poisson = FftPoisson::new(&domain);
    let mut sim = Simulation::builder()
        .domain(domain)
        .poisson_solver(poisson)
        .advector(SemiLagrangian::new())
        .integrator(StrangSplitting::new(1.0))
        .initial_conditions(snap)
        .time_final(50.0)
        .build()?;

    let exit = sim.run()?;
    exit.print_summary();
    Ok(())
}

Conservation-aware simulation with LoMaC

let mut sim = Simulation::builder()
    .domain(domain)
    .poisson_solver(FftPoisson::new(&domain))
    .advector(SemiLagrangian::new())
    .integrator(StrangSplitting::new(1.0))
    .initial_conditions(snap)
    .time_final(50.0)
    .lomac(true)  // enable local macroscopic conservation
    .build()?;

Architecture

Each solver component is a Rust trait; implementations are swapped independently:

Phase-space representations

Poisson solvers

Time integrators

The PhaseSpaceRepr trait

The central abstraction. All phase-space storage strategies implement this interface:

pub trait PhaseSpaceRepr: Send + Sync {
    fn compute_density(&self) -> DensityField;
    fn advect_x(&mut self, displacement: &DisplacementField, dt: f64);
    fn advect_v(&mut self, acceleration: &AccelerationField, dt: f64);
    fn moment(&self, position: &[f64; 3], order: usize) -> Tensor;
    fn total_mass(&self) -> f64;
    fn casimir_c2(&self) -> f64;
    fn entropy(&self) -> f64;
    fn stream_count(&self) -> StreamCountField;
    fn velocity_distribution(&self, position: &[f64; 3]) -> Vec<f64>;
    fn total_kinetic_energy(&self) -> f64;
    fn to_snapshot(&self, time: f64) -> PhaseSpaceSnapshot;
    fn load_snapshot(&mut self, snap: PhaseSpaceSnapshot);
    fn as_any(&self) -> &dyn Any;
}

Hierarchical Tucker (HT) tensor compression

A uniform 6D grid at N=64 per dimension requires 64⁶ ≈ 7×10¹⁰ cells. The HtTensor representation exploits the balanced binary tree structure of the x-v split to store f(x,v) in O(dNk² + dk³) memory, where k is the representation rank and N is the grid size per dimension.

Construction methods:

• HtTensor::from_full() — compress a full 6D array via hierarchical SVD (HSVD). O(N⁶).
• HtTensor::from_function_aca() — black-box construction via the HTACA algorithm (Ballani & Grasedyck 2013). Builds the HT decomposition by sampling O(dNk) fibers instead of evaluating all N⁶ grid points.

Operations (all in compressed format, never expanding to full):

• compute_density() — O(Nk²) velocity integration via tree contraction
• truncate(eps) — rank-adaptive recompression (orthogonalize + top-down SVD)
• add() — rank-concatenation addition, followed by truncate() to compress
• inner_product() / frobenius_norm() — O(dk⁴) via recursive Gram matrices
• evaluate() — single-point query in O(dk³)
• advect_x() / advect_v() — SLAR (Semi-Lagrangian Adaptive Rank) via HTACA reconstruction

Conservation framework (LoMaC)

The LoMaC (Local Macroscopic Conservation) framework restores exact conservation of mass, momentum, and energy after each time step. Enable via .lomac(true) on the simulation builder.

Components:

• KFVS — Kinetic Flux Vector Splitting macroscopic solver with half-Maxwellian fluxes
• Conservative SVD — Moment-preserving projection via Gram matrix inversion
• Rank-adaptive controller — Conservation-aware truncation tolerance with budget management

Initial conditions

Isolated equilibria (via sample_on_grid()):

• PlummerIC — Plummer sphere via analytic distribution function f(E)
• KingIC — King model via Poisson-Boltzmann ODE (RK4 integration)
• HernquistIC — Hernquist profile via closed-form f(E)
• NfwIC — NFW profile via numerical Eddington inversion

Cosmological:

• ZeldovichSingleMode — single-mode Zel'dovich pancake
• ZeldovichIC — multi-mode Zel'dovich ICs from Gaussian random field (Harrison-Zel'dovich spectrum, FFT-based, reproducible seeding)

Disk dynamics:

• DiskStabilityIC — exponential disk with Shu (1969) distribution function f(E, L_z), Toomre Q stability parameter, azimuthal perturbation modes (bars, spirals)

Multi-body and custom:

• MergerIC — two-body superposition f = f₁ + f₂ with offsets
• TidalIC — progenitor equilibrium model in an external host potential
• CustomIC / CustomICArray — user-provided callable or pre-computed array

Diagnostics

Conserved quantities monitored each timestep via GlobalDiagnostics:

• Total energy (kinetic + potential), momentum, angular momentum
• Casimir C₂, Boltzmann entropy
• Virial ratio, total mass in box
• Density profile (radial binning)

Analysis tools:

• L1/L2/L∞ field norms and error metrics
• ConservationSummary — energy, mass, momentum, C₂ drift tracking
• convergence_table — Richardson extrapolation and convergence order estimation
• CausticDetector — caustic surface detection, first caustic time

I/O and checkpointing

• Binary snapshot save/load (shape + time + data)
• CSV diagnostics (time series of all conserved quantities)
• JSON checkpoint (snapshot + diagnostics history)
• HT tensor checkpoint (tree structure without dense expansion)
• HDF5 (feature-gated): save_snapshot_hdf5, load_snapshot_hdf5, save_ht_checkpoint_hdf5, load_ht_checkpoint_hdf5

Exit conditions

Termination is configurable via the builder API:

Validation suite

163+ tests — run with cargo test --release -- --test-threads=1:

Physics validation

Convergence tests

Solver cross-validation

Plus integration tests, HT tensor/ACA tests (17), conservation framework tests (15), diagnostics tests (10), and solver-specific unit tests.

Equation-to-module mapping

Feature flags

[dependencies]
caustic = { version = "0.0.6", features = ["hdf5"] }

Performance

• Parallelism: rayon data parallelism across UniformGrid6D (density, advection), FftPoisson (FFT axes), TreePoisson (grid walk), SheetTracker (particle advection), Multigrid (residual, prolongation)
• Release profile: fat LTO, codegen-units = 1, target-cpu=native (via .cargo/config.toml)
• Benchmarks: criterion benchmarks (cargo bench), benchmark binary: solver_kernels
• Instrumentation: tracing::info_span! on all hot methods (zero overhead without a subscriber)
• Profiling profile: [profile.profiling] inherits release with debug symbols for perf/samply

Companion: phasma

phasma is a ratatui-based terminal UI that consumes caustic as a library dependency. It provides interactive parameter editing, live diagnostics rendering, density/phase-space heatmaps, energy conservation plots, radial profile charts, rank monitoring, and Poisson solver analysis — all from the terminal. phasma contains no solver logic; it delegates entirely to caustic.

Minimum supported Rust version

Rust edition 2024, targeting stable Rust 1.85+.

License

This project is licensed under the GNU General Public License v3.0. See LICENSE for details.

Citation

If you use caustic in academic work, please cite:

@software{caustic,
  title  = {caustic: A 6D Vlasov--Poisson solver framework},
  url    = {https://github.com/resonant-jovian/caustic},
  year   = {2026}
}