[petsc-users] Using SNESSHELL as a wrapper for a CFD solver.
Matthew Knepley
knepley at gmail.com
Wed May 31 17:21:14 CDT 2023
On Wed, May 31, 2023 at 5:51 PM Barry Smith <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>
> 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>
> *Date: *Wednesday, May 31, 2023 at 3:48 PM
> *To: *Kenneth C Hall <kenneth.c.hall at duke.edu>
> *Cc: *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 3:21 PM Kenneth C Hall <kenneth.c.hall at duke.edu>
> wrote:
>
> Hi,
>
> I am doing a number of problems using PETSc/SLEPc, but I also work on some
> non-PETSc/SLEPc flow solvers. I would like to use PETSc as a wrapper for
> this non-PETSc flow solver for compatibility, so I can use the tolerance
> monitoring, options, viewers, and for direct comparison to PETSc methods I
> am using.
>
> Here is what I am trying to do… I have a CFD solver that iterates with a
> nonlinear iterator of the form x := N(x). This can be expressed in a
> fortran routine of the form,
>
> SUBROUTINE MySolver(x)
> or
> SUBROUTINE MySolver(x,y)
>
> In the first case, x is over written. In the second, y = N(x). In any
> event, I want to do something like what is shown in the subroutine at the
> bottom of this email.
>
> The code below “works” in the sense that MySolver is called, but it is
> called exactly **once**. But MyMonitor and MyConverged are **not**
> called. Again, I want to iterate so MySolver should be called many times,
> as should MyMonitor and MyConverged.
>
>
> The SNESolve() method is called once per nonlinear solve, just as
> KSPSolve() is called once per linear solve. There may be iteration inside
> the method, but that is handled inside the particular implementation. For
> example, both Newton's method and Nonlinear Conjugate Gradient iterate, but
> the iteration is internal to both, and they both call the monitor and
> convergence test at each internal iterate.
>
> So, if your nonlinear solver should iterate, it should happen inside the
> SNESSolve call for the SNESSHELL object. Does this make sense? What does
> your solver do?
>
> Thanks,
>
> Matt
>
>
> The SNESView before and after SNESSolve looks like this:
>
>
> SNES Object: 1 MPI process
>
> type: shell
>
> SNES has not been set up so information may be incomplete
>
> maximum iterations=50, maximum function evaluations=10000
>
> tolerances: relative=1e-50, absolute=1e-10, solution=1e+06
>
> total number of function evaluations=0
>
> norm schedule ALWAYS
>
>
> SNES Object: 1 MPI process
>
> type: shell
>
> maximum iterations=50, maximum function evaluations=10000
>
> tolerances: relative=1e-50, absolute=1e-10, solution=1e+06
>
> total number of function evaluations=0
>
> norm schedule ALWAYS
>
> Any suggestions on how to do what I am trying to accomplish?
>
> Thanks.
> Kenneth Hall
>
>
>
> *#include <petsc/finclude/petsc.h>*
> *#include "macros.h"*
>
>
> *MODULE* SolveWithSNESShell_module
> *USE* MyPetscModule
> *CONTAINS*
> !
>
> !====================================================================================================
> *SUBROUTINE* MySolver(snes, x, ierr)
>
> !====================================================================================================
> !!
> !!
>
> !====================================================================================================
> !
> *USE* MyPetscModule
> *IMPLICIT* *NONE*
> !
> !.... declared passed variables
> SNES :: snes
> Vec :: x
> PetscErrorCode :: ierr
> !
> !.... code to find residual x := N(x)
> !.... (or alternatively y := N(x))
>
>
> *END* *SUBROUTINE* MySolver
> !
>
> !====================================================================================================
> *SUBROUTINE* MyMonitor(snes, its, rnorm, ierr)
>
> !====================================================================================================
> !!
> !!
>
> !====================================================================================================
> !
> *USE* MyPetscModule
> *IMPLICIT* *NONE*
> !
> !.... Declare passed variables
> SNES :: snes
> PetscInt :: its
> PetscReal :: rnorm
> PetscErrorCode :: ierr
> !
> !.... Code to print out convergence history
> !.... Code to print out convergence history
>
>
> *END* *SUBROUTINE* MyMonitor
>
>
>
> !====================================================================================================
> *SUBROUTINE* MyConverged(snes, it, xnorm, ynorm, znorm, reason, ierr)
> *USE* MyPetscModule
> *IMPLICIT* *NONE*
>
>
> SNES :: snes
> PetscInt :: it,ctx
> PetscReal :: xnorm, ynorm, znorm
> KSPConvergedReason :: reason
> PetscErrorCode :: ierr
>
>
> ! ... add convergence test here ...
> ! set reason to a positive value if convergence has been achieved
>
>
>
>
> *END* *SUBROUTINE* MyConverged
> *END* *MODULE* SolveWithSNESShell_module
> !
>
> !====================================================================================================
> *SUBROUTINE* SolveWithSNESShell
>
> !====================================================================================================
> !!
> !!
>
> !====================================================================================================
> !
> *USE* SolveWithSNESShell_module
> *IMPLICIT* *NONE*
> !
> !.... Declare passed variables
> *INTEGER* :: level_tmp
> !
> !.... Declare local variables
> *INTEGER* :: iz
> *INTEGER* :: imax
> *INTEGER* :: jmax
> *INTEGER* :: kmax
> SNES :: snes
> KSP :: ksp
> Vec :: x
> Vec :: y
> PetscViewer :: viewer
> PetscErrorCode :: ierr
> PetscReal :: rtol = 1.0D-10 !! relative tolerance
> PetscReal :: atol = 1.0D-50 !! absolute tolerance
> PetscReal :: dtol = 1.0D+06 !! divergence tolerance
> PetscInt :: maxits = 50
> PetscInt :: maxf = 10000
> *character*(*len*=1000):: args
> !
> !.... count the number of degrees of freedom.
> level = level_tmp
> n = 0
> *DO* iz = 1, hb(level)%nzone
> imax = hb(level)%zone(iz)%imax - 1
> jmax = hb(level)%zone(iz)%jmax - 1
> kmax = hb(level)%zone(iz)%kmax - 1
> n = n + imax * jmax * kmax
> *END* *DO*
> n = n * neqn
> !
> !.... Initialize PETSc
> PetscCall(PetscInitialize(PETSC_NULL_CHARACTER, ierr))
> !
> !.... Log
> PetscCall(PetscLogDefaultBegin(ierr))
> !
> !.... Hard-wired options.
> ! PetscCall(PetscOptionsInsertString(PETSC_NULL_OPTIONS, "command line
> style option here" , ierr))
> !
> !.... Command line options.
> *call* *GET_COMMAND*(args)
> PetscCall(PetscOptionsInsertString(PETSC_NULL_OPTIONS, args, ierr))
> !
> !.... view command line table
> PetscCall(PetscViewerASCIIOpen(PETSC_COMM_SELF,
> PETSC_VIEWER_STDOUT_SELF, viewer, ierr))
> PetscCall(PetscOptionsView(PETSC_NULL_OPTIONS, viewer, ierr))
> PetscCall(PetscViewerDestroy(viewer, ierr))
> !
> !.... Create PETSc vectors
> PetscCall(VecCreateSeq(PETSC_COMM_SELF, n, x, ierr))
> PetscCall(VecCreateSeq(PETSC_COMM_SELF, n, y, ierr))
> PetscCall(VecSet(x, 0.0d0, ierr))
> PetscCall(VecSet(y, 0.0d0, ierr))
>
>
> !.... SNES context
> PetscCall(SNESCreate(PETSC_COMM_SELF, snes, ierr))
> PetscCall(SNESSetType(snes, SNESSHELL, ierr))
> PetscCall(SNESShellSetSolve(snes, MySolver, ierr))
>
>
> !!!! PetscCall(SNESSetFunction(snes, x, MySolver, PETSC_NULL_INTEGER,
> ierr))
> !!!! this line causes a segmentation error if uncommented.
>
>
> PetscCall(SNESSetConvergenceTest(snes, MyConverged, 0,
> PETSC_NULL_FUNCTION, ierr))
> !
> !.... Set SNES options
> PetscCall(SNESSetFromOptions(snes, ierr))
>
>
> !.... Set tolerances
> PetscCall(SNESSetTolerances(snes, rtol, atol, dtol, maxits, maxf,
> ierr))
> !
> !.... SNES montior
> PetscCall(SNESMonitorSet(snes, MyMonitor, PETSC_NULL_INTEGER,
> PETSC_NULL_FUNCTION,ierr))
> !
> !.... Set the initial solution
> *CALL* HBToVecX(x)
> !
> !.... View snes context
> PetscCall(SNESView(snes, viewer, ierr))
> !
> !.... Solve SNES problem
> PetscCall(SNESSolve(snes, PETSC_NULL_VEC, x, ierr))
> !
> !.... View snes context
> PetscCall(SNESView(snes, viewer, ierr))
> !
> !.... dump the logs
> ! call PetscLogDump(ierr) ! Why does this cause error
> !
> !.... Destroy PETSc objects
> PetscCall(SNESDestroy(snes, ierr))
> PetscCall(VecDestroy(x, ierr))
> PetscCall(VecDestroy(y, ierr))
> PetscCall(PetscViewerDestroy(viewer, ierr))
> !
> !.... Finish
> PetscCall(PetscFinalize(ierr))
>
>
> *END* *SUBROUTINE* SolveWithSNESShell
>
>
>
>
> --
> What most experimenters take for granted before they begin their
> experiments is infinitely more interesting than any results to which their
> experiments lead.
> -- Norbert Wiener
>
> https://www.cse.buffalo.edu/~knepley/
> <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/>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20230531/8b6434ea/attachment-0001.html>
More information about the petsc-users
mailing list