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