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