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

Kenneth C Hall kenneth.c.hall at duke.edu
Wed May 31 17:42:39 CDT 2023


Matt,

Thanks for this.

Mach number preconditioning is as follows.  The Euler or Navier-Stokes equations are written as:

du/dt + M.(dF(u)/dx + dG(u)/dx) = 0

The matrix M is a preconditioning matrix which changes the wave speeds of the convection and the pressure characteristic speeds so that they are close to one another. The equations are no longer time accurate, but converge to a steady states faster. This is especially useful for very low speed flows where the (unpreconditioned) waves travel at very high speeds compared to the convection waves.

Kenneth

From: Matthew Knepley <knepley at gmail.com>
Date: Wednesday, May 31, 2023 at 6:21 PM
To: Barry Smith <bsmith at petsc.dev>
Cc: Kenneth C Hall <kenneth.c.hall at duke.edu>, petsc-users at mcs.anl.gov <petsc-users at mcs.anl.gov>
Subject: Re: [petsc-users] Using SNESSHELL as a wrapper for a CFD solver.
On Wed, May 31, 2023 at 5:51 PM Barry Smith <bsmith at petsc.dev<mailto:bsmith at petsc.dev>> wrote:

  Sorry, I wrote to quickly in my last email. You will need to create a SNESSHELL its solve simply calls your solver (for its one iteration) SNESNRICHARSON handles the rest.

Yes, this is a great suggestion. It should work with your current SNESSHELL.

I would also point out that we have the SNESMS, which is the multistage solver from Jameson. You can
use this as a residual smooth with FAS to do something similar to what you have, although I do not know
what Mach number preconditioning is (Jed probably knows).

  Thanks,

    Matt

On May 31, 2023, at 4:16 PM, Kenneth C Hall <kenneth.c.hall at duke.edu<mailto:kenneth.c.hall at 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.  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 <knepley at gmail.com<mailto:knepley at gmail.com>>
Date: Wednesday, May 31, 2023 at 3:48 PM
To: Kenneth C Hall <kenneth.c.hall at duke.edu<mailto:kenneth.c.hall at duke.edu>>
Cc: petsc-users at mcs.anl.gov<mailto:petsc-users at mcs.anl.gov> <petsc-users at mcs.anl.gov<mailto:petsc-users at 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.hall at duke.edu<mailto: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/<https://urldefense.com/v3/__http:/www.cse.buffalo.edu/*knepley/__;fg!!OToaGQ!sE2W1qI_dqcEHL1dOCnJ3Rdv9TATFVDBiBqx_tlQsOnjvGF7StDjsmVcm9Qkfe4XcTFOBtVjtFl5om07Rdjw$>



--
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!s-GnaelCE7pDC_xWpNI3mM5IDxDYowWhn_fSvgL-9JPgdORySPUEHjp-rhUcrlLnhv8LCx38Hr_TCa18TRGy$>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20230531/3cc04ed1/attachment-0001.html>


More information about the petsc-users mailing list