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
