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