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