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