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