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