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;  /* matrix used to construct the preconditioner of the Jacobian */
 57:   Mat            picard;        /* copy of jacobian_pre 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:   /* ---------------------------------- Testing --------------------------------- */
130:   PetscBool testFunc; // Test the function routine
131:   PetscBool testJac;  // Test the Jacobian routine

133:   /* ------------------------- Miscellaneous Information ------------------------ */

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

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

149:   PetscInt numLinearSolveFailures;
150:   PetscInt maxLinearSolveFailures;

152:   PetscBool domainerror;         /* set with SNESSetFunctionDomainError() */
153:   PetscBool jacobiandomainerror; /* set with SNESSetJacobianDomainError() */
154:   PetscBool checkjacdomainerror; /* if or not check Jacobian domain error after Jacobian evaluations */

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

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

163:   Vec     *vwork; /* more work vectors for Jacobian approx */
164:   PetscInt nvwork;

166:   PetscBool mf;          /* -snes_mf OR -snes_mf_operator was used on this snes */
167:   PetscBool mf_operator; /* -snes_mf_operator was used on this snes */
168:   PetscInt  mf_version;  /* The version of snes_mf used */

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

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

181: typedef struct _p_DMSNES  *DMSNES;
182: typedef struct _DMSNESOps *DMSNESOps;
183: struct _DMSNESOps {
184:   SNESFunctionFn *computefunction;
185:   SNESFunctionFn *computemffunction;
186:   SNESJacobianFn *computejacobian;

188:   /* objective */
189:   SNESObjectiveFn *computeobjective;

191:   /* Picard iteration functions */
192:   SNESFunctionFn *computepfunction;
193:   SNESJacobianFn *computepjacobian;

195:   /* User-defined smoother */
196:   SNESNGSFn *computegs;

198:   PetscErrorCode (*destroy)(DMSNES);
199:   PetscErrorCode (*duplicate)(DMSNES, DMSNES);
200: };

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

205:    Level: developer

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

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

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

218:    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`

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

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

242:   void *data;

244:   /* See developer note for DMSNES above */
245:   DM originaldm;
246: };
247: PETSC_EXTERN PetscErrorCode DMGetDMSNES(DM, DMSNES *);
248: PETSC_EXTERN PetscErrorCode DMSNESView(DMSNES, PetscViewer);
249: PETSC_EXTERN PetscErrorCode DMSNESLoad(DMSNES, PetscViewer);
250: PETSC_EXTERN PetscErrorCode DMGetDMSNESWrite(DM, DMSNES *);

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

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

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

295: PETSC_INTERN PetscErrorCode DMSNESUnsetFunctionContext_Internal(DM);
296: PETSC_EXTERN PetscErrorCode DMSNESUnsetJacobianContext_Internal(DM);

298: PETSC_EXTERN PetscLogEvent SNES_Solve;
299: PETSC_EXTERN PetscLogEvent SNES_SetUp;
300: PETSC_EXTERN PetscLogEvent SNES_FunctionEval;
301: PETSC_EXTERN PetscLogEvent SNES_JacobianEval;
302: PETSC_EXTERN PetscLogEvent SNES_NGSEval;
303: PETSC_EXTERN PetscLogEvent SNES_NGSFuncEval;
304: PETSC_EXTERN PetscLogEvent SNES_NewtonALEval;
305: PETSC_EXTERN PetscLogEvent SNES_NPCSolve;
306: PETSC_EXTERN PetscLogEvent SNES_ObjectiveEval;

308: PETSC_INTERN PetscBool  SNEScite;
309: PETSC_INTERN const char SNESCitation[];

311: /* Used by TAOBNK solvers */
312: PETSC_SINGLE_LIBRARY_VISIBILITY_INTERNAL PetscErrorCode KSPPostSolve_SNESEW(KSP, Vec, Vec, void *);
313: PETSC_SINGLE_LIBRARY_VISIBILITY_INTERNAL PetscErrorCode KSPPreSolve_SNESEW(KSP, Vec, Vec, void *);
314: PETSC_SINGLE_LIBRARY_VISIBILITY_INTERNAL PetscErrorCode SNESEWSetFromOptions_Private(SNESKSPEW *, PetscBool, MPI_Comm, const char *);

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

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

350: #define SNESCheckKSPSolve(snes) \
351:   do { \
352:     KSPConvergedReason kspreason; \
353:     PetscInt           lits; \
354:     PetscCall(KSPGetIterationNumber(snes->ksp, &lits)); \
355:     snes->linear_its += lits; \
356:     PetscCall(KSPGetConvergedReason(snes->ksp, &kspreason)); \
357:     if (kspreason < 0) { \
358:       if (kspreason == KSP_DIVERGED_NANORINF) { \
359:         PetscBool domainerror; \
360:         PetscCallMPI(MPIU_Allreduce(&snes->domainerror, &domainerror, 1, MPIU_BOOL, MPI_LOR, PetscObjectComm((PetscObject)snes))); \
361:         if (domainerror) snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; \
362:         else snes->reason = SNES_DIVERGED_LINEAR_SOLVE; \
363:         PetscFunctionReturn(PETSC_SUCCESS); \
364:       } else { \
365:         if (++snes->numLinearSolveFailures >= snes->maxLinearSolveFailures) { \
366:           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)); \
367:           snes->reason = SNES_DIVERGED_LINEAR_SOLVE; \
368:           PetscFunctionReturn(PETSC_SUCCESS); \
369:         } \
370:       } \
371:     } \
372:   } while (0)

374: #define SNESNeedNorm_Private(snes, iter) \
375:   (((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)) || \
376:    (snes)->normschedule == SNES_NORM_ALWAYS)