diff --git a/server/src/lib.rs b/server/src/lib.rs index f1620b6680..61b20c5247 100644 --- a/server/src/lib.rs +++ b/server/src/lib.rs @@ -31,6 +31,7 @@ pub mod settings; pub mod state_ext; pub mod sys; #[cfg(not(feature = "worldgen"))] mod test_world; +pub mod windsim; pub mod wiring; // Reexports diff --git a/server/src/sys/mod.rs b/server/src/sys/mod.rs index bdf98b6a22..07f2997ead 100644 --- a/server/src/sys/mod.rs +++ b/server/src/sys/mod.rs @@ -10,6 +10,7 @@ pub mod subscription; pub mod terrain; pub mod terrain_sync; pub mod waypoint; +pub mod windsim; pub mod wiring; use common_ecs::{dispatch, run_now, System}; @@ -32,6 +33,7 @@ pub fn add_server_systems(dispatch_builder: &mut DispatcherBuilder) { dispatch::(dispatch_builder, &[]); dispatch::(dispatch_builder, &[]); dispatch::(dispatch_builder, &[]); + dispatch::(dispatch_builder, &[]); } pub fn run_sync_systems(ecs: &mut specs::World) { diff --git a/server/src/sys/windsim.rs b/server/src/sys/windsim.rs new file mode 100644 index 0000000000..9c45e5e9aa --- /dev/null +++ b/server/src/sys/windsim.rs @@ -0,0 +1,37 @@ +#![allow(dead_code)] +use crate::windsim::WindSim; +use common::{ + comp::{Pos, Vel}, + resources::DeltaTime, +}; +use common_ecs::{Job, Origin, Phase, System}; +use specs::{Read, Write}; +use vek::*; + +/// This system updates the wind grid for the entire map +#[derive(Default)] +pub struct Sys; +impl<'a> System<'a> for Sys { + #[allow(clippy::type_complexity)] + type SystemData = (Read<'a, DeltaTime>, Write<'a, WindSim>); + + const NAME: &'static str = "windsim"; + const ORIGIN: Origin = Origin::Server; + const PHASE: Phase = Phase::Create; + + fn run(_job: &mut Job, (dt, mut windsim): Self::SystemData) { + let wind_sources: Vec<(Pos, Vel)> = vec![( + Pos(Vec3 { + x: 100.0, + y: 100.0, + z: 1.0, + }), + Vel(Vec3 { + x: 100.0, + y: 100.0, + z: 100.0, + }), + )]; + windsim.tick(wind_sources, &dt); + } +} diff --git a/server/src/windsim/fluid.rs b/server/src/windsim/fluid.rs new file mode 100644 index 0000000000..af1a392f51 --- /dev/null +++ b/server/src/windsim/fluid.rs @@ -0,0 +1,397 @@ +use super::types::*; +/// Macro for indexing into a 1D array using 3D coordinates. +macro_rules! IX { + ( $x: expr, $y: expr, $z: expr ) => {{ $x as usize + (X_SIZE + 2) * ($y as usize + (Y_SIZE + 2) * $z as usize) }}; +} + +/// +/// * `b` - The type of border. 1 for vertical vel walls, 2 for side vel walls, +/// 3 for back-front walls, 0 for dens. +fn set_borders(grid: &mut Box<[f32]>, b: u8) { + let max_x = X_SIZE + 1; + let max_y = Y_SIZE + 1; + let max_z = Z_SIZE + 1; + + for ii in 0..max_x { + for kk in 0..max_z { + let ix_top = IX!(ii, max_y, kk); + let ix_top_inset = IX!(ii, max_y - 1, kk); + let ix_bot = IX!(ii, 0, kk); + let ix_bot_inset = IX!(ii, 1, kk); + grid[ix_top] = { + if b == 2 { + -grid[ix_top_inset] + } else { + grid[ix_top_inset] + } + }; + grid[ix_bot] = { + if b == 2 { + -grid[ix_bot_inset] + } else { + grid[ix_bot_inset] + } + }; + } + } + // SIDE WALLS + for jj in 0..max_y { + for kk in 0..max_z { + let ix_left = IX!(0, jj, kk); + let ix_left_inset = IX!(1, jj, kk); + let ix_right = IX!(max_x, jj, kk); + let ix_right_inset = IX!(max_x - 1, jj, kk); + grid[ix_left] = { + if b == 1 { + -grid[ix_left_inset] + } else { + grid[ix_left_inset] + } + }; + grid[ix_right] = { + if b == 1 { + -grid[ix_right_inset] + } else { + grid[ix_right_inset] + } + }; + } + } + // BACK - FRONT WALLS + for ii in 0..max_x { + for jj in 0..max_y { + let ix_front = IX!(ii, jj, 0); + let ix_front_inset = IX!(ii, jj, 1); + let ix_back = IX!(ii, jj, max_z); + let ix_back_inset = IX!(ii, jj, max_z - 1); + grid[ix_front] = { + if b == 3 { + -grid[ix_front_inset] + } else { + grid[ix_front_inset] + } + }; + grid[ix_back] = { + if b == 3 { + -grid[ix_back_inset] + } else { + grid[ix_back_inset] + } + }; + } + } + + // For the 12 edges of the 3d grid + for ii in 1..max_x - 1 { + grid[IX!(ii, 0, 0)] = (grid[IX!(ii + 1, 0, 0)] + + grid[IX!(ii - 1, 0, 0)] + + grid[IX!(ii, 1, 0)] + + grid[IX!(ii, 0, 1)]) + / 4.0; + grid[IX!(ii, max_y, 0)] = (grid[IX!(ii + 1, max_y, 0)] + + grid[IX!(ii - 1, max_y, 0)] + + grid[IX!(ii, max_y - 1, 0)] + + grid[IX!(ii, max_y, 1)]) + / 4.0; + grid[IX!(ii, 0, max_z)] = (grid[IX!(ii + 1, 0, max_z)] + + grid[IX!(ii - 1, 0, max_z)] + + grid[IX!(ii, 0, max_z - 1)] + + grid[IX!(ii, 1, max_z)]) + / 4.0; + grid[IX!(ii, max_y, max_z)] = (grid[IX!(ii + 1, max_y, max_z)] + + grid[IX!(ii - 1, max_y, max_z)] + + grid[IX!(ii, max_y - 1, max_z)] + + grid[IX!(ii, max_y, max_z - 1)]) + / 4.0; + } + // The 8 vertices of the 3d grid + //grid[IX!(0, 0, 0)] = (grid[IX!(1, 0, Z_SIZE - 1)] + grid[IX!(0, 1, Z_SIZE + // - 1)]); +} + +fn linear_solver( + grid: &mut Box<[f32]>, + prev_grid: &Box<[f32]>, + diff_rate: f32, + denominator: f32, + borders: bool, + b: u8, +) { + // For each cell we get contributions from all 6 direct neighbors + let iterations = 1; + + for _ in 0..iterations { + for ii in 1..=X_SIZE { + for jj in 1..=Y_SIZE { + for kk in 1..=Z_SIZE { + let ix = IX!(ii, jj, kk); //Index of current cell + let ix_up = IX!(ii, jj - 1, kk); // 1 row up + let ix_down = IX!(ii, jj + 1, kk); // 1 row down + let ix_back = IX!(ii, jj, kk - 1); // 1 row back + let ix_front = IX!(ii, jj, kk + 1); // 1 row front + grid[ix] = (prev_grid[ix] + + diff_rate + * (grid[ix - 1] + + grid[ix + 1] + + grid[ix_up] + + grid[ix_down] + + grid[ix_back] + + grid[ix_front])) + / denominator; + } + } + } + if borders { + set_borders(grid, b) + } + } +} + +fn trace_backwards(pos: f32, distance_moved: f32, size: f32) -> (f32, f32, f32, f32) { + //estimated position in previous step (not on a cell center - cell centers are + // integers) + let mut prev_pos = pos - distance_moved; + + // is the estimated position past the edges? keep inside + if prev_pos < 0.5 { + prev_pos = 0.5 + } + if prev_pos > size + 0.5 { + prev_pos = size + 0.5 + } + + // find cell centers to left and right (or up and down, etc) + let p0 = prev_pos.floor(); + let p1 = prev_pos.ceil(); + + // how close to cell centers? 0 to 1 + let d1 = prev_pos - p0; + let d0 = 1.0 - d1; + + // return adjacent cell centers and distances to centers + (p0, p1, d0, d1) +} + +/// Diffuse values with neighbors. Depends on viscosity (non-viscosity). +fn diffuse( + grid: &mut Box<[f32]>, + prev_grid: &Box<[f32]>, + dt: f32, + viscosity: f32, + borders: bool, + b: u8, +) { + let diff_rate = dt * viscosity * grid.len() as f32; + if viscosity == 0.0 { + return; + } + linear_solver( + grid, + prev_grid, + diff_rate, + 1.0 + 6.0 * diff_rate, + borders, + b, + ) +} + +/// Move values accross velocity fields +#[allow(clippy::too_many_arguments)] +fn advect( + grid: &mut Box<[f32]>, + prev_grid: &Box<[f32]>, + vx_grid: &Box<[f32]>, + vy_grid: &Box<[f32]>, + vz_grid: &Box<[f32]>, + dt: f32, + borders: bool, + b: u8, +) { + let dt0x = dt * X_SIZE as f32; + let dt0y = dt * Y_SIZE as f32; + let dt0z = dt * Z_SIZE as f32; + + for ii in 1..=X_SIZE { + for jj in 1..=Y_SIZE { + for kk in 1..=Z_SIZE { + // current cell center + let ix = IX!(ii, jj, kk); + + // positions of and distance to adjacent cell centers to previous position of + // current cell center + let (x0, x1, r0, r1) = + trace_backwards(ii as f32, dt0x * vx_grid[ix], X_SIZE as f32); + let (y0, y1, s0, s1) = + trace_backwards(jj as f32, dt0y * vy_grid[ix], Y_SIZE as f32); + let (z0, z1, t0, t1) = + trace_backwards(kk as f32, dt0z * vz_grid[ix], Z_SIZE as f32); + + // all adjacent cell centers surrounding previous position (like vertices of a + // cube) + let ix000 = IX!(x0, y0, z0); + let ix010 = IX!(x0, y1, z0); + let ix100 = IX!(x1, y0, z0); + let ix110 = IX!(x1, y1, z0); + let ix001 = IX!(x0, y0, z1); + let ix011 = IX!(x0, y1, z1); + let ix111 = IX!(x1, y1, z1); + let ix101 = IX!(x1, y0, z1); + + // value of cell is weighed average of the values of the 8 cell centers + grid[ix] = r0 + * (s0 * (t0 * prev_grid[ix000] + t1 * prev_grid[ix001]) + + s1 * (t0 * prev_grid[ix010] + t1 * prev_grid[ix011])) + + r1 * (s0 * (t0 * prev_grid[ix100] + t1 * prev_grid[ix101]) + + s1 * (t0 * prev_grid[ix110] + t1 * prev_grid[ix111])); + } + } + } + if borders { + set_borders(grid, b); + } +} + +// Forces velocity to be mass conserving +#[allow(clippy::too_many_arguments)] +fn project( + vx_grid: &mut Box<[f32]>, + vy_grid: &mut Box<[f32]>, + vz_grid: &mut Box<[f32]>, + prev_x: &mut Box<[f32]>, + prev_y: &mut Box<[f32]>, + prev_z: &mut Box<[f32]>, + borders: bool, +) { + for ii in 1..=X_SIZE { + for jj in 1..=Y_SIZE { + for kk in 1..=Z_SIZE { + let ix = IX!(ii, jj, kk); + let ix_up = IX!(ii, jj - 1, kk); // 1 row up + let ix_down = IX!(ii, jj + 1, kk); // 1 row down + let ix_back = IX!(ii, jj, kk - 1); // 1 row back + let ix_front = IX!(ii, jj, kk + 1); // 1 row front + prev_y[ix] = -0.5 + * ((vx_grid[ix + 1] - vx_grid[ix - 1]) / X_SIZE as f32 + + (vy_grid[ix_down] - vy_grid[ix_up]) / Y_SIZE as f32 + + (vz_grid[ix_front] - vz_grid[ix_back]) / Z_SIZE as f32); + prev_x[ix] = 0.0; + } + } + } + if borders { + set_borders(prev_y, 1); + set_borders(prev_x, 2); + set_borders(prev_z, 3); + } + + //Gauss seidel to compute gradient field x with y (ignoring z right now) + linear_solver(prev_x, prev_y, 1.0, 6.0, false, 0); + + //Substract gradient field + for ii in 1..=X_SIZE { + for jj in 1..=Y_SIZE { + for kk in 1..=Z_SIZE { + let ix = IX!(ii, jj, kk); + let ix_up = IX!(ii, jj - 1, kk); // 1 row up + let ix_down = IX!(ii, jj + 1, kk); // 1 row down + let ix_back = IX!(ii, jj, kk - 1); // 1 row back + let ix_front = IX!(ii, jj, kk + 1); // 1 row front + vx_grid[ix] -= 0.5 * (prev_x[ix + 1] - prev_x[ix - 1]) * X_SIZE as f32; + vy_grid[ix] -= 0.5 * (prev_x[ix_down] - prev_x[ix_up]) * Y_SIZE as f32; + vz_grid[ix] -= 0.5 * (prev_x[ix_front] - prev_x[ix_back]) * Z_SIZE as f32; + } + } + } + if borders { + set_borders(vx_grid, 1); + set_borders(vy_grid, 2); + set_borders(vz_grid, 3); + } +} + +/// Step density +fn step_dens( + dens_grid: &mut Box<[f32]>, + vx_grid: &mut Box<[f32]>, + vy_grid: &mut Box<[f32]>, + vz_grid: &mut Box<[f32]>, + dt: f32, + diff: f32, + borders: bool, +) { + // Make a copy of the dens_grid + let prev_dens_grid = &mut dens_grid.clone(); + + // Swap + let (prev_dens_grid, dens_grid) = (dens_grid, prev_dens_grid); + + // Diffuse + diffuse(dens_grid, prev_dens_grid, dt, diff, borders, 0); + + // Swap + let (prev_dens_grid, dens_grid) = (dens_grid, prev_dens_grid); + + // Advect + advect( + dens_grid, + prev_dens_grid, + vx_grid, + vy_grid, + vz_grid, + dt, + borders, + 0, + ); +} + +/// Step velocity +fn step_vel( + vx_grid: &mut Box<[f32]>, + vy_grid: &mut Box<[f32]>, + vz_grid: &mut Box<[f32]>, + dt: f32, + diff: f32, + borders: bool, +) { + let prev_x = &mut vx_grid.clone(); + let prev_y = &mut vy_grid.clone(); + let prev_z = &mut vz_grid.clone(); + + // Swap grids + let (prev_x, vx_grid) = (vx_grid, prev_x); + let (prev_y, vy_grid) = (vy_grid, prev_y); + let (prev_z, vz_grid) = (vz_grid, prev_z); + + // Diffuse just like with density but with velocity instead + diffuse(vx_grid, prev_x, dt, diff, borders, 1); + diffuse(vy_grid, prev_y, dt, diff, borders, 2); + diffuse(vz_grid, prev_z, dt, diff, borders, 3); + + // For mass conservation before advect + project(vx_grid, vy_grid, vz_grid, prev_x, prev_y, prev_z, borders); + + // Swap grids + let (prev_x, vx_grid) = (vx_grid, prev_x); + let (prev_y, vy_grid) = (vy_grid, prev_y); + let (prev_z, vz_grid) = (vz_grid, prev_z); + + // Advect just like with density + advect(vx_grid, prev_x, prev_x, prev_y, prev_z, dt, borders, 1); + advect(vy_grid, prev_y, prev_x, prev_y, prev_z, dt, borders, 2); + advect(vz_grid, prev_z, prev_x, prev_y, prev_z, dt, borders, 2); + + project(vx_grid, vy_grid, vz_grid, prev_x, prev_y, prev_z, borders); +} + +pub fn step_fluid( + dens_grid: &mut Box<[f32]>, + vx_grid: &mut Box<[f32]>, + vy_grid: &mut Box<[f32]>, + vz_grid: &mut Box<[f32]>, + dt: f32, + viscosity: f32, + borders: bool, +) { + step_dens(dens_grid, vx_grid, vy_grid, vz_grid, dt, viscosity, borders); + step_vel(vx_grid, vy_grid, vz_grid, dt, viscosity, borders); +} diff --git a/server/src/windsim/mod.rs b/server/src/windsim/mod.rs new file mode 100644 index 0000000000..a279e97fe4 --- /dev/null +++ b/server/src/windsim/mod.rs @@ -0,0 +1,95 @@ +#![allow(dead_code)] +pub mod fluid; +mod types; +use common::{terrain::TerrainChunkSize, vol::RectVolSize}; +pub use fluid::step_fluid; +use types::{WindGrid, X_SIZE, Y_SIZE, Z_SIZE}; +use vek::*; +use world::sim::WorldSim; + +use common::{ + comp::{Pos, Vel}, + resources::DeltaTime, +}; +//use common_state::State; + +#[derive(Default)] +pub struct WindSim { + grid: WindGrid, + blocks_per_cell: Vec3, +} + +impl WindSim { + pub fn new(world_sim: &WorldSim) -> Self { + Self { + grid: WindGrid::default(), + blocks_per_cell: cell_size_in_blocks(world_sim), + } + } + + /// Converts world positions, to 3D grid positions. + /// Returns None if out of bounds, for example negative positions. + pub fn world_to_grid(&self, pos: Pos) -> Option> { + if pos + .0 + .map2(self.blocks_per_cell, |pi, si| { + pi >= 0.0 && pi <= (pi / si as f32) + }) + .reduce_and() + { + Some( + pos.0 + .map2(self.blocks_per_cell, |pi, si| pi as usize / si as usize), + ) + } else { + None + } + } + + pub fn tick(&mut self, sources: Vec<(Pos, Vel)>, dt: &DeltaTime) { + for (pos, vel) in sources { + self.grid + .add_velocity_source(pos.0.map(|e| e as usize), vel.0) + } + step_fluid( + &mut self.grid.density, + &mut self.grid.x_vel, + &mut self.grid.y_vel, + &mut self.grid.z_vel, + dt.0, + 0.1, + true, + ); + } +} + +fn cell_size_in_blocks(world_sim: &WorldSim) -> Vec3 { + let world_chunks: Vec2 = world_sim.get_size(); + let blocks_per_chunk: Vec2 = TerrainChunkSize::RECT_SIZE; + + let world_blocks: Vec2 = world_chunks.map2(blocks_per_chunk, |ai, bi| ai * bi); + + let grid_size = Vec3 { + x: X_SIZE as u32, + y: Y_SIZE as u32, + z: Z_SIZE as u32, + }; + let cell_size_xy: Vec2 = world_blocks.map2(grid_size.xy(), |ai, bi| ai / bi as u32); + Vec3 { + x: cell_size_xy.x, + y: cell_size_xy.y, + z: 500, + } +} + +// pub fn init(state: &mut State) { +// let mut grid = WindGrid { +// x_vel: vec![0_f32; SIZE_1D].into_boxed_slice(), +// y_vel: vec![0_f32; SIZE_1D].into_boxed_slice(), +// z_vel: vec![0_f32; SIZE_1D].into_boxed_slice(), +// density: vec![0_f32; SIZE_1D].into_boxed_slice(), +// }; +// +// state.ecs_mut().insert(grid); +// +// } diff --git a/server/src/windsim/types.rs b/server/src/windsim/types.rs new file mode 100644 index 0000000000..3fbb884bb9 --- /dev/null +++ b/server/src/windsim/types.rs @@ -0,0 +1,55 @@ +#![allow(dead_code)] +use vek::*; + +pub const X_SIZE: usize = 96; +pub const Y_SIZE: usize = 96; +pub const Z_SIZE: usize = 4; +pub const SIZE_1D: usize = (X_SIZE + 2) * (Y_SIZE + 2) * (Z_SIZE + 2); + +#[derive(Clone, Debug, PartialEq)] +pub struct WindGrid { + pub x_vel: Box<[f32]>, + pub y_vel: Box<[f32]>, + pub z_vel: Box<[f32]>, + pub density: Box<[f32]>, +} +impl Default for WindGrid { + fn default() -> Self { + WindGrid { + x_vel: vec![0_f32; SIZE_1D].into_boxed_slice(), + y_vel: vec![0_f32; SIZE_1D].into_boxed_slice(), + z_vel: vec![0_f32; SIZE_1D].into_boxed_slice(), + density: vec![0_f32; SIZE_1D].into_boxed_slice(), + } + } +} +impl WindGrid { + // Takes 3D grid position (not a world position) + pub fn add_velocity_source(&mut self, pos: Vec3, vel: Vec3) { + // Convert to 1D grid position + let index = Self::get_index(pos.x, pos.y, pos.z); + self.x_vel[index] = vel.x; + self.y_vel[index] = vel.y; + self.z_vel[index] = vel.z; + } + + // Takes 3D grid position (not a world position) + pub fn add_density_source(&mut self, pos: Vec3, dens: f32) { + let index = Self::get_index(pos.x, pos.y, pos.z); + self.density[index] = dens; + } + + pub fn get_index(x: usize, y: usize, z: usize) -> usize { + x + (X_SIZE + 2) * (y + (Y_SIZE + 2) * z) + } + + // Takes 3D grid position (not a world position) + pub fn get_velocity(&self, pos: Vec3) -> Vec3 { + let index = Self::get_index(pos.x, pos.y, pos.z); + let x = self.x_vel[index] as f32; + let y = self.y_vel[index] as f32; + let z = self.y_vel[index] as f32; + + Vec3 { x, y, z } + } +}