<div dir="ltr"><div class="gmail_extra"><div class="gmail_quote">On Wed, May 28, 2014 at 12:50 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"><div class="">On 28/05/14 18:24, Matthew Knepley wrote:<br>
> On Wed, May 28, 2014 at 12:06 PM, Lawrence Mitchell<br>
> <<a href="mailto:lawrence.mitchell@imperial.ac.uk">lawrence.mitchell@imperial.ac.uk</a><br>
</div><div><div class="h5">> <mailto:<a href="mailto:lawrence.mitchell@imperial.ac.uk">lawrence.mitchell@imperial.ac.uk</a>>> wrote:<br>
><br>
> 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>
><br>
><br>
> The abstraction that is breaking here is the (dof, offset) storage of<br>
> the global<br>
> indices. Here is my sketch of how this could be made to work. You could have<br>
> a replacement for DMCreateDefaultGlobalSection() which put valid global<br>
> indices in the field subsections, but the top-level global section would<br>
> not be<br>
> valid. These indices could mirror your MatNest indices.<br>
><br>
> The problem is that I do not know what other functions would break, since a<br>
> bunch of things use the global indices. If all you care about is FS,<br>
> then this<br>
> might work since FS pulls out subdms for the fields, and this operation<br>
> could<br>
> be made to respect the global indices in the subfields (I think).<br>
<br>
</div></div>Looking at DMCreateSubDM_Section_Private (which does the hard work of<br>
building the ISes) that looks possible, yes.<br>
<br>
So I only need this for FS, however I guess if you expose the API,<br>
someone else may want it for other things. I still think (reading the<br>
code at least) that MatNestFindIS would break for this use case if I<br>
have (say) a 3x3 block system and want to do -pc_fieldsplit_0_fields 0,1<br>
....<br>
<div class=""><br>
> However, to me, this sounds like an equivalent amount of work to writing a<br>
> map from the separated to the interlaced space using the l2g for your<br>
> submatrices.<br>
<br>
</div>So in this case, I would have a monolithic (interleaved) matrix which I<br>
assemble into by doing MatGetLocalSubMatrix to get a block. Rather than<br>
having a Nest for which I do MatNestGetSubMat to get a block. Is that<br>
right? So in this case I just need to arrange to create my monolithic<br>
matrices and the ISes to pass to MatGetLocalSubMatrix (note I can't let<br>
the DM create the Mat because in some cases I'm lying to it about the<br>
connectivity so the sparsity pattern it would give is wrong, even if the<br>
number of dofs and offsets are right). However, is it intended that<br>
this kind of setup ever work with DMSetMatType(..., MATNEST)? If so, I<br>
think I'll be back at square one with this approach.</blockquote><div><br></div><div>Yes, exactly. You would just use your ISes from the section to pull out the</div><div>submatrix, and it should then function exactly as the one from MatNestGetSubMat().</div>
<div>I don't think square one here, because the approach is more flexible. If the IS</div><div>matches the IS for a Nest block, then the submatrix can just be returned wholesale.</div><div>If the Mat is not a Nest, or the IS is interleaved, then the interface is the same.</div>
<div>I think this approach can do what you want.</div><div><br></div><div> Thanks,</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">
<span class="HOEnZb"><font color="#888888"><br>
Lawrence<br>
</font></span></blockquote></div><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>