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