<div dir="ltr"><div class="gmail_extra"><div class="gmail_quote">On Wed, May 28, 2014 at 12:06 PM, Lawrence Mitchell <span dir="ltr"><<a href="mailto:lawrence.mitchell@imperial.ac.uk" target="_blank">lawrence.mitchell@imperial.ac.uk</a>></span> wrote:<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">Dear all,<br>
<br>
a little context, I have a 3x3 block system that comes from a hybridised<br>
discretisation of a helmholtz operator, so I know both my "velocity" and<br>
"pressure" spaces are discontinuous and hence easy to invert. I'd like<br>
to try this by doing schur complements recursively, eliminating both<br>
variables onto the lagrange multiplier one after the other. Anyway,<br>
onto the problems:<br>
<br>
I assemble mixed systems by decomposing them into the individual pieces,<br>
building the individual operators and then gluing them together again in<br>
a MatNest (yes, yes, I know Jed says I should be using<br>
MatGetLocalSubMatrix and friends, but anyway). I have a DMPlex which<br>
I'm using to provide me with sections for the individual fields, so I'd<br>
like to be able to tell the DM about the glued together fields so I can<br>
hang it on a SNES for fieldsplit preconditioning. However, AFAICT,<br>
there seems to be no way for the default global section to deal with<br>
non-interleaved fields, which means that the ISes that the DM creates to<br>
describe the fields are interleaved, and hence don't look like the ises<br>
my MatNest has. Am I completely barking up the wrong tree here?<br></blockquote><div><br></div><div>The abstraction that is breaking here is the (dof, offset) storage of the global</div><div>indices. Here is my sketch of how this could be made to work. You could have</div>
<div>a replacement for DMCreateDefaultGlobalSection() which put valid global</div><div>indices in the field subsections, but the top-level global section would not be</div><div>valid. These indices could mirror your MatNest indices.</div>
<div><br></div><div>The problem is that I do not know what other functions would break, since a</div><div>bunch of things use the global indices. If all you care about is FS, then this</div><div>might work since FS pulls out subdms for the fields, and this operation could</div>
<div>be made to respect the global indices in the subfields (I think).</div><div><br></div><div>However, to me, this sounds like an equivalent amount of work to writing a</div><div>map from the separated to the interlaced space using the l2g for your</div>
<div>submatrices.</div><div><br></div><div> Matt</div><div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
I note in passing that it somehow seems unlikely that this will<br>
currently work for non 2x2 systems since MatNestFindIS can only cope<br>
with pulling out single blocks (rather than say a nest of a 2x2 block).<br>
<br>
This is a sketch of what I'm doing:<br>
<br>
import petsc4py<br>
import sys<br>
petsc4py.init(sys.argv)<br>
from petsc4py import PETSc<br>
<br>
<br>
bdy = PETSc.DMPlex().create()<br>
bdy.setDimension(1)<br>
bdy.createSquareBoundary([0, 0], [1, 1], [1, 1])<br>
<br>
dm = PETSc.DMPlex().generate(bdy)<br>
<br>
sec = dm.createSection(1, [1], [0, 0, 1])<br>
<br>
gsec = PETSc.Section().create()<br>
nfield = 2<br>
gsec.setNumFields(nfield)<br>
gsec.setChart(*sec.getChart())<br>
<br>
pStart, pEnd = gsec.getChart()<br>
offset = 0<br>
for p in range(pStart, pEnd):<br>
sdof = sec.getDof(p)<br>
if sdof > 0:<br>
for i in range(nfield):<br>
gsec.setFieldDof(p, i, sdof)<br>
# this needs the obvious change in petsc4py to work<br>
gsec.setFieldOffset(p, i, offset)<br>
offset += sdof<br>
gsec.setDof(p, sdof * nfield)<br>
<br>
for i in range(nfield):<br>
gsec.setFieldName(i, "%s" % i)<br>
<br>
dm.setDefaultSection(gsec)<br>
<br>
mat = dm.createMat()<br>
for i in range(2 * nfield):<br>
# mat.getSubMatrix<br>
mat.setValues(i, i, i/2 + 1)<br>
mat.assemble()<br>
<br>
# This works, since I got the Mat from the DM<br>
ksp = PETSc.KSP().create()<br>
ksp.setDM(dm)<br>
ksp.setDMActive(False)<br>
ksp.setOperators(mat)<br>
b = mat.createVecLeft()<br>
ksp.setUp()<br>
ksp.setFromOptions()<br>
for i in range(2 * nfield):<br>
b.setValues(i, (i/2 + 1))<br>
b.assemble()<br>
<br>
x = mat.createVecRight()<br>
<br>
ksp.solve(b, x)<br>
x.view()<br>
<br>
<br>
# This doesn't.<br>
m00 = PETSc.Mat().create()<br>
m00.setSizes((2, 2))<br>
m00.setUp()<br>
m00.setValues(0, 0, 1)<br>
m00.setValues(1, 1, 1)<br>
m00.assemble()<br>
<br>
zero = m00.duplicate()<br>
zero.assemble()<br>
m11 = m00.duplicate()<br>
m11.setValues(0, 0, 2)<br>
m11.setValues(1, 1, 2)<br>
m11.assemble()<br>
m22 = m00.duplicate()<br>
m22.setValues(0, 0, 3)<br>
m22.setValues(1, 1, 3)<br>
m22.assemble()<br>
<br>
<br>
if nfield == 2:<br>
mat = PETSc.Mat().createNest([[m00, zero.copy()],<br>
[zero.copy(), m11]])<br>
elif nfield == 3:<br>
mat = PETSc.Mat().createNest([[m00, zero.copy(), zero.copy()],<br>
[zero.copy(), m11, zero.copy()],<br>
[zero.copy(), zero.copy(), m22]])<br>
else:<br>
raise RuntimeError<br>
<br>
ksp.setOperators(mat)<br>
b = mat.createVecLeft()<br>
ksp.setUp()<br>
ksp.setFromOptions()<br>
for i in range(2 * nfield):<br>
b.setValues(i, (i/2 + 1))<br>
b.assemble()<br>
<br>
x = mat.createVecRight()<br>
<br>
ksp.solve(b, x)<br>
x.view()<br>
<br>
$ python foo.py -pc_type fieldsplit -pc_fieldsplit_type schur<br>
-pc_fieldsplit_schur_fact_type diag<br>
Vec Object: 1 MPI processes<br>
type: seq<br>
1<br>
1<br>
1<br>
1<br>
Traceback (most recent call last):<br>
File "petsc-bar.py", line 91, in <module><br>
ksp.setUp()<br>
File "KSP.pyx", line 336, in petsc4py.PETSc.KSP.setUp<br>
(src/petsc4py.PETSc.c:132758)<br>
petsc4py.PETSc.Error: error code 75<br>
[0] KSPSetUp() line 305 in<br>
/data/lmitche1/src/petsc/src/ksp/ksp/interface/itfunc.c<br>
[0] PCSetUp() line 902 in<br>
/data/lmitche1/src/petsc/src/ksp/pc/interface/precon.c<br>
[0] PCSetUp_FieldSplit() line 580 in<br>
/data/lmitche1/src/petsc/src/ksp/pc/impls/fieldsplit/fieldsplit.c<br>
[0] MatGetSubMatrix() line 7345 in<br>
/data/lmitche1/src/petsc/src/mat/interface/matrix.c<br>
[0] MatGetSubMatrix_Nest() line 387 in<br>
/data/lmitche1/src/petsc/src/mat/impls/nest/matnest.c<br>
[0] MatNestFindSubMat() line 371 in<br>
/data/lmitche1/src/petsc/src/mat/impls/nest/matnest.c<br>
[0] MatNestFindIS() line 310 in<br>
/data/lmitche1/src/petsc/src/mat/impls/nest/matnest.c<br>
[0] Arguments are incompatible<br>
[0] Could not find index set<br>
<br>
If I break in MatNestFindIS I see:<br>
<br>
(gdb) print ISView(is, PETSC_VIEWER_STDOUT_WORLD)<br>
IS Object: 1 MPI processes<br>
type: general<br>
Number of indices in set 2<br>
0 0<br>
1 2<br>
$1 = 0<br>
<br>
i.e. the DM has set up field 0 as living at dofs 0 and 2 (rather than<br>
what I want, 0 and 1).<br>
<br>
</blockquote></div><br><br clear="all"><div><br></div>-- <br>What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br>
-- Norbert Wiener
</div></div>