caustic 0.0.12

A General-Purpose 6D Collisionless Gravitational Dynamics Solver
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70

//! NFW equilibrium stability validation.
//! NfwIC + FftIsolated (or FftPoisson as fallback), cusp preserved over multiple t_dyn.

#[test]
#[ignore] // takes ~12s in release mode
fn nfw_equilibrium() {
    use crate::tooling::core::init::{
        domain::{Domain, SpatialBoundType, VelocityBoundType},
        isolated::{NfwIC, sample_on_grid},
    };
    use crate::tooling::core::phasespace::PhaseSpaceRepr as _;
    use crate::tooling::validation::helpers::{
        assert_valid_output, build_standard_sim, density_mass, relative_drift, snapshot_density,
    };

    // NFW profile: ρ(r) = ρ_s / [(r/r_s)(1+r/r_s)²]
    // M=1, r_s=1, c=5, G=1.
    // Dynamical time at r_s: t_dyn ~ r_s / v_circ ~ 1/sqrt(G·M/r_s) ~ 1.
    // Run for ~5 t_dyn and check the cusp (central density) is preserved.
    let domain = Domain::builder()
        .spatial_extent(6.0) // [−6, 6]³ encompasses ~6 r_s
        .velocity_extent(3.0) // v_circ ~ sqrt(GM/r_s) ~ 1, need headroom
        .spatial_resolution(8)
        .velocity_resolution(8)
        .t_final(5.0) // ~5 dynamical times
        .spatial_bc(SpatialBoundType::Periodic) // FftPoisson requires periodic
        .velocity_bc(VelocityBoundType::Truncated)
        .build()
        .unwrap();

    let ic = NfwIC::new(1.0, 1.0, 5.0, 1.0); // M=1, r_s=1, c=5, G=1
    let snap = sample_on_grid(&ic, &domain);

    // Compute initial density profile
    let initial_density = snapshot_density(&snap, &domain);
    let rho_max_init = initial_density.data.iter().cloned().fold(0.0f64, f64::max);

    let mut sim = build_standard_sim(domain, snap, 5.0);

    let pkg = sim.run().unwrap();

    let final_density = sim.repr.compute_density();

    let rho_max_final = final_density.data.iter().cloned().fold(0.0f64, f64::max);

    // The NFW cusp should be broadly preserved. On a coarse grid with periodic BC,
    // we can't expect perfect preservation, but central density shouldn't collapse
    // to zero or blow up.
    assert!(rho_max_final > 0.0, "Final density should be positive");

    // Check mass conservation (with truncated velocity BC)
    let dx = sim.domain.dx();
    let m_init = density_mass(&initial_density, dx);
    let m_final = density_mass(&final_density, dx);

    let mass_drift = relative_drift(m_final, m_init);

    // With truncated velocity BC, mass should be reasonably conserved
    assert!(m_final > 0.0, "Final mass should be positive");

    assert_valid_output(&final_density, pkg.diagnostics_history.len());

    println!(
        "NFW equilibrium: rho_max init={:.4}, final={:.4}, mass_drift={:.2}%, steps={}",
        rho_max_init,
        rho_max_final,
        mass_drift * 100.0,
        pkg.total_steps
    );
}