<div class="gmail_extra">On Thu, May 3, 2012 at 1:28 PM, Pavanakumar Mohanamuraly <span dir="ltr"><<a href="mailto:m.pavanakumar@gmail.com" target="_blank">m.pavanakumar@gmail.com</a>></span> wrote:<br><div class="gmail_quote">
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div class="gmail_quote"><div>Dear Matt,</div><div><br></div><div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
> Currently I am planning to use Petsc Krylov solver for implicit time<br>
> integration of my CFD code. I need some pointer on the implementation. I am<br>
> providing a detailed list of things I am not very clear with in both the<br>
> documentation and the mailing list archives.<br>
><br>
> --> The equations to be solved is of the form : [ I / \delta t + J ] \delta<br>
> U^n = - R[ U^n ]<br>
><br>
> --> The Jacobian J is hard to evaluate exactly and hence we use the finite<br>
> difference method to evaluate J. I want to use a matrix-free approach so<br>
> all I need is the action of this J on some vector v.<br>
><br>
> --> Jv = \frac{\partial R}{\partial U}v = ( R[ U^n + hv ] - R[ U^] ) / h,<br>
> h is some parameter. Since I is identity matrix and \delta t is fixed this<br>
> diagonal matrix is trivial to evaluate. Thus [ I / \delta t + J ] is<br>
> available as a user define function in my solver.<br>
><br>
> 1) Vector U^n is from an unstructured node distribution and is partitioned<br>
> outside of Petsc. I have both the global indexing and local indexing<br>
> information along with the ghost node information.<br>
> 2) I have already defined functions to move data across the ghost nodes<br>
> given the local variable address and memory stride/offset information.<br>
> 3) In Jacobian-vector Jv evaluation, one has to exchange the ghost cell<br>
> values of vector v with adjacent processor before the calculation is done.<br>
> 4) I can write the Jacobian-vector Jv evaluation function to perform ghost<br>
> node exchange before evaluation but I am not clear as to how I can<br>
> interface my local arrays and local / global indexing information into<br>
> Petsc.<br>
> 5) I understand that I have to create a Matrix Shell for the matrix free<br>
> method and define my user-defined function to evaluate the matrix-vector<br>
> product given an input vector. But it is not clear as to how Petsc<br>
> understands how I order my vectors locally and how that corresponds to the<br>
> global indexes.<br>
> 6) I cannot use the DA object as it is for structured/contiguous<br>
> partitioning .i.e., the local indices don't correspond to contiguous global<br>
> ordering/indexing.<br>
> 7) I also don't want to use the unstructured mesh and decomposition<br>
> routines in Petsc and I already have my own.<br>
><br>
> Can this be done in Petsc ?<br>
><br>
<br>
Yes, there are ??? parts:<br>
<br>
1) You want to manage unstructured data yourself. The easiest way to do<br>
this now is to implement a DM. You just need to provide<br>
routines for making global and local vectors (you have this), and<br>
routines for mapping between them (you have this). You will have<br>
to use petsc-dev, but this definitely looks like the easiest thing to<br>
do.<br>
<br></blockquote><div><br></div><div><br></div><div><div>Thanks for the reply. </div><div><br></div><div> </div></div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
a) If you consider using part of our infrastructure, you can build a<br>
PetscSection and a PetscSF, and everything else will be automatically<br>
calculated. These are very new things, so there is not a lot of<br>
documentation, and thus you might prefer to do the whole DM yourself.<br>
<br></blockquote><div><br></div><div><div><div>Most of the tutorial slides point out that DM is for structured arrays (maybe this mislead me). Can you point me to the functions that i would have to use to create </div></div>
</div></div></blockquote><div><br></div><div>No, the DA is for structured grids. It is a subclass of DM. There is also a DM for unstructured grids.</div><div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
<div class="gmail_quote"><div><div><div>1) DM array ghost cells</div></div></div></div></blockquote><div><br></div><div><a href="http://www.mcs.anl.gov/petsc/petsc-dev/docs/manualpages/DM/DMCreateGlobalVector.html">http://www.mcs.anl.gov/petsc/petsc-dev/docs/manualpages/DM/DMCreateGlobalVector.html</a></div>
<div><a href="http://www.mcs.anl.gov/petsc/petsc-dev/docs/manualpages/DM/DMCreateLocalVector.html">http://www.mcs.anl.gov/petsc/petsc-dev/docs/manualpages/DM/DMCreateLocalVector.html</a></div><div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
<div class="gmail_quote"><div><div><div>2) Local/Global indexing</div></div></div></div></blockquote><div><a href="http://www.mcs.anl.gov/petsc/petsc-dev/docs/manualpages/DM/DMGetLocalToGlobalMapping.html">http://www.mcs.anl.gov/petsc/petsc-dev/docs/manualpages/DM/DMGetLocalToGlobalMapping.html</a></div>
<div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div class="gmail_quote"><div><div><div>3) DM array exchange</div></div></div></div></blockquote><div><br></div>
<div><br class="Apple-interchange-newline"><a href="http://www.mcs.anl.gov/petsc/petsc-dev/docs/manualpages/DM/DMGlobalToLocalBegin.html">http://www.mcs.anl.gov/petsc/petsc-dev/docs/manualpages/DM/DMGlobalToLocalBegin.html</a></div>
<div><a href="http://www.mcs.anl.gov/petsc/petsc-dev/docs/manualpages/DM/DMGlobalToLocalEnd.html">http://www.mcs.anl.gov/petsc/petsc-dev/docs/manualpages/DM/DMGlobalToLocalEnd.html</a></div><div><a href="http://www.mcs.anl.gov/petsc/petsc-dev/docs/manualpages/DM/DMLocalToGlobalBegin.html">http://www.mcs.anl.gov/petsc/petsc-dev/docs/manualpages/DM/DMLocalToGlobalBegin.html</a></div>
<div><a href="http://www.mcs.anl.gov/petsc/petsc-dev/docs/manualpages/DM/DMLocalToGlobalEnd.html">http://www.mcs.anl.gov/petsc/petsc-dev/docs/manualpages/DM/DMLocalToGlobalEnd.html</a></div><div><br></div><div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
<div class="gmail_quote"><div><div><div> </div></div></div><div>It would be nice you can share some code snippet on writing our own DM's. </div></div></blockquote><div><br></div><div>You can look at any of the existing DMs in src/dm/impls</div>
<div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div class="gmail_quote"><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
2) You want to do matrix-free action of your Jacobian. This is easy using<br>
MatShell, however the real question is How will you precondition that<br>
equation?</blockquote><div><br></div><div>Currently I use the Jacobi diagonal preconditioner. I use the largest spectral radius of my Euler equations \lambda = U + a. BTW I actually have implemented my own parallel Krylov solver for my CFD code using the algorithm in Yousuf Saad's book. But my advisor urged me to use Petsc so that I would be able to use the various pre-conditioners in Petsc and the code would run more faster/optimal. I have not written my Krylov solver optimally.</div>
<div><br></div><div>So is there is no way I can use preconditioning in matrix-free context in Petsc ?</div></div></blockquote><div><br></div><div>It depends what kind of preconditioner you want to use. Not even Jacobi is matrix-free. You need the diagonal. You</div>
<div>can basically use Chebychev. What you should consider is providing a simplified matrix to use for building a</div><div>preconditioner.</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">
<div class="gmail_quote"><div>Regards,</div><div> </div></div>--<div> <br>Pavanakumar Mohanamuraly<br>
</div>
</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<br>
</div>