[petsc-dev] FieldSplit fixes in 3.3
Dmitry Karpeev
karpeev at mcs.anl.gov
Sun Aug 19 12:48:41 CDT 2012
The main reason for this being a patch in 3.3 is that recursive FieldSplit
is broken there (which is why I couldn't enable it in libMesh/Moose
correctly). More comments below.
On Sun, Aug 19, 2012 at 9:32 AM, Jed Brown <jedbrown at mcs.anl.gov> wrote:
> On Sun, Aug 19, 2012 at 2:31 AM, Dmitry Karpeev <karpeev at mcs.anl.gov>wrote:
>
>> I just pushed a couple of fixes to PCFIELDSPLIT in petsc-3.3 (
>> http://petsc.cs.iit.edu/petsc/releases/petsc-3.3/rev/943dde820f7f,
>
>
> Yikes:
>
> --- a/src/ksp/pc/impls/fieldsplit/fieldsplit.c
> +++ b/src/ksp/pc/impls/fieldsplit/fieldsplit.c
> @@ -14,6 +14,7 @@
> PetscInt *fields,*fields_col;
> VecScatter sctx;
> IS is,is_col;
> + DM dm;
> PC_FieldSplitLink next,previous;
> };
>
> @@ -263,6 +264,13 @@
> if(dms) {
> ierr = PetscInfo(pc,"Setting up physics based fieldsplit
> preconditioner using the embedded DM\n");CHKERRQ(ierr);
> for (ilink=jac->head,i=0; ilink; ilink=ilink->next,i++) {
> + /*
> + HACK:
> + keep a handle to DM here without increfing it: may need the DM
> to set up the Schur solvers;
> + can't rely on KSPGetDM() since it will create a DMShell if
> none was set;
> + can't use ilink->ksp->dm without creating a dependence on
> dmimpl.h.
> + */
> + ilink->dm = dms[i];
>
> I think your comment was intended to say "kspimpl.h", but you are creating
> a tangled web of barbed wire here. You are supposed to call KSPGetDM() and
> switch based on its *capability*, not *presence*. PETSc should never do
> something different based on existence of a DM.
>
I need to pass the DM from the A00 split to the Schur solvers, but not if
it's null. If it was PETSC_NULL in the A00 split,
KSPGetDM() will give me a DMShell, which will generate an error down the
line. I can't reliably get the DM from the split's KSP,
so I have to stash it somewhere. Otherwise recursive splits don't work. If
there is a better solution, I'm all for it.
>
> - /* set tabbing and options prefix of KSP inside the MatSchur */
> + /* set tabbing, options prefix and DM of KSP inside the MatSchur */
> ierr = MatSchurComplementGetKSP(jac->schur,&ksp);CHKERRQ(ierr);
> ierr =
> PetscObjectIncrementTabLevel((PetscObject)ksp,(PetscObject)pc,2);CHKERRQ(ierr);
> + {
> + PC pcinner;
> + ierr = KSPGetPC(ksp, &pcinner); CHKERRQ(ierr);
> + ierr =
> PetscObjectIncrementTabLevel((PetscObject)pcinner,(PetscObject)pc,2);CHKERRQ(ierr);
> + }
>
> There is KSPIncrementTabLevel() for a reason.
>
In petsc-dev, but not in petsc-3.3.
>
> - ierr = PetscSNPrintf(schurprefix,sizeof
> schurprefix,"%sfieldsplit_%s_",((PetscObject)pc)->prefix?((PetscObject)pc)->prefix:"",ilink->splitname);CHKERRQ(ierr);
> - ierr =
> KSPSetOptionsPrefix(jac->kspschur,schurprefix);CHKERRQ(ierr);
> /* really want setfromoptions called in
> PCSetFromOptions_FieldSplit(), but it is not ready yet */
> /* need to call this every time, since the jac->kspschur is freshly
> created, otherwise its options never get set */
> - ierr = KSPSetFromOptions(jac->kspschur);CHKERRQ(ierr);
> + ierr = PetscSNPrintf(schurprefix,sizeof
> schurprefix,"%sfieldsplit_%s_",((PetscObject)pc)->prefix?((PetscObject)pc)->prefix:"",ilink->splitname);CHKERRQ(ierr);
> + ierr = KSPSetOptionsPrefix(jac->kspschur,schurprefix);CHKERRQ(ierr);
> + ierr = KSPSetFromOptions(jac->kspschur); CHKERRQ(ierr);
> + ierr = KSPSetUp(jac->kspschur); CHKERRQ(ierr);
>
> Why this KSPSetUp()?
>
I guess this KSPSetUp() can go.
> You have changed the formatting here, but not explained why. (Also, can't
> you follow PETSc coding guidelines regarding placement of CHKERRQ?)
>
> - } else {
> - /* set up the individual PCs */
> - i = 0;
> - ilink = jac->head;
> - while (ilink) {
> - ierr =
> KSPSetOperators(ilink->ksp,jac->mat[i],jac->pmat[i],flag);CHKERRQ(ierr);
> - /* really want setfromoptions called in
> PCSetFromOptions_FieldSplit(), but it is not ready yet */
> - if (!jac->suboptionsset) {ierr =
> KSPSetFromOptions(ilink->ksp);CHKERRQ(ierr);}
> - ierr = KSPSetUp(ilink->ksp);CHKERRQ(ierr);
> - i++;
> - ilink = ilink->next;
> - }
> + }
> + /* set up the individual PCs */
> + i = 0;
> + ilink = jac->head;
> + while (ilink) {
> + ierr = KSPSetOperators(ilink->ksp,jac->mat[i],jac->pmat[i],flag);
> CHKERRQ(ierr);
> + /* really want setfromoptions called in
> PCSetFromOptions_FieldSplit(), but it is not ready yet */
> + if (!jac->suboptionsset) {ierr = KSPSetFromOptions(ilink->ksp);
> CHKERRQ(ierr);}
> + ierr = KSPSetUp(ilink->ksp);CHKERRQ(ierr);
> + i++;
> + ilink = ilink->next;
>
> Why do you want this stuff done even for Schur? At least one of these
> blocks is not used at all, so why are we setting it up?
>
In the previous version NONE of these were done for Schur. I suppose we
can handle these cases more carefully.
>
>
>> http://petsc.cs.iit.edu/petsc/releases/petsc-3.3/rev/5a9ebf885615)
>>
>
> /* HACK: Check for the constant null space */
> - ierr = MatGetNullSpace(pc->pmat, &sp);CHKERRQ(ierr);
> + /*
> + This replaces a previous test that checked whether the split's
> constants are injected into the ambient problem's nullspace.
> + That can create false negatives (declaring constants not in the
> nullspace, when they are) as in Stokes or A = [1, 1; 1, 0].
> + */
> if (sp) {
>
> So now you have stale comments and conditionals based on uninitialized
> memory (sp)? I suppose your intent is to always do this constant null space
> check regardless of whether the larger matrix has a null space set. I agree
> that the old code was incorrect because the restriction of the null space
> of the parent matrix is not comparable to the null space of the submatrix.
>
> ierr = MatGetNullSpace(jac->pmat[1], &sp);CHKERRQ(ierr);
> - if (sp) {ierr = MatSetNullSpace(jac->schur, sp);CHKERRQ(ierr);}
> - /* set tabbing, options prefix and DM of KSP inside the MatSchur */
> + /*
> + Do we really want to attach the A11-block's nullspace to S? Note
> that this may generate
> + "false positives", when the A11's nullspace isn't S's: Stokes or A
> = [1, 1; 1, 0].
> + Using a false nullspace may prevent KSP(S) from converging, since
> it might force an inconsistent rhs.
> + */
> + /* if (sp) {ierr = MatSetNullSpace(jac->schur, sp);CHKERRQ(ierr);}
> */
> + /* set tabbing, options prefix and DM of KSP inside the
> MatSchurComplement */
>
> I agree with this change, though the comment is incorrect. The typical
> problem is that S is nonsingular or has a small null space while A11 is has
> huge null space (e.g. the zero matrix). The true null space projection for
> the zero matrix prevents progress in any direction.
>
>
> I think we can repair the second patch, though we really need a general
> mechanism for providing Schur complement null spaces. Matt advocates
> attaching it to A11 (where the user knows full well that it's not the null
> space of A11 when they set it), making the old code "correct".
>
Okay, I'll back out that part of the second patch, since A11's nullspace is
assumed to be meant for S. I'm not sure what's incorrect about the
comment, though: if a vector is in the A11's kernel, but not in S's, I call
it a "false positive". This terminology may be wanting, but what's
incorrect about it?
> I think that's confusing and we need a specific API for it, but if people
> are using it that way, we shouldn't change it in 3.3.
>
> I'd prefer to backout the first patch, split it into separate pieces, and
> review each before pushing. I think it's a potentially big behavior change
> for 3.3.
>
Which pieces of the first patch do you think are too big for 3.3? The
splits' and Schur KSPs have to be set up (maybe not all in every situation)
in order for recursive splitting to work. Should we declare that
capability unavailable for 3.3? I can eliminate the unneeded KSPSetUp()
calls.
Dmitry.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-dev/attachments/20120819/4ba594a8/attachment.html>
More information about the petsc-dev
mailing list