Reg. BICGSTAB
 

Hi, i wrote a code for solving Pressure Poisson equation using BiCgSTAB iterative technique. My code is working fine except for the very fine grids(dx~1e-3) so i cant run some cases. Plz advice me reg. this. Is pre-conditioning of any use..?

Re: Reg. BICGSTAB
 

Hi,

It's known that BI-CGStab can breakdown, I found in a paper the reasons, there are 3, but I can't remember...

But basically, Pressure Poisson Equations matyrix is supposed to be symmetric and positive definite!!! And BiCGStab is designed for unsymmetric matrices!!! So why not use simply CG or CGS. With or whitout preconditionning, you will get very fast and accurate results.

If your matrice is not symmetric, you will probably have to think about GMRES. In this case you have to choose a fairly large restart but it should be ok. Not as fast as BICGStab but you'll always get good results.

Re: Reg. BICGSTAB
 

Hi, Thank you for your suggestions. If possible please give the details of that paper.

Re: Reg. BICGSTAB
 

The title of the paper is "The Breakdowns of BiCGStab" Link: http://www.springerlink.com/content/62w0t6nh2p8frala/

Re: Reg. BICGSTAB
 

If you're not preconditioning the matrix then it's likely that bicgstab will not converge (well not within a reasonable number of iterations anyway).

Things to check:

(1) make sure you are running in double precision and your convergence criteria is sensibly set.

(2) Try diagonal preconditioning to see if it improves matters.

(3) If (2) appears to work then try another simple preconditioner (Jacobi, Guass-Seidel, SOR or SSOR are all easy to program).

(4) If all the above fail try another solver such as GCR(k) or GMRES(k).

NOTE: most iterative methods have the potential to breakdown/stagnate (it's happened to me with both GCR and GMRES even with preconditioning).

Re: Reg. BICGSTAB
 

The most simple preconditionner you can use is the diagonal of the matrix. It's improving the simulation enough to see if you need to invesigate more the idea of the preconditionner or simply move to another method.

Re: Reg. BICGSTAB
 

hi, Could you tell us your grid size? I am using BiCGSTAB (without preconditioning) to solve pressure poisson equation with 550*450 grid points. It is converging and working very wel.. It is good if you once again check your implementation.

Re: Reg. BICGSTAB
 

I also used BICGStab for another kind of matrix, unsymmetric of course. Diverging for n=18000 if no preconditionning Diverging for n=32000 if diagonal preconditionning n being the size of the matrix of course. So it is a fact : BICGStab can breakdown for numerous reasons!!

Re: Reg. BICGSTAB
 

N=18000 or N=32000 is not important here, and it is not the reason it has broken down. And since it has broken down diagonal preconditioning is not gonna help because your matrix clearly is not diagonally dominant (not a M matrix). The reason for preconditioned BiCGStab is not very clear from what you have written. (I had made BiCGStab failed for even N=6).

If your case is too ill conditioned the fail safe method is LSQR.

Re: Reg. BICGSTAB
 

Do you actually mean diverging or do you mean "not converging" - they are different things?

The only time I've had "divergence" has been to programming error (or bad initialization resulting in a division by zero).

The not converging case occurs when the iteration increments either become zero (stagnation) or the residual oscillates without decreasing.

Re: Reg. BICGSTAB
 

The problem is coming when i put, say, 21 grids in 1/600cm(to capture some details in a case, i've used clustered grids.Total domain size was 5.5*0.1 mm with 550*62 grids),in short, when the dx or dy are of very low magnitude. Except for such kind of fine grids, it is working for any no. of grid points, so where could possibly the bug be..?. Please suggest how to modify the implementation as it is just being influenced by the grid dimension?

Re: Reg. BICGSTAB
 

(in continuation)

My code is Diverging for such grids.

Re: Reg. BICGSTAB
 

What I meant is that for my case, the BIGgstab was giving great results for "coarse" meshes and was breacking down for "finer" mesh. And my matrix is diagonal dominant... But currently I'm using SuperLU or GMRES, which is giving great satisfation, specially the first.

Re: Reg. BICGSTAB
 

BiCGStab with good preconditioner is very good. I give you an idea about how fast it could be : on a unstructured grid of 70 000 cells, and heat transfer equation with non linear source term. 1 iteration of it reduced the absolute error by factor of 10^-7. That is convergent in 1 iteration. The preconditioner was AMG. ;-)

So if the preconditioner is good, BiCGStab is also very good and robust.

Re: Reg. BICGSTAB
 

hi InddZen,

Is it true that the discretization of pressure Poisson equation always produces symmetric and positive definite matrix? Is this the case even for non-uniform grid? is CG faster than multigrid?

Re: Reg. BICGSTAB
 

No, afaik,For non-uniform/non-orthogonal grid the discretization will not produce symmetric matrix.

Re: Reg. BICGSTAB
 

Actually for incompressible case the matrix will be symmetrical and positive definite. Even for unstructured grids. CG could be used and is used usually, but i have found that BiCGStab (ILU preconditioner) converges much faster than CG. So I always use it.

Re: Reg. BICGSTAB
 

Ok, AMG is working fine for elliptic equations, what about other kind of matrices, for exemple the one obtained by a FEM discretisation of Navier Stokes equation (incompressible flow of course), specially for high Reynold Numbers... You should read again your papers... Actually I don't have time to spend on trying to developp the right AMG for my matrix, which is really special, probaly latter...

Re: Reg. BICGSTAB
 

"Actually for incompressible case the matrix will be symmetrical and positive definite."

No it won't - try discretizing the 1D problem on a nonuniform grid using finite differences!

Re: Reg. BICGSTAB
 

Try nondimensionalizing your equations and make sure you're using a relative (dimensionless) error estimate for convergence. Also make sure you're using double precision.

Re: Reg. BICGSTAB
 

Yes, it will. Try descretization by finite volume method. lol.

I will give you a quote from Peric's book:

"; moreover, its coefficient matrix is symmetric so special solvers for symmetric matrices can be used (e.g ICCG solver from the conjugate gradients family ...... )"

Chapter 8. Complex Geometries.

I do not need Peric to tell me that it will be symmetric, I have already written code for unstructured grids and I know it is symmetric for incompressible cases. It is unsymmetric for compressible because it has one more term due to density correction.

Re: Reg. BICGSTAB
 

You're wrong, try the suggested 1D example on a non-uniform grid - you'll end up with something of the form

a.P_i-1 + b.P_i + c.P_i+1 = f_i

with, in general, a /= c => matrix is not symmetric!

Re: Reg. BICGSTAB
 

Ok I am wrong , you won, happy now.

Re: Reg. BICGSTAB
 

You need to consider symmetry about the diagonal i.e. does A_i,j = A_j,i ?

An off-diagonal entry in row i column j of a symmetric matrix correspond with another in a different row. Your example is limited to terms in a single row of the matrix, so no conclusion on symmetry can be drawn from your example.

Infact, you're just considering if A_(i-b),i = A_(i+b),i. That's a very different property and is not generally true.

andy

Re: Reg. BICGSTAB
 

You need a(i+1) = c(i) for all the i's which in general is not true on a uniform grid - I just missed the indices off. If h_i is the spacing c(i) is calulated using h_(i-1) and h_i while a(i+1) will be calculated with h_i and h_(i+1). This is why it will not be symmetric.

Re: Reg. BICGSTAB
 

I will help you little bit:

http://en.wikipedia.org/wiki/Symmetric_matrix

Re: Reg. BICGSTAB
 

Ok... I'm still confused. Please can you point out the error in the following, as I cannot spot it myself!!

In general, for any node index k:

Rule 1) a_k references the node/cell to the left of node k and uses h_k and h_(k-1)

Rule 2) c_k references the node/cell to the right of node k and uses h_k and h_(k+1)

Now for symmetry, we require a_(j+1) = c_j, as you say.

By applying rule (1) above we find that a_(j+1) uses h_(j+1) and h_j.

By appying rule (2) above we find that c_j uses h_j and h_(j+1).

So ISTM we are using the same h's. This again indicates symmetry.

I hope we can get to the bottom of this!! I certainly thought this matrix should be symmetric - it's a Poisson equation so there's no upwind terms to break the symmetry, but I am open to evidence to the contary.

Best regards, Andy

Re: Reg. BICGSTAB
 

You should ask yourself if you're a better mathematician than me before posting such things!

Check my counter example more closely (the matrix is tridiagonal so symmetry requires a(i+1)=c(i) as I stated) or give a completely rigourous proof that the discrete form is always symmetric no matter what numerical technique is used! - hint you will not be able do this because there is a counter example.

Re: Reg. BICGSTAB
 

I pointed you to the wiki because you do not understand the basic definition of symmetric matrix which requires that for the matrix to be symmetric Aij = Aji .

That is the other element that should be equal to Aij does not even fall on same raw. (Aji is not present in the raw that contains Aij, in the matrix).

Now what you keep on writing and stubbornly arguing is the comparing the elements on the same raw and telling that since they are not same matrix is not symmetric.

I gave you a quote from Peric's book but you still think you are more intellegent than Peric.

From the wiki link I pointed you, an simple example of symmetric matrix:

| 1 2 3 |

| 2 4 -5 |

| -4 5 0 |

this is symmetric matrix and it is not symmetric according to your logic , since in raw 2 : 2 != -5

Learn basic maths first, will help you a lot in life.

correction
 

the matrix i quoted should read

| 1 2 3 |

| 2 4 -5 |

| 3 -5 6 |

Just a small typo.

Re: Reg. BICGSTAB
 

It's the fact that the three h's are not equal for a stretched grid that cause's the asymmetry. You really can't talk about the symmetry of the matrix without considering the finite difference grid.

The point is, for example,

c(i) = 2.0/( h(i)*(h(i) + h(i-1)) ) while a(i) = 2.0/( h(i-1)*(h(i) + h(i-1)) )

=> a(i+1) = 2.0/( h(i)*( h(i+1) + h(i)) ) /= c(i)

in general because the h's are not equal! More complex/accurate discretizations are even less likely to be symmetric. The above is an example of what can happen in the incompressible equations on a staggered grid.

Re: Reg. BICGSTAB
 

Thanks - that makes perfect sense now.

I've programmed this myself in the not so distant past, so I must have known this quite recently... I wonder why I forgot!! Good to revisit it I guess.

Best wishes, andy

Re: Reg. BICGSTAB
 

The pressure correction equation is poisson type equation, and it is not discretized the way you have written.

between i to j (the coefficient relating i to j) Aij is given as

Aij = Area / distance_ij

in 1D case area = 1 , so Aij = 1/distance_ij

for j relating to i, Aji = Area / distance_ji

Since distance_ij = distance_ji (irrespective of how other elements are speard out, distance between two points is what decide this coefficient. So your non uniformity is grid spacing does not effect this).

Aij = Aji and hence matrix is symmetric.

If you want to descretisize the way you have written feel free to do so, but softwares like Fluent StarCCM do not do what you have written. (But do as I mentioned above).

For incompressible cases, pressure matrix is symmetric. (finite volume methods).

Re: Reg. BICGSTAB
 

No the symmetry of the discretized equation depends on the grid and your choice of discretization - This is why when Ferziger and Peric discuss the pressure projection equation they point out that you can either use a defect correction form (thus lagging any asymmetry) or throw away some terms (namely the nonorthogonal ones).

The discretization I've giving is pretty much standard on staggered grids with grid strecthing when the grid changes are smooth. Convert the Laplacian operator into general nonorthogonal curvilinear coordinates and then convince yourself that it's still a symmetric operator.

Re: Reg. BICGSTAB
 

You want to say that what is written here is wrong:

http://www.cfd-online.com/Wiki/Discr...diffusion_term

And the non orthoginal terms go into the source part of equations and not in the matrix itself.

Clearly the matrix is symmetric. Of course there is no way you would agree with it, since you have decided to ignore all that is presented to you. Best of luck with your quest in CFD.

Re: Reg. BICGSTAB
 

And you're ignoring the facts (if you move the non-othorgonal terms, and all the pressure terms making up the extended stencilm, from the source term then then they contribute to the matrix and make it asymmetric!)

"Best of luck with your quest in CFD."

You really don't know anything about me - I'm not some "student" and have a job writing and developing CFD codes. You may also consider I've got a BSc, MSc and PhD in Mathematics + numerous years of work experience (on top of postdoctoral work).

My Last word on this pointless argument:

There is no one correct way to do CFD - if there was there would be no need to do research (hence your examples are useless - you need to demonstrate that "all possible discretizations" of the equations leads to a completely symmetric equation (no lagging terms in the iteration) which as I've shown, by counter example, is false).

Re: Reg. BICGSTAB
 

1. Where are non orthogonal terms in 1D case.

2. We are not talking about how you would do things, we are talking about how they are generally done. You could write second derivative using many neighbors and claim that it would be asymetric. But this is not the way it is done. It is done the way I have pointed you to the wiki page. And it comes out to be symmetric. To prove it otherwise you need to prove that what I pointed you is wrong. Go ahead and do it. I am waiting.

3. No matter how many degrees you have, it is the fact that the approach Fluent, StarCCM+ takes for pressure descretisation, the matrix comes out to be symmetric. (for incompressible cases). And this is what we been arguing. So your saying that pressure matrix is not symmetric stands wrong.

Re: Reg. BICGSTAB
 

All, BICGSTAB is good if your matrix is not symmetric pos. def. However, if your matrix is sym. pos. def. you are wasting computational resources by doing the second matrix-vector product dealing with the transpose. If you are not seeing oscillations in the relative residual, BICG may be faster if you are set on using such a method.

Furthermore, a good preconditioner is essential when using an of the sparse solvers. The main diagonal preconditioner is Jacobi preconditioning and really only works well when the matrix is really diagonally dominant (i.e., the matrix 1/d_ii is a good approximation of the inverse). Other good preconditioners are Incomplete Cholesky which should be used with sym. pos. def. matrices and Incomplete LU should be used otherwise. There are various implementations of these to consider (i.e., zero fill-in, Crout type, conservation of row and column sums). If your matrix is composed of blocks of smaller matrices then using a block Jacobi preconditioner can be a good option as well. Depending on the PDE of interest, there are many specialized preconditioners available that are designed for the sparsity pattern of your system. There are mathematical conferences dedicated to this subject.

From my experience, if your matrix is not sym. pos. def. then GMRES(m) and choosing (m<=50) minimizes the storage associated with the search directions in the Krylov subspace is a good method. If you use a good incomplete preconditoner (one that conserves row-sums) you can increase your convergence by an order of magnitude in some cases. However, choosing (m) is somewhat of an art since many matrices that are poorly conditioned can result in stagnation of the relative residual. The idea is to choose (m) to minimize storage and still get rapid convergence. In some instances you can get convergence using (m) restart whereas you do not with GMRES w/out restart. If you have a sym. pos. def. matrix then CG is one of the fastest (since it does a considerable less amount of computation), especially with a good preconditioner like ichol. You can determine how many iterations that CG should take with knowledge of the condition number and a little bit of algebra.

In general, the worse the condition number (~ the ratio of the largest to smallest eigenvalue), the more problems you will have finding a good combination of a solver and preconditioner to optimize your cpu effort.

There are several other iterative solvers that work well for large sparse systems. I recommend Y. Saad's text and C.T. Kelley's text on the subject. Both give good examples using common sparse systems we see in the discretization of PDEs. Saad's text also includes a detailed description of various types of preconditioners.

-- Regards,

DR