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