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

[Matthew Knepley]

On Wed, May 28, 2014 at 12:06 PM, Lawrence Mitchell <
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).

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.

   Matt


> I note in passing that it somehow seems unlikely that this will
> currently work for non 2x2 systems since MatNestFindIS can only cope
> with pulling out single blocks (rather than say a nest of a 2x2 block).
>
> This is a sketch of what I'm doing:
>
> import petsc4py
> import sys
> petsc4py.init(sys.argv)
> from petsc4py import PETSc
>
>
> bdy = PETSc.DMPlex().create()
> bdy.setDimension(1)
> bdy.createSquareBoundary([0, 0], [1, 1], [1, 1])
>
> dm = PETSc.DMPlex().generate(bdy)
>
> sec = dm.createSection(1, [1], [0, 0, 1])
>
> gsec = PETSc.Section().create()
> nfield = 2
> gsec.setNumFields(nfield)
> gsec.setChart(*sec.getChart())
>
> pStart, pEnd = gsec.getChart()
> offset = 0
> for p in range(pStart, pEnd):
>     sdof = sec.getDof(p)
>     if sdof > 0:
>         for i in range(nfield):
>             gsec.setFieldDof(p, i, sdof)
>             # this needs the obvious change in petsc4py to work
>             gsec.setFieldOffset(p, i, offset)
>             offset += sdof
>         gsec.setDof(p, sdof * nfield)
>
> for i in range(nfield):
>     gsec.setFieldName(i, "%s" % i)
>
> dm.setDefaultSection(gsec)
>
> mat = dm.createMat()
> for i in range(2 * nfield):
>     # mat.getSubMatrix
>     mat.setValues(i, i, i/2 + 1)
> mat.assemble()
>
> # This works, since I got the Mat from the DM
> ksp = PETSc.KSP().create()
> ksp.setDM(dm)
> ksp.setDMActive(False)
> ksp.setOperators(mat)
> b = mat.createVecLeft()
> ksp.setUp()
> ksp.setFromOptions()
> for i in range(2 * nfield):
>     b.setValues(i, (i/2 + 1))
> b.assemble()
>
> x = mat.createVecRight()
>
> ksp.solve(b, x)
> x.view()
>
>
> # This doesn't.
> m00 = PETSc.Mat().create()
> m00.setSizes((2, 2))
> m00.setUp()
> m00.setValues(0, 0, 1)
> m00.setValues(1, 1, 1)
> m00.assemble()
>
> zero = m00.duplicate()
> zero.assemble()
> m11 = m00.duplicate()
> m11.setValues(0, 0, 2)
> m11.setValues(1, 1, 2)
> m11.assemble()
> m22 = m00.duplicate()
> m22.setValues(0, 0, 3)
> m22.setValues(1, 1, 3)
> m22.assemble()
>
>
> if nfield == 2:
>     mat = PETSc.Mat().createNest([[m00, zero.copy()],
>                                   [zero.copy(), m11]])
> elif nfield == 3:
>     mat = PETSc.Mat().createNest([[m00, zero.copy(), zero.copy()],
>                                   [zero.copy(), m11, zero.copy()],
>                                   [zero.copy(), zero.copy(), m22]])
> else:
>     raise RuntimeError
>
> ksp.setOperators(mat)
> b = mat.createVecLeft()
> ksp.setUp()
> ksp.setFromOptions()
> for i in range(2 * nfield):
>     b.setValues(i, (i/2 + 1))
> b.assemble()
>
> x = mat.createVecRight()
>
> ksp.solve(b, x)
> x.view()
>
> $ python foo.py -pc_type fieldsplit -pc_fieldsplit_type schur
> -pc_fieldsplit_schur_fact_type diag
> Vec Object: 1 MPI processes
>   type: seq
> 1
> 1
> 1
> 1
> Traceback (most recent call last):
>   File "petsc-bar.py", line 91, in <module>
>     ksp.setUp()
>   File "KSP.pyx", line 336, in petsc4py.PETSc.KSP.setUp
> (src/petsc4py.PETSc.c:132758)
> petsc4py.PETSc.Error: error code 75
> [0] KSPSetUp() line 305 in
> /data/lmitche1/src/petsc/src/ksp/ksp/interface/itfunc.c
> [0] PCSetUp() line 902 in
> /data/lmitche1/src/petsc/src/ksp/pc/interface/precon.c
> [0] PCSetUp_FieldSplit() line 580 in
> /data/lmitche1/src/petsc/src/ksp/pc/impls/fieldsplit/fieldsplit.c
> [0] MatGetSubMatrix() line 7345 in
> /data/lmitche1/src/petsc/src/mat/interface/matrix.c
> [0] MatGetSubMatrix_Nest() line 387 in
> /data/lmitche1/src/petsc/src/mat/impls/nest/matnest.c
> [0] MatNestFindSubMat() line 371 in
> /data/lmitche1/src/petsc/src/mat/impls/nest/matnest.c
> [0] MatNestFindIS() line 310 in
> /data/lmitche1/src/petsc/src/mat/impls/nest/matnest.c
> [0] Arguments are incompatible
> [0] Could not find index set
>
> If I break in MatNestFindIS I see:
>
> (gdb) print ISView(is, PETSC_VIEWER_STDOUT_WORLD)
> IS Object: 1 MPI processes
>   type: general
> Number of indices in set 2
> 0 0
> 1 2
> $1 = 0
>
> i.e. the DM has set up field 0 as living at dofs 0 and 2 (rather than
> what I want, 0 and 1).
>
>


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