[petsc-users] Using SNESSHELL as a wrapper for a CFD solver.
Matthew Knepley
knepley at gmail.com
Wed May 31 14:48:00 CDT 2023
On Wed, May 31, 2023 at 3:21 PM Kenneth C Hall <kenneth.c.hall at 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/ <http://www.cse.buffalo.edu/~knepley/>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20230531/cf9319e8/attachment-0001.html>
More information about the petsc-users
mailing list