[petsc-dev] Seeking Advice on Petsc Preconditioners to Try

[Dave Nystrom]

Matthew Knepley writes:
 > On Wed, Nov 30, 2011 at 12:41 AM, Dave Nystrom <dnystrom1 at comcast.net>wrote:
 > 
 > > I have a linear system in a code that I have interfaced to petsc that is
 > > taking about 80 percent of the run time per timestep.  This linear system
 > > is a symmetric block banded matrix where the blocks are 2x2.  The matrix
 > > looks as follows:
 > >
 > >  1 2 3 4 5 6 7 8 9 0 1 2 3 4 5 6 7 8 9 0
 > > 1X X                     Y Y Y
 > > 2X X X                     Y Y Y
 > > 3  X X X                     Y Y Y
 > > 4    X X X                     Y Y Y
 > > 5      X X X                     Y Y Y
 > > 6        X X X                     Y Y Y
 > > 7          X X X                     Y Y Y
 > > 8            X X X                     Y Y Y
 > > 9              X X X                     Y Y Y
 > > 0                X X X                     Y Y Y
 > > 1                  X X X                     Y Y Y
 > > 2                    X X X                     Y Y Y
 > > 3Z                     X X X                     Y Y Y
 > > 4Z Z                     X X X                     Y Y Y
 > > 5Z Z Z                     X X X                     Y Y Y
 > > 6  Z Z Z                     X X X                     Y Y Y
 > > 7    Z Z Z                     X X X                     Y Y Y
 > > 8      Z Z Z                     X X X                     Y Y Y
 > > 9        Z Z Z                     X X X                     Y Y Y
 > > 0          Z Z Z                     X X X                     Y Y Y
 > >
 > > So in my diagram above, X, Y and Z are 2x2 blocks.  The symmetry of the
 > > matrix requires that X_ij = transpose(X_ji) and Y_ij = transpose(Z_ji).  So
 > > far, I have just input this matrix to petsc without indicating that it was
 > > block banded with 2x2 blocks.  I have also not told petsc that the matrix
 > > is symmetric.  And I have allowed petsc to decide the best way to store the
 > > matrix.
 > >
 > > I can solve this linear system over the course of a run using -ksp_type
 > > preonly -pc_type lu.  But that will not scale very well to larger problems
 > > that I want to solve.  I can also solve this system over the course of a
 > > run using -ksp_type cg -pc_type jacobi -vec_type cusp -mat_type aijcusp.
 > > However, over the course of a run, the iteration count ranges from 771 to
 > > 47300.  I have also tried sacusp, ainvcusp, sacusppoly, ilu(k) and icc(k)
 > > with k=0.  The sacusppoly preconditioner fails because of a thrust error
 > > related to an invalid device pointer, if I am remembering correctly.  I
 > > reported this problem to petsc-maint a while back and have also reported it
 > > for the cusp bugtracker.  But it does not appear that anyone has really
 > > looked into the bug.  For the other preconditioners of sacusp, ilu(k) and
 > > icc(k), they do not result in convergence to a solution and the runs fail.
 > 
 > All preconditioners are custom.  Have you done a literature search for PCs
 > known to work for this problem?  Can you say anything about the spectrum of
 > the operator? conditioning? what is the principal symbol (if its a PDE)?
 > The pattern is not enough to recommend a PC.

Hi Matt,

Thanks for your reply and sorry to take so long to get back to you.  I have
not done a literature search.  This is a time stepping code and so I have
different linear systems for each linear solve and time step.  I have taken
some time to study this particular system at a point in the nonlinear phase
of the problem although there is definitely more I could do in that area.  At
one point in the run, I have taken the matrix and generated a rhs using a
known solution vector of random numbers between 0 and 1.  Using some petsc
info that I believe Jed supplied to a previous question, I used petsc to
print out the condition number estimate for the preconditioned system for
various preconditioners.  The results seemed to be in the range of 10^5 to
10^7.  Also, when I computed error norms using a direct petsc LU solve, I was
only getting 11-12 digits of accuracy.  On my other linear systems, I would
get 15-16 digits of accuracy.  So the system is poorly conditioned.  However,
the operator is computed from fields that are highly non-uniform over the
spatial domain of the problem.  The problem is a formation problem for a
Field Reversed Configuration (FRC) magnetic fusion problem.

The matrix is symmetric but I am not sure if it is positive definite.  I can
check symmetry easily with petsc but am not sure if there is an easy way to
tell if the matrix is positive definite - if there is a good way to do this
with petsc, I'd love to hear about it.  I did not find anything in the online
help.  I wonder if petsc might be able to say something about the positive
definiteness of the system after a solution has been found.  Anyway, I'm not
up on the various ways to tell if a matrix is positive definite but I suppose
it would be expensive to solve for all the eigenvalues of the matrix just to
tell if it was positive definite.  If there are clever and inexpensive ways
to do this, I'd be interested in trying them.

However, I did try both minres and symmlq from petsc using jacobi
preconditioning and cuda.  The case where I used minres ran to completion and
used less than half the iterations that cg used so I was encouraged about
that.  However, there was an issue with the final solution that I need to
investigate.  The run using symmlq from petsc with jacobi preconditioning and
cuda had numerical problems that I have not investigated and did not run to
completion.

So, because of the physics involved, the linear system is difficult to
solve.  It seems poorly conditioned.  It is symmetric but might not be
positive definite.  The matrix comes from discretizing a pde on a regular
structured 2d mesh.

Does any of this additional info help?

Thanks,

Dave

 > Matt
 > 
 > > I'm wondering if there are suggestions of other preconditioners in petsc
 > > that I should try.  The only third party package that I have tried is the
 > > txpetscgpu package.  I have not tried hypre or any of the multigrid
 > > preconditioners yet.  I'm not sure how difficult it is to try those
 > > packages.  Anyway, so far I have not found a preconditioner available in
 > > petsc that provides a robust solution to this problem and would be
 > > interested in any suggestions that anyone might have of things to try.
 > >
 > > I'd be happy to provide additional info and am planning on packaging up a
 > > couple of examples of the matrix and rhs for people I am interacting with
 > > at Tech-X and EMPhotonics.  So I'd be happy to provide the matrix
 > > examples for this forum as well if anyone wants a copy.
 > >
 > > Thanks,
 > >
 > > Dave
 > >
 > > --
 > > Dave Nystrom
 > >
 > > phone: 505-661-9943 (home office)
 > >       505-662-6893 (home)
 > > skype: dave.nystrom76
 > > email: dnystrom1 at comcast.net
 > > smail: 219 Loma del Escolar
 > >       Los Alamos, NM 87544