pravash 1.2.0

Fluid dynamics simulation — SPH, Euler/Navier-Stokes, shallow water, buoyancy, drag, vortex dynamics
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
133
134
135
136
137
138
139
140
141
142
143
144
145
146

//! Vortex dynamics — vorticity, circulation, turbulence.

use std::f64::consts::PI;

/// Vorticity from a 2D velocity field (curl of velocity).
///
/// ω = ∂vy/∂x - ∂vx/∂y
#[inline]
#[must_use]
pub fn vorticity_2d(dvydx: f64, dvxdy: f64) -> f64 {
    dvydx - dvxdy
}

/// Lamb-Oseen vortex — decaying vortex with viscous core.
///
/// Returns tangential velocity at radius r from center.
/// Γ = circulation, ν = kinematic viscosity, t = time since creation.
#[inline]
#[must_use]
pub fn lamb_oseen_velocity(circulation: f64, radius: f64, viscosity: f64, time: f64) -> f64 {
    if radius <= 0.0 || time <= 0.0 || viscosity <= 0.0 {
        return 0.0;
    }
    let r2 = radius * radius;
    let core = 4.0 * viscosity * time;
    (circulation / (2.0 * PI * radius)) * (1.0 - (-r2 / core).exp())
}

/// Rankine vortex — solid body rotation inside core, irrotational outside.
///
/// v(r) = Γ·r/(2π·R²)  for r ≤ R (solid body)
/// v(r) = Γ/(2π·r)      for r > R (irrotational)
#[inline]
#[must_use]
pub fn rankine_velocity(circulation: f64, radius: f64, core_radius: f64) -> f64 {
    if radius <= 0.0 {
        return 0.0;
    }
    if radius <= core_radius {
        circulation * radius / (2.0 * PI * core_radius * core_radius)
    } else {
        circulation / (2.0 * PI * radius)
    }
}

/// Enstrophy — measure of total vorticity (integral of ω²).
///
/// For a discrete 2D field: Σ ω² × dx²
#[must_use]
pub fn enstrophy(vorticity_field: &[f64], dx: f64) -> f64 {
    let dx2 = dx * dx;
    vorticity_field.iter().map(|&w| w * w * dx2).sum()
}

/// Kolmogorov microscale — smallest turbulent eddy size.
///
/// η = (ν³/ε)^(1/4)
/// where ε = turbulent dissipation rate.
#[inline]
#[must_use]
pub fn kolmogorov_scale(viscosity: f64, dissipation_rate: f64) -> f64 {
    if viscosity <= 0.0 || dissipation_rate <= 0.0 {
        return f64::INFINITY;
    }
    // Use log form for numeric stability: η = exp(0.75·ln(ν) - 0.25·ln(ε))
    (0.75 * viscosity.ln() - 0.25 * dissipation_rate.ln()).exp()
}

#[cfg(test)]
mod tests {
    use super::*;

    const EPS: f64 = 1e-10;

    #[test]
    fn test_vorticity_2d() {
        assert!((vorticity_2d(1.0, 0.0) - 1.0).abs() < EPS);
        assert!((vorticity_2d(0.0, 1.0) - (-1.0)).abs() < EPS);
        assert!((vorticity_2d(2.0, 2.0)).abs() < EPS);
    }

    #[test]
    fn test_lamb_oseen_far_field() {
        let v = lamb_oseen_velocity(1.0, 10.0, 0.001, 1.0);
        // Far from core, approaches Γ/(2πr)
        let irrotational = 1.0 / (2.0 * PI * 10.0);
        assert!((v - irrotational).abs() < 0.01);
    }

    #[test]
    fn test_lamb_oseen_zero_radius() {
        assert!(lamb_oseen_velocity(1.0, 0.0, 0.001, 1.0).abs() < EPS);
    }

    #[test]
    fn test_lamb_oseen_zero_time() {
        assert!(lamb_oseen_velocity(1.0, 1.0, 0.001, 0.0).abs() < EPS);
    }

    #[test]
    fn test_rankine_inside_core() {
        let v = rankine_velocity(2.0 * PI, 0.5, 1.0);
        // Inside: v = Γr/(2πR²) = 2π·0.5/(2π·1) = 0.5
        assert!((v - 0.5).abs() < EPS);
    }

    #[test]
    fn test_rankine_outside_core() {
        let v = rankine_velocity(2.0 * PI, 2.0, 1.0);
        // Outside: v = Γ/(2πr) = 2π/(2π·2) = 0.5
        assert!((v - 0.5).abs() < EPS);
    }

    #[test]
    fn test_rankine_continuous_at_core() {
        let gamma = 5.0;
        let r_core = 1.0;
        let v_inside = rankine_velocity(gamma, r_core, r_core);
        let v_outside = rankine_velocity(gamma, r_core + EPS, r_core);
        assert!((v_inside - v_outside).abs() < 0.01);
    }

    #[test]
    fn test_enstrophy_zero() {
        assert!(enstrophy(&[0.0, 0.0, 0.0], 0.1).abs() < EPS);
    }

    #[test]
    fn test_enstrophy_positive() {
        let field = vec![1.0, -1.0, 2.0, 0.5];
        let e = enstrophy(&field, 0.1);
        assert!(e > 0.0);
    }

    #[test]
    fn test_kolmogorov_scale() {
        let eta = kolmogorov_scale(1e-6, 0.01);
        assert!(eta > 0.0);
        assert!(eta.is_finite());
    }

    #[test]
    fn test_kolmogorov_zero_dissipation() {
        assert!(kolmogorov_scale(1e-6, 0.0).is_infinite());
    }
}