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