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