Actual source code: snes.c
petsc-3.14.6 2021-03-30
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 : this takes an optional truth value (0/1/no/yes/true/false)
27: Level: intermediate
29: Notes:
30: Normally PETSc continues if a linear solver fails to converge, you can call SNESGetConvergedReason() after a SNESSolve()
31: to determine if it has converged.
33: .seealso: SNESGetErrorIfNotConverged(), KSPGetErrorIfNotConverged(), KSPSetErrorIFNotConverged()
34: @*/
35: PetscErrorCode SNESSetErrorIfNotConverged(SNES snes,PetscBool flg)
36: {
40: snes->errorifnotconverged = flg;
41: return(0);
42: }
44: /*@
45: SNESGetErrorIfNotConverged - Will SNESSolve() generate an error if the solver does not converge?
47: Not Collective
49: Input Parameter:
50: . snes - iterative context obtained from SNESCreate()
52: Output Parameter:
53: . flag - PETSC_TRUE if it will generate an error, else PETSC_FALSE
55: Level: intermediate
57: .seealso: SNESSetErrorIfNotConverged(), KSPGetErrorIfNotConverged(), KSPSetErrorIFNotConverged()
58: @*/
59: PetscErrorCode SNESGetErrorIfNotConverged(SNES snes,PetscBool *flag)
60: {
64: *flag = snes->errorifnotconverged;
65: return(0);
66: }
68: /*@
69: SNESSetAlwaysComputesFinalResidual - does the SNES always compute the residual at the final solution?
71: Logically Collective on SNES
73: Input Parameters:
74: + snes - the shell SNES
75: - flg - is the residual computed?
77: Level: advanced
79: .seealso: SNESGetAlwaysComputesFinalResidual()
80: @*/
81: PetscErrorCode SNESSetAlwaysComputesFinalResidual(SNES snes, PetscBool flg)
82: {
85: snes->alwayscomputesfinalresidual = flg;
86: return(0);
87: }
89: /*@
90: SNESGetAlwaysComputesFinalResidual - does the SNES always compute the residual at the final solution?
92: Logically Collective on SNES
94: Input Parameter:
95: . snes - the shell SNES
97: Output Parameter:
98: . flg - is the residual computed?
100: Level: advanced
102: .seealso: SNESSetAlwaysComputesFinalResidual()
103: @*/
104: PetscErrorCode SNESGetAlwaysComputesFinalResidual(SNES snes, PetscBool *flg)
105: {
108: *flg = snes->alwayscomputesfinalresidual;
109: return(0);
110: }
112: /*@
113: SNESSetFunctionDomainError - tells SNES that the input vector to your SNESFunction is not
114: in the functions domain. For example, negative pressure.
116: Logically Collective on SNES
118: Input Parameters:
119: . snes - the SNES context
121: Level: advanced
123: .seealso: SNESCreate(), SNESSetFunction(), SNESFunction
124: @*/
125: PetscErrorCode SNESSetFunctionDomainError(SNES snes)
126: {
129: if (snes->errorifnotconverged) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"User code indicates input vector is not in the function domain");
130: snes->domainerror = PETSC_TRUE;
131: return(0);
132: }
134: /*@
135: SNESSetJacobianDomainError - tells SNES that computeJacobian does not make sense any more. For example there is a negative element transformation.
137: Logically Collective on SNES
139: Input Parameters:
140: . snes - the SNES context
142: Level: advanced
144: .seealso: SNESCreate(), SNESSetFunction(), SNESFunction(), SNESSetFunctionDomainError()
145: @*/
146: PetscErrorCode SNESSetJacobianDomainError(SNES snes)
147: {
150: if (snes->errorifnotconverged) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"User code indicates computeJacobian does not make sense");
151: snes->jacobiandomainerror = PETSC_TRUE;
152: return(0);
153: }
155: /*@
156: SNESSetCheckJacobianDomainError - if or not to check jacobian domain error after each Jacobian evaluation. By default, we check Jacobian domain error
157: in the debug mode, and do not check it in the optimized mode.
159: Logically Collective on SNES
161: Input Parameters:
162: + snes - the SNES context
163: - flg - indicates if or not to check jacobian domain error after each Jacobian evaluation
165: Level: advanced
167: .seealso: SNESCreate(), SNESSetFunction(), SNESFunction(), SNESSetFunctionDomainError(), SNESGetCheckJacobianDomainError()
168: @*/
169: PetscErrorCode SNESSetCheckJacobianDomainError(SNES snes, PetscBool flg)
170: {
173: snes->checkjacdomainerror = flg;
174: return(0);
175: }
177: /*@
178: SNESGetCheckJacobianDomainError - Get an indicator whether or not we are checking Jacobian domain errors after each Jacobian evaluation.
180: Logically Collective on SNES
182: Input Parameters:
183: . snes - the SNES context
185: Output Parameters:
186: . flg - PETSC_FALSE indicates that we don't check jacobian domain errors after each Jacobian evaluation
188: Level: advanced
190: .seealso: SNESCreate(), SNESSetFunction(), SNESFunction(), SNESSetFunctionDomainError(), SNESSetCheckJacobianDomainError()
191: @*/
192: PetscErrorCode SNESGetCheckJacobianDomainError(SNES snes, PetscBool *flg)
193: {
197: *flg = snes->checkjacdomainerror;
198: return(0);
199: }
201: /*@
202: SNESGetFunctionDomainError - Gets the status of the domain error after a call to SNESComputeFunction;
204: Logically Collective on SNES
206: Input Parameters:
207: . snes - the SNES context
209: Output Parameters:
210: . domainerror - Set to PETSC_TRUE if there's a domain error; PETSC_FALSE otherwise.
212: Level: advanced
214: .seealso: SNESSetFunctionDomainError(), SNESComputeFunction()
215: @*/
216: PetscErrorCode SNESGetFunctionDomainError(SNES snes, PetscBool *domainerror)
217: {
221: *domainerror = snes->domainerror;
222: return(0);
223: }
225: /*@
226: SNESGetJacobianDomainError - Gets the status of the Jacobian domain error after a call to SNESComputeJacobian;
228: Logically Collective on SNES
230: Input Parameters:
231: . snes - the SNES context
233: Output Parameters:
234: . domainerror - Set to PETSC_TRUE if there's a jacobian domain error; PETSC_FALSE otherwise.
236: Level: advanced
238: .seealso: SNESSetFunctionDomainError(), SNESComputeFunction(),SNESGetFunctionDomainError()
239: @*/
240: PetscErrorCode SNESGetJacobianDomainError(SNES snes, PetscBool *domainerror)
241: {
245: *domainerror = snes->jacobiandomainerror;
246: return(0);
247: }
249: /*@C
250: SNESLoad - Loads a SNES that has been stored in binary with SNESView().
252: Collective on PetscViewer
254: Input Parameters:
255: + newdm - the newly loaded SNES, this needs to have been created with SNESCreate() or
256: some related function before a call to SNESLoad().
257: - viewer - binary file viewer, obtained from PetscViewerBinaryOpen()
259: Level: intermediate
261: Notes:
262: The type is determined by the data in the file, any type set into the SNES before this call is ignored.
264: Notes for advanced users:
265: Most users should not need to know the details of the binary storage
266: format, since SNESLoad() and TSView() completely hide these details.
267: But for anyone who's interested, the standard binary matrix storage
268: format is
269: .vb
270: has not yet been determined
271: .ve
273: .seealso: PetscViewerBinaryOpen(), SNESView(), MatLoad(), VecLoad()
274: @*/
275: PetscErrorCode SNESLoad(SNES snes, PetscViewer viewer)
276: {
278: PetscBool isbinary;
279: PetscInt classid;
280: char type[256];
281: KSP ksp;
282: DM dm;
283: DMSNES dmsnes;
288: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
289: if (!isbinary) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Invalid viewer; open viewer with PetscViewerBinaryOpen()");
291: PetscViewerBinaryRead(viewer,&classid,1,NULL,PETSC_INT);
292: if (classid != SNES_FILE_CLASSID) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONG,"Not SNES next in file");
293: PetscViewerBinaryRead(viewer,type,256,NULL,PETSC_CHAR);
294: SNESSetType(snes, type);
295: if (snes->ops->load) {
296: (*snes->ops->load)(snes,viewer);
297: }
298: SNESGetDM(snes,&dm);
299: DMGetDMSNES(dm,&dmsnes);
300: DMSNESLoad(dmsnes,viewer);
301: SNESGetKSP(snes,&ksp);
302: KSPLoad(ksp,viewer);
303: return(0);
304: }
306: #include <petscdraw.h>
307: #if defined(PETSC_HAVE_SAWS)
308: #include <petscviewersaws.h>
309: #endif
311: /*@C
312: SNESViewFromOptions - View from Options
314: Collective on SNES
316: Input Parameters:
317: + A - the application ordering context
318: . obj - Optional object
319: - name - command line option
321: Level: intermediate
322: .seealso: SNES, SNESView, PetscObjectViewFromOptions(), SNESCreate()
323: @*/
324: PetscErrorCode SNESViewFromOptions(SNES A,PetscObject obj,const char name[])
325: {
330: PetscObjectViewFromOptions((PetscObject)A,obj,name);
331: return(0);
332: }
334: PETSC_EXTERN PetscErrorCode SNESComputeJacobian_DMDA(SNES,Vec,Mat,Mat,void*);
336: /*@C
337: SNESView - Prints the SNES data structure.
339: Collective on SNES
341: Input Parameters:
342: + SNES - the SNES context
343: - viewer - visualization context
345: Options Database Key:
346: . -snes_view - Calls SNESView() at end of SNESSolve()
348: Notes:
349: The available visualization contexts include
350: + PETSC_VIEWER_STDOUT_SELF - standard output (default)
351: - PETSC_VIEWER_STDOUT_WORLD - synchronized standard
352: output where only the first processor opens
353: the file. All other processors send their
354: data to the first processor to print.
356: The user can open an alternative visualization context with
357: PetscViewerASCIIOpen() - output to a specified file.
359: In the debugger you can do "call SNESView(snes,0)" to display the SNES solver. (The same holds for any PETSc object viewer).
361: Level: beginner
363: .seealso: PetscViewerASCIIOpen()
364: @*/
365: PetscErrorCode SNESView(SNES snes,PetscViewer viewer)
366: {
367: SNESKSPEW *kctx;
369: KSP ksp;
370: SNESLineSearch linesearch;
371: PetscBool iascii,isstring,isbinary,isdraw;
372: DMSNES dmsnes;
373: #if defined(PETSC_HAVE_SAWS)
374: PetscBool issaws;
375: #endif
379: if (!viewer) {
380: PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)snes),&viewer);
381: }
385: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
386: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSTRING,&isstring);
387: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
388: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
389: #if defined(PETSC_HAVE_SAWS)
390: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSAWS,&issaws);
391: #endif
392: if (iascii) {
393: SNESNormSchedule normschedule;
394: DM dm;
395: PetscErrorCode (*cJ)(SNES,Vec,Mat,Mat,void*);
396: void *ctx;
397: const char *pre = "";
399: PetscObjectPrintClassNamePrefixType((PetscObject)snes,viewer);
400: if (!snes->setupcalled) {
401: PetscViewerASCIIPrintf(viewer," SNES has not been set up so information may be incomplete\n");
402: }
403: if (snes->ops->view) {
404: PetscViewerASCIIPushTab(viewer);
405: (*snes->ops->view)(snes,viewer);
406: PetscViewerASCIIPopTab(viewer);
407: }
408: PetscViewerASCIIPrintf(viewer," maximum iterations=%D, maximum function evaluations=%D\n",snes->max_its,snes->max_funcs);
409: PetscViewerASCIIPrintf(viewer," tolerances: relative=%g, absolute=%g, solution=%g\n",(double)snes->rtol,(double)snes->abstol,(double)snes->stol);
410: if (snes->usesksp) {
411: PetscViewerASCIIPrintf(viewer," total number of linear solver iterations=%D\n",snes->linear_its);
412: }
413: PetscViewerASCIIPrintf(viewer," total number of function evaluations=%D\n",snes->nfuncs);
414: SNESGetNormSchedule(snes, &normschedule);
415: if (normschedule > 0) {PetscViewerASCIIPrintf(viewer," norm schedule %s\n",SNESNormSchedules[normschedule]);}
416: if (snes->gridsequence) {
417: PetscViewerASCIIPrintf(viewer," total number of grid sequence refinements=%D\n",snes->gridsequence);
418: }
419: if (snes->ksp_ewconv) {
420: kctx = (SNESKSPEW*)snes->kspconvctx;
421: if (kctx) {
422: PetscViewerASCIIPrintf(viewer," Eisenstat-Walker computation of KSP relative tolerance (version %D)\n",kctx->version);
423: PetscViewerASCIIPrintf(viewer," rtol_0=%g, rtol_max=%g, threshold=%g\n",(double)kctx->rtol_0,(double)kctx->rtol_max,(double)kctx->threshold);
424: PetscViewerASCIIPrintf(viewer," gamma=%g, alpha=%g, alpha2=%g\n",(double)kctx->gamma,(double)kctx->alpha,(double)kctx->alpha2);
425: }
426: }
427: if (snes->lagpreconditioner == -1) {
428: PetscViewerASCIIPrintf(viewer," Preconditioned is never rebuilt\n");
429: } else if (snes->lagpreconditioner > 1) {
430: PetscViewerASCIIPrintf(viewer," Preconditioned is rebuilt every %D new Jacobians\n",snes->lagpreconditioner);
431: }
432: if (snes->lagjacobian == -1) {
433: PetscViewerASCIIPrintf(viewer," Jacobian is never rebuilt\n");
434: } else if (snes->lagjacobian > 1) {
435: PetscViewerASCIIPrintf(viewer," Jacobian is rebuilt every %D SNES iterations\n",snes->lagjacobian);
436: }
437: SNESGetDM(snes,&dm);
438: DMSNESGetJacobian(dm,&cJ,&ctx);
439: if (snes->mf_operator) {
440: PetscViewerASCIIPrintf(viewer," Jacobian is applied matrix-free with differencing\n");
441: pre = "Preconditioning ";
442: }
443: if (cJ == SNESComputeJacobianDefault) {
444: PetscViewerASCIIPrintf(viewer," %sJacobian is built using finite differences one column at a time\n",pre);
445: } else if (cJ == SNESComputeJacobianDefaultColor) {
446: PetscViewerASCIIPrintf(viewer," %sJacobian is built using finite differences with coloring\n",pre);
447: /* it slightly breaks data encapsulation for access the DMDA information directly */
448: } else if (cJ == SNESComputeJacobian_DMDA) {
449: MatFDColoring fdcoloring;
450: PetscObjectQuery((PetscObject)dm,"DMDASNES_FDCOLORING",(PetscObject*)&fdcoloring);
451: if (fdcoloring) {
452: PetscViewerASCIIPrintf(viewer," %sJacobian is built using colored finite differences on a DMDA\n",pre);
453: } else {
454: PetscViewerASCIIPrintf(viewer," %sJacobian is built using a DMDA local Jacobian\n",pre);
455: }
456: } else if (snes->mf) {
457: PetscViewerASCIIPrintf(viewer," Jacobian is applied matrix-free with differencing, no explict Jacobian\n");
458: }
459: } else if (isstring) {
460: const char *type;
461: SNESGetType(snes,&type);
462: PetscViewerStringSPrintf(viewer," SNESType: %-7.7s",type);
463: if (snes->ops->view) {(*snes->ops->view)(snes,viewer);}
464: } else if (isbinary) {
465: PetscInt classid = SNES_FILE_CLASSID;
466: MPI_Comm comm;
467: PetscMPIInt rank;
468: char type[256];
470: PetscObjectGetComm((PetscObject)snes,&comm);
471: MPI_Comm_rank(comm,&rank);
472: if (!rank) {
473: PetscViewerBinaryWrite(viewer,&classid,1,PETSC_INT);
474: PetscStrncpy(type,((PetscObject)snes)->type_name,sizeof(type));
475: PetscViewerBinaryWrite(viewer,type,sizeof(type),PETSC_CHAR);
476: }
477: if (snes->ops->view) {
478: (*snes->ops->view)(snes,viewer);
479: }
480: } else if (isdraw) {
481: PetscDraw draw;
482: char str[36];
483: PetscReal x,y,bottom,h;
485: PetscViewerDrawGetDraw(viewer,0,&draw);
486: PetscDrawGetCurrentPoint(draw,&x,&y);
487: PetscStrncpy(str,"SNES: ",sizeof(str));
488: PetscStrlcat(str,((PetscObject)snes)->type_name,sizeof(str));
489: PetscDrawStringBoxed(draw,x,y,PETSC_DRAW_BLUE,PETSC_DRAW_BLACK,str,NULL,&h);
490: bottom = y - h;
491: PetscDrawPushCurrentPoint(draw,x,bottom);
492: if (snes->ops->view) {
493: (*snes->ops->view)(snes,viewer);
494: }
495: #if defined(PETSC_HAVE_SAWS)
496: } else if (issaws) {
497: PetscMPIInt rank;
498: const char *name;
500: PetscObjectGetName((PetscObject)snes,&name);
501: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
502: if (!((PetscObject)snes)->amsmem && !rank) {
503: char dir[1024];
505: PetscObjectViewSAWs((PetscObject)snes,viewer);
506: PetscSNPrintf(dir,1024,"/PETSc/Objects/%s/its",name);
507: PetscStackCallSAWs(SAWs_Register,(dir,&snes->iter,1,SAWs_READ,SAWs_INT));
508: if (!snes->conv_hist) {
509: SNESSetConvergenceHistory(snes,NULL,NULL,PETSC_DECIDE,PETSC_TRUE);
510: }
511: PetscSNPrintf(dir,1024,"/PETSc/Objects/%s/conv_hist",name);
512: PetscStackCallSAWs(SAWs_Register,(dir,snes->conv_hist,10,SAWs_READ,SAWs_DOUBLE));
513: }
514: #endif
515: }
516: if (snes->linesearch) {
517: SNESGetLineSearch(snes, &linesearch);
518: PetscViewerASCIIPushTab(viewer);
519: SNESLineSearchView(linesearch, viewer);
520: PetscViewerASCIIPopTab(viewer);
521: }
522: if (snes->npc && snes->usesnpc) {
523: PetscViewerASCIIPushTab(viewer);
524: SNESView(snes->npc, viewer);
525: PetscViewerASCIIPopTab(viewer);
526: }
527: PetscViewerASCIIPushTab(viewer);
528: DMGetDMSNES(snes->dm,&dmsnes);
529: DMSNESView(dmsnes, viewer);
530: PetscViewerASCIIPopTab(viewer);
531: if (snes->usesksp) {
532: SNESGetKSP(snes,&ksp);
533: PetscViewerASCIIPushTab(viewer);
534: KSPView(ksp,viewer);
535: PetscViewerASCIIPopTab(viewer);
536: }
537: if (isdraw) {
538: PetscDraw draw;
539: PetscViewerDrawGetDraw(viewer,0,&draw);
540: PetscDrawPopCurrentPoint(draw);
541: }
542: return(0);
543: }
545: /*
546: We retain a list of functions that also take SNES command
547: line options. These are called at the end SNESSetFromOptions()
548: */
549: #define MAXSETFROMOPTIONS 5
550: static PetscInt numberofsetfromoptions;
551: static PetscErrorCode (*othersetfromoptions[MAXSETFROMOPTIONS])(SNES);
553: /*@C
554: SNESAddOptionsChecker - Adds an additional function to check for SNES options.
556: Not Collective
558: Input Parameter:
559: . snescheck - function that checks for options
561: Level: developer
563: .seealso: SNESSetFromOptions()
564: @*/
565: PetscErrorCode SNESAddOptionsChecker(PetscErrorCode (*snescheck)(SNES))
566: {
568: if (numberofsetfromoptions >= MAXSETFROMOPTIONS) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE, "Too many options checkers, only %D allowed", MAXSETFROMOPTIONS);
569: othersetfromoptions[numberofsetfromoptions++] = snescheck;
570: return(0);
571: }
573: PETSC_INTERN PetscErrorCode SNESDefaultMatrixFreeCreate2(SNES,Vec,Mat*);
575: static PetscErrorCode SNESSetUpMatrixFree_Private(SNES snes, PetscBool hasOperator, PetscInt version)
576: {
577: Mat J;
579: MatNullSpace nullsp;
584: if (!snes->vec_func && (snes->jacobian || snes->jacobian_pre)) {
585: Mat A = snes->jacobian, B = snes->jacobian_pre;
586: MatCreateVecs(A ? A : B, NULL,&snes->vec_func);
587: }
589: if (version == 1) {
590: MatCreateSNESMF(snes,&J);
591: MatMFFDSetOptionsPrefix(J,((PetscObject)snes)->prefix);
592: MatSetFromOptions(J);
593: } else if (version == 2) {
594: if (!snes->vec_func) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"SNESSetFunction() must be called first");
595: #if !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_REAL_SINGLE) && !defined(PETSC_USE_REAL___FLOAT128) && !defined(PETSC_USE_REAL___FP16)
596: SNESDefaultMatrixFreeCreate2(snes,snes->vec_func,&J);
597: #else
598: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP, "matrix-free operator routines (version 2)");
599: #endif
600: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE, "matrix-free operator routines, only version 1 and 2");
602: /* attach any user provided null space that was on Amat to the newly created matrix free matrix */
603: if (snes->jacobian) {
604: MatGetNullSpace(snes->jacobian,&nullsp);
605: if (nullsp) {
606: MatSetNullSpace(J,nullsp);
607: }
608: }
610: PetscInfo1(snes,"Setting default matrix-free operator routines (version %D)\n", version);
611: if (hasOperator) {
613: /* This version replaces the user provided Jacobian matrix with a
614: matrix-free version but still employs the user-provided preconditioner matrix. */
615: SNESSetJacobian(snes,J,NULL,NULL,NULL);
616: } else {
617: /* This version replaces both the user-provided Jacobian and the user-
618: provided preconditioner Jacobian with the default matrix free version. */
619: if ((snes->npcside== PC_LEFT) && snes->npc) {
620: if (!snes->jacobian){SNESSetJacobian(snes,J,NULL,NULL,NULL);}
621: } else {
622: KSP ksp;
623: PC pc;
624: PetscBool match;
626: SNESSetJacobian(snes,J,J,MatMFFDComputeJacobian,NULL);
627: /* Force no preconditioner */
628: SNESGetKSP(snes,&ksp);
629: KSPGetPC(ksp,&pc);
630: PetscObjectTypeCompare((PetscObject)pc,PCSHELL,&match);
631: if (!match) {
632: PetscInfo(snes,"Setting default matrix-free preconditioner routines\nThat is no preconditioner is being used\n");
633: PCSetType(pc,PCNONE);
634: }
635: }
636: }
637: MatDestroy(&J);
638: return(0);
639: }
641: static PetscErrorCode DMRestrictHook_SNESVecSol(DM dmfine,Mat Restrict,Vec Rscale,Mat Inject,DM dmcoarse,void *ctx)
642: {
643: SNES snes = (SNES)ctx;
645: Vec Xfine,Xfine_named = NULL,Xcoarse;
648: if (PetscLogPrintInfo) {
649: PetscInt finelevel,coarselevel,fineclevel,coarseclevel;
650: DMGetRefineLevel(dmfine,&finelevel);
651: DMGetCoarsenLevel(dmfine,&fineclevel);
652: DMGetRefineLevel(dmcoarse,&coarselevel);
653: DMGetCoarsenLevel(dmcoarse,&coarseclevel);
654: PetscInfo4(dmfine,"Restricting SNES solution vector from level %D-%D to level %D-%D\n",finelevel,fineclevel,coarselevel,coarseclevel);
655: }
656: if (dmfine == snes->dm) Xfine = snes->vec_sol;
657: else {
658: DMGetNamedGlobalVector(dmfine,"SNESVecSol",&Xfine_named);
659: Xfine = Xfine_named;
660: }
661: DMGetNamedGlobalVector(dmcoarse,"SNESVecSol",&Xcoarse);
662: if (Inject) {
663: MatRestrict(Inject,Xfine,Xcoarse);
664: } else {
665: MatRestrict(Restrict,Xfine,Xcoarse);
666: VecPointwiseMult(Xcoarse,Xcoarse,Rscale);
667: }
668: DMRestoreNamedGlobalVector(dmcoarse,"SNESVecSol",&Xcoarse);
669: if (Xfine_named) {DMRestoreNamedGlobalVector(dmfine,"SNESVecSol",&Xfine_named);}
670: return(0);
671: }
673: static PetscErrorCode DMCoarsenHook_SNESVecSol(DM dm,DM dmc,void *ctx)
674: {
678: DMCoarsenHookAdd(dmc,DMCoarsenHook_SNESVecSol,DMRestrictHook_SNESVecSol,ctx);
679: return(0);
680: }
682: /* This may be called to rediscretize the operator on levels of linear multigrid. The DM shuffle is so the user can
683: * safely call SNESGetDM() in their residual evaluation routine. */
684: static PetscErrorCode KSPComputeOperators_SNES(KSP ksp,Mat A,Mat B,void *ctx)
685: {
686: SNES snes = (SNES)ctx;
688: Vec X,Xnamed = NULL;
689: DM dmsave;
690: void *ctxsave;
691: PetscErrorCode (*jac)(SNES,Vec,Mat,Mat,void*) = NULL;
694: dmsave = snes->dm;
695: KSPGetDM(ksp,&snes->dm);
696: if (dmsave == snes->dm) X = snes->vec_sol; /* We are on the finest level */
697: else { /* We are on a coarser level, this vec was initialized using a DM restrict hook */
698: DMGetNamedGlobalVector(snes->dm,"SNESVecSol",&Xnamed);
699: X = Xnamed;
700: SNESGetJacobian(snes,NULL,NULL,&jac,&ctxsave);
701: /* If the DM's don't match up, the MatFDColoring context needed for the jacobian won't match up either -- fixit. */
702: if (jac == SNESComputeJacobianDefaultColor) {
703: SNESSetJacobian(snes,NULL,NULL,SNESComputeJacobianDefaultColor,NULL);
704: }
705: }
706: /* Make sure KSP DM has the Jacobian computation routine */
707: {
708: DMSNES sdm;
710: DMGetDMSNES(snes->dm, &sdm);
711: if (!sdm->ops->computejacobian) {
712: DMCopyDMSNES(dmsave, snes->dm);
713: }
714: }
715: /* Compute the operators */
716: SNESComputeJacobian(snes,X,A,B);
717: /* Put the previous context back */
718: if (snes->dm != dmsave && jac == SNESComputeJacobianDefaultColor) {
719: SNESSetJacobian(snes,NULL,NULL,jac,ctxsave);
720: }
722: if (Xnamed) {DMRestoreNamedGlobalVector(snes->dm,"SNESVecSol",&Xnamed);}
723: snes->dm = dmsave;
724: return(0);
725: }
727: /*@
728: SNESSetUpMatrices - ensures that matrices are available for SNES, to be called by SNESSetUp_XXX()
730: Collective
732: Input Arguments:
733: . snes - snes to configure
735: Level: developer
737: .seealso: SNESSetUp()
738: @*/
739: PetscErrorCode SNESSetUpMatrices(SNES snes)
740: {
742: DM dm;
743: DMSNES sdm;
746: SNESGetDM(snes,&dm);
747: DMGetDMSNES(dm,&sdm);
748: if (!snes->jacobian && snes->mf) {
749: Mat J;
750: void *functx;
751: MatCreateSNESMF(snes,&J);
752: MatMFFDSetOptionsPrefix(J,((PetscObject)snes)->prefix);
753: MatSetFromOptions(J);
754: SNESGetFunction(snes,NULL,NULL,&functx);
755: SNESSetJacobian(snes,J,J,NULL,NULL);
756: MatDestroy(&J);
757: } else if (snes->mf_operator && !snes->jacobian_pre && !snes->jacobian) {
758: Mat J,B;
759: MatCreateSNESMF(snes,&J);
760: MatMFFDSetOptionsPrefix(J,((PetscObject)snes)->prefix);
761: MatSetFromOptions(J);
762: DMCreateMatrix(snes->dm,&B);
763: /* sdm->computejacobian was already set to reach here */
764: SNESSetJacobian(snes,J,B,NULL,NULL);
765: MatDestroy(&J);
766: MatDestroy(&B);
767: } else if (!snes->jacobian_pre) {
768: PetscErrorCode (*nspconstr)(DM, PetscInt, PetscInt, MatNullSpace *);
769: PetscDS prob;
770: Mat J, B;
771: MatNullSpace nullspace = NULL;
772: PetscBool hasPrec = PETSC_FALSE;
773: PetscInt Nf;
775: J = snes->jacobian;
776: DMGetDS(dm, &prob);
777: if (prob) {PetscDSHasJacobianPreconditioner(prob, &hasPrec);}
778: if (J) {PetscObjectReference((PetscObject) J);}
779: else if (hasPrec) {DMCreateMatrix(snes->dm, &J);}
780: DMCreateMatrix(snes->dm, &B);
781: PetscDSGetNumFields(prob, &Nf);
782: DMGetNullSpaceConstructor(snes->dm, Nf, &nspconstr);
783: if (nspconstr) (*nspconstr)(snes->dm, Nf, Nf, &nullspace);
784: MatSetNullSpace(B, nullspace);
785: MatNullSpaceDestroy(&nullspace);
786: SNESSetJacobian(snes, J ? J : B, B, NULL, NULL);
787: MatDestroy(&J);
788: MatDestroy(&B);
789: }
790: {
791: KSP ksp;
792: SNESGetKSP(snes,&ksp);
793: KSPSetComputeOperators(ksp,KSPComputeOperators_SNES,snes);
794: DMCoarsenHookAdd(snes->dm,DMCoarsenHook_SNESVecSol,DMRestrictHook_SNESVecSol,snes);
795: }
796: return(0);
797: }
799: /*@C
800: SNESMonitorSetFromOptions - Sets a monitor function and viewer appropriate for the type indicated by the user
802: Collective on SNES
804: Input Parameters:
805: + snes - SNES object you wish to monitor
806: . name - the monitor type one is seeking
807: . help - message indicating what monitoring is done
808: . manual - manual page for the monitor
809: . monitor - the monitor function
810: - 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
812: Level: developer
814: .seealso: PetscOptionsGetViewer(), PetscOptionsGetReal(), PetscOptionsHasName(), PetscOptionsGetString(),
815: PetscOptionsGetIntArray(), PetscOptionsGetRealArray(), PetscOptionsBool()
816: PetscOptionsInt(), PetscOptionsString(), PetscOptionsReal(), PetscOptionsBool(),
817: PetscOptionsName(), PetscOptionsBegin(), PetscOptionsEnd(), PetscOptionsHead(),
818: PetscOptionsStringArray(),PetscOptionsRealArray(), PetscOptionsScalar(),
819: PetscOptionsBoolGroupBegin(), PetscOptionsBoolGroup(), PetscOptionsBoolGroupEnd(),
820: PetscOptionsFList(), PetscOptionsEList()
821: @*/
822: PetscErrorCode SNESMonitorSetFromOptions(SNES snes,const char name[],const char help[], const char manual[],PetscErrorCode (*monitor)(SNES,PetscInt,PetscReal,PetscViewerAndFormat*),PetscErrorCode (*monitorsetup)(SNES,PetscViewerAndFormat*))
823: {
824: PetscErrorCode ierr;
825: PetscViewer viewer;
826: PetscViewerFormat format;
827: PetscBool flg;
830: PetscOptionsGetViewer(PetscObjectComm((PetscObject)snes),((PetscObject) snes)->options,((PetscObject)snes)->prefix,name,&viewer,&format,&flg);
831: if (flg) {
832: PetscViewerAndFormat *vf;
833: PetscViewerAndFormatCreate(viewer,format,&vf);
834: PetscObjectDereference((PetscObject)viewer);
835: if (monitorsetup) {
836: (*monitorsetup)(snes,vf);
837: }
838: SNESMonitorSet(snes,(PetscErrorCode (*)(SNES,PetscInt,PetscReal,void*))monitor,vf,(PetscErrorCode (*)(void**))PetscViewerAndFormatDestroy);
839: }
840: return(0);
841: }
843: /*@
844: SNESSetFromOptions - Sets various SNES and KSP parameters from user options.
846: Collective on SNES
848: Input Parameter:
849: . snes - the SNES context
851: Options Database Keys:
852: + -snes_type <type> - newtonls, newtontr, ngmres, ncg, nrichardson, qn, vi, fas, SNESType for complete list
853: . -snes_stol - convergence tolerance in terms of the norm
854: of the change in the solution between steps
855: . -snes_atol <abstol> - absolute tolerance of residual norm
856: . -snes_rtol <rtol> - relative decrease in tolerance norm from initial
857: . -snes_divergence_tolerance <divtol> - if the residual goes above divtol*rnorm0, exit with divergence
858: . -snes_force_iteration <force> - force SNESSolve() to take at least one iteration
859: . -snes_max_it <max_it> - maximum number of iterations
860: . -snes_max_funcs <max_funcs> - maximum number of function evaluations
861: . -snes_max_fail <max_fail> - maximum number of line search failures allowed before stopping, default is none
862: . -snes_max_linear_solve_fail - number of linear solver failures before SNESSolve() stops
863: . -snes_lag_preconditioner <lag> - how often preconditioner is rebuilt (use -1 to never rebuild)
864: . -snes_lag_preconditioner_persists <true,false> - retains the -snes_lag_preconditioner information across multiple SNESSolve()
865: . -snes_lag_jacobian <lag> - how often Jacobian is rebuilt (use -1 to never rebuild)
866: . -snes_lag_jacobian_persists <true,false> - retains the -snes_lag_jacobian information across multiple SNESSolve()
867: . -snes_trtol <trtol> - trust region tolerance
868: . -snes_no_convergence_test - skip convergence test in nonlinear
869: solver; hence iterations will continue until max_it
870: or some other criterion is reached. Saves expense
871: of convergence test
872: . -snes_monitor [ascii][:filename][:viewer format] - prints residual norm at each iteration. if no filename given prints to stdout
873: . -snes_monitor_solution [ascii binary draw][:filename][:viewer format] - plots solution at each iteration
874: . -snes_monitor_residual [ascii binary draw][:filename][:viewer format] - plots residual (not its norm) at each iteration
875: . -snes_monitor_solution_update [ascii binary draw][:filename][:viewer format] - plots update to solution at each iteration
876: . -snes_monitor_lg_residualnorm - plots residual norm at each iteration
877: . -snes_monitor_lg_range - plots residual norm at each iteration
878: . -snes_fd - use finite differences to compute Jacobian; very slow, only for testing
879: . -snes_fd_color - use finite differences with coloring to compute Jacobian
880: . -snes_mf_ksp_monitor - if using matrix-free multiply then print h at each KSP iteration
881: . -snes_converged_reason - print the reason for convergence/divergence after each solve
882: . -npc_snes_type <type> - the SNES type to use as a nonlinear preconditioner
883: . -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.
884: - -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.
886: Options Database for Eisenstat-Walker method:
887: + -snes_ksp_ew - use Eisenstat-Walker method for determining linear system convergence
888: . -snes_ksp_ew_version ver - version of Eisenstat-Walker method
889: . -snes_ksp_ew_rtol0 <rtol0> - Sets rtol0
890: . -snes_ksp_ew_rtolmax <rtolmax> - Sets rtolmax
891: . -snes_ksp_ew_gamma <gamma> - Sets gamma
892: . -snes_ksp_ew_alpha <alpha> - Sets alpha
893: . -snes_ksp_ew_alpha2 <alpha2> - Sets alpha2
894: - -snes_ksp_ew_threshold <threshold> - Sets threshold
896: Notes:
897: To see all options, run your program with the -help option or consult the users manual
899: Notes:
900: SNES supports three approaches for computing (approximate) Jacobians: user provided via SNESSetJacobian(), matrix free, and computing explictly with
901: finite differences and coloring using MatFDColoring. It is also possible to use automatic differentiation and the MatFDColoring object.
903: Level: beginner
905: .seealso: SNESSetOptionsPrefix(), SNESResetFromOptions(), SNES, SNESCreate()
906: @*/
907: PetscErrorCode SNESSetFromOptions(SNES snes)
908: {
909: PetscBool flg,pcset,persist,set;
910: PetscInt i,indx,lag,grids;
911: const char *deft = SNESNEWTONLS;
912: const char *convtests[] = {"default","skip"};
913: SNESKSPEW *kctx = NULL;
914: char type[256], monfilename[PETSC_MAX_PATH_LEN];
916: PCSide pcside;
917: const char *optionsprefix;
921: SNESRegisterAll();
922: PetscObjectOptionsBegin((PetscObject)snes);
923: if (((PetscObject)snes)->type_name) deft = ((PetscObject)snes)->type_name;
924: PetscOptionsFList("-snes_type","Nonlinear solver method","SNESSetType",SNESList,deft,type,256,&flg);
925: if (flg) {
926: SNESSetType(snes,type);
927: } else if (!((PetscObject)snes)->type_name) {
928: SNESSetType(snes,deft);
929: }
930: PetscOptionsReal("-snes_stol","Stop if step length less than","SNESSetTolerances",snes->stol,&snes->stol,NULL);
931: PetscOptionsReal("-snes_atol","Stop if function norm less than","SNESSetTolerances",snes->abstol,&snes->abstol,NULL);
933: PetscOptionsReal("-snes_rtol","Stop if decrease in function norm less than","SNESSetTolerances",snes->rtol,&snes->rtol,NULL);
934: PetscOptionsReal("-snes_divergence_tolerance","Stop if residual norm increases by this factor","SNESSetDivergenceTolerance",snes->divtol,&snes->divtol,NULL);
935: PetscOptionsInt("-snes_max_it","Maximum iterations","SNESSetTolerances",snes->max_its,&snes->max_its,NULL);
936: PetscOptionsInt("-snes_max_funcs","Maximum function evaluations","SNESSetTolerances",snes->max_funcs,&snes->max_funcs,NULL);
937: PetscOptionsInt("-snes_max_fail","Maximum nonlinear step failures","SNESSetMaxNonlinearStepFailures",snes->maxFailures,&snes->maxFailures,NULL);
938: PetscOptionsInt("-snes_max_linear_solve_fail","Maximum failures in linear solves allowed","SNESSetMaxLinearSolveFailures",snes->maxLinearSolveFailures,&snes->maxLinearSolveFailures,NULL);
939: PetscOptionsBool("-snes_error_if_not_converged","Generate error if solver does not converge","SNESSetErrorIfNotConverged",snes->errorifnotconverged,&snes->errorifnotconverged,NULL);
940: PetscOptionsBool("-snes_force_iteration","Force SNESSolve() to take at least one iteration","SNESSetForceIteration",snes->forceiteration,&snes->forceiteration,NULL);
941: PetscOptionsBool("-snes_check_jacobian_domain_error","Check Jacobian domain error after Jacobian evaluation","SNESCheckJacobianDomainError",snes->checkjacdomainerror,&snes->checkjacdomainerror,NULL);
943: PetscOptionsInt("-snes_lag_preconditioner","How often to rebuild preconditioner","SNESSetLagPreconditioner",snes->lagpreconditioner,&lag,&flg);
944: if (flg) {
945: if (lag == -1) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_USER,"Cannot set the lag to -1 from the command line since the preconditioner must be built as least once, perhaps you mean -2");
946: SNESSetLagPreconditioner(snes,lag);
947: }
948: PetscOptionsBool("-snes_lag_preconditioner_persists","Preconditioner lagging through multiple SNES solves","SNESSetLagPreconditionerPersists",snes->lagjac_persist,&persist,&flg);
949: if (flg) {
950: SNESSetLagPreconditionerPersists(snes,persist);
951: }
952: PetscOptionsInt("-snes_lag_jacobian","How often to rebuild Jacobian","SNESSetLagJacobian",snes->lagjacobian,&lag,&flg);
953: if (flg) {
954: if (lag == -1) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_USER,"Cannot set the lag to -1 from the command line since the Jacobian must be built as least once, perhaps you mean -2");
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) {
959: SNESSetLagJacobianPersists(snes,persist);
960: }
962: PetscOptionsInt("-snes_grid_sequence","Use grid sequencing to generate initial guess","SNESSetGridSequence",snes->gridsequence,&grids,&flg);
963: if (flg) {
964: SNESSetGridSequence(snes,grids);
965: }
967: PetscOptionsEList("-snes_convergence_test","Convergence test","SNESSetConvergenceTest",convtests,2,"default",&indx,&flg);
968: if (flg) {
969: switch (indx) {
970: case 0: SNESSetConvergenceTest(snes,SNESConvergedDefault,NULL,NULL); break;
971: case 1: SNESSetConvergenceTest(snes,SNESConvergedSkip,NULL,NULL); break;
972: }
973: }
975: PetscOptionsEList("-snes_norm_schedule","SNES Norm schedule","SNESSetNormSchedule",SNESNormSchedules,5,"function",&indx,&flg);
976: if (flg) { SNESSetNormSchedule(snes,(SNESNormSchedule)indx); }
978: PetscOptionsEList("-snes_function_type","SNES Norm schedule","SNESSetFunctionType",SNESFunctionTypes,2,"unpreconditioned",&indx,&flg);
979: if (flg) { SNESSetFunctionType(snes,(SNESFunctionType)indx); }
981: kctx = (SNESKSPEW*)snes->kspconvctx;
983: PetscOptionsBool("-snes_ksp_ew","Use Eisentat-Walker linear system convergence test","SNESKSPSetUseEW",snes->ksp_ewconv,&snes->ksp_ewconv,NULL);
985: PetscOptionsInt("-snes_ksp_ew_version","Version 1, 2 or 3","SNESKSPSetParametersEW",kctx->version,&kctx->version,NULL);
986: PetscOptionsReal("-snes_ksp_ew_rtol0","0 <= rtol0 < 1","SNESKSPSetParametersEW",kctx->rtol_0,&kctx->rtol_0,NULL);
987: PetscOptionsReal("-snes_ksp_ew_rtolmax","0 <= rtolmax < 1","SNESKSPSetParametersEW",kctx->rtol_max,&kctx->rtol_max,NULL);
988: PetscOptionsReal("-snes_ksp_ew_gamma","0 <= gamma <= 1","SNESKSPSetParametersEW",kctx->gamma,&kctx->gamma,NULL);
989: PetscOptionsReal("-snes_ksp_ew_alpha","1 < alpha <= 2","SNESKSPSetParametersEW",kctx->alpha,&kctx->alpha,NULL);
990: PetscOptionsReal("-snes_ksp_ew_alpha2","alpha2","SNESKSPSetParametersEW",kctx->alpha2,&kctx->alpha2,NULL);
991: PetscOptionsReal("-snes_ksp_ew_threshold","0 < threshold < 1","SNESKSPSetParametersEW",kctx->threshold,&kctx->threshold,NULL);
993: flg = PETSC_FALSE;
994: PetscOptionsBool("-snes_monitor_cancel","Remove all monitors","SNESMonitorCancel",flg,&flg,&set);
995: if (set && flg) {SNESMonitorCancel(snes);}
997: SNESMonitorSetFromOptions(snes,"-snes_monitor","Monitor norm of function","SNESMonitorDefault",SNESMonitorDefault,NULL);
998: SNESMonitorSetFromOptions(snes,"-snes_monitor_short","Monitor norm of function with fewer digits","SNESMonitorDefaultShort",SNESMonitorDefaultShort,NULL);
999: SNESMonitorSetFromOptions(snes,"-snes_monitor_range","Monitor range of elements of function","SNESMonitorRange",SNESMonitorRange,NULL);
1001: SNESMonitorSetFromOptions(snes,"-snes_monitor_ratio","Monitor ratios of the norm of function for consecutive steps","SNESMonitorRatio",SNESMonitorRatio,SNESMonitorRatioSetUp);
1002: SNESMonitorSetFromOptions(snes,"-snes_monitor_field","Monitor norm of function (split into fields)","SNESMonitorDefaultField",SNESMonitorDefaultField,NULL);
1003: SNESMonitorSetFromOptions(snes,"-snes_monitor_solution","View solution at each iteration","SNESMonitorSolution",SNESMonitorSolution,NULL);
1004: SNESMonitorSetFromOptions(snes,"-snes_monitor_solution_update","View correction at each iteration","SNESMonitorSolutionUpdate",SNESMonitorSolutionUpdate,NULL);
1005: SNESMonitorSetFromOptions(snes,"-snes_monitor_residual","View residual at each iteration","SNESMonitorResidual",SNESMonitorResidual,NULL);
1006: SNESMonitorSetFromOptions(snes,"-snes_monitor_jacupdate_spectrum","Print the change in the spectrum of the Jacobian","SNESMonitorJacUpdateSpectrum",SNESMonitorJacUpdateSpectrum,NULL);
1007: SNESMonitorSetFromOptions(snes,"-snes_monitor_fields","Monitor norm of function per field","SNESMonitorSet",SNESMonitorFields,NULL);
1009: PetscOptionsString("-snes_monitor_python","Use Python function","SNESMonitorSet",NULL,monfilename,sizeof(monfilename),&flg);
1010: if (flg) {PetscPythonMonitorSet((PetscObject)snes,monfilename);}
1012: flg = PETSC_FALSE;
1013: PetscOptionsBool("-snes_monitor_lg_residualnorm","Plot function norm at each iteration","SNESMonitorLGResidualNorm",flg,&flg,NULL);
1014: if (flg) {
1015: PetscDrawLG ctx;
1017: SNESMonitorLGCreate(PetscObjectComm((PetscObject)snes),NULL,NULL,PETSC_DECIDE,PETSC_DECIDE,400,300,&ctx);
1018: SNESMonitorSet(snes,SNESMonitorLGResidualNorm,ctx,(PetscErrorCode (*)(void**))PetscDrawLGDestroy);
1019: }
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_fd","Use finite differences (slow) to compute Jacobian","SNESComputeJacobianDefault",flg,&flg,NULL);
1031: if (flg) {
1032: void *functx;
1033: DM dm;
1034: DMSNES sdm;
1035: SNESGetDM(snes,&dm);
1036: DMGetDMSNES(dm,&sdm);
1037: sdm->jacobianctx = NULL;
1038: SNESGetFunction(snes,NULL,NULL,&functx);
1039: SNESSetJacobian(snes,snes->jacobian,snes->jacobian_pre,SNESComputeJacobianDefault,functx);
1040: PetscInfo(snes,"Setting default finite difference Jacobian matrix\n");
1041: }
1043: flg = PETSC_FALSE;
1044: PetscOptionsBool("-snes_fd_function","Use finite differences (slow) to compute function from user objective","SNESObjectiveComputeFunctionDefaultFD",flg,&flg,NULL);
1045: if (flg) {
1046: SNESSetFunction(snes,NULL,SNESObjectiveComputeFunctionDefaultFD,NULL);
1047: }
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: DMSNES sdm;
1054: SNESGetDM(snes,&dm);
1055: DMGetDMSNES(dm,&sdm);
1056: sdm->jacobianctx = NULL;
1057: SNESSetJacobian(snes,snes->jacobian,snes->jacobian_pre,SNESComputeJacobianDefaultColor,NULL);
1058: PetscInfo(snes,"Setting default finite difference coloring Jacobian matrix\n");
1059: }
1061: flg = PETSC_FALSE;
1062: PetscOptionsBool("-snes_mf_operator","Use a Matrix-Free Jacobian with user-provided preconditioner matrix","SNESSetUseMatrixFree",PETSC_FALSE,&snes->mf_operator,&flg);
1063: if (flg && snes->mf_operator) {
1064: snes->mf_operator = PETSC_TRUE;
1065: snes->mf = PETSC_TRUE;
1066: }
1067: flg = PETSC_FALSE;
1068: PetscOptionsBool("-snes_mf","Use a Matrix-Free Jacobian with no preconditioner matrix","SNESSetUseMatrixFree",PETSC_FALSE,&snes->mf,&flg);
1069: if (!flg && snes->mf_operator) snes->mf = PETSC_TRUE;
1070: PetscOptionsInt("-snes_mf_version","Matrix-Free routines version 1 or 2","None",snes->mf_version,&snes->mf_version,NULL);
1072: flg = PETSC_FALSE;
1073: SNESGetNPCSide(snes,&pcside);
1074: PetscOptionsEnum("-snes_npc_side","SNES nonlinear preconditioner side","SNESSetNPCSide",PCSides,(PetscEnum)pcside,(PetscEnum*)&pcside,&flg);
1075: if (flg) {SNESSetNPCSide(snes,pcside);}
1077: #if defined(PETSC_HAVE_SAWS)
1078: /*
1079: Publish convergence information using SAWs
1080: */
1081: flg = PETSC_FALSE;
1082: PetscOptionsBool("-snes_monitor_saws","Publish SNES progress using SAWs","SNESMonitorSet",flg,&flg,NULL);
1083: if (flg) {
1084: void *ctx;
1085: SNESMonitorSAWsCreate(snes,&ctx);
1086: SNESMonitorSet(snes,SNESMonitorSAWs,ctx,SNESMonitorSAWsDestroy);
1087: }
1088: #endif
1089: #if defined(PETSC_HAVE_SAWS)
1090: {
1091: PetscBool set;
1092: flg = PETSC_FALSE;
1093: PetscOptionsBool("-snes_saws_block","Block for SAWs at end of SNESSolve","PetscObjectSAWsBlock",((PetscObject)snes)->amspublishblock,&flg,&set);
1094: if (set) {
1095: PetscObjectSAWsSetBlock((PetscObject)snes,flg);
1096: }
1097: }
1098: #endif
1100: for (i = 0; i < numberofsetfromoptions; i++) {
1101: (*othersetfromoptions[i])(snes);
1102: }
1104: if (snes->ops->setfromoptions) {
1105: (*snes->ops->setfromoptions)(PetscOptionsObject,snes);
1106: }
1108: /* process any options handlers added with PetscObjectAddOptionsHandler() */
1109: PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject)snes);
1110: PetscOptionsEnd();
1112: if (snes->linesearch) {
1113: SNESGetLineSearch(snes, &snes->linesearch);
1114: SNESLineSearchSetFromOptions(snes->linesearch);
1115: }
1117: if (snes->usesksp) {
1118: if (!snes->ksp) {SNESGetKSP(snes,&snes->ksp);}
1119: KSPSetOperators(snes->ksp,snes->jacobian,snes->jacobian_pre);
1120: KSPSetFromOptions(snes->ksp);
1121: }
1123: /* if user has set the SNES NPC type via options database, create it. */
1124: SNESGetOptionsPrefix(snes, &optionsprefix);
1125: PetscOptionsHasName(((PetscObject)snes)->options,optionsprefix, "-npc_snes_type", &pcset);
1126: if (pcset && (!snes->npc)) {
1127: SNESGetNPC(snes, &snes->npc);
1128: }
1129: if (snes->npc) {
1130: SNESSetFromOptions(snes->npc);
1131: }
1132: snes->setfromoptionscalled++;
1133: return(0);
1134: }
1136: /*@
1137: SNESResetFromOptions - Sets various SNES and KSP parameters from user options ONLY if the SNES was previously set from options
1139: Collective on SNES
1141: Input Parameter:
1142: . snes - the SNES context
1144: Level: beginner
1146: .seealso: SNESSetFromOptions(), SNESSetOptionsPrefix()
1147: @*/
1148: PetscErrorCode SNESResetFromOptions(SNES snes)
1149: {
1153: if (snes->setfromoptionscalled) {SNESSetFromOptions(snes);}
1154: return(0);
1155: }
1157: /*@C
1158: SNESSetComputeApplicationContext - Sets an optional function to compute a user-defined context for
1159: the nonlinear solvers.
1161: Logically Collective on SNES
1163: Input Parameters:
1164: + snes - the SNES context
1165: . compute - function to compute the context
1166: - destroy - function to destroy the context
1168: Level: intermediate
1170: Notes:
1171: This function is currently not available from Fortran.
1173: .seealso: SNESGetApplicationContext(), SNESSetComputeApplicationContext(), SNESGetApplicationContext()
1174: @*/
1175: PetscErrorCode SNESSetComputeApplicationContext(SNES snes,PetscErrorCode (*compute)(SNES,void**),PetscErrorCode (*destroy)(void**))
1176: {
1179: snes->ops->usercompute = compute;
1180: snes->ops->userdestroy = destroy;
1181: return(0);
1182: }
1184: /*@
1185: SNESSetApplicationContext - Sets the optional user-defined context for
1186: the nonlinear solvers.
1188: Logically Collective on SNES
1190: Input Parameters:
1191: + snes - the SNES context
1192: - usrP - optional user context
1194: Level: intermediate
1196: Fortran Notes:
1197: To use this from Fortran you must write a Fortran interface definition for this
1198: function that tells Fortran the Fortran derived data type that you are passing in as the ctx argument.
1200: .seealso: SNESGetApplicationContext()
1201: @*/
1202: PetscErrorCode SNESSetApplicationContext(SNES snes,void *usrP)
1203: {
1205: KSP ksp;
1209: SNESGetKSP(snes,&ksp);
1210: KSPSetApplicationContext(ksp,usrP);
1211: snes->user = usrP;
1212: return(0);
1213: }
1215: /*@
1216: SNESGetApplicationContext - Gets the user-defined context for the
1217: nonlinear solvers.
1219: Not Collective
1221: Input Parameter:
1222: . snes - SNES context
1224: Output Parameter:
1225: . usrP - user context
1227: Fortran Notes:
1228: To use this from Fortran you must write a Fortran interface definition for this
1229: function that tells Fortran the Fortran derived data type that you are passing in as the ctx argument.
1231: Level: intermediate
1233: .seealso: SNESSetApplicationContext()
1234: @*/
1235: PetscErrorCode SNESGetApplicationContext(SNES snes,void *usrP)
1236: {
1239: *(void**)usrP = snes->user;
1240: return(0);
1241: }
1243: /*@
1244: SNESSetUseMatrixFree - indicates that SNES should use matrix free finite difference matrix vector products internally to apply the Jacobian.
1246: Collective on SNES
1248: Input Parameters:
1249: + snes - SNES context
1250: . mf_operator - use matrix-free only for the Amat used by SNESSetJacobian(), this means the user provided Pmat will continue to be used
1251: - mf - use matrix-free for both the Amat and Pmat used by SNESSetJacobian(), both the Amat and Pmat set in SNESSetJacobian() will be ignored
1253: Options Database:
1254: + -snes_mf - use matrix free for both the mat and pmat operator
1255: . -snes_mf_operator - use matrix free only for the mat operator
1256: . -snes_fd_color - compute the Jacobian via coloring and finite differences.
1257: - -snes_fd - compute the Jacobian via finite differences (slow)
1259: Level: intermediate
1261: Notes:
1262: SNES supports three approaches for computing (approximate) Jacobians: user provided via SNESSetJacobian(), matrix free, and computing explictly with
1263: finite differences and coloring using MatFDColoring. It is also possible to use automatic differentiation and the MatFDColoring object.
1265: .seealso: SNESGetUseMatrixFree(), MatCreateSNESMF(), SNESComputeJacobianDefaultColor()
1266: @*/
1267: PetscErrorCode SNESSetUseMatrixFree(SNES snes,PetscBool mf_operator,PetscBool mf)
1268: {
1273: snes->mf = mf_operator ? PETSC_TRUE : mf;
1274: snes->mf_operator = mf_operator;
1275: return(0);
1276: }
1278: /*@
1279: SNESGetUseMatrixFree - indicates if the SNES uses matrix free finite difference matrix vector products to apply the Jacobian.
1281: Collective on SNES
1283: Input Parameter:
1284: . snes - SNES context
1286: Output Parameters:
1287: + mf_operator - use matrix-free only for the Amat used by SNESSetJacobian(), this means the user provided Pmat will continue to be used
1288: - mf - use matrix-free for both the Amat and Pmat used by SNESSetJacobian(), both the Amat and Pmat set in SNESSetJacobian() will be ignored
1290: Options Database:
1291: + -snes_mf - use matrix free for both the mat and pmat operator
1292: - -snes_mf_operator - use matrix free only for the mat operator
1294: Level: intermediate
1296: .seealso: SNESSetUseMatrixFree(), MatCreateSNESMF()
1297: @*/
1298: PetscErrorCode SNESGetUseMatrixFree(SNES snes,PetscBool *mf_operator,PetscBool *mf)
1299: {
1302: if (mf) *mf = snes->mf;
1303: if (mf_operator) *mf_operator = snes->mf_operator;
1304: return(0);
1305: }
1307: /*@
1308: SNESGetIterationNumber - Gets the number of nonlinear iterations completed
1309: at this time.
1311: Not Collective
1313: Input Parameter:
1314: . snes - SNES context
1316: Output Parameter:
1317: . iter - iteration number
1319: Notes:
1320: For example, during the computation of iteration 2 this would return 1.
1322: This is useful for using lagged Jacobians (where one does not recompute the
1323: Jacobian at each SNES iteration). For example, the code
1324: .vb
1325: SNESGetIterationNumber(snes,&it);
1326: if (!(it % 2)) {
1327: [compute Jacobian here]
1328: }
1329: .ve
1330: can be used in your ComputeJacobian() function to cause the Jacobian to be
1331: recomputed every second SNES iteration.
1333: After the SNES solve is complete this will return the number of nonlinear iterations used.
1335: Level: intermediate
1337: .seealso: SNESGetLinearSolveIterations()
1338: @*/
1339: PetscErrorCode SNESGetIterationNumber(SNES snes,PetscInt *iter)
1340: {
1344: *iter = snes->iter;
1345: return(0);
1346: }
1348: /*@
1349: SNESSetIterationNumber - Sets the current iteration number.
1351: Not Collective
1353: Input Parameter:
1354: + snes - SNES context
1355: - iter - iteration number
1357: Level: developer
1359: .seealso: SNESGetLinearSolveIterations()
1360: @*/
1361: PetscErrorCode SNESSetIterationNumber(SNES snes,PetscInt iter)
1362: {
1367: PetscObjectSAWsTakeAccess((PetscObject)snes);
1368: snes->iter = iter;
1369: PetscObjectSAWsGrantAccess((PetscObject)snes);
1370: return(0);
1371: }
1373: /*@
1374: SNESGetNonlinearStepFailures - Gets the number of unsuccessful steps
1375: attempted by the nonlinear solver.
1377: Not Collective
1379: Input Parameter:
1380: . snes - SNES context
1382: Output Parameter:
1383: . nfails - number of unsuccessful steps attempted
1385: Notes:
1386: This counter is reset to zero for each successive call to SNESSolve().
1388: Level: intermediate
1390: .seealso: SNESGetMaxLinearSolveFailures(), SNESGetLinearSolveIterations(), SNESSetMaxLinearSolveFailures(), SNESGetLinearSolveFailures(),
1391: SNESSetMaxNonlinearStepFailures(), SNESGetMaxNonlinearStepFailures()
1392: @*/
1393: PetscErrorCode SNESGetNonlinearStepFailures(SNES snes,PetscInt *nfails)
1394: {
1398: *nfails = snes->numFailures;
1399: return(0);
1400: }
1402: /*@
1403: SNESSetMaxNonlinearStepFailures - Sets the maximum number of unsuccessful steps
1404: attempted by the nonlinear solver before it gives up.
1406: Not Collective
1408: Input Parameters:
1409: + snes - SNES context
1410: - maxFails - maximum of unsuccessful steps
1412: Level: intermediate
1414: .seealso: SNESGetMaxLinearSolveFailures(), SNESGetLinearSolveIterations(), SNESSetMaxLinearSolveFailures(), SNESGetLinearSolveFailures(),
1415: SNESGetMaxNonlinearStepFailures(), SNESGetNonlinearStepFailures()
1416: @*/
1417: PetscErrorCode SNESSetMaxNonlinearStepFailures(SNES snes, PetscInt maxFails)
1418: {
1421: snes->maxFailures = maxFails;
1422: return(0);
1423: }
1425: /*@
1426: SNESGetMaxNonlinearStepFailures - Gets the maximum number of unsuccessful steps
1427: attempted by the nonlinear solver before it gives up.
1429: Not Collective
1431: Input Parameter:
1432: . snes - SNES context
1434: Output Parameter:
1435: . maxFails - maximum of unsuccessful steps
1437: Level: intermediate
1439: .seealso: SNESGetMaxLinearSolveFailures(), SNESGetLinearSolveIterations(), SNESSetMaxLinearSolveFailures(), SNESGetLinearSolveFailures(),
1440: SNESSetMaxNonlinearStepFailures(), SNESGetNonlinearStepFailures()
1442: @*/
1443: PetscErrorCode SNESGetMaxNonlinearStepFailures(SNES snes, PetscInt *maxFails)
1444: {
1448: *maxFails = snes->maxFailures;
1449: return(0);
1450: }
1452: /*@
1453: SNESGetNumberFunctionEvals - Gets the number of user provided function evaluations
1454: done by SNES.
1456: Not Collective
1458: Input Parameter:
1459: . snes - SNES context
1461: Output Parameter:
1462: . nfuncs - number of evaluations
1464: Level: intermediate
1466: Notes:
1467: Reset every time SNESSolve is called unless SNESSetCountersReset() is used.
1469: .seealso: SNESGetMaxLinearSolveFailures(), SNESGetLinearSolveIterations(), SNESSetMaxLinearSolveFailures(), SNESGetLinearSolveFailures(), SNESSetCountersReset()
1470: @*/
1471: PetscErrorCode SNESGetNumberFunctionEvals(SNES snes, PetscInt *nfuncs)
1472: {
1476: *nfuncs = snes->nfuncs;
1477: return(0);
1478: }
1480: /*@
1481: SNESGetLinearSolveFailures - Gets the number of failed (non-converged)
1482: linear solvers.
1484: Not Collective
1486: Input Parameter:
1487: . snes - SNES context
1489: Output Parameter:
1490: . nfails - number of failed solves
1492: Level: intermediate
1494: Options Database Keys:
1495: . -snes_max_linear_solve_fail <num> - The number of failures before the solve is terminated
1497: Notes:
1498: This counter is reset to zero for each successive call to SNESSolve().
1500: .seealso: SNESGetMaxLinearSolveFailures(), SNESGetLinearSolveIterations(), SNESSetMaxLinearSolveFailures()
1501: @*/
1502: PetscErrorCode SNESGetLinearSolveFailures(SNES snes,PetscInt *nfails)
1503: {
1507: *nfails = snes->numLinearSolveFailures;
1508: return(0);
1509: }
1511: /*@
1512: SNESSetMaxLinearSolveFailures - the number of failed linear solve attempts
1513: allowed before SNES returns with a diverged reason of SNES_DIVERGED_LINEAR_SOLVE
1515: Logically Collective on SNES
1517: Input Parameters:
1518: + snes - SNES context
1519: - maxFails - maximum allowed linear solve failures
1521: Level: intermediate
1523: Options Database Keys:
1524: . -snes_max_linear_solve_fail <num> - The number of failures before the solve is terminated
1526: Notes:
1527: By default this is 0; that is SNES returns on the first failed linear solve
1529: .seealso: SNESGetLinearSolveFailures(), SNESGetMaxLinearSolveFailures(), SNESGetLinearSolveIterations()
1530: @*/
1531: PetscErrorCode SNESSetMaxLinearSolveFailures(SNES snes, PetscInt maxFails)
1532: {
1536: snes->maxLinearSolveFailures = maxFails;
1537: return(0);
1538: }
1540: /*@
1541: SNESGetMaxLinearSolveFailures - gets the maximum number of linear solve failures that
1542: are allowed before SNES terminates
1544: Not Collective
1546: Input Parameter:
1547: . snes - SNES context
1549: Output Parameter:
1550: . maxFails - maximum of unsuccessful solves allowed
1552: Level: intermediate
1554: Notes:
1555: By default this is 1; that is SNES returns on the first failed linear solve
1557: .seealso: SNESGetLinearSolveFailures(), SNESGetLinearSolveIterations(), SNESSetMaxLinearSolveFailures(),
1558: @*/
1559: PetscErrorCode SNESGetMaxLinearSolveFailures(SNES snes, PetscInt *maxFails)
1560: {
1564: *maxFails = snes->maxLinearSolveFailures;
1565: return(0);
1566: }
1568: /*@
1569: SNESGetLinearSolveIterations - Gets the total number of linear iterations
1570: used by the nonlinear solver.
1572: Not Collective
1574: Input Parameter:
1575: . snes - SNES context
1577: Output Parameter:
1578: . lits - number of linear iterations
1580: Notes:
1581: This counter is reset to zero for each successive call to SNESSolve() unless SNESSetCountersReset() is used.
1583: 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
1584: then call KSPGetIterationNumber() after the failed solve.
1586: Level: intermediate
1588: .seealso: SNESGetIterationNumber(), SNESGetLinearSolveFailures(), SNESGetMaxLinearSolveFailures(), SNESSetCountersReset()
1589: @*/
1590: PetscErrorCode SNESGetLinearSolveIterations(SNES snes,PetscInt *lits)
1591: {
1595: *lits = snes->linear_its;
1596: return(0);
1597: }
1599: /*@
1600: SNESSetCountersReset - Sets whether or not the counters for linear iterations and function evaluations
1601: are reset every time SNESSolve() is called.
1603: Logically Collective on SNES
1605: Input Parameter:
1606: + snes - SNES context
1607: - reset - whether to reset the counters or not
1609: Notes:
1610: This defaults to PETSC_TRUE
1612: Level: developer
1614: .seealso: SNESGetNumberFunctionEvals(), SNESGetLinearSolveIterations(), SNESGetNPC()
1615: @*/
1616: PetscErrorCode SNESSetCountersReset(SNES snes,PetscBool reset)
1617: {
1621: snes->counters_reset = reset;
1622: return(0);
1623: }
1626: /*@
1627: SNESSetKSP - Sets a KSP context for the SNES object to use
1629: Not Collective, but the SNES and KSP objects must live on the same MPI_Comm
1631: Input Parameters:
1632: + snes - the SNES context
1633: - ksp - the KSP context
1635: Notes:
1636: The SNES object already has its KSP object, you can obtain with SNESGetKSP()
1637: so this routine is rarely needed.
1639: The KSP object that is already in the SNES object has its reference count
1640: decreased by one.
1642: Level: developer
1644: .seealso: KSPGetPC(), SNESCreate(), KSPCreate(), SNESSetKSP()
1645: @*/
1646: PetscErrorCode SNESSetKSP(SNES snes,KSP ksp)
1647: {
1654: PetscObjectReference((PetscObject)ksp);
1655: if (snes->ksp) {PetscObjectDereference((PetscObject)snes->ksp);}
1656: snes->ksp = ksp;
1657: return(0);
1658: }
1660: /* -----------------------------------------------------------*/
1661: /*@
1662: SNESCreate - Creates a nonlinear solver context.
1664: Collective
1666: Input Parameters:
1667: . comm - MPI communicator
1669: Output Parameter:
1670: . outsnes - the new SNES context
1672: Options Database Keys:
1673: + -snes_mf - Activates default matrix-free Jacobian-vector products,
1674: and no preconditioning matrix
1675: . -snes_mf_operator - Activates default matrix-free Jacobian-vector
1676: products, and a user-provided preconditioning matrix
1677: as set by SNESSetJacobian()
1678: - -snes_fd - Uses (slow!) finite differences to compute Jacobian
1680: Level: beginner
1682: Developer Notes:
1683: SNES always creates a KSP object even though many SNES methods do not use it. This is
1684: unfortunate and should be fixed at some point. The flag snes->usesksp indicates if the
1685: particular method does use KSP and regulates if the information about the KSP is printed
1686: in SNESView(). TSSetFromOptions() does call SNESSetFromOptions() which can lead to users being confused
1687: by help messages about meaningless SNES options.
1689: SNES always creates the snes->kspconvctx even though it is used by only one type. This should
1690: be fixed.
1692: .seealso: SNESSolve(), SNESDestroy(), SNES, SNESSetLagPreconditioner(), SNESSetLagJacobian()
1694: @*/
1695: PetscErrorCode SNESCreate(MPI_Comm comm,SNES *outsnes)
1696: {
1698: SNES snes;
1699: SNESKSPEW *kctx;
1703: *outsnes = NULL;
1704: SNESInitializePackage();
1706: PetscHeaderCreate(snes,SNES_CLASSID,"SNES","Nonlinear solver","SNES",comm,SNESDestroy,SNESView);
1708: snes->ops->converged = SNESConvergedDefault;
1709: snes->usesksp = PETSC_TRUE;
1710: snes->tolerancesset = PETSC_FALSE;
1711: snes->max_its = 50;
1712: snes->max_funcs = 10000;
1713: snes->norm = 0.0;
1714: snes->xnorm = 0.0;
1715: snes->ynorm = 0.0;
1716: snes->normschedule = SNES_NORM_ALWAYS;
1717: snes->functype = SNES_FUNCTION_DEFAULT;
1718: #if defined(PETSC_USE_REAL_SINGLE)
1719: snes->rtol = 1.e-5;
1720: #else
1721: snes->rtol = 1.e-8;
1722: #endif
1723: snes->ttol = 0.0;
1724: #if defined(PETSC_USE_REAL_SINGLE)
1725: snes->abstol = 1.e-25;
1726: #else
1727: snes->abstol = 1.e-50;
1728: #endif
1729: #if defined(PETSC_USE_REAL_SINGLE)
1730: snes->stol = 1.e-5;
1731: #else
1732: snes->stol = 1.e-8;
1733: #endif
1734: #if defined(PETSC_USE_REAL_SINGLE)
1735: snes->deltatol = 1.e-6;
1736: #else
1737: snes->deltatol = 1.e-12;
1738: #endif
1739: snes->divtol = 1.e4;
1740: snes->rnorm0 = 0;
1741: snes->nfuncs = 0;
1742: snes->numFailures = 0;
1743: snes->maxFailures = 1;
1744: snes->linear_its = 0;
1745: snes->lagjacobian = 1;
1746: snes->jac_iter = 0;
1747: snes->lagjac_persist = PETSC_FALSE;
1748: snes->lagpreconditioner = 1;
1749: snes->pre_iter = 0;
1750: snes->lagpre_persist = PETSC_FALSE;
1751: snes->numbermonitors = 0;
1752: snes->data = NULL;
1753: snes->setupcalled = PETSC_FALSE;
1754: snes->ksp_ewconv = PETSC_FALSE;
1755: snes->nwork = 0;
1756: snes->work = NULL;
1757: snes->nvwork = 0;
1758: snes->vwork = NULL;
1759: snes->conv_hist_len = 0;
1760: snes->conv_hist_max = 0;
1761: snes->conv_hist = NULL;
1762: snes->conv_hist_its = NULL;
1763: snes->conv_hist_reset = PETSC_TRUE;
1764: snes->counters_reset = PETSC_TRUE;
1765: snes->vec_func_init_set = PETSC_FALSE;
1766: snes->reason = SNES_CONVERGED_ITERATING;
1767: snes->npcside = PC_RIGHT;
1768: snes->setfromoptionscalled = 0;
1770: snes->mf = PETSC_FALSE;
1771: snes->mf_operator = PETSC_FALSE;
1772: snes->mf_version = 1;
1774: snes->numLinearSolveFailures = 0;
1775: snes->maxLinearSolveFailures = 1;
1777: snes->vizerotolerance = 1.e-8;
1778: snes->checkjacdomainerror = PetscDefined(USE_DEBUG) ? PETSC_TRUE : PETSC_FALSE;
1780: /* Set this to true if the implementation of SNESSolve_XXX does compute the residual at the final solution. */
1781: snes->alwayscomputesfinalresidual = PETSC_FALSE;
1783: /* Create context to compute Eisenstat-Walker relative tolerance for KSP */
1784: PetscNewLog(snes,&kctx);
1786: snes->kspconvctx = (void*)kctx;
1787: kctx->version = 2;
1788: kctx->rtol_0 = .3; /* Eisenstat and Walker suggest rtol_0=.5, but
1789: this was too large for some test cases */
1790: kctx->rtol_last = 0.0;
1791: kctx->rtol_max = .9;
1792: kctx->gamma = 1.0;
1793: kctx->alpha = .5*(1.0 + PetscSqrtReal(5.0));
1794: kctx->alpha2 = kctx->alpha;
1795: kctx->threshold = .1;
1796: kctx->lresid_last = 0.0;
1797: kctx->norm_last = 0.0;
1799: *outsnes = snes;
1800: return(0);
1801: }
1803: /*MC
1804: SNESFunction - Functional form used to convey the nonlinear function to be solved by SNES
1806: Synopsis:
1807: #include "petscsnes.h"
1808: PetscErrorCode SNESFunction(SNES snes,Vec x,Vec f,void *ctx);
1810: Collective on snes
1812: Input Parameters:
1813: + snes - the SNES context
1814: . x - state at which to evaluate residual
1815: - ctx - optional user-defined function context, passed in with SNESSetFunction()
1817: Output Parameter:
1818: . f - vector to put residual (function value)
1820: Level: intermediate
1822: .seealso: SNESSetFunction(), SNESGetFunction()
1823: M*/
1825: /*@C
1826: SNESSetFunction - Sets the function evaluation routine and function
1827: vector for use by the SNES routines in solving systems of nonlinear
1828: equations.
1830: Logically Collective on SNES
1832: Input Parameters:
1833: + snes - the SNES context
1834: . r - vector to store function value
1835: . f - function evaluation routine; see SNESFunction for calling sequence details
1836: - ctx - [optional] user-defined context for private data for the
1837: function evaluation routine (may be NULL)
1839: Notes:
1840: The Newton-like methods typically solve linear systems of the form
1841: $ f'(x) x = -f(x),
1842: where f'(x) denotes the Jacobian matrix and f(x) is the function.
1844: Level: beginner
1846: .seealso: SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESSetPicard(), SNESFunction
1847: @*/
1848: PetscErrorCode SNESSetFunction(SNES snes,Vec r,PetscErrorCode (*f)(SNES,Vec,Vec,void*),void *ctx)
1849: {
1851: DM dm;
1855: if (r) {
1858: PetscObjectReference((PetscObject)r);
1859: VecDestroy(&snes->vec_func);
1861: snes->vec_func = r;
1862: }
1863: SNESGetDM(snes,&dm);
1864: DMSNESSetFunction(dm,f,ctx);
1865: return(0);
1866: }
1869: /*@C
1870: SNESSetInitialFunction - Sets the function vector to be used as the
1871: function norm at the initialization of the method. In some
1872: instances, the user has precomputed the function before calling
1873: SNESSolve. This function allows one to avoid a redundant call
1874: to SNESComputeFunction in that case.
1876: Logically Collective on SNES
1878: Input Parameters:
1879: + snes - the SNES context
1880: - f - vector to store function value
1882: Notes:
1883: This should not be modified during the solution procedure.
1885: This is used extensively in the SNESFAS hierarchy and in nonlinear preconditioning.
1887: Level: developer
1889: .seealso: SNESSetFunction(), SNESComputeFunction(), SNESSetInitialFunctionNorm()
1890: @*/
1891: PetscErrorCode SNESSetInitialFunction(SNES snes, Vec f)
1892: {
1894: Vec vec_func;
1900: if (snes->npcside== PC_LEFT && snes->functype == SNES_FUNCTION_PRECONDITIONED) {
1901: snes->vec_func_init_set = PETSC_FALSE;
1902: return(0);
1903: }
1904: SNESGetFunction(snes,&vec_func,NULL,NULL);
1905: VecCopy(f, vec_func);
1907: snes->vec_func_init_set = PETSC_TRUE;
1908: return(0);
1909: }
1911: /*@
1912: SNESSetNormSchedule - Sets the SNESNormSchedule used in covergence and monitoring
1913: of the SNES method.
1915: Logically Collective on SNES
1917: Input Parameters:
1918: + snes - the SNES context
1919: - normschedule - the frequency of norm computation
1921: Options Database Key:
1922: . -snes_norm_schedule <none, always, initialonly, finalonly, initalfinalonly>
1924: Notes:
1925: Only certain SNES methods support certain SNESNormSchedules. Most require evaluation
1926: of the nonlinear function and the taking of its norm at every iteration to
1927: even ensure convergence at all. However, methods such as custom Gauss-Seidel methods
1928: (SNESNGS) and the like do not require the norm of the function to be computed, and therfore
1929: may either be monitored for convergence or not. As these are often used as nonlinear
1930: preconditioners, monitoring the norm of their error is not a useful enterprise within
1931: their solution.
1933: Level: developer
1935: .seealso: SNESGetNormSchedule(), SNESComputeFunction(), VecNorm(), SNESSetFunction(), SNESSetInitialFunction(), SNESNormSchedule
1936: @*/
1937: PetscErrorCode SNESSetNormSchedule(SNES snes, SNESNormSchedule normschedule)
1938: {
1941: snes->normschedule = normschedule;
1942: return(0);
1943: }
1946: /*@
1947: SNESGetNormSchedule - Gets the SNESNormSchedule used in covergence and monitoring
1948: of the SNES method.
1950: Logically Collective on SNES
1952: Input Parameters:
1953: + snes - the SNES context
1954: - normschedule - the type of the norm used
1956: Level: advanced
1958: .seealso: SNESSetNormSchedule(), SNESComputeFunction(), VecNorm(), SNESSetFunction(), SNESSetInitialFunction(), SNESNormSchedule
1959: @*/
1960: PetscErrorCode SNESGetNormSchedule(SNES snes, SNESNormSchedule *normschedule)
1961: {
1964: *normschedule = snes->normschedule;
1965: return(0);
1966: }
1969: /*@
1970: SNESSetFunctionNorm - Sets the last computed residual norm.
1972: Logically Collective on SNES
1974: Input Parameters:
1975: + snes - the SNES context
1977: - normschedule - the frequency of norm computation
1979: Level: developer
1981: .seealso: SNESGetNormSchedule(), SNESComputeFunction(), VecNorm(), SNESSetFunction(), SNESSetInitialFunction(), SNESNormSchedule
1982: @*/
1983: PetscErrorCode SNESSetFunctionNorm(SNES snes, PetscReal norm)
1984: {
1987: snes->norm = norm;
1988: return(0);
1989: }
1991: /*@
1992: SNESGetFunctionNorm - Gets the last computed norm of the residual
1994: Not Collective
1996: Input Parameter:
1997: . snes - the SNES context
1999: Output Parameter:
2000: . norm - the last computed residual norm
2002: Level: developer
2004: .seealso: SNESSetNormSchedule(), SNESComputeFunction(), VecNorm(), SNESSetFunction(), SNESSetInitialFunction(), SNESNormSchedule
2005: @*/
2006: PetscErrorCode SNESGetFunctionNorm(SNES snes, PetscReal *norm)
2007: {
2011: *norm = snes->norm;
2012: return(0);
2013: }
2015: /*@
2016: SNESGetUpdateNorm - Gets the last computed norm of the Newton update
2018: Not Collective
2020: Input Parameter:
2021: . snes - the SNES context
2023: Output Parameter:
2024: . ynorm - the last computed update norm
2026: Level: developer
2028: .seealso: SNESSetNormSchedule(), SNESComputeFunction(), SNESGetFunctionNorm()
2029: @*/
2030: PetscErrorCode SNESGetUpdateNorm(SNES snes, PetscReal *ynorm)
2031: {
2035: *ynorm = snes->ynorm;
2036: return(0);
2037: }
2039: /*@
2040: SNESGetSolutionNorm - Gets the last computed norm of the solution
2042: Not Collective
2044: Input Parameter:
2045: . snes - the SNES context
2047: Output Parameter:
2048: . xnorm - the last computed solution norm
2050: Level: developer
2052: .seealso: SNESSetNormSchedule(), SNESComputeFunction(), SNESGetFunctionNorm(), SNESGetUpdateNorm()
2053: @*/
2054: PetscErrorCode SNESGetSolutionNorm(SNES snes, PetscReal *xnorm)
2055: {
2059: *xnorm = snes->xnorm;
2060: return(0);
2061: }
2063: /*@C
2064: SNESSetFunctionType - Sets the SNESNormSchedule used in covergence and monitoring
2065: of the SNES method.
2067: Logically Collective on SNES
2069: Input Parameters:
2070: + snes - the SNES context
2071: - normschedule - the frequency of norm computation
2073: Notes:
2074: Only certain SNES methods support certain SNESNormSchedules. Most require evaluation
2075: of the nonlinear function and the taking of its norm at every iteration to
2076: even ensure convergence at all. However, methods such as custom Gauss-Seidel methods
2077: (SNESNGS) and the like do not require the norm of the function to be computed, and therfore
2078: may either be monitored for convergence or not. As these are often used as nonlinear
2079: preconditioners, monitoring the norm of their error is not a useful enterprise within
2080: their solution.
2082: Level: developer
2084: .seealso: SNESGetNormSchedule(), SNESComputeFunction(), VecNorm(), SNESSetFunction(), SNESSetInitialFunction(), SNESNormSchedule
2085: @*/
2086: PetscErrorCode SNESSetFunctionType(SNES snes, SNESFunctionType type)
2087: {
2090: snes->functype = type;
2091: return(0);
2092: }
2095: /*@C
2096: SNESGetFunctionType - Gets the SNESNormSchedule used in covergence and monitoring
2097: of the SNES method.
2099: Logically Collective on SNES
2101: Input Parameters:
2102: + snes - the SNES context
2103: - normschedule - the type of the norm used
2105: Level: advanced
2107: .seealso: SNESSetNormSchedule(), SNESComputeFunction(), VecNorm(), SNESSetFunction(), SNESSetInitialFunction(), SNESNormSchedule
2108: @*/
2109: PetscErrorCode SNESGetFunctionType(SNES snes, SNESFunctionType *type)
2110: {
2113: *type = snes->functype;
2114: return(0);
2115: }
2117: /*MC
2118: SNESNGSFunction - function used to convey a Gauss-Seidel sweep on the nonlinear function
2120: Synopsis:
2121: #include <petscsnes.h>
2122: $ SNESNGSFunction(SNES snes,Vec x,Vec b,void *ctx);
2124: Collective on snes
2126: Input Parameters:
2127: + X - solution vector
2128: . B - RHS vector
2129: - ctx - optional user-defined Gauss-Seidel context
2131: Output Parameter:
2132: . X - solution vector
2134: Level: intermediate
2136: .seealso: SNESSetNGS(), SNESGetNGS()
2137: M*/
2139: /*@C
2140: SNESSetNGS - Sets the user nonlinear Gauss-Seidel routine for
2141: use with composed nonlinear solvers.
2143: Input Parameters:
2144: + snes - the SNES context
2145: . f - function evaluation routine to apply Gauss-Seidel see SNESNGSFunction
2146: - ctx - [optional] user-defined context for private data for the
2147: smoother evaluation routine (may be NULL)
2149: Notes:
2150: The NGS routines are used by the composed nonlinear solver to generate
2151: a problem appropriate update to the solution, particularly FAS.
2153: Level: intermediate
2155: .seealso: SNESGetFunction(), SNESComputeNGS()
2156: @*/
2157: PetscErrorCode SNESSetNGS(SNES snes,PetscErrorCode (*f)(SNES,Vec,Vec,void*),void *ctx)
2158: {
2160: DM dm;
2164: SNESGetDM(snes,&dm);
2165: DMSNESSetNGS(dm,f,ctx);
2166: return(0);
2167: }
2169: PetscErrorCode SNESPicardComputeFunction(SNES snes,Vec x,Vec f,void *ctx)
2170: {
2172: DM dm;
2173: DMSNES sdm;
2176: SNESGetDM(snes,&dm);
2177: DMGetDMSNES(dm,&sdm);
2178: if (!sdm->ops->computepfunction) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE, "Must call SNESSetPicard() to provide Picard function.");
2179: if (!sdm->ops->computepjacobian) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE, "Must call SNESSetPicard() to provide Picard Jacobian.");
2180: /* A(x)*x - b(x) */
2181: PetscStackPush("SNES Picard user function");
2182: (*sdm->ops->computepfunction)(snes,x,f,sdm->pctx);
2183: PetscStackPop;
2184: PetscStackPush("SNES Picard user Jacobian");
2185: (*sdm->ops->computepjacobian)(snes,x,snes->jacobian,snes->jacobian_pre,sdm->pctx);
2186: PetscStackPop;
2187: VecScale(f,-1.0);
2188: MatMultAdd(snes->jacobian,x,f,f);
2189: return(0);
2190: }
2192: PetscErrorCode SNESPicardComputeJacobian(SNES snes,Vec x1,Mat J,Mat B,void *ctx)
2193: {
2195: /* the jacobian matrix should be pre-filled in SNESPicardComputeFunction */
2196: return(0);
2197: }
2199: /*@C
2200: SNESSetPicard - Use SNES to solve the semilinear-system A(x) x = b(x) via a Picard type iteration (Picard linearization)
2202: Logically Collective on SNES
2204: Input Parameters:
2205: + snes - the SNES context
2206: . r - vector to store function value
2207: . b - function evaluation routine
2208: . Amat - matrix with which A(x) x - b(x) is to be computed
2209: . Pmat - matrix from which preconditioner is computed (usually the same as Amat)
2210: . J - function to compute matrix value, see SNESJacobianFunction for details on its calling sequence
2211: - ctx - [optional] user-defined context for private data for the
2212: function evaluation routine (may be NULL)
2214: Notes:
2215: We do not recomemend using this routine. It is far better to provide the nonlinear function F() and some approximation to the Jacobian and use
2216: 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.
2218: One can call SNESSetPicard() or SNESSetFunction() (and possibly SNESSetJacobian()) but cannot call both
2220: $ Solves the equation A(x) x = b(x) via the defect correction algorithm A(x^{n}) (x^{n+1} - x^{n}) = b(x^{n}) - A(x^{n})x^{n}
2221: $ Note that when an exact solver is used this corresponds to the "classic" Picard A(x^{n}) x^{n+1} = b(x^{n}) iteration.
2223: Run with -snes_mf_operator to solve the system with Newton's method using A(x^{n}) to construct the preconditioner.
2225: We implement the defect correction form of the Picard iteration because it converges much more generally when inexact linear solvers are used then
2226: the direct Picard iteration A(x^n) x^{n+1} = b(x^n)
2228: 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
2229: 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
2230: different please contact us at petsc-dev@mcs.anl.gov and we'll have an entirely new argument :-).
2232: Level: intermediate
2234: .seealso: SNESGetFunction(), SNESSetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESGetPicard(), SNESLineSearchPreCheckPicard(), SNESJacobianFunction
2235: @*/
2236: PetscErrorCode SNESSetPicard(SNES snes,Vec r,PetscErrorCode (*b)(SNES,Vec,Vec,void*),Mat Amat, Mat Pmat, PetscErrorCode (*J)(SNES,Vec,Mat,Mat,void*),void *ctx)
2237: {
2239: DM dm;
2243: SNESGetDM(snes, &dm);
2244: DMSNESSetPicard(dm,b,J,ctx);
2245: SNESSetFunction(snes,r,SNESPicardComputeFunction,ctx);
2246: SNESSetJacobian(snes,Amat,Pmat,SNESPicardComputeJacobian,ctx);
2247: return(0);
2248: }
2250: /*@C
2251: SNESGetPicard - Returns the context for the Picard iteration
2253: Not Collective, but Vec is parallel if SNES is parallel. Collective if Vec is requested, but has not been created yet.
2255: Input Parameter:
2256: . snes - the SNES context
2258: Output Parameter:
2259: + r - the function (or NULL)
2260: . f - the function (or NULL); see SNESFunction for calling sequence details
2261: . Amat - the matrix used to defined the operation A(x) x - b(x) (or NULL)
2262: . Pmat - the matrix from which the preconditioner will be constructed (or NULL)
2263: . J - the function for matrix evaluation (or NULL); see SNESJacobianFunction for calling sequence details
2264: - ctx - the function context (or NULL)
2266: Level: advanced
2268: .seealso: SNESSetPicard(), SNESGetFunction(), SNESGetJacobian(), SNESGetDM(), SNESFunction, SNESJacobianFunction
2269: @*/
2270: 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)
2271: {
2273: DM dm;
2277: SNESGetFunction(snes,r,NULL,NULL);
2278: SNESGetJacobian(snes,Amat,Pmat,NULL,NULL);
2279: SNESGetDM(snes,&dm);
2280: DMSNESGetPicard(dm,f,J,ctx);
2281: return(0);
2282: }
2284: /*@C
2285: SNESSetComputeInitialGuess - Sets a routine used to compute an initial guess for the problem
2287: Logically Collective on SNES
2289: Input Parameters:
2290: + snes - the SNES context
2291: . func - function evaluation routine
2292: - ctx - [optional] user-defined context for private data for the
2293: function evaluation routine (may be NULL)
2295: Calling sequence of func:
2296: $ func (SNES snes,Vec x,void *ctx);
2298: . f - function vector
2299: - ctx - optional user-defined function context
2301: Level: intermediate
2303: .seealso: SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian()
2304: @*/
2305: PetscErrorCode SNESSetComputeInitialGuess(SNES snes,PetscErrorCode (*func)(SNES,Vec,void*),void *ctx)
2306: {
2309: if (func) snes->ops->computeinitialguess = func;
2310: if (ctx) snes->initialguessP = ctx;
2311: return(0);
2312: }
2314: /* --------------------------------------------------------------- */
2315: /*@C
2316: SNESGetRhs - Gets the vector for solving F(x) = rhs. If rhs is not set
2317: it assumes a zero right hand side.
2319: Logically Collective on SNES
2321: Input Parameter:
2322: . snes - the SNES context
2324: Output Parameter:
2325: . rhs - the right hand side vector or NULL if the right hand side vector is null
2327: Level: intermediate
2329: .seealso: SNESGetSolution(), SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESSetFunction()
2330: @*/
2331: PetscErrorCode SNESGetRhs(SNES snes,Vec *rhs)
2332: {
2336: *rhs = snes->vec_rhs;
2337: return(0);
2338: }
2340: /*@
2341: SNESComputeFunction - Calls the function that has been set with SNESSetFunction().
2343: Collective on SNES
2345: Input Parameters:
2346: + snes - the SNES context
2347: - x - input vector
2349: Output Parameter:
2350: . y - function vector, as set by SNESSetFunction()
2352: Notes:
2353: SNESComputeFunction() is typically used within nonlinear solvers
2354: implementations, so most users would not generally call this routine
2355: themselves.
2357: Level: developer
2359: .seealso: SNESSetFunction(), SNESGetFunction()
2360: @*/
2361: PetscErrorCode SNESComputeFunction(SNES snes,Vec x,Vec y)
2362: {
2364: DM dm;
2365: DMSNES sdm;
2373: VecValidValues(x,2,PETSC_TRUE);
2375: SNESGetDM(snes,&dm);
2376: DMGetDMSNES(dm,&sdm);
2377: if (sdm->ops->computefunction) {
2378: if (sdm->ops->computefunction != SNESObjectiveComputeFunctionDefaultFD) {
2379: PetscLogEventBegin(SNES_FunctionEval,snes,x,y,0);
2380: }
2381: VecLockReadPush(x);
2382: PetscStackPush("SNES user function");
2383: /* ensure domainerror is false prior to computefunction evaluation (may not have been reset) */
2384: snes->domainerror = PETSC_FALSE;
2385: (*sdm->ops->computefunction)(snes,x,y,sdm->functionctx);
2386: PetscStackPop;
2387: VecLockReadPop(x);
2388: if (sdm->ops->computefunction != SNESObjectiveComputeFunctionDefaultFD) {
2389: PetscLogEventEnd(SNES_FunctionEval,snes,x,y,0);
2390: }
2391: } else if (snes->vec_rhs) {
2392: MatMult(snes->jacobian, x, y);
2393: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE, "Must call SNESSetFunction() or SNESSetDM() before SNESComputeFunction(), likely called from SNESSolve().");
2394: if (snes->vec_rhs) {
2395: VecAXPY(y,-1.0,snes->vec_rhs);
2396: }
2397: snes->nfuncs++;
2398: /*
2399: domainerror might not be set on all processes; so we tag vector locally with Inf and the next inner product or norm will
2400: propagate the value to all processes
2401: */
2402: if (snes->domainerror) {
2403: VecSetInf(y);
2404: }
2405: return(0);
2406: }
2408: /*@
2409: SNESComputeNGS - Calls the Gauss-Seidel function that has been set with SNESSetNGS().
2411: Collective on SNES
2413: Input Parameters:
2414: + snes - the SNES context
2415: . x - input vector
2416: - b - rhs vector
2418: Output Parameter:
2419: . x - new solution vector
2421: Notes:
2422: SNESComputeNGS() is typically used within composed nonlinear solver
2423: implementations, so most users would not generally call this routine
2424: themselves.
2426: Level: developer
2428: .seealso: SNESSetNGS(), SNESComputeFunction()
2429: @*/
2430: PetscErrorCode SNESComputeNGS(SNES snes,Vec b,Vec x)
2431: {
2433: DM dm;
2434: DMSNES sdm;
2442: if (b) {VecValidValues(b,2,PETSC_TRUE);}
2443: PetscLogEventBegin(SNES_NGSEval,snes,x,b,0);
2444: SNESGetDM(snes,&dm);
2445: DMGetDMSNES(dm,&sdm);
2446: if (sdm->ops->computegs) {
2447: if (b) {VecLockReadPush(b);}
2448: PetscStackPush("SNES user NGS");
2449: (*sdm->ops->computegs)(snes,x,b,sdm->gsctx);
2450: PetscStackPop;
2451: if (b) {VecLockReadPop(b);}
2452: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE, "Must call SNESSetNGS() before SNESComputeNGS(), likely called from SNESSolve().");
2453: PetscLogEventEnd(SNES_NGSEval,snes,x,b,0);
2454: return(0);
2455: }
2457: PetscErrorCode SNESTestJacobian(SNES snes)
2458: {
2459: Mat A,B,C,D,jacobian;
2460: Vec x = snes->vec_sol,f = snes->vec_func;
2461: PetscErrorCode ierr;
2462: PetscReal nrm,gnorm;
2463: PetscReal threshold = 1.e-5;
2464: MatType mattype;
2465: PetscInt m,n,M,N;
2466: void *functx;
2467: PetscBool complete_print = PETSC_FALSE,threshold_print = PETSC_FALSE,test = PETSC_FALSE,flg,istranspose;
2468: PetscViewer viewer,mviewer;
2469: MPI_Comm comm;
2470: PetscInt tabs;
2471: static PetscBool directionsprinted = PETSC_FALSE;
2472: PetscViewerFormat format;
2475: PetscObjectOptionsBegin((PetscObject)snes);
2476: PetscOptionsName("-snes_test_jacobian","Compare hand-coded and finite difference Jacobians","None",&test);
2477: PetscOptionsReal("-snes_test_jacobian", "Threshold for element difference between hand-coded and finite difference being meaningful", "None", threshold, &threshold,NULL);
2478: PetscOptionsViewer("-snes_test_jacobian_view","View difference between hand-coded and finite difference Jacobians element entries","None",&mviewer,&format,&complete_print);
2479: if (!complete_print) {
2480: PetscOptionsDeprecated("-snes_test_jacobian_display","-snes_test_jacobian_view","3.13",NULL);
2481: PetscOptionsViewer("-snes_test_jacobian_display","Display difference between hand-coded and finite difference Jacobians","None",&mviewer,&format,&complete_print);
2482: }
2483: /* for compatibility with PETSc 3.9 and older. */
2484: PetscOptionsDeprecated("-snes_test_jacobian_display_threshold","-snes_test_jacobian","3.13","-snes_test_jacobian accepts an optional threshold (since v3.10)");
2485: PetscOptionsReal("-snes_test_jacobian_display_threshold", "Display difference between hand-coded and finite difference Jacobians which exceed input threshold", "None", threshold, &threshold, &threshold_print);
2486: PetscOptionsEnd();
2487: if (!test) return(0);
2489: PetscObjectGetComm((PetscObject)snes,&comm);
2490: PetscViewerASCIIGetStdout(comm,&viewer);
2491: PetscViewerASCIIGetTab(viewer, &tabs);
2492: PetscViewerASCIISetTab(viewer, ((PetscObject)snes)->tablevel);
2493: PetscViewerASCIIPrintf(viewer," ---------- Testing Jacobian -------------\n");
2494: if (!complete_print && !directionsprinted) {
2495: PetscViewerASCIIPrintf(viewer," Run with -snes_test_jacobian_view and optionally -snes_test_jacobian <threshold> to show difference\n");
2496: PetscViewerASCIIPrintf(viewer," of hand-coded and finite difference Jacobian entries greater than <threshold>.\n");
2497: }
2498: if (!directionsprinted) {
2499: PetscViewerASCIIPrintf(viewer," Testing hand-coded Jacobian, if (for double precision runs) ||J - Jfd||_F/||J||_F is\n");
2500: PetscViewerASCIIPrintf(viewer," O(1.e-8), the hand-coded Jacobian is probably correct.\n");
2501: directionsprinted = PETSC_TRUE;
2502: }
2503: if (complete_print) {
2504: PetscViewerPushFormat(mviewer,format);
2505: }
2507: PetscObjectTypeCompare((PetscObject)snes->jacobian,MATMFFD,&flg);
2508: if (!flg) jacobian = snes->jacobian;
2509: else jacobian = snes->jacobian_pre;
2511: if (!x) {
2512: MatCreateVecs(jacobian, &x, NULL);
2513: } else {
2514: PetscObjectReference((PetscObject) x);
2515: }
2516: if (!f) {
2517: VecDuplicate(x, &f);
2518: } else {
2519: PetscObjectReference((PetscObject) f);
2520: }
2521: /* evaluate the function at this point because SNESComputeJacobianDefault() assumes that the function has been evaluated and put into snes->vec_func */
2522: SNESComputeFunction(snes,x,f);
2523: VecDestroy(&f);
2524: PetscObjectTypeCompare((PetscObject)snes,SNESKSPTRANSPOSEONLY,&istranspose);
2525: while (jacobian) {
2526: Mat JT = NULL, Jsave = NULL;
2528: if (istranspose) {
2529: MatCreateTranspose(jacobian,&JT);
2530: Jsave = jacobian;
2531: jacobian = JT;
2532: }
2533: PetscObjectBaseTypeCompareAny((PetscObject)jacobian,&flg,MATSEQAIJ,MATMPIAIJ,MATSEQDENSE,MATMPIDENSE,MATSEQBAIJ,MATMPIBAIJ,MATSEQSBAIJ,MATMPISBAIJ,"");
2534: if (flg) {
2535: A = jacobian;
2536: PetscObjectReference((PetscObject)A);
2537: } else {
2538: MatComputeOperator(jacobian,MATAIJ,&A);
2539: }
2541: MatGetType(A,&mattype);
2542: MatGetSize(A,&M,&N);
2543: MatGetLocalSize(A,&m,&n);
2544: MatCreate(PetscObjectComm((PetscObject)A),&B);
2545: MatSetType(B,mattype);
2546: MatSetSizes(B,m,n,M,N);
2547: MatSetBlockSizesFromMats(B,A,A);
2548: MatSetUp(B);
2549: MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);
2551: SNESGetFunction(snes,NULL,NULL,&functx);
2552: SNESComputeJacobianDefault(snes,x,B,B,functx);
2554: MatDuplicate(B,MAT_COPY_VALUES,&D);
2555: MatAYPX(D,-1.0,A,DIFFERENT_NONZERO_PATTERN);
2556: MatNorm(D,NORM_FROBENIUS,&nrm);
2557: MatNorm(A,NORM_FROBENIUS,&gnorm);
2558: MatDestroy(&D);
2559: if (!gnorm) gnorm = 1; /* just in case */
2560: PetscViewerASCIIPrintf(viewer," ||J - Jfd||_F/||J||_F = %g, ||J - Jfd||_F = %g\n",(double)(nrm/gnorm),(double)nrm);
2562: if (complete_print) {
2563: PetscViewerASCIIPrintf(viewer," Hand-coded Jacobian ----------\n");
2564: MatView(A,mviewer);
2565: PetscViewerASCIIPrintf(viewer," Finite difference Jacobian ----------\n");
2566: MatView(B,mviewer);
2567: }
2569: if (threshold_print || complete_print) {
2570: PetscInt Istart, Iend, *ccols, bncols, cncols, j, row;
2571: PetscScalar *cvals;
2572: const PetscInt *bcols;
2573: const PetscScalar *bvals;
2575: MatCreate(PetscObjectComm((PetscObject)A),&C);
2576: MatSetType(C,mattype);
2577: MatSetSizes(C,m,n,M,N);
2578: MatSetBlockSizesFromMats(C,A,A);
2579: MatSetUp(C);
2580: MatSetOption(C,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);
2582: MatAYPX(B,-1.0,A,DIFFERENT_NONZERO_PATTERN);
2583: MatGetOwnershipRange(B,&Istart,&Iend);
2585: for (row = Istart; row < Iend; row++) {
2586: MatGetRow(B,row,&bncols,&bcols,&bvals);
2587: PetscMalloc2(bncols,&ccols,bncols,&cvals);
2588: for (j = 0, cncols = 0; j < bncols; j++) {
2589: if (PetscAbsScalar(bvals[j]) > threshold) {
2590: ccols[cncols] = bcols[j];
2591: cvals[cncols] = bvals[j];
2592: cncols += 1;
2593: }
2594: }
2595: if (cncols) {
2596: MatSetValues(C,1,&row,cncols,ccols,cvals,INSERT_VALUES);
2597: }
2598: MatRestoreRow(B,row,&bncols,&bcols,&bvals);
2599: PetscFree2(ccols,cvals);
2600: }
2601: MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
2602: MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
2603: PetscViewerASCIIPrintf(viewer," Hand-coded minus finite-difference Jacobian with tolerance %g ----------\n",(double)threshold);
2604: MatView(C,complete_print ? mviewer : viewer);
2605: MatDestroy(&C);
2606: }
2607: MatDestroy(&A);
2608: MatDestroy(&B);
2609: MatDestroy(&JT);
2610: if (Jsave) jacobian = Jsave;
2611: if (jacobian != snes->jacobian_pre) {
2612: jacobian = snes->jacobian_pre;
2613: PetscViewerASCIIPrintf(viewer," ---------- Testing Jacobian for preconditioner -------------\n");
2614: }
2615: else jacobian = NULL;
2616: }
2617: VecDestroy(&x);
2618: if (complete_print) {
2619: PetscViewerPopFormat(mviewer);
2620: }
2621: if (mviewer) { PetscViewerDestroy(&mviewer); }
2622: PetscViewerASCIISetTab(viewer,tabs);
2623: return(0);
2624: }
2626: /*@
2627: SNESComputeJacobian - Computes the Jacobian matrix that has been set with SNESSetJacobian().
2629: Collective on SNES
2631: Input Parameters:
2632: + snes - the SNES context
2633: - x - input vector
2635: Output Parameters:
2636: + A - Jacobian matrix
2637: - B - optional preconditioning matrix
2639: Options Database Keys:
2640: + -snes_lag_preconditioner <lag>
2641: . -snes_lag_jacobian <lag>
2642: . -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.
2643: . -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
2644: . -snes_compare_explicit - Compare the computed Jacobian to the finite difference Jacobian and output the differences
2645: . -snes_compare_explicit_draw - Compare the computed Jacobian to the finite difference Jacobian and draw the result
2646: . -snes_compare_explicit_contour - Compare the computed Jacobian to the finite difference Jacobian and draw a contour plot with the result
2647: . -snes_compare_operator - Make the comparison options above use the operator instead of the preconditioning matrix
2648: . -snes_compare_coloring - Compute the finite difference Jacobian using coloring and display norms of difference
2649: . -snes_compare_coloring_display - Compute the finite differece Jacobian using coloring and display verbose differences
2650: . -snes_compare_coloring_threshold - Display only those matrix entries that differ by more than a given threshold
2651: . -snes_compare_coloring_threshold_atol - Absolute tolerance for difference in matrix entries to be displayed by -snes_compare_coloring_threshold
2652: . -snes_compare_coloring_threshold_rtol - Relative tolerance for difference in matrix entries to be displayed by -snes_compare_coloring_threshold
2653: . -snes_compare_coloring_draw - Compute the finite differece Jacobian using coloring and draw differences
2654: - -snes_compare_coloring_draw_contour - Compute the finite differece Jacobian using coloring and show contours of matrices and differences
2657: Notes:
2658: Most users should not need to explicitly call this routine, as it
2659: is used internally within the nonlinear solvers.
2661: Developer Notes:
2662: 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
2663: for with the SNESType of test that has been removed.
2665: Level: developer
2667: .seealso: SNESSetJacobian(), KSPSetOperators(), MatStructure, SNESSetLagPreconditioner(), SNESSetLagJacobian()
2668: @*/
2669: PetscErrorCode SNESComputeJacobian(SNES snes,Vec X,Mat A,Mat B)
2670: {
2672: PetscBool flag;
2673: DM dm;
2674: DMSNES sdm;
2675: KSP ksp;
2681: VecValidValues(X,2,PETSC_TRUE);
2682: SNESGetDM(snes,&dm);
2683: DMGetDMSNES(dm,&sdm);
2685: if (!sdm->ops->computejacobian) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_USER,"Must call SNESSetJacobian(), DMSNESSetJacobian(), DMDASNESSetJacobianLocal(), etc");
2687: /* make sure that MatAssemblyBegin/End() is called on A matrix if it is matrix free */
2689: if (snes->lagjacobian == -2) {
2690: snes->lagjacobian = -1;
2692: PetscInfo(snes,"Recomputing Jacobian/preconditioner because lag is -2 (means compute Jacobian, but then never again) \n");
2693: } else if (snes->lagjacobian == -1) {
2694: PetscInfo(snes,"Reusing Jacobian/preconditioner because lag is -1\n");
2695: PetscObjectTypeCompare((PetscObject)A,MATMFFD,&flag);
2696: if (flag) {
2697: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
2698: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
2699: }
2700: return(0);
2701: } else if (snes->lagjacobian > 1 && (snes->iter + snes->jac_iter) % snes->lagjacobian) {
2702: PetscInfo2(snes,"Reusing Jacobian/preconditioner because lag is %D and SNES iteration is %D\n",snes->lagjacobian,snes->iter);
2703: PetscObjectTypeCompare((PetscObject)A,MATMFFD,&flag);
2704: if (flag) {
2705: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
2706: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
2707: }
2708: return(0);
2709: }
2710: if (snes->npc && snes->npcside== PC_LEFT) {
2711: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
2712: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
2713: return(0);
2714: }
2716: PetscLogEventBegin(SNES_JacobianEval,snes,X,A,B);
2717: VecLockReadPush(X);
2718: PetscStackPush("SNES user Jacobian function");
2719: (*sdm->ops->computejacobian)(snes,X,A,B,sdm->jacobianctx);
2720: PetscStackPop;
2721: VecLockReadPop(X);
2722: PetscLogEventEnd(SNES_JacobianEval,snes,X,A,B);
2724: /* attach latest linearization point to the preconditioning matrix */
2725: PetscObjectCompose((PetscObject)B,"__SNES_latest_X",(PetscObject)X);
2727: /* the next line ensures that snes->ksp exists */
2728: SNESGetKSP(snes,&ksp);
2729: if (snes->lagpreconditioner == -2) {
2730: PetscInfo(snes,"Rebuilding preconditioner exactly once since lag is -2\n");
2731: KSPSetReusePreconditioner(snes->ksp,PETSC_FALSE);
2732: snes->lagpreconditioner = -1;
2733: } else if (snes->lagpreconditioner == -1) {
2734: PetscInfo(snes,"Reusing preconditioner because lag is -1\n");
2735: KSPSetReusePreconditioner(snes->ksp,PETSC_TRUE);
2736: } else if (snes->lagpreconditioner > 1 && (snes->iter + snes->pre_iter) % snes->lagpreconditioner) {
2737: PetscInfo2(snes,"Reusing preconditioner because lag is %D and SNES iteration is %D\n",snes->lagpreconditioner,snes->iter);
2738: KSPSetReusePreconditioner(snes->ksp,PETSC_TRUE);
2739: } else {
2740: PetscInfo(snes,"Rebuilding preconditioner\n");
2741: KSPSetReusePreconditioner(snes->ksp,PETSC_FALSE);
2742: }
2744: SNESTestJacobian(snes);
2745: /* make sure user returned a correct Jacobian and preconditioner */
2748: {
2749: PetscBool flag = PETSC_FALSE,flag_draw = PETSC_FALSE,flag_contour = PETSC_FALSE,flag_operator = PETSC_FALSE;
2750: PetscOptionsGetViewer(PetscObjectComm((PetscObject)snes),((PetscObject) snes)->options,((PetscObject)snes)->prefix,"-snes_compare_explicit",NULL,NULL,&flag);
2751: PetscOptionsGetViewer(PetscObjectComm((PetscObject)snes),((PetscObject) snes)->options,((PetscObject)snes)->prefix,"-snes_compare_explicit_draw",NULL,NULL,&flag_draw);
2752: PetscOptionsGetViewer(PetscObjectComm((PetscObject)snes),((PetscObject) snes)->options,((PetscObject)snes)->prefix,"-snes_compare_explicit_draw_contour",NULL,NULL,&flag_contour);
2753: PetscOptionsGetViewer(PetscObjectComm((PetscObject)snes),((PetscObject) snes)->options,((PetscObject)snes)->prefix,"-snes_compare_operator",NULL,NULL,&flag_operator);
2754: if (flag || flag_draw || flag_contour) {
2755: Mat Bexp_mine = NULL,Bexp,FDexp;
2756: PetscViewer vdraw,vstdout;
2757: PetscBool flg;
2758: if (flag_operator) {
2759: MatComputeOperator(A,MATAIJ,&Bexp_mine);
2760: Bexp = Bexp_mine;
2761: } else {
2762: /* See if the preconditioning matrix can be viewed and added directly */
2763: PetscObjectBaseTypeCompareAny((PetscObject)B,&flg,MATSEQAIJ,MATMPIAIJ,MATSEQDENSE,MATMPIDENSE,MATSEQBAIJ,MATMPIBAIJ,MATSEQSBAIJ,MATMPIBAIJ,"");
2764: if (flg) Bexp = B;
2765: else {
2766: /* If the "preconditioning" matrix is itself MATSHELL or some other type without direct support */
2767: MatComputeOperator(B,MATAIJ,&Bexp_mine);
2768: Bexp = Bexp_mine;
2769: }
2770: }
2771: MatConvert(Bexp,MATSAME,MAT_INITIAL_MATRIX,&FDexp);
2772: SNESComputeJacobianDefault(snes,X,FDexp,FDexp,NULL);
2773: PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)snes),&vstdout);
2774: if (flag_draw || flag_contour) {
2775: PetscViewerDrawOpen(PetscObjectComm((PetscObject)snes),NULL,"Explicit Jacobians",PETSC_DECIDE,PETSC_DECIDE,300,300,&vdraw);
2776: if (flag_contour) {PetscViewerPushFormat(vdraw,PETSC_VIEWER_DRAW_CONTOUR);}
2777: } else vdraw = NULL;
2778: PetscViewerASCIIPrintf(vstdout,"Explicit %s\n",flag_operator ? "Jacobian" : "preconditioning Jacobian");
2779: if (flag) {MatView(Bexp,vstdout);}
2780: if (vdraw) {MatView(Bexp,vdraw);}
2781: PetscViewerASCIIPrintf(vstdout,"Finite difference Jacobian\n");
2782: if (flag) {MatView(FDexp,vstdout);}
2783: if (vdraw) {MatView(FDexp,vdraw);}
2784: MatAYPX(FDexp,-1.0,Bexp,SAME_NONZERO_PATTERN);
2785: PetscViewerASCIIPrintf(vstdout,"User-provided matrix minus finite difference Jacobian\n");
2786: if (flag) {MatView(FDexp,vstdout);}
2787: if (vdraw) { /* Always use contour for the difference */
2788: PetscViewerPushFormat(vdraw,PETSC_VIEWER_DRAW_CONTOUR);
2789: MatView(FDexp,vdraw);
2790: PetscViewerPopFormat(vdraw);
2791: }
2792: if (flag_contour) {PetscViewerPopFormat(vdraw);}
2793: PetscViewerDestroy(&vdraw);
2794: MatDestroy(&Bexp_mine);
2795: MatDestroy(&FDexp);
2796: }
2797: }
2798: {
2799: PetscBool flag = PETSC_FALSE,flag_display = PETSC_FALSE,flag_draw = PETSC_FALSE,flag_contour = PETSC_FALSE,flag_threshold = PETSC_FALSE;
2800: PetscReal threshold_atol = PETSC_SQRT_MACHINE_EPSILON,threshold_rtol = 10*PETSC_SQRT_MACHINE_EPSILON;
2801: PetscOptionsGetViewer(PetscObjectComm((PetscObject)snes),((PetscObject)snes)->options,((PetscObject)snes)->prefix,"-snes_compare_coloring",NULL,NULL,&flag);
2802: PetscOptionsGetViewer(PetscObjectComm((PetscObject)snes),((PetscObject)snes)->options,((PetscObject)snes)->prefix,"-snes_compare_coloring_display",NULL,NULL,&flag_display);
2803: PetscOptionsGetViewer(PetscObjectComm((PetscObject)snes),((PetscObject)snes)->options,((PetscObject)snes)->prefix,"-snes_compare_coloring_draw",NULL,NULL,&flag_draw);
2804: PetscOptionsGetViewer(PetscObjectComm((PetscObject)snes),((PetscObject)snes)->options,((PetscObject)snes)->prefix,"-snes_compare_coloring_draw_contour",NULL,NULL,&flag_contour);
2805: PetscOptionsGetViewer(PetscObjectComm((PetscObject)snes),((PetscObject)snes)->options,((PetscObject)snes)->prefix,"-snes_compare_coloring_threshold",NULL,NULL,&flag_threshold);
2806: if (flag_threshold) {
2807: PetscOptionsGetReal(((PetscObject)snes)->options,((PetscObject)snes)->prefix,"-snes_compare_coloring_threshold_rtol",&threshold_rtol,NULL);
2808: PetscOptionsGetReal(((PetscObject)snes)->options,((PetscObject)snes)->prefix,"-snes_compare_coloring_threshold_atol",&threshold_atol,NULL);
2809: }
2810: if (flag || flag_display || flag_draw || flag_contour || flag_threshold) {
2811: Mat Bfd;
2812: PetscViewer vdraw,vstdout;
2813: MatColoring coloring;
2814: ISColoring iscoloring;
2815: MatFDColoring matfdcoloring;
2816: PetscErrorCode (*func)(SNES,Vec,Vec,void*);
2817: void *funcctx;
2818: PetscReal norm1,norm2,normmax;
2820: MatDuplicate(B,MAT_DO_NOT_COPY_VALUES,&Bfd);
2821: MatColoringCreate(Bfd,&coloring);
2822: MatColoringSetType(coloring,MATCOLORINGSL);
2823: MatColoringSetFromOptions(coloring);
2824: MatColoringApply(coloring,&iscoloring);
2825: MatColoringDestroy(&coloring);
2826: MatFDColoringCreate(Bfd,iscoloring,&matfdcoloring);
2827: MatFDColoringSetFromOptions(matfdcoloring);
2828: MatFDColoringSetUp(Bfd,iscoloring,matfdcoloring);
2829: ISColoringDestroy(&iscoloring);
2831: /* This method of getting the function is currently unreliable since it doesn't work for DM local functions. */
2832: SNESGetFunction(snes,NULL,&func,&funcctx);
2833: MatFDColoringSetFunction(matfdcoloring,(PetscErrorCode (*)(void))func,funcctx);
2834: PetscObjectSetOptionsPrefix((PetscObject)matfdcoloring,((PetscObject)snes)->prefix);
2835: PetscObjectAppendOptionsPrefix((PetscObject)matfdcoloring,"coloring_");
2836: MatFDColoringSetFromOptions(matfdcoloring);
2837: MatFDColoringApply(Bfd,matfdcoloring,X,snes);
2838: MatFDColoringDestroy(&matfdcoloring);
2840: PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)snes),&vstdout);
2841: if (flag_draw || flag_contour) {
2842: PetscViewerDrawOpen(PetscObjectComm((PetscObject)snes),NULL,"Colored Jacobians",PETSC_DECIDE,PETSC_DECIDE,300,300,&vdraw);
2843: if (flag_contour) {PetscViewerPushFormat(vdraw,PETSC_VIEWER_DRAW_CONTOUR);}
2844: } else vdraw = NULL;
2845: PetscViewerASCIIPrintf(vstdout,"Explicit preconditioning Jacobian\n");
2846: if (flag_display) {MatView(B,vstdout);}
2847: if (vdraw) {MatView(B,vdraw);}
2848: PetscViewerASCIIPrintf(vstdout,"Colored Finite difference Jacobian\n");
2849: if (flag_display) {MatView(Bfd,vstdout);}
2850: if (vdraw) {MatView(Bfd,vdraw);}
2851: MatAYPX(Bfd,-1.0,B,SAME_NONZERO_PATTERN);
2852: MatNorm(Bfd,NORM_1,&norm1);
2853: MatNorm(Bfd,NORM_FROBENIUS,&norm2);
2854: MatNorm(Bfd,NORM_MAX,&normmax);
2855: PetscViewerASCIIPrintf(vstdout,"User-provided matrix minus finite difference Jacobian, norm1=%g normFrob=%g normmax=%g\n",(double)norm1,(double)norm2,(double)normmax);
2856: if (flag_display) {MatView(Bfd,vstdout);}
2857: if (vdraw) { /* Always use contour for the difference */
2858: PetscViewerPushFormat(vdraw,PETSC_VIEWER_DRAW_CONTOUR);
2859: MatView(Bfd,vdraw);
2860: PetscViewerPopFormat(vdraw);
2861: }
2862: if (flag_contour) {PetscViewerPopFormat(vdraw);}
2864: if (flag_threshold) {
2865: PetscInt bs,rstart,rend,i;
2866: MatGetBlockSize(B,&bs);
2867: MatGetOwnershipRange(B,&rstart,&rend);
2868: for (i=rstart; i<rend; i++) {
2869: const PetscScalar *ba,*ca;
2870: const PetscInt *bj,*cj;
2871: PetscInt bn,cn,j,maxentrycol = -1,maxdiffcol = -1,maxrdiffcol = -1;
2872: PetscReal maxentry = 0,maxdiff = 0,maxrdiff = 0;
2873: MatGetRow(B,i,&bn,&bj,&ba);
2874: MatGetRow(Bfd,i,&cn,&cj,&ca);
2875: if (bn != cn) SETERRQ(((PetscObject)A)->comm,PETSC_ERR_PLIB,"Unexpected different nonzero pattern in -snes_compare_coloring_threshold");
2876: for (j=0; j<bn; j++) {
2877: PetscReal rdiff = PetscAbsScalar(ca[j]) / (threshold_atol + threshold_rtol*PetscAbsScalar(ba[j]));
2878: if (PetscAbsScalar(ba[j]) > PetscAbs(maxentry)) {
2879: maxentrycol = bj[j];
2880: maxentry = PetscRealPart(ba[j]);
2881: }
2882: if (PetscAbsScalar(ca[j]) > PetscAbs(maxdiff)) {
2883: maxdiffcol = bj[j];
2884: maxdiff = PetscRealPart(ca[j]);
2885: }
2886: if (rdiff > maxrdiff) {
2887: maxrdiffcol = bj[j];
2888: maxrdiff = rdiff;
2889: }
2890: }
2891: if (maxrdiff > 1) {
2892: PetscViewerASCIIPrintf(vstdout,"row %D (maxentry=%g at %D, maxdiff=%g at %D, maxrdiff=%g at %D):",i,(double)maxentry,maxentrycol,(double)maxdiff,maxdiffcol,(double)maxrdiff,maxrdiffcol);
2893: for (j=0; j<bn; j++) {
2894: PetscReal rdiff;
2895: rdiff = PetscAbsScalar(ca[j]) / (threshold_atol + threshold_rtol*PetscAbsScalar(ba[j]));
2896: if (rdiff > 1) {
2897: PetscViewerASCIIPrintf(vstdout," (%D,%g:%g)",bj[j],(double)PetscRealPart(ba[j]),(double)PetscRealPart(ca[j]));
2898: }
2899: }
2900: PetscViewerASCIIPrintf(vstdout,"\n",i,maxentry,maxdiff,maxrdiff);
2901: }
2902: MatRestoreRow(B,i,&bn,&bj,&ba);
2903: MatRestoreRow(Bfd,i,&cn,&cj,&ca);
2904: }
2905: }
2906: PetscViewerDestroy(&vdraw);
2907: MatDestroy(&Bfd);
2908: }
2909: }
2910: return(0);
2911: }
2913: /*MC
2914: SNESJacobianFunction - Function used to convey the nonlinear Jacobian of the function to be solved by SNES
2916: Synopsis:
2917: #include "petscsnes.h"
2918: PetscErrorCode SNESJacobianFunction(SNES snes,Vec x,Mat Amat,Mat Pmat,void *ctx);
2920: Collective on snes
2922: Input Parameters:
2923: + x - input vector, the Jacobian is to be computed at this value
2924: - ctx - [optional] user-defined Jacobian context
2926: Output Parameters:
2927: + Amat - the matrix that defines the (approximate) Jacobian
2928: - Pmat - the matrix to be used in constructing the preconditioner, usually the same as Amat.
2930: Level: intermediate
2932: .seealso: SNESSetFunction(), SNESGetFunction(), SNESSetJacobian(), SNESGetJacobian()
2933: M*/
2935: /*@C
2936: SNESSetJacobian - Sets the function to compute Jacobian as well as the
2937: location to store the matrix.
2939: Logically Collective on SNES
2941: Input Parameters:
2942: + snes - the SNES context
2943: . Amat - the matrix that defines the (approximate) Jacobian
2944: . Pmat - the matrix to be used in constructing the preconditioner, usually the same as Amat.
2945: . J - Jacobian evaluation routine (if NULL then SNES retains any previously set value), see SNESJacobianFunction for details
2946: - ctx - [optional] user-defined context for private data for the
2947: Jacobian evaluation routine (may be NULL) (if NULL then SNES retains any previously set value)
2949: Notes:
2950: If the Amat matrix and Pmat matrix are different you must call MatAssemblyBegin/End() on
2951: each matrix.
2953: If you know the operator Amat has a null space you can use MatSetNullSpace() and MatSetTransposeNullSpace() to supply the null
2954: space to Amat and the KSP solvers will automatically use that null space as needed during the solution process.
2956: If using SNESComputeJacobianDefaultColor() to assemble a Jacobian, the ctx argument
2957: must be a MatFDColoring.
2959: Other defect-correction schemes can be used by computing a different matrix in place of the Jacobian. One common
2960: example is to use the "Picard linearization" which only differentiates through the highest order parts of each term.
2962: Level: beginner
2964: .seealso: KSPSetOperators(), SNESSetFunction(), MatMFFDComputeJacobian(), SNESComputeJacobianDefaultColor(), MatStructure, J,
2965: SNESSetPicard(), SNESJacobianFunction
2966: @*/
2967: PetscErrorCode SNESSetJacobian(SNES snes,Mat Amat,Mat Pmat,PetscErrorCode (*J)(SNES,Vec,Mat,Mat,void*),void *ctx)
2968: {
2970: DM dm;
2978: SNESGetDM(snes,&dm);
2979: DMSNESSetJacobian(dm,J,ctx);
2980: if (Amat) {
2981: PetscObjectReference((PetscObject)Amat);
2982: MatDestroy(&snes->jacobian);
2984: snes->jacobian = Amat;
2985: }
2986: if (Pmat) {
2987: PetscObjectReference((PetscObject)Pmat);
2988: MatDestroy(&snes->jacobian_pre);
2990: snes->jacobian_pre = Pmat;
2991: }
2992: return(0);
2993: }
2995: /*@C
2996: SNESGetJacobian - Returns the Jacobian matrix and optionally the user
2997: provided context for evaluating the Jacobian.
2999: Not Collective, but Mat object will be parallel if SNES object is
3001: Input Parameter:
3002: . snes - the nonlinear solver context
3004: Output Parameters:
3005: + Amat - location to stash (approximate) Jacobian matrix (or NULL)
3006: . Pmat - location to stash matrix used to compute the preconditioner (or NULL)
3007: . J - location to put Jacobian function (or NULL), see SNESJacobianFunction for details on its calling sequence
3008: - ctx - location to stash Jacobian ctx (or NULL)
3010: Level: advanced
3012: .seealso: SNESSetJacobian(), SNESComputeJacobian(), SNESJacobianFunction, SNESGetFunction()
3013: @*/
3014: PetscErrorCode SNESGetJacobian(SNES snes,Mat *Amat,Mat *Pmat,PetscErrorCode (**J)(SNES,Vec,Mat,Mat,void*),void **ctx)
3015: {
3017: DM dm;
3018: DMSNES sdm;
3022: if (Amat) *Amat = snes->jacobian;
3023: if (Pmat) *Pmat = snes->jacobian_pre;
3024: SNESGetDM(snes,&dm);
3025: DMGetDMSNES(dm,&sdm);
3026: if (J) *J = sdm->ops->computejacobian;
3027: if (ctx) *ctx = sdm->jacobianctx;
3028: return(0);
3029: }
3031: static PetscErrorCode SNESSetDefaultComputeJacobian(SNES snes)
3032: {
3034: DM dm;
3035: DMSNES sdm;
3038: SNESGetDM(snes,&dm);
3039: DMGetDMSNES(dm,&sdm);
3040: if (!sdm->ops->computejacobian && snes->jacobian_pre) {
3041: DM dm;
3042: PetscBool isdense,ismf;
3044: SNESGetDM(snes,&dm);
3045: PetscObjectTypeCompareAny((PetscObject)snes->jacobian_pre,&isdense,MATSEQDENSE,MATMPIDENSE,MATDENSE,NULL);
3046: PetscObjectTypeCompareAny((PetscObject)snes->jacobian_pre,&ismf,MATMFFD,MATSHELL,NULL);
3047: if (isdense) {
3048: DMSNESSetJacobian(dm,SNESComputeJacobianDefault,NULL);
3049: } else if (!ismf) {
3050: DMSNESSetJacobian(dm,SNESComputeJacobianDefaultColor,NULL);
3051: }
3052: }
3053: return(0);
3054: }
3056: /*@
3057: SNESSetUp - Sets up the internal data structures for the later use
3058: of a nonlinear solver.
3060: Collective on SNES
3062: Input Parameters:
3063: . snes - the SNES context
3065: Notes:
3066: For basic use of the SNES solvers the user need not explicitly call
3067: SNESSetUp(), since these actions will automatically occur during
3068: the call to SNESSolve(). However, if one wishes to control this
3069: phase separately, SNESSetUp() should be called after SNESCreate()
3070: and optional routines of the form SNESSetXXX(), but before SNESSolve().
3072: Level: advanced
3074: .seealso: SNESCreate(), SNESSolve(), SNESDestroy()
3075: @*/
3076: PetscErrorCode SNESSetUp(SNES snes)
3077: {
3079: DM dm;
3080: DMSNES sdm;
3081: SNESLineSearch linesearch, pclinesearch;
3082: void *lsprectx,*lspostctx;
3083: PetscErrorCode (*precheck)(SNESLineSearch,Vec,Vec,PetscBool*,void*);
3084: PetscErrorCode (*postcheck)(SNESLineSearch,Vec,Vec,Vec,PetscBool*,PetscBool*,void*);
3085: PetscErrorCode (*func)(SNES,Vec,Vec,void*);
3086: Vec f,fpc;
3087: void *funcctx;
3088: PetscErrorCode (*jac)(SNES,Vec,Mat,Mat,void*);
3089: void *jacctx,*appctx;
3090: Mat j,jpre;
3094: if (snes->setupcalled) return(0);
3095: PetscLogEventBegin(SNES_Setup,snes,0,0,0);
3097: if (!((PetscObject)snes)->type_name) {
3098: SNESSetType(snes,SNESNEWTONLS);
3099: }
3101: SNESGetFunction(snes,&snes->vec_func,NULL,NULL);
3103: SNESGetDM(snes,&dm);
3104: DMGetDMSNES(dm,&sdm);
3105: if (!sdm->ops->computefunction) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"Function never provided to SNES object");
3106: SNESSetDefaultComputeJacobian(snes);
3108: if (!snes->vec_func) {
3109: DMCreateGlobalVector(dm,&snes->vec_func);
3110: }
3112: if (!snes->ksp) {
3113: SNESGetKSP(snes, &snes->ksp);
3114: }
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) {
3164: SNESSetUpMatrixFree_Private(snes, snes->mf_operator, snes->mf_version);
3165: }
3166: if (snes->ops->usercompute && !snes->user) {
3167: (*snes->ops->usercompute)(snes,(void**)&snes->user);
3168: }
3170: snes->jac_iter = 0;
3171: snes->pre_iter = 0;
3173: if (snes->ops->setup) {
3174: (*snes->ops->setup)(snes);
3175: }
3177: SNESSetDefaultComputeJacobian(snes);
3179: if (snes->npc && (snes->npcside== PC_LEFT)) {
3180: if (snes->functype == SNES_FUNCTION_PRECONDITIONED) {
3181: if (snes->linesearch){
3182: SNESGetLineSearch(snes,&linesearch);
3183: SNESLineSearchSetFunction(linesearch,SNESComputeFunctionDefaultNPC);
3184: }
3185: }
3186: }
3187: PetscLogEventEnd(SNES_Setup,snes,0,0,0);
3188: snes->setupcalled = PETSC_TRUE;
3189: return(0);
3190: }
3192: /*@
3193: SNESReset - Resets a SNES context to the snessetupcalled = 0 state and removes any allocated Vecs and Mats
3195: Collective on SNES
3197: Input Parameter:
3198: . snes - iterative context obtained from SNESCreate()
3200: Level: intermediate
3202: Notes:
3203: Also calls the application context destroy routine set with SNESSetComputeApplicationContext()
3205: .seealso: SNESCreate(), SNESSetUp(), SNESSolve()
3206: @*/
3207: PetscErrorCode SNESReset(SNES snes)
3208: {
3213: if (snes->ops->userdestroy && snes->user) {
3214: (*snes->ops->userdestroy)((void**)&snes->user);
3215: snes->user = NULL;
3216: }
3217: if (snes->npc) {
3218: SNESReset(snes->npc);
3219: }
3221: if (snes->ops->reset) {
3222: (*snes->ops->reset)(snes);
3223: }
3224: if (snes->ksp) {
3225: KSPReset(snes->ksp);
3226: }
3228: if (snes->linesearch) {
3229: SNESLineSearchReset(snes->linesearch);
3230: }
3232: VecDestroy(&snes->vec_rhs);
3233: VecDestroy(&snes->vec_sol);
3234: VecDestroy(&snes->vec_sol_update);
3235: VecDestroy(&snes->vec_func);
3236: MatDestroy(&snes->jacobian);
3237: MatDestroy(&snes->jacobian_pre);
3238: VecDestroyVecs(snes->nwork,&snes->work);
3239: VecDestroyVecs(snes->nvwork,&snes->vwork);
3241: snes->alwayscomputesfinalresidual = PETSC_FALSE;
3243: snes->nwork = snes->nvwork = 0;
3244: snes->setupcalled = PETSC_FALSE;
3245: return(0);
3246: }
3248: /*@
3249: SNESDestroy - Destroys the nonlinear solver context that was created
3250: with SNESCreate().
3252: Collective on SNES
3254: Input Parameter:
3255: . snes - the SNES context
3257: Level: beginner
3259: .seealso: SNESCreate(), SNESSolve()
3260: @*/
3261: PetscErrorCode SNESDestroy(SNES *snes)
3262: {
3266: if (!*snes) return(0);
3268: if (--((PetscObject)(*snes))->refct > 0) {*snes = NULL; return(0);}
3270: SNESReset((*snes));
3271: SNESDestroy(&(*snes)->npc);
3273: /* if memory was published with SAWs then destroy it */
3274: PetscObjectSAWsViewOff((PetscObject)*snes);
3275: if ((*snes)->ops->destroy) {(*((*snes))->ops->destroy)((*snes));}
3277: if ((*snes)->dm) {DMCoarsenHookRemove((*snes)->dm,DMCoarsenHook_SNESVecSol,DMRestrictHook_SNESVecSol,*snes);}
3278: DMDestroy(&(*snes)->dm);
3279: KSPDestroy(&(*snes)->ksp);
3280: SNESLineSearchDestroy(&(*snes)->linesearch);
3282: PetscFree((*snes)->kspconvctx);
3283: if ((*snes)->ops->convergeddestroy) {
3284: (*(*snes)->ops->convergeddestroy)((*snes)->cnvP);
3285: }
3286: if ((*snes)->conv_hist_alloc) {
3287: PetscFree2((*snes)->conv_hist,(*snes)->conv_hist_its);
3288: }
3289: SNESMonitorCancel((*snes));
3290: PetscHeaderDestroy(snes);
3291: return(0);
3292: }
3294: /* ----------- Routines to set solver parameters ---------- */
3296: /*@
3297: SNESSetLagPreconditioner - Determines when the preconditioner is rebuilt in the nonlinear solve.
3299: Logically Collective on SNES
3301: Input Parameters:
3302: + snes - the SNES context
3303: - lag - -1 indicates NEVER rebuild, 1 means rebuild every time the Jacobian is computed within a single nonlinear solve, 2 means every second time
3304: the Jacobian is built etc. -2 indicates rebuild preconditioner at next chance but then never rebuild after that
3306: Options Database Keys:
3307: + -snes_lag_jacobian_persists <true,false> - sets the persistence
3308: . -snes_lag_jacobian <-2,1,2,...> - sets the lag
3309: . -snes_lag_preconditioner_persists <true,false> - sets the persistence
3310: - -snes_lag_preconditioner <-2,1,2,...> - sets the lag
3312: Notes:
3313: The default is 1
3314: The preconditioner is ALWAYS built in the first iteration of a nonlinear solve unless lag is -1 or SNESSetLagPreconditionerPersists() was called
3315: If -1 is used before the very first nonlinear solve the preconditioner is still built because there is no previous preconditioner to use
3317: Level: intermediate
3319: .seealso: SNESSetTrustRegionTolerance(), SNESGetLagPreconditioner(), SNESSetLagJacobian(), SNESGetLagJacobian(), SNESSetLagPreconditionerPersists(),
3320: SNESSetLagJacobianPersists()
3322: @*/
3323: PetscErrorCode SNESSetLagPreconditioner(SNES snes,PetscInt lag)
3324: {
3327: if (lag < -2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Lag must be -2, -1, 1 or greater");
3328: if (!lag) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Lag cannot be 0");
3330: snes->lagpreconditioner = lag;
3331: return(0);
3332: }
3334: /*@
3335: SNESSetGridSequence - sets the number of steps of grid sequencing that SNES does
3337: Logically Collective on SNES
3339: Input Parameters:
3340: + snes - the SNES context
3341: - steps - the number of refinements to do, defaults to 0
3343: Options Database Keys:
3344: . -snes_grid_sequence <steps>
3346: Level: intermediate
3348: Notes:
3349: Use SNESGetSolution() to extract the fine grid solution after grid sequencing.
3351: .seealso: SNESSetTrustRegionTolerance(), SNESGetLagPreconditioner(), SNESSetLagJacobian(), SNESGetLagJacobian(), SNESGetGridSequence()
3353: @*/
3354: PetscErrorCode SNESSetGridSequence(SNES snes,PetscInt steps)
3355: {
3359: snes->gridsequence = steps;
3360: return(0);
3361: }
3363: /*@
3364: SNESGetGridSequence - gets the number of steps of grid sequencing that SNES does
3366: Logically Collective on SNES
3368: Input Parameter:
3369: . snes - the SNES context
3371: Output Parameter:
3372: . steps - the number of refinements to do, defaults to 0
3374: Options Database Keys:
3375: . -snes_grid_sequence <steps>
3377: Level: intermediate
3379: Notes:
3380: Use SNESGetSolution() to extract the fine grid solution after grid sequencing.
3382: .seealso: SNESSetTrustRegionTolerance(), SNESGetLagPreconditioner(), SNESSetLagJacobian(), SNESGetLagJacobian(), SNESSetGridSequence()
3384: @*/
3385: PetscErrorCode SNESGetGridSequence(SNES snes,PetscInt *steps)
3386: {
3389: *steps = snes->gridsequence;
3390: return(0);
3391: }
3393: /*@
3394: SNESGetLagPreconditioner - Indicates how often the preconditioner is rebuilt
3396: Not Collective
3398: Input Parameter:
3399: . snes - the SNES context
3401: Output Parameter:
3402: . lag - -1 indicates NEVER rebuild, 1 means rebuild every time the Jacobian is computed within a single nonlinear solve, 2 means every second time
3403: the Jacobian is built etc. -2 indicates rebuild preconditioner at next chance but then never rebuild after that
3405: Options Database Keys:
3406: + -snes_lag_jacobian_persists <true,false> - sets the persistence
3407: . -snes_lag_jacobian <-2,1,2,...> - sets the lag
3408: . -snes_lag_preconditioner_persists <true,false> - sets the persistence
3409: - -snes_lag_preconditioner <-2,1,2,...> - sets the lag
3411: Notes:
3412: The default is 1
3413: The preconditioner is ALWAYS built in the first iteration of a nonlinear solve unless lag is -1
3415: Level: intermediate
3417: .seealso: SNESSetTrustRegionTolerance(), SNESSetLagPreconditioner(), SNESSetLagJacobianPersists(), SNESSetLagPreconditionerPersists()
3419: @*/
3420: PetscErrorCode SNESGetLagPreconditioner(SNES snes,PetscInt *lag)
3421: {
3424: *lag = snes->lagpreconditioner;
3425: return(0);
3426: }
3428: /*@
3429: SNESSetLagJacobian - Determines when the Jacobian is rebuilt in the nonlinear solve. See SNESSetLagPreconditioner() for determining how
3430: often the preconditioner is rebuilt.
3432: Logically Collective on SNES
3434: Input Parameters:
3435: + snes - the SNES context
3436: - lag - -1 indicates NEVER rebuild, 1 means rebuild every time the Jacobian is computed within a single nonlinear solve, 2 means every second time
3437: the Jacobian is built etc. -2 means rebuild at next chance but then never again
3439: Options Database Keys:
3440: + -snes_lag_jacobian_persists <true,false> - sets the persistence
3441: . -snes_lag_jacobian <-2,1,2,...> - sets the lag
3442: . -snes_lag_preconditioner_persists <true,false> - sets the persistence
3443: - -snes_lag_preconditioner <-2,1,2,...> - sets the lag.
3445: Notes:
3446: The default is 1
3447: The Jacobian is ALWAYS built in the first iteration of a nonlinear solve unless lag is -1
3448: 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
3449: at the next Newton step but never again (unless it is reset to another value)
3451: Level: intermediate
3453: .seealso: SNESSetTrustRegionTolerance(), SNESGetLagPreconditioner(), SNESSetLagPreconditioner(), SNESGetLagJacobianPersists(), SNESSetLagPreconditionerPersists()
3455: @*/
3456: PetscErrorCode SNESSetLagJacobian(SNES snes,PetscInt lag)
3457: {
3460: if (lag < -2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Lag must be -2, -1, 1 or greater");
3461: if (!lag) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Lag cannot be 0");
3463: snes->lagjacobian = lag;
3464: return(0);
3465: }
3467: /*@
3468: SNESGetLagJacobian - Indicates how often the Jacobian is rebuilt. See SNESGetLagPreconditioner() to determine when the preconditioner is rebuilt
3470: Not Collective
3472: Input Parameter:
3473: . snes - the SNES context
3475: Output Parameter:
3476: . lag - -1 indicates NEVER rebuild, 1 means rebuild every time the Jacobian is computed within a single nonlinear solve, 2 means every second time
3477: the Jacobian is built etc.
3479: Notes:
3480: The default is 1
3481: The jacobian is ALWAYS built in the first iteration of a nonlinear solve unless lag is -1 or SNESSetLagJacobianPersists() was called.
3483: Level: intermediate
3485: .seealso: SNESSetTrustRegionTolerance(), SNESSetLagJacobian(), SNESSetLagPreconditioner(), SNESGetLagPreconditioner(), SNESSetLagJacobianPersists(), SNESSetLagPreconditionerPersists()
3487: @*/
3488: PetscErrorCode SNESGetLagJacobian(SNES snes,PetscInt *lag)
3489: {
3492: *lag = snes->lagjacobian;
3493: return(0);
3494: }
3496: /*@
3497: SNESSetLagJacobianPersists - Set whether or not the Jacobian lagging persists through multiple solves
3499: Logically collective on SNES
3501: Input Parameter:
3502: + snes - the SNES context
3503: - flg - jacobian lagging persists if true
3505: Options Database Keys:
3506: + -snes_lag_jacobian_persists <true,false> - sets the persistence
3507: . -snes_lag_jacobian <-2,1,2,...> - sets the lag
3508: . -snes_lag_preconditioner_persists <true,false> - sets the persistence
3509: - -snes_lag_preconditioner <-2,1,2,...> - sets the lag
3512: Notes:
3513: This is useful both for nonlinear preconditioning, where it's appropriate to have the Jacobian be stale by
3514: several solves, and for implicit time-stepping, where Jacobian lagging in the inner nonlinear solve over several
3515: timesteps may present huge efficiency gains.
3517: Level: developer
3519: .seealso: SNESSetLagPreconditionerPersists(), SNESSetLagJacobian(), SNESGetLagJacobian(), SNESGetNPC(), SNESSetLagJacobianPersists()
3521: @*/
3522: PetscErrorCode SNESSetLagJacobianPersists(SNES snes,PetscBool flg)
3523: {
3527: snes->lagjac_persist = flg;
3528: return(0);
3529: }
3531: /*@
3532: SNESSetLagPreconditionerPersists - Set whether or not the preconditioner lagging persists through multiple solves
3534: Logically Collective on SNES
3536: Input Parameter:
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: This is useful both for nonlinear preconditioning, where it's appropriate to have the preconditioner be stale
3548: by several solves, and for implicit time-stepping, where preconditioner lagging in the inner nonlinear solve over
3549: several timesteps may present huge efficiency gains.
3551: Level: developer
3553: .seealso: SNESSetLagJacobianPersists(), SNESSetLagJacobian(), SNESGetLagJacobian(), SNESGetNPC(), SNESSetLagPreconditioner()
3555: @*/
3556: PetscErrorCode SNESSetLagPreconditionerPersists(SNES snes,PetscBool flg)
3557: {
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 Keys:
3575: . -snes_force_iteration <force> - Sets forcing an iteration
3577: Notes:
3578: This is used sometimes with TS to prevent TS from detecting a false steady state solution
3580: Level: intermediate
3582: .seealso: SNESSetTrustRegionTolerance(), SNESSetDivergenceTolerance()
3583: @*/
3584: PetscErrorCode SNESSetForceIteration(SNES snes,PetscBool force)
3585: {
3588: snes->forceiteration = force;
3589: return(0);
3590: }
3592: /*@
3593: SNESGetForceIteration - Whether or not to force SNESSolve() take at least one iteration regardless of the initial residual norm
3595: Logically Collective on SNES
3597: Input Parameters:
3598: . snes - the SNES context
3600: Output Parameter:
3601: . force - PETSC_TRUE requires at least one iteration.
3603: Level: intermediate
3605: .seealso: SNESSetForceIteration(), SNESSetTrustRegionTolerance(), SNESSetDivergenceTolerance()
3606: @*/
3607: PetscErrorCode SNESGetForceIteration(SNES snes,PetscBool *force)
3608: {
3611: *force = snes->forceiteration;
3612: return(0);
3613: }
3615: /*@
3616: SNESSetTolerances - Sets various parameters used in convergence tests.
3618: Logically Collective on SNES
3620: Input Parameters:
3621: + snes - the SNES context
3622: . abstol - absolute convergence tolerance
3623: . rtol - relative convergence tolerance
3624: . stol - convergence tolerance in terms of the norm of the change in the solution between steps, || delta x || < stol*|| x ||
3625: . maxit - maximum number of iterations
3626: - maxf - maximum number of function evaluations (-1 indicates no limit)
3628: Options Database Keys:
3629: + -snes_atol <abstol> - Sets abstol
3630: . -snes_rtol <rtol> - Sets rtol
3631: . -snes_stol <stol> - Sets stol
3632: . -snes_max_it <maxit> - Sets maxit
3633: - -snes_max_funcs <maxf> - Sets maxf
3635: Notes:
3636: The default maximum number of iterations is 50.
3637: The default maximum number of function evaluations is 1000.
3639: Level: intermediate
3641: .seealso: SNESSetTrustRegionTolerance(), SNESSetDivergenceTolerance(), SNESSetForceIteration()
3642: @*/
3643: PetscErrorCode SNESSetTolerances(SNES snes,PetscReal abstol,PetscReal rtol,PetscReal stol,PetscInt maxit,PetscInt maxf)
3644: {
3653: if (abstol != PETSC_DEFAULT) {
3654: if (abstol < 0.0) SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_OUTOFRANGE,"Absolute tolerance %g must be non-negative",(double)abstol);
3655: snes->abstol = abstol;
3656: }
3657: if (rtol != PETSC_DEFAULT) {
3658: if (rtol < 0.0 || 1.0 <= rtol) SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_OUTOFRANGE,"Relative tolerance %g must be non-negative and less than 1.0",(double)rtol);
3659: snes->rtol = rtol;
3660: }
3661: if (stol != PETSC_DEFAULT) {
3662: if (stol < 0.0) SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_OUTOFRANGE,"Step tolerance %g must be non-negative",(double)stol);
3663: snes->stol = stol;
3664: }
3665: if (maxit != PETSC_DEFAULT) {
3666: if (maxit < 0) SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_OUTOFRANGE,"Maximum number of iterations %D must be non-negative",maxit);
3667: snes->max_its = maxit;
3668: }
3669: if (maxf != PETSC_DEFAULT) {
3670: if (maxf < -1) SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_OUTOFRANGE,"Maximum number of function evaluations %D must be -1 or nonnegative",maxf);
3671: snes->max_funcs = maxf;
3672: }
3673: snes->tolerancesset = PETSC_TRUE;
3674: return(0);
3675: }
3677: /*@
3678: SNESSetDivergenceTolerance - Sets the divergence tolerance used for the SNES divergence test.
3680: Logically Collective on SNES
3682: Input Parameters:
3683: + snes - the SNES context
3684: - divtol - the divergence tolerance. Use -1 to deactivate the test.
3686: Options Database Keys:
3687: . -snes_divergence_tolerance <divtol> - Sets divtol
3689: Notes:
3690: The default divergence tolerance is 1e4.
3692: Level: intermediate
3694: .seealso: SNESSetTolerances(), SNESGetDivergenceTolerance
3695: @*/
3696: PetscErrorCode SNESSetDivergenceTolerance(SNES snes,PetscReal divtol)
3697: {
3702: if (divtol != PETSC_DEFAULT) {
3703: snes->divtol = divtol;
3704: }
3705: else {
3706: snes->divtol = 1.0e4;
3707: }
3708: return(0);
3709: }
3711: /*@
3712: SNESGetTolerances - Gets various parameters used in convergence tests.
3714: Not Collective
3716: Input Parameters:
3717: + snes - the SNES context
3718: . atol - absolute convergence tolerance
3719: . rtol - relative convergence tolerance
3720: . stol - convergence tolerance in terms of the norm
3721: of the change in the solution between steps
3722: . maxit - maximum number of iterations
3723: - maxf - maximum number of function evaluations
3725: Notes:
3726: The user can specify NULL for any parameter that is not needed.
3728: Level: intermediate
3730: .seealso: SNESSetTolerances()
3731: @*/
3732: PetscErrorCode SNESGetTolerances(SNES snes,PetscReal *atol,PetscReal *rtol,PetscReal *stol,PetscInt *maxit,PetscInt *maxf)
3733: {
3736: if (atol) *atol = snes->abstol;
3737: if (rtol) *rtol = snes->rtol;
3738: if (stol) *stol = snes->stol;
3739: if (maxit) *maxit = snes->max_its;
3740: if (maxf) *maxf = snes->max_funcs;
3741: return(0);
3742: }
3744: /*@
3745: SNESGetDivergenceTolerance - Gets divergence tolerance used in divergence test.
3747: Not Collective
3749: Input Parameters:
3750: + snes - the SNES context
3751: - divtol - divergence tolerance
3753: Level: intermediate
3755: .seealso: SNESSetDivergenceTolerance()
3756: @*/
3757: PetscErrorCode SNESGetDivergenceTolerance(SNES snes,PetscReal *divtol)
3758: {
3761: if (divtol) *divtol = snes->divtol;
3762: return(0);
3763: }
3765: /*@
3766: SNESSetTrustRegionTolerance - Sets the trust region parameter tolerance.
3768: Logically Collective on SNES
3770: Input Parameters:
3771: + snes - the SNES context
3772: - tol - tolerance
3774: Options Database Key:
3775: . -snes_trtol <tol> - Sets tol
3777: Level: intermediate
3779: .seealso: SNESSetTolerances()
3780: @*/
3781: PetscErrorCode SNESSetTrustRegionTolerance(SNES snes,PetscReal tol)
3782: {
3786: snes->deltatol = tol;
3787: return(0);
3788: }
3790: /*
3791: Duplicate the lg monitors for SNES from KSP; for some reason with
3792: dynamic libraries things don't work under Sun4 if we just use
3793: macros instead of functions
3794: */
3795: PetscErrorCode SNESMonitorLGResidualNorm(SNES snes,PetscInt it,PetscReal norm,void *ctx)
3796: {
3801: KSPMonitorLGResidualNorm((KSP)snes,it,norm,ctx);
3802: return(0);
3803: }
3805: PetscErrorCode SNESMonitorLGCreate(MPI_Comm comm,const char host[],const char label[],int x,int y,int m,int n,PetscDrawLG *lgctx)
3806: {
3810: KSPMonitorLGResidualNormCreate(comm,host,label,x,y,m,n,lgctx);
3811: return(0);
3812: }
3814: PETSC_INTERN PetscErrorCode SNESMonitorRange_Private(SNES,PetscInt,PetscReal*);
3816: PetscErrorCode SNESMonitorLGRange(SNES snes,PetscInt n,PetscReal rnorm,void *monctx)
3817: {
3818: PetscDrawLG lg;
3819: PetscErrorCode ierr;
3820: PetscReal x,y,per;
3821: PetscViewer v = (PetscViewer)monctx;
3822: static PetscReal prev; /* should be in the context */
3823: PetscDraw draw;
3827: PetscViewerDrawGetDrawLG(v,0,&lg);
3828: if (!n) {PetscDrawLGReset(lg);}
3829: PetscDrawLGGetDraw(lg,&draw);
3830: PetscDrawSetTitle(draw,"Residual norm");
3831: x = (PetscReal)n;
3832: if (rnorm > 0.0) y = PetscLog10Real(rnorm);
3833: else y = -15.0;
3834: PetscDrawLGAddPoint(lg,&x,&y);
3835: if (n < 20 || !(n % 5) || snes->reason) {
3836: PetscDrawLGDraw(lg);
3837: PetscDrawLGSave(lg);
3838: }
3840: PetscViewerDrawGetDrawLG(v,1,&lg);
3841: if (!n) {PetscDrawLGReset(lg);}
3842: PetscDrawLGGetDraw(lg,&draw);
3843: PetscDrawSetTitle(draw,"% elemts > .2*max elemt");
3844: SNESMonitorRange_Private(snes,n,&per);
3845: x = (PetscReal)n;
3846: y = 100.0*per;
3847: PetscDrawLGAddPoint(lg,&x,&y);
3848: if (n < 20 || !(n % 5) || snes->reason) {
3849: PetscDrawLGDraw(lg);
3850: PetscDrawLGSave(lg);
3851: }
3853: PetscViewerDrawGetDrawLG(v,2,&lg);
3854: if (!n) {prev = rnorm;PetscDrawLGReset(lg);}
3855: PetscDrawLGGetDraw(lg,&draw);
3856: PetscDrawSetTitle(draw,"(norm -oldnorm)/oldnorm");
3857: x = (PetscReal)n;
3858: y = (prev - rnorm)/prev;
3859: PetscDrawLGAddPoint(lg,&x,&y);
3860: if (n < 20 || !(n % 5) || snes->reason) {
3861: PetscDrawLGDraw(lg);
3862: PetscDrawLGSave(lg);
3863: }
3865: PetscViewerDrawGetDrawLG(v,3,&lg);
3866: if (!n) {PetscDrawLGReset(lg);}
3867: PetscDrawLGGetDraw(lg,&draw);
3868: PetscDrawSetTitle(draw,"(norm -oldnorm)/oldnorm*(% > .2 max)");
3869: x = (PetscReal)n;
3870: y = (prev - rnorm)/(prev*per);
3871: if (n > 2) { /*skip initial crazy value */
3872: PetscDrawLGAddPoint(lg,&x,&y);
3873: }
3874: if (n < 20 || !(n % 5) || snes->reason) {
3875: PetscDrawLGDraw(lg);
3876: PetscDrawLGSave(lg);
3877: }
3878: prev = rnorm;
3879: return(0);
3880: }
3882: /*@
3883: SNESMonitor - runs the user provided monitor routines, if they exist
3885: Collective on SNES
3887: Input Parameters:
3888: + snes - nonlinear solver context obtained from SNESCreate()
3889: . iter - iteration number
3890: - rnorm - relative norm of the residual
3892: Notes:
3893: This routine is called by the SNES implementations.
3894: It does not typically need to be called by the user.
3896: Level: developer
3898: .seealso: SNESMonitorSet()
3899: @*/
3900: PetscErrorCode SNESMonitor(SNES snes,PetscInt iter,PetscReal rnorm)
3901: {
3903: PetscInt i,n = snes->numbermonitors;
3906: VecLockReadPush(snes->vec_sol);
3907: for (i=0; i<n; i++) {
3908: (*snes->monitor[i])(snes,iter,rnorm,snes->monitorcontext[i]);
3909: }
3910: VecLockReadPop(snes->vec_sol);
3911: return(0);
3912: }
3914: /* ------------ Routines to set performance monitoring options ----------- */
3916: /*MC
3917: SNESMonitorFunction - functional form passed to SNESMonitorSet() to monitor convergence of nonlinear solver
3919: Synopsis:
3920: #include <petscsnes.h>
3921: $ PetscErrorCode SNESMonitorFunction(SNES snes,PetscInt its, PetscReal norm,void *mctx)
3923: Collective on snes
3925: Input Parameters:
3926: + snes - the SNES context
3927: . its - iteration number
3928: . norm - 2-norm function value (may be estimated)
3929: - mctx - [optional] monitoring context
3931: Level: advanced
3933: .seealso: SNESMonitorSet(), SNESMonitorGet()
3934: M*/
3936: /*@C
3937: SNESMonitorSet - Sets an ADDITIONAL function that is to be used at every
3938: iteration of the nonlinear solver to display the iteration's
3939: progress.
3941: Logically Collective on SNES
3943: Input Parameters:
3944: + snes - the SNES context
3945: . f - the monitor function, see SNESMonitorFunction for the calling sequence
3946: . mctx - [optional] user-defined context for private data for the
3947: monitor routine (use NULL if no context is desired)
3948: - monitordestroy - [optional] routine that frees monitor context
3949: (may be NULL)
3951: Options Database Keys:
3952: + -snes_monitor - sets SNESMonitorDefault()
3953: . -snes_monitor_lg_residualnorm - sets line graph monitor,
3954: uses SNESMonitorLGCreate()
3955: - -snes_monitor_cancel - cancels all monitors that have
3956: been hardwired into a code by
3957: calls to SNESMonitorSet(), but
3958: does not cancel those set via
3959: the options database.
3961: Notes:
3962: Several different monitoring routines may be set by calling
3963: SNESMonitorSet() multiple times; all will be called in the
3964: order in which they were set.
3966: Fortran Notes:
3967: Only a single monitor function can be set for each SNES object
3969: Level: intermediate
3971: .seealso: SNESMonitorDefault(), SNESMonitorCancel(), SNESMonitorFunction
3972: @*/
3973: PetscErrorCode SNESMonitorSet(SNES snes,PetscErrorCode (*f)(SNES,PetscInt,PetscReal,void*),void *mctx,PetscErrorCode (*monitordestroy)(void**))
3974: {
3975: PetscInt i;
3977: PetscBool identical;
3981: for (i=0; i<snes->numbermonitors;i++) {
3982: PetscMonitorCompare((PetscErrorCode (*)(void))f,mctx,monitordestroy,(PetscErrorCode (*)(void))snes->monitor[i],snes->monitorcontext[i],snes->monitordestroy[i],&identical);
3983: if (identical) return(0);
3984: }
3985: if (snes->numbermonitors >= MAXSNESMONITORS) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Too many monitors set");
3986: snes->monitor[snes->numbermonitors] = f;
3987: snes->monitordestroy[snes->numbermonitors] = monitordestroy;
3988: snes->monitorcontext[snes->numbermonitors++] = (void*)mctx;
3989: return(0);
3990: }
3992: /*@
3993: SNESMonitorCancel - Clears all the monitor functions for a SNES object.
3995: Logically Collective on SNES
3997: Input Parameters:
3998: . snes - the SNES context
4000: Options Database Key:
4001: . -snes_monitor_cancel - cancels all monitors that have been hardwired
4002: into a code by calls to SNESMonitorSet(), but does not cancel those
4003: set via the options database
4005: Notes:
4006: There is no way to clear one specific monitor from a SNES object.
4008: Level: intermediate
4010: .seealso: SNESMonitorDefault(), SNESMonitorSet()
4011: @*/
4012: PetscErrorCode SNESMonitorCancel(SNES snes)
4013: {
4015: PetscInt i;
4019: for (i=0; i<snes->numbermonitors; i++) {
4020: if (snes->monitordestroy[i]) {
4021: (*snes->monitordestroy[i])(&snes->monitorcontext[i]);
4022: }
4023: }
4024: snes->numbermonitors = 0;
4025: return(0);
4026: }
4028: /*MC
4029: SNESConvergenceTestFunction - functional form used for testing of convergence of nonlinear solver
4031: Synopsis:
4032: #include <petscsnes.h>
4033: $ PetscErrorCode SNESConvergenceTest(SNES snes,PetscInt it,PetscReal xnorm,PetscReal gnorm,PetscReal f,SNESConvergedReason *reason,void *cctx)
4035: Collective on snes
4037: Input Parameters:
4038: + snes - the SNES context
4039: . it - current iteration (0 is the first and is before any Newton step)
4040: . xnorm - 2-norm of current iterate
4041: . gnorm - 2-norm of current step
4042: . f - 2-norm of function
4043: - cctx - [optional] convergence context
4045: Output Parameter:
4046: . reason - reason for convergence/divergence, only needs to be set when convergence or divergence is detected
4048: Level: intermediate
4050: .seealso: SNESSetConvergenceTest(), SNESGetConvergenceTest()
4051: M*/
4053: /*@C
4054: SNESSetConvergenceTest - Sets the function that is to be used
4055: to test for convergence of the nonlinear iterative solution.
4057: Logically Collective on SNES
4059: Input Parameters:
4060: + snes - the SNES context
4061: . SNESConvergenceTestFunction - routine to test for convergence
4062: . cctx - [optional] context for private data for the convergence routine (may be NULL)
4063: - destroy - [optional] destructor for the context (may be NULL; PETSC_NULL_FUNCTION in Fortran)
4065: Level: advanced
4067: .seealso: SNESConvergedDefault(), SNESConvergedSkip(), SNESConvergenceTestFunction
4068: @*/
4069: PetscErrorCode SNESSetConvergenceTest(SNES snes,PetscErrorCode (*SNESConvergenceTestFunction)(SNES,PetscInt,PetscReal,PetscReal,PetscReal,SNESConvergedReason*,void*),void *cctx,PetscErrorCode (*destroy)(void*))
4070: {
4075: if (!SNESConvergenceTestFunction) SNESConvergenceTestFunction = SNESConvergedSkip;
4076: if (snes->ops->convergeddestroy) {
4077: (*snes->ops->convergeddestroy)(snes->cnvP);
4078: }
4079: snes->ops->converged = SNESConvergenceTestFunction;
4080: snes->ops->convergeddestroy = destroy;
4081: snes->cnvP = cctx;
4082: return(0);
4083: }
4085: /*@
4086: SNESGetConvergedReason - Gets the reason the SNES iteration was stopped.
4088: Not Collective
4090: Input Parameter:
4091: . snes - the SNES context
4093: Output Parameter:
4094: . reason - negative value indicates diverged, positive value converged, see SNESConvergedReason or the
4095: manual pages for the individual convergence tests for complete lists
4097: Options Database:
4098: . -snes_converged_reason - prints the reason to standard out
4100: Level: intermediate
4102: Notes:
4103: Should only be called after the call the SNESSolve() is complete, if it is called earlier it returns the value SNES__CONVERGED_ITERATING.
4105: .seealso: SNESSetConvergenceTest(), SNESSetConvergedReason(), SNESConvergedReason
4106: @*/
4107: PetscErrorCode SNESGetConvergedReason(SNES snes,SNESConvergedReason *reason)
4108: {
4112: *reason = snes->reason;
4113: return(0);
4114: }
4116: /*@
4117: SNESSetConvergedReason - Sets the reason the SNES iteration was stopped.
4119: Not Collective
4121: Input Parameters:
4122: + snes - the SNES context
4123: - reason - negative value indicates diverged, positive value converged, see SNESConvergedReason or the
4124: manual pages for the individual convergence tests for complete lists
4126: Level: intermediate
4128: .seealso: SNESGetConvergedReason(), SNESSetConvergenceTest(), SNESConvergedReason
4129: @*/
4130: PetscErrorCode SNESSetConvergedReason(SNES snes,SNESConvergedReason reason)
4131: {
4134: snes->reason = reason;
4135: return(0);
4136: }
4138: /*@
4139: SNESSetConvergenceHistory - Sets the array used to hold the convergence history.
4141: Logically Collective on SNES
4143: Input Parameters:
4144: + snes - iterative context obtained from SNESCreate()
4145: . a - array to hold history, this array will contain the function norms computed at each step
4146: . its - integer array holds the number of linear iterations for each solve.
4147: . na - size of a and its
4148: - reset - PETSC_TRUE indicates each new nonlinear solve resets the history counter to zero,
4149: else it continues storing new values for new nonlinear solves after the old ones
4151: Notes:
4152: If 'a' and 'its' are NULL then space is allocated for the history. If 'na' PETSC_DECIDE or PETSC_DEFAULT then a
4153: default array of length 10000 is allocated.
4155: This routine is useful, e.g., when running a code for purposes
4156: of accurate performance monitoring, when no I/O should be done
4157: during the section of code that is being timed.
4159: Level: intermediate
4161: .seealso: SNESGetConvergenceHistory()
4163: @*/
4164: PetscErrorCode SNESSetConvergenceHistory(SNES snes,PetscReal a[],PetscInt its[],PetscInt na,PetscBool reset)
4165: {
4172: if (!a) {
4173: if (na == PETSC_DECIDE || na == PETSC_DEFAULT) na = 1000;
4174: PetscCalloc2(na,&a,na,&its);
4175: snes->conv_hist_alloc = PETSC_TRUE;
4176: }
4177: snes->conv_hist = a;
4178: snes->conv_hist_its = its;
4179: snes->conv_hist_max = na;
4180: snes->conv_hist_len = 0;
4181: snes->conv_hist_reset = reset;
4182: return(0);
4183: }
4185: #if defined(PETSC_HAVE_MATLAB_ENGINE)
4186: #include <engine.h> /* MATLAB include file */
4187: #include <mex.h> /* MATLAB include file */
4189: PETSC_EXTERN mxArray *SNESGetConvergenceHistoryMatlab(SNES snes)
4190: {
4191: mxArray *mat;
4192: PetscInt i;
4193: PetscReal *ar;
4196: mat = mxCreateDoubleMatrix(snes->conv_hist_len,1,mxREAL);
4197: ar = (PetscReal*) mxGetData(mat);
4198: for (i=0; i<snes->conv_hist_len; i++) ar[i] = snes->conv_hist[i];
4199: PetscFunctionReturn(mat);
4200: }
4201: #endif
4203: /*@C
4204: SNESGetConvergenceHistory - Gets the array used to hold the convergence history.
4206: Not Collective
4208: Input Parameter:
4209: . snes - iterative context obtained from SNESCreate()
4211: Output Parameters:
4212: + a - array to hold history
4213: . its - integer array holds the number of linear iterations (or
4214: negative if not converged) for each solve.
4215: - na - size of a and its
4217: Notes:
4218: The calling sequence for this routine in Fortran is
4219: $ call SNESGetConvergenceHistory(SNES snes, integer na, integer ierr)
4221: This routine is useful, e.g., when running a code for purposes
4222: of accurate performance monitoring, when no I/O should be done
4223: during the section of code that is being timed.
4225: Level: intermediate
4227: .seealso: SNESSetConvergenceHistory()
4229: @*/
4230: PetscErrorCode SNESGetConvergenceHistory(SNES snes,PetscReal *a[],PetscInt *its[],PetscInt *na)
4231: {
4234: if (a) *a = snes->conv_hist;
4235: if (its) *its = snes->conv_hist_its;
4236: if (na) *na = snes->conv_hist_len;
4237: return(0);
4238: }
4240: /*@C
4241: SNESSetUpdate - Sets the general-purpose update function called
4242: at the beginning of every iteration of the nonlinear solve. Specifically
4243: it is called just before the Jacobian is "evaluated".
4245: Logically Collective on SNES
4247: Input Parameters:
4248: + snes - The nonlinear solver context
4249: - func - The function
4251: Calling sequence of func:
4252: $ func (SNES snes, PetscInt step);
4254: . step - The current step of the iteration
4256: Level: advanced
4258: Note: This is NOT what one uses to update the ghost points before a function evaluation, that should be done at the beginning of your FormFunction()
4259: This is not used by most users.
4261: .seealso SNESSetJacobian(), SNESSolve()
4262: @*/
4263: PetscErrorCode SNESSetUpdate(SNES snes, PetscErrorCode (*func)(SNES, PetscInt))
4264: {
4267: snes->ops->update = func;
4268: return(0);
4269: }
4271: /*
4272: SNESScaleStep_Private - Scales a step so that its length is less than the
4273: positive parameter delta.
4275: Input Parameters:
4276: + snes - the SNES context
4277: . y - approximate solution of linear system
4278: . fnorm - 2-norm of current function
4279: - delta - trust region size
4281: Output Parameters:
4282: + gpnorm - predicted function norm at the new point, assuming local
4283: linearization. The value is zero if the step lies within the trust
4284: region, and exceeds zero otherwise.
4285: - ynorm - 2-norm of the step
4287: Note:
4288: For non-trust region methods such as SNESNEWTONLS, the parameter delta
4289: is set to be the maximum allowable step size.
4291: */
4292: PetscErrorCode SNESScaleStep_Private(SNES snes,Vec y,PetscReal *fnorm,PetscReal *delta,PetscReal *gpnorm,PetscReal *ynorm)
4293: {
4294: PetscReal nrm;
4295: PetscScalar cnorm;
4303: VecNorm(y,NORM_2,&nrm);
4304: if (nrm > *delta) {
4305: nrm = *delta/nrm;
4306: *gpnorm = (1.0 - nrm)*(*fnorm);
4307: cnorm = nrm;
4308: VecScale(y,cnorm);
4309: *ynorm = *delta;
4310: } else {
4311: *gpnorm = 0.0;
4312: *ynorm = nrm;
4313: }
4314: return(0);
4315: }
4317: /*@C
4318: SNESConvergedReasonView - Displays the reason a SNES solve converged or diverged to a viewer
4320: Collective on SNES
4322: Parameter:
4323: + snes - iterative context obtained from SNESCreate()
4324: - viewer - the viewer to display the reason
4327: Options Database Keys:
4328: + -snes_converged_reason - print reason for converged or diverged, also prints number of iterations
4329: - -snes_converged_reason ::failed - only print reason and number of iterations when diverged
4331: Notes:
4332: To change the format of the output call PetscViewerPushFormat(viewer,format) before this call. Use PETSC_VIEWER_DEFAULT for the default,
4333: use PETSC_VIEWER_FAILED to only display a reason if it fails.
4335: Level: beginner
4337: .seealso: SNESCreate(), SNESSetUp(), SNESDestroy(), SNESSetTolerances(), SNESConvergedDefault(), SNESGetConvergedReason(), SNESConvergedReasonViewFromOptions(),
4338: PetscViewerPushFormat(), PetscViewerPopFormat()
4340: @*/
4341: PetscErrorCode SNESConvergedReasonView(SNES snes,PetscViewer viewer)
4342: {
4343: PetscViewerFormat format;
4344: PetscBool isAscii;
4345: PetscErrorCode ierr;
4348: if (!viewer) viewer = PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)snes));
4349: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isAscii);
4350: if (isAscii) {
4351: PetscViewerGetFormat(viewer, &format);
4352: PetscViewerASCIIAddTab(viewer,((PetscObject)snes)->tablevel);
4353: if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
4354: DM dm;
4355: Vec u;
4356: PetscDS prob;
4357: PetscInt Nf, f;
4358: PetscErrorCode (**exactSol)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar[], void *);
4359: void **exactCtx;
4360: PetscReal error;
4362: SNESGetDM(snes, &dm);
4363: SNESGetSolution(snes, &u);
4364: DMGetDS(dm, &prob);
4365: PetscDSGetNumFields(prob, &Nf);
4366: PetscMalloc2(Nf, &exactSol, Nf, &exactCtx);
4367: for (f = 0; f < Nf; ++f) {PetscDSGetExactSolution(prob, f, &exactSol[f], &exactCtx[f]);}
4368: DMComputeL2Diff(dm, 0.0, exactSol, exactCtx, u, &error);
4369: PetscFree2(exactSol, exactCtx);
4370: if (error < 1.0e-11) {PetscViewerASCIIPrintf(viewer, "L_2 Error: < 1.0e-11\n");}
4371: else {PetscViewerASCIIPrintf(viewer, "L_2 Error: %g\n", error);}
4372: }
4373: if (snes->reason > 0 && format != PETSC_VIEWER_FAILED) {
4374: if (((PetscObject) snes)->prefix) {
4375: PetscViewerASCIIPrintf(viewer,"Nonlinear %s solve converged due to %s iterations %D\n",((PetscObject) snes)->prefix,SNESConvergedReasons[snes->reason],snes->iter);
4376: } else {
4377: PetscViewerASCIIPrintf(viewer,"Nonlinear solve converged due to %s iterations %D\n",SNESConvergedReasons[snes->reason],snes->iter);
4378: }
4379: } else if (snes->reason <= 0) {
4380: if (((PetscObject) snes)->prefix) {
4381: PetscViewerASCIIPrintf(viewer,"Nonlinear %s solve did not converge due to %s iterations %D\n",((PetscObject) snes)->prefix,SNESConvergedReasons[snes->reason],snes->iter);
4382: } else {
4383: PetscViewerASCIIPrintf(viewer,"Nonlinear solve did not converge due to %s iterations %D\n",SNESConvergedReasons[snes->reason],snes->iter);
4384: }
4385: }
4386: PetscViewerASCIISubtractTab(viewer,((PetscObject)snes)->tablevel);
4387: }
4388: return(0);
4389: }
4391: /*@
4392: SNESConvergedReasonViewFromOptions - Processes command line options to determine if/how a SNESReason is to be viewed.
4394: Collective on SNES
4396: Input Parameters:
4397: . snes - the SNES object
4399: Level: intermediate
4401: .seealso: SNESCreate(), SNESSetUp(), SNESDestroy(), SNESSetTolerances(), SNESConvergedDefault(), SNESGetConvergedReason(), SNESConvergedReasonView()
4403: @*/
4404: PetscErrorCode SNESConvergedReasonViewFromOptions(SNES snes)
4405: {
4406: PetscErrorCode ierr;
4407: PetscViewer viewer;
4408: PetscBool flg;
4409: static PetscBool incall = PETSC_FALSE;
4410: PetscViewerFormat format;
4413: if (incall) return(0);
4414: incall = PETSC_TRUE;
4415: PetscOptionsGetViewer(PetscObjectComm((PetscObject)snes),((PetscObject)snes)->options,((PetscObject)snes)->prefix,"-snes_converged_reason",&viewer,&format,&flg);
4416: if (flg) {
4417: PetscViewerPushFormat(viewer,format);
4418: SNESConvergedReasonView(snes,viewer);
4419: PetscViewerPopFormat(viewer);
4420: PetscViewerDestroy(&viewer);
4421: }
4422: incall = PETSC_FALSE;
4423: return(0);
4424: }
4426: /*@
4427: SNESSolve - Solves a nonlinear system F(x) = b.
4428: Call SNESSolve() after calling SNESCreate() and optional routines of the form SNESSetXXX().
4430: Collective on SNES
4432: Input Parameters:
4433: + snes - the SNES context
4434: . b - the constant part of the equation F(x) = b, or NULL to use zero.
4435: - x - the solution vector.
4437: Notes:
4438: The user should initialize the vector,x, with the initial guess
4439: for the nonlinear solve prior to calling SNESSolve. In particular,
4440: to employ an initial guess of zero, the user should explicitly set
4441: this vector to zero by calling VecSet().
4443: Level: beginner
4445: .seealso: SNESCreate(), SNESDestroy(), SNESSetFunction(), SNESSetJacobian(), SNESSetGridSequence(), SNESGetSolution()
4446: @*/
4447: PetscErrorCode SNESSolve(SNES snes,Vec b,Vec x)
4448: {
4449: PetscErrorCode ierr;
4450: PetscBool flg;
4451: PetscInt grid;
4452: Vec xcreated = NULL;
4453: DM dm;
4462: /* High level operations using the nonlinear solver */
4463: {
4464: PetscViewer viewer;
4465: PetscViewerFormat format;
4466: PetscInt num;
4467: PetscBool flg;
4468: static PetscBool incall = PETSC_FALSE;
4470: if (!incall) {
4471: /* Estimate the convergence rate of the discretization */
4472: PetscOptionsGetViewer(PetscObjectComm((PetscObject) snes),((PetscObject)snes)->options, ((PetscObject) snes)->prefix, "-snes_convergence_estimate", &viewer, &format, &flg);
4473: if (flg) {
4474: PetscConvEst conv;
4475: DM dm;
4476: PetscReal *alpha; /* Convergence rate of the solution error for each field in the L_2 norm */
4477: PetscInt Nf;
4479: incall = PETSC_TRUE;
4480: SNESGetDM(snes, &dm);
4481: DMGetNumFields(dm, &Nf);
4482: PetscCalloc1(Nf, &alpha);
4483: PetscConvEstCreate(PetscObjectComm((PetscObject) snes), &conv);
4484: PetscConvEstSetSolver(conv, (PetscObject) snes);
4485: PetscConvEstSetFromOptions(conv);
4486: PetscConvEstSetUp(conv);
4487: PetscConvEstGetConvRate(conv, alpha);
4488: PetscViewerPushFormat(viewer, format);
4489: PetscConvEstRateView(conv, alpha, viewer);
4490: PetscViewerPopFormat(viewer);
4491: PetscViewerDestroy(&viewer);
4492: PetscConvEstDestroy(&conv);
4493: PetscFree(alpha);
4494: incall = PETSC_FALSE;
4495: }
4496: /* Adaptively refine the initial grid */
4497: num = 1;
4498: PetscOptionsGetInt(NULL, ((PetscObject) snes)->prefix, "-snes_adapt_initial", &num, &flg);
4499: if (flg) {
4500: DMAdaptor adaptor;
4502: incall = PETSC_TRUE;
4503: DMAdaptorCreate(PetscObjectComm((PetscObject)snes), &adaptor);
4504: DMAdaptorSetSolver(adaptor, snes);
4505: DMAdaptorSetSequenceLength(adaptor, num);
4506: DMAdaptorSetFromOptions(adaptor);
4507: DMAdaptorSetUp(adaptor);
4508: DMAdaptorAdapt(adaptor, x, DM_ADAPTATION_INITIAL, &dm, &x);
4509: DMAdaptorDestroy(&adaptor);
4510: incall = PETSC_FALSE;
4511: }
4512: /* Use grid sequencing to adapt */
4513: num = 0;
4514: PetscOptionsGetInt(NULL, ((PetscObject) snes)->prefix, "-snes_adapt_sequence", &num, NULL);
4515: if (num) {
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_SEQUENTIAL, &dm, &x);
4525: DMAdaptorDestroy(&adaptor);
4526: incall = PETSC_FALSE;
4527: }
4528: }
4529: }
4530: if (!x) {
4531: SNESGetDM(snes,&dm);
4532: DMCreateGlobalVector(dm,&xcreated);
4533: x = xcreated;
4534: }
4535: SNESViewFromOptions(snes,NULL,"-snes_view_pre");
4537: for (grid=0; grid<snes->gridsequence; grid++) {PetscViewerASCIIPushTab(PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)snes)));}
4538: for (grid=0; grid<snes->gridsequence+1; grid++) {
4540: /* set solution vector */
4541: if (!grid) {PetscObjectReference((PetscObject)x);}
4542: VecDestroy(&snes->vec_sol);
4543: snes->vec_sol = x;
4544: SNESGetDM(snes,&dm);
4546: /* set affine vector if provided */
4547: if (b) { PetscObjectReference((PetscObject)b); }
4548: VecDestroy(&snes->vec_rhs);
4549: snes->vec_rhs = b;
4551: if (snes->vec_rhs && (snes->vec_func == snes->vec_rhs)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_IDN,"Right hand side vector cannot be function vector");
4552: if (snes->vec_func == snes->vec_sol) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_IDN,"Solution vector cannot be function vector");
4553: if (snes->vec_rhs == snes->vec_sol) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_IDN,"Solution vector cannot be right hand side vector");
4554: if (!snes->vec_sol_update /* && snes->vec_sol */) {
4555: VecDuplicate(snes->vec_sol,&snes->vec_sol_update);
4556: PetscLogObjectParent((PetscObject)snes,(PetscObject)snes->vec_sol_update);
4557: }
4558: DMShellSetGlobalVector(dm,snes->vec_sol);
4559: SNESSetUp(snes);
4561: if (!grid) {
4562: if (snes->ops->computeinitialguess) {
4563: (*snes->ops->computeinitialguess)(snes,snes->vec_sol,snes->initialguessP);
4564: }
4565: }
4567: if (snes->conv_hist_reset) snes->conv_hist_len = 0;
4568: if (snes->counters_reset) {snes->nfuncs = 0; snes->linear_its = 0; snes->numFailures = 0;}
4570: PetscLogEventBegin(SNES_Solve,snes,0,0,0);
4571: (*snes->ops->solve)(snes);
4572: PetscLogEventEnd(SNES_Solve,snes,0,0,0);
4573: if (!snes->reason) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Internal error, solver returned without setting converged reason");
4574: snes->domainerror = PETSC_FALSE; /* clear the flag if it has been set */
4576: if (snes->lagjac_persist) snes->jac_iter += snes->iter;
4577: if (snes->lagpre_persist) snes->pre_iter += snes->iter;
4579: PetscOptionsGetViewer(PetscObjectComm((PetscObject)snes),((PetscObject)snes)->options,((PetscObject)snes)->prefix,"-snes_test_local_min",NULL,NULL,&flg);
4580: if (flg && !PetscPreLoadingOn) { SNESTestLocalMin(snes); }
4581: SNESConvergedReasonViewFromOptions(snes);
4583: if (snes->errorifnotconverged && snes->reason < 0) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_NOT_CONVERGED,"SNESSolve has not converged");
4584: if (snes->reason < 0) break;
4585: if (grid < snes->gridsequence) {
4586: DM fine;
4587: Vec xnew;
4588: Mat interp;
4590: DMRefine(snes->dm,PetscObjectComm((PetscObject)snes),&fine);
4591: if (!fine) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_INCOMP,"DMRefine() did not perform any refinement, cannot continue grid sequencing");
4592: DMCreateInterpolation(snes->dm,fine,&interp,NULL);
4593: DMCreateGlobalVector(fine,&xnew);
4594: MatInterpolate(interp,x,xnew);
4595: DMInterpolate(snes->dm,interp,fine);
4596: MatDestroy(&interp);
4597: x = xnew;
4599: SNESReset(snes);
4600: SNESSetDM(snes,fine);
4601: SNESResetFromOptions(snes);
4602: DMDestroy(&fine);
4603: PetscViewerASCIIPopTab(PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)snes)));
4604: }
4605: }
4606: SNESViewFromOptions(snes,NULL,"-snes_view");
4607: VecViewFromOptions(snes->vec_sol,(PetscObject)snes,"-snes_view_solution");
4608: DMMonitor(snes->dm);
4610: VecDestroy(&xcreated);
4611: PetscObjectSAWsBlock((PetscObject)snes);
4612: return(0);
4613: }
4615: /* --------- Internal routines for SNES Package --------- */
4617: /*@C
4618: SNESSetType - Sets the method for the nonlinear solver.
4620: Collective on SNES
4622: Input Parameters:
4623: + snes - the SNES context
4624: - type - a known method
4626: Options Database Key:
4627: . -snes_type <type> - Sets the method; use -help for a list
4628: of available methods (for instance, newtonls or newtontr)
4630: Notes:
4631: See "petsc/include/petscsnes.h" for available methods (for instance)
4632: + SNESNEWTONLS - Newton's method with line search
4633: (systems of nonlinear equations)
4634: - SNESNEWTONTR - Newton's method with trust region
4635: (systems of nonlinear equations)
4637: Normally, it is best to use the SNESSetFromOptions() command and then
4638: set the SNES solver type from the options database rather than by using
4639: this routine. Using the options database provides the user with
4640: maximum flexibility in evaluating the many nonlinear solvers.
4641: The SNESSetType() routine is provided for those situations where it
4642: is necessary to set the nonlinear solver independently of the command
4643: line or options database. This might be the case, for example, when
4644: the choice of solver changes during the execution of the program,
4645: and the user's application is taking responsibility for choosing the
4646: appropriate method.
4648: Developer Notes:
4649: SNESRegister() adds a constructor for a new SNESType to SNESList, SNESSetType() locates
4650: the constructor in that list and calls it to create the spexific object.
4652: Level: intermediate
4654: .seealso: SNESType, SNESCreate(), SNESDestroy(), SNESGetType(), SNESSetFromOptions()
4656: @*/
4657: PetscErrorCode SNESSetType(SNES snes,SNESType type)
4658: {
4659: PetscErrorCode ierr,(*r)(SNES);
4660: PetscBool match;
4666: PetscObjectTypeCompare((PetscObject)snes,type,&match);
4667: if (match) return(0);
4669: PetscFunctionListFind(SNESList,type,&r);
4670: if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unable to find requested SNES type %s",type);
4671: /* Destroy the previous private SNES context */
4672: if (snes->ops->destroy) {
4673: (*(snes)->ops->destroy)(snes);
4674: snes->ops->destroy = NULL;
4675: }
4676: /* Reinitialize function pointers in SNESOps structure */
4677: snes->ops->setup = NULL;
4678: snes->ops->solve = NULL;
4679: snes->ops->view = NULL;
4680: snes->ops->setfromoptions = NULL;
4681: snes->ops->destroy = NULL;
4683: /* It may happen the user has customized the line search before calling SNESSetType */
4684: if (((PetscObject)snes)->type_name) {
4685: SNESLineSearchDestroy(&snes->linesearch);
4686: }
4688: /* Call the SNESCreate_XXX routine for this particular Nonlinear solver */
4689: snes->setupcalled = PETSC_FALSE;
4691: PetscObjectChangeTypeName((PetscObject)snes,type);
4692: (*r)(snes);
4693: return(0);
4694: }
4696: /*@C
4697: SNESGetType - Gets the SNES method type and name (as a string).
4699: Not Collective
4701: Input Parameter:
4702: . snes - nonlinear solver context
4704: Output Parameter:
4705: . type - SNES method (a character string)
4707: Level: intermediate
4709: @*/
4710: PetscErrorCode SNESGetType(SNES snes,SNESType *type)
4711: {
4715: *type = ((PetscObject)snes)->type_name;
4716: return(0);
4717: }
4719: /*@
4720: SNESSetSolution - Sets the solution vector for use by the SNES routines.
4722: Logically Collective on SNES
4724: Input Parameters:
4725: + snes - the SNES context obtained from SNESCreate()
4726: - u - the solution vector
4728: Level: beginner
4730: @*/
4731: PetscErrorCode SNESSetSolution(SNES snes, Vec u)
4732: {
4733: DM dm;
4739: PetscObjectReference((PetscObject) u);
4740: VecDestroy(&snes->vec_sol);
4742: snes->vec_sol = u;
4744: SNESGetDM(snes, &dm);
4745: DMShellSetGlobalVector(dm, u);
4746: return(0);
4747: }
4749: /*@
4750: SNESGetSolution - Returns the vector where the approximate solution is
4751: stored. This is the fine grid solution when using SNESSetGridSequence().
4753: Not Collective, but Vec is parallel if SNES is parallel
4755: Input Parameter:
4756: . snes - the SNES context
4758: Output Parameter:
4759: . x - the solution
4761: Level: intermediate
4763: .seealso: SNESGetSolutionUpdate(), SNESGetFunction()
4764: @*/
4765: PetscErrorCode SNESGetSolution(SNES snes,Vec *x)
4766: {
4770: *x = snes->vec_sol;
4771: return(0);
4772: }
4774: /*@
4775: SNESGetSolutionUpdate - Returns the vector where the solution update is
4776: stored.
4778: Not Collective, but Vec is parallel if SNES is parallel
4780: Input Parameter:
4781: . snes - the SNES context
4783: Output Parameter:
4784: . x - the solution update
4786: Level: advanced
4788: .seealso: SNESGetSolution(), SNESGetFunction()
4789: @*/
4790: PetscErrorCode SNESGetSolutionUpdate(SNES snes,Vec *x)
4791: {
4795: *x = snes->vec_sol_update;
4796: return(0);
4797: }
4799: /*@C
4800: SNESGetFunction - Returns the vector where the function is stored.
4802: Not Collective, but Vec is parallel if SNES is parallel. Collective if Vec is requested, but has not been created yet.
4804: Input Parameter:
4805: . snes - the SNES context
4807: Output Parameter:
4808: + r - the vector that is used to store residuals (or NULL if you don't want it)
4809: . f - the function (or NULL if you don't want it); see SNESFunction for calling sequence details
4810: - ctx - the function context (or NULL if you don't want it)
4812: Level: advanced
4814: Notes: The vector r DOES NOT, in general contain the current value of the SNES nonlinear function
4816: .seealso: SNESSetFunction(), SNESGetSolution(), SNESFunction
4817: @*/
4818: PetscErrorCode SNESGetFunction(SNES snes,Vec *r,PetscErrorCode (**f)(SNES,Vec,Vec,void*),void **ctx)
4819: {
4821: DM dm;
4825: if (r) {
4826: if (!snes->vec_func) {
4827: if (snes->vec_rhs) {
4828: VecDuplicate(snes->vec_rhs,&snes->vec_func);
4829: } else if (snes->vec_sol) {
4830: VecDuplicate(snes->vec_sol,&snes->vec_func);
4831: } else if (snes->dm) {
4832: DMCreateGlobalVector(snes->dm,&snes->vec_func);
4833: }
4834: }
4835: *r = snes->vec_func;
4836: }
4837: SNESGetDM(snes,&dm);
4838: DMSNESGetFunction(dm,f,ctx);
4839: return(0);
4840: }
4842: /*@C
4843: SNESGetNGS - Returns the NGS function and context.
4845: Input Parameter:
4846: . snes - the SNES context
4848: Output Parameter:
4849: + f - the function (or NULL) see SNESNGSFunction for details
4850: - ctx - the function context (or NULL)
4852: Level: advanced
4854: .seealso: SNESSetNGS(), SNESGetFunction()
4855: @*/
4857: PetscErrorCode SNESGetNGS (SNES snes, PetscErrorCode (**f)(SNES, Vec, Vec, void*), void ** ctx)
4858: {
4860: DM dm;
4864: SNESGetDM(snes,&dm);
4865: DMSNESGetNGS(dm,f,ctx);
4866: return(0);
4867: }
4869: /*@C
4870: SNESSetOptionsPrefix - Sets the prefix used for searching for all
4871: SNES options in the database.
4873: Logically Collective on SNES
4875: Input Parameter:
4876: + snes - the SNES context
4877: - prefix - the prefix to prepend to all option names
4879: Notes:
4880: A hyphen (-) must NOT be given at the beginning of the prefix name.
4881: The first character of all runtime options is AUTOMATICALLY the hyphen.
4883: Level: advanced
4885: .seealso: SNESSetFromOptions()
4886: @*/
4887: PetscErrorCode SNESSetOptionsPrefix(SNES snes,const char prefix[])
4888: {
4893: PetscObjectSetOptionsPrefix((PetscObject)snes,prefix);
4894: if (!snes->ksp) {SNESGetKSP(snes,&snes->ksp);}
4895: if (snes->linesearch) {
4896: SNESGetLineSearch(snes,&snes->linesearch);
4897: PetscObjectSetOptionsPrefix((PetscObject)snes->linesearch,prefix);
4898: }
4899: KSPSetOptionsPrefix(snes->ksp,prefix);
4900: return(0);
4901: }
4903: /*@C
4904: SNESAppendOptionsPrefix - Appends to the prefix used for searching for all
4905: SNES options in the database.
4907: Logically Collective on SNES
4909: Input Parameters:
4910: + snes - the SNES context
4911: - prefix - the prefix to prepend to all option names
4913: Notes:
4914: A hyphen (-) must NOT be given at the beginning of the prefix name.
4915: The first character of all runtime options is AUTOMATICALLY the hyphen.
4917: Level: advanced
4919: .seealso: SNESGetOptionsPrefix()
4920: @*/
4921: PetscErrorCode SNESAppendOptionsPrefix(SNES snes,const char prefix[])
4922: {
4927: PetscObjectAppendOptionsPrefix((PetscObject)snes,prefix);
4928: if (!snes->ksp) {SNESGetKSP(snes,&snes->ksp);}
4929: if (snes->linesearch) {
4930: SNESGetLineSearch(snes,&snes->linesearch);
4931: PetscObjectAppendOptionsPrefix((PetscObject)snes->linesearch,prefix);
4932: }
4933: KSPAppendOptionsPrefix(snes->ksp,prefix);
4934: return(0);
4935: }
4937: /*@C
4938: SNESGetOptionsPrefix - Sets the prefix used for searching for all
4939: SNES options in the database.
4941: Not Collective
4943: Input Parameter:
4944: . snes - the SNES context
4946: Output Parameter:
4947: . prefix - pointer to the prefix string used
4949: Notes:
4950: On the fortran side, the user should pass in a string 'prefix' of
4951: sufficient length to hold the prefix.
4953: Level: advanced
4955: .seealso: SNESAppendOptionsPrefix()
4956: @*/
4957: PetscErrorCode SNESGetOptionsPrefix(SNES snes,const char *prefix[])
4958: {
4963: PetscObjectGetOptionsPrefix((PetscObject)snes,prefix);
4964: return(0);
4965: }
4968: /*@C
4969: SNESRegister - Adds a method to the nonlinear solver package.
4971: Not collective
4973: Input Parameters:
4974: + name_solver - name of a new user-defined solver
4975: - routine_create - routine to create method context
4977: Notes:
4978: SNESRegister() may be called multiple times to add several user-defined solvers.
4980: Sample usage:
4981: .vb
4982: SNESRegister("my_solver",MySolverCreate);
4983: .ve
4985: Then, your solver can be chosen with the procedural interface via
4986: $ SNESSetType(snes,"my_solver")
4987: or at runtime via the option
4988: $ -snes_type my_solver
4990: Level: advanced
4992: Note: If your function is not being put into a shared library then use SNESRegister() instead
4994: .seealso: SNESRegisterAll(), SNESRegisterDestroy()
4996: Level: advanced
4997: @*/
4998: PetscErrorCode SNESRegister(const char sname[],PetscErrorCode (*function)(SNES))
4999: {
5003: SNESInitializePackage();
5004: PetscFunctionListAdd(&SNESList,sname,function);
5005: return(0);
5006: }
5008: PetscErrorCode SNESTestLocalMin(SNES snes)
5009: {
5011: PetscInt N,i,j;
5012: Vec u,uh,fh;
5013: PetscScalar value;
5014: PetscReal norm;
5017: SNESGetSolution(snes,&u);
5018: VecDuplicate(u,&uh);
5019: VecDuplicate(u,&fh);
5021: /* currently only works for sequential */
5022: PetscPrintf(PetscObjectComm((PetscObject)snes),"Testing FormFunction() for local min\n");
5023: VecGetSize(u,&N);
5024: for (i=0; i<N; i++) {
5025: VecCopy(u,uh);
5026: PetscPrintf(PetscObjectComm((PetscObject)snes),"i = %D\n",i);
5027: for (j=-10; j<11; j++) {
5028: value = PetscSign(j)*PetscExpReal(PetscAbs(j)-10.0);
5029: VecSetValue(uh,i,value,ADD_VALUES);
5030: SNESComputeFunction(snes,uh,fh);
5031: VecNorm(fh,NORM_2,&norm);
5032: PetscPrintf(PetscObjectComm((PetscObject)snes)," j norm %D %18.16e\n",j,norm);
5033: value = -value;
5034: VecSetValue(uh,i,value,ADD_VALUES);
5035: }
5036: }
5037: VecDestroy(&uh);
5038: VecDestroy(&fh);
5039: return(0);
5040: }
5042: /*@
5043: SNESKSPSetUseEW - Sets SNES use Eisenstat-Walker method for
5044: computing relative tolerance for linear solvers within an inexact
5045: Newton method.
5047: Logically Collective on SNES
5049: Input Parameters:
5050: + snes - SNES context
5051: - flag - PETSC_TRUE or PETSC_FALSE
5053: Options Database:
5054: + -snes_ksp_ew - use Eisenstat-Walker method for determining linear system convergence
5055: . -snes_ksp_ew_version ver - version of Eisenstat-Walker method
5056: . -snes_ksp_ew_rtol0 <rtol0> - Sets rtol0
5057: . -snes_ksp_ew_rtolmax <rtolmax> - Sets rtolmax
5058: . -snes_ksp_ew_gamma <gamma> - Sets gamma
5059: . -snes_ksp_ew_alpha <alpha> - Sets alpha
5060: . -snes_ksp_ew_alpha2 <alpha2> - Sets alpha2
5061: - -snes_ksp_ew_threshold <threshold> - Sets threshold
5063: Notes:
5064: Currently, the default is to use a constant relative tolerance for
5065: the inner linear solvers. Alternatively, one can use the
5066: Eisenstat-Walker method, where the relative convergence tolerance
5067: is reset at each Newton iteration according progress of the nonlinear
5068: solver.
5070: Level: advanced
5072: Reference:
5073: S. C. Eisenstat and H. F. Walker, "Choosing the forcing terms in an
5074: inexact Newton method", SISC 17 (1), pp.16-32, 1996.
5076: .seealso: SNESKSPGetUseEW(), SNESKSPGetParametersEW(), SNESKSPSetParametersEW()
5077: @*/
5078: PetscErrorCode SNESKSPSetUseEW(SNES snes,PetscBool flag)
5079: {
5083: snes->ksp_ewconv = flag;
5084: return(0);
5085: }
5087: /*@
5088: SNESKSPGetUseEW - Gets if SNES is using Eisenstat-Walker method
5089: for computing relative tolerance for linear solvers within an
5090: inexact Newton method.
5092: Not Collective
5094: Input Parameter:
5095: . snes - SNES context
5097: Output Parameter:
5098: . flag - PETSC_TRUE or PETSC_FALSE
5100: Notes:
5101: Currently, the default is to use a constant relative tolerance for
5102: the inner linear solvers. Alternatively, one can use the
5103: Eisenstat-Walker method, where the relative convergence tolerance
5104: is reset at each Newton iteration according progress of the nonlinear
5105: solver.
5107: Level: advanced
5109: Reference:
5110: S. C. Eisenstat and H. F. Walker, "Choosing the forcing terms in an
5111: inexact Newton method", SISC 17 (1), pp.16-32, 1996.
5113: .seealso: SNESKSPSetUseEW(), SNESKSPGetParametersEW(), SNESKSPSetParametersEW()
5114: @*/
5115: PetscErrorCode SNESKSPGetUseEW(SNES snes, PetscBool *flag)
5116: {
5120: *flag = snes->ksp_ewconv;
5121: return(0);
5122: }
5124: /*@
5125: SNESKSPSetParametersEW - Sets parameters for Eisenstat-Walker
5126: convergence criteria for the linear solvers within an inexact
5127: Newton method.
5129: Logically Collective on SNES
5131: Input Parameters:
5132: + snes - SNES context
5133: . version - version 1, 2 (default is 2) or 3
5134: . rtol_0 - initial relative tolerance (0 <= rtol_0 < 1)
5135: . rtol_max - maximum relative tolerance (0 <= rtol_max < 1)
5136: . gamma - multiplicative factor for version 2 rtol computation
5137: (0 <= gamma2 <= 1)
5138: . alpha - power for version 2 rtol computation (1 < alpha <= 2)
5139: . alpha2 - power for safeguard
5140: - threshold - threshold for imposing safeguard (0 < threshold < 1)
5142: Note:
5143: Version 3 was contributed by Luis Chacon, June 2006.
5145: Use PETSC_DEFAULT to retain the default for any of the parameters.
5147: Level: advanced
5149: Reference:
5150: S. C. Eisenstat and H. F. Walker, "Choosing the forcing terms in an
5151: inexact Newton method", Utah State University Math. Stat. Dept. Res.
5152: Report 6/94/75, June, 1994, to appear in SIAM J. Sci. Comput.
5154: .seealso: SNESKSPSetUseEW(), SNESKSPGetUseEW(), SNESKSPGetParametersEW()
5155: @*/
5156: PetscErrorCode SNESKSPSetParametersEW(SNES snes,PetscInt version,PetscReal rtol_0,PetscReal rtol_max,PetscReal gamma,PetscReal alpha,PetscReal alpha2,PetscReal threshold)
5157: {
5158: SNESKSPEW *kctx;
5162: kctx = (SNESKSPEW*)snes->kspconvctx;
5163: if (!kctx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"No Eisenstat-Walker context existing");
5172: if (version != PETSC_DEFAULT) kctx->version = version;
5173: if (rtol_0 != PETSC_DEFAULT) kctx->rtol_0 = rtol_0;
5174: if (rtol_max != PETSC_DEFAULT) kctx->rtol_max = rtol_max;
5175: if (gamma != PETSC_DEFAULT) kctx->gamma = gamma;
5176: if (alpha != PETSC_DEFAULT) kctx->alpha = alpha;
5177: if (alpha2 != PETSC_DEFAULT) kctx->alpha2 = alpha2;
5178: if (threshold != PETSC_DEFAULT) kctx->threshold = threshold;
5180: if (kctx->version < 1 || kctx->version > 3) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Only versions 1, 2 and 3 are supported: %D",kctx->version);
5181: if (kctx->rtol_0 < 0.0 || kctx->rtol_0 >= 1.0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"0.0 <= rtol_0 < 1.0: %g",(double)kctx->rtol_0);
5182: if (kctx->rtol_max < 0.0 || kctx->rtol_max >= 1.0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"0.0 <= rtol_max (%g) < 1.0\n",(double)kctx->rtol_max);
5183: if (kctx->gamma < 0.0 || kctx->gamma > 1.0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"0.0 <= gamma (%g) <= 1.0\n",(double)kctx->gamma);
5184: if (kctx->alpha <= 1.0 || kctx->alpha > 2.0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"1.0 < alpha (%g) <= 2.0\n",(double)kctx->alpha);
5185: if (kctx->threshold <= 0.0 || kctx->threshold >= 1.0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"0.0 < threshold (%g) < 1.0\n",(double)kctx->threshold);
5186: return(0);
5187: }
5189: /*@
5190: SNESKSPGetParametersEW - Gets parameters for Eisenstat-Walker
5191: convergence criteria for the linear solvers within an inexact
5192: Newton method.
5194: Not Collective
5196: Input Parameters:
5197: snes - SNES context
5199: Output Parameters:
5200: + version - version 1, 2 (default is 2) or 3
5201: . rtol_0 - initial relative tolerance (0 <= rtol_0 < 1)
5202: . rtol_max - maximum relative tolerance (0 <= rtol_max < 1)
5203: . gamma - multiplicative factor for version 2 rtol computation (0 <= gamma2 <= 1)
5204: . alpha - power for version 2 rtol computation (1 < alpha <= 2)
5205: . alpha2 - power for safeguard
5206: - threshold - threshold for imposing safeguard (0 < threshold < 1)
5208: Level: advanced
5210: .seealso: SNESKSPSetUseEW(), SNESKSPGetUseEW(), SNESKSPSetParametersEW()
5211: @*/
5212: PetscErrorCode SNESKSPGetParametersEW(SNES snes,PetscInt *version,PetscReal *rtol_0,PetscReal *rtol_max,PetscReal *gamma,PetscReal *alpha,PetscReal *alpha2,PetscReal *threshold)
5213: {
5214: SNESKSPEW *kctx;
5218: kctx = (SNESKSPEW*)snes->kspconvctx;
5219: if (!kctx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"No Eisenstat-Walker context existing");
5220: if (version) *version = kctx->version;
5221: if (rtol_0) *rtol_0 = kctx->rtol_0;
5222: if (rtol_max) *rtol_max = kctx->rtol_max;
5223: if (gamma) *gamma = kctx->gamma;
5224: if (alpha) *alpha = kctx->alpha;
5225: if (alpha2) *alpha2 = kctx->alpha2;
5226: if (threshold) *threshold = kctx->threshold;
5227: return(0);
5228: }
5230: PetscErrorCode KSPPreSolve_SNESEW(KSP ksp, Vec b, Vec x, SNES snes)
5231: {
5233: SNESKSPEW *kctx = (SNESKSPEW*)snes->kspconvctx;
5234: PetscReal rtol = PETSC_DEFAULT,stol;
5237: if (!snes->ksp_ewconv) return(0);
5238: if (!snes->iter) {
5239: rtol = kctx->rtol_0; /* first time in, so use the original user rtol */
5240: VecNorm(snes->vec_func,NORM_2,&kctx->norm_first);
5241: }
5242: else {
5243: if (kctx->version == 1) {
5244: rtol = (snes->norm - kctx->lresid_last)/kctx->norm_last;
5245: if (rtol < 0.0) rtol = -rtol;
5246: stol = PetscPowReal(kctx->rtol_last,kctx->alpha2);
5247: if (stol > kctx->threshold) rtol = PetscMax(rtol,stol);
5248: } else if (kctx->version == 2) {
5249: rtol = kctx->gamma * PetscPowReal(snes->norm/kctx->norm_last,kctx->alpha);
5250: stol = kctx->gamma * PetscPowReal(kctx->rtol_last,kctx->alpha);
5251: if (stol > kctx->threshold) rtol = PetscMax(rtol,stol);
5252: } else if (kctx->version == 3) { /* contributed by Luis Chacon, June 2006. */
5253: rtol = kctx->gamma * PetscPowReal(snes->norm/kctx->norm_last,kctx->alpha);
5254: /* safeguard: avoid sharp decrease of rtol */
5255: stol = kctx->gamma*PetscPowReal(kctx->rtol_last,kctx->alpha);
5256: stol = PetscMax(rtol,stol);
5257: rtol = PetscMin(kctx->rtol_0,stol);
5258: /* safeguard: avoid oversolving */
5259: stol = kctx->gamma*(kctx->norm_first*snes->rtol)/snes->norm;
5260: stol = PetscMax(rtol,stol);
5261: rtol = PetscMin(kctx->rtol_0,stol);
5262: } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Only versions 1, 2 or 3 are supported: %D",kctx->version);
5263: }
5264: /* safeguard: avoid rtol greater than one */
5265: rtol = PetscMin(rtol,kctx->rtol_max);
5266: KSPSetTolerances(ksp,rtol,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
5267: PetscInfo3(snes,"iter %D, Eisenstat-Walker (version %D) KSP rtol=%g\n",snes->iter,kctx->version,(double)rtol);
5268: return(0);
5269: }
5271: PetscErrorCode KSPPostSolve_SNESEW(KSP ksp, Vec b, Vec x, SNES snes)
5272: {
5274: SNESKSPEW *kctx = (SNESKSPEW*)snes->kspconvctx;
5275: PCSide pcside;
5276: Vec lres;
5279: if (!snes->ksp_ewconv) return(0);
5280: KSPGetTolerances(ksp,&kctx->rtol_last,NULL,NULL,NULL);
5281: kctx->norm_last = snes->norm;
5282: if (kctx->version == 1) {
5283: PC pc;
5284: PetscBool isNone;
5286: KSPGetPC(ksp, &pc);
5287: PetscObjectTypeCompare((PetscObject) pc, PCNONE, &isNone);
5288: KSPGetPCSide(ksp,&pcside);
5289: if (pcside == PC_RIGHT || isNone) { /* XXX Should we also test KSP_UNPRECONDITIONED_NORM ? */
5290: /* KSP residual is true linear residual */
5291: KSPGetResidualNorm(ksp,&kctx->lresid_last);
5292: } else {
5293: /* KSP residual is preconditioned residual */
5294: /* compute true linear residual norm */
5295: VecDuplicate(b,&lres);
5296: MatMult(snes->jacobian,x,lres);
5297: VecAYPX(lres,-1.0,b);
5298: VecNorm(lres,NORM_2,&kctx->lresid_last);
5299: VecDestroy(&lres);
5300: }
5301: }
5302: return(0);
5303: }
5305: /*@
5306: SNESGetKSP - Returns the KSP context for a SNES solver.
5308: Not Collective, but if SNES object is parallel, then KSP object is parallel
5310: Input Parameter:
5311: . snes - the SNES context
5313: Output Parameter:
5314: . ksp - the KSP context
5316: Notes:
5317: The user can then directly manipulate the KSP context to set various
5318: options, etc. Likewise, the user can then extract and manipulate the
5319: PC contexts as well.
5321: Level: beginner
5323: .seealso: KSPGetPC(), SNESCreate(), KSPCreate(), SNESSetKSP()
5324: @*/
5325: PetscErrorCode SNESGetKSP(SNES snes,KSP *ksp)
5326: {
5333: if (!snes->ksp) {
5334: PetscBool monitor = PETSC_FALSE;
5336: KSPCreate(PetscObjectComm((PetscObject)snes),&snes->ksp);
5337: PetscObjectIncrementTabLevel((PetscObject)snes->ksp,(PetscObject)snes,1);
5338: PetscLogObjectParent((PetscObject)snes,(PetscObject)snes->ksp);
5340: KSPSetPreSolve(snes->ksp,(PetscErrorCode (*)(KSP,Vec,Vec,void*))KSPPreSolve_SNESEW,snes);
5341: KSPSetPostSolve(snes->ksp,(PetscErrorCode (*)(KSP,Vec,Vec,void*))KSPPostSolve_SNESEW,snes);
5343: PetscOptionsGetBool(((PetscObject)snes)->options,((PetscObject)snes)->prefix,"-ksp_monitor_snes",&monitor,NULL);
5344: if (monitor) {
5345: KSPMonitorSet(snes->ksp,KSPMonitorSNES,snes,NULL);
5346: }
5347: monitor = PETSC_FALSE;
5348: PetscOptionsGetBool(((PetscObject)snes)->options,((PetscObject)snes)->prefix,"-ksp_monitor_snes_lg",&monitor,NULL);
5349: if (monitor) {
5350: PetscObject *objs;
5351: KSPMonitorSNESLGResidualNormCreate(PetscObjectComm((PetscObject)snes),NULL,NULL,PETSC_DECIDE,PETSC_DECIDE,600,600,&objs);
5352: objs[0] = (PetscObject) snes;
5353: KSPMonitorSet(snes->ksp,(PetscErrorCode (*)(KSP,PetscInt,PetscReal,void*))KSPMonitorSNESLGResidualNorm,objs,(PetscErrorCode (*)(void**))KSPMonitorSNESLGResidualNormDestroy);
5354: }
5355: PetscObjectSetOptions((PetscObject)snes->ksp,((PetscObject)snes)->options);
5356: }
5357: *ksp = snes->ksp;
5358: return(0);
5359: }
5362: #include <petsc/private/dmimpl.h>
5363: /*@
5364: SNESSetDM - Sets the DM that may be used by some nonlinear solvers or their underlying preconditioners
5366: Logically Collective on SNES
5368: Input Parameters:
5369: + snes - the nonlinear solver context
5370: - dm - the dm, cannot be NULL
5372: Notes:
5373: A DM can only be used for solving one problem at a time because information about the problem is stored on the DM,
5374: even when not using interfaces like DMSNESSetFunction(). Use DMClone() to get a distinct DM when solving different
5375: problems using the same function space.
5377: Level: intermediate
5379: .seealso: SNESGetDM(), KSPSetDM(), KSPGetDM()
5380: @*/
5381: PetscErrorCode SNESSetDM(SNES snes,DM dm)
5382: {
5384: KSP ksp;
5385: DMSNES sdm;
5390: PetscObjectReference((PetscObject)dm);
5391: if (snes->dm) { /* Move the DMSNES context over to the new DM unless the new DM already has one */
5392: if (snes->dm->dmsnes && !dm->dmsnes) {
5393: DMCopyDMSNES(snes->dm,dm);
5394: DMGetDMSNES(snes->dm,&sdm);
5395: if (sdm->originaldm == snes->dm) sdm->originaldm = dm; /* Grant write privileges to the replacement DM */
5396: }
5397: DMCoarsenHookRemove(snes->dm,DMCoarsenHook_SNESVecSol,DMRestrictHook_SNESVecSol,snes);
5398: DMDestroy(&snes->dm);
5399: }
5400: snes->dm = dm;
5401: snes->dmAuto = PETSC_FALSE;
5403: SNESGetKSP(snes,&ksp);
5404: KSPSetDM(ksp,dm);
5405: KSPSetDMActive(ksp,PETSC_FALSE);
5406: if (snes->npc) {
5407: SNESSetDM(snes->npc, snes->dm);
5408: SNESSetNPCSide(snes,snes->npcside);
5409: }
5410: return(0);
5411: }
5413: /*@
5414: SNESGetDM - Gets the DM that may be used by some preconditioners
5416: Not Collective but DM obtained is parallel on SNES
5418: Input Parameter:
5419: . snes - the preconditioner context
5421: Output Parameter:
5422: . dm - the dm
5424: Level: intermediate
5426: .seealso: SNESSetDM(), KSPSetDM(), KSPGetDM()
5427: @*/
5428: PetscErrorCode SNESGetDM(SNES snes,DM *dm)
5429: {
5434: if (!snes->dm) {
5435: DMShellCreate(PetscObjectComm((PetscObject)snes),&snes->dm);
5436: snes->dmAuto = PETSC_TRUE;
5437: }
5438: *dm = snes->dm;
5439: return(0);
5440: }
5442: /*@
5443: SNESSetNPC - Sets the nonlinear preconditioner to be used.
5445: Collective on SNES
5447: Input Parameters:
5448: + snes - iterative context obtained from SNESCreate()
5449: - pc - the preconditioner object
5451: Notes:
5452: Use SNESGetNPC() to retrieve the preconditioner context (for example,
5453: to configure it using the API).
5455: Level: developer
5457: .seealso: SNESGetNPC(), SNESHasNPC()
5458: @*/
5459: PetscErrorCode SNESSetNPC(SNES snes, SNES pc)
5460: {
5467: PetscObjectReference((PetscObject) pc);
5468: SNESDestroy(&snes->npc);
5469: snes->npc = pc;
5470: PetscLogObjectParent((PetscObject)snes, (PetscObject)snes->npc);
5471: return(0);
5472: }
5474: /*@
5475: SNESGetNPC - Creates a nonlinear preconditioning solver (SNES) to be used to precondition the nonlinear solver.
5477: Not Collective; but any changes to the obtained SNES object must be applied collectively
5479: Input Parameter:
5480: . snes - iterative context obtained from SNESCreate()
5482: Output Parameter:
5483: . pc - preconditioner context
5485: Options Database:
5486: . -npc_snes_type <type> - set the type of the SNES to use as the nonlinear preconditioner
5488: Notes:
5489: If a SNES was previously set with SNESSetNPC() then that SNES is returned, otherwise a new SNES object is created.
5491: The (preconditioner) SNES returned automatically inherits the same nonlinear function and Jacobian supplied to the original
5492: SNES during SNESSetUp()
5494: Level: developer
5496: .seealso: SNESSetNPC(), SNESHasNPC(), SNES, SNESCreate()
5497: @*/
5498: PetscErrorCode SNESGetNPC(SNES snes, SNES *pc)
5499: {
5501: const char *optionsprefix;
5506: if (!snes->npc) {
5507: SNESCreate(PetscObjectComm((PetscObject)snes),&snes->npc);
5508: PetscObjectIncrementTabLevel((PetscObject)snes->npc,(PetscObject)snes,1);
5509: PetscLogObjectParent((PetscObject)snes,(PetscObject)snes->npc);
5510: SNESGetOptionsPrefix(snes,&optionsprefix);
5511: SNESSetOptionsPrefix(snes->npc,optionsprefix);
5512: SNESAppendOptionsPrefix(snes->npc,"npc_");
5513: SNESSetCountersReset(snes->npc,PETSC_FALSE);
5514: }
5515: *pc = snes->npc;
5516: return(0);
5517: }
5519: /*@
5520: SNESHasNPC - Returns whether a nonlinear preconditioner exists
5522: Not Collective
5524: Input Parameter:
5525: . snes - iterative context obtained from SNESCreate()
5527: Output Parameter:
5528: . has_npc - whether the SNES has an NPC or not
5530: Level: developer
5532: .seealso: SNESSetNPC(), SNESGetNPC()
5533: @*/
5534: PetscErrorCode SNESHasNPC(SNES snes, PetscBool *has_npc)
5535: {
5538: *has_npc = (PetscBool) (snes->npc ? PETSC_TRUE : PETSC_FALSE);
5539: return(0);
5540: }
5542: /*@
5543: SNESSetNPCSide - Sets the preconditioning side.
5545: Logically Collective on SNES
5547: Input Parameter:
5548: . snes - iterative context obtained from SNESCreate()
5550: Output Parameter:
5551: . side - the preconditioning side, where side is one of
5552: .vb
5553: PC_LEFT - left preconditioning
5554: PC_RIGHT - right preconditioning (default for most nonlinear solvers)
5555: .ve
5557: Options Database Keys:
5558: . -snes_pc_side <right,left>
5560: Notes:
5561: SNESNRICHARDSON and SNESNCG only support left preconditioning.
5563: Level: intermediate
5565: .seealso: SNESGetNPCSide(), KSPSetPCSide()
5566: @*/
5567: PetscErrorCode SNESSetNPCSide(SNES snes,PCSide side)
5568: {
5572: snes->npcside= side;
5573: return(0);
5574: }
5576: /*@
5577: SNESGetNPCSide - Gets the preconditioning side.
5579: Not Collective
5581: Input Parameter:
5582: . snes - iterative context obtained from SNESCreate()
5584: Output Parameter:
5585: . side - the preconditioning side, where side is one of
5586: .vb
5587: PC_LEFT - left preconditioning
5588: PC_RIGHT - right preconditioning (default for most nonlinear solvers)
5589: .ve
5591: Level: intermediate
5593: .seealso: SNESSetNPCSide(), KSPGetPCSide()
5594: @*/
5595: PetscErrorCode SNESGetNPCSide(SNES snes,PCSide *side)
5596: {
5600: *side = snes->npcside;
5601: return(0);
5602: }
5604: /*@
5605: SNESSetLineSearch - Sets the linesearch on the SNES instance.
5607: Collective on SNES
5609: Input Parameters:
5610: + snes - iterative context obtained from SNESCreate()
5611: - linesearch - the linesearch object
5613: Notes:
5614: Use SNESGetLineSearch() to retrieve the preconditioner context (for example,
5615: to configure it using the API).
5617: Level: developer
5619: .seealso: SNESGetLineSearch()
5620: @*/
5621: PetscErrorCode SNESSetLineSearch(SNES snes, SNESLineSearch linesearch)
5622: {
5629: PetscObjectReference((PetscObject) linesearch);
5630: SNESLineSearchDestroy(&snes->linesearch);
5632: snes->linesearch = linesearch;
5634: PetscLogObjectParent((PetscObject)snes, (PetscObject)snes->linesearch);
5635: return(0);
5636: }
5638: /*@
5639: SNESGetLineSearch - Returns a pointer to the line search context set with SNESSetLineSearch()
5640: or creates a default line search instance associated with the SNES and returns it.
5642: Not Collective
5644: Input Parameter:
5645: . snes - iterative context obtained from SNESCreate()
5647: Output Parameter:
5648: . linesearch - linesearch context
5650: Level: beginner
5652: .seealso: SNESSetLineSearch(), SNESLineSearchCreate()
5653: @*/
5654: PetscErrorCode SNESGetLineSearch(SNES snes, SNESLineSearch *linesearch)
5655: {
5657: const char *optionsprefix;
5662: if (!snes->linesearch) {
5663: SNESGetOptionsPrefix(snes, &optionsprefix);
5664: SNESLineSearchCreate(PetscObjectComm((PetscObject)snes), &snes->linesearch);
5665: SNESLineSearchSetSNES(snes->linesearch, snes);
5666: SNESLineSearchAppendOptionsPrefix(snes->linesearch, optionsprefix);
5667: PetscObjectIncrementTabLevel((PetscObject) snes->linesearch, (PetscObject) snes, 1);
5668: PetscLogObjectParent((PetscObject)snes, (PetscObject)snes->linesearch);
5669: }
5670: *linesearch = snes->linesearch;
5671: return(0);
5672: }