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 /*MC 1510 SNESFunction - function used to convey the nonlinear function to be solved by SNES 1511 1512 Synopsis: 1513 SNESFunction(SNES snes,Vec x,Vec f,void *ctx); 1514 1515 Input Parameters: 1516 + snes - the SNES context 1517 . x - state at which to evaluate residual 1518 - ctx - optional user-defined function context, passed in with SNESSetFunction() 1519 1520 Output Parameter: 1521 . f - vector to put residual (function value) 1522 1523 .seealso: SNESSetFunction(), SNESGetFunction() 1524 M*/ 1525 1526 #undef __FUNCT__ 1527 #define __FUNCT__ "SNESSetFunction" 1528 /*@C 1529 SNESSetFunction - Sets the function evaluation routine and function 1530 vector for use by the SNES routines in solving systems of nonlinear 1531 equations. 1532 1533 Logically Collective on SNES 1534 1535 Input Parameters: 1536 + snes - the SNES context 1537 . r - vector to store function value 1538 . SNESFunction - function evaluation routine 1539 - ctx - [optional] user-defined context for private data for the 1540 function evaluation routine (may be PETSC_NULL) 1541 1542 Notes: 1543 The Newton-like methods typically solve linear systems of the form 1544 $ f'(x) x = -f(x), 1545 where f'(x) denotes the Jacobian matrix and f(x) is the function. 1546 1547 Level: beginner 1548 1549 .keywords: SNES, nonlinear, set, function 1550 1551 .seealso: SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESSetPicard() 1552 @*/ 1553 PetscErrorCode SNESSetFunction(SNES snes,Vec r,PetscErrorCode (*SNESFunction)(SNES,Vec,Vec,void*),void *ctx) 1554 { 1555 PetscErrorCode ierr; 1556 DM dm; 1557 1558 PetscFunctionBegin; 1559 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1560 if (r) { 1561 PetscValidHeaderSpecific(r,VEC_CLASSID,2); 1562 PetscCheckSameComm(snes,1,r,2); 1563 ierr = PetscObjectReference((PetscObject)r);CHKERRQ(ierr); 1564 ierr = VecDestroy(&snes->vec_func);CHKERRQ(ierr); 1565 snes->vec_func = r; 1566 } 1567 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 1568 ierr = DMSNESSetFunction(dm,SNESFunction,ctx);CHKERRQ(ierr); 1569 PetscFunctionReturn(0); 1570 } 1571 1572 1573 #undef __FUNCT__ 1574 #define __FUNCT__ "SNESSetInitialFunction" 1575 /*@C 1576 SNESSetInitialFunction - Sets the function vector to be used as the 1577 function norm at the initialization of the method. In some 1578 instances, the user has precomputed the function before calling 1579 SNESSolve. This function allows one to avoid a redundant call 1580 to SNESComputeFunction in that case. 1581 1582 Logically Collective on SNES 1583 1584 Input Parameters: 1585 + snes - the SNES context 1586 - f - vector to store function value 1587 1588 Notes: 1589 This should not be modified during the solution procedure. 1590 1591 This is used extensively in the SNESFAS hierarchy and in nonlinear preconditioning. 1592 1593 Level: developer 1594 1595 .keywords: SNES, nonlinear, set, function 1596 1597 .seealso: SNESSetFunction(), SNESComputeFunction(), SNESSetInitialFunctionNorm() 1598 @*/ 1599 PetscErrorCode SNESSetInitialFunction(SNES snes, Vec f) 1600 { 1601 PetscErrorCode ierr; 1602 Vec vec_func; 1603 1604 PetscFunctionBegin; 1605 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1606 PetscValidHeaderSpecific(f,VEC_CLASSID,2); 1607 PetscCheckSameComm(snes,1,f,2); 1608 ierr = SNESGetFunction(snes,&vec_func,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr); 1609 ierr = VecCopy(f, vec_func);CHKERRQ(ierr); 1610 snes->vec_func_init_set = PETSC_TRUE; 1611 PetscFunctionReturn(0); 1612 } 1613 1614 1615 #undef __FUNCT__ 1616 #define __FUNCT__ "SNESSetInitialFunctionNorm" 1617 /*@C 1618 SNESSetInitialFunctionNorm - Sets the function norm to be used as the function norm 1619 at the initialization of the method. In some instances, the user has precomputed 1620 the function and its norm before calling SNESSolve. This function allows one to 1621 avoid a redundant call to SNESComputeFunction() and VecNorm() in that case. 1622 1623 Logically Collective on SNES 1624 1625 Input Parameters: 1626 + snes - the SNES context 1627 - fnorm - the norm of F as set by SNESSetInitialFunction() 1628 1629 This is used extensively in the SNESFAS hierarchy and in nonlinear preconditioning. 1630 1631 Level: developer 1632 1633 .keywords: SNES, nonlinear, set, function, norm 1634 1635 .seealso: SNESSetFunction(), SNESComputeFunction(), SNESSetInitialFunction() 1636 @*/ 1637 PetscErrorCode SNESSetInitialFunctionNorm(SNES snes, PetscReal fnorm) 1638 { 1639 PetscFunctionBegin; 1640 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1641 snes->norm_init = fnorm; 1642 snes->norm_init_set = PETSC_TRUE; 1643 PetscFunctionReturn(0); 1644 } 1645 1646 #undef __FUNCT__ 1647 #define __FUNCT__ "SNESSetNormType" 1648 /*@ 1649 SNESSetNormType - Sets the SNESNormType used in covergence and monitoring 1650 of the SNES method. 1651 1652 Logically Collective on SNES 1653 1654 Input Parameters: 1655 + snes - the SNES context 1656 - normtype - the type of the norm used 1657 1658 Notes: 1659 Only certain SNES methods support certain SNESNormTypes. Most require evaluation 1660 of the nonlinear function and the taking of its norm at every iteration to 1661 even ensure convergence at all. However, methods such as custom Gauss-Seidel methods 1662 (SNESGS) and the like do not require the norm of the function to be computed, and therfore 1663 may either be monitored for convergence or not. As these are often used as nonlinear 1664 preconditioners, monitoring the norm of their error is not a useful enterprise within 1665 their solution. 1666 1667 Level: developer 1668 1669 .keywords: SNES, nonlinear, set, function, norm, type 1670 1671 .seealso: SNESGetNormType(), SNESComputeFunction(), VecNorm(), SNESSetFunction(), SNESSetInitialFunction(), SNESNormType 1672 @*/ 1673 PetscErrorCode SNESSetNormType(SNES snes, SNESNormType normtype) 1674 { 1675 PetscFunctionBegin; 1676 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1677 snes->normtype = normtype; 1678 PetscFunctionReturn(0); 1679 } 1680 1681 1682 #undef __FUNCT__ 1683 #define __FUNCT__ "SNESGetNormType" 1684 /*@ 1685 SNESGetNormType - Gets the SNESNormType used in covergence and monitoring 1686 of the SNES method. 1687 1688 Logically Collective on SNES 1689 1690 Input Parameters: 1691 + snes - the SNES context 1692 - normtype - the type of the norm used 1693 1694 Level: advanced 1695 1696 .keywords: SNES, nonlinear, set, function, norm, type 1697 1698 .seealso: SNESSetNormType(), SNESComputeFunction(), VecNorm(), SNESSetFunction(), SNESSetInitialFunction(), SNESNormType 1699 @*/ 1700 PetscErrorCode SNESGetNormType(SNES snes, SNESNormType *normtype) 1701 { 1702 PetscFunctionBegin; 1703 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1704 *normtype = snes->normtype; 1705 PetscFunctionReturn(0); 1706 } 1707 1708 #undef __FUNCT__ 1709 #define __FUNCT__ "SNESSetGS" 1710 /*@C 1711 SNESSetGS - Sets the user nonlinear Gauss-Seidel routine for 1712 use with composed nonlinear solvers. 1713 1714 Input Parameters: 1715 + snes - the SNES context 1716 . gsfunc - function evaluation routine 1717 - ctx - [optional] user-defined context for private data for the 1718 smoother evaluation routine (may be PETSC_NULL) 1719 1720 Calling sequence of func: 1721 $ func (SNES snes,Vec x,Vec b,void *ctx); 1722 1723 + X - solution vector 1724 . B - RHS vector 1725 - ctx - optional user-defined Gauss-Seidel context 1726 1727 Notes: 1728 The GS routines are used by the composed nonlinear solver to generate 1729 a problem appropriate update to the solution, particularly FAS. 1730 1731 Level: intermediate 1732 1733 .keywords: SNES, nonlinear, set, Gauss-Seidel 1734 1735 .seealso: SNESGetFunction(), SNESComputeGS() 1736 @*/ 1737 PetscErrorCode SNESSetGS(SNES snes,PetscErrorCode (*gsfunc)(SNES,Vec,Vec,void*),void *ctx) 1738 { 1739 PetscErrorCode ierr; 1740 DM dm; 1741 1742 PetscFunctionBegin; 1743 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1744 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 1745 ierr = DMSNESSetGS(dm,gsfunc,ctx);CHKERRQ(ierr); 1746 PetscFunctionReturn(0); 1747 } 1748 1749 #undef __FUNCT__ 1750 #define __FUNCT__ "SNESSetGSSweeps" 1751 /*@ 1752 SNESSetGSSweeps - Sets the number of sweeps of GS to use. 1753 1754 Input Parameters: 1755 + snes - the SNES context 1756 - sweeps - the number of sweeps of GS to perform. 1757 1758 Level: intermediate 1759 1760 .keywords: SNES, nonlinear, set, Gauss-Siedel 1761 1762 .seealso: SNESSetGS(), SNESGetGS(), SNESSetPC(), SNESGetGSSweeps() 1763 @*/ 1764 1765 PetscErrorCode SNESSetGSSweeps(SNES snes, PetscInt sweeps) 1766 { 1767 PetscFunctionBegin; 1768 snes->gssweeps = sweeps; 1769 PetscFunctionReturn(0); 1770 } 1771 1772 1773 #undef __FUNCT__ 1774 #define __FUNCT__ "SNESGetGSSweeps" 1775 /*@ 1776 SNESGetGSSweeps - Gets the number of sweeps GS will use. 1777 1778 Input Parameters: 1779 . snes - the SNES context 1780 1781 Output Parameters: 1782 . sweeps - the number of sweeps of GS to perform. 1783 1784 Level: intermediate 1785 1786 .keywords: SNES, nonlinear, set, Gauss-Siedel 1787 1788 .seealso: SNESSetGS(), SNESGetGS(), SNESSetPC(), SNESSetGSSweeps() 1789 @*/ 1790 PetscErrorCode SNESGetGSSweeps(SNES snes, PetscInt * sweeps) 1791 { 1792 PetscFunctionBegin; 1793 *sweeps = snes->gssweeps; 1794 PetscFunctionReturn(0); 1795 } 1796 1797 #undef __FUNCT__ 1798 #define __FUNCT__ "SNESPicardComputeFunction" 1799 PetscErrorCode SNESPicardComputeFunction(SNES snes,Vec x,Vec f,void *ctx) 1800 { 1801 PetscErrorCode ierr; 1802 DM dm; 1803 DMSNES sdm; 1804 1805 PetscFunctionBegin; 1806 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 1807 ierr = DMGetDMSNES(dm,&sdm);CHKERRQ(ierr); 1808 /* A(x)*x - b(x) */ 1809 if (sdm->ops->computepfunction) { 1810 ierr = (*sdm->ops->computepfunction)(snes,x,f,sdm->pctx);CHKERRQ(ierr); 1811 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE, "Must call SNESSetPicard() to provide Picard function."); 1812 1813 if (sdm->ops->computepjacobian) { 1814 ierr = (*sdm->ops->computepjacobian)(snes,x,&snes->jacobian,&snes->jacobian_pre,&snes->matstruct,sdm->pctx);CHKERRQ(ierr); 1815 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE, "Must call SNESSetPicard() to provide Picard matrix."); 1816 ierr = VecScale(f,-1.0);CHKERRQ(ierr); 1817 ierr = MatMultAdd(snes->jacobian,x,f,f);CHKERRQ(ierr); 1818 PetscFunctionReturn(0); 1819 } 1820 1821 #undef __FUNCT__ 1822 #define __FUNCT__ "SNESPicardComputeJacobian" 1823 PetscErrorCode SNESPicardComputeJacobian(SNES snes,Vec x1,Mat *J,Mat *B,MatStructure *flag,void *ctx) 1824 { 1825 PetscFunctionBegin; 1826 /* the jacobian matrix should be pre-filled in SNESPicardComputeFunction */ 1827 *flag = snes->matstruct; 1828 PetscFunctionReturn(0); 1829 } 1830 1831 #undef __FUNCT__ 1832 #define __FUNCT__ "SNESSetPicard" 1833 /*@C 1834 SNESSetPicard - Use SNES to solve the semilinear-system A(x) x = b(x) via a Picard type iteration (Picard linearization) 1835 1836 Logically Collective on SNES 1837 1838 Input Parameters: 1839 + snes - the SNES context 1840 . r - vector to store function value 1841 . func - function evaluation routine 1842 . jmat - normally the same as mat but you can pass another matrix for which you compute A(x) x - b(x) (see jmat below) 1843 . mat - matrix to store A 1844 . mfunc - function to compute matrix value 1845 - ctx - [optional] user-defined context for private data for the 1846 function evaluation routine (may be PETSC_NULL) 1847 1848 Calling sequence of func: 1849 $ func (SNES snes,Vec x,Vec f,void *ctx); 1850 1851 + f - function vector 1852 - ctx - optional user-defined function context 1853 1854 Calling sequence of mfunc: 1855 $ mfunc (SNES snes,Vec x,Mat *jmat,Mat *mat,int *flag,void *ctx); 1856 1857 + x - input vector 1858 . 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(), 1859 normally just pass mat in this location 1860 . mat - form A(x) matrix 1861 . flag - flag indicating information about the preconditioner matrix 1862 structure (same as flag in KSPSetOperators()), one of SAME_NONZERO_PATTERN,DIFFERENT_NONZERO_PATTERN,SAME_PRECONDITIONER 1863 - ctx - [optional] user-defined Jacobian context 1864 1865 Notes: 1866 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 1867 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. 1868 1869 One can call SNESSetPicard() or SNESSetFunction() (and possibly SNESSetJacobian()) but cannot call both 1870 1871 $ 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} 1872 $ Note that when an exact solver is used this corresponds to the "classic" Picard A(x^{n}) x^{n+1} = b(x^{n}) iteration. 1873 1874 Run with -snes_mf_operator to solve the system with Newton's method using A(x^{n}) to construct the preconditioner. 1875 1876 We implement the defect correction form of the Picard iteration because it converges much more generally when inexact linear solvers are used then 1877 the direct Picard iteration A(x^n) x^{n+1} = b(x^n) 1878 1879 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 1880 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 1881 different please contact us at petsc-dev@mcs.anl.gov and we'll have an entirely new argument :-). 1882 1883 Level: intermediate 1884 1885 .keywords: SNES, nonlinear, set, function 1886 1887 .seealso: SNESGetFunction(), SNESSetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESGetPicard(), SNESLineSearchPreCheckPicard() 1888 @*/ 1889 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) 1890 { 1891 PetscErrorCode ierr; 1892 DM dm; 1893 1894 PetscFunctionBegin; 1895 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1896 ierr = SNESGetDM(snes, &dm);CHKERRQ(ierr); 1897 ierr = DMSNESSetPicard(dm,func,mfunc,ctx);CHKERRQ(ierr); 1898 ierr = SNESSetFunction(snes,r,SNESPicardComputeFunction,ctx);CHKERRQ(ierr); 1899 ierr = SNESSetJacobian(snes,jmat,mat,SNESPicardComputeJacobian,ctx);CHKERRQ(ierr); 1900 PetscFunctionReturn(0); 1901 } 1902 1903 1904 #undef __FUNCT__ 1905 #define __FUNCT__ "SNESGetPicard" 1906 /*@C 1907 SNESGetPicard - Returns the context for the Picard iteration 1908 1909 Not Collective, but Vec is parallel if SNES is parallel. Collective if Vec is requested, but has not been created yet. 1910 1911 Input Parameter: 1912 . snes - the SNES context 1913 1914 Output Parameter: 1915 + r - the function (or PETSC_NULL) 1916 . func - the function (or PETSC_NULL) 1917 . jmat - the picard matrix (or PETSC_NULL) 1918 . mat - the picard preconditioner matrix (or PETSC_NULL) 1919 . mfunc - the function for matrix evaluation (or PETSC_NULL) 1920 - ctx - the function context (or PETSC_NULL) 1921 1922 Level: advanced 1923 1924 .keywords: SNES, nonlinear, get, function 1925 1926 .seealso: SNESSetPicard, SNESGetFunction, SNESGetJacobian, SNESGetDM 1927 @*/ 1928 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) 1929 { 1930 PetscErrorCode ierr; 1931 DM dm; 1932 1933 PetscFunctionBegin; 1934 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1935 ierr = SNESGetFunction(snes,r,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr); 1936 ierr = SNESGetJacobian(snes,jmat,mat,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr); 1937 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 1938 ierr = DMSNESGetPicard(dm,func,mfunc,ctx);CHKERRQ(ierr); 1939 PetscFunctionReturn(0); 1940 } 1941 1942 #undef __FUNCT__ 1943 #define __FUNCT__ "SNESSetComputeInitialGuess" 1944 /*@C 1945 SNESSetComputeInitialGuess - Sets a routine used to compute an initial guess for the problem 1946 1947 Logically Collective on SNES 1948 1949 Input Parameters: 1950 + snes - the SNES context 1951 . func - function evaluation routine 1952 - ctx - [optional] user-defined context for private data for the 1953 function evaluation routine (may be PETSC_NULL) 1954 1955 Calling sequence of func: 1956 $ func (SNES snes,Vec x,void *ctx); 1957 1958 . f - function vector 1959 - ctx - optional user-defined function context 1960 1961 Level: intermediate 1962 1963 .keywords: SNES, nonlinear, set, function 1964 1965 .seealso: SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian() 1966 @*/ 1967 PetscErrorCode SNESSetComputeInitialGuess(SNES snes,PetscErrorCode (*func)(SNES,Vec,void*),void *ctx) 1968 { 1969 PetscFunctionBegin; 1970 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1971 if (func) snes->ops->computeinitialguess = func; 1972 if (ctx) snes->initialguessP = ctx; 1973 PetscFunctionReturn(0); 1974 } 1975 1976 /* --------------------------------------------------------------- */ 1977 #undef __FUNCT__ 1978 #define __FUNCT__ "SNESGetRhs" 1979 /*@C 1980 SNESGetRhs - Gets the vector for solving F(x) = rhs. If rhs is not set 1981 it assumes a zero right hand side. 1982 1983 Logically Collective on SNES 1984 1985 Input Parameter: 1986 . snes - the SNES context 1987 1988 Output Parameter: 1989 . rhs - the right hand side vector or PETSC_NULL if the right hand side vector is null 1990 1991 Level: intermediate 1992 1993 .keywords: SNES, nonlinear, get, function, right hand side 1994 1995 .seealso: SNESGetSolution(), SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESSetFunction() 1996 @*/ 1997 PetscErrorCode SNESGetRhs(SNES snes,Vec *rhs) 1998 { 1999 PetscFunctionBegin; 2000 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2001 PetscValidPointer(rhs,2); 2002 *rhs = snes->vec_rhs; 2003 PetscFunctionReturn(0); 2004 } 2005 2006 #undef __FUNCT__ 2007 #define __FUNCT__ "SNESComputeFunction" 2008 /*@ 2009 SNESComputeFunction - Calls the function that has been set with 2010 SNESSetFunction(). 2011 2012 Collective on SNES 2013 2014 Input Parameters: 2015 + snes - the SNES context 2016 - x - input vector 2017 2018 Output Parameter: 2019 . y - function vector, as set by SNESSetFunction() 2020 2021 Notes: 2022 SNESComputeFunction() is typically used within nonlinear solvers 2023 implementations, so most users would not generally call this routine 2024 themselves. 2025 2026 Level: developer 2027 2028 .keywords: SNES, nonlinear, compute, function 2029 2030 .seealso: SNESSetFunction(), SNESGetFunction() 2031 @*/ 2032 PetscErrorCode SNESComputeFunction(SNES snes,Vec x,Vec y) 2033 { 2034 PetscErrorCode ierr; 2035 DM dm; 2036 DMSNES sdm; 2037 2038 PetscFunctionBegin; 2039 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2040 PetscValidHeaderSpecific(x,VEC_CLASSID,2); 2041 PetscValidHeaderSpecific(y,VEC_CLASSID,3); 2042 PetscCheckSameComm(snes,1,x,2); 2043 PetscCheckSameComm(snes,1,y,3); 2044 VecValidValues(x,2,PETSC_TRUE); 2045 2046 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 2047 ierr = DMGetDMSNES(dm,&sdm);CHKERRQ(ierr); 2048 ierr = PetscLogEventBegin(SNES_FunctionEval,snes,x,y,0);CHKERRQ(ierr); 2049 if (snes->pc && snes->pcside == PC_LEFT) { 2050 ierr = VecCopy(x,y);CHKERRQ(ierr); 2051 ierr = SNESSolve(snes->pc,snes->vec_rhs,y);CHKERRQ(ierr); 2052 ierr = VecAYPX(y,-1.0,x);CHKERRQ(ierr); 2053 } else if (sdm->ops->computefunction) { 2054 PetscStackPush("SNES user function"); 2055 ierr = (*sdm->ops->computefunction)(snes,x,y,sdm->functionctx);CHKERRQ(ierr); 2056 PetscStackPop; 2057 } else if (snes->vec_rhs) { 2058 ierr = MatMult(snes->jacobian, x, y);CHKERRQ(ierr); 2059 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE, "Must call SNESSetFunction() or SNESSetDM() before SNESComputeFunction(), likely called from SNESSolve()."); 2060 if (snes->vec_rhs) { 2061 ierr = VecAXPY(y,-1.0,snes->vec_rhs);CHKERRQ(ierr); 2062 } 2063 snes->nfuncs++; 2064 ierr = PetscLogEventEnd(SNES_FunctionEval,snes,x,y,0);CHKERRQ(ierr); 2065 VecValidValues(y,3,PETSC_FALSE); 2066 PetscFunctionReturn(0); 2067 } 2068 2069 #undef __FUNCT__ 2070 #define __FUNCT__ "SNESComputeGS" 2071 /*@ 2072 SNESComputeGS - Calls the Gauss-Seidel function that has been set with 2073 SNESSetGS(). 2074 2075 Collective on SNES 2076 2077 Input Parameters: 2078 + snes - the SNES context 2079 . x - input vector 2080 - b - rhs vector 2081 2082 Output Parameter: 2083 . x - new solution vector 2084 2085 Notes: 2086 SNESComputeGS() is typically used within composed nonlinear solver 2087 implementations, so most users would not generally call this routine 2088 themselves. 2089 2090 Level: developer 2091 2092 .keywords: SNES, nonlinear, compute, function 2093 2094 .seealso: SNESSetGS(), SNESComputeFunction() 2095 @*/ 2096 PetscErrorCode SNESComputeGS(SNES snes,Vec b,Vec x) 2097 { 2098 PetscErrorCode ierr; 2099 PetscInt i; 2100 DM dm; 2101 DMSNES sdm; 2102 2103 PetscFunctionBegin; 2104 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2105 PetscValidHeaderSpecific(x,VEC_CLASSID,2); 2106 if (b) PetscValidHeaderSpecific(b,VEC_CLASSID,3); 2107 PetscCheckSameComm(snes,1,x,2); 2108 if (b) PetscCheckSameComm(snes,1,b,3); 2109 if (b) VecValidValues(b,2,PETSC_TRUE); 2110 ierr = PetscLogEventBegin(SNES_GSEval,snes,x,b,0);CHKERRQ(ierr); 2111 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 2112 ierr = DMGetDMSNES(dm,&sdm);CHKERRQ(ierr); 2113 if (sdm->ops->computegs) { 2114 for (i = 0; i < snes->gssweeps; i++) { 2115 PetscStackPush("SNES user GS"); 2116 ierr = (*sdm->ops->computegs)(snes,x,b,sdm->gsctx);CHKERRQ(ierr); 2117 PetscStackPop; 2118 } 2119 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE, "Must call SNESSetGS() before SNESComputeGS(), likely called from SNESSolve()."); 2120 ierr = PetscLogEventEnd(SNES_GSEval,snes,x,b,0);CHKERRQ(ierr); 2121 VecValidValues(x,3,PETSC_FALSE); 2122 PetscFunctionReturn(0); 2123 } 2124 2125 2126 #undef __FUNCT__ 2127 #define __FUNCT__ "SNESComputeJacobian" 2128 /*@ 2129 SNESComputeJacobian - Computes the Jacobian matrix that has been 2130 set with SNESSetJacobian(). 2131 2132 Collective on SNES and Mat 2133 2134 Input Parameters: 2135 + snes - the SNES context 2136 - x - input vector 2137 2138 Output Parameters: 2139 + A - Jacobian matrix 2140 . B - optional preconditioning matrix 2141 - flag - flag indicating matrix structure (one of, SAME_NONZERO_PATTERN,DIFFERENT_NONZERO_PATTERN,SAME_PRECONDITIONER) 2142 2143 Options Database Keys: 2144 + -snes_lag_preconditioner <lag> 2145 . -snes_lag_jacobian <lag> 2146 . -snes_compare_explicit - Compare the computed Jacobian to the finite difference Jacobian and output the differences 2147 . -snes_compare_explicit_draw - Compare the computed Jacobian to the finite difference Jacobian and draw the result 2148 . -snes_compare_explicit_contour - Compare the computed Jacobian to the finite difference Jacobian and draw a contour plot with the result 2149 . -snes_compare_operator - Make the comparison options above use the operator instead of the preconditioning matrix 2150 . -snes_compare_coloring - Compute the finite differece Jacobian using coloring and display norms of difference 2151 . -snes_compare_coloring_display - Compute the finite differece Jacobian using coloring and display verbose differences 2152 . -snes_compare_coloring_threshold - Display only those matrix entries that differ by more than a given threshold 2153 . -snes_compare_coloring_threshold_atol - Absolute tolerance for difference in matrix entries to be displayed by -snes_compare_coloring_threshold 2154 . -snes_compare_coloring_threshold_rtol - Relative tolerance for difference in matrix entries to be displayed by -snes_compare_coloring_threshold 2155 . -snes_compare_coloring_draw - Compute the finite differece Jacobian using coloring and draw differences 2156 - -snes_compare_coloring_draw_contour - Compute the finite differece Jacobian using coloring and show contours of matrices and differences 2157 2158 2159 Notes: 2160 Most users should not need to explicitly call this routine, as it 2161 is used internally within the nonlinear solvers. 2162 2163 See KSPSetOperators() for important information about setting the 2164 flag parameter. 2165 2166 Level: developer 2167 2168 .keywords: SNES, compute, Jacobian, matrix 2169 2170 .seealso: SNESSetJacobian(), KSPSetOperators(), MatStructure, SNESSetLagPreconditioner(), SNESSetLagJacobian() 2171 @*/ 2172 PetscErrorCode SNESComputeJacobian(SNES snes,Vec X,Mat *A,Mat *B,MatStructure *flg) 2173 { 2174 PetscErrorCode ierr; 2175 PetscBool flag; 2176 DM dm; 2177 DMSNES sdm; 2178 2179 PetscFunctionBegin; 2180 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2181 PetscValidHeaderSpecific(X,VEC_CLASSID,2); 2182 PetscValidPointer(flg,5); 2183 PetscCheckSameComm(snes,1,X,2); 2184 VecValidValues(X,2,PETSC_TRUE); 2185 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 2186 ierr = DMGetDMSNES(dm,&sdm);CHKERRQ(ierr); 2187 2188 if (!sdm->ops->computejacobian) SETERRQ(((PetscObject)snes)->comm,PETSC_ERR_USER,"Must call SNESSetJacobian(), DMSNESSetJacobian(), DMDASNESSetJacobianLocal(), etc"); 2189 2190 /* make sure that MatAssemblyBegin/End() is called on A matrix if it is matrix free */ 2191 2192 if (snes->lagjacobian == -2) { 2193 snes->lagjacobian = -1; 2194 ierr = PetscInfo(snes,"Recomputing Jacobian/preconditioner because lag is -2 (means compute Jacobian, but then never again) \n");CHKERRQ(ierr); 2195 } else if (snes->lagjacobian == -1) { 2196 *flg = SAME_PRECONDITIONER; 2197 ierr = PetscInfo(snes,"Reusing Jacobian/preconditioner because lag is -1\n");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 } else if (snes->lagjacobian > 1 && snes->iter % snes->lagjacobian) { 2205 *flg = SAME_PRECONDITIONER; 2206 ierr = PetscInfo2(snes,"Reusing Jacobian/preconditioner because lag is %D and SNES iteration is %D\n",snes->lagjacobian,snes->iter);CHKERRQ(ierr); 2207 ierr = PetscObjectTypeCompare((PetscObject)*A,MATMFFD,&flag);CHKERRQ(ierr); 2208 if (flag) { 2209 ierr = MatAssemblyBegin(*A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2210 ierr = MatAssemblyEnd(*A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2211 } 2212 PetscFunctionReturn(0); 2213 } 2214 2215 *flg = DIFFERENT_NONZERO_PATTERN; 2216 ierr = PetscLogEventBegin(SNES_JacobianEval,snes,X,*A,*B);CHKERRQ(ierr); 2217 2218 PetscStackPush("SNES user Jacobian function"); 2219 ierr = (*sdm->ops->computejacobian)(snes,X,A,B,flg,sdm->jacobianctx);CHKERRQ(ierr); 2220 PetscStackPop; 2221 2222 ierr = PetscLogEventEnd(SNES_JacobianEval,snes,X,*A,*B);CHKERRQ(ierr); 2223 2224 if (snes->lagpreconditioner == -2) { 2225 ierr = PetscInfo(snes,"Rebuilding preconditioner exactly once since lag is -2\n");CHKERRQ(ierr); 2226 snes->lagpreconditioner = -1; 2227 } else if (snes->lagpreconditioner == -1) { 2228 *flg = SAME_PRECONDITIONER; 2229 ierr = PetscInfo(snes,"Reusing preconditioner because lag is -1\n");CHKERRQ(ierr); 2230 } else if (snes->lagpreconditioner > 1 && snes->iter % snes->lagpreconditioner) { 2231 *flg = SAME_PRECONDITIONER; 2232 ierr = PetscInfo2(snes,"Reusing preconditioner because lag is %D and SNES iteration is %D\n",snes->lagpreconditioner,snes->iter);CHKERRQ(ierr); 2233 } 2234 2235 /* make sure user returned a correct Jacobian and preconditioner */ 2236 /* PetscValidHeaderSpecific(*A,MAT_CLASSID,3); 2237 PetscValidHeaderSpecific(*B,MAT_CLASSID,4); */ 2238 { 2239 PetscBool flag = PETSC_FALSE,flag_draw = PETSC_FALSE,flag_contour = PETSC_FALSE,flag_operator = PETSC_FALSE; 2240 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_explicit",&flag,PETSC_NULL);CHKERRQ(ierr); 2241 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_explicit_draw",&flag_draw,PETSC_NULL);CHKERRQ(ierr); 2242 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_explicit_draw_contour",&flag_contour,PETSC_NULL);CHKERRQ(ierr); 2243 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_operator",&flag_operator,PETSC_NULL);CHKERRQ(ierr); 2244 if (flag || flag_draw || flag_contour) { 2245 Mat Bexp_mine = PETSC_NULL,Bexp,FDexp; 2246 MatStructure mstruct; 2247 PetscViewer vdraw,vstdout; 2248 PetscBool flg; 2249 if (flag_operator) { 2250 ierr = MatComputeExplicitOperator(*A,&Bexp_mine);CHKERRQ(ierr); 2251 Bexp = Bexp_mine; 2252 } else { 2253 /* See if the preconditioning matrix can be viewed and added directly */ 2254 ierr = PetscObjectTypeCompareAny((PetscObject)*B,&flg,MATSEQAIJ,MATMPIAIJ,MATSEQDENSE,MATMPIDENSE,MATSEQBAIJ,MATMPIBAIJ,MATSEQSBAIJ,MATMPIBAIJ,"");CHKERRQ(ierr); 2255 if (flg) Bexp = *B; 2256 else { 2257 /* If the "preconditioning" matrix is itself MATSHELL or some other type without direct support */ 2258 ierr = MatComputeExplicitOperator(*B,&Bexp_mine);CHKERRQ(ierr); 2259 Bexp = Bexp_mine; 2260 } 2261 } 2262 ierr = MatConvert(Bexp,MATSAME,MAT_INITIAL_MATRIX,&FDexp);CHKERRQ(ierr); 2263 ierr = SNESDefaultComputeJacobian(snes,X,&FDexp,&FDexp,&mstruct,NULL);CHKERRQ(ierr); 2264 ierr = PetscViewerASCIIGetStdout(((PetscObject)snes)->comm,&vstdout);CHKERRQ(ierr); 2265 if (flag_draw || flag_contour) { 2266 ierr = PetscViewerDrawOpen(((PetscObject)snes)->comm,0,"Explicit Jacobians",PETSC_DECIDE,PETSC_DECIDE,300,300,&vdraw);CHKERRQ(ierr); 2267 if (flag_contour) {ierr = PetscViewerPushFormat(vdraw,PETSC_VIEWER_DRAW_CONTOUR);CHKERRQ(ierr);} 2268 } else vdraw = PETSC_NULL; 2269 ierr = PetscViewerASCIIPrintf(vstdout,"Explicit %s\n",flag_operator?"Jacobian":"preconditioning Jacobian");CHKERRQ(ierr); 2270 if (flag) {ierr = MatView(Bexp,vstdout);CHKERRQ(ierr);} 2271 if (vdraw) {ierr = MatView(Bexp,vdraw);CHKERRQ(ierr);} 2272 ierr = PetscViewerASCIIPrintf(vstdout,"Finite difference Jacobian\n");CHKERRQ(ierr); 2273 if (flag) {ierr = MatView(FDexp,vstdout);CHKERRQ(ierr);} 2274 if (vdraw) {ierr = MatView(FDexp,vdraw);CHKERRQ(ierr);} 2275 ierr = MatAYPX(FDexp,-1.0,Bexp,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 2276 ierr = PetscViewerASCIIPrintf(vstdout,"User-provided matrix minus finite difference Jacobian\n");CHKERRQ(ierr); 2277 if (flag) {ierr = MatView(FDexp,vstdout);CHKERRQ(ierr);} 2278 if (vdraw) { /* Always use contour for the difference */ 2279 ierr = PetscViewerPushFormat(vdraw,PETSC_VIEWER_DRAW_CONTOUR);CHKERRQ(ierr); 2280 ierr = MatView(FDexp,vdraw);CHKERRQ(ierr); 2281 ierr = PetscViewerPopFormat(vdraw);CHKERRQ(ierr); 2282 } 2283 if (flag_contour) {ierr = PetscViewerPopFormat(vdraw);CHKERRQ(ierr);} 2284 ierr = PetscViewerDestroy(&vdraw);CHKERRQ(ierr); 2285 ierr = MatDestroy(&Bexp_mine);CHKERRQ(ierr); 2286 ierr = MatDestroy(&FDexp);CHKERRQ(ierr); 2287 } 2288 } 2289 { 2290 PetscBool flag = PETSC_FALSE,flag_display = PETSC_FALSE,flag_draw = PETSC_FALSE,flag_contour = PETSC_FALSE,flag_threshold = PETSC_FALSE; 2291 PetscReal threshold_atol = PETSC_SQRT_MACHINE_EPSILON,threshold_rtol = 10*PETSC_SQRT_MACHINE_EPSILON; 2292 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_coloring",&flag,PETSC_NULL);CHKERRQ(ierr); 2293 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_coloring_display",&flag_display,PETSC_NULL);CHKERRQ(ierr); 2294 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_coloring_draw",&flag_draw,PETSC_NULL);CHKERRQ(ierr); 2295 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_coloring_draw_contour",&flag_contour,PETSC_NULL);CHKERRQ(ierr); 2296 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_coloring_threshold",&flag_threshold,PETSC_NULL);CHKERRQ(ierr); 2297 ierr = PetscOptionsGetReal(((PetscObject)snes)->prefix,"-snes_compare_coloring_threshold_rtol",&threshold_rtol,PETSC_NULL);CHKERRQ(ierr); 2298 ierr = PetscOptionsGetReal(((PetscObject)snes)->prefix,"-snes_compare_coloring_threshold_atol",&threshold_atol,PETSC_NULL);CHKERRQ(ierr); 2299 if (flag || flag_display || flag_draw || flag_contour || flag_threshold) { 2300 Mat Bfd; 2301 MatStructure mstruct; 2302 PetscViewer vdraw,vstdout; 2303 ISColoring iscoloring; 2304 MatFDColoring matfdcoloring; 2305 PetscErrorCode (*func)(SNES,Vec,Vec,void*); 2306 void *funcctx; 2307 PetscReal norm1,norm2,normmax; 2308 2309 ierr = MatDuplicate(*B,MAT_DO_NOT_COPY_VALUES,&Bfd);CHKERRQ(ierr); 2310 ierr = MatGetColoring(Bfd,MATCOLORINGSL,&iscoloring);CHKERRQ(ierr); 2311 ierr = MatFDColoringCreate(Bfd,iscoloring,&matfdcoloring);CHKERRQ(ierr); 2312 ierr = ISColoringDestroy(&iscoloring);CHKERRQ(ierr); 2313 2314 /* This method of getting the function is currently unreliable since it doesn't work for DM local functions. */ 2315 ierr = SNESGetFunction(snes,PETSC_NULL,&func,&funcctx);CHKERRQ(ierr); 2316 ierr = MatFDColoringSetFunction(matfdcoloring,(PetscErrorCode(*)(void))func,funcctx);CHKERRQ(ierr); 2317 ierr = PetscObjectSetOptionsPrefix((PetscObject)matfdcoloring,((PetscObject)snes)->prefix);CHKERRQ(ierr); 2318 ierr = PetscObjectAppendOptionsPrefix((PetscObject)matfdcoloring,"coloring_");CHKERRQ(ierr); 2319 ierr = MatFDColoringSetFromOptions(matfdcoloring);CHKERRQ(ierr); 2320 ierr = MatFDColoringApply(Bfd,matfdcoloring,X,&mstruct,snes);CHKERRQ(ierr); 2321 ierr = MatFDColoringDestroy(&matfdcoloring);CHKERRQ(ierr); 2322 2323 ierr = PetscViewerASCIIGetStdout(((PetscObject)snes)->comm,&vstdout);CHKERRQ(ierr); 2324 if (flag_draw || flag_contour) { 2325 ierr = PetscViewerDrawOpen(((PetscObject)snes)->comm,0,"Colored Jacobians",PETSC_DECIDE,PETSC_DECIDE,300,300,&vdraw);CHKERRQ(ierr); 2326 if (flag_contour) {ierr = PetscViewerPushFormat(vdraw,PETSC_VIEWER_DRAW_CONTOUR);CHKERRQ(ierr);} 2327 } else vdraw = PETSC_NULL; 2328 ierr = PetscViewerASCIIPrintf(vstdout,"Explicit preconditioning Jacobian\n");CHKERRQ(ierr); 2329 if (flag_display) {ierr = MatView(*B,vstdout);CHKERRQ(ierr);} 2330 if (vdraw) {ierr = MatView(*B,vdraw);CHKERRQ(ierr);} 2331 ierr = PetscViewerASCIIPrintf(vstdout,"Colored Finite difference Jacobian\n");CHKERRQ(ierr); 2332 if (flag_display) {ierr = MatView(Bfd,vstdout);CHKERRQ(ierr);} 2333 if (vdraw) {ierr = MatView(Bfd,vdraw);CHKERRQ(ierr);} 2334 ierr = MatAYPX(Bfd,-1.0,*B,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 2335 ierr = MatNorm(Bfd,NORM_1,&norm1);CHKERRQ(ierr); 2336 ierr = MatNorm(Bfd,NORM_FROBENIUS,&norm2);CHKERRQ(ierr); 2337 ierr = MatNorm(Bfd,NORM_MAX,&normmax);CHKERRQ(ierr); 2338 ierr = PetscViewerASCIIPrintf(vstdout,"User-provided matrix minus finite difference Jacobian, norm1=%G normFrob=%G normmax=%G\n",norm1,norm2,normmax);CHKERRQ(ierr); 2339 if (flag_display) {ierr = MatView(Bfd,vstdout);CHKERRQ(ierr);} 2340 if (vdraw) { /* Always use contour for the difference */ 2341 ierr = PetscViewerPushFormat(vdraw,PETSC_VIEWER_DRAW_CONTOUR);CHKERRQ(ierr); 2342 ierr = MatView(Bfd,vdraw);CHKERRQ(ierr); 2343 ierr = PetscViewerPopFormat(vdraw);CHKERRQ(ierr); 2344 } 2345 if (flag_contour) {ierr = PetscViewerPopFormat(vdraw);CHKERRQ(ierr);} 2346 2347 if (flag_threshold) { 2348 PetscInt bs,rstart,rend,i; 2349 ierr = MatGetBlockSize(*B,&bs);CHKERRQ(ierr); 2350 ierr = MatGetOwnershipRange(*B,&rstart,&rend);CHKERRQ(ierr); 2351 for (i=rstart; i<rend; i++) { 2352 const PetscScalar *ba,*ca; 2353 const PetscInt *bj,*cj; 2354 PetscInt bn,cn,j,maxentrycol = -1,maxdiffcol = -1,maxrdiffcol = -1; 2355 PetscReal maxentry = 0,maxdiff = 0,maxrdiff = 0; 2356 ierr = MatGetRow(*B,i,&bn,&bj,&ba);CHKERRQ(ierr); 2357 ierr = MatGetRow(Bfd,i,&cn,&cj,&ca);CHKERRQ(ierr); 2358 if (bn != cn) SETERRQ(((PetscObject)*A)->comm,PETSC_ERR_PLIB,"Unexpected different nonzero pattern in -snes_compare_coloring_threshold"); 2359 for (j=0; j<bn; j++) { 2360 PetscReal rdiff = PetscAbsScalar(ca[j]) / (threshold_atol + threshold_rtol*PetscAbsScalar(ba[j])); 2361 if (PetscAbsScalar(ba[j]) > PetscAbs(maxentry)) { 2362 maxentrycol = bj[j]; 2363 maxentry = PetscRealPart(ba[j]); 2364 } 2365 if (PetscAbsScalar(ca[j]) > PetscAbs(maxdiff)) { 2366 maxdiffcol = bj[j]; 2367 maxdiff = PetscRealPart(ca[j]); 2368 } 2369 if (rdiff > maxrdiff) { 2370 maxrdiffcol = bj[j]; 2371 maxrdiff = rdiff; 2372 } 2373 } 2374 if (maxrdiff > 1) { 2375 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); 2376 for (j=0; j<bn; j++) { 2377 PetscReal rdiff; 2378 rdiff = PetscAbsScalar(ca[j]) / (threshold_atol + threshold_rtol*PetscAbsScalar(ba[j])); 2379 if (rdiff > 1) { 2380 ierr = PetscViewerASCIIPrintf(vstdout," (%D,%G:%G)",bj[j],PetscRealPart(ba[j]),PetscRealPart(ca[j]));CHKERRQ(ierr); 2381 } 2382 } 2383 ierr = PetscViewerASCIIPrintf(vstdout,"\n",i,maxentry,maxdiff,maxrdiff);CHKERRQ(ierr); 2384 } 2385 ierr = MatRestoreRow(*B,i,&bn,&bj,&ba);CHKERRQ(ierr); 2386 ierr = MatRestoreRow(Bfd,i,&cn,&cj,&ca);CHKERRQ(ierr); 2387 } 2388 } 2389 ierr = PetscViewerDestroy(&vdraw);CHKERRQ(ierr); 2390 ierr = MatDestroy(&Bfd);CHKERRQ(ierr); 2391 } 2392 } 2393 PetscFunctionReturn(0); 2394 } 2395 2396 #undef __FUNCT__ 2397 #define __FUNCT__ "SNESSetJacobian" 2398 /*@C 2399 SNESSetJacobian - Sets the function to compute Jacobian as well as the 2400 location to store the matrix. 2401 2402 Logically Collective on SNES and Mat 2403 2404 Input Parameters: 2405 + snes - the SNES context 2406 . A - Jacobian matrix 2407 . B - preconditioner matrix (usually same as the Jacobian) 2408 . func - Jacobian evaluation routine (if PETSC_NULL then SNES retains any previously set value) 2409 - ctx - [optional] user-defined context for private data for the 2410 Jacobian evaluation routine (may be PETSC_NULL) (if PETSC_NULL then SNES retains any previously set value) 2411 2412 Calling sequence of func: 2413 $ func (SNES snes,Vec x,Mat *A,Mat *B,int *flag,void *ctx); 2414 2415 + x - input vector 2416 . A - Jacobian matrix 2417 . B - preconditioner matrix, usually the same as A 2418 . flag - flag indicating information about the preconditioner matrix 2419 structure (same as flag in KSPSetOperators()), one of SAME_NONZERO_PATTERN,DIFFERENT_NONZERO_PATTERN,SAME_PRECONDITIONER 2420 - ctx - [optional] user-defined Jacobian context 2421 2422 Notes: 2423 See KSPSetOperators() for important information about setting the flag 2424 output parameter in the routine func(). Be sure to read this information! 2425 2426 The routine func() takes Mat * as the matrix arguments rather than Mat. 2427 This allows the Jacobian evaluation routine to replace A and/or B with a 2428 completely new new matrix structure (not just different matrix elements) 2429 when appropriate, for instance, if the nonzero structure is changing 2430 throughout the global iterations. 2431 2432 If the A matrix and B matrix are different you must call MatAssemblyBegin/End() on 2433 each matrix. 2434 2435 If using SNESDefaultComputeJacobianColor() to assemble a Jacobian, the ctx argument 2436 must be a MatFDColoring. 2437 2438 Other defect-correction schemes can be used by computing a different matrix in place of the Jacobian. One common 2439 example is to use the "Picard linearization" which only differentiates through the highest order parts of each term. 2440 2441 Level: beginner 2442 2443 .keywords: SNES, nonlinear, set, Jacobian, matrix 2444 2445 .seealso: KSPSetOperators(), SNESSetFunction(), MatMFFDComputeJacobian(), SNESDefaultComputeJacobianColor(), MatStructure 2446 @*/ 2447 PetscErrorCode SNESSetJacobian(SNES snes,Mat A,Mat B,PetscErrorCode (*func)(SNES,Vec,Mat*,Mat*,MatStructure*,void*),void *ctx) 2448 { 2449 PetscErrorCode ierr; 2450 DM dm; 2451 2452 PetscFunctionBegin; 2453 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2454 if (A) PetscValidHeaderSpecific(A,MAT_CLASSID,2); 2455 if (B) PetscValidHeaderSpecific(B,MAT_CLASSID,3); 2456 if (A) PetscCheckSameComm(snes,1,A,2); 2457 if (B) PetscCheckSameComm(snes,1,B,3); 2458 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 2459 ierr = DMSNESSetJacobian(dm,func,ctx);CHKERRQ(ierr); 2460 if (A) { 2461 ierr = PetscObjectReference((PetscObject)A);CHKERRQ(ierr); 2462 ierr = MatDestroy(&snes->jacobian);CHKERRQ(ierr); 2463 snes->jacobian = A; 2464 } 2465 if (B) { 2466 ierr = PetscObjectReference((PetscObject)B);CHKERRQ(ierr); 2467 ierr = MatDestroy(&snes->jacobian_pre);CHKERRQ(ierr); 2468 snes->jacobian_pre = B; 2469 } 2470 PetscFunctionReturn(0); 2471 } 2472 2473 #undef __FUNCT__ 2474 #define __FUNCT__ "SNESGetJacobian" 2475 /*@C 2476 SNESGetJacobian - Returns the Jacobian matrix and optionally the user 2477 provided context for evaluating the Jacobian. 2478 2479 Not Collective, but Mat object will be parallel if SNES object is 2480 2481 Input Parameter: 2482 . snes - the nonlinear solver context 2483 2484 Output Parameters: 2485 + A - location to stash Jacobian matrix (or PETSC_NULL) 2486 . B - location to stash preconditioner matrix (or PETSC_NULL) 2487 . func - location to put Jacobian function (or PETSC_NULL) 2488 - ctx - location to stash Jacobian ctx (or PETSC_NULL) 2489 2490 Level: advanced 2491 2492 .seealso: SNESSetJacobian(), SNESComputeJacobian() 2493 @*/ 2494 PetscErrorCode SNESGetJacobian(SNES snes,Mat *A,Mat *B,PetscErrorCode (**func)(SNES,Vec,Mat*,Mat*,MatStructure*,void*),void **ctx) 2495 { 2496 PetscErrorCode ierr; 2497 DM dm; 2498 DMSNES sdm; 2499 2500 PetscFunctionBegin; 2501 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2502 if (A) *A = snes->jacobian; 2503 if (B) *B = snes->jacobian_pre; 2504 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 2505 ierr = DMGetDMSNES(dm,&sdm);CHKERRQ(ierr); 2506 if (func) *func = sdm->ops->computejacobian; 2507 if (ctx) *ctx = sdm->jacobianctx; 2508 PetscFunctionReturn(0); 2509 } 2510 2511 /* ----- Routines to initialize and destroy a nonlinear solver ---- */ 2512 2513 #undef __FUNCT__ 2514 #define __FUNCT__ "SNESSetUp" 2515 /*@ 2516 SNESSetUp - Sets up the internal data structures for the later use 2517 of a nonlinear solver. 2518 2519 Collective on SNES 2520 2521 Input Parameters: 2522 . snes - the SNES context 2523 2524 Notes: 2525 For basic use of the SNES solvers the user need not explicitly call 2526 SNESSetUp(), since these actions will automatically occur during 2527 the call to SNESSolve(). However, if one wishes to control this 2528 phase separately, SNESSetUp() should be called after SNESCreate() 2529 and optional routines of the form SNESSetXXX(), but before SNESSolve(). 2530 2531 Level: advanced 2532 2533 .keywords: SNES, nonlinear, setup 2534 2535 .seealso: SNESCreate(), SNESSolve(), SNESDestroy() 2536 @*/ 2537 PetscErrorCode SNESSetUp(SNES snes) 2538 { 2539 PetscErrorCode ierr; 2540 DM dm; 2541 DMSNES sdm; 2542 SNESLineSearch linesearch, pclinesearch; 2543 void *lsprectx,*lspostctx; 2544 SNESLineSearchPreCheckFunc lsprefunc; 2545 SNESLineSearchPostCheckFunc lspostfunc; 2546 PetscErrorCode (*func)(SNES,Vec,Vec,void*); 2547 Vec f,fpc; 2548 void *funcctx; 2549 PetscErrorCode (*jac)(SNES,Vec,Mat*,Mat*,MatStructure*,void*); 2550 void *jacctx,*appctx; 2551 Mat A,B; 2552 2553 PetscFunctionBegin; 2554 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2555 if (snes->setupcalled) PetscFunctionReturn(0); 2556 2557 if (!((PetscObject)snes)->type_name) { 2558 ierr = SNESSetType(snes,SNESNEWTONLS);CHKERRQ(ierr); 2559 } 2560 2561 ierr = SNESGetFunction(snes,&snes->vec_func,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr); 2562 if (snes->vec_func == snes->vec_sol) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_IDN,"Solution vector cannot be function vector"); 2563 if (snes->vec_rhs == snes->vec_sol) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_IDN,"Solution vector cannot be right hand side vector"); 2564 2565 if (!snes->vec_sol_update /* && snes->vec_sol */) { 2566 ierr = VecDuplicate(snes->vec_sol,&snes->vec_sol_update);CHKERRQ(ierr); 2567 ierr = PetscLogObjectParent(snes,snes->vec_sol_update);CHKERRQ(ierr); 2568 } 2569 2570 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 2571 ierr = DMShellSetGlobalVector(snes->dm,snes->vec_sol);CHKERRQ(ierr); 2572 ierr = DMGetDMSNES(dm,&sdm);CHKERRQ(ierr); 2573 if (!sdm->ops->computefunction) SETERRQ(((PetscObject)dm)->comm,PETSC_ERR_ARG_WRONGSTATE,"Function never provided to SNES object"); 2574 if (!sdm->ops->computejacobian) { 2575 ierr = DMSNESSetJacobian(dm,SNESDefaultComputeJacobianColor,PETSC_NULL);CHKERRQ(ierr); 2576 } 2577 if (!snes->vec_func) { 2578 ierr = DMCreateGlobalVector(dm,&snes->vec_func);CHKERRQ(ierr); 2579 } 2580 2581 if (!snes->ksp) { 2582 ierr = SNESGetKSP(snes, &snes->ksp);CHKERRQ(ierr); 2583 } 2584 2585 if (!snes->linesearch) { 2586 ierr = SNESGetSNESLineSearch(snes, &snes->linesearch);CHKERRQ(ierr); 2587 } 2588 2589 if (snes->pc && (snes->pcside == PC_LEFT)) { 2590 snes->mf = PETSC_TRUE; 2591 snes->mf_operator = PETSC_FALSE; 2592 } 2593 2594 if (snes->mf) { 2595 ierr = SNESSetUpMatrixFree_Private(snes, snes->mf_operator, snes->mf_version);CHKERRQ(ierr); 2596 } 2597 2598 if (snes->ops->usercompute && !snes->user) { 2599 ierr = (*snes->ops->usercompute)(snes,(void**)&snes->user);CHKERRQ(ierr); 2600 } 2601 2602 if (snes->pc) { 2603 /* copy the DM over */ 2604 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 2605 ierr = SNESSetDM(snes->pc,dm);CHKERRQ(ierr); 2606 2607 ierr = SNESGetFunction(snes,&f,&func,&funcctx);CHKERRQ(ierr); 2608 ierr = VecDuplicate(f,&fpc);CHKERRQ(ierr); 2609 ierr = SNESSetFunction(snes->pc,fpc,func,funcctx);CHKERRQ(ierr); 2610 ierr = SNESGetJacobian(snes,&A,&B,&jac,&jacctx);CHKERRQ(ierr); 2611 ierr = SNESSetJacobian(snes->pc,A,B,jac,jacctx);CHKERRQ(ierr); 2612 ierr = SNESGetApplicationContext(snes,&appctx);CHKERRQ(ierr); 2613 ierr = SNESSetApplicationContext(snes->pc,appctx);CHKERRQ(ierr); 2614 ierr = VecDestroy(&fpc);CHKERRQ(ierr); 2615 2616 /* copy the function pointers over */ 2617 ierr = PetscObjectCopyFortranFunctionPointers((PetscObject)snes,(PetscObject)snes->pc);CHKERRQ(ierr); 2618 2619 /* default to 1 iteration */ 2620 ierr = SNESSetTolerances(snes->pc, 0.0, 0.0, 0.0, 1, snes->pc->max_funcs);CHKERRQ(ierr); 2621 ierr = SNESSetNormType(snes->pc, SNES_NORM_FINAL_ONLY);CHKERRQ(ierr); 2622 ierr = SNESSetFromOptions(snes->pc);CHKERRQ(ierr); 2623 2624 /* copy the line search context over */ 2625 ierr = SNESGetSNESLineSearch(snes,&linesearch);CHKERRQ(ierr); 2626 ierr = SNESGetSNESLineSearch(snes->pc,&pclinesearch);CHKERRQ(ierr); 2627 ierr = SNESLineSearchGetPreCheck(linesearch,&lsprefunc,&lsprectx);CHKERRQ(ierr); 2628 ierr = SNESLineSearchGetPostCheck(linesearch,&lspostfunc,&lspostctx);CHKERRQ(ierr); 2629 ierr = SNESLineSearchSetPreCheck(pclinesearch,lsprefunc,lsprectx);CHKERRQ(ierr); 2630 ierr = SNESLineSearchSetPostCheck(pclinesearch,lspostfunc,lspostctx);CHKERRQ(ierr); 2631 ierr = PetscObjectCopyFortranFunctionPointers((PetscObject)linesearch, (PetscObject)pclinesearch);CHKERRQ(ierr); 2632 } 2633 2634 if (snes->ops->setup) { 2635 ierr = (*snes->ops->setup)(snes);CHKERRQ(ierr); 2636 } 2637 2638 snes->setupcalled = PETSC_TRUE; 2639 PetscFunctionReturn(0); 2640 } 2641 2642 #undef __FUNCT__ 2643 #define __FUNCT__ "SNESReset" 2644 /*@ 2645 SNESReset - Resets a SNES context to the snessetupcalled = 0 state and removes any allocated Vecs and Mats 2646 2647 Collective on SNES 2648 2649 Input Parameter: 2650 . snes - iterative context obtained from SNESCreate() 2651 2652 Level: intermediate 2653 2654 Notes: Also calls the application context destroy routine set with SNESSetComputeApplicationContext() 2655 2656 .keywords: SNES, destroy 2657 2658 .seealso: SNESCreate(), SNESSetUp(), SNESSolve() 2659 @*/ 2660 PetscErrorCode SNESReset(SNES snes) 2661 { 2662 PetscErrorCode ierr; 2663 2664 PetscFunctionBegin; 2665 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2666 if (snes->ops->userdestroy && snes->user) { 2667 ierr = (*snes->ops->userdestroy)((void**)&snes->user);CHKERRQ(ierr); 2668 snes->user = PETSC_NULL; 2669 } 2670 if (snes->pc) { 2671 ierr = SNESReset(snes->pc);CHKERRQ(ierr); 2672 } 2673 2674 if (snes->ops->reset) { 2675 ierr = (*snes->ops->reset)(snes);CHKERRQ(ierr); 2676 } 2677 if (snes->ksp) { 2678 ierr = KSPReset(snes->ksp);CHKERRQ(ierr); 2679 } 2680 2681 if (snes->linesearch) { 2682 ierr = SNESLineSearchReset(snes->linesearch);CHKERRQ(ierr); 2683 } 2684 2685 ierr = VecDestroy(&snes->vec_rhs);CHKERRQ(ierr); 2686 ierr = VecDestroy(&snes->vec_sol);CHKERRQ(ierr); 2687 ierr = VecDestroy(&snes->vec_sol_update);CHKERRQ(ierr); 2688 ierr = VecDestroy(&snes->vec_func);CHKERRQ(ierr); 2689 ierr = MatDestroy(&snes->jacobian);CHKERRQ(ierr); 2690 ierr = MatDestroy(&snes->jacobian_pre);CHKERRQ(ierr); 2691 ierr = VecDestroyVecs(snes->nwork,&snes->work);CHKERRQ(ierr); 2692 ierr = VecDestroyVecs(snes->nvwork,&snes->vwork);CHKERRQ(ierr); 2693 snes->nwork = snes->nvwork = 0; 2694 snes->setupcalled = PETSC_FALSE; 2695 PetscFunctionReturn(0); 2696 } 2697 2698 #undef __FUNCT__ 2699 #define __FUNCT__ "SNESDestroy" 2700 /*@ 2701 SNESDestroy - Destroys the nonlinear solver context that was created 2702 with SNESCreate(). 2703 2704 Collective on SNES 2705 2706 Input Parameter: 2707 . snes - the SNES context 2708 2709 Level: beginner 2710 2711 .keywords: SNES, nonlinear, destroy 2712 2713 .seealso: SNESCreate(), SNESSolve() 2714 @*/ 2715 PetscErrorCode SNESDestroy(SNES *snes) 2716 { 2717 PetscErrorCode ierr; 2718 2719 PetscFunctionBegin; 2720 if (!*snes) PetscFunctionReturn(0); 2721 PetscValidHeaderSpecific((*snes),SNES_CLASSID,1); 2722 if (--((PetscObject)(*snes))->refct > 0) {*snes = 0; PetscFunctionReturn(0);} 2723 2724 ierr = SNESReset((*snes));CHKERRQ(ierr); 2725 ierr = SNESDestroy(&(*snes)->pc);CHKERRQ(ierr); 2726 2727 /* if memory was published with AMS then destroy it */ 2728 ierr = PetscObjectDepublish((*snes));CHKERRQ(ierr); 2729 if ((*snes)->ops->destroy) {ierr = (*((*snes))->ops->destroy)((*snes));CHKERRQ(ierr);} 2730 2731 ierr = DMDestroy(&(*snes)->dm);CHKERRQ(ierr); 2732 ierr = KSPDestroy(&(*snes)->ksp);CHKERRQ(ierr); 2733 ierr = SNESLineSearchDestroy(&(*snes)->linesearch);CHKERRQ(ierr); 2734 2735 ierr = PetscFree((*snes)->kspconvctx);CHKERRQ(ierr); 2736 if ((*snes)->ops->convergeddestroy) { 2737 ierr = (*(*snes)->ops->convergeddestroy)((*snes)->cnvP);CHKERRQ(ierr); 2738 } 2739 if ((*snes)->conv_malloc) { 2740 ierr = PetscFree((*snes)->conv_hist);CHKERRQ(ierr); 2741 ierr = PetscFree((*snes)->conv_hist_its);CHKERRQ(ierr); 2742 } 2743 ierr = SNESMonitorCancel((*snes));CHKERRQ(ierr); 2744 ierr = PetscHeaderDestroy(snes);CHKERRQ(ierr); 2745 PetscFunctionReturn(0); 2746 } 2747 2748 /* ----------- Routines to set solver parameters ---------- */ 2749 2750 #undef __FUNCT__ 2751 #define __FUNCT__ "SNESSetLagPreconditioner" 2752 /*@ 2753 SNESSetLagPreconditioner - Determines when the preconditioner is rebuilt in the nonlinear solve. 2754 2755 Logically Collective on SNES 2756 2757 Input Parameters: 2758 + snes - the SNES context 2759 - lag - -1 indicates NEVER rebuild, 1 means rebuild every time the Jacobian is computed within a single nonlinear solve, 2 means every second time 2760 the Jacobian is built etc. -2 indicates rebuild preconditioner at next chance but then never rebuild after that 2761 2762 Options Database Keys: 2763 . -snes_lag_preconditioner <lag> 2764 2765 Notes: 2766 The default is 1 2767 The preconditioner is ALWAYS built in the first iteration of a nonlinear solve unless lag is -1 2768 If -1 is used before the very first nonlinear solve the preconditioner is still built because there is no previous preconditioner to use 2769 2770 Level: intermediate 2771 2772 .keywords: SNES, nonlinear, set, convergence, tolerances 2773 2774 .seealso: SNESSetTrustRegionTolerance(), SNESGetLagPreconditioner(), SNESSetLagJacobian(), SNESGetLagJacobian() 2775 2776 @*/ 2777 PetscErrorCode SNESSetLagPreconditioner(SNES snes,PetscInt lag) 2778 { 2779 PetscFunctionBegin; 2780 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2781 if (lag < -2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Lag must be -2, -1, 1 or greater"); 2782 if (!lag) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Lag cannot be 0"); 2783 PetscValidLogicalCollectiveInt(snes,lag,2); 2784 snes->lagpreconditioner = lag; 2785 PetscFunctionReturn(0); 2786 } 2787 2788 #undef __FUNCT__ 2789 #define __FUNCT__ "SNESSetGridSequence" 2790 /*@ 2791 SNESSetGridSequence - sets the number of steps of grid sequencing that SNES does 2792 2793 Logically Collective on SNES 2794 2795 Input Parameters: 2796 + snes - the SNES context 2797 - steps - the number of refinements to do, defaults to 0 2798 2799 Options Database Keys: 2800 . -snes_grid_sequence <steps> 2801 2802 Level: intermediate 2803 2804 Notes: 2805 Use SNESGetSolution() to extract the fine grid solution after grid sequencing. 2806 2807 .keywords: SNES, nonlinear, set, convergence, tolerances 2808 2809 .seealso: SNESSetTrustRegionTolerance(), SNESGetLagPreconditioner(), SNESSetLagJacobian(), SNESGetLagJacobian() 2810 2811 @*/ 2812 PetscErrorCode SNESSetGridSequence(SNES snes,PetscInt steps) 2813 { 2814 PetscFunctionBegin; 2815 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2816 PetscValidLogicalCollectiveInt(snes,steps,2); 2817 snes->gridsequence = steps; 2818 PetscFunctionReturn(0); 2819 } 2820 2821 #undef __FUNCT__ 2822 #define __FUNCT__ "SNESGetLagPreconditioner" 2823 /*@ 2824 SNESGetLagPreconditioner - Indicates how often the preconditioner is rebuilt 2825 2826 Not Collective 2827 2828 Input Parameter: 2829 . snes - the SNES context 2830 2831 Output Parameter: 2832 . lag - -1 indicates NEVER rebuild, 1 means rebuild every time the Jacobian is computed within a single nonlinear solve, 2 means every second time 2833 the Jacobian is built etc. -2 indicates rebuild preconditioner at next chance but then never rebuild after that 2834 2835 Options Database Keys: 2836 . -snes_lag_preconditioner <lag> 2837 2838 Notes: 2839 The default is 1 2840 The preconditioner is ALWAYS built in the first iteration of a nonlinear solve unless lag is -1 2841 2842 Level: intermediate 2843 2844 .keywords: SNES, nonlinear, set, convergence, tolerances 2845 2846 .seealso: SNESSetTrustRegionTolerance(), SNESSetLagPreconditioner() 2847 2848 @*/ 2849 PetscErrorCode SNESGetLagPreconditioner(SNES snes,PetscInt *lag) 2850 { 2851 PetscFunctionBegin; 2852 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2853 *lag = snes->lagpreconditioner; 2854 PetscFunctionReturn(0); 2855 } 2856 2857 #undef __FUNCT__ 2858 #define __FUNCT__ "SNESSetLagJacobian" 2859 /*@ 2860 SNESSetLagJacobian - Determines when the Jacobian is rebuilt in the nonlinear solve. See SNESSetLagPreconditioner() for determining how 2861 often the preconditioner is rebuilt. 2862 2863 Logically Collective on SNES 2864 2865 Input Parameters: 2866 + snes - the SNES context 2867 - lag - -1 indicates NEVER rebuild, 1 means rebuild every time the Jacobian is computed within a single nonlinear solve, 2 means every second time 2868 the Jacobian is built etc. -2 means rebuild at next chance but then never again 2869 2870 Options Database Keys: 2871 . -snes_lag_jacobian <lag> 2872 2873 Notes: 2874 The default is 1 2875 The Jacobian is ALWAYS built in the first iteration of a nonlinear solve unless lag is -1 2876 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 2877 at the next Newton step but never again (unless it is reset to another value) 2878 2879 Level: intermediate 2880 2881 .keywords: SNES, nonlinear, set, convergence, tolerances 2882 2883 .seealso: SNESSetTrustRegionTolerance(), SNESGetLagPreconditioner(), SNESSetLagPreconditioner(), SNESGetLagJacobian() 2884 2885 @*/ 2886 PetscErrorCode SNESSetLagJacobian(SNES snes,PetscInt lag) 2887 { 2888 PetscFunctionBegin; 2889 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2890 if (lag < -2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Lag must be -2, -1, 1 or greater"); 2891 if (!lag) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Lag cannot be 0"); 2892 PetscValidLogicalCollectiveInt(snes,lag,2); 2893 snes->lagjacobian = lag; 2894 PetscFunctionReturn(0); 2895 } 2896 2897 #undef __FUNCT__ 2898 #define __FUNCT__ "SNESGetLagJacobian" 2899 /*@ 2900 SNESGetLagJacobian - Indicates how often the Jacobian is rebuilt. See SNESGetLagPreconditioner() to determine when the preconditioner is rebuilt 2901 2902 Not Collective 2903 2904 Input Parameter: 2905 . snes - the SNES context 2906 2907 Output Parameter: 2908 . lag - -1 indicates NEVER rebuild, 1 means rebuild every time the Jacobian is computed within a single nonlinear solve, 2 means every second time 2909 the Jacobian is built etc. 2910 2911 Options Database Keys: 2912 . -snes_lag_jacobian <lag> 2913 2914 Notes: 2915 The default is 1 2916 The jacobian is ALWAYS built in the first iteration of a nonlinear solve unless lag is -1 2917 2918 Level: intermediate 2919 2920 .keywords: SNES, nonlinear, set, convergence, tolerances 2921 2922 .seealso: SNESSetTrustRegionTolerance(), SNESSetLagJacobian(), SNESSetLagPreconditioner(), SNESGetLagPreconditioner() 2923 2924 @*/ 2925 PetscErrorCode SNESGetLagJacobian(SNES snes,PetscInt *lag) 2926 { 2927 PetscFunctionBegin; 2928 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2929 *lag = snes->lagjacobian; 2930 PetscFunctionReturn(0); 2931 } 2932 2933 #undef __FUNCT__ 2934 #define __FUNCT__ "SNESSetTolerances" 2935 /*@ 2936 SNESSetTolerances - Sets various parameters used in convergence tests. 2937 2938 Logically Collective on SNES 2939 2940 Input Parameters: 2941 + snes - the SNES context 2942 . abstol - absolute convergence tolerance 2943 . rtol - relative convergence tolerance 2944 . stol - convergence tolerance in terms of the norm of the change in the solution between steps, || delta x || < stol*|| x || 2945 . maxit - maximum number of iterations 2946 - maxf - maximum number of function evaluations 2947 2948 Options Database Keys: 2949 + -snes_atol <abstol> - Sets abstol 2950 . -snes_rtol <rtol> - Sets rtol 2951 . -snes_stol <stol> - Sets stol 2952 . -snes_max_it <maxit> - Sets maxit 2953 - -snes_max_funcs <maxf> - Sets maxf 2954 2955 Notes: 2956 The default maximum number of iterations is 50. 2957 The default maximum number of function evaluations is 1000. 2958 2959 Level: intermediate 2960 2961 .keywords: SNES, nonlinear, set, convergence, tolerances 2962 2963 .seealso: SNESSetTrustRegionTolerance() 2964 @*/ 2965 PetscErrorCode SNESSetTolerances(SNES snes,PetscReal abstol,PetscReal rtol,PetscReal stol,PetscInt maxit,PetscInt maxf) 2966 { 2967 PetscFunctionBegin; 2968 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2969 PetscValidLogicalCollectiveReal(snes,abstol,2); 2970 PetscValidLogicalCollectiveReal(snes,rtol,3); 2971 PetscValidLogicalCollectiveReal(snes,stol,4); 2972 PetscValidLogicalCollectiveInt(snes,maxit,5); 2973 PetscValidLogicalCollectiveInt(snes,maxf,6); 2974 2975 if (abstol != PETSC_DEFAULT) { 2976 if (abstol < 0.0) SETERRQ1(((PetscObject)snes)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Absolute tolerance %G must be non-negative",abstol); 2977 snes->abstol = abstol; 2978 } 2979 if (rtol != PETSC_DEFAULT) { 2980 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); 2981 snes->rtol = rtol; 2982 } 2983 if (stol != PETSC_DEFAULT) { 2984 if (stol < 0.0) SETERRQ1(((PetscObject)snes)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Step tolerance %G must be non-negative",stol); 2985 snes->stol = stol; 2986 } 2987 if (maxit != PETSC_DEFAULT) { 2988 if (maxit < 0) SETERRQ1(((PetscObject)snes)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Maximum number of iterations %D must be non-negative",maxit); 2989 snes->max_its = maxit; 2990 } 2991 if (maxf != PETSC_DEFAULT) { 2992 if (maxf < 0) SETERRQ1(((PetscObject)snes)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Maximum number of function evaluations %D must be non-negative",maxf); 2993 snes->max_funcs = maxf; 2994 } 2995 snes->tolerancesset = PETSC_TRUE; 2996 PetscFunctionReturn(0); 2997 } 2998 2999 #undef __FUNCT__ 3000 #define __FUNCT__ "SNESGetTolerances" 3001 /*@ 3002 SNESGetTolerances - Gets various parameters used in convergence tests. 3003 3004 Not Collective 3005 3006 Input Parameters: 3007 + snes - the SNES context 3008 . atol - absolute convergence tolerance 3009 . rtol - relative convergence tolerance 3010 . stol - convergence tolerance in terms of the norm 3011 of the change in the solution between steps 3012 . maxit - maximum number of iterations 3013 - maxf - maximum number of function evaluations 3014 3015 Notes: 3016 The user can specify PETSC_NULL for any parameter that is not needed. 3017 3018 Level: intermediate 3019 3020 .keywords: SNES, nonlinear, get, convergence, tolerances 3021 3022 .seealso: SNESSetTolerances() 3023 @*/ 3024 PetscErrorCode SNESGetTolerances(SNES snes,PetscReal *atol,PetscReal *rtol,PetscReal *stol,PetscInt *maxit,PetscInt *maxf) 3025 { 3026 PetscFunctionBegin; 3027 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3028 if (atol) *atol = snes->abstol; 3029 if (rtol) *rtol = snes->rtol; 3030 if (stol) *stol = snes->stol; 3031 if (maxit) *maxit = snes->max_its; 3032 if (maxf) *maxf = snes->max_funcs; 3033 PetscFunctionReturn(0); 3034 } 3035 3036 #undef __FUNCT__ 3037 #define __FUNCT__ "SNESSetTrustRegionTolerance" 3038 /*@ 3039 SNESSetTrustRegionTolerance - Sets the trust region parameter tolerance. 3040 3041 Logically Collective on SNES 3042 3043 Input Parameters: 3044 + snes - the SNES context 3045 - tol - tolerance 3046 3047 Options Database Key: 3048 . -snes_trtol <tol> - Sets tol 3049 3050 Level: intermediate 3051 3052 .keywords: SNES, nonlinear, set, trust region, tolerance 3053 3054 .seealso: SNESSetTolerances() 3055 @*/ 3056 PetscErrorCode SNESSetTrustRegionTolerance(SNES snes,PetscReal tol) 3057 { 3058 PetscFunctionBegin; 3059 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3060 PetscValidLogicalCollectiveReal(snes,tol,2); 3061 snes->deltatol = tol; 3062 PetscFunctionReturn(0); 3063 } 3064 3065 /* 3066 Duplicate the lg monitors for SNES from KSP; for some reason with 3067 dynamic libraries things don't work under Sun4 if we just use 3068 macros instead of functions 3069 */ 3070 #undef __FUNCT__ 3071 #define __FUNCT__ "SNESMonitorLGResidualNorm" 3072 PetscErrorCode SNESMonitorLGResidualNorm(SNES snes,PetscInt it,PetscReal norm,void *ctx) 3073 { 3074 PetscErrorCode ierr; 3075 3076 PetscFunctionBegin; 3077 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3078 ierr = KSPMonitorLGResidualNorm((KSP)snes,it,norm,ctx);CHKERRQ(ierr); 3079 PetscFunctionReturn(0); 3080 } 3081 3082 #undef __FUNCT__ 3083 #define __FUNCT__ "SNESMonitorLGCreate" 3084 PetscErrorCode SNESMonitorLGCreate(const char host[],const char label[],int x,int y,int m,int n,PetscDrawLG *draw) 3085 { 3086 PetscErrorCode ierr; 3087 3088 PetscFunctionBegin; 3089 ierr = KSPMonitorLGResidualNormCreate(host,label,x,y,m,n,draw);CHKERRQ(ierr); 3090 PetscFunctionReturn(0); 3091 } 3092 3093 #undef __FUNCT__ 3094 #define __FUNCT__ "SNESMonitorLGDestroy" 3095 PetscErrorCode SNESMonitorLGDestroy(PetscDrawLG *draw) 3096 { 3097 PetscErrorCode ierr; 3098 3099 PetscFunctionBegin; 3100 ierr = KSPMonitorLGResidualNormDestroy(draw);CHKERRQ(ierr); 3101 PetscFunctionReturn(0); 3102 } 3103 3104 extern PetscErrorCode SNESMonitorRange_Private(SNES,PetscInt,PetscReal*); 3105 #undef __FUNCT__ 3106 #define __FUNCT__ "SNESMonitorLGRange" 3107 PetscErrorCode SNESMonitorLGRange(SNES snes,PetscInt n,PetscReal rnorm,void *monctx) 3108 { 3109 PetscDrawLG lg; 3110 PetscErrorCode ierr; 3111 PetscReal x,y,per; 3112 PetscViewer v = (PetscViewer)monctx; 3113 static PetscReal prev; /* should be in the context */ 3114 PetscDraw draw; 3115 3116 PetscFunctionBegin; 3117 ierr = PetscViewerDrawGetDrawLG(v,0,&lg);CHKERRQ(ierr); 3118 if (!n) {ierr = PetscDrawLGReset(lg);CHKERRQ(ierr);} 3119 ierr = PetscDrawLGGetDraw(lg,&draw);CHKERRQ(ierr); 3120 ierr = PetscDrawSetTitle(draw,"Residual norm");CHKERRQ(ierr); 3121 x = (PetscReal) n; 3122 if (rnorm > 0.0) y = log10(rnorm); else y = -15.0; 3123 ierr = PetscDrawLGAddPoint(lg,&x,&y);CHKERRQ(ierr); 3124 if (n < 20 || !(n % 5)) { 3125 ierr = PetscDrawLGDraw(lg);CHKERRQ(ierr); 3126 } 3127 3128 ierr = PetscViewerDrawGetDrawLG(v,1,&lg);CHKERRQ(ierr); 3129 if (!n) {ierr = PetscDrawLGReset(lg);CHKERRQ(ierr);} 3130 ierr = PetscDrawLGGetDraw(lg,&draw);CHKERRQ(ierr); 3131 ierr = PetscDrawSetTitle(draw,"% elemts > .2*max elemt");CHKERRQ(ierr); 3132 ierr = SNESMonitorRange_Private(snes,n,&per);CHKERRQ(ierr); 3133 x = (PetscReal) n; 3134 y = 100.0*per; 3135 ierr = PetscDrawLGAddPoint(lg,&x,&y);CHKERRQ(ierr); 3136 if (n < 20 || !(n % 5)) { 3137 ierr = PetscDrawLGDraw(lg);CHKERRQ(ierr); 3138 } 3139 3140 ierr = PetscViewerDrawGetDrawLG(v,2,&lg);CHKERRQ(ierr); 3141 if (!n) {prev = rnorm;ierr = PetscDrawLGReset(lg);CHKERRQ(ierr);} 3142 ierr = PetscDrawLGGetDraw(lg,&draw);CHKERRQ(ierr); 3143 ierr = PetscDrawSetTitle(draw,"(norm -oldnorm)/oldnorm");CHKERRQ(ierr); 3144 x = (PetscReal) n; 3145 y = (prev - rnorm)/prev; 3146 ierr = PetscDrawLGAddPoint(lg,&x,&y);CHKERRQ(ierr); 3147 if (n < 20 || !(n % 5)) { 3148 ierr = PetscDrawLGDraw(lg);CHKERRQ(ierr); 3149 } 3150 3151 ierr = PetscViewerDrawGetDrawLG(v,3,&lg);CHKERRQ(ierr); 3152 if (!n) {ierr = PetscDrawLGReset(lg);CHKERRQ(ierr);} 3153 ierr = PetscDrawLGGetDraw(lg,&draw);CHKERRQ(ierr); 3154 ierr = PetscDrawSetTitle(draw,"(norm -oldnorm)/oldnorm*(% > .2 max)");CHKERRQ(ierr); 3155 x = (PetscReal) n; 3156 y = (prev - rnorm)/(prev*per); 3157 if (n > 2) { /*skip initial crazy value */ 3158 ierr = PetscDrawLGAddPoint(lg,&x,&y);CHKERRQ(ierr); 3159 } 3160 if (n < 20 || !(n % 5)) { 3161 ierr = PetscDrawLGDraw(lg);CHKERRQ(ierr); 3162 } 3163 prev = rnorm; 3164 PetscFunctionReturn(0); 3165 } 3166 3167 #undef __FUNCT__ 3168 #define __FUNCT__ "SNESMonitor" 3169 /*@ 3170 SNESMonitor - runs the user provided monitor routines, if they exist 3171 3172 Collective on SNES 3173 3174 Input Parameters: 3175 + snes - nonlinear solver context obtained from SNESCreate() 3176 . iter - iteration number 3177 - rnorm - relative norm of the residual 3178 3179 Notes: 3180 This routine is called by the SNES implementations. 3181 It does not typically need to be called by the user. 3182 3183 Level: developer 3184 3185 .seealso: SNESMonitorSet() 3186 @*/ 3187 PetscErrorCode SNESMonitor(SNES snes,PetscInt iter,PetscReal rnorm) 3188 { 3189 PetscErrorCode ierr; 3190 PetscInt i,n = snes->numbermonitors; 3191 3192 PetscFunctionBegin; 3193 for (i=0; i<n; i++) { 3194 ierr = (*snes->monitor[i])(snes,iter,rnorm,snes->monitorcontext[i]);CHKERRQ(ierr); 3195 } 3196 PetscFunctionReturn(0); 3197 } 3198 3199 /* ------------ Routines to set performance monitoring options ----------- */ 3200 3201 #undef __FUNCT__ 3202 #define __FUNCT__ "SNESMonitorSet" 3203 /*@C 3204 SNESMonitorSet - Sets an ADDITIONAL function that is to be used at every 3205 iteration of the nonlinear solver to display the iteration's 3206 progress. 3207 3208 Logically Collective on SNES 3209 3210 Input Parameters: 3211 + snes - the SNES context 3212 . func - monitoring routine 3213 . mctx - [optional] user-defined context for private data for the 3214 monitor routine (use PETSC_NULL if no context is desired) 3215 - monitordestroy - [optional] routine that frees monitor context 3216 (may be PETSC_NULL) 3217 3218 Calling sequence of func: 3219 $ int func(SNES snes,PetscInt its, PetscReal norm,void *mctx) 3220 3221 + snes - the SNES context 3222 . its - iteration number 3223 . norm - 2-norm function value (may be estimated) 3224 - mctx - [optional] monitoring context 3225 3226 Options Database Keys: 3227 + -snes_monitor - sets SNESMonitorDefault() 3228 . -snes_monitor_lg_residualnorm - sets line graph monitor, 3229 uses SNESMonitorLGCreate() 3230 - -snes_monitor_cancel - cancels all monitors that have 3231 been hardwired into a code by 3232 calls to SNESMonitorSet(), but 3233 does not cancel those set via 3234 the options database. 3235 3236 Notes: 3237 Several different monitoring routines may be set by calling 3238 SNESMonitorSet() multiple times; all will be called in the 3239 order in which they were set. 3240 3241 Fortran notes: Only a single monitor function can be set for each SNES object 3242 3243 Level: intermediate 3244 3245 .keywords: SNES, nonlinear, set, monitor 3246 3247 .seealso: SNESMonitorDefault(), SNESMonitorCancel() 3248 @*/ 3249 PetscErrorCode SNESMonitorSet(SNES snes,PetscErrorCode (*monitor)(SNES,PetscInt,PetscReal,void*),void *mctx,PetscErrorCode (*monitordestroy)(void**)) 3250 { 3251 PetscInt i; 3252 PetscErrorCode ierr; 3253 3254 PetscFunctionBegin; 3255 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3256 if (snes->numbermonitors >= MAXSNESMONITORS) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Too many monitors set"); 3257 for (i=0; i<snes->numbermonitors;i++) { 3258 if (monitor == snes->monitor[i] && monitordestroy == snes->monitordestroy[i] && mctx == snes->monitorcontext[i]) { 3259 if (monitordestroy) { 3260 ierr = (*monitordestroy)(&mctx);CHKERRQ(ierr); 3261 } 3262 PetscFunctionReturn(0); 3263 } 3264 } 3265 snes->monitor[snes->numbermonitors] = monitor; 3266 snes->monitordestroy[snes->numbermonitors] = monitordestroy; 3267 snes->monitorcontext[snes->numbermonitors++] = (void*)mctx; 3268 PetscFunctionReturn(0); 3269 } 3270 3271 #undef __FUNCT__ 3272 #define __FUNCT__ "SNESMonitorCancel" 3273 /*@C 3274 SNESMonitorCancel - Clears all the monitor functions for a SNES object. 3275 3276 Logically Collective on SNES 3277 3278 Input Parameters: 3279 . snes - the SNES context 3280 3281 Options Database Key: 3282 . -snes_monitor_cancel - cancels all monitors that have been hardwired 3283 into a code by calls to SNESMonitorSet(), but does not cancel those 3284 set via the options database 3285 3286 Notes: 3287 There is no way to clear one specific monitor from a SNES object. 3288 3289 Level: intermediate 3290 3291 .keywords: SNES, nonlinear, set, monitor 3292 3293 .seealso: SNESMonitorDefault(), SNESMonitorSet() 3294 @*/ 3295 PetscErrorCode SNESMonitorCancel(SNES snes) 3296 { 3297 PetscErrorCode ierr; 3298 PetscInt i; 3299 3300 PetscFunctionBegin; 3301 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3302 for (i=0; i<snes->numbermonitors; i++) { 3303 if (snes->monitordestroy[i]) { 3304 ierr = (*snes->monitordestroy[i])(&snes->monitorcontext[i]);CHKERRQ(ierr); 3305 } 3306 } 3307 snes->numbermonitors = 0; 3308 PetscFunctionReturn(0); 3309 } 3310 3311 #undef __FUNCT__ 3312 #define __FUNCT__ "SNESSetConvergenceTest" 3313 /*@C 3314 SNESSetConvergenceTest - Sets the function that is to be used 3315 to test for convergence of the nonlinear iterative solution. 3316 3317 Logically Collective on SNES 3318 3319 Input Parameters: 3320 + snes - the SNES context 3321 . func - routine to test for convergence 3322 . cctx - [optional] context for private data for the convergence routine (may be PETSC_NULL) 3323 - destroy - [optional] destructor for the context (may be PETSC_NULL; PETSC_NULL_FUNCTION in Fortran) 3324 3325 Calling sequence of func: 3326 $ PetscErrorCode func (SNES snes,PetscInt it,PetscReal xnorm,PetscReal gnorm,PetscReal f,SNESConvergedReason *reason,void *cctx) 3327 3328 + snes - the SNES context 3329 . it - current iteration (0 is the first and is before any Newton step) 3330 . cctx - [optional] convergence context 3331 . reason - reason for convergence/divergence 3332 . xnorm - 2-norm of current iterate 3333 . gnorm - 2-norm of current step 3334 - f - 2-norm of function 3335 3336 Level: advanced 3337 3338 .keywords: SNES, nonlinear, set, convergence, test 3339 3340 .seealso: SNESDefaultConverged(), SNESSkipConverged() 3341 @*/ 3342 PetscErrorCode SNESSetConvergenceTest(SNES snes,PetscErrorCode (*func)(SNES,PetscInt,PetscReal,PetscReal,PetscReal,SNESConvergedReason*,void*),void *cctx,PetscErrorCode (*destroy)(void*)) 3343 { 3344 PetscErrorCode ierr; 3345 3346 PetscFunctionBegin; 3347 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3348 if (!func) func = SNESSkipConverged; 3349 if (snes->ops->convergeddestroy) { 3350 ierr = (*snes->ops->convergeddestroy)(snes->cnvP);CHKERRQ(ierr); 3351 } 3352 snes->ops->converged = func; 3353 snes->ops->convergeddestroy = destroy; 3354 snes->cnvP = cctx; 3355 PetscFunctionReturn(0); 3356 } 3357 3358 #undef __FUNCT__ 3359 #define __FUNCT__ "SNESGetConvergedReason" 3360 /*@ 3361 SNESGetConvergedReason - Gets the reason the SNES iteration was stopped. 3362 3363 Not Collective 3364 3365 Input Parameter: 3366 . snes - the SNES context 3367 3368 Output Parameter: 3369 . reason - negative value indicates diverged, positive value converged, see SNESConvergedReason or the 3370 manual pages for the individual convergence tests for complete lists 3371 3372 Level: intermediate 3373 3374 Notes: Can only be called after the call the SNESSolve() is complete. 3375 3376 .keywords: SNES, nonlinear, set, convergence, test 3377 3378 .seealso: SNESSetConvergenceTest(), SNESConvergedReason 3379 @*/ 3380 PetscErrorCode SNESGetConvergedReason(SNES snes,SNESConvergedReason *reason) 3381 { 3382 PetscFunctionBegin; 3383 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3384 PetscValidPointer(reason,2); 3385 *reason = snes->reason; 3386 PetscFunctionReturn(0); 3387 } 3388 3389 #undef __FUNCT__ 3390 #define __FUNCT__ "SNESSetConvergenceHistory" 3391 /*@ 3392 SNESSetConvergenceHistory - Sets the array used to hold the convergence history. 3393 3394 Logically Collective on SNES 3395 3396 Input Parameters: 3397 + snes - iterative context obtained from SNESCreate() 3398 . a - array to hold history, this array will contain the function norms computed at each step 3399 . its - integer array holds the number of linear iterations for each solve. 3400 . na - size of a and its 3401 - reset - PETSC_TRUE indicates each new nonlinear solve resets the history counter to zero, 3402 else it continues storing new values for new nonlinear solves after the old ones 3403 3404 Notes: 3405 If 'a' and 'its' are PETSC_NULL then space is allocated for the history. If 'na' PETSC_DECIDE or PETSC_DEFAULT then a 3406 default array of length 10000 is allocated. 3407 3408 This routine is useful, e.g., when running a code for purposes 3409 of accurate performance monitoring, when no I/O should be done 3410 during the section of code that is being timed. 3411 3412 Level: intermediate 3413 3414 .keywords: SNES, set, convergence, history 3415 3416 .seealso: SNESGetConvergenceHistory() 3417 3418 @*/ 3419 PetscErrorCode SNESSetConvergenceHistory(SNES snes,PetscReal a[],PetscInt its[],PetscInt na,PetscBool reset) 3420 { 3421 PetscErrorCode ierr; 3422 3423 PetscFunctionBegin; 3424 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3425 if (na) PetscValidScalarPointer(a,2); 3426 if (its) PetscValidIntPointer(its,3); 3427 if (na == PETSC_DECIDE || na == PETSC_DEFAULT || !a) { 3428 if (na == PETSC_DECIDE || na == PETSC_DEFAULT) na = 1000; 3429 ierr = PetscMalloc(na*sizeof(PetscReal),&a);CHKERRQ(ierr); 3430 ierr = PetscMalloc(na*sizeof(PetscInt),&its);CHKERRQ(ierr); 3431 snes->conv_malloc = PETSC_TRUE; 3432 } 3433 snes->conv_hist = a; 3434 snes->conv_hist_its = its; 3435 snes->conv_hist_max = na; 3436 snes->conv_hist_len = 0; 3437 snes->conv_hist_reset = reset; 3438 PetscFunctionReturn(0); 3439 } 3440 3441 #if defined(PETSC_HAVE_MATLAB_ENGINE) 3442 #include <engine.h> /* MATLAB include file */ 3443 #include <mex.h> /* MATLAB include file */ 3444 EXTERN_C_BEGIN 3445 #undef __FUNCT__ 3446 #define __FUNCT__ "SNESGetConvergenceHistoryMatlab" 3447 mxArray *SNESGetConvergenceHistoryMatlab(SNES snes) 3448 { 3449 mxArray *mat; 3450 PetscInt i; 3451 PetscReal *ar; 3452 3453 PetscFunctionBegin; 3454 mat = mxCreateDoubleMatrix(snes->conv_hist_len,1,mxREAL); 3455 ar = (PetscReal*) mxGetData(mat); 3456 for (i=0; i<snes->conv_hist_len; i++) { 3457 ar[i] = snes->conv_hist[i]; 3458 } 3459 PetscFunctionReturn(mat); 3460 } 3461 EXTERN_C_END 3462 #endif 3463 3464 3465 #undef __FUNCT__ 3466 #define __FUNCT__ "SNESGetConvergenceHistory" 3467 /*@C 3468 SNESGetConvergenceHistory - Gets the array used to hold the convergence history. 3469 3470 Not Collective 3471 3472 Input Parameter: 3473 . snes - iterative context obtained from SNESCreate() 3474 3475 Output Parameters: 3476 . a - array to hold history 3477 . its - integer array holds the number of linear iterations (or 3478 negative if not converged) for each solve. 3479 - na - size of a and its 3480 3481 Notes: 3482 The calling sequence for this routine in Fortran is 3483 $ call SNESGetConvergenceHistory(SNES snes, integer na, integer ierr) 3484 3485 This routine is useful, e.g., when running a code for purposes 3486 of accurate performance monitoring, when no I/O should be done 3487 during the section of code that is being timed. 3488 3489 Level: intermediate 3490 3491 .keywords: SNES, get, convergence, history 3492 3493 .seealso: SNESSetConvergencHistory() 3494 3495 @*/ 3496 PetscErrorCode SNESGetConvergenceHistory(SNES snes,PetscReal *a[],PetscInt *its[],PetscInt *na) 3497 { 3498 PetscFunctionBegin; 3499 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3500 if (a) *a = snes->conv_hist; 3501 if (its) *its = snes->conv_hist_its; 3502 if (na) *na = snes->conv_hist_len; 3503 PetscFunctionReturn(0); 3504 } 3505 3506 #undef __FUNCT__ 3507 #define __FUNCT__ "SNESSetUpdate" 3508 /*@C 3509 SNESSetUpdate - Sets the general-purpose update function called 3510 at the beginning of every iteration of the nonlinear solve. Specifically 3511 it is called just before the Jacobian is "evaluated". 3512 3513 Logically Collective on SNES 3514 3515 Input Parameters: 3516 . snes - The nonlinear solver context 3517 . func - The function 3518 3519 Calling sequence of func: 3520 . func (SNES snes, PetscInt step); 3521 3522 . step - The current step of the iteration 3523 3524 Level: advanced 3525 3526 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() 3527 This is not used by most users. 3528 3529 .keywords: SNES, update 3530 3531 .seealso SNESDefaultUpdate(), SNESSetJacobian(), SNESSolve() 3532 @*/ 3533 PetscErrorCode SNESSetUpdate(SNES snes, PetscErrorCode (*func)(SNES, PetscInt)) 3534 { 3535 PetscFunctionBegin; 3536 PetscValidHeaderSpecific(snes, SNES_CLASSID,1); 3537 snes->ops->update = func; 3538 PetscFunctionReturn(0); 3539 } 3540 3541 #undef __FUNCT__ 3542 #define __FUNCT__ "SNESDefaultUpdate" 3543 /*@ 3544 SNESDefaultUpdate - The default update function which does nothing. 3545 3546 Not collective 3547 3548 Input Parameters: 3549 . snes - The nonlinear solver context 3550 . step - The current step of the iteration 3551 3552 Level: intermediate 3553 3554 .keywords: SNES, update 3555 .seealso SNESSetUpdate(), SNESDefaultRhsBC(), SNESDefaultShortolutionBC() 3556 @*/ 3557 PetscErrorCode SNESDefaultUpdate(SNES snes, PetscInt step) 3558 { 3559 PetscFunctionBegin; 3560 PetscFunctionReturn(0); 3561 } 3562 3563 #undef __FUNCT__ 3564 #define __FUNCT__ "SNESScaleStep_Private" 3565 /* 3566 SNESScaleStep_Private - Scales a step so that its length is less than the 3567 positive parameter delta. 3568 3569 Input Parameters: 3570 + snes - the SNES context 3571 . y - approximate solution of linear system 3572 . fnorm - 2-norm of current function 3573 - delta - trust region size 3574 3575 Output Parameters: 3576 + gpnorm - predicted function norm at the new point, assuming local 3577 linearization. The value is zero if the step lies within the trust 3578 region, and exceeds zero otherwise. 3579 - ynorm - 2-norm of the step 3580 3581 Note: 3582 For non-trust region methods such as SNESNEWTONLS, the parameter delta 3583 is set to be the maximum allowable step size. 3584 3585 .keywords: SNES, nonlinear, scale, step 3586 */ 3587 PetscErrorCode SNESScaleStep_Private(SNES snes,Vec y,PetscReal *fnorm,PetscReal *delta,PetscReal *gpnorm,PetscReal *ynorm) 3588 { 3589 PetscReal nrm; 3590 PetscScalar cnorm; 3591 PetscErrorCode ierr; 3592 3593 PetscFunctionBegin; 3594 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3595 PetscValidHeaderSpecific(y,VEC_CLASSID,2); 3596 PetscCheckSameComm(snes,1,y,2); 3597 3598 ierr = VecNorm(y,NORM_2,&nrm);CHKERRQ(ierr); 3599 if (nrm > *delta) { 3600 nrm = *delta/nrm; 3601 *gpnorm = (1.0 - nrm)*(*fnorm); 3602 cnorm = nrm; 3603 ierr = VecScale(y,cnorm);CHKERRQ(ierr); 3604 *ynorm = *delta; 3605 } else { 3606 *gpnorm = 0.0; 3607 *ynorm = nrm; 3608 } 3609 PetscFunctionReturn(0); 3610 } 3611 3612 #undef __FUNCT__ 3613 #define __FUNCT__ "SNESSolve" 3614 /*@C 3615 SNESSolve - Solves a nonlinear system F(x) = b. 3616 Call SNESSolve() after calling SNESCreate() and optional routines of the form SNESSetXXX(). 3617 3618 Collective on SNES 3619 3620 Input Parameters: 3621 + snes - the SNES context 3622 . b - the constant part of the equation F(x) = b, or PETSC_NULL to use zero. 3623 - x - the solution vector. 3624 3625 Notes: 3626 The user should initialize the vector,x, with the initial guess 3627 for the nonlinear solve prior to calling SNESSolve. In particular, 3628 to employ an initial guess of zero, the user should explicitly set 3629 this vector to zero by calling VecSet(). 3630 3631 Level: beginner 3632 3633 .keywords: SNES, nonlinear, solve 3634 3635 .seealso: SNESCreate(), SNESDestroy(), SNESSetFunction(), SNESSetJacobian(), SNESSetGridSequence(), SNESGetSolution() 3636 @*/ 3637 PetscErrorCode SNESSolve(SNES snes,Vec b,Vec x) 3638 { 3639 PetscErrorCode ierr; 3640 PetscBool flg; 3641 PetscViewer viewer; 3642 PetscInt grid; 3643 Vec xcreated = PETSC_NULL; 3644 DM dm; 3645 PetscViewerFormat format; 3646 3647 PetscFunctionBegin; 3648 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3649 if (x) PetscValidHeaderSpecific(x,VEC_CLASSID,3); 3650 if (x) PetscCheckSameComm(snes,1,x,3); 3651 if (b) PetscValidHeaderSpecific(b,VEC_CLASSID,2); 3652 if (b) PetscCheckSameComm(snes,1,b,2); 3653 3654 if (!x) { 3655 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 3656 ierr = DMCreateGlobalVector(dm,&xcreated);CHKERRQ(ierr); 3657 x = xcreated; 3658 } 3659 3660 for (grid=0; grid<snes->gridsequence; grid++) {ierr = PetscViewerASCIIPushTab(PETSC_VIEWER_STDOUT_(((PetscObject)snes)->comm));CHKERRQ(ierr);} 3661 for (grid=0; grid<snes->gridsequence+1; grid++) { 3662 3663 /* set solution vector */ 3664 if (!grid) {ierr = PetscObjectReference((PetscObject)x);CHKERRQ(ierr);} 3665 ierr = VecDestroy(&snes->vec_sol);CHKERRQ(ierr); 3666 snes->vec_sol = x; 3667 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 3668 3669 /* set affine vector if provided */ 3670 if (b) { ierr = PetscObjectReference((PetscObject)b);CHKERRQ(ierr); } 3671 ierr = VecDestroy(&snes->vec_rhs);CHKERRQ(ierr); 3672 snes->vec_rhs = b; 3673 3674 ierr = SNESSetUp(snes);CHKERRQ(ierr); 3675 3676 if (!grid) { 3677 if (snes->ops->computeinitialguess) { 3678 ierr = (*snes->ops->computeinitialguess)(snes,snes->vec_sol,snes->initialguessP);CHKERRQ(ierr); 3679 } 3680 } 3681 3682 if (snes->conv_hist_reset) snes->conv_hist_len = 0; 3683 snes->nfuncs = 0; snes->linear_its = 0; snes->numFailures = 0; 3684 3685 ierr = PetscLogEventBegin(SNES_Solve,snes,0,0,0);CHKERRQ(ierr); 3686 ierr = (*snes->ops->solve)(snes);CHKERRQ(ierr); 3687 ierr = PetscLogEventEnd(SNES_Solve,snes,0,0,0);CHKERRQ(ierr); 3688 if (snes->domainerror){ 3689 snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; 3690 snes->domainerror = PETSC_FALSE; 3691 } 3692 if (!snes->reason) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Internal error, solver returned without setting converged reason"); 3693 3694 flg = PETSC_FALSE; 3695 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_test_local_min",&flg,PETSC_NULL);CHKERRQ(ierr); 3696 if (flg && !PetscPreLoadingOn) { ierr = SNESTestLocalMin(snes);CHKERRQ(ierr); } 3697 if (snes->printreason) { 3698 ierr = PetscViewerASCIIAddTab(PETSC_VIEWER_STDOUT_(((PetscObject)snes)->comm),((PetscObject)snes)->tablevel);CHKERRQ(ierr); 3699 if (snes->reason > 0) { 3700 ierr = PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_(((PetscObject)snes)->comm),"Nonlinear solve converged due to %s iterations %D\n",SNESConvergedReasons[snes->reason],snes->iter);CHKERRQ(ierr); 3701 } else { 3702 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); 3703 } 3704 ierr = PetscViewerASCIISubtractTab(PETSC_VIEWER_STDOUT_(((PetscObject)snes)->comm),((PetscObject)snes)->tablevel);CHKERRQ(ierr); 3705 } 3706 3707 if (snes->errorifnotconverged && snes->reason < 0) SETERRQ(((PetscObject)snes)->comm,PETSC_ERR_NOT_CONVERGED,"SNESSolve has not converged"); 3708 if (grid < snes->gridsequence) { 3709 DM fine; 3710 Vec xnew; 3711 Mat interp; 3712 3713 ierr = DMRefine(snes->dm,((PetscObject)snes)->comm,&fine);CHKERRQ(ierr); 3714 if (!fine) SETERRQ(((PetscObject)snes)->comm,PETSC_ERR_ARG_INCOMP,"DMRefine() did not perform any refinement, cannot continue grid sequencing"); 3715 ierr = DMCreateInterpolation(snes->dm,fine,&interp,PETSC_NULL);CHKERRQ(ierr); 3716 ierr = DMCreateGlobalVector(fine,&xnew);CHKERRQ(ierr); 3717 ierr = MatInterpolate(interp,x,xnew);CHKERRQ(ierr); 3718 ierr = DMInterpolate(snes->dm,interp,fine);CHKERRQ(ierr); 3719 ierr = MatDestroy(&interp);CHKERRQ(ierr); 3720 x = xnew; 3721 3722 ierr = SNESReset(snes);CHKERRQ(ierr); 3723 ierr = SNESSetDM(snes,fine);CHKERRQ(ierr); 3724 ierr = DMDestroy(&fine);CHKERRQ(ierr); 3725 ierr = PetscViewerASCIIPopTab(PETSC_VIEWER_STDOUT_(((PetscObject)snes)->comm));CHKERRQ(ierr); 3726 } 3727 } 3728 /* monitoring and viewing */ 3729 ierr = PetscOptionsGetViewer(((PetscObject)snes)->comm,((PetscObject)snes)->prefix,"-snes_view",&viewer,&format,&flg);CHKERRQ(ierr); 3730 if (flg && !PetscPreLoadingOn) { 3731 ierr = PetscViewerPushFormat(viewer,format);CHKERRQ(ierr); 3732 ierr = SNESView(snes,viewer);CHKERRQ(ierr); 3733 ierr = PetscViewerPopFormat(viewer);CHKERRQ(ierr); 3734 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 3735 } 3736 ierr = VecViewFromOptions(snes->vec_sol,"-snes_view_solution");CHKERRQ(ierr); 3737 3738 ierr = VecDestroy(&xcreated);CHKERRQ(ierr); 3739 PetscFunctionReturn(0); 3740 } 3741 3742 /* --------- Internal routines for SNES Package --------- */ 3743 3744 #undef __FUNCT__ 3745 #define __FUNCT__ "SNESSetType" 3746 /*@C 3747 SNESSetType - Sets the method for the nonlinear solver. 3748 3749 Collective on SNES 3750 3751 Input Parameters: 3752 + snes - the SNES context 3753 - type - a known method 3754 3755 Options Database Key: 3756 . -snes_type <type> - Sets the method; use -help for a list 3757 of available methods (for instance, newtonls or newtontr) 3758 3759 Notes: 3760 See "petsc/include/petscsnes.h" for available methods (for instance) 3761 + SNESNEWTONLS - Newton's method with line search 3762 (systems of nonlinear equations) 3763 . SNESNEWTONTR - Newton's method with trust region 3764 (systems of nonlinear equations) 3765 3766 Normally, it is best to use the SNESSetFromOptions() command and then 3767 set the SNES solver type from the options database rather than by using 3768 this routine. Using the options database provides the user with 3769 maximum flexibility in evaluating the many nonlinear solvers. 3770 The SNESSetType() routine is provided for those situations where it 3771 is necessary to set the nonlinear solver independently of the command 3772 line or options database. This might be the case, for example, when 3773 the choice of solver changes during the execution of the program, 3774 and the user's application is taking responsibility for choosing the 3775 appropriate method. 3776 3777 Level: intermediate 3778 3779 .keywords: SNES, set, type 3780 3781 .seealso: SNESType, SNESCreate() 3782 3783 @*/ 3784 PetscErrorCode SNESSetType(SNES snes,SNESType type) 3785 { 3786 PetscErrorCode ierr,(*r)(SNES); 3787 PetscBool match; 3788 3789 PetscFunctionBegin; 3790 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3791 PetscValidCharPointer(type,2); 3792 3793 ierr = PetscObjectTypeCompare((PetscObject)snes,type,&match);CHKERRQ(ierr); 3794 if (match) PetscFunctionReturn(0); 3795 3796 ierr = PetscFListFind(SNESList,((PetscObject)snes)->comm,type,PETSC_TRUE,(void (**)(void)) &r);CHKERRQ(ierr); 3797 if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unable to find requested SNES type %s",type); 3798 /* Destroy the previous private SNES context */ 3799 if (snes->ops->destroy) { 3800 ierr = (*(snes)->ops->destroy)(snes);CHKERRQ(ierr); 3801 snes->ops->destroy = PETSC_NULL; 3802 } 3803 /* Reinitialize function pointers in SNESOps structure */ 3804 snes->ops->setup = 0; 3805 snes->ops->solve = 0; 3806 snes->ops->view = 0; 3807 snes->ops->setfromoptions = 0; 3808 snes->ops->destroy = 0; 3809 /* Call the SNESCreate_XXX routine for this particular Nonlinear solver */ 3810 snes->setupcalled = PETSC_FALSE; 3811 ierr = PetscObjectChangeTypeName((PetscObject)snes,type);CHKERRQ(ierr); 3812 ierr = (*r)(snes);CHKERRQ(ierr); 3813 #if defined(PETSC_HAVE_AMS) 3814 if (PetscAMSPublishAll) { 3815 ierr = PetscObjectAMSPublish((PetscObject)snes);CHKERRQ(ierr); 3816 } 3817 #endif 3818 PetscFunctionReturn(0); 3819 } 3820 3821 3822 /* --------------------------------------------------------------------- */ 3823 #undef __FUNCT__ 3824 #define __FUNCT__ "SNESRegisterDestroy" 3825 /*@ 3826 SNESRegisterDestroy - Frees the list of nonlinear solvers that were 3827 registered by SNESRegisterDynamic(). 3828 3829 Not Collective 3830 3831 Level: advanced 3832 3833 .keywords: SNES, nonlinear, register, destroy 3834 3835 .seealso: SNESRegisterAll() 3836 @*/ 3837 PetscErrorCode SNESRegisterDestroy(void) 3838 { 3839 PetscErrorCode ierr; 3840 3841 PetscFunctionBegin; 3842 ierr = PetscFListDestroy(&SNESList);CHKERRQ(ierr); 3843 SNESRegisterAllCalled = PETSC_FALSE; 3844 PetscFunctionReturn(0); 3845 } 3846 3847 #undef __FUNCT__ 3848 #define __FUNCT__ "SNESGetType" 3849 /*@C 3850 SNESGetType - Gets the SNES method type and name (as a string). 3851 3852 Not Collective 3853 3854 Input Parameter: 3855 . snes - nonlinear solver context 3856 3857 Output Parameter: 3858 . type - SNES method (a character string) 3859 3860 Level: intermediate 3861 3862 .keywords: SNES, nonlinear, get, type, name 3863 @*/ 3864 PetscErrorCode SNESGetType(SNES snes,SNESType *type) 3865 { 3866 PetscFunctionBegin; 3867 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3868 PetscValidPointer(type,2); 3869 *type = ((PetscObject)snes)->type_name; 3870 PetscFunctionReturn(0); 3871 } 3872 3873 #undef __FUNCT__ 3874 #define __FUNCT__ "SNESGetSolution" 3875 /*@ 3876 SNESGetSolution - Returns the vector where the approximate solution is 3877 stored. This is the fine grid solution when using SNESSetGridSequence(). 3878 3879 Not Collective, but Vec is parallel if SNES is parallel 3880 3881 Input Parameter: 3882 . snes - the SNES context 3883 3884 Output Parameter: 3885 . x - the solution 3886 3887 Level: intermediate 3888 3889 .keywords: SNES, nonlinear, get, solution 3890 3891 .seealso: SNESGetSolutionUpdate(), SNESGetFunction() 3892 @*/ 3893 PetscErrorCode SNESGetSolution(SNES snes,Vec *x) 3894 { 3895 PetscFunctionBegin; 3896 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3897 PetscValidPointer(x,2); 3898 *x = snes->vec_sol; 3899 PetscFunctionReturn(0); 3900 } 3901 3902 #undef __FUNCT__ 3903 #define __FUNCT__ "SNESGetSolutionUpdate" 3904 /*@ 3905 SNESGetSolutionUpdate - Returns the vector where the solution update is 3906 stored. 3907 3908 Not Collective, but Vec is parallel if SNES is parallel 3909 3910 Input Parameter: 3911 . snes - the SNES context 3912 3913 Output Parameter: 3914 . x - the solution update 3915 3916 Level: advanced 3917 3918 .keywords: SNES, nonlinear, get, solution, update 3919 3920 .seealso: SNESGetSolution(), SNESGetFunction() 3921 @*/ 3922 PetscErrorCode SNESGetSolutionUpdate(SNES snes,Vec *x) 3923 { 3924 PetscFunctionBegin; 3925 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3926 PetscValidPointer(x,2); 3927 *x = snes->vec_sol_update; 3928 PetscFunctionReturn(0); 3929 } 3930 3931 #undef __FUNCT__ 3932 #define __FUNCT__ "SNESGetFunction" 3933 /*@C 3934 SNESGetFunction - Returns the vector where the function is stored. 3935 3936 Not Collective, but Vec is parallel if SNES is parallel. Collective if Vec is requested, but has not been created yet. 3937 3938 Input Parameter: 3939 . snes - the SNES context 3940 3941 Output Parameter: 3942 + r - the function (or PETSC_NULL) 3943 . func - the function (or PETSC_NULL) 3944 - ctx - the function context (or PETSC_NULL) 3945 3946 Level: advanced 3947 3948 .keywords: SNES, nonlinear, get, function 3949 3950 .seealso: SNESSetFunction(), SNESGetSolution() 3951 @*/ 3952 PetscErrorCode SNESGetFunction(SNES snes,Vec *r,PetscErrorCode (**func)(SNES,Vec,Vec,void*),void **ctx) 3953 { 3954 PetscErrorCode ierr; 3955 DM dm; 3956 3957 PetscFunctionBegin; 3958 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3959 if (r) { 3960 if (!snes->vec_func) { 3961 if (snes->vec_rhs) { 3962 ierr = VecDuplicate(snes->vec_rhs,&snes->vec_func);CHKERRQ(ierr); 3963 } else if (snes->vec_sol) { 3964 ierr = VecDuplicate(snes->vec_sol,&snes->vec_func);CHKERRQ(ierr); 3965 } else if (snes->dm) { 3966 ierr = DMCreateGlobalVector(snes->dm,&snes->vec_func);CHKERRQ(ierr); 3967 } 3968 } 3969 *r = snes->vec_func; 3970 } 3971 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 3972 ierr = DMSNESGetFunction(dm,func,ctx);CHKERRQ(ierr); 3973 PetscFunctionReturn(0); 3974 } 3975 3976 /*@C 3977 SNESGetGS - Returns the GS function and context. 3978 3979 Input Parameter: 3980 . snes - the SNES context 3981 3982 Output Parameter: 3983 + gsfunc - the function (or PETSC_NULL) 3984 - ctx - the function context (or PETSC_NULL) 3985 3986 Level: advanced 3987 3988 .keywords: SNES, nonlinear, get, function 3989 3990 .seealso: SNESSetGS(), SNESGetFunction() 3991 @*/ 3992 3993 #undef __FUNCT__ 3994 #define __FUNCT__ "SNESGetGS" 3995 PetscErrorCode SNESGetGS (SNES snes, PetscErrorCode(**func)(SNES, Vec, Vec, void*), void ** ctx) 3996 { 3997 PetscErrorCode ierr; 3998 DM dm; 3999 4000 PetscFunctionBegin; 4001 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4002 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 4003 ierr = DMSNESGetGS(dm,func,ctx);CHKERRQ(ierr); 4004 PetscFunctionReturn(0); 4005 } 4006 4007 #undef __FUNCT__ 4008 #define __FUNCT__ "SNESSetOptionsPrefix" 4009 /*@C 4010 SNESSetOptionsPrefix - Sets the prefix used for searching for all 4011 SNES options in the database. 4012 4013 Logically Collective on SNES 4014 4015 Input Parameter: 4016 + snes - the SNES context 4017 - prefix - the prefix to prepend to all option names 4018 4019 Notes: 4020 A hyphen (-) must NOT be given at the beginning of the prefix name. 4021 The first character of all runtime options is AUTOMATICALLY the hyphen. 4022 4023 Level: advanced 4024 4025 .keywords: SNES, set, options, prefix, database 4026 4027 .seealso: SNESSetFromOptions() 4028 @*/ 4029 PetscErrorCode SNESSetOptionsPrefix(SNES snes,const char prefix[]) 4030 { 4031 PetscErrorCode ierr; 4032 4033 PetscFunctionBegin; 4034 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4035 ierr = PetscObjectSetOptionsPrefix((PetscObject)snes,prefix);CHKERRQ(ierr); 4036 if (!snes->ksp) {ierr = SNESGetKSP(snes,&snes->ksp);CHKERRQ(ierr);} 4037 if (snes->linesearch) { 4038 ierr = SNESGetSNESLineSearch(snes,&snes->linesearch);CHKERRQ(ierr); 4039 ierr = PetscObjectSetOptionsPrefix((PetscObject)snes->linesearch,prefix);CHKERRQ(ierr); 4040 } 4041 ierr = KSPSetOptionsPrefix(snes->ksp,prefix);CHKERRQ(ierr); 4042 PetscFunctionReturn(0); 4043 } 4044 4045 #undef __FUNCT__ 4046 #define __FUNCT__ "SNESAppendOptionsPrefix" 4047 /*@C 4048 SNESAppendOptionsPrefix - Appends to the prefix used for searching for all 4049 SNES options in the database. 4050 4051 Logically Collective on SNES 4052 4053 Input Parameters: 4054 + snes - the SNES context 4055 - prefix - the prefix to prepend to all option names 4056 4057 Notes: 4058 A hyphen (-) must NOT be given at the beginning of the prefix name. 4059 The first character of all runtime options is AUTOMATICALLY the hyphen. 4060 4061 Level: advanced 4062 4063 .keywords: SNES, append, options, prefix, database 4064 4065 .seealso: SNESGetOptionsPrefix() 4066 @*/ 4067 PetscErrorCode SNESAppendOptionsPrefix(SNES snes,const char prefix[]) 4068 { 4069 PetscErrorCode ierr; 4070 4071 PetscFunctionBegin; 4072 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4073 ierr = PetscObjectAppendOptionsPrefix((PetscObject)snes,prefix);CHKERRQ(ierr); 4074 if (!snes->ksp) {ierr = SNESGetKSP(snes,&snes->ksp);CHKERRQ(ierr);} 4075 if (snes->linesearch) { 4076 ierr = SNESGetSNESLineSearch(snes,&snes->linesearch);CHKERRQ(ierr); 4077 ierr = PetscObjectAppendOptionsPrefix((PetscObject)snes->linesearch,prefix);CHKERRQ(ierr); 4078 } 4079 ierr = KSPAppendOptionsPrefix(snes->ksp,prefix);CHKERRQ(ierr); 4080 PetscFunctionReturn(0); 4081 } 4082 4083 #undef __FUNCT__ 4084 #define __FUNCT__ "SNESGetOptionsPrefix" 4085 /*@C 4086 SNESGetOptionsPrefix - Sets the prefix used for searching for all 4087 SNES options in the database. 4088 4089 Not Collective 4090 4091 Input Parameter: 4092 . snes - the SNES context 4093 4094 Output Parameter: 4095 . prefix - pointer to the prefix string used 4096 4097 Notes: On the fortran side, the user should pass in a string 'prefix' of 4098 sufficient length to hold the prefix. 4099 4100 Level: advanced 4101 4102 .keywords: SNES, get, options, prefix, database 4103 4104 .seealso: SNESAppendOptionsPrefix() 4105 @*/ 4106 PetscErrorCode SNESGetOptionsPrefix(SNES snes,const char *prefix[]) 4107 { 4108 PetscErrorCode ierr; 4109 4110 PetscFunctionBegin; 4111 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4112 ierr = PetscObjectGetOptionsPrefix((PetscObject)snes,prefix);CHKERRQ(ierr); 4113 PetscFunctionReturn(0); 4114 } 4115 4116 4117 #undef __FUNCT__ 4118 #define __FUNCT__ "SNESRegister" 4119 /*@C 4120 SNESRegister - See SNESRegisterDynamic() 4121 4122 Level: advanced 4123 @*/ 4124 PetscErrorCode SNESRegister(const char sname[],const char path[],const char name[],PetscErrorCode (*function)(SNES)) 4125 { 4126 char fullname[PETSC_MAX_PATH_LEN]; 4127 PetscErrorCode ierr; 4128 4129 PetscFunctionBegin; 4130 ierr = PetscFListConcat(path,name,fullname);CHKERRQ(ierr); 4131 ierr = PetscFListAdd(&SNESList,sname,fullname,(void (*)(void))function);CHKERRQ(ierr); 4132 PetscFunctionReturn(0); 4133 } 4134 4135 #undef __FUNCT__ 4136 #define __FUNCT__ "SNESTestLocalMin" 4137 PetscErrorCode SNESTestLocalMin(SNES snes) 4138 { 4139 PetscErrorCode ierr; 4140 PetscInt N,i,j; 4141 Vec u,uh,fh; 4142 PetscScalar value; 4143 PetscReal norm; 4144 4145 PetscFunctionBegin; 4146 ierr = SNESGetSolution(snes,&u);CHKERRQ(ierr); 4147 ierr = VecDuplicate(u,&uh);CHKERRQ(ierr); 4148 ierr = VecDuplicate(u,&fh);CHKERRQ(ierr); 4149 4150 /* currently only works for sequential */ 4151 ierr = PetscPrintf(PETSC_COMM_WORLD,"Testing FormFunction() for local min\n");CHKERRQ(ierr); 4152 ierr = VecGetSize(u,&N);CHKERRQ(ierr); 4153 for (i=0; i<N; i++) { 4154 ierr = VecCopy(u,uh);CHKERRQ(ierr); 4155 ierr = PetscPrintf(PETSC_COMM_WORLD,"i = %D\n",i);CHKERRQ(ierr); 4156 for (j=-10; j<11; j++) { 4157 value = PetscSign(j)*exp(PetscAbs(j)-10.0); 4158 ierr = VecSetValue(uh,i,value,ADD_VALUES);CHKERRQ(ierr); 4159 ierr = SNESComputeFunction(snes,uh,fh);CHKERRQ(ierr); 4160 ierr = VecNorm(fh,NORM_2,&norm);CHKERRQ(ierr); 4161 ierr = PetscPrintf(PETSC_COMM_WORLD," j norm %D %18.16e\n",j,norm);CHKERRQ(ierr); 4162 value = -value; 4163 ierr = VecSetValue(uh,i,value,ADD_VALUES);CHKERRQ(ierr); 4164 } 4165 } 4166 ierr = VecDestroy(&uh);CHKERRQ(ierr); 4167 ierr = VecDestroy(&fh);CHKERRQ(ierr); 4168 PetscFunctionReturn(0); 4169 } 4170 4171 #undef __FUNCT__ 4172 #define __FUNCT__ "SNESKSPSetUseEW" 4173 /*@ 4174 SNESKSPSetUseEW - Sets SNES use Eisenstat-Walker method for 4175 computing relative tolerance for linear solvers within an inexact 4176 Newton method. 4177 4178 Logically Collective on SNES 4179 4180 Input Parameters: 4181 + snes - SNES context 4182 - flag - PETSC_TRUE or PETSC_FALSE 4183 4184 Options Database: 4185 + -snes_ksp_ew - use Eisenstat-Walker method for determining linear system convergence 4186 . -snes_ksp_ew_version ver - version of Eisenstat-Walker method 4187 . -snes_ksp_ew_rtol0 <rtol0> - Sets rtol0 4188 . -snes_ksp_ew_rtolmax <rtolmax> - Sets rtolmax 4189 . -snes_ksp_ew_gamma <gamma> - Sets gamma 4190 . -snes_ksp_ew_alpha <alpha> - Sets alpha 4191 . -snes_ksp_ew_alpha2 <alpha2> - Sets alpha2 4192 - -snes_ksp_ew_threshold <threshold> - Sets threshold 4193 4194 Notes: 4195 Currently, the default is to use a constant relative tolerance for 4196 the inner linear solvers. Alternatively, one can use the 4197 Eisenstat-Walker method, where the relative convergence tolerance 4198 is reset at each Newton iteration according progress of the nonlinear 4199 solver. 4200 4201 Level: advanced 4202 4203 Reference: 4204 S. C. Eisenstat and H. F. Walker, "Choosing the forcing terms in an 4205 inexact Newton method", SISC 17 (1), pp.16-32, 1996. 4206 4207 .keywords: SNES, KSP, Eisenstat, Walker, convergence, test, inexact, Newton 4208 4209 .seealso: SNESKSPGetUseEW(), SNESKSPGetParametersEW(), SNESKSPSetParametersEW() 4210 @*/ 4211 PetscErrorCode SNESKSPSetUseEW(SNES snes,PetscBool flag) 4212 { 4213 PetscFunctionBegin; 4214 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4215 PetscValidLogicalCollectiveBool(snes,flag,2); 4216 snes->ksp_ewconv = flag; 4217 PetscFunctionReturn(0); 4218 } 4219 4220 #undef __FUNCT__ 4221 #define __FUNCT__ "SNESKSPGetUseEW" 4222 /*@ 4223 SNESKSPGetUseEW - Gets if SNES is using Eisenstat-Walker method 4224 for computing relative tolerance for linear solvers within an 4225 inexact Newton method. 4226 4227 Not Collective 4228 4229 Input Parameter: 4230 . snes - SNES context 4231 4232 Output Parameter: 4233 . flag - PETSC_TRUE or PETSC_FALSE 4234 4235 Notes: 4236 Currently, the default is to use a constant relative tolerance for 4237 the inner linear solvers. Alternatively, one can use the 4238 Eisenstat-Walker method, where the relative convergence tolerance 4239 is reset at each Newton iteration according progress of the nonlinear 4240 solver. 4241 4242 Level: advanced 4243 4244 Reference: 4245 S. C. Eisenstat and H. F. Walker, "Choosing the forcing terms in an 4246 inexact Newton method", SISC 17 (1), pp.16-32, 1996. 4247 4248 .keywords: SNES, KSP, Eisenstat, Walker, convergence, test, inexact, Newton 4249 4250 .seealso: SNESKSPSetUseEW(), SNESKSPGetParametersEW(), SNESKSPSetParametersEW() 4251 @*/ 4252 PetscErrorCode SNESKSPGetUseEW(SNES snes, PetscBool *flag) 4253 { 4254 PetscFunctionBegin; 4255 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4256 PetscValidPointer(flag,2); 4257 *flag = snes->ksp_ewconv; 4258 PetscFunctionReturn(0); 4259 } 4260 4261 #undef __FUNCT__ 4262 #define __FUNCT__ "SNESKSPSetParametersEW" 4263 /*@ 4264 SNESKSPSetParametersEW - Sets parameters for Eisenstat-Walker 4265 convergence criteria for the linear solvers within an inexact 4266 Newton method. 4267 4268 Logically Collective on SNES 4269 4270 Input Parameters: 4271 + snes - SNES context 4272 . version - version 1, 2 (default is 2) or 3 4273 . rtol_0 - initial relative tolerance (0 <= rtol_0 < 1) 4274 . rtol_max - maximum relative tolerance (0 <= rtol_max < 1) 4275 . gamma - multiplicative factor for version 2 rtol computation 4276 (0 <= gamma2 <= 1) 4277 . alpha - power for version 2 rtol computation (1 < alpha <= 2) 4278 . alpha2 - power for safeguard 4279 - threshold - threshold for imposing safeguard (0 < threshold < 1) 4280 4281 Note: 4282 Version 3 was contributed by Luis Chacon, June 2006. 4283 4284 Use PETSC_DEFAULT to retain the default for any of the parameters. 4285 4286 Level: advanced 4287 4288 Reference: 4289 S. C. Eisenstat and H. F. Walker, "Choosing the forcing terms in an 4290 inexact Newton method", Utah State University Math. Stat. Dept. Res. 4291 Report 6/94/75, June, 1994, to appear in SIAM J. Sci. Comput. 4292 4293 .keywords: SNES, KSP, Eisenstat, Walker, set, parameters 4294 4295 .seealso: SNESKSPSetUseEW(), SNESKSPGetUseEW(), SNESKSPGetParametersEW() 4296 @*/ 4297 PetscErrorCode SNESKSPSetParametersEW(SNES snes,PetscInt version,PetscReal rtol_0,PetscReal rtol_max, 4298 PetscReal gamma,PetscReal alpha,PetscReal alpha2,PetscReal threshold) 4299 { 4300 SNESKSPEW *kctx; 4301 PetscFunctionBegin; 4302 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4303 kctx = (SNESKSPEW*)snes->kspconvctx; 4304 if (!kctx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"No Eisenstat-Walker context existing"); 4305 PetscValidLogicalCollectiveInt(snes,version,2); 4306 PetscValidLogicalCollectiveReal(snes,rtol_0,3); 4307 PetscValidLogicalCollectiveReal(snes,rtol_max,4); 4308 PetscValidLogicalCollectiveReal(snes,gamma,5); 4309 PetscValidLogicalCollectiveReal(snes,alpha,6); 4310 PetscValidLogicalCollectiveReal(snes,alpha2,7); 4311 PetscValidLogicalCollectiveReal(snes,threshold,8); 4312 4313 if (version != PETSC_DEFAULT) kctx->version = version; 4314 if (rtol_0 != PETSC_DEFAULT) kctx->rtol_0 = rtol_0; 4315 if (rtol_max != PETSC_DEFAULT) kctx->rtol_max = rtol_max; 4316 if (gamma != PETSC_DEFAULT) kctx->gamma = gamma; 4317 if (alpha != PETSC_DEFAULT) kctx->alpha = alpha; 4318 if (alpha2 != PETSC_DEFAULT) kctx->alpha2 = alpha2; 4319 if (threshold != PETSC_DEFAULT) kctx->threshold = threshold; 4320 4321 if (kctx->version < 1 || kctx->version > 3) { 4322 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Only versions 1, 2 and 3 are supported: %D",kctx->version); 4323 } 4324 if (kctx->rtol_0 < 0.0 || kctx->rtol_0 >= 1.0) { 4325 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"0.0 <= rtol_0 < 1.0: %G",kctx->rtol_0); 4326 } 4327 if (kctx->rtol_max < 0.0 || kctx->rtol_max >= 1.0) { 4328 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"0.0 <= rtol_max (%G) < 1.0\n",kctx->rtol_max); 4329 } 4330 if (kctx->gamma < 0.0 || kctx->gamma > 1.0) { 4331 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"0.0 <= gamma (%G) <= 1.0\n",kctx->gamma); 4332 } 4333 if (kctx->alpha <= 1.0 || kctx->alpha > 2.0) { 4334 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"1.0 < alpha (%G) <= 2.0\n",kctx->alpha); 4335 } 4336 if (kctx->threshold <= 0.0 || kctx->threshold >= 1.0) { 4337 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"0.0 < threshold (%G) < 1.0\n",kctx->threshold); 4338 } 4339 PetscFunctionReturn(0); 4340 } 4341 4342 #undef __FUNCT__ 4343 #define __FUNCT__ "SNESKSPGetParametersEW" 4344 /*@ 4345 SNESKSPGetParametersEW - Gets parameters for Eisenstat-Walker 4346 convergence criteria for the linear solvers within an inexact 4347 Newton method. 4348 4349 Not Collective 4350 4351 Input Parameters: 4352 snes - SNES context 4353 4354 Output Parameters: 4355 + version - version 1, 2 (default is 2) or 3 4356 . rtol_0 - initial relative tolerance (0 <= rtol_0 < 1) 4357 . rtol_max - maximum relative tolerance (0 <= rtol_max < 1) 4358 . gamma - multiplicative factor for version 2 rtol computation 4359 (0 <= gamma2 <= 1) 4360 . alpha - power for version 2 rtol computation (1 < alpha <= 2) 4361 . alpha2 - power for safeguard 4362 - threshold - threshold for imposing safeguard (0 < threshold < 1) 4363 4364 Level: advanced 4365 4366 .keywords: SNES, KSP, Eisenstat, Walker, get, parameters 4367 4368 .seealso: SNESKSPSetUseEW(), SNESKSPGetUseEW(), SNESKSPSetParametersEW() 4369 @*/ 4370 PetscErrorCode SNESKSPGetParametersEW(SNES snes,PetscInt *version,PetscReal *rtol_0,PetscReal *rtol_max, 4371 PetscReal *gamma,PetscReal *alpha,PetscReal *alpha2,PetscReal *threshold) 4372 { 4373 SNESKSPEW *kctx; 4374 PetscFunctionBegin; 4375 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4376 kctx = (SNESKSPEW*)snes->kspconvctx; 4377 if (!kctx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"No Eisenstat-Walker context existing"); 4378 if (version) *version = kctx->version; 4379 if (rtol_0) *rtol_0 = kctx->rtol_0; 4380 if (rtol_max) *rtol_max = kctx->rtol_max; 4381 if (gamma) *gamma = kctx->gamma; 4382 if (alpha) *alpha = kctx->alpha; 4383 if (alpha2) *alpha2 = kctx->alpha2; 4384 if (threshold) *threshold = kctx->threshold; 4385 PetscFunctionReturn(0); 4386 } 4387 4388 #undef __FUNCT__ 4389 #define __FUNCT__ "SNESKSPEW_PreSolve" 4390 static PetscErrorCode SNESKSPEW_PreSolve(SNES snes, KSP ksp, Vec b, Vec x) 4391 { 4392 PetscErrorCode ierr; 4393 SNESKSPEW *kctx = (SNESKSPEW*)snes->kspconvctx; 4394 PetscReal rtol=PETSC_DEFAULT,stol; 4395 4396 PetscFunctionBegin; 4397 if (!kctx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"No Eisenstat-Walker context exists"); 4398 if (!snes->iter) { /* first time in, so use the original user rtol */ 4399 rtol = kctx->rtol_0; 4400 } else { 4401 if (kctx->version == 1) { 4402 rtol = (snes->norm - kctx->lresid_last)/kctx->norm_last; 4403 if (rtol < 0.0) rtol = -rtol; 4404 stol = pow(kctx->rtol_last,kctx->alpha2); 4405 if (stol > kctx->threshold) rtol = PetscMax(rtol,stol); 4406 } else if (kctx->version == 2) { 4407 rtol = kctx->gamma * pow(snes->norm/kctx->norm_last,kctx->alpha); 4408 stol = kctx->gamma * pow(kctx->rtol_last,kctx->alpha); 4409 if (stol > kctx->threshold) rtol = PetscMax(rtol,stol); 4410 } else if (kctx->version == 3) {/* contributed by Luis Chacon, June 2006. */ 4411 rtol = kctx->gamma * pow(snes->norm/kctx->norm_last,kctx->alpha); 4412 /* safeguard: avoid sharp decrease of rtol */ 4413 stol = kctx->gamma*pow(kctx->rtol_last,kctx->alpha); 4414 stol = PetscMax(rtol,stol); 4415 rtol = PetscMin(kctx->rtol_0,stol); 4416 /* safeguard: avoid oversolving */ 4417 stol = kctx->gamma*(snes->ttol)/snes->norm; 4418 stol = PetscMax(rtol,stol); 4419 rtol = PetscMin(kctx->rtol_0,stol); 4420 } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Only versions 1, 2 or 3 are supported: %D",kctx->version); 4421 } 4422 /* safeguard: avoid rtol greater than one */ 4423 rtol = PetscMin(rtol,kctx->rtol_max); 4424 ierr = KSPSetTolerances(ksp,rtol,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);CHKERRQ(ierr); 4425 ierr = PetscInfo3(snes,"iter %D, Eisenstat-Walker (version %D) KSP rtol=%G\n",snes->iter,kctx->version,rtol);CHKERRQ(ierr); 4426 PetscFunctionReturn(0); 4427 } 4428 4429 #undef __FUNCT__ 4430 #define __FUNCT__ "SNESKSPEW_PostSolve" 4431 static PetscErrorCode SNESKSPEW_PostSolve(SNES snes, KSP ksp, Vec b, Vec x) 4432 { 4433 PetscErrorCode ierr; 4434 SNESKSPEW *kctx = (SNESKSPEW*)snes->kspconvctx; 4435 PCSide pcside; 4436 Vec lres; 4437 4438 PetscFunctionBegin; 4439 if (!kctx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"No Eisenstat-Walker context exists"); 4440 ierr = KSPGetTolerances(ksp,&kctx->rtol_last,0,0,0);CHKERRQ(ierr); 4441 ierr = SNESGetFunctionNorm(snes,&kctx->norm_last);CHKERRQ(ierr); 4442 if (kctx->version == 1) { 4443 ierr = KSPGetPCSide(ksp,&pcside);CHKERRQ(ierr); 4444 if (pcside == PC_RIGHT) { /* XXX Should we also test KSP_UNPRECONDITIONED_NORM ? */ 4445 /* KSP residual is true linear residual */ 4446 ierr = KSPGetResidualNorm(ksp,&kctx->lresid_last);CHKERRQ(ierr); 4447 } else { 4448 /* KSP residual is preconditioned residual */ 4449 /* compute true linear residual norm */ 4450 ierr = VecDuplicate(b,&lres);CHKERRQ(ierr); 4451 ierr = MatMult(snes->jacobian,x,lres);CHKERRQ(ierr); 4452 ierr = VecAYPX(lres,-1.0,b);CHKERRQ(ierr); 4453 ierr = VecNorm(lres,NORM_2,&kctx->lresid_last);CHKERRQ(ierr); 4454 ierr = VecDestroy(&lres);CHKERRQ(ierr); 4455 } 4456 } 4457 PetscFunctionReturn(0); 4458 } 4459 4460 #undef __FUNCT__ 4461 #define __FUNCT__ "SNES_KSPSolve" 4462 PetscErrorCode SNES_KSPSolve(SNES snes, KSP ksp, Vec b, Vec x) 4463 { 4464 PetscErrorCode ierr; 4465 4466 PetscFunctionBegin; 4467 if (snes->ksp_ewconv) { ierr = SNESKSPEW_PreSolve(snes,ksp,b,x);CHKERRQ(ierr); } 4468 ierr = KSPSolve(ksp,b,x);CHKERRQ(ierr); 4469 if (snes->ksp_ewconv) { ierr = SNESKSPEW_PostSolve(snes,ksp,b,x);CHKERRQ(ierr); } 4470 PetscFunctionReturn(0); 4471 } 4472 4473 #include <petsc-private/dmimpl.h> 4474 #undef __FUNCT__ 4475 #define __FUNCT__ "SNESSetDM" 4476 /*@ 4477 SNESSetDM - Sets the DM that may be used by some preconditioners 4478 4479 Logically Collective on SNES 4480 4481 Input Parameters: 4482 + snes - the preconditioner context 4483 - dm - the dm 4484 4485 Level: intermediate 4486 4487 4488 .seealso: SNESGetDM(), KSPSetDM(), KSPGetDM() 4489 @*/ 4490 PetscErrorCode SNESSetDM(SNES snes,DM dm) 4491 { 4492 PetscErrorCode ierr; 4493 KSP ksp; 4494 DMSNES sdm; 4495 4496 PetscFunctionBegin; 4497 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4498 if (dm) {ierr = PetscObjectReference((PetscObject)dm);CHKERRQ(ierr);} 4499 if (snes->dm) { /* Move the DMSNES context over to the new DM unless the new DM already has one */ 4500 if (snes->dm->dmsnes && snes->dmAuto && !dm->dmsnes) { 4501 ierr = DMCopyDMSNES(snes->dm,dm);CHKERRQ(ierr); 4502 ierr = DMGetDMSNES(snes->dm,&sdm);CHKERRQ(ierr); 4503 if (sdm->originaldm == snes->dm) { /* Grant write privileges to the replacement DM */ 4504 sdm->originaldm = dm; 4505 } 4506 } 4507 ierr = DMDestroy(&snes->dm);CHKERRQ(ierr); 4508 } 4509 snes->dm = dm; 4510 snes->dmAuto = PETSC_FALSE; 4511 ierr = SNESGetKSP(snes,&ksp);CHKERRQ(ierr); 4512 ierr = KSPSetDM(ksp,dm);CHKERRQ(ierr); 4513 ierr = KSPSetDMActive(ksp,PETSC_FALSE);CHKERRQ(ierr); 4514 if (snes->pc) { 4515 ierr = SNESSetDM(snes->pc, snes->dm);CHKERRQ(ierr); 4516 ierr = SNESSetPCSide(snes,snes->pcside);CHKERRQ(ierr); 4517 } 4518 PetscFunctionReturn(0); 4519 } 4520 4521 #undef __FUNCT__ 4522 #define __FUNCT__ "SNESGetDM" 4523 /*@ 4524 SNESGetDM - Gets the DM that may be used by some preconditioners 4525 4526 Not Collective but DM obtained is parallel on SNES 4527 4528 Input Parameter: 4529 . snes - the preconditioner context 4530 4531 Output Parameter: 4532 . dm - the dm 4533 4534 Level: intermediate 4535 4536 4537 .seealso: SNESSetDM(), KSPSetDM(), KSPGetDM() 4538 @*/ 4539 PetscErrorCode SNESGetDM(SNES snes,DM *dm) 4540 { 4541 PetscErrorCode ierr; 4542 4543 PetscFunctionBegin; 4544 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4545 if (!snes->dm) { 4546 ierr = DMShellCreate(((PetscObject)snes)->comm,&snes->dm);CHKERRQ(ierr); 4547 snes->dmAuto = PETSC_TRUE; 4548 } 4549 *dm = snes->dm; 4550 PetscFunctionReturn(0); 4551 } 4552 4553 #undef __FUNCT__ 4554 #define __FUNCT__ "SNESSetPC" 4555 /*@ 4556 SNESSetPC - Sets the nonlinear preconditioner to be used. 4557 4558 Collective on SNES 4559 4560 Input Parameters: 4561 + snes - iterative context obtained from SNESCreate() 4562 - pc - the preconditioner object 4563 4564 Notes: 4565 Use SNESGetPC() to retrieve the preconditioner context (for example, 4566 to configure it using the API). 4567 4568 Level: developer 4569 4570 .keywords: SNES, set, precondition 4571 .seealso: SNESGetPC() 4572 @*/ 4573 PetscErrorCode SNESSetPC(SNES snes, SNES pc) 4574 { 4575 PetscErrorCode ierr; 4576 4577 PetscFunctionBegin; 4578 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); 4579 PetscValidHeaderSpecific(pc, SNES_CLASSID, 2); 4580 PetscCheckSameComm(snes, 1, pc, 2); 4581 ierr = PetscObjectReference((PetscObject) pc);CHKERRQ(ierr); 4582 ierr = SNESDestroy(&snes->pc);CHKERRQ(ierr); 4583 snes->pc = pc; 4584 ierr = PetscLogObjectParent(snes, snes->pc);CHKERRQ(ierr); 4585 PetscFunctionReturn(0); 4586 } 4587 4588 #undef __FUNCT__ 4589 #define __FUNCT__ "SNESGetPC" 4590 /*@ 4591 SNESGetPC - Returns a pointer to the nonlinear preconditioning context set with SNESSetPC(). 4592 4593 Not Collective 4594 4595 Input Parameter: 4596 . snes - iterative context obtained from SNESCreate() 4597 4598 Output Parameter: 4599 . pc - preconditioner context 4600 4601 Level: developer 4602 4603 .keywords: SNES, get, preconditioner 4604 .seealso: SNESSetPC() 4605 @*/ 4606 PetscErrorCode SNESGetPC(SNES snes, SNES *pc) 4607 { 4608 PetscErrorCode ierr; 4609 const char *optionsprefix; 4610 4611 PetscFunctionBegin; 4612 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); 4613 PetscValidPointer(pc, 2); 4614 if (!snes->pc) { 4615 ierr = SNESCreate(((PetscObject) snes)->comm,&snes->pc);CHKERRQ(ierr); 4616 ierr = PetscObjectIncrementTabLevel((PetscObject)snes->pc,(PetscObject)snes,1);CHKERRQ(ierr); 4617 ierr = PetscLogObjectParent(snes,snes->pc);CHKERRQ(ierr); 4618 ierr = SNESGetOptionsPrefix(snes,&optionsprefix);CHKERRQ(ierr); 4619 ierr = SNESSetOptionsPrefix(snes->pc,optionsprefix);CHKERRQ(ierr); 4620 ierr = SNESAppendOptionsPrefix(snes->pc,"npc_");CHKERRQ(ierr); 4621 } 4622 *pc = snes->pc; 4623 PetscFunctionReturn(0); 4624 } 4625 4626 4627 #undef __FUNCT__ 4628 #define __FUNCT__ "SNESSetPCSide" 4629 /*@ 4630 SNESSetPCSide - Sets the preconditioning side. 4631 4632 Logically Collective on SNES 4633 4634 Input Parameter: 4635 . snes - iterative context obtained from SNESCreate() 4636 4637 Output Parameter: 4638 . side - the preconditioning side, where side is one of 4639 .vb 4640 PC_LEFT - left preconditioning (default) 4641 PC_RIGHT - right preconditioning 4642 .ve 4643 4644 Options Database Keys: 4645 . -snes_pc_side <right,left> 4646 4647 Level: intermediate 4648 4649 .keywords: SNES, set, right, left, side, preconditioner, flag 4650 4651 .seealso: SNESGetPCSide(), KSPSetPCSide() 4652 @*/ 4653 PetscErrorCode SNESSetPCSide(SNES snes,PCSide side) 4654 { 4655 PetscFunctionBegin; 4656 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4657 PetscValidLogicalCollectiveEnum(snes,side,2); 4658 snes->pcside = side; 4659 PetscFunctionReturn(0); 4660 } 4661 4662 #undef __FUNCT__ 4663 #define __FUNCT__ "SNESGetPCSide" 4664 /*@ 4665 SNESGetPCSide - Gets the preconditioning side. 4666 4667 Not Collective 4668 4669 Input Parameter: 4670 . snes - iterative context obtained from SNESCreate() 4671 4672 Output Parameter: 4673 . side - the preconditioning side, where side is one of 4674 .vb 4675 PC_LEFT - left preconditioning (default) 4676 PC_RIGHT - right preconditioning 4677 .ve 4678 4679 Level: intermediate 4680 4681 .keywords: SNES, get, right, left, side, preconditioner, flag 4682 4683 .seealso: SNESSetPCSide(), KSPGetPCSide() 4684 @*/ 4685 PetscErrorCode SNESGetPCSide(SNES snes,PCSide *side) 4686 { 4687 PetscFunctionBegin; 4688 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4689 PetscValidPointer(side,2); 4690 *side = snes->pcside; 4691 PetscFunctionReturn(0); 4692 } 4693 4694 #undef __FUNCT__ 4695 #define __FUNCT__ "SNESSetSNESLineSearch" 4696 /*@ 4697 SNESSetSNESLineSearch - Sets the linesearch on the SNES instance. 4698 4699 Collective on SNES 4700 4701 Input Parameters: 4702 + snes - iterative context obtained from SNESCreate() 4703 - linesearch - the linesearch object 4704 4705 Notes: 4706 Use SNESGetSNESLineSearch() to retrieve the preconditioner context (for example, 4707 to configure it using the API). 4708 4709 Level: developer 4710 4711 .keywords: SNES, set, linesearch 4712 .seealso: SNESGetSNESLineSearch() 4713 @*/ 4714 PetscErrorCode SNESSetSNESLineSearch(SNES snes, SNESLineSearch linesearch) 4715 { 4716 PetscErrorCode ierr; 4717 4718 PetscFunctionBegin; 4719 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); 4720 PetscValidHeaderSpecific(linesearch, SNESLINESEARCH_CLASSID, 2); 4721 PetscCheckSameComm(snes, 1, linesearch, 2); 4722 ierr = PetscObjectReference((PetscObject) linesearch);CHKERRQ(ierr); 4723 ierr = SNESLineSearchDestroy(&snes->linesearch);CHKERRQ(ierr); 4724 snes->linesearch = linesearch; 4725 ierr = PetscLogObjectParent(snes, snes->linesearch);CHKERRQ(ierr); 4726 PetscFunctionReturn(0); 4727 } 4728 4729 #undef __FUNCT__ 4730 #define __FUNCT__ "SNESGetSNESLineSearch" 4731 /*@C 4732 SNESGetSNESLineSearch - Returns a pointer to the line search context set with SNESSetLineSearch() 4733 or creates a default line search instance associated with the SNES and returns it. 4734 4735 Not Collective 4736 4737 Input Parameter: 4738 . snes - iterative context obtained from SNESCreate() 4739 4740 Output Parameter: 4741 . linesearch - linesearch context 4742 4743 Level: developer 4744 4745 .keywords: SNES, get, linesearch 4746 .seealso: SNESSetSNESLineSearch() 4747 @*/ 4748 PetscErrorCode SNESGetSNESLineSearch(SNES snes, SNESLineSearch *linesearch) 4749 { 4750 PetscErrorCode ierr; 4751 const char *optionsprefix; 4752 4753 PetscFunctionBegin; 4754 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); 4755 PetscValidPointer(linesearch, 2); 4756 if (!snes->linesearch) { 4757 ierr = SNESGetOptionsPrefix(snes, &optionsprefix);CHKERRQ(ierr); 4758 ierr = SNESLineSearchCreate(((PetscObject) snes)->comm, &snes->linesearch);CHKERRQ(ierr); 4759 ierr = SNESLineSearchSetSNES(snes->linesearch, snes);CHKERRQ(ierr); 4760 ierr = SNESLineSearchAppendOptionsPrefix(snes->linesearch, optionsprefix);CHKERRQ(ierr); 4761 ierr = PetscObjectIncrementTabLevel((PetscObject) snes->linesearch, (PetscObject) snes, 1);CHKERRQ(ierr); 4762 ierr = PetscLogObjectParent(snes, snes->linesearch);CHKERRQ(ierr); 4763 } 4764 *linesearch = snes->linesearch; 4765 PetscFunctionReturn(0); 4766 } 4767 4768 #if defined(PETSC_HAVE_MATLAB_ENGINE) 4769 #include <mex.h> 4770 4771 typedef struct {char *funcname; mxArray *ctx;} SNESMatlabContext; 4772 4773 #undef __FUNCT__ 4774 #define __FUNCT__ "SNESComputeFunction_Matlab" 4775 /* 4776 SNESComputeFunction_Matlab - Calls the function that has been set with 4777 SNESSetFunctionMatlab(). 4778 4779 Collective on SNES 4780 4781 Input Parameters: 4782 + snes - the SNES context 4783 - x - input vector 4784 4785 Output Parameter: 4786 . y - function vector, as set by SNESSetFunction() 4787 4788 Notes: 4789 SNESComputeFunction() is typically used within nonlinear solvers 4790 implementations, so most users would not generally call this routine 4791 themselves. 4792 4793 Level: developer 4794 4795 .keywords: SNES, nonlinear, compute, function 4796 4797 .seealso: SNESSetFunction(), SNESGetFunction() 4798 */ 4799 PetscErrorCode SNESComputeFunction_Matlab(SNES snes,Vec x,Vec y, void *ctx) 4800 { 4801 PetscErrorCode ierr; 4802 SNESMatlabContext *sctx = (SNESMatlabContext *)ctx; 4803 int nlhs = 1,nrhs = 5; 4804 mxArray *plhs[1],*prhs[5]; 4805 long long int lx = 0,ly = 0,ls = 0; 4806 4807 PetscFunctionBegin; 4808 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4809 PetscValidHeaderSpecific(x,VEC_CLASSID,2); 4810 PetscValidHeaderSpecific(y,VEC_CLASSID,3); 4811 PetscCheckSameComm(snes,1,x,2); 4812 PetscCheckSameComm(snes,1,y,3); 4813 4814 /* call Matlab function in ctx with arguments x and y */ 4815 4816 ierr = PetscMemcpy(&ls,&snes,sizeof(snes));CHKERRQ(ierr); 4817 ierr = PetscMemcpy(&lx,&x,sizeof(x));CHKERRQ(ierr); 4818 ierr = PetscMemcpy(&ly,&y,sizeof(x));CHKERRQ(ierr); 4819 prhs[0] = mxCreateDoubleScalar((double)ls); 4820 prhs[1] = mxCreateDoubleScalar((double)lx); 4821 prhs[2] = mxCreateDoubleScalar((double)ly); 4822 prhs[3] = mxCreateString(sctx->funcname); 4823 prhs[4] = sctx->ctx; 4824 ierr = mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscSNESComputeFunctionInternal");CHKERRQ(ierr); 4825 ierr = mxGetScalar(plhs[0]);CHKERRQ(ierr); 4826 mxDestroyArray(prhs[0]); 4827 mxDestroyArray(prhs[1]); 4828 mxDestroyArray(prhs[2]); 4829 mxDestroyArray(prhs[3]); 4830 mxDestroyArray(plhs[0]); 4831 PetscFunctionReturn(0); 4832 } 4833 4834 4835 #undef __FUNCT__ 4836 #define __FUNCT__ "SNESSetFunctionMatlab" 4837 /* 4838 SNESSetFunctionMatlab - Sets the function evaluation routine and function 4839 vector for use by the SNES routines in solving systems of nonlinear 4840 equations from MATLAB. Here the function is a string containing the name of a MATLAB function 4841 4842 Logically Collective on SNES 4843 4844 Input Parameters: 4845 + snes - the SNES context 4846 . r - vector to store function value 4847 - func - function evaluation routine 4848 4849 Calling sequence of func: 4850 $ func (SNES snes,Vec x,Vec f,void *ctx); 4851 4852 4853 Notes: 4854 The Newton-like methods typically solve linear systems of the form 4855 $ f'(x) x = -f(x), 4856 where f'(x) denotes the Jacobian matrix and f(x) is the function. 4857 4858 Level: beginner 4859 4860 .keywords: SNES, nonlinear, set, function 4861 4862 .seealso: SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESSetFunction() 4863 */ 4864 PetscErrorCode SNESSetFunctionMatlab(SNES snes,Vec r,const char *func,mxArray *ctx) 4865 { 4866 PetscErrorCode ierr; 4867 SNESMatlabContext *sctx; 4868 4869 PetscFunctionBegin; 4870 /* currently sctx is memory bleed */ 4871 ierr = PetscMalloc(sizeof(SNESMatlabContext),&sctx);CHKERRQ(ierr); 4872 ierr = PetscStrallocpy(func,&sctx->funcname);CHKERRQ(ierr); 4873 /* 4874 This should work, but it doesn't 4875 sctx->ctx = ctx; 4876 mexMakeArrayPersistent(sctx->ctx); 4877 */ 4878 sctx->ctx = mxDuplicateArray(ctx); 4879 ierr = SNESSetFunction(snes,r,SNESComputeFunction_Matlab,sctx);CHKERRQ(ierr); 4880 PetscFunctionReturn(0); 4881 } 4882 4883 #undef __FUNCT__ 4884 #define __FUNCT__ "SNESComputeJacobian_Matlab" 4885 /* 4886 SNESComputeJacobian_Matlab - Calls the function that has been set with 4887 SNESSetJacobianMatlab(). 4888 4889 Collective on SNES 4890 4891 Input Parameters: 4892 + snes - the SNES context 4893 . x - input vector 4894 . A, B - the matrices 4895 - ctx - user context 4896 4897 Output Parameter: 4898 . flag - structure of the matrix 4899 4900 Level: developer 4901 4902 .keywords: SNES, nonlinear, compute, function 4903 4904 .seealso: SNESSetFunction(), SNESGetFunction() 4905 @*/ 4906 PetscErrorCode SNESComputeJacobian_Matlab(SNES snes,Vec x,Mat *A,Mat *B,MatStructure *flag, void *ctx) 4907 { 4908 PetscErrorCode ierr; 4909 SNESMatlabContext *sctx = (SNESMatlabContext *)ctx; 4910 int nlhs = 2,nrhs = 6; 4911 mxArray *plhs[2],*prhs[6]; 4912 long long int lx = 0,lA = 0,ls = 0, lB = 0; 4913 4914 PetscFunctionBegin; 4915 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4916 PetscValidHeaderSpecific(x,VEC_CLASSID,2); 4917 4918 /* call Matlab function in ctx with arguments x and y */ 4919 4920 ierr = PetscMemcpy(&ls,&snes,sizeof(snes));CHKERRQ(ierr); 4921 ierr = PetscMemcpy(&lx,&x,sizeof(x));CHKERRQ(ierr); 4922 ierr = PetscMemcpy(&lA,A,sizeof(x));CHKERRQ(ierr); 4923 ierr = PetscMemcpy(&lB,B,sizeof(x));CHKERRQ(ierr); 4924 prhs[0] = mxCreateDoubleScalar((double)ls); 4925 prhs[1] = mxCreateDoubleScalar((double)lx); 4926 prhs[2] = mxCreateDoubleScalar((double)lA); 4927 prhs[3] = mxCreateDoubleScalar((double)lB); 4928 prhs[4] = mxCreateString(sctx->funcname); 4929 prhs[5] = sctx->ctx; 4930 ierr = mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscSNESComputeJacobianInternal");CHKERRQ(ierr); 4931 ierr = mxGetScalar(plhs[0]);CHKERRQ(ierr); 4932 *flag = (MatStructure) mxGetScalar(plhs[1]);CHKERRQ(ierr); 4933 mxDestroyArray(prhs[0]); 4934 mxDestroyArray(prhs[1]); 4935 mxDestroyArray(prhs[2]); 4936 mxDestroyArray(prhs[3]); 4937 mxDestroyArray(prhs[4]); 4938 mxDestroyArray(plhs[0]); 4939 mxDestroyArray(plhs[1]); 4940 PetscFunctionReturn(0); 4941 } 4942 4943 4944 #undef __FUNCT__ 4945 #define __FUNCT__ "SNESSetJacobianMatlab" 4946 /* 4947 SNESSetJacobianMatlab - Sets the Jacobian function evaluation routine and two empty Jacobian matrices 4948 vector for use by the SNES routines in solving systems of nonlinear 4949 equations from MATLAB. Here the function is a string containing the name of a MATLAB function 4950 4951 Logically Collective on SNES 4952 4953 Input Parameters: 4954 + snes - the SNES context 4955 . A,B - Jacobian matrices 4956 . func - function evaluation routine 4957 - ctx - user context 4958 4959 Calling sequence of func: 4960 $ flag = func (SNES snes,Vec x,Mat A,Mat B,void *ctx); 4961 4962 4963 Level: developer 4964 4965 .keywords: SNES, nonlinear, set, function 4966 4967 .seealso: SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESSetFunction() 4968 */ 4969 PetscErrorCode SNESSetJacobianMatlab(SNES snes,Mat A,Mat B,const char *func,mxArray *ctx) 4970 { 4971 PetscErrorCode ierr; 4972 SNESMatlabContext *sctx; 4973 4974 PetscFunctionBegin; 4975 /* currently sctx is memory bleed */ 4976 ierr = PetscMalloc(sizeof(SNESMatlabContext),&sctx);CHKERRQ(ierr); 4977 ierr = PetscStrallocpy(func,&sctx->funcname);CHKERRQ(ierr); 4978 /* 4979 This should work, but it doesn't 4980 sctx->ctx = ctx; 4981 mexMakeArrayPersistent(sctx->ctx); 4982 */ 4983 sctx->ctx = mxDuplicateArray(ctx); 4984 ierr = SNESSetJacobian(snes,A,B,SNESComputeJacobian_Matlab,sctx);CHKERRQ(ierr); 4985 PetscFunctionReturn(0); 4986 } 4987 4988 #undef __FUNCT__ 4989 #define __FUNCT__ "SNESMonitor_Matlab" 4990 /* 4991 SNESMonitor_Matlab - Calls the function that has been set with SNESMonitorSetMatlab(). 4992 4993 Collective on SNES 4994 4995 .seealso: SNESSetFunction(), SNESGetFunction() 4996 @*/ 4997 PetscErrorCode SNESMonitor_Matlab(SNES snes,PetscInt it, PetscReal fnorm, void *ctx) 4998 { 4999 PetscErrorCode ierr; 5000 SNESMatlabContext *sctx = (SNESMatlabContext *)ctx; 5001 int nlhs = 1,nrhs = 6; 5002 mxArray *plhs[1],*prhs[6]; 5003 long long int lx = 0,ls = 0; 5004 Vec x=snes->vec_sol; 5005 5006 PetscFunctionBegin; 5007 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 5008 5009 ierr = PetscMemcpy(&ls,&snes,sizeof(snes));CHKERRQ(ierr); 5010 ierr = PetscMemcpy(&lx,&x,sizeof(x));CHKERRQ(ierr); 5011 prhs[0] = mxCreateDoubleScalar((double)ls); 5012 prhs[1] = mxCreateDoubleScalar((double)it); 5013 prhs[2] = mxCreateDoubleScalar((double)fnorm); 5014 prhs[3] = mxCreateDoubleScalar((double)lx); 5015 prhs[4] = mxCreateString(sctx->funcname); 5016 prhs[5] = sctx->ctx; 5017 ierr = mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscSNESMonitorInternal");CHKERRQ(ierr); 5018 ierr = mxGetScalar(plhs[0]);CHKERRQ(ierr); 5019 mxDestroyArray(prhs[0]); 5020 mxDestroyArray(prhs[1]); 5021 mxDestroyArray(prhs[2]); 5022 mxDestroyArray(prhs[3]); 5023 mxDestroyArray(prhs[4]); 5024 mxDestroyArray(plhs[0]); 5025 PetscFunctionReturn(0); 5026 } 5027 5028 5029 #undef __FUNCT__ 5030 #define __FUNCT__ "SNESMonitorSetMatlab" 5031 /* 5032 SNESMonitorSetMatlab - Sets the monitor function from MATLAB 5033 5034 Level: developer 5035 5036 .keywords: SNES, nonlinear, set, function 5037 5038 .seealso: SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESSetFunction() 5039 */ 5040 PetscErrorCode SNESMonitorSetMatlab(SNES snes,const char *func,mxArray *ctx) 5041 { 5042 PetscErrorCode ierr; 5043 SNESMatlabContext *sctx; 5044 5045 PetscFunctionBegin; 5046 /* currently sctx is memory bleed */ 5047 ierr = PetscMalloc(sizeof(SNESMatlabContext),&sctx);CHKERRQ(ierr); 5048 ierr = PetscStrallocpy(func,&sctx->funcname);CHKERRQ(ierr); 5049 /* 5050 This should work, but it doesn't 5051 sctx->ctx = ctx; 5052 mexMakeArrayPersistent(sctx->ctx); 5053 */ 5054 sctx->ctx = mxDuplicateArray(ctx); 5055 ierr = SNESMonitorSet(snes,SNESMonitor_Matlab,sctx,PETSC_NULL);CHKERRQ(ierr); 5056 PetscFunctionReturn(0); 5057 } 5058 5059 #endif 5060