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; 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->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->computepfunction) { 1762 ierr = (*sdm->computepfunction)(snes,x,f,sdm->pctx);CHKERRQ(ierr); 1763 } else { 1764 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE, "Must call SNESSetPicard() to provide Picard function."); 1765 } 1766 1767 if (sdm->computepjacobian) { 1768 ierr = (*sdm->computepjacobian)(snes,x,&snes->jacobian,&snes->jacobian_pre,&snes->matstruct,sdm->pctx);CHKERRQ(ierr); 1769 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE, "Must call SNESSetPicard() to provide Picard matrix."); 1770 1771 ierr = VecView(x,PETSC_VIEWER_BINARY_WORLD);CHKERRQ(ierr); 1772 ierr = VecView(f,PETSC_VIEWER_BINARY_WORLD);CHKERRQ(ierr); 1773 ierr = VecScale(f,-1.0);CHKERRQ(ierr); 1774 ierr = MatMultAdd(snes->jacobian_pre,x,f,f);CHKERRQ(ierr); 1775 PetscFunctionReturn(0); 1776 } 1777 1778 #undef __FUNCT__ 1779 #define __FUNCT__ "SNESPicardComputeJacobian" 1780 PetscErrorCode SNESPicardComputeJacobian(SNES snes,Vec x1,Mat *J,Mat *B,MatStructure *flag,void *ctx) 1781 { 1782 PetscFunctionBegin; 1783 /* the jacobian matrix should be pre-filled in SNESPicardComputeFunction */ 1784 *flag = snes->matstruct; 1785 PetscFunctionReturn(0); 1786 } 1787 1788 #undef __FUNCT__ 1789 #define __FUNCT__ "SNESSetPicard" 1790 /*@C 1791 SNESSetPicard - Use SNES to solve the semilinear-system A(x) x = b(x) via a Picard type iteration (Picard linearization) 1792 1793 Logically Collective on SNES 1794 1795 Input Parameters: 1796 + snes - the SNES context 1797 . r - vector to store function value 1798 . func - function evaluation routine 1799 . 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) 1800 . mat - matrix to store A 1801 . mfunc - function to compute matrix value 1802 - ctx - [optional] user-defined context for private data for the 1803 function evaluation routine (may be PETSC_NULL) 1804 1805 Calling sequence of func: 1806 $ func (SNES snes,Vec x,Vec f,void *ctx); 1807 1808 + f - function vector 1809 - ctx - optional user-defined function context 1810 1811 Calling sequence of mfunc: 1812 $ mfunc (SNES snes,Vec x,Mat *jmat,Mat *mat,int *flag,void *ctx); 1813 1814 + x - input vector 1815 . 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(), 1816 normally just pass mat in this location 1817 . mat - form A(x) matrix 1818 . flag - flag indicating information about the preconditioner matrix 1819 structure (same as flag in KSPSetOperators()), one of SAME_NONZERO_PATTERN,DIFFERENT_NONZERO_PATTERN,SAME_PRECONDITIONER 1820 - ctx - [optional] user-defined Jacobian context 1821 1822 Notes: 1823 One can call SNESSetPicard() or SNESSetFunction() (and possibly SNESSetJacobian()) but cannot call both 1824 1825 $ 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} 1826 $ Note that when an exact solver is used this corresponds to the "classic" Picard A(x^{n}) x^{n+1} = b(x^{n}) iteration. 1827 1828 Run with -snes_mf_operator to solve the system with Newton's method using A(x^{n}) to construct the preconditioner. 1829 1830 We implement the defect correction form of the Picard iteration because it converges much more generally when inexact linear solvers are used then 1831 the direct Picard iteration A(x^n) x^{n+1} = b(x^n) 1832 1833 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 1834 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 1835 different please contact us at petsc-dev@mcs.anl.gov and we'll have an entirely new argument :-). 1836 1837 Level: beginner 1838 1839 .keywords: SNES, nonlinear, set, function 1840 1841 .seealso: SNESGetFunction(), SNESSetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESGetPicard(), SNESLineSearchPreCheckPicard() 1842 @*/ 1843 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) 1844 { 1845 PetscErrorCode ierr; 1846 DM dm; 1847 1848 PetscFunctionBegin; 1849 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1850 ierr = SNESGetDM(snes, &dm);CHKERRQ(ierr); 1851 ierr = DMSNESSetPicard(dm,func,mfunc,ctx);CHKERRQ(ierr); 1852 ierr = SNESSetFunction(snes,r,SNESPicardComputeFunction,ctx);CHKERRQ(ierr); 1853 ierr = SNESSetJacobian(snes,jmat,mat,SNESPicardComputeJacobian,ctx);CHKERRQ(ierr); 1854 PetscFunctionReturn(0); 1855 } 1856 1857 1858 #undef __FUNCT__ 1859 #define __FUNCT__ "SNESGetPicard" 1860 /*@C 1861 SNESGetPicard - Returns the context for the Picard iteration 1862 1863 Not Collective, but Vec is parallel if SNES is parallel. Collective if Vec is requested, but has not been created yet. 1864 1865 Input Parameter: 1866 . snes - the SNES context 1867 1868 Output Parameter: 1869 + r - the function (or PETSC_NULL) 1870 . func - the function (or PETSC_NULL) 1871 . jmat - the picard matrix (or PETSC_NULL) 1872 . mat - the picard preconditioner matrix (or PETSC_NULL) 1873 . mfunc - the function for matrix evaluation (or PETSC_NULL) 1874 - ctx - the function context (or PETSC_NULL) 1875 1876 Level: advanced 1877 1878 .keywords: SNES, nonlinear, get, function 1879 1880 .seealso: SNESSetPicard, SNESGetFunction, SNESGetJacobian, SNESGetDM 1881 @*/ 1882 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) 1883 { 1884 PetscErrorCode ierr; 1885 DM dm; 1886 1887 PetscFunctionBegin; 1888 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1889 ierr = SNESGetFunction(snes,r,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr); 1890 ierr = SNESGetJacobian(snes,jmat,mat,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr); 1891 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 1892 ierr = DMSNESGetPicard(dm,func,mfunc,ctx);CHKERRQ(ierr); 1893 PetscFunctionReturn(0); 1894 } 1895 1896 #undef __FUNCT__ 1897 #define __FUNCT__ "SNESSetComputeInitialGuess" 1898 /*@C 1899 SNESSetComputeInitialGuess - Sets a routine used to compute an initial guess for the problem 1900 1901 Logically Collective on SNES 1902 1903 Input Parameters: 1904 + snes - the SNES context 1905 . func - function evaluation routine 1906 - ctx - [optional] user-defined context for private data for the 1907 function evaluation routine (may be PETSC_NULL) 1908 1909 Calling sequence of func: 1910 $ func (SNES snes,Vec x,void *ctx); 1911 1912 . f - function vector 1913 - ctx - optional user-defined function context 1914 1915 Level: intermediate 1916 1917 .keywords: SNES, nonlinear, set, function 1918 1919 .seealso: SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian() 1920 @*/ 1921 PetscErrorCode SNESSetComputeInitialGuess(SNES snes,PetscErrorCode (*func)(SNES,Vec,void*),void *ctx) 1922 { 1923 PetscFunctionBegin; 1924 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1925 if (func) snes->ops->computeinitialguess = func; 1926 if (ctx) snes->initialguessP = ctx; 1927 PetscFunctionReturn(0); 1928 } 1929 1930 /* --------------------------------------------------------------- */ 1931 #undef __FUNCT__ 1932 #define __FUNCT__ "SNESGetRhs" 1933 /*@C 1934 SNESGetRhs - Gets the vector for solving F(x) = rhs. If rhs is not set 1935 it assumes a zero right hand side. 1936 1937 Logically Collective on SNES 1938 1939 Input Parameter: 1940 . snes - the SNES context 1941 1942 Output Parameter: 1943 . rhs - the right hand side vector or PETSC_NULL if the right hand side vector is null 1944 1945 Level: intermediate 1946 1947 .keywords: SNES, nonlinear, get, function, right hand side 1948 1949 .seealso: SNESGetSolution(), SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESSetFunction() 1950 @*/ 1951 PetscErrorCode SNESGetRhs(SNES snes,Vec *rhs) 1952 { 1953 PetscFunctionBegin; 1954 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1955 PetscValidPointer(rhs,2); 1956 *rhs = snes->vec_rhs; 1957 PetscFunctionReturn(0); 1958 } 1959 1960 #undef __FUNCT__ 1961 #define __FUNCT__ "SNESComputeFunction" 1962 /*@ 1963 SNESComputeFunction - Calls the function that has been set with 1964 SNESSetFunction(). 1965 1966 Collective on SNES 1967 1968 Input Parameters: 1969 + snes - the SNES context 1970 - x - input vector 1971 1972 Output Parameter: 1973 . y - function vector, as set by SNESSetFunction() 1974 1975 Notes: 1976 SNESComputeFunction() is typically used within nonlinear solvers 1977 implementations, so most users would not generally call this routine 1978 themselves. 1979 1980 Level: developer 1981 1982 .keywords: SNES, nonlinear, compute, function 1983 1984 .seealso: SNESSetFunction(), SNESGetFunction() 1985 @*/ 1986 PetscErrorCode SNESComputeFunction(SNES snes,Vec x,Vec y) 1987 { 1988 PetscErrorCode ierr; 1989 DM dm; 1990 DMSNES sdm; 1991 1992 PetscFunctionBegin; 1993 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1994 PetscValidHeaderSpecific(x,VEC_CLASSID,2); 1995 PetscValidHeaderSpecific(y,VEC_CLASSID,3); 1996 PetscCheckSameComm(snes,1,x,2); 1997 PetscCheckSameComm(snes,1,y,3); 1998 VecValidValues(x,2,PETSC_TRUE); 1999 2000 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 2001 ierr = DMGetDMSNES(dm,&sdm);CHKERRQ(ierr); 2002 ierr = PetscLogEventBegin(SNES_FunctionEval,snes,x,y,0);CHKERRQ(ierr); 2003 if (snes->pc && snes->pcside == PC_LEFT) { 2004 ierr = VecCopy(x,y);CHKERRQ(ierr); 2005 ierr = SNESSolve(snes->pc,snes->vec_rhs,y);CHKERRQ(ierr); 2006 ierr = VecAYPX(y,-1.0,x);CHKERRQ(ierr); 2007 } else if (sdm->computefunction) { 2008 PetscStackPush("SNES user function"); 2009 ierr = (*sdm->computefunction)(snes,x,y,sdm->functionctx);CHKERRQ(ierr); 2010 PetscStackPop; 2011 } else if (snes->vec_rhs) { 2012 ierr = MatMult(snes->jacobian, x, y);CHKERRQ(ierr); 2013 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE, "Must call SNESSetFunction() or SNESSetDM() before SNESComputeFunction(), likely called from SNESSolve()."); 2014 if (snes->vec_rhs) { 2015 ierr = VecAXPY(y,-1.0,snes->vec_rhs);CHKERRQ(ierr); 2016 } 2017 snes->nfuncs++; 2018 ierr = PetscLogEventEnd(SNES_FunctionEval,snes,x,y,0);CHKERRQ(ierr); 2019 VecValidValues(y,3,PETSC_FALSE); 2020 PetscFunctionReturn(0); 2021 } 2022 2023 #undef __FUNCT__ 2024 #define __FUNCT__ "SNESComputeGS" 2025 /*@ 2026 SNESComputeGS - Calls the Gauss-Seidel function that has been set with 2027 SNESSetGS(). 2028 2029 Collective on SNES 2030 2031 Input Parameters: 2032 + snes - the SNES context 2033 . x - input vector 2034 - b - rhs vector 2035 2036 Output Parameter: 2037 . x - new solution vector 2038 2039 Notes: 2040 SNESComputeGS() is typically used within composed nonlinear solver 2041 implementations, so most users would not generally call this routine 2042 themselves. 2043 2044 Level: developer 2045 2046 .keywords: SNES, nonlinear, compute, function 2047 2048 .seealso: SNESSetGS(), SNESComputeFunction() 2049 @*/ 2050 PetscErrorCode SNESComputeGS(SNES snes,Vec b,Vec x) 2051 { 2052 PetscErrorCode ierr; 2053 PetscInt i; 2054 DM dm; 2055 DMSNES sdm; 2056 2057 PetscFunctionBegin; 2058 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2059 PetscValidHeaderSpecific(x,VEC_CLASSID,2); 2060 if (b) PetscValidHeaderSpecific(b,VEC_CLASSID,3); 2061 PetscCheckSameComm(snes,1,x,2); 2062 if (b) PetscCheckSameComm(snes,1,b,3); 2063 if (b) VecValidValues(b,2,PETSC_TRUE); 2064 ierr = PetscLogEventBegin(SNES_GSEval,snes,x,b,0);CHKERRQ(ierr); 2065 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 2066 ierr = DMGetDMSNES(dm,&sdm);CHKERRQ(ierr); 2067 if (sdm->computegs) { 2068 for (i = 0; i < snes->gssweeps; i++) { 2069 PetscStackPush("SNES user GS"); 2070 ierr = (*sdm->computegs)(snes,x,b,sdm->gsctx);CHKERRQ(ierr); 2071 PetscStackPop; 2072 } 2073 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE, "Must call SNESSetGS() before SNESComputeGS(), likely called from SNESSolve()."); 2074 ierr = PetscLogEventEnd(SNES_GSEval,snes,x,b,0);CHKERRQ(ierr); 2075 VecValidValues(x,3,PETSC_FALSE); 2076 PetscFunctionReturn(0); 2077 } 2078 2079 2080 #undef __FUNCT__ 2081 #define __FUNCT__ "SNESComputeJacobian" 2082 /*@ 2083 SNESComputeJacobian - Computes the Jacobian matrix that has been 2084 set with SNESSetJacobian(). 2085 2086 Collective on SNES and Mat 2087 2088 Input Parameters: 2089 + snes - the SNES context 2090 - x - input vector 2091 2092 Output Parameters: 2093 + A - Jacobian matrix 2094 . B - optional preconditioning matrix 2095 - flag - flag indicating matrix structure (one of, SAME_NONZERO_PATTERN,DIFFERENT_NONZERO_PATTERN,SAME_PRECONDITIONER) 2096 2097 Options Database Keys: 2098 + -snes_lag_preconditioner <lag> 2099 . -snes_lag_jacobian <lag> 2100 . -snes_compare_explicit - Compare the computed Jacobian to the finite difference Jacobian and output the differences 2101 . -snes_compare_explicit_draw - Compare the computed Jacobian to the finite difference Jacobian and draw the result 2102 . -snes_compare_explicit_contour - Compare the computed Jacobian to the finite difference Jacobian and draw a contour plot with the result 2103 . -snes_compare_operator - Make the comparison options above use the operator instead of the preconditioning matrix 2104 . -snes_compare_coloring - Compute the finite differece Jacobian using coloring and display norms of difference 2105 . -snes_compare_coloring_display - Compute the finite differece Jacobian using coloring and display verbose differences 2106 . -snes_compare_coloring_threshold - Display only those matrix entries that differ by more than a given threshold 2107 . -snes_compare_coloring_threshold_atol - Absolute tolerance for difference in matrix entries to be displayed by -snes_compare_coloring_threshold 2108 . -snes_compare_coloring_threshold_rtol - Relative tolerance for difference in matrix entries to be displayed by -snes_compare_coloring_threshold 2109 . -snes_compare_coloring_draw - Compute the finite differece Jacobian using coloring and draw differences 2110 - -snes_compare_coloring_draw_contour - Compute the finite differece Jacobian using coloring and show contours of matrices and differences 2111 2112 2113 Notes: 2114 Most users should not need to explicitly call this routine, as it 2115 is used internally within the nonlinear solvers. 2116 2117 See KSPSetOperators() for important information about setting the 2118 flag parameter. 2119 2120 Level: developer 2121 2122 .keywords: SNES, compute, Jacobian, matrix 2123 2124 .seealso: SNESSetJacobian(), KSPSetOperators(), MatStructure, SNESSetLagPreconditioner(), SNESSetLagJacobian() 2125 @*/ 2126 PetscErrorCode SNESComputeJacobian(SNES snes,Vec X,Mat *A,Mat *B,MatStructure *flg) 2127 { 2128 PetscErrorCode ierr; 2129 PetscBool flag; 2130 DM dm; 2131 DMSNES sdm; 2132 2133 PetscFunctionBegin; 2134 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2135 PetscValidHeaderSpecific(X,VEC_CLASSID,2); 2136 PetscValidPointer(flg,5); 2137 PetscCheckSameComm(snes,1,X,2); 2138 VecValidValues(X,2,PETSC_TRUE); 2139 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 2140 ierr = DMGetDMSNES(dm,&sdm);CHKERRQ(ierr); 2141 2142 if (!sdm->computejacobian) SETERRQ(((PetscObject)snes)->comm,PETSC_ERR_USER,"Must call SNESSetJacobian(), DMSNESSetJacobian(), DMDASNESSetJacobianLocal(), etc"); 2143 2144 /* make sure that MatAssemblyBegin/End() is called on A matrix if it is matrix free */ 2145 2146 if (snes->lagjacobian == -2) { 2147 snes->lagjacobian = -1; 2148 ierr = PetscInfo(snes,"Recomputing Jacobian/preconditioner because lag is -2 (means compute Jacobian, but then never again) \n");CHKERRQ(ierr); 2149 } else if (snes->lagjacobian == -1) { 2150 *flg = SAME_PRECONDITIONER; 2151 ierr = PetscInfo(snes,"Reusing Jacobian/preconditioner because lag is -1\n");CHKERRQ(ierr); 2152 ierr = PetscObjectTypeCompare((PetscObject)*A,MATMFFD,&flag);CHKERRQ(ierr); 2153 if (flag) { 2154 ierr = MatAssemblyBegin(*A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2155 ierr = MatAssemblyEnd(*A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2156 } 2157 PetscFunctionReturn(0); 2158 } else if (snes->lagjacobian > 1 && snes->iter % snes->lagjacobian) { 2159 *flg = SAME_PRECONDITIONER; 2160 ierr = PetscInfo2(snes,"Reusing Jacobian/preconditioner because lag is %D and SNES iteration is %D\n",snes->lagjacobian,snes->iter);CHKERRQ(ierr); 2161 ierr = PetscObjectTypeCompare((PetscObject)*A,MATMFFD,&flag);CHKERRQ(ierr); 2162 if (flag) { 2163 ierr = MatAssemblyBegin(*A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2164 ierr = MatAssemblyEnd(*A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2165 } 2166 PetscFunctionReturn(0); 2167 } 2168 2169 *flg = DIFFERENT_NONZERO_PATTERN; 2170 ierr = PetscLogEventBegin(SNES_JacobianEval,snes,X,*A,*B);CHKERRQ(ierr); 2171 2172 PetscStackPush("SNES user Jacobian function"); 2173 ierr = (*sdm->computejacobian)(snes,X,A,B,flg,sdm->jacobianctx);CHKERRQ(ierr); 2174 PetscStackPop; 2175 2176 ierr = PetscLogEventEnd(SNES_JacobianEval,snes,X,*A,*B);CHKERRQ(ierr); 2177 2178 if (snes->lagpreconditioner == -2) { 2179 ierr = PetscInfo(snes,"Rebuilding preconditioner exactly once since lag is -2\n");CHKERRQ(ierr); 2180 snes->lagpreconditioner = -1; 2181 } else if (snes->lagpreconditioner == -1) { 2182 *flg = SAME_PRECONDITIONER; 2183 ierr = PetscInfo(snes,"Reusing preconditioner because lag is -1\n");CHKERRQ(ierr); 2184 } else if (snes->lagpreconditioner > 1 && snes->iter % snes->lagpreconditioner) { 2185 *flg = SAME_PRECONDITIONER; 2186 ierr = PetscInfo2(snes,"Reusing preconditioner because lag is %D and SNES iteration is %D\n",snes->lagpreconditioner,snes->iter);CHKERRQ(ierr); 2187 } 2188 2189 /* make sure user returned a correct Jacobian and preconditioner */ 2190 /* PetscValidHeaderSpecific(*A,MAT_CLASSID,3); 2191 PetscValidHeaderSpecific(*B,MAT_CLASSID,4); */ 2192 { 2193 PetscBool flag = PETSC_FALSE,flag_draw = PETSC_FALSE,flag_contour = PETSC_FALSE,flag_operator = PETSC_FALSE; 2194 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_explicit",&flag,PETSC_NULL);CHKERRQ(ierr); 2195 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_explicit_draw",&flag_draw,PETSC_NULL);CHKERRQ(ierr); 2196 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_explicit_draw_contour",&flag_contour,PETSC_NULL);CHKERRQ(ierr); 2197 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_operator",&flag_operator,PETSC_NULL);CHKERRQ(ierr); 2198 if (flag || flag_draw || flag_contour) { 2199 Mat Bexp_mine = PETSC_NULL,Bexp,FDexp; 2200 MatStructure mstruct; 2201 PetscViewer vdraw,vstdout; 2202 PetscBool flg; 2203 if (flag_operator) { 2204 ierr = MatComputeExplicitOperator(*A,&Bexp_mine);CHKERRQ(ierr); 2205 Bexp = Bexp_mine; 2206 } else { 2207 /* See if the preconditioning matrix can be viewed and added directly */ 2208 ierr = PetscObjectTypeCompareAny((PetscObject)*B,&flg,MATSEQAIJ,MATMPIAIJ,MATSEQDENSE,MATMPIDENSE,MATSEQBAIJ,MATMPIBAIJ,MATSEQSBAIJ,MATMPIBAIJ,"");CHKERRQ(ierr); 2209 if (flg) Bexp = *B; 2210 else { 2211 /* If the "preconditioning" matrix is itself MATSHELL or some other type without direct support */ 2212 ierr = MatComputeExplicitOperator(*B,&Bexp_mine);CHKERRQ(ierr); 2213 Bexp = Bexp_mine; 2214 } 2215 } 2216 ierr = MatConvert(Bexp,MATSAME,MAT_INITIAL_MATRIX,&FDexp);CHKERRQ(ierr); 2217 ierr = SNESDefaultComputeJacobian(snes,X,&FDexp,&FDexp,&mstruct,NULL);CHKERRQ(ierr); 2218 ierr = PetscViewerASCIIGetStdout(((PetscObject)snes)->comm,&vstdout);CHKERRQ(ierr); 2219 if (flag_draw || flag_contour) { 2220 ierr = PetscViewerDrawOpen(((PetscObject)snes)->comm,0,"Explicit Jacobians",PETSC_DECIDE,PETSC_DECIDE,300,300,&vdraw);CHKERRQ(ierr); 2221 if (flag_contour) {ierr = PetscViewerPushFormat(vdraw,PETSC_VIEWER_DRAW_CONTOUR);CHKERRQ(ierr);} 2222 } else vdraw = PETSC_NULL; 2223 ierr = PetscViewerASCIIPrintf(vstdout,"Explicit %s\n",flag_operator?"Jacobian":"preconditioning Jacobian");CHKERRQ(ierr); 2224 if (flag) {ierr = MatView(Bexp,vstdout);CHKERRQ(ierr);} 2225 if (vdraw) {ierr = MatView(Bexp,vdraw);CHKERRQ(ierr);} 2226 ierr = PetscViewerASCIIPrintf(vstdout,"Finite difference Jacobian\n");CHKERRQ(ierr); 2227 if (flag) {ierr = MatView(FDexp,vstdout);CHKERRQ(ierr);} 2228 if (vdraw) {ierr = MatView(FDexp,vdraw);CHKERRQ(ierr);} 2229 ierr = MatAYPX(FDexp,-1.0,Bexp,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 2230 ierr = PetscViewerASCIIPrintf(vstdout,"User-provided matrix minus finite difference Jacobian\n");CHKERRQ(ierr); 2231 if (flag) {ierr = MatView(FDexp,vstdout);CHKERRQ(ierr);} 2232 if (vdraw) { /* Always use contour for the difference */ 2233 ierr = PetscViewerPushFormat(vdraw,PETSC_VIEWER_DRAW_CONTOUR);CHKERRQ(ierr); 2234 ierr = MatView(FDexp,vdraw);CHKERRQ(ierr); 2235 ierr = PetscViewerPopFormat(vdraw);CHKERRQ(ierr); 2236 } 2237 if (flag_contour) {ierr = PetscViewerPopFormat(vdraw);CHKERRQ(ierr);} 2238 ierr = PetscViewerDestroy(&vdraw);CHKERRQ(ierr); 2239 ierr = MatDestroy(&Bexp_mine);CHKERRQ(ierr); 2240 ierr = MatDestroy(&FDexp);CHKERRQ(ierr); 2241 } 2242 } 2243 { 2244 PetscBool flag = PETSC_FALSE,flag_display = PETSC_FALSE,flag_draw = PETSC_FALSE,flag_contour = PETSC_FALSE,flag_threshold = PETSC_FALSE; 2245 PetscReal threshold_atol = PETSC_SQRT_MACHINE_EPSILON,threshold_rtol = 10*PETSC_SQRT_MACHINE_EPSILON; 2246 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_coloring",&flag,PETSC_NULL);CHKERRQ(ierr); 2247 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_coloring_display",&flag_display,PETSC_NULL);CHKERRQ(ierr); 2248 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_coloring_draw",&flag_draw,PETSC_NULL);CHKERRQ(ierr); 2249 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_coloring_draw_contour",&flag_contour,PETSC_NULL);CHKERRQ(ierr); 2250 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_coloring_threshold",&flag_threshold,PETSC_NULL);CHKERRQ(ierr); 2251 ierr = PetscOptionsGetReal(((PetscObject)snes)->prefix,"-snes_compare_coloring_threshold_rtol",&threshold_rtol,PETSC_NULL);CHKERRQ(ierr); 2252 ierr = PetscOptionsGetReal(((PetscObject)snes)->prefix,"-snes_compare_coloring_threshold_atol",&threshold_atol,PETSC_NULL);CHKERRQ(ierr); 2253 if (flag || flag_display || flag_draw || flag_contour || flag_threshold) { 2254 Mat Bfd; 2255 MatStructure mstruct; 2256 PetscViewer vdraw,vstdout; 2257 ISColoring iscoloring; 2258 MatFDColoring matfdcoloring; 2259 PetscErrorCode (*func)(SNES,Vec,Vec,void*); 2260 void *funcctx; 2261 PetscReal norm1,norm2,normmax; 2262 2263 ierr = MatDuplicate(*B,MAT_DO_NOT_COPY_VALUES,&Bfd);CHKERRQ(ierr); 2264 ierr = MatGetColoring(Bfd,MATCOLORINGSL,&iscoloring);CHKERRQ(ierr); 2265 ierr = MatFDColoringCreate(Bfd,iscoloring,&matfdcoloring);CHKERRQ(ierr); 2266 ierr = ISColoringDestroy(&iscoloring);CHKERRQ(ierr); 2267 2268 /* This method of getting the function is currently unreliable since it doesn't work for DM local functions. */ 2269 ierr = SNESGetFunction(snes,PETSC_NULL,&func,&funcctx);CHKERRQ(ierr); 2270 ierr = MatFDColoringSetFunction(matfdcoloring,(PetscErrorCode(*)(void))func,funcctx);CHKERRQ(ierr); 2271 ierr = PetscObjectSetOptionsPrefix((PetscObject)matfdcoloring,((PetscObject)snes)->prefix);CHKERRQ(ierr); 2272 ierr = PetscObjectAppendOptionsPrefix((PetscObject)matfdcoloring,"coloring_");CHKERRQ(ierr); 2273 ierr = MatFDColoringSetFromOptions(matfdcoloring);CHKERRQ(ierr); 2274 ierr = MatFDColoringApply(Bfd,matfdcoloring,X,&mstruct,snes);CHKERRQ(ierr); 2275 ierr = MatFDColoringDestroy(&matfdcoloring);CHKERRQ(ierr); 2276 2277 ierr = PetscViewerASCIIGetStdout(((PetscObject)snes)->comm,&vstdout);CHKERRQ(ierr); 2278 if (flag_draw || flag_contour) { 2279 ierr = PetscViewerDrawOpen(((PetscObject)snes)->comm,0,"Colored Jacobians",PETSC_DECIDE,PETSC_DECIDE,300,300,&vdraw);CHKERRQ(ierr); 2280 if (flag_contour) {ierr = PetscViewerPushFormat(vdraw,PETSC_VIEWER_DRAW_CONTOUR);CHKERRQ(ierr);} 2281 } else vdraw = PETSC_NULL; 2282 ierr = PetscViewerASCIIPrintf(vstdout,"Explicit preconditioning Jacobian\n");CHKERRQ(ierr); 2283 if (flag_display) {ierr = MatView(*B,vstdout);CHKERRQ(ierr);} 2284 if (vdraw) {ierr = MatView(*B,vdraw);CHKERRQ(ierr);} 2285 ierr = PetscViewerASCIIPrintf(vstdout,"Colored Finite difference Jacobian\n");CHKERRQ(ierr); 2286 if (flag_display) {ierr = MatView(Bfd,vstdout);CHKERRQ(ierr);} 2287 if (vdraw) {ierr = MatView(Bfd,vdraw);CHKERRQ(ierr);} 2288 ierr = MatAYPX(Bfd,-1.0,*B,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 2289 ierr = MatNorm(Bfd,NORM_1,&norm1);CHKERRQ(ierr); 2290 ierr = MatNorm(Bfd,NORM_FROBENIUS,&norm2);CHKERRQ(ierr); 2291 ierr = MatNorm(Bfd,NORM_MAX,&normmax);CHKERRQ(ierr); 2292 ierr = PetscViewerASCIIPrintf(vstdout,"User-provided matrix minus finite difference Jacobian, norm1=%G normFrob=%G normmax=%G\n",norm1,norm2,normmax);CHKERRQ(ierr); 2293 if (flag_display) {ierr = MatView(Bfd,vstdout);CHKERRQ(ierr);} 2294 if (vdraw) { /* Always use contour for the difference */ 2295 ierr = PetscViewerPushFormat(vdraw,PETSC_VIEWER_DRAW_CONTOUR);CHKERRQ(ierr); 2296 ierr = MatView(Bfd,vdraw);CHKERRQ(ierr); 2297 ierr = PetscViewerPopFormat(vdraw);CHKERRQ(ierr); 2298 } 2299 if (flag_contour) {ierr = PetscViewerPopFormat(vdraw);CHKERRQ(ierr);} 2300 2301 if (flag_threshold) { 2302 PetscInt bs,rstart,rend,i; 2303 ierr = MatGetBlockSize(*B,&bs);CHKERRQ(ierr); 2304 ierr = MatGetOwnershipRange(*B,&rstart,&rend);CHKERRQ(ierr); 2305 for (i=rstart; i<rend; i++) { 2306 const PetscScalar *ba,*ca; 2307 const PetscInt *bj,*cj; 2308 PetscInt bn,cn,j,maxentrycol = -1,maxdiffcol = -1,maxrdiffcol = -1; 2309 PetscReal maxentry = 0,maxdiff = 0,maxrdiff = 0; 2310 ierr = MatGetRow(*B,i,&bn,&bj,&ba);CHKERRQ(ierr); 2311 ierr = MatGetRow(Bfd,i,&cn,&cj,&ca);CHKERRQ(ierr); 2312 if (bn != cn) SETERRQ(((PetscObject)*A)->comm,PETSC_ERR_PLIB,"Unexpected different nonzero pattern in -snes_compare_coloring_threshold"); 2313 for (j=0; j<bn; j++) { 2314 PetscReal rdiff = PetscAbsScalar(ca[j]) / (threshold_atol + threshold_rtol*PetscAbsScalar(ba[j])); 2315 if (PetscAbsScalar(ba[j]) > PetscAbs(maxentry)) { 2316 maxentrycol = bj[j]; 2317 maxentry = PetscRealPart(ba[j]); 2318 } 2319 if (PetscAbsScalar(ca[j]) > PetscAbs(maxdiff)) { 2320 maxdiffcol = bj[j]; 2321 maxdiff = PetscRealPart(ca[j]); 2322 } 2323 if (rdiff > maxrdiff) { 2324 maxrdiffcol = bj[j]; 2325 maxrdiff = rdiff; 2326 } 2327 } 2328 if (maxrdiff > 1) { 2329 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); 2330 for (j=0; j<bn; j++) { 2331 PetscReal rdiff; 2332 rdiff = PetscAbsScalar(ca[j]) / (threshold_atol + threshold_rtol*PetscAbsScalar(ba[j])); 2333 if (rdiff > 1) { 2334 ierr = PetscViewerASCIIPrintf(vstdout," (%D,%G:%G)",bj[j],PetscRealPart(ba[j]),PetscRealPart(ca[j]));CHKERRQ(ierr); 2335 } 2336 } 2337 ierr = PetscViewerASCIIPrintf(vstdout,"\n",i,maxentry,maxdiff,maxrdiff);CHKERRQ(ierr); 2338 } 2339 ierr = MatRestoreRow(*B,i,&bn,&bj,&ba);CHKERRQ(ierr); 2340 ierr = MatRestoreRow(Bfd,i,&cn,&cj,&ca);CHKERRQ(ierr); 2341 } 2342 } 2343 ierr = PetscViewerDestroy(&vdraw);CHKERRQ(ierr); 2344 ierr = MatDestroy(&Bfd);CHKERRQ(ierr); 2345 } 2346 } 2347 PetscFunctionReturn(0); 2348 } 2349 2350 #undef __FUNCT__ 2351 #define __FUNCT__ "SNESSetJacobian" 2352 /*@C 2353 SNESSetJacobian - Sets the function to compute Jacobian as well as the 2354 location to store the matrix. 2355 2356 Logically Collective on SNES and Mat 2357 2358 Input Parameters: 2359 + snes - the SNES context 2360 . A - Jacobian matrix 2361 . B - preconditioner matrix (usually same as the Jacobian) 2362 . func - Jacobian evaluation routine (if PETSC_NULL then SNES retains any previously set value) 2363 - ctx - [optional] user-defined context for private data for the 2364 Jacobian evaluation routine (may be PETSC_NULL) (if PETSC_NULL then SNES retains any previously set value) 2365 2366 Calling sequence of func: 2367 $ func (SNES snes,Vec x,Mat *A,Mat *B,int *flag,void *ctx); 2368 2369 + x - input vector 2370 . A - Jacobian matrix 2371 . B - preconditioner matrix, usually the same as A 2372 . flag - flag indicating information about the preconditioner matrix 2373 structure (same as flag in KSPSetOperators()), one of SAME_NONZERO_PATTERN,DIFFERENT_NONZERO_PATTERN,SAME_PRECONDITIONER 2374 - ctx - [optional] user-defined Jacobian context 2375 2376 Notes: 2377 See KSPSetOperators() for important information about setting the flag 2378 output parameter in the routine func(). Be sure to read this information! 2379 2380 The routine func() takes Mat * as the matrix arguments rather than Mat. 2381 This allows the Jacobian evaluation routine to replace A and/or B with a 2382 completely new new matrix structure (not just different matrix elements) 2383 when appropriate, for instance, if the nonzero structure is changing 2384 throughout the global iterations. 2385 2386 If the A matrix and B matrix are different you must call MatAssemblyBegin/End() on 2387 each matrix. 2388 2389 If using SNESDefaultComputeJacobianColor() to assemble a Jacobian, the ctx argument 2390 must be a MatFDColoring. 2391 2392 Other defect-correction schemes can be used by computing a different matrix in place of the Jacobian. One common 2393 example is to use the "Picard linearization" which only differentiates through the highest order parts of each term. 2394 2395 Level: beginner 2396 2397 .keywords: SNES, nonlinear, set, Jacobian, matrix 2398 2399 .seealso: KSPSetOperators(), SNESSetFunction(), MatMFFDComputeJacobian(), SNESDefaultComputeJacobianColor(), MatStructure 2400 @*/ 2401 PetscErrorCode SNESSetJacobian(SNES snes,Mat A,Mat B,PetscErrorCode (*func)(SNES,Vec,Mat*,Mat*,MatStructure*,void*),void *ctx) 2402 { 2403 PetscErrorCode ierr; 2404 DM dm; 2405 2406 PetscFunctionBegin; 2407 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2408 if (A) PetscValidHeaderSpecific(A,MAT_CLASSID,2); 2409 if (B) PetscValidHeaderSpecific(B,MAT_CLASSID,3); 2410 if (A) PetscCheckSameComm(snes,1,A,2); 2411 if (B) PetscCheckSameComm(snes,1,B,3); 2412 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 2413 ierr = DMSNESSetJacobian(dm,func,ctx);CHKERRQ(ierr); 2414 if (A) { 2415 ierr = PetscObjectReference((PetscObject)A);CHKERRQ(ierr); 2416 ierr = MatDestroy(&snes->jacobian);CHKERRQ(ierr); 2417 snes->jacobian = A; 2418 } 2419 if (B) { 2420 ierr = PetscObjectReference((PetscObject)B);CHKERRQ(ierr); 2421 ierr = MatDestroy(&snes->jacobian_pre);CHKERRQ(ierr); 2422 snes->jacobian_pre = B; 2423 } 2424 PetscFunctionReturn(0); 2425 } 2426 2427 #undef __FUNCT__ 2428 #define __FUNCT__ "SNESGetJacobian" 2429 /*@C 2430 SNESGetJacobian - Returns the Jacobian matrix and optionally the user 2431 provided context for evaluating the Jacobian. 2432 2433 Not Collective, but Mat object will be parallel if SNES object is 2434 2435 Input Parameter: 2436 . snes - the nonlinear solver context 2437 2438 Output Parameters: 2439 + A - location to stash Jacobian matrix (or PETSC_NULL) 2440 . B - location to stash preconditioner matrix (or PETSC_NULL) 2441 . func - location to put Jacobian function (or PETSC_NULL) 2442 - ctx - location to stash Jacobian ctx (or PETSC_NULL) 2443 2444 Level: advanced 2445 2446 .seealso: SNESSetJacobian(), SNESComputeJacobian() 2447 @*/ 2448 PetscErrorCode SNESGetJacobian(SNES snes,Mat *A,Mat *B,PetscErrorCode (**func)(SNES,Vec,Mat*,Mat*,MatStructure*,void*),void **ctx) 2449 { 2450 PetscErrorCode ierr; 2451 DM dm; 2452 DMSNES sdm; 2453 2454 PetscFunctionBegin; 2455 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2456 if (A) *A = snes->jacobian; 2457 if (B) *B = snes->jacobian_pre; 2458 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 2459 ierr = DMGetDMSNES(dm,&sdm);CHKERRQ(ierr); 2460 if (func) *func = sdm->computejacobian; 2461 if (ctx) *ctx = sdm->jacobianctx; 2462 PetscFunctionReturn(0); 2463 } 2464 2465 /* ----- Routines to initialize and destroy a nonlinear solver ---- */ 2466 2467 #undef __FUNCT__ 2468 #define __FUNCT__ "SNESSetUp" 2469 /*@ 2470 SNESSetUp - Sets up the internal data structures for the later use 2471 of a nonlinear solver. 2472 2473 Collective on SNES 2474 2475 Input Parameters: 2476 . snes - the SNES context 2477 2478 Notes: 2479 For basic use of the SNES solvers the user need not explicitly call 2480 SNESSetUp(), since these actions will automatically occur during 2481 the call to SNESSolve(). However, if one wishes to control this 2482 phase separately, SNESSetUp() should be called after SNESCreate() 2483 and optional routines of the form SNESSetXXX(), but before SNESSolve(). 2484 2485 Level: advanced 2486 2487 .keywords: SNES, nonlinear, setup 2488 2489 .seealso: SNESCreate(), SNESSolve(), SNESDestroy() 2490 @*/ 2491 PetscErrorCode SNESSetUp(SNES snes) 2492 { 2493 PetscErrorCode ierr; 2494 DM dm; 2495 DMSNES sdm; 2496 SNESLineSearch linesearch, pclinesearch; 2497 void *lsprectx,*lspostctx; 2498 SNESLineSearchPreCheckFunc lsprefunc; 2499 SNESLineSearchPostCheckFunc lspostfunc; 2500 PetscErrorCode (*func)(SNES,Vec,Vec,void*); 2501 Vec f,fpc; 2502 void *funcctx; 2503 PetscErrorCode (*jac)(SNES,Vec,Mat*,Mat*,MatStructure*,void*); 2504 void *jacctx,*appctx; 2505 Mat A,B; 2506 2507 PetscFunctionBegin; 2508 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2509 if (snes->setupcalled) PetscFunctionReturn(0); 2510 2511 if (!((PetscObject)snes)->type_name) { 2512 ierr = SNESSetType(snes,SNESLS);CHKERRQ(ierr); 2513 } 2514 2515 ierr = SNESGetFunction(snes,&snes->vec_func,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr); 2516 if (snes->vec_func == snes->vec_sol) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_IDN,"Solution vector cannot be function vector"); 2517 if (snes->vec_rhs == snes->vec_sol) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_IDN,"Solution vector cannot be right hand side vector"); 2518 2519 if (!snes->vec_sol_update /* && snes->vec_sol */) { 2520 ierr = VecDuplicate(snes->vec_sol,&snes->vec_sol_update);CHKERRQ(ierr); 2521 ierr = PetscLogObjectParent(snes,snes->vec_sol_update);CHKERRQ(ierr); 2522 } 2523 2524 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 2525 ierr = DMShellSetGlobalVector(snes->dm,snes->vec_sol);CHKERRQ(ierr); 2526 ierr = DMGetDMSNES(dm,&sdm);CHKERRQ(ierr); 2527 if (!sdm->computefunction) SETERRQ(((PetscObject)dm)->comm,PETSC_ERR_ARG_WRONGSTATE,"Function never provided to SNES object"); 2528 if (!sdm->computejacobian) { 2529 ierr = DMSNESSetJacobian(dm,SNESDefaultComputeJacobianColor,PETSC_NULL);CHKERRQ(ierr); 2530 } 2531 if (!snes->vec_func) { 2532 ierr = DMCreateGlobalVector(dm,&snes->vec_func);CHKERRQ(ierr); 2533 } 2534 2535 if (!snes->ksp) {ierr = SNESGetKSP(snes, &snes->ksp);CHKERRQ(ierr);} 2536 2537 if (!snes->linesearch) {ierr = SNESGetSNESLineSearch(snes, &snes->linesearch);} 2538 2539 if (snes->pc && (snes->pcside == PC_LEFT)) snes->mf = PETSC_TRUE; 2540 2541 if (snes->mf) { 2542 ierr = SNESSetUpMatrixFree_Private(snes, snes->mf_operator, snes->mf_version);CHKERRQ(ierr); 2543 } 2544 2545 2546 if (snes->ops->usercompute && !snes->user) { 2547 ierr = (*snes->ops->usercompute)(snes,(void**)&snes->user);CHKERRQ(ierr); 2548 } 2549 2550 if (snes->pc) { 2551 /* copy the DM over */ 2552 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 2553 ierr = SNESSetDM(snes->pc,dm);CHKERRQ(ierr); 2554 2555 ierr = SNESGetFunction(snes,&f,&func,&funcctx);CHKERRQ(ierr); 2556 ierr = VecDuplicate(f,&fpc);CHKERRQ(ierr); 2557 ierr = SNESSetFunction(snes->pc,fpc,func,funcctx);CHKERRQ(ierr); 2558 ierr = SNESGetJacobian(snes,&A,&B,&jac,&jacctx);CHKERRQ(ierr); 2559 ierr = SNESSetJacobian(snes->pc,A,B,jac,jacctx);CHKERRQ(ierr); 2560 ierr = SNESGetApplicationContext(snes,&appctx);CHKERRQ(ierr); 2561 ierr = SNESSetApplicationContext(snes->pc,appctx);CHKERRQ(ierr); 2562 ierr = VecDestroy(&fpc);CHKERRQ(ierr); 2563 2564 /* copy the function pointers over */ 2565 ierr = PetscObjectCopyFortranFunctionPointers((PetscObject)snes,(PetscObject)snes->pc);CHKERRQ(ierr); 2566 2567 /* default to 1 iteration */ 2568 ierr = SNESSetTolerances(snes->pc, 0.0, 0.0, 0.0, 1, snes->pc->max_funcs);CHKERRQ(ierr); 2569 ierr = SNESSetNormType(snes->pc, SNES_NORM_FINAL_ONLY);CHKERRQ(ierr); 2570 ierr = SNESSetFromOptions(snes->pc);CHKERRQ(ierr); 2571 2572 /* copy the line search context over */ 2573 ierr = SNESGetSNESLineSearch(snes,&linesearch);CHKERRQ(ierr); 2574 ierr = SNESGetSNESLineSearch(snes->pc,&pclinesearch);CHKERRQ(ierr); 2575 ierr = SNESLineSearchGetPreCheck(linesearch,&lsprefunc,&lsprectx);CHKERRQ(ierr); 2576 ierr = SNESLineSearchGetPostCheck(linesearch,&lspostfunc,&lspostctx);CHKERRQ(ierr); 2577 ierr = SNESLineSearchSetPreCheck(pclinesearch,lsprefunc,lsprectx);CHKERRQ(ierr); 2578 ierr = SNESLineSearchSetPostCheck(pclinesearch,lspostfunc,lspostctx);CHKERRQ(ierr); 2579 ierr = PetscObjectCopyFortranFunctionPointers((PetscObject)linesearch, (PetscObject)pclinesearch);CHKERRQ(ierr); 2580 } 2581 2582 if (snes->ops->setup) { 2583 ierr = (*snes->ops->setup)(snes);CHKERRQ(ierr); 2584 } 2585 2586 snes->setupcalled = PETSC_TRUE; 2587 PetscFunctionReturn(0); 2588 } 2589 2590 #undef __FUNCT__ 2591 #define __FUNCT__ "SNESReset" 2592 /*@ 2593 SNESReset - Resets a SNES context to the snessetupcalled = 0 state and removes any allocated Vecs and Mats 2594 2595 Collective on SNES 2596 2597 Input Parameter: 2598 . snes - iterative context obtained from SNESCreate() 2599 2600 Level: intermediate 2601 2602 Notes: Also calls the application context destroy routine set with SNESSetComputeApplicationContext() 2603 2604 .keywords: SNES, destroy 2605 2606 .seealso: SNESCreate(), SNESSetUp(), SNESSolve() 2607 @*/ 2608 PetscErrorCode SNESReset(SNES snes) 2609 { 2610 PetscErrorCode ierr; 2611 2612 PetscFunctionBegin; 2613 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2614 if (snes->ops->userdestroy && snes->user) { 2615 ierr = (*snes->ops->userdestroy)((void**)&snes->user);CHKERRQ(ierr); 2616 snes->user = PETSC_NULL; 2617 } 2618 if (snes->pc) { 2619 ierr = SNESReset(snes->pc);CHKERRQ(ierr); 2620 } 2621 2622 if (snes->ops->reset) { 2623 ierr = (*snes->ops->reset)(snes);CHKERRQ(ierr); 2624 } 2625 if (snes->ksp) { 2626 ierr = KSPReset(snes->ksp);CHKERRQ(ierr); 2627 } 2628 2629 if (snes->linesearch) { 2630 ierr = SNESLineSearchReset(snes->linesearch);CHKERRQ(ierr); 2631 } 2632 2633 ierr = VecDestroy(&snes->vec_rhs);CHKERRQ(ierr); 2634 ierr = VecDestroy(&snes->vec_sol);CHKERRQ(ierr); 2635 ierr = VecDestroy(&snes->vec_sol_update);CHKERRQ(ierr); 2636 ierr = VecDestroy(&snes->vec_func);CHKERRQ(ierr); 2637 ierr = MatDestroy(&snes->jacobian);CHKERRQ(ierr); 2638 ierr = MatDestroy(&snes->jacobian_pre);CHKERRQ(ierr); 2639 ierr = VecDestroyVecs(snes->nwork,&snes->work);CHKERRQ(ierr); 2640 ierr = VecDestroyVecs(snes->nvwork,&snes->vwork);CHKERRQ(ierr); 2641 snes->nwork = snes->nvwork = 0; 2642 snes->setupcalled = PETSC_FALSE; 2643 PetscFunctionReturn(0); 2644 } 2645 2646 #undef __FUNCT__ 2647 #define __FUNCT__ "SNESDestroy" 2648 /*@ 2649 SNESDestroy - Destroys the nonlinear solver context that was created 2650 with SNESCreate(). 2651 2652 Collective on SNES 2653 2654 Input Parameter: 2655 . snes - the SNES context 2656 2657 Level: beginner 2658 2659 .keywords: SNES, nonlinear, destroy 2660 2661 .seealso: SNESCreate(), SNESSolve() 2662 @*/ 2663 PetscErrorCode SNESDestroy(SNES *snes) 2664 { 2665 PetscErrorCode ierr; 2666 2667 PetscFunctionBegin; 2668 if (!*snes) PetscFunctionReturn(0); 2669 PetscValidHeaderSpecific((*snes),SNES_CLASSID,1); 2670 if (--((PetscObject)(*snes))->refct > 0) {*snes = 0; PetscFunctionReturn(0);} 2671 2672 ierr = SNESReset((*snes));CHKERRQ(ierr); 2673 ierr = SNESDestroy(&(*snes)->pc);CHKERRQ(ierr); 2674 2675 /* if memory was published with AMS then destroy it */ 2676 ierr = PetscObjectDepublish((*snes));CHKERRQ(ierr); 2677 if ((*snes)->ops->destroy) {ierr = (*((*snes))->ops->destroy)((*snes));CHKERRQ(ierr);} 2678 2679 ierr = DMDestroy(&(*snes)->dm);CHKERRQ(ierr); 2680 ierr = KSPDestroy(&(*snes)->ksp);CHKERRQ(ierr); 2681 ierr = SNESLineSearchDestroy(&(*snes)->linesearch);CHKERRQ(ierr); 2682 2683 ierr = PetscFree((*snes)->kspconvctx);CHKERRQ(ierr); 2684 if ((*snes)->ops->convergeddestroy) { 2685 ierr = (*(*snes)->ops->convergeddestroy)((*snes)->cnvP);CHKERRQ(ierr); 2686 } 2687 if ((*snes)->conv_malloc) { 2688 ierr = PetscFree((*snes)->conv_hist);CHKERRQ(ierr); 2689 ierr = PetscFree((*snes)->conv_hist_its);CHKERRQ(ierr); 2690 } 2691 ierr = SNESMonitorCancel((*snes));CHKERRQ(ierr); 2692 ierr = PetscHeaderDestroy(snes);CHKERRQ(ierr); 2693 PetscFunctionReturn(0); 2694 } 2695 2696 /* ----------- Routines to set solver parameters ---------- */ 2697 2698 #undef __FUNCT__ 2699 #define __FUNCT__ "SNESSetLagPreconditioner" 2700 /*@ 2701 SNESSetLagPreconditioner - Determines when the preconditioner is rebuilt in the nonlinear solve. 2702 2703 Logically Collective on SNES 2704 2705 Input Parameters: 2706 + snes - the SNES context 2707 - lag - -1 indicates NEVER rebuild, 1 means rebuild every time the Jacobian is computed within a single nonlinear solve, 2 means every second time 2708 the Jacobian is built etc. -2 indicates rebuild preconditioner at next chance but then never rebuild after that 2709 2710 Options Database Keys: 2711 . -snes_lag_preconditioner <lag> 2712 2713 Notes: 2714 The default is 1 2715 The preconditioner is ALWAYS built in the first iteration of a nonlinear solve unless lag is -1 2716 If -1 is used before the very first nonlinear solve the preconditioner is still built because there is no previous preconditioner to use 2717 2718 Level: intermediate 2719 2720 .keywords: SNES, nonlinear, set, convergence, tolerances 2721 2722 .seealso: SNESSetTrustRegionTolerance(), SNESGetLagPreconditioner(), SNESSetLagJacobian(), SNESGetLagJacobian() 2723 2724 @*/ 2725 PetscErrorCode SNESSetLagPreconditioner(SNES snes,PetscInt lag) 2726 { 2727 PetscFunctionBegin; 2728 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2729 if (lag < -2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Lag must be -2, -1, 1 or greater"); 2730 if (!lag) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Lag cannot be 0"); 2731 PetscValidLogicalCollectiveInt(snes,lag,2); 2732 snes->lagpreconditioner = lag; 2733 PetscFunctionReturn(0); 2734 } 2735 2736 #undef __FUNCT__ 2737 #define __FUNCT__ "SNESSetGridSequence" 2738 /*@ 2739 SNESSetGridSequence - sets the number of steps of grid sequencing that SNES does 2740 2741 Logically Collective on SNES 2742 2743 Input Parameters: 2744 + snes - the SNES context 2745 - steps - the number of refinements to do, defaults to 0 2746 2747 Options Database Keys: 2748 . -snes_grid_sequence <steps> 2749 2750 Level: intermediate 2751 2752 Notes: 2753 Use SNESGetSolution() to extract the fine grid solution after grid sequencing. 2754 2755 .keywords: SNES, nonlinear, set, convergence, tolerances 2756 2757 .seealso: SNESSetTrustRegionTolerance(), SNESGetLagPreconditioner(), SNESSetLagJacobian(), SNESGetLagJacobian() 2758 2759 @*/ 2760 PetscErrorCode SNESSetGridSequence(SNES snes,PetscInt steps) 2761 { 2762 PetscFunctionBegin; 2763 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2764 PetscValidLogicalCollectiveInt(snes,steps,2); 2765 snes->gridsequence = steps; 2766 PetscFunctionReturn(0); 2767 } 2768 2769 #undef __FUNCT__ 2770 #define __FUNCT__ "SNESGetLagPreconditioner" 2771 /*@ 2772 SNESGetLagPreconditioner - Indicates how often the preconditioner is rebuilt 2773 2774 Not Collective 2775 2776 Input Parameter: 2777 . snes - the SNES context 2778 2779 Output Parameter: 2780 . lag - -1 indicates NEVER rebuild, 1 means rebuild every time the Jacobian is computed within a single nonlinear solve, 2 means every second time 2781 the Jacobian is built etc. -2 indicates rebuild preconditioner at next chance but then never rebuild after that 2782 2783 Options Database Keys: 2784 . -snes_lag_preconditioner <lag> 2785 2786 Notes: 2787 The default is 1 2788 The preconditioner is ALWAYS built in the first iteration of a nonlinear solve unless lag is -1 2789 2790 Level: intermediate 2791 2792 .keywords: SNES, nonlinear, set, convergence, tolerances 2793 2794 .seealso: SNESSetTrustRegionTolerance(), SNESSetLagPreconditioner() 2795 2796 @*/ 2797 PetscErrorCode SNESGetLagPreconditioner(SNES snes,PetscInt *lag) 2798 { 2799 PetscFunctionBegin; 2800 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2801 *lag = snes->lagpreconditioner; 2802 PetscFunctionReturn(0); 2803 } 2804 2805 #undef __FUNCT__ 2806 #define __FUNCT__ "SNESSetLagJacobian" 2807 /*@ 2808 SNESSetLagJacobian - Determines when the Jacobian is rebuilt in the nonlinear solve. See SNESSetLagPreconditioner() for determining how 2809 often the preconditioner is rebuilt. 2810 2811 Logically Collective on SNES 2812 2813 Input Parameters: 2814 + snes - the SNES context 2815 - lag - -1 indicates NEVER rebuild, 1 means rebuild every time the Jacobian is computed within a single nonlinear solve, 2 means every second time 2816 the Jacobian is built etc. -2 means rebuild at next chance but then never again 2817 2818 Options Database Keys: 2819 . -snes_lag_jacobian <lag> 2820 2821 Notes: 2822 The default is 1 2823 The Jacobian is ALWAYS built in the first iteration of a nonlinear solve unless lag is -1 2824 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 2825 at the next Newton step but never again (unless it is reset to another value) 2826 2827 Level: intermediate 2828 2829 .keywords: SNES, nonlinear, set, convergence, tolerances 2830 2831 .seealso: SNESSetTrustRegionTolerance(), SNESGetLagPreconditioner(), SNESSetLagPreconditioner(), SNESGetLagJacobian() 2832 2833 @*/ 2834 PetscErrorCode SNESSetLagJacobian(SNES snes,PetscInt lag) 2835 { 2836 PetscFunctionBegin; 2837 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2838 if (lag < -2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Lag must be -2, -1, 1 or greater"); 2839 if (!lag) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Lag cannot be 0"); 2840 PetscValidLogicalCollectiveInt(snes,lag,2); 2841 snes->lagjacobian = lag; 2842 PetscFunctionReturn(0); 2843 } 2844 2845 #undef __FUNCT__ 2846 #define __FUNCT__ "SNESGetLagJacobian" 2847 /*@ 2848 SNESGetLagJacobian - Indicates how often the Jacobian is rebuilt. See SNESGetLagPreconditioner() to determine when the preconditioner is rebuilt 2849 2850 Not Collective 2851 2852 Input Parameter: 2853 . snes - the SNES context 2854 2855 Output Parameter: 2856 . lag - -1 indicates NEVER rebuild, 1 means rebuild every time the Jacobian is computed within a single nonlinear solve, 2 means every second time 2857 the Jacobian is built etc. 2858 2859 Options Database Keys: 2860 . -snes_lag_jacobian <lag> 2861 2862 Notes: 2863 The default is 1 2864 The jacobian is ALWAYS built in the first iteration of a nonlinear solve unless lag is -1 2865 2866 Level: intermediate 2867 2868 .keywords: SNES, nonlinear, set, convergence, tolerances 2869 2870 .seealso: SNESSetTrustRegionTolerance(), SNESSetLagJacobian(), SNESSetLagPreconditioner(), SNESGetLagPreconditioner() 2871 2872 @*/ 2873 PetscErrorCode SNESGetLagJacobian(SNES snes,PetscInt *lag) 2874 { 2875 PetscFunctionBegin; 2876 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2877 *lag = snes->lagjacobian; 2878 PetscFunctionReturn(0); 2879 } 2880 2881 #undef __FUNCT__ 2882 #define __FUNCT__ "SNESSetTolerances" 2883 /*@ 2884 SNESSetTolerances - Sets various parameters used in convergence tests. 2885 2886 Logically Collective on SNES 2887 2888 Input Parameters: 2889 + snes - the SNES context 2890 . abstol - absolute convergence tolerance 2891 . rtol - relative convergence tolerance 2892 . stol - convergence tolerance in terms of the norm 2893 of the change in the solution between steps 2894 . maxit - maximum number of iterations 2895 - maxf - maximum number of function evaluations 2896 2897 Options Database Keys: 2898 + -snes_atol <abstol> - Sets abstol 2899 . -snes_rtol <rtol> - Sets rtol 2900 . -snes_stol <stol> - Sets stol 2901 . -snes_max_it <maxit> - Sets maxit 2902 - -snes_max_funcs <maxf> - Sets maxf 2903 2904 Notes: 2905 The default maximum number of iterations is 50. 2906 The default maximum number of function evaluations is 1000. 2907 2908 Level: intermediate 2909 2910 .keywords: SNES, nonlinear, set, convergence, tolerances 2911 2912 .seealso: SNESSetTrustRegionTolerance() 2913 @*/ 2914 PetscErrorCode SNESSetTolerances(SNES snes,PetscReal abstol,PetscReal rtol,PetscReal stol,PetscInt maxit,PetscInt maxf) 2915 { 2916 PetscFunctionBegin; 2917 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2918 PetscValidLogicalCollectiveReal(snes,abstol,2); 2919 PetscValidLogicalCollectiveReal(snes,rtol,3); 2920 PetscValidLogicalCollectiveReal(snes,stol,4); 2921 PetscValidLogicalCollectiveInt(snes,maxit,5); 2922 PetscValidLogicalCollectiveInt(snes,maxf,6); 2923 2924 if (abstol != PETSC_DEFAULT) { 2925 if (abstol < 0.0) SETERRQ1(((PetscObject)snes)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Absolute tolerance %G must be non-negative",abstol); 2926 snes->abstol = abstol; 2927 } 2928 if (rtol != PETSC_DEFAULT) { 2929 if (rtol < 0.0 || 1.0 <= rtol) SETERRQ1(((PetscObject)snes)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Relative tolerance %G must be non-negative and less than 1.0",rtol); 2930 snes->rtol = rtol; 2931 } 2932 if (stol != PETSC_DEFAULT) { 2933 if (stol < 0.0) SETERRQ1(((PetscObject)snes)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Step tolerance %G must be non-negative",stol); 2934 snes->stol = stol; 2935 } 2936 if (maxit != PETSC_DEFAULT) { 2937 if (maxit < 0) SETERRQ1(((PetscObject)snes)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Maximum number of iterations %D must be non-negative",maxit); 2938 snes->max_its = maxit; 2939 } 2940 if (maxf != PETSC_DEFAULT) { 2941 if (maxf < 0) SETERRQ1(((PetscObject)snes)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Maximum number of function evaluations %D must be non-negative",maxf); 2942 snes->max_funcs = maxf; 2943 } 2944 snes->tolerancesset = PETSC_TRUE; 2945 PetscFunctionReturn(0); 2946 } 2947 2948 #undef __FUNCT__ 2949 #define __FUNCT__ "SNESGetTolerances" 2950 /*@ 2951 SNESGetTolerances - Gets various parameters used in convergence tests. 2952 2953 Not Collective 2954 2955 Input Parameters: 2956 + snes - the SNES context 2957 . atol - absolute convergence tolerance 2958 . rtol - relative convergence tolerance 2959 . stol - convergence tolerance in terms of the norm 2960 of the change in the solution between steps 2961 . maxit - maximum number of iterations 2962 - maxf - maximum number of function evaluations 2963 2964 Notes: 2965 The user can specify PETSC_NULL for any parameter that is not needed. 2966 2967 Level: intermediate 2968 2969 .keywords: SNES, nonlinear, get, convergence, tolerances 2970 2971 .seealso: SNESSetTolerances() 2972 @*/ 2973 PetscErrorCode SNESGetTolerances(SNES snes,PetscReal *atol,PetscReal *rtol,PetscReal *stol,PetscInt *maxit,PetscInt *maxf) 2974 { 2975 PetscFunctionBegin; 2976 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2977 if (atol) *atol = snes->abstol; 2978 if (rtol) *rtol = snes->rtol; 2979 if (stol) *stol = snes->stol; 2980 if (maxit) *maxit = snes->max_its; 2981 if (maxf) *maxf = snes->max_funcs; 2982 PetscFunctionReturn(0); 2983 } 2984 2985 #undef __FUNCT__ 2986 #define __FUNCT__ "SNESSetTrustRegionTolerance" 2987 /*@ 2988 SNESSetTrustRegionTolerance - Sets the trust region parameter tolerance. 2989 2990 Logically Collective on SNES 2991 2992 Input Parameters: 2993 + snes - the SNES context 2994 - tol - tolerance 2995 2996 Options Database Key: 2997 . -snes_trtol <tol> - Sets tol 2998 2999 Level: intermediate 3000 3001 .keywords: SNES, nonlinear, set, trust region, tolerance 3002 3003 .seealso: SNESSetTolerances() 3004 @*/ 3005 PetscErrorCode SNESSetTrustRegionTolerance(SNES snes,PetscReal tol) 3006 { 3007 PetscFunctionBegin; 3008 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3009 PetscValidLogicalCollectiveReal(snes,tol,2); 3010 snes->deltatol = tol; 3011 PetscFunctionReturn(0); 3012 } 3013 3014 /* 3015 Duplicate the lg monitors for SNES from KSP; for some reason with 3016 dynamic libraries things don't work under Sun4 if we just use 3017 macros instead of functions 3018 */ 3019 #undef __FUNCT__ 3020 #define __FUNCT__ "SNESMonitorLGResidualNorm" 3021 PetscErrorCode SNESMonitorLGResidualNorm(SNES snes,PetscInt it,PetscReal norm,void *ctx) 3022 { 3023 PetscErrorCode ierr; 3024 3025 PetscFunctionBegin; 3026 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3027 ierr = KSPMonitorLGResidualNorm((KSP)snes,it,norm,ctx);CHKERRQ(ierr); 3028 PetscFunctionReturn(0); 3029 } 3030 3031 #undef __FUNCT__ 3032 #define __FUNCT__ "SNESMonitorLGCreate" 3033 PetscErrorCode SNESMonitorLGCreate(const char host[],const char label[],int x,int y,int m,int n,PetscDrawLG *draw) 3034 { 3035 PetscErrorCode ierr; 3036 3037 PetscFunctionBegin; 3038 ierr = KSPMonitorLGResidualNormCreate(host,label,x,y,m,n,draw);CHKERRQ(ierr); 3039 PetscFunctionReturn(0); 3040 } 3041 3042 #undef __FUNCT__ 3043 #define __FUNCT__ "SNESMonitorLGDestroy" 3044 PetscErrorCode SNESMonitorLGDestroy(PetscDrawLG *draw) 3045 { 3046 PetscErrorCode ierr; 3047 3048 PetscFunctionBegin; 3049 ierr = KSPMonitorLGResidualNormDestroy(draw);CHKERRQ(ierr); 3050 PetscFunctionReturn(0); 3051 } 3052 3053 extern PetscErrorCode SNESMonitorRange_Private(SNES,PetscInt,PetscReal*); 3054 #undef __FUNCT__ 3055 #define __FUNCT__ "SNESMonitorLGRange" 3056 PetscErrorCode SNESMonitorLGRange(SNES snes,PetscInt n,PetscReal rnorm,void *monctx) 3057 { 3058 PetscDrawLG lg; 3059 PetscErrorCode ierr; 3060 PetscReal x,y,per; 3061 PetscViewer v = (PetscViewer)monctx; 3062 static PetscReal prev; /* should be in the context */ 3063 PetscDraw draw; 3064 3065 PetscFunctionBegin; 3066 ierr = PetscViewerDrawGetDrawLG(v,0,&lg);CHKERRQ(ierr); 3067 if (!n) {ierr = PetscDrawLGReset(lg);CHKERRQ(ierr);} 3068 ierr = PetscDrawLGGetDraw(lg,&draw);CHKERRQ(ierr); 3069 ierr = PetscDrawSetTitle(draw,"Residual norm");CHKERRQ(ierr); 3070 x = (PetscReal) n; 3071 if (rnorm > 0.0) y = log10(rnorm); else y = -15.0; 3072 ierr = PetscDrawLGAddPoint(lg,&x,&y);CHKERRQ(ierr); 3073 if (n < 20 || !(n % 5)) { 3074 ierr = PetscDrawLGDraw(lg);CHKERRQ(ierr); 3075 } 3076 3077 ierr = PetscViewerDrawGetDrawLG(v,1,&lg);CHKERRQ(ierr); 3078 if (!n) {ierr = PetscDrawLGReset(lg);CHKERRQ(ierr);} 3079 ierr = PetscDrawLGGetDraw(lg,&draw);CHKERRQ(ierr); 3080 ierr = PetscDrawSetTitle(draw,"% elemts > .2*max elemt");CHKERRQ(ierr); 3081 ierr = SNESMonitorRange_Private(snes,n,&per);CHKERRQ(ierr); 3082 x = (PetscReal) n; 3083 y = 100.0*per; 3084 ierr = PetscDrawLGAddPoint(lg,&x,&y);CHKERRQ(ierr); 3085 if (n < 20 || !(n % 5)) { 3086 ierr = PetscDrawLGDraw(lg);CHKERRQ(ierr); 3087 } 3088 3089 ierr = PetscViewerDrawGetDrawLG(v,2,&lg);CHKERRQ(ierr); 3090 if (!n) {prev = rnorm;ierr = PetscDrawLGReset(lg);CHKERRQ(ierr);} 3091 ierr = PetscDrawLGGetDraw(lg,&draw);CHKERRQ(ierr); 3092 ierr = PetscDrawSetTitle(draw,"(norm -oldnorm)/oldnorm");CHKERRQ(ierr); 3093 x = (PetscReal) n; 3094 y = (prev - rnorm)/prev; 3095 ierr = PetscDrawLGAddPoint(lg,&x,&y);CHKERRQ(ierr); 3096 if (n < 20 || !(n % 5)) { 3097 ierr = PetscDrawLGDraw(lg);CHKERRQ(ierr); 3098 } 3099 3100 ierr = PetscViewerDrawGetDrawLG(v,3,&lg);CHKERRQ(ierr); 3101 if (!n) {ierr = PetscDrawLGReset(lg);CHKERRQ(ierr);} 3102 ierr = PetscDrawLGGetDraw(lg,&draw);CHKERRQ(ierr); 3103 ierr = PetscDrawSetTitle(draw,"(norm -oldnorm)/oldnorm*(% > .2 max)");CHKERRQ(ierr); 3104 x = (PetscReal) n; 3105 y = (prev - rnorm)/(prev*per); 3106 if (n > 2) { /*skip initial crazy value */ 3107 ierr = PetscDrawLGAddPoint(lg,&x,&y);CHKERRQ(ierr); 3108 } 3109 if (n < 20 || !(n % 5)) { 3110 ierr = PetscDrawLGDraw(lg);CHKERRQ(ierr); 3111 } 3112 prev = rnorm; 3113 PetscFunctionReturn(0); 3114 } 3115 3116 #undef __FUNCT__ 3117 #define __FUNCT__ "SNESMonitor" 3118 /*@ 3119 SNESMonitor - runs the user provided monitor routines, if they exist 3120 3121 Collective on SNES 3122 3123 Input Parameters: 3124 + snes - nonlinear solver context obtained from SNESCreate() 3125 . iter - iteration number 3126 - rnorm - relative norm of the residual 3127 3128 Notes: 3129 This routine is called by the SNES implementations. 3130 It does not typically need to be called by the user. 3131 3132 Level: developer 3133 3134 .seealso: SNESMonitorSet() 3135 @*/ 3136 PetscErrorCode SNESMonitor(SNES snes,PetscInt iter,PetscReal rnorm) 3137 { 3138 PetscErrorCode ierr; 3139 PetscInt i,n = snes->numbermonitors; 3140 3141 PetscFunctionBegin; 3142 for (i=0; i<n; i++) { 3143 ierr = (*snes->monitor[i])(snes,iter,rnorm,snes->monitorcontext[i]);CHKERRQ(ierr); 3144 } 3145 PetscFunctionReturn(0); 3146 } 3147 3148 /* ------------ Routines to set performance monitoring options ----------- */ 3149 3150 #undef __FUNCT__ 3151 #define __FUNCT__ "SNESMonitorSet" 3152 /*@C 3153 SNESMonitorSet - Sets an ADDITIONAL function that is to be used at every 3154 iteration of the nonlinear solver to display the iteration's 3155 progress. 3156 3157 Logically Collective on SNES 3158 3159 Input Parameters: 3160 + snes - the SNES context 3161 . func - monitoring routine 3162 . mctx - [optional] user-defined context for private data for the 3163 monitor routine (use PETSC_NULL if no context is desired) 3164 - monitordestroy - [optional] routine that frees monitor context 3165 (may be PETSC_NULL) 3166 3167 Calling sequence of func: 3168 $ int func(SNES snes,PetscInt its, PetscReal norm,void *mctx) 3169 3170 + snes - the SNES context 3171 . its - iteration number 3172 . norm - 2-norm function value (may be estimated) 3173 - mctx - [optional] monitoring context 3174 3175 Options Database Keys: 3176 + -snes_monitor - sets SNESMonitorDefault() 3177 . -snes_monitor_lg_residualnorm - sets line graph monitor, 3178 uses SNESMonitorLGCreate() 3179 - -snes_monitor_cancel - cancels all monitors that have 3180 been hardwired into a code by 3181 calls to SNESMonitorSet(), but 3182 does not cancel those set via 3183 the options database. 3184 3185 Notes: 3186 Several different monitoring routines may be set by calling 3187 SNESMonitorSet() multiple times; all will be called in the 3188 order in which they were set. 3189 3190 Fortran notes: Only a single monitor function can be set for each SNES object 3191 3192 Level: intermediate 3193 3194 .keywords: SNES, nonlinear, set, monitor 3195 3196 .seealso: SNESMonitorDefault(), SNESMonitorCancel() 3197 @*/ 3198 PetscErrorCode SNESMonitorSet(SNES snes,PetscErrorCode (*monitor)(SNES,PetscInt,PetscReal,void*),void *mctx,PetscErrorCode (*monitordestroy)(void**)) 3199 { 3200 PetscInt i; 3201 PetscErrorCode ierr; 3202 3203 PetscFunctionBegin; 3204 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3205 if (snes->numbermonitors >= MAXSNESMONITORS) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Too many monitors set"); 3206 for (i=0; i<snes->numbermonitors;i++) { 3207 if (monitor == snes->monitor[i] && monitordestroy == snes->monitordestroy[i] && mctx == snes->monitorcontext[i]) { 3208 if (monitordestroy) { 3209 ierr = (*monitordestroy)(&mctx);CHKERRQ(ierr); 3210 } 3211 PetscFunctionReturn(0); 3212 } 3213 } 3214 snes->monitor[snes->numbermonitors] = monitor; 3215 snes->monitordestroy[snes->numbermonitors] = monitordestroy; 3216 snes->monitorcontext[snes->numbermonitors++] = (void*)mctx; 3217 PetscFunctionReturn(0); 3218 } 3219 3220 #undef __FUNCT__ 3221 #define __FUNCT__ "SNESMonitorCancel" 3222 /*@C 3223 SNESMonitorCancel - Clears all the monitor functions for a SNES object. 3224 3225 Logically Collective on SNES 3226 3227 Input Parameters: 3228 . snes - the SNES context 3229 3230 Options Database Key: 3231 . -snes_monitor_cancel - cancels all monitors that have been hardwired 3232 into a code by calls to SNESMonitorSet(), but does not cancel those 3233 set via the options database 3234 3235 Notes: 3236 There is no way to clear one specific monitor from a SNES object. 3237 3238 Level: intermediate 3239 3240 .keywords: SNES, nonlinear, set, monitor 3241 3242 .seealso: SNESMonitorDefault(), SNESMonitorSet() 3243 @*/ 3244 PetscErrorCode SNESMonitorCancel(SNES snes) 3245 { 3246 PetscErrorCode ierr; 3247 PetscInt i; 3248 3249 PetscFunctionBegin; 3250 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3251 for (i=0; i<snes->numbermonitors; i++) { 3252 if (snes->monitordestroy[i]) { 3253 ierr = (*snes->monitordestroy[i])(&snes->monitorcontext[i]);CHKERRQ(ierr); 3254 } 3255 } 3256 snes->numbermonitors = 0; 3257 PetscFunctionReturn(0); 3258 } 3259 3260 #undef __FUNCT__ 3261 #define __FUNCT__ "SNESSetConvergenceTest" 3262 /*@C 3263 SNESSetConvergenceTest - Sets the function that is to be used 3264 to test for convergence of the nonlinear iterative solution. 3265 3266 Logically Collective on SNES 3267 3268 Input Parameters: 3269 + snes - the SNES context 3270 . func - routine to test for convergence 3271 . cctx - [optional] context for private data for the convergence routine (may be PETSC_NULL) 3272 - destroy - [optional] destructor for the context (may be PETSC_NULL; PETSC_NULL_FUNCTION in Fortran) 3273 3274 Calling sequence of func: 3275 $ PetscErrorCode func (SNES snes,PetscInt it,PetscReal xnorm,PetscReal gnorm,PetscReal f,SNESConvergedReason *reason,void *cctx) 3276 3277 + snes - the SNES context 3278 . it - current iteration (0 is the first and is before any Newton step) 3279 . cctx - [optional] convergence context 3280 . reason - reason for convergence/divergence 3281 . xnorm - 2-norm of current iterate 3282 . gnorm - 2-norm of current step 3283 - f - 2-norm of function 3284 3285 Level: advanced 3286 3287 .keywords: SNES, nonlinear, set, convergence, test 3288 3289 .seealso: SNESDefaultConverged(), SNESSkipConverged() 3290 @*/ 3291 PetscErrorCode SNESSetConvergenceTest(SNES snes,PetscErrorCode (*func)(SNES,PetscInt,PetscReal,PetscReal,PetscReal,SNESConvergedReason*,void*),void *cctx,PetscErrorCode (*destroy)(void*)) 3292 { 3293 PetscErrorCode ierr; 3294 3295 PetscFunctionBegin; 3296 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3297 if (!func) func = SNESSkipConverged; 3298 if (snes->ops->convergeddestroy) { 3299 ierr = (*snes->ops->convergeddestroy)(snes->cnvP);CHKERRQ(ierr); 3300 } 3301 snes->ops->converged = func; 3302 snes->ops->convergeddestroy = destroy; 3303 snes->cnvP = cctx; 3304 PetscFunctionReturn(0); 3305 } 3306 3307 #undef __FUNCT__ 3308 #define __FUNCT__ "SNESGetConvergedReason" 3309 /*@ 3310 SNESGetConvergedReason - Gets the reason the SNES iteration was stopped. 3311 3312 Not Collective 3313 3314 Input Parameter: 3315 . snes - the SNES context 3316 3317 Output Parameter: 3318 . reason - negative value indicates diverged, positive value converged, see SNESConvergedReason or the 3319 manual pages for the individual convergence tests for complete lists 3320 3321 Level: intermediate 3322 3323 Notes: Can only be called after the call the SNESSolve() is complete. 3324 3325 .keywords: SNES, nonlinear, set, convergence, test 3326 3327 .seealso: SNESSetConvergenceTest(), SNESConvergedReason 3328 @*/ 3329 PetscErrorCode SNESGetConvergedReason(SNES snes,SNESConvergedReason *reason) 3330 { 3331 PetscFunctionBegin; 3332 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3333 PetscValidPointer(reason,2); 3334 *reason = snes->reason; 3335 PetscFunctionReturn(0); 3336 } 3337 3338 #undef __FUNCT__ 3339 #define __FUNCT__ "SNESSetConvergenceHistory" 3340 /*@ 3341 SNESSetConvergenceHistory - Sets the array used to hold the convergence history. 3342 3343 Logically Collective on SNES 3344 3345 Input Parameters: 3346 + snes - iterative context obtained from SNESCreate() 3347 . a - array to hold history, this array will contain the function norms computed at each step 3348 . its - integer array holds the number of linear iterations for each solve. 3349 . na - size of a and its 3350 - reset - PETSC_TRUE indicates each new nonlinear solve resets the history counter to zero, 3351 else it continues storing new values for new nonlinear solves after the old ones 3352 3353 Notes: 3354 If 'a' and 'its' are PETSC_NULL then space is allocated for the history. If 'na' PETSC_DECIDE or PETSC_DEFAULT then a 3355 default array of length 10000 is allocated. 3356 3357 This routine is useful, e.g., when running a code for purposes 3358 of accurate performance monitoring, when no I/O should be done 3359 during the section of code that is being timed. 3360 3361 Level: intermediate 3362 3363 .keywords: SNES, set, convergence, history 3364 3365 .seealso: SNESGetConvergenceHistory() 3366 3367 @*/ 3368 PetscErrorCode SNESSetConvergenceHistory(SNES snes,PetscReal a[],PetscInt its[],PetscInt na,PetscBool reset) 3369 { 3370 PetscErrorCode ierr; 3371 3372 PetscFunctionBegin; 3373 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3374 if (na) PetscValidScalarPointer(a,2); 3375 if (its) PetscValidIntPointer(its,3); 3376 if (na == PETSC_DECIDE || na == PETSC_DEFAULT || !a) { 3377 if (na == PETSC_DECIDE || na == PETSC_DEFAULT) na = 1000; 3378 ierr = PetscMalloc(na*sizeof(PetscReal),&a);CHKERRQ(ierr); 3379 ierr = PetscMalloc(na*sizeof(PetscInt),&its);CHKERRQ(ierr); 3380 snes->conv_malloc = PETSC_TRUE; 3381 } 3382 snes->conv_hist = a; 3383 snes->conv_hist_its = its; 3384 snes->conv_hist_max = na; 3385 snes->conv_hist_len = 0; 3386 snes->conv_hist_reset = reset; 3387 PetscFunctionReturn(0); 3388 } 3389 3390 #if defined(PETSC_HAVE_MATLAB_ENGINE) 3391 #include <engine.h> /* MATLAB include file */ 3392 #include <mex.h> /* MATLAB include file */ 3393 EXTERN_C_BEGIN 3394 #undef __FUNCT__ 3395 #define __FUNCT__ "SNESGetConvergenceHistoryMatlab" 3396 mxArray *SNESGetConvergenceHistoryMatlab(SNES snes) 3397 { 3398 mxArray *mat; 3399 PetscInt i; 3400 PetscReal *ar; 3401 3402 PetscFunctionBegin; 3403 mat = mxCreateDoubleMatrix(snes->conv_hist_len,1,mxREAL); 3404 ar = (PetscReal*) mxGetData(mat); 3405 for (i=0; i<snes->conv_hist_len; i++) { 3406 ar[i] = snes->conv_hist[i]; 3407 } 3408 PetscFunctionReturn(mat); 3409 } 3410 EXTERN_C_END 3411 #endif 3412 3413 3414 #undef __FUNCT__ 3415 #define __FUNCT__ "SNESGetConvergenceHistory" 3416 /*@C 3417 SNESGetConvergenceHistory - Gets the array used to hold the convergence history. 3418 3419 Not Collective 3420 3421 Input Parameter: 3422 . snes - iterative context obtained from SNESCreate() 3423 3424 Output Parameters: 3425 . a - array to hold history 3426 . its - integer array holds the number of linear iterations (or 3427 negative if not converged) for each solve. 3428 - na - size of a and its 3429 3430 Notes: 3431 The calling sequence for this routine in Fortran is 3432 $ call SNESGetConvergenceHistory(SNES snes, integer na, integer ierr) 3433 3434 This routine is useful, e.g., when running a code for purposes 3435 of accurate performance monitoring, when no I/O should be done 3436 during the section of code that is being timed. 3437 3438 Level: intermediate 3439 3440 .keywords: SNES, get, convergence, history 3441 3442 .seealso: SNESSetConvergencHistory() 3443 3444 @*/ 3445 PetscErrorCode SNESGetConvergenceHistory(SNES snes,PetscReal *a[],PetscInt *its[],PetscInt *na) 3446 { 3447 PetscFunctionBegin; 3448 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3449 if (a) *a = snes->conv_hist; 3450 if (its) *its = snes->conv_hist_its; 3451 if (na) *na = snes->conv_hist_len; 3452 PetscFunctionReturn(0); 3453 } 3454 3455 #undef __FUNCT__ 3456 #define __FUNCT__ "SNESSetUpdate" 3457 /*@C 3458 SNESSetUpdate - Sets the general-purpose update function called 3459 at the beginning of every iteration of the nonlinear solve. Specifically 3460 it is called just before the Jacobian is "evaluated". 3461 3462 Logically Collective on SNES 3463 3464 Input Parameters: 3465 . snes - The nonlinear solver context 3466 . func - The function 3467 3468 Calling sequence of func: 3469 . func (SNES snes, PetscInt step); 3470 3471 . step - The current step of the iteration 3472 3473 Level: advanced 3474 3475 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() 3476 This is not used by most users. 3477 3478 .keywords: SNES, update 3479 3480 .seealso SNESDefaultUpdate(), SNESSetJacobian(), SNESSolve() 3481 @*/ 3482 PetscErrorCode SNESSetUpdate(SNES snes, PetscErrorCode (*func)(SNES, PetscInt)) 3483 { 3484 PetscFunctionBegin; 3485 PetscValidHeaderSpecific(snes, SNES_CLASSID,1); 3486 snes->ops->update = func; 3487 PetscFunctionReturn(0); 3488 } 3489 3490 #undef __FUNCT__ 3491 #define __FUNCT__ "SNESDefaultUpdate" 3492 /*@ 3493 SNESDefaultUpdate - The default update function which does nothing. 3494 3495 Not collective 3496 3497 Input Parameters: 3498 . snes - The nonlinear solver context 3499 . step - The current step of the iteration 3500 3501 Level: intermediate 3502 3503 .keywords: SNES, update 3504 .seealso SNESSetUpdate(), SNESDefaultRhsBC(), SNESDefaultShortolutionBC() 3505 @*/ 3506 PetscErrorCode SNESDefaultUpdate(SNES snes, PetscInt step) 3507 { 3508 PetscFunctionBegin; 3509 PetscFunctionReturn(0); 3510 } 3511 3512 #undef __FUNCT__ 3513 #define __FUNCT__ "SNESScaleStep_Private" 3514 /* 3515 SNESScaleStep_Private - Scales a step so that its length is less than the 3516 positive parameter delta. 3517 3518 Input Parameters: 3519 + snes - the SNES context 3520 . y - approximate solution of linear system 3521 . fnorm - 2-norm of current function 3522 - delta - trust region size 3523 3524 Output Parameters: 3525 + gpnorm - predicted function norm at the new point, assuming local 3526 linearization. The value is zero if the step lies within the trust 3527 region, and exceeds zero otherwise. 3528 - ynorm - 2-norm of the step 3529 3530 Note: 3531 For non-trust region methods such as SNESLS, the parameter delta 3532 is set to be the maximum allowable step size. 3533 3534 .keywords: SNES, nonlinear, scale, step 3535 */ 3536 PetscErrorCode SNESScaleStep_Private(SNES snes,Vec y,PetscReal *fnorm,PetscReal *delta,PetscReal *gpnorm,PetscReal *ynorm) 3537 { 3538 PetscReal nrm; 3539 PetscScalar cnorm; 3540 PetscErrorCode ierr; 3541 3542 PetscFunctionBegin; 3543 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3544 PetscValidHeaderSpecific(y,VEC_CLASSID,2); 3545 PetscCheckSameComm(snes,1,y,2); 3546 3547 ierr = VecNorm(y,NORM_2,&nrm);CHKERRQ(ierr); 3548 if (nrm > *delta) { 3549 nrm = *delta/nrm; 3550 *gpnorm = (1.0 - nrm)*(*fnorm); 3551 cnorm = nrm; 3552 ierr = VecScale(y,cnorm);CHKERRQ(ierr); 3553 *ynorm = *delta; 3554 } else { 3555 *gpnorm = 0.0; 3556 *ynorm = nrm; 3557 } 3558 PetscFunctionReturn(0); 3559 } 3560 3561 #undef __FUNCT__ 3562 #define __FUNCT__ "SNESSolve" 3563 /*@C 3564 SNESSolve - Solves a nonlinear system F(x) = b. 3565 Call SNESSolve() after calling SNESCreate() and optional routines of the form SNESSetXXX(). 3566 3567 Collective on SNES 3568 3569 Input Parameters: 3570 + snes - the SNES context 3571 . b - the constant part of the equation F(x) = b, or PETSC_NULL to use zero. 3572 - x - the solution vector. 3573 3574 Notes: 3575 The user should initialize the vector,x, with the initial guess 3576 for the nonlinear solve prior to calling SNESSolve. In particular, 3577 to employ an initial guess of zero, the user should explicitly set 3578 this vector to zero by calling VecSet(). 3579 3580 Level: beginner 3581 3582 .keywords: SNES, nonlinear, solve 3583 3584 .seealso: SNESCreate(), SNESDestroy(), SNESSetFunction(), SNESSetJacobian(), SNESSetGridSequence(), SNESGetSolution() 3585 @*/ 3586 PetscErrorCode SNESSolve(SNES snes,Vec b,Vec x) 3587 { 3588 PetscErrorCode ierr; 3589 PetscBool flg; 3590 char filename[PETSC_MAX_PATH_LEN]; 3591 PetscViewer viewer; 3592 PetscInt grid; 3593 Vec xcreated = PETSC_NULL; 3594 DM dm; 3595 3596 PetscFunctionBegin; 3597 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3598 if (x) PetscValidHeaderSpecific(x,VEC_CLASSID,3); 3599 if (x) PetscCheckSameComm(snes,1,x,3); 3600 if (b) PetscValidHeaderSpecific(b,VEC_CLASSID,2); 3601 if (b) PetscCheckSameComm(snes,1,b,2); 3602 3603 if (!x) { 3604 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 3605 ierr = DMCreateGlobalVector(dm,&xcreated);CHKERRQ(ierr); 3606 x = xcreated; 3607 } 3608 3609 for (grid=0; grid<snes->gridsequence; grid++) {ierr = PetscViewerASCIIPushTab(PETSC_VIEWER_STDOUT_(((PetscObject)snes)->comm));CHKERRQ(ierr);} 3610 for (grid=0; grid<snes->gridsequence+1; grid++) { 3611 3612 /* set solution vector */ 3613 if (!grid) {ierr = PetscObjectReference((PetscObject)x);CHKERRQ(ierr);} 3614 ierr = VecDestroy(&snes->vec_sol);CHKERRQ(ierr); 3615 snes->vec_sol = x; 3616 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 3617 3618 /* set affine vector if provided */ 3619 if (b) { ierr = PetscObjectReference((PetscObject)b);CHKERRQ(ierr); } 3620 ierr = VecDestroy(&snes->vec_rhs);CHKERRQ(ierr); 3621 snes->vec_rhs = b; 3622 3623 ierr = SNESSetUp(snes);CHKERRQ(ierr); 3624 3625 if (!grid) { 3626 if (snes->ops->computeinitialguess) { 3627 ierr = (*snes->ops->computeinitialguess)(snes,snes->vec_sol,snes->initialguessP);CHKERRQ(ierr); 3628 } 3629 } 3630 3631 if (snes->conv_hist_reset) snes->conv_hist_len = 0; 3632 snes->nfuncs = 0; snes->linear_its = 0; snes->numFailures = 0; 3633 3634 ierr = PetscLogEventBegin(SNES_Solve,snes,0,0,0);CHKERRQ(ierr); 3635 ierr = (*snes->ops->solve)(snes);CHKERRQ(ierr); 3636 ierr = PetscLogEventEnd(SNES_Solve,snes,0,0,0);CHKERRQ(ierr); 3637 if (snes->domainerror){ 3638 snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; 3639 snes->domainerror = PETSC_FALSE; 3640 } 3641 if (!snes->reason) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Internal error, solver returned without setting converged reason"); 3642 3643 flg = PETSC_FALSE; 3644 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_test_local_min",&flg,PETSC_NULL);CHKERRQ(ierr); 3645 if (flg && !PetscPreLoadingOn) { ierr = SNESTestLocalMin(snes);CHKERRQ(ierr); } 3646 if (snes->printreason) { 3647 ierr = PetscViewerASCIIAddTab(PETSC_VIEWER_STDOUT_(((PetscObject)snes)->comm),((PetscObject)snes)->tablevel);CHKERRQ(ierr); 3648 if (snes->reason > 0) { 3649 ierr = PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_(((PetscObject)snes)->comm),"Nonlinear solve converged due to %s iterations %D\n",SNESConvergedReasons[snes->reason],snes->iter);CHKERRQ(ierr); 3650 } else { 3651 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); 3652 } 3653 ierr = PetscViewerASCIISubtractTab(PETSC_VIEWER_STDOUT_(((PetscObject)snes)->comm),((PetscObject)snes)->tablevel);CHKERRQ(ierr); 3654 } 3655 3656 if (snes->errorifnotconverged && snes->reason < 0) SETERRQ(((PetscObject)snes)->comm,PETSC_ERR_NOT_CONVERGED,"SNESSolve has not converged"); 3657 if (grid < snes->gridsequence) { 3658 DM fine; 3659 Vec xnew; 3660 Mat interp; 3661 3662 ierr = DMRefine(snes->dm,((PetscObject)snes)->comm,&fine);CHKERRQ(ierr); 3663 if (!fine) SETERRQ(((PetscObject)snes)->comm,PETSC_ERR_ARG_INCOMP,"DMRefine() did not perform any refinement, cannot continue grid sequencing"); 3664 ierr = DMCreateInterpolation(snes->dm,fine,&interp,PETSC_NULL);CHKERRQ(ierr); 3665 ierr = DMCreateGlobalVector(fine,&xnew);CHKERRQ(ierr); 3666 ierr = MatInterpolate(interp,x,xnew);CHKERRQ(ierr); 3667 ierr = DMInterpolate(snes->dm,interp,fine);CHKERRQ(ierr); 3668 ierr = MatDestroy(&interp);CHKERRQ(ierr); 3669 x = xnew; 3670 3671 ierr = SNESReset(snes);CHKERRQ(ierr); 3672 ierr = SNESSetDM(snes,fine);CHKERRQ(ierr); 3673 ierr = DMDestroy(&fine);CHKERRQ(ierr); 3674 ierr = PetscViewerASCIIPopTab(PETSC_VIEWER_STDOUT_(((PetscObject)snes)->comm));CHKERRQ(ierr); 3675 } 3676 } 3677 /* monitoring and viewing */ 3678 flg = PETSC_FALSE; 3679 ierr = PetscOptionsGetString(((PetscObject)snes)->prefix,"-snes_view",filename,PETSC_MAX_PATH_LEN,&flg);CHKERRQ(ierr); 3680 if (flg && !PetscPreLoadingOn) { 3681 ierr = PetscViewerASCIIOpen(((PetscObject)snes)->comm,filename,&viewer);CHKERRQ(ierr); 3682 ierr = SNESView(snes,viewer);CHKERRQ(ierr); 3683 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 3684 } 3685 3686 flg = PETSC_FALSE; 3687 ierr = PetscOptionsGetString(((PetscObject)snes)->prefix,"-snes_view_solution_vtk",filename,PETSC_MAX_PATH_LEN,&flg);CHKERRQ(ierr); 3688 if (flg) { 3689 PetscViewer viewer; 3690 ierr = PetscViewerCreate(((PetscObject)snes)->comm,&viewer);CHKERRQ(ierr); 3691 ierr = PetscViewerSetType(viewer,PETSCVIEWERVTK);CHKERRQ(ierr); 3692 ierr = PetscViewerFileSetName(viewer,filename);CHKERRQ(ierr); 3693 ierr = VecView(snes->vec_sol,viewer);CHKERRQ(ierr); 3694 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 3695 } 3696 3697 ierr = VecDestroy(&xcreated);CHKERRQ(ierr); 3698 PetscFunctionReturn(0); 3699 } 3700 3701 /* --------- Internal routines for SNES Package --------- */ 3702 3703 #undef __FUNCT__ 3704 #define __FUNCT__ "SNESSetType" 3705 /*@C 3706 SNESSetType - Sets the method for the nonlinear solver. 3707 3708 Collective on SNES 3709 3710 Input Parameters: 3711 + snes - the SNES context 3712 - type - a known method 3713 3714 Options Database Key: 3715 . -snes_type <type> - Sets the method; use -help for a list 3716 of available methods (for instance, ls or tr) 3717 3718 Notes: 3719 See "petsc/include/petscsnes.h" for available methods (for instance) 3720 + SNESLS - Newton's method with line search 3721 (systems of nonlinear equations) 3722 . SNESTR - Newton's method with trust region 3723 (systems of nonlinear equations) 3724 3725 Normally, it is best to use the SNESSetFromOptions() command and then 3726 set the SNES solver type from the options database rather than by using 3727 this routine. Using the options database provides the user with 3728 maximum flexibility in evaluating the many nonlinear solvers. 3729 The SNESSetType() routine is provided for those situations where it 3730 is necessary to set the nonlinear solver independently of the command 3731 line or options database. This might be the case, for example, when 3732 the choice of solver changes during the execution of the program, 3733 and the user's application is taking responsibility for choosing the 3734 appropriate method. 3735 3736 Level: intermediate 3737 3738 .keywords: SNES, set, type 3739 3740 .seealso: SNESType, SNESCreate() 3741 3742 @*/ 3743 PetscErrorCode SNESSetType(SNES snes,SNESType type) 3744 { 3745 PetscErrorCode ierr,(*r)(SNES); 3746 PetscBool match; 3747 3748 PetscFunctionBegin; 3749 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3750 PetscValidCharPointer(type,2); 3751 3752 ierr = PetscObjectTypeCompare((PetscObject)snes,type,&match);CHKERRQ(ierr); 3753 if (match) PetscFunctionReturn(0); 3754 3755 ierr = PetscFListFind(SNESList,((PetscObject)snes)->comm,type,PETSC_TRUE,(void (**)(void)) &r);CHKERRQ(ierr); 3756 if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unable to find requested SNES type %s",type); 3757 /* Destroy the previous private SNES context */ 3758 if (snes->ops->destroy) { 3759 ierr = (*(snes)->ops->destroy)(snes);CHKERRQ(ierr); 3760 snes->ops->destroy = PETSC_NULL; 3761 } 3762 /* Reinitialize function pointers in SNESOps structure */ 3763 snes->ops->setup = 0; 3764 snes->ops->solve = 0; 3765 snes->ops->view = 0; 3766 snes->ops->setfromoptions = 0; 3767 snes->ops->destroy = 0; 3768 /* Call the SNESCreate_XXX routine for this particular Nonlinear solver */ 3769 snes->setupcalled = PETSC_FALSE; 3770 ierr = PetscObjectChangeTypeName((PetscObject)snes,type);CHKERRQ(ierr); 3771 ierr = (*r)(snes);CHKERRQ(ierr); 3772 #if defined(PETSC_HAVE_AMS) 3773 if (PetscAMSPublishAll) { 3774 ierr = PetscObjectAMSPublish((PetscObject)snes);CHKERRQ(ierr); 3775 } 3776 #endif 3777 PetscFunctionReturn(0); 3778 } 3779 3780 3781 /* --------------------------------------------------------------------- */ 3782 #undef __FUNCT__ 3783 #define __FUNCT__ "SNESRegisterDestroy" 3784 /*@ 3785 SNESRegisterDestroy - Frees the list of nonlinear solvers that were 3786 registered by SNESRegisterDynamic(). 3787 3788 Not Collective 3789 3790 Level: advanced 3791 3792 .keywords: SNES, nonlinear, register, destroy 3793 3794 .seealso: SNESRegisterAll(), SNESRegisterAll() 3795 @*/ 3796 PetscErrorCode SNESRegisterDestroy(void) 3797 { 3798 PetscErrorCode ierr; 3799 3800 PetscFunctionBegin; 3801 ierr = PetscFListDestroy(&SNESList);CHKERRQ(ierr); 3802 SNESRegisterAllCalled = PETSC_FALSE; 3803 PetscFunctionReturn(0); 3804 } 3805 3806 #undef __FUNCT__ 3807 #define __FUNCT__ "SNESGetType" 3808 /*@C 3809 SNESGetType - Gets the SNES method type and name (as a string). 3810 3811 Not Collective 3812 3813 Input Parameter: 3814 . snes - nonlinear solver context 3815 3816 Output Parameter: 3817 . type - SNES method (a character string) 3818 3819 Level: intermediate 3820 3821 .keywords: SNES, nonlinear, get, type, name 3822 @*/ 3823 PetscErrorCode SNESGetType(SNES snes,SNESType *type) 3824 { 3825 PetscFunctionBegin; 3826 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3827 PetscValidPointer(type,2); 3828 *type = ((PetscObject)snes)->type_name; 3829 PetscFunctionReturn(0); 3830 } 3831 3832 #undef __FUNCT__ 3833 #define __FUNCT__ "SNESGetSolution" 3834 /*@ 3835 SNESGetSolution - Returns the vector where the approximate solution is 3836 stored. This is the fine grid solution when using SNESSetGridSequence(). 3837 3838 Not Collective, but Vec is parallel if SNES is parallel 3839 3840 Input Parameter: 3841 . snes - the SNES context 3842 3843 Output Parameter: 3844 . x - the solution 3845 3846 Level: intermediate 3847 3848 .keywords: SNES, nonlinear, get, solution 3849 3850 .seealso: SNESGetSolutionUpdate(), SNESGetFunction() 3851 @*/ 3852 PetscErrorCode SNESGetSolution(SNES snes,Vec *x) 3853 { 3854 PetscFunctionBegin; 3855 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3856 PetscValidPointer(x,2); 3857 *x = snes->vec_sol; 3858 PetscFunctionReturn(0); 3859 } 3860 3861 #undef __FUNCT__ 3862 #define __FUNCT__ "SNESGetSolutionUpdate" 3863 /*@ 3864 SNESGetSolutionUpdate - Returns the vector where the solution update is 3865 stored. 3866 3867 Not Collective, but Vec is parallel if SNES is parallel 3868 3869 Input Parameter: 3870 . snes - the SNES context 3871 3872 Output Parameter: 3873 . x - the solution update 3874 3875 Level: advanced 3876 3877 .keywords: SNES, nonlinear, get, solution, update 3878 3879 .seealso: SNESGetSolution(), SNESGetFunction() 3880 @*/ 3881 PetscErrorCode SNESGetSolutionUpdate(SNES snes,Vec *x) 3882 { 3883 PetscFunctionBegin; 3884 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3885 PetscValidPointer(x,2); 3886 *x = snes->vec_sol_update; 3887 PetscFunctionReturn(0); 3888 } 3889 3890 #undef __FUNCT__ 3891 #define __FUNCT__ "SNESGetFunction" 3892 /*@C 3893 SNESGetFunction - Returns the vector where the function is stored. 3894 3895 Not Collective, but Vec is parallel if SNES is parallel. Collective if Vec is requested, but has not been created yet. 3896 3897 Input Parameter: 3898 . snes - the SNES context 3899 3900 Output Parameter: 3901 + r - the function (or PETSC_NULL) 3902 . func - the function (or PETSC_NULL) 3903 - ctx - the function context (or PETSC_NULL) 3904 3905 Level: advanced 3906 3907 .keywords: SNES, nonlinear, get, function 3908 3909 .seealso: SNESSetFunction(), SNESGetSolution() 3910 @*/ 3911 PetscErrorCode SNESGetFunction(SNES snes,Vec *r,PetscErrorCode (**func)(SNES,Vec,Vec,void*),void **ctx) 3912 { 3913 PetscErrorCode ierr; 3914 DM dm; 3915 3916 PetscFunctionBegin; 3917 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3918 if (r) { 3919 if (!snes->vec_func) { 3920 if (snes->vec_rhs) { 3921 ierr = VecDuplicate(snes->vec_rhs,&snes->vec_func);CHKERRQ(ierr); 3922 } else if (snes->vec_sol) { 3923 ierr = VecDuplicate(snes->vec_sol,&snes->vec_func);CHKERRQ(ierr); 3924 } else if (snes->dm) { 3925 ierr = DMCreateGlobalVector(snes->dm,&snes->vec_func);CHKERRQ(ierr); 3926 } 3927 } 3928 *r = snes->vec_func; 3929 } 3930 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 3931 ierr = DMSNESGetFunction(dm,func,ctx);CHKERRQ(ierr); 3932 PetscFunctionReturn(0); 3933 } 3934 3935 /*@C 3936 SNESGetGS - Returns the GS function and context. 3937 3938 Input Parameter: 3939 . snes - the SNES context 3940 3941 Output Parameter: 3942 + gsfunc - the function (or PETSC_NULL) 3943 - ctx - the function context (or PETSC_NULL) 3944 3945 Level: advanced 3946 3947 .keywords: SNES, nonlinear, get, function 3948 3949 .seealso: SNESSetGS(), SNESGetFunction() 3950 @*/ 3951 3952 #undef __FUNCT__ 3953 #define __FUNCT__ "SNESGetGS" 3954 PetscErrorCode SNESGetGS (SNES snes, PetscErrorCode(**func)(SNES, Vec, Vec, void*), void ** ctx) 3955 { 3956 PetscErrorCode ierr; 3957 DM dm; 3958 3959 PetscFunctionBegin; 3960 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3961 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 3962 ierr = DMSNESGetGS(dm,func,ctx);CHKERRQ(ierr); 3963 PetscFunctionReturn(0); 3964 } 3965 3966 #undef __FUNCT__ 3967 #define __FUNCT__ "SNESSetOptionsPrefix" 3968 /*@C 3969 SNESSetOptionsPrefix - Sets the prefix used for searching for all 3970 SNES options in the database. 3971 3972 Logically Collective on SNES 3973 3974 Input Parameter: 3975 + snes - the SNES context 3976 - prefix - the prefix to prepend to all option names 3977 3978 Notes: 3979 A hyphen (-) must NOT be given at the beginning of the prefix name. 3980 The first character of all runtime options is AUTOMATICALLY the hyphen. 3981 3982 Level: advanced 3983 3984 .keywords: SNES, set, options, prefix, database 3985 3986 .seealso: SNESSetFromOptions() 3987 @*/ 3988 PetscErrorCode SNESSetOptionsPrefix(SNES snes,const char prefix[]) 3989 { 3990 PetscErrorCode ierr; 3991 3992 PetscFunctionBegin; 3993 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3994 ierr = PetscObjectSetOptionsPrefix((PetscObject)snes,prefix);CHKERRQ(ierr); 3995 if (!snes->ksp) {ierr = SNESGetKSP(snes,&snes->ksp);CHKERRQ(ierr);} 3996 if (snes->linesearch) { 3997 ierr = SNESGetSNESLineSearch(snes,&snes->linesearch);CHKERRQ(ierr); 3998 ierr = PetscObjectSetOptionsPrefix((PetscObject)snes->linesearch,prefix);CHKERRQ(ierr); 3999 } 4000 ierr = KSPSetOptionsPrefix(snes->ksp,prefix);CHKERRQ(ierr); 4001 PetscFunctionReturn(0); 4002 } 4003 4004 #undef __FUNCT__ 4005 #define __FUNCT__ "SNESAppendOptionsPrefix" 4006 /*@C 4007 SNESAppendOptionsPrefix - Appends to the prefix used for searching for all 4008 SNES options in the database. 4009 4010 Logically Collective on SNES 4011 4012 Input Parameters: 4013 + snes - the SNES context 4014 - prefix - the prefix to prepend to all option names 4015 4016 Notes: 4017 A hyphen (-) must NOT be given at the beginning of the prefix name. 4018 The first character of all runtime options is AUTOMATICALLY the hyphen. 4019 4020 Level: advanced 4021 4022 .keywords: SNES, append, options, prefix, database 4023 4024 .seealso: SNESGetOptionsPrefix() 4025 @*/ 4026 PetscErrorCode SNESAppendOptionsPrefix(SNES snes,const char prefix[]) 4027 { 4028 PetscErrorCode ierr; 4029 4030 PetscFunctionBegin; 4031 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4032 ierr = PetscObjectAppendOptionsPrefix((PetscObject)snes,prefix);CHKERRQ(ierr); 4033 if (!snes->ksp) {ierr = SNESGetKSP(snes,&snes->ksp);CHKERRQ(ierr);} 4034 if (snes->linesearch) { 4035 ierr = SNESGetSNESLineSearch(snes,&snes->linesearch);CHKERRQ(ierr); 4036 ierr = PetscObjectAppendOptionsPrefix((PetscObject)snes->linesearch,prefix);CHKERRQ(ierr); 4037 } 4038 ierr = KSPAppendOptionsPrefix(snes->ksp,prefix);CHKERRQ(ierr); 4039 PetscFunctionReturn(0); 4040 } 4041 4042 #undef __FUNCT__ 4043 #define __FUNCT__ "SNESGetOptionsPrefix" 4044 /*@C 4045 SNESGetOptionsPrefix - Sets the prefix used for searching for all 4046 SNES options in the database. 4047 4048 Not Collective 4049 4050 Input Parameter: 4051 . snes - the SNES context 4052 4053 Output Parameter: 4054 . prefix - pointer to the prefix string used 4055 4056 Notes: On the fortran side, the user should pass in a string 'prefix' of 4057 sufficient length to hold the prefix. 4058 4059 Level: advanced 4060 4061 .keywords: SNES, get, options, prefix, database 4062 4063 .seealso: SNESAppendOptionsPrefix() 4064 @*/ 4065 PetscErrorCode SNESGetOptionsPrefix(SNES snes,const char *prefix[]) 4066 { 4067 PetscErrorCode ierr; 4068 4069 PetscFunctionBegin; 4070 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4071 ierr = PetscObjectGetOptionsPrefix((PetscObject)snes,prefix);CHKERRQ(ierr); 4072 PetscFunctionReturn(0); 4073 } 4074 4075 4076 #undef __FUNCT__ 4077 #define __FUNCT__ "SNESRegister" 4078 /*@C 4079 SNESRegister - See SNESRegisterDynamic() 4080 4081 Level: advanced 4082 @*/ 4083 PetscErrorCode SNESRegister(const char sname[],const char path[],const char name[],PetscErrorCode (*function)(SNES)) 4084 { 4085 char fullname[PETSC_MAX_PATH_LEN]; 4086 PetscErrorCode ierr; 4087 4088 PetscFunctionBegin; 4089 ierr = PetscFListConcat(path,name,fullname);CHKERRQ(ierr); 4090 ierr = PetscFListAdd(&SNESList,sname,fullname,(void (*)(void))function);CHKERRQ(ierr); 4091 PetscFunctionReturn(0); 4092 } 4093 4094 #undef __FUNCT__ 4095 #define __FUNCT__ "SNESTestLocalMin" 4096 PetscErrorCode SNESTestLocalMin(SNES snes) 4097 { 4098 PetscErrorCode ierr; 4099 PetscInt N,i,j; 4100 Vec u,uh,fh; 4101 PetscScalar value; 4102 PetscReal norm; 4103 4104 PetscFunctionBegin; 4105 ierr = SNESGetSolution(snes,&u);CHKERRQ(ierr); 4106 ierr = VecDuplicate(u,&uh);CHKERRQ(ierr); 4107 ierr = VecDuplicate(u,&fh);CHKERRQ(ierr); 4108 4109 /* currently only works for sequential */ 4110 ierr = PetscPrintf(PETSC_COMM_WORLD,"Testing FormFunction() for local min\n"); 4111 ierr = VecGetSize(u,&N);CHKERRQ(ierr); 4112 for (i=0; i<N; i++) { 4113 ierr = VecCopy(u,uh);CHKERRQ(ierr); 4114 ierr = PetscPrintf(PETSC_COMM_WORLD,"i = %D\n",i);CHKERRQ(ierr); 4115 for (j=-10; j<11; j++) { 4116 value = PetscSign(j)*exp(PetscAbs(j)-10.0); 4117 ierr = VecSetValue(uh,i,value,ADD_VALUES);CHKERRQ(ierr); 4118 ierr = SNESComputeFunction(snes,uh,fh);CHKERRQ(ierr); 4119 ierr = VecNorm(fh,NORM_2,&norm);CHKERRQ(ierr); 4120 ierr = PetscPrintf(PETSC_COMM_WORLD," j norm %D %18.16e\n",j,norm);CHKERRQ(ierr); 4121 value = -value; 4122 ierr = VecSetValue(uh,i,value,ADD_VALUES);CHKERRQ(ierr); 4123 } 4124 } 4125 ierr = VecDestroy(&uh);CHKERRQ(ierr); 4126 ierr = VecDestroy(&fh);CHKERRQ(ierr); 4127 PetscFunctionReturn(0); 4128 } 4129 4130 #undef __FUNCT__ 4131 #define __FUNCT__ "SNESKSPSetUseEW" 4132 /*@ 4133 SNESKSPSetUseEW - Sets SNES use Eisenstat-Walker method for 4134 computing relative tolerance for linear solvers within an inexact 4135 Newton method. 4136 4137 Logically Collective on SNES 4138 4139 Input Parameters: 4140 + snes - SNES context 4141 - flag - PETSC_TRUE or PETSC_FALSE 4142 4143 Options Database: 4144 + -snes_ksp_ew - use Eisenstat-Walker method for determining linear system convergence 4145 . -snes_ksp_ew_version ver - version of Eisenstat-Walker method 4146 . -snes_ksp_ew_rtol0 <rtol0> - Sets rtol0 4147 . -snes_ksp_ew_rtolmax <rtolmax> - Sets rtolmax 4148 . -snes_ksp_ew_gamma <gamma> - Sets gamma 4149 . -snes_ksp_ew_alpha <alpha> - Sets alpha 4150 . -snes_ksp_ew_alpha2 <alpha2> - Sets alpha2 4151 - -snes_ksp_ew_threshold <threshold> - Sets threshold 4152 4153 Notes: 4154 Currently, the default is to use a constant relative tolerance for 4155 the inner linear solvers. Alternatively, one can use the 4156 Eisenstat-Walker method, where the relative convergence tolerance 4157 is reset at each Newton iteration according progress of the nonlinear 4158 solver. 4159 4160 Level: advanced 4161 4162 Reference: 4163 S. C. Eisenstat and H. F. Walker, "Choosing the forcing terms in an 4164 inexact Newton method", SISC 17 (1), pp.16-32, 1996. 4165 4166 .keywords: SNES, KSP, Eisenstat, Walker, convergence, test, inexact, Newton 4167 4168 .seealso: SNESKSPGetUseEW(), SNESKSPGetParametersEW(), SNESKSPSetParametersEW() 4169 @*/ 4170 PetscErrorCode SNESKSPSetUseEW(SNES snes,PetscBool flag) 4171 { 4172 PetscFunctionBegin; 4173 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4174 PetscValidLogicalCollectiveBool(snes,flag,2); 4175 snes->ksp_ewconv = flag; 4176 PetscFunctionReturn(0); 4177 } 4178 4179 #undef __FUNCT__ 4180 #define __FUNCT__ "SNESKSPGetUseEW" 4181 /*@ 4182 SNESKSPGetUseEW - Gets if SNES is using Eisenstat-Walker method 4183 for computing relative tolerance for linear solvers within an 4184 inexact Newton method. 4185 4186 Not Collective 4187 4188 Input Parameter: 4189 . snes - SNES context 4190 4191 Output Parameter: 4192 . flag - PETSC_TRUE or PETSC_FALSE 4193 4194 Notes: 4195 Currently, the default is to use a constant relative tolerance for 4196 the inner linear solvers. Alternatively, one can use the 4197 Eisenstat-Walker method, where the relative convergence tolerance 4198 is reset at each Newton iteration according progress of the nonlinear 4199 solver. 4200 4201 Level: advanced 4202 4203 Reference: 4204 S. C. Eisenstat and H. F. Walker, "Choosing the forcing terms in an 4205 inexact Newton method", SISC 17 (1), pp.16-32, 1996. 4206 4207 .keywords: SNES, KSP, Eisenstat, Walker, convergence, test, inexact, Newton 4208 4209 .seealso: SNESKSPSetUseEW(), SNESKSPGetParametersEW(), SNESKSPSetParametersEW() 4210 @*/ 4211 PetscErrorCode SNESKSPGetUseEW(SNES snes, PetscBool *flag) 4212 { 4213 PetscFunctionBegin; 4214 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4215 PetscValidPointer(flag,2); 4216 *flag = snes->ksp_ewconv; 4217 PetscFunctionReturn(0); 4218 } 4219 4220 #undef __FUNCT__ 4221 #define __FUNCT__ "SNESKSPSetParametersEW" 4222 /*@ 4223 SNESKSPSetParametersEW - Sets parameters for Eisenstat-Walker 4224 convergence criteria for the linear solvers within an inexact 4225 Newton method. 4226 4227 Logically Collective on SNES 4228 4229 Input Parameters: 4230 + snes - SNES context 4231 . version - version 1, 2 (default is 2) or 3 4232 . rtol_0 - initial relative tolerance (0 <= rtol_0 < 1) 4233 . rtol_max - maximum relative tolerance (0 <= rtol_max < 1) 4234 . gamma - multiplicative factor for version 2 rtol computation 4235 (0 <= gamma2 <= 1) 4236 . alpha - power for version 2 rtol computation (1 < alpha <= 2) 4237 . alpha2 - power for safeguard 4238 - threshold - threshold for imposing safeguard (0 < threshold < 1) 4239 4240 Note: 4241 Version 3 was contributed by Luis Chacon, June 2006. 4242 4243 Use PETSC_DEFAULT to retain the default for any of the parameters. 4244 4245 Level: advanced 4246 4247 Reference: 4248 S. C. Eisenstat and H. F. Walker, "Choosing the forcing terms in an 4249 inexact Newton method", Utah State University Math. Stat. Dept. Res. 4250 Report 6/94/75, June, 1994, to appear in SIAM J. Sci. Comput. 4251 4252 .keywords: SNES, KSP, Eisenstat, Walker, set, parameters 4253 4254 .seealso: SNESKSPSetUseEW(), SNESKSPGetUseEW(), SNESKSPGetParametersEW() 4255 @*/ 4256 PetscErrorCode SNESKSPSetParametersEW(SNES snes,PetscInt version,PetscReal rtol_0,PetscReal rtol_max, 4257 PetscReal gamma,PetscReal alpha,PetscReal alpha2,PetscReal threshold) 4258 { 4259 SNESKSPEW *kctx; 4260 PetscFunctionBegin; 4261 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4262 kctx = (SNESKSPEW*)snes->kspconvctx; 4263 if (!kctx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"No Eisenstat-Walker context existing"); 4264 PetscValidLogicalCollectiveInt(snes,version,2); 4265 PetscValidLogicalCollectiveReal(snes,rtol_0,3); 4266 PetscValidLogicalCollectiveReal(snes,rtol_max,4); 4267 PetscValidLogicalCollectiveReal(snes,gamma,5); 4268 PetscValidLogicalCollectiveReal(snes,alpha,6); 4269 PetscValidLogicalCollectiveReal(snes,alpha2,7); 4270 PetscValidLogicalCollectiveReal(snes,threshold,8); 4271 4272 if (version != PETSC_DEFAULT) kctx->version = version; 4273 if (rtol_0 != PETSC_DEFAULT) kctx->rtol_0 = rtol_0; 4274 if (rtol_max != PETSC_DEFAULT) kctx->rtol_max = rtol_max; 4275 if (gamma != PETSC_DEFAULT) kctx->gamma = gamma; 4276 if (alpha != PETSC_DEFAULT) kctx->alpha = alpha; 4277 if (alpha2 != PETSC_DEFAULT) kctx->alpha2 = alpha2; 4278 if (threshold != PETSC_DEFAULT) kctx->threshold = threshold; 4279 4280 if (kctx->version < 1 || kctx->version > 3) { 4281 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Only versions 1, 2 and 3 are supported: %D",kctx->version); 4282 } 4283 if (kctx->rtol_0 < 0.0 || kctx->rtol_0 >= 1.0) { 4284 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"0.0 <= rtol_0 < 1.0: %G",kctx->rtol_0); 4285 } 4286 if (kctx->rtol_max < 0.0 || kctx->rtol_max >= 1.0) { 4287 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"0.0 <= rtol_max (%G) < 1.0\n",kctx->rtol_max); 4288 } 4289 if (kctx->gamma < 0.0 || kctx->gamma > 1.0) { 4290 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"0.0 <= gamma (%G) <= 1.0\n",kctx->gamma); 4291 } 4292 if (kctx->alpha <= 1.0 || kctx->alpha > 2.0) { 4293 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"1.0 < alpha (%G) <= 2.0\n",kctx->alpha); 4294 } 4295 if (kctx->threshold <= 0.0 || kctx->threshold >= 1.0) { 4296 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"0.0 < threshold (%G) < 1.0\n",kctx->threshold); 4297 } 4298 PetscFunctionReturn(0); 4299 } 4300 4301 #undef __FUNCT__ 4302 #define __FUNCT__ "SNESKSPGetParametersEW" 4303 /*@ 4304 SNESKSPGetParametersEW - Gets parameters for Eisenstat-Walker 4305 convergence criteria for the linear solvers within an inexact 4306 Newton method. 4307 4308 Not Collective 4309 4310 Input Parameters: 4311 snes - SNES context 4312 4313 Output Parameters: 4314 + version - version 1, 2 (default is 2) or 3 4315 . rtol_0 - initial relative tolerance (0 <= rtol_0 < 1) 4316 . rtol_max - maximum relative tolerance (0 <= rtol_max < 1) 4317 . gamma - multiplicative factor for version 2 rtol computation 4318 (0 <= gamma2 <= 1) 4319 . alpha - power for version 2 rtol computation (1 < alpha <= 2) 4320 . alpha2 - power for safeguard 4321 - threshold - threshold for imposing safeguard (0 < threshold < 1) 4322 4323 Level: advanced 4324 4325 .keywords: SNES, KSP, Eisenstat, Walker, get, parameters 4326 4327 .seealso: SNESKSPSetUseEW(), SNESKSPGetUseEW(), SNESKSPSetParametersEW() 4328 @*/ 4329 PetscErrorCode SNESKSPGetParametersEW(SNES snes,PetscInt *version,PetscReal *rtol_0,PetscReal *rtol_max, 4330 PetscReal *gamma,PetscReal *alpha,PetscReal *alpha2,PetscReal *threshold) 4331 { 4332 SNESKSPEW *kctx; 4333 PetscFunctionBegin; 4334 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4335 kctx = (SNESKSPEW*)snes->kspconvctx; 4336 if (!kctx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"No Eisenstat-Walker context existing"); 4337 if (version) *version = kctx->version; 4338 if (rtol_0) *rtol_0 = kctx->rtol_0; 4339 if (rtol_max) *rtol_max = kctx->rtol_max; 4340 if (gamma) *gamma = kctx->gamma; 4341 if (alpha) *alpha = kctx->alpha; 4342 if (alpha2) *alpha2 = kctx->alpha2; 4343 if (threshold) *threshold = kctx->threshold; 4344 PetscFunctionReturn(0); 4345 } 4346 4347 #undef __FUNCT__ 4348 #define __FUNCT__ "SNESKSPEW_PreSolve" 4349 static PetscErrorCode SNESKSPEW_PreSolve(SNES snes, KSP ksp, Vec b, Vec x) 4350 { 4351 PetscErrorCode ierr; 4352 SNESKSPEW *kctx = (SNESKSPEW*)snes->kspconvctx; 4353 PetscReal rtol=PETSC_DEFAULT,stol; 4354 4355 PetscFunctionBegin; 4356 if (!kctx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"No Eisenstat-Walker context exists"); 4357 if (!snes->iter) { /* first time in, so use the original user rtol */ 4358 rtol = kctx->rtol_0; 4359 } else { 4360 if (kctx->version == 1) { 4361 rtol = (snes->norm - kctx->lresid_last)/kctx->norm_last; 4362 if (rtol < 0.0) rtol = -rtol; 4363 stol = pow(kctx->rtol_last,kctx->alpha2); 4364 if (stol > kctx->threshold) rtol = PetscMax(rtol,stol); 4365 } else if (kctx->version == 2) { 4366 rtol = kctx->gamma * pow(snes->norm/kctx->norm_last,kctx->alpha); 4367 stol = kctx->gamma * pow(kctx->rtol_last,kctx->alpha); 4368 if (stol > kctx->threshold) rtol = PetscMax(rtol,stol); 4369 } else if (kctx->version == 3) {/* contributed by Luis Chacon, June 2006. */ 4370 rtol = kctx->gamma * pow(snes->norm/kctx->norm_last,kctx->alpha); 4371 /* safeguard: avoid sharp decrease of rtol */ 4372 stol = kctx->gamma*pow(kctx->rtol_last,kctx->alpha); 4373 stol = PetscMax(rtol,stol); 4374 rtol = PetscMin(kctx->rtol_0,stol); 4375 /* safeguard: avoid oversolving */ 4376 stol = kctx->gamma*(snes->ttol)/snes->norm; 4377 stol = PetscMax(rtol,stol); 4378 rtol = PetscMin(kctx->rtol_0,stol); 4379 } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Only versions 1, 2 or 3 are supported: %D",kctx->version); 4380 } 4381 /* safeguard: avoid rtol greater than one */ 4382 rtol = PetscMin(rtol,kctx->rtol_max); 4383 ierr = KSPSetTolerances(ksp,rtol,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);CHKERRQ(ierr); 4384 ierr = PetscInfo3(snes,"iter %D, Eisenstat-Walker (version %D) KSP rtol=%G\n",snes->iter,kctx->version,rtol);CHKERRQ(ierr); 4385 PetscFunctionReturn(0); 4386 } 4387 4388 #undef __FUNCT__ 4389 #define __FUNCT__ "SNESKSPEW_PostSolve" 4390 static PetscErrorCode SNESKSPEW_PostSolve(SNES snes, KSP ksp, Vec b, Vec x) 4391 { 4392 PetscErrorCode ierr; 4393 SNESKSPEW *kctx = (SNESKSPEW*)snes->kspconvctx; 4394 PCSide pcside; 4395 Vec lres; 4396 4397 PetscFunctionBegin; 4398 if (!kctx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"No Eisenstat-Walker context exists"); 4399 ierr = KSPGetTolerances(ksp,&kctx->rtol_last,0,0,0);CHKERRQ(ierr); 4400 ierr = SNESGetFunctionNorm(snes,&kctx->norm_last);CHKERRQ(ierr); 4401 if (kctx->version == 1) { 4402 ierr = KSPGetPCSide(ksp,&pcside);CHKERRQ(ierr); 4403 if (pcside == PC_RIGHT) { /* XXX Should we also test KSP_UNPRECONDITIONED_NORM ? */ 4404 /* KSP residual is true linear residual */ 4405 ierr = KSPGetResidualNorm(ksp,&kctx->lresid_last);CHKERRQ(ierr); 4406 } else { 4407 /* KSP residual is preconditioned residual */ 4408 /* compute true linear residual norm */ 4409 ierr = VecDuplicate(b,&lres);CHKERRQ(ierr); 4410 ierr = MatMult(snes->jacobian,x,lres);CHKERRQ(ierr); 4411 ierr = VecAYPX(lres,-1.0,b);CHKERRQ(ierr); 4412 ierr = VecNorm(lres,NORM_2,&kctx->lresid_last);CHKERRQ(ierr); 4413 ierr = VecDestroy(&lres);CHKERRQ(ierr); 4414 } 4415 } 4416 PetscFunctionReturn(0); 4417 } 4418 4419 #undef __FUNCT__ 4420 #define __FUNCT__ "SNES_KSPSolve" 4421 PetscErrorCode SNES_KSPSolve(SNES snes, KSP ksp, Vec b, Vec x) 4422 { 4423 PetscErrorCode ierr; 4424 4425 PetscFunctionBegin; 4426 if (snes->ksp_ewconv) { ierr = SNESKSPEW_PreSolve(snes,ksp,b,x);CHKERRQ(ierr); } 4427 ierr = KSPSolve(ksp,b,x);CHKERRQ(ierr); 4428 if (snes->ksp_ewconv) { ierr = SNESKSPEW_PostSolve(snes,ksp,b,x);CHKERRQ(ierr); } 4429 PetscFunctionReturn(0); 4430 } 4431 4432 #undef __FUNCT__ 4433 #define __FUNCT__ "SNESSetDM" 4434 /*@ 4435 SNESSetDM - Sets the DM that may be used by some preconditioners 4436 4437 Logically Collective on SNES 4438 4439 Input Parameters: 4440 + snes - the preconditioner context 4441 - dm - the dm 4442 4443 Level: intermediate 4444 4445 4446 .seealso: SNESGetDM(), KSPSetDM(), KSPGetDM() 4447 @*/ 4448 PetscErrorCode SNESSetDM(SNES snes,DM dm) 4449 { 4450 PetscErrorCode ierr; 4451 KSP ksp; 4452 DMSNES sdm; 4453 4454 PetscFunctionBegin; 4455 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4456 if (dm) {ierr = PetscObjectReference((PetscObject)dm);CHKERRQ(ierr);} 4457 if (snes->dm) { /* Move the DMSNES context over to the new DM unless the new DM already has one */ 4458 PetscContainer oldcontainer,container; 4459 ierr = PetscObjectQuery((PetscObject)snes->dm,"DMSNES",(PetscObject*)&oldcontainer);CHKERRQ(ierr); 4460 ierr = PetscObjectQuery((PetscObject)dm,"DMSNES",(PetscObject*)&container);CHKERRQ(ierr); 4461 if (oldcontainer && snes->dmAuto && !container) { 4462 ierr = DMCopyDMSNES(snes->dm,dm);CHKERRQ(ierr); 4463 ierr = DMGetDMSNES(snes->dm,&sdm);CHKERRQ(ierr); 4464 if (sdm->originaldm == snes->dm) { /* Grant write privileges to the replacement DM */ 4465 sdm->originaldm = dm; 4466 } 4467 } 4468 ierr = DMDestroy(&snes->dm);CHKERRQ(ierr); 4469 } 4470 snes->dm = dm; 4471 snes->dmAuto = PETSC_FALSE; 4472 ierr = SNESGetKSP(snes,&ksp);CHKERRQ(ierr); 4473 ierr = KSPSetDM(ksp,dm);CHKERRQ(ierr); 4474 ierr = KSPSetDMActive(ksp,PETSC_FALSE);CHKERRQ(ierr); 4475 if (snes->pc) { 4476 ierr = SNESSetDM(snes->pc, snes->dm);CHKERRQ(ierr); 4477 ierr = SNESSetPCSide(snes,snes->pcside);CHKERRQ(ierr); 4478 } 4479 PetscFunctionReturn(0); 4480 } 4481 4482 #undef __FUNCT__ 4483 #define __FUNCT__ "SNESGetDM" 4484 /*@ 4485 SNESGetDM - Gets the DM that may be used by some preconditioners 4486 4487 Not Collective but DM obtained is parallel on SNES 4488 4489 Input Parameter: 4490 . snes - the preconditioner context 4491 4492 Output Parameter: 4493 . dm - the dm 4494 4495 Level: intermediate 4496 4497 4498 .seealso: SNESSetDM(), KSPSetDM(), KSPGetDM() 4499 @*/ 4500 PetscErrorCode SNESGetDM(SNES snes,DM *dm) 4501 { 4502 PetscErrorCode ierr; 4503 4504 PetscFunctionBegin; 4505 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4506 if (!snes->dm) { 4507 ierr = DMShellCreate(((PetscObject)snes)->comm,&snes->dm);CHKERRQ(ierr); 4508 snes->dmAuto = PETSC_TRUE; 4509 } 4510 *dm = snes->dm; 4511 PetscFunctionReturn(0); 4512 } 4513 4514 #undef __FUNCT__ 4515 #define __FUNCT__ "SNESSetPC" 4516 /*@ 4517 SNESSetPC - Sets the nonlinear preconditioner to be used. 4518 4519 Collective on SNES 4520 4521 Input Parameters: 4522 + snes - iterative context obtained from SNESCreate() 4523 - pc - the preconditioner object 4524 4525 Notes: 4526 Use SNESGetPC() to retrieve the preconditioner context (for example, 4527 to configure it using the API). 4528 4529 Level: developer 4530 4531 .keywords: SNES, set, precondition 4532 .seealso: SNESGetPC() 4533 @*/ 4534 PetscErrorCode SNESSetPC(SNES snes, SNES pc) 4535 { 4536 PetscErrorCode ierr; 4537 4538 PetscFunctionBegin; 4539 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); 4540 PetscValidHeaderSpecific(pc, SNES_CLASSID, 2); 4541 PetscCheckSameComm(snes, 1, pc, 2); 4542 ierr = PetscObjectReference((PetscObject) pc);CHKERRQ(ierr); 4543 ierr = SNESDestroy(&snes->pc);CHKERRQ(ierr); 4544 snes->pc = pc; 4545 ierr = PetscLogObjectParent(snes, snes->pc);CHKERRQ(ierr); 4546 PetscFunctionReturn(0); 4547 } 4548 4549 #undef __FUNCT__ 4550 #define __FUNCT__ "SNESGetPC" 4551 /*@ 4552 SNESGetPC - Returns a pointer to the nonlinear preconditioning context set with SNESSetPC(). 4553 4554 Not Collective 4555 4556 Input Parameter: 4557 . snes - iterative context obtained from SNESCreate() 4558 4559 Output Parameter: 4560 . pc - preconditioner context 4561 4562 Level: developer 4563 4564 .keywords: SNES, get, preconditioner 4565 .seealso: SNESSetPC() 4566 @*/ 4567 PetscErrorCode SNESGetPC(SNES snes, SNES *pc) 4568 { 4569 PetscErrorCode ierr; 4570 const char *optionsprefix; 4571 4572 PetscFunctionBegin; 4573 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); 4574 PetscValidPointer(pc, 2); 4575 if (!snes->pc) { 4576 ierr = SNESCreate(((PetscObject) snes)->comm,&snes->pc);CHKERRQ(ierr); 4577 ierr = PetscObjectIncrementTabLevel((PetscObject)snes->pc,(PetscObject)snes,1);CHKERRQ(ierr); 4578 ierr = PetscLogObjectParent(snes,snes->pc);CHKERRQ(ierr); 4579 ierr = SNESGetOptionsPrefix(snes,&optionsprefix);CHKERRQ(ierr); 4580 ierr = SNESSetOptionsPrefix(snes->pc,optionsprefix);CHKERRQ(ierr); 4581 ierr = SNESAppendOptionsPrefix(snes->pc,"npc_");CHKERRQ(ierr); 4582 } 4583 *pc = snes->pc; 4584 PetscFunctionReturn(0); 4585 } 4586 4587 4588 #undef __FUNCT__ 4589 #define __FUNCT__ "SNESSetPCSide" 4590 /*@ 4591 SNESSetPCSide - Sets the preconditioning side. 4592 4593 Logically Collective on SNES 4594 4595 Input Parameter: 4596 . snes - iterative context obtained from SNESCreate() 4597 4598 Output Parameter: 4599 . side - the preconditioning side, where side is one of 4600 .vb 4601 PC_LEFT - left preconditioning (default) 4602 PC_RIGHT - right preconditioning 4603 .ve 4604 4605 Options Database Keys: 4606 . -snes_pc_side <right,left> 4607 4608 Level: intermediate 4609 4610 .keywords: SNES, set, right, left, side, preconditioner, flag 4611 4612 .seealso: SNESGetPCSide(), KSPSetPCSide() 4613 @*/ 4614 PetscErrorCode SNESSetPCSide(SNES snes,PCSide side) 4615 { 4616 PetscFunctionBegin; 4617 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4618 PetscValidLogicalCollectiveEnum(snes,side,2); 4619 snes->pcside = side; 4620 PetscFunctionReturn(0); 4621 } 4622 4623 #undef __FUNCT__ 4624 #define __FUNCT__ "SNESGetPCSide" 4625 /*@ 4626 SNESGetPCSide - Gets the preconditioning side. 4627 4628 Not Collective 4629 4630 Input Parameter: 4631 . snes - iterative context obtained from SNESCreate() 4632 4633 Output Parameter: 4634 . side - the preconditioning side, where side is one of 4635 .vb 4636 PC_LEFT - left preconditioning (default) 4637 PC_RIGHT - right preconditioning 4638 .ve 4639 4640 Level: intermediate 4641 4642 .keywords: SNES, get, right, left, side, preconditioner, flag 4643 4644 .seealso: SNESSetPCSide(), KSPGetPCSide() 4645 @*/ 4646 PetscErrorCode SNESGetPCSide(SNES snes,PCSide *side) 4647 { 4648 PetscFunctionBegin; 4649 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4650 PetscValidPointer(side,2); 4651 *side = snes->pcside; 4652 PetscFunctionReturn(0); 4653 } 4654 4655 #undef __FUNCT__ 4656 #define __FUNCT__ "SNESSetSNESLineSearch" 4657 /*@ 4658 SNESSetSNESLineSearch - Sets the linesearch on the SNES instance. 4659 4660 Collective on SNES 4661 4662 Input Parameters: 4663 + snes - iterative context obtained from SNESCreate() 4664 - linesearch - the linesearch object 4665 4666 Notes: 4667 Use SNESGetSNESLineSearch() to retrieve the preconditioner context (for example, 4668 to configure it using the API). 4669 4670 Level: developer 4671 4672 .keywords: SNES, set, linesearch 4673 .seealso: SNESGetSNESLineSearch() 4674 @*/ 4675 PetscErrorCode SNESSetSNESLineSearch(SNES snes, SNESLineSearch linesearch) 4676 { 4677 PetscErrorCode ierr; 4678 4679 PetscFunctionBegin; 4680 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); 4681 PetscValidHeaderSpecific(linesearch, SNESLINESEARCH_CLASSID, 2); 4682 PetscCheckSameComm(snes, 1, linesearch, 2); 4683 ierr = PetscObjectReference((PetscObject) linesearch);CHKERRQ(ierr); 4684 ierr = SNESLineSearchDestroy(&snes->linesearch);CHKERRQ(ierr); 4685 snes->linesearch = linesearch; 4686 ierr = PetscLogObjectParent(snes, snes->linesearch);CHKERRQ(ierr); 4687 PetscFunctionReturn(0); 4688 } 4689 4690 #undef __FUNCT__ 4691 #define __FUNCT__ "SNESGetSNESLineSearch" 4692 /*@C 4693 SNESGetSNESLineSearch - Returns a pointer to the line search context set with SNESSetLineSearch() 4694 or creates a default line search instance associated with the SNES and returns it. 4695 4696 Not Collective 4697 4698 Input Parameter: 4699 . snes - iterative context obtained from SNESCreate() 4700 4701 Output Parameter: 4702 . linesearch - linesearch context 4703 4704 Level: developer 4705 4706 .keywords: SNES, get, linesearch 4707 .seealso: SNESSetSNESLineSearch() 4708 @*/ 4709 PetscErrorCode SNESGetSNESLineSearch(SNES snes, SNESLineSearch *linesearch) 4710 { 4711 PetscErrorCode ierr; 4712 const char *optionsprefix; 4713 4714 PetscFunctionBegin; 4715 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); 4716 PetscValidPointer(linesearch, 2); 4717 if (!snes->linesearch) { 4718 ierr = SNESGetOptionsPrefix(snes, &optionsprefix);CHKERRQ(ierr); 4719 ierr = SNESLineSearchCreate(((PetscObject) snes)->comm, &snes->linesearch);CHKERRQ(ierr); 4720 ierr = SNESLineSearchSetSNES(snes->linesearch, snes);CHKERRQ(ierr); 4721 ierr = SNESLineSearchAppendOptionsPrefix(snes->linesearch, optionsprefix);CHKERRQ(ierr); 4722 ierr = PetscObjectIncrementTabLevel((PetscObject) snes->linesearch, (PetscObject) snes, 1);CHKERRQ(ierr); 4723 ierr = PetscLogObjectParent(snes, snes->linesearch);CHKERRQ(ierr); 4724 } 4725 *linesearch = snes->linesearch; 4726 PetscFunctionReturn(0); 4727 } 4728 4729 #if defined(PETSC_HAVE_MATLAB_ENGINE) 4730 #include <mex.h> 4731 4732 typedef struct {char *funcname; mxArray *ctx;} SNESMatlabContext; 4733 4734 #undef __FUNCT__ 4735 #define __FUNCT__ "SNESComputeFunction_Matlab" 4736 /* 4737 SNESComputeFunction_Matlab - Calls the function that has been set with 4738 SNESSetFunctionMatlab(). 4739 4740 Collective on SNES 4741 4742 Input Parameters: 4743 + snes - the SNES context 4744 - x - input vector 4745 4746 Output Parameter: 4747 . y - function vector, as set by SNESSetFunction() 4748 4749 Notes: 4750 SNESComputeFunction() is typically used within nonlinear solvers 4751 implementations, so most users would not generally call this routine 4752 themselves. 4753 4754 Level: developer 4755 4756 .keywords: SNES, nonlinear, compute, function 4757 4758 .seealso: SNESSetFunction(), SNESGetFunction() 4759 */ 4760 PetscErrorCode SNESComputeFunction_Matlab(SNES snes,Vec x,Vec y, void *ctx) 4761 { 4762 PetscErrorCode ierr; 4763 SNESMatlabContext *sctx = (SNESMatlabContext *)ctx; 4764 int nlhs = 1,nrhs = 5; 4765 mxArray *plhs[1],*prhs[5]; 4766 long long int lx = 0,ly = 0,ls = 0; 4767 4768 PetscFunctionBegin; 4769 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4770 PetscValidHeaderSpecific(x,VEC_CLASSID,2); 4771 PetscValidHeaderSpecific(y,VEC_CLASSID,3); 4772 PetscCheckSameComm(snes,1,x,2); 4773 PetscCheckSameComm(snes,1,y,3); 4774 4775 /* call Matlab function in ctx with arguments x and y */ 4776 4777 ierr = PetscMemcpy(&ls,&snes,sizeof(snes));CHKERRQ(ierr); 4778 ierr = PetscMemcpy(&lx,&x,sizeof(x));CHKERRQ(ierr); 4779 ierr = PetscMemcpy(&ly,&y,sizeof(x));CHKERRQ(ierr); 4780 prhs[0] = mxCreateDoubleScalar((double)ls); 4781 prhs[1] = mxCreateDoubleScalar((double)lx); 4782 prhs[2] = mxCreateDoubleScalar((double)ly); 4783 prhs[3] = mxCreateString(sctx->funcname); 4784 prhs[4] = sctx->ctx; 4785 ierr = mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscSNESComputeFunctionInternal");CHKERRQ(ierr); 4786 ierr = mxGetScalar(plhs[0]);CHKERRQ(ierr); 4787 mxDestroyArray(prhs[0]); 4788 mxDestroyArray(prhs[1]); 4789 mxDestroyArray(prhs[2]); 4790 mxDestroyArray(prhs[3]); 4791 mxDestroyArray(plhs[0]); 4792 PetscFunctionReturn(0); 4793 } 4794 4795 4796 #undef __FUNCT__ 4797 #define __FUNCT__ "SNESSetFunctionMatlab" 4798 /* 4799 SNESSetFunctionMatlab - Sets the function evaluation routine and function 4800 vector for use by the SNES routines in solving systems of nonlinear 4801 equations from MATLAB. Here the function is a string containing the name of a MATLAB function 4802 4803 Logically Collective on SNES 4804 4805 Input Parameters: 4806 + snes - the SNES context 4807 . r - vector to store function value 4808 - func - function evaluation routine 4809 4810 Calling sequence of func: 4811 $ func (SNES snes,Vec x,Vec f,void *ctx); 4812 4813 4814 Notes: 4815 The Newton-like methods typically solve linear systems of the form 4816 $ f'(x) x = -f(x), 4817 where f'(x) denotes the Jacobian matrix and f(x) is the function. 4818 4819 Level: beginner 4820 4821 .keywords: SNES, nonlinear, set, function 4822 4823 .seealso: SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESSetFunction() 4824 */ 4825 PetscErrorCode SNESSetFunctionMatlab(SNES snes,Vec r,const char *func,mxArray *ctx) 4826 { 4827 PetscErrorCode ierr; 4828 SNESMatlabContext *sctx; 4829 4830 PetscFunctionBegin; 4831 /* currently sctx is memory bleed */ 4832 ierr = PetscMalloc(sizeof(SNESMatlabContext),&sctx);CHKERRQ(ierr); 4833 ierr = PetscStrallocpy(func,&sctx->funcname);CHKERRQ(ierr); 4834 /* 4835 This should work, but it doesn't 4836 sctx->ctx = ctx; 4837 mexMakeArrayPersistent(sctx->ctx); 4838 */ 4839 sctx->ctx = mxDuplicateArray(ctx); 4840 ierr = SNESSetFunction(snes,r,SNESComputeFunction_Matlab,sctx);CHKERRQ(ierr); 4841 PetscFunctionReturn(0); 4842 } 4843 4844 #undef __FUNCT__ 4845 #define __FUNCT__ "SNESComputeJacobian_Matlab" 4846 /* 4847 SNESComputeJacobian_Matlab - Calls the function that has been set with 4848 SNESSetJacobianMatlab(). 4849 4850 Collective on SNES 4851 4852 Input Parameters: 4853 + snes - the SNES context 4854 . x - input vector 4855 . A, B - the matrices 4856 - ctx - user context 4857 4858 Output Parameter: 4859 . flag - structure of the matrix 4860 4861 Level: developer 4862 4863 .keywords: SNES, nonlinear, compute, function 4864 4865 .seealso: SNESSetFunction(), SNESGetFunction() 4866 @*/ 4867 PetscErrorCode SNESComputeJacobian_Matlab(SNES snes,Vec x,Mat *A,Mat *B,MatStructure *flag, void *ctx) 4868 { 4869 PetscErrorCode ierr; 4870 SNESMatlabContext *sctx = (SNESMatlabContext *)ctx; 4871 int nlhs = 2,nrhs = 6; 4872 mxArray *plhs[2],*prhs[6]; 4873 long long int lx = 0,lA = 0,ls = 0, lB = 0; 4874 4875 PetscFunctionBegin; 4876 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4877 PetscValidHeaderSpecific(x,VEC_CLASSID,2); 4878 4879 /* call Matlab function in ctx with arguments x and y */ 4880 4881 ierr = PetscMemcpy(&ls,&snes,sizeof(snes));CHKERRQ(ierr); 4882 ierr = PetscMemcpy(&lx,&x,sizeof(x));CHKERRQ(ierr); 4883 ierr = PetscMemcpy(&lA,A,sizeof(x));CHKERRQ(ierr); 4884 ierr = PetscMemcpy(&lB,B,sizeof(x));CHKERRQ(ierr); 4885 prhs[0] = mxCreateDoubleScalar((double)ls); 4886 prhs[1] = mxCreateDoubleScalar((double)lx); 4887 prhs[2] = mxCreateDoubleScalar((double)lA); 4888 prhs[3] = mxCreateDoubleScalar((double)lB); 4889 prhs[4] = mxCreateString(sctx->funcname); 4890 prhs[5] = sctx->ctx; 4891 ierr = mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscSNESComputeJacobianInternal");CHKERRQ(ierr); 4892 ierr = mxGetScalar(plhs[0]);CHKERRQ(ierr); 4893 *flag = (MatStructure) mxGetScalar(plhs[1]);CHKERRQ(ierr); 4894 mxDestroyArray(prhs[0]); 4895 mxDestroyArray(prhs[1]); 4896 mxDestroyArray(prhs[2]); 4897 mxDestroyArray(prhs[3]); 4898 mxDestroyArray(prhs[4]); 4899 mxDestroyArray(plhs[0]); 4900 mxDestroyArray(plhs[1]); 4901 PetscFunctionReturn(0); 4902 } 4903 4904 4905 #undef __FUNCT__ 4906 #define __FUNCT__ "SNESSetJacobianMatlab" 4907 /* 4908 SNESSetJacobianMatlab - Sets the Jacobian function evaluation routine and two empty Jacobian matrices 4909 vector for use by the SNES routines in solving systems of nonlinear 4910 equations from MATLAB. Here the function is a string containing the name of a MATLAB function 4911 4912 Logically Collective on SNES 4913 4914 Input Parameters: 4915 + snes - the SNES context 4916 . A,B - Jacobian matrices 4917 . func - function evaluation routine 4918 - ctx - user context 4919 4920 Calling sequence of func: 4921 $ flag = func (SNES snes,Vec x,Mat A,Mat B,void *ctx); 4922 4923 4924 Level: developer 4925 4926 .keywords: SNES, nonlinear, set, function 4927 4928 .seealso: SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESSetFunction() 4929 */ 4930 PetscErrorCode SNESSetJacobianMatlab(SNES snes,Mat A,Mat B,const char *func,mxArray *ctx) 4931 { 4932 PetscErrorCode ierr; 4933 SNESMatlabContext *sctx; 4934 4935 PetscFunctionBegin; 4936 /* currently sctx is memory bleed */ 4937 ierr = PetscMalloc(sizeof(SNESMatlabContext),&sctx);CHKERRQ(ierr); 4938 ierr = PetscStrallocpy(func,&sctx->funcname);CHKERRQ(ierr); 4939 /* 4940 This should work, but it doesn't 4941 sctx->ctx = ctx; 4942 mexMakeArrayPersistent(sctx->ctx); 4943 */ 4944 sctx->ctx = mxDuplicateArray(ctx); 4945 ierr = SNESSetJacobian(snes,A,B,SNESComputeJacobian_Matlab,sctx);CHKERRQ(ierr); 4946 PetscFunctionReturn(0); 4947 } 4948 4949 #undef __FUNCT__ 4950 #define __FUNCT__ "SNESMonitor_Matlab" 4951 /* 4952 SNESMonitor_Matlab - Calls the function that has been set with SNESMonitorSetMatlab(). 4953 4954 Collective on SNES 4955 4956 .seealso: SNESSetFunction(), SNESGetFunction() 4957 @*/ 4958 PetscErrorCode SNESMonitor_Matlab(SNES snes,PetscInt it, PetscReal fnorm, void *ctx) 4959 { 4960 PetscErrorCode ierr; 4961 SNESMatlabContext *sctx = (SNESMatlabContext *)ctx; 4962 int nlhs = 1,nrhs = 6; 4963 mxArray *plhs[1],*prhs[6]; 4964 long long int lx = 0,ls = 0; 4965 Vec x=snes->vec_sol; 4966 4967 PetscFunctionBegin; 4968 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4969 4970 ierr = PetscMemcpy(&ls,&snes,sizeof(snes));CHKERRQ(ierr); 4971 ierr = PetscMemcpy(&lx,&x,sizeof(x));CHKERRQ(ierr); 4972 prhs[0] = mxCreateDoubleScalar((double)ls); 4973 prhs[1] = mxCreateDoubleScalar((double)it); 4974 prhs[2] = mxCreateDoubleScalar((double)fnorm); 4975 prhs[3] = mxCreateDoubleScalar((double)lx); 4976 prhs[4] = mxCreateString(sctx->funcname); 4977 prhs[5] = sctx->ctx; 4978 ierr = mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscSNESMonitorInternal");CHKERRQ(ierr); 4979 ierr = mxGetScalar(plhs[0]);CHKERRQ(ierr); 4980 mxDestroyArray(prhs[0]); 4981 mxDestroyArray(prhs[1]); 4982 mxDestroyArray(prhs[2]); 4983 mxDestroyArray(prhs[3]); 4984 mxDestroyArray(prhs[4]); 4985 mxDestroyArray(plhs[0]); 4986 PetscFunctionReturn(0); 4987 } 4988 4989 4990 #undef __FUNCT__ 4991 #define __FUNCT__ "SNESMonitorSetMatlab" 4992 /* 4993 SNESMonitorSetMatlab - Sets the monitor function from MATLAB 4994 4995 Level: developer 4996 4997 .keywords: SNES, nonlinear, set, function 4998 4999 .seealso: SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESSetFunction() 5000 */ 5001 PetscErrorCode SNESMonitorSetMatlab(SNES snes,const char *func,mxArray *ctx) 5002 { 5003 PetscErrorCode ierr; 5004 SNESMatlabContext *sctx; 5005 5006 PetscFunctionBegin; 5007 /* currently sctx is memory bleed */ 5008 ierr = PetscMalloc(sizeof(SNESMatlabContext),&sctx);CHKERRQ(ierr); 5009 ierr = PetscStrallocpy(func,&sctx->funcname);CHKERRQ(ierr); 5010 /* 5011 This should work, but it doesn't 5012 sctx->ctx = ctx; 5013 mexMakeArrayPersistent(sctx->ctx); 5014 */ 5015 sctx->ctx = mxDuplicateArray(ctx); 5016 ierr = SNESMonitorSet(snes,SNESMonitor_Matlab,sctx,PETSC_NULL);CHKERRQ(ierr); 5017 PetscFunctionReturn(0); 5018 } 5019 5020 #endif 5021