Sign in to follow this  
Winegums

Jacobi Iterative Fluid Solver

Recommended Posts

Hi, I'm having some trouble writing an iterative fluid solver. Based on the Jos Stam stuff (http://www.dgp.toronto.edu/people/stam/reality/Research/pdf/GDC03.pdf) I identified the point in the Projection routine where he uses a gauss-seidel solver, and I've tried putting in a jacobi one instead, but it doesn't work. I've looked at two different sources which explain a code based implementation using shaders (which should transfer over into a C++ implementation easily enough), but no dice. the method is as follows:
void JacobiSolver(float* x, float* b)
{

	float *x0 = new float[size];

	float alpha = -1;

	for(int i = 0; i < size ; i++)
	{	
		x0[i] = x[i];
	}

	for(int k = 0; k < 20 ; k++)		//Rather than converging, just repeat
	{
		for (int i=1 ; i<=N ; i++ )	//for x
		{
			for (int j=1 ; j<=N ; j++ )	//for y
			{
				x[IX(i, j)] = (x0[IX(i-1,j)] + x0[IX(i+1,j)] + x0[IX(i,j-1)] + x0[IX(i,j+1)] + (alpha * b[IX(i,j)]))/4;
			}
		}

		for(int i = 0; i < size ; i++)
		{
			x0[i] = x[i];
		}	
		
	}

	delete x0;
}

I noticed that in this paper alpha is defined as one grid step squared, whereas here it is defined as -1. I'm not sure if this is an implementation difference, but either way it doesn't work. The other possibility I've considered is that how I'm using the data is wrong. Does anything change when swapping from gauss-Seidel to Jacobi as far as the rest of the helmholtz-hodge approach is concerned?

Share this post


Link to post
Share on other sites
The only difference between Jacobi and Gauss-Seidel is that in Gauss-Seidel you replace x0 with x, and eliminate the step where you copy the x values into x0. This usually just speeds up the process because you're always using the latest values of x and you eliminate the updating step at the end of every iteration.

Alpha is part of the system you're solving, not part of the algorithm...

And a GPU based solution might use Jacobi's method instead of Gauss-Seidel because it may be better suited to a GPU, where all the calculations for an iteration are done in parallel. If you're doing sequential updates on each element, just use Gauss-Seidel.

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