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) {ierr = SNESGetSNESLineSearch(snes,&snes->linesearch);CHKERRQ(ierr);} 3898 ierr = KSPSetOptionsPrefix(snes->ksp,prefix);CHKERRQ(ierr); 3899 ierr = PetscObjectSetOptionsPrefix((PetscObject)snes->linesearch,prefix);CHKERRQ(ierr); 3900 PetscFunctionReturn(0); 3901 } 3902 3903 #undef __FUNCT__ 3904 #define __FUNCT__ "SNESAppendOptionsPrefix" 3905 /*@C 3906 SNESAppendOptionsPrefix - Appends to the prefix used for searching for all 3907 SNES options in the database. 3908 3909 Logically Collective on SNES 3910 3911 Input Parameters: 3912 + snes - the SNES context 3913 - prefix - the prefix to prepend to all option names 3914 3915 Notes: 3916 A hyphen (-) must NOT be given at the beginning of the prefix name. 3917 The first character of all runtime options is AUTOMATICALLY the hyphen. 3918 3919 Level: advanced 3920 3921 .keywords: SNES, append, options, prefix, database 3922 3923 .seealso: SNESGetOptionsPrefix() 3924 @*/ 3925 PetscErrorCode SNESAppendOptionsPrefix(SNES snes,const char prefix[]) 3926 { 3927 PetscErrorCode ierr; 3928 3929 PetscFunctionBegin; 3930 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3931 ierr = PetscObjectAppendOptionsPrefix((PetscObject)snes,prefix);CHKERRQ(ierr); 3932 if (!snes->ksp) {ierr = SNESGetKSP(snes,&snes->ksp);CHKERRQ(ierr);} 3933 if (!snes->linesearch) {ierr = SNESGetSNESLineSearch(snes,&snes->linesearch);CHKERRQ(ierr);} 3934 ierr = KSPAppendOptionsPrefix(snes->ksp,prefix);CHKERRQ(ierr); 3935 ierr = PetscObjectAppendOptionsPrefix((PetscObject)snes->linesearch,prefix);CHKERRQ(ierr); 3936 PetscFunctionReturn(0); 3937 } 3938 3939 #undef __FUNCT__ 3940 #define __FUNCT__ "SNESGetOptionsPrefix" 3941 /*@C 3942 SNESGetOptionsPrefix - Sets the prefix used for searching for all 3943 SNES options in the database. 3944 3945 Not Collective 3946 3947 Input Parameter: 3948 . snes - the SNES context 3949 3950 Output Parameter: 3951 . prefix - pointer to the prefix string used 3952 3953 Notes: On the fortran side, the user should pass in a string 'prefix' of 3954 sufficient length to hold the prefix. 3955 3956 Level: advanced 3957 3958 .keywords: SNES, get, options, prefix, database 3959 3960 .seealso: SNESAppendOptionsPrefix() 3961 @*/ 3962 PetscErrorCode SNESGetOptionsPrefix(SNES snes,const char *prefix[]) 3963 { 3964 PetscErrorCode ierr; 3965 3966 PetscFunctionBegin; 3967 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3968 ierr = PetscObjectGetOptionsPrefix((PetscObject)snes,prefix);CHKERRQ(ierr); 3969 PetscFunctionReturn(0); 3970 } 3971 3972 3973 #undef __FUNCT__ 3974 #define __FUNCT__ "SNESRegister" 3975 /*@C 3976 SNESRegister - See SNESRegisterDynamic() 3977 3978 Level: advanced 3979 @*/ 3980 PetscErrorCode SNESRegister(const char sname[],const char path[],const char name[],PetscErrorCode (*function)(SNES)) 3981 { 3982 char fullname[PETSC_MAX_PATH_LEN]; 3983 PetscErrorCode ierr; 3984 3985 PetscFunctionBegin; 3986 ierr = PetscFListConcat(path,name,fullname);CHKERRQ(ierr); 3987 ierr = PetscFListAdd(&SNESList,sname,fullname,(void (*)(void))function);CHKERRQ(ierr); 3988 PetscFunctionReturn(0); 3989 } 3990 3991 #undef __FUNCT__ 3992 #define __FUNCT__ "SNESTestLocalMin" 3993 PetscErrorCode SNESTestLocalMin(SNES snes) 3994 { 3995 PetscErrorCode ierr; 3996 PetscInt N,i,j; 3997 Vec u,uh,fh; 3998 PetscScalar value; 3999 PetscReal norm; 4000 4001 PetscFunctionBegin; 4002 ierr = SNESGetSolution(snes,&u);CHKERRQ(ierr); 4003 ierr = VecDuplicate(u,&uh);CHKERRQ(ierr); 4004 ierr = VecDuplicate(u,&fh);CHKERRQ(ierr); 4005 4006 /* currently only works for sequential */ 4007 ierr = PetscPrintf(PETSC_COMM_WORLD,"Testing FormFunction() for local min\n"); 4008 ierr = VecGetSize(u,&N);CHKERRQ(ierr); 4009 for (i=0; i<N; i++) { 4010 ierr = VecCopy(u,uh);CHKERRQ(ierr); 4011 ierr = PetscPrintf(PETSC_COMM_WORLD,"i = %D\n",i);CHKERRQ(ierr); 4012 for (j=-10; j<11; j++) { 4013 value = PetscSign(j)*exp(PetscAbs(j)-10.0); 4014 ierr = VecSetValue(uh,i,value,ADD_VALUES);CHKERRQ(ierr); 4015 ierr = SNESComputeFunction(snes,uh,fh);CHKERRQ(ierr); 4016 ierr = VecNorm(fh,NORM_2,&norm);CHKERRQ(ierr); 4017 ierr = PetscPrintf(PETSC_COMM_WORLD," j norm %D %18.16e\n",j,norm);CHKERRQ(ierr); 4018 value = -value; 4019 ierr = VecSetValue(uh,i,value,ADD_VALUES);CHKERRQ(ierr); 4020 } 4021 } 4022 ierr = VecDestroy(&uh);CHKERRQ(ierr); 4023 ierr = VecDestroy(&fh);CHKERRQ(ierr); 4024 PetscFunctionReturn(0); 4025 } 4026 4027 #undef __FUNCT__ 4028 #define __FUNCT__ "SNESKSPSetUseEW" 4029 /*@ 4030 SNESKSPSetUseEW - Sets SNES use Eisenstat-Walker method for 4031 computing relative tolerance for linear solvers within an inexact 4032 Newton method. 4033 4034 Logically Collective on SNES 4035 4036 Input Parameters: 4037 + snes - SNES context 4038 - flag - PETSC_TRUE or PETSC_FALSE 4039 4040 Options Database: 4041 + -snes_ksp_ew - use Eisenstat-Walker method for determining linear system convergence 4042 . -snes_ksp_ew_version ver - version of Eisenstat-Walker method 4043 . -snes_ksp_ew_rtol0 <rtol0> - Sets rtol0 4044 . -snes_ksp_ew_rtolmax <rtolmax> - Sets rtolmax 4045 . -snes_ksp_ew_gamma <gamma> - Sets gamma 4046 . -snes_ksp_ew_alpha <alpha> - Sets alpha 4047 . -snes_ksp_ew_alpha2 <alpha2> - Sets alpha2 4048 - -snes_ksp_ew_threshold <threshold> - Sets threshold 4049 4050 Notes: 4051 Currently, the default is to use a constant relative tolerance for 4052 the inner linear solvers. Alternatively, one can use the 4053 Eisenstat-Walker method, where the relative convergence tolerance 4054 is reset at each Newton iteration according progress of the nonlinear 4055 solver. 4056 4057 Level: advanced 4058 4059 Reference: 4060 S. C. Eisenstat and H. F. Walker, "Choosing the forcing terms in an 4061 inexact Newton method", SISC 17 (1), pp.16-32, 1996. 4062 4063 .keywords: SNES, KSP, Eisenstat, Walker, convergence, test, inexact, Newton 4064 4065 .seealso: SNESKSPGetUseEW(), SNESKSPGetParametersEW(), SNESKSPSetParametersEW() 4066 @*/ 4067 PetscErrorCode SNESKSPSetUseEW(SNES snes,PetscBool flag) 4068 { 4069 PetscFunctionBegin; 4070 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4071 PetscValidLogicalCollectiveBool(snes,flag,2); 4072 snes->ksp_ewconv = flag; 4073 PetscFunctionReturn(0); 4074 } 4075 4076 #undef __FUNCT__ 4077 #define __FUNCT__ "SNESKSPGetUseEW" 4078 /*@ 4079 SNESKSPGetUseEW - Gets if SNES is using Eisenstat-Walker method 4080 for computing relative tolerance for linear solvers within an 4081 inexact Newton method. 4082 4083 Not Collective 4084 4085 Input Parameter: 4086 . snes - SNES context 4087 4088 Output Parameter: 4089 . flag - PETSC_TRUE or PETSC_FALSE 4090 4091 Notes: 4092 Currently, the default is to use a constant relative tolerance for 4093 the inner linear solvers. Alternatively, one can use the 4094 Eisenstat-Walker method, where the relative convergence tolerance 4095 is reset at each Newton iteration according progress of the nonlinear 4096 solver. 4097 4098 Level: advanced 4099 4100 Reference: 4101 S. C. Eisenstat and H. F. Walker, "Choosing the forcing terms in an 4102 inexact Newton method", SISC 17 (1), pp.16-32, 1996. 4103 4104 .keywords: SNES, KSP, Eisenstat, Walker, convergence, test, inexact, Newton 4105 4106 .seealso: SNESKSPSetUseEW(), SNESKSPGetParametersEW(), SNESKSPSetParametersEW() 4107 @*/ 4108 PetscErrorCode SNESKSPGetUseEW(SNES snes, PetscBool *flag) 4109 { 4110 PetscFunctionBegin; 4111 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4112 PetscValidPointer(flag,2); 4113 *flag = snes->ksp_ewconv; 4114 PetscFunctionReturn(0); 4115 } 4116 4117 #undef __FUNCT__ 4118 #define __FUNCT__ "SNESKSPSetParametersEW" 4119 /*@ 4120 SNESKSPSetParametersEW - Sets parameters for Eisenstat-Walker 4121 convergence criteria for the linear solvers within an inexact 4122 Newton method. 4123 4124 Logically Collective on SNES 4125 4126 Input Parameters: 4127 + snes - SNES context 4128 . version - version 1, 2 (default is 2) or 3 4129 . rtol_0 - initial relative tolerance (0 <= rtol_0 < 1) 4130 . rtol_max - maximum relative tolerance (0 <= rtol_max < 1) 4131 . gamma - multiplicative factor for version 2 rtol computation 4132 (0 <= gamma2 <= 1) 4133 . alpha - power for version 2 rtol computation (1 < alpha <= 2) 4134 . alpha2 - power for safeguard 4135 - threshold - threshold for imposing safeguard (0 < threshold < 1) 4136 4137 Note: 4138 Version 3 was contributed by Luis Chacon, June 2006. 4139 4140 Use PETSC_DEFAULT to retain the default for any of the parameters. 4141 4142 Level: advanced 4143 4144 Reference: 4145 S. C. Eisenstat and H. F. Walker, "Choosing the forcing terms in an 4146 inexact Newton method", Utah State University Math. Stat. Dept. Res. 4147 Report 6/94/75, June, 1994, to appear in SIAM J. Sci. Comput. 4148 4149 .keywords: SNES, KSP, Eisenstat, Walker, set, parameters 4150 4151 .seealso: SNESKSPSetUseEW(), SNESKSPGetUseEW(), SNESKSPGetParametersEW() 4152 @*/ 4153 PetscErrorCode SNESKSPSetParametersEW(SNES snes,PetscInt version,PetscReal rtol_0,PetscReal rtol_max, 4154 PetscReal gamma,PetscReal alpha,PetscReal alpha2,PetscReal threshold) 4155 { 4156 SNESKSPEW *kctx; 4157 PetscFunctionBegin; 4158 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4159 kctx = (SNESKSPEW*)snes->kspconvctx; 4160 if (!kctx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"No Eisenstat-Walker context existing"); 4161 PetscValidLogicalCollectiveInt(snes,version,2); 4162 PetscValidLogicalCollectiveReal(snes,rtol_0,3); 4163 PetscValidLogicalCollectiveReal(snes,rtol_max,4); 4164 PetscValidLogicalCollectiveReal(snes,gamma,5); 4165 PetscValidLogicalCollectiveReal(snes,alpha,6); 4166 PetscValidLogicalCollectiveReal(snes,alpha2,7); 4167 PetscValidLogicalCollectiveReal(snes,threshold,8); 4168 4169 if (version != PETSC_DEFAULT) kctx->version = version; 4170 if (rtol_0 != PETSC_DEFAULT) kctx->rtol_0 = rtol_0; 4171 if (rtol_max != PETSC_DEFAULT) kctx->rtol_max = rtol_max; 4172 if (gamma != PETSC_DEFAULT) kctx->gamma = gamma; 4173 if (alpha != PETSC_DEFAULT) kctx->alpha = alpha; 4174 if (alpha2 != PETSC_DEFAULT) kctx->alpha2 = alpha2; 4175 if (threshold != PETSC_DEFAULT) kctx->threshold = threshold; 4176 4177 if (kctx->version < 1 || kctx->version > 3) { 4178 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Only versions 1, 2 and 3 are supported: %D",kctx->version); 4179 } 4180 if (kctx->rtol_0 < 0.0 || kctx->rtol_0 >= 1.0) { 4181 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"0.0 <= rtol_0 < 1.0: %G",kctx->rtol_0); 4182 } 4183 if (kctx->rtol_max < 0.0 || kctx->rtol_max >= 1.0) { 4184 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"0.0 <= rtol_max (%G) < 1.0\n",kctx->rtol_max); 4185 } 4186 if (kctx->gamma < 0.0 || kctx->gamma > 1.0) { 4187 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"0.0 <= gamma (%G) <= 1.0\n",kctx->gamma); 4188 } 4189 if (kctx->alpha <= 1.0 || kctx->alpha > 2.0) { 4190 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"1.0 < alpha (%G) <= 2.0\n",kctx->alpha); 4191 } 4192 if (kctx->threshold <= 0.0 || kctx->threshold >= 1.0) { 4193 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"0.0 < threshold (%G) < 1.0\n",kctx->threshold); 4194 } 4195 PetscFunctionReturn(0); 4196 } 4197 4198 #undef __FUNCT__ 4199 #define __FUNCT__ "SNESKSPGetParametersEW" 4200 /*@ 4201 SNESKSPGetParametersEW - Gets parameters for Eisenstat-Walker 4202 convergence criteria for the linear solvers within an inexact 4203 Newton method. 4204 4205 Not Collective 4206 4207 Input Parameters: 4208 snes - SNES context 4209 4210 Output Parameters: 4211 + version - version 1, 2 (default is 2) or 3 4212 . rtol_0 - initial relative tolerance (0 <= rtol_0 < 1) 4213 . rtol_max - maximum relative tolerance (0 <= rtol_max < 1) 4214 . gamma - multiplicative factor for version 2 rtol computation 4215 (0 <= gamma2 <= 1) 4216 . alpha - power for version 2 rtol computation (1 < alpha <= 2) 4217 . alpha2 - power for safeguard 4218 - threshold - threshold for imposing safeguard (0 < threshold < 1) 4219 4220 Level: advanced 4221 4222 .keywords: SNES, KSP, Eisenstat, Walker, get, parameters 4223 4224 .seealso: SNESKSPSetUseEW(), SNESKSPGetUseEW(), SNESKSPSetParametersEW() 4225 @*/ 4226 PetscErrorCode SNESKSPGetParametersEW(SNES snes,PetscInt *version,PetscReal *rtol_0,PetscReal *rtol_max, 4227 PetscReal *gamma,PetscReal *alpha,PetscReal *alpha2,PetscReal *threshold) 4228 { 4229 SNESKSPEW *kctx; 4230 PetscFunctionBegin; 4231 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4232 kctx = (SNESKSPEW*)snes->kspconvctx; 4233 if (!kctx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"No Eisenstat-Walker context existing"); 4234 if(version) *version = kctx->version; 4235 if(rtol_0) *rtol_0 = kctx->rtol_0; 4236 if(rtol_max) *rtol_max = kctx->rtol_max; 4237 if(gamma) *gamma = kctx->gamma; 4238 if(alpha) *alpha = kctx->alpha; 4239 if(alpha2) *alpha2 = kctx->alpha2; 4240 if(threshold) *threshold = kctx->threshold; 4241 PetscFunctionReturn(0); 4242 } 4243 4244 #undef __FUNCT__ 4245 #define __FUNCT__ "SNESKSPEW_PreSolve" 4246 static PetscErrorCode SNESKSPEW_PreSolve(SNES snes, KSP ksp, Vec b, Vec x) 4247 { 4248 PetscErrorCode ierr; 4249 SNESKSPEW *kctx = (SNESKSPEW*)snes->kspconvctx; 4250 PetscReal rtol=PETSC_DEFAULT,stol; 4251 4252 PetscFunctionBegin; 4253 if (!kctx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"No Eisenstat-Walker context exists"); 4254 if (!snes->iter) { /* first time in, so use the original user rtol */ 4255 rtol = kctx->rtol_0; 4256 } else { 4257 if (kctx->version == 1) { 4258 rtol = (snes->norm - kctx->lresid_last)/kctx->norm_last; 4259 if (rtol < 0.0) rtol = -rtol; 4260 stol = pow(kctx->rtol_last,kctx->alpha2); 4261 if (stol > kctx->threshold) rtol = PetscMax(rtol,stol); 4262 } else if (kctx->version == 2) { 4263 rtol = kctx->gamma * pow(snes->norm/kctx->norm_last,kctx->alpha); 4264 stol = kctx->gamma * pow(kctx->rtol_last,kctx->alpha); 4265 if (stol > kctx->threshold) rtol = PetscMax(rtol,stol); 4266 } else if (kctx->version == 3) {/* contributed by Luis Chacon, June 2006. */ 4267 rtol = kctx->gamma * pow(snes->norm/kctx->norm_last,kctx->alpha); 4268 /* safeguard: avoid sharp decrease of rtol */ 4269 stol = kctx->gamma*pow(kctx->rtol_last,kctx->alpha); 4270 stol = PetscMax(rtol,stol); 4271 rtol = PetscMin(kctx->rtol_0,stol); 4272 /* safeguard: avoid oversolving */ 4273 stol = kctx->gamma*(snes->ttol)/snes->norm; 4274 stol = PetscMax(rtol,stol); 4275 rtol = PetscMin(kctx->rtol_0,stol); 4276 } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Only versions 1, 2 or 3 are supported: %D",kctx->version); 4277 } 4278 /* safeguard: avoid rtol greater than one */ 4279 rtol = PetscMin(rtol,kctx->rtol_max); 4280 ierr = KSPSetTolerances(ksp,rtol,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);CHKERRQ(ierr); 4281 ierr = PetscInfo3(snes,"iter %D, Eisenstat-Walker (version %D) KSP rtol=%G\n",snes->iter,kctx->version,rtol);CHKERRQ(ierr); 4282 PetscFunctionReturn(0); 4283 } 4284 4285 #undef __FUNCT__ 4286 #define __FUNCT__ "SNESKSPEW_PostSolve" 4287 static PetscErrorCode SNESKSPEW_PostSolve(SNES snes, KSP ksp, Vec b, Vec x) 4288 { 4289 PetscErrorCode ierr; 4290 SNESKSPEW *kctx = (SNESKSPEW*)snes->kspconvctx; 4291 PCSide pcside; 4292 Vec lres; 4293 4294 PetscFunctionBegin; 4295 if (!kctx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"No Eisenstat-Walker context exists"); 4296 ierr = KSPGetTolerances(ksp,&kctx->rtol_last,0,0,0);CHKERRQ(ierr); 4297 ierr = SNESGetFunctionNorm(snes,&kctx->norm_last);CHKERRQ(ierr); 4298 if (kctx->version == 1) { 4299 ierr = KSPGetPCSide(ksp,&pcside);CHKERRQ(ierr); 4300 if (pcside == PC_RIGHT) { /* XXX Should we also test KSP_UNPRECONDITIONED_NORM ? */ 4301 /* KSP residual is true linear residual */ 4302 ierr = KSPGetResidualNorm(ksp,&kctx->lresid_last);CHKERRQ(ierr); 4303 } else { 4304 /* KSP residual is preconditioned residual */ 4305 /* compute true linear residual norm */ 4306 ierr = VecDuplicate(b,&lres);CHKERRQ(ierr); 4307 ierr = MatMult(snes->jacobian,x,lres);CHKERRQ(ierr); 4308 ierr = VecAYPX(lres,-1.0,b);CHKERRQ(ierr); 4309 ierr = VecNorm(lres,NORM_2,&kctx->lresid_last);CHKERRQ(ierr); 4310 ierr = VecDestroy(&lres);CHKERRQ(ierr); 4311 } 4312 } 4313 PetscFunctionReturn(0); 4314 } 4315 4316 #undef __FUNCT__ 4317 #define __FUNCT__ "SNES_KSPSolve" 4318 PetscErrorCode SNES_KSPSolve(SNES snes, KSP ksp, Vec b, Vec x) 4319 { 4320 PetscErrorCode ierr; 4321 4322 PetscFunctionBegin; 4323 if (snes->ksp_ewconv) { ierr = SNESKSPEW_PreSolve(snes,ksp,b,x);CHKERRQ(ierr); } 4324 ierr = KSPSolve(ksp,b,x);CHKERRQ(ierr); 4325 if (snes->ksp_ewconv) { ierr = SNESKSPEW_PostSolve(snes,ksp,b,x);CHKERRQ(ierr); } 4326 PetscFunctionReturn(0); 4327 } 4328 4329 #undef __FUNCT__ 4330 #define __FUNCT__ "SNESSetDM" 4331 /*@ 4332 SNESSetDM - Sets the DM that may be used by some preconditioners 4333 4334 Logically Collective on SNES 4335 4336 Input Parameters: 4337 + snes - the preconditioner context 4338 - dm - the dm 4339 4340 Level: intermediate 4341 4342 4343 .seealso: SNESGetDM(), KSPSetDM(), KSPGetDM() 4344 @*/ 4345 PetscErrorCode SNESSetDM(SNES snes,DM dm) 4346 { 4347 PetscErrorCode ierr; 4348 KSP ksp; 4349 SNESDM sdm; 4350 4351 PetscFunctionBegin; 4352 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4353 if (dm) {ierr = PetscObjectReference((PetscObject)dm);CHKERRQ(ierr);} 4354 if (snes->dm) { /* Move the SNESDM context over to the new DM unless the new DM already has one */ 4355 PetscContainer oldcontainer,container; 4356 ierr = PetscObjectQuery((PetscObject)snes->dm,"SNESDM",(PetscObject*)&oldcontainer);CHKERRQ(ierr); 4357 ierr = PetscObjectQuery((PetscObject)dm,"SNESDM",(PetscObject*)&container);CHKERRQ(ierr); 4358 if (oldcontainer && !container) { 4359 ierr = DMSNESCopyContext(snes->dm,dm);CHKERRQ(ierr); 4360 ierr = DMSNESGetContext(snes->dm,&sdm);CHKERRQ(ierr); 4361 if (sdm->originaldm == snes->dm) { /* Grant write privileges to the replacement DM */ 4362 sdm->originaldm = dm; 4363 } 4364 } 4365 ierr = DMDestroy(&snes->dm);CHKERRQ(ierr); 4366 } 4367 snes->dm = dm; 4368 ierr = SNESGetKSP(snes,&ksp);CHKERRQ(ierr); 4369 ierr = KSPSetDM(ksp,dm);CHKERRQ(ierr); 4370 ierr = KSPSetDMActive(ksp,PETSC_FALSE);CHKERRQ(ierr); 4371 if (snes->pc) { 4372 ierr = SNESSetDM(snes->pc, snes->dm);CHKERRQ(ierr); 4373 } 4374 PetscFunctionReturn(0); 4375 } 4376 4377 #undef __FUNCT__ 4378 #define __FUNCT__ "SNESGetDM" 4379 /*@ 4380 SNESGetDM - Gets the DM that may be used by some preconditioners 4381 4382 Not Collective but DM obtained is parallel on SNES 4383 4384 Input Parameter: 4385 . snes - the preconditioner context 4386 4387 Output Parameter: 4388 . dm - the dm 4389 4390 Level: intermediate 4391 4392 4393 .seealso: SNESSetDM(), KSPSetDM(), KSPGetDM() 4394 @*/ 4395 PetscErrorCode SNESGetDM(SNES snes,DM *dm) 4396 { 4397 PetscErrorCode ierr; 4398 4399 PetscFunctionBegin; 4400 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4401 if (!snes->dm) { 4402 ierr = DMShellCreate(((PetscObject)snes)->comm,&snes->dm);CHKERRQ(ierr); 4403 } 4404 *dm = snes->dm; 4405 PetscFunctionReturn(0); 4406 } 4407 4408 #undef __FUNCT__ 4409 #define __FUNCT__ "SNESSetPC" 4410 /*@ 4411 SNESSetPC - Sets the nonlinear preconditioner to be used. 4412 4413 Collective on SNES 4414 4415 Input Parameters: 4416 + snes - iterative context obtained from SNESCreate() 4417 - pc - the preconditioner object 4418 4419 Notes: 4420 Use SNESGetPC() to retrieve the preconditioner context (for example, 4421 to configure it using the API). 4422 4423 Level: developer 4424 4425 .keywords: SNES, set, precondition 4426 .seealso: SNESGetPC() 4427 @*/ 4428 PetscErrorCode SNESSetPC(SNES snes, SNES pc) 4429 { 4430 PetscErrorCode ierr; 4431 4432 PetscFunctionBegin; 4433 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); 4434 PetscValidHeaderSpecific(pc, SNES_CLASSID, 2); 4435 PetscCheckSameComm(snes, 1, pc, 2); 4436 ierr = PetscObjectReference((PetscObject) pc);CHKERRQ(ierr); 4437 ierr = SNESDestroy(&snes->pc);CHKERRQ(ierr); 4438 snes->pc = pc; 4439 ierr = PetscLogObjectParent(snes, snes->pc);CHKERRQ(ierr); 4440 PetscFunctionReturn(0); 4441 } 4442 4443 #undef __FUNCT__ 4444 #define __FUNCT__ "SNESGetPC" 4445 /*@ 4446 SNESGetPC - Returns a pointer to the nonlinear preconditioning context set with SNESSetPC(). 4447 4448 Not Collective 4449 4450 Input Parameter: 4451 . snes - iterative context obtained from SNESCreate() 4452 4453 Output Parameter: 4454 . pc - preconditioner context 4455 4456 Level: developer 4457 4458 .keywords: SNES, get, preconditioner 4459 .seealso: SNESSetPC() 4460 @*/ 4461 PetscErrorCode SNESGetPC(SNES snes, SNES *pc) 4462 { 4463 PetscErrorCode ierr; 4464 const char *optionsprefix; 4465 4466 PetscFunctionBegin; 4467 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); 4468 PetscValidPointer(pc, 2); 4469 if (!snes->pc) { 4470 ierr = SNESCreate(((PetscObject) snes)->comm,&snes->pc);CHKERRQ(ierr); 4471 ierr = PetscObjectIncrementTabLevel((PetscObject)snes->pc,(PetscObject)snes,1);CHKERRQ(ierr); 4472 ierr = PetscLogObjectParent(snes,snes->pc);CHKERRQ(ierr); 4473 ierr = SNESGetOptionsPrefix(snes,&optionsprefix);CHKERRQ(ierr); 4474 ierr = SNESSetOptionsPrefix(snes->pc,optionsprefix);CHKERRQ(ierr); 4475 ierr = SNESAppendOptionsPrefix(snes->pc,"npc_");CHKERRQ(ierr); 4476 } 4477 *pc = snes->pc; 4478 PetscFunctionReturn(0); 4479 } 4480 4481 #undef __FUNCT__ 4482 #define __FUNCT__ "SNESSetSNESLineSearch" 4483 /*@ 4484 SNESSetSNESLineSearch - Sets the linesearch on the SNES instance. 4485 4486 Collective on SNES 4487 4488 Input Parameters: 4489 + snes - iterative context obtained from SNESCreate() 4490 - linesearch - the linesearch object 4491 4492 Notes: 4493 Use SNESGetSNESLineSearch() to retrieve the preconditioner context (for example, 4494 to configure it using the API). 4495 4496 Level: developer 4497 4498 .keywords: SNES, set, linesearch 4499 .seealso: SNESGetSNESLineSearch() 4500 @*/ 4501 PetscErrorCode SNESSetSNESLineSearch(SNES snes, SNESLineSearch linesearch) 4502 { 4503 PetscErrorCode ierr; 4504 4505 PetscFunctionBegin; 4506 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); 4507 PetscValidHeaderSpecific(linesearch, SNESLINESEARCH_CLASSID, 2); 4508 PetscCheckSameComm(snes, 1, linesearch, 2); 4509 ierr = PetscObjectReference((PetscObject) linesearch);CHKERRQ(ierr); 4510 ierr = SNESLineSearchDestroy(&snes->linesearch);CHKERRQ(ierr); 4511 snes->linesearch = linesearch; 4512 ierr = PetscLogObjectParent(snes, snes->linesearch);CHKERRQ(ierr); 4513 PetscFunctionReturn(0); 4514 } 4515 4516 #undef __FUNCT__ 4517 #define __FUNCT__ "SNESGetSNESLineSearch" 4518 /*@C 4519 SNESGetSNESLineSearch - Returns a pointer to the line search context set with SNESSetLineSearch() 4520 or creates a default line search instance associated with the SNES and returns it. 4521 4522 Not Collective 4523 4524 Input Parameter: 4525 . snes - iterative context obtained from SNESCreate() 4526 4527 Output Parameter: 4528 . linesearch - linesearch context 4529 4530 Level: developer 4531 4532 .keywords: SNES, get, linesearch 4533 .seealso: SNESSetSNESLineSearch() 4534 @*/ 4535 PetscErrorCode SNESGetSNESLineSearch(SNES snes, SNESLineSearch *linesearch) 4536 { 4537 PetscErrorCode ierr; 4538 const char *optionsprefix; 4539 4540 PetscFunctionBegin; 4541 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); 4542 PetscValidPointer(linesearch, 2); 4543 if (!snes->linesearch) { 4544 ierr = SNESGetOptionsPrefix(snes, &optionsprefix);CHKERRQ(ierr); 4545 ierr = SNESLineSearchCreate(((PetscObject) snes)->comm, &snes->linesearch);CHKERRQ(ierr); 4546 ierr = SNESLineSearchSetSNES(snes->linesearch, snes);CHKERRQ(ierr); 4547 ierr = SNESLineSearchAppendOptionsPrefix(snes->linesearch, optionsprefix);CHKERRQ(ierr); 4548 ierr = PetscObjectIncrementTabLevel((PetscObject) snes->linesearch, (PetscObject) snes, 1);CHKERRQ(ierr); 4549 ierr = PetscLogObjectParent(snes, snes->linesearch);CHKERRQ(ierr); 4550 } 4551 *linesearch = snes->linesearch; 4552 PetscFunctionReturn(0); 4553 } 4554 4555 #if defined(PETSC_HAVE_MATLAB_ENGINE) 4556 #include <mex.h> 4557 4558 typedef struct {char *funcname; mxArray *ctx;} SNESMatlabContext; 4559 4560 #undef __FUNCT__ 4561 #define __FUNCT__ "SNESComputeFunction_Matlab" 4562 /* 4563 SNESComputeFunction_Matlab - Calls the function that has been set with 4564 SNESSetFunctionMatlab(). 4565 4566 Collective on SNES 4567 4568 Input Parameters: 4569 + snes - the SNES context 4570 - x - input vector 4571 4572 Output Parameter: 4573 . y - function vector, as set by SNESSetFunction() 4574 4575 Notes: 4576 SNESComputeFunction() is typically used within nonlinear solvers 4577 implementations, so most users would not generally call this routine 4578 themselves. 4579 4580 Level: developer 4581 4582 .keywords: SNES, nonlinear, compute, function 4583 4584 .seealso: SNESSetFunction(), SNESGetFunction() 4585 */ 4586 PetscErrorCode SNESComputeFunction_Matlab(SNES snes,Vec x,Vec y, void *ctx) 4587 { 4588 PetscErrorCode ierr; 4589 SNESMatlabContext *sctx = (SNESMatlabContext *)ctx; 4590 int nlhs = 1,nrhs = 5; 4591 mxArray *plhs[1],*prhs[5]; 4592 long long int lx = 0,ly = 0,ls = 0; 4593 4594 PetscFunctionBegin; 4595 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4596 PetscValidHeaderSpecific(x,VEC_CLASSID,2); 4597 PetscValidHeaderSpecific(y,VEC_CLASSID,3); 4598 PetscCheckSameComm(snes,1,x,2); 4599 PetscCheckSameComm(snes,1,y,3); 4600 4601 /* call Matlab function in ctx with arguments x and y */ 4602 4603 ierr = PetscMemcpy(&ls,&snes,sizeof(snes));CHKERRQ(ierr); 4604 ierr = PetscMemcpy(&lx,&x,sizeof(x));CHKERRQ(ierr); 4605 ierr = PetscMemcpy(&ly,&y,sizeof(x));CHKERRQ(ierr); 4606 prhs[0] = mxCreateDoubleScalar((double)ls); 4607 prhs[1] = mxCreateDoubleScalar((double)lx); 4608 prhs[2] = mxCreateDoubleScalar((double)ly); 4609 prhs[3] = mxCreateString(sctx->funcname); 4610 prhs[4] = sctx->ctx; 4611 ierr = mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscSNESComputeFunctionInternal");CHKERRQ(ierr); 4612 ierr = mxGetScalar(plhs[0]);CHKERRQ(ierr); 4613 mxDestroyArray(prhs[0]); 4614 mxDestroyArray(prhs[1]); 4615 mxDestroyArray(prhs[2]); 4616 mxDestroyArray(prhs[3]); 4617 mxDestroyArray(plhs[0]); 4618 PetscFunctionReturn(0); 4619 } 4620 4621 4622 #undef __FUNCT__ 4623 #define __FUNCT__ "SNESSetFunctionMatlab" 4624 /* 4625 SNESSetFunctionMatlab - Sets the function evaluation routine and function 4626 vector for use by the SNES routines in solving systems of nonlinear 4627 equations from MATLAB. Here the function is a string containing the name of a MATLAB function 4628 4629 Logically Collective on SNES 4630 4631 Input Parameters: 4632 + snes - the SNES context 4633 . r - vector to store function value 4634 - func - function evaluation routine 4635 4636 Calling sequence of func: 4637 $ func (SNES snes,Vec x,Vec f,void *ctx); 4638 4639 4640 Notes: 4641 The Newton-like methods typically solve linear systems of the form 4642 $ f'(x) x = -f(x), 4643 where f'(x) denotes the Jacobian matrix and f(x) is the function. 4644 4645 Level: beginner 4646 4647 .keywords: SNES, nonlinear, set, function 4648 4649 .seealso: SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESSetFunction() 4650 */ 4651 PetscErrorCode SNESSetFunctionMatlab(SNES snes,Vec r,const char *func,mxArray *ctx) 4652 { 4653 PetscErrorCode ierr; 4654 SNESMatlabContext *sctx; 4655 4656 PetscFunctionBegin; 4657 /* currently sctx is memory bleed */ 4658 ierr = PetscMalloc(sizeof(SNESMatlabContext),&sctx);CHKERRQ(ierr); 4659 ierr = PetscStrallocpy(func,&sctx->funcname);CHKERRQ(ierr); 4660 /* 4661 This should work, but it doesn't 4662 sctx->ctx = ctx; 4663 mexMakeArrayPersistent(sctx->ctx); 4664 */ 4665 sctx->ctx = mxDuplicateArray(ctx); 4666 ierr = SNESSetFunction(snes,r,SNESComputeFunction_Matlab,sctx);CHKERRQ(ierr); 4667 PetscFunctionReturn(0); 4668 } 4669 4670 #undef __FUNCT__ 4671 #define __FUNCT__ "SNESComputeJacobian_Matlab" 4672 /* 4673 SNESComputeJacobian_Matlab - Calls the function that has been set with 4674 SNESSetJacobianMatlab(). 4675 4676 Collective on SNES 4677 4678 Input Parameters: 4679 + snes - the SNES context 4680 . x - input vector 4681 . A, B - the matrices 4682 - ctx - user context 4683 4684 Output Parameter: 4685 . flag - structure of the matrix 4686 4687 Level: developer 4688 4689 .keywords: SNES, nonlinear, compute, function 4690 4691 .seealso: SNESSetFunction(), SNESGetFunction() 4692 @*/ 4693 PetscErrorCode SNESComputeJacobian_Matlab(SNES snes,Vec x,Mat *A,Mat *B,MatStructure *flag, void *ctx) 4694 { 4695 PetscErrorCode ierr; 4696 SNESMatlabContext *sctx = (SNESMatlabContext *)ctx; 4697 int nlhs = 2,nrhs = 6; 4698 mxArray *plhs[2],*prhs[6]; 4699 long long int lx = 0,lA = 0,ls = 0, lB = 0; 4700 4701 PetscFunctionBegin; 4702 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4703 PetscValidHeaderSpecific(x,VEC_CLASSID,2); 4704 4705 /* call Matlab function in ctx with arguments x and y */ 4706 4707 ierr = PetscMemcpy(&ls,&snes,sizeof(snes));CHKERRQ(ierr); 4708 ierr = PetscMemcpy(&lx,&x,sizeof(x));CHKERRQ(ierr); 4709 ierr = PetscMemcpy(&lA,A,sizeof(x));CHKERRQ(ierr); 4710 ierr = PetscMemcpy(&lB,B,sizeof(x));CHKERRQ(ierr); 4711 prhs[0] = mxCreateDoubleScalar((double)ls); 4712 prhs[1] = mxCreateDoubleScalar((double)lx); 4713 prhs[2] = mxCreateDoubleScalar((double)lA); 4714 prhs[3] = mxCreateDoubleScalar((double)lB); 4715 prhs[4] = mxCreateString(sctx->funcname); 4716 prhs[5] = sctx->ctx; 4717 ierr = mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscSNESComputeJacobianInternal");CHKERRQ(ierr); 4718 ierr = mxGetScalar(plhs[0]);CHKERRQ(ierr); 4719 *flag = (MatStructure) mxGetScalar(plhs[1]);CHKERRQ(ierr); 4720 mxDestroyArray(prhs[0]); 4721 mxDestroyArray(prhs[1]); 4722 mxDestroyArray(prhs[2]); 4723 mxDestroyArray(prhs[3]); 4724 mxDestroyArray(prhs[4]); 4725 mxDestroyArray(plhs[0]); 4726 mxDestroyArray(plhs[1]); 4727 PetscFunctionReturn(0); 4728 } 4729 4730 4731 #undef __FUNCT__ 4732 #define __FUNCT__ "SNESSetJacobianMatlab" 4733 /* 4734 SNESSetJacobianMatlab - Sets the Jacobian function evaluation routine and two empty Jacobian matrices 4735 vector for use by the SNES routines in solving systems of nonlinear 4736 equations from MATLAB. Here the function is a string containing the name of a MATLAB function 4737 4738 Logically Collective on SNES 4739 4740 Input Parameters: 4741 + snes - the SNES context 4742 . A,B - Jacobian matrices 4743 . func - function evaluation routine 4744 - ctx - user context 4745 4746 Calling sequence of func: 4747 $ flag = func (SNES snes,Vec x,Mat A,Mat B,void *ctx); 4748 4749 4750 Level: developer 4751 4752 .keywords: SNES, nonlinear, set, function 4753 4754 .seealso: SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESSetFunction() 4755 */ 4756 PetscErrorCode SNESSetJacobianMatlab(SNES snes,Mat A,Mat B,const char *func,mxArray *ctx) 4757 { 4758 PetscErrorCode ierr; 4759 SNESMatlabContext *sctx; 4760 4761 PetscFunctionBegin; 4762 /* currently sctx is memory bleed */ 4763 ierr = PetscMalloc(sizeof(SNESMatlabContext),&sctx);CHKERRQ(ierr); 4764 ierr = PetscStrallocpy(func,&sctx->funcname);CHKERRQ(ierr); 4765 /* 4766 This should work, but it doesn't 4767 sctx->ctx = ctx; 4768 mexMakeArrayPersistent(sctx->ctx); 4769 */ 4770 sctx->ctx = mxDuplicateArray(ctx); 4771 ierr = SNESSetJacobian(snes,A,B,SNESComputeJacobian_Matlab,sctx);CHKERRQ(ierr); 4772 PetscFunctionReturn(0); 4773 } 4774 4775 #undef __FUNCT__ 4776 #define __FUNCT__ "SNESMonitor_Matlab" 4777 /* 4778 SNESMonitor_Matlab - Calls the function that has been set with SNESMonitorSetMatlab(). 4779 4780 Collective on SNES 4781 4782 .seealso: SNESSetFunction(), SNESGetFunction() 4783 @*/ 4784 PetscErrorCode SNESMonitor_Matlab(SNES snes,PetscInt it, PetscReal fnorm, void *ctx) 4785 { 4786 PetscErrorCode ierr; 4787 SNESMatlabContext *sctx = (SNESMatlabContext *)ctx; 4788 int nlhs = 1,nrhs = 6; 4789 mxArray *plhs[1],*prhs[6]; 4790 long long int lx = 0,ls = 0; 4791 Vec x=snes->vec_sol; 4792 4793 PetscFunctionBegin; 4794 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4795 4796 ierr = PetscMemcpy(&ls,&snes,sizeof(snes));CHKERRQ(ierr); 4797 ierr = PetscMemcpy(&lx,&x,sizeof(x));CHKERRQ(ierr); 4798 prhs[0] = mxCreateDoubleScalar((double)ls); 4799 prhs[1] = mxCreateDoubleScalar((double)it); 4800 prhs[2] = mxCreateDoubleScalar((double)fnorm); 4801 prhs[3] = mxCreateDoubleScalar((double)lx); 4802 prhs[4] = mxCreateString(sctx->funcname); 4803 prhs[5] = sctx->ctx; 4804 ierr = mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscSNESMonitorInternal");CHKERRQ(ierr); 4805 ierr = mxGetScalar(plhs[0]);CHKERRQ(ierr); 4806 mxDestroyArray(prhs[0]); 4807 mxDestroyArray(prhs[1]); 4808 mxDestroyArray(prhs[2]); 4809 mxDestroyArray(prhs[3]); 4810 mxDestroyArray(prhs[4]); 4811 mxDestroyArray(plhs[0]); 4812 PetscFunctionReturn(0); 4813 } 4814 4815 4816 #undef __FUNCT__ 4817 #define __FUNCT__ "SNESMonitorSetMatlab" 4818 /* 4819 SNESMonitorSetMatlab - Sets the monitor function from MATLAB 4820 4821 Level: developer 4822 4823 .keywords: SNES, nonlinear, set, function 4824 4825 .seealso: SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESSetFunction() 4826 */ 4827 PetscErrorCode SNESMonitorSetMatlab(SNES snes,const char *func,mxArray *ctx) 4828 { 4829 PetscErrorCode ierr; 4830 SNESMatlabContext *sctx; 4831 4832 PetscFunctionBegin; 4833 /* currently sctx is memory bleed */ 4834 ierr = PetscMalloc(sizeof(SNESMatlabContext),&sctx);CHKERRQ(ierr); 4835 ierr = PetscStrallocpy(func,&sctx->funcname);CHKERRQ(ierr); 4836 /* 4837 This should work, but it doesn't 4838 sctx->ctx = ctx; 4839 mexMakeArrayPersistent(sctx->ctx); 4840 */ 4841 sctx->ctx = mxDuplicateArray(ctx); 4842 ierr = SNESMonitorSet(snes,SNESMonitor_Matlab,sctx,PETSC_NULL);CHKERRQ(ierr); 4843 PetscFunctionReturn(0); 4844 } 4845 4846 #endif 4847