[petsc-dev] DM field interfaces

[Barry Smith]

On Feb 25, 2012, at 3:50 PM, Matthew Knepley wrote:

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

   There is nothing "wrong" with doing it this way. It is really only a question of design:

   1) a) Do we have two different ways of indicating subsets of DOFs, sometimes an IS and sometimes a list of integers representing a set of canonical fields or 
       b) do we ALWAYS use an IS and have more IS subclasses for efficient representations of canonical fields. and
   2) Since the above approach only allows defining fields in terms of canonical fields we would also need DMFieldSplittingAddFieldIS(...,IS); for the more general case.
        Note that I also don't like having both MatZeroRows() and MatZeroRowsIS() and think there should just be one and better ISs

   Note that in some sense we are just dealing with C's lack of function overloading.

    All 4 bullets below are true whether we go with "always use IS to indicate fields" approach of mine or the "you can also indicate fields with in terms of canonical fields 
    without creating a canonical IS" model of you and Jed. So what you seem to prefer is pretty much just a question on syntax.

    Note also that the 4 bullets below are what Jed, Dmitry and I have been advocating. 

  Barry


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