[petsc-dev] Seeking Advice on Petsc Preconditioners to Try
Matthew Knepley
knepley at gmail.com
Fri Dec 2 18:16:22 CST 2011
On Fri, Dec 2, 2011 at 5:46 PM, Dave Nystrom
<Dave.Nystrom at tachyonlogic.com>wrote:
> 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.
>
No, this does not really help. For instance, from a Stokes problem we would
expect to hear, I am solving a Laplacian in the momentum equation, and
enforcing a differential constraint with the divergence. That allows us to
understand and build a preconditioner.
There is a simply staggeringly large literature on that problem. The time
has
come for a literature search.
It is very likely that sacusp was ineffective because you have a vector
problem,
but did not indicate this separation to the solver.
Matt
> 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
>
--
What most experimenters take for granted before they begin their
experiments is infinitely more interesting than any results to which their
experiments lead.
-- Norbert Wiener
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-dev/attachments/20111202/279042ae/attachment.html>
More information about the petsc-dev
mailing list