<div dir="ltr"><div class="gmail_extra"><div class="gmail_quote">On Fri, Jan 26, 2018 at 8:39 PM, Stephen Wornom <span dir="ltr"><<a href="mailto:stephen.wornom@inria.fr" target="_blank">stephen.wornom@inria.fr</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div><div style="font-family:"times new roman","new york",times,serif;font-size:18pt;color:rgb(0,0,0)"><div>I am disparate to implement the petsc library gmres in an CFD code and need a little help to get going.<br></div><div><br></div><div>23 years ago, gmres was implemented/written in the code by Y. Saad, modified by A. Malevsky, version February 1, 1995.</div><div><br></div><div>petsc, most likely started with the same version, added new solvers, ...etc and maintained and corrected bugs reported by users<br></div><div>and I would like to use the latest petsc version of gmres<br></div><div><br></div><div>Why?<br></div><div>Normally we have no convergence problems using gmres.<br></div><div>However sometimes gmres diverges instantaneously, no warning, just a negative density.<br></div><div>Our code is parallel mpi, unstructured, uses a low Mach preconditioner, FV Roe scheme.<br></div><div><br></div><div>Thus I would like to add an petsc option.<br></div><div>Would someone outline, in detail, the steps that I need to follow?<br></div><div>I need all the encourage that I can get to add the option in the code.<br></div></div></div></blockquote><div><br></div><div>I normally divide this process into four steps:</div><div><br></div><div> 1) Call PetscInitialize() and PetscFinalize() once in the code.</div><div><br></div><div> We usually do this at the beginning and the end of the code, but you can do it anywhere you want as long as</div><div> Init precedes any PETSc calls, and Fin follows them.</div><div><br></div><div> 2) Convert the input vectors to PETSc Vecs</div><div><br></div><div> This should be easy. We have wrapper constructors</div><div><br></div><div> <a href="http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Vec/VecCreateSeqWithArray.html">http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Vec/VecCreateSeqWithArray.html</a></div><div> <a href="http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Vec/VecCreateMPIWithArray.html#VecCreateMPIWithArray">http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Vec/VecCreateMPIWithArray.html#VecCreateMPIWithArray</a><br></div><div><br></div><div>There is an example of using these in Fortran90 to create a solver interface that takes raw arrays</div><div><br></div><div> <a href="http://www.mcs.anl.gov/petsc/petsc-current/src/ksp/ksp/examples/tutorials/ex13f90.F90.html">http://www.mcs.anl.gov/petsc/petsc-current/src/ksp/ksp/examples/tutorials/ex13f90.F90.html</a></div><div><br></div><div> 3) Convert the system matrix to a PETSc Mat</div><div><br></div><div> If you do not want to use PETSc preconditioners, you can easily do this by wrapping just the action of the matrix:</div><div><br></div><div> <a href="http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MatShellSetOperation.html">http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MatShellSetOperation.html</a></div><div><br></div><div> and there are Fortran examples listed at the bottom of that page. If instead you want to enable preconditioners, then</div><div> you must put the values into a PETSc Mat structure, and you would use MatSetValues() in the naive way for that.</div><div><br></div><div> 4) Create a solver, pass it the system matrix, rhs, and solution</div><div><br></div><div> There are many examples of this. For instance</div><div><br></div><div> <a href="http://www.mcs.anl.gov/petsc/petsc-current/src/ksp/ksp/examples/tutorials/ex1f.F90.html">http://www.mcs.anl.gov/petsc/petsc-current/src/ksp/ksp/examples/tutorials/ex1f.F90.html</a></div><div><br></div><div> Thanks,</div><div><br></div><div> Matt</div><div> </div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div><div style="font-family:"times new roman","new york",times,serif;font-size:18pt;color:rgb(0,0,0)"><div></div><div>petsc is installed on our computers<br></div><div>Thanks in advance,<br></div><div>Stephen<br></div><div><br></div><div><br></div><div> IF ( petsc .EQ. 1 ) THEN<br>c petsc CALL MatAssembleBegin<br>c petsc CALL MatAssembleEnd</div><div>... add steps<br> ELSE<br> CALL GMRESASR ! from our code<br> ENDIF<br><div><br></div></div><div><br></div><div><br></div><div>This is from our code which may be useful information<br></div><div>c-----------------------------<wbr>------------------------------<wbr>------------<br>c flexible GMRES routine. This is a version of GMRES which allows a <br>c a variable preconditioner. Implemented with a reverse communication <br>c protocole for flexibility -<br>c DISTRIBUTED VERSION (USES DISTDOT FOR DDOT) <br>c explicit (exact) residual norms for restarts <br>c written by Y. Saad, modified by A. Malevsky, version February 1, 1995<br>c-----------------------------<wbr>------------------------------<wbr>------------<br>c This Is A Reverse Communication Implementation. <br>c-----------------------------<wbr>-------------------- <br>c USAGE: (see also comments for icode below). CGMRES<br>c should be put in a loop and the loop should be active for as<br>c long as icode is not equal to 0. On RETURN fgmres will<br>c 1) either be requesting the new preconditioned vector applied<br>c to wk1 in case icode.eq.1 (result should be put in wk2) <br>c 2) or be requesting the product of A applied to the vector wk1<br>c in case icode.eq.2 (result should be put in wk2) <br>c 3) or be terminated in case icode .eq. 0. <br>c on entry always set icode = 0. So icode should be set back to zero<br>c upon convergence.<br>c-----------------------------<wbr>------------------------------<wbr>------------<br>c Here is a typical way of running fgmres: <br>c<br>c icode = 0<br>c 1 continue<br>c CALL fgmres (n,im,rhs,sol,i,vv,w,wk1, wk2,eps,maxits,iout,icode)<br>c<br>c if (icode .eq. 1) then<br>c CALL precon(n, wk1, wk2) <--- user's variable preconditioning<br>c goto 1<br>c else if (icode .ge. 2) then<br>c CALL matvec (n,wk1, wk2) <--- user's matrix vector product. <br>c goto 1<br>c else <br>c ----- done ---- <br>c .........<br> </div></div></div></blockquote></div><br><br clear="all"><div><br></div>-- <br><div class="gmail_signature"><div dir="ltr"><div><div dir="ltr"><div>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><br></div><div><a href="http://www.caam.rice.edu/~mk51/" target="_blank">https://www.cse.buffalo.edu/~knepley/</a><br></div></div></div></div></div>
</div></div>