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 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 2146 *flg = DIFFERENT_NONZERO_PATTERN; 2147 ierr = PetscLogEventBegin(SNES_JacobianEval,snes,X,*A,*B);CHKERRQ(ierr); 2148 2149 PetscStackPush("SNES user Jacobian function"); 2150 ierr = (*sdm->ops->computejacobian)(snes,X,A,B,flg,sdm->jacobianctx);CHKERRQ(ierr); 2151 PetscStackPop; 2152 2153 ierr = PetscLogEventEnd(SNES_JacobianEval,snes,X,*A,*B);CHKERRQ(ierr); 2154 2155 if (snes->lagpreconditioner == -2) { 2156 ierr = PetscInfo(snes,"Rebuilding preconditioner exactly once since lag is -2\n");CHKERRQ(ierr); 2157 2158 snes->lagpreconditioner = -1; 2159 } else if (snes->lagpreconditioner == -1) { 2160 *flg = SAME_PRECONDITIONER; 2161 ierr = PetscInfo(snes,"Reusing preconditioner because lag is -1\n");CHKERRQ(ierr); 2162 } else if (snes->lagpreconditioner > 1 && snes->iter % snes->lagpreconditioner) { 2163 *flg = SAME_PRECONDITIONER; 2164 ierr = PetscInfo2(snes,"Reusing preconditioner because lag is %D and SNES iteration is %D\n",snes->lagpreconditioner,snes->iter);CHKERRQ(ierr); 2165 } 2166 2167 /* make sure user returned a correct Jacobian and preconditioner */ 2168 /* PetscValidHeaderSpecific(*A,MAT_CLASSID,3); 2169 PetscValidHeaderSpecific(*B,MAT_CLASSID,4); */ 2170 { 2171 PetscBool flag = PETSC_FALSE,flag_draw = PETSC_FALSE,flag_contour = PETSC_FALSE,flag_operator = PETSC_FALSE; 2172 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_explicit",&flag,NULL);CHKERRQ(ierr); 2173 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_explicit_draw",&flag_draw,NULL);CHKERRQ(ierr); 2174 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_explicit_draw_contour",&flag_contour,NULL);CHKERRQ(ierr); 2175 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_operator",&flag_operator,NULL);CHKERRQ(ierr); 2176 if (flag || flag_draw || flag_contour) { 2177 Mat Bexp_mine = NULL,Bexp,FDexp; 2178 MatStructure mstruct; 2179 PetscViewer vdraw,vstdout; 2180 PetscBool flg; 2181 if (flag_operator) { 2182 ierr = MatComputeExplicitOperator(*A,&Bexp_mine);CHKERRQ(ierr); 2183 Bexp = Bexp_mine; 2184 } else { 2185 /* See if the preconditioning matrix can be viewed and added directly */ 2186 ierr = PetscObjectTypeCompareAny((PetscObject)*B,&flg,MATSEQAIJ,MATMPIAIJ,MATSEQDENSE,MATMPIDENSE,MATSEQBAIJ,MATMPIBAIJ,MATSEQSBAIJ,MATMPIBAIJ,"");CHKERRQ(ierr); 2187 if (flg) Bexp = *B; 2188 else { 2189 /* If the "preconditioning" matrix is itself MATSHELL or some other type without direct support */ 2190 ierr = MatComputeExplicitOperator(*B,&Bexp_mine);CHKERRQ(ierr); 2191 Bexp = Bexp_mine; 2192 } 2193 } 2194 ierr = MatConvert(Bexp,MATSAME,MAT_INITIAL_MATRIX,&FDexp);CHKERRQ(ierr); 2195 ierr = SNESDefaultComputeJacobian(snes,X,&FDexp,&FDexp,&mstruct,NULL);CHKERRQ(ierr); 2196 ierr = PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)snes),&vstdout);CHKERRQ(ierr); 2197 if (flag_draw || flag_contour) { 2198 ierr = PetscViewerDrawOpen(PetscObjectComm((PetscObject)snes),0,"Explicit Jacobians",PETSC_DECIDE,PETSC_DECIDE,300,300,&vdraw);CHKERRQ(ierr); 2199 if (flag_contour) {ierr = PetscViewerPushFormat(vdraw,PETSC_VIEWER_DRAW_CONTOUR);CHKERRQ(ierr);} 2200 } else vdraw = NULL; 2201 ierr = PetscViewerASCIIPrintf(vstdout,"Explicit %s\n",flag_operator ? "Jacobian" : "preconditioning Jacobian");CHKERRQ(ierr); 2202 if (flag) {ierr = MatView(Bexp,vstdout);CHKERRQ(ierr);} 2203 if (vdraw) {ierr = MatView(Bexp,vdraw);CHKERRQ(ierr);} 2204 ierr = PetscViewerASCIIPrintf(vstdout,"Finite difference Jacobian\n");CHKERRQ(ierr); 2205 if (flag) {ierr = MatView(FDexp,vstdout);CHKERRQ(ierr);} 2206 if (vdraw) {ierr = MatView(FDexp,vdraw);CHKERRQ(ierr);} 2207 ierr = MatAYPX(FDexp,-1.0,Bexp,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 2208 ierr = PetscViewerASCIIPrintf(vstdout,"User-provided matrix minus finite difference Jacobian\n");CHKERRQ(ierr); 2209 if (flag) {ierr = MatView(FDexp,vstdout);CHKERRQ(ierr);} 2210 if (vdraw) { /* Always use contour for the difference */ 2211 ierr = PetscViewerPushFormat(vdraw,PETSC_VIEWER_DRAW_CONTOUR);CHKERRQ(ierr); 2212 ierr = MatView(FDexp,vdraw);CHKERRQ(ierr); 2213 ierr = PetscViewerPopFormat(vdraw);CHKERRQ(ierr); 2214 } 2215 if (flag_contour) {ierr = PetscViewerPopFormat(vdraw);CHKERRQ(ierr);} 2216 ierr = PetscViewerDestroy(&vdraw);CHKERRQ(ierr); 2217 ierr = MatDestroy(&Bexp_mine);CHKERRQ(ierr); 2218 ierr = MatDestroy(&FDexp);CHKERRQ(ierr); 2219 } 2220 } 2221 { 2222 PetscBool flag = PETSC_FALSE,flag_display = PETSC_FALSE,flag_draw = PETSC_FALSE,flag_contour = PETSC_FALSE,flag_threshold = PETSC_FALSE; 2223 PetscReal threshold_atol = PETSC_SQRT_MACHINE_EPSILON,threshold_rtol = 10*PETSC_SQRT_MACHINE_EPSILON; 2224 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_coloring",&flag,NULL);CHKERRQ(ierr); 2225 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_coloring_display",&flag_display,NULL);CHKERRQ(ierr); 2226 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_coloring_draw",&flag_draw,NULL);CHKERRQ(ierr); 2227 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_coloring_draw_contour",&flag_contour,NULL);CHKERRQ(ierr); 2228 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_coloring_threshold",&flag_threshold,NULL);CHKERRQ(ierr); 2229 ierr = PetscOptionsGetReal(((PetscObject)snes)->prefix,"-snes_compare_coloring_threshold_rtol",&threshold_rtol,NULL);CHKERRQ(ierr); 2230 ierr = PetscOptionsGetReal(((PetscObject)snes)->prefix,"-snes_compare_coloring_threshold_atol",&threshold_atol,NULL);CHKERRQ(ierr); 2231 if (flag || flag_display || flag_draw || flag_contour || flag_threshold) { 2232 Mat Bfd; 2233 MatStructure mstruct; 2234 PetscViewer vdraw,vstdout; 2235 ISColoring iscoloring; 2236 MatFDColoring matfdcoloring; 2237 PetscErrorCode (*func)(SNES,Vec,Vec,void*); 2238 void *funcctx; 2239 PetscReal norm1,norm2,normmax; 2240 2241 ierr = MatDuplicate(*B,MAT_DO_NOT_COPY_VALUES,&Bfd);CHKERRQ(ierr); 2242 ierr = MatGetColoring(Bfd,MATCOLORINGSL,&iscoloring);CHKERRQ(ierr); 2243 ierr = MatFDColoringCreate(Bfd,iscoloring,&matfdcoloring);CHKERRQ(ierr); 2244 ierr = ISColoringDestroy(&iscoloring);CHKERRQ(ierr); 2245 2246 /* This method of getting the function is currently unreliable since it doesn't work for DM local functions. */ 2247 ierr = SNESGetFunction(snes,NULL,&func,&funcctx);CHKERRQ(ierr); 2248 ierr = MatFDColoringSetFunction(matfdcoloring,(PetscErrorCode (*)(void))func,funcctx);CHKERRQ(ierr); 2249 ierr = PetscObjectSetOptionsPrefix((PetscObject)matfdcoloring,((PetscObject)snes)->prefix);CHKERRQ(ierr); 2250 ierr = PetscObjectAppendOptionsPrefix((PetscObject)matfdcoloring,"coloring_");CHKERRQ(ierr); 2251 ierr = MatFDColoringSetFromOptions(matfdcoloring);CHKERRQ(ierr); 2252 ierr = MatFDColoringApply(Bfd,matfdcoloring,X,&mstruct,snes);CHKERRQ(ierr); 2253 ierr = MatFDColoringDestroy(&matfdcoloring);CHKERRQ(ierr); 2254 2255 ierr = PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)snes),&vstdout);CHKERRQ(ierr); 2256 if (flag_draw || flag_contour) { 2257 ierr = PetscViewerDrawOpen(PetscObjectComm((PetscObject)snes),0,"Colored Jacobians",PETSC_DECIDE,PETSC_DECIDE,300,300,&vdraw);CHKERRQ(ierr); 2258 if (flag_contour) {ierr = PetscViewerPushFormat(vdraw,PETSC_VIEWER_DRAW_CONTOUR);CHKERRQ(ierr);} 2259 } else vdraw = NULL; 2260 ierr = PetscViewerASCIIPrintf(vstdout,"Explicit preconditioning Jacobian\n");CHKERRQ(ierr); 2261 if (flag_display) {ierr = MatView(*B,vstdout);CHKERRQ(ierr);} 2262 if (vdraw) {ierr = MatView(*B,vdraw);CHKERRQ(ierr);} 2263 ierr = PetscViewerASCIIPrintf(vstdout,"Colored Finite difference Jacobian\n");CHKERRQ(ierr); 2264 if (flag_display) {ierr = MatView(Bfd,vstdout);CHKERRQ(ierr);} 2265 if (vdraw) {ierr = MatView(Bfd,vdraw);CHKERRQ(ierr);} 2266 ierr = MatAYPX(Bfd,-1.0,*B,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 2267 ierr = MatNorm(Bfd,NORM_1,&norm1);CHKERRQ(ierr); 2268 ierr = MatNorm(Bfd,NORM_FROBENIUS,&norm2);CHKERRQ(ierr); 2269 ierr = MatNorm(Bfd,NORM_MAX,&normmax);CHKERRQ(ierr); 2270 ierr = PetscViewerASCIIPrintf(vstdout,"User-provided matrix minus finite difference Jacobian, norm1=%G normFrob=%G normmax=%G\n",norm1,norm2,normmax);CHKERRQ(ierr); 2271 if (flag_display) {ierr = MatView(Bfd,vstdout);CHKERRQ(ierr);} 2272 if (vdraw) { /* Always use contour for the difference */ 2273 ierr = PetscViewerPushFormat(vdraw,PETSC_VIEWER_DRAW_CONTOUR);CHKERRQ(ierr); 2274 ierr = MatView(Bfd,vdraw);CHKERRQ(ierr); 2275 ierr = PetscViewerPopFormat(vdraw);CHKERRQ(ierr); 2276 } 2277 if (flag_contour) {ierr = PetscViewerPopFormat(vdraw);CHKERRQ(ierr);} 2278 2279 if (flag_threshold) { 2280 PetscInt bs,rstart,rend,i; 2281 ierr = MatGetBlockSize(*B,&bs);CHKERRQ(ierr); 2282 ierr = MatGetOwnershipRange(*B,&rstart,&rend);CHKERRQ(ierr); 2283 for (i=rstart; i<rend; i++) { 2284 const PetscScalar *ba,*ca; 2285 const PetscInt *bj,*cj; 2286 PetscInt bn,cn,j,maxentrycol = -1,maxdiffcol = -1,maxrdiffcol = -1; 2287 PetscReal maxentry = 0,maxdiff = 0,maxrdiff = 0; 2288 ierr = MatGetRow(*B,i,&bn,&bj,&ba);CHKERRQ(ierr); 2289 ierr = MatGetRow(Bfd,i,&cn,&cj,&ca);CHKERRQ(ierr); 2290 if (bn != cn) SETERRQ(((PetscObject)*A)->comm,PETSC_ERR_PLIB,"Unexpected different nonzero pattern in -snes_compare_coloring_threshold"); 2291 for (j=0; j<bn; j++) { 2292 PetscReal rdiff = PetscAbsScalar(ca[j]) / (threshold_atol + threshold_rtol*PetscAbsScalar(ba[j])); 2293 if (PetscAbsScalar(ba[j]) > PetscAbs(maxentry)) { 2294 maxentrycol = bj[j]; 2295 maxentry = PetscRealPart(ba[j]); 2296 } 2297 if (PetscAbsScalar(ca[j]) > PetscAbs(maxdiff)) { 2298 maxdiffcol = bj[j]; 2299 maxdiff = PetscRealPart(ca[j]); 2300 } 2301 if (rdiff > maxrdiff) { 2302 maxrdiffcol = bj[j]; 2303 maxrdiff = rdiff; 2304 } 2305 } 2306 if (maxrdiff > 1) { 2307 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); 2308 for (j=0; j<bn; j++) { 2309 PetscReal rdiff; 2310 rdiff = PetscAbsScalar(ca[j]) / (threshold_atol + threshold_rtol*PetscAbsScalar(ba[j])); 2311 if (rdiff > 1) { 2312 ierr = PetscViewerASCIIPrintf(vstdout," (%D,%G:%G)",bj[j],PetscRealPart(ba[j]),PetscRealPart(ca[j]));CHKERRQ(ierr); 2313 } 2314 } 2315 ierr = PetscViewerASCIIPrintf(vstdout,"\n",i,maxentry,maxdiff,maxrdiff);CHKERRQ(ierr); 2316 } 2317 ierr = MatRestoreRow(*B,i,&bn,&bj,&ba);CHKERRQ(ierr); 2318 ierr = MatRestoreRow(Bfd,i,&cn,&cj,&ca);CHKERRQ(ierr); 2319 } 2320 } 2321 ierr = PetscViewerDestroy(&vdraw);CHKERRQ(ierr); 2322 ierr = MatDestroy(&Bfd);CHKERRQ(ierr); 2323 } 2324 } 2325 PetscFunctionReturn(0); 2326 } 2327 2328 /*MC 2329 SNESJacobianFunction - function used to convey the nonlinear Jacobian of the function to be solved by SNES 2330 2331 Synopsis: 2332 #include "petscsnes.h" 2333 $ SNESJacobianFunction(SNES snes,Vec x,Mat *A,Mat *B,int *flag,void *ctx); 2334 2335 + x - input vector 2336 . A - Jacobian matrix 2337 . B - preconditioner matrix, usually the same as A 2338 . flag - flag indicating information about the preconditioner matrix 2339 structure (same as flag in KSPSetOperators()), one of SAME_NONZERO_PATTERN,DIFFERENT_NONZERO_PATTERN,SAME_PRECONDITIONER 2340 - ctx - [optional] user-defined Jacobian context 2341 2342 .seealso: SNESSetFunction(), SNESGetFunction(), SNESSetJacobian(), SNESGetJacobian() 2343 M*/ 2344 2345 #undef __FUNCT__ 2346 #define __FUNCT__ "SNESSetJacobian" 2347 /*@C 2348 SNESSetJacobian - Sets the function to compute Jacobian as well as the 2349 location to store the matrix. 2350 2351 Logically Collective on SNES and Mat 2352 2353 Input Parameters: 2354 + snes - the SNES context 2355 . A - Jacobian matrix 2356 . B - preconditioner matrix (usually same as the Jacobian) 2357 . SNESJacobianFunction - Jacobian evaluation routine (if NULL then SNES retains any previously set value) 2358 - ctx - [optional] user-defined context for private data for the 2359 Jacobian evaluation routine (may be NULL) (if NULL then SNES retains any previously set value) 2360 2361 Notes: 2362 See KSPSetOperators() for important information about setting the flag 2363 output parameter in the routine func(). Be sure to read this information! 2364 2365 The routine func() takes Mat * as the matrix arguments rather than Mat. 2366 This allows the Jacobian evaluation routine to replace A and/or B with a 2367 completely new new matrix structure (not just different matrix elements) 2368 when appropriate, for instance, if the nonzero structure is changing 2369 throughout the global iterations. 2370 2371 If the A matrix and B matrix are different you must call MatAssemblyBegin/End() on 2372 each matrix. 2373 2374 If using SNESDefaultComputeJacobianColor() to assemble a Jacobian, the ctx argument 2375 must be a MatFDColoring. 2376 2377 Other defect-correction schemes can be used by computing a different matrix in place of the Jacobian. One common 2378 example is to use the "Picard linearization" which only differentiates through the highest order parts of each term. 2379 2380 Level: beginner 2381 2382 .keywords: SNES, nonlinear, set, Jacobian, matrix 2383 2384 .seealso: KSPSetOperators(), SNESSetFunction(), MatMFFDComputeJacobian(), SNESDefaultComputeJacobianColor(), MatStructure, SNESJacobianFunction 2385 @*/ 2386 PetscErrorCode SNESSetJacobian(SNES snes,Mat A,Mat B,PetscErrorCode (*SNESJacobianFunction)(SNES,Vec,Mat*,Mat*,MatStructure*,void*),void *ctx) 2387 { 2388 PetscErrorCode ierr; 2389 DM dm; 2390 2391 PetscFunctionBegin; 2392 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2393 if (A) PetscValidHeaderSpecific(A,MAT_CLASSID,2); 2394 if (B) PetscValidHeaderSpecific(B,MAT_CLASSID,3); 2395 if (A) PetscCheckSameComm(snes,1,A,2); 2396 if (B) PetscCheckSameComm(snes,1,B,3); 2397 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 2398 ierr = DMSNESSetJacobian(dm,SNESJacobianFunction,ctx);CHKERRQ(ierr); 2399 if (A) { 2400 ierr = PetscObjectReference((PetscObject)A);CHKERRQ(ierr); 2401 ierr = MatDestroy(&snes->jacobian);CHKERRQ(ierr); 2402 2403 snes->jacobian = A; 2404 } 2405 if (B) { 2406 ierr = PetscObjectReference((PetscObject)B);CHKERRQ(ierr); 2407 ierr = MatDestroy(&snes->jacobian_pre);CHKERRQ(ierr); 2408 2409 snes->jacobian_pre = B; 2410 } 2411 PetscFunctionReturn(0); 2412 } 2413 2414 #undef __FUNCT__ 2415 #define __FUNCT__ "SNESGetJacobian" 2416 /*@C 2417 SNESGetJacobian - Returns the Jacobian matrix and optionally the user 2418 provided context for evaluating the Jacobian. 2419 2420 Not Collective, but Mat object will be parallel if SNES object is 2421 2422 Input Parameter: 2423 . snes - the nonlinear solver context 2424 2425 Output Parameters: 2426 + A - location to stash Jacobian matrix (or NULL) 2427 . B - location to stash preconditioner matrix (or NULL) 2428 . SNESJacobianFunction - location to put Jacobian function (or NULL) 2429 - ctx - location to stash Jacobian ctx (or NULL) 2430 2431 Level: advanced 2432 2433 .seealso: SNESSetJacobian(), SNESComputeJacobian() 2434 @*/ 2435 PetscErrorCode SNESGetJacobian(SNES snes,Mat *A,Mat *B,PetscErrorCode (**SNESJacobianFunction)(SNES,Vec,Mat*,Mat*,MatStructure*,void*),void **ctx) 2436 { 2437 PetscErrorCode ierr; 2438 DM dm; 2439 DMSNES sdm; 2440 2441 PetscFunctionBegin; 2442 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2443 if (A) *A = snes->jacobian; 2444 if (B) *B = snes->jacobian_pre; 2445 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 2446 ierr = DMGetDMSNES(dm,&sdm);CHKERRQ(ierr); 2447 if (SNESJacobianFunction) *SNESJacobianFunction = sdm->ops->computejacobian; 2448 if (ctx) *ctx = sdm->jacobianctx; 2449 PetscFunctionReturn(0); 2450 } 2451 2452 #undef __FUNCT__ 2453 #define __FUNCT__ "SNESSetUp" 2454 /*@ 2455 SNESSetUp - Sets up the internal data structures for the later use 2456 of a nonlinear solver. 2457 2458 Collective on SNES 2459 2460 Input Parameters: 2461 . snes - the SNES context 2462 2463 Notes: 2464 For basic use of the SNES solvers the user need not explicitly call 2465 SNESSetUp(), since these actions will automatically occur during 2466 the call to SNESSolve(). However, if one wishes to control this 2467 phase separately, SNESSetUp() should be called after SNESCreate() 2468 and optional routines of the form SNESSetXXX(), but before SNESSolve(). 2469 2470 Level: advanced 2471 2472 .keywords: SNES, nonlinear, setup 2473 2474 .seealso: SNESCreate(), SNESSolve(), SNESDestroy() 2475 @*/ 2476 PetscErrorCode SNESSetUp(SNES snes) 2477 { 2478 PetscErrorCode ierr; 2479 DM dm; 2480 DMSNES sdm; 2481 SNESLineSearch linesearch, pclinesearch; 2482 void *lsprectx,*lspostctx; 2483 PetscErrorCode (*precheck)(SNESLineSearch,Vec,Vec,PetscBool*,void*); 2484 PetscErrorCode (*postcheck)(SNESLineSearch,Vec,Vec,Vec,PetscBool*,PetscBool*,void*); 2485 PetscErrorCode (*func)(SNES,Vec,Vec,void*); 2486 Vec f,fpc; 2487 void *funcctx; 2488 PetscErrorCode (*jac)(SNES,Vec,Mat*,Mat*,MatStructure*,void*); 2489 void *jacctx,*appctx; 2490 Mat A,B; 2491 2492 PetscFunctionBegin; 2493 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2494 if (snes->setupcalled) PetscFunctionReturn(0); 2495 2496 if (!((PetscObject)snes)->type_name) { 2497 ierr = SNESSetType(snes,SNESNEWTONLS);CHKERRQ(ierr); 2498 } 2499 2500 ierr = SNESGetFunction(snes,&snes->vec_func,NULL,NULL);CHKERRQ(ierr); 2501 if (snes->vec_func == snes->vec_sol) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_IDN,"Solution vector cannot be function vector"); 2502 if (snes->vec_rhs == snes->vec_sol) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_IDN,"Solution vector cannot be right hand side vector"); 2503 2504 if (!snes->vec_sol_update /* && snes->vec_sol */) { 2505 ierr = VecDuplicate(snes->vec_sol,&snes->vec_sol_update);CHKERRQ(ierr); 2506 ierr = PetscLogObjectParent(snes,snes->vec_sol_update);CHKERRQ(ierr); 2507 } 2508 2509 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 2510 ierr = DMShellSetGlobalVector(snes->dm,snes->vec_sol);CHKERRQ(ierr); 2511 ierr = DMGetDMSNES(dm,&sdm);CHKERRQ(ierr); 2512 if (!sdm->ops->computefunction) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"Function never provided to SNES object"); 2513 if (!sdm->ops->computejacobian) { 2514 ierr = DMSNESSetJacobian(dm,SNESDefaultComputeJacobianColor,NULL);CHKERRQ(ierr); 2515 } 2516 if (!snes->vec_func) { 2517 ierr = DMCreateGlobalVector(dm,&snes->vec_func);CHKERRQ(ierr); 2518 } 2519 2520 if (!snes->ksp) { 2521 ierr = SNESGetKSP(snes, &snes->ksp);CHKERRQ(ierr); 2522 } 2523 2524 if (!snes->linesearch) { 2525 ierr = SNESGetSNESLineSearch(snes, &snes->linesearch);CHKERRQ(ierr); 2526 } 2527 2528 if (snes->pc && (snes->pcside == PC_LEFT)) { 2529 snes->mf = PETSC_TRUE; 2530 snes->mf_operator = PETSC_FALSE; 2531 } 2532 2533 if (snes->mf) { 2534 ierr = SNESSetUpMatrixFree_Private(snes, snes->mf_operator, snes->mf_version);CHKERRQ(ierr); 2535 } 2536 2537 if (snes->ops->usercompute && !snes->user) { 2538 ierr = (*snes->ops->usercompute)(snes,(void**)&snes->user);CHKERRQ(ierr); 2539 } 2540 2541 if (snes->pc) { 2542 /* copy the DM over */ 2543 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 2544 ierr = SNESSetDM(snes->pc,dm);CHKERRQ(ierr); 2545 2546 ierr = SNESGetFunction(snes,&f,&func,&funcctx);CHKERRQ(ierr); 2547 ierr = VecDuplicate(f,&fpc);CHKERRQ(ierr); 2548 ierr = SNESSetFunction(snes->pc,fpc,func,funcctx);CHKERRQ(ierr); 2549 ierr = SNESGetJacobian(snes,&A,&B,&jac,&jacctx);CHKERRQ(ierr); 2550 ierr = SNESSetJacobian(snes->pc,A,B,jac,jacctx);CHKERRQ(ierr); 2551 ierr = SNESGetApplicationContext(snes,&appctx);CHKERRQ(ierr); 2552 ierr = SNESSetApplicationContext(snes->pc,appctx);CHKERRQ(ierr); 2553 ierr = VecDestroy(&fpc);CHKERRQ(ierr); 2554 2555 /* copy the function pointers over */ 2556 ierr = PetscObjectCopyFortranFunctionPointers((PetscObject)snes,(PetscObject)snes->pc);CHKERRQ(ierr); 2557 2558 /* default to 1 iteration */ 2559 ierr = SNESSetTolerances(snes->pc,0.0,0.0,0.0,1,snes->pc->max_funcs);CHKERRQ(ierr); 2560 if (snes->pcside==PC_RIGHT) { 2561 ierr = SNESSetNormType(snes->pc,SNES_NORM_FINAL_ONLY);CHKERRQ(ierr); 2562 } else { 2563 ierr = SNESSetNormType(snes->pc,SNES_NORM_NONE);CHKERRQ(ierr); 2564 } 2565 ierr = SNESSetFromOptions(snes->pc);CHKERRQ(ierr); 2566 2567 /* copy the line search context over */ 2568 ierr = SNESGetSNESLineSearch(snes,&linesearch);CHKERRQ(ierr); 2569 ierr = SNESGetSNESLineSearch(snes->pc,&pclinesearch);CHKERRQ(ierr); 2570 ierr = SNESLineSearchGetPreCheck(linesearch,&precheck,&lsprectx);CHKERRQ(ierr); 2571 ierr = SNESLineSearchGetPostCheck(linesearch,&postcheck,&lspostctx);CHKERRQ(ierr); 2572 ierr = SNESLineSearchSetPreCheck(pclinesearch,precheck,lsprectx);CHKERRQ(ierr); 2573 ierr = SNESLineSearchSetPostCheck(pclinesearch,postcheck,lspostctx);CHKERRQ(ierr); 2574 ierr = PetscObjectCopyFortranFunctionPointers((PetscObject)linesearch, (PetscObject)pclinesearch);CHKERRQ(ierr); 2575 } 2576 2577 if (snes->ops->setup) { 2578 ierr = (*snes->ops->setup)(snes);CHKERRQ(ierr); 2579 } 2580 2581 snes->setupcalled = PETSC_TRUE; 2582 PetscFunctionReturn(0); 2583 } 2584 2585 #undef __FUNCT__ 2586 #define __FUNCT__ "SNESReset" 2587 /*@ 2588 SNESReset - Resets a SNES context to the snessetupcalled = 0 state and removes any allocated Vecs and Mats 2589 2590 Collective on SNES 2591 2592 Input Parameter: 2593 . snes - iterative context obtained from SNESCreate() 2594 2595 Level: intermediate 2596 2597 Notes: Also calls the application context destroy routine set with SNESSetComputeApplicationContext() 2598 2599 .keywords: SNES, destroy 2600 2601 .seealso: SNESCreate(), SNESSetUp(), SNESSolve() 2602 @*/ 2603 PetscErrorCode SNESReset(SNES snes) 2604 { 2605 PetscErrorCode ierr; 2606 2607 PetscFunctionBegin; 2608 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2609 if (snes->ops->userdestroy && snes->user) { 2610 ierr = (*snes->ops->userdestroy)((void**)&snes->user);CHKERRQ(ierr); 2611 snes->user = NULL; 2612 } 2613 if (snes->pc) { 2614 ierr = SNESReset(snes->pc);CHKERRQ(ierr); 2615 } 2616 2617 if (snes->ops->reset) { 2618 ierr = (*snes->ops->reset)(snes);CHKERRQ(ierr); 2619 } 2620 if (snes->ksp) { 2621 ierr = KSPReset(snes->ksp);CHKERRQ(ierr); 2622 } 2623 2624 if (snes->linesearch) { 2625 ierr = SNESLineSearchReset(snes->linesearch);CHKERRQ(ierr); 2626 } 2627 2628 ierr = VecDestroy(&snes->vec_rhs);CHKERRQ(ierr); 2629 ierr = VecDestroy(&snes->vec_sol);CHKERRQ(ierr); 2630 ierr = VecDestroy(&snes->vec_sol_update);CHKERRQ(ierr); 2631 ierr = VecDestroy(&snes->vec_func);CHKERRQ(ierr); 2632 ierr = MatDestroy(&snes->jacobian);CHKERRQ(ierr); 2633 ierr = MatDestroy(&snes->jacobian_pre);CHKERRQ(ierr); 2634 ierr = VecDestroyVecs(snes->nwork,&snes->work);CHKERRQ(ierr); 2635 ierr = VecDestroyVecs(snes->nvwork,&snes->vwork);CHKERRQ(ierr); 2636 2637 snes->nwork = snes->nvwork = 0; 2638 snes->setupcalled = PETSC_FALSE; 2639 PetscFunctionReturn(0); 2640 } 2641 2642 #undef __FUNCT__ 2643 #define __FUNCT__ "SNESDestroy" 2644 /*@ 2645 SNESDestroy - Destroys the nonlinear solver context that was created 2646 with SNESCreate(). 2647 2648 Collective on SNES 2649 2650 Input Parameter: 2651 . snes - the SNES context 2652 2653 Level: beginner 2654 2655 .keywords: SNES, nonlinear, destroy 2656 2657 .seealso: SNESCreate(), SNESSolve() 2658 @*/ 2659 PetscErrorCode SNESDestroy(SNES *snes) 2660 { 2661 PetscErrorCode ierr; 2662 2663 PetscFunctionBegin; 2664 if (!*snes) PetscFunctionReturn(0); 2665 PetscValidHeaderSpecific((*snes),SNES_CLASSID,1); 2666 if (--((PetscObject)(*snes))->refct > 0) {*snes = 0; PetscFunctionReturn(0);} 2667 2668 ierr = SNESReset((*snes));CHKERRQ(ierr); 2669 ierr = SNESDestroy(&(*snes)->pc);CHKERRQ(ierr); 2670 2671 /* if memory was published with AMS then destroy it */ 2672 ierr = PetscObjectDepublish((*snes));CHKERRQ(ierr); 2673 if ((*snes)->ops->destroy) {ierr = (*((*snes))->ops->destroy)((*snes));CHKERRQ(ierr);} 2674 2675 ierr = DMDestroy(&(*snes)->dm);CHKERRQ(ierr); 2676 ierr = KSPDestroy(&(*snes)->ksp);CHKERRQ(ierr); 2677 ierr = SNESLineSearchDestroy(&(*snes)->linesearch);CHKERRQ(ierr); 2678 2679 ierr = PetscFree((*snes)->kspconvctx);CHKERRQ(ierr); 2680 if ((*snes)->ops->convergeddestroy) { 2681 ierr = (*(*snes)->ops->convergeddestroy)((*snes)->cnvP);CHKERRQ(ierr); 2682 } 2683 if ((*snes)->conv_malloc) { 2684 ierr = PetscFree((*snes)->conv_hist);CHKERRQ(ierr); 2685 ierr = PetscFree((*snes)->conv_hist_its);CHKERRQ(ierr); 2686 } 2687 ierr = SNESMonitorCancel((*snes));CHKERRQ(ierr); 2688 ierr = PetscHeaderDestroy(snes);CHKERRQ(ierr); 2689 PetscFunctionReturn(0); 2690 } 2691 2692 /* ----------- Routines to set solver parameters ---------- */ 2693 2694 #undef __FUNCT__ 2695 #define __FUNCT__ "SNESSetLagPreconditioner" 2696 /*@ 2697 SNESSetLagPreconditioner - Determines when the preconditioner is rebuilt in the nonlinear solve. 2698 2699 Logically Collective on SNES 2700 2701 Input Parameters: 2702 + snes - the SNES context 2703 - lag - -1 indicates NEVER rebuild, 1 means rebuild every time the Jacobian is computed within a single nonlinear solve, 2 means every second time 2704 the Jacobian is built etc. -2 indicates rebuild preconditioner at next chance but then never rebuild after that 2705 2706 Options Database Keys: 2707 . -snes_lag_preconditioner <lag> 2708 2709 Notes: 2710 The default is 1 2711 The preconditioner is ALWAYS built in the first iteration of a nonlinear solve unless lag is -1 2712 If -1 is used before the very first nonlinear solve the preconditioner is still built because there is no previous preconditioner to use 2713 2714 Level: intermediate 2715 2716 .keywords: SNES, nonlinear, set, convergence, tolerances 2717 2718 .seealso: SNESSetTrustRegionTolerance(), SNESGetLagPreconditioner(), SNESSetLagJacobian(), SNESGetLagJacobian() 2719 2720 @*/ 2721 PetscErrorCode SNESSetLagPreconditioner(SNES snes,PetscInt lag) 2722 { 2723 PetscFunctionBegin; 2724 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2725 if (lag < -2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Lag must be -2, -1, 1 or greater"); 2726 if (!lag) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Lag cannot be 0"); 2727 PetscValidLogicalCollectiveInt(snes,lag,2); 2728 snes->lagpreconditioner = lag; 2729 PetscFunctionReturn(0); 2730 } 2731 2732 #undef __FUNCT__ 2733 #define __FUNCT__ "SNESSetGridSequence" 2734 /*@ 2735 SNESSetGridSequence - sets the number of steps of grid sequencing that SNES does 2736 2737 Logically Collective on SNES 2738 2739 Input Parameters: 2740 + snes - the SNES context 2741 - steps - the number of refinements to do, defaults to 0 2742 2743 Options Database Keys: 2744 . -snes_grid_sequence <steps> 2745 2746 Level: intermediate 2747 2748 Notes: 2749 Use SNESGetSolution() to extract the fine grid solution after grid sequencing. 2750 2751 .keywords: SNES, nonlinear, set, convergence, tolerances 2752 2753 .seealso: SNESSetTrustRegionTolerance(), SNESGetLagPreconditioner(), SNESSetLagJacobian(), SNESGetLagJacobian() 2754 2755 @*/ 2756 PetscErrorCode SNESSetGridSequence(SNES snes,PetscInt steps) 2757 { 2758 PetscFunctionBegin; 2759 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2760 PetscValidLogicalCollectiveInt(snes,steps,2); 2761 snes->gridsequence = steps; 2762 PetscFunctionReturn(0); 2763 } 2764 2765 #undef __FUNCT__ 2766 #define __FUNCT__ "SNESGetLagPreconditioner" 2767 /*@ 2768 SNESGetLagPreconditioner - Indicates how often the preconditioner is rebuilt 2769 2770 Not Collective 2771 2772 Input Parameter: 2773 . snes - the SNES context 2774 2775 Output Parameter: 2776 . lag - -1 indicates NEVER rebuild, 1 means rebuild every time the Jacobian is computed within a single nonlinear solve, 2 means every second time 2777 the Jacobian is built etc. -2 indicates rebuild preconditioner at next chance but then never rebuild after that 2778 2779 Options Database Keys: 2780 . -snes_lag_preconditioner <lag> 2781 2782 Notes: 2783 The default is 1 2784 The preconditioner is ALWAYS built in the first iteration of a nonlinear solve unless lag is -1 2785 2786 Level: intermediate 2787 2788 .keywords: SNES, nonlinear, set, convergence, tolerances 2789 2790 .seealso: SNESSetTrustRegionTolerance(), SNESSetLagPreconditioner() 2791 2792 @*/ 2793 PetscErrorCode SNESGetLagPreconditioner(SNES snes,PetscInt *lag) 2794 { 2795 PetscFunctionBegin; 2796 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2797 *lag = snes->lagpreconditioner; 2798 PetscFunctionReturn(0); 2799 } 2800 2801 #undef __FUNCT__ 2802 #define __FUNCT__ "SNESSetLagJacobian" 2803 /*@ 2804 SNESSetLagJacobian - Determines when the Jacobian is rebuilt in the nonlinear solve. See SNESSetLagPreconditioner() for determining how 2805 often the preconditioner is rebuilt. 2806 2807 Logically Collective on SNES 2808 2809 Input Parameters: 2810 + snes - the SNES context 2811 - lag - -1 indicates NEVER rebuild, 1 means rebuild every time the Jacobian is computed within a single nonlinear solve, 2 means every second time 2812 the Jacobian is built etc. -2 means rebuild at next chance but then never again 2813 2814 Options Database Keys: 2815 . -snes_lag_jacobian <lag> 2816 2817 Notes: 2818 The default is 1 2819 The Jacobian is ALWAYS built in the first iteration of a nonlinear solve unless lag is -1 2820 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 2821 at the next Newton step but never again (unless it is reset to another value) 2822 2823 Level: intermediate 2824 2825 .keywords: SNES, nonlinear, set, convergence, tolerances 2826 2827 .seealso: SNESSetTrustRegionTolerance(), SNESGetLagPreconditioner(), SNESSetLagPreconditioner(), SNESGetLagJacobian() 2828 2829 @*/ 2830 PetscErrorCode SNESSetLagJacobian(SNES snes,PetscInt lag) 2831 { 2832 PetscFunctionBegin; 2833 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2834 if (lag < -2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Lag must be -2, -1, 1 or greater"); 2835 if (!lag) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Lag cannot be 0"); 2836 PetscValidLogicalCollectiveInt(snes,lag,2); 2837 snes->lagjacobian = lag; 2838 PetscFunctionReturn(0); 2839 } 2840 2841 #undef __FUNCT__ 2842 #define __FUNCT__ "SNESGetLagJacobian" 2843 /*@ 2844 SNESGetLagJacobian - Indicates how often the Jacobian is rebuilt. See SNESGetLagPreconditioner() to determine when the preconditioner is rebuilt 2845 2846 Not Collective 2847 2848 Input Parameter: 2849 . snes - the SNES context 2850 2851 Output Parameter: 2852 . lag - -1 indicates NEVER rebuild, 1 means rebuild every time the Jacobian is computed within a single nonlinear solve, 2 means every second time 2853 the Jacobian is built etc. 2854 2855 Options Database Keys: 2856 . -snes_lag_jacobian <lag> 2857 2858 Notes: 2859 The default is 1 2860 The jacobian is ALWAYS built in the first iteration of a nonlinear solve unless lag is -1 2861 2862 Level: intermediate 2863 2864 .keywords: SNES, nonlinear, set, convergence, tolerances 2865 2866 .seealso: SNESSetTrustRegionTolerance(), SNESSetLagJacobian(), SNESSetLagPreconditioner(), SNESGetLagPreconditioner() 2867 2868 @*/ 2869 PetscErrorCode SNESGetLagJacobian(SNES snes,PetscInt *lag) 2870 { 2871 PetscFunctionBegin; 2872 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2873 *lag = snes->lagjacobian; 2874 PetscFunctionReturn(0); 2875 } 2876 2877 #undef __FUNCT__ 2878 #define __FUNCT__ "SNESSetTolerances" 2879 /*@ 2880 SNESSetTolerances - Sets various parameters used in convergence tests. 2881 2882 Logically Collective on SNES 2883 2884 Input Parameters: 2885 + snes - the SNES context 2886 . abstol - absolute convergence tolerance 2887 . rtol - relative convergence tolerance 2888 . stol - convergence tolerance in terms of the norm of the change in the solution between steps, || delta x || < stol*|| x || 2889 . maxit - maximum number of iterations 2890 - maxf - maximum number of function evaluations 2891 2892 Options Database Keys: 2893 + -snes_atol <abstol> - Sets abstol 2894 . -snes_rtol <rtol> - Sets rtol 2895 . -snes_stol <stol> - Sets stol 2896 . -snes_max_it <maxit> - Sets maxit 2897 - -snes_max_funcs <maxf> - Sets maxf 2898 2899 Notes: 2900 The default maximum number of iterations is 50. 2901 The default maximum number of function evaluations is 1000. 2902 2903 Level: intermediate 2904 2905 .keywords: SNES, nonlinear, set, convergence, tolerances 2906 2907 .seealso: SNESSetTrustRegionTolerance() 2908 @*/ 2909 PetscErrorCode SNESSetTolerances(SNES snes,PetscReal abstol,PetscReal rtol,PetscReal stol,PetscInt maxit,PetscInt maxf) 2910 { 2911 PetscFunctionBegin; 2912 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2913 PetscValidLogicalCollectiveReal(snes,abstol,2); 2914 PetscValidLogicalCollectiveReal(snes,rtol,3); 2915 PetscValidLogicalCollectiveReal(snes,stol,4); 2916 PetscValidLogicalCollectiveInt(snes,maxit,5); 2917 PetscValidLogicalCollectiveInt(snes,maxf,6); 2918 2919 if (abstol != PETSC_DEFAULT) { 2920 if (abstol < 0.0) SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_OUTOFRANGE,"Absolute tolerance %G must be non-negative",abstol); 2921 snes->abstol = abstol; 2922 } 2923 if (rtol != PETSC_DEFAULT) { 2924 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); 2925 snes->rtol = rtol; 2926 } 2927 if (stol != PETSC_DEFAULT) { 2928 if (stol < 0.0) SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_OUTOFRANGE,"Step tolerance %G must be non-negative",stol); 2929 snes->stol = stol; 2930 } 2931 if (maxit != PETSC_DEFAULT) { 2932 if (maxit < 0) SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_OUTOFRANGE,"Maximum number of iterations %D must be non-negative",maxit); 2933 snes->max_its = maxit; 2934 } 2935 if (maxf != PETSC_DEFAULT) { 2936 if (maxf < 0) SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_OUTOFRANGE,"Maximum number of function evaluations %D must be non-negative",maxf); 2937 snes->max_funcs = maxf; 2938 } 2939 snes->tolerancesset = PETSC_TRUE; 2940 PetscFunctionReturn(0); 2941 } 2942 2943 #undef __FUNCT__ 2944 #define __FUNCT__ "SNESGetTolerances" 2945 /*@ 2946 SNESGetTolerances - Gets various parameters used in convergence tests. 2947 2948 Not Collective 2949 2950 Input Parameters: 2951 + snes - the SNES context 2952 . atol - absolute convergence tolerance 2953 . rtol - relative convergence tolerance 2954 . stol - convergence tolerance in terms of the norm 2955 of the change in the solution between steps 2956 . maxit - maximum number of iterations 2957 - maxf - maximum number of function evaluations 2958 2959 Notes: 2960 The user can specify NULL for any parameter that is not needed. 2961 2962 Level: intermediate 2963 2964 .keywords: SNES, nonlinear, get, convergence, tolerances 2965 2966 .seealso: SNESSetTolerances() 2967 @*/ 2968 PetscErrorCode SNESGetTolerances(SNES snes,PetscReal *atol,PetscReal *rtol,PetscReal *stol,PetscInt *maxit,PetscInt *maxf) 2969 { 2970 PetscFunctionBegin; 2971 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2972 if (atol) *atol = snes->abstol; 2973 if (rtol) *rtol = snes->rtol; 2974 if (stol) *stol = snes->stol; 2975 if (maxit) *maxit = snes->max_its; 2976 if (maxf) *maxf = snes->max_funcs; 2977 PetscFunctionReturn(0); 2978 } 2979 2980 #undef __FUNCT__ 2981 #define __FUNCT__ "SNESSetTrustRegionTolerance" 2982 /*@ 2983 SNESSetTrustRegionTolerance - Sets the trust region parameter tolerance. 2984 2985 Logically Collective on SNES 2986 2987 Input Parameters: 2988 + snes - the SNES context 2989 - tol - tolerance 2990 2991 Options Database Key: 2992 . -snes_trtol <tol> - Sets tol 2993 2994 Level: intermediate 2995 2996 .keywords: SNES, nonlinear, set, trust region, tolerance 2997 2998 .seealso: SNESSetTolerances() 2999 @*/ 3000 PetscErrorCode SNESSetTrustRegionTolerance(SNES snes,PetscReal tol) 3001 { 3002 PetscFunctionBegin; 3003 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3004 PetscValidLogicalCollectiveReal(snes,tol,2); 3005 snes->deltatol = tol; 3006 PetscFunctionReturn(0); 3007 } 3008 3009 /* 3010 Duplicate the lg monitors for SNES from KSP; for some reason with 3011 dynamic libraries things don't work under Sun4 if we just use 3012 macros instead of functions 3013 */ 3014 #undef __FUNCT__ 3015 #define __FUNCT__ "SNESMonitorLGResidualNorm" 3016 PetscErrorCode SNESMonitorLGResidualNorm(SNES snes,PetscInt it,PetscReal norm,void *ctx) 3017 { 3018 PetscErrorCode ierr; 3019 3020 PetscFunctionBegin; 3021 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3022 ierr = KSPMonitorLGResidualNorm((KSP)snes,it,norm,ctx);CHKERRQ(ierr); 3023 PetscFunctionReturn(0); 3024 } 3025 3026 #undef __FUNCT__ 3027 #define __FUNCT__ "SNESMonitorLGCreate" 3028 PetscErrorCode SNESMonitorLGCreate(const char host[],const char label[],int x,int y,int m,int n,PetscDrawLG *draw) 3029 { 3030 PetscErrorCode ierr; 3031 3032 PetscFunctionBegin; 3033 ierr = KSPMonitorLGResidualNormCreate(host,label,x,y,m,n,draw);CHKERRQ(ierr); 3034 PetscFunctionReturn(0); 3035 } 3036 3037 #undef __FUNCT__ 3038 #define __FUNCT__ "SNESMonitorLGDestroy" 3039 PetscErrorCode SNESMonitorLGDestroy(PetscDrawLG *draw) 3040 { 3041 PetscErrorCode ierr; 3042 3043 PetscFunctionBegin; 3044 ierr = KSPMonitorLGResidualNormDestroy(draw);CHKERRQ(ierr); 3045 PetscFunctionReturn(0); 3046 } 3047 3048 extern PetscErrorCode SNESMonitorRange_Private(SNES,PetscInt,PetscReal*); 3049 #undef __FUNCT__ 3050 #define __FUNCT__ "SNESMonitorLGRange" 3051 PetscErrorCode SNESMonitorLGRange(SNES snes,PetscInt n,PetscReal rnorm,void *monctx) 3052 { 3053 PetscDrawLG lg; 3054 PetscErrorCode ierr; 3055 PetscReal x,y,per; 3056 PetscViewer v = (PetscViewer)monctx; 3057 static PetscReal prev; /* should be in the context */ 3058 PetscDraw draw; 3059 3060 PetscFunctionBegin; 3061 ierr = PetscViewerDrawGetDrawLG(v,0,&lg);CHKERRQ(ierr); 3062 if (!n) {ierr = PetscDrawLGReset(lg);CHKERRQ(ierr);} 3063 ierr = PetscDrawLGGetDraw(lg,&draw);CHKERRQ(ierr); 3064 ierr = PetscDrawSetTitle(draw,"Residual norm");CHKERRQ(ierr); 3065 x = (PetscReal)n; 3066 if (rnorm > 0.0) y = log10(rnorm); 3067 else y = -15.0; 3068 ierr = PetscDrawLGAddPoint(lg,&x,&y);CHKERRQ(ierr); 3069 if (n < 20 || !(n % 5)) { 3070 ierr = PetscDrawLGDraw(lg);CHKERRQ(ierr); 3071 } 3072 3073 ierr = PetscViewerDrawGetDrawLG(v,1,&lg);CHKERRQ(ierr); 3074 if (!n) {ierr = PetscDrawLGReset(lg);CHKERRQ(ierr);} 3075 ierr = PetscDrawLGGetDraw(lg,&draw);CHKERRQ(ierr); 3076 ierr = PetscDrawSetTitle(draw,"% elemts > .2*max elemt");CHKERRQ(ierr); 3077 ierr = SNESMonitorRange_Private(snes,n,&per);CHKERRQ(ierr); 3078 x = (PetscReal)n; 3079 y = 100.0*per; 3080 ierr = PetscDrawLGAddPoint(lg,&x,&y);CHKERRQ(ierr); 3081 if (n < 20 || !(n % 5)) { 3082 ierr = PetscDrawLGDraw(lg);CHKERRQ(ierr); 3083 } 3084 3085 ierr = PetscViewerDrawGetDrawLG(v,2,&lg);CHKERRQ(ierr); 3086 if (!n) {prev = rnorm;ierr = PetscDrawLGReset(lg);CHKERRQ(ierr);} 3087 ierr = PetscDrawLGGetDraw(lg,&draw);CHKERRQ(ierr); 3088 ierr = PetscDrawSetTitle(draw,"(norm -oldnorm)/oldnorm");CHKERRQ(ierr); 3089 x = (PetscReal)n; 3090 y = (prev - rnorm)/prev; 3091 ierr = PetscDrawLGAddPoint(lg,&x,&y);CHKERRQ(ierr); 3092 if (n < 20 || !(n % 5)) { 3093 ierr = PetscDrawLGDraw(lg);CHKERRQ(ierr); 3094 } 3095 3096 ierr = PetscViewerDrawGetDrawLG(v,3,&lg);CHKERRQ(ierr); 3097 if (!n) {ierr = PetscDrawLGReset(lg);CHKERRQ(ierr);} 3098 ierr = PetscDrawLGGetDraw(lg,&draw);CHKERRQ(ierr); 3099 ierr = PetscDrawSetTitle(draw,"(norm -oldnorm)/oldnorm*(% > .2 max)");CHKERRQ(ierr); 3100 x = (PetscReal)n; 3101 y = (prev - rnorm)/(prev*per); 3102 if (n > 2) { /*skip initial crazy value */ 3103 ierr = PetscDrawLGAddPoint(lg,&x,&y);CHKERRQ(ierr); 3104 } 3105 if (n < 20 || !(n % 5)) { 3106 ierr = PetscDrawLGDraw(lg);CHKERRQ(ierr); 3107 } 3108 prev = rnorm; 3109 PetscFunctionReturn(0); 3110 } 3111 3112 #undef __FUNCT__ 3113 #define __FUNCT__ "SNESMonitor" 3114 /*@ 3115 SNESMonitor - runs the user provided monitor routines, if they exist 3116 3117 Collective on SNES 3118 3119 Input Parameters: 3120 + snes - nonlinear solver context obtained from SNESCreate() 3121 . iter - iteration number 3122 - rnorm - relative norm of the residual 3123 3124 Notes: 3125 This routine is called by the SNES implementations. 3126 It does not typically need to be called by the user. 3127 3128 Level: developer 3129 3130 .seealso: SNESMonitorSet() 3131 @*/ 3132 PetscErrorCode SNESMonitor(SNES snes,PetscInt iter,PetscReal rnorm) 3133 { 3134 PetscErrorCode ierr; 3135 PetscInt i,n = snes->numbermonitors; 3136 3137 PetscFunctionBegin; 3138 for (i=0; i<n; i++) { 3139 ierr = (*snes->monitor[i])(snes,iter,rnorm,snes->monitorcontext[i]);CHKERRQ(ierr); 3140 } 3141 PetscFunctionReturn(0); 3142 } 3143 3144 /* ------------ Routines to set performance monitoring options ----------- */ 3145 3146 /*MC 3147 SNESMonitorFunction - functional form passed to SNESMonitorSet() to monitor convergence of nonlinear solver 3148 3149 Synopsis: 3150 #include "petscsnes.h" 3151 $ PetscErrorCode SNESMonitorFunction(SNES snes,PetscInt its, PetscReal norm,void *mctx) 3152 3153 + snes - the SNES context 3154 . its - iteration number 3155 . norm - 2-norm function value (may be estimated) 3156 - mctx - [optional] monitoring context 3157 3158 .seealso: SNESMonitorSet(), SNESMonitorGet() 3159 M*/ 3160 3161 #undef __FUNCT__ 3162 #define __FUNCT__ "SNESMonitorSet" 3163 /*@C 3164 SNESMonitorSet - Sets an ADDITIONAL function that is to be used at every 3165 iteration of the nonlinear solver to display the iteration's 3166 progress. 3167 3168 Logically Collective on SNES 3169 3170 Input Parameters: 3171 + snes - the SNES context 3172 . SNESMonitorFunction - monitoring routine 3173 . mctx - [optional] user-defined context for private data for the 3174 monitor routine (use NULL if no context is desired) 3175 - monitordestroy - [optional] routine that frees monitor context 3176 (may be NULL) 3177 3178 Options Database Keys: 3179 + -snes_monitor - sets SNESMonitorDefault() 3180 . -snes_monitor_lg_residualnorm - sets line graph monitor, 3181 uses SNESMonitorLGCreate() 3182 - -snes_monitor_cancel - cancels all monitors that have 3183 been hardwired into a code by 3184 calls to SNESMonitorSet(), but 3185 does not cancel those set via 3186 the options database. 3187 3188 Notes: 3189 Several different monitoring routines may be set by calling 3190 SNESMonitorSet() multiple times; all will be called in the 3191 order in which they were set. 3192 3193 Fortran notes: Only a single monitor function can be set for each SNES object 3194 3195 Level: intermediate 3196 3197 .keywords: SNES, nonlinear, set, monitor 3198 3199 .seealso: SNESMonitorDefault(), SNESMonitorCancel(), SNESMonitorFunction 3200 @*/ 3201 PetscErrorCode SNESMonitorSet(SNES snes,PetscErrorCode (*SNESMonitorFunction)(SNES,PetscInt,PetscReal,void*),void *mctx,PetscErrorCode (*monitordestroy)(void**)) 3202 { 3203 PetscInt i; 3204 PetscErrorCode ierr; 3205 3206 PetscFunctionBegin; 3207 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3208 if (snes->numbermonitors >= MAXSNESMONITORS) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Too many monitors set"); 3209 for (i=0; i<snes->numbermonitors;i++) { 3210 if (SNESMonitorFunction == snes->monitor[i] && monitordestroy == snes->monitordestroy[i] && mctx == snes->monitorcontext[i]) { 3211 if (monitordestroy) { 3212 ierr = (*monitordestroy)(&mctx);CHKERRQ(ierr); 3213 } 3214 PetscFunctionReturn(0); 3215 } 3216 } 3217 snes->monitor[snes->numbermonitors] = SNESMonitorFunction; 3218 snes->monitordestroy[snes->numbermonitors] = monitordestroy; 3219 snes->monitorcontext[snes->numbermonitors++] = (void*)mctx; 3220 PetscFunctionReturn(0); 3221 } 3222 3223 #undef __FUNCT__ 3224 #define __FUNCT__ "SNESMonitorCancel" 3225 /*@C 3226 SNESMonitorCancel - Clears all the monitor functions for a SNES object. 3227 3228 Logically Collective on SNES 3229 3230 Input Parameters: 3231 . snes - the SNES context 3232 3233 Options Database Key: 3234 . -snes_monitor_cancel - cancels all monitors that have been hardwired 3235 into a code by calls to SNESMonitorSet(), but does not cancel those 3236 set via the options database 3237 3238 Notes: 3239 There is no way to clear one specific monitor from a SNES object. 3240 3241 Level: intermediate 3242 3243 .keywords: SNES, nonlinear, set, monitor 3244 3245 .seealso: SNESMonitorDefault(), SNESMonitorSet() 3246 @*/ 3247 PetscErrorCode SNESMonitorCancel(SNES snes) 3248 { 3249 PetscErrorCode ierr; 3250 PetscInt i; 3251 3252 PetscFunctionBegin; 3253 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3254 for (i=0; i<snes->numbermonitors; i++) { 3255 if (snes->monitordestroy[i]) { 3256 ierr = (*snes->monitordestroy[i])(&snes->monitorcontext[i]);CHKERRQ(ierr); 3257 } 3258 } 3259 snes->numbermonitors = 0; 3260 PetscFunctionReturn(0); 3261 } 3262 3263 /*MC 3264 SNESConvergenceTestFunction - functional form used for testing of convergence of nonlinear solver 3265 3266 Synopsis: 3267 #include "petscsnes.h" 3268 $ PetscErrorCode SNESConvergenceTest(SNES snes,PetscInt it,PetscReal xnorm,PetscReal gnorm,PetscReal f,SNESConvergedReason *reason,void *cctx) 3269 3270 + snes - the SNES context 3271 . it - current iteration (0 is the first and is before any Newton step) 3272 . cctx - [optional] convergence context 3273 . reason - reason for convergence/divergence 3274 . xnorm - 2-norm of current iterate 3275 . gnorm - 2-norm of current step 3276 - f - 2-norm of function 3277 3278 3279 .seealso: SNESSetConvergenceTest(), SNESGetConvergenceTest() 3280 M*/ 3281 3282 #undef __FUNCT__ 3283 #define __FUNCT__ "SNESSetConvergenceTest" 3284 /*@C 3285 SNESSetConvergenceTest - Sets the function that is to be used 3286 to test for convergence of the nonlinear iterative solution. 3287 3288 Logically Collective on SNES 3289 3290 Input Parameters: 3291 + snes - the SNES context 3292 . SNESConvergenceTestFunction - routine to test for convergence 3293 . cctx - [optional] context for private data for the convergence routine (may be NULL) 3294 - destroy - [optional] destructor for the context (may be NULL; NULL_FUNCTION in Fortran) 3295 3296 Level: advanced 3297 3298 .keywords: SNES, nonlinear, set, convergence, test 3299 3300 .seealso: SNESDefaultConverged(), SNESSkipConverged(), SNESConvergenceTestFunction 3301 @*/ 3302 PetscErrorCode SNESSetConvergenceTest(SNES snes,PetscErrorCode (*SNESConvergenceTestFunction)(SNES,PetscInt,PetscReal,PetscReal,PetscReal,SNESConvergedReason*,void*),void *cctx,PetscErrorCode (*destroy)(void*)) 3303 { 3304 PetscErrorCode ierr; 3305 3306 PetscFunctionBegin; 3307 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3308 if (!SNESConvergenceTestFunction) SNESConvergenceTestFunction = SNESSkipConverged; 3309 if (snes->ops->convergeddestroy) { 3310 ierr = (*snes->ops->convergeddestroy)(snes->cnvP);CHKERRQ(ierr); 3311 } 3312 snes->ops->converged = SNESConvergenceTestFunction; 3313 snes->ops->convergeddestroy = destroy; 3314 snes->cnvP = cctx; 3315 PetscFunctionReturn(0); 3316 } 3317 3318 #undef __FUNCT__ 3319 #define __FUNCT__ "SNESGetConvergedReason" 3320 /*@ 3321 SNESGetConvergedReason - Gets the reason the SNES iteration was stopped. 3322 3323 Not Collective 3324 3325 Input Parameter: 3326 . snes - the SNES context 3327 3328 Output Parameter: 3329 . reason - negative value indicates diverged, positive value converged, see SNESConvergedReason or the 3330 manual pages for the individual convergence tests for complete lists 3331 3332 Level: intermediate 3333 3334 Notes: Can only be called after the call the SNESSolve() is complete. 3335 3336 .keywords: SNES, nonlinear, set, convergence, test 3337 3338 .seealso: SNESSetConvergenceTest(), SNESConvergedReason 3339 @*/ 3340 PetscErrorCode SNESGetConvergedReason(SNES snes,SNESConvergedReason *reason) 3341 { 3342 PetscFunctionBegin; 3343 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3344 PetscValidPointer(reason,2); 3345 *reason = snes->reason; 3346 PetscFunctionReturn(0); 3347 } 3348 3349 #undef __FUNCT__ 3350 #define __FUNCT__ "SNESSetConvergenceHistory" 3351 /*@ 3352 SNESSetConvergenceHistory - Sets the array used to hold the convergence history. 3353 3354 Logically Collective on SNES 3355 3356 Input Parameters: 3357 + snes - iterative context obtained from SNESCreate() 3358 . a - array to hold history, this array will contain the function norms computed at each step 3359 . its - integer array holds the number of linear iterations for each solve. 3360 . na - size of a and its 3361 - reset - PETSC_TRUE indicates each new nonlinear solve resets the history counter to zero, 3362 else it continues storing new values for new nonlinear solves after the old ones 3363 3364 Notes: 3365 If 'a' and 'its' are NULL then space is allocated for the history. If 'na' PETSC_DECIDE or PETSC_DEFAULT then a 3366 default array of length 10000 is allocated. 3367 3368 This routine is useful, e.g., when running a code for purposes 3369 of accurate performance monitoring, when no I/O should be done 3370 during the section of code that is being timed. 3371 3372 Level: intermediate 3373 3374 .keywords: SNES, set, convergence, history 3375 3376 .seealso: SNESGetConvergenceHistory() 3377 3378 @*/ 3379 PetscErrorCode SNESSetConvergenceHistory(SNES snes,PetscReal a[],PetscInt its[],PetscInt na,PetscBool reset) 3380 { 3381 PetscErrorCode ierr; 3382 3383 PetscFunctionBegin; 3384 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3385 if (na) PetscValidScalarPointer(a,2); 3386 if (its) PetscValidIntPointer(its,3); 3387 if (na == PETSC_DECIDE || na == PETSC_DEFAULT || !a) { 3388 if (na == PETSC_DECIDE || na == PETSC_DEFAULT) na = 1000; 3389 ierr = PetscMalloc(na*sizeof(PetscReal),&a);CHKERRQ(ierr); 3390 ierr = PetscMalloc(na*sizeof(PetscInt),&its);CHKERRQ(ierr); 3391 3392 snes->conv_malloc = PETSC_TRUE; 3393 } 3394 snes->conv_hist = a; 3395 snes->conv_hist_its = its; 3396 snes->conv_hist_max = na; 3397 snes->conv_hist_len = 0; 3398 snes->conv_hist_reset = reset; 3399 PetscFunctionReturn(0); 3400 } 3401 3402 #if defined(PETSC_HAVE_MATLAB_ENGINE) 3403 #include <engine.h> /* MATLAB include file */ 3404 #include <mex.h> /* MATLAB include file */ 3405 EXTERN_C_BEGIN 3406 #undef __FUNCT__ 3407 #define __FUNCT__ "SNESGetConvergenceHistoryMatlab" 3408 mxArray *SNESGetConvergenceHistoryMatlab(SNES snes) 3409 { 3410 mxArray *mat; 3411 PetscInt i; 3412 PetscReal *ar; 3413 3414 PetscFunctionBegin; 3415 mat = mxCreateDoubleMatrix(snes->conv_hist_len,1,mxREAL); 3416 ar = (PetscReal*) mxGetData(mat); 3417 for (i=0; i<snes->conv_hist_len; i++) ar[i] = snes->conv_hist[i]; 3418 PetscFunctionReturn(mat); 3419 } 3420 EXTERN_C_END 3421 #endif 3422 3423 #undef __FUNCT__ 3424 #define __FUNCT__ "SNESGetConvergenceHistory" 3425 /*@C 3426 SNESGetConvergenceHistory - Gets the array used to hold the convergence history. 3427 3428 Not Collective 3429 3430 Input Parameter: 3431 . snes - iterative context obtained from SNESCreate() 3432 3433 Output Parameters: 3434 . a - array to hold history 3435 . its - integer array holds the number of linear iterations (or 3436 negative if not converged) for each solve. 3437 - na - size of a and its 3438 3439 Notes: 3440 The calling sequence for this routine in Fortran is 3441 $ call SNESGetConvergenceHistory(SNES snes, integer na, integer ierr) 3442 3443 This routine is useful, e.g., when running a code for purposes 3444 of accurate performance monitoring, when no I/O should be done 3445 during the section of code that is being timed. 3446 3447 Level: intermediate 3448 3449 .keywords: SNES, get, convergence, history 3450 3451 .seealso: SNESSetConvergencHistory() 3452 3453 @*/ 3454 PetscErrorCode SNESGetConvergenceHistory(SNES snes,PetscReal *a[],PetscInt *its[],PetscInt *na) 3455 { 3456 PetscFunctionBegin; 3457 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3458 if (a) *a = snes->conv_hist; 3459 if (its) *its = snes->conv_hist_its; 3460 if (na) *na = snes->conv_hist_len; 3461 PetscFunctionReturn(0); 3462 } 3463 3464 #undef __FUNCT__ 3465 #define __FUNCT__ "SNESSetUpdate" 3466 /*@C 3467 SNESSetUpdate - Sets the general-purpose update function called 3468 at the beginning of every iteration of the nonlinear solve. Specifically 3469 it is called just before the Jacobian is "evaluated". 3470 3471 Logically Collective on SNES 3472 3473 Input Parameters: 3474 . snes - The nonlinear solver context 3475 . func - The function 3476 3477 Calling sequence of func: 3478 . func (SNES snes, PetscInt step); 3479 3480 . step - The current step of the iteration 3481 3482 Level: advanced 3483 3484 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() 3485 This is not used by most users. 3486 3487 .keywords: SNES, update 3488 3489 .seealso SNESDefaultUpdate(), SNESSetJacobian(), SNESSolve() 3490 @*/ 3491 PetscErrorCode SNESSetUpdate(SNES snes, PetscErrorCode (*func)(SNES, PetscInt)) 3492 { 3493 PetscFunctionBegin; 3494 PetscValidHeaderSpecific(snes, SNES_CLASSID,1); 3495 snes->ops->update = func; 3496 PetscFunctionReturn(0); 3497 } 3498 3499 #undef __FUNCT__ 3500 #define __FUNCT__ "SNESDefaultUpdate" 3501 /*@ 3502 SNESDefaultUpdate - The default update function which does nothing. 3503 3504 Not collective 3505 3506 Input Parameters: 3507 . snes - The nonlinear solver context 3508 . step - The current step of the iteration 3509 3510 Level: intermediate 3511 3512 .keywords: SNES, update 3513 .seealso SNESSetUpdate(), SNESDefaultRhsBC(), SNESDefaultShortolutionBC() 3514 @*/ 3515 PetscErrorCode SNESDefaultUpdate(SNES snes, PetscInt step) 3516 { 3517 PetscFunctionBegin; 3518 PetscFunctionReturn(0); 3519 } 3520 3521 #undef __FUNCT__ 3522 #define __FUNCT__ "SNESScaleStep_Private" 3523 /* 3524 SNESScaleStep_Private - Scales a step so that its length is less than the 3525 positive parameter delta. 3526 3527 Input Parameters: 3528 + snes - the SNES context 3529 . y - approximate solution of linear system 3530 . fnorm - 2-norm of current function 3531 - delta - trust region size 3532 3533 Output Parameters: 3534 + gpnorm - predicted function norm at the new point, assuming local 3535 linearization. The value is zero if the step lies within the trust 3536 region, and exceeds zero otherwise. 3537 - ynorm - 2-norm of the step 3538 3539 Note: 3540 For non-trust region methods such as SNESNEWTONLS, the parameter delta 3541 is set to be the maximum allowable step size. 3542 3543 .keywords: SNES, nonlinear, scale, step 3544 */ 3545 PetscErrorCode SNESScaleStep_Private(SNES snes,Vec y,PetscReal *fnorm,PetscReal *delta,PetscReal *gpnorm,PetscReal *ynorm) 3546 { 3547 PetscReal nrm; 3548 PetscScalar cnorm; 3549 PetscErrorCode ierr; 3550 3551 PetscFunctionBegin; 3552 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3553 PetscValidHeaderSpecific(y,VEC_CLASSID,2); 3554 PetscCheckSameComm(snes,1,y,2); 3555 3556 ierr = VecNorm(y,NORM_2,&nrm);CHKERRQ(ierr); 3557 if (nrm > *delta) { 3558 nrm = *delta/nrm; 3559 *gpnorm = (1.0 - nrm)*(*fnorm); 3560 cnorm = nrm; 3561 ierr = VecScale(y,cnorm);CHKERRQ(ierr); 3562 *ynorm = *delta; 3563 } else { 3564 *gpnorm = 0.0; 3565 *ynorm = nrm; 3566 } 3567 PetscFunctionReturn(0); 3568 } 3569 3570 #undef __FUNCT__ 3571 #define __FUNCT__ "SNESSolve" 3572 /*@C 3573 SNESSolve - Solves a nonlinear system F(x) = b. 3574 Call SNESSolve() after calling SNESCreate() and optional routines of the form SNESSetXXX(). 3575 3576 Collective on SNES 3577 3578 Input Parameters: 3579 + snes - the SNES context 3580 . b - the constant part of the equation F(x) = b, or NULL to use zero. 3581 - x - the solution vector. 3582 3583 Notes: 3584 The user should initialize the vector,x, with the initial guess 3585 for the nonlinear solve prior to calling SNESSolve. In particular, 3586 to employ an initial guess of zero, the user should explicitly set 3587 this vector to zero by calling VecSet(). 3588 3589 Level: beginner 3590 3591 .keywords: SNES, nonlinear, solve 3592 3593 .seealso: SNESCreate(), SNESDestroy(), SNESSetFunction(), SNESSetJacobian(), SNESSetGridSequence(), SNESGetSolution() 3594 @*/ 3595 PetscErrorCode SNESSolve(SNES snes,Vec b,Vec x) 3596 { 3597 PetscErrorCode ierr; 3598 PetscBool flg; 3599 PetscViewer viewer; 3600 PetscInt grid; 3601 Vec xcreated = NULL; 3602 DM dm; 3603 PetscViewerFormat format; 3604 3605 PetscFunctionBegin; 3606 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3607 if (x) PetscValidHeaderSpecific(x,VEC_CLASSID,3); 3608 if (x) PetscCheckSameComm(snes,1,x,3); 3609 if (b) PetscValidHeaderSpecific(b,VEC_CLASSID,2); 3610 if (b) PetscCheckSameComm(snes,1,b,2); 3611 3612 if (!x) { 3613 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 3614 ierr = DMCreateGlobalVector(dm,&xcreated);CHKERRQ(ierr); 3615 x = xcreated; 3616 } 3617 3618 for (grid=0; grid<snes->gridsequence; grid++) {ierr = PetscViewerASCIIPushTab(PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)snes)));CHKERRQ(ierr);} 3619 for (grid=0; grid<snes->gridsequence+1; grid++) { 3620 3621 /* set solution vector */ 3622 if (!grid) {ierr = PetscObjectReference((PetscObject)x);CHKERRQ(ierr);} 3623 ierr = VecDestroy(&snes->vec_sol);CHKERRQ(ierr); 3624 snes->vec_sol = x; 3625 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 3626 3627 /* set affine vector if provided */ 3628 if (b) { ierr = PetscObjectReference((PetscObject)b);CHKERRQ(ierr); } 3629 ierr = VecDestroy(&snes->vec_rhs);CHKERRQ(ierr); 3630 snes->vec_rhs = b; 3631 3632 ierr = SNESSetUp(snes);CHKERRQ(ierr); 3633 3634 if (!grid) { 3635 if (snes->ops->computeinitialguess) { 3636 ierr = (*snes->ops->computeinitialguess)(snes,snes->vec_sol,snes->initialguessP);CHKERRQ(ierr); 3637 } 3638 } 3639 3640 if (snes->conv_hist_reset) snes->conv_hist_len = 0; 3641 snes->nfuncs = 0; snes->linear_its = 0; snes->numFailures = 0; 3642 3643 ierr = PetscLogEventBegin(SNES_Solve,snes,0,0,0);CHKERRQ(ierr); 3644 ierr = (*snes->ops->solve)(snes);CHKERRQ(ierr); 3645 ierr = PetscLogEventEnd(SNES_Solve,snes,0,0,0);CHKERRQ(ierr); 3646 if (snes->domainerror) { 3647 snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; 3648 snes->domainerror = PETSC_FALSE; 3649 } 3650 if (!snes->reason) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Internal error, solver returned without setting converged reason"); 3651 3652 flg = PETSC_FALSE; 3653 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_test_local_min",&flg,NULL);CHKERRQ(ierr); 3654 if (flg && !PetscPreLoadingOn) { ierr = SNESTestLocalMin(snes);CHKERRQ(ierr); } 3655 if (snes->printreason) { 3656 ierr = PetscViewerASCIIAddTab(PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)snes)),((PetscObject)snes)->tablevel);CHKERRQ(ierr); 3657 if (snes->reason > 0) { 3658 ierr = PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)snes)),"Nonlinear solve converged due to %s iterations %D\n",SNESConvergedReasons[snes->reason],snes->iter);CHKERRQ(ierr); 3659 } else { 3660 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); 3661 } 3662 ierr = PetscViewerASCIISubtractTab(PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)snes)),((PetscObject)snes)->tablevel);CHKERRQ(ierr); 3663 } 3664 3665 if (snes->errorifnotconverged && snes->reason < 0) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_NOT_CONVERGED,"SNESSolve has not converged"); 3666 if (grid < snes->gridsequence) { 3667 DM fine; 3668 Vec xnew; 3669 Mat interp; 3670 3671 ierr = DMRefine(snes->dm,PetscObjectComm((PetscObject)snes),&fine);CHKERRQ(ierr); 3672 if (!fine) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_INCOMP,"DMRefine() did not perform any refinement, cannot continue grid sequencing"); 3673 ierr = DMCreateInterpolation(snes->dm,fine,&interp,NULL);CHKERRQ(ierr); 3674 ierr = DMCreateGlobalVector(fine,&xnew);CHKERRQ(ierr); 3675 ierr = MatInterpolate(interp,x,xnew);CHKERRQ(ierr); 3676 ierr = DMInterpolate(snes->dm,interp,fine);CHKERRQ(ierr); 3677 ierr = MatDestroy(&interp);CHKERRQ(ierr); 3678 x = xnew; 3679 3680 ierr = SNESReset(snes);CHKERRQ(ierr); 3681 ierr = SNESSetDM(snes,fine);CHKERRQ(ierr); 3682 ierr = DMDestroy(&fine);CHKERRQ(ierr); 3683 ierr = PetscViewerASCIIPopTab(PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)snes)));CHKERRQ(ierr); 3684 } 3685 } 3686 /* monitoring and viewing */ 3687 ierr = PetscOptionsGetViewer(PetscObjectComm((PetscObject)snes),((PetscObject)snes)->prefix,"-snes_view",&viewer,&format,&flg);CHKERRQ(ierr); 3688 if (flg && !PetscPreLoadingOn) { 3689 ierr = PetscViewerPushFormat(viewer,format);CHKERRQ(ierr); 3690 ierr = SNESView(snes,viewer);CHKERRQ(ierr); 3691 ierr = PetscViewerPopFormat(viewer);CHKERRQ(ierr); 3692 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 3693 } 3694 ierr = VecViewFromOptions(snes->vec_sol,"-snes_view_solution");CHKERRQ(ierr); 3695 3696 ierr = VecDestroy(&xcreated);CHKERRQ(ierr); 3697 PetscFunctionReturn(0); 3698 } 3699 3700 /* --------- Internal routines for SNES Package --------- */ 3701 3702 #undef __FUNCT__ 3703 #define __FUNCT__ "SNESSetType" 3704 /*@C 3705 SNESSetType - Sets the method for the nonlinear solver. 3706 3707 Collective on SNES 3708 3709 Input Parameters: 3710 + snes - the SNES context 3711 - type - a known method 3712 3713 Options Database Key: 3714 . -snes_type <type> - Sets the method; use -help for a list 3715 of available methods (for instance, newtonls or newtontr) 3716 3717 Notes: 3718 See "petsc/include/petscsnes.h" for available methods (for instance) 3719 + SNESNEWTONLS - Newton's method with line search 3720 (systems of nonlinear equations) 3721 . SNESNEWTONTR - Newton's method with trust region 3722 (systems of nonlinear equations) 3723 3724 Normally, it is best to use the SNESSetFromOptions() command and then 3725 set the SNES solver type from the options database rather than by using 3726 this routine. Using the options database provides the user with 3727 maximum flexibility in evaluating the many nonlinear solvers. 3728 The SNESSetType() routine is provided for those situations where it 3729 is necessary to set the nonlinear solver independently of the command 3730 line or options database. This might be the case, for example, when 3731 the choice of solver changes during the execution of the program, 3732 and the user's application is taking responsibility for choosing the 3733 appropriate method. 3734 3735 Level: intermediate 3736 3737 .keywords: SNES, set, type 3738 3739 .seealso: SNESType, SNESCreate() 3740 3741 @*/ 3742 PetscErrorCode SNESSetType(SNES snes,SNESType type) 3743 { 3744 PetscErrorCode ierr,(*r)(SNES); 3745 PetscBool match; 3746 3747 PetscFunctionBegin; 3748 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3749 PetscValidCharPointer(type,2); 3750 3751 ierr = PetscObjectTypeCompare((PetscObject)snes,type,&match);CHKERRQ(ierr); 3752 if (match) PetscFunctionReturn(0); 3753 3754 ierr = PetscFunctionListFind(PetscObjectComm((PetscObject)snes),SNESList,type,PETSC_TRUE,(void (**)(void)) &r);CHKERRQ(ierr); 3755 if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unable to find requested SNES type %s",type); 3756 /* Destroy the previous private SNES context */ 3757 if (snes->ops->destroy) { 3758 ierr = (*(snes)->ops->destroy)(snes);CHKERRQ(ierr); 3759 snes->ops->destroy = NULL; 3760 } 3761 /* Reinitialize function pointers in SNESOps structure */ 3762 snes->ops->setup = 0; 3763 snes->ops->solve = 0; 3764 snes->ops->view = 0; 3765 snes->ops->setfromoptions = 0; 3766 snes->ops->destroy = 0; 3767 /* Call the SNESCreate_XXX routine for this particular Nonlinear solver */ 3768 snes->setupcalled = PETSC_FALSE; 3769 3770 ierr = PetscObjectChangeTypeName((PetscObject)snes,type);CHKERRQ(ierr); 3771 ierr = (*r)(snes);CHKERRQ(ierr); 3772 #if defined(PETSC_HAVE_AMS) 3773 if (PetscAMSPublishAll) { 3774 ierr = PetscObjectAMSPublish((PetscObject)snes);CHKERRQ(ierr); 3775 } 3776 #endif 3777 PetscFunctionReturn(0); 3778 } 3779 3780 /* --------------------------------------------------------------------- */ 3781 #undef __FUNCT__ 3782 #define __FUNCT__ "SNESRegisterDestroy" 3783 /*@ 3784 SNESRegisterDestroy - Frees the list of nonlinear solvers that were 3785 registered by SNESRegisterDynamic(). 3786 3787 Not Collective 3788 3789 Level: advanced 3790 3791 .keywords: SNES, nonlinear, register, destroy 3792 3793 .seealso: SNESRegisterAll() 3794 @*/ 3795 PetscErrorCode SNESRegisterDestroy(void) 3796 { 3797 PetscErrorCode ierr; 3798 3799 PetscFunctionBegin; 3800 ierr = PetscFunctionListDestroy(&SNESList);CHKERRQ(ierr); 3801 3802 SNESRegisterAllCalled = PETSC_FALSE; 3803 PetscFunctionReturn(0); 3804 } 3805 3806 #undef __FUNCT__ 3807 #define __FUNCT__ "SNESGetType" 3808 /*@C 3809 SNESGetType - Gets the SNES method type and name (as a string). 3810 3811 Not Collective 3812 3813 Input Parameter: 3814 . snes - nonlinear solver context 3815 3816 Output Parameter: 3817 . type - SNES method (a character string) 3818 3819 Level: intermediate 3820 3821 .keywords: SNES, nonlinear, get, type, name 3822 @*/ 3823 PetscErrorCode SNESGetType(SNES snes,SNESType *type) 3824 { 3825 PetscFunctionBegin; 3826 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3827 PetscValidPointer(type,2); 3828 *type = ((PetscObject)snes)->type_name; 3829 PetscFunctionReturn(0); 3830 } 3831 3832 #undef __FUNCT__ 3833 #define __FUNCT__ "SNESGetSolution" 3834 /*@ 3835 SNESGetSolution - Returns the vector where the approximate solution is 3836 stored. This is the fine grid solution when using SNESSetGridSequence(). 3837 3838 Not Collective, but Vec is parallel if SNES is parallel 3839 3840 Input Parameter: 3841 . snes - the SNES context 3842 3843 Output Parameter: 3844 . x - the solution 3845 3846 Level: intermediate 3847 3848 .keywords: SNES, nonlinear, get, solution 3849 3850 .seealso: SNESGetSolutionUpdate(), SNESGetFunction() 3851 @*/ 3852 PetscErrorCode SNESGetSolution(SNES snes,Vec *x) 3853 { 3854 PetscFunctionBegin; 3855 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3856 PetscValidPointer(x,2); 3857 *x = snes->vec_sol; 3858 PetscFunctionReturn(0); 3859 } 3860 3861 #undef __FUNCT__ 3862 #define __FUNCT__ "SNESGetSolutionUpdate" 3863 /*@ 3864 SNESGetSolutionUpdate - Returns the vector where the solution update is 3865 stored. 3866 3867 Not Collective, but Vec is parallel if SNES is parallel 3868 3869 Input Parameter: 3870 . snes - the SNES context 3871 3872 Output Parameter: 3873 . x - the solution update 3874 3875 Level: advanced 3876 3877 .keywords: SNES, nonlinear, get, solution, update 3878 3879 .seealso: SNESGetSolution(), SNESGetFunction() 3880 @*/ 3881 PetscErrorCode SNESGetSolutionUpdate(SNES snes,Vec *x) 3882 { 3883 PetscFunctionBegin; 3884 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3885 PetscValidPointer(x,2); 3886 *x = snes->vec_sol_update; 3887 PetscFunctionReturn(0); 3888 } 3889 3890 #undef __FUNCT__ 3891 #define __FUNCT__ "SNESGetFunction" 3892 /*@C 3893 SNESGetFunction - Returns the vector where the function is stored. 3894 3895 Not Collective, but Vec is parallel if SNES is parallel. Collective if Vec is requested, but has not been created yet. 3896 3897 Input Parameter: 3898 . snes - the SNES context 3899 3900 Output Parameter: 3901 + r - the vector that is used to store residuals (or NULL if you don't want it) 3902 . SNESFunction- the function (or NULL if you don't want it) 3903 - ctx - the function context (or NULL if you don't want it) 3904 3905 Level: advanced 3906 3907 .keywords: SNES, nonlinear, get, function 3908 3909 .seealso: SNESSetFunction(), SNESGetSolution(), SNESFunction 3910 @*/ 3911 PetscErrorCode SNESGetFunction(SNES snes,Vec *r,PetscErrorCode (**SNESFunction)(SNES,Vec,Vec,void*),void **ctx) 3912 { 3913 PetscErrorCode ierr; 3914 DM dm; 3915 3916 PetscFunctionBegin; 3917 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3918 if (r) { 3919 if (!snes->vec_func) { 3920 if (snes->vec_rhs) { 3921 ierr = VecDuplicate(snes->vec_rhs,&snes->vec_func);CHKERRQ(ierr); 3922 } else if (snes->vec_sol) { 3923 ierr = VecDuplicate(snes->vec_sol,&snes->vec_func);CHKERRQ(ierr); 3924 } else if (snes->dm) { 3925 ierr = DMCreateGlobalVector(snes->dm,&snes->vec_func);CHKERRQ(ierr); 3926 } 3927 } 3928 *r = snes->vec_func; 3929 } 3930 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 3931 ierr = DMSNESGetFunction(dm,SNESFunction,ctx);CHKERRQ(ierr); 3932 PetscFunctionReturn(0); 3933 } 3934 3935 /*@C 3936 SNESGetGS - Returns the GS function and context. 3937 3938 Input Parameter: 3939 . snes - the SNES context 3940 3941 Output Parameter: 3942 + SNESGSFunction - the function (or NULL) 3943 - ctx - the function context (or NULL) 3944 3945 Level: advanced 3946 3947 .keywords: SNES, nonlinear, get, function 3948 3949 .seealso: SNESSetGS(), SNESGetFunction() 3950 @*/ 3951 3952 #undef __FUNCT__ 3953 #define __FUNCT__ "SNESGetGS" 3954 PetscErrorCode SNESGetGS (SNES snes, PetscErrorCode (**SNESGSFunction)(SNES, Vec, Vec, void*), void ** ctx) 3955 { 3956 PetscErrorCode ierr; 3957 DM dm; 3958 3959 PetscFunctionBegin; 3960 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3961 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 3962 ierr = DMSNESGetGS(dm,SNESGSFunction,ctx);CHKERRQ(ierr); 3963 PetscFunctionReturn(0); 3964 } 3965 3966 #undef __FUNCT__ 3967 #define __FUNCT__ "SNESSetOptionsPrefix" 3968 /*@C 3969 SNESSetOptionsPrefix - Sets the prefix used for searching for all 3970 SNES options in the database. 3971 3972 Logically Collective on SNES 3973 3974 Input Parameter: 3975 + snes - the SNES context 3976 - prefix - the prefix to prepend to all option names 3977 3978 Notes: 3979 A hyphen (-) must NOT be given at the beginning of the prefix name. 3980 The first character of all runtime options is AUTOMATICALLY the hyphen. 3981 3982 Level: advanced 3983 3984 .keywords: SNES, set, options, prefix, database 3985 3986 .seealso: SNESSetFromOptions() 3987 @*/ 3988 PetscErrorCode SNESSetOptionsPrefix(SNES snes,const char prefix[]) 3989 { 3990 PetscErrorCode ierr; 3991 3992 PetscFunctionBegin; 3993 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3994 ierr = PetscObjectSetOptionsPrefix((PetscObject)snes,prefix);CHKERRQ(ierr); 3995 if (!snes->ksp) {ierr = SNESGetKSP(snes,&snes->ksp);CHKERRQ(ierr);} 3996 if (snes->linesearch) { 3997 ierr = SNESGetSNESLineSearch(snes,&snes->linesearch);CHKERRQ(ierr); 3998 ierr = PetscObjectSetOptionsPrefix((PetscObject)snes->linesearch,prefix);CHKERRQ(ierr); 3999 } 4000 ierr = KSPSetOptionsPrefix(snes->ksp,prefix);CHKERRQ(ierr); 4001 PetscFunctionReturn(0); 4002 } 4003 4004 #undef __FUNCT__ 4005 #define __FUNCT__ "SNESAppendOptionsPrefix" 4006 /*@C 4007 SNESAppendOptionsPrefix - Appends to the prefix used for searching for all 4008 SNES options in the database. 4009 4010 Logically Collective on SNES 4011 4012 Input Parameters: 4013 + snes - the SNES context 4014 - prefix - the prefix to prepend to all option names 4015 4016 Notes: 4017 A hyphen (-) must NOT be given at the beginning of the prefix name. 4018 The first character of all runtime options is AUTOMATICALLY the hyphen. 4019 4020 Level: advanced 4021 4022 .keywords: SNES, append, options, prefix, database 4023 4024 .seealso: SNESGetOptionsPrefix() 4025 @*/ 4026 PetscErrorCode SNESAppendOptionsPrefix(SNES snes,const char prefix[]) 4027 { 4028 PetscErrorCode ierr; 4029 4030 PetscFunctionBegin; 4031 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4032 ierr = PetscObjectAppendOptionsPrefix((PetscObject)snes,prefix);CHKERRQ(ierr); 4033 if (!snes->ksp) {ierr = SNESGetKSP(snes,&snes->ksp);CHKERRQ(ierr);} 4034 if (snes->linesearch) { 4035 ierr = SNESGetSNESLineSearch(snes,&snes->linesearch);CHKERRQ(ierr); 4036 ierr = PetscObjectAppendOptionsPrefix((PetscObject)snes->linesearch,prefix);CHKERRQ(ierr); 4037 } 4038 ierr = KSPAppendOptionsPrefix(snes->ksp,prefix);CHKERRQ(ierr); 4039 PetscFunctionReturn(0); 4040 } 4041 4042 #undef __FUNCT__ 4043 #define __FUNCT__ "SNESGetOptionsPrefix" 4044 /*@C 4045 SNESGetOptionsPrefix - Sets the prefix used for searching for all 4046 SNES options in the database. 4047 4048 Not Collective 4049 4050 Input Parameter: 4051 . snes - the SNES context 4052 4053 Output Parameter: 4054 . prefix - pointer to the prefix string used 4055 4056 Notes: On the fortran side, the user should pass in a string 'prefix' of 4057 sufficient length to hold the prefix. 4058 4059 Level: advanced 4060 4061 .keywords: SNES, get, options, prefix, database 4062 4063 .seealso: SNESAppendOptionsPrefix() 4064 @*/ 4065 PetscErrorCode SNESGetOptionsPrefix(SNES snes,const char *prefix[]) 4066 { 4067 PetscErrorCode ierr; 4068 4069 PetscFunctionBegin; 4070 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4071 ierr = PetscObjectGetOptionsPrefix((PetscObject)snes,prefix);CHKERRQ(ierr); 4072 PetscFunctionReturn(0); 4073 } 4074 4075 4076 #undef __FUNCT__ 4077 #define __FUNCT__ "SNESRegister" 4078 /*@C 4079 SNESRegister - See SNESRegisterDynamic() 4080 4081 Level: advanced 4082 @*/ 4083 PetscErrorCode SNESRegister(const char sname[],const char path[],const char name[],PetscErrorCode (*function)(SNES)) 4084 { 4085 char fullname[PETSC_MAX_PATH_LEN]; 4086 PetscErrorCode ierr; 4087 4088 PetscFunctionBegin; 4089 ierr = PetscFunctionListConcat(path,name,fullname);CHKERRQ(ierr); 4090 ierr = PetscFunctionListAdd(PETSC_COMM_WORLD,&SNESList,sname,fullname,(void (*)(void))function);CHKERRQ(ierr); 4091 PetscFunctionReturn(0); 4092 } 4093 4094 #undef __FUNCT__ 4095 #define __FUNCT__ "SNESTestLocalMin" 4096 PetscErrorCode SNESTestLocalMin(SNES snes) 4097 { 4098 PetscErrorCode ierr; 4099 PetscInt N,i,j; 4100 Vec u,uh,fh; 4101 PetscScalar value; 4102 PetscReal norm; 4103 4104 PetscFunctionBegin; 4105 ierr = SNESGetSolution(snes,&u);CHKERRQ(ierr); 4106 ierr = VecDuplicate(u,&uh);CHKERRQ(ierr); 4107 ierr = VecDuplicate(u,&fh);CHKERRQ(ierr); 4108 4109 /* currently only works for sequential */ 4110 ierr = PetscPrintf(PETSC_COMM_WORLD,"Testing FormFunction() for local min\n");CHKERRQ(ierr); 4111 ierr = VecGetSize(u,&N);CHKERRQ(ierr); 4112 for (i=0; i<N; i++) { 4113 ierr = VecCopy(u,uh);CHKERRQ(ierr); 4114 ierr = PetscPrintf(PETSC_COMM_WORLD,"i = %D\n",i);CHKERRQ(ierr); 4115 for (j=-10; j<11; j++) { 4116 value = PetscSign(j)*exp(PetscAbs(j)-10.0); 4117 ierr = VecSetValue(uh,i,value,ADD_VALUES);CHKERRQ(ierr); 4118 ierr = SNESComputeFunction(snes,uh,fh);CHKERRQ(ierr); 4119 ierr = VecNorm(fh,NORM_2,&norm);CHKERRQ(ierr); 4120 ierr = PetscPrintf(PETSC_COMM_WORLD," j norm %D %18.16e\n",j,norm);CHKERRQ(ierr); 4121 value = -value; 4122 ierr = VecSetValue(uh,i,value,ADD_VALUES);CHKERRQ(ierr); 4123 } 4124 } 4125 ierr = VecDestroy(&uh);CHKERRQ(ierr); 4126 ierr = VecDestroy(&fh);CHKERRQ(ierr); 4127 PetscFunctionReturn(0); 4128 } 4129 4130 #undef __FUNCT__ 4131 #define __FUNCT__ "SNESKSPSetUseEW" 4132 /*@ 4133 SNESKSPSetUseEW - Sets SNES use Eisenstat-Walker method for 4134 computing relative tolerance for linear solvers within an inexact 4135 Newton method. 4136 4137 Logically Collective on SNES 4138 4139 Input Parameters: 4140 + snes - SNES context 4141 - flag - PETSC_TRUE or PETSC_FALSE 4142 4143 Options Database: 4144 + -snes_ksp_ew - use Eisenstat-Walker method for determining linear system convergence 4145 . -snes_ksp_ew_version ver - version of Eisenstat-Walker method 4146 . -snes_ksp_ew_rtol0 <rtol0> - Sets rtol0 4147 . -snes_ksp_ew_rtolmax <rtolmax> - Sets rtolmax 4148 . -snes_ksp_ew_gamma <gamma> - Sets gamma 4149 . -snes_ksp_ew_alpha <alpha> - Sets alpha 4150 . -snes_ksp_ew_alpha2 <alpha2> - Sets alpha2 4151 - -snes_ksp_ew_threshold <threshold> - Sets threshold 4152 4153 Notes: 4154 Currently, the default is to use a constant relative tolerance for 4155 the inner linear solvers. Alternatively, one can use the 4156 Eisenstat-Walker method, where the relative convergence tolerance 4157 is reset at each Newton iteration according progress of the nonlinear 4158 solver. 4159 4160 Level: advanced 4161 4162 Reference: 4163 S. C. Eisenstat and H. F. Walker, "Choosing the forcing terms in an 4164 inexact Newton method", SISC 17 (1), pp.16-32, 1996. 4165 4166 .keywords: SNES, KSP, Eisenstat, Walker, convergence, test, inexact, Newton 4167 4168 .seealso: SNESKSPGetUseEW(), SNESKSPGetParametersEW(), SNESKSPSetParametersEW() 4169 @*/ 4170 PetscErrorCode SNESKSPSetUseEW(SNES snes,PetscBool flag) 4171 { 4172 PetscFunctionBegin; 4173 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4174 PetscValidLogicalCollectiveBool(snes,flag,2); 4175 snes->ksp_ewconv = flag; 4176 PetscFunctionReturn(0); 4177 } 4178 4179 #undef __FUNCT__ 4180 #define __FUNCT__ "SNESKSPGetUseEW" 4181 /*@ 4182 SNESKSPGetUseEW - Gets if SNES is using Eisenstat-Walker method 4183 for computing relative tolerance for linear solvers within an 4184 inexact Newton method. 4185 4186 Not Collective 4187 4188 Input Parameter: 4189 . snes - SNES context 4190 4191 Output Parameter: 4192 . flag - PETSC_TRUE or PETSC_FALSE 4193 4194 Notes: 4195 Currently, the default is to use a constant relative tolerance for 4196 the inner linear solvers. Alternatively, one can use the 4197 Eisenstat-Walker method, where the relative convergence tolerance 4198 is reset at each Newton iteration according progress of the nonlinear 4199 solver. 4200 4201 Level: advanced 4202 4203 Reference: 4204 S. C. Eisenstat and H. F. Walker, "Choosing the forcing terms in an 4205 inexact Newton method", SISC 17 (1), pp.16-32, 1996. 4206 4207 .keywords: SNES, KSP, Eisenstat, Walker, convergence, test, inexact, Newton 4208 4209 .seealso: SNESKSPSetUseEW(), SNESKSPGetParametersEW(), SNESKSPSetParametersEW() 4210 @*/ 4211 PetscErrorCode SNESKSPGetUseEW(SNES snes, PetscBool *flag) 4212 { 4213 PetscFunctionBegin; 4214 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4215 PetscValidPointer(flag,2); 4216 *flag = snes->ksp_ewconv; 4217 PetscFunctionReturn(0); 4218 } 4219 4220 #undef __FUNCT__ 4221 #define __FUNCT__ "SNESKSPSetParametersEW" 4222 /*@ 4223 SNESKSPSetParametersEW - Sets parameters for Eisenstat-Walker 4224 convergence criteria for the linear solvers within an inexact 4225 Newton method. 4226 4227 Logically Collective on SNES 4228 4229 Input Parameters: 4230 + snes - SNES context 4231 . version - version 1, 2 (default is 2) or 3 4232 . rtol_0 - initial relative tolerance (0 <= rtol_0 < 1) 4233 . rtol_max - maximum relative tolerance (0 <= rtol_max < 1) 4234 . gamma - multiplicative factor for version 2 rtol computation 4235 (0 <= gamma2 <= 1) 4236 . alpha - power for version 2 rtol computation (1 < alpha <= 2) 4237 . alpha2 - power for safeguard 4238 - threshold - threshold for imposing safeguard (0 < threshold < 1) 4239 4240 Note: 4241 Version 3 was contributed by Luis Chacon, June 2006. 4242 4243 Use PETSC_DEFAULT to retain the default for any of the parameters. 4244 4245 Level: advanced 4246 4247 Reference: 4248 S. C. Eisenstat and H. F. Walker, "Choosing the forcing terms in an 4249 inexact Newton method", Utah State University Math. Stat. Dept. Res. 4250 Report 6/94/75, June, 1994, to appear in SIAM J. Sci. Comput. 4251 4252 .keywords: SNES, KSP, Eisenstat, Walker, set, parameters 4253 4254 .seealso: SNESKSPSetUseEW(), SNESKSPGetUseEW(), SNESKSPGetParametersEW() 4255 @*/ 4256 PetscErrorCode SNESKSPSetParametersEW(SNES snes,PetscInt version,PetscReal rtol_0,PetscReal rtol_max,PetscReal gamma,PetscReal alpha,PetscReal alpha2,PetscReal threshold) 4257 { 4258 SNESKSPEW *kctx; 4259 4260 PetscFunctionBegin; 4261 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4262 kctx = (SNESKSPEW*)snes->kspconvctx; 4263 if (!kctx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"No Eisenstat-Walker context existing"); 4264 PetscValidLogicalCollectiveInt(snes,version,2); 4265 PetscValidLogicalCollectiveReal(snes,rtol_0,3); 4266 PetscValidLogicalCollectiveReal(snes,rtol_max,4); 4267 PetscValidLogicalCollectiveReal(snes,gamma,5); 4268 PetscValidLogicalCollectiveReal(snes,alpha,6); 4269 PetscValidLogicalCollectiveReal(snes,alpha2,7); 4270 PetscValidLogicalCollectiveReal(snes,threshold,8); 4271 4272 if (version != PETSC_DEFAULT) kctx->version = version; 4273 if (rtol_0 != PETSC_DEFAULT) kctx->rtol_0 = rtol_0; 4274 if (rtol_max != PETSC_DEFAULT) kctx->rtol_max = rtol_max; 4275 if (gamma != PETSC_DEFAULT) kctx->gamma = gamma; 4276 if (alpha != PETSC_DEFAULT) kctx->alpha = alpha; 4277 if (alpha2 != PETSC_DEFAULT) kctx->alpha2 = alpha2; 4278 if (threshold != PETSC_DEFAULT) kctx->threshold = threshold; 4279 4280 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); 4281 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); 4282 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); 4283 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); 4284 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); 4285 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); 4286 PetscFunctionReturn(0); 4287 } 4288 4289 #undef __FUNCT__ 4290 #define __FUNCT__ "SNESKSPGetParametersEW" 4291 /*@ 4292 SNESKSPGetParametersEW - Gets parameters for Eisenstat-Walker 4293 convergence criteria for the linear solvers within an inexact 4294 Newton method. 4295 4296 Not Collective 4297 4298 Input Parameters: 4299 snes - SNES context 4300 4301 Output Parameters: 4302 + version - version 1, 2 (default is 2) or 3 4303 . rtol_0 - initial relative tolerance (0 <= rtol_0 < 1) 4304 . rtol_max - maximum relative tolerance (0 <= rtol_max < 1) 4305 . gamma - multiplicative factor for version 2 rtol computation (0 <= gamma2 <= 1) 4306 . alpha - power for version 2 rtol computation (1 < alpha <= 2) 4307 . alpha2 - power for safeguard 4308 - threshold - threshold for imposing safeguard (0 < threshold < 1) 4309 4310 Level: advanced 4311 4312 .keywords: SNES, KSP, Eisenstat, Walker, get, parameters 4313 4314 .seealso: SNESKSPSetUseEW(), SNESKSPGetUseEW(), SNESKSPSetParametersEW() 4315 @*/ 4316 PetscErrorCode SNESKSPGetParametersEW(SNES snes,PetscInt *version,PetscReal *rtol_0,PetscReal *rtol_max,PetscReal *gamma,PetscReal *alpha,PetscReal *alpha2,PetscReal *threshold) 4317 { 4318 SNESKSPEW *kctx; 4319 4320 PetscFunctionBegin; 4321 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4322 kctx = (SNESKSPEW*)snes->kspconvctx; 4323 if (!kctx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"No Eisenstat-Walker context existing"); 4324 if (version) *version = kctx->version; 4325 if (rtol_0) *rtol_0 = kctx->rtol_0; 4326 if (rtol_max) *rtol_max = kctx->rtol_max; 4327 if (gamma) *gamma = kctx->gamma; 4328 if (alpha) *alpha = kctx->alpha; 4329 if (alpha2) *alpha2 = kctx->alpha2; 4330 if (threshold) *threshold = kctx->threshold; 4331 PetscFunctionReturn(0); 4332 } 4333 4334 #undef __FUNCT__ 4335 #define __FUNCT__ "SNESKSPEW_PreSolve" 4336 static PetscErrorCode SNESKSPEW_PreSolve(SNES snes, KSP ksp, Vec b, Vec x) 4337 { 4338 PetscErrorCode ierr; 4339 SNESKSPEW *kctx = (SNESKSPEW*)snes->kspconvctx; 4340 PetscReal rtol = PETSC_DEFAULT,stol; 4341 4342 PetscFunctionBegin; 4343 if (!kctx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"No Eisenstat-Walker context exists"); 4344 if (!snes->iter) rtol = kctx->rtol_0; /* first time in, so use the original user rtol */ 4345 else { 4346 if (kctx->version == 1) { 4347 rtol = (snes->norm - kctx->lresid_last)/kctx->norm_last; 4348 if (rtol < 0.0) rtol = -rtol; 4349 stol = pow(kctx->rtol_last,kctx->alpha2); 4350 if (stol > kctx->threshold) rtol = PetscMax(rtol,stol); 4351 } else if (kctx->version == 2) { 4352 rtol = kctx->gamma * pow(snes->norm/kctx->norm_last,kctx->alpha); 4353 stol = kctx->gamma * pow(kctx->rtol_last,kctx->alpha); 4354 if (stol > kctx->threshold) rtol = PetscMax(rtol,stol); 4355 } else if (kctx->version == 3) { /* contributed by Luis Chacon, June 2006. */ 4356 rtol = kctx->gamma * pow(snes->norm/kctx->norm_last,kctx->alpha); 4357 /* safeguard: avoid sharp decrease of rtol */ 4358 stol = kctx->gamma*pow(kctx->rtol_last,kctx->alpha); 4359 stol = PetscMax(rtol,stol); 4360 rtol = PetscMin(kctx->rtol_0,stol); 4361 /* safeguard: avoid oversolving */ 4362 stol = kctx->gamma*(snes->ttol)/snes->norm; 4363 stol = PetscMax(rtol,stol); 4364 rtol = PetscMin(kctx->rtol_0,stol); 4365 } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Only versions 1, 2 or 3 are supported: %D",kctx->version); 4366 } 4367 /* safeguard: avoid rtol greater than one */ 4368 rtol = PetscMin(rtol,kctx->rtol_max); 4369 ierr = KSPSetTolerances(ksp,rtol,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);CHKERRQ(ierr); 4370 ierr = PetscInfo3(snes,"iter %D, Eisenstat-Walker (version %D) KSP rtol=%G\n",snes->iter,kctx->version,rtol);CHKERRQ(ierr); 4371 PetscFunctionReturn(0); 4372 } 4373 4374 #undef __FUNCT__ 4375 #define __FUNCT__ "SNESKSPEW_PostSolve" 4376 static PetscErrorCode SNESKSPEW_PostSolve(SNES snes, KSP ksp, Vec b, Vec x) 4377 { 4378 PetscErrorCode ierr; 4379 SNESKSPEW *kctx = (SNESKSPEW*)snes->kspconvctx; 4380 PCSide pcside; 4381 Vec lres; 4382 4383 PetscFunctionBegin; 4384 if (!kctx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"No Eisenstat-Walker context exists"); 4385 ierr = KSPGetTolerances(ksp,&kctx->rtol_last,0,0,0);CHKERRQ(ierr); 4386 ierr = SNESGetFunctionNorm(snes,&kctx->norm_last);CHKERRQ(ierr); 4387 if (kctx->version == 1) { 4388 ierr = KSPGetPCSide(ksp,&pcside);CHKERRQ(ierr); 4389 if (pcside == PC_RIGHT) { /* XXX Should we also test KSP_UNPRECONDITIONED_NORM ? */ 4390 /* KSP residual is true linear residual */ 4391 ierr = KSPGetResidualNorm(ksp,&kctx->lresid_last);CHKERRQ(ierr); 4392 } else { 4393 /* KSP residual is preconditioned residual */ 4394 /* compute true linear residual norm */ 4395 ierr = VecDuplicate(b,&lres);CHKERRQ(ierr); 4396 ierr = MatMult(snes->jacobian,x,lres);CHKERRQ(ierr); 4397 ierr = VecAYPX(lres,-1.0,b);CHKERRQ(ierr); 4398 ierr = VecNorm(lres,NORM_2,&kctx->lresid_last);CHKERRQ(ierr); 4399 ierr = VecDestroy(&lres);CHKERRQ(ierr); 4400 } 4401 } 4402 PetscFunctionReturn(0); 4403 } 4404 4405 #undef __FUNCT__ 4406 #define __FUNCT__ "SNES_KSPSolve" 4407 PetscErrorCode SNES_KSPSolve(SNES snes, KSP ksp, Vec b, Vec x) 4408 { 4409 PetscErrorCode ierr; 4410 4411 PetscFunctionBegin; 4412 if (snes->ksp_ewconv) { ierr = SNESKSPEW_PreSolve(snes,ksp,b,x);CHKERRQ(ierr); } 4413 ierr = KSPSolve(ksp,b,x);CHKERRQ(ierr); 4414 if (snes->ksp_ewconv) { ierr = SNESKSPEW_PostSolve(snes,ksp,b,x);CHKERRQ(ierr); } 4415 PetscFunctionReturn(0); 4416 } 4417 4418 #include <petsc-private/dmimpl.h> 4419 #undef __FUNCT__ 4420 #define __FUNCT__ "SNESSetDM" 4421 /*@ 4422 SNESSetDM - Sets the DM that may be used by some preconditioners 4423 4424 Logically Collective on SNES 4425 4426 Input Parameters: 4427 + snes - the preconditioner context 4428 - dm - the dm 4429 4430 Level: intermediate 4431 4432 .seealso: SNESGetDM(), KSPSetDM(), KSPGetDM() 4433 @*/ 4434 PetscErrorCode SNESSetDM(SNES snes,DM dm) 4435 { 4436 PetscErrorCode ierr; 4437 KSP ksp; 4438 DMSNES sdm; 4439 4440 PetscFunctionBegin; 4441 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4442 if (dm) {ierr = PetscObjectReference((PetscObject)dm);CHKERRQ(ierr);} 4443 if (snes->dm) { /* Move the DMSNES context over to the new DM unless the new DM already has one */ 4444 if (snes->dm->dmsnes && snes->dmAuto && !dm->dmsnes) { 4445 ierr = DMCopyDMSNES(snes->dm,dm);CHKERRQ(ierr); 4446 ierr = DMGetDMSNES(snes->dm,&sdm);CHKERRQ(ierr); 4447 if (sdm->originaldm == snes->dm) sdm->originaldm = dm; /* Grant write privileges to the replacement DM */ 4448 } 4449 ierr = DMDestroy(&snes->dm);CHKERRQ(ierr); 4450 } 4451 snes->dm = dm; 4452 snes->dmAuto = PETSC_FALSE; 4453 4454 ierr = SNESGetKSP(snes,&ksp);CHKERRQ(ierr); 4455 ierr = KSPSetDM(ksp,dm);CHKERRQ(ierr); 4456 ierr = KSPSetDMActive(ksp,PETSC_FALSE);CHKERRQ(ierr); 4457 if (snes->pc) { 4458 ierr = SNESSetDM(snes->pc, snes->dm);CHKERRQ(ierr); 4459 ierr = SNESSetPCSide(snes,snes->pcside);CHKERRQ(ierr); 4460 } 4461 PetscFunctionReturn(0); 4462 } 4463 4464 #undef __FUNCT__ 4465 #define __FUNCT__ "SNESGetDM" 4466 /*@ 4467 SNESGetDM - Gets the DM that may be used by some preconditioners 4468 4469 Not Collective but DM obtained is parallel on SNES 4470 4471 Input Parameter: 4472 . snes - the preconditioner context 4473 4474 Output Parameter: 4475 . dm - the dm 4476 4477 Level: intermediate 4478 4479 .seealso: SNESSetDM(), KSPSetDM(), KSPGetDM() 4480 @*/ 4481 PetscErrorCode SNESGetDM(SNES snes,DM *dm) 4482 { 4483 PetscErrorCode ierr; 4484 4485 PetscFunctionBegin; 4486 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4487 if (!snes->dm) { 4488 ierr = DMShellCreate(PetscObjectComm((PetscObject)snes),&snes->dm);CHKERRQ(ierr); 4489 snes->dmAuto = PETSC_TRUE; 4490 } 4491 *dm = snes->dm; 4492 PetscFunctionReturn(0); 4493 } 4494 4495 #undef __FUNCT__ 4496 #define __FUNCT__ "SNESSetPC" 4497 /*@ 4498 SNESSetPC - Sets the nonlinear preconditioner to be used. 4499 4500 Collective on SNES 4501 4502 Input Parameters: 4503 + snes - iterative context obtained from SNESCreate() 4504 - pc - the preconditioner object 4505 4506 Notes: 4507 Use SNESGetPC() to retrieve the preconditioner context (for example, 4508 to configure it using the API). 4509 4510 Level: developer 4511 4512 .keywords: SNES, set, precondition 4513 .seealso: SNESGetPC() 4514 @*/ 4515 PetscErrorCode SNESSetPC(SNES snes, SNES pc) 4516 { 4517 PetscErrorCode ierr; 4518 4519 PetscFunctionBegin; 4520 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); 4521 PetscValidHeaderSpecific(pc, SNES_CLASSID, 2); 4522 PetscCheckSameComm(snes, 1, pc, 2); 4523 ierr = PetscObjectReference((PetscObject) pc);CHKERRQ(ierr); 4524 ierr = SNESDestroy(&snes->pc);CHKERRQ(ierr); 4525 snes->pc = pc; 4526 ierr = PetscLogObjectParent(snes, snes->pc);CHKERRQ(ierr); 4527 PetscFunctionReturn(0); 4528 } 4529 4530 #undef __FUNCT__ 4531 #define __FUNCT__ "SNESGetPC" 4532 /*@ 4533 SNESGetPC - Returns a pointer to the nonlinear preconditioning context set with SNESSetPC(). 4534 4535 Not Collective 4536 4537 Input Parameter: 4538 . snes - iterative context obtained from SNESCreate() 4539 4540 Output Parameter: 4541 . pc - preconditioner context 4542 4543 Level: developer 4544 4545 .keywords: SNES, get, preconditioner 4546 .seealso: SNESSetPC() 4547 @*/ 4548 PetscErrorCode SNESGetPC(SNES snes, SNES *pc) 4549 { 4550 PetscErrorCode ierr; 4551 const char *optionsprefix; 4552 4553 PetscFunctionBegin; 4554 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); 4555 PetscValidPointer(pc, 2); 4556 if (!snes->pc) { 4557 ierr = SNESCreate(PetscObjectComm((PetscObject)snes),&snes->pc);CHKERRQ(ierr); 4558 ierr = PetscObjectIncrementTabLevel((PetscObject)snes->pc,(PetscObject)snes,1);CHKERRQ(ierr); 4559 ierr = PetscLogObjectParent(snes,snes->pc);CHKERRQ(ierr); 4560 ierr = SNESGetOptionsPrefix(snes,&optionsprefix);CHKERRQ(ierr); 4561 ierr = SNESSetOptionsPrefix(snes->pc,optionsprefix);CHKERRQ(ierr); 4562 ierr = SNESAppendOptionsPrefix(snes->pc,"npc_");CHKERRQ(ierr); 4563 } 4564 *pc = snes->pc; 4565 PetscFunctionReturn(0); 4566 } 4567 4568 #undef __FUNCT__ 4569 #define __FUNCT__ "SNESSetPCSide" 4570 /*@ 4571 SNESSetPCSide - Sets the preconditioning side. 4572 4573 Logically Collective on SNES 4574 4575 Input Parameter: 4576 . snes - iterative context obtained from SNESCreate() 4577 4578 Output Parameter: 4579 . side - the preconditioning side, where side is one of 4580 .vb 4581 PC_LEFT - left preconditioning (default) 4582 PC_RIGHT - right preconditioning 4583 .ve 4584 4585 Options Database Keys: 4586 . -snes_pc_side <right,left> 4587 4588 Level: intermediate 4589 4590 .keywords: SNES, set, right, left, side, preconditioner, flag 4591 4592 .seealso: SNESGetPCSide(), KSPSetPCSide() 4593 @*/ 4594 PetscErrorCode SNESSetPCSide(SNES snes,PCSide side) 4595 { 4596 PetscFunctionBegin; 4597 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4598 PetscValidLogicalCollectiveEnum(snes,side,2); 4599 snes->pcside = side; 4600 PetscFunctionReturn(0); 4601 } 4602 4603 #undef __FUNCT__ 4604 #define __FUNCT__ "SNESGetPCSide" 4605 /*@ 4606 SNESGetPCSide - Gets the preconditioning side. 4607 4608 Not Collective 4609 4610 Input Parameter: 4611 . snes - iterative context obtained from SNESCreate() 4612 4613 Output Parameter: 4614 . side - the preconditioning side, where side is one of 4615 .vb 4616 PC_LEFT - left preconditioning (default) 4617 PC_RIGHT - right preconditioning 4618 .ve 4619 4620 Level: intermediate 4621 4622 .keywords: SNES, get, right, left, side, preconditioner, flag 4623 4624 .seealso: SNESSetPCSide(), KSPGetPCSide() 4625 @*/ 4626 PetscErrorCode SNESGetPCSide(SNES snes,PCSide *side) 4627 { 4628 PetscFunctionBegin; 4629 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4630 PetscValidPointer(side,2); 4631 *side = snes->pcside; 4632 PetscFunctionReturn(0); 4633 } 4634 4635 #undef __FUNCT__ 4636 #define __FUNCT__ "SNESSetSNESLineSearch" 4637 /*@ 4638 SNESSetSNESLineSearch - Sets the linesearch on the SNES instance. 4639 4640 Collective on SNES 4641 4642 Input Parameters: 4643 + snes - iterative context obtained from SNESCreate() 4644 - linesearch - the linesearch object 4645 4646 Notes: 4647 Use SNESGetSNESLineSearch() to retrieve the preconditioner context (for example, 4648 to configure it using the API). 4649 4650 Level: developer 4651 4652 .keywords: SNES, set, linesearch 4653 .seealso: SNESGetSNESLineSearch() 4654 @*/ 4655 PetscErrorCode SNESSetSNESLineSearch(SNES snes, SNESLineSearch linesearch) 4656 { 4657 PetscErrorCode ierr; 4658 4659 PetscFunctionBegin; 4660 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); 4661 PetscValidHeaderSpecific(linesearch, SNESLINESEARCH_CLASSID, 2); 4662 PetscCheckSameComm(snes, 1, linesearch, 2); 4663 ierr = PetscObjectReference((PetscObject) linesearch);CHKERRQ(ierr); 4664 ierr = SNESLineSearchDestroy(&snes->linesearch);CHKERRQ(ierr); 4665 4666 snes->linesearch = linesearch; 4667 4668 ierr = PetscLogObjectParent(snes, snes->linesearch);CHKERRQ(ierr); 4669 PetscFunctionReturn(0); 4670 } 4671 4672 #undef __FUNCT__ 4673 #define __FUNCT__ "SNESGetSNESLineSearch" 4674 /*@C 4675 SNESGetSNESLineSearch - Returns a pointer to the line search context set with SNESSetLineSearch() 4676 or creates a default line search instance associated with the SNES and returns it. 4677 4678 Not Collective 4679 4680 Input Parameter: 4681 . snes - iterative context obtained from SNESCreate() 4682 4683 Output Parameter: 4684 . linesearch - linesearch context 4685 4686 Level: developer 4687 4688 .keywords: SNES, get, linesearch 4689 .seealso: SNESSetSNESLineSearch() 4690 @*/ 4691 PetscErrorCode SNESGetSNESLineSearch(SNES snes, SNESLineSearch *linesearch) 4692 { 4693 PetscErrorCode ierr; 4694 const char *optionsprefix; 4695 4696 PetscFunctionBegin; 4697 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); 4698 PetscValidPointer(linesearch, 2); 4699 if (!snes->linesearch) { 4700 ierr = SNESGetOptionsPrefix(snes, &optionsprefix);CHKERRQ(ierr); 4701 ierr = SNESLineSearchCreate(PetscObjectComm((PetscObject)snes), &snes->linesearch);CHKERRQ(ierr); 4702 ierr = SNESLineSearchSetSNES(snes->linesearch, snes);CHKERRQ(ierr); 4703 ierr = SNESLineSearchAppendOptionsPrefix(snes->linesearch, optionsprefix);CHKERRQ(ierr); 4704 ierr = PetscObjectIncrementTabLevel((PetscObject) snes->linesearch, (PetscObject) snes, 1);CHKERRQ(ierr); 4705 ierr = PetscLogObjectParent(snes, snes->linesearch);CHKERRQ(ierr); 4706 } 4707 *linesearch = snes->linesearch; 4708 PetscFunctionReturn(0); 4709 } 4710 4711 #if defined(PETSC_HAVE_MATLAB_ENGINE) 4712 #include <mex.h> 4713 4714 typedef struct {char *funcname; mxArray *ctx;} SNESMatlabContext; 4715 4716 #undef __FUNCT__ 4717 #define __FUNCT__ "SNESComputeFunction_Matlab" 4718 /* 4719 SNESComputeFunction_Matlab - Calls the function that has been set with SNESSetFunctionMatlab(). 4720 4721 Collective on SNES 4722 4723 Input Parameters: 4724 + snes - the SNES context 4725 - x - input vector 4726 4727 Output Parameter: 4728 . y - function vector, as set by SNESSetFunction() 4729 4730 Notes: 4731 SNESComputeFunction() is typically used within nonlinear solvers 4732 implementations, so most users would not generally call this routine 4733 themselves. 4734 4735 Level: developer 4736 4737 .keywords: SNES, nonlinear, compute, function 4738 4739 .seealso: SNESSetFunction(), SNESGetFunction() 4740 */ 4741 PetscErrorCode SNESComputeFunction_Matlab(SNES snes,Vec x,Vec y, void *ctx) 4742 { 4743 PetscErrorCode ierr; 4744 SNESMatlabContext *sctx = (SNESMatlabContext*)ctx; 4745 int nlhs = 1,nrhs = 5; 4746 mxArray *plhs[1],*prhs[5]; 4747 long long int lx = 0,ly = 0,ls = 0; 4748 4749 PetscFunctionBegin; 4750 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4751 PetscValidHeaderSpecific(x,VEC_CLASSID,2); 4752 PetscValidHeaderSpecific(y,VEC_CLASSID,3); 4753 PetscCheckSameComm(snes,1,x,2); 4754 PetscCheckSameComm(snes,1,y,3); 4755 4756 /* call Matlab function in ctx with arguments x and y */ 4757 4758 ierr = PetscMemcpy(&ls,&snes,sizeof(snes));CHKERRQ(ierr); 4759 ierr = PetscMemcpy(&lx,&x,sizeof(x));CHKERRQ(ierr); 4760 ierr = PetscMemcpy(&ly,&y,sizeof(x));CHKERRQ(ierr); 4761 prhs[0] = mxCreateDoubleScalar((double)ls); 4762 prhs[1] = mxCreateDoubleScalar((double)lx); 4763 prhs[2] = mxCreateDoubleScalar((double)ly); 4764 prhs[3] = mxCreateString(sctx->funcname); 4765 prhs[4] = sctx->ctx; 4766 ierr = mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscSNESComputeFunctionInternal");CHKERRQ(ierr); 4767 ierr = mxGetScalar(plhs[0]);CHKERRQ(ierr); 4768 mxDestroyArray(prhs[0]); 4769 mxDestroyArray(prhs[1]); 4770 mxDestroyArray(prhs[2]); 4771 mxDestroyArray(prhs[3]); 4772 mxDestroyArray(plhs[0]); 4773 PetscFunctionReturn(0); 4774 } 4775 4776 #undef __FUNCT__ 4777 #define __FUNCT__ "SNESSetFunctionMatlab" 4778 /* 4779 SNESSetFunctionMatlab - Sets the function evaluation routine and function 4780 vector for use by the SNES routines in solving systems of nonlinear 4781 equations from MATLAB. Here the function is a string containing the name of a MATLAB function 4782 4783 Logically Collective on SNES 4784 4785 Input Parameters: 4786 + snes - the SNES context 4787 . r - vector to store function value 4788 - SNESFunction - function evaluation routine 4789 4790 Notes: 4791 The Newton-like methods typically solve linear systems of the form 4792 $ f'(x) x = -f(x), 4793 where f'(x) denotes the Jacobian matrix and f(x) is the function. 4794 4795 Level: beginner 4796 4797 .keywords: SNES, nonlinear, set, function 4798 4799 .seealso: SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESSetFunction() 4800 */ 4801 PetscErrorCode SNESSetFunctionMatlab(SNES snes,Vec r,const char *SNESFunction,mxArray *ctx) 4802 { 4803 PetscErrorCode ierr; 4804 SNESMatlabContext *sctx; 4805 4806 PetscFunctionBegin; 4807 /* currently sctx is memory bleed */ 4808 ierr = PetscMalloc(sizeof(SNESMatlabContext),&sctx);CHKERRQ(ierr); 4809 ierr = PetscStrallocpy(SNESFunction,&sctx->funcname);CHKERRQ(ierr); 4810 /* 4811 This should work, but it doesn't 4812 sctx->ctx = ctx; 4813 mexMakeArrayPersistent(sctx->ctx); 4814 */ 4815 sctx->ctx = mxDuplicateArray(ctx); 4816 ierr = SNESSetFunction(snes,r,SNESComputeFunction_Matlab,sctx);CHKERRQ(ierr); 4817 PetscFunctionReturn(0); 4818 } 4819 4820 #undef __FUNCT__ 4821 #define __FUNCT__ "SNESComputeJacobian_Matlab" 4822 /* 4823 SNESComputeJacobian_Matlab - Calls the function that has been set with SNESSetJacobianMatlab(). 4824 4825 Collective on SNES 4826 4827 Input Parameters: 4828 + snes - the SNES context 4829 . x - input vector 4830 . A, B - the matrices 4831 - ctx - user context 4832 4833 Output Parameter: 4834 . flag - structure of the matrix 4835 4836 Level: developer 4837 4838 .keywords: SNES, nonlinear, compute, function 4839 4840 .seealso: SNESSetFunction(), SNESGetFunction() 4841 @*/ 4842 PetscErrorCode SNESComputeJacobian_Matlab(SNES snes,Vec x,Mat *A,Mat *B,MatStructure *flag, void *ctx) 4843 { 4844 PetscErrorCode ierr; 4845 SNESMatlabContext *sctx = (SNESMatlabContext*)ctx; 4846 int nlhs = 2,nrhs = 6; 4847 mxArray *plhs[2],*prhs[6]; 4848 long long int lx = 0,lA = 0,ls = 0, lB = 0; 4849 4850 PetscFunctionBegin; 4851 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4852 PetscValidHeaderSpecific(x,VEC_CLASSID,2); 4853 4854 /* call Matlab function in ctx with arguments x and y */ 4855 4856 ierr = PetscMemcpy(&ls,&snes,sizeof(snes));CHKERRQ(ierr); 4857 ierr = PetscMemcpy(&lx,&x,sizeof(x));CHKERRQ(ierr); 4858 ierr = PetscMemcpy(&lA,A,sizeof(x));CHKERRQ(ierr); 4859 ierr = PetscMemcpy(&lB,B,sizeof(x));CHKERRQ(ierr); 4860 prhs[0] = mxCreateDoubleScalar((double)ls); 4861 prhs[1] = mxCreateDoubleScalar((double)lx); 4862 prhs[2] = mxCreateDoubleScalar((double)lA); 4863 prhs[3] = mxCreateDoubleScalar((double)lB); 4864 prhs[4] = mxCreateString(sctx->funcname); 4865 prhs[5] = sctx->ctx; 4866 ierr = mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscSNESComputeJacobianInternal");CHKERRQ(ierr); 4867 ierr = mxGetScalar(plhs[0]);CHKERRQ(ierr); 4868 *flag = (MatStructure) mxGetScalar(plhs[1]);CHKERRQ(ierr); 4869 mxDestroyArray(prhs[0]); 4870 mxDestroyArray(prhs[1]); 4871 mxDestroyArray(prhs[2]); 4872 mxDestroyArray(prhs[3]); 4873 mxDestroyArray(prhs[4]); 4874 mxDestroyArray(plhs[0]); 4875 mxDestroyArray(plhs[1]); 4876 PetscFunctionReturn(0); 4877 } 4878 4879 #undef __FUNCT__ 4880 #define __FUNCT__ "SNESSetJacobianMatlab" 4881 /* 4882 SNESSetJacobianMatlab - Sets the Jacobian function evaluation routine and two empty Jacobian matrices 4883 vector for use by the SNES routines in solving systems of nonlinear 4884 equations from MATLAB. Here the function is a string containing the name of a MATLAB function 4885 4886 Logically Collective on SNES 4887 4888 Input Parameters: 4889 + snes - the SNES context 4890 . A,B - Jacobian matrices 4891 . SNESJacobianFunction - function evaluation routine 4892 - ctx - user context 4893 4894 Level: developer 4895 4896 .keywords: SNES, nonlinear, set, function 4897 4898 .seealso: SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESSetFunction(), SNESJacobianFunction 4899 */ 4900 PetscErrorCode SNESSetJacobianMatlab(SNES snes,Mat A,Mat B,const char *SNESJacobianFunction,mxArray *ctx) 4901 { 4902 PetscErrorCode ierr; 4903 SNESMatlabContext *sctx; 4904 4905 PetscFunctionBegin; 4906 /* currently sctx is memory bleed */ 4907 ierr = PetscMalloc(sizeof(SNESMatlabContext),&sctx);CHKERRQ(ierr); 4908 ierr = PetscStrallocpy(SNESJacobianFunction,&sctx->funcname);CHKERRQ(ierr); 4909 /* 4910 This should work, but it doesn't 4911 sctx->ctx = ctx; 4912 mexMakeArrayPersistent(sctx->ctx); 4913 */ 4914 sctx->ctx = mxDuplicateArray(ctx); 4915 ierr = SNESSetJacobian(snes,A,B,SNESComputeJacobian_Matlab,sctx);CHKERRQ(ierr); 4916 PetscFunctionReturn(0); 4917 } 4918 4919 #undef __FUNCT__ 4920 #define __FUNCT__ "SNESMonitor_Matlab" 4921 /* 4922 SNESMonitor_Matlab - Calls the function that has been set with SNESMonitorSetMatlab(). 4923 4924 Collective on SNES 4925 4926 .seealso: SNESSetFunction(), SNESGetFunction() 4927 @*/ 4928 PetscErrorCode SNESMonitor_Matlab(SNES snes,PetscInt it, PetscReal fnorm, void *ctx) 4929 { 4930 PetscErrorCode ierr; 4931 SNESMatlabContext *sctx = (SNESMatlabContext*)ctx; 4932 int nlhs = 1,nrhs = 6; 4933 mxArray *plhs[1],*prhs[6]; 4934 long long int lx = 0,ls = 0; 4935 Vec x = snes->vec_sol; 4936 4937 PetscFunctionBegin; 4938 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4939 4940 ierr = PetscMemcpy(&ls,&snes,sizeof(snes));CHKERRQ(ierr); 4941 ierr = PetscMemcpy(&lx,&x,sizeof(x));CHKERRQ(ierr); 4942 prhs[0] = mxCreateDoubleScalar((double)ls); 4943 prhs[1] = mxCreateDoubleScalar((double)it); 4944 prhs[2] = mxCreateDoubleScalar((double)fnorm); 4945 prhs[3] = mxCreateDoubleScalar((double)lx); 4946 prhs[4] = mxCreateString(sctx->funcname); 4947 prhs[5] = sctx->ctx; 4948 ierr = mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscSNESMonitorInternal");CHKERRQ(ierr); 4949 ierr = mxGetScalar(plhs[0]);CHKERRQ(ierr); 4950 mxDestroyArray(prhs[0]); 4951 mxDestroyArray(prhs[1]); 4952 mxDestroyArray(prhs[2]); 4953 mxDestroyArray(prhs[3]); 4954 mxDestroyArray(prhs[4]); 4955 mxDestroyArray(plhs[0]); 4956 PetscFunctionReturn(0); 4957 } 4958 4959 #undef __FUNCT__ 4960 #define __FUNCT__ "SNESMonitorSetMatlab" 4961 /* 4962 SNESMonitorSetMatlab - Sets the monitor function from MATLAB 4963 4964 Level: developer 4965 4966 .keywords: SNES, nonlinear, set, function 4967 4968 .seealso: SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESSetFunction() 4969 */ 4970 PetscErrorCode SNESMonitorSetMatlab(SNES snes,const char *SNESMonitorFunction,mxArray *ctx) 4971 { 4972 PetscErrorCode ierr; 4973 SNESMatlabContext *sctx; 4974 4975 PetscFunctionBegin; 4976 /* currently sctx is memory bleed */ 4977 ierr = PetscMalloc(sizeof(SNESMatlabContext),&sctx);CHKERRQ(ierr); 4978 ierr = PetscStrallocpy(SNESMonitorFunction,&sctx->funcname);CHKERRQ(ierr); 4979 /* 4980 This should work, but it doesn't 4981 sctx->ctx = ctx; 4982 mexMakeArrayPersistent(sctx->ctx); 4983 */ 4984 sctx->ctx = mxDuplicateArray(ctx); 4985 ierr = SNESMonitorSet(snes,SNESMonitor_Matlab,sctx,NULL);CHKERRQ(ierr); 4986 PetscFunctionReturn(0); 4987 } 4988 4989 #endif 4990