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