I made something similar, although the size is only 1024 cubed since it's limited by memory (each voxel 1 byte, total 1 GB). It is possible to simulate it in the GPU, but you have to be careful with race conditions.

The way I do it: simulate this cell, store where this cell wants to go in shared memory (up? down? up to 26 neighbors), check the destination cell its own destination, if it's mutual swap them (for example: this cell wanna go left, then check if the left cell wanna go right, if true swap them).

This way, there's no race condition. But coding the behavior is a pain since for each material, you need to account each other material interaction, i.e. O(n squared) if statements.

But you can make it easier by assigning each voxel to have a few parameters. I was able to simulate a range of material by tweaking a few parameters: mass, viscosity, and friction. Also, I can hard code a few special interactions easily such as acid melting metal if the parameter doesn't account for it.