# extending stable fluids to 3d

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

## Recommended Posts

hi everyone I've been trying to extend Stam's 2d solver I found on the web to 3 dimensions. I'm doing this "on sight" as I am not familiar with cfd, but I do feel I have a general understanding of the principles explained in "Stable fluids" and "real time fluid dynamics for games". My problem is the 3d solver is not stable thus exploding from time to time. The 2d version works fine though. The obvious place to look for a bug is the projection step. I did some research, annoyed some people with experience in cfd, but no one could tell me why it explodes, so I figured its time to bother the gamedev forums :) stam projects his vector fields like this: (using a cell-centered grid)

#define IX(i,j) ((i)+(N+2)*(j))

void project ( int N, float * u, float * v, float * p, float * div )
{
int i, j;

FOR_EACH_CELL
div[IX(i,j)] = -0.5f*(u[IX(i+1,j)]-u[IX(i-1,j)]
+v[IX(i,j+1)]-v[IX(i,j-1)])/N;
p[IX(i,j)] = 0;
END_FOR
set_bnd ( N, 0, div ); set_bnd ( N, 0, p );

lin_solve ( N, 0, p, div, 1, 4 );

FOR_EACH_CELL
u[IX(i,j)] -= 0.5f*N*(p[IX(i+1,j)]-p[IX(i-1,j)]);
v[IX(i,j)] -= 0.5f*N*(p[IX(i,j+1)]-p[IX(i,j-1)]);
END_FOR
set_bnd ( N, 1, u ); set_bnd ( N, 2, v );
}


I extended this to:
#define IX(i,j,k) ((i) + (N+2)*(j) + ((N*N)+(4*N)+(4))*(k))

void project_3d ( int N, float * u, float * v, float * w, float * p, float * div )
{
int i, j, k;

FOR_EACH_CELL
div[IX(i,j,k)] = - 0.5f*(u[IX(i+1,j,k)]-u[IX(i-1,j,k)]
+v[IX(i,j+1,k)]-v[IX(i,j-1,k)]
+w[IX(i,j,k+1)]-w[IX(i,j,k-1)])/N;
p[IX(i,j,k)] = 0;
END_FOR
set_bnd_3d( N, 0, div ); set_bnd_3d ( N, 0, p );

lin_solve_3d ( N, 0, p, div, 1, 6 );

FOR_EACH_CELL

u[IX(i,j,k)] -= 0.5f*N*(p[IX(i+1,j,k)]-p[IX(i-1,j,k)]);
v[IX(i,j,k)] -= 0.5f*N*(p[IX(i,j+1,k)]-p[IX(i,j-1,k)]);
w[IX(i,j,k)] -= 0.5f*N*(p[IX(i,j,k+1)]-p[IX(i,j,k-1)]);

END_FOR
set_bnd_3d ( N, 1, u ); set_bnd_3d ( N, 2, v ); set_bnd_3d ( N, 3, w );

}


The lin_solve function uses a gauss-seidel relaxation method to solve a linear system. I considered implementing my own conjugate gradient method, but i'd rather get this method to work first. After all, gauss seidel should work, regardless of the number of dimensions right? anyways, here's the code:
void lin_solve ( int N, int b, float * x, float * x0, float a, float c )
{
int i, j, k;

for ( k=0 ; k<20 ; k++ ) {
FOR_EACH_CELL
x[IX(i,j)] = (x0[IX(i,j)] +
a*( x[IX(i-1,j)]+x[IX(i+1,j)]
+x[IX(i,j-1)]+x[IX(i,j+1)]))/c;
END_FOR
set_bnd ( N, b, x );
}
}


void lin_solve_3d ( int N, int b, float * x, float * x0, float a, float c )
{
int i, j, k, count;

int sz = (N+2)*(N+2)*(N+2);
float * temp = (float*) malloc (sz*sizeof(float) );

for ( count=0 ; count<100 ; count++)
{
FOR_EACH_CELL
temp[IX(i,j,k)] = ( x0[IX(i,j,k)]
+ a*(x[IX(i-1,j,k)]+x[IX(i+1,j,k)]
+x[IX(i,j-1,k)]+x[IX(i,j+1,k)]
+x[IX(i,j,k-1)]+x[IX(i,j,k+1)])
) / c;
END_FOR

FOR_EACH_CELL
x[IX(i,j,k)] = temp[IX(i,j,k)];
END_FOR

set_bnd_3d ( N, b, x );
}

free( temp );
}

The boundaries used are "walls", maybe there's a mistake there?

anyways, thanks to everyone getting to read this post this far,
many more thanks to people who might be able to help :)
I really look forward to be able to play with a 3d solver.

Greets ,
bassman flash



##### Share on other sites
sorry, forget to [/source] the last chunck 'o code ^^,
so here's the outro:

The boundaries used are "walls", maybe there's a mistake there?

anyways, thanks to everyone getting to read this post this far,
many more thanks to people who might be able to help :)
I really look forward to be able to play with a 3d solver.

Greets ,
bassman flash

##### Share on other sites
i also implemented that paper, one of the best ive ever read.

in any case, your 3d gauss-seidel and 2d version are not identical. why the use of a temp array in the 3d case? not sure if thats the cause of instability, but its not going to produce the precise same results, thats for sure.

##### Share on other sites
the introduction of temp was one of my many (near desperate)
throws to a solution - didn't help. I saw it in some guys implementation
of the 2d version. It seemed logical, cause the calculation is
more symmetrical this way. But the solver blows up with and withouth it.

bassman

1. 1
Rutin
46
2. 2
3. 3
4. 4
5. 5
JoeJ
18

• 13
• 10
• 12
• 10
• 13
• ### Forum Statistics

• Total Topics
632998
• Total Posts
3009802
• ### Who's Online (See full list)

There are no registered users currently online

×