Using SNESSHELL as a wrapper for a CFD solver.
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 !.... 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
On Wed, May 31, 2023 at 3:21 PM Kenneth C Hall <[email protected]> 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/>
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 <[email protected]> Date: Wednesday, May 31, 2023 at 3:48 PM To: Kenneth C Hall <[email protected]> Cc: [email protected] <[email protected]> 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 <[email protected]<mailto:[email protected]>> 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$>
On May 31, 2023, at 4:16 PM, Kenneth C Hall <[email protected]> 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 <[email protected] <mailto:[email protected]>> Date: Wednesday, May 31, 2023 at 3:48 PM To: Kenneth C Hall <[email protected] <mailto:[email protected]>> Cc: [email protected] <mailto:[email protected]> <[email protected] <mailto:[email protected]>> 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 <[email protected] <mailto:[email protected]>> 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$>
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.
On May 31, 2023, at 4:16 PM, Kenneth C Hall <[email protected]> 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 <[email protected] <mailto:[email protected]>> Date: Wednesday, May 31, 2023 at 3:48 PM To: Kenneth C Hall <[email protected] <mailto:[email protected]>> Cc: [email protected] <mailto:[email protected]> <[email protected] <mailto:[email protected]>> 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 <[email protected] <mailto:[email protected]>> 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$>
On Wed, May 31, 2023 at 5:51 PM Barry Smith <[email protected]> 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 <[email protected]> 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 <[email protected]> *Date: *Wednesday, May 31, 2023 at 3:48 PM *To: *Kenneth C Hall <[email protected]> *Cc: *[email protected] <[email protected]> *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 <[email protected]> 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/ <http://www.cse.buffalo.edu/~knepley/>
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 <[email protected]> Date: Wednesday, May 31, 2023 at 6:21 PM To: Barry Smith <[email protected]> Cc: Kenneth C Hall <[email protected]>, [email protected] <[email protected]> Subject: Re: [petsc-users] Using SNESSHELL as a wrapper for a CFD solver. On Wed, May 31, 2023 at 5:51 PM Barry Smith <[email protected]<mailto:[email protected]>> 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 <[email protected]<mailto:[email protected]>> 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 <[email protected]<mailto:[email protected]>> Date: Wednesday, May 31, 2023 at 3:48 PM To: Kenneth C Hall <[email protected]<mailto:[email protected]>> Cc: [email protected]<mailto:[email protected]> <[email protected]<mailto:[email protected]>> 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 <[email protected]<mailto:[email protected]>> 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$>
Barry, Thanks for your reply. Yes, for SNESNGMRES, I use N(x) – x as the residuals, with basic line search. This converges in half the “iterations” as the original CFD solver, but each SNESNGMRES iteration requires two residual calls, so ends up converging in almost exactly the same amount of time as the CFD solver alone. I will try the SNESSetNPC/SNESNRICHARSON/SNESSHELL combination for the “bare CFD” version. Thanks for your helpful comments and reply. Kenneth From: Barry Smith <[email protected]> Date: Wednesday, May 31, 2023 at 5:50 PM To: Kenneth C Hall <[email protected]> Cc: [email protected] <[email protected]> Subject: Re: [petsc-users] Using SNESSHELL as a wrapper for a CFD solver. 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. On May 31, 2023, at 4:16 PM, Kenneth C Hall <[email protected]> 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 <[email protected]<mailto:[email protected]>> Date: Wednesday, May 31, 2023 at 3:48 PM To: Kenneth C Hall <[email protected]<mailto:[email protected]>> Cc: [email protected]<mailto:[email protected]> <[email protected]<mailto:[email protected]>> 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 <[email protected]<mailto:[email protected]>> 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$>
You should also try the SNESNGRES with your SNESSHELL as the SNESSetPC().
On May 31, 2023, at 6:33 PM, Kenneth C Hall <[email protected]> wrote:
Barry,
Thanks for your reply. Yes, for SNESNGMRES, I use N(x) – x as the residuals, with basic line search. This converges in half the “iterations” as the original CFD solver, but each SNESNGMRES iteration requires two residual calls, so ends up converging in almost exactly the same amount of time as the CFD solver alone.
I will try the SNESSetNPC/SNESNRICHARSON/SNESSHELL combination for the “bare CFD” version.
Thanks for your helpful comments and reply.
Kenneth
From: Barry Smith <[email protected] <mailto:[email protected]>> Date: Wednesday, May 31, 2023 at 5:50 PM To: Kenneth C Hall <[email protected] <mailto:[email protected]>> Cc: [email protected] <mailto:[email protected]> <[email protected] <mailto:[email protected]>> Subject: Re: [petsc-users] Using SNESSHELL as a wrapper for a CFD solver.
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.
On May 31, 2023, at 4:16 PM, Kenneth C Hall <[email protected] <mailto:[email protected]>> 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 <[email protected] <mailto:[email protected]>> Date: Wednesday, May 31, 2023 at 3:48 PM To: Kenneth C Hall <[email protected] <mailto:[email protected]>> Cc: [email protected] <mailto:[email protected]> <[email protected] <mailto:[email protected]>> 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 <[email protected] <mailto:[email protected]>> 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$>
participants (3)
-
Barry Smith -
Kenneth C Hall -
Matthew Knepley