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 :

Click to run the simulation (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] = 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

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 :

11 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. 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)
    Please help.

  9. Jānis Elsts says:

    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?

Leave a Reply