[petsc-users] move from KSPSetNullSpace to MatSetNullSpace

[Barry Smith]

   If we move the location of the nullspace used in removal from the pmat to the mat would that completely resolve the problem for you?

  Barry


  On Jun 23, 2015, at 7:42 AM, Stephan Kramer <s.kramer at imperial.ac.uk> wrote:
> 
>>> 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