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

[Barry Smith]

> On Feb 23, 2016, at 11:35 PM, Mohammad Mirzadeh <mirzadeh at gmail.com> wrote:
> 
> Dear all,
> 
> 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.

  How come it is not purely symmetric? The usual finite elements with pure Neumann or periodic bc will give a completely symmetric matrix.

  Barry

> 
> 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.
> 
> Now, here are my questions (sorry if they are too many!):
>  
> 1) Is there any efficient way of calculating nullspace of A^T in this case? Is SVD the only way?
> 
> 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?
> 
> 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.

> 
> 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
> 
>