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