On Fri, Dec 2, 2011 at 5:46 PM, Dave Nystrom <span dir="ltr"><<a href="mailto:Dave.Nystrom@tachyonlogic.com">Dave.Nystrom@tachyonlogic.com</a>></span> wrote:<br><div class="gmail_quote"><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex;">
Matthew Knepley writes:<br>
> On Wed, Nov 30, 2011 at 12:41 AM, Dave Nystrom <<a href="mailto:dnystrom1@comcast.net">dnystrom1@comcast.net</a>>wrote:<br>
><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<br>
> > is a symmetric block banded matrix where the blocks are 2x2. The matrix<br>
> > looks 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<br>
> > is 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<br>
> > run 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 you say anything about the spectrum of<br>
> the operator? conditioning? what is the principal symbol (if its a PDE)?<br>
> The pattern is not enough to recommend a PC.<br>
<br>
Hi Matt,<br>
<br>
Thanks for your reply and sorry to take so long to get back to you. I have<br>
not done a literature search. This is a time stepping code and so I have<br>
different linear systems for each linear solve and time step. I have taken<br>
some time to study this particular system at a point in the nonlinear phase<br>
of the problem although there is definitely more I could do in that area. At<br>
one point in the run, I have taken the matrix and generated a rhs using a<br>
known solution vector of random numbers between 0 and 1. Using some petsc<br>
info that I believe Jed supplied to a previous question, I used petsc to<br>
print out the condition number estimate for the preconditioned system for<br>
various preconditioners. The results seemed to be in the range of 10^5 to<br>
10^7. Also, when I computed error norms using a direct petsc LU solve, I was<br>
only getting 11-12 digits of accuracy. On my other linear systems, I would<br>
get 15-16 digits of accuracy. So the system is poorly conditioned. However,<br>
the operator is computed from fields that are highly non-uniform over the<br>
spatial domain of the problem. The problem is a formation problem for a<br>
Field Reversed Configuration (FRC) magnetic fusion problem.<br></blockquote><div><br></div><div>No, this does not really help. For instance, from a Stokes problem we would</div><div>expect to hear, I am solving a Laplacian in the momentum equation, and</div>
<div>enforcing a differential constraint with the divergence. That allows us to</div><div>understand and build a preconditioner.</div><div><br></div><div>There is a simply staggeringly large literature on that problem. The time has</div>
<div>come for a literature search.</div><div><br></div><div>It is very likely that sacusp was ineffective because you have a vector problem,</div><div>but did not indicate this separation to the solver.</div><div><br></div>
<div> Matt</div><div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex;">
The matrix is symmetric but I am not sure if it is positive definite. I can<br>
check symmetry easily with petsc but am not sure if there is an easy way to<br>
tell if the matrix is positive definite - if there is a good way to do this<br>
with petsc, I'd love to hear about it. I did not find anything in the online<br>
help. I wonder if petsc might be able to say something about the positive<br>
definiteness of the system after a solution has been found. Anyway, I'm not<br>
up on the various ways to tell if a matrix is positive definite but I suppose<br>
it would be expensive to solve for all the eigenvalues of the matrix just to<br>
tell if it was positive definite. If there are clever and inexpensive ways<br>
to do this, I'd be interested in trying them.<br>
<br>
However, I did try both minres and symmlq from petsc using jacobi<br>
preconditioning and cuda. The case where I used minres ran to completion and<br>
used less than half the iterations that cg used so I was encouraged about<br>
that. However, there was an issue with the final solution that I need to<br>
investigate. The run using symmlq from petsc with jacobi preconditioning and<br>
cuda had numerical problems that I have not investigated and did not run to<br>
completion.<br>
<br>
So, because of the physics involved, the linear system is difficult to<br>
solve. It seems poorly conditioned. It is symmetric but might not be<br>
positive definite. The matrix comes from discretizing a pde on a regular<br>
structured 2d mesh.<br>
<br>
Does any of this additional info help?<br>
<br>
Thanks,<br>
<br>
Dave<br>
<br>
> Matt<br>
><br>
> > I'm wondering if there are suggestions of other preconditioners in petsc<br>
> > that 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<br>
> > interested 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<br>
> > at Tech-X and EMPhotonics. So I'd be happy to provide the matrix<br>
> > examples for 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>
</blockquote></div><br><br clear="all"><div><br></div>-- <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>