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
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132

//! Two-stream instability validation.
//! IC: f(x,v) ∝ v² exp(-v²/2) (1 + ε cos(kx)) — two counter-streaming beams.
//! The perturbation should grow and saturate at a finite amplitude.

#[test]
fn two_stream_instability() {
    use crate::tooling::core::algos::lagrangian::SemiLagrangian;
    use crate::tooling::core::algos::uniform::UniformGrid6D;
    use crate::tooling::core::init::domain::{Domain, SpatialBoundType, VelocityBoundType};
    use crate::tooling::core::integrator::TimeIntegrator as _;
    use crate::tooling::core::phasespace::PhaseSpaceRepr as _;
    use crate::tooling::core::poisson::fft::FftPoisson;
    use crate::tooling::core::time::strang::StrangSplitting;

    // Two-stream: f(x,v) = C * v₁² * exp(-v²/2) * (1 + ε cos(k x₁))
    // The v² factor creates two beams moving in ±x₁ direction.
    // With G=1, this is gravitationally unstable for appropriate k.
    let lx = std::f64::consts::PI;
    let lv = 4.0f64;
    let sigma = 1.0f64;
    let g = 1.0f64;
    let epsilon = 0.05f64;
    let k = std::f64::consts::PI / lx; // fundamental mode

    let domain = Domain::builder()
        .spatial_extent(lx)
        .velocity_extent(lv)
        .spatial_resolution(8)
        .velocity_resolution(8)
        .t_final(2.0)
        .spatial_bc(SpatialBoundType::Periodic)
        .velocity_bc(VelocityBoundType::Open)
        .build()
        .unwrap();

    let mut grid = UniformGrid6D::new(domain.clone());
    let dx = domain.dx();
    let dv = domain.dv();
    let [nx1, nx2, nx3, nv1, nv2, nv3] = grid.sizes();

    // Compute normalization for v₁² exp(-v²/2σ²) distribution
    let mut s_norm = 0.0f64;
    for iv1 in 0..nv1 {
        let v1 = -lv + (iv1 as f64 + 0.5) * dv[0];
        for iv2 in 0..nv2 {
            let v2 = -lv + (iv2 as f64 + 0.5) * dv[1];
            for iv3 in 0..nv3 {
                let v3 = -lv + (iv3 as f64 + 0.5) * dv[2];
                let v2sq = v1 * v1 + v2 * v2 + v3 * v3;
                s_norm += v1 * v1 * (-v2sq / (2.0 * sigma * sigma)).exp() * dv[0] * dv[1] * dv[2];
            }
        }
    }
    let c = 1.0 / s_norm.max(1e-30);

    for ix1 in 0..nx1 {
        let x1 = -lx + (ix1 as f64 + 0.5) * dx[0];
        let perturb = 1.0 + epsilon * (k * x1).cos();
        for ix2 in 0..nx2 {
            for ix3 in 0..nx3 {
                for iv1 in 0..nv1 {
                    let v1 = -lv + (iv1 as f64 + 0.5) * dv[0];
                    for iv2 in 0..nv2 {
                        let v2 = -lv + (iv2 as f64 + 0.5) * dv[1];
                        for iv3 in 0..nv3 {
                            let v3 = -lv + (iv3 as f64 + 0.5) * dv[2];
                            let v2sq = v1 * v1 + v2 * v2 + v3 * v3;
                            let f = c * v1 * v1 * (-v2sq / (2.0 * sigma * sigma)).exp() * perturb;
                            let idx = grid.index([ix1, ix2, ix3], [iv1, iv2, iv3]);
                            grid.data[idx] = f.max(0.0);
                        }
                    }
                }
            }
        }
    }

    let rho_init = grid.compute_density();
    let rho_max_init = rho_init.data.iter().cloned().fold(0.0f64, f64::max);
    let rho_min_init = rho_init.data.iter().cloned().fold(f64::MAX, f64::min);
    let amp_init = (rho_max_init - rho_min_init).max(1e-30);

    let poisson = FftPoisson::new(&domain);
    let advector = SemiLagrangian::new();
    let mut integrator = StrangSplitting::new(g);
    let dt = 0.05f64;

    // Track amplitude over time for saturation detection
    let mut max_amp = amp_init;
    for step in 0..40 {
        integrator
            .advance(&mut grid, &poisson, &advector, dt)
            .unwrap();

        if (step + 1) % 10 == 0 {
            let rho = grid.compute_density();
            let rho_max = rho.data.iter().cloned().fold(0.0f64, f64::max);
            let rho_min = rho.data.iter().cloned().fold(f64::MAX, f64::min);
            let amp = rho_max - rho_min;
            if amp > max_amp {
                max_amp = amp;
            }
        }
    }

    let rho_final = grid.compute_density();
    assert!(
        !rho_final.data.iter().any(|x| x.is_nan()),
        "Density contains NaN after two-stream evolution"
    );

    let rho_max_final = rho_final.data.iter().cloned().fold(0.0f64, f64::max);
    let rho_min_final = rho_final.data.iter().cloned().fold(f64::MAX, f64::min);
    let amp_final = rho_max_final - rho_min_final;

    // The two-stream instability should cause the perturbation to grow
    // (at least initially, before possible saturation)
    assert!(
        max_amp > amp_init,
        "Two-stream: perturbation should grow. amp_init={:.4e}, max_amp={:.4e}",
        amp_init,
        max_amp
    );

    println!(
        "Two-stream instability: amp_init={:.4e}, amp_final={:.4e}, max_amp={:.4e}, ratio={:.2}",
        amp_init,
        amp_final,
        max_amp,
        max_amp / amp_init
    );
}