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