1 #ifndef __SNESIMPL_H
2 #define __SNESIMPL_H
3
4 #include <petscsnes.h>
5 #include <petsc/private/petscimpl.h>
6
7 PETSC_EXTERN PetscBool SNESRegisterAllCalled;
8 PETSC_EXTERN PetscErrorCode SNESRegisterAll(void);
9
10 typedef struct _SNESOps *SNESOps;
11
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 };
31
32 /*
33 Nonlinear solver context
34 */
35 #define MAXSNESMONITORS 5
36
37 struct _p_SNES {
38 PETSCHEADER(struct _SNESOps);
39 DM dm;
40 PetscBool dmAuto; /* SNES created currently used DM automatically */
41 SNES npc;
42 PCSide npcside;
43 PetscBool usesnpc; /* type can use a nonlinear preconditioner */
44
45 /* ------------------------ User-provided stuff -------------------------------*/
46 void *user; /* user-defined context */
47
48 Vec vec_rhs; /* If non-null, solve F(x) = rhs */
49 Vec vec_sol; /* pointer to solution */
50
51 Vec vec_func; /* pointer to function */
52
53 Mat jacobian; /* Jacobian matrix */
54 Mat jacobian_pre; /* preconditioner matrix */
55 void *initialguessP; /* user-defined initial guess context */
56 KSP ksp; /* linear solver context */
57 SNESLineSearch linesearch; /* line search context */
58 PetscBool usesksp;
59 MatStructure matstruct; /* Used by Picard solver */
60
61 Vec vec_sol_update; /* pointer to solution update */
62
63 Vec scaling; /* scaling vector */
64 void *scaP; /* scaling context */
65
66 PetscReal precheck_picard_angle; /* For use with SNESLineSearchPreCheckPicard */
67
68 /* ------------------------Time stepping hooks-----------------------------------*/
69
70 /* ---------------- PETSc-provided (or user-provided) stuff ---------------------*/
71
72 PetscErrorCode (*monitor[MAXSNESMONITORS])(SNES,PetscInt,PetscReal,void*); /* monitor routine */
73 PetscErrorCode (*monitordestroy[MAXSNESMONITORS])(void**); /* monitor context destroy routine */
74 void *monitorcontext[MAXSNESMONITORS]; /* monitor context */
75 PetscInt numbermonitors; /* number of monitors */
76 void *cnvP; /* convergence context */
77 SNESConvergedReason reason;
78 PetscBool errorifnotconverged;
79
80 /* --- Routines and data that are unique to each particular solver --- */
81
82 PetscBool setupcalled; /* true if setup has been called */
83 void *data; /* implementation-specific data */
84
85 /* -------------------------- Parameters -------------------------------------- */
86
87 PetscInt max_its; /* max number of iterations */
88 PetscInt max_funcs; /* max number of function evals */
89 PetscInt nfuncs; /* number of function evaluations */
90 PetscInt iter; /* global iteration number */
91 PetscInt linear_its; /* total number of linear solver iterations */
92 PetscReal norm; /* residual norm of current iterate */
93 PetscReal ynorm; /* update norm of current iterate */
94 PetscReal xnorm; /* solution norm of current iterate */
95 PetscReal rtol; /* relative tolerance */
96 PetscReal divtol; /* relative divergence tolerance */
97 PetscReal abstol; /* absolute tolerance */
98 PetscReal stol; /* step length tolerance*/
99 PetscReal deltatol; /* trust region convergence tolerance */
100 PetscBool forceiteration; /* Force SNES to take at least one iteration regardless of the initial residual norm */
101 PetscInt lagpreconditioner; /* SNESSetLagPreconditioner() */
102 PetscInt lagjacobian; /* SNESSetLagJacobian() */
103 PetscInt jac_iter; /* The present iteration of the Jacobian lagging */
104 PetscBool lagjac_persist; /* The jac_iter persists until reset */
105 PetscInt pre_iter; /* The present iteration of the Preconditioner lagging */
106 PetscBool lagpre_persist; /* The pre_iter persists until reset */
107 PetscInt gridsequence; /* number of grid sequence steps to take; defaults to zero */
108
109 PetscBool tolerancesset; /* SNESSetTolerances() called and tolerances should persist through SNESCreate_XXX()*/
110
111 PetscBool vec_func_init_set; /* the initial function has been set */
112
113 SNESNormSchedule normschedule; /* Norm computation type for SNES instance */
114 SNESFunctionType functype; /* Function type for the SNES instance */
115
116 /* ------------------------ Default work-area management ---------------------- */
117
118 PetscInt nwork;
119 Vec *work;
120
121 /* ------------------------- Miscellaneous Information ------------------------ */
122
123 PetscInt setfromoptionscalled;
124 PetscReal *conv_hist; /* If !0, stores function norm (or
125 gradient norm) at each iteration */
126 PetscInt *conv_hist_its; /* linear iterations for each Newton step */
127 PetscInt conv_hist_len; /* size of convergence history array */
128 PetscInt conv_hist_max; /* actual amount of data in conv_history */
129 PetscBool conv_hist_reset; /* reset counter for each new SNES solve */
130 PetscBool conv_hist_alloc;
131 PetscBool counters_reset; /* reset counter for each new SNES solve */
132
133 /* the next two are used for failures in the line search; they should be put elsewhere */
134 PetscInt numFailures; /* number of unsuccessful step attempts */
135 PetscInt maxFailures; /* maximum number of unsuccessful step attempts */
136
137 PetscInt numLinearSolveFailures;
138 PetscInt maxLinearSolveFailures;
139
140 PetscBool domainerror; /* set with SNESSetFunctionDomainError() */
141 PetscBool jacobiandomainerror; /* set with SNESSetJacobianDomainError() */
142 PetscBool checkjacdomainerror; /* if or not check Jacobian domain error after Jacobian evaluations */
143
144 PetscBool ksp_ewconv; /* flag indicating use of Eisenstat-Walker KSP convergence criteria */
145 void *kspconvctx; /* Eisenstat-Walker KSP convergence context */
146
147 /* SNESConvergedDefault context: split it off into a separate var/struct to be passed as context to SNESConvergedDefault? */
148 PetscReal ttol; /* rtol*initial_residual_norm */
149 PetscReal rnorm0; /* initial residual norm (used for divergence testing) */
150
151 Vec *vwork; /* more work vectors for Jacobian approx */
152 PetscInt nvwork;
153
154 PetscBool mf; /* -snes_mf was used on this snes */
155 PetscBool mf_operator; /* -snes_mf_operator was used on this snes */
156 PetscInt mf_version; /* The version of snes_mf used */
157
158 PetscReal vizerotolerance; /* tolerance for considering an x[] value to be on the bound */
159 Vec xl,xu; /* upper and lower bounds for box constrained VI problems */
160 PetscInt ntruebounds; /* number of non-infinite bounds set for VI box constraints */
161 PetscBool usersetbounds; /* bounds have been set via SNESVISetVariableBounds(), rather than via computevariablebounds() callback. */
162
163 PetscBool alwayscomputesfinalresidual; /* Does SNESSolve_XXX always compute the value of the residual at the final
164 * solution and put it in vec_func? Used inside SNESSolve_FAS to determine
165 * if the final residual must be computed before restricting or prolonging
166 * it. */
167
168 };
169
170 typedef struct _p_DMSNES *DMSNES;
171 typedef struct _DMSNESOps *DMSNESOps;
172 struct _DMSNESOps {
173 PetscErrorCode (*computefunction)(SNES,Vec,Vec,void*);
174 PetscErrorCode (*computejacobian)(SNES,Vec,Mat,Mat,void*);
175
176 /* objective */
177 PetscErrorCode (*computeobjective)(SNES,Vec,PetscReal*,void*);
178
179 /* Picard iteration functions */
180 PetscErrorCode (*computepfunction)(SNES,Vec,Vec,void*);
181 PetscErrorCode (*computepjacobian)(SNES,Vec,Mat,Mat,void*);
182
183 /* User-defined smoother */
184 PetscErrorCode (*computegs)(SNES,Vec,Vec,void*);
185
186 PetscErrorCode (*destroy)(DMSNES);
187 PetscErrorCode (*duplicate)(DMSNES,DMSNES);
188 };
189
190 struct _p_DMSNES {
191 PETSCHEADER(struct _DMSNESOps);
192 void *functionctx;
193 void *gsctx;
194 void *pctx;
195 void *jacobianctx;
196 void *objectivectx;
197
198 void *data;
199
200 /* This is NOT reference counted. The DM on which this context was first created is cached here to implement one-way
201 * copy-on-write. When DMGetDMSNESWrite() sees a request using a different DM, it makes a copy. Thus, if a user
202 * only interacts directly with one level, e.g., using SNESSetFunction(), then SNESSetUp_FAS() is called to build
203 * coarse levels, then the user changes the routine with another call to SNESSetFunction(), it automatically
204 * propagates to all the levels. If instead, they get out a specific level and set the function on that level,
205 * subsequent changes to the original level will no longer propagate to that level.
206 */
207 DM originaldm;
208 };
209 PETSC_EXTERN PetscErrorCode DMGetDMSNES(DM,DMSNES*);
210 PETSC_EXTERN PetscErrorCode DMSNESView(DMSNES,PetscViewer);
211 PETSC_EXTERN PetscErrorCode DMSNESLoad(DMSNES,PetscViewer);
212 PETSC_EXTERN PetscErrorCode DMGetDMSNESWrite(DM,DMSNES*);
213
214
215 /* Context for Eisenstat-Walker convergence criteria for KSP solvers */
216 typedef struct {
217 PetscInt version; /* flag indicating version 1 or 2 of test */
218 PetscReal rtol_0; /* initial rtol */
219 PetscReal rtol_last; /* last rtol */
220 PetscReal rtol_max; /* maximum rtol */
221 PetscReal gamma; /* mult. factor for version 2 rtol computation */
222 PetscReal alpha; /* power for version 2 rtol computation */
223 PetscReal alpha2; /* power for safeguard */
224 PetscReal threshold; /* threshold for imposing safeguard */
225 PetscReal lresid_last; /* linear residual from last iteration */
226 PetscReal norm_last; /* function norm from last iteration */
227 PetscReal norm_first; /* function norm from the beginning of the first iteration. */
228 } SNESKSPEW;
229
SNESLogConvergenceHistory(SNES snes,PetscReal res,PetscInt its)230 PETSC_STATIC_INLINE PetscErrorCode SNESLogConvergenceHistory(SNES snes,PetscReal res,PetscInt its)
231 {
232 PetscErrorCode ierr;
233
234 PetscFunctionBegin;
235 ierr = PetscObjectSAWsTakeAccess((PetscObject)snes);CHKERRQ(ierr);
236 if (snes->conv_hist && snes->conv_hist_max > snes->conv_hist_len) {
237 if (snes->conv_hist) snes->conv_hist[snes->conv_hist_len] = res;
238 if (snes->conv_hist_its) snes->conv_hist_its[snes->conv_hist_len] = its;
239 snes->conv_hist_len++;
240 }
241 ierr = PetscObjectSAWsGrantAccess((PetscObject)snes);CHKERRQ(ierr);
242 PetscFunctionReturn(0);
243 }
244
245 PETSC_EXTERN PetscErrorCode SNESVIProjectOntoBounds(SNES,Vec);
246 PETSC_INTERN PetscErrorCode SNESVICheckLocalMin_Private(SNES,Mat,Vec,Vec,PetscReal,PetscBool*);
247 PETSC_INTERN PetscErrorCode SNESReset_VI(SNES);
248 PETSC_INTERN PetscErrorCode SNESDestroy_VI(SNES);
249 PETSC_INTERN PetscErrorCode SNESView_VI(SNES,PetscViewer);
250 PETSC_INTERN PetscErrorCode SNESSetFromOptions_VI(PetscOptionItems*,SNES);
251 PETSC_INTERN PetscErrorCode SNESSetUp_VI(SNES);
252 PETSC_EXTERN_TYPEDEF typedef PetscErrorCode (*SNESVIComputeVariableBoundsFunction)(SNES,Vec,Vec);
253 PETSC_INTERN PetscErrorCode SNESVISetComputeVariableBounds_VI(SNES,SNESVIComputeVariableBoundsFunction);
254 PETSC_INTERN PetscErrorCode SNESVISetVariableBounds_VI(SNES,Vec,Vec);
255 PETSC_INTERN PetscErrorCode SNESConvergedDefault_VI(SNES,PetscInt,PetscReal,PetscReal,PetscReal,SNESConvergedReason*,void*);
256
257 PetscErrorCode SNESScaleStep_Private(SNES,Vec,PetscReal*,PetscReal*,PetscReal*,PetscReal*);
258 PETSC_EXTERN PetscErrorCode DMSNESCheck_Internal(SNES,DM,Vec);
259
260 PETSC_EXTERN PetscLogEvent SNES_Solve;
261 PETSC_EXTERN PetscLogEvent SNES_Setup;
262 PETSC_EXTERN PetscLogEvent SNES_LineSearch;
263 PETSC_EXTERN PetscLogEvent SNES_FunctionEval;
264 PETSC_EXTERN PetscLogEvent SNES_JacobianEval;
265 PETSC_EXTERN PetscLogEvent SNES_NGSEval;
266 PETSC_EXTERN PetscLogEvent SNES_NGSFuncEval;
267 PETSC_EXTERN PetscLogEvent SNES_NPCSolve;
268 PETSC_EXTERN PetscLogEvent SNES_ObjectiveEval;
269
270 PETSC_INTERN PetscBool SNEScite;
271 PETSC_INTERN const char SNESCitation[];
272
273 /*
274 Either generate an error or mark as diverged when a real from a SNES function norm is Nan or Inf.
275 domainerror is reset here, once reason is set, to allow subsequent iterations to be feasible (e.g. line search).
276 */
277 #define SNESCheckFunctionNorm(snes,beta) do { \
278 if (PetscIsInfOrNanReal(beta)) {\
279 if (snes->errorifnotconverged) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_NOT_CONVERGED,"SNESSolve has not converged due to Nan or Inf norm");\
280 else {\
281 PetscBool domainerror;\
282 PetscErrorCode ierr = MPIU_Allreduce(&snes->domainerror,&domainerror,1,MPIU_BOOL,MPI_LOR,PetscObjectComm((PetscObject)snes));CHKERRQ(ierr);\
283 if (domainerror) {\
284 snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN;\
285 snes->domainerror = PETSC_FALSE;\
286 } else snes->reason = SNES_DIVERGED_FNORM_NAN;\
287 PetscFunctionReturn(0);\
288 }\
289 } } while (0)
290
291 #define SNESCheckJacobianDomainerror(snes) do { \
292 if (snes->checkjacdomainerror) {\
293 PetscBool domainerror;\
294 PetscErrorCode ierr = MPIU_Allreduce(&snes->jacobiandomainerror,&domainerror,1,MPIU_BOOL,MPI_LOR,PetscObjectComm((PetscObject)snes));CHKERRQ(ierr);\
295 if (domainerror) {\
296 snes->reason = SNES_DIVERGED_JACOBIAN_DOMAIN;\
297 if (snes->errorifnotconverged) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_NOT_CONVERGED,"SNESSolve has not converged due to Jacobian domain error");\
298 PetscFunctionReturn(0);\
299 }\
300 } } while (0)
301
302 #define SNESCheckKSPSolve(snes)\
303 do {\
304 KSPConvergedReason kspreason; \
305 PetscErrorCode ierr; \
306 PetscInt lits; \
307 ierr = KSPGetIterationNumber(snes->ksp,&lits);CHKERRQ(ierr); \
308 snes->linear_its += lits; \
309 ierr = KSPGetConvergedReason(snes->ksp,&kspreason);CHKERRQ(ierr);\
310 if (kspreason < 0) {\
311 if (kspreason == KSP_DIVERGED_NANORINF) {\
312 PetscBool domainerror;\
313 ierr = MPIU_Allreduce(&snes->domainerror,&domainerror,1,MPIU_BOOL,MPI_LOR,PetscObjectComm((PetscObject)snes));CHKERRQ(ierr);\
314 if (domainerror) snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN;\
315 else snes->reason = SNES_DIVERGED_LINEAR_SOLVE; \
316 PetscFunctionReturn(0);\
317 } else {\
318 if (++snes->numLinearSolveFailures >= snes->maxLinearSolveFailures) {\
319 ierr = PetscInfo3(snes,"iter=%D, number linear solve failures %D greater than current SNES allowed %D, stopping solve\n",snes->iter,snes->numLinearSolveFailures,snes->maxLinearSolveFailures);CHKERRQ(ierr);\
320 snes->reason = SNES_DIVERGED_LINEAR_SOLVE;\
321 PetscFunctionReturn(0);\
322 }\
323 }\
324 }\
325 } while (0)
326
327 #endif
328