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