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 PetscFunctionReturn(0); 1789 } 1790 1791 #undef __FUNCT__ 1792 #define __FUNCT__ "SNESPicardComputeJacobian" 1793 PetscErrorCode SNESPicardComputeJacobian(SNES snes,Vec x1,Mat *J,Mat *B,MatStructure *flag,void *ctx) 1794 { 1795 PetscFunctionBegin; 1796 /* the jacobian matrix should be pre-filled in SNESPicardComputeFunction */ 1797 *flag = snes->matstruct; 1798 PetscFunctionReturn(0); 1799 } 1800 1801 #undef __FUNCT__ 1802 #define __FUNCT__ "SNESSetPicard" 1803 /*@C 1804 SNESSetPicard - Use SNES to solve the semilinear-system A(x) x = b(x) via a Picard type iteration (Picard linearization) 1805 1806 Logically Collective on SNES 1807 1808 Input Parameters: 1809 + snes - the SNES context 1810 . r - vector to store function value 1811 . func - function evaluation routine 1812 . 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) 1813 . mat - matrix to store A 1814 . mfunc - function to compute matrix value 1815 - ctx - [optional] user-defined context for private data for the 1816 function evaluation routine (may be PETSC_NULL) 1817 1818 Calling sequence of func: 1819 $ func (SNES snes,Vec x,Vec f,void *ctx); 1820 1821 + f - function vector 1822 - ctx - optional user-defined function context 1823 1824 Calling sequence of mfunc: 1825 $ mfunc (SNES snes,Vec x,Mat *jmat,Mat *mat,int *flag,void *ctx); 1826 1827 + x - input vector 1828 . 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(), 1829 normally just pass mat in this location 1830 . mat - form A(x) matrix 1831 . flag - flag indicating information about the preconditioner matrix 1832 structure (same as flag in KSPSetOperators()), one of SAME_NONZERO_PATTERN,DIFFERENT_NONZERO_PATTERN,SAME_PRECONDITIONER 1833 - ctx - [optional] user-defined Jacobian context 1834 1835 Notes: 1836 One can call SNESSetPicard() or SNESSetFunction() (and possibly SNESSetJacobian()) but cannot call both 1837 1838 $ 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} 1839 $ Note that when an exact solver is used this corresponds to the "classic" Picard A(x^{n}) x^{n+1} = b(x^{n}) iteration. 1840 1841 Run with -snes_mf_operator to solve the system with Newton's method using A(x^{n}) to construct the preconditioner. 1842 1843 We implement the defect correction form of the Picard iteration because it converges much more generally when inexact linear solvers are used then 1844 the direct Picard iteration A(x^n) x^{n+1} = b(x^n) 1845 1846 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 1847 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 1848 different please contact us at petsc-dev@mcs.anl.gov and we'll have an entirely new argument :-). 1849 1850 Level: beginner 1851 1852 .keywords: SNES, nonlinear, set, function 1853 1854 .seealso: SNESGetFunction(), SNESSetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESGetPicard(), SNESLineSearchPreCheckPicard() 1855 @*/ 1856 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) 1857 { 1858 PetscErrorCode ierr; 1859 DM dm; 1860 1861 PetscFunctionBegin; 1862 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1863 ierr = SNESGetDM(snes, &dm);CHKERRQ(ierr); 1864 ierr = DMSNESSetPicard(dm,func,mfunc,ctx);CHKERRQ(ierr); 1865 ierr = SNESSetFunction(snes,r,SNESPicardComputeFunction,ctx);CHKERRQ(ierr); 1866 ierr = SNESSetJacobian(snes,jmat,mat,SNESPicardComputeJacobian,ctx);CHKERRQ(ierr); 1867 PetscFunctionReturn(0); 1868 } 1869 1870 1871 #undef __FUNCT__ 1872 #define __FUNCT__ "SNESGetPicard" 1873 /*@C 1874 SNESGetPicard - Returns the context for the Picard iteration 1875 1876 Not Collective, but Vec is parallel if SNES is parallel. Collective if Vec is requested, but has not been created yet. 1877 1878 Input Parameter: 1879 . snes - the SNES context 1880 1881 Output Parameter: 1882 + r - the function (or PETSC_NULL) 1883 . func - the function (or PETSC_NULL) 1884 . jmat - the picard matrix (or PETSC_NULL) 1885 . mat - the picard preconditioner matrix (or PETSC_NULL) 1886 . mfunc - the function for matrix evaluation (or PETSC_NULL) 1887 - ctx - the function context (or PETSC_NULL) 1888 1889 Level: advanced 1890 1891 .keywords: SNES, nonlinear, get, function 1892 1893 .seealso: SNESSetPicard, SNESGetFunction, SNESGetJacobian, SNESGetDM 1894 @*/ 1895 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) 1896 { 1897 PetscErrorCode ierr; 1898 DM dm; 1899 1900 PetscFunctionBegin; 1901 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1902 ierr = SNESGetFunction(snes,r,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr); 1903 ierr = SNESGetJacobian(snes,jmat,mat,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr); 1904 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 1905 ierr = DMSNESGetPicard(dm,func,mfunc,ctx);CHKERRQ(ierr); 1906 PetscFunctionReturn(0); 1907 } 1908 1909 #undef __FUNCT__ 1910 #define __FUNCT__ "SNESSetComputeInitialGuess" 1911 /*@C 1912 SNESSetComputeInitialGuess - Sets a routine used to compute an initial guess for the problem 1913 1914 Logically Collective on SNES 1915 1916 Input Parameters: 1917 + snes - the SNES context 1918 . func - function evaluation routine 1919 - ctx - [optional] user-defined context for private data for the 1920 function evaluation routine (may be PETSC_NULL) 1921 1922 Calling sequence of func: 1923 $ func (SNES snes,Vec x,void *ctx); 1924 1925 . f - function vector 1926 - ctx - optional user-defined function context 1927 1928 Level: intermediate 1929 1930 .keywords: SNES, nonlinear, set, function 1931 1932 .seealso: SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian() 1933 @*/ 1934 PetscErrorCode SNESSetComputeInitialGuess(SNES snes,PetscErrorCode (*func)(SNES,Vec,void*),void *ctx) 1935 { 1936 PetscFunctionBegin; 1937 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1938 if (func) snes->ops->computeinitialguess = func; 1939 if (ctx) snes->initialguessP = ctx; 1940 PetscFunctionReturn(0); 1941 } 1942 1943 /* --------------------------------------------------------------- */ 1944 #undef __FUNCT__ 1945 #define __FUNCT__ "SNESGetRhs" 1946 /*@C 1947 SNESGetRhs - Gets the vector for solving F(x) = rhs. If rhs is not set 1948 it assumes a zero right hand side. 1949 1950 Logically Collective on SNES 1951 1952 Input Parameter: 1953 . snes - the SNES context 1954 1955 Output Parameter: 1956 . rhs - the right hand side vector or PETSC_NULL if the right hand side vector is null 1957 1958 Level: intermediate 1959 1960 .keywords: SNES, nonlinear, get, function, right hand side 1961 1962 .seealso: SNESGetSolution(), SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESSetFunction() 1963 @*/ 1964 PetscErrorCode SNESGetRhs(SNES snes,Vec *rhs) 1965 { 1966 PetscFunctionBegin; 1967 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1968 PetscValidPointer(rhs,2); 1969 *rhs = snes->vec_rhs; 1970 PetscFunctionReturn(0); 1971 } 1972 1973 #undef __FUNCT__ 1974 #define __FUNCT__ "SNESComputeFunction" 1975 /*@ 1976 SNESComputeFunction - Calls the function that has been set with 1977 SNESSetFunction(). 1978 1979 Collective on SNES 1980 1981 Input Parameters: 1982 + snes - the SNES context 1983 - x - input vector 1984 1985 Output Parameter: 1986 . y - function vector, as set by SNESSetFunction() 1987 1988 Notes: 1989 SNESComputeFunction() is typically used within nonlinear solvers 1990 implementations, so most users would not generally call this routine 1991 themselves. 1992 1993 Level: developer 1994 1995 .keywords: SNES, nonlinear, compute, function 1996 1997 .seealso: SNESSetFunction(), SNESGetFunction() 1998 @*/ 1999 PetscErrorCode SNESComputeFunction(SNES snes,Vec x,Vec y) 2000 { 2001 PetscErrorCode ierr; 2002 DM dm; 2003 DMSNES sdm; 2004 2005 PetscFunctionBegin; 2006 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2007 PetscValidHeaderSpecific(x,VEC_CLASSID,2); 2008 PetscValidHeaderSpecific(y,VEC_CLASSID,3); 2009 PetscCheckSameComm(snes,1,x,2); 2010 PetscCheckSameComm(snes,1,y,3); 2011 VecValidValues(x,2,PETSC_TRUE); 2012 2013 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 2014 ierr = DMGetDMSNES(dm,&sdm);CHKERRQ(ierr); 2015 ierr = PetscLogEventBegin(SNES_FunctionEval,snes,x,y,0);CHKERRQ(ierr); 2016 if (snes->pc && snes->pcside == PC_LEFT) { 2017 ierr = VecCopy(x,y);CHKERRQ(ierr); 2018 ierr = SNESSolve(snes->pc,snes->vec_rhs,y);CHKERRQ(ierr); 2019 ierr = VecAYPX(y,-1.0,x);CHKERRQ(ierr); 2020 } else if (sdm->ops->computefunction) { 2021 PetscStackPush("SNES user function"); 2022 ierr = (*sdm->ops->computefunction)(snes,x,y,sdm->functionctx);CHKERRQ(ierr); 2023 PetscStackPop; 2024 } else if (snes->vec_rhs) { 2025 ierr = MatMult(snes->jacobian, x, y);CHKERRQ(ierr); 2026 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE, "Must call SNESSetFunction() or SNESSetDM() before SNESComputeFunction(), likely called from SNESSolve()."); 2027 if (snes->vec_rhs) { 2028 ierr = VecAXPY(y,-1.0,snes->vec_rhs);CHKERRQ(ierr); 2029 } 2030 snes->nfuncs++; 2031 ierr = PetscLogEventEnd(SNES_FunctionEval,snes,x,y,0);CHKERRQ(ierr); 2032 VecValidValues(y,3,PETSC_FALSE); 2033 PetscFunctionReturn(0); 2034 } 2035 2036 #undef __FUNCT__ 2037 #define __FUNCT__ "SNESComputeGS" 2038 /*@ 2039 SNESComputeGS - Calls the Gauss-Seidel function that has been set with 2040 SNESSetGS(). 2041 2042 Collective on SNES 2043 2044 Input Parameters: 2045 + snes - the SNES context 2046 . x - input vector 2047 - b - rhs vector 2048 2049 Output Parameter: 2050 . x - new solution vector 2051 2052 Notes: 2053 SNESComputeGS() is typically used within composed nonlinear solver 2054 implementations, so most users would not generally call this routine 2055 themselves. 2056 2057 Level: developer 2058 2059 .keywords: SNES, nonlinear, compute, function 2060 2061 .seealso: SNESSetGS(), SNESComputeFunction() 2062 @*/ 2063 PetscErrorCode SNESComputeGS(SNES snes,Vec b,Vec x) 2064 { 2065 PetscErrorCode ierr; 2066 PetscInt i; 2067 DM dm; 2068 DMSNES sdm; 2069 2070 PetscFunctionBegin; 2071 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2072 PetscValidHeaderSpecific(x,VEC_CLASSID,2); 2073 if (b) PetscValidHeaderSpecific(b,VEC_CLASSID,3); 2074 PetscCheckSameComm(snes,1,x,2); 2075 if (b) PetscCheckSameComm(snes,1,b,3); 2076 if (b) VecValidValues(b,2,PETSC_TRUE); 2077 ierr = PetscLogEventBegin(SNES_GSEval,snes,x,b,0);CHKERRQ(ierr); 2078 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 2079 ierr = DMGetDMSNES(dm,&sdm);CHKERRQ(ierr); 2080 if (sdm->ops->computegs) { 2081 for (i = 0; i < snes->gssweeps; i++) { 2082 PetscStackPush("SNES user GS"); 2083 ierr = (*sdm->ops->computegs)(snes,x,b,sdm->gsctx);CHKERRQ(ierr); 2084 PetscStackPop; 2085 } 2086 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE, "Must call SNESSetGS() before SNESComputeGS(), likely called from SNESSolve()."); 2087 ierr = PetscLogEventEnd(SNES_GSEval,snes,x,b,0);CHKERRQ(ierr); 2088 VecValidValues(x,3,PETSC_FALSE); 2089 PetscFunctionReturn(0); 2090 } 2091 2092 2093 #undef __FUNCT__ 2094 #define __FUNCT__ "SNESComputeJacobian" 2095 /*@ 2096 SNESComputeJacobian - Computes the Jacobian matrix that has been 2097 set with SNESSetJacobian(). 2098 2099 Collective on SNES and Mat 2100 2101 Input Parameters: 2102 + snes - the SNES context 2103 - x - input vector 2104 2105 Output Parameters: 2106 + A - Jacobian matrix 2107 . B - optional preconditioning matrix 2108 - flag - flag indicating matrix structure (one of, SAME_NONZERO_PATTERN,DIFFERENT_NONZERO_PATTERN,SAME_PRECONDITIONER) 2109 2110 Options Database Keys: 2111 + -snes_lag_preconditioner <lag> 2112 . -snes_lag_jacobian <lag> 2113 . -snes_compare_explicit - Compare the computed Jacobian to the finite difference Jacobian and output the differences 2114 . -snes_compare_explicit_draw - Compare the computed Jacobian to the finite difference Jacobian and draw the result 2115 . -snes_compare_explicit_contour - Compare the computed Jacobian to the finite difference Jacobian and draw a contour plot with the result 2116 . -snes_compare_operator - Make the comparison options above use the operator instead of the preconditioning matrix 2117 . -snes_compare_coloring - Compute the finite differece Jacobian using coloring and display norms of difference 2118 . -snes_compare_coloring_display - Compute the finite differece Jacobian using coloring and display verbose differences 2119 . -snes_compare_coloring_threshold - Display only those matrix entries that differ by more than a given threshold 2120 . -snes_compare_coloring_threshold_atol - Absolute tolerance for difference in matrix entries to be displayed by -snes_compare_coloring_threshold 2121 . -snes_compare_coloring_threshold_rtol - Relative tolerance for difference in matrix entries to be displayed by -snes_compare_coloring_threshold 2122 . -snes_compare_coloring_draw - Compute the finite differece Jacobian using coloring and draw differences 2123 - -snes_compare_coloring_draw_contour - Compute the finite differece Jacobian using coloring and show contours of matrices and differences 2124 2125 2126 Notes: 2127 Most users should not need to explicitly call this routine, as it 2128 is used internally within the nonlinear solvers. 2129 2130 See KSPSetOperators() for important information about setting the 2131 flag parameter. 2132 2133 Level: developer 2134 2135 .keywords: SNES, compute, Jacobian, matrix 2136 2137 .seealso: SNESSetJacobian(), KSPSetOperators(), MatStructure, SNESSetLagPreconditioner(), SNESSetLagJacobian() 2138 @*/ 2139 PetscErrorCode SNESComputeJacobian(SNES snes,Vec X,Mat *A,Mat *B,MatStructure *flg) 2140 { 2141 PetscErrorCode ierr; 2142 PetscBool flag; 2143 DM dm; 2144 DMSNES sdm; 2145 2146 PetscFunctionBegin; 2147 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2148 PetscValidHeaderSpecific(X,VEC_CLASSID,2); 2149 PetscValidPointer(flg,5); 2150 PetscCheckSameComm(snes,1,X,2); 2151 VecValidValues(X,2,PETSC_TRUE); 2152 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 2153 ierr = DMGetDMSNES(dm,&sdm);CHKERRQ(ierr); 2154 2155 if (!sdm->ops->computejacobian) SETERRQ(((PetscObject)snes)->comm,PETSC_ERR_USER,"Must call SNESSetJacobian(), DMSNESSetJacobian(), DMDASNESSetJacobianLocal(), etc"); 2156 2157 /* make sure that MatAssemblyBegin/End() is called on A matrix if it is matrix free */ 2158 2159 if (snes->lagjacobian == -2) { 2160 snes->lagjacobian = -1; 2161 ierr = PetscInfo(snes,"Recomputing Jacobian/preconditioner because lag is -2 (means compute Jacobian, but then never again) \n");CHKERRQ(ierr); 2162 } else if (snes->lagjacobian == -1) { 2163 *flg = SAME_PRECONDITIONER; 2164 ierr = PetscInfo(snes,"Reusing Jacobian/preconditioner because lag is -1\n");CHKERRQ(ierr); 2165 ierr = PetscObjectTypeCompare((PetscObject)*A,MATMFFD,&flag);CHKERRQ(ierr); 2166 if (flag) { 2167 ierr = MatAssemblyBegin(*A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2168 ierr = MatAssemblyEnd(*A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2169 } 2170 PetscFunctionReturn(0); 2171 } else if (snes->lagjacobian > 1 && snes->iter % snes->lagjacobian) { 2172 *flg = SAME_PRECONDITIONER; 2173 ierr = PetscInfo2(snes,"Reusing Jacobian/preconditioner because lag is %D and SNES iteration is %D\n",snes->lagjacobian,snes->iter);CHKERRQ(ierr); 2174 ierr = PetscObjectTypeCompare((PetscObject)*A,MATMFFD,&flag);CHKERRQ(ierr); 2175 if (flag) { 2176 ierr = MatAssemblyBegin(*A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2177 ierr = MatAssemblyEnd(*A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2178 } 2179 PetscFunctionReturn(0); 2180 } 2181 2182 *flg = DIFFERENT_NONZERO_PATTERN; 2183 ierr = PetscLogEventBegin(SNES_JacobianEval,snes,X,*A,*B);CHKERRQ(ierr); 2184 2185 PetscStackPush("SNES user Jacobian function"); 2186 ierr = (*sdm->ops->computejacobian)(snes,X,A,B,flg,sdm->jacobianctx);CHKERRQ(ierr); 2187 PetscStackPop; 2188 2189 ierr = PetscLogEventEnd(SNES_JacobianEval,snes,X,*A,*B);CHKERRQ(ierr); 2190 2191 if (snes->lagpreconditioner == -2) { 2192 ierr = PetscInfo(snes,"Rebuilding preconditioner exactly once since lag is -2\n");CHKERRQ(ierr); 2193 snes->lagpreconditioner = -1; 2194 } else if (snes->lagpreconditioner == -1) { 2195 *flg = SAME_PRECONDITIONER; 2196 ierr = PetscInfo(snes,"Reusing preconditioner because lag is -1\n");CHKERRQ(ierr); 2197 } else if (snes->lagpreconditioner > 1 && snes->iter % snes->lagpreconditioner) { 2198 *flg = SAME_PRECONDITIONER; 2199 ierr = PetscInfo2(snes,"Reusing preconditioner because lag is %D and SNES iteration is %D\n",snes->lagpreconditioner,snes->iter);CHKERRQ(ierr); 2200 } 2201 2202 /* make sure user returned a correct Jacobian and preconditioner */ 2203 /* PetscValidHeaderSpecific(*A,MAT_CLASSID,3); 2204 PetscValidHeaderSpecific(*B,MAT_CLASSID,4); */ 2205 { 2206 PetscBool flag = PETSC_FALSE,flag_draw = PETSC_FALSE,flag_contour = PETSC_FALSE,flag_operator = PETSC_FALSE; 2207 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_explicit",&flag,PETSC_NULL);CHKERRQ(ierr); 2208 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_explicit_draw",&flag_draw,PETSC_NULL);CHKERRQ(ierr); 2209 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_explicit_draw_contour",&flag_contour,PETSC_NULL);CHKERRQ(ierr); 2210 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_operator",&flag_operator,PETSC_NULL);CHKERRQ(ierr); 2211 if (flag || flag_draw || flag_contour) { 2212 Mat Bexp_mine = PETSC_NULL,Bexp,FDexp; 2213 MatStructure mstruct; 2214 PetscViewer vdraw,vstdout; 2215 PetscBool flg; 2216 if (flag_operator) { 2217 ierr = MatComputeExplicitOperator(*A,&Bexp_mine);CHKERRQ(ierr); 2218 Bexp = Bexp_mine; 2219 } else { 2220 /* See if the preconditioning matrix can be viewed and added directly */ 2221 ierr = PetscObjectTypeCompareAny((PetscObject)*B,&flg,MATSEQAIJ,MATMPIAIJ,MATSEQDENSE,MATMPIDENSE,MATSEQBAIJ,MATMPIBAIJ,MATSEQSBAIJ,MATMPIBAIJ,"");CHKERRQ(ierr); 2222 if (flg) Bexp = *B; 2223 else { 2224 /* If the "preconditioning" matrix is itself MATSHELL or some other type without direct support */ 2225 ierr = MatComputeExplicitOperator(*B,&Bexp_mine);CHKERRQ(ierr); 2226 Bexp = Bexp_mine; 2227 } 2228 } 2229 ierr = MatConvert(Bexp,MATSAME,MAT_INITIAL_MATRIX,&FDexp);CHKERRQ(ierr); 2230 ierr = SNESDefaultComputeJacobian(snes,X,&FDexp,&FDexp,&mstruct,NULL);CHKERRQ(ierr); 2231 ierr = PetscViewerASCIIGetStdout(((PetscObject)snes)->comm,&vstdout);CHKERRQ(ierr); 2232 if (flag_draw || flag_contour) { 2233 ierr = PetscViewerDrawOpen(((PetscObject)snes)->comm,0,"Explicit Jacobians",PETSC_DECIDE,PETSC_DECIDE,300,300,&vdraw);CHKERRQ(ierr); 2234 if (flag_contour) {ierr = PetscViewerPushFormat(vdraw,PETSC_VIEWER_DRAW_CONTOUR);CHKERRQ(ierr);} 2235 } else vdraw = PETSC_NULL; 2236 ierr = PetscViewerASCIIPrintf(vstdout,"Explicit %s\n",flag_operator?"Jacobian":"preconditioning Jacobian");CHKERRQ(ierr); 2237 if (flag) {ierr = MatView(Bexp,vstdout);CHKERRQ(ierr);} 2238 if (vdraw) {ierr = MatView(Bexp,vdraw);CHKERRQ(ierr);} 2239 ierr = PetscViewerASCIIPrintf(vstdout,"Finite difference Jacobian\n");CHKERRQ(ierr); 2240 if (flag) {ierr = MatView(FDexp,vstdout);CHKERRQ(ierr);} 2241 if (vdraw) {ierr = MatView(FDexp,vdraw);CHKERRQ(ierr);} 2242 ierr = MatAYPX(FDexp,-1.0,Bexp,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 2243 ierr = PetscViewerASCIIPrintf(vstdout,"User-provided matrix minus finite difference Jacobian\n");CHKERRQ(ierr); 2244 if (flag) {ierr = MatView(FDexp,vstdout);CHKERRQ(ierr);} 2245 if (vdraw) { /* Always use contour for the difference */ 2246 ierr = PetscViewerPushFormat(vdraw,PETSC_VIEWER_DRAW_CONTOUR);CHKERRQ(ierr); 2247 ierr = MatView(FDexp,vdraw);CHKERRQ(ierr); 2248 ierr = PetscViewerPopFormat(vdraw);CHKERRQ(ierr); 2249 } 2250 if (flag_contour) {ierr = PetscViewerPopFormat(vdraw);CHKERRQ(ierr);} 2251 ierr = PetscViewerDestroy(&vdraw);CHKERRQ(ierr); 2252 ierr = MatDestroy(&Bexp_mine);CHKERRQ(ierr); 2253 ierr = MatDestroy(&FDexp);CHKERRQ(ierr); 2254 } 2255 } 2256 { 2257 PetscBool flag = PETSC_FALSE,flag_display = PETSC_FALSE,flag_draw = PETSC_FALSE,flag_contour = PETSC_FALSE,flag_threshold = PETSC_FALSE; 2258 PetscReal threshold_atol = PETSC_SQRT_MACHINE_EPSILON,threshold_rtol = 10*PETSC_SQRT_MACHINE_EPSILON; 2259 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_coloring",&flag,PETSC_NULL);CHKERRQ(ierr); 2260 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_coloring_display",&flag_display,PETSC_NULL);CHKERRQ(ierr); 2261 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_coloring_draw",&flag_draw,PETSC_NULL);CHKERRQ(ierr); 2262 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_coloring_draw_contour",&flag_contour,PETSC_NULL);CHKERRQ(ierr); 2263 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_coloring_threshold",&flag_threshold,PETSC_NULL);CHKERRQ(ierr); 2264 ierr = PetscOptionsGetReal(((PetscObject)snes)->prefix,"-snes_compare_coloring_threshold_rtol",&threshold_rtol,PETSC_NULL);CHKERRQ(ierr); 2265 ierr = PetscOptionsGetReal(((PetscObject)snes)->prefix,"-snes_compare_coloring_threshold_atol",&threshold_atol,PETSC_NULL);CHKERRQ(ierr); 2266 if (flag || flag_display || flag_draw || flag_contour || flag_threshold) { 2267 Mat Bfd; 2268 MatStructure mstruct; 2269 PetscViewer vdraw,vstdout; 2270 ISColoring iscoloring; 2271 MatFDColoring matfdcoloring; 2272 PetscErrorCode (*func)(SNES,Vec,Vec,void*); 2273 void *funcctx; 2274 PetscReal norm1,norm2,normmax; 2275 2276 ierr = MatDuplicate(*B,MAT_DO_NOT_COPY_VALUES,&Bfd);CHKERRQ(ierr); 2277 ierr = MatGetColoring(Bfd,MATCOLORINGSL,&iscoloring);CHKERRQ(ierr); 2278 ierr = MatFDColoringCreate(Bfd,iscoloring,&matfdcoloring);CHKERRQ(ierr); 2279 ierr = ISColoringDestroy(&iscoloring);CHKERRQ(ierr); 2280 2281 /* This method of getting the function is currently unreliable since it doesn't work for DM local functions. */ 2282 ierr = SNESGetFunction(snes,PETSC_NULL,&func,&funcctx);CHKERRQ(ierr); 2283 ierr = MatFDColoringSetFunction(matfdcoloring,(PetscErrorCode(*)(void))func,funcctx);CHKERRQ(ierr); 2284 ierr = PetscObjectSetOptionsPrefix((PetscObject)matfdcoloring,((PetscObject)snes)->prefix);CHKERRQ(ierr); 2285 ierr = PetscObjectAppendOptionsPrefix((PetscObject)matfdcoloring,"coloring_");CHKERRQ(ierr); 2286 ierr = MatFDColoringSetFromOptions(matfdcoloring);CHKERRQ(ierr); 2287 ierr = MatFDColoringApply(Bfd,matfdcoloring,X,&mstruct,snes);CHKERRQ(ierr); 2288 ierr = MatFDColoringDestroy(&matfdcoloring);CHKERRQ(ierr); 2289 2290 ierr = PetscViewerASCIIGetStdout(((PetscObject)snes)->comm,&vstdout);CHKERRQ(ierr); 2291 if (flag_draw || flag_contour) { 2292 ierr = PetscViewerDrawOpen(((PetscObject)snes)->comm,0,"Colored Jacobians",PETSC_DECIDE,PETSC_DECIDE,300,300,&vdraw);CHKERRQ(ierr); 2293 if (flag_contour) {ierr = PetscViewerPushFormat(vdraw,PETSC_VIEWER_DRAW_CONTOUR);CHKERRQ(ierr);} 2294 } else vdraw = PETSC_NULL; 2295 ierr = PetscViewerASCIIPrintf(vstdout,"Explicit preconditioning Jacobian\n");CHKERRQ(ierr); 2296 if (flag_display) {ierr = MatView(*B,vstdout);CHKERRQ(ierr);} 2297 if (vdraw) {ierr = MatView(*B,vdraw);CHKERRQ(ierr);} 2298 ierr = PetscViewerASCIIPrintf(vstdout,"Colored Finite difference Jacobian\n");CHKERRQ(ierr); 2299 if (flag_display) {ierr = MatView(Bfd,vstdout);CHKERRQ(ierr);} 2300 if (vdraw) {ierr = MatView(Bfd,vdraw);CHKERRQ(ierr);} 2301 ierr = MatAYPX(Bfd,-1.0,*B,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 2302 ierr = MatNorm(Bfd,NORM_1,&norm1);CHKERRQ(ierr); 2303 ierr = MatNorm(Bfd,NORM_FROBENIUS,&norm2);CHKERRQ(ierr); 2304 ierr = MatNorm(Bfd,NORM_MAX,&normmax);CHKERRQ(ierr); 2305 ierr = PetscViewerASCIIPrintf(vstdout,"User-provided matrix minus finite difference Jacobian, norm1=%G normFrob=%G normmax=%G\n",norm1,norm2,normmax);CHKERRQ(ierr); 2306 if (flag_display) {ierr = MatView(Bfd,vstdout);CHKERRQ(ierr);} 2307 if (vdraw) { /* Always use contour for the difference */ 2308 ierr = PetscViewerPushFormat(vdraw,PETSC_VIEWER_DRAW_CONTOUR);CHKERRQ(ierr); 2309 ierr = MatView(Bfd,vdraw);CHKERRQ(ierr); 2310 ierr = PetscViewerPopFormat(vdraw);CHKERRQ(ierr); 2311 } 2312 if (flag_contour) {ierr = PetscViewerPopFormat(vdraw);CHKERRQ(ierr);} 2313 2314 if (flag_threshold) { 2315 PetscInt bs,rstart,rend,i; 2316 ierr = MatGetBlockSize(*B,&bs);CHKERRQ(ierr); 2317 ierr = MatGetOwnershipRange(*B,&rstart,&rend);CHKERRQ(ierr); 2318 for (i=rstart; i<rend; i++) { 2319 const PetscScalar *ba,*ca; 2320 const PetscInt *bj,*cj; 2321 PetscInt bn,cn,j,maxentrycol = -1,maxdiffcol = -1,maxrdiffcol = -1; 2322 PetscReal maxentry = 0,maxdiff = 0,maxrdiff = 0; 2323 ierr = MatGetRow(*B,i,&bn,&bj,&ba);CHKERRQ(ierr); 2324 ierr = MatGetRow(Bfd,i,&cn,&cj,&ca);CHKERRQ(ierr); 2325 if (bn != cn) SETERRQ(((PetscObject)*A)->comm,PETSC_ERR_PLIB,"Unexpected different nonzero pattern in -snes_compare_coloring_threshold"); 2326 for (j=0; j<bn; j++) { 2327 PetscReal rdiff = PetscAbsScalar(ca[j]) / (threshold_atol + threshold_rtol*PetscAbsScalar(ba[j])); 2328 if (PetscAbsScalar(ba[j]) > PetscAbs(maxentry)) { 2329 maxentrycol = bj[j]; 2330 maxentry = PetscRealPart(ba[j]); 2331 } 2332 if (PetscAbsScalar(ca[j]) > PetscAbs(maxdiff)) { 2333 maxdiffcol = bj[j]; 2334 maxdiff = PetscRealPart(ca[j]); 2335 } 2336 if (rdiff > maxrdiff) { 2337 maxrdiffcol = bj[j]; 2338 maxrdiff = rdiff; 2339 } 2340 } 2341 if (maxrdiff > 1) { 2342 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); 2343 for (j=0; j<bn; j++) { 2344 PetscReal rdiff; 2345 rdiff = PetscAbsScalar(ca[j]) / (threshold_atol + threshold_rtol*PetscAbsScalar(ba[j])); 2346 if (rdiff > 1) { 2347 ierr = PetscViewerASCIIPrintf(vstdout," (%D,%G:%G)",bj[j],PetscRealPart(ba[j]),PetscRealPart(ca[j]));CHKERRQ(ierr); 2348 } 2349 } 2350 ierr = PetscViewerASCIIPrintf(vstdout,"\n",i,maxentry,maxdiff,maxrdiff);CHKERRQ(ierr); 2351 } 2352 ierr = MatRestoreRow(*B,i,&bn,&bj,&ba);CHKERRQ(ierr); 2353 ierr = MatRestoreRow(Bfd,i,&cn,&cj,&ca);CHKERRQ(ierr); 2354 } 2355 } 2356 ierr = PetscViewerDestroy(&vdraw);CHKERRQ(ierr); 2357 ierr = MatDestroy(&Bfd);CHKERRQ(ierr); 2358 } 2359 } 2360 PetscFunctionReturn(0); 2361 } 2362 2363 #undef __FUNCT__ 2364 #define __FUNCT__ "SNESSetJacobian" 2365 /*@C 2366 SNESSetJacobian - Sets the function to compute Jacobian as well as the 2367 location to store the matrix. 2368 2369 Logically Collective on SNES and Mat 2370 2371 Input Parameters: 2372 + snes - the SNES context 2373 . A - Jacobian matrix 2374 . B - preconditioner matrix (usually same as the Jacobian) 2375 . func - Jacobian evaluation routine (if PETSC_NULL then SNES retains any previously set value) 2376 - ctx - [optional] user-defined context for private data for the 2377 Jacobian evaluation routine (may be PETSC_NULL) (if PETSC_NULL then SNES retains any previously set value) 2378 2379 Calling sequence of func: 2380 $ func (SNES snes,Vec x,Mat *A,Mat *B,int *flag,void *ctx); 2381 2382 + x - input vector 2383 . A - Jacobian matrix 2384 . B - preconditioner matrix, usually the same as A 2385 . flag - flag indicating information about the preconditioner matrix 2386 structure (same as flag in KSPSetOperators()), one of SAME_NONZERO_PATTERN,DIFFERENT_NONZERO_PATTERN,SAME_PRECONDITIONER 2387 - ctx - [optional] user-defined Jacobian context 2388 2389 Notes: 2390 See KSPSetOperators() for important information about setting the flag 2391 output parameter in the routine func(). Be sure to read this information! 2392 2393 The routine func() takes Mat * as the matrix arguments rather than Mat. 2394 This allows the Jacobian evaluation routine to replace A and/or B with a 2395 completely new new matrix structure (not just different matrix elements) 2396 when appropriate, for instance, if the nonzero structure is changing 2397 throughout the global iterations. 2398 2399 If the A matrix and B matrix are different you must call MatAssemblyBegin/End() on 2400 each matrix. 2401 2402 If using SNESDefaultComputeJacobianColor() to assemble a Jacobian, the ctx argument 2403 must be a MatFDColoring. 2404 2405 Other defect-correction schemes can be used by computing a different matrix in place of the Jacobian. One common 2406 example is to use the "Picard linearization" which only differentiates through the highest order parts of each term. 2407 2408 Level: beginner 2409 2410 .keywords: SNES, nonlinear, set, Jacobian, matrix 2411 2412 .seealso: KSPSetOperators(), SNESSetFunction(), MatMFFDComputeJacobian(), SNESDefaultComputeJacobianColor(), MatStructure 2413 @*/ 2414 PetscErrorCode SNESSetJacobian(SNES snes,Mat A,Mat B,PetscErrorCode (*func)(SNES,Vec,Mat*,Mat*,MatStructure*,void*),void *ctx) 2415 { 2416 PetscErrorCode ierr; 2417 DM dm; 2418 2419 PetscFunctionBegin; 2420 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2421 if (A) PetscValidHeaderSpecific(A,MAT_CLASSID,2); 2422 if (B) PetscValidHeaderSpecific(B,MAT_CLASSID,3); 2423 if (A) PetscCheckSameComm(snes,1,A,2); 2424 if (B) PetscCheckSameComm(snes,1,B,3); 2425 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 2426 ierr = DMSNESSetJacobian(dm,func,ctx);CHKERRQ(ierr); 2427 if (A) { 2428 ierr = PetscObjectReference((PetscObject)A);CHKERRQ(ierr); 2429 ierr = MatDestroy(&snes->jacobian);CHKERRQ(ierr); 2430 snes->jacobian = A; 2431 } 2432 if (B) { 2433 ierr = PetscObjectReference((PetscObject)B);CHKERRQ(ierr); 2434 ierr = MatDestroy(&snes->jacobian_pre);CHKERRQ(ierr); 2435 snes->jacobian_pre = B; 2436 } 2437 PetscFunctionReturn(0); 2438 } 2439 2440 #undef __FUNCT__ 2441 #define __FUNCT__ "SNESGetJacobian" 2442 /*@C 2443 SNESGetJacobian - Returns the Jacobian matrix and optionally the user 2444 provided context for evaluating the Jacobian. 2445 2446 Not Collective, but Mat object will be parallel if SNES object is 2447 2448 Input Parameter: 2449 . snes - the nonlinear solver context 2450 2451 Output Parameters: 2452 + A - location to stash Jacobian matrix (or PETSC_NULL) 2453 . B - location to stash preconditioner matrix (or PETSC_NULL) 2454 . func - location to put Jacobian function (or PETSC_NULL) 2455 - ctx - location to stash Jacobian ctx (or PETSC_NULL) 2456 2457 Level: advanced 2458 2459 .seealso: SNESSetJacobian(), SNESComputeJacobian() 2460 @*/ 2461 PetscErrorCode SNESGetJacobian(SNES snes,Mat *A,Mat *B,PetscErrorCode (**func)(SNES,Vec,Mat*,Mat*,MatStructure*,void*),void **ctx) 2462 { 2463 PetscErrorCode ierr; 2464 DM dm; 2465 DMSNES sdm; 2466 2467 PetscFunctionBegin; 2468 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2469 if (A) *A = snes->jacobian; 2470 if (B) *B = snes->jacobian_pre; 2471 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 2472 ierr = DMGetDMSNES(dm,&sdm);CHKERRQ(ierr); 2473 if (func) *func = sdm->ops->computejacobian; 2474 if (ctx) *ctx = sdm->jacobianctx; 2475 PetscFunctionReturn(0); 2476 } 2477 2478 /* ----- Routines to initialize and destroy a nonlinear solver ---- */ 2479 2480 #undef __FUNCT__ 2481 #define __FUNCT__ "SNESSetUp" 2482 /*@ 2483 SNESSetUp - Sets up the internal data structures for the later use 2484 of a nonlinear solver. 2485 2486 Collective on SNES 2487 2488 Input Parameters: 2489 . snes - the SNES context 2490 2491 Notes: 2492 For basic use of the SNES solvers the user need not explicitly call 2493 SNESSetUp(), since these actions will automatically occur during 2494 the call to SNESSolve(). However, if one wishes to control this 2495 phase separately, SNESSetUp() should be called after SNESCreate() 2496 and optional routines of the form SNESSetXXX(), but before SNESSolve(). 2497 2498 Level: advanced 2499 2500 .keywords: SNES, nonlinear, setup 2501 2502 .seealso: SNESCreate(), SNESSolve(), SNESDestroy() 2503 @*/ 2504 PetscErrorCode SNESSetUp(SNES snes) 2505 { 2506 PetscErrorCode ierr; 2507 DM dm; 2508 DMSNES sdm; 2509 SNESLineSearch linesearch, pclinesearch; 2510 void *lsprectx,*lspostctx; 2511 SNESLineSearchPreCheckFunc lsprefunc; 2512 SNESLineSearchPostCheckFunc lspostfunc; 2513 PetscErrorCode (*func)(SNES,Vec,Vec,void*); 2514 Vec f,fpc; 2515 void *funcctx; 2516 PetscErrorCode (*jac)(SNES,Vec,Mat*,Mat*,MatStructure*,void*); 2517 void *jacctx,*appctx; 2518 Mat A,B; 2519 2520 PetscFunctionBegin; 2521 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2522 if (snes->setupcalled) PetscFunctionReturn(0); 2523 2524 if (!((PetscObject)snes)->type_name) { 2525 ierr = SNESSetType(snes,SNESLS);CHKERRQ(ierr); 2526 } 2527 2528 ierr = SNESGetFunction(snes,&snes->vec_func,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr); 2529 if (snes->vec_func == snes->vec_sol) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_IDN,"Solution vector cannot be function vector"); 2530 if (snes->vec_rhs == snes->vec_sol) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_IDN,"Solution vector cannot be right hand side vector"); 2531 2532 if (!snes->vec_sol_update /* && snes->vec_sol */) { 2533 ierr = VecDuplicate(snes->vec_sol,&snes->vec_sol_update);CHKERRQ(ierr); 2534 ierr = PetscLogObjectParent(snes,snes->vec_sol_update);CHKERRQ(ierr); 2535 } 2536 2537 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 2538 ierr = DMShellSetGlobalVector(snes->dm,snes->vec_sol);CHKERRQ(ierr); 2539 ierr = DMGetDMSNES(dm,&sdm);CHKERRQ(ierr); 2540 if (!sdm->ops->computefunction) SETERRQ(((PetscObject)dm)->comm,PETSC_ERR_ARG_WRONGSTATE,"Function never provided to SNES object"); 2541 if (!sdm->ops->computejacobian) { 2542 ierr = DMSNESSetJacobian(dm,SNESDefaultComputeJacobianColor,PETSC_NULL);CHKERRQ(ierr); 2543 } 2544 if (!snes->vec_func) { 2545 ierr = DMCreateGlobalVector(dm,&snes->vec_func);CHKERRQ(ierr); 2546 } 2547 2548 if (!snes->ksp) {ierr = SNESGetKSP(snes, &snes->ksp);CHKERRQ(ierr);} 2549 2550 if (!snes->linesearch) {ierr = SNESGetSNESLineSearch(snes, &snes->linesearch);} 2551 2552 if (snes->pc && (snes->pcside == PC_LEFT)) snes->mf = PETSC_TRUE; 2553 2554 if (snes->mf) { 2555 ierr = SNESSetUpMatrixFree_Private(snes, snes->mf_operator, snes->mf_version);CHKERRQ(ierr); 2556 } 2557 2558 2559 if (snes->ops->usercompute && !snes->user) { 2560 ierr = (*snes->ops->usercompute)(snes,(void**)&snes->user);CHKERRQ(ierr); 2561 } 2562 2563 if (snes->pc) { 2564 /* copy the DM over */ 2565 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 2566 ierr = SNESSetDM(snes->pc,dm);CHKERRQ(ierr); 2567 2568 ierr = SNESGetFunction(snes,&f,&func,&funcctx);CHKERRQ(ierr); 2569 ierr = VecDuplicate(f,&fpc);CHKERRQ(ierr); 2570 ierr = SNESSetFunction(snes->pc,fpc,func,funcctx);CHKERRQ(ierr); 2571 ierr = SNESGetJacobian(snes,&A,&B,&jac,&jacctx);CHKERRQ(ierr); 2572 ierr = SNESSetJacobian(snes->pc,A,B,jac,jacctx);CHKERRQ(ierr); 2573 ierr = SNESGetApplicationContext(snes,&appctx);CHKERRQ(ierr); 2574 ierr = SNESSetApplicationContext(snes->pc,appctx);CHKERRQ(ierr); 2575 ierr = VecDestroy(&fpc);CHKERRQ(ierr); 2576 2577 /* copy the function pointers over */ 2578 ierr = PetscObjectCopyFortranFunctionPointers((PetscObject)snes,(PetscObject)snes->pc);CHKERRQ(ierr); 2579 2580 /* default to 1 iteration */ 2581 ierr = SNESSetTolerances(snes->pc, 0.0, 0.0, 0.0, 1, snes->pc->max_funcs);CHKERRQ(ierr); 2582 ierr = SNESSetNormType(snes->pc, SNES_NORM_FINAL_ONLY);CHKERRQ(ierr); 2583 ierr = SNESSetFromOptions(snes->pc);CHKERRQ(ierr); 2584 2585 /* copy the line search context over */ 2586 ierr = SNESGetSNESLineSearch(snes,&linesearch);CHKERRQ(ierr); 2587 ierr = SNESGetSNESLineSearch(snes->pc,&pclinesearch);CHKERRQ(ierr); 2588 ierr = SNESLineSearchGetPreCheck(linesearch,&lsprefunc,&lsprectx);CHKERRQ(ierr); 2589 ierr = SNESLineSearchGetPostCheck(linesearch,&lspostfunc,&lspostctx);CHKERRQ(ierr); 2590 ierr = SNESLineSearchSetPreCheck(pclinesearch,lsprefunc,lsprectx);CHKERRQ(ierr); 2591 ierr = SNESLineSearchSetPostCheck(pclinesearch,lspostfunc,lspostctx);CHKERRQ(ierr); 2592 ierr = PetscObjectCopyFortranFunctionPointers((PetscObject)linesearch, (PetscObject)pclinesearch);CHKERRQ(ierr); 2593 } 2594 2595 if (snes->ops->setup) { 2596 ierr = (*snes->ops->setup)(snes);CHKERRQ(ierr); 2597 } 2598 2599 snes->setupcalled = PETSC_TRUE; 2600 PetscFunctionReturn(0); 2601 } 2602 2603 #undef __FUNCT__ 2604 #define __FUNCT__ "SNESReset" 2605 /*@ 2606 SNESReset - Resets a SNES context to the snessetupcalled = 0 state and removes any allocated Vecs and Mats 2607 2608 Collective on SNES 2609 2610 Input Parameter: 2611 . snes - iterative context obtained from SNESCreate() 2612 2613 Level: intermediate 2614 2615 Notes: Also calls the application context destroy routine set with SNESSetComputeApplicationContext() 2616 2617 .keywords: SNES, destroy 2618 2619 .seealso: SNESCreate(), SNESSetUp(), SNESSolve() 2620 @*/ 2621 PetscErrorCode SNESReset(SNES snes) 2622 { 2623 PetscErrorCode ierr; 2624 2625 PetscFunctionBegin; 2626 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2627 if (snes->ops->userdestroy && snes->user) { 2628 ierr = (*snes->ops->userdestroy)((void**)&snes->user);CHKERRQ(ierr); 2629 snes->user = PETSC_NULL; 2630 } 2631 if (snes->pc) { 2632 ierr = SNESReset(snes->pc);CHKERRQ(ierr); 2633 } 2634 2635 if (snes->ops->reset) { 2636 ierr = (*snes->ops->reset)(snes);CHKERRQ(ierr); 2637 } 2638 if (snes->ksp) { 2639 ierr = KSPReset(snes->ksp);CHKERRQ(ierr); 2640 } 2641 2642 if (snes->linesearch) { 2643 ierr = SNESLineSearchReset(snes->linesearch);CHKERRQ(ierr); 2644 } 2645 2646 ierr = VecDestroy(&snes->vec_rhs);CHKERRQ(ierr); 2647 ierr = VecDestroy(&snes->vec_sol);CHKERRQ(ierr); 2648 ierr = VecDestroy(&snes->vec_sol_update);CHKERRQ(ierr); 2649 ierr = VecDestroy(&snes->vec_func);CHKERRQ(ierr); 2650 ierr = MatDestroy(&snes->jacobian);CHKERRQ(ierr); 2651 ierr = MatDestroy(&snes->jacobian_pre);CHKERRQ(ierr); 2652 ierr = VecDestroyVecs(snes->nwork,&snes->work);CHKERRQ(ierr); 2653 ierr = VecDestroyVecs(snes->nvwork,&snes->vwork);CHKERRQ(ierr); 2654 snes->nwork = snes->nvwork = 0; 2655 snes->setupcalled = PETSC_FALSE; 2656 PetscFunctionReturn(0); 2657 } 2658 2659 #undef __FUNCT__ 2660 #define __FUNCT__ "SNESDestroy" 2661 /*@ 2662 SNESDestroy - Destroys the nonlinear solver context that was created 2663 with SNESCreate(). 2664 2665 Collective on SNES 2666 2667 Input Parameter: 2668 . snes - the SNES context 2669 2670 Level: beginner 2671 2672 .keywords: SNES, nonlinear, destroy 2673 2674 .seealso: SNESCreate(), SNESSolve() 2675 @*/ 2676 PetscErrorCode SNESDestroy(SNES *snes) 2677 { 2678 PetscErrorCode ierr; 2679 2680 PetscFunctionBegin; 2681 if (!*snes) PetscFunctionReturn(0); 2682 PetscValidHeaderSpecific((*snes),SNES_CLASSID,1); 2683 if (--((PetscObject)(*snes))->refct > 0) {*snes = 0; PetscFunctionReturn(0);} 2684 2685 ierr = SNESReset((*snes));CHKERRQ(ierr); 2686 ierr = SNESDestroy(&(*snes)->pc);CHKERRQ(ierr); 2687 2688 /* if memory was published with AMS then destroy it */ 2689 ierr = PetscObjectDepublish((*snes));CHKERRQ(ierr); 2690 if ((*snes)->ops->destroy) {ierr = (*((*snes))->ops->destroy)((*snes));CHKERRQ(ierr);} 2691 2692 ierr = DMDestroy(&(*snes)->dm);CHKERRQ(ierr); 2693 ierr = KSPDestroy(&(*snes)->ksp);CHKERRQ(ierr); 2694 ierr = SNESLineSearchDestroy(&(*snes)->linesearch);CHKERRQ(ierr); 2695 2696 ierr = PetscFree((*snes)->kspconvctx);CHKERRQ(ierr); 2697 if ((*snes)->ops->convergeddestroy) { 2698 ierr = (*(*snes)->ops->convergeddestroy)((*snes)->cnvP);CHKERRQ(ierr); 2699 } 2700 if ((*snes)->conv_malloc) { 2701 ierr = PetscFree((*snes)->conv_hist);CHKERRQ(ierr); 2702 ierr = PetscFree((*snes)->conv_hist_its);CHKERRQ(ierr); 2703 } 2704 ierr = SNESMonitorCancel((*snes));CHKERRQ(ierr); 2705 ierr = PetscHeaderDestroy(snes);CHKERRQ(ierr); 2706 PetscFunctionReturn(0); 2707 } 2708 2709 /* ----------- Routines to set solver parameters ---------- */ 2710 2711 #undef __FUNCT__ 2712 #define __FUNCT__ "SNESSetLagPreconditioner" 2713 /*@ 2714 SNESSetLagPreconditioner - Determines when the preconditioner is rebuilt in the nonlinear solve. 2715 2716 Logically Collective on SNES 2717 2718 Input Parameters: 2719 + snes - the SNES context 2720 - lag - -1 indicates NEVER rebuild, 1 means rebuild every time the Jacobian is computed within a single nonlinear solve, 2 means every second time 2721 the Jacobian is built etc. -2 indicates rebuild preconditioner at next chance but then never rebuild after that 2722 2723 Options Database Keys: 2724 . -snes_lag_preconditioner <lag> 2725 2726 Notes: 2727 The default is 1 2728 The preconditioner is ALWAYS built in the first iteration of a nonlinear solve unless lag is -1 2729 If -1 is used before the very first nonlinear solve the preconditioner is still built because there is no previous preconditioner to use 2730 2731 Level: intermediate 2732 2733 .keywords: SNES, nonlinear, set, convergence, tolerances 2734 2735 .seealso: SNESSetTrustRegionTolerance(), SNESGetLagPreconditioner(), SNESSetLagJacobian(), SNESGetLagJacobian() 2736 2737 @*/ 2738 PetscErrorCode SNESSetLagPreconditioner(SNES snes,PetscInt lag) 2739 { 2740 PetscFunctionBegin; 2741 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2742 if (lag < -2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Lag must be -2, -1, 1 or greater"); 2743 if (!lag) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Lag cannot be 0"); 2744 PetscValidLogicalCollectiveInt(snes,lag,2); 2745 snes->lagpreconditioner = lag; 2746 PetscFunctionReturn(0); 2747 } 2748 2749 #undef __FUNCT__ 2750 #define __FUNCT__ "SNESSetGridSequence" 2751 /*@ 2752 SNESSetGridSequence - sets the number of steps of grid sequencing that SNES does 2753 2754 Logically Collective on SNES 2755 2756 Input Parameters: 2757 + snes - the SNES context 2758 - steps - the number of refinements to do, defaults to 0 2759 2760 Options Database Keys: 2761 . -snes_grid_sequence <steps> 2762 2763 Level: intermediate 2764 2765 Notes: 2766 Use SNESGetSolution() to extract the fine grid solution after grid sequencing. 2767 2768 .keywords: SNES, nonlinear, set, convergence, tolerances 2769 2770 .seealso: SNESSetTrustRegionTolerance(), SNESGetLagPreconditioner(), SNESSetLagJacobian(), SNESGetLagJacobian() 2771 2772 @*/ 2773 PetscErrorCode SNESSetGridSequence(SNES snes,PetscInt steps) 2774 { 2775 PetscFunctionBegin; 2776 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2777 PetscValidLogicalCollectiveInt(snes,steps,2); 2778 snes->gridsequence = steps; 2779 PetscFunctionReturn(0); 2780 } 2781 2782 #undef __FUNCT__ 2783 #define __FUNCT__ "SNESGetLagPreconditioner" 2784 /*@ 2785 SNESGetLagPreconditioner - Indicates how often the preconditioner is rebuilt 2786 2787 Not Collective 2788 2789 Input Parameter: 2790 . snes - the SNES context 2791 2792 Output Parameter: 2793 . lag - -1 indicates NEVER rebuild, 1 means rebuild every time the Jacobian is computed within a single nonlinear solve, 2 means every second time 2794 the Jacobian is built etc. -2 indicates rebuild preconditioner at next chance but then never rebuild after that 2795 2796 Options Database Keys: 2797 . -snes_lag_preconditioner <lag> 2798 2799 Notes: 2800 The default is 1 2801 The preconditioner is ALWAYS built in the first iteration of a nonlinear solve unless lag is -1 2802 2803 Level: intermediate 2804 2805 .keywords: SNES, nonlinear, set, convergence, tolerances 2806 2807 .seealso: SNESSetTrustRegionTolerance(), SNESSetLagPreconditioner() 2808 2809 @*/ 2810 PetscErrorCode SNESGetLagPreconditioner(SNES snes,PetscInt *lag) 2811 { 2812 PetscFunctionBegin; 2813 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2814 *lag = snes->lagpreconditioner; 2815 PetscFunctionReturn(0); 2816 } 2817 2818 #undef __FUNCT__ 2819 #define __FUNCT__ "SNESSetLagJacobian" 2820 /*@ 2821 SNESSetLagJacobian - Determines when the Jacobian is rebuilt in the nonlinear solve. See SNESSetLagPreconditioner() for determining how 2822 often the preconditioner is rebuilt. 2823 2824 Logically Collective on SNES 2825 2826 Input Parameters: 2827 + snes - the SNES context 2828 - lag - -1 indicates NEVER rebuild, 1 means rebuild every time the Jacobian is computed within a single nonlinear solve, 2 means every second time 2829 the Jacobian is built etc. -2 means rebuild at next chance but then never again 2830 2831 Options Database Keys: 2832 . -snes_lag_jacobian <lag> 2833 2834 Notes: 2835 The default is 1 2836 The Jacobian is ALWAYS built in the first iteration of a nonlinear solve unless lag is -1 2837 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 2838 at the next Newton step but never again (unless it is reset to another value) 2839 2840 Level: intermediate 2841 2842 .keywords: SNES, nonlinear, set, convergence, tolerances 2843 2844 .seealso: SNESSetTrustRegionTolerance(), SNESGetLagPreconditioner(), SNESSetLagPreconditioner(), SNESGetLagJacobian() 2845 2846 @*/ 2847 PetscErrorCode SNESSetLagJacobian(SNES snes,PetscInt lag) 2848 { 2849 PetscFunctionBegin; 2850 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2851 if (lag < -2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Lag must be -2, -1, 1 or greater"); 2852 if (!lag) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Lag cannot be 0"); 2853 PetscValidLogicalCollectiveInt(snes,lag,2); 2854 snes->lagjacobian = lag; 2855 PetscFunctionReturn(0); 2856 } 2857 2858 #undef __FUNCT__ 2859 #define __FUNCT__ "SNESGetLagJacobian" 2860 /*@ 2861 SNESGetLagJacobian - Indicates how often the Jacobian is rebuilt. See SNESGetLagPreconditioner() to determine when the preconditioner is rebuilt 2862 2863 Not Collective 2864 2865 Input Parameter: 2866 . snes - the SNES context 2867 2868 Output Parameter: 2869 . lag - -1 indicates NEVER rebuild, 1 means rebuild every time the Jacobian is computed within a single nonlinear solve, 2 means every second time 2870 the Jacobian is built etc. 2871 2872 Options Database Keys: 2873 . -snes_lag_jacobian <lag> 2874 2875 Notes: 2876 The default is 1 2877 The jacobian is ALWAYS built in the first iteration of a nonlinear solve unless lag is -1 2878 2879 Level: intermediate 2880 2881 .keywords: SNES, nonlinear, set, convergence, tolerances 2882 2883 .seealso: SNESSetTrustRegionTolerance(), SNESSetLagJacobian(), SNESSetLagPreconditioner(), SNESGetLagPreconditioner() 2884 2885 @*/ 2886 PetscErrorCode SNESGetLagJacobian(SNES snes,PetscInt *lag) 2887 { 2888 PetscFunctionBegin; 2889 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2890 *lag = snes->lagjacobian; 2891 PetscFunctionReturn(0); 2892 } 2893 2894 #undef __FUNCT__ 2895 #define __FUNCT__ "SNESSetTolerances" 2896 /*@ 2897 SNESSetTolerances - Sets various parameters used in convergence tests. 2898 2899 Logically Collective on SNES 2900 2901 Input Parameters: 2902 + snes - the SNES context 2903 . abstol - absolute convergence tolerance 2904 . rtol - relative convergence tolerance 2905 . stol - convergence tolerance in terms of the norm of the change in the solution between steps, || delta x || < stol*|| x || 2906 . maxit - maximum number of iterations 2907 - maxf - maximum number of function evaluations 2908 2909 Options Database Keys: 2910 + -snes_atol <abstol> - Sets abstol 2911 . -snes_rtol <rtol> - Sets rtol 2912 . -snes_stol <stol> - Sets stol 2913 . -snes_max_it <maxit> - Sets maxit 2914 - -snes_max_funcs <maxf> - Sets maxf 2915 2916 Notes: 2917 The default maximum number of iterations is 50. 2918 The default maximum number of function evaluations is 1000. 2919 2920 Level: intermediate 2921 2922 .keywords: SNES, nonlinear, set, convergence, tolerances 2923 2924 .seealso: SNESSetTrustRegionTolerance() 2925 @*/ 2926 PetscErrorCode SNESSetTolerances(SNES snes,PetscReal abstol,PetscReal rtol,PetscReal stol,PetscInt maxit,PetscInt maxf) 2927 { 2928 PetscFunctionBegin; 2929 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2930 PetscValidLogicalCollectiveReal(snes,abstol,2); 2931 PetscValidLogicalCollectiveReal(snes,rtol,3); 2932 PetscValidLogicalCollectiveReal(snes,stol,4); 2933 PetscValidLogicalCollectiveInt(snes,maxit,5); 2934 PetscValidLogicalCollectiveInt(snes,maxf,6); 2935 2936 if (abstol != PETSC_DEFAULT) { 2937 if (abstol < 0.0) SETERRQ1(((PetscObject)snes)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Absolute tolerance %G must be non-negative",abstol); 2938 snes->abstol = abstol; 2939 } 2940 if (rtol != PETSC_DEFAULT) { 2941 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); 2942 snes->rtol = rtol; 2943 } 2944 if (stol != PETSC_DEFAULT) { 2945 if (stol < 0.0) SETERRQ1(((PetscObject)snes)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Step tolerance %G must be non-negative",stol); 2946 snes->stol = stol; 2947 } 2948 if (maxit != PETSC_DEFAULT) { 2949 if (maxit < 0) SETERRQ1(((PetscObject)snes)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Maximum number of iterations %D must be non-negative",maxit); 2950 snes->max_its = maxit; 2951 } 2952 if (maxf != PETSC_DEFAULT) { 2953 if (maxf < 0) SETERRQ1(((PetscObject)snes)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Maximum number of function evaluations %D must be non-negative",maxf); 2954 snes->max_funcs = maxf; 2955 } 2956 snes->tolerancesset = PETSC_TRUE; 2957 PetscFunctionReturn(0); 2958 } 2959 2960 #undef __FUNCT__ 2961 #define __FUNCT__ "SNESGetTolerances" 2962 /*@ 2963 SNESGetTolerances - Gets various parameters used in convergence tests. 2964 2965 Not Collective 2966 2967 Input Parameters: 2968 + snes - the SNES context 2969 . atol - absolute convergence tolerance 2970 . rtol - relative convergence tolerance 2971 . stol - convergence tolerance in terms of the norm 2972 of the change in the solution between steps 2973 . maxit - maximum number of iterations 2974 - maxf - maximum number of function evaluations 2975 2976 Notes: 2977 The user can specify PETSC_NULL for any parameter that is not needed. 2978 2979 Level: intermediate 2980 2981 .keywords: SNES, nonlinear, get, convergence, tolerances 2982 2983 .seealso: SNESSetTolerances() 2984 @*/ 2985 PetscErrorCode SNESGetTolerances(SNES snes,PetscReal *atol,PetscReal *rtol,PetscReal *stol,PetscInt *maxit,PetscInt *maxf) 2986 { 2987 PetscFunctionBegin; 2988 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2989 if (atol) *atol = snes->abstol; 2990 if (rtol) *rtol = snes->rtol; 2991 if (stol) *stol = snes->stol; 2992 if (maxit) *maxit = snes->max_its; 2993 if (maxf) *maxf = snes->max_funcs; 2994 PetscFunctionReturn(0); 2995 } 2996 2997 #undef __FUNCT__ 2998 #define __FUNCT__ "SNESSetTrustRegionTolerance" 2999 /*@ 3000 SNESSetTrustRegionTolerance - Sets the trust region parameter tolerance. 3001 3002 Logically Collective on SNES 3003 3004 Input Parameters: 3005 + snes - the SNES context 3006 - tol - tolerance 3007 3008 Options Database Key: 3009 . -snes_trtol <tol> - Sets tol 3010 3011 Level: intermediate 3012 3013 .keywords: SNES, nonlinear, set, trust region, tolerance 3014 3015 .seealso: SNESSetTolerances() 3016 @*/ 3017 PetscErrorCode SNESSetTrustRegionTolerance(SNES snes,PetscReal tol) 3018 { 3019 PetscFunctionBegin; 3020 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3021 PetscValidLogicalCollectiveReal(snes,tol,2); 3022 snes->deltatol = tol; 3023 PetscFunctionReturn(0); 3024 } 3025 3026 /* 3027 Duplicate the lg monitors for SNES from KSP; for some reason with 3028 dynamic libraries things don't work under Sun4 if we just use 3029 macros instead of functions 3030 */ 3031 #undef __FUNCT__ 3032 #define __FUNCT__ "SNESMonitorLGResidualNorm" 3033 PetscErrorCode SNESMonitorLGResidualNorm(SNES snes,PetscInt it,PetscReal norm,void *ctx) 3034 { 3035 PetscErrorCode ierr; 3036 3037 PetscFunctionBegin; 3038 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3039 ierr = KSPMonitorLGResidualNorm((KSP)snes,it,norm,ctx);CHKERRQ(ierr); 3040 PetscFunctionReturn(0); 3041 } 3042 3043 #undef __FUNCT__ 3044 #define __FUNCT__ "SNESMonitorLGCreate" 3045 PetscErrorCode SNESMonitorLGCreate(const char host[],const char label[],int x,int y,int m,int n,PetscDrawLG *draw) 3046 { 3047 PetscErrorCode ierr; 3048 3049 PetscFunctionBegin; 3050 ierr = KSPMonitorLGResidualNormCreate(host,label,x,y,m,n,draw);CHKERRQ(ierr); 3051 PetscFunctionReturn(0); 3052 } 3053 3054 #undef __FUNCT__ 3055 #define __FUNCT__ "SNESMonitorLGDestroy" 3056 PetscErrorCode SNESMonitorLGDestroy(PetscDrawLG *draw) 3057 { 3058 PetscErrorCode ierr; 3059 3060 PetscFunctionBegin; 3061 ierr = KSPMonitorLGResidualNormDestroy(draw);CHKERRQ(ierr); 3062 PetscFunctionReturn(0); 3063 } 3064 3065 extern PetscErrorCode SNESMonitorRange_Private(SNES,PetscInt,PetscReal*); 3066 #undef __FUNCT__ 3067 #define __FUNCT__ "SNESMonitorLGRange" 3068 PetscErrorCode SNESMonitorLGRange(SNES snes,PetscInt n,PetscReal rnorm,void *monctx) 3069 { 3070 PetscDrawLG lg; 3071 PetscErrorCode ierr; 3072 PetscReal x,y,per; 3073 PetscViewer v = (PetscViewer)monctx; 3074 static PetscReal prev; /* should be in the context */ 3075 PetscDraw draw; 3076 3077 PetscFunctionBegin; 3078 ierr = PetscViewerDrawGetDrawLG(v,0,&lg);CHKERRQ(ierr); 3079 if (!n) {ierr = PetscDrawLGReset(lg);CHKERRQ(ierr);} 3080 ierr = PetscDrawLGGetDraw(lg,&draw);CHKERRQ(ierr); 3081 ierr = PetscDrawSetTitle(draw,"Residual norm");CHKERRQ(ierr); 3082 x = (PetscReal) n; 3083 if (rnorm > 0.0) y = log10(rnorm); else y = -15.0; 3084 ierr = PetscDrawLGAddPoint(lg,&x,&y);CHKERRQ(ierr); 3085 if (n < 20 || !(n % 5)) { 3086 ierr = PetscDrawLGDraw(lg);CHKERRQ(ierr); 3087 } 3088 3089 ierr = PetscViewerDrawGetDrawLG(v,1,&lg);CHKERRQ(ierr); 3090 if (!n) {ierr = PetscDrawLGReset(lg);CHKERRQ(ierr);} 3091 ierr = PetscDrawLGGetDraw(lg,&draw);CHKERRQ(ierr); 3092 ierr = PetscDrawSetTitle(draw,"% elemts > .2*max elemt");CHKERRQ(ierr); 3093 ierr = SNESMonitorRange_Private(snes,n,&per);CHKERRQ(ierr); 3094 x = (PetscReal) n; 3095 y = 100.0*per; 3096 ierr = PetscDrawLGAddPoint(lg,&x,&y);CHKERRQ(ierr); 3097 if (n < 20 || !(n % 5)) { 3098 ierr = PetscDrawLGDraw(lg);CHKERRQ(ierr); 3099 } 3100 3101 ierr = PetscViewerDrawGetDrawLG(v,2,&lg);CHKERRQ(ierr); 3102 if (!n) {prev = rnorm;ierr = PetscDrawLGReset(lg);CHKERRQ(ierr);} 3103 ierr = PetscDrawLGGetDraw(lg,&draw);CHKERRQ(ierr); 3104 ierr = PetscDrawSetTitle(draw,"(norm -oldnorm)/oldnorm");CHKERRQ(ierr); 3105 x = (PetscReal) n; 3106 y = (prev - rnorm)/prev; 3107 ierr = PetscDrawLGAddPoint(lg,&x,&y);CHKERRQ(ierr); 3108 if (n < 20 || !(n % 5)) { 3109 ierr = PetscDrawLGDraw(lg);CHKERRQ(ierr); 3110 } 3111 3112 ierr = PetscViewerDrawGetDrawLG(v,3,&lg);CHKERRQ(ierr); 3113 if (!n) {ierr = PetscDrawLGReset(lg);CHKERRQ(ierr);} 3114 ierr = PetscDrawLGGetDraw(lg,&draw);CHKERRQ(ierr); 3115 ierr = PetscDrawSetTitle(draw,"(norm -oldnorm)/oldnorm*(% > .2 max)");CHKERRQ(ierr); 3116 x = (PetscReal) n; 3117 y = (prev - rnorm)/(prev*per); 3118 if (n > 2) { /*skip initial crazy value */ 3119 ierr = PetscDrawLGAddPoint(lg,&x,&y);CHKERRQ(ierr); 3120 } 3121 if (n < 20 || !(n % 5)) { 3122 ierr = PetscDrawLGDraw(lg);CHKERRQ(ierr); 3123 } 3124 prev = rnorm; 3125 PetscFunctionReturn(0); 3126 } 3127 3128 #undef __FUNCT__ 3129 #define __FUNCT__ "SNESMonitor" 3130 /*@ 3131 SNESMonitor - runs the user provided monitor routines, if they exist 3132 3133 Collective on SNES 3134 3135 Input Parameters: 3136 + snes - nonlinear solver context obtained from SNESCreate() 3137 . iter - iteration number 3138 - rnorm - relative norm of the residual 3139 3140 Notes: 3141 This routine is called by the SNES implementations. 3142 It does not typically need to be called by the user. 3143 3144 Level: developer 3145 3146 .seealso: SNESMonitorSet() 3147 @*/ 3148 PetscErrorCode SNESMonitor(SNES snes,PetscInt iter,PetscReal rnorm) 3149 { 3150 PetscErrorCode ierr; 3151 PetscInt i,n = snes->numbermonitors; 3152 3153 PetscFunctionBegin; 3154 for (i=0; i<n; i++) { 3155 ierr = (*snes->monitor[i])(snes,iter,rnorm,snes->monitorcontext[i]);CHKERRQ(ierr); 3156 } 3157 PetscFunctionReturn(0); 3158 } 3159 3160 /* ------------ Routines to set performance monitoring options ----------- */ 3161 3162 #undef __FUNCT__ 3163 #define __FUNCT__ "SNESMonitorSet" 3164 /*@C 3165 SNESMonitorSet - Sets an ADDITIONAL function that is to be used at every 3166 iteration of the nonlinear solver to display the iteration's 3167 progress. 3168 3169 Logically Collective on SNES 3170 3171 Input Parameters: 3172 + snes - the SNES context 3173 . func - monitoring routine 3174 . mctx - [optional] user-defined context for private data for the 3175 monitor routine (use PETSC_NULL if no context is desired) 3176 - monitordestroy - [optional] routine that frees monitor context 3177 (may be PETSC_NULL) 3178 3179 Calling sequence of func: 3180 $ int func(SNES snes,PetscInt its, PetscReal norm,void *mctx) 3181 3182 + snes - the SNES context 3183 . its - iteration number 3184 . norm - 2-norm function value (may be estimated) 3185 - mctx - [optional] monitoring context 3186 3187 Options Database Keys: 3188 + -snes_monitor - sets SNESMonitorDefault() 3189 . -snes_monitor_lg_residualnorm - sets line graph monitor, 3190 uses SNESMonitorLGCreate() 3191 - -snes_monitor_cancel - cancels all monitors that have 3192 been hardwired into a code by 3193 calls to SNESMonitorSet(), but 3194 does not cancel those set via 3195 the options database. 3196 3197 Notes: 3198 Several different monitoring routines may be set by calling 3199 SNESMonitorSet() multiple times; all will be called in the 3200 order in which they were set. 3201 3202 Fortran notes: Only a single monitor function can be set for each SNES object 3203 3204 Level: intermediate 3205 3206 .keywords: SNES, nonlinear, set, monitor 3207 3208 .seealso: SNESMonitorDefault(), SNESMonitorCancel() 3209 @*/ 3210 PetscErrorCode SNESMonitorSet(SNES snes,PetscErrorCode (*monitor)(SNES,PetscInt,PetscReal,void*),void *mctx,PetscErrorCode (*monitordestroy)(void**)) 3211 { 3212 PetscInt i; 3213 PetscErrorCode ierr; 3214 3215 PetscFunctionBegin; 3216 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3217 if (snes->numbermonitors >= MAXSNESMONITORS) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Too many monitors set"); 3218 for (i=0; i<snes->numbermonitors;i++) { 3219 if (monitor == snes->monitor[i] && monitordestroy == snes->monitordestroy[i] && mctx == snes->monitorcontext[i]) { 3220 if (monitordestroy) { 3221 ierr = (*monitordestroy)(&mctx);CHKERRQ(ierr); 3222 } 3223 PetscFunctionReturn(0); 3224 } 3225 } 3226 snes->monitor[snes->numbermonitors] = monitor; 3227 snes->monitordestroy[snes->numbermonitors] = monitordestroy; 3228 snes->monitorcontext[snes->numbermonitors++] = (void*)mctx; 3229 PetscFunctionReturn(0); 3230 } 3231 3232 #undef __FUNCT__ 3233 #define __FUNCT__ "SNESMonitorCancel" 3234 /*@C 3235 SNESMonitorCancel - Clears all the monitor functions for a SNES object. 3236 3237 Logically Collective on SNES 3238 3239 Input Parameters: 3240 . snes - the SNES context 3241 3242 Options Database Key: 3243 . -snes_monitor_cancel - cancels all monitors that have been hardwired 3244 into a code by calls to SNESMonitorSet(), but does not cancel those 3245 set via the options database 3246 3247 Notes: 3248 There is no way to clear one specific monitor from a SNES object. 3249 3250 Level: intermediate 3251 3252 .keywords: SNES, nonlinear, set, monitor 3253 3254 .seealso: SNESMonitorDefault(), SNESMonitorSet() 3255 @*/ 3256 PetscErrorCode SNESMonitorCancel(SNES snes) 3257 { 3258 PetscErrorCode ierr; 3259 PetscInt i; 3260 3261 PetscFunctionBegin; 3262 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3263 for (i=0; i<snes->numbermonitors; i++) { 3264 if (snes->monitordestroy[i]) { 3265 ierr = (*snes->monitordestroy[i])(&snes->monitorcontext[i]);CHKERRQ(ierr); 3266 } 3267 } 3268 snes->numbermonitors = 0; 3269 PetscFunctionReturn(0); 3270 } 3271 3272 #undef __FUNCT__ 3273 #define __FUNCT__ "SNESSetConvergenceTest" 3274 /*@C 3275 SNESSetConvergenceTest - Sets the function that is to be used 3276 to test for convergence of the nonlinear iterative solution. 3277 3278 Logically Collective on SNES 3279 3280 Input Parameters: 3281 + snes - the SNES context 3282 . func - routine to test for convergence 3283 . cctx - [optional] context for private data for the convergence routine (may be PETSC_NULL) 3284 - destroy - [optional] destructor for the context (may be PETSC_NULL; PETSC_NULL_FUNCTION in Fortran) 3285 3286 Calling sequence of func: 3287 $ PetscErrorCode func (SNES snes,PetscInt it,PetscReal xnorm,PetscReal gnorm,PetscReal f,SNESConvergedReason *reason,void *cctx) 3288 3289 + snes - the SNES context 3290 . it - current iteration (0 is the first and is before any Newton step) 3291 . cctx - [optional] convergence context 3292 . reason - reason for convergence/divergence 3293 . xnorm - 2-norm of current iterate 3294 . gnorm - 2-norm of current step 3295 - f - 2-norm of function 3296 3297 Level: advanced 3298 3299 .keywords: SNES, nonlinear, set, convergence, test 3300 3301 .seealso: SNESDefaultConverged(), SNESSkipConverged() 3302 @*/ 3303 PetscErrorCode SNESSetConvergenceTest(SNES snes,PetscErrorCode (*func)(SNES,PetscInt,PetscReal,PetscReal,PetscReal,SNESConvergedReason*,void*),void *cctx,PetscErrorCode (*destroy)(void*)) 3304 { 3305 PetscErrorCode ierr; 3306 3307 PetscFunctionBegin; 3308 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3309 if (!func) func = SNESSkipConverged; 3310 if (snes->ops->convergeddestroy) { 3311 ierr = (*snes->ops->convergeddestroy)(snes->cnvP);CHKERRQ(ierr); 3312 } 3313 snes->ops->converged = func; 3314 snes->ops->convergeddestroy = destroy; 3315 snes->cnvP = cctx; 3316 PetscFunctionReturn(0); 3317 } 3318 3319 #undef __FUNCT__ 3320 #define __FUNCT__ "SNESGetConvergedReason" 3321 /*@ 3322 SNESGetConvergedReason - Gets the reason the SNES iteration was stopped. 3323 3324 Not Collective 3325 3326 Input Parameter: 3327 . snes - the SNES context 3328 3329 Output Parameter: 3330 . reason - negative value indicates diverged, positive value converged, see SNESConvergedReason or the 3331 manual pages for the individual convergence tests for complete lists 3332 3333 Level: intermediate 3334 3335 Notes: Can only be called after the call the SNESSolve() is complete. 3336 3337 .keywords: SNES, nonlinear, set, convergence, test 3338 3339 .seealso: SNESSetConvergenceTest(), SNESConvergedReason 3340 @*/ 3341 PetscErrorCode SNESGetConvergedReason(SNES snes,SNESConvergedReason *reason) 3342 { 3343 PetscFunctionBegin; 3344 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3345 PetscValidPointer(reason,2); 3346 *reason = snes->reason; 3347 PetscFunctionReturn(0); 3348 } 3349 3350 #undef __FUNCT__ 3351 #define __FUNCT__ "SNESSetConvergenceHistory" 3352 /*@ 3353 SNESSetConvergenceHistory - Sets the array used to hold the convergence history. 3354 3355 Logically Collective on SNES 3356 3357 Input Parameters: 3358 + snes - iterative context obtained from SNESCreate() 3359 . a - array to hold history, this array will contain the function norms computed at each step 3360 . its - integer array holds the number of linear iterations for each solve. 3361 . na - size of a and its 3362 - reset - PETSC_TRUE indicates each new nonlinear solve resets the history counter to zero, 3363 else it continues storing new values for new nonlinear solves after the old ones 3364 3365 Notes: 3366 If 'a' and 'its' are PETSC_NULL then space is allocated for the history. If 'na' PETSC_DECIDE or PETSC_DEFAULT then a 3367 default array of length 10000 is allocated. 3368 3369 This routine is useful, e.g., when running a code for purposes 3370 of accurate performance monitoring, when no I/O should be done 3371 during the section of code that is being timed. 3372 3373 Level: intermediate 3374 3375 .keywords: SNES, set, convergence, history 3376 3377 .seealso: SNESGetConvergenceHistory() 3378 3379 @*/ 3380 PetscErrorCode SNESSetConvergenceHistory(SNES snes,PetscReal a[],PetscInt its[],PetscInt na,PetscBool reset) 3381 { 3382 PetscErrorCode ierr; 3383 3384 PetscFunctionBegin; 3385 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3386 if (na) PetscValidScalarPointer(a,2); 3387 if (its) PetscValidIntPointer(its,3); 3388 if (na == PETSC_DECIDE || na == PETSC_DEFAULT || !a) { 3389 if (na == PETSC_DECIDE || na == PETSC_DEFAULT) na = 1000; 3390 ierr = PetscMalloc(na*sizeof(PetscReal),&a);CHKERRQ(ierr); 3391 ierr = PetscMalloc(na*sizeof(PetscInt),&its);CHKERRQ(ierr); 3392 snes->conv_malloc = PETSC_TRUE; 3393 } 3394 snes->conv_hist = a; 3395 snes->conv_hist_its = its; 3396 snes->conv_hist_max = na; 3397 snes->conv_hist_len = 0; 3398 snes->conv_hist_reset = reset; 3399 PetscFunctionReturn(0); 3400 } 3401 3402 #if defined(PETSC_HAVE_MATLAB_ENGINE) 3403 #include <engine.h> /* MATLAB include file */ 3404 #include <mex.h> /* MATLAB include file */ 3405 EXTERN_C_BEGIN 3406 #undef __FUNCT__ 3407 #define __FUNCT__ "SNESGetConvergenceHistoryMatlab" 3408 mxArray *SNESGetConvergenceHistoryMatlab(SNES snes) 3409 { 3410 mxArray *mat; 3411 PetscInt i; 3412 PetscReal *ar; 3413 3414 PetscFunctionBegin; 3415 mat = mxCreateDoubleMatrix(snes->conv_hist_len,1,mxREAL); 3416 ar = (PetscReal*) mxGetData(mat); 3417 for (i=0; i<snes->conv_hist_len; i++) { 3418 ar[i] = snes->conv_hist[i]; 3419 } 3420 PetscFunctionReturn(mat); 3421 } 3422 EXTERN_C_END 3423 #endif 3424 3425 3426 #undef __FUNCT__ 3427 #define __FUNCT__ "SNESGetConvergenceHistory" 3428 /*@C 3429 SNESGetConvergenceHistory - Gets the array used to hold the convergence history. 3430 3431 Not Collective 3432 3433 Input Parameter: 3434 . snes - iterative context obtained from SNESCreate() 3435 3436 Output Parameters: 3437 . a - array to hold history 3438 . its - integer array holds the number of linear iterations (or 3439 negative if not converged) for each solve. 3440 - na - size of a and its 3441 3442 Notes: 3443 The calling sequence for this routine in Fortran is 3444 $ call SNESGetConvergenceHistory(SNES snes, integer na, integer ierr) 3445 3446 This routine is useful, e.g., when running a code for purposes 3447 of accurate performance monitoring, when no I/O should be done 3448 during the section of code that is being timed. 3449 3450 Level: intermediate 3451 3452 .keywords: SNES, get, convergence, history 3453 3454 .seealso: SNESSetConvergencHistory() 3455 3456 @*/ 3457 PetscErrorCode SNESGetConvergenceHistory(SNES snes,PetscReal *a[],PetscInt *its[],PetscInt *na) 3458 { 3459 PetscFunctionBegin; 3460 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3461 if (a) *a = snes->conv_hist; 3462 if (its) *its = snes->conv_hist_its; 3463 if (na) *na = snes->conv_hist_len; 3464 PetscFunctionReturn(0); 3465 } 3466 3467 #undef __FUNCT__ 3468 #define __FUNCT__ "SNESSetUpdate" 3469 /*@C 3470 SNESSetUpdate - Sets the general-purpose update function called 3471 at the beginning of every iteration of the nonlinear solve. Specifically 3472 it is called just before the Jacobian is "evaluated". 3473 3474 Logically Collective on SNES 3475 3476 Input Parameters: 3477 . snes - The nonlinear solver context 3478 . func - The function 3479 3480 Calling sequence of func: 3481 . func (SNES snes, PetscInt step); 3482 3483 . step - The current step of the iteration 3484 3485 Level: advanced 3486 3487 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() 3488 This is not used by most users. 3489 3490 .keywords: SNES, update 3491 3492 .seealso SNESDefaultUpdate(), SNESSetJacobian(), SNESSolve() 3493 @*/ 3494 PetscErrorCode SNESSetUpdate(SNES snes, PetscErrorCode (*func)(SNES, PetscInt)) 3495 { 3496 PetscFunctionBegin; 3497 PetscValidHeaderSpecific(snes, SNES_CLASSID,1); 3498 snes->ops->update = func; 3499 PetscFunctionReturn(0); 3500 } 3501 3502 #undef __FUNCT__ 3503 #define __FUNCT__ "SNESDefaultUpdate" 3504 /*@ 3505 SNESDefaultUpdate - The default update function which does nothing. 3506 3507 Not collective 3508 3509 Input Parameters: 3510 . snes - The nonlinear solver context 3511 . step - The current step of the iteration 3512 3513 Level: intermediate 3514 3515 .keywords: SNES, update 3516 .seealso SNESSetUpdate(), SNESDefaultRhsBC(), SNESDefaultShortolutionBC() 3517 @*/ 3518 PetscErrorCode SNESDefaultUpdate(SNES snes, PetscInt step) 3519 { 3520 PetscFunctionBegin; 3521 PetscFunctionReturn(0); 3522 } 3523 3524 #undef __FUNCT__ 3525 #define __FUNCT__ "SNESScaleStep_Private" 3526 /* 3527 SNESScaleStep_Private - Scales a step so that its length is less than the 3528 positive parameter delta. 3529 3530 Input Parameters: 3531 + snes - the SNES context 3532 . y - approximate solution of linear system 3533 . fnorm - 2-norm of current function 3534 - delta - trust region size 3535 3536 Output Parameters: 3537 + gpnorm - predicted function norm at the new point, assuming local 3538 linearization. The value is zero if the step lies within the trust 3539 region, and exceeds zero otherwise. 3540 - ynorm - 2-norm of the step 3541 3542 Note: 3543 For non-trust region methods such as SNESLS, the parameter delta 3544 is set to be the maximum allowable step size. 3545 3546 .keywords: SNES, nonlinear, scale, step 3547 */ 3548 PetscErrorCode SNESScaleStep_Private(SNES snes,Vec y,PetscReal *fnorm,PetscReal *delta,PetscReal *gpnorm,PetscReal *ynorm) 3549 { 3550 PetscReal nrm; 3551 PetscScalar cnorm; 3552 PetscErrorCode ierr; 3553 3554 PetscFunctionBegin; 3555 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3556 PetscValidHeaderSpecific(y,VEC_CLASSID,2); 3557 PetscCheckSameComm(snes,1,y,2); 3558 3559 ierr = VecNorm(y,NORM_2,&nrm);CHKERRQ(ierr); 3560 if (nrm > *delta) { 3561 nrm = *delta/nrm; 3562 *gpnorm = (1.0 - nrm)*(*fnorm); 3563 cnorm = nrm; 3564 ierr = VecScale(y,cnorm);CHKERRQ(ierr); 3565 *ynorm = *delta; 3566 } else { 3567 *gpnorm = 0.0; 3568 *ynorm = nrm; 3569 } 3570 PetscFunctionReturn(0); 3571 } 3572 3573 #undef __FUNCT__ 3574 #define __FUNCT__ "SNESSolve" 3575 /*@C 3576 SNESSolve - Solves a nonlinear system F(x) = b. 3577 Call SNESSolve() after calling SNESCreate() and optional routines of the form SNESSetXXX(). 3578 3579 Collective on SNES 3580 3581 Input Parameters: 3582 + snes - the SNES context 3583 . b - the constant part of the equation F(x) = b, or PETSC_NULL to use zero. 3584 - x - the solution vector. 3585 3586 Notes: 3587 The user should initialize the vector,x, with the initial guess 3588 for the nonlinear solve prior to calling SNESSolve. In particular, 3589 to employ an initial guess of zero, the user should explicitly set 3590 this vector to zero by calling VecSet(). 3591 3592 Level: beginner 3593 3594 .keywords: SNES, nonlinear, solve 3595 3596 .seealso: SNESCreate(), SNESDestroy(), SNESSetFunction(), SNESSetJacobian(), SNESSetGridSequence(), SNESGetSolution() 3597 @*/ 3598 PetscErrorCode SNESSolve(SNES snes,Vec b,Vec x) 3599 { 3600 PetscErrorCode ierr; 3601 PetscBool flg; 3602 char filename[PETSC_MAX_PATH_LEN]; 3603 PetscViewer viewer; 3604 PetscInt grid; 3605 Vec xcreated = PETSC_NULL; 3606 DM dm; 3607 3608 PetscFunctionBegin; 3609 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3610 if (x) PetscValidHeaderSpecific(x,VEC_CLASSID,3); 3611 if (x) PetscCheckSameComm(snes,1,x,3); 3612 if (b) PetscValidHeaderSpecific(b,VEC_CLASSID,2); 3613 if (b) PetscCheckSameComm(snes,1,b,2); 3614 3615 if (!x) { 3616 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 3617 ierr = DMCreateGlobalVector(dm,&xcreated);CHKERRQ(ierr); 3618 x = xcreated; 3619 } 3620 3621 for (grid=0; grid<snes->gridsequence; grid++) {ierr = PetscViewerASCIIPushTab(PETSC_VIEWER_STDOUT_(((PetscObject)snes)->comm));CHKERRQ(ierr);} 3622 for (grid=0; grid<snes->gridsequence+1; grid++) { 3623 3624 /* set solution vector */ 3625 if (!grid) {ierr = PetscObjectReference((PetscObject)x);CHKERRQ(ierr);} 3626 ierr = VecDestroy(&snes->vec_sol);CHKERRQ(ierr); 3627 snes->vec_sol = x; 3628 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 3629 3630 /* set affine vector if provided */ 3631 if (b) { ierr = PetscObjectReference((PetscObject)b);CHKERRQ(ierr); } 3632 ierr = VecDestroy(&snes->vec_rhs);CHKERRQ(ierr); 3633 snes->vec_rhs = b; 3634 3635 ierr = SNESSetUp(snes);CHKERRQ(ierr); 3636 3637 if (!grid) { 3638 if (snes->ops->computeinitialguess) { 3639 ierr = (*snes->ops->computeinitialguess)(snes,snes->vec_sol,snes->initialguessP);CHKERRQ(ierr); 3640 } 3641 } 3642 3643 if (snes->conv_hist_reset) snes->conv_hist_len = 0; 3644 snes->nfuncs = 0; snes->linear_its = 0; snes->numFailures = 0; 3645 3646 ierr = PetscLogEventBegin(SNES_Solve,snes,0,0,0);CHKERRQ(ierr); 3647 ierr = (*snes->ops->solve)(snes);CHKERRQ(ierr); 3648 ierr = PetscLogEventEnd(SNES_Solve,snes,0,0,0);CHKERRQ(ierr); 3649 if (snes->domainerror){ 3650 snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; 3651 snes->domainerror = PETSC_FALSE; 3652 } 3653 if (!snes->reason) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Internal error, solver returned without setting converged reason"); 3654 3655 flg = PETSC_FALSE; 3656 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_test_local_min",&flg,PETSC_NULL);CHKERRQ(ierr); 3657 if (flg && !PetscPreLoadingOn) { ierr = SNESTestLocalMin(snes);CHKERRQ(ierr); } 3658 if (snes->printreason) { 3659 ierr = PetscViewerASCIIAddTab(PETSC_VIEWER_STDOUT_(((PetscObject)snes)->comm),((PetscObject)snes)->tablevel);CHKERRQ(ierr); 3660 if (snes->reason > 0) { 3661 ierr = PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_(((PetscObject)snes)->comm),"Nonlinear solve converged due to %s iterations %D\n",SNESConvergedReasons[snes->reason],snes->iter);CHKERRQ(ierr); 3662 } else { 3663 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); 3664 } 3665 ierr = PetscViewerASCIISubtractTab(PETSC_VIEWER_STDOUT_(((PetscObject)snes)->comm),((PetscObject)snes)->tablevel);CHKERRQ(ierr); 3666 } 3667 3668 if (snes->errorifnotconverged && snes->reason < 0) SETERRQ(((PetscObject)snes)->comm,PETSC_ERR_NOT_CONVERGED,"SNESSolve has not converged"); 3669 if (grid < snes->gridsequence) { 3670 DM fine; 3671 Vec xnew; 3672 Mat interp; 3673 3674 ierr = DMRefine(snes->dm,((PetscObject)snes)->comm,&fine);CHKERRQ(ierr); 3675 if (!fine) SETERRQ(((PetscObject)snes)->comm,PETSC_ERR_ARG_INCOMP,"DMRefine() did not perform any refinement, cannot continue grid sequencing"); 3676 ierr = DMCreateInterpolation(snes->dm,fine,&interp,PETSC_NULL);CHKERRQ(ierr); 3677 ierr = DMCreateGlobalVector(fine,&xnew);CHKERRQ(ierr); 3678 ierr = MatInterpolate(interp,x,xnew);CHKERRQ(ierr); 3679 ierr = DMInterpolate(snes->dm,interp,fine);CHKERRQ(ierr); 3680 ierr = MatDestroy(&interp);CHKERRQ(ierr); 3681 x = xnew; 3682 3683 ierr = SNESReset(snes);CHKERRQ(ierr); 3684 ierr = SNESSetDM(snes,fine);CHKERRQ(ierr); 3685 ierr = DMDestroy(&fine);CHKERRQ(ierr); 3686 ierr = PetscViewerASCIIPopTab(PETSC_VIEWER_STDOUT_(((PetscObject)snes)->comm));CHKERRQ(ierr); 3687 } 3688 } 3689 /* monitoring and viewing */ 3690 flg = PETSC_FALSE; 3691 ierr = PetscOptionsGetString(((PetscObject)snes)->prefix,"-snes_view",filename,PETSC_MAX_PATH_LEN,&flg);CHKERRQ(ierr); 3692 if (flg && !PetscPreLoadingOn) { 3693 ierr = PetscViewerASCIIOpen(((PetscObject)snes)->comm,filename,&viewer);CHKERRQ(ierr); 3694 ierr = SNESView(snes,viewer);CHKERRQ(ierr); 3695 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 3696 } 3697 3698 flg = PETSC_FALSE; 3699 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_view_draw",&flg,PETSC_NULL);CHKERRQ(ierr); 3700 if (flg) { 3701 ierr = PetscViewerDrawOpen(((PetscObject)snes)->comm,PETSC_NULL,"SNES Solver",0,0,600,600,&viewer);CHKERRQ(ierr); 3702 ierr = SNESView(snes,viewer);CHKERRQ(ierr); 3703 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 3704 } 3705 3706 flg = PETSC_FALSE; 3707 ierr = PetscOptionsGetString(((PetscObject)snes)->prefix,"-snes_view_solution_vtk",filename,PETSC_MAX_PATH_LEN,&flg);CHKERRQ(ierr); 3708 if (flg) { 3709 PetscViewer viewer; 3710 ierr = PetscViewerCreate(((PetscObject)snes)->comm,&viewer);CHKERRQ(ierr); 3711 ierr = PetscViewerSetType(viewer,PETSCVIEWERVTK);CHKERRQ(ierr); 3712 ierr = PetscViewerFileSetName(viewer,filename);CHKERRQ(ierr); 3713 ierr = VecView(snes->vec_sol,viewer);CHKERRQ(ierr); 3714 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 3715 } 3716 3717 ierr = VecDestroy(&xcreated);CHKERRQ(ierr); 3718 PetscFunctionReturn(0); 3719 } 3720 3721 /* --------- Internal routines for SNES Package --------- */ 3722 3723 #undef __FUNCT__ 3724 #define __FUNCT__ "SNESSetType" 3725 /*@C 3726 SNESSetType - Sets the method for the nonlinear solver. 3727 3728 Collective on SNES 3729 3730 Input Parameters: 3731 + snes - the SNES context 3732 - type - a known method 3733 3734 Options Database Key: 3735 . -snes_type <type> - Sets the method; use -help for a list 3736 of available methods (for instance, ls or tr) 3737 3738 Notes: 3739 See "petsc/include/petscsnes.h" for available methods (for instance) 3740 + SNESLS - Newton's method with line search 3741 (systems of nonlinear equations) 3742 . SNESTR - Newton's method with trust region 3743 (systems of nonlinear equations) 3744 3745 Normally, it is best to use the SNESSetFromOptions() command and then 3746 set the SNES solver type from the options database rather than by using 3747 this routine. Using the options database provides the user with 3748 maximum flexibility in evaluating the many nonlinear solvers. 3749 The SNESSetType() routine is provided for those situations where it 3750 is necessary to set the nonlinear solver independently of the command 3751 line or options database. This might be the case, for example, when 3752 the choice of solver changes during the execution of the program, 3753 and the user's application is taking responsibility for choosing the 3754 appropriate method. 3755 3756 Level: intermediate 3757 3758 .keywords: SNES, set, type 3759 3760 .seealso: SNESType, SNESCreate() 3761 3762 @*/ 3763 PetscErrorCode SNESSetType(SNES snes,SNESType type) 3764 { 3765 PetscErrorCode ierr,(*r)(SNES); 3766 PetscBool match; 3767 3768 PetscFunctionBegin; 3769 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3770 PetscValidCharPointer(type,2); 3771 3772 ierr = PetscObjectTypeCompare((PetscObject)snes,type,&match);CHKERRQ(ierr); 3773 if (match) PetscFunctionReturn(0); 3774 3775 ierr = PetscFListFind(SNESList,((PetscObject)snes)->comm,type,PETSC_TRUE,(void (**)(void)) &r);CHKERRQ(ierr); 3776 if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unable to find requested SNES type %s",type); 3777 /* Destroy the previous private SNES context */ 3778 if (snes->ops->destroy) { 3779 ierr = (*(snes)->ops->destroy)(snes);CHKERRQ(ierr); 3780 snes->ops->destroy = PETSC_NULL; 3781 } 3782 /* Reinitialize function pointers in SNESOps structure */ 3783 snes->ops->setup = 0; 3784 snes->ops->solve = 0; 3785 snes->ops->view = 0; 3786 snes->ops->setfromoptions = 0; 3787 snes->ops->destroy = 0; 3788 /* Call the SNESCreate_XXX routine for this particular Nonlinear solver */ 3789 snes->setupcalled = PETSC_FALSE; 3790 ierr = PetscObjectChangeTypeName((PetscObject)snes,type);CHKERRQ(ierr); 3791 ierr = (*r)(snes);CHKERRQ(ierr); 3792 #if defined(PETSC_HAVE_AMS) 3793 if (PetscAMSPublishAll) { 3794 ierr = PetscObjectAMSPublish((PetscObject)snes);CHKERRQ(ierr); 3795 } 3796 #endif 3797 PetscFunctionReturn(0); 3798 } 3799 3800 3801 /* --------------------------------------------------------------------- */ 3802 #undef __FUNCT__ 3803 #define __FUNCT__ "SNESRegisterDestroy" 3804 /*@ 3805 SNESRegisterDestroy - Frees the list of nonlinear solvers that were 3806 registered by SNESRegisterDynamic(). 3807 3808 Not Collective 3809 3810 Level: advanced 3811 3812 .keywords: SNES, nonlinear, register, destroy 3813 3814 .seealso: SNESRegisterAll(), SNESRegisterAll() 3815 @*/ 3816 PetscErrorCode SNESRegisterDestroy(void) 3817 { 3818 PetscErrorCode ierr; 3819 3820 PetscFunctionBegin; 3821 ierr = PetscFListDestroy(&SNESList);CHKERRQ(ierr); 3822 SNESRegisterAllCalled = PETSC_FALSE; 3823 PetscFunctionReturn(0); 3824 } 3825 3826 #undef __FUNCT__ 3827 #define __FUNCT__ "SNESGetType" 3828 /*@C 3829 SNESGetType - Gets the SNES method type and name (as a string). 3830 3831 Not Collective 3832 3833 Input Parameter: 3834 . snes - nonlinear solver context 3835 3836 Output Parameter: 3837 . type - SNES method (a character string) 3838 3839 Level: intermediate 3840 3841 .keywords: SNES, nonlinear, get, type, name 3842 @*/ 3843 PetscErrorCode SNESGetType(SNES snes,SNESType *type) 3844 { 3845 PetscFunctionBegin; 3846 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3847 PetscValidPointer(type,2); 3848 *type = ((PetscObject)snes)->type_name; 3849 PetscFunctionReturn(0); 3850 } 3851 3852 #undef __FUNCT__ 3853 #define __FUNCT__ "SNESGetSolution" 3854 /*@ 3855 SNESGetSolution - Returns the vector where the approximate solution is 3856 stored. This is the fine grid solution when using SNESSetGridSequence(). 3857 3858 Not Collective, but Vec is parallel if SNES is parallel 3859 3860 Input Parameter: 3861 . snes - the SNES context 3862 3863 Output Parameter: 3864 . x - the solution 3865 3866 Level: intermediate 3867 3868 .keywords: SNES, nonlinear, get, solution 3869 3870 .seealso: SNESGetSolutionUpdate(), SNESGetFunction() 3871 @*/ 3872 PetscErrorCode SNESGetSolution(SNES snes,Vec *x) 3873 { 3874 PetscFunctionBegin; 3875 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3876 PetscValidPointer(x,2); 3877 *x = snes->vec_sol; 3878 PetscFunctionReturn(0); 3879 } 3880 3881 #undef __FUNCT__ 3882 #define __FUNCT__ "SNESGetSolutionUpdate" 3883 /*@ 3884 SNESGetSolutionUpdate - Returns the vector where the solution update is 3885 stored. 3886 3887 Not Collective, but Vec is parallel if SNES is parallel 3888 3889 Input Parameter: 3890 . snes - the SNES context 3891 3892 Output Parameter: 3893 . x - the solution update 3894 3895 Level: advanced 3896 3897 .keywords: SNES, nonlinear, get, solution, update 3898 3899 .seealso: SNESGetSolution(), SNESGetFunction() 3900 @*/ 3901 PetscErrorCode SNESGetSolutionUpdate(SNES snes,Vec *x) 3902 { 3903 PetscFunctionBegin; 3904 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3905 PetscValidPointer(x,2); 3906 *x = snes->vec_sol_update; 3907 PetscFunctionReturn(0); 3908 } 3909 3910 #undef __FUNCT__ 3911 #define __FUNCT__ "SNESGetFunction" 3912 /*@C 3913 SNESGetFunction - Returns the vector where the function is stored. 3914 3915 Not Collective, but Vec is parallel if SNES is parallel. Collective if Vec is requested, but has not been created yet. 3916 3917 Input Parameter: 3918 . snes - the SNES context 3919 3920 Output Parameter: 3921 + r - the function (or PETSC_NULL) 3922 . func - the function (or PETSC_NULL) 3923 - ctx - the function context (or PETSC_NULL) 3924 3925 Level: advanced 3926 3927 .keywords: SNES, nonlinear, get, function 3928 3929 .seealso: SNESSetFunction(), SNESGetSolution() 3930 @*/ 3931 PetscErrorCode SNESGetFunction(SNES snes,Vec *r,PetscErrorCode (**func)(SNES,Vec,Vec,void*),void **ctx) 3932 { 3933 PetscErrorCode ierr; 3934 DM dm; 3935 3936 PetscFunctionBegin; 3937 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3938 if (r) { 3939 if (!snes->vec_func) { 3940 if (snes->vec_rhs) { 3941 ierr = VecDuplicate(snes->vec_rhs,&snes->vec_func);CHKERRQ(ierr); 3942 } else if (snes->vec_sol) { 3943 ierr = VecDuplicate(snes->vec_sol,&snes->vec_func);CHKERRQ(ierr); 3944 } else if (snes->dm) { 3945 ierr = DMCreateGlobalVector(snes->dm,&snes->vec_func);CHKERRQ(ierr); 3946 } 3947 } 3948 *r = snes->vec_func; 3949 } 3950 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 3951 ierr = DMSNESGetFunction(dm,func,ctx);CHKERRQ(ierr); 3952 PetscFunctionReturn(0); 3953 } 3954 3955 /*@C 3956 SNESGetGS - Returns the GS function and context. 3957 3958 Input Parameter: 3959 . snes - the SNES context 3960 3961 Output Parameter: 3962 + gsfunc - the function (or PETSC_NULL) 3963 - ctx - the function context (or PETSC_NULL) 3964 3965 Level: advanced 3966 3967 .keywords: SNES, nonlinear, get, function 3968 3969 .seealso: SNESSetGS(), SNESGetFunction() 3970 @*/ 3971 3972 #undef __FUNCT__ 3973 #define __FUNCT__ "SNESGetGS" 3974 PetscErrorCode SNESGetGS (SNES snes, PetscErrorCode(**func)(SNES, Vec, Vec, void*), void ** ctx) 3975 { 3976 PetscErrorCode ierr; 3977 DM dm; 3978 3979 PetscFunctionBegin; 3980 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3981 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 3982 ierr = DMSNESGetGS(dm,func,ctx);CHKERRQ(ierr); 3983 PetscFunctionReturn(0); 3984 } 3985 3986 #undef __FUNCT__ 3987 #define __FUNCT__ "SNESSetOptionsPrefix" 3988 /*@C 3989 SNESSetOptionsPrefix - Sets the prefix used for searching for all 3990 SNES options in the database. 3991 3992 Logically Collective on SNES 3993 3994 Input Parameter: 3995 + snes - the SNES context 3996 - prefix - the prefix to prepend to all option names 3997 3998 Notes: 3999 A hyphen (-) must NOT be given at the beginning of the prefix name. 4000 The first character of all runtime options is AUTOMATICALLY the hyphen. 4001 4002 Level: advanced 4003 4004 .keywords: SNES, set, options, prefix, database 4005 4006 .seealso: SNESSetFromOptions() 4007 @*/ 4008 PetscErrorCode SNESSetOptionsPrefix(SNES snes,const char prefix[]) 4009 { 4010 PetscErrorCode ierr; 4011 4012 PetscFunctionBegin; 4013 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4014 ierr = PetscObjectSetOptionsPrefix((PetscObject)snes,prefix);CHKERRQ(ierr); 4015 if (!snes->ksp) {ierr = SNESGetKSP(snes,&snes->ksp);CHKERRQ(ierr);} 4016 if (snes->linesearch) { 4017 ierr = SNESGetSNESLineSearch(snes,&snes->linesearch);CHKERRQ(ierr); 4018 ierr = PetscObjectSetOptionsPrefix((PetscObject)snes->linesearch,prefix);CHKERRQ(ierr); 4019 } 4020 ierr = KSPSetOptionsPrefix(snes->ksp,prefix);CHKERRQ(ierr); 4021 PetscFunctionReturn(0); 4022 } 4023 4024 #undef __FUNCT__ 4025 #define __FUNCT__ "SNESAppendOptionsPrefix" 4026 /*@C 4027 SNESAppendOptionsPrefix - Appends to the prefix used for searching for all 4028 SNES options in the database. 4029 4030 Logically Collective on SNES 4031 4032 Input Parameters: 4033 + snes - the SNES context 4034 - prefix - the prefix to prepend to all option names 4035 4036 Notes: 4037 A hyphen (-) must NOT be given at the beginning of the prefix name. 4038 The first character of all runtime options is AUTOMATICALLY the hyphen. 4039 4040 Level: advanced 4041 4042 .keywords: SNES, append, options, prefix, database 4043 4044 .seealso: SNESGetOptionsPrefix() 4045 @*/ 4046 PetscErrorCode SNESAppendOptionsPrefix(SNES snes,const char prefix[]) 4047 { 4048 PetscErrorCode ierr; 4049 4050 PetscFunctionBegin; 4051 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4052 ierr = PetscObjectAppendOptionsPrefix((PetscObject)snes,prefix);CHKERRQ(ierr); 4053 if (!snes->ksp) {ierr = SNESGetKSP(snes,&snes->ksp);CHKERRQ(ierr);} 4054 if (snes->linesearch) { 4055 ierr = SNESGetSNESLineSearch(snes,&snes->linesearch);CHKERRQ(ierr); 4056 ierr = PetscObjectAppendOptionsPrefix((PetscObject)snes->linesearch,prefix);CHKERRQ(ierr); 4057 } 4058 ierr = KSPAppendOptionsPrefix(snes->ksp,prefix);CHKERRQ(ierr); 4059 PetscFunctionReturn(0); 4060 } 4061 4062 #undef __FUNCT__ 4063 #define __FUNCT__ "SNESGetOptionsPrefix" 4064 /*@C 4065 SNESGetOptionsPrefix - Sets the prefix used for searching for all 4066 SNES options in the database. 4067 4068 Not Collective 4069 4070 Input Parameter: 4071 . snes - the SNES context 4072 4073 Output Parameter: 4074 . prefix - pointer to the prefix string used 4075 4076 Notes: On the fortran side, the user should pass in a string 'prefix' of 4077 sufficient length to hold the prefix. 4078 4079 Level: advanced 4080 4081 .keywords: SNES, get, options, prefix, database 4082 4083 .seealso: SNESAppendOptionsPrefix() 4084 @*/ 4085 PetscErrorCode SNESGetOptionsPrefix(SNES snes,const char *prefix[]) 4086 { 4087 PetscErrorCode ierr; 4088 4089 PetscFunctionBegin; 4090 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4091 ierr = PetscObjectGetOptionsPrefix((PetscObject)snes,prefix);CHKERRQ(ierr); 4092 PetscFunctionReturn(0); 4093 } 4094 4095 4096 #undef __FUNCT__ 4097 #define __FUNCT__ "SNESRegister" 4098 /*@C 4099 SNESRegister - See SNESRegisterDynamic() 4100 4101 Level: advanced 4102 @*/ 4103 PetscErrorCode SNESRegister(const char sname[],const char path[],const char name[],PetscErrorCode (*function)(SNES)) 4104 { 4105 char fullname[PETSC_MAX_PATH_LEN]; 4106 PetscErrorCode ierr; 4107 4108 PetscFunctionBegin; 4109 ierr = PetscFListConcat(path,name,fullname);CHKERRQ(ierr); 4110 ierr = PetscFListAdd(&SNESList,sname,fullname,(void (*)(void))function);CHKERRQ(ierr); 4111 PetscFunctionReturn(0); 4112 } 4113 4114 #undef __FUNCT__ 4115 #define __FUNCT__ "SNESTestLocalMin" 4116 PetscErrorCode SNESTestLocalMin(SNES snes) 4117 { 4118 PetscErrorCode ierr; 4119 PetscInt N,i,j; 4120 Vec u,uh,fh; 4121 PetscScalar value; 4122 PetscReal norm; 4123 4124 PetscFunctionBegin; 4125 ierr = SNESGetSolution(snes,&u);CHKERRQ(ierr); 4126 ierr = VecDuplicate(u,&uh);CHKERRQ(ierr); 4127 ierr = VecDuplicate(u,&fh);CHKERRQ(ierr); 4128 4129 /* currently only works for sequential */ 4130 ierr = PetscPrintf(PETSC_COMM_WORLD,"Testing FormFunction() for local min\n"); 4131 ierr = VecGetSize(u,&N);CHKERRQ(ierr); 4132 for (i=0; i<N; i++) { 4133 ierr = VecCopy(u,uh);CHKERRQ(ierr); 4134 ierr = PetscPrintf(PETSC_COMM_WORLD,"i = %D\n",i);CHKERRQ(ierr); 4135 for (j=-10; j<11; j++) { 4136 value = PetscSign(j)*exp(PetscAbs(j)-10.0); 4137 ierr = VecSetValue(uh,i,value,ADD_VALUES);CHKERRQ(ierr); 4138 ierr = SNESComputeFunction(snes,uh,fh);CHKERRQ(ierr); 4139 ierr = VecNorm(fh,NORM_2,&norm);CHKERRQ(ierr); 4140 ierr = PetscPrintf(PETSC_COMM_WORLD," j norm %D %18.16e\n",j,norm);CHKERRQ(ierr); 4141 value = -value; 4142 ierr = VecSetValue(uh,i,value,ADD_VALUES);CHKERRQ(ierr); 4143 } 4144 } 4145 ierr = VecDestroy(&uh);CHKERRQ(ierr); 4146 ierr = VecDestroy(&fh);CHKERRQ(ierr); 4147 PetscFunctionReturn(0); 4148 } 4149 4150 #undef __FUNCT__ 4151 #define __FUNCT__ "SNESKSPSetUseEW" 4152 /*@ 4153 SNESKSPSetUseEW - Sets SNES use Eisenstat-Walker method for 4154 computing relative tolerance for linear solvers within an inexact 4155 Newton method. 4156 4157 Logically Collective on SNES 4158 4159 Input Parameters: 4160 + snes - SNES context 4161 - flag - PETSC_TRUE or PETSC_FALSE 4162 4163 Options Database: 4164 + -snes_ksp_ew - use Eisenstat-Walker method for determining linear system convergence 4165 . -snes_ksp_ew_version ver - version of Eisenstat-Walker method 4166 . -snes_ksp_ew_rtol0 <rtol0> - Sets rtol0 4167 . -snes_ksp_ew_rtolmax <rtolmax> - Sets rtolmax 4168 . -snes_ksp_ew_gamma <gamma> - Sets gamma 4169 . -snes_ksp_ew_alpha <alpha> - Sets alpha 4170 . -snes_ksp_ew_alpha2 <alpha2> - Sets alpha2 4171 - -snes_ksp_ew_threshold <threshold> - Sets threshold 4172 4173 Notes: 4174 Currently, the default is to use a constant relative tolerance for 4175 the inner linear solvers. Alternatively, one can use the 4176 Eisenstat-Walker method, where the relative convergence tolerance 4177 is reset at each Newton iteration according progress of the nonlinear 4178 solver. 4179 4180 Level: advanced 4181 4182 Reference: 4183 S. C. Eisenstat and H. F. Walker, "Choosing the forcing terms in an 4184 inexact Newton method", SISC 17 (1), pp.16-32, 1996. 4185 4186 .keywords: SNES, KSP, Eisenstat, Walker, convergence, test, inexact, Newton 4187 4188 .seealso: SNESKSPGetUseEW(), SNESKSPGetParametersEW(), SNESKSPSetParametersEW() 4189 @*/ 4190 PetscErrorCode SNESKSPSetUseEW(SNES snes,PetscBool flag) 4191 { 4192 PetscFunctionBegin; 4193 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4194 PetscValidLogicalCollectiveBool(snes,flag,2); 4195 snes->ksp_ewconv = flag; 4196 PetscFunctionReturn(0); 4197 } 4198 4199 #undef __FUNCT__ 4200 #define __FUNCT__ "SNESKSPGetUseEW" 4201 /*@ 4202 SNESKSPGetUseEW - Gets if SNES is using Eisenstat-Walker method 4203 for computing relative tolerance for linear solvers within an 4204 inexact Newton method. 4205 4206 Not Collective 4207 4208 Input Parameter: 4209 . snes - SNES context 4210 4211 Output Parameter: 4212 . flag - PETSC_TRUE or PETSC_FALSE 4213 4214 Notes: 4215 Currently, the default is to use a constant relative tolerance for 4216 the inner linear solvers. Alternatively, one can use the 4217 Eisenstat-Walker method, where the relative convergence tolerance 4218 is reset at each Newton iteration according progress of the nonlinear 4219 solver. 4220 4221 Level: advanced 4222 4223 Reference: 4224 S. C. Eisenstat and H. F. Walker, "Choosing the forcing terms in an 4225 inexact Newton method", SISC 17 (1), pp.16-32, 1996. 4226 4227 .keywords: SNES, KSP, Eisenstat, Walker, convergence, test, inexact, Newton 4228 4229 .seealso: SNESKSPSetUseEW(), SNESKSPGetParametersEW(), SNESKSPSetParametersEW() 4230 @*/ 4231 PetscErrorCode SNESKSPGetUseEW(SNES snes, PetscBool *flag) 4232 { 4233 PetscFunctionBegin; 4234 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4235 PetscValidPointer(flag,2); 4236 *flag = snes->ksp_ewconv; 4237 PetscFunctionReturn(0); 4238 } 4239 4240 #undef __FUNCT__ 4241 #define __FUNCT__ "SNESKSPSetParametersEW" 4242 /*@ 4243 SNESKSPSetParametersEW - Sets parameters for Eisenstat-Walker 4244 convergence criteria for the linear solvers within an inexact 4245 Newton method. 4246 4247 Logically Collective on SNES 4248 4249 Input Parameters: 4250 + snes - SNES context 4251 . version - version 1, 2 (default is 2) or 3 4252 . rtol_0 - initial relative tolerance (0 <= rtol_0 < 1) 4253 . rtol_max - maximum relative tolerance (0 <= rtol_max < 1) 4254 . gamma - multiplicative factor for version 2 rtol computation 4255 (0 <= gamma2 <= 1) 4256 . alpha - power for version 2 rtol computation (1 < alpha <= 2) 4257 . alpha2 - power for safeguard 4258 - threshold - threshold for imposing safeguard (0 < threshold < 1) 4259 4260 Note: 4261 Version 3 was contributed by Luis Chacon, June 2006. 4262 4263 Use PETSC_DEFAULT to retain the default for any of the parameters. 4264 4265 Level: advanced 4266 4267 Reference: 4268 S. C. Eisenstat and H. F. Walker, "Choosing the forcing terms in an 4269 inexact Newton method", Utah State University Math. Stat. Dept. Res. 4270 Report 6/94/75, June, 1994, to appear in SIAM J. Sci. Comput. 4271 4272 .keywords: SNES, KSP, Eisenstat, Walker, set, parameters 4273 4274 .seealso: SNESKSPSetUseEW(), SNESKSPGetUseEW(), SNESKSPGetParametersEW() 4275 @*/ 4276 PetscErrorCode SNESKSPSetParametersEW(SNES snes,PetscInt version,PetscReal rtol_0,PetscReal rtol_max, 4277 PetscReal gamma,PetscReal alpha,PetscReal alpha2,PetscReal threshold) 4278 { 4279 SNESKSPEW *kctx; 4280 PetscFunctionBegin; 4281 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4282 kctx = (SNESKSPEW*)snes->kspconvctx; 4283 if (!kctx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"No Eisenstat-Walker context existing"); 4284 PetscValidLogicalCollectiveInt(snes,version,2); 4285 PetscValidLogicalCollectiveReal(snes,rtol_0,3); 4286 PetscValidLogicalCollectiveReal(snes,rtol_max,4); 4287 PetscValidLogicalCollectiveReal(snes,gamma,5); 4288 PetscValidLogicalCollectiveReal(snes,alpha,6); 4289 PetscValidLogicalCollectiveReal(snes,alpha2,7); 4290 PetscValidLogicalCollectiveReal(snes,threshold,8); 4291 4292 if (version != PETSC_DEFAULT) kctx->version = version; 4293 if (rtol_0 != PETSC_DEFAULT) kctx->rtol_0 = rtol_0; 4294 if (rtol_max != PETSC_DEFAULT) kctx->rtol_max = rtol_max; 4295 if (gamma != PETSC_DEFAULT) kctx->gamma = gamma; 4296 if (alpha != PETSC_DEFAULT) kctx->alpha = alpha; 4297 if (alpha2 != PETSC_DEFAULT) kctx->alpha2 = alpha2; 4298 if (threshold != PETSC_DEFAULT) kctx->threshold = threshold; 4299 4300 if (kctx->version < 1 || kctx->version > 3) { 4301 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Only versions 1, 2 and 3 are supported: %D",kctx->version); 4302 } 4303 if (kctx->rtol_0 < 0.0 || kctx->rtol_0 >= 1.0) { 4304 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"0.0 <= rtol_0 < 1.0: %G",kctx->rtol_0); 4305 } 4306 if (kctx->rtol_max < 0.0 || kctx->rtol_max >= 1.0) { 4307 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"0.0 <= rtol_max (%G) < 1.0\n",kctx->rtol_max); 4308 } 4309 if (kctx->gamma < 0.0 || kctx->gamma > 1.0) { 4310 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"0.0 <= gamma (%G) <= 1.0\n",kctx->gamma); 4311 } 4312 if (kctx->alpha <= 1.0 || kctx->alpha > 2.0) { 4313 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"1.0 < alpha (%G) <= 2.0\n",kctx->alpha); 4314 } 4315 if (kctx->threshold <= 0.0 || kctx->threshold >= 1.0) { 4316 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"0.0 < threshold (%G) < 1.0\n",kctx->threshold); 4317 } 4318 PetscFunctionReturn(0); 4319 } 4320 4321 #undef __FUNCT__ 4322 #define __FUNCT__ "SNESKSPGetParametersEW" 4323 /*@ 4324 SNESKSPGetParametersEW - Gets parameters for Eisenstat-Walker 4325 convergence criteria for the linear solvers within an inexact 4326 Newton method. 4327 4328 Not Collective 4329 4330 Input Parameters: 4331 snes - SNES context 4332 4333 Output Parameters: 4334 + version - version 1, 2 (default is 2) or 3 4335 . rtol_0 - initial relative tolerance (0 <= rtol_0 < 1) 4336 . rtol_max - maximum relative tolerance (0 <= rtol_max < 1) 4337 . gamma - multiplicative factor for version 2 rtol computation 4338 (0 <= gamma2 <= 1) 4339 . alpha - power for version 2 rtol computation (1 < alpha <= 2) 4340 . alpha2 - power for safeguard 4341 - threshold - threshold for imposing safeguard (0 < threshold < 1) 4342 4343 Level: advanced 4344 4345 .keywords: SNES, KSP, Eisenstat, Walker, get, parameters 4346 4347 .seealso: SNESKSPSetUseEW(), SNESKSPGetUseEW(), SNESKSPSetParametersEW() 4348 @*/ 4349 PetscErrorCode SNESKSPGetParametersEW(SNES snes,PetscInt *version,PetscReal *rtol_0,PetscReal *rtol_max, 4350 PetscReal *gamma,PetscReal *alpha,PetscReal *alpha2,PetscReal *threshold) 4351 { 4352 SNESKSPEW *kctx; 4353 PetscFunctionBegin; 4354 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4355 kctx = (SNESKSPEW*)snes->kspconvctx; 4356 if (!kctx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"No Eisenstat-Walker context existing"); 4357 if (version) *version = kctx->version; 4358 if (rtol_0) *rtol_0 = kctx->rtol_0; 4359 if (rtol_max) *rtol_max = kctx->rtol_max; 4360 if (gamma) *gamma = kctx->gamma; 4361 if (alpha) *alpha = kctx->alpha; 4362 if (alpha2) *alpha2 = kctx->alpha2; 4363 if (threshold) *threshold = kctx->threshold; 4364 PetscFunctionReturn(0); 4365 } 4366 4367 #undef __FUNCT__ 4368 #define __FUNCT__ "SNESKSPEW_PreSolve" 4369 static PetscErrorCode SNESKSPEW_PreSolve(SNES snes, KSP ksp, Vec b, Vec x) 4370 { 4371 PetscErrorCode ierr; 4372 SNESKSPEW *kctx = (SNESKSPEW*)snes->kspconvctx; 4373 PetscReal rtol=PETSC_DEFAULT,stol; 4374 4375 PetscFunctionBegin; 4376 if (!kctx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"No Eisenstat-Walker context exists"); 4377 if (!snes->iter) { /* first time in, so use the original user rtol */ 4378 rtol = kctx->rtol_0; 4379 } else { 4380 if (kctx->version == 1) { 4381 rtol = (snes->norm - kctx->lresid_last)/kctx->norm_last; 4382 if (rtol < 0.0) rtol = -rtol; 4383 stol = pow(kctx->rtol_last,kctx->alpha2); 4384 if (stol > kctx->threshold) rtol = PetscMax(rtol,stol); 4385 } else if (kctx->version == 2) { 4386 rtol = kctx->gamma * pow(snes->norm/kctx->norm_last,kctx->alpha); 4387 stol = kctx->gamma * pow(kctx->rtol_last,kctx->alpha); 4388 if (stol > kctx->threshold) rtol = PetscMax(rtol,stol); 4389 } else if (kctx->version == 3) {/* contributed by Luis Chacon, June 2006. */ 4390 rtol = kctx->gamma * pow(snes->norm/kctx->norm_last,kctx->alpha); 4391 /* safeguard: avoid sharp decrease of rtol */ 4392 stol = kctx->gamma*pow(kctx->rtol_last,kctx->alpha); 4393 stol = PetscMax(rtol,stol); 4394 rtol = PetscMin(kctx->rtol_0,stol); 4395 /* safeguard: avoid oversolving */ 4396 stol = kctx->gamma*(snes->ttol)/snes->norm; 4397 stol = PetscMax(rtol,stol); 4398 rtol = PetscMin(kctx->rtol_0,stol); 4399 } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Only versions 1, 2 or 3 are supported: %D",kctx->version); 4400 } 4401 /* safeguard: avoid rtol greater than one */ 4402 rtol = PetscMin(rtol,kctx->rtol_max); 4403 ierr = KSPSetTolerances(ksp,rtol,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);CHKERRQ(ierr); 4404 ierr = PetscInfo3(snes,"iter %D, Eisenstat-Walker (version %D) KSP rtol=%G\n",snes->iter,kctx->version,rtol);CHKERRQ(ierr); 4405 PetscFunctionReturn(0); 4406 } 4407 4408 #undef __FUNCT__ 4409 #define __FUNCT__ "SNESKSPEW_PostSolve" 4410 static PetscErrorCode SNESKSPEW_PostSolve(SNES snes, KSP ksp, Vec b, Vec x) 4411 { 4412 PetscErrorCode ierr; 4413 SNESKSPEW *kctx = (SNESKSPEW*)snes->kspconvctx; 4414 PCSide pcside; 4415 Vec lres; 4416 4417 PetscFunctionBegin; 4418 if (!kctx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"No Eisenstat-Walker context exists"); 4419 ierr = KSPGetTolerances(ksp,&kctx->rtol_last,0,0,0);CHKERRQ(ierr); 4420 ierr = SNESGetFunctionNorm(snes,&kctx->norm_last);CHKERRQ(ierr); 4421 if (kctx->version == 1) { 4422 ierr = KSPGetPCSide(ksp,&pcside);CHKERRQ(ierr); 4423 if (pcside == PC_RIGHT) { /* XXX Should we also test KSP_UNPRECONDITIONED_NORM ? */ 4424 /* KSP residual is true linear residual */ 4425 ierr = KSPGetResidualNorm(ksp,&kctx->lresid_last);CHKERRQ(ierr); 4426 } else { 4427 /* KSP residual is preconditioned residual */ 4428 /* compute true linear residual norm */ 4429 ierr = VecDuplicate(b,&lres);CHKERRQ(ierr); 4430 ierr = MatMult(snes->jacobian,x,lres);CHKERRQ(ierr); 4431 ierr = VecAYPX(lres,-1.0,b);CHKERRQ(ierr); 4432 ierr = VecNorm(lres,NORM_2,&kctx->lresid_last);CHKERRQ(ierr); 4433 ierr = VecDestroy(&lres);CHKERRQ(ierr); 4434 } 4435 } 4436 PetscFunctionReturn(0); 4437 } 4438 4439 #undef __FUNCT__ 4440 #define __FUNCT__ "SNES_KSPSolve" 4441 PetscErrorCode SNES_KSPSolve(SNES snes, KSP ksp, Vec b, Vec x) 4442 { 4443 PetscErrorCode ierr; 4444 4445 PetscFunctionBegin; 4446 if (snes->ksp_ewconv) { ierr = SNESKSPEW_PreSolve(snes,ksp,b,x);CHKERRQ(ierr); } 4447 ierr = KSPSolve(ksp,b,x);CHKERRQ(ierr); 4448 if (snes->ksp_ewconv) { ierr = SNESKSPEW_PostSolve(snes,ksp,b,x);CHKERRQ(ierr); } 4449 PetscFunctionReturn(0); 4450 } 4451 4452 #include <petsc-private/dmimpl.h> 4453 #undef __FUNCT__ 4454 #define __FUNCT__ "SNESSetDM" 4455 /*@ 4456 SNESSetDM - Sets the DM that may be used by some preconditioners 4457 4458 Logically Collective on SNES 4459 4460 Input Parameters: 4461 + snes - the preconditioner context 4462 - dm - the dm 4463 4464 Level: intermediate 4465 4466 4467 .seealso: SNESGetDM(), KSPSetDM(), KSPGetDM() 4468 @*/ 4469 PetscErrorCode SNESSetDM(SNES snes,DM dm) 4470 { 4471 PetscErrorCode ierr; 4472 KSP ksp; 4473 DMSNES sdm; 4474 4475 PetscFunctionBegin; 4476 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4477 if (dm) {ierr = PetscObjectReference((PetscObject)dm);CHKERRQ(ierr);} 4478 if (snes->dm) { /* Move the DMSNES context over to the new DM unless the new DM already has one */ 4479 if (snes->dm->dmsnes && snes->dmAuto && !dm->dmsnes) { 4480 ierr = DMCopyDMSNES(snes->dm,dm);CHKERRQ(ierr); 4481 ierr = DMGetDMSNES(snes->dm,&sdm);CHKERRQ(ierr); 4482 if (sdm->originaldm == snes->dm) { /* Grant write privileges to the replacement DM */ 4483 sdm->originaldm = dm; 4484 } 4485 } 4486 ierr = DMDestroy(&snes->dm);CHKERRQ(ierr); 4487 } 4488 snes->dm = dm; 4489 snes->dmAuto = PETSC_FALSE; 4490 ierr = SNESGetKSP(snes,&ksp);CHKERRQ(ierr); 4491 ierr = KSPSetDM(ksp,dm);CHKERRQ(ierr); 4492 ierr = KSPSetDMActive(ksp,PETSC_FALSE);CHKERRQ(ierr); 4493 if (snes->pc) { 4494 ierr = SNESSetDM(snes->pc, snes->dm);CHKERRQ(ierr); 4495 ierr = SNESSetPCSide(snes,snes->pcside);CHKERRQ(ierr); 4496 } 4497 PetscFunctionReturn(0); 4498 } 4499 4500 #undef __FUNCT__ 4501 #define __FUNCT__ "SNESGetDM" 4502 /*@ 4503 SNESGetDM - Gets the DM that may be used by some preconditioners 4504 4505 Not Collective but DM obtained is parallel on SNES 4506 4507 Input Parameter: 4508 . snes - the preconditioner context 4509 4510 Output Parameter: 4511 . dm - the dm 4512 4513 Level: intermediate 4514 4515 4516 .seealso: SNESSetDM(), KSPSetDM(), KSPGetDM() 4517 @*/ 4518 PetscErrorCode SNESGetDM(SNES snes,DM *dm) 4519 { 4520 PetscErrorCode ierr; 4521 4522 PetscFunctionBegin; 4523 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4524 if (!snes->dm) { 4525 ierr = DMShellCreate(((PetscObject)snes)->comm,&snes->dm);CHKERRQ(ierr); 4526 snes->dmAuto = PETSC_TRUE; 4527 } 4528 *dm = snes->dm; 4529 PetscFunctionReturn(0); 4530 } 4531 4532 #undef __FUNCT__ 4533 #define __FUNCT__ "SNESSetPC" 4534 /*@ 4535 SNESSetPC - Sets the nonlinear preconditioner to be used. 4536 4537 Collective on SNES 4538 4539 Input Parameters: 4540 + snes - iterative context obtained from SNESCreate() 4541 - pc - the preconditioner object 4542 4543 Notes: 4544 Use SNESGetPC() to retrieve the preconditioner context (for example, 4545 to configure it using the API). 4546 4547 Level: developer 4548 4549 .keywords: SNES, set, precondition 4550 .seealso: SNESGetPC() 4551 @*/ 4552 PetscErrorCode SNESSetPC(SNES snes, SNES pc) 4553 { 4554 PetscErrorCode ierr; 4555 4556 PetscFunctionBegin; 4557 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); 4558 PetscValidHeaderSpecific(pc, SNES_CLASSID, 2); 4559 PetscCheckSameComm(snes, 1, pc, 2); 4560 ierr = PetscObjectReference((PetscObject) pc);CHKERRQ(ierr); 4561 ierr = SNESDestroy(&snes->pc);CHKERRQ(ierr); 4562 snes->pc = pc; 4563 ierr = PetscLogObjectParent(snes, snes->pc);CHKERRQ(ierr); 4564 PetscFunctionReturn(0); 4565 } 4566 4567 #undef __FUNCT__ 4568 #define __FUNCT__ "SNESGetPC" 4569 /*@ 4570 SNESGetPC - Returns a pointer to the nonlinear preconditioning context set with SNESSetPC(). 4571 4572 Not Collective 4573 4574 Input Parameter: 4575 . snes - iterative context obtained from SNESCreate() 4576 4577 Output Parameter: 4578 . pc - preconditioner context 4579 4580 Level: developer 4581 4582 .keywords: SNES, get, preconditioner 4583 .seealso: SNESSetPC() 4584 @*/ 4585 PetscErrorCode SNESGetPC(SNES snes, SNES *pc) 4586 { 4587 PetscErrorCode ierr; 4588 const char *optionsprefix; 4589 4590 PetscFunctionBegin; 4591 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); 4592 PetscValidPointer(pc, 2); 4593 if (!snes->pc) { 4594 ierr = SNESCreate(((PetscObject) snes)->comm,&snes->pc);CHKERRQ(ierr); 4595 ierr = PetscObjectIncrementTabLevel((PetscObject)snes->pc,(PetscObject)snes,1);CHKERRQ(ierr); 4596 ierr = PetscLogObjectParent(snes,snes->pc);CHKERRQ(ierr); 4597 ierr = SNESGetOptionsPrefix(snes,&optionsprefix);CHKERRQ(ierr); 4598 ierr = SNESSetOptionsPrefix(snes->pc,optionsprefix);CHKERRQ(ierr); 4599 ierr = SNESAppendOptionsPrefix(snes->pc,"npc_");CHKERRQ(ierr); 4600 } 4601 *pc = snes->pc; 4602 PetscFunctionReturn(0); 4603 } 4604 4605 4606 #undef __FUNCT__ 4607 #define __FUNCT__ "SNESSetPCSide" 4608 /*@ 4609 SNESSetPCSide - Sets the preconditioning side. 4610 4611 Logically Collective on SNES 4612 4613 Input Parameter: 4614 . snes - iterative context obtained from SNESCreate() 4615 4616 Output Parameter: 4617 . side - the preconditioning side, where side is one of 4618 .vb 4619 PC_LEFT - left preconditioning (default) 4620 PC_RIGHT - right preconditioning 4621 .ve 4622 4623 Options Database Keys: 4624 . -snes_pc_side <right,left> 4625 4626 Level: intermediate 4627 4628 .keywords: SNES, set, right, left, side, preconditioner, flag 4629 4630 .seealso: SNESGetPCSide(), KSPSetPCSide() 4631 @*/ 4632 PetscErrorCode SNESSetPCSide(SNES snes,PCSide side) 4633 { 4634 PetscFunctionBegin; 4635 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4636 PetscValidLogicalCollectiveEnum(snes,side,2); 4637 snes->pcside = side; 4638 PetscFunctionReturn(0); 4639 } 4640 4641 #undef __FUNCT__ 4642 #define __FUNCT__ "SNESGetPCSide" 4643 /*@ 4644 SNESGetPCSide - Gets the preconditioning side. 4645 4646 Not Collective 4647 4648 Input Parameter: 4649 . snes - iterative context obtained from SNESCreate() 4650 4651 Output Parameter: 4652 . side - the preconditioning side, where side is one of 4653 .vb 4654 PC_LEFT - left preconditioning (default) 4655 PC_RIGHT - right preconditioning 4656 .ve 4657 4658 Level: intermediate 4659 4660 .keywords: SNES, get, right, left, side, preconditioner, flag 4661 4662 .seealso: SNESSetPCSide(), KSPGetPCSide() 4663 @*/ 4664 PetscErrorCode SNESGetPCSide(SNES snes,PCSide *side) 4665 { 4666 PetscFunctionBegin; 4667 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4668 PetscValidPointer(side,2); 4669 *side = snes->pcside; 4670 PetscFunctionReturn(0); 4671 } 4672 4673 #undef __FUNCT__ 4674 #define __FUNCT__ "SNESSetSNESLineSearch" 4675 /*@ 4676 SNESSetSNESLineSearch - Sets the linesearch on the SNES instance. 4677 4678 Collective on SNES 4679 4680 Input Parameters: 4681 + snes - iterative context obtained from SNESCreate() 4682 - linesearch - the linesearch object 4683 4684 Notes: 4685 Use SNESGetSNESLineSearch() to retrieve the preconditioner context (for example, 4686 to configure it using the API). 4687 4688 Level: developer 4689 4690 .keywords: SNES, set, linesearch 4691 .seealso: SNESGetSNESLineSearch() 4692 @*/ 4693 PetscErrorCode SNESSetSNESLineSearch(SNES snes, SNESLineSearch linesearch) 4694 { 4695 PetscErrorCode ierr; 4696 4697 PetscFunctionBegin; 4698 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); 4699 PetscValidHeaderSpecific(linesearch, SNESLINESEARCH_CLASSID, 2); 4700 PetscCheckSameComm(snes, 1, linesearch, 2); 4701 ierr = PetscObjectReference((PetscObject) linesearch);CHKERRQ(ierr); 4702 ierr = SNESLineSearchDestroy(&snes->linesearch);CHKERRQ(ierr); 4703 snes->linesearch = linesearch; 4704 ierr = PetscLogObjectParent(snes, snes->linesearch);CHKERRQ(ierr); 4705 PetscFunctionReturn(0); 4706 } 4707 4708 #undef __FUNCT__ 4709 #define __FUNCT__ "SNESGetSNESLineSearch" 4710 /*@C 4711 SNESGetSNESLineSearch - Returns a pointer to the line search context set with SNESSetLineSearch() 4712 or creates a default line search instance associated with the SNES and returns it. 4713 4714 Not Collective 4715 4716 Input Parameter: 4717 . snes - iterative context obtained from SNESCreate() 4718 4719 Output Parameter: 4720 . linesearch - linesearch context 4721 4722 Level: developer 4723 4724 .keywords: SNES, get, linesearch 4725 .seealso: SNESSetSNESLineSearch() 4726 @*/ 4727 PetscErrorCode SNESGetSNESLineSearch(SNES snes, SNESLineSearch *linesearch) 4728 { 4729 PetscErrorCode ierr; 4730 const char *optionsprefix; 4731 4732 PetscFunctionBegin; 4733 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); 4734 PetscValidPointer(linesearch, 2); 4735 if (!snes->linesearch) { 4736 ierr = SNESGetOptionsPrefix(snes, &optionsprefix);CHKERRQ(ierr); 4737 ierr = SNESLineSearchCreate(((PetscObject) snes)->comm, &snes->linesearch);CHKERRQ(ierr); 4738 ierr = SNESLineSearchSetSNES(snes->linesearch, snes);CHKERRQ(ierr); 4739 ierr = SNESLineSearchAppendOptionsPrefix(snes->linesearch, optionsprefix);CHKERRQ(ierr); 4740 ierr = PetscObjectIncrementTabLevel((PetscObject) snes->linesearch, (PetscObject) snes, 1);CHKERRQ(ierr); 4741 ierr = PetscLogObjectParent(snes, snes->linesearch);CHKERRQ(ierr); 4742 } 4743 *linesearch = snes->linesearch; 4744 PetscFunctionReturn(0); 4745 } 4746 4747 #if defined(PETSC_HAVE_MATLAB_ENGINE) 4748 #include <mex.h> 4749 4750 typedef struct {char *funcname; mxArray *ctx;} SNESMatlabContext; 4751 4752 #undef __FUNCT__ 4753 #define __FUNCT__ "SNESComputeFunction_Matlab" 4754 /* 4755 SNESComputeFunction_Matlab - Calls the function that has been set with 4756 SNESSetFunctionMatlab(). 4757 4758 Collective on SNES 4759 4760 Input Parameters: 4761 + snes - the SNES context 4762 - x - input vector 4763 4764 Output Parameter: 4765 . y - function vector, as set by SNESSetFunction() 4766 4767 Notes: 4768 SNESComputeFunction() is typically used within nonlinear solvers 4769 implementations, so most users would not generally call this routine 4770 themselves. 4771 4772 Level: developer 4773 4774 .keywords: SNES, nonlinear, compute, function 4775 4776 .seealso: SNESSetFunction(), SNESGetFunction() 4777 */ 4778 PetscErrorCode SNESComputeFunction_Matlab(SNES snes,Vec x,Vec y, void *ctx) 4779 { 4780 PetscErrorCode ierr; 4781 SNESMatlabContext *sctx = (SNESMatlabContext *)ctx; 4782 int nlhs = 1,nrhs = 5; 4783 mxArray *plhs[1],*prhs[5]; 4784 long long int lx = 0,ly = 0,ls = 0; 4785 4786 PetscFunctionBegin; 4787 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4788 PetscValidHeaderSpecific(x,VEC_CLASSID,2); 4789 PetscValidHeaderSpecific(y,VEC_CLASSID,3); 4790 PetscCheckSameComm(snes,1,x,2); 4791 PetscCheckSameComm(snes,1,y,3); 4792 4793 /* call Matlab function in ctx with arguments x and y */ 4794 4795 ierr = PetscMemcpy(&ls,&snes,sizeof(snes));CHKERRQ(ierr); 4796 ierr = PetscMemcpy(&lx,&x,sizeof(x));CHKERRQ(ierr); 4797 ierr = PetscMemcpy(&ly,&y,sizeof(x));CHKERRQ(ierr); 4798 prhs[0] = mxCreateDoubleScalar((double)ls); 4799 prhs[1] = mxCreateDoubleScalar((double)lx); 4800 prhs[2] = mxCreateDoubleScalar((double)ly); 4801 prhs[3] = mxCreateString(sctx->funcname); 4802 prhs[4] = sctx->ctx; 4803 ierr = mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscSNESComputeFunctionInternal");CHKERRQ(ierr); 4804 ierr = mxGetScalar(plhs[0]);CHKERRQ(ierr); 4805 mxDestroyArray(prhs[0]); 4806 mxDestroyArray(prhs[1]); 4807 mxDestroyArray(prhs[2]); 4808 mxDestroyArray(prhs[3]); 4809 mxDestroyArray(plhs[0]); 4810 PetscFunctionReturn(0); 4811 } 4812 4813 4814 #undef __FUNCT__ 4815 #define __FUNCT__ "SNESSetFunctionMatlab" 4816 /* 4817 SNESSetFunctionMatlab - Sets the function evaluation routine and function 4818 vector for use by the SNES routines in solving systems of nonlinear 4819 equations from MATLAB. Here the function is a string containing the name of a MATLAB function 4820 4821 Logically Collective on SNES 4822 4823 Input Parameters: 4824 + snes - the SNES context 4825 . r - vector to store function value 4826 - func - function evaluation routine 4827 4828 Calling sequence of func: 4829 $ func (SNES snes,Vec x,Vec f,void *ctx); 4830 4831 4832 Notes: 4833 The Newton-like methods typically solve linear systems of the form 4834 $ f'(x) x = -f(x), 4835 where f'(x) denotes the Jacobian matrix and f(x) is the function. 4836 4837 Level: beginner 4838 4839 .keywords: SNES, nonlinear, set, function 4840 4841 .seealso: SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESSetFunction() 4842 */ 4843 PetscErrorCode SNESSetFunctionMatlab(SNES snes,Vec r,const char *func,mxArray *ctx) 4844 { 4845 PetscErrorCode ierr; 4846 SNESMatlabContext *sctx; 4847 4848 PetscFunctionBegin; 4849 /* currently sctx is memory bleed */ 4850 ierr = PetscMalloc(sizeof(SNESMatlabContext),&sctx);CHKERRQ(ierr); 4851 ierr = PetscStrallocpy(func,&sctx->funcname);CHKERRQ(ierr); 4852 /* 4853 This should work, but it doesn't 4854 sctx->ctx = ctx; 4855 mexMakeArrayPersistent(sctx->ctx); 4856 */ 4857 sctx->ctx = mxDuplicateArray(ctx); 4858 ierr = SNESSetFunction(snes,r,SNESComputeFunction_Matlab,sctx);CHKERRQ(ierr); 4859 PetscFunctionReturn(0); 4860 } 4861 4862 #undef __FUNCT__ 4863 #define __FUNCT__ "SNESComputeJacobian_Matlab" 4864 /* 4865 SNESComputeJacobian_Matlab - Calls the function that has been set with 4866 SNESSetJacobianMatlab(). 4867 4868 Collective on SNES 4869 4870 Input Parameters: 4871 + snes - the SNES context 4872 . x - input vector 4873 . A, B - the matrices 4874 - ctx - user context 4875 4876 Output Parameter: 4877 . flag - structure of the matrix 4878 4879 Level: developer 4880 4881 .keywords: SNES, nonlinear, compute, function 4882 4883 .seealso: SNESSetFunction(), SNESGetFunction() 4884 @*/ 4885 PetscErrorCode SNESComputeJacobian_Matlab(SNES snes,Vec x,Mat *A,Mat *B,MatStructure *flag, void *ctx) 4886 { 4887 PetscErrorCode ierr; 4888 SNESMatlabContext *sctx = (SNESMatlabContext *)ctx; 4889 int nlhs = 2,nrhs = 6; 4890 mxArray *plhs[2],*prhs[6]; 4891 long long int lx = 0,lA = 0,ls = 0, lB = 0; 4892 4893 PetscFunctionBegin; 4894 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4895 PetscValidHeaderSpecific(x,VEC_CLASSID,2); 4896 4897 /* call Matlab function in ctx with arguments x and y */ 4898 4899 ierr = PetscMemcpy(&ls,&snes,sizeof(snes));CHKERRQ(ierr); 4900 ierr = PetscMemcpy(&lx,&x,sizeof(x));CHKERRQ(ierr); 4901 ierr = PetscMemcpy(&lA,A,sizeof(x));CHKERRQ(ierr); 4902 ierr = PetscMemcpy(&lB,B,sizeof(x));CHKERRQ(ierr); 4903 prhs[0] = mxCreateDoubleScalar((double)ls); 4904 prhs[1] = mxCreateDoubleScalar((double)lx); 4905 prhs[2] = mxCreateDoubleScalar((double)lA); 4906 prhs[3] = mxCreateDoubleScalar((double)lB); 4907 prhs[4] = mxCreateString(sctx->funcname); 4908 prhs[5] = sctx->ctx; 4909 ierr = mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscSNESComputeJacobianInternal");CHKERRQ(ierr); 4910 ierr = mxGetScalar(plhs[0]);CHKERRQ(ierr); 4911 *flag = (MatStructure) mxGetScalar(plhs[1]);CHKERRQ(ierr); 4912 mxDestroyArray(prhs[0]); 4913 mxDestroyArray(prhs[1]); 4914 mxDestroyArray(prhs[2]); 4915 mxDestroyArray(prhs[3]); 4916 mxDestroyArray(prhs[4]); 4917 mxDestroyArray(plhs[0]); 4918 mxDestroyArray(plhs[1]); 4919 PetscFunctionReturn(0); 4920 } 4921 4922 4923 #undef __FUNCT__ 4924 #define __FUNCT__ "SNESSetJacobianMatlab" 4925 /* 4926 SNESSetJacobianMatlab - Sets the Jacobian function evaluation routine and two empty Jacobian matrices 4927 vector for use by the SNES routines in solving systems of nonlinear 4928 equations from MATLAB. Here the function is a string containing the name of a MATLAB function 4929 4930 Logically Collective on SNES 4931 4932 Input Parameters: 4933 + snes - the SNES context 4934 . A,B - Jacobian matrices 4935 . func - function evaluation routine 4936 - ctx - user context 4937 4938 Calling sequence of func: 4939 $ flag = func (SNES snes,Vec x,Mat A,Mat B,void *ctx); 4940 4941 4942 Level: developer 4943 4944 .keywords: SNES, nonlinear, set, function 4945 4946 .seealso: SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESSetFunction() 4947 */ 4948 PetscErrorCode SNESSetJacobianMatlab(SNES snes,Mat A,Mat B,const char *func,mxArray *ctx) 4949 { 4950 PetscErrorCode ierr; 4951 SNESMatlabContext *sctx; 4952 4953 PetscFunctionBegin; 4954 /* currently sctx is memory bleed */ 4955 ierr = PetscMalloc(sizeof(SNESMatlabContext),&sctx);CHKERRQ(ierr); 4956 ierr = PetscStrallocpy(func,&sctx->funcname);CHKERRQ(ierr); 4957 /* 4958 This should work, but it doesn't 4959 sctx->ctx = ctx; 4960 mexMakeArrayPersistent(sctx->ctx); 4961 */ 4962 sctx->ctx = mxDuplicateArray(ctx); 4963 ierr = SNESSetJacobian(snes,A,B,SNESComputeJacobian_Matlab,sctx);CHKERRQ(ierr); 4964 PetscFunctionReturn(0); 4965 } 4966 4967 #undef __FUNCT__ 4968 #define __FUNCT__ "SNESMonitor_Matlab" 4969 /* 4970 SNESMonitor_Matlab - Calls the function that has been set with SNESMonitorSetMatlab(). 4971 4972 Collective on SNES 4973 4974 .seealso: SNESSetFunction(), SNESGetFunction() 4975 @*/ 4976 PetscErrorCode SNESMonitor_Matlab(SNES snes,PetscInt it, PetscReal fnorm, void *ctx) 4977 { 4978 PetscErrorCode ierr; 4979 SNESMatlabContext *sctx = (SNESMatlabContext *)ctx; 4980 int nlhs = 1,nrhs = 6; 4981 mxArray *plhs[1],*prhs[6]; 4982 long long int lx = 0,ls = 0; 4983 Vec x=snes->vec_sol; 4984 4985 PetscFunctionBegin; 4986 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4987 4988 ierr = PetscMemcpy(&ls,&snes,sizeof(snes));CHKERRQ(ierr); 4989 ierr = PetscMemcpy(&lx,&x,sizeof(x));CHKERRQ(ierr); 4990 prhs[0] = mxCreateDoubleScalar((double)ls); 4991 prhs[1] = mxCreateDoubleScalar((double)it); 4992 prhs[2] = mxCreateDoubleScalar((double)fnorm); 4993 prhs[3] = mxCreateDoubleScalar((double)lx); 4994 prhs[4] = mxCreateString(sctx->funcname); 4995 prhs[5] = sctx->ctx; 4996 ierr = mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscSNESMonitorInternal");CHKERRQ(ierr); 4997 ierr = mxGetScalar(plhs[0]);CHKERRQ(ierr); 4998 mxDestroyArray(prhs[0]); 4999 mxDestroyArray(prhs[1]); 5000 mxDestroyArray(prhs[2]); 5001 mxDestroyArray(prhs[3]); 5002 mxDestroyArray(prhs[4]); 5003 mxDestroyArray(plhs[0]); 5004 PetscFunctionReturn(0); 5005 } 5006 5007 5008 #undef __FUNCT__ 5009 #define __FUNCT__ "SNESMonitorSetMatlab" 5010 /* 5011 SNESMonitorSetMatlab - Sets the monitor function from MATLAB 5012 5013 Level: developer 5014 5015 .keywords: SNES, nonlinear, set, function 5016 5017 .seealso: SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESSetFunction() 5018 */ 5019 PetscErrorCode SNESMonitorSetMatlab(SNES snes,const char *func,mxArray *ctx) 5020 { 5021 PetscErrorCode ierr; 5022 SNESMatlabContext *sctx; 5023 5024 PetscFunctionBegin; 5025 /* currently sctx is memory bleed */ 5026 ierr = PetscMalloc(sizeof(SNESMatlabContext),&sctx);CHKERRQ(ierr); 5027 ierr = PetscStrallocpy(func,&sctx->funcname);CHKERRQ(ierr); 5028 /* 5029 This should work, but it doesn't 5030 sctx->ctx = ctx; 5031 mexMakeArrayPersistent(sctx->ctx); 5032 */ 5033 sctx->ctx = mxDuplicateArray(ctx); 5034 ierr = SNESMonitorSet(snes,SNESMonitor_Matlab,sctx,PETSC_NULL);CHKERRQ(ierr); 5035 PetscFunctionReturn(0); 5036 } 5037 5038 #endif 5039