From 60c55aee550695d358d7860f7cdbe0a693015ffb Mon Sep 17 00:00:00 2001 From: Lol3rrr Date: Thu, 25 Dec 2025 22:53:36 +0100 Subject: [PATCH] Add more support --- references/fluidsim.html | 779 ++++++++++++++++++++++++++++++++++++++ src/lib.rs | 2 + src/main.rs | 108 +++--- src/ten_minute_physics.rs | 276 ++++++++++++++ 4 files changed, 1101 insertions(+), 64 deletions(-) create mode 100644 references/fluidsim.html create mode 100644 src/ten_minute_physics.rs 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; + } +}