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
#![allow(dead_code)]
use approx::AbsDiffEq;
use num::Zero;
use crate::math::Real;
use crate::object::Fluid;
pub struct TimestepManager {
cfl_coeff: Real,
min_num_substeps: u32,
max_num_substeps: u32,
dt: Real,
inv_dt: Real,
total_step_size: Real,
remaining_time: Real,
particle_radius: Real,
}
impl TimestepManager {
pub fn new(particle_radius: Real) -> Self {
Self {
cfl_coeff: na::convert::<_, Real>(0.4),
min_num_substeps: 1,
max_num_substeps: 10,
particle_radius,
dt: na::zero::<Real>(),
inv_dt: na::zero::<Real>(),
total_step_size: na::zero::<Real>(),
remaining_time: na::zero::<Real>(),
}
}
fn max_substep(&self, fluids: &[Fluid]) -> Real {
let mut max_sq_vel = na::zero::<Real>();
for (v, a) in fluids
.iter()
.flat_map(|f| f.velocities.iter().zip(f.accelerations.iter()))
{
max_sq_vel = max_sq_vel.max((v + a * self.remaining_time).norm_squared());
}
self.particle_radius * na::convert::<_, Real>(2.0) / max_sq_vel.sqrt() * self.cfl_coeff
}
pub fn reset(&mut self, total_step_size: Real) {
self.total_step_size = total_step_size;
self.remaining_time = total_step_size;
}
#[inline]
pub fn is_done(&self) -> bool {
self.remaining_time <= Real::default_epsilon()
}
#[inline]
pub fn dt(&self) -> Real {
self.dt
}
#[inline]
pub fn inv_dt(&self) -> Real {
self.inv_dt
}
#[inline]
pub fn advance(&mut self, fluids: &[Fluid]) {
let substep = self.compute_substep(fluids);
self.dt = substep;
self.inv_dt = if substep.is_zero() {
na::zero::<Real>()
} else {
na::one::<Real>() / substep
};
self.remaining_time -= self.dt;
}
fn compute_substep(&self, _fluids: &[Fluid]) -> Real {
return self.total_step_size;
}
}