Jump to content
  • Advertisement
Sign in to follow this  
Gon

Multigrid for Jos Stam Stable fluids help needed!

This topic is 2532 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

To anyone reading this ahoy there,

Now back to serious mode, I'm working in my old stable fluids code and I was trying to implement the Multigrid method (not full MG)
I wanted to implement the V-cycle and latter the W-cycle.
However after some searching on the internet (a lot of searching) I managed to program something which is not giving the expected results (code will follow briefly).
I know this because I have Jacobi, Gauss-Seidel, SOR, and simple Conjugate Gradient working all giving the same results yet MG V-cycle is not.
So anyone familiar with this that can find any issue in the code some help would be most welcomed. I'm just testing for a simple 8x8x8 matrix. Also I perform 13 iters each smooth (Using Jacobi solver). And Jacobi alone gets there for the test case in less than the 13.


#define _LEVELS 2

float CLINEAR_SOLVERS::laplacian(long i,long j,long k,float *x,float a,float iter)
{
return((x[_IX(i,j,k)]*iter)-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)]));
}

void CLINEAR_SOLVERS::_mg(CGRID *g,float *x,float *x0,float a,float iter,long iters)
{
long lvl;
float *r=(float*) calloc(g->T_cells, sizeof(float));

if(r!=NULL) {
for(lvl=0;lvl<(_LEVELS-1);lvl++){
smoth(g,x,x0,a,iter,iters,lvl*2);
update_r(g,r,x,x0,a,iter,lvl*2);
restrict(g,r,x,x0,lvl*2);
}

smoth(g,x,x0,a,iter,iters,(_LEVELS-1)*2); //^-1 ??????

for(lvl=(_LEVELS-2);lvl>=0;lvl--){
prolong(g,x,lvl*2);
smoth(g,x,x0,a,iter,iters,lvl*2);
}

free(r);
}
}

void CLINEAR_SOLVERS::smoth(CGRID *g,float *x,float *x0,float a,float iter, long iters,long inc)
{
long i,j,k,it;
float *aux=(float*) calloc(g->T_cells, sizeof(float));

if(inc==0) inc++;

if(aux) {
for(it=0;it<iters;it++) {
FOR_EACH_CELL1
aux[_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)]))/iter;
END_FOR_EACH_CELL

FOR_EACH_CELL1
#ifdef _SOR
x[_IX(i,j,k)]=(1-_W)*x[_IX(i,j,k)]+_W*aux[_IX(i,j,k)];
#else
x[_IX(i,j,k)]=aux[_IX(i,j,k)];
#endif
END_FOR_EACH_CELL
}

free(aux);
}
}

void CLINEAR_SOLVERS::update_r(CGRID *g,float *r,float *x,float *x0,float a,float iter,long inc)
{
long i,j,k;

if(inc==0) inc++;

FOR_EACH_CELL1
r[_IX(i,j,k)]=x0[_IX(i,j,k)]-laplacian(i,j,k,x,a,iter);
END_FOR_EACH_CELL
}

void CLINEAR_SOLVERS::restrict(CGRID *g,float *r,float *x,float *x0,long inc)
{
long i,j,k;

if(inc==0) inc++;

// 64 points stencil full weighted
// | |
// ||1 2 1||2 4 2||1 2 1||
//1/64||2 4 2||4 8 4||2 4 2||
// ||1 2 1||2 4 2||1 2 1||
// | |

// 32 points stencil half weighted
// | |
// ||0 1 0||1 2 1||0 1 0||
//1/32||1 2 1||2 8 2||1 2 1||
// ||0 1 0||1 2 1||0 1 0||
// | |

// 16 points stencil
// | |
// ||0 0 0||0 1 0||0 0 0||
//1/16||0 2 0||1 8 1||0 2 0||
// ||0 0 0||0 1 0||0 0 0||
// | |

FOR_EACH_CELL1
x0[_IX(i,j,k)]=((0.5*(
r[_IX(i-1,j,k)]+r[_IX(i+1,j,k)]+
r[_IX(i,j-1,k)]+r[_IX(i,j+1,k)]+
r[_IX(i,j,k-1)]+r[_IX(i,j,k+1)]))+(0.25*(
r[_IX(i-1,j,k-1)]+r[_IX(i+1,j,k-1)]+
r[_IX(i,j-1,k-1)]+r[_IX(i,j+1,k-1)]+
r[_IX(i-1,j,k+1)]+r[_IX(i+1,j,k+1)]+
r[_IX(i,j-1,k+1)]+r[_IX(i,j+1,k+1)]+
r[_IX(i-1,j-1,k)]+r[_IX(i-1,j+1,k)]+
r[_IX(i+1,j-1,k)]+r[_IX(i+1,j+1,k)]))+(0.125*(
r[_IX(i-1,j-1,k-1)]+r[_IX(i-1,j+1,k-1)]+
r[_IX(i-1,j-1,k+1)]+r[_IX(i-1,j+1,k+1)]+
r[_IX(i+1,j-1,k-1)]+r[_IX(i+1,j+1,k-1)]+
r[_IX(i+1,j-1,k+1)]+r[_IX(i+1,j+1,k+1)]))+
r[_IX(i,j,k)]);


//performing 64 point stencil
//x0[_IX(i,j,k)]=1/64*((4*(
// r[_IX(i-1,j,k)]+r[_IX(i+1,j,k)]+
// r[_IX(i,j-1,k)]+r[_IX(i,j+1,k)]+
// r[_IX(i,j,k-1)]+r[_IX(i,j,k+1)]))+(2*(
// r[_IX(i-1,j,k-1)]+r[_IX(i+1,j,k-1)]+
// r[_IX(i,j-1,k-1)]+r[_IX(i,j+1,k-1)]+
// r[_IX(i-1,j,k+1)]+r[_IX(i+1,j,k+1)]+
// r[_IX(i,j-1,k+1)]+r[_IX(i,j+1,k+1)]+
// r[_IX(i-1,j-1,k)]+r[_IX(i-1,j+1,k)]+
// r[_IX(i+1,j-1,k)]+r[_IX(i+1,j+1,k)]))+
// r[_IX(i-1,j-1,k-1)]+r[_IX(i-1,j+1,k-1)]+
// r[_IX(i-1,j-1,k+1)]+r[_IX(i-1,j+1,k+1)]+
// r[_IX(i+1,j-1,k-1)]+r[_IX(i+1,j+1,k-1)]+
// r[_IX(i+1,j-1,k+1)]+r[_IX(i+1,j+1,k+1)]+
// (8*(r[_IX(i,j,k)])));

//performing 32 point stencil
//x0[_IX(i,j,k)]=1/32*((2*(
// r[_IX(i-1,j,k)]+r[_IX(i+1,j,k)]+
// r[_IX(i,j-1,k)]+r[_IX(i,j+1,k)]+
// r[_IX(i,j,k-1)]+r[_IX(i,j,k+1)]))+
// r[_IX(i-1,j,k-1)]+r[_IX(i+1,j,k-1)]+
// r[_IX(i,j-1,k-1)]+r[_IX(i,j+1,k-1)]+
// r[_IX(i-1,j,k+1)]+r[_IX(i+1,j,k+1)]+
// r[_IX(i,j-1,k+1)]+r[_IX(i,j+1,k+1)]+
// r[_IX(i-1,j-1,k)]+r[_IX(i-1,j+1,k)]+
// r[_IX(i+1,j-1,k)]+r[_IX(i+1,j+1,k)]+
// (8*(r[_IX(i,j,k)])));

//performing 16 point stencil
//x0[_IX(i,j,k)]=1/16*(
// r[_IX(i-1,j,k)]+r[_IX(i+1,j,k)]+
// r[_IX(i,j-1,k)]+r[_IX(i,j+1,k)]+
// (2*(r[_IX(i,j,k-1)]+r[_IX(i,j,k+1)]))+8*(r[_IX(i,j,k)]));
x[_IX(i,j,k)]=0;
END_FOR_EACH_CELL
}

void CLINEAR_SOLVERS::prolong(CGRID *g,float *x,long inc)
{
long i,j,k;

if(inc==0) inc++;

FOR_EACH_CELL1
x[_IX(i,j,k)]+=x[_IX(i,j,k)];
x[_IX(i-1,j,k)]+=(0.5f*x[_IX(i,j,k)]);
x[_IX(i+1,j,k)]+=(0.5f*x[_IX(i,j,k)]);
x[_IX(i,j-1,k)]+=(0.5f*x[_IX(i,j,k)]);
x[_IX(i,j+1,k)]+=(0.5f*x[_IX(i,j,k)]);
x[_IX(i,j,k-1)]+=(0.5f*x[_IX(i,j,k)]);
x[_IX(i,j,k+1)]+=(0.5f*x[_IX(i,j,k)]);

if((k!=1)&&(k!=((g->NZ)-2))){
x[_IX(i-1,j,k-1)]+=(0.25f*x[_IX(i,j,k)]);
x[_IX(i+1,j,k-1)]+=(0.25f*x[_IX(i,j,k)]);
x[_IX(i,j-1,k-1)]+=(0.25f*x[_IX(i,j,k)]);
x[_IX(i,j+1,k-1)]+=(0.25f*x[_IX(i,j,k)]);

x[_IX(i-1,j,k+1)]+=(0.25f*x[_IX(i,j,k)]);
x[_IX(i+1,j,k+1)]+=(0.25f*x[_IX(i,j,k)]);
x[_IX(i,j-1,k+1)]+=(0.25f*x[_IX(i,j,k)]);
x[_IX(i,j+1,k+1)]+=(0.25f*x[_IX(i,j,k)]);

x[_IX(i-1,j-1,k-1)]+=(0.125f*x[_IX(i,j,k)]);
x[_IX(i-1,j+1,k-1)]+=(0.125f*x[_IX(i,j,k)]);
x[_IX(i+1,j-1,k-1)]+=(0.125f*x[_IX(i,j,k)]);
x[_IX(i+1,j+1,k-1)]+=(0.125f*x[_IX(i,j,k)]);

x[_IX(i-1,j-1,k+1)]+=(0.125f*x[_IX(i,j,k)]);
x[_IX(i-1,j+1,k+1)]+=(0.125f*x[_IX(i,j,k)]);
x[_IX(i+1,j-1,k+1)]+=(0.125f*x[_IX(i,j,k)]);
x[_IX(i+1,j+1,k+1)]+=(0.125f*x[_IX(i,j,k)]);
}
END_FOR_EACH_CELL
}

Share this post


Link to post
Share on other sites
Advertisement
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!