[petsc-dev] DMComplex and PCFieldSplit

Barry Smith bsmith at mcs.anl.gov
Fri Feb 24 16:16:15 CST 2012 

> I have DMComplexGetDefaultSection(), which is the first thing I call when constructing the FieldSplit info.

  Found it. So you want PCFieldSplit to troll through those inner loops of 

       /* Define fields */
        PetscSection section, sectionGlobal;
        PetscInt    *fieldSizes, **fieldIndices;
        PetscInt     nF, f, pStart, pEnd, p;

        ierr = DMComplexGetDefaultSection(pc->dm, &section);CHKERRQ(ierr);
        ierr = DMComplexGetDefaultGlobalSection(pc->dm, &sectionGlobal);CHKERRQ(ierr);
        ierr = PetscSectionGetNumFields(section, &nF);CHKERRQ(ierr);
        ierr = PetscMalloc2(nF,PetscInt,&fieldSizes,nF,PetscInt *,&fieldIndices);CHKERRQ(ierr);
        ierr = PetscSectionGetChart(sectionGlobal, &pStart, &pEnd);CHKERRQ(ierr);
        for(f = 0; f < nF; ++f) {
          fieldSizes[f] = 0;
        }
        for(p = pStart; p < pEnd; ++p) {
          PetscInt gdof;

          ierr = PetscSectionGetDof(sectionGlobal, p, &gdof);CHKERRQ(ierr);
          if (gdof > 0) {
            for(f = 0; f < nF; ++f) {
              PetscInt fcdof;

              ierr = PetscSectionGetFieldConstraintDof(section, p, f, &fcdof);CHKERRQ(ierr);
              fieldSizes[f] += fcdof;
            }
          }
        }
        for(f = 0; f < nF; ++f) {
          ierr = PetscMalloc(fieldSizes[f] * sizeof(PetscInt), &fieldIndices[f]);CHKERRQ(ierr);
          fieldSizes[f] = 0;
        }
        for(p = pStart; p < pEnd; ++p) {
          PetscInt gdof, goff;

          ierr = PetscSectionGetDof(sectionGlobal, p, &gdof);CHKERRQ(ierr);
          if (gdof > 0) {
            ierr = PetscSectionGetOffset(sectionGlobal, p, &goff);CHKERRQ(ierr);
            for(f = 0; f < nF; ++f) {
              PetscInt fcdof, fc;

              ierr = PetscSectionGetFieldConstraintDof(section, p, f, &fcdof);CHKERRQ(ierr);
              for(fc = 0; fc < fcdof; ++fc, ++fieldSizes[f]) {
                fieldIndices[f][fieldSizes[f]] = goff++;
              }
            }
          }
        }
        for(f = 0; f < nF; ++f) {
          IS          field;
          const char *fieldname;

          ierr = PetscSectionGetFieldName(section, f, &fieldname);CHKERRQ(ierr);
          ierr = ISCreateGeneral(((PetscObject) pc)->comm, fieldSizes[f], fieldIndices[f], PETSC_OWN_POINTER, &field);CHKERRQ(ierr);
          ierr = PetscPrintf(((PetscObject) pc)->comm, "Field %s\n", fieldname);CHKERRQ(ierr);
          ierr = ISView(field, PETSC_NULL);CHKERRQ(ierr);
          ierr = PCFieldSplitSetIS(pc, fieldname, field);CHKERRQ(ierr);
          ierr = ISDestroy(&field);CHKERRQ(ierr);
        }
        ierr = PetscFree2(fieldSizes,fieldIndices);CHKERRQ(ierr);
      }

every time instead of call DMGetFieldsIS(dm&is);   That's insane, just take the crud above and call it DMGetFieldsIS_Complex() and we'll have PCFieldSplit independent of any DM guts.

I don't understand what we are fighting about.


On Feb 24, 2012, at 3:49 PM, Matthew Knepley wrote:

> 
> It is the correct level of granularity, especially fr FEM which PETSc has almost none of right now
> because we make it unnaturally hard due to our ridiculous infatuation with Finite Difference. That kind
> of thinking colors everything.

   I never said that PetscSection doesn't belong anywhere. It just doesn't seem to belong anywhere in PC/KSP land, it is unnecessarily low for that. Sure it may be ok for the finite element assembly etc.

   Barry


>  
>  Regardless, say the concept of PetscSection as a way to keep the information about fields is fine. PetscSection is a better representation of an array of ISs.
> 
> 1)  Imbedding these constraint thingys in it seems completely wrong. Information about the constraint things should just be in another PetscSection not inside the PetscSection that tells you the field division.
> 
> I disagree since you need this information in every situation where you ask about sizes, etc. Again, there is only ONE
> problem in PETSc that does this, written by me. It is often useful to try something before trying to figure out what is the
> best way to do it.
>  
> 2) One could then have DMGetPetscSection(DM,section) that tells you the decomposition, but since the solvers need DMs for the pieces (in PCFieldSplit) it would be DMGetSections(DM,section,DM *subdms) and not so different from my first email with PetscSection replacing my array of ISs. So why not fix IS to be a bit more dynamic and build there is no real need for PetscSection except as an array of ISs.
> 
> I have DMComplexGetDefaultSection(), which is the first thing I call when constructing the FieldSplit info.
> 
>    Matt
>  
>  Jed is laughing his pants off at the reaction to his side comment in that email :-)
> 
> 
>   Barry
> 
> On Feb 24, 2012, at 3:32 PM, Matthew Knepley wrote:
> 
> > On Fri, Feb 24, 2012 at 3:26 PM, Barry Smith <bsmith at mcs.anl.gov> wrote:
> >
> > On Feb 24, 2012, at 3:13 PM, Matthew Knepley wrote:
> >
> > > On Fri, Feb 24, 2012 at 3:05 PM, Barry Smith <bsmith at mcs.anl.gov> wrote:
> > >
> > > On Feb 24, 2012, at 2:41 PM, Matthew Knepley wrote:
> > >
> > > > On Fri, Feb 24, 2012 at 2:39 PM, Barry Smith <bsmith at mcs.anl.gov> wrote:
> > > >
> > > > On Feb 24, 2012, at 2:23 PM, Jed Brown wrote:
> > > >
> > > > > We need a richer field interface. Some way for a client of DM to discover vector and tensor structure.
> > > >
> > > >   No idea what a tensor is, but DMGetSubDMS(dm,PetscInt *cnt,IS *ises,DM *dms);   and DMPushSubDMType(DM,"name of a type of splitting like u,v,p or velocity,p"); DMPopSubDMType(); and DMGetAvailableSubDMTypes(DM,char **);
> > > >
> > > >  Do you need more than that?
> > > >
> > > > Yuck Yuck Yuck. Why are we persisting in shoving all of the index stuff into the DM interface. I think this is so
> > > > unwieldy and complicated. PetscSection is small, powerful, and does everything suggested and more. I see
> > > > no argument in support of stuff like the above.
> > >
> > >  I have no idea what a PetscSection is and
> > >
> > > PetscSection - This is a mapping from DMMESH points to sets of values, which is
> > >  our presentation of a fibre bundle.
> > >
> > > doesn't help explain it.   PetscSection is small? It has more methods than DM, that is not small.
> > >
> > > Don't make me return the count. It has fewer, and definitely fewer concepts. Most of those are
> > > simple get/set.
> > >
> > > PetscSectionGetFieldComponents() seems to be intuitive, cool. By why does it return an array of integers? The universal way in PETSc would have it returning an IS. Wait a minute, it returns one int, the number of components? Then why is called GetFieldComponents and not GetNumFieldComponents()
> > >
> > > It could be changed to GetNumFieldComponents.
> > >
> > > And what is all this constraint stuff? PetscSectionGetConstraintDof() Huh? It looks inside the "section" and gets an array of ints from a section inside the section?
> > >
> > > Yes. It keeps track of constrained dofs. People, who through laziness or willful blindness, refuse to acknowledge that some
> > > problems make more sense when constrained dofs ar eliminated.
> > >
> > > Here's my guess: a PetscSection is way of "chopping up" a canonical vector into a bunch of fields. Fine, but then why does it have a bunch of other stuff in it (like constraints, what the hell are constraints)? And why isn't a PetscSection just an array of IS that define the fields?
> > >
> > > Constraints are just that. IS is designed to be immutable. Building on it internally is not a good idea for this reason.
> > >
> > > The reason to shove all the index stuff into the DM interface is because generally one doesn't just want to know the entries of a field but one also wants to create vectors for those field variables only, create matrices, know something about the mesh for that field so one can generate interpolations and coarsen for multigrid. One cannot do any of those things with PetscSection.
> > >
> > > Yes you can.
> >
> >   Hmm, I guess I haven't pulled recently. I don't see a PetscSectionGetInterpolation() a PetscSectionCoarsen() a PetscSectionGetGlobalVector() etc etc.
> >
> > >
> > >  Please explain in freshman linear algebra what a PetscSection is.
> > >
> > > Exactly what it says in the help. It is a mapping from mesh points (vertices, edges, faces, cells) to dofs. How can this be easier?
> >
> >   Say I want an IS that represents all the cell center dofs, how do I obtain that from a PetscSection?  I don't see the functionality for that but it seems the most basic kind of thing needed.
> >
> > DMComplexGetHeight(dm, 0, &cStart, &cEnd); /* Get the range of mesh points which are cells, height is like co-dimension */
> > numCellDof = 0;
> > for(c = cStart; c < cEnd; ++c) {
> >   PetscSectionGetDof(section, c, &dof);
> >   numCellDof += dof;
> > }
> > PetscMalloc(numCellDof * sizeof(PetscInt), &indices);
> > numCellDof = 0;
> > for(c = cStart; c < cEnd; ++c) {
> >   PetscSectionGetDof(section, c, &dof);
> >   PetscSectionGetOffset(section, c, &off);
> >   for(d = 0; d < dof; ++d, ++numCellDof) {
> >     indices[numCellDof] = off+d;
> >   }
> > }
> >
> >    Matt
> >
> >
> >   Barry
> >
> > >
> > >    Matt
> > >
> > >
> > >    Barry
> > >
> > >
> > >
> > >
> > >
> > >
> > > >
> > > >    Matt
> > > >
> > > >    Barry
> > > >
> > > > Note that current DMDA's default subdms served up would be a DM for each field consisting of one of the degrees of freedom per node and it could also serve up  DMs that represent any collections of those for example dof 0 and 1 per node, then dof 2 per node. We'll need to come up with some naming convention for these.
> > > >
> > > >
> > > > > Possibly also (eventually) a way to solve equations of state so that we can formulate reduced problems in non-conservative variables.
> > > > >
> > > > > On Feb 24, 2012 2:19 PM, "Matthew Knepley" <knepley at gmail.com> wrote:
> > > > > On Fri, Feb 24, 2012 at 1:58 PM, Jed Brown <jedbrown at mcs.anl.gov> wrote:
> > > > > http://petsc.cs.iit.edu/petsc/petsc-dev/rev/518ff70e8a0a
> > > > >
> > > > > Matt, I want to get rid of these conditionals, not add more. We should have a DM base interface for getting fields on which to split, that common interface should _replace_ the DMComposite specialization.
> > > > >
> > > > > I agree. However, I won't put anything in until I do it by hand once.
> > > > >
> > > > > Also, I was tempted to just promote DMGetGlobalISes(), but its not quite right. I am planning
> > > > > to put the IS method in PetscSection. I am hoping eventually DMDA uses PetscSection for
> > > > > layout.
> > > > >
> > > > >    Matt
> > > > >
> > > > > --
> > > > > 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
> > > >
> > > >
> > > >
> > > >
> > > > --
> > > > 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
> > >
> > >
> > >
> > >
> > > --
> > > 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
> >
> >
> >
> >
> > --
> > 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
> 
> 
> 
> 
> -- 
> 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

More information about the petsc-dev mailing list