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