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 <[email protected]> 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
