<div dir="ltr"><div class="gmail_extra"><div class="gmail_quote">On Tue, Jun 23, 2015 at 7:42 AM, Stephan Kramer <span dir="ltr"><<a href="mailto:s.kramer@imperial.ac.uk" target="_blank">s.kramer@imperial.ac.uk</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
On Mon, Jun 22, 2015 at 1:13 PM, Stephan Kramer <s.kramer at<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
<a href="http://imperial.ac.uk" rel="noreferrer" target="_blank">imperial.ac.uk</a>><br>
wrote:<br>
<br>
 Dear petsc devs<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
<br>
I've been trying to move our code from using KSPSetNullSpace to use<br>
MatSetNullSpace instead. Although I appreciate the convenience of the<br>
nullspace being propagated automatically through the solver hierarchy,<br>
I'm<br>
still a little confused on how to deal with the case that mat/=pmat in a<br>
ksp.<br>
<br>
If I read the code correctly I need to call MatSetNullSpace on the pmat<br>
in<br>
order for a ksp to project that nullspace out of the residual during the<br>
krylov iteration. However the nullspaces of mat and pmat are not<br>
necessarily the same. For instance, what should I do if the system that I<br>
want to solve has a nullspace (i.e. the `mat` has a null vector), but the<br>
preconditioner matrix does not.<br>
<br>
As an example of existing setups that we have for which this is the case:<br>
take a standard Stokes velocity, pressure system - where we want to solve<br>
for pressure through a Schur complement. Suppose the boundary conditions<br>
are Dirichlet for velocity on all boundaries, then the pressure equation<br>
has the standard, constant nullspace. A frequently used preconditioner is<br>
the "scaled pressure mass matrix" where G^T K^{-1} G is approximated by a<br>
pressure mass matrix that is scaled by the value of the viscosity. So in<br>
this case the actual system has a nullspace, but the preconditioner<br>
matrix<br>
does not.<br>
<br>
The way we previously used to solve this is by setting the nullspace on<br>
the ksp of the Schur system, where the mat comes from<br>
MatCreateSchurComplement and the pmat is the scaled mass matrix. We then<br>
set the pctype to PCKSP and do not set a nullspace on its associated<br>
ksp. I<br>
don't really see how I can do this using only MatSetNullSpace - unless I<br>
create a copy of the mass matrix and have one copy with and one copy<br>
without a nullspace so that I could use the one with the nullspace for<br>
the<br>
pmat of the Schur ksp (and thus the mat of the pcksp) and the copy<br>
without<br>
the nullspace as the pmat for the pcksp.<br>
<br>
We would very much appreciate some guidance on what the correct way to<br>
deal with this kind of situation is<br>
<br>
<br>
</blockquote>
I can understand that this is inconsistent (we have not really thought of<br>
a<br>
good way to make this consistent). However, does<br>
this produce the wrong results? If A has a null space, won't you get the<br>
same answer by attaching it to P, or am I missing the<br>
import of the above example?<br>
<br>
</blockquote>
<br>
Attaching the nullspace to P means it will be applied in the PCKSP solve as<br>
well - which in this case is just a mass matrix solve which doesn't<br>
have a nullspace. I was under the impression that removing a nullspace in<br>
a solve<br>
where the matrix doesn't having one, would lead to trouble - but I'm happy<br>
to<br>
be told wrong.<br>
<br>
</blockquote>
<br>
I guess I was saying that this solve is only applied after a system with<br>
that null space,<br>
so I did not think it would change the answer.<br>
<br>
  Thanks,<br>
<br>
     Matt<br>
</blockquote>
<br>
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:<br>
<br>
Say the original PCKSP solve (that doesn't apply the nullspace) gets passed some residual r from the outer solve, so it solves:<br>
<br>
  Pz=r<br>
<br>
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:<br>
<br>
  NM^{-1}P z = NM^{-1}r<br>
<br>
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:<br>
<br>
  z -> z + P^{-1}M n<br>
<br>
for arbitrary n in the nullspace.<br>
<br>
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.<br></blockquote><div><br></div><div>This analysis does not make sense to me. The whole point of having all this nullspace stuff is so that we do not get</div><div>a component of the null space in the solution. The way we avoid these in Krylov methods is to project them out, </div><div>which is what happens when you attach it. Thus, your PCKSP solution 'z' will not have any 'n' component. Your</div><div>outer solve will also not have any 'n' component, so I do not see where the inconsistency comes in.</div><div><br></div><div>  Thanks,</div><div><br></div><div>    Matt</div><div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
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.<br>
<br>
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.<br>
<br>
Cheers<span class="HOEnZb"><font color="#888888"><br>
Stephan<br>
</font></span></blockquote></div><br><br clear="all"><div><br></div>-- <br><div class="gmail_signature">What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br>-- Norbert Wiener</div>
</div></div>