[petsc-users] petsc-users Digest, Vol 41, Issue 17

[Matthew Knepley]

On Thu, May 3, 2012 at 1:28 PM, Pavanakumar Mohanamuraly <
m.pavanakumar at gmail.com> wrote:

> Dear Matt,
>
>
>
>> > Currently I am planning to use Petsc Krylov solver for implicit time
>> > integration of my CFD code. I need some pointer on the implementation.
>> I am
>> > providing a detailed list of things I am not very clear with in both the
>> > documentation and the mailing list archives.
>> >
>> > --> The equations to be solved is of the form :  [ I / \delta t + J ]
>> \delta
>> > U^n = - R[ U^n ]
>> >
>> > --> The Jacobian J is hard to evaluate exactly and hence we use the
>> finite
>> > difference method to evaluate J. I want to use a matrix-free approach so
>> > all I need is the action of this J on some vector v.
>> >
>> > --> Jv = \frac{\partial R}{\partial U}v = ( R[ U^n + hv ] - R[ U^] ) /
>> h,
>> > h is some parameter. Since I is identity matrix and \delta t is fixed
>> this
>> > diagonal matrix is trivial to evaluate. Thus [ I / \delta t + J ] is
>> > available as a user define function in my solver.
>> >
>> > 1) Vector U^n is from an unstructured node distribution and is
>> partitioned
>> > outside of Petsc. I have both the global indexing and local indexing
>> > information along with the ghost node information.
>> > 2) I have already defined functions to move data across the ghost nodes
>> > given the local variable address and memory stride/offset information.
>> > 3) In Jacobian-vector Jv evaluation, one has to exchange the ghost cell
>> > values of vector v with adjacent processor before the calculation is
>> done.
>> > 4) I can write the Jacobian-vector Jv evaluation function to perform
>> ghost
>> > node exchange before evaluation but I am not clear as to how I can
>> > interface my local arrays and local / global indexing information into
>> > Petsc.
>> > 5) I understand that I have to create a Matrix Shell for the matrix free
>> > method and define my user-defined function to evaluate the matrix-vector
>> > product given an input vector. But it is not clear as to how Petsc
>> > understands how I order my vectors locally and how that corresponds to
>> the
>> > global indexes.
>> > 6) I cannot use the DA object as it is for structured/contiguous
>> > partitioning .i.e., the local indices don't correspond to contiguous
>> global
>> > ordering/indexing.
>> > 7) I also don't want to use the unstructured mesh and decomposition
>> > routines in Petsc and I already have my own.
>> >
>> > Can this be done in Petsc ?
>> >
>>
>> Yes, there are ??? parts:
>>
>>  1) You want to manage unstructured data yourself. The easiest way to do
>> this now is to implement a DM. You just need to provide
>>      routines for making global and local vectors (you have this), and
>> routines for mapping between them (you have this). You will have
>>      to use petsc-dev, but this definitely looks like the easiest thing to
>> do.
>>
>>
>
> Thanks for the reply.
>
>
>
>>     a) If you consider using part of our infrastructure, you can build a
>> PetscSection and a PetscSF, and everything else will be automatically
>>         calculated. These are very new things, so there is not a lot of
>> documentation, and thus you might prefer to do the whole DM yourself.
>>
>>
> 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
>

No, the DA is for structured grids. It is a subclass of DM. There is also a
DM for unstructured grids.


> 1) DM array ghost cells
>

http://www.mcs.anl.gov/petsc/petsc-dev/docs/manualpages/DM/DMCreateGlobalVector.html
http://www.mcs.anl.gov/petsc/petsc-dev/docs/manualpages/DM/DMCreateLocalVector.html


> 2) Local/Global indexing
>
http://www.mcs.anl.gov/petsc/petsc-dev/docs/manualpages/DM/DMGetLocalToGlobalMapping.html


> 3) DM array exchange
>


http://www.mcs.anl.gov/petsc/petsc-dev/docs/manualpages/DM/DMGlobalToLocalBegin.html
http://www.mcs.anl.gov/petsc/petsc-dev/docs/manualpages/DM/DMGlobalToLocalEnd.html
http://www.mcs.anl.gov/petsc/petsc-dev/docs/manualpages/DM/DMLocalToGlobalBegin.html
http://www.mcs.anl.gov/petsc/petsc-dev/docs/manualpages/DM/DMLocalToGlobalEnd.html



>
> It would be nice you can share some code snippet on writing our own DM's.
>

You can look at any of the existing DMs in src/dm/impls


>  2) You want to do matrix-free action of your Jacobian. This is easy using
>> MatShell, however the real question is How will you precondition that
>>       equation?
>
>
> 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.
>
> So is there is no way I can use preconditioning in matrix-free context in
> Petsc ?
>

It depends what kind of preconditioner you want to use. Not even Jacobi is
matrix-free. You need the diagonal. You
can basically use Chebychev. What you should consider is providing a
simplified matrix to use for building a
preconditioner.

   Matt


> Regards,
>
> --
>
> Pavanakumar Mohanamuraly
>



-- 
What most experimenters take for granted before they begin their
experiments is infinitely more interesting than any results to which their
experiments lead.
-- Norbert Wiener
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20120503/41ee1280/attachment.htm>