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