Further question about PC with Jaocbi Row Sum

[Shi Jin]

Hi there,

First of all, I want to thank Matt  for his previous suggestion on the use of -pc_jacobi_rowsum 1 option. Now I have a more theoretical question I hope you may help although it does not necessarily connects to PETsc directly.

Basically, I am trying to speed up the solution of a finite element mass matrix, which is constructed on a second order tetrahedral element. My idea is to use the lumped mass matrix to precondition it. However, it does not seem to work originally but I guessed it may has something to do with the fact that it is second order since the theory for second order lumped mass matrix is not so clear, at least to me. So I decided to work out the linear element first, since everything there is mathematically well established.

OK. I first constructed a mass matrix based on linear elements. I can construct its lumped mass matrix by three methods:
1. sum of each row
2. scale the diagonal entry
3. use a nodal quadrature rule to construct a diagonal matrix
These three methods turned out to produce identical diagonal matrices as the lumped mass matrix, just as the theory predicts. So perfect!
I can further test that the lumped mass matrix has similar eigenvalues to the original consistent mass matrix, although different. And the solutions of a linear system is quite similar too. So I understand the reason lots of people replace solving the consistent mass matrix with solving the lumped one to achieve much improved efficiency but without losing much of the accuracy. So far, it is all making sense.

I naturally think it could be very helpful if I can use this lumped mass matrix  as the prediction  matrix  for  the  consistent  mass matrix solver,  where  the  consistent  matrix  is  kept  to  have the  better accuracy. However, my tests show that it does not help at all. I tried several ways to do it within PETSc, such as setting
KSPSetOperators( solM, M, lumpedM SAME_PRECONDITIONER);
or directly use the -pc_type  jacobi -pc_jacobi_rowsum 1 to the built-in method.
These two methods turned out to be equivalent but they both produce less efficient solutions: it actually took twice more steps to converge than without these options.
This is quite puzzling to me. Although I have to admit that I have seen a lot of replacing the consistent mass matrix with the lumped one in the literature but have not seen much of using the lumped mass matrix as a preconditioner. Maybe using the lumped matrix for preconditioning simply does not work? I would love to hear some comment here.

If that's all, I don't feel too bad. Then I came back to the second order elements since that's what I want to use. Accidentally  I decided to try to solve the  second  order  consistent mass matrix with the   -pc_type  jacobi -pc_jacobi_rowsum 1 option. Bang! It converges almost three times faster.  For a particular system, it usually converges in 9-10 iterations and now it converges in 2-3 iterations.  This is amazing! But I don't know why it is so. 

If that's all, I would be just happy. Then I ran my single particle sedimentation code with -pc_type  jacobi -pc_jacobi_rowsum 1 and it does run a lot faster. However, the results I got are slightly different from what I used to, which is weird since the only thing changed is the preconditioner while the same linear system was solved. I tried several -pc_type options, and they are all consistent with the old one. So I am a little bit hesitant adapting this new speed up method. What troubles me most is that the new simulation results are actually even closer to the experiments we are comparing, which may suggest that the row sum PC is even better. But this is just one test case I would rather believe it happens to cause errors in the direction to compensate  other simulation errors. If I had other well established test case which unfortunately I don't, I would imagine it may work differently.

So my strongest puzzle is that how could a change in the pre-conditioner make such an observable change in the solutions. I understand different PCs produce different solutions but they should be numerically very close and non-detectable on a physical quantity level plot, right? Is there something particular about this rowsum method?

I apologize about this lengthy email but I do hope to have some in depth scientific discussion.
Thank you very much.

Shi Jin, PhD




      ____________________________________________________________________________________
You rock. That's why Blockbuster's offering you one month of Blockbuster Total Access, No Cost.  
http://tc.deals.yahoo.com/tc/blockbuster/text5.com