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