[petsc-users] Neumann BC with non-symmetric matrix

Thu Feb 25 16:18:34 CST 2016

Barry,

Thanks a lot for the detailed discussion. Things make much more sense now,
especially I was confused why the manual says to provide call both
'MatSetNullSpace' and 'MatSetTransposeNullSpace'. I have couple of more
questions and some observations I have made since yesterday.

1) Is there a systematic way to tell KSP to stop when it stagnates? I am
_not_ using SNES.

2) Once KSP converges to the least-square solution, the residual must be in
the nullspace of A^T because otherwise it would have been reduced by the
solver. So is it (at lest theoretically) possible to use the residual
vector as an approximate basis for the n(A^T)? In general this wouldn't be
enough but I'm thinking since the nullspace is one-dimensional, maybe I
could use the residual itself to manually project solution onto range of A
after calling KSPSolve?

3) Are preconditoners applied to the left by default? If not which one are
and which one are not?

4) So then if I provide the nullspace of A, KSP residual should converge,
correct? I have made a few observations:

   4-1) When I provide the nullspace of A through MatSetNullSpace, I
(generally) do see the residual (preconditioned) become very small (~1e-12
or so) but then it sometimes stagnates (say around 1e-10). Is this normal?

   4-2) Another observation is that the true residual stagnates way earlier
which I assume is an indication that the RHS is in fact not in the range of
A.

   4-3) Also, I've seen that the choice of KSP and PC have considerable
effects on the solver. For instance, by default I use hypre+bcgs and I have
noticed I need to change coarse-grid relaxation from Gaussian Elimination
to symmetric-sor/Jacobi otherwise hypre has issues. I assume this is
because the AMG preconditioner is also singular on the coarse level?

   4-4) Along the same lines, I tried a couple of other PCs such as
{jacobi, sor, gamg, ilu} and none of them were able to converge with bcgs
as the KSP. However, with gmres, almost all of them converge with the
exception of gamg.

   4-5) If I use gmres instead of bcgs, and for any of {jacobi, sor, ilu},
the true residual seems to be generally 2-3 orders of magnitude smaller
(1e-6 vs 1e-3). I suppose this is just because gmres is more robust?

   4-6) With hypre, the true residual is always larger (~1e-3) than other
PCs no matter if I use bcgs or gmres.

Sorry for the long discussion but this has turned out to be quite
educational for me!

Thanks,
Mohammad

On Wed, Feb 24, 2016 at 2:21 PM, Barry Smith <bsmith at mcs.anl.gov> wrote:

>
> > On Feb 24, 2016, at 9:07 AM, Mohammad Mirzadeh <mirzadeh at gmail.com>
> wrote:
> >
> > Barry,
> >
> > On Wednesday, February 24, 2016, Barry Smith <bsmith at mcs.anl.gov> wrote:
> >
> > > On Feb 24, 2016, at 12:07 AM, Mohammad Mirzadeh <mirzadeh at gmail.com>
> wrote:
> >
> >
> > At the PDE level the compatibility condition is satisfied. However I
> suspect that at the discrete level the rhs is not not exactly in the range.
> The reason for my suspicion is that I know for a fact that my
> discretization is not conservative at the hanging nodes.
>
>    Could be.
>
> >
> > Thanks for the suggestion. I'll give it a try. Howerver, does GMRES
> fundamentally behave differently than BiCGSTAB for these systems? I have
> seen people saying that it can keep the solution in the range if rhs is
> already in the range but I thought any KSP would do the same?
>
>
>     Here is the deal.  There are two issues here
>
> 1)  Making sure that b is the in the range of A.  If b is not in the range
> of A then it is obviously impossible to find an x such that A x = b. If we
> divide b into a part in the range of A called br and the rest call it bp
> then b = br + bp   and one can solve Ax = br   and the left over residual
> is bp. Normally if you run GMRES with an inconsistent right hand side (that
> is bp != 0) it will solve Ax = br automatically and thus give you a "least
> squares" answer which is likely what you want. These means in some sense
> you don't really need to worry about making sure that b is in the range of
> A. But the residuals of GMRES will stagnate, which makes sense because it
> cannot get rid of the bp part. In the least squares sense however GMRES has
> converged. If you provide MatSetTransposeNullSpace() then KSP automatically
> removes this space from b when it starts the Krylov method, this is nice
> because the GMRES residuals will go to zero.
>
> 2) Making sure that huge multiples of the null space of A do not get into
> the computed solution.
>
> With left preconditioning Krylov methods construct the solution in the
> space {Bb, BABb, BABABb, ..} if the range of B contains entries in the null
> space of A then the Krylov space will contain vectors in the null space of
> A. What can happen is that in each iteration of the Krylov space those
> entries grow and you end up with a solution x + xn where xn is in the null
> space of A and very large, thus it looks like GMRES is not converging since
> the solution "jump" each iteration. If you force the range of B to exclude
> the null space of A, that is project out the null space of A after applying
> B then nothing in the null space ever gets into the Krylov space and you
> get the "minimum norm" solution to the least squares problem which is what
> you almost always want. When you provide MatSetNullSpace() the KSP solvers
> automatically remove the null space after applying B.
>
>  All this stuff applies for any Krylov method.
>
>   So based on my understanding, you should just provide the null space
> that you can and forget out the transpose null space and use left
> preconditioning with GMRES (forget about what I said about trying with
> right preconditioning).  Let GMRES iterate until the residual norm has
> stabilized and use the resulting solution which is the least squares
> solution. If you are using KSP inside SNES you may need to set
> SNESSetMaxLinearSolveFailures() to a large value so it doesn't stop when it
> thinks the GMRES has failed.
>
>   Barry
>
> Matt and Jed, please check my logic I often flip my ranges/null spaces and
> may have incorrect presentation above.
> >
> >
> > > >
> > > > 2) I have tried fixing the solution at an arbitrary point, and while
> it generally works, for some problems I get numerical artifacts, e.g.
> slight asymmetry in the solution and/or increased error close to the point
> where I fix the solution. Is this, more or less, expected as a known
> artifact?
> >
> >   Yeah this approach is not good. We never recommend it.
> > > >
> > > > 3) An alternative to 2 is to enforce some global constraint on the
> solution, e.g. to require that the average be zero. My question here is
> two-fold:
> > >
> > >   Requiring the average be zero is exactly the same as providing a
> null space of the constant function. Saying the average is zero is the same
> as saying the solution is orthogonal to the constant function. I don't see
> any reason to introduce the Lagrange multiplier and all its complications
> inside of just providing the constant null space.
> > >
> > > Is this also true at the discrete level when the matrix is
> non-symmetric? I have always viewed this as just a constraint that could
> really be anything.
> > >
> > > >
> > > > 3-1) Is this generally any better than solution 2, in terms of not
> messing too much with the condition number of the matrix?
> > > >
> > > > 3-2) I don't quite know how to implement this using PETSc. Generally
> speaking I'd like to solve
> > > >
> > > > | A        U |   | X |   | B |
> > > > |            | * |   | = |   |
> > > > | U^T      0 |   | s |   | 0 |
> > > >
> > > >
> > > > where U is a constant vector (of ones) and s is effectively a
> Lagrange multiplier. I suspect I need to use MatCreateSchurComplement and
> pass that to the KSP? Or do I need create my own matrix type from scratch
> through MatCreateShell?
> > > >
> > > > Any help is appreciated!
> > > >
> > > > Thanks,
> > > > Mohammad
> > > >
> > > >
> > >
> > >
> > >
> > > --
> > > Sent from Gmail Mobile
> >
> >
> >
> > --
> > Sent from Gmail Mobile
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20160225/f8ecdf2b/attachment.html>