caustic/tooling/core/algos/
lagrangian.rs

1//! Semi-Lagrangian advector. Traces characteristics backwards from each grid point and
2//! interpolates. Eliminates CFL constraint; Δt limited only by accuracy requirements.
3
4use super::super::{advecator::Advector, phasespace::PhaseSpaceRepr, types::*};
5
6/// Semi-Lagrangian advector with cubic spline interpolation.
7pub struct SemiLagrangian {
8    pub interpolation_order: usize,
9}
10
11impl Default for SemiLagrangian {
12    fn default() -> Self {
13        Self::new()
14    }
15}
16
17impl SemiLagrangian {
18    pub fn new() -> Self {
19        Self {
20            interpolation_order: 3,
21        }
22    }
23}
24
25impl Advector for SemiLagrangian {
26    fn drift(&self, repr: &mut dyn PhaseSpaceRepr, dt: f64) {
27        // UniformGrid6D ignores the displacement field; it computes shifts from its velocity grid.
28        let dummy = DisplacementField {
29            dx: vec![],
30            dy: vec![],
31            dz: vec![],
32            shape: [0, 0, 0],
33        };
34        repr.advect_x(&dummy, dt);
35    }
36
37    fn kick(&self, repr: &mut dyn PhaseSpaceRepr, acceleration: &AccelerationField, dt: f64) {
38        repr.advect_v(acceleration, dt);
39    }
40
41    fn step(&self, repr: &mut dyn PhaseSpaceRepr, acceleration: &AccelerationField, dt: f64) {
42        self.kick(repr, acceleration, dt);
43        self.drift(repr, dt);
44    }
45}
46
47/// 4-point Catmull-Rom cubic spline interpolation at fractional position `x` in
48/// periodic array of length `n`. `x` is a fractional cell index in `[0, n)`.
49pub fn cubic_spline_interpolate(data: &[f64], x: f64, n: usize) -> f64 {
50    let i = x.floor() as isize;
51    let t = x - x.floor();
52    let wrap = |j: isize| (j.rem_euclid(n as isize)) as usize;
53    let p0 = data[wrap(i - 1)];
54    let p1 = data[wrap(i)];
55    let p2 = data[wrap(i + 1)];
56    let p3 = data[wrap(i + 2)];
57    let a0 = p1;
58    let a1 = 0.5 * (-p0 + p2);
59    let a2 = p0 - 2.5 * p1 + 2.0 * p2 - 0.5 * p3;
60    let a3 = -0.5 * p0 + 1.5 * p1 - 1.5 * p2 + 0.5 * p3;
61    a0 + a1 * t + a2 * t * t + a3 * t * t * t
62}
63
64/// 4-point Catmull-Rom cubic spline interpolation with clamped (open) boundary.
65fn cubic_spline_interpolate_open(data: &[f64], x: f64, n: usize) -> f64 {
66    let i = x.floor() as isize;
67    let t = x - x.floor();
68    let clamp = |j: isize| j.clamp(0, n as isize - 1) as usize;
69    let p0 = data[clamp(i - 1)];
70    let p1 = data[clamp(i)];
71    let p2 = data[clamp(i + 1)];
72    let p3 = data[clamp(i + 2)];
73    let a0 = p1;
74    let a1 = 0.5 * (-p0 + p2);
75    let a2 = p0 - 2.5 * p1 + 2.0 * p2 - 0.5 * p3;
76    let a3 = -0.5 * p0 + 1.5 * p1 - 1.5 * p2 + 0.5 * p3;
77    a0 + a1 * t + a2 * t * t + a3 * t * t * t
78}
79
80/// 1D semi-Lagrangian shift of a grid line.
81///
82/// - `data`: array of `n` values on a periodic/open 1D grid spanning `[−l, l]`
83/// - `disp`: displacement (physical units) — the whole line shifts by this amount
84/// - `cell_size`: grid spacing (= 2l/n)
85/// - `n`: number of cells
86/// - `l`: half-extent of the domain
87/// - `periodic`: true → wrap-around; false → absorbing (out-of-bounds → 0)
88pub fn sl_shift_1d(
89    data: &[f64],
90    disp: f64,
91    cell_size: f64,
92    n: usize,
93    l: f64,
94    periodic: bool,
95) -> Vec<f64> {
96    (0..n)
97        .map(|i| {
98            // Center of output cell i: -l + (i + 0.5) * cell_size
99            // Departure point (back-trace): center_i - disp
100            let center_i = -l + (i as f64 + 0.5) * cell_size;
101            let departure_phys = center_i - disp;
102            // Fractional cell index of departure point
103            let departure_idx = (departure_phys + l) / cell_size - 0.5;
104            if periodic {
105                cubic_spline_interpolate(data, departure_idx, n)
106            } else {
107                // Absorbing BC: anything outside [−0.5, n−0.5) maps to zero
108                if departure_idx < -0.5 || departure_idx >= n as f64 - 0.5 {
109                    0.0
110                } else {
111                    let clamped = departure_idx.clamp(0.0, n as f64 - 1.0 - 1e-10);
112                    cubic_spline_interpolate_open(data, clamped, n)
113                }
114            }
115        })
116        .collect()
117}
caustic/tooling/core/algos/lagrangian.rs

caustic/tooling/core/algos/
lagrangian.rs
