On Tue, Sep 10, 2013 at 6:53 PM, Barry Smith <[email protected]> wrote:
> > On Sep 10, 2013, at 6:11 PM, Matthew Knepley <[email protected]> wrote: > > > On Tue, Sep 10, 2013 at 5:14 PM, Barry Smith <[email protected]> wrote: > > > > We had this discussion last November. Now that DMPlex has further > matured and is likely to play a big role in some major applications I'd > like to revisit the issue and see if PetscSection can be refactored to be > easier to understand and more useful and more "compatible" with how DMPlex > is developing. > > > > PetscSection seems to have two distinct roles: > > > > Role 1 (I'll call this object XX) Provides a way of indexing into > an array of "items" where different items can be of different sizes and > there can be spaces in the memory between items. > > > > --------- > > struct _n_PetscUniformSection atlasLayout; /* Layout for the atlas */ > > PetscInt *atlasDof; /* Describes layout of > storage, point --> # of values */ > > PetscInt *atlasOff; /* Describes layout of > storage, point --> offset into storage */ > > PetscInt maxDof; /* Maximum dof on any > point */ > > PetscBool setup; > > > > So if we have a regular C array pointer to type int, double etc we do > addr = array + p == &array[p] to get the address of the pth item (where p > may only be valid between pstart and pend.) With PetscSection instead one > would do XXGetOffset(xx,p,&offset); addr = (type*) (array + offset); while > XXGetDOF(xx,p,dof); gives you the "size" of the memory at location addr. > > > > For historical reasons the offsets are in terms of int* instead of > char*; (this should be fixed) > > > > I do not understand this. The offsets are integers, and I think they > should be. > > I was unclear. The offset value should be in bytes from the beginning > of the array, currently it is in ints from the beginning of the array. Ints > is just plan confusing, Shri has to have this strange size sizeof(his > struct)/sizeof(PetscInt) to indicate how big his chunk is. No, its has nothing to do with Int. Where did this come from? In fact, I use it all the time to mean Scalar offsets. It is interpreter by the creator of the Section. I told him to put in bytes for the sizes. Then you do not need any of that crap. > > > > > > This is a useful thing but in itself not a killer app. What makes > it a killer app are > > > > DMPlexDistributeField(DM dm, PetscSF pointSF, XX originalSection, Vec > originalVec, XX newSection, Vec newVec) > > DMPlexDistributeData(DM dm, PetscSF pointSF, XX originalSection, > MPI_Datatype datatype, void *originalData, XX newSection, void **newData) > > > > Note that these do not actually even use the DM (the code should be > fixed by moving it out of DMPlex into the "lower level" code.) Note also > these really create newSection internally (the code should be fixed by > making these take PetscSection *newSection). > > > > I can move them over to vsectionis.c. I do not agree that the signature > should change. I thought we were going to > > the VecLoad() paradigm where an empty object comes in. > > > > What this means is that if I have a set of "points" distributed > across MPI processes and for each point a "bunch of data" (in an "array" > offsetted by a XX) and I redistribute the points across MPI processes > including adding ghost points then I can redistribute all the "bunch of > data" (including that needed for ghost points) trivially and conveniently > access it using the newSection. > > > > Yep, this should be very useful, an extension of the very useful > VecScatter. > > > > > > Role 2 (I'll call this object CC) Provides a container that tells > you which parts of "bunches of data" in an array or Vec are associated with > different "physical simulation fields" and are possibly constrained (for > example by boundary conditions). Uses XX to indicate each field etc. > > > > -------- > > PetscSection bc; /* Describes constraints, > point --> # local dofs which are constrained */ > > PetscInt *bcIndices; /* Local indices for > constrained dofs */ > > > > > > PetscInt numFields; /* The number of fields > making up the degrees of freedom */ > > const char **fieldNames; /* The field names */ > > PetscInt *numFieldComponents; /* The number of > components in each field */ > > PetscSection *field; /* A section describing > the layout and constraints for each field */ > > > > > > Suggestion: > > > > We split PetscSection into XX and CC reorganize and clean them > up. > > > > Okay. The reason I did it this way was that any FEM code is going to > need CC, and it needs to contain XX, so > > I just made XX the container. We can make CC the container and have the > split that you want. It makes sense. > > > > Relationship of XX to IS. > > > > An IS defines a subset of entries in a Vec or possibly in an > array of type int. Generally individual entries in an IS do not have > particular meaning relative to other entries. That is one does not iterate > over entries in an IS, one just uses the IS to get a new Vec and then does > whatever on that Vec. > > > > An XX is something you often iterate over, for example > > > > for (p = pStart; p < pEnd; ++p) { > > PetscInt dof, off, s; > > ierr = PetscSectionGetDof(mesh->supportSection, p, > &dof);CHKERRQ(ierr); > > ierr = PetscSectionGetOffset(mesh->supportSection, p, > &off);CHKERRQ(ierr); > > do something > > > > that is you get a single entry of XX and do something with it (though > you can also get a bunch of entries and do something on the whole bunch). > > > > We do have these kind of offset structures in PETSc (MatAIJ), but they > are hardcoded in. > > > > Question: > > > > Is XX a powerful and useful thing? Do we have the correct API for > it? Can it be used for other purposes? > > > > I think its all yes. > > > > Can it be made more powerful: for example the "chunk" that it gives > you is always a contiqous set of bytes in the "array". Would be be better > off with a YY object that allowed them to be non-contiquous? Or should YY > be a higher level object implemented with XX? > > > > Let me take this statement apart. So for any point p, the Section gives > you back (offset, dof), so that you can use it to refer to > > > > array[offset]...array[offset+dof] > > > > which is the "chunk" that Barry is talking about. It is contiguous by > design. That is how we achieve compression in the representation. I claim > you > > can get the YY object with two XX objects. The first XX object says how > many "points" in the second XX correspond to each chunk. Each one of > > the subchunks in the second XX is contiguous but they create a > non-contiguous chunk for the first XX. In general, any non-contiguous thing > is > > represented by a set of contiguous things. This is not a good > representation is you mostly have scatter, scalar data. However, I think we > mostly have > > scattered blocks of data. > > Is this second thing, built on two XX worthing of its own abstraction? > For example when I way give me the coordinates of the nodes of a triangular > in a mesh? I am not sure I understand what you want here. I don't think this YY thing is all that useful since I can't see the generality. Matt > > > > Looking for a vigorous discussion and then code that I can better > understand. > > > > Me too, so "Jed is probably all wrong about this" :) > > > > Matt > > > > > > Barry > > > > > > > > On Nov 8, 2012, at 10:18 PM, Jed Brown <[email protected]> wrote: > > > > > This is inevitable, but may still not resolve anything. I'm sure Matt > will correct whatever below is incorrect. The main question is whether to > promote PetscSection to a PetscObject or to hide/replace it with something > else. > > > > > > So PetscSection is an unloved (by most) beastie that is sort of > intermediate between an IS and a DM. Its core functionality is to provide > indirect indexing with variable block sizes. It consists of two parts, the > first of which is less than a PetscLayout. > > > > > > struct _n_PetscUniformSection { > > > MPI_Comm comm; > > > PetscInt pStart, pEnd; /* The chart: all points are contained in > [pStart, pEnd) */ > > > PetscInt numDof; /* Describes layout of storage, point --> > (constant # of values, (p - pStart)*constant # of values) */ > > > }; > > > > > > numDof is always 1 and PetscUniformSection is never referenced from a > *.c file so this is really just a range of valid "point" values > [pStart...pEnd). > > > > > > struct _n_PetscSection { > > > struct _n_PetscUniformSection atlasLayout; /* Layout for the atlas > */ > > > PetscInt *atlasDof; /* Describes layout of > storage, point --> # of values */ > > > PetscInt *atlasOff; /* Describes layout of > storage, point --> offset into storage */ > > > > > > ^^ These two fields are the essential part. Each index in > [pStart,pEnd) is mapped to a block of size atlasDof[p-pStart] located at > atlasOff[p-pStart]. Matt composes these indirect maps for storing meshes > and other indexing tasks. This "unsorted offset mapping" is indeed used > frequently in PETSc, though outside of Matt's code, it's just raw arrays. > > > > > > User's primarily encounter PetscSection when accessing unstructured > meshes or any time the number of degrees of freedom or locations of each > "point" (topological entity in the mesh) is not uniform. DMs do a similar > thing (provide more logical access to a Vec), but DMDAVecGetArray() is of > course only natural for a structured grid. > > > > > > At the user interface level, provided the user doesn't need low-level > access to topology, I think there are three possibilities. > > > > > > 1. To compute a residual at a cell, user does > > > > > > VecGetArray(F,&f); // normal access to the Vec > > > DMComplexGetHeightStratum(dm,0,&cstart,&cend); // "height 0" is like > codimension 0 with respect to the maximum topological dimension. Matt and I > have have a philosophical debate over whether users prefer to think in > "strata" (that can skip dimensions, so "height 1" could be a face, an edge, > or a vertex depending on what is stored in the mesh) or in "topological > dimensions" (where codim=1 is a "face" for a 3D mesh and an edge for a 2D > mesh). > > > > > > DMGetDefaultSection(dm,§ion); > > > for (c=cstart; c<cend; c++) { > > > PetscInt off; > > > PetscSectionGetOffset(section,c,&off); > > > f[off+0] = residual of first equation at this cell; > > > f[off+1] = residual of second equation; > > > // brushing under the rug how we access neighbor values > > > } > > > > > > The above is basically the current model. Matt has > DMComplexGet/SetClosure(), but that interface carries a sleeper reference > to the Vec's array and he says he plans to remove it. > > > > > > 2. Add DMComplex-specific access routines so the user does not need to > see the PetscSection. Presumably this would be something like > > > DMComplexGetPointOffset(dm,c,&offset); // offset into owned part of > global vector? > > > DMComplexGetPointOffsetLocal(dm,c,&loffset); // offset into local > vector > > > > > > 3. Access cursors (another object). > > > DMComplexGetCursors(dm,Vec U,3,&uface,&ucellL,&ucellR); > > > // similar for F > > > DMComplexGetHeightStratum(dm,1,&fstart,&fend); > > > for (f=fstart; f<fend; f++) { > > > const PetscInt *c; > > > const PetscScalar *uL,*uR,*uF; > > > DMComplexGetSupport(dm,f,&c); // left and right cells > > > DMCursorGetValues(uface,f,&uF); // only used for staggered/mixed > discretization > > > DMCursorGetValues(ucellL,c[0],&uL); > > > DMCursorGetValues(ucellR,c[1],&uR); > > > // Solve Riemann problem with "face" values uF[] (if using a > staggered discretization), uL[], and uR[]. > > > DMCursorRestoreValues(uface,f,&uF); // ... > > > > > > The "advantage" of the cursors here is that DMComplex doesn't need to > hold the context manager itself (perhaps for multiple threads), the code > dealing with the number of concurrent accesses is outside the inner loop, > and the implementations of these accessors are trivial (offset). I think > it's still not compelling for the sketched discretization above, but > consider also FEM methods where we need to access all degrees of freedom on > the closure of an element. The following packs a buffer with the closure > (because it's not contiguous in memory). > > > > > > DMCursorGetClosure(cell,c[0],&uc); > > > > > > > > > 4. Alternative suggestions? > > > > > > > > > I'm so sick of the name DMComplex. Let's just rename to DMPlex. It's > shorter and doesn't have a double-meaning. Or pretty much any other name > but "Complex". > > > > > > For completeness, the last bit of _n_PetscSection: > > > > > > PetscInt maxDof; /* Maximum dof on any > point */ > > > PetscSection bc; /* Describes > constraints, point --> # local dofs which are constrained */ > > > PetscInt *bcIndices; /* Local indices for > constrained dofs */ > > > > > > These BC fields are for doing masked updates, such as > VecSetValuesSection() which ignores constrained dofs. > > > > > > PetscInt refcnt; /* Vecs obtained with > VecDuplicate() and from MatGetVecs() reuse map of input object */ > > > PetscBool setup; > > > > > > PetscInt numFields; /* The number of fields > making up the degrees of freedom */ > > > const char **fieldNames; /* The field names */ > > > PetscInt *numFieldComponents; /* The number of > components in each field */ > > > PetscSection *field; /* A section describing > the layout and constraints for each field */ > > > }; > > > > > > > > > > > > > -- > > 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
