[petsc-dev] PetscSection

[Barry Smith]

On Nov 9, 2012, at 7:17 PM, Jed Brown <jedbrown at mcs.anl.gov> wrote:

> On Fri, Nov 9, 2012 at 5:12 PM, Barry Smith <bsmith at mcs.anl.gov> wrote:
> 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.
> 
> Matt and I have been debating a similar point independently. I think we concluded that we could make do with index sets, DMImpl-based vector access, and and sub-DMs (for splitting out a piece of a larger problem, e.g., to solve an auxiliary problem along a fault).
>  
> 
>       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.
> 
> The key feature is accessing a range by block number. In the case of Section, the linear scalar indices (from concatenating the scalar entries for each block) is not used directly.
>  
> 
>     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?
> 
> Matt does not explicitly store the vector indices with support on the element. Instead he has the "closure" operation which, given an cell number, gathers all the lower-dimensional topological "points" (faces, edges, and vertices) using a breadth-first traversal. (This is something that might make sense to cache explicitly for each element. It uses more memory, but it makes the access pattern simpler so might pay off. Most FEM codes do that.)

   Just an implementation issue. To me the correct abstract model is "indices of an element" etc and whether this is done via a "closure" or listing them etc is just an implementation issue. 

   I realize this is different than Matt's abstract model. I am not saying my model is better than Matt's (from an abstract mathematicians Matt's is probably better) but I believe my model is much more sellable to the masses (remember very few abstract mathematicians use PETSc :-)).

   So we are set? We design and implement a better, more general IS concept?


   Barry

> 
> Each of those "points" in the closure has zero or more dofs in the Vec, as indicated by the "default" Section.
> 
> for each element e:
>     PetscScalar *cu,*cf;
>     // need workspace for cf
>     DMPlexVecGetClosure(dm,U,e,&csize,&cu); // I'm just going to pretend that he has abandoned the colliding name
>     ComputeElementResidual(csize,cu,cf); // geometry and function space conveniently tucked under the rug
>     DMPlexVecRestoreClosure(dm,U,e,&csize,&cu);
>     DMPlexVecSetClosure(dm,F,e,csize,cf,ADD_VALUES);
> 
> 
> 
> Note that there is no PetscSection above because its use is tucked inside the convenience routines which internally rely on DMComplexGetTransitiveClosure() and then figure out what part of the Vec to access using PetscSectionGetOffset() and PetscSectionGetDof().