Actual source code: snesimpl.h

  1: #ifndef __SNESIMPL_H

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

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

 10: typedef struct _SNESOps *SNESOps;

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

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

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

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

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

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

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

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

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

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

 70:   /* ------------------------Time stepping hooks-----------------------------------*/

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

 74:   PetscErrorCode (*monitor[MAXSNESMONITORS])(SNES, PetscInt, PetscReal, void *); /* monitor routine */
 75:   PetscErrorCode (*monitordestroy[MAXSNESMONITORS])(void **);                    /* monitor context destroy routine */
 76:   void               *monitorcontext[MAXSNESMONITORS];                           /* monitor context */
 77:   PetscInt            numbermonitors;                                            /* number of monitors */
 78:   PetscBool           pauseFinal;                                                /* pause all drawing monitor at the final iterate */
 79:   void               *cnvP;                                                      /* convergence context */
 80:   SNESConvergedReason reason;                                                    /* converged reason */
 81:   PetscErrorCode (*reasonview[MAXSNESREASONVIEWS])(SNES, void *);                /* snes converged reason view */
 82:   PetscErrorCode (*reasonviewdestroy[MAXSNESREASONVIEWS])(void **);              /* reason view context destroy routine */
 83:   void     *reasonviewcontext[MAXSNESREASONVIEWS];                               /* reason view context */
 84:   PetscInt  numberreasonviews;                                                   /* number of reason views */
 85:   PetscBool errorifnotconverged;

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

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

 92:   /* --------------------------  Parameters -------------------------------------- */

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

116:   PetscBool tolerancesset; /* SNESSetTolerances() called and tolerances should persist through SNESCreate_XXX()*/

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

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

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

125:   PetscInt nwork;
126:   Vec     *work;

128:   /* ------------------------- Miscellaneous Information ------------------------ */

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

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

144:   PetscInt numLinearSolveFailures;
145:   PetscInt maxLinearSolveFailures;

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

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

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

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

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

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

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

176: typedef struct _p_DMSNES  *DMSNES;
177: typedef struct _DMSNESOps *DMSNESOps;
178: struct _DMSNESOps {
179:   PetscErrorCode (*computefunction)(SNES, Vec, Vec, void *);
180:   PetscErrorCode (*computemffunction)(SNES, Vec, Vec, void *);
181:   PetscErrorCode (*computejacobian)(SNES, Vec, Mat, Mat, void *);

183:   /* objective */
184:   PetscErrorCode (*computeobjective)(SNES, Vec, PetscReal *, void *);

186:   /* Picard iteration functions */
187:   PetscErrorCode (*computepfunction)(SNES, Vec, Vec, void *);
188:   PetscErrorCode (*computepjacobian)(SNES, Vec, Mat, Mat, void *);

190:   /* User-defined smoother */
191:   PetscErrorCode (*computegs)(SNES, Vec, Vec, void *);

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

197: struct _p_DMSNES {
198:   PETSCHEADER(struct _DMSNESOps);
199:   PetscContainer functionctxcontainer;
200:   PetscContainer jacobianctxcontainer;
201:   void          *mffunctionctx;
202:   void          *gsctx;
203:   void          *pctx;
204:   void          *objectivectx;

206:   void *data;

208:   /* This is NOT reference counted. The DM on which this context was first created is cached here to implement one-way
209:    * copy-on-write. When DMGetDMSNESWrite() sees a request using a different DM, it makes a copy. Thus, if a user
210:    * only interacts directly with one level, e.g., using SNESSetFunction(), then SNESSetUp_FAS() is called to build
211:    * coarse levels, then the user changes the routine with another call to SNESSetFunction(), it automatically
212:    * propagates to all the levels. If instead, they get out a specific level and set the function on that level,
213:    * subsequent changes to the original level will no longer propagate to that level.
214:    */
215:   DM originaldm;
216: };
217: PETSC_EXTERN PetscErrorCode DMGetDMSNES(DM, DMSNES *);
218: PETSC_EXTERN PetscErrorCode DMSNESView(DMSNES, PetscViewer);
219: PETSC_EXTERN PetscErrorCode DMSNESLoad(DMSNES, PetscViewer);
220: PETSC_EXTERN PetscErrorCode DMGetDMSNESWrite(DM, DMSNES *);

222: /* Context for Eisenstat-Walker convergence criteria for KSP solvers */
223: typedef struct {
224:   PetscInt  version;     /* flag indicating version (1,2,3 or 4) */
225:   PetscReal rtol_0;      /* initial rtol */
226:   PetscReal rtol_last;   /* last rtol */
227:   PetscReal rtol_max;    /* maximum rtol */
228:   PetscReal gamma;       /* mult. factor for version 2 rtol computation */
229:   PetscReal alpha;       /* power for version 2 rtol computation */
230:   PetscReal alpha2;      /* power for safeguard */
231:   PetscReal threshold;   /* threshold for imposing safeguard */
232:   PetscReal lresid_last; /* linear residual from last iteration */
233:   PetscReal norm_last;   /* function norm from last iteration */
234:   PetscReal norm_first;  /* function norm from the beginning of the first iteration. */
235:   PetscReal rtol_last_2, rk_last, rk_last_2;
236:   PetscReal v4_p1, v4_p2, v4_p3, v4_m1, v4_m2, v4_m3, v4_m4;
237: } SNESKSPEW;

239: static inline PetscErrorCode SNESLogConvergenceHistory(SNES snes, PetscReal res, PetscInt its)
240: {
241:   PetscObjectSAWsTakeAccess((PetscObject)snes);
242:   if (snes->conv_hist && snes->conv_hist_max > snes->conv_hist_len) {
243:     if (snes->conv_hist) snes->conv_hist[snes->conv_hist_len] = res;
244:     if (snes->conv_hist_its) snes->conv_hist_its[snes->conv_hist_len] = its;
245:     snes->conv_hist_len++;
246:   }
247:   PetscObjectSAWsGrantAccess((PetscObject)snes);
248:   return 0;
249: }

251: PETSC_EXTERN PetscErrorCode SNESVIProjectOntoBounds(SNES, Vec);
252: PETSC_INTERN PetscErrorCode SNESVICheckLocalMin_Private(SNES, Mat, Vec, Vec, PetscReal, PetscBool *);
253: PETSC_INTERN PetscErrorCode SNESReset_VI(SNES);
254: PETSC_INTERN PetscErrorCode SNESDestroy_VI(SNES);
255: PETSC_INTERN PetscErrorCode SNESView_VI(SNES, PetscViewer);
256: PETSC_INTERN PetscErrorCode SNESSetFromOptions_VI(SNES, PetscOptionItems *);
257: PETSC_INTERN PetscErrorCode SNESSetUp_VI(SNES);
258: PETSC_EXTERN_TYPEDEF typedef PetscErrorCode (*SNESVIComputeVariableBoundsFunction)(SNES, Vec, Vec);
259: PETSC_INTERN PetscErrorCode SNESVISetComputeVariableBounds_VI(SNES, SNESVIComputeVariableBoundsFunction);
260: PETSC_INTERN PetscErrorCode SNESVISetVariableBounds_VI(SNES, Vec, Vec);
261: PETSC_INTERN PetscErrorCode SNESConvergedDefault_VI(SNES, PetscInt, PetscReal, PetscReal, PetscReal, SNESConvergedReason *, void *);

263: PetscErrorCode              SNESScaleStep_Private(SNES, Vec, PetscReal *, PetscReal *, PetscReal *, PetscReal *);
264: PETSC_EXTERN PetscErrorCode DMSNESUnsetFunctionContext_Internal(DM);
265: PETSC_EXTERN PetscErrorCode DMSNESUnsetJacobianContext_Internal(DM);
266: PETSC_EXTERN PetscErrorCode DMSNESCheck_Internal(SNES, DM, Vec);

268: PETSC_EXTERN PetscLogEvent SNES_Solve;
269: PETSC_EXTERN PetscLogEvent SNES_Setup;
270: PETSC_EXTERN PetscLogEvent SNES_LineSearch;
271: PETSC_EXTERN PetscLogEvent SNES_FunctionEval;
272: PETSC_EXTERN PetscLogEvent SNES_JacobianEval;
273: PETSC_EXTERN PetscLogEvent SNES_NGSEval;
274: PETSC_EXTERN PetscLogEvent SNES_NGSFuncEval;
275: PETSC_EXTERN PetscLogEvent SNES_NPCSolve;
276: PETSC_EXTERN PetscLogEvent SNES_ObjectiveEval;

278: PETSC_INTERN PetscBool  SNEScite;
279: PETSC_INTERN const char SNESCitation[];

281: /* Used by TAOBNK solvers */
282: PETSC_EXTERN PetscErrorCode KSPPostSolve_SNESEW(KSP, Vec, Vec, SNES);
283: PETSC_EXTERN PetscErrorCode KSPPreSolve_SNESEW(KSP, Vec, Vec, SNES);
284: PETSC_EXTERN PetscErrorCode SNESEWSetFromOptions_Private(SNESKSPEW *, MPI_Comm, const char *);

286: /*
287:     Either generate an error or mark as diverged when a real from a SNES function norm is Nan or Inf.
288:     domainerror is reset here, once reason is set, to allow subsequent iterations to be feasible (e.g. line search).
289: */
290: #define SNESCheckFunctionNorm(snes, beta) \
291:   do { \
292:     if (PetscIsInfOrNanReal(beta)) { \
294:       { \
295:         PetscBool domainerror; \
296:         MPIU_Allreduce(&snes->domainerror, &domainerror, 1, MPIU_BOOL, MPI_LOR, PetscObjectComm((PetscObject)snes)); \
297:         if (domainerror) { \
298:           snes->reason      = SNES_DIVERGED_FUNCTION_DOMAIN; \
299:           snes->domainerror = PETSC_FALSE; \
300:         } else snes->reason = SNES_DIVERGED_FNORM_NAN; \
301:         return 0; \
302:       } \
303:     } \
304:   } while (0)

306: #define SNESCheckJacobianDomainerror(snes) \
307:   do { \
308:     if (snes->checkjacdomainerror) { \
309:       PetscBool domainerror; \
310:       MPIU_Allreduce(&snes->jacobiandomainerror, &domainerror, 1, MPIU_BOOL, MPI_LOR, PetscObjectComm((PetscObject)snes)); \
311:       if (domainerror) { \
312:         snes->reason = SNES_DIVERGED_JACOBIAN_DOMAIN; \
314:         return 0; \
315:       } \
316:     } \
317:   } while (0)

319: #define SNESCheckKSPSolve(snes) \
320:   do { \
321:     KSPConvergedReason kspreason; \
322:     PetscInt           lits; \
323:     KSPGetIterationNumber(snes->ksp, &lits); \
324:     snes->linear_its += lits; \
325:     KSPGetConvergedReason(snes->ksp, &kspreason); \
326:     if (kspreason < 0) { \
327:       if (kspreason == KSP_DIVERGED_NANORINF) { \
328:         PetscBool domainerror; \
329:         MPIU_Allreduce(&snes->domainerror, &domainerror, 1, MPIU_BOOL, MPI_LOR, PetscObjectComm((PetscObject)snes)); \
330:         if (domainerror) snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; \
331:         else snes->reason = SNES_DIVERGED_LINEAR_SOLVE; \
332:         return 0; \
333:       } else { \
334:         if (++snes->numLinearSolveFailures >= snes->maxLinearSolveFailures) { \
335:           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); \
336:           snes->reason = SNES_DIVERGED_LINEAR_SOLVE; \
337:           return 0; \
338:         } \
339:       } \
340:     } \
341:   } while (0)

343: #endif