# Cell Based Simulation

I'm working on a cell based simulation in a pixel shader (webgl, shadertoy.com specifically), but am having some challenges trying to make the simulation consistent.

Each cell (pixel) has the following format:

• Red Component = Cell Type
• Green, Blue, Alpha Component = Cell data for that specific cell type

Cell types include:

• empty space
• stone
• sand
• water
• lots of others

Logic rules are like:

• empty space does nothing
• stone does nothing
• if the cell below sand is empty, it will move there. Else it will try moving down left and down right.
• water has the same rules as sand, but if it wasn't able to make a move using that logic, it will try to flow to the left (simplification of the actual logic)
• more complex logic, like a seed might absorb neighboring water cells, and sprout a plant after absorbing 3 water cells.

A specific problem I'm trying to solve is that when two particles want to move to the same cell, they both do, but the first one is over-written by the second one and so is destroyed. Check the screenshot below: Black = empty space, white = stone, blue = water.

The problem is that the water on the left wants to move down right, and the water on the right wants to move down left. They both do that, and one of the water cells is effectively destroyed since they both cleared out where they were moving from.

The pixel shader program runs for each cell in the simulation, but of course only has write access to the cell that it's currently running for (since the shader can only emit one pixel - the pixel that it's being run for).

The program can do extra texture reads though to get knowledge of it's neighbors.

Right now I have the logic like this (simplified pseudo code):

if current cell is water
{
if cell below is empty
{
set cell below to current cell
set current cell to empty
}
else if cell to lower left is empty
{
set cell at lower left to current cell
set current cell to empty
}
else if cell to lower right is empty
{
set cell at lower right to current cell
set current cell to empty
}
}


To simplify logic, I have this code run for the current pixel, as well as all 8 neighboring pixels. This way, I don't have to write "source logic" and "destination logic", I can just write the natural logic, and both sides get picked up by virtue of neighbors being simulated, and the fact that there is only write access to the pixel we are actually running the program on.

I do make sure and give each simulation clean versions of the data (ie neighbor simulations changing data don't affect each other, they ALL are based only on the previous state of the cells last frame, and write the new state into a different buffer for the next frame).

To solve this, I thought about simulating not just the neighbors that are 1 cell away, but simulating the neighbors that are 2 cells away too, to be able to detect this case. In the 5x5 simulation code i would say "once a cell has been modified by any of the simulations in the 5x5 grid, don't allow it to be modified by other simulation operations in this same 5x5 grid".

I believe this would solve this immediate problem, but at the same time, since each cell's future state is now tied to neighbor cells up to 2 pixels away, that i can't just consider the 5x5 grid to know for sure what each cell in that grid will do, but would instead need to know the state of their neighbors that are up to 2 away, and the state of THEIR neighbors up to 2 away etc.

Does anyone know the correct way to deal with a simulation of this nature, to get decent / consistent results?

Edit: Here's the current source code of my simulation shader. It reads from last frame's state and outputs the current frame's state.

/*
============================================================
This buffer represents the "next state" of the simulation.

Buf A copies this data at the beginning of each frame.

Image renders this data.
============================================================
*/

// variables
const vec2 txState = vec2(0.0, 0.0);
// xyzw = dust cell value to spawn
// x = dust type.  yzw = dust params

// Grid Cell Format:
// x = dust type.  0..255 mapped to 0..1
// yzw = dust params

//============================================================

float isInside( vec2 p, vec2 c ) { vec2 d = abs(p-0.5-c) - 0.5; return -max(d.x,d.y); }
float isInside( vec2 p, vec4 c ) { vec2 d = abs(p-0.5-c.xy-c.zw*0.5) - 0.5*c.zw - 0.5; return -max(d.x,d.y); }

vec4 loadValue( in vec2 re )
{
return texture2D( iChannel0, (0.5+re) / iChannelResolution.xy, -100.0 );
}

void storeValue( in vec2 re, in vec4 va, inout vec4 fragColor, in vec2 fragCoord )
{
fragColor = ( isInside(fragCoord,re) > 0.0 ) ? va : fragColor;
}

void storeValue( in vec4 re, in vec4 va, inout vec4 fragColor, in vec2 fragCoord )
{
fragColor = ( isInside(fragCoord,re) > 0.0 ) ? va : fragColor;
}

//============================================================
vec4 GetCellData (in vec2 cell)
{
// Bottom row of pixels is special since it holds game state variables.
// Always return stone for those cells.
// Also return stone for out of bound reads.
if (cell.x >= 0.0 && cell.y >= 1.0 && cell.x < iChannelResolution.x && cell.y < iChannelResolution.y)
return texture2D( iChannel0, (cell + 0.5) / iChannelResolution.xy, -100.0 );
else
return vec4(1.0 / 255.0, 0.0, 0.0, 0.0);
}

//============================================================
void WriteCell (vec2 cell, vec4 cellData, in vec2 progCell, inout vec4 progCellData, inout bool progCellWritten)
{
if (progCellWritten || cell != progCell)
return;

progCellData = cellData;
progCellWritten = true;
}

//============================================================
float rand(vec2 co)
{
return fract(sin(dot(co.xy ,vec2(12.9898,78.233))) * 43758.5453);
}

//============================================================
void DoSim (
in vec2 cell,        // the cell coordinates we are running the sim for
in vec4 cellData02,
in vec4 cellData12,
in vec4 cellData22,
in vec4 cellData01,
in vec4 cellData11,  // the sim cell
in vec4 cellData21,
in vec4 cellData00,
in vec4 cellData10,
in vec4 cellData20,
in vec2 progCell,       // the cell that is running the simulation
inout vec4 progCellData,// the cell data for the cell running the simulation
inout bool progCellWritten
)
{
float dustType = floor(cellData11.x * 255.0);

// water logic
if (dustType == 2.0)
{
// if the tile below us is empty, move there
if (floor(cellData10.x*255.0) == 0.0)
{
WriteCell(cell + vec2( 0.0, -1.0), cellData11, progCell, progCellData, progCellWritten);
WriteCell(cell                   , vec4(0.0) , progCell, progCellData, progCellWritten);
}
// else if the tile to the bottom left is empty, move there
else if (floor(cellData00.x*255.0) == 0.0)
{
WriteCell(cell + vec2(-1.0, -1.0), cellData11, progCell, progCellData, progCellWritten);
WriteCell(cell                   , vec4(0.0) , progCell, progCellData, progCellWritten);
}
// else if the tile to the bottom right is empty, move there
else if (floor(cellData20.x*255.0) == 0.0)
{
WriteCell(cell + vec2( 1.0, -1.0), cellData11, progCell, progCellData, progCellWritten);
WriteCell(cell                   , vec4(0.0) , progCell, progCellData, progCellWritten);
}
// else if there's anything underneath us, do flow logic
else if (floor(cellData10.x*255.0) != 0.0)
{
// if flowing left
if (cellData11.y < 0.5)
{
// if the left is empty, move there
if (floor(cellData01.x*255.0) == 0.0)
{
WriteCell(cell + vec2(-1.0,  0.0), cellData11, progCell, progCellData, progCellWritten);
WriteCell(cell                   , vec4(0.0) , progCell, progCellData, progCellWritten);
}
// else flip flow direction
else
{
cellData11.y = 1.0;
WriteCell(cell + vec2( 0.0,  0.0), cellData11, progCell, progCellData, progCellWritten);
}
}
// else if flowing right
else
{
// if the right is empty, move there
if (floor(cellData21.x*255.0) == 0.0)
{
WriteCell(cell + vec2( 1.0,  0.0), cellData11, progCell, progCellData, progCellWritten);
WriteCell(cell                   , vec4(0.0) , progCell, progCellData, progCellWritten);
}
// else flip flow directions
else
{
cellData11.y = 0.0;
WriteCell(cell + vec2( 0.0,  0.0), cellData11, progCell, progCellData, progCellWritten);
}
}
}
}
// Sand Logic
else if (dustType == 3.0)
{
// if the tile below us is empty, move there
if (floor(cellData10.x*255.0) == 0.0)
{
WriteCell(cell + vec2( 0.0, -1.0), cellData11, progCell, progCellData, progCellWritten);
WriteCell(cell                   , vec4(0.0) , progCell, progCellData, progCellWritten);
}
// else if the tile to the bottom left is empty, move there
else if (floor(cellData00.x*255.0) == 0.0)
{
WriteCell(cell + vec2(-1.0, -1.0), cellData11, progCell, progCellData, progCellWritten);
WriteCell(cell                   , vec4(0.0) , progCell, progCellData, progCellWritten);
}
// else if the tile to the bottom right is empty, move there
else if (floor(cellData20.x*255.0) == 0.0)
{
WriteCell(cell + vec2( 1.0, -1.0), cellData11, progCell, progCellData, progCellWritten);
WriteCell(cell                   , vec4(0.0) , progCell, progCellData, progCellWritten);
}
}
}

//============================================================
void mainImage( out vec4 fragColor, in vec2 fragCoord )
{

// get the 5x5 grid of cells surrounding this cell.
// so we have the neighbors, and their neighbors.
// needed so we can sim ourselves and all neighbors.
// this ensures everyone is in agreement, and lets us write simple sim code.
vec2 cell       = floor(fragCoord);
vec4 cellData04 = GetCellData(cell + vec2(-2.0,  2.0));
vec4 cellData03 = GetCellData(cell + vec2(-2.0,  1.0));
vec4 cellData02 = GetCellData(cell + vec2(-2.0,  0.0));
vec4 cellData01 = GetCellData(cell + vec2(-2.0, -1.0));
vec4 cellData00 = GetCellData(cell + vec2(-2.0, -2.0));

vec4 cellData14 = GetCellData(cell + vec2(-1.0,  2.0));
vec4 cellData13 = GetCellData(cell + vec2(-1.0,  1.0));
vec4 cellData12 = GetCellData(cell + vec2(-1.0,  0.0));
vec4 cellData11 = GetCellData(cell + vec2(-1.0, -1.0));
vec4 cellData10 = GetCellData(cell + vec2(-1.0, -2.0));

vec4 cellData24 = GetCellData(cell + vec2( 0.0,  2.0));
vec4 cellData23 = GetCellData(cell + vec2( 0.0,  1.0));
vec4 cellData22 = GetCellData(cell + vec2( 0.0,  0.0));
vec4 cellData21 = GetCellData(cell + vec2( 0.0, -1.0));
vec4 cellData20 = GetCellData(cell + vec2( 0.0, -2.0));

vec4 cellData34 = GetCellData(cell + vec2( 1.0,  2.0));
vec4 cellData33 = GetCellData(cell + vec2( 1.0,  1.0));
vec4 cellData32 = GetCellData(cell + vec2( 1.0,  0.0));
vec4 cellData31 = GetCellData(cell + vec2( 1.0, -1.0));
vec4 cellData30 = GetCellData(cell + vec2( 1.0, -2.0));

vec4 cellData44 = GetCellData(cell + vec2( 2.0,  2.0));
vec4 cellData43 = GetCellData(cell + vec2( 2.0,  1.0));
vec4 cellData42 = GetCellData(cell + vec2( 2.0,  0.0));
vec4 cellData41 = GetCellData(cell + vec2( 2.0, -1.0));
vec4 cellData40 = GetCellData(cell + vec2( 2.0, -2.0));

// save off a copy, so that our mutable cell data isn't used in logic
vec4 progCellData = cellData22;

//----- Simulation -----
// Run simulation for this cell and all neighbors.
// This is to handle the fact that we or a neighbor might want to modify this particular cell.
// Doing it this way makes it so we write simpler, but more solid (less error prone) code.

//    Sim Cell Location       |               Top Row                |             Middle Row             |             Bottom Row              |  Program Cell Info

bool progCellWritten = false;
DoSim(cell + vec2(-1.0, -1.0),   cellData02, cellData12, cellData22,   cellData01, cellData11, cellData21,   cellData00, cellData10, cellData20,   cell, progCellData, progCellWritten);
DoSim(cell + vec2( 0.0, -1.0),   cellData12, cellData22, cellData32,   cellData11, cellData21, cellData31,   cellData10, cellData20, cellData30,   cell, progCellData, progCellWritten);
DoSim(cell + vec2( 1.0, -1.0),   cellData22, cellData32, cellData42,   cellData21, cellData31, cellData41,   cellData20, cellData30, cellData40,   cell, progCellData, progCellWritten);

DoSim(cell + vec2(-1.0,  0.0),   cellData03, cellData13, cellData23,   cellData02, cellData12, cellData22,   cellData01, cellData11, cellData21,   cell, progCellData, progCellWritten);
DoSim(cell + vec2( 0.0,  0.0),   cellData13, cellData23, cellData33,   cellData12, cellData22, cellData32,   cellData11, cellData21, cellData31,   cell, progCellData, progCellWritten);
DoSim(cell + vec2( 1.0,  0.0),   cellData23, cellData33, cellData43,   cellData22, cellData32, cellData42,   cellData21, cellData31, cellData41,   cell, progCellData, progCellWritten);

DoSim(cell + vec2(-1.0,  1.0),   cellData04, cellData14, cellData24,   cellData03, cellData13, cellData23,   cellData02, cellData12, cellData22,   cell, progCellData, progCellWritten);
DoSim(cell + vec2( 0.0,  1.0),   cellData14, cellData24, cellData34,   cellData13, cellData23, cellData33,   cellData12, cellData22, cellData32,   cell, progCellData, progCellWritten);
DoSim(cell + vec2( 1.0,  1.0),   cellData24, cellData34, cellData44,   cellData23, cellData33, cellData43,   cellData22, cellData32, cellData42,   cell, progCellData, progCellWritten);

//----- Save State -----
fragColor = vec4(0.0);
storeValue(txState, state   , fragColor, fragCoord);
if (cell.y > 0.0)
fragColor = progCellData;
}

• You are overwriting the neighbors. Then you are asking, "Why are my neighbors being overwritten?". Think about it. – MickLH Jan 22 '16 at 19:02
• Yeah that's the problem that I'm trying to solve. The issue is that my simulation is run for each cell, and has write access only to the current cell. I'm wondering about whether there's some techniques for simulating well within the restrictions I have imposed on me by the environment the simulation is running in. – Alan Wolfe Jan 22 '16 at 19:07
• Transform the logic so that instead of writing to a neighbor cell B while considering "center" cell A. You check cell A while considering cell B. – MickLH Jan 22 '16 at 19:09
• Please dont take this as disrespect, but your answers smell a bit of not understanding the problem sufficiently - or I'm just not getting what you are saying. Can you give an example of what you mean, and how it solves the specific problem I'm hitting, or any similar issues? – Alan Wolfe Jan 22 '16 at 19:21
• Sorry, but you're just not getting what I said at all. This problem is immensely simple, if I went any further, I would literally be writing the code for you. – MickLH Jan 22 '16 at 19:31