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