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 3628 while (1) ; 3629 ierr = PetscLogEventEnd(SNES_Solve,snes,0,0,0);CHKERRQ(ierr); 3630 if (snes->domainerror) { 3631 snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; 3632 snes->domainerror = PETSC_FALSE; 3633 } 3634 if (!snes->reason) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Internal error, solver returned without setting converged reason"); 3635 3636 flg = PETSC_FALSE; 3637 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_test_local_min",&flg,NULL);CHKERRQ(ierr); 3638 if (flg && !PetscPreLoadingOn) { ierr = SNESTestLocalMin(snes);CHKERRQ(ierr); } 3639 if (snes->printreason) { 3640 ierr = PetscViewerASCIIAddTab(PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)snes)),((PetscObject)snes)->tablevel);CHKERRQ(ierr); 3641 if (snes->reason > 0) { 3642 ierr = PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)snes)),"Nonlinear solve converged due to %s iterations %D\n",SNESConvergedReasons[snes->reason],snes->iter);CHKERRQ(ierr); 3643 } else { 3644 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); 3645 } 3646 ierr = PetscViewerASCIISubtractTab(PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)snes)),((PetscObject)snes)->tablevel);CHKERRQ(ierr); 3647 } 3648 3649 if (snes->errorifnotconverged && snes->reason < 0) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_NOT_CONVERGED,"SNESSolve has not converged"); 3650 if (grid < snes->gridsequence) { 3651 DM fine; 3652 Vec xnew; 3653 Mat interp; 3654 3655 ierr = DMRefine(snes->dm,PetscObjectComm((PetscObject)snes),&fine);CHKERRQ(ierr); 3656 if (!fine) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_INCOMP,"DMRefine() did not perform any refinement, cannot continue grid sequencing"); 3657 ierr = DMCreateInterpolation(snes->dm,fine,&interp,NULL);CHKERRQ(ierr); 3658 ierr = DMCreateGlobalVector(fine,&xnew);CHKERRQ(ierr); 3659 ierr = MatInterpolate(interp,x,xnew);CHKERRQ(ierr); 3660 ierr = DMInterpolate(snes->dm,interp,fine);CHKERRQ(ierr); 3661 ierr = MatDestroy(&interp);CHKERRQ(ierr); 3662 x = xnew; 3663 3664 ierr = SNESReset(snes);CHKERRQ(ierr); 3665 ierr = SNESSetDM(snes,fine);CHKERRQ(ierr); 3666 ierr = DMDestroy(&fine);CHKERRQ(ierr); 3667 ierr = PetscViewerASCIIPopTab(PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)snes)));CHKERRQ(ierr); 3668 } 3669 } 3670 /* monitoring and viewing */ 3671 ierr = PetscOptionsGetViewer(PetscObjectComm((PetscObject)snes),((PetscObject)snes)->prefix,"-snes_view",&viewer,&format,&flg);CHKERRQ(ierr); 3672 if (flg && !PetscPreLoadingOn) { 3673 ierr = PetscViewerPushFormat(viewer,format);CHKERRQ(ierr); 3674 ierr = SNESView(snes,viewer);CHKERRQ(ierr); 3675 ierr = PetscViewerPopFormat(viewer);CHKERRQ(ierr); 3676 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 3677 } 3678 ierr = VecViewFromOptions(snes->vec_sol,"-snes_view_solution");CHKERRQ(ierr); 3679 3680 ierr = VecDestroy(&xcreated);CHKERRQ(ierr); 3681 PetscFunctionReturn(0); 3682 } 3683 3684 /* --------- Internal routines for SNES Package --------- */ 3685 3686 #undef __FUNCT__ 3687 #define __FUNCT__ "SNESSetType" 3688 /*@C 3689 SNESSetType - Sets the method for the nonlinear solver. 3690 3691 Collective on SNES 3692 3693 Input Parameters: 3694 + snes - the SNES context 3695 - type - a known method 3696 3697 Options Database Key: 3698 . -snes_type <type> - Sets the method; use -help for a list 3699 of available methods (for instance, newtonls or newtontr) 3700 3701 Notes: 3702 See "petsc/include/petscsnes.h" for available methods (for instance) 3703 + SNESNEWTONLS - Newton's method with line search 3704 (systems of nonlinear equations) 3705 . SNESNEWTONTR - Newton's method with trust region 3706 (systems of nonlinear equations) 3707 3708 Normally, it is best to use the SNESSetFromOptions() command and then 3709 set the SNES solver type from the options database rather than by using 3710 this routine. Using the options database provides the user with 3711 maximum flexibility in evaluating the many nonlinear solvers. 3712 The SNESSetType() routine is provided for those situations where it 3713 is necessary to set the nonlinear solver independently of the command 3714 line or options database. This might be the case, for example, when 3715 the choice of solver changes during the execution of the program, 3716 and the user's application is taking responsibility for choosing the 3717 appropriate method. 3718 3719 Level: intermediate 3720 3721 .keywords: SNES, set, type 3722 3723 .seealso: SNESType, SNESCreate() 3724 3725 @*/ 3726 PetscErrorCode SNESSetType(SNES snes,SNESType type) 3727 { 3728 PetscErrorCode ierr,(*r)(SNES); 3729 PetscBool match; 3730 3731 PetscFunctionBegin; 3732 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3733 PetscValidCharPointer(type,2); 3734 3735 ierr = PetscObjectTypeCompare((PetscObject)snes,type,&match);CHKERRQ(ierr); 3736 if (match) PetscFunctionReturn(0); 3737 3738 ierr = PetscFunctionListFind(PetscObjectComm((PetscObject)snes),SNESList,type,PETSC_TRUE,(void (**)(void)) &r);CHKERRQ(ierr); 3739 if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unable to find requested SNES type %s",type); 3740 /* Destroy the previous private SNES context */ 3741 if (snes->ops->destroy) { 3742 ierr = (*(snes)->ops->destroy)(snes);CHKERRQ(ierr); 3743 snes->ops->destroy = NULL; 3744 } 3745 /* Reinitialize function pointers in SNESOps structure */ 3746 snes->ops->setup = 0; 3747 snes->ops->solve = 0; 3748 snes->ops->view = 0; 3749 snes->ops->setfromoptions = 0; 3750 snes->ops->destroy = 0; 3751 /* Call the SNESCreate_XXX routine for this particular Nonlinear solver */ 3752 snes->setupcalled = PETSC_FALSE; 3753 3754 ierr = PetscObjectChangeTypeName((PetscObject)snes,type);CHKERRQ(ierr); 3755 ierr = (*r)(snes);CHKERRQ(ierr); 3756 #if defined(PETSC_HAVE_AMS) 3757 if (PetscAMSPublishAll) { 3758 ierr = PetscObjectAMSPublish((PetscObject)snes);CHKERRQ(ierr); 3759 } 3760 #endif 3761 PetscFunctionReturn(0); 3762 } 3763 3764 /* --------------------------------------------------------------------- */ 3765 #undef __FUNCT__ 3766 #define __FUNCT__ "SNESRegisterDestroy" 3767 /*@ 3768 SNESRegisterDestroy - Frees the list of nonlinear solvers that were 3769 registered by SNESRegisterDynamic(). 3770 3771 Not Collective 3772 3773 Level: advanced 3774 3775 .keywords: SNES, nonlinear, register, destroy 3776 3777 .seealso: SNESRegisterAll() 3778 @*/ 3779 PetscErrorCode SNESRegisterDestroy(void) 3780 { 3781 PetscErrorCode ierr; 3782 3783 PetscFunctionBegin; 3784 ierr = PetscFunctionListDestroy(&SNESList);CHKERRQ(ierr); 3785 3786 SNESRegisterAllCalled = PETSC_FALSE; 3787 PetscFunctionReturn(0); 3788 } 3789 3790 #undef __FUNCT__ 3791 #define __FUNCT__ "SNESGetType" 3792 /*@C 3793 SNESGetType - Gets the SNES method type and name (as a string). 3794 3795 Not Collective 3796 3797 Input Parameter: 3798 . snes - nonlinear solver context 3799 3800 Output Parameter: 3801 . type - SNES method (a character string) 3802 3803 Level: intermediate 3804 3805 .keywords: SNES, nonlinear, get, type, name 3806 @*/ 3807 PetscErrorCode SNESGetType(SNES snes,SNESType *type) 3808 { 3809 PetscFunctionBegin; 3810 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3811 PetscValidPointer(type,2); 3812 *type = ((PetscObject)snes)->type_name; 3813 PetscFunctionReturn(0); 3814 } 3815 3816 #undef __FUNCT__ 3817 #define __FUNCT__ "SNESGetSolution" 3818 /*@ 3819 SNESGetSolution - Returns the vector where the approximate solution is 3820 stored. This is the fine grid solution when using SNESSetGridSequence(). 3821 3822 Not Collective, but Vec is parallel if SNES is parallel 3823 3824 Input Parameter: 3825 . snes - the SNES context 3826 3827 Output Parameter: 3828 . x - the solution 3829 3830 Level: intermediate 3831 3832 .keywords: SNES, nonlinear, get, solution 3833 3834 .seealso: SNESGetSolutionUpdate(), SNESGetFunction() 3835 @*/ 3836 PetscErrorCode SNESGetSolution(SNES snes,Vec *x) 3837 { 3838 PetscFunctionBegin; 3839 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3840 PetscValidPointer(x,2); 3841 *x = snes->vec_sol; 3842 PetscFunctionReturn(0); 3843 } 3844 3845 #undef __FUNCT__ 3846 #define __FUNCT__ "SNESGetSolutionUpdate" 3847 /*@ 3848 SNESGetSolutionUpdate - Returns the vector where the solution update is 3849 stored. 3850 3851 Not Collective, but Vec is parallel if SNES is parallel 3852 3853 Input Parameter: 3854 . snes - the SNES context 3855 3856 Output Parameter: 3857 . x - the solution update 3858 3859 Level: advanced 3860 3861 .keywords: SNES, nonlinear, get, solution, update 3862 3863 .seealso: SNESGetSolution(), SNESGetFunction() 3864 @*/ 3865 PetscErrorCode SNESGetSolutionUpdate(SNES snes,Vec *x) 3866 { 3867 PetscFunctionBegin; 3868 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3869 PetscValidPointer(x,2); 3870 *x = snes->vec_sol_update; 3871 PetscFunctionReturn(0); 3872 } 3873 3874 #undef __FUNCT__ 3875 #define __FUNCT__ "SNESGetFunction" 3876 /*@C 3877 SNESGetFunction - Returns the vector where the function is stored. 3878 3879 Not Collective, but Vec is parallel if SNES is parallel. Collective if Vec is requested, but has not been created yet. 3880 3881 Input Parameter: 3882 . snes - the SNES context 3883 3884 Output Parameter: 3885 + r - the vector that is used to store residuals (or NULL if you don't want it) 3886 . SNESFunction- the function (or NULL if you don't want it) 3887 - ctx - the function context (or NULL if you don't want it) 3888 3889 Level: advanced 3890 3891 .keywords: SNES, nonlinear, get, function 3892 3893 .seealso: SNESSetFunction(), SNESGetSolution(), SNESFunction 3894 @*/ 3895 PetscErrorCode SNESGetFunction(SNES snes,Vec *r,PetscErrorCode (**SNESFunction)(SNES,Vec,Vec,void*),void **ctx) 3896 { 3897 PetscErrorCode ierr; 3898 DM dm; 3899 3900 PetscFunctionBegin; 3901 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3902 if (r) { 3903 if (!snes->vec_func) { 3904 if (snes->vec_rhs) { 3905 ierr = VecDuplicate(snes->vec_rhs,&snes->vec_func);CHKERRQ(ierr); 3906 } else if (snes->vec_sol) { 3907 ierr = VecDuplicate(snes->vec_sol,&snes->vec_func);CHKERRQ(ierr); 3908 } else if (snes->dm) { 3909 ierr = DMCreateGlobalVector(snes->dm,&snes->vec_func);CHKERRQ(ierr); 3910 } 3911 } 3912 *r = snes->vec_func; 3913 } 3914 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 3915 ierr = DMSNESGetFunction(dm,SNESFunction,ctx);CHKERRQ(ierr); 3916 PetscFunctionReturn(0); 3917 } 3918 3919 /*@C 3920 SNESGetGS - Returns the GS function and context. 3921 3922 Input Parameter: 3923 . snes - the SNES context 3924 3925 Output Parameter: 3926 + SNESGSFunction - the function (or NULL) 3927 - ctx - the function context (or NULL) 3928 3929 Level: advanced 3930 3931 .keywords: SNES, nonlinear, get, function 3932 3933 .seealso: SNESSetGS(), SNESGetFunction() 3934 @*/ 3935 3936 #undef __FUNCT__ 3937 #define __FUNCT__ "SNESGetGS" 3938 PetscErrorCode SNESGetGS (SNES snes, PetscErrorCode (**SNESGSFunction)(SNES, Vec, Vec, void*), void ** ctx) 3939 { 3940 PetscErrorCode ierr; 3941 DM dm; 3942 3943 PetscFunctionBegin; 3944 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3945 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 3946 ierr = DMSNESGetGS(dm,SNESGSFunction,ctx);CHKERRQ(ierr); 3947 PetscFunctionReturn(0); 3948 } 3949 3950 #undef __FUNCT__ 3951 #define __FUNCT__ "SNESSetOptionsPrefix" 3952 /*@C 3953 SNESSetOptionsPrefix - Sets the prefix used for searching for all 3954 SNES options in the database. 3955 3956 Logically Collective on SNES 3957 3958 Input Parameter: 3959 + snes - the SNES context 3960 - prefix - the prefix to prepend to all option names 3961 3962 Notes: 3963 A hyphen (-) must NOT be given at the beginning of the prefix name. 3964 The first character of all runtime options is AUTOMATICALLY the hyphen. 3965 3966 Level: advanced 3967 3968 .keywords: SNES, set, options, prefix, database 3969 3970 .seealso: SNESSetFromOptions() 3971 @*/ 3972 PetscErrorCode SNESSetOptionsPrefix(SNES snes,const char prefix[]) 3973 { 3974 PetscErrorCode ierr; 3975 3976 PetscFunctionBegin; 3977 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3978 ierr = PetscObjectSetOptionsPrefix((PetscObject)snes,prefix);CHKERRQ(ierr); 3979 if (!snes->ksp) {ierr = SNESGetKSP(snes,&snes->ksp);CHKERRQ(ierr);} 3980 if (snes->linesearch) { 3981 ierr = SNESGetSNESLineSearch(snes,&snes->linesearch);CHKERRQ(ierr); 3982 ierr = PetscObjectSetOptionsPrefix((PetscObject)snes->linesearch,prefix);CHKERRQ(ierr); 3983 } 3984 ierr = KSPSetOptionsPrefix(snes->ksp,prefix);CHKERRQ(ierr); 3985 PetscFunctionReturn(0); 3986 } 3987 3988 #undef __FUNCT__ 3989 #define __FUNCT__ "SNESAppendOptionsPrefix" 3990 /*@C 3991 SNESAppendOptionsPrefix - Appends to the prefix used for searching for all 3992 SNES options in the database. 3993 3994 Logically Collective on SNES 3995 3996 Input Parameters: 3997 + snes - the SNES context 3998 - prefix - the prefix to prepend to all option names 3999 4000 Notes: 4001 A hyphen (-) must NOT be given at the beginning of the prefix name. 4002 The first character of all runtime options is AUTOMATICALLY the hyphen. 4003 4004 Level: advanced 4005 4006 .keywords: SNES, append, options, prefix, database 4007 4008 .seealso: SNESGetOptionsPrefix() 4009 @*/ 4010 PetscErrorCode SNESAppendOptionsPrefix(SNES snes,const char prefix[]) 4011 { 4012 PetscErrorCode ierr; 4013 4014 PetscFunctionBegin; 4015 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4016 ierr = PetscObjectAppendOptionsPrefix((PetscObject)snes,prefix);CHKERRQ(ierr); 4017 if (!snes->ksp) {ierr = SNESGetKSP(snes,&snes->ksp);CHKERRQ(ierr);} 4018 if (snes->linesearch) { 4019 ierr = SNESGetSNESLineSearch(snes,&snes->linesearch);CHKERRQ(ierr); 4020 ierr = PetscObjectAppendOptionsPrefix((PetscObject)snes->linesearch,prefix);CHKERRQ(ierr); 4021 } 4022 ierr = KSPAppendOptionsPrefix(snes->ksp,prefix);CHKERRQ(ierr); 4023 PetscFunctionReturn(0); 4024 } 4025 4026 #undef __FUNCT__ 4027 #define __FUNCT__ "SNESGetOptionsPrefix" 4028 /*@C 4029 SNESGetOptionsPrefix - Sets the prefix used for searching for all 4030 SNES options in the database. 4031 4032 Not Collective 4033 4034 Input Parameter: 4035 . snes - the SNES context 4036 4037 Output Parameter: 4038 . prefix - pointer to the prefix string used 4039 4040 Notes: On the fortran side, the user should pass in a string 'prefix' of 4041 sufficient length to hold the prefix. 4042 4043 Level: advanced 4044 4045 .keywords: SNES, get, options, prefix, database 4046 4047 .seealso: SNESAppendOptionsPrefix() 4048 @*/ 4049 PetscErrorCode SNESGetOptionsPrefix(SNES snes,const char *prefix[]) 4050 { 4051 PetscErrorCode ierr; 4052 4053 PetscFunctionBegin; 4054 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4055 ierr = PetscObjectGetOptionsPrefix((PetscObject)snes,prefix);CHKERRQ(ierr); 4056 PetscFunctionReturn(0); 4057 } 4058 4059 4060 #undef __FUNCT__ 4061 #define __FUNCT__ "SNESRegister" 4062 /*@C 4063 SNESRegister - See SNESRegisterDynamic() 4064 4065 Level: advanced 4066 @*/ 4067 PetscErrorCode SNESRegister(const char sname[],const char path[],const char name[],PetscErrorCode (*function)(SNES)) 4068 { 4069 char fullname[PETSC_MAX_PATH_LEN]; 4070 PetscErrorCode ierr; 4071 4072 PetscFunctionBegin; 4073 ierr = PetscFunctionListConcat(path,name,fullname);CHKERRQ(ierr); 4074 ierr = PetscFunctionListAdd(PETSC_COMM_WORLD,&SNESList,sname,fullname,(void (*)(void))function);CHKERRQ(ierr); 4075 PetscFunctionReturn(0); 4076 } 4077 4078 #undef __FUNCT__ 4079 #define __FUNCT__ "SNESTestLocalMin" 4080 PetscErrorCode SNESTestLocalMin(SNES snes) 4081 { 4082 PetscErrorCode ierr; 4083 PetscInt N,i,j; 4084 Vec u,uh,fh; 4085 PetscScalar value; 4086 PetscReal norm; 4087 4088 PetscFunctionBegin; 4089 ierr = SNESGetSolution(snes,&u);CHKERRQ(ierr); 4090 ierr = VecDuplicate(u,&uh);CHKERRQ(ierr); 4091 ierr = VecDuplicate(u,&fh);CHKERRQ(ierr); 4092 4093 /* currently only works for sequential */ 4094 ierr = PetscPrintf(PETSC_COMM_WORLD,"Testing FormFunction() for local min\n");CHKERRQ(ierr); 4095 ierr = VecGetSize(u,&N);CHKERRQ(ierr); 4096 for (i=0; i<N; i++) { 4097 ierr = VecCopy(u,uh);CHKERRQ(ierr); 4098 ierr = PetscPrintf(PETSC_COMM_WORLD,"i = %D\n",i);CHKERRQ(ierr); 4099 for (j=-10; j<11; j++) { 4100 value = PetscSign(j)*exp(PetscAbs(j)-10.0); 4101 ierr = VecSetValue(uh,i,value,ADD_VALUES);CHKERRQ(ierr); 4102 ierr = SNESComputeFunction(snes,uh,fh);CHKERRQ(ierr); 4103 ierr = VecNorm(fh,NORM_2,&norm);CHKERRQ(ierr); 4104 ierr = PetscPrintf(PETSC_COMM_WORLD," j norm %D %18.16e\n",j,norm);CHKERRQ(ierr); 4105 value = -value; 4106 ierr = VecSetValue(uh,i,value,ADD_VALUES);CHKERRQ(ierr); 4107 } 4108 } 4109 ierr = VecDestroy(&uh);CHKERRQ(ierr); 4110 ierr = VecDestroy(&fh);CHKERRQ(ierr); 4111 PetscFunctionReturn(0); 4112 } 4113 4114 #undef __FUNCT__ 4115 #define __FUNCT__ "SNESKSPSetUseEW" 4116 /*@ 4117 SNESKSPSetUseEW - Sets SNES use Eisenstat-Walker method for 4118 computing relative tolerance for linear solvers within an inexact 4119 Newton method. 4120 4121 Logically Collective on SNES 4122 4123 Input Parameters: 4124 + snes - SNES context 4125 - flag - PETSC_TRUE or PETSC_FALSE 4126 4127 Options Database: 4128 + -snes_ksp_ew - use Eisenstat-Walker method for determining linear system convergence 4129 . -snes_ksp_ew_version ver - version of Eisenstat-Walker method 4130 . -snes_ksp_ew_rtol0 <rtol0> - Sets rtol0 4131 . -snes_ksp_ew_rtolmax <rtolmax> - Sets rtolmax 4132 . -snes_ksp_ew_gamma <gamma> - Sets gamma 4133 . -snes_ksp_ew_alpha <alpha> - Sets alpha 4134 . -snes_ksp_ew_alpha2 <alpha2> - Sets alpha2 4135 - -snes_ksp_ew_threshold <threshold> - Sets threshold 4136 4137 Notes: 4138 Currently, the default is to use a constant relative tolerance for 4139 the inner linear solvers. Alternatively, one can use the 4140 Eisenstat-Walker method, where the relative convergence tolerance 4141 is reset at each Newton iteration according progress of the nonlinear 4142 solver. 4143 4144 Level: advanced 4145 4146 Reference: 4147 S. C. Eisenstat and H. F. Walker, "Choosing the forcing terms in an 4148 inexact Newton method", SISC 17 (1), pp.16-32, 1996. 4149 4150 .keywords: SNES, KSP, Eisenstat, Walker, convergence, test, inexact, Newton 4151 4152 .seealso: SNESKSPGetUseEW(), SNESKSPGetParametersEW(), SNESKSPSetParametersEW() 4153 @*/ 4154 PetscErrorCode SNESKSPSetUseEW(SNES snes,PetscBool flag) 4155 { 4156 PetscFunctionBegin; 4157 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4158 PetscValidLogicalCollectiveBool(snes,flag,2); 4159 snes->ksp_ewconv = flag; 4160 PetscFunctionReturn(0); 4161 } 4162 4163 #undef __FUNCT__ 4164 #define __FUNCT__ "SNESKSPGetUseEW" 4165 /*@ 4166 SNESKSPGetUseEW - Gets if SNES is using Eisenstat-Walker method 4167 for computing relative tolerance for linear solvers within an 4168 inexact Newton method. 4169 4170 Not Collective 4171 4172 Input Parameter: 4173 . snes - SNES context 4174 4175 Output Parameter: 4176 . flag - PETSC_TRUE or PETSC_FALSE 4177 4178 Notes: 4179 Currently, the default is to use a constant relative tolerance for 4180 the inner linear solvers. Alternatively, one can use the 4181 Eisenstat-Walker method, where the relative convergence tolerance 4182 is reset at each Newton iteration according progress of the nonlinear 4183 solver. 4184 4185 Level: advanced 4186 4187 Reference: 4188 S. C. Eisenstat and H. F. Walker, "Choosing the forcing terms in an 4189 inexact Newton method", SISC 17 (1), pp.16-32, 1996. 4190 4191 .keywords: SNES, KSP, Eisenstat, Walker, convergence, test, inexact, Newton 4192 4193 .seealso: SNESKSPSetUseEW(), SNESKSPGetParametersEW(), SNESKSPSetParametersEW() 4194 @*/ 4195 PetscErrorCode SNESKSPGetUseEW(SNES snes, PetscBool *flag) 4196 { 4197 PetscFunctionBegin; 4198 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4199 PetscValidPointer(flag,2); 4200 *flag = snes->ksp_ewconv; 4201 PetscFunctionReturn(0); 4202 } 4203 4204 #undef __FUNCT__ 4205 #define __FUNCT__ "SNESKSPSetParametersEW" 4206 /*@ 4207 SNESKSPSetParametersEW - Sets parameters for Eisenstat-Walker 4208 convergence criteria for the linear solvers within an inexact 4209 Newton method. 4210 4211 Logically Collective on SNES 4212 4213 Input Parameters: 4214 + snes - SNES context 4215 . version - version 1, 2 (default is 2) or 3 4216 . rtol_0 - initial relative tolerance (0 <= rtol_0 < 1) 4217 . rtol_max - maximum relative tolerance (0 <= rtol_max < 1) 4218 . gamma - multiplicative factor for version 2 rtol computation 4219 (0 <= gamma2 <= 1) 4220 . alpha - power for version 2 rtol computation (1 < alpha <= 2) 4221 . alpha2 - power for safeguard 4222 - threshold - threshold for imposing safeguard (0 < threshold < 1) 4223 4224 Note: 4225 Version 3 was contributed by Luis Chacon, June 2006. 4226 4227 Use PETSC_DEFAULT to retain the default for any of the parameters. 4228 4229 Level: advanced 4230 4231 Reference: 4232 S. C. Eisenstat and H. F. Walker, "Choosing the forcing terms in an 4233 inexact Newton method", Utah State University Math. Stat. Dept. Res. 4234 Report 6/94/75, June, 1994, to appear in SIAM J. Sci. Comput. 4235 4236 .keywords: SNES, KSP, Eisenstat, Walker, set, parameters 4237 4238 .seealso: SNESKSPSetUseEW(), SNESKSPGetUseEW(), SNESKSPGetParametersEW() 4239 @*/ 4240 PetscErrorCode SNESKSPSetParametersEW(SNES snes,PetscInt version,PetscReal rtol_0,PetscReal rtol_max,PetscReal gamma,PetscReal alpha,PetscReal alpha2,PetscReal threshold) 4241 { 4242 SNESKSPEW *kctx; 4243 4244 PetscFunctionBegin; 4245 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4246 kctx = (SNESKSPEW*)snes->kspconvctx; 4247 if (!kctx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"No Eisenstat-Walker context existing"); 4248 PetscValidLogicalCollectiveInt(snes,version,2); 4249 PetscValidLogicalCollectiveReal(snes,rtol_0,3); 4250 PetscValidLogicalCollectiveReal(snes,rtol_max,4); 4251 PetscValidLogicalCollectiveReal(snes,gamma,5); 4252 PetscValidLogicalCollectiveReal(snes,alpha,6); 4253 PetscValidLogicalCollectiveReal(snes,alpha2,7); 4254 PetscValidLogicalCollectiveReal(snes,threshold,8); 4255 4256 if (version != PETSC_DEFAULT) kctx->version = version; 4257 if (rtol_0 != PETSC_DEFAULT) kctx->rtol_0 = rtol_0; 4258 if (rtol_max != PETSC_DEFAULT) kctx->rtol_max = rtol_max; 4259 if (gamma != PETSC_DEFAULT) kctx->gamma = gamma; 4260 if (alpha != PETSC_DEFAULT) kctx->alpha = alpha; 4261 if (alpha2 != PETSC_DEFAULT) kctx->alpha2 = alpha2; 4262 if (threshold != PETSC_DEFAULT) kctx->threshold = threshold; 4263 4264 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); 4265 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); 4266 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); 4267 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); 4268 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); 4269 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); 4270 PetscFunctionReturn(0); 4271 } 4272 4273 #undef __FUNCT__ 4274 #define __FUNCT__ "SNESKSPGetParametersEW" 4275 /*@ 4276 SNESKSPGetParametersEW - Gets parameters for Eisenstat-Walker 4277 convergence criteria for the linear solvers within an inexact 4278 Newton method. 4279 4280 Not Collective 4281 4282 Input Parameters: 4283 snes - SNES context 4284 4285 Output Parameters: 4286 + version - version 1, 2 (default is 2) or 3 4287 . rtol_0 - initial relative tolerance (0 <= rtol_0 < 1) 4288 . rtol_max - maximum relative tolerance (0 <= rtol_max < 1) 4289 . gamma - multiplicative factor for version 2 rtol computation (0 <= gamma2 <= 1) 4290 . alpha - power for version 2 rtol computation (1 < alpha <= 2) 4291 . alpha2 - power for safeguard 4292 - threshold - threshold for imposing safeguard (0 < threshold < 1) 4293 4294 Level: advanced 4295 4296 .keywords: SNES, KSP, Eisenstat, Walker, get, parameters 4297 4298 .seealso: SNESKSPSetUseEW(), SNESKSPGetUseEW(), SNESKSPSetParametersEW() 4299 @*/ 4300 PetscErrorCode SNESKSPGetParametersEW(SNES snes,PetscInt *version,PetscReal *rtol_0,PetscReal *rtol_max,PetscReal *gamma,PetscReal *alpha,PetscReal *alpha2,PetscReal *threshold) 4301 { 4302 SNESKSPEW *kctx; 4303 4304 PetscFunctionBegin; 4305 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4306 kctx = (SNESKSPEW*)snes->kspconvctx; 4307 if (!kctx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"No Eisenstat-Walker context existing"); 4308 if (version) *version = kctx->version; 4309 if (rtol_0) *rtol_0 = kctx->rtol_0; 4310 if (rtol_max) *rtol_max = kctx->rtol_max; 4311 if (gamma) *gamma = kctx->gamma; 4312 if (alpha) *alpha = kctx->alpha; 4313 if (alpha2) *alpha2 = kctx->alpha2; 4314 if (threshold) *threshold = kctx->threshold; 4315 PetscFunctionReturn(0); 4316 } 4317 4318 #undef __FUNCT__ 4319 #define __FUNCT__ "SNESKSPEW_PreSolve" 4320 static PetscErrorCode SNESKSPEW_PreSolve(SNES snes, KSP ksp, Vec b, Vec x) 4321 { 4322 PetscErrorCode ierr; 4323 SNESKSPEW *kctx = (SNESKSPEW*)snes->kspconvctx; 4324 PetscReal rtol = PETSC_DEFAULT,stol; 4325 4326 PetscFunctionBegin; 4327 if (!kctx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"No Eisenstat-Walker context exists"); 4328 if (!snes->iter) rtol = kctx->rtol_0; /* first time in, so use the original user rtol */ 4329 else { 4330 if (kctx->version == 1) { 4331 rtol = (snes->norm - kctx->lresid_last)/kctx->norm_last; 4332 if (rtol < 0.0) rtol = -rtol; 4333 stol = PetscPowReal(kctx->rtol_last,kctx->alpha2); 4334 if (stol > kctx->threshold) rtol = PetscMax(rtol,stol); 4335 } else if (kctx->version == 2) { 4336 rtol = kctx->gamma * PetscPowReal(snes->norm/kctx->norm_last,kctx->alpha); 4337 stol = kctx->gamma * PetscPowReal(kctx->rtol_last,kctx->alpha); 4338 if (stol > kctx->threshold) rtol = PetscMax(rtol,stol); 4339 } else if (kctx->version == 3) { /* contributed by Luis Chacon, June 2006. */ 4340 rtol = kctx->gamma * PetscPowReal(snes->norm/kctx->norm_last,kctx->alpha); 4341 /* safeguard: avoid sharp decrease of rtol */ 4342 stol = kctx->gamma*PetscPowReal(kctx->rtol_last,kctx->alpha); 4343 stol = PetscMax(rtol,stol); 4344 rtol = PetscMin(kctx->rtol_0,stol); 4345 /* safeguard: avoid oversolving */ 4346 stol = kctx->gamma*(snes->ttol)/snes->norm; 4347 stol = PetscMax(rtol,stol); 4348 rtol = PetscMin(kctx->rtol_0,stol); 4349 } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Only versions 1, 2 or 3 are supported: %D",kctx->version); 4350 } 4351 /* safeguard: avoid rtol greater than one */ 4352 rtol = PetscMin(rtol,kctx->rtol_max); 4353 ierr = KSPSetTolerances(ksp,rtol,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);CHKERRQ(ierr); 4354 ierr = PetscInfo3(snes,"iter %D, Eisenstat-Walker (version %D) KSP rtol=%G\n",snes->iter,kctx->version,rtol);CHKERRQ(ierr); 4355 PetscFunctionReturn(0); 4356 } 4357 4358 #undef __FUNCT__ 4359 #define __FUNCT__ "SNESKSPEW_PostSolve" 4360 static PetscErrorCode SNESKSPEW_PostSolve(SNES snes, KSP ksp, Vec b, Vec x) 4361 { 4362 PetscErrorCode ierr; 4363 SNESKSPEW *kctx = (SNESKSPEW*)snes->kspconvctx; 4364 PCSide pcside; 4365 Vec lres; 4366 4367 PetscFunctionBegin; 4368 if (!kctx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"No Eisenstat-Walker context exists"); 4369 ierr = KSPGetTolerances(ksp,&kctx->rtol_last,0,0,0);CHKERRQ(ierr); 4370 ierr = SNESGetFunctionNorm(snes,&kctx->norm_last);CHKERRQ(ierr); 4371 if (kctx->version == 1) { 4372 ierr = KSPGetPCSide(ksp,&pcside);CHKERRQ(ierr); 4373 if (pcside == PC_RIGHT) { /* XXX Should we also test KSP_UNPRECONDITIONED_NORM ? */ 4374 /* KSP residual is true linear residual */ 4375 ierr = KSPGetResidualNorm(ksp,&kctx->lresid_last);CHKERRQ(ierr); 4376 } else { 4377 /* KSP residual is preconditioned residual */ 4378 /* compute true linear residual norm */ 4379 ierr = VecDuplicate(b,&lres);CHKERRQ(ierr); 4380 ierr = MatMult(snes->jacobian,x,lres);CHKERRQ(ierr); 4381 ierr = VecAYPX(lres,-1.0,b);CHKERRQ(ierr); 4382 ierr = VecNorm(lres,NORM_2,&kctx->lresid_last);CHKERRQ(ierr); 4383 ierr = VecDestroy(&lres);CHKERRQ(ierr); 4384 } 4385 } 4386 PetscFunctionReturn(0); 4387 } 4388 4389 #undef __FUNCT__ 4390 #define __FUNCT__ "SNES_KSPSolve" 4391 PetscErrorCode SNES_KSPSolve(SNES snes, KSP ksp, Vec b, Vec x) 4392 { 4393 PetscErrorCode ierr; 4394 4395 PetscFunctionBegin; 4396 if (snes->ksp_ewconv) { ierr = SNESKSPEW_PreSolve(snes,ksp,b,x);CHKERRQ(ierr); } 4397 ierr = KSPSolve(ksp,b,x);CHKERRQ(ierr); 4398 if (snes->ksp_ewconv) { ierr = SNESKSPEW_PostSolve(snes,ksp,b,x);CHKERRQ(ierr); } 4399 PetscFunctionReturn(0); 4400 } 4401 4402 #include <petsc-private/dmimpl.h> 4403 #undef __FUNCT__ 4404 #define __FUNCT__ "SNESSetDM" 4405 /*@ 4406 SNESSetDM - Sets the DM that may be used by some preconditioners 4407 4408 Logically Collective on SNES 4409 4410 Input Parameters: 4411 + snes - the preconditioner context 4412 - dm - the dm 4413 4414 Level: intermediate 4415 4416 .seealso: SNESGetDM(), KSPSetDM(), KSPGetDM() 4417 @*/ 4418 PetscErrorCode SNESSetDM(SNES snes,DM dm) 4419 { 4420 PetscErrorCode ierr; 4421 KSP ksp; 4422 DMSNES sdm; 4423 4424 PetscFunctionBegin; 4425 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4426 if (dm) {ierr = PetscObjectReference((PetscObject)dm);CHKERRQ(ierr);} 4427 if (snes->dm) { /* Move the DMSNES context over to the new DM unless the new DM already has one */ 4428 if (snes->dm->dmsnes && snes->dmAuto && !dm->dmsnes) { 4429 ierr = DMCopyDMSNES(snes->dm,dm);CHKERRQ(ierr); 4430 ierr = DMGetDMSNES(snes->dm,&sdm);CHKERRQ(ierr); 4431 if (sdm->originaldm == snes->dm) sdm->originaldm = dm; /* Grant write privileges to the replacement DM */ 4432 } 4433 ierr = DMDestroy(&snes->dm);CHKERRQ(ierr); 4434 } 4435 snes->dm = dm; 4436 snes->dmAuto = PETSC_FALSE; 4437 4438 ierr = SNESGetKSP(snes,&ksp);CHKERRQ(ierr); 4439 ierr = KSPSetDM(ksp,dm);CHKERRQ(ierr); 4440 ierr = KSPSetDMActive(ksp,PETSC_FALSE);CHKERRQ(ierr); 4441 if (snes->pc) { 4442 ierr = SNESSetDM(snes->pc, snes->dm);CHKERRQ(ierr); 4443 ierr = SNESSetPCSide(snes,snes->pcside);CHKERRQ(ierr); 4444 } 4445 PetscFunctionReturn(0); 4446 } 4447 4448 #undef __FUNCT__ 4449 #define __FUNCT__ "SNESGetDM" 4450 /*@ 4451 SNESGetDM - Gets the DM that may be used by some preconditioners 4452 4453 Not Collective but DM obtained is parallel on SNES 4454 4455 Input Parameter: 4456 . snes - the preconditioner context 4457 4458 Output Parameter: 4459 . dm - the dm 4460 4461 Level: intermediate 4462 4463 .seealso: SNESSetDM(), KSPSetDM(), KSPGetDM() 4464 @*/ 4465 PetscErrorCode SNESGetDM(SNES snes,DM *dm) 4466 { 4467 PetscErrorCode ierr; 4468 4469 PetscFunctionBegin; 4470 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4471 if (!snes->dm) { 4472 ierr = DMShellCreate(PetscObjectComm((PetscObject)snes),&snes->dm);CHKERRQ(ierr); 4473 snes->dmAuto = PETSC_TRUE; 4474 } 4475 *dm = snes->dm; 4476 PetscFunctionReturn(0); 4477 } 4478 4479 #undef __FUNCT__ 4480 #define __FUNCT__ "SNESSetPC" 4481 /*@ 4482 SNESSetPC - Sets the nonlinear preconditioner to be used. 4483 4484 Collective on SNES 4485 4486 Input Parameters: 4487 + snes - iterative context obtained from SNESCreate() 4488 - pc - the preconditioner object 4489 4490 Notes: 4491 Use SNESGetPC() to retrieve the preconditioner context (for example, 4492 to configure it using the API). 4493 4494 Level: developer 4495 4496 .keywords: SNES, set, precondition 4497 .seealso: SNESGetPC() 4498 @*/ 4499 PetscErrorCode SNESSetPC(SNES snes, SNES pc) 4500 { 4501 PetscErrorCode ierr; 4502 4503 PetscFunctionBegin; 4504 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); 4505 PetscValidHeaderSpecific(pc, SNES_CLASSID, 2); 4506 PetscCheckSameComm(snes, 1, pc, 2); 4507 ierr = PetscObjectReference((PetscObject) pc);CHKERRQ(ierr); 4508 ierr = SNESDestroy(&snes->pc);CHKERRQ(ierr); 4509 snes->pc = pc; 4510 ierr = PetscLogObjectParent(snes, snes->pc);CHKERRQ(ierr); 4511 PetscFunctionReturn(0); 4512 } 4513 4514 #undef __FUNCT__ 4515 #define __FUNCT__ "SNESGetPC" 4516 /*@ 4517 SNESGetPC - Returns a pointer to the nonlinear preconditioning context set with SNESSetPC(). 4518 4519 Not Collective 4520 4521 Input Parameter: 4522 . snes - iterative context obtained from SNESCreate() 4523 4524 Output Parameter: 4525 . pc - preconditioner context 4526 4527 Level: developer 4528 4529 .keywords: SNES, get, preconditioner 4530 .seealso: SNESSetPC() 4531 @*/ 4532 PetscErrorCode SNESGetPC(SNES snes, SNES *pc) 4533 { 4534 PetscErrorCode ierr; 4535 const char *optionsprefix; 4536 4537 PetscFunctionBegin; 4538 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); 4539 PetscValidPointer(pc, 2); 4540 if (!snes->pc) { 4541 ierr = SNESCreate(PetscObjectComm((PetscObject)snes),&snes->pc);CHKERRQ(ierr); 4542 ierr = PetscObjectIncrementTabLevel((PetscObject)snes->pc,(PetscObject)snes,1);CHKERRQ(ierr); 4543 ierr = PetscLogObjectParent(snes,snes->pc);CHKERRQ(ierr); 4544 ierr = SNESGetOptionsPrefix(snes,&optionsprefix);CHKERRQ(ierr); 4545 ierr = SNESSetOptionsPrefix(snes->pc,optionsprefix);CHKERRQ(ierr); 4546 ierr = SNESAppendOptionsPrefix(snes->pc,"npc_");CHKERRQ(ierr); 4547 } 4548 *pc = snes->pc; 4549 PetscFunctionReturn(0); 4550 } 4551 4552 #undef __FUNCT__ 4553 #define __FUNCT__ "SNESSetPCSide" 4554 /*@ 4555 SNESSetPCSide - Sets the preconditioning side. 4556 4557 Logically Collective on SNES 4558 4559 Input Parameter: 4560 . snes - iterative context obtained from SNESCreate() 4561 4562 Output Parameter: 4563 . side - the preconditioning side, where side is one of 4564 .vb 4565 PC_LEFT - left preconditioning (default) 4566 PC_RIGHT - right preconditioning 4567 .ve 4568 4569 Options Database Keys: 4570 . -snes_pc_side <right,left> 4571 4572 Level: intermediate 4573 4574 .keywords: SNES, set, right, left, side, preconditioner, flag 4575 4576 .seealso: SNESGetPCSide(), KSPSetPCSide() 4577 @*/ 4578 PetscErrorCode SNESSetPCSide(SNES snes,PCSide side) 4579 { 4580 PetscFunctionBegin; 4581 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4582 PetscValidLogicalCollectiveEnum(snes,side,2); 4583 snes->pcside = side; 4584 PetscFunctionReturn(0); 4585 } 4586 4587 #undef __FUNCT__ 4588 #define __FUNCT__ "SNESGetPCSide" 4589 /*@ 4590 SNESGetPCSide - Gets the preconditioning side. 4591 4592 Not Collective 4593 4594 Input Parameter: 4595 . snes - iterative context obtained from SNESCreate() 4596 4597 Output Parameter: 4598 . side - the preconditioning side, where side is one of 4599 .vb 4600 PC_LEFT - left preconditioning (default) 4601 PC_RIGHT - right preconditioning 4602 .ve 4603 4604 Level: intermediate 4605 4606 .keywords: SNES, get, right, left, side, preconditioner, flag 4607 4608 .seealso: SNESSetPCSide(), KSPGetPCSide() 4609 @*/ 4610 PetscErrorCode SNESGetPCSide(SNES snes,PCSide *side) 4611 { 4612 PetscFunctionBegin; 4613 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4614 PetscValidPointer(side,2); 4615 *side = snes->pcside; 4616 PetscFunctionReturn(0); 4617 } 4618 4619 #undef __FUNCT__ 4620 #define __FUNCT__ "SNESSetSNESLineSearch" 4621 /*@ 4622 SNESSetSNESLineSearch - Sets the linesearch on the SNES instance. 4623 4624 Collective on SNES 4625 4626 Input Parameters: 4627 + snes - iterative context obtained from SNESCreate() 4628 - linesearch - the linesearch object 4629 4630 Notes: 4631 Use SNESGetSNESLineSearch() to retrieve the preconditioner context (for example, 4632 to configure it using the API). 4633 4634 Level: developer 4635 4636 .keywords: SNES, set, linesearch 4637 .seealso: SNESGetSNESLineSearch() 4638 @*/ 4639 PetscErrorCode SNESSetSNESLineSearch(SNES snes, SNESLineSearch linesearch) 4640 { 4641 PetscErrorCode ierr; 4642 4643 PetscFunctionBegin; 4644 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); 4645 PetscValidHeaderSpecific(linesearch, SNESLINESEARCH_CLASSID, 2); 4646 PetscCheckSameComm(snes, 1, linesearch, 2); 4647 ierr = PetscObjectReference((PetscObject) linesearch);CHKERRQ(ierr); 4648 ierr = SNESLineSearchDestroy(&snes->linesearch);CHKERRQ(ierr); 4649 4650 snes->linesearch = linesearch; 4651 4652 ierr = PetscLogObjectParent(snes, snes->linesearch);CHKERRQ(ierr); 4653 PetscFunctionReturn(0); 4654 } 4655 4656 #undef __FUNCT__ 4657 #define __FUNCT__ "SNESGetSNESLineSearch" 4658 /*@C 4659 SNESGetSNESLineSearch - Returns a pointer to the line search context set with SNESSetLineSearch() 4660 or creates a default line search instance associated with the SNES and returns it. 4661 4662 Not Collective 4663 4664 Input Parameter: 4665 . snes - iterative context obtained from SNESCreate() 4666 4667 Output Parameter: 4668 . linesearch - linesearch context 4669 4670 Level: developer 4671 4672 .keywords: SNES, get, linesearch 4673 .seealso: SNESSetSNESLineSearch() 4674 @*/ 4675 PetscErrorCode SNESGetSNESLineSearch(SNES snes, SNESLineSearch *linesearch) 4676 { 4677 PetscErrorCode ierr; 4678 const char *optionsprefix; 4679 4680 PetscFunctionBegin; 4681 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); 4682 PetscValidPointer(linesearch, 2); 4683 if (!snes->linesearch) { 4684 ierr = SNESGetOptionsPrefix(snes, &optionsprefix);CHKERRQ(ierr); 4685 ierr = SNESLineSearchCreate(PetscObjectComm((PetscObject)snes), &snes->linesearch);CHKERRQ(ierr); 4686 ierr = SNESLineSearchSetSNES(snes->linesearch, snes);CHKERRQ(ierr); 4687 ierr = SNESLineSearchAppendOptionsPrefix(snes->linesearch, optionsprefix);CHKERRQ(ierr); 4688 ierr = PetscObjectIncrementTabLevel((PetscObject) snes->linesearch, (PetscObject) snes, 1);CHKERRQ(ierr); 4689 ierr = PetscLogObjectParent(snes, snes->linesearch);CHKERRQ(ierr); 4690 } 4691 *linesearch = snes->linesearch; 4692 PetscFunctionReturn(0); 4693 } 4694 4695 #if defined(PETSC_HAVE_MATLAB_ENGINE) 4696 #include <mex.h> 4697 4698 typedef struct {char *funcname; mxArray *ctx;} SNESMatlabContext; 4699 4700 #undef __FUNCT__ 4701 #define __FUNCT__ "SNESComputeFunction_Matlab" 4702 /* 4703 SNESComputeFunction_Matlab - Calls the function that has been set with SNESSetFunctionMatlab(). 4704 4705 Collective on SNES 4706 4707 Input Parameters: 4708 + snes - the SNES context 4709 - x - input vector 4710 4711 Output Parameter: 4712 . y - function vector, as set by SNESSetFunction() 4713 4714 Notes: 4715 SNESComputeFunction() is typically used within nonlinear solvers 4716 implementations, so most users would not generally call this routine 4717 themselves. 4718 4719 Level: developer 4720 4721 .keywords: SNES, nonlinear, compute, function 4722 4723 .seealso: SNESSetFunction(), SNESGetFunction() 4724 */ 4725 PetscErrorCode SNESComputeFunction_Matlab(SNES snes,Vec x,Vec y, void *ctx) 4726 { 4727 PetscErrorCode ierr; 4728 SNESMatlabContext *sctx = (SNESMatlabContext*)ctx; 4729 int nlhs = 1,nrhs = 5; 4730 mxArray *plhs[1],*prhs[5]; 4731 long long int lx = 0,ly = 0,ls = 0; 4732 4733 PetscFunctionBegin; 4734 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4735 PetscValidHeaderSpecific(x,VEC_CLASSID,2); 4736 PetscValidHeaderSpecific(y,VEC_CLASSID,3); 4737 PetscCheckSameComm(snes,1,x,2); 4738 PetscCheckSameComm(snes,1,y,3); 4739 4740 /* call Matlab function in ctx with arguments x and y */ 4741 4742 ierr = PetscMemcpy(&ls,&snes,sizeof(snes));CHKERRQ(ierr); 4743 ierr = PetscMemcpy(&lx,&x,sizeof(x));CHKERRQ(ierr); 4744 ierr = PetscMemcpy(&ly,&y,sizeof(x));CHKERRQ(ierr); 4745 prhs[0] = mxCreateDoubleScalar((double)ls); 4746 prhs[1] = mxCreateDoubleScalar((double)lx); 4747 prhs[2] = mxCreateDoubleScalar((double)ly); 4748 prhs[3] = mxCreateString(sctx->funcname); 4749 prhs[4] = sctx->ctx; 4750 ierr = mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscSNESComputeFunctionInternal");CHKERRQ(ierr); 4751 ierr = mxGetScalar(plhs[0]);CHKERRQ(ierr); 4752 mxDestroyArray(prhs[0]); 4753 mxDestroyArray(prhs[1]); 4754 mxDestroyArray(prhs[2]); 4755 mxDestroyArray(prhs[3]); 4756 mxDestroyArray(plhs[0]); 4757 PetscFunctionReturn(0); 4758 } 4759 4760 #undef __FUNCT__ 4761 #define __FUNCT__ "SNESSetFunctionMatlab" 4762 /* 4763 SNESSetFunctionMatlab - Sets the function evaluation routine and function 4764 vector for use by the SNES routines in solving systems of nonlinear 4765 equations from MATLAB. Here the function is a string containing the name of a MATLAB function 4766 4767 Logically Collective on SNES 4768 4769 Input Parameters: 4770 + snes - the SNES context 4771 . r - vector to store function value 4772 - SNESFunction - function evaluation routine 4773 4774 Notes: 4775 The Newton-like methods typically solve linear systems of the form 4776 $ f'(x) x = -f(x), 4777 where f'(x) denotes the Jacobian matrix and f(x) is the function. 4778 4779 Level: beginner 4780 4781 .keywords: SNES, nonlinear, set, function 4782 4783 .seealso: SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESSetFunction() 4784 */ 4785 PetscErrorCode SNESSetFunctionMatlab(SNES snes,Vec r,const char *SNESFunction,mxArray *ctx) 4786 { 4787 PetscErrorCode ierr; 4788 SNESMatlabContext *sctx; 4789 4790 PetscFunctionBegin; 4791 /* currently sctx is memory bleed */ 4792 ierr = PetscMalloc(sizeof(SNESMatlabContext),&sctx);CHKERRQ(ierr); 4793 ierr = PetscStrallocpy(SNESFunction,&sctx->funcname);CHKERRQ(ierr); 4794 /* 4795 This should work, but it doesn't 4796 sctx->ctx = ctx; 4797 mexMakeArrayPersistent(sctx->ctx); 4798 */ 4799 sctx->ctx = mxDuplicateArray(ctx); 4800 ierr = SNESSetFunction(snes,r,SNESComputeFunction_Matlab,sctx);CHKERRQ(ierr); 4801 PetscFunctionReturn(0); 4802 } 4803 4804 #undef __FUNCT__ 4805 #define __FUNCT__ "SNESComputeJacobian_Matlab" 4806 /* 4807 SNESComputeJacobian_Matlab - Calls the function that has been set with SNESSetJacobianMatlab(). 4808 4809 Collective on SNES 4810 4811 Input Parameters: 4812 + snes - the SNES context 4813 . x - input vector 4814 . A, B - the matrices 4815 - ctx - user context 4816 4817 Output Parameter: 4818 . flag - structure of the matrix 4819 4820 Level: developer 4821 4822 .keywords: SNES, nonlinear, compute, function 4823 4824 .seealso: SNESSetFunction(), SNESGetFunction() 4825 @*/ 4826 PetscErrorCode SNESComputeJacobian_Matlab(SNES snes,Vec x,Mat *A,Mat *B,MatStructure *flag, void *ctx) 4827 { 4828 PetscErrorCode ierr; 4829 SNESMatlabContext *sctx = (SNESMatlabContext*)ctx; 4830 int nlhs = 2,nrhs = 6; 4831 mxArray *plhs[2],*prhs[6]; 4832 long long int lx = 0,lA = 0,ls = 0, lB = 0; 4833 4834 PetscFunctionBegin; 4835 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4836 PetscValidHeaderSpecific(x,VEC_CLASSID,2); 4837 4838 /* call Matlab function in ctx with arguments x and y */ 4839 4840 ierr = PetscMemcpy(&ls,&snes,sizeof(snes));CHKERRQ(ierr); 4841 ierr = PetscMemcpy(&lx,&x,sizeof(x));CHKERRQ(ierr); 4842 ierr = PetscMemcpy(&lA,A,sizeof(x));CHKERRQ(ierr); 4843 ierr = PetscMemcpy(&lB,B,sizeof(x));CHKERRQ(ierr); 4844 prhs[0] = mxCreateDoubleScalar((double)ls); 4845 prhs[1] = mxCreateDoubleScalar((double)lx); 4846 prhs[2] = mxCreateDoubleScalar((double)lA); 4847 prhs[3] = mxCreateDoubleScalar((double)lB); 4848 prhs[4] = mxCreateString(sctx->funcname); 4849 prhs[5] = sctx->ctx; 4850 ierr = mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscSNESComputeJacobianInternal");CHKERRQ(ierr); 4851 ierr = mxGetScalar(plhs[0]);CHKERRQ(ierr); 4852 *flag = (MatStructure) mxGetScalar(plhs[1]);CHKERRQ(ierr); 4853 mxDestroyArray(prhs[0]); 4854 mxDestroyArray(prhs[1]); 4855 mxDestroyArray(prhs[2]); 4856 mxDestroyArray(prhs[3]); 4857 mxDestroyArray(prhs[4]); 4858 mxDestroyArray(plhs[0]); 4859 mxDestroyArray(plhs[1]); 4860 PetscFunctionReturn(0); 4861 } 4862 4863 #undef __FUNCT__ 4864 #define __FUNCT__ "SNESSetJacobianMatlab" 4865 /* 4866 SNESSetJacobianMatlab - Sets the Jacobian function evaluation routine and two empty Jacobian matrices 4867 vector for use by the SNES routines in solving systems of nonlinear 4868 equations from MATLAB. Here the function is a string containing the name of a MATLAB function 4869 4870 Logically Collective on SNES 4871 4872 Input Parameters: 4873 + snes - the SNES context 4874 . A,B - Jacobian matrices 4875 . SNESJacobianFunction - function evaluation routine 4876 - ctx - user context 4877 4878 Level: developer 4879 4880 .keywords: SNES, nonlinear, set, function 4881 4882 .seealso: SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESSetFunction(), SNESJacobianFunction 4883 */ 4884 PetscErrorCode SNESSetJacobianMatlab(SNES snes,Mat A,Mat B,const char *SNESJacobianFunction,mxArray *ctx) 4885 { 4886 PetscErrorCode ierr; 4887 SNESMatlabContext *sctx; 4888 4889 PetscFunctionBegin; 4890 /* currently sctx is memory bleed */ 4891 ierr = PetscMalloc(sizeof(SNESMatlabContext),&sctx);CHKERRQ(ierr); 4892 ierr = PetscStrallocpy(SNESJacobianFunction,&sctx->funcname);CHKERRQ(ierr); 4893 /* 4894 This should work, but it doesn't 4895 sctx->ctx = ctx; 4896 mexMakeArrayPersistent(sctx->ctx); 4897 */ 4898 sctx->ctx = mxDuplicateArray(ctx); 4899 ierr = SNESSetJacobian(snes,A,B,SNESComputeJacobian_Matlab,sctx);CHKERRQ(ierr); 4900 PetscFunctionReturn(0); 4901 } 4902 4903 #undef __FUNCT__ 4904 #define __FUNCT__ "SNESMonitor_Matlab" 4905 /* 4906 SNESMonitor_Matlab - Calls the function that has been set with SNESMonitorSetMatlab(). 4907 4908 Collective on SNES 4909 4910 .seealso: SNESSetFunction(), SNESGetFunction() 4911 @*/ 4912 PetscErrorCode SNESMonitor_Matlab(SNES snes,PetscInt it, PetscReal fnorm, void *ctx) 4913 { 4914 PetscErrorCode ierr; 4915 SNESMatlabContext *sctx = (SNESMatlabContext*)ctx; 4916 int nlhs = 1,nrhs = 6; 4917 mxArray *plhs[1],*prhs[6]; 4918 long long int lx = 0,ls = 0; 4919 Vec x = snes->vec_sol; 4920 4921 PetscFunctionBegin; 4922 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4923 4924 ierr = PetscMemcpy(&ls,&snes,sizeof(snes));CHKERRQ(ierr); 4925 ierr = PetscMemcpy(&lx,&x,sizeof(x));CHKERRQ(ierr); 4926 prhs[0] = mxCreateDoubleScalar((double)ls); 4927 prhs[1] = mxCreateDoubleScalar((double)it); 4928 prhs[2] = mxCreateDoubleScalar((double)fnorm); 4929 prhs[3] = mxCreateDoubleScalar((double)lx); 4930 prhs[4] = mxCreateString(sctx->funcname); 4931 prhs[5] = sctx->ctx; 4932 ierr = mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscSNESMonitorInternal");CHKERRQ(ierr); 4933 ierr = mxGetScalar(plhs[0]);CHKERRQ(ierr); 4934 mxDestroyArray(prhs[0]); 4935 mxDestroyArray(prhs[1]); 4936 mxDestroyArray(prhs[2]); 4937 mxDestroyArray(prhs[3]); 4938 mxDestroyArray(prhs[4]); 4939 mxDestroyArray(plhs[0]); 4940 PetscFunctionReturn(0); 4941 } 4942 4943 #undef __FUNCT__ 4944 #define __FUNCT__ "SNESMonitorSetMatlab" 4945 /* 4946 SNESMonitorSetMatlab - Sets the monitor function from MATLAB 4947 4948 Level: developer 4949 4950 .keywords: SNES, nonlinear, set, function 4951 4952 .seealso: SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESSetFunction() 4953 */ 4954 PetscErrorCode SNESMonitorSetMatlab(SNES snes,const char *SNESMonitorFunction,mxArray *ctx) 4955 { 4956 PetscErrorCode ierr; 4957 SNESMatlabContext *sctx; 4958 4959 PetscFunctionBegin; 4960 /* currently sctx is memory bleed */ 4961 ierr = PetscMalloc(sizeof(SNESMatlabContext),&sctx);CHKERRQ(ierr); 4962 ierr = PetscStrallocpy(SNESMonitorFunction,&sctx->funcname);CHKERRQ(ierr); 4963 /* 4964 This should work, but it doesn't 4965 sctx->ctx = ctx; 4966 mexMakeArrayPersistent(sctx->ctx); 4967 */ 4968 sctx->ctx = mxDuplicateArray(ctx); 4969 ierr = SNESMonitorSet(snes,SNESMonitor_Matlab,sctx,NULL);CHKERRQ(ierr); 4970 PetscFunctionReturn(0); 4971 } 4972 4973 #endif 4974