Further question about PC with Jaocbi Row Sum

[Barry Smith]

On Apr 8, 2008, at 7:16 PM, Shi Jin wrote:
> 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.

  The lumped mass matrix being a good replacement for the mass matrix  
is a question about approximation. How good is each to the continuous  
L_2 norm? It isn't
really a question about how close each is to the other.

Being a good preconditioner is a linear algebra question, what is the  
(complex) relationship between the eigenvalues and eigenvectors of the  
two matrices?  (what happens to the eigenvalues of B^{-1} M?)

I think these two questions are distinct, intuitively they
seem to be related, but mathematically I don't think there is a direct  
relationship so I am not surprised by your observations.


    By the way, I have seen cases where using the lumped mass matrix  
resulted in BETTER approximation to the continuous solution then using  
the
"true" mass matrix; again this is counter intuitive but there is  
nothing mathematically that says it shouldn't be.

>
>
> 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?

    Yes, so long as you use a tight enough convergence tolerance with  
KSPSetTolerances(). Also by default with most KSP solvers the  
"preconditioned" residual
norm is used to determine convergence, thus in some way the  
preconditioner helps determines when the KSP stops.

> Is there something particular about this rowsum method?

    No. If you use a -ksp_rtol of 1.e-12 and still get different  
answers, this needs to be investigated.

    Barry

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