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