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