## Simple Fluid Simulation With Cellular Automata

Last week I couldn’t use my regular dev. machine (broken graphics card), so all my WordPress-related plans were on hold. To pass the time, I built a simple water simulation in Processing. Today I’m going to show you this little application and explain how it works. Online demo and source code are included.

### Big Words, or The Theoretical Part

Fluid dynamics is a complex topic. If you so much as look it up in Wikipedia, you will be immediately assaulted by integrals, differentials and other agents of mathematical insanity. To accurately model a fluid you would need a quite an in-depth understanding of physics and calculus. And while that’s certainly not something one couldn’t figure out given enough time, “proper” fluid simulation is still extremely rare in non-academic applications and games because it’s **slow**.

One way to create a faster (though less accurate) simulation is by using cellular automata to represent the water. A cellular automaton is basically a grid where each cell can be in one of a finite number of states, plus a set of rules that determine how a cell can change from one state to another. The rules are typically local – that is, they consider only the current cell and it’s direct neighbours when determining the new state. During each step of the simulation, you simply loop through the entire grid and apply the rules to each cell.

When simulating water or other fluids it can more appropriate to store a continuous value (e.g. water mass) in each cell, instead of using discrete states.

### CA-Based Water Simulation

Before we go into details, lets take a look at the promised online demo. It will give you a good idea about what you can achieve with this approach :

*(klicken Sie auf das Bild; Java required)*

Alright, on to the code stuff. Lets try simulating the most common fluid – water.

**The basic data structure**

//Map dimensions final int map_width = 16; final int map_height = 16; //Block types final int AIR = 0; final int GROUND = 1; final int WATER = 2; //Data structures int[][] blocks = new int[map_width+2][map_height+2]; float[][] mass = new float[map_width+2][map_height+2], new_mass = new float[map_width+2][map_height+2];

I use two two-dimensional arrays to represent the simulation world. The `blocks`

array defines a basic “map” where each cell can contain ground (solid, stops water flow), air (an empty cell, water can flow in) or water. Each water cell also has a corresponding entry in the `mass`

array that defines how much water it contains. All water cells/blocks start out with one unit of water.

There’s also a third array – `new_mass`

that is used to store intermediate mass values when running the simulation. Operating directly on the `mass`

array is unadvisable because you would get various sideffects like water spreading at different speeds depending on the order in which you update the cells.

**Overall algorithm**

The main idea is to treat water as a slightly compressible liquid, as described in this article. In terms of implementation, this means that if there are two or more water cells stacked vertically, the bottom cells will be able to hold slightly more water than normal. This way we don’t need to explicitly track pressure to make the water level equalize in communicating vessels – we can just look at how much excess water a cell has, and move it upwards if required.

Lets set up the exact fluid properties as constants :

//Water properties final float MaxMass = 1.0; //The normal, un-pressurized mass of a full water cell final float MaxCompress = 0.02; //How much excess water a cell can store, compared to the cell above it final float MinMass = 0.0001; //Ignore cells that are almost dry

For example, if one cell contains 1.0 units of water, the cell below it should contain up to 1.02 units, the cell below that one should contain 1.04, the next one 1.06, and so on. There is one however special case that we need to consider to get a plausible simulation – what if the if the top cell contains less than `MaxMass`

units of water? In my implementation the bottom cell will contain a proportionally smaller excess amount, not the usual `mass_of_the_cell_above + MaxCompression`

units.

Based on the above, we can simulate water movement by looping over the `mass`

array and applying these cellular automaton rules to each cell :

- Take the mass of the current cell and the cell
**below**it and figure out how much water the bottom cell should contain. If it has less than that, remove the corresponding amount from the current cell and add it to the bottom cell. - Check the cell to the
**left**of this one. If it has less water, move over enough water to make both cells contain the same amount. - Do the same thing for the
**right**neighbour. - Do the same thing as in step 1., but for the cell
**above**the current one.

While doing this you also need to keep track of how much water the current cell still has remaining, or you could end up with negative-mass water blocks.

After the new mass of each cell has been calculated, you can safely copy over the new values to the `mass`

array. You also need to update the `blocks`

array so that cells that are now dry get marked as empty and vice versa.

Here’s the simulation code :

void simulate_compression(){ float Flow = 0; float remaining_mass; //Calculate and apply flow for each block for (int x = 1; x <= map_width; x++){ for(int y = 1; y <= map_height; y++){ //Skip inert ground blocks if ( blocks[x][y] == GROUND) continue; //Custom push-only flow Flow = 0; remaining_mass = mass[x][y]; if ( remaining_mass <= 0 ) continue; //The block below this one if ( (blocks[x][y-1] != GROUND) ){ Flow = get_stable_state_b( remaining_mass + mass[x][y-1] ) - mass[x][y-1]; if ( Flow > MinFlow ){ Flow *= 0.5; //leads to smoother flow } Flow = constrain( Flow, 0, min(MaxSpeed, remaining_mass) ); new_mass[x][y] -= Flow; new_mass[x][y-1] += Flow; remaining_mass -= Flow; } if ( remaining_mass <= 0 ) continue; //Left if ( blocks[x-1][y] != GROUND ){ //Equalize the amount of water in this block and it's neighbour Flow = (mass[x][y] - mass[x-1][y])/4; if ( Flow > MinFlow ){ Flow *= 0.5; } Flow = constrain(Flow, 0, remaining_mass); new_mass[x][y] -= Flow; new_mass[x-1][y] += Flow; remaining_mass -= Flow; } if ( remaining_mass <= 0 ) continue; //Right if ( blocks[x+1][y] != GROUND ){ //Equalize the amount of water in this block and it's neighbour Flow = (mass[x][y] - mass[x+1][y])/4; if ( Flow > MinFlow ){ Flow *= 0.5; } Flow = constrain(Flow, 0, remaining_mass); new_mass[x][y] -= Flow; new_mass[x+1][y] += Flow; remaining_mass -= Flow; } if ( remaining_mass <= 0 ) continue; //Up. Only compressed water flows upwards. if ( blocks[x][y+1] != GROUND ){ Flow = remaining_mass - get_stable_state_b( remaining_mass + mass[x][y+1] ); if ( Flow > MinFlow ){ Flow *= 0.5; } Flow = constrain( Flow, 0, min(MaxSpeed, remaining_mass) ); new_mass[x][y] -= Flow; new_mass[x][y+1] += Flow; remaining_mass -= Flow; } } } //Copy the new mass values to the mass array for (int x = 0; x < map_width + 2; x++){ for (int y = 0; y < map_height + 2; y++){ mass[x][y] = new_mass[x][y]; } } for (int x = 1; x <= map_width; x++){ for(int y = 1; y <= map_height; y++){ //Skip ground blocks if(blocks[x][y] == GROUND) continue; //Flag/unflag water blocks if (mass[x][y] > MinMass){ blocks[x][y] = WATER; } else { blocks[x][y] = AIR; } } } //Remove any water that has left the map for (int x =0; x < map_width+2; x++){ mass[x][0] = 0; mass[x][map_height+1] = 0; } for (int y = 1; y < map_height+1; y++){ mass[0][y] = 0; mass[map_width+1][y] = 0; } }

And here’s the function that calculates how to distribute a given amount of water between two vertically adjacent cells :

//Returns the amount of water that should be in the bottom cell. float get_stable_state_b ( float total_mass ){ if ( total_mass <= 1 ){ return 1; } else if ( total_mass < 2*MaxMass + MaxCompress ){ return (MaxMass*MaxMass + total_mass*MaxCompress)/(MaxMass + MaxCompress); } else { return (total_mass + MaxCompress)/2; } }

**Full source code**

**Download the source**(23 KB)- You will need the Java-based Processing programming language to run the code.

I didn’t discuss how to render the water, or how to handle user interaction, but the source code also includes a simplistic implementation of those things.

### Conclusion

CA-based algorithms have both advantages and disadvantages. On the one hand, they’re relatively fast and the simulation quality is often good enough for games and software toys. On the other, they can produce behaviour that is very unrealistic – for example, water falling into a basin will form a little “hill” and spread out slowly, instead of near-instantly as real water would. Whether this is an acceptable downside depends on your application.

For a different take on simulating fluid dynamics with cellular automata, check out this interview with the author of Dwarf Fortress.

**Related posts :**

[…] Source code and a live demo are also available. This tutorial is along the same vein as my previous fluid simulation post that dealt with pseudo-compressible fluids and pressure, but the specific algorithm is much simpler […]

Thank you for this instructive tutorial, and final link to Dwarf fortress’ system. It’s the only simple tutorial I found, and it’s pretty neat. I didn’t actually run your code (live demo triggers errors with Chrome, and I didn’t bother installing Processing), but it was very instructive, thank you !

In water_sim.pde lign 440 you used map_height instead of map_width.

Seems good, but where do you define “constrain”?

constrain is pre-defined library function in the Processing environment.

As for the riche issues, louis vuitton handbags’s ready to wear is generally hidden away in the upper floors of its stores. The louis vuitton handbags brand obtained the bronze medal at the 1867 World’s Fair and a gold

grill worn on the bottom of the bag on the market.

So congratulations to louis vuitton handbags designer Marc Jacobs backstage after

his louis vuitton handbags show on Wednesday, July 13, 2011 at 2223 North West Shore Boulevard

in Tampa at 6:39 pm.

Could you explain the branches of get_stable_state, I don’t really understand that part ?

Other than that, thanks this is exactly what I’ve been looking for.

It’s been a long time since I wrote that. Lets see.. here’s what I have in my notes:

I’m trying to convert this script to 3d, which is not hard, but i also want bytes,not floats, so this is as far as i got.

void simulate_compression(){

int Flow = 0;

int MinFlow = 1;

int MaxSpeed=8;

int remaining_mass;

int[,,] mass = new int[16,16,16];;

int[,,] new_mass= new int[16,16,16];

//Calculate and apply flow for each block

for (int x = 0; x < 16; x++) {

for (int y = 0; y < 16; y++) {

for (int z = 0; z < 16; z++) {

mass [x, y, z] = thischunk.Fluids [x, y, z].water;

new_mass [x, y, z] = 0;

}

}

}

for (int x = 0; x < 16; x++) {

for (int y = 0; y < 16; y++) {

for (int z = 0; z < 16; z++) {

//Skip inert ground blocks

if (thischunk.blocks [x, y, z].GetBlockInfo ().watercapacity == 0) {

continue;

}

//Custom push-only flow

Flow = 0;

remaining_mass = mass [x,y,z];

if (remaining_mass MinFlow) {

// Flow = Mathf.FloorToInt(Flow* 0.5f); //leads to smoother flow

}

Flow = Mathf.Clamp (Flow, 0, Mathf.Min (MaxSpeed, remaining_mass));

new_mass [x, y, z] -= Flow;

new_mass [x, y – 1, z] += Flow;

remaining_mass -= Flow;

}

}

if (remaining_mass MinFlow) {

Flow = Mathf.RoundToInt(Flow* 0.5f); //leads to smoother flow

}

//Flow = Mathf.Clamp (Flow, 0, remaining_mass);

new_mass [x,y,z] -= Flow;

new_mass [x – 1,y,z] += Flow;

remaining_mass -= Flow;

}

}

if (remaining_mass MinFlow) {

Flow = Mathf.RoundToInt(Flow* 0.5f); //leads to smoother flow

}

//Flow = Mathf.Clamp (Flow, 0, remaining_mass);

new_mass [x,y,z] -= Flow;

new_mass [x,y,z – 1] += Flow;

remaining_mass -= Flow;

}

}

if (remaining_mass MinFlow) {

Flow = Mathf.RoundToInt(Flow* 0.5f); //leads to smoother flow

}

//Flow = Mathf.Clamp (Flow, 0, remaining_mass);

new_mass [x,y,z] -= Flow;

new_mass [x,y,z + 1] += Flow;

remaining_mass -= Flow;

}

}

if (remaining_mass MinFlow) {

Flow = Mathf.RoundToInt(Flow* 0.5f); //leads to smoother flow

}

//Flow = Mathf.Clamp (Flow, 0, remaining_mass);

new_mass [x,y,z] -= Flow;

new_mass [x + 1,y,z] += Flow;

remaining_mass -= Flow;

}

}

if (remaining_mass MinFlow) {

// Flow = Mathf.FloorToInt(Flow* 0.5f); //leads to smoother flow

}

Flow = Mathf.Clamp (Flow, 0, Mathf.Min (MaxSpeed, remaining_mass));

new_mass [x,y,z] -= Flow;

new_mass [x,y + 1,z] += Flow;

remaining_mass -= Flow;

}

}

}

}

}

//Copy the new mass values to the mass array

for (int x = 0; x < 16; x++) {

for (int y = 0; y < 16; y++) {

for (int z = 0; z < 16; z++) {

if (new_mass [x,y,z]<=0) {

thischunk.Fluids [x, y, z].water = 0;

continue;

}

thischunk.Fluids [x, y, z].water = (byte)new_mass [x,y,z];

}

}

}

}

//Returns the amount of water that should be in the bottom cell.

private int get_stable_state_b ( int total_mass ){

if ( total_mass <= 8 ){

return 8;

} else if ( total_mass < 2*MaxMass + MaxCompress ){

return Mathf.RoundToInt((MaxMass*MaxMass + total_mass*MaxCompress)/(MaxMass + MaxCompress));

} else {

return Mathf.RoundToInt(total_mass + MaxCompress)/2;

}

}

I can't get the water to level out because of the water next to it is only one less, and not less then one. anyideas how i could fix this problem.

a block holds 8 water normal and 255max(cause bytes)

Please help.

I haven’t done any fluid simulation in year. I’m afraid I won’t be able to help much.

That said, if I was trying to adapt this algorithm to bytes, I would probably set the normal amount of water for a “full” block to 128 instead of 8. Since bytes have a much smaller number of possible values than floats, moving the baseline closer to the middle of the range gives you more flexibility and makes precision errors less visible.

I would also ignore differences of one or less because bytes don’t have enough precision to handle them.

Maybe something like that would work?