Actual source code: snes.c

  1: #include <petsc/private/snesimpl.h>
  2: #include <petscdmshell.h>
  3: #include <petscdraw.h>
  4: #include <petscds.h>
  5: #include <petscdmadaptor.h>
  6: #include <petscconvest.h>

  8: PetscBool         SNESRegisterAllCalled = PETSC_FALSE;
  9: PetscFunctionList SNESList              = NULL;

 11: /* Logging support */
 12: PetscClassId  SNES_CLASSID, DMSNES_CLASSID;
 13: PetscLogEvent SNES_Solve, SNES_Setup, SNES_FunctionEval, SNES_JacobianEval, SNES_NGSEval, SNES_NGSFuncEval, SNES_NPCSolve, SNES_ObjectiveEval;

 15: /*@
 16:    SNESSetErrorIfNotConverged - Causes `SNESSolve()` to generate an error if the solver has not converged.

 18:    Logically Collective on snes

 20:    Input Parameters:
 21: +  snes - iterative context obtained from `SNESCreate()`
 22: -  flg - `PETSC_TRUE` indicates you want the error generated

 24:    Options database keys:
 25: .  -snes_error_if_not_converged <true,false> - cause an immediate error condition and stop the program if the solver does not converge

 27:    Level: intermediate

 29:    Note:
 30:    Normally PETSc continues if a solver fails to converge, you can call `SNESGetConvergedReason()` after a `SNESSolve()`
 31:    to determine if it has converged. Otherwise the solution may be inaccurate or wrong

 33: .seealso: `SNES`, `SNESGetErrorIfNotConverged()`, `KSPGetErrorIfNotConverged()`, `KSPSetErrorIfNotConverged()`
 34: @*/
 35: PetscErrorCode SNESSetErrorIfNotConverged(SNES snes, PetscBool flg)
 36: {
 39:   snes->errorifnotconverged = flg;
 40:   return 0;
 41: }

 43: /*@
 44:    SNESGetErrorIfNotConverged - Indicates if `SNESSolve()` will generate an error if the solver does not converge?

 46:    Not Collective

 48:    Input Parameter:
 49: .  snes - iterative context obtained from `SNESCreate()`

 51:    Output Parameter:
 52: .  flag - `PETSC_TRUE` if it will generate an error, else `PETSC_FALSE`

 54:    Level: intermediate

 56: .seealso: `SNES`, `SNESSolve()`, `SNESSetErrorIfNotConverged()`, `KSPGetErrorIfNotConverged()`, `KSPSetErrorIfNotConverged()`
 57: @*/
 58: PetscErrorCode SNESGetErrorIfNotConverged(SNES snes, PetscBool *flag)
 59: {
 62:   *flag = snes->errorifnotconverged;
 63:   return 0;
 64: }

 66: /*@
 67:     SNESSetAlwaysComputesFinalResidual - tells the `SNES` to always compute the residual at the final solution

 69:    Logically Collective on snes

 71:     Input Parameters:
 72: +   snes - the shell `SNES`
 73: -   flg - `PETSC_TRUE` to always compute the residual

 75:    Level: advanced

 77:    Note:
 78:    Some solvers (such as smoothers in a `SNESFAS`) do not need the residual computed at the final solution so skip computing it
 79:    to save time.

 81: .seealso: `SNES`, `SNESSolve()`, `SNESGetAlwaysComputesFinalResidual()`
 82: @*/
 83: PetscErrorCode SNESSetAlwaysComputesFinalResidual(SNES snes, PetscBool flg)
 84: {
 86:   snes->alwayscomputesfinalresidual = flg;
 87:   return 0;
 88: }

 90: /*@
 91:     SNESGetAlwaysComputesFinalResidual - checks if the `SNES` always computes the residual at the final solution

 93:    Logically Collective on snes

 95:     Input Parameter:
 96: .   snes - the `SNES` context

 98:     Output Parameter:
 99: .   flg - `PETSC_TRUE` if the residual is computed

101:    Level: advanced

103: .seealso: `SNES`, `SNESSolve()`, `SNESSetAlwaysComputesFinalResidual()`
104: @*/
105: PetscErrorCode SNESGetAlwaysComputesFinalResidual(SNES snes, PetscBool *flg)
106: {
108:   *flg = snes->alwayscomputesfinalresidual;
109:   return 0;
110: }

112: /*@
113:    SNESSetFunctionDomainError - tells `SNES` that the input vector, a proposed new solution, to your function you provided to `SNESSetFunction()` is not
114:      in the functions domain. For example, a step with negative pressure.

116:    Logically Collective on snes

118:    Input Parameters:
119: .  snes - the `SNES` context

121:    Level: advanced

123:    Note:
124:    You can direct `SNES` to avoid certain steps by using `SNESVISetVariableBounds()`, `SNESVISetComputeVariableBounds()` or
125:    `SNESLineSearchSetPreCheck()`, `SNESLineSearchSetPostCheck()`

127: .seealso: `SNESCreate()`, `SNESSetFunction()`, `SNESFunction`, `SNESSetJacobianDomainError()`, `SNESVISetVariableBounds()`,
128:           `SNESVISetComputeVariableBounds()`, `SNESLineSearchSetPreCheck()`, `SNESLineSearchSetPostCheck()`
129: @*/
130: PetscErrorCode SNESSetFunctionDomainError(SNES snes)
131: {
134:   snes->domainerror = PETSC_TRUE;
135:   return 0;
136: }

138: /*@
139:    SNESSetJacobianDomainError - tells `SNES` that the function you provided to `SNESSetJacobian()` at the proposed step. For example there is a negative element transformation.

141:    Logically Collective on snes

143:    Input Parameters:
144: .  snes - the `SNES` context

146:    Level: advanced

148:    Note:
149:    You can direct `SNES` to avoid certain steps by using `SNESVISetVariableBounds()`, `SNESVISetComputeVariableBounds()` or
150:    `SNESLineSearchSetPreCheck()`, `SNESLineSearchSetPostCheck()`

152: .seealso: `SNESCreate()`, `SNESSetFunction()`, `SNESFunction()`, `SNESSetFunctionDomainError()`, `SNESVISetVariableBounds()`,
153:           `SNESVISetComputeVariableBounds()`, `SNESLineSearchSetPreCheck()`, `SNESLineSearchSetPostCheck()`
154: @*/
155: PetscErrorCode SNESSetJacobianDomainError(SNES snes)
156: {
159:   snes->jacobiandomainerror = PETSC_TRUE;
160:   return 0;
161: }

163: /*@
164:    SNESSetCheckJacobianDomainError - tells `SNESSolve()` whether to check if the user called `SNESSetJacobianDomainError()` Jacobian domain error after
165:    each Jacobian evaluation. By default, we check Jacobian domain error in the debug mode, and do not check it in the optimized mode.

167:    Logically Collective on snes

169:    Input Parameters:
170: +  snes - the SNES context
171: -  flg  - indicates if or not to check Jacobian domain error after each Jacobian evaluation

173:    Level: advanced

175:    Note:
176:    Checks require one extra parallel synchronization for each Jacobian evaluation

178: .seealso: `SNESCreate()`, `SNESSetFunction()`, `SNESFunction()`, `SNESSetFunctionDomainError()`, `SNESGetCheckJacobianDomainError()`
179: @*/
180: PetscErrorCode SNESSetCheckJacobianDomainError(SNES snes, PetscBool flg)
181: {
183:   snes->checkjacdomainerror = flg;
184:   return 0;
185: }

187: /*@
188:    SNESGetCheckJacobianDomainError - Get an indicator whether or not we are checking Jacobian domain errors after each Jacobian evaluation.

190:    Logically Collective on snes

192:    Input Parameters:
193: .  snes - the `SNES` context

195:    Output Parameters:
196: .  flg  - `PETSC_FALSE` indicates that we don't check Jacobian domain errors after each Jacobian evaluation

198:    Level: advanced

200: .seealso: `SNESCreate()`, `SNESSetFunction()`, `SNESFunction()`, `SNESSetFunctionDomainError()`, `SNESSetCheckJacobianDomainError()`
201: @*/
202: PetscErrorCode SNESGetCheckJacobianDomainError(SNES snes, PetscBool *flg)
203: {
206:   *flg = snes->checkjacdomainerror;
207:   return 0;
208: }

210: /*@
211:    SNESGetFunctionDomainError - Gets the status of the domain error after a call to `SNESComputeFunction()`;

213:    Logically Collective

215:    Input Parameters:
216: .  snes - the `SNES` context

218:    Output Parameters:
219: .  domainerror - Set to `PETSC_TRUE` if there's a domain error; `PETSC_FALSE` otherwise.

221:    Level: developer

223: .seealso: `SNESSetFunctionDomainError()`, `SNESComputeFunction()`
224: @*/
225: PetscErrorCode SNESGetFunctionDomainError(SNES snes, PetscBool *domainerror)
226: {
229:   *domainerror = snes->domainerror;
230:   return 0;
231: }

233: /*@
234:    SNESGetJacobianDomainError - Gets the status of the Jacobian domain error after a call to `SNESComputeJacobian()`;

236:    Logically Collective on snes

238:    Input Parameters:
239: .  snes - the `SNES` context

241:    Output Parameters:
242: .  domainerror - Set to `PETSC_TRUE` if there's a Jacobian domain error; `PETSC_FALSE` otherwise.

244:    Level: advanced

246: .seealso: `SNESSetFunctionDomainError()`, `SNESComputeFunction()`, `SNESGetFunctionDomainError()`
247: @*/
248: PetscErrorCode SNESGetJacobianDomainError(SNES snes, PetscBool *domainerror)
249: {
252:   *domainerror = snes->jacobiandomainerror;
253:   return 0;
254: }

256: /*@C
257:   SNESLoad - Loads a `SNES` that has been stored in `PETSCVIEWERBINARY` with `SNESView()`.

259:   Collective on snes

261:   Input Parameters:
262: + newdm - the newly loaded `SNES`, this needs to have been created with `SNESCreate()` or
263:            some related function before a call to `SNESLoad()`.
264: - viewer - binary file viewer, obtained from `PetscViewerBinaryOpen()`

266:    Level: intermediate

268:   Note:
269:    The type is determined by the data in the file, any type set into the `SNES` before this call is ignored.

271: .seealso: `SNESCreate()`, `SNESType`, `PetscViewerBinaryOpen()`, `SNESView()`, `MatLoad()`, `VecLoad()`
272: @*/
273: PetscErrorCode SNESLoad(SNES snes, PetscViewer viewer)
274: {
275:   PetscBool isbinary;
276:   PetscInt  classid;
277:   char      type[256];
278:   KSP       ksp;
279:   DM        dm;
280:   DMSNES    dmsnes;

284:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary);

287:   PetscViewerBinaryRead(viewer, &classid, 1, NULL, PETSC_INT);
289:   PetscViewerBinaryRead(viewer, type, 256, NULL, PETSC_CHAR);
290:   SNESSetType(snes, type);
291:   PetscTryTypeMethod(snes, load, viewer);
292:   SNESGetDM(snes, &dm);
293:   DMGetDMSNES(dm, &dmsnes);
294:   DMSNESLoad(dmsnes, viewer);
295:   SNESGetKSP(snes, &ksp);
296:   KSPLoad(ksp, viewer);
297:   return 0;
298: }

300: #include <petscdraw.h>
301: #if defined(PETSC_HAVE_SAWS)
302: #include <petscviewersaws.h>
303: #endif

305: /*@C
306:    SNESViewFromOptions - View a `SNES` based on the options database

308:    Collective on snes

310:    Input Parameters:
311: +  A - the `SNES` context
312: .  obj - Optional object
313: -  name - command line option

315:    Level: intermediate

317: .seealso: `SNES`, `SNESView`, `PetscObjectViewFromOptions()`, `SNESCreate()`
318: @*/
319: PetscErrorCode SNESViewFromOptions(SNES A, PetscObject obj, const char name[])
320: {
322:   PetscObjectViewFromOptions((PetscObject)A, obj, name);
323:   return 0;
324: }

326: PETSC_EXTERN PetscErrorCode SNESComputeJacobian_DMDA(SNES, Vec, Mat, Mat, void *);

328: /*@C
329:    SNESView - Prints the `SNES` data structure.

331:    Collective on snes

333:    Input Parameters:
334: +  snes - the `SNES` context
335: -  viewer - the `PetscViewer`

337:    Options Database Key:
338: .  -snes_view - Calls `SNESView()` at end of `SNESSolve()`

340:    Notes:
341:    The available visualization contexts include
342: +     `PETSC_VIEWER_STDOUT_SELF` - standard output (default)
343: -     `PETSC_VIEWER_STDOUT_WORLD` - synchronized standard
344:          output where only the first processor opens
345:          the file.  All other processors send their
346:          data to the first processor to print.

348:    The available formats include
349: +     `PETSC_VIEWER_DEFAULT` - standard output (default)
350: -     `PETSC_VIEWER_ASCII_INFO_DETAIL` - more verbose output for `SNESNASM`

352:    The user can open an alternative visualization context with
353:    `PetscViewerASCIIOpen()` - output to a specified file.

355:   In the debugger you can do "call `SNESView`(snes,0)" to display the `SNES` solver. (The same holds for any PETSc object viewer).

357:    Level: beginner

359: .seealso: `SNES`, `SNESLoad()`, `SNESCreate()`, `PetscViewerASCIIOpen()`
360: @*/
361: PetscErrorCode SNESView(SNES snes, PetscViewer viewer)
362: {
363:   SNESKSPEW     *kctx;
364:   KSP            ksp;
365:   SNESLineSearch linesearch;
366:   PetscBool      iascii, isstring, isbinary, isdraw;
367:   DMSNES         dmsnes;
368: #if defined(PETSC_HAVE_SAWS)
369:   PetscBool issaws;
370: #endif

373:   if (!viewer) PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)snes), &viewer);

377:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
378:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERSTRING, &isstring);
379:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary);
380:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw);
381: #if defined(PETSC_HAVE_SAWS)
382:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERSAWS, &issaws);
383: #endif
384:   if (iascii) {
385:     SNESNormSchedule normschedule;
386:     DM               dm;
387:     PetscErrorCode (*cJ)(SNES, Vec, Mat, Mat, void *);
388:     void       *ctx;
389:     const char *pre = "";

391:     PetscObjectPrintClassNamePrefixType((PetscObject)snes, viewer);
392:     if (!snes->setupcalled) PetscViewerASCIIPrintf(viewer, "  SNES has not been set up so information may be incomplete\n");
393:     if (snes->ops->view) {
394:       PetscViewerASCIIPushTab(viewer);
395:       PetscUseTypeMethod(snes, view, viewer);
396:       PetscViewerASCIIPopTab(viewer);
397:     }
398:     PetscViewerASCIIPrintf(viewer, "  maximum iterations=%" PetscInt_FMT ", maximum function evaluations=%" PetscInt_FMT "\n", snes->max_its, snes->max_funcs);
399:     PetscViewerASCIIPrintf(viewer, "  tolerances: relative=%g, absolute=%g, solution=%g\n", (double)snes->rtol, (double)snes->abstol, (double)snes->stol);
400:     if (snes->usesksp) PetscViewerASCIIPrintf(viewer, "  total number of linear solver iterations=%" PetscInt_FMT "\n", snes->linear_its);
401:     PetscViewerASCIIPrintf(viewer, "  total number of function evaluations=%" PetscInt_FMT "\n", snes->nfuncs);
402:     SNESGetNormSchedule(snes, &normschedule);
403:     if (normschedule > 0) PetscViewerASCIIPrintf(viewer, "  norm schedule %s\n", SNESNormSchedules[normschedule]);
404:     if (snes->gridsequence) PetscViewerASCIIPrintf(viewer, "  total number of grid sequence refinements=%" PetscInt_FMT "\n", snes->gridsequence);
405:     if (snes->ksp_ewconv) {
406:       kctx = (SNESKSPEW *)snes->kspconvctx;
407:       if (kctx) {
408:         PetscViewerASCIIPrintf(viewer, "  Eisenstat-Walker computation of KSP relative tolerance (version %" PetscInt_FMT ")\n", kctx->version);
409:         PetscViewerASCIIPrintf(viewer, "    rtol_0=%g, rtol_max=%g, threshold=%g\n", (double)kctx->rtol_0, (double)kctx->rtol_max, (double)kctx->threshold);
410:         PetscViewerASCIIPrintf(viewer, "    gamma=%g, alpha=%g, alpha2=%g\n", (double)kctx->gamma, (double)kctx->alpha, (double)kctx->alpha2);
411:       }
412:     }
413:     if (snes->lagpreconditioner == -1) {
414:       PetscViewerASCIIPrintf(viewer, "  Preconditioned is never rebuilt\n");
415:     } else if (snes->lagpreconditioner > 1) {
416:       PetscViewerASCIIPrintf(viewer, "  Preconditioned is rebuilt every %" PetscInt_FMT " new Jacobians\n", snes->lagpreconditioner);
417:     }
418:     if (snes->lagjacobian == -1) {
419:       PetscViewerASCIIPrintf(viewer, "  Jacobian is never rebuilt\n");
420:     } else if (snes->lagjacobian > 1) {
421:       PetscViewerASCIIPrintf(viewer, "  Jacobian is rebuilt every %" PetscInt_FMT " SNES iterations\n", snes->lagjacobian);
422:     }
423:     SNESGetDM(snes, &dm);
424:     DMSNESGetJacobian(dm, &cJ, &ctx);
425:     if (snes->mf_operator) {
426:       PetscViewerASCIIPrintf(viewer, "  Jacobian is applied matrix-free with differencing\n");
427:       pre = "Preconditioning ";
428:     }
429:     if (cJ == SNESComputeJacobianDefault) {
430:       PetscViewerASCIIPrintf(viewer, "  %sJacobian is built using finite differences one column at a time\n", pre);
431:     } else if (cJ == SNESComputeJacobianDefaultColor) {
432:       PetscViewerASCIIPrintf(viewer, "  %sJacobian is built using finite differences with coloring\n", pre);
433:       /* it slightly breaks data encapsulation for access the DMDA information directly */
434:     } else if (cJ == SNESComputeJacobian_DMDA) {
435:       MatFDColoring fdcoloring;
436:       PetscObjectQuery((PetscObject)dm, "DMDASNES_FDCOLORING", (PetscObject *)&fdcoloring);
437:       if (fdcoloring) {
438:         PetscViewerASCIIPrintf(viewer, "  %sJacobian is built using colored finite differences on a DMDA\n", pre);
439:       } else {
440:         PetscViewerASCIIPrintf(viewer, "  %sJacobian is built using a DMDA local Jacobian\n", pre);
441:       }
442:     } else if (snes->mf) {
443:       PetscViewerASCIIPrintf(viewer, "  Jacobian is applied matrix-free with differencing, no explicit Jacobian\n");
444:     }
445:   } else if (isstring) {
446:     const char *type;
447:     SNESGetType(snes, &type);
448:     PetscViewerStringSPrintf(viewer, " SNESType: %-7.7s", type);
449:     PetscTryTypeMethod(snes, view, viewer);
450:   } else if (isbinary) {
451:     PetscInt    classid = SNES_FILE_CLASSID;
452:     MPI_Comm    comm;
453:     PetscMPIInt rank;
454:     char        type[256];

456:     PetscObjectGetComm((PetscObject)snes, &comm);
457:     MPI_Comm_rank(comm, &rank);
458:     if (rank == 0) {
459:       PetscViewerBinaryWrite(viewer, &classid, 1, PETSC_INT);
460:       PetscStrncpy(type, ((PetscObject)snes)->type_name, sizeof(type));
461:       PetscViewerBinaryWrite(viewer, type, sizeof(type), PETSC_CHAR);
462:     }
463:     PetscTryTypeMethod(snes, view, viewer);
464:   } else if (isdraw) {
465:     PetscDraw draw;
466:     char      str[36];
467:     PetscReal x, y, bottom, h;

469:     PetscViewerDrawGetDraw(viewer, 0, &draw);
470:     PetscDrawGetCurrentPoint(draw, &x, &y);
471:     PetscStrncpy(str, "SNES: ", sizeof(str));
472:     PetscStrlcat(str, ((PetscObject)snes)->type_name, sizeof(str));
473:     PetscDrawStringBoxed(draw, x, y, PETSC_DRAW_BLUE, PETSC_DRAW_BLACK, str, NULL, &h);
474:     bottom = y - h;
475:     PetscDrawPushCurrentPoint(draw, x, bottom);
476:     PetscTryTypeMethod(snes, view, viewer);
477: #if defined(PETSC_HAVE_SAWS)
478:   } else if (issaws) {
479:     PetscMPIInt rank;
480:     const char *name;

482:     PetscObjectGetName((PetscObject)snes, &name);
483:     MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
484:     if (!((PetscObject)snes)->amsmem && rank == 0) {
485:       char dir[1024];

487:       PetscObjectViewSAWs((PetscObject)snes, viewer);
488:       PetscSNPrintf(dir, 1024, "/PETSc/Objects/%s/its", name);
489:       SAWs_Register, (dir, &snes->iter, 1, SAWs_READ, SAWs_INT);
490:       if (!snes->conv_hist) SNESSetConvergenceHistory(snes, NULL, NULL, PETSC_DECIDE, PETSC_TRUE);
491:       PetscSNPrintf(dir, 1024, "/PETSc/Objects/%s/conv_hist", name);
492:       SAWs_Register, (dir, snes->conv_hist, 10, SAWs_READ, SAWs_DOUBLE);
493:     }
494: #endif
495:   }
496:   if (snes->linesearch) {
497:     SNESGetLineSearch(snes, &linesearch);
498:     PetscViewerASCIIPushTab(viewer);
499:     SNESLineSearchView(linesearch, viewer);
500:     PetscViewerASCIIPopTab(viewer);
501:   }
502:   if (snes->npc && snes->usesnpc) {
503:     PetscViewerASCIIPushTab(viewer);
504:     SNESView(snes->npc, viewer);
505:     PetscViewerASCIIPopTab(viewer);
506:   }
507:   PetscViewerASCIIPushTab(viewer);
508:   DMGetDMSNES(snes->dm, &dmsnes);
509:   DMSNESView(dmsnes, viewer);
510:   PetscViewerASCIIPopTab(viewer);
511:   if (snes->usesksp) {
512:     SNESGetKSP(snes, &ksp);
513:     PetscViewerASCIIPushTab(viewer);
514:     KSPView(ksp, viewer);
515:     PetscViewerASCIIPopTab(viewer);
516:   }
517:   if (isdraw) {
518:     PetscDraw draw;
519:     PetscViewerDrawGetDraw(viewer, 0, &draw);
520:     PetscDrawPopCurrentPoint(draw);
521:   }
522:   return 0;
523: }

525: /*
526:   We retain a list of functions that also take SNES command
527:   line options. These are called at the end SNESSetFromOptions()
528: */
529: #define MAXSETFROMOPTIONS 5
530: static PetscInt numberofsetfromoptions;
531: static PetscErrorCode (*othersetfromoptions[MAXSETFROMOPTIONS])(SNES);

533: /*@C
534:   SNESAddOptionsChecker - Adds an additional function to check for `SNES` options.

536:   Not Collective

538:   Input Parameter:
539: . snescheck - function that checks for options

541:   Level: developer

543: .seealso: `SNES`, `SNESSetFromOptions()`
544: @*/
545: PetscErrorCode SNESAddOptionsChecker(PetscErrorCode (*snescheck)(SNES))
546: {
548:   othersetfromoptions[numberofsetfromoptions++] = snescheck;
549:   return 0;
550: }

552: static PetscErrorCode SNESSetUpMatrixFree_Private(SNES snes, PetscBool hasOperator, PetscInt version)
553: {
554:   Mat          J;
555:   MatNullSpace nullsp;


559:   if (!snes->vec_func && (snes->jacobian || snes->jacobian_pre)) {
560:     Mat A = snes->jacobian, B = snes->jacobian_pre;
561:     MatCreateVecs(A ? A : B, NULL, &snes->vec_func);
562:   }

564:   if (version == 1) {
565:     MatCreateSNESMF(snes, &J);
566:     MatMFFDSetOptionsPrefix(J, ((PetscObject)snes)->prefix);
567:     MatSetFromOptions(J);
568:     /* TODO: the version 2 code should be merged into the MatCreateSNESMF() and MatCreateMFFD() infrastructure and then removed */
569:   } else if (version == 2) {
571: #if !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_REAL_SINGLE) && !defined(PETSC_USE_REAL___FLOAT128) && !defined(PETSC_USE_REAL___FP16)
572:     MatCreateSNESMFMore(snes, snes->vec_func, &J);
573: #else
574:     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "matrix-free operator routines (version 2)");
575: #endif
576:   } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "matrix-free operator routines, only version 1 and 2");

578:   /* attach any user provided null space that was on Amat to the newly created matrix free matrix */
579:   if (snes->jacobian) {
580:     MatGetNullSpace(snes->jacobian, &nullsp);
581:     if (nullsp) MatSetNullSpace(J, nullsp);
582:   }

584:   PetscInfo(snes, "Setting default matrix-free operator routines (version %" PetscInt_FMT ")\n", version);
585:   if (hasOperator) {
586:     /* This version replaces the user provided Jacobian matrix with a
587:        matrix-free version but still employs the user-provided preconditioner matrix. */
588:     SNESSetJacobian(snes, J, NULL, NULL, NULL);
589:   } else {
590:     /* This version replaces both the user-provided Jacobian and the user-
591:      provided preconditioner Jacobian with the default matrix free version. */
592:     if (snes->npcside == PC_LEFT && snes->npc) {
593:       if (!snes->jacobian) SNESSetJacobian(snes, J, NULL, NULL, NULL);
594:     } else {
595:       KSP       ksp;
596:       PC        pc;
597:       PetscBool match;

599:       SNESSetJacobian(snes, J, J, MatMFFDComputeJacobian, NULL);
600:       /* Force no preconditioner */
601:       SNESGetKSP(snes, &ksp);
602:       KSPGetPC(ksp, &pc);
603:       PetscObjectTypeCompareAny((PetscObject)pc, &match, PCSHELL, PCH2OPUS, "");
604:       if (!match) {
605:         PetscInfo(snes, "Setting default matrix-free preconditioner routines\nThat is no preconditioner is being used\n");
606:         PCSetType(pc, PCNONE);
607:       }
608:     }
609:   }
610:   MatDestroy(&J);
611:   return 0;
612: }

614: static PetscErrorCode DMRestrictHook_SNESVecSol(DM dmfine, Mat Restrict, Vec Rscale, Mat Inject, DM dmcoarse, void *ctx)
615: {
616:   SNES snes = (SNES)ctx;
617:   Vec  Xfine, Xfine_named = NULL, Xcoarse;

619:   if (PetscLogPrintInfo) {
620:     PetscInt finelevel, coarselevel, fineclevel, coarseclevel;
621:     DMGetRefineLevel(dmfine, &finelevel);
622:     DMGetCoarsenLevel(dmfine, &fineclevel);
623:     DMGetRefineLevel(dmcoarse, &coarselevel);
624:     DMGetCoarsenLevel(dmcoarse, &coarseclevel);
625:     PetscInfo(dmfine, "Restricting SNES solution vector from level %" PetscInt_FMT "-%" PetscInt_FMT " to level %" PetscInt_FMT "-%" PetscInt_FMT "\n", finelevel, fineclevel, coarselevel, coarseclevel);
626:   }
627:   if (dmfine == snes->dm) Xfine = snes->vec_sol;
628:   else {
629:     DMGetNamedGlobalVector(dmfine, "SNESVecSol", &Xfine_named);
630:     Xfine = Xfine_named;
631:   }
632:   DMGetNamedGlobalVector(dmcoarse, "SNESVecSol", &Xcoarse);
633:   if (Inject) {
634:     MatRestrict(Inject, Xfine, Xcoarse);
635:   } else {
636:     MatRestrict(Restrict, Xfine, Xcoarse);
637:     VecPointwiseMult(Xcoarse, Xcoarse, Rscale);
638:   }
639:   DMRestoreNamedGlobalVector(dmcoarse, "SNESVecSol", &Xcoarse);
640:   if (Xfine_named) DMRestoreNamedGlobalVector(dmfine, "SNESVecSol", &Xfine_named);
641:   return 0;
642: }

644: static PetscErrorCode DMCoarsenHook_SNESVecSol(DM dm, DM dmc, void *ctx)
645: {
646:   DMCoarsenHookAdd(dmc, DMCoarsenHook_SNESVecSol, DMRestrictHook_SNESVecSol, ctx);
647:   return 0;
648: }

650: /* This may be called to rediscretize the operator on levels of linear multigrid. The DM shuffle is so the user can
651:  * safely call SNESGetDM() in their residual evaluation routine. */
652: static PetscErrorCode KSPComputeOperators_SNES(KSP ksp, Mat A, Mat B, void *ctx)
653: {
654:   SNES  snes = (SNES)ctx;
655:   Vec   X, Xnamed = NULL;
656:   DM    dmsave;
657:   void *ctxsave;
658:   PetscErrorCode (*jac)(SNES, Vec, Mat, Mat, void *) = NULL;

660:   dmsave = snes->dm;
661:   KSPGetDM(ksp, &snes->dm);
662:   if (dmsave == snes->dm) X = snes->vec_sol; /* We are on the finest level */
663:   else { /* We are on a coarser level, this vec was initialized using a DM restrict hook */ DMGetNamedGlobalVector(snes->dm, "SNESVecSol", &Xnamed);
664:     X = Xnamed;
665:     SNESGetJacobian(snes, NULL, NULL, &jac, &ctxsave);
666:     /* If the DM's don't match up, the MatFDColoring context needed for the jacobian won't match up either -- fixit. */
667:     if (jac == SNESComputeJacobianDefaultColor) SNESSetJacobian(snes, NULL, NULL, SNESComputeJacobianDefaultColor, NULL);
668:   }
669:   /* Make sure KSP DM has the Jacobian computation routine */
670:   {
671:     DMSNES sdm;

673:     DMGetDMSNES(snes->dm, &sdm);
674:     if (!sdm->ops->computejacobian) DMCopyDMSNES(dmsave, snes->dm);
675:   }
676:   /* Compute the operators */
677:   SNESComputeJacobian(snes, X, A, B);
678:   /* Put the previous context back */
679:   if (snes->dm != dmsave && jac == SNESComputeJacobianDefaultColor) SNESSetJacobian(snes, NULL, NULL, jac, ctxsave);

681:   if (Xnamed) DMRestoreNamedGlobalVector(snes->dm, "SNESVecSol", &Xnamed);
682:   snes->dm = dmsave;
683:   return 0;
684: }

686: /*@
687:    SNESSetUpMatrices - ensures that matrices are available for `SNES`, this is called by `SNESSetUp_XXX()`

689:    Collective

691:    Input Parameter:
692: .  snes - snes to configure

694:    Level: developer

696: .seealso: `SNES`, `SNESSetUp()`
697: @*/
698: PetscErrorCode SNESSetUpMatrices(SNES snes)
699: {
700:   DM     dm;
701:   DMSNES sdm;

703:   SNESGetDM(snes, &dm);
704:   DMGetDMSNES(dm, &sdm);
705:   if (!snes->jacobian && snes->mf) {
706:     Mat   J;
707:     void *functx;
708:     MatCreateSNESMF(snes, &J);
709:     MatMFFDSetOptionsPrefix(J, ((PetscObject)snes)->prefix);
710:     MatSetFromOptions(J);
711:     SNESGetFunction(snes, NULL, NULL, &functx);
712:     SNESSetJacobian(snes, J, J, NULL, NULL);
713:     MatDestroy(&J);
714:   } else if (snes->mf_operator && !snes->jacobian_pre && !snes->jacobian) {
715:     Mat J, B;
716:     MatCreateSNESMF(snes, &J);
717:     MatMFFDSetOptionsPrefix(J, ((PetscObject)snes)->prefix);
718:     MatSetFromOptions(J);
719:     DMCreateMatrix(snes->dm, &B);
720:     /* sdm->computejacobian was already set to reach here */
721:     SNESSetJacobian(snes, J, B, NULL, NULL);
722:     MatDestroy(&J);
723:     MatDestroy(&B);
724:   } else if (!snes->jacobian_pre) {
725:     PetscDS   prob;
726:     Mat       J, B;
727:     PetscBool hasPrec = PETSC_FALSE;

729:     J = snes->jacobian;
730:     DMGetDS(dm, &prob);
731:     if (prob) PetscDSHasJacobianPreconditioner(prob, &hasPrec);
732:     if (J) PetscObjectReference((PetscObject)J);
733:     else if (hasPrec) DMCreateMatrix(snes->dm, &J);
734:     DMCreateMatrix(snes->dm, &B);
735:     SNESSetJacobian(snes, J ? J : B, B, NULL, NULL);
736:     MatDestroy(&J);
737:     MatDestroy(&B);
738:   }
739:   {
740:     KSP ksp;
741:     SNESGetKSP(snes, &ksp);
742:     KSPSetComputeOperators(ksp, KSPComputeOperators_SNES, snes);
743:     DMCoarsenHookAdd(snes->dm, DMCoarsenHook_SNESVecSol, DMRestrictHook_SNESVecSol, snes);
744:   }
745:   return 0;
746: }

748: static PetscErrorCode SNESMonitorPauseFinal_Internal(SNES snes)
749: {
750:   PetscInt i;

752:   if (!snes->pauseFinal) return 0;
753:   for (i = 0; i < snes->numbermonitors; ++i) {
754:     PetscViewerAndFormat *vf = (PetscViewerAndFormat *)snes->monitorcontext[i];
755:     PetscDraw             draw;
756:     PetscReal             lpause;

758:     if (!vf) continue;
759:     if (vf->lg) {
761:       if (((PetscObject)vf->lg)->classid != PETSC_DRAWLG_CLASSID) continue;
762:       PetscDrawLGGetDraw(vf->lg, &draw);
763:       PetscDrawGetPause(draw, &lpause);
764:       PetscDrawSetPause(draw, -1.0);
765:       PetscDrawPause(draw);
766:       PetscDrawSetPause(draw, lpause);
767:     } else {
768:       PetscBool isdraw;

771:       if (((PetscObject)vf->viewer)->classid != PETSC_VIEWER_CLASSID) continue;
772:       PetscObjectTypeCompare((PetscObject)vf->viewer, PETSCVIEWERDRAW, &isdraw);
773:       if (!isdraw) continue;
774:       PetscViewerDrawGetDraw(vf->viewer, 0, &draw);
775:       PetscDrawGetPause(draw, &lpause);
776:       PetscDrawSetPause(draw, -1.0);
777:       PetscDrawPause(draw);
778:       PetscDrawSetPause(draw, lpause);
779:     }
780:   }
781:   return 0;
782: }

784: /*@C
785:    SNESMonitorSetFromOptions - Sets a monitor function and viewer appropriate for the type indicated by the user

787:    Collective on snes

789:    Input Parameters:
790: +  snes - SNES object you wish to monitor
791: .  name - the monitor type one is seeking
792: .  help - message indicating what monitoring is done
793: .  manual - manual page for the monitor
794: .  monitor - the monitor function
795: -  monitorsetup - a function that is called once ONLY if the user selected this monitor that may set additional features of the `SNES` or `PetscViewer` objects

797:    Options Database Key:
798: .  -name - trigger the use of this monitor in `SNESSetFromOptions()`

800:    Level: advanced

802: .seealso: `PetscOptionsGetViewer()`, `PetscOptionsGetReal()`, `PetscOptionsHasName()`, `PetscOptionsGetString()`,
803:           `PetscOptionsGetIntArray()`, `PetscOptionsGetRealArray()`, `PetscOptionsBool()`
804:           `PetscOptionsInt()`, `PetscOptionsString()`, `PetscOptionsReal()`, `PetscOptionsBool()`,
805:           `PetscOptionsName()`, `PetscOptionsBegin()`, `PetscOptionsEnd()`, `PetscOptionsHeadBegin()`,
806:           `PetscOptionsStringArray()`, `PetscOptionsRealArray()`, `PetscOptionsScalar()`,
807:           `PetscOptionsBoolGroupBegin()`, `PetscOptionsBoolGroup()`, `PetscOptionsBoolGroupEnd()`,
808:           `PetscOptionsFList()`, `PetscOptionsEList()`
809: @*/
810: PetscErrorCode SNESMonitorSetFromOptions(SNES snes, const char name[], const char help[], const char manual[], PetscErrorCode (*monitor)(SNES, PetscInt, PetscReal, PetscViewerAndFormat *), PetscErrorCode (*monitorsetup)(SNES, PetscViewerAndFormat *))
811: {
812:   PetscViewer       viewer;
813:   PetscViewerFormat format;
814:   PetscBool         flg;

816:   PetscOptionsGetViewer(PetscObjectComm((PetscObject)snes), ((PetscObject)snes)->options, ((PetscObject)snes)->prefix, name, &viewer, &format, &flg);
817:   if (flg) {
818:     PetscViewerAndFormat *vf;
819:     PetscViewerAndFormatCreate(viewer, format, &vf);
820:     PetscObjectDereference((PetscObject)viewer);
821:     if (monitorsetup) (*monitorsetup)(snes, vf);
822:     SNESMonitorSet(snes, (PetscErrorCode(*)(SNES, PetscInt, PetscReal, void *))monitor, vf, (PetscErrorCode(*)(void **))PetscViewerAndFormatDestroy);
823:   }
824:   return 0;
825: }

827: PetscErrorCode SNESEWSetFromOptions_Private(SNESKSPEW *kctx, MPI_Comm comm, const char *prefix)
828: {
829:   PetscOptionsBegin(comm, prefix, "Eisenstat and Walker type forcing options", "KSP");
830:   PetscOptionsInt("-ksp_ew_version", "Version 1, 2 or 3", NULL, kctx->version, &kctx->version, NULL);
831:   PetscOptionsReal("-ksp_ew_rtol0", "0 <= rtol0 < 1", NULL, kctx->rtol_0, &kctx->rtol_0, NULL);
832:   PetscOptionsReal("-ksp_ew_rtolmax", "0 <= rtolmax < 1", NULL, kctx->rtol_max, &kctx->rtol_max, NULL);
833:   PetscOptionsReal("-ksp_ew_gamma", "0 <= gamma <= 1", NULL, kctx->gamma, &kctx->gamma, NULL);
834:   PetscOptionsReal("-ksp_ew_alpha", "1 < alpha <= 2", NULL, kctx->alpha, &kctx->alpha, NULL);
835:   PetscOptionsReal("-ksp_ew_alpha2", "alpha2", NULL, kctx->alpha2, &kctx->alpha2, NULL);
836:   PetscOptionsReal("-ksp_ew_threshold", "0 < threshold < 1", NULL, kctx->threshold, &kctx->threshold, NULL);
837:   PetscOptionsReal("-ksp_ew_v4_p1", "p1", NULL, kctx->v4_p1, &kctx->v4_p1, NULL);
838:   PetscOptionsReal("-ksp_ew_v4_p2", "p2", NULL, kctx->v4_p2, &kctx->v4_p2, NULL);
839:   PetscOptionsReal("-ksp_ew_v4_p3", "p3", NULL, kctx->v4_p3, &kctx->v4_p3, NULL);
840:   PetscOptionsReal("-ksp_ew_v4_m1", "Scaling when rk-1 in [p2,p3)", NULL, kctx->v4_m1, &kctx->v4_m1, NULL);
841:   PetscOptionsReal("-ksp_ew_v4_m2", "Scaling when rk-1 in [p3,+infty)", NULL, kctx->v4_m2, &kctx->v4_m2, NULL);
842:   PetscOptionsReal("-ksp_ew_v4_m3", "Threshold for successive rtol (0.1 in Eq.7)", NULL, kctx->v4_m3, &kctx->v4_m3, NULL);
843:   PetscOptionsReal("-ksp_ew_v4_m4", "Adaptation scaling (0.5 in Eq.7)", NULL, kctx->v4_m4, &kctx->v4_m4, NULL);
844:   PetscOptionsEnd();
845:   return 0;
846: }

848: /*@
849:    SNESSetFromOptions - Sets various `SNES` and `KSP` parameters from user options.

851:    Collective on snes

853:    Input Parameter:
854: .  snes - the `SNES` context

856:    Options Database Keys:
857: +  -snes_type <type> - newtonls, newtontr, ngmres, ncg, nrichardson, qn, vi, fas, `SNESType` for complete list
858: .  -snes_stol - convergence tolerance in terms of the norm
859:                 of the change in the solution between steps
860: .  -snes_atol <abstol> - absolute tolerance of residual norm
861: .  -snes_rtol <rtol> - relative decrease in tolerance norm from initial
862: .  -snes_divergence_tolerance <divtol> - if the residual goes above divtol*rnorm0, exit with divergence
863: .  -snes_force_iteration <force> - force SNESSolve() to take at least one iteration
864: .  -snes_max_it <max_it> - maximum number of iterations
865: .  -snes_max_funcs <max_funcs> - maximum number of function evaluations
866: .  -snes_max_fail <max_fail> - maximum number of line search failures allowed before stopping, default is none
867: .  -snes_max_linear_solve_fail - number of linear solver failures before SNESSolve() stops
868: .  -snes_lag_preconditioner <lag> - how often preconditioner is rebuilt (use -1 to never rebuild)
869: .  -snes_lag_preconditioner_persists <true,false> - retains the -snes_lag_preconditioner information across multiple SNESSolve()
870: .  -snes_lag_jacobian <lag> - how often Jacobian is rebuilt (use -1 to never rebuild)
871: .  -snes_lag_jacobian_persists <true,false> - retains the -snes_lag_jacobian information across multiple SNESSolve()
872: .  -snes_trtol <trtol> - trust region tolerance
873: .  -snes_convergence_test - <default,skip,correct_pressure> convergence test in nonlinear solver.
874:                                default `SNESConvergedDefault()`. skip `SNESConvergedSkip()` means continue iterating until max_it or some other criterion is reached, saving expense
875:                                of convergence test. correct_pressure S`NESConvergedCorrectPressure()` has special handling of a pressure null space.
876: .  -snes_monitor [ascii][:filename][:viewer format] - prints residual norm at each iteration. if no filename given prints to stdout
877: .  -snes_monitor_solution [ascii binary draw][:filename][:viewer format] - plots solution at each iteration
878: .  -snes_monitor_residual [ascii binary draw][:filename][:viewer format] - plots residual (not its norm) at each iteration
879: .  -snes_monitor_solution_update [ascii binary draw][:filename][:viewer format] - plots update to solution at each iteration
880: .  -snes_monitor_lg_residualnorm - plots residual norm at each iteration
881: .  -snes_monitor_lg_range - plots residual norm at each iteration
882: .  -snes_monitor_pause_final - Pauses all monitor drawing after the solver ends
883: .  -snes_fd - use finite differences to compute Jacobian; very slow, only for testing
884: .  -snes_fd_color - use finite differences with coloring to compute Jacobian
885: .  -snes_mf_ksp_monitor - if using matrix-free multiply then print h at each KSP iteration
886: .  -snes_converged_reason - print the reason for convergence/divergence after each solve
887: .  -npc_snes_type <type> - the SNES type to use as a nonlinear preconditioner
888: .   -snes_test_jacobian <optional threshold> - compare the user provided Jacobian with one computed via finite differences to check for errors.  If a threshold is given, display only those entries whose difference is greater than the threshold.
889: -   -snes_test_jacobian_view - display the user provided Jacobian, the finite difference Jacobian and the difference between them to help users detect the location of errors in the user provided Jacobian.

891:     Options Database Keys for Eisenstat-Walker method:
892: +  -snes_ksp_ew - use Eisenstat-Walker method for determining linear system convergence
893: .  -snes_ksp_ew_version ver - version of  Eisenstat-Walker method
894: .  -snes_ksp_ew_rtol0 <rtol0> - Sets rtol0
895: .  -snes_ksp_ew_rtolmax <rtolmax> - Sets rtolmax
896: .  -snes_ksp_ew_gamma <gamma> - Sets gamma
897: .  -snes_ksp_ew_alpha <alpha> - Sets alpha
898: .  -snes_ksp_ew_alpha2 <alpha2> - Sets alpha2
899: -  -snes_ksp_ew_threshold <threshold> - Sets threshold

901:    Notes:
902:    To see all options, run your program with the -help option or consult the users manual

904:    `SNES` supports three approaches for computing (approximate) Jacobians: user provided via `SNESSetJacobian()`, matrix free, and computing explicitly with
905:    finite differences and coloring using `MatFDColoring`. It is also possible to use automatic differentiation and the `MatFDColoring` object.

907:    Level: beginner

909: .seealso: `SNESType`, `SNESSetOptionsPrefix()`, `SNESResetFromOptions()`, `SNES`, `SNESCreate()`
910: @*/
911: PetscErrorCode SNESSetFromOptions(SNES snes)
912: {
913:   PetscBool   flg, pcset, persist, set;
914:   PetscInt    i, indx, lag, grids;
915:   const char *deft        = SNESNEWTONLS;
916:   const char *convtests[] = {"default", "skip", "correct_pressure"};
917:   SNESKSPEW  *kctx        = NULL;
918:   char        type[256], monfilename[PETSC_MAX_PATH_LEN], ewprefix[256];
919:   PCSide      pcside;
920:   const char *optionsprefix;

923:   SNESRegisterAll();
924:   PetscObjectOptionsBegin((PetscObject)snes);
925:   if (((PetscObject)snes)->type_name) deft = ((PetscObject)snes)->type_name;
926:   PetscOptionsFList("-snes_type", "Nonlinear solver method", "SNESSetType", SNESList, deft, type, 256, &flg);
927:   if (flg) {
928:     SNESSetType(snes, type);
929:   } else if (!((PetscObject)snes)->type_name) {
930:     SNESSetType(snes, deft);
931:   }
932:   PetscOptionsReal("-snes_stol", "Stop if step length less than", "SNESSetTolerances", snes->stol, &snes->stol, NULL);
933:   PetscOptionsReal("-snes_atol", "Stop if function norm less than", "SNESSetTolerances", snes->abstol, &snes->abstol, NULL);

935:   PetscOptionsReal("-snes_rtol", "Stop if decrease in function norm less than", "SNESSetTolerances", snes->rtol, &snes->rtol, NULL);
936:   PetscOptionsReal("-snes_divergence_tolerance", "Stop if residual norm increases by this factor", "SNESSetDivergenceTolerance", snes->divtol, &snes->divtol, NULL);
937:   PetscOptionsInt("-snes_max_it", "Maximum iterations", "SNESSetTolerances", snes->max_its, &snes->max_its, NULL);
938:   PetscOptionsInt("-snes_max_funcs", "Maximum function evaluations", "SNESSetTolerances", snes->max_funcs, &snes->max_funcs, NULL);
939:   PetscOptionsInt("-snes_max_fail", "Maximum nonlinear step failures", "SNESSetMaxNonlinearStepFailures", snes->maxFailures, &snes->maxFailures, NULL);
940:   PetscOptionsInt("-snes_max_linear_solve_fail", "Maximum failures in linear solves allowed", "SNESSetMaxLinearSolveFailures", snes->maxLinearSolveFailures, &snes->maxLinearSolveFailures, NULL);
941:   PetscOptionsBool("-snes_error_if_not_converged", "Generate error if solver does not converge", "SNESSetErrorIfNotConverged", snes->errorifnotconverged, &snes->errorifnotconverged, NULL);
942:   PetscOptionsBool("-snes_force_iteration", "Force SNESSolve() to take at least one iteration", "SNESSetForceIteration", snes->forceiteration, &snes->forceiteration, NULL);
943:   PetscOptionsBool("-snes_check_jacobian_domain_error", "Check Jacobian domain error after Jacobian evaluation", "SNESCheckJacobianDomainError", snes->checkjacdomainerror, &snes->checkjacdomainerror, NULL);

945:   PetscOptionsInt("-snes_lag_preconditioner", "How often to rebuild preconditioner", "SNESSetLagPreconditioner", snes->lagpreconditioner, &lag, &flg);
946:   if (flg) {
948:     SNESSetLagPreconditioner(snes, lag);
949:   }
950:   PetscOptionsBool("-snes_lag_preconditioner_persists", "Preconditioner lagging through multiple SNES solves", "SNESSetLagPreconditionerPersists", snes->lagjac_persist, &persist, &flg);
951:   if (flg) SNESSetLagPreconditionerPersists(snes, persist);
952:   PetscOptionsInt("-snes_lag_jacobian", "How often to rebuild Jacobian", "SNESSetLagJacobian", snes->lagjacobian, &lag, &flg);
953:   if (flg) {
955:     SNESSetLagJacobian(snes, lag);
956:   }
957:   PetscOptionsBool("-snes_lag_jacobian_persists", "Jacobian lagging through multiple SNES solves", "SNESSetLagJacobianPersists", snes->lagjac_persist, &persist, &flg);
958:   if (flg) SNESSetLagJacobianPersists(snes, persist);

960:   PetscOptionsInt("-snes_grid_sequence", "Use grid sequencing to generate initial guess", "SNESSetGridSequence", snes->gridsequence, &grids, &flg);
961:   if (flg) SNESSetGridSequence(snes, grids);

963:   PetscOptionsEList("-snes_convergence_test", "Convergence test", "SNESSetConvergenceTest", convtests, sizeof(convtests) / sizeof(char *), "default", &indx, &flg);
964:   if (flg) {
965:     switch (indx) {
966:     case 0:
967:       SNESSetConvergenceTest(snes, SNESConvergedDefault, NULL, NULL);
968:       break;
969:     case 1:
970:       SNESSetConvergenceTest(snes, SNESConvergedSkip, NULL, NULL);
971:       break;
972:     case 2:
973:       SNESSetConvergenceTest(snes, SNESConvergedCorrectPressure, NULL, NULL);
974:       break;
975:     }
976:   }

978:   PetscOptionsEList("-snes_norm_schedule", "SNES Norm schedule", "SNESSetNormSchedule", SNESNormSchedules, 5, "function", &indx, &flg);
979:   if (flg) SNESSetNormSchedule(snes, (SNESNormSchedule)indx);

981:   PetscOptionsEList("-snes_function_type", "SNES Norm schedule", "SNESSetFunctionType", SNESFunctionTypes, 2, "unpreconditioned", &indx, &flg);
982:   if (flg) SNESSetFunctionType(snes, (SNESFunctionType)indx);

984:   kctx = (SNESKSPEW *)snes->kspconvctx;

986:   PetscOptionsBool("-snes_ksp_ew", "Use Eisentat-Walker linear system convergence test", "SNESKSPSetUseEW", snes->ksp_ewconv, &snes->ksp_ewconv, NULL);

988:   SNESGetOptionsPrefix(snes, &optionsprefix);
989:   PetscSNPrintf(ewprefix, sizeof(ewprefix), "%s%s", optionsprefix ? optionsprefix : "", "snes_");
990:   SNESEWSetFromOptions_Private(kctx, PetscObjectComm((PetscObject)snes), ewprefix);

992:   PetscOptionsInt("-snes_ksp_ew_version", "Version 1, 2 or 3", "SNESKSPSetParametersEW", kctx->version, &kctx->version, NULL);
993:   PetscOptionsReal("-snes_ksp_ew_rtol0", "0 <= rtol0 < 1", "SNESKSPSetParametersEW", kctx->rtol_0, &kctx->rtol_0, NULL);
994:   PetscOptionsReal("-snes_ksp_ew_rtolmax", "0 <= rtolmax < 1", "SNESKSPSetParametersEW", kctx->rtol_max, &kctx->rtol_max, NULL);
995:   PetscOptionsReal("-snes_ksp_ew_gamma", "0 <= gamma <= 1", "SNESKSPSetParametersEW", kctx->gamma, &kctx->gamma, NULL);
996:   PetscOptionsReal("-snes_ksp_ew_alpha", "1 < alpha <= 2", "SNESKSPSetParametersEW", kctx->alpha, &kctx->alpha, NULL);
997:   PetscOptionsReal("-snes_ksp_ew_alpha2", "alpha2", "SNESKSPSetParametersEW", kctx->alpha2, &kctx->alpha2, NULL);
998:   PetscOptionsReal("-snes_ksp_ew_threshold", "0 < threshold < 1", "SNESKSPSetParametersEW", kctx->threshold, &kctx->threshold, NULL);

1000:   flg = PETSC_FALSE;
1001:   PetscOptionsBool("-snes_monitor_cancel", "Remove all monitors", "SNESMonitorCancel", flg, &flg, &set);
1002:   if (set && flg) SNESMonitorCancel(snes);

1004:   SNESMonitorSetFromOptions(snes, "-snes_monitor", "Monitor norm of function", "SNESMonitorDefault", SNESMonitorDefault, SNESMonitorDefaultSetUp);
1005:   SNESMonitorSetFromOptions(snes, "-snes_monitor_short", "Monitor norm of function with fewer digits", "SNESMonitorDefaultShort", SNESMonitorDefaultShort, NULL);
1006:   SNESMonitorSetFromOptions(snes, "-snes_monitor_range", "Monitor range of elements of function", "SNESMonitorRange", SNESMonitorRange, NULL);

1008:   SNESMonitorSetFromOptions(snes, "-snes_monitor_ratio", "Monitor ratios of the norm of function for consecutive steps", "SNESMonitorRatio", SNESMonitorRatio, SNESMonitorRatioSetUp);
1009:   SNESMonitorSetFromOptions(snes, "-snes_monitor_field", "Monitor norm of function (split into fields)", "SNESMonitorDefaultField", SNESMonitorDefaultField, NULL);
1010:   SNESMonitorSetFromOptions(snes, "-snes_monitor_solution", "View solution at each iteration", "SNESMonitorSolution", SNESMonitorSolution, NULL);
1011:   SNESMonitorSetFromOptions(snes, "-snes_monitor_solution_update", "View correction at each iteration", "SNESMonitorSolutionUpdate", SNESMonitorSolutionUpdate, NULL);
1012:   SNESMonitorSetFromOptions(snes, "-snes_monitor_residual", "View residual at each iteration", "SNESMonitorResidual", SNESMonitorResidual, NULL);
1013:   SNESMonitorSetFromOptions(snes, "-snes_monitor_jacupdate_spectrum", "Print the change in the spectrum of the Jacobian", "SNESMonitorJacUpdateSpectrum", SNESMonitorJacUpdateSpectrum, NULL);
1014:   SNESMonitorSetFromOptions(snes, "-snes_monitor_fields", "Monitor norm of function per field", "SNESMonitorSet", SNESMonitorFields, NULL);
1015:   PetscOptionsBool("-snes_monitor_pause_final", "Pauses all draw monitors at the final iterate", "SNESMonitorPauseFinal_Internal", PETSC_FALSE, &snes->pauseFinal, NULL);

1017:   PetscOptionsString("-snes_monitor_python", "Use Python function", "SNESMonitorSet", NULL, monfilename, sizeof(monfilename), &flg);
1018:   if (flg) PetscPythonMonitorSet((PetscObject)snes, monfilename);

1020:   flg = PETSC_FALSE;
1021:   PetscOptionsBool("-snes_monitor_lg_range", "Plot function range at each iteration", "SNESMonitorLGRange", flg, &flg, NULL);
1022:   if (flg) {
1023:     PetscViewer ctx;

1025:     PetscViewerDrawOpen(PetscObjectComm((PetscObject)snes), NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 400, 300, &ctx);
1026:     SNESMonitorSet(snes, SNESMonitorLGRange, ctx, (PetscErrorCode(*)(void **))PetscViewerDestroy);
1027:   }

1029:   flg = PETSC_FALSE;
1030:   PetscOptionsBool("-snes_converged_reason_view_cancel", "Remove all converged reason viewers", "SNESConvergedReasonViewCancel", flg, &flg, &set);
1031:   if (set && flg) SNESConvergedReasonViewCancel(snes);

1033:   flg = PETSC_FALSE;
1034:   PetscOptionsBool("-snes_fd", "Use finite differences (slow) to compute Jacobian", "SNESComputeJacobianDefault", flg, &flg, NULL);
1035:   if (flg) {
1036:     void *functx;
1037:     DM    dm;
1038:     SNESGetDM(snes, &dm);
1039:     DMSNESUnsetJacobianContext_Internal(dm);
1040:     SNESGetFunction(snes, NULL, NULL, &functx);
1041:     SNESSetJacobian(snes, snes->jacobian, snes->jacobian_pre, SNESComputeJacobianDefault, functx);
1042:     PetscInfo(snes, "Setting default finite difference Jacobian matrix\n");
1043:   }

1045:   flg = PETSC_FALSE;
1046:   PetscOptionsBool("-snes_fd_function", "Use finite differences (slow) to compute function from user objective", "SNESObjectiveComputeFunctionDefaultFD", flg, &flg, NULL);
1047:   if (flg) SNESSetFunction(snes, NULL, SNESObjectiveComputeFunctionDefaultFD, NULL);

1049:   flg = PETSC_FALSE;
1050:   PetscOptionsBool("-snes_fd_color", "Use finite differences with coloring to compute Jacobian", "SNESComputeJacobianDefaultColor", flg, &flg, NULL);
1051:   if (flg) {
1052:     DM dm;
1053:     SNESGetDM(snes, &dm);
1054:     DMSNESUnsetJacobianContext_Internal(dm);
1055:     SNESSetJacobian(snes, snes->jacobian, snes->jacobian_pre, SNESComputeJacobianDefaultColor, NULL);
1056:     PetscInfo(snes, "Setting default finite difference coloring Jacobian matrix\n");
1057:   }

1059:   flg = PETSC_FALSE;
1060:   PetscOptionsBool("-snes_mf_operator", "Use a Matrix-Free Jacobian with user-provided preconditioner matrix", "SNESSetUseMatrixFree", PETSC_FALSE, &snes->mf_operator, &flg);
1061:   if (flg && snes->mf_operator) {
1062:     snes->mf_operator = PETSC_TRUE;
1063:     snes->mf          = PETSC_TRUE;
1064:   }
1065:   flg = PETSC_FALSE;
1066:   PetscOptionsBool("-snes_mf", "Use a Matrix-Free Jacobian with no preconditioner matrix", "SNESSetUseMatrixFree", PETSC_FALSE, &snes->mf, &flg);
1067:   if (!flg && snes->mf_operator) snes->mf = PETSC_TRUE;
1068:   PetscOptionsInt("-snes_mf_version", "Matrix-Free routines version 1 or 2", "None", snes->mf_version, &snes->mf_version, NULL);

1070:   flg = PETSC_FALSE;
1071:   SNESGetNPCSide(snes, &pcside);
1072:   PetscOptionsEnum("-snes_npc_side", "SNES nonlinear preconditioner side", "SNESSetNPCSide", PCSides, (PetscEnum)pcside, (PetscEnum *)&pcside, &flg);
1073:   if (flg) SNESSetNPCSide(snes, pcside);

1075: #if defined(PETSC_HAVE_SAWS)
1076:   /*
1077:     Publish convergence information using SAWs
1078:   */
1079:   flg = PETSC_FALSE;
1080:   PetscOptionsBool("-snes_monitor_saws", "Publish SNES progress using SAWs", "SNESMonitorSet", flg, &flg, NULL);
1081:   if (flg) {
1082:     void *ctx;
1083:     SNESMonitorSAWsCreate(snes, &ctx);
1084:     SNESMonitorSet(snes, SNESMonitorSAWs, ctx, SNESMonitorSAWsDestroy);
1085:   }
1086: #endif
1087: #if defined(PETSC_HAVE_SAWS)
1088:   {
1089:     PetscBool set;
1090:     flg = PETSC_FALSE;
1091:     PetscOptionsBool("-snes_saws_block", "Block for SAWs at end of SNESSolve", "PetscObjectSAWsBlock", ((PetscObject)snes)->amspublishblock, &flg, &set);
1092:     if (set) PetscObjectSAWsSetBlock((PetscObject)snes, flg);
1093:   }
1094: #endif

1096:   for (i = 0; i < numberofsetfromoptions; i++) (*othersetfromoptions[i])(snes);

1098:   PetscTryTypeMethod(snes, setfromoptions, PetscOptionsObject);

1100:   /* process any options handlers added with PetscObjectAddOptionsHandler() */
1101:   PetscObjectProcessOptionsHandlers((PetscObject)snes, PetscOptionsObject);
1102:   PetscOptionsEnd();

1104:   if (snes->linesearch) {
1105:     SNESGetLineSearch(snes, &snes->linesearch);
1106:     SNESLineSearchSetFromOptions(snes->linesearch);
1107:   }

1109:   if (snes->usesksp) {
1110:     if (!snes->ksp) SNESGetKSP(snes, &snes->ksp);
1111:     KSPSetOperators(snes->ksp, snes->jacobian, snes->jacobian_pre);
1112:     KSPSetFromOptions(snes->ksp);
1113:   }

1115:   /* if user has set the SNES NPC type via options database, create it. */
1116:   SNESGetOptionsPrefix(snes, &optionsprefix);
1117:   PetscOptionsHasName(((PetscObject)snes)->options, optionsprefix, "-npc_snes_type", &pcset);
1118:   if (pcset && (!snes->npc)) SNESGetNPC(snes, &snes->npc);
1119:   if (snes->npc) SNESSetFromOptions(snes->npc);
1120:   snes->setfromoptionscalled++;
1121:   return 0;
1122: }

1124: /*@
1125:    SNESResetFromOptions - Sets various `SNES` and `KSP` parameters from user options ONLY if the `SNESSetFromOptions()` was previously set from options

1127:    Collective on snes

1129:    Input Parameter:
1130: .  snes - the `SNES` context

1132:    Level: beginner

1134: .seealso: `SNES`, `SNESSetFromOptions()`, `SNESSetOptionsPrefix()`
1135: @*/
1136: PetscErrorCode SNESResetFromOptions(SNES snes)
1137: {
1138:   if (snes->setfromoptionscalled) SNESSetFromOptions(snes);
1139:   return 0;
1140: }

1142: /*@C
1143:    SNESSetComputeApplicationContext - Sets an optional function to compute a user-defined context for
1144:    the nonlinear solvers.

1146:    Logically Collective on snes

1148:    Input Parameters:
1149: +  snes - the `SNES` context
1150: .  compute - function to compute the context
1151: -  destroy - function to destroy the context

1153:    Level: intermediate

1155:    Note:
1156:    This routine is useful if you are performing grid sequencing or using `SNESFAS` and need the appropriate context generated for each level.

1158:    Use `SNESSetApplicationContext()` to see the context immediately

1160:    Fortran Note:
1161:    This function is currently not available from Fortran.

1163: .seealso: `SNESGetApplicationContext()`, `SNESSetComputeApplicationContext()`, `SNESSetApplicationContext()`
1164: @*/
1165: PetscErrorCode SNESSetComputeApplicationContext(SNES snes, PetscErrorCode (*compute)(SNES, void **), PetscErrorCode (*destroy)(void **))
1166: {
1168:   snes->ops->usercompute = compute;
1169:   snes->ops->userdestroy = destroy;
1170:   return 0;
1171: }

1173: /*@
1174:    SNESSetApplicationContext - Sets the optional user-defined context for the nonlinear solvers.

1176:    Logically Collective on snes

1178:    Input Parameters:
1179: +  snes - the `SNES` context
1180: -  usrP - optional user context

1182:    Level: intermediate

1184:    Notes:
1185:    Users can provide a context when constructing the `SNES` options and then access it inside their function, Jacobian, or other evaluation function
1186:    with `SNESGetApplicationContext()`

1188:    To provide a function that computes the context for you use `SNESSetComputeApplicationContext()`

1190:    Fortran Note:
1191:     To use this from Fortran you must write a Fortran interface definition for this
1192:     function that tells Fortran the Fortran derived data type that you are passing in as the ctx argument.

1194: .seealso: `SNES`, `SNESSetComputeApplicationContext()`, `SNESGetApplicationContext()`
1195: @*/
1196: PetscErrorCode SNESSetApplicationContext(SNES snes, void *usrP)
1197: {
1198:   KSP ksp;

1201:   SNESGetKSP(snes, &ksp);
1202:   KSPSetApplicationContext(ksp, usrP);
1203:   snes->user = usrP;
1204:   return 0;
1205: }

1207: /*@
1208:    SNESGetApplicationContext - Gets the user-defined context for the
1209:    nonlinear solvers set with `SNESGetApplicationContext()` or with `SNESSetComputeApplicationContext()`

1211:    Not Collective

1213:    Input Parameter:
1214: .  snes - `SNES` context

1216:    Output Parameter:
1217: .  usrP - user context

1219:    Fortran Note:
1220:     To use this from Fortran you must write a Fortran interface definition for this
1221:     function that tells Fortran the Fortran derived data type that you are passing in as the ctx argument.

1223:    Level: intermediate

1225: .seealso: `SNESSetApplicationContext()`
1226: @*/
1227: PetscErrorCode SNESGetApplicationContext(SNES snes, void *usrP)
1228: {
1230:   *(void **)usrP = snes->user;
1231:   return 0;
1232: }

1234: /*@
1235:    SNESSetUseMatrixFree - indicates that `SNES` should use matrix free finite difference matrix vector products to apply the Jacobian.

1237:    Logically Collective on snes, the values must be the same on all MPI ranks

1239:    Input Parameters:
1240: +  snes - `SNES` context
1241: .  mf_operator - use matrix-free only for the Amat used by `SNESSetJacobian()`, this means the user provided Pmat will continue to be used
1242: -  mf - use matrix-free for both the Amat and Pmat used by `SNESSetJacobian()`, both the Amat and Pmat set in `SNESSetJacobian()` will be ignored. With
1243:    this option no matrix element based preconditioners can be used in the linear solve since the matrix won't be explicitly available

1245:    Options Database Keys:
1246: + -snes_mf_operator - use matrix free only for the mat operator
1247: . -snes_mf - use matrix-free for both the mat and pmat operator
1248: . -snes_fd_color - compute the Jacobian via coloring and finite differences.
1249: - -snes_fd - compute the Jacobian via finite differences (slow)

1251:    Level: intermediate

1253:    Note:
1254:    SNES supports three approaches for computing (approximate) Jacobians: user provided via `SNESSetJacobian()`, matrix-free, and computing explicitly with
1255:    finite differences and coloring using `MatFDColoring`. It is also possible to use automatic differentiation and the `MatFDColoring` object.

1257: .seealso: `SNES`, `SNESGetUseMatrixFree()`, `MatCreateSNESMF()`, `SNESComputeJacobianDefaultColor()`
1258: @*/
1259: PetscErrorCode SNESSetUseMatrixFree(SNES snes, PetscBool mf_operator, PetscBool mf)
1260: {
1264:   snes->mf          = mf_operator ? PETSC_TRUE : mf;
1265:   snes->mf_operator = mf_operator;
1266:   return 0;
1267: }

1269: /*@
1270:    SNESGetUseMatrixFree - indicates if the SNES uses matrix-free finite difference matrix vector products to apply the Jacobian.

1272:    Not Collective, but the resulting flags will be the same on all MPI ranks

1274:    Input Parameter:
1275: .  snes - `SNES` context

1277:    Output Parameters:
1278: +  mf_operator - use matrix-free only for the Amat used by `SNESSetJacobian()`, this means the user provided Pmat will continue to be used
1279: -  mf - use matrix-free for both the Amat and Pmat used by `SNESSetJacobian()`, both the Amat and Pmat set in `SNESSetJacobian()` will be ignored

1281:    Level: intermediate

1283: .seealso: `SNES`, `SNESSetUseMatrixFree()`, `MatCreateSNESMF()`
1284: @*/
1285: PetscErrorCode SNESGetUseMatrixFree(SNES snes, PetscBool *mf_operator, PetscBool *mf)
1286: {
1288:   if (mf) *mf = snes->mf;
1289:   if (mf_operator) *mf_operator = snes->mf_operator;
1290:   return 0;
1291: }

1293: /*@
1294:    SNESGetIterationNumber - Gets the number of nonlinear iterations completed
1295:    at this time.

1297:    Not Collective

1299:    Input Parameter:
1300: .  snes - `SNES` context

1302:    Output Parameter:
1303: .  iter - iteration number

1305:    Notes:
1306:    For example, during the computation of iteration 2 this would return 1.

1308:    This is useful for using lagged Jacobians (where one does not recompute the
1309:    Jacobian at each `SNES` iteration). For example, the code
1310: .vb
1311:       SNESGetIterationNumber(snes,&it);
1312:       if (!(it % 2)) {
1313:         [compute Jacobian here]
1314:       }
1315: .ve
1316:    can be used in your function that computes the Jacobian to cause the Jacobian to be
1317:    recomputed every second `SNES` iteration. See also `SNESSetLagJacobian()`

1319:    After the `SNES` solve is complete this will return the number of nonlinear iterations used.

1321:    Level: intermediate

1323: .seealso: `SNES`, `SNESSolve()`, `SNESSetLagJacobian()`, `SNESGetLinearSolveIterations()`
1324: @*/
1325: PetscErrorCode SNESGetIterationNumber(SNES snes, PetscInt *iter)
1326: {
1329:   *iter = snes->iter;
1330:   return 0;
1331: }

1333: /*@
1334:    SNESSetIterationNumber - Sets the current iteration number.

1336:    Not Collective

1338:    Input Parameters:
1339: +  snes - `SNES` context
1340: -  iter - iteration number

1342:    Level: developer

1344: .seealso: `SNESGetLinearSolveIterations()`
1345: @*/
1346: PetscErrorCode SNESSetIterationNumber(SNES snes, PetscInt iter)
1347: {
1349:   PetscObjectSAWsTakeAccess((PetscObject)snes);
1350:   snes->iter = iter;
1351:   PetscObjectSAWsGrantAccess((PetscObject)snes);
1352:   return 0;
1353: }

1355: /*@
1356:    SNESGetNonlinearStepFailures - Gets the number of unsuccessful steps
1357:    attempted by the nonlinear solver.

1359:    Not Collective

1361:    Input Parameter:
1362: .  snes - `SNES` context

1364:    Output Parameter:
1365: .  nfails - number of unsuccessful steps attempted

1367:    Note:
1368:    This counter is reset to zero for each successive call to `SNESSolve()`.

1370:    Level: intermediate

1372: .seealso: `SNES`, `SNESGetMaxLinearSolveFailures()`, `SNESGetLinearSolveIterations()`, `SNESSetMaxLinearSolveFailures()`, `SNESGetLinearSolveFailures()`,
1373:           `SNESSetMaxNonlinearStepFailures()`, `SNESGetMaxNonlinearStepFailures()`
1374: @*/
1375: PetscErrorCode SNESGetNonlinearStepFailures(SNES snes, PetscInt *nfails)
1376: {
1379:   *nfails = snes->numFailures;
1380:   return 0;
1381: }

1383: /*@
1384:    SNESSetMaxNonlinearStepFailures - Sets the maximum number of unsuccessful steps
1385:    attempted by the nonlinear solver before it gives up and generates an error

1387:    Not Collective

1389:    Input Parameters:
1390: +  snes     - `SNES` context
1391: -  maxFails - maximum of unsuccessful steps

1393:    Level: intermediate

1395: .seealso: `SNESSetErrorIfNotConverged()`, `SNESGetMaxLinearSolveFailures()`, `SNESGetLinearSolveIterations()`, `SNESSetMaxLinearSolveFailures()`, `SNESGetLinearSolveFailures()`,
1396:           `SNESGetMaxNonlinearStepFailures()`, `SNESGetNonlinearStepFailures()`
1397: @*/
1398: PetscErrorCode SNESSetMaxNonlinearStepFailures(SNES snes, PetscInt maxFails)
1399: {
1401:   snes->maxFailures = maxFails;
1402:   return 0;
1403: }

1405: /*@
1406:    SNESGetMaxNonlinearStepFailures - Gets the maximum number of unsuccessful steps
1407:    attempted by the nonlinear solver before it gives up and generates an error

1409:    Not Collective

1411:    Input Parameter:
1412: .  snes     - SNES context

1414:    Output Parameter:
1415: .  maxFails - maximum of unsuccessful steps

1417:    Level: intermediate

1419: .seealso: `SNESSetErrorIfNotConverged()`, `SNESGetMaxLinearSolveFailures()`, `SNESGetLinearSolveIterations()`, `SNESSetMaxLinearSolveFailures()`, `SNESGetLinearSolveFailures()`,
1420:           `SNESSetMaxNonlinearStepFailures()`, `SNESGetNonlinearStepFailures()`
1421: @*/
1422: PetscErrorCode SNESGetMaxNonlinearStepFailures(SNES snes, PetscInt *maxFails)
1423: {
1426:   *maxFails = snes->maxFailures;
1427:   return 0;
1428: }

1430: /*@
1431:    SNESGetNumberFunctionEvals - Gets the number of user provided function evaluations
1432:      done by the `SNES` object

1434:    Not Collective

1436:    Input Parameter:
1437: .  snes     - `SNES` context

1439:    Output Parameter:
1440: .  nfuncs - number of evaluations

1442:    Level: intermediate

1444:    Note:
1445:     Reset every time `SNESSolve()` is called unless `SNESSetCountersReset()` is used.

1447: .seealso: `SNES`, `SNESGetMaxLinearSolveFailures()`, `SNESGetLinearSolveIterations()`, `SNESSetMaxLinearSolveFailures()`, `SNESGetLinearSolveFailures()`, `SNESSetCountersReset()`
1448: @*/
1449: PetscErrorCode SNESGetNumberFunctionEvals(SNES snes, PetscInt *nfuncs)
1450: {
1453:   *nfuncs = snes->nfuncs;
1454:   return 0;
1455: }

1457: /*@
1458:    SNESGetLinearSolveFailures - Gets the number of failed (non-converged)
1459:    linear solvers.

1461:    Not Collective

1463:    Input Parameter:
1464: .  snes - `SNES` context

1466:    Output Parameter:
1467: .  nfails - number of failed solves

1469:    Options Database Key:
1470: . -snes_max_linear_solve_fail <num> - The number of failures before the solve is terminated

1472:    Level: intermediate

1474:    Note:
1475:    This counter is reset to zero for each successive call to `SNESSolve()`.

1477: .seealso: `SNESGetMaxLinearSolveFailures()`, `SNESGetLinearSolveIterations()`, `SNESSetMaxLinearSolveFailures()`
1478: @*/
1479: PetscErrorCode SNESGetLinearSolveFailures(SNES snes, PetscInt *nfails)
1480: {
1483:   *nfails = snes->numLinearSolveFailures;
1484:   return 0;
1485: }

1487: /*@
1488:    SNESSetMaxLinearSolveFailures - the number of failed linear solve attempts
1489:    allowed before `SNES` returns with a diverged reason of `SNES_DIVERGED_LINEAR_SOLVE`

1491:    Logically Collective on snes

1493:    Input Parameters:
1494: +  snes     - `SNES` context
1495: -  maxFails - maximum allowed linear solve failures

1497:    Level: intermediate

1499:    Options Database Key:
1500: . -snes_max_linear_solve_fail <num> - The number of failures before the solve is terminated

1502:    Note:
1503:     By default this is 0; that is `SNES` returns on the first failed linear solve

1505: .seealso: `SNESSetErrorIfNotConverged()`, `SNESGetLinearSolveFailures()`, `SNESGetMaxLinearSolveFailures()`, `SNESGetLinearSolveIterations()`
1506: @*/
1507: PetscErrorCode SNESSetMaxLinearSolveFailures(SNES snes, PetscInt maxFails)
1508: {
1511:   snes->maxLinearSolveFailures = maxFails;
1512:   return 0;
1513: }

1515: /*@
1516:    SNESGetMaxLinearSolveFailures - gets the maximum number of linear solve failures that
1517:      are allowed before `SNES` returns as unsuccessful

1519:    Not Collective

1521:    Input Parameter:
1522: .  snes     - `SNES` context

1524:    Output Parameter:
1525: .  maxFails - maximum of unsuccessful solves allowed

1527:    Level: intermediate

1529:    Note:
1530:     By default this is 1; that is `SNES` returns on the first failed linear solve

1532: .seealso: `SNESSetErrorIfNotConverged()`, `SNESGetLinearSolveFailures()`, `SNESGetLinearSolveIterations()`, `SNESSetMaxLinearSolveFailures()`,
1533: @*/
1534: PetscErrorCode SNESGetMaxLinearSolveFailures(SNES snes, PetscInt *maxFails)
1535: {
1538:   *maxFails = snes->maxLinearSolveFailures;
1539:   return 0;
1540: }

1542: /*@
1543:    SNESGetLinearSolveIterations - Gets the total number of linear iterations
1544:    used by the nonlinear solver.

1546:    Not Collective

1548:    Input Parameter:
1549: .  snes - `SNES` context

1551:    Output Parameter:
1552: .  lits - number of linear iterations

1554:    Notes:
1555:    This counter is reset to zero for each successive call to `SNESSolve()` unless `SNESSetCountersReset()` is used.

1557:    If the linear solver fails inside the `SNESSolve()` the iterations for that call to the linear solver are not included. If you wish to count them
1558:    then call `KSPGetIterationNumber()` after the failed solve.

1560:    Level: intermediate

1562: .seealso: `SNES`, `SNESGetIterationNumber()`, `SNESGetLinearSolveFailures()`, `SNESGetMaxLinearSolveFailures()`, `SNESSetCountersReset()`
1563: @*/
1564: PetscErrorCode SNESGetLinearSolveIterations(SNES snes, PetscInt *lits)
1565: {
1568:   *lits = snes->linear_its;
1569:   return 0;
1570: }

1572: /*@
1573:    SNESSetCountersReset - Sets whether or not the counters for linear iterations and function evaluations
1574:    are reset every time `SNESSolve()` is called.

1576:    Logically Collective on snes

1578:    Input Parameters:
1579: +  snes - `SNES` context
1580: -  reset - whether to reset the counters or not, defaults to `PETSC_TRUE`

1582:    Level: developer

1584: .seealso: `SNESGetNumberFunctionEvals()`, `SNESGetLinearSolveIterations()`, `SNESGetNPC()`
1585: @*/
1586: PetscErrorCode SNESSetCountersReset(SNES snes, PetscBool reset)
1587: {
1590:   snes->counters_reset = reset;
1591:   return 0;
1592: }

1594: /*@
1595:    SNESSetKSP - Sets a `KSP` context for the `SNES` object to use

1597:    Not Collective, but the `SNES` and `KSP` objects must live on the same MPI_Comm

1599:    Input Parameters:
1600: +  snes - the `SNES` context
1601: -  ksp - the `KSP` context

1603:    Notes:
1604:    The `SNES` object already has its `KSP` object, you can obtain with `SNESGetKSP()`
1605:    so this routine is rarely needed.

1607:    The `KSP` object that is already in the `SNES` object has its reference count
1608:    decreased by one.

1610:    Level: developer

1612: .seealso: `SNES`, `KSP`, `KSPGetPC()`, `SNESCreate()`, `KSPCreate()`, `SNESSetKSP()`
1613: @*/
1614: PetscErrorCode SNESSetKSP(SNES snes, KSP ksp)
1615: {
1619:   PetscObjectReference((PetscObject)ksp);
1620:   if (snes->ksp) PetscObjectDereference((PetscObject)snes->ksp);
1621:   snes->ksp = ksp;
1622:   return 0;
1623: }

1625: /*@
1626:    SNESCreate - Creates a nonlinear solver context.

1628:    Collective

1630:    Input Parameter:
1631: .  comm - MPI communicator

1633:    Output Parameter:
1634: .  outsnes - the new SNES context

1636:    Options Database Keys:
1637: +   -snes_mf - Activates default matrix-free Jacobian-vector products,
1638:                and no preconditioning matrix
1639: .   -snes_mf_operator - Activates default matrix-free Jacobian-vector
1640:                products, and a user-provided preconditioning matrix
1641:                as set by SNESSetJacobian()
1642: -   -snes_fd - Uses (slow!) finite differences to compute Jacobian

1644:    Level: beginner

1646:    Developer Notes:
1647:    `SNES` always creates a `KSP` object even though many `SNES` methods do not use it. This is
1648:    unfortunate and should be fixed at some point. The flag snes->usesksp indicates if the
1649:    particular method does use `KSP` and regulates if the information about the `KSP` is printed
1650:    in `SNESView()`.

1652:    `TSSetFromOptions()` does call `SNESSetFromOptions()` which can lead to users being confused
1653:    by help messages about meaningless `SNES` options.

1655:    `SNES` always creates the snes->kspconvctx even though it is used by only one type. This should
1656:                     be fixed.

1658: .seealso: `SNES`, `SNESSolve()`, `SNESDestroy()`, `SNES`, `SNESSetLagPreconditioner()`, `SNESSetLagJacobian()`
1659: @*/
1660: PetscErrorCode SNESCreate(MPI_Comm comm, SNES *outsnes)
1661: {
1662:   SNES       snes;
1663:   SNESKSPEW *kctx;

1666:   *outsnes = NULL;
1667:   SNESInitializePackage();

1669:   PetscHeaderCreate(snes, SNES_CLASSID, "SNES", "Nonlinear solver", "SNES", comm, SNESDestroy, SNESView);

1671:   snes->ops->converged = SNESConvergedDefault;
1672:   snes->usesksp        = PETSC_TRUE;
1673:   snes->tolerancesset  = PETSC_FALSE;
1674:   snes->max_its        = 50;
1675:   snes->max_funcs      = 10000;
1676:   snes->norm           = 0.0;
1677:   snes->xnorm          = 0.0;
1678:   snes->ynorm          = 0.0;
1679:   snes->normschedule   = SNES_NORM_ALWAYS;
1680:   snes->functype       = SNES_FUNCTION_DEFAULT;
1681: #if defined(PETSC_USE_REAL_SINGLE)
1682:   snes->rtol = 1.e-5;
1683: #else
1684:   snes->rtol = 1.e-8;
1685: #endif
1686:   snes->ttol = 0.0;
1687: #if defined(PETSC_USE_REAL_SINGLE)
1688:   snes->abstol = 1.e-25;
1689: #else
1690:   snes->abstol = 1.e-50;
1691: #endif
1692: #if defined(PETSC_USE_REAL_SINGLE)
1693:   snes->stol = 1.e-5;
1694: #else
1695:   snes->stol = 1.e-8;
1696: #endif
1697: #if defined(PETSC_USE_REAL_SINGLE)
1698:   snes->deltatol = 1.e-6;
1699: #else
1700:   snes->deltatol = 1.e-12;
1701: #endif
1702:   snes->divtol               = 1.e4;
1703:   snes->rnorm0               = 0;
1704:   snes->nfuncs               = 0;
1705:   snes->numFailures          = 0;
1706:   snes->maxFailures          = 1;
1707:   snes->linear_its           = 0;
1708:   snes->lagjacobian          = 1;
1709:   snes->jac_iter             = 0;
1710:   snes->lagjac_persist       = PETSC_FALSE;
1711:   snes->lagpreconditioner    = 1;
1712:   snes->pre_iter             = 0;
1713:   snes->lagpre_persist       = PETSC_FALSE;
1714:   snes->numbermonitors       = 0;
1715:   snes->numberreasonviews    = 0;
1716:   snes->data                 = NULL;
1717:   snes->setupcalled          = PETSC_FALSE;
1718:   snes->ksp_ewconv           = PETSC_FALSE;
1719:   snes->nwork                = 0;
1720:   snes->work                 = NULL;
1721:   snes->nvwork               = 0;
1722:   snes->vwork                = NULL;
1723:   snes->conv_hist_len        = 0;
1724:   snes->conv_hist_max        = 0;
1725:   snes->conv_hist            = NULL;
1726:   snes->conv_hist_its        = NULL;
1727:   snes->conv_hist_reset      = PETSC_TRUE;
1728:   snes->counters_reset       = PETSC_TRUE;
1729:   snes->vec_func_init_set    = PETSC_FALSE;
1730:   snes->reason               = SNES_CONVERGED_ITERATING;
1731:   snes->npcside              = PC_RIGHT;
1732:   snes->setfromoptionscalled = 0;

1734:   snes->mf          = PETSC_FALSE;
1735:   snes->mf_operator = PETSC_FALSE;
1736:   snes->mf_version  = 1;

1738:   snes->numLinearSolveFailures = 0;
1739:   snes->maxLinearSolveFailures = 1;

1741:   snes->vizerotolerance     = 1.e-8;
1742:   snes->checkjacdomainerror = PetscDefined(USE_DEBUG) ? PETSC_TRUE : PETSC_FALSE;

1744:   /* Set this to true if the implementation of SNESSolve_XXX does compute the residual at the final solution. */
1745:   snes->alwayscomputesfinalresidual = PETSC_FALSE;

1747:   /* Create context to compute Eisenstat-Walker relative tolerance for KSP */
1748:   PetscNew(&kctx);

1750:   snes->kspconvctx  = (void *)kctx;
1751:   kctx->version     = 2;
1752:   kctx->rtol_0      = 0.3; /* Eisenstat and Walker suggest rtol_0=.5, but
1753:                              this was too large for some test cases */
1754:   kctx->rtol_last   = 0.0;
1755:   kctx->rtol_max    = 0.9;
1756:   kctx->gamma       = 1.0;
1757:   kctx->alpha       = 0.5 * (1.0 + PetscSqrtReal(5.0));
1758:   kctx->alpha2      = kctx->alpha;
1759:   kctx->threshold   = 0.1;
1760:   kctx->lresid_last = 0.0;
1761:   kctx->norm_last   = 0.0;

1763:   kctx->rk_last     = 0.0;
1764:   kctx->rk_last_2   = 0.0;
1765:   kctx->rtol_last_2 = 0.0;
1766:   kctx->v4_p1       = 0.1;
1767:   kctx->v4_p2       = 0.4;
1768:   kctx->v4_p3       = 0.7;
1769:   kctx->v4_m1       = 0.8;
1770:   kctx->v4_m2       = 0.5;
1771:   kctx->v4_m3       = 0.1;
1772:   kctx->v4_m4       = 0.5;

1774:   *outsnes = snes;
1775:   return 0;
1776: }

1778: /*MC
1779:     SNESFunction - Functional form used to convey the nonlinear function to `SNES` in `SNESSetFunction()`

1781:      Synopsis:
1782:      #include "petscsnes.h"
1783:      PetscErrorCode SNESFunction(SNES snes,Vec x,Vec f,void *ctx);

1785:      Collective on snes

1787:      Input Parameters:
1788: +     snes - the `SNES` context
1789: .     x    - state at which to evaluate residual
1790: -     ctx     - optional user-defined function context, passed in with `SNESSetFunction()`

1792:      Output Parameter:
1793: .     f  - vector to put residual (function value)

1795:    Level: intermediate

1797: .seealso: `SNESSetFunction()`, `SNESGetFunction()`
1798: M*/

1800: /*@C
1801:    SNESSetFunction - Sets the function evaluation routine and function
1802:    vector for use by the `SNES` routines in solving systems of nonlinear
1803:    equations.

1805:    Logically Collective on snes

1807:    Input Parameters:
1808: +  snes - the `SNES` context
1809: .  r - vector to store function values, may be NULL
1810: .  f - function evaluation routine; see `SNESFunction` for calling sequence details
1811: -  ctx - [optional] user-defined context for private data for the
1812:          function evaluation routine (may be NULL)

1814:    Level: beginner

1816: .seealso: `SNES`, `SNESGetFunction()`, `SNESComputeFunction()`, `SNESSetJacobian()`, `SNESSetPicard()`, `SNESFunction`
1817: @*/
1818: PetscErrorCode SNESSetFunction(SNES snes, Vec r, PetscErrorCode (*f)(SNES, Vec, Vec, void *), void *ctx)
1819: {
1820:   DM dm;

1823:   if (r) {
1826:     PetscObjectReference((PetscObject)r);
1827:     VecDestroy(&snes->vec_func);
1828:     snes->vec_func = r;
1829:   }
1830:   SNESGetDM(snes, &dm);
1831:   DMSNESSetFunction(dm, f, ctx);
1832:   if (f == SNESPicardComputeFunction) DMSNESSetMFFunction(dm, SNESPicardComputeMFFunction, ctx);
1833:   return 0;
1834: }

1836: /*@C
1837:    SNESSetInitialFunction - Sets the function vector to be used as the
1838:    initial function value at the initialization of the method.  In some
1839:    instances, the user has precomputed the function before calling
1840:    `SNESSolve()`.  This function allows one to avoid a redundant call
1841:    to `SNESComputeFunction()` in that case.

1843:    Logically Collective on snes

1845:    Input Parameters:
1846: +  snes - the `SNES` context
1847: -  f - vector to store function value

1849:    Notes:
1850:    This should not be modified during the solution procedure.

1852:    This is used extensively in the `SNESFAS` hierarchy and in nonlinear preconditioning.

1854:    Level: developer

1856: .seealso: `SNES`, `SNESFAS`, `SNESSetFunction()`, `SNESComputeFunction()`, `SNESSetInitialFunctionNorm()`
1857: @*/
1858: PetscErrorCode SNESSetInitialFunction(SNES snes, Vec f)
1859: {
1860:   Vec vec_func;

1865:   if (snes->npcside == PC_LEFT && snes->functype == SNES_FUNCTION_PRECONDITIONED) {
1866:     snes->vec_func_init_set = PETSC_FALSE;
1867:     return 0;
1868:   }
1869:   SNESGetFunction(snes, &vec_func, NULL, NULL);
1870:   VecCopy(f, vec_func);

1872:   snes->vec_func_init_set = PETSC_TRUE;
1873:   return 0;
1874: }

1876: /*@
1877:    SNESSetNormSchedule - Sets the `SNESNormSchedule` used in convergence and monitoring
1878:    of the `SNES` method, when norms are computed in the solving process

1880:    Logically Collective on snes

1882:    Input Parameters:
1883: +  snes - the `SNES` context
1884: -  normschedule - the frequency of norm computation

1886:    Options Database Key:
1887: .  -snes_norm_schedule <none, always, initialonly, finalonly, initialfinalonly> - set the schedule

1889:    Notes:
1890:    Only certain `SNES` methods support certain `SNESNormSchedules`.  Most require evaluation
1891:    of the nonlinear function and the taking of its norm at every iteration to
1892:    even ensure convergence at all.  However, methods such as custom Gauss-Seidel methods
1893:    `SNESNGS` and the like do not require the norm of the function to be computed, and therefore
1894:    may either be monitored for convergence or not.  As these are often used as nonlinear
1895:    preconditioners, monitoring the norm of their error is not a useful enterprise within
1896:    their solution.

1898:    Level: advanced

1900: .seealso: `SNESNormSchedule`, `SNESGetNormSchedule()`, `SNESComputeFunction()`, `VecNorm()`, `SNESSetFunction()`, `SNESSetInitialFunction()`, `SNESNormSchedule`
1901: @*/
1902: PetscErrorCode SNESSetNormSchedule(SNES snes, SNESNormSchedule normschedule)
1903: {
1905:   snes->normschedule = normschedule;
1906:   return 0;
1907: }

1909: /*@
1910:    SNESGetNormSchedule - Gets the `SNESNormSchedule` used in convergence and monitoring
1911:    of the `SNES` method.

1913:    Logically Collective on snes

1915:    Input Parameters:
1916: +  snes - the `SNES` context
1917: -  normschedule - the type of the norm used

1919:    Level: advanced

1921: .seealso: `SNES`, `SNESSetNormSchedule()`, `SNESComputeFunction()`, `VecNorm()`, `SNESSetFunction()`, `SNESSetInitialFunction()`, `SNESNormSchedule`
1922: @*/
1923: PetscErrorCode SNESGetNormSchedule(SNES snes, SNESNormSchedule *normschedule)
1924: {
1926:   *normschedule = snes->normschedule;
1927:   return 0;
1928: }

1930: /*@
1931:   SNESSetFunctionNorm - Sets the last computed residual norm.

1933:   Logically Collective on snes

1935:   Input Parameters:
1936: +  snes - the `SNES` context
1937: -  norm - the value of the norm

1939:   Level: developer

1941: .seealso: `SNES`, `SNESGetNormSchedule()`, `SNESComputeFunction()`, `VecNorm()`, `SNESSetFunction()`, `SNESSetInitialFunction()`, `SNESNormSchedule`
1942: @*/
1943: PetscErrorCode SNESSetFunctionNorm(SNES snes, PetscReal norm)
1944: {
1946:   snes->norm = norm;
1947:   return 0;
1948: }

1950: /*@
1951:   SNESGetFunctionNorm - Gets the last computed norm of the residual

1953:   Not Collective

1955:   Input Parameter:
1956: . snes - the `SNES` context

1958:   Output Parameter:
1959: . norm - the last computed residual norm

1961:   Level: developer

1963: .seealso: `SNES`, `SNESSetNormSchedule()`, `SNESComputeFunction()`, `VecNorm()`, `SNESSetFunction()`, `SNESSetInitialFunction()`, `SNESNormSchedule`
1964: @*/
1965: PetscErrorCode SNESGetFunctionNorm(SNES snes, PetscReal *norm)
1966: {
1969:   *norm = snes->norm;
1970:   return 0;
1971: }

1973: /*@
1974:   SNESGetUpdateNorm - Gets the last computed norm of the solution update

1976:   Not Collective

1978:   Input Parameter:
1979: . snes - the `SNES` context

1981:   Output Parameter:
1982: . ynorm - the last computed update norm

1984:   Level: developer

1986:   Note:
1987:   The new solution is the current solution plus the update, so this norm is an indication of the size of the update

1989: .seealso: `SNES`, `SNESSetNormSchedule()`, `SNESComputeFunction()`, `SNESGetFunctionNorm()`
1990: @*/
1991: PetscErrorCode SNESGetUpdateNorm(SNES snes, PetscReal *ynorm)
1992: {
1995:   *ynorm = snes->ynorm;
1996:   return 0;
1997: }

1999: /*@
2000:   SNESGetSolutionNorm - Gets the last computed norm of the solution

2002:   Not Collective

2004:   Input Parameter:
2005: . snes - the `SNES` context

2007:   Output Parameter:
2008: . xnorm - the last computed solution norm

2010:   Level: developer

2012: .seealso: `SNES`, `SNESSetNormSchedule()`, `SNESComputeFunction()`, `SNESGetFunctionNorm()`, `SNESGetUpdateNorm()`
2013: @*/
2014: PetscErrorCode SNESGetSolutionNorm(SNES snes, PetscReal *xnorm)
2015: {
2018:   *xnorm = snes->xnorm;
2019:   return 0;
2020: }

2022: /*@C
2023:    SNESSetFunctionType - Sets the `SNESFunctionType`
2024:    of the `SNES` method.

2026:    Logically Collective on snes

2028:    Input Parameters:
2029: +  snes - the `SNES` context
2030: -  type - the function type

2032:    Level: developer

2034:    Notes:
2035:    Possible values of the function type
2036: +  `SNES_FUNCTION_DEFAULT` - the default for the given `SNESType`
2037: .  `SNES_FUNCTION_UNPRECONDITIONED` - an unpreconditioned function evaluation (this is the function provided with `SNESSetFunction()`
2038: -  `SNES_FUNCTION_PRECONDITIONED` - a transformation of the function provided with `SNESSetFunction()`

2040:    Different `SNESType`s use this value in different ways

2042: .seealso: `SNES`, `SNESFunctionType`, `SNESGetNormSchedule()`, `SNESComputeFunction()`, `VecNorm()`, `SNESSetFunction()`, `SNESSetInitialFunction()`, `SNESNormSchedule`
2043: @*/
2044: PetscErrorCode SNESSetFunctionType(SNES snes, SNESFunctionType type)
2045: {
2047:   snes->functype = type;
2048:   return 0;
2049: }

2051: /*@C
2052:    SNESGetFunctionType - Gets the `SNESFunctionType` used in convergence and monitoring set with `SNESSetFunctionType()`
2053:    of the SNES method.

2055:    Logically Collective on snes

2057:    Input Parameters:
2058: +  snes - the `SNES` context
2059: -  type - the type of the function evaluation, see `SNESSetFunctionType()`

2061:    Level: advanced

2063: .seealso: `SNESSetFunctionType()`, `SNESFunctionType`, `SNESSetNormSchedule()`, `SNESComputeFunction()`, `VecNorm()`, `SNESSetFunction()`, `SNESSetInitialFunction()`, `SNESNormSchedule`
2064: @*/
2065: PetscErrorCode SNESGetFunctionType(SNES snes, SNESFunctionType *type)
2066: {
2068:   *type = snes->functype;
2069:   return 0;
2070: }

2072: /*MC
2073:     SNESNGSFunction - function used to apply a Gauss-Seidel sweep on the nonlinear function

2075:      Synopsis:
2076: #include <petscsnes.h>
2077: $    SNESNGSFunction(SNES snes,Vec x,Vec b,void *ctx);

2079:      Collective on snes

2081:      Input Parameters:
2082: +  X   - solution vector
2083: .  B   - RHS vector
2084: -  ctx - optional user-defined Gauss-Seidel context

2086:      Output Parameter:
2087: .  X   - solution vector

2089:    Level: intermediate

2091: .seealso: `SNESNGS`, `SNESSetNGS()`, `SNESGetNGS()`
2092: M*/

2094: /*@C
2095:    SNESSetNGS - Sets the user nonlinear Gauss-Seidel routine for
2096:    use with composed nonlinear solvers.

2098:    Input Parameters:
2099: +  snes   - the SNES context
2100: .  f - function evaluation routine to apply Gauss-Seidel see `SNESNGSFunction`
2101: -  ctx    - [optional] user-defined context for private data for the
2102:             smoother evaluation routine (may be NULL)

2104:    Calling sequence of f:
2105: $  PetscErrorCode f(SNES snes,Vec X,Vec B,void *ctx);

2107:    Arguments of f:
2108: +  snes - the `SNES` context
2109: .  X - the current solution
2110: .  B - the right hand side vector (which may be NULL)
2111: -  ctx - a user provided context

2113:    Note:
2114:    The `SNESNGS` routines are used by the composed nonlinear solver to generate
2115:     a problem appropriate update to the solution, particularly `SNESFAS`.

2117:    Level: intermediate

2119: .seealso: `SNESGetNGS()`, `SNESNGSFunction`, `SNESNCG`, `SNESGetFunction()`, `SNESComputeNGS()`
2120: @*/
2121: PetscErrorCode SNESSetNGS(SNES snes, PetscErrorCode (*f)(SNES, Vec, Vec, void *), void *ctx)
2122: {
2123:   DM dm;

2126:   SNESGetDM(snes, &dm);
2127:   DMSNESSetNGS(dm, f, ctx);
2128:   return 0;
2129: }

2131: /*
2132:      This is used for -snes_mf_operator; it uses a duplicate of snes->jacobian_pre because snes->jacobian_pre cannot be
2133:    changed during the KSPSolve()
2134: */
2135: PetscErrorCode SNESPicardComputeMFFunction(SNES snes, Vec x, Vec f, void *ctx)
2136: {
2137:   DM     dm;
2138:   DMSNES sdm;

2140:   SNESGetDM(snes, &dm);
2141:   DMGetDMSNES(dm, &sdm);
2142:   /*  A(x)*x - b(x) */
2143:   if (sdm->ops->computepfunction) {
2144:     PetscCallBack("SNES Picard callback function", (*sdm->ops->computepfunction)(snes, x, f, sdm->pctx));
2145:     VecScale(f, -1.0);
2146:     /* Cannot share nonzero pattern because of the possible use of SNESComputeJacobianDefault() */
2147:     if (!snes->picard) MatDuplicate(snes->jacobian_pre, MAT_DO_NOT_COPY_VALUES, &snes->picard);
2148:     PetscCallBack("SNES Picard callback Jacobian", (*sdm->ops->computepjacobian)(snes, x, snes->picard, snes->picard, sdm->pctx));
2149:     MatMultAdd(snes->picard, x, f, f);
2150:   } else {
2151:     PetscCallBack("SNES Picard callback Jacobian", (*sdm->ops->computepjacobian)(snes, x, snes->picard, snes->picard, sdm->pctx));
2152:     MatMult(snes->picard, x, f);
2153:   }
2154:   return 0;
2155: }

2157: PetscErrorCode SNESPicardComputeFunction(SNES snes, Vec x, Vec f, void *ctx)
2158: {
2159:   DM     dm;
2160:   DMSNES sdm;

2162:   SNESGetDM(snes, &dm);
2163:   DMGetDMSNES(dm, &sdm);
2164:   /*  A(x)*x - b(x) */
2165:   if (sdm->ops->computepfunction) {
2166:     PetscCallBack("SNES Picard callback function", (*sdm->ops->computepfunction)(snes, x, f, sdm->pctx));
2167:     VecScale(f, -1.0);
2168:     PetscCallBack("SNES Picard callback Jacobian", (*sdm->ops->computepjacobian)(snes, x, snes->jacobian, snes->jacobian_pre, sdm->pctx));
2169:     MatMultAdd(snes->jacobian_pre, x, f, f);
2170:   } else {
2171:     PetscCallBack("SNES Picard callback Jacobian", (*sdm->ops->computepjacobian)(snes, x, snes->jacobian, snes->jacobian_pre, sdm->pctx));
2172:     MatMult(snes->jacobian_pre, x, f);
2173:   }
2174:   return 0;
2175: }

2177: PetscErrorCode SNESPicardComputeJacobian(SNES snes, Vec x1, Mat J, Mat B, void *ctx)
2178: {
2179:   /* the jacobian matrix should be pre-filled in SNESPicardComputeFunction */
2180:   /* must assembly if matrix-free to get the last SNES solution */
2181:   MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY);
2182:   MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY);
2183:   return 0;
2184: }

2186: /*@C
2187:    SNESSetPicard - Use `SNES` to solve the system A(x) x = bp(x) + b via a Picard type iteration (Picard linearization)

2189:    Logically Collective on snes

2191:    Input Parameters:
2192: +  snes - the `SNES` context
2193: .  r - vector to store function values, may be NULL
2194: .  bp - function evaluation routine, may be NULL
2195: .  Amat - matrix with which A(x) x - bp(x) - b is to be computed
2196: .  Pmat - matrix from which preconditioner is computed (usually the same as Amat)
2197: .  J  - function to compute matrix values, see SNESJacobianFunction() for details on its calling sequence
2198: -  ctx - [optional] user-defined context for private data for the function evaluation routine (may be NULL)

2200:    Notes:
2201:     It is often better to provide the nonlinear function F() and some approximation to its Jacobian directly and use
2202:     an approximate Newton solver. This interface is provided to allow porting/testing a previous Picard based code in PETSc before converting it to approximate Newton.

2204:     One can call `SNESSetPicard()` or `SNESSetFunction()` (and possibly `SNESSetJacobian()`) but cannot call both

2206: $     Solves the equation A(x) x = bp(x) - b via the defect correction algorithm A(x^{n}) (x^{n+1} - x^{n}) = bp(x^{n}) + b - A(x^{n})x^{n}
2207: $     Note that when an exact solver is used this corresponds to the "classic" Picard A(x^{n}) x^{n+1} = bp(x^{n}) + b iteration.

2209:      Run with -snes_mf_operator to solve the system with Newton's method using A(x^{n}) to construct the preconditioner.

2211:    We implement the defect correction form of the Picard iteration because it converges much more generally when inexact linear solvers are used then
2212:    the direct Picard iteration A(x^n) x^{n+1} = bp(x^n) + b

2214:    There is some controversity over the definition of a Picard iteration for nonlinear systems but almost everyone agrees that it involves a linear solve and some
2215:    believe it is the iteration  A(x^{n}) x^{n+1} = b(x^{n}) hence we use the name Picard. If anyone has an authoritative  reference that defines the Picard iteration
2216:    different please contact us at petsc-dev@mcs.anl.gov and we'll have an entirely new argument :-).

2218:    When used with -snes_mf_operator this will run matrix-free Newton's method where the matrix-vector product is of the true Jacobian of A(x)x - bp(x) -b and
2219:     A(x^{n}) is used to build the preconditioner

2221:    When used with -snes_fd this will compute the true Jacobian (very slowly one column at at time) and thus represent Newton's method.

2223:    When used with -snes_fd_coloring this will compute the Jacobian via coloring and thus represent a faster implementation of Newton's method. But the
2224:    the nonzero structure of the Jacobian is, in general larger than that of the Picard matrix A so you must provide in A the needed nonzero structure for the correct
2225:    coloring. When using `DMDA` this may mean creating the matrix A with `DMCreateMatrix()` using a wider stencil than strictly needed for A or with a `DMDA_STENCIL_BOX`.
2226:    See the commment in src/snes/tutorials/ex15.c.

2228:    Level: intermediate

2230: .seealso: `SNES`, `SNESGetFunction()`, `SNESSetFunction()`, `SNESComputeFunction()`, `SNESSetJacobian()`, `SNESGetPicard()`, `SNESLineSearchPreCheckPicard()`, `SNESJacobianFunction`
2231: @*/
2232: PetscErrorCode SNESSetPicard(SNES snes, Vec r, PetscErrorCode (*bp)(SNES, Vec, Vec, void *), Mat Amat, Mat Pmat, PetscErrorCode (*J)(SNES, Vec, Mat, Mat, void *), void *ctx)
2233: {
2234:   DM dm;

2237:   SNESGetDM(snes, &dm);
2238:   DMSNESSetPicard(dm, bp, J, ctx);
2239:   DMSNESSetMFFunction(dm, SNESPicardComputeMFFunction, ctx);
2240:   SNESSetFunction(snes, r, SNESPicardComputeFunction, ctx);
2241:   SNESSetJacobian(snes, Amat, Pmat, SNESPicardComputeJacobian, ctx);
2242:   return 0;
2243: }

2245: /*@C
2246:    SNESGetPicard - Returns the context for the Picard iteration

2248:    Not Collective, but `Vec` is parallel if `SNES` is parallel. Collective if `Vec` is requested, but has not been created yet.

2250:    Input Parameter:
2251: .  snes - the `SNES` context

2253:    Output Parameters:
2254: +  r - the function (or NULL)
2255: .  f - the function (or NULL); see `SNESFunction` for calling sequence details
2256: .  Amat - the matrix used to defined the operation A(x) x - b(x) (or NULL)
2257: .  Pmat  - the matrix from which the preconditioner will be constructed (or NULL)
2258: .  J - the function for matrix evaluation (or NULL); see `SNESJacobianFunction` for calling sequence details
2259: -  ctx - the function context (or NULL)

2261:    Level: advanced

2263: .seealso: `SNESSetFunction()`, `SNESSetPicard()`, `SNESGetFunction()`, `SNESGetJacobian()`, `SNESGetDM()`, `SNESFunction`, `SNESJacobianFunction`
2264: @*/
2265: PetscErrorCode SNESGetPicard(SNES snes, Vec *r, PetscErrorCode (**f)(SNES, Vec, Vec, void *), Mat *Amat, Mat *Pmat, PetscErrorCode (**J)(SNES, Vec, Mat, Mat, void *), void **ctx)
2266: {
2267:   DM dm;

2270:   SNESGetFunction(snes, r, NULL, NULL);
2271:   SNESGetJacobian(snes, Amat, Pmat, NULL, NULL);
2272:   SNESGetDM(snes, &dm);
2273:   DMSNESGetPicard(dm, f, J, ctx);
2274:   return 0;
2275: }

2277: /*@C
2278:    SNESSetComputeInitialGuess - Sets a routine used to compute an initial guess for the problem

2280:    Logically Collective on snes

2282:    Input Parameters:
2283: +  snes - the `SNES` context
2284: .  func - function evaluation routine
2285: -  ctx - [optional] user-defined context for private data for the
2286:          function evaluation routine (may be NULL)

2288:    Calling sequence of func:
2289: $    func (SNES snes,Vec x,void *ctx);

2291: .  f - function vector
2292: -  ctx - optional user-defined function context

2294:    Level: intermediate

2296: .seealso: `SNES`, `SNESSolve()`, `SNESGetFunction()`, `SNESComputeFunction()`, `SNESSetJacobian()`
2297: @*/
2298: PetscErrorCode SNESSetComputeInitialGuess(SNES snes, PetscErrorCode (*func)(SNES, Vec, void *), void *ctx)
2299: {
2301:   if (func) snes->ops->computeinitialguess = func;
2302:   if (ctx) snes->initialguessP = ctx;
2303:   return 0;
2304: }

2306: /*@C
2307:    SNESGetRhs - Gets the vector for solving F(x) = rhs. If rhs is not set
2308:    it assumes a zero right hand side.

2310:    Logically Collective on snes

2312:    Input Parameter:
2313: .  snes - the `SNES` context

2315:    Output Parameter:
2316: .  rhs - the right hand side vector or NULL if the right hand side vector is null

2318:    Level: intermediate

2320: .seealso: `SNES`, `SNESGetSolution()`, `SNESGetFunction()`, `SNESComputeFunction()`, `SNESSetJacobian()`, `SNESSetFunction()`
2321: @*/
2322: PetscErrorCode SNESGetRhs(SNES snes, Vec *rhs)
2323: {
2326:   *rhs = snes->vec_rhs;
2327:   return 0;
2328: }

2330: /*@
2331:    SNESComputeFunction - Calls the function that has been set with `SNESSetFunction()`.

2333:    Collective on snes

2335:    Input Parameters:
2336: +  snes - the `SNES` context
2337: -  x - input vector

2339:    Output Parameter:
2340: .  y - function vector, as set by `SNESSetFunction()`

2342:    Note:
2343:    `SNESComputeFunction()` is typically used within nonlinear solvers
2344:    implementations, so users would not generally call this routine themselves.

2346:    Level: developer

2348: .seealso: `SNES`, `SNESSetFunction()`, `SNESGetFunction()`, `SNESComputeMFFunction()`
2349: @*/
2350: PetscErrorCode SNESComputeFunction(SNES snes, Vec x, Vec y)
2351: {
2352:   DM     dm;
2353:   DMSNES sdm;

2360:   VecValidValues_Internal(x, 2, PETSC_TRUE);

2362:   SNESGetDM(snes, &dm);
2363:   DMGetDMSNES(dm, &sdm);
2364:   if (sdm->ops->computefunction) {
2365:     if (sdm->ops->computefunction != SNESObjectiveComputeFunctionDefaultFD) PetscLogEventBegin(SNES_FunctionEval, snes, x, y, 0);
2366:     VecLockReadPush(x);
2367:     /* ensure domainerror is false prior to computefunction evaluation (may not have been reset) */
2368:     snes->domainerror = PETSC_FALSE;
2369:     {
2370:       void *ctx;
2371:       PetscErrorCode (*computefunction)(SNES, Vec, Vec, void *);
2372:       DMSNESGetFunction(dm, &computefunction, &ctx);
2373:       PetscCallBack("SNES callback function", (*computefunction)(snes, x, y, ctx));
2374:     }
2375:     VecLockReadPop(x);
2376:     if (sdm->ops->computefunction != SNESObjectiveComputeFunctionDefaultFD) PetscLogEventEnd(SNES_FunctionEval, snes, x, y, 0);
2377:   } else if (snes->vec_rhs) {
2378:     MatMult(snes->jacobian, x, y);
2379:   } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must call SNESSetFunction() or SNESSetDM() before SNESComputeFunction(), likely called from SNESSolve().");
2380:   if (snes->vec_rhs) VecAXPY(y, -1.0, snes->vec_rhs);
2381:   snes->nfuncs++;
2382:   /*
2383:      domainerror might not be set on all processes; so we tag vector locally with Inf and the next inner product or norm will
2384:      propagate the value to all processes
2385:   */
2386:   if (snes->domainerror) VecSetInf(y);
2387:   return 0;
2388: }

2390: /*@
2391:    SNESComputeMFFunction - Calls the function that has been set with `SNESSetMFFunction()`.

2393:    Collective on snes

2395:    Input Parameters:
2396: +  snes - the `SNES` context
2397: -  x - input vector

2399:    Output Parameter:
2400: .  y - function vector, as set by `SNESSetMFFunction()`

2402:    Notes:
2403:    `SNESComputeMFFunction()` is used within the matrix vector products called by the matrix created with `MatCreateSNESMF()`
2404:    so users would not generally call this routine themselves.

2406:     Since this function is intended for use with finite differencing it does not subtract the right hand side vector provided with `SNESSolve()`
2407:     while `SNESComputeFunction()` does. As such, this routine cannot be used with  `MatMFFDSetBase()` with a provided F function value even if it applies the
2408:     same function as `SNESComputeFunction()` if a `SNESSolve()` right hand side vector is use because the two functions difference would include this right hand side function.

2410:    Level: developer

2412: .seealso: `SNES`, `SNESSetFunction()`, `SNESGetFunction()`, `SNESComputeFunction()`, `MatCreateSNESMF`
2413: @*/
2414: PetscErrorCode SNESComputeMFFunction(SNES snes, Vec x, Vec y)
2415: {
2416:   DM     dm;
2417:   DMSNES sdm;

2424:   VecValidValues_Internal(x, 2, PETSC_TRUE);

2426:   SNESGetDM(snes, &dm);
2427:   DMGetDMSNES(dm, &sdm);
2428:   PetscLogEventBegin(SNES_FunctionEval, snes, x, y, 0);
2429:   VecLockReadPush(x);
2430:   /* ensure domainerror is false prior to computefunction evaluation (may not have been reset) */
2431:   snes->domainerror = PETSC_FALSE;
2432:   PetscCallBack("SNES callback function", (*sdm->ops->computemffunction)(snes, x, y, sdm->mffunctionctx));
2433:   VecLockReadPop(x);
2434:   PetscLogEventEnd(SNES_FunctionEval, snes, x, y, 0);
2435:   snes->nfuncs++;
2436:   /*
2437:      domainerror might not be set on all processes; so we tag vector locally with Inf and the next inner product or norm will
2438:      propagate the value to all processes
2439:   */
2440:   if (snes->domainerror) VecSetInf(y);
2441:   return 0;
2442: }

2444: /*@
2445:    SNESComputeNGS - Calls the Gauss-Seidel function that has been set with  `SNESSetNGS()`.

2447:    Collective on snes

2449:    Input Parameters:
2450: +  snes - the `SNES` context
2451: .  x - input vector
2452: -  b - rhs vector

2454:    Output Parameter:
2455: .  x - new solution vector

2457:    Note:
2458:    `SNESComputeNGS()` is typically used within composed nonlinear solver
2459:    implementations, so most users would not generally call this routine
2460:    themselves.

2462:    Level: developer

2464: .seealso: `SNESNGS`, `SNESSetNGS()`, `SNESComputeFunction()`
2465: @*/
2466: PetscErrorCode SNESComputeNGS(SNES snes, Vec b, Vec x)
2467: {
2468:   DM     dm;
2469:   DMSNES sdm;

2476:   if (b) VecValidValues_Internal(b, 2, PETSC_TRUE);
2477:   PetscLogEventBegin(SNES_NGSEval, snes, x, b, 0);
2478:   SNESGetDM(snes, &dm);
2479:   DMGetDMSNES(dm, &sdm);
2480:   if (sdm->ops->computegs) {
2481:     if (b) VecLockReadPush(b);
2482:     PetscCallBack("SNES callback NGS", (*sdm->ops->computegs)(snes, x, b, sdm->gsctx));
2483:     if (b) VecLockReadPop(b);
2484:   } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must call SNESSetNGS() before SNESComputeNGS(), likely called from SNESSolve().");
2485:   PetscLogEventEnd(SNES_NGSEval, snes, x, b, 0);
2486:   return 0;
2487: }

2489: PetscErrorCode SNESTestJacobian(SNES snes)
2490: {
2491:   Mat               A, B, C, D, jacobian;
2492:   Vec               x = snes->vec_sol, f = snes->vec_func;
2493:   PetscReal         nrm, gnorm;
2494:   PetscReal         threshold = 1.e-5;
2495:   MatType           mattype;
2496:   PetscInt          m, n, M, N;
2497:   void             *functx;
2498:   PetscBool         complete_print = PETSC_FALSE, threshold_print = PETSC_FALSE, test = PETSC_FALSE, flg, istranspose;
2499:   PetscViewer       viewer, mviewer;
2500:   MPI_Comm          comm;
2501:   PetscInt          tabs;
2502:   static PetscBool  directionsprinted = PETSC_FALSE;
2503:   PetscViewerFormat format;

2505:   PetscObjectOptionsBegin((PetscObject)snes);
2506:   PetscOptionsName("-snes_test_jacobian", "Compare hand-coded and finite difference Jacobians", "None", &test);
2507:   PetscOptionsReal("-snes_test_jacobian", "Threshold for element difference between hand-coded and finite difference being meaningful", "None", threshold, &threshold, NULL);
2508:   PetscOptionsViewer("-snes_test_jacobian_view", "View difference between hand-coded and finite difference Jacobians element entries", "None", &mviewer, &format, &complete_print);
2509:   if (!complete_print) {
2510:     PetscOptionsDeprecated("-snes_test_jacobian_display", "-snes_test_jacobian_view", "3.13", NULL);
2511:     PetscOptionsViewer("-snes_test_jacobian_display", "Display difference between hand-coded and finite difference Jacobians", "None", &mviewer, &format, &complete_print);
2512:   }
2513:   /* for compatibility with PETSc 3.9 and older. */
2514:   PetscOptionsDeprecated("-snes_test_jacobian_display_threshold", "-snes_test_jacobian", "3.13", "-snes_test_jacobian accepts an optional threshold (since v3.10)");
2515:   PetscOptionsReal("-snes_test_jacobian_display_threshold", "Display difference between hand-coded and finite difference Jacobians which exceed input threshold", "None", threshold, &threshold, &threshold_print);
2516:   PetscOptionsEnd();
2517:   if (!test) return 0;

2519:   PetscObjectGetComm((PetscObject)snes, &comm);
2520:   PetscViewerASCIIGetStdout(comm, &viewer);
2521:   PetscViewerASCIIGetTab(viewer, &tabs);
2522:   PetscViewerASCIISetTab(viewer, ((PetscObject)snes)->tablevel);
2523:   PetscViewerASCIIPrintf(viewer, "  ---------- Testing Jacobian -------------\n");
2524:   if (!complete_print && !directionsprinted) {
2525:     PetscViewerASCIIPrintf(viewer, "  Run with -snes_test_jacobian_view and optionally -snes_test_jacobian <threshold> to show difference\n");
2526:     PetscViewerASCIIPrintf(viewer, "    of hand-coded and finite difference Jacobian entries greater than <threshold>.\n");
2527:   }
2528:   if (!directionsprinted) {
2529:     PetscViewerASCIIPrintf(viewer, "  Testing hand-coded Jacobian, if (for double precision runs) ||J - Jfd||_F/||J||_F is\n");
2530:     PetscViewerASCIIPrintf(viewer, "    O(1.e-8), the hand-coded Jacobian is probably correct.\n");
2531:     directionsprinted = PETSC_TRUE;
2532:   }
2533:   if (complete_print) PetscViewerPushFormat(mviewer, format);

2535:   PetscObjectTypeCompare((PetscObject)snes->jacobian, MATMFFD, &flg);
2536:   if (!flg) jacobian = snes->jacobian;
2537:   else jacobian = snes->jacobian_pre;

2539:   if (!x) {
2540:     MatCreateVecs(jacobian, &x, NULL);
2541:   } else {
2542:     PetscObjectReference((PetscObject)x);
2543:   }
2544:   if (!f) {
2545:     VecDuplicate(x, &f);
2546:   } else {
2547:     PetscObjectReference((PetscObject)f);
2548:   }
2549:   /* evaluate the function at this point because SNESComputeJacobianDefault() assumes that the function has been evaluated and put into snes->vec_func */
2550:   SNESComputeFunction(snes, x, f);
2551:   VecDestroy(&f);
2552:   PetscObjectTypeCompare((PetscObject)snes, SNESKSPTRANSPOSEONLY, &istranspose);
2553:   while (jacobian) {
2554:     Mat JT = NULL, Jsave = NULL;

2556:     if (istranspose) {
2557:       MatCreateTranspose(jacobian, &JT);
2558:       Jsave    = jacobian;
2559:       jacobian = JT;
2560:     }
2561:     PetscObjectBaseTypeCompareAny((PetscObject)jacobian, &flg, MATSEQAIJ, MATMPIAIJ, MATSEQDENSE, MATMPIDENSE, MATSEQBAIJ, MATMPIBAIJ, MATSEQSBAIJ, MATMPISBAIJ, "");
2562:     if (flg) {
2563:       A = jacobian;
2564:       PetscObjectReference((PetscObject)A);
2565:     } else {
2566:       MatComputeOperator(jacobian, MATAIJ, &A);
2567:     }

2569:     MatGetType(A, &mattype);
2570:     MatGetSize(A, &M, &N);
2571:     MatGetLocalSize(A, &m, &n);
2572:     MatCreate(PetscObjectComm((PetscObject)A), &B);
2573:     MatSetType(B, mattype);
2574:     MatSetSizes(B, m, n, M, N);
2575:     MatSetBlockSizesFromMats(B, A, A);
2576:     MatSetUp(B);
2577:     MatSetOption(B, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);

2579:     SNESGetFunction(snes, NULL, NULL, &functx);
2580:     SNESComputeJacobianDefault(snes, x, B, B, functx);

2582:     MatDuplicate(B, MAT_COPY_VALUES, &D);
2583:     MatAYPX(D, -1.0, A, DIFFERENT_NONZERO_PATTERN);
2584:     MatNorm(D, NORM_FROBENIUS, &nrm);
2585:     MatNorm(A, NORM_FROBENIUS, &gnorm);
2586:     MatDestroy(&D);
2587:     if (!gnorm) gnorm = 1; /* just in case */
2588:     PetscViewerASCIIPrintf(viewer, "  ||J - Jfd||_F/||J||_F = %g, ||J - Jfd||_F = %g\n", (double)(nrm / gnorm), (double)nrm);

2590:     if (complete_print) {
2591:       PetscViewerASCIIPrintf(viewer, "  Hand-coded Jacobian ----------\n");
2592:       MatView(A, mviewer);
2593:       PetscViewerASCIIPrintf(viewer, "  Finite difference Jacobian ----------\n");
2594:       MatView(B, mviewer);
2595:     }

2597:     if (threshold_print || complete_print) {
2598:       PetscInt           Istart, Iend, *ccols, bncols, cncols, j, row;
2599:       PetscScalar       *cvals;
2600:       const PetscInt    *bcols;
2601:       const PetscScalar *bvals;

2603:       MatCreate(PetscObjectComm((PetscObject)A), &C);
2604:       MatSetType(C, mattype);
2605:       MatSetSizes(C, m, n, M, N);
2606:       MatSetBlockSizesFromMats(C, A, A);
2607:       MatSetUp(C);
2608:       MatSetOption(C, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);

2610:       MatAYPX(B, -1.0, A, DIFFERENT_NONZERO_PATTERN);
2611:       MatGetOwnershipRange(B, &Istart, &Iend);

2613:       for (row = Istart; row < Iend; row++) {
2614:         MatGetRow(B, row, &bncols, &bcols, &bvals);
2615:         PetscMalloc2(bncols, &ccols, bncols, &cvals);
2616:         for (j = 0, cncols = 0; j < bncols; j++) {
2617:           if (PetscAbsScalar(bvals[j]) > threshold) {
2618:             ccols[cncols] = bcols[j];
2619:             cvals[cncols] = bvals[j];
2620:             cncols += 1;
2621:           }
2622:         }
2623:         if (cncols) MatSetValues(C, 1, &row, cncols, ccols, cvals, INSERT_VALUES);
2624:         MatRestoreRow(B, row, &bncols, &bcols, &bvals);
2625:         PetscFree2(ccols, cvals);
2626:       }
2627:       MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY);
2628:       MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY);
2629:       PetscViewerASCIIPrintf(viewer, "  Hand-coded minus finite-difference Jacobian with tolerance %g ----------\n", (double)threshold);
2630:       MatView(C, complete_print ? mviewer : viewer);
2631:       MatDestroy(&C);
2632:     }
2633:     MatDestroy(&A);
2634:     MatDestroy(&B);
2635:     MatDestroy(&JT);
2636:     if (Jsave) jacobian = Jsave;
2637:     if (jacobian != snes->jacobian_pre) {
2638:       jacobian = snes->jacobian_pre;
2639:       PetscViewerASCIIPrintf(viewer, "  ---------- Testing Jacobian for preconditioner -------------\n");
2640:     } else jacobian = NULL;
2641:   }
2642:   VecDestroy(&x);
2643:   if (complete_print) PetscViewerPopFormat(mviewer);
2644:   if (mviewer) PetscViewerDestroy(&mviewer);
2645:   PetscViewerASCIISetTab(viewer, tabs);
2646:   return 0;
2647: }

2649: /*@
2650:    SNESComputeJacobian - Computes the Jacobian matrix that has been set with `SNESSetJacobian()`.

2652:    Collective on snes

2654:    Input Parameters:
2655: +  snes - the `SNES` context
2656: -  x - input vector

2658:    Output Parameters:
2659: +  A - Jacobian matrix
2660: -  B - optional matrix for building the preconditioner

2662:   Options Database Keys:
2663: +    -snes_lag_preconditioner <lag> - how often to rebuild preconditioner
2664: .    -snes_lag_jacobian <lag> - how often to rebuild Jacobian
2665: .    -snes_test_jacobian <optional threshold> - compare the user provided Jacobian with one compute via finite differences to check for errors.  If a threshold is given, display only those entries whose difference is greater than the threshold.
2666: .    -snes_test_jacobian_view - display the user provided Jacobian, the finite difference Jacobian and the difference between them to help users detect the location of errors in the user provided Jacobian
2667: .    -snes_compare_explicit - Compare the computed Jacobian to the finite difference Jacobian and output the differences
2668: .    -snes_compare_explicit_draw  - Compare the computed Jacobian to the finite difference Jacobian and draw the result
2669: .    -snes_compare_explicit_contour  - Compare the computed Jacobian to the finite difference Jacobian and draw a contour plot with the result
2670: .    -snes_compare_operator  - Make the comparison options above use the operator instead of the preconditioning matrix
2671: .    -snes_compare_coloring - Compute the finite difference Jacobian using coloring and display norms of difference
2672: .    -snes_compare_coloring_display - Compute the finite difference Jacobian using coloring and display verbose differences
2673: .    -snes_compare_coloring_threshold - Display only those matrix entries that differ by more than a given threshold
2674: .    -snes_compare_coloring_threshold_atol - Absolute tolerance for difference in matrix entries to be displayed by -snes_compare_coloring_threshold
2675: .    -snes_compare_coloring_threshold_rtol - Relative tolerance for difference in matrix entries to be displayed by -snes_compare_coloring_threshold
2676: .    -snes_compare_coloring_draw - Compute the finite difference Jacobian using coloring and draw differences
2677: -    -snes_compare_coloring_draw_contour - Compute the finite difference Jacobian using coloring and show contours of matrices and differences

2679:    Note:
2680:    Most users should not need to explicitly call this routine, as it
2681:    is used internally within the nonlinear solvers.

2683:    Developer Note:
2684:     This has duplicative ways of checking the accuracy of the user provided Jacobian (see the options above). This is for historical reasons, the routine SNESTestJacobian() use to used
2685:       for with the SNESType of test that has been removed.

2687:    Level: developer

2689: .seealso: `SNESSetJacobian()`, `KSPSetOperators()`, `MatStructure`, `SNESSetLagPreconditioner()`, `SNESSetLagJacobian()`
2690: @*/
2691: PetscErrorCode SNESComputeJacobian(SNES snes, Vec X, Mat A, Mat B)
2692: {
2693:   PetscBool flag;
2694:   DM        dm;
2695:   DMSNES    sdm;
2696:   KSP       ksp;

2701:   VecValidValues_Internal(X, 2, PETSC_TRUE);
2702:   SNESGetDM(snes, &dm);
2703:   DMGetDMSNES(dm, &sdm);

2705:   /* make sure that MatAssemblyBegin/End() is called on A matrix if it is matrix free */
2706:   if (snes->lagjacobian == -2) {
2707:     snes->lagjacobian = -1;

2709:     PetscInfo(snes, "Recomputing Jacobian/preconditioner because lag is -2 (means compute Jacobian, but then never again) \n");
2710:   } else if (snes->lagjacobian == -1) {
2711:     PetscInfo(snes, "Reusing Jacobian/preconditioner because lag is -1\n");
2712:     PetscObjectTypeCompare((PetscObject)A, MATMFFD, &flag);
2713:     if (flag) {
2714:       MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
2715:       MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
2716:     }
2717:     return 0;
2718:   } else if (snes->lagjacobian > 1 && (snes->iter + snes->jac_iter) % snes->lagjacobian) {
2719:     PetscInfo(snes, "Reusing Jacobian/preconditioner because lag is %" PetscInt_FMT " and SNES iteration is %" PetscInt_FMT "\n", snes->lagjacobian, snes->iter);
2720:     PetscObjectTypeCompare((PetscObject)A, MATMFFD, &flag);
2721:     if (flag) {
2722:       MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
2723:       MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
2724:     }
2725:     return 0;
2726:   }
2727:   if (snes->npc && snes->npcside == PC_LEFT) {
2728:     MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
2729:     MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
2730:     return 0;
2731:   }

2733:   PetscLogEventBegin(SNES_JacobianEval, snes, X, A, B);
2734:   VecLockReadPush(X);
2735:   {
2736:     void *ctx;
2737:     PetscErrorCode (*J)(SNES, Vec, Mat, Mat, void *);
2738:     DMSNESGetJacobian(dm, &J, &ctx);
2739:     PetscCallBack("SNES callback Jacobian", (*J)(snes, X, A, B, ctx));
2740:   }
2741:   VecLockReadPop(X);
2742:   PetscLogEventEnd(SNES_JacobianEval, snes, X, A, B);

2744:   /* attach latest linearization point to the preconditioning matrix */
2745:   PetscObjectCompose((PetscObject)B, "__SNES_latest_X", (PetscObject)X);

2747:   /* the next line ensures that snes->ksp exists */
2748:   SNESGetKSP(snes, &ksp);
2749:   if (snes->lagpreconditioner == -2) {
2750:     PetscInfo(snes, "Rebuilding preconditioner exactly once since lag is -2\n");
2751:     KSPSetReusePreconditioner(snes->ksp, PETSC_FALSE);
2752:     snes->lagpreconditioner = -1;
2753:   } else if (snes->lagpreconditioner == -1) {
2754:     PetscInfo(snes, "Reusing preconditioner because lag is -1\n");
2755:     KSPSetReusePreconditioner(snes->ksp, PETSC_TRUE);
2756:   } else if (snes->lagpreconditioner > 1 && (snes->iter + snes->pre_iter) % snes->lagpreconditioner) {
2757:     PetscInfo(snes, "Reusing preconditioner because lag is %" PetscInt_FMT " and SNES iteration is %" PetscInt_FMT "\n", snes->lagpreconditioner, snes->iter);
2758:     KSPSetReusePreconditioner(snes->ksp, PETSC_TRUE);
2759:   } else {
2760:     PetscInfo(snes, "Rebuilding preconditioner\n");
2761:     KSPSetReusePreconditioner(snes->ksp, PETSC_FALSE);
2762:   }

2764:   SNESTestJacobian(snes);
2765:   /* make sure user returned a correct Jacobian and preconditioner */
2768:   {
2769:     PetscBool flag = PETSC_FALSE, flag_draw = PETSC_FALSE, flag_contour = PETSC_FALSE, flag_operator = PETSC_FALSE;
2770:     PetscOptionsGetViewer(PetscObjectComm((PetscObject)snes), ((PetscObject)snes)->options, ((PetscObject)snes)->prefix, "-snes_compare_explicit", NULL, NULL, &flag);
2771:     PetscOptionsGetViewer(PetscObjectComm((PetscObject)snes), ((PetscObject)snes)->options, ((PetscObject)snes)->prefix, "-snes_compare_explicit_draw", NULL, NULL, &flag_draw);
2772:     PetscOptionsGetViewer(PetscObjectComm((PetscObject)snes), ((PetscObject)snes)->options, ((PetscObject)snes)->prefix, "-snes_compare_explicit_draw_contour", NULL, NULL, &flag_contour);
2773:     PetscOptionsGetViewer(PetscObjectComm((PetscObject)snes), ((PetscObject)snes)->options, ((PetscObject)snes)->prefix, "-snes_compare_operator", NULL, NULL, &flag_operator);
2774:     if (flag || flag_draw || flag_contour) {
2775:       Mat         Bexp_mine = NULL, Bexp, FDexp;
2776:       PetscViewer vdraw, vstdout;
2777:       PetscBool   flg;
2778:       if (flag_operator) {
2779:         MatComputeOperator(A, MATAIJ, &Bexp_mine);
2780:         Bexp = Bexp_mine;
2781:       } else {
2782:         /* See if the preconditioning matrix can be viewed and added directly */
2783:         PetscObjectBaseTypeCompareAny((PetscObject)B, &flg, MATSEQAIJ, MATMPIAIJ, MATSEQDENSE, MATMPIDENSE, MATSEQBAIJ, MATMPIBAIJ, MATSEQSBAIJ, MATMPIBAIJ, "");
2784:         if (flg) Bexp = B;
2785:         else {
2786:           /* If the "preconditioning" matrix is itself MATSHELL or some other type without direct support */
2787:           MatComputeOperator(B, MATAIJ, &Bexp_mine);
2788:           Bexp = Bexp_mine;
2789:         }
2790:       }
2791:       MatConvert(Bexp, MATSAME, MAT_INITIAL_MATRIX, &FDexp);
2792:       SNESComputeJacobianDefault(snes, X, FDexp, FDexp, NULL);
2793:       PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)snes), &vstdout);
2794:       if (flag_draw || flag_contour) {
2795:         PetscViewerDrawOpen(PetscObjectComm((PetscObject)snes), NULL, "Explicit Jacobians", PETSC_DECIDE, PETSC_DECIDE, 300, 300, &vdraw);
2796:         if (flag_contour) PetscViewerPushFormat(vdraw, PETSC_VIEWER_DRAW_CONTOUR);
2797:       } else vdraw = NULL;
2798:       PetscViewerASCIIPrintf(vstdout, "Explicit %s\n", flag_operator ? "Jacobian" : "preconditioning Jacobian");
2799:       if (flag) MatView(Bexp, vstdout);
2800:       if (vdraw) MatView(Bexp, vdraw);
2801:       PetscViewerASCIIPrintf(vstdout, "Finite difference Jacobian\n");
2802:       if (flag) MatView(FDexp, vstdout);
2803:       if (vdraw) MatView(FDexp, vdraw);
2804:       MatAYPX(FDexp, -1.0, Bexp, SAME_NONZERO_PATTERN);
2805:       PetscViewerASCIIPrintf(vstdout, "User-provided matrix minus finite difference Jacobian\n");
2806:       if (flag) MatView(FDexp, vstdout);
2807:       if (vdraw) { /* Always use contour for the difference */
2808:         PetscViewerPushFormat(vdraw, PETSC_VIEWER_DRAW_CONTOUR);
2809:         MatView(FDexp, vdraw);
2810:         PetscViewerPopFormat(vdraw);
2811:       }
2812:       if (flag_contour) PetscViewerPopFormat(vdraw);
2813:       PetscViewerDestroy(&vdraw);
2814:       MatDestroy(&Bexp_mine);
2815:       MatDestroy(&FDexp);
2816:     }
2817:   }
2818:   {
2819:     PetscBool flag = PETSC_FALSE, flag_display = PETSC_FALSE, flag_draw = PETSC_FALSE, flag_contour = PETSC_FALSE, flag_threshold = PETSC_FALSE;
2820:     PetscReal threshold_atol = PETSC_SQRT_MACHINE_EPSILON, threshold_rtol = 10 * PETSC_SQRT_MACHINE_EPSILON;
2821:     PetscOptionsGetViewer(PetscObjectComm((PetscObject)snes), ((PetscObject)snes)->options, ((PetscObject)snes)->prefix, "-snes_compare_coloring", NULL, NULL, &flag);
2822:     PetscOptionsGetViewer(PetscObjectComm((PetscObject)snes), ((PetscObject)snes)->options, ((PetscObject)snes)->prefix, "-snes_compare_coloring_display", NULL, NULL, &flag_display);
2823:     PetscOptionsGetViewer(PetscObjectComm((PetscObject)snes), ((PetscObject)snes)->options, ((PetscObject)snes)->prefix, "-snes_compare_coloring_draw", NULL, NULL, &flag_draw);
2824:     PetscOptionsGetViewer(PetscObjectComm((PetscObject)snes), ((PetscObject)snes)->options, ((PetscObject)snes)->prefix, "-snes_compare_coloring_draw_contour", NULL, NULL, &flag_contour);
2825:     PetscOptionsGetViewer(PetscObjectComm((PetscObject)snes), ((PetscObject)snes)->options, ((PetscObject)snes)->prefix, "-snes_compare_coloring_threshold", NULL, NULL, &flag_threshold);
2826:     if (flag_threshold) {
2827:       PetscOptionsGetReal(((PetscObject)snes)->options, ((PetscObject)snes)->prefix, "-snes_compare_coloring_threshold_rtol", &threshold_rtol, NULL);
2828:       PetscOptionsGetReal(((PetscObject)snes)->options, ((PetscObject)snes)->prefix, "-snes_compare_coloring_threshold_atol", &threshold_atol, NULL);
2829:     }
2830:     if (flag || flag_display || flag_draw || flag_contour || flag_threshold) {
2831:       Mat           Bfd;
2832:       PetscViewer   vdraw, vstdout;
2833:       MatColoring   coloring;
2834:       ISColoring    iscoloring;
2835:       MatFDColoring matfdcoloring;
2836:       PetscErrorCode (*func)(SNES, Vec, Vec, void *);
2837:       void     *funcctx;
2838:       PetscReal norm1, norm2, normmax;

2840:       MatDuplicate(B, MAT_DO_NOT_COPY_VALUES, &Bfd);
2841:       MatColoringCreate(Bfd, &coloring);
2842:       MatColoringSetType(coloring, MATCOLORINGSL);
2843:       MatColoringSetFromOptions(coloring);
2844:       MatColoringApply(coloring, &iscoloring);
2845:       MatColoringDestroy(&coloring);
2846:       MatFDColoringCreate(Bfd, iscoloring, &matfdcoloring);
2847:       MatFDColoringSetFromOptions(matfdcoloring);
2848:       MatFDColoringSetUp(Bfd, iscoloring, matfdcoloring);
2849:       ISColoringDestroy(&iscoloring);

2851:       /* This method of getting the function is currently unreliable since it doesn't work for DM local functions. */
2852:       SNESGetFunction(snes, NULL, &func, &funcctx);
2853:       MatFDColoringSetFunction(matfdcoloring, (PetscErrorCode(*)(void))func, funcctx);
2854:       PetscObjectSetOptionsPrefix((PetscObject)matfdcoloring, ((PetscObject)snes)->prefix);
2855:       PetscObjectAppendOptionsPrefix((PetscObject)matfdcoloring, "coloring_");
2856:       MatFDColoringSetFromOptions(matfdcoloring);
2857:       MatFDColoringApply(Bfd, matfdcoloring, X, snes);
2858:       MatFDColoringDestroy(&matfdcoloring);

2860:       PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)snes), &vstdout);
2861:       if (flag_draw || flag_contour) {
2862:         PetscViewerDrawOpen(PetscObjectComm((PetscObject)snes), NULL, "Colored Jacobians", PETSC_DECIDE, PETSC_DECIDE, 300, 300, &vdraw);
2863:         if (flag_contour) PetscViewerPushFormat(vdraw, PETSC_VIEWER_DRAW_CONTOUR);
2864:       } else vdraw = NULL;
2865:       PetscViewerASCIIPrintf(vstdout, "Explicit preconditioning Jacobian\n");
2866:       if (flag_display) MatView(B, vstdout);
2867:       if (vdraw) MatView(B, vdraw);
2868:       PetscViewerASCIIPrintf(vstdout, "Colored Finite difference Jacobian\n");
2869:       if (flag_display) MatView(Bfd, vstdout);
2870:       if (vdraw) MatView(Bfd, vdraw);
2871:       MatAYPX(Bfd, -1.0, B, SAME_NONZERO_PATTERN);
2872:       MatNorm(Bfd, NORM_1, &norm1);
2873:       MatNorm(Bfd, NORM_FROBENIUS, &norm2);
2874:       MatNorm(Bfd, NORM_MAX, &normmax);
2875:       PetscViewerASCIIPrintf(vstdout, "User-provided matrix minus finite difference Jacobian, norm1=%g normFrob=%g normmax=%g\n", (double)norm1, (double)norm2, (double)normmax);
2876:       if (flag_display) MatView(Bfd, vstdout);
2877:       if (vdraw) { /* Always use contour for the difference */
2878:         PetscViewerPushFormat(vdraw, PETSC_VIEWER_DRAW_CONTOUR);
2879:         MatView(Bfd, vdraw);
2880:         PetscViewerPopFormat(vdraw);
2881:       }
2882:       if (flag_contour) PetscViewerPopFormat(vdraw);

2884:       if (flag_threshold) {
2885:         PetscInt bs, rstart, rend, i;
2886:         MatGetBlockSize(B, &bs);
2887:         MatGetOwnershipRange(B, &rstart, &rend);
2888:         for (i = rstart; i < rend; i++) {
2889:           const PetscScalar *ba, *ca;
2890:           const PetscInt    *bj, *cj;
2891:           PetscInt           bn, cn, j, maxentrycol = -1, maxdiffcol = -1, maxrdiffcol = -1;
2892:           PetscReal          maxentry = 0, maxdiff = 0, maxrdiff = 0;
2893:           MatGetRow(B, i, &bn, &bj, &ba);
2894:           MatGetRow(Bfd, i, &cn, &cj, &ca);
2896:           for (j = 0; j < bn; j++) {
2897:             PetscReal rdiff = PetscAbsScalar(ca[j]) / (threshold_atol + threshold_rtol * PetscAbsScalar(ba[j]));
2898:             if (PetscAbsScalar(ba[j]) > PetscAbs(maxentry)) {
2899:               maxentrycol = bj[j];
2900:               maxentry    = PetscRealPart(ba[j]);
2901:             }
2902:             if (PetscAbsScalar(ca[j]) > PetscAbs(maxdiff)) {
2903:               maxdiffcol = bj[j];
2904:               maxdiff    = PetscRealPart(ca[j]);
2905:             }
2906:             if (rdiff > maxrdiff) {
2907:               maxrdiffcol = bj[j];
2908:               maxrdiff    = rdiff;
2909:             }
2910:           }
2911:           if (maxrdiff > 1) {
2912:             PetscViewerASCIIPrintf(vstdout, "row %" PetscInt_FMT " (maxentry=%g at %" PetscInt_FMT ", maxdiff=%g at %" PetscInt_FMT ", maxrdiff=%g at %" PetscInt_FMT "):", i, (double)maxentry, maxentrycol, (double)maxdiff, maxdiffcol, (double)maxrdiff, maxrdiffcol);
2913:             for (j = 0; j < bn; j++) {
2914:               PetscReal rdiff;
2915:               rdiff = PetscAbsScalar(ca[j]) / (threshold_atol + threshold_rtol * PetscAbsScalar(ba[j]));
2916:               if (rdiff > 1) PetscViewerASCIIPrintf(vstdout, " (%" PetscInt_FMT ",%g:%g)", bj[j], (double)PetscRealPart(ba[j]), (double)PetscRealPart(ca[j]));
2917:             }
2918:             PetscViewerASCIIPrintf(vstdout, "\n");
2919:           }
2920:           MatRestoreRow(B, i, &bn, &bj, &ba);
2921:           MatRestoreRow(Bfd, i, &cn, &cj, &ca);
2922:         }
2923:       }
2924:       PetscViewerDestroy(&vdraw);
2925:       MatDestroy(&Bfd);
2926:     }
2927:   }
2928:   return 0;
2929: }

2931: /*MC
2932:     SNESJacobianFunction - Function used by `SNES` to compute the nonlinear Jacobian of the function to be solved by `SNES`

2934:      Synopsis:
2935:      #include "petscsnes.h"
2936:      PetscErrorCode SNESJacobianFunction(SNES snes,Vec x,Mat Amat,Mat Pmat,void *ctx);

2938:      Collective on snes

2940:     Input Parameters:
2941: +  x - input vector, the Jacobian is to be computed at this value
2942: -  ctx - [optional] user-defined Jacobian context

2944:     Output Parameters:
2945: +  Amat - the matrix that defines the (approximate) Jacobian
2946: -  Pmat - the matrix to be used in constructing the preconditioner, usually the same as Amat.

2948:    Level: intermediate

2950: .seealso: `SNES`, `SNESSetFunction()`, `SNESGetFunction()`, `SNESSetJacobian()`, `SNESGetJacobian()`
2951: M*/

2953: /*@C
2954:    SNESSetJacobian - Sets the function to compute Jacobian as well as the
2955:    location to store the matrix.

2957:    Logically Collective on snes

2959:    Input Parameters:
2960: +  snes - the `SNES` context
2961: .  Amat - the matrix that defines the (approximate) Jacobian
2962: .  Pmat - the matrix to be used in constructing the preconditioner, usually the same as Amat.
2963: .  J - Jacobian evaluation routine (if NULL then `SNES` retains any previously set value), see `SNESJacobianFunction` for details
2964: -  ctx - [optional] user-defined context for private data for the
2965:          Jacobian evaluation routine (may be NULL) (if NULL then SNES retains any previously set value)

2967:    Notes:
2968:    If the Amat matrix and Pmat matrix are different you must call `MatAssemblyBegin()`/`MatAssemblyEnd()` on
2969:    each matrix.

2971:    If you know the operator Amat has a null space you can use `MatSetNullSpace()` and `MatSetTransposeNullSpace()` to supply the null
2972:    space to Amat and the KSP solvers will automatically use that null space as needed during the solution process.

2974:    If using `SNESComputeJacobianDefaultColor()` to assemble a Jacobian, the ctx argument
2975:    must be a `MatFDColoring`.

2977:    Other defect-correction schemes can be used by computing a different matrix in place of the Jacobian.  One common
2978:    example is to use the "Picard linearization" which only differentiates through the highest order parts of each term using `SNESSetPicard()`

2980:    Level: beginner

2982: .seealso: `SNES`, `KSPSetOperators()`, `SNESSetFunction()`, `MatMFFDComputeJacobian()`, `SNESComputeJacobianDefaultColor()`, `MatStructure`, `J`,
2983:           `SNESSetPicard()`, `SNESJacobianFunction`
2984: @*/
2985: PetscErrorCode SNESSetJacobian(SNES snes, Mat Amat, Mat Pmat, PetscErrorCode (*J)(SNES, Vec, Mat, Mat, void *), void *ctx)
2986: {
2987:   DM dm;

2994:   SNESGetDM(snes, &dm);
2995:   DMSNESSetJacobian(dm, J, ctx);
2996:   if (Amat) {
2997:     PetscObjectReference((PetscObject)Amat);
2998:     MatDestroy(&snes->jacobian);

3000:     snes->jacobian = Amat;
3001:   }
3002:   if (Pmat) {
3003:     PetscObjectReference((PetscObject)Pmat);
3004:     MatDestroy(&snes->jacobian_pre);

3006:     snes->jacobian_pre = Pmat;
3007:   }
3008:   return 0;
3009: }

3011: /*@C
3012:    SNESGetJacobian - Returns the Jacobian matrix and optionally the user
3013:    provided context for evaluating the Jacobian.

3015:    Not Collective, but `Mat` object will be parallel if `SNES` object is

3017:    Input Parameter:
3018: .  snes - the nonlinear solver context

3020:    Output Parameters:
3021: +  Amat - location to stash (approximate) Jacobian matrix (or NULL)
3022: .  Pmat - location to stash matrix used to compute the preconditioner (or NULL)
3023: .  J - location to put Jacobian function (or NULL), see SNESJacobianFunction for details on its calling sequence
3024: -  ctx - location to stash Jacobian ctx (or NULL)

3026:    Level: advanced

3028: .seealso: `SNES`, `Mat`, `SNESSetJacobian()`, `SNESComputeJacobian()`, `SNESJacobianFunction`, `SNESGetFunction()`
3029: @*/
3030: PetscErrorCode SNESGetJacobian(SNES snes, Mat *Amat, Mat *Pmat, PetscErrorCode (**J)(SNES, Vec, Mat, Mat, void *), void **ctx)
3031: {
3032:   DM dm;

3035:   if (Amat) *Amat = snes->jacobian;
3036:   if (Pmat) *Pmat = snes->jacobian_pre;
3037:   SNESGetDM(snes, &dm);
3038:   DMSNESGetJacobian(dm, J, ctx);
3039:   return 0;
3040: }

3042: static PetscErrorCode SNESSetDefaultComputeJacobian(SNES snes)
3043: {
3044:   DM     dm;
3045:   DMSNES sdm;

3047:   SNESGetDM(snes, &dm);
3048:   DMGetDMSNES(dm, &sdm);
3049:   if (!sdm->ops->computejacobian && snes->jacobian_pre) {
3050:     DM        dm;
3051:     PetscBool isdense, ismf;

3053:     SNESGetDM(snes, &dm);
3054:     PetscObjectTypeCompareAny((PetscObject)snes->jacobian_pre, &isdense, MATSEQDENSE, MATMPIDENSE, MATDENSE, NULL);
3055:     PetscObjectTypeCompareAny((PetscObject)snes->jacobian_pre, &ismf, MATMFFD, MATSHELL, NULL);
3056:     if (isdense) {
3057:       DMSNESSetJacobian(dm, SNESComputeJacobianDefault, NULL);
3058:     } else if (!ismf) {
3059:       DMSNESSetJacobian(dm, SNESComputeJacobianDefaultColor, NULL);
3060:     }
3061:   }
3062:   return 0;
3063: }

3065: /*@
3066:    SNESSetUp - Sets up the internal data structures for the later use
3067:    of a nonlinear solver.

3069:    Collective on snes

3071:    Input Parameters:
3072: .  snes - the `SNES` context

3074:    Note:
3075:    For basic use of the `SNES` solvers the user need not explicitly call
3076:    `SNESSetUp()`, since these actions will automatically occur during
3077:    the call to `SNESSolve()`.  However, if one wishes to control this
3078:    phase separately, `SNESSetUp()` should be called after `SNESCreate()`
3079:    and optional routines of the form SNESSetXXX(), but before `SNESSolve()`.

3081:    Level: advanced

3083: .seealso: `SNES`, `SNESCreate()`, `SNESSolve()`, `SNESDestroy()`
3084: @*/
3085: PetscErrorCode SNESSetUp(SNES snes)
3086: {
3087:   DM             dm;
3088:   DMSNES         sdm;
3089:   SNESLineSearch linesearch, pclinesearch;
3090:   void          *lsprectx, *lspostctx;
3091:   PetscErrorCode (*precheck)(SNESLineSearch, Vec, Vec, PetscBool *, void *);
3092:   PetscErrorCode (*postcheck)(SNESLineSearch, Vec, Vec, Vec, PetscBool *, PetscBool *, void *);
3093:   PetscErrorCode (*func)(SNES, Vec, Vec, void *);
3094:   Vec   f, fpc;
3095:   void *funcctx;
3096:   PetscErrorCode (*jac)(SNES, Vec, Mat, Mat, void *);
3097:   void *jacctx, *appctx;
3098:   Mat   j, jpre;

3101:   if (snes->setupcalled) return 0;
3102:   PetscLogEventBegin(SNES_Setup, snes, 0, 0, 0);

3104:   if (!((PetscObject)snes)->type_name) SNESSetType(snes, SNESNEWTONLS);

3106:   SNESGetFunction(snes, &snes->vec_func, NULL, NULL);

3108:   SNESGetDM(snes, &dm);
3109:   DMGetDMSNES(dm, &sdm);
3110:   SNESSetDefaultComputeJacobian(snes);

3112:   if (!snes->vec_func) DMCreateGlobalVector(dm, &snes->vec_func);

3114:   if (!snes->ksp) SNESGetKSP(snes, &snes->ksp);

3116:   if (snes->linesearch) {
3117:     SNESGetLineSearch(snes, &snes->linesearch);
3118:     SNESLineSearchSetFunction(snes->linesearch, SNESComputeFunction);
3119:   }

3121:   if (snes->npc && snes->npcside == PC_LEFT) {
3122:     snes->mf          = PETSC_TRUE;
3123:     snes->mf_operator = PETSC_FALSE;
3124:   }

3126:   if (snes->npc) {
3127:     /* copy the DM over */
3128:     SNESGetDM(snes, &dm);
3129:     SNESSetDM(snes->npc, dm);

3131:     SNESGetFunction(snes, &f, &func, &funcctx);
3132:     VecDuplicate(f, &fpc);
3133:     SNESSetFunction(snes->npc, fpc, func, funcctx);
3134:     SNESGetJacobian(snes, &j, &jpre, &jac, &jacctx);
3135:     SNESSetJacobian(snes->npc, j, jpre, jac, jacctx);
3136:     SNESGetApplicationContext(snes, &appctx);
3137:     SNESSetApplicationContext(snes->npc, appctx);
3138:     VecDestroy(&fpc);

3140:     /* copy the function pointers over */
3141:     PetscObjectCopyFortranFunctionPointers((PetscObject)snes, (PetscObject)snes->npc);

3143:     /* default to 1 iteration */
3144:     SNESSetTolerances(snes->npc, 0.0, 0.0, 0.0, 1, snes->npc->max_funcs);
3145:     if (snes->npcside == PC_RIGHT) {
3146:       SNESSetNormSchedule(snes->npc, SNES_NORM_FINAL_ONLY);
3147:     } else {
3148:       SNESSetNormSchedule(snes->npc, SNES_NORM_NONE);
3149:     }
3150:     SNESSetFromOptions(snes->npc);

3152:     /* copy the line search context over */
3153:     if (snes->linesearch && snes->npc->linesearch) {
3154:       SNESGetLineSearch(snes, &linesearch);
3155:       SNESGetLineSearch(snes->npc, &pclinesearch);
3156:       SNESLineSearchGetPreCheck(linesearch, &precheck, &lsprectx);
3157:       SNESLineSearchGetPostCheck(linesearch, &postcheck, &lspostctx);
3158:       SNESLineSearchSetPreCheck(pclinesearch, precheck, lsprectx);
3159:       SNESLineSearchSetPostCheck(pclinesearch, postcheck, lspostctx);
3160:       PetscObjectCopyFortranFunctionPointers((PetscObject)linesearch, (PetscObject)pclinesearch);
3161:     }
3162:   }
3163:   if (snes->mf) SNESSetUpMatrixFree_Private(snes, snes->mf_operator, snes->mf_version);
3164:   if (snes->ops->usercompute && !snes->user) (*snes->ops->usercompute)(snes, (void **)&snes->user);

3166:   snes->jac_iter = 0;
3167:   snes->pre_iter = 0;

3169:   PetscTryTypeMethod(snes, setup);

3171:   SNESSetDefaultComputeJacobian(snes);

3173:   if (snes->npc && snes->npcside == PC_LEFT) {
3174:     if (snes->functype == SNES_FUNCTION_PRECONDITIONED) {
3175:       if (snes->linesearch) {
3176:         SNESGetLineSearch(snes, &linesearch);
3177:         SNESLineSearchSetFunction(linesearch, SNESComputeFunctionDefaultNPC);
3178:       }
3179:     }
3180:   }
3181:   PetscLogEventEnd(SNES_Setup, snes, 0, 0, 0);
3182:   snes->setupcalled = PETSC_TRUE;
3183:   return 0;
3184: }

3186: /*@
3187:    SNESReset - Resets a `SNES` context to the snessetupcalled = 0 state and removes any allocated `Vec`s and `Mat`s

3189:    Collective on snes

3191:    Input Parameter:
3192: .  snes - iterative context obtained from `SNESCreate()`

3194:    Level: intermediate

3196:    Notes:
3197:    Call this if you wish to reuse a `SNES` but with different size vectors

3199:    Also calls the application context destroy routine set with `SNESSetComputeApplicationContext()`

3201: .seealso: `SNES`, `SNESDestroy()`, `SNESCreate()`, `SNESSetUp()`, `SNESSolve()`
3202: @*/
3203: PetscErrorCode SNESReset(SNES snes)
3204: {
3206:   if (snes->ops->userdestroy && snes->user) {
3207:     (*snes->ops->userdestroy)((void **)&snes->user);
3208:     snes->user = NULL;
3209:   }
3210:   if (snes->npc) SNESReset(snes->npc);

3212:   PetscTryTypeMethod(snes, reset);
3213:   if (snes->ksp) KSPReset(snes->ksp);

3215:   if (snes->linesearch) SNESLineSearchReset(snes->linesearch);

3217:   VecDestroy(&snes->vec_rhs);
3218:   VecDestroy(&snes->vec_sol);
3219:   VecDestroy(&snes->vec_sol_update);
3220:   VecDestroy(&snes->vec_func);
3221:   MatDestroy(&snes->jacobian);
3222:   MatDestroy(&snes->jacobian_pre);
3223:   MatDestroy(&snes->picard);
3224:   VecDestroyVecs(snes->nwork, &snes->work);
3225:   VecDestroyVecs(snes->nvwork, &snes->vwork);

3227:   snes->alwayscomputesfinalresidual = PETSC_FALSE;

3229:   snes->nwork = snes->nvwork = 0;
3230:   snes->setupcalled          = PETSC_FALSE;
3231:   return 0;
3232: }

3234: /*@
3235:    SNESConvergedReasonViewCancel - Clears all the reason view functions for a `SNES` object.

3237:    Collective on snes

3239:    Input Parameter:
3240: .  snes - iterative context obtained from `SNESCreate()`

3242:    Level: intermediate

3244: .seealso: `SNES`, `SNESCreate()`, `SNESDestroy()`, `SNESReset()`
3245: @*/
3246: PetscErrorCode SNESConvergedReasonViewCancel(SNES snes)
3247: {
3248:   PetscInt i;

3251:   for (i = 0; i < snes->numberreasonviews; i++) {
3252:     if (snes->reasonviewdestroy[i]) (*snes->reasonviewdestroy[i])(&snes->reasonviewcontext[i]);
3253:   }
3254:   snes->numberreasonviews = 0;
3255:   return 0;
3256: }

3258: /*@C
3259:    SNESDestroy - Destroys the nonlinear solver context that was created
3260:    with `SNESCreate()`.

3262:    Collective on snes

3264:    Input Parameter:
3265: .  snes - the `SNES` context

3267:    Level: beginner

3269: .seealso: `SNES`, `SNESCreate()`, `SNESSolve()`
3270: @*/
3271: PetscErrorCode SNESDestroy(SNES *snes)
3272: {
3273:   if (!*snes) return 0;
3275:   if (--((PetscObject)(*snes))->refct > 0) {
3276:     *snes = NULL;
3277:     return 0;
3278:   }

3280:   SNESReset((*snes));
3281:   SNESDestroy(&(*snes)->npc);

3283:   /* if memory was published with SAWs then destroy it */
3284:   PetscObjectSAWsViewOff((PetscObject)*snes);
3285:   PetscTryTypeMethod((*snes), destroy);

3287:   if ((*snes)->dm) DMCoarsenHookRemove((*snes)->dm, DMCoarsenHook_SNESVecSol, DMRestrictHook_SNESVecSol, *snes);
3288:   DMDestroy(&(*snes)->dm);
3289:   KSPDestroy(&(*snes)->ksp);
3290:   SNESLineSearchDestroy(&(*snes)->linesearch);

3292:   PetscFree((*snes)->kspconvctx);
3293:   if ((*snes)->ops->convergeddestroy) (*(*snes)->ops->convergeddestroy)((*snes)->cnvP);
3294:   if ((*snes)->conv_hist_alloc) PetscFree2((*snes)->conv_hist, (*snes)->conv_hist_its);
3295:   SNESMonitorCancel((*snes));
3296:   SNESConvergedReasonViewCancel((*snes));
3297:   PetscHeaderDestroy(snes);
3298:   return 0;
3299: }

3301: /* ----------- Routines to set solver parameters ---------- */

3303: /*@
3304:    SNESSetLagPreconditioner - Determines when the preconditioner is rebuilt in the nonlinear solve.

3306:    Logically Collective on snes

3308:    Input Parameters:
3309: +  snes - the `SNES` context
3310: -  lag - 1 means rebuild every time the Jacobian is computed within a single nonlinear solve, 2 means every second time
3311:          the Jacobian is built etc. -2 indicates rebuild preconditioner at next chance but then never rebuild after that

3313:    Options Database Keys:
3314: +    -snes_lag_jacobian_persists <true,false> - sets the persistence
3315: .    -snes_lag_jacobian <-2,1,2,...> - sets the lag
3316: .    -snes_lag_preconditioner_persists <true,false> - sets the persistence
3317: -    -snes_lag_preconditioner <-2,1,2,...> - sets the lag

3319:    Notes:
3320:    The default is 1
3321:    The preconditioner is ALWAYS built in the first iteration of a nonlinear solve unless lag is -1 or `SNESSetLagPreconditionerPersists()` was called

3323:    `SNESSetLagPreconditionerPersists()` allows using the same uniform lagging (for example every second linear solve) across multiple nonlinear solves.

3325:    Level: intermediate

3327: .seealso: `SNESSetTrustRegionTolerance()`, `SNESGetLagPreconditioner()`, `SNESSetLagJacobian()`, `SNESGetLagJacobian()`, `SNESSetLagPreconditionerPersists()`,
3328:           `SNESSetLagJacobianPersists()`, `SNES`, `SNESSolve()`
3329: @*/
3330: PetscErrorCode SNESSetLagPreconditioner(SNES snes, PetscInt lag)
3331: {
3336:   snes->lagpreconditioner = lag;
3337:   return 0;
3338: }

3340: /*@
3341:    SNESSetGridSequence - sets the number of steps of grid sequencing that `SNES` will do

3343:    Logically Collective on snes

3345:    Input Parameters:
3346: +  snes - the `SNES` context
3347: -  steps - the number of refinements to do, defaults to 0

3349:    Options Database Key:
3350: .    -snes_grid_sequence <steps> - Use grid sequencing to generate initial guess

3352:    Level: intermediate

3354:    Note:
3355:    Use `SNESGetSolution()` to extract the fine grid solution after grid sequencing.

3357: .seealso: `SNES`, `SNESSetTrustRegionTolerance()`, `SNESGetLagPreconditioner()`, `SNESSetLagJacobian()`, `SNESGetLagJacobian()`, `SNESGetGridSequence()`
3358: @*/
3359: PetscErrorCode SNESSetGridSequence(SNES snes, PetscInt steps)
3360: {
3363:   snes->gridsequence = steps;
3364:   return 0;
3365: }

3367: /*@
3368:    SNESGetGridSequence - gets the number of steps of grid sequencing that `SNES` will do

3370:    Logically Collective on snes

3372:    Input Parameter:
3373: .  snes - the `SNES` context

3375:    Output Parameter:
3376: .  steps - the number of refinements to do, defaults to 0

3378:    Options Database Key:
3379: .    -snes_grid_sequence <steps> - set number of refinements

3381:    Level: intermediate

3383:    Note:
3384:    Use `SNESGetSolution()` to extract the fine grid solution after grid sequencing.

3386: .seealso: `SNESSetTrustRegionTolerance()`, `SNESGetLagPreconditioner()`, `SNESSetLagJacobian()`, `SNESGetLagJacobian()`, `SNESSetGridSequence()`
3387: @*/
3388: PetscErrorCode SNESGetGridSequence(SNES snes, PetscInt *steps)
3389: {
3391:   *steps = snes->gridsequence;
3392:   return 0;
3393: }

3395: /*@
3396:    SNESGetLagPreconditioner - Return how often the preconditioner is rebuilt

3398:    Not Collective

3400:    Input Parameter:
3401: .  snes - the `SNES` context

3403:    Output Parameter:
3404: .   lag - -1 indicates NEVER rebuild, 1 means rebuild every time the Jacobian is computed within a single nonlinear solve, 2 means every second time
3405:          the Jacobian is built etc. -2 indicates rebuild preconditioner at next chance but then never rebuild after that

3407:    Options Database Keys:
3408: +    -snes_lag_jacobian_persists <true,false> - sets the persistence
3409: .    -snes_lag_jacobian <-2,1,2,...> - sets the lag
3410: .    -snes_lag_preconditioner_persists <true,false> - sets the persistence
3411: -    -snes_lag_preconditioner <-2,1,2,...> - sets the lag

3413:    Notes:
3414:    The default is 1

3416:    The preconditioner is ALWAYS built in the first iteration of a nonlinear solve unless lag is -1

3418:    Level: intermediate

3420: .seealso: `SNES`, `SNESSetTrustRegionTolerance()`, `SNESSetLagPreconditioner()`, `SNESSetLagJacobianPersists()`, `SNESSetLagPreconditionerPersists()`
3421: @*/
3422: PetscErrorCode SNESGetLagPreconditioner(SNES snes, PetscInt *lag)
3423: {
3425:   *lag = snes->lagpreconditioner;
3426:   return 0;
3427: }

3429: /*@
3430:    SNESSetLagJacobian - Set when the Jacobian is rebuilt in the nonlinear solve. See `SNESSetLagPreconditioner()` for determining how
3431:      often the preconditioner is rebuilt.

3433:    Logically Collective on snes

3435:    Input Parameters:
3436: +  snes - the `SNES` context
3437: -  lag - -1 indicates NEVER rebuild, 1 means rebuild every time the Jacobian is computed within a single nonlinear solve, 2 means every second time
3438:          the Jacobian is built etc. -2 means rebuild at next chance but then never again

3440:    Options Database Keys:
3441: +    -snes_lag_jacobian_persists <true,false> - sets the persistence
3442: .    -snes_lag_jacobian <-2,1,2,...> - sets the lag
3443: .    -snes_lag_preconditioner_persists <true,false> - sets the persistence
3444: -    -snes_lag_preconditioner <-2,1,2,...> - sets the lag.

3446:    Notes:
3447:    The default is 1

3449:    The Jacobian is ALWAYS built in the first iteration of a nonlinear solve unless lag is -1

3451:    If  -1 is used before the very first nonlinear solve the CODE WILL FAIL! because no Jacobian is used, use -2 to indicate you want it recomputed
3452:    at the next Newton step but never again (unless it is reset to another value)

3454:    Level: intermediate

3456: .seealso: `SNES`, `SNESSetTrustRegionTolerance()`, `SNESGetLagPreconditioner()`, `SNESSetLagPreconditioner()`, `SNESGetLagJacobianPersists()`, `SNESSetLagPreconditionerPersists()`
3457: @*/
3458: PetscErrorCode SNESSetLagJacobian(SNES snes, PetscInt lag)
3459: {
3464:   snes->lagjacobian = lag;
3465:   return 0;
3466: }

3468: /*@
3469:    SNESGetLagJacobian - Get how often the Jacobian is rebuilt. See `SNESGetLagPreconditioner()` to determine when the preconditioner is rebuilt

3471:    Not Collective

3473:    Input Parameter:
3474: .  snes - the `SNES` context

3476:    Output Parameter:
3477: .   lag - -1 indicates NEVER rebuild, 1 means rebuild every time the Jacobian is computed within a single nonlinear solve, 2 means every second time
3478:          the Jacobian is built etc.

3480:    Notes:
3481:    The default is 1

3483:    The jacobian is ALWAYS built in the first iteration of a nonlinear solve unless lag is -1 or `SNESSetLagJacobianPersists()` was called.

3485:    Level: intermediate

3487: .seealso: `SNES`, `SNESSetTrustRegionTolerance()`, `SNESSetLagJacobian()`, `SNESSetLagPreconditioner()`, `SNESGetLagPreconditioner()`, `SNESSetLagJacobianPersists()`, `SNESSetLagPreconditionerPersists()`

3489: @*/
3490: PetscErrorCode SNESGetLagJacobian(SNES snes, PetscInt *lag)
3491: {
3493:   *lag = snes->lagjacobian;
3494:   return 0;
3495: }

3497: /*@
3498:    SNESSetLagJacobianPersists - Set whether or not the Jacobian lagging persists through multiple nonlinear solves

3500:    Logically collective on snes

3502:    Input Parameters:
3503: +  snes - the `SNES` context
3504: -   flg - jacobian lagging persists if true

3506:    Options Database Keys:
3507: +    -snes_lag_jacobian_persists <true,false> - sets the persistence
3508: .    -snes_lag_jacobian <-2,1,2,...> - sets the lag
3509: .    -snes_lag_preconditioner_persists <true,false> - sets the persistence
3510: -    -snes_lag_preconditioner <-2,1,2,...> - sets the lag

3512:    Notes:
3513:     Normally when `SNESSetLagJacobian()` is used, the Jacobian is always rebuilt at the beginning of each new nonlinear solve, this removes that.

3515:     This is useful both for nonlinear preconditioning, where it's appropriate to have the Jacobian be stale by
3516:    several solves, and for implicit time-stepping, where Jacobian lagging in the inner nonlinear solve over several
3517:    timesteps may present huge efficiency gains.

3519:    Level: advanced

3521: .seealso: `SNES, `SNESSetLagPreconditionerPersists()`, `SNESSetLagJacobian()`, `SNESGetLagJacobian()`, `SNESGetNPC()`, `SNESSetLagJacobianPersists()`
3522: @*/
3523: PetscErrorCode SNESSetLagJacobianPersists(SNES snes, PetscBool flg)
3524: {
3527:   snes->lagjac_persist = flg;
3528:   return 0;
3529: }

3531: /*@
3532:    SNESSetLagPreconditionerPersists - Set whether or not the preconditioner lagging persists through multiple nonlinear solves

3534:    Logically Collective on snes

3536:    Input Parameters:
3537: +  snes - the `SNES` context
3538: -   flg - preconditioner lagging persists if true

3540:    Options Database Keys:
3541: +    -snes_lag_jacobian_persists <true,false> - sets the persistence
3542: .    -snes_lag_jacobian <-2,1,2,...> - sets the lag
3543: .    -snes_lag_preconditioner_persists <true,false> - sets the persistence
3544: -    -snes_lag_preconditioner <-2,1,2,...> - sets the lag

3546:    Notes:
3547:     Normally when `SNESSetLagPreconditioner()` is used, the preconditioner is always rebuilt at the beginning of each new nonlinear solve, this removes that.

3549:    This is useful both for nonlinear preconditioning, where it's appropriate to have the preconditioner be stale
3550:    by several solves, and for implicit time-stepping, where preconditioner lagging in the inner nonlinear solve over
3551:    several timesteps may present huge efficiency gains.

3553:    Level: developer

3555: .seealso: `SNES`, `SNESSetLagJacobianPersists()`, `SNESSetLagJacobian()`, `SNESGetLagJacobian()`, `SNESGetNPC()`, `SNESSetLagPreconditioner()`
3556: @*/
3557: PetscErrorCode SNESSetLagPreconditionerPersists(SNES snes, PetscBool flg)
3558: {
3561:   snes->lagpre_persist = flg;
3562:   return 0;
3563: }

3565: /*@
3566:    SNESSetForceIteration - force `SNESSolve()` to take at least one iteration regardless of the initial residual norm

3568:    Logically Collective on snes

3570:    Input Parameters:
3571: +  snes - the `SNES` context
3572: -  force - `PETSC_TRUE` require at least one iteration

3574:    Options Database Key:
3575: .    -snes_force_iteration <force> - Sets forcing an iteration

3577:    Note:
3578:    This is used sometimes with `TS` to prevent `TS` from detecting a false steady state solution

3580:    Level: intermediate

3582: .seealso: `SNES`, `TS`, `SNESSetTrustRegionTolerance()`, `SNESSetDivergenceTolerance()`
3583: @*/
3584: PetscErrorCode SNESSetForceIteration(SNES snes, PetscBool force)
3585: {
3587:   snes->forceiteration = force;
3588:   return 0;
3589: }

3591: /*@
3592:    SNESGetForceIteration - Check whether or not `SNESSolve()` take at least one iteration regardless of the initial residual norm

3594:    Logically Collective on snes

3596:    Input Parameters:
3597: .  snes - the `SNES` context

3599:    Output Parameter:
3600: .  force - PETSC_TRUE requires at least one iteration.

3602:    Level: intermediate

3604: .seealso: `SNES`, `SNESSetForceIteration()`, `SNESSetTrustRegionTolerance()`, `SNESSetDivergenceTolerance()`
3605: @*/
3606: PetscErrorCode SNESGetForceIteration(SNES snes, PetscBool *force)
3607: {
3609:   *force = snes->forceiteration;
3610:   return 0;
3611: }

3613: /*@
3614:    SNESSetTolerances - Sets `SNES` various parameters used in convergence tests.

3616:    Logically Collective on snes

3618:    Input Parameters:
3619: +  snes - the `SNES` context
3620: .  abstol - absolute convergence tolerance
3621: .  rtol - relative convergence tolerance
3622: .  stol -  convergence tolerance in terms of the norm of the change in the solution between steps,  || delta x || < stol*|| x ||
3623: .  maxit - maximum number of iterations, default 50.
3624: -  maxf - maximum number of function evaluations (-1 indicates no limit), default 1000

3626:    Options Database Keys:
3627: +    -snes_atol <abstol> - Sets abstol
3628: .    -snes_rtol <rtol> - Sets rtol
3629: .    -snes_stol <stol> - Sets stol
3630: .    -snes_max_it <maxit> - Sets maxit
3631: -    -snes_max_funcs <maxf> - Sets maxf

3633:    Level: intermediate

3635: .seealso: `SNESolve()`, `SNES`, `SNESSetTrustRegionTolerance()`, `SNESSetDivergenceTolerance()`, `SNESSetForceIteration()`
3636: @*/
3637: PetscErrorCode SNESSetTolerances(SNES snes, PetscReal abstol, PetscReal rtol, PetscReal stol, PetscInt maxit, PetscInt maxf)
3638: {

3646:   if (abstol != PETSC_DEFAULT) {
3648:     snes->abstol = abstol;
3649:   }
3650:   if (rtol != PETSC_DEFAULT) {
3652:     snes->rtol = rtol;
3653:   }
3654:   if (stol != PETSC_DEFAULT) {
3656:     snes->stol = stol;
3657:   }
3658:   if (maxit != PETSC_DEFAULT) {
3660:     snes->max_its = maxit;
3661:   }
3662:   if (maxf != PETSC_DEFAULT) {
3664:     snes->max_funcs = maxf;
3665:   }
3666:   snes->tolerancesset = PETSC_TRUE;
3667:   return 0;
3668: }

3670: /*@
3671:    SNESSetDivergenceTolerance - Sets the divergence tolerance used for the `SNES` divergence test.

3673:    Logically Collective on snes

3675:    Input Parameters:
3676: +  snes - the `SNES` context
3677: -  divtol - the divergence tolerance. Use -1 to deactivate the test, default is 1e4

3679:    Options Database Key:
3680: .    -snes_divergence_tolerance <divtol> - Sets divtol

3682:    Level: intermediate

3684: .seealso: `SNES`, `SNESSolve()`, `SNESSetTolerances()`, `SNESGetDivergenceTolerance`
3685: @*/
3686: PetscErrorCode SNESSetDivergenceTolerance(SNES snes, PetscReal divtol)
3687: {

3691:   if (divtol != PETSC_DEFAULT) {
3692:     snes->divtol = divtol;
3693:   } else {
3694:     snes->divtol = 1.0e4;
3695:   }
3696:   return 0;
3697: }

3699: /*@
3700:    SNESGetTolerances - Gets various parameters used in convergence tests.

3702:    Not Collective

3704:    Input Parameters:
3705: +  snes - the `SNES` context
3706: .  atol - absolute convergence tolerance
3707: .  rtol - relative convergence tolerance
3708: .  stol -  convergence tolerance in terms of the norm
3709:            of the change in the solution between steps
3710: .  maxit - maximum number of iterations
3711: -  maxf - maximum number of function evaluations

3713:    Note:
3714:    The user can specify NULL for any parameter that is not needed.

3716:    Level: intermediate

3718: .seealso: `SNES`, `SNESSetTolerances()`
3719: @*/
3720: PetscErrorCode SNESGetTolerances(SNES snes, PetscReal *atol, PetscReal *rtol, PetscReal *stol, PetscInt *maxit, PetscInt *maxf)
3721: {
3723:   if (atol) *atol = snes->abstol;
3724:   if (rtol) *rtol = snes->rtol;
3725:   if (stol) *stol = snes->stol;
3726:   if (maxit) *maxit = snes->max_its;
3727:   if (maxf) *maxf = snes->max_funcs;
3728:   return 0;
3729: }

3731: /*@
3732:    SNESGetDivergenceTolerance - Gets divergence tolerance used in divergence test.

3734:    Not Collective

3736:    Input Parameters:
3737: +  snes - the `SNES` context
3738: -  divtol - divergence tolerance

3740:    Level: intermediate

3742: .seealso: `SNES`, `SNESSetDivergenceTolerance()`
3743: @*/
3744: PetscErrorCode SNESGetDivergenceTolerance(SNES snes, PetscReal *divtol)
3745: {
3747:   if (divtol) *divtol = snes->divtol;
3748:   return 0;
3749: }

3751: /*@
3752:    SNESSetTrustRegionTolerance - Sets the trust region parameter tolerance.

3754:    Logically Collective on snes

3756:    Input Parameters:
3757: +  snes - the `SNES` context
3758: -  tol - tolerance

3760:    Options Database Key:
3761: .  -snes_trtol <tol> - Sets tol

3763:    Level: intermediate

3765: .seealso: `SNES`, `SNESNEWTONTRDC`, `SNESSetTolerances()`
3766: @*/
3767: PetscErrorCode SNESSetTrustRegionTolerance(SNES snes, PetscReal tol)
3768: {
3771:   snes->deltatol = tol;
3772:   return 0;
3773: }

3775: PETSC_INTERN PetscErrorCode SNESMonitorRange_Private(SNES, PetscInt, PetscReal *);

3777: PetscErrorCode SNESMonitorLGRange(SNES snes, PetscInt n, PetscReal rnorm, void *monctx)
3778: {
3779:   PetscDrawLG      lg;
3780:   PetscReal        x, y, per;
3781:   PetscViewer      v = (PetscViewer)monctx;
3782:   static PetscReal prev; /* should be in the context */
3783:   PetscDraw        draw;

3786:   PetscViewerDrawGetDrawLG(v, 0, &lg);
3787:   if (!n) PetscDrawLGReset(lg);
3788:   PetscDrawLGGetDraw(lg, &draw);
3789:   PetscDrawSetTitle(draw, "Residual norm");
3790:   x = (PetscReal)n;
3791:   if (rnorm > 0.0) y = PetscLog10Real(rnorm);
3792:   else y = -15.0;
3793:   PetscDrawLGAddPoint(lg, &x, &y);
3794:   if (n < 20 || !(n % 5) || snes->reason) {
3795:     PetscDrawLGDraw(lg);
3796:     PetscDrawLGSave(lg);
3797:   }

3799:   PetscViewerDrawGetDrawLG(v, 1, &lg);
3800:   if (!n) PetscDrawLGReset(lg);
3801:   PetscDrawLGGetDraw(lg, &draw);
3802:   PetscDrawSetTitle(draw, "% elemts > .2*max elemt");
3803:   SNESMonitorRange_Private(snes, n, &per);
3804:   x = (PetscReal)n;
3805:   y = 100.0 * per;
3806:   PetscDrawLGAddPoint(lg, &x, &y);
3807:   if (n < 20 || !(n % 5) || snes->reason) {
3808:     PetscDrawLGDraw(lg);
3809:     PetscDrawLGSave(lg);
3810:   }

3812:   PetscViewerDrawGetDrawLG(v, 2, &lg);
3813:   if (!n) {
3814:     prev = rnorm;
3815:     PetscDrawLGReset(lg);
3816:   }
3817:   PetscDrawLGGetDraw(lg, &draw);
3818:   PetscDrawSetTitle(draw, "(norm -oldnorm)/oldnorm");
3819:   x = (PetscReal)n;
3820:   y = (prev - rnorm) / prev;
3821:   PetscDrawLGAddPoint(lg, &x, &y);
3822:   if (n < 20 || !(n % 5) || snes->reason) {
3823:     PetscDrawLGDraw(lg);
3824:     PetscDrawLGSave(lg);
3825:   }

3827:   PetscViewerDrawGetDrawLG(v, 3, &lg);
3828:   if (!n) PetscDrawLGReset(lg);
3829:   PetscDrawLGGetDraw(lg, &draw);
3830:   PetscDrawSetTitle(draw, "(norm -oldnorm)/oldnorm*(% > .2 max)");
3831:   x = (PetscReal)n;
3832:   y = (prev - rnorm) / (prev * per);
3833:   if (n > 2) { /*skip initial crazy value */
3834:     PetscDrawLGAddPoint(lg, &x, &y);
3835:   }
3836:   if (n < 20 || !(n % 5) || snes->reason) {
3837:     PetscDrawLGDraw(lg);
3838:     PetscDrawLGSave(lg);
3839:   }
3840:   prev = rnorm;
3841:   return 0;
3842: }

3844: /*@
3845:    SNESMonitor - runs the user provided monitor routines, if they exist

3847:    Collective on snes

3849:    Input Parameters:
3850: +  snes - nonlinear solver context obtained from `SNESCreate()`
3851: .  iter - iteration number
3852: -  rnorm - relative norm of the residual

3854:    Note:
3855:    This routine is called by the `SNES` implementations.
3856:    It does not typically need to be called by the user.

3858:    Level: developer

3860: .seealso: `SNES`, `SNESMonitorSet()`
3861: @*/
3862: PetscErrorCode SNESMonitor(SNES snes, PetscInt iter, PetscReal rnorm)
3863: {
3864:   PetscInt i, n = snes->numbermonitors;

3866:   VecLockReadPush(snes->vec_sol);
3867:   for (i = 0; i < n; i++) (*snes->monitor[i])(snes, iter, rnorm, snes->monitorcontext[i]);
3868:   VecLockReadPop(snes->vec_sol);
3869:   return 0;
3870: }

3872: /* ------------ Routines to set performance monitoring options ----------- */

3874: /*MC
3875:     SNESMonitorFunction - functional form passed to `SNESMonitorSet()` to monitor convergence of nonlinear solver

3877:      Synopsis:
3878: #include <petscsnes.h>
3879: $    PetscErrorCode SNESMonitorFunction(SNES snes,PetscInt its, PetscReal norm,void *mctx)

3881:      Collective on snes

3883:     Input Parameters:
3884: +    snes - the `SNES` context
3885: .    its - iteration number
3886: .    norm - 2-norm function value (may be estimated)
3887: -    mctx - [optional] monitoring context

3889:    Level: advanced

3891: .seealso: `SNESMonitorSet()`, `SNESMonitorSet()`, `SNESMonitorGet()`
3892: M*/

3894: /*@C
3895:    SNESMonitorSet - Sets an ADDITIONAL function that is to be used at every
3896:    iteration of the nonlinear solver to display the iteration's
3897:    progress.

3899:    Logically Collective on snes

3901:    Input Parameters:
3902: +  snes - the `SNES` context
3903: .  f - the monitor function, see `SNESMonitorFunction` for the calling sequence
3904: .  mctx - [optional] user-defined context for private data for the
3905:           monitor routine (use NULL if no context is desired)
3906: -  monitordestroy - [optional] routine that frees monitor context
3907:           (may be NULL)

3909:    Options Database Keys:
3910: +    -snes_monitor        - sets `SNESMonitorDefault()`
3911: .    -snes_monitor draw::draw_lg - sets line graph monitor,
3912: -    -snes_monitor_cancel - cancels all monitors that have
3913:                             been hardwired into a code by
3914:                             calls to SNESMonitorSet(), but
3915:                             does not cancel those set via
3916:                             the options database.

3918:    Note:
3919:    Several different monitoring routines may be set by calling
3920:    `SNESMonitorSet()` multiple times; all will be called in the
3921:    order in which they were set.

3923:    Fortran Note:
3924:    Only a single monitor function can be set for each `SNES` object

3926:    Level: intermediate

3928: .seealso: `SNES`, `SNESSolve()`, `SNESMonitorDefault()`, `SNESMonitorCancel()`, `SNESMonitorFunction`
3929: @*/
3930: PetscErrorCode SNESMonitorSet(SNES snes, PetscErrorCode (*f)(SNES, PetscInt, PetscReal, void *), void *mctx, PetscErrorCode (*monitordestroy)(void **))
3931: {
3932:   PetscInt  i;
3933:   PetscBool identical;

3936:   for (i = 0; i < snes->numbermonitors; i++) {
3937:     PetscMonitorCompare((PetscErrorCode(*)(void))f, mctx, monitordestroy, (PetscErrorCode(*)(void))snes->monitor[i], snes->monitorcontext[i], snes->monitordestroy[i], &identical);
3938:     if (identical) return 0;
3939:   }
3941:   snes->monitor[snes->numbermonitors]          = f;
3942:   snes->monitordestroy[snes->numbermonitors]   = monitordestroy;
3943:   snes->monitorcontext[snes->numbermonitors++] = (void *)mctx;
3944:   return 0;
3945: }

3947: /*@
3948:    SNESMonitorCancel - Clears all the monitor functions for a `SNES` object.

3950:    Logically Collective on snes

3952:    Input Parameters:
3953: .  snes - the `SNES` context

3955:    Options Database Key:
3956: .  -snes_monitor_cancel - cancels all monitors that have been hardwired
3957:     into a code by calls to SNESMonitorSet(), but does not cancel those
3958:     set via the options database

3960:    Note:
3961:    There is no way to clear one specific monitor from a `SNES` object.

3963:    Level: intermediate

3965: .seealso: `SNES`, `SNESMonitorGet()`, `SNESMonitorDefault()`, `SNESMonitorSet()`
3966: @*/
3967: PetscErrorCode SNESMonitorCancel(SNES snes)
3968: {
3969:   PetscInt i;

3972:   for (i = 0; i < snes->numbermonitors; i++) {
3973:     if (snes->monitordestroy[i]) (*snes->monitordestroy[i])(&snes->monitorcontext[i]);
3974:   }
3975:   snes->numbermonitors = 0;
3976:   return 0;
3977: }

3979: /*MC
3980:     SNESConvergenceTestFunction - functional form used for testing of convergence of nonlinear solver

3982:      Synopsis:
3983: #include <petscsnes.h>
3984: $     PetscErrorCode SNESConvergenceTest(SNES snes,PetscInt it,PetscReal xnorm,PetscReal gnorm,PetscReal f,SNESConvergedReason *reason,void *cctx)

3986:      Collective on snes

3988:     Input Parameters:
3989: +    snes - the `SNES` context
3990: .    it - current iteration (0 is the first and is before any Newton step)
3991: .    xnorm - 2-norm of current iterate
3992: .    gnorm - 2-norm of current step
3993: .    f - 2-norm of function
3994: -    cctx - [optional] convergence context

3996:     Output Parameter:
3997: .    reason - reason for convergence/divergence, only needs to be set when convergence or divergence is detected

3999:    Level: intermediate

4001: .seealso: `SNES`, `SNESSolve`, `SNESSetConvergenceTest()`, `SNESGetConvergenceTest()`
4002: M*/

4004: /*@C
4005:    SNESSetConvergenceTest - Sets the function that is to be used
4006:    to test for convergence of the nonlinear iterative solution.

4008:    Logically Collective on snes

4010:    Input Parameters:
4011: +  snes - the `SNES` context
4012: .  `SNESConvergenceTestFunction` - routine to test for convergence
4013: .  cctx - [optional] context for private data for the convergence routine  (may be NULL)
4014: -  destroy - [optional] destructor for the context (may be NULL; PETSC_NULL_FUNCTION in Fortran)

4016:    Level: advanced

4018: .seealso: `SNES`, `SNESConvergedDefault()`, `SNESConvergedSkip()`, `SNESConvergenceTestFunction`
4019: @*/
4020: PetscErrorCode SNESSetConvergenceTest(SNES snes, PetscErrorCode (*SNESConvergenceTestFunction)(SNES, PetscInt, PetscReal, PetscReal, PetscReal, SNESConvergedReason *, void *), void *cctx, PetscErrorCode (*destroy)(void *))
4021: {
4023:   if (!SNESConvergenceTestFunction) SNESConvergenceTestFunction = SNESConvergedSkip;
4024:   if (snes->ops->convergeddestroy) (*snes->ops->convergeddestroy)(snes->cnvP);
4025:   snes->ops->converged        = SNESConvergenceTestFunction;
4026:   snes->ops->convergeddestroy = destroy;
4027:   snes->cnvP                  = cctx;
4028:   return 0;
4029: }

4031: /*@
4032:    SNESGetConvergedReason - Gets the reason the `SNES` iteration was stopped.

4034:    Not Collective

4036:    Input Parameter:
4037: .  snes - the `SNES` context

4039:    Output Parameter:
4040: .  reason - negative value indicates diverged, positive value converged, see `SNESConvergedReason` for the individual convergence tests for complete lists

4042:    Options Database Key:
4043: .   -snes_converged_reason - prints the reason to standard out

4045:    Level: intermediate

4047:    Note:
4048:     Should only be called after the call the `SNESSolve()` is complete, if it is called earlier it returns the value `SNES__CONVERGED_ITERATING`.

4050: .seealso: `SNESSolve()`, `SNESSetConvergenceTest()`, `SNESSetConvergedReason()`, `SNESConvergedReason`, `SNESGetConvergedReasonString()`
4051: @*/
4052: PetscErrorCode SNESGetConvergedReason(SNES snes, SNESConvergedReason *reason)
4053: {
4056:   *reason = snes->reason;
4057:   return 0;
4058: }

4060: /*@C
4061:    SNESGetConvergedReasonString - Return a human readable string for `SNESConvergedReason`

4063:    Not Collective

4065:    Input Parameter:
4066: .  snes - the `SNES` context

4068:    Output Parameter:
4069: .  strreason - a human readable string that describes SNES converged reason

4071:    Level: beginner

4073: .seealso: `SNES`, `SNESGetConvergedReason()`
4074: @*/
4075: PetscErrorCode SNESGetConvergedReasonString(SNES snes, const char **strreason)
4076: {
4079:   *strreason = SNESConvergedReasons[snes->reason];
4080:   return 0;
4081: }

4083: /*@
4084:    SNESSetConvergedReason - Sets the reason the `SNES` iteration was stopped.

4086:    Not Collective

4088:    Input Parameters:
4089: +  snes - the `SNES` context
4090: -  reason - negative value indicates diverged, positive value converged, see `SNESConvergedReason` or the
4091:             manual pages for the individual convergence tests for complete lists

4093:    Level: developer

4095:    Developer Note:
4096:    Called inside the various `SNESSolve()` implementations

4098: .seealso: `SNESGetConvergedReason()`, `SNESSetConvergenceTest()`, `SNESConvergedReason`
4099: @*/
4100: PetscErrorCode SNESSetConvergedReason(SNES snes, SNESConvergedReason reason)
4101: {
4103:   snes->reason = reason;
4104:   return 0;
4105: }

4107: /*@
4108:    SNESSetConvergenceHistory - Sets the array used to hold the convergence history.

4110:    Logically Collective on snes

4112:    Input Parameters:
4113: +  snes - iterative context obtained from `SNESCreate()`
4114: .  a   - array to hold history, this array will contain the function norms computed at each step
4115: .  its - integer array holds the number of linear iterations for each solve.
4116: .  na  - size of a and its
4117: -  reset - `PETSC_TRUE` indicates each new nonlinear solve resets the history counter to zero,
4118:            else it continues storing new values for new nonlinear solves after the old ones

4120:    Notes:
4121:    If 'a' and 'its' are NULL then space is allocated for the history. If 'na' `PETSC_DECIDE` or `PETSC_DEFAULT` then a
4122:    default array of length 10000 is allocated.

4124:    This routine is useful, e.g., when running a code for purposes
4125:    of accurate performance monitoring, when no I/O should be done
4126:    during the section of code that is being timed.

4128:    Level: intermediate

4130: .seealso: `SNES`, `SNESSolve()`, `SNESGetConvergenceHistory()`
4131: @*/
4132: PetscErrorCode SNESSetConvergenceHistory(SNES snes, PetscReal a[], PetscInt its[], PetscInt na, PetscBool reset)
4133: {
4137:   if (!a) {
4138:     if (na == PETSC_DECIDE || na == PETSC_DEFAULT) na = 1000;
4139:     PetscCalloc2(na, &a, na, &its);
4140:     snes->conv_hist_alloc = PETSC_TRUE;
4141:   }
4142:   snes->conv_hist       = a;
4143:   snes->conv_hist_its   = its;
4144:   snes->conv_hist_max   = (size_t)na;
4145:   snes->conv_hist_len   = 0;
4146:   snes->conv_hist_reset = reset;
4147:   return 0;
4148: }

4150: #if defined(PETSC_HAVE_MATLAB)
4151:   #include <engine.h> /* MATLAB include file */
4152:   #include <mex.h>    /* MATLAB include file */

4154: PETSC_EXTERN mxArray *SNESGetConvergenceHistoryMatlab(SNES snes)
4155: {
4156:   mxArray   *mat;
4157:   PetscInt   i;
4158:   PetscReal *ar;

4160:   mat = mxCreateDoubleMatrix(snes->conv_hist_len, 1, mxREAL);
4161:   ar  = (PetscReal *)mxGetData(mat);
4162:   for (i = 0; i < snes->conv_hist_len; i++) ar[i] = snes->conv_hist[i];
4163:   return mat;
4164: }
4165: #endif

4167: /*@C
4168:    SNESGetConvergenceHistory - Gets the array used to hold the convergence history.

4170:    Not Collective

4172:    Input Parameter:
4173: .  snes - iterative context obtained from `SNESCreate()`

4175:    Output Parameters:
4176: +  a   - array to hold history, usually was set with `SNESSetConvergenceHistory()`
4177: .  its - integer array holds the number of linear iterations (or
4178:          negative if not converged) for each solve.
4179: -  na  - size of a and its

4181:    Notes:
4182:     The calling sequence for this routine in Fortran is
4183: $   call SNESGetConvergenceHistory(SNES snes, integer na, integer ierr)

4185:    This routine is useful, e.g., when running a code for purposes
4186:    of accurate performance monitoring, when no I/O should be done
4187:    during the section of code that is being timed.

4189:    Level: intermediate

4191: .seealso: `SNES`, `SNESSolve()`, `SNESSetConvergenceHistory()`
4192: @*/
4193: PetscErrorCode SNESGetConvergenceHistory(SNES snes, PetscReal *a[], PetscInt *its[], PetscInt *na)
4194: {
4196:   if (a) *a = snes->conv_hist;
4197:   if (its) *its = snes->conv_hist_its;
4198:   if (na) *na = (PetscInt)snes->conv_hist_len;
4199:   return 0;
4200: }

4202: /*@C
4203:   SNESSetUpdate - Sets the general-purpose update function called
4204:   at the beginning of every iteration of the nonlinear solve. Specifically
4205:   it is called just before the Jacobian is "evaluated".

4207:   Logically Collective on snes

4209:   Input Parameters:
4210: + snes - The nonlinear solver context
4211: - func - The function

4213:   Calling sequence of func:
4214: $ func (SNES snes, PetscInt step);

4216: . step - The current step of the iteration

4218:   Level: advanced

4220:   Note:
4221:      This is NOT what one uses to update the ghost points before a function evaluation, that should be done at the beginning of your function provided
4222:      to `SNESSetFunction()`, or `SNESSetPicard()`
4223:      This is not used by most users.

4225:      There are a varity of function hooks one many set that are called at different stages of the nonlinear solution process, see the functions listed below.

4227: .seealso: `SNES`, `SNESSolve()`, `SNESSetJacobian()`, `SNESSolve()`, `SNESLineSearchSetPreCheck()`, `SNESLineSearchSetPostCheck()`, `SNESNewtonTRSetPreCheck()`, `SNESNewtonTRSetPostCheck()`,
4228:          `SNESMonitorSet()`, `SNESSetDivergenceTest()`
4229: @*/
4230: PetscErrorCode SNESSetUpdate(SNES snes, PetscErrorCode (*func)(SNES, PetscInt))
4231: {
4233:   snes->ops->update = func;
4234:   return 0;
4235: }

4237: /*
4238:    SNESScaleStep_Private - Scales a step so that its length is less than the
4239:    positive parameter delta.

4241:     Input Parameters:
4242: +   snes - the `SNES` context
4243: .   y - approximate solution of linear system
4244: .   fnorm - 2-norm of current function
4245: -   delta - trust region size

4247:     Output Parameters:
4248: +   gpnorm - predicted function norm at the new point, assuming local
4249:     linearization.  The value is zero if the step lies within the trust
4250:     region, and exceeds zero otherwise.
4251: -   ynorm - 2-norm of the step

4253:     Level: developer

4255:     Note:
4256:     For non-trust region methods such as `SNESNEWTONLS`, the parameter delta
4257:     is set to be the maximum allowable step size.
4258: */
4259: PetscErrorCode SNESScaleStep_Private(SNES snes, Vec y, PetscReal *fnorm, PetscReal *delta, PetscReal *gpnorm, PetscReal *ynorm)
4260: {
4261:   PetscReal   nrm;
4262:   PetscScalar cnorm;


4268:   VecNorm(y, NORM_2, &nrm);
4269:   if (nrm > *delta) {
4270:     nrm     = *delta / nrm;
4271:     *gpnorm = (1.0 - nrm) * (*fnorm);
4272:     cnorm   = nrm;
4273:     VecScale(y, cnorm);
4274:     *ynorm = *delta;
4275:   } else {
4276:     *gpnorm = 0.0;
4277:     *ynorm  = nrm;
4278:   }
4279:   return 0;
4280: }

4282: /*@C
4283:    SNESConvergedReasonView - Displays the reason a `SNES` solve converged or diverged to a viewer

4285:    Collective on snes

4287:    Parameter:
4288: +  snes - iterative context obtained from `SNESCreate()`
4289: -  viewer - the viewer to display the reason

4291:    Options Database Keys:
4292: +  -snes_converged_reason - print reason for converged or diverged, also prints number of iterations
4293: -  -snes_converged_reason ::failed - only print reason and number of iterations when diverged

4295:   Note:
4296:      To change the format of the output call `PetscViewerPushFormat`(viewer,format) before this call. Use `PETSC_VIEWER_DEFAULT` for the default,
4297:      use `PETSC_VIEWER_FAILED` to only display a reason if it fails.

4299:    Level: beginner

4301: .seealso: `SNESConvergedReason`, `PetscViewer`, `SNES`,
4302:           `SNESCreate()`, `SNESSetUp()`, `SNESDestroy()`, `SNESSetTolerances()`, `SNESConvergedDefault()`, `SNESGetConvergedReason()`,
4303:           `SNESConvergedReasonViewFromOptions()`,
4304:           `PetscViewerPushFormat()`, `PetscViewerPopFormat()`
4305: @*/
4306: PetscErrorCode SNESConvergedReasonView(SNES snes, PetscViewer viewer)
4307: {
4308:   PetscViewerFormat format;
4309:   PetscBool         isAscii;

4311:   if (!viewer) viewer = PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)snes));
4312:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isAscii);
4313:   if (isAscii) {
4314:     PetscViewerGetFormat(viewer, &format);
4315:     PetscViewerASCIIAddTab(viewer, ((PetscObject)snes)->tablevel);
4316:     if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
4317:       DM       dm;
4318:       Vec      u;
4319:       PetscDS  prob;
4320:       PetscInt Nf, f;
4321:       PetscErrorCode (**exactSol)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar[], void *);
4322:       void    **exactCtx;
4323:       PetscReal error;

4325:       SNESGetDM(snes, &dm);
4326:       SNESGetSolution(snes, &u);
4327:       DMGetDS(dm, &prob);
4328:       PetscDSGetNumFields(prob, &Nf);
4329:       PetscMalloc2(Nf, &exactSol, Nf, &exactCtx);
4330:       for (f = 0; f < Nf; ++f) PetscDSGetExactSolution(prob, f, &exactSol[f], &exactCtx[f]);
4331:       DMComputeL2Diff(dm, 0.0, exactSol, exactCtx, u, &error);
4332:       PetscFree2(exactSol, exactCtx);
4333:       if (error < 1.0e-11) PetscViewerASCIIPrintf(viewer, "L_2 Error: < 1.0e-11\n");
4334:       else PetscViewerASCIIPrintf(viewer, "L_2 Error: %g\n", (double)error);
4335:     }
4336:     if (snes->reason > 0 && format != PETSC_VIEWER_FAILED) {
4337:       if (((PetscObject)snes)->prefix) {
4338:         PetscViewerASCIIPrintf(viewer, "Nonlinear %s solve converged due to %s iterations %" PetscInt_FMT "\n", ((PetscObject)snes)->prefix, SNESConvergedReasons[snes->reason], snes->iter);
4339:       } else {
4340:         PetscViewerASCIIPrintf(viewer, "Nonlinear solve converged due to %s iterations %" PetscInt_FMT "\n", SNESConvergedReasons[snes->reason], snes->iter);
4341:       }
4342:     } else if (snes->reason <= 0) {
4343:       if (((PetscObject)snes)->prefix) {
4344:         PetscViewerASCIIPrintf(viewer, "Nonlinear %s solve did not converge due to %s iterations %" PetscInt_FMT "\n", ((PetscObject)snes)->prefix, SNESConvergedReasons[snes->reason], snes->iter);
4345:       } else {
4346:         PetscViewerASCIIPrintf(viewer, "Nonlinear solve did not converge due to %s iterations %" PetscInt_FMT "\n", SNESConvergedReasons[snes->reason], snes->iter);
4347:       }
4348:     }
4349:     PetscViewerASCIISubtractTab(viewer, ((PetscObject)snes)->tablevel);
4350:   }
4351:   return 0;
4352: }

4354: /*@C
4355:    SNESConvergedReasonViewSet - Sets an ADDITIONAL function that is to be used at the
4356:     end of the nonlinear solver to display the conver reason of the nonlinear solver.

4358:    Logically Collective on snes

4360:    Input Parameters:
4361: +  snes - the `SNES` context
4362: .  f - the snes converged reason view function
4363: .  vctx - [optional] user-defined context for private data for the
4364:           snes converged reason view routine (use NULL if no context is desired)
4365: -  reasonviewdestroy - [optional] routine that frees reasonview context
4366:           (may be NULL)

4368:    Options Database Keys:
4369: +    -snes_converged_reason        - sets a default `SNESConvergedReasonView()`
4370: -    -snes_converged_reason_view_cancel - cancels all converged reason viewers that have
4371:                             been hardwired into a code by
4372:                             calls to `SNESConvergedReasonViewSet()`, but
4373:                             does not cancel those set via
4374:                             the options database.

4376:    Note:
4377:    Several different converged reason view routines may be set by calling
4378:    `SNESConvergedReasonViewSet()` multiple times; all will be called in the
4379:    order in which they were set.

4381:    Level: intermediate

4383: .seealso: `SNES`, `SNESSolve()`, `SNESConvergedReason`, `SNESGetConvergedReason()`, `SNESConvergedReasonView()`, `SNESConvergedReasonViewCancel()`
4384: @*/
4385: PetscErrorCode SNESConvergedReasonViewSet(SNES snes, PetscErrorCode (*f)(SNES, void *), void *vctx, PetscErrorCode (*reasonviewdestroy)(void **))
4386: {
4387:   PetscInt  i;
4388:   PetscBool identical;

4391:   for (i = 0; i < snes->numberreasonviews; i++) {
4392:     PetscMonitorCompare((PetscErrorCode(*)(void))f, vctx, reasonviewdestroy, (PetscErrorCode(*)(void))snes->reasonview[i], snes->reasonviewcontext[i], snes->reasonviewdestroy[i], &identical);
4393:     if (identical) return 0;
4394:   }
4396:   snes->reasonview[snes->numberreasonviews]          = f;
4397:   snes->reasonviewdestroy[snes->numberreasonviews]   = reasonviewdestroy;
4398:   snes->reasonviewcontext[snes->numberreasonviews++] = (void *)vctx;
4399:   return 0;
4400: }

4402: /*@
4403:   SNESConvergedReasonViewFromOptions - Processes command line options to determine if/how a `SNESConvergedReason` is to be viewed.
4404:                                        All the user-provided convergedReasonView routines will be involved as well, if they exist.

4406:   Collective on snes

4408:   Input Parameters:
4409: . snes   - the `SNES` object

4411:   Level: advanced

4413: .seealso: `SNES`, `SNESConvergedReason`, `SNESConvergedReasonViewSet()`, `SNESCreate()`, `SNESSetUp()`, `SNESDestroy()`,
4414:           `SNESSetTolerances()`, `SNESConvergedDefault()`, `SNESGetConvergedReason()`, `SNESConvergedReasonView()`
4415: @*/
4416: PetscErrorCode SNESConvergedReasonViewFromOptions(SNES snes)
4417: {
4418:   PetscViewer       viewer;
4419:   PetscBool         flg;
4420:   static PetscBool  incall = PETSC_FALSE;
4421:   PetscViewerFormat format;
4422:   PetscInt          i;

4424:   if (incall) return 0;
4425:   incall = PETSC_TRUE;

4427:   /* All user-provided viewers are called first, if they exist. */
4428:   for (i = 0; i < snes->numberreasonviews; i++) (*snes->reasonview[i])(snes, snes->reasonviewcontext[i]);

4430:   /* Call PETSc default routine if users ask for it */
4431:   PetscOptionsGetViewer(PetscObjectComm((PetscObject)snes), ((PetscObject)snes)->options, ((PetscObject)snes)->prefix, "-snes_converged_reason", &viewer, &format, &flg);
4432:   if (flg) {
4433:     PetscViewerPushFormat(viewer, format);
4434:     SNESConvergedReasonView(snes, viewer);
4435:     PetscViewerPopFormat(viewer);
4436:     PetscViewerDestroy(&viewer);
4437:   }
4438:   incall = PETSC_FALSE;
4439:   return 0;
4440: }

4442: /*@
4443:    SNESSolve - Solves a nonlinear system F(x) = b.
4444:    Call `SNESSolve()` after calling `SNESCreate()` and optional routines of the form `SNESSetXXX()`.

4446:    Collective on snes

4448:    Input Parameters:
4449: +  snes - the `SNES` context
4450: .  b - the constant part of the equation F(x) = b, or NULL to use zero.
4451: -  x - the solution vector.

4453:    Note:
4454:    The user should initialize the vector,x, with the initial guess
4455:    for the nonlinear solve prior to calling `SNESSolve()`.  In particular,
4456:    to employ an initial guess of zero, the user should explicitly set
4457:    this vector to zero by calling `VecSet()`.

4459:    Level: beginner

4461: .seealso: `SNES`, `SNESCreate()`, `SNESDestroy()`, `SNESSetFunction()`, `SNESSetJacobian()`, `SNESSetGridSequence()`, `SNESGetSolution()`,
4462:           `SNESNewtonTRSetPreCheck()`, `SNESNewtonTRGetPreCheck()`, `SNESNewtonTRSetPostCheck()`, `SNESNewtonTRGetPostCheck()`,
4463:           `SNESLineSearchSetPostCheck()`, `SNESLineSearchGetPostCheck()`, `SNESLineSearchSetPreCheck()`, `SNESLineSearchGetPreCheck()`
4464: @*/
4465: PetscErrorCode SNESSolve(SNES snes, Vec b, Vec x)
4466: {
4467:   PetscBool flg;
4468:   PetscInt  grid;
4469:   Vec       xcreated = NULL;
4470:   DM        dm;


4478:   /* High level operations using the nonlinear solver */
4479:   {
4480:     PetscViewer       viewer;
4481:     PetscViewerFormat format;
4482:     PetscInt          num;
4483:     PetscBool         flg;
4484:     static PetscBool  incall = PETSC_FALSE;

4486:     if (!incall) {
4487:       /* Estimate the convergence rate of the discretization */
4488:       PetscOptionsGetViewer(PetscObjectComm((PetscObject)snes), ((PetscObject)snes)->options, ((PetscObject)snes)->prefix, "-snes_convergence_estimate", &viewer, &format, &flg);
4489:       if (flg) {
4490:         PetscConvEst conv;
4491:         DM           dm;
4492:         PetscReal   *alpha; /* Convergence rate of the solution error for each field in the L_2 norm */
4493:         PetscInt     Nf;

4495:         incall = PETSC_TRUE;
4496:         SNESGetDM(snes, &dm);
4497:         DMGetNumFields(dm, &Nf);
4498:         PetscCalloc1(Nf, &alpha);
4499:         PetscConvEstCreate(PetscObjectComm((PetscObject)snes), &conv);
4500:         PetscConvEstSetSolver(conv, (PetscObject)snes);
4501:         PetscConvEstSetFromOptions(conv);
4502:         PetscConvEstSetUp(conv);
4503:         PetscConvEstGetConvRate(conv, alpha);
4504:         PetscViewerPushFormat(viewer, format);
4505:         PetscConvEstRateView(conv, alpha, viewer);
4506:         PetscViewerPopFormat(viewer);
4507:         PetscViewerDestroy(&viewer);
4508:         PetscConvEstDestroy(&conv);
4509:         PetscFree(alpha);
4510:         incall = PETSC_FALSE;
4511:       }
4512:       /* Adaptively refine the initial grid */
4513:       num = 1;
4514:       PetscOptionsGetInt(NULL, ((PetscObject)snes)->prefix, "-snes_adapt_initial", &num, &flg);
4515:       if (flg) {
4516:         DMAdaptor adaptor;

4518:         incall = PETSC_TRUE;
4519:         DMAdaptorCreate(PetscObjectComm((PetscObject)snes), &adaptor);
4520:         DMAdaptorSetSolver(adaptor, snes);
4521:         DMAdaptorSetSequenceLength(adaptor, num);
4522:         DMAdaptorSetFromOptions(adaptor);
4523:         DMAdaptorSetUp(adaptor);
4524:         DMAdaptorAdapt(adaptor, x, DM_ADAPTATION_INITIAL, &dm, &x);
4525:         DMAdaptorDestroy(&adaptor);
4526:         incall = PETSC_FALSE;
4527:       }
4528:       /* Use grid sequencing to adapt */
4529:       num = 0;
4530:       PetscOptionsGetInt(NULL, ((PetscObject)snes)->prefix, "-snes_adapt_sequence", &num, NULL);
4531:       if (num) {
4532:         DMAdaptor adaptor;

4534:         incall = PETSC_TRUE;
4535:         DMAdaptorCreate(PetscObjectComm((PetscObject)snes), &adaptor);
4536:         DMAdaptorSetSolver(adaptor, snes);
4537:         DMAdaptorSetSequenceLength(adaptor, num);
4538:         DMAdaptorSetFromOptions(adaptor);
4539:         DMAdaptorSetUp(adaptor);
4540:         DMAdaptorAdapt(adaptor, x, DM_ADAPTATION_SEQUENTIAL, &dm, &x);
4541:         DMAdaptorDestroy(&adaptor);
4542:         incall = PETSC_FALSE;
4543:       }
4544:     }
4545:   }
4546:   if (!x) x = snes->vec_sol;
4547:   if (!x) {
4548:     SNESGetDM(snes, &dm);
4549:     DMCreateGlobalVector(dm, &xcreated);
4550:     x = xcreated;
4551:   }
4552:   SNESViewFromOptions(snes, NULL, "-snes_view_pre");

4554:   for (grid = 0; grid < snes->gridsequence; grid++) PetscViewerASCIIPushTab(PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)snes)));
4555:   for (grid = 0; grid < snes->gridsequence + 1; grid++) {
4556:     /* set solution vector */
4557:     if (!grid) PetscObjectReference((PetscObject)x);
4558:     VecDestroy(&snes->vec_sol);
4559:     snes->vec_sol = x;
4560:     SNESGetDM(snes, &dm);

4562:     /* set affine vector if provided */
4563:     if (b) PetscObjectReference((PetscObject)b);
4564:     VecDestroy(&snes->vec_rhs);
4565:     snes->vec_rhs = b;

4570:     if (!snes->vec_sol_update /* && snes->vec_sol */) { VecDuplicate(snes->vec_sol, &snes->vec_sol_update); }
4571:     DMShellSetGlobalVector(dm, snes->vec_sol);
4572:     SNESSetUp(snes);

4574:     if (!grid) {
4575:       if (snes->ops->computeinitialguess) PetscCallBack("SNES callback initial guess", (*snes->ops->computeinitialguess)(snes, snes->vec_sol, snes->initialguessP));
4576:     }

4578:     if (snes->conv_hist_reset) snes->conv_hist_len = 0;
4579:     if (snes->counters_reset) {
4580:       snes->nfuncs      = 0;
4581:       snes->linear_its  = 0;
4582:       snes->numFailures = 0;
4583:     }

4585:     PetscLogEventBegin(SNES_Solve, snes, 0, 0, 0);
4586:     PetscUseTypeMethod(snes, solve);
4587:     PetscLogEventEnd(SNES_Solve, snes, 0, 0, 0);
4589:     snes->domainerror = PETSC_FALSE; /* clear the flag if it has been set */

4591:     if (snes->lagjac_persist) snes->jac_iter += snes->iter;
4592:     if (snes->lagpre_persist) snes->pre_iter += snes->iter;

4594:     PetscOptionsGetViewer(PetscObjectComm((PetscObject)snes), ((PetscObject)snes)->options, ((PetscObject)snes)->prefix, "-snes_test_local_min", NULL, NULL, &flg);
4595:     if (flg && !PetscPreLoadingOn) SNESTestLocalMin(snes);
4596:     /* Call converged reason views. This may involve user-provided viewers as well */
4597:     SNESConvergedReasonViewFromOptions(snes);

4600:     if (snes->reason < 0) break;
4601:     if (grid < snes->gridsequence) {
4602:       DM  fine;
4603:       Vec xnew;
4604:       Mat interp;

4606:       DMRefine(snes->dm, PetscObjectComm((PetscObject)snes), &fine);
4608:       DMCreateInterpolation(snes->dm, fine, &interp, NULL);
4609:       DMCreateGlobalVector(fine, &xnew);
4610:       MatInterpolate(interp, x, xnew);
4611:       DMInterpolate(snes->dm, interp, fine);
4612:       MatDestroy(&interp);
4613:       x = xnew;

4615:       SNESReset(snes);
4616:       SNESSetDM(snes, fine);
4617:       SNESResetFromOptions(snes);
4618:       DMDestroy(&fine);
4619:       PetscViewerASCIIPopTab(PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)snes)));
4620:     }
4621:   }
4622:   SNESViewFromOptions(snes, NULL, "-snes_view");
4623:   VecViewFromOptions(snes->vec_sol, (PetscObject)snes, "-snes_view_solution");
4624:   DMMonitor(snes->dm);
4625:   SNESMonitorPauseFinal_Internal(snes);

4627:   VecDestroy(&xcreated);
4628:   PetscObjectSAWsBlock((PetscObject)snes);
4629:   return 0;
4630: }

4632: /* --------- Internal routines for SNES Package --------- */

4634: /*@C
4635:    SNESSetType - Sets the method for the nonlinear solver.

4637:    Collective on snes

4639:    Input Parameters:
4640: +  snes - the `SNES` context
4641: -  type - a known method

4643:    Options Database Key:
4644: .  -snes_type <type> - Sets the method; use -help for a list
4645:    of available methods (for instance, newtonls or newtontr)

4647:    Notes:
4648:    See "petsc/include/petscsnes.h" for available methods (for instance)
4649: +    `SNESNEWTONLS` - Newton's method with line search
4650:      (systems of nonlinear equations)
4651: -    `SNESNEWTONTRDC` - Newton's method with trust region
4652:      (systems of nonlinear equations)

4654:   Normally, it is best to use the `SNESSetFromOptions()` command and then
4655:   set the `SNES` solver type from the options database rather than by using
4656:   this routine.  Using the options database provides the user with
4657:   maximum flexibility in evaluating the many nonlinear solvers.
4658:   The `SNESSetType()` routine is provided for those situations where it
4659:   is necessary to set the nonlinear solver independently of the command
4660:   line or options database.  This might be the case, for example, when
4661:   the choice of solver changes during the execution of the program,
4662:   and the user's application is taking responsibility for choosing the
4663:   appropriate method.

4665:     Developer Note:
4666:     `SNESRegister()` adds a constructor for a new `SNESType` to `SNESList`, `SNESSetType()` locates
4667:     the constructor in that list and calls it to create the specific object.

4669:   Level: intermediate

4671: .seealso: `SNES`, `SNESSolve()`, `SNESType`, `SNESCreate()`, `SNESDestroy()`, `SNESGetType()`, `SNESSetFromOptions()`
4672: @*/
4673: PetscErrorCode SNESSetType(SNES snes, SNESType type)
4674: {
4675:   PetscBool match;
4676:   PetscErrorCode (*r)(SNES);


4681:   PetscObjectTypeCompare((PetscObject)snes, type, &match);
4682:   if (match) return 0;

4684:   PetscFunctionListFind(SNESList, type, &r);
4686:   /* Destroy the previous private SNES context */
4687:   PetscTryTypeMethod(snes, destroy);
4688:   /* Reinitialize function pointers in SNESOps structure */
4689:   snes->ops->setup          = NULL;
4690:   snes->ops->solve          = NULL;
4691:   snes->ops->view           = NULL;
4692:   snes->ops->setfromoptions = NULL;
4693:   snes->ops->destroy        = NULL;

4695:   /* It may happen the user has customized the line search before calling SNESSetType */
4696:   if (((PetscObject)snes)->type_name) SNESLineSearchDestroy(&snes->linesearch);

4698:   /* Call the SNESCreate_XXX routine for this particular Nonlinear solver */
4699:   snes->setupcalled = PETSC_FALSE;

4701:   PetscObjectChangeTypeName((PetscObject)snes, type);
4702:   (*r)(snes);
4703:   return 0;
4704: }

4706: /*@C
4707:    SNESGetType - Gets the `SNES` method type and name (as a string).

4709:    Not Collective

4711:    Input Parameter:
4712: .  snes - nonlinear solver context

4714:    Output Parameter:
4715: .  type - `SNES` method (a character string)

4717:    Level: intermediate

4719: .seealso: `SNESSetType()`, `SNESType`, `SNESSetFromOptions()`, `SNES`
4720: @*/
4721: PetscErrorCode SNESGetType(SNES snes, SNESType *type)
4722: {
4725:   *type = ((PetscObject)snes)->type_name;
4726:   return 0;
4727: }

4729: /*@
4730:   SNESSetSolution - Sets the solution vector for use by the `SNES` routines.

4732:   Logically Collective on snes

4734:   Input Parameters:
4735: + snes - the `SNES` context obtained from `SNESCreate()`
4736: - u    - the solution vector

4738:   Level: beginner

4740: .seealso: `SNES`, `SNESSolve()`, `SNESGetSolution()`, `Vec`
4741: @*/
4742: PetscErrorCode SNESSetSolution(SNES snes, Vec u)
4743: {
4744:   DM dm;

4748:   PetscObjectReference((PetscObject)u);
4749:   VecDestroy(&snes->vec_sol);

4751:   snes->vec_sol = u;

4753:   SNESGetDM(snes, &dm);
4754:   DMShellSetGlobalVector(dm, u);
4755:   return 0;
4756: }

4758: /*@
4759:    SNESGetSolution - Returns the vector where the approximate solution is
4760:    stored. This is the fine grid solution when using `SNESSetGridSequence()`.

4762:    Not Collective, but x is parallel if snes is parallel

4764:    Input Parameter:
4765: .  snes - the `SNES` context

4767:    Output Parameter:
4768: .  x - the solution

4770:    Level: intermediate

4772: .seealso: `SNESSetSolution()`, `SNESSolve()`, `SNES`, `SNESGetSolutionUpdate()`, `SNESGetFunction()`
4773: @*/
4774: PetscErrorCode SNESGetSolution(SNES snes, Vec *x)
4775: {
4778:   *x = snes->vec_sol;
4779:   return 0;
4780: }

4782: /*@
4783:    SNESGetSolutionUpdate - Returns the vector where the solution update is
4784:    stored.

4786:    Not Collective, but x is parallel if snes is parallel

4788:    Input Parameter:
4789: .  snes - the `SNES` context

4791:    Output Parameter:
4792: .  x - the solution update

4794:    Level: advanced

4796: .seealso: `SNES`, `SNESGetSolution()`, `SNESGetFunction()`
4797: @*/
4798: PetscErrorCode SNESGetSolutionUpdate(SNES snes, Vec *x)
4799: {
4802:   *x = snes->vec_sol_update;
4803:   return 0;
4804: }

4806: /*@C
4807:    SNESGetFunction - Returns the function that defines the nonlinear system set with `SNESSetFunction()`

4809:    Not Collective, but r is parallel if snes is parallel. Collective if r is requested, but has not been created yet.

4811:    Input Parameter:
4812: .  snes - the `SNES` context

4814:    Output Parameters:
4815: +  r - the vector that is used to store residuals (or NULL if you don't want it)
4816: .  f - the function (or NULL if you don't want it); see `SNESFunction` for calling sequence details
4817: -  ctx - the function context (or NULL if you don't want it)

4819:    Level: advanced

4821:     Note:
4822:    The vector r DOES NOT, in general, contain the current value of the `SNES` nonlinear function

4824: .seealso: `SNES, `SNESSolve()`, `SNESSetFunction()`, `SNESGetSolution()`, `SNESFunction`
4825: @*/
4826: PetscErrorCode SNESGetFunction(SNES snes, Vec *r, PetscErrorCode (**f)(SNES, Vec, Vec, void *), void **ctx)
4827: {
4828:   DM dm;

4831:   if (r) {
4832:     if (!snes->vec_func) {
4833:       if (snes->vec_rhs) {
4834:         VecDuplicate(snes->vec_rhs, &snes->vec_func);
4835:       } else if (snes->vec_sol) {
4836:         VecDuplicate(snes->vec_sol, &snes->vec_func);
4837:       } else if (snes->dm) {
4838:         DMCreateGlobalVector(snes->dm, &snes->vec_func);
4839:       }
4840:     }
4841:     *r = snes->vec_func;
4842:   }
4843:   SNESGetDM(snes, &dm);
4844:   DMSNESGetFunction(dm, f, ctx);
4845:   return 0;
4846: }

4848: /*@C
4849:    SNESGetNGS - Returns the `SNESNGS` function and context set with `SNESSetNGS()`

4851:    Input Parameter:
4852: .  snes - the `SNES` context

4854:    Output Parameters:
4855: +  f - the function (or NULL) see `SNESNGSFunction` for details
4856: -  ctx    - the function context (or NULL)

4858:    Level: advanced

4860: .seealso: `SNESSetNGS()`, `SNESGetFunction()`
4861: @*/

4863: PetscErrorCode SNESGetNGS(SNES snes, PetscErrorCode (**f)(SNES, Vec, Vec, void *), void **ctx)
4864: {
4865:   DM dm;

4868:   SNESGetDM(snes, &dm);
4869:   DMSNESGetNGS(dm, f, ctx);
4870:   return 0;
4871: }

4873: /*@C
4874:    SNESSetOptionsPrefix - Sets the prefix used for searching for all
4875:    `SNES` options in the database.

4877:    Logically Collective on snes

4879:    Input Parameters:
4880: +  snes - the `SNES` context
4881: -  prefix - the prefix to prepend to all option names

4883:    Note:
4884:    A hyphen (-) must NOT be given at the beginning of the prefix name.
4885:    The first character of all runtime options is AUTOMATICALLY the hyphen.

4887:    Level: advanced

4889: .seealso: `SNES`, `SNESSetFromOptions()`, `SNESAppendOptionsPrefix()`
4890: @*/
4891: PetscErrorCode SNESSetOptionsPrefix(SNES snes, const char prefix[])
4892: {
4894:   PetscObjectSetOptionsPrefix((PetscObject)snes, prefix);
4895:   if (!snes->ksp) SNESGetKSP(snes, &snes->ksp);
4896:   if (snes->linesearch) {
4897:     SNESGetLineSearch(snes, &snes->linesearch);
4898:     PetscObjectSetOptionsPrefix((PetscObject)snes->linesearch, prefix);
4899:   }
4900:   KSPSetOptionsPrefix(snes->ksp, prefix);
4901:   return 0;
4902: }

4904: /*@C
4905:    SNESAppendOptionsPrefix - Appends to the prefix used for searching for all
4906:    `SNES` options in the database.

4908:    Logically Collective on snes

4910:    Input Parameters:
4911: +  snes - the `SNES` context
4912: -  prefix - the prefix to prepend to all option names

4914:    Note:
4915:    A hyphen (-) must NOT be given at the beginning of the prefix name.
4916:    The first character of all runtime options is AUTOMATICALLY the hyphen.

4918:    Level: advanced

4920: .seealso: `SNESGetOptionsPrefix()`, `SNESSetOptionsPrefix()`
4921: @*/
4922: PetscErrorCode SNESAppendOptionsPrefix(SNES snes, const char prefix[])
4923: {
4925:   PetscObjectAppendOptionsPrefix((PetscObject)snes, prefix);
4926:   if (!snes->ksp) SNESGetKSP(snes, &snes->ksp);
4927:   if (snes->linesearch) {
4928:     SNESGetLineSearch(snes, &snes->linesearch);
4929:     PetscObjectAppendOptionsPrefix((PetscObject)snes->linesearch, prefix);
4930:   }
4931:   KSPAppendOptionsPrefix(snes->ksp, prefix);
4932:   return 0;
4933: }

4935: /*@C
4936:    SNESGetOptionsPrefix - Gets the prefix used for searching for all
4937:    `SNES` options in the database.

4939:    Not Collective

4941:    Input Parameter:
4942: .  snes - the `SNES` context

4944:    Output Parameter:
4945: .  prefix - pointer to the prefix string used

4947:    Fortran Note:
4948:     On the fortran side, the user should pass in a string 'prefix' of
4949:    sufficient length to hold the prefix.

4951:    Level: advanced

4953: .seealso: `SNES`, `SNESSetOptionsPrefix()`, `SNESAppendOptionsPrefix()`
4954: @*/
4955: PetscErrorCode SNESGetOptionsPrefix(SNES snes, const char *prefix[])
4956: {
4958:   PetscObjectGetOptionsPrefix((PetscObject)snes, prefix);
4959:   return 0;
4960: }

4962: /*@C
4963:   SNESRegister - Adds a method to the nonlinear solver package.

4965:    Not collective

4967:    Input Parameters:
4968: +  name_solver - name of a new user-defined solver
4969: -  routine_create - routine to create method context

4971:    Note:
4972:    `SNESRegister()` may be called multiple times to add several user-defined solvers.

4974:    Sample usage:
4975: .vb
4976:    SNESRegister("my_solver",MySolverCreate);
4977: .ve

4979:    Then, your solver can be chosen with the procedural interface via
4980: $     SNESSetType(snes,"my_solver")
4981:    or at runtime via the option
4982: $     -snes_type my_solver

4984:    Level: advanced

4986: .seealso: `SNESRegisterAll()`, `SNESRegisterDestroy()`
4987: @*/
4988: PetscErrorCode SNESRegister(const char sname[], PetscErrorCode (*function)(SNES))
4989: {
4990:   SNESInitializePackage();
4991:   PetscFunctionListAdd(&SNESList, sname, function);
4992:   return 0;
4993: }

4995: PetscErrorCode SNESTestLocalMin(SNES snes)
4996: {
4997:   PetscInt    N, i, j;
4998:   Vec         u, uh, fh;
4999:   PetscScalar value;
5000:   PetscReal   norm;

5002:   SNESGetSolution(snes, &u);
5003:   VecDuplicate(u, &uh);
5004:   VecDuplicate(u, &fh);

5006:   /* currently only works for sequential */
5007:   PetscPrintf(PetscObjectComm((PetscObject)snes), "Testing FormFunction() for local min\n");
5008:   VecGetSize(u, &N);
5009:   for (i = 0; i < N; i++) {
5010:     VecCopy(u, uh);
5011:     PetscPrintf(PetscObjectComm((PetscObject)snes), "i = %" PetscInt_FMT "\n", i);
5012:     for (j = -10; j < 11; j++) {
5013:       value = PetscSign(j) * PetscExpReal(PetscAbs(j) - 10.0);
5014:       VecSetValue(uh, i, value, ADD_VALUES);
5015:       SNESComputeFunction(snes, uh, fh);
5016:       VecNorm(fh, NORM_2, &norm);
5017:       PetscPrintf(PetscObjectComm((PetscObject)snes), "       j norm %" PetscInt_FMT " %18.16e\n", j, (double)norm);
5018:       value = -value;
5019:       VecSetValue(uh, i, value, ADD_VALUES);
5020:     }
5021:   }
5022:   VecDestroy(&uh);
5023:   VecDestroy(&fh);
5024:   return 0;
5025: }

5027: /*@
5028:    SNESKSPSetUseEW - Sets `SNES` to the use Eisenstat-Walker method for
5029:    computing relative tolerance for linear solvers within an inexact
5030:    Newton method.

5032:    Logically Collective on snes

5034:    Input Parameters:
5035: +  snes - `SNES` context
5036: -  flag - `PETSC_TRUE` or `PETSC_FALSE`

5038:     Options Database Keys:
5039: +  -snes_ksp_ew - use Eisenstat-Walker method for determining linear system convergence
5040: .  -snes_ksp_ew_version ver - version of  Eisenstat-Walker method
5041: .  -snes_ksp_ew_rtol0 <rtol0> - Sets rtol0
5042: .  -snes_ksp_ew_rtolmax <rtolmax> - Sets rtolmax
5043: .  -snes_ksp_ew_gamma <gamma> - Sets gamma
5044: .  -snes_ksp_ew_alpha <alpha> - Sets alpha
5045: .  -snes_ksp_ew_alpha2 <alpha2> - Sets alpha2
5046: -  -snes_ksp_ew_threshold <threshold> - Sets threshold

5048:    Note:
5049:    The default is to use a constant relative tolerance for
5050:    the inner linear solvers.  Alternatively, one can use the
5051:    Eisenstat-Walker method, where the relative convergence tolerance
5052:    is reset at each Newton iteration according progress of the nonlinear
5053:    solver.

5055:    Level: advanced

5057:    Reference:
5058: .  - * S. C. Eisenstat and H. F. Walker, "Choosing the forcing terms in an inexact Newton method", SISC 17 (1), pp.16-32, 1996.

5060: .seealso: `KSP`, `SNES`, `SNESKSPGetUseEW()`, `SNESKSPGetParametersEW()`, `SNESKSPSetParametersEW()`
5061: @*/
5062: PetscErrorCode SNESKSPSetUseEW(SNES snes, PetscBool flag)
5063: {
5066:   snes->ksp_ewconv = flag;
5067:   return 0;
5068: }

5070: /*@
5071:    SNESKSPGetUseEW - Gets if `SNES` is using Eisenstat-Walker method
5072:    for computing relative tolerance for linear solvers within an
5073:    inexact Newton method.

5075:    Not Collective

5077:    Input Parameter:
5078: .  snes - `SNES` context

5080:    Output Parameter:
5081: .  flag - `PETSC_TRUE` or `PETSC_FALSE`

5083:    Level: advanced

5085: .seealso: `SNESKSPSetUseEW()`, `SNESKSPGetParametersEW()`, `SNESKSPSetParametersEW()`
5086: @*/
5087: PetscErrorCode SNESKSPGetUseEW(SNES snes, PetscBool *flag)
5088: {
5091:   *flag = snes->ksp_ewconv;
5092:   return 0;
5093: }

5095: /*@
5096:    SNESKSPSetParametersEW - Sets parameters for Eisenstat-Walker
5097:    convergence criteria for the linear solvers within an inexact
5098:    Newton method.

5100:    Logically Collective on snes

5102:    Input Parameters:
5103: +    snes - `SNES` context
5104: .    version - version 1, 2 (default is 2), 3 or 4
5105: .    rtol_0 - initial relative tolerance (0 <= rtol_0 < 1)
5106: .    rtol_max - maximum relative tolerance (0 <= rtol_max < 1)
5107: .    gamma - multiplicative factor for version 2 rtol computation
5108:              (0 <= gamma2 <= 1)
5109: .    alpha - power for version 2 rtol computation (1 < alpha <= 2)
5110: .    alpha2 - power for safeguard
5111: -    threshold - threshold for imposing safeguard (0 < threshold < 1)

5113:    Notes:
5114:    Version 3 was contributed by Luis Chacon, June 2006.

5116:    Use `PETSC_DEFAULT` to retain the default for any of the parameters.

5118:    Level: advanced

5120: .seealso: `SNES`, `SNESKSPSetUseEW()`, `SNESKSPGetUseEW()`, `SNESKSPGetParametersEW()`
5121: @*/
5122: PetscErrorCode SNESKSPSetParametersEW(SNES snes, PetscInt version, PetscReal rtol_0, PetscReal rtol_max, PetscReal gamma, PetscReal alpha, PetscReal alpha2, PetscReal threshold)
5123: {
5124:   SNESKSPEW *kctx;

5127:   kctx = (SNESKSPEW *)snes->kspconvctx;

5137:   if (version != PETSC_DEFAULT) kctx->version = version;
5138:   if (rtol_0 != PETSC_DEFAULT) kctx->rtol_0 = rtol_0;
5139:   if (rtol_max != PETSC_DEFAULT) kctx->rtol_max = rtol_max;
5140:   if (gamma != PETSC_DEFAULT) kctx->gamma = gamma;
5141:   if (alpha != PETSC_DEFAULT) kctx->alpha = alpha;
5142:   if (alpha2 != PETSC_DEFAULT) kctx->alpha2 = alpha2;
5143:   if (threshold != PETSC_DEFAULT) kctx->threshold = threshold;

5151:   return 0;
5152: }

5154: /*@
5155:    SNESKSPGetParametersEW - Gets parameters for Eisenstat-Walker
5156:    convergence criteria for the linear solvers within an inexact
5157:    Newton method.

5159:    Not Collective

5161:    Input Parameter:
5162: .    snes - `SNES` context

5164:    Output Parameters:
5165: +    version - version 1, 2 (default is 2), 3 or 4
5166: .    rtol_0 - initial relative tolerance (0 <= rtol_0 < 1)
5167: .    rtol_max - maximum relative tolerance (0 <= rtol_max < 1)
5168: .    gamma - multiplicative factor for version 2 rtol computation (0 <= gamma2 <= 1)
5169: .    alpha - power for version 2 rtol computation (1 < alpha <= 2)
5170: .    alpha2 - power for safeguard
5171: -    threshold - threshold for imposing safeguard (0 < threshold < 1)

5173:    Level: advanced

5175: .seealso: `SNES`, `SNESKSPSetUseEW()`, `SNESKSPGetUseEW()`, `SNESKSPSetParametersEW()`
5176: @*/
5177: PetscErrorCode SNESKSPGetParametersEW(SNES snes, PetscInt *version, PetscReal *rtol_0, PetscReal *rtol_max, PetscReal *gamma, PetscReal *alpha, PetscReal *alpha2, PetscReal *threshold)
5178: {
5179:   SNESKSPEW *kctx;

5182:   kctx = (SNESKSPEW *)snes->kspconvctx;
5184:   if (version) *version = kctx->version;
5185:   if (rtol_0) *rtol_0 = kctx->rtol_0;
5186:   if (rtol_max) *rtol_max = kctx->rtol_max;
5187:   if (gamma) *gamma = kctx->gamma;
5188:   if (alpha) *alpha = kctx->alpha;
5189:   if (alpha2) *alpha2 = kctx->alpha2;
5190:   if (threshold) *threshold = kctx->threshold;
5191:   return 0;
5192: }

5194: PetscErrorCode KSPPreSolve_SNESEW(KSP ksp, Vec b, Vec x, SNES snes)
5195: {
5196:   SNESKSPEW *kctx = (SNESKSPEW *)snes->kspconvctx;
5197:   PetscReal  rtol = PETSC_DEFAULT, stol;

5199:   if (!snes->ksp_ewconv) return 0;
5200:   if (!snes->iter) {
5201:     rtol = kctx->rtol_0; /* first time in, so use the original user rtol */
5202:     VecNorm(snes->vec_func, NORM_2, &kctx->norm_first);
5203:   } else {
5204:     if (kctx->version == 1) {
5205:       rtol = PetscAbsReal(snes->norm - kctx->lresid_last) / kctx->norm_last;
5206:       stol = PetscPowReal(kctx->rtol_last, kctx->alpha2);
5207:       if (stol > kctx->threshold) rtol = PetscMax(rtol, stol);
5208:     } else if (kctx->version == 2) {
5209:       rtol = kctx->gamma * PetscPowReal(snes->norm / kctx->norm_last, kctx->alpha);
5210:       stol = kctx->gamma * PetscPowReal(kctx->rtol_last, kctx->alpha);
5211:       if (stol > kctx->threshold) rtol = PetscMax(rtol, stol);
5212:     } else if (kctx->version == 3) { /* contributed by Luis Chacon, June 2006. */
5213:       rtol = kctx->gamma * PetscPowReal(snes->norm / kctx->norm_last, kctx->alpha);
5214:       /* safeguard: avoid sharp decrease of rtol */
5215:       stol = kctx->gamma * PetscPowReal(kctx->rtol_last, kctx->alpha);
5216:       stol = PetscMax(rtol, stol);
5217:       rtol = PetscMin(kctx->rtol_0, stol);
5218:       /* safeguard: avoid oversolving */
5219:       stol = kctx->gamma * (kctx->norm_first * snes->rtol) / snes->norm;
5220:       stol = PetscMax(rtol, stol);
5221:       rtol = PetscMin(kctx->rtol_0, stol);
5222:     } else if (kctx->version == 4) { /* H.-B. An et al. Journal of Computational and Applied Mathematics 200 (2007) 47-60 */
5223:       PetscReal ared = PetscAbsReal(kctx->norm_last - snes->norm);
5224:       PetscReal pred = PetscAbsReal(kctx->norm_last - kctx->lresid_last);
5225:       PetscReal rk   = ared / pred;
5226:       if (rk < kctx->v4_p1) rtol = 1. - 2. * kctx->v4_p1;
5227:       else if (rk < kctx->v4_p2) rtol = kctx->rtol_last;
5228:       else if (rk < kctx->v4_p3) rtol = kctx->v4_m1 * kctx->rtol_last;
5229:       else rtol = kctx->v4_m2 * kctx->rtol_last;

5231:       if (kctx->rtol_last_2 > kctx->v4_m3 && kctx->rtol_last > kctx->v4_m3 && kctx->rk_last_2 < kctx->v4_p1 && kctx->rk_last < kctx->v4_p1) {
5232:         rtol = kctx->v4_m4 * kctx->rtol_last;
5233:         //printf("iter %" PetscInt_FMT ", Eisenstat-Walker (version %" PetscInt_FMT ") KSP rtol=%g (rk %g ps %g %g %g) (AD)\n",snes->iter,kctx->version,(double)rtol,rk,kctx->v4_p1,kctx->v4_p2,kctx->v4_p3);
5234:       } else {
5235:         //printf("iter %" PetscInt_FMT ", Eisenstat-Walker (version %" PetscInt_FMT ") KSP rtol=%g (rk %g ps %g %g %g)\n",snes->iter,kctx->version,(double)rtol,rk,kctx->v4_p1,kctx->v4_p2,kctx->v4_p3);
5236:       }
5237:       kctx->rtol_last_2 = kctx->rtol_last;
5238:       kctx->rk_last_2   = kctx->rk_last;
5239:       kctx->rk_last     = rk;
5240:     } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Only versions 1-4 are supported: %" PetscInt_FMT, kctx->version);
5241:   }
5242:   /* safeguard: avoid rtol greater than rtol_max */
5243:   rtol = PetscMin(rtol, kctx->rtol_max);
5244:   KSPSetTolerances(ksp, rtol, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT);
5245:   PetscInfo(snes, "iter %" PetscInt_FMT ", Eisenstat-Walker (version %" PetscInt_FMT ") KSP rtol=%g\n", snes->iter, kctx->version, (double)rtol);
5246:   return 0;
5247: }

5249: PetscErrorCode KSPPostSolve_SNESEW(KSP ksp, Vec b, Vec x, SNES snes)
5250: {
5251:   SNESKSPEW *kctx = (SNESKSPEW *)snes->kspconvctx;
5252:   PCSide     pcside;
5253:   Vec        lres;

5255:   if (!snes->ksp_ewconv) return 0;
5256:   KSPGetTolerances(ksp, &kctx->rtol_last, NULL, NULL, NULL);
5257:   kctx->norm_last = snes->norm;
5258:   if (kctx->version == 1 || kctx->version == 4) {
5259:     PC        pc;
5260:     PetscBool getRes;

5262:     KSPGetPC(ksp, &pc);
5263:     PetscObjectTypeCompare((PetscObject)pc, PCNONE, &getRes);
5264:     if (!getRes) {
5265:       KSPNormType normtype;

5267:       KSPGetNormType(ksp, &normtype);
5268:       getRes = (PetscBool)(normtype == KSP_NORM_UNPRECONDITIONED);
5269:     }
5270:     KSPGetPCSide(ksp, &pcside);
5271:     if (pcside == PC_RIGHT || getRes) { /* KSP residual is true linear residual */
5272:       KSPGetResidualNorm(ksp, &kctx->lresid_last);
5273:     } else {
5274:       /* KSP residual is preconditioned residual */
5275:       /* compute true linear residual norm */
5276:       Mat J;
5277:       KSPGetOperators(ksp, &J, NULL);
5278:       VecDuplicate(b, &lres);
5279:       MatMult(J, x, lres);
5280:       VecAYPX(lres, -1.0, b);
5281:       VecNorm(lres, NORM_2, &kctx->lresid_last);
5282:       VecDestroy(&lres);
5283:     }
5284:   }
5285:   return 0;
5286: }

5288: /*@
5289:    SNESGetKSP - Returns the `KSP` context for a `SNES` solver.

5291:    Not Collective, but if snes is parallel, then ksp is parallel

5293:    Input Parameter:
5294: .  snes - the `SNES` context

5296:    Output Parameter:
5297: .  ksp - the `KSP` context

5299:    Notes:
5300:    The user can then directly manipulate the `KSP` context to set various
5301:    options, etc.  Likewise, the user can then extract and manipulate the
5302:    `PC` contexts as well.

5304:    Some `SNESType`s do not use a `KSP` but a `KSP` is still returned by this function

5306:    Level: beginner

5308: .seealso: `SNES`, `KSP`, `PC`, `KSPGetPC()`, `SNESCreate()`, `KSPCreate()`, `SNESSetKSP()`
5309: @*/
5310: PetscErrorCode SNESGetKSP(SNES snes, KSP *ksp)
5311: {

5315:   if (!snes->ksp) {
5316:     KSPCreate(PetscObjectComm((PetscObject)snes), &snes->ksp);
5317:     PetscObjectIncrementTabLevel((PetscObject)snes->ksp, (PetscObject)snes, 1);

5319:     KSPSetPreSolve(snes->ksp, (PetscErrorCode(*)(KSP, Vec, Vec, void *))KSPPreSolve_SNESEW, snes);
5320:     KSPSetPostSolve(snes->ksp, (PetscErrorCode(*)(KSP, Vec, Vec, void *))KSPPostSolve_SNESEW, snes);

5322:     KSPMonitorSetFromOptions(snes->ksp, "-snes_monitor_ksp", "snes_preconditioned_residual", snes);
5323:     PetscObjectSetOptions((PetscObject)snes->ksp, ((PetscObject)snes)->options);
5324:   }
5325:   *ksp = snes->ksp;
5326:   return 0;
5327: }

5329: #include <petsc/private/dmimpl.h>
5330: /*@
5331:    SNESSetDM - Sets the `DM` that may be used by some nonlinear solvers or their underlying preconditioners

5333:    Logically Collective on snes

5335:    Input Parameters:
5336: +  snes - the nonlinear solver context
5337: -  dm - the dm, cannot be NULL

5339:    Note:
5340:    A `DM` can only be used for solving one problem at a time because information about the problem is stored on the `DM`,
5341:    even when not using interfaces like `DMSNESSetFunction()`.  Use `DMClone()` to get a distinct `DM` when solving different
5342:    problems using the same function space.

5344:    Level: intermediate

5346: .seealso: `DM`, `SNESGetDM()`, `KSPSetDM()`, `KSPGetDM()`
5347: @*/
5348: PetscErrorCode SNESSetDM(SNES snes, DM dm)
5349: {
5350:   KSP    ksp;
5351:   DMSNES sdm;

5355:   PetscObjectReference((PetscObject)dm);
5356:   if (snes->dm) { /* Move the DMSNES context over to the new DM unless the new DM already has one */
5357:     if (snes->dm->dmsnes && !dm->dmsnes) {
5358:       DMCopyDMSNES(snes->dm, dm);
5359:       DMGetDMSNES(snes->dm, &sdm);
5360:       if (sdm->originaldm == snes->dm) sdm->originaldm = dm; /* Grant write privileges to the replacement DM */
5361:     }
5362:     DMCoarsenHookRemove(snes->dm, DMCoarsenHook_SNESVecSol, DMRestrictHook_SNESVecSol, snes);
5363:     DMDestroy(&snes->dm);
5364:   }
5365:   snes->dm     = dm;
5366:   snes->dmAuto = PETSC_FALSE;

5368:   SNESGetKSP(snes, &ksp);
5369:   KSPSetDM(ksp, dm);
5370:   KSPSetDMActive(ksp, PETSC_FALSE);
5371:   if (snes->npc) {
5372:     SNESSetDM(snes->npc, snes->dm);
5373:     SNESSetNPCSide(snes, snes->npcside);
5374:   }
5375:   return 0;
5376: }

5378: /*@
5379:    SNESGetDM - Gets the `DM` that may be used by some preconditioners

5381:    Not Collective but dm obtained is parallel on snes

5383:    Input Parameter:
5384: . snes - the preconditioner context

5386:    Output Parameter:
5387: .  dm - the dm

5389:    Level: intermediate

5391: .seealso: `DM`, `SNESSetDM()`, `KSPSetDM()`, `KSPGetDM()`
5392: @*/
5393: PetscErrorCode SNESGetDM(SNES snes, DM *dm)
5394: {
5396:   if (!snes->dm) {
5397:     DMShellCreate(PetscObjectComm((PetscObject)snes), &snes->dm);
5398:     snes->dmAuto = PETSC_TRUE;
5399:   }
5400:   *dm = snes->dm;
5401:   return 0;
5402: }

5404: /*@
5405:   SNESSetNPC - Sets the nonlinear preconditioner to be used.

5407:   Collective on snes

5409:   Input Parameters:
5410: + snes - iterative context obtained from `SNESCreate()`
5411: - npc   - the preconditioner object

5413:   Notes:
5414:   Use `SNESGetNPC()` to retrieve the preconditioner context (for example,
5415:   to configure it using the API).

5417:   Only some `SNESType` can use a nonlinear preconditioner

5419:   Level: developer

5421: .seealso: `SNESNGS`, `SNESFAS`, `SNESGetNPC()`, `SNESHasNPC()`
5422: @*/
5423: PetscErrorCode SNESSetNPC(SNES snes, SNES npc)
5424: {
5428:   PetscObjectReference((PetscObject)npc);
5429:   SNESDestroy(&snes->npc);
5430:   snes->npc = npc;
5431:   return 0;
5432: }

5434: /*@
5435:   SNESGetNPC - Gets a nonlinear preconditioning solver SNES` to be used to precondition the original nonlinear solver.

5437:   Not Collective; but any changes to the obtained the npc object must be applied collectively

5439:   Input Parameter:
5440: . snes - iterative context obtained from `SNESCreate()`

5442:   Output Parameter:
5443: . npc - preconditioner context

5445:   Options Database Key:
5446: . -npc_snes_type <type> - set the type of the `SNES` to use as the nonlinear preconditioner

5448:   Notes:
5449:     If a `SNES` was previously set with `SNESSetNPC()` then that value is returned, otherwise a new `SNES` object is created.

5451:     The (preconditioner) `SNES` returned automatically inherits the same nonlinear function and Jacobian supplied to the original
5452:     `SNES`

5454:   Level: developer

5456: .seealso: `SNESSetNPC()`, `SNESHasNPC()`, `SNES`, `SNESCreate()`
5457: @*/
5458: PetscErrorCode SNESGetNPC(SNES snes, SNES *pc)
5459: {
5460:   const char *optionsprefix;

5464:   if (!snes->npc) {
5465:     SNESCreate(PetscObjectComm((PetscObject)snes), &snes->npc);
5466:     PetscObjectIncrementTabLevel((PetscObject)snes->npc, (PetscObject)snes, 1);
5467:     SNESGetOptionsPrefix(snes, &optionsprefix);
5468:     SNESSetOptionsPrefix(snes->npc, optionsprefix);
5469:     SNESAppendOptionsPrefix(snes->npc, "npc_");
5470:     SNESSetCountersReset(snes->npc, PETSC_FALSE);
5471:   }
5472:   *pc = snes->npc;
5473:   return 0;
5474: }

5476: /*@
5477:   SNESHasNPC - Returns whether a nonlinear preconditioner exists

5479:   Not Collective

5481:   Input Parameter:
5482: . snes - iterative context obtained from `SNESCreate()`

5484:   Output Parameter:
5485: . has_npc - whether the `SNES` has an NPC or not

5487:   Level: developer

5489: .seealso: `SNESSetNPC()`, `SNESGetNPC()`
5490: @*/
5491: PetscErrorCode SNESHasNPC(SNES snes, PetscBool *has_npc)
5492: {
5494:   *has_npc = (PetscBool)(snes->npc ? PETSC_TRUE : PETSC_FALSE);
5495:   return 0;
5496: }

5498: /*@
5499:     SNESSetNPCSide - Sets the preconditioning side.

5501:     Logically Collective on snes

5503:     Input Parameter:
5504: .   snes - iterative context obtained from `SNESCreate()`

5506:     Output Parameter:
5507: .   side - the preconditioning side, where side is one of
5508: .vb
5509:       PC_LEFT - left preconditioning
5510:       PC_RIGHT - right preconditioning (default for most nonlinear solvers)
5511: .ve

5513:     Options Database Key:
5514: .   -snes_npc_side <right,left> - nonlinear preconditioner side

5516:     Note:
5517:     `SNESNRICHARDSON` and `SNESNCG` only support left preconditioning.

5519:     Level: intermediate

5521: .seealso: `SNESType`, `SNESGetNPCSide()`, `KSPSetPCSide()`
5522: @*/
5523: PetscErrorCode SNESSetNPCSide(SNES snes, PCSide side)
5524: {
5527:   if (side == PC_SIDE_DEFAULT) side = PC_RIGHT;
5529:   snes->npcside = side;
5530:   return 0;
5531: }

5533: /*@
5534:     SNESGetNPCSide - Gets the preconditioning side.

5536:     Not Collective

5538:     Input Parameter:
5539: .   snes - iterative context obtained from `SNESCreate()`

5541:     Output Parameter:
5542: .   side - the preconditioning side, where side is one of
5543: .vb
5544:       `PC_LEFT` - left preconditioning
5545:       `PC_RIGHT` - right preconditioning (default for most nonlinear solvers)
5546: .ve

5548:     Level: intermediate

5550: .seealso: `SNES`, `SNESSetNPCSide()`, `KSPGetPCSide()`
5551: @*/
5552: PetscErrorCode SNESGetNPCSide(SNES snes, PCSide *side)
5553: {
5556:   *side = snes->npcside;
5557:   return 0;
5558: }

5560: /*@
5561:   SNESSetLineSearch - Sets the linesearch on the `SNES` instance.

5563:   Collective on snes

5565:   Input Parameters:
5566: + snes - iterative context obtained from `SNESCreate()`
5567: - linesearch   - the linesearch object

5569:   Note:
5570:   Use `SNESGetLineSearch()` to retrieve the preconditioner context (for example,
5571:   to configure it using the API).

5573:   Level: developer

5575: .seealso: `SNESGetLineSearch()`
5576: @*/
5577: PetscErrorCode SNESSetLineSearch(SNES snes, SNESLineSearch linesearch)
5578: {
5582:   PetscObjectReference((PetscObject)linesearch);
5583:   SNESLineSearchDestroy(&snes->linesearch);

5585:   snes->linesearch = linesearch;

5587:   return 0;
5588: }

5590: /*@
5591:   SNESGetLineSearch - Returns a pointer to the line search context set with `SNESSetLineSearch()`
5592:   or creates a default line search instance associated with the `SNES` and returns it.

5594:   Not Collective

5596:   Input Parameter:
5597: . snes - iterative context obtained from `SNESCreate()`

5599:   Output Parameter:
5600: . linesearch - linesearch context

5602:   Level: beginner

5604: .seealso: `SNESLineSearch`, `SNESSetLineSearch()`, `SNESLineSearchCreate()`
5605: @*/
5606: PetscErrorCode SNESGetLineSearch(SNES snes, SNESLineSearch *linesearch)
5607: {
5608:   const char *optionsprefix;

5612:   if (!snes->linesearch) {
5613:     SNESGetOptionsPrefix(snes, &optionsprefix);
5614:     SNESLineSearchCreate(PetscObjectComm((PetscObject)snes), &snes->linesearch);
5615:     SNESLineSearchSetSNES(snes->linesearch, snes);
5616:     SNESLineSearchAppendOptionsPrefix(snes->linesearch, optionsprefix);
5617:     PetscObjectIncrementTabLevel((PetscObject)snes->linesearch, (PetscObject)snes, 1);
5618:   }
5619:   *linesearch = snes->linesearch;
5620:   return 0;
5621: }