[petsc-users] Using SNESSHELL as a wrapper for a CFD solver.

31 May 2023

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

   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

[petsc-users] Using SNESSHELL as a wrapper for a CFD solver.

Kenneth C Hall