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