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