Actual source code: snes.c
petsc-3.11.4 2019-09-28
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_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: .keywords: SNES, set, initial guess, nonzero
35: .seealso: SNESGetErrorIfNotConverged(), KSPGetErrorIfNotConverged(), KSPSetErrorIFNotConverged()
36: @*/
37: PetscErrorCode SNESSetErrorIfNotConverged(SNES snes,PetscBool flg)
38: {
42: snes->errorifnotconverged = flg;
43: return(0);
44: }
46: /*@
47: SNESGetErrorIfNotConverged - Will SNESSolve() generate an error if the solver does not converge?
49: Not Collective
51: Input Parameter:
52: . snes - iterative context obtained from SNESCreate()
54: Output Parameter:
55: . flag - PETSC_TRUE if it will generate an error, else PETSC_FALSE
57: Level: intermediate
59: .keywords: SNES, set, initial guess, nonzero
61: .seealso: SNESSetErrorIfNotConverged(), KSPGetErrorIfNotConverged(), KSPSetErrorIFNotConverged()
62: @*/
63: PetscErrorCode SNESGetErrorIfNotConverged(SNES snes,PetscBool *flag)
64: {
68: *flag = snes->errorifnotconverged;
69: return(0);
70: }
72: /*@
73: SNESSetAlwaysComputesFinalResidual - does the SNES always compute the residual at the final solution?
75: Logically Collective on SNES
77: Input Parameters:
78: + snes - the shell SNES
79: - flg - is the residual computed?
81: Level: advanced
83: .seealso: SNESGetAlwaysComputesFinalResidual()
84: @*/
85: PetscErrorCode SNESSetAlwaysComputesFinalResidual(SNES snes, PetscBool flg)
86: {
89: snes->alwayscomputesfinalresidual = flg;
90: return(0);
91: }
93: /*@
94: SNESGetAlwaysComputesFinalResidual - does the SNES always compute the residual at the final solution?
96: Logically Collective on SNES
98: Input Parameter:
99: . snes - the shell SNES
101: Output Parameter:
102: . flg - is the residual computed?
104: Level: advanced
106: .seealso: SNESSetAlwaysComputesFinalResidual()
107: @*/
108: PetscErrorCode SNESGetAlwaysComputesFinalResidual(SNES snes, PetscBool *flg)
109: {
112: *flg = snes->alwayscomputesfinalresidual;
113: return(0);
114: }
116: /*@
117: SNESSetFunctionDomainError - tells SNES that the input vector to your SNESFunction is not
118: in the functions domain. For example, negative pressure.
120: Logically Collective on SNES
122: Input Parameters:
123: . snes - the SNES context
125: Level: advanced
127: .keywords: SNES, view
129: .seealso: SNESCreate(), SNESSetFunction(), SNESFunction
130: @*/
131: PetscErrorCode SNESSetFunctionDomainError(SNES snes)
132: {
135: if (snes->errorifnotconverged) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"User code indicates input vector is not in the function domain");
136: snes->domainerror = PETSC_TRUE;
137: return(0);
138: }
140: /*@
141: SNESSetJacobianDomainError - tells SNES that computeJacobian does not make sense any more. For example there is a negative element transformation.
143: Logically Collective on SNES
145: Input Parameters:
146: . snes - the SNES context
148: Level: advanced
150: .keywords: SNES, view
152: .seealso: SNESCreate(), SNESSetFunction(), SNESFunction(), SNESSetFunctionDomainError()
153: @*/
154: PetscErrorCode SNESSetJacobianDomainError(SNES snes)
155: {
158: if (snes->errorifnotconverged) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"User code indicates computeJacobian does not make sense");
159: snes->jacobiandomainerror = PETSC_TRUE;
160: return(0);
161: }
163: /*@
164: SNESSetCheckJacobianDomainError - if or not to check jacobian domain error after each Jacobian evaluation. By default, we check Jacobian domain error
165: in the debug mode, and do not check it in the optimized mode.
167: Logically Collective on SNES
169: Input Parameters:
170: . snes - the SNES context
171: . flg - indicates if or not to check jacobian domain error after each Jacobian evaluation
173: Level: advanced
175: .keywords: SNES, view
177: .seealso: SNESCreate(), SNESSetFunction(), SNESFunction(), SNESSetFunctionDomainError(), SNESGetCheckJacobianDomainError()
178: @*/
179: PetscErrorCode SNESSetCheckJacobianDomainError(SNES snes, PetscBool flg)
180: {
183: snes->checkjacdomainerror = flg;
184: return(0);
185: }
187: /*@
188: SNESGetCheckJacobianDomainError - Get an indicator whether or not we are checking Jacobian domain errors after each Jacobian evaluation.
190: Logically Collective on SNES
192: Input Parameters:
193: . snes - the SNES context
195: Output Parameters:
196: . flg - PETSC_FALSE indicates that we don't check jacobian domain errors after each Jacobian evaluation
198: Level: advanced
200: .keywords: SNES, view
202: .seealso: SNESCreate(), SNESSetFunction(), SNESFunction(), SNESSetFunctionDomainError(), SNESSetCheckJacobianDomainError()
203: @*/
204: PetscErrorCode SNESGetCheckJacobianDomainError(SNES snes, PetscBool *flg)
205: {
209: *flg = snes->checkjacdomainerror;
210: return(0);
211: }
213: /*@
214: SNESGetFunctionDomainError - Gets the status of the domain error after a call to SNESComputeFunction;
216: Logically Collective on SNES
218: Input Parameters:
219: . snes - the SNES context
221: Output Parameters:
222: . domainerror - Set to PETSC_TRUE if there's a domain error; PETSC_FALSE otherwise.
224: Level: advanced
226: .keywords: SNES, view
228: .seealso: SNESSetFunctionDomainError(), SNESComputeFunction()
229: @*/
230: PetscErrorCode SNESGetFunctionDomainError(SNES snes, PetscBool *domainerror)
231: {
235: *domainerror = snes->domainerror;
236: return(0);
237: }
239: /*@
240: SNESGetJacobianDomainError - Gets the status of the Jacobian domain error after a call to SNESComputeJacobian;
242: Logically Collective on SNES
244: Input Parameters:
245: . snes - the SNES context
247: Output Parameters:
248: . domainerror - Set to PETSC_TRUE if there's a jacobian domain error; PETSC_FALSE otherwise.
250: Level: advanced
252: .keywords: SNES, view
254: .seealso: SNESSetFunctionDomainError(), SNESComputeFunction(),SNESGetFunctionDomainError()
255: @*/
256: PetscErrorCode SNESGetJacobianDomainError(SNES snes, PetscBool *domainerror)
257: {
261: *domainerror = snes->jacobiandomainerror;
262: return(0);
263: }
265: /*@C
266: SNESLoad - Loads a SNES that has been stored in binary with SNESView().
268: Collective on PetscViewer
270: Input Parameters:
271: + newdm - the newly loaded SNES, this needs to have been created with SNESCreate() or
272: some related function before a call to SNESLoad().
273: - viewer - binary file viewer, obtained from PetscViewerBinaryOpen()
275: Level: intermediate
277: Notes:
278: The type is determined by the data in the file, any type set into the SNES before this call is ignored.
280: Notes for advanced users:
281: Most users should not need to know the details of the binary storage
282: format, since SNESLoad() and TSView() completely hide these details.
283: But for anyone who's interested, the standard binary matrix storage
284: format is
285: .vb
286: has not yet been determined
287: .ve
289: .seealso: PetscViewerBinaryOpen(), SNESView(), MatLoad(), VecLoad()
290: @*/
291: PetscErrorCode SNESLoad(SNES snes, PetscViewer viewer)
292: {
294: PetscBool isbinary;
295: PetscInt classid;
296: char type[256];
297: KSP ksp;
298: DM dm;
299: DMSNES dmsnes;
304: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
305: if (!isbinary) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Invalid viewer; open viewer with PetscViewerBinaryOpen()");
307: PetscViewerBinaryRead(viewer,&classid,1,NULL,PETSC_INT);
308: if (classid != SNES_FILE_CLASSID) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONG,"Not SNES next in file");
309: PetscViewerBinaryRead(viewer,type,256,NULL,PETSC_CHAR);
310: SNESSetType(snes, type);
311: if (snes->ops->load) {
312: (*snes->ops->load)(snes,viewer);
313: }
314: SNESGetDM(snes,&dm);
315: DMGetDMSNES(dm,&dmsnes);
316: DMSNESLoad(dmsnes,viewer);
317: SNESGetKSP(snes,&ksp);
318: KSPLoad(ksp,viewer);
319: return(0);
320: }
322: #include <petscdraw.h>
323: #if defined(PETSC_HAVE_SAWS)
324: #include <petscviewersaws.h>
325: #endif
327: /*@C
328: SNESView - Prints the SNES data structure.
330: Collective on SNES
332: Input Parameters:
333: + SNES - the SNES context
334: - viewer - visualization context
336: Options Database Key:
337: . -snes_view - Calls SNESView() at end of SNESSolve()
339: Notes:
340: The available visualization contexts include
341: + PETSC_VIEWER_STDOUT_SELF - standard output (default)
342: - PETSC_VIEWER_STDOUT_WORLD - synchronized standard
343: output where only the first processor opens
344: the file. All other processors send their
345: data to the first processor to print.
347: The user can open an alternative visualization context with
348: PetscViewerASCIIOpen() - output to a specified file.
350: Level: beginner
352: .keywords: SNES, view
354: .seealso: PetscViewerASCIIOpen()
355: @*/
356: PetscErrorCode SNESView(SNES snes,PetscViewer viewer)
357: {
358: SNESKSPEW *kctx;
360: KSP ksp;
361: SNESLineSearch linesearch;
362: PetscBool iascii,isstring,isbinary,isdraw;
363: DMSNES dmsnes;
364: #if defined(PETSC_HAVE_SAWS)
365: PetscBool issaws;
366: #endif
370: if (!viewer) {
371: PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)snes),&viewer);
372: }
376: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
377: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSTRING,&isstring);
378: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
379: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
380: #if defined(PETSC_HAVE_SAWS)
381: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSAWS,&issaws);
382: #endif
383: if (iascii) {
384: SNESNormSchedule normschedule;
385: DM dm;
386: PetscErrorCode (*cJ)(SNES,Vec,Mat,Mat,void*);
387: void *ctx;
389: PetscObjectPrintClassNamePrefixType((PetscObject)snes,viewer);
390: if (!snes->setupcalled) {
391: PetscViewerASCIIPrintf(viewer," SNES has not been set up so information may be incomplete\n");
392: }
393: if (snes->ops->view) {
394: PetscViewerASCIIPushTab(viewer);
395: (*snes->ops->view)(snes,viewer);
396: PetscViewerASCIIPopTab(viewer);
397: }
398: PetscViewerASCIIPrintf(viewer," maximum iterations=%D, maximum function evaluations=%D\n",snes->max_its,snes->max_funcs);
399: PetscViewerASCIIPrintf(viewer," tolerances: relative=%g, absolute=%g, solution=%g\n",(double)snes->rtol,(double)snes->abstol,(double)snes->stol);
400: if (snes->usesksp) {
401: PetscViewerASCIIPrintf(viewer," total number of linear solver iterations=%D\n",snes->linear_its);
402: }
403: PetscViewerASCIIPrintf(viewer," total number of function evaluations=%D\n",snes->nfuncs);
404: SNESGetNormSchedule(snes, &normschedule);
405: if (normschedule > 0) {PetscViewerASCIIPrintf(viewer," norm schedule %s\n",SNESNormSchedules[normschedule]);}
406: if (snes->gridsequence) {
407: PetscViewerASCIIPrintf(viewer," total number of grid sequence refinements=%D\n",snes->gridsequence);
408: }
409: if (snes->ksp_ewconv) {
410: kctx = (SNESKSPEW*)snes->kspconvctx;
411: if (kctx) {
412: PetscViewerASCIIPrintf(viewer," Eisenstat-Walker computation of KSP relative tolerance (version %D)\n",kctx->version);
413: PetscViewerASCIIPrintf(viewer," rtol_0=%g, rtol_max=%g, threshold=%g\n",(double)kctx->rtol_0,(double)kctx->rtol_max,(double)kctx->threshold);
414: PetscViewerASCIIPrintf(viewer," gamma=%g, alpha=%g, alpha2=%g\n",(double)kctx->gamma,(double)kctx->alpha,(double)kctx->alpha2);
415: }
416: }
417: if (snes->lagpreconditioner == -1) {
418: PetscViewerASCIIPrintf(viewer," Preconditioned is never rebuilt\n");
419: } else if (snes->lagpreconditioner > 1) {
420: PetscViewerASCIIPrintf(viewer," Preconditioned is rebuilt every %D new Jacobians\n",snes->lagpreconditioner);
421: }
422: if (snes->lagjacobian == -1) {
423: PetscViewerASCIIPrintf(viewer," Jacobian is never rebuilt\n");
424: } else if (snes->lagjacobian > 1) {
425: PetscViewerASCIIPrintf(viewer," Jacobian is rebuilt every %D SNES iterations\n",snes->lagjacobian);
426: }
427: SNESGetDM(snes,&dm);
428: DMSNESGetJacobian(dm,&cJ,&ctx);
429: if (cJ == SNESComputeJacobianDefault) {
430: PetscViewerASCIIPrintf(viewer," Jacobian is built using finite differences one column at a time\n");
431: } else if (cJ == SNESComputeJacobianDefaultColor) {
432: PetscViewerASCIIPrintf(viewer," Jacobian is built using finite differences with coloring\n");
433: }
434: } else if (isstring) {
435: const char *type;
436: SNESGetType(snes,&type);
437: PetscViewerStringSPrintf(viewer," %-3.3s",type);
438: } else if (isbinary) {
439: PetscInt classid = SNES_FILE_CLASSID;
440: MPI_Comm comm;
441: PetscMPIInt rank;
442: char type[256];
444: PetscObjectGetComm((PetscObject)snes,&comm);
445: MPI_Comm_rank(comm,&rank);
446: if (!rank) {
447: PetscViewerBinaryWrite(viewer,&classid,1,PETSC_INT,PETSC_FALSE);
448: PetscStrncpy(type,((PetscObject)snes)->type_name,sizeof(type));
449: PetscViewerBinaryWrite(viewer,type,sizeof(type),PETSC_CHAR,PETSC_FALSE);
450: }
451: if (snes->ops->view) {
452: (*snes->ops->view)(snes,viewer);
453: }
454: } else if (isdraw) {
455: PetscDraw draw;
456: char str[36];
457: PetscReal x,y,bottom,h;
459: PetscViewerDrawGetDraw(viewer,0,&draw);
460: PetscDrawGetCurrentPoint(draw,&x,&y);
461: PetscStrncpy(str,"SNES: ",sizeof(str));
462: PetscStrlcat(str,((PetscObject)snes)->type_name,sizeof(str));
463: PetscDrawStringBoxed(draw,x,y,PETSC_DRAW_BLUE,PETSC_DRAW_BLACK,str,NULL,&h);
464: bottom = y - h;
465: PetscDrawPushCurrentPoint(draw,x,bottom);
466: if (snes->ops->view) {
467: (*snes->ops->view)(snes,viewer);
468: }
469: #if defined(PETSC_HAVE_SAWS)
470: } else if (issaws) {
471: PetscMPIInt rank;
472: const char *name;
474: PetscObjectGetName((PetscObject)snes,&name);
475: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
476: if (!((PetscObject)snes)->amsmem && !rank) {
477: char dir[1024];
479: PetscObjectViewSAWs((PetscObject)snes,viewer);
480: PetscSNPrintf(dir,1024,"/PETSc/Objects/%s/its",name);
481: PetscStackCallSAWs(SAWs_Register,(dir,&snes->iter,1,SAWs_READ,SAWs_INT));
482: if (!snes->conv_hist) {
483: SNESSetConvergenceHistory(snes,NULL,NULL,PETSC_DECIDE,PETSC_TRUE);
484: }
485: PetscSNPrintf(dir,1024,"/PETSc/Objects/%s/conv_hist",name);
486: PetscStackCallSAWs(SAWs_Register,(dir,snes->conv_hist,10,SAWs_READ,SAWs_DOUBLE));
487: }
488: #endif
489: }
490: if (snes->linesearch) {
491: SNESGetLineSearch(snes, &linesearch);
492: PetscViewerASCIIPushTab(viewer);
493: SNESLineSearchView(linesearch, viewer);
494: PetscViewerASCIIPopTab(viewer);
495: }
496: if (snes->npc && snes->usesnpc) {
497: PetscViewerASCIIPushTab(viewer);
498: SNESView(snes->npc, viewer);
499: PetscViewerASCIIPopTab(viewer);
500: }
501: PetscViewerASCIIPushTab(viewer);
502: DMGetDMSNES(snes->dm,&dmsnes);
503: DMSNESView(dmsnes, viewer);
504: PetscViewerASCIIPopTab(viewer);
505: if (snes->usesksp) {
506: SNESGetKSP(snes,&ksp);
507: PetscViewerASCIIPushTab(viewer);
508: KSPView(ksp,viewer);
509: PetscViewerASCIIPopTab(viewer);
510: }
511: if (isdraw) {
512: PetscDraw draw;
513: PetscViewerDrawGetDraw(viewer,0,&draw);
514: PetscDrawPopCurrentPoint(draw);
515: }
516: return(0);
517: }
519: /*
520: We retain a list of functions that also take SNES command
521: line options. These are called at the end SNESSetFromOptions()
522: */
523: #define MAXSETFROMOPTIONS 5
524: static PetscInt numberofsetfromoptions;
525: static PetscErrorCode (*othersetfromoptions[MAXSETFROMOPTIONS])(SNES);
527: /*@C
528: SNESAddOptionsChecker - Adds an additional function to check for SNES options.
530: Not Collective
532: Input Parameter:
533: . snescheck - function that checks for options
535: Level: developer
537: .seealso: SNESSetFromOptions()
538: @*/
539: PetscErrorCode SNESAddOptionsChecker(PetscErrorCode (*snescheck)(SNES))
540: {
542: if (numberofsetfromoptions >= MAXSETFROMOPTIONS) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE, "Too many options checkers, only %D allowed", MAXSETFROMOPTIONS);
543: othersetfromoptions[numberofsetfromoptions++] = snescheck;
544: return(0);
545: }
547: PETSC_INTERN PetscErrorCode SNESDefaultMatrixFreeCreate2(SNES,Vec,Mat*);
549: static PetscErrorCode SNESSetUpMatrixFree_Private(SNES snes, PetscBool hasOperator, PetscInt version)
550: {
551: Mat J;
552: KSP ksp;
553: PC pc;
554: PetscBool match;
556: MatNullSpace nullsp;
561: if (!snes->vec_func && (snes->jacobian || snes->jacobian_pre)) {
562: Mat A = snes->jacobian, B = snes->jacobian_pre;
563: MatCreateVecs(A ? A : B, NULL,&snes->vec_func);
564: }
566: if (version == 1) {
567: MatCreateSNESMF(snes,&J);
568: MatMFFDSetOptionsPrefix(J,((PetscObject)snes)->prefix);
569: MatSetFromOptions(J);
570: } else if (version == 2) {
571: if (!snes->vec_func) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"SNESSetFunction() must be called first");
572: #if !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_REAL_SINGLE) && !defined(PETSC_USE_REAL___FLOAT128) && !defined(PETSC_USE_REAL___FP16)
573: SNESDefaultMatrixFreeCreate2(snes,snes->vec_func,&J);
574: #else
575: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP, "matrix-free operator rutines (version 2)");
576: #endif
577: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE, "matrix-free operator rutines, only version 1 and 2");
579: /* attach any user provided null space that was on Amat to the newly created matrix free matrix */
580: if (snes->jacobian) {
581: MatGetNullSpace(snes->jacobian,&nullsp);
582: if (nullsp) {
583: MatSetNullSpace(J,nullsp);
584: }
585: }
587: PetscInfo1(snes,"Setting default matrix-free operator routines (version %D)\n", version);
588: if (hasOperator) {
590: /* This version replaces the user provided Jacobian matrix with a
591: matrix-free version but still employs the user-provided preconditioner matrix. */
592: SNESSetJacobian(snes,J,0,0,0);
593: } else {
594: /* This version replaces both the user-provided Jacobian and the user-
595: provided preconditioner Jacobian with the default matrix free version. */
596: if ((snes->npcside== PC_LEFT) && snes->npc) {
597: if (!snes->jacobian){SNESSetJacobian(snes,J,0,0,0);}
598: } else {
599: SNESSetJacobian(snes,J,J,MatMFFDComputeJacobian,0);
600: }
601: /* Force no preconditioner */
602: SNESGetKSP(snes,&ksp);
603: KSPGetPC(ksp,&pc);
604: PetscObjectTypeCompare((PetscObject)pc,PCSHELL,&match);
605: if (!match) {
606: PetscInfo(snes,"Setting default matrix-free preconditioner routines\nThat is no preconditioner is being used\n");
607: PCSetType(pc,PCNONE);
608: }
609: }
610: MatDestroy(&J);
611: return(0);
612: }
614: static PetscErrorCode DMRestrictHook_SNESVecSol(DM dmfine,Mat Restrict,Vec Rscale,Mat Inject,DM dmcoarse,void *ctx)
615: {
616: SNES snes = (SNES)ctx;
618: Vec Xfine,Xfine_named = NULL,Xcoarse;
621: if (PetscLogPrintInfo) {
622: PetscInt finelevel,coarselevel,fineclevel,coarseclevel;
623: DMGetRefineLevel(dmfine,&finelevel);
624: DMGetCoarsenLevel(dmfine,&fineclevel);
625: DMGetRefineLevel(dmcoarse,&coarselevel);
626: DMGetCoarsenLevel(dmcoarse,&coarseclevel);
627: PetscInfo4(dmfine,"Restricting SNES solution vector from level %D-%D to level %D-%D\n",finelevel,fineclevel,coarselevel,coarseclevel);
628: }
629: if (dmfine == snes->dm) Xfine = snes->vec_sol;
630: else {
631: DMGetNamedGlobalVector(dmfine,"SNESVecSol",&Xfine_named);
632: Xfine = Xfine_named;
633: }
634: DMGetNamedGlobalVector(dmcoarse,"SNESVecSol",&Xcoarse);
635: if (Inject) {
636: MatRestrict(Inject,Xfine,Xcoarse);
637: } else {
638: MatRestrict(Restrict,Xfine,Xcoarse);
639: VecPointwiseMult(Xcoarse,Xcoarse,Rscale);
640: }
641: DMRestoreNamedGlobalVector(dmcoarse,"SNESVecSol",&Xcoarse);
642: if (Xfine_named) {DMRestoreNamedGlobalVector(dmfine,"SNESVecSol",&Xfine_named);}
643: return(0);
644: }
646: static PetscErrorCode DMCoarsenHook_SNESVecSol(DM dm,DM dmc,void *ctx)
647: {
651: DMCoarsenHookAdd(dmc,DMCoarsenHook_SNESVecSol,DMRestrictHook_SNESVecSol,ctx);
652: return(0);
653: }
655: /* This may be called to rediscretize the operator on levels of linear multigrid. The DM shuffle is so the user can
656: * safely call SNESGetDM() in their residual evaluation routine. */
657: static PetscErrorCode KSPComputeOperators_SNES(KSP ksp,Mat A,Mat B,void *ctx)
658: {
659: SNES snes = (SNES)ctx;
661: Vec X,Xnamed = NULL;
662: DM dmsave;
663: void *ctxsave;
664: PetscErrorCode (*jac)(SNES,Vec,Mat,Mat,void*) = NULL;
667: dmsave = snes->dm;
668: KSPGetDM(ksp,&snes->dm);
669: if (dmsave == snes->dm) X = snes->vec_sol; /* We are on the finest level */
670: else { /* We are on a coarser level, this vec was initialized using a DM restrict hook */
671: DMGetNamedGlobalVector(snes->dm,"SNESVecSol",&Xnamed);
672: X = Xnamed;
673: SNESGetJacobian(snes,NULL,NULL,&jac,&ctxsave);
674: /* If the DM's don't match up, the MatFDColoring context needed for the jacobian won't match up either -- fixit. */
675: if (jac == SNESComputeJacobianDefaultColor) {
676: SNESSetJacobian(snes,NULL,NULL,SNESComputeJacobianDefaultColor,0);
677: }
678: }
679: /* Make sure KSP DM has the Jacobian computation routine */
680: {
681: DMSNES sdm;
683: DMGetDMSNES(snes->dm, &sdm);
684: if (!sdm->ops->computejacobian) {
685: DMCopyDMSNES(dmsave, snes->dm);
686: }
687: }
688: /* Compute the operators */
689: SNESComputeJacobian(snes,X,A,B);
690: /* Put the previous context back */
691: if (snes->dm != dmsave && jac == SNESComputeJacobianDefaultColor) {
692: SNESSetJacobian(snes,NULL,NULL,jac,ctxsave);
693: }
695: if (Xnamed) {DMRestoreNamedGlobalVector(snes->dm,"SNESVecSol",&Xnamed);}
696: snes->dm = dmsave;
697: return(0);
698: }
700: /*@
701: SNESSetUpMatrices - ensures that matrices are available for SNES, to be called by SNESSetUp_XXX()
703: Collective
705: Input Arguments:
706: . snes - snes to configure
708: Level: developer
710: .seealso: SNESSetUp()
711: @*/
712: PetscErrorCode SNESSetUpMatrices(SNES snes)
713: {
715: DM dm;
716: DMSNES sdm;
719: SNESGetDM(snes,&dm);
720: DMGetDMSNES(dm,&sdm);
721: if (!sdm->ops->computejacobian) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_PLIB,"DMSNES not properly configured");
722: else if (!snes->jacobian && snes->mf) {
723: Mat J;
724: void *functx;
725: MatCreateSNESMF(snes,&J);
726: MatMFFDSetOptionsPrefix(J,((PetscObject)snes)->prefix);
727: MatSetFromOptions(J);
728: SNESGetFunction(snes,NULL,NULL,&functx);
729: SNESSetJacobian(snes,J,J,0,0);
730: MatDestroy(&J);
731: } else if (snes->mf_operator && !snes->jacobian_pre && !snes->jacobian) {
732: Mat J,B;
733: MatCreateSNESMF(snes,&J);
734: MatMFFDSetOptionsPrefix(J,((PetscObject)snes)->prefix);
735: MatSetFromOptions(J);
736: DMCreateMatrix(snes->dm,&B);
737: /* sdm->computejacobian was already set to reach here */
738: SNESSetJacobian(snes,J,B,NULL,NULL);
739: MatDestroy(&J);
740: MatDestroy(&B);
741: } else if (!snes->jacobian_pre) {
742: PetscErrorCode (*nspconstr)(DM, PetscInt, MatNullSpace *);
743: PetscDS prob;
744: Mat J, B;
745: MatNullSpace nullspace = NULL;
746: PetscBool hasPrec = PETSC_FALSE;
747: PetscInt Nf;
749: J = snes->jacobian;
750: DMGetDS(dm, &prob);
751: if (prob) {PetscDSHasJacobianPreconditioner(prob, &hasPrec);}
752: if (J) {PetscObjectReference((PetscObject) J);}
753: else if (hasPrec) {DMCreateMatrix(snes->dm, &J);}
754: DMCreateMatrix(snes->dm, &B);
755: PetscDSGetNumFields(prob, &Nf);
756: DMGetNullSpaceConstructor(snes->dm, Nf, &nspconstr);
757: if (nspconstr) (*nspconstr)(snes->dm, -1, &nullspace);
758: MatSetNullSpace(B, nullspace);
759: MatNullSpaceDestroy(&nullspace);
760: SNESSetJacobian(snes, J ? J : B, B, NULL, NULL);
761: MatDestroy(&J);
762: MatDestroy(&B);
763: }
764: {
765: KSP ksp;
766: SNESGetKSP(snes,&ksp);
767: KSPSetComputeOperators(ksp,KSPComputeOperators_SNES,snes);
768: DMCoarsenHookAdd(snes->dm,DMCoarsenHook_SNESVecSol,DMRestrictHook_SNESVecSol,snes);
769: }
770: return(0);
771: }
773: /*@C
774: SNESMonitorSetFromOptions - Sets a monitor function and viewer appropriate for the type indicated by the user
776: Collective on SNES
778: Input Parameters:
779: + snes - SNES object you wish to monitor
780: . name - the monitor type one is seeking
781: . help - message indicating what monitoring is done
782: . manual - manual page for the monitor
783: . monitor - the monitor function
784: - 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
786: Level: developer
788: .seealso: PetscOptionsGetViewer(), PetscOptionsGetReal(), PetscOptionsHasName(), PetscOptionsGetString(),
789: PetscOptionsGetIntArray(), PetscOptionsGetRealArray(), PetscOptionsBool()
790: PetscOptionsInt(), PetscOptionsString(), PetscOptionsReal(), PetscOptionsBool(),
791: PetscOptionsName(), PetscOptionsBegin(), PetscOptionsEnd(), PetscOptionsHead(),
792: PetscOptionsStringArray(),PetscOptionsRealArray(), PetscOptionsScalar(),
793: PetscOptionsBoolGroupBegin(), PetscOptionsBoolGroup(), PetscOptionsBoolGroupEnd(),
794: PetscOptionsFList(), PetscOptionsEList()
795: @*/
796: PetscErrorCode SNESMonitorSetFromOptions(SNES snes,const char name[],const char help[], const char manual[],PetscErrorCode (*monitor)(SNES,PetscInt,PetscReal,PetscViewerAndFormat*),PetscErrorCode (*monitorsetup)(SNES,PetscViewerAndFormat*))
797: {
798: PetscErrorCode ierr;
799: PetscViewer viewer;
800: PetscViewerFormat format;
801: PetscBool flg;
804: PetscOptionsGetViewer(PetscObjectComm((PetscObject)snes),((PetscObject) snes)->options,((PetscObject)snes)->prefix,name,&viewer,&format,&flg);
805: if (flg) {
806: PetscViewerAndFormat *vf;
807: PetscViewerAndFormatCreate(viewer,format,&vf);
808: PetscObjectDereference((PetscObject)viewer);
809: if (monitorsetup) {
810: (*monitorsetup)(snes,vf);
811: }
812: SNESMonitorSet(snes,(PetscErrorCode (*)(SNES,PetscInt,PetscReal,void*))monitor,vf,(PetscErrorCode (*)(void**))PetscViewerAndFormatDestroy);
813: }
814: return(0);
815: }
817: /*@
818: SNESSetFromOptions - Sets various SNES and KSP parameters from user options.
820: Collective on SNES
822: Input Parameter:
823: . snes - the SNES context
825: Options Database Keys:
826: + -snes_type <type> - newtonls, newtontr, ngmres, ncg, nrichardson, qn, vi, fas, SNESType for complete list
827: . -snes_stol - convergence tolerance in terms of the norm
828: of the change in the solution between steps
829: . -snes_atol <abstol> - absolute tolerance of residual norm
830: . -snes_rtol <rtol> - relative decrease in tolerance norm from initial
831: . -snes_divergence_tolerance <divtol> - if the residual goes above divtol*rnorm0, exit with divergence
832: . -snes_force_iteration <force> - force SNESSolve() to take at least one iteration
833: . -snes_max_it <max_it> - maximum number of iterations
834: . -snes_max_funcs <max_funcs> - maximum number of function evaluations
835: . -snes_max_fail <max_fail> - maximum number of line search failures allowed before stopping, default is none
836: . -snes_max_linear_solve_fail - number of linear solver failures before SNESSolve() stops
837: . -snes_lag_preconditioner <lag> - how often preconditioner is rebuilt (use -1 to never rebuild)
838: . -snes_lag_jacobian <lag> - how often Jacobian is rebuilt (use -1 to never rebuild)
839: . -snes_trtol <trtol> - trust region tolerance
840: . -snes_no_convergence_test - skip convergence test in nonlinear
841: solver; hence iterations will continue until max_it
842: or some other criterion is reached. Saves expense
843: of convergence test
844: . -snes_monitor [ascii][:filename][:viewer format] - prints residual norm at each iteration. if no filename given prints to stdout
845: . -snes_monitor_solution [ascii binary draw][:filename][:viewer format] - plots solution at each iteration
846: . -snes_monitor_residual [ascii binary draw][:filename][:viewer format] - plots residual (not its norm) at each iteration
847: . -snes_monitor_solution_update [ascii binary draw][:filename][:viewer format] - plots update to solution at each iteration
848: . -snes_monitor_lg_residualnorm - plots residual norm at each iteration
849: . -snes_monitor_lg_range - plots residual norm at each iteration
850: . -snes_fd - use finite differences to compute Jacobian; very slow, only for testing
851: . -snes_fd_color - use finite differences with coloring to compute Jacobian
852: . -snes_mf_ksp_monitor - if using matrix-free multiply then print h at each KSP iteration
853: - -snes_converged_reason - print the reason for convergence/divergence after each solve
855: Options Database for Eisenstat-Walker method:
856: + -snes_ksp_ew - use Eisenstat-Walker method for determining linear system convergence
857: . -snes_ksp_ew_version ver - version of Eisenstat-Walker method
858: . -snes_ksp_ew_rtol0 <rtol0> - Sets rtol0
859: . -snes_ksp_ew_rtolmax <rtolmax> - Sets rtolmax
860: . -snes_ksp_ew_gamma <gamma> - Sets gamma
861: . -snes_ksp_ew_alpha <alpha> - Sets alpha
862: . -snes_ksp_ew_alpha2 <alpha2> - Sets alpha2
863: - -snes_ksp_ew_threshold <threshold> - Sets threshold
865: Notes:
866: To see all options, run your program with the -help option or consult
867: Users-Manual: ch_snes
869: Level: beginner
871: .keywords: SNES, nonlinear, set, options, database
873: .seealso: SNESSetOptionsPrefix(), SNESResetFromOptions()
874: @*/
875: PetscErrorCode SNESSetFromOptions(SNES snes)
876: {
877: PetscBool flg,pcset,persist,set;
878: PetscInt i,indx,lag,grids;
879: const char *deft = SNESNEWTONLS;
880: const char *convtests[] = {"default","skip"};
881: SNESKSPEW *kctx = NULL;
882: char type[256], monfilename[PETSC_MAX_PATH_LEN];
884: PCSide pcside;
885: const char *optionsprefix;
889: SNESRegisterAll();
890: PetscObjectOptionsBegin((PetscObject)snes);
891: if (((PetscObject)snes)->type_name) deft = ((PetscObject)snes)->type_name;
892: PetscOptionsFList("-snes_type","Nonlinear solver method","SNESSetType",SNESList,deft,type,256,&flg);
893: if (flg) {
894: SNESSetType(snes,type);
895: } else if (!((PetscObject)snes)->type_name) {
896: SNESSetType(snes,deft);
897: }
898: PetscOptionsReal("-snes_stol","Stop if step length less than","SNESSetTolerances",snes->stol,&snes->stol,NULL);
899: PetscOptionsReal("-snes_atol","Stop if function norm less than","SNESSetTolerances",snes->abstol,&snes->abstol,NULL);
901: PetscOptionsReal("-snes_rtol","Stop if decrease in function norm less than","SNESSetTolerances",snes->rtol,&snes->rtol,NULL);
902: PetscOptionsReal("-snes_divergence_tolerance","Stop if residual norm increases by this factor","SNESSetDivergenceTolerance",snes->divtol,&snes->divtol,NULL);
903: PetscOptionsInt("-snes_max_it","Maximum iterations","SNESSetTolerances",snes->max_its,&snes->max_its,NULL);
904: PetscOptionsInt("-snes_max_funcs","Maximum function evaluations","SNESSetTolerances",snes->max_funcs,&snes->max_funcs,NULL);
905: PetscOptionsInt("-snes_max_fail","Maximum nonlinear step failures","SNESSetMaxNonlinearStepFailures",snes->maxFailures,&snes->maxFailures,NULL);
906: PetscOptionsInt("-snes_max_linear_solve_fail","Maximum failures in linear solves allowed","SNESSetMaxLinearSolveFailures",snes->maxLinearSolveFailures,&snes->maxLinearSolveFailures,NULL);
907: PetscOptionsBool("-snes_error_if_not_converged","Generate error if solver does not converge","SNESSetErrorIfNotConverged",snes->errorifnotconverged,&snes->errorifnotconverged,NULL);
908: PetscOptionsBool("-snes_force_iteration","Force SNESSolve() to take at least one iteration","SNESSetForceIteration",snes->forceiteration,&snes->forceiteration,NULL);
909: PetscOptionsBool("-snes_check_jacobian_domain_error","Check Jacobian domain error after Jacobian evaluation","SNESCheckJacobianDomainError",snes->checkjacdomainerror,&snes->checkjacdomainerror,NULL);
911: PetscOptionsInt("-snes_lag_preconditioner","How often to rebuild preconditioner","SNESSetLagPreconditioner",snes->lagpreconditioner,&lag,&flg);
912: if (flg) {
913: SNESSetLagPreconditioner(snes,lag);
914: }
915: PetscOptionsBool("-snes_lag_preconditioner_persists","Preconditioner lagging through multiple solves","SNESSetLagPreconditionerPersists",snes->lagjac_persist,&persist,&flg);
916: if (flg) {
917: SNESSetLagPreconditionerPersists(snes,persist);
918: }
919: PetscOptionsInt("-snes_lag_jacobian","How often to rebuild Jacobian","SNESSetLagJacobian",snes->lagjacobian,&lag,&flg);
920: if (flg) {
921: SNESSetLagJacobian(snes,lag);
922: }
923: PetscOptionsBool("-snes_lag_jacobian_persists","Jacobian lagging through multiple solves","SNESSetLagJacobianPersists",snes->lagjac_persist,&persist,&flg);
924: if (flg) {
925: SNESSetLagJacobianPersists(snes,persist);
926: }
928: PetscOptionsInt("-snes_grid_sequence","Use grid sequencing to generate initial guess","SNESSetGridSequence",snes->gridsequence,&grids,&flg);
929: if (flg) {
930: SNESSetGridSequence(snes,grids);
931: }
933: PetscOptionsEList("-snes_convergence_test","Convergence test","SNESSetConvergenceTest",convtests,2,"default",&indx,&flg);
934: if (flg) {
935: switch (indx) {
936: case 0: SNESSetConvergenceTest(snes,SNESConvergedDefault,NULL,NULL); break;
937: case 1: SNESSetConvergenceTest(snes,SNESConvergedSkip,NULL,NULL); break;
938: }
939: }
941: PetscOptionsEList("-snes_norm_schedule","SNES Norm schedule","SNESSetNormSchedule",SNESNormSchedules,5,"function",&indx,&flg);
942: if (flg) { SNESSetNormSchedule(snes,(SNESNormSchedule)indx); }
944: PetscOptionsEList("-snes_function_type","SNES Norm schedule","SNESSetFunctionType",SNESFunctionTypes,2,"unpreconditioned",&indx,&flg);
945: if (flg) { SNESSetFunctionType(snes,(SNESFunctionType)indx); }
947: kctx = (SNESKSPEW*)snes->kspconvctx;
949: PetscOptionsBool("-snes_ksp_ew","Use Eisentat-Walker linear system convergence test","SNESKSPSetUseEW",snes->ksp_ewconv,&snes->ksp_ewconv,NULL);
951: PetscOptionsInt("-snes_ksp_ew_version","Version 1, 2 or 3","SNESKSPSetParametersEW",kctx->version,&kctx->version,NULL);
952: PetscOptionsReal("-snes_ksp_ew_rtol0","0 <= rtol0 < 1","SNESKSPSetParametersEW",kctx->rtol_0,&kctx->rtol_0,NULL);
953: PetscOptionsReal("-snes_ksp_ew_rtolmax","0 <= rtolmax < 1","SNESKSPSetParametersEW",kctx->rtol_max,&kctx->rtol_max,NULL);
954: PetscOptionsReal("-snes_ksp_ew_gamma","0 <= gamma <= 1","SNESKSPSetParametersEW",kctx->gamma,&kctx->gamma,NULL);
955: PetscOptionsReal("-snes_ksp_ew_alpha","1 < alpha <= 2","SNESKSPSetParametersEW",kctx->alpha,&kctx->alpha,NULL);
956: PetscOptionsReal("-snes_ksp_ew_alpha2","alpha2","SNESKSPSetParametersEW",kctx->alpha2,&kctx->alpha2,NULL);
957: PetscOptionsReal("-snes_ksp_ew_threshold","0 < threshold < 1","SNESKSPSetParametersEW",kctx->threshold,&kctx->threshold,NULL);
959: flg = PETSC_FALSE;
960: PetscOptionsBool("-snes_monitor_cancel","Remove all monitors","SNESMonitorCancel",flg,&flg,&set);
961: if (set && flg) {SNESMonitorCancel(snes);}
963: SNESMonitorSetFromOptions(snes,"-snes_monitor","Monitor norm of function","SNESMonitorDefault",SNESMonitorDefault,NULL);
964: SNESMonitorSetFromOptions(snes,"-snes_monitor_short","Monitor norm of function with fewer digits","SNESMonitorDefaultShort",SNESMonitorDefaultShort,NULL);
965: SNESMonitorSetFromOptions(snes,"-snes_monitor_range","Monitor range of elements of function","SNESMonitorRange",SNESMonitorRange,NULL);
967: SNESMonitorSetFromOptions(snes,"-snes_monitor_ratio","Monitor ratios of the norm of function for consecutive steps","SNESMonitorRatio",SNESMonitorRatio,SNESMonitorRatioSetUp);
968: SNESMonitorSetFromOptions(snes,"-snes_monitor_field","Monitor norm of function (split into fields)","SNESMonitorDefaultField",SNESMonitorDefaultField,NULL);
969: SNESMonitorSetFromOptions(snes,"-snes_monitor_solution","View solution at each iteration","SNESMonitorSolution",SNESMonitorSolution,NULL);
970: SNESMonitorSetFromOptions(snes,"-snes_monitor_solution_update","View correction at each iteration","SNESMonitorSolutionUpdate",SNESMonitorSolutionUpdate,NULL);
971: SNESMonitorSetFromOptions(snes,"-snes_monitor_residual","View residual at each iteration","SNESMonitorResidual",SNESMonitorResidual,NULL);
972: SNESMonitorSetFromOptions(snes,"-snes_monitor_jacupdate_spectrum","Print the change in the spectrum of the Jacobian","SNESMonitorJacUpdateSpectrum",SNESMonitorJacUpdateSpectrum,NULL);
973: SNESMonitorSetFromOptions(snes,"-snes_monitor_fields","Monitor norm of function per field","SNESMonitorSet",SNESMonitorFields,NULL);
975: PetscOptionsString("-snes_monitor_python","Use Python function","SNESMonitorSet",0,monfilename,PETSC_MAX_PATH_LEN,&flg);
976: if (flg) {PetscPythonMonitorSet((PetscObject)snes,monfilename);}
979: flg = PETSC_FALSE;
980: PetscOptionsBool("-snes_monitor_lg_residualnorm","Plot function norm at each iteration","SNESMonitorLGResidualNorm",flg,&flg,NULL);
981: if (flg) {
982: PetscDrawLG ctx;
984: SNESMonitorLGCreate(PetscObjectComm((PetscObject)snes),NULL,NULL,PETSC_DECIDE,PETSC_DECIDE,400,300,&ctx);
985: SNESMonitorSet(snes,SNESMonitorLGResidualNorm,ctx,(PetscErrorCode (*)(void**))PetscDrawLGDestroy);
986: }
987: flg = PETSC_FALSE;
988: PetscOptionsBool("-snes_monitor_lg_range","Plot function range at each iteration","SNESMonitorLGRange",flg,&flg,NULL);
989: if (flg) {
990: PetscViewer ctx;
992: PetscViewerDrawOpen(PetscObjectComm((PetscObject)snes),NULL,NULL,PETSC_DECIDE,PETSC_DECIDE,400,300,&ctx);
993: SNESMonitorSet(snes,SNESMonitorLGRange,ctx,(PetscErrorCode (*)(void**))PetscViewerDestroy);
994: }
998: flg = PETSC_FALSE;
999: PetscOptionsBool("-snes_fd","Use finite differences (slow) to compute Jacobian","SNESComputeJacobianDefault",flg,&flg,NULL);
1000: if (flg) {
1001: void *functx;
1002: DM dm;
1003: DMSNES sdm;
1004: SNESGetDM(snes,&dm);
1005: DMGetDMSNES(dm,&sdm);
1006: sdm->jacobianctx = NULL;
1007: SNESGetFunction(snes,NULL,NULL,&functx);
1008: SNESSetJacobian(snes,snes->jacobian,snes->jacobian_pre,SNESComputeJacobianDefault,functx);
1009: PetscInfo(snes,"Setting default finite difference Jacobian matrix\n");
1010: }
1012: flg = PETSC_FALSE;
1013: PetscOptionsBool("-snes_fd_function","Use finite differences (slow) to compute function from user objective","SNESObjectiveComputeFunctionDefaultFD",flg,&flg,NULL);
1014: if (flg) {
1015: SNESSetFunction(snes,NULL,SNESObjectiveComputeFunctionDefaultFD,NULL);
1016: }
1018: flg = PETSC_FALSE;
1019: PetscOptionsBool("-snes_fd_color","Use finite differences with coloring to compute Jacobian","SNESComputeJacobianDefaultColor",flg,&flg,NULL);
1020: if (flg) {
1021: DM dm;
1022: DMSNES sdm;
1023: SNESGetDM(snes,&dm);
1024: DMGetDMSNES(dm,&sdm);
1025: sdm->jacobianctx = NULL;
1026: SNESSetJacobian(snes,snes->jacobian,snes->jacobian_pre,SNESComputeJacobianDefaultColor,0);
1027: PetscInfo(snes,"Setting default finite difference coloring Jacobian matrix\n");
1028: }
1030: flg = PETSC_FALSE;
1031: PetscOptionsBool("-snes_mf_operator","Use a Matrix-Free Jacobian with user-provided preconditioner matrix","SNESSetUseMatrixFree",PETSC_FALSE,&snes->mf_operator,&flg);
1032: if (flg && snes->mf_operator) {
1033: snes->mf_operator = PETSC_TRUE;
1034: snes->mf = PETSC_TRUE;
1035: }
1036: flg = PETSC_FALSE;
1037: PetscOptionsBool("-snes_mf","Use a Matrix-Free Jacobian with no preconditioner matrix","SNESSetUseMatrixFree",PETSC_FALSE,&snes->mf,&flg);
1038: if (!flg && snes->mf_operator) snes->mf = PETSC_TRUE;
1039: PetscOptionsInt("-snes_mf_version","Matrix-Free routines version 1 or 2","None",snes->mf_version,&snes->mf_version,0);
1041: flg = PETSC_FALSE;
1042: SNESGetNPCSide(snes,&pcside);
1043: PetscOptionsEnum("-snes_npc_side","SNES nonlinear preconditioner side","SNESSetNPCSide",PCSides,(PetscEnum)pcside,(PetscEnum*)&pcside,&flg);
1044: if (flg) {SNESSetNPCSide(snes,pcside);}
1046: #if defined(PETSC_HAVE_SAWS)
1047: /*
1048: Publish convergence information using SAWs
1049: */
1050: flg = PETSC_FALSE;
1051: PetscOptionsBool("-snes_monitor_saws","Publish SNES progress using SAWs","SNESMonitorSet",flg,&flg,NULL);
1052: if (flg) {
1053: void *ctx;
1054: SNESMonitorSAWsCreate(snes,&ctx);
1055: SNESMonitorSet(snes,SNESMonitorSAWs,ctx,SNESMonitorSAWsDestroy);
1056: }
1057: #endif
1058: #if defined(PETSC_HAVE_SAWS)
1059: {
1060: PetscBool set;
1061: flg = PETSC_FALSE;
1062: PetscOptionsBool("-snes_saws_block","Block for SAWs at end of SNESSolve","PetscObjectSAWsBlock",((PetscObject)snes)->amspublishblock,&flg,&set);
1063: if (set) {
1064: PetscObjectSAWsSetBlock((PetscObject)snes,flg);
1065: }
1066: }
1067: #endif
1069: for (i = 0; i < numberofsetfromoptions; i++) {
1070: (*othersetfromoptions[i])(snes);
1071: }
1073: if (snes->ops->setfromoptions) {
1074: (*snes->ops->setfromoptions)(PetscOptionsObject,snes);
1075: }
1077: /* process any options handlers added with PetscObjectAddOptionsHandler() */
1078: PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject)snes);
1079: PetscOptionsEnd();
1081: if (!snes->linesearch) {
1082: SNESGetLineSearch(snes, &snes->linesearch);
1083: }
1084: SNESLineSearchSetFromOptions(snes->linesearch);
1086: if (!snes->ksp) {SNESGetKSP(snes,&snes->ksp);}
1087: KSPSetOperators(snes->ksp,snes->jacobian,snes->jacobian_pre);
1088: KSPSetFromOptions(snes->ksp);
1090: /* if someone has set the SNES NPC type, create it. */
1091: SNESGetOptionsPrefix(snes, &optionsprefix);
1092: PetscOptionsHasName(((PetscObject)snes)->options,optionsprefix, "-npc_snes_type", &pcset);
1093: if (pcset && (!snes->npc)) {
1094: SNESGetNPC(snes, &snes->npc);
1095: }
1096: snes->setfromoptionscalled++;
1097: return(0);
1098: }
1100: /*@
1101: SNESResetFromOptions - Sets various SNES and KSP parameters from user options ONLY if the SNES was previously set from options
1103: Collective on SNES
1105: Input Parameter:
1106: . snes - the SNES context
1108: Level: beginner
1110: .keywords: SNES, nonlinear, set, options, database
1112: .seealso: SNESSetFromOptions(), SNESSetOptionsPrefix()
1113: @*/
1114: PetscErrorCode SNESResetFromOptions(SNES snes)
1115: {
1119: if (snes->setfromoptionscalled) {SNESSetFromOptions(snes);}
1120: return(0);
1121: }
1123: /*@C
1124: SNESSetComputeApplicationContext - Sets an optional function to compute a user-defined context for
1125: the nonlinear solvers.
1127: Logically Collective on SNES
1129: Input Parameters:
1130: + snes - the SNES context
1131: . compute - function to compute the context
1132: - destroy - function to destroy the context
1134: Level: intermediate
1136: Notes:
1137: This function is currently not available from Fortran.
1139: .keywords: SNES, nonlinear, set, application, context
1141: .seealso: SNESGetApplicationContext(), SNESSetComputeApplicationContext(), SNESGetApplicationContext()
1142: @*/
1143: PetscErrorCode SNESSetComputeApplicationContext(SNES snes,PetscErrorCode (*compute)(SNES,void**),PetscErrorCode (*destroy)(void**))
1144: {
1147: snes->ops->usercompute = compute;
1148: snes->ops->userdestroy = destroy;
1149: return(0);
1150: }
1152: /*@
1153: SNESSetApplicationContext - Sets the optional user-defined context for
1154: the nonlinear solvers.
1156: Logically Collective on SNES
1158: Input Parameters:
1159: + snes - the SNES context
1160: - usrP - optional user context
1162: Level: intermediate
1164: Fortran Notes:
1165: To use this from Fortran you must write a Fortran interface definition for this
1166: function that tells Fortran the Fortran derived data type that you are passing in as the ctx argument.
1168: .keywords: SNES, nonlinear, set, application, context
1170: .seealso: SNESGetApplicationContext()
1171: @*/
1172: PetscErrorCode SNESSetApplicationContext(SNES snes,void *usrP)
1173: {
1175: KSP ksp;
1179: SNESGetKSP(snes,&ksp);
1180: KSPSetApplicationContext(ksp,usrP);
1181: snes->user = usrP;
1182: return(0);
1183: }
1185: /*@
1186: SNESGetApplicationContext - Gets the user-defined context for the
1187: nonlinear solvers.
1189: Not Collective
1191: Input Parameter:
1192: . snes - SNES context
1194: Output Parameter:
1195: . usrP - user context
1197: Fortran Notes:
1198: To use this from Fortran you must write a Fortran interface definition for this
1199: function that tells Fortran the Fortran derived data type that you are passing in as the ctx argument.
1201: Level: intermediate
1203: .keywords: SNES, nonlinear, get, application, context
1205: .seealso: SNESSetApplicationContext()
1206: @*/
1207: PetscErrorCode SNESGetApplicationContext(SNES snes,void *usrP)
1208: {
1211: *(void**)usrP = snes->user;
1212: return(0);
1213: }
1215: /*@
1216: SNESSetUseMatrixFree - indicates that SNES should use matrix free finite difference matrix vector products internally to apply
1217: the Jacobian.
1219: Collective on SNES
1221: Input Parameters:
1222: + snes - SNES context
1223: . mf - use matrix-free for both the Amat and Pmat used by SNESSetJacobian(), both the Amat and Pmat set in SNESSetJacobian() will be ignored
1224: - mf_operator - use matrix-free only for the Amat used by SNESSetJacobian(), this means the user provided Pmat will continue to be used
1226: Options Database:
1227: + -snes_mf - use matrix free for both the mat and pmat operator
1228: - -snes_mf_operator - use matrix free only for the mat operator
1230: Level: intermediate
1232: .keywords: SNES, nonlinear, get, iteration, number,
1234: .seealso: SNESGetUseMatrixFree(), MatCreateSNESMF()
1235: @*/
1236: PetscErrorCode SNESSetUseMatrixFree(SNES snes,PetscBool mf_operator,PetscBool mf)
1237: {
1242: if (mf && !mf_operator) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_INCOMP,"If using mf must also use mf_operator");
1243: snes->mf = mf;
1244: snes->mf_operator = mf_operator;
1245: return(0);
1246: }
1248: /*@
1249: SNESGetUseMatrixFree - indicates if the SNES uses matrix free finite difference matrix vector products to apply
1250: the Jacobian.
1252: Collective on SNES
1254: Input Parameter:
1255: . snes - SNES context
1257: Output Parameters:
1258: + mf - use matrix-free for both the Amat and Pmat used by SNESSetJacobian(), both the Amat and Pmat set in SNESSetJacobian() will be ignored
1259: - mf_operator - use matrix-free only for the Amat used by SNESSetJacobian(), this means the user provided Pmat will continue to be used
1261: Options Database:
1262: + -snes_mf - use matrix free for both the mat and pmat operator
1263: - -snes_mf_operator - use matrix free only for the mat operator
1265: Level: intermediate
1267: .keywords: SNES, nonlinear, get, iteration, number,
1269: .seealso: SNESSetUseMatrixFree(), MatCreateSNESMF()
1270: @*/
1271: PetscErrorCode SNESGetUseMatrixFree(SNES snes,PetscBool *mf_operator,PetscBool *mf)
1272: {
1275: if (mf) *mf = snes->mf;
1276: if (mf_operator) *mf_operator = snes->mf_operator;
1277: return(0);
1278: }
1280: /*@
1281: SNESGetIterationNumber - Gets the number of nonlinear iterations completed
1282: at this time.
1284: Not Collective
1286: Input Parameter:
1287: . snes - SNES context
1289: Output Parameter:
1290: . iter - iteration number
1292: Notes:
1293: For example, during the computation of iteration 2 this would return 1.
1295: This is useful for using lagged Jacobians (where one does not recompute the
1296: Jacobian at each SNES iteration). For example, the code
1297: .vb
1298: SNESGetIterationNumber(snes,&it);
1299: if (!(it % 2)) {
1300: [compute Jacobian here]
1301: }
1302: .ve
1303: can be used in your ComputeJacobian() function to cause the Jacobian to be
1304: recomputed every second SNES iteration.
1306: After the SNES solve is complete this will return the number of nonlinear iterations used.
1308: Level: intermediate
1310: .keywords: SNES, nonlinear, get, iteration, number,
1312: .seealso: SNESGetLinearSolveIterations()
1313: @*/
1314: PetscErrorCode SNESGetIterationNumber(SNES snes,PetscInt *iter)
1315: {
1319: *iter = snes->iter;
1320: return(0);
1321: }
1323: /*@
1324: SNESSetIterationNumber - Sets the current iteration number.
1326: Not Collective
1328: Input Parameter:
1329: . snes - SNES context
1330: . iter - iteration number
1332: Level: developer
1334: .keywords: SNES, nonlinear, set, iteration, number,
1336: .seealso: SNESGetLinearSolveIterations()
1337: @*/
1338: PetscErrorCode SNESSetIterationNumber(SNES snes,PetscInt iter)
1339: {
1344: PetscObjectSAWsTakeAccess((PetscObject)snes);
1345: snes->iter = iter;
1346: PetscObjectSAWsGrantAccess((PetscObject)snes);
1347: return(0);
1348: }
1350: /*@
1351: SNESGetNonlinearStepFailures - Gets the number of unsuccessful steps
1352: attempted by the nonlinear solver.
1354: Not Collective
1356: Input Parameter:
1357: . snes - SNES context
1359: Output Parameter:
1360: . nfails - number of unsuccessful steps attempted
1362: Notes:
1363: This counter is reset to zero for each successive call to SNESSolve().
1365: Level: intermediate
1367: .keywords: SNES, nonlinear, get, number, unsuccessful, steps
1369: .seealso: SNESGetMaxLinearSolveFailures(), SNESGetLinearSolveIterations(), SNESSetMaxLinearSolveFailures(), SNESGetLinearSolveFailures(),
1370: SNESSetMaxNonlinearStepFailures(), SNESGetMaxNonlinearStepFailures()
1371: @*/
1372: PetscErrorCode SNESGetNonlinearStepFailures(SNES snes,PetscInt *nfails)
1373: {
1377: *nfails = snes->numFailures;
1378: return(0);
1379: }
1381: /*@
1382: SNESSetMaxNonlinearStepFailures - Sets the maximum number of unsuccessful steps
1383: attempted by the nonlinear solver before it gives up.
1385: Not Collective
1387: Input Parameters:
1388: + snes - SNES context
1389: - maxFails - maximum of unsuccessful steps
1391: Level: intermediate
1393: .keywords: SNES, nonlinear, set, maximum, unsuccessful, steps
1395: .seealso: SNESGetMaxLinearSolveFailures(), SNESGetLinearSolveIterations(), SNESSetMaxLinearSolveFailures(), SNESGetLinearSolveFailures(),
1396: SNESGetMaxNonlinearStepFailures(), SNESGetNonlinearStepFailures()
1397: @*/
1398: PetscErrorCode SNESSetMaxNonlinearStepFailures(SNES snes, PetscInt maxFails)
1399: {
1402: snes->maxFailures = maxFails;
1403: return(0);
1404: }
1406: /*@
1407: SNESGetMaxNonlinearStepFailures - Gets the maximum number of unsuccessful steps
1408: attempted by the nonlinear solver before it gives up.
1410: Not Collective
1412: Input Parameter:
1413: . snes - SNES context
1415: Output Parameter:
1416: . maxFails - maximum of unsuccessful steps
1418: Level: intermediate
1420: .keywords: SNES, nonlinear, get, maximum, unsuccessful, steps
1422: .seealso: SNESGetMaxLinearSolveFailures(), SNESGetLinearSolveIterations(), SNESSetMaxLinearSolveFailures(), SNESGetLinearSolveFailures(),
1423: SNESSetMaxNonlinearStepFailures(), SNESGetNonlinearStepFailures()
1425: @*/
1426: PetscErrorCode SNESGetMaxNonlinearStepFailures(SNES snes, PetscInt *maxFails)
1427: {
1431: *maxFails = snes->maxFailures;
1432: return(0);
1433: }
1435: /*@
1436: SNESGetNumberFunctionEvals - Gets the number of user provided function evaluations
1437: done by SNES.
1439: Not Collective
1441: Input Parameter:
1442: . snes - SNES context
1444: Output Parameter:
1445: . nfuncs - number of evaluations
1447: Level: intermediate
1449: Notes:
1450: Reset every time SNESSolve is called unless SNESSetCountersReset() is used.
1452: .keywords: SNES, nonlinear, get, maximum, unsuccessful, steps
1454: .seealso: SNESGetMaxLinearSolveFailures(), SNESGetLinearSolveIterations(), SNESSetMaxLinearSolveFailures(), SNESGetLinearSolveFailures(), SNESSetCountersReset()
1455: @*/
1456: PetscErrorCode SNESGetNumberFunctionEvals(SNES snes, PetscInt *nfuncs)
1457: {
1461: *nfuncs = snes->nfuncs;
1462: return(0);
1463: }
1465: /*@
1466: SNESGetLinearSolveFailures - Gets the number of failed (non-converged)
1467: linear solvers.
1469: Not Collective
1471: Input Parameter:
1472: . snes - SNES context
1474: Output Parameter:
1475: . nfails - number of failed solves
1477: Level: intermediate
1479: Options Database Keys:
1480: . -snes_max_linear_solve_fail <num> - The number of failures before the solve is terminated
1482: Notes:
1483: This counter is reset to zero for each successive call to SNESSolve().
1485: .keywords: SNES, nonlinear, get, number, unsuccessful, steps
1487: .seealso: SNESGetMaxLinearSolveFailures(), SNESGetLinearSolveIterations(), SNESSetMaxLinearSolveFailures()
1488: @*/
1489: PetscErrorCode SNESGetLinearSolveFailures(SNES snes,PetscInt *nfails)
1490: {
1494: *nfails = snes->numLinearSolveFailures;
1495: return(0);
1496: }
1498: /*@
1499: SNESSetMaxLinearSolveFailures - the number of failed linear solve attempts
1500: allowed before SNES returns with a diverged reason of SNES_DIVERGED_LINEAR_SOLVE
1502: Logically Collective on SNES
1504: Input Parameters:
1505: + snes - SNES context
1506: - maxFails - maximum allowed linear solve failures
1508: Level: intermediate
1510: Options Database Keys:
1511: . -snes_max_linear_solve_fail <num> - The number of failures before the solve is terminated
1513: Notes:
1514: By default this is 0; that is SNES returns on the first failed linear solve
1516: .keywords: SNES, nonlinear, set, maximum, unsuccessful, steps
1518: .seealso: SNESGetLinearSolveFailures(), SNESGetMaxLinearSolveFailures(), SNESGetLinearSolveIterations()
1519: @*/
1520: PetscErrorCode SNESSetMaxLinearSolveFailures(SNES snes, PetscInt maxFails)
1521: {
1525: snes->maxLinearSolveFailures = maxFails;
1526: return(0);
1527: }
1529: /*@
1530: SNESGetMaxLinearSolveFailures - gets the maximum number of linear solve failures that
1531: are allowed before SNES terminates
1533: Not Collective
1535: Input Parameter:
1536: . snes - SNES context
1538: Output Parameter:
1539: . maxFails - maximum of unsuccessful solves allowed
1541: Level: intermediate
1543: Notes:
1544: By default this is 1; that is SNES returns on the first failed linear solve
1546: .keywords: SNES, nonlinear, get, maximum, unsuccessful, steps
1548: .seealso: SNESGetLinearSolveFailures(), SNESGetLinearSolveIterations(), SNESSetMaxLinearSolveFailures(),
1549: @*/
1550: PetscErrorCode SNESGetMaxLinearSolveFailures(SNES snes, PetscInt *maxFails)
1551: {
1555: *maxFails = snes->maxLinearSolveFailures;
1556: return(0);
1557: }
1559: /*@
1560: SNESGetLinearSolveIterations - Gets the total number of linear iterations
1561: used by the nonlinear solver.
1563: Not Collective
1565: Input Parameter:
1566: . snes - SNES context
1568: Output Parameter:
1569: . lits - number of linear iterations
1571: Notes:
1572: This counter is reset to zero for each successive call to SNESSolve() unless SNESSetCountersReset() is used.
1574: 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
1575: then call KSPGetIterationNumber() after the failed solve.
1577: Level: intermediate
1579: .keywords: SNES, nonlinear, get, number, linear, iterations
1581: .seealso: SNESGetIterationNumber(), SNESGetLinearSolveFailures(), SNESGetMaxLinearSolveFailures(), SNESSetCountersReset()
1582: @*/
1583: PetscErrorCode SNESGetLinearSolveIterations(SNES snes,PetscInt *lits)
1584: {
1588: *lits = snes->linear_its;
1589: return(0);
1590: }
1592: /*@
1593: SNESSetCountersReset - Sets whether or not the counters for linear iterations and function evaluations
1594: are reset every time SNESSolve() is called.
1596: Logically Collective on SNES
1598: Input Parameter:
1599: + snes - SNES context
1600: - reset - whether to reset the counters or not
1602: Notes:
1603: This defaults to PETSC_TRUE
1605: Level: developer
1607: .keywords: SNES, nonlinear, set, reset, number, linear, iterations
1609: .seealso: SNESGetNumberFunctionEvals(), SNESGetLinearSolveIterations(), SNESGetNPC()
1610: @*/
1611: PetscErrorCode SNESSetCountersReset(SNES snes,PetscBool reset)
1612: {
1616: snes->counters_reset = reset;
1617: return(0);
1618: }
1621: /*@
1622: SNESSetKSP - Sets a KSP context for the SNES object to use
1624: Not Collective, but the SNES and KSP objects must live on the same MPI_Comm
1626: Input Parameters:
1627: + snes - the SNES context
1628: - ksp - the KSP context
1630: Notes:
1631: The SNES object already has its KSP object, you can obtain with SNESGetKSP()
1632: so this routine is rarely needed.
1634: The KSP object that is already in the SNES object has its reference count
1635: decreased by one.
1637: Level: developer
1639: .keywords: SNES, nonlinear, get, KSP, context
1641: .seealso: KSPGetPC(), SNESCreate(), KSPCreate(), SNESSetKSP()
1642: @*/
1643: PetscErrorCode SNESSetKSP(SNES snes,KSP ksp)
1644: {
1651: PetscObjectReference((PetscObject)ksp);
1652: if (snes->ksp) {PetscObjectDereference((PetscObject)snes->ksp);}
1653: snes->ksp = ksp;
1654: return(0);
1655: }
1657: /* -----------------------------------------------------------*/
1658: /*@
1659: SNESCreate - Creates a nonlinear solver context.
1661: Collective on MPI_Comm
1663: Input Parameters:
1664: . comm - MPI communicator
1666: Output Parameter:
1667: . outsnes - the new SNES context
1669: Options Database Keys:
1670: + -snes_mf - Activates default matrix-free Jacobian-vector products,
1671: and no preconditioning matrix
1672: . -snes_mf_operator - Activates default matrix-free Jacobian-vector
1673: products, and a user-provided preconditioning matrix
1674: as set by SNESSetJacobian()
1675: - -snes_fd - Uses (slow!) finite differences to compute Jacobian
1677: Level: beginner
1679: Developer Notes:
1680: SNES always creates a KSP object even though many SNES methods do not use it. This is
1681: unfortunate and should be fixed at some point. The flag snes->usesksp indicates if the
1682: particular method does use KSP and regulates if the information about the KSP is printed
1683: in SNESView(). TSSetFromOptions() does call SNESSetFromOptions() which can lead to users being confused
1684: by help messages about meaningless SNES options.
1686: SNES always creates the snes->kspconvctx even though it is used by only one type. This should
1687: be fixed.
1689: .keywords: SNES, nonlinear, create, context
1691: .seealso: SNESSolve(), SNESDestroy(), SNES, SNESSetLagPreconditioner()
1693: @*/
1694: PetscErrorCode SNESCreate(MPI_Comm comm,SNES *outsnes)
1695: {
1697: SNES snes;
1698: SNESKSPEW *kctx;
1702: *outsnes = NULL;
1703: SNESInitializePackage();
1705: PetscHeaderCreate(snes,SNES_CLASSID,"SNES","Nonlinear solver","SNES",comm,SNESDestroy,SNESView);
1707: snes->ops->converged = SNESConvergedDefault;
1708: snes->usesksp = PETSC_TRUE;
1709: snes->tolerancesset = PETSC_FALSE;
1710: snes->max_its = 50;
1711: snes->max_funcs = 10000;
1712: snes->norm = 0.0;
1713: snes->xnorm = 0.0;
1714: snes->ynorm = 0.0;
1715: snes->normschedule = SNES_NORM_ALWAYS;
1716: snes->functype = SNES_FUNCTION_DEFAULT;
1717: #if defined(PETSC_USE_REAL_SINGLE)
1718: snes->rtol = 1.e-5;
1719: #else
1720: snes->rtol = 1.e-8;
1721: #endif
1722: snes->ttol = 0.0;
1723: #if defined(PETSC_USE_REAL_SINGLE)
1724: snes->abstol = 1.e-25;
1725: #else
1726: snes->abstol = 1.e-50;
1727: #endif
1728: #if defined(PETSC_USE_REAL_SINGLE)
1729: snes->stol = 1.e-5;
1730: #else
1731: snes->stol = 1.e-8;
1732: #endif
1733: #if defined(PETSC_USE_REAL_SINGLE)
1734: snes->deltatol = 1.e-6;
1735: #else
1736: snes->deltatol = 1.e-12;
1737: #endif
1738: snes->divtol = 1.e4;
1739: snes->rnorm0 = 0;
1740: snes->nfuncs = 0;
1741: snes->numFailures = 0;
1742: snes->maxFailures = 1;
1743: snes->linear_its = 0;
1744: snes->lagjacobian = 1;
1745: snes->jac_iter = 0;
1746: snes->lagjac_persist = PETSC_FALSE;
1747: snes->lagpreconditioner = 1;
1748: snes->pre_iter = 0;
1749: snes->lagpre_persist = PETSC_FALSE;
1750: snes->numbermonitors = 0;
1751: snes->data = 0;
1752: snes->setupcalled = PETSC_FALSE;
1753: snes->ksp_ewconv = PETSC_FALSE;
1754: snes->nwork = 0;
1755: snes->work = 0;
1756: snes->nvwork = 0;
1757: snes->vwork = 0;
1758: snes->conv_hist_len = 0;
1759: snes->conv_hist_max = 0;
1760: snes->conv_hist = NULL;
1761: snes->conv_hist_its = NULL;
1762: snes->conv_hist_reset = PETSC_TRUE;
1763: snes->counters_reset = PETSC_TRUE;
1764: snes->vec_func_init_set = PETSC_FALSE;
1765: snes->reason = SNES_CONVERGED_ITERATING;
1766: snes->npcside = PC_RIGHT;
1767: snes->setfromoptionscalled = 0;
1769: snes->mf = PETSC_FALSE;
1770: snes->mf_operator = PETSC_FALSE;
1771: snes->mf_version = 1;
1773: snes->numLinearSolveFailures = 0;
1774: snes->maxLinearSolveFailures = 1;
1776: snes->vizerotolerance = 1.e-8;
1777: #if defined(PETSC_USE_DEBUG)
1778: snes->checkjacdomainerror = PETSC_TRUE;
1779: #else
1780: snes->checkjacdomainerror = PETSC_FALSE;
1781: #endif
1783: /* Set this to true if the implementation of SNESSolve_XXX does compute the residual at the final solution. */
1784: snes->alwayscomputesfinalresidual = PETSC_FALSE;
1786: /* Create context to compute Eisenstat-Walker relative tolerance for KSP */
1787: PetscNewLog(snes,&kctx);
1789: snes->kspconvctx = (void*)kctx;
1790: kctx->version = 2;
1791: kctx->rtol_0 = .3; /* Eisenstat and Walker suggest rtol_0=.5, but
1792: this was too large for some test cases */
1793: kctx->rtol_last = 0.0;
1794: kctx->rtol_max = .9;
1795: kctx->gamma = 1.0;
1796: kctx->alpha = .5*(1.0 + PetscSqrtReal(5.0));
1797: kctx->alpha2 = kctx->alpha;
1798: kctx->threshold = .1;
1799: kctx->lresid_last = 0.0;
1800: kctx->norm_last = 0.0;
1802: *outsnes = snes;
1803: return(0);
1804: }
1806: /*MC
1807: SNESFunction - Functional form used to convey the nonlinear function to be solved by SNES
1809: Synopsis:
1810: #include "petscsnes.h"
1811: PetscErrorCode SNESFunction(SNES snes,Vec x,Vec f,void *ctx);
1813: Input Parameters:
1814: + snes - the SNES context
1815: . x - state at which to evaluate residual
1816: - ctx - optional user-defined function context, passed in with SNESSetFunction()
1818: Output Parameter:
1819: . f - vector to put residual (function value)
1821: Level: intermediate
1823: .seealso: SNESSetFunction(), SNESGetFunction()
1824: M*/
1826: /*@C
1827: SNESSetFunction - Sets the function evaluation routine and function
1828: vector for use by the SNES routines in solving systems of nonlinear
1829: equations.
1831: Logically Collective on SNES
1833: Input Parameters:
1834: + snes - the SNES context
1835: . r - vector to store function value
1836: . f - function evaluation routine; see SNESFunction for calling sequence details
1837: - ctx - [optional] user-defined context for private data for the
1838: function evaluation routine (may be NULL)
1840: Notes:
1841: The Newton-like methods typically solve linear systems of the form
1842: $ f'(x) x = -f(x),
1843: where f'(x) denotes the Jacobian matrix and f(x) is the function.
1845: Level: beginner
1847: .keywords: SNES, nonlinear, set, function
1849: .seealso: SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESSetPicard(), SNESFunction
1850: @*/
1851: PetscErrorCode SNESSetFunction(SNES snes,Vec r,PetscErrorCode (*f)(SNES,Vec,Vec,void*),void *ctx)
1852: {
1854: DM dm;
1858: if (r) {
1861: PetscObjectReference((PetscObject)r);
1862: VecDestroy(&snes->vec_func);
1864: snes->vec_func = r;
1865: }
1866: SNESGetDM(snes,&dm);
1867: DMSNESSetFunction(dm,f,ctx);
1868: return(0);
1869: }
1872: /*@C
1873: SNESSetInitialFunction - Sets the function vector to be used as the
1874: function norm at the initialization of the method. In some
1875: instances, the user has precomputed the function before calling
1876: SNESSolve. This function allows one to avoid a redundant call
1877: to SNESComputeFunction in that case.
1879: Logically Collective on SNES
1881: Input Parameters:
1882: + snes - the SNES context
1883: - f - vector to store function value
1885: Notes:
1886: This should not be modified during the solution procedure.
1888: This is used extensively in the SNESFAS hierarchy and in nonlinear preconditioning.
1890: Level: developer
1892: .keywords: SNES, nonlinear, set, function
1894: .seealso: SNESSetFunction(), SNESComputeFunction(), SNESSetInitialFunctionNorm()
1895: @*/
1896: PetscErrorCode SNESSetInitialFunction(SNES snes, Vec f)
1897: {
1899: Vec vec_func;
1905: if (snes->npcside== PC_LEFT && snes->functype == SNES_FUNCTION_PRECONDITIONED) {
1906: snes->vec_func_init_set = PETSC_FALSE;
1907: return(0);
1908: }
1909: SNESGetFunction(snes,&vec_func,NULL,NULL);
1910: VecCopy(f, vec_func);
1912: snes->vec_func_init_set = PETSC_TRUE;
1913: return(0);
1914: }
1916: /*@
1917: SNESSetNormSchedule - Sets the SNESNormSchedule used in covergence and monitoring
1918: of the SNES method.
1920: Logically Collective on SNES
1922: Input Parameters:
1923: + snes - the SNES context
1924: - normschedule - the frequency of norm computation
1926: Options Database Key:
1927: . -snes_norm_schedule <none, always, initialonly, finalonly, initalfinalonly>
1929: Notes:
1930: Only certain SNES methods support certain SNESNormSchedules. Most require evaluation
1931: of the nonlinear function and the taking of its norm at every iteration to
1932: even ensure convergence at all. However, methods such as custom Gauss-Seidel methods
1933: (SNESNGS) and the like do not require the norm of the function to be computed, and therfore
1934: may either be monitored for convergence or not. As these are often used as nonlinear
1935: preconditioners, monitoring the norm of their error is not a useful enterprise within
1936: their solution.
1938: Level: developer
1940: .keywords: SNES, nonlinear, set, function, norm, type
1942: .seealso: SNESGetNormSchedule(), SNESComputeFunction(), VecNorm(), SNESSetFunction(), SNESSetInitialFunction(), SNESNormSchedule
1943: @*/
1944: PetscErrorCode SNESSetNormSchedule(SNES snes, SNESNormSchedule normschedule)
1945: {
1948: snes->normschedule = normschedule;
1949: return(0);
1950: }
1953: /*@
1954: SNESGetNormSchedule - Gets the SNESNormSchedule used in covergence and monitoring
1955: of the SNES method.
1957: Logically Collective on SNES
1959: Input Parameters:
1960: + snes - the SNES context
1961: - normschedule - the type of the norm used
1963: Level: advanced
1965: .keywords: SNES, nonlinear, set, function, norm, type
1967: .seealso: SNESSetNormSchedule(), SNESComputeFunction(), VecNorm(), SNESSetFunction(), SNESSetInitialFunction(), SNESNormSchedule
1968: @*/
1969: PetscErrorCode SNESGetNormSchedule(SNES snes, SNESNormSchedule *normschedule)
1970: {
1973: *normschedule = snes->normschedule;
1974: return(0);
1975: }
1978: /*@
1979: SNESSetFunctionNorm - Sets the last computed residual norm.
1981: Logically Collective on SNES
1983: Input Parameters:
1984: + snes - the SNES context
1986: - normschedule - the frequency of norm computation
1988: Level: developer
1990: .keywords: SNES, nonlinear, set, function, norm, type
1991: .seealso: SNESGetNormSchedule(), SNESComputeFunction(), VecNorm(), SNESSetFunction(), SNESSetInitialFunction(), SNESNormSchedule
1992: @*/
1993: PetscErrorCode SNESSetFunctionNorm(SNES snes, PetscReal norm)
1994: {
1997: snes->norm = norm;
1998: return(0);
1999: }
2001: /*@
2002: SNESGetFunctionNorm - Gets the last computed norm of the residual
2004: Not Collective
2006: Input Parameter:
2007: . snes - the SNES context
2009: Output Parameter:
2010: . norm - the last computed residual norm
2012: Level: developer
2014: .keywords: SNES, nonlinear, set, function, norm, type
2015: .seealso: SNESSetNormSchedule(), SNESComputeFunction(), VecNorm(), SNESSetFunction(), SNESSetInitialFunction(), SNESNormSchedule
2016: @*/
2017: PetscErrorCode SNESGetFunctionNorm(SNES snes, PetscReal *norm)
2018: {
2022: *norm = snes->norm;
2023: return(0);
2024: }
2026: /*@
2027: SNESGetUpdateNorm - Gets the last computed norm of the Newton update
2029: Not Collective
2031: Input Parameter:
2032: . snes - the SNES context
2034: Output Parameter:
2035: . ynorm - the last computed update norm
2037: Level: developer
2039: .keywords: SNES, nonlinear, set, function, norm, type
2040: .seealso: SNESSetNormSchedule(), SNESComputeFunction(), SNESGetFunctionNorm()
2041: @*/
2042: PetscErrorCode SNESGetUpdateNorm(SNES snes, PetscReal *ynorm)
2043: {
2047: *ynorm = snes->ynorm;
2048: return(0);
2049: }
2051: /*@
2052: SNESGetSolutionNorm - Gets the last computed norm of the solution
2054: Not Collective
2056: Input Parameter:
2057: . snes - the SNES context
2059: Output Parameter:
2060: . xnorm - the last computed solution norm
2062: Level: developer
2064: .keywords: SNES, nonlinear, set, function, norm, type
2065: .seealso: SNESSetNormSchedule(), SNESComputeFunction(), SNESGetFunctionNorm(), SNESGetUpdateNorm()
2066: @*/
2067: PetscErrorCode SNESGetSolutionNorm(SNES snes, PetscReal *xnorm)
2068: {
2072: *xnorm = snes->xnorm;
2073: return(0);
2074: }
2076: /*@C
2077: SNESSetFunctionType - Sets the SNESNormSchedule used in covergence and monitoring
2078: of the SNES method.
2080: Logically Collective on SNES
2082: Input Parameters:
2083: + snes - the SNES context
2084: - normschedule - the frequency of norm computation
2086: Notes:
2087: Only certain SNES methods support certain SNESNormSchedules. Most require evaluation
2088: of the nonlinear function and the taking of its norm at every iteration to
2089: even ensure convergence at all. However, methods such as custom Gauss-Seidel methods
2090: (SNESNGS) and the like do not require the norm of the function to be computed, and therfore
2091: may either be monitored for convergence or not. As these are often used as nonlinear
2092: preconditioners, monitoring the norm of their error is not a useful enterprise within
2093: their solution.
2095: Level: developer
2097: .keywords: SNES, nonlinear, set, function, norm, type
2099: .seealso: SNESGetNormSchedule(), SNESComputeFunction(), VecNorm(), SNESSetFunction(), SNESSetInitialFunction(), SNESNormSchedule
2100: @*/
2101: PetscErrorCode SNESSetFunctionType(SNES snes, SNESFunctionType type)
2102: {
2105: snes->functype = type;
2106: return(0);
2107: }
2110: /*@C
2111: SNESGetFunctionType - Gets the SNESNormSchedule used in covergence and monitoring
2112: of the SNES method.
2114: Logically Collective on SNES
2116: Input Parameters:
2117: + snes - the SNES context
2118: - normschedule - the type of the norm used
2120: Level: advanced
2122: .keywords: SNES, nonlinear, set, function, norm, type
2124: .seealso: SNESSetNormSchedule(), SNESComputeFunction(), VecNorm(), SNESSetFunction(), SNESSetInitialFunction(), SNESNormSchedule
2125: @*/
2126: PetscErrorCode SNESGetFunctionType(SNES snes, SNESFunctionType *type)
2127: {
2130: *type = snes->functype;
2131: return(0);
2132: }
2134: /*MC
2135: SNESNGSFunction - function used to convey a Gauss-Seidel sweep on the nonlinear function
2137: Synopsis:
2138: #include <petscsnes.h>
2139: $ SNESNGSFunction(SNES snes,Vec x,Vec b,void *ctx);
2141: + X - solution vector
2142: . B - RHS vector
2143: - ctx - optional user-defined Gauss-Seidel context
2145: Level: intermediate
2147: .seealso: SNESSetNGS(), SNESGetNGS()
2148: M*/
2150: /*@C
2151: SNESSetNGS - Sets the user nonlinear Gauss-Seidel routine for
2152: use with composed nonlinear solvers.
2154: Input Parameters:
2155: + snes - the SNES context
2156: . f - function evaluation routine to apply Gauss-Seidel see SNESNGSFunction
2157: - ctx - [optional] user-defined context for private data for the
2158: smoother evaluation routine (may be NULL)
2160: Notes:
2161: The NGS routines are used by the composed nonlinear solver to generate
2162: a problem appropriate update to the solution, particularly FAS.
2164: Level: intermediate
2166: .keywords: SNES, nonlinear, set, Gauss-Seidel
2168: .seealso: SNESGetFunction(), SNESComputeNGS()
2169: @*/
2170: PetscErrorCode SNESSetNGS(SNES snes,PetscErrorCode (*f)(SNES,Vec,Vec,void*),void *ctx)
2171: {
2173: DM dm;
2177: SNESGetDM(snes,&dm);
2178: DMSNESSetNGS(dm,f,ctx);
2179: return(0);
2180: }
2182: PetscErrorCode SNESPicardComputeFunction(SNES snes,Vec x,Vec f,void *ctx)
2183: {
2185: DM dm;
2186: DMSNES sdm;
2189: SNESGetDM(snes,&dm);
2190: DMGetDMSNES(dm,&sdm);
2191: if (!sdm->ops->computepfunction) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE, "Must call SNESSetPicard() to provide Picard function.");
2192: if (!sdm->ops->computepjacobian) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE, "Must call SNESSetPicard() to provide Picard Jacobian.");
2193: /* A(x)*x - b(x) */
2194: PetscStackPush("SNES Picard user function");
2195: (*sdm->ops->computepfunction)(snes,x,f,sdm->pctx);
2196: PetscStackPop;
2197: PetscStackPush("SNES Picard user Jacobian");
2198: (*sdm->ops->computepjacobian)(snes,x,snes->jacobian,snes->jacobian_pre,sdm->pctx);
2199: PetscStackPop;
2200: VecScale(f,-1.0);
2201: MatMultAdd(snes->jacobian,x,f,f);
2202: return(0);
2203: }
2205: PetscErrorCode SNESPicardComputeJacobian(SNES snes,Vec x1,Mat J,Mat B,void *ctx)
2206: {
2208: /* the jacobian matrix should be pre-filled in SNESPicardComputeFunction */
2209: return(0);
2210: }
2212: /*@C
2213: SNESSetPicard - Use SNES to solve the semilinear-system A(x) x = b(x) via a Picard type iteration (Picard linearization)
2215: Logically Collective on SNES
2217: Input Parameters:
2218: + snes - the SNES context
2219: . r - vector to store function value
2220: . b - function evaluation routine
2221: . Amat - matrix with which A(x) x - b(x) is to be computed
2222: . Pmat - matrix from which preconditioner is computed (usually the same as Amat)
2223: . J - function to compute matrix value, see SNESJacobianFunction for details on its calling sequence
2224: - ctx - [optional] user-defined context for private data for the
2225: function evaluation routine (may be NULL)
2227: Notes:
2228: 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
2229: 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.
2231: One can call SNESSetPicard() or SNESSetFunction() (and possibly SNESSetJacobian()) but cannot call both
2233: $ 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}
2234: $ Note that when an exact solver is used this corresponds to the "classic" Picard A(x^{n}) x^{n+1} = b(x^{n}) iteration.
2236: Run with -snes_mf_operator to solve the system with Newton's method using A(x^{n}) to construct the preconditioner.
2238: We implement the defect correction form of the Picard iteration because it converges much more generally when inexact linear solvers are used then
2239: the direct Picard iteration A(x^n) x^{n+1} = b(x^n)
2241: 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
2242: 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
2243: different please contact us at petsc-dev@mcs.anl.gov and we'll have an entirely new argument :-).
2245: Level: intermediate
2247: .keywords: SNES, nonlinear, set, function
2249: .seealso: SNESGetFunction(), SNESSetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESGetPicard(), SNESLineSearchPreCheckPicard(), SNESJacobianFunction
2250: @*/
2251: 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)
2252: {
2254: DM dm;
2258: SNESGetDM(snes, &dm);
2259: DMSNESSetPicard(dm,b,J,ctx);
2260: SNESSetFunction(snes,r,SNESPicardComputeFunction,ctx);
2261: SNESSetJacobian(snes,Amat,Pmat,SNESPicardComputeJacobian,ctx);
2262: return(0);
2263: }
2265: /*@C
2266: SNESGetPicard - Returns the context for the Picard iteration
2268: Not Collective, but Vec is parallel if SNES is parallel. Collective if Vec is requested, but has not been created yet.
2270: Input Parameter:
2271: . snes - the SNES context
2273: Output Parameter:
2274: + r - the function (or NULL)
2275: . f - the function (or NULL); see SNESFunction for calling sequence details
2276: . Amat - the matrix used to defined the operation A(x) x - b(x) (or NULL)
2277: . Pmat - the matrix from which the preconditioner will be constructed (or NULL)
2278: . J - the function for matrix evaluation (or NULL); see SNESJacobianFunction for calling sequence details
2279: - ctx - the function context (or NULL)
2281: Level: advanced
2283: .keywords: SNES, nonlinear, get, function
2285: .seealso: SNESSetPicard(), SNESGetFunction(), SNESGetJacobian(), SNESGetDM(), SNESFunction, SNESJacobianFunction
2286: @*/
2287: 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)
2288: {
2290: DM dm;
2294: SNESGetFunction(snes,r,NULL,NULL);
2295: SNESGetJacobian(snes,Amat,Pmat,NULL,NULL);
2296: SNESGetDM(snes,&dm);
2297: DMSNESGetPicard(dm,f,J,ctx);
2298: return(0);
2299: }
2301: /*@C
2302: SNESSetComputeInitialGuess - Sets a routine used to compute an initial guess for the problem
2304: Logically Collective on SNES
2306: Input Parameters:
2307: + snes - the SNES context
2308: . func - function evaluation routine
2309: - ctx - [optional] user-defined context for private data for the
2310: function evaluation routine (may be NULL)
2312: Calling sequence of func:
2313: $ func (SNES snes,Vec x,void *ctx);
2315: . f - function vector
2316: - ctx - optional user-defined function context
2318: Level: intermediate
2320: .keywords: SNES, nonlinear, set, function
2322: .seealso: SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian()
2323: @*/
2324: PetscErrorCode SNESSetComputeInitialGuess(SNES snes,PetscErrorCode (*func)(SNES,Vec,void*),void *ctx)
2325: {
2328: if (func) snes->ops->computeinitialguess = func;
2329: if (ctx) snes->initialguessP = ctx;
2330: return(0);
2331: }
2333: /* --------------------------------------------------------------- */
2334: /*@C
2335: SNESGetRhs - Gets the vector for solving F(x) = rhs. If rhs is not set
2336: it assumes a zero right hand side.
2338: Logically Collective on SNES
2340: Input Parameter:
2341: . snes - the SNES context
2343: Output Parameter:
2344: . rhs - the right hand side vector or NULL if the right hand side vector is null
2346: Level: intermediate
2348: .keywords: SNES, nonlinear, get, function, right hand side
2350: .seealso: SNESGetSolution(), SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESSetFunction()
2351: @*/
2352: PetscErrorCode SNESGetRhs(SNES snes,Vec *rhs)
2353: {
2357: *rhs = snes->vec_rhs;
2358: return(0);
2359: }
2361: /*@
2362: SNESComputeFunction - Calls the function that has been set with SNESSetFunction().
2364: Collective on SNES
2366: Input Parameters:
2367: + snes - the SNES context
2368: - x - input vector
2370: Output Parameter:
2371: . y - function vector, as set by SNESSetFunction()
2373: Notes:
2374: SNESComputeFunction() is typically used within nonlinear solvers
2375: implementations, so most users would not generally call this routine
2376: themselves.
2378: Level: developer
2380: .keywords: SNES, nonlinear, compute, function
2382: .seealso: SNESSetFunction(), SNESGetFunction()
2383: @*/
2384: PetscErrorCode SNESComputeFunction(SNES snes,Vec x,Vec y)
2385: {
2387: DM dm;
2388: DMSNES sdm;
2396: VecValidValues(x,2,PETSC_TRUE);
2398: SNESGetDM(snes,&dm);
2399: DMGetDMSNES(dm,&sdm);
2400: if (sdm->ops->computefunction) {
2401: if (sdm->ops->computefunction != SNESObjectiveComputeFunctionDefaultFD) {
2402: PetscLogEventBegin(SNES_FunctionEval,snes,x,y,0);
2403: }
2404: VecLockReadPush(x);
2405: PetscStackPush("SNES user function");
2406: (*sdm->ops->computefunction)(snes,x,y,sdm->functionctx);
2407: PetscStackPop;
2408: VecLockReadPop(x);
2409: if (sdm->ops->computefunction != SNESObjectiveComputeFunctionDefaultFD) {
2410: PetscLogEventEnd(SNES_FunctionEval,snes,x,y,0);
2411: }
2412: } else if (snes->vec_rhs) {
2413: MatMult(snes->jacobian, x, y);
2414: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE, "Must call SNESSetFunction() or SNESSetDM() before SNESComputeFunction(), likely called from SNESSolve().");
2415: if (snes->vec_rhs) {
2416: VecAXPY(y,-1.0,snes->vec_rhs);
2417: }
2418: snes->nfuncs++;
2419: /*
2420: domainerror might not be set on all processes; so we tag vector locally with Inf and the next inner product or norm will
2421: propagate the value to all processes
2422: */
2423: if (snes->domainerror) {
2424: VecSetInf(y);
2425: }
2426: return(0);
2427: }
2429: /*@
2430: SNESComputeNGS - Calls the Gauss-Seidel function that has been set with SNESSetNGS().
2432: Collective on SNES
2434: Input Parameters:
2435: + snes - the SNES context
2436: . x - input vector
2437: - b - rhs vector
2439: Output Parameter:
2440: . x - new solution vector
2442: Notes:
2443: SNESComputeNGS() is typically used within composed nonlinear solver
2444: implementations, so most users would not generally call this routine
2445: themselves.
2447: Level: developer
2449: .keywords: SNES, nonlinear, compute, function
2451: .seealso: SNESSetNGS(), SNESComputeFunction()
2452: @*/
2453: PetscErrorCode SNESComputeNGS(SNES snes,Vec b,Vec x)
2454: {
2456: DM dm;
2457: DMSNES sdm;
2465: if (b) {VecValidValues(b,2,PETSC_TRUE);}
2466: PetscLogEventBegin(SNES_NGSEval,snes,x,b,0);
2467: SNESGetDM(snes,&dm);
2468: DMGetDMSNES(dm,&sdm);
2469: if (sdm->ops->computegs) {
2470: if (b) {VecLockReadPush(b);}
2471: PetscStackPush("SNES user NGS");
2472: (*sdm->ops->computegs)(snes,x,b,sdm->gsctx);
2473: PetscStackPop;
2474: if (b) {VecLockReadPop(b);}
2475: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE, "Must call SNESSetNGS() before SNESComputeNGS(), likely called from SNESSolve().");
2476: PetscLogEventEnd(SNES_NGSEval,snes,x,b,0);
2477: return(0);
2478: }
2480: PetscErrorCode SNESTestJacobian(SNES snes)
2481: {
2482: Mat A,B,C,D,jacobian;
2483: Vec x = snes->vec_sol,f = snes->vec_func;
2484: PetscErrorCode ierr;
2485: PetscReal nrm,gnorm;
2486: PetscReal threshold = 1.e-5;
2487: PetscInt m,n,M,N;
2488: void *functx;
2489: PetscBool complete_print = PETSC_FALSE,threshold_print = PETSC_FALSE,test = PETSC_FALSE,flg;
2490: PetscViewer viewer,mviewer;
2491: MPI_Comm comm;
2492: PetscInt tabs;
2493: static PetscBool directionsprinted = PETSC_FALSE;
2494: PetscViewerFormat format;
2497: PetscObjectOptionsBegin((PetscObject)snes);
2498: PetscOptionsName("-snes_test_jacobian","Compare hand-coded and finite difference Jacobians","None",&test);
2499: PetscOptionsReal("-snes_test_jacobian", "Threshold for element difference between hand-coded and finite difference being meaningful", "None", threshold, &threshold,NULL);
2500: PetscOptionsViewer("-snes_test_jacobian_view","View difference between hand-coded and finite difference Jacobians element entries","None",&mviewer,&format,&complete_print);
2501: if (!complete_print) {
2502: PetscOptionsViewer("-snes_test_jacobian_display","Display difference between hand-coded and finite difference Jacobians","None",&mviewer,&format,&complete_print);
2503: }
2504: /* for compatibility with PETSc 3.9 and older. */
2505: PetscOptionsReal("-snes_test_jacobian_display_threshold", "Display difference between hand-coded and finite difference Jacobians which exceed input threshold", "None", threshold, &threshold, &threshold_print);
2506: PetscOptionsEnd();
2507: if (!test) return(0);
2509: PetscObjectGetComm((PetscObject)snes,&comm);
2510: PetscViewerASCIIGetStdout(comm,&viewer);
2511: PetscViewerASCIIGetTab(viewer, &tabs);
2512: PetscViewerASCIISetTab(viewer, ((PetscObject)snes)->tablevel);
2513: PetscViewerASCIIPrintf(viewer," ---------- Testing Jacobian -------------\n");
2514: if (!complete_print && !directionsprinted) {
2515: PetscViewerASCIIPrintf(viewer," Run with -snes_test_jacobian_view and optionally -snes_test_jacobian <threshold> to show difference\n");
2516: PetscViewerASCIIPrintf(viewer," of hand-coded and finite difference Jacobian entries greater than <threshold>.\n");
2517: }
2518: if (!directionsprinted) {
2519: PetscViewerASCIIPrintf(viewer," Testing hand-coded Jacobian, if (for double precision runs) ||J - Jfd||_F/||J||_F is\n");
2520: PetscViewerASCIIPrintf(viewer," O(1.e-8), the hand-coded Jacobian is probably correct.\n");
2521: directionsprinted = PETSC_TRUE;
2522: }
2523: if (complete_print) {
2524: PetscViewerPushFormat(mviewer,format);
2525: }
2527: /* evaluate the function at this point because SNESComputeJacobianDefault() assumes that the function has been evaluated and put into snes->vec_func */
2528: SNESComputeFunction(snes,x,f);
2530: PetscObjectTypeCompare((PetscObject)snes->jacobian,MATMFFD,&flg);
2531: if (!flg) jacobian = snes->jacobian;
2532: else jacobian = snes->jacobian_pre;
2534: while (jacobian) {
2535: PetscObjectTypeCompareAny((PetscObject)jacobian,&flg,MATSEQAIJ,MATMPIAIJ,MATSEQDENSE,MATMPIDENSE,MATSEQBAIJ,MATMPIBAIJ,MATSEQSBAIJ,MATMPIBAIJ,"");
2536: if (flg) {
2537: A = jacobian;
2538: PetscObjectReference((PetscObject)A);
2539: } else {
2540: MatComputeExplicitOperator(jacobian,&A);
2541: }
2543: MatCreate(PetscObjectComm((PetscObject)A),&B);
2544: MatGetSize(A,&M,&N);
2545: MatGetLocalSize(A,&m,&n);
2546: MatSetSizes(B,m,n,M,N);
2547: MatSetType(B,((PetscObject)A)->type_name);
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(jacobian,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: MatAYPX(B,-1.0,A,DIFFERENT_NONZERO_PATTERN);
2576: MatCreate(PetscObjectComm((PetscObject)A),&C);
2577: MatSetSizes(C,m,n,M,N);
2578: MatSetType(C,((PetscObject)A)->type_name);
2579: MatSetUp(C);
2580: MatSetOption(C,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);
2581: MatGetOwnershipRange(B,&Istart,&Iend);
2583: for (row = Istart; row < Iend; row++) {
2584: MatGetRow(B,row,&bncols,&bcols,&bvals);
2585: PetscMalloc2(bncols,&ccols,bncols,&cvals);
2586: for (j = 0, cncols = 0; j < bncols; j++) {
2587: if (PetscAbsScalar(bvals[j]) > threshold) {
2588: ccols[cncols] = bcols[j];
2589: cvals[cncols] = bvals[j];
2590: cncols += 1;
2591: }
2592: }
2593: if (cncols) {
2594: MatSetValues(C,1,&row,cncols,ccols,cvals,INSERT_VALUES);
2595: }
2596: MatRestoreRow(B,row,&bncols,&bcols,&bvals);
2597: PetscFree2(ccols,cvals);
2598: }
2599: MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
2600: MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
2601: PetscViewerASCIIPrintf(viewer," Hand-coded minus finite-difference Jacobian with tolerance %g ----------\n",(double)threshold);
2602: MatView(C,complete_print ? mviewer : viewer);
2603: MatDestroy(&C);
2604: }
2605: MatDestroy(&A);
2606: MatDestroy(&B);
2608: if (jacobian != snes->jacobian_pre) {
2609: jacobian = snes->jacobian_pre;
2610: PetscViewerASCIIPrintf(viewer," ---------- Testing Jacobian for preconditioner -------------\n");
2611: }
2612: else jacobian = NULL;
2613: }
2614: if (complete_print) {
2615: PetscViewerPopFormat(mviewer);
2616: }
2617: if (mviewer) { PetscViewerDestroy(&mviewer); }
2618: PetscViewerASCIISetTab(viewer,tabs);
2619: return(0);
2620: }
2622: /*@
2623: SNESComputeJacobian - Computes the Jacobian matrix that has been set with SNESSetJacobian().
2625: Collective on SNES and Mat
2627: Input Parameters:
2628: + snes - the SNES context
2629: - x - input vector
2631: Output Parameters:
2632: + A - Jacobian matrix
2633: - B - optional preconditioning matrix
2635: Options Database Keys:
2636: + -snes_lag_preconditioner <lag>
2637: . -snes_lag_jacobian <lag>
2638: . -snes_test_jacobian - compare the user provided Jacobian with one compute via finite differences to check for errors
2639: . -snes_test_jacobian_display - 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
2640: . -snes_test_jacobian_display_threshold <numerical value> - display entries in the difference between the user provided Jacobian and finite difference Jacobian that are greater than a certain value to help users detect errors
2641: . -snes_compare_explicit - Compare the computed Jacobian to the finite difference Jacobian and output the differences
2642: . -snes_compare_explicit_draw - Compare the computed Jacobian to the finite difference Jacobian and draw the result
2643: . -snes_compare_explicit_contour - Compare the computed Jacobian to the finite difference Jacobian and draw a contour plot with the result
2644: . -snes_compare_operator - Make the comparison options above use the operator instead of the preconditioning matrix
2645: . -snes_compare_coloring - Compute the finite difference Jacobian using coloring and display norms of difference
2646: . -snes_compare_coloring_display - Compute the finite differece Jacobian using coloring and display verbose differences
2647: . -snes_compare_coloring_threshold - Display only those matrix entries that differ by more than a given threshold
2648: . -snes_compare_coloring_threshold_atol - Absolute tolerance for difference in matrix entries to be displayed by -snes_compare_coloring_threshold
2649: . -snes_compare_coloring_threshold_rtol - Relative tolerance for difference in matrix entries to be displayed by -snes_compare_coloring_threshold
2650: . -snes_compare_coloring_draw - Compute the finite differece Jacobian using coloring and draw differences
2651: - -snes_compare_coloring_draw_contour - Compute the finite differece Jacobian using coloring and show contours of matrices and differences
2654: Notes:
2655: Most users should not need to explicitly call this routine, as it
2656: is used internally within the nonlinear solvers.
2658: Developer Notes:
2659: 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
2660: for with the SNESType of test that has been removed.
2662: Level: developer
2664: .keywords: SNES, compute, Jacobian, matrix
2666: .seealso: SNESSetJacobian(), KSPSetOperators(), MatStructure, SNESSetLagPreconditioner(), SNESSetLagJacobian()
2667: @*/
2668: PetscErrorCode SNESComputeJacobian(SNES snes,Vec X,Mat A,Mat B)
2669: {
2671: PetscBool flag;
2672: DM dm;
2673: DMSNES sdm;
2674: KSP ksp;
2680: VecValidValues(X,2,PETSC_TRUE);
2681: SNESGetDM(snes,&dm);
2682: DMGetDMSNES(dm,&sdm);
2684: if (!sdm->ops->computejacobian) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_USER,"Must call SNESSetJacobian(), DMSNESSetJacobian(), DMDASNESSetJacobianLocal(), etc");
2686: /* make sure that MatAssemblyBegin/End() is called on A matrix if it is matrix free */
2688: if (snes->lagjacobian == -2) {
2689: snes->lagjacobian = -1;
2691: PetscInfo(snes,"Recomputing Jacobian/preconditioner because lag is -2 (means compute Jacobian, but then never again) \n");
2692: } else if (snes->lagjacobian == -1) {
2693: PetscInfo(snes,"Reusing Jacobian/preconditioner because lag is -1\n");
2694: PetscObjectTypeCompare((PetscObject)A,MATMFFD,&flag);
2695: if (flag) {
2696: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
2697: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
2698: }
2699: return(0);
2700: } else if (snes->lagjacobian > 1 && (snes->iter + snes->jac_iter) % snes->lagjacobian) {
2701: PetscInfo2(snes,"Reusing Jacobian/preconditioner because lag is %D and SNES iteration is %D\n",snes->lagjacobian,snes->iter);
2702: PetscObjectTypeCompare((PetscObject)A,MATMFFD,&flag);
2703: if (flag) {
2704: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
2705: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
2706: }
2707: return(0);
2708: }
2709: if (snes->npc && snes->npcside== PC_LEFT) {
2710: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
2711: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
2712: return(0);
2713: }
2715: PetscLogEventBegin(SNES_JacobianEval,snes,X,A,B);
2716: VecLockReadPush(X);
2717: PetscStackPush("SNES user Jacobian function");
2718: (*sdm->ops->computejacobian)(snes,X,A,B,sdm->jacobianctx);
2719: PetscStackPop;
2720: VecLockReadPop(X);
2721: PetscLogEventEnd(SNES_JacobianEval,snes,X,A,B);
2723: /* the next line ensures that snes->ksp exists */
2724: SNESGetKSP(snes,&ksp);
2725: if (snes->lagpreconditioner == -2) {
2726: PetscInfo(snes,"Rebuilding preconditioner exactly once since lag is -2\n");
2727: KSPSetReusePreconditioner(snes->ksp,PETSC_FALSE);
2728: snes->lagpreconditioner = -1;
2729: } else if (snes->lagpreconditioner == -1) {
2730: PetscInfo(snes,"Reusing preconditioner because lag is -1\n");
2731: KSPSetReusePreconditioner(snes->ksp,PETSC_TRUE);
2732: } else if (snes->lagpreconditioner > 1 && (snes->iter + snes->pre_iter) % snes->lagpreconditioner) {
2733: PetscInfo2(snes,"Reusing preconditioner because lag is %D and SNES iteration is %D\n",snes->lagpreconditioner,snes->iter);
2734: KSPSetReusePreconditioner(snes->ksp,PETSC_TRUE);
2735: } else {
2736: PetscInfo(snes,"Rebuilding preconditioner\n");
2737: KSPSetReusePreconditioner(snes->ksp,PETSC_FALSE);
2738: }
2740: SNESTestJacobian(snes);
2741: /* make sure user returned a correct Jacobian and preconditioner */
2744: {
2745: PetscBool flag = PETSC_FALSE,flag_draw = PETSC_FALSE,flag_contour = PETSC_FALSE,flag_operator = PETSC_FALSE;
2746: PetscOptionsGetViewer(PetscObjectComm((PetscObject)snes),((PetscObject) snes)->options,((PetscObject)snes)->prefix,"-snes_compare_explicit",NULL,NULL,&flag);
2747: PetscOptionsGetViewer(PetscObjectComm((PetscObject)snes),((PetscObject) snes)->options,((PetscObject)snes)->prefix,"-snes_compare_explicit_draw",NULL,NULL,&flag_draw);
2748: PetscOptionsGetViewer(PetscObjectComm((PetscObject)snes),((PetscObject) snes)->options,((PetscObject)snes)->prefix,"-snes_compare_explicit_draw_contour",NULL,NULL,&flag_contour);
2749: PetscOptionsGetViewer(PetscObjectComm((PetscObject)snes),((PetscObject) snes)->options,((PetscObject)snes)->prefix,"-snes_compare_operator",NULL,NULL,&flag_operator);
2750: if (flag || flag_draw || flag_contour) {
2751: Mat Bexp_mine = NULL,Bexp,FDexp;
2752: PetscViewer vdraw,vstdout;
2753: PetscBool flg;
2754: if (flag_operator) {
2755: MatComputeExplicitOperator(A,&Bexp_mine);
2756: Bexp = Bexp_mine;
2757: } else {
2758: /* See if the preconditioning matrix can be viewed and added directly */
2759: PetscObjectTypeCompareAny((PetscObject)B,&flg,MATSEQAIJ,MATMPIAIJ,MATSEQDENSE,MATMPIDENSE,MATSEQBAIJ,MATMPIBAIJ,MATSEQSBAIJ,MATMPIBAIJ,"");
2760: if (flg) Bexp = B;
2761: else {
2762: /* If the "preconditioning" matrix is itself MATSHELL or some other type without direct support */
2763: MatComputeExplicitOperator(B,&Bexp_mine);
2764: Bexp = Bexp_mine;
2765: }
2766: }
2767: MatConvert(Bexp,MATSAME,MAT_INITIAL_MATRIX,&FDexp);
2768: SNESComputeJacobianDefault(snes,X,FDexp,FDexp,NULL);
2769: PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)snes),&vstdout);
2770: if (flag_draw || flag_contour) {
2771: PetscViewerDrawOpen(PetscObjectComm((PetscObject)snes),0,"Explicit Jacobians",PETSC_DECIDE,PETSC_DECIDE,300,300,&vdraw);
2772: if (flag_contour) {PetscViewerPushFormat(vdraw,PETSC_VIEWER_DRAW_CONTOUR);}
2773: } else vdraw = NULL;
2774: PetscViewerASCIIPrintf(vstdout,"Explicit %s\n",flag_operator ? "Jacobian" : "preconditioning Jacobian");
2775: if (flag) {MatView(Bexp,vstdout);}
2776: if (vdraw) {MatView(Bexp,vdraw);}
2777: PetscViewerASCIIPrintf(vstdout,"Finite difference Jacobian\n");
2778: if (flag) {MatView(FDexp,vstdout);}
2779: if (vdraw) {MatView(FDexp,vdraw);}
2780: MatAYPX(FDexp,-1.0,Bexp,SAME_NONZERO_PATTERN);
2781: PetscViewerASCIIPrintf(vstdout,"User-provided matrix minus finite difference Jacobian\n");
2782: if (flag) {MatView(FDexp,vstdout);}
2783: if (vdraw) { /* Always use contour for the difference */
2784: PetscViewerPushFormat(vdraw,PETSC_VIEWER_DRAW_CONTOUR);
2785: MatView(FDexp,vdraw);
2786: PetscViewerPopFormat(vdraw);
2787: }
2788: if (flag_contour) {PetscViewerPopFormat(vdraw);}
2789: PetscViewerDestroy(&vdraw);
2790: MatDestroy(&Bexp_mine);
2791: MatDestroy(&FDexp);
2792: }
2793: }
2794: {
2795: PetscBool flag = PETSC_FALSE,flag_display = PETSC_FALSE,flag_draw = PETSC_FALSE,flag_contour = PETSC_FALSE,flag_threshold = PETSC_FALSE;
2796: PetscReal threshold_atol = PETSC_SQRT_MACHINE_EPSILON,threshold_rtol = 10*PETSC_SQRT_MACHINE_EPSILON;
2797: PetscOptionsGetViewer(PetscObjectComm((PetscObject)snes),((PetscObject)snes)->options,((PetscObject)snes)->prefix,"-snes_compare_coloring",NULL,NULL,&flag);
2798: PetscOptionsGetViewer(PetscObjectComm((PetscObject)snes),((PetscObject)snes)->options,((PetscObject)snes)->prefix,"-snes_compare_coloring_display",NULL,NULL,&flag_display);
2799: PetscOptionsGetViewer(PetscObjectComm((PetscObject)snes),((PetscObject)snes)->options,((PetscObject)snes)->prefix,"-snes_compare_coloring_draw",NULL,NULL,&flag_draw);
2800: PetscOptionsGetViewer(PetscObjectComm((PetscObject)snes),((PetscObject)snes)->options,((PetscObject)snes)->prefix,"-snes_compare_coloring_draw_contour",NULL,NULL,&flag_contour);
2801: PetscOptionsGetViewer(PetscObjectComm((PetscObject)snes),((PetscObject)snes)->options,((PetscObject)snes)->prefix,"-snes_compare_coloring_threshold",NULL,NULL,&flag_threshold);
2802: if (flag_threshold) {
2803: PetscOptionsGetReal(((PetscObject)snes)->options,((PetscObject)snes)->prefix,"-snes_compare_coloring_threshold_rtol",&threshold_rtol,NULL);
2804: PetscOptionsGetReal(((PetscObject)snes)->options,((PetscObject)snes)->prefix,"-snes_compare_coloring_threshold_atol",&threshold_atol,NULL);
2805: }
2806: if (flag || flag_display || flag_draw || flag_contour || flag_threshold) {
2807: Mat Bfd;
2808: PetscViewer vdraw,vstdout;
2809: MatColoring coloring;
2810: ISColoring iscoloring;
2811: MatFDColoring matfdcoloring;
2812: PetscErrorCode (*func)(SNES,Vec,Vec,void*);
2813: void *funcctx;
2814: PetscReal norm1,norm2,normmax;
2816: MatDuplicate(B,MAT_DO_NOT_COPY_VALUES,&Bfd);
2817: MatColoringCreate(Bfd,&coloring);
2818: MatColoringSetType(coloring,MATCOLORINGSL);
2819: MatColoringSetFromOptions(coloring);
2820: MatColoringApply(coloring,&iscoloring);
2821: MatColoringDestroy(&coloring);
2822: MatFDColoringCreate(Bfd,iscoloring,&matfdcoloring);
2823: MatFDColoringSetFromOptions(matfdcoloring);
2824: MatFDColoringSetUp(Bfd,iscoloring,matfdcoloring);
2825: ISColoringDestroy(&iscoloring);
2827: /* This method of getting the function is currently unreliable since it doesn't work for DM local functions. */
2828: SNESGetFunction(snes,NULL,&func,&funcctx);
2829: MatFDColoringSetFunction(matfdcoloring,(PetscErrorCode (*)(void))func,funcctx);
2830: PetscObjectSetOptionsPrefix((PetscObject)matfdcoloring,((PetscObject)snes)->prefix);
2831: PetscObjectAppendOptionsPrefix((PetscObject)matfdcoloring,"coloring_");
2832: MatFDColoringSetFromOptions(matfdcoloring);
2833: MatFDColoringApply(Bfd,matfdcoloring,X,snes);
2834: MatFDColoringDestroy(&matfdcoloring);
2836: PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)snes),&vstdout);
2837: if (flag_draw || flag_contour) {
2838: PetscViewerDrawOpen(PetscObjectComm((PetscObject)snes),0,"Colored Jacobians",PETSC_DECIDE,PETSC_DECIDE,300,300,&vdraw);
2839: if (flag_contour) {PetscViewerPushFormat(vdraw,PETSC_VIEWER_DRAW_CONTOUR);}
2840: } else vdraw = NULL;
2841: PetscViewerASCIIPrintf(vstdout,"Explicit preconditioning Jacobian\n");
2842: if (flag_display) {MatView(B,vstdout);}
2843: if (vdraw) {MatView(B,vdraw);}
2844: PetscViewerASCIIPrintf(vstdout,"Colored Finite difference Jacobian\n");
2845: if (flag_display) {MatView(Bfd,vstdout);}
2846: if (vdraw) {MatView(Bfd,vdraw);}
2847: MatAYPX(Bfd,-1.0,B,SAME_NONZERO_PATTERN);
2848: MatNorm(Bfd,NORM_1,&norm1);
2849: MatNorm(Bfd,NORM_FROBENIUS,&norm2);
2850: MatNorm(Bfd,NORM_MAX,&normmax);
2851: PetscViewerASCIIPrintf(vstdout,"User-provided matrix minus finite difference Jacobian, norm1=%g normFrob=%g normmax=%g\n",(double)norm1,(double)norm2,(double)normmax);
2852: if (flag_display) {MatView(Bfd,vstdout);}
2853: if (vdraw) { /* Always use contour for the difference */
2854: PetscViewerPushFormat(vdraw,PETSC_VIEWER_DRAW_CONTOUR);
2855: MatView(Bfd,vdraw);
2856: PetscViewerPopFormat(vdraw);
2857: }
2858: if (flag_contour) {PetscViewerPopFormat(vdraw);}
2860: if (flag_threshold) {
2861: PetscInt bs,rstart,rend,i;
2862: MatGetBlockSize(B,&bs);
2863: MatGetOwnershipRange(B,&rstart,&rend);
2864: for (i=rstart; i<rend; i++) {
2865: const PetscScalar *ba,*ca;
2866: const PetscInt *bj,*cj;
2867: PetscInt bn,cn,j,maxentrycol = -1,maxdiffcol = -1,maxrdiffcol = -1;
2868: PetscReal maxentry = 0,maxdiff = 0,maxrdiff = 0;
2869: MatGetRow(B,i,&bn,&bj,&ba);
2870: MatGetRow(Bfd,i,&cn,&cj,&ca);
2871: if (bn != cn) SETERRQ(((PetscObject)A)->comm,PETSC_ERR_PLIB,"Unexpected different nonzero pattern in -snes_compare_coloring_threshold");
2872: for (j=0; j<bn; j++) {
2873: PetscReal rdiff = PetscAbsScalar(ca[j]) / (threshold_atol + threshold_rtol*PetscAbsScalar(ba[j]));
2874: if (PetscAbsScalar(ba[j]) > PetscAbs(maxentry)) {
2875: maxentrycol = bj[j];
2876: maxentry = PetscRealPart(ba[j]);
2877: }
2878: if (PetscAbsScalar(ca[j]) > PetscAbs(maxdiff)) {
2879: maxdiffcol = bj[j];
2880: maxdiff = PetscRealPart(ca[j]);
2881: }
2882: if (rdiff > maxrdiff) {
2883: maxrdiffcol = bj[j];
2884: maxrdiff = rdiff;
2885: }
2886: }
2887: if (maxrdiff > 1) {
2888: 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);
2889: for (j=0; j<bn; j++) {
2890: PetscReal rdiff;
2891: rdiff = PetscAbsScalar(ca[j]) / (threshold_atol + threshold_rtol*PetscAbsScalar(ba[j]));
2892: if (rdiff > 1) {
2893: PetscViewerASCIIPrintf(vstdout," (%D,%g:%g)",bj[j],(double)PetscRealPart(ba[j]),(double)PetscRealPart(ca[j]));
2894: }
2895: }
2896: PetscViewerASCIIPrintf(vstdout,"\n",i,maxentry,maxdiff,maxrdiff);
2897: }
2898: MatRestoreRow(B,i,&bn,&bj,&ba);
2899: MatRestoreRow(Bfd,i,&cn,&cj,&ca);
2900: }
2901: }
2902: PetscViewerDestroy(&vdraw);
2903: MatDestroy(&Bfd);
2904: }
2905: }
2906: return(0);
2907: }
2909: /*MC
2910: SNESJacobianFunction - Function used to convey the nonlinear Jacobian of the function to be solved by SNES
2912: Synopsis:
2913: #include "petscsnes.h"
2914: PetscErrorCode SNESJacobianFunction(SNES snes,Vec x,Mat Amat,Mat Pmat,void *ctx);
2916: + x - input vector
2917: . Amat - the matrix that defines the (approximate) Jacobian
2918: . Pmat - the matrix to be used in constructing the preconditioner, usually the same as Amat.
2919: - ctx - [optional] user-defined Jacobian context
2921: Level: intermediate
2923: .seealso: SNESSetFunction(), SNESGetFunction(), SNESSetJacobian(), SNESGetJacobian()
2924: M*/
2926: /*@C
2927: SNESSetJacobian - Sets the function to compute Jacobian as well as the
2928: location to store the matrix.
2930: Logically Collective on SNES and Mat
2932: Input Parameters:
2933: + snes - the SNES context
2934: . Amat - the matrix that defines the (approximate) Jacobian
2935: . Pmat - the matrix to be used in constructing the preconditioner, usually the same as Amat.
2936: . J - Jacobian evaluation routine (if NULL then SNES retains any previously set value), see SNESJacobianFunction for details
2937: - ctx - [optional] user-defined context for private data for the
2938: Jacobian evaluation routine (may be NULL) (if NULL then SNES retains any previously set value)
2940: Notes:
2941: If the Amat matrix and Pmat matrix are different you must call MatAssemblyBegin/End() on
2942: each matrix.
2944: If you know the operator Amat has a null space you can use MatSetNullSpace() and MatSetTransposeNullSpace() to supply the null
2945: space to Amat and the KSP solvers will automatically use that null space as needed during the solution process.
2947: If using SNESComputeJacobianDefaultColor() to assemble a Jacobian, the ctx argument
2948: must be a MatFDColoring.
2950: Other defect-correction schemes can be used by computing a different matrix in place of the Jacobian. One common
2951: example is to use the "Picard linearization" which only differentiates through the highest order parts of each term.
2953: Level: beginner
2955: .keywords: SNES, nonlinear, set, Jacobian, matrix
2957: .seealso: KSPSetOperators(), SNESSetFunction(), MatMFFDComputeJacobian(), SNESComputeJacobianDefaultColor(), MatStructure, J,
2958: SNESSetPicard(), SNESJacobianFunction
2959: @*/
2960: PetscErrorCode SNESSetJacobian(SNES snes,Mat Amat,Mat Pmat,PetscErrorCode (*J)(SNES,Vec,Mat,Mat,void*),void *ctx)
2961: {
2963: DM dm;
2971: SNESGetDM(snes,&dm);
2972: DMSNESSetJacobian(dm,J,ctx);
2973: if (Amat) {
2974: PetscObjectReference((PetscObject)Amat);
2975: MatDestroy(&snes->jacobian);
2977: snes->jacobian = Amat;
2978: }
2979: if (Pmat) {
2980: PetscObjectReference((PetscObject)Pmat);
2981: MatDestroy(&snes->jacobian_pre);
2983: snes->jacobian_pre = Pmat;
2984: }
2985: return(0);
2986: }
2988: /*@C
2989: SNESGetJacobian - Returns the Jacobian matrix and optionally the user
2990: provided context for evaluating the Jacobian.
2992: Not Collective, but Mat object will be parallel if SNES object is
2994: Input Parameter:
2995: . snes - the nonlinear solver context
2997: Output Parameters:
2998: + Amat - location to stash (approximate) Jacobian matrix (or NULL)
2999: . Pmat - location to stash matrix used to compute the preconditioner (or NULL)
3000: . J - location to put Jacobian function (or NULL), see SNESJacobianFunction for details on its calling sequence
3001: - ctx - location to stash Jacobian ctx (or NULL)
3003: Level: advanced
3005: .seealso: SNESSetJacobian(), SNESComputeJacobian(), SNESJacobianFunction, SNESGetFunction()
3006: @*/
3007: PetscErrorCode SNESGetJacobian(SNES snes,Mat *Amat,Mat *Pmat,PetscErrorCode (**J)(SNES,Vec,Mat,Mat,void*),void **ctx)
3008: {
3010: DM dm;
3011: DMSNES sdm;
3015: if (Amat) *Amat = snes->jacobian;
3016: if (Pmat) *Pmat = snes->jacobian_pre;
3017: SNESGetDM(snes,&dm);
3018: DMGetDMSNES(dm,&sdm);
3019: if (J) *J = sdm->ops->computejacobian;
3020: if (ctx) *ctx = sdm->jacobianctx;
3021: return(0);
3022: }
3024: /*@
3025: SNESSetUp - Sets up the internal data structures for the later use
3026: of a nonlinear solver.
3028: Collective on SNES
3030: Input Parameters:
3031: . snes - the SNES context
3033: Notes:
3034: For basic use of the SNES solvers the user need not explicitly call
3035: SNESSetUp(), since these actions will automatically occur during
3036: the call to SNESSolve(). However, if one wishes to control this
3037: phase separately, SNESSetUp() should be called after SNESCreate()
3038: and optional routines of the form SNESSetXXX(), but before SNESSolve().
3040: Level: advanced
3042: .keywords: SNES, nonlinear, setup
3044: .seealso: SNESCreate(), SNESSolve(), SNESDestroy()
3045: @*/
3046: PetscErrorCode SNESSetUp(SNES snes)
3047: {
3049: DM dm;
3050: DMSNES sdm;
3051: SNESLineSearch linesearch, pclinesearch;
3052: void *lsprectx,*lspostctx;
3053: PetscErrorCode (*precheck)(SNESLineSearch,Vec,Vec,PetscBool*,void*);
3054: PetscErrorCode (*postcheck)(SNESLineSearch,Vec,Vec,Vec,PetscBool*,PetscBool*,void*);
3055: PetscErrorCode (*func)(SNES,Vec,Vec,void*);
3056: Vec f,fpc;
3057: void *funcctx;
3058: PetscErrorCode (*jac)(SNES,Vec,Mat,Mat,void*);
3059: void *jacctx,*appctx;
3060: Mat j,jpre;
3064: if (snes->setupcalled) return(0);
3066: if (!((PetscObject)snes)->type_name) {
3067: SNESSetType(snes,SNESNEWTONLS);
3068: }
3070: SNESGetFunction(snes,&snes->vec_func,NULL,NULL);
3072: SNESGetDM(snes,&dm);
3073: DMGetDMSNES(dm,&sdm);
3074: if (!sdm->ops->computefunction) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"Function never provided to SNES object");
3075: if (!sdm->ops->computejacobian) {
3076: DMSNESSetJacobian(dm,SNESComputeJacobianDefaultColor,NULL);
3077: }
3078: if (!snes->vec_func) {
3079: DMCreateGlobalVector(dm,&snes->vec_func);
3080: }
3082: if (!snes->ksp) {
3083: SNESGetKSP(snes, &snes->ksp);
3084: }
3086: if (!snes->linesearch) {
3087: SNESGetLineSearch(snes, &snes->linesearch);
3088: }
3089: SNESLineSearchSetFunction(snes->linesearch,SNESComputeFunction);
3091: if (snes->npc && (snes->npcside== PC_LEFT)) {
3092: snes->mf = PETSC_TRUE;
3093: snes->mf_operator = PETSC_FALSE;
3094: }
3096: if (snes->npc) {
3097: /* copy the DM over */
3098: SNESGetDM(snes,&dm);
3099: SNESSetDM(snes->npc,dm);
3101: SNESGetFunction(snes,&f,&func,&funcctx);
3102: VecDuplicate(f,&fpc);
3103: SNESSetFunction(snes->npc,fpc,func,funcctx);
3104: SNESGetJacobian(snes,&j,&jpre,&jac,&jacctx);
3105: SNESSetJacobian(snes->npc,j,jpre,jac,jacctx);
3106: SNESGetApplicationContext(snes,&appctx);
3107: SNESSetApplicationContext(snes->npc,appctx);
3108: VecDestroy(&fpc);
3110: /* copy the function pointers over */
3111: PetscObjectCopyFortranFunctionPointers((PetscObject)snes,(PetscObject)snes->npc);
3113: /* default to 1 iteration */
3114: SNESSetTolerances(snes->npc,0.0,0.0,0.0,1,snes->npc->max_funcs);
3115: if (snes->npcside==PC_RIGHT) {
3116: SNESSetNormSchedule(snes->npc,SNES_NORM_FINAL_ONLY);
3117: } else {
3118: SNESSetNormSchedule(snes->npc,SNES_NORM_NONE);
3119: }
3120: SNESSetFromOptions(snes->npc);
3122: /* copy the line search context over */
3123: SNESGetLineSearch(snes,&linesearch);
3124: SNESGetLineSearch(snes->npc,&pclinesearch);
3125: SNESLineSearchGetPreCheck(linesearch,&precheck,&lsprectx);
3126: SNESLineSearchGetPostCheck(linesearch,&postcheck,&lspostctx);
3127: SNESLineSearchSetPreCheck(pclinesearch,precheck,lsprectx);
3128: SNESLineSearchSetPostCheck(pclinesearch,postcheck,lspostctx);
3129: PetscObjectCopyFortranFunctionPointers((PetscObject)linesearch, (PetscObject)pclinesearch);
3130: }
3131: if (snes->mf) {
3132: SNESSetUpMatrixFree_Private(snes, snes->mf_operator, snes->mf_version);
3133: }
3134: if (snes->ops->usercompute && !snes->user) {
3135: (*snes->ops->usercompute)(snes,(void**)&snes->user);
3136: }
3138: snes->jac_iter = 0;
3139: snes->pre_iter = 0;
3141: if (snes->ops->setup) {
3142: (*snes->ops->setup)(snes);
3143: }
3145: if (snes->npc && (snes->npcside== PC_LEFT)) {
3146: if (snes->functype == SNES_FUNCTION_PRECONDITIONED) {
3147: SNESGetLineSearch(snes,&linesearch);
3148: SNESLineSearchSetFunction(linesearch,SNESComputeFunctionDefaultNPC);
3149: }
3150: }
3152: snes->setupcalled = PETSC_TRUE;
3153: return(0);
3154: }
3156: /*@
3157: SNESReset - Resets a SNES context to the snessetupcalled = 0 state and removes any allocated Vecs and Mats
3159: Collective on SNES
3161: Input Parameter:
3162: . snes - iterative context obtained from SNESCreate()
3164: Level: intermediate
3166: Notes:
3167: Also calls the application context destroy routine set with SNESSetComputeApplicationContext()
3169: .keywords: SNES, destroy
3171: .seealso: SNESCreate(), SNESSetUp(), SNESSolve()
3172: @*/
3173: PetscErrorCode SNESReset(SNES snes)
3174: {
3179: if (snes->ops->userdestroy && snes->user) {
3180: (*snes->ops->userdestroy)((void**)&snes->user);
3181: snes->user = NULL;
3182: }
3183: if (snes->npc) {
3184: SNESReset(snes->npc);
3185: }
3187: if (snes->ops->reset) {
3188: (*snes->ops->reset)(snes);
3189: }
3190: if (snes->ksp) {
3191: KSPReset(snes->ksp);
3192: }
3194: if (snes->linesearch) {
3195: SNESLineSearchReset(snes->linesearch);
3196: }
3198: VecDestroy(&snes->vec_rhs);
3199: VecDestroy(&snes->vec_sol);
3200: VecDestroy(&snes->vec_sol_update);
3201: VecDestroy(&snes->vec_func);
3202: MatDestroy(&snes->jacobian);
3203: MatDestroy(&snes->jacobian_pre);
3204: VecDestroyVecs(snes->nwork,&snes->work);
3205: VecDestroyVecs(snes->nvwork,&snes->vwork);
3207: snes->alwayscomputesfinalresidual = PETSC_FALSE;
3209: snes->nwork = snes->nvwork = 0;
3210: snes->setupcalled = PETSC_FALSE;
3211: return(0);
3212: }
3214: /*@
3215: SNESDestroy - Destroys the nonlinear solver context that was created
3216: with SNESCreate().
3218: Collective on SNES
3220: Input Parameter:
3221: . snes - the SNES context
3223: Level: beginner
3225: .keywords: SNES, nonlinear, destroy
3227: .seealso: SNESCreate(), SNESSolve()
3228: @*/
3229: PetscErrorCode SNESDestroy(SNES *snes)
3230: {
3234: if (!*snes) return(0);
3236: if (--((PetscObject)(*snes))->refct > 0) {*snes = 0; return(0);}
3238: SNESReset((*snes));
3239: SNESDestroy(&(*snes)->npc);
3241: /* if memory was published with SAWs then destroy it */
3242: PetscObjectSAWsViewOff((PetscObject)*snes);
3243: if ((*snes)->ops->destroy) {(*((*snes))->ops->destroy)((*snes));}
3245: if ((*snes)->dm) {DMCoarsenHookRemove((*snes)->dm,DMCoarsenHook_SNESVecSol,DMRestrictHook_SNESVecSol,*snes);}
3246: DMDestroy(&(*snes)->dm);
3247: KSPDestroy(&(*snes)->ksp);
3248: SNESLineSearchDestroy(&(*snes)->linesearch);
3250: PetscFree((*snes)->kspconvctx);
3251: if ((*snes)->ops->convergeddestroy) {
3252: (*(*snes)->ops->convergeddestroy)((*snes)->cnvP);
3253: }
3254: if ((*snes)->conv_malloc) {
3255: PetscFree((*snes)->conv_hist);
3256: PetscFree((*snes)->conv_hist_its);
3257: }
3258: SNESMonitorCancel((*snes));
3259: PetscHeaderDestroy(snes);
3260: return(0);
3261: }
3263: /* ----------- Routines to set solver parameters ---------- */
3265: /*@
3266: SNESSetLagPreconditioner - Determines when the preconditioner is rebuilt in the nonlinear solve.
3268: Logically Collective on SNES
3270: Input Parameters:
3271: + snes - the SNES context
3272: - lag - -1 indicates NEVER rebuild, 1 means rebuild every time the Jacobian is computed within a single nonlinear solve, 2 means every second time
3273: the Jacobian is built etc. -2 indicates rebuild preconditioner at next chance but then never rebuild after that
3275: Options Database Keys:
3276: . -snes_lag_preconditioner <lag>
3278: Notes:
3279: The default is 1
3280: The preconditioner is ALWAYS built in the first iteration of a nonlinear solve unless lag is -1
3281: If -1 is used before the very first nonlinear solve the preconditioner is still built because there is no previous preconditioner to use
3283: Level: intermediate
3285: .keywords: SNES, nonlinear, set, convergence, tolerances
3287: .seealso: SNESSetTrustRegionTolerance(), SNESGetLagPreconditioner(), SNESSetLagJacobian(), SNESGetLagJacobian()
3289: @*/
3290: PetscErrorCode SNESSetLagPreconditioner(SNES snes,PetscInt lag)
3291: {
3294: if (lag < -2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Lag must be -2, -1, 1 or greater");
3295: if (!lag) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Lag cannot be 0");
3297: snes->lagpreconditioner = lag;
3298: return(0);
3299: }
3301: /*@
3302: SNESSetGridSequence - sets the number of steps of grid sequencing that SNES does
3304: Logically Collective on SNES
3306: Input Parameters:
3307: + snes - the SNES context
3308: - steps - the number of refinements to do, defaults to 0
3310: Options Database Keys:
3311: . -snes_grid_sequence <steps>
3313: Level: intermediate
3315: Notes:
3316: Use SNESGetSolution() to extract the fine grid solution after grid sequencing.
3318: .keywords: SNES, nonlinear, set, convergence, tolerances
3320: .seealso: SNESSetTrustRegionTolerance(), SNESGetLagPreconditioner(), SNESSetLagJacobian(), SNESGetLagJacobian(), SNESGetGridSequence()
3322: @*/
3323: PetscErrorCode SNESSetGridSequence(SNES snes,PetscInt steps)
3324: {
3328: snes->gridsequence = steps;
3329: return(0);
3330: }
3332: /*@
3333: SNESGetGridSequence - gets the number of steps of grid sequencing that SNES does
3335: Logically Collective on SNES
3337: Input Parameter:
3338: . snes - the SNES context
3340: Output Parameter:
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: .keywords: SNES, nonlinear, set, convergence, tolerances
3353: .seealso: SNESSetTrustRegionTolerance(), SNESGetLagPreconditioner(), SNESSetLagJacobian(), SNESGetLagJacobian(), SNESSetGridSequence()
3355: @*/
3356: PetscErrorCode SNESGetGridSequence(SNES snes,PetscInt *steps)
3357: {
3360: *steps = snes->gridsequence;
3361: return(0);
3362: }
3364: /*@
3365: SNESGetLagPreconditioner - Indicates how often the preconditioner is rebuilt
3367: Not Collective
3369: Input Parameter:
3370: . snes - the SNES context
3372: Output Parameter:
3373: . lag - -1 indicates NEVER rebuild, 1 means rebuild every time the Jacobian is computed within a single nonlinear solve, 2 means every second time
3374: the Jacobian is built etc. -2 indicates rebuild preconditioner at next chance but then never rebuild after that
3376: Options Database Keys:
3377: . -snes_lag_preconditioner <lag>
3379: Notes:
3380: The default is 1
3381: The preconditioner is ALWAYS built in the first iteration of a nonlinear solve unless lag is -1
3383: Level: intermediate
3385: .keywords: SNES, nonlinear, set, convergence, tolerances
3387: .seealso: SNESSetTrustRegionTolerance(), SNESSetLagPreconditioner()
3389: @*/
3390: PetscErrorCode SNESGetLagPreconditioner(SNES snes,PetscInt *lag)
3391: {
3394: *lag = snes->lagpreconditioner;
3395: return(0);
3396: }
3398: /*@
3399: SNESSetLagJacobian - Determines when the Jacobian is rebuilt in the nonlinear solve. See SNESSetLagPreconditioner() for determining how
3400: often the preconditioner is rebuilt.
3402: Logically Collective on SNES
3404: Input Parameters:
3405: + snes - the SNES context
3406: - lag - -1 indicates NEVER rebuild, 1 means rebuild every time the Jacobian is computed within a single nonlinear solve, 2 means every second time
3407: the Jacobian is built etc. -2 means rebuild at next chance but then never again
3409: Options Database Keys:
3410: . -snes_lag_jacobian <lag>
3412: Notes:
3413: The default is 1
3414: The Jacobian is ALWAYS built in the first iteration of a nonlinear solve unless lag is -1
3415: 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
3416: at the next Newton step but never again (unless it is reset to another value)
3418: Level: intermediate
3420: .keywords: SNES, nonlinear, set, convergence, tolerances
3422: .seealso: SNESSetTrustRegionTolerance(), SNESGetLagPreconditioner(), SNESSetLagPreconditioner(), SNESGetLagJacobian()
3424: @*/
3425: PetscErrorCode SNESSetLagJacobian(SNES snes,PetscInt lag)
3426: {
3429: if (lag < -2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Lag must be -2, -1, 1 or greater");
3430: if (!lag) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Lag cannot be 0");
3432: snes->lagjacobian = lag;
3433: return(0);
3434: }
3436: /*@
3437: SNESGetLagJacobian - Indicates how often the Jacobian is rebuilt. See SNESGetLagPreconditioner() to determine when the preconditioner is rebuilt
3439: Not Collective
3441: Input Parameter:
3442: . snes - the SNES context
3444: Output Parameter:
3445: . lag - -1 indicates NEVER rebuild, 1 means rebuild every time the Jacobian is computed within a single nonlinear solve, 2 means every second time
3446: the Jacobian is built etc.
3448: Options Database Keys:
3449: . -snes_lag_jacobian <lag>
3451: Notes:
3452: The default is 1
3453: The jacobian is ALWAYS built in the first iteration of a nonlinear solve unless lag is -1
3455: Level: intermediate
3457: .keywords: SNES, nonlinear, set, convergence, tolerances
3459: .seealso: SNESSetTrustRegionTolerance(), SNESSetLagJacobian(), SNESSetLagPreconditioner(), SNESGetLagPreconditioner()
3461: @*/
3462: PetscErrorCode SNESGetLagJacobian(SNES snes,PetscInt *lag)
3463: {
3466: *lag = snes->lagjacobian;
3467: return(0);
3468: }
3470: /*@
3471: SNESSetLagJacobianPersists - Set whether or not the Jacobian lagging persists through multiple solves
3473: Logically collective on SNES
3475: Input Parameter:
3476: + snes - the SNES context
3477: - flg - jacobian lagging persists if true
3479: Options Database Keys:
3480: . -snes_lag_jacobian_persists <flg>
3482: Notes:
3483: This is useful both for nonlinear preconditioning, where it's appropriate to have the Jacobian be stale by
3484: several solves, and for implicit time-stepping, where Jacobian lagging in the inner nonlinear solve over several
3485: timesteps may present huge efficiency gains.
3487: Level: developer
3489: .keywords: SNES, nonlinear, lag
3491: .seealso: SNESSetLagPreconditionerPersists(), SNESSetLagJacobian(), SNESGetLagJacobian(), SNESGetNPC()
3493: @*/
3494: PetscErrorCode SNESSetLagJacobianPersists(SNES snes,PetscBool flg)
3495: {
3499: snes->lagjac_persist = flg;
3500: return(0);
3501: }
3503: /*@
3504: SNESSetLagPreconditionerPersists - Set whether or not the preconditioner lagging persists through multiple solves
3506: Logically Collective on SNES
3508: Input Parameter:
3509: + snes - the SNES context
3510: - flg - preconditioner lagging persists if true
3512: Options Database Keys:
3513: . -snes_lag_jacobian_persists <flg>
3515: Notes:
3516: This is useful both for nonlinear preconditioning, where it's appropriate to have the preconditioner be stale
3517: by several solves, and for implicit time-stepping, where preconditioner lagging in the inner nonlinear solve over
3518: several timesteps may present huge efficiency gains.
3520: Level: developer
3522: .keywords: SNES, nonlinear, lag
3524: .seealso: SNESSetLagJacobianPersists(), SNESSetLagJacobian(), SNESGetLagJacobian(), SNESGetNPC()
3526: @*/
3527: PetscErrorCode SNESSetLagPreconditionerPersists(SNES snes,PetscBool flg)
3528: {
3532: snes->lagpre_persist = flg;
3533: return(0);
3534: }
3536: /*@
3537: SNESSetForceIteration - force SNESSolve() to take at least one iteration regardless of the initial residual norm
3539: Logically Collective on SNES
3541: Input Parameters:
3542: + snes - the SNES context
3543: - force - PETSC_TRUE require at least one iteration
3545: Options Database Keys:
3546: . -snes_force_iteration <force> - Sets forcing an iteration
3548: Notes:
3549: This is used sometimes with TS to prevent TS from detecting a false steady state solution
3551: Level: intermediate
3553: .keywords: SNES, nonlinear, set, convergence, tolerances
3555: .seealso: SNESSetTrustRegionTolerance(), SNESSetDivergenceTolerance()
3556: @*/
3557: PetscErrorCode SNESSetForceIteration(SNES snes,PetscBool force)
3558: {
3561: snes->forceiteration = force;
3562: return(0);
3563: }
3565: /*@
3566: SNESGetForceIteration - Whether or not to force SNESSolve() take at least one iteration regardless of the initial residual norm
3568: Logically Collective on SNES
3570: Input Parameters:
3571: . snes - the SNES context
3573: Output Parameter:
3574: . force - PETSC_TRUE requires at least one iteration.
3576: .keywords: SNES, nonlinear, set, convergence, tolerances
3578: Level: intermediate
3580: .seealso: SNESSetForceIteration(), SNESSetTrustRegionTolerance(), SNESSetDivergenceTolerance()
3581: @*/
3582: PetscErrorCode SNESGetForceIteration(SNES snes,PetscBool *force)
3583: {
3586: *force = snes->forceiteration;
3587: return(0);
3588: }
3590: /*@
3591: SNESSetTolerances - Sets various parameters used in convergence tests.
3593: Logically Collective on SNES
3595: Input Parameters:
3596: + snes - the SNES context
3597: . abstol - absolute convergence tolerance
3598: . rtol - relative convergence tolerance
3599: . stol - convergence tolerance in terms of the norm of the change in the solution between steps, || delta x || < stol*|| x ||
3600: . maxit - maximum number of iterations
3601: - maxf - maximum number of function evaluations (-1 indicates no limit)
3603: Options Database Keys:
3604: + -snes_atol <abstol> - Sets abstol
3605: . -snes_rtol <rtol> - Sets rtol
3606: . -snes_stol <stol> - Sets stol
3607: . -snes_max_it <maxit> - Sets maxit
3608: - -snes_max_funcs <maxf> - Sets maxf
3610: Notes:
3611: The default maximum number of iterations is 50.
3612: The default maximum number of function evaluations is 1000.
3614: Level: intermediate
3616: .keywords: SNES, nonlinear, set, convergence, tolerances
3618: .seealso: SNESSetTrustRegionTolerance(), SNESSetDivergenceTolerance(), SNESSetForceIteration()
3619: @*/
3620: PetscErrorCode SNESSetTolerances(SNES snes,PetscReal abstol,PetscReal rtol,PetscReal stol,PetscInt maxit,PetscInt maxf)
3621: {
3630: if (abstol != PETSC_DEFAULT) {
3631: if (abstol < 0.0) SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_OUTOFRANGE,"Absolute tolerance %g must be non-negative",(double)abstol);
3632: snes->abstol = abstol;
3633: }
3634: if (rtol != PETSC_DEFAULT) {
3635: 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);
3636: snes->rtol = rtol;
3637: }
3638: if (stol != PETSC_DEFAULT) {
3639: if (stol < 0.0) SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_OUTOFRANGE,"Step tolerance %g must be non-negative",(double)stol);
3640: snes->stol = stol;
3641: }
3642: if (maxit != PETSC_DEFAULT) {
3643: if (maxit < 0) SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_OUTOFRANGE,"Maximum number of iterations %D must be non-negative",maxit);
3644: snes->max_its = maxit;
3645: }
3646: if (maxf != PETSC_DEFAULT) {
3647: if (maxf < -1) SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_OUTOFRANGE,"Maximum number of function evaluations %D must be -1 or nonnegative",maxf);
3648: snes->max_funcs = maxf;
3649: }
3650: snes->tolerancesset = PETSC_TRUE;
3651: return(0);
3652: }
3654: /*@
3655: SNESSetDivergenceTolerance - Sets the divergence tolerance used for the SNES divergence test.
3657: Logically Collective on SNES
3659: Input Parameters:
3660: + snes - the SNES context
3661: - divtol - the divergence tolerance. Use -1 to deactivate the test.
3663: Options Database Keys:
3664: + -snes_divergence_tolerance <divtol> - Sets divtol
3666: Notes:
3667: The default divergence tolerance is 1e4.
3669: Level: intermediate
3671: .keywords: SNES, nonlinear, set, divergence, tolerance
3673: .seealso: SNESSetTolerances(), SNESGetDivergenceTolerance
3674: @*/
3675: PetscErrorCode SNESSetDivergenceTolerance(SNES snes,PetscReal divtol)
3676: {
3681: if (divtol != PETSC_DEFAULT) {
3682: snes->divtol = divtol;
3683: }
3684: else {
3685: snes->divtol = 1.0e4;
3686: }
3687: return(0);
3688: }
3690: /*@
3691: SNESGetTolerances - Gets various parameters used in convergence tests.
3693: Not Collective
3695: Input Parameters:
3696: + snes - the SNES context
3697: . atol - absolute convergence tolerance
3698: . rtol - relative convergence tolerance
3699: . stol - convergence tolerance in terms of the norm
3700: of the change in the solution between steps
3701: . maxit - maximum number of iterations
3702: - maxf - maximum number of function evaluations
3704: Notes:
3705: The user can specify NULL for any parameter that is not needed.
3707: Level: intermediate
3709: .keywords: SNES, nonlinear, get, convergence, tolerances
3711: .seealso: SNESSetTolerances()
3712: @*/
3713: PetscErrorCode SNESGetTolerances(SNES snes,PetscReal *atol,PetscReal *rtol,PetscReal *stol,PetscInt *maxit,PetscInt *maxf)
3714: {
3717: if (atol) *atol = snes->abstol;
3718: if (rtol) *rtol = snes->rtol;
3719: if (stol) *stol = snes->stol;
3720: if (maxit) *maxit = snes->max_its;
3721: if (maxf) *maxf = snes->max_funcs;
3722: return(0);
3723: }
3725: /*@
3726: SNESGetDivergenceTolerance - Gets divergence tolerance used in divergence test.
3728: Not Collective
3730: Input Parameters:
3731: + snes - the SNES context
3732: - divtol - divergence tolerance
3734: Level: intermediate
3736: .keywords: SNES, nonlinear, get, divergence, tolerance
3738: .seealso: SNESSetDivergenceTolerance()
3739: @*/
3740: PetscErrorCode SNESGetDivergenceTolerance(SNES snes,PetscReal *divtol)
3741: {
3744: if (divtol) *divtol = snes->divtol;
3745: return(0);
3746: }
3748: /*@
3749: SNESSetTrustRegionTolerance - Sets the trust region parameter tolerance.
3751: Logically Collective on SNES
3753: Input Parameters:
3754: + snes - the SNES context
3755: - tol - tolerance
3757: Options Database Key:
3758: . -snes_trtol <tol> - Sets tol
3760: Level: intermediate
3762: .keywords: SNES, nonlinear, set, trust region, tolerance
3764: .seealso: SNESSetTolerances()
3765: @*/
3766: PetscErrorCode SNESSetTrustRegionTolerance(SNES snes,PetscReal tol)
3767: {
3771: snes->deltatol = tol;
3772: return(0);
3773: }
3775: /*
3776: Duplicate the lg monitors for SNES from KSP; for some reason with
3777: dynamic libraries things don't work under Sun4 if we just use
3778: macros instead of functions
3779: */
3780: PetscErrorCode SNESMonitorLGResidualNorm(SNES snes,PetscInt it,PetscReal norm,void *ctx)
3781: {
3786: KSPMonitorLGResidualNorm((KSP)snes,it,norm,ctx);
3787: return(0);
3788: }
3790: PetscErrorCode SNESMonitorLGCreate(MPI_Comm comm,const char host[],const char label[],int x,int y,int m,int n,PetscDrawLG *lgctx)
3791: {
3795: KSPMonitorLGResidualNormCreate(comm,host,label,x,y,m,n,lgctx);
3796: return(0);
3797: }
3799: PETSC_INTERN PetscErrorCode SNESMonitorRange_Private(SNES,PetscInt,PetscReal*);
3801: PetscErrorCode SNESMonitorLGRange(SNES snes,PetscInt n,PetscReal rnorm,void *monctx)
3802: {
3803: PetscDrawLG lg;
3804: PetscErrorCode ierr;
3805: PetscReal x,y,per;
3806: PetscViewer v = (PetscViewer)monctx;
3807: static PetscReal prev; /* should be in the context */
3808: PetscDraw draw;
3812: PetscViewerDrawGetDrawLG(v,0,&lg);
3813: if (!n) {PetscDrawLGReset(lg);}
3814: PetscDrawLGGetDraw(lg,&draw);
3815: PetscDrawSetTitle(draw,"Residual norm");
3816: x = (PetscReal)n;
3817: if (rnorm > 0.0) y = PetscLog10Real(rnorm);
3818: else y = -15.0;
3819: PetscDrawLGAddPoint(lg,&x,&y);
3820: if (n < 20 || !(n % 5) || snes->reason) {
3821: PetscDrawLGDraw(lg);
3822: PetscDrawLGSave(lg);
3823: }
3825: PetscViewerDrawGetDrawLG(v,1,&lg);
3826: if (!n) {PetscDrawLGReset(lg);}
3827: PetscDrawLGGetDraw(lg,&draw);
3828: PetscDrawSetTitle(draw,"% elemts > .2*max elemt");
3829: SNESMonitorRange_Private(snes,n,&per);
3830: x = (PetscReal)n;
3831: y = 100.0*per;
3832: PetscDrawLGAddPoint(lg,&x,&y);
3833: if (n < 20 || !(n % 5) || snes->reason) {
3834: PetscDrawLGDraw(lg);
3835: PetscDrawLGSave(lg);
3836: }
3838: PetscViewerDrawGetDrawLG(v,2,&lg);
3839: if (!n) {prev = rnorm;PetscDrawLGReset(lg);}
3840: PetscDrawLGGetDraw(lg,&draw);
3841: PetscDrawSetTitle(draw,"(norm -oldnorm)/oldnorm");
3842: x = (PetscReal)n;
3843: y = (prev - rnorm)/prev;
3844: PetscDrawLGAddPoint(lg,&x,&y);
3845: if (n < 20 || !(n % 5) || snes->reason) {
3846: PetscDrawLGDraw(lg);
3847: PetscDrawLGSave(lg);
3848: }
3850: PetscViewerDrawGetDrawLG(v,3,&lg);
3851: if (!n) {PetscDrawLGReset(lg);}
3852: PetscDrawLGGetDraw(lg,&draw);
3853: PetscDrawSetTitle(draw,"(norm -oldnorm)/oldnorm*(% > .2 max)");
3854: x = (PetscReal)n;
3855: y = (prev - rnorm)/(prev*per);
3856: if (n > 2) { /*skip initial crazy value */
3857: PetscDrawLGAddPoint(lg,&x,&y);
3858: }
3859: if (n < 20 || !(n % 5) || snes->reason) {
3860: PetscDrawLGDraw(lg);
3861: PetscDrawLGSave(lg);
3862: }
3863: prev = rnorm;
3864: return(0);
3865: }
3867: /*@
3868: SNESMonitor - runs the user provided monitor routines, if they exist
3870: Collective on SNES
3872: Input Parameters:
3873: + snes - nonlinear solver context obtained from SNESCreate()
3874: . iter - iteration number
3875: - rnorm - relative norm of the residual
3877: Notes:
3878: This routine is called by the SNES implementations.
3879: It does not typically need to be called by the user.
3881: Level: developer
3883: .seealso: SNESMonitorSet()
3884: @*/
3885: PetscErrorCode SNESMonitor(SNES snes,PetscInt iter,PetscReal rnorm)
3886: {
3888: PetscInt i,n = snes->numbermonitors;
3891: VecLockReadPush(snes->vec_sol);
3892: for (i=0; i<n; i++) {
3893: (*snes->monitor[i])(snes,iter,rnorm,snes->monitorcontext[i]);
3894: }
3895: VecLockReadPop(snes->vec_sol);
3896: return(0);
3897: }
3899: /* ------------ Routines to set performance monitoring options ----------- */
3901: /*MC
3902: SNESMonitorFunction - functional form passed to SNESMonitorSet() to monitor convergence of nonlinear solver
3904: Synopsis:
3905: #include <petscsnes.h>
3906: $ PetscErrorCode SNESMonitorFunction(SNES snes,PetscInt its, PetscReal norm,void *mctx)
3908: + snes - the SNES context
3909: . its - iteration number
3910: . norm - 2-norm function value (may be estimated)
3911: - mctx - [optional] monitoring context
3913: Level: advanced
3915: .seealso: SNESMonitorSet(), SNESMonitorGet()
3916: M*/
3918: /*@C
3919: SNESMonitorSet - Sets an ADDITIONAL function that is to be used at every
3920: iteration of the nonlinear solver to display the iteration's
3921: progress.
3923: Logically Collective on SNES
3925: Input Parameters:
3926: + snes - the SNES context
3927: . f - the monitor function, see SNESMonitorFunction for the calling sequence
3928: . mctx - [optional] user-defined context for private data for the
3929: monitor routine (use NULL if no context is desired)
3930: - monitordestroy - [optional] routine that frees monitor context
3931: (may be NULL)
3933: Options Database Keys:
3934: + -snes_monitor - sets SNESMonitorDefault()
3935: . -snes_monitor_lg_residualnorm - sets line graph monitor,
3936: uses SNESMonitorLGCreate()
3937: - -snes_monitor_cancel - cancels all monitors that have
3938: been hardwired into a code by
3939: calls to SNESMonitorSet(), but
3940: does not cancel those set via
3941: the options database.
3943: Notes:
3944: Several different monitoring routines may be set by calling
3945: SNESMonitorSet() multiple times; all will be called in the
3946: order in which they were set.
3948: Fortran Notes:
3949: Only a single monitor function can be set for each SNES object
3951: Level: intermediate
3953: .keywords: SNES, nonlinear, set, monitor
3955: .seealso: SNESMonitorDefault(), SNESMonitorCancel(), SNESMonitorFunction
3956: @*/
3957: PetscErrorCode SNESMonitorSet(SNES snes,PetscErrorCode (*f)(SNES,PetscInt,PetscReal,void*),void *mctx,PetscErrorCode (*monitordestroy)(void**))
3958: {
3959: PetscInt i;
3961: PetscBool identical;
3965: for (i=0; i<snes->numbermonitors;i++) {
3966: PetscMonitorCompare((PetscErrorCode (*)(void))f,mctx,monitordestroy,(PetscErrorCode (*)(void))snes->monitor[i],snes->monitorcontext[i],snes->monitordestroy[i],&identical);
3967: if (identical) return(0);
3968: }
3969: if (snes->numbermonitors >= MAXSNESMONITORS) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Too many monitors set");
3970: snes->monitor[snes->numbermonitors] = f;
3971: snes->monitordestroy[snes->numbermonitors] = monitordestroy;
3972: snes->monitorcontext[snes->numbermonitors++] = (void*)mctx;
3973: return(0);
3974: }
3976: /*@
3977: SNESMonitorCancel - Clears all the monitor functions for a SNES object.
3979: Logically Collective on SNES
3981: Input Parameters:
3982: . snes - the SNES context
3984: Options Database Key:
3985: . -snes_monitor_cancel - cancels all monitors that have been hardwired
3986: into a code by calls to SNESMonitorSet(), but does not cancel those
3987: set via the options database
3989: Notes:
3990: There is no way to clear one specific monitor from a SNES object.
3992: Level: intermediate
3994: .keywords: SNES, nonlinear, set, monitor
3996: .seealso: SNESMonitorDefault(), SNESMonitorSet()
3997: @*/
3998: PetscErrorCode SNESMonitorCancel(SNES snes)
3999: {
4001: PetscInt i;
4005: for (i=0; i<snes->numbermonitors; i++) {
4006: if (snes->monitordestroy[i]) {
4007: (*snes->monitordestroy[i])(&snes->monitorcontext[i]);
4008: }
4009: }
4010: snes->numbermonitors = 0;
4011: return(0);
4012: }
4014: /*MC
4015: SNESConvergenceTestFunction - functional form used for testing of convergence of nonlinear solver
4017: Synopsis:
4018: #include <petscsnes.h>
4019: $ PetscErrorCode SNESConvergenceTest(SNES snes,PetscInt it,PetscReal xnorm,PetscReal gnorm,PetscReal f,SNESConvergedReason *reason,void *cctx)
4021: + snes - the SNES context
4022: . it - current iteration (0 is the first and is before any Newton step)
4023: . cctx - [optional] convergence context
4024: . reason - reason for convergence/divergence
4025: . xnorm - 2-norm of current iterate
4026: . gnorm - 2-norm of current step
4027: - f - 2-norm of function
4029: Level: intermediate
4031: .seealso: SNESSetConvergenceTest(), SNESGetConvergenceTest()
4032: M*/
4034: /*@C
4035: SNESSetConvergenceTest - Sets the function that is to be used
4036: to test for convergence of the nonlinear iterative solution.
4038: Logically Collective on SNES
4040: Input Parameters:
4041: + snes - the SNES context
4042: . SNESConvergenceTestFunction - routine to test for convergence
4043: . cctx - [optional] context for private data for the convergence routine (may be NULL)
4044: - destroy - [optional] destructor for the context (may be NULL; PETSC_NULL_FUNCTION in Fortran)
4046: Level: advanced
4048: .keywords: SNES, nonlinear, set, convergence, test
4050: .seealso: SNESConvergedDefault(), SNESConvergedSkip(), SNESConvergenceTestFunction
4051: @*/
4052: PetscErrorCode SNESSetConvergenceTest(SNES snes,PetscErrorCode (*SNESConvergenceTestFunction)(SNES,PetscInt,PetscReal,PetscReal,PetscReal,SNESConvergedReason*,void*),void *cctx,PetscErrorCode (*destroy)(void*))
4053: {
4058: if (!SNESConvergenceTestFunction) SNESConvergenceTestFunction = SNESConvergedSkip;
4059: if (snes->ops->convergeddestroy) {
4060: (*snes->ops->convergeddestroy)(snes->cnvP);
4061: }
4062: snes->ops->converged = SNESConvergenceTestFunction;
4063: snes->ops->convergeddestroy = destroy;
4064: snes->cnvP = cctx;
4065: return(0);
4066: }
4068: /*@
4069: SNESGetConvergedReason - Gets the reason the SNES iteration was stopped.
4071: Not Collective
4073: Input Parameter:
4074: . snes - the SNES context
4076: Output Parameter:
4077: . reason - negative value indicates diverged, positive value converged, see SNESConvergedReason or the
4078: manual pages for the individual convergence tests for complete lists
4080: Options Database:
4081: . -snes_converged_reason - prints the reason to standard out
4083: Level: intermediate
4085: Notes:
4086: Should only be called after the call the SNESSolve() is complete, if it is called earlier it returns the value SNES__CONVERGED_ITERATING.
4088: .keywords: SNES, nonlinear, set, convergence, test
4090: .seealso: SNESSetConvergenceTest(), SNESSetConvergedReason(), SNESConvergedReason
4091: @*/
4092: PetscErrorCode SNESGetConvergedReason(SNES snes,SNESConvergedReason *reason)
4093: {
4097: *reason = snes->reason;
4098: return(0);
4099: }
4101: /*@
4102: SNESSetConvergedReason - Sets the reason the SNES iteration was stopped.
4104: Not Collective
4106: Input Parameters:
4107: + snes - the SNES context
4108: - reason - negative value indicates diverged, positive value converged, see SNESConvergedReason or the
4109: manual pages for the individual convergence tests for complete lists
4111: Level: intermediate
4113: .keywords: SNES, nonlinear, set, convergence, test
4114: .seealso: SNESGetConvergedReason(), SNESSetConvergenceTest(), SNESConvergedReason
4115: @*/
4116: PetscErrorCode SNESSetConvergedReason(SNES snes,SNESConvergedReason reason)
4117: {
4120: snes->reason = reason;
4121: return(0);
4122: }
4124: /*@
4125: SNESSetConvergenceHistory - Sets the array used to hold the convergence history.
4127: Logically Collective on SNES
4129: Input Parameters:
4130: + snes - iterative context obtained from SNESCreate()
4131: . a - array to hold history, this array will contain the function norms computed at each step
4132: . its - integer array holds the number of linear iterations for each solve.
4133: . na - size of a and its
4134: - reset - PETSC_TRUE indicates each new nonlinear solve resets the history counter to zero,
4135: else it continues storing new values for new nonlinear solves after the old ones
4137: Notes:
4138: If 'a' and 'its' are NULL then space is allocated for the history. If 'na' PETSC_DECIDE or PETSC_DEFAULT then a
4139: default array of length 10000 is allocated.
4141: This routine is useful, e.g., when running a code for purposes
4142: of accurate performance monitoring, when no I/O should be done
4143: during the section of code that is being timed.
4145: Level: intermediate
4147: .keywords: SNES, set, convergence, history
4149: .seealso: SNESGetConvergenceHistory()
4151: @*/
4152: PetscErrorCode SNESSetConvergenceHistory(SNES snes,PetscReal a[],PetscInt its[],PetscInt na,PetscBool reset)
4153: {
4160: if (!a) {
4161: if (na == PETSC_DECIDE || na == PETSC_DEFAULT) na = 1000;
4162: PetscCalloc1(na,&a);
4163: PetscCalloc1(na,&its);
4165: snes->conv_malloc = PETSC_TRUE;
4166: }
4167: snes->conv_hist = a;
4168: snes->conv_hist_its = its;
4169: snes->conv_hist_max = na;
4170: snes->conv_hist_len = 0;
4171: snes->conv_hist_reset = reset;
4172: return(0);
4173: }
4175: #if defined(PETSC_HAVE_MATLAB_ENGINE)
4176: #include <engine.h> /* MATLAB include file */
4177: #include <mex.h> /* MATLAB include file */
4179: PETSC_EXTERN mxArray *SNESGetConvergenceHistoryMatlab(SNES snes)
4180: {
4181: mxArray *mat;
4182: PetscInt i;
4183: PetscReal *ar;
4186: mat = mxCreateDoubleMatrix(snes->conv_hist_len,1,mxREAL);
4187: ar = (PetscReal*) mxGetData(mat);
4188: for (i=0; i<snes->conv_hist_len; i++) ar[i] = snes->conv_hist[i];
4189: PetscFunctionReturn(mat);
4190: }
4191: #endif
4193: /*@C
4194: SNESGetConvergenceHistory - Gets the array used to hold the convergence history.
4196: Not Collective
4198: Input Parameter:
4199: . snes - iterative context obtained from SNESCreate()
4201: Output Parameters:
4202: . a - array to hold history
4203: . its - integer array holds the number of linear iterations (or
4204: negative if not converged) for each solve.
4205: - na - size of a and its
4207: Notes:
4208: The calling sequence for this routine in Fortran is
4209: $ call SNESGetConvergenceHistory(SNES snes, integer na, integer ierr)
4211: This routine is useful, e.g., when running a code for purposes
4212: of accurate performance monitoring, when no I/O should be done
4213: during the section of code that is being timed.
4215: Level: intermediate
4217: .keywords: SNES, get, convergence, history
4219: .seealso: SNESSetConvergencHistory()
4221: @*/
4222: PetscErrorCode SNESGetConvergenceHistory(SNES snes,PetscReal *a[],PetscInt *its[],PetscInt *na)
4223: {
4226: if (a) *a = snes->conv_hist;
4227: if (its) *its = snes->conv_hist_its;
4228: if (na) *na = snes->conv_hist_len;
4229: return(0);
4230: }
4232: /*@C
4233: SNESSetUpdate - Sets the general-purpose update function called
4234: at the beginning of every iteration of the nonlinear solve. Specifically
4235: it is called just before the Jacobian is "evaluated".
4237: Logically Collective on SNES
4239: Input Parameters:
4240: . snes - The nonlinear solver context
4241: . func - The function
4243: Calling sequence of func:
4244: . func (SNES snes, PetscInt step);
4246: . step - The current step of the iteration
4248: Level: advanced
4250: 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()
4251: This is not used by most users.
4253: .keywords: SNES, update
4255: .seealso SNESSetJacobian(), SNESSolve()
4256: @*/
4257: PetscErrorCode SNESSetUpdate(SNES snes, PetscErrorCode (*func)(SNES, PetscInt))
4258: {
4261: snes->ops->update = func;
4262: return(0);
4263: }
4265: /*
4266: SNESScaleStep_Private - Scales a step so that its length is less than the
4267: positive parameter delta.
4269: Input Parameters:
4270: + snes - the SNES context
4271: . y - approximate solution of linear system
4272: . fnorm - 2-norm of current function
4273: - delta - trust region size
4275: Output Parameters:
4276: + gpnorm - predicted function norm at the new point, assuming local
4277: linearization. The value is zero if the step lies within the trust
4278: region, and exceeds zero otherwise.
4279: - ynorm - 2-norm of the step
4281: Note:
4282: For non-trust region methods such as SNESNEWTONLS, the parameter delta
4283: is set to be the maximum allowable step size.
4285: .keywords: SNES, nonlinear, scale, step
4286: */
4287: PetscErrorCode SNESScaleStep_Private(SNES snes,Vec y,PetscReal *fnorm,PetscReal *delta,PetscReal *gpnorm,PetscReal *ynorm)
4288: {
4289: PetscReal nrm;
4290: PetscScalar cnorm;
4298: VecNorm(y,NORM_2,&nrm);
4299: if (nrm > *delta) {
4300: nrm = *delta/nrm;
4301: *gpnorm = (1.0 - nrm)*(*fnorm);
4302: cnorm = nrm;
4303: VecScale(y,cnorm);
4304: *ynorm = *delta;
4305: } else {
4306: *gpnorm = 0.0;
4307: *ynorm = nrm;
4308: }
4309: return(0);
4310: }
4312: /*@
4313: SNESReasonView - Displays the reason a SNES solve converged or diverged to a viewer
4315: Collective on SNES
4317: Parameter:
4318: + snes - iterative context obtained from SNESCreate()
4319: - viewer - the viewer to display the reason
4322: Options Database Keys:
4323: . -snes_converged_reason - print reason for converged or diverged, also prints number of iterations
4325: Level: beginner
4327: .keywords: SNES, solve, linear system
4329: .seealso: SNESCreate(), SNESSetUp(), SNESDestroy(), SNESSetTolerances(), SNESConvergedDefault()
4331: @*/
4332: PetscErrorCode SNESReasonView(SNES snes,PetscViewer viewer)
4333: {
4334: PetscViewerFormat format;
4335: PetscBool isAscii;
4336: PetscErrorCode ierr;
4339: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isAscii);
4340: if (isAscii) {
4341: PetscViewerGetFormat(viewer, &format);
4342: PetscViewerASCIIAddTab(viewer,((PetscObject)snes)->tablevel);
4343: if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
4344: DM dm;
4345: Vec u;
4346: PetscDS prob;
4347: PetscInt Nf, f;
4348: PetscErrorCode (**exactFuncs)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar[], void *);
4349: PetscReal error;
4351: SNESGetDM(snes, &dm);
4352: SNESGetSolution(snes, &u);
4353: DMGetDS(dm, &prob);
4354: PetscDSGetNumFields(prob, &Nf);
4355: PetscMalloc1(Nf, &exactFuncs);
4356: for (f = 0; f < Nf; ++f) {PetscDSGetExactSolution(prob, f, &exactFuncs[f]);}
4357: DMComputeL2Diff(dm, 0.0, exactFuncs, NULL, u, &error);
4358: PetscFree(exactFuncs);
4359: if (error < 1.0e-11) {PetscViewerASCIIPrintf(viewer, "L_2 Error: < 1.0e-11\n");}
4360: else {PetscViewerASCIIPrintf(viewer, "L_2 Error: %g\n", error);}
4361: }
4362: if (snes->reason > 0) {
4363: if (((PetscObject) snes)->prefix) {
4364: PetscViewerASCIIPrintf(viewer,"Nonlinear %s solve converged due to %s iterations %D\n",((PetscObject) snes)->prefix,SNESConvergedReasons[snes->reason],snes->iter);
4365: } else {
4366: PetscViewerASCIIPrintf(viewer,"Nonlinear solve converged due to %s iterations %D\n",SNESConvergedReasons[snes->reason],snes->iter);
4367: }
4368: } else {
4369: if (((PetscObject) snes)->prefix) {
4370: PetscViewerASCIIPrintf(viewer,"Nonlinear %s solve did not converge due to %s iterations %D\n",((PetscObject) snes)->prefix,SNESConvergedReasons[snes->reason],snes->iter);
4371: } else {
4372: PetscViewerASCIIPrintf(viewer,"Nonlinear solve did not converge due to %s iterations %D\n",SNESConvergedReasons[snes->reason],snes->iter);
4373: }
4374: }
4375: PetscViewerASCIISubtractTab(viewer,((PetscObject)snes)->tablevel);
4376: }
4377: return(0);
4378: }
4380: /*@C
4381: SNESReasonViewFromOptions - Processes command line options to determine if/how a SNESReason is to be viewed.
4383: Collective on SNES
4385: Input Parameters:
4386: . snes - the SNES object
4388: Level: intermediate
4390: @*/
4391: PetscErrorCode SNESReasonViewFromOptions(SNES snes)
4392: {
4393: PetscErrorCode ierr;
4394: PetscViewer viewer;
4395: PetscBool flg;
4396: static PetscBool incall = PETSC_FALSE;
4397: PetscViewerFormat format;
4400: if (incall) return(0);
4401: incall = PETSC_TRUE;
4402: PetscOptionsGetViewer(PetscObjectComm((PetscObject)snes),((PetscObject)snes)->options,((PetscObject)snes)->prefix,"-snes_converged_reason",&viewer,&format,&flg);
4403: if (flg) {
4404: PetscViewerPushFormat(viewer,format);
4405: SNESReasonView(snes,viewer);
4406: PetscViewerPopFormat(viewer);
4407: PetscViewerDestroy(&viewer);
4408: }
4409: incall = PETSC_FALSE;
4410: return(0);
4411: }
4413: /*@
4414: SNESSolve - Solves a nonlinear system F(x) = b.
4415: Call SNESSolve() after calling SNESCreate() and optional routines of the form SNESSetXXX().
4417: Collective on SNES
4419: Input Parameters:
4420: + snes - the SNES context
4421: . b - the constant part of the equation F(x) = b, or NULL to use zero.
4422: - x - the solution vector.
4424: Notes:
4425: The user should initialize the vector,x, with the initial guess
4426: for the nonlinear solve prior to calling SNESSolve. In particular,
4427: to employ an initial guess of zero, the user should explicitly set
4428: this vector to zero by calling VecSet().
4430: Level: beginner
4432: .keywords: SNES, nonlinear, solve
4434: .seealso: SNESCreate(), SNESDestroy(), SNESSetFunction(), SNESSetJacobian(), SNESSetGridSequence(), SNESGetSolution()
4435: @*/
4436: PetscErrorCode SNESSolve(SNES snes,Vec b,Vec x)
4437: {
4438: PetscErrorCode ierr;
4439: PetscBool flg;
4440: PetscInt grid;
4441: Vec xcreated = NULL;
4442: DM dm;
4451: /* High level operations using the nonlinear solver */
4452: {
4453: PetscViewer viewer;
4454: PetscViewerFormat format;
4455: PetscInt num;
4456: PetscBool flg;
4457: static PetscBool incall = PETSC_FALSE;
4459: if (!incall) {
4460: /* Estimate the convergence rate of the discretization */
4461: PetscOptionsGetViewer(PetscObjectComm((PetscObject) snes),((PetscObject)snes)->options, ((PetscObject) snes)->prefix, "-snes_convergence_estimate", &viewer, &format, &flg);
4462: if (flg) {
4463: PetscConvEst conv;
4464: DM dm;
4465: PetscReal *alpha; /* Convergence rate of the solution error for each field in the L_2 norm */
4466: PetscInt Nf;
4468: incall = PETSC_TRUE;
4469: SNESGetDM(snes, &dm);
4470: DMGetNumFields(dm, &Nf);
4471: PetscMalloc1(Nf, &alpha);
4472: PetscConvEstCreate(PetscObjectComm((PetscObject) snes), &conv);
4473: PetscConvEstSetSolver(conv, snes);
4474: PetscConvEstSetFromOptions(conv);
4475: PetscConvEstSetUp(conv);
4476: PetscConvEstGetConvRate(conv, alpha);
4477: PetscViewerPushFormat(viewer, format);
4478: PetscConvEstRateView(conv, alpha, viewer);
4479: PetscViewerPopFormat(viewer);
4480: PetscViewerDestroy(&viewer);
4481: PetscConvEstDestroy(&conv);
4482: PetscFree(alpha);
4483: incall = PETSC_FALSE;
4484: }
4485: /* Adaptively refine the initial grid */
4486: num = 1;
4487: PetscOptionsGetInt(NULL, ((PetscObject) snes)->prefix, "-snes_adapt_initial", &num, &flg);
4488: if (flg) {
4489: DMAdaptor adaptor;
4491: incall = PETSC_TRUE;
4492: DMAdaptorCreate(PETSC_COMM_WORLD, &adaptor);
4493: DMAdaptorSetSolver(adaptor, snes);
4494: DMAdaptorSetSequenceLength(adaptor, num);
4495: DMAdaptorSetFromOptions(adaptor);
4496: DMAdaptorSetUp(adaptor);
4497: DMAdaptorAdapt(adaptor, x, DM_ADAPTATION_INITIAL, &dm, &x);
4498: DMAdaptorDestroy(&adaptor);
4499: incall = PETSC_FALSE;
4500: }
4501: /* Use grid sequencing to adapt */
4502: num = 0;
4503: PetscOptionsGetInt(NULL, ((PetscObject) snes)->prefix, "-snes_adapt_sequence", &num, NULL);
4504: if (num) {
4505: DMAdaptor adaptor;
4507: incall = PETSC_TRUE;
4508: DMAdaptorCreate(PETSC_COMM_WORLD, &adaptor);
4509: DMAdaptorSetSolver(adaptor, snes);
4510: DMAdaptorSetSequenceLength(adaptor, num);
4511: DMAdaptorSetFromOptions(adaptor);
4512: DMAdaptorSetUp(adaptor);
4513: DMAdaptorAdapt(adaptor, x, DM_ADAPTATION_SEQUENTIAL, &dm, &x);
4514: DMAdaptorDestroy(&adaptor);
4515: incall = PETSC_FALSE;
4516: }
4517: }
4518: }
4519: if (!x) {
4520: SNESGetDM(snes,&dm);
4521: DMCreateGlobalVector(dm,&xcreated);
4522: x = xcreated;
4523: }
4524: SNESViewFromOptions(snes,NULL,"-snes_view_pre");
4526: for (grid=0; grid<snes->gridsequence; grid++) {PetscViewerASCIIPushTab(PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)snes)));}
4527: for (grid=0; grid<snes->gridsequence+1; grid++) {
4529: /* set solution vector */
4530: if (!grid) {PetscObjectReference((PetscObject)x);}
4531: VecDestroy(&snes->vec_sol);
4532: snes->vec_sol = x;
4533: SNESGetDM(snes,&dm);
4535: /* set affine vector if provided */
4536: if (b) { PetscObjectReference((PetscObject)b); }
4537: VecDestroy(&snes->vec_rhs);
4538: snes->vec_rhs = b;
4540: 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");
4541: if (snes->vec_func == snes->vec_sol) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_IDN,"Solution vector cannot be function vector");
4542: if (snes->vec_rhs == snes->vec_sol) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_IDN,"Solution vector cannot be right hand side vector");
4543: if (!snes->vec_sol_update /* && snes->vec_sol */) {
4544: VecDuplicate(snes->vec_sol,&snes->vec_sol_update);
4545: PetscLogObjectParent((PetscObject)snes,(PetscObject)snes->vec_sol_update);
4546: }
4547: DMShellSetGlobalVector(dm,snes->vec_sol);
4548: SNESSetUp(snes);
4550: if (!grid) {
4551: if (snes->ops->computeinitialguess) {
4552: (*snes->ops->computeinitialguess)(snes,snes->vec_sol,snes->initialguessP);
4553: }
4554: }
4556: if (snes->conv_hist_reset) snes->conv_hist_len = 0;
4557: if (snes->counters_reset) {snes->nfuncs = 0; snes->linear_its = 0; snes->numFailures = 0;}
4559: PetscLogEventBegin(SNES_Solve,snes,0,0,0);
4560: (*snes->ops->solve)(snes);
4561: PetscLogEventEnd(SNES_Solve,snes,0,0,0);
4562: if (!snes->reason) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Internal error, solver returned without setting converged reason");
4563: snes->domainerror = PETSC_FALSE; /* clear the flag if it has been set */
4565: if (snes->lagjac_persist) snes->jac_iter += snes->iter;
4566: if (snes->lagpre_persist) snes->pre_iter += snes->iter;
4568: PetscOptionsGetViewer(PetscObjectComm((PetscObject)snes),((PetscObject)snes)->options,((PetscObject)snes)->prefix,"-snes_test_local_min",NULL,NULL,&flg);
4569: if (flg && !PetscPreLoadingOn) { SNESTestLocalMin(snes); }
4570: SNESReasonViewFromOptions(snes);
4572: if (snes->errorifnotconverged && snes->reason < 0) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_NOT_CONVERGED,"SNESSolve has not converged");
4573: if (snes->reason < 0) break;
4574: if (grid < snes->gridsequence) {
4575: DM fine;
4576: Vec xnew;
4577: Mat interp;
4579: DMRefine(snes->dm,PetscObjectComm((PetscObject)snes),&fine);
4580: if (!fine) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_INCOMP,"DMRefine() did not perform any refinement, cannot continue grid sequencing");
4581: DMCreateInterpolation(snes->dm,fine,&interp,NULL);
4582: DMCreateGlobalVector(fine,&xnew);
4583: MatInterpolate(interp,x,xnew);
4584: DMInterpolate(snes->dm,interp,fine);
4585: MatDestroy(&interp);
4586: x = xnew;
4588: SNESReset(snes);
4589: SNESSetDM(snes,fine);
4590: SNESResetFromOptions(snes);
4591: DMDestroy(&fine);
4592: PetscViewerASCIIPopTab(PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)snes)));
4593: }
4594: }
4595: SNESViewFromOptions(snes,NULL,"-snes_view");
4596: VecViewFromOptions(snes->vec_sol,(PetscObject)snes,"-snes_view_solution");
4598: VecDestroy(&xcreated);
4599: PetscObjectSAWsBlock((PetscObject)snes);
4600: return(0);
4601: }
4603: /* --------- Internal routines for SNES Package --------- */
4605: /*@C
4606: SNESSetType - Sets the method for the nonlinear solver.
4608: Collective on SNES
4610: Input Parameters:
4611: + snes - the SNES context
4612: - type - a known method
4614: Options Database Key:
4615: . -snes_type <type> - Sets the method; use -help for a list
4616: of available methods (for instance, newtonls or newtontr)
4618: Notes:
4619: See "petsc/include/petscsnes.h" for available methods (for instance)
4620: + SNESNEWTONLS - Newton's method with line search
4621: (systems of nonlinear equations)
4622: . SNESNEWTONTR - Newton's method with trust region
4623: (systems of nonlinear equations)
4625: Normally, it is best to use the SNESSetFromOptions() command and then
4626: set the SNES solver type from the options database rather than by using
4627: this routine. Using the options database provides the user with
4628: maximum flexibility in evaluating the many nonlinear solvers.
4629: The SNESSetType() routine is provided for those situations where it
4630: is necessary to set the nonlinear solver independently of the command
4631: line or options database. This might be the case, for example, when
4632: the choice of solver changes during the execution of the program,
4633: and the user's application is taking responsibility for choosing the
4634: appropriate method.
4636: Developer Notes:
4637: SNESRegister() adds a constructor for a new SNESType to SNESList, SNESSetType() locates
4638: the constructor in that list and calls it to create the spexific object.
4640: Level: intermediate
4642: .keywords: SNES, set, type
4644: .seealso: SNESType, SNESCreate(), SNESDestroy(), SNESGetType(), SNESSetFromOptions()
4646: @*/
4647: PetscErrorCode SNESSetType(SNES snes,SNESType type)
4648: {
4649: PetscErrorCode ierr,(*r)(SNES);
4650: PetscBool match;
4656: PetscObjectTypeCompare((PetscObject)snes,type,&match);
4657: if (match) return(0);
4659: PetscFunctionListFind(SNESList,type,&r);
4660: if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unable to find requested SNES type %s",type);
4661: /* Destroy the previous private SNES context */
4662: if (snes->ops->destroy) {
4663: (*(snes)->ops->destroy)(snes);
4664: snes->ops->destroy = NULL;
4665: }
4666: /* Reinitialize function pointers in SNESOps structure */
4667: snes->ops->setup = 0;
4668: snes->ops->solve = 0;
4669: snes->ops->view = 0;
4670: snes->ops->setfromoptions = 0;
4671: snes->ops->destroy = 0;
4672: /* Call the SNESCreate_XXX routine for this particular Nonlinear solver */
4673: snes->setupcalled = PETSC_FALSE;
4675: PetscObjectChangeTypeName((PetscObject)snes,type);
4676: (*r)(snes);
4677: return(0);
4678: }
4680: /*@C
4681: SNESGetType - Gets the SNES method type and name (as a string).
4683: Not Collective
4685: Input Parameter:
4686: . snes - nonlinear solver context
4688: Output Parameter:
4689: . type - SNES method (a character string)
4691: Level: intermediate
4693: .keywords: SNES, nonlinear, get, type, name
4694: @*/
4695: PetscErrorCode SNESGetType(SNES snes,SNESType *type)
4696: {
4700: *type = ((PetscObject)snes)->type_name;
4701: return(0);
4702: }
4704: /*@
4705: SNESSetSolution - Sets the solution vector for use by the SNES routines.
4707: Logically Collective on SNES and Vec
4709: Input Parameters:
4710: + snes - the SNES context obtained from SNESCreate()
4711: - u - the solution vector
4713: Level: beginner
4715: .keywords: SNES, set, solution
4716: @*/
4717: PetscErrorCode SNESSetSolution(SNES snes, Vec u)
4718: {
4719: DM dm;
4725: PetscObjectReference((PetscObject) u);
4726: VecDestroy(&snes->vec_sol);
4728: snes->vec_sol = u;
4730: SNESGetDM(snes, &dm);
4731: DMShellSetGlobalVector(dm, u);
4732: return(0);
4733: }
4735: /*@
4736: SNESGetSolution - Returns the vector where the approximate solution is
4737: stored. This is the fine grid solution when using SNESSetGridSequence().
4739: Not Collective, but Vec is parallel if SNES is parallel
4741: Input Parameter:
4742: . snes - the SNES context
4744: Output Parameter:
4745: . x - the solution
4747: Level: intermediate
4749: .keywords: SNES, nonlinear, get, solution
4751: .seealso: SNESGetSolutionUpdate(), SNESGetFunction()
4752: @*/
4753: PetscErrorCode SNESGetSolution(SNES snes,Vec *x)
4754: {
4758: *x = snes->vec_sol;
4759: return(0);
4760: }
4762: /*@
4763: SNESGetSolutionUpdate - Returns the vector where the solution update is
4764: stored.
4766: Not Collective, but Vec is parallel if SNES is parallel
4768: Input Parameter:
4769: . snes - the SNES context
4771: Output Parameter:
4772: . x - the solution update
4774: Level: advanced
4776: .keywords: SNES, nonlinear, get, solution, update
4778: .seealso: SNESGetSolution(), SNESGetFunction()
4779: @*/
4780: PetscErrorCode SNESGetSolutionUpdate(SNES snes,Vec *x)
4781: {
4785: *x = snes->vec_sol_update;
4786: return(0);
4787: }
4789: /*@C
4790: SNESGetFunction - Returns the vector where the function is stored.
4792: Not Collective, but Vec is parallel if SNES is parallel. Collective if Vec is requested, but has not been created yet.
4794: Input Parameter:
4795: . snes - the SNES context
4797: Output Parameter:
4798: + r - the vector that is used to store residuals (or NULL if you don't want it)
4799: . f - the function (or NULL if you don't want it); see SNESFunction for calling sequence details
4800: - ctx - the function context (or NULL if you don't want it)
4802: Level: advanced
4804: Notes: The vector r DOES NOT, in general contain the current value of the SNES nonlinear function
4806: .keywords: SNES, nonlinear, get, function
4808: .seealso: SNESSetFunction(), SNESGetSolution(), SNESFunction
4809: @*/
4810: PetscErrorCode SNESGetFunction(SNES snes,Vec *r,PetscErrorCode (**f)(SNES,Vec,Vec,void*),void **ctx)
4811: {
4813: DM dm;
4817: if (r) {
4818: if (!snes->vec_func) {
4819: if (snes->vec_rhs) {
4820: VecDuplicate(snes->vec_rhs,&snes->vec_func);
4821: } else if (snes->vec_sol) {
4822: VecDuplicate(snes->vec_sol,&snes->vec_func);
4823: } else if (snes->dm) {
4824: DMCreateGlobalVector(snes->dm,&snes->vec_func);
4825: }
4826: }
4827: *r = snes->vec_func;
4828: }
4829: SNESGetDM(snes,&dm);
4830: DMSNESGetFunction(dm,f,ctx);
4831: return(0);
4832: }
4834: /*@C
4835: SNESGetNGS - Returns the NGS function and context.
4837: Input Parameter:
4838: . snes - the SNES context
4840: Output Parameter:
4841: + f - the function (or NULL) see SNESNGSFunction for details
4842: - ctx - the function context (or NULL)
4844: Level: advanced
4846: .keywords: SNES, nonlinear, get, function
4848: .seealso: SNESSetNGS(), SNESGetFunction()
4849: @*/
4851: PetscErrorCode SNESGetNGS (SNES snes, PetscErrorCode (**f)(SNES, Vec, Vec, void*), void ** ctx)
4852: {
4854: DM dm;
4858: SNESGetDM(snes,&dm);
4859: DMSNESGetNGS(dm,f,ctx);
4860: return(0);
4861: }
4863: /*@C
4864: SNESSetOptionsPrefix - Sets the prefix used for searching for all
4865: SNES options in the database.
4867: Logically Collective on SNES
4869: Input Parameter:
4870: + snes - the SNES context
4871: - prefix - the prefix to prepend to all option names
4873: Notes:
4874: A hyphen (-) must NOT be given at the beginning of the prefix name.
4875: The first character of all runtime options is AUTOMATICALLY the hyphen.
4877: Level: advanced
4879: .keywords: SNES, set, options, prefix, database
4881: .seealso: SNESSetFromOptions()
4882: @*/
4883: PetscErrorCode SNESSetOptionsPrefix(SNES snes,const char prefix[])
4884: {
4889: PetscObjectSetOptionsPrefix((PetscObject)snes,prefix);
4890: if (!snes->ksp) {SNESGetKSP(snes,&snes->ksp);}
4891: if (snes->linesearch) {
4892: SNESGetLineSearch(snes,&snes->linesearch);
4893: PetscObjectSetOptionsPrefix((PetscObject)snes->linesearch,prefix);
4894: }
4895: KSPSetOptionsPrefix(snes->ksp,prefix);
4896: return(0);
4897: }
4899: /*@C
4900: SNESAppendOptionsPrefix - Appends to the prefix used for searching for all
4901: SNES options in the database.
4903: Logically Collective on SNES
4905: Input Parameters:
4906: + snes - the SNES context
4907: - prefix - the prefix to prepend to all option names
4909: Notes:
4910: A hyphen (-) must NOT be given at the beginning of the prefix name.
4911: The first character of all runtime options is AUTOMATICALLY the hyphen.
4913: Level: advanced
4915: .keywords: SNES, append, options, prefix, database
4917: .seealso: SNESGetOptionsPrefix()
4918: @*/
4919: PetscErrorCode SNESAppendOptionsPrefix(SNES snes,const char prefix[])
4920: {
4925: PetscObjectAppendOptionsPrefix((PetscObject)snes,prefix);
4926: if (!snes->ksp) {SNESGetKSP(snes,&snes->ksp);}
4927: if (snes->linesearch) {
4928: SNESGetLineSearch(snes,&snes->linesearch);
4929: PetscObjectAppendOptionsPrefix((PetscObject)snes->linesearch,prefix);
4930: }
4931: KSPAppendOptionsPrefix(snes->ksp,prefix);
4932: return(0);
4933: }
4935: /*@C
4936: SNESGetOptionsPrefix - Sets the prefix used for searching for all
4937: SNES options in the database.
4939: Not Collective
4941: Input Parameter:
4942: . snes - the SNES context
4944: Output Parameter:
4945: . prefix - pointer to the prefix string used
4947: Notes:
4948: On the fortran side, the user should pass in a string 'prefix' of
4949: sufficient length to hold the prefix.
4951: Level: advanced
4953: .keywords: SNES, get, options, prefix, database
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: .keywords: SNES, nonlinear, register
4996: .seealso: SNESRegisterAll(), SNESRegisterDestroy()
4998: Level: advanced
4999: @*/
5000: PetscErrorCode SNESRegister(const char sname[],PetscErrorCode (*function)(SNES))
5001: {
5005: SNESInitializePackage();
5006: PetscFunctionListAdd(&SNESList,sname,function);
5007: return(0);
5008: }
5010: PetscErrorCode SNESTestLocalMin(SNES snes)
5011: {
5013: PetscInt N,i,j;
5014: Vec u,uh,fh;
5015: PetscScalar value;
5016: PetscReal norm;
5019: SNESGetSolution(snes,&u);
5020: VecDuplicate(u,&uh);
5021: VecDuplicate(u,&fh);
5023: /* currently only works for sequential */
5024: PetscPrintf(PETSC_COMM_WORLD,"Testing FormFunction() for local min\n");
5025: VecGetSize(u,&N);
5026: for (i=0; i<N; i++) {
5027: VecCopy(u,uh);
5028: PetscPrintf(PETSC_COMM_WORLD,"i = %D\n",i);
5029: for (j=-10; j<11; j++) {
5030: value = PetscSign(j)*PetscExpReal(PetscAbs(j)-10.0);
5031: VecSetValue(uh,i,value,ADD_VALUES);
5032: SNESComputeFunction(snes,uh,fh);
5033: VecNorm(fh,NORM_2,&norm);
5034: PetscPrintf(PETSC_COMM_WORLD," j norm %D %18.16e\n",j,norm);
5035: value = -value;
5036: VecSetValue(uh,i,value,ADD_VALUES);
5037: }
5038: }
5039: VecDestroy(&uh);
5040: VecDestroy(&fh);
5041: return(0);
5042: }
5044: /*@
5045: SNESKSPSetUseEW - Sets SNES use Eisenstat-Walker method for
5046: computing relative tolerance for linear solvers within an inexact
5047: Newton method.
5049: Logically Collective on SNES
5051: Input Parameters:
5052: + snes - SNES context
5053: - flag - PETSC_TRUE or PETSC_FALSE
5055: Options Database:
5056: + -snes_ksp_ew - use Eisenstat-Walker method for determining linear system convergence
5057: . -snes_ksp_ew_version ver - version of Eisenstat-Walker method
5058: . -snes_ksp_ew_rtol0 <rtol0> - Sets rtol0
5059: . -snes_ksp_ew_rtolmax <rtolmax> - Sets rtolmax
5060: . -snes_ksp_ew_gamma <gamma> - Sets gamma
5061: . -snes_ksp_ew_alpha <alpha> - Sets alpha
5062: . -snes_ksp_ew_alpha2 <alpha2> - Sets alpha2
5063: - -snes_ksp_ew_threshold <threshold> - Sets threshold
5065: Notes:
5066: Currently, the default is to use a constant relative tolerance for
5067: the inner linear solvers. Alternatively, one can use the
5068: Eisenstat-Walker method, where the relative convergence tolerance
5069: is reset at each Newton iteration according progress of the nonlinear
5070: solver.
5072: Level: advanced
5074: Reference:
5075: S. C. Eisenstat and H. F. Walker, "Choosing the forcing terms in an
5076: inexact Newton method", SISC 17 (1), pp.16-32, 1996.
5078: .keywords: SNES, KSP, Eisenstat, Walker, convergence, test, inexact, Newton
5080: .seealso: SNESKSPGetUseEW(), SNESKSPGetParametersEW(), SNESKSPSetParametersEW()
5081: @*/
5082: PetscErrorCode SNESKSPSetUseEW(SNES snes,PetscBool flag)
5083: {
5087: snes->ksp_ewconv = flag;
5088: return(0);
5089: }
5091: /*@
5092: SNESKSPGetUseEW - Gets if SNES is using Eisenstat-Walker method
5093: for computing relative tolerance for linear solvers within an
5094: inexact Newton method.
5096: Not Collective
5098: Input Parameter:
5099: . snes - SNES context
5101: Output Parameter:
5102: . flag - PETSC_TRUE or PETSC_FALSE
5104: Notes:
5105: Currently, the default is to use a constant relative tolerance for
5106: the inner linear solvers. Alternatively, one can use the
5107: Eisenstat-Walker method, where the relative convergence tolerance
5108: is reset at each Newton iteration according progress of the nonlinear
5109: solver.
5111: Level: advanced
5113: Reference:
5114: S. C. Eisenstat and H. F. Walker, "Choosing the forcing terms in an
5115: inexact Newton method", SISC 17 (1), pp.16-32, 1996.
5117: .keywords: SNES, KSP, Eisenstat, Walker, convergence, test, inexact, Newton
5119: .seealso: SNESKSPSetUseEW(), SNESKSPGetParametersEW(), SNESKSPSetParametersEW()
5120: @*/
5121: PetscErrorCode SNESKSPGetUseEW(SNES snes, PetscBool *flag)
5122: {
5126: *flag = snes->ksp_ewconv;
5127: return(0);
5128: }
5130: /*@
5131: SNESKSPSetParametersEW - Sets parameters for Eisenstat-Walker
5132: convergence criteria for the linear solvers within an inexact
5133: Newton method.
5135: Logically Collective on SNES
5137: Input Parameters:
5138: + snes - SNES context
5139: . version - version 1, 2 (default is 2) or 3
5140: . rtol_0 - initial relative tolerance (0 <= rtol_0 < 1)
5141: . rtol_max - maximum relative tolerance (0 <= rtol_max < 1)
5142: . gamma - multiplicative factor for version 2 rtol computation
5143: (0 <= gamma2 <= 1)
5144: . alpha - power for version 2 rtol computation (1 < alpha <= 2)
5145: . alpha2 - power for safeguard
5146: - threshold - threshold for imposing safeguard (0 < threshold < 1)
5148: Note:
5149: Version 3 was contributed by Luis Chacon, June 2006.
5151: Use PETSC_DEFAULT to retain the default for any of the parameters.
5153: Level: advanced
5155: Reference:
5156: S. C. Eisenstat and H. F. Walker, "Choosing the forcing terms in an
5157: inexact Newton method", Utah State University Math. Stat. Dept. Res.
5158: Report 6/94/75, June, 1994, to appear in SIAM J. Sci. Comput.
5160: .keywords: SNES, KSP, Eisenstat, Walker, set, parameters
5162: .seealso: SNESKSPSetUseEW(), SNESKSPGetUseEW(), SNESKSPGetParametersEW()
5163: @*/
5164: PetscErrorCode SNESKSPSetParametersEW(SNES snes,PetscInt version,PetscReal rtol_0,PetscReal rtol_max,PetscReal gamma,PetscReal alpha,PetscReal alpha2,PetscReal threshold)
5165: {
5166: SNESKSPEW *kctx;
5170: kctx = (SNESKSPEW*)snes->kspconvctx;
5171: if (!kctx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"No Eisenstat-Walker context existing");
5180: if (version != PETSC_DEFAULT) kctx->version = version;
5181: if (rtol_0 != PETSC_DEFAULT) kctx->rtol_0 = rtol_0;
5182: if (rtol_max != PETSC_DEFAULT) kctx->rtol_max = rtol_max;
5183: if (gamma != PETSC_DEFAULT) kctx->gamma = gamma;
5184: if (alpha != PETSC_DEFAULT) kctx->alpha = alpha;
5185: if (alpha2 != PETSC_DEFAULT) kctx->alpha2 = alpha2;
5186: if (threshold != PETSC_DEFAULT) kctx->threshold = threshold;
5188: 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);
5189: 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);
5190: 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);
5191: 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);
5192: 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);
5193: 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);
5194: return(0);
5195: }
5197: /*@
5198: SNESKSPGetParametersEW - Gets parameters for Eisenstat-Walker
5199: convergence criteria for the linear solvers within an inexact
5200: Newton method.
5202: Not Collective
5204: Input Parameters:
5205: snes - SNES context
5207: Output Parameters:
5208: + version - version 1, 2 (default is 2) or 3
5209: . rtol_0 - initial relative tolerance (0 <= rtol_0 < 1)
5210: . rtol_max - maximum relative tolerance (0 <= rtol_max < 1)
5211: . gamma - multiplicative factor for version 2 rtol computation (0 <= gamma2 <= 1)
5212: . alpha - power for version 2 rtol computation (1 < alpha <= 2)
5213: . alpha2 - power for safeguard
5214: - threshold - threshold for imposing safeguard (0 < threshold < 1)
5216: Level: advanced
5218: .keywords: SNES, KSP, Eisenstat, Walker, get, parameters
5220: .seealso: SNESKSPSetUseEW(), SNESKSPGetUseEW(), SNESKSPSetParametersEW()
5221: @*/
5222: PetscErrorCode SNESKSPGetParametersEW(SNES snes,PetscInt *version,PetscReal *rtol_0,PetscReal *rtol_max,PetscReal *gamma,PetscReal *alpha,PetscReal *alpha2,PetscReal *threshold)
5223: {
5224: SNESKSPEW *kctx;
5228: kctx = (SNESKSPEW*)snes->kspconvctx;
5229: if (!kctx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"No Eisenstat-Walker context existing");
5230: if (version) *version = kctx->version;
5231: if (rtol_0) *rtol_0 = kctx->rtol_0;
5232: if (rtol_max) *rtol_max = kctx->rtol_max;
5233: if (gamma) *gamma = kctx->gamma;
5234: if (alpha) *alpha = kctx->alpha;
5235: if (alpha2) *alpha2 = kctx->alpha2;
5236: if (threshold) *threshold = kctx->threshold;
5237: return(0);
5238: }
5240: PetscErrorCode KSPPreSolve_SNESEW(KSP ksp, Vec b, Vec x, SNES snes)
5241: {
5243: SNESKSPEW *kctx = (SNESKSPEW*)snes->kspconvctx;
5244: PetscReal rtol = PETSC_DEFAULT,stol;
5247: if (!snes->ksp_ewconv) return(0);
5248: if (!snes->iter) {
5249: rtol = kctx->rtol_0; /* first time in, so use the original user rtol */
5250: VecNorm(snes->vec_func,NORM_2,&kctx->norm_first);
5251: }
5252: else {
5253: if (kctx->version == 1) {
5254: rtol = (snes->norm - kctx->lresid_last)/kctx->norm_last;
5255: if (rtol < 0.0) rtol = -rtol;
5256: stol = PetscPowReal(kctx->rtol_last,kctx->alpha2);
5257: if (stol > kctx->threshold) rtol = PetscMax(rtol,stol);
5258: } else if (kctx->version == 2) {
5259: rtol = kctx->gamma * PetscPowReal(snes->norm/kctx->norm_last,kctx->alpha);
5260: stol = kctx->gamma * PetscPowReal(kctx->rtol_last,kctx->alpha);
5261: if (stol > kctx->threshold) rtol = PetscMax(rtol,stol);
5262: } else if (kctx->version == 3) { /* contributed by Luis Chacon, June 2006. */
5263: rtol = kctx->gamma * PetscPowReal(snes->norm/kctx->norm_last,kctx->alpha);
5264: /* safeguard: avoid sharp decrease of rtol */
5265: stol = kctx->gamma*PetscPowReal(kctx->rtol_last,kctx->alpha);
5266: stol = PetscMax(rtol,stol);
5267: rtol = PetscMin(kctx->rtol_0,stol);
5268: /* safeguard: avoid oversolving */
5269: stol = kctx->gamma*(kctx->norm_first*snes->rtol)/snes->norm;
5270: stol = PetscMax(rtol,stol);
5271: rtol = PetscMin(kctx->rtol_0,stol);
5272: } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Only versions 1, 2 or 3 are supported: %D",kctx->version);
5273: }
5274: /* safeguard: avoid rtol greater than one */
5275: rtol = PetscMin(rtol,kctx->rtol_max);
5276: KSPSetTolerances(ksp,rtol,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
5277: PetscInfo3(snes,"iter %D, Eisenstat-Walker (version %D) KSP rtol=%g\n",snes->iter,kctx->version,(double)rtol);
5278: return(0);
5279: }
5281: PetscErrorCode KSPPostSolve_SNESEW(KSP ksp, Vec b, Vec x, SNES snes)
5282: {
5284: SNESKSPEW *kctx = (SNESKSPEW*)snes->kspconvctx;
5285: PCSide pcside;
5286: Vec lres;
5289: if (!snes->ksp_ewconv) return(0);
5290: KSPGetTolerances(ksp,&kctx->rtol_last,0,0,0);
5291: kctx->norm_last = snes->norm;
5292: if (kctx->version == 1) {
5293: PC pc;
5294: PetscBool isNone;
5296: KSPGetPC(ksp, &pc);
5297: PetscObjectTypeCompare((PetscObject) pc, PCNONE, &isNone);
5298: KSPGetPCSide(ksp,&pcside);
5299: if (pcside == PC_RIGHT || isNone) { /* XXX Should we also test KSP_UNPRECONDITIONED_NORM ? */
5300: /* KSP residual is true linear residual */
5301: KSPGetResidualNorm(ksp,&kctx->lresid_last);
5302: } else {
5303: /* KSP residual is preconditioned residual */
5304: /* compute true linear residual norm */
5305: VecDuplicate(b,&lres);
5306: MatMult(snes->jacobian,x,lres);
5307: VecAYPX(lres,-1.0,b);
5308: VecNorm(lres,NORM_2,&kctx->lresid_last);
5309: VecDestroy(&lres);
5310: }
5311: }
5312: return(0);
5313: }
5315: /*@
5316: SNESGetKSP - Returns the KSP context for a SNES solver.
5318: Not Collective, but if SNES object is parallel, then KSP object is parallel
5320: Input Parameter:
5321: . snes - the SNES context
5323: Output Parameter:
5324: . ksp - the KSP context
5326: Notes:
5327: The user can then directly manipulate the KSP context to set various
5328: options, etc. Likewise, the user can then extract and manipulate the
5329: PC contexts as well.
5331: Level: beginner
5333: .keywords: SNES, nonlinear, get, KSP, context
5335: .seealso: KSPGetPC(), SNESCreate(), KSPCreate(), SNESSetKSP()
5336: @*/
5337: PetscErrorCode SNESGetKSP(SNES snes,KSP *ksp)
5338: {
5345: if (!snes->ksp) {
5346: PetscBool monitor = PETSC_FALSE;
5348: KSPCreate(PetscObjectComm((PetscObject)snes),&snes->ksp);
5349: PetscObjectIncrementTabLevel((PetscObject)snes->ksp,(PetscObject)snes,1);
5350: PetscLogObjectParent((PetscObject)snes,(PetscObject)snes->ksp);
5352: KSPSetPreSolve(snes->ksp,(PetscErrorCode (*)(KSP,Vec,Vec,void*))KSPPreSolve_SNESEW,snes);
5353: KSPSetPostSolve(snes->ksp,(PetscErrorCode (*)(KSP,Vec,Vec,void*))KSPPostSolve_SNESEW,snes);
5355: PetscOptionsGetBool(((PetscObject)snes)->options,((PetscObject)snes)->prefix,"-ksp_monitor_snes",&monitor,NULL);
5356: if (monitor) {
5357: KSPMonitorSet(snes->ksp,KSPMonitorSNES,snes,NULL);
5358: }
5359: monitor = PETSC_FALSE;
5360: PetscOptionsGetBool(((PetscObject)snes)->options,((PetscObject)snes)->prefix,"-ksp_monitor_snes_lg",&monitor,NULL);
5361: if (monitor) {
5362: PetscObject *objs;
5363: KSPMonitorSNESLGResidualNormCreate(PetscObjectComm((PetscObject)snes),NULL,NULL,PETSC_DECIDE,PETSC_DECIDE,600,600,&objs);
5364: objs[0] = (PetscObject) snes;
5365: KSPMonitorSet(snes->ksp,(PetscErrorCode (*)(KSP,PetscInt,PetscReal,void*))KSPMonitorSNESLGResidualNorm,objs,(PetscErrorCode (*)(void**))KSPMonitorSNESLGResidualNormDestroy);
5366: }
5367: PetscObjectSetOptions((PetscObject)snes->ksp,((PetscObject)snes)->options);
5368: }
5369: *ksp = snes->ksp;
5370: return(0);
5371: }
5374: #include <petsc/private/dmimpl.h>
5375: /*@
5376: SNESSetDM - Sets the DM that may be used by some nonlinear solvers or their underlying preconditioners
5378: Logically Collective on SNES
5380: Input Parameters:
5381: + snes - the nonlinear solver context
5382: - dm - the dm, cannot be NULL
5384: Notes:
5385: A DM can only be used for solving one problem at a time because information about the problem is stored on the DM,
5386: even when not using interfaces like DMSNESSetFunction(). Use DMClone() to get a distinct DM when solving different
5387: problems using the same function space.
5389: Level: intermediate
5391: .seealso: SNESGetDM(), KSPSetDM(), KSPGetDM()
5392: @*/
5393: PetscErrorCode SNESSetDM(SNES snes,DM dm)
5394: {
5396: KSP ksp;
5397: DMSNES sdm;
5402: PetscObjectReference((PetscObject)dm);
5403: if (snes->dm) { /* Move the DMSNES context over to the new DM unless the new DM already has one */
5404: if (snes->dm->dmsnes && !dm->dmsnes) {
5405: DMCopyDMSNES(snes->dm,dm);
5406: DMGetDMSNES(snes->dm,&sdm);
5407: if (sdm->originaldm == snes->dm) sdm->originaldm = dm; /* Grant write privileges to the replacement DM */
5408: }
5409: DMCoarsenHookRemove(snes->dm,DMCoarsenHook_SNESVecSol,DMRestrictHook_SNESVecSol,snes);
5410: DMDestroy(&snes->dm);
5411: }
5412: snes->dm = dm;
5413: snes->dmAuto = PETSC_FALSE;
5415: SNESGetKSP(snes,&ksp);
5416: KSPSetDM(ksp,dm);
5417: KSPSetDMActive(ksp,PETSC_FALSE);
5418: if (snes->npc) {
5419: SNESSetDM(snes->npc, snes->dm);
5420: SNESSetNPCSide(snes,snes->npcside);
5421: }
5422: return(0);
5423: }
5425: /*@
5426: SNESGetDM - Gets the DM that may be used by some preconditioners
5428: Not Collective but DM obtained is parallel on SNES
5430: Input Parameter:
5431: . snes - the preconditioner context
5433: Output Parameter:
5434: . dm - the dm
5436: Level: intermediate
5438: .seealso: SNESSetDM(), KSPSetDM(), KSPGetDM()
5439: @*/
5440: PetscErrorCode SNESGetDM(SNES snes,DM *dm)
5441: {
5446: if (!snes->dm) {
5447: DMShellCreate(PetscObjectComm((PetscObject)snes),&snes->dm);
5448: snes->dmAuto = PETSC_TRUE;
5449: }
5450: *dm = snes->dm;
5451: return(0);
5452: }
5454: /*@
5455: SNESSetNPC - Sets the nonlinear preconditioner to be used.
5457: Collective on SNES
5459: Input Parameters:
5460: + snes - iterative context obtained from SNESCreate()
5461: - pc - the preconditioner object
5463: Notes:
5464: Use SNESGetNPC() to retrieve the preconditioner context (for example,
5465: to configure it using the API).
5467: Level: developer
5469: .keywords: SNES, set, precondition
5470: .seealso: SNESGetNPC(), SNESHasNPC()
5471: @*/
5472: PetscErrorCode SNESSetNPC(SNES snes, SNES pc)
5473: {
5480: PetscObjectReference((PetscObject) pc);
5481: SNESDestroy(&snes->npc);
5482: snes->npc = pc;
5483: PetscLogObjectParent((PetscObject)snes, (PetscObject)snes->npc);
5484: return(0);
5485: }
5487: /*@
5488: SNESGetNPC - Creates a nonlinear preconditioning solver (SNES) to be used to precondition the nonlinear solver.
5490: Not Collective; but any changes to the obtained SNES object must be applied collectively
5492: Input Parameter:
5493: . snes - iterative context obtained from SNESCreate()
5495: Output Parameter:
5496: . pc - preconditioner context
5498: Notes:
5499: If a SNES was previously set with SNESSetNPC() then that SNES is returned.
5501: The (preconditioner) SNES returned automatically inherits the same nonlinear function and Jacobian supplied to the original
5502: SNES during SNESSetUp()
5504: Level: developer
5506: .keywords: SNES, get, preconditioner
5507: .seealso: SNESSetNPC(), SNESHasNPC(), SNES, SNESCreate()
5508: @*/
5509: PetscErrorCode SNESGetNPC(SNES snes, SNES *pc)
5510: {
5512: const char *optionsprefix;
5517: if (!snes->npc) {
5518: SNESCreate(PetscObjectComm((PetscObject)snes),&snes->npc);
5519: PetscObjectIncrementTabLevel((PetscObject)snes->npc,(PetscObject)snes,1);
5520: PetscLogObjectParent((PetscObject)snes,(PetscObject)snes->npc);
5521: SNESGetOptionsPrefix(snes,&optionsprefix);
5522: SNESSetOptionsPrefix(snes->npc,optionsprefix);
5523: SNESAppendOptionsPrefix(snes->npc,"npc_");
5524: SNESSetCountersReset(snes->npc,PETSC_FALSE);
5525: }
5526: *pc = snes->npc;
5527: return(0);
5528: }
5530: /*@
5531: SNESHasNPC - Returns whether a nonlinear preconditioner exists
5533: Not Collective
5535: Input Parameter:
5536: . snes - iterative context obtained from SNESCreate()
5538: Output Parameter:
5539: . has_npc - whether the SNES has an NPC or not
5541: Level: developer
5543: .keywords: SNES, has, preconditioner
5544: .seealso: SNESSetNPC(), SNESGetNPC()
5545: @*/
5546: PetscErrorCode SNESHasNPC(SNES snes, PetscBool *has_npc)
5547: {
5550: *has_npc = (PetscBool) (snes->npc ? PETSC_TRUE : PETSC_FALSE);
5551: return(0);
5552: }
5554: /*@
5555: SNESSetNPCSide - Sets the preconditioning side.
5557: Logically Collective on SNES
5559: Input Parameter:
5560: . snes - iterative context obtained from SNESCreate()
5562: Output Parameter:
5563: . side - the preconditioning side, where side is one of
5564: .vb
5565: PC_LEFT - left preconditioning
5566: PC_RIGHT - right preconditioning (default for most nonlinear solvers)
5567: .ve
5569: Options Database Keys:
5570: . -snes_pc_side <right,left>
5572: Notes:
5573: SNESNRICHARDSON and SNESNCG only support left preconditioning.
5575: Level: intermediate
5577: .keywords: SNES, set, right, left, side, preconditioner, flag
5579: .seealso: SNESGetNPCSide(), KSPSetPCSide()
5580: @*/
5581: PetscErrorCode SNESSetNPCSide(SNES snes,PCSide side)
5582: {
5586: snes->npcside= side;
5587: return(0);
5588: }
5590: /*@
5591: SNESGetNPCSide - Gets the preconditioning side.
5593: Not Collective
5595: Input Parameter:
5596: . snes - iterative context obtained from SNESCreate()
5598: Output Parameter:
5599: . side - the preconditioning side, where side is one of
5600: .vb
5601: PC_LEFT - left preconditioning
5602: PC_RIGHT - right preconditioning (default for most nonlinear solvers)
5603: .ve
5605: Level: intermediate
5607: .keywords: SNES, get, right, left, side, preconditioner, flag
5609: .seealso: SNESSetNPCSide(), KSPGetPCSide()
5610: @*/
5611: PetscErrorCode SNESGetNPCSide(SNES snes,PCSide *side)
5612: {
5616: *side = snes->npcside;
5617: return(0);
5618: }
5620: /*@
5621: SNESSetLineSearch - Sets the linesearch on the SNES instance.
5623: Collective on SNES
5625: Input Parameters:
5626: + snes - iterative context obtained from SNESCreate()
5627: - linesearch - the linesearch object
5629: Notes:
5630: Use SNESGetLineSearch() to retrieve the preconditioner context (for example,
5631: to configure it using the API).
5633: Level: developer
5635: .keywords: SNES, set, linesearch
5636: .seealso: SNESGetLineSearch()
5637: @*/
5638: PetscErrorCode SNESSetLineSearch(SNES snes, SNESLineSearch linesearch)
5639: {
5646: PetscObjectReference((PetscObject) linesearch);
5647: SNESLineSearchDestroy(&snes->linesearch);
5649: snes->linesearch = linesearch;
5651: PetscLogObjectParent((PetscObject)snes, (PetscObject)snes->linesearch);
5652: return(0);
5653: }
5655: /*@
5656: SNESGetLineSearch - Returns a pointer to the line search context set with SNESSetLineSearch()
5657: or creates a default line search instance associated with the SNES and returns it.
5659: Not Collective
5661: Input Parameter:
5662: . snes - iterative context obtained from SNESCreate()
5664: Output Parameter:
5665: . linesearch - linesearch context
5667: Level: beginner
5669: .keywords: SNES, get, linesearch
5670: .seealso: SNESSetLineSearch(), SNESLineSearchCreate()
5671: @*/
5672: PetscErrorCode SNESGetLineSearch(SNES snes, SNESLineSearch *linesearch)
5673: {
5675: const char *optionsprefix;
5680: if (!snes->linesearch) {
5681: SNESGetOptionsPrefix(snes, &optionsprefix);
5682: SNESLineSearchCreate(PetscObjectComm((PetscObject)snes), &snes->linesearch);
5683: SNESLineSearchSetSNES(snes->linesearch, snes);
5684: SNESLineSearchAppendOptionsPrefix(snes->linesearch, optionsprefix);
5685: PetscObjectIncrementTabLevel((PetscObject) snes->linesearch, (PetscObject) snes, 1);
5686: PetscLogObjectParent((PetscObject)snes, (PetscObject)snes->linesearch);
5687: }
5688: *linesearch = snes->linesearch;
5689: return(0);
5690: }
5692: #if defined(PETSC_HAVE_MATLAB_ENGINE)
5693: #include <mex.h>
5695: typedef struct {char *funcname; mxArray *ctx;} SNESMatlabContext;
5697: /*
5698: SNESComputeFunction_Matlab - Calls the function that has been set with SNESSetFunctionMatlab().
5700: Collective on SNES
5702: Input Parameters:
5703: + snes - the SNES context
5704: - x - input vector
5706: Output Parameter:
5707: . y - function vector, as set by SNESSetFunction()
5709: Notes:
5710: SNESComputeFunction() is typically used within nonlinear solvers
5711: implementations, so most users would not generally call this routine
5712: themselves.
5714: Level: developer
5716: .keywords: SNES, nonlinear, compute, function
5718: .seealso: SNESSetFunction(), SNESGetFunction()
5719: */
5720: PetscErrorCode SNESComputeFunction_Matlab(SNES snes,Vec x,Vec y, void *ctx)
5721: {
5722: PetscErrorCode ierr;
5723: SNESMatlabContext *sctx = (SNESMatlabContext*)ctx;
5724: int nlhs = 1,nrhs = 5;
5725: mxArray *plhs[1],*prhs[5];
5726: long long int lx = 0,ly = 0,ls = 0;
5735: /* call Matlab function in ctx with arguments x and y */
5737: PetscMemcpy(&ls,&snes,sizeof(snes));
5738: PetscMemcpy(&lx,&x,sizeof(x));
5739: PetscMemcpy(&ly,&y,sizeof(x));
5740: prhs[0] = mxCreateDoubleScalar((double)ls);
5741: prhs[1] = mxCreateDoubleScalar((double)lx);
5742: prhs[2] = mxCreateDoubleScalar((double)ly);
5743: prhs[3] = mxCreateString(sctx->funcname);
5744: prhs[4] = sctx->ctx;
5745: mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscSNESComputeFunctionInternal");
5746: mxGetScalar(plhs[0]);
5747: mxDestroyArray(prhs[0]);
5748: mxDestroyArray(prhs[1]);
5749: mxDestroyArray(prhs[2]);
5750: mxDestroyArray(prhs[3]);
5751: mxDestroyArray(plhs[0]);
5752: return(0);
5753: }
5755: /*
5756: SNESSetFunctionMatlab - Sets the function evaluation routine and function
5757: vector for use by the SNES routines in solving systems of nonlinear
5758: equations from MATLAB. Here the function is a string containing the name of a MATLAB function
5760: Logically Collective on SNES
5762: Input Parameters:
5763: + snes - the SNES context
5764: . r - vector to store function value
5765: - f - function evaluation routine
5767: Notes:
5768: The Newton-like methods typically solve linear systems of the form
5769: $ f'(x) x = -f(x),
5770: where f'(x) denotes the Jacobian matrix and f(x) is the function.
5772: Level: beginner
5774: Developer Note: This bleeds the allocated memory SNESMatlabContext *sctx;
5776: .keywords: SNES, nonlinear, set, function
5778: .seealso: SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESSetFunction()
5779: */
5780: PetscErrorCode SNESSetFunctionMatlab(SNES snes,Vec r,const char *f,mxArray *ctx)
5781: {
5782: PetscErrorCode ierr;
5783: SNESMatlabContext *sctx;
5786: /* currently sctx is memory bleed */
5787: PetscNew(&sctx);
5788: PetscStrallocpy(f,&sctx->funcname);
5789: /*
5790: This should work, but it doesn't
5791: sctx->ctx = ctx;
5792: mexMakeArrayPersistent(sctx->ctx);
5793: */
5794: sctx->ctx = mxDuplicateArray(ctx);
5795: SNESSetFunction(snes,r,SNESComputeFunction_Matlab,sctx);
5796: return(0);
5797: }
5799: /*
5800: SNESComputeJacobian_Matlab - Calls the function that has been set with SNESSetJacobianMatlab().
5802: Collective on SNES
5804: Input Parameters:
5805: + snes - the SNES context
5806: . x - input vector
5807: . A, B - the matrices
5808: - ctx - user context
5810: Level: developer
5812: .keywords: SNES, nonlinear, compute, function
5814: .seealso: SNESSetFunction(), SNESGetFunction()
5815: @*/
5816: PetscErrorCode SNESComputeJacobian_Matlab(SNES snes,Vec x,Mat A,Mat B,void *ctx)
5817: {
5818: PetscErrorCode ierr;
5819: SNESMatlabContext *sctx = (SNESMatlabContext*)ctx;
5820: int nlhs = 2,nrhs = 6;
5821: mxArray *plhs[2],*prhs[6];
5822: long long int lx = 0,lA = 0,ls = 0, lB = 0;
5828: /* call Matlab function in ctx with arguments x and y */
5830: PetscMemcpy(&ls,&snes,sizeof(snes));
5831: PetscMemcpy(&lx,&x,sizeof(x));
5832: PetscMemcpy(&lA,A,sizeof(x));
5833: PetscMemcpy(&lB,B,sizeof(x));
5834: prhs[0] = mxCreateDoubleScalar((double)ls);
5835: prhs[1] = mxCreateDoubleScalar((double)lx);
5836: prhs[2] = mxCreateDoubleScalar((double)lA);
5837: prhs[3] = mxCreateDoubleScalar((double)lB);
5838: prhs[4] = mxCreateString(sctx->funcname);
5839: prhs[5] = sctx->ctx;
5840: mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscSNESComputeJacobianInternal");
5841: mxGetScalar(plhs[0]);
5842: mxDestroyArray(prhs[0]);
5843: mxDestroyArray(prhs[1]);
5844: mxDestroyArray(prhs[2]);
5845: mxDestroyArray(prhs[3]);
5846: mxDestroyArray(prhs[4]);
5847: mxDestroyArray(plhs[0]);
5848: mxDestroyArray(plhs[1]);
5849: return(0);
5850: }
5852: /*
5853: SNESSetJacobianMatlab - Sets the Jacobian function evaluation routine and two empty Jacobian matrices
5854: vector for use by the SNES routines in solving systems of nonlinear
5855: equations from MATLAB. Here the function is a string containing the name of a MATLAB function
5857: Logically Collective on SNES
5859: Input Parameters:
5860: + snes - the SNES context
5861: . A,B - Jacobian matrices
5862: . J - function evaluation routine
5863: - ctx - user context
5865: Level: developer
5867: Developer Note: This bleeds the allocated memory SNESMatlabContext *sctx;
5869: .keywords: SNES, nonlinear, set, function
5871: .seealso: SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESSetFunction(), J
5872: */
5873: PetscErrorCode SNESSetJacobianMatlab(SNES snes,Mat A,Mat B,const char *J,mxArray *ctx)
5874: {
5875: PetscErrorCode ierr;
5876: SNESMatlabContext *sctx;
5879: /* currently sctx is memory bleed */
5880: PetscNew(&sctx);
5881: PetscStrallocpy(J,&sctx->funcname);
5882: /*
5883: This should work, but it doesn't
5884: sctx->ctx = ctx;
5885: mexMakeArrayPersistent(sctx->ctx);
5886: */
5887: sctx->ctx = mxDuplicateArray(ctx);
5888: SNESSetJacobian(snes,A,B,SNESComputeJacobian_Matlab,sctx);
5889: return(0);
5890: }
5892: /*
5893: SNESMonitor_Matlab - Calls the function that has been set with SNESMonitorSetMatlab().
5895: Collective on SNES
5897: .seealso: SNESSetFunction(), SNESGetFunction()
5898: @*/
5899: PetscErrorCode SNESMonitor_Matlab(SNES snes,PetscInt it, PetscReal fnorm, void *ctx)
5900: {
5901: PetscErrorCode ierr;
5902: SNESMatlabContext *sctx = (SNESMatlabContext*)ctx;
5903: int nlhs = 1,nrhs = 6;
5904: mxArray *plhs[1],*prhs[6];
5905: long long int lx = 0,ls = 0;
5906: Vec x = snes->vec_sol;
5911: PetscMemcpy(&ls,&snes,sizeof(snes));
5912: PetscMemcpy(&lx,&x,sizeof(x));
5913: prhs[0] = mxCreateDoubleScalar((double)ls);
5914: prhs[1] = mxCreateDoubleScalar((double)it);
5915: prhs[2] = mxCreateDoubleScalar((double)fnorm);
5916: prhs[3] = mxCreateDoubleScalar((double)lx);
5917: prhs[4] = mxCreateString(sctx->funcname);
5918: prhs[5] = sctx->ctx;
5919: mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscSNESMonitorInternal");
5920: mxGetScalar(plhs[0]);
5921: mxDestroyArray(prhs[0]);
5922: mxDestroyArray(prhs[1]);
5923: mxDestroyArray(prhs[2]);
5924: mxDestroyArray(prhs[3]);
5925: mxDestroyArray(prhs[4]);
5926: mxDestroyArray(plhs[0]);
5927: return(0);
5928: }
5930: /*
5931: SNESMonitorSetMatlab - Sets the monitor function from MATLAB
5933: Level: developer
5935: Developer Note: This bleeds the allocated memory SNESMatlabContext *sctx;
5937: .keywords: SNES, nonlinear, set, function
5939: .seealso: SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESSetFunction()
5940: */
5941: PetscErrorCode SNESMonitorSetMatlab(SNES snes,const char *f,mxArray *ctx)
5942: {
5943: PetscErrorCode ierr;
5944: SNESMatlabContext *sctx;
5947: /* currently sctx is memory bleed */
5948: PetscNew(&sctx);
5949: PetscStrallocpy(f,&sctx->funcname);
5950: /*
5951: This should work, but it doesn't
5952: sctx->ctx = ctx;
5953: mexMakeArrayPersistent(sctx->ctx);
5954: */
5955: sctx->ctx = mxDuplicateArray(ctx);
5956: SNESMonitorSet(snes,SNESMonitor_Matlab,sctx,NULL);
5957: return(0);
5958: }
5960: #endif