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)(PetscOptionItems*,SNES); /* 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: PetscInt conv_hist_len; /* size of convergence history array */
135: PetscInt 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: void *functionctx;
201: void *mffunctionctx;
202: void *gsctx;
203: void *pctx;
204: void *jacobianctx;
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 or 2 of test */
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: } SNESKSPEW;
238: PETSC_STATIC_INLINE PetscErrorCode SNESLogConvergenceHistory(SNES snes,PetscReal res,PetscInt its)
239: {
243: PetscObjectSAWsTakeAccess((PetscObject)snes);
244: if (snes->conv_hist && snes->conv_hist_max > snes->conv_hist_len) {
245: if (snes->conv_hist) snes->conv_hist[snes->conv_hist_len] = res;
246: if (snes->conv_hist_its) snes->conv_hist_its[snes->conv_hist_len] = its;
247: snes->conv_hist_len++;
248: }
249: PetscObjectSAWsGrantAccess((PetscObject)snes);
250: return(0);
251: }
253: PETSC_EXTERN PetscErrorCode SNESVIProjectOntoBounds(SNES,Vec);
254: PETSC_INTERN PetscErrorCode SNESVICheckLocalMin_Private(SNES,Mat,Vec,Vec,PetscReal,PetscBool*);
255: PETSC_INTERN PetscErrorCode SNESReset_VI(SNES);
256: PETSC_INTERN PetscErrorCode SNESDestroy_VI(SNES);
257: PETSC_INTERN PetscErrorCode SNESView_VI(SNES,PetscViewer);
258: PETSC_INTERN PetscErrorCode SNESSetFromOptions_VI(PetscOptionItems*,SNES);
259: PETSC_INTERN PetscErrorCode SNESSetUp_VI(SNES);
260: PETSC_EXTERN_TYPEDEF typedef PetscErrorCode (*SNESVIComputeVariableBoundsFunction)(SNES,Vec,Vec);
261: PETSC_INTERN PetscErrorCode SNESVISetComputeVariableBounds_VI(SNES,SNESVIComputeVariableBoundsFunction);
262: PETSC_INTERN PetscErrorCode SNESVISetVariableBounds_VI(SNES,Vec,Vec);
263: PETSC_INTERN PetscErrorCode SNESConvergedDefault_VI(SNES,PetscInt,PetscReal,PetscReal,PetscReal,SNESConvergedReason*,void*);
265: PetscErrorCode SNESScaleStep_Private(SNES,Vec,PetscReal*,PetscReal*,PetscReal*,PetscReal*);
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: /*
282: Either generate an error or mark as diverged when a real from a SNES function norm is Nan or Inf.
283: domainerror is reset here, once reason is set, to allow subsequent iterations to be feasible (e.g. line search).
284: */
285: #define SNESCheckFunctionNorm(snes,beta) do { \
286: if (PetscIsInfOrNanReal(beta)) {\
287: if (snes->errorifnotconverged) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_NOT_CONVERGED,"SNESSolve has not converged due to Nan or Inf norm");\
288: else {\
289: PetscBool domainerror;\
290: PetscErrorCode MPIU_Allreduce(&snes->domainerror,&domainerror,1,MPIU_BOOL,MPI_LOR,PetscObjectComm((PetscObject)snes));\
291: if (domainerror) {\
292: snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN;\
293: snes->domainerror = PETSC_FALSE;\
294: } else snes->reason = SNES_DIVERGED_FNORM_NAN;\
295: return(0);\
296: }\
297: } } while (0)
299: #define SNESCheckJacobianDomainerror(snes) do { \
300: if (snes->checkjacdomainerror) {\
301: PetscBool domainerror;\
302: PetscErrorCode MPIU_Allreduce(&snes->jacobiandomainerror,&domainerror,1,MPIU_BOOL,MPI_LOR,PetscObjectComm((PetscObject)snes));\
303: if (domainerror) {\
304: snes->reason = SNES_DIVERGED_JACOBIAN_DOMAIN;\
305: if (snes->errorifnotconverged) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_NOT_CONVERGED,"SNESSolve has not converged due to Jacobian domain error");\
306: return(0);\
307: }\
308: } } while (0)
310: #define SNESCheckKSPSolve(snes)\
311: do {\
312: KSPConvergedReason kspreason; \
314: PetscInt lits; \
315: KSPGetIterationNumber(snes->ksp,&lits); \
316: snes->linear_its += lits; \
317: KSPGetConvergedReason(snes->ksp,&kspreason);\
318: if (kspreason < 0) {\
319: if (kspreason == KSP_DIVERGED_NANORINF) {\
320: PetscBool domainerror;\
321: MPIU_Allreduce(&snes->domainerror,&domainerror,1,MPIU_BOOL,MPI_LOR,PetscObjectComm((PetscObject)snes));\
322: if (domainerror) snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN;\
323: else snes->reason = SNES_DIVERGED_LINEAR_SOLVE; \
324: return(0);\
325: } else {\
326: if (++snes->numLinearSolveFailures >= snes->maxLinearSolveFailures) {\
327: PetscInfo3(snes,"iter=%D, number linear solve failures %D greater than current SNES allowed %D, stopping solve\n",snes->iter,snes->numLinearSolveFailures,snes->maxLinearSolveFailures);\
328: snes->reason = SNES_DIVERGED_LINEAR_SOLVE;\
329: return(0);\
330: }\
331: }\
332: }\
333: } while (0)
335: #endif