<br><br><div class="gmail_quote">On Sun, Aug 19, 2012 at 1:17 PM, Jed Brown <span dir="ltr"><<a href="mailto:jedbrown@mcs.anl.gov" target="_blank">jedbrown@mcs.anl.gov</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
<div class="gmail_quote"><div>On Sun, Aug 19, 2012 at 12:48 PM, Dmitry Karpeev <span dir="ltr"><<a href="mailto:karpeev@mcs.anl.gov" target="_blank">karpeev@mcs.anl.gov</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
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).</blockquote></div><div><br>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. </div>
</div></blockquote><div>I don't think the present fix is about the API or working through side effects (nullspace stuff excluded and partially backed out here: <a href="http://petsc.cs.iit.edu/petsc/releases/petsc-3.3/rev/1723d4624521">http://petsc.cs.iit.edu/petsc/releases/petsc-3.3/rev/1723d4624521</a>)</div>
<div>The current interface assumes that the A00 solver and the Schur inner solver are to be set up identically (something we relaxed in petsc-dev). In order to use fieldsplit on these solvers the DM for the corresponding split must be forwarded to both of these A00 solvers. It was being set only on the "outer" A00 solver. Partly this is an artifact of our "split" setup process: some of it occurs in PCFieldSplitSetDefaults(), some later in PCSetUp_FieldSplit(). Some objects obtained in PCFieldSplitSetDefaults() need to be cached until PCSetUp_FieldSplit(), and the splits' DMs are among those. I now put in proper reference counting for them here:</div>
<div><a href="http://petsc.cs.iit.edu/petsc/releases/petsc-3.3/rev/188af9799779">http://petsc.cs.iit.edu/petsc/releases/petsc-3.3/rev/188af9799779</a></div><div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
<div class="gmail_quote"><div>Maybe we can reach some agreement on the proper way to do things.</div>
</div></blockquote><div>That's fine with me. The reason I wanted to fix this particular problem now is that some dependent packages (e.g., Moose) only rely on release versions of petsc and would not be able to use this functionality properly for quite some time. As far as I see it, there are </div>
<div>other problems with fieldsplit that can (and should) wait for petsc-dev to be resolved (e.g., the proper handling PCReset()).</div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
<div class="gmail_quote"><div><div>
<div><br></div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"> More comments below.<br><br><div class="gmail_quote"><div>On Sun, Aug 19, 2012 at 9:32 AM, Jed Brown <span dir="ltr"><<a href="mailto:jedbrown@mcs.anl.gov" target="_blank">jedbrown@mcs.anl.gov</a>></span> wrote:<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div class="gmail_quote"><div>On Sun, Aug 19, 2012 at 2:31 AM, Dmitry Karpeev <span dir="ltr"><<a href="mailto:karpeev@mcs.anl.gov" target="_blank">karpeev@mcs.anl.gov</a>></span> wrote:<br>
</div><div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
I just pushed a couple of fixes to PCFIELDSPLIT in petsc-3.3 (<a href="http://petsc.cs.iit.edu/petsc/releases/petsc-3.3/rev/943dde820f7f" target="_blank">http://petsc.cs.iit.edu/petsc/releases/petsc-3.3/rev/943dde820f7f</a>, </blockquote>
<div><br></div></div><div>Yikes:</div><div><br></div><div><div><font face="courier new, monospace">--- a/src/ksp/pc/impls/fieldsplit/fieldsplit.c</font></div><div><font face="courier new, monospace">+++ b/src/ksp/pc/impls/fieldsplit/fieldsplit.c</font></div>
<div><font face="courier new, monospace">@@ -14,6 +14,7 @@</font></div><div><font face="courier new, monospace"> PetscInt *fields,*fields_col;</font></div><div><font face="courier new, monospace"> VecScatter sctx;</font></div>
<div><font face="courier new, monospace"> IS is,is_col;</font></div><div><font face="courier new, monospace">+ DM dm;</font></div><div><font face="courier new, monospace"> PC_FieldSplitLink next,previous;</font></div>
<div><font face="courier new, monospace"> };</font></div><div><font face="courier new, monospace"> </font></div><div><font face="courier new, monospace">@@ -263,6 +264,13 @@</font></div><div><font face="courier new, monospace"> if(dms) {</font></div>
<div><font face="courier new, monospace"> ierr = PetscInfo(pc,"Setting up physics based fieldsplit preconditioner using the embedded DM\n");CHKERRQ(ierr);</font></div><div><font face="courier new, monospace"> for (ilink=jac->head,i=0; ilink; ilink=ilink->next,i++) {</font></div>
<div><font face="courier new, monospace">+ /* </font></div><div><font face="courier new, monospace">+ HACK: </font></div><div><font face="courier new, monospace">+ keep a handle to DM here without increfing it: may need the DM to set up the Schur solvers; </font></div>
<div><font face="courier new, monospace">+ can't rely on KSPGetDM() since it will create a DMShell if none was set;</font></div><div><font face="courier new, monospace">+ can't use ilink->ksp->dm without creating a dependence on dmimpl.h.</font></div>
<div><font face="courier new, monospace">+ */</font></div><div><font face="courier new, monospace">+ ilink->dm = dms[i]; </font></div></div><div><br></div><div>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 <b>capability</b>, not <b>presence</b>. PETSc should never do something different based on existence of a DM.</div>
</div></blockquote></div><div>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, </div><div>KSPGetDM() will give me a DMShell, which will generate an error down the line. </div>
</div></blockquote><div><br></div></div></div><div>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.</div>
</div></blockquote><div>Currently the inner and outer A00 solvers are essentially identified (duplicated), so they should be set up identically, including the DM.</div><div>Especially when that DM defines a recursive split. This is in part because there is no way to configure those two differently programmatically or from the command line. So if the outer A00 is using -fieldsplit_0_pc_type_fieldsplit, the inner A00 will as well.</div>
<div>However, the DM on the inner A00 will be rather different (or absent) and produce splits incompatible with the rest of the options (e.g., a single default split, while -fieldsplit_0_pc_fieldsplit_type schur). </div>
<div><br></div><div>Incidentally, this raises this question down the road: in petsc-dev the inner and outer A00 can be configure separately using different prefixes, but how should the inner A00 DM be configured?</div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
<div class="gmail_quote"><div><div>
<div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div class="gmail_quote"><div> I can't reliably get the DM from the split's KSP,</div>
<div>so I have to stash it somewhere. Otherwise recursive splits don't work. If there is a better solution, I'm all for it.</div><div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
<div class="gmail_quote">
<div><br></div><div><div><font face="courier new, monospace">- /* set tabbing and options prefix of KSP inside the MatSchur */</font></div><div><font face="courier new, monospace">+ /* set tabbing, options prefix and DM of KSP inside the MatSchur */</font></div>
<div><font face="courier new, monospace"> ierr = MatSchurComplementGetKSP(jac->schur,&ksp);CHKERRQ(ierr);</font></div><div><font face="courier new, monospace"> ierr = PetscObjectIncrementTabLevel((PetscObject)ksp,(PetscObject)pc,2);CHKERRQ(ierr);</font></div>
<div><font face="courier new, monospace">+ {</font></div><div><font face="courier new, monospace">+ PC pcinner;</font></div><div><font face="courier new, monospace">+ ierr = KSPGetPC(ksp, &pcinner); CHKERRQ(ierr);</font></div>
<div><font face="courier new, monospace">+ ierr = PetscObjectIncrementTabLevel((PetscObject)pcinner,(PetscObject)pc,2);CHKERRQ(ierr);</font></div><div><font face="courier new, monospace">+ }</font></div>
</div><div><br></div><div>There is KSPIncrementTabLevel() for a reason. </div></div></blockquote></div><div>In petsc-dev, but not in petsc-3.3. </div><div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
<div class="gmail_quote"><div><br></div><div><div><font face="courier new, monospace">- ierr = PetscSNPrintf(schurprefix,sizeof schurprefix,"%sfieldsplit_%s_",((PetscObject)pc)->prefix?((PetscObject)pc)->prefix:"",ilink->splitname);CHKERRQ(ierr);</font></div>
<div><font face="courier new, monospace">- ierr = KSPSetOptionsPrefix(jac->kspschur,schurprefix);CHKERRQ(ierr);</font></div><div><font face="courier new, monospace"> /* really want setfromoptions called in PCSetFromOptions_FieldSplit(), but it is not ready yet */</font></div>
<div><font face="courier new, monospace"> /* need to call this every time, since the jac->kspschur is freshly created, otherwise its options never get set */</font></div><div><font face="courier new, monospace">- ierr = KSPSetFromOptions(jac->kspschur);CHKERRQ(ierr);</font></div>
<div><font face="courier new, monospace">+ ierr = PetscSNPrintf(schurprefix,sizeof schurprefix,"%sfieldsplit_%s_",((PetscObject)pc)->prefix?((PetscObject)pc)->prefix:"",ilink->splitname);CHKERRQ(ierr);</font></div>
<div><font face="courier new, monospace">+ ierr = KSPSetOptionsPrefix(jac->kspschur,schurprefix);CHKERRQ(ierr);</font></div><div><font face="courier new, monospace">+ ierr = KSPSetFromOptions(jac->kspschur); CHKERRQ(ierr);</font></div>
<div><font face="courier new, monospace">+ ierr = KSPSetUp(jac->kspschur); CHKERRQ(ierr);</font></div></div><div><br></div><div>Why this KSPSetUp()? </div></div></blockquote></div><div>I guess this KSPSetUp() can go. </div>
<div><div>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div class="gmail_quote"><div>You have changed the formatting here, but not explained why. (Also, can't you follow PETSc coding guidelines regarding placement of CHKERRQ?)</div>
<div><br></div><div><div><font face="courier new, monospace">- } else {</font></div><div><font face="courier new, monospace">- /* set up the individual PCs */</font></div><div><font face="courier new, monospace">- i = 0;</font></div>
<div><font face="courier new, monospace">- ilink = jac->head;</font></div><div><font face="courier new, monospace">- while (ilink) {</font></div><div><font face="courier new, monospace">- ierr = KSPSetOperators(ilink->ksp,jac->mat[i],jac->pmat[i],flag);CHKERRQ(ierr);</font></div>
<div><font face="courier new, monospace">- /* really want setfromoptions called in PCSetFromOptions_FieldSplit(), but it is not ready yet */</font></div><div><font face="courier new, monospace">- if (!jac->suboptionsset) {ierr = KSPSetFromOptions(ilink->ksp);CHKERRQ(ierr);}</font></div>
<div><font face="courier new, monospace">- ierr = KSPSetUp(ilink->ksp);CHKERRQ(ierr);</font></div><div><font face="courier new, monospace">- i++;</font></div><div><font face="courier new, monospace">- ilink = ilink->next;</font></div>
<div><font face="courier new, monospace">- }</font></div><div><font face="courier new, monospace">+ } </font></div><div><font face="courier new, monospace">+ /* set up the individual PCs */</font></div><div><font face="courier new, monospace">+ i = 0;</font></div>
<div><font face="courier new, monospace">+ ilink = jac->head;</font></div><div><font face="courier new, monospace">+ while (ilink) {</font></div><div><font face="courier new, monospace">+ ierr = KSPSetOperators(ilink->ksp,jac->mat[i],jac->pmat[i],flag); CHKERRQ(ierr);</font></div>
<div><font face="courier new, monospace">+ /* really want setfromoptions called in PCSetFromOptions_FieldSplit(), but it is not ready yet */</font></div><div><font face="courier new, monospace">+ if (!jac->suboptionsset) {ierr = KSPSetFromOptions(ilink->ksp); CHKERRQ(ierr);}</font></div>
<div><font face="courier new, monospace">+ ierr = KSPSetUp(ilink->ksp);CHKERRQ(ierr);</font></div><div><font face="courier new, monospace">+ i++;</font></div><div><font face="courier new, monospace">+ ilink = ilink->next;</font></div>
</div><div><br></div><div>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?</div></div></blockquote></div></div><div>In the previous version NONE of these were done for Schur. I suppose we can handle these cases more carefully.</div>
</div></blockquote><div><br></div></div></div><div>You haven't explained why that was wrong. You are certainly setting up more than necessary here.</div></div></blockquote><div>I backed out all of KSPSetUp() calls -- they occur inside KSPSolve() anyway, and only when necessary.</div>
<div>KSPSetFromOptions() also occurs more sparingly now:</div><div><a href="http://petsc.cs.iit.edu/petsc/releases/petsc-3.3/rev/3c0d43fb5911">http://petsc.cs.iit.edu/petsc/releases/petsc-3.3/rev/3c0d43fb5911</a></div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
<div class="gmail_quote"><div>
<div><br></div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
<div class="gmail_quote"><div>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div class="gmail_quote"><div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
<div><a href="http://petsc.cs.iit.edu/petsc/releases/petsc-3.3/rev/5a9ebf885615" target="_blank">http://petsc.cs.iit.edu/petsc/releases/petsc-3.3/rev/5a9ebf885615</a>)</div></blockquote><div><br></div><div><font face="courier new, monospace"> /* HACK: Check for the constant null space */</font></div>
<div><font face="courier new, monospace">- ierr = MatGetNullSpace(pc->pmat, &sp);CHKERRQ(ierr);</font></div><div><div><font face="courier new, monospace">+ /* </font></div><div><font face="courier new, monospace">+ This replaces a previous test that checked whether the split's constants are injected into the ambient problem's nullspace.</font></div>
<div><font face="courier new, monospace">+ That can create false negatives (declaring constants not in the nullspace, when they are) as in Stokes or A = [1, 1; 1, 0].</font></div><div><font face="courier new, monospace">+ */</font></div>
<div><font face="courier new, monospace"> if (sp) {</font></div></div><div><br></div><div>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.</div>
<div><br></div><div><div><font face="courier new, monospace"> ierr = MatGetNullSpace(jac->pmat[1], &sp);CHKERRQ(ierr);</font></div><div><font face="courier new, monospace">- if (sp) {ierr = MatSetNullSpace(jac->schur, sp);CHKERRQ(ierr);}</font></div>
<div><font face="courier new, monospace">- /* set tabbing, options prefix and DM of KSP inside the MatSchur */</font></div><div><font face="courier new, monospace">+ /* </font></div><div><font face="courier new, monospace">+ Do we really want to attach the A11-block's nullspace to S? Note that this may generate</font></div>
<div><font face="courier new, monospace">+ "false positives", when the A11's nullspace isn't S's: Stokes or A = [1, 1; 1, 0].</font></div><div><font face="courier new, monospace">+ Using a false nullspace may prevent KSP(S) from converging, since it might force an inconsistent rhs. </font></div>
<div><font face="courier new, monospace">+ */</font></div><div><font face="courier new, monospace">+ /* if (sp) {ierr = MatSetNullSpace(jac->schur, sp);CHKERRQ(ierr);} */</font></div><div><font face="courier new, monospace">+ /* set tabbing, options prefix and DM of KSP inside the MatSchurComplement */</font></div>
</div><div><br></div><div>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.</div>
<div><br></div><div><br></div><div>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". </div>
</div></blockquote></div><div>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?</div>
</div></blockquote><div><br></div></div><div>"force an inconsistent rhs"</div></div></blockquote><div>The KSP solve will fail only if the search space X is contracted to the point where the rhs b is not in the range AX, hence,</div>
<div>inconsistent. I backed out this hunk, however.</div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div class="gmail_quote"><div><br></div><div>S may well be nonsingular.</div>
<div><div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
<div class="gmail_quote"><div>
<div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div class="gmail_quote"><div>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.</div>
<div><br></div><div> 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.</div></div></blockquote></div><div>
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.</div>
</div></blockquote><div><br></div></div><div>1. I hate those hacky dm reference-non-references.</div><div><br></div><div>2. I don't want unused KSPs to be set up. (We should eventually fix the model so they don't exist...)<br>
<br>3. Two of the three IncrementTabLevels that you introduced are not necessary because the KSP's tab level was incremented before getting the PC. </div></div></blockquote><div>Where? </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
<div class="gmail_quote"><div>This unnecessary code is ugly and confusing.</div></div></blockquote><div>The current code is equivalent to KSPIncrementTabLevel(), introduced in petsc-dev for this very reason.</div><div>My recursive fieldsplit test cases produce what to me is very readable -ksp_monitor output with these tabbing options. </div>
<div>It wasn't the case before.</div><div><br></div><div>Dmitry.</div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div class="gmail_quote"><div><br></div><div>
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.</div></div>
</blockquote></div><br>