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