[petsc-dev] PetscSection

[Matthew Knepley]

On Thu, Nov 8, 2012 at 11: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).

This allows blocks, but as Jed said I have not done anything with blocks yet.

> 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,&section);
> 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.

This has been redone to remove Jed's objection (see new manual section). It now
uses DM temp arrays,

for(c = cStart; c < cEnd; ++c) {
   PetscInt numVals;
   PetscScalar *vals;

   DMComplexVecGetClosure(dm, section, vec, c, &numVals, &vals);
   /* Compute residual */
   DMComplexVecRestoreClosure(dm, section, vec, c, &numVals, &vals);
   DMComplexVecSetClosure(dm, section, resvec, c, vals, ADD_VALUES);
}

> 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

This is cumbersome because you must make another DM for every PetscSection
you want to use. However, I am not completely against this now because
DMComplexClone()
is easy and accomplishes this.

> 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);

I hate cursors. I had the same experience with them no matter what
they are called
(iterators, etc.) You need so much information, that you end up with
the whole object
in this "external" thing. I think they never pay off.

> 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.

They are also for creating reduced systems.

>   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 */
> };

These are for multi-field splitting/

     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