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

     83:   PetscViewer       convergedreasonviewer;
     84:   PetscViewerFormat convergedreasonformat;
     85:   PetscErrorCode (*reasonview[MAXSNESREASONVIEWS])(SNES, void *);   /* snes converged reason view */
     86:   PetscErrorCode (*reasonviewdestroy[MAXSNESREASONVIEWS])(void **); /* reason view context destroy routine */
     87:   void     *reasonviewcontext[MAXSNESREASONVIEWS];                  /* reason view context */
     88:   PetscInt  numberreasonviews;                                      /* number of reason views */
     89:   PetscBool errorifnotconverged;

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

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

     96:   /* --------------------------  Parameters -------------------------------------- */

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

    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_NewtonALEval;
    302: PETSC_EXTERN PetscLogEvent SNES_NPCSolve;
    303: PETSC_EXTERN PetscLogEvent SNES_ObjectiveEval;

    305: PETSC_INTERN PetscBool  SNEScite;
    306: PETSC_INTERN const char SNESCitation[];

    308: /* Used by TAOBNK solvers */
    309: PETSC_EXTERN PetscErrorCode KSPPostSolve_SNESEW(KSP, Vec, Vec, void *);
    310: PETSC_EXTERN PetscErrorCode KSPPreSolve_SNESEW(KSP, Vec, Vec, void *);
    311: PETSC_EXTERN PetscErrorCode SNESEWSetFromOptions_Private(SNESKSPEW *, PetscBool, MPI_Comm, const char *);

    313: /*
    314:     Either generate an error or mark as diverged when a real from a SNES function norm is Nan or Inf.
    315:     domainerror is reset here, once reason is set, to allow subsequent iterations to be feasible (e.g. line search).
    316: */
    317: #define SNESCheckFunctionNorm(snes, beta) \
    318:   do { \
    319:     if (PetscIsInfOrNanReal(beta)) { \
    320:       PetscCheck(!snes->errorifnotconverged, PetscObjectComm((PetscObject)snes), PETSC_ERR_NOT_CONVERGED, "SNESSolve has not converged due to Nan or Inf norm"); \
    321:       { \
    322:         PetscBool domainerror; \
    323:         PetscCallMPI(MPIU_Allreduce(&snes->domainerror, &domainerror, 1, MPIU_BOOL, MPI_LOR, PetscObjectComm((PetscObject)snes))); \
    324:         if (domainerror) { \
    325:           snes->reason      = SNES_DIVERGED_FUNCTION_DOMAIN; \
    326:           snes->domainerror = PETSC_FALSE; \
    327:         } else snes->reason = SNES_DIVERGED_FNORM_NAN; \
    328:         PetscFunctionReturn(PETSC_SUCCESS); \
    329:       } \
    330:     } \
    331:   } while (0)

    333: #define SNESCheckJacobianDomainerror(snes) \
    334:   do { \
    335:     if (snes->checkjacdomainerror) { \
    336:       PetscBool domainerror; \
    337:       PetscCallMPI(MPIU_Allreduce(&snes->jacobiandomainerror, &domainerror, 1, MPIU_BOOL, MPI_LOR, PetscObjectComm((PetscObject)snes))); \
    338:       if (domainerror) { \
    339:         snes->reason              = SNES_DIVERGED_JACOBIAN_DOMAIN; \
    340:         snes->jacobiandomainerror = PETSC_FALSE; \
    341:         PetscCheck(!snes->errorifnotconverged, PetscObjectComm((PetscObject)snes), PETSC_ERR_NOT_CONVERGED, "SNESSolve has not converged due to Jacobian domain error"); \
    342:         PetscFunctionReturn(PETSC_SUCCESS); \
    343:       } \
    344:     } \
    345:   } while (0)

    347: #define SNESCheckKSPSolve(snes) \
    348:   do { \
    349:     KSPConvergedReason kspreason; \
    350:     PetscInt           lits; \
    351:     PetscCall(KSPGetIterationNumber(snes->ksp, &lits)); \
    352:     snes->linear_its += lits; \
    353:     PetscCall(KSPGetConvergedReason(snes->ksp, &kspreason)); \
    354:     if (kspreason < 0) { \
    355:       if (kspreason == KSP_DIVERGED_NANORINF) { \
    356:         PetscBool domainerror; \
    357:         PetscCallMPI(MPIU_Allreduce(&snes->domainerror, &domainerror, 1, MPIU_BOOL, MPI_LOR, PetscObjectComm((PetscObject)snes))); \
    358:         if (domainerror) snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; \
    359:         else snes->reason = SNES_DIVERGED_LINEAR_SOLVE; \
    360:         PetscFunctionReturn(PETSC_SUCCESS); \
    361:       } else { \
    362:         if (++snes->numLinearSolveFailures >= snes->maxLinearSolveFailures) { \
    363:           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)); \
    364:           snes->reason = SNES_DIVERGED_LINEAR_SOLVE; \
    365:           PetscFunctionReturn(PETSC_SUCCESS); \
    366:         } \
    367:       } \
    368:     } \
    369:   } while (0)

    371: #define SNESNeedNorm_Private(snes, iter) \
    372:   (((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)) || \
    373:    (snes)->normschedule == SNES_NORM_ALWAYS)