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 ierr = (*newdm->ops->load)(newdm,viewer);CHKERRQ(ierr); 176 ierr = SNESGetKSP(newdm,&ksp);CHKERRQ(ierr); 177 ierr = KSPLoad(ksp,viewer);CHKERRQ(ierr); 178 PetscFunctionReturn(0); 179 } 180 181 #undef __FUNCT__ 182 #define __FUNCT__ "SNESView" 183 /*@C 184 SNESView - Prints the SNES data structure. 185 186 Collective on SNES 187 188 Input Parameters: 189 + SNES - the SNES context 190 - viewer - visualization context 191 192 Options Database Key: 193 . -snes_view - Calls SNESView() at end of SNESSolve() 194 195 Notes: 196 The available visualization contexts include 197 + PETSC_VIEWER_STDOUT_SELF - standard output (default) 198 - PETSC_VIEWER_STDOUT_WORLD - synchronized standard 199 output where only the first processor opens 200 the file. All other processors send their 201 data to the first processor to print. 202 203 The user can open an alternative visualization context with 204 PetscViewerASCIIOpen() - output to a specified file. 205 206 Level: beginner 207 208 .keywords: SNES, view 209 210 .seealso: PetscViewerASCIIOpen() 211 @*/ 212 PetscErrorCode SNESView(SNES snes,PetscViewer viewer) 213 { 214 SNESKSPEW *kctx; 215 PetscErrorCode ierr; 216 KSP ksp; 217 SNESLineSearch linesearch; 218 PetscBool iascii,isstring,isbinary; 219 220 PetscFunctionBegin; 221 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 222 if (!viewer) { 223 ierr = PetscViewerASCIIGetStdout(((PetscObject)snes)->comm,&viewer);CHKERRQ(ierr); 224 } 225 PetscValidHeaderSpecific(viewer,PETSC_VIEWER_CLASSID,2); 226 PetscCheckSameComm(snes,1,viewer,2); 227 228 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 229 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSTRING,&isstring);CHKERRQ(ierr); 230 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr); 231 if (iascii) { 232 ierr = PetscObjectPrintClassNamePrefixType((PetscObject)snes,viewer,"SNES Object");CHKERRQ(ierr); 233 if (snes->ops->view) { 234 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 235 ierr = (*snes->ops->view)(snes,viewer);CHKERRQ(ierr); 236 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 237 } 238 ierr = PetscViewerASCIIPrintf(viewer," maximum iterations=%D, maximum function evaluations=%D\n",snes->max_its,snes->max_funcs);CHKERRQ(ierr); 239 ierr = PetscViewerASCIIPrintf(viewer," tolerances: relative=%G, absolute=%G, solution=%G\n", 240 snes->rtol,snes->abstol,snes->stol);CHKERRQ(ierr); 241 ierr = PetscViewerASCIIPrintf(viewer," total number of linear solver iterations=%D\n",snes->linear_its);CHKERRQ(ierr); 242 ierr = PetscViewerASCIIPrintf(viewer," total number of function evaluations=%D\n",snes->nfuncs);CHKERRQ(ierr); 243 if (snes->gridsequence) { 244 ierr = PetscViewerASCIIPrintf(viewer," total number of grid sequence refinements=%D\n",snes->gridsequence);CHKERRQ(ierr); 245 } 246 if (snes->ksp_ewconv) { 247 kctx = (SNESKSPEW *)snes->kspconvctx; 248 if (kctx) { 249 ierr = PetscViewerASCIIPrintf(viewer," Eisenstat-Walker computation of KSP relative tolerance (version %D)\n",kctx->version);CHKERRQ(ierr); 250 ierr = PetscViewerASCIIPrintf(viewer," rtol_0=%G, rtol_max=%G, threshold=%G\n",kctx->rtol_0,kctx->rtol_max,kctx->threshold);CHKERRQ(ierr); 251 ierr = PetscViewerASCIIPrintf(viewer," gamma=%G, alpha=%G, alpha2=%G\n",kctx->gamma,kctx->alpha,kctx->alpha2);CHKERRQ(ierr); 252 } 253 } 254 if (snes->lagpreconditioner == -1) { 255 ierr = PetscViewerASCIIPrintf(viewer," Preconditioned is never rebuilt\n");CHKERRQ(ierr); 256 } else if (snes->lagpreconditioner > 1) { 257 ierr = PetscViewerASCIIPrintf(viewer," Preconditioned is rebuilt every %D new Jacobians\n",snes->lagpreconditioner);CHKERRQ(ierr); 258 } 259 if (snes->lagjacobian == -1) { 260 ierr = PetscViewerASCIIPrintf(viewer," Jacobian is never rebuilt\n");CHKERRQ(ierr); 261 } else if (snes->lagjacobian > 1) { 262 ierr = PetscViewerASCIIPrintf(viewer," Jacobian is rebuilt every %D SNES iterations\n",snes->lagjacobian);CHKERRQ(ierr); 263 } 264 } else if (isstring) { 265 const char *type; 266 ierr = SNESGetType(snes,&type);CHKERRQ(ierr); 267 ierr = PetscViewerStringSPrintf(viewer," %-3.3s",type);CHKERRQ(ierr); 268 } else if (isbinary) { 269 PetscInt classid = SNES_FILE_CLASSID; 270 MPI_Comm comm; 271 PetscMPIInt rank; 272 char type[256]; 273 274 ierr = PetscObjectGetComm((PetscObject)snes,&comm);CHKERRQ(ierr); 275 ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr); 276 if (!rank) { 277 ierr = PetscViewerBinaryWrite(viewer,&classid,1,PETSC_INT,PETSC_FALSE);CHKERRQ(ierr); 278 ierr = PetscStrncpy(type,((PetscObject)snes)->type_name,256);CHKERRQ(ierr); 279 ierr = PetscViewerBinaryWrite(viewer,type,256,PETSC_CHAR,PETSC_FALSE);CHKERRQ(ierr); 280 } 281 if (snes->ops->view) { 282 ierr = (*snes->ops->view)(snes,viewer);CHKERRQ(ierr); 283 } 284 } 285 if (snes->pc && snes->usespc) { 286 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 287 ierr = SNESView(snes->pc, viewer);CHKERRQ(ierr); 288 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 289 } 290 if (snes->usesksp) { 291 ierr = SNESGetKSP(snes,&ksp);CHKERRQ(ierr); 292 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 293 ierr = KSPView(ksp,viewer);CHKERRQ(ierr); 294 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 295 } 296 if (snes->linesearch) { 297 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 298 ierr = SNESGetSNESLineSearch(snes, &linesearch);CHKERRQ(ierr); 299 ierr = SNESLineSearchView(linesearch, viewer);CHKERRQ(ierr); 300 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 301 } 302 PetscFunctionReturn(0); 303 } 304 305 /* 306 We retain a list of functions that also take SNES command 307 line options. These are called at the end SNESSetFromOptions() 308 */ 309 #define MAXSETFROMOPTIONS 5 310 static PetscInt numberofsetfromoptions; 311 static PetscErrorCode (*othersetfromoptions[MAXSETFROMOPTIONS])(SNES); 312 313 #undef __FUNCT__ 314 #define __FUNCT__ "SNESAddOptionsChecker" 315 /*@C 316 SNESAddOptionsChecker - Adds an additional function to check for SNES options. 317 318 Not Collective 319 320 Input Parameter: 321 . snescheck - function that checks for options 322 323 Level: developer 324 325 .seealso: SNESSetFromOptions() 326 @*/ 327 PetscErrorCode SNESAddOptionsChecker(PetscErrorCode (*snescheck)(SNES)) 328 { 329 PetscFunctionBegin; 330 if (numberofsetfromoptions >= MAXSETFROMOPTIONS) { 331 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE, "Too many options checkers, only %D allowed", MAXSETFROMOPTIONS); 332 } 333 othersetfromoptions[numberofsetfromoptions++] = snescheck; 334 PetscFunctionReturn(0); 335 } 336 337 extern PetscErrorCode SNESDefaultMatrixFreeCreate2(SNES,Vec,Mat*); 338 339 #undef __FUNCT__ 340 #define __FUNCT__ "SNESSetUpMatrixFree_Private" 341 static PetscErrorCode SNESSetUpMatrixFree_Private(SNES snes, PetscBool hasOperator, PetscInt version) 342 { 343 Mat J; 344 KSP ksp; 345 PC pc; 346 PetscBool match; 347 PetscErrorCode ierr; 348 349 PetscFunctionBegin; 350 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 351 352 if (!snes->vec_func && (snes->jacobian || snes->jacobian_pre)) { 353 Mat A = snes->jacobian, B = snes->jacobian_pre; 354 ierr = MatGetVecs(A ? A : B, PETSC_NULL,&snes->vec_func);CHKERRQ(ierr); 355 } 356 357 if (version == 1) { 358 ierr = MatCreateSNESMF(snes,&J);CHKERRQ(ierr); 359 ierr = MatMFFDSetOptionsPrefix(J,((PetscObject)snes)->prefix);CHKERRQ(ierr); 360 ierr = MatSetFromOptions(J);CHKERRQ(ierr); 361 } else if (version == 2) { 362 if (!snes->vec_func) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"SNESSetFunction() must be called first"); 363 #if !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_REAL_SINGLE) && !defined(PETSC_USE_REAL___FLOAT128) 364 ierr = SNESDefaultMatrixFreeCreate2(snes,snes->vec_func,&J);CHKERRQ(ierr); 365 #else 366 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP, "matrix-free operator rutines (version 2)"); 367 #endif 368 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE, "matrix-free operator rutines, only version 1 and 2"); 369 370 ierr = PetscInfo1(snes,"Setting default matrix-free operator routines (version %D)\n", version);CHKERRQ(ierr); 371 if (hasOperator) { 372 373 /* This version replaces the user provided Jacobian matrix with a 374 matrix-free version but still employs the user-provided preconditioner matrix. */ 375 ierr = SNESSetJacobian(snes,J,0,0,0);CHKERRQ(ierr); 376 } else { 377 /* This version replaces both the user-provided Jacobian and the user- 378 provided preconditioner Jacobian with the default matrix free version. */ 379 380 ierr = SNESSetJacobian(snes,J,J,0,0);CHKERRQ(ierr); 381 /* Force no preconditioner */ 382 ierr = SNESGetKSP(snes,&ksp);CHKERRQ(ierr); 383 ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr); 384 ierr = PetscObjectTypeCompare((PetscObject)pc,PCSHELL,&match);CHKERRQ(ierr); 385 if (!match) { 386 ierr = PetscInfo(snes,"Setting default matrix-free preconditioner routines\nThat is no preconditioner is being used\n");CHKERRQ(ierr); 387 ierr = PCSetType(pc,PCNONE);CHKERRQ(ierr); 388 } 389 } 390 ierr = MatDestroy(&J);CHKERRQ(ierr); 391 PetscFunctionReturn(0); 392 } 393 394 #undef __FUNCT__ 395 #define __FUNCT__ "DMRestrictHook_SNESVecSol" 396 static PetscErrorCode DMRestrictHook_SNESVecSol(DM dmfine,Mat Restrict,Vec Rscale,Mat Inject,DM dmcoarse,void *ctx) 397 { 398 SNES snes = (SNES)ctx; 399 PetscErrorCode ierr; 400 Vec Xfine,Xfine_named = PETSC_NULL,Xcoarse; 401 402 PetscFunctionBegin; 403 if (PetscLogPrintInfo) { 404 PetscInt finelevel,coarselevel,fineclevel,coarseclevel; 405 ierr = DMGetRefineLevel(dmfine,&finelevel);CHKERRQ(ierr); 406 ierr = DMGetCoarsenLevel(dmfine,&fineclevel);CHKERRQ(ierr); 407 ierr = DMGetRefineLevel(dmcoarse,&coarselevel);CHKERRQ(ierr); 408 ierr = DMGetCoarsenLevel(dmcoarse,&coarseclevel);CHKERRQ(ierr); 409 ierr = PetscInfo4(dmfine,"Restricting SNES solution vector from level %D-%D to level %D-%D\n",finelevel,fineclevel,coarselevel,coarseclevel);CHKERRQ(ierr); 410 } 411 if (dmfine == snes->dm) Xfine = snes->vec_sol; 412 else { 413 ierr = DMGetNamedGlobalVector(dmfine,"SNESVecSol",&Xfine_named);CHKERRQ(ierr); 414 Xfine = Xfine_named; 415 } 416 ierr = DMGetNamedGlobalVector(dmcoarse,"SNESVecSol",&Xcoarse);CHKERRQ(ierr); 417 ierr = MatRestrict(Restrict,Xfine,Xcoarse);CHKERRQ(ierr); 418 ierr = VecPointwiseMult(Xcoarse,Xcoarse,Rscale);CHKERRQ(ierr); 419 ierr = DMRestoreNamedGlobalVector(dmcoarse,"SNESVecSol",&Xcoarse);CHKERRQ(ierr); 420 if (Xfine_named) {ierr = DMRestoreNamedGlobalVector(dmfine,"SNESVecSol",&Xfine_named);CHKERRQ(ierr);} 421 PetscFunctionReturn(0); 422 } 423 424 #undef __FUNCT__ 425 #define __FUNCT__ "DMCoarsenHook_SNESVecSol" 426 static PetscErrorCode DMCoarsenHook_SNESVecSol(DM dm,DM dmc,void *ctx) 427 { 428 PetscErrorCode ierr; 429 430 PetscFunctionBegin; 431 ierr = DMCoarsenHookAdd(dmc,DMCoarsenHook_SNESVecSol,DMRestrictHook_SNESVecSol,ctx);CHKERRQ(ierr); 432 PetscFunctionReturn(0); 433 } 434 435 #undef __FUNCT__ 436 #define __FUNCT__ "KSPComputeOperators_SNES" 437 /* This may be called to rediscretize the operator on levels of linear multigrid. The DM shuffle is so the user can 438 * safely call SNESGetDM() in their residual evaluation routine. */ 439 static PetscErrorCode KSPComputeOperators_SNES(KSP ksp,Mat A,Mat B,MatStructure *mstruct,void *ctx) 440 { 441 SNES snes = (SNES)ctx; 442 PetscErrorCode ierr; 443 Mat Asave = A,Bsave = B; 444 Vec X,Xnamed = PETSC_NULL; 445 DM dmsave; 446 void *ctxsave; 447 PetscErrorCode (*jac)(SNES,Vec,Mat*,Mat*,MatStructure*,void*); 448 449 PetscFunctionBegin; 450 dmsave = snes->dm; 451 ierr = KSPGetDM(ksp,&snes->dm);CHKERRQ(ierr); 452 if (dmsave == snes->dm) X = snes->vec_sol; /* We are on the finest level */ 453 else { /* We are on a coarser level, this vec was initialized using a DM restrict hook */ 454 ierr = DMGetNamedGlobalVector(snes->dm,"SNESVecSol",&Xnamed);CHKERRQ(ierr); 455 X = Xnamed; 456 ierr = SNESGetJacobian(snes,PETSC_NULL,PETSC_NULL,&jac,&ctxsave);CHKERRQ(ierr); 457 /* If the DM's don't match up, the MatFDColoring context needed for the jacobian won't match up either -- fixit. */ 458 if (jac == SNESDefaultComputeJacobianColor) { 459 ierr = SNESSetJacobian(snes,PETSC_NULL,PETSC_NULL,SNESDefaultComputeJacobianColor,0);CHKERRQ(ierr); 460 } 461 } 462 /* put the previous context back */ 463 464 ierr = SNESComputeJacobian(snes,X,&A,&B,mstruct);CHKERRQ(ierr); 465 if (snes->dm != dmsave && jac == SNESDefaultComputeJacobianColor) { 466 ierr = SNESSetJacobian(snes,PETSC_NULL,PETSC_NULL,jac,ctxsave);CHKERRQ(ierr); 467 } 468 469 if (A != Asave || B != Bsave) SETERRQ(((PetscObject)snes)->comm,PETSC_ERR_SUP,"No support for changing matrices at this time"); 470 if (Xnamed) { 471 ierr = DMRestoreNamedGlobalVector(snes->dm,"SNESVecSol",&Xnamed);CHKERRQ(ierr); 472 } 473 snes->dm = dmsave; 474 PetscFunctionReturn(0); 475 } 476 477 #undef __FUNCT__ 478 #define __FUNCT__ "SNESSetUpMatrices" 479 /*@ 480 SNESSetUpMatrices - ensures that matrices are available for SNES, to be called by SNESSetUp_XXX() 481 482 Collective 483 484 Input Arguments: 485 . snes - snes to configure 486 487 Level: developer 488 489 .seealso: SNESSetUp() 490 @*/ 491 PetscErrorCode SNESSetUpMatrices(SNES snes) 492 { 493 PetscErrorCode ierr; 494 DM dm; 495 SNESDM sdm; 496 497 PetscFunctionBegin; 498 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 499 ierr = DMSNESGetContext(dm,&sdm);CHKERRQ(ierr); 500 if (!sdm->computejacobian) { 501 SETERRQ(((PetscObject)snes)->comm,PETSC_ERR_PLIB,"SNESDM not properly configured"); 502 } else if (!snes->jacobian && snes->mf) { 503 Mat J; 504 void *functx; 505 ierr = MatCreateSNESMF(snes,&J);CHKERRQ(ierr); 506 ierr = MatMFFDSetOptionsPrefix(J,((PetscObject)snes)->prefix);CHKERRQ(ierr); 507 ierr = MatSetFromOptions(J);CHKERRQ(ierr); 508 ierr = SNESGetFunction(snes,PETSC_NULL,PETSC_NULL,&functx);CHKERRQ(ierr); 509 ierr = SNESSetJacobian(snes,J,J,0,0);CHKERRQ(ierr); 510 ierr = MatDestroy(&J);CHKERRQ(ierr); 511 } else if (snes->mf_operator && !snes->jacobian_pre && !snes->jacobian) { 512 Mat J,B; 513 ierr = MatCreateSNESMF(snes,&J);CHKERRQ(ierr); 514 ierr = MatMFFDSetOptionsPrefix(J,((PetscObject)snes)->prefix);CHKERRQ(ierr); 515 ierr = MatSetFromOptions(J);CHKERRQ(ierr); 516 ierr = DMCreateMatrix(snes->dm,MATAIJ,&B);CHKERRQ(ierr); 517 /* sdm->computejacobian was already set to reach here */ 518 ierr = SNESSetJacobian(snes,J,B,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr); 519 ierr = MatDestroy(&J);CHKERRQ(ierr); 520 ierr = MatDestroy(&B);CHKERRQ(ierr); 521 } else if (!snes->jacobian_pre) { 522 Mat J,B; 523 J = snes->jacobian; 524 ierr = DMCreateMatrix(snes->dm,MATAIJ,&B);CHKERRQ(ierr); 525 ierr = SNESSetJacobian(snes,J?J:B,B,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr); 526 ierr = MatDestroy(&B);CHKERRQ(ierr); 527 } 528 { 529 KSP ksp; 530 ierr = SNESGetKSP(snes,&ksp);CHKERRQ(ierr); 531 ierr = KSPSetComputeOperators(ksp,KSPComputeOperators_SNES,snes);CHKERRQ(ierr); 532 ierr = DMCoarsenHookAdd(snes->dm,DMCoarsenHook_SNESVecSol,DMRestrictHook_SNESVecSol,snes);CHKERRQ(ierr); 533 } 534 PetscFunctionReturn(0); 535 } 536 537 #undef __FUNCT__ 538 #define __FUNCT__ "SNESSetFromOptions" 539 /*@ 540 SNESSetFromOptions - Sets various SNES and KSP parameters from user options. 541 542 Collective on SNES 543 544 Input Parameter: 545 . snes - the SNES context 546 547 Options Database Keys: 548 + -snes_type <type> - ls, tr, ngmres, ncg, richardson, qn, vi, fas 549 . -snes_stol - convergence tolerance in terms of the norm 550 of the change in the solution between steps 551 . -snes_atol <abstol> - absolute tolerance of residual norm 552 . -snes_rtol <rtol> - relative decrease in tolerance norm from initial 553 . -snes_max_it <max_it> - maximum number of iterations 554 . -snes_max_funcs <max_funcs> - maximum number of function evaluations 555 . -snes_max_fail <max_fail> - maximum number of line search failures allowed before stopping, default is none 556 . -snes_max_linear_solve_fail - number of linear solver failures before SNESSolve() stops 557 . -snes_lag_preconditioner <lag> - how often preconditioner is rebuilt (use -1 to never rebuild) 558 . -snes_lag_jacobian <lag> - how often Jacobian is rebuilt (use -1 to never rebuild) 559 . -snes_trtol <trtol> - trust region tolerance 560 . -snes_no_convergence_test - skip convergence test in nonlinear 561 solver; hence iterations will continue until max_it 562 or some other criterion is reached. Saves expense 563 of convergence test 564 . -snes_monitor <optional filename> - prints residual norm at each iteration. if no 565 filename given prints to stdout 566 . -snes_monitor_solution - plots solution at each iteration 567 . -snes_monitor_residual - plots residual (not its norm) at each iteration 568 . -snes_monitor_solution_update - plots update to solution at each iteration 569 . -snes_monitor_lg_residualnorm - plots residual norm at each iteration 570 . -snes_monitor_lg_range - plots residual norm at each iteration 571 . -snes_fd - use finite differences to compute Jacobian; very slow, only for testing 572 . -snes_fd_color - use finite differences with coloring to compute Jacobian 573 . -snes_mf_ksp_monitor - if using matrix-free multiply then print h at each KSP iteration 574 - -snes_converged_reason - print the reason for convergence/divergence after each solve 575 576 Options Database for Eisenstat-Walker method: 577 + -snes_ksp_ew - use Eisenstat-Walker method for determining linear system convergence 578 . -snes_ksp_ew_version ver - version of Eisenstat-Walker method 579 . -snes_ksp_ew_rtol0 <rtol0> - Sets rtol0 580 . -snes_ksp_ew_rtolmax <rtolmax> - Sets rtolmax 581 . -snes_ksp_ew_gamma <gamma> - Sets gamma 582 . -snes_ksp_ew_alpha <alpha> - Sets alpha 583 . -snes_ksp_ew_alpha2 <alpha2> - Sets alpha2 584 - -snes_ksp_ew_threshold <threshold> - Sets threshold 585 586 Notes: 587 To see all options, run your program with the -help option or consult 588 the <A href="../../docs/manual.pdf#nameddest=ch_snes">SNES chapter of the users manual</A>. 589 590 Level: beginner 591 592 .keywords: SNES, nonlinear, set, options, database 593 594 .seealso: SNESSetOptionsPrefix() 595 @*/ 596 PetscErrorCode SNESSetFromOptions(SNES snes) 597 { 598 PetscBool flg,pcset; 599 PetscInt i,indx,lag,grids; 600 MatStructure matflag; 601 const char *deft = SNESLS; 602 const char *convtests[] = {"default","skip"}; 603 SNESKSPEW *kctx = NULL; 604 char type[256], monfilename[PETSC_MAX_PATH_LEN]; 605 PetscViewer monviewer; 606 PetscErrorCode ierr; 607 PCSide pcside; 608 const char *optionsprefix; 609 610 PetscFunctionBegin; 611 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 612 613 if (!SNESRegisterAllCalled) {ierr = SNESRegisterAll(PETSC_NULL);CHKERRQ(ierr);} 614 ierr = PetscObjectOptionsBegin((PetscObject)snes);CHKERRQ(ierr); 615 if (((PetscObject)snes)->type_name) { deft = ((PetscObject)snes)->type_name; } 616 ierr = PetscOptionsList("-snes_type","Nonlinear solver method","SNESSetType",SNESList,deft,type,256,&flg);CHKERRQ(ierr); 617 if (flg) { 618 ierr = SNESSetType(snes,type);CHKERRQ(ierr); 619 } else if (!((PetscObject)snes)->type_name) { 620 ierr = SNESSetType(snes,deft);CHKERRQ(ierr); 621 } 622 /* not used here, but called so will go into help messaage */ 623 ierr = PetscOptionsName("-snes_view","Print detailed information on solver used","SNESView",0);CHKERRQ(ierr); 624 625 ierr = PetscOptionsReal("-snes_stol","Stop if step length less than","SNESSetTolerances",snes->stol,&snes->stol,0);CHKERRQ(ierr); 626 ierr = PetscOptionsReal("-snes_atol","Stop if function norm less than","SNESSetTolerances",snes->abstol,&snes->abstol,0);CHKERRQ(ierr); 627 628 ierr = PetscOptionsReal("-snes_rtol","Stop if decrease in function norm less than","SNESSetTolerances",snes->rtol,&snes->rtol,0);CHKERRQ(ierr); 629 ierr = PetscOptionsInt("-snes_max_it","Maximum iterations","SNESSetTolerances",snes->max_its,&snes->max_its,PETSC_NULL);CHKERRQ(ierr); 630 ierr = PetscOptionsInt("-snes_max_funcs","Maximum function evaluations","SNESSetTolerances",snes->max_funcs,&snes->max_funcs,PETSC_NULL);CHKERRQ(ierr); 631 ierr = PetscOptionsInt("-snes_max_fail","Maximum nonlinear step failures","SNESSetMaxNonlinearStepFailures",snes->maxFailures,&snes->maxFailures,PETSC_NULL);CHKERRQ(ierr); 632 ierr = PetscOptionsInt("-snes_max_linear_solve_fail","Maximum failures in linear solves allowed","SNESSetMaxLinearSolveFailures",snes->maxLinearSolveFailures,&snes->maxLinearSolveFailures,PETSC_NULL);CHKERRQ(ierr); 633 ierr = PetscOptionsBool("-snes_error_if_not_converged","Generate error if solver does not converge","SNESSetErrorIfNotConverged",snes->errorifnotconverged,&snes->errorifnotconverged,PETSC_NULL);CHKERRQ(ierr); 634 635 ierr = PetscOptionsInt("-snes_lag_preconditioner","How often to rebuild preconditioner","SNESSetLagPreconditioner",snes->lagpreconditioner,&lag,&flg);CHKERRQ(ierr); 636 if (flg) { 637 ierr = SNESSetLagPreconditioner(snes,lag);CHKERRQ(ierr); 638 } 639 ierr = PetscOptionsInt("-snes_lag_jacobian","How often to rebuild Jacobian","SNESSetLagJacobian",snes->lagjacobian,&lag,&flg);CHKERRQ(ierr); 640 if (flg) { 641 ierr = SNESSetLagJacobian(snes,lag);CHKERRQ(ierr); 642 } 643 ierr = PetscOptionsInt("-snes_grid_sequence","Use grid sequencing to generate initial guess","SNESSetGridSequence",snes->gridsequence,&grids,&flg);CHKERRQ(ierr); 644 if (flg) { 645 ierr = SNESSetGridSequence(snes,grids);CHKERRQ(ierr); 646 } 647 648 ierr = PetscOptionsEList("-snes_convergence_test","Convergence test","SNESSetConvergenceTest",convtests,2,"default",&indx,&flg);CHKERRQ(ierr); 649 if (flg) { 650 switch (indx) { 651 case 0: ierr = SNESSetConvergenceTest(snes,SNESDefaultConverged,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr); break; 652 case 1: ierr = SNESSetConvergenceTest(snes,SNESSkipConverged,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr); break; 653 } 654 } 655 656 ierr = PetscOptionsBool("-snes_converged_reason","Print reason for converged or diverged","SNESSolve",snes->printreason,&snes->printreason,PETSC_NULL);CHKERRQ(ierr); 657 658 ierr = PetscOptionsEList("-snes_norm_type","SNES Norm type","SNESSetNormType",SNESNormTypes,5,"function",&indx,&flg);CHKERRQ(ierr); 659 if (flg) { ierr = SNESSetNormType(snes,(SNESNormType)indx);CHKERRQ(ierr); } 660 661 kctx = (SNESKSPEW *)snes->kspconvctx; 662 663 ierr = PetscOptionsBool("-snes_ksp_ew","Use Eisentat-Walker linear system convergence test","SNESKSPSetUseEW",snes->ksp_ewconv,&snes->ksp_ewconv,PETSC_NULL);CHKERRQ(ierr); 664 665 ierr = PetscOptionsInt("-snes_ksp_ew_version","Version 1, 2 or 3","SNESKSPSetParametersEW",kctx->version,&kctx->version,0);CHKERRQ(ierr); 666 ierr = PetscOptionsReal("-snes_ksp_ew_rtol0","0 <= rtol0 < 1","SNESKSPSetParametersEW",kctx->rtol_0,&kctx->rtol_0,0);CHKERRQ(ierr); 667 ierr = PetscOptionsReal("-snes_ksp_ew_rtolmax","0 <= rtolmax < 1","SNESKSPSetParametersEW",kctx->rtol_max,&kctx->rtol_max,0);CHKERRQ(ierr); 668 ierr = PetscOptionsReal("-snes_ksp_ew_gamma","0 <= gamma <= 1","SNESKSPSetParametersEW",kctx->gamma,&kctx->gamma,0);CHKERRQ(ierr); 669 ierr = PetscOptionsReal("-snes_ksp_ew_alpha","1 < alpha <= 2","SNESKSPSetParametersEW",kctx->alpha,&kctx->alpha,0);CHKERRQ(ierr); 670 ierr = PetscOptionsReal("-snes_ksp_ew_alpha2","alpha2","SNESKSPSetParametersEW",kctx->alpha2,&kctx->alpha2,0);CHKERRQ(ierr); 671 ierr = PetscOptionsReal("-snes_ksp_ew_threshold","0 < threshold < 1","SNESKSPSetParametersEW",kctx->threshold,&kctx->threshold,0);CHKERRQ(ierr); 672 673 flg = PETSC_FALSE; 674 ierr = PetscOptionsBool("-snes_monitor_cancel","Remove all monitors","SNESMonitorCancel",flg,&flg,PETSC_NULL);CHKERRQ(ierr); 675 if (flg) {ierr = SNESMonitorCancel(snes);CHKERRQ(ierr);} 676 677 ierr = PetscOptionsString("-snes_monitor","Monitor norm of function","SNESMonitorSet","stdout",monfilename,PETSC_MAX_PATH_LEN,&flg);CHKERRQ(ierr); 678 if (flg) { 679 ierr = PetscViewerASCIIOpen(((PetscObject)snes)->comm,monfilename,&monviewer);CHKERRQ(ierr); 680 ierr = SNESMonitorSet(snes,SNESMonitorDefault,monviewer,(PetscErrorCode (*)(void**))PetscViewerDestroy);CHKERRQ(ierr); 681 } 682 683 ierr = PetscOptionsString("-snes_monitor_range","Monitor range of elements of function","SNESMonitorSet","stdout",monfilename,PETSC_MAX_PATH_LEN,&flg);CHKERRQ(ierr); 684 if (flg) { 685 ierr = SNESMonitorSet(snes,SNESMonitorRange,0,0);CHKERRQ(ierr); 686 } 687 688 ierr = PetscOptionsString("-snes_ratiomonitor","Monitor ratios of norms of function","SNESMonitorSetRatio","stdout",monfilename,PETSC_MAX_PATH_LEN,&flg);CHKERRQ(ierr); 689 if (flg) { 690 ierr = PetscViewerASCIIOpen(((PetscObject)snes)->comm,monfilename,&monviewer);CHKERRQ(ierr); 691 ierr = SNESMonitorSetRatio(snes,monviewer);CHKERRQ(ierr); 692 } 693 694 ierr = PetscOptionsString("-snes_monitor_short","Monitor norm of function (fewer digits)","SNESMonitorSet","stdout",monfilename,PETSC_MAX_PATH_LEN,&flg);CHKERRQ(ierr); 695 if (flg) { 696 ierr = PetscViewerASCIIOpen(((PetscObject)snes)->comm,monfilename,&monviewer);CHKERRQ(ierr); 697 ierr = SNESMonitorSet(snes,SNESMonitorDefaultShort,monviewer,(PetscErrorCode (*)(void**))PetscViewerDestroy);CHKERRQ(ierr); 698 } 699 700 ierr = PetscOptionsString("-snes_monitor_python","Use Python function","SNESMonitorSet",0,monfilename,PETSC_MAX_PATH_LEN,&flg);CHKERRQ(ierr); 701 if (flg) {ierr = PetscPythonMonitorSet((PetscObject)snes,monfilename);CHKERRQ(ierr);} 702 703 flg = PETSC_FALSE; 704 ierr = PetscOptionsBool("-snes_monitor_solution","Plot solution at each iteration","SNESMonitorSolution",flg,&flg,PETSC_NULL);CHKERRQ(ierr); 705 if (flg) {ierr = SNESMonitorSet(snes,SNESMonitorSolution,0,0);CHKERRQ(ierr);} 706 flg = PETSC_FALSE; 707 ierr = PetscOptionsBool("-snes_monitor_solution_update","Plot correction at each iteration","SNESMonitorSolutionUpdate",flg,&flg,PETSC_NULL);CHKERRQ(ierr); 708 if (flg) {ierr = SNESMonitorSet(snes,SNESMonitorSolutionUpdate,0,0);CHKERRQ(ierr);} 709 flg = PETSC_FALSE; 710 ierr = PetscOptionsBool("-snes_monitor_residual","Plot residual at each iteration","SNESMonitorResidual",flg,&flg,PETSC_NULL);CHKERRQ(ierr); 711 if (flg) {ierr = SNESMonitorSet(snes,SNESMonitorResidual,0,0);CHKERRQ(ierr);} 712 flg = PETSC_FALSE; 713 ierr = PetscOptionsBool("-snes_monitor_lg_residualnorm","Plot function norm at each iteration","SNESMonitorLGResidualNorm",flg,&flg,PETSC_NULL);CHKERRQ(ierr); 714 if (flg) { 715 PetscDrawLG ctx; 716 717 ierr = SNESMonitorLGCreate(0,0,PETSC_DECIDE,PETSC_DECIDE,300,300,&ctx);CHKERRQ(ierr); 718 ierr = SNESMonitorSet(snes,SNESMonitorLGResidualNorm,ctx,(PetscErrorCode (*)(void**))SNESMonitorLGDestroy);CHKERRQ(ierr); 719 } 720 flg = PETSC_FALSE; 721 ierr = PetscOptionsBool("-snes_monitor_lg_range","Plot function range at each iteration","SNESMonitorLGRange",flg,&flg,PETSC_NULL);CHKERRQ(ierr); 722 if (flg) { 723 PetscViewer ctx; 724 725 ierr = PetscViewerDrawOpen(((PetscObject)snes)->comm,0,0,PETSC_DECIDE,PETSC_DECIDE,300,300,&ctx);CHKERRQ(ierr); 726 ierr = SNESMonitorSet(snes,SNESMonitorLGRange,ctx,(PetscErrorCode (*)(void**))PetscViewerDestroy);CHKERRQ(ierr); 727 } 728 729 flg = PETSC_FALSE; 730 ierr = PetscOptionsBool("-snes_monitor_jacupdate_spectrum","Print the change in the spectrum of the Jacobian","SNESMonitorJacUpdateSpectrum",flg,&flg,PETSC_NULL);CHKERRQ(ierr); 731 if (flg) {ierr = SNESMonitorSet(snes,SNESMonitorJacUpdateSpectrum,0,0);CHKERRQ(ierr);} 732 733 flg = PETSC_FALSE; 734 ierr = PetscOptionsBool("-snes_fd","Use finite differences (slow) to compute Jacobian","SNESDefaultComputeJacobian",flg,&flg,PETSC_NULL);CHKERRQ(ierr); 735 if (flg) { 736 void *functx; 737 ierr = SNESGetFunction(snes,PETSC_NULL,PETSC_NULL,&functx);CHKERRQ(ierr); 738 ierr = SNESSetJacobian(snes,snes->jacobian,snes->jacobian_pre,SNESDefaultComputeJacobian,functx);CHKERRQ(ierr); 739 ierr = PetscInfo(snes,"Setting default finite difference Jacobian matrix\n");CHKERRQ(ierr); 740 } 741 742 flg = PETSC_FALSE; 743 ierr = PetscOptionsBool("-snes_fd_function","Use finite differences (slow) to compute function from user objective","SNESDefaultObjectiveComputeFunctionFD",flg,&flg,PETSC_NULL);CHKERRQ(ierr); 744 if (flg) { 745 ierr = SNESSetFunction(snes,PETSC_NULL,SNESDefaultObjectiveComputeFunctionFD,PETSC_NULL);CHKERRQ(ierr); 746 } 747 748 flg = PETSC_FALSE; 749 ierr = PetscOptionsBool("-snes_fd_color","Use finite differences with coloring to compute Jacobian","SNESDefaultComputeJacobianColor",flg,&flg,PETSC_NULL);CHKERRQ(ierr); 750 if (flg) { 751 void *functx; 752 ierr = SNESGetFunction(snes,PETSC_NULL,PETSC_NULL,&functx);CHKERRQ(ierr); 753 ierr = SNESSetJacobian(snes,snes->jacobian,snes->jacobian_pre,SNESDefaultComputeJacobianColor,0);CHKERRQ(ierr); 754 ierr = PetscInfo(snes,"Setting default finite difference coloring Jacobian matrix\n");CHKERRQ(ierr); 755 } 756 757 flg = PETSC_FALSE; 758 ierr = PetscOptionsBool("-snes_mf_operator","Use a Matrix-Free Jacobian with user-provided preconditioner matrix","MatCreateSNESMF",PETSC_FALSE,&snes->mf_operator,&flg);CHKERRQ(ierr); 759 if (flg && snes->mf_operator) { 760 snes->mf_operator = PETSC_TRUE; 761 snes->mf = PETSC_TRUE; 762 } 763 flg = PETSC_FALSE; 764 ierr = PetscOptionsBool("-snes_mf","Use a Matrix-Free Jacobian with no preconditioner matrix","MatCreateSNESMF",PETSC_FALSE,&snes->mf,&flg);CHKERRQ(ierr); 765 if (!flg && snes->mf_operator) snes->mf = PETSC_TRUE; 766 ierr = PetscOptionsInt("-snes_mf_version","Matrix-Free routines version 1 or 2","None",snes->mf_version,&snes->mf_version,0);CHKERRQ(ierr); 767 768 flg = PETSC_FALSE; 769 ierr = SNESGetPCSide(snes,&pcside); 770 ierr = PetscOptionsEnum("-snes_npc_side","SNES nonlinear preconditioner side","SNESSetPCSide",PCSides,(PetscEnum)pcside,(PetscEnum*)&pcside,&flg);CHKERRQ(ierr); 771 if (flg) {ierr = SNESSetPCSide(snes,pcside);CHKERRQ(ierr);} 772 773 /* GS Options */ 774 ierr = PetscOptionsInt("-snes_gs_sweeps","Number of sweeps of GS to apply","SNESComputeGS",snes->gssweeps,&snes->gssweeps,PETSC_NULL);CHKERRQ(ierr); 775 776 for (i = 0; i < numberofsetfromoptions; i++) { 777 ierr = (*othersetfromoptions[i])(snes);CHKERRQ(ierr); 778 } 779 780 if (snes->ops->setfromoptions) { 781 ierr = (*snes->ops->setfromoptions)(snes);CHKERRQ(ierr); 782 } 783 784 /* process any options handlers added with PetscObjectAddOptionsHandler() */ 785 ierr = PetscObjectProcessOptionsHandlers((PetscObject)snes);CHKERRQ(ierr); 786 ierr = PetscOptionsEnd();CHKERRQ(ierr); 787 788 if (!snes->ksp) {ierr = SNESGetKSP(snes,&snes->ksp);CHKERRQ(ierr);} 789 ierr = KSPGetOperators(snes->ksp,PETSC_NULL,PETSC_NULL,&matflag); 790 ierr = KSPSetOperators(snes->ksp,snes->jacobian,snes->jacobian_pre,matflag);CHKERRQ(ierr); 791 ierr = KSPSetFromOptions(snes->ksp);CHKERRQ(ierr); 792 793 if (!snes->linesearch) { 794 ierr = SNESGetSNESLineSearch(snes, &snes->linesearch);CHKERRQ(ierr); 795 } 796 ierr = SNESLineSearchSetFromOptions(snes->linesearch);CHKERRQ(ierr); 797 798 /* if someone has set the SNES PC type, create it. */ 799 ierr = SNESGetOptionsPrefix(snes, &optionsprefix);CHKERRQ(ierr); 800 ierr = PetscOptionsHasName(optionsprefix, "-npc_snes_type", &pcset);CHKERRQ(ierr); 801 if (pcset && (!snes->pc)) { 802 ierr = SNESGetPC(snes, &snes->pc);CHKERRQ(ierr); 803 } 804 PetscFunctionReturn(0); 805 } 806 807 #undef __FUNCT__ 808 #define __FUNCT__ "SNESSetComputeApplicationContext" 809 /*@ 810 SNESSetComputeApplicationContext - Sets an optional function to compute a user-defined context for 811 the nonlinear solvers. 812 813 Logically Collective on SNES 814 815 Input Parameters: 816 + snes - the SNES context 817 . compute - function to compute the context 818 - destroy - function to destroy the context 819 820 Level: intermediate 821 822 .keywords: SNES, nonlinear, set, application, context 823 824 .seealso: SNESGetApplicationContext(), SNESSetComputeApplicationContext(), SNESGetApplicationContext() 825 @*/ 826 PetscErrorCode SNESSetComputeApplicationContext(SNES snes,PetscErrorCode (*compute)(SNES,void**),PetscErrorCode (*destroy)(void**)) 827 { 828 PetscFunctionBegin; 829 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 830 snes->ops->usercompute = compute; 831 snes->ops->userdestroy = destroy; 832 PetscFunctionReturn(0); 833 } 834 835 #undef __FUNCT__ 836 #define __FUNCT__ "SNESSetApplicationContext" 837 /*@ 838 SNESSetApplicationContext - Sets the optional user-defined context for 839 the nonlinear solvers. 840 841 Logically Collective on SNES 842 843 Input Parameters: 844 + snes - the SNES context 845 - usrP - optional user context 846 847 Level: intermediate 848 849 .keywords: SNES, nonlinear, set, application, context 850 851 .seealso: SNESGetApplicationContext() 852 @*/ 853 PetscErrorCode SNESSetApplicationContext(SNES snes,void *usrP) 854 { 855 PetscErrorCode ierr; 856 KSP ksp; 857 858 PetscFunctionBegin; 859 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 860 ierr = SNESGetKSP(snes,&ksp);CHKERRQ(ierr); 861 ierr = KSPSetApplicationContext(ksp,usrP);CHKERRQ(ierr); 862 snes->user = usrP; 863 PetscFunctionReturn(0); 864 } 865 866 #undef __FUNCT__ 867 #define __FUNCT__ "SNESGetApplicationContext" 868 /*@ 869 SNESGetApplicationContext - Gets the user-defined context for the 870 nonlinear solvers. 871 872 Not Collective 873 874 Input Parameter: 875 . snes - SNES context 876 877 Output Parameter: 878 . usrP - user context 879 880 Level: intermediate 881 882 .keywords: SNES, nonlinear, get, application, context 883 884 .seealso: SNESSetApplicationContext() 885 @*/ 886 PetscErrorCode SNESGetApplicationContext(SNES snes,void *usrP) 887 { 888 PetscFunctionBegin; 889 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 890 *(void**)usrP = snes->user; 891 PetscFunctionReturn(0); 892 } 893 894 #undef __FUNCT__ 895 #define __FUNCT__ "SNESGetIterationNumber" 896 /*@ 897 SNESGetIterationNumber - Gets the number of nonlinear iterations completed 898 at this time. 899 900 Not Collective 901 902 Input Parameter: 903 . snes - SNES context 904 905 Output Parameter: 906 . iter - iteration number 907 908 Notes: 909 For example, during the computation of iteration 2 this would return 1. 910 911 This is useful for using lagged Jacobians (where one does not recompute the 912 Jacobian at each SNES iteration). For example, the code 913 .vb 914 ierr = SNESGetIterationNumber(snes,&it); 915 if (!(it % 2)) { 916 [compute Jacobian here] 917 } 918 .ve 919 can be used in your ComputeJacobian() function to cause the Jacobian to be 920 recomputed every second SNES iteration. 921 922 Level: intermediate 923 924 .keywords: SNES, nonlinear, get, iteration, number, 925 926 .seealso: SNESGetFunctionNorm(), SNESGetLinearSolveIterations() 927 @*/ 928 PetscErrorCode SNESGetIterationNumber(SNES snes,PetscInt* iter) 929 { 930 PetscFunctionBegin; 931 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 932 PetscValidIntPointer(iter,2); 933 *iter = snes->iter; 934 PetscFunctionReturn(0); 935 } 936 937 #undef __FUNCT__ 938 #define __FUNCT__ "SNESSetIterationNumber" 939 /*@ 940 SNESSetIterationNumber - Sets the current iteration number. 941 942 Not Collective 943 944 Input Parameter: 945 . snes - SNES context 946 . iter - iteration number 947 948 Level: developer 949 950 .keywords: SNES, nonlinear, set, iteration, number, 951 952 .seealso: SNESGetFunctionNorm(), SNESGetLinearSolveIterations() 953 @*/ 954 PetscErrorCode SNESSetIterationNumber(SNES snes,PetscInt iter) 955 { 956 PetscErrorCode ierr; 957 958 PetscFunctionBegin; 959 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 960 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 961 snes->iter = iter; 962 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 963 PetscFunctionReturn(0); 964 } 965 966 #undef __FUNCT__ 967 #define __FUNCT__ "SNESGetFunctionNorm" 968 /*@ 969 SNESGetFunctionNorm - Gets the norm of the current function that was set 970 with SNESSSetFunction(). 971 972 Collective on SNES 973 974 Input Parameter: 975 . snes - SNES context 976 977 Output Parameter: 978 . fnorm - 2-norm of function 979 980 Level: intermediate 981 982 .keywords: SNES, nonlinear, get, function, norm 983 984 .seealso: SNESGetFunction(), SNESGetIterationNumber(), SNESGetLinearSolveIterations() 985 @*/ 986 PetscErrorCode SNESGetFunctionNorm(SNES snes,PetscReal *fnorm) 987 { 988 PetscFunctionBegin; 989 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 990 PetscValidScalarPointer(fnorm,2); 991 *fnorm = snes->norm; 992 PetscFunctionReturn(0); 993 } 994 995 996 #undef __FUNCT__ 997 #define __FUNCT__ "SNESSetFunctionNorm" 998 /*@ 999 SNESSetFunctionNorm - Sets the 2-norm of the current function computed using VecNorm(). 1000 1001 Collective on SNES 1002 1003 Input Parameter: 1004 . snes - SNES context 1005 . fnorm - 2-norm of function 1006 1007 Level: developer 1008 1009 .keywords: SNES, nonlinear, set, function, norm 1010 1011 .seealso: SNESSetFunction(), SNESSetIterationNumber(), VecNorm(). 1012 @*/ 1013 PetscErrorCode SNESSetFunctionNorm(SNES snes,PetscReal fnorm) 1014 { 1015 1016 PetscErrorCode ierr; 1017 1018 PetscFunctionBegin; 1019 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1020 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 1021 snes->norm = fnorm; 1022 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 1023 PetscFunctionReturn(0); 1024 } 1025 1026 #undef __FUNCT__ 1027 #define __FUNCT__ "SNESGetNonlinearStepFailures" 1028 /*@ 1029 SNESGetNonlinearStepFailures - Gets the number of unsuccessful steps 1030 attempted by the nonlinear solver. 1031 1032 Not Collective 1033 1034 Input Parameter: 1035 . snes - SNES context 1036 1037 Output Parameter: 1038 . nfails - number of unsuccessful steps attempted 1039 1040 Notes: 1041 This counter is reset to zero for each successive call to SNESSolve(). 1042 1043 Level: intermediate 1044 1045 .keywords: SNES, nonlinear, get, number, unsuccessful, steps 1046 1047 .seealso: SNESGetMaxLinearSolveFailures(), SNESGetLinearSolveIterations(), SNESSetMaxLinearSolveFailures(), SNESGetLinearSolveFailures(), 1048 SNESSetMaxNonlinearStepFailures(), SNESGetMaxNonlinearStepFailures() 1049 @*/ 1050 PetscErrorCode SNESGetNonlinearStepFailures(SNES snes,PetscInt* nfails) 1051 { 1052 PetscFunctionBegin; 1053 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1054 PetscValidIntPointer(nfails,2); 1055 *nfails = snes->numFailures; 1056 PetscFunctionReturn(0); 1057 } 1058 1059 #undef __FUNCT__ 1060 #define __FUNCT__ "SNESSetMaxNonlinearStepFailures" 1061 /*@ 1062 SNESSetMaxNonlinearStepFailures - Sets the maximum number of unsuccessful steps 1063 attempted by the nonlinear solver before it gives up. 1064 1065 Not Collective 1066 1067 Input Parameters: 1068 + snes - SNES context 1069 - maxFails - maximum of unsuccessful steps 1070 1071 Level: intermediate 1072 1073 .keywords: SNES, nonlinear, set, maximum, unsuccessful, steps 1074 1075 .seealso: SNESGetMaxLinearSolveFailures(), SNESGetLinearSolveIterations(), SNESSetMaxLinearSolveFailures(), SNESGetLinearSolveFailures(), 1076 SNESGetMaxNonlinearStepFailures(), SNESGetNonlinearStepFailures() 1077 @*/ 1078 PetscErrorCode SNESSetMaxNonlinearStepFailures(SNES snes, PetscInt maxFails) 1079 { 1080 PetscFunctionBegin; 1081 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1082 snes->maxFailures = maxFails; 1083 PetscFunctionReturn(0); 1084 } 1085 1086 #undef __FUNCT__ 1087 #define __FUNCT__ "SNESGetMaxNonlinearStepFailures" 1088 /*@ 1089 SNESGetMaxNonlinearStepFailures - Gets the maximum number of unsuccessful steps 1090 attempted by the nonlinear solver before it gives up. 1091 1092 Not Collective 1093 1094 Input Parameter: 1095 . snes - SNES context 1096 1097 Output Parameter: 1098 . maxFails - maximum of unsuccessful steps 1099 1100 Level: intermediate 1101 1102 .keywords: SNES, nonlinear, get, maximum, unsuccessful, steps 1103 1104 .seealso: SNESGetMaxLinearSolveFailures(), SNESGetLinearSolveIterations(), SNESSetMaxLinearSolveFailures(), SNESGetLinearSolveFailures(), 1105 SNESSetMaxNonlinearStepFailures(), SNESGetNonlinearStepFailures() 1106 1107 @*/ 1108 PetscErrorCode SNESGetMaxNonlinearStepFailures(SNES snes, PetscInt *maxFails) 1109 { 1110 PetscFunctionBegin; 1111 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1112 PetscValidIntPointer(maxFails,2); 1113 *maxFails = snes->maxFailures; 1114 PetscFunctionReturn(0); 1115 } 1116 1117 #undef __FUNCT__ 1118 #define __FUNCT__ "SNESGetNumberFunctionEvals" 1119 /*@ 1120 SNESGetNumberFunctionEvals - Gets the number of user provided function evaluations 1121 done by SNES. 1122 1123 Not Collective 1124 1125 Input Parameter: 1126 . snes - SNES context 1127 1128 Output Parameter: 1129 . nfuncs - number of evaluations 1130 1131 Level: intermediate 1132 1133 .keywords: SNES, nonlinear, get, maximum, unsuccessful, steps 1134 1135 .seealso: SNESGetMaxLinearSolveFailures(), SNESGetLinearSolveIterations(), SNESSetMaxLinearSolveFailures(), SNESGetLinearSolveFailures() 1136 @*/ 1137 PetscErrorCode SNESGetNumberFunctionEvals(SNES snes, PetscInt *nfuncs) 1138 { 1139 PetscFunctionBegin; 1140 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1141 PetscValidIntPointer(nfuncs,2); 1142 *nfuncs = snes->nfuncs; 1143 PetscFunctionReturn(0); 1144 } 1145 1146 #undef __FUNCT__ 1147 #define __FUNCT__ "SNESGetLinearSolveFailures" 1148 /*@ 1149 SNESGetLinearSolveFailures - Gets the number of failed (non-converged) 1150 linear solvers. 1151 1152 Not Collective 1153 1154 Input Parameter: 1155 . snes - SNES context 1156 1157 Output Parameter: 1158 . nfails - number of failed solves 1159 1160 Notes: 1161 This counter is reset to zero for each successive call to SNESSolve(). 1162 1163 Level: intermediate 1164 1165 .keywords: SNES, nonlinear, get, number, unsuccessful, steps 1166 1167 .seealso: SNESGetMaxLinearSolveFailures(), SNESGetLinearSolveIterations(), SNESSetMaxLinearSolveFailures() 1168 @*/ 1169 PetscErrorCode SNESGetLinearSolveFailures(SNES snes,PetscInt* nfails) 1170 { 1171 PetscFunctionBegin; 1172 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1173 PetscValidIntPointer(nfails,2); 1174 *nfails = snes->numLinearSolveFailures; 1175 PetscFunctionReturn(0); 1176 } 1177 1178 #undef __FUNCT__ 1179 #define __FUNCT__ "SNESSetMaxLinearSolveFailures" 1180 /*@ 1181 SNESSetMaxLinearSolveFailures - the number of failed linear solve attempts 1182 allowed before SNES returns with a diverged reason of SNES_DIVERGED_LINEAR_SOLVE 1183 1184 Logically Collective on SNES 1185 1186 Input Parameters: 1187 + snes - SNES context 1188 - maxFails - maximum allowed linear solve failures 1189 1190 Level: intermediate 1191 1192 Notes: By default this is 0; that is SNES returns on the first failed linear solve 1193 1194 .keywords: SNES, nonlinear, set, maximum, unsuccessful, steps 1195 1196 .seealso: SNESGetLinearSolveFailures(), SNESGetMaxLinearSolveFailures(), SNESGetLinearSolveIterations() 1197 @*/ 1198 PetscErrorCode SNESSetMaxLinearSolveFailures(SNES snes, PetscInt maxFails) 1199 { 1200 PetscFunctionBegin; 1201 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1202 PetscValidLogicalCollectiveInt(snes,maxFails,2); 1203 snes->maxLinearSolveFailures = maxFails; 1204 PetscFunctionReturn(0); 1205 } 1206 1207 #undef __FUNCT__ 1208 #define __FUNCT__ "SNESGetMaxLinearSolveFailures" 1209 /*@ 1210 SNESGetMaxLinearSolveFailures - gets the maximum number of linear solve failures that 1211 are allowed before SNES terminates 1212 1213 Not Collective 1214 1215 Input Parameter: 1216 . snes - SNES context 1217 1218 Output Parameter: 1219 . maxFails - maximum of unsuccessful solves allowed 1220 1221 Level: intermediate 1222 1223 Notes: By default this is 1; that is SNES returns on the first failed linear solve 1224 1225 .keywords: SNES, nonlinear, get, maximum, unsuccessful, steps 1226 1227 .seealso: SNESGetLinearSolveFailures(), SNESGetLinearSolveIterations(), SNESSetMaxLinearSolveFailures(), 1228 @*/ 1229 PetscErrorCode SNESGetMaxLinearSolveFailures(SNES snes, PetscInt *maxFails) 1230 { 1231 PetscFunctionBegin; 1232 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1233 PetscValidIntPointer(maxFails,2); 1234 *maxFails = snes->maxLinearSolveFailures; 1235 PetscFunctionReturn(0); 1236 } 1237 1238 #undef __FUNCT__ 1239 #define __FUNCT__ "SNESGetLinearSolveIterations" 1240 /*@ 1241 SNESGetLinearSolveIterations - Gets the total number of linear iterations 1242 used by the nonlinear solver. 1243 1244 Not Collective 1245 1246 Input Parameter: 1247 . snes - SNES context 1248 1249 Output Parameter: 1250 . lits - number of linear iterations 1251 1252 Notes: 1253 This counter is reset to zero for each successive call to SNESSolve(). 1254 1255 Level: intermediate 1256 1257 .keywords: SNES, nonlinear, get, number, linear, iterations 1258 1259 .seealso: SNESGetIterationNumber(), SNESGetFunctionNorm(), SNESGetLinearSolveFailures(), SNESGetMaxLinearSolveFailures() 1260 @*/ 1261 PetscErrorCode SNESGetLinearSolveIterations(SNES snes,PetscInt* lits) 1262 { 1263 PetscFunctionBegin; 1264 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1265 PetscValidIntPointer(lits,2); 1266 *lits = snes->linear_its; 1267 PetscFunctionReturn(0); 1268 } 1269 1270 #undef __FUNCT__ 1271 #define __FUNCT__ "SNESGetKSP" 1272 /*@ 1273 SNESGetKSP - Returns the KSP context for a SNES solver. 1274 1275 Not Collective, but if SNES object is parallel, then KSP object is parallel 1276 1277 Input Parameter: 1278 . snes - the SNES context 1279 1280 Output Parameter: 1281 . ksp - the KSP context 1282 1283 Notes: 1284 The user can then directly manipulate the KSP context to set various 1285 options, etc. Likewise, the user can then extract and manipulate the 1286 PC contexts as well. 1287 1288 Level: beginner 1289 1290 .keywords: SNES, nonlinear, get, KSP, context 1291 1292 .seealso: KSPGetPC(), SNESCreate(), KSPCreate(), SNESSetKSP() 1293 @*/ 1294 PetscErrorCode SNESGetKSP(SNES snes,KSP *ksp) 1295 { 1296 PetscErrorCode ierr; 1297 1298 PetscFunctionBegin; 1299 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1300 PetscValidPointer(ksp,2); 1301 1302 if (!snes->ksp) { 1303 ierr = KSPCreate(((PetscObject)snes)->comm,&snes->ksp);CHKERRQ(ierr); 1304 ierr = PetscObjectIncrementTabLevel((PetscObject)snes->ksp,(PetscObject)snes,1);CHKERRQ(ierr); 1305 ierr = PetscLogObjectParent(snes,snes->ksp);CHKERRQ(ierr); 1306 } 1307 *ksp = snes->ksp; 1308 PetscFunctionReturn(0); 1309 } 1310 1311 #undef __FUNCT__ 1312 #define __FUNCT__ "SNESSetKSP" 1313 /*@ 1314 SNESSetKSP - Sets a KSP context for the SNES object to use 1315 1316 Not Collective, but the SNES and KSP objects must live on the same MPI_Comm 1317 1318 Input Parameters: 1319 + snes - the SNES context 1320 - ksp - the KSP context 1321 1322 Notes: 1323 The SNES object already has its KSP object, you can obtain with SNESGetKSP() 1324 so this routine is rarely needed. 1325 1326 The KSP object that is already in the SNES object has its reference count 1327 decreased by one. 1328 1329 Level: developer 1330 1331 .keywords: SNES, nonlinear, get, KSP, context 1332 1333 .seealso: KSPGetPC(), SNESCreate(), KSPCreate(), SNESSetKSP() 1334 @*/ 1335 PetscErrorCode SNESSetKSP(SNES snes,KSP ksp) 1336 { 1337 PetscErrorCode ierr; 1338 1339 PetscFunctionBegin; 1340 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1341 PetscValidHeaderSpecific(ksp,KSP_CLASSID,2); 1342 PetscCheckSameComm(snes,1,ksp,2); 1343 ierr = PetscObjectReference((PetscObject)ksp);CHKERRQ(ierr); 1344 if (snes->ksp) {ierr = PetscObjectDereference((PetscObject)snes->ksp);CHKERRQ(ierr);} 1345 snes->ksp = ksp; 1346 PetscFunctionReturn(0); 1347 } 1348 1349 #if 0 1350 #undef __FUNCT__ 1351 #define __FUNCT__ "SNESPublish_Petsc" 1352 static PetscErrorCode SNESPublish_Petsc(PetscObject obj) 1353 { 1354 PetscFunctionBegin; 1355 PetscFunctionReturn(0); 1356 } 1357 #endif 1358 1359 /* -----------------------------------------------------------*/ 1360 #undef __FUNCT__ 1361 #define __FUNCT__ "SNESCreate" 1362 /*@ 1363 SNESCreate - Creates a nonlinear solver context. 1364 1365 Collective on MPI_Comm 1366 1367 Input Parameters: 1368 . comm - MPI communicator 1369 1370 Output Parameter: 1371 . outsnes - the new SNES context 1372 1373 Options Database Keys: 1374 + -snes_mf - Activates default matrix-free Jacobian-vector products, 1375 and no preconditioning matrix 1376 . -snes_mf_operator - Activates default matrix-free Jacobian-vector 1377 products, and a user-provided preconditioning matrix 1378 as set by SNESSetJacobian() 1379 - -snes_fd - Uses (slow!) finite differences to compute Jacobian 1380 1381 Level: beginner 1382 1383 .keywords: SNES, nonlinear, create, context 1384 1385 .seealso: SNESSolve(), SNESDestroy(), SNES, SNESSetLagPreconditioner() 1386 1387 @*/ 1388 PetscErrorCode SNESCreate(MPI_Comm comm,SNES *outsnes) 1389 { 1390 PetscErrorCode ierr; 1391 SNES snes; 1392 SNESKSPEW *kctx; 1393 1394 PetscFunctionBegin; 1395 PetscValidPointer(outsnes,2); 1396 *outsnes = PETSC_NULL; 1397 #ifndef PETSC_USE_DYNAMIC_LIBRARIES 1398 ierr = SNESInitializePackage(PETSC_NULL);CHKERRQ(ierr); 1399 #endif 1400 1401 ierr = PetscHeaderCreate(snes,_p_SNES,struct _SNESOps,SNES_CLASSID,0,"SNES","Nonlinear solver","SNES",comm,SNESDestroy,SNESView);CHKERRQ(ierr); 1402 1403 snes->ops->converged = SNESDefaultConverged; 1404 snes->usesksp = PETSC_TRUE; 1405 snes->tolerancesset = PETSC_FALSE; 1406 snes->max_its = 50; 1407 snes->max_funcs = 10000; 1408 snes->norm = 0.0; 1409 snes->normtype = SNES_NORM_FUNCTION; 1410 snes->rtol = 1.e-8; 1411 snes->ttol = 0.0; 1412 snes->abstol = 1.e-50; 1413 snes->stol = 1.e-8; 1414 snes->deltatol = 1.e-12; 1415 snes->nfuncs = 0; 1416 snes->numFailures = 0; 1417 snes->maxFailures = 1; 1418 snes->linear_its = 0; 1419 snes->lagjacobian = 1; 1420 snes->lagpreconditioner = 1; 1421 snes->numbermonitors = 0; 1422 snes->data = 0; 1423 snes->setupcalled = PETSC_FALSE; 1424 snes->ksp_ewconv = PETSC_FALSE; 1425 snes->nwork = 0; 1426 snes->work = 0; 1427 snes->nvwork = 0; 1428 snes->vwork = 0; 1429 snes->conv_hist_len = 0; 1430 snes->conv_hist_max = 0; 1431 snes->conv_hist = PETSC_NULL; 1432 snes->conv_hist_its = PETSC_NULL; 1433 snes->conv_hist_reset = PETSC_TRUE; 1434 snes->vec_func_init_set = PETSC_FALSE; 1435 snes->norm_init = 0.; 1436 snes->norm_init_set = PETSC_FALSE; 1437 snes->reason = SNES_CONVERGED_ITERATING; 1438 snes->gssweeps = 1; 1439 1440 snes->pcside = PC_RIGHT; 1441 1442 snes->mf = PETSC_FALSE; 1443 snes->mf_operator = PETSC_FALSE; 1444 snes->mf_version = 1; 1445 1446 snes->numLinearSolveFailures = 0; 1447 snes->maxLinearSolveFailures = 1; 1448 1449 /* Create context to compute Eisenstat-Walker relative tolerance for KSP */ 1450 ierr = PetscNewLog(snes,SNESKSPEW,&kctx);CHKERRQ(ierr); 1451 snes->kspconvctx = (void*)kctx; 1452 kctx->version = 2; 1453 kctx->rtol_0 = .3; /* Eisenstat and Walker suggest rtol_0=.5, but 1454 this was too large for some test cases */ 1455 kctx->rtol_last = 0.0; 1456 kctx->rtol_max = .9; 1457 kctx->gamma = 1.0; 1458 kctx->alpha = .5*(1.0 + PetscSqrtReal(5.0)); 1459 kctx->alpha2 = kctx->alpha; 1460 kctx->threshold = .1; 1461 kctx->lresid_last = 0.0; 1462 kctx->norm_last = 0.0; 1463 1464 *outsnes = snes; 1465 PetscFunctionReturn(0); 1466 } 1467 1468 #undef __FUNCT__ 1469 #define __FUNCT__ "SNESSetFunction" 1470 /*@C 1471 SNESSetFunction - Sets the function evaluation routine and function 1472 vector for use by the SNES routines in solving systems of nonlinear 1473 equations. 1474 1475 Logically Collective on SNES 1476 1477 Input Parameters: 1478 + snes - the SNES context 1479 . r - vector to store function value 1480 . func - function evaluation routine 1481 - ctx - [optional] user-defined context for private data for the 1482 function evaluation routine (may be PETSC_NULL) 1483 1484 Calling sequence of func: 1485 $ func (SNES snes,Vec x,Vec f,void *ctx); 1486 1487 + snes - the SNES context 1488 . x - state at which to evaluate residual 1489 . f - vector to put residual 1490 - ctx - optional user-defined function context 1491 1492 Notes: 1493 The Newton-like methods typically solve linear systems of the form 1494 $ f'(x) x = -f(x), 1495 where f'(x) denotes the Jacobian matrix and f(x) is the function. 1496 1497 Level: beginner 1498 1499 .keywords: SNES, nonlinear, set, function 1500 1501 .seealso: SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESSetPicard() 1502 @*/ 1503 PetscErrorCode SNESSetFunction(SNES snes,Vec r,PetscErrorCode (*func)(SNES,Vec,Vec,void*),void *ctx) 1504 { 1505 PetscErrorCode ierr; 1506 DM dm; 1507 1508 PetscFunctionBegin; 1509 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1510 if (r) { 1511 PetscValidHeaderSpecific(r,VEC_CLASSID,2); 1512 PetscCheckSameComm(snes,1,r,2); 1513 ierr = PetscObjectReference((PetscObject)r);CHKERRQ(ierr); 1514 ierr = VecDestroy(&snes->vec_func);CHKERRQ(ierr); 1515 snes->vec_func = r; 1516 } 1517 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 1518 ierr = DMSNESSetFunction(dm,func,ctx);CHKERRQ(ierr); 1519 PetscFunctionReturn(0); 1520 } 1521 1522 1523 #undef __FUNCT__ 1524 #define __FUNCT__ "SNESSetInitialFunction" 1525 /*@C 1526 SNESSetInitialFunction - Sets the function vector to be used as the 1527 function norm at the initialization of the method. In some 1528 instances, the user has precomputed the function before calling 1529 SNESSolve. This function allows one to avoid a redundant call 1530 to SNESComputeFunction in that case. 1531 1532 Logically Collective on SNES 1533 1534 Input Parameters: 1535 + snes - the SNES context 1536 - f - vector to store function value 1537 1538 Notes: 1539 This should not be modified during the solution procedure. 1540 1541 This is used extensively in the SNESFAS hierarchy and in nonlinear preconditioning. 1542 1543 Level: developer 1544 1545 .keywords: SNES, nonlinear, set, function 1546 1547 .seealso: SNESSetFunction(), SNESComputeFunction(), SNESSetInitialFunctionNorm() 1548 @*/ 1549 PetscErrorCode SNESSetInitialFunction(SNES snes, Vec f) 1550 { 1551 PetscErrorCode ierr; 1552 Vec vec_func; 1553 1554 PetscFunctionBegin; 1555 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1556 PetscValidHeaderSpecific(f,VEC_CLASSID,2); 1557 PetscCheckSameComm(snes,1,f,2); 1558 ierr = SNESGetFunction(snes,&vec_func,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr); 1559 ierr = VecCopy(f, vec_func);CHKERRQ(ierr); 1560 snes->vec_func_init_set = PETSC_TRUE; 1561 PetscFunctionReturn(0); 1562 } 1563 1564 1565 #undef __FUNCT__ 1566 #define __FUNCT__ "SNESSetInitialFunctionNorm" 1567 /*@C 1568 SNESSetInitialFunctionNorm - Sets the function norm to be used as the function norm 1569 at the initialization of the method. In some instances, the user has precomputed 1570 the function and its norm before calling SNESSolve. This function allows one to 1571 avoid a redundant call to SNESComputeFunction() and VecNorm() in that case. 1572 1573 Logically Collective on SNES 1574 1575 Input Parameters: 1576 + snes - the SNES context 1577 - fnorm - the norm of F as set by SNESSetInitialFunction() 1578 1579 This is used extensively in the SNESFAS hierarchy and in nonlinear preconditioning. 1580 1581 Level: developer 1582 1583 .keywords: SNES, nonlinear, set, function, norm 1584 1585 .seealso: SNESSetFunction(), SNESComputeFunction(), SNESSetInitialFunction() 1586 @*/ 1587 PetscErrorCode SNESSetInitialFunctionNorm(SNES snes, PetscReal fnorm) 1588 { 1589 PetscFunctionBegin; 1590 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1591 snes->norm_init = fnorm; 1592 snes->norm_init_set = PETSC_TRUE; 1593 PetscFunctionReturn(0); 1594 } 1595 1596 #undef __FUNCT__ 1597 #define __FUNCT__ "SNESSetNormType" 1598 /*@ 1599 SNESSetNormType - Sets the SNESNormType used in covergence and monitoring 1600 of the SNES method. 1601 1602 Logically Collective on SNES 1603 1604 Input Parameters: 1605 + snes - the SNES context 1606 - normtype - the type of the norm used 1607 1608 Notes: 1609 Only certain SNES methods support certain SNESNormTypes. Most require evaluation 1610 of the nonlinear function and the taking of its norm at every iteration to 1611 even ensure convergence at all. However, methods such as custom Gauss-Seidel methods 1612 (SNESGS) and the like do not require the norm of the function to be computed, and therfore 1613 may either be monitored for convergence or not. As these are often used as nonlinear 1614 preconditioners, monitoring the norm of their error is not a useful enterprise within 1615 their solution. 1616 1617 Level: developer 1618 1619 .keywords: SNES, nonlinear, set, function, norm, type 1620 1621 .seealso: SNESGetNormType(), SNESComputeFunction(), VecNorm(), SNESSetFunction(), SNESSetInitialFunction(), SNESNormType 1622 @*/ 1623 PetscErrorCode SNESSetNormType(SNES snes, SNESNormType normtype) 1624 { 1625 PetscFunctionBegin; 1626 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1627 snes->normtype = normtype; 1628 PetscFunctionReturn(0); 1629 } 1630 1631 1632 #undef __FUNCT__ 1633 #define __FUNCT__ "SNESGetNormType" 1634 /*@ 1635 SNESGetNormType - Gets the SNESNormType used in covergence and monitoring 1636 of the SNES method. 1637 1638 Logically Collective on SNES 1639 1640 Input Parameters: 1641 + snes - the SNES context 1642 - normtype - the type of the norm used 1643 1644 Level: advanced 1645 1646 .keywords: SNES, nonlinear, set, function, norm, type 1647 1648 .seealso: SNESSetNormType(), SNESComputeFunction(), VecNorm(), SNESSetFunction(), SNESSetInitialFunction(), SNESNormType 1649 @*/ 1650 PetscErrorCode SNESGetNormType(SNES snes, SNESNormType *normtype) 1651 { 1652 PetscFunctionBegin; 1653 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1654 *normtype = snes->normtype; 1655 PetscFunctionReturn(0); 1656 } 1657 1658 #undef __FUNCT__ 1659 #define __FUNCT__ "SNESSetGS" 1660 /*@C 1661 SNESSetGS - Sets the user nonlinear Gauss-Seidel routine for 1662 use with composed nonlinear solvers. 1663 1664 Input Parameters: 1665 + snes - the SNES context 1666 . gsfunc - function evaluation routine 1667 - ctx - [optional] user-defined context for private data for the 1668 smoother evaluation routine (may be PETSC_NULL) 1669 1670 Calling sequence of func: 1671 $ func (SNES snes,Vec x,Vec b,void *ctx); 1672 1673 + X - solution vector 1674 . B - RHS vector 1675 - ctx - optional user-defined Gauss-Seidel context 1676 1677 Notes: 1678 The GS routines are used by the composed nonlinear solver to generate 1679 a problem appropriate update to the solution, particularly FAS. 1680 1681 Level: intermediate 1682 1683 .keywords: SNES, nonlinear, set, Gauss-Seidel 1684 1685 .seealso: SNESGetFunction(), SNESComputeGS() 1686 @*/ 1687 PetscErrorCode SNESSetGS(SNES snes,PetscErrorCode (*gsfunc)(SNES,Vec,Vec,void*),void *ctx) 1688 { 1689 PetscErrorCode ierr; 1690 DM dm; 1691 1692 PetscFunctionBegin; 1693 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1694 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 1695 ierr = DMSNESSetGS(dm,gsfunc,ctx);CHKERRQ(ierr); 1696 PetscFunctionReturn(0); 1697 } 1698 1699 #undef __FUNCT__ 1700 #define __FUNCT__ "SNESSetGSSweeps" 1701 /*@ 1702 SNESSetGSSweeps - Sets the number of sweeps of GS to use. 1703 1704 Input Parameters: 1705 + snes - the SNES context 1706 - sweeps - the number of sweeps of GS to perform. 1707 1708 Level: intermediate 1709 1710 .keywords: SNES, nonlinear, set, Gauss-Siedel 1711 1712 .seealso: SNESSetGS(), SNESGetGS(), SNESSetPC(), SNESGetGSSweeps() 1713 @*/ 1714 1715 PetscErrorCode SNESSetGSSweeps(SNES snes, PetscInt sweeps) 1716 { 1717 PetscFunctionBegin; 1718 snes->gssweeps = sweeps; 1719 PetscFunctionReturn(0); 1720 } 1721 1722 1723 #undef __FUNCT__ 1724 #define __FUNCT__ "SNESGetGSSweeps" 1725 /*@ 1726 SNESGetGSSweeps - Gets the number of sweeps GS will use. 1727 1728 Input Parameters: 1729 . snes - the SNES context 1730 1731 Output Parameters: 1732 . sweeps - the number of sweeps of GS to perform. 1733 1734 Level: intermediate 1735 1736 .keywords: SNES, nonlinear, set, Gauss-Siedel 1737 1738 .seealso: SNESSetGS(), SNESGetGS(), SNESSetPC(), SNESSetGSSweeps() 1739 @*/ 1740 PetscErrorCode SNESGetGSSweeps(SNES snes, PetscInt * sweeps) 1741 { 1742 PetscFunctionBegin; 1743 *sweeps = snes->gssweeps; 1744 PetscFunctionReturn(0); 1745 } 1746 1747 #undef __FUNCT__ 1748 #define __FUNCT__ "SNESPicardComputeFunction" 1749 PetscErrorCode SNESPicardComputeFunction(SNES snes,Vec x,Vec f,void *ctx) 1750 { 1751 PetscErrorCode ierr; 1752 DM dm; 1753 SNESDM sdm; 1754 1755 PetscFunctionBegin; 1756 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 1757 ierr = DMSNESGetContext(dm,&sdm);CHKERRQ(ierr); 1758 /* A(x)*x - b(x) */ 1759 if (sdm->computepfunction) { 1760 ierr = (*sdm->computepfunction)(snes,x,f,sdm->pctx);CHKERRQ(ierr); 1761 } else if (snes->dm) { 1762 ierr = DMComputeFunction(snes->dm,x,f);CHKERRQ(ierr); 1763 } else { 1764 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE, "Must call SNESSetPicard() or SNESSetDM() before SNESPicardComputeFunction to provide Picard function."); 1765 } 1766 1767 if (sdm->computepjacobian) { 1768 ierr = (*sdm->computepjacobian)(snes,x,&snes->jacobian,&snes->jacobian_pre,&snes->matstruct,sdm->pctx);CHKERRQ(ierr); 1769 } else if (snes->dm) { 1770 ierr = DMComputeJacobian(snes->dm,x,snes->jacobian,snes->jacobian_pre,&snes->matstruct);CHKERRQ(ierr); 1771 } else { 1772 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE, "Must call SNESSetPicard() or SNESSetDM() before SNESPicardComputeFunction to provide Picard matrix."); 1773 } 1774 1775 ierr = VecView(x,PETSC_VIEWER_BINARY_WORLD);CHKERRQ(ierr); 1776 ierr = VecView(f,PETSC_VIEWER_BINARY_WORLD);CHKERRQ(ierr); 1777 ierr = VecScale(f,-1.0);CHKERRQ(ierr); 1778 ierr = MatMultAdd(snes->jacobian_pre,x,f,f);CHKERRQ(ierr); 1779 PetscFunctionReturn(0); 1780 } 1781 1782 #undef __FUNCT__ 1783 #define __FUNCT__ "SNESPicardComputeJacobian" 1784 PetscErrorCode SNESPicardComputeJacobian(SNES snes,Vec x1,Mat *J,Mat *B,MatStructure *flag,void *ctx) 1785 { 1786 PetscFunctionBegin; 1787 /* the jacobian matrix should be pre-filled in SNESPicardComputeFunction */ 1788 *flag = snes->matstruct; 1789 PetscFunctionReturn(0); 1790 } 1791 1792 #undef __FUNCT__ 1793 #define __FUNCT__ "SNESSetPicard" 1794 /*@C 1795 SNESSetPicard - Use SNES to solve the semilinear-system A(x) x = b(x) via a Picard type iteration (Picard linearization) 1796 1797 Logically Collective on SNES 1798 1799 Input Parameters: 1800 + snes - the SNES context 1801 . r - vector to store function value 1802 . func - function evaluation routine 1803 . jmat - normally the same as mat but you can pass another matrix for which you compute the Jacobian of A(x) x - b(x) (see jmat below) 1804 . mat - matrix to store A 1805 . mfunc - function to compute matrix value 1806 - ctx - [optional] user-defined context for private data for the 1807 function evaluation routine (may be PETSC_NULL) 1808 1809 Calling sequence of func: 1810 $ func (SNES snes,Vec x,Vec f,void *ctx); 1811 1812 + f - function vector 1813 - ctx - optional user-defined function context 1814 1815 Calling sequence of mfunc: 1816 $ mfunc (SNES snes,Vec x,Mat *jmat,Mat *mat,int *flag,void *ctx); 1817 1818 + x - input vector 1819 . jmat - Form Jacobian matrix of A(x) x - b(x) if available, not there is really no reason to use it in this way since then you can just use SNESSetJacobian(), 1820 normally just pass mat in this location 1821 . mat - form A(x) matrix 1822 . flag - flag indicating information about the preconditioner matrix 1823 structure (same as flag in KSPSetOperators()), one of SAME_NONZERO_PATTERN,DIFFERENT_NONZERO_PATTERN,SAME_PRECONDITIONER 1824 - ctx - [optional] user-defined Jacobian context 1825 1826 Notes: 1827 One can call SNESSetPicard() or SNESSetFunction() (and possibly SNESSetJacobian()) but cannot call both 1828 1829 $ Solves the equation A(x) x = b(x) via the defect correction algorithm A(x^{n}) (x^{n+1} - x^{n}) = b(x^{n}) - A(x^{n})x^{n} 1830 $ Note that when an exact solver is used this corresponds to the "classic" Picard A(x^{n}) x^{n+1} = b(x^{n}) iteration. 1831 1832 Run with -snes_mf_operator to solve the system with Newton's method using A(x^{n}) to construct the preconditioner. 1833 1834 We implement the defect correction form of the Picard iteration because it converges much more generally when inexact linear solvers are used then 1835 the direct Picard iteration A(x^n) x^{n+1} = b(x^n) 1836 1837 There is some controversity over the definition of a Picard iteration for nonlinear systems but almost everyone agrees that it involves a linear solve and some 1838 believe it is the iteration A(x^{n}) x^{n+1} = b(x^{n}) hence we use the name Picard. If anyone has an authoritative reference that defines the Picard iteration 1839 different please contact us at petsc-dev@mcs.anl.gov and we'll have an entirely new argument :-). 1840 1841 Level: beginner 1842 1843 .keywords: SNES, nonlinear, set, function 1844 1845 .seealso: SNESGetFunction(), SNESSetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESGetPicard(), SNESLineSearchPreCheckPicard() 1846 @*/ 1847 PetscErrorCode SNESSetPicard(SNES snes,Vec r,PetscErrorCode (*func)(SNES,Vec,Vec,void*),Mat jmat, Mat mat, PetscErrorCode (*mfunc)(SNES,Vec,Mat*,Mat*,MatStructure*,void*),void *ctx) 1848 { 1849 PetscErrorCode ierr; 1850 DM dm; 1851 1852 PetscFunctionBegin; 1853 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1854 ierr = SNESGetDM(snes, &dm);CHKERRQ(ierr); 1855 ierr = DMSNESSetPicard(dm,func,mfunc,ctx);CHKERRQ(ierr); 1856 ierr = SNESSetFunction(snes,r,SNESPicardComputeFunction,ctx);CHKERRQ(ierr); 1857 ierr = SNESSetJacobian(snes,jmat,mat,SNESPicardComputeJacobian,ctx);CHKERRQ(ierr); 1858 PetscFunctionReturn(0); 1859 } 1860 1861 1862 #undef __FUNCT__ 1863 #define __FUNCT__ "SNESGetPicard" 1864 /*@C 1865 SNESGetPicard - Returns the context for the Picard iteration 1866 1867 Not Collective, but Vec is parallel if SNES is parallel. Collective if Vec is requested, but has not been created yet. 1868 1869 Input Parameter: 1870 . snes - the SNES context 1871 1872 Output Parameter: 1873 + r - the function (or PETSC_NULL) 1874 . func - the function (or PETSC_NULL) 1875 . jmat - the picard matrix (or PETSC_NULL) 1876 . mat - the picard preconditioner matrix (or PETSC_NULL) 1877 . mfunc - the function for matrix evaluation (or PETSC_NULL) 1878 - ctx - the function context (or PETSC_NULL) 1879 1880 Level: advanced 1881 1882 .keywords: SNES, nonlinear, get, function 1883 1884 .seealso: SNESSetPicard, SNESGetFunction, SNESGetJacobian, SNESGetDM 1885 @*/ 1886 PetscErrorCode SNESGetPicard(SNES snes,Vec *r,PetscErrorCode (**func)(SNES,Vec,Vec,void*),Mat *jmat, Mat *mat, PetscErrorCode (**mfunc)(SNES,Vec,Mat*,Mat*,MatStructure*,void*),void **ctx) 1887 { 1888 PetscErrorCode ierr; 1889 DM dm; 1890 1891 PetscFunctionBegin; 1892 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1893 ierr = SNESGetFunction(snes,r,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr); 1894 ierr = SNESGetJacobian(snes,jmat,mat,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr); 1895 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 1896 ierr = DMSNESGetPicard(dm,func,mfunc,ctx);CHKERRQ(ierr); 1897 PetscFunctionReturn(0); 1898 } 1899 1900 #undef __FUNCT__ 1901 #define __FUNCT__ "SNESSetComputeInitialGuess" 1902 /*@C 1903 SNESSetComputeInitialGuess - Sets a routine used to compute an initial guess for the problem 1904 1905 Logically Collective on SNES 1906 1907 Input Parameters: 1908 + snes - the SNES context 1909 . func - function evaluation routine 1910 - ctx - [optional] user-defined context for private data for the 1911 function evaluation routine (may be PETSC_NULL) 1912 1913 Calling sequence of func: 1914 $ func (SNES snes,Vec x,void *ctx); 1915 1916 . f - function vector 1917 - ctx - optional user-defined function context 1918 1919 Level: intermediate 1920 1921 .keywords: SNES, nonlinear, set, function 1922 1923 .seealso: SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian() 1924 @*/ 1925 PetscErrorCode SNESSetComputeInitialGuess(SNES snes,PetscErrorCode (*func)(SNES,Vec,void*),void *ctx) 1926 { 1927 PetscFunctionBegin; 1928 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1929 if (func) snes->ops->computeinitialguess = func; 1930 if (ctx) snes->initialguessP = ctx; 1931 PetscFunctionReturn(0); 1932 } 1933 1934 /* --------------------------------------------------------------- */ 1935 #undef __FUNCT__ 1936 #define __FUNCT__ "SNESGetRhs" 1937 /*@C 1938 SNESGetRhs - Gets the vector for solving F(x) = rhs. If rhs is not set 1939 it assumes a zero right hand side. 1940 1941 Logically Collective on SNES 1942 1943 Input Parameter: 1944 . snes - the SNES context 1945 1946 Output Parameter: 1947 . rhs - the right hand side vector or PETSC_NULL if the right hand side vector is null 1948 1949 Level: intermediate 1950 1951 .keywords: SNES, nonlinear, get, function, right hand side 1952 1953 .seealso: SNESGetSolution(), SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESSetFunction() 1954 @*/ 1955 PetscErrorCode SNESGetRhs(SNES snes,Vec *rhs) 1956 { 1957 PetscFunctionBegin; 1958 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1959 PetscValidPointer(rhs,2); 1960 *rhs = snes->vec_rhs; 1961 PetscFunctionReturn(0); 1962 } 1963 1964 #undef __FUNCT__ 1965 #define __FUNCT__ "SNESComputeFunction" 1966 /*@ 1967 SNESComputeFunction - Calls the function that has been set with 1968 SNESSetFunction(). 1969 1970 Collective on SNES 1971 1972 Input Parameters: 1973 + snes - the SNES context 1974 - x - input vector 1975 1976 Output Parameter: 1977 . y - function vector, as set by SNESSetFunction() 1978 1979 Notes: 1980 SNESComputeFunction() is typically used within nonlinear solvers 1981 implementations, so most users would not generally call this routine 1982 themselves. 1983 1984 Level: developer 1985 1986 .keywords: SNES, nonlinear, compute, function 1987 1988 .seealso: SNESSetFunction(), SNESGetFunction() 1989 @*/ 1990 PetscErrorCode SNESComputeFunction(SNES snes,Vec x,Vec y) 1991 { 1992 PetscErrorCode ierr; 1993 DM dm; 1994 SNESDM sdm; 1995 1996 PetscFunctionBegin; 1997 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1998 PetscValidHeaderSpecific(x,VEC_CLASSID,2); 1999 PetscValidHeaderSpecific(y,VEC_CLASSID,3); 2000 PetscCheckSameComm(snes,1,x,2); 2001 PetscCheckSameComm(snes,1,y,3); 2002 VecValidValues(x,2,PETSC_TRUE); 2003 2004 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 2005 ierr = DMSNESGetContext(dm,&sdm);CHKERRQ(ierr); 2006 ierr = PetscLogEventBegin(SNES_FunctionEval,snes,x,y,0);CHKERRQ(ierr); 2007 if (snes->pc && snes->pcside == PC_LEFT) { 2008 ierr = VecCopy(x,y);CHKERRQ(ierr); 2009 ierr = SNESSolve(snes->pc,snes->vec_rhs,y);CHKERRQ(ierr); 2010 ierr = VecAYPX(y,-1.0,x);CHKERRQ(ierr); 2011 } else if (sdm->computefunction) { 2012 PetscStackPush("SNES user function"); 2013 ierr = (*sdm->computefunction)(snes,x,y,sdm->functionctx);CHKERRQ(ierr); 2014 PetscStackPop; 2015 } else if (snes->dm) { 2016 ierr = DMComputeFunction(snes->dm,x,y);CHKERRQ(ierr); 2017 } else if (snes->vec_rhs) { 2018 ierr = MatMult(snes->jacobian, x, y);CHKERRQ(ierr); 2019 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE, "Must call SNESSetFunction() or SNESSetDM() before SNESComputeFunction(), likely called from SNESSolve()."); 2020 if (snes->vec_rhs) { 2021 ierr = VecAXPY(y,-1.0,snes->vec_rhs);CHKERRQ(ierr); 2022 } 2023 snes->nfuncs++; 2024 ierr = PetscLogEventEnd(SNES_FunctionEval,snes,x,y,0);CHKERRQ(ierr); 2025 VecValidValues(y,3,PETSC_FALSE); 2026 PetscFunctionReturn(0); 2027 } 2028 2029 #undef __FUNCT__ 2030 #define __FUNCT__ "SNESComputeGS" 2031 /*@ 2032 SNESComputeGS - Calls the Gauss-Seidel function that has been set with 2033 SNESSetGS(). 2034 2035 Collective on SNES 2036 2037 Input Parameters: 2038 + snes - the SNES context 2039 . x - input vector 2040 - b - rhs vector 2041 2042 Output Parameter: 2043 . x - new solution vector 2044 2045 Notes: 2046 SNESComputeGS() is typically used within composed nonlinear solver 2047 implementations, so most users would not generally call this routine 2048 themselves. 2049 2050 Level: developer 2051 2052 .keywords: SNES, nonlinear, compute, function 2053 2054 .seealso: SNESSetGS(), SNESComputeFunction() 2055 @*/ 2056 PetscErrorCode SNESComputeGS(SNES snes,Vec b,Vec x) 2057 { 2058 PetscErrorCode ierr; 2059 PetscInt i; 2060 DM dm; 2061 SNESDM sdm; 2062 2063 PetscFunctionBegin; 2064 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2065 PetscValidHeaderSpecific(x,VEC_CLASSID,2); 2066 if (b) PetscValidHeaderSpecific(b,VEC_CLASSID,3); 2067 PetscCheckSameComm(snes,1,x,2); 2068 if (b) PetscCheckSameComm(snes,1,b,3); 2069 if (b) VecValidValues(b,2,PETSC_TRUE); 2070 ierr = PetscLogEventBegin(SNES_GSEval,snes,x,b,0);CHKERRQ(ierr); 2071 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 2072 ierr = DMSNESGetContext(dm,&sdm);CHKERRQ(ierr); 2073 if (sdm->computegs) { 2074 for (i = 0; i < snes->gssweeps; i++) { 2075 PetscStackPush("SNES user GS"); 2076 ierr = (*sdm->computegs)(snes,x,b,sdm->gsctx);CHKERRQ(ierr); 2077 PetscStackPop; 2078 } 2079 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE, "Must call SNESSetGS() before SNESComputeGS(), likely called from SNESSolve()."); 2080 ierr = PetscLogEventEnd(SNES_GSEval,snes,x,b,0);CHKERRQ(ierr); 2081 VecValidValues(x,3,PETSC_FALSE); 2082 PetscFunctionReturn(0); 2083 } 2084 2085 2086 #undef __FUNCT__ 2087 #define __FUNCT__ "SNESComputeJacobian" 2088 /*@ 2089 SNESComputeJacobian - Computes the Jacobian matrix that has been 2090 set with SNESSetJacobian(). 2091 2092 Collective on SNES and Mat 2093 2094 Input Parameters: 2095 + snes - the SNES context 2096 - x - input vector 2097 2098 Output Parameters: 2099 + A - Jacobian matrix 2100 . B - optional preconditioning matrix 2101 - flag - flag indicating matrix structure (one of, SAME_NONZERO_PATTERN,DIFFERENT_NONZERO_PATTERN,SAME_PRECONDITIONER) 2102 2103 Options Database Keys: 2104 + -snes_lag_preconditioner <lag> 2105 . -snes_lag_jacobian <lag> 2106 . -snes_compare_explicit - Compare the computed Jacobian to the finite difference Jacobian and output the differences 2107 . -snes_compare_explicit_draw - Compare the computed Jacobian to the finite difference Jacobian and draw the result 2108 . -snes_compare_explicit_contour - Compare the computed Jacobian to the finite difference Jacobian and draw a contour plot with the result 2109 . -snes_compare_operator - Make the comparison options above use the operator instead of the preconditioning matrix 2110 . -snes_compare_coloring - Compute the finite differece Jacobian using coloring and display norms of difference 2111 . -snes_compare_coloring_display - Compute the finite differece Jacobian using coloring and display verbose differences 2112 . -snes_compare_coloring_threshold - Display only those matrix entries that differ by more than a given threshold 2113 . -snes_compare_coloring_threshold_atol - Absolute tolerance for difference in matrix entries to be displayed by -snes_compare_coloring_threshold 2114 . -snes_compare_coloring_threshold_rtol - Relative tolerance for difference in matrix entries to be displayed by -snes_compare_coloring_threshold 2115 . -snes_compare_coloring_draw - Compute the finite differece Jacobian using coloring and draw differences 2116 - -snes_compare_coloring_draw_contour - Compute the finite differece Jacobian using coloring and show contours of matrices and differences 2117 2118 2119 Notes: 2120 Most users should not need to explicitly call this routine, as it 2121 is used internally within the nonlinear solvers. 2122 2123 See KSPSetOperators() for important information about setting the 2124 flag parameter. 2125 2126 Level: developer 2127 2128 .keywords: SNES, compute, Jacobian, matrix 2129 2130 .seealso: SNESSetJacobian(), KSPSetOperators(), MatStructure, SNESSetLagPreconditioner(), SNESSetLagJacobian() 2131 @*/ 2132 PetscErrorCode SNESComputeJacobian(SNES snes,Vec X,Mat *A,Mat *B,MatStructure *flg) 2133 { 2134 PetscErrorCode ierr; 2135 PetscBool flag; 2136 DM dm; 2137 SNESDM sdm; 2138 2139 PetscFunctionBegin; 2140 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2141 PetscValidHeaderSpecific(X,VEC_CLASSID,2); 2142 PetscValidPointer(flg,5); 2143 PetscCheckSameComm(snes,1,X,2); 2144 VecValidValues(X,2,PETSC_TRUE); 2145 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 2146 ierr = DMSNESGetContext(dm,&sdm);CHKERRQ(ierr); 2147 2148 if (!sdm->computejacobian) SETERRQ(((PetscObject)snes)->comm,PETSC_ERR_USER,"Must call SNESSetJacobian(), DMSNESSetJacobian(), DMDASNESSetJacobianLocal(), etc"); 2149 2150 /* make sure that MatAssemblyBegin/End() is called on A matrix if it is matrix free */ 2151 2152 if (snes->lagjacobian == -2) { 2153 snes->lagjacobian = -1; 2154 ierr = PetscInfo(snes,"Recomputing Jacobian/preconditioner because lag is -2 (means compute Jacobian, but then never again) \n");CHKERRQ(ierr); 2155 } else if (snes->lagjacobian == -1) { 2156 *flg = SAME_PRECONDITIONER; 2157 ierr = PetscInfo(snes,"Reusing Jacobian/preconditioner because lag is -1\n");CHKERRQ(ierr); 2158 ierr = PetscObjectTypeCompare((PetscObject)*A,MATMFFD,&flag);CHKERRQ(ierr); 2159 if (flag) { 2160 ierr = MatAssemblyBegin(*A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2161 ierr = MatAssemblyEnd(*A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2162 } 2163 PetscFunctionReturn(0); 2164 } else if (snes->lagjacobian > 1 && snes->iter % snes->lagjacobian) { 2165 *flg = SAME_PRECONDITIONER; 2166 ierr = PetscInfo2(snes,"Reusing Jacobian/preconditioner because lag is %D and SNES iteration is %D\n",snes->lagjacobian,snes->iter);CHKERRQ(ierr); 2167 ierr = PetscObjectTypeCompare((PetscObject)*A,MATMFFD,&flag);CHKERRQ(ierr); 2168 if (flag) { 2169 ierr = MatAssemblyBegin(*A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2170 ierr = MatAssemblyEnd(*A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2171 } 2172 PetscFunctionReturn(0); 2173 } 2174 2175 *flg = DIFFERENT_NONZERO_PATTERN; 2176 ierr = PetscLogEventBegin(SNES_JacobianEval,snes,X,*A,*B);CHKERRQ(ierr); 2177 2178 PetscStackPush("SNES user Jacobian function"); 2179 if (snes->mf && !snes->mf_operator) { 2180 ierr = MatMFFDComputeJacobian(snes,X,A,B,flg,sdm->jacobianctx);CHKERRQ(ierr); 2181 } else { 2182 ierr = (*sdm->computejacobian)(snes,X,A,B,flg,sdm->jacobianctx);CHKERRQ(ierr); 2183 } 2184 PetscStackPop; 2185 2186 ierr = PetscLogEventEnd(SNES_JacobianEval,snes,X,*A,*B);CHKERRQ(ierr); 2187 2188 if (snes->lagpreconditioner == -2) { 2189 ierr = PetscInfo(snes,"Rebuilding preconditioner exactly once since lag is -2\n");CHKERRQ(ierr); 2190 snes->lagpreconditioner = -1; 2191 } else if (snes->lagpreconditioner == -1) { 2192 *flg = SAME_PRECONDITIONER; 2193 ierr = PetscInfo(snes,"Reusing preconditioner because lag is -1\n");CHKERRQ(ierr); 2194 } else if (snes->lagpreconditioner > 1 && snes->iter % snes->lagpreconditioner) { 2195 *flg = SAME_PRECONDITIONER; 2196 ierr = PetscInfo2(snes,"Reusing preconditioner because lag is %D and SNES iteration is %D\n",snes->lagpreconditioner,snes->iter);CHKERRQ(ierr); 2197 } 2198 2199 /* make sure user returned a correct Jacobian and preconditioner */ 2200 /* PetscValidHeaderSpecific(*A,MAT_CLASSID,3); 2201 PetscValidHeaderSpecific(*B,MAT_CLASSID,4); */ 2202 { 2203 PetscBool flag = PETSC_FALSE,flag_draw = PETSC_FALSE,flag_contour = PETSC_FALSE,flag_operator = PETSC_FALSE; 2204 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_explicit",&flag,PETSC_NULL);CHKERRQ(ierr); 2205 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_explicit_draw",&flag_draw,PETSC_NULL);CHKERRQ(ierr); 2206 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_explicit_draw_contour",&flag_contour,PETSC_NULL);CHKERRQ(ierr); 2207 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_operator",&flag_operator,PETSC_NULL);CHKERRQ(ierr); 2208 if (flag || flag_draw || flag_contour) { 2209 Mat Bexp_mine = PETSC_NULL,Bexp,FDexp; 2210 MatStructure mstruct; 2211 PetscViewer vdraw,vstdout; 2212 PetscBool flg; 2213 if (flag_operator) { 2214 ierr = MatComputeExplicitOperator(*A,&Bexp_mine);CHKERRQ(ierr); 2215 Bexp = Bexp_mine; 2216 } else { 2217 /* See if the preconditioning matrix can be viewed and added directly */ 2218 ierr = PetscObjectTypeCompareAny((PetscObject)*B,&flg,MATSEQAIJ,MATMPIAIJ,MATSEQDENSE,MATMPIDENSE,MATSEQBAIJ,MATMPIBAIJ,MATSEQSBAIJ,MATMPIBAIJ,"");CHKERRQ(ierr); 2219 if (flg) Bexp = *B; 2220 else { 2221 /* If the "preconditioning" matrix is itself MATSHELL or some other type without direct support */ 2222 ierr = MatComputeExplicitOperator(*B,&Bexp_mine);CHKERRQ(ierr); 2223 Bexp = Bexp_mine; 2224 } 2225 } 2226 ierr = MatConvert(Bexp,MATSAME,MAT_INITIAL_MATRIX,&FDexp);CHKERRQ(ierr); 2227 ierr = SNESDefaultComputeJacobian(snes,X,&FDexp,&FDexp,&mstruct,NULL);CHKERRQ(ierr); 2228 ierr = PetscViewerASCIIGetStdout(((PetscObject)snes)->comm,&vstdout);CHKERRQ(ierr); 2229 if (flag_draw || flag_contour) { 2230 ierr = PetscViewerDrawOpen(((PetscObject)snes)->comm,0,"Explicit Jacobians",PETSC_DECIDE,PETSC_DECIDE,300,300,&vdraw);CHKERRQ(ierr); 2231 if (flag_contour) {ierr = PetscViewerPushFormat(vdraw,PETSC_VIEWER_DRAW_CONTOUR);CHKERRQ(ierr);} 2232 } else vdraw = PETSC_NULL; 2233 ierr = PetscViewerASCIIPrintf(vstdout,"Explicit %s\n",flag_operator?"Jacobian":"preconditioning Jacobian");CHKERRQ(ierr); 2234 if (flag) {ierr = MatView(Bexp,vstdout);CHKERRQ(ierr);} 2235 if (vdraw) {ierr = MatView(Bexp,vdraw);CHKERRQ(ierr);} 2236 ierr = PetscViewerASCIIPrintf(vstdout,"Finite difference Jacobian\n");CHKERRQ(ierr); 2237 if (flag) {ierr = MatView(FDexp,vstdout);CHKERRQ(ierr);} 2238 if (vdraw) {ierr = MatView(FDexp,vdraw);CHKERRQ(ierr);} 2239 ierr = MatAYPX(FDexp,-1.0,Bexp,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 2240 ierr = PetscViewerASCIIPrintf(vstdout,"User-provided matrix minus finite difference Jacobian\n");CHKERRQ(ierr); 2241 if (flag) {ierr = MatView(FDexp,vstdout);CHKERRQ(ierr);} 2242 if (vdraw) { /* Always use contour for the difference */ 2243 ierr = PetscViewerPushFormat(vdraw,PETSC_VIEWER_DRAW_CONTOUR);CHKERRQ(ierr); 2244 ierr = MatView(FDexp,vdraw);CHKERRQ(ierr); 2245 ierr = PetscViewerPopFormat(vdraw);CHKERRQ(ierr); 2246 } 2247 if (flag_contour) {ierr = PetscViewerPopFormat(vdraw);CHKERRQ(ierr);} 2248 ierr = PetscViewerDestroy(&vdraw);CHKERRQ(ierr); 2249 ierr = MatDestroy(&Bexp_mine);CHKERRQ(ierr); 2250 ierr = MatDestroy(&FDexp);CHKERRQ(ierr); 2251 } 2252 } 2253 { 2254 PetscBool flag = PETSC_FALSE,flag_display = PETSC_FALSE,flag_draw = PETSC_FALSE,flag_contour = PETSC_FALSE,flag_threshold = PETSC_FALSE; 2255 PetscReal threshold_atol = PETSC_SQRT_MACHINE_EPSILON,threshold_rtol = 10*PETSC_SQRT_MACHINE_EPSILON; 2256 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_coloring",&flag,PETSC_NULL);CHKERRQ(ierr); 2257 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_coloring_display",&flag_display,PETSC_NULL);CHKERRQ(ierr); 2258 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_coloring_draw",&flag_draw,PETSC_NULL);CHKERRQ(ierr); 2259 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_coloring_draw_contour",&flag_contour,PETSC_NULL);CHKERRQ(ierr); 2260 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_coloring_threshold",&flag_threshold,PETSC_NULL);CHKERRQ(ierr); 2261 ierr = PetscOptionsGetReal(((PetscObject)snes)->prefix,"-snes_compare_coloring_threshold_rtol",&threshold_rtol,PETSC_NULL);CHKERRQ(ierr); 2262 ierr = PetscOptionsGetReal(((PetscObject)snes)->prefix,"-snes_compare_coloring_threshold_atol",&threshold_atol,PETSC_NULL);CHKERRQ(ierr); 2263 if (flag || flag_display || flag_draw || flag_contour || flag_threshold) { 2264 Mat Bfd; 2265 MatStructure mstruct; 2266 PetscViewer vdraw,vstdout; 2267 ISColoring iscoloring; 2268 MatFDColoring matfdcoloring; 2269 PetscErrorCode (*func)(SNES,Vec,Vec,void*); 2270 void *funcctx; 2271 PetscReal norm1,norm2,normmax; 2272 2273 ierr = MatDuplicate(*B,MAT_DO_NOT_COPY_VALUES,&Bfd);CHKERRQ(ierr); 2274 ierr = MatGetColoring(Bfd,MATCOLORINGSL,&iscoloring);CHKERRQ(ierr); 2275 ierr = MatFDColoringCreate(Bfd,iscoloring,&matfdcoloring);CHKERRQ(ierr); 2276 ierr = ISColoringDestroy(&iscoloring);CHKERRQ(ierr); 2277 2278 /* This method of getting the function is currently unreliable since it doesn't work for DM local functions. */ 2279 ierr = SNESGetFunction(snes,PETSC_NULL,&func,&funcctx);CHKERRQ(ierr); 2280 ierr = MatFDColoringSetFunction(matfdcoloring,(PetscErrorCode(*)(void))func,funcctx);CHKERRQ(ierr); 2281 ierr = PetscObjectSetOptionsPrefix((PetscObject)matfdcoloring,((PetscObject)snes)->prefix);CHKERRQ(ierr); 2282 ierr = PetscObjectAppendOptionsPrefix((PetscObject)matfdcoloring,"coloring_");CHKERRQ(ierr); 2283 ierr = MatFDColoringSetFromOptions(matfdcoloring);CHKERRQ(ierr); 2284 ierr = MatFDColoringApply(Bfd,matfdcoloring,X,&mstruct,snes);CHKERRQ(ierr); 2285 ierr = MatFDColoringDestroy(&matfdcoloring);CHKERRQ(ierr); 2286 2287 ierr = PetscViewerASCIIGetStdout(((PetscObject)snes)->comm,&vstdout);CHKERRQ(ierr); 2288 if (flag_draw || flag_contour) { 2289 ierr = PetscViewerDrawOpen(((PetscObject)snes)->comm,0,"Colored Jacobians",PETSC_DECIDE,PETSC_DECIDE,300,300,&vdraw);CHKERRQ(ierr); 2290 if (flag_contour) {ierr = PetscViewerPushFormat(vdraw,PETSC_VIEWER_DRAW_CONTOUR);CHKERRQ(ierr);} 2291 } else vdraw = PETSC_NULL; 2292 ierr = PetscViewerASCIIPrintf(vstdout,"Explicit preconditioning Jacobian\n");CHKERRQ(ierr); 2293 if (flag_display) {ierr = MatView(*B,vstdout);CHKERRQ(ierr);} 2294 if (vdraw) {ierr = MatView(*B,vdraw);CHKERRQ(ierr);} 2295 ierr = PetscViewerASCIIPrintf(vstdout,"Colored Finite difference Jacobian\n");CHKERRQ(ierr); 2296 if (flag_display) {ierr = MatView(Bfd,vstdout);CHKERRQ(ierr);} 2297 if (vdraw) {ierr = MatView(Bfd,vdraw);CHKERRQ(ierr);} 2298 ierr = MatAYPX(Bfd,-1.0,*B,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 2299 ierr = MatNorm(Bfd,NORM_1,&norm1);CHKERRQ(ierr); 2300 ierr = MatNorm(Bfd,NORM_FROBENIUS,&norm2);CHKERRQ(ierr); 2301 ierr = MatNorm(Bfd,NORM_MAX,&normmax);CHKERRQ(ierr); 2302 ierr = PetscViewerASCIIPrintf(vstdout,"User-provided matrix minus finite difference Jacobian, norm1=%G normFrob=%G normmax=%G\n",norm1,norm2,normmax);CHKERRQ(ierr); 2303 if (flag_display) {ierr = MatView(Bfd,vstdout);CHKERRQ(ierr);} 2304 if (vdraw) { /* Always use contour for the difference */ 2305 ierr = PetscViewerPushFormat(vdraw,PETSC_VIEWER_DRAW_CONTOUR);CHKERRQ(ierr); 2306 ierr = MatView(Bfd,vdraw);CHKERRQ(ierr); 2307 ierr = PetscViewerPopFormat(vdraw);CHKERRQ(ierr); 2308 } 2309 if (flag_contour) {ierr = PetscViewerPopFormat(vdraw);CHKERRQ(ierr);} 2310 2311 if (flag_threshold) { 2312 PetscInt bs,rstart,rend,i; 2313 ierr = MatGetBlockSize(*B,&bs);CHKERRQ(ierr); 2314 ierr = MatGetOwnershipRange(*B,&rstart,&rend);CHKERRQ(ierr); 2315 for (i=rstart; i<rend; i++) { 2316 const PetscScalar *ba,*ca; 2317 const PetscInt *bj,*cj; 2318 PetscInt bn,cn,j,maxentrycol = -1,maxdiffcol = -1,maxrdiffcol = -1; 2319 PetscReal maxentry = 0,maxdiff = 0,maxrdiff = 0; 2320 ierr = MatGetRow(*B,i,&bn,&bj,&ba);CHKERRQ(ierr); 2321 ierr = MatGetRow(Bfd,i,&cn,&cj,&ca);CHKERRQ(ierr); 2322 if (bn != cn) SETERRQ(((PetscObject)*A)->comm,PETSC_ERR_PLIB,"Unexpected different nonzero pattern in -snes_compare_coloring_threshold"); 2323 for (j=0; j<bn; j++) { 2324 PetscReal rdiff = PetscAbsScalar(ca[j]) / (threshold_atol + threshold_rtol*PetscAbsScalar(ba[j])); 2325 if (PetscAbsScalar(ba[j]) > PetscAbs(maxentry)) { 2326 maxentrycol = bj[j]; 2327 maxentry = PetscRealPart(ba[j]); 2328 } 2329 if (PetscAbsScalar(ca[j]) > PetscAbs(maxdiff)) { 2330 maxdiffcol = bj[j]; 2331 maxdiff = PetscRealPart(ca[j]); 2332 } 2333 if (rdiff > maxrdiff) { 2334 maxrdiffcol = bj[j]; 2335 maxrdiff = rdiff; 2336 } 2337 } 2338 if (maxrdiff > 1) { 2339 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); 2340 for (j=0; j<bn; j++) { 2341 PetscReal rdiff; 2342 rdiff = PetscAbsScalar(ca[j]) / (threshold_atol + threshold_rtol*PetscAbsScalar(ba[j])); 2343 if (rdiff > 1) { 2344 ierr = PetscViewerASCIIPrintf(vstdout," (%D,%G:%G)",bj[j],PetscRealPart(ba[j]),PetscRealPart(ca[j]));CHKERRQ(ierr); 2345 } 2346 } 2347 ierr = PetscViewerASCIIPrintf(vstdout,"\n",i,maxentry,maxdiff,maxrdiff);CHKERRQ(ierr); 2348 } 2349 ierr = MatRestoreRow(*B,i,&bn,&bj,&ba);CHKERRQ(ierr); 2350 ierr = MatRestoreRow(Bfd,i,&cn,&cj,&ca);CHKERRQ(ierr); 2351 } 2352 } 2353 ierr = PetscViewerDestroy(&vdraw);CHKERRQ(ierr); 2354 ierr = MatDestroy(&Bfd);CHKERRQ(ierr); 2355 } 2356 } 2357 PetscFunctionReturn(0); 2358 } 2359 2360 #undef __FUNCT__ 2361 #define __FUNCT__ "SNESSetJacobian" 2362 /*@C 2363 SNESSetJacobian - Sets the function to compute Jacobian as well as the 2364 location to store the matrix. 2365 2366 Logically Collective on SNES and Mat 2367 2368 Input Parameters: 2369 + snes - the SNES context 2370 . A - Jacobian matrix 2371 . B - preconditioner matrix (usually same as the Jacobian) 2372 . func - Jacobian evaluation routine (if PETSC_NULL then SNES retains any previously set value) 2373 - ctx - [optional] user-defined context for private data for the 2374 Jacobian evaluation routine (may be PETSC_NULL) (if PETSC_NULL then SNES retains any previously set value) 2375 2376 Calling sequence of func: 2377 $ func (SNES snes,Vec x,Mat *A,Mat *B,int *flag,void *ctx); 2378 2379 + x - input vector 2380 . A - Jacobian matrix 2381 . B - preconditioner matrix, usually the same as A 2382 . flag - flag indicating information about the preconditioner matrix 2383 structure (same as flag in KSPSetOperators()), one of SAME_NONZERO_PATTERN,DIFFERENT_NONZERO_PATTERN,SAME_PRECONDITIONER 2384 - ctx - [optional] user-defined Jacobian context 2385 2386 Notes: 2387 See KSPSetOperators() for important information about setting the flag 2388 output parameter in the routine func(). Be sure to read this information! 2389 2390 The routine func() takes Mat * as the matrix arguments rather than Mat. 2391 This allows the Jacobian evaluation routine to replace A and/or B with a 2392 completely new new matrix structure (not just different matrix elements) 2393 when appropriate, for instance, if the nonzero structure is changing 2394 throughout the global iterations. 2395 2396 If the A matrix and B matrix are different you must call MatAssemblyBegin/End() on 2397 each matrix. 2398 2399 If using SNESDefaultComputeJacobianColor() to assemble a Jacobian, the ctx argument 2400 must be a MatFDColoring. 2401 2402 Other defect-correction schemes can be used by computing a different matrix in place of the Jacobian. One common 2403 example is to use the "Picard linearization" which only differentiates through the highest order parts of each term. 2404 2405 Level: beginner 2406 2407 .keywords: SNES, nonlinear, set, Jacobian, matrix 2408 2409 .seealso: KSPSetOperators(), SNESSetFunction(), MatMFFDComputeJacobian(), SNESDefaultComputeJacobianColor(), MatStructure 2410 @*/ 2411 PetscErrorCode SNESSetJacobian(SNES snes,Mat A,Mat B,PetscErrorCode (*func)(SNES,Vec,Mat*,Mat*,MatStructure*,void*),void *ctx) 2412 { 2413 PetscErrorCode ierr; 2414 DM dm; 2415 2416 PetscFunctionBegin; 2417 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2418 if (A) PetscValidHeaderSpecific(A,MAT_CLASSID,2); 2419 if (B) PetscValidHeaderSpecific(B,MAT_CLASSID,3); 2420 if (A) PetscCheckSameComm(snes,1,A,2); 2421 if (B) PetscCheckSameComm(snes,1,B,3); 2422 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 2423 ierr = DMSNESSetJacobian(dm,func,ctx);CHKERRQ(ierr); 2424 if (A) { 2425 ierr = PetscObjectReference((PetscObject)A);CHKERRQ(ierr); 2426 ierr = MatDestroy(&snes->jacobian);CHKERRQ(ierr); 2427 snes->jacobian = A; 2428 } 2429 if (B) { 2430 ierr = PetscObjectReference((PetscObject)B);CHKERRQ(ierr); 2431 ierr = MatDestroy(&snes->jacobian_pre);CHKERRQ(ierr); 2432 snes->jacobian_pre = B; 2433 } 2434 PetscFunctionReturn(0); 2435 } 2436 2437 #undef __FUNCT__ 2438 #define __FUNCT__ "SNESGetJacobian" 2439 /*@C 2440 SNESGetJacobian - Returns the Jacobian matrix and optionally the user 2441 provided context for evaluating the Jacobian. 2442 2443 Not Collective, but Mat object will be parallel if SNES object is 2444 2445 Input Parameter: 2446 . snes - the nonlinear solver context 2447 2448 Output Parameters: 2449 + A - location to stash Jacobian matrix (or PETSC_NULL) 2450 . B - location to stash preconditioner matrix (or PETSC_NULL) 2451 . func - location to put Jacobian function (or PETSC_NULL) 2452 - ctx - location to stash Jacobian ctx (or PETSC_NULL) 2453 2454 Level: advanced 2455 2456 .seealso: SNESSetJacobian(), SNESComputeJacobian() 2457 @*/ 2458 PetscErrorCode SNESGetJacobian(SNES snes,Mat *A,Mat *B,PetscErrorCode (**func)(SNES,Vec,Mat*,Mat*,MatStructure*,void*),void **ctx) 2459 { 2460 PetscErrorCode ierr; 2461 DM dm; 2462 SNESDM sdm; 2463 2464 PetscFunctionBegin; 2465 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2466 if (A) *A = snes->jacobian; 2467 if (B) *B = snes->jacobian_pre; 2468 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 2469 ierr = DMSNESGetContext(dm,&sdm);CHKERRQ(ierr); 2470 if (func) *func = sdm->computejacobian; 2471 if (ctx) *ctx = sdm->jacobianctx; 2472 PetscFunctionReturn(0); 2473 } 2474 2475 /* ----- Routines to initialize and destroy a nonlinear solver ---- */ 2476 2477 #undef __FUNCT__ 2478 #define __FUNCT__ "SNESSetUp" 2479 /*@ 2480 SNESSetUp - Sets up the internal data structures for the later use 2481 of a nonlinear solver. 2482 2483 Collective on SNES 2484 2485 Input Parameters: 2486 . snes - the SNES context 2487 2488 Notes: 2489 For basic use of the SNES solvers the user need not explicitly call 2490 SNESSetUp(), since these actions will automatically occur during 2491 the call to SNESSolve(). However, if one wishes to control this 2492 phase separately, SNESSetUp() should be called after SNESCreate() 2493 and optional routines of the form SNESSetXXX(), but before SNESSolve(). 2494 2495 Level: advanced 2496 2497 .keywords: SNES, nonlinear, setup 2498 2499 .seealso: SNESCreate(), SNESSolve(), SNESDestroy() 2500 @*/ 2501 PetscErrorCode SNESSetUp(SNES snes) 2502 { 2503 PetscErrorCode ierr; 2504 DM dm; 2505 SNESDM sdm; 2506 SNESLineSearch linesearch, pclinesearch; 2507 void *lsprectx,*lspostctx; 2508 SNESLineSearchPreCheckFunc lsprefunc; 2509 SNESLineSearchPostCheckFunc lspostfunc; 2510 PetscErrorCode (*func)(SNES,Vec,Vec,void*); 2511 Vec f,fpc; 2512 void *funcctx; 2513 PetscErrorCode (*jac)(SNES,Vec,Mat*,Mat*,MatStructure*,void*); 2514 void *jacctx,*appctx; 2515 Mat A,B; 2516 2517 PetscFunctionBegin; 2518 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2519 if (snes->setupcalled) PetscFunctionReturn(0); 2520 2521 if (!((PetscObject)snes)->type_name) { 2522 ierr = SNESSetType(snes,SNESLS);CHKERRQ(ierr); 2523 } 2524 2525 ierr = SNESGetFunction(snes,&snes->vec_func,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr); 2526 if (snes->vec_func == snes->vec_sol) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_IDN,"Solution vector cannot be function vector"); 2527 if (snes->vec_rhs == snes->vec_sol) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_IDN,"Solution vector cannot be right hand side vector"); 2528 2529 if (!snes->vec_sol_update /* && snes->vec_sol */) { 2530 ierr = VecDuplicate(snes->vec_sol,&snes->vec_sol_update);CHKERRQ(ierr); 2531 ierr = PetscLogObjectParent(snes,snes->vec_sol_update);CHKERRQ(ierr); 2532 } 2533 2534 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 2535 ierr = DMShellSetGlobalVector(snes->dm,snes->vec_sol);CHKERRQ(ierr); 2536 ierr = DMSNESSetUpLegacy(dm);CHKERRQ(ierr); /* To be removed when function routines are taken out of the DM package */ 2537 ierr = DMSNESGetContext(dm,&sdm);CHKERRQ(ierr); 2538 if (!sdm->computefunction) SETERRQ(((PetscObject)snes)->comm,PETSC_ERR_ARG_WRONGSTATE,"Must provide a residual function with SNESSetFunction(), DMSNESSetFunction(), DMDASNESSetFunctionLocal(), etc"); 2539 if (!snes->vec_func) { 2540 ierr = DMCreateGlobalVector(dm,&snes->vec_func);CHKERRQ(ierr); 2541 } 2542 2543 if (!snes->ksp) {ierr = SNESGetKSP(snes, &snes->ksp);CHKERRQ(ierr);} 2544 2545 if (!snes->linesearch) {ierr = SNESGetSNESLineSearch(snes, &snes->linesearch);} 2546 2547 if (snes->pc && (snes->pcside == PC_LEFT)) snes->mf = PETSC_TRUE; 2548 2549 if (snes->mf) { ierr = SNESSetUpMatrixFree_Private(snes, snes->mf_operator, snes->mf_version);CHKERRQ(ierr); } 2550 2551 2552 if (snes->ops->usercompute && !snes->user) { 2553 ierr = (*snes->ops->usercompute)(snes,(void**)&snes->user);CHKERRQ(ierr); 2554 } 2555 2556 if (snes->pc) { 2557 /* copy the DM over */ 2558 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 2559 ierr = SNESSetDM(snes->pc,dm);CHKERRQ(ierr); 2560 2561 /* copy the legacy SNES context not related to the DM over*/ 2562 ierr = SNESGetFunction(snes,&f,&func,&funcctx);CHKERRQ(ierr); 2563 ierr = VecDuplicate(f,&fpc);CHKERRQ(ierr); 2564 ierr = SNESSetFunction(snes->pc,fpc,func,funcctx);CHKERRQ(ierr); 2565 ierr = SNESGetJacobian(snes,&A,&B,&jac,&jacctx);CHKERRQ(ierr); 2566 ierr = SNESSetJacobian(snes->pc,A,B,jac,jacctx);CHKERRQ(ierr); 2567 ierr = SNESGetApplicationContext(snes,&appctx);CHKERRQ(ierr); 2568 ierr = SNESSetApplicationContext(snes->pc,appctx);CHKERRQ(ierr); 2569 ierr = VecDestroy(&fpc);CHKERRQ(ierr); 2570 2571 /* copy the function pointers over */ 2572 ierr = PetscObjectCopyFortranFunctionPointers((PetscObject)snes,(PetscObject)snes->pc);CHKERRQ(ierr); 2573 2574 /* default to 1 iteration */ 2575 ierr = SNESSetTolerances(snes->pc, 0.0, 0.0, 0.0, 1, snes->pc->max_funcs);CHKERRQ(ierr); 2576 ierr = SNESSetNormType(snes->pc, SNES_NORM_FINAL_ONLY);CHKERRQ(ierr); 2577 ierr = SNESSetFromOptions(snes->pc);CHKERRQ(ierr); 2578 2579 /* copy the line search context over */ 2580 ierr = SNESGetSNESLineSearch(snes,&linesearch);CHKERRQ(ierr); 2581 ierr = SNESGetSNESLineSearch(snes->pc,&pclinesearch);CHKERRQ(ierr); 2582 ierr = SNESLineSearchGetPreCheck(linesearch,&lsprefunc,&lsprectx);CHKERRQ(ierr); 2583 ierr = SNESLineSearchGetPostCheck(linesearch,&lspostfunc,&lspostctx);CHKERRQ(ierr); 2584 ierr = SNESLineSearchSetPreCheck(pclinesearch,lsprefunc,lsprectx);CHKERRQ(ierr); 2585 ierr = SNESLineSearchSetPostCheck(pclinesearch,lspostfunc,lspostctx);CHKERRQ(ierr); 2586 ierr = PetscObjectCopyFortranFunctionPointers((PetscObject)linesearch, (PetscObject)pclinesearch);CHKERRQ(ierr); 2587 } 2588 2589 if (snes->ops->setup) { 2590 ierr = (*snes->ops->setup)(snes);CHKERRQ(ierr); 2591 } 2592 2593 snes->setupcalled = PETSC_TRUE; 2594 PetscFunctionReturn(0); 2595 } 2596 2597 #undef __FUNCT__ 2598 #define __FUNCT__ "SNESReset" 2599 /*@ 2600 SNESReset - Resets a SNES context to the snessetupcalled = 0 state and removes any allocated Vecs and Mats 2601 2602 Collective on SNES 2603 2604 Input Parameter: 2605 . snes - iterative context obtained from SNESCreate() 2606 2607 Level: intermediate 2608 2609 Notes: Also calls the application context destroy routine set with SNESSetComputeApplicationContext() 2610 2611 .keywords: SNES, destroy 2612 2613 .seealso: SNESCreate(), SNESSetUp(), SNESSolve() 2614 @*/ 2615 PetscErrorCode SNESReset(SNES snes) 2616 { 2617 PetscErrorCode ierr; 2618 2619 PetscFunctionBegin; 2620 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2621 if (snes->ops->userdestroy && snes->user) { 2622 ierr = (*snes->ops->userdestroy)((void**)&snes->user);CHKERRQ(ierr); 2623 snes->user = PETSC_NULL; 2624 } 2625 if (snes->pc) { 2626 ierr = SNESReset(snes->pc);CHKERRQ(ierr); 2627 } 2628 2629 if (snes->ops->reset) { 2630 ierr = (*snes->ops->reset)(snes);CHKERRQ(ierr); 2631 } 2632 if (snes->ksp) { 2633 ierr = KSPReset(snes->ksp);CHKERRQ(ierr); 2634 } 2635 2636 if (snes->linesearch) { 2637 ierr = SNESLineSearchReset(snes->linesearch);CHKERRQ(ierr); 2638 } 2639 2640 ierr = VecDestroy(&snes->vec_rhs);CHKERRQ(ierr); 2641 ierr = VecDestroy(&snes->vec_sol);CHKERRQ(ierr); 2642 ierr = VecDestroy(&snes->vec_sol_update);CHKERRQ(ierr); 2643 ierr = VecDestroy(&snes->vec_func);CHKERRQ(ierr); 2644 ierr = MatDestroy(&snes->jacobian);CHKERRQ(ierr); 2645 ierr = MatDestroy(&snes->jacobian_pre);CHKERRQ(ierr); 2646 ierr = VecDestroyVecs(snes->nwork,&snes->work);CHKERRQ(ierr); 2647 ierr = VecDestroyVecs(snes->nvwork,&snes->vwork);CHKERRQ(ierr); 2648 snes->nwork = snes->nvwork = 0; 2649 snes->setupcalled = PETSC_FALSE; 2650 PetscFunctionReturn(0); 2651 } 2652 2653 #undef __FUNCT__ 2654 #define __FUNCT__ "SNESDestroy" 2655 /*@ 2656 SNESDestroy - Destroys the nonlinear solver context that was created 2657 with SNESCreate(). 2658 2659 Collective on SNES 2660 2661 Input Parameter: 2662 . snes - the SNES context 2663 2664 Level: beginner 2665 2666 .keywords: SNES, nonlinear, destroy 2667 2668 .seealso: SNESCreate(), SNESSolve() 2669 @*/ 2670 PetscErrorCode SNESDestroy(SNES *snes) 2671 { 2672 PetscErrorCode ierr; 2673 2674 PetscFunctionBegin; 2675 if (!*snes) PetscFunctionReturn(0); 2676 PetscValidHeaderSpecific((*snes),SNES_CLASSID,1); 2677 if (--((PetscObject)(*snes))->refct > 0) {*snes = 0; PetscFunctionReturn(0);} 2678 2679 ierr = SNESReset((*snes));CHKERRQ(ierr); 2680 ierr = SNESDestroy(&(*snes)->pc);CHKERRQ(ierr); 2681 2682 /* if memory was published with AMS then destroy it */ 2683 ierr = PetscObjectDepublish((*snes));CHKERRQ(ierr); 2684 if ((*snes)->ops->destroy) {ierr = (*((*snes))->ops->destroy)((*snes));CHKERRQ(ierr);} 2685 2686 ierr = DMDestroy(&(*snes)->dm);CHKERRQ(ierr); 2687 ierr = KSPDestroy(&(*snes)->ksp);CHKERRQ(ierr); 2688 ierr = SNESLineSearchDestroy(&(*snes)->linesearch);CHKERRQ(ierr); 2689 2690 ierr = PetscFree((*snes)->kspconvctx);CHKERRQ(ierr); 2691 if ((*snes)->ops->convergeddestroy) { 2692 ierr = (*(*snes)->ops->convergeddestroy)((*snes)->cnvP);CHKERRQ(ierr); 2693 } 2694 if ((*snes)->conv_malloc) { 2695 ierr = PetscFree((*snes)->conv_hist);CHKERRQ(ierr); 2696 ierr = PetscFree((*snes)->conv_hist_its);CHKERRQ(ierr); 2697 } 2698 ierr = SNESMonitorCancel((*snes));CHKERRQ(ierr); 2699 ierr = PetscHeaderDestroy(snes);CHKERRQ(ierr); 2700 PetscFunctionReturn(0); 2701 } 2702 2703 /* ----------- Routines to set solver parameters ---------- */ 2704 2705 #undef __FUNCT__ 2706 #define __FUNCT__ "SNESSetLagPreconditioner" 2707 /*@ 2708 SNESSetLagPreconditioner - Determines when the preconditioner is rebuilt in the nonlinear solve. 2709 2710 Logically Collective on SNES 2711 2712 Input Parameters: 2713 + snes - the SNES context 2714 - lag - -1 indicates NEVER rebuild, 1 means rebuild every time the Jacobian is computed within a single nonlinear solve, 2 means every second time 2715 the Jacobian is built etc. -2 indicates rebuild preconditioner at next chance but then never rebuild after that 2716 2717 Options Database Keys: 2718 . -snes_lag_preconditioner <lag> 2719 2720 Notes: 2721 The default is 1 2722 The preconditioner is ALWAYS built in the first iteration of a nonlinear solve unless lag is -1 2723 If -1 is used before the very first nonlinear solve the preconditioner is still built because there is no previous preconditioner to use 2724 2725 Level: intermediate 2726 2727 .keywords: SNES, nonlinear, set, convergence, tolerances 2728 2729 .seealso: SNESSetTrustRegionTolerance(), SNESGetLagPreconditioner(), SNESSetLagJacobian(), SNESGetLagJacobian() 2730 2731 @*/ 2732 PetscErrorCode SNESSetLagPreconditioner(SNES snes,PetscInt lag) 2733 { 2734 PetscFunctionBegin; 2735 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2736 if (lag < -2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Lag must be -2, -1, 1 or greater"); 2737 if (!lag) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Lag cannot be 0"); 2738 PetscValidLogicalCollectiveInt(snes,lag,2); 2739 snes->lagpreconditioner = lag; 2740 PetscFunctionReturn(0); 2741 } 2742 2743 #undef __FUNCT__ 2744 #define __FUNCT__ "SNESSetGridSequence" 2745 /*@ 2746 SNESSetGridSequence - sets the number of steps of grid sequencing that SNES does 2747 2748 Logically Collective on SNES 2749 2750 Input Parameters: 2751 + snes - the SNES context 2752 - steps - the number of refinements to do, defaults to 0 2753 2754 Options Database Keys: 2755 . -snes_grid_sequence <steps> 2756 2757 Level: intermediate 2758 2759 Notes: 2760 Use SNESGetSolution() to extract the fine grid solution after grid sequencing. 2761 2762 .keywords: SNES, nonlinear, set, convergence, tolerances 2763 2764 .seealso: SNESSetTrustRegionTolerance(), SNESGetLagPreconditioner(), SNESSetLagJacobian(), SNESGetLagJacobian() 2765 2766 @*/ 2767 PetscErrorCode SNESSetGridSequence(SNES snes,PetscInt steps) 2768 { 2769 PetscFunctionBegin; 2770 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2771 PetscValidLogicalCollectiveInt(snes,steps,2); 2772 snes->gridsequence = steps; 2773 PetscFunctionReturn(0); 2774 } 2775 2776 #undef __FUNCT__ 2777 #define __FUNCT__ "SNESGetLagPreconditioner" 2778 /*@ 2779 SNESGetLagPreconditioner - Indicates how often the preconditioner is rebuilt 2780 2781 Not Collective 2782 2783 Input Parameter: 2784 . snes - the SNES context 2785 2786 Output Parameter: 2787 . lag - -1 indicates NEVER rebuild, 1 means rebuild every time the Jacobian is computed within a single nonlinear solve, 2 means every second time 2788 the Jacobian is built etc. -2 indicates rebuild preconditioner at next chance but then never rebuild after that 2789 2790 Options Database Keys: 2791 . -snes_lag_preconditioner <lag> 2792 2793 Notes: 2794 The default is 1 2795 The preconditioner is ALWAYS built in the first iteration of a nonlinear solve unless lag is -1 2796 2797 Level: intermediate 2798 2799 .keywords: SNES, nonlinear, set, convergence, tolerances 2800 2801 .seealso: SNESSetTrustRegionTolerance(), SNESSetLagPreconditioner() 2802 2803 @*/ 2804 PetscErrorCode SNESGetLagPreconditioner(SNES snes,PetscInt *lag) 2805 { 2806 PetscFunctionBegin; 2807 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2808 *lag = snes->lagpreconditioner; 2809 PetscFunctionReturn(0); 2810 } 2811 2812 #undef __FUNCT__ 2813 #define __FUNCT__ "SNESSetLagJacobian" 2814 /*@ 2815 SNESSetLagJacobian - Determines when the Jacobian is rebuilt in the nonlinear solve. See SNESSetLagPreconditioner() for determining how 2816 often the preconditioner is rebuilt. 2817 2818 Logically Collective on SNES 2819 2820 Input Parameters: 2821 + snes - the SNES context 2822 - lag - -1 indicates NEVER rebuild, 1 means rebuild every time the Jacobian is computed within a single nonlinear solve, 2 means every second time 2823 the Jacobian is built etc. -2 means rebuild at next chance but then never again 2824 2825 Options Database Keys: 2826 . -snes_lag_jacobian <lag> 2827 2828 Notes: 2829 The default is 1 2830 The Jacobian is ALWAYS built in the first iteration of a nonlinear solve unless lag is -1 2831 If -1 is used before the very first nonlinear solve the CODE WILL FAIL! because no Jacobian is used, use -2 to indicate you want it recomputed 2832 at the next Newton step but never again (unless it is reset to another value) 2833 2834 Level: intermediate 2835 2836 .keywords: SNES, nonlinear, set, convergence, tolerances 2837 2838 .seealso: SNESSetTrustRegionTolerance(), SNESGetLagPreconditioner(), SNESSetLagPreconditioner(), SNESGetLagJacobian() 2839 2840 @*/ 2841 PetscErrorCode SNESSetLagJacobian(SNES snes,PetscInt lag) 2842 { 2843 PetscFunctionBegin; 2844 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2845 if (lag < -2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Lag must be -2, -1, 1 or greater"); 2846 if (!lag) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Lag cannot be 0"); 2847 PetscValidLogicalCollectiveInt(snes,lag,2); 2848 snes->lagjacobian = lag; 2849 PetscFunctionReturn(0); 2850 } 2851 2852 #undef __FUNCT__ 2853 #define __FUNCT__ "SNESGetLagJacobian" 2854 /*@ 2855 SNESGetLagJacobian - Indicates how often the Jacobian is rebuilt. See SNESGetLagPreconditioner() to determine when the preconditioner is rebuilt 2856 2857 Not Collective 2858 2859 Input Parameter: 2860 . snes - the SNES context 2861 2862 Output Parameter: 2863 . lag - -1 indicates NEVER rebuild, 1 means rebuild every time the Jacobian is computed within a single nonlinear solve, 2 means every second time 2864 the Jacobian is built etc. 2865 2866 Options Database Keys: 2867 . -snes_lag_jacobian <lag> 2868 2869 Notes: 2870 The default is 1 2871 The jacobian is ALWAYS built in the first iteration of a nonlinear solve unless lag is -1 2872 2873 Level: intermediate 2874 2875 .keywords: SNES, nonlinear, set, convergence, tolerances 2876 2877 .seealso: SNESSetTrustRegionTolerance(), SNESSetLagJacobian(), SNESSetLagPreconditioner(), SNESGetLagPreconditioner() 2878 2879 @*/ 2880 PetscErrorCode SNESGetLagJacobian(SNES snes,PetscInt *lag) 2881 { 2882 PetscFunctionBegin; 2883 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2884 *lag = snes->lagjacobian; 2885 PetscFunctionReturn(0); 2886 } 2887 2888 #undef __FUNCT__ 2889 #define __FUNCT__ "SNESSetTolerances" 2890 /*@ 2891 SNESSetTolerances - Sets various parameters used in convergence tests. 2892 2893 Logically Collective on SNES 2894 2895 Input Parameters: 2896 + snes - the SNES context 2897 . abstol - absolute convergence tolerance 2898 . rtol - relative convergence tolerance 2899 . stol - convergence tolerance in terms of the norm 2900 of the change in the solution between steps 2901 . maxit - maximum number of iterations 2902 - maxf - maximum number of function evaluations 2903 2904 Options Database Keys: 2905 + -snes_atol <abstol> - Sets abstol 2906 . -snes_rtol <rtol> - Sets rtol 2907 . -snes_stol <stol> - Sets stol 2908 . -snes_max_it <maxit> - Sets maxit 2909 - -snes_max_funcs <maxf> - Sets maxf 2910 2911 Notes: 2912 The default maximum number of iterations is 50. 2913 The default maximum number of function evaluations is 1000. 2914 2915 Level: intermediate 2916 2917 .keywords: SNES, nonlinear, set, convergence, tolerances 2918 2919 .seealso: SNESSetTrustRegionTolerance() 2920 @*/ 2921 PetscErrorCode SNESSetTolerances(SNES snes,PetscReal abstol,PetscReal rtol,PetscReal stol,PetscInt maxit,PetscInt maxf) 2922 { 2923 PetscFunctionBegin; 2924 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2925 PetscValidLogicalCollectiveReal(snes,abstol,2); 2926 PetscValidLogicalCollectiveReal(snes,rtol,3); 2927 PetscValidLogicalCollectiveReal(snes,stol,4); 2928 PetscValidLogicalCollectiveInt(snes,maxit,5); 2929 PetscValidLogicalCollectiveInt(snes,maxf,6); 2930 2931 if (abstol != PETSC_DEFAULT) { 2932 if (abstol < 0.0) SETERRQ1(((PetscObject)snes)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Absolute tolerance %G must be non-negative",abstol); 2933 snes->abstol = abstol; 2934 } 2935 if (rtol != PETSC_DEFAULT) { 2936 if (rtol < 0.0 || 1.0 <= rtol) SETERRQ1(((PetscObject)snes)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Relative tolerance %G must be non-negative and less than 1.0",rtol); 2937 snes->rtol = rtol; 2938 } 2939 if (stol != PETSC_DEFAULT) { 2940 if (stol < 0.0) SETERRQ1(((PetscObject)snes)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Step tolerance %G must be non-negative",stol); 2941 snes->stol = stol; 2942 } 2943 if (maxit != PETSC_DEFAULT) { 2944 if (maxit < 0) SETERRQ1(((PetscObject)snes)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Maximum number of iterations %D must be non-negative",maxit); 2945 snes->max_its = maxit; 2946 } 2947 if (maxf != PETSC_DEFAULT) { 2948 if (maxf < 0) SETERRQ1(((PetscObject)snes)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Maximum number of function evaluations %D must be non-negative",maxf); 2949 snes->max_funcs = maxf; 2950 } 2951 snes->tolerancesset = PETSC_TRUE; 2952 PetscFunctionReturn(0); 2953 } 2954 2955 #undef __FUNCT__ 2956 #define __FUNCT__ "SNESGetTolerances" 2957 /*@ 2958 SNESGetTolerances - Gets various parameters used in convergence tests. 2959 2960 Not Collective 2961 2962 Input Parameters: 2963 + snes - the SNES context 2964 . atol - absolute convergence tolerance 2965 . rtol - relative convergence tolerance 2966 . stol - convergence tolerance in terms of the norm 2967 of the change in the solution between steps 2968 . maxit - maximum number of iterations 2969 - maxf - maximum number of function evaluations 2970 2971 Notes: 2972 The user can specify PETSC_NULL for any parameter that is not needed. 2973 2974 Level: intermediate 2975 2976 .keywords: SNES, nonlinear, get, convergence, tolerances 2977 2978 .seealso: SNESSetTolerances() 2979 @*/ 2980 PetscErrorCode SNESGetTolerances(SNES snes,PetscReal *atol,PetscReal *rtol,PetscReal *stol,PetscInt *maxit,PetscInt *maxf) 2981 { 2982 PetscFunctionBegin; 2983 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2984 if (atol) *atol = snes->abstol; 2985 if (rtol) *rtol = snes->rtol; 2986 if (stol) *stol = snes->stol; 2987 if (maxit) *maxit = snes->max_its; 2988 if (maxf) *maxf = snes->max_funcs; 2989 PetscFunctionReturn(0); 2990 } 2991 2992 #undef __FUNCT__ 2993 #define __FUNCT__ "SNESSetTrustRegionTolerance" 2994 /*@ 2995 SNESSetTrustRegionTolerance - Sets the trust region parameter tolerance. 2996 2997 Logically Collective on SNES 2998 2999 Input Parameters: 3000 + snes - the SNES context 3001 - tol - tolerance 3002 3003 Options Database Key: 3004 . -snes_trtol <tol> - Sets tol 3005 3006 Level: intermediate 3007 3008 .keywords: SNES, nonlinear, set, trust region, tolerance 3009 3010 .seealso: SNESSetTolerances() 3011 @*/ 3012 PetscErrorCode SNESSetTrustRegionTolerance(SNES snes,PetscReal tol) 3013 { 3014 PetscFunctionBegin; 3015 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3016 PetscValidLogicalCollectiveReal(snes,tol,2); 3017 snes->deltatol = tol; 3018 PetscFunctionReturn(0); 3019 } 3020 3021 /* 3022 Duplicate the lg monitors for SNES from KSP; for some reason with 3023 dynamic libraries things don't work under Sun4 if we just use 3024 macros instead of functions 3025 */ 3026 #undef __FUNCT__ 3027 #define __FUNCT__ "SNESMonitorLGResidualNorm" 3028 PetscErrorCode SNESMonitorLGResidualNorm(SNES snes,PetscInt it,PetscReal norm,void *ctx) 3029 { 3030 PetscErrorCode ierr; 3031 3032 PetscFunctionBegin; 3033 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3034 ierr = KSPMonitorLGResidualNorm((KSP)snes,it,norm,ctx);CHKERRQ(ierr); 3035 PetscFunctionReturn(0); 3036 } 3037 3038 #undef __FUNCT__ 3039 #define __FUNCT__ "SNESMonitorLGCreate" 3040 PetscErrorCode SNESMonitorLGCreate(const char host[],const char label[],int x,int y,int m,int n,PetscDrawLG *draw) 3041 { 3042 PetscErrorCode ierr; 3043 3044 PetscFunctionBegin; 3045 ierr = KSPMonitorLGResidualNormCreate(host,label,x,y,m,n,draw);CHKERRQ(ierr); 3046 PetscFunctionReturn(0); 3047 } 3048 3049 #undef __FUNCT__ 3050 #define __FUNCT__ "SNESMonitorLGDestroy" 3051 PetscErrorCode SNESMonitorLGDestroy(PetscDrawLG *draw) 3052 { 3053 PetscErrorCode ierr; 3054 3055 PetscFunctionBegin; 3056 ierr = KSPMonitorLGResidualNormDestroy(draw);CHKERRQ(ierr); 3057 PetscFunctionReturn(0); 3058 } 3059 3060 extern PetscErrorCode SNESMonitorRange_Private(SNES,PetscInt,PetscReal*); 3061 #undef __FUNCT__ 3062 #define __FUNCT__ "SNESMonitorLGRange" 3063 PetscErrorCode SNESMonitorLGRange(SNES snes,PetscInt n,PetscReal rnorm,void *monctx) 3064 { 3065 PetscDrawLG lg; 3066 PetscErrorCode ierr; 3067 PetscReal x,y,per; 3068 PetscViewer v = (PetscViewer)monctx; 3069 static PetscReal prev; /* should be in the context */ 3070 PetscDraw draw; 3071 3072 PetscFunctionBegin; 3073 ierr = PetscViewerDrawGetDrawLG(v,0,&lg);CHKERRQ(ierr); 3074 if (!n) {ierr = PetscDrawLGReset(lg);CHKERRQ(ierr);} 3075 ierr = PetscDrawLGGetDraw(lg,&draw);CHKERRQ(ierr); 3076 ierr = PetscDrawSetTitle(draw,"Residual norm");CHKERRQ(ierr); 3077 x = (PetscReal) n; 3078 if (rnorm > 0.0) y = log10(rnorm); else y = -15.0; 3079 ierr = PetscDrawLGAddPoint(lg,&x,&y);CHKERRQ(ierr); 3080 if (n < 20 || !(n % 5)) { 3081 ierr = PetscDrawLGDraw(lg);CHKERRQ(ierr); 3082 } 3083 3084 ierr = PetscViewerDrawGetDrawLG(v,1,&lg);CHKERRQ(ierr); 3085 if (!n) {ierr = PetscDrawLGReset(lg);CHKERRQ(ierr);} 3086 ierr = PetscDrawLGGetDraw(lg,&draw);CHKERRQ(ierr); 3087 ierr = PetscDrawSetTitle(draw,"% elemts > .2*max elemt");CHKERRQ(ierr); 3088 ierr = SNESMonitorRange_Private(snes,n,&per);CHKERRQ(ierr); 3089 x = (PetscReal) n; 3090 y = 100.0*per; 3091 ierr = PetscDrawLGAddPoint(lg,&x,&y);CHKERRQ(ierr); 3092 if (n < 20 || !(n % 5)) { 3093 ierr = PetscDrawLGDraw(lg);CHKERRQ(ierr); 3094 } 3095 3096 ierr = PetscViewerDrawGetDrawLG(v,2,&lg);CHKERRQ(ierr); 3097 if (!n) {prev = rnorm;ierr = PetscDrawLGReset(lg);CHKERRQ(ierr);} 3098 ierr = PetscDrawLGGetDraw(lg,&draw);CHKERRQ(ierr); 3099 ierr = PetscDrawSetTitle(draw,"(norm -oldnorm)/oldnorm");CHKERRQ(ierr); 3100 x = (PetscReal) n; 3101 y = (prev - rnorm)/prev; 3102 ierr = PetscDrawLGAddPoint(lg,&x,&y);CHKERRQ(ierr); 3103 if (n < 20 || !(n % 5)) { 3104 ierr = PetscDrawLGDraw(lg);CHKERRQ(ierr); 3105 } 3106 3107 ierr = PetscViewerDrawGetDrawLG(v,3,&lg);CHKERRQ(ierr); 3108 if (!n) {ierr = PetscDrawLGReset(lg);CHKERRQ(ierr);} 3109 ierr = PetscDrawLGGetDraw(lg,&draw);CHKERRQ(ierr); 3110 ierr = PetscDrawSetTitle(draw,"(norm -oldnorm)/oldnorm*(% > .2 max)");CHKERRQ(ierr); 3111 x = (PetscReal) n; 3112 y = (prev - rnorm)/(prev*per); 3113 if (n > 2) { /*skip initial crazy value */ 3114 ierr = PetscDrawLGAddPoint(lg,&x,&y);CHKERRQ(ierr); 3115 } 3116 if (n < 20 || !(n % 5)) { 3117 ierr = PetscDrawLGDraw(lg);CHKERRQ(ierr); 3118 } 3119 prev = rnorm; 3120 PetscFunctionReturn(0); 3121 } 3122 3123 #undef __FUNCT__ 3124 #define __FUNCT__ "SNESMonitor" 3125 /*@ 3126 SNESMonitor - runs the user provided monitor routines, if they exist 3127 3128 Collective on SNES 3129 3130 Input Parameters: 3131 + snes - nonlinear solver context obtained from SNESCreate() 3132 . iter - iteration number 3133 - rnorm - relative norm of the residual 3134 3135 Notes: 3136 This routine is called by the SNES implementations. 3137 It does not typically need to be called by the user. 3138 3139 Level: developer 3140 3141 .seealso: SNESMonitorSet() 3142 @*/ 3143 PetscErrorCode SNESMonitor(SNES snes,PetscInt iter,PetscReal rnorm) 3144 { 3145 PetscErrorCode ierr; 3146 PetscInt i,n = snes->numbermonitors; 3147 3148 PetscFunctionBegin; 3149 for (i=0; i<n; i++) { 3150 ierr = (*snes->monitor[i])(snes,iter,rnorm,snes->monitorcontext[i]);CHKERRQ(ierr); 3151 } 3152 PetscFunctionReturn(0); 3153 } 3154 3155 /* ------------ Routines to set performance monitoring options ----------- */ 3156 3157 #undef __FUNCT__ 3158 #define __FUNCT__ "SNESMonitorSet" 3159 /*@C 3160 SNESMonitorSet - Sets an ADDITIONAL function that is to be used at every 3161 iteration of the nonlinear solver to display the iteration's 3162 progress. 3163 3164 Logically Collective on SNES 3165 3166 Input Parameters: 3167 + snes - the SNES context 3168 . func - monitoring routine 3169 . mctx - [optional] user-defined context for private data for the 3170 monitor routine (use PETSC_NULL if no context is desired) 3171 - monitordestroy - [optional] routine that frees monitor context 3172 (may be PETSC_NULL) 3173 3174 Calling sequence of func: 3175 $ int func(SNES snes,PetscInt its, PetscReal norm,void *mctx) 3176 3177 + snes - the SNES context 3178 . its - iteration number 3179 . norm - 2-norm function value (may be estimated) 3180 - mctx - [optional] monitoring context 3181 3182 Options Database Keys: 3183 + -snes_monitor - sets SNESMonitorDefault() 3184 . -snes_monitor_lg_residualnorm - sets line graph monitor, 3185 uses SNESMonitorLGCreate() 3186 - -snes_monitor_cancel - cancels all monitors that have 3187 been hardwired into a code by 3188 calls to SNESMonitorSet(), but 3189 does not cancel those set via 3190 the options database. 3191 3192 Notes: 3193 Several different monitoring routines may be set by calling 3194 SNESMonitorSet() multiple times; all will be called in the 3195 order in which they were set. 3196 3197 Fortran notes: Only a single monitor function can be set for each SNES object 3198 3199 Level: intermediate 3200 3201 .keywords: SNES, nonlinear, set, monitor 3202 3203 .seealso: SNESMonitorDefault(), SNESMonitorCancel() 3204 @*/ 3205 PetscErrorCode SNESMonitorSet(SNES snes,PetscErrorCode (*monitor)(SNES,PetscInt,PetscReal,void*),void *mctx,PetscErrorCode (*monitordestroy)(void**)) 3206 { 3207 PetscInt i; 3208 PetscErrorCode ierr; 3209 3210 PetscFunctionBegin; 3211 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3212 if (snes->numbermonitors >= MAXSNESMONITORS) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Too many monitors set"); 3213 for (i=0; i<snes->numbermonitors;i++) { 3214 if (monitor == snes->monitor[i] && monitordestroy == snes->monitordestroy[i] && mctx == snes->monitorcontext[i]) { 3215 if (monitordestroy) { 3216 ierr = (*monitordestroy)(&mctx);CHKERRQ(ierr); 3217 } 3218 PetscFunctionReturn(0); 3219 } 3220 } 3221 snes->monitor[snes->numbermonitors] = monitor; 3222 snes->monitordestroy[snes->numbermonitors] = monitordestroy; 3223 snes->monitorcontext[snes->numbermonitors++] = (void*)mctx; 3224 PetscFunctionReturn(0); 3225 } 3226 3227 #undef __FUNCT__ 3228 #define __FUNCT__ "SNESMonitorCancel" 3229 /*@C 3230 SNESMonitorCancel - Clears all the monitor functions for a SNES object. 3231 3232 Logically Collective on SNES 3233 3234 Input Parameters: 3235 . snes - the SNES context 3236 3237 Options Database Key: 3238 . -snes_monitor_cancel - cancels all monitors that have been hardwired 3239 into a code by calls to SNESMonitorSet(), but does not cancel those 3240 set via the options database 3241 3242 Notes: 3243 There is no way to clear one specific monitor from a SNES object. 3244 3245 Level: intermediate 3246 3247 .keywords: SNES, nonlinear, set, monitor 3248 3249 .seealso: SNESMonitorDefault(), SNESMonitorSet() 3250 @*/ 3251 PetscErrorCode SNESMonitorCancel(SNES snes) 3252 { 3253 PetscErrorCode ierr; 3254 PetscInt i; 3255 3256 PetscFunctionBegin; 3257 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3258 for (i=0; i<snes->numbermonitors; i++) { 3259 if (snes->monitordestroy[i]) { 3260 ierr = (*snes->monitordestroy[i])(&snes->monitorcontext[i]);CHKERRQ(ierr); 3261 } 3262 } 3263 snes->numbermonitors = 0; 3264 PetscFunctionReturn(0); 3265 } 3266 3267 #undef __FUNCT__ 3268 #define __FUNCT__ "SNESSetConvergenceTest" 3269 /*@C 3270 SNESSetConvergenceTest - Sets the function that is to be used 3271 to test for convergence of the nonlinear iterative solution. 3272 3273 Logically Collective on SNES 3274 3275 Input Parameters: 3276 + snes - the SNES context 3277 . func - routine to test for convergence 3278 . cctx - [optional] context for private data for the convergence routine (may be PETSC_NULL) 3279 - destroy - [optional] destructor for the context (may be PETSC_NULL; PETSC_NULL_FUNCTION in Fortran) 3280 3281 Calling sequence of func: 3282 $ PetscErrorCode func (SNES snes,PetscInt it,PetscReal xnorm,PetscReal gnorm,PetscReal f,SNESConvergedReason *reason,void *cctx) 3283 3284 + snes - the SNES context 3285 . it - current iteration (0 is the first and is before any Newton step) 3286 . cctx - [optional] convergence context 3287 . reason - reason for convergence/divergence 3288 . xnorm - 2-norm of current iterate 3289 . gnorm - 2-norm of current step 3290 - f - 2-norm of function 3291 3292 Level: advanced 3293 3294 .keywords: SNES, nonlinear, set, convergence, test 3295 3296 .seealso: SNESDefaultConverged(), SNESSkipConverged() 3297 @*/ 3298 PetscErrorCode SNESSetConvergenceTest(SNES snes,PetscErrorCode (*func)(SNES,PetscInt,PetscReal,PetscReal,PetscReal,SNESConvergedReason*,void*),void *cctx,PetscErrorCode (*destroy)(void*)) 3299 { 3300 PetscErrorCode ierr; 3301 3302 PetscFunctionBegin; 3303 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3304 if (!func) func = SNESSkipConverged; 3305 if (snes->ops->convergeddestroy) { 3306 ierr = (*snes->ops->convergeddestroy)(snes->cnvP);CHKERRQ(ierr); 3307 } 3308 snes->ops->converged = func; 3309 snes->ops->convergeddestroy = destroy; 3310 snes->cnvP = cctx; 3311 PetscFunctionReturn(0); 3312 } 3313 3314 #undef __FUNCT__ 3315 #define __FUNCT__ "SNESGetConvergedReason" 3316 /*@ 3317 SNESGetConvergedReason - Gets the reason the SNES iteration was stopped. 3318 3319 Not Collective 3320 3321 Input Parameter: 3322 . snes - the SNES context 3323 3324 Output Parameter: 3325 . reason - negative value indicates diverged, positive value converged, see SNESConvergedReason or the 3326 manual pages for the individual convergence tests for complete lists 3327 3328 Level: intermediate 3329 3330 Notes: Can only be called after the call the SNESSolve() is complete. 3331 3332 .keywords: SNES, nonlinear, set, convergence, test 3333 3334 .seealso: SNESSetConvergenceTest(), SNESConvergedReason 3335 @*/ 3336 PetscErrorCode SNESGetConvergedReason(SNES snes,SNESConvergedReason *reason) 3337 { 3338 PetscFunctionBegin; 3339 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3340 PetscValidPointer(reason,2); 3341 *reason = snes->reason; 3342 PetscFunctionReturn(0); 3343 } 3344 3345 #undef __FUNCT__ 3346 #define __FUNCT__ "SNESSetConvergenceHistory" 3347 /*@ 3348 SNESSetConvergenceHistory - Sets the array used to hold the convergence history. 3349 3350 Logically Collective on SNES 3351 3352 Input Parameters: 3353 + snes - iterative context obtained from SNESCreate() 3354 . a - array to hold history, this array will contain the function norms computed at each step 3355 . its - integer array holds the number of linear iterations for each solve. 3356 . na - size of a and its 3357 - reset - PETSC_TRUE indicates each new nonlinear solve resets the history counter to zero, 3358 else it continues storing new values for new nonlinear solves after the old ones 3359 3360 Notes: 3361 If 'a' and 'its' are PETSC_NULL then space is allocated for the history. If 'na' PETSC_DECIDE or PETSC_DEFAULT then a 3362 default array of length 10000 is allocated. 3363 3364 This routine is useful, e.g., when running a code for purposes 3365 of accurate performance monitoring, when no I/O should be done 3366 during the section of code that is being timed. 3367 3368 Level: intermediate 3369 3370 .keywords: SNES, set, convergence, history 3371 3372 .seealso: SNESGetConvergenceHistory() 3373 3374 @*/ 3375 PetscErrorCode SNESSetConvergenceHistory(SNES snes,PetscReal a[],PetscInt its[],PetscInt na,PetscBool reset) 3376 { 3377 PetscErrorCode ierr; 3378 3379 PetscFunctionBegin; 3380 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3381 if (na) PetscValidScalarPointer(a,2); 3382 if (its) PetscValidIntPointer(its,3); 3383 if (na == PETSC_DECIDE || na == PETSC_DEFAULT || !a) { 3384 if (na == PETSC_DECIDE || na == PETSC_DEFAULT) na = 1000; 3385 ierr = PetscMalloc(na*sizeof(PetscReal),&a);CHKERRQ(ierr); 3386 ierr = PetscMalloc(na*sizeof(PetscInt),&its);CHKERRQ(ierr); 3387 snes->conv_malloc = PETSC_TRUE; 3388 } 3389 snes->conv_hist = a; 3390 snes->conv_hist_its = its; 3391 snes->conv_hist_max = na; 3392 snes->conv_hist_len = 0; 3393 snes->conv_hist_reset = reset; 3394 PetscFunctionReturn(0); 3395 } 3396 3397 #if defined(PETSC_HAVE_MATLAB_ENGINE) 3398 #include <engine.h> /* MATLAB include file */ 3399 #include <mex.h> /* MATLAB include file */ 3400 EXTERN_C_BEGIN 3401 #undef __FUNCT__ 3402 #define __FUNCT__ "SNESGetConvergenceHistoryMatlab" 3403 mxArray *SNESGetConvergenceHistoryMatlab(SNES snes) 3404 { 3405 mxArray *mat; 3406 PetscInt i; 3407 PetscReal *ar; 3408 3409 PetscFunctionBegin; 3410 mat = mxCreateDoubleMatrix(snes->conv_hist_len,1,mxREAL); 3411 ar = (PetscReal*) mxGetData(mat); 3412 for (i=0; i<snes->conv_hist_len; i++) { 3413 ar[i] = snes->conv_hist[i]; 3414 } 3415 PetscFunctionReturn(mat); 3416 } 3417 EXTERN_C_END 3418 #endif 3419 3420 3421 #undef __FUNCT__ 3422 #define __FUNCT__ "SNESGetConvergenceHistory" 3423 /*@C 3424 SNESGetConvergenceHistory - Gets the array used to hold the convergence history. 3425 3426 Not Collective 3427 3428 Input Parameter: 3429 . snes - iterative context obtained from SNESCreate() 3430 3431 Output Parameters: 3432 . a - array to hold history 3433 . its - integer array holds the number of linear iterations (or 3434 negative if not converged) for each solve. 3435 - na - size of a and its 3436 3437 Notes: 3438 The calling sequence for this routine in Fortran is 3439 $ call SNESGetConvergenceHistory(SNES snes, integer na, integer ierr) 3440 3441 This routine is useful, e.g., when running a code for purposes 3442 of accurate performance monitoring, when no I/O should be done 3443 during the section of code that is being timed. 3444 3445 Level: intermediate 3446 3447 .keywords: SNES, get, convergence, history 3448 3449 .seealso: SNESSetConvergencHistory() 3450 3451 @*/ 3452 PetscErrorCode SNESGetConvergenceHistory(SNES snes,PetscReal *a[],PetscInt *its[],PetscInt *na) 3453 { 3454 PetscFunctionBegin; 3455 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3456 if (a) *a = snes->conv_hist; 3457 if (its) *its = snes->conv_hist_its; 3458 if (na) *na = snes->conv_hist_len; 3459 PetscFunctionReturn(0); 3460 } 3461 3462 #undef __FUNCT__ 3463 #define __FUNCT__ "SNESSetUpdate" 3464 /*@C 3465 SNESSetUpdate - Sets the general-purpose update function called 3466 at the beginning of every iteration of the nonlinear solve. Specifically 3467 it is called just before the Jacobian is "evaluated". 3468 3469 Logically Collective on SNES 3470 3471 Input Parameters: 3472 . snes - The nonlinear solver context 3473 . func - The function 3474 3475 Calling sequence of func: 3476 . func (SNES snes, PetscInt step); 3477 3478 . step - The current step of the iteration 3479 3480 Level: advanced 3481 3482 Note: This is NOT what one uses to update the ghost points before a function evaluation, that should be done at the beginning of your FormFunction() 3483 This is not used by most users. 3484 3485 .keywords: SNES, update 3486 3487 .seealso SNESDefaultUpdate(), SNESSetJacobian(), SNESSolve() 3488 @*/ 3489 PetscErrorCode SNESSetUpdate(SNES snes, PetscErrorCode (*func)(SNES, PetscInt)) 3490 { 3491 PetscFunctionBegin; 3492 PetscValidHeaderSpecific(snes, SNES_CLASSID,1); 3493 snes->ops->update = func; 3494 PetscFunctionReturn(0); 3495 } 3496 3497 #undef __FUNCT__ 3498 #define __FUNCT__ "SNESDefaultUpdate" 3499 /*@ 3500 SNESDefaultUpdate - The default update function which does nothing. 3501 3502 Not collective 3503 3504 Input Parameters: 3505 . snes - The nonlinear solver context 3506 . step - The current step of the iteration 3507 3508 Level: intermediate 3509 3510 .keywords: SNES, update 3511 .seealso SNESSetUpdate(), SNESDefaultRhsBC(), SNESDefaultShortolutionBC() 3512 @*/ 3513 PetscErrorCode SNESDefaultUpdate(SNES snes, PetscInt step) 3514 { 3515 PetscFunctionBegin; 3516 PetscFunctionReturn(0); 3517 } 3518 3519 #undef __FUNCT__ 3520 #define __FUNCT__ "SNESScaleStep_Private" 3521 /* 3522 SNESScaleStep_Private - Scales a step so that its length is less than the 3523 positive parameter delta. 3524 3525 Input Parameters: 3526 + snes - the SNES context 3527 . y - approximate solution of linear system 3528 . fnorm - 2-norm of current function 3529 - delta - trust region size 3530 3531 Output Parameters: 3532 + gpnorm - predicted function norm at the new point, assuming local 3533 linearization. The value is zero if the step lies within the trust 3534 region, and exceeds zero otherwise. 3535 - ynorm - 2-norm of the step 3536 3537 Note: 3538 For non-trust region methods such as SNESLS, the parameter delta 3539 is set to be the maximum allowable step size. 3540 3541 .keywords: SNES, nonlinear, scale, step 3542 */ 3543 PetscErrorCode SNESScaleStep_Private(SNES snes,Vec y,PetscReal *fnorm,PetscReal *delta,PetscReal *gpnorm,PetscReal *ynorm) 3544 { 3545 PetscReal nrm; 3546 PetscScalar cnorm; 3547 PetscErrorCode ierr; 3548 3549 PetscFunctionBegin; 3550 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3551 PetscValidHeaderSpecific(y,VEC_CLASSID,2); 3552 PetscCheckSameComm(snes,1,y,2); 3553 3554 ierr = VecNorm(y,NORM_2,&nrm);CHKERRQ(ierr); 3555 if (nrm > *delta) { 3556 nrm = *delta/nrm; 3557 *gpnorm = (1.0 - nrm)*(*fnorm); 3558 cnorm = nrm; 3559 ierr = VecScale(y,cnorm);CHKERRQ(ierr); 3560 *ynorm = *delta; 3561 } else { 3562 *gpnorm = 0.0; 3563 *ynorm = nrm; 3564 } 3565 PetscFunctionReturn(0); 3566 } 3567 3568 #undef __FUNCT__ 3569 #define __FUNCT__ "SNESSolve" 3570 /*@C 3571 SNESSolve - Solves a nonlinear system F(x) = b. 3572 Call SNESSolve() after calling SNESCreate() and optional routines of the form SNESSetXXX(). 3573 3574 Collective on SNES 3575 3576 Input Parameters: 3577 + snes - the SNES context 3578 . b - the constant part of the equation F(x) = b, or PETSC_NULL to use zero. 3579 - x - the solution vector. 3580 3581 Notes: 3582 The user should initialize the vector,x, with the initial guess 3583 for the nonlinear solve prior to calling SNESSolve. In particular, 3584 to employ an initial guess of zero, the user should explicitly set 3585 this vector to zero by calling VecSet(). 3586 3587 Level: beginner 3588 3589 .keywords: SNES, nonlinear, solve 3590 3591 .seealso: SNESCreate(), SNESDestroy(), SNESSetFunction(), SNESSetJacobian(), SNESSetGridSequence(), SNESGetSolution() 3592 @*/ 3593 PetscErrorCode SNESSolve(SNES snes,Vec b,Vec x) 3594 { 3595 PetscErrorCode ierr; 3596 PetscBool flg; 3597 char filename[PETSC_MAX_PATH_LEN]; 3598 PetscViewer viewer; 3599 PetscInt grid; 3600 Vec xcreated = PETSC_NULL; 3601 DM dm; 3602 3603 PetscFunctionBegin; 3604 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3605 if (x) PetscValidHeaderSpecific(x,VEC_CLASSID,3); 3606 if (x) PetscCheckSameComm(snes,1,x,3); 3607 if (b) PetscValidHeaderSpecific(b,VEC_CLASSID,2); 3608 if (b) PetscCheckSameComm(snes,1,b,2); 3609 3610 if (!x) { 3611 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 3612 ierr = DMCreateGlobalVector(dm,&xcreated);CHKERRQ(ierr); 3613 x = xcreated; 3614 } 3615 3616 for (grid=0; grid<snes->gridsequence; grid++) {ierr = PetscViewerASCIIPushTab(PETSC_VIEWER_STDOUT_(((PetscObject)snes)->comm));CHKERRQ(ierr);} 3617 for (grid=0; grid<snes->gridsequence+1; grid++) { 3618 3619 /* set solution vector */ 3620 if (!grid) {ierr = PetscObjectReference((PetscObject)x);CHKERRQ(ierr);} 3621 ierr = VecDestroy(&snes->vec_sol);CHKERRQ(ierr); 3622 snes->vec_sol = x; 3623 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 3624 3625 /* set affine vector if provided */ 3626 if (b) { ierr = PetscObjectReference((PetscObject)b);CHKERRQ(ierr); } 3627 ierr = VecDestroy(&snes->vec_rhs);CHKERRQ(ierr); 3628 snes->vec_rhs = b; 3629 3630 ierr = SNESSetUp(snes);CHKERRQ(ierr); 3631 3632 if (!grid) { 3633 if (snes->ops->computeinitialguess) { 3634 ierr = (*snes->ops->computeinitialguess)(snes,snes->vec_sol,snes->initialguessP);CHKERRQ(ierr); 3635 } else if (snes->dm) { 3636 PetscBool ig; 3637 ierr = DMHasInitialGuess(snes->dm,&ig);CHKERRQ(ierr); 3638 if (ig) { 3639 ierr = DMComputeInitialGuess(snes->dm,snes->vec_sol);CHKERRQ(ierr); 3640 } 3641 } 3642 } 3643 3644 if (snes->conv_hist_reset) snes->conv_hist_len = 0; 3645 snes->nfuncs = 0; snes->linear_its = 0; snes->numFailures = 0; 3646 3647 ierr = PetscLogEventBegin(SNES_Solve,snes,0,0,0);CHKERRQ(ierr); 3648 ierr = (*snes->ops->solve)(snes);CHKERRQ(ierr); 3649 ierr = PetscLogEventEnd(SNES_Solve,snes,0,0,0);CHKERRQ(ierr); 3650 if (snes->domainerror){ 3651 snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; 3652 snes->domainerror = PETSC_FALSE; 3653 } 3654 if (!snes->reason) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Internal error, solver returned without setting converged reason"); 3655 3656 flg = PETSC_FALSE; 3657 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_test_local_min",&flg,PETSC_NULL);CHKERRQ(ierr); 3658 if (flg && !PetscPreLoadingOn) { ierr = SNESTestLocalMin(snes);CHKERRQ(ierr); } 3659 if (snes->printreason) { 3660 ierr = PetscViewerASCIIAddTab(PETSC_VIEWER_STDOUT_(((PetscObject)snes)->comm),((PetscObject)snes)->tablevel);CHKERRQ(ierr); 3661 if (snes->reason > 0) { 3662 ierr = PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_(((PetscObject)snes)->comm),"Nonlinear solve converged due to %s iterations %D\n",SNESConvergedReasons[snes->reason],snes->iter);CHKERRQ(ierr); 3663 } else { 3664 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); 3665 } 3666 ierr = PetscViewerASCIISubtractTab(PETSC_VIEWER_STDOUT_(((PetscObject)snes)->comm),((PetscObject)snes)->tablevel);CHKERRQ(ierr); 3667 } 3668 3669 if (snes->errorifnotconverged && snes->reason < 0) SETERRQ(((PetscObject)snes)->comm,PETSC_ERR_NOT_CONVERGED,"SNESSolve has not converged"); 3670 if (grid < snes->gridsequence) { 3671 DM fine; 3672 Vec xnew; 3673 Mat interp; 3674 3675 ierr = DMRefine(snes->dm,((PetscObject)snes)->comm,&fine);CHKERRQ(ierr); 3676 if (!fine) SETERRQ(((PetscObject)snes)->comm,PETSC_ERR_ARG_INCOMP,"DMRefine() did not perform any refinement, cannot continue grid sequencing"); 3677 ierr = DMCreateInterpolation(snes->dm,fine,&interp,PETSC_NULL);CHKERRQ(ierr); 3678 ierr = DMCreateGlobalVector(fine,&xnew);CHKERRQ(ierr); 3679 ierr = MatInterpolate(interp,x,xnew);CHKERRQ(ierr); 3680 ierr = DMInterpolate(snes->dm,interp,fine);CHKERRQ(ierr); 3681 ierr = MatDestroy(&interp);CHKERRQ(ierr); 3682 x = xnew; 3683 3684 ierr = SNESReset(snes);CHKERRQ(ierr); 3685 ierr = SNESSetDM(snes,fine);CHKERRQ(ierr); 3686 ierr = DMDestroy(&fine);CHKERRQ(ierr); 3687 ierr = PetscViewerASCIIPopTab(PETSC_VIEWER_STDOUT_(((PetscObject)snes)->comm));CHKERRQ(ierr); 3688 } 3689 } 3690 /* monitoring and viewing */ 3691 flg = PETSC_FALSE; 3692 ierr = PetscOptionsGetString(((PetscObject)snes)->prefix,"-snes_view",filename,PETSC_MAX_PATH_LEN,&flg);CHKERRQ(ierr); 3693 if (flg && !PetscPreLoadingOn) { 3694 ierr = PetscViewerASCIIOpen(((PetscObject)snes)->comm,filename,&viewer);CHKERRQ(ierr); 3695 ierr = SNESView(snes,viewer);CHKERRQ(ierr); 3696 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 3697 } 3698 3699 flg = PETSC_FALSE; 3700 ierr = PetscOptionsGetString(((PetscObject)snes)->prefix,"-snes_view_solution_vtk",filename,PETSC_MAX_PATH_LEN,&flg);CHKERRQ(ierr); 3701 if (flg) { 3702 PetscViewer viewer; 3703 ierr = PetscViewerCreate(((PetscObject)snes)->comm,&viewer);CHKERRQ(ierr); 3704 ierr = PetscViewerSetType(viewer,PETSCVIEWERVTK);CHKERRQ(ierr); 3705 ierr = PetscViewerFileSetName(viewer,filename);CHKERRQ(ierr); 3706 ierr = VecView(snes->vec_sol,viewer);CHKERRQ(ierr); 3707 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 3708 } 3709 3710 ierr = VecDestroy(&xcreated);CHKERRQ(ierr); 3711 PetscFunctionReturn(0); 3712 } 3713 3714 /* --------- Internal routines for SNES Package --------- */ 3715 3716 #undef __FUNCT__ 3717 #define __FUNCT__ "SNESSetType" 3718 /*@C 3719 SNESSetType - Sets the method for the nonlinear solver. 3720 3721 Collective on SNES 3722 3723 Input Parameters: 3724 + snes - the SNES context 3725 - type - a known method 3726 3727 Options Database Key: 3728 . -snes_type <type> - Sets the method; use -help for a list 3729 of available methods (for instance, ls or tr) 3730 3731 Notes: 3732 See "petsc/include/petscsnes.h" for available methods (for instance) 3733 + SNESLS - Newton's method with line search 3734 (systems of nonlinear equations) 3735 . SNESTR - Newton's method with trust region 3736 (systems of nonlinear equations) 3737 3738 Normally, it is best to use the SNESSetFromOptions() command and then 3739 set the SNES solver type from the options database rather than by using 3740 this routine. Using the options database provides the user with 3741 maximum flexibility in evaluating the many nonlinear solvers. 3742 The SNESSetType() routine is provided for those situations where it 3743 is necessary to set the nonlinear solver independently of the command 3744 line or options database. This might be the case, for example, when 3745 the choice of solver changes during the execution of the program, 3746 and the user's application is taking responsibility for choosing the 3747 appropriate method. 3748 3749 Level: intermediate 3750 3751 .keywords: SNES, set, type 3752 3753 .seealso: SNESType, SNESCreate() 3754 3755 @*/ 3756 PetscErrorCode SNESSetType(SNES snes,SNESType type) 3757 { 3758 PetscErrorCode ierr,(*r)(SNES); 3759 PetscBool match; 3760 3761 PetscFunctionBegin; 3762 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3763 PetscValidCharPointer(type,2); 3764 3765 ierr = PetscObjectTypeCompare((PetscObject)snes,type,&match);CHKERRQ(ierr); 3766 if (match) PetscFunctionReturn(0); 3767 3768 ierr = PetscFListFind(SNESList,((PetscObject)snes)->comm,type,PETSC_TRUE,(void (**)(void)) &r);CHKERRQ(ierr); 3769 if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unable to find requested SNES type %s",type); 3770 /* Destroy the previous private SNES context */ 3771 if (snes->ops->destroy) { 3772 ierr = (*(snes)->ops->destroy)(snes);CHKERRQ(ierr); 3773 snes->ops->destroy = PETSC_NULL; 3774 } 3775 /* Reinitialize function pointers in SNESOps structure */ 3776 snes->ops->setup = 0; 3777 snes->ops->solve = 0; 3778 snes->ops->view = 0; 3779 snes->ops->setfromoptions = 0; 3780 snes->ops->destroy = 0; 3781 /* Call the SNESCreate_XXX routine for this particular Nonlinear solver */ 3782 snes->setupcalled = PETSC_FALSE; 3783 ierr = PetscObjectChangeTypeName((PetscObject)snes,type);CHKERRQ(ierr); 3784 ierr = (*r)(snes);CHKERRQ(ierr); 3785 #if defined(PETSC_HAVE_AMS) 3786 if (PetscAMSPublishAll) { 3787 ierr = PetscObjectAMSPublish((PetscObject)snes);CHKERRQ(ierr); 3788 } 3789 #endif 3790 PetscFunctionReturn(0); 3791 } 3792 3793 3794 /* --------------------------------------------------------------------- */ 3795 #undef __FUNCT__ 3796 #define __FUNCT__ "SNESRegisterDestroy" 3797 /*@ 3798 SNESRegisterDestroy - Frees the list of nonlinear solvers that were 3799 registered by SNESRegisterDynamic(). 3800 3801 Not Collective 3802 3803 Level: advanced 3804 3805 .keywords: SNES, nonlinear, register, destroy 3806 3807 .seealso: SNESRegisterAll(), SNESRegisterAll() 3808 @*/ 3809 PetscErrorCode SNESRegisterDestroy(void) 3810 { 3811 PetscErrorCode ierr; 3812 3813 PetscFunctionBegin; 3814 ierr = PetscFListDestroy(&SNESList);CHKERRQ(ierr); 3815 SNESRegisterAllCalled = PETSC_FALSE; 3816 PetscFunctionReturn(0); 3817 } 3818 3819 #undef __FUNCT__ 3820 #define __FUNCT__ "SNESGetType" 3821 /*@C 3822 SNESGetType - Gets the SNES method type and name (as a string). 3823 3824 Not Collective 3825 3826 Input Parameter: 3827 . snes - nonlinear solver context 3828 3829 Output Parameter: 3830 . type - SNES method (a character string) 3831 3832 Level: intermediate 3833 3834 .keywords: SNES, nonlinear, get, type, name 3835 @*/ 3836 PetscErrorCode SNESGetType(SNES snes,SNESType *type) 3837 { 3838 PetscFunctionBegin; 3839 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3840 PetscValidPointer(type,2); 3841 *type = ((PetscObject)snes)->type_name; 3842 PetscFunctionReturn(0); 3843 } 3844 3845 #undef __FUNCT__ 3846 #define __FUNCT__ "SNESGetSolution" 3847 /*@ 3848 SNESGetSolution - Returns the vector where the approximate solution is 3849 stored. This is the fine grid solution when using SNESSetGridSequence(). 3850 3851 Not Collective, but Vec is parallel if SNES is parallel 3852 3853 Input Parameter: 3854 . snes - the SNES context 3855 3856 Output Parameter: 3857 . x - the solution 3858 3859 Level: intermediate 3860 3861 .keywords: SNES, nonlinear, get, solution 3862 3863 .seealso: SNESGetSolutionUpdate(), SNESGetFunction() 3864 @*/ 3865 PetscErrorCode SNESGetSolution(SNES snes,Vec *x) 3866 { 3867 PetscFunctionBegin; 3868 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3869 PetscValidPointer(x,2); 3870 *x = snes->vec_sol; 3871 PetscFunctionReturn(0); 3872 } 3873 3874 #undef __FUNCT__ 3875 #define __FUNCT__ "SNESGetSolutionUpdate" 3876 /*@ 3877 SNESGetSolutionUpdate - Returns the vector where the solution update is 3878 stored. 3879 3880 Not Collective, but Vec is parallel if SNES is parallel 3881 3882 Input Parameter: 3883 . snes - the SNES context 3884 3885 Output Parameter: 3886 . x - the solution update 3887 3888 Level: advanced 3889 3890 .keywords: SNES, nonlinear, get, solution, update 3891 3892 .seealso: SNESGetSolution(), SNESGetFunction() 3893 @*/ 3894 PetscErrorCode SNESGetSolutionUpdate(SNES snes,Vec *x) 3895 { 3896 PetscFunctionBegin; 3897 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3898 PetscValidPointer(x,2); 3899 *x = snes->vec_sol_update; 3900 PetscFunctionReturn(0); 3901 } 3902 3903 #undef __FUNCT__ 3904 #define __FUNCT__ "SNESGetFunction" 3905 /*@C 3906 SNESGetFunction - Returns the vector where the function is stored. 3907 3908 Not Collective, but Vec is parallel if SNES is parallel. Collective if Vec is requested, but has not been created yet. 3909 3910 Input Parameter: 3911 . snes - the SNES context 3912 3913 Output Parameter: 3914 + r - the function (or PETSC_NULL) 3915 . func - the function (or PETSC_NULL) 3916 - ctx - the function context (or PETSC_NULL) 3917 3918 Level: advanced 3919 3920 .keywords: SNES, nonlinear, get, function 3921 3922 .seealso: SNESSetFunction(), SNESGetSolution() 3923 @*/ 3924 PetscErrorCode SNESGetFunction(SNES snes,Vec *r,PetscErrorCode (**func)(SNES,Vec,Vec,void*),void **ctx) 3925 { 3926 PetscErrorCode ierr; 3927 DM dm; 3928 3929 PetscFunctionBegin; 3930 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3931 if (r) { 3932 if (!snes->vec_func) { 3933 if (snes->vec_rhs) { 3934 ierr = VecDuplicate(snes->vec_rhs,&snes->vec_func);CHKERRQ(ierr); 3935 } else if (snes->vec_sol) { 3936 ierr = VecDuplicate(snes->vec_sol,&snes->vec_func);CHKERRQ(ierr); 3937 } else if (snes->dm) { 3938 ierr = DMCreateGlobalVector(snes->dm,&snes->vec_func);CHKERRQ(ierr); 3939 } 3940 } 3941 *r = snes->vec_func; 3942 } 3943 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 3944 ierr = DMSNESGetFunction(dm,func,ctx);CHKERRQ(ierr); 3945 PetscFunctionReturn(0); 3946 } 3947 3948 /*@C 3949 SNESGetGS - Returns the GS function and context. 3950 3951 Input Parameter: 3952 . snes - the SNES context 3953 3954 Output Parameter: 3955 + gsfunc - the function (or PETSC_NULL) 3956 - ctx - the function context (or PETSC_NULL) 3957 3958 Level: advanced 3959 3960 .keywords: SNES, nonlinear, get, function 3961 3962 .seealso: SNESSetGS(), SNESGetFunction() 3963 @*/ 3964 3965 #undef __FUNCT__ 3966 #define __FUNCT__ "SNESGetGS" 3967 PetscErrorCode SNESGetGS (SNES snes, PetscErrorCode(**func)(SNES, Vec, Vec, void*), void ** ctx) 3968 { 3969 PetscErrorCode ierr; 3970 DM dm; 3971 3972 PetscFunctionBegin; 3973 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3974 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 3975 ierr = DMSNESGetGS(dm,func,ctx);CHKERRQ(ierr); 3976 PetscFunctionReturn(0); 3977 } 3978 3979 #undef __FUNCT__ 3980 #define __FUNCT__ "SNESSetOptionsPrefix" 3981 /*@C 3982 SNESSetOptionsPrefix - Sets the prefix used for searching for all 3983 SNES options in the database. 3984 3985 Logically Collective on SNES 3986 3987 Input Parameter: 3988 + snes - the SNES context 3989 - prefix - the prefix to prepend to all option names 3990 3991 Notes: 3992 A hyphen (-) must NOT be given at the beginning of the prefix name. 3993 The first character of all runtime options is AUTOMATICALLY the hyphen. 3994 3995 Level: advanced 3996 3997 .keywords: SNES, set, options, prefix, database 3998 3999 .seealso: SNESSetFromOptions() 4000 @*/ 4001 PetscErrorCode SNESSetOptionsPrefix(SNES snes,const char prefix[]) 4002 { 4003 PetscErrorCode ierr; 4004 4005 PetscFunctionBegin; 4006 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4007 ierr = PetscObjectSetOptionsPrefix((PetscObject)snes,prefix);CHKERRQ(ierr); 4008 if (!snes->ksp) {ierr = SNESGetKSP(snes,&snes->ksp);CHKERRQ(ierr);} 4009 if (snes->linesearch) { 4010 ierr = SNESGetSNESLineSearch(snes,&snes->linesearch);CHKERRQ(ierr); 4011 ierr = PetscObjectSetOptionsPrefix((PetscObject)snes->linesearch,prefix);CHKERRQ(ierr); 4012 } 4013 ierr = KSPSetOptionsPrefix(snes->ksp,prefix);CHKERRQ(ierr); 4014 PetscFunctionReturn(0); 4015 } 4016 4017 #undef __FUNCT__ 4018 #define __FUNCT__ "SNESAppendOptionsPrefix" 4019 /*@C 4020 SNESAppendOptionsPrefix - Appends to the prefix used for searching for all 4021 SNES options in the database. 4022 4023 Logically Collective on SNES 4024 4025 Input Parameters: 4026 + snes - the SNES context 4027 - prefix - the prefix to prepend to all option names 4028 4029 Notes: 4030 A hyphen (-) must NOT be given at the beginning of the prefix name. 4031 The first character of all runtime options is AUTOMATICALLY the hyphen. 4032 4033 Level: advanced 4034 4035 .keywords: SNES, append, options, prefix, database 4036 4037 .seealso: SNESGetOptionsPrefix() 4038 @*/ 4039 PetscErrorCode SNESAppendOptionsPrefix(SNES snes,const char prefix[]) 4040 { 4041 PetscErrorCode ierr; 4042 4043 PetscFunctionBegin; 4044 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4045 ierr = PetscObjectAppendOptionsPrefix((PetscObject)snes,prefix);CHKERRQ(ierr); 4046 if (!snes->ksp) {ierr = SNESGetKSP(snes,&snes->ksp);CHKERRQ(ierr);} 4047 if (snes->linesearch) { 4048 ierr = SNESGetSNESLineSearch(snes,&snes->linesearch);CHKERRQ(ierr); 4049 ierr = PetscObjectAppendOptionsPrefix((PetscObject)snes->linesearch,prefix);CHKERRQ(ierr); 4050 } 4051 ierr = KSPAppendOptionsPrefix(snes->ksp,prefix);CHKERRQ(ierr); 4052 PetscFunctionReturn(0); 4053 } 4054 4055 #undef __FUNCT__ 4056 #define __FUNCT__ "SNESGetOptionsPrefix" 4057 /*@C 4058 SNESGetOptionsPrefix - Sets the prefix used for searching for all 4059 SNES options in the database. 4060 4061 Not Collective 4062 4063 Input Parameter: 4064 . snes - the SNES context 4065 4066 Output Parameter: 4067 . prefix - pointer to the prefix string used 4068 4069 Notes: On the fortran side, the user should pass in a string 'prefix' of 4070 sufficient length to hold the prefix. 4071 4072 Level: advanced 4073 4074 .keywords: SNES, get, options, prefix, database 4075 4076 .seealso: SNESAppendOptionsPrefix() 4077 @*/ 4078 PetscErrorCode SNESGetOptionsPrefix(SNES snes,const char *prefix[]) 4079 { 4080 PetscErrorCode ierr; 4081 4082 PetscFunctionBegin; 4083 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4084 ierr = PetscObjectGetOptionsPrefix((PetscObject)snes,prefix);CHKERRQ(ierr); 4085 PetscFunctionReturn(0); 4086 } 4087 4088 4089 #undef __FUNCT__ 4090 #define __FUNCT__ "SNESRegister" 4091 /*@C 4092 SNESRegister - See SNESRegisterDynamic() 4093 4094 Level: advanced 4095 @*/ 4096 PetscErrorCode SNESRegister(const char sname[],const char path[],const char name[],PetscErrorCode (*function)(SNES)) 4097 { 4098 char fullname[PETSC_MAX_PATH_LEN]; 4099 PetscErrorCode ierr; 4100 4101 PetscFunctionBegin; 4102 ierr = PetscFListConcat(path,name,fullname);CHKERRQ(ierr); 4103 ierr = PetscFListAdd(&SNESList,sname,fullname,(void (*)(void))function);CHKERRQ(ierr); 4104 PetscFunctionReturn(0); 4105 } 4106 4107 #undef __FUNCT__ 4108 #define __FUNCT__ "SNESTestLocalMin" 4109 PetscErrorCode SNESTestLocalMin(SNES snes) 4110 { 4111 PetscErrorCode ierr; 4112 PetscInt N,i,j; 4113 Vec u,uh,fh; 4114 PetscScalar value; 4115 PetscReal norm; 4116 4117 PetscFunctionBegin; 4118 ierr = SNESGetSolution(snes,&u);CHKERRQ(ierr); 4119 ierr = VecDuplicate(u,&uh);CHKERRQ(ierr); 4120 ierr = VecDuplicate(u,&fh);CHKERRQ(ierr); 4121 4122 /* currently only works for sequential */ 4123 ierr = PetscPrintf(PETSC_COMM_WORLD,"Testing FormFunction() for local min\n"); 4124 ierr = VecGetSize(u,&N);CHKERRQ(ierr); 4125 for (i=0; i<N; i++) { 4126 ierr = VecCopy(u,uh);CHKERRQ(ierr); 4127 ierr = PetscPrintf(PETSC_COMM_WORLD,"i = %D\n",i);CHKERRQ(ierr); 4128 for (j=-10; j<11; j++) { 4129 value = PetscSign(j)*exp(PetscAbs(j)-10.0); 4130 ierr = VecSetValue(uh,i,value,ADD_VALUES);CHKERRQ(ierr); 4131 ierr = SNESComputeFunction(snes,uh,fh);CHKERRQ(ierr); 4132 ierr = VecNorm(fh,NORM_2,&norm);CHKERRQ(ierr); 4133 ierr = PetscPrintf(PETSC_COMM_WORLD," j norm %D %18.16e\n",j,norm);CHKERRQ(ierr); 4134 value = -value; 4135 ierr = VecSetValue(uh,i,value,ADD_VALUES);CHKERRQ(ierr); 4136 } 4137 } 4138 ierr = VecDestroy(&uh);CHKERRQ(ierr); 4139 ierr = VecDestroy(&fh);CHKERRQ(ierr); 4140 PetscFunctionReturn(0); 4141 } 4142 4143 #undef __FUNCT__ 4144 #define __FUNCT__ "SNESKSPSetUseEW" 4145 /*@ 4146 SNESKSPSetUseEW - Sets SNES use Eisenstat-Walker method for 4147 computing relative tolerance for linear solvers within an inexact 4148 Newton method. 4149 4150 Logically Collective on SNES 4151 4152 Input Parameters: 4153 + snes - SNES context 4154 - flag - PETSC_TRUE or PETSC_FALSE 4155 4156 Options Database: 4157 + -snes_ksp_ew - use Eisenstat-Walker method for determining linear system convergence 4158 . -snes_ksp_ew_version ver - version of Eisenstat-Walker method 4159 . -snes_ksp_ew_rtol0 <rtol0> - Sets rtol0 4160 . -snes_ksp_ew_rtolmax <rtolmax> - Sets rtolmax 4161 . -snes_ksp_ew_gamma <gamma> - Sets gamma 4162 . -snes_ksp_ew_alpha <alpha> - Sets alpha 4163 . -snes_ksp_ew_alpha2 <alpha2> - Sets alpha2 4164 - -snes_ksp_ew_threshold <threshold> - Sets threshold 4165 4166 Notes: 4167 Currently, the default is to use a constant relative tolerance for 4168 the inner linear solvers. Alternatively, one can use the 4169 Eisenstat-Walker method, where the relative convergence tolerance 4170 is reset at each Newton iteration according progress of the nonlinear 4171 solver. 4172 4173 Level: advanced 4174 4175 Reference: 4176 S. C. Eisenstat and H. F. Walker, "Choosing the forcing terms in an 4177 inexact Newton method", SISC 17 (1), pp.16-32, 1996. 4178 4179 .keywords: SNES, KSP, Eisenstat, Walker, convergence, test, inexact, Newton 4180 4181 .seealso: SNESKSPGetUseEW(), SNESKSPGetParametersEW(), SNESKSPSetParametersEW() 4182 @*/ 4183 PetscErrorCode SNESKSPSetUseEW(SNES snes,PetscBool flag) 4184 { 4185 PetscFunctionBegin; 4186 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4187 PetscValidLogicalCollectiveBool(snes,flag,2); 4188 snes->ksp_ewconv = flag; 4189 PetscFunctionReturn(0); 4190 } 4191 4192 #undef __FUNCT__ 4193 #define __FUNCT__ "SNESKSPGetUseEW" 4194 /*@ 4195 SNESKSPGetUseEW - Gets if SNES is using Eisenstat-Walker method 4196 for computing relative tolerance for linear solvers within an 4197 inexact Newton method. 4198 4199 Not Collective 4200 4201 Input Parameter: 4202 . snes - SNES context 4203 4204 Output Parameter: 4205 . flag - PETSC_TRUE or PETSC_FALSE 4206 4207 Notes: 4208 Currently, the default is to use a constant relative tolerance for 4209 the inner linear solvers. Alternatively, one can use the 4210 Eisenstat-Walker method, where the relative convergence tolerance 4211 is reset at each Newton iteration according progress of the nonlinear 4212 solver. 4213 4214 Level: advanced 4215 4216 Reference: 4217 S. C. Eisenstat and H. F. Walker, "Choosing the forcing terms in an 4218 inexact Newton method", SISC 17 (1), pp.16-32, 1996. 4219 4220 .keywords: SNES, KSP, Eisenstat, Walker, convergence, test, inexact, Newton 4221 4222 .seealso: SNESKSPSetUseEW(), SNESKSPGetParametersEW(), SNESKSPSetParametersEW() 4223 @*/ 4224 PetscErrorCode SNESKSPGetUseEW(SNES snes, PetscBool *flag) 4225 { 4226 PetscFunctionBegin; 4227 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4228 PetscValidPointer(flag,2); 4229 *flag = snes->ksp_ewconv; 4230 PetscFunctionReturn(0); 4231 } 4232 4233 #undef __FUNCT__ 4234 #define __FUNCT__ "SNESKSPSetParametersEW" 4235 /*@ 4236 SNESKSPSetParametersEW - Sets parameters for Eisenstat-Walker 4237 convergence criteria for the linear solvers within an inexact 4238 Newton method. 4239 4240 Logically Collective on SNES 4241 4242 Input Parameters: 4243 + snes - SNES context 4244 . version - version 1, 2 (default is 2) or 3 4245 . rtol_0 - initial relative tolerance (0 <= rtol_0 < 1) 4246 . rtol_max - maximum relative tolerance (0 <= rtol_max < 1) 4247 . gamma - multiplicative factor for version 2 rtol computation 4248 (0 <= gamma2 <= 1) 4249 . alpha - power for version 2 rtol computation (1 < alpha <= 2) 4250 . alpha2 - power for safeguard 4251 - threshold - threshold for imposing safeguard (0 < threshold < 1) 4252 4253 Note: 4254 Version 3 was contributed by Luis Chacon, June 2006. 4255 4256 Use PETSC_DEFAULT to retain the default for any of the parameters. 4257 4258 Level: advanced 4259 4260 Reference: 4261 S. C. Eisenstat and H. F. Walker, "Choosing the forcing terms in an 4262 inexact Newton method", Utah State University Math. Stat. Dept. Res. 4263 Report 6/94/75, June, 1994, to appear in SIAM J. Sci. Comput. 4264 4265 .keywords: SNES, KSP, Eisenstat, Walker, set, parameters 4266 4267 .seealso: SNESKSPSetUseEW(), SNESKSPGetUseEW(), SNESKSPGetParametersEW() 4268 @*/ 4269 PetscErrorCode SNESKSPSetParametersEW(SNES snes,PetscInt version,PetscReal rtol_0,PetscReal rtol_max, 4270 PetscReal gamma,PetscReal alpha,PetscReal alpha2,PetscReal threshold) 4271 { 4272 SNESKSPEW *kctx; 4273 PetscFunctionBegin; 4274 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4275 kctx = (SNESKSPEW*)snes->kspconvctx; 4276 if (!kctx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"No Eisenstat-Walker context existing"); 4277 PetscValidLogicalCollectiveInt(snes,version,2); 4278 PetscValidLogicalCollectiveReal(snes,rtol_0,3); 4279 PetscValidLogicalCollectiveReal(snes,rtol_max,4); 4280 PetscValidLogicalCollectiveReal(snes,gamma,5); 4281 PetscValidLogicalCollectiveReal(snes,alpha,6); 4282 PetscValidLogicalCollectiveReal(snes,alpha2,7); 4283 PetscValidLogicalCollectiveReal(snes,threshold,8); 4284 4285 if (version != PETSC_DEFAULT) kctx->version = version; 4286 if (rtol_0 != PETSC_DEFAULT) kctx->rtol_0 = rtol_0; 4287 if (rtol_max != PETSC_DEFAULT) kctx->rtol_max = rtol_max; 4288 if (gamma != PETSC_DEFAULT) kctx->gamma = gamma; 4289 if (alpha != PETSC_DEFAULT) kctx->alpha = alpha; 4290 if (alpha2 != PETSC_DEFAULT) kctx->alpha2 = alpha2; 4291 if (threshold != PETSC_DEFAULT) kctx->threshold = threshold; 4292 4293 if (kctx->version < 1 || kctx->version > 3) { 4294 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Only versions 1, 2 and 3 are supported: %D",kctx->version); 4295 } 4296 if (kctx->rtol_0 < 0.0 || kctx->rtol_0 >= 1.0) { 4297 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"0.0 <= rtol_0 < 1.0: %G",kctx->rtol_0); 4298 } 4299 if (kctx->rtol_max < 0.0 || kctx->rtol_max >= 1.0) { 4300 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"0.0 <= rtol_max (%G) < 1.0\n",kctx->rtol_max); 4301 } 4302 if (kctx->gamma < 0.0 || kctx->gamma > 1.0) { 4303 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"0.0 <= gamma (%G) <= 1.0\n",kctx->gamma); 4304 } 4305 if (kctx->alpha <= 1.0 || kctx->alpha > 2.0) { 4306 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"1.0 < alpha (%G) <= 2.0\n",kctx->alpha); 4307 } 4308 if (kctx->threshold <= 0.0 || kctx->threshold >= 1.0) { 4309 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"0.0 < threshold (%G) < 1.0\n",kctx->threshold); 4310 } 4311 PetscFunctionReturn(0); 4312 } 4313 4314 #undef __FUNCT__ 4315 #define __FUNCT__ "SNESKSPGetParametersEW" 4316 /*@ 4317 SNESKSPGetParametersEW - Gets parameters for Eisenstat-Walker 4318 convergence criteria for the linear solvers within an inexact 4319 Newton method. 4320 4321 Not Collective 4322 4323 Input Parameters: 4324 snes - SNES context 4325 4326 Output Parameters: 4327 + version - version 1, 2 (default is 2) or 3 4328 . rtol_0 - initial relative tolerance (0 <= rtol_0 < 1) 4329 . rtol_max - maximum relative tolerance (0 <= rtol_max < 1) 4330 . gamma - multiplicative factor for version 2 rtol computation 4331 (0 <= gamma2 <= 1) 4332 . alpha - power for version 2 rtol computation (1 < alpha <= 2) 4333 . alpha2 - power for safeguard 4334 - threshold - threshold for imposing safeguard (0 < threshold < 1) 4335 4336 Level: advanced 4337 4338 .keywords: SNES, KSP, Eisenstat, Walker, get, parameters 4339 4340 .seealso: SNESKSPSetUseEW(), SNESKSPGetUseEW(), SNESKSPSetParametersEW() 4341 @*/ 4342 PetscErrorCode SNESKSPGetParametersEW(SNES snes,PetscInt *version,PetscReal *rtol_0,PetscReal *rtol_max, 4343 PetscReal *gamma,PetscReal *alpha,PetscReal *alpha2,PetscReal *threshold) 4344 { 4345 SNESKSPEW *kctx; 4346 PetscFunctionBegin; 4347 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4348 kctx = (SNESKSPEW*)snes->kspconvctx; 4349 if (!kctx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"No Eisenstat-Walker context existing"); 4350 if (version) *version = kctx->version; 4351 if (rtol_0) *rtol_0 = kctx->rtol_0; 4352 if (rtol_max) *rtol_max = kctx->rtol_max; 4353 if (gamma) *gamma = kctx->gamma; 4354 if (alpha) *alpha = kctx->alpha; 4355 if (alpha2) *alpha2 = kctx->alpha2; 4356 if (threshold) *threshold = kctx->threshold; 4357 PetscFunctionReturn(0); 4358 } 4359 4360 #undef __FUNCT__ 4361 #define __FUNCT__ "SNESKSPEW_PreSolve" 4362 static PetscErrorCode SNESKSPEW_PreSolve(SNES snes, KSP ksp, Vec b, Vec x) 4363 { 4364 PetscErrorCode ierr; 4365 SNESKSPEW *kctx = (SNESKSPEW*)snes->kspconvctx; 4366 PetscReal rtol=PETSC_DEFAULT,stol; 4367 4368 PetscFunctionBegin; 4369 if (!kctx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"No Eisenstat-Walker context exists"); 4370 if (!snes->iter) { /* first time in, so use the original user rtol */ 4371 rtol = kctx->rtol_0; 4372 } else { 4373 if (kctx->version == 1) { 4374 rtol = (snes->norm - kctx->lresid_last)/kctx->norm_last; 4375 if (rtol < 0.0) rtol = -rtol; 4376 stol = pow(kctx->rtol_last,kctx->alpha2); 4377 if (stol > kctx->threshold) rtol = PetscMax(rtol,stol); 4378 } else if (kctx->version == 2) { 4379 rtol = kctx->gamma * pow(snes->norm/kctx->norm_last,kctx->alpha); 4380 stol = kctx->gamma * pow(kctx->rtol_last,kctx->alpha); 4381 if (stol > kctx->threshold) rtol = PetscMax(rtol,stol); 4382 } else if (kctx->version == 3) {/* contributed by Luis Chacon, June 2006. */ 4383 rtol = kctx->gamma * pow(snes->norm/kctx->norm_last,kctx->alpha); 4384 /* safeguard: avoid sharp decrease of rtol */ 4385 stol = kctx->gamma*pow(kctx->rtol_last,kctx->alpha); 4386 stol = PetscMax(rtol,stol); 4387 rtol = PetscMin(kctx->rtol_0,stol); 4388 /* safeguard: avoid oversolving */ 4389 stol = kctx->gamma*(snes->ttol)/snes->norm; 4390 stol = PetscMax(rtol,stol); 4391 rtol = PetscMin(kctx->rtol_0,stol); 4392 } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Only versions 1, 2 or 3 are supported: %D",kctx->version); 4393 } 4394 /* safeguard: avoid rtol greater than one */ 4395 rtol = PetscMin(rtol,kctx->rtol_max); 4396 ierr = KSPSetTolerances(ksp,rtol,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);CHKERRQ(ierr); 4397 ierr = PetscInfo3(snes,"iter %D, Eisenstat-Walker (version %D) KSP rtol=%G\n",snes->iter,kctx->version,rtol);CHKERRQ(ierr); 4398 PetscFunctionReturn(0); 4399 } 4400 4401 #undef __FUNCT__ 4402 #define __FUNCT__ "SNESKSPEW_PostSolve" 4403 static PetscErrorCode SNESKSPEW_PostSolve(SNES snes, KSP ksp, Vec b, Vec x) 4404 { 4405 PetscErrorCode ierr; 4406 SNESKSPEW *kctx = (SNESKSPEW*)snes->kspconvctx; 4407 PCSide pcside; 4408 Vec lres; 4409 4410 PetscFunctionBegin; 4411 if (!kctx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"No Eisenstat-Walker context exists"); 4412 ierr = KSPGetTolerances(ksp,&kctx->rtol_last,0,0,0);CHKERRQ(ierr); 4413 ierr = SNESGetFunctionNorm(snes,&kctx->norm_last);CHKERRQ(ierr); 4414 if (kctx->version == 1) { 4415 ierr = KSPGetPCSide(ksp,&pcside);CHKERRQ(ierr); 4416 if (pcside == PC_RIGHT) { /* XXX Should we also test KSP_UNPRECONDITIONED_NORM ? */ 4417 /* KSP residual is true linear residual */ 4418 ierr = KSPGetResidualNorm(ksp,&kctx->lresid_last);CHKERRQ(ierr); 4419 } else { 4420 /* KSP residual is preconditioned residual */ 4421 /* compute true linear residual norm */ 4422 ierr = VecDuplicate(b,&lres);CHKERRQ(ierr); 4423 ierr = MatMult(snes->jacobian,x,lres);CHKERRQ(ierr); 4424 ierr = VecAYPX(lres,-1.0,b);CHKERRQ(ierr); 4425 ierr = VecNorm(lres,NORM_2,&kctx->lresid_last);CHKERRQ(ierr); 4426 ierr = VecDestroy(&lres);CHKERRQ(ierr); 4427 } 4428 } 4429 PetscFunctionReturn(0); 4430 } 4431 4432 #undef __FUNCT__ 4433 #define __FUNCT__ "SNES_KSPSolve" 4434 PetscErrorCode SNES_KSPSolve(SNES snes, KSP ksp, Vec b, Vec x) 4435 { 4436 PetscErrorCode ierr; 4437 4438 PetscFunctionBegin; 4439 if (snes->ksp_ewconv) { ierr = SNESKSPEW_PreSolve(snes,ksp,b,x);CHKERRQ(ierr); } 4440 ierr = KSPSolve(ksp,b,x);CHKERRQ(ierr); 4441 if (snes->ksp_ewconv) { ierr = SNESKSPEW_PostSolve(snes,ksp,b,x);CHKERRQ(ierr); } 4442 PetscFunctionReturn(0); 4443 } 4444 4445 #undef __FUNCT__ 4446 #define __FUNCT__ "SNESSetDM" 4447 /*@ 4448 SNESSetDM - Sets the DM that may be used by some preconditioners 4449 4450 Logically Collective on SNES 4451 4452 Input Parameters: 4453 + snes - the preconditioner context 4454 - dm - the dm 4455 4456 Level: intermediate 4457 4458 4459 .seealso: SNESGetDM(), KSPSetDM(), KSPGetDM() 4460 @*/ 4461 PetscErrorCode SNESSetDM(SNES snes,DM dm) 4462 { 4463 PetscErrorCode ierr; 4464 KSP ksp; 4465 SNESDM sdm; 4466 4467 PetscFunctionBegin; 4468 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4469 if (dm) {ierr = PetscObjectReference((PetscObject)dm);CHKERRQ(ierr);} 4470 if (snes->dm) { /* Move the SNESDM context over to the new DM unless the new DM already has one */ 4471 PetscContainer oldcontainer,container; 4472 ierr = PetscObjectQuery((PetscObject)snes->dm,"SNESDM",(PetscObject*)&oldcontainer);CHKERRQ(ierr); 4473 ierr = PetscObjectQuery((PetscObject)dm,"SNESDM",(PetscObject*)&container);CHKERRQ(ierr); 4474 if (oldcontainer && snes->dmAuto && !container) { 4475 ierr = DMSNESCopyContext(snes->dm,dm);CHKERRQ(ierr); 4476 ierr = DMSNESGetContext(snes->dm,&sdm);CHKERRQ(ierr); 4477 if (sdm->originaldm == snes->dm) { /* Grant write privileges to the replacement DM */ 4478 sdm->originaldm = dm; 4479 } 4480 } 4481 ierr = DMDestroy(&snes->dm);CHKERRQ(ierr); 4482 } 4483 snes->dm = dm; 4484 snes->dmAuto = PETSC_FALSE; 4485 ierr = SNESGetKSP(snes,&ksp);CHKERRQ(ierr); 4486 ierr = KSPSetDM(ksp,dm);CHKERRQ(ierr); 4487 ierr = KSPSetDMActive(ksp,PETSC_FALSE);CHKERRQ(ierr); 4488 if (snes->pc) { 4489 ierr = SNESSetDM(snes->pc, snes->dm);CHKERRQ(ierr); 4490 ierr = SNESSetPCSide(snes,snes->pcside);CHKERRQ(ierr); 4491 } 4492 PetscFunctionReturn(0); 4493 } 4494 4495 #undef __FUNCT__ 4496 #define __FUNCT__ "SNESGetDM" 4497 /*@ 4498 SNESGetDM - Gets the DM that may be used by some preconditioners 4499 4500 Not Collective but DM obtained is parallel on SNES 4501 4502 Input Parameter: 4503 . snes - the preconditioner context 4504 4505 Output Parameter: 4506 . dm - the dm 4507 4508 Level: intermediate 4509 4510 4511 .seealso: SNESSetDM(), KSPSetDM(), KSPGetDM() 4512 @*/ 4513 PetscErrorCode SNESGetDM(SNES snes,DM *dm) 4514 { 4515 PetscErrorCode ierr; 4516 4517 PetscFunctionBegin; 4518 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4519 if (!snes->dm) { 4520 ierr = DMShellCreate(((PetscObject)snes)->comm,&snes->dm);CHKERRQ(ierr); 4521 snes->dmAuto = PETSC_TRUE; 4522 } 4523 *dm = snes->dm; 4524 PetscFunctionReturn(0); 4525 } 4526 4527 #undef __FUNCT__ 4528 #define __FUNCT__ "SNESSetPC" 4529 /*@ 4530 SNESSetPC - Sets the nonlinear preconditioner to be used. 4531 4532 Collective on SNES 4533 4534 Input Parameters: 4535 + snes - iterative context obtained from SNESCreate() 4536 - pc - the preconditioner object 4537 4538 Notes: 4539 Use SNESGetPC() to retrieve the preconditioner context (for example, 4540 to configure it using the API). 4541 4542 Level: developer 4543 4544 .keywords: SNES, set, precondition 4545 .seealso: SNESGetPC() 4546 @*/ 4547 PetscErrorCode SNESSetPC(SNES snes, SNES pc) 4548 { 4549 PetscErrorCode ierr; 4550 4551 PetscFunctionBegin; 4552 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); 4553 PetscValidHeaderSpecific(pc, SNES_CLASSID, 2); 4554 PetscCheckSameComm(snes, 1, pc, 2); 4555 ierr = PetscObjectReference((PetscObject) pc);CHKERRQ(ierr); 4556 ierr = SNESDestroy(&snes->pc);CHKERRQ(ierr); 4557 snes->pc = pc; 4558 ierr = PetscLogObjectParent(snes, snes->pc);CHKERRQ(ierr); 4559 PetscFunctionReturn(0); 4560 } 4561 4562 #undef __FUNCT__ 4563 #define __FUNCT__ "SNESGetPC" 4564 /*@ 4565 SNESGetPC - Returns a pointer to the nonlinear preconditioning context set with SNESSetPC(). 4566 4567 Not Collective 4568 4569 Input Parameter: 4570 . snes - iterative context obtained from SNESCreate() 4571 4572 Output Parameter: 4573 . pc - preconditioner context 4574 4575 Level: developer 4576 4577 .keywords: SNES, get, preconditioner 4578 .seealso: SNESSetPC() 4579 @*/ 4580 PetscErrorCode SNESGetPC(SNES snes, SNES *pc) 4581 { 4582 PetscErrorCode ierr; 4583 const char *optionsprefix; 4584 4585 PetscFunctionBegin; 4586 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); 4587 PetscValidPointer(pc, 2); 4588 if (!snes->pc) { 4589 ierr = SNESCreate(((PetscObject) snes)->comm,&snes->pc);CHKERRQ(ierr); 4590 ierr = PetscObjectIncrementTabLevel((PetscObject)snes->pc,(PetscObject)snes,1);CHKERRQ(ierr); 4591 ierr = PetscLogObjectParent(snes,snes->pc);CHKERRQ(ierr); 4592 ierr = SNESGetOptionsPrefix(snes,&optionsprefix);CHKERRQ(ierr); 4593 ierr = SNESSetOptionsPrefix(snes->pc,optionsprefix);CHKERRQ(ierr); 4594 ierr = SNESAppendOptionsPrefix(snes->pc,"npc_");CHKERRQ(ierr); 4595 } 4596 *pc = snes->pc; 4597 PetscFunctionReturn(0); 4598 } 4599 4600 4601 #undef __FUNCT__ 4602 #define __FUNCT__ "SNESSetPCSide" 4603 /*@ 4604 SNESSetPCSide - Sets the preconditioning side. 4605 4606 Logically Collective on SNES 4607 4608 Input Parameter: 4609 . snes - iterative context obtained from SNESCreate() 4610 4611 Output Parameter: 4612 . side - the preconditioning side, where side is one of 4613 .vb 4614 PC_LEFT - left preconditioning (default) 4615 PC_RIGHT - right preconditioning 4616 .ve 4617 4618 Options Database Keys: 4619 . -snes_pc_side <right,left> 4620 4621 Level: intermediate 4622 4623 .keywords: SNES, set, right, left, side, preconditioner, flag 4624 4625 .seealso: SNESGetPCSide(), KSPSetPCSide() 4626 @*/ 4627 PetscErrorCode SNESSetPCSide(SNES snes,PCSide side) 4628 { 4629 PetscFunctionBegin; 4630 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4631 PetscValidLogicalCollectiveEnum(snes,side,2); 4632 snes->pcside = side; 4633 PetscFunctionReturn(0); 4634 } 4635 4636 #undef __FUNCT__ 4637 #define __FUNCT__ "SNESGetPCSide" 4638 /*@ 4639 SNESGetPCSide - Gets the preconditioning side. 4640 4641 Not Collective 4642 4643 Input Parameter: 4644 . snes - iterative context obtained from SNESCreate() 4645 4646 Output Parameter: 4647 . side - the preconditioning side, where side is one of 4648 .vb 4649 PC_LEFT - left preconditioning (default) 4650 PC_RIGHT - right preconditioning 4651 .ve 4652 4653 Level: intermediate 4654 4655 .keywords: SNES, get, right, left, side, preconditioner, flag 4656 4657 .seealso: SNESSetPCSide(), KSPGetPCSide() 4658 @*/ 4659 PetscErrorCode SNESGetPCSide(SNES snes,PCSide *side) 4660 { 4661 PetscFunctionBegin; 4662 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4663 PetscValidPointer(side,2); 4664 *side = snes->pcside; 4665 PetscFunctionReturn(0); 4666 } 4667 4668 #undef __FUNCT__ 4669 #define __FUNCT__ "SNESSetSNESLineSearch" 4670 /*@ 4671 SNESSetSNESLineSearch - Sets the linesearch on the SNES instance. 4672 4673 Collective on SNES 4674 4675 Input Parameters: 4676 + snes - iterative context obtained from SNESCreate() 4677 - linesearch - the linesearch object 4678 4679 Notes: 4680 Use SNESGetSNESLineSearch() to retrieve the preconditioner context (for example, 4681 to configure it using the API). 4682 4683 Level: developer 4684 4685 .keywords: SNES, set, linesearch 4686 .seealso: SNESGetSNESLineSearch() 4687 @*/ 4688 PetscErrorCode SNESSetSNESLineSearch(SNES snes, SNESLineSearch linesearch) 4689 { 4690 PetscErrorCode ierr; 4691 4692 PetscFunctionBegin; 4693 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); 4694 PetscValidHeaderSpecific(linesearch, SNESLINESEARCH_CLASSID, 2); 4695 PetscCheckSameComm(snes, 1, linesearch, 2); 4696 ierr = PetscObjectReference((PetscObject) linesearch);CHKERRQ(ierr); 4697 ierr = SNESLineSearchDestroy(&snes->linesearch);CHKERRQ(ierr); 4698 snes->linesearch = linesearch; 4699 ierr = PetscLogObjectParent(snes, snes->linesearch);CHKERRQ(ierr); 4700 PetscFunctionReturn(0); 4701 } 4702 4703 #undef __FUNCT__ 4704 #define __FUNCT__ "SNESGetSNESLineSearch" 4705 /*@C 4706 SNESGetSNESLineSearch - Returns a pointer to the line search context set with SNESSetLineSearch() 4707 or creates a default line search instance associated with the SNES and returns it. 4708 4709 Not Collective 4710 4711 Input Parameter: 4712 . snes - iterative context obtained from SNESCreate() 4713 4714 Output Parameter: 4715 . linesearch - linesearch context 4716 4717 Level: developer 4718 4719 .keywords: SNES, get, linesearch 4720 .seealso: SNESSetSNESLineSearch() 4721 @*/ 4722 PetscErrorCode SNESGetSNESLineSearch(SNES snes, SNESLineSearch *linesearch) 4723 { 4724 PetscErrorCode ierr; 4725 const char *optionsprefix; 4726 4727 PetscFunctionBegin; 4728 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); 4729 PetscValidPointer(linesearch, 2); 4730 if (!snes->linesearch) { 4731 ierr = SNESGetOptionsPrefix(snes, &optionsprefix);CHKERRQ(ierr); 4732 ierr = SNESLineSearchCreate(((PetscObject) snes)->comm, &snes->linesearch);CHKERRQ(ierr); 4733 ierr = SNESLineSearchSetSNES(snes->linesearch, snes);CHKERRQ(ierr); 4734 ierr = SNESLineSearchAppendOptionsPrefix(snes->linesearch, optionsprefix);CHKERRQ(ierr); 4735 ierr = PetscObjectIncrementTabLevel((PetscObject) snes->linesearch, (PetscObject) snes, 1);CHKERRQ(ierr); 4736 ierr = PetscLogObjectParent(snes, snes->linesearch);CHKERRQ(ierr); 4737 } 4738 *linesearch = snes->linesearch; 4739 PetscFunctionReturn(0); 4740 } 4741 4742 #if defined(PETSC_HAVE_MATLAB_ENGINE) 4743 #include <mex.h> 4744 4745 typedef struct {char *funcname; mxArray *ctx;} SNESMatlabContext; 4746 4747 #undef __FUNCT__ 4748 #define __FUNCT__ "SNESComputeFunction_Matlab" 4749 /* 4750 SNESComputeFunction_Matlab - Calls the function that has been set with 4751 SNESSetFunctionMatlab(). 4752 4753 Collective on SNES 4754 4755 Input Parameters: 4756 + snes - the SNES context 4757 - x - input vector 4758 4759 Output Parameter: 4760 . y - function vector, as set by SNESSetFunction() 4761 4762 Notes: 4763 SNESComputeFunction() is typically used within nonlinear solvers 4764 implementations, so most users would not generally call this routine 4765 themselves. 4766 4767 Level: developer 4768 4769 .keywords: SNES, nonlinear, compute, function 4770 4771 .seealso: SNESSetFunction(), SNESGetFunction() 4772 */ 4773 PetscErrorCode SNESComputeFunction_Matlab(SNES snes,Vec x,Vec y, void *ctx) 4774 { 4775 PetscErrorCode ierr; 4776 SNESMatlabContext *sctx = (SNESMatlabContext *)ctx; 4777 int nlhs = 1,nrhs = 5; 4778 mxArray *plhs[1],*prhs[5]; 4779 long long int lx = 0,ly = 0,ls = 0; 4780 4781 PetscFunctionBegin; 4782 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4783 PetscValidHeaderSpecific(x,VEC_CLASSID,2); 4784 PetscValidHeaderSpecific(y,VEC_CLASSID,3); 4785 PetscCheckSameComm(snes,1,x,2); 4786 PetscCheckSameComm(snes,1,y,3); 4787 4788 /* call Matlab function in ctx with arguments x and y */ 4789 4790 ierr = PetscMemcpy(&ls,&snes,sizeof(snes));CHKERRQ(ierr); 4791 ierr = PetscMemcpy(&lx,&x,sizeof(x));CHKERRQ(ierr); 4792 ierr = PetscMemcpy(&ly,&y,sizeof(x));CHKERRQ(ierr); 4793 prhs[0] = mxCreateDoubleScalar((double)ls); 4794 prhs[1] = mxCreateDoubleScalar((double)lx); 4795 prhs[2] = mxCreateDoubleScalar((double)ly); 4796 prhs[3] = mxCreateString(sctx->funcname); 4797 prhs[4] = sctx->ctx; 4798 ierr = mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscSNESComputeFunctionInternal");CHKERRQ(ierr); 4799 ierr = mxGetScalar(plhs[0]);CHKERRQ(ierr); 4800 mxDestroyArray(prhs[0]); 4801 mxDestroyArray(prhs[1]); 4802 mxDestroyArray(prhs[2]); 4803 mxDestroyArray(prhs[3]); 4804 mxDestroyArray(plhs[0]); 4805 PetscFunctionReturn(0); 4806 } 4807 4808 4809 #undef __FUNCT__ 4810 #define __FUNCT__ "SNESSetFunctionMatlab" 4811 /* 4812 SNESSetFunctionMatlab - Sets the function evaluation routine and function 4813 vector for use by the SNES routines in solving systems of nonlinear 4814 equations from MATLAB. Here the function is a string containing the name of a MATLAB function 4815 4816 Logically Collective on SNES 4817 4818 Input Parameters: 4819 + snes - the SNES context 4820 . r - vector to store function value 4821 - func - function evaluation routine 4822 4823 Calling sequence of func: 4824 $ func (SNES snes,Vec x,Vec f,void *ctx); 4825 4826 4827 Notes: 4828 The Newton-like methods typically solve linear systems of the form 4829 $ f'(x) x = -f(x), 4830 where f'(x) denotes the Jacobian matrix and f(x) is the function. 4831 4832 Level: beginner 4833 4834 .keywords: SNES, nonlinear, set, function 4835 4836 .seealso: SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESSetFunction() 4837 */ 4838 PetscErrorCode SNESSetFunctionMatlab(SNES snes,Vec r,const char *func,mxArray *ctx) 4839 { 4840 PetscErrorCode ierr; 4841 SNESMatlabContext *sctx; 4842 4843 PetscFunctionBegin; 4844 /* currently sctx is memory bleed */ 4845 ierr = PetscMalloc(sizeof(SNESMatlabContext),&sctx);CHKERRQ(ierr); 4846 ierr = PetscStrallocpy(func,&sctx->funcname);CHKERRQ(ierr); 4847 /* 4848 This should work, but it doesn't 4849 sctx->ctx = ctx; 4850 mexMakeArrayPersistent(sctx->ctx); 4851 */ 4852 sctx->ctx = mxDuplicateArray(ctx); 4853 ierr = SNESSetFunction(snes,r,SNESComputeFunction_Matlab,sctx);CHKERRQ(ierr); 4854 PetscFunctionReturn(0); 4855 } 4856 4857 #undef __FUNCT__ 4858 #define __FUNCT__ "SNESComputeJacobian_Matlab" 4859 /* 4860 SNESComputeJacobian_Matlab - Calls the function that has been set with 4861 SNESSetJacobianMatlab(). 4862 4863 Collective on SNES 4864 4865 Input Parameters: 4866 + snes - the SNES context 4867 . x - input vector 4868 . A, B - the matrices 4869 - ctx - user context 4870 4871 Output Parameter: 4872 . flag - structure of the matrix 4873 4874 Level: developer 4875 4876 .keywords: SNES, nonlinear, compute, function 4877 4878 .seealso: SNESSetFunction(), SNESGetFunction() 4879 @*/ 4880 PetscErrorCode SNESComputeJacobian_Matlab(SNES snes,Vec x,Mat *A,Mat *B,MatStructure *flag, void *ctx) 4881 { 4882 PetscErrorCode ierr; 4883 SNESMatlabContext *sctx = (SNESMatlabContext *)ctx; 4884 int nlhs = 2,nrhs = 6; 4885 mxArray *plhs[2],*prhs[6]; 4886 long long int lx = 0,lA = 0,ls = 0, lB = 0; 4887 4888 PetscFunctionBegin; 4889 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4890 PetscValidHeaderSpecific(x,VEC_CLASSID,2); 4891 4892 /* call Matlab function in ctx with arguments x and y */ 4893 4894 ierr = PetscMemcpy(&ls,&snes,sizeof(snes));CHKERRQ(ierr); 4895 ierr = PetscMemcpy(&lx,&x,sizeof(x));CHKERRQ(ierr); 4896 ierr = PetscMemcpy(&lA,A,sizeof(x));CHKERRQ(ierr); 4897 ierr = PetscMemcpy(&lB,B,sizeof(x));CHKERRQ(ierr); 4898 prhs[0] = mxCreateDoubleScalar((double)ls); 4899 prhs[1] = mxCreateDoubleScalar((double)lx); 4900 prhs[2] = mxCreateDoubleScalar((double)lA); 4901 prhs[3] = mxCreateDoubleScalar((double)lB); 4902 prhs[4] = mxCreateString(sctx->funcname); 4903 prhs[5] = sctx->ctx; 4904 ierr = mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscSNESComputeJacobianInternal");CHKERRQ(ierr); 4905 ierr = mxGetScalar(plhs[0]);CHKERRQ(ierr); 4906 *flag = (MatStructure) mxGetScalar(plhs[1]);CHKERRQ(ierr); 4907 mxDestroyArray(prhs[0]); 4908 mxDestroyArray(prhs[1]); 4909 mxDestroyArray(prhs[2]); 4910 mxDestroyArray(prhs[3]); 4911 mxDestroyArray(prhs[4]); 4912 mxDestroyArray(plhs[0]); 4913 mxDestroyArray(plhs[1]); 4914 PetscFunctionReturn(0); 4915 } 4916 4917 4918 #undef __FUNCT__ 4919 #define __FUNCT__ "SNESSetJacobianMatlab" 4920 /* 4921 SNESSetJacobianMatlab - Sets the Jacobian function evaluation routine and two empty Jacobian matrices 4922 vector for use by the SNES routines in solving systems of nonlinear 4923 equations from MATLAB. Here the function is a string containing the name of a MATLAB function 4924 4925 Logically Collective on SNES 4926 4927 Input Parameters: 4928 + snes - the SNES context 4929 . A,B - Jacobian matrices 4930 . func - function evaluation routine 4931 - ctx - user context 4932 4933 Calling sequence of func: 4934 $ flag = func (SNES snes,Vec x,Mat A,Mat B,void *ctx); 4935 4936 4937 Level: developer 4938 4939 .keywords: SNES, nonlinear, set, function 4940 4941 .seealso: SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESSetFunction() 4942 */ 4943 PetscErrorCode SNESSetJacobianMatlab(SNES snes,Mat A,Mat B,const char *func,mxArray *ctx) 4944 { 4945 PetscErrorCode ierr; 4946 SNESMatlabContext *sctx; 4947 4948 PetscFunctionBegin; 4949 /* currently sctx is memory bleed */ 4950 ierr = PetscMalloc(sizeof(SNESMatlabContext),&sctx);CHKERRQ(ierr); 4951 ierr = PetscStrallocpy(func,&sctx->funcname);CHKERRQ(ierr); 4952 /* 4953 This should work, but it doesn't 4954 sctx->ctx = ctx; 4955 mexMakeArrayPersistent(sctx->ctx); 4956 */ 4957 sctx->ctx = mxDuplicateArray(ctx); 4958 ierr = SNESSetJacobian(snes,A,B,SNESComputeJacobian_Matlab,sctx);CHKERRQ(ierr); 4959 PetscFunctionReturn(0); 4960 } 4961 4962 #undef __FUNCT__ 4963 #define __FUNCT__ "SNESMonitor_Matlab" 4964 /* 4965 SNESMonitor_Matlab - Calls the function that has been set with SNESMonitorSetMatlab(). 4966 4967 Collective on SNES 4968 4969 .seealso: SNESSetFunction(), SNESGetFunction() 4970 @*/ 4971 PetscErrorCode SNESMonitor_Matlab(SNES snes,PetscInt it, PetscReal fnorm, void *ctx) 4972 { 4973 PetscErrorCode ierr; 4974 SNESMatlabContext *sctx = (SNESMatlabContext *)ctx; 4975 int nlhs = 1,nrhs = 6; 4976 mxArray *plhs[1],*prhs[6]; 4977 long long int lx = 0,ls = 0; 4978 Vec x=snes->vec_sol; 4979 4980 PetscFunctionBegin; 4981 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4982 4983 ierr = PetscMemcpy(&ls,&snes,sizeof(snes));CHKERRQ(ierr); 4984 ierr = PetscMemcpy(&lx,&x,sizeof(x));CHKERRQ(ierr); 4985 prhs[0] = mxCreateDoubleScalar((double)ls); 4986 prhs[1] = mxCreateDoubleScalar((double)it); 4987 prhs[2] = mxCreateDoubleScalar((double)fnorm); 4988 prhs[3] = mxCreateDoubleScalar((double)lx); 4989 prhs[4] = mxCreateString(sctx->funcname); 4990 prhs[5] = sctx->ctx; 4991 ierr = mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscSNESMonitorInternal");CHKERRQ(ierr); 4992 ierr = mxGetScalar(plhs[0]);CHKERRQ(ierr); 4993 mxDestroyArray(prhs[0]); 4994 mxDestroyArray(prhs[1]); 4995 mxDestroyArray(prhs[2]); 4996 mxDestroyArray(prhs[3]); 4997 mxDestroyArray(prhs[4]); 4998 mxDestroyArray(plhs[0]); 4999 PetscFunctionReturn(0); 5000 } 5001 5002 5003 #undef __FUNCT__ 5004 #define __FUNCT__ "SNESMonitorSetMatlab" 5005 /* 5006 SNESMonitorSetMatlab - Sets the monitor function from MATLAB 5007 5008 Level: developer 5009 5010 .keywords: SNES, nonlinear, set, function 5011 5012 .seealso: SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESSetFunction() 5013 */ 5014 PetscErrorCode SNESMonitorSetMatlab(SNES snes,const char *func,mxArray *ctx) 5015 { 5016 PetscErrorCode ierr; 5017 SNESMatlabContext *sctx; 5018 5019 PetscFunctionBegin; 5020 /* currently sctx is memory bleed */ 5021 ierr = PetscMalloc(sizeof(SNESMatlabContext),&sctx);CHKERRQ(ierr); 5022 ierr = PetscStrallocpy(func,&sctx->funcname);CHKERRQ(ierr); 5023 /* 5024 This should work, but it doesn't 5025 sctx->ctx = ctx; 5026 mexMakeArrayPersistent(sctx->ctx); 5027 */ 5028 sctx->ctx = mxDuplicateArray(ctx); 5029 ierr = SNESMonitorSet(snes,SNESMonitor_Matlab,sctx,PETSC_NULL);CHKERRQ(ierr); 5030 PetscFunctionReturn(0); 5031 } 5032 5033 #endif 5034