1 2 #include <petsc-private/snesimpl.h> /*I "petscsnes.h" I*/ 3 #include <petscdmshell.h> 4 5 PetscBool SNESRegisterAllCalled = PETSC_FALSE; 6 PetscFunctionList SNESList = NULL; 7 8 /* Logging support */ 9 PetscClassId SNES_CLASSID, DMSNES_CLASSID; 10 PetscLogEvent SNES_Solve, SNES_FunctionEval, SNES_JacobianEval, SNES_GSEval, SNES_NPCSolve; 11 12 #undef __FUNCT__ 13 #define __FUNCT__ "SNESSetErrorIfNotConverged" 14 /*@ 15 SNESSetErrorIfNotConverged - Causes SNESSolve() to generate an error if the solver has not converged. 16 17 Logically Collective on SNES 18 19 Input Parameters: 20 + snes - iterative context obtained from SNESCreate() 21 - flg - PETSC_TRUE indicates you want the error generated 22 23 Options database keys: 24 . -snes_error_if_not_converged : this takes an optional truth value (0/1/no/yes/true/false) 25 26 Level: intermediate 27 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. 31 32 .keywords: SNES, set, initial guess, nonzero 33 34 .seealso: SNESGetErrorIfNotConverged(), KSPGetErrorIfNotConverged(), KSPSetErrorIFNotConverged() 35 @*/ 36 PetscErrorCode SNESSetErrorIfNotConverged(SNES snes,PetscBool flg) 37 { 38 PetscFunctionBegin; 39 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 40 PetscValidLogicalCollectiveBool(snes,flg,2); 41 snes->errorifnotconverged = flg; 42 PetscFunctionReturn(0); 43 } 44 45 #undef __FUNCT__ 46 #define __FUNCT__ "SNESGetErrorIfNotConverged" 47 /*@ 48 SNESGetErrorIfNotConverged - Will SNESSolve() generate an error if the solver does not converge? 49 50 Not Collective 51 52 Input Parameter: 53 . snes - iterative context obtained from SNESCreate() 54 55 Output Parameter: 56 . flag - PETSC_TRUE if it will generate an error, else PETSC_FALSE 57 58 Level: intermediate 59 60 .keywords: SNES, set, initial guess, nonzero 61 62 .seealso: SNESSetErrorIfNotConverged(), KSPGetErrorIfNotConverged(), KSPSetErrorIFNotConverged() 63 @*/ 64 PetscErrorCode SNESGetErrorIfNotConverged(SNES snes,PetscBool *flag) 65 { 66 PetscFunctionBegin; 67 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 68 PetscValidPointer(flag,2); 69 *flag = snes->errorifnotconverged; 70 PetscFunctionReturn(0); 71 } 72 73 #undef __FUNCT__ 74 #define __FUNCT__ "SNESSetFunctionDomainError" 75 /*@ 76 SNESSetFunctionDomainError - tells SNES that the input vector to your SNESFunction is not 77 in the functions domain. For example, negative pressure. 78 79 Logically Collective on SNES 80 81 Input Parameters: 82 . snes - the SNES context 83 84 Level: advanced 85 86 .keywords: SNES, view 87 88 .seealso: SNESCreate(), SNESSetFunction(), SNESFunction 89 @*/ 90 PetscErrorCode SNESSetFunctionDomainError(SNES snes) 91 { 92 PetscFunctionBegin; 93 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 94 snes->domainerror = PETSC_TRUE; 95 PetscFunctionReturn(0); 96 } 97 98 #undef __FUNCT__ 99 #define __FUNCT__ "SNESGetFunctionDomainError" 100 /*@ 101 SNESGetFunctionDomainError - Gets the status of the domain error after a call to SNESComputeFunction; 102 103 Logically Collective on SNES 104 105 Input Parameters: 106 . snes - the SNES context 107 108 Output Parameters: 109 . domainerror - Set to PETSC_TRUE if there's a domain error; PETSC_FALSE otherwise. 110 111 Level: advanced 112 113 .keywords: SNES, view 114 115 .seealso: SNESSetFunctionDomainError(), SNESComputeFunction() 116 @*/ 117 PetscErrorCode SNESGetFunctionDomainError(SNES snes, PetscBool *domainerror) 118 { 119 PetscFunctionBegin; 120 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 121 PetscValidPointer(domainerror, 2); 122 *domainerror = snes->domainerror; 123 PetscFunctionReturn(0); 124 } 125 126 #undef __FUNCT__ 127 #define __FUNCT__ "SNESLoad" 128 /*@C 129 SNESLoad - Loads a SNES that has been stored in binary with SNESView(). 130 131 Collective on PetscViewer 132 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() 137 138 Level: intermediate 139 140 Notes: 141 The type is determined by the data in the file, any type set into the SNES before this call is ignored. 142 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 151 152 .seealso: PetscViewerBinaryOpen(), SNESView(), MatLoad(), VecLoad() 153 @*/ 154 PetscErrorCode SNESLoad(SNES snes, PetscViewer viewer) 155 { 156 PetscErrorCode ierr; 157 PetscBool isbinary; 158 PetscInt classid; 159 char type[256]; 160 KSP ksp; 161 DM dm; 162 DMSNES dmsnes; 163 164 PetscFunctionBegin; 165 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 166 PetscValidHeaderSpecific(viewer,PETSC_VIEWER_CLASSID,2); 167 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr); 168 if (!isbinary) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Invalid viewer; open viewer with PetscViewerBinaryOpen()"); 169 170 ierr = PetscViewerBinaryRead(viewer,&classid,1,PETSC_INT);CHKERRQ(ierr); 171 if (classid != SNES_FILE_CLASSID) SETERRQ(((PetscObject)snes)->comm,PETSC_ERR_ARG_WRONG,"Not SNES next in file"); 172 ierr = PetscViewerBinaryRead(viewer,type,256,PETSC_CHAR);CHKERRQ(ierr); 173 ierr = SNESSetType(snes, type);CHKERRQ(ierr); 174 if (snes->ops->load) { 175 ierr = (*snes->ops->load)(snes,viewer);CHKERRQ(ierr); 176 } 177 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 178 ierr = DMGetDMSNES(dm,&dmsnes);CHKERRQ(ierr); 179 ierr = DMSNESLoad(dmsnes,viewer);CHKERRQ(ierr); 180 ierr = SNESGetKSP(snes,&ksp);CHKERRQ(ierr); 181 ierr = KSPLoad(ksp,viewer);CHKERRQ(ierr); 182 PetscFunctionReturn(0); 183 } 184 185 #undef __FUNCT__ 186 #define __FUNCT__ "SNESView" 187 /*@C 188 SNESView - Prints the SNES data structure. 189 190 Collective on SNES 191 192 Input Parameters: 193 + SNES - the SNES context 194 - viewer - visualization context 195 196 Options Database Key: 197 . -snes_view - Calls SNESView() at end of SNESSolve() 198 199 Notes: 200 The available visualization contexts include 201 + PETSC_VIEWER_STDOUT_SELF - standard output (default) 202 - PETSC_VIEWER_STDOUT_WORLD - synchronized standard 203 output where only the first processor opens 204 the file. All other processors send their 205 data to the first processor to print. 206 207 The user can open an alternative visualization context with 208 PetscViewerASCIIOpen() - output to a specified file. 209 210 Level: beginner 211 212 .keywords: SNES, view 213 214 .seealso: PetscViewerASCIIOpen() 215 @*/ 216 PetscErrorCode SNESView(SNES snes,PetscViewer viewer) 217 { 218 SNESKSPEW *kctx; 219 PetscErrorCode ierr; 220 KSP ksp; 221 SNESLineSearch linesearch; 222 PetscBool iascii,isstring,isbinary,isdraw; 223 DMSNES dmsnes; 224 225 PetscFunctionBegin; 226 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 227 if (!viewer) { 228 ierr = PetscViewerASCIIGetStdout(((PetscObject)snes)->comm,&viewer);CHKERRQ(ierr); 229 } 230 PetscValidHeaderSpecific(viewer,PETSC_VIEWER_CLASSID,2); 231 PetscCheckSameComm(snes,1,viewer,2); 232 233 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 234 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSTRING,&isstring);CHKERRQ(ierr); 235 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr); 236 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);CHKERRQ(ierr); 237 if (iascii) { 238 ierr = PetscObjectPrintClassNamePrefixType((PetscObject)snes,viewer,"SNES Object");CHKERRQ(ierr); 239 if (snes->ops->view) { 240 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 241 ierr = (*snes->ops->view)(snes,viewer);CHKERRQ(ierr); 242 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 243 } 244 ierr = PetscViewerASCIIPrintf(viewer," maximum iterations=%D, maximum function evaluations=%D\n",snes->max_its,snes->max_funcs);CHKERRQ(ierr); 245 ierr = PetscViewerASCIIPrintf(viewer," tolerances: relative=%G, absolute=%G, solution=%G\n",snes->rtol,snes->abstol,snes->stol);CHKERRQ(ierr); 246 ierr = PetscViewerASCIIPrintf(viewer," total number of linear solver iterations=%D\n",snes->linear_its);CHKERRQ(ierr); 247 ierr = PetscViewerASCIIPrintf(viewer," total number of function evaluations=%D\n",snes->nfuncs);CHKERRQ(ierr); 248 if (snes->gridsequence) { 249 ierr = PetscViewerASCIIPrintf(viewer," total number of grid sequence refinements=%D\n",snes->gridsequence);CHKERRQ(ierr); 250 } 251 if (snes->ksp_ewconv) { 252 kctx = (SNESKSPEW*)snes->kspconvctx; 253 if (kctx) { 254 ierr = PetscViewerASCIIPrintf(viewer," Eisenstat-Walker computation of KSP relative tolerance (version %D)\n",kctx->version);CHKERRQ(ierr); 255 ierr = PetscViewerASCIIPrintf(viewer," rtol_0=%G, rtol_max=%G, threshold=%G\n",kctx->rtol_0,kctx->rtol_max,kctx->threshold);CHKERRQ(ierr); 256 ierr = PetscViewerASCIIPrintf(viewer," gamma=%G, alpha=%G, alpha2=%G\n",kctx->gamma,kctx->alpha,kctx->alpha2);CHKERRQ(ierr); 257 } 258 } 259 if (snes->lagpreconditioner == -1) { 260 ierr = PetscViewerASCIIPrintf(viewer," Preconditioned is never rebuilt\n");CHKERRQ(ierr); 261 } else if (snes->lagpreconditioner > 1) { 262 ierr = PetscViewerASCIIPrintf(viewer," Preconditioned is rebuilt every %D new Jacobians\n",snes->lagpreconditioner);CHKERRQ(ierr); 263 } 264 if (snes->lagjacobian == -1) { 265 ierr = PetscViewerASCIIPrintf(viewer," Jacobian is never rebuilt\n");CHKERRQ(ierr); 266 } else if (snes->lagjacobian > 1) { 267 ierr = PetscViewerASCIIPrintf(viewer," Jacobian is rebuilt every %D SNES iterations\n",snes->lagjacobian);CHKERRQ(ierr); 268 } 269 } else if (isstring) { 270 const char *type; 271 ierr = SNESGetType(snes,&type);CHKERRQ(ierr); 272 ierr = PetscViewerStringSPrintf(viewer," %-3.3s",type);CHKERRQ(ierr); 273 } else if (isbinary) { 274 PetscInt classid = SNES_FILE_CLASSID; 275 MPI_Comm comm; 276 PetscMPIInt rank; 277 char type[256]; 278 279 ierr = PetscObjectGetComm((PetscObject)snes,&comm);CHKERRQ(ierr); 280 ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr); 281 if (!rank) { 282 ierr = PetscViewerBinaryWrite(viewer,&classid,1,PETSC_INT,PETSC_FALSE);CHKERRQ(ierr); 283 ierr = PetscStrncpy(type,((PetscObject)snes)->type_name,256);CHKERRQ(ierr); 284 ierr = PetscViewerBinaryWrite(viewer,type,256,PETSC_CHAR,PETSC_FALSE);CHKERRQ(ierr); 285 } 286 if (snes->ops->view) { 287 ierr = (*snes->ops->view)(snes,viewer);CHKERRQ(ierr); 288 } 289 } else if (isdraw) { 290 PetscDraw draw; 291 char str[36]; 292 PetscReal x,y,bottom,h; 293 294 ierr = PetscViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr); 295 ierr = PetscDrawGetCurrentPoint(draw,&x,&y);CHKERRQ(ierr); 296 ierr = PetscStrcpy(str,"SNES: ");CHKERRQ(ierr); 297 ierr = PetscStrcat(str,((PetscObject)snes)->type_name);CHKERRQ(ierr); 298 ierr = PetscDrawBoxedString(draw,x,y,PETSC_DRAW_BLUE,PETSC_DRAW_BLACK,str,NULL,&h);CHKERRQ(ierr); 299 bottom = y - h; 300 ierr = PetscDrawPushCurrentPoint(draw,x,bottom);CHKERRQ(ierr); 301 if (snes->ops->view) { 302 ierr = (*snes->ops->view)(snes,viewer);CHKERRQ(ierr); 303 } 304 } 305 if (snes->linesearch) { 306 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 307 ierr = SNESGetSNESLineSearch(snes, &linesearch);CHKERRQ(ierr); 308 ierr = SNESLineSearchView(linesearch, viewer);CHKERRQ(ierr); 309 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 310 } 311 if (snes->pc && snes->usespc) { 312 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 313 ierr = SNESView(snes->pc, viewer);CHKERRQ(ierr); 314 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 315 } 316 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 317 ierr = DMGetDMSNES(snes->dm,&dmsnes);CHKERRQ(ierr); 318 ierr = DMSNESView(dmsnes, viewer);CHKERRQ(ierr); 319 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 320 if (snes->usesksp) { 321 ierr = SNESGetKSP(snes,&ksp);CHKERRQ(ierr); 322 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 323 ierr = KSPView(ksp,viewer);CHKERRQ(ierr); 324 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 325 } 326 if (isdraw) { 327 PetscDraw draw; 328 ierr = PetscViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr); 329 ierr = PetscDrawPopCurrentPoint(draw);CHKERRQ(ierr); 330 } 331 PetscFunctionReturn(0); 332 } 333 334 /* 335 We retain a list of functions that also take SNES command 336 line options. These are called at the end SNESSetFromOptions() 337 */ 338 #define MAXSETFROMOPTIONS 5 339 static PetscInt numberofsetfromoptions; 340 static PetscErrorCode (*othersetfromoptions[MAXSETFROMOPTIONS])(SNES); 341 342 #undef __FUNCT__ 343 #define __FUNCT__ "SNESAddOptionsChecker" 344 /*@C 345 SNESAddOptionsChecker - Adds an additional function to check for SNES options. 346 347 Not Collective 348 349 Input Parameter: 350 . snescheck - function that checks for options 351 352 Level: developer 353 354 .seealso: SNESSetFromOptions() 355 @*/ 356 PetscErrorCode SNESAddOptionsChecker(PetscErrorCode (*snescheck)(SNES)) 357 { 358 PetscFunctionBegin; 359 if (numberofsetfromoptions >= MAXSETFROMOPTIONS) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE, "Too many options checkers, only %D allowed", MAXSETFROMOPTIONS); 360 othersetfromoptions[numberofsetfromoptions++] = snescheck; 361 PetscFunctionReturn(0); 362 } 363 364 extern PetscErrorCode SNESDefaultMatrixFreeCreate2(SNES,Vec,Mat*); 365 366 #undef __FUNCT__ 367 #define __FUNCT__ "SNESSetUpMatrixFree_Private" 368 static PetscErrorCode SNESSetUpMatrixFree_Private(SNES snes, PetscBool hasOperator, PetscInt version) 369 { 370 Mat J; 371 KSP ksp; 372 PC pc; 373 PetscBool match; 374 PetscErrorCode ierr; 375 376 PetscFunctionBegin; 377 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 378 379 if (!snes->vec_func && (snes->jacobian || snes->jacobian_pre)) { 380 Mat A = snes->jacobian, B = snes->jacobian_pre; 381 ierr = MatGetVecs(A ? A : B, NULL,&snes->vec_func);CHKERRQ(ierr); 382 } 383 384 if (version == 1) { 385 ierr = MatCreateSNESMF(snes,&J);CHKERRQ(ierr); 386 ierr = MatMFFDSetOptionsPrefix(J,((PetscObject)snes)->prefix);CHKERRQ(ierr); 387 ierr = MatSetFromOptions(J);CHKERRQ(ierr); 388 } else if (version == 2) { 389 if (!snes->vec_func) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"SNESSetFunction() must be called first"); 390 #if !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_REAL_SINGLE) && !defined(PETSC_USE_REAL___FLOAT128) 391 ierr = SNESDefaultMatrixFreeCreate2(snes,snes->vec_func,&J);CHKERRQ(ierr); 392 #else 393 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP, "matrix-free operator rutines (version 2)"); 394 #endif 395 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE, "matrix-free operator rutines, only version 1 and 2"); 396 397 ierr = PetscInfo1(snes,"Setting default matrix-free operator routines (version %D)\n", version);CHKERRQ(ierr); 398 if (hasOperator) { 399 400 /* This version replaces the user provided Jacobian matrix with a 401 matrix-free version but still employs the user-provided preconditioner matrix. */ 402 ierr = SNESSetJacobian(snes,J,0,0,0);CHKERRQ(ierr); 403 } else { 404 /* This version replaces both the user-provided Jacobian and the user- 405 provided preconditioner Jacobian with the default matrix free version. */ 406 if ((snes->pcside == PC_LEFT) && snes->pc) { 407 ierr = SNESSetJacobian(snes,J,0,0,0);CHKERRQ(ierr); 408 } else { 409 ierr = SNESSetJacobian(snes,J,J,MatMFFDComputeJacobian,0);CHKERRQ(ierr); 410 } 411 /* Force no preconditioner */ 412 ierr = SNESGetKSP(snes,&ksp);CHKERRQ(ierr); 413 ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr); 414 ierr = PetscObjectTypeCompare((PetscObject)pc,PCSHELL,&match);CHKERRQ(ierr); 415 if (!match) { 416 ierr = PetscInfo(snes,"Setting default matrix-free preconditioner routines\nThat is no preconditioner is being used\n");CHKERRQ(ierr); 417 ierr = PCSetType(pc,PCNONE);CHKERRQ(ierr); 418 } 419 } 420 ierr = MatDestroy(&J);CHKERRQ(ierr); 421 PetscFunctionReturn(0); 422 } 423 424 #undef __FUNCT__ 425 #define __FUNCT__ "DMRestrictHook_SNESVecSol" 426 static PetscErrorCode DMRestrictHook_SNESVecSol(DM dmfine,Mat Restrict,Vec Rscale,Mat Inject,DM dmcoarse,void *ctx) 427 { 428 SNES snes = (SNES)ctx; 429 PetscErrorCode ierr; 430 Vec Xfine,Xfine_named = NULL,Xcoarse; 431 432 PetscFunctionBegin; 433 if (PetscLogPrintInfo) { 434 PetscInt finelevel,coarselevel,fineclevel,coarseclevel; 435 ierr = DMGetRefineLevel(dmfine,&finelevel);CHKERRQ(ierr); 436 ierr = DMGetCoarsenLevel(dmfine,&fineclevel);CHKERRQ(ierr); 437 ierr = DMGetRefineLevel(dmcoarse,&coarselevel);CHKERRQ(ierr); 438 ierr = DMGetCoarsenLevel(dmcoarse,&coarseclevel);CHKERRQ(ierr); 439 ierr = PetscInfo4(dmfine,"Restricting SNES solution vector from level %D-%D to level %D-%D\n",finelevel,fineclevel,coarselevel,coarseclevel);CHKERRQ(ierr); 440 } 441 if (dmfine == snes->dm) Xfine = snes->vec_sol; 442 else { 443 ierr = DMGetNamedGlobalVector(dmfine,"SNESVecSol",&Xfine_named);CHKERRQ(ierr); 444 Xfine = Xfine_named; 445 } 446 ierr = DMGetNamedGlobalVector(dmcoarse,"SNESVecSol",&Xcoarse);CHKERRQ(ierr); 447 ierr = MatRestrict(Restrict,Xfine,Xcoarse);CHKERRQ(ierr); 448 ierr = VecPointwiseMult(Xcoarse,Xcoarse,Rscale);CHKERRQ(ierr); 449 ierr = DMRestoreNamedGlobalVector(dmcoarse,"SNESVecSol",&Xcoarse);CHKERRQ(ierr); 450 if (Xfine_named) {ierr = DMRestoreNamedGlobalVector(dmfine,"SNESVecSol",&Xfine_named);CHKERRQ(ierr);} 451 PetscFunctionReturn(0); 452 } 453 454 #undef __FUNCT__ 455 #define __FUNCT__ "DMCoarsenHook_SNESVecSol" 456 static PetscErrorCode DMCoarsenHook_SNESVecSol(DM dm,DM dmc,void *ctx) 457 { 458 PetscErrorCode ierr; 459 460 PetscFunctionBegin; 461 ierr = DMCoarsenHookAdd(dmc,DMCoarsenHook_SNESVecSol,DMRestrictHook_SNESVecSol,ctx);CHKERRQ(ierr); 462 PetscFunctionReturn(0); 463 } 464 465 #undef __FUNCT__ 466 #define __FUNCT__ "KSPComputeOperators_SNES" 467 /* This may be called to rediscretize the operator on levels of linear multigrid. The DM shuffle is so the user can 468 * safely call SNESGetDM() in their residual evaluation routine. */ 469 static PetscErrorCode KSPComputeOperators_SNES(KSP ksp,Mat A,Mat B,MatStructure *mstruct,void *ctx) 470 { 471 SNES snes = (SNES)ctx; 472 PetscErrorCode ierr; 473 Mat Asave = A,Bsave = B; 474 Vec X,Xnamed = NULL; 475 DM dmsave; 476 void *ctxsave; 477 PetscErrorCode (*jac)(SNES,Vec,Mat*,Mat*,MatStructure*,void*); 478 479 PetscFunctionBegin; 480 dmsave = snes->dm; 481 ierr = KSPGetDM(ksp,&snes->dm);CHKERRQ(ierr); 482 if (dmsave == snes->dm) X = snes->vec_sol; /* We are on the finest level */ 483 else { /* We are on a coarser level, this vec was initialized using a DM restrict hook */ 484 ierr = DMGetNamedGlobalVector(snes->dm,"SNESVecSol",&Xnamed);CHKERRQ(ierr); 485 X = Xnamed; 486 ierr = SNESGetJacobian(snes,NULL,NULL,&jac,&ctxsave);CHKERRQ(ierr); 487 /* If the DM's don't match up, the MatFDColoring context needed for the jacobian won't match up either -- fixit. */ 488 if (jac == SNESDefaultComputeJacobianColor) { 489 ierr = SNESSetJacobian(snes,NULL,NULL,SNESDefaultComputeJacobianColor,0);CHKERRQ(ierr); 490 } 491 } 492 /* put the previous context back */ 493 494 ierr = SNESComputeJacobian(snes,X,&A,&B,mstruct);CHKERRQ(ierr); 495 if (snes->dm != dmsave && jac == SNESDefaultComputeJacobianColor) { 496 ierr = SNESSetJacobian(snes,NULL,NULL,jac,ctxsave);CHKERRQ(ierr); 497 } 498 499 if (A != Asave || B != Bsave) SETERRQ(((PetscObject)snes)->comm,PETSC_ERR_SUP,"No support for changing matrices at this time"); 500 if (Xnamed) { 501 ierr = DMRestoreNamedGlobalVector(snes->dm,"SNESVecSol",&Xnamed);CHKERRQ(ierr); 502 } 503 snes->dm = dmsave; 504 PetscFunctionReturn(0); 505 } 506 507 #undef __FUNCT__ 508 #define __FUNCT__ "SNESSetUpMatrices" 509 /*@ 510 SNESSetUpMatrices - ensures that matrices are available for SNES, to be called by SNESSetUp_XXX() 511 512 Collective 513 514 Input Arguments: 515 . snes - snes to configure 516 517 Level: developer 518 519 .seealso: SNESSetUp() 520 @*/ 521 PetscErrorCode SNESSetUpMatrices(SNES snes) 522 { 523 PetscErrorCode ierr; 524 DM dm; 525 DMSNES sdm; 526 527 PetscFunctionBegin; 528 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 529 ierr = DMGetDMSNES(dm,&sdm);CHKERRQ(ierr); 530 if (!sdm->ops->computejacobian) SETERRQ(((PetscObject)snes)->comm,PETSC_ERR_PLIB,"DMSNES not properly configured"); 531 else if (!snes->jacobian && snes->mf) { 532 Mat J; 533 void *functx; 534 ierr = MatCreateSNESMF(snes,&J);CHKERRQ(ierr); 535 ierr = MatMFFDSetOptionsPrefix(J,((PetscObject)snes)->prefix);CHKERRQ(ierr); 536 ierr = MatSetFromOptions(J);CHKERRQ(ierr); 537 ierr = SNESGetFunction(snes,NULL,NULL,&functx);CHKERRQ(ierr); 538 ierr = SNESSetJacobian(snes,J,J,0,0);CHKERRQ(ierr); 539 ierr = MatDestroy(&J);CHKERRQ(ierr); 540 } else if (snes->mf_operator && !snes->jacobian_pre && !snes->jacobian) { 541 Mat J,B; 542 ierr = MatCreateSNESMF(snes,&J);CHKERRQ(ierr); 543 ierr = MatMFFDSetOptionsPrefix(J,((PetscObject)snes)->prefix);CHKERRQ(ierr); 544 ierr = MatSetFromOptions(J);CHKERRQ(ierr); 545 ierr = DMCreateMatrix(snes->dm,MATAIJ,&B);CHKERRQ(ierr); 546 /* sdm->computejacobian was already set to reach here */ 547 ierr = SNESSetJacobian(snes,J,B,NULL,NULL);CHKERRQ(ierr); 548 ierr = MatDestroy(&J);CHKERRQ(ierr); 549 ierr = MatDestroy(&B);CHKERRQ(ierr); 550 } else if (!snes->jacobian_pre) { 551 Mat J,B; 552 J = snes->jacobian; 553 ierr = DMCreateMatrix(snes->dm,MATAIJ,&B);CHKERRQ(ierr); 554 ierr = SNESSetJacobian(snes,J ? J : B,B,NULL,NULL);CHKERRQ(ierr); 555 ierr = MatDestroy(&B);CHKERRQ(ierr); 556 } 557 { 558 KSP ksp; 559 ierr = SNESGetKSP(snes,&ksp);CHKERRQ(ierr); 560 ierr = KSPSetComputeOperators(ksp,KSPComputeOperators_SNES,snes);CHKERRQ(ierr); 561 ierr = DMCoarsenHookAdd(snes->dm,DMCoarsenHook_SNESVecSol,DMRestrictHook_SNESVecSol,snes);CHKERRQ(ierr); 562 } 563 PetscFunctionReturn(0); 564 } 565 566 #undef __FUNCT__ 567 #define __FUNCT__ "SNESSetFromOptions" 568 /*@ 569 SNESSetFromOptions - Sets various SNES and KSP parameters from user options. 570 571 Collective on SNES 572 573 Input Parameter: 574 . snes - the SNES context 575 576 Options Database Keys: 577 + -snes_type <type> - ls, tr, ngmres, ncg, richardson, qn, vi, fas 578 . -snes_stol - convergence tolerance in terms of the norm 579 of the change in the solution between steps 580 . -snes_atol <abstol> - absolute tolerance of residual norm 581 . -snes_rtol <rtol> - relative decrease in tolerance norm from initial 582 . -snes_max_it <max_it> - maximum number of iterations 583 . -snes_max_funcs <max_funcs> - maximum number of function evaluations 584 . -snes_max_fail <max_fail> - maximum number of line search failures allowed before stopping, default is none 585 . -snes_max_linear_solve_fail - number of linear solver failures before SNESSolve() stops 586 . -snes_lag_preconditioner <lag> - how often preconditioner is rebuilt (use -1 to never rebuild) 587 . -snes_lag_jacobian <lag> - how often Jacobian is rebuilt (use -1 to never rebuild) 588 . -snes_trtol <trtol> - trust region tolerance 589 . -snes_no_convergence_test - skip convergence test in nonlinear 590 solver; hence iterations will continue until max_it 591 or some other criterion is reached. Saves expense 592 of convergence test 593 . -snes_monitor <optional filename> - prints residual norm at each iteration. if no 594 filename given prints to stdout 595 . -snes_monitor_solution - plots solution at each iteration 596 . -snes_monitor_residual - plots residual (not its norm) at each iteration 597 . -snes_monitor_solution_update - plots update to solution at each iteration 598 . -snes_monitor_lg_residualnorm - plots residual norm at each iteration 599 . -snes_monitor_lg_range - plots residual norm at each iteration 600 . -snes_fd - use finite differences to compute Jacobian; very slow, only for testing 601 . -snes_fd_color - use finite differences with coloring to compute Jacobian 602 . -snes_mf_ksp_monitor - if using matrix-free multiply then print h at each KSP iteration 603 - -snes_converged_reason - print the reason for convergence/divergence after each solve 604 605 Options Database for Eisenstat-Walker method: 606 + -snes_ksp_ew - use Eisenstat-Walker method for determining linear system convergence 607 . -snes_ksp_ew_version ver - version of Eisenstat-Walker method 608 . -snes_ksp_ew_rtol0 <rtol0> - Sets rtol0 609 . -snes_ksp_ew_rtolmax <rtolmax> - Sets rtolmax 610 . -snes_ksp_ew_gamma <gamma> - Sets gamma 611 . -snes_ksp_ew_alpha <alpha> - Sets alpha 612 . -snes_ksp_ew_alpha2 <alpha2> - Sets alpha2 613 - -snes_ksp_ew_threshold <threshold> - Sets threshold 614 615 Notes: 616 To see all options, run your program with the -help option or consult 617 the <A href="../../docs/manual.pdf#nameddest=ch_snes">SNES chapter of the users manual</A>. 618 619 Level: beginner 620 621 .keywords: SNES, nonlinear, set, options, database 622 623 .seealso: SNESSetOptionsPrefix() 624 @*/ 625 PetscErrorCode SNESSetFromOptions(SNES snes) 626 { 627 PetscBool flg,pcset; 628 PetscInt i,indx,lag,grids; 629 MatStructure matflag; 630 const char *deft = SNESNEWTONLS; 631 const char *convtests[] = {"default","skip"}; 632 SNESKSPEW *kctx = NULL; 633 char type[256], monfilename[PETSC_MAX_PATH_LEN]; 634 PetscViewer monviewer; 635 PetscErrorCode ierr; 636 PCSide pcside; 637 const char *optionsprefix; 638 639 PetscFunctionBegin; 640 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 641 642 if (!SNESRegisterAllCalled) {ierr = SNESRegisterAll(NULL);CHKERRQ(ierr);} 643 ierr = PetscObjectOptionsBegin((PetscObject)snes);CHKERRQ(ierr); 644 if (((PetscObject)snes)->type_name) deft = ((PetscObject)snes)->type_name; 645 ierr = PetscOptionsList("-snes_type","Nonlinear solver method","SNESSetType",SNESList,deft,type,256,&flg);CHKERRQ(ierr); 646 if (flg) { 647 ierr = SNESSetType(snes,type);CHKERRQ(ierr); 648 } else if (!((PetscObject)snes)->type_name) { 649 ierr = SNESSetType(snes,deft);CHKERRQ(ierr); 650 } 651 /* not used here, but called so will go into help messaage */ 652 ierr = PetscOptionsName("-snes_view","Print detailed information on solver used","SNESView",0);CHKERRQ(ierr); 653 654 ierr = PetscOptionsReal("-snes_stol","Stop if step length less than","SNESSetTolerances",snes->stol,&snes->stol,0);CHKERRQ(ierr); 655 ierr = PetscOptionsReal("-snes_atol","Stop if function norm less than","SNESSetTolerances",snes->abstol,&snes->abstol,0);CHKERRQ(ierr); 656 657 ierr = PetscOptionsReal("-snes_rtol","Stop if decrease in function norm less than","SNESSetTolerances",snes->rtol,&snes->rtol,0);CHKERRQ(ierr); 658 ierr = PetscOptionsInt("-snes_max_it","Maximum iterations","SNESSetTolerances",snes->max_its,&snes->max_its,NULL);CHKERRQ(ierr); 659 ierr = PetscOptionsInt("-snes_max_funcs","Maximum function evaluations","SNESSetTolerances",snes->max_funcs,&snes->max_funcs,NULL);CHKERRQ(ierr); 660 ierr = PetscOptionsInt("-snes_max_fail","Maximum nonlinear step failures","SNESSetMaxNonlinearStepFailures",snes->maxFailures,&snes->maxFailures,NULL);CHKERRQ(ierr); 661 ierr = PetscOptionsInt("-snes_max_linear_solve_fail","Maximum failures in linear solves allowed","SNESSetMaxLinearSolveFailures",snes->maxLinearSolveFailures,&snes->maxLinearSolveFailures,NULL);CHKERRQ(ierr); 662 ierr = PetscOptionsBool("-snes_error_if_not_converged","Generate error if solver does not converge","SNESSetErrorIfNotConverged",snes->errorifnotconverged,&snes->errorifnotconverged,NULL);CHKERRQ(ierr); 663 664 ierr = PetscOptionsInt("-snes_lag_preconditioner","How often to rebuild preconditioner","SNESSetLagPreconditioner",snes->lagpreconditioner,&lag,&flg);CHKERRQ(ierr); 665 if (flg) { 666 ierr = SNESSetLagPreconditioner(snes,lag);CHKERRQ(ierr); 667 } 668 ierr = PetscOptionsInt("-snes_lag_jacobian","How often to rebuild Jacobian","SNESSetLagJacobian",snes->lagjacobian,&lag,&flg);CHKERRQ(ierr); 669 if (flg) { 670 ierr = SNESSetLagJacobian(snes,lag);CHKERRQ(ierr); 671 } 672 ierr = PetscOptionsInt("-snes_grid_sequence","Use grid sequencing to generate initial guess","SNESSetGridSequence",snes->gridsequence,&grids,&flg);CHKERRQ(ierr); 673 if (flg) { 674 ierr = SNESSetGridSequence(snes,grids);CHKERRQ(ierr); 675 } 676 677 ierr = PetscOptionsEList("-snes_convergence_test","Convergence test","SNESSetConvergenceTest",convtests,2,"default",&indx,&flg);CHKERRQ(ierr); 678 if (flg) { 679 switch (indx) { 680 case 0: ierr = SNESSetConvergenceTest(snes,SNESDefaultConverged,NULL,NULL);CHKERRQ(ierr); break; 681 case 1: ierr = SNESSetConvergenceTest(snes,SNESSkipConverged,NULL,NULL);CHKERRQ(ierr); break; 682 } 683 } 684 685 ierr = PetscOptionsBool("-snes_converged_reason","Print reason for converged or diverged","SNESSolve",snes->printreason,&snes->printreason,NULL);CHKERRQ(ierr); 686 687 ierr = PetscOptionsEList("-snes_norm_type","SNES Norm type","SNESSetNormType",SNESNormTypes,5,"function",&indx,&flg);CHKERRQ(ierr); 688 if (flg) { ierr = SNESSetNormType(snes,(SNESNormType)indx);CHKERRQ(ierr); } 689 690 kctx = (SNESKSPEW*)snes->kspconvctx; 691 692 ierr = PetscOptionsBool("-snes_ksp_ew","Use Eisentat-Walker linear system convergence test","SNESKSPSetUseEW",snes->ksp_ewconv,&snes->ksp_ewconv,NULL);CHKERRQ(ierr); 693 694 ierr = PetscOptionsInt("-snes_ksp_ew_version","Version 1, 2 or 3","SNESKSPSetParametersEW",kctx->version,&kctx->version,0);CHKERRQ(ierr); 695 ierr = PetscOptionsReal("-snes_ksp_ew_rtol0","0 <= rtol0 < 1","SNESKSPSetParametersEW",kctx->rtol_0,&kctx->rtol_0,0);CHKERRQ(ierr); 696 ierr = PetscOptionsReal("-snes_ksp_ew_rtolmax","0 <= rtolmax < 1","SNESKSPSetParametersEW",kctx->rtol_max,&kctx->rtol_max,0);CHKERRQ(ierr); 697 ierr = PetscOptionsReal("-snes_ksp_ew_gamma","0 <= gamma <= 1","SNESKSPSetParametersEW",kctx->gamma,&kctx->gamma,0);CHKERRQ(ierr); 698 ierr = PetscOptionsReal("-snes_ksp_ew_alpha","1 < alpha <= 2","SNESKSPSetParametersEW",kctx->alpha,&kctx->alpha,0);CHKERRQ(ierr); 699 ierr = PetscOptionsReal("-snes_ksp_ew_alpha2","alpha2","SNESKSPSetParametersEW",kctx->alpha2,&kctx->alpha2,0);CHKERRQ(ierr); 700 ierr = PetscOptionsReal("-snes_ksp_ew_threshold","0 < threshold < 1","SNESKSPSetParametersEW",kctx->threshold,&kctx->threshold,0);CHKERRQ(ierr); 701 702 flg = PETSC_FALSE; 703 ierr = PetscOptionsBool("-snes_check_jacobian","Check each Jacobian with a differenced one","SNESUpdateCheckJacobian",flg,&flg,NULL);CHKERRQ(ierr); 704 if (flg) { 705 ierr = SNESSetUpdate(snes,SNESUpdateCheckJacobian);CHKERRQ(ierr); 706 } 707 708 flg = PETSC_FALSE; 709 ierr = PetscOptionsBool("-snes_monitor_cancel","Remove all monitors","SNESMonitorCancel",flg,&flg,NULL);CHKERRQ(ierr); 710 if (flg) {ierr = SNESMonitorCancel(snes);CHKERRQ(ierr);} 711 712 ierr = PetscOptionsString("-snes_monitor","Monitor norm of function","SNESMonitorSet","stdout",monfilename,PETSC_MAX_PATH_LEN,&flg);CHKERRQ(ierr); 713 if (flg) { 714 ierr = PetscViewerASCIIOpen(((PetscObject)snes)->comm,monfilename,&monviewer);CHKERRQ(ierr); 715 ierr = SNESMonitorSet(snes,SNESMonitorDefault,monviewer,(PetscErrorCode (*)(void**))PetscViewerDestroy);CHKERRQ(ierr); 716 } 717 718 ierr = PetscOptionsString("-snes_monitor_range","Monitor range of elements of function","SNESMonitorSet","stdout",monfilename,PETSC_MAX_PATH_LEN,&flg);CHKERRQ(ierr); 719 if (flg) { 720 ierr = SNESMonitorSet(snes,SNESMonitorRange,0,0);CHKERRQ(ierr); 721 } 722 723 ierr = PetscOptionsString("-snes_ratiomonitor","Monitor ratios of norms of function","SNESMonitorSetRatio","stdout",monfilename,PETSC_MAX_PATH_LEN,&flg);CHKERRQ(ierr); 724 if (flg) { 725 ierr = PetscViewerASCIIOpen(((PetscObject)snes)->comm,monfilename,&monviewer);CHKERRQ(ierr); 726 ierr = SNESMonitorSetRatio(snes,monviewer);CHKERRQ(ierr); 727 } 728 729 ierr = PetscOptionsString("-snes_monitor_short","Monitor norm of function (fewer digits)","SNESMonitorSet","stdout",monfilename,PETSC_MAX_PATH_LEN,&flg);CHKERRQ(ierr); 730 if (flg) { 731 ierr = PetscViewerASCIIOpen(((PetscObject)snes)->comm,monfilename,&monviewer);CHKERRQ(ierr); 732 ierr = SNESMonitorSet(snes,SNESMonitorDefaultShort,monviewer,(PetscErrorCode (*)(void**))PetscViewerDestroy);CHKERRQ(ierr); 733 } 734 735 ierr = PetscOptionsString("-snes_monitor_python","Use Python function","SNESMonitorSet",0,monfilename,PETSC_MAX_PATH_LEN,&flg);CHKERRQ(ierr); 736 if (flg) {ierr = PetscPythonMonitorSet((PetscObject)snes,monfilename);CHKERRQ(ierr);} 737 738 flg = PETSC_FALSE; 739 ierr = PetscOptionsBool("-snes_monitor_solution","Plot solution at each iteration","SNESMonitorSolution",flg,&flg,NULL);CHKERRQ(ierr); 740 if (flg) {ierr = SNESMonitorSet(snes,SNESMonitorSolution,0,0);CHKERRQ(ierr);} 741 flg = PETSC_FALSE; 742 ierr = PetscOptionsBool("-snes_monitor_solution_update","Plot correction at each iteration","SNESMonitorSolutionUpdate",flg,&flg,NULL);CHKERRQ(ierr); 743 if (flg) {ierr = SNESMonitorSet(snes,SNESMonitorSolutionUpdate,0,0);CHKERRQ(ierr);} 744 flg = PETSC_FALSE; 745 ierr = PetscOptionsBool("-snes_monitor_residual","Plot residual at each iteration","SNESMonitorResidual",flg,&flg,NULL);CHKERRQ(ierr); 746 if (flg) {ierr = SNESMonitorSet(snes,SNESMonitorResidual,0,0);CHKERRQ(ierr);} 747 flg = PETSC_FALSE; 748 ierr = PetscOptionsBool("-snes_monitor_lg_residualnorm","Plot function norm at each iteration","SNESMonitorLGResidualNorm",flg,&flg,NULL);CHKERRQ(ierr); 749 if (flg) { 750 PetscDrawLG ctx; 751 752 ierr = SNESMonitorLGCreate(0,0,PETSC_DECIDE,PETSC_DECIDE,300,300,&ctx);CHKERRQ(ierr); 753 ierr = SNESMonitorSet(snes,SNESMonitorLGResidualNorm,ctx,(PetscErrorCode (*)(void**))SNESMonitorLGDestroy);CHKERRQ(ierr); 754 } 755 flg = PETSC_FALSE; 756 ierr = PetscOptionsBool("-snes_monitor_lg_range","Plot function range at each iteration","SNESMonitorLGRange",flg,&flg,NULL);CHKERRQ(ierr); 757 if (flg) { 758 PetscViewer ctx; 759 760 ierr = PetscViewerDrawOpen(((PetscObject)snes)->comm,0,0,PETSC_DECIDE,PETSC_DECIDE,300,300,&ctx);CHKERRQ(ierr); 761 ierr = SNESMonitorSet(snes,SNESMonitorLGRange,ctx,(PetscErrorCode (*)(void**))PetscViewerDestroy);CHKERRQ(ierr); 762 } 763 764 flg = PETSC_FALSE; 765 ierr = PetscOptionsBool("-snes_monitor_jacupdate_spectrum","Print the change in the spectrum of the Jacobian","SNESMonitorJacUpdateSpectrum",flg,&flg,NULL);CHKERRQ(ierr); 766 if (flg) {ierr = SNESMonitorSet(snes,SNESMonitorJacUpdateSpectrum,0,0);CHKERRQ(ierr);} 767 768 flg = PETSC_FALSE; 769 ierr = PetscOptionsBool("-snes_fd","Use finite differences (slow) to compute Jacobian","SNESDefaultComputeJacobian",flg,&flg,NULL);CHKERRQ(ierr); 770 if (flg) { 771 void *functx; 772 ierr = SNESGetFunction(snes,NULL,NULL,&functx);CHKERRQ(ierr); 773 ierr = SNESSetJacobian(snes,snes->jacobian,snes->jacobian_pre,SNESDefaultComputeJacobian,functx);CHKERRQ(ierr); 774 ierr = PetscInfo(snes,"Setting default finite difference Jacobian matrix\n");CHKERRQ(ierr); 775 } 776 777 flg = PETSC_FALSE; 778 ierr = PetscOptionsBool("-snes_fd_function","Use finite differences (slow) to compute function from user objective","SNESDefaultObjectiveComputeFunctionFD",flg,&flg,NULL);CHKERRQ(ierr); 779 if (flg) { 780 ierr = SNESSetFunction(snes,NULL,SNESDefaultObjectiveComputeFunctionFD,NULL);CHKERRQ(ierr); 781 } 782 783 flg = PETSC_FALSE; 784 ierr = PetscOptionsBool("-snes_fd_color","Use finite differences with coloring to compute Jacobian","SNESDefaultComputeJacobianColor",flg,&flg,NULL);CHKERRQ(ierr); 785 if (flg) { 786 void *functx; 787 ierr = SNESGetFunction(snes,NULL,NULL,&functx);CHKERRQ(ierr); 788 ierr = SNESSetJacobian(snes,snes->jacobian,snes->jacobian_pre,SNESDefaultComputeJacobianColor,0);CHKERRQ(ierr); 789 ierr = PetscInfo(snes,"Setting default finite difference coloring Jacobian matrix\n");CHKERRQ(ierr); 790 } 791 792 flg = PETSC_FALSE; 793 ierr = PetscOptionsBool("-snes_mf_operator","Use a Matrix-Free Jacobian with user-provided preconditioner matrix","MatCreateSNESMF",PETSC_FALSE,&snes->mf_operator,&flg);CHKERRQ(ierr); 794 if (flg && snes->mf_operator) { 795 snes->mf_operator = PETSC_TRUE; 796 snes->mf = PETSC_TRUE; 797 } 798 flg = PETSC_FALSE; 799 ierr = PetscOptionsBool("-snes_mf","Use a Matrix-Free Jacobian with no preconditioner matrix","MatCreateSNESMF",PETSC_FALSE,&snes->mf,&flg);CHKERRQ(ierr); 800 if (!flg && snes->mf_operator) snes->mf = PETSC_TRUE; 801 ierr = PetscOptionsInt("-snes_mf_version","Matrix-Free routines version 1 or 2","None",snes->mf_version,&snes->mf_version,0);CHKERRQ(ierr); 802 803 flg = PETSC_FALSE; 804 ierr = SNESGetPCSide(snes,&pcside);CHKERRQ(ierr); 805 ierr = PetscOptionsEnum("-snes_npc_side","SNES nonlinear preconditioner side","SNESSetPCSide",PCSides,(PetscEnum)pcside,(PetscEnum*)&pcside,&flg);CHKERRQ(ierr); 806 if (flg) {ierr = SNESSetPCSide(snes,pcside);CHKERRQ(ierr);} 807 808 for (i = 0; i < numberofsetfromoptions; i++) { 809 ierr = (*othersetfromoptions[i])(snes);CHKERRQ(ierr); 810 } 811 812 if (snes->ops->setfromoptions) { 813 ierr = (*snes->ops->setfromoptions)(snes);CHKERRQ(ierr); 814 } 815 816 /* process any options handlers added with PetscObjectAddOptionsHandler() */ 817 ierr = PetscObjectProcessOptionsHandlers((PetscObject)snes);CHKERRQ(ierr); 818 ierr = PetscOptionsEnd();CHKERRQ(ierr); 819 820 if (!snes->ksp) {ierr = SNESGetKSP(snes,&snes->ksp);CHKERRQ(ierr);} 821 ierr = KSPGetOperators(snes->ksp,NULL,NULL,&matflag);CHKERRQ(ierr); 822 ierr = KSPSetOperators(snes->ksp,snes->jacobian,snes->jacobian_pre,matflag);CHKERRQ(ierr); 823 ierr = KSPSetFromOptions(snes->ksp);CHKERRQ(ierr); 824 825 if (!snes->linesearch) { 826 ierr = SNESGetSNESLineSearch(snes, &snes->linesearch);CHKERRQ(ierr); 827 } 828 ierr = SNESLineSearchSetFromOptions(snes->linesearch);CHKERRQ(ierr); 829 830 /* if someone has set the SNES PC type, create it. */ 831 ierr = SNESGetOptionsPrefix(snes, &optionsprefix);CHKERRQ(ierr); 832 ierr = PetscOptionsHasName(optionsprefix, "-npc_snes_type", &pcset);CHKERRQ(ierr); 833 if (pcset && (!snes->pc)) { 834 ierr = SNESGetPC(snes, &snes->pc);CHKERRQ(ierr); 835 } 836 PetscFunctionReturn(0); 837 } 838 839 #undef __FUNCT__ 840 #define __FUNCT__ "SNESSetComputeApplicationContext" 841 /*@ 842 SNESSetComputeApplicationContext - Sets an optional function to compute a user-defined context for 843 the nonlinear solvers. 844 845 Logically Collective on SNES 846 847 Input Parameters: 848 + snes - the SNES context 849 . compute - function to compute the context 850 - destroy - function to destroy the context 851 852 Level: intermediate 853 854 .keywords: SNES, nonlinear, set, application, context 855 856 .seealso: SNESGetApplicationContext(), SNESSetComputeApplicationContext(), SNESGetApplicationContext() 857 @*/ 858 PetscErrorCode SNESSetComputeApplicationContext(SNES snes,PetscErrorCode (*compute)(SNES,void**),PetscErrorCode (*destroy)(void**)) 859 { 860 PetscFunctionBegin; 861 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 862 snes->ops->usercompute = compute; 863 snes->ops->userdestroy = destroy; 864 PetscFunctionReturn(0); 865 } 866 867 #undef __FUNCT__ 868 #define __FUNCT__ "SNESSetApplicationContext" 869 /*@ 870 SNESSetApplicationContext - Sets the optional user-defined context for 871 the nonlinear solvers. 872 873 Logically Collective on SNES 874 875 Input Parameters: 876 + snes - the SNES context 877 - usrP - optional user context 878 879 Level: intermediate 880 881 .keywords: SNES, nonlinear, set, application, context 882 883 .seealso: SNESGetApplicationContext() 884 @*/ 885 PetscErrorCode SNESSetApplicationContext(SNES snes,void *usrP) 886 { 887 PetscErrorCode ierr; 888 KSP ksp; 889 890 PetscFunctionBegin; 891 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 892 ierr = SNESGetKSP(snes,&ksp);CHKERRQ(ierr); 893 ierr = KSPSetApplicationContext(ksp,usrP);CHKERRQ(ierr); 894 snes->user = usrP; 895 PetscFunctionReturn(0); 896 } 897 898 #undef __FUNCT__ 899 #define __FUNCT__ "SNESGetApplicationContext" 900 /*@ 901 SNESGetApplicationContext - Gets the user-defined context for the 902 nonlinear solvers. 903 904 Not Collective 905 906 Input Parameter: 907 . snes - SNES context 908 909 Output Parameter: 910 . usrP - user context 911 912 Level: intermediate 913 914 .keywords: SNES, nonlinear, get, application, context 915 916 .seealso: SNESSetApplicationContext() 917 @*/ 918 PetscErrorCode SNESGetApplicationContext(SNES snes,void *usrP) 919 { 920 PetscFunctionBegin; 921 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 922 *(void**)usrP = snes->user; 923 PetscFunctionReturn(0); 924 } 925 926 #undef __FUNCT__ 927 #define __FUNCT__ "SNESGetIterationNumber" 928 /*@ 929 SNESGetIterationNumber - Gets the number of nonlinear iterations completed 930 at this time. 931 932 Not Collective 933 934 Input Parameter: 935 . snes - SNES context 936 937 Output Parameter: 938 . iter - iteration number 939 940 Notes: 941 For example, during the computation of iteration 2 this would return 1. 942 943 This is useful for using lagged Jacobians (where one does not recompute the 944 Jacobian at each SNES iteration). For example, the code 945 .vb 946 ierr = SNESGetIterationNumber(snes,&it); 947 if (!(it % 2)) { 948 [compute Jacobian here] 949 } 950 .ve 951 can be used in your ComputeJacobian() function to cause the Jacobian to be 952 recomputed every second SNES iteration. 953 954 Level: intermediate 955 956 .keywords: SNES, nonlinear, get, iteration, number, 957 958 .seealso: SNESGetFunctionNorm(), SNESGetLinearSolveIterations() 959 @*/ 960 PetscErrorCode SNESGetIterationNumber(SNES snes,PetscInt *iter) 961 { 962 PetscFunctionBegin; 963 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 964 PetscValidIntPointer(iter,2); 965 *iter = snes->iter; 966 PetscFunctionReturn(0); 967 } 968 969 #undef __FUNCT__ 970 #define __FUNCT__ "SNESSetIterationNumber" 971 /*@ 972 SNESSetIterationNumber - Sets the current iteration number. 973 974 Not Collective 975 976 Input Parameter: 977 . snes - SNES context 978 . iter - iteration number 979 980 Level: developer 981 982 .keywords: SNES, nonlinear, set, iteration, number, 983 984 .seealso: SNESGetFunctionNorm(), SNESGetLinearSolveIterations() 985 @*/ 986 PetscErrorCode SNESSetIterationNumber(SNES snes,PetscInt iter) 987 { 988 PetscErrorCode ierr; 989 990 PetscFunctionBegin; 991 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 992 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 993 snes->iter = iter; 994 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 995 PetscFunctionReturn(0); 996 } 997 998 #undef __FUNCT__ 999 #define __FUNCT__ "SNESGetFunctionNorm" 1000 /*@ 1001 SNESGetFunctionNorm - Gets the norm of the current function that was set 1002 with SNESSSetFunction(). 1003 1004 Collective on SNES 1005 1006 Input Parameter: 1007 . snes - SNES context 1008 1009 Output Parameter: 1010 . fnorm - 2-norm of function 1011 1012 Level: intermediate 1013 1014 .keywords: SNES, nonlinear, get, function, norm 1015 1016 .seealso: SNESGetFunction(), SNESGetIterationNumber(), SNESGetLinearSolveIterations() 1017 @*/ 1018 PetscErrorCode SNESGetFunctionNorm(SNES snes,PetscReal *fnorm) 1019 { 1020 PetscFunctionBegin; 1021 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1022 PetscValidScalarPointer(fnorm,2); 1023 *fnorm = snes->norm; 1024 PetscFunctionReturn(0); 1025 } 1026 1027 1028 #undef __FUNCT__ 1029 #define __FUNCT__ "SNESSetFunctionNorm" 1030 /*@ 1031 SNESSetFunctionNorm - Sets the 2-norm of the current function computed using VecNorm(). 1032 1033 Collective on SNES 1034 1035 Input Parameter: 1036 . snes - SNES context 1037 . fnorm - 2-norm of function 1038 1039 Level: developer 1040 1041 .keywords: SNES, nonlinear, set, function, norm 1042 1043 .seealso: SNESSetFunction(), SNESSetIterationNumber(), VecNorm(). 1044 @*/ 1045 PetscErrorCode SNESSetFunctionNorm(SNES snes,PetscReal fnorm) 1046 { 1047 1048 PetscErrorCode ierr; 1049 1050 PetscFunctionBegin; 1051 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1052 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 1053 snes->norm = fnorm; 1054 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 1055 PetscFunctionReturn(0); 1056 } 1057 1058 #undef __FUNCT__ 1059 #define __FUNCT__ "SNESGetNonlinearStepFailures" 1060 /*@ 1061 SNESGetNonlinearStepFailures - Gets the number of unsuccessful steps 1062 attempted by the nonlinear solver. 1063 1064 Not Collective 1065 1066 Input Parameter: 1067 . snes - SNES context 1068 1069 Output Parameter: 1070 . nfails - number of unsuccessful steps attempted 1071 1072 Notes: 1073 This counter is reset to zero for each successive call to SNESSolve(). 1074 1075 Level: intermediate 1076 1077 .keywords: SNES, nonlinear, get, number, unsuccessful, steps 1078 1079 .seealso: SNESGetMaxLinearSolveFailures(), SNESGetLinearSolveIterations(), SNESSetMaxLinearSolveFailures(), SNESGetLinearSolveFailures(), 1080 SNESSetMaxNonlinearStepFailures(), SNESGetMaxNonlinearStepFailures() 1081 @*/ 1082 PetscErrorCode SNESGetNonlinearStepFailures(SNES snes,PetscInt *nfails) 1083 { 1084 PetscFunctionBegin; 1085 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1086 PetscValidIntPointer(nfails,2); 1087 *nfails = snes->numFailures; 1088 PetscFunctionReturn(0); 1089 } 1090 1091 #undef __FUNCT__ 1092 #define __FUNCT__ "SNESSetMaxNonlinearStepFailures" 1093 /*@ 1094 SNESSetMaxNonlinearStepFailures - Sets the maximum number of unsuccessful steps 1095 attempted by the nonlinear solver before it gives up. 1096 1097 Not Collective 1098 1099 Input Parameters: 1100 + snes - SNES context 1101 - maxFails - maximum of unsuccessful steps 1102 1103 Level: intermediate 1104 1105 .keywords: SNES, nonlinear, set, maximum, unsuccessful, steps 1106 1107 .seealso: SNESGetMaxLinearSolveFailures(), SNESGetLinearSolveIterations(), SNESSetMaxLinearSolveFailures(), SNESGetLinearSolveFailures(), 1108 SNESGetMaxNonlinearStepFailures(), SNESGetNonlinearStepFailures() 1109 @*/ 1110 PetscErrorCode SNESSetMaxNonlinearStepFailures(SNES snes, PetscInt maxFails) 1111 { 1112 PetscFunctionBegin; 1113 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1114 snes->maxFailures = maxFails; 1115 PetscFunctionReturn(0); 1116 } 1117 1118 #undef __FUNCT__ 1119 #define __FUNCT__ "SNESGetMaxNonlinearStepFailures" 1120 /*@ 1121 SNESGetMaxNonlinearStepFailures - Gets the maximum number of unsuccessful steps 1122 attempted by the nonlinear solver before it gives up. 1123 1124 Not Collective 1125 1126 Input Parameter: 1127 . snes - SNES context 1128 1129 Output Parameter: 1130 . maxFails - maximum of unsuccessful steps 1131 1132 Level: intermediate 1133 1134 .keywords: SNES, nonlinear, get, maximum, unsuccessful, steps 1135 1136 .seealso: SNESGetMaxLinearSolveFailures(), SNESGetLinearSolveIterations(), SNESSetMaxLinearSolveFailures(), SNESGetLinearSolveFailures(), 1137 SNESSetMaxNonlinearStepFailures(), SNESGetNonlinearStepFailures() 1138 1139 @*/ 1140 PetscErrorCode SNESGetMaxNonlinearStepFailures(SNES snes, PetscInt *maxFails) 1141 { 1142 PetscFunctionBegin; 1143 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1144 PetscValidIntPointer(maxFails,2); 1145 *maxFails = snes->maxFailures; 1146 PetscFunctionReturn(0); 1147 } 1148 1149 #undef __FUNCT__ 1150 #define __FUNCT__ "SNESGetNumberFunctionEvals" 1151 /*@ 1152 SNESGetNumberFunctionEvals - Gets the number of user provided function evaluations 1153 done by SNES. 1154 1155 Not Collective 1156 1157 Input Parameter: 1158 . snes - SNES context 1159 1160 Output Parameter: 1161 . nfuncs - number of evaluations 1162 1163 Level: intermediate 1164 1165 .keywords: SNES, nonlinear, get, maximum, unsuccessful, steps 1166 1167 .seealso: SNESGetMaxLinearSolveFailures(), SNESGetLinearSolveIterations(), SNESSetMaxLinearSolveFailures(), SNESGetLinearSolveFailures() 1168 @*/ 1169 PetscErrorCode SNESGetNumberFunctionEvals(SNES snes, PetscInt *nfuncs) 1170 { 1171 PetscFunctionBegin; 1172 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1173 PetscValidIntPointer(nfuncs,2); 1174 *nfuncs = snes->nfuncs; 1175 PetscFunctionReturn(0); 1176 } 1177 1178 #undef __FUNCT__ 1179 #define __FUNCT__ "SNESGetLinearSolveFailures" 1180 /*@ 1181 SNESGetLinearSolveFailures - Gets the number of failed (non-converged) 1182 linear solvers. 1183 1184 Not Collective 1185 1186 Input Parameter: 1187 . snes - SNES context 1188 1189 Output Parameter: 1190 . nfails - number of failed solves 1191 1192 Notes: 1193 This counter is reset to zero for each successive call to SNESSolve(). 1194 1195 Level: intermediate 1196 1197 .keywords: SNES, nonlinear, get, number, unsuccessful, steps 1198 1199 .seealso: SNESGetMaxLinearSolveFailures(), SNESGetLinearSolveIterations(), SNESSetMaxLinearSolveFailures() 1200 @*/ 1201 PetscErrorCode SNESGetLinearSolveFailures(SNES snes,PetscInt *nfails) 1202 { 1203 PetscFunctionBegin; 1204 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1205 PetscValidIntPointer(nfails,2); 1206 *nfails = snes->numLinearSolveFailures; 1207 PetscFunctionReturn(0); 1208 } 1209 1210 #undef __FUNCT__ 1211 #define __FUNCT__ "SNESSetMaxLinearSolveFailures" 1212 /*@ 1213 SNESSetMaxLinearSolveFailures - the number of failed linear solve attempts 1214 allowed before SNES returns with a diverged reason of SNES_DIVERGED_LINEAR_SOLVE 1215 1216 Logically Collective on SNES 1217 1218 Input Parameters: 1219 + snes - SNES context 1220 - maxFails - maximum allowed linear solve failures 1221 1222 Level: intermediate 1223 1224 Notes: By default this is 0; that is SNES returns on the first failed linear solve 1225 1226 .keywords: SNES, nonlinear, set, maximum, unsuccessful, steps 1227 1228 .seealso: SNESGetLinearSolveFailures(), SNESGetMaxLinearSolveFailures(), SNESGetLinearSolveIterations() 1229 @*/ 1230 PetscErrorCode SNESSetMaxLinearSolveFailures(SNES snes, PetscInt maxFails) 1231 { 1232 PetscFunctionBegin; 1233 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1234 PetscValidLogicalCollectiveInt(snes,maxFails,2); 1235 snes->maxLinearSolveFailures = maxFails; 1236 PetscFunctionReturn(0); 1237 } 1238 1239 #undef __FUNCT__ 1240 #define __FUNCT__ "SNESGetMaxLinearSolveFailures" 1241 /*@ 1242 SNESGetMaxLinearSolveFailures - gets the maximum number of linear solve failures that 1243 are allowed before SNES terminates 1244 1245 Not Collective 1246 1247 Input Parameter: 1248 . snes - SNES context 1249 1250 Output Parameter: 1251 . maxFails - maximum of unsuccessful solves allowed 1252 1253 Level: intermediate 1254 1255 Notes: By default this is 1; that is SNES returns on the first failed linear solve 1256 1257 .keywords: SNES, nonlinear, get, maximum, unsuccessful, steps 1258 1259 .seealso: SNESGetLinearSolveFailures(), SNESGetLinearSolveIterations(), SNESSetMaxLinearSolveFailures(), 1260 @*/ 1261 PetscErrorCode SNESGetMaxLinearSolveFailures(SNES snes, PetscInt *maxFails) 1262 { 1263 PetscFunctionBegin; 1264 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1265 PetscValidIntPointer(maxFails,2); 1266 *maxFails = snes->maxLinearSolveFailures; 1267 PetscFunctionReturn(0); 1268 } 1269 1270 #undef __FUNCT__ 1271 #define __FUNCT__ "SNESGetLinearSolveIterations" 1272 /*@ 1273 SNESGetLinearSolveIterations - Gets the total number of linear iterations 1274 used by the nonlinear solver. 1275 1276 Not Collective 1277 1278 Input Parameter: 1279 . snes - SNES context 1280 1281 Output Parameter: 1282 . lits - number of linear iterations 1283 1284 Notes: 1285 This counter is reset to zero for each successive call to SNESSolve(). 1286 1287 Level: intermediate 1288 1289 .keywords: SNES, nonlinear, get, number, linear, iterations 1290 1291 .seealso: SNESGetIterationNumber(), SNESGetFunctionNorm(), SNESGetLinearSolveFailures(), SNESGetMaxLinearSolveFailures() 1292 @*/ 1293 PetscErrorCode SNESGetLinearSolveIterations(SNES snes,PetscInt *lits) 1294 { 1295 PetscFunctionBegin; 1296 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1297 PetscValidIntPointer(lits,2); 1298 *lits = snes->linear_its; 1299 PetscFunctionReturn(0); 1300 } 1301 1302 #undef __FUNCT__ 1303 #define __FUNCT__ "SNESGetKSP" 1304 /*@ 1305 SNESGetKSP - Returns the KSP context for a SNES solver. 1306 1307 Not Collective, but if SNES object is parallel, then KSP object is parallel 1308 1309 Input Parameter: 1310 . snes - the SNES context 1311 1312 Output Parameter: 1313 . ksp - the KSP context 1314 1315 Notes: 1316 The user can then directly manipulate the KSP context to set various 1317 options, etc. Likewise, the user can then extract and manipulate the 1318 PC contexts as well. 1319 1320 Level: beginner 1321 1322 .keywords: SNES, nonlinear, get, KSP, context 1323 1324 .seealso: KSPGetPC(), SNESCreate(), KSPCreate(), SNESSetKSP() 1325 @*/ 1326 PetscErrorCode SNESGetKSP(SNES snes,KSP *ksp) 1327 { 1328 PetscErrorCode ierr; 1329 1330 PetscFunctionBegin; 1331 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1332 PetscValidPointer(ksp,2); 1333 1334 if (!snes->ksp) { 1335 ierr = KSPCreate(((PetscObject)snes)->comm,&snes->ksp);CHKERRQ(ierr); 1336 ierr = PetscObjectIncrementTabLevel((PetscObject)snes->ksp,(PetscObject)snes,1);CHKERRQ(ierr); 1337 ierr = PetscLogObjectParent(snes,snes->ksp);CHKERRQ(ierr); 1338 } 1339 *ksp = snes->ksp; 1340 PetscFunctionReturn(0); 1341 } 1342 1343 #undef __FUNCT__ 1344 #define __FUNCT__ "SNESSetKSP" 1345 /*@ 1346 SNESSetKSP - Sets a KSP context for the SNES object to use 1347 1348 Not Collective, but the SNES and KSP objects must live on the same MPI_Comm 1349 1350 Input Parameters: 1351 + snes - the SNES context 1352 - ksp - the KSP context 1353 1354 Notes: 1355 The SNES object already has its KSP object, you can obtain with SNESGetKSP() 1356 so this routine is rarely needed. 1357 1358 The KSP object that is already in the SNES object has its reference count 1359 decreased by one. 1360 1361 Level: developer 1362 1363 .keywords: SNES, nonlinear, get, KSP, context 1364 1365 .seealso: KSPGetPC(), SNESCreate(), KSPCreate(), SNESSetKSP() 1366 @*/ 1367 PetscErrorCode SNESSetKSP(SNES snes,KSP ksp) 1368 { 1369 PetscErrorCode ierr; 1370 1371 PetscFunctionBegin; 1372 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1373 PetscValidHeaderSpecific(ksp,KSP_CLASSID,2); 1374 PetscCheckSameComm(snes,1,ksp,2); 1375 ierr = PetscObjectReference((PetscObject)ksp);CHKERRQ(ierr); 1376 if (snes->ksp) {ierr = PetscObjectDereference((PetscObject)snes->ksp);CHKERRQ(ierr);} 1377 snes->ksp = ksp; 1378 PetscFunctionReturn(0); 1379 } 1380 1381 #if 0 1382 #undef __FUNCT__ 1383 #define __FUNCT__ "SNESPublish_Petsc" 1384 static PetscErrorCode SNESPublish_Petsc(PetscObject obj) 1385 { 1386 PetscFunctionBegin; 1387 PetscFunctionReturn(0); 1388 } 1389 #endif 1390 1391 /* -----------------------------------------------------------*/ 1392 #undef __FUNCT__ 1393 #define __FUNCT__ "SNESCreate" 1394 /*@ 1395 SNESCreate - Creates a nonlinear solver context. 1396 1397 Collective on MPI_Comm 1398 1399 Input Parameters: 1400 . comm - MPI communicator 1401 1402 Output Parameter: 1403 . outsnes - the new SNES context 1404 1405 Options Database Keys: 1406 + -snes_mf - Activates default matrix-free Jacobian-vector products, 1407 and no preconditioning matrix 1408 . -snes_mf_operator - Activates default matrix-free Jacobian-vector 1409 products, and a user-provided preconditioning matrix 1410 as set by SNESSetJacobian() 1411 - -snes_fd - Uses (slow!) finite differences to compute Jacobian 1412 1413 Level: beginner 1414 1415 .keywords: SNES, nonlinear, create, context 1416 1417 .seealso: SNESSolve(), SNESDestroy(), SNES, SNESSetLagPreconditioner() 1418 1419 @*/ 1420 PetscErrorCode SNESCreate(MPI_Comm comm,SNES *outsnes) 1421 { 1422 PetscErrorCode ierr; 1423 SNES snes; 1424 SNESKSPEW *kctx; 1425 1426 PetscFunctionBegin; 1427 PetscValidPointer(outsnes,2); 1428 *outsnes = NULL; 1429 #if !defined(PETSC_USE_DYNAMIC_LIBRARIES) 1430 ierr = SNESInitializePackage(NULL);CHKERRQ(ierr); 1431 #endif 1432 1433 ierr = PetscHeaderCreate(snes,_p_SNES,struct _SNESOps,SNES_CLASSID,0,"SNES","Nonlinear solver","SNES",comm,SNESDestroy,SNESView);CHKERRQ(ierr); 1434 1435 snes->ops->converged = SNESDefaultConverged; 1436 snes->usesksp = PETSC_TRUE; 1437 snes->tolerancesset = PETSC_FALSE; 1438 snes->max_its = 50; 1439 snes->max_funcs = 10000; 1440 snes->norm = 0.0; 1441 snes->normtype = SNES_NORM_FUNCTION; 1442 snes->rtol = 1.e-8; 1443 snes->ttol = 0.0; 1444 snes->abstol = 1.e-50; 1445 snes->stol = 1.e-8; 1446 snes->deltatol = 1.e-12; 1447 snes->nfuncs = 0; 1448 snes->numFailures = 0; 1449 snes->maxFailures = 1; 1450 snes->linear_its = 0; 1451 snes->lagjacobian = 1; 1452 snes->lagpreconditioner = 1; 1453 snes->numbermonitors = 0; 1454 snes->data = 0; 1455 snes->setupcalled = PETSC_FALSE; 1456 snes->ksp_ewconv = PETSC_FALSE; 1457 snes->nwork = 0; 1458 snes->work = 0; 1459 snes->nvwork = 0; 1460 snes->vwork = 0; 1461 snes->conv_hist_len = 0; 1462 snes->conv_hist_max = 0; 1463 snes->conv_hist = NULL; 1464 snes->conv_hist_its = NULL; 1465 snes->conv_hist_reset = PETSC_TRUE; 1466 snes->vec_func_init_set = PETSC_FALSE; 1467 snes->norm_init = 0.; 1468 snes->norm_init_set = PETSC_FALSE; 1469 snes->reason = SNES_CONVERGED_ITERATING; 1470 snes->pcside = PC_RIGHT; 1471 1472 snes->mf = PETSC_FALSE; 1473 snes->mf_operator = PETSC_FALSE; 1474 snes->mf_version = 1; 1475 1476 snes->numLinearSolveFailures = 0; 1477 snes->maxLinearSolveFailures = 1; 1478 1479 /* Create context to compute Eisenstat-Walker relative tolerance for KSP */ 1480 ierr = PetscNewLog(snes,SNESKSPEW,&kctx);CHKERRQ(ierr); 1481 1482 snes->kspconvctx = (void*)kctx; 1483 kctx->version = 2; 1484 kctx->rtol_0 = .3; /* Eisenstat and Walker suggest rtol_0=.5, but 1485 this was too large for some test cases */ 1486 kctx->rtol_last = 0.0; 1487 kctx->rtol_max = .9; 1488 kctx->gamma = 1.0; 1489 kctx->alpha = .5*(1.0 + PetscSqrtReal(5.0)); 1490 kctx->alpha2 = kctx->alpha; 1491 kctx->threshold = .1; 1492 kctx->lresid_last = 0.0; 1493 kctx->norm_last = 0.0; 1494 1495 *outsnes = snes; 1496 PetscFunctionReturn(0); 1497 } 1498 1499 /*MC 1500 SNESFunction - functional form used to convey the nonlinear function to be solved by SNES 1501 1502 Synopsis: 1503 #include "petscsnes.h" 1504 SNESFunction(SNES snes,Vec x,Vec f,void *ctx); 1505 1506 Input Parameters: 1507 + snes - the SNES context 1508 . x - state at which to evaluate residual 1509 - ctx - optional user-defined function context, passed in with SNESSetFunction() 1510 1511 Output Parameter: 1512 . f - vector to put residual (function value) 1513 1514 .seealso: SNESSetFunction(), SNESGetFunction() 1515 M*/ 1516 1517 #undef __FUNCT__ 1518 #define __FUNCT__ "SNESSetFunction" 1519 /*@C 1520 SNESSetFunction - Sets the function evaluation routine and function 1521 vector for use by the SNES routines in solving systems of nonlinear 1522 equations. 1523 1524 Logically Collective on SNES 1525 1526 Input Parameters: 1527 + snes - the SNES context 1528 . r - vector to store function value 1529 . SNESFunction - function evaluation routine 1530 - ctx - [optional] user-defined context for private data for the 1531 function evaluation routine (may be NULL) 1532 1533 Notes: 1534 The Newton-like methods typically solve linear systems of the form 1535 $ f'(x) x = -f(x), 1536 where f'(x) denotes the Jacobian matrix and f(x) is the function. 1537 1538 Level: beginner 1539 1540 .keywords: SNES, nonlinear, set, function 1541 1542 .seealso: SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESSetPicard(), SNESFunction 1543 @*/ 1544 PetscErrorCode SNESSetFunction(SNES snes,Vec r,PetscErrorCode (*SNESFunction)(SNES,Vec,Vec,void*),void *ctx) 1545 { 1546 PetscErrorCode ierr; 1547 DM dm; 1548 1549 PetscFunctionBegin; 1550 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1551 if (r) { 1552 PetscValidHeaderSpecific(r,VEC_CLASSID,2); 1553 PetscCheckSameComm(snes,1,r,2); 1554 ierr = PetscObjectReference((PetscObject)r);CHKERRQ(ierr); 1555 ierr = VecDestroy(&snes->vec_func);CHKERRQ(ierr); 1556 1557 snes->vec_func = r; 1558 } 1559 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 1560 ierr = DMSNESSetFunction(dm,SNESFunction,ctx);CHKERRQ(ierr); 1561 PetscFunctionReturn(0); 1562 } 1563 1564 1565 #undef __FUNCT__ 1566 #define __FUNCT__ "SNESSetInitialFunction" 1567 /*@C 1568 SNESSetInitialFunction - Sets the function vector to be used as the 1569 function norm at the initialization of the method. In some 1570 instances, the user has precomputed the function before calling 1571 SNESSolve. This function allows one to avoid a redundant call 1572 to SNESComputeFunction in that case. 1573 1574 Logically Collective on SNES 1575 1576 Input Parameters: 1577 + snes - the SNES context 1578 - f - vector to store function value 1579 1580 Notes: 1581 This should not be modified during the solution procedure. 1582 1583 This is used extensively in the SNESFAS hierarchy and in nonlinear preconditioning. 1584 1585 Level: developer 1586 1587 .keywords: SNES, nonlinear, set, function 1588 1589 .seealso: SNESSetFunction(), SNESComputeFunction(), SNESSetInitialFunctionNorm() 1590 @*/ 1591 PetscErrorCode SNESSetInitialFunction(SNES snes, Vec f) 1592 { 1593 PetscErrorCode ierr; 1594 Vec vec_func; 1595 1596 PetscFunctionBegin; 1597 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1598 PetscValidHeaderSpecific(f,VEC_CLASSID,2); 1599 PetscCheckSameComm(snes,1,f,2); 1600 ierr = SNESGetFunction(snes,&vec_func,NULL,NULL);CHKERRQ(ierr); 1601 ierr = VecCopy(f, vec_func);CHKERRQ(ierr); 1602 1603 snes->vec_func_init_set = PETSC_TRUE; 1604 PetscFunctionReturn(0); 1605 } 1606 1607 1608 #undef __FUNCT__ 1609 #define __FUNCT__ "SNESSetInitialFunctionNorm" 1610 /*@C 1611 SNESSetInitialFunctionNorm - Sets the function norm to be used as the function norm 1612 at the initialization of the method. In some instances, the user has precomputed 1613 the function and its norm before calling SNESSolve. This function allows one to 1614 avoid a redundant call to SNESComputeFunction() and VecNorm() in that case. 1615 1616 Logically Collective on SNES 1617 1618 Input Parameters: 1619 + snes - the SNES context 1620 - fnorm - the norm of F as set by SNESSetInitialFunction() 1621 1622 This is used extensively in the SNESFAS hierarchy and in nonlinear preconditioning. 1623 1624 Level: developer 1625 1626 .keywords: SNES, nonlinear, set, function, norm 1627 1628 .seealso: SNESSetFunction(), SNESComputeFunction(), SNESSetInitialFunction() 1629 @*/ 1630 PetscErrorCode SNESSetInitialFunctionNorm(SNES snes, PetscReal fnorm) 1631 { 1632 PetscFunctionBegin; 1633 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1634 snes->norm_init = fnorm; 1635 snes->norm_init_set = PETSC_TRUE; 1636 PetscFunctionReturn(0); 1637 } 1638 1639 #undef __FUNCT__ 1640 #define __FUNCT__ "SNESSetNormType" 1641 /*@ 1642 SNESSetNormType - Sets the SNESNormType used in covergence and monitoring 1643 of the SNES method. 1644 1645 Logically Collective on SNES 1646 1647 Input Parameters: 1648 + snes - the SNES context 1649 - normtype - the type of the norm used 1650 1651 Notes: 1652 Only certain SNES methods support certain SNESNormTypes. Most require evaluation 1653 of the nonlinear function and the taking of its norm at every iteration to 1654 even ensure convergence at all. However, methods such as custom Gauss-Seidel methods 1655 (SNESGS) and the like do not require the norm of the function to be computed, and therfore 1656 may either be monitored for convergence or not. As these are often used as nonlinear 1657 preconditioners, monitoring the norm of their error is not a useful enterprise within 1658 their solution. 1659 1660 Level: developer 1661 1662 .keywords: SNES, nonlinear, set, function, norm, type 1663 1664 .seealso: SNESGetNormType(), SNESComputeFunction(), VecNorm(), SNESSetFunction(), SNESSetInitialFunction(), SNESNormType 1665 @*/ 1666 PetscErrorCode SNESSetNormType(SNES snes, SNESNormType normtype) 1667 { 1668 PetscFunctionBegin; 1669 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1670 snes->normtype = normtype; 1671 PetscFunctionReturn(0); 1672 } 1673 1674 1675 #undef __FUNCT__ 1676 #define __FUNCT__ "SNESGetNormType" 1677 /*@ 1678 SNESGetNormType - Gets the SNESNormType used in covergence and monitoring 1679 of the SNES method. 1680 1681 Logically Collective on SNES 1682 1683 Input Parameters: 1684 + snes - the SNES context 1685 - normtype - the type of the norm used 1686 1687 Level: advanced 1688 1689 .keywords: SNES, nonlinear, set, function, norm, type 1690 1691 .seealso: SNESSetNormType(), SNESComputeFunction(), VecNorm(), SNESSetFunction(), SNESSetInitialFunction(), SNESNormType 1692 @*/ 1693 PetscErrorCode SNESGetNormType(SNES snes, SNESNormType *normtype) 1694 { 1695 PetscFunctionBegin; 1696 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1697 *normtype = snes->normtype; 1698 PetscFunctionReturn(0); 1699 } 1700 1701 /*MC 1702 SNESGSFunction - function used to convey a Gauss-Seidel sweep on the nonlinear function 1703 1704 Synopsis: 1705 #include "petscsnes.h" 1706 $ SNESGSFunction(SNES snes,Vec x,Vec b,void *ctx); 1707 1708 + X - solution vector 1709 . B - RHS vector 1710 - ctx - optional user-defined Gauss-Seidel context 1711 1712 .seealso: SNESSetGS(), SNESGetGS() 1713 M*/ 1714 1715 #undef __FUNCT__ 1716 #define __FUNCT__ "SNESSetGS" 1717 /*@C 1718 SNESSetGS - Sets the user nonlinear Gauss-Seidel routine for 1719 use with composed nonlinear solvers. 1720 1721 Input Parameters: 1722 + snes - the SNES context 1723 . SNESGSFunction - function evaluation routine 1724 - ctx - [optional] user-defined context for private data for the 1725 smoother evaluation routine (may be NULL) 1726 1727 Notes: 1728 The GS routines are used by the composed nonlinear solver to generate 1729 a problem appropriate update to the solution, particularly FAS. 1730 1731 Level: intermediate 1732 1733 .keywords: SNES, nonlinear, set, Gauss-Seidel 1734 1735 .seealso: SNESGetFunction(), SNESComputeGS() 1736 @*/ 1737 PetscErrorCode SNESSetGS(SNES snes,PetscErrorCode (*SNESGSFunction)(SNES,Vec,Vec,void*),void *ctx) 1738 { 1739 PetscErrorCode ierr; 1740 DM dm; 1741 1742 PetscFunctionBegin; 1743 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1744 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 1745 ierr = DMSNESSetGS(dm,SNESGSFunction,ctx);CHKERRQ(ierr); 1746 PetscFunctionReturn(0); 1747 } 1748 1749 #undef __FUNCT__ 1750 #define __FUNCT__ "SNESPicardComputeFunction" 1751 PetscErrorCode SNESPicardComputeFunction(SNES snes,Vec x,Vec f,void *ctx) 1752 { 1753 PetscErrorCode ierr; 1754 DM dm; 1755 DMSNES sdm; 1756 1757 PetscFunctionBegin; 1758 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 1759 ierr = DMGetDMSNES(dm,&sdm);CHKERRQ(ierr); 1760 /* A(x)*x - b(x) */ 1761 if (sdm->ops->computepfunction) { 1762 ierr = (*sdm->ops->computepfunction)(snes,x,f,sdm->pctx);CHKERRQ(ierr); 1763 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE, "Must call SNESSetPicard() to provide Picard function."); 1764 1765 if (sdm->ops->computepjacobian) { 1766 ierr = (*sdm->ops->computepjacobian)(snes,x,&snes->jacobian,&snes->jacobian_pre,&snes->matstruct,sdm->pctx);CHKERRQ(ierr); 1767 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE, "Must call SNESSetPicard() to provide Picard matrix."); 1768 ierr = VecScale(f,-1.0);CHKERRQ(ierr); 1769 ierr = MatMultAdd(snes->jacobian,x,f,f);CHKERRQ(ierr); 1770 PetscFunctionReturn(0); 1771 } 1772 1773 #undef __FUNCT__ 1774 #define __FUNCT__ "SNESPicardComputeJacobian" 1775 PetscErrorCode SNESPicardComputeJacobian(SNES snes,Vec x1,Mat *J,Mat *B,MatStructure *flag,void *ctx) 1776 { 1777 PetscFunctionBegin; 1778 /* the jacobian matrix should be pre-filled in SNESPicardComputeFunction */ 1779 *flag = snes->matstruct; 1780 PetscFunctionReturn(0); 1781 } 1782 1783 #undef __FUNCT__ 1784 #define __FUNCT__ "SNESSetPicard" 1785 /*@C 1786 SNESSetPicard - Use SNES to solve the semilinear-system A(x) x = b(x) via a Picard type iteration (Picard linearization) 1787 1788 Logically Collective on SNES 1789 1790 Input Parameters: 1791 + snes - the SNES context 1792 . r - vector to store function value 1793 . SNESFunction - function evaluation routine 1794 . jmat - normally the same as mat but you can pass another matrix for which you compute A(x) x - b(x) (see jmat below) 1795 . mat - matrix to store A 1796 . SNESJacobianFunction - function to compute matrix value 1797 - ctx - [optional] user-defined context for private data for the 1798 function evaluation routine (may be NULL) 1799 1800 Notes: 1801 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 1802 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. 1803 1804 One can call SNESSetPicard() or SNESSetFunction() (and possibly SNESSetJacobian()) but cannot call both 1805 1806 $ 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} 1807 $ Note that when an exact solver is used this corresponds to the "classic" Picard A(x^{n}) x^{n+1} = b(x^{n}) iteration. 1808 1809 Run with -snes_mf_operator to solve the system with Newton's method using A(x^{n}) to construct the preconditioner. 1810 1811 We implement the defect correction form of the Picard iteration because it converges much more generally when inexact linear solvers are used then 1812 the direct Picard iteration A(x^n) x^{n+1} = b(x^n) 1813 1814 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 1815 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 1816 different please contact us at petsc-dev@mcs.anl.gov and we'll have an entirely new argument :-). 1817 1818 Level: intermediate 1819 1820 .keywords: SNES, nonlinear, set, function 1821 1822 .seealso: SNESGetFunction(), SNESSetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESGetPicard(), SNESLineSearchPreCheckPicard() 1823 @*/ 1824 PetscErrorCode SNESSetPicard(SNES snes,Vec r,PetscErrorCode (*SNESFunction)(SNES,Vec,Vec,void*),Mat jmat, Mat mat, PetscErrorCode (*SNESJacobianFunction)(SNES,Vec,Mat*,Mat*,MatStructure*,void*),void *ctx) 1825 { 1826 PetscErrorCode ierr; 1827 DM dm; 1828 1829 PetscFunctionBegin; 1830 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1831 ierr = SNESGetDM(snes, &dm);CHKERRQ(ierr); 1832 ierr = DMSNESSetPicard(dm,SNESFunction,SNESJacobianFunction,ctx);CHKERRQ(ierr); 1833 ierr = SNESSetFunction(snes,r,SNESPicardComputeFunction,ctx);CHKERRQ(ierr); 1834 ierr = SNESSetJacobian(snes,jmat,mat,SNESPicardComputeJacobian,ctx);CHKERRQ(ierr); 1835 PetscFunctionReturn(0); 1836 } 1837 1838 #undef __FUNCT__ 1839 #define __FUNCT__ "SNESGetPicard" 1840 /*@C 1841 SNESGetPicard - Returns the context for the Picard iteration 1842 1843 Not Collective, but Vec is parallel if SNES is parallel. Collective if Vec is requested, but has not been created yet. 1844 1845 Input Parameter: 1846 . snes - the SNES context 1847 1848 Output Parameter: 1849 + r - the function (or NULL) 1850 . SNESFunction - the function (or NULL) 1851 . jmat - the picard matrix (or NULL) 1852 . mat - the picard preconditioner matrix (or NULL) 1853 . SNESJacobianFunction - the function for matrix evaluation (or NULL) 1854 - ctx - the function context (or NULL) 1855 1856 Level: advanced 1857 1858 .keywords: SNES, nonlinear, get, function 1859 1860 .seealso: SNESSetPicard, SNESGetFunction, SNESGetJacobian, SNESGetDM, SNESFunction, SNESJacobianFunction 1861 @*/ 1862 PetscErrorCode SNESGetPicard(SNES snes,Vec *r,PetscErrorCode (**SNESFunction)(SNES,Vec,Vec,void*),Mat *jmat, Mat *mat, PetscErrorCode (**SNESJacobianFunction)(SNES,Vec,Mat*,Mat*,MatStructure*,void*),void **ctx) 1863 { 1864 PetscErrorCode ierr; 1865 DM dm; 1866 1867 PetscFunctionBegin; 1868 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1869 ierr = SNESGetFunction(snes,r,NULL,NULL);CHKERRQ(ierr); 1870 ierr = SNESGetJacobian(snes,jmat,mat,NULL,NULL);CHKERRQ(ierr); 1871 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 1872 ierr = DMSNESGetPicard(dm,SNESFunction,SNESJacobianFunction,ctx);CHKERRQ(ierr); 1873 PetscFunctionReturn(0); 1874 } 1875 1876 #undef __FUNCT__ 1877 #define __FUNCT__ "SNESSetComputeInitialGuess" 1878 /*@C 1879 SNESSetComputeInitialGuess - Sets a routine used to compute an initial guess for the problem 1880 1881 Logically Collective on SNES 1882 1883 Input Parameters: 1884 + snes - the SNES context 1885 . func - function evaluation routine 1886 - ctx - [optional] user-defined context for private data for the 1887 function evaluation routine (may be NULL) 1888 1889 Calling sequence of func: 1890 $ func (SNES snes,Vec x,void *ctx); 1891 1892 . f - function vector 1893 - ctx - optional user-defined function context 1894 1895 Level: intermediate 1896 1897 .keywords: SNES, nonlinear, set, function 1898 1899 .seealso: SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian() 1900 @*/ 1901 PetscErrorCode SNESSetComputeInitialGuess(SNES snes,PetscErrorCode (*func)(SNES,Vec,void*),void *ctx) 1902 { 1903 PetscFunctionBegin; 1904 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1905 if (func) snes->ops->computeinitialguess = func; 1906 if (ctx) snes->initialguessP = ctx; 1907 PetscFunctionReturn(0); 1908 } 1909 1910 /* --------------------------------------------------------------- */ 1911 #undef __FUNCT__ 1912 #define __FUNCT__ "SNESGetRhs" 1913 /*@C 1914 SNESGetRhs - Gets the vector for solving F(x) = rhs. If rhs is not set 1915 it assumes a zero right hand side. 1916 1917 Logically Collective on SNES 1918 1919 Input Parameter: 1920 . snes - the SNES context 1921 1922 Output Parameter: 1923 . rhs - the right hand side vector or NULL if the right hand side vector is null 1924 1925 Level: intermediate 1926 1927 .keywords: SNES, nonlinear, get, function, right hand side 1928 1929 .seealso: SNESGetSolution(), SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESSetFunction() 1930 @*/ 1931 PetscErrorCode SNESGetRhs(SNES snes,Vec *rhs) 1932 { 1933 PetscFunctionBegin; 1934 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1935 PetscValidPointer(rhs,2); 1936 *rhs = snes->vec_rhs; 1937 PetscFunctionReturn(0); 1938 } 1939 1940 #undef __FUNCT__ 1941 #define __FUNCT__ "SNESComputeFunction" 1942 /*@ 1943 SNESComputeFunction - Calls the function that has been set with SNESSetFunction(). 1944 1945 Collective on SNES 1946 1947 Input Parameters: 1948 + snes - the SNES context 1949 - x - input vector 1950 1951 Output Parameter: 1952 . y - function vector, as set by SNESSetFunction() 1953 1954 Notes: 1955 SNESComputeFunction() is typically used within nonlinear solvers 1956 implementations, so most users would not generally call this routine 1957 themselves. 1958 1959 Level: developer 1960 1961 .keywords: SNES, nonlinear, compute, function 1962 1963 .seealso: SNESSetFunction(), SNESGetFunction() 1964 @*/ 1965 PetscErrorCode SNESComputeFunction(SNES snes,Vec x,Vec y) 1966 { 1967 PetscErrorCode ierr; 1968 DM dm; 1969 DMSNES sdm; 1970 1971 PetscFunctionBegin; 1972 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1973 PetscValidHeaderSpecific(x,VEC_CLASSID,2); 1974 PetscValidHeaderSpecific(y,VEC_CLASSID,3); 1975 PetscCheckSameComm(snes,1,x,2); 1976 PetscCheckSameComm(snes,1,y,3); 1977 VecValidValues(x,2,PETSC_TRUE); 1978 1979 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 1980 ierr = DMGetDMSNES(dm,&sdm);CHKERRQ(ierr); 1981 if (snes->pc && snes->pcside == PC_LEFT) { 1982 ierr = VecCopy(x,y);CHKERRQ(ierr); 1983 ierr = PetscLogEventBegin(SNES_NPCSolve,snes->pc,x,y,0);CHKERRQ(ierr); 1984 ierr = SNESSolve(snes->pc,snes->vec_rhs,y);CHKERRQ(ierr); 1985 ierr = PetscLogEventEnd(SNES_NPCSolve,snes->pc,x,y,0);CHKERRQ(ierr); 1986 ierr = VecAYPX(y,-1.0,x);CHKERRQ(ierr); 1987 } else if (sdm->ops->computefunction) { 1988 ierr = PetscLogEventBegin(SNES_FunctionEval,snes,x,y,0);CHKERRQ(ierr); 1989 PetscStackPush("SNES user function"); 1990 ierr = (*sdm->ops->computefunction)(snes,x,y,sdm->functionctx);CHKERRQ(ierr); 1991 PetscStackPop; 1992 ierr = PetscLogEventEnd(SNES_FunctionEval,snes,x,y,0);CHKERRQ(ierr); 1993 } else if (snes->vec_rhs) { 1994 ierr = MatMult(snes->jacobian, x, y);CHKERRQ(ierr); 1995 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE, "Must call SNESSetFunction() or SNESSetDM() before SNESComputeFunction(), likely called from SNESSolve()."); 1996 if (snes->vec_rhs) { 1997 ierr = VecAXPY(y,-1.0,snes->vec_rhs);CHKERRQ(ierr); 1998 } 1999 snes->nfuncs++; 2000 VecValidValues(y,3,PETSC_FALSE); 2001 PetscFunctionReturn(0); 2002 } 2003 2004 #undef __FUNCT__ 2005 #define __FUNCT__ "SNESComputeGS" 2006 /*@ 2007 SNESComputeGS - Calls the Gauss-Seidel function that has been set with SNESSetGS(). 2008 2009 Collective on SNES 2010 2011 Input Parameters: 2012 + snes - the SNES context 2013 . x - input vector 2014 - b - rhs vector 2015 2016 Output Parameter: 2017 . x - new solution vector 2018 2019 Notes: 2020 SNESComputeGS() is typically used within composed nonlinear solver 2021 implementations, so most users would not generally call this routine 2022 themselves. 2023 2024 Level: developer 2025 2026 .keywords: SNES, nonlinear, compute, function 2027 2028 .seealso: SNESSetGS(), SNESComputeFunction() 2029 @*/ 2030 PetscErrorCode SNESComputeGS(SNES snes,Vec b,Vec x) 2031 { 2032 PetscErrorCode ierr; 2033 DM dm; 2034 DMSNES sdm; 2035 2036 PetscFunctionBegin; 2037 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2038 PetscValidHeaderSpecific(x,VEC_CLASSID,2); 2039 if (b) PetscValidHeaderSpecific(b,VEC_CLASSID,3); 2040 PetscCheckSameComm(snes,1,x,2); 2041 if (b) PetscCheckSameComm(snes,1,b,3); 2042 if (b) VecValidValues(b,2,PETSC_TRUE); 2043 ierr = PetscLogEventBegin(SNES_GSEval,snes,x,b,0);CHKERRQ(ierr); 2044 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 2045 ierr = DMGetDMSNES(dm,&sdm);CHKERRQ(ierr); 2046 if (sdm->ops->computegs) { 2047 PetscStackPush("SNES user GS"); 2048 ierr = (*sdm->ops->computegs)(snes,x,b,sdm->gsctx);CHKERRQ(ierr); 2049 PetscStackPop; 2050 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE, "Must call SNESSetGS() before SNESComputeGS(), likely called from SNESSolve()."); 2051 ierr = PetscLogEventEnd(SNES_GSEval,snes,x,b,0);CHKERRQ(ierr); 2052 VecValidValues(x,3,PETSC_FALSE); 2053 PetscFunctionReturn(0); 2054 } 2055 2056 #undef __FUNCT__ 2057 #define __FUNCT__ "SNESComputeJacobian" 2058 /*@ 2059 SNESComputeJacobian - Computes the Jacobian matrix that has been set with SNESSetJacobian(). 2060 2061 Collective on SNES and Mat 2062 2063 Input Parameters: 2064 + snes - the SNES context 2065 - x - input vector 2066 2067 Output Parameters: 2068 + A - Jacobian matrix 2069 . B - optional preconditioning matrix 2070 - flag - flag indicating matrix structure (one of, SAME_NONZERO_PATTERN,DIFFERENT_NONZERO_PATTERN,SAME_PRECONDITIONER) 2071 2072 Options Database Keys: 2073 + -snes_lag_preconditioner <lag> 2074 . -snes_lag_jacobian <lag> 2075 . -snes_compare_explicit - Compare the computed Jacobian to the finite difference Jacobian and output the differences 2076 . -snes_compare_explicit_draw - Compare the computed Jacobian to the finite difference Jacobian and draw the result 2077 . -snes_compare_explicit_contour - Compare the computed Jacobian to the finite difference Jacobian and draw a contour plot with the result 2078 . -snes_compare_operator - Make the comparison options above use the operator instead of the preconditioning matrix 2079 . -snes_compare_coloring - Compute the finite differece Jacobian using coloring and display norms of difference 2080 . -snes_compare_coloring_display - Compute the finite differece Jacobian using coloring and display verbose differences 2081 . -snes_compare_coloring_threshold - Display only those matrix entries that differ by more than a given threshold 2082 . -snes_compare_coloring_threshold_atol - Absolute tolerance for difference in matrix entries to be displayed by -snes_compare_coloring_threshold 2083 . -snes_compare_coloring_threshold_rtol - Relative tolerance for difference in matrix entries to be displayed by -snes_compare_coloring_threshold 2084 . -snes_compare_coloring_draw - Compute the finite differece Jacobian using coloring and draw differences 2085 - -snes_compare_coloring_draw_contour - Compute the finite differece Jacobian using coloring and show contours of matrices and differences 2086 2087 2088 Notes: 2089 Most users should not need to explicitly call this routine, as it 2090 is used internally within the nonlinear solvers. 2091 2092 See KSPSetOperators() for important information about setting the 2093 flag parameter. 2094 2095 Level: developer 2096 2097 .keywords: SNES, compute, Jacobian, matrix 2098 2099 .seealso: SNESSetJacobian(), KSPSetOperators(), MatStructure, SNESSetLagPreconditioner(), SNESSetLagJacobian() 2100 @*/ 2101 PetscErrorCode SNESComputeJacobian(SNES snes,Vec X,Mat *A,Mat *B,MatStructure *flg) 2102 { 2103 PetscErrorCode ierr; 2104 PetscBool flag; 2105 DM dm; 2106 DMSNES sdm; 2107 2108 PetscFunctionBegin; 2109 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2110 PetscValidHeaderSpecific(X,VEC_CLASSID,2); 2111 PetscValidPointer(flg,5); 2112 PetscCheckSameComm(snes,1,X,2); 2113 VecValidValues(X,2,PETSC_TRUE); 2114 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 2115 ierr = DMGetDMSNES(dm,&sdm);CHKERRQ(ierr); 2116 2117 if (!sdm->ops->computejacobian) SETERRQ(((PetscObject)snes)->comm,PETSC_ERR_USER,"Must call SNESSetJacobian(), DMSNESSetJacobian(), DMDASNESSetJacobianLocal(), etc"); 2118 2119 /* make sure that MatAssemblyBegin/End() is called on A matrix if it is matrix free */ 2120 2121 if (snes->lagjacobian == -2) { 2122 snes->lagjacobian = -1; 2123 2124 ierr = PetscInfo(snes,"Recomputing Jacobian/preconditioner because lag is -2 (means compute Jacobian, but then never again) \n");CHKERRQ(ierr); 2125 } else if (snes->lagjacobian == -1) { 2126 *flg = SAME_PRECONDITIONER; 2127 ierr = PetscInfo(snes,"Reusing Jacobian/preconditioner because lag is -1\n");CHKERRQ(ierr); 2128 ierr = PetscObjectTypeCompare((PetscObject)*A,MATMFFD,&flag);CHKERRQ(ierr); 2129 if (flag) { 2130 ierr = MatAssemblyBegin(*A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2131 ierr = MatAssemblyEnd(*A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2132 } 2133 PetscFunctionReturn(0); 2134 } else if (snes->lagjacobian > 1 && snes->iter % snes->lagjacobian) { 2135 *flg = SAME_PRECONDITIONER; 2136 ierr = PetscInfo2(snes,"Reusing Jacobian/preconditioner because lag is %D and SNES iteration is %D\n",snes->lagjacobian,snes->iter);CHKERRQ(ierr); 2137 ierr = PetscObjectTypeCompare((PetscObject)*A,MATMFFD,&flag);CHKERRQ(ierr); 2138 if (flag) { 2139 ierr = MatAssemblyBegin(*A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2140 ierr = MatAssemblyEnd(*A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2141 } 2142 PetscFunctionReturn(0); 2143 } 2144 2145 *flg = DIFFERENT_NONZERO_PATTERN; 2146 ierr = PetscLogEventBegin(SNES_JacobianEval,snes,X,*A,*B);CHKERRQ(ierr); 2147 2148 PetscStackPush("SNES user Jacobian function"); 2149 ierr = (*sdm->ops->computejacobian)(snes,X,A,B,flg,sdm->jacobianctx);CHKERRQ(ierr); 2150 PetscStackPop; 2151 2152 ierr = PetscLogEventEnd(SNES_JacobianEval,snes,X,*A,*B);CHKERRQ(ierr); 2153 2154 if (snes->lagpreconditioner == -2) { 2155 ierr = PetscInfo(snes,"Rebuilding preconditioner exactly once since lag is -2\n");CHKERRQ(ierr); 2156 2157 snes->lagpreconditioner = -1; 2158 } else if (snes->lagpreconditioner == -1) { 2159 *flg = SAME_PRECONDITIONER; 2160 ierr = PetscInfo(snes,"Reusing preconditioner because lag is -1\n");CHKERRQ(ierr); 2161 } else if (snes->lagpreconditioner > 1 && snes->iter % snes->lagpreconditioner) { 2162 *flg = SAME_PRECONDITIONER; 2163 ierr = PetscInfo2(snes,"Reusing preconditioner because lag is %D and SNES iteration is %D\n",snes->lagpreconditioner,snes->iter);CHKERRQ(ierr); 2164 } 2165 2166 /* make sure user returned a correct Jacobian and preconditioner */ 2167 /* PetscValidHeaderSpecific(*A,MAT_CLASSID,3); 2168 PetscValidHeaderSpecific(*B,MAT_CLASSID,4); */ 2169 { 2170 PetscBool flag = PETSC_FALSE,flag_draw = PETSC_FALSE,flag_contour = PETSC_FALSE,flag_operator = PETSC_FALSE; 2171 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_explicit",&flag,NULL);CHKERRQ(ierr); 2172 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_explicit_draw",&flag_draw,NULL);CHKERRQ(ierr); 2173 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_explicit_draw_contour",&flag_contour,NULL);CHKERRQ(ierr); 2174 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_operator",&flag_operator,NULL);CHKERRQ(ierr); 2175 if (flag || flag_draw || flag_contour) { 2176 Mat Bexp_mine = NULL,Bexp,FDexp; 2177 MatStructure mstruct; 2178 PetscViewer vdraw,vstdout; 2179 PetscBool flg; 2180 if (flag_operator) { 2181 ierr = MatComputeExplicitOperator(*A,&Bexp_mine);CHKERRQ(ierr); 2182 Bexp = Bexp_mine; 2183 } else { 2184 /* See if the preconditioning matrix can be viewed and added directly */ 2185 ierr = PetscObjectTypeCompareAny((PetscObject)*B,&flg,MATSEQAIJ,MATMPIAIJ,MATSEQDENSE,MATMPIDENSE,MATSEQBAIJ,MATMPIBAIJ,MATSEQSBAIJ,MATMPIBAIJ,"");CHKERRQ(ierr); 2186 if (flg) Bexp = *B; 2187 else { 2188 /* If the "preconditioning" matrix is itself MATSHELL or some other type without direct support */ 2189 ierr = MatComputeExplicitOperator(*B,&Bexp_mine);CHKERRQ(ierr); 2190 Bexp = Bexp_mine; 2191 } 2192 } 2193 ierr = MatConvert(Bexp,MATSAME,MAT_INITIAL_MATRIX,&FDexp);CHKERRQ(ierr); 2194 ierr = SNESDefaultComputeJacobian(snes,X,&FDexp,&FDexp,&mstruct,NULL);CHKERRQ(ierr); 2195 ierr = PetscViewerASCIIGetStdout(((PetscObject)snes)->comm,&vstdout);CHKERRQ(ierr); 2196 if (flag_draw || flag_contour) { 2197 ierr = PetscViewerDrawOpen(((PetscObject)snes)->comm,0,"Explicit Jacobians",PETSC_DECIDE,PETSC_DECIDE,300,300,&vdraw);CHKERRQ(ierr); 2198 if (flag_contour) {ierr = PetscViewerPushFormat(vdraw,PETSC_VIEWER_DRAW_CONTOUR);CHKERRQ(ierr);} 2199 } else vdraw = NULL; 2200 ierr = PetscViewerASCIIPrintf(vstdout,"Explicit %s\n",flag_operator ? "Jacobian" : "preconditioning Jacobian");CHKERRQ(ierr); 2201 if (flag) {ierr = MatView(Bexp,vstdout);CHKERRQ(ierr);} 2202 if (vdraw) {ierr = MatView(Bexp,vdraw);CHKERRQ(ierr);} 2203 ierr = PetscViewerASCIIPrintf(vstdout,"Finite difference Jacobian\n");CHKERRQ(ierr); 2204 if (flag) {ierr = MatView(FDexp,vstdout);CHKERRQ(ierr);} 2205 if (vdraw) {ierr = MatView(FDexp,vdraw);CHKERRQ(ierr);} 2206 ierr = MatAYPX(FDexp,-1.0,Bexp,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 2207 ierr = PetscViewerASCIIPrintf(vstdout,"User-provided matrix minus finite difference Jacobian\n");CHKERRQ(ierr); 2208 if (flag) {ierr = MatView(FDexp,vstdout);CHKERRQ(ierr);} 2209 if (vdraw) { /* Always use contour for the difference */ 2210 ierr = PetscViewerPushFormat(vdraw,PETSC_VIEWER_DRAW_CONTOUR);CHKERRQ(ierr); 2211 ierr = MatView(FDexp,vdraw);CHKERRQ(ierr); 2212 ierr = PetscViewerPopFormat(vdraw);CHKERRQ(ierr); 2213 } 2214 if (flag_contour) {ierr = PetscViewerPopFormat(vdraw);CHKERRQ(ierr);} 2215 ierr = PetscViewerDestroy(&vdraw);CHKERRQ(ierr); 2216 ierr = MatDestroy(&Bexp_mine);CHKERRQ(ierr); 2217 ierr = MatDestroy(&FDexp);CHKERRQ(ierr); 2218 } 2219 } 2220 { 2221 PetscBool flag = PETSC_FALSE,flag_display = PETSC_FALSE,flag_draw = PETSC_FALSE,flag_contour = PETSC_FALSE,flag_threshold = PETSC_FALSE; 2222 PetscReal threshold_atol = PETSC_SQRT_MACHINE_EPSILON,threshold_rtol = 10*PETSC_SQRT_MACHINE_EPSILON; 2223 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_coloring",&flag,NULL);CHKERRQ(ierr); 2224 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_coloring_display",&flag_display,NULL);CHKERRQ(ierr); 2225 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_coloring_draw",&flag_draw,NULL);CHKERRQ(ierr); 2226 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_coloring_draw_contour",&flag_contour,NULL);CHKERRQ(ierr); 2227 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_coloring_threshold",&flag_threshold,NULL);CHKERRQ(ierr); 2228 ierr = PetscOptionsGetReal(((PetscObject)snes)->prefix,"-snes_compare_coloring_threshold_rtol",&threshold_rtol,NULL);CHKERRQ(ierr); 2229 ierr = PetscOptionsGetReal(((PetscObject)snes)->prefix,"-snes_compare_coloring_threshold_atol",&threshold_atol,NULL);CHKERRQ(ierr); 2230 if (flag || flag_display || flag_draw || flag_contour || flag_threshold) { 2231 Mat Bfd; 2232 MatStructure mstruct; 2233 PetscViewer vdraw,vstdout; 2234 ISColoring iscoloring; 2235 MatFDColoring matfdcoloring; 2236 PetscErrorCode (*func)(SNES,Vec,Vec,void*); 2237 void *funcctx; 2238 PetscReal norm1,norm2,normmax; 2239 2240 ierr = MatDuplicate(*B,MAT_DO_NOT_COPY_VALUES,&Bfd);CHKERRQ(ierr); 2241 ierr = MatGetColoring(Bfd,MATCOLORINGSL,&iscoloring);CHKERRQ(ierr); 2242 ierr = MatFDColoringCreate(Bfd,iscoloring,&matfdcoloring);CHKERRQ(ierr); 2243 ierr = ISColoringDestroy(&iscoloring);CHKERRQ(ierr); 2244 2245 /* This method of getting the function is currently unreliable since it doesn't work for DM local functions. */ 2246 ierr = SNESGetFunction(snes,NULL,&func,&funcctx);CHKERRQ(ierr); 2247 ierr = MatFDColoringSetFunction(matfdcoloring,(PetscErrorCode (*)(void))func,funcctx);CHKERRQ(ierr); 2248 ierr = PetscObjectSetOptionsPrefix((PetscObject)matfdcoloring,((PetscObject)snes)->prefix);CHKERRQ(ierr); 2249 ierr = PetscObjectAppendOptionsPrefix((PetscObject)matfdcoloring,"coloring_");CHKERRQ(ierr); 2250 ierr = MatFDColoringSetFromOptions(matfdcoloring);CHKERRQ(ierr); 2251 ierr = MatFDColoringApply(Bfd,matfdcoloring,X,&mstruct,snes);CHKERRQ(ierr); 2252 ierr = MatFDColoringDestroy(&matfdcoloring);CHKERRQ(ierr); 2253 2254 ierr = PetscViewerASCIIGetStdout(((PetscObject)snes)->comm,&vstdout);CHKERRQ(ierr); 2255 if (flag_draw || flag_contour) { 2256 ierr = PetscViewerDrawOpen(((PetscObject)snes)->comm,0,"Colored Jacobians",PETSC_DECIDE,PETSC_DECIDE,300,300,&vdraw);CHKERRQ(ierr); 2257 if (flag_contour) {ierr = PetscViewerPushFormat(vdraw,PETSC_VIEWER_DRAW_CONTOUR);CHKERRQ(ierr);} 2258 } else vdraw = NULL; 2259 ierr = PetscViewerASCIIPrintf(vstdout,"Explicit preconditioning Jacobian\n");CHKERRQ(ierr); 2260 if (flag_display) {ierr = MatView(*B,vstdout);CHKERRQ(ierr);} 2261 if (vdraw) {ierr = MatView(*B,vdraw);CHKERRQ(ierr);} 2262 ierr = PetscViewerASCIIPrintf(vstdout,"Colored Finite difference Jacobian\n");CHKERRQ(ierr); 2263 if (flag_display) {ierr = MatView(Bfd,vstdout);CHKERRQ(ierr);} 2264 if (vdraw) {ierr = MatView(Bfd,vdraw);CHKERRQ(ierr);} 2265 ierr = MatAYPX(Bfd,-1.0,*B,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 2266 ierr = MatNorm(Bfd,NORM_1,&norm1);CHKERRQ(ierr); 2267 ierr = MatNorm(Bfd,NORM_FROBENIUS,&norm2);CHKERRQ(ierr); 2268 ierr = MatNorm(Bfd,NORM_MAX,&normmax);CHKERRQ(ierr); 2269 ierr = PetscViewerASCIIPrintf(vstdout,"User-provided matrix minus finite difference Jacobian, norm1=%G normFrob=%G normmax=%G\n",norm1,norm2,normmax);CHKERRQ(ierr); 2270 if (flag_display) {ierr = MatView(Bfd,vstdout);CHKERRQ(ierr);} 2271 if (vdraw) { /* Always use contour for the difference */ 2272 ierr = PetscViewerPushFormat(vdraw,PETSC_VIEWER_DRAW_CONTOUR);CHKERRQ(ierr); 2273 ierr = MatView(Bfd,vdraw);CHKERRQ(ierr); 2274 ierr = PetscViewerPopFormat(vdraw);CHKERRQ(ierr); 2275 } 2276 if (flag_contour) {ierr = PetscViewerPopFormat(vdraw);CHKERRQ(ierr);} 2277 2278 if (flag_threshold) { 2279 PetscInt bs,rstart,rend,i; 2280 ierr = MatGetBlockSize(*B,&bs);CHKERRQ(ierr); 2281 ierr = MatGetOwnershipRange(*B,&rstart,&rend);CHKERRQ(ierr); 2282 for (i=rstart; i<rend; i++) { 2283 const PetscScalar *ba,*ca; 2284 const PetscInt *bj,*cj; 2285 PetscInt bn,cn,j,maxentrycol = -1,maxdiffcol = -1,maxrdiffcol = -1; 2286 PetscReal maxentry = 0,maxdiff = 0,maxrdiff = 0; 2287 ierr = MatGetRow(*B,i,&bn,&bj,&ba);CHKERRQ(ierr); 2288 ierr = MatGetRow(Bfd,i,&cn,&cj,&ca);CHKERRQ(ierr); 2289 if (bn != cn) SETERRQ(((PetscObject)*A)->comm,PETSC_ERR_PLIB,"Unexpected different nonzero pattern in -snes_compare_coloring_threshold"); 2290 for (j=0; j<bn; j++) { 2291 PetscReal rdiff = PetscAbsScalar(ca[j]) / (threshold_atol + threshold_rtol*PetscAbsScalar(ba[j])); 2292 if (PetscAbsScalar(ba[j]) > PetscAbs(maxentry)) { 2293 maxentrycol = bj[j]; 2294 maxentry = PetscRealPart(ba[j]); 2295 } 2296 if (PetscAbsScalar(ca[j]) > PetscAbs(maxdiff)) { 2297 maxdiffcol = bj[j]; 2298 maxdiff = PetscRealPart(ca[j]); 2299 } 2300 if (rdiff > maxrdiff) { 2301 maxrdiffcol = bj[j]; 2302 maxrdiff = rdiff; 2303 } 2304 } 2305 if (maxrdiff > 1) { 2306 ierr = PetscViewerASCIIPrintf(vstdout,"row %D (maxentry=%G at %D, maxdiff=%G at %D, maxrdiff=%G at %D):",i,maxentry,maxentrycol,maxdiff,maxdiffcol,maxrdiff,maxrdiffcol);CHKERRQ(ierr); 2307 for (j=0; j<bn; j++) { 2308 PetscReal rdiff; 2309 rdiff = PetscAbsScalar(ca[j]) / (threshold_atol + threshold_rtol*PetscAbsScalar(ba[j])); 2310 if (rdiff > 1) { 2311 ierr = PetscViewerASCIIPrintf(vstdout," (%D,%G:%G)",bj[j],PetscRealPart(ba[j]),PetscRealPart(ca[j]));CHKERRQ(ierr); 2312 } 2313 } 2314 ierr = PetscViewerASCIIPrintf(vstdout,"\n",i,maxentry,maxdiff,maxrdiff);CHKERRQ(ierr); 2315 } 2316 ierr = MatRestoreRow(*B,i,&bn,&bj,&ba);CHKERRQ(ierr); 2317 ierr = MatRestoreRow(Bfd,i,&cn,&cj,&ca);CHKERRQ(ierr); 2318 } 2319 } 2320 ierr = PetscViewerDestroy(&vdraw);CHKERRQ(ierr); 2321 ierr = MatDestroy(&Bfd);CHKERRQ(ierr); 2322 } 2323 } 2324 PetscFunctionReturn(0); 2325 } 2326 2327 /*MC 2328 SNESJacobianFunction - function used to convey the nonlinear Jacobian of the function to be solved by SNES 2329 2330 Synopsis: 2331 #include "petscsnes.h" 2332 $ SNESJacobianFunction(SNES snes,Vec x,Mat *A,Mat *B,int *flag,void *ctx); 2333 2334 + x - input vector 2335 . A - Jacobian matrix 2336 . B - preconditioner matrix, usually the same as A 2337 . flag - flag indicating information about the preconditioner matrix 2338 structure (same as flag in KSPSetOperators()), one of SAME_NONZERO_PATTERN,DIFFERENT_NONZERO_PATTERN,SAME_PRECONDITIONER 2339 - ctx - [optional] user-defined Jacobian context 2340 2341 .seealso: SNESSetFunction(), SNESGetFunction(), SNESSetJacobian(), SNESGetJacobian() 2342 M*/ 2343 2344 #undef __FUNCT__ 2345 #define __FUNCT__ "SNESSetJacobian" 2346 /*@C 2347 SNESSetJacobian - Sets the function to compute Jacobian as well as the 2348 location to store the matrix. 2349 2350 Logically Collective on SNES and Mat 2351 2352 Input Parameters: 2353 + snes - the SNES context 2354 . A - Jacobian matrix 2355 . B - preconditioner matrix (usually same as the Jacobian) 2356 . SNESJacobianFunction - Jacobian evaluation routine (if NULL then SNES retains any previously set value) 2357 - ctx - [optional] user-defined context for private data for the 2358 Jacobian evaluation routine (may be NULL) (if NULL then SNES retains any previously set value) 2359 2360 Notes: 2361 See KSPSetOperators() for important information about setting the flag 2362 output parameter in the routine func(). Be sure to read this information! 2363 2364 The routine func() takes Mat * as the matrix arguments rather than Mat. 2365 This allows the Jacobian evaluation routine to replace A and/or B with a 2366 completely new new matrix structure (not just different matrix elements) 2367 when appropriate, for instance, if the nonzero structure is changing 2368 throughout the global iterations. 2369 2370 If the A matrix and B matrix are different you must call MatAssemblyBegin/End() on 2371 each matrix. 2372 2373 If using SNESDefaultComputeJacobianColor() to assemble a Jacobian, the ctx argument 2374 must be a MatFDColoring. 2375 2376 Other defect-correction schemes can be used by computing a different matrix in place of the Jacobian. One common 2377 example is to use the "Picard linearization" which only differentiates through the highest order parts of each term. 2378 2379 Level: beginner 2380 2381 .keywords: SNES, nonlinear, set, Jacobian, matrix 2382 2383 .seealso: KSPSetOperators(), SNESSetFunction(), MatMFFDComputeJacobian(), SNESDefaultComputeJacobianColor(), MatStructure, SNESJacobianFunction 2384 @*/ 2385 PetscErrorCode SNESSetJacobian(SNES snes,Mat A,Mat B,PetscErrorCode (*SNESJacobianFunction)(SNES,Vec,Mat*,Mat*,MatStructure*,void*),void *ctx) 2386 { 2387 PetscErrorCode ierr; 2388 DM dm; 2389 2390 PetscFunctionBegin; 2391 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2392 if (A) PetscValidHeaderSpecific(A,MAT_CLASSID,2); 2393 if (B) PetscValidHeaderSpecific(B,MAT_CLASSID,3); 2394 if (A) PetscCheckSameComm(snes,1,A,2); 2395 if (B) PetscCheckSameComm(snes,1,B,3); 2396 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 2397 ierr = DMSNESSetJacobian(dm,SNESJacobianFunction,ctx);CHKERRQ(ierr); 2398 if (A) { 2399 ierr = PetscObjectReference((PetscObject)A);CHKERRQ(ierr); 2400 ierr = MatDestroy(&snes->jacobian);CHKERRQ(ierr); 2401 2402 snes->jacobian = A; 2403 } 2404 if (B) { 2405 ierr = PetscObjectReference((PetscObject)B);CHKERRQ(ierr); 2406 ierr = MatDestroy(&snes->jacobian_pre);CHKERRQ(ierr); 2407 2408 snes->jacobian_pre = B; 2409 } 2410 PetscFunctionReturn(0); 2411 } 2412 2413 #undef __FUNCT__ 2414 #define __FUNCT__ "SNESGetJacobian" 2415 /*@C 2416 SNESGetJacobian - Returns the Jacobian matrix and optionally the user 2417 provided context for evaluating the Jacobian. 2418 2419 Not Collective, but Mat object will be parallel if SNES object is 2420 2421 Input Parameter: 2422 . snes - the nonlinear solver context 2423 2424 Output Parameters: 2425 + A - location to stash Jacobian matrix (or NULL) 2426 . B - location to stash preconditioner matrix (or NULL) 2427 . SNESJacobianFunction - location to put Jacobian function (or NULL) 2428 - ctx - location to stash Jacobian ctx (or NULL) 2429 2430 Level: advanced 2431 2432 .seealso: SNESSetJacobian(), SNESComputeJacobian() 2433 @*/ 2434 PetscErrorCode SNESGetJacobian(SNES snes,Mat *A,Mat *B,PetscErrorCode (**SNESJacobianFunction)(SNES,Vec,Mat*,Mat*,MatStructure*,void*),void **ctx) 2435 { 2436 PetscErrorCode ierr; 2437 DM dm; 2438 DMSNES sdm; 2439 2440 PetscFunctionBegin; 2441 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2442 if (A) *A = snes->jacobian; 2443 if (B) *B = snes->jacobian_pre; 2444 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 2445 ierr = DMGetDMSNES(dm,&sdm);CHKERRQ(ierr); 2446 if (SNESJacobianFunction) *SNESJacobianFunction = sdm->ops->computejacobian; 2447 if (ctx) *ctx = sdm->jacobianctx; 2448 PetscFunctionReturn(0); 2449 } 2450 2451 #undef __FUNCT__ 2452 #define __FUNCT__ "SNESSetUp" 2453 /*@ 2454 SNESSetUp - Sets up the internal data structures for the later use 2455 of a nonlinear solver. 2456 2457 Collective on SNES 2458 2459 Input Parameters: 2460 . snes - the SNES context 2461 2462 Notes: 2463 For basic use of the SNES solvers the user need not explicitly call 2464 SNESSetUp(), since these actions will automatically occur during 2465 the call to SNESSolve(). However, if one wishes to control this 2466 phase separately, SNESSetUp() should be called after SNESCreate() 2467 and optional routines of the form SNESSetXXX(), but before SNESSolve(). 2468 2469 Level: advanced 2470 2471 .keywords: SNES, nonlinear, setup 2472 2473 .seealso: SNESCreate(), SNESSolve(), SNESDestroy() 2474 @*/ 2475 PetscErrorCode SNESSetUp(SNES snes) 2476 { 2477 PetscErrorCode ierr; 2478 DM dm; 2479 DMSNES sdm; 2480 SNESLineSearch linesearch, pclinesearch; 2481 void *lsprectx,*lspostctx; 2482 PetscErrorCode (*precheck)(SNESLineSearch,Vec,Vec,PetscBool*,void*); 2483 PetscErrorCode (*postcheck)(SNESLineSearch,Vec,Vec,Vec,PetscBool*,PetscBool*,void*); 2484 PetscErrorCode (*func)(SNES,Vec,Vec,void*); 2485 Vec f,fpc; 2486 void *funcctx; 2487 PetscErrorCode (*jac)(SNES,Vec,Mat*,Mat*,MatStructure*,void*); 2488 void *jacctx,*appctx; 2489 Mat A,B; 2490 2491 PetscFunctionBegin; 2492 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2493 if (snes->setupcalled) PetscFunctionReturn(0); 2494 2495 if (!((PetscObject)snes)->type_name) { 2496 ierr = SNESSetType(snes,SNESNEWTONLS);CHKERRQ(ierr); 2497 } 2498 2499 ierr = SNESGetFunction(snes,&snes->vec_func,NULL,NULL);CHKERRQ(ierr); 2500 if (snes->vec_func == snes->vec_sol) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_IDN,"Solution vector cannot be function vector"); 2501 if (snes->vec_rhs == snes->vec_sol) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_IDN,"Solution vector cannot be right hand side vector"); 2502 2503 if (!snes->vec_sol_update /* && snes->vec_sol */) { 2504 ierr = VecDuplicate(snes->vec_sol,&snes->vec_sol_update);CHKERRQ(ierr); 2505 ierr = PetscLogObjectParent(snes,snes->vec_sol_update);CHKERRQ(ierr); 2506 } 2507 2508 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 2509 ierr = DMShellSetGlobalVector(snes->dm,snes->vec_sol);CHKERRQ(ierr); 2510 ierr = DMGetDMSNES(dm,&sdm);CHKERRQ(ierr); 2511 if (!sdm->ops->computefunction) SETERRQ(((PetscObject)dm)->comm,PETSC_ERR_ARG_WRONGSTATE,"Function never provided to SNES object"); 2512 if (!sdm->ops->computejacobian) { 2513 ierr = DMSNESSetJacobian(dm,SNESDefaultComputeJacobianColor,NULL);CHKERRQ(ierr); 2514 } 2515 if (!snes->vec_func) { 2516 ierr = DMCreateGlobalVector(dm,&snes->vec_func);CHKERRQ(ierr); 2517 } 2518 2519 if (!snes->ksp) { 2520 ierr = SNESGetKSP(snes, &snes->ksp);CHKERRQ(ierr); 2521 } 2522 2523 if (!snes->linesearch) { 2524 ierr = SNESGetSNESLineSearch(snes, &snes->linesearch);CHKERRQ(ierr); 2525 } 2526 2527 if (snes->pc && (snes->pcside == PC_LEFT)) { 2528 snes->mf = PETSC_TRUE; 2529 snes->mf_operator = PETSC_FALSE; 2530 } 2531 2532 if (snes->mf) { 2533 ierr = SNESSetUpMatrixFree_Private(snes, snes->mf_operator, snes->mf_version);CHKERRQ(ierr); 2534 } 2535 2536 if (snes->ops->usercompute && !snes->user) { 2537 ierr = (*snes->ops->usercompute)(snes,(void**)&snes->user);CHKERRQ(ierr); 2538 } 2539 2540 if (snes->pc) { 2541 /* copy the DM over */ 2542 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 2543 ierr = SNESSetDM(snes->pc,dm);CHKERRQ(ierr); 2544 2545 ierr = SNESGetFunction(snes,&f,&func,&funcctx);CHKERRQ(ierr); 2546 ierr = VecDuplicate(f,&fpc);CHKERRQ(ierr); 2547 ierr = SNESSetFunction(snes->pc,fpc,func,funcctx);CHKERRQ(ierr); 2548 ierr = SNESGetJacobian(snes,&A,&B,&jac,&jacctx);CHKERRQ(ierr); 2549 ierr = SNESSetJacobian(snes->pc,A,B,jac,jacctx);CHKERRQ(ierr); 2550 ierr = SNESGetApplicationContext(snes,&appctx);CHKERRQ(ierr); 2551 ierr = SNESSetApplicationContext(snes->pc,appctx);CHKERRQ(ierr); 2552 ierr = VecDestroy(&fpc);CHKERRQ(ierr); 2553 2554 /* copy the function pointers over */ 2555 ierr = PetscObjectCopyFortranFunctionPointers((PetscObject)snes,(PetscObject)snes->pc);CHKERRQ(ierr); 2556 2557 /* default to 1 iteration */ 2558 ierr = SNESSetTolerances(snes->pc,0.0,0.0,0.0,1,snes->pc->max_funcs);CHKERRQ(ierr); 2559 if (snes->pcside==PC_RIGHT) { 2560 ierr = SNESSetNormType(snes->pc,SNES_NORM_FINAL_ONLY);CHKERRQ(ierr); 2561 } else { 2562 ierr = SNESSetNormType(snes->pc,SNES_NORM_NONE);CHKERRQ(ierr); 2563 } 2564 ierr = SNESSetFromOptions(snes->pc);CHKERRQ(ierr); 2565 2566 /* copy the line search context over */ 2567 ierr = SNESGetSNESLineSearch(snes,&linesearch);CHKERRQ(ierr); 2568 ierr = SNESGetSNESLineSearch(snes->pc,&pclinesearch);CHKERRQ(ierr); 2569 ierr = SNESLineSearchGetPreCheck(linesearch,&precheck,&lsprectx);CHKERRQ(ierr); 2570 ierr = SNESLineSearchGetPostCheck(linesearch,&postcheck,&lspostctx);CHKERRQ(ierr); 2571 ierr = SNESLineSearchSetPreCheck(pclinesearch,precheck,lsprectx);CHKERRQ(ierr); 2572 ierr = SNESLineSearchSetPostCheck(pclinesearch,postcheck,lspostctx);CHKERRQ(ierr); 2573 ierr = PetscObjectCopyFortranFunctionPointers((PetscObject)linesearch, (PetscObject)pclinesearch);CHKERRQ(ierr); 2574 } 2575 2576 if (snes->ops->setup) { 2577 ierr = (*snes->ops->setup)(snes);CHKERRQ(ierr); 2578 } 2579 2580 snes->setupcalled = PETSC_TRUE; 2581 PetscFunctionReturn(0); 2582 } 2583 2584 #undef __FUNCT__ 2585 #define __FUNCT__ "SNESReset" 2586 /*@ 2587 SNESReset - Resets a SNES context to the snessetupcalled = 0 state and removes any allocated Vecs and Mats 2588 2589 Collective on SNES 2590 2591 Input Parameter: 2592 . snes - iterative context obtained from SNESCreate() 2593 2594 Level: intermediate 2595 2596 Notes: Also calls the application context destroy routine set with SNESSetComputeApplicationContext() 2597 2598 .keywords: SNES, destroy 2599 2600 .seealso: SNESCreate(), SNESSetUp(), SNESSolve() 2601 @*/ 2602 PetscErrorCode SNESReset(SNES snes) 2603 { 2604 PetscErrorCode ierr; 2605 2606 PetscFunctionBegin; 2607 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2608 if (snes->ops->userdestroy && snes->user) { 2609 ierr = (*snes->ops->userdestroy)((void**)&snes->user);CHKERRQ(ierr); 2610 snes->user = NULL; 2611 } 2612 if (snes->pc) { 2613 ierr = SNESReset(snes->pc);CHKERRQ(ierr); 2614 } 2615 2616 if (snes->ops->reset) { 2617 ierr = (*snes->ops->reset)(snes);CHKERRQ(ierr); 2618 } 2619 if (snes->ksp) { 2620 ierr = KSPReset(snes->ksp);CHKERRQ(ierr); 2621 } 2622 2623 if (snes->linesearch) { 2624 ierr = SNESLineSearchReset(snes->linesearch);CHKERRQ(ierr); 2625 } 2626 2627 ierr = VecDestroy(&snes->vec_rhs);CHKERRQ(ierr); 2628 ierr = VecDestroy(&snes->vec_sol);CHKERRQ(ierr); 2629 ierr = VecDestroy(&snes->vec_sol_update);CHKERRQ(ierr); 2630 ierr = VecDestroy(&snes->vec_func);CHKERRQ(ierr); 2631 ierr = MatDestroy(&snes->jacobian);CHKERRQ(ierr); 2632 ierr = MatDestroy(&snes->jacobian_pre);CHKERRQ(ierr); 2633 ierr = VecDestroyVecs(snes->nwork,&snes->work);CHKERRQ(ierr); 2634 ierr = VecDestroyVecs(snes->nvwork,&snes->vwork);CHKERRQ(ierr); 2635 2636 snes->nwork = snes->nvwork = 0; 2637 snes->setupcalled = PETSC_FALSE; 2638 PetscFunctionReturn(0); 2639 } 2640 2641 #undef __FUNCT__ 2642 #define __FUNCT__ "SNESDestroy" 2643 /*@ 2644 SNESDestroy - Destroys the nonlinear solver context that was created 2645 with SNESCreate(). 2646 2647 Collective on SNES 2648 2649 Input Parameter: 2650 . snes - the SNES context 2651 2652 Level: beginner 2653 2654 .keywords: SNES, nonlinear, destroy 2655 2656 .seealso: SNESCreate(), SNESSolve() 2657 @*/ 2658 PetscErrorCode SNESDestroy(SNES *snes) 2659 { 2660 PetscErrorCode ierr; 2661 2662 PetscFunctionBegin; 2663 if (!*snes) PetscFunctionReturn(0); 2664 PetscValidHeaderSpecific((*snes),SNES_CLASSID,1); 2665 if (--((PetscObject)(*snes))->refct > 0) {*snes = 0; PetscFunctionReturn(0);} 2666 2667 ierr = SNESReset((*snes));CHKERRQ(ierr); 2668 ierr = SNESDestroy(&(*snes)->pc);CHKERRQ(ierr); 2669 2670 /* if memory was published with AMS then destroy it */ 2671 ierr = PetscObjectDepublish((*snes));CHKERRQ(ierr); 2672 if ((*snes)->ops->destroy) {ierr = (*((*snes))->ops->destroy)((*snes));CHKERRQ(ierr);} 2673 2674 ierr = DMDestroy(&(*snes)->dm);CHKERRQ(ierr); 2675 ierr = KSPDestroy(&(*snes)->ksp);CHKERRQ(ierr); 2676 ierr = SNESLineSearchDestroy(&(*snes)->linesearch);CHKERRQ(ierr); 2677 2678 ierr = PetscFree((*snes)->kspconvctx);CHKERRQ(ierr); 2679 if ((*snes)->ops->convergeddestroy) { 2680 ierr = (*(*snes)->ops->convergeddestroy)((*snes)->cnvP);CHKERRQ(ierr); 2681 } 2682 if ((*snes)->conv_malloc) { 2683 ierr = PetscFree((*snes)->conv_hist);CHKERRQ(ierr); 2684 ierr = PetscFree((*snes)->conv_hist_its);CHKERRQ(ierr); 2685 } 2686 ierr = SNESMonitorCancel((*snes));CHKERRQ(ierr); 2687 ierr = PetscHeaderDestroy(snes);CHKERRQ(ierr); 2688 PetscFunctionReturn(0); 2689 } 2690 2691 /* ----------- Routines to set solver parameters ---------- */ 2692 2693 #undef __FUNCT__ 2694 #define __FUNCT__ "SNESSetLagPreconditioner" 2695 /*@ 2696 SNESSetLagPreconditioner - Determines when the preconditioner is rebuilt in the nonlinear solve. 2697 2698 Logically Collective on SNES 2699 2700 Input Parameters: 2701 + snes - the SNES context 2702 - lag - -1 indicates NEVER rebuild, 1 means rebuild every time the Jacobian is computed within a single nonlinear solve, 2 means every second time 2703 the Jacobian is built etc. -2 indicates rebuild preconditioner at next chance but then never rebuild after that 2704 2705 Options Database Keys: 2706 . -snes_lag_preconditioner <lag> 2707 2708 Notes: 2709 The default is 1 2710 The preconditioner is ALWAYS built in the first iteration of a nonlinear solve unless lag is -1 2711 If -1 is used before the very first nonlinear solve the preconditioner is still built because there is no previous preconditioner to use 2712 2713 Level: intermediate 2714 2715 .keywords: SNES, nonlinear, set, convergence, tolerances 2716 2717 .seealso: SNESSetTrustRegionTolerance(), SNESGetLagPreconditioner(), SNESSetLagJacobian(), SNESGetLagJacobian() 2718 2719 @*/ 2720 PetscErrorCode SNESSetLagPreconditioner(SNES snes,PetscInt lag) 2721 { 2722 PetscFunctionBegin; 2723 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2724 if (lag < -2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Lag must be -2, -1, 1 or greater"); 2725 if (!lag) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Lag cannot be 0"); 2726 PetscValidLogicalCollectiveInt(snes,lag,2); 2727 snes->lagpreconditioner = lag; 2728 PetscFunctionReturn(0); 2729 } 2730 2731 #undef __FUNCT__ 2732 #define __FUNCT__ "SNESSetGridSequence" 2733 /*@ 2734 SNESSetGridSequence - sets the number of steps of grid sequencing that SNES does 2735 2736 Logically Collective on SNES 2737 2738 Input Parameters: 2739 + snes - the SNES context 2740 - steps - the number of refinements to do, defaults to 0 2741 2742 Options Database Keys: 2743 . -snes_grid_sequence <steps> 2744 2745 Level: intermediate 2746 2747 Notes: 2748 Use SNESGetSolution() to extract the fine grid solution after grid sequencing. 2749 2750 .keywords: SNES, nonlinear, set, convergence, tolerances 2751 2752 .seealso: SNESSetTrustRegionTolerance(), SNESGetLagPreconditioner(), SNESSetLagJacobian(), SNESGetLagJacobian() 2753 2754 @*/ 2755 PetscErrorCode SNESSetGridSequence(SNES snes,PetscInt steps) 2756 { 2757 PetscFunctionBegin; 2758 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2759 PetscValidLogicalCollectiveInt(snes,steps,2); 2760 snes->gridsequence = steps; 2761 PetscFunctionReturn(0); 2762 } 2763 2764 #undef __FUNCT__ 2765 #define __FUNCT__ "SNESGetLagPreconditioner" 2766 /*@ 2767 SNESGetLagPreconditioner - Indicates how often the preconditioner is rebuilt 2768 2769 Not Collective 2770 2771 Input Parameter: 2772 . snes - the SNES context 2773 2774 Output Parameter: 2775 . lag - -1 indicates NEVER rebuild, 1 means rebuild every time the Jacobian is computed within a single nonlinear solve, 2 means every second time 2776 the Jacobian is built etc. -2 indicates rebuild preconditioner at next chance but then never rebuild after that 2777 2778 Options Database Keys: 2779 . -snes_lag_preconditioner <lag> 2780 2781 Notes: 2782 The default is 1 2783 The preconditioner is ALWAYS built in the first iteration of a nonlinear solve unless lag is -1 2784 2785 Level: intermediate 2786 2787 .keywords: SNES, nonlinear, set, convergence, tolerances 2788 2789 .seealso: SNESSetTrustRegionTolerance(), SNESSetLagPreconditioner() 2790 2791 @*/ 2792 PetscErrorCode SNESGetLagPreconditioner(SNES snes,PetscInt *lag) 2793 { 2794 PetscFunctionBegin; 2795 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2796 *lag = snes->lagpreconditioner; 2797 PetscFunctionReturn(0); 2798 } 2799 2800 #undef __FUNCT__ 2801 #define __FUNCT__ "SNESSetLagJacobian" 2802 /*@ 2803 SNESSetLagJacobian - Determines when the Jacobian is rebuilt in the nonlinear solve. See SNESSetLagPreconditioner() for determining how 2804 often the preconditioner is rebuilt. 2805 2806 Logically Collective on SNES 2807 2808 Input Parameters: 2809 + snes - the SNES context 2810 - lag - -1 indicates NEVER rebuild, 1 means rebuild every time the Jacobian is computed within a single nonlinear solve, 2 means every second time 2811 the Jacobian is built etc. -2 means rebuild at next chance but then never again 2812 2813 Options Database Keys: 2814 . -snes_lag_jacobian <lag> 2815 2816 Notes: 2817 The default is 1 2818 The Jacobian is ALWAYS built in the first iteration of a nonlinear solve unless lag is -1 2819 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 2820 at the next Newton step but never again (unless it is reset to another value) 2821 2822 Level: intermediate 2823 2824 .keywords: SNES, nonlinear, set, convergence, tolerances 2825 2826 .seealso: SNESSetTrustRegionTolerance(), SNESGetLagPreconditioner(), SNESSetLagPreconditioner(), SNESGetLagJacobian() 2827 2828 @*/ 2829 PetscErrorCode SNESSetLagJacobian(SNES snes,PetscInt lag) 2830 { 2831 PetscFunctionBegin; 2832 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2833 if (lag < -2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Lag must be -2, -1, 1 or greater"); 2834 if (!lag) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Lag cannot be 0"); 2835 PetscValidLogicalCollectiveInt(snes,lag,2); 2836 snes->lagjacobian = lag; 2837 PetscFunctionReturn(0); 2838 } 2839 2840 #undef __FUNCT__ 2841 #define __FUNCT__ "SNESGetLagJacobian" 2842 /*@ 2843 SNESGetLagJacobian - Indicates how often the Jacobian is rebuilt. See SNESGetLagPreconditioner() to determine when the preconditioner is rebuilt 2844 2845 Not Collective 2846 2847 Input Parameter: 2848 . snes - the SNES context 2849 2850 Output Parameter: 2851 . lag - -1 indicates NEVER rebuild, 1 means rebuild every time the Jacobian is computed within a single nonlinear solve, 2 means every second time 2852 the Jacobian is built etc. 2853 2854 Options Database Keys: 2855 . -snes_lag_jacobian <lag> 2856 2857 Notes: 2858 The default is 1 2859 The jacobian is ALWAYS built in the first iteration of a nonlinear solve unless lag is -1 2860 2861 Level: intermediate 2862 2863 .keywords: SNES, nonlinear, set, convergence, tolerances 2864 2865 .seealso: SNESSetTrustRegionTolerance(), SNESSetLagJacobian(), SNESSetLagPreconditioner(), SNESGetLagPreconditioner() 2866 2867 @*/ 2868 PetscErrorCode SNESGetLagJacobian(SNES snes,PetscInt *lag) 2869 { 2870 PetscFunctionBegin; 2871 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2872 *lag = snes->lagjacobian; 2873 PetscFunctionReturn(0); 2874 } 2875 2876 #undef __FUNCT__ 2877 #define __FUNCT__ "SNESSetTolerances" 2878 /*@ 2879 SNESSetTolerances - Sets various parameters used in convergence tests. 2880 2881 Logically Collective on SNES 2882 2883 Input Parameters: 2884 + snes - the SNES context 2885 . abstol - absolute convergence tolerance 2886 . rtol - relative convergence tolerance 2887 . stol - convergence tolerance in terms of the norm of the change in the solution between steps, || delta x || < stol*|| x || 2888 . maxit - maximum number of iterations 2889 - maxf - maximum number of function evaluations 2890 2891 Options Database Keys: 2892 + -snes_atol <abstol> - Sets abstol 2893 . -snes_rtol <rtol> - Sets rtol 2894 . -snes_stol <stol> - Sets stol 2895 . -snes_max_it <maxit> - Sets maxit 2896 - -snes_max_funcs <maxf> - Sets maxf 2897 2898 Notes: 2899 The default maximum number of iterations is 50. 2900 The default maximum number of function evaluations is 1000. 2901 2902 Level: intermediate 2903 2904 .keywords: SNES, nonlinear, set, convergence, tolerances 2905 2906 .seealso: SNESSetTrustRegionTolerance() 2907 @*/ 2908 PetscErrorCode SNESSetTolerances(SNES snes,PetscReal abstol,PetscReal rtol,PetscReal stol,PetscInt maxit,PetscInt maxf) 2909 { 2910 PetscFunctionBegin; 2911 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2912 PetscValidLogicalCollectiveReal(snes,abstol,2); 2913 PetscValidLogicalCollectiveReal(snes,rtol,3); 2914 PetscValidLogicalCollectiveReal(snes,stol,4); 2915 PetscValidLogicalCollectiveInt(snes,maxit,5); 2916 PetscValidLogicalCollectiveInt(snes,maxf,6); 2917 2918 if (abstol != PETSC_DEFAULT) { 2919 if (abstol < 0.0) SETERRQ1(((PetscObject)snes)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Absolute tolerance %G must be non-negative",abstol); 2920 snes->abstol = abstol; 2921 } 2922 if (rtol != PETSC_DEFAULT) { 2923 if (rtol < 0.0 || 1.0 <= rtol) SETERRQ1(((PetscObject)snes)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Relative tolerance %G must be non-negative and less than 1.0",rtol); 2924 snes->rtol = rtol; 2925 } 2926 if (stol != PETSC_DEFAULT) { 2927 if (stol < 0.0) SETERRQ1(((PetscObject)snes)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Step tolerance %G must be non-negative",stol); 2928 snes->stol = stol; 2929 } 2930 if (maxit != PETSC_DEFAULT) { 2931 if (maxit < 0) SETERRQ1(((PetscObject)snes)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Maximum number of iterations %D must be non-negative",maxit); 2932 snes->max_its = maxit; 2933 } 2934 if (maxf != PETSC_DEFAULT) { 2935 if (maxf < 0) SETERRQ1(((PetscObject)snes)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Maximum number of function evaluations %D must be non-negative",maxf); 2936 snes->max_funcs = maxf; 2937 } 2938 snes->tolerancesset = PETSC_TRUE; 2939 PetscFunctionReturn(0); 2940 } 2941 2942 #undef __FUNCT__ 2943 #define __FUNCT__ "SNESGetTolerances" 2944 /*@ 2945 SNESGetTolerances - Gets various parameters used in convergence tests. 2946 2947 Not Collective 2948 2949 Input Parameters: 2950 + snes - the SNES context 2951 . atol - absolute convergence tolerance 2952 . rtol - relative convergence tolerance 2953 . stol - convergence tolerance in terms of the norm 2954 of the change in the solution between steps 2955 . maxit - maximum number of iterations 2956 - maxf - maximum number of function evaluations 2957 2958 Notes: 2959 The user can specify NULL for any parameter that is not needed. 2960 2961 Level: intermediate 2962 2963 .keywords: SNES, nonlinear, get, convergence, tolerances 2964 2965 .seealso: SNESSetTolerances() 2966 @*/ 2967 PetscErrorCode SNESGetTolerances(SNES snes,PetscReal *atol,PetscReal *rtol,PetscReal *stol,PetscInt *maxit,PetscInt *maxf) 2968 { 2969 PetscFunctionBegin; 2970 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2971 if (atol) *atol = snes->abstol; 2972 if (rtol) *rtol = snes->rtol; 2973 if (stol) *stol = snes->stol; 2974 if (maxit) *maxit = snes->max_its; 2975 if (maxf) *maxf = snes->max_funcs; 2976 PetscFunctionReturn(0); 2977 } 2978 2979 #undef __FUNCT__ 2980 #define __FUNCT__ "SNESSetTrustRegionTolerance" 2981 /*@ 2982 SNESSetTrustRegionTolerance - Sets the trust region parameter tolerance. 2983 2984 Logically Collective on SNES 2985 2986 Input Parameters: 2987 + snes - the SNES context 2988 - tol - tolerance 2989 2990 Options Database Key: 2991 . -snes_trtol <tol> - Sets tol 2992 2993 Level: intermediate 2994 2995 .keywords: SNES, nonlinear, set, trust region, tolerance 2996 2997 .seealso: SNESSetTolerances() 2998 @*/ 2999 PetscErrorCode SNESSetTrustRegionTolerance(SNES snes,PetscReal tol) 3000 { 3001 PetscFunctionBegin; 3002 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3003 PetscValidLogicalCollectiveReal(snes,tol,2); 3004 snes->deltatol = tol; 3005 PetscFunctionReturn(0); 3006 } 3007 3008 /* 3009 Duplicate the lg monitors for SNES from KSP; for some reason with 3010 dynamic libraries things don't work under Sun4 if we just use 3011 macros instead of functions 3012 */ 3013 #undef __FUNCT__ 3014 #define __FUNCT__ "SNESMonitorLGResidualNorm" 3015 PetscErrorCode SNESMonitorLGResidualNorm(SNES snes,PetscInt it,PetscReal norm,void *ctx) 3016 { 3017 PetscErrorCode ierr; 3018 3019 PetscFunctionBegin; 3020 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3021 ierr = KSPMonitorLGResidualNorm((KSP)snes,it,norm,ctx);CHKERRQ(ierr); 3022 PetscFunctionReturn(0); 3023 } 3024 3025 #undef __FUNCT__ 3026 #define __FUNCT__ "SNESMonitorLGCreate" 3027 PetscErrorCode SNESMonitorLGCreate(const char host[],const char label[],int x,int y,int m,int n,PetscDrawLG *draw) 3028 { 3029 PetscErrorCode ierr; 3030 3031 PetscFunctionBegin; 3032 ierr = KSPMonitorLGResidualNormCreate(host,label,x,y,m,n,draw);CHKERRQ(ierr); 3033 PetscFunctionReturn(0); 3034 } 3035 3036 #undef __FUNCT__ 3037 #define __FUNCT__ "SNESMonitorLGDestroy" 3038 PetscErrorCode SNESMonitorLGDestroy(PetscDrawLG *draw) 3039 { 3040 PetscErrorCode ierr; 3041 3042 PetscFunctionBegin; 3043 ierr = KSPMonitorLGResidualNormDestroy(draw);CHKERRQ(ierr); 3044 PetscFunctionReturn(0); 3045 } 3046 3047 extern PetscErrorCode SNESMonitorRange_Private(SNES,PetscInt,PetscReal*); 3048 #undef __FUNCT__ 3049 #define __FUNCT__ "SNESMonitorLGRange" 3050 PetscErrorCode SNESMonitorLGRange(SNES snes,PetscInt n,PetscReal rnorm,void *monctx) 3051 { 3052 PetscDrawLG lg; 3053 PetscErrorCode ierr; 3054 PetscReal x,y,per; 3055 PetscViewer v = (PetscViewer)monctx; 3056 static PetscReal prev; /* should be in the context */ 3057 PetscDraw draw; 3058 3059 PetscFunctionBegin; 3060 ierr = PetscViewerDrawGetDrawLG(v,0,&lg);CHKERRQ(ierr); 3061 if (!n) {ierr = PetscDrawLGReset(lg);CHKERRQ(ierr);} 3062 ierr = PetscDrawLGGetDraw(lg,&draw);CHKERRQ(ierr); 3063 ierr = PetscDrawSetTitle(draw,"Residual norm");CHKERRQ(ierr); 3064 x = (PetscReal)n; 3065 if (rnorm > 0.0) y = log10(rnorm); 3066 else y = -15.0; 3067 ierr = PetscDrawLGAddPoint(lg,&x,&y);CHKERRQ(ierr); 3068 if (n < 20 || !(n % 5)) { 3069 ierr = PetscDrawLGDraw(lg);CHKERRQ(ierr); 3070 } 3071 3072 ierr = PetscViewerDrawGetDrawLG(v,1,&lg);CHKERRQ(ierr); 3073 if (!n) {ierr = PetscDrawLGReset(lg);CHKERRQ(ierr);} 3074 ierr = PetscDrawLGGetDraw(lg,&draw);CHKERRQ(ierr); 3075 ierr = PetscDrawSetTitle(draw,"% elemts > .2*max elemt");CHKERRQ(ierr); 3076 ierr = SNESMonitorRange_Private(snes,n,&per);CHKERRQ(ierr); 3077 x = (PetscReal)n; 3078 y = 100.0*per; 3079 ierr = PetscDrawLGAddPoint(lg,&x,&y);CHKERRQ(ierr); 3080 if (n < 20 || !(n % 5)) { 3081 ierr = PetscDrawLGDraw(lg);CHKERRQ(ierr); 3082 } 3083 3084 ierr = PetscViewerDrawGetDrawLG(v,2,&lg);CHKERRQ(ierr); 3085 if (!n) {prev = rnorm;ierr = PetscDrawLGReset(lg);CHKERRQ(ierr);} 3086 ierr = PetscDrawLGGetDraw(lg,&draw);CHKERRQ(ierr); 3087 ierr = PetscDrawSetTitle(draw,"(norm -oldnorm)/oldnorm");CHKERRQ(ierr); 3088 x = (PetscReal)n; 3089 y = (prev - rnorm)/prev; 3090 ierr = PetscDrawLGAddPoint(lg,&x,&y);CHKERRQ(ierr); 3091 if (n < 20 || !(n % 5)) { 3092 ierr = PetscDrawLGDraw(lg);CHKERRQ(ierr); 3093 } 3094 3095 ierr = PetscViewerDrawGetDrawLG(v,3,&lg);CHKERRQ(ierr); 3096 if (!n) {ierr = PetscDrawLGReset(lg);CHKERRQ(ierr);} 3097 ierr = PetscDrawLGGetDraw(lg,&draw);CHKERRQ(ierr); 3098 ierr = PetscDrawSetTitle(draw,"(norm -oldnorm)/oldnorm*(% > .2 max)");CHKERRQ(ierr); 3099 x = (PetscReal)n; 3100 y = (prev - rnorm)/(prev*per); 3101 if (n > 2) { /*skip initial crazy value */ 3102 ierr = PetscDrawLGAddPoint(lg,&x,&y);CHKERRQ(ierr); 3103 } 3104 if (n < 20 || !(n % 5)) { 3105 ierr = PetscDrawLGDraw(lg);CHKERRQ(ierr); 3106 } 3107 prev = rnorm; 3108 PetscFunctionReturn(0); 3109 } 3110 3111 #undef __FUNCT__ 3112 #define __FUNCT__ "SNESMonitor" 3113 /*@ 3114 SNESMonitor - runs the user provided monitor routines, if they exist 3115 3116 Collective on SNES 3117 3118 Input Parameters: 3119 + snes - nonlinear solver context obtained from SNESCreate() 3120 . iter - iteration number 3121 - rnorm - relative norm of the residual 3122 3123 Notes: 3124 This routine is called by the SNES implementations. 3125 It does not typically need to be called by the user. 3126 3127 Level: developer 3128 3129 .seealso: SNESMonitorSet() 3130 @*/ 3131 PetscErrorCode SNESMonitor(SNES snes,PetscInt iter,PetscReal rnorm) 3132 { 3133 PetscErrorCode ierr; 3134 PetscInt i,n = snes->numbermonitors; 3135 3136 PetscFunctionBegin; 3137 for (i=0; i<n; i++) { 3138 ierr = (*snes->monitor[i])(snes,iter,rnorm,snes->monitorcontext[i]);CHKERRQ(ierr); 3139 } 3140 PetscFunctionReturn(0); 3141 } 3142 3143 /* ------------ Routines to set performance monitoring options ----------- */ 3144 3145 /*MC 3146 SNESMonitorFunction - functional form passed to SNESMonitorSet() to monitor convergence of nonlinear solver 3147 3148 Synopsis: 3149 #include "petscsnes.h" 3150 $ PetscErrorCode SNESMonitorFunction(SNES snes,PetscInt its, PetscReal norm,void *mctx) 3151 3152 + snes - the SNES context 3153 . its - iteration number 3154 . norm - 2-norm function value (may be estimated) 3155 - mctx - [optional] monitoring context 3156 3157 .seealso: SNESMonitorSet(), SNESMonitorGet() 3158 M*/ 3159 3160 #undef __FUNCT__ 3161 #define __FUNCT__ "SNESMonitorSet" 3162 /*@C 3163 SNESMonitorSet - Sets an ADDITIONAL function that is to be used at every 3164 iteration of the nonlinear solver to display the iteration's 3165 progress. 3166 3167 Logically Collective on SNES 3168 3169 Input Parameters: 3170 + snes - the SNES context 3171 . SNESMonitorFunction - monitoring routine 3172 . mctx - [optional] user-defined context for private data for the 3173 monitor routine (use NULL if no context is desired) 3174 - monitordestroy - [optional] routine that frees monitor context 3175 (may be NULL) 3176 3177 Options Database Keys: 3178 + -snes_monitor - sets SNESMonitorDefault() 3179 . -snes_monitor_lg_residualnorm - sets line graph monitor, 3180 uses SNESMonitorLGCreate() 3181 - -snes_monitor_cancel - cancels all monitors that have 3182 been hardwired into a code by 3183 calls to SNESMonitorSet(), but 3184 does not cancel those set via 3185 the options database. 3186 3187 Notes: 3188 Several different monitoring routines may be set by calling 3189 SNESMonitorSet() multiple times; all will be called in the 3190 order in which they were set. 3191 3192 Fortran notes: Only a single monitor function can be set for each SNES object 3193 3194 Level: intermediate 3195 3196 .keywords: SNES, nonlinear, set, monitor 3197 3198 .seealso: SNESMonitorDefault(), SNESMonitorCancel(), SNESMonitorFunction 3199 @*/ 3200 PetscErrorCode SNESMonitorSet(SNES snes,PetscErrorCode (*SNESMonitorFunction)(SNES,PetscInt,PetscReal,void*),void *mctx,PetscErrorCode (*monitordestroy)(void**)) 3201 { 3202 PetscInt i; 3203 PetscErrorCode ierr; 3204 3205 PetscFunctionBegin; 3206 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3207 if (snes->numbermonitors >= MAXSNESMONITORS) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Too many monitors set"); 3208 for (i=0; i<snes->numbermonitors;i++) { 3209 if (SNESMonitorFunction == snes->monitor[i] && monitordestroy == snes->monitordestroy[i] && mctx == snes->monitorcontext[i]) { 3210 if (monitordestroy) { 3211 ierr = (*monitordestroy)(&mctx);CHKERRQ(ierr); 3212 } 3213 PetscFunctionReturn(0); 3214 } 3215 } 3216 snes->monitor[snes->numbermonitors] = SNESMonitorFunction; 3217 snes->monitordestroy[snes->numbermonitors] = monitordestroy; 3218 snes->monitorcontext[snes->numbermonitors++] = (void*)mctx; 3219 PetscFunctionReturn(0); 3220 } 3221 3222 #undef __FUNCT__ 3223 #define __FUNCT__ "SNESMonitorCancel" 3224 /*@C 3225 SNESMonitorCancel - Clears all the monitor functions for a SNES object. 3226 3227 Logically Collective on SNES 3228 3229 Input Parameters: 3230 . snes - the SNES context 3231 3232 Options Database Key: 3233 . -snes_monitor_cancel - cancels all monitors that have been hardwired 3234 into a code by calls to SNESMonitorSet(), but does not cancel those 3235 set via the options database 3236 3237 Notes: 3238 There is no way to clear one specific monitor from a SNES object. 3239 3240 Level: intermediate 3241 3242 .keywords: SNES, nonlinear, set, monitor 3243 3244 .seealso: SNESMonitorDefault(), SNESMonitorSet() 3245 @*/ 3246 PetscErrorCode SNESMonitorCancel(SNES snes) 3247 { 3248 PetscErrorCode ierr; 3249 PetscInt i; 3250 3251 PetscFunctionBegin; 3252 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3253 for (i=0; i<snes->numbermonitors; i++) { 3254 if (snes->monitordestroy[i]) { 3255 ierr = (*snes->monitordestroy[i])(&snes->monitorcontext[i]);CHKERRQ(ierr); 3256 } 3257 } 3258 snes->numbermonitors = 0; 3259 PetscFunctionReturn(0); 3260 } 3261 3262 /*MC 3263 SNESConvergenceTestFunction - functional form used for testing of convergence of nonlinear solver 3264 3265 Synopsis: 3266 #include "petscsnes.h" 3267 $ PetscErrorCode SNESConvergenceTest(SNES snes,PetscInt it,PetscReal xnorm,PetscReal gnorm,PetscReal f,SNESConvergedReason *reason,void *cctx) 3268 3269 + snes - the SNES context 3270 . it - current iteration (0 is the first and is before any Newton step) 3271 . cctx - [optional] convergence context 3272 . reason - reason for convergence/divergence 3273 . xnorm - 2-norm of current iterate 3274 . gnorm - 2-norm of current step 3275 - f - 2-norm of function 3276 3277 3278 .seealso: SNESSetConvergenceTest(), SNESGetConvergenceTest() 3279 M*/ 3280 3281 #undef __FUNCT__ 3282 #define __FUNCT__ "SNESSetConvergenceTest" 3283 /*@C 3284 SNESSetConvergenceTest - Sets the function that is to be used 3285 to test for convergence of the nonlinear iterative solution. 3286 3287 Logically Collective on SNES 3288 3289 Input Parameters: 3290 + snes - the SNES context 3291 . SNESConvergenceTestFunction - routine to test for convergence 3292 . cctx - [optional] context for private data for the convergence routine (may be NULL) 3293 - destroy - [optional] destructor for the context (may be NULL; NULL_FUNCTION in Fortran) 3294 3295 Level: advanced 3296 3297 .keywords: SNES, nonlinear, set, convergence, test 3298 3299 .seealso: SNESDefaultConverged(), SNESSkipConverged(), SNESConvergenceTestFunction 3300 @*/ 3301 PetscErrorCode SNESSetConvergenceTest(SNES snes,PetscErrorCode (*SNESConvergenceTestFunction)(SNES,PetscInt,PetscReal,PetscReal,PetscReal,SNESConvergedReason*,void*),void *cctx,PetscErrorCode (*destroy)(void*)) 3302 { 3303 PetscErrorCode ierr; 3304 3305 PetscFunctionBegin; 3306 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3307 if (!SNESConvergenceTestFunction) SNESConvergenceTestFunction = SNESSkipConverged; 3308 if (snes->ops->convergeddestroy) { 3309 ierr = (*snes->ops->convergeddestroy)(snes->cnvP);CHKERRQ(ierr); 3310 } 3311 snes->ops->converged = SNESConvergenceTestFunction; 3312 snes->ops->convergeddestroy = destroy; 3313 snes->cnvP = cctx; 3314 PetscFunctionReturn(0); 3315 } 3316 3317 #undef __FUNCT__ 3318 #define __FUNCT__ "SNESGetConvergedReason" 3319 /*@ 3320 SNESGetConvergedReason - Gets the reason the SNES iteration was stopped. 3321 3322 Not Collective 3323 3324 Input Parameter: 3325 . snes - the SNES context 3326 3327 Output Parameter: 3328 . reason - negative value indicates diverged, positive value converged, see SNESConvergedReason or the 3329 manual pages for the individual convergence tests for complete lists 3330 3331 Level: intermediate 3332 3333 Notes: Can only be called after the call the SNESSolve() is complete. 3334 3335 .keywords: SNES, nonlinear, set, convergence, test 3336 3337 .seealso: SNESSetConvergenceTest(), SNESConvergedReason 3338 @*/ 3339 PetscErrorCode SNESGetConvergedReason(SNES snes,SNESConvergedReason *reason) 3340 { 3341 PetscFunctionBegin; 3342 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3343 PetscValidPointer(reason,2); 3344 *reason = snes->reason; 3345 PetscFunctionReturn(0); 3346 } 3347 3348 #undef __FUNCT__ 3349 #define __FUNCT__ "SNESSetConvergenceHistory" 3350 /*@ 3351 SNESSetConvergenceHistory - Sets the array used to hold the convergence history. 3352 3353 Logically Collective on SNES 3354 3355 Input Parameters: 3356 + snes - iterative context obtained from SNESCreate() 3357 . a - array to hold history, this array will contain the function norms computed at each step 3358 . its - integer array holds the number of linear iterations for each solve. 3359 . na - size of a and its 3360 - reset - PETSC_TRUE indicates each new nonlinear solve resets the history counter to zero, 3361 else it continues storing new values for new nonlinear solves after the old ones 3362 3363 Notes: 3364 If 'a' and 'its' are NULL then space is allocated for the history. If 'na' PETSC_DECIDE or PETSC_DEFAULT then a 3365 default array of length 10000 is allocated. 3366 3367 This routine is useful, e.g., when running a code for purposes 3368 of accurate performance monitoring, when no I/O should be done 3369 during the section of code that is being timed. 3370 3371 Level: intermediate 3372 3373 .keywords: SNES, set, convergence, history 3374 3375 .seealso: SNESGetConvergenceHistory() 3376 3377 @*/ 3378 PetscErrorCode SNESSetConvergenceHistory(SNES snes,PetscReal a[],PetscInt its[],PetscInt na,PetscBool reset) 3379 { 3380 PetscErrorCode ierr; 3381 3382 PetscFunctionBegin; 3383 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3384 if (na) PetscValidScalarPointer(a,2); 3385 if (its) PetscValidIntPointer(its,3); 3386 if (na == PETSC_DECIDE || na == PETSC_DEFAULT || !a) { 3387 if (na == PETSC_DECIDE || na == PETSC_DEFAULT) na = 1000; 3388 ierr = PetscMalloc(na*sizeof(PetscReal),&a);CHKERRQ(ierr); 3389 ierr = PetscMalloc(na*sizeof(PetscInt),&its);CHKERRQ(ierr); 3390 3391 snes->conv_malloc = PETSC_TRUE; 3392 } 3393 snes->conv_hist = a; 3394 snes->conv_hist_its = its; 3395 snes->conv_hist_max = na; 3396 snes->conv_hist_len = 0; 3397 snes->conv_hist_reset = reset; 3398 PetscFunctionReturn(0); 3399 } 3400 3401 #if defined(PETSC_HAVE_MATLAB_ENGINE) 3402 #include <engine.h> /* MATLAB include file */ 3403 #include <mex.h> /* MATLAB include file */ 3404 EXTERN_C_BEGIN 3405 #undef __FUNCT__ 3406 #define __FUNCT__ "SNESGetConvergenceHistoryMatlab" 3407 mxArray *SNESGetConvergenceHistoryMatlab(SNES snes) 3408 { 3409 mxArray *mat; 3410 PetscInt i; 3411 PetscReal *ar; 3412 3413 PetscFunctionBegin; 3414 mat = mxCreateDoubleMatrix(snes->conv_hist_len,1,mxREAL); 3415 ar = (PetscReal*) mxGetData(mat); 3416 for (i=0; i<snes->conv_hist_len; i++) ar[i] = snes->conv_hist[i]; 3417 PetscFunctionReturn(mat); 3418 } 3419 EXTERN_C_END 3420 #endif 3421 3422 #undef __FUNCT__ 3423 #define __FUNCT__ "SNESGetConvergenceHistory" 3424 /*@C 3425 SNESGetConvergenceHistory - Gets the array used to hold the convergence history. 3426 3427 Not Collective 3428 3429 Input Parameter: 3430 . snes - iterative context obtained from SNESCreate() 3431 3432 Output Parameters: 3433 . a - array to hold history 3434 . its - integer array holds the number of linear iterations (or 3435 negative if not converged) for each solve. 3436 - na - size of a and its 3437 3438 Notes: 3439 The calling sequence for this routine in Fortran is 3440 $ call SNESGetConvergenceHistory(SNES snes, integer na, integer ierr) 3441 3442 This routine is useful, e.g., when running a code for purposes 3443 of accurate performance monitoring, when no I/O should be done 3444 during the section of code that is being timed. 3445 3446 Level: intermediate 3447 3448 .keywords: SNES, get, convergence, history 3449 3450 .seealso: SNESSetConvergencHistory() 3451 3452 @*/ 3453 PetscErrorCode SNESGetConvergenceHistory(SNES snes,PetscReal *a[],PetscInt *its[],PetscInt *na) 3454 { 3455 PetscFunctionBegin; 3456 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3457 if (a) *a = snes->conv_hist; 3458 if (its) *its = snes->conv_hist_its; 3459 if (na) *na = snes->conv_hist_len; 3460 PetscFunctionReturn(0); 3461 } 3462 3463 #undef __FUNCT__ 3464 #define __FUNCT__ "SNESSetUpdate" 3465 /*@C 3466 SNESSetUpdate - Sets the general-purpose update function called 3467 at the beginning of every iteration of the nonlinear solve. Specifically 3468 it is called just before the Jacobian is "evaluated". 3469 3470 Logically Collective on SNES 3471 3472 Input Parameters: 3473 . snes - The nonlinear solver context 3474 . func - The function 3475 3476 Calling sequence of func: 3477 . func (SNES snes, PetscInt step); 3478 3479 . step - The current step of the iteration 3480 3481 Level: advanced 3482 3483 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() 3484 This is not used by most users. 3485 3486 .keywords: SNES, update 3487 3488 .seealso SNESDefaultUpdate(), SNESSetJacobian(), SNESSolve() 3489 @*/ 3490 PetscErrorCode SNESSetUpdate(SNES snes, PetscErrorCode (*func)(SNES, PetscInt)) 3491 { 3492 PetscFunctionBegin; 3493 PetscValidHeaderSpecific(snes, SNES_CLASSID,1); 3494 snes->ops->update = func; 3495 PetscFunctionReturn(0); 3496 } 3497 3498 #undef __FUNCT__ 3499 #define __FUNCT__ "SNESDefaultUpdate" 3500 /*@ 3501 SNESDefaultUpdate - The default update function which does nothing. 3502 3503 Not collective 3504 3505 Input Parameters: 3506 . snes - The nonlinear solver context 3507 . step - The current step of the iteration 3508 3509 Level: intermediate 3510 3511 .keywords: SNES, update 3512 .seealso SNESSetUpdate(), SNESDefaultRhsBC(), SNESDefaultShortolutionBC() 3513 @*/ 3514 PetscErrorCode SNESDefaultUpdate(SNES snes, PetscInt step) 3515 { 3516 PetscFunctionBegin; 3517 PetscFunctionReturn(0); 3518 } 3519 3520 #undef __FUNCT__ 3521 #define __FUNCT__ "SNESScaleStep_Private" 3522 /* 3523 SNESScaleStep_Private - Scales a step so that its length is less than the 3524 positive parameter delta. 3525 3526 Input Parameters: 3527 + snes - the SNES context 3528 . y - approximate solution of linear system 3529 . fnorm - 2-norm of current function 3530 - delta - trust region size 3531 3532 Output Parameters: 3533 + gpnorm - predicted function norm at the new point, assuming local 3534 linearization. The value is zero if the step lies within the trust 3535 region, and exceeds zero otherwise. 3536 - ynorm - 2-norm of the step 3537 3538 Note: 3539 For non-trust region methods such as SNESNEWTONLS, the parameter delta 3540 is set to be the maximum allowable step size. 3541 3542 .keywords: SNES, nonlinear, scale, step 3543 */ 3544 PetscErrorCode SNESScaleStep_Private(SNES snes,Vec y,PetscReal *fnorm,PetscReal *delta,PetscReal *gpnorm,PetscReal *ynorm) 3545 { 3546 PetscReal nrm; 3547 PetscScalar cnorm; 3548 PetscErrorCode ierr; 3549 3550 PetscFunctionBegin; 3551 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3552 PetscValidHeaderSpecific(y,VEC_CLASSID,2); 3553 PetscCheckSameComm(snes,1,y,2); 3554 3555 ierr = VecNorm(y,NORM_2,&nrm);CHKERRQ(ierr); 3556 if (nrm > *delta) { 3557 nrm = *delta/nrm; 3558 *gpnorm = (1.0 - nrm)*(*fnorm); 3559 cnorm = nrm; 3560 ierr = VecScale(y,cnorm);CHKERRQ(ierr); 3561 *ynorm = *delta; 3562 } else { 3563 *gpnorm = 0.0; 3564 *ynorm = nrm; 3565 } 3566 PetscFunctionReturn(0); 3567 } 3568 3569 #undef __FUNCT__ 3570 #define __FUNCT__ "SNESSolve" 3571 /*@C 3572 SNESSolve - Solves a nonlinear system F(x) = b. 3573 Call SNESSolve() after calling SNESCreate() and optional routines of the form SNESSetXXX(). 3574 3575 Collective on SNES 3576 3577 Input Parameters: 3578 + snes - the SNES context 3579 . b - the constant part of the equation F(x) = b, or NULL to use zero. 3580 - x - the solution vector. 3581 3582 Notes: 3583 The user should initialize the vector,x, with the initial guess 3584 for the nonlinear solve prior to calling SNESSolve. In particular, 3585 to employ an initial guess of zero, the user should explicitly set 3586 this vector to zero by calling VecSet(). 3587 3588 Level: beginner 3589 3590 .keywords: SNES, nonlinear, solve 3591 3592 .seealso: SNESCreate(), SNESDestroy(), SNESSetFunction(), SNESSetJacobian(), SNESSetGridSequence(), SNESGetSolution() 3593 @*/ 3594 PetscErrorCode SNESSolve(SNES snes,Vec b,Vec x) 3595 { 3596 PetscErrorCode ierr; 3597 PetscBool flg; 3598 PetscViewer viewer; 3599 PetscInt grid; 3600 Vec xcreated = NULL; 3601 DM dm; 3602 PetscViewerFormat format; 3603 3604 PetscFunctionBegin; 3605 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3606 if (x) PetscValidHeaderSpecific(x,VEC_CLASSID,3); 3607 if (x) PetscCheckSameComm(snes,1,x,3); 3608 if (b) PetscValidHeaderSpecific(b,VEC_CLASSID,2); 3609 if (b) PetscCheckSameComm(snes,1,b,2); 3610 3611 if (!x) { 3612 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 3613 ierr = DMCreateGlobalVector(dm,&xcreated);CHKERRQ(ierr); 3614 x = xcreated; 3615 } 3616 3617 for (grid=0; grid<snes->gridsequence; grid++) {ierr = PetscViewerASCIIPushTab(PETSC_VIEWER_STDOUT_(((PetscObject)snes)->comm));CHKERRQ(ierr);} 3618 for (grid=0; grid<snes->gridsequence+1; grid++) { 3619 3620 /* set solution vector */ 3621 if (!grid) {ierr = PetscObjectReference((PetscObject)x);CHKERRQ(ierr);} 3622 ierr = VecDestroy(&snes->vec_sol);CHKERRQ(ierr); 3623 snes->vec_sol = x; 3624 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 3625 3626 /* set affine vector if provided */ 3627 if (b) { ierr = PetscObjectReference((PetscObject)b);CHKERRQ(ierr); } 3628 ierr = VecDestroy(&snes->vec_rhs);CHKERRQ(ierr); 3629 snes->vec_rhs = b; 3630 3631 ierr = SNESSetUp(snes);CHKERRQ(ierr); 3632 3633 if (!grid) { 3634 if (snes->ops->computeinitialguess) { 3635 ierr = (*snes->ops->computeinitialguess)(snes,snes->vec_sol,snes->initialguessP);CHKERRQ(ierr); 3636 } 3637 } 3638 3639 if (snes->conv_hist_reset) snes->conv_hist_len = 0; 3640 snes->nfuncs = 0; snes->linear_its = 0; snes->numFailures = 0; 3641 3642 ierr = PetscLogEventBegin(SNES_Solve,snes,0,0,0);CHKERRQ(ierr); 3643 ierr = (*snes->ops->solve)(snes);CHKERRQ(ierr); 3644 ierr = PetscLogEventEnd(SNES_Solve,snes,0,0,0);CHKERRQ(ierr); 3645 if (snes->domainerror) { 3646 snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; 3647 snes->domainerror = PETSC_FALSE; 3648 } 3649 if (!snes->reason) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Internal error, solver returned without setting converged reason"); 3650 3651 flg = PETSC_FALSE; 3652 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_test_local_min",&flg,NULL);CHKERRQ(ierr); 3653 if (flg && !PetscPreLoadingOn) { ierr = SNESTestLocalMin(snes);CHKERRQ(ierr); } 3654 if (snes->printreason) { 3655 ierr = PetscViewerASCIIAddTab(PETSC_VIEWER_STDOUT_(((PetscObject)snes)->comm),((PetscObject)snes)->tablevel);CHKERRQ(ierr); 3656 if (snes->reason > 0) { 3657 ierr = PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_(((PetscObject)snes)->comm),"Nonlinear solve converged due to %s iterations %D\n",SNESConvergedReasons[snes->reason],snes->iter);CHKERRQ(ierr); 3658 } else { 3659 ierr = PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_(((PetscObject)snes)->comm),"Nonlinear solve did not converge due to %s iterations %D\n",SNESConvergedReasons[snes->reason],snes->iter);CHKERRQ(ierr); 3660 } 3661 ierr = PetscViewerASCIISubtractTab(PETSC_VIEWER_STDOUT_(((PetscObject)snes)->comm),((PetscObject)snes)->tablevel);CHKERRQ(ierr); 3662 } 3663 3664 if (snes->errorifnotconverged && snes->reason < 0) SETERRQ(((PetscObject)snes)->comm,PETSC_ERR_NOT_CONVERGED,"SNESSolve has not converged"); 3665 if (grid < snes->gridsequence) { 3666 DM fine; 3667 Vec xnew; 3668 Mat interp; 3669 3670 ierr = DMRefine(snes->dm,((PetscObject)snes)->comm,&fine);CHKERRQ(ierr); 3671 if (!fine) SETERRQ(((PetscObject)snes)->comm,PETSC_ERR_ARG_INCOMP,"DMRefine() did not perform any refinement, cannot continue grid sequencing"); 3672 ierr = DMCreateInterpolation(snes->dm,fine,&interp,NULL);CHKERRQ(ierr); 3673 ierr = DMCreateGlobalVector(fine,&xnew);CHKERRQ(ierr); 3674 ierr = MatInterpolate(interp,x,xnew);CHKERRQ(ierr); 3675 ierr = DMInterpolate(snes->dm,interp,fine);CHKERRQ(ierr); 3676 ierr = MatDestroy(&interp);CHKERRQ(ierr); 3677 x = xnew; 3678 3679 ierr = SNESReset(snes);CHKERRQ(ierr); 3680 ierr = SNESSetDM(snes,fine);CHKERRQ(ierr); 3681 ierr = DMDestroy(&fine);CHKERRQ(ierr); 3682 ierr = PetscViewerASCIIPopTab(PETSC_VIEWER_STDOUT_(((PetscObject)snes)->comm));CHKERRQ(ierr); 3683 } 3684 } 3685 /* monitoring and viewing */ 3686 ierr = PetscOptionsGetViewer(((PetscObject)snes)->comm,((PetscObject)snes)->prefix,"-snes_view",&viewer,&format,&flg);CHKERRQ(ierr); 3687 if (flg && !PetscPreLoadingOn) { 3688 ierr = PetscViewerPushFormat(viewer,format);CHKERRQ(ierr); 3689 ierr = SNESView(snes,viewer);CHKERRQ(ierr); 3690 ierr = PetscViewerPopFormat(viewer);CHKERRQ(ierr); 3691 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 3692 } 3693 ierr = VecViewFromOptions(snes->vec_sol,"-snes_view_solution");CHKERRQ(ierr); 3694 3695 ierr = VecDestroy(&xcreated);CHKERRQ(ierr); 3696 PetscFunctionReturn(0); 3697 } 3698 3699 /* --------- Internal routines for SNES Package --------- */ 3700 3701 #undef __FUNCT__ 3702 #define __FUNCT__ "SNESSetType" 3703 /*@C 3704 SNESSetType - Sets the method for the nonlinear solver. 3705 3706 Collective on SNES 3707 3708 Input Parameters: 3709 + snes - the SNES context 3710 - type - a known method 3711 3712 Options Database Key: 3713 . -snes_type <type> - Sets the method; use -help for a list 3714 of available methods (for instance, newtonls or newtontr) 3715 3716 Notes: 3717 See "petsc/include/petscsnes.h" for available methods (for instance) 3718 + SNESNEWTONLS - Newton's method with line search 3719 (systems of nonlinear equations) 3720 . SNESNEWTONTR - Newton's method with trust region 3721 (systems of nonlinear equations) 3722 3723 Normally, it is best to use the SNESSetFromOptions() command and then 3724 set the SNES solver type from the options database rather than by using 3725 this routine. Using the options database provides the user with 3726 maximum flexibility in evaluating the many nonlinear solvers. 3727 The SNESSetType() routine is provided for those situations where it 3728 is necessary to set the nonlinear solver independently of the command 3729 line or options database. This might be the case, for example, when 3730 the choice of solver changes during the execution of the program, 3731 and the user's application is taking responsibility for choosing the 3732 appropriate method. 3733 3734 Level: intermediate 3735 3736 .keywords: SNES, set, type 3737 3738 .seealso: SNESType, SNESCreate() 3739 3740 @*/ 3741 PetscErrorCode SNESSetType(SNES snes,SNESType type) 3742 { 3743 PetscErrorCode ierr,(*r)(SNES); 3744 PetscBool match; 3745 3746 PetscFunctionBegin; 3747 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3748 PetscValidCharPointer(type,2); 3749 3750 ierr = PetscObjectTypeCompare((PetscObject)snes,type,&match);CHKERRQ(ierr); 3751 if (match) PetscFunctionReturn(0); 3752 3753 ierr = PetscFunctionListFind(((PetscObject)snes)->comm,SNESList,type,PETSC_TRUE,(void (**)(void)) &r);CHKERRQ(ierr); 3754 if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unable to find requested SNES type %s",type); 3755 /* Destroy the previous private SNES context */ 3756 if (snes->ops->destroy) { 3757 ierr = (*(snes)->ops->destroy)(snes);CHKERRQ(ierr); 3758 snes->ops->destroy = NULL; 3759 } 3760 /* Reinitialize function pointers in SNESOps structure */ 3761 snes->ops->setup = 0; 3762 snes->ops->solve = 0; 3763 snes->ops->view = 0; 3764 snes->ops->setfromoptions = 0; 3765 snes->ops->destroy = 0; 3766 /* Call the SNESCreate_XXX routine for this particular Nonlinear solver */ 3767 snes->setupcalled = PETSC_FALSE; 3768 3769 ierr = PetscObjectChangeTypeName((PetscObject)snes,type);CHKERRQ(ierr); 3770 ierr = (*r)(snes);CHKERRQ(ierr); 3771 #if defined(PETSC_HAVE_AMS) 3772 if (PetscAMSPublishAll) { 3773 ierr = PetscObjectAMSPublish((PetscObject)snes);CHKERRQ(ierr); 3774 } 3775 #endif 3776 PetscFunctionReturn(0); 3777 } 3778 3779 /* --------------------------------------------------------------------- */ 3780 #undef __FUNCT__ 3781 #define __FUNCT__ "SNESRegisterDestroy" 3782 /*@ 3783 SNESRegisterDestroy - Frees the list of nonlinear solvers that were 3784 registered by SNESRegisterDynamic(). 3785 3786 Not Collective 3787 3788 Level: advanced 3789 3790 .keywords: SNES, nonlinear, register, destroy 3791 3792 .seealso: SNESRegisterAll() 3793 @*/ 3794 PetscErrorCode SNESRegisterDestroy(void) 3795 { 3796 PetscErrorCode ierr; 3797 3798 PetscFunctionBegin; 3799 ierr = PetscFunctionListDestroy(&SNESList);CHKERRQ(ierr); 3800 3801 SNESRegisterAllCalled = PETSC_FALSE; 3802 PetscFunctionReturn(0); 3803 } 3804 3805 #undef __FUNCT__ 3806 #define __FUNCT__ "SNESGetType" 3807 /*@C 3808 SNESGetType - Gets the SNES method type and name (as a string). 3809 3810 Not Collective 3811 3812 Input Parameter: 3813 . snes - nonlinear solver context 3814 3815 Output Parameter: 3816 . type - SNES method (a character string) 3817 3818 Level: intermediate 3819 3820 .keywords: SNES, nonlinear, get, type, name 3821 @*/ 3822 PetscErrorCode SNESGetType(SNES snes,SNESType *type) 3823 { 3824 PetscFunctionBegin; 3825 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3826 PetscValidPointer(type,2); 3827 *type = ((PetscObject)snes)->type_name; 3828 PetscFunctionReturn(0); 3829 } 3830 3831 #undef __FUNCT__ 3832 #define __FUNCT__ "SNESGetSolution" 3833 /*@ 3834 SNESGetSolution - Returns the vector where the approximate solution is 3835 stored. This is the fine grid solution when using SNESSetGridSequence(). 3836 3837 Not Collective, but Vec is parallel if SNES is parallel 3838 3839 Input Parameter: 3840 . snes - the SNES context 3841 3842 Output Parameter: 3843 . x - the solution 3844 3845 Level: intermediate 3846 3847 .keywords: SNES, nonlinear, get, solution 3848 3849 .seealso: SNESGetSolutionUpdate(), SNESGetFunction() 3850 @*/ 3851 PetscErrorCode SNESGetSolution(SNES snes,Vec *x) 3852 { 3853 PetscFunctionBegin; 3854 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3855 PetscValidPointer(x,2); 3856 *x = snes->vec_sol; 3857 PetscFunctionReturn(0); 3858 } 3859 3860 #undef __FUNCT__ 3861 #define __FUNCT__ "SNESGetSolutionUpdate" 3862 /*@ 3863 SNESGetSolutionUpdate - Returns the vector where the solution update is 3864 stored. 3865 3866 Not Collective, but Vec is parallel if SNES is parallel 3867 3868 Input Parameter: 3869 . snes - the SNES context 3870 3871 Output Parameter: 3872 . x - the solution update 3873 3874 Level: advanced 3875 3876 .keywords: SNES, nonlinear, get, solution, update 3877 3878 .seealso: SNESGetSolution(), SNESGetFunction() 3879 @*/ 3880 PetscErrorCode SNESGetSolutionUpdate(SNES snes,Vec *x) 3881 { 3882 PetscFunctionBegin; 3883 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3884 PetscValidPointer(x,2); 3885 *x = snes->vec_sol_update; 3886 PetscFunctionReturn(0); 3887 } 3888 3889 #undef __FUNCT__ 3890 #define __FUNCT__ "SNESGetFunction" 3891 /*@C 3892 SNESGetFunction - Returns the vector where the function is stored. 3893 3894 Not Collective, but Vec is parallel if SNES is parallel. Collective if Vec is requested, but has not been created yet. 3895 3896 Input Parameter: 3897 . snes - the SNES context 3898 3899 Output Parameter: 3900 + r - the vector that is used to store residuals (or NULL if you don't want it) 3901 . SNESFunction- the function (or NULL if you don't want it) 3902 - ctx - the function context (or NULL if you don't want it) 3903 3904 Level: advanced 3905 3906 .keywords: SNES, nonlinear, get, function 3907 3908 .seealso: SNESSetFunction(), SNESGetSolution(), SNESFunction 3909 @*/ 3910 PetscErrorCode SNESGetFunction(SNES snes,Vec *r,PetscErrorCode (**SNESFunction)(SNES,Vec,Vec,void*),void **ctx) 3911 { 3912 PetscErrorCode ierr; 3913 DM dm; 3914 3915 PetscFunctionBegin; 3916 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3917 if (r) { 3918 if (!snes->vec_func) { 3919 if (snes->vec_rhs) { 3920 ierr = VecDuplicate(snes->vec_rhs,&snes->vec_func);CHKERRQ(ierr); 3921 } else if (snes->vec_sol) { 3922 ierr = VecDuplicate(snes->vec_sol,&snes->vec_func);CHKERRQ(ierr); 3923 } else if (snes->dm) { 3924 ierr = DMCreateGlobalVector(snes->dm,&snes->vec_func);CHKERRQ(ierr); 3925 } 3926 } 3927 *r = snes->vec_func; 3928 } 3929 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 3930 ierr = DMSNESGetFunction(dm,SNESFunction,ctx);CHKERRQ(ierr); 3931 PetscFunctionReturn(0); 3932 } 3933 3934 /*@C 3935 SNESGetGS - Returns the GS function and context. 3936 3937 Input Parameter: 3938 . snes - the SNES context 3939 3940 Output Parameter: 3941 + SNESGSFunction - the function (or NULL) 3942 - ctx - the function context (or NULL) 3943 3944 Level: advanced 3945 3946 .keywords: SNES, nonlinear, get, function 3947 3948 .seealso: SNESSetGS(), SNESGetFunction() 3949 @*/ 3950 3951 #undef __FUNCT__ 3952 #define __FUNCT__ "SNESGetGS" 3953 PetscErrorCode SNESGetGS (SNES snes, PetscErrorCode (**SNESGSFunction)(SNES, Vec, Vec, void*), void ** ctx) 3954 { 3955 PetscErrorCode ierr; 3956 DM dm; 3957 3958 PetscFunctionBegin; 3959 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3960 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 3961 ierr = DMSNESGetGS(dm,SNESGSFunction,ctx);CHKERRQ(ierr); 3962 PetscFunctionReturn(0); 3963 } 3964 3965 #undef __FUNCT__ 3966 #define __FUNCT__ "SNESSetOptionsPrefix" 3967 /*@C 3968 SNESSetOptionsPrefix - Sets the prefix used for searching for all 3969 SNES options in the database. 3970 3971 Logically Collective on SNES 3972 3973 Input Parameter: 3974 + snes - the SNES context 3975 - prefix - the prefix to prepend to all option names 3976 3977 Notes: 3978 A hyphen (-) must NOT be given at the beginning of the prefix name. 3979 The first character of all runtime options is AUTOMATICALLY the hyphen. 3980 3981 Level: advanced 3982 3983 .keywords: SNES, set, options, prefix, database 3984 3985 .seealso: SNESSetFromOptions() 3986 @*/ 3987 PetscErrorCode SNESSetOptionsPrefix(SNES snes,const char prefix[]) 3988 { 3989 PetscErrorCode ierr; 3990 3991 PetscFunctionBegin; 3992 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3993 ierr = PetscObjectSetOptionsPrefix((PetscObject)snes,prefix);CHKERRQ(ierr); 3994 if (!snes->ksp) {ierr = SNESGetKSP(snes,&snes->ksp);CHKERRQ(ierr);} 3995 if (snes->linesearch) { 3996 ierr = SNESGetSNESLineSearch(snes,&snes->linesearch);CHKERRQ(ierr); 3997 ierr = PetscObjectSetOptionsPrefix((PetscObject)snes->linesearch,prefix);CHKERRQ(ierr); 3998 } 3999 ierr = KSPSetOptionsPrefix(snes->ksp,prefix);CHKERRQ(ierr); 4000 PetscFunctionReturn(0); 4001 } 4002 4003 #undef __FUNCT__ 4004 #define __FUNCT__ "SNESAppendOptionsPrefix" 4005 /*@C 4006 SNESAppendOptionsPrefix - Appends to the prefix used for searching for all 4007 SNES options in the database. 4008 4009 Logically Collective on SNES 4010 4011 Input Parameters: 4012 + snes - the SNES context 4013 - prefix - the prefix to prepend to all option names 4014 4015 Notes: 4016 A hyphen (-) must NOT be given at the beginning of the prefix name. 4017 The first character of all runtime options is AUTOMATICALLY the hyphen. 4018 4019 Level: advanced 4020 4021 .keywords: SNES, append, options, prefix, database 4022 4023 .seealso: SNESGetOptionsPrefix() 4024 @*/ 4025 PetscErrorCode SNESAppendOptionsPrefix(SNES snes,const char prefix[]) 4026 { 4027 PetscErrorCode ierr; 4028 4029 PetscFunctionBegin; 4030 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4031 ierr = PetscObjectAppendOptionsPrefix((PetscObject)snes,prefix);CHKERRQ(ierr); 4032 if (!snes->ksp) {ierr = SNESGetKSP(snes,&snes->ksp);CHKERRQ(ierr);} 4033 if (snes->linesearch) { 4034 ierr = SNESGetSNESLineSearch(snes,&snes->linesearch);CHKERRQ(ierr); 4035 ierr = PetscObjectAppendOptionsPrefix((PetscObject)snes->linesearch,prefix);CHKERRQ(ierr); 4036 } 4037 ierr = KSPAppendOptionsPrefix(snes->ksp,prefix);CHKERRQ(ierr); 4038 PetscFunctionReturn(0); 4039 } 4040 4041 #undef __FUNCT__ 4042 #define __FUNCT__ "SNESGetOptionsPrefix" 4043 /*@C 4044 SNESGetOptionsPrefix - Sets the prefix used for searching for all 4045 SNES options in the database. 4046 4047 Not Collective 4048 4049 Input Parameter: 4050 . snes - the SNES context 4051 4052 Output Parameter: 4053 . prefix - pointer to the prefix string used 4054 4055 Notes: On the fortran side, the user should pass in a string 'prefix' of 4056 sufficient length to hold the prefix. 4057 4058 Level: advanced 4059 4060 .keywords: SNES, get, options, prefix, database 4061 4062 .seealso: SNESAppendOptionsPrefix() 4063 @*/ 4064 PetscErrorCode SNESGetOptionsPrefix(SNES snes,const char *prefix[]) 4065 { 4066 PetscErrorCode ierr; 4067 4068 PetscFunctionBegin; 4069 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4070 ierr = PetscObjectGetOptionsPrefix((PetscObject)snes,prefix);CHKERRQ(ierr); 4071 PetscFunctionReturn(0); 4072 } 4073 4074 4075 #undef __FUNCT__ 4076 #define __FUNCT__ "SNESRegister" 4077 /*@C 4078 SNESRegister - See SNESRegisterDynamic() 4079 4080 Level: advanced 4081 @*/ 4082 PetscErrorCode SNESRegister(const char sname[],const char path[],const char name[],PetscErrorCode (*function)(SNES)) 4083 { 4084 char fullname[PETSC_MAX_PATH_LEN]; 4085 PetscErrorCode ierr; 4086 4087 PetscFunctionBegin; 4088 ierr = PetscFunctionListConcat(path,name,fullname);CHKERRQ(ierr); 4089 ierr = PetscFunctionListAdd(PETSC_COMM_WORLD,&SNESList,sname,fullname,(void (*)(void))function);CHKERRQ(ierr); 4090 PetscFunctionReturn(0); 4091 } 4092 4093 #undef __FUNCT__ 4094 #define __FUNCT__ "SNESTestLocalMin" 4095 PetscErrorCode SNESTestLocalMin(SNES snes) 4096 { 4097 PetscErrorCode ierr; 4098 PetscInt N,i,j; 4099 Vec u,uh,fh; 4100 PetscScalar value; 4101 PetscReal norm; 4102 4103 PetscFunctionBegin; 4104 ierr = SNESGetSolution(snes,&u);CHKERRQ(ierr); 4105 ierr = VecDuplicate(u,&uh);CHKERRQ(ierr); 4106 ierr = VecDuplicate(u,&fh);CHKERRQ(ierr); 4107 4108 /* currently only works for sequential */ 4109 ierr = PetscPrintf(PETSC_COMM_WORLD,"Testing FormFunction() for local min\n");CHKERRQ(ierr); 4110 ierr = VecGetSize(u,&N);CHKERRQ(ierr); 4111 for (i=0; i<N; i++) { 4112 ierr = VecCopy(u,uh);CHKERRQ(ierr); 4113 ierr = PetscPrintf(PETSC_COMM_WORLD,"i = %D\n",i);CHKERRQ(ierr); 4114 for (j=-10; j<11; j++) { 4115 value = PetscSign(j)*exp(PetscAbs(j)-10.0); 4116 ierr = VecSetValue(uh,i,value,ADD_VALUES);CHKERRQ(ierr); 4117 ierr = SNESComputeFunction(snes,uh,fh);CHKERRQ(ierr); 4118 ierr = VecNorm(fh,NORM_2,&norm);CHKERRQ(ierr); 4119 ierr = PetscPrintf(PETSC_COMM_WORLD," j norm %D %18.16e\n",j,norm);CHKERRQ(ierr); 4120 value = -value; 4121 ierr = VecSetValue(uh,i,value,ADD_VALUES);CHKERRQ(ierr); 4122 } 4123 } 4124 ierr = VecDestroy(&uh);CHKERRQ(ierr); 4125 ierr = VecDestroy(&fh);CHKERRQ(ierr); 4126 PetscFunctionReturn(0); 4127 } 4128 4129 #undef __FUNCT__ 4130 #define __FUNCT__ "SNESKSPSetUseEW" 4131 /*@ 4132 SNESKSPSetUseEW - Sets SNES use Eisenstat-Walker method for 4133 computing relative tolerance for linear solvers within an inexact 4134 Newton method. 4135 4136 Logically Collective on SNES 4137 4138 Input Parameters: 4139 + snes - SNES context 4140 - flag - PETSC_TRUE or PETSC_FALSE 4141 4142 Options Database: 4143 + -snes_ksp_ew - use Eisenstat-Walker method for determining linear system convergence 4144 . -snes_ksp_ew_version ver - version of Eisenstat-Walker method 4145 . -snes_ksp_ew_rtol0 <rtol0> - Sets rtol0 4146 . -snes_ksp_ew_rtolmax <rtolmax> - Sets rtolmax 4147 . -snes_ksp_ew_gamma <gamma> - Sets gamma 4148 . -snes_ksp_ew_alpha <alpha> - Sets alpha 4149 . -snes_ksp_ew_alpha2 <alpha2> - Sets alpha2 4150 - -snes_ksp_ew_threshold <threshold> - Sets threshold 4151 4152 Notes: 4153 Currently, the default is to use a constant relative tolerance for 4154 the inner linear solvers. Alternatively, one can use the 4155 Eisenstat-Walker method, where the relative convergence tolerance 4156 is reset at each Newton iteration according progress of the nonlinear 4157 solver. 4158 4159 Level: advanced 4160 4161 Reference: 4162 S. C. Eisenstat and H. F. Walker, "Choosing the forcing terms in an 4163 inexact Newton method", SISC 17 (1), pp.16-32, 1996. 4164 4165 .keywords: SNES, KSP, Eisenstat, Walker, convergence, test, inexact, Newton 4166 4167 .seealso: SNESKSPGetUseEW(), SNESKSPGetParametersEW(), SNESKSPSetParametersEW() 4168 @*/ 4169 PetscErrorCode SNESKSPSetUseEW(SNES snes,PetscBool flag) 4170 { 4171 PetscFunctionBegin; 4172 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4173 PetscValidLogicalCollectiveBool(snes,flag,2); 4174 snes->ksp_ewconv = flag; 4175 PetscFunctionReturn(0); 4176 } 4177 4178 #undef __FUNCT__ 4179 #define __FUNCT__ "SNESKSPGetUseEW" 4180 /*@ 4181 SNESKSPGetUseEW - Gets if SNES is using Eisenstat-Walker method 4182 for computing relative tolerance for linear solvers within an 4183 inexact Newton method. 4184 4185 Not Collective 4186 4187 Input Parameter: 4188 . snes - SNES context 4189 4190 Output Parameter: 4191 . flag - PETSC_TRUE or PETSC_FALSE 4192 4193 Notes: 4194 Currently, the default is to use a constant relative tolerance for 4195 the inner linear solvers. Alternatively, one can use the 4196 Eisenstat-Walker method, where the relative convergence tolerance 4197 is reset at each Newton iteration according progress of the nonlinear 4198 solver. 4199 4200 Level: advanced 4201 4202 Reference: 4203 S. C. Eisenstat and H. F. Walker, "Choosing the forcing terms in an 4204 inexact Newton method", SISC 17 (1), pp.16-32, 1996. 4205 4206 .keywords: SNES, KSP, Eisenstat, Walker, convergence, test, inexact, Newton 4207 4208 .seealso: SNESKSPSetUseEW(), SNESKSPGetParametersEW(), SNESKSPSetParametersEW() 4209 @*/ 4210 PetscErrorCode SNESKSPGetUseEW(SNES snes, PetscBool *flag) 4211 { 4212 PetscFunctionBegin; 4213 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4214 PetscValidPointer(flag,2); 4215 *flag = snes->ksp_ewconv; 4216 PetscFunctionReturn(0); 4217 } 4218 4219 #undef __FUNCT__ 4220 #define __FUNCT__ "SNESKSPSetParametersEW" 4221 /*@ 4222 SNESKSPSetParametersEW - Sets parameters for Eisenstat-Walker 4223 convergence criteria for the linear solvers within an inexact 4224 Newton method. 4225 4226 Logically Collective on SNES 4227 4228 Input Parameters: 4229 + snes - SNES context 4230 . version - version 1, 2 (default is 2) or 3 4231 . rtol_0 - initial relative tolerance (0 <= rtol_0 < 1) 4232 . rtol_max - maximum relative tolerance (0 <= rtol_max < 1) 4233 . gamma - multiplicative factor for version 2 rtol computation 4234 (0 <= gamma2 <= 1) 4235 . alpha - power for version 2 rtol computation (1 < alpha <= 2) 4236 . alpha2 - power for safeguard 4237 - threshold - threshold for imposing safeguard (0 < threshold < 1) 4238 4239 Note: 4240 Version 3 was contributed by Luis Chacon, June 2006. 4241 4242 Use PETSC_DEFAULT to retain the default for any of the parameters. 4243 4244 Level: advanced 4245 4246 Reference: 4247 S. C. Eisenstat and H. F. Walker, "Choosing the forcing terms in an 4248 inexact Newton method", Utah State University Math. Stat. Dept. Res. 4249 Report 6/94/75, June, 1994, to appear in SIAM J. Sci. Comput. 4250 4251 .keywords: SNES, KSP, Eisenstat, Walker, set, parameters 4252 4253 .seealso: SNESKSPSetUseEW(), SNESKSPGetUseEW(), SNESKSPGetParametersEW() 4254 @*/ 4255 PetscErrorCode SNESKSPSetParametersEW(SNES snes,PetscInt version,PetscReal rtol_0,PetscReal rtol_max,PetscReal gamma,PetscReal alpha,PetscReal alpha2,PetscReal threshold) 4256 { 4257 SNESKSPEW *kctx; 4258 4259 PetscFunctionBegin; 4260 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4261 kctx = (SNESKSPEW*)snes->kspconvctx; 4262 if (!kctx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"No Eisenstat-Walker context existing"); 4263 PetscValidLogicalCollectiveInt(snes,version,2); 4264 PetscValidLogicalCollectiveReal(snes,rtol_0,3); 4265 PetscValidLogicalCollectiveReal(snes,rtol_max,4); 4266 PetscValidLogicalCollectiveReal(snes,gamma,5); 4267 PetscValidLogicalCollectiveReal(snes,alpha,6); 4268 PetscValidLogicalCollectiveReal(snes,alpha2,7); 4269 PetscValidLogicalCollectiveReal(snes,threshold,8); 4270 4271 if (version != PETSC_DEFAULT) kctx->version = version; 4272 if (rtol_0 != PETSC_DEFAULT) kctx->rtol_0 = rtol_0; 4273 if (rtol_max != PETSC_DEFAULT) kctx->rtol_max = rtol_max; 4274 if (gamma != PETSC_DEFAULT) kctx->gamma = gamma; 4275 if (alpha != PETSC_DEFAULT) kctx->alpha = alpha; 4276 if (alpha2 != PETSC_DEFAULT) kctx->alpha2 = alpha2; 4277 if (threshold != PETSC_DEFAULT) kctx->threshold = threshold; 4278 4279 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); 4280 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); 4281 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); 4282 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); 4283 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); 4284 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); 4285 PetscFunctionReturn(0); 4286 } 4287 4288 #undef __FUNCT__ 4289 #define __FUNCT__ "SNESKSPGetParametersEW" 4290 /*@ 4291 SNESKSPGetParametersEW - Gets parameters for Eisenstat-Walker 4292 convergence criteria for the linear solvers within an inexact 4293 Newton method. 4294 4295 Not Collective 4296 4297 Input Parameters: 4298 snes - SNES context 4299 4300 Output Parameters: 4301 + version - version 1, 2 (default is 2) or 3 4302 . rtol_0 - initial relative tolerance (0 <= rtol_0 < 1) 4303 . rtol_max - maximum relative tolerance (0 <= rtol_max < 1) 4304 . gamma - multiplicative factor for version 2 rtol computation (0 <= gamma2 <= 1) 4305 . alpha - power for version 2 rtol computation (1 < alpha <= 2) 4306 . alpha2 - power for safeguard 4307 - threshold - threshold for imposing safeguard (0 < threshold < 1) 4308 4309 Level: advanced 4310 4311 .keywords: SNES, KSP, Eisenstat, Walker, get, parameters 4312 4313 .seealso: SNESKSPSetUseEW(), SNESKSPGetUseEW(), SNESKSPSetParametersEW() 4314 @*/ 4315 PetscErrorCode SNESKSPGetParametersEW(SNES snes,PetscInt *version,PetscReal *rtol_0,PetscReal *rtol_max,PetscReal *gamma,PetscReal *alpha,PetscReal *alpha2,PetscReal *threshold) 4316 { 4317 SNESKSPEW *kctx; 4318 4319 PetscFunctionBegin; 4320 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4321 kctx = (SNESKSPEW*)snes->kspconvctx; 4322 if (!kctx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"No Eisenstat-Walker context existing"); 4323 if (version) *version = kctx->version; 4324 if (rtol_0) *rtol_0 = kctx->rtol_0; 4325 if (rtol_max) *rtol_max = kctx->rtol_max; 4326 if (gamma) *gamma = kctx->gamma; 4327 if (alpha) *alpha = kctx->alpha; 4328 if (alpha2) *alpha2 = kctx->alpha2; 4329 if (threshold) *threshold = kctx->threshold; 4330 PetscFunctionReturn(0); 4331 } 4332 4333 #undef __FUNCT__ 4334 #define __FUNCT__ "SNESKSPEW_PreSolve" 4335 static PetscErrorCode SNESKSPEW_PreSolve(SNES snes, KSP ksp, Vec b, Vec x) 4336 { 4337 PetscErrorCode ierr; 4338 SNESKSPEW *kctx = (SNESKSPEW*)snes->kspconvctx; 4339 PetscReal rtol = PETSC_DEFAULT,stol; 4340 4341 PetscFunctionBegin; 4342 if (!kctx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"No Eisenstat-Walker context exists"); 4343 if (!snes->iter) rtol = kctx->rtol_0; /* first time in, so use the original user rtol */ 4344 else { 4345 if (kctx->version == 1) { 4346 rtol = (snes->norm - kctx->lresid_last)/kctx->norm_last; 4347 if (rtol < 0.0) rtol = -rtol; 4348 stol = pow(kctx->rtol_last,kctx->alpha2); 4349 if (stol > kctx->threshold) rtol = PetscMax(rtol,stol); 4350 } else if (kctx->version == 2) { 4351 rtol = kctx->gamma * pow(snes->norm/kctx->norm_last,kctx->alpha); 4352 stol = kctx->gamma * pow(kctx->rtol_last,kctx->alpha); 4353 if (stol > kctx->threshold) rtol = PetscMax(rtol,stol); 4354 } else if (kctx->version == 3) { /* contributed by Luis Chacon, June 2006. */ 4355 rtol = kctx->gamma * pow(snes->norm/kctx->norm_last,kctx->alpha); 4356 /* safeguard: avoid sharp decrease of rtol */ 4357 stol = kctx->gamma*pow(kctx->rtol_last,kctx->alpha); 4358 stol = PetscMax(rtol,stol); 4359 rtol = PetscMin(kctx->rtol_0,stol); 4360 /* safeguard: avoid oversolving */ 4361 stol = kctx->gamma*(snes->ttol)/snes->norm; 4362 stol = PetscMax(rtol,stol); 4363 rtol = PetscMin(kctx->rtol_0,stol); 4364 } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Only versions 1, 2 or 3 are supported: %D",kctx->version); 4365 } 4366 /* safeguard: avoid rtol greater than one */ 4367 rtol = PetscMin(rtol,kctx->rtol_max); 4368 ierr = KSPSetTolerances(ksp,rtol,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);CHKERRQ(ierr); 4369 ierr = PetscInfo3(snes,"iter %D, Eisenstat-Walker (version %D) KSP rtol=%G\n",snes->iter,kctx->version,rtol);CHKERRQ(ierr); 4370 PetscFunctionReturn(0); 4371 } 4372 4373 #undef __FUNCT__ 4374 #define __FUNCT__ "SNESKSPEW_PostSolve" 4375 static PetscErrorCode SNESKSPEW_PostSolve(SNES snes, KSP ksp, Vec b, Vec x) 4376 { 4377 PetscErrorCode ierr; 4378 SNESKSPEW *kctx = (SNESKSPEW*)snes->kspconvctx; 4379 PCSide pcside; 4380 Vec lres; 4381 4382 PetscFunctionBegin; 4383 if (!kctx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"No Eisenstat-Walker context exists"); 4384 ierr = KSPGetTolerances(ksp,&kctx->rtol_last,0,0,0);CHKERRQ(ierr); 4385 ierr = SNESGetFunctionNorm(snes,&kctx->norm_last);CHKERRQ(ierr); 4386 if (kctx->version == 1) { 4387 ierr = KSPGetPCSide(ksp,&pcside);CHKERRQ(ierr); 4388 if (pcside == PC_RIGHT) { /* XXX Should we also test KSP_UNPRECONDITIONED_NORM ? */ 4389 /* KSP residual is true linear residual */ 4390 ierr = KSPGetResidualNorm(ksp,&kctx->lresid_last);CHKERRQ(ierr); 4391 } else { 4392 /* KSP residual is preconditioned residual */ 4393 /* compute true linear residual norm */ 4394 ierr = VecDuplicate(b,&lres);CHKERRQ(ierr); 4395 ierr = MatMult(snes->jacobian,x,lres);CHKERRQ(ierr); 4396 ierr = VecAYPX(lres,-1.0,b);CHKERRQ(ierr); 4397 ierr = VecNorm(lres,NORM_2,&kctx->lresid_last);CHKERRQ(ierr); 4398 ierr = VecDestroy(&lres);CHKERRQ(ierr); 4399 } 4400 } 4401 PetscFunctionReturn(0); 4402 } 4403 4404 #undef __FUNCT__ 4405 #define __FUNCT__ "SNES_KSPSolve" 4406 PetscErrorCode SNES_KSPSolve(SNES snes, KSP ksp, Vec b, Vec x) 4407 { 4408 PetscErrorCode ierr; 4409 4410 PetscFunctionBegin; 4411 if (snes->ksp_ewconv) { ierr = SNESKSPEW_PreSolve(snes,ksp,b,x);CHKERRQ(ierr); } 4412 ierr = KSPSolve(ksp,b,x);CHKERRQ(ierr); 4413 if (snes->ksp_ewconv) { ierr = SNESKSPEW_PostSolve(snes,ksp,b,x);CHKERRQ(ierr); } 4414 PetscFunctionReturn(0); 4415 } 4416 4417 #include <petsc-private/dmimpl.h> 4418 #undef __FUNCT__ 4419 #define __FUNCT__ "SNESSetDM" 4420 /*@ 4421 SNESSetDM - Sets the DM that may be used by some preconditioners 4422 4423 Logically Collective on SNES 4424 4425 Input Parameters: 4426 + snes - the preconditioner context 4427 - dm - the dm 4428 4429 Level: intermediate 4430 4431 .seealso: SNESGetDM(), KSPSetDM(), KSPGetDM() 4432 @*/ 4433 PetscErrorCode SNESSetDM(SNES snes,DM dm) 4434 { 4435 PetscErrorCode ierr; 4436 KSP ksp; 4437 DMSNES sdm; 4438 4439 PetscFunctionBegin; 4440 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4441 if (dm) {ierr = PetscObjectReference((PetscObject)dm);CHKERRQ(ierr);} 4442 if (snes->dm) { /* Move the DMSNES context over to the new DM unless the new DM already has one */ 4443 if (snes->dm->dmsnes && snes->dmAuto && !dm->dmsnes) { 4444 ierr = DMCopyDMSNES(snes->dm,dm);CHKERRQ(ierr); 4445 ierr = DMGetDMSNES(snes->dm,&sdm);CHKERRQ(ierr); 4446 if (sdm->originaldm == snes->dm) sdm->originaldm = dm; /* Grant write privileges to the replacement DM */ 4447 } 4448 ierr = DMDestroy(&snes->dm);CHKERRQ(ierr); 4449 } 4450 snes->dm = dm; 4451 snes->dmAuto = PETSC_FALSE; 4452 4453 ierr = SNESGetKSP(snes,&ksp);CHKERRQ(ierr); 4454 ierr = KSPSetDM(ksp,dm);CHKERRQ(ierr); 4455 ierr = KSPSetDMActive(ksp,PETSC_FALSE);CHKERRQ(ierr); 4456 if (snes->pc) { 4457 ierr = SNESSetDM(snes->pc, snes->dm);CHKERRQ(ierr); 4458 ierr = SNESSetPCSide(snes,snes->pcside);CHKERRQ(ierr); 4459 } 4460 PetscFunctionReturn(0); 4461 } 4462 4463 #undef __FUNCT__ 4464 #define __FUNCT__ "SNESGetDM" 4465 /*@ 4466 SNESGetDM - Gets the DM that may be used by some preconditioners 4467 4468 Not Collective but DM obtained is parallel on SNES 4469 4470 Input Parameter: 4471 . snes - the preconditioner context 4472 4473 Output Parameter: 4474 . dm - the dm 4475 4476 Level: intermediate 4477 4478 .seealso: SNESSetDM(), KSPSetDM(), KSPGetDM() 4479 @*/ 4480 PetscErrorCode SNESGetDM(SNES snes,DM *dm) 4481 { 4482 PetscErrorCode ierr; 4483 4484 PetscFunctionBegin; 4485 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4486 if (!snes->dm) { 4487 ierr = DMShellCreate(((PetscObject)snes)->comm,&snes->dm);CHKERRQ(ierr); 4488 snes->dmAuto = PETSC_TRUE; 4489 } 4490 *dm = snes->dm; 4491 PetscFunctionReturn(0); 4492 } 4493 4494 #undef __FUNCT__ 4495 #define __FUNCT__ "SNESSetPC" 4496 /*@ 4497 SNESSetPC - Sets the nonlinear preconditioner to be used. 4498 4499 Collective on SNES 4500 4501 Input Parameters: 4502 + snes - iterative context obtained from SNESCreate() 4503 - pc - the preconditioner object 4504 4505 Notes: 4506 Use SNESGetPC() to retrieve the preconditioner context (for example, 4507 to configure it using the API). 4508 4509 Level: developer 4510 4511 .keywords: SNES, set, precondition 4512 .seealso: SNESGetPC() 4513 @*/ 4514 PetscErrorCode SNESSetPC(SNES snes, SNES pc) 4515 { 4516 PetscErrorCode ierr; 4517 4518 PetscFunctionBegin; 4519 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); 4520 PetscValidHeaderSpecific(pc, SNES_CLASSID, 2); 4521 PetscCheckSameComm(snes, 1, pc, 2); 4522 ierr = PetscObjectReference((PetscObject) pc);CHKERRQ(ierr); 4523 ierr = SNESDestroy(&snes->pc);CHKERRQ(ierr); 4524 snes->pc = pc; 4525 ierr = PetscLogObjectParent(snes, snes->pc);CHKERRQ(ierr); 4526 PetscFunctionReturn(0); 4527 } 4528 4529 #undef __FUNCT__ 4530 #define __FUNCT__ "SNESGetPC" 4531 /*@ 4532 SNESGetPC - Returns a pointer to the nonlinear preconditioning context set with SNESSetPC(). 4533 4534 Not Collective 4535 4536 Input Parameter: 4537 . snes - iterative context obtained from SNESCreate() 4538 4539 Output Parameter: 4540 . pc - preconditioner context 4541 4542 Level: developer 4543 4544 .keywords: SNES, get, preconditioner 4545 .seealso: SNESSetPC() 4546 @*/ 4547 PetscErrorCode SNESGetPC(SNES snes, SNES *pc) 4548 { 4549 PetscErrorCode ierr; 4550 const char *optionsprefix; 4551 4552 PetscFunctionBegin; 4553 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); 4554 PetscValidPointer(pc, 2); 4555 if (!snes->pc) { 4556 ierr = SNESCreate(((PetscObject) snes)->comm,&snes->pc);CHKERRQ(ierr); 4557 ierr = PetscObjectIncrementTabLevel((PetscObject)snes->pc,(PetscObject)snes,1);CHKERRQ(ierr); 4558 ierr = PetscLogObjectParent(snes,snes->pc);CHKERRQ(ierr); 4559 ierr = SNESGetOptionsPrefix(snes,&optionsprefix);CHKERRQ(ierr); 4560 ierr = SNESSetOptionsPrefix(snes->pc,optionsprefix);CHKERRQ(ierr); 4561 ierr = SNESAppendOptionsPrefix(snes->pc,"npc_");CHKERRQ(ierr); 4562 } 4563 *pc = snes->pc; 4564 PetscFunctionReturn(0); 4565 } 4566 4567 #undef __FUNCT__ 4568 #define __FUNCT__ "SNESSetPCSide" 4569 /*@ 4570 SNESSetPCSide - Sets the preconditioning side. 4571 4572 Logically Collective on SNES 4573 4574 Input Parameter: 4575 . snes - iterative context obtained from SNESCreate() 4576 4577 Output Parameter: 4578 . side - the preconditioning side, where side is one of 4579 .vb 4580 PC_LEFT - left preconditioning (default) 4581 PC_RIGHT - right preconditioning 4582 .ve 4583 4584 Options Database Keys: 4585 . -snes_pc_side <right,left> 4586 4587 Level: intermediate 4588 4589 .keywords: SNES, set, right, left, side, preconditioner, flag 4590 4591 .seealso: SNESGetPCSide(), KSPSetPCSide() 4592 @*/ 4593 PetscErrorCode SNESSetPCSide(SNES snes,PCSide side) 4594 { 4595 PetscFunctionBegin; 4596 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4597 PetscValidLogicalCollectiveEnum(snes,side,2); 4598 snes->pcside = side; 4599 PetscFunctionReturn(0); 4600 } 4601 4602 #undef __FUNCT__ 4603 #define __FUNCT__ "SNESGetPCSide" 4604 /*@ 4605 SNESGetPCSide - Gets the preconditioning side. 4606 4607 Not Collective 4608 4609 Input Parameter: 4610 . snes - iterative context obtained from SNESCreate() 4611 4612 Output Parameter: 4613 . side - the preconditioning side, where side is one of 4614 .vb 4615 PC_LEFT - left preconditioning (default) 4616 PC_RIGHT - right preconditioning 4617 .ve 4618 4619 Level: intermediate 4620 4621 .keywords: SNES, get, right, left, side, preconditioner, flag 4622 4623 .seealso: SNESSetPCSide(), KSPGetPCSide() 4624 @*/ 4625 PetscErrorCode SNESGetPCSide(SNES snes,PCSide *side) 4626 { 4627 PetscFunctionBegin; 4628 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4629 PetscValidPointer(side,2); 4630 *side = snes->pcside; 4631 PetscFunctionReturn(0); 4632 } 4633 4634 #undef __FUNCT__ 4635 #define __FUNCT__ "SNESSetSNESLineSearch" 4636 /*@ 4637 SNESSetSNESLineSearch - Sets the linesearch on the SNES instance. 4638 4639 Collective on SNES 4640 4641 Input Parameters: 4642 + snes - iterative context obtained from SNESCreate() 4643 - linesearch - the linesearch object 4644 4645 Notes: 4646 Use SNESGetSNESLineSearch() to retrieve the preconditioner context (for example, 4647 to configure it using the API). 4648 4649 Level: developer 4650 4651 .keywords: SNES, set, linesearch 4652 .seealso: SNESGetSNESLineSearch() 4653 @*/ 4654 PetscErrorCode SNESSetSNESLineSearch(SNES snes, SNESLineSearch linesearch) 4655 { 4656 PetscErrorCode ierr; 4657 4658 PetscFunctionBegin; 4659 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); 4660 PetscValidHeaderSpecific(linesearch, SNESLINESEARCH_CLASSID, 2); 4661 PetscCheckSameComm(snes, 1, linesearch, 2); 4662 ierr = PetscObjectReference((PetscObject) linesearch);CHKERRQ(ierr); 4663 ierr = SNESLineSearchDestroy(&snes->linesearch);CHKERRQ(ierr); 4664 4665 snes->linesearch = linesearch; 4666 4667 ierr = PetscLogObjectParent(snes, snes->linesearch);CHKERRQ(ierr); 4668 PetscFunctionReturn(0); 4669 } 4670 4671 #undef __FUNCT__ 4672 #define __FUNCT__ "SNESGetSNESLineSearch" 4673 /*@C 4674 SNESGetSNESLineSearch - Returns a pointer to the line search context set with SNESSetLineSearch() 4675 or creates a default line search instance associated with the SNES and returns it. 4676 4677 Not Collective 4678 4679 Input Parameter: 4680 . snes - iterative context obtained from SNESCreate() 4681 4682 Output Parameter: 4683 . linesearch - linesearch context 4684 4685 Level: developer 4686 4687 .keywords: SNES, get, linesearch 4688 .seealso: SNESSetSNESLineSearch() 4689 @*/ 4690 PetscErrorCode SNESGetSNESLineSearch(SNES snes, SNESLineSearch *linesearch) 4691 { 4692 PetscErrorCode ierr; 4693 const char *optionsprefix; 4694 4695 PetscFunctionBegin; 4696 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); 4697 PetscValidPointer(linesearch, 2); 4698 if (!snes->linesearch) { 4699 ierr = SNESGetOptionsPrefix(snes, &optionsprefix);CHKERRQ(ierr); 4700 ierr = SNESLineSearchCreate(((PetscObject) snes)->comm, &snes->linesearch);CHKERRQ(ierr); 4701 ierr = SNESLineSearchSetSNES(snes->linesearch, snes);CHKERRQ(ierr); 4702 ierr = SNESLineSearchAppendOptionsPrefix(snes->linesearch, optionsprefix);CHKERRQ(ierr); 4703 ierr = PetscObjectIncrementTabLevel((PetscObject) snes->linesearch, (PetscObject) snes, 1);CHKERRQ(ierr); 4704 ierr = PetscLogObjectParent(snes, snes->linesearch);CHKERRQ(ierr); 4705 } 4706 *linesearch = snes->linesearch; 4707 PetscFunctionReturn(0); 4708 } 4709 4710 #if defined(PETSC_HAVE_MATLAB_ENGINE) 4711 #include <mex.h> 4712 4713 typedef struct {char *funcname; mxArray *ctx;} SNESMatlabContext; 4714 4715 #undef __FUNCT__ 4716 #define __FUNCT__ "SNESComputeFunction_Matlab" 4717 /* 4718 SNESComputeFunction_Matlab - Calls the function that has been set with SNESSetFunctionMatlab(). 4719 4720 Collective on SNES 4721 4722 Input Parameters: 4723 + snes - the SNES context 4724 - x - input vector 4725 4726 Output Parameter: 4727 . y - function vector, as set by SNESSetFunction() 4728 4729 Notes: 4730 SNESComputeFunction() is typically used within nonlinear solvers 4731 implementations, so most users would not generally call this routine 4732 themselves. 4733 4734 Level: developer 4735 4736 .keywords: SNES, nonlinear, compute, function 4737 4738 .seealso: SNESSetFunction(), SNESGetFunction() 4739 */ 4740 PetscErrorCode SNESComputeFunction_Matlab(SNES snes,Vec x,Vec y, void *ctx) 4741 { 4742 PetscErrorCode ierr; 4743 SNESMatlabContext *sctx = (SNESMatlabContext*)ctx; 4744 int nlhs = 1,nrhs = 5; 4745 mxArray *plhs[1],*prhs[5]; 4746 long long int lx = 0,ly = 0,ls = 0; 4747 4748 PetscFunctionBegin; 4749 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4750 PetscValidHeaderSpecific(x,VEC_CLASSID,2); 4751 PetscValidHeaderSpecific(y,VEC_CLASSID,3); 4752 PetscCheckSameComm(snes,1,x,2); 4753 PetscCheckSameComm(snes,1,y,3); 4754 4755 /* call Matlab function in ctx with arguments x and y */ 4756 4757 ierr = PetscMemcpy(&ls,&snes,sizeof(snes));CHKERRQ(ierr); 4758 ierr = PetscMemcpy(&lx,&x,sizeof(x));CHKERRQ(ierr); 4759 ierr = PetscMemcpy(&ly,&y,sizeof(x));CHKERRQ(ierr); 4760 prhs[0] = mxCreateDoubleScalar((double)ls); 4761 prhs[1] = mxCreateDoubleScalar((double)lx); 4762 prhs[2] = mxCreateDoubleScalar((double)ly); 4763 prhs[3] = mxCreateString(sctx->funcname); 4764 prhs[4] = sctx->ctx; 4765 ierr = mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscSNESComputeFunctionInternal");CHKERRQ(ierr); 4766 ierr = mxGetScalar(plhs[0]);CHKERRQ(ierr); 4767 mxDestroyArray(prhs[0]); 4768 mxDestroyArray(prhs[1]); 4769 mxDestroyArray(prhs[2]); 4770 mxDestroyArray(prhs[3]); 4771 mxDestroyArray(plhs[0]); 4772 PetscFunctionReturn(0); 4773 } 4774 4775 #undef __FUNCT__ 4776 #define __FUNCT__ "SNESSetFunctionMatlab" 4777 /* 4778 SNESSetFunctionMatlab - Sets the function evaluation routine and function 4779 vector for use by the SNES routines in solving systems of nonlinear 4780 equations from MATLAB. Here the function is a string containing the name of a MATLAB function 4781 4782 Logically Collective on SNES 4783 4784 Input Parameters: 4785 + snes - the SNES context 4786 . r - vector to store function value 4787 - SNESFunction - function evaluation routine 4788 4789 Notes: 4790 The Newton-like methods typically solve linear systems of the form 4791 $ f'(x) x = -f(x), 4792 where f'(x) denotes the Jacobian matrix and f(x) is the function. 4793 4794 Level: beginner 4795 4796 .keywords: SNES, nonlinear, set, function 4797 4798 .seealso: SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESSetFunction() 4799 */ 4800 PetscErrorCode SNESSetFunctionMatlab(SNES snes,Vec r,const char *SNESFunction,mxArray *ctx) 4801 { 4802 PetscErrorCode ierr; 4803 SNESMatlabContext *sctx; 4804 4805 PetscFunctionBegin; 4806 /* currently sctx is memory bleed */ 4807 ierr = PetscMalloc(sizeof(SNESMatlabContext),&sctx);CHKERRQ(ierr); 4808 ierr = PetscStrallocpy(SNESFunction,&sctx->funcname);CHKERRQ(ierr); 4809 /* 4810 This should work, but it doesn't 4811 sctx->ctx = ctx; 4812 mexMakeArrayPersistent(sctx->ctx); 4813 */ 4814 sctx->ctx = mxDuplicateArray(ctx); 4815 ierr = SNESSetFunction(snes,r,SNESComputeFunction_Matlab,sctx);CHKERRQ(ierr); 4816 PetscFunctionReturn(0); 4817 } 4818 4819 #undef __FUNCT__ 4820 #define __FUNCT__ "SNESComputeJacobian_Matlab" 4821 /* 4822 SNESComputeJacobian_Matlab - Calls the function that has been set with SNESSetJacobianMatlab(). 4823 4824 Collective on SNES 4825 4826 Input Parameters: 4827 + snes - the SNES context 4828 . x - input vector 4829 . A, B - the matrices 4830 - ctx - user context 4831 4832 Output Parameter: 4833 . flag - structure of the matrix 4834 4835 Level: developer 4836 4837 .keywords: SNES, nonlinear, compute, function 4838 4839 .seealso: SNESSetFunction(), SNESGetFunction() 4840 @*/ 4841 PetscErrorCode SNESComputeJacobian_Matlab(SNES snes,Vec x,Mat *A,Mat *B,MatStructure *flag, void *ctx) 4842 { 4843 PetscErrorCode ierr; 4844 SNESMatlabContext *sctx = (SNESMatlabContext*)ctx; 4845 int nlhs = 2,nrhs = 6; 4846 mxArray *plhs[2],*prhs[6]; 4847 long long int lx = 0,lA = 0,ls = 0, lB = 0; 4848 4849 PetscFunctionBegin; 4850 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4851 PetscValidHeaderSpecific(x,VEC_CLASSID,2); 4852 4853 /* call Matlab function in ctx with arguments x and y */ 4854 4855 ierr = PetscMemcpy(&ls,&snes,sizeof(snes));CHKERRQ(ierr); 4856 ierr = PetscMemcpy(&lx,&x,sizeof(x));CHKERRQ(ierr); 4857 ierr = PetscMemcpy(&lA,A,sizeof(x));CHKERRQ(ierr); 4858 ierr = PetscMemcpy(&lB,B,sizeof(x));CHKERRQ(ierr); 4859 prhs[0] = mxCreateDoubleScalar((double)ls); 4860 prhs[1] = mxCreateDoubleScalar((double)lx); 4861 prhs[2] = mxCreateDoubleScalar((double)lA); 4862 prhs[3] = mxCreateDoubleScalar((double)lB); 4863 prhs[4] = mxCreateString(sctx->funcname); 4864 prhs[5] = sctx->ctx; 4865 ierr = mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscSNESComputeJacobianInternal");CHKERRQ(ierr); 4866 ierr = mxGetScalar(plhs[0]);CHKERRQ(ierr); 4867 *flag = (MatStructure) mxGetScalar(plhs[1]);CHKERRQ(ierr); 4868 mxDestroyArray(prhs[0]); 4869 mxDestroyArray(prhs[1]); 4870 mxDestroyArray(prhs[2]); 4871 mxDestroyArray(prhs[3]); 4872 mxDestroyArray(prhs[4]); 4873 mxDestroyArray(plhs[0]); 4874 mxDestroyArray(plhs[1]); 4875 PetscFunctionReturn(0); 4876 } 4877 4878 #undef __FUNCT__ 4879 #define __FUNCT__ "SNESSetJacobianMatlab" 4880 /* 4881 SNESSetJacobianMatlab - Sets the Jacobian function evaluation routine and two empty Jacobian matrices 4882 vector for use by the SNES routines in solving systems of nonlinear 4883 equations from MATLAB. Here the function is a string containing the name of a MATLAB function 4884 4885 Logically Collective on SNES 4886 4887 Input Parameters: 4888 + snes - the SNES context 4889 . A,B - Jacobian matrices 4890 . SNESJacobianFunction - function evaluation routine 4891 - ctx - user context 4892 4893 Level: developer 4894 4895 .keywords: SNES, nonlinear, set, function 4896 4897 .seealso: SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESSetFunction(), SNESJacobianFunction 4898 */ 4899 PetscErrorCode SNESSetJacobianMatlab(SNES snes,Mat A,Mat B,const char *SNESJacobianFunction,mxArray *ctx) 4900 { 4901 PetscErrorCode ierr; 4902 SNESMatlabContext *sctx; 4903 4904 PetscFunctionBegin; 4905 /* currently sctx is memory bleed */ 4906 ierr = PetscMalloc(sizeof(SNESMatlabContext),&sctx);CHKERRQ(ierr); 4907 ierr = PetscStrallocpy(SNESJacobianFunction,&sctx->funcname);CHKERRQ(ierr); 4908 /* 4909 This should work, but it doesn't 4910 sctx->ctx = ctx; 4911 mexMakeArrayPersistent(sctx->ctx); 4912 */ 4913 sctx->ctx = mxDuplicateArray(ctx); 4914 ierr = SNESSetJacobian(snes,A,B,SNESComputeJacobian_Matlab,sctx);CHKERRQ(ierr); 4915 PetscFunctionReturn(0); 4916 } 4917 4918 #undef __FUNCT__ 4919 #define __FUNCT__ "SNESMonitor_Matlab" 4920 /* 4921 SNESMonitor_Matlab - Calls the function that has been set with SNESMonitorSetMatlab(). 4922 4923 Collective on SNES 4924 4925 .seealso: SNESSetFunction(), SNESGetFunction() 4926 @*/ 4927 PetscErrorCode SNESMonitor_Matlab(SNES snes,PetscInt it, PetscReal fnorm, void *ctx) 4928 { 4929 PetscErrorCode ierr; 4930 SNESMatlabContext *sctx = (SNESMatlabContext*)ctx; 4931 int nlhs = 1,nrhs = 6; 4932 mxArray *plhs[1],*prhs[6]; 4933 long long int lx = 0,ls = 0; 4934 Vec x = snes->vec_sol; 4935 4936 PetscFunctionBegin; 4937 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4938 4939 ierr = PetscMemcpy(&ls,&snes,sizeof(snes));CHKERRQ(ierr); 4940 ierr = PetscMemcpy(&lx,&x,sizeof(x));CHKERRQ(ierr); 4941 prhs[0] = mxCreateDoubleScalar((double)ls); 4942 prhs[1] = mxCreateDoubleScalar((double)it); 4943 prhs[2] = mxCreateDoubleScalar((double)fnorm); 4944 prhs[3] = mxCreateDoubleScalar((double)lx); 4945 prhs[4] = mxCreateString(sctx->funcname); 4946 prhs[5] = sctx->ctx; 4947 ierr = mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscSNESMonitorInternal");CHKERRQ(ierr); 4948 ierr = mxGetScalar(plhs[0]);CHKERRQ(ierr); 4949 mxDestroyArray(prhs[0]); 4950 mxDestroyArray(prhs[1]); 4951 mxDestroyArray(prhs[2]); 4952 mxDestroyArray(prhs[3]); 4953 mxDestroyArray(prhs[4]); 4954 mxDestroyArray(plhs[0]); 4955 PetscFunctionReturn(0); 4956 } 4957 4958 #undef __FUNCT__ 4959 #define __FUNCT__ "SNESMonitorSetMatlab" 4960 /* 4961 SNESMonitorSetMatlab - Sets the monitor function from MATLAB 4962 4963 Level: developer 4964 4965 .keywords: SNES, nonlinear, set, function 4966 4967 .seealso: SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESSetFunction() 4968 */ 4969 PetscErrorCode SNESMonitorSetMatlab(SNES snes,const char *SNESMonitorFunction,mxArray *ctx) 4970 { 4971 PetscErrorCode ierr; 4972 SNESMatlabContext *sctx; 4973 4974 PetscFunctionBegin; 4975 /* currently sctx is memory bleed */ 4976 ierr = PetscMalloc(sizeof(SNESMatlabContext),&sctx);CHKERRQ(ierr); 4977 ierr = PetscStrallocpy(SNESMonitorFunction,&sctx->funcname);CHKERRQ(ierr); 4978 /* 4979 This should work, but it doesn't 4980 sctx->ctx = ctx; 4981 mexMakeArrayPersistent(sctx->ctx); 4982 */ 4983 sctx->ctx = mxDuplicateArray(ctx); 4984 ierr = SNESMonitorSet(snes,SNESMonitor_Matlab,sctx,NULL);CHKERRQ(ierr); 4985 PetscFunctionReturn(0); 4986 } 4987 4988 #endif 4989