Jump to content
  • Advertisement
Sign in to follow this  
bassman

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.

If you intended to correct an error in the post then please contact us.

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

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
Advertisement
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
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
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
Sign in to follow this  

  • Advertisement
×

Important Information

By using GameDev.net, you agree to our community Guidelines, Terms of Use, and Privacy Policy.

We are the game development community.

Whether you are an indie, hobbyist, AAA developer, or just trying to learn, GameDev.net is the place for you to learn, share, and connect with the games industry. Learn more About Us or sign up!

Sign me up!