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

//! Shared helper functions and macros for time integrators.
//!
//! Eliminates duplication of the dynamical-time CFL estimate, the
//! density→potential→acceleration Poisson pipeline, timing instrumentation,
//! progress reporting, and SpectralV hypercollision dispatch.

use super::super::{
    phasespace::PhaseSpaceRepr,
    solver::PoissonSolver,
    types::{AccelerationField, DensityField, PotentialField},
};

/// Compute the maximum stable timestep from the dynamical time.
///
/// dt = cfl_factor / sqrt(G * rho_max)
///
/// Returns a large sentinel (1e10) if the density or G is non-positive.
/// Used by all splitting integrators (Strang, Yoshida, Lie, BM4, RKN6, etc.).
/// The unsplit integrator uses a different CFL based on spatial/velocity grids.
pub fn dynamical_timestep(repr: &dyn PhaseSpaceRepr, g: f64, cfl_factor: f64) -> f64 {
    let density = repr.compute_density();
    let rho_max = density.data.iter().cloned().fold(0.0_f64, f64::max);
    if rho_max <= 0.0 || g <= 0.0 {
        return 1e10;
    }
    let t_dyn = 1.0 / (g * rho_max).sqrt();
    cfl_factor * t_dyn
}

/// Solve the full Poisson pipeline: density → potential → acceleration.
///
/// This is the most common 3-line sequence in time integrators. Returns all
/// three fields so the caller can use them for diagnostics and conservation.
pub fn solve_poisson(
    repr: &dyn PhaseSpaceRepr,
    solver: &dyn PoissonSolver,
    g: f64,
) -> (DensityField, PotentialField, AccelerationField) {
    let density = repr.compute_density();
    let potential = solver.solve(&density, g);
    let acceleration = solver.compute_acceleration(&potential);
    (density, potential, acceleration)
}

/// Apply SpectralV hypercollision damping if the representation is `SpectralV`.
///
/// This is a no-op for all other representation types. Used after each kick
/// sub-step in Strang, Yoshida, and Lie splitting to suppress Gibbs ringing.
pub fn apply_hypercollision_if_spectral(repr: &mut dyn PhaseSpaceRepr, dt: f64) {
    if let Some(spectral) = repr
        .as_any_mut()
        .downcast_mut::<super::super::algos::spectral::SpectralV>()
    {
        spectral.apply_hypercollision(dt);
    }
}

/// Time a block and accumulate milliseconds into a `StepTimings` field.
///
/// Usage: `time_ms!(timings, drift_ms, { advector.drift(repr, dt) });`
macro_rules! time_ms {
    ($timings:expr, $field:ident, $body:expr) => {{
        let _t0 = std::time::Instant::now();
        let _result = $body;
        $timings.$field += _t0.elapsed().as_secs_f64() * 1000.0;
        _result
    }};
}

/// Report integrator phase and sub-step progress to the TUI.
///
/// No-op if progress is `None`. Usage:
/// `report_phase!(self.progress, StepPhase::DriftHalf1, 0, 5);`
macro_rules! report_phase {
    ($progress:expr, $phase:expr, $step:expr, $total:expr) => {
        if let Some(p) = &$progress {
            p.set_phase($phase);
            p.set_sub_step($step, $total);
        }
    };
}

#[allow(unused_imports)]
pub(crate) use report_phase;
#[allow(unused_imports)]
pub(crate) use time_ms;