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