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