[petsc-dev] PetscSection

[Barry Smith]

On Nov 8, 2012, at 10:51 PM, Matthew Knepley <knepley at gmail.com> wrote:

> On Thu, Nov 8, 2012 at 11:46 PM, Barry Smith <bsmith at mcs.anl.gov> wrote:
>> 
>>  Why isn't it
>>> 
>>> struct _n_PetscSection {
>>>  PetscInt  firstIndex,lastIndex;   /* and why not just have this start at zero and go to n? */
>>>  PetscInt *dofs;
>>>  PetscInt *offsets;
>> };
>> 
>>   all that other clutter in there makes it mighty confusing as to the purpose of the beast.
> 
> I would not fight that.
> 
>>   Also why is the name PetscSection?    If I have an array of numbers this just seems to organize them  into a series of blocks (and allow traversals through the blocks), for example,
>> 
>>    0 1 2 3 4 5 6 7 8 9 10 11 12 13 14
>>    X X    Y Y Y      Z Z   Z
>> 
>>   Now a block IS allows traversals through a series of blocks (but all blocks must be the same size).  It is important to note that  though you can think of a block IS as simply an implementation of a IS (that is a collection of integers) that is more efficient than listing all the integers in an ISGeneral, it actually has additional information (each block as its own fundamental unit).
>> 
>>  With this point of view, how is a PetscSection more than a block IS with variable block size? And why not implement it that way?
> 
> Because you can do more general things. The offsets do not have to be
> prefix sums of the sizes. This allows
> me to make views into larger vectors, make "global" sections with
> negative offsets for off-process values,
> etc. IS is not flexible enough even for the things I can dream up, let
> alone users.

   Matt,

      You are "down defining" IS so you can say it doesn't support doing these things. I am proposing "up defining" IS to be able to do these things. 

      I consider IS to be _THE_  PETSc tool for "indexing  stuff" (and want it to be the only tool, having two tools is redundant and confusing). When IS doesn't have a capability for some "kind of indexing" we wish to do then I think the right thing to do is generalize IS as needed, not introduce a different class. 

      Consider the normal scenario for finite elements:  (simple case)
      
          Loop over elements:
               Get the vector entries Ue from U for the element
               Compute with Ue producing Fe
               Put the Fe into the output vector F.

        this is an example of indexing.  Plain old ISGeneral doesn't directly provide a way to do this,  Having a concept of a collection of groups of indices (of varying sizes for different element types) makes this kind of finite element scenario straight forward.  Your PetscSection seems like yet a different type of useful indexing.

    What I propose is we "start at the beginning" using all the ideas we've gained over the years and design a new IS (which we can call PetscIS (or something else)) using good abstractions of "indexing stuff". Once we have this good design
IS will quickly be removed and replaced with the new IS construct. I am open to any idea on what the model (set of abstractions, whatever ..) we use for the new IS. 

    Though you use fancy abstract words with PetscSection like Section and Atlas you actually have only one concrete specific implementation that is merely block size and ranges. Is that really all we want/need for a very general and useful concept of efficient indexing?  Perhaps it is, but let's explore it. For example in pseudo code like above show me how you would do the normal scenario for finite elements using your PetscSection?

  Barry




> 
>   Matt
> 
>>   My hatred of more than a tiny number of abstract classes in PETSc has served me well over the years, and though I cheat sometimes to maintain this small number I believe that having excessive classes merely shows you don't understand your domain well enough.
>> 
>>   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,&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.
>>> 
>>> 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