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