August 2012 - petsc-dev - mailman3.cels.anl.gov

Re: [petsc-dev] FieldSplit fixes in 3.3
by Jed Brown 19 Aug '12

19 Aug '12

On Sun, Aug 19, 2012 at 12:48 PM, Dmitry Karpeev <karpeev(a)mcs.anl.gov>wrote: > 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). The current interface doesn't explicitly support all the ways to get information into the split, but several people have worked around the limitations to get the necessary information in there. Having you and Matt changing things to rely on mutually incompatible side-effects does not help. Maybe we can reach some agreement on the proper way to do things. More comments below. > > On Sun, Aug 19, 2012 at 9:32 AM, Jed Brown <jedbrown(a)mcs.anl.gov> wrote: > >> On Sun, Aug 19, 2012 at 2:31 AM, Dmitry Karpeev <karpeev(a)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. > Let's fix that "down the line". Setting a DM should never force it to be used or cause an error due to unsupported operation in a case where not having a DM is also acceptable. > 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. > You haven't explained why that was wrong. You are certainly setting up more than necessary here. >> >>> 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? > "force an inconsistent rhs" S may well be nonsingular. > > >> 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. > 1. I hate those hacky dm reference-non-references. 2. I don't want unused KSPs to be set up. (We should eventually fix the model so they don't exist...) 3. Two of the three IncrementTabLevels that you introduced are not necessary because the KSP's tab level was incremented before getting the PC. This unnecessary code is ugly and confusing. If we get rid of these three things, the patch reduces to fixing the tab level of the PC residing inside the MatSchurComplement, which I have no problem with.

1 0

Re: [petsc-dev] FieldSplit fixes in 3.3
by Dmitry Karpeev 19 Aug '12

19 Aug '12

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(a)mcs.anl.gov> wrote: > On Sun, Aug 19, 2012 at 2:31 AM, Dmitry Karpeev <karpeev(a)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.

1 0

Re: [petsc-dev] FieldSplit fixes in 3.3
by Jed Brown 19 Aug '12

19 Aug '12

On Sun, Aug 19, 2012 at 2:31 AM, Dmitry Karpeev <karpeev(a)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. - /* 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. - 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()? 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? > 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". 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.

1 0

Re: [petsc-dev] FieldSplit fixes in 3.3
by Matthew Knepley 19 Aug '12

19 Aug '12

On Sun, Aug 19, 2012 at 2:31 AM, Dmitry Karpeev <karpeev(a)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, > http://petsc.cs.iit.edu/petsc/releases/petsc-3.3/rev/5a9ebf885615) In > particular, the second changeset > modifies the way the splits' nullspaces are configured. In my view the > old behavior was not quite correct. > Here are the relevant comments: > http://petsc.cs.iit.edu/petsc/releases/petsc-3.3/rev/5a9ebf885615#l1.9 > and > http://petsc.cs.iit.edu/petsc/releases/petsc-3.3/rev/5a9ebf885615#l1.37 > > However, I don't have use cases where the old way of handling block > nullspaces > leads to a failure. If there is justifiable legitimate usage for the old > ways, the > changeset can be easily backed out. > 1) I fixed the Schur prefixes in dev, so do not merge this over it 2) It was my understanding that the Schur nullspace should be attached to A_11, but if we want to find another way to do it, that is fine. Matt > Dmitry. > -- 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

1 0

Re: [petsc-dev] removed PLAPACK interface Removed generic MatGetArray/MatRestoreArray()
by Jed Brown 19 Aug '12

19 Aug '12

On Sun, Aug 19, 2012 at 6:22 AM, Alexander Grayver <agrayver(a)gfz-potsdam.de>wrote: > I actually used them to implement Hadamard product since there is no > MatMatMultPointwise. > Is there now a different way to do that? > There is MatSeqAIJGetArray() and MatDenseGetArray().

1 0

Re: [petsc-dev] removed PLAPACK interface Removed generic MatGetArray/MatRestoreArray()
by Alexander Grayver 19 Aug '12

19 Aug '12

On 19.08.2012 04:43, Barry Smith wrote: > removed PLAPACK interface > Removed generic MatGetArray/MatRestoreArray() I actually used them to implement Hadamard product since there is no MatMatMultPointwise. Is there now a different way to do that? -- Regards, Alexander

1 0

FieldSplit fixes in 3.3
by Dmitry Karpeev 19 Aug '12

19 Aug '12

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, http://petsc.cs.iit.edu/petsc/releases/petsc-3.3/rev/5a9ebf885615) In particular, the second changeset modifies the way the splits' nullspaces are configured. In my view the old behavior was not quite correct. Here are the relevant comments: http://petsc.cs.iit.edu/petsc/releases/petsc-3.3/rev/5a9ebf885615#l1.9 and http://petsc.cs.iit.edu/petsc/releases/petsc-3.3/rev/5a9ebf885615#l1.37 However, I don't have use cases where the old way of handling block nullspaces leads to a failure. If there is justifiable legitimate usage for the old ways, the changeset can be easily backed out. Dmitry.

1 0

removed PLAPACK interface Removed generic MatGetArray/MatRestoreArray()
by Barry Smith 18 Aug '12

18 Aug '12

removed PLAPACK interface Removed generic MatGetArray/MatRestoreArray()

1 0

Re: [petsc-dev] Ghost values in sieve
by Matthew Knepley 18 Aug '12

18 Aug '12

On Sat, Aug 18, 2012 at 4:11 PM, Chris Eldred <chris.eldred(a)gmail.com>wrote: > Yes I think that makes sense. How do you access the associated PetscSF > for the DMComplex? Is it possible to define multiple PetscSF's for a > given DMComplex, each associated with different PetscSections? > 1) DMGetDefaultSF() for the one for the default Section 2) Yes, but you need to create them yourself using DMCreateDefaultSF() > My current understanding of using DMComplex is (for the serial case): > > 1) Create a DMComplex object, set its cones/supports, etc. > A good guide is to look at complexcreate.c, where I create a few default meshes. > 2) Create PetscSections that describe the layout of data across the > DMComplex > 3) Set a PetscSection as the default section > 4) Use DMCreateLocalVector to create the Vec objects needed for the > computation > > Is this correct? > Yes > In the parallel case, what else needs to be done? In parallel, a) I call DMComplexDistribute(). The alternative is to create your mesh directly in parallel, but this does not work with any mesh generator I know, and you would have to create the PetscSF for points that this function makes automatically. b) Now the section you create is local c) You can use the local section and pointSF to create a global section, which knows about overlap d) You can use the local and global section to make the dofSF c) and d) are done automatically for the default section (data layout). Matt > > -Chris > > On Fri, Aug 17, 2012 at 1:24 PM, Matthew Knepley <knepley(a)gmail.com> > wrote: > > On Fri, Aug 17, 2012 at 1:21 PM, Chris Eldred <chris.eldred(a)gmail.com> > > wrote: > >> > >> I was wondering how ghost values were handled when creating a > >> distributed (global) vector from a specific PetscSection over a > >> DMComplex object- ie what defines the stencils used to figure out > >> which values need to be ghosted? What about local vs. global indexes? > > > > > > This answer will have to be long because the question is mixing levels in > > the interface. > > > > A PetscSection is only an offset structure, so it has no designated > fields > > for overlap information. What I do for DM is to create 2 sections: one > > describing the layout of the local vector, and one describing the global > > vector and ghosts. The global section can be created automatically given > > the local section and a PetscSF (successor to VecScatter). > > > > The PetscSF describes the parallel overlap, and also has scatters and > > reductions. > > This is created by the DMComplex instead of a scatter since it can handle > > arbitrary types, and has a better query interface. DMComplex creates a > > (non-overlapping) partition of some stratum of points, usually cells for > > FEM. > > Then the closure(p) U star(p) for each point is added to that partition, > and > > the dofs on any shared point are ghosted. > > > > This is all now done using the vanilla DM interface, so DMGlobalToLocal() > > and LocalToGlobal() and cousins work as normal. > > > > Does that make sense? > > > > Matt > > > >> > >> > >> -Chris > >> > >> -- > >> Chris Eldred > >> DOE Computational Science Graduate Fellow > >> Graduate Student, Atmospheric Science, Colorado State University > >> B.S. Applied Computational Physics, Carnegie Mellon University, 2009 > >> chris.eldred(a)gmail.com > > > > > > > > > > -- > > 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 > > > > -- > Chris Eldred > DOE Computational Science Graduate Fellow > Graduate Student, Atmospheric Science, Colorado State University > B.S. Applied Computational Physics, Carnegie Mellon University, 2009 > chris.eldred(a)gmail.com > -- 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

1 0

Re: [petsc-dev] Ghost values in sieve
by Chris Eldred 18 Aug '12

18 Aug '12

Yes I think that makes sense. How do you access the associated PetscSF for the DMComplex? Is it possible to define multiple PetscSF's for a given DMComplex, each associated with different PetscSections? My current understanding of using DMComplex is (for the serial case): 1) Create a DMComplex object, set its cones/supports, etc. 2) Create PetscSections that describe the layout of data across the DMComplex 3) Set a PetscSection as the default section 4) Use DMCreateLocalVector to create the Vec objects needed for the computation Is this correct? In the parallel case, what else needs to be done? -Chris On Fri, Aug 17, 2012 at 1:24 PM, Matthew Knepley <knepley(a)gmail.com> wrote: > On Fri, Aug 17, 2012 at 1:21 PM, Chris Eldred <chris.eldred(a)gmail.com> > wrote: >> >> I was wondering how ghost values were handled when creating a >> distributed (global) vector from a specific PetscSection over a >> DMComplex object- ie what defines the stencils used to figure out >> which values need to be ghosted? What about local vs. global indexes? > > > This answer will have to be long because the question is mixing levels in > the interface. > > A PetscSection is only an offset structure, so it has no designated fields > for overlap information. What I do for DM is to create 2 sections: one > describing the layout of the local vector, and one describing the global > vector and ghosts. The global section can be created automatically given > the local section and a PetscSF (successor to VecScatter). > > The PetscSF describes the parallel overlap, and also has scatters and > reductions. > This is created by the DMComplex instead of a scatter since it can handle > arbitrary types, and has a better query interface. DMComplex creates a > (non-overlapping) partition of some stratum of points, usually cells for > FEM. > Then the closure(p) U star(p) for each point is added to that partition, and > the dofs on any shared point are ghosted. > > This is all now done using the vanilla DM interface, so DMGlobalToLocal() > and LocalToGlobal() and cousins work as normal. > > Does that make sense? > > Matt > >> >> >> -Chris >> >> -- >> Chris Eldred >> DOE Computational Science Graduate Fellow >> Graduate Student, Atmospheric Science, Colorado State University >> B.S. Applied Computational Physics, Carnegie Mellon University, 2009 >> chris.eldred(a)gmail.com > > > > > -- > 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 -- Chris Eldred DOE Computational Science Graduate Fellow Graduate Student, Atmospheric Science, Colorado State University B.S. Applied Computational Physics, Carnegie Mellon University, 2009 chris.eldred(a)gmail.com

1 0