<br><br><div class="gmail_quote">On Fri, Feb 10, 2012 at 7:10 PM, Barry Smith <span dir="ltr"><<a href="mailto:bsmith@mcs.anl.gov">bsmith@mcs.anl.gov</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">

<div class="im"><br>
On Feb 10, 2012, at 7:04 PM, Dmitry Karpeev wrote:<br>
<br>
><br>
><br>
> On Fri, Feb 10, 2012 at 7:02 PM, Barry Smith <<a href="mailto:bsmith@mcs.anl.gov">bsmith@mcs.anl.gov</a>> wrote:<br>
><br>
> On Feb 10, 2012, at 6:45 PM, Dmitry Karpeev wrote:<br>
><br>
> ><br>
> ><br>
> > On Fri, Feb 10, 2012 at 5:24 PM, Barry Smith <<a href="mailto:bsmith@mcs.anl.gov">bsmith@mcs.anl.gov</a>> wrote:<br>
> ><br>
> > On Feb 10, 2012, at 5:14 PM, Jed Brown wrote:<br>
> ><br>
> > > On Fri, Feb 10, 2012 at 17:08, Dmitry Karpeev <<a href="mailto:karpeev@mcs.anl.gov">karpeev@mcs.anl.gov</a>> wrote:<br>
> > > I think it's completely natural for a DM to assemble two operators -- the discretizations for the two are likely to be related anyway -- as soon as we decide that it's natural for KSP to take in two matrices and, more importantly, for the callback set with DMSetJacobian() to compute two matrices: if a DM knows how to compute  two "Jacobians",<br>


> ><br>
> >    DM doesn't/shouldn't know how to compute two Jacobians. Where did you get that from?<br>
> ><br>
> >  Here's the current declaration of DMSetJacobian:<br>
> >  PetscErrorCode  DMSetJacobian(DM dm,PetscErrorCode (*f)(DM,Vec,Mat,Mat,MatStructure*))<br>
> > This attaches (to a single DM) a callback that computes *two* matrices, just like SNES/KSP would want.<br>
><br>
>   You are right this is odd now and eventually needs some kind of fixing. I kept the two Mats to match the SNESSetJacobian(). Of course we cannot really use it since a DM can only give back one matrix. Shouldn't this be changed to the function f taking two DMs?<br>


><br>
> I think SNES should make two callbacks -- one for each DM -- to have the "user" (or an FEM library) compute J and J_pre.<br>
<br>
</div>   That is what I was thinking originally. But Jed pointed out that sometimes the work between computing the two Jacobians may be mostly common so there should be some way to allow the user to provide code that does both together and not have the two callbacks as you list above.<br>


<br>
    We'll never have a perfect design.<br></blockquote><div><br></div><div>How can there be a single callback, when the two DMs could be of different types?</div><div>I guess it could be done using the double-dispatch mechanism: </div>

<div>DMComputeOperators(DM dmJ, DM dmP, Mat *J, Mat *P) would dispatch based on the types of the two DMs, </div><div>and, if no such method exists, defaulting internally to calling</div><div>DMComputeMatrix(DM dm, Mat *A) for both J and P.</div>

<div>I still think that DMComputeMatrix() should be different from DMCreateMatrix() -- the latter just preallocates --</div><div>eliminating the need for DMSetPreallocateOnly().</div><div><br></div><div>If the two DMs are of the same type and, say, share a mesh, then the dual dispatch can take advantage of</div>

<div>a corresponding implementation.</div><div><br></div><div><br></div><div>Dmitry.</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>
   Barry<br>
</font></span><div class="HOEnZb"><div class="h5">><br>
> Dmitry.<br>
><br>
>   Barry<br>
><br>
> > Dmitry.<br>
> ><br>
> > > why wouldn't it know how to create/preallocate the two corresponding matrices?<br>
> ><br>
> ><br>
> >   The reason I don't like having the single DM doing this is that one could use very different beasties to do the true Jacobian and preconditioner (not just a simple stencil change) and shoving that stuff into a single DM is unnatural. For example true problem is on an unstructured grid, preconditioner problem on a simplier structured grid.<br>


> ><br>
> ><br>
> > ><br>
> > > Lots of functions get messy if the DM has multiple ways to do something. Should DMCreateLocalVector() use the full stencil or the preconditioning stencil, what should DMGlobalToLocalBegin() be updating, etc.<br>


> > ><br>
> > > Barry's solution of having separate DMs sounds cleaner to me, at least modulo needing conventions about which DM on which to PetscObjectCompose() things needed by certain callbacks (e.g. in the FAS with TS stuff I'm doing).<br>


> ><br>
> >   Clean that up and provide a formal way to do those callbacks that don't use kludgy PetscObjectCompose() or PetscObjectComposeFunction()<br>
> ><br>
> >   Barry<br>
> ><br>
> > Yes, it will be slightly painful to pass two DM into the KSP and PC and track them in multigrid but we do do that with the two operators already so using two DM seems very natural.<br>
> ><br>
><br>
><br>
<br>
</div></div></blockquote></div><br>