On Wed, May 31, 2023 at 5:51 PM Barry Smith <bsm...@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.h...@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 <knep...@gmail.com> > *Date: *Wednesday, May 31, 2023 at 3:48 PM > *To: *Kenneth C Hall <kenneth.c.h...@duke.edu> > *Cc: *petsc-users@mcs.anl.gov <petsc-users@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.h...@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/>