Actual source code: snes.c
petsc-3.4.5 2014-06-29
2: #include <petsc-private/snesimpl.h> /*I "petscsnes.h" I*/
3: #include <petscdmshell.h>
5: PetscBool SNESRegisterAllCalled = PETSC_FALSE;
6: PetscFunctionList SNESList = NULL;
8: /* Logging support */
9: PetscClassId SNES_CLASSID, DMSNES_CLASSID;
10: PetscLogEvent SNES_Solve, SNES_FunctionEval, SNES_JacobianEval, SNES_GSEval, SNES_NPCSolve;
14: /*@
15: SNESSetErrorIfNotConverged - Causes SNESSolve() to generate an error if the solver has not converged.
17: Logically Collective on SNES
19: Input Parameters:
20: + snes - iterative context obtained from SNESCreate()
21: - flg - PETSC_TRUE indicates you want the error generated
23: Options database keys:
24: . -snes_error_if_not_converged : this takes an optional truth value (0/1/no/yes/true/false)
26: Level: intermediate
28: Notes:
29: Normally PETSc continues if a linear solver fails to converge, you can call SNESGetConvergedReason() after a SNESSolve()
30: to determine if it has converged.
32: .keywords: SNES, set, initial guess, nonzero
34: .seealso: SNESGetErrorIfNotConverged(), KSPGetErrorIfNotConverged(), KSPSetErrorIFNotConverged()
35: @*/
36: PetscErrorCode SNESSetErrorIfNotConverged(SNES snes,PetscBool flg)
37: {
41: snes->errorifnotconverged = flg;
42: return(0);
43: }
47: /*@
48: SNESGetErrorIfNotConverged - Will SNESSolve() generate an error if the solver does not converge?
50: Not Collective
52: Input Parameter:
53: . snes - iterative context obtained from SNESCreate()
55: Output Parameter:
56: . flag - PETSC_TRUE if it will generate an error, else PETSC_FALSE
58: Level: intermediate
60: .keywords: SNES, set, initial guess, nonzero
62: .seealso: SNESSetErrorIfNotConverged(), KSPGetErrorIfNotConverged(), KSPSetErrorIFNotConverged()
63: @*/
64: PetscErrorCode SNESGetErrorIfNotConverged(SNES snes,PetscBool *flag)
65: {
69: *flag = snes->errorifnotconverged;
70: return(0);
71: }
75: /*@
76: SNESSetFunctionDomainError - tells SNES that the input vector to your SNESFunction is not
77: in the functions domain. For example, negative pressure.
79: Logically Collective on SNES
81: Input Parameters:
82: . snes - the SNES context
84: Level: advanced
86: .keywords: SNES, view
88: .seealso: SNESCreate(), SNESSetFunction(), SNESFunction
89: @*/
90: PetscErrorCode SNESSetFunctionDomainError(SNES snes)
91: {
94: snes->domainerror = PETSC_TRUE;
95: return(0);
96: }
100: /*@
101: SNESGetFunctionDomainError - Gets the status of the domain error after a call to SNESComputeFunction;
103: Logically Collective on SNES
105: Input Parameters:
106: . snes - the SNES context
108: Output Parameters:
109: . domainerror - Set to PETSC_TRUE if there's a domain error; PETSC_FALSE otherwise.
111: Level: advanced
113: .keywords: SNES, view
115: .seealso: SNESSetFunctionDomainError(), SNESComputeFunction()
116: @*/
117: PetscErrorCode SNESGetFunctionDomainError(SNES snes, PetscBool *domainerror)
118: {
122: *domainerror = snes->domainerror;
123: return(0);
124: }
128: /*@C
129: SNESLoad - Loads a SNES that has been stored in binary with SNESView().
131: Collective on PetscViewer
133: Input Parameters:
134: + newdm - the newly loaded SNES, this needs to have been created with SNESCreate() or
135: some related function before a call to SNESLoad().
136: - viewer - binary file viewer, obtained from PetscViewerBinaryOpen()
138: Level: intermediate
140: Notes:
141: The type is determined by the data in the file, any type set into the SNES before this call is ignored.
143: Notes for advanced users:
144: Most users should not need to know the details of the binary storage
145: format, since SNESLoad() and TSView() completely hide these details.
146: But for anyone who's interested, the standard binary matrix storage
147: format is
148: .vb
149: has not yet been determined
150: .ve
152: .seealso: PetscViewerBinaryOpen(), SNESView(), MatLoad(), VecLoad()
153: @*/
154: PetscErrorCode SNESLoad(SNES snes, PetscViewer viewer)
155: {
157: PetscBool isbinary;
158: PetscInt classid;
159: char type[256];
160: KSP ksp;
161: DM dm;
162: DMSNES dmsnes;
167: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
168: if (!isbinary) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Invalid viewer; open viewer with PetscViewerBinaryOpen()");
170: PetscViewerBinaryRead(viewer,&classid,1,PETSC_INT);
171: if (classid != SNES_FILE_CLASSID) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONG,"Not SNES next in file");
172: PetscViewerBinaryRead(viewer,type,256,PETSC_CHAR);
173: SNESSetType(snes, type);
174: if (snes->ops->load) {
175: (*snes->ops->load)(snes,viewer);
176: }
177: SNESGetDM(snes,&dm);
178: DMGetDMSNES(dm,&dmsnes);
179: DMSNESLoad(dmsnes,viewer);
180: SNESGetKSP(snes,&ksp);
181: KSPLoad(ksp,viewer);
182: return(0);
183: }
185: #include <petscdraw.h>
186: #if defined(PETSC_HAVE_AMS)
187: #include <petscviewerams.h>
188: #endif
191: /*@C
192: SNESView - Prints the SNES data structure.
194: Collective on SNES
196: Input Parameters:
197: + SNES - the SNES context
198: - viewer - visualization context
200: Options Database Key:
201: . -snes_view - Calls SNESView() at end of SNESSolve()
203: Notes:
204: The available visualization contexts include
205: + PETSC_VIEWER_STDOUT_SELF - standard output (default)
206: - PETSC_VIEWER_STDOUT_WORLD - synchronized standard
207: output where only the first processor opens
208: the file. All other processors send their
209: data to the first processor to print.
211: The user can open an alternative visualization context with
212: PetscViewerASCIIOpen() - output to a specified file.
214: Level: beginner
216: .keywords: SNES, view
218: .seealso: PetscViewerASCIIOpen()
219: @*/
220: PetscErrorCode SNESView(SNES snes,PetscViewer viewer)
221: {
222: SNESKSPEW *kctx;
224: KSP ksp;
225: SNESLineSearch linesearch;
226: PetscBool iascii,isstring,isbinary,isdraw;
227: DMSNES dmsnes;
228: #if defined(PETSC_HAVE_AMS)
229: PetscBool isams;
230: #endif
234: if (!viewer) {
235: PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)snes),&viewer);
236: }
240: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
241: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSTRING,&isstring);
242: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
243: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
244: #if defined(PETSC_HAVE_AMS)
245: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERAMS,&isams);
246: #endif
247: if (iascii) {
248: PetscObjectPrintClassNamePrefixType((PetscObject)snes,viewer,"SNES Object");
249: if (!snes->setupcalled) {
250: PetscViewerASCIIPrintf(viewer," SNES has not been set up so information may be incomplete\n");
251: }
252: if (snes->ops->view) {
253: PetscViewerASCIIPushTab(viewer);
254: (*snes->ops->view)(snes,viewer);
255: PetscViewerASCIIPopTab(viewer);
256: }
257: PetscViewerASCIIPrintf(viewer," maximum iterations=%D, maximum function evaluations=%D\n",snes->max_its,snes->max_funcs);
258: PetscViewerASCIIPrintf(viewer," tolerances: relative=%G, absolute=%G, solution=%G\n",snes->rtol,snes->abstol,snes->stol);
259: PetscViewerASCIIPrintf(viewer," total number of linear solver iterations=%D\n",snes->linear_its);
260: PetscViewerASCIIPrintf(viewer," total number of function evaluations=%D\n",snes->nfuncs);
261: if (snes->gridsequence) {
262: PetscViewerASCIIPrintf(viewer," total number of grid sequence refinements=%D\n",snes->gridsequence);
263: }
264: if (snes->ksp_ewconv) {
265: kctx = (SNESKSPEW*)snes->kspconvctx;
266: if (kctx) {
267: PetscViewerASCIIPrintf(viewer," Eisenstat-Walker computation of KSP relative tolerance (version %D)\n",kctx->version);
268: PetscViewerASCIIPrintf(viewer," rtol_0=%G, rtol_max=%G, threshold=%G\n",kctx->rtol_0,kctx->rtol_max,kctx->threshold);
269: PetscViewerASCIIPrintf(viewer," gamma=%G, alpha=%G, alpha2=%G\n",kctx->gamma,kctx->alpha,kctx->alpha2);
270: }
271: }
272: if (snes->lagpreconditioner == -1) {
273: PetscViewerASCIIPrintf(viewer," Preconditioned is never rebuilt\n");
274: } else if (snes->lagpreconditioner > 1) {
275: PetscViewerASCIIPrintf(viewer," Preconditioned is rebuilt every %D new Jacobians\n",snes->lagpreconditioner);
276: }
277: if (snes->lagjacobian == -1) {
278: PetscViewerASCIIPrintf(viewer," Jacobian is never rebuilt\n");
279: } else if (snes->lagjacobian > 1) {
280: PetscViewerASCIIPrintf(viewer," Jacobian is rebuilt every %D SNES iterations\n",snes->lagjacobian);
281: }
282: } else if (isstring) {
283: const char *type;
284: SNESGetType(snes,&type);
285: PetscViewerStringSPrintf(viewer," %-3.3s",type);
286: } else if (isbinary) {
287: PetscInt classid = SNES_FILE_CLASSID;
288: MPI_Comm comm;
289: PetscMPIInt rank;
290: char type[256];
292: PetscObjectGetComm((PetscObject)snes,&comm);
293: MPI_Comm_rank(comm,&rank);
294: if (!rank) {
295: PetscViewerBinaryWrite(viewer,&classid,1,PETSC_INT,PETSC_FALSE);
296: PetscStrncpy(type,((PetscObject)snes)->type_name,256);
297: PetscViewerBinaryWrite(viewer,type,256,PETSC_CHAR,PETSC_FALSE);
298: }
299: if (snes->ops->view) {
300: (*snes->ops->view)(snes,viewer);
301: }
302: } else if (isdraw) {
303: PetscDraw draw;
304: char str[36];
305: PetscReal x,y,bottom,h;
307: PetscViewerDrawGetDraw(viewer,0,&draw);
308: PetscDrawGetCurrentPoint(draw,&x,&y);
309: PetscStrcpy(str,"SNES: ");
310: PetscStrcat(str,((PetscObject)snes)->type_name);
311: PetscDrawBoxedString(draw,x,y,PETSC_DRAW_BLUE,PETSC_DRAW_BLACK,str,NULL,&h);
312: bottom = y - h;
313: PetscDrawPushCurrentPoint(draw,x,bottom);
314: if (snes->ops->view) {
315: (*snes->ops->view)(snes,viewer);
316: }
317: #if defined(PETSC_HAVE_AMS)
318: } else if (isams) {
319: if (((PetscObject)snes)->amsmem == -1) {
320: PetscObjectViewAMS((PetscObject)snes,viewer);
321: PetscStackCallAMS(AMS_Memory_take_access,(((PetscObject)snes)->amsmem));
322: PetscStackCallAMS(AMS_Memory_add_field,(((PetscObject)snes)->amsmem,"its",&snes->iter,1,AMS_INT,AMS_READ,AMS_COMMON,AMS_REDUCT_UNDEF));
323: if (!snes->conv_hist) {
324: SNESSetConvergenceHistory(snes,NULL,NULL,PETSC_DECIDE,PETSC_FALSE);
325: }
326: PetscStackCallAMS(AMS_Memory_add_field,(((PetscObject)snes)->amsmem,"conv_hist",snes->conv_hist,10,AMS_DOUBLE,AMS_READ,AMS_COMMON,AMS_REDUCT_UNDEF));
327: PetscStackCallAMS(AMS_Memory_grant_access,(((PetscObject)snes)->amsmem));
328: }
329: #endif
330: }
331: if (snes->linesearch) {
332: PetscViewerASCIIPushTab(viewer);
333: SNESGetLineSearch(snes, &linesearch);
334: SNESLineSearchView(linesearch, viewer);
335: PetscViewerASCIIPopTab(viewer);
336: }
337: if (snes->pc && snes->usespc) {
338: PetscViewerASCIIPushTab(viewer);
339: SNESView(snes->pc, viewer);
340: PetscViewerASCIIPopTab(viewer);
341: }
342: PetscViewerASCIIPushTab(viewer);
343: DMGetDMSNES(snes->dm,&dmsnes);
344: DMSNESView(dmsnes, viewer);
345: PetscViewerASCIIPopTab(viewer);
346: if (snes->usesksp) {
347: SNESGetKSP(snes,&ksp);
348: PetscViewerASCIIPushTab(viewer);
349: KSPView(ksp,viewer);
350: PetscViewerASCIIPopTab(viewer);
351: }
352: if (isdraw) {
353: PetscDraw draw;
354: PetscViewerDrawGetDraw(viewer,0,&draw);
355: PetscDrawPopCurrentPoint(draw);
356: }
357: return(0);
358: }
360: /*
361: We retain a list of functions that also take SNES command
362: line options. These are called at the end SNESSetFromOptions()
363: */
364: #define MAXSETFROMOPTIONS 5
365: static PetscInt numberofsetfromoptions;
366: static PetscErrorCode (*othersetfromoptions[MAXSETFROMOPTIONS])(SNES);
370: /*@C
371: SNESAddOptionsChecker - Adds an additional function to check for SNES options.
373: Not Collective
375: Input Parameter:
376: . snescheck - function that checks for options
378: Level: developer
380: .seealso: SNESSetFromOptions()
381: @*/
382: PetscErrorCode SNESAddOptionsChecker(PetscErrorCode (*snescheck)(SNES))
383: {
385: if (numberofsetfromoptions >= MAXSETFROMOPTIONS) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE, "Too many options checkers, only %D allowed", MAXSETFROMOPTIONS);
386: othersetfromoptions[numberofsetfromoptions++] = snescheck;
387: return(0);
388: }
390: extern PetscErrorCode SNESDefaultMatrixFreeCreate2(SNES,Vec,Mat*);
394: static PetscErrorCode SNESSetUpMatrixFree_Private(SNES snes, PetscBool hasOperator, PetscInt version)
395: {
396: Mat J;
397: KSP ksp;
398: PC pc;
399: PetscBool match;
405: if (!snes->vec_func && (snes->jacobian || snes->jacobian_pre)) {
406: Mat A = snes->jacobian, B = snes->jacobian_pre;
407: MatGetVecs(A ? A : B, NULL,&snes->vec_func);
408: }
410: if (version == 1) {
411: MatCreateSNESMF(snes,&J);
412: MatMFFDSetOptionsPrefix(J,((PetscObject)snes)->prefix);
413: MatSetFromOptions(J);
414: } else if (version == 2) {
415: if (!snes->vec_func) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"SNESSetFunction() must be called first");
416: #if !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_REAL_SINGLE) && !defined(PETSC_USE_REAL___FLOAT128)
417: SNESDefaultMatrixFreeCreate2(snes,snes->vec_func,&J);
418: #else
419: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP, "matrix-free operator rutines (version 2)");
420: #endif
421: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE, "matrix-free operator rutines, only version 1 and 2");
423: PetscInfo1(snes,"Setting default matrix-free operator routines (version %D)\n", version);
424: if (hasOperator) {
426: /* This version replaces the user provided Jacobian matrix with a
427: matrix-free version but still employs the user-provided preconditioner matrix. */
428: SNESSetJacobian(snes,J,0,0,0);
429: } else {
430: /* This version replaces both the user-provided Jacobian and the user-
431: provided preconditioner Jacobian with the default matrix free version. */
432: if ((snes->pcside == PC_LEFT) && snes->pc) {
433: if (!snes->jacobian){SNESSetJacobian(snes,J,0,0,0);}
434: } else {
435: SNESSetJacobian(snes,J,J,MatMFFDComputeJacobian,0);
436: }
437: /* Force no preconditioner */
438: SNESGetKSP(snes,&ksp);
439: KSPGetPC(ksp,&pc);
440: PetscObjectTypeCompare((PetscObject)pc,PCSHELL,&match);
441: if (!match) {
442: PetscInfo(snes,"Setting default matrix-free preconditioner routines\nThat is no preconditioner is being used\n");
443: PCSetType(pc,PCNONE);
444: }
445: }
446: MatDestroy(&J);
447: return(0);
448: }
452: static PetscErrorCode DMRestrictHook_SNESVecSol(DM dmfine,Mat Restrict,Vec Rscale,Mat Inject,DM dmcoarse,void *ctx)
453: {
454: SNES snes = (SNES)ctx;
456: Vec Xfine,Xfine_named = NULL,Xcoarse;
459: if (PetscLogPrintInfo) {
460: PetscInt finelevel,coarselevel,fineclevel,coarseclevel;
461: DMGetRefineLevel(dmfine,&finelevel);
462: DMGetCoarsenLevel(dmfine,&fineclevel);
463: DMGetRefineLevel(dmcoarse,&coarselevel);
464: DMGetCoarsenLevel(dmcoarse,&coarseclevel);
465: PetscInfo4(dmfine,"Restricting SNES solution vector from level %D-%D to level %D-%D\n",finelevel,fineclevel,coarselevel,coarseclevel);
466: }
467: if (dmfine == snes->dm) Xfine = snes->vec_sol;
468: else {
469: DMGetNamedGlobalVector(dmfine,"SNESVecSol",&Xfine_named);
470: Xfine = Xfine_named;
471: }
472: DMGetNamedGlobalVector(dmcoarse,"SNESVecSol",&Xcoarse);
473: MatRestrict(Restrict,Xfine,Xcoarse);
474: VecPointwiseMult(Xcoarse,Xcoarse,Rscale);
475: DMRestoreNamedGlobalVector(dmcoarse,"SNESVecSol",&Xcoarse);
476: if (Xfine_named) {DMRestoreNamedGlobalVector(dmfine,"SNESVecSol",&Xfine_named);}
477: return(0);
478: }
482: static PetscErrorCode DMCoarsenHook_SNESVecSol(DM dm,DM dmc,void *ctx)
483: {
487: DMCoarsenHookAdd(dmc,DMCoarsenHook_SNESVecSol,DMRestrictHook_SNESVecSol,ctx);
488: return(0);
489: }
493: /* This may be called to rediscretize the operator on levels of linear multigrid. The DM shuffle is so the user can
494: * safely call SNESGetDM() in their residual evaluation routine. */
495: static PetscErrorCode KSPComputeOperators_SNES(KSP ksp,Mat A,Mat B,MatStructure *mstruct,void *ctx)
496: {
497: SNES snes = (SNES)ctx;
499: Mat Asave = A,Bsave = B;
500: Vec X,Xnamed = NULL;
501: DM dmsave;
502: void *ctxsave;
503: PetscErrorCode (*jac)(SNES,Vec,Mat*,Mat*,MatStructure*,void*);
506: dmsave = snes->dm;
507: KSPGetDM(ksp,&snes->dm);
508: if (dmsave == snes->dm) X = snes->vec_sol; /* We are on the finest level */
509: else { /* We are on a coarser level, this vec was initialized using a DM restrict hook */
510: DMGetNamedGlobalVector(snes->dm,"SNESVecSol",&Xnamed);
511: X = Xnamed;
512: SNESGetJacobian(snes,NULL,NULL,&jac,&ctxsave);
513: /* If the DM's don't match up, the MatFDColoring context needed for the jacobian won't match up either -- fixit. */
514: if (jac == SNESComputeJacobianDefaultColor) {
515: SNESSetJacobian(snes,NULL,NULL,SNESComputeJacobianDefaultColor,0);
516: }
517: }
518: /* put the previous context back */
520: SNESComputeJacobian(snes,X,&A,&B,mstruct);
521: if (snes->dm != dmsave && jac == SNESComputeJacobianDefaultColor) {
522: SNESSetJacobian(snes,NULL,NULL,jac,ctxsave);
523: }
525: if (A != Asave || B != Bsave) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_SUP,"No support for changing matrices at this time");
526: if (Xnamed) {
527: DMRestoreNamedGlobalVector(snes->dm,"SNESVecSol",&Xnamed);
528: }
529: snes->dm = dmsave;
530: return(0);
531: }
535: /*@
536: SNESSetUpMatrices - ensures that matrices are available for SNES, to be called by SNESSetUp_XXX()
538: Collective
540: Input Arguments:
541: . snes - snes to configure
543: Level: developer
545: .seealso: SNESSetUp()
546: @*/
547: PetscErrorCode SNESSetUpMatrices(SNES snes)
548: {
550: DM dm;
551: DMSNES sdm;
554: SNESGetDM(snes,&dm);
555: DMGetDMSNES(dm,&sdm);
556: if (!sdm->ops->computejacobian) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_PLIB,"DMSNES not properly configured");
557: else if (!snes->jacobian && snes->mf) {
558: Mat J;
559: void *functx;
560: MatCreateSNESMF(snes,&J);
561: MatMFFDSetOptionsPrefix(J,((PetscObject)snes)->prefix);
562: MatSetFromOptions(J);
563: SNESGetFunction(snes,NULL,NULL,&functx);
564: SNESSetJacobian(snes,J,J,0,0);
565: MatDestroy(&J);
566: } else if (snes->mf_operator && !snes->jacobian_pre && !snes->jacobian) {
567: Mat J,B;
568: MatCreateSNESMF(snes,&J);
569: MatMFFDSetOptionsPrefix(J,((PetscObject)snes)->prefix);
570: MatSetFromOptions(J);
571: DMCreateMatrix(snes->dm,MATAIJ,&B);
572: /* sdm->computejacobian was already set to reach here */
573: SNESSetJacobian(snes,J,B,NULL,NULL);
574: MatDestroy(&J);
575: MatDestroy(&B);
576: } else if (!snes->jacobian_pre) {
577: Mat J,B;
578: J = snes->jacobian;
579: DMCreateMatrix(snes->dm,MATAIJ,&B);
580: SNESSetJacobian(snes,J ? J : B,B,NULL,NULL);
581: MatDestroy(&B);
582: }
583: {
584: KSP ksp;
585: SNESGetKSP(snes,&ksp);
586: KSPSetComputeOperators(ksp,KSPComputeOperators_SNES,snes);
587: DMCoarsenHookAdd(snes->dm,DMCoarsenHook_SNESVecSol,DMRestrictHook_SNESVecSol,snes);
588: }
589: return(0);
590: }
594: /*@
595: SNESSetFromOptions - Sets various SNES and KSP parameters from user options.
597: Collective on SNES
599: Input Parameter:
600: . snes - the SNES context
602: Options Database Keys:
603: + -snes_type <type> - ls, tr, ngmres, ncg, richardson, qn, vi, fas
604: . -snes_stol - convergence tolerance in terms of the norm
605: of the change in the solution between steps
606: . -snes_atol <abstol> - absolute tolerance of residual norm
607: . -snes_rtol <rtol> - relative decrease in tolerance norm from initial
608: . -snes_max_it <max_it> - maximum number of iterations
609: . -snes_max_funcs <max_funcs> - maximum number of function evaluations
610: . -snes_max_fail <max_fail> - maximum number of line search failures allowed before stopping, default is none
611: . -snes_max_linear_solve_fail - number of linear solver failures before SNESSolve() stops
612: . -snes_lag_preconditioner <lag> - how often preconditioner is rebuilt (use -1 to never rebuild)
613: . -snes_lag_jacobian <lag> - how often Jacobian is rebuilt (use -1 to never rebuild)
614: . -snes_trtol <trtol> - trust region tolerance
615: . -snes_no_convergence_test - skip convergence test in nonlinear
616: solver; hence iterations will continue until max_it
617: or some other criterion is reached. Saves expense
618: of convergence test
619: . -snes_monitor <optional filename> - prints residual norm at each iteration. if no
620: filename given prints to stdout
621: . -snes_monitor_solution - plots solution at each iteration
622: . -snes_monitor_residual - plots residual (not its norm) at each iteration
623: . -snes_monitor_solution_update - plots update to solution at each iteration
624: . -snes_monitor_lg_residualnorm - plots residual norm at each iteration
625: . -snes_monitor_lg_range - plots residual norm at each iteration
626: . -snes_fd - use finite differences to compute Jacobian; very slow, only for testing
627: . -snes_fd_color - use finite differences with coloring to compute Jacobian
628: . -snes_mf_ksp_monitor - if using matrix-free multiply then print h at each KSP iteration
629: - -snes_converged_reason - print the reason for convergence/divergence after each solve
631: Options Database for Eisenstat-Walker method:
632: + -snes_ksp_ew - use Eisenstat-Walker method for determining linear system convergence
633: . -snes_ksp_ew_version ver - version of Eisenstat-Walker method
634: . -snes_ksp_ew_rtol0 <rtol0> - Sets rtol0
635: . -snes_ksp_ew_rtolmax <rtolmax> - Sets rtolmax
636: . -snes_ksp_ew_gamma <gamma> - Sets gamma
637: . -snes_ksp_ew_alpha <alpha> - Sets alpha
638: . -snes_ksp_ew_alpha2 <alpha2> - Sets alpha2
639: - -snes_ksp_ew_threshold <threshold> - Sets threshold
641: Notes:
642: To see all options, run your program with the -help option or consult
643: the <A href="../../docs/manual.pdf#nameddest=ch_snes">SNES chapter of the users manual</A>.
645: Level: beginner
647: .keywords: SNES, nonlinear, set, options, database
649: .seealso: SNESSetOptionsPrefix()
650: @*/
651: PetscErrorCode SNESSetFromOptions(SNES snes)
652: {
653: PetscBool flg,pcset;
654: PetscInt i,indx,lag,grids;
655: MatStructure matflag;
656: const char *deft = SNESNEWTONLS;
657: const char *convtests[] = {"default","skip"};
658: SNESKSPEW *kctx = NULL;
659: char type[256], monfilename[PETSC_MAX_PATH_LEN];
660: PetscViewer monviewer;
662: PCSide pcside;
663: const char *optionsprefix;
667: if (!SNESRegisterAllCalled) {SNESRegisterAll();}
668: PetscObjectOptionsBegin((PetscObject)snes);
669: if (((PetscObject)snes)->type_name) deft = ((PetscObject)snes)->type_name;
670: PetscOptionsList("-snes_type","Nonlinear solver method","SNESSetType",SNESList,deft,type,256,&flg);
671: if (flg) {
672: SNESSetType(snes,type);
673: } else if (!((PetscObject)snes)->type_name) {
674: SNESSetType(snes,deft);
675: }
676: /* not used here, but called so will go into help messaage */
677: PetscOptionsName("-snes_view","Print detailed information on solver used","SNESView",0);
679: PetscOptionsReal("-snes_stol","Stop if step length less than","SNESSetTolerances",snes->stol,&snes->stol,0);
680: PetscOptionsReal("-snes_atol","Stop if function norm less than","SNESSetTolerances",snes->abstol,&snes->abstol,0);
682: PetscOptionsReal("-snes_rtol","Stop if decrease in function norm less than","SNESSetTolerances",snes->rtol,&snes->rtol,0);
683: PetscOptionsInt("-snes_max_it","Maximum iterations","SNESSetTolerances",snes->max_its,&snes->max_its,NULL);
684: PetscOptionsInt("-snes_max_funcs","Maximum function evaluations","SNESSetTolerances",snes->max_funcs,&snes->max_funcs,NULL);
685: PetscOptionsInt("-snes_max_fail","Maximum nonlinear step failures","SNESSetMaxNonlinearStepFailures",snes->maxFailures,&snes->maxFailures,NULL);
686: PetscOptionsInt("-snes_max_linear_solve_fail","Maximum failures in linear solves allowed","SNESSetMaxLinearSolveFailures",snes->maxLinearSolveFailures,&snes->maxLinearSolveFailures,NULL);
687: PetscOptionsBool("-snes_error_if_not_converged","Generate error if solver does not converge","SNESSetErrorIfNotConverged",snes->errorifnotconverged,&snes->errorifnotconverged,NULL);
689: PetscOptionsInt("-snes_lag_preconditioner","How often to rebuild preconditioner","SNESSetLagPreconditioner",snes->lagpreconditioner,&lag,&flg);
690: if (flg) {
691: SNESSetLagPreconditioner(snes,lag);
692: }
693: PetscOptionsInt("-snes_lag_jacobian","How often to rebuild Jacobian","SNESSetLagJacobian",snes->lagjacobian,&lag,&flg);
694: if (flg) {
695: SNESSetLagJacobian(snes,lag);
696: }
697: PetscOptionsInt("-snes_grid_sequence","Use grid sequencing to generate initial guess","SNESSetGridSequence",snes->gridsequence,&grids,&flg);
698: if (flg) {
699: SNESSetGridSequence(snes,grids);
700: }
702: PetscOptionsEList("-snes_convergence_test","Convergence test","SNESSetConvergenceTest",convtests,2,"default",&indx,&flg);
703: if (flg) {
704: switch (indx) {
705: case 0: SNESSetConvergenceTest(snes,SNESConvergedDefault,NULL,NULL); break;
706: case 1: SNESSetConvergenceTest(snes,SNESSkipConverged,NULL,NULL); break;
707: }
708: }
710: PetscOptionsBool("-snes_converged_reason","Print reason for converged or diverged","SNESSolve",snes->printreason,&snes->printreason,NULL);
712: PetscOptionsEList("-snes_norm_type","SNES Norm type","SNESSetNormType",SNESNormTypes,5,"function",&indx,&flg);
713: if (flg) { SNESSetNormType(snes,(SNESNormType)indx); }
715: kctx = (SNESKSPEW*)snes->kspconvctx;
717: PetscOptionsBool("-snes_ksp_ew","Use Eisentat-Walker linear system convergence test","SNESKSPSetUseEW",snes->ksp_ewconv,&snes->ksp_ewconv,NULL);
719: PetscOptionsInt("-snes_ksp_ew_version","Version 1, 2 or 3","SNESKSPSetParametersEW",kctx->version,&kctx->version,0);
720: PetscOptionsReal("-snes_ksp_ew_rtol0","0 <= rtol0 < 1","SNESKSPSetParametersEW",kctx->rtol_0,&kctx->rtol_0,0);
721: PetscOptionsReal("-snes_ksp_ew_rtolmax","0 <= rtolmax < 1","SNESKSPSetParametersEW",kctx->rtol_max,&kctx->rtol_max,0);
722: PetscOptionsReal("-snes_ksp_ew_gamma","0 <= gamma <= 1","SNESKSPSetParametersEW",kctx->gamma,&kctx->gamma,0);
723: PetscOptionsReal("-snes_ksp_ew_alpha","1 < alpha <= 2","SNESKSPSetParametersEW",kctx->alpha,&kctx->alpha,0);
724: PetscOptionsReal("-snes_ksp_ew_alpha2","alpha2","SNESKSPSetParametersEW",kctx->alpha2,&kctx->alpha2,0);
725: PetscOptionsReal("-snes_ksp_ew_threshold","0 < threshold < 1","SNESKSPSetParametersEW",kctx->threshold,&kctx->threshold,0);
727: flg = PETSC_FALSE;
728: PetscOptionsBool("-snes_check_jacobian","Check each Jacobian with a differenced one","SNESUpdateCheckJacobian",flg,&flg,NULL);
729: if (flg) {
730: SNESSetUpdate(snes,SNESUpdateCheckJacobian);
731: }
733: flg = PETSC_FALSE;
734: PetscOptionsBool("-snes_monitor_cancel","Remove all monitors","SNESMonitorCancel",flg,&flg,NULL);
735: if (flg) {SNESMonitorCancel(snes);}
737: PetscOptionsString("-snes_monitor","Monitor norm of function","SNESMonitorSet","stdout",monfilename,PETSC_MAX_PATH_LEN,&flg);
738: if (flg) {
739: PetscViewerASCIIOpen(PetscObjectComm((PetscObject)snes),monfilename,&monviewer);
740: SNESMonitorSet(snes,SNESMonitorDefault,monviewer,(PetscErrorCode (*)(void**))PetscViewerDestroy);
741: }
743: PetscOptionsString("-snes_monitor_range","Monitor range of elements of function","SNESMonitorSet","stdout",monfilename,PETSC_MAX_PATH_LEN,&flg);
744: if (flg) {
745: SNESMonitorSet(snes,SNESMonitorRange,0,0);
746: }
748: PetscOptionsString("-snes_ratiomonitor","Monitor ratios of norms of function","SNESMonitorSetRatio","stdout",monfilename,PETSC_MAX_PATH_LEN,&flg);
749: if (flg) {
750: PetscViewerASCIIOpen(PetscObjectComm((PetscObject)snes),monfilename,&monviewer);
751: SNESMonitorSetRatio(snes,monviewer);
752: }
754: PetscOptionsString("-snes_monitor_short","Monitor norm of function (fewer digits)","SNESMonitorSet","stdout",monfilename,PETSC_MAX_PATH_LEN,&flg);
755: if (flg) {
756: PetscViewerASCIIOpen(PetscObjectComm((PetscObject)snes),monfilename,&monviewer);
757: SNESMonitorSet(snes,SNESMonitorDefaultShort,monviewer,(PetscErrorCode (*)(void**))PetscViewerDestroy);
758: }
760: PetscOptionsString("-snes_monitor_python","Use Python function","SNESMonitorSet",0,monfilename,PETSC_MAX_PATH_LEN,&flg);
761: if (flg) {PetscPythonMonitorSet((PetscObject)snes,monfilename);}
763: flg = PETSC_FALSE;
764: PetscOptionsBool("-snes_monitor_solution","Plot solution at each iteration","SNESMonitorSolution",flg,&flg,NULL);
765: if (flg) {SNESMonitorSet(snes,SNESMonitorSolution,0,0);}
766: flg = PETSC_FALSE;
767: PetscOptionsBool("-snes_monitor_solution_update","Plot correction at each iteration","SNESMonitorSolutionUpdate",flg,&flg,NULL);
768: if (flg) {SNESMonitorSet(snes,SNESMonitorSolutionUpdate,0,0);}
769: flg = PETSC_FALSE;
770: PetscOptionsBool("-snes_monitor_residual","Plot residual at each iteration","SNESMonitorResidual",flg,&flg,NULL);
771: if (flg) {SNESMonitorSet(snes,SNESMonitorResidual,0,0);}
772: flg = PETSC_FALSE;
773: PetscOptionsBool("-snes_monitor_lg_residualnorm","Plot function norm at each iteration","SNESMonitorLGResidualNorm",flg,&flg,NULL);
774: if (flg) {
775: PetscDrawLG ctx;
777: SNESMonitorLGCreate(0,0,PETSC_DECIDE,PETSC_DECIDE,300,300,&ctx);
778: SNESMonitorSet(snes,SNESMonitorLGResidualNorm,ctx,(PetscErrorCode (*)(void**))SNESMonitorLGDestroy);
779: }
780: flg = PETSC_FALSE;
781: PetscOptionsBool("-snes_monitor_lg_range","Plot function range at each iteration","SNESMonitorLGRange",flg,&flg,NULL);
782: if (flg) {
783: PetscViewer ctx;
785: PetscViewerDrawOpen(PetscObjectComm((PetscObject)snes),0,0,PETSC_DECIDE,PETSC_DECIDE,300,300,&ctx);
786: SNESMonitorSet(snes,SNESMonitorLGRange,ctx,(PetscErrorCode (*)(void**))PetscViewerDestroy);
787: }
789: flg = PETSC_FALSE;
790: PetscOptionsBool("-snes_monitor_jacupdate_spectrum","Print the change in the spectrum of the Jacobian","SNESMonitorJacUpdateSpectrum",flg,&flg,NULL);
791: if (flg) {SNESMonitorSet(snes,SNESMonitorJacUpdateSpectrum,0,0);}
793: flg = PETSC_FALSE;
794: PetscOptionsBool("-snes_fd","Use finite differences (slow) to compute Jacobian","SNESComputeJacobianDefault",flg,&flg,NULL);
795: if (flg) {
796: void *functx;
797: SNESGetFunction(snes,NULL,NULL,&functx);
798: SNESSetJacobian(snes,snes->jacobian,snes->jacobian_pre,SNESComputeJacobianDefault,functx);
799: PetscInfo(snes,"Setting default finite difference Jacobian matrix\n");
800: }
802: flg = PETSC_FALSE;
803: PetscOptionsBool("-snes_fd_function","Use finite differences (slow) to compute function from user objective","SNESObjectiveComputeFunctionDefaultFD",flg,&flg,NULL);
804: if (flg) {
805: SNESSetFunction(snes,NULL,SNESObjectiveComputeFunctionDefaultFD,NULL);
806: }
808: flg = PETSC_FALSE;
809: PetscOptionsBool("-snes_fd_color","Use finite differences with coloring to compute Jacobian","SNESComputeJacobianDefaultColor",flg,&flg,NULL);
810: if (flg) {
811: DM dm;
812: DMSNES sdm;
813: SNESGetDM(snes,&dm);
814: DMGetDMSNES(dm,&sdm);
815: sdm->jacobianctx = NULL;
816: SNESSetJacobian(snes,snes->jacobian,snes->jacobian_pre,SNESComputeJacobianDefaultColor,0);
817: PetscInfo(snes,"Setting default finite difference coloring Jacobian matrix\n");
818: }
820: flg = PETSC_FALSE;
821: PetscOptionsBool("-snes_mf_operator","Use a Matrix-Free Jacobian with user-provided preconditioner matrix","MatCreateSNESMF",PETSC_FALSE,&snes->mf_operator,&flg);
822: if (flg && snes->mf_operator) {
823: snes->mf_operator = PETSC_TRUE;
824: snes->mf = PETSC_TRUE;
825: }
826: flg = PETSC_FALSE;
827: PetscOptionsBool("-snes_mf","Use a Matrix-Free Jacobian with no preconditioner matrix","MatCreateSNESMF",PETSC_FALSE,&snes->mf,&flg);
828: if (!flg && snes->mf_operator) snes->mf = PETSC_TRUE;
829: PetscOptionsInt("-snes_mf_version","Matrix-Free routines version 1 or 2","None",snes->mf_version,&snes->mf_version,0);
831: flg = PETSC_FALSE;
832: SNESGetPCSide(snes,&pcside);
833: PetscOptionsEnum("-snes_npc_side","SNES nonlinear preconditioner side","SNESSetPCSide",PCSides,(PetscEnum)pcside,(PetscEnum*)&pcside,&flg);
834: if (flg) {SNESSetPCSide(snes,pcside);}
836: #if defined(PETSC_HAVE_AMS)
837: {
838: PetscBool set;
839: flg = PETSC_FALSE;
840: PetscOptionsBool("-snes_ams_block","Block for AMS memory snooper at end of SNESSolve","PetscObjectAMSBlock",((PetscObject)snes)->amspublishblock,&flg,&set);
841: if (set) {
842: PetscObjectAMSSetBlock((PetscObject)snes,flg);
843: }
844: }
845: #endif
847: for (i = 0; i < numberofsetfromoptions; i++) {
848: (*othersetfromoptions[i])(snes);
849: }
851: if (snes->ops->setfromoptions) {
852: (*snes->ops->setfromoptions)(snes);
853: }
855: /* process any options handlers added with PetscObjectAddOptionsHandler() */
856: PetscObjectProcessOptionsHandlers((PetscObject)snes);
857: PetscOptionsEnd();
859: if (!snes->ksp) {SNESGetKSP(snes,&snes->ksp);}
860: KSPGetOperators(snes->ksp,NULL,NULL,&matflag);
861: KSPSetOperators(snes->ksp,snes->jacobian,snes->jacobian_pre,matflag);
862: KSPSetFromOptions(snes->ksp);
864: if (!snes->linesearch) {
865: SNESGetLineSearch(snes, &snes->linesearch);
866: }
867: SNESLineSearchSetFromOptions(snes->linesearch);
869: /* if someone has set the SNES PC type, create it. */
870: SNESGetOptionsPrefix(snes, &optionsprefix);
871: PetscOptionsHasName(optionsprefix, "-npc_snes_type", &pcset);
872: if (pcset && (!snes->pc)) {
873: SNESGetPC(snes, &snes->pc);
874: }
875: return(0);
876: }
880: /*@C
881: SNESSetComputeApplicationContext - Sets an optional function to compute a user-defined context for
882: the nonlinear solvers.
884: Logically Collective on SNES
886: Input Parameters:
887: + snes - the SNES context
888: . compute - function to compute the context
889: - destroy - function to destroy the context
891: Level: intermediate
893: Notes:
894: This function is currently not available from Fortran.
896: .keywords: SNES, nonlinear, set, application, context
898: .seealso: SNESGetApplicationContext(), SNESSetComputeApplicationContext(), SNESGetApplicationContext()
899: @*/
900: PetscErrorCode SNESSetComputeApplicationContext(SNES snes,PetscErrorCode (*compute)(SNES,void**),PetscErrorCode (*destroy)(void**))
901: {
904: snes->ops->usercompute = compute;
905: snes->ops->userdestroy = destroy;
906: return(0);
907: }
911: /*@
912: SNESSetApplicationContext - Sets the optional user-defined context for
913: the nonlinear solvers.
915: Logically Collective on SNES
917: Input Parameters:
918: + snes - the SNES context
919: - usrP - optional user context
921: Level: intermediate
923: .keywords: SNES, nonlinear, set, application, context
925: .seealso: SNESGetApplicationContext()
926: @*/
927: PetscErrorCode SNESSetApplicationContext(SNES snes,void *usrP)
928: {
930: KSP ksp;
934: SNESGetKSP(snes,&ksp);
935: KSPSetApplicationContext(ksp,usrP);
936: snes->user = usrP;
937: return(0);
938: }
942: /*@
943: SNESGetApplicationContext - Gets the user-defined context for the
944: nonlinear solvers.
946: Not Collective
948: Input Parameter:
949: . snes - SNES context
951: Output Parameter:
952: . usrP - user context
954: Level: intermediate
956: .keywords: SNES, nonlinear, get, application, context
958: .seealso: SNESSetApplicationContext()
959: @*/
960: PetscErrorCode SNESGetApplicationContext(SNES snes,void *usrP)
961: {
964: *(void**)usrP = snes->user;
965: return(0);
966: }
970: /*@
971: SNESGetIterationNumber - Gets the number of nonlinear iterations completed
972: at this time.
974: Not Collective
976: Input Parameter:
977: . snes - SNES context
979: Output Parameter:
980: . iter - iteration number
982: Notes:
983: For example, during the computation of iteration 2 this would return 1.
985: This is useful for using lagged Jacobians (where one does not recompute the
986: Jacobian at each SNES iteration). For example, the code
987: .vb
988: SNESGetIterationNumber(snes,&it);
989: if (!(it % 2)) {
990: [compute Jacobian here]
991: }
992: .ve
993: can be used in your ComputeJacobian() function to cause the Jacobian to be
994: recomputed every second SNES iteration.
996: Level: intermediate
998: .keywords: SNES, nonlinear, get, iteration, number,
1000: .seealso: SNESGetFunctionNorm(), SNESGetLinearSolveIterations()
1001: @*/
1002: PetscErrorCode SNESGetIterationNumber(SNES snes,PetscInt *iter)
1003: {
1007: *iter = snes->iter;
1008: return(0);
1009: }
1013: /*@
1014: SNESSetIterationNumber - Sets the current iteration number.
1016: Not Collective
1018: Input Parameter:
1019: . snes - SNES context
1020: . iter - iteration number
1022: Level: developer
1024: .keywords: SNES, nonlinear, set, iteration, number,
1026: .seealso: SNESGetFunctionNorm(), SNESGetLinearSolveIterations()
1027: @*/
1028: PetscErrorCode SNESSetIterationNumber(SNES snes,PetscInt iter)
1029: {
1034: PetscObjectAMSTakeAccess((PetscObject)snes);
1035: snes->iter = iter;
1036: PetscObjectAMSGrantAccess((PetscObject)snes);
1037: return(0);
1038: }
1042: /*@
1043: SNESGetFunctionNorm - Gets the norm of the current function that was set
1044: with SNESSSetFunction().
1046: Collective on SNES
1048: Input Parameter:
1049: . snes - SNES context
1051: Output Parameter:
1052: . fnorm - 2-norm of function
1054: Level: intermediate
1056: .keywords: SNES, nonlinear, get, function, norm
1058: .seealso: SNESGetFunction(), SNESGetIterationNumber(), SNESGetLinearSolveIterations()
1059: @*/
1060: PetscErrorCode SNESGetFunctionNorm(SNES snes,PetscReal *fnorm)
1061: {
1065: *fnorm = snes->norm;
1066: return(0);
1067: }
1072: /*@
1073: SNESSetFunctionNorm - Sets the 2-norm of the current function computed using VecNorm().
1075: Collective on SNES
1077: Input Parameter:
1078: . snes - SNES context
1079: . fnorm - 2-norm of function
1081: Level: developer
1083: .keywords: SNES, nonlinear, set, function, norm
1085: .seealso: SNESSetFunction(), SNESSetIterationNumber(), VecNorm().
1086: @*/
1087: PetscErrorCode SNESSetFunctionNorm(SNES snes,PetscReal fnorm)
1088: {
1094: PetscObjectAMSTakeAccess((PetscObject)snes);
1095: snes->norm = fnorm;
1096: PetscObjectAMSGrantAccess((PetscObject)snes);
1097: return(0);
1098: }
1102: /*@
1103: SNESGetNonlinearStepFailures - Gets the number of unsuccessful steps
1104: attempted by the nonlinear solver.
1106: Not Collective
1108: Input Parameter:
1109: . snes - SNES context
1111: Output Parameter:
1112: . nfails - number of unsuccessful steps attempted
1114: Notes:
1115: This counter is reset to zero for each successive call to SNESSolve().
1117: Level: intermediate
1119: .keywords: SNES, nonlinear, get, number, unsuccessful, steps
1121: .seealso: SNESGetMaxLinearSolveFailures(), SNESGetLinearSolveIterations(), SNESSetMaxLinearSolveFailures(), SNESGetLinearSolveFailures(),
1122: SNESSetMaxNonlinearStepFailures(), SNESGetMaxNonlinearStepFailures()
1123: @*/
1124: PetscErrorCode SNESGetNonlinearStepFailures(SNES snes,PetscInt *nfails)
1125: {
1129: *nfails = snes->numFailures;
1130: return(0);
1131: }
1135: /*@
1136: SNESSetMaxNonlinearStepFailures - Sets the maximum number of unsuccessful steps
1137: attempted by the nonlinear solver before it gives up.
1139: Not Collective
1141: Input Parameters:
1142: + snes - SNES context
1143: - maxFails - maximum of unsuccessful steps
1145: Level: intermediate
1147: .keywords: SNES, nonlinear, set, maximum, unsuccessful, steps
1149: .seealso: SNESGetMaxLinearSolveFailures(), SNESGetLinearSolveIterations(), SNESSetMaxLinearSolveFailures(), SNESGetLinearSolveFailures(),
1150: SNESGetMaxNonlinearStepFailures(), SNESGetNonlinearStepFailures()
1151: @*/
1152: PetscErrorCode SNESSetMaxNonlinearStepFailures(SNES snes, PetscInt maxFails)
1153: {
1156: snes->maxFailures = maxFails;
1157: return(0);
1158: }
1162: /*@
1163: SNESGetMaxNonlinearStepFailures - Gets the maximum number of unsuccessful steps
1164: attempted by the nonlinear solver before it gives up.
1166: Not Collective
1168: Input Parameter:
1169: . snes - SNES context
1171: Output Parameter:
1172: . maxFails - maximum of unsuccessful steps
1174: Level: intermediate
1176: .keywords: SNES, nonlinear, get, maximum, unsuccessful, steps
1178: .seealso: SNESGetMaxLinearSolveFailures(), SNESGetLinearSolveIterations(), SNESSetMaxLinearSolveFailures(), SNESGetLinearSolveFailures(),
1179: SNESSetMaxNonlinearStepFailures(), SNESGetNonlinearStepFailures()
1181: @*/
1182: PetscErrorCode SNESGetMaxNonlinearStepFailures(SNES snes, PetscInt *maxFails)
1183: {
1187: *maxFails = snes->maxFailures;
1188: return(0);
1189: }
1193: /*@
1194: SNESGetNumberFunctionEvals - Gets the number of user provided function evaluations
1195: done by SNES.
1197: Not Collective
1199: Input Parameter:
1200: . snes - SNES context
1202: Output Parameter:
1203: . nfuncs - number of evaluations
1205: Level: intermediate
1207: .keywords: SNES, nonlinear, get, maximum, unsuccessful, steps
1209: .seealso: SNESGetMaxLinearSolveFailures(), SNESGetLinearSolveIterations(), SNESSetMaxLinearSolveFailures(), SNESGetLinearSolveFailures()
1210: @*/
1211: PetscErrorCode SNESGetNumberFunctionEvals(SNES snes, PetscInt *nfuncs)
1212: {
1216: *nfuncs = snes->nfuncs;
1217: return(0);
1218: }
1222: /*@
1223: SNESGetLinearSolveFailures - Gets the number of failed (non-converged)
1224: linear solvers.
1226: Not Collective
1228: Input Parameter:
1229: . snes - SNES context
1231: Output Parameter:
1232: . nfails - number of failed solves
1234: Notes:
1235: This counter is reset to zero for each successive call to SNESSolve().
1237: Level: intermediate
1239: .keywords: SNES, nonlinear, get, number, unsuccessful, steps
1241: .seealso: SNESGetMaxLinearSolveFailures(), SNESGetLinearSolveIterations(), SNESSetMaxLinearSolveFailures()
1242: @*/
1243: PetscErrorCode SNESGetLinearSolveFailures(SNES snes,PetscInt *nfails)
1244: {
1248: *nfails = snes->numLinearSolveFailures;
1249: return(0);
1250: }
1254: /*@
1255: SNESSetMaxLinearSolveFailures - the number of failed linear solve attempts
1256: allowed before SNES returns with a diverged reason of SNES_DIVERGED_LINEAR_SOLVE
1258: Logically Collective on SNES
1260: Input Parameters:
1261: + snes - SNES context
1262: - maxFails - maximum allowed linear solve failures
1264: Level: intermediate
1266: Notes: By default this is 0; that is SNES returns on the first failed linear solve
1268: .keywords: SNES, nonlinear, set, maximum, unsuccessful, steps
1270: .seealso: SNESGetLinearSolveFailures(), SNESGetMaxLinearSolveFailures(), SNESGetLinearSolveIterations()
1271: @*/
1272: PetscErrorCode SNESSetMaxLinearSolveFailures(SNES snes, PetscInt maxFails)
1273: {
1277: snes->maxLinearSolveFailures = maxFails;
1278: return(0);
1279: }
1283: /*@
1284: SNESGetMaxLinearSolveFailures - gets the maximum number of linear solve failures that
1285: are allowed before SNES terminates
1287: Not Collective
1289: Input Parameter:
1290: . snes - SNES context
1292: Output Parameter:
1293: . maxFails - maximum of unsuccessful solves allowed
1295: Level: intermediate
1297: Notes: By default this is 1; that is SNES returns on the first failed linear solve
1299: .keywords: SNES, nonlinear, get, maximum, unsuccessful, steps
1301: .seealso: SNESGetLinearSolveFailures(), SNESGetLinearSolveIterations(), SNESSetMaxLinearSolveFailures(),
1302: @*/
1303: PetscErrorCode SNESGetMaxLinearSolveFailures(SNES snes, PetscInt *maxFails)
1304: {
1308: *maxFails = snes->maxLinearSolveFailures;
1309: return(0);
1310: }
1314: /*@
1315: SNESGetLinearSolveIterations - Gets the total number of linear iterations
1316: used by the nonlinear solver.
1318: Not Collective
1320: Input Parameter:
1321: . snes - SNES context
1323: Output Parameter:
1324: . lits - number of linear iterations
1326: Notes:
1327: This counter is reset to zero for each successive call to SNESSolve().
1329: Level: intermediate
1331: .keywords: SNES, nonlinear, get, number, linear, iterations
1333: .seealso: SNESGetIterationNumber(), SNESGetFunctionNorm(), SNESGetLinearSolveFailures(), SNESGetMaxLinearSolveFailures()
1334: @*/
1335: PetscErrorCode SNESGetLinearSolveIterations(SNES snes,PetscInt *lits)
1336: {
1340: *lits = snes->linear_its;
1341: return(0);
1342: }
1347: /*@
1348: SNESSetKSP - Sets a KSP context for the SNES object to use
1350: Not Collective, but the SNES and KSP objects must live on the same MPI_Comm
1352: Input Parameters:
1353: + snes - the SNES context
1354: - ksp - the KSP context
1356: Notes:
1357: The SNES object already has its KSP object, you can obtain with SNESGetKSP()
1358: so this routine is rarely needed.
1360: The KSP object that is already in the SNES object has its reference count
1361: decreased by one.
1363: Level: developer
1365: .keywords: SNES, nonlinear, get, KSP, context
1367: .seealso: KSPGetPC(), SNESCreate(), KSPCreate(), SNESSetKSP()
1368: @*/
1369: PetscErrorCode SNESSetKSP(SNES snes,KSP ksp)
1370: {
1377: PetscObjectReference((PetscObject)ksp);
1378: if (snes->ksp) {PetscObjectDereference((PetscObject)snes->ksp);}
1379: snes->ksp = ksp;
1380: return(0);
1381: }
1383: /* -----------------------------------------------------------*/
1386: /*@
1387: SNESCreate - Creates a nonlinear solver context.
1389: Collective on MPI_Comm
1391: Input Parameters:
1392: . comm - MPI communicator
1394: Output Parameter:
1395: . outsnes - the new SNES context
1397: Options Database Keys:
1398: + -snes_mf - Activates default matrix-free Jacobian-vector products,
1399: and no preconditioning matrix
1400: . -snes_mf_operator - Activates default matrix-free Jacobian-vector
1401: products, and a user-provided preconditioning matrix
1402: as set by SNESSetJacobian()
1403: - -snes_fd - Uses (slow!) finite differences to compute Jacobian
1405: Level: beginner
1407: .keywords: SNES, nonlinear, create, context
1409: .seealso: SNESSolve(), SNESDestroy(), SNES, SNESSetLagPreconditioner()
1411: @*/
1412: PetscErrorCode SNESCreate(MPI_Comm comm,SNES *outsnes)
1413: {
1415: SNES snes;
1416: SNESKSPEW *kctx;
1420: *outsnes = NULL;
1421: #if !defined(PETSC_USE_DYNAMIC_LIBRARIES)
1422: SNESInitializePackage();
1423: #endif
1425: PetscHeaderCreate(snes,_p_SNES,struct _SNESOps,SNES_CLASSID,"SNES","Nonlinear solver","SNES",comm,SNESDestroy,SNESView);
1427: snes->ops->converged = SNESConvergedDefault;
1428: snes->usesksp = PETSC_TRUE;
1429: snes->tolerancesset = PETSC_FALSE;
1430: snes->max_its = 50;
1431: snes->max_funcs = 10000;
1432: snes->norm = 0.0;
1433: snes->normtype = SNES_NORM_FUNCTION;
1434: snes->rtol = 1.e-8;
1435: snes->ttol = 0.0;
1436: snes->abstol = 1.e-50;
1437: snes->stol = 1.e-8;
1438: snes->deltatol = 1.e-12;
1439: snes->nfuncs = 0;
1440: snes->numFailures = 0;
1441: snes->maxFailures = 1;
1442: snes->linear_its = 0;
1443: snes->lagjacobian = 1;
1444: snes->lagpreconditioner = 1;
1445: snes->numbermonitors = 0;
1446: snes->data = 0;
1447: snes->setupcalled = PETSC_FALSE;
1448: snes->ksp_ewconv = PETSC_FALSE;
1449: snes->nwork = 0;
1450: snes->work = 0;
1451: snes->nvwork = 0;
1452: snes->vwork = 0;
1453: snes->conv_hist_len = 0;
1454: snes->conv_hist_max = 0;
1455: snes->conv_hist = NULL;
1456: snes->conv_hist_its = NULL;
1457: snes->conv_hist_reset = PETSC_TRUE;
1458: snes->vec_func_init_set = PETSC_FALSE;
1459: snes->norm_init = 0.;
1460: snes->norm_init_set = PETSC_FALSE;
1461: snes->reason = SNES_CONVERGED_ITERATING;
1462: snes->pcside = PC_RIGHT;
1464: snes->mf = PETSC_FALSE;
1465: snes->mf_operator = PETSC_FALSE;
1466: snes->mf_version = 1;
1468: snes->numLinearSolveFailures = 0;
1469: snes->maxLinearSolveFailures = 1;
1471: /* Create context to compute Eisenstat-Walker relative tolerance for KSP */
1472: PetscNewLog(snes,SNESKSPEW,&kctx);
1474: snes->kspconvctx = (void*)kctx;
1475: kctx->version = 2;
1476: kctx->rtol_0 = .3; /* Eisenstat and Walker suggest rtol_0=.5, but
1477: this was too large for some test cases */
1478: kctx->rtol_last = 0.0;
1479: kctx->rtol_max = .9;
1480: kctx->gamma = 1.0;
1481: kctx->alpha = .5*(1.0 + PetscSqrtReal(5.0));
1482: kctx->alpha2 = kctx->alpha;
1483: kctx->threshold = .1;
1484: kctx->lresid_last = 0.0;
1485: kctx->norm_last = 0.0;
1487: *outsnes = snes;
1488: return(0);
1489: }
1491: /*MC
1492: SNESFunction - functional form used to convey the nonlinear function to be solved by SNES
1494: Synopsis:
1495: #include "petscsnes.h"
1496: SNESFunction(SNES snes,Vec x,Vec f,void *ctx);
1498: Input Parameters:
1499: + snes - the SNES context
1500: . x - state at which to evaluate residual
1501: - ctx - optional user-defined function context, passed in with SNESSetFunction()
1503: Output Parameter:
1504: . f - vector to put residual (function value)
1506: Level: intermediate
1508: .seealso: SNESSetFunction(), SNESGetFunction()
1509: M*/
1513: /*@C
1514: SNESSetFunction - Sets the function evaluation routine and function
1515: vector for use by the SNES routines in solving systems of nonlinear
1516: equations.
1518: Logically Collective on SNES
1520: Input Parameters:
1521: + snes - the SNES context
1522: . r - vector to store function value
1523: . SNESFunction - function evaluation routine
1524: - ctx - [optional] user-defined context for private data for the
1525: function evaluation routine (may be NULL)
1527: Notes:
1528: The Newton-like methods typically solve linear systems of the form
1529: $ f'(x) x = -f(x),
1530: where f'(x) denotes the Jacobian matrix and f(x) is the function.
1532: Level: beginner
1534: .keywords: SNES, nonlinear, set, function
1536: .seealso: SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESSetPicard(), SNESFunction
1537: @*/
1538: PetscErrorCode SNESSetFunction(SNES snes,Vec r,PetscErrorCode (*SNESFunction)(SNES,Vec,Vec,void*),void *ctx)
1539: {
1541: DM dm;
1545: if (r) {
1548: PetscObjectReference((PetscObject)r);
1549: VecDestroy(&snes->vec_func);
1551: snes->vec_func = r;
1552: }
1553: SNESGetDM(snes,&dm);
1554: DMSNESSetFunction(dm,SNESFunction,ctx);
1555: return(0);
1556: }
1561: /*@C
1562: SNESSetInitialFunction - Sets the function vector to be used as the
1563: function norm at the initialization of the method. In some
1564: instances, the user has precomputed the function before calling
1565: SNESSolve. This function allows one to avoid a redundant call
1566: to SNESComputeFunction in that case.
1568: Logically Collective on SNES
1570: Input Parameters:
1571: + snes - the SNES context
1572: - f - vector to store function value
1574: Notes:
1575: This should not be modified during the solution procedure.
1577: This is used extensively in the SNESFAS hierarchy and in nonlinear preconditioning.
1579: Level: developer
1581: .keywords: SNES, nonlinear, set, function
1583: .seealso: SNESSetFunction(), SNESComputeFunction(), SNESSetInitialFunctionNorm()
1584: @*/
1585: PetscErrorCode SNESSetInitialFunction(SNES snes, Vec f)
1586: {
1588: Vec vec_func;
1594: SNESGetFunction(snes,&vec_func,NULL,NULL);
1595: VecCopy(f, vec_func);
1597: snes->vec_func_init_set = PETSC_TRUE;
1598: return(0);
1599: }
1604: /*@C
1605: SNESSetInitialFunctionNorm - Sets the function norm to be used as the function norm
1606: at the initialization of the method. In some instances, the user has precomputed
1607: the function and its norm before calling SNESSolve. This function allows one to
1608: avoid a redundant call to SNESComputeFunction() and VecNorm() in that case.
1610: Logically Collective on SNES
1612: Input Parameters:
1613: + snes - the SNES context
1614: - fnorm - the norm of F as set by SNESSetInitialFunction()
1616: This is used extensively in the SNESFAS hierarchy and in nonlinear preconditioning.
1618: Level: developer
1620: .keywords: SNES, nonlinear, set, function, norm
1622: .seealso: SNESSetFunction(), SNESComputeFunction(), SNESSetInitialFunction()
1623: @*/
1624: PetscErrorCode SNESSetInitialFunctionNorm(SNES snes, PetscReal fnorm)
1625: {
1628: snes->norm_init = fnorm;
1629: snes->norm_init_set = PETSC_TRUE;
1630: return(0);
1631: }
1635: /*@
1636: SNESSetNormType - Sets the SNESNormType used in covergence and monitoring
1637: of the SNES method.
1639: Logically Collective on SNES
1641: Input Parameters:
1642: + snes - the SNES context
1643: - normtype - the type of the norm used
1645: Notes:
1646: Only certain SNES methods support certain SNESNormTypes. Most require evaluation
1647: of the nonlinear function and the taking of its norm at every iteration to
1648: even ensure convergence at all. However, methods such as custom Gauss-Seidel methods
1649: (SNESGS) and the like do not require the norm of the function to be computed, and therfore
1650: may either be monitored for convergence or not. As these are often used as nonlinear
1651: preconditioners, monitoring the norm of their error is not a useful enterprise within
1652: their solution.
1654: Level: developer
1656: .keywords: SNES, nonlinear, set, function, norm, type
1658: .seealso: SNESGetNormType(), SNESComputeFunction(), VecNorm(), SNESSetFunction(), SNESSetInitialFunction(), SNESNormType
1659: @*/
1660: PetscErrorCode SNESSetNormType(SNES snes, SNESNormType normtype)
1661: {
1664: snes->normtype = normtype;
1665: return(0);
1666: }
1671: /*@
1672: SNESGetNormType - Gets the SNESNormType used in covergence and monitoring
1673: of the SNES method.
1675: Logically Collective on SNES
1677: Input Parameters:
1678: + snes - the SNES context
1679: - normtype - the type of the norm used
1681: Level: advanced
1683: .keywords: SNES, nonlinear, set, function, norm, type
1685: .seealso: SNESSetNormType(), SNESComputeFunction(), VecNorm(), SNESSetFunction(), SNESSetInitialFunction(), SNESNormType
1686: @*/
1687: PetscErrorCode SNESGetNormType(SNES snes, SNESNormType *normtype)
1688: {
1691: *normtype = snes->normtype;
1692: return(0);
1693: }
1695: /*MC
1696: SNESGSFunction - function used to convey a Gauss-Seidel sweep on the nonlinear function
1698: Synopsis:
1699: #include "petscsnes.h"
1700: $ SNESGSFunction(SNES snes,Vec x,Vec b,void *ctx);
1702: + X - solution vector
1703: . B - RHS vector
1704: - ctx - optional user-defined Gauss-Seidel context
1706: Level: intermediate
1708: .seealso: SNESSetGS(), SNESGetGS()
1709: M*/
1713: /*@C
1714: SNESSetGS - Sets the user nonlinear Gauss-Seidel routine for
1715: use with composed nonlinear solvers.
1717: Input Parameters:
1718: + snes - the SNES context
1719: . SNESGSFunction - function evaluation routine
1720: - ctx - [optional] user-defined context for private data for the
1721: smoother evaluation routine (may be NULL)
1723: Notes:
1724: The GS routines are used by the composed nonlinear solver to generate
1725: a problem appropriate update to the solution, particularly FAS.
1727: Level: intermediate
1729: .keywords: SNES, nonlinear, set, Gauss-Seidel
1731: .seealso: SNESGetFunction(), SNESComputeGS()
1732: @*/
1733: PetscErrorCode SNESSetGS(SNES snes,PetscErrorCode (*SNESGSFunction)(SNES,Vec,Vec,void*),void *ctx)
1734: {
1736: DM dm;
1740: SNESGetDM(snes,&dm);
1741: DMSNESSetGS(dm,SNESGSFunction,ctx);
1742: return(0);
1743: }
1747: PETSC_EXTERN PetscErrorCode SNESPicardComputeFunction(SNES snes,Vec x,Vec f,void *ctx)
1748: {
1750: DM dm;
1751: DMSNES sdm;
1754: SNESGetDM(snes,&dm);
1755: DMGetDMSNES(dm,&sdm);
1756: /* A(x)*x - b(x) */
1757: if (sdm->ops->computepfunction) {
1758: (*sdm->ops->computepfunction)(snes,x,f,sdm->pctx);
1759: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE, "Must call SNESSetPicard() to provide Picard function.");
1761: if (sdm->ops->computepjacobian) {
1762: (*sdm->ops->computepjacobian)(snes,x,&snes->jacobian,&snes->jacobian_pre,&snes->matstruct,sdm->pctx);
1763: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE, "Must call SNESSetPicard() to provide Picard matrix.");
1764: VecScale(f,-1.0);
1765: MatMultAdd(snes->jacobian,x,f,f);
1766: return(0);
1767: }
1771: PETSC_EXTERN PetscErrorCode SNESPicardComputeJacobian(SNES snes,Vec x1,Mat *J,Mat *B,MatStructure *flag,void *ctx)
1772: {
1774: /* the jacobian matrix should be pre-filled in SNESPicardComputeFunction */
1775: *flag = snes->matstruct;
1776: return(0);
1777: }
1781: /*@C
1782: SNESSetPicard - Use SNES to solve the semilinear-system A(x) x = b(x) via a Picard type iteration (Picard linearization)
1784: Logically Collective on SNES
1786: Input Parameters:
1787: + snes - the SNES context
1788: . r - vector to store function value
1789: . SNESFunction - function evaluation routine
1790: . Amat - matrix with which A(x) x - b(x) is to be computed
1791: . Pmat - matrix from which preconditioner is computed (usually the same as Amat)
1792: . SNESJacobianFunction - function to compute matrix value
1793: - ctx - [optional] user-defined context for private data for the
1794: function evaluation routine (may be NULL)
1796: Notes:
1797: 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
1798: 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.
1800: One can call SNESSetPicard() or SNESSetFunction() (and possibly SNESSetJacobian()) but cannot call both
1802: $ 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}
1803: $ Note that when an exact solver is used this corresponds to the "classic" Picard A(x^{n}) x^{n+1} = b(x^{n}) iteration.
1805: Run with -snes_mf_operator to solve the system with Newton's method using A(x^{n}) to construct the preconditioner.
1807: We implement the defect correction form of the Picard iteration because it converges much more generally when inexact linear solvers are used then
1808: the direct Picard iteration A(x^n) x^{n+1} = b(x^n)
1810: 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
1811: 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
1812: different please contact us at petsc-dev@mcs.anl.gov and we'll have an entirely new argument :-).
1814: Level: intermediate
1816: .keywords: SNES, nonlinear, set, function
1818: .seealso: SNESGetFunction(), SNESSetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESGetPicard(), SNESLineSearchPreCheckPicard()
1819: @*/
1820: PetscErrorCode SNESSetPicard(SNES snes,Vec r,PetscErrorCode (*SNESFunction)(SNES,Vec,Vec,void*),Mat Amat, Mat Pmat, PetscErrorCode (*SNESJacobianFunction)(SNES,Vec,Mat*,Mat*,MatStructure*,void*),void *ctx)
1821: {
1823: DM dm;
1827: SNESGetDM(snes, &dm);
1828: DMSNESSetPicard(dm,SNESFunction,SNESJacobianFunction,ctx);
1829: SNESSetFunction(snes,r,SNESPicardComputeFunction,ctx);
1830: SNESSetJacobian(snes,Amat,Pmat,SNESPicardComputeJacobian,ctx);
1831: return(0);
1832: }
1836: /*@C
1837: SNESGetPicard - Returns the context for the Picard iteration
1839: Not Collective, but Vec is parallel if SNES is parallel. Collective if Vec is requested, but has not been created yet.
1841: Input Parameter:
1842: . snes - the SNES context
1844: Output Parameter:
1845: + r - the function (or NULL)
1846: . SNESFunction - the function (or NULL)
1847: . Amat - the matrix used to defined the operation A(x) x - b(x) (or NULL)
1848: . Pmat - the matrix from which the preconditioner will be constructed (or NULL)
1849: . SNESJacobianFunction - the function for matrix evaluation (or NULL)
1850: - ctx - the function context (or NULL)
1852: Level: advanced
1854: .keywords: SNES, nonlinear, get, function
1856: .seealso: SNESSetPicard(), SNESGetFunction(), SNESGetJacobian(), SNESGetDM(), SNESFunction, SNESJacobianFunction
1857: @*/
1858: PetscErrorCode SNESGetPicard(SNES snes,Vec *r,PetscErrorCode (**SNESFunction)(SNES,Vec,Vec,void*),Mat *Amat, Mat *Pmat, PetscErrorCode (**SNESJacobianFunction)(SNES,Vec,Mat*,Mat*,MatStructure*,void*),void **ctx)
1859: {
1861: DM dm;
1865: SNESGetFunction(snes,r,NULL,NULL);
1866: SNESGetJacobian(snes,Amat,Pmat,NULL,NULL);
1867: SNESGetDM(snes,&dm);
1868: DMSNESGetPicard(dm,SNESFunction,SNESJacobianFunction,ctx);
1869: return(0);
1870: }
1874: /*@C
1875: SNESSetComputeInitialGuess - Sets a routine used to compute an initial guess for the problem
1877: Logically Collective on SNES
1879: Input Parameters:
1880: + snes - the SNES context
1881: . func - function evaluation routine
1882: - ctx - [optional] user-defined context for private data for the
1883: function evaluation routine (may be NULL)
1885: Calling sequence of func:
1886: $ func (SNES snes,Vec x,void *ctx);
1888: . f - function vector
1889: - ctx - optional user-defined function context
1891: Level: intermediate
1893: .keywords: SNES, nonlinear, set, function
1895: .seealso: SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian()
1896: @*/
1897: PetscErrorCode SNESSetComputeInitialGuess(SNES snes,PetscErrorCode (*func)(SNES,Vec,void*),void *ctx)
1898: {
1901: if (func) snes->ops->computeinitialguess = func;
1902: if (ctx) snes->initialguessP = ctx;
1903: return(0);
1904: }
1906: /* --------------------------------------------------------------- */
1909: /*@C
1910: SNESGetRhs - Gets the vector for solving F(x) = rhs. If rhs is not set
1911: it assumes a zero right hand side.
1913: Logically Collective on SNES
1915: Input Parameter:
1916: . snes - the SNES context
1918: Output Parameter:
1919: . rhs - the right hand side vector or NULL if the right hand side vector is null
1921: Level: intermediate
1923: .keywords: SNES, nonlinear, get, function, right hand side
1925: .seealso: SNESGetSolution(), SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESSetFunction()
1926: @*/
1927: PetscErrorCode SNESGetRhs(SNES snes,Vec *rhs)
1928: {
1932: *rhs = snes->vec_rhs;
1933: return(0);
1934: }
1938: /*@
1939: SNESComputeFunction - Calls the function that has been set with SNESSetFunction().
1941: Collective on SNES
1943: Input Parameters:
1944: + snes - the SNES context
1945: - x - input vector
1947: Output Parameter:
1948: . y - function vector, as set by SNESSetFunction()
1950: Notes:
1951: SNESComputeFunction() is typically used within nonlinear solvers
1952: implementations, so most users would not generally call this routine
1953: themselves.
1955: Level: developer
1957: .keywords: SNES, nonlinear, compute, function
1959: .seealso: SNESSetFunction(), SNESGetFunction()
1960: @*/
1961: PetscErrorCode SNESComputeFunction(SNES snes,Vec x,Vec y)
1962: {
1964: DM dm;
1965: DMSNES sdm;
1973: VecValidValues(x,2,PETSC_TRUE);
1975: SNESGetDM(snes,&dm);
1976: DMGetDMSNES(dm,&sdm);
1977: if (snes->pc && snes->pcside == PC_LEFT) {
1978: VecCopy(x,y);
1979: PetscLogEventBegin(SNES_NPCSolve,snes->pc,x,y,0);
1980: SNESSolve(snes->pc,snes->vec_rhs,y);
1981: PetscLogEventEnd(SNES_NPCSolve,snes->pc,x,y,0);
1982: VecAYPX(y,-1.0,x);
1983: VecValidValues(y,3,PETSC_FALSE);
1984: return(0);
1985: } else if (sdm->ops->computefunction) {
1986: PetscLogEventBegin(SNES_FunctionEval,snes,x,y,0);
1987: PetscStackPush("SNES user function");
1988: (*sdm->ops->computefunction)(snes,x,y,sdm->functionctx);
1989: PetscStackPop;
1990: PetscLogEventEnd(SNES_FunctionEval,snes,x,y,0);
1991: } else if (snes->vec_rhs) {
1992: MatMult(snes->jacobian, x, y);
1993: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE, "Must call SNESSetFunction() or SNESSetDM() before SNESComputeFunction(), likely called from SNESSolve().");
1994: if (snes->vec_rhs) {
1995: VecAXPY(y,-1.0,snes->vec_rhs);
1996: }
1997: snes->nfuncs++;
1998: VecValidValues(y,3,PETSC_FALSE);
1999: return(0);
2000: }
2004: /*@
2005: SNESComputeGS - Calls the Gauss-Seidel function that has been set with SNESSetGS().
2007: Collective on SNES
2009: Input Parameters:
2010: + snes - the SNES context
2011: . x - input vector
2012: - b - rhs vector
2014: Output Parameter:
2015: . x - new solution vector
2017: Notes:
2018: SNESComputeGS() is typically used within composed nonlinear solver
2019: implementations, so most users would not generally call this routine
2020: themselves.
2022: Level: developer
2024: .keywords: SNES, nonlinear, compute, function
2026: .seealso: SNESSetGS(), SNESComputeFunction()
2027: @*/
2028: PetscErrorCode SNESComputeGS(SNES snes,Vec b,Vec x)
2029: {
2031: DM dm;
2032: DMSNES sdm;
2040: if (b) {VecValidValues(b,2,PETSC_TRUE);}
2041: PetscLogEventBegin(SNES_GSEval,snes,x,b,0);
2042: SNESGetDM(snes,&dm);
2043: DMGetDMSNES(dm,&sdm);
2044: if (sdm->ops->computegs) {
2045: PetscStackPush("SNES user GS");
2046: (*sdm->ops->computegs)(snes,x,b,sdm->gsctx);
2047: PetscStackPop;
2048: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE, "Must call SNESSetGS() before SNESComputeGS(), likely called from SNESSolve().");
2049: PetscLogEventEnd(SNES_GSEval,snes,x,b,0);
2050: VecValidValues(x,3,PETSC_FALSE);
2051: return(0);
2052: }
2056: /*@
2057: SNESComputeJacobian - Computes the Jacobian matrix that has been set with SNESSetJacobian().
2059: Collective on SNES and Mat
2061: Input Parameters:
2062: + snes - the SNES context
2063: - x - input vector
2065: Output Parameters:
2066: + A - Jacobian matrix
2067: . B - optional preconditioning matrix
2068: - flag - flag indicating matrix structure (one of, SAME_NONZERO_PATTERN,DIFFERENT_NONZERO_PATTERN,SAME_PRECONDITIONER)
2070: Options Database Keys:
2071: + -snes_lag_preconditioner <lag>
2072: . -snes_lag_jacobian <lag>
2073: . -snes_compare_explicit - Compare the computed Jacobian to the finite difference Jacobian and output the differences
2074: . -snes_compare_explicit_draw - Compare the computed Jacobian to the finite difference Jacobian and draw the result
2075: . -snes_compare_explicit_contour - Compare the computed Jacobian to the finite difference Jacobian and draw a contour plot with the result
2076: . -snes_compare_operator - Make the comparison options above use the operator instead of the preconditioning matrix
2077: . -snes_compare_coloring - Compute the finite differece Jacobian using coloring and display norms of difference
2078: . -snes_compare_coloring_display - Compute the finite differece Jacobian using coloring and display verbose differences
2079: . -snes_compare_coloring_threshold - Display only those matrix entries that differ by more than a given threshold
2080: . -snes_compare_coloring_threshold_atol - Absolute tolerance for difference in matrix entries to be displayed by -snes_compare_coloring_threshold
2081: . -snes_compare_coloring_threshold_rtol - Relative tolerance for difference in matrix entries to be displayed by -snes_compare_coloring_threshold
2082: . -snes_compare_coloring_draw - Compute the finite differece Jacobian using coloring and draw differences
2083: - -snes_compare_coloring_draw_contour - Compute the finite differece Jacobian using coloring and show contours of matrices and differences
2086: Notes:
2087: Most users should not need to explicitly call this routine, as it
2088: is used internally within the nonlinear solvers.
2090: See KSPSetOperators() for important information about setting the
2091: flag parameter.
2093: Level: developer
2095: .keywords: SNES, compute, Jacobian, matrix
2097: .seealso: SNESSetJacobian(), KSPSetOperators(), MatStructure, SNESSetLagPreconditioner(), SNESSetLagJacobian()
2098: @*/
2099: PetscErrorCode SNESComputeJacobian(SNES snes,Vec X,Mat *A,Mat *B,MatStructure *flg)
2100: {
2102: PetscBool flag;
2103: DM dm;
2104: DMSNES sdm;
2111: VecValidValues(X,2,PETSC_TRUE);
2112: SNESGetDM(snes,&dm);
2113: DMGetDMSNES(dm,&sdm);
2115: if (!sdm->ops->computejacobian) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_USER,"Must call SNESSetJacobian(), DMSNESSetJacobian(), DMDASNESSetJacobianLocal(), etc");
2117: /* make sure that MatAssemblyBegin/End() is called on A matrix if it is matrix free */
2119: if (snes->lagjacobian == -2) {
2120: snes->lagjacobian = -1;
2122: PetscInfo(snes,"Recomputing Jacobian/preconditioner because lag is -2 (means compute Jacobian, but then never again) \n");
2123: } else if (snes->lagjacobian == -1) {
2124: *flg = SAME_PRECONDITIONER;
2125: PetscInfo(snes,"Reusing Jacobian/preconditioner because lag is -1\n");
2126: PetscObjectTypeCompare((PetscObject)*A,MATMFFD,&flag);
2127: if (flag) {
2128: MatAssemblyBegin(*A,MAT_FINAL_ASSEMBLY);
2129: MatAssemblyEnd(*A,MAT_FINAL_ASSEMBLY);
2130: }
2131: return(0);
2132: } else if (snes->lagjacobian > 1 && snes->iter % snes->lagjacobian) {
2133: *flg = SAME_PRECONDITIONER;
2134: PetscInfo2(snes,"Reusing Jacobian/preconditioner because lag is %D and SNES iteration is %D\n",snes->lagjacobian,snes->iter);
2135: PetscObjectTypeCompare((PetscObject)*A,MATMFFD,&flag);
2136: if (flag) {
2137: MatAssemblyBegin(*A,MAT_FINAL_ASSEMBLY);
2138: MatAssemblyEnd(*A,MAT_FINAL_ASSEMBLY);
2139: }
2140: return(0);
2141: }
2142: if (snes->pc && snes->pcside == PC_LEFT) {
2143: MatAssemblyBegin(*A,MAT_FINAL_ASSEMBLY);
2144: MatAssemblyEnd(*A,MAT_FINAL_ASSEMBLY);
2145: return(0);
2146: }
2148: *flg = DIFFERENT_NONZERO_PATTERN;
2149: PetscLogEventBegin(SNES_JacobianEval,snes,X,*A,*B);
2151: PetscStackPush("SNES user Jacobian function");
2152: (*sdm->ops->computejacobian)(snes,X,A,B,flg,sdm->jacobianctx);
2153: PetscStackPop;
2155: PetscLogEventEnd(SNES_JacobianEval,snes,X,*A,*B);
2157: if (snes->lagpreconditioner == -2) {
2158: PetscInfo(snes,"Rebuilding preconditioner exactly once since lag is -2\n");
2160: snes->lagpreconditioner = -1;
2161: } else if (snes->lagpreconditioner == -1) {
2162: *flg = SAME_PRECONDITIONER;
2163: PetscInfo(snes,"Reusing preconditioner because lag is -1\n");
2164: } else if (snes->lagpreconditioner > 1 && snes->iter % snes->lagpreconditioner) {
2165: *flg = SAME_PRECONDITIONER;
2166: PetscInfo2(snes,"Reusing preconditioner because lag is %D and SNES iteration is %D\n",snes->lagpreconditioner,snes->iter);
2167: }
2169: /* make sure user returned a correct Jacobian and preconditioner */
2172: {
2173: PetscBool flag = PETSC_FALSE,flag_draw = PETSC_FALSE,flag_contour = PETSC_FALSE,flag_operator = PETSC_FALSE;
2174: PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_explicit",&flag,NULL);
2175: PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_explicit_draw",&flag_draw,NULL);
2176: PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_explicit_draw_contour",&flag_contour,NULL);
2177: PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_operator",&flag_operator,NULL);
2178: if (flag || flag_draw || flag_contour) {
2179: Mat Bexp_mine = NULL,Bexp,FDexp;
2180: MatStructure mstruct;
2181: PetscViewer vdraw,vstdout;
2182: PetscBool flg;
2183: if (flag_operator) {
2184: MatComputeExplicitOperator(*A,&Bexp_mine);
2185: Bexp = Bexp_mine;
2186: } else {
2187: /* See if the preconditioning matrix can be viewed and added directly */
2188: PetscObjectTypeCompareAny((PetscObject)*B,&flg,MATSEQAIJ,MATMPIAIJ,MATSEQDENSE,MATMPIDENSE,MATSEQBAIJ,MATMPIBAIJ,MATSEQSBAIJ,MATMPIBAIJ,"");
2189: if (flg) Bexp = *B;
2190: else {
2191: /* If the "preconditioning" matrix is itself MATSHELL or some other type without direct support */
2192: MatComputeExplicitOperator(*B,&Bexp_mine);
2193: Bexp = Bexp_mine;
2194: }
2195: }
2196: MatConvert(Bexp,MATSAME,MAT_INITIAL_MATRIX,&FDexp);
2197: SNESComputeJacobianDefault(snes,X,&FDexp,&FDexp,&mstruct,NULL);
2198: PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)snes),&vstdout);
2199: if (flag_draw || flag_contour) {
2200: PetscViewerDrawOpen(PetscObjectComm((PetscObject)snes),0,"Explicit Jacobians",PETSC_DECIDE,PETSC_DECIDE,300,300,&vdraw);
2201: if (flag_contour) {PetscViewerPushFormat(vdraw,PETSC_VIEWER_DRAW_CONTOUR);}
2202: } else vdraw = NULL;
2203: PetscViewerASCIIPrintf(vstdout,"Explicit %s\n",flag_operator ? "Jacobian" : "preconditioning Jacobian");
2204: if (flag) {MatView(Bexp,vstdout);}
2205: if (vdraw) {MatView(Bexp,vdraw);}
2206: PetscViewerASCIIPrintf(vstdout,"Finite difference Jacobian\n");
2207: if (flag) {MatView(FDexp,vstdout);}
2208: if (vdraw) {MatView(FDexp,vdraw);}
2209: MatAYPX(FDexp,-1.0,Bexp,SAME_NONZERO_PATTERN);
2210: PetscViewerASCIIPrintf(vstdout,"User-provided matrix minus finite difference Jacobian\n");
2211: if (flag) {MatView(FDexp,vstdout);}
2212: if (vdraw) { /* Always use contour for the difference */
2213: PetscViewerPushFormat(vdraw,PETSC_VIEWER_DRAW_CONTOUR);
2214: MatView(FDexp,vdraw);
2215: PetscViewerPopFormat(vdraw);
2216: }
2217: if (flag_contour) {PetscViewerPopFormat(vdraw);}
2218: PetscViewerDestroy(&vdraw);
2219: MatDestroy(&Bexp_mine);
2220: MatDestroy(&FDexp);
2221: }
2222: }
2223: {
2224: PetscBool flag = PETSC_FALSE,flag_display = PETSC_FALSE,flag_draw = PETSC_FALSE,flag_contour = PETSC_FALSE,flag_threshold = PETSC_FALSE;
2225: PetscReal threshold_atol = PETSC_SQRT_MACHINE_EPSILON,threshold_rtol = 10*PETSC_SQRT_MACHINE_EPSILON;
2226: PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_coloring",&flag,NULL);
2227: PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_coloring_display",&flag_display,NULL);
2228: PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_coloring_draw",&flag_draw,NULL);
2229: PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_coloring_draw_contour",&flag_contour,NULL);
2230: PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_coloring_threshold",&flag_threshold,NULL);
2231: PetscOptionsGetReal(((PetscObject)snes)->prefix,"-snes_compare_coloring_threshold_rtol",&threshold_rtol,NULL);
2232: PetscOptionsGetReal(((PetscObject)snes)->prefix,"-snes_compare_coloring_threshold_atol",&threshold_atol,NULL);
2233: if (flag || flag_display || flag_draw || flag_contour || flag_threshold) {
2234: Mat Bfd;
2235: MatStructure mstruct;
2236: PetscViewer vdraw,vstdout;
2237: ISColoring iscoloring;
2238: MatFDColoring matfdcoloring;
2239: PetscErrorCode (*func)(SNES,Vec,Vec,void*);
2240: void *funcctx;
2241: PetscReal norm1,norm2,normmax;
2243: MatDuplicate(*B,MAT_DO_NOT_COPY_VALUES,&Bfd);
2244: MatGetColoring(Bfd,MATCOLORINGSL,&iscoloring);
2245: MatFDColoringCreate(Bfd,iscoloring,&matfdcoloring);
2246: ISColoringDestroy(&iscoloring);
2248: /* This method of getting the function is currently unreliable since it doesn't work for DM local functions. */
2249: SNESGetFunction(snes,NULL,&func,&funcctx);
2250: MatFDColoringSetFunction(matfdcoloring,(PetscErrorCode (*)(void))func,funcctx);
2251: PetscObjectSetOptionsPrefix((PetscObject)matfdcoloring,((PetscObject)snes)->prefix);
2252: PetscObjectAppendOptionsPrefix((PetscObject)matfdcoloring,"coloring_");
2253: MatFDColoringSetFromOptions(matfdcoloring);
2254: MatFDColoringApply(Bfd,matfdcoloring,X,&mstruct,snes);
2255: MatFDColoringDestroy(&matfdcoloring);
2257: PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)snes),&vstdout);
2258: if (flag_draw || flag_contour) {
2259: PetscViewerDrawOpen(PetscObjectComm((PetscObject)snes),0,"Colored Jacobians",PETSC_DECIDE,PETSC_DECIDE,300,300,&vdraw);
2260: if (flag_contour) {PetscViewerPushFormat(vdraw,PETSC_VIEWER_DRAW_CONTOUR);}
2261: } else vdraw = NULL;
2262: PetscViewerASCIIPrintf(vstdout,"Explicit preconditioning Jacobian\n");
2263: if (flag_display) {MatView(*B,vstdout);}
2264: if (vdraw) {MatView(*B,vdraw);}
2265: PetscViewerASCIIPrintf(vstdout,"Colored Finite difference Jacobian\n");
2266: if (flag_display) {MatView(Bfd,vstdout);}
2267: if (vdraw) {MatView(Bfd,vdraw);}
2268: MatAYPX(Bfd,-1.0,*B,SAME_NONZERO_PATTERN);
2269: MatNorm(Bfd,NORM_1,&norm1);
2270: MatNorm(Bfd,NORM_FROBENIUS,&norm2);
2271: MatNorm(Bfd,NORM_MAX,&normmax);
2272: PetscViewerASCIIPrintf(vstdout,"User-provided matrix minus finite difference Jacobian, norm1=%G normFrob=%G normmax=%G\n",norm1,norm2,normmax);
2273: if (flag_display) {MatView(Bfd,vstdout);}
2274: if (vdraw) { /* Always use contour for the difference */
2275: PetscViewerPushFormat(vdraw,PETSC_VIEWER_DRAW_CONTOUR);
2276: MatView(Bfd,vdraw);
2277: PetscViewerPopFormat(vdraw);
2278: }
2279: if (flag_contour) {PetscViewerPopFormat(vdraw);}
2281: if (flag_threshold) {
2282: PetscInt bs,rstart,rend,i;
2283: MatGetBlockSize(*B,&bs);
2284: MatGetOwnershipRange(*B,&rstart,&rend);
2285: for (i=rstart; i<rend; i++) {
2286: const PetscScalar *ba,*ca;
2287: const PetscInt *bj,*cj;
2288: PetscInt bn,cn,j,maxentrycol = -1,maxdiffcol = -1,maxrdiffcol = -1;
2289: PetscReal maxentry = 0,maxdiff = 0,maxrdiff = 0;
2290: MatGetRow(*B,i,&bn,&bj,&ba);
2291: MatGetRow(Bfd,i,&cn,&cj,&ca);
2292: if (bn != cn) SETERRQ(((PetscObject)*A)->comm,PETSC_ERR_PLIB,"Unexpected different nonzero pattern in -snes_compare_coloring_threshold");
2293: for (j=0; j<bn; j++) {
2294: PetscReal rdiff = PetscAbsScalar(ca[j]) / (threshold_atol + threshold_rtol*PetscAbsScalar(ba[j]));
2295: if (PetscAbsScalar(ba[j]) > PetscAbs(maxentry)) {
2296: maxentrycol = bj[j];
2297: maxentry = PetscRealPart(ba[j]);
2298: }
2299: if (PetscAbsScalar(ca[j]) > PetscAbs(maxdiff)) {
2300: maxdiffcol = bj[j];
2301: maxdiff = PetscRealPart(ca[j]);
2302: }
2303: if (rdiff > maxrdiff) {
2304: maxrdiffcol = bj[j];
2305: maxrdiff = rdiff;
2306: }
2307: }
2308: if (maxrdiff > 1) {
2309: PetscViewerASCIIPrintf(vstdout,"row %D (maxentry=%G at %D, maxdiff=%G at %D, maxrdiff=%G at %D):",i,maxentry,maxentrycol,maxdiff,maxdiffcol,maxrdiff,maxrdiffcol);
2310: for (j=0; j<bn; j++) {
2311: PetscReal rdiff;
2312: rdiff = PetscAbsScalar(ca[j]) / (threshold_atol + threshold_rtol*PetscAbsScalar(ba[j]));
2313: if (rdiff > 1) {
2314: PetscViewerASCIIPrintf(vstdout," (%D,%G:%G)",bj[j],PetscRealPart(ba[j]),PetscRealPart(ca[j]));
2315: }
2316: }
2317: PetscViewerASCIIPrintf(vstdout,"\n",i,maxentry,maxdiff,maxrdiff);
2318: }
2319: MatRestoreRow(*B,i,&bn,&bj,&ba);
2320: MatRestoreRow(Bfd,i,&cn,&cj,&ca);
2321: }
2322: }
2323: PetscViewerDestroy(&vdraw);
2324: MatDestroy(&Bfd);
2325: }
2326: }
2327: return(0);
2328: }
2330: /*MC
2331: SNESJacobianFunction - function used to convey the nonlinear Jacobian of the function to be solved by SNES
2333: Synopsis:
2334: #include "petscsnes.h"
2335: $ SNESJacobianFunction(SNES snes,Vec x,Mat *Amat,Mat *Pmat,int *flag,void *ctx);
2337: + x - input vector
2338: . Amat - the matrix that defines the (approximate) Jacobian
2339: . Pmat - the matrix to be used in constructing the preconditioner, usually the same as Amat.
2340: . flag - flag indicating information about the preconditioner matrix
2341: structure (same as flag in KSPSetOperators()), one of SAME_NONZERO_PATTERN,DIFFERENT_NONZERO_PATTERN,SAME_PRECONDITIONER
2342: - ctx - [optional] user-defined Jacobian context
2344: Level: intermediate
2346: .seealso: SNESSetFunction(), SNESGetFunction(), SNESSetJacobian(), SNESGetJacobian()
2347: M*/
2351: /*@C
2352: SNESSetJacobian - Sets the function to compute Jacobian as well as the
2353: location to store the matrix.
2355: Logically Collective on SNES and Mat
2357: Input Parameters:
2358: + snes - the SNES context
2359: . Amat - the matrix that defines the (approximate) Jacobian
2360: . Pmat - the matrix to be used in constructing the preconditioner, usually the same as Amat.
2361: . SNESJacobianFunction - Jacobian evaluation routine (if NULL then SNES retains any previously set value)
2362: - ctx - [optional] user-defined context for private data for the
2363: Jacobian evaluation routine (may be NULL) (if NULL then SNES retains any previously set value)
2365: Notes:
2366: See KSPSetOperators() for important information about setting the flag
2367: output parameter in the routine func(). Be sure to read this information!
2369: The routine func() takes Mat * as the matrix arguments rather than Mat.
2370: This allows the Jacobian evaluation routine to replace A and/or B with a
2371: completely new new matrix structure (not just different matrix elements)
2372: when appropriate, for instance, if the nonzero structure is changing
2373: throughout the global iterations.
2375: If the Amat matrix and Pmat matrix are different you must call MatAssemblyBegin/End() on
2376: each matrix.
2378: If using SNESComputeJacobianDefaultColor() to assemble a Jacobian, the ctx argument
2379: must be a MatFDColoring.
2381: Other defect-correction schemes can be used by computing a different matrix in place of the Jacobian. One common
2382: example is to use the "Picard linearization" which only differentiates through the highest order parts of each term.
2384: Level: beginner
2386: .keywords: SNES, nonlinear, set, Jacobian, matrix
2388: .seealso: KSPSetOperators(), SNESSetFunction(), MatMFFDComputeJacobian(), SNESComputeJacobianDefaultColor(), MatStructure, SNESJacobianFunction, SNESSetPicard()
2389: @*/
2390: PetscErrorCode SNESSetJacobian(SNES snes,Mat Amat,Mat Pmat,PetscErrorCode (*SNESJacobianFunction)(SNES,Vec,Mat*,Mat*,MatStructure*,void*),void *ctx)
2391: {
2393: DM dm;
2401: SNESGetDM(snes,&dm);
2402: DMSNESSetJacobian(dm,SNESJacobianFunction,ctx);
2403: if (Amat) {
2404: PetscObjectReference((PetscObject)Amat);
2405: MatDestroy(&snes->jacobian);
2407: snes->jacobian = Amat;
2408: }
2409: if (Pmat) {
2410: PetscObjectReference((PetscObject)Pmat);
2411: MatDestroy(&snes->jacobian_pre);
2413: snes->jacobian_pre = Pmat;
2414: }
2415: return(0);
2416: }
2420: /*@C
2421: SNESGetJacobian - Returns the Jacobian matrix and optionally the user
2422: provided context for evaluating the Jacobian.
2424: Not Collective, but Mat object will be parallel if SNES object is
2426: Input Parameter:
2427: . snes - the nonlinear solver context
2429: Output Parameters:
2430: + Amat - location to stash (approximate) Jacobian matrix (or NULL)
2431: . Pmat - location to stash matrix used to compute the preconditioner (or NULL)
2432: . SNESJacobianFunction - location to put Jacobian function (or NULL)
2433: - ctx - location to stash Jacobian ctx (or NULL)
2435: Level: advanced
2437: .seealso: SNESSetJacobian(), SNESComputeJacobian()
2438: @*/
2439: PetscErrorCode SNESGetJacobian(SNES snes,Mat *Amat,Mat *Pmat,PetscErrorCode (**SNESJacobianFunction)(SNES,Vec,Mat*,Mat*,MatStructure*,void*),void **ctx)
2440: {
2442: DM dm;
2443: DMSNES sdm;
2447: if (Amat) *Amat = snes->jacobian;
2448: if (Pmat) *Pmat = snes->jacobian_pre;
2449: SNESGetDM(snes,&dm);
2450: DMGetDMSNES(dm,&sdm);
2451: if (SNESJacobianFunction) *SNESJacobianFunction = sdm->ops->computejacobian;
2452: if (ctx) *ctx = sdm->jacobianctx;
2453: return(0);
2454: }
2458: /*@
2459: SNESSetUp - Sets up the internal data structures for the later use
2460: of a nonlinear solver.
2462: Collective on SNES
2464: Input Parameters:
2465: . snes - the SNES context
2467: Notes:
2468: For basic use of the SNES solvers the user need not explicitly call
2469: SNESSetUp(), since these actions will automatically occur during
2470: the call to SNESSolve(). However, if one wishes to control this
2471: phase separately, SNESSetUp() should be called after SNESCreate()
2472: and optional routines of the form SNESSetXXX(), but before SNESSolve().
2474: Level: advanced
2476: .keywords: SNES, nonlinear, setup
2478: .seealso: SNESCreate(), SNESSolve(), SNESDestroy()
2479: @*/
2480: PetscErrorCode SNESSetUp(SNES snes)
2481: {
2483: DM dm;
2484: DMSNES sdm;
2485: SNESLineSearch linesearch, pclinesearch;
2486: void *lsprectx,*lspostctx;
2487: PetscErrorCode (*precheck)(SNESLineSearch,Vec,Vec,PetscBool*,void*);
2488: PetscErrorCode (*postcheck)(SNESLineSearch,Vec,Vec,Vec,PetscBool*,PetscBool*,void*);
2489: PetscErrorCode (*func)(SNES,Vec,Vec,void*);
2490: Vec f,fpc;
2491: void *funcctx;
2492: PetscErrorCode (*jac)(SNES,Vec,Mat*,Mat*,MatStructure*,void*);
2493: void *jacctx,*appctx;
2494: Mat A,B;
2498: if (snes->setupcalled) return(0);
2500: if (!((PetscObject)snes)->type_name) {
2501: SNESSetType(snes,SNESNEWTONLS);
2502: }
2504: SNESGetFunction(snes,&snes->vec_func,NULL,NULL);
2505: if (snes->vec_func == snes->vec_sol) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_IDN,"Solution vector cannot be function vector");
2506: if (snes->vec_rhs == snes->vec_sol) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_IDN,"Solution vector cannot be right hand side vector");
2508: if (!snes->vec_sol_update /* && snes->vec_sol */) {
2509: VecDuplicate(snes->vec_sol,&snes->vec_sol_update);
2510: PetscLogObjectParent(snes,snes->vec_sol_update);
2511: }
2513: SNESGetDM(snes,&dm);
2514: DMShellSetGlobalVector(snes->dm,snes->vec_sol);
2515: DMGetDMSNES(dm,&sdm);
2516: if (!sdm->ops->computefunction) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"Function never provided to SNES object");
2517: if (!sdm->ops->computejacobian) {
2518: DMSNESSetJacobian(dm,SNESComputeJacobianDefaultColor,NULL);
2519: }
2520: if (!snes->vec_func) {
2521: DMCreateGlobalVector(dm,&snes->vec_func);
2522: }
2524: if (!snes->ksp) {
2525: SNESGetKSP(snes, &snes->ksp);
2526: }
2528: if (!snes->linesearch) {
2529: SNESGetLineSearch(snes, &snes->linesearch);
2530: }
2532: if (snes->pc && (snes->pcside == PC_LEFT)) {
2533: snes->mf = PETSC_TRUE;
2534: snes->mf_operator = PETSC_FALSE;
2535: }
2537: if (snes->mf) {
2538: SNESSetUpMatrixFree_Private(snes, snes->mf_operator, snes->mf_version);
2539: }
2541: if (snes->ops->usercompute && !snes->user) {
2542: (*snes->ops->usercompute)(snes,(void**)&snes->user);
2543: }
2545: if (snes->pc) {
2546: /* copy the DM over */
2547: SNESGetDM(snes,&dm);
2548: SNESSetDM(snes->pc,dm);
2550: SNESGetFunction(snes,&f,&func,&funcctx);
2551: VecDuplicate(f,&fpc);
2552: SNESSetFunction(snes->pc,fpc,func,funcctx);
2553: SNESGetJacobian(snes,&A,&B,&jac,&jacctx);
2554: SNESSetJacobian(snes->pc,A,B,jac,jacctx);
2555: SNESGetApplicationContext(snes,&appctx);
2556: SNESSetApplicationContext(snes->pc,appctx);
2557: VecDestroy(&fpc);
2559: /* copy the function pointers over */
2560: PetscObjectCopyFortranFunctionPointers((PetscObject)snes,(PetscObject)snes->pc);
2562: /* default to 1 iteration */
2563: SNESSetTolerances(snes->pc,0.0,0.0,0.0,1,snes->pc->max_funcs);
2564: if (snes->pcside==PC_RIGHT) {
2565: SNESSetNormType(snes->pc,SNES_NORM_FINAL_ONLY);
2566: } else {
2567: SNESSetNormType(snes->pc,SNES_NORM_NONE);
2568: }
2569: SNESSetFromOptions(snes->pc);
2571: /* copy the line search context over */
2572: SNESGetLineSearch(snes,&linesearch);
2573: SNESGetLineSearch(snes->pc,&pclinesearch);
2574: SNESLineSearchGetPreCheck(linesearch,&precheck,&lsprectx);
2575: SNESLineSearchGetPostCheck(linesearch,&postcheck,&lspostctx);
2576: SNESLineSearchSetPreCheck(pclinesearch,precheck,lsprectx);
2577: SNESLineSearchSetPostCheck(pclinesearch,postcheck,lspostctx);
2578: PetscObjectCopyFortranFunctionPointers((PetscObject)linesearch, (PetscObject)pclinesearch);
2579: }
2581: if (snes->ops->setup) {
2582: (*snes->ops->setup)(snes);
2583: }
2585: snes->setupcalled = PETSC_TRUE;
2586: return(0);
2587: }
2591: /*@
2592: SNESReset - Resets a SNES context to the snessetupcalled = 0 state and removes any allocated Vecs and Mats
2594: Collective on SNES
2596: Input Parameter:
2597: . snes - iterative context obtained from SNESCreate()
2599: Level: intermediate
2601: Notes: Also calls the application context destroy routine set with SNESSetComputeApplicationContext()
2603: .keywords: SNES, destroy
2605: .seealso: SNESCreate(), SNESSetUp(), SNESSolve()
2606: @*/
2607: PetscErrorCode SNESReset(SNES snes)
2608: {
2613: if (snes->ops->userdestroy && snes->user) {
2614: (*snes->ops->userdestroy)((void**)&snes->user);
2615: snes->user = NULL;
2616: }
2617: if (snes->pc) {
2618: SNESReset(snes->pc);
2619: }
2621: if (snes->ops->reset) {
2622: (*snes->ops->reset)(snes);
2623: }
2624: if (snes->ksp) {
2625: KSPReset(snes->ksp);
2626: }
2628: if (snes->linesearch) {
2629: SNESLineSearchReset(snes->linesearch);
2630: }
2632: VecDestroy(&snes->vec_rhs);
2633: VecDestroy(&snes->vec_sol);
2634: VecDestroy(&snes->vec_sol_update);
2635: VecDestroy(&snes->vec_func);
2636: MatDestroy(&snes->jacobian);
2637: MatDestroy(&snes->jacobian_pre);
2638: VecDestroyVecs(snes->nwork,&snes->work);
2639: VecDestroyVecs(snes->nvwork,&snes->vwork);
2641: snes->nwork = snes->nvwork = 0;
2642: snes->setupcalled = PETSC_FALSE;
2643: return(0);
2644: }
2648: /*@
2649: SNESDestroy - Destroys the nonlinear solver context that was created
2650: with SNESCreate().
2652: Collective on SNES
2654: Input Parameter:
2655: . snes - the SNES context
2657: Level: beginner
2659: .keywords: SNES, nonlinear, destroy
2661: .seealso: SNESCreate(), SNESSolve()
2662: @*/
2663: PetscErrorCode SNESDestroy(SNES *snes)
2664: {
2668: if (!*snes) return(0);
2670: if (--((PetscObject)(*snes))->refct > 0) {*snes = 0; return(0);}
2672: SNESReset((*snes));
2673: SNESDestroy(&(*snes)->pc);
2675: /* if memory was published with AMS then destroy it */
2676: PetscObjectAMSViewOff((PetscObject)*snes);
2677: if ((*snes)->ops->destroy) {(*((*snes))->ops->destroy)((*snes));}
2679: DMDestroy(&(*snes)->dm);
2680: KSPDestroy(&(*snes)->ksp);
2681: SNESLineSearchDestroy(&(*snes)->linesearch);
2683: PetscFree((*snes)->kspconvctx);
2684: if ((*snes)->ops->convergeddestroy) {
2685: (*(*snes)->ops->convergeddestroy)((*snes)->cnvP);
2686: }
2687: if ((*snes)->conv_malloc) {
2688: PetscFree((*snes)->conv_hist);
2689: PetscFree((*snes)->conv_hist_its);
2690: }
2691: SNESMonitorCancel((*snes));
2692: PetscHeaderDestroy(snes);
2693: return(0);
2694: }
2696: /* ----------- Routines to set solver parameters ---------- */
2700: /*@
2701: SNESSetLagPreconditioner - Determines when the preconditioner is rebuilt in the nonlinear solve.
2703: Logically Collective on SNES
2705: Input Parameters:
2706: + snes - the SNES context
2707: - lag - -1 indicates NEVER rebuild, 1 means rebuild every time the Jacobian is computed within a single nonlinear solve, 2 means every second time
2708: the Jacobian is built etc. -2 indicates rebuild preconditioner at next chance but then never rebuild after that
2710: Options Database Keys:
2711: . -snes_lag_preconditioner <lag>
2713: Notes:
2714: The default is 1
2715: The preconditioner is ALWAYS built in the first iteration of a nonlinear solve unless lag is -1
2716: If -1 is used before the very first nonlinear solve the preconditioner is still built because there is no previous preconditioner to use
2718: Level: intermediate
2720: .keywords: SNES, nonlinear, set, convergence, tolerances
2722: .seealso: SNESSetTrustRegionTolerance(), SNESGetLagPreconditioner(), SNESSetLagJacobian(), SNESGetLagJacobian()
2724: @*/
2725: PetscErrorCode SNESSetLagPreconditioner(SNES snes,PetscInt lag)
2726: {
2729: if (lag < -2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Lag must be -2, -1, 1 or greater");
2730: if (!lag) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Lag cannot be 0");
2732: snes->lagpreconditioner = lag;
2733: return(0);
2734: }
2738: /*@
2739: SNESSetGridSequence - sets the number of steps of grid sequencing that SNES does
2741: Logically Collective on SNES
2743: Input Parameters:
2744: + snes - the SNES context
2745: - steps - the number of refinements to do, defaults to 0
2747: Options Database Keys:
2748: . -snes_grid_sequence <steps>
2750: Level: intermediate
2752: Notes:
2753: Use SNESGetSolution() to extract the fine grid solution after grid sequencing.
2755: .keywords: SNES, nonlinear, set, convergence, tolerances
2757: .seealso: SNESSetTrustRegionTolerance(), SNESGetLagPreconditioner(), SNESSetLagJacobian(), SNESGetLagJacobian()
2759: @*/
2760: PetscErrorCode SNESSetGridSequence(SNES snes,PetscInt steps)
2761: {
2765: snes->gridsequence = steps;
2766: return(0);
2767: }
2771: /*@
2772: SNESGetLagPreconditioner - Indicates how often the preconditioner is rebuilt
2774: Not Collective
2776: Input Parameter:
2777: . snes - the SNES context
2779: Output Parameter:
2780: . lag - -1 indicates NEVER rebuild, 1 means rebuild every time the Jacobian is computed within a single nonlinear solve, 2 means every second time
2781: the Jacobian is built etc. -2 indicates rebuild preconditioner at next chance but then never rebuild after that
2783: Options Database Keys:
2784: . -snes_lag_preconditioner <lag>
2786: Notes:
2787: The default is 1
2788: The preconditioner is ALWAYS built in the first iteration of a nonlinear solve unless lag is -1
2790: Level: intermediate
2792: .keywords: SNES, nonlinear, set, convergence, tolerances
2794: .seealso: SNESSetTrustRegionTolerance(), SNESSetLagPreconditioner()
2796: @*/
2797: PetscErrorCode SNESGetLagPreconditioner(SNES snes,PetscInt *lag)
2798: {
2801: *lag = snes->lagpreconditioner;
2802: return(0);
2803: }
2807: /*@
2808: SNESSetLagJacobian - Determines when the Jacobian is rebuilt in the nonlinear solve. See SNESSetLagPreconditioner() for determining how
2809: often the preconditioner is rebuilt.
2811: Logically Collective on SNES
2813: Input Parameters:
2814: + snes - the SNES context
2815: - lag - -1 indicates NEVER rebuild, 1 means rebuild every time the Jacobian is computed within a single nonlinear solve, 2 means every second time
2816: the Jacobian is built etc. -2 means rebuild at next chance but then never again
2818: Options Database Keys:
2819: . -snes_lag_jacobian <lag>
2821: Notes:
2822: The default is 1
2823: The Jacobian is ALWAYS built in the first iteration of a nonlinear solve unless lag is -1
2824: 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
2825: at the next Newton step but never again (unless it is reset to another value)
2827: Level: intermediate
2829: .keywords: SNES, nonlinear, set, convergence, tolerances
2831: .seealso: SNESSetTrustRegionTolerance(), SNESGetLagPreconditioner(), SNESSetLagPreconditioner(), SNESGetLagJacobian()
2833: @*/
2834: PetscErrorCode SNESSetLagJacobian(SNES snes,PetscInt lag)
2835: {
2838: if (lag < -2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Lag must be -2, -1, 1 or greater");
2839: if (!lag) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Lag cannot be 0");
2841: snes->lagjacobian = lag;
2842: return(0);
2843: }
2847: /*@
2848: SNESGetLagJacobian - Indicates how often the Jacobian is rebuilt. See SNESGetLagPreconditioner() to determine when the preconditioner is rebuilt
2850: Not Collective
2852: Input Parameter:
2853: . snes - the SNES context
2855: Output Parameter:
2856: . lag - -1 indicates NEVER rebuild, 1 means rebuild every time the Jacobian is computed within a single nonlinear solve, 2 means every second time
2857: the Jacobian is built etc.
2859: Options Database Keys:
2860: . -snes_lag_jacobian <lag>
2862: Notes:
2863: The default is 1
2864: The jacobian is ALWAYS built in the first iteration of a nonlinear solve unless lag is -1
2866: Level: intermediate
2868: .keywords: SNES, nonlinear, set, convergence, tolerances
2870: .seealso: SNESSetTrustRegionTolerance(), SNESSetLagJacobian(), SNESSetLagPreconditioner(), SNESGetLagPreconditioner()
2872: @*/
2873: PetscErrorCode SNESGetLagJacobian(SNES snes,PetscInt *lag)
2874: {
2877: *lag = snes->lagjacobian;
2878: return(0);
2879: }
2883: /*@
2884: SNESSetTolerances - Sets various parameters used in convergence tests.
2886: Logically Collective on SNES
2888: Input Parameters:
2889: + snes - the SNES context
2890: . abstol - absolute convergence tolerance
2891: . rtol - relative convergence tolerance
2892: . stol - convergence tolerance in terms of the norm of the change in the solution between steps, || delta x || < stol*|| x ||
2893: . maxit - maximum number of iterations
2894: - maxf - maximum number of function evaluations
2896: Options Database Keys:
2897: + -snes_atol <abstol> - Sets abstol
2898: . -snes_rtol <rtol> - Sets rtol
2899: . -snes_stol <stol> - Sets stol
2900: . -snes_max_it <maxit> - Sets maxit
2901: - -snes_max_funcs <maxf> - Sets maxf
2903: Notes:
2904: The default maximum number of iterations is 50.
2905: The default maximum number of function evaluations is 1000.
2907: Level: intermediate
2909: .keywords: SNES, nonlinear, set, convergence, tolerances
2911: .seealso: SNESSetTrustRegionTolerance()
2912: @*/
2913: PetscErrorCode SNESSetTolerances(SNES snes,PetscReal abstol,PetscReal rtol,PetscReal stol,PetscInt maxit,PetscInt maxf)
2914: {
2923: if (abstol != PETSC_DEFAULT) {
2924: if (abstol < 0.0) SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_OUTOFRANGE,"Absolute tolerance %G must be non-negative",abstol);
2925: snes->abstol = abstol;
2926: }
2927: if (rtol != PETSC_DEFAULT) {
2928: 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",rtol);
2929: snes->rtol = rtol;
2930: }
2931: if (stol != PETSC_DEFAULT) {
2932: if (stol < 0.0) SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_OUTOFRANGE,"Step tolerance %G must be non-negative",stol);
2933: snes->stol = stol;
2934: }
2935: if (maxit != PETSC_DEFAULT) {
2936: if (maxit < 0) SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_OUTOFRANGE,"Maximum number of iterations %D must be non-negative",maxit);
2937: snes->max_its = maxit;
2938: }
2939: if (maxf != PETSC_DEFAULT) {
2940: if (maxf < 0) SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_OUTOFRANGE,"Maximum number of function evaluations %D must be non-negative",maxf);
2941: snes->max_funcs = maxf;
2942: }
2943: snes->tolerancesset = PETSC_TRUE;
2944: return(0);
2945: }
2949: /*@
2950: SNESGetTolerances - Gets various parameters used in convergence tests.
2952: Not Collective
2954: Input Parameters:
2955: + snes - the SNES context
2956: . atol - absolute convergence tolerance
2957: . rtol - relative convergence tolerance
2958: . stol - convergence tolerance in terms of the norm
2959: of the change in the solution between steps
2960: . maxit - maximum number of iterations
2961: - maxf - maximum number of function evaluations
2963: Notes:
2964: The user can specify NULL for any parameter that is not needed.
2966: Level: intermediate
2968: .keywords: SNES, nonlinear, get, convergence, tolerances
2970: .seealso: SNESSetTolerances()
2971: @*/
2972: PetscErrorCode SNESGetTolerances(SNES snes,PetscReal *atol,PetscReal *rtol,PetscReal *stol,PetscInt *maxit,PetscInt *maxf)
2973: {
2976: if (atol) *atol = snes->abstol;
2977: if (rtol) *rtol = snes->rtol;
2978: if (stol) *stol = snes->stol;
2979: if (maxit) *maxit = snes->max_its;
2980: if (maxf) *maxf = snes->max_funcs;
2981: return(0);
2982: }
2986: /*@
2987: SNESSetTrustRegionTolerance - Sets the trust region parameter tolerance.
2989: Logically Collective on SNES
2991: Input Parameters:
2992: + snes - the SNES context
2993: - tol - tolerance
2995: Options Database Key:
2996: . -snes_trtol <tol> - Sets tol
2998: Level: intermediate
3000: .keywords: SNES, nonlinear, set, trust region, tolerance
3002: .seealso: SNESSetTolerances()
3003: @*/
3004: PetscErrorCode SNESSetTrustRegionTolerance(SNES snes,PetscReal tol)
3005: {
3009: snes->deltatol = tol;
3010: return(0);
3011: }
3013: /*
3014: Duplicate the lg monitors for SNES from KSP; for some reason with
3015: dynamic libraries things don't work under Sun4 if we just use
3016: macros instead of functions
3017: */
3020: PetscErrorCode SNESMonitorLGResidualNorm(SNES snes,PetscInt it,PetscReal norm,void *ctx)
3021: {
3026: KSPMonitorLGResidualNorm((KSP)snes,it,norm,ctx);
3027: return(0);
3028: }
3032: PetscErrorCode SNESMonitorLGCreate(const char host[],const char label[],int x,int y,int m,int n,PetscDrawLG *draw)
3033: {
3037: KSPMonitorLGResidualNormCreate(host,label,x,y,m,n,draw);
3038: return(0);
3039: }
3043: PetscErrorCode SNESMonitorLGDestroy(PetscDrawLG *draw)
3044: {
3048: KSPMonitorLGResidualNormDestroy(draw);
3049: return(0);
3050: }
3052: extern PetscErrorCode SNESMonitorRange_Private(SNES,PetscInt,PetscReal*);
3055: PetscErrorCode SNESMonitorLGRange(SNES snes,PetscInt n,PetscReal rnorm,void *monctx)
3056: {
3057: PetscDrawLG lg;
3058: PetscErrorCode ierr;
3059: PetscReal x,y,per;
3060: PetscViewer v = (PetscViewer)monctx;
3061: static PetscReal prev; /* should be in the context */
3062: PetscDraw draw;
3065: PetscViewerDrawGetDrawLG(v,0,&lg);
3066: if (!n) {PetscDrawLGReset(lg);}
3067: PetscDrawLGGetDraw(lg,&draw);
3068: PetscDrawSetTitle(draw,"Residual norm");
3069: x = (PetscReal)n;
3070: if (rnorm > 0.0) y = log10(rnorm);
3071: else y = -15.0;
3072: PetscDrawLGAddPoint(lg,&x,&y);
3073: if (n < 20 || !(n % 5)) {
3074: PetscDrawLGDraw(lg);
3075: }
3077: PetscViewerDrawGetDrawLG(v,1,&lg);
3078: if (!n) {PetscDrawLGReset(lg);}
3079: PetscDrawLGGetDraw(lg,&draw);
3080: PetscDrawSetTitle(draw,"% elemts > .2*max elemt");
3081: SNESMonitorRange_Private(snes,n,&per);
3082: x = (PetscReal)n;
3083: y = 100.0*per;
3084: PetscDrawLGAddPoint(lg,&x,&y);
3085: if (n < 20 || !(n % 5)) {
3086: PetscDrawLGDraw(lg);
3087: }
3089: PetscViewerDrawGetDrawLG(v,2,&lg);
3090: if (!n) {prev = rnorm;PetscDrawLGReset(lg);}
3091: PetscDrawLGGetDraw(lg,&draw);
3092: PetscDrawSetTitle(draw,"(norm -oldnorm)/oldnorm");
3093: x = (PetscReal)n;
3094: y = (prev - rnorm)/prev;
3095: PetscDrawLGAddPoint(lg,&x,&y);
3096: if (n < 20 || !(n % 5)) {
3097: PetscDrawLGDraw(lg);
3098: }
3100: PetscViewerDrawGetDrawLG(v,3,&lg);
3101: if (!n) {PetscDrawLGReset(lg);}
3102: PetscDrawLGGetDraw(lg,&draw);
3103: PetscDrawSetTitle(draw,"(norm -oldnorm)/oldnorm*(% > .2 max)");
3104: x = (PetscReal)n;
3105: y = (prev - rnorm)/(prev*per);
3106: if (n > 2) { /*skip initial crazy value */
3107: PetscDrawLGAddPoint(lg,&x,&y);
3108: }
3109: if (n < 20 || !(n % 5)) {
3110: PetscDrawLGDraw(lg);
3111: }
3112: prev = rnorm;
3113: return(0);
3114: }
3118: /*@
3119: SNESMonitor - runs the user provided monitor routines, if they exist
3121: Collective on SNES
3123: Input Parameters:
3124: + snes - nonlinear solver context obtained from SNESCreate()
3125: . iter - iteration number
3126: - rnorm - relative norm of the residual
3128: Notes:
3129: This routine is called by the SNES implementations.
3130: It does not typically need to be called by the user.
3132: Level: developer
3134: .seealso: SNESMonitorSet()
3135: @*/
3136: PetscErrorCode SNESMonitor(SNES snes,PetscInt iter,PetscReal rnorm)
3137: {
3139: PetscInt i,n = snes->numbermonitors;
3142: for (i=0; i<n; i++) {
3143: (*snes->monitor[i])(snes,iter,rnorm,snes->monitorcontext[i]);
3144: }
3145: return(0);
3146: }
3148: /* ------------ Routines to set performance monitoring options ----------- */
3150: /*MC
3151: SNESMonitorFunction - functional form passed to SNESMonitorSet() to monitor convergence of nonlinear solver
3153: Synopsis:
3154: #include "petscsnes.h"
3155: $ PetscErrorCode SNESMonitorFunction(SNES snes,PetscInt its, PetscReal norm,void *mctx)
3157: + snes - the SNES context
3158: . its - iteration number
3159: . norm - 2-norm function value (may be estimated)
3160: - mctx - [optional] monitoring context
3162: Level: advanced
3164: .seealso: SNESMonitorSet(), SNESMonitorGet()
3165: M*/
3169: /*@C
3170: SNESMonitorSet - Sets an ADDITIONAL function that is to be used at every
3171: iteration of the nonlinear solver to display the iteration's
3172: progress.
3174: Logically Collective on SNES
3176: Input Parameters:
3177: + snes - the SNES context
3178: . SNESMonitorFunction - monitoring routine
3179: . mctx - [optional] user-defined context for private data for the
3180: monitor routine (use NULL if no context is desired)
3181: - monitordestroy - [optional] routine that frees monitor context
3182: (may be NULL)
3184: Options Database Keys:
3185: + -snes_monitor - sets SNESMonitorDefault()
3186: . -snes_monitor_lg_residualnorm - sets line graph monitor,
3187: uses SNESMonitorLGCreate()
3188: - -snes_monitor_cancel - cancels all monitors that have
3189: been hardwired into a code by
3190: calls to SNESMonitorSet(), but
3191: does not cancel those set via
3192: the options database.
3194: Notes:
3195: Several different monitoring routines may be set by calling
3196: SNESMonitorSet() multiple times; all will be called in the
3197: order in which they were set.
3199: Fortran notes: Only a single monitor function can be set for each SNES object
3201: Level: intermediate
3203: .keywords: SNES, nonlinear, set, monitor
3205: .seealso: SNESMonitorDefault(), SNESMonitorCancel(), SNESMonitorFunction
3206: @*/
3207: PetscErrorCode SNESMonitorSet(SNES snes,PetscErrorCode (*SNESMonitorFunction)(SNES,PetscInt,PetscReal,void*),void *mctx,PetscErrorCode (*monitordestroy)(void**))
3208: {
3209: PetscInt i;
3214: if (snes->numbermonitors >= MAXSNESMONITORS) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Too many monitors set");
3215: for (i=0; i<snes->numbermonitors;i++) {
3216: if (SNESMonitorFunction == snes->monitor[i] && monitordestroy == snes->monitordestroy[i] && mctx == snes->monitorcontext[i]) {
3217: if (monitordestroy) {
3218: (*monitordestroy)(&mctx);
3219: }
3220: return(0);
3221: }
3222: }
3223: snes->monitor[snes->numbermonitors] = SNESMonitorFunction;
3224: snes->monitordestroy[snes->numbermonitors] = monitordestroy;
3225: snes->monitorcontext[snes->numbermonitors++] = (void*)mctx;
3226: return(0);
3227: }
3231: /*@
3232: SNESMonitorCancel - Clears all the monitor functions for a SNES object.
3234: Logically Collective on SNES
3236: Input Parameters:
3237: . snes - the SNES context
3239: Options Database Key:
3240: . -snes_monitor_cancel - cancels all monitors that have been hardwired
3241: into a code by calls to SNESMonitorSet(), but does not cancel those
3242: set via the options database
3244: Notes:
3245: There is no way to clear one specific monitor from a SNES object.
3247: Level: intermediate
3249: .keywords: SNES, nonlinear, set, monitor
3251: .seealso: SNESMonitorDefault(), SNESMonitorSet()
3252: @*/
3253: PetscErrorCode SNESMonitorCancel(SNES snes)
3254: {
3256: PetscInt i;
3260: for (i=0; i<snes->numbermonitors; i++) {
3261: if (snes->monitordestroy[i]) {
3262: (*snes->monitordestroy[i])(&snes->monitorcontext[i]);
3263: }
3264: }
3265: snes->numbermonitors = 0;
3266: return(0);
3267: }
3269: /*MC
3270: SNESConvergenceTestFunction - functional form used for testing of convergence of nonlinear solver
3272: Synopsis:
3273: #include "petscsnes.h"
3274: $ PetscErrorCode SNESConvergenceTest(SNES snes,PetscInt it,PetscReal xnorm,PetscReal gnorm,PetscReal f,SNESConvergedReason *reason,void *cctx)
3276: + snes - the SNES context
3277: . it - current iteration (0 is the first and is before any Newton step)
3278: . cctx - [optional] convergence context
3279: . reason - reason for convergence/divergence
3280: . xnorm - 2-norm of current iterate
3281: . gnorm - 2-norm of current step
3282: - f - 2-norm of function
3284: Level: intermediate
3286: .seealso: SNESSetConvergenceTest(), SNESGetConvergenceTest()
3287: M*/
3291: /*@C
3292: SNESSetConvergenceTest - Sets the function that is to be used
3293: to test for convergence of the nonlinear iterative solution.
3295: Logically Collective on SNES
3297: Input Parameters:
3298: + snes - the SNES context
3299: . SNESConvergenceTestFunction - routine to test for convergence
3300: . cctx - [optional] context for private data for the convergence routine (may be NULL)
3301: - destroy - [optional] destructor for the context (may be NULL; NULL_FUNCTION in Fortran)
3303: Level: advanced
3305: .keywords: SNES, nonlinear, set, convergence, test
3307: .seealso: SNESConvergedDefault(), SNESSkipConverged(), SNESConvergenceTestFunction
3308: @*/
3309: PetscErrorCode SNESSetConvergenceTest(SNES snes,PetscErrorCode (*SNESConvergenceTestFunction)(SNES,PetscInt,PetscReal,PetscReal,PetscReal,SNESConvergedReason*,void*),void *cctx,PetscErrorCode (*destroy)(void*))
3310: {
3315: if (!SNESConvergenceTestFunction) SNESConvergenceTestFunction = SNESSkipConverged;
3316: if (snes->ops->convergeddestroy) {
3317: (*snes->ops->convergeddestroy)(snes->cnvP);
3318: }
3319: snes->ops->converged = SNESConvergenceTestFunction;
3320: snes->ops->convergeddestroy = destroy;
3321: snes->cnvP = cctx;
3322: return(0);
3323: }
3327: /*@
3328: SNESGetConvergedReason - Gets the reason the SNES iteration was stopped.
3330: Not Collective
3332: Input Parameter:
3333: . snes - the SNES context
3335: Output Parameter:
3336: . reason - negative value indicates diverged, positive value converged, see SNESConvergedReason or the
3337: manual pages for the individual convergence tests for complete lists
3339: Level: intermediate
3341: Notes: Can only be called after the call the SNESSolve() is complete.
3343: .keywords: SNES, nonlinear, set, convergence, test
3345: .seealso: SNESSetConvergenceTest(), SNESConvergedReason
3346: @*/
3347: PetscErrorCode SNESGetConvergedReason(SNES snes,SNESConvergedReason *reason)
3348: {
3352: *reason = snes->reason;
3353: return(0);
3354: }
3358: /*@
3359: SNESSetConvergenceHistory - Sets the array used to hold the convergence history.
3361: Logically Collective on SNES
3363: Input Parameters:
3364: + snes - iterative context obtained from SNESCreate()
3365: . a - array to hold history, this array will contain the function norms computed at each step
3366: . its - integer array holds the number of linear iterations for each solve.
3367: . na - size of a and its
3368: - reset - PETSC_TRUE indicates each new nonlinear solve resets the history counter to zero,
3369: else it continues storing new values for new nonlinear solves after the old ones
3371: Notes:
3372: If 'a' and 'its' are NULL then space is allocated for the history. If 'na' PETSC_DECIDE or PETSC_DEFAULT then a
3373: default array of length 10000 is allocated.
3375: This routine is useful, e.g., when running a code for purposes
3376: of accurate performance monitoring, when no I/O should be done
3377: during the section of code that is being timed.
3379: Level: intermediate
3381: .keywords: SNES, set, convergence, history
3383: .seealso: SNESGetConvergenceHistory()
3385: @*/
3386: PetscErrorCode SNESSetConvergenceHistory(SNES snes,PetscReal a[],PetscInt its[],PetscInt na,PetscBool reset)
3387: {
3394: if (!a) {
3395: if (na == PETSC_DECIDE || na == PETSC_DEFAULT) na = 1000;
3396: PetscMalloc(na*sizeof(PetscReal),&a);
3397: PetscMalloc(na*sizeof(PetscInt),&its);
3399: snes->conv_malloc = PETSC_TRUE;
3400: }
3401: snes->conv_hist = a;
3402: snes->conv_hist_its = its;
3403: snes->conv_hist_max = na;
3404: snes->conv_hist_len = 0;
3405: snes->conv_hist_reset = reset;
3406: return(0);
3407: }
3409: #if defined(PETSC_HAVE_MATLAB_ENGINE)
3410: #include <engine.h> /* MATLAB include file */
3411: #include <mex.h> /* MATLAB include file */
3415: PETSC_EXTERN mxArray *SNESGetConvergenceHistoryMatlab(SNES snes)
3416: {
3417: mxArray *mat;
3418: PetscInt i;
3419: PetscReal *ar;
3422: mat = mxCreateDoubleMatrix(snes->conv_hist_len,1,mxREAL);
3423: ar = (PetscReal*) mxGetData(mat);
3424: for (i=0; i<snes->conv_hist_len; i++) ar[i] = snes->conv_hist[i];
3425: PetscFunctionReturn(mat);
3426: }
3427: #endif
3431: /*@C
3432: SNESGetConvergenceHistory - Gets the array used to hold the convergence history.
3434: Not Collective
3436: Input Parameter:
3437: . snes - iterative context obtained from SNESCreate()
3439: Output Parameters:
3440: . a - array to hold history
3441: . its - integer array holds the number of linear iterations (or
3442: negative if not converged) for each solve.
3443: - na - size of a and its
3445: Notes:
3446: The calling sequence for this routine in Fortran is
3447: $ call SNESGetConvergenceHistory(SNES snes, integer na, integer ierr)
3449: This routine is useful, e.g., when running a code for purposes
3450: of accurate performance monitoring, when no I/O should be done
3451: during the section of code that is being timed.
3453: Level: intermediate
3455: .keywords: SNES, get, convergence, history
3457: .seealso: SNESSetConvergencHistory()
3459: @*/
3460: PetscErrorCode SNESGetConvergenceHistory(SNES snes,PetscReal *a[],PetscInt *its[],PetscInt *na)
3461: {
3464: if (a) *a = snes->conv_hist;
3465: if (its) *its = snes->conv_hist_its;
3466: if (na) *na = snes->conv_hist_len;
3467: return(0);
3468: }
3472: /*@C
3473: SNESSetUpdate - Sets the general-purpose update function called
3474: at the beginning of every iteration of the nonlinear solve. Specifically
3475: it is called just before the Jacobian is "evaluated".
3477: Logically Collective on SNES
3479: Input Parameters:
3480: . snes - The nonlinear solver context
3481: . func - The function
3483: Calling sequence of func:
3484: . func (SNES snes, PetscInt step);
3486: . step - The current step of the iteration
3488: Level: advanced
3490: 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()
3491: This is not used by most users.
3493: .keywords: SNES, update
3495: .seealso SNESSetJacobian(), SNESSolve()
3496: @*/
3497: PetscErrorCode SNESSetUpdate(SNES snes, PetscErrorCode (*func)(SNES, PetscInt))
3498: {
3501: snes->ops->update = func;
3502: return(0);
3503: }
3507: /*
3508: SNESScaleStep_Private - Scales a step so that its length is less than the
3509: positive parameter delta.
3511: Input Parameters:
3512: + snes - the SNES context
3513: . y - approximate solution of linear system
3514: . fnorm - 2-norm of current function
3515: - delta - trust region size
3517: Output Parameters:
3518: + gpnorm - predicted function norm at the new point, assuming local
3519: linearization. The value is zero if the step lies within the trust
3520: region, and exceeds zero otherwise.
3521: - ynorm - 2-norm of the step
3523: Note:
3524: For non-trust region methods such as SNESNEWTONLS, the parameter delta
3525: is set to be the maximum allowable step size.
3527: .keywords: SNES, nonlinear, scale, step
3528: */
3529: PetscErrorCode SNESScaleStep_Private(SNES snes,Vec y,PetscReal *fnorm,PetscReal *delta,PetscReal *gpnorm,PetscReal *ynorm)
3530: {
3531: PetscReal nrm;
3532: PetscScalar cnorm;
3540: VecNorm(y,NORM_2,&nrm);
3541: if (nrm > *delta) {
3542: nrm = *delta/nrm;
3543: *gpnorm = (1.0 - nrm)*(*fnorm);
3544: cnorm = nrm;
3545: VecScale(y,cnorm);
3546: *ynorm = *delta;
3547: } else {
3548: *gpnorm = 0.0;
3549: *ynorm = nrm;
3550: }
3551: return(0);
3552: }
3556: /*@C
3557: SNESSolve - Solves a nonlinear system F(x) = b.
3558: Call SNESSolve() after calling SNESCreate() and optional routines of the form SNESSetXXX().
3560: Collective on SNES
3562: Input Parameters:
3563: + snes - the SNES context
3564: . b - the constant part of the equation F(x) = b, or NULL to use zero.
3565: - x - the solution vector.
3567: Notes:
3568: The user should initialize the vector,x, with the initial guess
3569: for the nonlinear solve prior to calling SNESSolve. In particular,
3570: to employ an initial guess of zero, the user should explicitly set
3571: this vector to zero by calling VecSet().
3573: Level: beginner
3575: .keywords: SNES, nonlinear, solve
3577: .seealso: SNESCreate(), SNESDestroy(), SNESSetFunction(), SNESSetJacobian(), SNESSetGridSequence(), SNESGetSolution()
3578: @*/
3579: PetscErrorCode SNESSolve(SNES snes,Vec b,Vec x)
3580: {
3581: PetscErrorCode ierr;
3582: PetscBool flg;
3583: PetscViewer viewer;
3584: PetscInt grid;
3585: Vec xcreated = NULL;
3586: DM dm;
3587: PetscViewerFormat format;
3596: if (!x) {
3597: SNESGetDM(snes,&dm);
3598: DMCreateGlobalVector(dm,&xcreated);
3599: x = xcreated;
3600: }
3602: PetscOptionsGetViewer(PetscObjectComm((PetscObject)snes),((PetscObject)snes)->prefix,"-snes_view_pre",&viewer,&format,&flg);
3603: if (flg && !PetscPreLoadingOn) {
3604: PetscViewerPushFormat(viewer,format);
3605: SNESView(snes,viewer);
3606: PetscViewerPopFormat(viewer);
3607: PetscViewerDestroy(&viewer);
3608: }
3610: for (grid=0; grid<snes->gridsequence; grid++) {PetscViewerASCIIPushTab(PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)snes)));}
3611: for (grid=0; grid<snes->gridsequence+1; grid++) {
3613: /* set solution vector */
3614: if (!grid) {PetscObjectReference((PetscObject)x);}
3615: VecDestroy(&snes->vec_sol);
3616: snes->vec_sol = x;
3617: SNESGetDM(snes,&dm);
3619: /* set affine vector if provided */
3620: if (b) { PetscObjectReference((PetscObject)b); }
3621: VecDestroy(&snes->vec_rhs);
3622: snes->vec_rhs = b;
3624: SNESSetUp(snes);
3626: if (!grid) {
3627: if (snes->ops->computeinitialguess) {
3628: (*snes->ops->computeinitialguess)(snes,snes->vec_sol,snes->initialguessP);
3629: }
3630: }
3632: if (snes->conv_hist_reset) snes->conv_hist_len = 0;
3633: snes->nfuncs = 0; snes->linear_its = 0; snes->numFailures = 0;
3635: PetscLogEventBegin(SNES_Solve,snes,0,0,0);
3636: (*snes->ops->solve)(snes);
3637: PetscLogEventEnd(SNES_Solve,snes,0,0,0);
3638: if (snes->domainerror) {
3639: snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN;
3640: snes->domainerror = PETSC_FALSE;
3641: }
3642: if (!snes->reason) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Internal error, solver returned without setting converged reason");
3644: flg = PETSC_FALSE;
3645: PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_test_local_min",&flg,NULL);
3646: if (flg && !PetscPreLoadingOn) { SNESTestLocalMin(snes); }
3647: if (snes->printreason) {
3648: PetscViewerASCIIAddTab(PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)snes)),((PetscObject)snes)->tablevel);
3649: if (snes->reason > 0) {
3650: PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)snes)),"Nonlinear solve converged due to %s iterations %D\n",SNESConvergedReasons[snes->reason],snes->iter);
3651: } else {
3652: PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)snes)),"Nonlinear solve did not converge due to %s iterations %D\n",SNESConvergedReasons[snes->reason],snes->iter);
3653: }
3654: PetscViewerASCIISubtractTab(PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)snes)),((PetscObject)snes)->tablevel);
3655: }
3657: if (snes->errorifnotconverged && snes->reason < 0) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_NOT_CONVERGED,"SNESSolve has not converged");
3658: if (grid < snes->gridsequence) {
3659: DM fine;
3660: Vec xnew;
3661: Mat interp;
3663: DMRefine(snes->dm,PetscObjectComm((PetscObject)snes),&fine);
3664: if (!fine) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_INCOMP,"DMRefine() did not perform any refinement, cannot continue grid sequencing");
3665: DMCreateInterpolation(snes->dm,fine,&interp,NULL);
3666: DMCreateGlobalVector(fine,&xnew);
3667: MatInterpolate(interp,x,xnew);
3668: DMInterpolate(snes->dm,interp,fine);
3669: MatDestroy(&interp);
3670: x = xnew;
3672: SNESReset(snes);
3673: SNESSetDM(snes,fine);
3674: DMDestroy(&fine);
3675: PetscViewerASCIIPopTab(PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)snes)));
3676: }
3677: }
3678: /* monitoring and viewing */
3679: PetscOptionsGetViewer(PetscObjectComm((PetscObject)snes),((PetscObject)snes)->prefix,"-snes_view",&viewer,&format,&flg);
3680: if (flg && !PetscPreLoadingOn) {
3681: PetscViewerPushFormat(viewer,format);
3682: SNESView(snes,viewer);
3683: PetscViewerPopFormat(viewer);
3684: PetscViewerDestroy(&viewer);
3685: }
3686: VecViewFromOptions(snes->vec_sol,((PetscObject)snes)->prefix,"-snes_view_solution");
3688: VecDestroy(&xcreated);
3689: PetscObjectAMSBlock((PetscObject)snes);
3690: return(0);
3691: }
3693: /* --------- Internal routines for SNES Package --------- */
3697: /*@C
3698: SNESSetType - Sets the method for the nonlinear solver.
3700: Collective on SNES
3702: Input Parameters:
3703: + snes - the SNES context
3704: - type - a known method
3706: Options Database Key:
3707: . -snes_type <type> - Sets the method; use -help for a list
3708: of available methods (for instance, newtonls or newtontr)
3710: Notes:
3711: See "petsc/include/petscsnes.h" for available methods (for instance)
3712: + SNESNEWTONLS - Newton's method with line search
3713: (systems of nonlinear equations)
3714: . SNESNEWTONTR - Newton's method with trust region
3715: (systems of nonlinear equations)
3717: Normally, it is best to use the SNESSetFromOptions() command and then
3718: set the SNES solver type from the options database rather than by using
3719: this routine. Using the options database provides the user with
3720: maximum flexibility in evaluating the many nonlinear solvers.
3721: The SNESSetType() routine is provided for those situations where it
3722: is necessary to set the nonlinear solver independently of the command
3723: line or options database. This might be the case, for example, when
3724: the choice of solver changes during the execution of the program,
3725: and the user's application is taking responsibility for choosing the
3726: appropriate method.
3728: Developer Notes: SNESRegister() adds a constructor for a new SNESType to SNESList, SNESSetType() locates
3729: the constructor in that list and calls it to create the spexific object.
3731: Level: intermediate
3733: .keywords: SNES, set, type
3735: .seealso: SNESType, SNESCreate(), SNESDestroy(), SNESGetType(), SNESSetFromOptions()
3737: @*/
3738: PetscErrorCode SNESSetType(SNES snes,SNESType type)
3739: {
3740: PetscErrorCode ierr,(*r)(SNES);
3741: PetscBool match;
3747: PetscObjectTypeCompare((PetscObject)snes,type,&match);
3748: if (match) return(0);
3750: PetscFunctionListFind(SNESList,type,&r);
3751: if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unable to find requested SNES type %s",type);
3752: /* Destroy the previous private SNES context */
3753: if (snes->ops->destroy) {
3754: (*(snes)->ops->destroy)(snes);
3755: snes->ops->destroy = NULL;
3756: }
3757: /* Reinitialize function pointers in SNESOps structure */
3758: snes->ops->setup = 0;
3759: snes->ops->solve = 0;
3760: snes->ops->view = 0;
3761: snes->ops->setfromoptions = 0;
3762: snes->ops->destroy = 0;
3763: /* Call the SNESCreate_XXX routine for this particular Nonlinear solver */
3764: snes->setupcalled = PETSC_FALSE;
3766: PetscObjectChangeTypeName((PetscObject)snes,type);
3767: (*r)(snes);
3768: return(0);
3769: }
3773: /*@C
3774: SNESGetType - Gets the SNES method type and name (as a string).
3776: Not Collective
3778: Input Parameter:
3779: . snes - nonlinear solver context
3781: Output Parameter:
3782: . type - SNES method (a character string)
3784: Level: intermediate
3786: .keywords: SNES, nonlinear, get, type, name
3787: @*/
3788: PetscErrorCode SNESGetType(SNES snes,SNESType *type)
3789: {
3793: *type = ((PetscObject)snes)->type_name;
3794: return(0);
3795: }
3799: /*@
3800: SNESGetSolution - Returns the vector where the approximate solution is
3801: stored. This is the fine grid solution when using SNESSetGridSequence().
3803: Not Collective, but Vec is parallel if SNES is parallel
3805: Input Parameter:
3806: . snes - the SNES context
3808: Output Parameter:
3809: . x - the solution
3811: Level: intermediate
3813: .keywords: SNES, nonlinear, get, solution
3815: .seealso: SNESGetSolutionUpdate(), SNESGetFunction()
3816: @*/
3817: PetscErrorCode SNESGetSolution(SNES snes,Vec *x)
3818: {
3822: *x = snes->vec_sol;
3823: return(0);
3824: }
3828: /*@
3829: SNESGetSolutionUpdate - Returns the vector where the solution update is
3830: stored.
3832: Not Collective, but Vec is parallel if SNES is parallel
3834: Input Parameter:
3835: . snes - the SNES context
3837: Output Parameter:
3838: . x - the solution update
3840: Level: advanced
3842: .keywords: SNES, nonlinear, get, solution, update
3844: .seealso: SNESGetSolution(), SNESGetFunction()
3845: @*/
3846: PetscErrorCode SNESGetSolutionUpdate(SNES snes,Vec *x)
3847: {
3851: *x = snes->vec_sol_update;
3852: return(0);
3853: }
3857: /*@C
3858: SNESGetFunction - Returns the vector where the function is stored.
3860: Not Collective, but Vec is parallel if SNES is parallel. Collective if Vec is requested, but has not been created yet.
3862: Input Parameter:
3863: . snes - the SNES context
3865: Output Parameter:
3866: + r - the vector that is used to store residuals (or NULL if you don't want it)
3867: . SNESFunction- the function (or NULL if you don't want it)
3868: - ctx - the function context (or NULL if you don't want it)
3870: Level: advanced
3872: .keywords: SNES, nonlinear, get, function
3874: .seealso: SNESSetFunction(), SNESGetSolution(), SNESFunction
3875: @*/
3876: PetscErrorCode SNESGetFunction(SNES snes,Vec *r,PetscErrorCode (**SNESFunction)(SNES,Vec,Vec,void*),void **ctx)
3877: {
3879: DM dm;
3883: if (r) {
3884: if (!snes->vec_func) {
3885: if (snes->vec_rhs) {
3886: VecDuplicate(snes->vec_rhs,&snes->vec_func);
3887: } else if (snes->vec_sol) {
3888: VecDuplicate(snes->vec_sol,&snes->vec_func);
3889: } else if (snes->dm) {
3890: DMCreateGlobalVector(snes->dm,&snes->vec_func);
3891: }
3892: }
3893: *r = snes->vec_func;
3894: }
3895: SNESGetDM(snes,&dm);
3896: DMSNESGetFunction(dm,SNESFunction,ctx);
3897: return(0);
3898: }
3900: /*@C
3901: SNESGetGS - Returns the GS function and context.
3903: Input Parameter:
3904: . snes - the SNES context
3906: Output Parameter:
3907: + SNESGSFunction - the function (or NULL)
3908: - ctx - the function context (or NULL)
3910: Level: advanced
3912: .keywords: SNES, nonlinear, get, function
3914: .seealso: SNESSetGS(), SNESGetFunction()
3915: @*/
3919: PetscErrorCode SNESGetGS (SNES snes, PetscErrorCode (**SNESGSFunction)(SNES, Vec, Vec, void*), void ** ctx)
3920: {
3922: DM dm;
3926: SNESGetDM(snes,&dm);
3927: DMSNESGetGS(dm,SNESGSFunction,ctx);
3928: return(0);
3929: }
3933: /*@C
3934: SNESSetOptionsPrefix - Sets the prefix used for searching for all
3935: SNES options in the database.
3937: Logically Collective on SNES
3939: Input Parameter:
3940: + snes - the SNES context
3941: - prefix - the prefix to prepend to all option names
3943: Notes:
3944: A hyphen (-) must NOT be given at the beginning of the prefix name.
3945: The first character of all runtime options is AUTOMATICALLY the hyphen.
3947: Level: advanced
3949: .keywords: SNES, set, options, prefix, database
3951: .seealso: SNESSetFromOptions()
3952: @*/
3953: PetscErrorCode SNESSetOptionsPrefix(SNES snes,const char prefix[])
3954: {
3959: PetscObjectSetOptionsPrefix((PetscObject)snes,prefix);
3960: if (!snes->ksp) {SNESGetKSP(snes,&snes->ksp);}
3961: if (snes->linesearch) {
3962: SNESGetLineSearch(snes,&snes->linesearch);
3963: PetscObjectSetOptionsPrefix((PetscObject)snes->linesearch,prefix);
3964: }
3965: KSPSetOptionsPrefix(snes->ksp,prefix);
3966: return(0);
3967: }
3971: /*@C
3972: SNESAppendOptionsPrefix - Appends to the prefix used for searching for all
3973: SNES options in the database.
3975: Logically Collective on SNES
3977: Input Parameters:
3978: + snes - the SNES context
3979: - prefix - the prefix to prepend to all option names
3981: Notes:
3982: A hyphen (-) must NOT be given at the beginning of the prefix name.
3983: The first character of all runtime options is AUTOMATICALLY the hyphen.
3985: Level: advanced
3987: .keywords: SNES, append, options, prefix, database
3989: .seealso: SNESGetOptionsPrefix()
3990: @*/
3991: PetscErrorCode SNESAppendOptionsPrefix(SNES snes,const char prefix[])
3992: {
3997: PetscObjectAppendOptionsPrefix((PetscObject)snes,prefix);
3998: if (!snes->ksp) {SNESGetKSP(snes,&snes->ksp);}
3999: if (snes->linesearch) {
4000: SNESGetLineSearch(snes,&snes->linesearch);
4001: PetscObjectAppendOptionsPrefix((PetscObject)snes->linesearch,prefix);
4002: }
4003: KSPAppendOptionsPrefix(snes->ksp,prefix);
4004: return(0);
4005: }
4009: /*@C
4010: SNESGetOptionsPrefix - Sets the prefix used for searching for all
4011: SNES options in the database.
4013: Not Collective
4015: Input Parameter:
4016: . snes - the SNES context
4018: Output Parameter:
4019: . prefix - pointer to the prefix string used
4021: Notes: On the fortran side, the user should pass in a string 'prefix' of
4022: sufficient length to hold the prefix.
4024: Level: advanced
4026: .keywords: SNES, get, options, prefix, database
4028: .seealso: SNESAppendOptionsPrefix()
4029: @*/
4030: PetscErrorCode SNESGetOptionsPrefix(SNES snes,const char *prefix[])
4031: {
4036: PetscObjectGetOptionsPrefix((PetscObject)snes,prefix);
4037: return(0);
4038: }
4043: /*@C
4044: SNESRegister - Adds a method to the nonlinear solver package.
4046: Not collective
4048: Input Parameters:
4049: + name_solver - name of a new user-defined solver
4050: - routine_create - routine to create method context
4052: Notes:
4053: SNESRegister() may be called multiple times to add several user-defined solvers.
4055: Sample usage:
4056: .vb
4057: SNESRegister("my_solver",MySolverCreate);
4058: .ve
4060: Then, your solver can be chosen with the procedural interface via
4061: $ SNESSetType(snes,"my_solver")
4062: or at runtime via the option
4063: $ -snes_type my_solver
4065: Level: advanced
4067: Note: If your function is not being put into a shared library then use SNESRegister() instead
4069: .keywords: SNES, nonlinear, register
4071: .seealso: SNESRegisterAll(), SNESRegisterDestroy()
4073: Level: advanced
4074: @*/
4075: PetscErrorCode SNESRegister(const char sname[],PetscErrorCode (*function)(SNES))
4076: {
4080: PetscFunctionListAdd(&SNESList,sname,function);
4081: return(0);
4082: }
4086: PetscErrorCode SNESTestLocalMin(SNES snes)
4087: {
4089: PetscInt N,i,j;
4090: Vec u,uh,fh;
4091: PetscScalar value;
4092: PetscReal norm;
4095: SNESGetSolution(snes,&u);
4096: VecDuplicate(u,&uh);
4097: VecDuplicate(u,&fh);
4099: /* currently only works for sequential */
4100: PetscPrintf(PETSC_COMM_WORLD,"Testing FormFunction() for local min\n");
4101: VecGetSize(u,&N);
4102: for (i=0; i<N; i++) {
4103: VecCopy(u,uh);
4104: PetscPrintf(PETSC_COMM_WORLD,"i = %D\n",i);
4105: for (j=-10; j<11; j++) {
4106: value = PetscSign(j)*exp(PetscAbs(j)-10.0);
4107: VecSetValue(uh,i,value,ADD_VALUES);
4108: SNESComputeFunction(snes,uh,fh);
4109: VecNorm(fh,NORM_2,&norm);
4110: PetscPrintf(PETSC_COMM_WORLD," j norm %D %18.16e\n",j,norm);
4111: value = -value;
4112: VecSetValue(uh,i,value,ADD_VALUES);
4113: }
4114: }
4115: VecDestroy(&uh);
4116: VecDestroy(&fh);
4117: return(0);
4118: }
4122: /*@
4123: SNESKSPSetUseEW - Sets SNES use Eisenstat-Walker method for
4124: computing relative tolerance for linear solvers within an inexact
4125: Newton method.
4127: Logically Collective on SNES
4129: Input Parameters:
4130: + snes - SNES context
4131: - flag - PETSC_TRUE or PETSC_FALSE
4133: Options Database:
4134: + -snes_ksp_ew - use Eisenstat-Walker method for determining linear system convergence
4135: . -snes_ksp_ew_version ver - version of Eisenstat-Walker method
4136: . -snes_ksp_ew_rtol0 <rtol0> - Sets rtol0
4137: . -snes_ksp_ew_rtolmax <rtolmax> - Sets rtolmax
4138: . -snes_ksp_ew_gamma <gamma> - Sets gamma
4139: . -snes_ksp_ew_alpha <alpha> - Sets alpha
4140: . -snes_ksp_ew_alpha2 <alpha2> - Sets alpha2
4141: - -snes_ksp_ew_threshold <threshold> - Sets threshold
4143: Notes:
4144: Currently, the default is to use a constant relative tolerance for
4145: the inner linear solvers. Alternatively, one can use the
4146: Eisenstat-Walker method, where the relative convergence tolerance
4147: is reset at each Newton iteration according progress of the nonlinear
4148: solver.
4150: Level: advanced
4152: Reference:
4153: S. C. Eisenstat and H. F. Walker, "Choosing the forcing terms in an
4154: inexact Newton method", SISC 17 (1), pp.16-32, 1996.
4156: .keywords: SNES, KSP, Eisenstat, Walker, convergence, test, inexact, Newton
4158: .seealso: SNESKSPGetUseEW(), SNESKSPGetParametersEW(), SNESKSPSetParametersEW()
4159: @*/
4160: PetscErrorCode SNESKSPSetUseEW(SNES snes,PetscBool flag)
4161: {
4165: snes->ksp_ewconv = flag;
4166: return(0);
4167: }
4171: /*@
4172: SNESKSPGetUseEW - Gets if SNES is using Eisenstat-Walker method
4173: for computing relative tolerance for linear solvers within an
4174: inexact Newton method.
4176: Not Collective
4178: Input Parameter:
4179: . snes - SNES context
4181: Output Parameter:
4182: . flag - PETSC_TRUE or PETSC_FALSE
4184: Notes:
4185: Currently, the default is to use a constant relative tolerance for
4186: the inner linear solvers. Alternatively, one can use the
4187: Eisenstat-Walker method, where the relative convergence tolerance
4188: is reset at each Newton iteration according progress of the nonlinear
4189: solver.
4191: Level: advanced
4193: Reference:
4194: S. C. Eisenstat and H. F. Walker, "Choosing the forcing terms in an
4195: inexact Newton method", SISC 17 (1), pp.16-32, 1996.
4197: .keywords: SNES, KSP, Eisenstat, Walker, convergence, test, inexact, Newton
4199: .seealso: SNESKSPSetUseEW(), SNESKSPGetParametersEW(), SNESKSPSetParametersEW()
4200: @*/
4201: PetscErrorCode SNESKSPGetUseEW(SNES snes, PetscBool *flag)
4202: {
4206: *flag = snes->ksp_ewconv;
4207: return(0);
4208: }
4212: /*@
4213: SNESKSPSetParametersEW - Sets parameters for Eisenstat-Walker
4214: convergence criteria for the linear solvers within an inexact
4215: Newton method.
4217: Logically Collective on SNES
4219: Input Parameters:
4220: + snes - SNES context
4221: . version - version 1, 2 (default is 2) or 3
4222: . rtol_0 - initial relative tolerance (0 <= rtol_0 < 1)
4223: . rtol_max - maximum relative tolerance (0 <= rtol_max < 1)
4224: . gamma - multiplicative factor for version 2 rtol computation
4225: (0 <= gamma2 <= 1)
4226: . alpha - power for version 2 rtol computation (1 < alpha <= 2)
4227: . alpha2 - power for safeguard
4228: - threshold - threshold for imposing safeguard (0 < threshold < 1)
4230: Note:
4231: Version 3 was contributed by Luis Chacon, June 2006.
4233: Use PETSC_DEFAULT to retain the default for any of the parameters.
4235: Level: advanced
4237: Reference:
4238: S. C. Eisenstat and H. F. Walker, "Choosing the forcing terms in an
4239: inexact Newton method", Utah State University Math. Stat. Dept. Res.
4240: Report 6/94/75, June, 1994, to appear in SIAM J. Sci. Comput.
4242: .keywords: SNES, KSP, Eisenstat, Walker, set, parameters
4244: .seealso: SNESKSPSetUseEW(), SNESKSPGetUseEW(), SNESKSPGetParametersEW()
4245: @*/
4246: PetscErrorCode SNESKSPSetParametersEW(SNES snes,PetscInt version,PetscReal rtol_0,PetscReal rtol_max,PetscReal gamma,PetscReal alpha,PetscReal alpha2,PetscReal threshold)
4247: {
4248: SNESKSPEW *kctx;
4252: kctx = (SNESKSPEW*)snes->kspconvctx;
4253: if (!kctx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"No Eisenstat-Walker context existing");
4262: if (version != PETSC_DEFAULT) kctx->version = version;
4263: if (rtol_0 != PETSC_DEFAULT) kctx->rtol_0 = rtol_0;
4264: if (rtol_max != PETSC_DEFAULT) kctx->rtol_max = rtol_max;
4265: if (gamma != PETSC_DEFAULT) kctx->gamma = gamma;
4266: if (alpha != PETSC_DEFAULT) kctx->alpha = alpha;
4267: if (alpha2 != PETSC_DEFAULT) kctx->alpha2 = alpha2;
4268: if (threshold != PETSC_DEFAULT) kctx->threshold = threshold;
4270: 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);
4271: 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",kctx->rtol_0);
4272: 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",kctx->rtol_max);
4273: if (kctx->gamma < 0.0 || kctx->gamma > 1.0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"0.0 <= gamma (%G) <= 1.0\n",kctx->gamma);
4274: if (kctx->alpha <= 1.0 || kctx->alpha > 2.0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"1.0 < alpha (%G) <= 2.0\n",kctx->alpha);
4275: if (kctx->threshold <= 0.0 || kctx->threshold >= 1.0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"0.0 < threshold (%G) < 1.0\n",kctx->threshold);
4276: return(0);
4277: }
4281: /*@
4282: SNESKSPGetParametersEW - Gets parameters for Eisenstat-Walker
4283: convergence criteria for the linear solvers within an inexact
4284: Newton method.
4286: Not Collective
4288: Input Parameters:
4289: snes - SNES context
4291: Output Parameters:
4292: + version - version 1, 2 (default is 2) or 3
4293: . rtol_0 - initial relative tolerance (0 <= rtol_0 < 1)
4294: . rtol_max - maximum relative tolerance (0 <= rtol_max < 1)
4295: . gamma - multiplicative factor for version 2 rtol computation (0 <= gamma2 <= 1)
4296: . alpha - power for version 2 rtol computation (1 < alpha <= 2)
4297: . alpha2 - power for safeguard
4298: - threshold - threshold for imposing safeguard (0 < threshold < 1)
4300: Level: advanced
4302: .keywords: SNES, KSP, Eisenstat, Walker, get, parameters
4304: .seealso: SNESKSPSetUseEW(), SNESKSPGetUseEW(), SNESKSPSetParametersEW()
4305: @*/
4306: PetscErrorCode SNESKSPGetParametersEW(SNES snes,PetscInt *version,PetscReal *rtol_0,PetscReal *rtol_max,PetscReal *gamma,PetscReal *alpha,PetscReal *alpha2,PetscReal *threshold)
4307: {
4308: SNESKSPEW *kctx;
4312: kctx = (SNESKSPEW*)snes->kspconvctx;
4313: if (!kctx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"No Eisenstat-Walker context existing");
4314: if (version) *version = kctx->version;
4315: if (rtol_0) *rtol_0 = kctx->rtol_0;
4316: if (rtol_max) *rtol_max = kctx->rtol_max;
4317: if (gamma) *gamma = kctx->gamma;
4318: if (alpha) *alpha = kctx->alpha;
4319: if (alpha2) *alpha2 = kctx->alpha2;
4320: if (threshold) *threshold = kctx->threshold;
4321: return(0);
4322: }
4326: PetscErrorCode SNESKSPEW_PreSolve(KSP ksp, Vec b, Vec x, SNES snes)
4327: {
4329: SNESKSPEW *kctx = (SNESKSPEW*)snes->kspconvctx;
4330: PetscReal rtol = PETSC_DEFAULT,stol;
4333: if (!snes->ksp_ewconv) return(0);
4334: if (!snes->iter) rtol = kctx->rtol_0; /* first time in, so use the original user rtol */
4335: else {
4336: if (kctx->version == 1) {
4337: rtol = (snes->norm - kctx->lresid_last)/kctx->norm_last;
4338: if (rtol < 0.0) rtol = -rtol;
4339: stol = PetscPowReal(kctx->rtol_last,kctx->alpha2);
4340: if (stol > kctx->threshold) rtol = PetscMax(rtol,stol);
4341: } else if (kctx->version == 2) {
4342: rtol = kctx->gamma * PetscPowReal(snes->norm/kctx->norm_last,kctx->alpha);
4343: stol = kctx->gamma * PetscPowReal(kctx->rtol_last,kctx->alpha);
4344: if (stol > kctx->threshold) rtol = PetscMax(rtol,stol);
4345: } else if (kctx->version == 3) { /* contributed by Luis Chacon, June 2006. */
4346: rtol = kctx->gamma * PetscPowReal(snes->norm/kctx->norm_last,kctx->alpha);
4347: /* safeguard: avoid sharp decrease of rtol */
4348: stol = kctx->gamma*PetscPowReal(kctx->rtol_last,kctx->alpha);
4349: stol = PetscMax(rtol,stol);
4350: rtol = PetscMin(kctx->rtol_0,stol);
4351: /* safeguard: avoid oversolving */
4352: stol = kctx->gamma*(snes->ttol)/snes->norm;
4353: stol = PetscMax(rtol,stol);
4354: rtol = PetscMin(kctx->rtol_0,stol);
4355: } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Only versions 1, 2 or 3 are supported: %D",kctx->version);
4356: }
4357: /* safeguard: avoid rtol greater than one */
4358: rtol = PetscMin(rtol,kctx->rtol_max);
4359: KSPSetTolerances(ksp,rtol,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
4360: PetscInfo3(snes,"iter %D, Eisenstat-Walker (version %D) KSP rtol=%G\n",snes->iter,kctx->version,rtol);
4361: return(0);
4362: }
4366: PetscErrorCode SNESKSPEW_PostSolve(KSP ksp, Vec b, Vec x, SNES snes)
4367: {
4369: SNESKSPEW *kctx = (SNESKSPEW*)snes->kspconvctx;
4370: PCSide pcside;
4371: Vec lres;
4374: if (!snes->ksp_ewconv) return(0);
4375: KSPGetTolerances(ksp,&kctx->rtol_last,0,0,0);
4376: SNESGetFunctionNorm(snes,&kctx->norm_last);
4377: if (kctx->version == 1) {
4378: KSPGetPCSide(ksp,&pcside);
4379: if (pcside == PC_RIGHT) { /* XXX Should we also test KSP_UNPRECONDITIONED_NORM ? */
4380: /* KSP residual is true linear residual */
4381: KSPGetResidualNorm(ksp,&kctx->lresid_last);
4382: } else {
4383: /* KSP residual is preconditioned residual */
4384: /* compute true linear residual norm */
4385: VecDuplicate(b,&lres);
4386: MatMult(snes->jacobian,x,lres);
4387: VecAYPX(lres,-1.0,b);
4388: VecNorm(lres,NORM_2,&kctx->lresid_last);
4389: VecDestroy(&lres);
4390: }
4391: }
4392: return(0);
4393: }
4397: /*@
4398: SNESGetKSP - Returns the KSP context for a SNES solver.
4400: Not Collective, but if SNES object is parallel, then KSP object is parallel
4402: Input Parameter:
4403: . snes - the SNES context
4405: Output Parameter:
4406: . ksp - the KSP context
4408: Notes:
4409: The user can then directly manipulate the KSP context to set various
4410: options, etc. Likewise, the user can then extract and manipulate the
4411: PC contexts as well.
4413: Level: beginner
4415: .keywords: SNES, nonlinear, get, KSP, context
4417: .seealso: KSPGetPC(), SNESCreate(), KSPCreate(), SNESSetKSP()
4418: @*/
4419: PetscErrorCode SNESGetKSP(SNES snes,KSP *ksp)
4420: {
4427: if (!snes->ksp) {
4428: KSPCreate(PetscObjectComm((PetscObject)snes),&snes->ksp);
4429: PetscObjectIncrementTabLevel((PetscObject)snes->ksp,(PetscObject)snes,1);
4430: PetscLogObjectParent(snes,snes->ksp);
4432: KSPSetPreSolve(snes->ksp,(PetscErrorCode (*)(KSP,Vec,Vec,void*))SNESKSPEW_PreSolve,snes);
4433: KSPSetPostSolve(snes->ksp,(PetscErrorCode (*)(KSP,Vec,Vec,void*))SNESKSPEW_PostSolve,snes);
4434: }
4435: *ksp = snes->ksp;
4436: return(0);
4437: }
4440: #include <petsc-private/dmimpl.h>
4443: /*@
4444: SNESSetDM - Sets the DM that may be used by some preconditioners
4446: Logically Collective on SNES
4448: Input Parameters:
4449: + snes - the preconditioner context
4450: - dm - the dm
4452: Level: intermediate
4454: .seealso: SNESGetDM(), KSPSetDM(), KSPGetDM()
4455: @*/
4456: PetscErrorCode SNESSetDM(SNES snes,DM dm)
4457: {
4459: KSP ksp;
4460: DMSNES sdm;
4464: if (dm) {PetscObjectReference((PetscObject)dm);}
4465: if (snes->dm) { /* Move the DMSNES context over to the new DM unless the new DM already has one */
4466: if (snes->dm->dmsnes && snes->dmAuto && !dm->dmsnes) {
4467: DMCopyDMSNES(snes->dm,dm);
4468: DMGetDMSNES(snes->dm,&sdm);
4469: if (sdm->originaldm == snes->dm) sdm->originaldm = dm; /* Grant write privileges to the replacement DM */
4470: }
4471: DMDestroy(&snes->dm);
4472: }
4473: snes->dm = dm;
4474: snes->dmAuto = PETSC_FALSE;
4476: SNESGetKSP(snes,&ksp);
4477: KSPSetDM(ksp,dm);
4478: KSPSetDMActive(ksp,PETSC_FALSE);
4479: if (snes->pc) {
4480: SNESSetDM(snes->pc, snes->dm);
4481: SNESSetPCSide(snes,snes->pcside);
4482: }
4483: return(0);
4484: }
4488: /*@
4489: SNESGetDM - Gets the DM that may be used by some preconditioners
4491: Not Collective but DM obtained is parallel on SNES
4493: Input Parameter:
4494: . snes - the preconditioner context
4496: Output Parameter:
4497: . dm - the dm
4499: Level: intermediate
4501: .seealso: SNESSetDM(), KSPSetDM(), KSPGetDM()
4502: @*/
4503: PetscErrorCode SNESGetDM(SNES snes,DM *dm)
4504: {
4509: if (!snes->dm) {
4510: DMShellCreate(PetscObjectComm((PetscObject)snes),&snes->dm);
4511: snes->dmAuto = PETSC_TRUE;
4512: }
4513: *dm = snes->dm;
4514: return(0);
4515: }
4519: /*@
4520: SNESSetPC - Sets the nonlinear preconditioner to be used.
4522: Collective on SNES
4524: Input Parameters:
4525: + snes - iterative context obtained from SNESCreate()
4526: - pc - the preconditioner object
4528: Notes:
4529: Use SNESGetPC() to retrieve the preconditioner context (for example,
4530: to configure it using the API).
4532: Level: developer
4534: .keywords: SNES, set, precondition
4535: .seealso: SNESGetPC()
4536: @*/
4537: PetscErrorCode SNESSetPC(SNES snes, SNES pc)
4538: {
4545: PetscObjectReference((PetscObject) pc);
4546: SNESDestroy(&snes->pc);
4547: snes->pc = pc;
4548: PetscLogObjectParent(snes, snes->pc);
4549: return(0);
4550: }
4554: /*@
4555: SNESGetPC - Returns a pointer to the nonlinear preconditioning context set with SNESSetPC().
4557: Not Collective
4559: Input Parameter:
4560: . snes - iterative context obtained from SNESCreate()
4562: Output Parameter:
4563: . pc - preconditioner context
4565: Level: developer
4567: .keywords: SNES, get, preconditioner
4568: .seealso: SNESSetPC()
4569: @*/
4570: PetscErrorCode SNESGetPC(SNES snes, SNES *pc)
4571: {
4573: const char *optionsprefix;
4578: if (!snes->pc) {
4579: SNESCreate(PetscObjectComm((PetscObject)snes),&snes->pc);
4580: PetscObjectIncrementTabLevel((PetscObject)snes->pc,(PetscObject)snes,1);
4581: PetscLogObjectParent(snes,snes->pc);
4582: SNESGetOptionsPrefix(snes,&optionsprefix);
4583: SNESSetOptionsPrefix(snes->pc,optionsprefix);
4584: SNESAppendOptionsPrefix(snes->pc,"npc_");
4585: }
4586: *pc = snes->pc;
4587: return(0);
4588: }
4592: /*@
4593: SNESSetPCSide - Sets the preconditioning side.
4595: Logically Collective on SNES
4597: Input Parameter:
4598: . snes - iterative context obtained from SNESCreate()
4600: Output Parameter:
4601: . side - the preconditioning side, where side is one of
4602: .vb
4603: PC_LEFT - left preconditioning (default)
4604: PC_RIGHT - right preconditioning
4605: .ve
4607: Options Database Keys:
4608: . -snes_pc_side <right,left>
4610: Level: intermediate
4612: .keywords: SNES, set, right, left, side, preconditioner, flag
4614: .seealso: SNESGetPCSide(), KSPSetPCSide()
4615: @*/
4616: PetscErrorCode SNESSetPCSide(SNES snes,PCSide side)
4617: {
4621: snes->pcside = side;
4622: return(0);
4623: }
4627: /*@
4628: SNESGetPCSide - Gets the preconditioning side.
4630: Not Collective
4632: Input Parameter:
4633: . snes - iterative context obtained from SNESCreate()
4635: Output Parameter:
4636: . side - the preconditioning side, where side is one of
4637: .vb
4638: PC_LEFT - left preconditioning (default)
4639: PC_RIGHT - right preconditioning
4640: .ve
4642: Level: intermediate
4644: .keywords: SNES, get, right, left, side, preconditioner, flag
4646: .seealso: SNESSetPCSide(), KSPGetPCSide()
4647: @*/
4648: PetscErrorCode SNESGetPCSide(SNES snes,PCSide *side)
4649: {
4653: *side = snes->pcside;
4654: return(0);
4655: }
4659: /*@
4660: SNESSetLineSearch - Sets the linesearch on the SNES instance.
4662: Collective on SNES
4664: Input Parameters:
4665: + snes - iterative context obtained from SNESCreate()
4666: - linesearch - the linesearch object
4668: Notes:
4669: Use SNESGetLineSearch() to retrieve the preconditioner context (for example,
4670: to configure it using the API).
4672: Level: developer
4674: .keywords: SNES, set, linesearch
4675: .seealso: SNESGetLineSearch()
4676: @*/
4677: PetscErrorCode SNESSetLineSearch(SNES snes, SNESLineSearch linesearch)
4678: {
4685: PetscObjectReference((PetscObject) linesearch);
4686: SNESLineSearchDestroy(&snes->linesearch);
4688: snes->linesearch = linesearch;
4690: PetscLogObjectParent(snes, snes->linesearch);
4691: return(0);
4692: }
4696: /*@
4697: SNESGetLineSearch - Returns a pointer to the line search context set with SNESSetLineSearch()
4698: or creates a default line search instance associated with the SNES and returns it.
4700: Not Collective
4702: Input Parameter:
4703: . snes - iterative context obtained from SNESCreate()
4705: Output Parameter:
4706: . linesearch - linesearch context
4708: Level: developer
4710: .keywords: SNES, get, linesearch
4711: .seealso: SNESSetLineSearch()
4712: @*/
4713: PetscErrorCode SNESGetLineSearch(SNES snes, SNESLineSearch *linesearch)
4714: {
4716: const char *optionsprefix;
4721: if (!snes->linesearch) {
4722: SNESGetOptionsPrefix(snes, &optionsprefix);
4723: SNESLineSearchCreate(PetscObjectComm((PetscObject)snes), &snes->linesearch);
4724: SNESLineSearchSetSNES(snes->linesearch, snes);
4725: SNESLineSearchAppendOptionsPrefix(snes->linesearch, optionsprefix);
4726: PetscObjectIncrementTabLevel((PetscObject) snes->linesearch, (PetscObject) snes, 1);
4727: PetscLogObjectParent(snes, snes->linesearch);
4728: }
4729: *linesearch = snes->linesearch;
4730: return(0);
4731: }
4733: #if defined(PETSC_HAVE_MATLAB_ENGINE)
4734: #include <mex.h>
4736: typedef struct {char *funcname; mxArray *ctx;} SNESMatlabContext;
4740: /*
4741: SNESComputeFunction_Matlab - Calls the function that has been set with SNESSetFunctionMatlab().
4743: Collective on SNES
4745: Input Parameters:
4746: + snes - the SNES context
4747: - x - input vector
4749: Output Parameter:
4750: . y - function vector, as set by SNESSetFunction()
4752: Notes:
4753: SNESComputeFunction() is typically used within nonlinear solvers
4754: implementations, so most users would not generally call this routine
4755: themselves.
4757: Level: developer
4759: .keywords: SNES, nonlinear, compute, function
4761: .seealso: SNESSetFunction(), SNESGetFunction()
4762: */
4763: PetscErrorCode SNESComputeFunction_Matlab(SNES snes,Vec x,Vec y, void *ctx)
4764: {
4765: PetscErrorCode ierr;
4766: SNESMatlabContext *sctx = (SNESMatlabContext*)ctx;
4767: int nlhs = 1,nrhs = 5;
4768: mxArray *plhs[1],*prhs[5];
4769: long long int lx = 0,ly = 0,ls = 0;
4778: /* call Matlab function in ctx with arguments x and y */
4780: PetscMemcpy(&ls,&snes,sizeof(snes));
4781: PetscMemcpy(&lx,&x,sizeof(x));
4782: PetscMemcpy(&ly,&y,sizeof(x));
4783: prhs[0] = mxCreateDoubleScalar((double)ls);
4784: prhs[1] = mxCreateDoubleScalar((double)lx);
4785: prhs[2] = mxCreateDoubleScalar((double)ly);
4786: prhs[3] = mxCreateString(sctx->funcname);
4787: prhs[4] = sctx->ctx;
4788: mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscSNESComputeFunctionInternal");
4789: mxGetScalar(plhs[0]);
4790: mxDestroyArray(prhs[0]);
4791: mxDestroyArray(prhs[1]);
4792: mxDestroyArray(prhs[2]);
4793: mxDestroyArray(prhs[3]);
4794: mxDestroyArray(plhs[0]);
4795: return(0);
4796: }
4800: /*
4801: SNESSetFunctionMatlab - Sets the function evaluation routine and function
4802: vector for use by the SNES routines in solving systems of nonlinear
4803: equations from MATLAB. Here the function is a string containing the name of a MATLAB function
4805: Logically Collective on SNES
4807: Input Parameters:
4808: + snes - the SNES context
4809: . r - vector to store function value
4810: - SNESFunction - function evaluation routine
4812: Notes:
4813: The Newton-like methods typically solve linear systems of the form
4814: $ f'(x) x = -f(x),
4815: where f'(x) denotes the Jacobian matrix and f(x) is the function.
4817: Level: beginner
4819: Developer Note: This bleeds the allocated memory SNESMatlabContext *sctx;
4821: .keywords: SNES, nonlinear, set, function
4823: .seealso: SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESSetFunction()
4824: */
4825: PetscErrorCode SNESSetFunctionMatlab(SNES snes,Vec r,const char *SNESFunction,mxArray *ctx)
4826: {
4827: PetscErrorCode ierr;
4828: SNESMatlabContext *sctx;
4831: /* currently sctx is memory bleed */
4832: PetscMalloc(sizeof(SNESMatlabContext),&sctx);
4833: PetscStrallocpy(SNESFunction,&sctx->funcname);
4834: /*
4835: This should work, but it doesn't
4836: sctx->ctx = ctx;
4837: mexMakeArrayPersistent(sctx->ctx);
4838: */
4839: sctx->ctx = mxDuplicateArray(ctx);
4840: SNESSetFunction(snes,r,SNESComputeFunction_Matlab,sctx);
4841: return(0);
4842: }
4846: /*
4847: SNESComputeJacobian_Matlab - Calls the function that has been set with SNESSetJacobianMatlab().
4849: Collective on SNES
4851: Input Parameters:
4852: + snes - the SNES context
4853: . x - input vector
4854: . A, B - the matrices
4855: - ctx - user context
4857: Output Parameter:
4858: . flag - structure of the matrix
4860: Level: developer
4862: .keywords: SNES, nonlinear, compute, function
4864: .seealso: SNESSetFunction(), SNESGetFunction()
4865: @*/
4866: PetscErrorCode SNESComputeJacobian_Matlab(SNES snes,Vec x,Mat *A,Mat *B,MatStructure *flag, void *ctx)
4867: {
4868: PetscErrorCode ierr;
4869: SNESMatlabContext *sctx = (SNESMatlabContext*)ctx;
4870: int nlhs = 2,nrhs = 6;
4871: mxArray *plhs[2],*prhs[6];
4872: long long int lx = 0,lA = 0,ls = 0, lB = 0;
4878: /* call Matlab function in ctx with arguments x and y */
4880: PetscMemcpy(&ls,&snes,sizeof(snes));
4881: PetscMemcpy(&lx,&x,sizeof(x));
4882: PetscMemcpy(&lA,A,sizeof(x));
4883: PetscMemcpy(&lB,B,sizeof(x));
4884: prhs[0] = mxCreateDoubleScalar((double)ls);
4885: prhs[1] = mxCreateDoubleScalar((double)lx);
4886: prhs[2] = mxCreateDoubleScalar((double)lA);
4887: prhs[3] = mxCreateDoubleScalar((double)lB);
4888: prhs[4] = mxCreateString(sctx->funcname);
4889: prhs[5] = sctx->ctx;
4890: mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscSNESComputeJacobianInternal");
4891: mxGetScalar(plhs[0]);
4892: *flag = (MatStructure) mxGetScalar(plhs[1]);
4893: mxDestroyArray(prhs[0]);
4894: mxDestroyArray(prhs[1]);
4895: mxDestroyArray(prhs[2]);
4896: mxDestroyArray(prhs[3]);
4897: mxDestroyArray(prhs[4]);
4898: mxDestroyArray(plhs[0]);
4899: mxDestroyArray(plhs[1]);
4900: return(0);
4901: }
4905: /*
4906: SNESSetJacobianMatlab - Sets the Jacobian function evaluation routine and two empty Jacobian matrices
4907: vector for use by the SNES routines in solving systems of nonlinear
4908: equations from MATLAB. Here the function is a string containing the name of a MATLAB function
4910: Logically Collective on SNES
4912: Input Parameters:
4913: + snes - the SNES context
4914: . A,B - Jacobian matrices
4915: . SNESJacobianFunction - function evaluation routine
4916: - ctx - user context
4918: Level: developer
4920: Developer Note: This bleeds the allocated memory SNESMatlabContext *sctx;
4922: .keywords: SNES, nonlinear, set, function
4924: .seealso: SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESSetFunction(), SNESJacobianFunction
4925: */
4926: PetscErrorCode SNESSetJacobianMatlab(SNES snes,Mat A,Mat B,const char *SNESJacobianFunction,mxArray *ctx)
4927: {
4928: PetscErrorCode ierr;
4929: SNESMatlabContext *sctx;
4932: /* currently sctx is memory bleed */
4933: PetscMalloc(sizeof(SNESMatlabContext),&sctx);
4934: PetscStrallocpy(SNESJacobianFunction,&sctx->funcname);
4935: /*
4936: This should work, but it doesn't
4937: sctx->ctx = ctx;
4938: mexMakeArrayPersistent(sctx->ctx);
4939: */
4940: sctx->ctx = mxDuplicateArray(ctx);
4941: SNESSetJacobian(snes,A,B,SNESComputeJacobian_Matlab,sctx);
4942: return(0);
4943: }
4947: /*
4948: SNESMonitor_Matlab - Calls the function that has been set with SNESMonitorSetMatlab().
4950: Collective on SNES
4952: .seealso: SNESSetFunction(), SNESGetFunction()
4953: @*/
4954: PetscErrorCode SNESMonitor_Matlab(SNES snes,PetscInt it, PetscReal fnorm, void *ctx)
4955: {
4956: PetscErrorCode ierr;
4957: SNESMatlabContext *sctx = (SNESMatlabContext*)ctx;
4958: int nlhs = 1,nrhs = 6;
4959: mxArray *plhs[1],*prhs[6];
4960: long long int lx = 0,ls = 0;
4961: Vec x = snes->vec_sol;
4966: PetscMemcpy(&ls,&snes,sizeof(snes));
4967: PetscMemcpy(&lx,&x,sizeof(x));
4968: prhs[0] = mxCreateDoubleScalar((double)ls);
4969: prhs[1] = mxCreateDoubleScalar((double)it);
4970: prhs[2] = mxCreateDoubleScalar((double)fnorm);
4971: prhs[3] = mxCreateDoubleScalar((double)lx);
4972: prhs[4] = mxCreateString(sctx->funcname);
4973: prhs[5] = sctx->ctx;
4974: mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscSNESMonitorInternal");
4975: mxGetScalar(plhs[0]);
4976: mxDestroyArray(prhs[0]);
4977: mxDestroyArray(prhs[1]);
4978: mxDestroyArray(prhs[2]);
4979: mxDestroyArray(prhs[3]);
4980: mxDestroyArray(prhs[4]);
4981: mxDestroyArray(plhs[0]);
4982: return(0);
4983: }
4987: /*
4988: SNESMonitorSetMatlab - Sets the monitor function from MATLAB
4990: Level: developer
4992: Developer Note: This bleeds the allocated memory SNESMatlabContext *sctx;
4994: .keywords: SNES, nonlinear, set, function
4996: .seealso: SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESSetFunction()
4997: */
4998: PetscErrorCode SNESMonitorSetMatlab(SNES snes,const char *SNESMonitorFunction,mxArray *ctx)
4999: {
5000: PetscErrorCode ierr;
5001: SNESMatlabContext *sctx;
5004: /* currently sctx is memory bleed */
5005: PetscMalloc(sizeof(SNESMatlabContext),&sctx);
5006: PetscStrallocpy(SNESMonitorFunction,&sctx->funcname);
5007: /*
5008: This should work, but it doesn't
5009: sctx->ctx = ctx;
5010: mexMakeArrayPersistent(sctx->ctx);
5011: */
5012: sctx->ctx = mxDuplicateArray(ctx);
5013: SNESMonitorSet(snes,SNESMonitor_Matlab,sctx,NULL);
5014: return(0);
5015: }
5017: #endif