[petsc-dev] DMPlex sections with non-interleaved fields

[Matthew Knepley]

On Wed, May 28, 2014 at 12:50 PM, Lawrence Mitchell <
lawrence.mitchell at imperial.ac.uk> wrote:

> On 28/05/14 18:24, Matthew Knepley wrote:
> > On Wed, May 28, 2014 at 12:06 PM, Lawrence Mitchell
> > <lawrence.mitchell at imperial.ac.uk
> > <mailto:lawrence.mitchell at imperial.ac.uk>> wrote:
> >
> >     Dear all,
> >
> >     a little context, I have a 3x3 block system that comes from a
> hybridised
> >     discretisation of a helmholtz operator, so I know both my "velocity"
> and
> >     "pressure" spaces are discontinuous and hence easy to invert.  I'd
> like
> >     to try this by doing schur complements recursively, eliminating both
> >     variables onto the lagrange multiplier one after the other.  Anyway,
> >     onto the problems:
> >
> >     I assemble mixed systems by decomposing them into the individual
> pieces,
> >     building the individual operators and then gluing them together
> again in
> >     a MatNest (yes, yes, I know Jed says I should be using
> >     MatGetLocalSubMatrix and friends, but anyway).  I have a DMPlex which
> >     I'm using to provide me with sections for the individual fields, so
> I'd
> >     like to be able to tell the DM about the glued together fields so I
> can
> >     hang it on a SNES for fieldsplit preconditioning.  However, AFAICT,
> >     there seems to be no way for the default global section to deal with
> >     non-interleaved fields, which means that the ISes that the DM
> creates to
> >     describe the fields are interleaved, and hence don't look like the
> ises
> >     my MatNest has.  Am I completely barking up the wrong tree here?
> >
> >
> > The abstraction that is breaking here is the (dof, offset) storage of
> > the global
> > indices. Here is my sketch of how this could be made to work. You could
> have
> > a replacement for DMCreateDefaultGlobalSection() which put valid global
> > indices in the field subsections, but the top-level global section would
> > not be
> > valid. These indices could mirror your MatNest indices.
> >
> > The problem is that I do not know what other functions would break,
> since a
> > bunch of things use the global indices. If all you care about is FS,
> > then this
> > might work since FS pulls out subdms for the fields, and this operation
> > could
> > be made to respect the global indices in the subfields (I think).
>
> Looking at DMCreateSubDM_Section_Private (which does the hard work of
> building the ISes) that looks possible, yes.
>
> So I only need this for FS, however I guess if you expose the API,
> someone else may want it for other things.  I still think (reading the
> code at least) that MatNestFindIS would break for this use case if I
> have (say) a 3x3 block system and want to do -pc_fieldsplit_0_fields 0,1
> ....
>
> > However, to me, this sounds like an equivalent amount of work to writing
> a
> > map from the separated to the interlaced space using the l2g for your
> > submatrices.
>
> So in this case, I would have a monolithic (interleaved) matrix which I
> assemble into by doing MatGetLocalSubMatrix to get a block.  Rather than
> having a Nest for which I do MatNestGetSubMat to get a block.  Is that
> right?  So in this case I just need to arrange to create my monolithic
> matrices and the ISes to pass to MatGetLocalSubMatrix (note I can't let
> the DM create the Mat because in some cases I'm lying to it about the
> connectivity so the sparsity pattern it would give is wrong, even if the
> number of dofs and offsets are right).  However, is it intended that
> this kind of setup ever work with DMSetMatType(..., MATNEST)?  If so, I
> think I'll be back at square one with this approach.


Yes, exactly. You would just use your ISes from the section to pull out the
submatrix, and it should then function exactly as the one from
MatNestGetSubMat().
I don't think square one here, because the approach is more flexible. If
the IS
matches the IS for a Nest block, then the submatrix can just be returned
wholesale.
If the Mat is not a Nest, or the IS is interleaved, then the interface is
the same.
I think this approach can do what you want.

  Thanks,

     Matt


>
> Lawrence
>

-- 
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/20140528/7641db06/attachment.html>