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