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 } else if (snes->dm) { 3638 PetscBool ig; 3639 ierr = DMHasInitialGuess(snes->dm,&ig);CHKERRQ(ierr); 3640 if (ig) { 3641 ierr = DMComputeInitialGuess(snes->dm,snes->vec_sol);CHKERRQ(ierr); 3642 } 3643 } 3644 } 3645 3646 if (snes->conv_hist_reset) snes->conv_hist_len = 0; 3647 snes->nfuncs = 0; snes->linear_its = 0; snes->numFailures = 0; 3648 3649 ierr = PetscLogEventBegin(SNES_Solve,snes,0,0,0);CHKERRQ(ierr); 3650 ierr = (*snes->ops->solve)(snes);CHKERRQ(ierr); 3651 ierr = PetscLogEventEnd(SNES_Solve,snes,0,0,0);CHKERRQ(ierr); 3652 if (snes->domainerror){ 3653 snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; 3654 snes->domainerror = PETSC_FALSE; 3655 } 3656 if (!snes->reason) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Internal error, solver returned without setting converged reason"); 3657 3658 flg = PETSC_FALSE; 3659 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_test_local_min",&flg,PETSC_NULL);CHKERRQ(ierr); 3660 if (flg && !PetscPreLoadingOn) { ierr = SNESTestLocalMin(snes);CHKERRQ(ierr); } 3661 if (snes->printreason) { 3662 ierr = PetscViewerASCIIAddTab(PETSC_VIEWER_STDOUT_(((PetscObject)snes)->comm),((PetscObject)snes)->tablevel);CHKERRQ(ierr); 3663 if (snes->reason > 0) { 3664 ierr = PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_(((PetscObject)snes)->comm),"Nonlinear solve converged due to %s iterations %D\n",SNESConvergedReasons[snes->reason],snes->iter);CHKERRQ(ierr); 3665 } else { 3666 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); 3667 } 3668 ierr = PetscViewerASCIISubtractTab(PETSC_VIEWER_STDOUT_(((PetscObject)snes)->comm),((PetscObject)snes)->tablevel);CHKERRQ(ierr); 3669 } 3670 3671 if (snes->errorifnotconverged && snes->reason < 0) SETERRQ(((PetscObject)snes)->comm,PETSC_ERR_NOT_CONVERGED,"SNESSolve has not converged"); 3672 if (grid < snes->gridsequence) { 3673 DM fine; 3674 Vec xnew; 3675 Mat interp; 3676 3677 ierr = DMRefine(snes->dm,((PetscObject)snes)->comm,&fine);CHKERRQ(ierr); 3678 if (!fine) SETERRQ(((PetscObject)snes)->comm,PETSC_ERR_ARG_INCOMP,"DMRefine() did not perform any refinement, cannot continue grid sequencing"); 3679 ierr = DMCreateInterpolation(snes->dm,fine,&interp,PETSC_NULL);CHKERRQ(ierr); 3680 ierr = DMCreateGlobalVector(fine,&xnew);CHKERRQ(ierr); 3681 ierr = MatInterpolate(interp,x,xnew);CHKERRQ(ierr); 3682 ierr = DMInterpolate(snes->dm,interp,fine);CHKERRQ(ierr); 3683 ierr = MatDestroy(&interp);CHKERRQ(ierr); 3684 x = xnew; 3685 3686 ierr = SNESReset(snes);CHKERRQ(ierr); 3687 ierr = SNESSetDM(snes,fine);CHKERRQ(ierr); 3688 ierr = DMDestroy(&fine);CHKERRQ(ierr); 3689 ierr = PetscViewerASCIIPopTab(PETSC_VIEWER_STDOUT_(((PetscObject)snes)->comm));CHKERRQ(ierr); 3690 } 3691 } 3692 /* monitoring and viewing */ 3693 flg = PETSC_FALSE; 3694 ierr = PetscOptionsGetString(((PetscObject)snes)->prefix,"-snes_view",filename,PETSC_MAX_PATH_LEN,&flg);CHKERRQ(ierr); 3695 if (flg && !PetscPreLoadingOn) { 3696 ierr = PetscViewerASCIIOpen(((PetscObject)snes)->comm,filename,&viewer);CHKERRQ(ierr); 3697 ierr = SNESView(snes,viewer);CHKERRQ(ierr); 3698 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 3699 } 3700 3701 flg = PETSC_FALSE; 3702 ierr = PetscOptionsGetString(((PetscObject)snes)->prefix,"-snes_view_solution_vtk",filename,PETSC_MAX_PATH_LEN,&flg);CHKERRQ(ierr); 3703 if (flg) { 3704 PetscViewer viewer; 3705 ierr = PetscViewerCreate(((PetscObject)snes)->comm,&viewer);CHKERRQ(ierr); 3706 ierr = PetscViewerSetType(viewer,PETSCVIEWERVTK);CHKERRQ(ierr); 3707 ierr = PetscViewerFileSetName(viewer,filename);CHKERRQ(ierr); 3708 ierr = VecView(snes->vec_sol,viewer);CHKERRQ(ierr); 3709 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 3710 } 3711 3712 ierr = VecDestroy(&xcreated);CHKERRQ(ierr); 3713 PetscFunctionReturn(0); 3714 } 3715 3716 /* --------- Internal routines for SNES Package --------- */ 3717 3718 #undef __FUNCT__ 3719 #define __FUNCT__ "SNESSetType" 3720 /*@C 3721 SNESSetType - Sets the method for the nonlinear solver. 3722 3723 Collective on SNES 3724 3725 Input Parameters: 3726 + snes - the SNES context 3727 - type - a known method 3728 3729 Options Database Key: 3730 . -snes_type <type> - Sets the method; use -help for a list 3731 of available methods (for instance, ls or tr) 3732 3733 Notes: 3734 See "petsc/include/petscsnes.h" for available methods (for instance) 3735 + SNESLS - Newton's method with line search 3736 (systems of nonlinear equations) 3737 . SNESTR - Newton's method with trust region 3738 (systems of nonlinear equations) 3739 3740 Normally, it is best to use the SNESSetFromOptions() command and then 3741 set the SNES solver type from the options database rather than by using 3742 this routine. Using the options database provides the user with 3743 maximum flexibility in evaluating the many nonlinear solvers. 3744 The SNESSetType() routine is provided for those situations where it 3745 is necessary to set the nonlinear solver independently of the command 3746 line or options database. This might be the case, for example, when 3747 the choice of solver changes during the execution of the program, 3748 and the user's application is taking responsibility for choosing the 3749 appropriate method. 3750 3751 Level: intermediate 3752 3753 .keywords: SNES, set, type 3754 3755 .seealso: SNESType, SNESCreate() 3756 3757 @*/ 3758 PetscErrorCode SNESSetType(SNES snes,SNESType type) 3759 { 3760 PetscErrorCode ierr,(*r)(SNES); 3761 PetscBool match; 3762 3763 PetscFunctionBegin; 3764 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3765 PetscValidCharPointer(type,2); 3766 3767 ierr = PetscObjectTypeCompare((PetscObject)snes,type,&match);CHKERRQ(ierr); 3768 if (match) PetscFunctionReturn(0); 3769 3770 ierr = PetscFListFind(SNESList,((PetscObject)snes)->comm,type,PETSC_TRUE,(void (**)(void)) &r);CHKERRQ(ierr); 3771 if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unable to find requested SNES type %s",type); 3772 /* Destroy the previous private SNES context */ 3773 if (snes->ops->destroy) { 3774 ierr = (*(snes)->ops->destroy)(snes);CHKERRQ(ierr); 3775 snes->ops->destroy = PETSC_NULL; 3776 } 3777 /* Reinitialize function pointers in SNESOps structure */ 3778 snes->ops->setup = 0; 3779 snes->ops->solve = 0; 3780 snes->ops->view = 0; 3781 snes->ops->setfromoptions = 0; 3782 snes->ops->destroy = 0; 3783 /* Call the SNESCreate_XXX routine for this particular Nonlinear solver */ 3784 snes->setupcalled = PETSC_FALSE; 3785 ierr = PetscObjectChangeTypeName((PetscObject)snes,type);CHKERRQ(ierr); 3786 ierr = (*r)(snes);CHKERRQ(ierr); 3787 #if defined(PETSC_HAVE_AMS) 3788 if (PetscAMSPublishAll) { 3789 ierr = PetscObjectAMSPublish((PetscObject)snes);CHKERRQ(ierr); 3790 } 3791 #endif 3792 PetscFunctionReturn(0); 3793 } 3794 3795 3796 /* --------------------------------------------------------------------- */ 3797 #undef __FUNCT__ 3798 #define __FUNCT__ "SNESRegisterDestroy" 3799 /*@ 3800 SNESRegisterDestroy - Frees the list of nonlinear solvers that were 3801 registered by SNESRegisterDynamic(). 3802 3803 Not Collective 3804 3805 Level: advanced 3806 3807 .keywords: SNES, nonlinear, register, destroy 3808 3809 .seealso: SNESRegisterAll(), SNESRegisterAll() 3810 @*/ 3811 PetscErrorCode SNESRegisterDestroy(void) 3812 { 3813 PetscErrorCode ierr; 3814 3815 PetscFunctionBegin; 3816 ierr = PetscFListDestroy(&SNESList);CHKERRQ(ierr); 3817 SNESRegisterAllCalled = PETSC_FALSE; 3818 PetscFunctionReturn(0); 3819 } 3820 3821 #undef __FUNCT__ 3822 #define __FUNCT__ "SNESGetType" 3823 /*@C 3824 SNESGetType - Gets the SNES method type and name (as a string). 3825 3826 Not Collective 3827 3828 Input Parameter: 3829 . snes - nonlinear solver context 3830 3831 Output Parameter: 3832 . type - SNES method (a character string) 3833 3834 Level: intermediate 3835 3836 .keywords: SNES, nonlinear, get, type, name 3837 @*/ 3838 PetscErrorCode SNESGetType(SNES snes,SNESType *type) 3839 { 3840 PetscFunctionBegin; 3841 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3842 PetscValidPointer(type,2); 3843 *type = ((PetscObject)snes)->type_name; 3844 PetscFunctionReturn(0); 3845 } 3846 3847 #undef __FUNCT__ 3848 #define __FUNCT__ "SNESGetSolution" 3849 /*@ 3850 SNESGetSolution - Returns the vector where the approximate solution is 3851 stored. This is the fine grid solution when using SNESSetGridSequence(). 3852 3853 Not Collective, but Vec is parallel if SNES is parallel 3854 3855 Input Parameter: 3856 . snes - the SNES context 3857 3858 Output Parameter: 3859 . x - the solution 3860 3861 Level: intermediate 3862 3863 .keywords: SNES, nonlinear, get, solution 3864 3865 .seealso: SNESGetSolutionUpdate(), SNESGetFunction() 3866 @*/ 3867 PetscErrorCode SNESGetSolution(SNES snes,Vec *x) 3868 { 3869 PetscFunctionBegin; 3870 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3871 PetscValidPointer(x,2); 3872 *x = snes->vec_sol; 3873 PetscFunctionReturn(0); 3874 } 3875 3876 #undef __FUNCT__ 3877 #define __FUNCT__ "SNESGetSolutionUpdate" 3878 /*@ 3879 SNESGetSolutionUpdate - Returns the vector where the solution update is 3880 stored. 3881 3882 Not Collective, but Vec is parallel if SNES is parallel 3883 3884 Input Parameter: 3885 . snes - the SNES context 3886 3887 Output Parameter: 3888 . x - the solution update 3889 3890 Level: advanced 3891 3892 .keywords: SNES, nonlinear, get, solution, update 3893 3894 .seealso: SNESGetSolution(), SNESGetFunction() 3895 @*/ 3896 PetscErrorCode SNESGetSolutionUpdate(SNES snes,Vec *x) 3897 { 3898 PetscFunctionBegin; 3899 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3900 PetscValidPointer(x,2); 3901 *x = snes->vec_sol_update; 3902 PetscFunctionReturn(0); 3903 } 3904 3905 #undef __FUNCT__ 3906 #define __FUNCT__ "SNESGetFunction" 3907 /*@C 3908 SNESGetFunction - Returns the vector where the function is stored. 3909 3910 Not Collective, but Vec is parallel if SNES is parallel. Collective if Vec is requested, but has not been created yet. 3911 3912 Input Parameter: 3913 . snes - the SNES context 3914 3915 Output Parameter: 3916 + r - the function (or PETSC_NULL) 3917 . func - the function (or PETSC_NULL) 3918 - ctx - the function context (or PETSC_NULL) 3919 3920 Level: advanced 3921 3922 .keywords: SNES, nonlinear, get, function 3923 3924 .seealso: SNESSetFunction(), SNESGetSolution() 3925 @*/ 3926 PetscErrorCode SNESGetFunction(SNES snes,Vec *r,PetscErrorCode (**func)(SNES,Vec,Vec,void*),void **ctx) 3927 { 3928 PetscErrorCode ierr; 3929 DM dm; 3930 3931 PetscFunctionBegin; 3932 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3933 if (r) { 3934 if (!snes->vec_func) { 3935 if (snes->vec_rhs) { 3936 ierr = VecDuplicate(snes->vec_rhs,&snes->vec_func);CHKERRQ(ierr); 3937 } else if (snes->vec_sol) { 3938 ierr = VecDuplicate(snes->vec_sol,&snes->vec_func);CHKERRQ(ierr); 3939 } else if (snes->dm) { 3940 ierr = DMCreateGlobalVector(snes->dm,&snes->vec_func);CHKERRQ(ierr); 3941 } 3942 } 3943 *r = snes->vec_func; 3944 } 3945 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 3946 ierr = DMSNESGetFunction(dm,func,ctx);CHKERRQ(ierr); 3947 PetscFunctionReturn(0); 3948 } 3949 3950 /*@C 3951 SNESGetGS - Returns the GS function and context. 3952 3953 Input Parameter: 3954 . snes - the SNES context 3955 3956 Output Parameter: 3957 + gsfunc - the function (or PETSC_NULL) 3958 - ctx - the function context (or PETSC_NULL) 3959 3960 Level: advanced 3961 3962 .keywords: SNES, nonlinear, get, function 3963 3964 .seealso: SNESSetGS(), SNESGetFunction() 3965 @*/ 3966 3967 #undef __FUNCT__ 3968 #define __FUNCT__ "SNESGetGS" 3969 PetscErrorCode SNESGetGS (SNES snes, PetscErrorCode(**func)(SNES, Vec, Vec, void*), void ** ctx) 3970 { 3971 PetscErrorCode ierr; 3972 DM dm; 3973 3974 PetscFunctionBegin; 3975 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3976 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 3977 ierr = DMSNESGetGS(dm,func,ctx);CHKERRQ(ierr); 3978 PetscFunctionReturn(0); 3979 } 3980 3981 #undef __FUNCT__ 3982 #define __FUNCT__ "SNESSetOptionsPrefix" 3983 /*@C 3984 SNESSetOptionsPrefix - Sets the prefix used for searching for all 3985 SNES options in the database. 3986 3987 Logically Collective on SNES 3988 3989 Input Parameter: 3990 + snes - the SNES context 3991 - prefix - the prefix to prepend to all option names 3992 3993 Notes: 3994 A hyphen (-) must NOT be given at the beginning of the prefix name. 3995 The first character of all runtime options is AUTOMATICALLY the hyphen. 3996 3997 Level: advanced 3998 3999 .keywords: SNES, set, options, prefix, database 4000 4001 .seealso: SNESSetFromOptions() 4002 @*/ 4003 PetscErrorCode SNESSetOptionsPrefix(SNES snes,const char prefix[]) 4004 { 4005 PetscErrorCode ierr; 4006 4007 PetscFunctionBegin; 4008 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4009 ierr = PetscObjectSetOptionsPrefix((PetscObject)snes,prefix);CHKERRQ(ierr); 4010 if (!snes->ksp) {ierr = SNESGetKSP(snes,&snes->ksp);CHKERRQ(ierr);} 4011 if (snes->linesearch) { 4012 ierr = SNESGetSNESLineSearch(snes,&snes->linesearch);CHKERRQ(ierr); 4013 ierr = PetscObjectSetOptionsPrefix((PetscObject)snes->linesearch,prefix);CHKERRQ(ierr); 4014 } 4015 ierr = KSPSetOptionsPrefix(snes->ksp,prefix);CHKERRQ(ierr); 4016 PetscFunctionReturn(0); 4017 } 4018 4019 #undef __FUNCT__ 4020 #define __FUNCT__ "SNESAppendOptionsPrefix" 4021 /*@C 4022 SNESAppendOptionsPrefix - Appends to the prefix used for searching for all 4023 SNES options in the database. 4024 4025 Logically Collective on SNES 4026 4027 Input Parameters: 4028 + snes - the SNES context 4029 - prefix - the prefix to prepend to all option names 4030 4031 Notes: 4032 A hyphen (-) must NOT be given at the beginning of the prefix name. 4033 The first character of all runtime options is AUTOMATICALLY the hyphen. 4034 4035 Level: advanced 4036 4037 .keywords: SNES, append, options, prefix, database 4038 4039 .seealso: SNESGetOptionsPrefix() 4040 @*/ 4041 PetscErrorCode SNESAppendOptionsPrefix(SNES snes,const char prefix[]) 4042 { 4043 PetscErrorCode ierr; 4044 4045 PetscFunctionBegin; 4046 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4047 ierr = PetscObjectAppendOptionsPrefix((PetscObject)snes,prefix);CHKERRQ(ierr); 4048 if (!snes->ksp) {ierr = SNESGetKSP(snes,&snes->ksp);CHKERRQ(ierr);} 4049 if (snes->linesearch) { 4050 ierr = SNESGetSNESLineSearch(snes,&snes->linesearch);CHKERRQ(ierr); 4051 ierr = PetscObjectAppendOptionsPrefix((PetscObject)snes->linesearch,prefix);CHKERRQ(ierr); 4052 } 4053 ierr = KSPAppendOptionsPrefix(snes->ksp,prefix);CHKERRQ(ierr); 4054 PetscFunctionReturn(0); 4055 } 4056 4057 #undef __FUNCT__ 4058 #define __FUNCT__ "SNESGetOptionsPrefix" 4059 /*@C 4060 SNESGetOptionsPrefix - Sets the prefix used for searching for all 4061 SNES options in the database. 4062 4063 Not Collective 4064 4065 Input Parameter: 4066 . snes - the SNES context 4067 4068 Output Parameter: 4069 . prefix - pointer to the prefix string used 4070 4071 Notes: On the fortran side, the user should pass in a string 'prefix' of 4072 sufficient length to hold the prefix. 4073 4074 Level: advanced 4075 4076 .keywords: SNES, get, options, prefix, database 4077 4078 .seealso: SNESAppendOptionsPrefix() 4079 @*/ 4080 PetscErrorCode SNESGetOptionsPrefix(SNES snes,const char *prefix[]) 4081 { 4082 PetscErrorCode ierr; 4083 4084 PetscFunctionBegin; 4085 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4086 ierr = PetscObjectGetOptionsPrefix((PetscObject)snes,prefix);CHKERRQ(ierr); 4087 PetscFunctionReturn(0); 4088 } 4089 4090 4091 #undef __FUNCT__ 4092 #define __FUNCT__ "SNESRegister" 4093 /*@C 4094 SNESRegister - See SNESRegisterDynamic() 4095 4096 Level: advanced 4097 @*/ 4098 PetscErrorCode SNESRegister(const char sname[],const char path[],const char name[],PetscErrorCode (*function)(SNES)) 4099 { 4100 char fullname[PETSC_MAX_PATH_LEN]; 4101 PetscErrorCode ierr; 4102 4103 PetscFunctionBegin; 4104 ierr = PetscFListConcat(path,name,fullname);CHKERRQ(ierr); 4105 ierr = PetscFListAdd(&SNESList,sname,fullname,(void (*)(void))function);CHKERRQ(ierr); 4106 PetscFunctionReturn(0); 4107 } 4108 4109 #undef __FUNCT__ 4110 #define __FUNCT__ "SNESTestLocalMin" 4111 PetscErrorCode SNESTestLocalMin(SNES snes) 4112 { 4113 PetscErrorCode ierr; 4114 PetscInt N,i,j; 4115 Vec u,uh,fh; 4116 PetscScalar value; 4117 PetscReal norm; 4118 4119 PetscFunctionBegin; 4120 ierr = SNESGetSolution(snes,&u);CHKERRQ(ierr); 4121 ierr = VecDuplicate(u,&uh);CHKERRQ(ierr); 4122 ierr = VecDuplicate(u,&fh);CHKERRQ(ierr); 4123 4124 /* currently only works for sequential */ 4125 ierr = PetscPrintf(PETSC_COMM_WORLD,"Testing FormFunction() for local min\n"); 4126 ierr = VecGetSize(u,&N);CHKERRQ(ierr); 4127 for (i=0; i<N; i++) { 4128 ierr = VecCopy(u,uh);CHKERRQ(ierr); 4129 ierr = PetscPrintf(PETSC_COMM_WORLD,"i = %D\n",i);CHKERRQ(ierr); 4130 for (j=-10; j<11; j++) { 4131 value = PetscSign(j)*exp(PetscAbs(j)-10.0); 4132 ierr = VecSetValue(uh,i,value,ADD_VALUES);CHKERRQ(ierr); 4133 ierr = SNESComputeFunction(snes,uh,fh);CHKERRQ(ierr); 4134 ierr = VecNorm(fh,NORM_2,&norm);CHKERRQ(ierr); 4135 ierr = PetscPrintf(PETSC_COMM_WORLD," j norm %D %18.16e\n",j,norm);CHKERRQ(ierr); 4136 value = -value; 4137 ierr = VecSetValue(uh,i,value,ADD_VALUES);CHKERRQ(ierr); 4138 } 4139 } 4140 ierr = VecDestroy(&uh);CHKERRQ(ierr); 4141 ierr = VecDestroy(&fh);CHKERRQ(ierr); 4142 PetscFunctionReturn(0); 4143 } 4144 4145 #undef __FUNCT__ 4146 #define __FUNCT__ "SNESKSPSetUseEW" 4147 /*@ 4148 SNESKSPSetUseEW - Sets SNES use Eisenstat-Walker method for 4149 computing relative tolerance for linear solvers within an inexact 4150 Newton method. 4151 4152 Logically Collective on SNES 4153 4154 Input Parameters: 4155 + snes - SNES context 4156 - flag - PETSC_TRUE or PETSC_FALSE 4157 4158 Options Database: 4159 + -snes_ksp_ew - use Eisenstat-Walker method for determining linear system convergence 4160 . -snes_ksp_ew_version ver - version of Eisenstat-Walker method 4161 . -snes_ksp_ew_rtol0 <rtol0> - Sets rtol0 4162 . -snes_ksp_ew_rtolmax <rtolmax> - Sets rtolmax 4163 . -snes_ksp_ew_gamma <gamma> - Sets gamma 4164 . -snes_ksp_ew_alpha <alpha> - Sets alpha 4165 . -snes_ksp_ew_alpha2 <alpha2> - Sets alpha2 4166 - -snes_ksp_ew_threshold <threshold> - Sets threshold 4167 4168 Notes: 4169 Currently, the default is to use a constant relative tolerance for 4170 the inner linear solvers. Alternatively, one can use the 4171 Eisenstat-Walker method, where the relative convergence tolerance 4172 is reset at each Newton iteration according progress of the nonlinear 4173 solver. 4174 4175 Level: advanced 4176 4177 Reference: 4178 S. C. Eisenstat and H. F. Walker, "Choosing the forcing terms in an 4179 inexact Newton method", SISC 17 (1), pp.16-32, 1996. 4180 4181 .keywords: SNES, KSP, Eisenstat, Walker, convergence, test, inexact, Newton 4182 4183 .seealso: SNESKSPGetUseEW(), SNESKSPGetParametersEW(), SNESKSPSetParametersEW() 4184 @*/ 4185 PetscErrorCode SNESKSPSetUseEW(SNES snes,PetscBool flag) 4186 { 4187 PetscFunctionBegin; 4188 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4189 PetscValidLogicalCollectiveBool(snes,flag,2); 4190 snes->ksp_ewconv = flag; 4191 PetscFunctionReturn(0); 4192 } 4193 4194 #undef __FUNCT__ 4195 #define __FUNCT__ "SNESKSPGetUseEW" 4196 /*@ 4197 SNESKSPGetUseEW - Gets if SNES is using Eisenstat-Walker method 4198 for computing relative tolerance for linear solvers within an 4199 inexact Newton method. 4200 4201 Not Collective 4202 4203 Input Parameter: 4204 . snes - SNES context 4205 4206 Output Parameter: 4207 . flag - PETSC_TRUE or PETSC_FALSE 4208 4209 Notes: 4210 Currently, the default is to use a constant relative tolerance for 4211 the inner linear solvers. Alternatively, one can use the 4212 Eisenstat-Walker method, where the relative convergence tolerance 4213 is reset at each Newton iteration according progress of the nonlinear 4214 solver. 4215 4216 Level: advanced 4217 4218 Reference: 4219 S. C. Eisenstat and H. F. Walker, "Choosing the forcing terms in an 4220 inexact Newton method", SISC 17 (1), pp.16-32, 1996. 4221 4222 .keywords: SNES, KSP, Eisenstat, Walker, convergence, test, inexact, Newton 4223 4224 .seealso: SNESKSPSetUseEW(), SNESKSPGetParametersEW(), SNESKSPSetParametersEW() 4225 @*/ 4226 PetscErrorCode SNESKSPGetUseEW(SNES snes, PetscBool *flag) 4227 { 4228 PetscFunctionBegin; 4229 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4230 PetscValidPointer(flag,2); 4231 *flag = snes->ksp_ewconv; 4232 PetscFunctionReturn(0); 4233 } 4234 4235 #undef __FUNCT__ 4236 #define __FUNCT__ "SNESKSPSetParametersEW" 4237 /*@ 4238 SNESKSPSetParametersEW - Sets parameters for Eisenstat-Walker 4239 convergence criteria for the linear solvers within an inexact 4240 Newton method. 4241 4242 Logically Collective on SNES 4243 4244 Input Parameters: 4245 + snes - SNES context 4246 . version - version 1, 2 (default is 2) or 3 4247 . rtol_0 - initial relative tolerance (0 <= rtol_0 < 1) 4248 . rtol_max - maximum relative tolerance (0 <= rtol_max < 1) 4249 . gamma - multiplicative factor for version 2 rtol computation 4250 (0 <= gamma2 <= 1) 4251 . alpha - power for version 2 rtol computation (1 < alpha <= 2) 4252 . alpha2 - power for safeguard 4253 - threshold - threshold for imposing safeguard (0 < threshold < 1) 4254 4255 Note: 4256 Version 3 was contributed by Luis Chacon, June 2006. 4257 4258 Use PETSC_DEFAULT to retain the default for any of the parameters. 4259 4260 Level: advanced 4261 4262 Reference: 4263 S. C. Eisenstat and H. F. Walker, "Choosing the forcing terms in an 4264 inexact Newton method", Utah State University Math. Stat. Dept. Res. 4265 Report 6/94/75, June, 1994, to appear in SIAM J. Sci. Comput. 4266 4267 .keywords: SNES, KSP, Eisenstat, Walker, set, parameters 4268 4269 .seealso: SNESKSPSetUseEW(), SNESKSPGetUseEW(), SNESKSPGetParametersEW() 4270 @*/ 4271 PetscErrorCode SNESKSPSetParametersEW(SNES snes,PetscInt version,PetscReal rtol_0,PetscReal rtol_max, 4272 PetscReal gamma,PetscReal alpha,PetscReal alpha2,PetscReal threshold) 4273 { 4274 SNESKSPEW *kctx; 4275 PetscFunctionBegin; 4276 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4277 kctx = (SNESKSPEW*)snes->kspconvctx; 4278 if (!kctx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"No Eisenstat-Walker context existing"); 4279 PetscValidLogicalCollectiveInt(snes,version,2); 4280 PetscValidLogicalCollectiveReal(snes,rtol_0,3); 4281 PetscValidLogicalCollectiveReal(snes,rtol_max,4); 4282 PetscValidLogicalCollectiveReal(snes,gamma,5); 4283 PetscValidLogicalCollectiveReal(snes,alpha,6); 4284 PetscValidLogicalCollectiveReal(snes,alpha2,7); 4285 PetscValidLogicalCollectiveReal(snes,threshold,8); 4286 4287 if (version != PETSC_DEFAULT) kctx->version = version; 4288 if (rtol_0 != PETSC_DEFAULT) kctx->rtol_0 = rtol_0; 4289 if (rtol_max != PETSC_DEFAULT) kctx->rtol_max = rtol_max; 4290 if (gamma != PETSC_DEFAULT) kctx->gamma = gamma; 4291 if (alpha != PETSC_DEFAULT) kctx->alpha = alpha; 4292 if (alpha2 != PETSC_DEFAULT) kctx->alpha2 = alpha2; 4293 if (threshold != PETSC_DEFAULT) kctx->threshold = threshold; 4294 4295 if (kctx->version < 1 || kctx->version > 3) { 4296 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Only versions 1, 2 and 3 are supported: %D",kctx->version); 4297 } 4298 if (kctx->rtol_0 < 0.0 || kctx->rtol_0 >= 1.0) { 4299 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"0.0 <= rtol_0 < 1.0: %G",kctx->rtol_0); 4300 } 4301 if (kctx->rtol_max < 0.0 || kctx->rtol_max >= 1.0) { 4302 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"0.0 <= rtol_max (%G) < 1.0\n",kctx->rtol_max); 4303 } 4304 if (kctx->gamma < 0.0 || kctx->gamma > 1.0) { 4305 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"0.0 <= gamma (%G) <= 1.0\n",kctx->gamma); 4306 } 4307 if (kctx->alpha <= 1.0 || kctx->alpha > 2.0) { 4308 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"1.0 < alpha (%G) <= 2.0\n",kctx->alpha); 4309 } 4310 if (kctx->threshold <= 0.0 || kctx->threshold >= 1.0) { 4311 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"0.0 < threshold (%G) < 1.0\n",kctx->threshold); 4312 } 4313 PetscFunctionReturn(0); 4314 } 4315 4316 #undef __FUNCT__ 4317 #define __FUNCT__ "SNESKSPGetParametersEW" 4318 /*@ 4319 SNESKSPGetParametersEW - Gets parameters for Eisenstat-Walker 4320 convergence criteria for the linear solvers within an inexact 4321 Newton method. 4322 4323 Not Collective 4324 4325 Input Parameters: 4326 snes - SNES context 4327 4328 Output Parameters: 4329 + version - version 1, 2 (default is 2) or 3 4330 . rtol_0 - initial relative tolerance (0 <= rtol_0 < 1) 4331 . rtol_max - maximum relative tolerance (0 <= rtol_max < 1) 4332 . gamma - multiplicative factor for version 2 rtol computation 4333 (0 <= gamma2 <= 1) 4334 . alpha - power for version 2 rtol computation (1 < alpha <= 2) 4335 . alpha2 - power for safeguard 4336 - threshold - threshold for imposing safeguard (0 < threshold < 1) 4337 4338 Level: advanced 4339 4340 .keywords: SNES, KSP, Eisenstat, Walker, get, parameters 4341 4342 .seealso: SNESKSPSetUseEW(), SNESKSPGetUseEW(), SNESKSPSetParametersEW() 4343 @*/ 4344 PetscErrorCode SNESKSPGetParametersEW(SNES snes,PetscInt *version,PetscReal *rtol_0,PetscReal *rtol_max, 4345 PetscReal *gamma,PetscReal *alpha,PetscReal *alpha2,PetscReal *threshold) 4346 { 4347 SNESKSPEW *kctx; 4348 PetscFunctionBegin; 4349 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4350 kctx = (SNESKSPEW*)snes->kspconvctx; 4351 if (!kctx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"No Eisenstat-Walker context existing"); 4352 if (version) *version = kctx->version; 4353 if (rtol_0) *rtol_0 = kctx->rtol_0; 4354 if (rtol_max) *rtol_max = kctx->rtol_max; 4355 if (gamma) *gamma = kctx->gamma; 4356 if (alpha) *alpha = kctx->alpha; 4357 if (alpha2) *alpha2 = kctx->alpha2; 4358 if (threshold) *threshold = kctx->threshold; 4359 PetscFunctionReturn(0); 4360 } 4361 4362 #undef __FUNCT__ 4363 #define __FUNCT__ "SNESKSPEW_PreSolve" 4364 static PetscErrorCode SNESKSPEW_PreSolve(SNES snes, KSP ksp, Vec b, Vec x) 4365 { 4366 PetscErrorCode ierr; 4367 SNESKSPEW *kctx = (SNESKSPEW*)snes->kspconvctx; 4368 PetscReal rtol=PETSC_DEFAULT,stol; 4369 4370 PetscFunctionBegin; 4371 if (!kctx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"No Eisenstat-Walker context exists"); 4372 if (!snes->iter) { /* first time in, so use the original user rtol */ 4373 rtol = kctx->rtol_0; 4374 } else { 4375 if (kctx->version == 1) { 4376 rtol = (snes->norm - kctx->lresid_last)/kctx->norm_last; 4377 if (rtol < 0.0) rtol = -rtol; 4378 stol = pow(kctx->rtol_last,kctx->alpha2); 4379 if (stol > kctx->threshold) rtol = PetscMax(rtol,stol); 4380 } else if (kctx->version == 2) { 4381 rtol = kctx->gamma * pow(snes->norm/kctx->norm_last,kctx->alpha); 4382 stol = kctx->gamma * pow(kctx->rtol_last,kctx->alpha); 4383 if (stol > kctx->threshold) rtol = PetscMax(rtol,stol); 4384 } else if (kctx->version == 3) {/* contributed by Luis Chacon, June 2006. */ 4385 rtol = kctx->gamma * pow(snes->norm/kctx->norm_last,kctx->alpha); 4386 /* safeguard: avoid sharp decrease of rtol */ 4387 stol = kctx->gamma*pow(kctx->rtol_last,kctx->alpha); 4388 stol = PetscMax(rtol,stol); 4389 rtol = PetscMin(kctx->rtol_0,stol); 4390 /* safeguard: avoid oversolving */ 4391 stol = kctx->gamma*(snes->ttol)/snes->norm; 4392 stol = PetscMax(rtol,stol); 4393 rtol = PetscMin(kctx->rtol_0,stol); 4394 } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Only versions 1, 2 or 3 are supported: %D",kctx->version); 4395 } 4396 /* safeguard: avoid rtol greater than one */ 4397 rtol = PetscMin(rtol,kctx->rtol_max); 4398 ierr = KSPSetTolerances(ksp,rtol,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);CHKERRQ(ierr); 4399 ierr = PetscInfo3(snes,"iter %D, Eisenstat-Walker (version %D) KSP rtol=%G\n",snes->iter,kctx->version,rtol);CHKERRQ(ierr); 4400 PetscFunctionReturn(0); 4401 } 4402 4403 #undef __FUNCT__ 4404 #define __FUNCT__ "SNESKSPEW_PostSolve" 4405 static PetscErrorCode SNESKSPEW_PostSolve(SNES snes, KSP ksp, Vec b, Vec x) 4406 { 4407 PetscErrorCode ierr; 4408 SNESKSPEW *kctx = (SNESKSPEW*)snes->kspconvctx; 4409 PCSide pcside; 4410 Vec lres; 4411 4412 PetscFunctionBegin; 4413 if (!kctx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"No Eisenstat-Walker context exists"); 4414 ierr = KSPGetTolerances(ksp,&kctx->rtol_last,0,0,0);CHKERRQ(ierr); 4415 ierr = SNESGetFunctionNorm(snes,&kctx->norm_last);CHKERRQ(ierr); 4416 if (kctx->version == 1) { 4417 ierr = KSPGetPCSide(ksp,&pcside);CHKERRQ(ierr); 4418 if (pcside == PC_RIGHT) { /* XXX Should we also test KSP_UNPRECONDITIONED_NORM ? */ 4419 /* KSP residual is true linear residual */ 4420 ierr = KSPGetResidualNorm(ksp,&kctx->lresid_last);CHKERRQ(ierr); 4421 } else { 4422 /* KSP residual is preconditioned residual */ 4423 /* compute true linear residual norm */ 4424 ierr = VecDuplicate(b,&lres);CHKERRQ(ierr); 4425 ierr = MatMult(snes->jacobian,x,lres);CHKERRQ(ierr); 4426 ierr = VecAYPX(lres,-1.0,b);CHKERRQ(ierr); 4427 ierr = VecNorm(lres,NORM_2,&kctx->lresid_last);CHKERRQ(ierr); 4428 ierr = VecDestroy(&lres);CHKERRQ(ierr); 4429 } 4430 } 4431 PetscFunctionReturn(0); 4432 } 4433 4434 #undef __FUNCT__ 4435 #define __FUNCT__ "SNES_KSPSolve" 4436 PetscErrorCode SNES_KSPSolve(SNES snes, KSP ksp, Vec b, Vec x) 4437 { 4438 PetscErrorCode ierr; 4439 4440 PetscFunctionBegin; 4441 if (snes->ksp_ewconv) { ierr = SNESKSPEW_PreSolve(snes,ksp,b,x);CHKERRQ(ierr); } 4442 ierr = KSPSolve(ksp,b,x);CHKERRQ(ierr); 4443 if (snes->ksp_ewconv) { ierr = SNESKSPEW_PostSolve(snes,ksp,b,x);CHKERRQ(ierr); } 4444 PetscFunctionReturn(0); 4445 } 4446 4447 #undef __FUNCT__ 4448 #define __FUNCT__ "SNESSetDM" 4449 /*@ 4450 SNESSetDM - Sets the DM that may be used by some preconditioners 4451 4452 Logically Collective on SNES 4453 4454 Input Parameters: 4455 + snes - the preconditioner context 4456 - dm - the dm 4457 4458 Level: intermediate 4459 4460 4461 .seealso: SNESGetDM(), KSPSetDM(), KSPGetDM() 4462 @*/ 4463 PetscErrorCode SNESSetDM(SNES snes,DM dm) 4464 { 4465 PetscErrorCode ierr; 4466 KSP ksp; 4467 SNESDM sdm; 4468 4469 PetscFunctionBegin; 4470 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4471 if (dm) {ierr = PetscObjectReference((PetscObject)dm);CHKERRQ(ierr);} 4472 if (snes->dm) { /* Move the SNESDM context over to the new DM unless the new DM already has one */ 4473 PetscContainer oldcontainer,container; 4474 ierr = PetscObjectQuery((PetscObject)snes->dm,"SNESDM",(PetscObject*)&oldcontainer);CHKERRQ(ierr); 4475 ierr = PetscObjectQuery((PetscObject)dm,"SNESDM",(PetscObject*)&container);CHKERRQ(ierr); 4476 if (oldcontainer && snes->dmAuto && !container) { 4477 ierr = DMSNESCopyContext(snes->dm,dm);CHKERRQ(ierr); 4478 ierr = DMSNESGetContext(snes->dm,&sdm);CHKERRQ(ierr); 4479 if (sdm->originaldm == snes->dm) { /* Grant write privileges to the replacement DM */ 4480 sdm->originaldm = dm; 4481 } 4482 } 4483 ierr = DMDestroy(&snes->dm);CHKERRQ(ierr); 4484 } 4485 snes->dm = dm; 4486 snes->dmAuto = PETSC_FALSE; 4487 ierr = SNESGetKSP(snes,&ksp);CHKERRQ(ierr); 4488 ierr = KSPSetDM(ksp,dm);CHKERRQ(ierr); 4489 ierr = KSPSetDMActive(ksp,PETSC_FALSE);CHKERRQ(ierr); 4490 if (snes->pc) { 4491 ierr = SNESSetDM(snes->pc, snes->dm);CHKERRQ(ierr); 4492 ierr = SNESSetPCSide(snes,snes->pcside);CHKERRQ(ierr); 4493 } 4494 PetscFunctionReturn(0); 4495 } 4496 4497 #undef __FUNCT__ 4498 #define __FUNCT__ "SNESGetDM" 4499 /*@ 4500 SNESGetDM - Gets the DM that may be used by some preconditioners 4501 4502 Not Collective but DM obtained is parallel on SNES 4503 4504 Input Parameter: 4505 . snes - the preconditioner context 4506 4507 Output Parameter: 4508 . dm - the dm 4509 4510 Level: intermediate 4511 4512 4513 .seealso: SNESSetDM(), KSPSetDM(), KSPGetDM() 4514 @*/ 4515 PetscErrorCode SNESGetDM(SNES snes,DM *dm) 4516 { 4517 PetscErrorCode ierr; 4518 4519 PetscFunctionBegin; 4520 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4521 if (!snes->dm) { 4522 ierr = DMShellCreate(((PetscObject)snes)->comm,&snes->dm);CHKERRQ(ierr); 4523 snes->dmAuto = PETSC_TRUE; 4524 } 4525 *dm = snes->dm; 4526 PetscFunctionReturn(0); 4527 } 4528 4529 #undef __FUNCT__ 4530 #define __FUNCT__ "SNESSetPC" 4531 /*@ 4532 SNESSetPC - Sets the nonlinear preconditioner to be used. 4533 4534 Collective on SNES 4535 4536 Input Parameters: 4537 + snes - iterative context obtained from SNESCreate() 4538 - pc - the preconditioner object 4539 4540 Notes: 4541 Use SNESGetPC() to retrieve the preconditioner context (for example, 4542 to configure it using the API). 4543 4544 Level: developer 4545 4546 .keywords: SNES, set, precondition 4547 .seealso: SNESGetPC() 4548 @*/ 4549 PetscErrorCode SNESSetPC(SNES snes, SNES pc) 4550 { 4551 PetscErrorCode ierr; 4552 4553 PetscFunctionBegin; 4554 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); 4555 PetscValidHeaderSpecific(pc, SNES_CLASSID, 2); 4556 PetscCheckSameComm(snes, 1, pc, 2); 4557 ierr = PetscObjectReference((PetscObject) pc);CHKERRQ(ierr); 4558 ierr = SNESDestroy(&snes->pc);CHKERRQ(ierr); 4559 snes->pc = pc; 4560 ierr = PetscLogObjectParent(snes, snes->pc);CHKERRQ(ierr); 4561 PetscFunctionReturn(0); 4562 } 4563 4564 #undef __FUNCT__ 4565 #define __FUNCT__ "SNESGetPC" 4566 /*@ 4567 SNESGetPC - Returns a pointer to the nonlinear preconditioning context set with SNESSetPC(). 4568 4569 Not Collective 4570 4571 Input Parameter: 4572 . snes - iterative context obtained from SNESCreate() 4573 4574 Output Parameter: 4575 . pc - preconditioner context 4576 4577 Level: developer 4578 4579 .keywords: SNES, get, preconditioner 4580 .seealso: SNESSetPC() 4581 @*/ 4582 PetscErrorCode SNESGetPC(SNES snes, SNES *pc) 4583 { 4584 PetscErrorCode ierr; 4585 const char *optionsprefix; 4586 4587 PetscFunctionBegin; 4588 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); 4589 PetscValidPointer(pc, 2); 4590 if (!snes->pc) { 4591 ierr = SNESCreate(((PetscObject) snes)->comm,&snes->pc);CHKERRQ(ierr); 4592 ierr = PetscObjectIncrementTabLevel((PetscObject)snes->pc,(PetscObject)snes,1);CHKERRQ(ierr); 4593 ierr = PetscLogObjectParent(snes,snes->pc);CHKERRQ(ierr); 4594 ierr = SNESGetOptionsPrefix(snes,&optionsprefix);CHKERRQ(ierr); 4595 ierr = SNESSetOptionsPrefix(snes->pc,optionsprefix);CHKERRQ(ierr); 4596 ierr = SNESAppendOptionsPrefix(snes->pc,"npc_");CHKERRQ(ierr); 4597 } 4598 *pc = snes->pc; 4599 PetscFunctionReturn(0); 4600 } 4601 4602 4603 #undef __FUNCT__ 4604 #define __FUNCT__ "SNESSetPCSide" 4605 /*@ 4606 SNESSetPCSide - Sets the preconditioning side. 4607 4608 Logically Collective on SNES 4609 4610 Input Parameter: 4611 . snes - iterative context obtained from SNESCreate() 4612 4613 Output Parameter: 4614 . side - the preconditioning side, where side is one of 4615 .vb 4616 PC_LEFT - left preconditioning (default) 4617 PC_RIGHT - right preconditioning 4618 .ve 4619 4620 Options Database Keys: 4621 . -snes_pc_side <right,left> 4622 4623 Level: intermediate 4624 4625 .keywords: SNES, set, right, left, side, preconditioner, flag 4626 4627 .seealso: SNESGetPCSide(), KSPSetPCSide() 4628 @*/ 4629 PetscErrorCode SNESSetPCSide(SNES snes,PCSide side) 4630 { 4631 PetscFunctionBegin; 4632 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4633 PetscValidLogicalCollectiveEnum(snes,side,2); 4634 snes->pcside = side; 4635 PetscFunctionReturn(0); 4636 } 4637 4638 #undef __FUNCT__ 4639 #define __FUNCT__ "SNESGetPCSide" 4640 /*@ 4641 SNESGetPCSide - Gets the preconditioning side. 4642 4643 Not Collective 4644 4645 Input Parameter: 4646 . snes - iterative context obtained from SNESCreate() 4647 4648 Output Parameter: 4649 . side - the preconditioning side, where side is one of 4650 .vb 4651 PC_LEFT - left preconditioning (default) 4652 PC_RIGHT - right preconditioning 4653 .ve 4654 4655 Level: intermediate 4656 4657 .keywords: SNES, get, right, left, side, preconditioner, flag 4658 4659 .seealso: SNESSetPCSide(), KSPGetPCSide() 4660 @*/ 4661 PetscErrorCode SNESGetPCSide(SNES snes,PCSide *side) 4662 { 4663 PetscFunctionBegin; 4664 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4665 PetscValidPointer(side,2); 4666 *side = snes->pcside; 4667 PetscFunctionReturn(0); 4668 } 4669 4670 #undef __FUNCT__ 4671 #define __FUNCT__ "SNESSetSNESLineSearch" 4672 /*@ 4673 SNESSetSNESLineSearch - Sets the linesearch on the SNES instance. 4674 4675 Collective on SNES 4676 4677 Input Parameters: 4678 + snes - iterative context obtained from SNESCreate() 4679 - linesearch - the linesearch object 4680 4681 Notes: 4682 Use SNESGetSNESLineSearch() to retrieve the preconditioner context (for example, 4683 to configure it using the API). 4684 4685 Level: developer 4686 4687 .keywords: SNES, set, linesearch 4688 .seealso: SNESGetSNESLineSearch() 4689 @*/ 4690 PetscErrorCode SNESSetSNESLineSearch(SNES snes, SNESLineSearch linesearch) 4691 { 4692 PetscErrorCode ierr; 4693 4694 PetscFunctionBegin; 4695 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); 4696 PetscValidHeaderSpecific(linesearch, SNESLINESEARCH_CLASSID, 2); 4697 PetscCheckSameComm(snes, 1, linesearch, 2); 4698 ierr = PetscObjectReference((PetscObject) linesearch);CHKERRQ(ierr); 4699 ierr = SNESLineSearchDestroy(&snes->linesearch);CHKERRQ(ierr); 4700 snes->linesearch = linesearch; 4701 ierr = PetscLogObjectParent(snes, snes->linesearch);CHKERRQ(ierr); 4702 PetscFunctionReturn(0); 4703 } 4704 4705 #undef __FUNCT__ 4706 #define __FUNCT__ "SNESGetSNESLineSearch" 4707 /*@C 4708 SNESGetSNESLineSearch - Returns a pointer to the line search context set with SNESSetLineSearch() 4709 or creates a default line search instance associated with the SNES and returns it. 4710 4711 Not Collective 4712 4713 Input Parameter: 4714 . snes - iterative context obtained from SNESCreate() 4715 4716 Output Parameter: 4717 . linesearch - linesearch context 4718 4719 Level: developer 4720 4721 .keywords: SNES, get, linesearch 4722 .seealso: SNESSetSNESLineSearch() 4723 @*/ 4724 PetscErrorCode SNESGetSNESLineSearch(SNES snes, SNESLineSearch *linesearch) 4725 { 4726 PetscErrorCode ierr; 4727 const char *optionsprefix; 4728 4729 PetscFunctionBegin; 4730 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); 4731 PetscValidPointer(linesearch, 2); 4732 if (!snes->linesearch) { 4733 ierr = SNESGetOptionsPrefix(snes, &optionsprefix);CHKERRQ(ierr); 4734 ierr = SNESLineSearchCreate(((PetscObject) snes)->comm, &snes->linesearch);CHKERRQ(ierr); 4735 ierr = SNESLineSearchSetSNES(snes->linesearch, snes);CHKERRQ(ierr); 4736 ierr = SNESLineSearchAppendOptionsPrefix(snes->linesearch, optionsprefix);CHKERRQ(ierr); 4737 ierr = PetscObjectIncrementTabLevel((PetscObject) snes->linesearch, (PetscObject) snes, 1);CHKERRQ(ierr); 4738 ierr = PetscLogObjectParent(snes, snes->linesearch);CHKERRQ(ierr); 4739 } 4740 *linesearch = snes->linesearch; 4741 PetscFunctionReturn(0); 4742 } 4743 4744 #if defined(PETSC_HAVE_MATLAB_ENGINE) 4745 #include <mex.h> 4746 4747 typedef struct {char *funcname; mxArray *ctx;} SNESMatlabContext; 4748 4749 #undef __FUNCT__ 4750 #define __FUNCT__ "SNESComputeFunction_Matlab" 4751 /* 4752 SNESComputeFunction_Matlab - Calls the function that has been set with 4753 SNESSetFunctionMatlab(). 4754 4755 Collective on SNES 4756 4757 Input Parameters: 4758 + snes - the SNES context 4759 - x - input vector 4760 4761 Output Parameter: 4762 . y - function vector, as set by SNESSetFunction() 4763 4764 Notes: 4765 SNESComputeFunction() is typically used within nonlinear solvers 4766 implementations, so most users would not generally call this routine 4767 themselves. 4768 4769 Level: developer 4770 4771 .keywords: SNES, nonlinear, compute, function 4772 4773 .seealso: SNESSetFunction(), SNESGetFunction() 4774 */ 4775 PetscErrorCode SNESComputeFunction_Matlab(SNES snes,Vec x,Vec y, void *ctx) 4776 { 4777 PetscErrorCode ierr; 4778 SNESMatlabContext *sctx = (SNESMatlabContext *)ctx; 4779 int nlhs = 1,nrhs = 5; 4780 mxArray *plhs[1],*prhs[5]; 4781 long long int lx = 0,ly = 0,ls = 0; 4782 4783 PetscFunctionBegin; 4784 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4785 PetscValidHeaderSpecific(x,VEC_CLASSID,2); 4786 PetscValidHeaderSpecific(y,VEC_CLASSID,3); 4787 PetscCheckSameComm(snes,1,x,2); 4788 PetscCheckSameComm(snes,1,y,3); 4789 4790 /* call Matlab function in ctx with arguments x and y */ 4791 4792 ierr = PetscMemcpy(&ls,&snes,sizeof(snes));CHKERRQ(ierr); 4793 ierr = PetscMemcpy(&lx,&x,sizeof(x));CHKERRQ(ierr); 4794 ierr = PetscMemcpy(&ly,&y,sizeof(x));CHKERRQ(ierr); 4795 prhs[0] = mxCreateDoubleScalar((double)ls); 4796 prhs[1] = mxCreateDoubleScalar((double)lx); 4797 prhs[2] = mxCreateDoubleScalar((double)ly); 4798 prhs[3] = mxCreateString(sctx->funcname); 4799 prhs[4] = sctx->ctx; 4800 ierr = mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscSNESComputeFunctionInternal");CHKERRQ(ierr); 4801 ierr = mxGetScalar(plhs[0]);CHKERRQ(ierr); 4802 mxDestroyArray(prhs[0]); 4803 mxDestroyArray(prhs[1]); 4804 mxDestroyArray(prhs[2]); 4805 mxDestroyArray(prhs[3]); 4806 mxDestroyArray(plhs[0]); 4807 PetscFunctionReturn(0); 4808 } 4809 4810 4811 #undef __FUNCT__ 4812 #define __FUNCT__ "SNESSetFunctionMatlab" 4813 /* 4814 SNESSetFunctionMatlab - Sets the function evaluation routine and function 4815 vector for use by the SNES routines in solving systems of nonlinear 4816 equations from MATLAB. Here the function is a string containing the name of a MATLAB function 4817 4818 Logically Collective on SNES 4819 4820 Input Parameters: 4821 + snes - the SNES context 4822 . r - vector to store function value 4823 - func - function evaluation routine 4824 4825 Calling sequence of func: 4826 $ func (SNES snes,Vec x,Vec f,void *ctx); 4827 4828 4829 Notes: 4830 The Newton-like methods typically solve linear systems of the form 4831 $ f'(x) x = -f(x), 4832 where f'(x) denotes the Jacobian matrix and f(x) is the function. 4833 4834 Level: beginner 4835 4836 .keywords: SNES, nonlinear, set, function 4837 4838 .seealso: SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESSetFunction() 4839 */ 4840 PetscErrorCode SNESSetFunctionMatlab(SNES snes,Vec r,const char *func,mxArray *ctx) 4841 { 4842 PetscErrorCode ierr; 4843 SNESMatlabContext *sctx; 4844 4845 PetscFunctionBegin; 4846 /* currently sctx is memory bleed */ 4847 ierr = PetscMalloc(sizeof(SNESMatlabContext),&sctx);CHKERRQ(ierr); 4848 ierr = PetscStrallocpy(func,&sctx->funcname);CHKERRQ(ierr); 4849 /* 4850 This should work, but it doesn't 4851 sctx->ctx = ctx; 4852 mexMakeArrayPersistent(sctx->ctx); 4853 */ 4854 sctx->ctx = mxDuplicateArray(ctx); 4855 ierr = SNESSetFunction(snes,r,SNESComputeFunction_Matlab,sctx);CHKERRQ(ierr); 4856 PetscFunctionReturn(0); 4857 } 4858 4859 #undef __FUNCT__ 4860 #define __FUNCT__ "SNESComputeJacobian_Matlab" 4861 /* 4862 SNESComputeJacobian_Matlab - Calls the function that has been set with 4863 SNESSetJacobianMatlab(). 4864 4865 Collective on SNES 4866 4867 Input Parameters: 4868 + snes - the SNES context 4869 . x - input vector 4870 . A, B - the matrices 4871 - ctx - user context 4872 4873 Output Parameter: 4874 . flag - structure of the matrix 4875 4876 Level: developer 4877 4878 .keywords: SNES, nonlinear, compute, function 4879 4880 .seealso: SNESSetFunction(), SNESGetFunction() 4881 @*/ 4882 PetscErrorCode SNESComputeJacobian_Matlab(SNES snes,Vec x,Mat *A,Mat *B,MatStructure *flag, void *ctx) 4883 { 4884 PetscErrorCode ierr; 4885 SNESMatlabContext *sctx = (SNESMatlabContext *)ctx; 4886 int nlhs = 2,nrhs = 6; 4887 mxArray *plhs[2],*prhs[6]; 4888 long long int lx = 0,lA = 0,ls = 0, lB = 0; 4889 4890 PetscFunctionBegin; 4891 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4892 PetscValidHeaderSpecific(x,VEC_CLASSID,2); 4893 4894 /* call Matlab function in ctx with arguments x and y */ 4895 4896 ierr = PetscMemcpy(&ls,&snes,sizeof(snes));CHKERRQ(ierr); 4897 ierr = PetscMemcpy(&lx,&x,sizeof(x));CHKERRQ(ierr); 4898 ierr = PetscMemcpy(&lA,A,sizeof(x));CHKERRQ(ierr); 4899 ierr = PetscMemcpy(&lB,B,sizeof(x));CHKERRQ(ierr); 4900 prhs[0] = mxCreateDoubleScalar((double)ls); 4901 prhs[1] = mxCreateDoubleScalar((double)lx); 4902 prhs[2] = mxCreateDoubleScalar((double)lA); 4903 prhs[3] = mxCreateDoubleScalar((double)lB); 4904 prhs[4] = mxCreateString(sctx->funcname); 4905 prhs[5] = sctx->ctx; 4906 ierr = mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscSNESComputeJacobianInternal");CHKERRQ(ierr); 4907 ierr = mxGetScalar(plhs[0]);CHKERRQ(ierr); 4908 *flag = (MatStructure) mxGetScalar(plhs[1]);CHKERRQ(ierr); 4909 mxDestroyArray(prhs[0]); 4910 mxDestroyArray(prhs[1]); 4911 mxDestroyArray(prhs[2]); 4912 mxDestroyArray(prhs[3]); 4913 mxDestroyArray(prhs[4]); 4914 mxDestroyArray(plhs[0]); 4915 mxDestroyArray(plhs[1]); 4916 PetscFunctionReturn(0); 4917 } 4918 4919 4920 #undef __FUNCT__ 4921 #define __FUNCT__ "SNESSetJacobianMatlab" 4922 /* 4923 SNESSetJacobianMatlab - Sets the Jacobian function evaluation routine and two empty Jacobian matrices 4924 vector for use by the SNES routines in solving systems of nonlinear 4925 equations from MATLAB. Here the function is a string containing the name of a MATLAB function 4926 4927 Logically Collective on SNES 4928 4929 Input Parameters: 4930 + snes - the SNES context 4931 . A,B - Jacobian matrices 4932 . func - function evaluation routine 4933 - ctx - user context 4934 4935 Calling sequence of func: 4936 $ flag = func (SNES snes,Vec x,Mat A,Mat B,void *ctx); 4937 4938 4939 Level: developer 4940 4941 .keywords: SNES, nonlinear, set, function 4942 4943 .seealso: SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESSetFunction() 4944 */ 4945 PetscErrorCode SNESSetJacobianMatlab(SNES snes,Mat A,Mat B,const char *func,mxArray *ctx) 4946 { 4947 PetscErrorCode ierr; 4948 SNESMatlabContext *sctx; 4949 4950 PetscFunctionBegin; 4951 /* currently sctx is memory bleed */ 4952 ierr = PetscMalloc(sizeof(SNESMatlabContext),&sctx);CHKERRQ(ierr); 4953 ierr = PetscStrallocpy(func,&sctx->funcname);CHKERRQ(ierr); 4954 /* 4955 This should work, but it doesn't 4956 sctx->ctx = ctx; 4957 mexMakeArrayPersistent(sctx->ctx); 4958 */ 4959 sctx->ctx = mxDuplicateArray(ctx); 4960 ierr = SNESSetJacobian(snes,A,B,SNESComputeJacobian_Matlab,sctx);CHKERRQ(ierr); 4961 PetscFunctionReturn(0); 4962 } 4963 4964 #undef __FUNCT__ 4965 #define __FUNCT__ "SNESMonitor_Matlab" 4966 /* 4967 SNESMonitor_Matlab - Calls the function that has been set with SNESMonitorSetMatlab(). 4968 4969 Collective on SNES 4970 4971 .seealso: SNESSetFunction(), SNESGetFunction() 4972 @*/ 4973 PetscErrorCode SNESMonitor_Matlab(SNES snes,PetscInt it, PetscReal fnorm, void *ctx) 4974 { 4975 PetscErrorCode ierr; 4976 SNESMatlabContext *sctx = (SNESMatlabContext *)ctx; 4977 int nlhs = 1,nrhs = 6; 4978 mxArray *plhs[1],*prhs[6]; 4979 long long int lx = 0,ls = 0; 4980 Vec x=snes->vec_sol; 4981 4982 PetscFunctionBegin; 4983 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4984 4985 ierr = PetscMemcpy(&ls,&snes,sizeof(snes));CHKERRQ(ierr); 4986 ierr = PetscMemcpy(&lx,&x,sizeof(x));CHKERRQ(ierr); 4987 prhs[0] = mxCreateDoubleScalar((double)ls); 4988 prhs[1] = mxCreateDoubleScalar((double)it); 4989 prhs[2] = mxCreateDoubleScalar((double)fnorm); 4990 prhs[3] = mxCreateDoubleScalar((double)lx); 4991 prhs[4] = mxCreateString(sctx->funcname); 4992 prhs[5] = sctx->ctx; 4993 ierr = mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscSNESMonitorInternal");CHKERRQ(ierr); 4994 ierr = mxGetScalar(plhs[0]);CHKERRQ(ierr); 4995 mxDestroyArray(prhs[0]); 4996 mxDestroyArray(prhs[1]); 4997 mxDestroyArray(prhs[2]); 4998 mxDestroyArray(prhs[3]); 4999 mxDestroyArray(prhs[4]); 5000 mxDestroyArray(plhs[0]); 5001 PetscFunctionReturn(0); 5002 } 5003 5004 5005 #undef __FUNCT__ 5006 #define __FUNCT__ "SNESMonitorSetMatlab" 5007 /* 5008 SNESMonitorSetMatlab - Sets the monitor function from MATLAB 5009 5010 Level: developer 5011 5012 .keywords: SNES, nonlinear, set, function 5013 5014 .seealso: SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESSetFunction() 5015 */ 5016 PetscErrorCode SNESMonitorSetMatlab(SNES snes,const char *func,mxArray *ctx) 5017 { 5018 PetscErrorCode ierr; 5019 SNESMatlabContext *sctx; 5020 5021 PetscFunctionBegin; 5022 /* currently sctx is memory bleed */ 5023 ierr = PetscMalloc(sizeof(SNESMatlabContext),&sctx);CHKERRQ(ierr); 5024 ierr = PetscStrallocpy(func,&sctx->funcname);CHKERRQ(ierr); 5025 /* 5026 This should work, but it doesn't 5027 sctx->ctx = ctx; 5028 mexMakeArrayPersistent(sctx->ctx); 5029 */ 5030 sctx->ctx = mxDuplicateArray(ctx); 5031 ierr = SNESMonitorSet(snes,SNESMonitor_Matlab,sctx,PETSC_NULL);CHKERRQ(ierr); 5032 PetscFunctionReturn(0); 5033 } 5034 5035 #endif 5036