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