FLIP fluid simulation example

[keptsecret]

Trial from the free task list for contributing to Nabla.

Performs a real-time FLIP fluid simulation, using compute shaders.

[devshgraphicsprogramming]

again, 3D dispatches!

what you're doing right now, is using the worst way to map 1D index to texel addresses causing cache misses when sampling neighbours

You don't benefit from workgroup touching neighbouring cells and caching, because you're doing 128 voxels in a straigh-ish line, which means you're tapping 1170 distinct cell values, a ratio of over 9.14

If you were using a 512 workgroup then you'd tap 4626 distinct cells giving you a ratio of 9.04

However if you used a 8x8x8 workgroup in a 3D dispatch, you'd only tap 1000 distinct cells, giving you 350% reduction in cache pressure!

[devshgraphicsprogramming]

Also then you could you know.. preload everything to groupshared [10][10][10] and sample from there

[devshgraphicsprogramming]

this has somewhat bad load balancing, there is another way to do this, ATOMICS

Basically you'd need to clear the velocity buffer to 0 in some previous dispatch and in updateFluidCells (since you're going through particles anyway) you'd do atomicAdd of the weighted velocity to the cell.

why this works?

• you're supposed to split the velocity buffer into mono-channel anyway (3 images one per component)
• emulated float atomic ops can be done via CAS-loops on R32_UINT underlying
• there are also extensions for atomic float and atomic float16 (would have to check if we list them in SFeatures or SLimits of the physical device) if you feel like specializing with device traits

E.g. CAS loop atomic

float32_t oldValue = 0.f;
do
{
   float32_t expectedValue = oldValue;
   float32_t newValue = oldValue+weightedVelocity;
   uint32_t oldBitpattern;
   InterlockedCompareExchange(velocity[component][cIdx],asuint(expectedValue),asuint(newValue),oldBitpattern);
   oldValue = asfloat(oldBitpattern);
} while (oldValue!=expectedValue);

[devshgraphicsprogramming]

then if all the weighted accumulation happensin updateFluidCells this dispatch only needs to enforce the boundary condition so can be kernel-fusioned with the updateNeighbour dispatch

[devshgraphicsprogramming]

Yep there should be device_traits for shaderBufferFloat16AtomicAdd, shaderBufferFloat32AtomicAdd, shaderImageFloat32AtomicAdd if you compile your code with ILogicalDevice

[devshgraphicsprogramming]

your local histogram was so small that you could have prefix summed it in shared memory with a workgroup right away just like our counting sort example.

anyway with a Global sort over the whole image, that won't be possible anymore

[devshgraphicsprogramming]

this could have been kernel fusioned with updateNeighbor

[devshgraphicsprogramming]

this could be kernel fused, if you used 8x8x8 workgroup, preloaded 10x10x10 (1 cell skirt) updated the axis cell in smem and then did the neighbor axis update on the inner 8x8x8 of the 10x10x10

See https://youtu.be/Ol_sHFVXvC0?si=RrP91pPpMzODMEDi&t=1585

[devshgraphicsprogramming]

I'm pretty sure this can be kernel fused.

Also why do you have so many copies of velocity laying around?
I can count at least 4, do they all need to be live at the same time!?

[devshgraphicsprogramming]

this can be kernel fused with the updateVelocity and diffusion pipeline if you preaload a 1 cell border around 8x8x8 into shared memory and do some redundant work

[devshgraphicsprogramming]

a 8x8x8 dispatch has 512 threads.

Recent Nvidia GPUs can have 2,3 or 4 of those resident on an SM at the same time.

Vulkan requires implementations supports at least 16kb of shared memory per workgroup, but in practice all desktop GPUs let you use 32kb.

It just so happens that the old NV GPU arch (Pascal) that can have max 2 workgroups co-resident has at least 64kb of smem, so enough to support 2 workgroups using 32kb each, and newer GPUs have even more.

You are accessing:

• divergence scalar
• pressure scalar
• solid testing

32kb is enough to store 5461 cells of this data if you use float16_t, and 3276 if you use float32_t for the scalars.

These mean 17^3 or 14^3 grids are possible to keep in shared memory. *

Meaning you can do (17-8)/2=4 or (14-8)/2=3 iterations of pressure solving in a single dispatch.

• How is it possible to do only have one copy of your buffers without ping-ponging? Red-Black Ordering!
  https://www.youtube.com/watch?app=desktop&v=giTZ89q-Bpk

[devshgraphicsprogramming]

What I want you do do:

1. Share Descriptor Sets and Layouts between related dispatches (esp the ones I said could be kernel fused) so there's less of them
2. Don't transition your 3D images back and forth between GENERAL and READ_ONLY
3. Make sure everything that's a "grid" is a 3D Image (one image per vector component, RGBA formats waste alpha) in a descriptor array binding and not a Buffer
4. Use 3D dispatches and reduce the 1D -> 3D address conversions (especially ones with integer divisions and modulo shouldn't happen, you will never need to go cell -> particle)
5. Use float atomics (you can require the feature) to accumulate weighted particle velocities and get rid of the particle sorts and per grid cell particle lists
6. Move from SSBO and Dynamic UBO to BDA and make sure to use SoA for particles
7. Cut down/elimanate useless velocity buffer copies (you have like 4 floating about) which don't need to be live at the same time
8. See if you can leverage red-black ordering for some of your grid updates so you don't need an in and out buffer (esp the pressure solver)

The shared memory and kernel fusion improvements we can leave for someone else as their recruitment task in the future. Unless you find some the fusions as trivial as I do.

Nsight profiling might be nice to do before/after changes to see if we're winning and how much.

P.S. The only time you should use a 1D dispatch is when you're going over the particle list, so you can map particleIx = gl_GlobalInvocationID.x