# Multigrid for Jos Stam Stable fluids help needed!

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

## 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 } 

• ### What is your GameDev Story?

In 2019 we are celebrating 20 years of GameDev.net! Share your GameDev Story with us.

• 28
• 16
• 10
• 10
• 11
• ### Forum Statistics

• Total Topics
634111
• Total Posts
3015563
×