We are changing the code to always use a random right hand side for computing the Chebyshev estimates instead of the given right hand side so hopefully this problem will go away in the near future.
Barry > On Sep 2, 2015, at 1:52 AM, Filippo Leonardi <[email protected]> wrote: > > So I spent a lot of time debugging my code for this one, it turns out that in > my debug run, where I solve the same system multiple times, reusing the > preconditioner if possible, sometimes, it happens that the rhs is 0. It has > nothing to to with scaling. Just that when I did my debug run I just had > reworked the code to use nonuniform meshes, so I assumed incorrectly, hey, > there is an issue with scaling. > > I do not perform the checks for incompatible rhs in this case, despite > knowing that PETSc breaks in this case (I knew that, but I forgot). I now > introduced a check in my debug build that check if a rhs is zero or not. > > Long story short: > I have the suspect that PETSc could be improved in this respect: I am not an > expert but somehow it shouldn't break or should warn the user (even if he is > doing somthing "nonsensical" like solving Ax = 0). Could it be that if I set > PETSc to solve a system with rhs 0 it breaks (maybe if combined with a > multiple solve?)? What I do is the following: > > i approximate the Leray projection operator with Laplace, so i solve Ax = b > multiple times, reusing the matrix/precond. When I debug the code, sometimes > I have a b = 0. > > (In fact, this was my first misktake: > http://lists.mcs.anl.gov/pipermail/petsc-users/2014-November/023478.html > ) > > Maybe because the relative error cannot possibly work? If so, can I suggest > to insert some sort of check for idiots like me :) (I repeated this mistake > twice knowing what the issue was)? > > On Tuesday 01 September 2015 12:13:34 Matthew Knepley wrote: > > On Tue, Sep 1, 2015 at 6:53 AM, Filippo Leonardi <[email protected]> > > > > wrote: > > > On Tuesday 01 September 2015 06:32:38 you wrote: > > > > > > 1) The KSP view does not say it is shifting. Are you using the > > > > > > latest > > > > > > > > > > > > > > > > > > > > > > > > release? > > > > > > > > > > yes, 3.6. Does PETSc warn for that even if I set the nullspace? I can > > > > > > also > > > > > > > > check MUMPS or something else. > > > > > > > > I am not sure what you think PETSc does here. If shifting were enabled, > > > > > > it > > > > > > > would add some > > > > > > > > diagonal matrix to the input matrix and continue the factorization. This > > > > > > > > would mean that the > > > > > > > > factors were not, in fact, a factorization of the input matrix, and you > > > > > > > > would not get the exact > > > > > > > > solution in one iterate. > > > > > > I though PETSc would've replace my pivots with small eps, which is > > > actually not a problem in my case > > > > > > > > > 2) If it shifted, it would not solve in a single iterate. > > > > > > > > > > even with preonly? > > > > > > > > You would have a large residual. Do you? > > > > > > Actually, I get a perfect solution. > > > > > > > > > 3) Your GAMG results imply that something is wrong with the coarse > > > > > > > > > > > > solve. > > > > > > > > > > > > > > > > > > > > > > > > This is exactly what would happen if > > > > > > > > > > > > > > > > > > > > > > > > that problem was not solved accurately (its off by > 10 orders of > > > > > > > > > > > > > > > > > > > > > > > > magnitude). > > > > > > > > > > yes, but GAMG builds is own coarse solvers so either the problem is > > > > > > > > > > already in the definition of A and b (likely) or it is a bug in gamg. > > > > > > > > Yes. GAMG uses the constants to build the basis, on the assumption that > > > > > > > > they are in the (near) nullspace of the > > > > > > > > operator with no boundary conditions. Since this is far off, I think > > > > this > > > > > > > > must not be true for your A. > > > > > > > > > > It sounds like your operator is not singular, and its not the > > > > > > Laplacian > > > > > > > > > since it does not look like the Neumann version > > > > > > > > > > > > > > > > > > > > > > > > has constants as a null space. > > > > > > > > > > I'm using periodic boundaries, and constants are in kern(A) > > > > > > > > Did you check? > > > > > > Checked with VecSet and MatMult just in case, I get a machine eps constant > > > vector. > > > > Okay, it seems that we have the best chance of figuring out the problem by > > looking at option 2) and 4). > > > > First, is your rhs consistent, meaning is it orthogonal to your nullspace? > > > > Matt > > > > > > Matt > > > > > > > > > Thanks a lot for your time. > > > > > > > > > > > Thanks,
