W-Shadow.com

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] = 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.