[petsc-dev] DM field interfaces

[Matthew Knepley]

On Sat, Feb 25, 2012 at 3:27 PM, Barry Smith <bsmith at mcs.anl.gov> wrote:

>
>   Let's back up a step before everyone gets totally confused (and worry
> less about exact syntax to use instead focus on concepts)
>
> Proposal 1: Jed
>
>    // Conservative gas dynamics with fields ordered as [rho, rho*u, rho*v,
> rho*w, E]
>   DMGetFieldSplitting(dm,"conservative-vector",&fscons); // creates if it
> doesn't exist
>   DMFieldSplittingAddField(dm,fscons,"rho",{0});
>   DMFieldSplittingAddField(dm,fscons,"momentum",{1,2,3});
>
>     This seems to be premised on having an underlying "canonical" set of
> fields that one can define unions of to make new fields. For example with
> DMDA the underlying
>    fields are simply the 0th, 1st, 2nd etc DOF at each node. For a simple
> staggered grid 0 canonical may be the pressure centers, 1 the x velocities,
> 2 the y velocities
>
> Proposal 2: Barry
>
>     Use IS's to define the fields, not unions of canonical fields. Reason:
> more general than proposal 1
>
> Proposal 3: Jed
>
>    Observation: Using IS does not give all desired generality since the
> new field may be a linear combination of DOFs. Or, yikes, a nonlinear
> transformation of DOF.
>    Thus we may need something like DMTransformVec() to do these
> transformations.
>
>
> Comment:  Barry
>
>     I very much like the universality in PETSc of using IS __ALWAYS__ to
> grab a set of indexed things and not using another indexing technique like
> 0 to indicate the "first field"
>     and 1 to indicate the "second field".
>
>       1) Matt seems to like this too since he very recently lobbied for
> simplifying PCFIELDSPLIT to ONLY use the IS for indicating parts and to
> remove
>            the   (somewhat redundant) PetscInt          nfields;  PetscInt
>          *fields; which we will do.
>

It is true, I agree with Barry here.


>        2) We currently have VecStrideXXXXX(vec,start,....) that index
> using the 0th field, 1st field etc Does not match the paradigm
>
>        3) VecGetSubVec(Vec,IS,Vec*) uses IS to pull out pieces so matches
> the paradigm. Way to go Jed!
>
>        4) There may be other places that do not match the paradigm
>
>  Because Jed's proposal 1 violates the use of IS for indexing paradigm and
> is less general than Proposal 2 I'd like to drop that Proposal 1; but it
> does bring up an issue. I think we need another IS implementation that is
> more "general" than stride for handling sets of DOFs of freedom associated
> with those canonical sets of fields. So it is easy to access canonical
> fields that are not simply stride. For example
> ISCreateCanonical(MPI_Comm,PetscInt ncanonical,PetscInt *fields), we can
> even have "predefined" ones available for use such as
> IS_CANONICAL_(MPI_Comm,3,{1,2,3}) Thus we can preserve some of the
> simplicity of proposal 1.  Jed's (again not worrying about syntax)
>

I would prefer not proliferating IS types here.


> DMFieldSplittingAddField(dm,fscons,"momentum",{1,2,3});  becomes
>

Here is how I understand this call. It says "I would like a field for FS
(meaning a set of indices) from
the thing this DM calls {1,2,3}". I do not see any advantage of shoving
{1,2,3} into an IS wrapper. What
is wrong with this sequence

  DMFieldSplittingAddField(dm, "momentum", [1,2,3])
  DMFieldSplittingAddField(dm, "mass", [0])
  DMFieldSplittingAddField(dm, "energy", [5])
  DMCreateFieldIS(dm, &numFields, &fieldNames, &fieldISes);

1) The ISes here are exactly what we expect, sets of integers, so FS never
has to deal with anything else.

2) We can move all the field splitting stuff to specific DM classes so it
is easy and expressive

3) We just need each DM to respond to CreateFieldIS() properly

4) I acknowledge that it does not solver the more general problem, but I
don't think its necessary right now

   Matt



> DMFieldSplittingAddField(dm,fscons,"momentum",IS_CANONICAL_(comm,3,{1,2,3})
>
> Note that for a Vec with a given block size the canonical fields are
> simply the stride fields so IS_CANONICAL_(comm,1,{start}) is pretty much
> ISStride().
>
> We could push this even further and assume that the underlying object
> being indexed in has "named" fields, we could still use those named fields
> to index via the IS if we had
> and ISCreateNamed(comm,"Name") and helper IS_NAMED_(comm,"Name") so one
> could do VecGetSubVec(vec,IS_NAMED_(comm,"T"),&subvec); etc.
>
> Ok, so this solves the original issue of indicating sets of indices for
> DM, useful for standard use of block methods for matrices :-). It allows
> general IS, canonical fields and named fields all under the same interface.
>
> ----------------
>
> Now to the hard stuff Jed dumped on us.
>
> Proposal 3: Jed
>
>    Observation: Using IS does not give all desired generality since the
> new field may be a linear combination of DOFs. Or, yikes, a nonlinear
> transformation of DOF.
>    Thus we may need something like DMTransformVec() to do these
> transformations.
>
>
> So this moves us beyond the traditional block linear algebra methods but
> if we could solve it we'd have some very nice new powerful capabilities on
> our hands.
>
> Question 1:  Do we need to solve this for the current generation of
> PCFieldsplit?
> Question 2:  Do we want to solve it in PCFieldSplit or should that come
> later in a different PC, PCBasisChangeThingy or something?
> Question 3: Can we go ahead with deciding on exact syntax etc for Proposal
> 2 and implementing before we resolve the Proposal 3 issue?
>
>
> When I see DMTransformVec() I get worried, we already have
> DMGetInterpolation/Restriction() that maps between two DMs, cannot we not
> somehow extend that instead of putting in more methods? (I know they return
> Mat which is a linear operator and if I had listened to Matt many years ago
> everything in PETSc would be a nonlinear operator so everything would be
> fine).  If at all possibly I'd to reuse/fix up those old interfaces to
> handle this need.
>
>  Barry
>
>
>
>
>
>
> On Feb 25, 2012, at 1:33 PM, Jed Brown wrote:
>
> > On Sat, Feb 25, 2012 at 13:10, Dmitry Karpeev <karpeev at mcs.anl.gov>
> wrote:
> > I wrote this before the thread moved to petsc-dev, so there is some
> overlap with the below.
> > A given problem with the same "topology" -- a mesh and a total index
> layout -- might define several useful splits -- collections of fields. Each
> field might have further useful splits. I have useful examples from my
> interaction with Moose/libMesh, although these are rather general.  A field
> is now represented by an IS.  We could try to put netsted split definitions
> into ISs (e.g., ISNest), but I think it's much more natural to keep all
> nestedness inside DMs -- we already have examples of this in the
> coarsening, refinement for MG and a pretty well understood interaction of
> this sort of DM hierarchy. with a hierarchy of PCs.
> >
> > Thus, I think it would be useful for a DM to return a split as a list of
> DMs, which can be further split, if necessary:
> >        DMGetFieldDMs(DM dm, PetscInt *fieldcount, DM **fielddms);
> > The field names can simply be the corresponding DM names, etc.  The
> fielddms can be lightweight, sharing the underlying
> > topology.  That's up to the implementation to arrange, though.  This
> addresses nested splits.
> >
> > There are two concepts, defining what the space for a given split _is_
> and how to get the part of the current state is needed to get a vector
> living in that space. The IS is a rudimentary way to define what part is
> needed. (You can use VecGetSubVector(X,isfield,&Xfield) or other methods to
> reference that part of the state.) Unfortunately, it's not sufficient even
> for a linear change of basis. An IS combined with a PF is sufficient for
> non-mixed spaces. Mixed/staggered spaces need even more because the
> transformation is defined through quadrature and projection. Maybe the
> transformation should be entirely handled by the DM instead of giving the
> user objects that apply the transformation.
> >
> > But is the interface going to be DMTransformVec(dm,X,Y,dmt,Yt) that
> transforms the perturbation Y to the state X into transformed variables
> according to dmt? Or should we be getting back an object that does this
> sort of transformation, so that we can amortize some setup coming from
> matching dmt relative to dm? (Note that we also need an inverse
> transformation.)
> >
> >
> > We could use the same mehanism to implement alternative splits on the
> same topology  -- these are different DMs sharing the topology internally.
>  Instead of pushing and popping splits, they can be retrieved by name as
> >        DMGetSplitDM(DM dm, const char *splitname, DM *splitdm);
> > Then calling DMGetFieldDMs on the splitdm retrieves the corresponding
> fields.
> > Again, splitdm and dm can share common data structures, as arranged  by
> the particular DM implementation.
> > The advantage of this approach over push/pop is that we can
> simultaneously use different splits of the same topology.
> >
> > It is rather easy to apply the same scheme and have DMs serve up
> (overlapping) domain decompositions useful for (G)ASM.
> > I'm now implementing both of these schemes (fieldsplit and subdomains)
> on top of libMesh, but the code resides in the libMesh source tree to avoid
> a circular dependence: as designed right now libMesh depends on PETSc.  The
> API has to live in DM, though.
> >
>
>


-- 
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/20120225/afc3c9bd/attachment.html>