flux: use nearest sampling to stabilize fluid for short timesteps

flux/shader/adjust_advection.comp.wgsl

@@ -25,11 +25,11 @@ fn main(
   let velocity = textureLoad(velocity_texture, global_id.xy, 0).xy;
 
   let size = vec2<f32>(textureDimensions(velocity_texture));
-  let sample_position = vec2<f32>(global_id.xy);
+  let advected_position  = (vec2<f32>(global_id.xy) + 1.0) - uniforms.timestep * velocity;
   // NOTE: Using floor here produces very different results from the GL version.
-  let pos2 = (vec2f(floor(sample_position.x + 1.0), floor(sample_position.y + 1.0)));
-  let adv = floor(sample_position + 1.0) - uniforms.timestep * velocity;
-  let min_max_sampling_position = (0.5 + (adv)) / size;
+  // Using floor individually on each dimensions also seems to produce different results!
+  // floor produces a distinct diagonal bias, moving fluid in waves to/from the upper right corner.
+  let min_max_sampling_position = (0.5 + round(advected_position)) / size;
   let l = textureSampleLevel(velocity_texture, linear_sampler, min_max_sampling_position, 0.0, vec2<i32>(-1, 0)).xy;
   let r = textureSampleLevel(velocity_texture, linear_sampler, min_max_sampling_position, 0.0, vec2<i32>(1, 0)).xy;
   let b = textureSampleLevel(velocity_texture, linear_sampler, min_max_sampling_position, 0.0, vec2<i32>(0, -1)).xy;

flux/shader/advect.comp.wgsl

@@ -35,7 +35,7 @@ fn main(
   let size = vec2<f32>(textureDimensions(velocity_texture));
   let sample_position = vec2<f32>(global_id.xy);
 
-  let advected_position = ((sample_position + 0.5) - direction.direction * uniforms.timestep * velocity) / size;
+  let advected_position = ((sample_position + 0.5) + direction.direction * uniforms.timestep * velocity) / size;
   let decay = 1.0 + uniforms.dissipation * uniforms.timestep;
   let new_velocity = textureSampleLevel(velocity_texture, linear_sampler, advected_position, 0.0).xy / decay;
   textureStore(out_texture, global_id.xy, vec4<f32>(new_velocity, 0.0, 0.0));

flux/shader/diffuse.comp.wgsl

@@ -10,6 +10,7 @@ struct FluidUniforms {
 
 @group(0) @binding(0) var<uniform> uniforms: FluidUniforms;
 @group(0) @binding(1) var linear_sampler: sampler;
+@group(0) @binding(1) var nearest_sampler: sampler;
 
 @group(1) @binding(0) var velocity_texture: texture_2d<f32>;
 @group(1) @binding(1) var out_texture: texture_storage_2d<rg32float, write>;
@@ -23,10 +24,10 @@ fn main(
 
   let size = textureDimensions(velocity_texture, 0);
   let sample_position = (vec2<f32>(global_id.xy) + 0.0) / vec2<f32>(size);
-  let l = textureSampleLevel(velocity_texture, linear_sampler, sample_position, 0.0, vec2<i32>(-1, 0)).xy;
-  let r = textureSampleLevel(velocity_texture, linear_sampler, sample_position, 0.0, vec2<i32>(1, 0)).xy;
-  let b = textureSampleLevel(velocity_texture, linear_sampler, sample_position, 0.0, vec2<i32>(0, -1)).xy;
-  let t = textureSampleLevel(velocity_texture, linear_sampler, sample_position, 0.0, vec2<i32>(0, 1)).xy;
+  let l = textureSampleLevel(velocity_texture, nearest_sampler, sample_position, 0.0, vec2<i32>(-1, 0)).xy;
+  let r = textureSampleLevel(velocity_texture, nearest_sampler, sample_position, 0.0, vec2<i32>(1, 0)).xy;
+  let b = textureSampleLevel(velocity_texture, nearest_sampler, sample_position, 0.0, vec2<i32>(0, -1)).xy;
+  let t = textureSampleLevel(velocity_texture, nearest_sampler, sample_position, 0.0, vec2<i32>(0, 1)).xy;
 
   let new_velocity = uniforms.stencil_factor * (l + r + b + t + uniforms.center_factor * velocity);
 

flux/shader/divergence.comp.wgsl

@@ -1,4 +1,4 @@
-@group(0) @binding(0) var linear_sampler: sampler;
+@group(0) @binding(0) var nearest_sampler: sampler;
 @group(0) @binding(1) var out_divergence_texture: texture_storage_2d<r32float, write>;
 
 @group(1) @binding(0) var velocity_texture: texture_2d<f32>;
@@ -12,10 +12,10 @@ fn main(
   let size = (textureDimensions(velocity_texture));
   let sample_position = (vec2<f32>(global_id.xy) + 0.0) / vec2<f32>(size);
 
-  let l = textureSampleLevel(velocity_texture, linear_sampler, sample_position, 0.0, vec2<i32>(-1, 0)).x;
-  let r = textureSampleLevel(velocity_texture, linear_sampler, sample_position, 0.0, vec2<i32>(1, 0)).x;
-  let t = textureSampleLevel(velocity_texture, linear_sampler, sample_position, 0.0, vec2<i32>(0, 1)).y;
-  let b = textureSampleLevel(velocity_texture, linear_sampler, sample_position, 0.0, vec2<i32>(0, -1)).y;
+  let l = textureSampleLevel(velocity_texture, nearest_sampler, sample_position, 0.0, vec2<i32>(-1, 0)).x;
+  let r = textureSampleLevel(velocity_texture, nearest_sampler, sample_position, 0.0, vec2<i32>(1, 0)).x;
+  let t = textureSampleLevel(velocity_texture, nearest_sampler, sample_position, 0.0, vec2<i32>(0, 1)).y;
+  let b = textureSampleLevel(velocity_texture, nearest_sampler, sample_position, 0.0, vec2<i32>(0, -1)).y;
 
   let new_divergence = 0.5 * ((r - l) + (t - b));
 

flux/shader/generate_noise.comp.wgsl

@@ -114,7 +114,7 @@ fn main(
 ) {
   let size = vec2<f32>(textureDimensions(out_texture));
   let texel_position = vec2<f32>(global_id.xy) / size;
-  
+
   var noise = vec2(0.0);
   for (var i = 0; i < 3; i = i + 1) {
     let channel = global.channels[i];

flux/shader/solve_pressure.comp.wgsl

@@ -10,6 +10,7 @@ struct FluidUniforms {
 
 @group(0) @binding(0) var<uniform> uniforms: FluidUniforms;
 @group(0) @binding(1) var linear_sampler: sampler;
+@group(0) @binding(2) var nearest_sampler: sampler;
 
 @group(1) @binding(0) var divergence_texture: texture_2d<f32>;
 
@@ -27,10 +28,10 @@ fn main(
   let pressure = textureLoad(pressure_texture, global_id.xy, 0).x;
   let divergence = textureLoad(divergence_texture, global_id.xy, 0).x;
 
-  var l = textureSampleLevel(pressure_texture, linear_sampler, sample_position, 0.0, vec2<i32>(-1, 0)).x;
-  var r = textureSampleLevel(pressure_texture, linear_sampler, sample_position, 0.0, vec2<i32>(1, 0)).x;
-  var b = textureSampleLevel(pressure_texture, linear_sampler, sample_position, 0.0, vec2<i32>(0, -1)).x;
-  var t = textureSampleLevel(pressure_texture, linear_sampler, sample_position, 0.0, vec2<i32>(0, 1)).x;
+  var l = textureSampleLevel(pressure_texture, nearest_sampler, sample_position, 0.0, vec2<i32>(-1, 0)).x;
+  var r = textureSampleLevel(pressure_texture, nearest_sampler, sample_position, 0.0, vec2<i32>(1, 0)).x;
+  var b = textureSampleLevel(pressure_texture, nearest_sampler, sample_position, 0.0, vec2<i32>(0, -1)).x;
+  var t = textureSampleLevel(pressure_texture, nearest_sampler, sample_position, 0.0, vec2<i32>(0, 1)).x;
 
   if (global_id.x == 0u) {
     l = pressure;

flux/shader/subtract_gradient.comp.wgsl

@@ -10,6 +10,7 @@ struct FluidUniforms {
 
 @group(0) @binding(0) var<uniform> uniforms: FluidUniforms;
 @group(0) @binding(1) var linear_sampler: sampler;
+@group(0) @binding(2) var nearest_sampler: sampler;
 
 @group(1) @binding(0) var pressure_texture: texture_2d<f32>;
 @group(1) @binding(1) var out_pressure_texture: texture_storage_2d<r32float, write>;
@@ -27,10 +28,10 @@ fn main(
 
   let pressure = textureLoad(pressure_texture, global_id.xy, 0).x;
 
-  var l = textureSampleLevel(pressure_texture, linear_sampler, sample_position, 0.0, vec2<i32>(-1, 0)).x;
-  var r = textureSampleLevel(pressure_texture, linear_sampler, sample_position, 0.0, vec2<i32>(1, 0)).x;
-  var b = textureSampleLevel(pressure_texture, linear_sampler, sample_position, 0.0, vec2<i32>(0, -1)).x;
-  var t = textureSampleLevel(pressure_texture, linear_sampler, sample_position, 0.0, vec2<i32>(0, 1)).x;
+  var l = textureSampleLevel(pressure_texture, nearest_sampler, sample_position, 0.0, vec2<i32>(-1, 0)).x;
+  var r = textureSampleLevel(pressure_texture, nearest_sampler, sample_position, 0.0, vec2<i32>(1, 0)).x;
+  var b = textureSampleLevel(pressure_texture, nearest_sampler, sample_position, 0.0, vec2<i32>(0, -1)).x;
+  var t = textureSampleLevel(pressure_texture, nearest_sampler, sample_position, 0.0, vec2<i32>(0, 1)).x;
 
   // Enforce the following boundary conditions:
   //

flux/src/render/fluid.rs

@@ -297,7 +297,15 @@ impl Context {
             label: Some("sampler:linear"),
             mag_filter: wgpu::FilterMode::Linear,
             min_filter: wgpu::FilterMode::Linear,
-            mipmap_filter: wgpu::FilterMode::Linear,
+            address_mode_u: wgpu::AddressMode::ClampToEdge,
+            address_mode_v: wgpu::AddressMode::ClampToEdge,
+            ..Default::default()
+        });
+
+        let nearest_sampler = device.create_sampler(&wgpu::SamplerDescriptor {
+            label: Some("sampler:nearest"),
+            mag_filter: wgpu::FilterMode::Nearest,
+            min_filter: wgpu::FilterMode::Nearest,
             address_mode_u: wgpu::AddressMode::ClampToEdge,
             address_mode_v: wgpu::AddressMode::ClampToEdge,
             ..Default::default()
@@ -387,6 +395,13 @@ impl Context {
                         ty: wgpu::BindingType::Sampler(wgpu::SamplerBindingType::Filtering),
                         count: None,
                     },
+                    // nearest_sampler
+                    wgpu::BindGroupLayoutEntry {
+                        binding: 2,
+                        visibility: wgpu::ShaderStages::COMPUTE,
+                        ty: wgpu::BindingType::Sampler(wgpu::SamplerBindingType::Filtering),
+                        count: None,
+                    },
                 ],
             });
 
@@ -439,6 +454,10 @@ impl Context {
                     binding: 1,
                     resource: wgpu::BindingResource::Sampler(&linear_sampler),
                 },
+                wgpu::BindGroupEntry {
+                    binding: 2,
+                    resource: wgpu::BindingResource::Sampler(&nearest_sampler),
+                },
             ],
         });
 
@@ -659,7 +678,7 @@ impl Context {
             entries: &[
                 wgpu::BindGroupEntry {
                     binding: 0,
-                    resource: wgpu::BindingResource::Sampler(&linear_sampler),
+                    resource: wgpu::BindingResource::Sampler(&nearest_sampler),
                 },
                 wgpu::BindGroupEntry {
                     binding: 1,
@@ -671,7 +690,7 @@ impl Context {
         let divergence_pipeline_layout =
             device.create_pipeline_layout(&wgpu::PipelineLayoutDescriptor {
                 label: Some("pipeline_layout:divergence"),
-                bind_group_layouts: &[&&divergence_bind_group_layout, &velocity_bind_group_layout],
+                bind_group_layouts: &[&divergence_bind_group_layout, &velocity_bind_group_layout],
                 push_constant_ranges: &[],
             });