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