Barry,<br><br>On Wednesday, February 24, 2016, Barry Smith <<a href="mailto:bsmith@mcs.anl.gov">bsmith@mcs.anl.gov</a>> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><br>
> On Feb 24, 2016, at 12:07 AM, Mohammad Mirzadeh <<a href="javascript:;" onclick="_e(event, 'cvml', 'mirzadeh@gmail.com')">mirzadeh@gmail.com</a>> wrote:<br>
><br>
><br>
> Barry,<br>
> On Wednesday, February 24, 2016, Barry Smith <<a href="javascript:;" onclick="_e(event, 'cvml', 'bsmith@mcs.anl.gov')">bsmith@mcs.anl.gov</a>> wrote:<br>
><br>
> > On Feb 23, 2016, at 11:35 PM, Mohammad Mirzadeh <<a href="javascript:;" onclick="_e(event, 'cvml', 'mirzadeh@gmail.com')">mirzadeh@gmail.com</a>> wrote:<br>
> ><br>
> > Dear all,<br>
> ><br>
> > I am dealing with a situation I was hoping to get some suggestions here. Suppose after discretizing a poisson equation with purely neumann (or periodic) bc I end up with a matrix that is *almost* symmetric, i.e. it is symmetric for almost all grid points with the exception of a few points.<br>
><br>
> How come it is not purely symmetric? The usual finite elements with pure Neumann or periodic bc will give a completely symmetric matrix.<br>
><br>
> Barry<br>
><br>
><br>
> So this is a finite difference discretization on adaptive Cartesian grids. It turns out that the discretization is non-symmetric at the corse-fine interface. It's actually not because of the BC itself.<br>
<br>
Oh yes, I remember that issue vaguely now.<br>
<br>
> ><br>
> > The correct way of handling this problem is by specifying the nullspace to MatSetNullSpace. However, since the matrix is non-symmetric in general I would need to pass the nullspace of A^T. Now it turns out that if A is *sufficiently close to being symmetric*, I can get away with the constant vector, which is the nullspace of A and not A^T, but obviously this does not always work. Sometimes the KSP converges and in other situations the residual stagnates which is to be expected.<br>
> ><br>
> > Now, here are my questions (sorry if they are too many!):<br>
> ><br>
> > 1) Is there any efficient way of calculating nullspace of A^T in this case?<br>
<br>
Not that I know of.<br>
<br>
> Is SVD the only way?<br>
<br>
I think if you make sure you only use right preconditioning such as with -ksp_type gmres -ksp_pc_side right AND you know that the right hand side is automatically in the range of A then the null space of A^T is never needed in the computation. Is your right hand side in the range of A?<br>
<br>
Barry</blockquote><div><br></div><div>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. </div><div><br></div><div>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?<span></span></div><div><br></div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
<br>
> ><br>
> > 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?<br>
<br>
Yeah this approach is not good. We never recommend it.<br>
> ><br>
> > 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:<br>
><br>
> 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.<br>
><br>
> 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.<br>
><br>
> ><br>
> > 3-1) Is this generally any better than solution 2, in terms of not messing too much with the condition number of the matrix?<br>
> ><br>
> > 3-2) I don't quite know how to implement this using PETSc. Generally speaking I'd like to solve<br>
> ><br>
> > | A U | | X | | B |<br>
> > | | * | | = | |<br>
> > | U^T 0 | | s | | 0 |<br>
> ><br>
> ><br>
> > 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?<br>
> ><br>
> > Any help is appreciated!<br>
> ><br>
> > Thanks,<br>
> > Mohammad<br>
> ><br>
> ><br>
><br>
><br>
><br>
> --<br>
> Sent from Gmail Mobile<br>
<br>
</blockquote><br><br>-- <br>Sent from Gmail Mobile<br>