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 of the change in the solution between steps, || delta x || < stol*|| x || 2893 . maxit - maximum number of iterations 2894 - maxf - maximum number of function evaluations 2895 2896 Options Database Keys: 2897 + -snes_atol <abstol> - Sets abstol 2898 . -snes_rtol <rtol> - Sets rtol 2899 . -snes_stol <stol> - Sets stol 2900 . -snes_max_it <maxit> - Sets maxit 2901 - -snes_max_funcs <maxf> - Sets maxf 2902 2903 Notes: 2904 The default maximum number of iterations is 50. 2905 The default maximum number of function evaluations is 1000. 2906 2907 Level: intermediate 2908 2909 .keywords: SNES, nonlinear, set, convergence, tolerances 2910 2911 .seealso: SNESSetTrustRegionTolerance() 2912 @*/ 2913 PetscErrorCode SNESSetTolerances(SNES snes,PetscReal abstol,PetscReal rtol,PetscReal stol,PetscInt maxit,PetscInt maxf) 2914 { 2915 PetscFunctionBegin; 2916 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2917 PetscValidLogicalCollectiveReal(snes,abstol,2); 2918 PetscValidLogicalCollectiveReal(snes,rtol,3); 2919 PetscValidLogicalCollectiveReal(snes,stol,4); 2920 PetscValidLogicalCollectiveInt(snes,maxit,5); 2921 PetscValidLogicalCollectiveInt(snes,maxf,6); 2922 2923 if (abstol != PETSC_DEFAULT) { 2924 if (abstol < 0.0) SETERRQ1(((PetscObject)snes)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Absolute tolerance %G must be non-negative",abstol); 2925 snes->abstol = abstol; 2926 } 2927 if (rtol != PETSC_DEFAULT) { 2928 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); 2929 snes->rtol = rtol; 2930 } 2931 if (stol != PETSC_DEFAULT) { 2932 if (stol < 0.0) SETERRQ1(((PetscObject)snes)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Step tolerance %G must be non-negative",stol); 2933 snes->stol = stol; 2934 } 2935 if (maxit != PETSC_DEFAULT) { 2936 if (maxit < 0) SETERRQ1(((PetscObject)snes)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Maximum number of iterations %D must be non-negative",maxit); 2937 snes->max_its = maxit; 2938 } 2939 if (maxf != PETSC_DEFAULT) { 2940 if (maxf < 0) SETERRQ1(((PetscObject)snes)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Maximum number of function evaluations %D must be non-negative",maxf); 2941 snes->max_funcs = maxf; 2942 } 2943 snes->tolerancesset = PETSC_TRUE; 2944 PetscFunctionReturn(0); 2945 } 2946 2947 #undef __FUNCT__ 2948 #define __FUNCT__ "SNESGetTolerances" 2949 /*@ 2950 SNESGetTolerances - Gets various parameters used in convergence tests. 2951 2952 Not Collective 2953 2954 Input Parameters: 2955 + snes - the SNES context 2956 . atol - absolute convergence tolerance 2957 . rtol - relative convergence tolerance 2958 . stol - convergence tolerance in terms of the norm 2959 of the change in the solution between steps 2960 . maxit - maximum number of iterations 2961 - maxf - maximum number of function evaluations 2962 2963 Notes: 2964 The user can specify PETSC_NULL for any parameter that is not needed. 2965 2966 Level: intermediate 2967 2968 .keywords: SNES, nonlinear, get, convergence, tolerances 2969 2970 .seealso: SNESSetTolerances() 2971 @*/ 2972 PetscErrorCode SNESGetTolerances(SNES snes,PetscReal *atol,PetscReal *rtol,PetscReal *stol,PetscInt *maxit,PetscInt *maxf) 2973 { 2974 PetscFunctionBegin; 2975 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2976 if (atol) *atol = snes->abstol; 2977 if (rtol) *rtol = snes->rtol; 2978 if (stol) *stol = snes->stol; 2979 if (maxit) *maxit = snes->max_its; 2980 if (maxf) *maxf = snes->max_funcs; 2981 PetscFunctionReturn(0); 2982 } 2983 2984 #undef __FUNCT__ 2985 #define __FUNCT__ "SNESSetTrustRegionTolerance" 2986 /*@ 2987 SNESSetTrustRegionTolerance - Sets the trust region parameter tolerance. 2988 2989 Logically Collective on SNES 2990 2991 Input Parameters: 2992 + snes - the SNES context 2993 - tol - tolerance 2994 2995 Options Database Key: 2996 . -snes_trtol <tol> - Sets tol 2997 2998 Level: intermediate 2999 3000 .keywords: SNES, nonlinear, set, trust region, tolerance 3001 3002 .seealso: SNESSetTolerances() 3003 @*/ 3004 PetscErrorCode SNESSetTrustRegionTolerance(SNES snes,PetscReal tol) 3005 { 3006 PetscFunctionBegin; 3007 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3008 PetscValidLogicalCollectiveReal(snes,tol,2); 3009 snes->deltatol = tol; 3010 PetscFunctionReturn(0); 3011 } 3012 3013 /* 3014 Duplicate the lg monitors for SNES from KSP; for some reason with 3015 dynamic libraries things don't work under Sun4 if we just use 3016 macros instead of functions 3017 */ 3018 #undef __FUNCT__ 3019 #define __FUNCT__ "SNESMonitorLGResidualNorm" 3020 PetscErrorCode SNESMonitorLGResidualNorm(SNES snes,PetscInt it,PetscReal norm,void *ctx) 3021 { 3022 PetscErrorCode ierr; 3023 3024 PetscFunctionBegin; 3025 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3026 ierr = KSPMonitorLGResidualNorm((KSP)snes,it,norm,ctx);CHKERRQ(ierr); 3027 PetscFunctionReturn(0); 3028 } 3029 3030 #undef __FUNCT__ 3031 #define __FUNCT__ "SNESMonitorLGCreate" 3032 PetscErrorCode SNESMonitorLGCreate(const char host[],const char label[],int x,int y,int m,int n,PetscDrawLG *draw) 3033 { 3034 PetscErrorCode ierr; 3035 3036 PetscFunctionBegin; 3037 ierr = KSPMonitorLGResidualNormCreate(host,label,x,y,m,n,draw);CHKERRQ(ierr); 3038 PetscFunctionReturn(0); 3039 } 3040 3041 #undef __FUNCT__ 3042 #define __FUNCT__ "SNESMonitorLGDestroy" 3043 PetscErrorCode SNESMonitorLGDestroy(PetscDrawLG *draw) 3044 { 3045 PetscErrorCode ierr; 3046 3047 PetscFunctionBegin; 3048 ierr = KSPMonitorLGResidualNormDestroy(draw);CHKERRQ(ierr); 3049 PetscFunctionReturn(0); 3050 } 3051 3052 extern PetscErrorCode SNESMonitorRange_Private(SNES,PetscInt,PetscReal*); 3053 #undef __FUNCT__ 3054 #define __FUNCT__ "SNESMonitorLGRange" 3055 PetscErrorCode SNESMonitorLGRange(SNES snes,PetscInt n,PetscReal rnorm,void *monctx) 3056 { 3057 PetscDrawLG lg; 3058 PetscErrorCode ierr; 3059 PetscReal x,y,per; 3060 PetscViewer v = (PetscViewer)monctx; 3061 static PetscReal prev; /* should be in the context */ 3062 PetscDraw draw; 3063 3064 PetscFunctionBegin; 3065 ierr = PetscViewerDrawGetDrawLG(v,0,&lg);CHKERRQ(ierr); 3066 if (!n) {ierr = PetscDrawLGReset(lg);CHKERRQ(ierr);} 3067 ierr = PetscDrawLGGetDraw(lg,&draw);CHKERRQ(ierr); 3068 ierr = PetscDrawSetTitle(draw,"Residual norm");CHKERRQ(ierr); 3069 x = (PetscReal) n; 3070 if (rnorm > 0.0) y = log10(rnorm); else y = -15.0; 3071 ierr = PetscDrawLGAddPoint(lg,&x,&y);CHKERRQ(ierr); 3072 if (n < 20 || !(n % 5)) { 3073 ierr = PetscDrawLGDraw(lg);CHKERRQ(ierr); 3074 } 3075 3076 ierr = PetscViewerDrawGetDrawLG(v,1,&lg);CHKERRQ(ierr); 3077 if (!n) {ierr = PetscDrawLGReset(lg);CHKERRQ(ierr);} 3078 ierr = PetscDrawLGGetDraw(lg,&draw);CHKERRQ(ierr); 3079 ierr = PetscDrawSetTitle(draw,"% elemts > .2*max elemt");CHKERRQ(ierr); 3080 ierr = SNESMonitorRange_Private(snes,n,&per);CHKERRQ(ierr); 3081 x = (PetscReal) n; 3082 y = 100.0*per; 3083 ierr = PetscDrawLGAddPoint(lg,&x,&y);CHKERRQ(ierr); 3084 if (n < 20 || !(n % 5)) { 3085 ierr = PetscDrawLGDraw(lg);CHKERRQ(ierr); 3086 } 3087 3088 ierr = PetscViewerDrawGetDrawLG(v,2,&lg);CHKERRQ(ierr); 3089 if (!n) {prev = rnorm;ierr = PetscDrawLGReset(lg);CHKERRQ(ierr);} 3090 ierr = PetscDrawLGGetDraw(lg,&draw);CHKERRQ(ierr); 3091 ierr = PetscDrawSetTitle(draw,"(norm -oldnorm)/oldnorm");CHKERRQ(ierr); 3092 x = (PetscReal) n; 3093 y = (prev - rnorm)/prev; 3094 ierr = PetscDrawLGAddPoint(lg,&x,&y);CHKERRQ(ierr); 3095 if (n < 20 || !(n % 5)) { 3096 ierr = PetscDrawLGDraw(lg);CHKERRQ(ierr); 3097 } 3098 3099 ierr = PetscViewerDrawGetDrawLG(v,3,&lg);CHKERRQ(ierr); 3100 if (!n) {ierr = PetscDrawLGReset(lg);CHKERRQ(ierr);} 3101 ierr = PetscDrawLGGetDraw(lg,&draw);CHKERRQ(ierr); 3102 ierr = PetscDrawSetTitle(draw,"(norm -oldnorm)/oldnorm*(% > .2 max)");CHKERRQ(ierr); 3103 x = (PetscReal) n; 3104 y = (prev - rnorm)/(prev*per); 3105 if (n > 2) { /*skip initial crazy value */ 3106 ierr = PetscDrawLGAddPoint(lg,&x,&y);CHKERRQ(ierr); 3107 } 3108 if (n < 20 || !(n % 5)) { 3109 ierr = PetscDrawLGDraw(lg);CHKERRQ(ierr); 3110 } 3111 prev = rnorm; 3112 PetscFunctionReturn(0); 3113 } 3114 3115 #undef __FUNCT__ 3116 #define __FUNCT__ "SNESMonitor" 3117 /*@ 3118 SNESMonitor - runs the user provided monitor routines, if they exist 3119 3120 Collective on SNES 3121 3122 Input Parameters: 3123 + snes - nonlinear solver context obtained from SNESCreate() 3124 . iter - iteration number 3125 - rnorm - relative norm of the residual 3126 3127 Notes: 3128 This routine is called by the SNES implementations. 3129 It does not typically need to be called by the user. 3130 3131 Level: developer 3132 3133 .seealso: SNESMonitorSet() 3134 @*/ 3135 PetscErrorCode SNESMonitor(SNES snes,PetscInt iter,PetscReal rnorm) 3136 { 3137 PetscErrorCode ierr; 3138 PetscInt i,n = snes->numbermonitors; 3139 3140 PetscFunctionBegin; 3141 for (i=0; i<n; i++) { 3142 ierr = (*snes->monitor[i])(snes,iter,rnorm,snes->monitorcontext[i]);CHKERRQ(ierr); 3143 } 3144 PetscFunctionReturn(0); 3145 } 3146 3147 /* ------------ Routines to set performance monitoring options ----------- */ 3148 3149 #undef __FUNCT__ 3150 #define __FUNCT__ "SNESMonitorSet" 3151 /*@C 3152 SNESMonitorSet - Sets an ADDITIONAL function that is to be used at every 3153 iteration of the nonlinear solver to display the iteration's 3154 progress. 3155 3156 Logically Collective on SNES 3157 3158 Input Parameters: 3159 + snes - the SNES context 3160 . func - monitoring routine 3161 . mctx - [optional] user-defined context for private data for the 3162 monitor routine (use PETSC_NULL if no context is desired) 3163 - monitordestroy - [optional] routine that frees monitor context 3164 (may be PETSC_NULL) 3165 3166 Calling sequence of func: 3167 $ int func(SNES snes,PetscInt its, PetscReal norm,void *mctx) 3168 3169 + snes - the SNES context 3170 . its - iteration number 3171 . norm - 2-norm function value (may be estimated) 3172 - mctx - [optional] monitoring context 3173 3174 Options Database Keys: 3175 + -snes_monitor - sets SNESMonitorDefault() 3176 . -snes_monitor_lg_residualnorm - sets line graph monitor, 3177 uses SNESMonitorLGCreate() 3178 - -snes_monitor_cancel - cancels all monitors that have 3179 been hardwired into a code by 3180 calls to SNESMonitorSet(), but 3181 does not cancel those set via 3182 the options database. 3183 3184 Notes: 3185 Several different monitoring routines may be set by calling 3186 SNESMonitorSet() multiple times; all will be called in the 3187 order in which they were set. 3188 3189 Fortran notes: Only a single monitor function can be set for each SNES object 3190 3191 Level: intermediate 3192 3193 .keywords: SNES, nonlinear, set, monitor 3194 3195 .seealso: SNESMonitorDefault(), SNESMonitorCancel() 3196 @*/ 3197 PetscErrorCode SNESMonitorSet(SNES snes,PetscErrorCode (*monitor)(SNES,PetscInt,PetscReal,void*),void *mctx,PetscErrorCode (*monitordestroy)(void**)) 3198 { 3199 PetscInt i; 3200 PetscErrorCode ierr; 3201 3202 PetscFunctionBegin; 3203 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3204 if (snes->numbermonitors >= MAXSNESMONITORS) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Too many monitors set"); 3205 for (i=0; i<snes->numbermonitors;i++) { 3206 if (monitor == snes->monitor[i] && monitordestroy == snes->monitordestroy[i] && mctx == snes->monitorcontext[i]) { 3207 if (monitordestroy) { 3208 ierr = (*monitordestroy)(&mctx);CHKERRQ(ierr); 3209 } 3210 PetscFunctionReturn(0); 3211 } 3212 } 3213 snes->monitor[snes->numbermonitors] = monitor; 3214 snes->monitordestroy[snes->numbermonitors] = monitordestroy; 3215 snes->monitorcontext[snes->numbermonitors++] = (void*)mctx; 3216 PetscFunctionReturn(0); 3217 } 3218 3219 #undef __FUNCT__ 3220 #define __FUNCT__ "SNESMonitorCancel" 3221 /*@C 3222 SNESMonitorCancel - Clears all the monitor functions for a SNES object. 3223 3224 Logically Collective on SNES 3225 3226 Input Parameters: 3227 . snes - the SNES context 3228 3229 Options Database Key: 3230 . -snes_monitor_cancel - cancels all monitors that have been hardwired 3231 into a code by calls to SNESMonitorSet(), but does not cancel those 3232 set via the options database 3233 3234 Notes: 3235 There is no way to clear one specific monitor from a SNES object. 3236 3237 Level: intermediate 3238 3239 .keywords: SNES, nonlinear, set, monitor 3240 3241 .seealso: SNESMonitorDefault(), SNESMonitorSet() 3242 @*/ 3243 PetscErrorCode SNESMonitorCancel(SNES snes) 3244 { 3245 PetscErrorCode ierr; 3246 PetscInt i; 3247 3248 PetscFunctionBegin; 3249 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3250 for (i=0; i<snes->numbermonitors; i++) { 3251 if (snes->monitordestroy[i]) { 3252 ierr = (*snes->monitordestroy[i])(&snes->monitorcontext[i]);CHKERRQ(ierr); 3253 } 3254 } 3255 snes->numbermonitors = 0; 3256 PetscFunctionReturn(0); 3257 } 3258 3259 #undef __FUNCT__ 3260 #define __FUNCT__ "SNESSetConvergenceTest" 3261 /*@C 3262 SNESSetConvergenceTest - Sets the function that is to be used 3263 to test for convergence of the nonlinear iterative solution. 3264 3265 Logically Collective on SNES 3266 3267 Input Parameters: 3268 + snes - the SNES context 3269 . func - routine to test for convergence 3270 . cctx - [optional] context for private data for the convergence routine (may be PETSC_NULL) 3271 - destroy - [optional] destructor for the context (may be PETSC_NULL; PETSC_NULL_FUNCTION in Fortran) 3272 3273 Calling sequence of func: 3274 $ PetscErrorCode func (SNES snes,PetscInt it,PetscReal xnorm,PetscReal gnorm,PetscReal f,SNESConvergedReason *reason,void *cctx) 3275 3276 + snes - the SNES context 3277 . it - current iteration (0 is the first and is before any Newton step) 3278 . cctx - [optional] convergence context 3279 . reason - reason for convergence/divergence 3280 . xnorm - 2-norm of current iterate 3281 . gnorm - 2-norm of current step 3282 - f - 2-norm of function 3283 3284 Level: advanced 3285 3286 .keywords: SNES, nonlinear, set, convergence, test 3287 3288 .seealso: SNESDefaultConverged(), SNESSkipConverged() 3289 @*/ 3290 PetscErrorCode SNESSetConvergenceTest(SNES snes,PetscErrorCode (*func)(SNES,PetscInt,PetscReal,PetscReal,PetscReal,SNESConvergedReason*,void*),void *cctx,PetscErrorCode (*destroy)(void*)) 3291 { 3292 PetscErrorCode ierr; 3293 3294 PetscFunctionBegin; 3295 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3296 if (!func) func = SNESSkipConverged; 3297 if (snes->ops->convergeddestroy) { 3298 ierr = (*snes->ops->convergeddestroy)(snes->cnvP);CHKERRQ(ierr); 3299 } 3300 snes->ops->converged = func; 3301 snes->ops->convergeddestroy = destroy; 3302 snes->cnvP = cctx; 3303 PetscFunctionReturn(0); 3304 } 3305 3306 #undef __FUNCT__ 3307 #define __FUNCT__ "SNESGetConvergedReason" 3308 /*@ 3309 SNESGetConvergedReason - Gets the reason the SNES iteration was stopped. 3310 3311 Not Collective 3312 3313 Input Parameter: 3314 . snes - the SNES context 3315 3316 Output Parameter: 3317 . reason - negative value indicates diverged, positive value converged, see SNESConvergedReason or the 3318 manual pages for the individual convergence tests for complete lists 3319 3320 Level: intermediate 3321 3322 Notes: Can only be called after the call the SNESSolve() is complete. 3323 3324 .keywords: SNES, nonlinear, set, convergence, test 3325 3326 .seealso: SNESSetConvergenceTest(), SNESConvergedReason 3327 @*/ 3328 PetscErrorCode SNESGetConvergedReason(SNES snes,SNESConvergedReason *reason) 3329 { 3330 PetscFunctionBegin; 3331 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3332 PetscValidPointer(reason,2); 3333 *reason = snes->reason; 3334 PetscFunctionReturn(0); 3335 } 3336 3337 #undef __FUNCT__ 3338 #define __FUNCT__ "SNESSetConvergenceHistory" 3339 /*@ 3340 SNESSetConvergenceHistory - Sets the array used to hold the convergence history. 3341 3342 Logically Collective on SNES 3343 3344 Input Parameters: 3345 + snes - iterative context obtained from SNESCreate() 3346 . a - array to hold history, this array will contain the function norms computed at each step 3347 . its - integer array holds the number of linear iterations for each solve. 3348 . na - size of a and its 3349 - reset - PETSC_TRUE indicates each new nonlinear solve resets the history counter to zero, 3350 else it continues storing new values for new nonlinear solves after the old ones 3351 3352 Notes: 3353 If 'a' and 'its' are PETSC_NULL then space is allocated for the history. If 'na' PETSC_DECIDE or PETSC_DEFAULT then a 3354 default array of length 10000 is allocated. 3355 3356 This routine is useful, e.g., when running a code for purposes 3357 of accurate performance monitoring, when no I/O should be done 3358 during the section of code that is being timed. 3359 3360 Level: intermediate 3361 3362 .keywords: SNES, set, convergence, history 3363 3364 .seealso: SNESGetConvergenceHistory() 3365 3366 @*/ 3367 PetscErrorCode SNESSetConvergenceHistory(SNES snes,PetscReal a[],PetscInt its[],PetscInt na,PetscBool reset) 3368 { 3369 PetscErrorCode ierr; 3370 3371 PetscFunctionBegin; 3372 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3373 if (na) PetscValidScalarPointer(a,2); 3374 if (its) PetscValidIntPointer(its,3); 3375 if (na == PETSC_DECIDE || na == PETSC_DEFAULT || !a) { 3376 if (na == PETSC_DECIDE || na == PETSC_DEFAULT) na = 1000; 3377 ierr = PetscMalloc(na*sizeof(PetscReal),&a);CHKERRQ(ierr); 3378 ierr = PetscMalloc(na*sizeof(PetscInt),&its);CHKERRQ(ierr); 3379 snes->conv_malloc = PETSC_TRUE; 3380 } 3381 snes->conv_hist = a; 3382 snes->conv_hist_its = its; 3383 snes->conv_hist_max = na; 3384 snes->conv_hist_len = 0; 3385 snes->conv_hist_reset = reset; 3386 PetscFunctionReturn(0); 3387 } 3388 3389 #if defined(PETSC_HAVE_MATLAB_ENGINE) 3390 #include <engine.h> /* MATLAB include file */ 3391 #include <mex.h> /* MATLAB include file */ 3392 EXTERN_C_BEGIN 3393 #undef __FUNCT__ 3394 #define __FUNCT__ "SNESGetConvergenceHistoryMatlab" 3395 mxArray *SNESGetConvergenceHistoryMatlab(SNES snes) 3396 { 3397 mxArray *mat; 3398 PetscInt i; 3399 PetscReal *ar; 3400 3401 PetscFunctionBegin; 3402 mat = mxCreateDoubleMatrix(snes->conv_hist_len,1,mxREAL); 3403 ar = (PetscReal*) mxGetData(mat); 3404 for (i=0; i<snes->conv_hist_len; i++) { 3405 ar[i] = snes->conv_hist[i]; 3406 } 3407 PetscFunctionReturn(mat); 3408 } 3409 EXTERN_C_END 3410 #endif 3411 3412 3413 #undef __FUNCT__ 3414 #define __FUNCT__ "SNESGetConvergenceHistory" 3415 /*@C 3416 SNESGetConvergenceHistory - Gets the array used to hold the convergence history. 3417 3418 Not Collective 3419 3420 Input Parameter: 3421 . snes - iterative context obtained from SNESCreate() 3422 3423 Output Parameters: 3424 . a - array to hold history 3425 . its - integer array holds the number of linear iterations (or 3426 negative if not converged) for each solve. 3427 - na - size of a and its 3428 3429 Notes: 3430 The calling sequence for this routine in Fortran is 3431 $ call SNESGetConvergenceHistory(SNES snes, integer na, integer ierr) 3432 3433 This routine is useful, e.g., when running a code for purposes 3434 of accurate performance monitoring, when no I/O should be done 3435 during the section of code that is being timed. 3436 3437 Level: intermediate 3438 3439 .keywords: SNES, get, convergence, history 3440 3441 .seealso: SNESSetConvergencHistory() 3442 3443 @*/ 3444 PetscErrorCode SNESGetConvergenceHistory(SNES snes,PetscReal *a[],PetscInt *its[],PetscInt *na) 3445 { 3446 PetscFunctionBegin; 3447 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3448 if (a) *a = snes->conv_hist; 3449 if (its) *its = snes->conv_hist_its; 3450 if (na) *na = snes->conv_hist_len; 3451 PetscFunctionReturn(0); 3452 } 3453 3454 #undef __FUNCT__ 3455 #define __FUNCT__ "SNESSetUpdate" 3456 /*@C 3457 SNESSetUpdate - Sets the general-purpose update function called 3458 at the beginning of every iteration of the nonlinear solve. Specifically 3459 it is called just before the Jacobian is "evaluated". 3460 3461 Logically Collective on SNES 3462 3463 Input Parameters: 3464 . snes - The nonlinear solver context 3465 . func - The function 3466 3467 Calling sequence of func: 3468 . func (SNES snes, PetscInt step); 3469 3470 . step - The current step of the iteration 3471 3472 Level: advanced 3473 3474 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() 3475 This is not used by most users. 3476 3477 .keywords: SNES, update 3478 3479 .seealso SNESDefaultUpdate(), SNESSetJacobian(), SNESSolve() 3480 @*/ 3481 PetscErrorCode SNESSetUpdate(SNES snes, PetscErrorCode (*func)(SNES, PetscInt)) 3482 { 3483 PetscFunctionBegin; 3484 PetscValidHeaderSpecific(snes, SNES_CLASSID,1); 3485 snes->ops->update = func; 3486 PetscFunctionReturn(0); 3487 } 3488 3489 #undef __FUNCT__ 3490 #define __FUNCT__ "SNESDefaultUpdate" 3491 /*@ 3492 SNESDefaultUpdate - The default update function which does nothing. 3493 3494 Not collective 3495 3496 Input Parameters: 3497 . snes - The nonlinear solver context 3498 . step - The current step of the iteration 3499 3500 Level: intermediate 3501 3502 .keywords: SNES, update 3503 .seealso SNESSetUpdate(), SNESDefaultRhsBC(), SNESDefaultShortolutionBC() 3504 @*/ 3505 PetscErrorCode SNESDefaultUpdate(SNES snes, PetscInt step) 3506 { 3507 PetscFunctionBegin; 3508 PetscFunctionReturn(0); 3509 } 3510 3511 #undef __FUNCT__ 3512 #define __FUNCT__ "SNESScaleStep_Private" 3513 /* 3514 SNESScaleStep_Private - Scales a step so that its length is less than the 3515 positive parameter delta. 3516 3517 Input Parameters: 3518 + snes - the SNES context 3519 . y - approximate solution of linear system 3520 . fnorm - 2-norm of current function 3521 - delta - trust region size 3522 3523 Output Parameters: 3524 + gpnorm - predicted function norm at the new point, assuming local 3525 linearization. The value is zero if the step lies within the trust 3526 region, and exceeds zero otherwise. 3527 - ynorm - 2-norm of the step 3528 3529 Note: 3530 For non-trust region methods such as SNESLS, the parameter delta 3531 is set to be the maximum allowable step size. 3532 3533 .keywords: SNES, nonlinear, scale, step 3534 */ 3535 PetscErrorCode SNESScaleStep_Private(SNES snes,Vec y,PetscReal *fnorm,PetscReal *delta,PetscReal *gpnorm,PetscReal *ynorm) 3536 { 3537 PetscReal nrm; 3538 PetscScalar cnorm; 3539 PetscErrorCode ierr; 3540 3541 PetscFunctionBegin; 3542 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3543 PetscValidHeaderSpecific(y,VEC_CLASSID,2); 3544 PetscCheckSameComm(snes,1,y,2); 3545 3546 ierr = VecNorm(y,NORM_2,&nrm);CHKERRQ(ierr); 3547 if (nrm > *delta) { 3548 nrm = *delta/nrm; 3549 *gpnorm = (1.0 - nrm)*(*fnorm); 3550 cnorm = nrm; 3551 ierr = VecScale(y,cnorm);CHKERRQ(ierr); 3552 *ynorm = *delta; 3553 } else { 3554 *gpnorm = 0.0; 3555 *ynorm = nrm; 3556 } 3557 PetscFunctionReturn(0); 3558 } 3559 3560 #undef __FUNCT__ 3561 #define __FUNCT__ "SNESSolve" 3562 /*@C 3563 SNESSolve - Solves a nonlinear system F(x) = b. 3564 Call SNESSolve() after calling SNESCreate() and optional routines of the form SNESSetXXX(). 3565 3566 Collective on SNES 3567 3568 Input Parameters: 3569 + snes - the SNES context 3570 . b - the constant part of the equation F(x) = b, or PETSC_NULL to use zero. 3571 - x - the solution vector. 3572 3573 Notes: 3574 The user should initialize the vector,x, with the initial guess 3575 for the nonlinear solve prior to calling SNESSolve. In particular, 3576 to employ an initial guess of zero, the user should explicitly set 3577 this vector to zero by calling VecSet(). 3578 3579 Level: beginner 3580 3581 .keywords: SNES, nonlinear, solve 3582 3583 .seealso: SNESCreate(), SNESDestroy(), SNESSetFunction(), SNESSetJacobian(), SNESSetGridSequence(), SNESGetSolution() 3584 @*/ 3585 PetscErrorCode SNESSolve(SNES snes,Vec b,Vec x) 3586 { 3587 PetscErrorCode ierr; 3588 PetscBool flg; 3589 char filename[PETSC_MAX_PATH_LEN]; 3590 PetscViewer viewer; 3591 PetscInt grid; 3592 Vec xcreated = PETSC_NULL; 3593 DM dm; 3594 3595 PetscFunctionBegin; 3596 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3597 if (x) PetscValidHeaderSpecific(x,VEC_CLASSID,3); 3598 if (x) PetscCheckSameComm(snes,1,x,3); 3599 if (b) PetscValidHeaderSpecific(b,VEC_CLASSID,2); 3600 if (b) PetscCheckSameComm(snes,1,b,2); 3601 3602 if (!x) { 3603 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 3604 ierr = DMCreateGlobalVector(dm,&xcreated);CHKERRQ(ierr); 3605 x = xcreated; 3606 } 3607 3608 for (grid=0; grid<snes->gridsequence; grid++) {ierr = PetscViewerASCIIPushTab(PETSC_VIEWER_STDOUT_(((PetscObject)snes)->comm));CHKERRQ(ierr);} 3609 for (grid=0; grid<snes->gridsequence+1; grid++) { 3610 3611 /* set solution vector */ 3612 if (!grid) {ierr = PetscObjectReference((PetscObject)x);CHKERRQ(ierr);} 3613 ierr = VecDestroy(&snes->vec_sol);CHKERRQ(ierr); 3614 snes->vec_sol = x; 3615 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 3616 3617 /* set affine vector if provided */ 3618 if (b) { ierr = PetscObjectReference((PetscObject)b);CHKERRQ(ierr); } 3619 ierr = VecDestroy(&snes->vec_rhs);CHKERRQ(ierr); 3620 snes->vec_rhs = b; 3621 3622 ierr = SNESSetUp(snes);CHKERRQ(ierr); 3623 3624 if (!grid) { 3625 if (snes->ops->computeinitialguess) { 3626 ierr = (*snes->ops->computeinitialguess)(snes,snes->vec_sol,snes->initialguessP);CHKERRQ(ierr); 3627 } 3628 } 3629 3630 if (snes->conv_hist_reset) snes->conv_hist_len = 0; 3631 snes->nfuncs = 0; snes->linear_its = 0; snes->numFailures = 0; 3632 3633 ierr = PetscLogEventBegin(SNES_Solve,snes,0,0,0);CHKERRQ(ierr); 3634 ierr = (*snes->ops->solve)(snes);CHKERRQ(ierr); 3635 ierr = PetscLogEventEnd(SNES_Solve,snes,0,0,0);CHKERRQ(ierr); 3636 if (snes->domainerror){ 3637 snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; 3638 snes->domainerror = PETSC_FALSE; 3639 } 3640 if (!snes->reason) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Internal error, solver returned without setting converged reason"); 3641 3642 flg = PETSC_FALSE; 3643 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_test_local_min",&flg,PETSC_NULL);CHKERRQ(ierr); 3644 if (flg && !PetscPreLoadingOn) { ierr = SNESTestLocalMin(snes);CHKERRQ(ierr); } 3645 if (snes->printreason) { 3646 ierr = PetscViewerASCIIAddTab(PETSC_VIEWER_STDOUT_(((PetscObject)snes)->comm),((PetscObject)snes)->tablevel);CHKERRQ(ierr); 3647 if (snes->reason > 0) { 3648 ierr = PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_(((PetscObject)snes)->comm),"Nonlinear solve converged due to %s iterations %D\n",SNESConvergedReasons[snes->reason],snes->iter);CHKERRQ(ierr); 3649 } else { 3650 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); 3651 } 3652 ierr = PetscViewerASCIISubtractTab(PETSC_VIEWER_STDOUT_(((PetscObject)snes)->comm),((PetscObject)snes)->tablevel);CHKERRQ(ierr); 3653 } 3654 3655 if (snes->errorifnotconverged && snes->reason < 0) SETERRQ(((PetscObject)snes)->comm,PETSC_ERR_NOT_CONVERGED,"SNESSolve has not converged"); 3656 if (grid < snes->gridsequence) { 3657 DM fine; 3658 Vec xnew; 3659 Mat interp; 3660 3661 ierr = DMRefine(snes->dm,((PetscObject)snes)->comm,&fine);CHKERRQ(ierr); 3662 if (!fine) SETERRQ(((PetscObject)snes)->comm,PETSC_ERR_ARG_INCOMP,"DMRefine() did not perform any refinement, cannot continue grid sequencing"); 3663 ierr = DMCreateInterpolation(snes->dm,fine,&interp,PETSC_NULL);CHKERRQ(ierr); 3664 ierr = DMCreateGlobalVector(fine,&xnew);CHKERRQ(ierr); 3665 ierr = MatInterpolate(interp,x,xnew);CHKERRQ(ierr); 3666 ierr = DMInterpolate(snes->dm,interp,fine);CHKERRQ(ierr); 3667 ierr = MatDestroy(&interp);CHKERRQ(ierr); 3668 x = xnew; 3669 3670 ierr = SNESReset(snes);CHKERRQ(ierr); 3671 ierr = SNESSetDM(snes,fine);CHKERRQ(ierr); 3672 ierr = DMDestroy(&fine);CHKERRQ(ierr); 3673 ierr = PetscViewerASCIIPopTab(PETSC_VIEWER_STDOUT_(((PetscObject)snes)->comm));CHKERRQ(ierr); 3674 } 3675 } 3676 /* monitoring and viewing */ 3677 flg = PETSC_FALSE; 3678 ierr = PetscOptionsGetString(((PetscObject)snes)->prefix,"-snes_view",filename,PETSC_MAX_PATH_LEN,&flg);CHKERRQ(ierr); 3679 if (flg && !PetscPreLoadingOn) { 3680 ierr = PetscViewerASCIIOpen(((PetscObject)snes)->comm,filename,&viewer);CHKERRQ(ierr); 3681 ierr = SNESView(snes,viewer);CHKERRQ(ierr); 3682 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 3683 } 3684 3685 flg = PETSC_FALSE; 3686 ierr = PetscOptionsGetString(((PetscObject)snes)->prefix,"-snes_view_solution_vtk",filename,PETSC_MAX_PATH_LEN,&flg);CHKERRQ(ierr); 3687 if (flg) { 3688 PetscViewer viewer; 3689 ierr = PetscViewerCreate(((PetscObject)snes)->comm,&viewer);CHKERRQ(ierr); 3690 ierr = PetscViewerSetType(viewer,PETSCVIEWERVTK);CHKERRQ(ierr); 3691 ierr = PetscViewerFileSetName(viewer,filename);CHKERRQ(ierr); 3692 ierr = VecView(snes->vec_sol,viewer);CHKERRQ(ierr); 3693 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 3694 } 3695 3696 ierr = VecDestroy(&xcreated);CHKERRQ(ierr); 3697 PetscFunctionReturn(0); 3698 } 3699 3700 /* --------- Internal routines for SNES Package --------- */ 3701 3702 #undef __FUNCT__ 3703 #define __FUNCT__ "SNESSetType" 3704 /*@C 3705 SNESSetType - Sets the method for the nonlinear solver. 3706 3707 Collective on SNES 3708 3709 Input Parameters: 3710 + snes - the SNES context 3711 - type - a known method 3712 3713 Options Database Key: 3714 . -snes_type <type> - Sets the method; use -help for a list 3715 of available methods (for instance, ls or tr) 3716 3717 Notes: 3718 See "petsc/include/petscsnes.h" for available methods (for instance) 3719 + SNESLS - Newton's method with line search 3720 (systems of nonlinear equations) 3721 . SNESTR - Newton's method with trust region 3722 (systems of nonlinear equations) 3723 3724 Normally, it is best to use the SNESSetFromOptions() command and then 3725 set the SNES solver type from the options database rather than by using 3726 this routine. Using the options database provides the user with 3727 maximum flexibility in evaluating the many nonlinear solvers. 3728 The SNESSetType() routine is provided for those situations where it 3729 is necessary to set the nonlinear solver independently of the command 3730 line or options database. This might be the case, for example, when 3731 the choice of solver changes during the execution of the program, 3732 and the user's application is taking responsibility for choosing the 3733 appropriate method. 3734 3735 Level: intermediate 3736 3737 .keywords: SNES, set, type 3738 3739 .seealso: SNESType, SNESCreate() 3740 3741 @*/ 3742 PetscErrorCode SNESSetType(SNES snes,SNESType type) 3743 { 3744 PetscErrorCode ierr,(*r)(SNES); 3745 PetscBool match; 3746 3747 PetscFunctionBegin; 3748 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3749 PetscValidCharPointer(type,2); 3750 3751 ierr = PetscObjectTypeCompare((PetscObject)snes,type,&match);CHKERRQ(ierr); 3752 if (match) PetscFunctionReturn(0); 3753 3754 ierr = PetscFListFind(SNESList,((PetscObject)snes)->comm,type,PETSC_TRUE,(void (**)(void)) &r);CHKERRQ(ierr); 3755 if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unable to find requested SNES type %s",type); 3756 /* Destroy the previous private SNES context */ 3757 if (snes->ops->destroy) { 3758 ierr = (*(snes)->ops->destroy)(snes);CHKERRQ(ierr); 3759 snes->ops->destroy = PETSC_NULL; 3760 } 3761 /* Reinitialize function pointers in SNESOps structure */ 3762 snes->ops->setup = 0; 3763 snes->ops->solve = 0; 3764 snes->ops->view = 0; 3765 snes->ops->setfromoptions = 0; 3766 snes->ops->destroy = 0; 3767 /* Call the SNESCreate_XXX routine for this particular Nonlinear solver */ 3768 snes->setupcalled = PETSC_FALSE; 3769 ierr = PetscObjectChangeTypeName((PetscObject)snes,type);CHKERRQ(ierr); 3770 ierr = (*r)(snes);CHKERRQ(ierr); 3771 #if defined(PETSC_HAVE_AMS) 3772 if (PetscAMSPublishAll) { 3773 ierr = PetscObjectAMSPublish((PetscObject)snes);CHKERRQ(ierr); 3774 } 3775 #endif 3776 PetscFunctionReturn(0); 3777 } 3778 3779 3780 /* --------------------------------------------------------------------- */ 3781 #undef __FUNCT__ 3782 #define __FUNCT__ "SNESRegisterDestroy" 3783 /*@ 3784 SNESRegisterDestroy - Frees the list of nonlinear solvers that were 3785 registered by SNESRegisterDynamic(). 3786 3787 Not Collective 3788 3789 Level: advanced 3790 3791 .keywords: SNES, nonlinear, register, destroy 3792 3793 .seealso: SNESRegisterAll(), SNESRegisterAll() 3794 @*/ 3795 PetscErrorCode SNESRegisterDestroy(void) 3796 { 3797 PetscErrorCode ierr; 3798 3799 PetscFunctionBegin; 3800 ierr = PetscFListDestroy(&SNESList);CHKERRQ(ierr); 3801 SNESRegisterAllCalled = PETSC_FALSE; 3802 PetscFunctionReturn(0); 3803 } 3804 3805 #undef __FUNCT__ 3806 #define __FUNCT__ "SNESGetType" 3807 /*@C 3808 SNESGetType - Gets the SNES method type and name (as a string). 3809 3810 Not Collective 3811 3812 Input Parameter: 3813 . snes - nonlinear solver context 3814 3815 Output Parameter: 3816 . type - SNES method (a character string) 3817 3818 Level: intermediate 3819 3820 .keywords: SNES, nonlinear, get, type, name 3821 @*/ 3822 PetscErrorCode SNESGetType(SNES snes,SNESType *type) 3823 { 3824 PetscFunctionBegin; 3825 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3826 PetscValidPointer(type,2); 3827 *type = ((PetscObject)snes)->type_name; 3828 PetscFunctionReturn(0); 3829 } 3830 3831 #undef __FUNCT__ 3832 #define __FUNCT__ "SNESGetSolution" 3833 /*@ 3834 SNESGetSolution - Returns the vector where the approximate solution is 3835 stored. This is the fine grid solution when using SNESSetGridSequence(). 3836 3837 Not Collective, but Vec is parallel if SNES is parallel 3838 3839 Input Parameter: 3840 . snes - the SNES context 3841 3842 Output Parameter: 3843 . x - the solution 3844 3845 Level: intermediate 3846 3847 .keywords: SNES, nonlinear, get, solution 3848 3849 .seealso: SNESGetSolutionUpdate(), SNESGetFunction() 3850 @*/ 3851 PetscErrorCode SNESGetSolution(SNES snes,Vec *x) 3852 { 3853 PetscFunctionBegin; 3854 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3855 PetscValidPointer(x,2); 3856 *x = snes->vec_sol; 3857 PetscFunctionReturn(0); 3858 } 3859 3860 #undef __FUNCT__ 3861 #define __FUNCT__ "SNESGetSolutionUpdate" 3862 /*@ 3863 SNESGetSolutionUpdate - Returns the vector where the solution update is 3864 stored. 3865 3866 Not Collective, but Vec is parallel if SNES is parallel 3867 3868 Input Parameter: 3869 . snes - the SNES context 3870 3871 Output Parameter: 3872 . x - the solution update 3873 3874 Level: advanced 3875 3876 .keywords: SNES, nonlinear, get, solution, update 3877 3878 .seealso: SNESGetSolution(), SNESGetFunction() 3879 @*/ 3880 PetscErrorCode SNESGetSolutionUpdate(SNES snes,Vec *x) 3881 { 3882 PetscFunctionBegin; 3883 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3884 PetscValidPointer(x,2); 3885 *x = snes->vec_sol_update; 3886 PetscFunctionReturn(0); 3887 } 3888 3889 #undef __FUNCT__ 3890 #define __FUNCT__ "SNESGetFunction" 3891 /*@C 3892 SNESGetFunction - Returns the vector where the function is stored. 3893 3894 Not Collective, but Vec is parallel if SNES is parallel. Collective if Vec is requested, but has not been created yet. 3895 3896 Input Parameter: 3897 . snes - the SNES context 3898 3899 Output Parameter: 3900 + r - the function (or PETSC_NULL) 3901 . func - the function (or PETSC_NULL) 3902 - ctx - the function context (or PETSC_NULL) 3903 3904 Level: advanced 3905 3906 .keywords: SNES, nonlinear, get, function 3907 3908 .seealso: SNESSetFunction(), SNESGetSolution() 3909 @*/ 3910 PetscErrorCode SNESGetFunction(SNES snes,Vec *r,PetscErrorCode (**func)(SNES,Vec,Vec,void*),void **ctx) 3911 { 3912 PetscErrorCode ierr; 3913 DM dm; 3914 3915 PetscFunctionBegin; 3916 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3917 if (r) { 3918 if (!snes->vec_func) { 3919 if (snes->vec_rhs) { 3920 ierr = VecDuplicate(snes->vec_rhs,&snes->vec_func);CHKERRQ(ierr); 3921 } else if (snes->vec_sol) { 3922 ierr = VecDuplicate(snes->vec_sol,&snes->vec_func);CHKERRQ(ierr); 3923 } else if (snes->dm) { 3924 ierr = DMCreateGlobalVector(snes->dm,&snes->vec_func);CHKERRQ(ierr); 3925 } 3926 } 3927 *r = snes->vec_func; 3928 } 3929 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 3930 ierr = DMSNESGetFunction(dm,func,ctx);CHKERRQ(ierr); 3931 PetscFunctionReturn(0); 3932 } 3933 3934 /*@C 3935 SNESGetGS - Returns the GS function and context. 3936 3937 Input Parameter: 3938 . snes - the SNES context 3939 3940 Output Parameter: 3941 + gsfunc - the function (or PETSC_NULL) 3942 - ctx - the function context (or PETSC_NULL) 3943 3944 Level: advanced 3945 3946 .keywords: SNES, nonlinear, get, function 3947 3948 .seealso: SNESSetGS(), SNESGetFunction() 3949 @*/ 3950 3951 #undef __FUNCT__ 3952 #define __FUNCT__ "SNESGetGS" 3953 PetscErrorCode SNESGetGS (SNES snes, PetscErrorCode(**func)(SNES, Vec, Vec, void*), void ** ctx) 3954 { 3955 PetscErrorCode ierr; 3956 DM dm; 3957 3958 PetscFunctionBegin; 3959 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3960 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 3961 ierr = DMSNESGetGS(dm,func,ctx);CHKERRQ(ierr); 3962 PetscFunctionReturn(0); 3963 } 3964 3965 #undef __FUNCT__ 3966 #define __FUNCT__ "SNESSetOptionsPrefix" 3967 /*@C 3968 SNESSetOptionsPrefix - Sets the prefix used for searching for all 3969 SNES options in the database. 3970 3971 Logically Collective on SNES 3972 3973 Input Parameter: 3974 + snes - the SNES context 3975 - prefix - the prefix to prepend to all option names 3976 3977 Notes: 3978 A hyphen (-) must NOT be given at the beginning of the prefix name. 3979 The first character of all runtime options is AUTOMATICALLY the hyphen. 3980 3981 Level: advanced 3982 3983 .keywords: SNES, set, options, prefix, database 3984 3985 .seealso: SNESSetFromOptions() 3986 @*/ 3987 PetscErrorCode SNESSetOptionsPrefix(SNES snes,const char prefix[]) 3988 { 3989 PetscErrorCode ierr; 3990 3991 PetscFunctionBegin; 3992 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3993 ierr = PetscObjectSetOptionsPrefix((PetscObject)snes,prefix);CHKERRQ(ierr); 3994 if (!snes->ksp) {ierr = SNESGetKSP(snes,&snes->ksp);CHKERRQ(ierr);} 3995 if (snes->linesearch) { 3996 ierr = SNESGetSNESLineSearch(snes,&snes->linesearch);CHKERRQ(ierr); 3997 ierr = PetscObjectSetOptionsPrefix((PetscObject)snes->linesearch,prefix);CHKERRQ(ierr); 3998 } 3999 ierr = KSPSetOptionsPrefix(snes->ksp,prefix);CHKERRQ(ierr); 4000 PetscFunctionReturn(0); 4001 } 4002 4003 #undef __FUNCT__ 4004 #define __FUNCT__ "SNESAppendOptionsPrefix" 4005 /*@C 4006 SNESAppendOptionsPrefix - Appends to the prefix used for searching for all 4007 SNES options in the database. 4008 4009 Logically Collective on SNES 4010 4011 Input Parameters: 4012 + snes - the SNES context 4013 - prefix - the prefix to prepend to all option names 4014 4015 Notes: 4016 A hyphen (-) must NOT be given at the beginning of the prefix name. 4017 The first character of all runtime options is AUTOMATICALLY the hyphen. 4018 4019 Level: advanced 4020 4021 .keywords: SNES, append, options, prefix, database 4022 4023 .seealso: SNESGetOptionsPrefix() 4024 @*/ 4025 PetscErrorCode SNESAppendOptionsPrefix(SNES snes,const char prefix[]) 4026 { 4027 PetscErrorCode ierr; 4028 4029 PetscFunctionBegin; 4030 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4031 ierr = PetscObjectAppendOptionsPrefix((PetscObject)snes,prefix);CHKERRQ(ierr); 4032 if (!snes->ksp) {ierr = SNESGetKSP(snes,&snes->ksp);CHKERRQ(ierr);} 4033 if (snes->linesearch) { 4034 ierr = SNESGetSNESLineSearch(snes,&snes->linesearch);CHKERRQ(ierr); 4035 ierr = PetscObjectAppendOptionsPrefix((PetscObject)snes->linesearch,prefix);CHKERRQ(ierr); 4036 } 4037 ierr = KSPAppendOptionsPrefix(snes->ksp,prefix);CHKERRQ(ierr); 4038 PetscFunctionReturn(0); 4039 } 4040 4041 #undef __FUNCT__ 4042 #define __FUNCT__ "SNESGetOptionsPrefix" 4043 /*@C 4044 SNESGetOptionsPrefix - Sets the prefix used for searching for all 4045 SNES options in the database. 4046 4047 Not Collective 4048 4049 Input Parameter: 4050 . snes - the SNES context 4051 4052 Output Parameter: 4053 . prefix - pointer to the prefix string used 4054 4055 Notes: On the fortran side, the user should pass in a string 'prefix' of 4056 sufficient length to hold the prefix. 4057 4058 Level: advanced 4059 4060 .keywords: SNES, get, options, prefix, database 4061 4062 .seealso: SNESAppendOptionsPrefix() 4063 @*/ 4064 PetscErrorCode SNESGetOptionsPrefix(SNES snes,const char *prefix[]) 4065 { 4066 PetscErrorCode ierr; 4067 4068 PetscFunctionBegin; 4069 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4070 ierr = PetscObjectGetOptionsPrefix((PetscObject)snes,prefix);CHKERRQ(ierr); 4071 PetscFunctionReturn(0); 4072 } 4073 4074 4075 #undef __FUNCT__ 4076 #define __FUNCT__ "SNESRegister" 4077 /*@C 4078 SNESRegister - See SNESRegisterDynamic() 4079 4080 Level: advanced 4081 @*/ 4082 PetscErrorCode SNESRegister(const char sname[],const char path[],const char name[],PetscErrorCode (*function)(SNES)) 4083 { 4084 char fullname[PETSC_MAX_PATH_LEN]; 4085 PetscErrorCode ierr; 4086 4087 PetscFunctionBegin; 4088 ierr = PetscFListConcat(path,name,fullname);CHKERRQ(ierr); 4089 ierr = PetscFListAdd(&SNESList,sname,fullname,(void (*)(void))function);CHKERRQ(ierr); 4090 PetscFunctionReturn(0); 4091 } 4092 4093 #undef __FUNCT__ 4094 #define __FUNCT__ "SNESTestLocalMin" 4095 PetscErrorCode SNESTestLocalMin(SNES snes) 4096 { 4097 PetscErrorCode ierr; 4098 PetscInt N,i,j; 4099 Vec u,uh,fh; 4100 PetscScalar value; 4101 PetscReal norm; 4102 4103 PetscFunctionBegin; 4104 ierr = SNESGetSolution(snes,&u);CHKERRQ(ierr); 4105 ierr = VecDuplicate(u,&uh);CHKERRQ(ierr); 4106 ierr = VecDuplicate(u,&fh);CHKERRQ(ierr); 4107 4108 /* currently only works for sequential */ 4109 ierr = PetscPrintf(PETSC_COMM_WORLD,"Testing FormFunction() for local min\n"); 4110 ierr = VecGetSize(u,&N);CHKERRQ(ierr); 4111 for (i=0; i<N; i++) { 4112 ierr = VecCopy(u,uh);CHKERRQ(ierr); 4113 ierr = PetscPrintf(PETSC_COMM_WORLD,"i = %D\n",i);CHKERRQ(ierr); 4114 for (j=-10; j<11; j++) { 4115 value = PetscSign(j)*exp(PetscAbs(j)-10.0); 4116 ierr = VecSetValue(uh,i,value,ADD_VALUES);CHKERRQ(ierr); 4117 ierr = SNESComputeFunction(snes,uh,fh);CHKERRQ(ierr); 4118 ierr = VecNorm(fh,NORM_2,&norm);CHKERRQ(ierr); 4119 ierr = PetscPrintf(PETSC_COMM_WORLD," j norm %D %18.16e\n",j,norm);CHKERRQ(ierr); 4120 value = -value; 4121 ierr = VecSetValue(uh,i,value,ADD_VALUES);CHKERRQ(ierr); 4122 } 4123 } 4124 ierr = VecDestroy(&uh);CHKERRQ(ierr); 4125 ierr = VecDestroy(&fh);CHKERRQ(ierr); 4126 PetscFunctionReturn(0); 4127 } 4128 4129 #undef __FUNCT__ 4130 #define __FUNCT__ "SNESKSPSetUseEW" 4131 /*@ 4132 SNESKSPSetUseEW - Sets SNES use Eisenstat-Walker method for 4133 computing relative tolerance for linear solvers within an inexact 4134 Newton method. 4135 4136 Logically Collective on SNES 4137 4138 Input Parameters: 4139 + snes - SNES context 4140 - flag - PETSC_TRUE or PETSC_FALSE 4141 4142 Options Database: 4143 + -snes_ksp_ew - use Eisenstat-Walker method for determining linear system convergence 4144 . -snes_ksp_ew_version ver - version of Eisenstat-Walker method 4145 . -snes_ksp_ew_rtol0 <rtol0> - Sets rtol0 4146 . -snes_ksp_ew_rtolmax <rtolmax> - Sets rtolmax 4147 . -snes_ksp_ew_gamma <gamma> - Sets gamma 4148 . -snes_ksp_ew_alpha <alpha> - Sets alpha 4149 . -snes_ksp_ew_alpha2 <alpha2> - Sets alpha2 4150 - -snes_ksp_ew_threshold <threshold> - Sets threshold 4151 4152 Notes: 4153 Currently, the default is to use a constant relative tolerance for 4154 the inner linear solvers. Alternatively, one can use the 4155 Eisenstat-Walker method, where the relative convergence tolerance 4156 is reset at each Newton iteration according progress of the nonlinear 4157 solver. 4158 4159 Level: advanced 4160 4161 Reference: 4162 S. C. Eisenstat and H. F. Walker, "Choosing the forcing terms in an 4163 inexact Newton method", SISC 17 (1), pp.16-32, 1996. 4164 4165 .keywords: SNES, KSP, Eisenstat, Walker, convergence, test, inexact, Newton 4166 4167 .seealso: SNESKSPGetUseEW(), SNESKSPGetParametersEW(), SNESKSPSetParametersEW() 4168 @*/ 4169 PetscErrorCode SNESKSPSetUseEW(SNES snes,PetscBool flag) 4170 { 4171 PetscFunctionBegin; 4172 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4173 PetscValidLogicalCollectiveBool(snes,flag,2); 4174 snes->ksp_ewconv = flag; 4175 PetscFunctionReturn(0); 4176 } 4177 4178 #undef __FUNCT__ 4179 #define __FUNCT__ "SNESKSPGetUseEW" 4180 /*@ 4181 SNESKSPGetUseEW - Gets if SNES is using Eisenstat-Walker method 4182 for computing relative tolerance for linear solvers within an 4183 inexact Newton method. 4184 4185 Not Collective 4186 4187 Input Parameter: 4188 . snes - SNES context 4189 4190 Output Parameter: 4191 . flag - PETSC_TRUE or PETSC_FALSE 4192 4193 Notes: 4194 Currently, the default is to use a constant relative tolerance for 4195 the inner linear solvers. Alternatively, one can use the 4196 Eisenstat-Walker method, where the relative convergence tolerance 4197 is reset at each Newton iteration according progress of the nonlinear 4198 solver. 4199 4200 Level: advanced 4201 4202 Reference: 4203 S. C. Eisenstat and H. F. Walker, "Choosing the forcing terms in an 4204 inexact Newton method", SISC 17 (1), pp.16-32, 1996. 4205 4206 .keywords: SNES, KSP, Eisenstat, Walker, convergence, test, inexact, Newton 4207 4208 .seealso: SNESKSPSetUseEW(), SNESKSPGetParametersEW(), SNESKSPSetParametersEW() 4209 @*/ 4210 PetscErrorCode SNESKSPGetUseEW(SNES snes, PetscBool *flag) 4211 { 4212 PetscFunctionBegin; 4213 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4214 PetscValidPointer(flag,2); 4215 *flag = snes->ksp_ewconv; 4216 PetscFunctionReturn(0); 4217 } 4218 4219 #undef __FUNCT__ 4220 #define __FUNCT__ "SNESKSPSetParametersEW" 4221 /*@ 4222 SNESKSPSetParametersEW - Sets parameters for Eisenstat-Walker 4223 convergence criteria for the linear solvers within an inexact 4224 Newton method. 4225 4226 Logically Collective on SNES 4227 4228 Input Parameters: 4229 + snes - SNES context 4230 . version - version 1, 2 (default is 2) or 3 4231 . rtol_0 - initial relative tolerance (0 <= rtol_0 < 1) 4232 . rtol_max - maximum relative tolerance (0 <= rtol_max < 1) 4233 . gamma - multiplicative factor for version 2 rtol computation 4234 (0 <= gamma2 <= 1) 4235 . alpha - power for version 2 rtol computation (1 < alpha <= 2) 4236 . alpha2 - power for safeguard 4237 - threshold - threshold for imposing safeguard (0 < threshold < 1) 4238 4239 Note: 4240 Version 3 was contributed by Luis Chacon, June 2006. 4241 4242 Use PETSC_DEFAULT to retain the default for any of the parameters. 4243 4244 Level: advanced 4245 4246 Reference: 4247 S. C. Eisenstat and H. F. Walker, "Choosing the forcing terms in an 4248 inexact Newton method", Utah State University Math. Stat. Dept. Res. 4249 Report 6/94/75, June, 1994, to appear in SIAM J. Sci. Comput. 4250 4251 .keywords: SNES, KSP, Eisenstat, Walker, set, parameters 4252 4253 .seealso: SNESKSPSetUseEW(), SNESKSPGetUseEW(), SNESKSPGetParametersEW() 4254 @*/ 4255 PetscErrorCode SNESKSPSetParametersEW(SNES snes,PetscInt version,PetscReal rtol_0,PetscReal rtol_max, 4256 PetscReal gamma,PetscReal alpha,PetscReal alpha2,PetscReal threshold) 4257 { 4258 SNESKSPEW *kctx; 4259 PetscFunctionBegin; 4260 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4261 kctx = (SNESKSPEW*)snes->kspconvctx; 4262 if (!kctx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"No Eisenstat-Walker context existing"); 4263 PetscValidLogicalCollectiveInt(snes,version,2); 4264 PetscValidLogicalCollectiveReal(snes,rtol_0,3); 4265 PetscValidLogicalCollectiveReal(snes,rtol_max,4); 4266 PetscValidLogicalCollectiveReal(snes,gamma,5); 4267 PetscValidLogicalCollectiveReal(snes,alpha,6); 4268 PetscValidLogicalCollectiveReal(snes,alpha2,7); 4269 PetscValidLogicalCollectiveReal(snes,threshold,8); 4270 4271 if (version != PETSC_DEFAULT) kctx->version = version; 4272 if (rtol_0 != PETSC_DEFAULT) kctx->rtol_0 = rtol_0; 4273 if (rtol_max != PETSC_DEFAULT) kctx->rtol_max = rtol_max; 4274 if (gamma != PETSC_DEFAULT) kctx->gamma = gamma; 4275 if (alpha != PETSC_DEFAULT) kctx->alpha = alpha; 4276 if (alpha2 != PETSC_DEFAULT) kctx->alpha2 = alpha2; 4277 if (threshold != PETSC_DEFAULT) kctx->threshold = threshold; 4278 4279 if (kctx->version < 1 || kctx->version > 3) { 4280 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Only versions 1, 2 and 3 are supported: %D",kctx->version); 4281 } 4282 if (kctx->rtol_0 < 0.0 || kctx->rtol_0 >= 1.0) { 4283 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"0.0 <= rtol_0 < 1.0: %G",kctx->rtol_0); 4284 } 4285 if (kctx->rtol_max < 0.0 || kctx->rtol_max >= 1.0) { 4286 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"0.0 <= rtol_max (%G) < 1.0\n",kctx->rtol_max); 4287 } 4288 if (kctx->gamma < 0.0 || kctx->gamma > 1.0) { 4289 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"0.0 <= gamma (%G) <= 1.0\n",kctx->gamma); 4290 } 4291 if (kctx->alpha <= 1.0 || kctx->alpha > 2.0) { 4292 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"1.0 < alpha (%G) <= 2.0\n",kctx->alpha); 4293 } 4294 if (kctx->threshold <= 0.0 || kctx->threshold >= 1.0) { 4295 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"0.0 < threshold (%G) < 1.0\n",kctx->threshold); 4296 } 4297 PetscFunctionReturn(0); 4298 } 4299 4300 #undef __FUNCT__ 4301 #define __FUNCT__ "SNESKSPGetParametersEW" 4302 /*@ 4303 SNESKSPGetParametersEW - Gets parameters for Eisenstat-Walker 4304 convergence criteria for the linear solvers within an inexact 4305 Newton method. 4306 4307 Not Collective 4308 4309 Input Parameters: 4310 snes - SNES context 4311 4312 Output Parameters: 4313 + version - version 1, 2 (default is 2) or 3 4314 . rtol_0 - initial relative tolerance (0 <= rtol_0 < 1) 4315 . rtol_max - maximum relative tolerance (0 <= rtol_max < 1) 4316 . gamma - multiplicative factor for version 2 rtol computation 4317 (0 <= gamma2 <= 1) 4318 . alpha - power for version 2 rtol computation (1 < alpha <= 2) 4319 . alpha2 - power for safeguard 4320 - threshold - threshold for imposing safeguard (0 < threshold < 1) 4321 4322 Level: advanced 4323 4324 .keywords: SNES, KSP, Eisenstat, Walker, get, parameters 4325 4326 .seealso: SNESKSPSetUseEW(), SNESKSPGetUseEW(), SNESKSPSetParametersEW() 4327 @*/ 4328 PetscErrorCode SNESKSPGetParametersEW(SNES snes,PetscInt *version,PetscReal *rtol_0,PetscReal *rtol_max, 4329 PetscReal *gamma,PetscReal *alpha,PetscReal *alpha2,PetscReal *threshold) 4330 { 4331 SNESKSPEW *kctx; 4332 PetscFunctionBegin; 4333 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4334 kctx = (SNESKSPEW*)snes->kspconvctx; 4335 if (!kctx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"No Eisenstat-Walker context existing"); 4336 if (version) *version = kctx->version; 4337 if (rtol_0) *rtol_0 = kctx->rtol_0; 4338 if (rtol_max) *rtol_max = kctx->rtol_max; 4339 if (gamma) *gamma = kctx->gamma; 4340 if (alpha) *alpha = kctx->alpha; 4341 if (alpha2) *alpha2 = kctx->alpha2; 4342 if (threshold) *threshold = kctx->threshold; 4343 PetscFunctionReturn(0); 4344 } 4345 4346 #undef __FUNCT__ 4347 #define __FUNCT__ "SNESKSPEW_PreSolve" 4348 static PetscErrorCode SNESKSPEW_PreSolve(SNES snes, KSP ksp, Vec b, Vec x) 4349 { 4350 PetscErrorCode ierr; 4351 SNESKSPEW *kctx = (SNESKSPEW*)snes->kspconvctx; 4352 PetscReal rtol=PETSC_DEFAULT,stol; 4353 4354 PetscFunctionBegin; 4355 if (!kctx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"No Eisenstat-Walker context exists"); 4356 if (!snes->iter) { /* first time in, so use the original user rtol */ 4357 rtol = kctx->rtol_0; 4358 } else { 4359 if (kctx->version == 1) { 4360 rtol = (snes->norm - kctx->lresid_last)/kctx->norm_last; 4361 if (rtol < 0.0) rtol = -rtol; 4362 stol = pow(kctx->rtol_last,kctx->alpha2); 4363 if (stol > kctx->threshold) rtol = PetscMax(rtol,stol); 4364 } else if (kctx->version == 2) { 4365 rtol = kctx->gamma * pow(snes->norm/kctx->norm_last,kctx->alpha); 4366 stol = kctx->gamma * pow(kctx->rtol_last,kctx->alpha); 4367 if (stol > kctx->threshold) rtol = PetscMax(rtol,stol); 4368 } else if (kctx->version == 3) {/* contributed by Luis Chacon, June 2006. */ 4369 rtol = kctx->gamma * pow(snes->norm/kctx->norm_last,kctx->alpha); 4370 /* safeguard: avoid sharp decrease of rtol */ 4371 stol = kctx->gamma*pow(kctx->rtol_last,kctx->alpha); 4372 stol = PetscMax(rtol,stol); 4373 rtol = PetscMin(kctx->rtol_0,stol); 4374 /* safeguard: avoid oversolving */ 4375 stol = kctx->gamma*(snes->ttol)/snes->norm; 4376 stol = PetscMax(rtol,stol); 4377 rtol = PetscMin(kctx->rtol_0,stol); 4378 } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Only versions 1, 2 or 3 are supported: %D",kctx->version); 4379 } 4380 /* safeguard: avoid rtol greater than one */ 4381 rtol = PetscMin(rtol,kctx->rtol_max); 4382 ierr = KSPSetTolerances(ksp,rtol,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);CHKERRQ(ierr); 4383 ierr = PetscInfo3(snes,"iter %D, Eisenstat-Walker (version %D) KSP rtol=%G\n",snes->iter,kctx->version,rtol);CHKERRQ(ierr); 4384 PetscFunctionReturn(0); 4385 } 4386 4387 #undef __FUNCT__ 4388 #define __FUNCT__ "SNESKSPEW_PostSolve" 4389 static PetscErrorCode SNESKSPEW_PostSolve(SNES snes, KSP ksp, Vec b, Vec x) 4390 { 4391 PetscErrorCode ierr; 4392 SNESKSPEW *kctx = (SNESKSPEW*)snes->kspconvctx; 4393 PCSide pcside; 4394 Vec lres; 4395 4396 PetscFunctionBegin; 4397 if (!kctx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"No Eisenstat-Walker context exists"); 4398 ierr = KSPGetTolerances(ksp,&kctx->rtol_last,0,0,0);CHKERRQ(ierr); 4399 ierr = SNESGetFunctionNorm(snes,&kctx->norm_last);CHKERRQ(ierr); 4400 if (kctx->version == 1) { 4401 ierr = KSPGetPCSide(ksp,&pcside);CHKERRQ(ierr); 4402 if (pcside == PC_RIGHT) { /* XXX Should we also test KSP_UNPRECONDITIONED_NORM ? */ 4403 /* KSP residual is true linear residual */ 4404 ierr = KSPGetResidualNorm(ksp,&kctx->lresid_last);CHKERRQ(ierr); 4405 } else { 4406 /* KSP residual is preconditioned residual */ 4407 /* compute true linear residual norm */ 4408 ierr = VecDuplicate(b,&lres);CHKERRQ(ierr); 4409 ierr = MatMult(snes->jacobian,x,lres);CHKERRQ(ierr); 4410 ierr = VecAYPX(lres,-1.0,b);CHKERRQ(ierr); 4411 ierr = VecNorm(lres,NORM_2,&kctx->lresid_last);CHKERRQ(ierr); 4412 ierr = VecDestroy(&lres);CHKERRQ(ierr); 4413 } 4414 } 4415 PetscFunctionReturn(0); 4416 } 4417 4418 #undef __FUNCT__ 4419 #define __FUNCT__ "SNES_KSPSolve" 4420 PetscErrorCode SNES_KSPSolve(SNES snes, KSP ksp, Vec b, Vec x) 4421 { 4422 PetscErrorCode ierr; 4423 4424 PetscFunctionBegin; 4425 if (snes->ksp_ewconv) { ierr = SNESKSPEW_PreSolve(snes,ksp,b,x);CHKERRQ(ierr); } 4426 ierr = KSPSolve(ksp,b,x);CHKERRQ(ierr); 4427 if (snes->ksp_ewconv) { ierr = SNESKSPEW_PostSolve(snes,ksp,b,x);CHKERRQ(ierr); } 4428 PetscFunctionReturn(0); 4429 } 4430 4431 #undef __FUNCT__ 4432 #define __FUNCT__ "SNESSetDM" 4433 /*@ 4434 SNESSetDM - Sets the DM that may be used by some preconditioners 4435 4436 Logically Collective on SNES 4437 4438 Input Parameters: 4439 + snes - the preconditioner context 4440 - dm - the dm 4441 4442 Level: intermediate 4443 4444 4445 .seealso: SNESGetDM(), KSPSetDM(), KSPGetDM() 4446 @*/ 4447 PetscErrorCode SNESSetDM(SNES snes,DM dm) 4448 { 4449 PetscErrorCode ierr; 4450 KSP ksp; 4451 DMSNES sdm; 4452 4453 PetscFunctionBegin; 4454 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4455 if (dm) {ierr = PetscObjectReference((PetscObject)dm);CHKERRQ(ierr);} 4456 if (snes->dm) { /* Move the DMSNES context over to the new DM unless the new DM already has one */ 4457 PetscContainer oldcontainer,container; 4458 ierr = PetscObjectQuery((PetscObject)snes->dm,"DMSNES",(PetscObject*)&oldcontainer);CHKERRQ(ierr); 4459 ierr = PetscObjectQuery((PetscObject)dm,"DMSNES",(PetscObject*)&container);CHKERRQ(ierr); 4460 if (oldcontainer && snes->dmAuto && !container) { 4461 ierr = DMCopyDMSNES(snes->dm,dm);CHKERRQ(ierr); 4462 ierr = DMGetDMSNES(snes->dm,&sdm);CHKERRQ(ierr); 4463 if (sdm->originaldm == snes->dm) { /* Grant write privileges to the replacement DM */ 4464 sdm->originaldm = dm; 4465 } 4466 } 4467 ierr = DMDestroy(&snes->dm);CHKERRQ(ierr); 4468 } 4469 snes->dm = dm; 4470 snes->dmAuto = PETSC_FALSE; 4471 ierr = SNESGetKSP(snes,&ksp);CHKERRQ(ierr); 4472 ierr = KSPSetDM(ksp,dm);CHKERRQ(ierr); 4473 ierr = KSPSetDMActive(ksp,PETSC_FALSE);CHKERRQ(ierr); 4474 if (snes->pc) { 4475 ierr = SNESSetDM(snes->pc, snes->dm);CHKERRQ(ierr); 4476 ierr = SNESSetPCSide(snes,snes->pcside);CHKERRQ(ierr); 4477 } 4478 PetscFunctionReturn(0); 4479 } 4480 4481 #undef __FUNCT__ 4482 #define __FUNCT__ "SNESGetDM" 4483 /*@ 4484 SNESGetDM - Gets the DM that may be used by some preconditioners 4485 4486 Not Collective but DM obtained is parallel on SNES 4487 4488 Input Parameter: 4489 . snes - the preconditioner context 4490 4491 Output Parameter: 4492 . dm - the dm 4493 4494 Level: intermediate 4495 4496 4497 .seealso: SNESSetDM(), KSPSetDM(), KSPGetDM() 4498 @*/ 4499 PetscErrorCode SNESGetDM(SNES snes,DM *dm) 4500 { 4501 PetscErrorCode ierr; 4502 4503 PetscFunctionBegin; 4504 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4505 if (!snes->dm) { 4506 ierr = DMShellCreate(((PetscObject)snes)->comm,&snes->dm);CHKERRQ(ierr); 4507 snes->dmAuto = PETSC_TRUE; 4508 } 4509 *dm = snes->dm; 4510 PetscFunctionReturn(0); 4511 } 4512 4513 #undef __FUNCT__ 4514 #define __FUNCT__ "SNESSetPC" 4515 /*@ 4516 SNESSetPC - Sets the nonlinear preconditioner to be used. 4517 4518 Collective on SNES 4519 4520 Input Parameters: 4521 + snes - iterative context obtained from SNESCreate() 4522 - pc - the preconditioner object 4523 4524 Notes: 4525 Use SNESGetPC() to retrieve the preconditioner context (for example, 4526 to configure it using the API). 4527 4528 Level: developer 4529 4530 .keywords: SNES, set, precondition 4531 .seealso: SNESGetPC() 4532 @*/ 4533 PetscErrorCode SNESSetPC(SNES snes, SNES pc) 4534 { 4535 PetscErrorCode ierr; 4536 4537 PetscFunctionBegin; 4538 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); 4539 PetscValidHeaderSpecific(pc, SNES_CLASSID, 2); 4540 PetscCheckSameComm(snes, 1, pc, 2); 4541 ierr = PetscObjectReference((PetscObject) pc);CHKERRQ(ierr); 4542 ierr = SNESDestroy(&snes->pc);CHKERRQ(ierr); 4543 snes->pc = pc; 4544 ierr = PetscLogObjectParent(snes, snes->pc);CHKERRQ(ierr); 4545 PetscFunctionReturn(0); 4546 } 4547 4548 #undef __FUNCT__ 4549 #define __FUNCT__ "SNESGetPC" 4550 /*@ 4551 SNESGetPC - Returns a pointer to the nonlinear preconditioning context set with SNESSetPC(). 4552 4553 Not Collective 4554 4555 Input Parameter: 4556 . snes - iterative context obtained from SNESCreate() 4557 4558 Output Parameter: 4559 . pc - preconditioner context 4560 4561 Level: developer 4562 4563 .keywords: SNES, get, preconditioner 4564 .seealso: SNESSetPC() 4565 @*/ 4566 PetscErrorCode SNESGetPC(SNES snes, SNES *pc) 4567 { 4568 PetscErrorCode ierr; 4569 const char *optionsprefix; 4570 4571 PetscFunctionBegin; 4572 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); 4573 PetscValidPointer(pc, 2); 4574 if (!snes->pc) { 4575 ierr = SNESCreate(((PetscObject) snes)->comm,&snes->pc);CHKERRQ(ierr); 4576 ierr = PetscObjectIncrementTabLevel((PetscObject)snes->pc,(PetscObject)snes,1);CHKERRQ(ierr); 4577 ierr = PetscLogObjectParent(snes,snes->pc);CHKERRQ(ierr); 4578 ierr = SNESGetOptionsPrefix(snes,&optionsprefix);CHKERRQ(ierr); 4579 ierr = SNESSetOptionsPrefix(snes->pc,optionsprefix);CHKERRQ(ierr); 4580 ierr = SNESAppendOptionsPrefix(snes->pc,"npc_");CHKERRQ(ierr); 4581 } 4582 *pc = snes->pc; 4583 PetscFunctionReturn(0); 4584 } 4585 4586 4587 #undef __FUNCT__ 4588 #define __FUNCT__ "SNESSetPCSide" 4589 /*@ 4590 SNESSetPCSide - Sets the preconditioning side. 4591 4592 Logically Collective on SNES 4593 4594 Input Parameter: 4595 . snes - iterative context obtained from SNESCreate() 4596 4597 Output Parameter: 4598 . side - the preconditioning side, where side is one of 4599 .vb 4600 PC_LEFT - left preconditioning (default) 4601 PC_RIGHT - right preconditioning 4602 .ve 4603 4604 Options Database Keys: 4605 . -snes_pc_side <right,left> 4606 4607 Level: intermediate 4608 4609 .keywords: SNES, set, right, left, side, preconditioner, flag 4610 4611 .seealso: SNESGetPCSide(), KSPSetPCSide() 4612 @*/ 4613 PetscErrorCode SNESSetPCSide(SNES snes,PCSide side) 4614 { 4615 PetscFunctionBegin; 4616 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4617 PetscValidLogicalCollectiveEnum(snes,side,2); 4618 snes->pcside = side; 4619 PetscFunctionReturn(0); 4620 } 4621 4622 #undef __FUNCT__ 4623 #define __FUNCT__ "SNESGetPCSide" 4624 /*@ 4625 SNESGetPCSide - Gets the preconditioning side. 4626 4627 Not Collective 4628 4629 Input Parameter: 4630 . snes - iterative context obtained from SNESCreate() 4631 4632 Output Parameter: 4633 . side - the preconditioning side, where side is one of 4634 .vb 4635 PC_LEFT - left preconditioning (default) 4636 PC_RIGHT - right preconditioning 4637 .ve 4638 4639 Level: intermediate 4640 4641 .keywords: SNES, get, right, left, side, preconditioner, flag 4642 4643 .seealso: SNESSetPCSide(), KSPGetPCSide() 4644 @*/ 4645 PetscErrorCode SNESGetPCSide(SNES snes,PCSide *side) 4646 { 4647 PetscFunctionBegin; 4648 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4649 PetscValidPointer(side,2); 4650 *side = snes->pcside; 4651 PetscFunctionReturn(0); 4652 } 4653 4654 #undef __FUNCT__ 4655 #define __FUNCT__ "SNESSetSNESLineSearch" 4656 /*@ 4657 SNESSetSNESLineSearch - Sets the linesearch on the SNES instance. 4658 4659 Collective on SNES 4660 4661 Input Parameters: 4662 + snes - iterative context obtained from SNESCreate() 4663 - linesearch - the linesearch object 4664 4665 Notes: 4666 Use SNESGetSNESLineSearch() to retrieve the preconditioner context (for example, 4667 to configure it using the API). 4668 4669 Level: developer 4670 4671 .keywords: SNES, set, linesearch 4672 .seealso: SNESGetSNESLineSearch() 4673 @*/ 4674 PetscErrorCode SNESSetSNESLineSearch(SNES snes, SNESLineSearch linesearch) 4675 { 4676 PetscErrorCode ierr; 4677 4678 PetscFunctionBegin; 4679 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); 4680 PetscValidHeaderSpecific(linesearch, SNESLINESEARCH_CLASSID, 2); 4681 PetscCheckSameComm(snes, 1, linesearch, 2); 4682 ierr = PetscObjectReference((PetscObject) linesearch);CHKERRQ(ierr); 4683 ierr = SNESLineSearchDestroy(&snes->linesearch);CHKERRQ(ierr); 4684 snes->linesearch = linesearch; 4685 ierr = PetscLogObjectParent(snes, snes->linesearch);CHKERRQ(ierr); 4686 PetscFunctionReturn(0); 4687 } 4688 4689 #undef __FUNCT__ 4690 #define __FUNCT__ "SNESGetSNESLineSearch" 4691 /*@C 4692 SNESGetSNESLineSearch - Returns a pointer to the line search context set with SNESSetLineSearch() 4693 or creates a default line search instance associated with the SNES and returns it. 4694 4695 Not Collective 4696 4697 Input Parameter: 4698 . snes - iterative context obtained from SNESCreate() 4699 4700 Output Parameter: 4701 . linesearch - linesearch context 4702 4703 Level: developer 4704 4705 .keywords: SNES, get, linesearch 4706 .seealso: SNESSetSNESLineSearch() 4707 @*/ 4708 PetscErrorCode SNESGetSNESLineSearch(SNES snes, SNESLineSearch *linesearch) 4709 { 4710 PetscErrorCode ierr; 4711 const char *optionsprefix; 4712 4713 PetscFunctionBegin; 4714 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); 4715 PetscValidPointer(linesearch, 2); 4716 if (!snes->linesearch) { 4717 ierr = SNESGetOptionsPrefix(snes, &optionsprefix);CHKERRQ(ierr); 4718 ierr = SNESLineSearchCreate(((PetscObject) snes)->comm, &snes->linesearch);CHKERRQ(ierr); 4719 ierr = SNESLineSearchSetSNES(snes->linesearch, snes);CHKERRQ(ierr); 4720 ierr = SNESLineSearchAppendOptionsPrefix(snes->linesearch, optionsprefix);CHKERRQ(ierr); 4721 ierr = PetscObjectIncrementTabLevel((PetscObject) snes->linesearch, (PetscObject) snes, 1);CHKERRQ(ierr); 4722 ierr = PetscLogObjectParent(snes, snes->linesearch);CHKERRQ(ierr); 4723 } 4724 *linesearch = snes->linesearch; 4725 PetscFunctionReturn(0); 4726 } 4727 4728 #if defined(PETSC_HAVE_MATLAB_ENGINE) 4729 #include <mex.h> 4730 4731 typedef struct {char *funcname; mxArray *ctx;} SNESMatlabContext; 4732 4733 #undef __FUNCT__ 4734 #define __FUNCT__ "SNESComputeFunction_Matlab" 4735 /* 4736 SNESComputeFunction_Matlab - Calls the function that has been set with 4737 SNESSetFunctionMatlab(). 4738 4739 Collective on SNES 4740 4741 Input Parameters: 4742 + snes - the SNES context 4743 - x - input vector 4744 4745 Output Parameter: 4746 . y - function vector, as set by SNESSetFunction() 4747 4748 Notes: 4749 SNESComputeFunction() is typically used within nonlinear solvers 4750 implementations, so most users would not generally call this routine 4751 themselves. 4752 4753 Level: developer 4754 4755 .keywords: SNES, nonlinear, compute, function 4756 4757 .seealso: SNESSetFunction(), SNESGetFunction() 4758 */ 4759 PetscErrorCode SNESComputeFunction_Matlab(SNES snes,Vec x,Vec y, void *ctx) 4760 { 4761 PetscErrorCode ierr; 4762 SNESMatlabContext *sctx = (SNESMatlabContext *)ctx; 4763 int nlhs = 1,nrhs = 5; 4764 mxArray *plhs[1],*prhs[5]; 4765 long long int lx = 0,ly = 0,ls = 0; 4766 4767 PetscFunctionBegin; 4768 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4769 PetscValidHeaderSpecific(x,VEC_CLASSID,2); 4770 PetscValidHeaderSpecific(y,VEC_CLASSID,3); 4771 PetscCheckSameComm(snes,1,x,2); 4772 PetscCheckSameComm(snes,1,y,3); 4773 4774 /* call Matlab function in ctx with arguments x and y */ 4775 4776 ierr = PetscMemcpy(&ls,&snes,sizeof(snes));CHKERRQ(ierr); 4777 ierr = PetscMemcpy(&lx,&x,sizeof(x));CHKERRQ(ierr); 4778 ierr = PetscMemcpy(&ly,&y,sizeof(x));CHKERRQ(ierr); 4779 prhs[0] = mxCreateDoubleScalar((double)ls); 4780 prhs[1] = mxCreateDoubleScalar((double)lx); 4781 prhs[2] = mxCreateDoubleScalar((double)ly); 4782 prhs[3] = mxCreateString(sctx->funcname); 4783 prhs[4] = sctx->ctx; 4784 ierr = mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscSNESComputeFunctionInternal");CHKERRQ(ierr); 4785 ierr = mxGetScalar(plhs[0]);CHKERRQ(ierr); 4786 mxDestroyArray(prhs[0]); 4787 mxDestroyArray(prhs[1]); 4788 mxDestroyArray(prhs[2]); 4789 mxDestroyArray(prhs[3]); 4790 mxDestroyArray(plhs[0]); 4791 PetscFunctionReturn(0); 4792 } 4793 4794 4795 #undef __FUNCT__ 4796 #define __FUNCT__ "SNESSetFunctionMatlab" 4797 /* 4798 SNESSetFunctionMatlab - Sets the function evaluation routine and function 4799 vector for use by the SNES routines in solving systems of nonlinear 4800 equations from MATLAB. Here the function is a string containing the name of a MATLAB function 4801 4802 Logically Collective on SNES 4803 4804 Input Parameters: 4805 + snes - the SNES context 4806 . r - vector to store function value 4807 - func - function evaluation routine 4808 4809 Calling sequence of func: 4810 $ func (SNES snes,Vec x,Vec f,void *ctx); 4811 4812 4813 Notes: 4814 The Newton-like methods typically solve linear systems of the form 4815 $ f'(x) x = -f(x), 4816 where f'(x) denotes the Jacobian matrix and f(x) is the function. 4817 4818 Level: beginner 4819 4820 .keywords: SNES, nonlinear, set, function 4821 4822 .seealso: SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESSetFunction() 4823 */ 4824 PetscErrorCode SNESSetFunctionMatlab(SNES snes,Vec r,const char *func,mxArray *ctx) 4825 { 4826 PetscErrorCode ierr; 4827 SNESMatlabContext *sctx; 4828 4829 PetscFunctionBegin; 4830 /* currently sctx is memory bleed */ 4831 ierr = PetscMalloc(sizeof(SNESMatlabContext),&sctx);CHKERRQ(ierr); 4832 ierr = PetscStrallocpy(func,&sctx->funcname);CHKERRQ(ierr); 4833 /* 4834 This should work, but it doesn't 4835 sctx->ctx = ctx; 4836 mexMakeArrayPersistent(sctx->ctx); 4837 */ 4838 sctx->ctx = mxDuplicateArray(ctx); 4839 ierr = SNESSetFunction(snes,r,SNESComputeFunction_Matlab,sctx);CHKERRQ(ierr); 4840 PetscFunctionReturn(0); 4841 } 4842 4843 #undef __FUNCT__ 4844 #define __FUNCT__ "SNESComputeJacobian_Matlab" 4845 /* 4846 SNESComputeJacobian_Matlab - Calls the function that has been set with 4847 SNESSetJacobianMatlab(). 4848 4849 Collective on SNES 4850 4851 Input Parameters: 4852 + snes - the SNES context 4853 . x - input vector 4854 . A, B - the matrices 4855 - ctx - user context 4856 4857 Output Parameter: 4858 . flag - structure of the matrix 4859 4860 Level: developer 4861 4862 .keywords: SNES, nonlinear, compute, function 4863 4864 .seealso: SNESSetFunction(), SNESGetFunction() 4865 @*/ 4866 PetscErrorCode SNESComputeJacobian_Matlab(SNES snes,Vec x,Mat *A,Mat *B,MatStructure *flag, void *ctx) 4867 { 4868 PetscErrorCode ierr; 4869 SNESMatlabContext *sctx = (SNESMatlabContext *)ctx; 4870 int nlhs = 2,nrhs = 6; 4871 mxArray *plhs[2],*prhs[6]; 4872 long long int lx = 0,lA = 0,ls = 0, lB = 0; 4873 4874 PetscFunctionBegin; 4875 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4876 PetscValidHeaderSpecific(x,VEC_CLASSID,2); 4877 4878 /* call Matlab function in ctx with arguments x and y */ 4879 4880 ierr = PetscMemcpy(&ls,&snes,sizeof(snes));CHKERRQ(ierr); 4881 ierr = PetscMemcpy(&lx,&x,sizeof(x));CHKERRQ(ierr); 4882 ierr = PetscMemcpy(&lA,A,sizeof(x));CHKERRQ(ierr); 4883 ierr = PetscMemcpy(&lB,B,sizeof(x));CHKERRQ(ierr); 4884 prhs[0] = mxCreateDoubleScalar((double)ls); 4885 prhs[1] = mxCreateDoubleScalar((double)lx); 4886 prhs[2] = mxCreateDoubleScalar((double)lA); 4887 prhs[3] = mxCreateDoubleScalar((double)lB); 4888 prhs[4] = mxCreateString(sctx->funcname); 4889 prhs[5] = sctx->ctx; 4890 ierr = mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscSNESComputeJacobianInternal");CHKERRQ(ierr); 4891 ierr = mxGetScalar(plhs[0]);CHKERRQ(ierr); 4892 *flag = (MatStructure) mxGetScalar(plhs[1]);CHKERRQ(ierr); 4893 mxDestroyArray(prhs[0]); 4894 mxDestroyArray(prhs[1]); 4895 mxDestroyArray(prhs[2]); 4896 mxDestroyArray(prhs[3]); 4897 mxDestroyArray(prhs[4]); 4898 mxDestroyArray(plhs[0]); 4899 mxDestroyArray(plhs[1]); 4900 PetscFunctionReturn(0); 4901 } 4902 4903 4904 #undef __FUNCT__ 4905 #define __FUNCT__ "SNESSetJacobianMatlab" 4906 /* 4907 SNESSetJacobianMatlab - Sets the Jacobian function evaluation routine and two empty Jacobian matrices 4908 vector for use by the SNES routines in solving systems of nonlinear 4909 equations from MATLAB. Here the function is a string containing the name of a MATLAB function 4910 4911 Logically Collective on SNES 4912 4913 Input Parameters: 4914 + snes - the SNES context 4915 . A,B - Jacobian matrices 4916 . func - function evaluation routine 4917 - ctx - user context 4918 4919 Calling sequence of func: 4920 $ flag = func (SNES snes,Vec x,Mat A,Mat B,void *ctx); 4921 4922 4923 Level: developer 4924 4925 .keywords: SNES, nonlinear, set, function 4926 4927 .seealso: SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESSetFunction() 4928 */ 4929 PetscErrorCode SNESSetJacobianMatlab(SNES snes,Mat A,Mat B,const char *func,mxArray *ctx) 4930 { 4931 PetscErrorCode ierr; 4932 SNESMatlabContext *sctx; 4933 4934 PetscFunctionBegin; 4935 /* currently sctx is memory bleed */ 4936 ierr = PetscMalloc(sizeof(SNESMatlabContext),&sctx);CHKERRQ(ierr); 4937 ierr = PetscStrallocpy(func,&sctx->funcname);CHKERRQ(ierr); 4938 /* 4939 This should work, but it doesn't 4940 sctx->ctx = ctx; 4941 mexMakeArrayPersistent(sctx->ctx); 4942 */ 4943 sctx->ctx = mxDuplicateArray(ctx); 4944 ierr = SNESSetJacobian(snes,A,B,SNESComputeJacobian_Matlab,sctx);CHKERRQ(ierr); 4945 PetscFunctionReturn(0); 4946 } 4947 4948 #undef __FUNCT__ 4949 #define __FUNCT__ "SNESMonitor_Matlab" 4950 /* 4951 SNESMonitor_Matlab - Calls the function that has been set with SNESMonitorSetMatlab(). 4952 4953 Collective on SNES 4954 4955 .seealso: SNESSetFunction(), SNESGetFunction() 4956 @*/ 4957 PetscErrorCode SNESMonitor_Matlab(SNES snes,PetscInt it, PetscReal fnorm, void *ctx) 4958 { 4959 PetscErrorCode ierr; 4960 SNESMatlabContext *sctx = (SNESMatlabContext *)ctx; 4961 int nlhs = 1,nrhs = 6; 4962 mxArray *plhs[1],*prhs[6]; 4963 long long int lx = 0,ls = 0; 4964 Vec x=snes->vec_sol; 4965 4966 PetscFunctionBegin; 4967 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4968 4969 ierr = PetscMemcpy(&ls,&snes,sizeof(snes));CHKERRQ(ierr); 4970 ierr = PetscMemcpy(&lx,&x,sizeof(x));CHKERRQ(ierr); 4971 prhs[0] = mxCreateDoubleScalar((double)ls); 4972 prhs[1] = mxCreateDoubleScalar((double)it); 4973 prhs[2] = mxCreateDoubleScalar((double)fnorm); 4974 prhs[3] = mxCreateDoubleScalar((double)lx); 4975 prhs[4] = mxCreateString(sctx->funcname); 4976 prhs[5] = sctx->ctx; 4977 ierr = mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscSNESMonitorInternal");CHKERRQ(ierr); 4978 ierr = mxGetScalar(plhs[0]);CHKERRQ(ierr); 4979 mxDestroyArray(prhs[0]); 4980 mxDestroyArray(prhs[1]); 4981 mxDestroyArray(prhs[2]); 4982 mxDestroyArray(prhs[3]); 4983 mxDestroyArray(prhs[4]); 4984 mxDestroyArray(plhs[0]); 4985 PetscFunctionReturn(0); 4986 } 4987 4988 4989 #undef __FUNCT__ 4990 #define __FUNCT__ "SNESMonitorSetMatlab" 4991 /* 4992 SNESMonitorSetMatlab - Sets the monitor function from MATLAB 4993 4994 Level: developer 4995 4996 .keywords: SNES, nonlinear, set, function 4997 4998 .seealso: SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESSetFunction() 4999 */ 5000 PetscErrorCode SNESMonitorSetMatlab(SNES snes,const char *func,mxArray *ctx) 5001 { 5002 PetscErrorCode ierr; 5003 SNESMatlabContext *sctx; 5004 5005 PetscFunctionBegin; 5006 /* currently sctx is memory bleed */ 5007 ierr = PetscMalloc(sizeof(SNESMatlabContext),&sctx);CHKERRQ(ierr); 5008 ierr = PetscStrallocpy(func,&sctx->funcname);CHKERRQ(ierr); 5009 /* 5010 This should work, but it doesn't 5011 sctx->ctx = ctx; 5012 mexMakeArrayPersistent(sctx->ctx); 5013 */ 5014 sctx->ctx = mxDuplicateArray(ctx); 5015 ierr = SNESMonitorSet(snes,SNESMonitor_Matlab,sctx,PETSC_NULL);CHKERRQ(ierr); 5016 PetscFunctionReturn(0); 5017 } 5018 5019 #endif 5020