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

[Lawrence Mitchell]

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?

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