[petsc-users] move from KSPSetNullSpace to MatSetNullSpace

[Matthew Knepley]

On Tue, 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.
>

This analysis does not make sense to me. The whole point of having all this
nullspace stuff is so that we do not get
a component of the null space in the solution. The way we avoid these in
Krylov methods is to project them out,
which is what happens when you attach it. Thus, your PCKSP solution 'z'
will not have any 'n' component. Your
outer solve will also not have any 'n' component, so I do not see where the
inconsistency comes in.

  Thanks,

    Matt


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



-- 
What most experimenters take for granted before they begin their
experiments is infinitely more interesting than any results to which their
experiments lead.
-- Norbert Wiener
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20150623/2fc14469/attachment.html>