[petsc-users] move from KSPSetNullSpace to MatSetNullSpace

[Stephan Kramer]

>> On Mon, Jun 22, 2015 at 1:13 PM, Stephan Kramer <s.kramer at
>>> imperial.ac.uk>
>>> wrote:
>>>
>>>  Dear petsc devs
>>>>
>>>> I've been trying to move our code from using KSPSetNullSpace to use
>>>> MatSetNullSpace instead. Although I appreciate the convenience of the
>>>> nullspace being propagated automatically through the solver hierarchy,
>>>> I'm
>>>> still a little confused on how to deal with the case that mat/=pmat in a
>>>> ksp.
>>>>
>>>> If I read the code correctly I need to call MatSetNullSpace on the pmat
>>>> in
>>>> order for a ksp to project that nullspace out of the residual during the
>>>> krylov iteration. However the nullspaces of mat and pmat are not
>>>> necessarily the same. For instance, what should I do if the system that I
>>>> want to solve has a nullspace (i.e. the `mat` has a null vector), but the
>>>> preconditioner matrix does not.
>>>>
>>>> As an example of existing setups that we have for which this is the case:
>>>> take a standard Stokes velocity, pressure system - where we want to solve
>>>> for pressure through a Schur complement. Suppose the boundary conditions
>>>> are Dirichlet for velocity on all boundaries, then the pressure equation
>>>> has the standard, constant nullspace. A frequently used preconditioner is
>>>> the "scaled pressure mass matrix" where G^T K^{-1} G is approximated by a
>>>> pressure mass matrix that is scaled by the value of the viscosity. So in
>>>> this case the actual system has a nullspace, but the preconditioner
>>>> matrix
>>>> does not.
>>>>
>>>> The way we previously used to solve this is by setting the nullspace on
>>>> the ksp of the Schur system, where the mat comes from
>>>> MatCreateSchurComplement and the pmat is the scaled mass matrix. We then
>>>> set the pctype to PCKSP and do not set a nullspace on its associated
>>>> ksp. I
>>>> don't really see how I can do this using only MatSetNullSpace - unless I
>>>> create a copy of the mass matrix and have one copy with and one copy
>>>> without a nullspace so that I could use the one with the nullspace for
>>>> the
>>>> pmat of the Schur ksp (and thus the mat of the pcksp) and the copy
>>>> without
>>>> the nullspace as the pmat for the pcksp.
>>>>
>>>> We would very much appreciate some guidance on what the correct way to
>>>> deal with this kind of situation is
>>>>
>>>>
>>> I can understand that this is inconsistent (we have not really thought of
>>> a
>>> good way to make this consistent). However, does
>>> this produce the wrong results? If A has a null space, won't you get the
>>> same answer by attaching it to P, or am I missing the
>>> import of the above example?
>>>
>>
>> Attaching the nullspace to P means it will be applied in the PCKSP solve as
>> well - which in this case is just a mass matrix solve which doesn't
>> have a nullspace. I was under the impression that removing a nullspace in
>> a solve
>> where the matrix doesn't having one, would lead to trouble - but I'm happy
>> to
>> be told wrong.
>>
>
> I guess I was saying that this solve is only applied after a system with
> that null space,
> so I did not think it would change the answer.
>
>   Thanks,
>
>      Matt

Sorry if we're talking completely cross-purposes here: but the pcksp 
solve (that doesn't actually have a nullspace) is inside a solve that 
does have a nullspace. If what you are saying is that applying the 
nullspace inside the pcksp solve does not affect the outer solve, I can 
only see that if the difference in outcome of the pcksp solve (with or 
without the nullspace) would be inside the nullspace, because such a 
difference would be projected out anyway straight after the pcksp solve. 
I don't see that this is true however. My reasoning is the following:

Say the original PCKSP solve (that doesn't apply the nullspace) gets 
passed some residual r from the outer solve, so it solves:

   Pz=r

where P is our mass matrix approximation to the system in the outer 
solve, and P is full rank. Now if we do remove the nullspace during the 
pcksp solve, effectively we change the system to the following:

   NM^{-1}P z = NM^{-1}r

where M^{-1} is the preconditioner inside the pcksp solve (assuming 
left-preconditiong here), and N is the projection operator that projects 
out the nullspace. This system is now rank-deficient, where I can add:

   z -> z + P^{-1}M n

for arbitrary n in the nullspace.

So not only is the possible difference between solving with or without a 
nullspace not found in that nullspace, but worse, I've ended up with a 
preconditioner that's rank deficient in the outer solve.

I've experimented with the example I described before with the Schur 
complement of a Stokes velocity,pressure system on the outside using 
fgmres and the viscosit scaled pressure mass matrix as preconditioner 
which is solved in the pcksp solve with cg+sor. With petsc 3.5 (so I can 
still use KSPSetNullSpace) if I set a nullspace only on the outside ksp, 
or if I set it both on the outside ksp and the pcksp sub_ksp, I get 
different answers (and also a different a number of iterations). The 
difference in answer visually looks like some grid scale noise that 
indeed could very well be of the form of a constant vector multiplied by 
an inverse mass matrix.

If you agree that applying the nullspace inside the pcksp indeed gives 
incorrect answers, I'm still not clear how I can set this up using 
MatSetNullSpace only.

Cheers
Stephan