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.
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)
DMFieldSplittingAddField(dm,fscons,"momentum",{1,2,3}); becomes
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.
>
