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
