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