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