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