[petsc-dev] PetscSection

[Jed Brown]

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

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().
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-dev/attachments/20121109/c7aa3631/attachment.html>