diff --git a/references/fluidsim.html b/references/fluidsim.html
new file mode 100644
index 0000000..40b2cdb
--- /dev/null
+++ b/references/fluidsim.html
@@ -0,0 +1,779 @@
+
+
+
+
+
+
+
+ Euler Fluid
+
+
+
+
+
+
+
+
+
+ Streamlines
+ Velocities
+ Pressure
+ Smoke
+ Overrelax
+
+
+
+
+
+
diff --git a/src/lib.rs b/src/lib.rs
index 7f822c3..f27147e 100644
--- a/src/lib.rs
+++ b/src/lib.rs
@@ -1,3 +1,5 @@
+pub mod ten_minute_physics;
+
#[derive(Debug, Clone)]
pub struct Domain {
pub width: usize,
diff --git a/src/main.rs b/src/main.rs
index 5ab91bc..3830b9d 100644
--- a/src/main.rs
+++ b/src/main.rs
@@ -4,68 +4,58 @@ use ratatui::{text::Text, Frame, widgets::canvas::Canvas, widgets::canvas::Point
use fluidsim::{Domain, simulate};
fn main() {
- let width = 100;
+ let width = 140;
let height = 60;
- let mut domain = Domain {
- width: width + 2,
- height: height + 2,
- cells: core::iter::repeat_n((0.0, 0.0), (width+2)*(height+2)).collect(),
- enabled: core::iter::repeat_n(1.0, (width+2)*(height+2)).collect(),
- };
+ let mut t_domain = fluidsim::ten_minute_physics::Domain::new(1000.0, width, height);
+ t_domain.set_s_with(|x, y| {
+ 1.0
+ });
+ t_domain.set_u_with(|x, y| {
+ if x == 0 {
+ 2.0
+ } else {
+ 0.0
+ }
+ });
- let boundary = 10;
- for y in boundary..domain.height-boundary {
- let cell = domain.cells.get_mut(y*domain.width + 1).unwrap();
- cell.0 = 1.0;
- }
-
- for y in 0..domain.height {
- *domain.enabled.get_mut(y*domain.width).unwrap() = 0.0;
- *domain.enabled.get_mut(y*domain.width + domain.width-1).unwrap() = 0.0;
- }
- for x in 0..domain.width {
- *domain.enabled.get_mut(0*domain.width + x).unwrap() = 0.0;
- *domain.enabled.get_mut((domain.height-1)*domain.width + x).unwrap() = 0.0;
- }
+ t_domain.set_s_with(|x, y| {
+ if x > 40 && x < 50 && y > 20 && y < 30 {
+ return 0.0;
+ }
- let (domain_tx, domain_rx) = std::sync::mpsc::channel();
+ 1.0
+ });
+
+ let domain_mtx_1 = std::sync::Arc::new(std::sync::Mutex::new(t_domain.clone()));
+ let domain_mtx = domain_mtx_1.clone();
std::thread::spawn(move || {
- let mut start = domain;
-
- domain_tx.send(Domain {
- width: start.width,
- height: start.height,
- cells: start.cells.clone(),
- enabled: start.enabled.clone(),
- }).unwrap();
+ let mut start = t_domain;
for i in 0.. {
// println!("Iteration {}", i);
- simulate(&mut start);
+ start.reset_pressure();
- if i % 100 == 0 {
- domain_tx.send(Domain {
- width: start.width,
- height: start.height,
- cells: start.cells.clone(),
- enabled: start.enabled.clone(),
- }).unwrap();
+ start.solve_incompressibility(100, 0.01);
+ start.extrapolate();
+ start.advect_vel(0.01);
+ start.advect_smoke(0.01);
+
+ if let Ok(mut mtx) = domain_mtx.try_lock() {
+ *mtx = start.clone();
}
- }
- std::thread::sleep(std::time::Duration::from_millis(10));
+ std::thread::sleep(std::time::Duration::from_millis(10));
+ }
});
- let mut domain = domain_rx.recv().unwrap();
+ let domain_mtx = domain_mtx_1;
+ let mut domain = domain_mtx.lock().unwrap().clone();
let mut terminal = ratatui::init();
loop {
- domain = match domain_rx.try_recv() {
- Ok(d) => d,
- Err(_) => domain,
- };
+ domain = domain_mtx.lock().unwrap().clone();
terminal.draw(|frame| {
draw(frame, &domain)
@@ -80,33 +70,23 @@ fn main() {
ratatui::restore();
}
-fn draw(frame: &mut Frame, domain: &Domain) {
- let canvas = Canvas::default().x_bounds([0.0, domain.width as f64 - 2.0]).y_bounds([0.0, domain.height as f64 - 2.0]).paint(|ctx| {
- for y in 1..domain.height-1 {
- for x in 1..domain.width-1 {
- let own_cell_idx = y * domain.width + x;
- let top_cell_idx = (y-1) * domain.width + x;
- let right_cell_idx = y * domain.width + x + 1;
+fn draw(frame: &mut Frame, domain: &fluidsim::ten_minute_physics::Domain) {
+ let (width, height) = domain.inner_dims();
- let own_cell = domain.cells.get(own_cell_idx).unwrap();
- let top_cell = domain.cells.get(top_cell_idx).unwrap();
- let right_cell = domain.cells.get(right_cell_idx).unwrap();
+ let canvas = Canvas::default().x_bounds([0.0, width as f64]).y_bounds([0.0, height as f64]).paint(|ctx| {
+ let (min_v, max_v) = domain.smoke_iter().fold((f32::MAX, f32::MIN), |(mip, map), v| {
+ (mip.min(v.2), map.max(v.2))
+ });
- let d = right_cell.0 - own_cell.0 + top_cell.1 - own_cell.1;
- // let text = format!("{:.4}", d);
- // ctx.print(-0.5 + x as f64, -0.5 + y as f64, text);
- let color = match d.total_cmp(&0.0) {
- core::cmp::Ordering::Less if d.abs() > 0.0001 => ratatui::style::Color::Rgb(255, 0, 0),
- core::cmp::Ordering::Greater if d.abs() > 0.0001 => ratatui::style::Color::Rgb(0, 255, 0),
- _ => ratatui::style::Color::White,
- };
+ for (x, y, d) in domain.smoke_iter() {
+ let inter = ((d - min_v) / (max_v - min_v)).sqrt();
+ let color = ratatui::style::Color::Rgb((255.0*inter) as u8, (255.0*inter) as u8, 0);
ctx.draw(&Points {
coords: &[(-0.5 + x as f64, -0.5 + y as f64)],
color: color,
})
- }
}
});
frame.render_widget(canvas, frame.area());
diff --git a/src/ten_minute_physics.rs b/src/ten_minute_physics.rs
new file mode 100644
index 0000000..11e3adc
--- /dev/null
+++ b/src/ten_minute_physics.rs
@@ -0,0 +1,276 @@
+#[derive(Clone)]
+pub struct Domain {
+ density: f32,
+ /// The inner width (without the extra padding cells at the edge)
+ width: usize,
+ /// The inner height (without the extra padding cells at the edge)
+ height: usize,
+ u: Vec,
+ v: Vec,
+ s: Vec,
+ p: Vec,
+ /// Smoke field
+ m: Vec,
+}
+
+enum Field {
+ U,
+ V,
+ Smoke,
+}
+
+impl Domain {
+ pub fn new(density:f32, width: usize, height: usize) -> Self {
+ let num_cells = (width + 2) * (height + 2);
+
+ Self {
+ density,
+ width,
+ height,
+ u: vec![0.0; num_cells],
+ v: vec![0.0; num_cells],
+ s: vec![0.0; num_cells],
+ p: vec![0.0; num_cells],
+ m: vec![1.0; num_cells],
+ }
+ }
+
+ pub fn reset_pressure(&mut self) {
+ self.p.fill(0.0);
+ }
+ pub fn set_s_with(&mut self, func: F) where F: Fn(usize, usize) -> f32 {
+ for i in 1..self.width+1 {
+ for j in 1..self.height+1 {
+ self.s[i*(self.height+2)+j] = func(i-1, j-1);
+ }
+ }
+ }
+
+ pub fn set_u_with(&mut self, func: F) where F: Fn(usize, usize) -> f32 {
+ for i in 1..self.width+1 {
+ for j in 1..self.height+1 {
+ self.u[i*(self.height+2)+j] = func(i-1, j-1);
+ }
+ }
+ }
+ pub fn set_v_with(&mut self, func: F) where F: Fn(usize, usize) -> f32 {
+ for i in 1..self.width+1 {
+ for j in 1..self.height+1 {
+ self.v[i*(self.height+2)+j] = func(i-1, j-1);
+ }
+ }
+ }
+
+ pub fn inner_dims(&self) -> (usize, usize) {
+ (self.width, self.height)
+ }
+
+ pub fn div_iter(&self) -> impl Iterator- {
+ let n = self.height + 2;
+
+ (1..self.width+1)
+ .flat_map(|x| core::iter::repeat(x).zip(1..self.height+1))
+ .map(move |(x, y)| {
+ let div = self.u[(x+1)*n + y] - self.u[x*n + y] + self.v[x*n + y + 1] - self.v[x*n+y];
+
+ (x - 1, y - 1, div)
+ })
+ }
+
+ pub fn vel_iter(&self) -> impl Iterator
- {
+ let n = self.height + 2;
+
+ (1..self.width+1)
+ .flat_map(|x| core::iter::repeat(x).zip(1..self.height+1))
+ .map(move |(x, y)| {
+ let u = self.u[x*n+y];
+ let v = self.v[x*n+y];
+
+ (x - 1, y - 1, (u*u + v*v).sqrt())
+ })
+ }
+
+ pub fn pressure_iter(&self) -> impl Iterator
- {
+ let n = self.height + 2;
+
+ (1..self.width+1)
+ .flat_map(|x| core::iter::repeat(x).zip(1..self.height+1))
+ .map(move |(i, j)| {
+ (i - 1, j - 1, self.p[i*n+j])
+ })
+ }
+
+ pub fn smoke_iter(&self) -> impl Iterator
- {
+ let n = self.height + 2;
+
+ (1..self.width+1)
+ .flat_map(|x| core::iter::repeat(x).zip(1..self.height+1))
+ .map(move |(i, j)| {
+ (i - 1, j - 1, self.m[i*n+j])
+ })
+ }
+
+ pub fn solve_incompressibility(&mut self, num_iters: usize, dT: f32) {
+ let n = self.height + 2;
+ let cp = self.density * 1.0 / dT;
+
+ for _ in 0..num_iters {
+ for i in 1..self.width+1 {
+ for j in 1..self.height+1 {
+ if self.s[i*n + j] == 0.0 {
+ continue;
+ }
+
+ let s = self.s[i*n + j];
+ let sx0 = self.s[(i-1)*n + j];
+ let sx1 = self.s[(i+1)*n + j];
+ let sy0 = self.s[i*n + j - 1];
+ let sy1 = self.s[i*n + j + 1];
+
+ let s = sx0 + sx1 + sy0 + sy1;
+ if s == 0.0 {
+ continue;
+ }
+
+ let div = self.u[(i+1)*n + j] - self.u[i*n + j] + self.v[i*n + j + 1] - self.v[i*n+j];
+
+ let p = -div / s;
+ let p = p * 1.9;
+ self.p[i*n+j] += cp * p;
+
+ self.u[i*n + j] -= sx0 * p;
+ self.u[(i+1)*n + j] += sx1 * p;
+ self.v[i*n+j] -= sy0 * p;
+ self.v[i*n + j + 1] += sy1 * p;
+ }
+ }
+ }
+ }
+
+ pub fn extrapolate(&mut self) {
+ let n = self.height + 2;
+
+ for i in 0..self.width+2 {
+ self.u[i*n + 0] = self.u[i*n + 1];
+ self.u[i*n + self.height+1] = self.u[i*n + self.height];
+ }
+ for j in 0..self.height+2 {
+ self.v[0*n + j] = self.u[1*n + j];
+ self.v[(self.width+1)*n + j] = self.u[self.width*n + j];
+ }
+ }
+
+ fn sample_field(&self, x: f32, y: f32, field: Field) -> f32 {
+ let n = self.height + 2;
+ let h = 1.0;
+
+ let h1 = 1.0/h;
+ let h2 = 0.5*h;
+
+ let x = (x.min((self.width+2) as f32 * h)).max(h);
+ let y = (y.min((self.height+2) as f32 * h)).max(h);
+
+ let (f, dx, dy) = match field {
+ Field::U => (&self.u, 0.0, h2),
+ Field::V => (&self.v, h2, 0.0),
+ Field::Smoke => (&self.m, h2, h2),
+ };
+
+ let x0 = ((x-dx)*h1).floor().min(self.width as f32 + 1.0);
+ let tx = ((x-dx) - x0*h) * h1;
+ let x1 = (x0 + 1.0).min(self.width as f32 + 1.0);
+
+ let x0 = x0 as usize;
+ let x1 = x1 as usize;
+
+ let y0 = ((y-dy)*h1).floor().min(self.height as f32 + 1.0);
+ let ty = ((y-dy) - y0*h) * h1;
+ let y1 = (y0 + 1.0).min(self.height as f32 + 1.0);
+
+ let y0 = y0 as usize;
+ let y1 = y1 as usize;
+
+ let sx = 1.0 - tx;
+ let sy = 1.0 - ty;
+
+ let val = sx * sy * f[x0*n + y0]
+ + tx*sy * f[x1*n + y0]
+ + tx*ty * f[x1*n + y1]
+ + sx*ty * f[x0*n + y1];
+
+ return val;
+ }
+
+ pub fn avg_u(&self, i: usize, j: usize) -> f32 {
+ let n = self.height + 2;
+ let u = (self.u[i*n + j-1] + self.u[i*n+j] + self.u[(i+1)*n + j-1] + self.u[(i+1)*n + j]) * 0.25;
+ return u;
+ }
+
+ pub fn avg_v(&self, i: usize, j: usize) -> f32 {
+ let n = self.height + 2;
+ let v = (self.v[(i-1)*n + j] + self.v[i*n+j] + self.v[(i-1)*n + j + 1] + self.v[i*n + j*1]) * 0.25;
+ return v;
+ }
+
+ pub fn advect_vel(&mut self, dt: f32) {
+ let mut newU = self.u.clone();
+ let mut newV = self.v.clone();
+
+ let n = self.height + 2;
+ let h = 1.0;
+ let h2 = 0.5 * h;
+
+ for i in 1..self.width+2 {
+ for j in 1..self.height+2 {
+ // u component
+ if self.s[i*n + j] != 0.0 && self.s[(i-1)*n + j] != 0.0 && j < self.height+1 {
+ let x = i as f32 * h;
+ let y = j as f32 *h + h2;
+ let u = self.u[i*n + j];
+ let v = self.avg_v(i, j);
+ let x = x - dt*u;
+ let y = y - dt*v;
+ let u = self.sample_field(x, y, Field::U);
+ newU[i*n + j] = u;
+ }
+
+ // v component
+ if self.s[i*n + j] != 0.0 && self.s[i*n + j - 1] != 0.0 && i < self.width + 1 {
+ let x = i as f32 * h;
+ let y = j as f32 *h + h2;
+ let u = self.avg_u(i, j);
+ let v = self.v[i*n+j];
+ let x = x - dt*u;
+ let y = y - dt*v;
+ let v = self.sample_field(x, y, Field::V);
+ newV[i*n + j] = v;
+ }
+ }
+ }
+
+ self.u = newU;
+ self.v = newV;
+ }
+
+ pub fn advect_smoke(&mut self, dt: f32) {
+ let mut newM = self.m.clone();
+
+ let n = self.height + 2;
+ let h = 1.0;
+ let h2 = 0.5*h;
+
+ for i in 1..self.width+1 {
+ for j in 1..self.height+1 {
+ let u = (self.u[i*n+j] + self.u[(i+1)*n + j]) * 0.5;
+ let v = (self.v[i*n+j] + self.v[i*n + j+1]) * 0.5;
+ let x = i as f32 *h + h2 - dt*u;
+ let y = j as f32 *h + h2 - dt*v;
+
+ newM[i*n+j] = self.sample_field(x, y, Field::Smoke);
+ }
+ }
+
+ self.m = newM;
+ }
+}