# Fluid Simulations - Projection/Mass Conservation

This topic is 3678 days old which is more than the 365 day threshold we allow for new replies. Please post a new topic.

## Recommended Posts

Hi, I've been reading and implimenting a method presented in this paper: http://www.dgp.toronto.edu/people/stam/reality/Research/pdf/GDC03.pdf One of the key functions of this method is the Projection() function, which is presented below:

void project ( int N, float * u, float * v, float * p, float * div )
{
int i, j, k;
float h;
h = 1.0/N;
for ( i=1 ; i<=N ; i++ )
{
for ( j=1 ; j<=N ; j++ )
{
div[IX(i,j)] = -0.5*h*(u[IX(i+1,j)]-u[IX(i-1,j)]+
v[IX(i,j+1)]-v[IX(i,j-1)]);
p[IX(i,j)] = 0;
}
}
set_bnd ( N, 0, div ); set_bnd ( N, 0, p );
for ( k=0 ; k<20 ; k++ ) {
for ( i=1 ; i<=N ; i++ ) {
for ( j=1 ; j<=N ; j++ ) {
p[IX(i,j)] = (div[IX(i,j)]+p[IX(i-1,j)]+p[IX(i+1,j)]+
p[IX(i,j-1)]+p[IX(i,j+1)])/4;
}
}
set_bnd ( N, 0, p );
}
for ( i=1 ; i<=N ; i++ ) {
for ( j=1 ; j<=N ; j++ ) {
u[IX(i,j)] -= 0.5*(p[IX(i+1,j)]-p[IX(i-1,j)])/h;
v[IX(i,j)] -= 0.5*(p[IX(i,j+1)]-p[IX(i,j-1)])/h;
}
}
set_bnd ( N, 1, u ); set_bnd ( N, 2, v );
}


In particular the first nested for loop confuses me.
for ( i=1 ; i<=N ; i++ )
{
for ( j=1 ; j<=N ; j++ )
{
div[IX(i,j)] = -0.5*h*(u[IX(i+1,j)]-u[IX(i-1,j)]+
v[IX(i,j+1)]-v[IX(i,j-1)]);
p[IX(i,j)] = 0;
}
}


Where is the equation for this coming from? Stam makes reference to Hodge decomposition, Poisson equation and Gauss-Seidel relaxation but I can't join the dots up.

##### Share on other sites
I think the first loop is just setting up some initial conditions for the second loop, which solves the poisson equation. The div array contains the net velocity flowing into (or from) each of the cells, and p will contain the solution to the poisson equation.

Edit: I read some more about it, and it isn't that complicated. As described in the paper, the velocity flows can be decomposed into the mass conserving, swirly flows and the not swirly gradient flows. The not so swirly part is the part that causes the net flow into a cell to be non-zero. The first step, performed by the first loop, is to calculate the net flow into each cell. Then in the second loop, the poisson equation is solved through the iterative relaxation method. This results in p containing the "heightmap" as described in the paper. This heightmap's derivative is the flows that account for all the non-zero flows into the cells. The projection step is the subtraction of these flows from the original flows, which projects the original flows onto the space of mass-preserving flows, performed in the third loop. This results in nice swirly vortices.

Hodge decomposition = the flows can be decomposed into the mass-preserving (swirly) and gradient (not swirly) components
Poisson equation = the form of the differential equation that is solved as a step in finding the gradient component
Gauss-Seidel relaxation = the method used (in the second loop) to solve the poisson equation
Projection = Projecting the flows into the space of purely mass-preserving flows, accomplished by subtracting out the gradient flows. Just like projecting a vector v onto a plane with normal n with the formula (v - v*(v*n)), except that the part being subtracted out is found through the solution of a difference equation instead of a dot product.

[Edited by - Vorpy on January 28, 2008 1:24:46 AM]

##### Share on other sites
OK I think I get the theory side.

We enter this function knowing that our fluid is made up of two fields, an incompressible field and a gradient field. Gradient field is what we want to get rid of.

We use the first loop to find the net flow into each grid cell.
The second loop uses Poisson equation to find the gradient field.

The third loop then subtracts the gradient field from the fluid.

Yeah?

If so, this has led me to more questions...

i) Whats up with the first loop, still? Surley the net flow into the cell is..
(absoluteValue(grid[n+1][n] - grid[n-1][n]) +
absoluteValue(grid[n][n+1] - grid[n][n-1])) / 2

why is there multiplication by cell height and -0.5 (rather than just 0.5)?

ii) Is there a good place to find out about the poisson equation? The paper skims over it but I think I'll need to understand it reasonably well. I checked wikipedia, but I didn't really get it.

EDIT: Just noticed that at the start things are multiplied by h and at the end they are divided by h, thats something else I don't get..

Thanks a lot for replying Vorpy. If I could rate you twice I would...