> On May 31, 2023, at 4:16 PM, Kenneth C Hall <kenneth.c.h...@duke.edu> wrote: > > Matt, > > Thanks for your quick reply. I think what you say makes sense. > > You asked what my code does. The MySolver program performs one iteration of a > CFD iteration. The CFD scheme is an explicit scheme that uses multigrid, Mach > number preconditioning, and residual smoothing. Typically, I have to call > MySolver on the order of 40 to 100 times to get acceptable convergence. > > And in fact, I have another version of this Petsc code that uses SNESNGMRES > to solve the problem with MySolver providing the residuals as R = N(x) - x.
So you use SNESSetFunction( a function you wrote that computes R = N(x) - x )? Have you considered using the SNESSetNPC() mechanism? That is, use your solver as a nonlinear preconditioner for any of the PETSc solvers? In particular, you could use it with SNESNRICHARDSON which provides a lineseach, monitor, and convergence testing in a loop. I think this is exactly the loop structure you want, so no need for SNESSHELL. > But I would like a version where I am using just MySolver, without any other > operations applied to it. So I am trying to plug MySolver into the PETSc > system to provide monitoring and other features, and for consistency and > comparison to these other (more appropriate!) uses of PETSc. > > Thanks. > Kenneth > > > From: Matthew Knepley <knep...@gmail.com <mailto:knep...@gmail.com>> > Date: Wednesday, May 31, 2023 at 3:48 PM > To: Kenneth C Hall <kenneth.c.h...@duke.edu <mailto:kenneth.c.h...@duke.edu>> > Cc: petsc-users@mcs.anl.gov <mailto:petsc-users@mcs.anl.gov> > <petsc-users@mcs.anl.gov <mailto:petsc-users@mcs.anl.gov>> > Subject: Re: [petsc-users] Using SNESSHELL as a wrapper for a CFD solver. > > On Wed, May 31, 2023 at 3:21 PM Kenneth C Hall <kenneth.c.h...@duke.edu > <mailto:kenneth.c.h...@duke.edu>> wrote: > Hi, > > I am doing a number of problems using PETSc/SLEPc, but I also work on some > non-PETSc/SLEPc flow solvers. I would like to use PETSc as a wrapper for > this non-PETSc flow solver for compatibility, so I can use the tolerance > monitoring, options, viewers, and for direct comparison to PETSc methods I am > using. > > Here is what I am trying to do… I have a CFD solver that iterates with a > nonlinear iterator of the form x := N(x). This can be expressed in a fortran > routine of the form, > > SUBROUTINE MySolver(x) > or > SUBROUTINE MySolver(x,y) > > In the first case, x is over written. In the second, y = N(x). In any event, > I want to do something like what is shown in the subroutine at the bottom of > this email. > > The code below “works” in the sense that MySolver is called, but it is called > exactly *once*. But MyMonitor and MyConverged are *not* called. Again, I want > to iterate so MySolver should be called many times, as should MyMonitor and > MyConverged. > > The SNESolve() method is called once per nonlinear solve, just as KSPSolve() > is called once per linear solve. There may be iteration inside the method, > but that is handled inside the particular implementation. For example, both > Newton's method and Nonlinear Conjugate Gradient iterate, but the iteration > is internal to both, and they both call the monitor and convergence test at > each internal iterate. > > So, if your nonlinear solver should iterate, it should happen inside the > SNESSolve call for the SNESSHELL object. Does this make sense? What does your > solver do? > > Thanks, > > Matt > > The SNESView before and after SNESSolve looks like this: > > SNES Object: 1 MPI process > type: shell > SNES has not been set up so information may be incomplete > maximum iterations=50, maximum function evaluations=10000 > tolerances: relative=1e-50, absolute=1e-10, solution=1e+06 > total number of function evaluations=0 > norm schedule ALWAYS > > SNES Object: 1 MPI process > type: shell > maximum iterations=50, maximum function evaluations=10000 > tolerances: relative=1e-50, absolute=1e-10, solution=1e+06 > total number of function evaluations=0 > norm schedule ALWAYS > > Any suggestions on how to do what I am trying to accomplish? > > Thanks. > Kenneth Hall > > > > #include <petsc/finclude/petsc.h> > #include "macros.h" > > MODULE SolveWithSNESShell_module > USE MyPetscModule > CONTAINS > ! > !==================================================================================================== > SUBROUTINE MySolver(snes, x, ierr) > !==================================================================================================== > !! > !! > !==================================================================================================== > ! > USE MyPetscModule > IMPLICIT NONE > ! > !.... declared passed variables > SNES :: snes > Vec :: x > PetscErrorCode :: ierr > ! > !.... code to find residual x := N(x) > !.... (or alternatively y := N(x)) > > END SUBROUTINE MySolver > ! > !==================================================================================================== > SUBROUTINE MyMonitor(snes, its, rnorm, ierr) > !==================================================================================================== > !! > !! > !==================================================================================================== > ! > USE MyPetscModule > IMPLICIT NONE > ! > !.... Declare passed variables > SNES :: snes > PetscInt :: its > PetscReal :: rnorm > PetscErrorCode :: ierr > ! > !.... Code to print out convergence history > !.... Code to print out convergence history > > END SUBROUTINE MyMonitor > > !==================================================================================================== > SUBROUTINE MyConverged(snes, it, xnorm, ynorm, znorm, reason, ierr) > USE MyPetscModule > IMPLICIT NONE > > SNES :: snes > PetscInt :: it,ctx > PetscReal :: xnorm, ynorm, znorm > KSPConvergedReason :: reason > PetscErrorCode :: ierr > > ! ... add convergence test here ... > ! set reason to a positive value if convergence has been achieved > > > END SUBROUTINE MyConverged > END MODULE SolveWithSNESShell_module > ! > !==================================================================================================== > SUBROUTINE SolveWithSNESShell > !==================================================================================================== > !! > !! > !==================================================================================================== > ! > USE SolveWithSNESShell_module > IMPLICIT NONE > ! > !.... Declare passed variables > INTEGER :: level_tmp > ! > !.... Declare local variables > INTEGER :: iz > INTEGER :: imax > INTEGER :: jmax > INTEGER :: kmax > SNES :: snes > KSP :: ksp > Vec :: x > Vec :: y > PetscViewer :: viewer > PetscErrorCode :: ierr > PetscReal :: rtol = 1.0D-10 !! relative tolerance > PetscReal :: atol = 1.0D-50 !! absolute tolerance > PetscReal :: dtol = 1.0D+06 !! divergence tolerance > PetscInt :: maxits = 50 > PetscInt :: maxf = 10000 > character(len=1000):: args > ! > !.... count the number of degrees of freedom. > level = level_tmp > n = 0 > DO iz = 1, hb(level)%nzone > imax = hb(level)%zone(iz)%imax - 1 > jmax = hb(level)%zone(iz)%jmax - 1 > kmax = hb(level)%zone(iz)%kmax - 1 > n = n + imax * jmax * kmax > END DO > n = n * neqn > ! > !.... Initialize PETSc > PetscCall(PetscInitialize(PETSC_NULL_CHARACTER, ierr)) > ! > !.... Log > PetscCall(PetscLogDefaultBegin(ierr)) > ! > !.... Hard-wired options. > ! PetscCall(PetscOptionsInsertString(PETSC_NULL_OPTIONS, "command line > style option here" , ierr)) > ! > !.... Command line options. > call GET_COMMAND(args) > PetscCall(PetscOptionsInsertString(PETSC_NULL_OPTIONS, args, ierr)) > ! > !.... view command line table > PetscCall(PetscViewerASCIIOpen(PETSC_COMM_SELF, > PETSC_VIEWER_STDOUT_SELF, viewer, ierr)) > PetscCall(PetscOptionsView(PETSC_NULL_OPTIONS, viewer, ierr)) > PetscCall(PetscViewerDestroy(viewer, ierr)) > ! > !.... Create PETSc vectors > PetscCall(VecCreateSeq(PETSC_COMM_SELF, n, x, ierr)) > PetscCall(VecCreateSeq(PETSC_COMM_SELF, n, y, ierr)) > PetscCall(VecSet(x, 0.0d0, ierr)) > PetscCall(VecSet(y, 0.0d0, ierr)) > > !.... SNES context > PetscCall(SNESCreate(PETSC_COMM_SELF, snes, ierr)) > PetscCall(SNESSetType(snes, SNESSHELL, ierr)) > PetscCall(SNESShellSetSolve(snes, MySolver, ierr)) > > !!!! PetscCall(SNESSetFunction(snes, x, MySolver, PETSC_NULL_INTEGER, ierr)) > !!!! this line causes a segmentation error if uncommented. > > PetscCall(SNESSetConvergenceTest(snes, MyConverged, 0, > PETSC_NULL_FUNCTION, ierr)) > ! > !.... Set SNES options > PetscCall(SNESSetFromOptions(snes, ierr)) > > !.... Set tolerances > PetscCall(SNESSetTolerances(snes, rtol, atol, dtol, maxits, maxf, ierr)) > ! > !.... SNES montior > PetscCall(SNESMonitorSet(snes, MyMonitor, PETSC_NULL_INTEGER, > PETSC_NULL_FUNCTION,ierr)) > ! > !.... Set the initial solution > CALL HBToVecX(x) > ! > !.... View snes context > PetscCall(SNESView(snes, viewer, ierr)) > ! > !.... Solve SNES problem > PetscCall(SNESSolve(snes, PETSC_NULL_VEC, x, ierr)) > ! > !.... View snes context > PetscCall(SNESView(snes, viewer, ierr)) > ! > !.... dump the logs > ! call PetscLogDump(ierr) ! Why does this cause error > ! > !.... Destroy PETSc objects > PetscCall(SNESDestroy(snes, ierr)) > PetscCall(VecDestroy(x, ierr)) > PetscCall(VecDestroy(y, ierr)) > PetscCall(PetscViewerDestroy(viewer, ierr)) > ! > !.... Finish > PetscCall(PetscFinalize(ierr)) > > END SUBROUTINE SolveWithSNESShell > > > > -- > 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/ > <https://urldefense.com/v3/__http:/www.cse.buffalo.edu/*knepley/__;fg!!OToaGQ!sE2W1qI_dqcEHL1dOCnJ3Rdv9TATFVDBiBqx_tlQsOnjvGF7StDjsmVcm9Qkfe4XcTFOBtVjtFl5om07Rdjw$>
