[arotester]

Hello,

I have implemented BiCGSTAB without any preconditioning and it works fine with the simple unsteady diffusion equation but when I use it to solve the pressure correction equation in NS(semi-Explicit formulation) equ. it shows the alpha and omega values in O(500) which is not correct and the solution diverges. I am trying a test case as a simple 2D channel flow.

I have already tested the flow through the channel with GSSOR solver and it is working well but as soon as I replace GSSOR with the BiCGSTAB the solution diverges.
I am not able to debug the issue, can anyone please give any insights? the code is as follows,

// 		-------------------------------------Pressure Correction calculations--------------------------------------
	for(i=1;i<=nx;i++)
		{
			pc[i][ny+1]=pc[i][ny];
			pc[i][0]=pc[i][1];
		}
	for(j=1;j<=ny;j++)
		{
			pc[0][j]=pc[1][j];
			pc[nx+1][j]=-pc[nx][j];
		}
		
		
	                                ap=dt*(2.0*dy/dx+2.0*dx/dy);	
					aw=-dt*dy/dx;
					ae=aw;
					as=-dt*dx/dy;
					an=as; 
	alp_num=0.0;
	for(i=1;i<=nx;i++)
		{
		for(j=1;j<=ny;j++)
			{
				Res[i][j]=t0[i][j]-(ap*t[i][j]+aw*t[i-1][j]+ae*t[i+1][j]+an*t[i][j+1]+as*t[i][j-1]);
				direc[i][j]=Res[i][j];
				Res_star[i][j]=Res[i][j];		
				alp_num=alp_num+Res[i][j]*Res_star[i][j];
			}
		}		
//    Next iterations
do
{

// Matrix multiplication [A]*[Res_0*] for alp_denominator
alp_den=0.0;
	for(i=1;i<=nx;i++)
		{
		for(j=1;j<=ny;j++)
			{
				
				if(i==1 && j==1) // SW
					AP[i][j]=ap*direc[i][j]+ae*direc[i+1][j]+an*direc[i][j+1];
				if(i==nx && j==1) // SE
					AP[i][j]=ap*direc[i][j]+aw*direc[i-1][j]+an*direc[i][j+1];	
				if(i==nx && j==ny) // NE 
					AP[i][j]=ap*direc[i][j]+aw*direc[i-1][j]+as*direc[i][j-1];
				if(i==1 && j==ny) // NW
					AP[i][j]=ap*direc[i][j]+ae*direc[i+1][j]+as*direc[i][j-1];			
				if(i==1 && (j>1 && j<ny)) //West
					AP[i][j]=ap*direc[i][j]+ae*direc[i+1][j]+an*direc[i][j+1]+as*direc[i][j-1];	
				if(i==nx && (j>1 && j<ny)) //East
					AP[i][j]=ap*direc[i][j]+aw*direc[i-1][j]+an*direc[i][j+1]+as*direc[i][j-1];		
				if(j==1 && (i>1 && i<nx)) //South
					AP[i][j]=ap*direc[i][j]+aw*direc[i-1][j]+ae*direc[i+1][j]+an*direc[i][j+1];			
				if(j==ny && (i>1 && i<nx)) //North
					AP[i][j]=ap*direc[i][j]+aw*direc[i-1][j]+ae*direc[i+1][j]+as*direc[i][j-1];
				if((1<i<nx) && (1<j<ny)) //Interior
					AP[i][j]=ap*direc[i][j]+aw*direc[i-1][j]+ae*direc[i+1][j]+as*direc[i][j-1]+an*direc[i][j+1];
							
				alp_den=alp_den+Res_star[i][j]*AP[i][j];				
			}
		}
// Finding alpha
Relaxation=alp_num/alp_den;

// Shadow element 2 calculation 
for(i=1;i<=nx;i++)
		{
		for(j=1;j<=ny;j++)
			{
				shad2[i][j]=Res[i][j]-Relaxation*AP[i][j];	// Updated pressure correction value
			}
		}	

// Matrix multiplication [A]*[s] for alp_denominator
om_num=0.0;
om_den=0.0;
	for(i=1;i<=nx;i++)
		{
		for(j=1;j<=ny;j++)
			{
				
				if(i==1 && j==1) // SW
					APshad[i][j]=ap*(shad2[i][j])+ae*(shad2[i+1][j])+an*(shad2[i][j+1]);
				if(i==nx && j==1) // SE
					APshad[i][j]=ap*(shad2[i][j])+aw*(shad2[i-1][j])+an*(shad2[i][j+1]);	
				if(i==nx && j==ny) // NE 
					APshad[i][j]=ap*(shad2[i][j])+aw*(shad2[i-1][j])+as*(shad2[i][j-1]);
				if(i==1 && j==ny) // NW
					APshad[i][j]=ap*(shad2[i][j])+ae*(shad2[i+1][j])+as*(shad2[i][j-1]);			
				if(i==1 && (j>1 && j<ny)) //West
					APshad[i][j]=ap*(shad2[i][j])+ae*(shad2[i+1][j])+an*(shad2[i][j+1])+as*(shad2[i][j-1]);	
				if(i==nx && (j>1 && j<ny)) //East
					APshad[i][j]=ap*(shad2[i][j])+aw*(shad2[i-1][j])+an*(shad2[i][j+1])+as*(shad2[i][j-1]);		
				if(j==1 && (i>1 && i<nx)) //South
					APshad[i][j]=ap*(shad2[i][j])+aw*(shad2[i-1][j])+ae*(shad2[i+1][j])+an*(shad2[i][j+1]);			
				if(j==ny && (i>1 && i<nx)) //North
					APshad[i][j]=ap*(shad2[i][j])+aw*(shad2[i-1][j])+ae*(shad2[i+1][j])+as*(shad2[i][j-1]);
				if((1<i<nx) && (1<j<ny)) //Interior
					APshad[i][j]=ap*(shad2[i][j])+aw*(shad2[i-1][j])+ae*(shad2[i+1][j])+as*(shad2[i][j-1])+an*(shad2[i][j+1]);
				om_num=om_num+APshad[i][j]*shad2[i][j];
				om_den=om_den+APshad[i][j]*APshad[i][j];
			}
		}
//Calculating omega
omega=om_num/om_den;


for(i=1;i<=nx;i++)
		{
		for(j=1;j<=ny;j++)
			{
				pc[i][j]=pc[i][j]+Relaxation*direc[i][j]+omega*(shad2[i][j]);	// Updated pressure correction value
				Res[i][j]=shad2[i][j]-omega*APshad[i][j];	// New Residue
			}
		}	
b_den=alp_num;
alp_num=0.0;
for(i=1;i<=nx;i++)
		{
		for(j=1;j<=ny;j++)
			{
				alp_num=alp_num+Res[i][j]*Res_star[i][j];		// updated value for next iteration
			}
		}		
direc_fac=(alp_num/b_den)*(Relaxation/omega);

for(i=1;i<=nx;i++)
		{
		for(j=1;j<=ny;j++)
			{
				
				direc[i][j]=Res[i][j]+direc_fac*(direc[i][j]-omega*AP[i][j]);	// updated value for the next iteration
			}
		}	
gRMS=sqrt(fabs(alp_num/(nx*ny)));
gcount++;
//if(gcount%3000==0)
printf("gRMS=%f\t\tgcount=%d\talp=%f\n",gRMS,gcount,Relaxation);
time=time+dt;	

}while(gRMS>0.000001);

Regards,
Ashish

[FMDenaro]

What about the BCs. for the pressure equation?

[arotester]

Originally Posted by FMDenaro

What about the BCs. for the pressure equation?

I have given homogenous Neumann at the inlet and the walls whereas the Dirichlet(zero) at the outlet.

[arotester]

Update:

I fixed the reference pressure near the outlet close to zero and applied the homogenous Neumann at all boundaries. This gave the correct answer with GSSOR, CGS and BiCGSTAB but the values of alpha and omega seem to be in the order of 200 to 300 which I believe is not correct.

Can someone please give an answer to this problem that why it is giving such a high value (it was showing such high values even before)?


Regards,
Ashish

[arotester]

Hello,

Does anyone have any feedback on the code????

[AliE]

Hi, unfortunately I have no time to check you code but honestly speaking... Do you really need to do your own BiCGSTAB? Moreover a preconditioner is essential to speed up the calculations!

There are several libraries that you can you "just" to solve a linear system. For example, I use PETSc which is simple and efficient. If you look online maybe you can find even something easier. In this case you will forget every detail related to the linear solver and you wil have more focus on your own code. Think about it

[arotester]

Originally Posted by AliE

Hi, unfortunately I have no time to check you code but honestly speaking... Do you really need to do your own BiCGSTAB? Moreover a preconditioner is essential to speed up the calculations!

There are several libraries that you can you "just" to solve a linear system. For example, I use PETSc which is simple and efficient. If you look online maybe you can find even something easier. In this case you will forget every detail related to the linear solver and you wil have more focus on your own code. Think about it

Thank you for the suggestion.
I will take a look at PETSc and its packages.
You are right that I can focus more on my own code if I use the solver packages but I have already written the code and it's working but there are some bugs in it. It will take me some time to debug the code and same time will be taken to learn to adopt PETSc.

So, debugging this code will be a priority and I am still trying for it.

Thank you,

Regards