## 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 :

1. 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.
2. 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.
3. Do the same thing for the right neighbour.
4. 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;
mass[x][map_height+1] = 0;
}
for (int y = 1; y < map_height+1; y++){
mass[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

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 :

### 9 Responses to “Simple Fluid Simulation With Cellular Automata”

1. […] 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 […]

2. Corentin says:

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.

3. Dave says:

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

4. Jānis Elsts says:

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

5. Louis Vuirron Outlet says:

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.

6. John says:

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.

7. Jānis Elsts says:

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

```if x <= 1, all water goes to the lower cell
* a = 0
* b = 1

if x > 1 & x < 2*MaxMass + MaxCompress, the lower cell should have MaxMass + (upper_cell/MaxMass) * MaxCompress
b = MaxMass + (a/MaxMass)*MaxCompress
a = x - b

->

b = MaxMass + ((x - b)/MaxMass)*MaxCompress ->
b = MaxMass + (x*MaxCompress - b*MaxCompress)/MaxMass
b*MaxMass = MaxMass^2 + (x*MaxCompress - b*MaxCompress)
b*(MaxMass + MaxCompress) = MaxMass*MaxMass + x*MaxCompress

* b = (MaxMass*MaxMass + x*MaxCompress)/(MaxMass + MaxCompress)
* a = x - b;

if x >= 2 * MaxMass + MaxCompress, the lower cell should have upper+MaxCompress

b = a + MaxCompress
a = x - b

->

b = x - b + MaxCompress ->
2b = x + MaxCompress ->

* b = (x + MaxCompress)/2
* a = x - b
```
8. Ronnie says:

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)
9. Jānis Elsts says: