Actual source code: snesimpl.h

  1: #pragma once

  3: #include <petscsnes.h>
  4: #include <petsc/private/petscimpl.h>

  6: /* SUBMANSEC = SNES */

  8: PETSC_EXTERN PetscBool      SNESRegisterAllCalled;
  9: PETSC_EXTERN PetscErrorCode SNESRegisterAll(void);

 11: typedef struct _SNESOps *SNESOps;

 13: struct _SNESOps {
 14:   PetscErrorCode (*computeinitialguess)(SNES, Vec, void *);
 15:   PetscErrorCode (*computescaling)(Vec, Vec, void *);
 16:   PetscErrorCode (*update)(SNES, PetscInt); /* General purpose function for update */
 17:   PetscErrorCode (*converged)(SNES, PetscInt, PetscReal, PetscReal, PetscReal, SNESConvergedReason *, void *);
 18:   PetscErrorCode (*convergeddestroy)(void *);
 19:   PetscErrorCode (*setup)(SNES); /* routine to set up the nonlinear solver */
 20:   PetscErrorCode (*solve)(SNES); /* actual nonlinear solver */
 21:   PetscErrorCode (*view)(SNES, PetscViewer);
 22:   PetscErrorCode (*setfromoptions)(SNES, PetscOptionItems *); /* sets options from database */
 23:   PetscErrorCode (*destroy)(SNES);
 24:   PetscErrorCode (*reset)(SNES);
 25:   PetscErrorCode (*usercompute)(SNES, void **);
 26:   PetscCtxDestroyFn *ctxdestroy;
 27:   PetscErrorCode (*computevariablebounds)(SNES, Vec, Vec); /* user provided routine to set box constrained variable bounds */
 28:   PetscErrorCode (*computepfunction)(SNES, Vec, Vec, void *);
 29:   PetscErrorCode (*computepjacobian)(SNES, Vec, Mat, Mat, void *);
 30:   PetscErrorCode (*load)(SNES, PetscViewer);
 31: };

 33: /*
 34:    Nonlinear solver context
 35:  */
 36: #define MAXSNESMONITORS    5
 37: #define MAXSNESREASONVIEWS 5

 39: struct _p_SNES {
 40:   PETSCHEADER(struct _SNESOps);
 41:   DM        dm;
 42:   PetscBool dmAuto; /* SNES created currently used DM automatically */
 43:   SNES      npc;
 44:   PCSide    npcside;
 45:   PetscBool usesnpc; /* type can use a nonlinear preconditioner */

 47:   /*  ------------------------ User-provided stuff -------------------------------*/
 48:   void *ctx; /* user-defined context */

 50:   Vec vec_rhs; /* If non-null, solve F(x) = rhs */
 51:   Vec vec_sol; /* pointer to solution */

 53:   Vec vec_func; /* pointer to function */

 55:   Mat            jacobian;      /* Jacobian matrix */
 56:   Mat            jacobian_pre;  /* preconditioner matrix */
 57:   Mat            picard;        /* copy of preconditioner matrix needed for Picard with -snes_mf_operator */
 58:   void          *initialguessP; /* user-defined initial guess context */
 59:   KSP            ksp;           /* linear solver context */
 60:   SNESLineSearch linesearch;    /* line search context */
 61:   PetscBool      usesksp;
 62:   MatStructure   matstruct; /* Used by Picard solver */

 64:   Vec vec_sol_update; /* pointer to solution update */

 66:   Vec   scaling; /* scaling vector */
 67:   void *scaP;    /* scaling context */

 69:   PetscReal precheck_picard_angle; /* For use with SNESLineSearchPreCheckPicard */

 71:   /* ------------------------Time stepping hooks-----------------------------------*/

 73:   /* ---------------- PETSc-provided (or user-provided) stuff ---------------------*/

 75:   PetscErrorCode (*monitor[MAXSNESMONITORS])(SNES, PetscInt, PetscReal, void *); /* monitor routine */
 76:   PetscCtxDestroyFn  *monitordestroy[MAXSNESMONITORS];                           /* monitor context destroy routine */
 77:   void               *monitorcontext[MAXSNESMONITORS];                           /* monitor context */
 78:   PetscInt            numbermonitors;                                            /* number of monitors */
 79:   PetscBool           pauseFinal;                                                /* pause all drawing monitor at the final iterate */
 80:   void               *cnvP;                                                      /* convergence context */
 81:   SNESConvergedReason reason;                                                    /* converged reason */

 83:   PetscViewer       convergedreasonviewer;
 84:   PetscViewerFormat convergedreasonformat;
 85:   PetscErrorCode (*reasonview[MAXSNESREASONVIEWS])(SNES, void *); /* snes converged reason view */
 86:   PetscCtxDestroyFn *reasonviewdestroy[MAXSNESREASONVIEWS];       /* reason view context destroy routine */
 87:   void              *reasonviewcontext[MAXSNESREASONVIEWS];       /* reason view context */
 88:   PetscInt           numberreasonviews;                           /* number of reason views */
 89:   PetscBool          errorifnotconverged;

 91:   /* --- Routines and data that are unique to each particular solver --- */

 93:   PetscBool setupcalled; /* true if setup has been called */
 94:   void     *data;        /* implementation-specific data */

 96:   /* --------------------------  Parameters -------------------------------------- */

 98:   PetscInt  nfuncs;                                 /* number of function evaluations */
 99:   PetscInt  iter;                                   /* global iteration number */
100:   PetscInt  linear_its;                             /* total number of linear solver iterations */
101:   PetscReal norm;                                   /* residual norm of current iterate */
102:   PetscReal ynorm;                                  /* update norm of current iterate */
103:   PetscReal xnorm;                                  /* solution norm of current iterate */
104:   PetscBool forceiteration;                         /* Force SNES to take at least one iteration regardless of the initial residual norm */
105:   PetscInt  lagpreconditioner;                      /* SNESSetLagPreconditioner() */
106:   PetscInt  lagjacobian;                            /* SNESSetLagJacobian() */
107:   PetscInt  jac_iter;                               /* The present iteration of the Jacobian lagging */
108:   PetscBool lagjac_persist;                         /* The jac_iter persists until reset */
109:   PetscInt  pre_iter;                               /* The present iteration of the Preconditioner lagging */
110:   PetscBool lagpre_persist;                         /* The pre_iter persists until reset */
111:   PetscInt  gridsequence;                           /* number of grid sequence steps to take; defaults to zero */
112:   PetscObjectParameterDeclare(PetscInt, max_its);   /* max number of iterations */
113:   PetscObjectParameterDeclare(PetscInt, max_funcs); /* max number of function evals */
114:   PetscObjectParameterDeclare(PetscReal, rtol);     /* relative tolerance */
115:   PetscObjectParameterDeclare(PetscReal, divtol);   /* relative divergence tolerance */
116:   PetscObjectParameterDeclare(PetscReal, abstol);   /* absolute tolerance */
117:   PetscObjectParameterDeclare(PetscReal, stol);     /* step length tolerance*/

119:   PetscBool vec_func_init_set; /* the initial function has been set */

121:   SNESNormSchedule normschedule; /* Norm computation type for SNES instance */
122:   SNESFunctionType functype;     /* Function type for the SNES instance */

124:   /* ------------------------ Default work-area management ---------------------- */

126:   PetscInt nwork;
127:   Vec     *work;

129:   /* ------------------------- Miscellaneous Information ------------------------ */

131:   PetscInt   setfromoptionscalled;
132:   PetscReal *conv_hist;       /* If !0, stores function norm (or
133:                                     gradient norm) at each iteration */
134:   PetscInt  *conv_hist_its;   /* linear iterations for each Newton step */
135:   size_t     conv_hist_len;   /* size of convergence history array */
136:   size_t     conv_hist_max;   /* actual amount of data in conv_history */
137:   PetscBool  conv_hist_reset; /* reset counter for each new SNES solve */
138:   PetscBool  conv_hist_alloc;
139:   PetscBool  counters_reset; /* reset counter for each new SNES solve */

141:   /* the next two are used for failures in the line search; they should be put elsewhere */
142:   PetscInt numFailures; /* number of unsuccessful step attempts */
143:   PetscInt maxFailures; /* maximum number of unsuccessful step attempts */

145:   PetscInt numLinearSolveFailures;
146:   PetscInt maxLinearSolveFailures;

148:   PetscBool domainerror;         /* set with SNESSetFunctionDomainError() */
149:   PetscBool jacobiandomainerror; /* set with SNESSetJacobianDomainError() */
150:   PetscBool checkjacdomainerror; /* if or not check Jacobian domain error after Jacobian evaluations */

152:   PetscBool ksp_ewconv; /* flag indicating use of Eisenstat-Walker KSP convergence criteria */
153:   void     *kspconvctx; /* Eisenstat-Walker KSP convergence context */

155:   /* SNESConvergedDefault context: split it off into a separate var/struct to be passed as context to SNESConvergedDefault? */
156:   PetscReal ttol;   /* rtol*initial_residual_norm */
157:   PetscReal rnorm0; /* initial residual norm (used for divergence testing) */

159:   Vec     *vwork; /* more work vectors for Jacobian approx */
160:   PetscInt nvwork;

162:   PetscBool mf;          /* -snes_mf OR -snes_mf_operator was used on this snes */
163:   PetscBool mf_operator; /* -snes_mf_operator was used on this snes */
164:   PetscInt  mf_version;  /* The version of snes_mf used */

166:   PetscReal vizerotolerance; /* tolerance for considering an x[] value to be on the bound */
167:   Vec       xl, xu;          /* upper and lower bounds for box constrained VI problems */
168:   PetscInt  ntruebounds;     /* number of non-infinite bounds set for VI box constraints */
169:   PetscBool usersetbounds;   /* bounds have been set via SNESVISetVariableBounds(), rather than via computevariablebounds() callback. */

171:   PetscBool alwayscomputesfinalresidual; /* Does SNESSolve_XXX always compute the value of the residual at the final
172:                                              * solution and put it in vec_func?  Used inside SNESSolve_FAS to determine
173:                                              * if the final residual must be computed before restricting or prolonging
174:                                              * it. */
175: };

177: typedef struct _p_DMSNES  *DMSNES;
178: typedef struct _DMSNESOps *DMSNESOps;
179: struct _DMSNESOps {
180:   SNESFunctionFn *computefunction;
181:   SNESFunctionFn *computemffunction;
182:   SNESJacobianFn *computejacobian;

184:   /* objective */
185:   SNESObjectiveFn *computeobjective;

187:   /* Picard iteration functions */
188:   SNESFunctionFn *computepfunction;
189:   SNESJacobianFn *computepjacobian;

191:   /* User-defined smoother */
192:   SNESNGSFn *computegs;

194:   PetscErrorCode (*destroy)(DMSNES);
195:   PetscErrorCode (*duplicate)(DMSNES, DMSNES);
196: };

198: /*S
199:    DMSNES - Object held by a `DM` that contains all the callback functions and their contexts needed by a `SNES`

201:    Level: developer

203:    Notes:
204:    Users provides callback functions and their contexts to `SNES` using, for example, `SNESSetFunction()`. These values are stored
205:    in a `DMSNES` that is contained in the `DM` associated with the `SNES`. If no `DM` was provided by
206:    the user with `SNESSetDM()` it is automatically created by `SNESGetDM()` with `DMShellCreate()`.

208:    Users very rarely need to worked directly with the `DMSNES` object, rather they work with the `SNES` and the `DM` they created

210:    Multiple `DM` can share a single `DMSNES`, often each `DM` is associated with
211:    a grid refinement level. `DMGetDMSNES()` returns the `DMSNES` associated with a `DM`. `DMGetDMSNESWrite()` returns a unique
212:    `DMSNES` that is only associated with the current `DM`, making a copy of the shared `DMSNES` if needed (copy-on-write).

214:    See `DMKSP` for details on why there is a needed for `DMSNES` instead of simply storing the user callbacks directly in the `DM` or the `TS`

216:    Developer Note:
217:    The `originaldm` inside the `DMSNES` is NOT reference counted (to prevent a reference count loop between a `DM` and a `DMSNES`).
218:    The `DM` on which this context was first created is cached here to implement one-way
219:    copy-on-write. When `DMGetDMSNESWrite()` sees a request using a different `DM`, it makes a copy of the `TSDM`. Thus, if a user
220:    only interacts directly with one level, e.g., using `TSSetIFunction()`, then coarse levels of a multilevel item
221:    integrator are built, then the user changes the routine with another call to `TSSetIFunction()`, it automatically
222:    propagates to all the levels. If instead, they get out a specific level and set the function on that level,
223:    subsequent changes to the original level will no longer propagate to that level.

225: .seealso: [](ch_ts), `SNES`, `SNESCreate()`, `DM`, `DMGetDMSNESWrite()`, `DMGetDMSNES()`,  `DMKSP`, `DMTS`, `DMSNESSetFunction()`, `DMSNESGetFunction()`,
226:           `DMSNESSetFunctionContextDestroy()`, `DMSNESSetMFFunction()`, `DMSNESSetNGS()`, `DMSNESGetNGS()`, `DMSNESSetJacobian()`, `DMSNESGetJacobian()`,
227:           `DMSNESSetJacobianContextDestroy()`, `DMSNESSetPicard()`, `DMSNESGetPicard()`, `DMSNESSetObjective()`, `DMSNESGetObjective()`, `DMCopyDMSNES()`
228: S*/
229: struct _p_DMSNES {
230:   PETSCHEADER(struct _DMSNESOps);
231:   PetscContainer functionctxcontainer;
232:   PetscContainer jacobianctxcontainer;
233:   void          *mffunctionctx;
234:   void          *gsctx;
235:   void          *pctx;
236:   void          *objectivectx;

238:   void *data;

240:   /* See developer note for DMSNES above */
241:   DM originaldm;
242: };
243: PETSC_EXTERN PetscErrorCode DMGetDMSNES(DM, DMSNES *);
244: PETSC_EXTERN PetscErrorCode DMSNESView(DMSNES, PetscViewer);
245: PETSC_EXTERN PetscErrorCode DMSNESLoad(DMSNES, PetscViewer);
246: PETSC_EXTERN PetscErrorCode DMGetDMSNESWrite(DM, DMSNES *);

248: /* Context for Eisenstat-Walker convergence criteria for KSP solvers */
249: typedef struct {
250:   PetscInt  version;     /* flag indicating version (1,2,3 or 4) */
251:   PetscReal rtol_0;      /* initial rtol */
252:   PetscReal rtol_last;   /* last rtol */
253:   PetscReal rtol_max;    /* maximum rtol */
254:   PetscReal gamma;       /* mult. factor for version 2 rtol computation */
255:   PetscReal alpha;       /* power for version 2 rtol computation */
256:   PetscReal alpha2;      /* power for safeguard */
257:   PetscReal threshold;   /* threshold for imposing safeguard */
258:   PetscReal lresid_last; /* linear residual from last iteration */
259:   PetscReal norm_last;   /* function norm from last iteration */
260:   PetscReal norm_first;  /* function norm from the beginning of the first iteration. */
261:   PetscReal rtol_last_2, rk_last, rk_last_2;
262:   PetscReal v4_p1, v4_p2, v4_p3, v4_m1, v4_m2, v4_m3, v4_m4;
263: } SNESKSPEW;

265: static inline PetscErrorCode SNESLogConvergenceHistory(SNES snes, PetscReal res, PetscInt its)
266: {
267:   PetscFunctionBegin;
268:   PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
269:   if (snes->conv_hist && snes->conv_hist_max > snes->conv_hist_len) {
270:     if (snes->conv_hist) snes->conv_hist[snes->conv_hist_len] = res;
271:     if (snes->conv_hist_its) snes->conv_hist_its[snes->conv_hist_len] = its;
272:     snes->conv_hist_len++;
273:   }
274:   PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
275:   PetscFunctionReturn(PETSC_SUCCESS);
276: }

278: PETSC_EXTERN PetscErrorCode SNESVIProjectOntoBounds(SNES, Vec);
279: PETSC_INTERN PetscErrorCode SNESVICheckLocalMin_Private(SNES, Mat, Vec, Vec, PetscReal, PetscBool *);
280: PETSC_INTERN PetscErrorCode SNESReset_VI(SNES);
281: PETSC_INTERN PetscErrorCode SNESDestroy_VI(SNES);
282: PETSC_INTERN PetscErrorCode SNESView_VI(SNES, PetscViewer);
283: PETSC_INTERN PetscErrorCode SNESSetFromOptions_VI(SNES, PetscOptionItems *);
284: PETSC_INTERN PetscErrorCode SNESSetUp_VI(SNES);
285: PETSC_EXTERN_TYPEDEF typedef PetscErrorCode(SNESVIComputeVariableBoundsFn)(SNES, Vec, Vec);
286: PETSC_EXTERN_TYPEDEF typedef SNESVIComputeVariableBoundsFn *SNESVIComputeVariableBoundsFunction; // deprecated version
287: PETSC_INTERN PetscErrorCode                                 SNESVISetComputeVariableBounds_VI(SNES, SNESVIComputeVariableBoundsFn);
288: PETSC_INTERN PetscErrorCode                                 SNESVISetVariableBounds_VI(SNES, Vec, Vec);
289: PETSC_INTERN PetscErrorCode                                 SNESConvergedDefault_VI(SNES, PetscInt, PetscReal, PetscReal, PetscReal, SNESConvergedReason *, void *);

291: PETSC_EXTERN PetscErrorCode DMSNESUnsetFunctionContext_Internal(DM);
292: PETSC_EXTERN PetscErrorCode DMSNESUnsetJacobianContext_Internal(DM);
293: PETSC_EXTERN PetscErrorCode DMSNESCheck_Internal(SNES, DM, Vec);

295: PETSC_EXTERN PetscLogEvent SNES_Solve;
296: PETSC_EXTERN PetscLogEvent SNES_SetUp;
297: PETSC_EXTERN PetscLogEvent SNES_FunctionEval;
298: PETSC_EXTERN PetscLogEvent SNES_JacobianEval;
299: PETSC_EXTERN PetscLogEvent SNES_NGSEval;
300: PETSC_EXTERN PetscLogEvent SNES_NGSFuncEval;
301: PETSC_EXTERN PetscLogEvent SNES_NewtonALEval;
302: PETSC_EXTERN PetscLogEvent SNES_NPCSolve;
303: PETSC_EXTERN PetscLogEvent SNES_ObjectiveEval;

305: PETSC_INTERN PetscBool  SNEScite;
306: PETSC_INTERN const char SNESCitation[];

308: /* Used by TAOBNK solvers */
309: PETSC_EXTERN PetscErrorCode KSPPostSolve_SNESEW(KSP, Vec, Vec, void *);
310: PETSC_EXTERN PetscErrorCode KSPPreSolve_SNESEW(KSP, Vec, Vec, void *);
311: PETSC_EXTERN PetscErrorCode SNESEWSetFromOptions_Private(SNESKSPEW *, PetscBool, MPI_Comm, const char *);

313: /*
314:     Either generate an error or mark as diverged when a real from a SNES function norm is Nan or Inf.
315:     domainerror is reset here, once reason is set, to allow subsequent iterations to be feasible (e.g. line search).
316: */
317: #define SNESCheckFunctionNorm(snes, beta) \
318:   do { \
319:     if (PetscIsInfOrNanReal(beta)) { \
320:       PetscCheck(!snes->errorifnotconverged, PetscObjectComm((PetscObject)snes), PETSC_ERR_NOT_CONVERGED, "SNESSolve has not converged due to Nan or Inf norm"); \
321:       { \
322:         PetscBool domainerror; \
323:         PetscCallMPI(MPIU_Allreduce(&snes->domainerror, &domainerror, 1, MPIU_BOOL, MPI_LOR, PetscObjectComm((PetscObject)snes))); \
324:         if (domainerror) { \
325:           snes->reason      = SNES_DIVERGED_FUNCTION_DOMAIN; \
326:           snes->domainerror = PETSC_FALSE; \
327:         } else snes->reason = SNES_DIVERGED_FNORM_NAN; \
328:         PetscFunctionReturn(PETSC_SUCCESS); \
329:       } \
330:     } \
331:   } while (0)

333: #define SNESCheckJacobianDomainerror(snes) \
334:   do { \
335:     if (snes->checkjacdomainerror) { \
336:       PetscBool domainerror; \
337:       PetscCallMPI(MPIU_Allreduce(&snes->jacobiandomainerror, &domainerror, 1, MPIU_BOOL, MPI_LOR, PetscObjectComm((PetscObject)snes))); \
338:       if (domainerror) { \
339:         snes->reason              = SNES_DIVERGED_JACOBIAN_DOMAIN; \
340:         snes->jacobiandomainerror = PETSC_FALSE; \
341:         PetscCheck(!snes->errorifnotconverged, PetscObjectComm((PetscObject)snes), PETSC_ERR_NOT_CONVERGED, "SNESSolve has not converged due to Jacobian domain error"); \
342:         PetscFunctionReturn(PETSC_SUCCESS); \
343:       } \
344:     } \
345:   } while (0)

347: #define SNESCheckKSPSolve(snes) \
348:   do { \
349:     KSPConvergedReason kspreason; \
350:     PetscInt           lits; \
351:     PetscCall(KSPGetIterationNumber(snes->ksp, &lits)); \
352:     snes->linear_its += lits; \
353:     PetscCall(KSPGetConvergedReason(snes->ksp, &kspreason)); \
354:     if (kspreason < 0) { \
355:       if (kspreason == KSP_DIVERGED_NANORINF) { \
356:         PetscBool domainerror; \
357:         PetscCallMPI(MPIU_Allreduce(&snes->domainerror, &domainerror, 1, MPIU_BOOL, MPI_LOR, PetscObjectComm((PetscObject)snes))); \
358:         if (domainerror) snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; \
359:         else snes->reason = SNES_DIVERGED_LINEAR_SOLVE; \
360:         PetscFunctionReturn(PETSC_SUCCESS); \
361:       } else { \
362:         if (++snes->numLinearSolveFailures >= snes->maxLinearSolveFailures) { \
363:           PetscCall(PetscInfo(snes, "iter=%" PetscInt_FMT ", number linear solve failures %" PetscInt_FMT " greater than current SNES allowed %" PetscInt_FMT ", stopping solve\n", snes->iter, snes->numLinearSolveFailures, snes->maxLinearSolveFailures)); \
364:           snes->reason = SNES_DIVERGED_LINEAR_SOLVE; \
365:           PetscFunctionReturn(PETSC_SUCCESS); \
366:         } \
367:       } \
368:     } \
369:   } while (0)

371: #define SNESNeedNorm_Private(snes, iter) \
372:   (((iter) == (snes)->max_its && ((snes)->normschedule == SNES_NORM_FINAL_ONLY || (snes)->normschedule == SNES_NORM_INITIAL_FINAL_ONLY)) || ((iter) == 0 && ((snes)->normschedule == SNES_NORM_INITIAL_ONLY || (snes)->normschedule == SNES_NORM_INITIAL_FINAL_ONLY)) || \
373:    (snes)->normschedule == SNES_NORM_ALWAYS)