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
#![allow(dead_code)]
use na::RealField;
use crate::object::Fluid;
pub struct TimestepManager<N: RealField> {
cfl_coeff: N,
min_num_substeps: u32,
max_num_substeps: u32,
dt: N,
inv_dt: N,
total_step_size: N,
remaining_time: N,
particle_radius: N,
}
impl<N: RealField> TimestepManager<N> {
pub fn new(particle_radius: N) -> Self {
Self {
cfl_coeff: na::convert(0.4),
min_num_substeps: 1,
max_num_substeps: 10,
particle_radius,
dt: N::zero(),
inv_dt: N::zero(),
total_step_size: N::zero(),
remaining_time: N::zero(),
}
}
fn max_substep(&self, fluids: &[Fluid<N>]) -> N {
let mut max_sq_vel = N::zero();
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(2.0) / max_sq_vel.sqrt() * self.cfl_coeff
}
pub fn reset(&mut self, total_step_size: N) {
self.total_step_size = total_step_size;
self.remaining_time = total_step_size;
}
#[inline]
pub fn is_done(&self) -> bool {
self.remaining_time <= N::default_epsilon()
}
#[inline]
pub fn dt(&self) -> N {
self.dt
}
#[inline]
pub fn inv_dt(&self) -> N {
self.inv_dt
}
#[inline]
pub fn advance(&mut self, fluids: &[Fluid<N>]) {
let substep = self.compute_substep(fluids);
self.dt = substep;
self.inv_dt = if substep.is_zero() {
N::zero()
} else {
N::one() / substep
};
self.remaining_time -= self.dt;
}
fn compute_substep(&self, _fluids: &[Fluid<N>]) -> N {
return self.total_step_size;
}
}