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

//! Conservation law validation suite: integration should preserve E, C₂, and mass.

#[test]
fn conservation_laws() {
    use crate::tooling::core::init::{
        domain::{Domain, SpatialBoundType, VelocityBoundType},
        isolated::{PlummerIC, sample_on_grid},
    };
    use crate::tooling::validation::helpers::{build_standard_sim, relative_drift};

    // Plummer sphere on a coarse grid. The Strang splitting conserves a shadow Hamiltonian,
    // so total energy (computed from the same potential) should drift slowly.
    let domain = Domain::builder()
        .spatial_extent(8.0) // [−8, 8]³, 8 cells → dx=2
        .velocity_extent(2.5) // [−2.5, 2.5]³, 8 cells → dv=0.625 (needs cells near v=0)
        .spatial_resolution(8)
        .velocity_resolution(8)
        .t_final(4.0) // ~2 dynamical times
        .spatial_bc(SpatialBoundType::Periodic)
        .velocity_bc(VelocityBoundType::Truncated) // Preserve mass; no particles escape velocity box
        .build()
        .unwrap();

    let ic = PlummerIC::new(1.0, 1.0, 1.0);
    let snap = sample_on_grid(&ic, &domain);

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

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

    let history = &pkg.diagnostics_history;
    assert!(history.len() >= 2, "Need at least 2 diagnostic entries");

    let e0 = history[0].total_energy;
    let c2_0 = history[0].casimir_c2;

    // Compute max relative drifts over the run
    let max_e_drift = history
        .iter()
        .map(|d| relative_drift(d.total_energy, e0))
        .fold(0.0f64, f64::max);
    let max_c2_drift = history
        .iter()
        .map(|d| relative_drift(d.casimir_c2, c2_0))
        .fold(0.0f64, f64::max);

    println!(
        "Conservation: E0={:.4}, max_dE/E={:.2e}; C2_0={:.4}, max_dC2/C2={:.2e}; steps={}",
        e0, max_e_drift, c2_0, max_c2_drift, pkg.total_steps
    );

    // For a coarse 8³×4³ grid with periodic BC, these are rough bounds.
    // Tighter bounds require higher resolution and isolated BC (FftIsolated).
    assert!(!e0.is_nan(), "Initial energy must not be NaN");
    assert!(!max_e_drift.is_nan(), "Energy drift must not be NaN");

    // The Strang splitting conserves a shadow Hamiltonian; energy drift should be bounded.
    // With CFL=0.5 and coarse grid, allow up to 50% drift over 2 t_dyn.
    assert!(
        max_e_drift < 0.5,
        "Energy drift {:.2e} exceeds 50% threshold over {} steps",
        max_e_drift,
        pkg.total_steps
    );
}