<div class="gmail_quote">On Tue, Dec 20, 2011 at 09:33, Dave Nystrom <span dir="ltr"><<a href="mailto:dnystrom1@comcast.net">dnystrom1@comcast.net</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
Hi Mark,<br>
<br>
I would like to try GAMG on some of my linear solves.  Could you suggest how<br>
to get started?  Is it more complicated than something like:<br>
<br>
-ksp_type cg -pc_type gamg<br></blockquote><div><br></div><div>Sure, start with this.</div><div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
<br>
I'm guessing I should first try it on one of my easier linear solves.  I have<br>
5 of them that would have a block size of 1.  Are the other GAMG option<br>
defaults good to start with or should I be trying to configure them as well?<br>
If so, I'm not familiar enough with multigrid to know off hand how to do<br>
that.<br></blockquote><div><br></div><div>What sort of problems? The defaults should be reasonable for elliptic problems.</div><div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">

<br>
Thanks,<br>
<font color="#888888"><br>
Dave<br>
</font><div><div></div><div class="h5"><br>
Mark F. Adams writes:<br>
 ><br>
 > On Dec 2, 2011, at 6:06 PM, Dave Nystrom wrote:<br>
 ><br>
 > > Mark F. Adams writes:<br>
 > >> It sounds like you have a symmetric positive definite systems like du/dt -<br>
 > >> div(alpha(x) grad)u.  The du/dt term makes the systems easier to solve.<br>
 > >> I'm guessing your hard system does not have this mass term and so is<br>
 > >> purely elliptic.  Multigrid is well suited for this type of problem, but<br>
 > >> the vector nature requires some thought.  You could use PETSc AMG -pc_type<br>
 > >> gamg but you need to tell it that you have a system of two dof/vertex.<br>
 > >> You can do that with something like:<br>
 > >><br>
 > >> ierr = MatSetBlockSize( mat, 2 );      CHKERRQ(ierr);<br>
 > >><br>
 > >> For the best results from GAMG you need to give it null space information<br>
 > >> but we can worry about that later.<br>
 > ><br>
 > > Hi Mark,<br>
 > ><br>
 > > I have been interested in trying some of the multigrid capabilities in<br>
 > > petsc.  I'm not sure I remember seeing GAMG so I guess I should go look for<br>
 > > that.<br>
 ><br>
 > GAMG is pretty new.<br>
 ><br>
 > > I have tried sacusp and sacusppoly but did not get good results on<br>
 > > this particular linear system.<br>
 > > In particular, sacusppoly seems broken.  I<br>
 > > can't get it to work even with the petsc src/ksp/ksp/examples/tutorials/ex2.c<br>
 > > example.  Thrust complains about an invalid device pointer I believe.<br>
 > > Anyway, I can get the other preconditioners to work just fine on this petsc<br>
 > > example problem.  When I try sacusp on this matrix for the case of generating<br>
 > > a rhs from a known solution vector, the computed solution seems to diverge<br>
 > > from the exact solution.  We also have an interface to an external agmg<br>
 > > package which is not able to solve this problem<br>
 > > but works well on the other 5<br>
 > > linear solves.  So I'd like to try more from the multigrid toolbox but do not<br>
 > > know much about how to supply the extra stuff that these packages often need.<br>
 > ><br>
 > > So, it sounds like you are suggesting that I try gamg and that I could at<br>
 > > least try it out without having to initially supply lots of additional info.<br>
 > > So I will take a look at gamg.<br>
 > ><br>
 ><br>
 > There are many things that can break a solver but most probably want to know that its a system so if you can set the block size and try gamg then that would be a good start.<br>
 ><br>
 > Mark<br>
 ><br>
 > > Thanks,<br>
 > ><br>
 > > Dave<br>
 > ><br>
 > >> Mark<br>
 > >><br>
 > >> On Nov 30, 2011, at 8:15 AM, Matthew Knepley wrote:<br>
 > >><br>
 > >>> On Wed, Nov 30, 2011 at 12:41 AM, Dave Nystrom <<a href="mailto:dnystrom1@comcast.net">dnystrom1@comcast.net</a>> wrote:<br>
 > >>> I have a linear system in a code that I have interfaced to petsc that is<br>
 > >>> taking about 80 percent of the run time per timestep.  This linear system is<br>
 > >>> a symmetric block banded matrix where the blocks are 2x2.  The matrix looks<br>
 > >>> as follows:<br>
 > >>><br>
 > >>> 1 2 3 4 5 6 7 8 9 0 1 2 3 4 5 6 7 8 9 0<br>
 > >>> 1X X                     Y Y Y<br>
 > >>> 2X X X                     Y Y Y<br>
 > >>> 3  X X X                     Y Y Y<br>
 > >>> 4    X X X                     Y Y Y<br>
 > >>> 5      X X X                     Y Y Y<br>
 > >>> 6        X X X                     Y Y Y<br>
 > >>> 7          X X X                     Y Y Y<br>
 > >>> 8            X X X                     Y Y Y<br>
 > >>> 9              X X X                     Y Y Y<br>
 > >>> 0                X X X                     Y Y Y<br>
 > >>> 1                  X X X                     Y Y Y<br>
 > >>> 2                    X X X                     Y Y Y<br>
 > >>> 3Z                     X X X                     Y Y Y<br>
 > >>> 4Z Z                     X X X                     Y Y Y<br>
 > >>> 5Z Z Z                     X X X                     Y Y Y<br>
 > >>> 6  Z Z Z                     X X X                     Y Y Y<br>
 > >>> 7    Z Z Z                     X X X                     Y Y Y<br>
 > >>> 8      Z Z Z                     X X X                     Y Y Y<br>
 > >>> 9        Z Z Z                     X X X                     Y Y Y<br>
 > >>> 0          Z Z Z                     X X X                     Y Y Y<br>
 > >>><br>
 > >>> So in my diagram above, X, Y and Z are 2x2 blocks.  The symmetry of the<br>
 > >>> matrix requires that X_ij = transpose(X_ji) and Y_ij = transpose(Z_ji).  So<br>
 > >>> far, I have just input this matrix to petsc without indicating that it was<br>
 > >>> block banded with 2x2 blocks.  I have also not told petsc that the matrix is<br>
 > >>> symmetric.  And I have allowed petsc to decide the best way to store the<br>
 > >>> matrix.<br>
 > >>><br>
 > >>> I can solve this linear system over the course of a run using -ksp_type<br>
 > >>> preonly -pc_type lu.  But that will not scale very well to larger problems<br>
 > >>> that I want to solve.  I can also solve this system over the course of a run<br>
 > >>> using -ksp_type cg -pc_type jacobi -vec_type cusp -mat_type aijcusp.<br>
 > >>> However, over the course of a run, the iteration count ranges from 771 to<br>
 > >>> 47300.  I have also tried sacusp, ainvcusp, sacusppoly, ilu(k) and icc(k)<br>
 > >>> with k=0.  The sacusppoly preconditioner fails because of a thrust error<br>
 > >>> related to an invalid device pointer, if I am remembering correctly.  I<br>
 > >>> reported this problem to petsc-maint a while back and have also reported it<br>
 > >>> for the cusp bugtracker.  But it does not appear that anyone has really<br>
 > >>> looked into the bug.  For the other preconditioners of sacusp, ilu(k) and<br>
 > >>> icc(k), they do not result in convergence to a solution and the runs fail.<br>
 > >>><br>
 > >>> All preconditioners are custom. Have you done a literature search for PCs<br>
 > >>> known to work for this problem? Can yu say anything about the spectrum of the<br>
 > >>> operator? conditioning? what is the principal symbol (if its a PDE)? The pattern<br>
 > >>> is not enough to recommend a PC.<br>
 > >>><br>
 > >>>   Matt<br>
 > >>><br>
 > >>> I'm wondering if there are suggestions of other preconditioners in petsc that<br>
 > >>> I should try.  The only third party package that I have tried is the<br>
 > >>> txpetscgpu package.  I have not tried hypre or any of the multigrid<br>
 > >>> preconditioners yet.  I'm not sure how difficult it is to try those<br>
 > >>> packages.  Anyway, so far I have not found a preconditioner available in<br>
 > >>> petsc that provides a robust solution to this problem and would be interested<br>
 > >>> in any suggestions that anyone might have of things to try.<br>
 > >>><br>
 > >>> I'd be happy to provide additional info and am planning on packaging up a<br>
 > >>> couple of examples of the matrix and rhs for people I am interacting with at<br>
 > >>> Tech-X and EMPhotonics.  So I'd be happy to provide the matrix examples for<br>
 > >>> this forum as well if anyone wants a copy.<br>
 > >>><br>
 > >>> Thanks,<br>
 > >>><br>
 > >>> Dave<br>
 > >>><br>
 > >>> --<br>
 > >>> Dave Nystrom<br>
 > >>><br>
 > >>> phone: <a href="tel:505-661-9943" value="+15056619943">505-661-9943</a> (home office)<br>
 > >>>      <a href="tel:505-662-6893" value="+15056626893">505-662-6893</a> (home)<br>
 > >>> skype: dave.nystrom76<br>
 > >>> email: <a href="mailto:dnystrom1@comcast.net">dnystrom1@comcast.net</a><br>
 > >>> smail: 219 Loma del Escolar<br>
 > >>>      Los Alamos, NM 87544<br>
 > >>><br>
 > >>><br>
 > >>><br>
 > >>> --<br>
 > >>> What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br>
 > >>> -- Norbert Wiener<br>
 > >><br>
 > ><br>
 ><br>
</div></div></blockquote></div><br>