[petsc-dev] DMComplex and PCFieldSplit
Matthew Knepley
knepley at gmail.com
Fri Feb 24 16:23:48 CST 2012
On Fri, Feb 24, 2012 at 4:16 PM, Barry Smith <bsmith at mcs.anl.gov> wrote:
>
> > 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, §ion);CHKERRQ(ierr);
> ierr = DMComplexGetDefaultGlobalSection(pc->dm,
> §ionGlobal);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.
Yes, I agree with that. I will move it.
Matt
>
> 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
>
>
--
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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-dev/attachments/20120224/3d97b051/attachment.html>
More information about the petsc-dev
mailing list