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: PetscErrorCode (*userdestroy)(void **);
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 *user; /* 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; /* preconditioner matrix */
57: Mat picard; /* copy of preconditioner matrix 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: PetscErrorCode (*monitordestroy[MAXSNESMONITORS])(void **); /* 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 */
82: PetscErrorCode (*reasonview[MAXSNESREASONVIEWS])(SNES, void *); /* snes converged reason view */
83: PetscErrorCode (*reasonviewdestroy[MAXSNESREASONVIEWS])(void **); /* reason view context destroy routine */
84: void *reasonviewcontext[MAXSNESREASONVIEWS]; /* reason view context */
85: PetscInt numberreasonviews; /* number of reason views */
86: PetscBool errorifnotconverged;
88: /* --- Routines and data that are unique to each particular solver --- */
90: PetscBool setupcalled; /* true if setup has been called */
91: void *data; /* implementation-specific data */
93: /* -------------------------- Parameters -------------------------------------- */
95: PetscInt max_its; /* max number of iterations */
96: PetscInt max_funcs; /* max number of function evals */
97: PetscInt nfuncs; /* number of function evaluations */
98: PetscInt iter; /* global iteration number */
99: PetscInt linear_its; /* total number of linear solver iterations */
100: PetscReal norm; /* residual norm of current iterate */
101: PetscReal ynorm; /* update norm of current iterate */
102: PetscReal xnorm; /* solution norm of current iterate */
103: PetscReal rtol; /* relative tolerance */
104: PetscReal divtol; /* relative divergence tolerance */
105: PetscReal abstol; /* absolute tolerance */
106: PetscReal stol; /* step length tolerance*/
107: PetscReal deltatol; /* trust region convergence tolerance */
108: PetscBool forceiteration; /* Force SNES to take at least one iteration regardless of the initial residual norm */
109: PetscInt lagpreconditioner; /* SNESSetLagPreconditioner() */
110: PetscInt lagjacobian; /* SNESSetLagJacobian() */
111: PetscInt jac_iter; /* The present iteration of the Jacobian lagging */
112: PetscBool lagjac_persist; /* The jac_iter persists until reset */
113: PetscInt pre_iter; /* The present iteration of the Preconditioner lagging */
114: PetscBool lagpre_persist; /* The pre_iter persists until reset */
115: PetscInt gridsequence; /* number of grid sequence steps to take; defaults to zero */
117: PetscBool tolerancesset; /* SNESSetTolerances() called and tolerances should persist through SNESCreate_XXX()*/
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: /* ------------------------- Miscellaneous Information ------------------------ */
131: PetscInt setfromoptionscalled;
132: PetscReal *conv_hist; /* If !0, stores function norm (or
133: gradient norm) at each iteration */
134: PetscInt *conv_hist_its; /* linear iterations for each Newton step */
135: size_t conv_hist_len; /* size of convergence history array */
136: size_t conv_hist_max; /* actual amount of data in conv_history */
137: PetscBool conv_hist_reset; /* reset counter for each new SNES solve */
138: PetscBool conv_hist_alloc;
139: PetscBool counters_reset; /* reset counter for each new SNES solve */
141: /* the next two are used for failures in the line search; they should be put elsewhere */
142: PetscInt numFailures; /* number of unsuccessful step attempts */
143: PetscInt maxFailures; /* maximum number of unsuccessful step attempts */
145: PetscInt numLinearSolveFailures;
146: PetscInt maxLinearSolveFailures;
148: PetscBool domainerror; /* set with SNESSetFunctionDomainError() */
149: PetscBool jacobiandomainerror; /* set with SNESSetJacobianDomainError() */
150: PetscBool checkjacdomainerror; /* if or not check Jacobian domain error after Jacobian evaluations */
152: PetscBool ksp_ewconv; /* flag indicating use of Eisenstat-Walker KSP convergence criteria */
153: void *kspconvctx; /* Eisenstat-Walker KSP convergence context */
155: /* SNESConvergedDefault context: split it off into a separate var/struct to be passed as context to SNESConvergedDefault? */
156: PetscReal ttol; /* rtol*initial_residual_norm */
157: PetscReal rnorm0; /* initial residual norm (used for divergence testing) */
159: Vec *vwork; /* more work vectors for Jacobian approx */
160: PetscInt nvwork;
162: PetscBool mf; /* -snes_mf OR -snes_mf_operator was used on this snes */
163: PetscBool mf_operator; /* -snes_mf_operator was used on this snes */
164: PetscInt mf_version; /* The version of snes_mf used */
166: PetscReal vizerotolerance; /* tolerance for considering an x[] value to be on the bound */
167: Vec xl, xu; /* upper and lower bounds for box constrained VI problems */
168: PetscInt ntruebounds; /* number of non-infinite bounds set for VI box constraints */
169: PetscBool usersetbounds; /* bounds have been set via SNESVISetVariableBounds(), rather than via computevariablebounds() callback. */
171: PetscBool alwayscomputesfinalresidual; /* Does SNESSolve_XXX always compute the value of the residual at the final
172: * solution and put it in vec_func? Used inside SNESSolve_FAS to determine
173: * if the final residual must be computed before restricting or prolonging
174: * it. */
175: };
177: typedef struct _p_DMSNES *DMSNES;
178: typedef struct _DMSNESOps *DMSNESOps;
179: struct _DMSNESOps {
180: SNESFunctionFn *computefunction;
181: SNESFunctionFn *computemffunction;
182: SNESJacobianFn *computejacobian;
184: /* objective */
185: SNESObjectiveFn *computeobjective;
187: /* Picard iteration functions */
188: SNESFunctionFn *computepfunction;
189: SNESJacobianFn *computepjacobian;
191: /* User-defined smoother */
192: SNESNGSFn *computegs;
194: PetscErrorCode (*destroy)(DMSNES);
195: PetscErrorCode (*duplicate)(DMSNES, DMSNES);
196: };
198: /*S
199: DMSNES - Object held by a `DM` that contains all the callback functions and their contexts needed by a `SNES`
201: Level: developer
203: Notes:
204: Users provides callback functions and their contexts to `SNES` using, for example, `SNESSetFunction()`. These values are stored
205: in a `DMSNES` that is contained in the `DM` associated with the `SNES`. If no `DM` was provided by
206: the user with `SNESSetDM()` it is automatically created by `SNESGetDM()` with `DMShellCreate()`.
208: Users very rarely need to worked directly with the `DMSNES` object, rather they work with the `SNES` and the `DM` they created
210: Multiple `DM` can share a single `DMSNES`, often each `DM` is associated with
211: a grid refinement level. `DMGetDMSNES()` returns the `DMSNES` associated with a `DM`. `DMGetDMSNESWrite()` returns a unique
212: `DMSNES` that is only associated with the current `DM`, making a copy of the shared `DMSNES` if needed (copy-on-write).
214: 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`
216: Developer Note:
217: The `originaldm` inside the `DMSNES` is NOT reference counted (to prevent a reference count loop between a `DM` and a `DMSNES`).
218: The `DM` on which this context was first created is cached here to implement one-way
219: copy-on-write. When `DMGetDMSNESWrite()` sees a request using a different `DM`, it makes a copy of the `TSDM`. Thus, if a user
220: only interacts directly with one level, e.g., using `TSSetIFunction()`, then coarse levels of a multilevel item
221: integrator are built, then the user changes the routine with another call to `TSSetIFunction()`, it automatically
222: propagates to all the levels. If instead, they get out a specific level and set the function on that level,
223: subsequent changes to the original level will no longer propagate to that level.
225: .seealso: [](ch_ts), `SNES`, `SNESCreate()`, `DM`, `DMGetDMSNESWrite()`, `DMGetDMSNES()`, `DMKSP`, `DMTS`, `DMSNESSetFunction()`, `DMSNESGetFunction()`,
226: `DMSNESSetFunctionContextDestroy()`, `DMSNESSetMFFunction()`, `DMSNESSetNGS()`, `DMSNESGetNGS()`, `DMSNESSetJacobian()`, `DMSNESGetJacobian()`,
227: `DMSNESSetJacobianContextDestroy()`, `DMSNESSetPicard()`, `DMSNESGetPicard()`, `DMSNESSetObjective()`, `DMSNESGetObjective()`, `DMCopyDMSNES()`
228: S*/
229: struct _p_DMSNES {
230: PETSCHEADER(struct _DMSNESOps);
231: PetscContainer functionctxcontainer;
232: PetscContainer jacobianctxcontainer;
233: void *mffunctionctx;
234: void *gsctx;
235: void *pctx;
236: void *objectivectx;
238: void *data;
240: /* See developer note for DMSNES above */
241: DM originaldm;
242: };
243: PETSC_EXTERN PetscErrorCode DMGetDMSNES(DM, DMSNES *);
244: PETSC_EXTERN PetscErrorCode DMSNESView(DMSNES, PetscViewer);
245: PETSC_EXTERN PetscErrorCode DMSNESLoad(DMSNES, PetscViewer);
246: PETSC_EXTERN PetscErrorCode DMGetDMSNESWrite(DM, DMSNES *);
248: /* Context for Eisenstat-Walker convergence criteria for KSP solvers */
249: typedef struct {
250: PetscInt version; /* flag indicating version (1,2,3 or 4) */
251: PetscReal rtol_0; /* initial rtol */
252: PetscReal rtol_last; /* last rtol */
253: PetscReal rtol_max; /* maximum rtol */
254: PetscReal gamma; /* mult. factor for version 2 rtol computation */
255: PetscReal alpha; /* power for version 2 rtol computation */
256: PetscReal alpha2; /* power for safeguard */
257: PetscReal threshold; /* threshold for imposing safeguard */
258: PetscReal lresid_last; /* linear residual from last iteration */
259: PetscReal norm_last; /* function norm from last iteration */
260: PetscReal norm_first; /* function norm from the beginning of the first iteration. */
261: PetscReal rtol_last_2, rk_last, rk_last_2;
262: PetscReal v4_p1, v4_p2, v4_p3, v4_m1, v4_m2, v4_m3, v4_m4;
263: } SNESKSPEW;
265: static inline PetscErrorCode SNESLogConvergenceHistory(SNES snes, PetscReal res, PetscInt its)
266: {
267: PetscFunctionBegin;
268: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
269: if (snes->conv_hist && snes->conv_hist_max > snes->conv_hist_len) {
270: if (snes->conv_hist) snes->conv_hist[snes->conv_hist_len] = res;
271: if (snes->conv_hist_its) snes->conv_hist_its[snes->conv_hist_len] = its;
272: snes->conv_hist_len++;
273: }
274: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
275: PetscFunctionReturn(PETSC_SUCCESS);
276: }
278: PETSC_EXTERN PetscErrorCode SNESVIProjectOntoBounds(SNES, Vec);
279: PETSC_INTERN PetscErrorCode SNESVICheckLocalMin_Private(SNES, Mat, Vec, Vec, PetscReal, PetscBool *);
280: PETSC_INTERN PetscErrorCode SNESReset_VI(SNES);
281: PETSC_INTERN PetscErrorCode SNESDestroy_VI(SNES);
282: PETSC_INTERN PetscErrorCode SNESView_VI(SNES, PetscViewer);
283: PETSC_INTERN PetscErrorCode SNESSetFromOptions_VI(SNES, PetscOptionItems *);
284: PETSC_INTERN PetscErrorCode SNESSetUp_VI(SNES);
285: PETSC_EXTERN_TYPEDEF typedef PetscErrorCode(SNESVIComputeVariableBoundsFn)(SNES, Vec, Vec);
286: PETSC_EXTERN_TYPEDEF typedef SNESVIComputeVariableBoundsFn *SNESVIComputeVariableBoundsFunction; // deprecated version
287: PETSC_INTERN PetscErrorCode SNESVISetComputeVariableBounds_VI(SNES, SNESVIComputeVariableBoundsFn);
288: PETSC_INTERN PetscErrorCode SNESVISetVariableBounds_VI(SNES, Vec, Vec);
289: PETSC_INTERN PetscErrorCode SNESConvergedDefault_VI(SNES, PetscInt, PetscReal, PetscReal, PetscReal, SNESConvergedReason *, void *);
291: PETSC_EXTERN PetscErrorCode DMSNESUnsetFunctionContext_Internal(DM);
292: PETSC_EXTERN PetscErrorCode DMSNESUnsetJacobianContext_Internal(DM);
293: PETSC_EXTERN PetscErrorCode DMSNESCheck_Internal(SNES, DM, Vec);
295: PETSC_EXTERN PetscLogEvent SNES_Solve;
296: PETSC_EXTERN PetscLogEvent SNES_SetUp;
297: PETSC_EXTERN PetscLogEvent SNES_FunctionEval;
298: PETSC_EXTERN PetscLogEvent SNES_JacobianEval;
299: PETSC_EXTERN PetscLogEvent SNES_NGSEval;
300: PETSC_EXTERN PetscLogEvent SNES_NGSFuncEval;
301: PETSC_EXTERN PetscLogEvent SNES_NPCSolve;
302: PETSC_EXTERN PetscLogEvent SNES_ObjectiveEval;
304: PETSC_INTERN PetscBool SNEScite;
305: PETSC_INTERN const char SNESCitation[];
307: /* Used by TAOBNK solvers */
308: PETSC_EXTERN PetscErrorCode KSPPostSolve_SNESEW(KSP, Vec, Vec, void *);
309: PETSC_EXTERN PetscErrorCode KSPPreSolve_SNESEW(KSP, Vec, Vec, void *);
310: PETSC_EXTERN PetscErrorCode SNESEWSetFromOptions_Private(SNESKSPEW *, PetscBool, MPI_Comm, const char *);
312: /*
313: Either generate an error or mark as diverged when a real from a SNES function norm is Nan or Inf.
314: domainerror is reset here, once reason is set, to allow subsequent iterations to be feasible (e.g. line search).
315: */
316: #define SNESCheckFunctionNorm(snes, beta) \
317: do { \
318: if (PetscIsInfOrNanReal(beta)) { \
319: PetscCheck(!snes->errorifnotconverged, PetscObjectComm((PetscObject)snes), PETSC_ERR_NOT_CONVERGED, "SNESSolve has not converged due to Nan or Inf norm"); \
320: { \
321: PetscBool domainerror; \
322: PetscCall(MPIU_Allreduce(&snes->domainerror, &domainerror, 1, MPIU_BOOL, MPI_LOR, PetscObjectComm((PetscObject)snes))); \
323: if (domainerror) { \
324: snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; \
325: snes->domainerror = PETSC_FALSE; \
326: } else snes->reason = SNES_DIVERGED_FNORM_NAN; \
327: PetscFunctionReturn(PETSC_SUCCESS); \
328: } \
329: } \
330: } while (0)
332: #define SNESCheckJacobianDomainerror(snes) \
333: do { \
334: if (snes->checkjacdomainerror) { \
335: PetscBool domainerror; \
336: PetscCall(MPIU_Allreduce(&snes->jacobiandomainerror, &domainerror, 1, MPIU_BOOL, MPI_LOR, PetscObjectComm((PetscObject)snes))); \
337: if (domainerror) { \
338: snes->reason = SNES_DIVERGED_JACOBIAN_DOMAIN; \
339: snes->jacobiandomainerror = PETSC_FALSE; \
340: PetscCheck(!snes->errorifnotconverged, PetscObjectComm((PetscObject)snes), PETSC_ERR_NOT_CONVERGED, "SNESSolve has not converged due to Jacobian domain error"); \
341: PetscFunctionReturn(PETSC_SUCCESS); \
342: } \
343: } \
344: } while (0)
346: #define SNESCheckKSPSolve(snes) \
347: do { \
348: KSPConvergedReason kspreason; \
349: PetscInt lits; \
350: PetscCall(KSPGetIterationNumber(snes->ksp, &lits)); \
351: snes->linear_its += lits; \
352: PetscCall(KSPGetConvergedReason(snes->ksp, &kspreason)); \
353: if (kspreason < 0) { \
354: if (kspreason == KSP_DIVERGED_NANORINF) { \
355: PetscBool domainerror; \
356: PetscCall(MPIU_Allreduce(&snes->domainerror, &domainerror, 1, MPIU_BOOL, MPI_LOR, PetscObjectComm((PetscObject)snes))); \
357: if (domainerror) snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; \
358: else snes->reason = SNES_DIVERGED_LINEAR_SOLVE; \
359: PetscFunctionReturn(PETSC_SUCCESS); \
360: } else { \
361: if (++snes->numLinearSolveFailures >= snes->maxLinearSolveFailures) { \
362: 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)); \
363: snes->reason = SNES_DIVERGED_LINEAR_SOLVE; \
364: PetscFunctionReturn(PETSC_SUCCESS); \
365: } \
366: } \
367: } \
368: } while (0)
370: #define SNESNeedNorm_Private(snes, iter) \
371: (((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)) || \
372: (snes)->normschedule == SNES_NORM_ALWAYS)