[petsc-users] Novice: How to implement the petsc Gmres libraries in my CFD code

[Matthew Knepley]

On Fri, Jan 26, 2018 at 8:39 PM, Stephen Wornom <stephen.wornom at inria.fr>
wrote:

> I am disparate to implement  the petsc library gmres in an CFD code and
> need a little help to get going.
>
> 23 years ago, gmres was implemented/written in the code by Y. Saad,
> modified by A. Malevsky, version February 1, 1995.
>
> petsc, most likely started with the same version, added new solvers,
> ...etc and maintained and corrected bugs reported by users
> and I would like to use the latest petsc version of gmres
>
> Why?
> Normally we have no convergence problems using gmres.
> However sometimes gmres diverges instantaneously, no warning, just a
> negative density.
> Our code is parallel mpi, unstructured, uses a low Mach preconditioner, FV
> Roe scheme.
>
> Thus I would like to add an petsc option.
> Would  someone outline, in detail, the steps that I need to follow?
> I need all the encourage that I can get to add the option in the code.
>

I normally divide this process into four steps:

  1) Call PetscInitialize() and PetscFinalize() once in the code.

  We usually do this at the beginning and the end of the code, but you can
do it anywhere you want as long as
  Init precedes any PETSc calls, and Fin follows them.

  2) Convert the input vectors to PETSc Vecs

  This should be easy. We have wrapper constructors


http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Vec/VecCreateSeqWithArray.html

http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Vec/VecCreateMPIWithArray.html#VecCreateMPIWithArray

There is an example of using these in Fortran90 to create a solver
interface that takes raw arrays


http://www.mcs.anl.gov/petsc/petsc-current/src/ksp/ksp/examples/tutorials/ex13f90.F90.html

  3) Convert the system matrix to a PETSc Mat

  If you do not want to use PETSc preconditioners, you can easily do this
by wrapping just the action of the matrix:


http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MatShellSetOperation.html

  and there are Fortran examples listed at the bottom of that page. If
instead you want to enable preconditioners, then
  you must put the values into a PETSc Mat structure, and you would use
MatSetValues() in the naive way for that.

  4) Create a solver, pass it the system matrix, rhs, and solution

  There are many examples of this. For instance


http://www.mcs.anl.gov/petsc/petsc-current/src/ksp/ksp/examples/tutorials/ex1f.F90.html

  Thanks,

     Matt


> petsc is installed on our computers
> Thanks in advance,
> Stephen
>
>
>             IF ( petsc .EQ. 1 ) THEN
> c petsc        CALL MatAssembleBegin
> c petsc        CALL MatAssembleEnd
> ... add steps
>             ELSE
>                CALL GMRESASR ! from our code
>             ENDIF
>
>
>
> This is from our code which may be useful information
> c-----------------------------------------------------------------------
> c flexible GMRES routine. This is a version of GMRES which allows a
> c a variable preconditioner. Implemented with a reverse communication
> c protocole for flexibility -
> c DISTRIBUTED VERSION (USES DISTDOT FOR DDOT)
> c explicit (exact) residual norms for restarts
> c written by Y. Saad, modified by A. Malevsky, version February 1, 1995
> c-----------------------------------------------------------------------
> c This Is A Reverse Communication Implementation.
> c-------------------------------------------------
> c USAGE: (see also comments for icode below). CGMRES
> c should be put in a loop and the loop should be active for as
> c long as icode is not equal to 0. On RETURN fgmres will
> c    1) either be requesting the new preconditioned vector applied
> c       to wk1 in case icode.eq.1 (result should be put in wk2)
> c    2) or be requesting the product of A applied to the vector wk1
> c       in case icode.eq.2 (result should be put in wk2)
> c    3) or be terminated in case icode .eq. 0.
> c on entry always set icode = 0. So icode should be set back to zero
> c upon convergence.
> c-----------------------------------------------------------------------
> c Here is a typical way of running fgmres:
> c
> c      icode = 0
> c 1    continue
> c      CALL fgmres (n,im,rhs,sol,i,vv,w,wk1, wk2,eps,maxits,iout,icode)
> c
> c      if (icode .eq. 1) then
> c         CALL  precon(n, wk1, wk2)    <--- user's variable preconditioning
> c         goto 1
> c      else if (icode .ge. 2) then
> c         CALL  matvec (n,wk1, wk2)    <--- user's matrix vector product.
> c         goto 1
> c      else
> c         ----- done ----
> c         .........
>
>



-- 
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

https://www.cse.buffalo.edu/~knepley/ <http://www.caam.rice.edu/~mk51/>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20180127/8ca004b2/attachment-0001.html>