Sign in to follow this  
bassman

extending stable fluids to 3d

Recommended Posts

bassman    122
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 );
	}
}

from which I made:
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 this post


Link to post
Share on other sites
bassman    122
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 this post


Link to post
Share on other sites
Eelco    301
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 this post


Link to post
Share on other sites
bassman    122
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.

thanks for your reaction,
bassman

Share this post


Link to post
Share on other sites

Create an account or sign in to comment

You need to be a member in order to leave a comment

Create an account

Sign up for a new account in our community. It's easy!

Register a new account

Sign in

Already have an account? Sign in here.

Sign In Now

Sign in to follow this