[petsc-dev] PetscSection
Shri
abhyshr at mcs.anl.gov
Tue Sep 10 20:36:58 CDT 2013
On Sep 10, 2013, at 6:57 PM, Matthew Knepley wrote:
> On Tue, Sep 10, 2013 at 6:53 PM, Barry Smith <bsmith at mcs.anl.gov> wrote:
>
> On Sep 10, 2013, at 6:11 PM, Matthew Knepley <knepley at gmail.com> wrote:
>
> > On Tue, Sep 10, 2013 at 5:14 PM, Barry Smith <bsmith at mcs.anl.gov> 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.
After speaking with Barry today, I realized that I was abusing PetscSection and doing a roundabout indexing to access the desired data. I was using four sections to distribute data for four different structures. I understand now that I need to
i) Create only one section to manage the parameters and use PetscSectionAddDof() to add the dofs (sizes of structs) for the different structs required for residual evaluation.
ii) Allocate a chunk of memory equal to the size returned from PetscSectionGetStorageSize(). Note that the user can do (i) and (ii) directly while reading the circuit parameters file. I am doing this after the file has been read and the different structs have been created, hence the need for allocating the chunk of memory.
iii) Loop over the edges, vertices and copy over the parameters using PetscSectionGetOffset().
iv) Use DMPlexDistributeData() with this section and the chunk of memory.
I tried to use MPI_BYTE as the MPI_DataType in DMPlexDistributeData() but got the following error. Jed explained that it had something to do with how the data is packed and suggested to divide the sizes of the structs in (ii) by PetscInt and use MPI_INT for DMPlexDistributeData() as a workaround.
[0]PETSC ERROR: --------------------- Error Message ------------------------------------
[0]PETSC ERROR: No support for this operation for this object type!
[0]PETSC ERROR: No support for type size not divisible by 4!
[0]PETSC ERROR: ------------------------------------------------------------------------
[0]PETSC ERROR: Petsc Development GIT revision: d45619dec29bfb59cf96225a84e0a74106da50ca GIT Date: 2013-08-17 23:45:05 -0500
[1]PETSC ERROR: --------------------- Error Message ------------------------------------
[1]PETSC ERROR: No support for this operation for this object type!
[1]PETSC ERROR: No support for type size not divisible by 4!
[1]PETSC ERROR: ------------------------------------------------------------------------
[1]PETSC ERROR: Petsc Development GIT revision: d45619dec29bfb59cf96225a84e0a74106da50ca GIT Date: 2013-08-17 23:45:05 -0500
[0]PETSC ERROR: See docs/changes/index.html for recent updates.
[0]PETSC ERROR: See docs/faq.html for hints about trouble shooting.
[0]PETSC ERROR: See docs/index.html for manual pages.
[0]PETSC ERROR: ------------------------------------------------------------------------
[0]PETSC ERROR: ./PF2 on a debug-mode named Shrirangs-MacBook-Pro.local by Shri Tue Sep 10 20:22:44 2013
[0]PETSC ERROR: Libraries linked from /Users/Shri/Documents/petsc/debug-mode/lib
[0]PETSC ERROR: Configure run at Sat Aug 31 09:27:36 2013
[0]PETSC ERROR: Configure options --download-chaco --download-metis --download-mpich --download-mumps --download-parmetis --download-scalapack --download-superlu_dist COPTFLAGS=-g CLFAGS=-g FLAGS=-g FOPTFLAGS=-g PETSC_ARCH=debug-mode
[1]PETSC ERROR: See docs/changes/index.html for recent updates.
[1]PETSC ERROR: See docs/faq.html for hints about trouble shooting.
[1]PETSC ERROR: See docs/index.html for manual pages.
[1]PETSC ERROR: ------------------------------------------------------------------------
[0]PETSC ERROR: ------------------------------------------------------------------------
[0]PETSC ERROR: PetscSFBasicPackTypeSetup() line 386 in /Users/Shri/Documents/petsc/src/vec/is/sf/impls/basic/sfbasic.c
[0]PETSC ERROR: PetscSFBasicGetPack() line 509 in /Users/Shri/Documents/petsc/src/vec/is/sf/impls/basic/sfbasic.c
[0]PETSC ERROR: [1]PETSC ERROR: ./PF2 on a debug-mode named Shrirangs-MacBook-Pro.local by Shri Tue Sep 10 20:22:44 2013
[1]PETSC ERROR: Libraries linked from /Users/Shri/Documents/petsc/debug-mode/lib
[1]PETSC ERROR: Configure run at Sat Aug 31 09:27:36 2013
[1]PETSC ERROR: Configure options --download-chaco --download-metis --download-mpich --download-mumps --download-parmetis --download-scalapack --download-superlu_dist COPTFLAGS=-g CLFAGS=-g FLAGS=-g FOPTFLAGS=-g PETSC_ARCH=debug-mode
[1]PETSC ERROR: ------------------------------------------------------------------------
[1]PETSC ERROR: PetscSFBcastBegin_Basic() line 645 in /Users/Shri/Documents/petsc/src/vec/is/sf/impls/basic/sfbasic.c
[0]PETSC ERROR: PetscSFBcastBegin() line 918 in /Users/Shri/Documents/petsc/src/vec/is/sf/interface/sf.c
[0]PETSC ERROR: DMPlexDistributeData() line 2796 in /Users/Shri/Documents/petsc/src/dm/impls/plex/plex.c
[0]PETSC ERROR: PetscSFBasicPackTypeSetup() line 386 in /Users/Shri/Documents/petsc/src/vec/is/sf/impls/basic/sfbasic.c
[1]PETSC ERROR: PetscSFBasicGetPack() line 509 in /Users/Shri/Documents/petsc/src/vec/is/sf/impls/basic/sfbasic.c
[1]PETSC ERROR: PetscSFBcastBegin_Basic() line 645 in /Users/Shri/Documents/petsc/src/vec/is/sf/impls/basic/sfbasic.c
main() line 582 in pf2.c
[1]PETSC ERROR: PetscSFBcastBegin() line 918 in /Users/Shri/Documents/petsc/src/vec/is/sf/interface/sf.c
[1]PETSC ERROR: DMPlexDistributeData() line 2796 in /Users/Shri/Documents/petsc/src/dm/impls/plex/plex.c
[1]PETSC ERROR: main() line 582 in pf2.c
application called MPI_Abort(MPI_COMM_WORLD, 56) - process 0
[cli_0]: aborting job:
application called MPI_Abort(MPI_COMM_WORLD, 56) - process 0
application called MPI_Abort(MPI_COMM_WORLD, 56) - process 1
[cli_1]: aborting job:
application called MPI_Abort(MPI_COMM_WORLD, 56) - process 1
===================================================================================
= BAD TERMINATION OF ONE OF YOUR APPLICATION PROCESSES
= EXIT CODE: 56
= CLEANING UP REMAINING PROCESSES
= YOU CAN IGNORE THE BELOW CLEANUP MESSAGES
===================================================================================
Btw, DMPlexDistributeData() is missing an extern declaration in petscdmplex.h.
>
>
> >
> >
> > 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 <jedbrown at mcs.anl.gov> 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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-dev/attachments/20130910/7f14cad2/attachment.html>
More information about the petsc-dev
mailing list