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