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. */

175: };

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

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

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

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

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

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

207:   void *data;

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

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

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

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

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

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

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

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

287: /*
288:     Either generate an error or mark as diverged when a real from a SNES function norm is Nan or Inf.
289:     domainerror is reset here, once reason is set, to allow subsequent iterations to be feasible (e.g. line search).
290: */
291: #define SNESCheckFunctionNorm(snes,beta) 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:   } } while (0)

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

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

340: #endif