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 } else if (!snes->vec_func && snes->dm) { 1370 ierr = DMCreateGlobalVector(snes->dm,&snes->vec_func);CHKERRQ(ierr); 1371 } 1372 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 1373 ierr = DMSNESSetFunction(dm,func,ctx);CHKERRQ(ierr); 1374 PetscFunctionReturn(0); 1375 } 1376 1377 1378 #undef __FUNCT__ 1379 #define __FUNCT__ "SNESSetGS" 1380 /*@C 1381 SNESSetGS - Sets the user nonlinear Gauss-Seidel routine for 1382 use with composed nonlinear solvers. 1383 1384 Input Parameters: 1385 + snes - the SNES context 1386 . gsfunc - function evaluation routine 1387 - ctx - [optional] user-defined context for private data for the 1388 smoother evaluation routine (may be PETSC_NULL) 1389 1390 Calling sequence of func: 1391 $ func (SNES snes,Vec x,Vec b,void *ctx); 1392 1393 + X - solution vector 1394 . B - RHS vector 1395 - ctx - optional user-defined Gauss-Seidel context 1396 1397 Notes: 1398 The GS routines are used by the composed nonlinear solver to generate 1399 a problem appropriate update to the solution, particularly FAS. 1400 1401 Level: intermediate 1402 1403 .keywords: SNES, nonlinear, set, Gauss-Seidel 1404 1405 .seealso: SNESGetFunction(), SNESComputeGS() 1406 @*/ 1407 PetscErrorCode SNESSetGS(SNES snes,PetscErrorCode (*gsfunc)(SNES,Vec,Vec,void*),void *ctx) 1408 { 1409 PetscErrorCode ierr; 1410 DM dm; 1411 1412 PetscFunctionBegin; 1413 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1414 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 1415 ierr = DMSNESSetGS(dm,gsfunc,ctx);CHKERRQ(ierr); 1416 PetscFunctionReturn(0); 1417 } 1418 1419 #undef __FUNCT__ 1420 #define __FUNCT__ "SNESSetGSSweeps" 1421 /*@ 1422 SNESSetGSSweeps - Sets the number of sweeps of GS to use. 1423 1424 Input Parameters: 1425 + snes - the SNES context 1426 - sweeps - the number of sweeps of GS to perform. 1427 1428 Level: intermediate 1429 1430 .keywords: SNES, nonlinear, set, Gauss-Siedel 1431 1432 .seealso: SNESSetGS(), SNESGetGS(), SNESSetPC(), SNESGetGSSweeps() 1433 @*/ 1434 1435 PetscErrorCode SNESSetGSSweeps(SNES snes, PetscInt sweeps) { 1436 PetscFunctionBegin; 1437 snes->gssweeps = sweeps; 1438 PetscFunctionReturn(0); 1439 } 1440 1441 1442 #undef __FUNCT__ 1443 #define __FUNCT__ "SNESGetGSSweeps" 1444 /*@ 1445 SNESGetGSSweeps - Gets the number of sweeps GS will use. 1446 1447 Input Parameters: 1448 . snes - the SNES context 1449 1450 Output Parameters: 1451 . sweeps - the number of sweeps of GS to perform. 1452 1453 Level: intermediate 1454 1455 .keywords: SNES, nonlinear, set, Gauss-Siedel 1456 1457 .seealso: SNESSetGS(), SNESGetGS(), SNESSetPC(), SNESSetGSSweeps() 1458 @*/ 1459 PetscErrorCode SNESGetGSSweeps(SNES snes, PetscInt * sweeps) { 1460 PetscFunctionBegin; 1461 *sweeps = snes->gssweeps; 1462 PetscFunctionReturn(0); 1463 } 1464 1465 #undef __FUNCT__ 1466 #define __FUNCT__ "SNESPicardComputeFunction" 1467 PetscErrorCode SNESPicardComputeFunction(SNES snes,Vec x,Vec f,void *ctx) 1468 { 1469 PetscErrorCode ierr; 1470 void *functx,*jacctx; 1471 1472 PetscFunctionBegin; 1473 ierr = SNESGetFunction(snes,PETSC_NULL,PETSC_NULL,&functx);CHKERRQ(ierr); 1474 ierr = SNESGetJacobian(snes,PETSC_NULL,PETSC_NULL,PETSC_NULL,&jacctx);CHKERRQ(ierr); 1475 /* A(x)*x - b(x) */ 1476 ierr = (*snes->ops->computepjacobian)(snes,x,&snes->jacobian,&snes->jacobian_pre,&snes->matstruct,jacctx);CHKERRQ(ierr); 1477 ierr = (*snes->ops->computepfunction)(snes,x,f,functx);CHKERRQ(ierr); 1478 ierr = VecView(x,PETSC_VIEWER_BINARY_WORLD);CHKERRQ(ierr); 1479 ierr = VecView(f,PETSC_VIEWER_BINARY_WORLD);CHKERRQ(ierr); 1480 ierr = VecScale(f,-1.0);CHKERRQ(ierr); 1481 ierr = MatMultAdd(snes->jacobian_pre,x,f,f);CHKERRQ(ierr); 1482 PetscFunctionReturn(0); 1483 } 1484 1485 #undef __FUNCT__ 1486 #define __FUNCT__ "SNESPicardComputeJacobian" 1487 PetscErrorCode SNESPicardComputeJacobian(SNES snes,Vec x1,Mat *J,Mat *B,MatStructure *flag,void *ctx) 1488 { 1489 PetscFunctionBegin; 1490 *flag = snes->matstruct; 1491 PetscFunctionReturn(0); 1492 } 1493 1494 #undef __FUNCT__ 1495 #define __FUNCT__ "SNESSetPicard" 1496 /*@C 1497 SNESSetPicard - Use SNES to solve the semilinear-system A(x) x = b(x) via a Picard type iteration (Picard linearization) 1498 1499 Logically Collective on SNES 1500 1501 Input Parameters: 1502 + snes - the SNES context 1503 . r - vector to store function value 1504 . func - function evaluation routine 1505 . 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) 1506 . mat - matrix to store A 1507 . mfunc - function to compute matrix value 1508 - ctx - [optional] user-defined context for private data for the 1509 function evaluation routine (may be PETSC_NULL) 1510 1511 Calling sequence of func: 1512 $ func (SNES snes,Vec x,Vec f,void *ctx); 1513 1514 + f - function vector 1515 - ctx - optional user-defined function context 1516 1517 Calling sequence of mfunc: 1518 $ mfunc (SNES snes,Vec x,Mat *jmat,Mat *mat,int *flag,void *ctx); 1519 1520 + x - input vector 1521 . 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(), 1522 normally just pass mat in this location 1523 . mat - form A(x) matrix 1524 . flag - flag indicating information about the preconditioner matrix 1525 structure (same as flag in KSPSetOperators()), one of SAME_NONZERO_PATTERN,DIFFERENT_NONZERO_PATTERN,SAME_PRECONDITIONER 1526 - ctx - [optional] user-defined Jacobian context 1527 1528 Notes: 1529 One can call SNESSetPicard() or SNESSetFunction() (and possibly SNESSetJacobian()) but cannot call both 1530 1531 $ 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} 1532 $ Note that when an exact solver is used this corresponds to the "classic" Picard A(x^{n}) x^{n+1} = b(x^{n}) iteration. 1533 1534 Run with -snes_mf_operator to solve the system with Newton's method using A(x^{n}) to construct the preconditioner. 1535 1536 We implement the defect correction form of the Picard iteration because it converges much more generally when inexact linear solvers are used then 1537 the direct Picard iteration A(x^n) x^{n+1} = b(x^n) 1538 1539 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 1540 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 1541 different please contact us at petsc-dev@mcs.anl.gov and we'll have an entirely new argument :-). 1542 1543 Level: beginner 1544 1545 .keywords: SNES, nonlinear, set, function 1546 1547 .seealso: SNESGetFunction(), SNESSetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESGetPicard(), SNESLineSearchPreCheckPicard() 1548 @*/ 1549 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) 1550 { 1551 PetscErrorCode ierr; 1552 PetscFunctionBegin; 1553 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1554 snes->ops->computepfunction = func; 1555 snes->ops->computepjacobian = mfunc; 1556 ierr = SNESSetFunction(snes,r,SNESPicardComputeFunction,ctx);CHKERRQ(ierr); 1557 ierr = SNESSetJacobian(snes,jmat,mat,SNESPicardComputeJacobian,ctx);CHKERRQ(ierr); 1558 PetscFunctionReturn(0); 1559 } 1560 1561 #undef __FUNCT__ 1562 #define __FUNCT__ "SNESSetComputeInitialGuess" 1563 /*@C 1564 SNESSetComputeInitialGuess - Sets a routine used to compute an initial guess for the problem 1565 1566 Logically Collective on SNES 1567 1568 Input Parameters: 1569 + snes - the SNES context 1570 . func - function evaluation routine 1571 - ctx - [optional] user-defined context for private data for the 1572 function evaluation routine (may be PETSC_NULL) 1573 1574 Calling sequence of func: 1575 $ func (SNES snes,Vec x,void *ctx); 1576 1577 . f - function vector 1578 - ctx - optional user-defined function context 1579 1580 Level: intermediate 1581 1582 .keywords: SNES, nonlinear, set, function 1583 1584 .seealso: SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian() 1585 @*/ 1586 PetscErrorCode SNESSetComputeInitialGuess(SNES snes,PetscErrorCode (*func)(SNES,Vec,void*),void *ctx) 1587 { 1588 PetscFunctionBegin; 1589 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1590 if (func) snes->ops->computeinitialguess = func; 1591 if (ctx) snes->initialguessP = ctx; 1592 PetscFunctionReturn(0); 1593 } 1594 1595 /* --------------------------------------------------------------- */ 1596 #undef __FUNCT__ 1597 #define __FUNCT__ "SNESGetRhs" 1598 /*@C 1599 SNESGetRhs - Gets the vector for solving F(x) = rhs. If rhs is not set 1600 it assumes a zero right hand side. 1601 1602 Logically Collective on SNES 1603 1604 Input Parameter: 1605 . snes - the SNES context 1606 1607 Output Parameter: 1608 . rhs - the right hand side vector or PETSC_NULL if the right hand side vector is null 1609 1610 Level: intermediate 1611 1612 .keywords: SNES, nonlinear, get, function, right hand side 1613 1614 .seealso: SNESGetSolution(), SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESSetFunction() 1615 @*/ 1616 PetscErrorCode SNESGetRhs(SNES snes,Vec *rhs) 1617 { 1618 PetscFunctionBegin; 1619 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1620 PetscValidPointer(rhs,2); 1621 *rhs = snes->vec_rhs; 1622 PetscFunctionReturn(0); 1623 } 1624 1625 #undef __FUNCT__ 1626 #define __FUNCT__ "SNESComputeFunction" 1627 /*@ 1628 SNESComputeFunction - Calls the function that has been set with 1629 SNESSetFunction(). 1630 1631 Collective on SNES 1632 1633 Input Parameters: 1634 + snes - the SNES context 1635 - x - input vector 1636 1637 Output Parameter: 1638 . y - function vector, as set by SNESSetFunction() 1639 1640 Notes: 1641 SNESComputeFunction() is typically used within nonlinear solvers 1642 implementations, so most users would not generally call this routine 1643 themselves. 1644 1645 Level: developer 1646 1647 .keywords: SNES, nonlinear, compute, function 1648 1649 .seealso: SNESSetFunction(), SNESGetFunction() 1650 @*/ 1651 PetscErrorCode SNESComputeFunction(SNES snes,Vec x,Vec y) 1652 { 1653 PetscErrorCode ierr; 1654 DM dm; 1655 SNESDM sdm; 1656 1657 PetscFunctionBegin; 1658 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1659 PetscValidHeaderSpecific(x,VEC_CLASSID,2); 1660 PetscValidHeaderSpecific(y,VEC_CLASSID,3); 1661 PetscCheckSameComm(snes,1,x,2); 1662 PetscCheckSameComm(snes,1,y,3); 1663 VecValidValues(x,2,PETSC_TRUE); 1664 1665 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 1666 ierr = DMSNESGetContext(dm,&sdm);CHKERRQ(ierr); 1667 ierr = PetscLogEventBegin(SNES_FunctionEval,snes,x,y,0);CHKERRQ(ierr); 1668 if (sdm->computefunction) { 1669 PetscStackPush("SNES user function"); 1670 ierr = (*sdm->computefunction)(snes,x,y,sdm->functionctx);CHKERRQ(ierr); 1671 PetscStackPop; 1672 } else if (snes->dm) { 1673 ierr = DMComputeFunction(snes->dm,x,y);CHKERRQ(ierr); 1674 } else if (snes->vec_rhs) { 1675 ierr = MatMult(snes->jacobian, x, y);CHKERRQ(ierr); 1676 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE, "Must call SNESSetFunction() or SNESSetDM() before SNESComputeFunction(), likely called from SNESSolve()."); 1677 if (snes->vec_rhs) { 1678 ierr = VecAXPY(y,-1.0,snes->vec_rhs);CHKERRQ(ierr); 1679 } 1680 snes->nfuncs++; 1681 ierr = PetscLogEventEnd(SNES_FunctionEval,snes,x,y,0);CHKERRQ(ierr); 1682 VecValidValues(y,3,PETSC_FALSE); 1683 PetscFunctionReturn(0); 1684 } 1685 1686 #undef __FUNCT__ 1687 #define __FUNCT__ "SNESComputeGS" 1688 /*@ 1689 SNESComputeGS - Calls the Gauss-Seidel function that has been set with 1690 SNESSetGS(). 1691 1692 Collective on SNES 1693 1694 Input Parameters: 1695 + snes - the SNES context 1696 . x - input vector 1697 - b - rhs vector 1698 1699 Output Parameter: 1700 . x - new solution vector 1701 1702 Notes: 1703 SNESComputeGS() is typically used within composed nonlinear solver 1704 implementations, so most users would not generally call this routine 1705 themselves. 1706 1707 Level: developer 1708 1709 .keywords: SNES, nonlinear, compute, function 1710 1711 .seealso: SNESSetGS(), SNESComputeFunction() 1712 @*/ 1713 PetscErrorCode SNESComputeGS(SNES snes,Vec b,Vec x) 1714 { 1715 PetscErrorCode ierr; 1716 PetscInt i; 1717 DM dm; 1718 SNESDM sdm; 1719 1720 PetscFunctionBegin; 1721 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1722 PetscValidHeaderSpecific(x,VEC_CLASSID,2); 1723 if (b) PetscValidHeaderSpecific(b,VEC_CLASSID,3); 1724 PetscCheckSameComm(snes,1,x,2); 1725 if(b) PetscCheckSameComm(snes,1,b,3); 1726 if (b) VecValidValues(b,2,PETSC_TRUE); 1727 ierr = PetscLogEventBegin(SNES_GSEval,snes,x,b,0);CHKERRQ(ierr); 1728 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 1729 ierr = DMSNESGetContext(dm,&sdm);CHKERRQ(ierr); 1730 if (sdm->computegs) { 1731 for (i = 0; i < snes->gssweeps; i++) { 1732 PetscStackPush("SNES user GS"); 1733 ierr = (*sdm->computegs)(snes,x,b,sdm->gsctx);CHKERRQ(ierr); 1734 PetscStackPop; 1735 } 1736 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE, "Must call SNESSetGS() before SNESComputeGS(), likely called from SNESSolve()."); 1737 ierr = PetscLogEventEnd(SNES_GSEval,snes,x,b,0);CHKERRQ(ierr); 1738 VecValidValues(x,3,PETSC_FALSE); 1739 PetscFunctionReturn(0); 1740 } 1741 1742 1743 #undef __FUNCT__ 1744 #define __FUNCT__ "SNESComputeJacobian" 1745 /*@ 1746 SNESComputeJacobian - Computes the Jacobian matrix that has been 1747 set with SNESSetJacobian(). 1748 1749 Collective on SNES and Mat 1750 1751 Input Parameters: 1752 + snes - the SNES context 1753 - x - input vector 1754 1755 Output Parameters: 1756 + A - Jacobian matrix 1757 . B - optional preconditioning matrix 1758 - flag - flag indicating matrix structure (one of, SAME_NONZERO_PATTERN,DIFFERENT_NONZERO_PATTERN,SAME_PRECONDITIONER) 1759 1760 Options Database Keys: 1761 + -snes_lag_preconditioner <lag> 1762 . -snes_lag_jacobian <lag> 1763 . -snes_compare_explicit - Compare the computed Jacobian to the finite difference Jacobian and output the differences 1764 . -snes_compare_explicit_draw - Compare the computed Jacobian to the finite difference Jacobian and draw the result 1765 . -snes_compare_explicit_contour - Compare the computed Jacobian to the finite difference Jacobian and draw a contour plot with the result 1766 . -snes_compare_operator - Make the comparison options above use the operator instead of the preconditioning matrix 1767 . -snes_compare_coloring - Compute the finite differece Jacobian using coloring and display norms of difference 1768 . -snes_compare_coloring_display - Compute the finite differece Jacobian using coloring and display verbose differences 1769 . -snes_compare_coloring_threshold - Display only those matrix entries that differ by more than a given threshold 1770 . -snes_compare_coloring_threshold_atol - Absolute tolerance for difference in matrix entries to be displayed by -snes_compare_coloring_threshold 1771 . -snes_compare_coloring_threshold_rtol - Relative tolerance for difference in matrix entries to be displayed by -snes_compare_coloring_threshold 1772 . -snes_compare_coloring_draw - Compute the finite differece Jacobian using coloring and draw differences 1773 - -snes_compare_coloring_draw_contour - Compute the finite differece Jacobian using coloring and show contours of matrices and differences 1774 1775 1776 Notes: 1777 Most users should not need to explicitly call this routine, as it 1778 is used internally within the nonlinear solvers. 1779 1780 See KSPSetOperators() for important information about setting the 1781 flag parameter. 1782 1783 Level: developer 1784 1785 .keywords: SNES, compute, Jacobian, matrix 1786 1787 .seealso: SNESSetJacobian(), KSPSetOperators(), MatStructure, SNESSetLagPreconditioner(), SNESSetLagJacobian() 1788 @*/ 1789 PetscErrorCode SNESComputeJacobian(SNES snes,Vec X,Mat *A,Mat *B,MatStructure *flg) 1790 { 1791 PetscErrorCode ierr; 1792 PetscBool flag; 1793 DM dm; 1794 SNESDM sdm; 1795 1796 PetscFunctionBegin; 1797 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1798 PetscValidHeaderSpecific(X,VEC_CLASSID,2); 1799 PetscValidPointer(flg,5); 1800 PetscCheckSameComm(snes,1,X,2); 1801 VecValidValues(X,2,PETSC_TRUE); 1802 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 1803 ierr = DMSNESGetContext(dm,&sdm);CHKERRQ(ierr); 1804 if (!sdm->computejacobian) SETERRQ(((PetscObject)snes)->comm,PETSC_ERR_USER,"Must call SNESSetJacobian(), DMSNESSetJacobian(), DMDASNESSetJacobianLocal(), etc"); 1805 1806 /* make sure that MatAssemblyBegin/End() is called on A matrix if it is matrix free */ 1807 1808 if (snes->lagjacobian == -2) { 1809 snes->lagjacobian = -1; 1810 ierr = PetscInfo(snes,"Recomputing Jacobian/preconditioner because lag is -2 (means compute Jacobian, but then never again) \n");CHKERRQ(ierr); 1811 } else if (snes->lagjacobian == -1) { 1812 *flg = SAME_PRECONDITIONER; 1813 ierr = PetscInfo(snes,"Reusing Jacobian/preconditioner because lag is -1\n");CHKERRQ(ierr); 1814 ierr = PetscTypeCompare((PetscObject)*A,MATMFFD,&flag);CHKERRQ(ierr); 1815 if (flag) { 1816 ierr = MatAssemblyBegin(*A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1817 ierr = MatAssemblyEnd(*A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1818 } 1819 PetscFunctionReturn(0); 1820 } else if (snes->lagjacobian > 1 && snes->iter % snes->lagjacobian) { 1821 *flg = SAME_PRECONDITIONER; 1822 ierr = PetscInfo2(snes,"Reusing Jacobian/preconditioner because lag is %D and SNES iteration is %D\n",snes->lagjacobian,snes->iter);CHKERRQ(ierr); 1823 ierr = PetscTypeCompare((PetscObject)*A,MATMFFD,&flag);CHKERRQ(ierr); 1824 if (flag) { 1825 ierr = MatAssemblyBegin(*A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1826 ierr = MatAssemblyEnd(*A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1827 } 1828 PetscFunctionReturn(0); 1829 } 1830 1831 *flg = DIFFERENT_NONZERO_PATTERN; 1832 ierr = PetscLogEventBegin(SNES_JacobianEval,snes,X,*A,*B);CHKERRQ(ierr); 1833 PetscStackPush("SNES user Jacobian function"); 1834 ierr = (*sdm->computejacobian)(snes,X,A,B,flg,sdm->jacobianctx);CHKERRQ(ierr); 1835 PetscStackPop; 1836 ierr = PetscLogEventEnd(SNES_JacobianEval,snes,X,*A,*B);CHKERRQ(ierr); 1837 1838 if (snes->lagpreconditioner == -2) { 1839 ierr = PetscInfo(snes,"Rebuilding preconditioner exactly once since lag is -2\n");CHKERRQ(ierr); 1840 snes->lagpreconditioner = -1; 1841 } else if (snes->lagpreconditioner == -1) { 1842 *flg = SAME_PRECONDITIONER; 1843 ierr = PetscInfo(snes,"Reusing preconditioner because lag is -1\n");CHKERRQ(ierr); 1844 } else if (snes->lagpreconditioner > 1 && snes->iter % snes->lagpreconditioner) { 1845 *flg = SAME_PRECONDITIONER; 1846 ierr = PetscInfo2(snes,"Reusing preconditioner because lag is %D and SNES iteration is %D\n",snes->lagpreconditioner,snes->iter);CHKERRQ(ierr); 1847 } 1848 1849 /* make sure user returned a correct Jacobian and preconditioner */ 1850 /* PetscValidHeaderSpecific(*A,MAT_CLASSID,3); 1851 PetscValidHeaderSpecific(*B,MAT_CLASSID,4); */ 1852 { 1853 PetscBool flag = PETSC_FALSE,flag_draw = PETSC_FALSE,flag_contour = PETSC_FALSE,flag_operator = PETSC_FALSE; 1854 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_explicit",&flag,PETSC_NULL);CHKERRQ(ierr); 1855 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_explicit_draw",&flag_draw,PETSC_NULL);CHKERRQ(ierr); 1856 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_explicit_draw_contour",&flag_contour,PETSC_NULL);CHKERRQ(ierr); 1857 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_operator",&flag_operator,PETSC_NULL);CHKERRQ(ierr); 1858 if (flag || flag_draw || flag_contour) { 1859 Mat Bexp_mine = PETSC_NULL,Bexp,FDexp; 1860 MatStructure mstruct; 1861 PetscViewer vdraw,vstdout; 1862 PetscBool flg; 1863 if (flag_operator) { 1864 ierr = MatComputeExplicitOperator(*A,&Bexp_mine);CHKERRQ(ierr); 1865 Bexp = Bexp_mine; 1866 } else { 1867 /* See if the preconditioning matrix can be viewed and added directly */ 1868 ierr = PetscTypeCompareAny((PetscObject)*B,&flg,MATSEQAIJ,MATMPIAIJ,MATSEQDENSE,MATMPIDENSE,MATSEQBAIJ,MATMPIBAIJ,MATSEQSBAIJ,MATMPIBAIJ,"");CHKERRQ(ierr); 1869 if (flg) Bexp = *B; 1870 else { 1871 /* If the "preconditioning" matrix is itself MATSHELL or some other type without direct support */ 1872 ierr = MatComputeExplicitOperator(*B,&Bexp_mine);CHKERRQ(ierr); 1873 Bexp = Bexp_mine; 1874 } 1875 } 1876 ierr = MatConvert(Bexp,MATSAME,MAT_INITIAL_MATRIX,&FDexp);CHKERRQ(ierr); 1877 ierr = SNESDefaultComputeJacobian(snes,X,&FDexp,&FDexp,&mstruct,NULL);CHKERRQ(ierr); 1878 ierr = PetscViewerASCIIGetStdout(((PetscObject)snes)->comm,&vstdout);CHKERRQ(ierr); 1879 if (flag_draw || flag_contour) { 1880 ierr = PetscViewerDrawOpen(((PetscObject)snes)->comm,0,"Explicit Jacobians",PETSC_DECIDE,PETSC_DECIDE,300,300,&vdraw);CHKERRQ(ierr); 1881 if (flag_contour) {ierr = PetscViewerPushFormat(vdraw,PETSC_VIEWER_DRAW_CONTOUR);CHKERRQ(ierr);} 1882 } else vdraw = PETSC_NULL; 1883 ierr = PetscViewerASCIIPrintf(vstdout,"Explicit %s\n",flag_operator?"Jacobian":"preconditioning Jacobian");CHKERRQ(ierr); 1884 if (flag) {ierr = MatView(Bexp,vstdout);CHKERRQ(ierr);} 1885 if (vdraw) {ierr = MatView(Bexp,vdraw);CHKERRQ(ierr);} 1886 ierr = PetscViewerASCIIPrintf(vstdout,"Finite difference Jacobian\n");CHKERRQ(ierr); 1887 if (flag) {ierr = MatView(FDexp,vstdout);CHKERRQ(ierr);} 1888 if (vdraw) {ierr = MatView(FDexp,vdraw);CHKERRQ(ierr);} 1889 ierr = MatAYPX(FDexp,-1.0,Bexp,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 1890 ierr = PetscViewerASCIIPrintf(vstdout,"User-provided matrix minus finite difference Jacobian\n");CHKERRQ(ierr); 1891 if (flag) {ierr = MatView(FDexp,vstdout);CHKERRQ(ierr);} 1892 if (vdraw) { /* Always use contour for the difference */ 1893 ierr = PetscViewerPushFormat(vdraw,PETSC_VIEWER_DRAW_CONTOUR);CHKERRQ(ierr); 1894 ierr = MatView(FDexp,vdraw);CHKERRQ(ierr); 1895 ierr = PetscViewerPopFormat(vdraw);CHKERRQ(ierr); 1896 } 1897 if (flag_contour) {ierr = PetscViewerPopFormat(vdraw);CHKERRQ(ierr);} 1898 ierr = PetscViewerDestroy(&vdraw);CHKERRQ(ierr); 1899 ierr = MatDestroy(&Bexp_mine);CHKERRQ(ierr); 1900 ierr = MatDestroy(&FDexp);CHKERRQ(ierr); 1901 } 1902 } 1903 { 1904 PetscBool flag = PETSC_FALSE,flag_display = PETSC_FALSE,flag_draw = PETSC_FALSE,flag_contour = PETSC_FALSE,flag_threshold = PETSC_FALSE; 1905 PetscReal threshold_atol = PETSC_SQRT_MACHINE_EPSILON,threshold_rtol = 10*PETSC_SQRT_MACHINE_EPSILON; 1906 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_coloring",&flag,PETSC_NULL);CHKERRQ(ierr); 1907 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_coloring_display",&flag_display,PETSC_NULL);CHKERRQ(ierr); 1908 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_coloring_draw",&flag_draw,PETSC_NULL);CHKERRQ(ierr); 1909 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_coloring_draw_contour",&flag_contour,PETSC_NULL);CHKERRQ(ierr); 1910 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_coloring_threshold",&flag_threshold,PETSC_NULL);CHKERRQ(ierr); 1911 ierr = PetscOptionsGetReal(((PetscObject)snes)->prefix,"-snes_compare_coloring_threshold_rtol",&threshold_rtol,PETSC_NULL);CHKERRQ(ierr); 1912 ierr = PetscOptionsGetReal(((PetscObject)snes)->prefix,"-snes_compare_coloring_threshold_atol",&threshold_atol,PETSC_NULL);CHKERRQ(ierr); 1913 if (flag || flag_display || flag_draw || flag_contour || flag_threshold) { 1914 Mat Bfd; 1915 MatStructure mstruct; 1916 PetscViewer vdraw,vstdout; 1917 ISColoring iscoloring; 1918 MatFDColoring matfdcoloring; 1919 PetscErrorCode (*func)(SNES,Vec,Vec,void*); 1920 void *funcctx; 1921 PetscReal norm1,norm2,normmax; 1922 1923 ierr = MatDuplicate(*B,MAT_DO_NOT_COPY_VALUES,&Bfd);CHKERRQ(ierr); 1924 ierr = MatGetColoring(Bfd,MATCOLORINGSL,&iscoloring);CHKERRQ(ierr); 1925 ierr = MatFDColoringCreate(Bfd,iscoloring,&matfdcoloring);CHKERRQ(ierr); 1926 ierr = ISColoringDestroy(&iscoloring);CHKERRQ(ierr); 1927 1928 /* This method of getting the function is currently unreliable since it doesn't work for DM local functions. */ 1929 ierr = SNESGetFunction(snes,PETSC_NULL,&func,&funcctx);CHKERRQ(ierr); 1930 ierr = MatFDColoringSetFunction(matfdcoloring,(PetscErrorCode(*)(void))func,funcctx);CHKERRQ(ierr); 1931 ierr = PetscObjectSetOptionsPrefix((PetscObject)matfdcoloring,((PetscObject)snes)->prefix);CHKERRQ(ierr); 1932 ierr = PetscObjectAppendOptionsPrefix((PetscObject)matfdcoloring,"coloring_");CHKERRQ(ierr); 1933 ierr = MatFDColoringSetFromOptions(matfdcoloring);CHKERRQ(ierr); 1934 ierr = MatFDColoringApply(Bfd,matfdcoloring,X,&mstruct,snes);CHKERRQ(ierr); 1935 ierr = MatFDColoringDestroy(&matfdcoloring);CHKERRQ(ierr); 1936 1937 ierr = PetscViewerASCIIGetStdout(((PetscObject)snes)->comm,&vstdout);CHKERRQ(ierr); 1938 if (flag_draw || flag_contour) { 1939 ierr = PetscViewerDrawOpen(((PetscObject)snes)->comm,0,"Colored Jacobians",PETSC_DECIDE,PETSC_DECIDE,300,300,&vdraw);CHKERRQ(ierr); 1940 if (flag_contour) {ierr = PetscViewerPushFormat(vdraw,PETSC_VIEWER_DRAW_CONTOUR);CHKERRQ(ierr);} 1941 } else vdraw = PETSC_NULL; 1942 ierr = PetscViewerASCIIPrintf(vstdout,"Explicit preconditioning Jacobian\n");CHKERRQ(ierr); 1943 if (flag_display) {ierr = MatView(*B,vstdout);CHKERRQ(ierr);} 1944 if (vdraw) {ierr = MatView(*B,vdraw);CHKERRQ(ierr);} 1945 ierr = PetscViewerASCIIPrintf(vstdout,"Colored Finite difference Jacobian\n");CHKERRQ(ierr); 1946 if (flag_display) {ierr = MatView(Bfd,vstdout);CHKERRQ(ierr);} 1947 if (vdraw) {ierr = MatView(Bfd,vdraw);CHKERRQ(ierr);} 1948 ierr = MatAYPX(Bfd,-1.0,*B,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 1949 ierr = MatNorm(Bfd,NORM_1,&norm1);CHKERRQ(ierr); 1950 ierr = MatNorm(Bfd,NORM_FROBENIUS,&norm2);CHKERRQ(ierr); 1951 ierr = MatNorm(Bfd,NORM_MAX,&normmax);CHKERRQ(ierr); 1952 ierr = PetscViewerASCIIPrintf(vstdout,"User-provided matrix minus finite difference Jacobian, norm1=%G normFrob=%G normmax=%G\n",norm1,norm2,normmax);CHKERRQ(ierr); 1953 if (flag_display) {ierr = MatView(Bfd,vstdout);CHKERRQ(ierr);} 1954 if (vdraw) { /* Always use contour for the difference */ 1955 ierr = PetscViewerPushFormat(vdraw,PETSC_VIEWER_DRAW_CONTOUR);CHKERRQ(ierr); 1956 ierr = MatView(Bfd,vdraw);CHKERRQ(ierr); 1957 ierr = PetscViewerPopFormat(vdraw);CHKERRQ(ierr); 1958 } 1959 if (flag_contour) {ierr = PetscViewerPopFormat(vdraw);CHKERRQ(ierr);} 1960 1961 if (flag_threshold) { 1962 PetscInt bs,rstart,rend,i; 1963 ierr = MatGetBlockSize(*B,&bs);CHKERRQ(ierr); 1964 ierr = MatGetOwnershipRange(*B,&rstart,&rend);CHKERRQ(ierr); 1965 for (i=rstart; i<rend; i++) { 1966 const PetscScalar *ba,*ca; 1967 const PetscInt *bj,*cj; 1968 PetscInt bn,cn,j,maxentrycol = -1,maxdiffcol = -1,maxrdiffcol = -1; 1969 PetscReal maxentry = 0,maxdiff = 0,maxrdiff = 0; 1970 ierr = MatGetRow(*B,i,&bn,&bj,&ba);CHKERRQ(ierr); 1971 ierr = MatGetRow(Bfd,i,&cn,&cj,&ca);CHKERRQ(ierr); 1972 if (bn != cn) SETERRQ(((PetscObject)*A)->comm,PETSC_ERR_PLIB,"Unexpected different nonzero pattern in -snes_compare_coloring_threshold"); 1973 for (j=0; j<bn; j++) { 1974 PetscReal rdiff = PetscAbsScalar(ca[j]) / (threshold_atol + threshold_rtol*PetscAbsScalar(ba[j])); 1975 if (PetscAbsScalar(ba[j]) > PetscAbs(maxentry)) { 1976 maxentrycol = bj[j]; 1977 maxentry = PetscRealPart(ba[j]); 1978 } 1979 if (PetscAbsScalar(ca[j]) > PetscAbs(maxdiff)) { 1980 maxdiffcol = bj[j]; 1981 maxdiff = PetscRealPart(ca[j]); 1982 } 1983 if (rdiff > maxrdiff) { 1984 maxrdiffcol = bj[j]; 1985 maxrdiff = rdiff; 1986 } 1987 } 1988 if (maxrdiff > 1) { 1989 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); 1990 for (j=0; j<bn; j++) { 1991 PetscReal rdiff; 1992 rdiff = PetscAbsScalar(ca[j]) / (threshold_atol + threshold_rtol*PetscAbsScalar(ba[j])); 1993 if (rdiff > 1) { 1994 ierr = PetscViewerASCIIPrintf(vstdout," (%D,%G:%G)",bj[j],PetscRealPart(ba[j]),PetscRealPart(ca[j]));CHKERRQ(ierr); 1995 } 1996 } 1997 ierr = PetscViewerASCIIPrintf(vstdout,"\n",i,maxentry,maxdiff,maxrdiff);CHKERRQ(ierr); 1998 } 1999 ierr = MatRestoreRow(*B,i,&bn,&bj,&ba);CHKERRQ(ierr); 2000 ierr = MatRestoreRow(Bfd,i,&cn,&cj,&ca);CHKERRQ(ierr); 2001 } 2002 } 2003 ierr = PetscViewerDestroy(&vdraw);CHKERRQ(ierr); 2004 ierr = MatDestroy(&Bfd);CHKERRQ(ierr); 2005 } 2006 } 2007 PetscFunctionReturn(0); 2008 } 2009 2010 #undef __FUNCT__ 2011 #define __FUNCT__ "SNESSetJacobian" 2012 /*@C 2013 SNESSetJacobian - Sets the function to compute Jacobian as well as the 2014 location to store the matrix. 2015 2016 Logically Collective on SNES and Mat 2017 2018 Input Parameters: 2019 + snes - the SNES context 2020 . A - Jacobian matrix 2021 . B - preconditioner matrix (usually same as the Jacobian) 2022 . func - Jacobian evaluation routine (if PETSC_NULL then SNES retains any previously set value) 2023 - ctx - [optional] user-defined context for private data for the 2024 Jacobian evaluation routine (may be PETSC_NULL) (if PETSC_NULL then SNES retains any previously set value) 2025 2026 Calling sequence of func: 2027 $ func (SNES snes,Vec x,Mat *A,Mat *B,int *flag,void *ctx); 2028 2029 + x - input vector 2030 . A - Jacobian matrix 2031 . B - preconditioner matrix, usually the same as A 2032 . flag - flag indicating information about the preconditioner matrix 2033 structure (same as flag in KSPSetOperators()), one of SAME_NONZERO_PATTERN,DIFFERENT_NONZERO_PATTERN,SAME_PRECONDITIONER 2034 - ctx - [optional] user-defined Jacobian context 2035 2036 Notes: 2037 See KSPSetOperators() for important information about setting the flag 2038 output parameter in the routine func(). Be sure to read this information! 2039 2040 The routine func() takes Mat * as the matrix arguments rather than Mat. 2041 This allows the Jacobian evaluation routine to replace A and/or B with a 2042 completely new new matrix structure (not just different matrix elements) 2043 when appropriate, for instance, if the nonzero structure is changing 2044 throughout the global iterations. 2045 2046 If the A matrix and B matrix are different you must call MatAssemblyBegin/End() on 2047 each matrix. 2048 2049 If using SNESDefaultComputeJacobianColor() to assemble a Jacobian, the ctx argument 2050 must be a MatFDColoring. 2051 2052 Other defect-correction schemes can be used by computing a different matrix in place of the Jacobian. One common 2053 example is to use the "Picard linearization" which only differentiates through the highest order parts of each term. 2054 2055 Level: beginner 2056 2057 .keywords: SNES, nonlinear, set, Jacobian, matrix 2058 2059 .seealso: KSPSetOperators(), SNESSetFunction(), MatMFFDComputeJacobian(), SNESDefaultComputeJacobianColor(), MatStructure 2060 @*/ 2061 PetscErrorCode SNESSetJacobian(SNES snes,Mat A,Mat B,PetscErrorCode (*func)(SNES,Vec,Mat*,Mat*,MatStructure*,void*),void *ctx) 2062 { 2063 PetscErrorCode ierr; 2064 DM dm; 2065 2066 PetscFunctionBegin; 2067 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2068 if (A) PetscValidHeaderSpecific(A,MAT_CLASSID,2); 2069 if (B) PetscValidHeaderSpecific(B,MAT_CLASSID,3); 2070 if (A) PetscCheckSameComm(snes,1,A,2); 2071 if (B) PetscCheckSameComm(snes,1,B,3); 2072 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 2073 ierr = DMSNESSetJacobian(dm,func,ctx);CHKERRQ(ierr); 2074 if (A) { 2075 ierr = PetscObjectReference((PetscObject)A);CHKERRQ(ierr); 2076 ierr = MatDestroy(&snes->jacobian);CHKERRQ(ierr); 2077 snes->jacobian = A; 2078 } 2079 if (B) { 2080 ierr = PetscObjectReference((PetscObject)B);CHKERRQ(ierr); 2081 ierr = MatDestroy(&snes->jacobian_pre);CHKERRQ(ierr); 2082 snes->jacobian_pre = B; 2083 } 2084 PetscFunctionReturn(0); 2085 } 2086 2087 #undef __FUNCT__ 2088 #define __FUNCT__ "SNESGetJacobian" 2089 /*@C 2090 SNESGetJacobian - Returns the Jacobian matrix and optionally the user 2091 provided context for evaluating the Jacobian. 2092 2093 Not Collective, but Mat object will be parallel if SNES object is 2094 2095 Input Parameter: 2096 . snes - the nonlinear solver context 2097 2098 Output Parameters: 2099 + A - location to stash Jacobian matrix (or PETSC_NULL) 2100 . B - location to stash preconditioner matrix (or PETSC_NULL) 2101 . func - location to put Jacobian function (or PETSC_NULL) 2102 - ctx - location to stash Jacobian ctx (or PETSC_NULL) 2103 2104 Level: advanced 2105 2106 .seealso: SNESSetJacobian(), SNESComputeJacobian() 2107 @*/ 2108 PetscErrorCode SNESGetJacobian(SNES snes,Mat *A,Mat *B,PetscErrorCode (**func)(SNES,Vec,Mat*,Mat*,MatStructure*,void*),void **ctx) 2109 { 2110 PetscErrorCode ierr; 2111 DM dm; 2112 SNESDM sdm; 2113 2114 PetscFunctionBegin; 2115 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2116 if (A) *A = snes->jacobian; 2117 if (B) *B = snes->jacobian_pre; 2118 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 2119 ierr = DMSNESGetContext(dm,&sdm);CHKERRQ(ierr); 2120 if (func) *func = sdm->computejacobian; 2121 if (ctx) *ctx = sdm->jacobianctx; 2122 PetscFunctionReturn(0); 2123 } 2124 2125 /* ----- Routines to initialize and destroy a nonlinear solver ---- */ 2126 2127 #undef __FUNCT__ 2128 #define __FUNCT__ "SNESSetUp" 2129 /*@ 2130 SNESSetUp - Sets up the internal data structures for the later use 2131 of a nonlinear solver. 2132 2133 Collective on SNES 2134 2135 Input Parameters: 2136 . snes - the SNES context 2137 2138 Notes: 2139 For basic use of the SNES solvers the user need not explicitly call 2140 SNESSetUp(), since these actions will automatically occur during 2141 the call to SNESSolve(). However, if one wishes to control this 2142 phase separately, SNESSetUp() should be called after SNESCreate() 2143 and optional routines of the form SNESSetXXX(), but before SNESSolve(). 2144 2145 Level: advanced 2146 2147 .keywords: SNES, nonlinear, setup 2148 2149 .seealso: SNESCreate(), SNESSolve(), SNESDestroy() 2150 @*/ 2151 PetscErrorCode SNESSetUp(SNES snes) 2152 { 2153 PetscErrorCode ierr; 2154 DM dm; 2155 SNESDM sdm; 2156 2157 PetscFunctionBegin; 2158 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2159 if (snes->setupcalled) PetscFunctionReturn(0); 2160 2161 if (!((PetscObject)snes)->type_name) { 2162 ierr = SNESSetType(snes,SNESLS);CHKERRQ(ierr); 2163 } 2164 2165 ierr = SNESGetFunction(snes,&snes->vec_func,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr); 2166 if (snes->vec_func == snes->vec_sol) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_IDN,"Solution vector cannot be function vector"); 2167 if (snes->vec_rhs == snes->vec_sol) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_IDN,"Solution vector cannot be right hand side vector"); 2168 2169 if (!snes->vec_sol_update /* && snes->vec_sol */) { 2170 ierr = VecDuplicate(snes->vec_sol,&snes->vec_sol_update);CHKERRQ(ierr); 2171 ierr = PetscLogObjectParent(snes,snes->vec_sol_update);CHKERRQ(ierr); 2172 } 2173 2174 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 2175 ierr = DMSNESSetUpLegacy(dm);CHKERRQ(ierr); /* To be removed when function routines are taken out of the DM package */ 2176 ierr = DMSNESGetContext(dm,&sdm);CHKERRQ(ierr); 2177 if (!sdm->computefunction) SETERRQ(((PetscObject)snes)->comm,PETSC_ERR_ARG_WRONGSTATE,"Must provide a residual function with SNESSetFunction(), DMSNESSetFunction(), DMDASNESSetFunctionLocal(), etc"); 2178 if (!snes->vec_func) { 2179 ierr = DMCreateGlobalVector(dm,&snes->vec_func);CHKERRQ(ierr); 2180 } 2181 2182 if (!snes->ksp) {ierr = SNESGetKSP(snes, &snes->ksp);CHKERRQ(ierr);} 2183 2184 if (snes->ops->usercompute && !snes->user) { 2185 ierr = (*snes->ops->usercompute)(snes,(void**)&snes->user);CHKERRQ(ierr); 2186 } 2187 2188 if (snes->ops->setup) { 2189 ierr = (*snes->ops->setup)(snes);CHKERRQ(ierr); 2190 } 2191 2192 snes->setupcalled = PETSC_TRUE; 2193 PetscFunctionReturn(0); 2194 } 2195 2196 #undef __FUNCT__ 2197 #define __FUNCT__ "SNESReset" 2198 /*@ 2199 SNESReset - Resets a SNES context to the snessetupcalled = 0 state and removes any allocated Vecs and Mats 2200 2201 Collective on SNES 2202 2203 Input Parameter: 2204 . snes - iterative context obtained from SNESCreate() 2205 2206 Level: intermediate 2207 2208 Notes: Also calls the application context destroy routine set with SNESSetComputeApplicationContext() 2209 2210 .keywords: SNES, destroy 2211 2212 .seealso: SNESCreate(), SNESSetUp(), SNESSolve() 2213 @*/ 2214 PetscErrorCode SNESReset(SNES snes) 2215 { 2216 PetscErrorCode ierr; 2217 2218 PetscFunctionBegin; 2219 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2220 if (snes->ops->userdestroy && snes->user) { 2221 ierr = (*snes->ops->userdestroy)((void**)&snes->user);CHKERRQ(ierr); 2222 snes->user = PETSC_NULL; 2223 } 2224 if (snes->pc) { 2225 ierr = SNESReset(snes->pc);CHKERRQ(ierr); 2226 } 2227 2228 if (snes->ops->reset) { 2229 ierr = (*snes->ops->reset)(snes);CHKERRQ(ierr); 2230 } 2231 if (snes->ksp) {ierr = KSPReset(snes->ksp);CHKERRQ(ierr);} 2232 ierr = VecDestroy(&snes->vec_rhs);CHKERRQ(ierr); 2233 ierr = VecDestroy(&snes->vec_sol);CHKERRQ(ierr); 2234 ierr = VecDestroy(&snes->vec_sol_update);CHKERRQ(ierr); 2235 ierr = VecDestroy(&snes->vec_func);CHKERRQ(ierr); 2236 ierr = MatDestroy(&snes->jacobian);CHKERRQ(ierr); 2237 ierr = MatDestroy(&snes->jacobian_pre);CHKERRQ(ierr); 2238 ierr = VecDestroyVecs(snes->nwork,&snes->work);CHKERRQ(ierr); 2239 ierr = VecDestroyVecs(snes->nvwork,&snes->vwork);CHKERRQ(ierr); 2240 snes->nwork = snes->nvwork = 0; 2241 snes->setupcalled = PETSC_FALSE; 2242 PetscFunctionReturn(0); 2243 } 2244 2245 #undef __FUNCT__ 2246 #define __FUNCT__ "SNESDestroy" 2247 /*@ 2248 SNESDestroy - Destroys the nonlinear solver context that was created 2249 with SNESCreate(). 2250 2251 Collective on SNES 2252 2253 Input Parameter: 2254 . snes - the SNES context 2255 2256 Level: beginner 2257 2258 .keywords: SNES, nonlinear, destroy 2259 2260 .seealso: SNESCreate(), SNESSolve() 2261 @*/ 2262 PetscErrorCode SNESDestroy(SNES *snes) 2263 { 2264 PetscErrorCode ierr; 2265 2266 PetscFunctionBegin; 2267 if (!*snes) PetscFunctionReturn(0); 2268 PetscValidHeaderSpecific((*snes),SNES_CLASSID,1); 2269 if (--((PetscObject)(*snes))->refct > 0) {*snes = 0; PetscFunctionReturn(0);} 2270 2271 ierr = SNESReset((*snes));CHKERRQ(ierr); 2272 ierr = SNESDestroy(&(*snes)->pc);CHKERRQ(ierr); 2273 2274 /* if memory was published with AMS then destroy it */ 2275 ierr = PetscObjectDepublish((*snes));CHKERRQ(ierr); 2276 if ((*snes)->ops->destroy) {ierr = (*((*snes))->ops->destroy)((*snes));CHKERRQ(ierr);} 2277 2278 ierr = DMDestroy(&(*snes)->dm);CHKERRQ(ierr); 2279 ierr = KSPDestroy(&(*snes)->ksp);CHKERRQ(ierr); 2280 2281 ierr = PetscFree((*snes)->kspconvctx);CHKERRQ(ierr); 2282 if ((*snes)->ops->convergeddestroy) { 2283 ierr = (*(*snes)->ops->convergeddestroy)((*snes)->cnvP);CHKERRQ(ierr); 2284 } 2285 if ((*snes)->conv_malloc) { 2286 ierr = PetscFree((*snes)->conv_hist);CHKERRQ(ierr); 2287 ierr = PetscFree((*snes)->conv_hist_its);CHKERRQ(ierr); 2288 } 2289 ierr = PetscViewerDestroy(&(*snes)->ls_monitor);CHKERRQ(ierr); 2290 ierr = SNESMonitorCancel((*snes));CHKERRQ(ierr); 2291 ierr = PetscHeaderDestroy(snes);CHKERRQ(ierr); 2292 PetscFunctionReturn(0); 2293 } 2294 2295 /* ----------- Routines to set solver parameters ---------- */ 2296 2297 #undef __FUNCT__ 2298 #define __FUNCT__ "SNESSetLagPreconditioner" 2299 /*@ 2300 SNESSetLagPreconditioner - Determines when the preconditioner is rebuilt in the nonlinear solve. 2301 2302 Logically Collective on SNES 2303 2304 Input Parameters: 2305 + snes - the SNES context 2306 - lag - -1 indicates NEVER rebuild, 1 means rebuild every time the Jacobian is computed within a single nonlinear solve, 2 means every second time 2307 the Jacobian is built etc. -2 indicates rebuild preconditioner at next chance but then never rebuild after that 2308 2309 Options Database Keys: 2310 . -snes_lag_preconditioner <lag> 2311 2312 Notes: 2313 The default is 1 2314 The preconditioner is ALWAYS built in the first iteration of a nonlinear solve unless lag is -1 2315 If -1 is used before the very first nonlinear solve the preconditioner is still built because there is no previous preconditioner to use 2316 2317 Level: intermediate 2318 2319 .keywords: SNES, nonlinear, set, convergence, tolerances 2320 2321 .seealso: SNESSetTrustRegionTolerance(), SNESGetLagPreconditioner(), SNESSetLagJacobian(), SNESGetLagJacobian() 2322 2323 @*/ 2324 PetscErrorCode SNESSetLagPreconditioner(SNES snes,PetscInt lag) 2325 { 2326 PetscFunctionBegin; 2327 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2328 if (lag < -2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Lag must be -2, -1, 1 or greater"); 2329 if (!lag) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Lag cannot be 0"); 2330 PetscValidLogicalCollectiveInt(snes,lag,2); 2331 snes->lagpreconditioner = lag; 2332 PetscFunctionReturn(0); 2333 } 2334 2335 #undef __FUNCT__ 2336 #define __FUNCT__ "SNESSetGridSequence" 2337 /*@ 2338 SNESSetGridSequence - sets the number of steps of grid sequencing that SNES does 2339 2340 Logically Collective on SNES 2341 2342 Input Parameters: 2343 + snes - the SNES context 2344 - steps - the number of refinements to do, defaults to 0 2345 2346 Options Database Keys: 2347 . -snes_grid_sequence <steps> 2348 2349 Level: intermediate 2350 2351 Notes: 2352 Use SNESGetSolution() to extract the fine grid solution after grid sequencing. 2353 2354 .keywords: SNES, nonlinear, set, convergence, tolerances 2355 2356 .seealso: SNESSetTrustRegionTolerance(), SNESGetLagPreconditioner(), SNESSetLagJacobian(), SNESGetLagJacobian() 2357 2358 @*/ 2359 PetscErrorCode SNESSetGridSequence(SNES snes,PetscInt steps) 2360 { 2361 PetscFunctionBegin; 2362 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2363 PetscValidLogicalCollectiveInt(snes,steps,2); 2364 snes->gridsequence = steps; 2365 PetscFunctionReturn(0); 2366 } 2367 2368 #undef __FUNCT__ 2369 #define __FUNCT__ "SNESGetLagPreconditioner" 2370 /*@ 2371 SNESGetLagPreconditioner - Indicates how often the preconditioner is rebuilt 2372 2373 Not Collective 2374 2375 Input Parameter: 2376 . snes - the SNES context 2377 2378 Output Parameter: 2379 . lag - -1 indicates NEVER rebuild, 1 means rebuild every time the Jacobian is computed within a single nonlinear solve, 2 means every second time 2380 the Jacobian is built etc. -2 indicates rebuild preconditioner at next chance but then never rebuild after that 2381 2382 Options Database Keys: 2383 . -snes_lag_preconditioner <lag> 2384 2385 Notes: 2386 The default is 1 2387 The preconditioner is ALWAYS built in the first iteration of a nonlinear solve unless lag is -1 2388 2389 Level: intermediate 2390 2391 .keywords: SNES, nonlinear, set, convergence, tolerances 2392 2393 .seealso: SNESSetTrustRegionTolerance(), SNESSetLagPreconditioner() 2394 2395 @*/ 2396 PetscErrorCode SNESGetLagPreconditioner(SNES snes,PetscInt *lag) 2397 { 2398 PetscFunctionBegin; 2399 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2400 *lag = snes->lagpreconditioner; 2401 PetscFunctionReturn(0); 2402 } 2403 2404 #undef __FUNCT__ 2405 #define __FUNCT__ "SNESSetLagJacobian" 2406 /*@ 2407 SNESSetLagJacobian - Determines when the Jacobian is rebuilt in the nonlinear solve. See SNESSetLagPreconditioner() for determining how 2408 often the preconditioner is rebuilt. 2409 2410 Logically Collective on SNES 2411 2412 Input Parameters: 2413 + snes - the SNES context 2414 - lag - -1 indicates NEVER rebuild, 1 means rebuild every time the Jacobian is computed within a single nonlinear solve, 2 means every second time 2415 the Jacobian is built etc. -2 means rebuild at next chance but then never again 2416 2417 Options Database Keys: 2418 . -snes_lag_jacobian <lag> 2419 2420 Notes: 2421 The default is 1 2422 The Jacobian is ALWAYS built in the first iteration of a nonlinear solve unless lag is -1 2423 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 2424 at the next Newton step but never again (unless it is reset to another value) 2425 2426 Level: intermediate 2427 2428 .keywords: SNES, nonlinear, set, convergence, tolerances 2429 2430 .seealso: SNESSetTrustRegionTolerance(), SNESGetLagPreconditioner(), SNESSetLagPreconditioner(), SNESGetLagJacobian() 2431 2432 @*/ 2433 PetscErrorCode SNESSetLagJacobian(SNES snes,PetscInt lag) 2434 { 2435 PetscFunctionBegin; 2436 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2437 if (lag < -2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Lag must be -2, -1, 1 or greater"); 2438 if (!lag) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Lag cannot be 0"); 2439 PetscValidLogicalCollectiveInt(snes,lag,2); 2440 snes->lagjacobian = lag; 2441 PetscFunctionReturn(0); 2442 } 2443 2444 #undef __FUNCT__ 2445 #define __FUNCT__ "SNESGetLagJacobian" 2446 /*@ 2447 SNESGetLagJacobian - Indicates how often the Jacobian is rebuilt. See SNESGetLagPreconditioner() to determine when the preconditioner is rebuilt 2448 2449 Not Collective 2450 2451 Input Parameter: 2452 . snes - the SNES context 2453 2454 Output Parameter: 2455 . lag - -1 indicates NEVER rebuild, 1 means rebuild every time the Jacobian is computed within a single nonlinear solve, 2 means every second time 2456 the Jacobian is built etc. 2457 2458 Options Database Keys: 2459 . -snes_lag_jacobian <lag> 2460 2461 Notes: 2462 The default is 1 2463 The jacobian is ALWAYS built in the first iteration of a nonlinear solve unless lag is -1 2464 2465 Level: intermediate 2466 2467 .keywords: SNES, nonlinear, set, convergence, tolerances 2468 2469 .seealso: SNESSetTrustRegionTolerance(), SNESSetLagJacobian(), SNESSetLagPreconditioner(), SNESGetLagPreconditioner() 2470 2471 @*/ 2472 PetscErrorCode SNESGetLagJacobian(SNES snes,PetscInt *lag) 2473 { 2474 PetscFunctionBegin; 2475 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2476 *lag = snes->lagjacobian; 2477 PetscFunctionReturn(0); 2478 } 2479 2480 #undef __FUNCT__ 2481 #define __FUNCT__ "SNESSetTolerances" 2482 /*@ 2483 SNESSetTolerances - Sets various parameters used in convergence tests. 2484 2485 Logically Collective on SNES 2486 2487 Input Parameters: 2488 + snes - the SNES context 2489 . abstol - absolute convergence tolerance 2490 . rtol - relative convergence tolerance 2491 . stol - convergence tolerance in terms of the norm 2492 of the change in the solution between steps 2493 . maxit - maximum number of iterations 2494 - maxf - maximum number of function evaluations 2495 2496 Options Database Keys: 2497 + -snes_atol <abstol> - Sets abstol 2498 . -snes_rtol <rtol> - Sets rtol 2499 . -snes_stol <stol> - Sets stol 2500 . -snes_max_it <maxit> - Sets maxit 2501 - -snes_max_funcs <maxf> - Sets maxf 2502 2503 Notes: 2504 The default maximum number of iterations is 50. 2505 The default maximum number of function evaluations is 1000. 2506 2507 Level: intermediate 2508 2509 .keywords: SNES, nonlinear, set, convergence, tolerances 2510 2511 .seealso: SNESSetTrustRegionTolerance() 2512 @*/ 2513 PetscErrorCode SNESSetTolerances(SNES snes,PetscReal abstol,PetscReal rtol,PetscReal stol,PetscInt maxit,PetscInt maxf) 2514 { 2515 PetscFunctionBegin; 2516 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2517 PetscValidLogicalCollectiveReal(snes,abstol,2); 2518 PetscValidLogicalCollectiveReal(snes,rtol,3); 2519 PetscValidLogicalCollectiveReal(snes,stol,4); 2520 PetscValidLogicalCollectiveInt(snes,maxit,5); 2521 PetscValidLogicalCollectiveInt(snes,maxf,6); 2522 2523 if (abstol != PETSC_DEFAULT) { 2524 if (abstol < 0.0) SETERRQ1(((PetscObject)snes)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Absolute tolerance %G must be non-negative",abstol); 2525 snes->abstol = abstol; 2526 } 2527 if (rtol != PETSC_DEFAULT) { 2528 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); 2529 snes->rtol = rtol; 2530 } 2531 if (stol != PETSC_DEFAULT) { 2532 if (stol < 0.0) SETERRQ1(((PetscObject)snes)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Step tolerance %G must be non-negative",stol); 2533 snes->xtol = stol; 2534 } 2535 if (maxit != PETSC_DEFAULT) { 2536 if (maxit < 0) SETERRQ1(((PetscObject)snes)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Maximum number of iterations %D must be non-negative",maxit); 2537 snes->max_its = maxit; 2538 } 2539 if (maxf != PETSC_DEFAULT) { 2540 if (maxf < 0) SETERRQ1(((PetscObject)snes)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Maximum number of function evaluations %D must be non-negative",maxf); 2541 snes->max_funcs = maxf; 2542 } 2543 PetscFunctionReturn(0); 2544 } 2545 2546 #undef __FUNCT__ 2547 #define __FUNCT__ "SNESGetTolerances" 2548 /*@ 2549 SNESGetTolerances - Gets various parameters used in convergence tests. 2550 2551 Not Collective 2552 2553 Input Parameters: 2554 + snes - the SNES context 2555 . atol - absolute convergence tolerance 2556 . rtol - relative convergence tolerance 2557 . stol - convergence tolerance in terms of the norm 2558 of the change in the solution between steps 2559 . maxit - maximum number of iterations 2560 - maxf - maximum number of function evaluations 2561 2562 Notes: 2563 The user can specify PETSC_NULL for any parameter that is not needed. 2564 2565 Level: intermediate 2566 2567 .keywords: SNES, nonlinear, get, convergence, tolerances 2568 2569 .seealso: SNESSetTolerances() 2570 @*/ 2571 PetscErrorCode SNESGetTolerances(SNES snes,PetscReal *atol,PetscReal *rtol,PetscReal *stol,PetscInt *maxit,PetscInt *maxf) 2572 { 2573 PetscFunctionBegin; 2574 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2575 if (atol) *atol = snes->abstol; 2576 if (rtol) *rtol = snes->rtol; 2577 if (stol) *stol = snes->xtol; 2578 if (maxit) *maxit = snes->max_its; 2579 if (maxf) *maxf = snes->max_funcs; 2580 PetscFunctionReturn(0); 2581 } 2582 2583 #undef __FUNCT__ 2584 #define __FUNCT__ "SNESSetTrustRegionTolerance" 2585 /*@ 2586 SNESSetTrustRegionTolerance - Sets the trust region parameter tolerance. 2587 2588 Logically Collective on SNES 2589 2590 Input Parameters: 2591 + snes - the SNES context 2592 - tol - tolerance 2593 2594 Options Database Key: 2595 . -snes_trtol <tol> - Sets tol 2596 2597 Level: intermediate 2598 2599 .keywords: SNES, nonlinear, set, trust region, tolerance 2600 2601 .seealso: SNESSetTolerances() 2602 @*/ 2603 PetscErrorCode SNESSetTrustRegionTolerance(SNES snes,PetscReal tol) 2604 { 2605 PetscFunctionBegin; 2606 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2607 PetscValidLogicalCollectiveReal(snes,tol,2); 2608 snes->deltatol = tol; 2609 PetscFunctionReturn(0); 2610 } 2611 2612 /* 2613 Duplicate the lg monitors for SNES from KSP; for some reason with 2614 dynamic libraries things don't work under Sun4 if we just use 2615 macros instead of functions 2616 */ 2617 #undef __FUNCT__ 2618 #define __FUNCT__ "SNESMonitorLG" 2619 PetscErrorCode SNESMonitorLG(SNES snes,PetscInt it,PetscReal norm,void *ctx) 2620 { 2621 PetscErrorCode ierr; 2622 2623 PetscFunctionBegin; 2624 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2625 ierr = KSPMonitorLG((KSP)snes,it,norm,ctx);CHKERRQ(ierr); 2626 PetscFunctionReturn(0); 2627 } 2628 2629 #undef __FUNCT__ 2630 #define __FUNCT__ "SNESMonitorLGCreate" 2631 PetscErrorCode SNESMonitorLGCreate(const char host[],const char label[],int x,int y,int m,int n,PetscDrawLG *draw) 2632 { 2633 PetscErrorCode ierr; 2634 2635 PetscFunctionBegin; 2636 ierr = KSPMonitorLGCreate(host,label,x,y,m,n,draw);CHKERRQ(ierr); 2637 PetscFunctionReturn(0); 2638 } 2639 2640 #undef __FUNCT__ 2641 #define __FUNCT__ "SNESMonitorLGDestroy" 2642 PetscErrorCode SNESMonitorLGDestroy(PetscDrawLG *draw) 2643 { 2644 PetscErrorCode ierr; 2645 2646 PetscFunctionBegin; 2647 ierr = KSPMonitorLGDestroy(draw);CHKERRQ(ierr); 2648 PetscFunctionReturn(0); 2649 } 2650 2651 extern PetscErrorCode SNESMonitorRange_Private(SNES,PetscInt,PetscReal*); 2652 #undef __FUNCT__ 2653 #define __FUNCT__ "SNESMonitorLGRange" 2654 PetscErrorCode SNESMonitorLGRange(SNES snes,PetscInt n,PetscReal rnorm,void *monctx) 2655 { 2656 PetscDrawLG lg; 2657 PetscErrorCode ierr; 2658 PetscReal x,y,per; 2659 PetscViewer v = (PetscViewer)monctx; 2660 static PetscReal prev; /* should be in the context */ 2661 PetscDraw draw; 2662 PetscFunctionBegin; 2663 if (!monctx) { 2664 MPI_Comm comm; 2665 2666 ierr = PetscObjectGetComm((PetscObject)snes,&comm);CHKERRQ(ierr); 2667 v = PETSC_VIEWER_DRAW_(comm); 2668 } 2669 ierr = PetscViewerDrawGetDrawLG(v,0,&lg);CHKERRQ(ierr); 2670 if (!n) {ierr = PetscDrawLGReset(lg);CHKERRQ(ierr);} 2671 ierr = PetscDrawLGGetDraw(lg,&draw);CHKERRQ(ierr); 2672 ierr = PetscDrawSetTitle(draw,"Residual norm");CHKERRQ(ierr); 2673 x = (PetscReal) n; 2674 if (rnorm > 0.0) y = log10(rnorm); else y = -15.0; 2675 ierr = PetscDrawLGAddPoint(lg,&x,&y);CHKERRQ(ierr); 2676 if (n < 20 || !(n % 5)) { 2677 ierr = PetscDrawLGDraw(lg);CHKERRQ(ierr); 2678 } 2679 2680 ierr = PetscViewerDrawGetDrawLG(v,1,&lg);CHKERRQ(ierr); 2681 if (!n) {ierr = PetscDrawLGReset(lg);CHKERRQ(ierr);} 2682 ierr = PetscDrawLGGetDraw(lg,&draw);CHKERRQ(ierr); 2683 ierr = PetscDrawSetTitle(draw,"% elemts > .2*max elemt");CHKERRQ(ierr); 2684 ierr = SNESMonitorRange_Private(snes,n,&per);CHKERRQ(ierr); 2685 x = (PetscReal) n; 2686 y = 100.0*per; 2687 ierr = PetscDrawLGAddPoint(lg,&x,&y);CHKERRQ(ierr); 2688 if (n < 20 || !(n % 5)) { 2689 ierr = PetscDrawLGDraw(lg);CHKERRQ(ierr); 2690 } 2691 2692 ierr = PetscViewerDrawGetDrawLG(v,2,&lg);CHKERRQ(ierr); 2693 if (!n) {prev = rnorm;ierr = PetscDrawLGReset(lg);CHKERRQ(ierr);} 2694 ierr = PetscDrawLGGetDraw(lg,&draw);CHKERRQ(ierr); 2695 ierr = PetscDrawSetTitle(draw,"(norm -oldnorm)/oldnorm");CHKERRQ(ierr); 2696 x = (PetscReal) n; 2697 y = (prev - rnorm)/prev; 2698 ierr = PetscDrawLGAddPoint(lg,&x,&y);CHKERRQ(ierr); 2699 if (n < 20 || !(n % 5)) { 2700 ierr = PetscDrawLGDraw(lg);CHKERRQ(ierr); 2701 } 2702 2703 ierr = PetscViewerDrawGetDrawLG(v,3,&lg);CHKERRQ(ierr); 2704 if (!n) {ierr = PetscDrawLGReset(lg);CHKERRQ(ierr);} 2705 ierr = PetscDrawLGGetDraw(lg,&draw);CHKERRQ(ierr); 2706 ierr = PetscDrawSetTitle(draw,"(norm -oldnorm)/oldnorm*(% > .2 max)");CHKERRQ(ierr); 2707 x = (PetscReal) n; 2708 y = (prev - rnorm)/(prev*per); 2709 if (n > 2) { /*skip initial crazy value */ 2710 ierr = PetscDrawLGAddPoint(lg,&x,&y);CHKERRQ(ierr); 2711 } 2712 if (n < 20 || !(n % 5)) { 2713 ierr = PetscDrawLGDraw(lg);CHKERRQ(ierr); 2714 } 2715 prev = rnorm; 2716 PetscFunctionReturn(0); 2717 } 2718 2719 #undef __FUNCT__ 2720 #define __FUNCT__ "SNESMonitorLGRangeCreate" 2721 PetscErrorCode SNESMonitorLGRangeCreate(const char host[],const char label[],int x,int y,int m,int n,PetscDrawLG *draw) 2722 { 2723 PetscErrorCode ierr; 2724 2725 PetscFunctionBegin; 2726 ierr = KSPMonitorLGCreate(host,label,x,y,m,n,draw);CHKERRQ(ierr); 2727 PetscFunctionReturn(0); 2728 } 2729 2730 #undef __FUNCT__ 2731 #define __FUNCT__ "SNESMonitorLGRangeDestroy" 2732 PetscErrorCode SNESMonitorLGRangeDestroy(PetscDrawLG *draw) 2733 { 2734 PetscErrorCode ierr; 2735 2736 PetscFunctionBegin; 2737 ierr = KSPMonitorLGDestroy(draw);CHKERRQ(ierr); 2738 PetscFunctionReturn(0); 2739 } 2740 2741 #undef __FUNCT__ 2742 #define __FUNCT__ "SNESMonitor" 2743 /*@ 2744 SNESMonitor - runs the user provided monitor routines, if they exist 2745 2746 Collective on SNES 2747 2748 Input Parameters: 2749 + snes - nonlinear solver context obtained from SNESCreate() 2750 . iter - iteration number 2751 - rnorm - relative norm of the residual 2752 2753 Notes: 2754 This routine is called by the SNES implementations. 2755 It does not typically need to be called by the user. 2756 2757 Level: developer 2758 2759 .seealso: SNESMonitorSet() 2760 @*/ 2761 PetscErrorCode SNESMonitor(SNES snes,PetscInt iter,PetscReal rnorm) 2762 { 2763 PetscErrorCode ierr; 2764 PetscInt i,n = snes->numbermonitors; 2765 2766 PetscFunctionBegin; 2767 for (i=0; i<n; i++) { 2768 ierr = (*snes->monitor[i])(snes,iter,rnorm,snes->monitorcontext[i]);CHKERRQ(ierr); 2769 } 2770 PetscFunctionReturn(0); 2771 } 2772 2773 /* ------------ Routines to set performance monitoring options ----------- */ 2774 2775 #undef __FUNCT__ 2776 #define __FUNCT__ "SNESMonitorSet" 2777 /*@C 2778 SNESMonitorSet - Sets an ADDITIONAL function that is to be used at every 2779 iteration of the nonlinear solver to display the iteration's 2780 progress. 2781 2782 Logically Collective on SNES 2783 2784 Input Parameters: 2785 + snes - the SNES context 2786 . func - monitoring routine 2787 . mctx - [optional] user-defined context for private data for the 2788 monitor routine (use PETSC_NULL if no context is desired) 2789 - monitordestroy - [optional] routine that frees monitor context 2790 (may be PETSC_NULL) 2791 2792 Calling sequence of func: 2793 $ int func(SNES snes,PetscInt its, PetscReal norm,void *mctx) 2794 2795 + snes - the SNES context 2796 . its - iteration number 2797 . norm - 2-norm function value (may be estimated) 2798 - mctx - [optional] monitoring context 2799 2800 Options Database Keys: 2801 + -snes_monitor - sets SNESMonitorDefault() 2802 . -snes_monitor_draw - sets line graph monitor, 2803 uses SNESMonitorLGCreate() 2804 - -snes_monitor_cancel - cancels all monitors that have 2805 been hardwired into a code by 2806 calls to SNESMonitorSet(), but 2807 does not cancel those set via 2808 the options database. 2809 2810 Notes: 2811 Several different monitoring routines may be set by calling 2812 SNESMonitorSet() multiple times; all will be called in the 2813 order in which they were set. 2814 2815 Fortran notes: Only a single monitor function can be set for each SNES object 2816 2817 Level: intermediate 2818 2819 .keywords: SNES, nonlinear, set, monitor 2820 2821 .seealso: SNESMonitorDefault(), SNESMonitorCancel() 2822 @*/ 2823 PetscErrorCode SNESMonitorSet(SNES snes,PetscErrorCode (*monitor)(SNES,PetscInt,PetscReal,void*),void *mctx,PetscErrorCode (*monitordestroy)(void**)) 2824 { 2825 PetscInt i; 2826 PetscErrorCode ierr; 2827 2828 PetscFunctionBegin; 2829 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2830 if (snes->numbermonitors >= MAXSNESMONITORS) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Too many monitors set"); 2831 for (i=0; i<snes->numbermonitors;i++) { 2832 if (monitor == snes->monitor[i] && monitordestroy == snes->monitordestroy[i] && mctx == snes->monitorcontext[i]) { 2833 if (monitordestroy) { 2834 ierr = (*monitordestroy)(&mctx);CHKERRQ(ierr); 2835 } 2836 PetscFunctionReturn(0); 2837 } 2838 } 2839 snes->monitor[snes->numbermonitors] = monitor; 2840 snes->monitordestroy[snes->numbermonitors] = monitordestroy; 2841 snes->monitorcontext[snes->numbermonitors++] = (void*)mctx; 2842 PetscFunctionReturn(0); 2843 } 2844 2845 #undef __FUNCT__ 2846 #define __FUNCT__ "SNESMonitorCancel" 2847 /*@C 2848 SNESMonitorCancel - Clears all the monitor functions for a SNES object. 2849 2850 Logically Collective on SNES 2851 2852 Input Parameters: 2853 . snes - the SNES context 2854 2855 Options Database Key: 2856 . -snes_monitor_cancel - cancels all monitors that have been hardwired 2857 into a code by calls to SNESMonitorSet(), but does not cancel those 2858 set via the options database 2859 2860 Notes: 2861 There is no way to clear one specific monitor from a SNES object. 2862 2863 Level: intermediate 2864 2865 .keywords: SNES, nonlinear, set, monitor 2866 2867 .seealso: SNESMonitorDefault(), SNESMonitorSet() 2868 @*/ 2869 PetscErrorCode SNESMonitorCancel(SNES snes) 2870 { 2871 PetscErrorCode ierr; 2872 PetscInt i; 2873 2874 PetscFunctionBegin; 2875 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2876 for (i=0; i<snes->numbermonitors; i++) { 2877 if (snes->monitordestroy[i]) { 2878 ierr = (*snes->monitordestroy[i])(&snes->monitorcontext[i]);CHKERRQ(ierr); 2879 } 2880 } 2881 snes->numbermonitors = 0; 2882 PetscFunctionReturn(0); 2883 } 2884 2885 #undef __FUNCT__ 2886 #define __FUNCT__ "SNESSetConvergenceTest" 2887 /*@C 2888 SNESSetConvergenceTest - Sets the function that is to be used 2889 to test for convergence of the nonlinear iterative solution. 2890 2891 Logically Collective on SNES 2892 2893 Input Parameters: 2894 + snes - the SNES context 2895 . func - routine to test for convergence 2896 . cctx - [optional] context for private data for the convergence routine (may be PETSC_NULL) 2897 - destroy - [optional] destructor for the context (may be PETSC_NULL; PETSC_NULL_FUNCTION in Fortran) 2898 2899 Calling sequence of func: 2900 $ PetscErrorCode func (SNES snes,PetscInt it,PetscReal xnorm,PetscReal gnorm,PetscReal f,SNESConvergedReason *reason,void *cctx) 2901 2902 + snes - the SNES context 2903 . it - current iteration (0 is the first and is before any Newton step) 2904 . cctx - [optional] convergence context 2905 . reason - reason for convergence/divergence 2906 . xnorm - 2-norm of current iterate 2907 . gnorm - 2-norm of current step 2908 - f - 2-norm of function 2909 2910 Level: advanced 2911 2912 .keywords: SNES, nonlinear, set, convergence, test 2913 2914 .seealso: SNESDefaultConverged(), SNESSkipConverged() 2915 @*/ 2916 PetscErrorCode SNESSetConvergenceTest(SNES snes,PetscErrorCode (*func)(SNES,PetscInt,PetscReal,PetscReal,PetscReal,SNESConvergedReason*,void*),void *cctx,PetscErrorCode (*destroy)(void*)) 2917 { 2918 PetscErrorCode ierr; 2919 2920 PetscFunctionBegin; 2921 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2922 if (!func) func = SNESSkipConverged; 2923 if (snes->ops->convergeddestroy) { 2924 ierr = (*snes->ops->convergeddestroy)(snes->cnvP);CHKERRQ(ierr); 2925 } 2926 snes->ops->converged = func; 2927 snes->ops->convergeddestroy = destroy; 2928 snes->cnvP = cctx; 2929 PetscFunctionReturn(0); 2930 } 2931 2932 #undef __FUNCT__ 2933 #define __FUNCT__ "SNESGetConvergedReason" 2934 /*@ 2935 SNESGetConvergedReason - Gets the reason the SNES iteration was stopped. 2936 2937 Not Collective 2938 2939 Input Parameter: 2940 . snes - the SNES context 2941 2942 Output Parameter: 2943 . reason - negative value indicates diverged, positive value converged, see SNESConvergedReason or the 2944 manual pages for the individual convergence tests for complete lists 2945 2946 Level: intermediate 2947 2948 Notes: Can only be called after the call the SNESSolve() is complete. 2949 2950 .keywords: SNES, nonlinear, set, convergence, test 2951 2952 .seealso: SNESSetConvergenceTest(), SNESConvergedReason 2953 @*/ 2954 PetscErrorCode SNESGetConvergedReason(SNES snes,SNESConvergedReason *reason) 2955 { 2956 PetscFunctionBegin; 2957 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2958 PetscValidPointer(reason,2); 2959 *reason = snes->reason; 2960 PetscFunctionReturn(0); 2961 } 2962 2963 #undef __FUNCT__ 2964 #define __FUNCT__ "SNESSetConvergenceHistory" 2965 /*@ 2966 SNESSetConvergenceHistory - Sets the array used to hold the convergence history. 2967 2968 Logically Collective on SNES 2969 2970 Input Parameters: 2971 + snes - iterative context obtained from SNESCreate() 2972 . a - array to hold history, this array will contain the function norms computed at each step 2973 . its - integer array holds the number of linear iterations for each solve. 2974 . na - size of a and its 2975 - reset - PETSC_TRUE indicates each new nonlinear solve resets the history counter to zero, 2976 else it continues storing new values for new nonlinear solves after the old ones 2977 2978 Notes: 2979 If 'a' and 'its' are PETSC_NULL then space is allocated for the history. If 'na' PETSC_DECIDE or PETSC_DEFAULT then a 2980 default array of length 10000 is allocated. 2981 2982 This routine is useful, e.g., when running a code for purposes 2983 of accurate performance monitoring, when no I/O should be done 2984 during the section of code that is being timed. 2985 2986 Level: intermediate 2987 2988 .keywords: SNES, set, convergence, history 2989 2990 .seealso: SNESGetConvergenceHistory() 2991 2992 @*/ 2993 PetscErrorCode SNESSetConvergenceHistory(SNES snes,PetscReal a[],PetscInt its[],PetscInt na,PetscBool reset) 2994 { 2995 PetscErrorCode ierr; 2996 2997 PetscFunctionBegin; 2998 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2999 if (na) PetscValidScalarPointer(a,2); 3000 if (its) PetscValidIntPointer(its,3); 3001 if (na == PETSC_DECIDE || na == PETSC_DEFAULT || !a) { 3002 if (na == PETSC_DECIDE || na == PETSC_DEFAULT) na = 1000; 3003 ierr = PetscMalloc(na*sizeof(PetscReal),&a);CHKERRQ(ierr); 3004 ierr = PetscMalloc(na*sizeof(PetscInt),&its);CHKERRQ(ierr); 3005 snes->conv_malloc = PETSC_TRUE; 3006 } 3007 snes->conv_hist = a; 3008 snes->conv_hist_its = its; 3009 snes->conv_hist_max = na; 3010 snes->conv_hist_len = 0; 3011 snes->conv_hist_reset = reset; 3012 PetscFunctionReturn(0); 3013 } 3014 3015 #if defined(PETSC_HAVE_MATLAB_ENGINE) 3016 #include <engine.h> /* MATLAB include file */ 3017 #include <mex.h> /* MATLAB include file */ 3018 EXTERN_C_BEGIN 3019 #undef __FUNCT__ 3020 #define __FUNCT__ "SNESGetConvergenceHistoryMatlab" 3021 mxArray *SNESGetConvergenceHistoryMatlab(SNES snes) 3022 { 3023 mxArray *mat; 3024 PetscInt i; 3025 PetscReal *ar; 3026 3027 PetscFunctionBegin; 3028 mat = mxCreateDoubleMatrix(snes->conv_hist_len,1,mxREAL); 3029 ar = (PetscReal*) mxGetData(mat); 3030 for (i=0; i<snes->conv_hist_len; i++) { 3031 ar[i] = snes->conv_hist[i]; 3032 } 3033 PetscFunctionReturn(mat); 3034 } 3035 EXTERN_C_END 3036 #endif 3037 3038 3039 #undef __FUNCT__ 3040 #define __FUNCT__ "SNESGetConvergenceHistory" 3041 /*@C 3042 SNESGetConvergenceHistory - Gets the array used to hold the convergence history. 3043 3044 Not Collective 3045 3046 Input Parameter: 3047 . snes - iterative context obtained from SNESCreate() 3048 3049 Output Parameters: 3050 . a - array to hold history 3051 . its - integer array holds the number of linear iterations (or 3052 negative if not converged) for each solve. 3053 - na - size of a and its 3054 3055 Notes: 3056 The calling sequence for this routine in Fortran is 3057 $ call SNESGetConvergenceHistory(SNES snes, integer na, integer ierr) 3058 3059 This routine is useful, e.g., when running a code for purposes 3060 of accurate performance monitoring, when no I/O should be done 3061 during the section of code that is being timed. 3062 3063 Level: intermediate 3064 3065 .keywords: SNES, get, convergence, history 3066 3067 .seealso: SNESSetConvergencHistory() 3068 3069 @*/ 3070 PetscErrorCode SNESGetConvergenceHistory(SNES snes,PetscReal *a[],PetscInt *its[],PetscInt *na) 3071 { 3072 PetscFunctionBegin; 3073 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3074 if (a) *a = snes->conv_hist; 3075 if (its) *its = snes->conv_hist_its; 3076 if (na) *na = snes->conv_hist_len; 3077 PetscFunctionReturn(0); 3078 } 3079 3080 #undef __FUNCT__ 3081 #define __FUNCT__ "SNESSetUpdate" 3082 /*@C 3083 SNESSetUpdate - Sets the general-purpose update function called 3084 at the beginning of every iteration of the nonlinear solve. Specifically 3085 it is called just before the Jacobian is "evaluated". 3086 3087 Logically Collective on SNES 3088 3089 Input Parameters: 3090 . snes - The nonlinear solver context 3091 . func - The function 3092 3093 Calling sequence of func: 3094 . func (SNES snes, PetscInt step); 3095 3096 . step - The current step of the iteration 3097 3098 Level: advanced 3099 3100 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() 3101 This is not used by most users. 3102 3103 .keywords: SNES, update 3104 3105 .seealso SNESDefaultUpdate(), SNESSetJacobian(), SNESSolve() 3106 @*/ 3107 PetscErrorCode SNESSetUpdate(SNES snes, PetscErrorCode (*func)(SNES, PetscInt)) 3108 { 3109 PetscFunctionBegin; 3110 PetscValidHeaderSpecific(snes, SNES_CLASSID,1); 3111 snes->ops->update = func; 3112 PetscFunctionReturn(0); 3113 } 3114 3115 #undef __FUNCT__ 3116 #define __FUNCT__ "SNESDefaultUpdate" 3117 /*@ 3118 SNESDefaultUpdate - The default update function which does nothing. 3119 3120 Not collective 3121 3122 Input Parameters: 3123 . snes - The nonlinear solver context 3124 . step - The current step of the iteration 3125 3126 Level: intermediate 3127 3128 .keywords: SNES, update 3129 .seealso SNESSetUpdate(), SNESDefaultRhsBC(), SNESDefaultShortolutionBC() 3130 @*/ 3131 PetscErrorCode SNESDefaultUpdate(SNES snes, PetscInt step) 3132 { 3133 PetscFunctionBegin; 3134 PetscFunctionReturn(0); 3135 } 3136 3137 #undef __FUNCT__ 3138 #define __FUNCT__ "SNESScaleStep_Private" 3139 /* 3140 SNESScaleStep_Private - Scales a step so that its length is less than the 3141 positive parameter delta. 3142 3143 Input Parameters: 3144 + snes - the SNES context 3145 . y - approximate solution of linear system 3146 . fnorm - 2-norm of current function 3147 - delta - trust region size 3148 3149 Output Parameters: 3150 + gpnorm - predicted function norm at the new point, assuming local 3151 linearization. The value is zero if the step lies within the trust 3152 region, and exceeds zero otherwise. 3153 - ynorm - 2-norm of the step 3154 3155 Note: 3156 For non-trust region methods such as SNESLS, the parameter delta 3157 is set to be the maximum allowable step size. 3158 3159 .keywords: SNES, nonlinear, scale, step 3160 */ 3161 PetscErrorCode SNESScaleStep_Private(SNES snes,Vec y,PetscReal *fnorm,PetscReal *delta,PetscReal *gpnorm,PetscReal *ynorm) 3162 { 3163 PetscReal nrm; 3164 PetscScalar cnorm; 3165 PetscErrorCode ierr; 3166 3167 PetscFunctionBegin; 3168 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3169 PetscValidHeaderSpecific(y,VEC_CLASSID,2); 3170 PetscCheckSameComm(snes,1,y,2); 3171 3172 ierr = VecNorm(y,NORM_2,&nrm);CHKERRQ(ierr); 3173 if (nrm > *delta) { 3174 nrm = *delta/nrm; 3175 *gpnorm = (1.0 - nrm)*(*fnorm); 3176 cnorm = nrm; 3177 ierr = VecScale(y,cnorm);CHKERRQ(ierr); 3178 *ynorm = *delta; 3179 } else { 3180 *gpnorm = 0.0; 3181 *ynorm = nrm; 3182 } 3183 PetscFunctionReturn(0); 3184 } 3185 3186 #undef __FUNCT__ 3187 #define __FUNCT__ "SNESSolve" 3188 /*@C 3189 SNESSolve - Solves a nonlinear system F(x) = b. 3190 Call SNESSolve() after calling SNESCreate() and optional routines of the form SNESSetXXX(). 3191 3192 Collective on SNES 3193 3194 Input Parameters: 3195 + snes - the SNES context 3196 . b - the constant part of the equation F(x) = b, or PETSC_NULL to use zero. 3197 - x - the solution vector. 3198 3199 Notes: 3200 The user should initialize the vector,x, with the initial guess 3201 for the nonlinear solve prior to calling SNESSolve. In particular, 3202 to employ an initial guess of zero, the user should explicitly set 3203 this vector to zero by calling VecSet(). 3204 3205 Level: beginner 3206 3207 .keywords: SNES, nonlinear, solve 3208 3209 .seealso: SNESCreate(), SNESDestroy(), SNESSetFunction(), SNESSetJacobian(), SNESSetGridSequence(), SNESGetSolution() 3210 @*/ 3211 PetscErrorCode SNESSolve(SNES snes,Vec b,Vec x) 3212 { 3213 PetscErrorCode ierr; 3214 PetscBool flg; 3215 char filename[PETSC_MAX_PATH_LEN]; 3216 PetscViewer viewer; 3217 PetscInt grid; 3218 Vec xcreated = PETSC_NULL; 3219 3220 PetscFunctionBegin; 3221 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3222 if (x) PetscValidHeaderSpecific(x,VEC_CLASSID,3); 3223 if (x) PetscCheckSameComm(snes,1,x,3); 3224 if (b) PetscValidHeaderSpecific(b,VEC_CLASSID,2); 3225 if (b) PetscCheckSameComm(snes,1,b,2); 3226 3227 if (!x && snes->dm) { 3228 ierr = DMCreateGlobalVector(snes->dm,&xcreated);CHKERRQ(ierr); 3229 x = xcreated; 3230 } 3231 3232 for (grid=0; grid<snes->gridsequence; grid++) {ierr = PetscViewerASCIIPushTab(PETSC_VIEWER_STDOUT_(((PetscObject)snes)->comm));CHKERRQ(ierr);} 3233 for (grid=0; grid<snes->gridsequence+1; grid++) { 3234 3235 /* set solution vector */ 3236 if (!grid) {ierr = PetscObjectReference((PetscObject)x);CHKERRQ(ierr);} 3237 ierr = VecDestroy(&snes->vec_sol);CHKERRQ(ierr); 3238 snes->vec_sol = x; 3239 /* set afine vector if provided */ 3240 if (b) { ierr = PetscObjectReference((PetscObject)b);CHKERRQ(ierr); } 3241 ierr = VecDestroy(&snes->vec_rhs);CHKERRQ(ierr); 3242 snes->vec_rhs = b; 3243 3244 ierr = SNESSetUp(snes);CHKERRQ(ierr); 3245 3246 if (!grid) { 3247 if (snes->ops->computeinitialguess) { 3248 ierr = (*snes->ops->computeinitialguess)(snes,snes->vec_sol,snes->initialguessP);CHKERRQ(ierr); 3249 } else if (snes->dm) { 3250 PetscBool ig; 3251 ierr = DMHasInitialGuess(snes->dm,&ig);CHKERRQ(ierr); 3252 if (ig) { 3253 ierr = DMComputeInitialGuess(snes->dm,snes->vec_sol);CHKERRQ(ierr); 3254 } 3255 } 3256 } 3257 3258 if (snes->conv_hist_reset) snes->conv_hist_len = 0; 3259 snes->nfuncs = 0; snes->linear_its = 0; snes->numFailures = 0; 3260 3261 ierr = PetscLogEventBegin(SNES_Solve,snes,0,0,0);CHKERRQ(ierr); 3262 ierr = (*snes->ops->solve)(snes);CHKERRQ(ierr); 3263 ierr = PetscLogEventEnd(SNES_Solve,snes,0,0,0);CHKERRQ(ierr); 3264 if (snes->domainerror){ 3265 snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; 3266 snes->domainerror = PETSC_FALSE; 3267 } 3268 if (!snes->reason) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Internal error, solver returned without setting converged reason"); 3269 3270 ierr = PetscOptionsGetString(((PetscObject)snes)->prefix,"-snes_view",filename,PETSC_MAX_PATH_LEN,&flg);CHKERRQ(ierr); 3271 if (flg && !PetscPreLoadingOn) { 3272 ierr = PetscViewerASCIIOpen(((PetscObject)snes)->comm,filename,&viewer);CHKERRQ(ierr); 3273 ierr = SNESView(snes,viewer);CHKERRQ(ierr); 3274 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 3275 } 3276 3277 flg = PETSC_FALSE; 3278 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_test_local_min",&flg,PETSC_NULL);CHKERRQ(ierr); 3279 if (flg && !PetscPreLoadingOn) { ierr = SNESTestLocalMin(snes);CHKERRQ(ierr); } 3280 if (snes->printreason) { 3281 ierr = PetscViewerASCIIAddTab(PETSC_VIEWER_STDOUT_(((PetscObject)snes)->comm),((PetscObject)snes)->tablevel);CHKERRQ(ierr); 3282 if (snes->reason > 0) { 3283 ierr = PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_(((PetscObject)snes)->comm),"Nonlinear solve converged due to %s\n",SNESConvergedReasons[snes->reason]);CHKERRQ(ierr); 3284 } else { 3285 ierr = PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_(((PetscObject)snes)->comm),"Nonlinear solve did not converge due to %s\n",SNESConvergedReasons[snes->reason]);CHKERRQ(ierr); 3286 } 3287 ierr = PetscViewerASCIISubtractTab(PETSC_VIEWER_STDOUT_(((PetscObject)snes)->comm),((PetscObject)snes)->tablevel);CHKERRQ(ierr); 3288 } 3289 3290 flg = PETSC_FALSE; 3291 ierr = PetscOptionsGetString(((PetscObject)snes)->prefix,"-snes_view_solution_vtk",filename,PETSC_MAX_PATH_LEN,&flg);CHKERRQ(ierr); 3292 if (flg) { 3293 PetscViewer viewer; 3294 ierr = PetscViewerCreate(((PetscObject)snes)->comm,&viewer);CHKERRQ(ierr); 3295 ierr = PetscViewerSetType(viewer,PETSCVIEWERVTK);CHKERRQ(ierr); 3296 ierr = PetscViewerFileSetName(viewer,filename);CHKERRQ(ierr); 3297 ierr = VecView(snes->vec_sol,viewer);CHKERRQ(ierr); 3298 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 3299 } 3300 3301 if (snes->errorifnotconverged && snes->reason < 0) SETERRQ(((PetscObject)snes)->comm,PETSC_ERR_NOT_CONVERGED,"SNESSolve has not converged"); 3302 if (grid < snes->gridsequence) { 3303 DM fine; 3304 Vec xnew; 3305 Mat interp; 3306 3307 ierr = DMRefine(snes->dm,((PetscObject)snes)->comm,&fine);CHKERRQ(ierr); 3308 ierr = DMCreateInterpolation(snes->dm,fine,&interp,PETSC_NULL);CHKERRQ(ierr); 3309 ierr = DMCreateGlobalVector(fine,&xnew);CHKERRQ(ierr); 3310 ierr = MatInterpolate(interp,x,xnew);CHKERRQ(ierr); 3311 ierr = MatDestroy(&interp);CHKERRQ(ierr); 3312 x = xnew; 3313 3314 ierr = SNESReset(snes);CHKERRQ(ierr); 3315 ierr = SNESSetDM(snes,fine);CHKERRQ(ierr); 3316 ierr = DMDestroy(&fine);CHKERRQ(ierr); 3317 ierr = PetscViewerASCIIPopTab(PETSC_VIEWER_STDOUT_(((PetscObject)snes)->comm));CHKERRQ(ierr); 3318 } 3319 } 3320 ierr = VecDestroy(&xcreated);CHKERRQ(ierr); 3321 PetscFunctionReturn(0); 3322 } 3323 3324 /* --------- Internal routines for SNES Package --------- */ 3325 3326 #undef __FUNCT__ 3327 #define __FUNCT__ "SNESSetType" 3328 /*@C 3329 SNESSetType - Sets the method for the nonlinear solver. 3330 3331 Collective on SNES 3332 3333 Input Parameters: 3334 + snes - the SNES context 3335 - type - a known method 3336 3337 Options Database Key: 3338 . -snes_type <type> - Sets the method; use -help for a list 3339 of available methods (for instance, ls or tr) 3340 3341 Notes: 3342 See "petsc/include/petscsnes.h" for available methods (for instance) 3343 + SNESLS - Newton's method with line search 3344 (systems of nonlinear equations) 3345 . SNESTR - Newton's method with trust region 3346 (systems of nonlinear equations) 3347 3348 Normally, it is best to use the SNESSetFromOptions() command and then 3349 set the SNES solver type from the options database rather than by using 3350 this routine. Using the options database provides the user with 3351 maximum flexibility in evaluating the many nonlinear solvers. 3352 The SNESSetType() routine is provided for those situations where it 3353 is necessary to set the nonlinear solver independently of the command 3354 line or options database. This might be the case, for example, when 3355 the choice of solver changes during the execution of the program, 3356 and the user's application is taking responsibility for choosing the 3357 appropriate method. 3358 3359 Level: intermediate 3360 3361 .keywords: SNES, set, type 3362 3363 .seealso: SNESType, SNESCreate() 3364 3365 @*/ 3366 PetscErrorCode SNESSetType(SNES snes,const SNESType type) 3367 { 3368 PetscErrorCode ierr,(*r)(SNES); 3369 PetscBool match; 3370 3371 PetscFunctionBegin; 3372 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3373 PetscValidCharPointer(type,2); 3374 3375 ierr = PetscTypeCompare((PetscObject)snes,type,&match);CHKERRQ(ierr); 3376 if (match) PetscFunctionReturn(0); 3377 3378 ierr = PetscFListFind(SNESList,((PetscObject)snes)->comm,type,PETSC_TRUE,(void (**)(void)) &r);CHKERRQ(ierr); 3379 if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unable to find requested SNES type %s",type); 3380 /* Destroy the previous private SNES context */ 3381 if (snes->ops->destroy) { 3382 ierr = (*(snes)->ops->destroy)(snes);CHKERRQ(ierr); 3383 snes->ops->destroy = PETSC_NULL; 3384 } 3385 /* Reinitialize function pointers in SNESOps structure */ 3386 snes->ops->setup = 0; 3387 snes->ops->solve = 0; 3388 snes->ops->view = 0; 3389 snes->ops->setfromoptions = 0; 3390 snes->ops->destroy = 0; 3391 /* Call the SNESCreate_XXX routine for this particular Nonlinear solver */ 3392 snes->setupcalled = PETSC_FALSE; 3393 ierr = PetscObjectChangeTypeName((PetscObject)snes,type);CHKERRQ(ierr); 3394 ierr = (*r)(snes);CHKERRQ(ierr); 3395 #if defined(PETSC_HAVE_AMS) 3396 if (PetscAMSPublishAll) { 3397 ierr = PetscObjectAMSPublish((PetscObject)snes);CHKERRQ(ierr); 3398 } 3399 #endif 3400 PetscFunctionReturn(0); 3401 } 3402 3403 3404 /* --------------------------------------------------------------------- */ 3405 #undef __FUNCT__ 3406 #define __FUNCT__ "SNESRegisterDestroy" 3407 /*@ 3408 SNESRegisterDestroy - Frees the list of nonlinear solvers that were 3409 registered by SNESRegisterDynamic(). 3410 3411 Not Collective 3412 3413 Level: advanced 3414 3415 .keywords: SNES, nonlinear, register, destroy 3416 3417 .seealso: SNESRegisterAll(), SNESRegisterAll() 3418 @*/ 3419 PetscErrorCode SNESRegisterDestroy(void) 3420 { 3421 PetscErrorCode ierr; 3422 3423 PetscFunctionBegin; 3424 ierr = PetscFListDestroy(&SNESList);CHKERRQ(ierr); 3425 SNESRegisterAllCalled = PETSC_FALSE; 3426 PetscFunctionReturn(0); 3427 } 3428 3429 #undef __FUNCT__ 3430 #define __FUNCT__ "SNESGetType" 3431 /*@C 3432 SNESGetType - Gets the SNES method type and name (as a string). 3433 3434 Not Collective 3435 3436 Input Parameter: 3437 . snes - nonlinear solver context 3438 3439 Output Parameter: 3440 . type - SNES method (a character string) 3441 3442 Level: intermediate 3443 3444 .keywords: SNES, nonlinear, get, type, name 3445 @*/ 3446 PetscErrorCode SNESGetType(SNES snes,const SNESType *type) 3447 { 3448 PetscFunctionBegin; 3449 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3450 PetscValidPointer(type,2); 3451 *type = ((PetscObject)snes)->type_name; 3452 PetscFunctionReturn(0); 3453 } 3454 3455 #undef __FUNCT__ 3456 #define __FUNCT__ "SNESGetSolution" 3457 /*@ 3458 SNESGetSolution - Returns the vector where the approximate solution is 3459 stored. This is the fine grid solution when using SNESSetGridSequence(). 3460 3461 Not Collective, but Vec is parallel if SNES is parallel 3462 3463 Input Parameter: 3464 . snes - the SNES context 3465 3466 Output Parameter: 3467 . x - the solution 3468 3469 Level: intermediate 3470 3471 .keywords: SNES, nonlinear, get, solution 3472 3473 .seealso: SNESGetSolutionUpdate(), SNESGetFunction() 3474 @*/ 3475 PetscErrorCode SNESGetSolution(SNES snes,Vec *x) 3476 { 3477 PetscFunctionBegin; 3478 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3479 PetscValidPointer(x,2); 3480 *x = snes->vec_sol; 3481 PetscFunctionReturn(0); 3482 } 3483 3484 #undef __FUNCT__ 3485 #define __FUNCT__ "SNESGetSolutionUpdate" 3486 /*@ 3487 SNESGetSolutionUpdate - Returns the vector where the solution update is 3488 stored. 3489 3490 Not Collective, but Vec is parallel if SNES is parallel 3491 3492 Input Parameter: 3493 . snes - the SNES context 3494 3495 Output Parameter: 3496 . x - the solution update 3497 3498 Level: advanced 3499 3500 .keywords: SNES, nonlinear, get, solution, update 3501 3502 .seealso: SNESGetSolution(), SNESGetFunction() 3503 @*/ 3504 PetscErrorCode SNESGetSolutionUpdate(SNES snes,Vec *x) 3505 { 3506 PetscFunctionBegin; 3507 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3508 PetscValidPointer(x,2); 3509 *x = snes->vec_sol_update; 3510 PetscFunctionReturn(0); 3511 } 3512 3513 #undef __FUNCT__ 3514 #define __FUNCT__ "SNESGetFunction" 3515 /*@C 3516 SNESGetFunction - Returns the vector where the function is stored. 3517 3518 Not Collective, but Vec is parallel if SNES is parallel. Collective if Vec is requested, but has not been created yet. 3519 3520 Input Parameter: 3521 . snes - the SNES context 3522 3523 Output Parameter: 3524 + r - the function (or PETSC_NULL) 3525 . func - the function (or PETSC_NULL) 3526 - ctx - the function context (or PETSC_NULL) 3527 3528 Level: advanced 3529 3530 .keywords: SNES, nonlinear, get, function 3531 3532 .seealso: SNESSetFunction(), SNESGetSolution() 3533 @*/ 3534 PetscErrorCode SNESGetFunction(SNES snes,Vec *r,PetscErrorCode (**func)(SNES,Vec,Vec,void*),void **ctx) 3535 { 3536 PetscErrorCode ierr; 3537 DM dm; 3538 3539 PetscFunctionBegin; 3540 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3541 if (r) { 3542 if (!snes->vec_func) { 3543 if (snes->vec_rhs) { 3544 ierr = VecDuplicate(snes->vec_rhs,&snes->vec_func);CHKERRQ(ierr); 3545 } else if (snes->vec_sol) { 3546 ierr = VecDuplicate(snes->vec_sol,&snes->vec_func);CHKERRQ(ierr); 3547 } else if (snes->dm) { 3548 ierr = DMCreateGlobalVector(snes->dm,&snes->vec_func);CHKERRQ(ierr); 3549 } 3550 } 3551 *r = snes->vec_func; 3552 } 3553 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 3554 ierr = DMSNESGetFunction(dm,func,ctx);CHKERRQ(ierr); 3555 PetscFunctionReturn(0); 3556 } 3557 3558 /*@C 3559 SNESGetGS - Returns the GS function and context. 3560 3561 Input Parameter: 3562 . snes - the SNES context 3563 3564 Output Parameter: 3565 + gsfunc - the function (or PETSC_NULL) 3566 - ctx - the function context (or PETSC_NULL) 3567 3568 Level: advanced 3569 3570 .keywords: SNES, nonlinear, get, function 3571 3572 .seealso: SNESSetGS(), SNESGetFunction() 3573 @*/ 3574 3575 #undef __FUNCT__ 3576 #define __FUNCT__ "SNESGetGS" 3577 PetscErrorCode SNESGetGS (SNES snes, PetscErrorCode(**func)(SNES, Vec, Vec, void*), void ** ctx) 3578 { 3579 PetscErrorCode ierr; 3580 DM dm; 3581 3582 PetscFunctionBegin; 3583 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3584 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 3585 ierr = DMSNESGetGS(dm,func,ctx);CHKERRQ(ierr); 3586 PetscFunctionReturn(0); 3587 } 3588 3589 #undef __FUNCT__ 3590 #define __FUNCT__ "SNESSetOptionsPrefix" 3591 /*@C 3592 SNESSetOptionsPrefix - Sets the prefix used for searching for all 3593 SNES options in the database. 3594 3595 Logically Collective on SNES 3596 3597 Input Parameter: 3598 + snes - the SNES context 3599 - prefix - the prefix to prepend to all option names 3600 3601 Notes: 3602 A hyphen (-) must NOT be given at the beginning of the prefix name. 3603 The first character of all runtime options is AUTOMATICALLY the hyphen. 3604 3605 Level: advanced 3606 3607 .keywords: SNES, set, options, prefix, database 3608 3609 .seealso: SNESSetFromOptions() 3610 @*/ 3611 PetscErrorCode SNESSetOptionsPrefix(SNES snes,const char prefix[]) 3612 { 3613 PetscErrorCode ierr; 3614 3615 PetscFunctionBegin; 3616 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3617 ierr = PetscObjectSetOptionsPrefix((PetscObject)snes,prefix);CHKERRQ(ierr); 3618 if (!snes->ksp) {ierr = SNESGetKSP(snes,&snes->ksp);CHKERRQ(ierr);} 3619 ierr = KSPSetOptionsPrefix(snes->ksp,prefix);CHKERRQ(ierr); 3620 PetscFunctionReturn(0); 3621 } 3622 3623 #undef __FUNCT__ 3624 #define __FUNCT__ "SNESAppendOptionsPrefix" 3625 /*@C 3626 SNESAppendOptionsPrefix - Appends to the prefix used for searching for all 3627 SNES options in the database. 3628 3629 Logically Collective on SNES 3630 3631 Input Parameters: 3632 + snes - the SNES context 3633 - prefix - the prefix to prepend to all option names 3634 3635 Notes: 3636 A hyphen (-) must NOT be given at the beginning of the prefix name. 3637 The first character of all runtime options is AUTOMATICALLY the hyphen. 3638 3639 Level: advanced 3640 3641 .keywords: SNES, append, options, prefix, database 3642 3643 .seealso: SNESGetOptionsPrefix() 3644 @*/ 3645 PetscErrorCode SNESAppendOptionsPrefix(SNES snes,const char prefix[]) 3646 { 3647 PetscErrorCode ierr; 3648 3649 PetscFunctionBegin; 3650 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3651 ierr = PetscObjectAppendOptionsPrefix((PetscObject)snes,prefix);CHKERRQ(ierr); 3652 if (!snes->ksp) {ierr = SNESGetKSP(snes,&snes->ksp);CHKERRQ(ierr);} 3653 ierr = KSPAppendOptionsPrefix(snes->ksp,prefix);CHKERRQ(ierr); 3654 PetscFunctionReturn(0); 3655 } 3656 3657 #undef __FUNCT__ 3658 #define __FUNCT__ "SNESGetOptionsPrefix" 3659 /*@C 3660 SNESGetOptionsPrefix - Sets the prefix used for searching for all 3661 SNES options in the database. 3662 3663 Not Collective 3664 3665 Input Parameter: 3666 . snes - the SNES context 3667 3668 Output Parameter: 3669 . prefix - pointer to the prefix string used 3670 3671 Notes: On the fortran side, the user should pass in a string 'prefix' of 3672 sufficient length to hold the prefix. 3673 3674 Level: advanced 3675 3676 .keywords: SNES, get, options, prefix, database 3677 3678 .seealso: SNESAppendOptionsPrefix() 3679 @*/ 3680 PetscErrorCode SNESGetOptionsPrefix(SNES snes,const char *prefix[]) 3681 { 3682 PetscErrorCode ierr; 3683 3684 PetscFunctionBegin; 3685 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3686 ierr = PetscObjectGetOptionsPrefix((PetscObject)snes,prefix);CHKERRQ(ierr); 3687 PetscFunctionReturn(0); 3688 } 3689 3690 3691 #undef __FUNCT__ 3692 #define __FUNCT__ "SNESRegister" 3693 /*@C 3694 SNESRegister - See SNESRegisterDynamic() 3695 3696 Level: advanced 3697 @*/ 3698 PetscErrorCode SNESRegister(const char sname[],const char path[],const char name[],PetscErrorCode (*function)(SNES)) 3699 { 3700 char fullname[PETSC_MAX_PATH_LEN]; 3701 PetscErrorCode ierr; 3702 3703 PetscFunctionBegin; 3704 ierr = PetscFListConcat(path,name,fullname);CHKERRQ(ierr); 3705 ierr = PetscFListAdd(&SNESList,sname,fullname,(void (*)(void))function);CHKERRQ(ierr); 3706 PetscFunctionReturn(0); 3707 } 3708 3709 #undef __FUNCT__ 3710 #define __FUNCT__ "SNESTestLocalMin" 3711 PetscErrorCode SNESTestLocalMin(SNES snes) 3712 { 3713 PetscErrorCode ierr; 3714 PetscInt N,i,j; 3715 Vec u,uh,fh; 3716 PetscScalar value; 3717 PetscReal norm; 3718 3719 PetscFunctionBegin; 3720 ierr = SNESGetSolution(snes,&u);CHKERRQ(ierr); 3721 ierr = VecDuplicate(u,&uh);CHKERRQ(ierr); 3722 ierr = VecDuplicate(u,&fh);CHKERRQ(ierr); 3723 3724 /* currently only works for sequential */ 3725 ierr = PetscPrintf(PETSC_COMM_WORLD,"Testing FormFunction() for local min\n"); 3726 ierr = VecGetSize(u,&N);CHKERRQ(ierr); 3727 for (i=0; i<N; i++) { 3728 ierr = VecCopy(u,uh);CHKERRQ(ierr); 3729 ierr = PetscPrintf(PETSC_COMM_WORLD,"i = %D\n",i);CHKERRQ(ierr); 3730 for (j=-10; j<11; j++) { 3731 value = PetscSign(j)*exp(PetscAbs(j)-10.0); 3732 ierr = VecSetValue(uh,i,value,ADD_VALUES);CHKERRQ(ierr); 3733 ierr = SNESComputeFunction(snes,uh,fh);CHKERRQ(ierr); 3734 ierr = VecNorm(fh,NORM_2,&norm);CHKERRQ(ierr); 3735 ierr = PetscPrintf(PETSC_COMM_WORLD," j norm %D %18.16e\n",j,norm);CHKERRQ(ierr); 3736 value = -value; 3737 ierr = VecSetValue(uh,i,value,ADD_VALUES);CHKERRQ(ierr); 3738 } 3739 } 3740 ierr = VecDestroy(&uh);CHKERRQ(ierr); 3741 ierr = VecDestroy(&fh);CHKERRQ(ierr); 3742 PetscFunctionReturn(0); 3743 } 3744 3745 #undef __FUNCT__ 3746 #define __FUNCT__ "SNESKSPSetUseEW" 3747 /*@ 3748 SNESKSPSetUseEW - Sets SNES use Eisenstat-Walker method for 3749 computing relative tolerance for linear solvers within an inexact 3750 Newton method. 3751 3752 Logically Collective on SNES 3753 3754 Input Parameters: 3755 + snes - SNES context 3756 - flag - PETSC_TRUE or PETSC_FALSE 3757 3758 Options Database: 3759 + -snes_ksp_ew - use Eisenstat-Walker method for determining linear system convergence 3760 . -snes_ksp_ew_version ver - version of Eisenstat-Walker method 3761 . -snes_ksp_ew_rtol0 <rtol0> - Sets rtol0 3762 . -snes_ksp_ew_rtolmax <rtolmax> - Sets rtolmax 3763 . -snes_ksp_ew_gamma <gamma> - Sets gamma 3764 . -snes_ksp_ew_alpha <alpha> - Sets alpha 3765 . -snes_ksp_ew_alpha2 <alpha2> - Sets alpha2 3766 - -snes_ksp_ew_threshold <threshold> - Sets threshold 3767 3768 Notes: 3769 Currently, the default is to use a constant relative tolerance for 3770 the inner linear solvers. Alternatively, one can use the 3771 Eisenstat-Walker method, where the relative convergence tolerance 3772 is reset at each Newton iteration according progress of the nonlinear 3773 solver. 3774 3775 Level: advanced 3776 3777 Reference: 3778 S. C. Eisenstat and H. F. Walker, "Choosing the forcing terms in an 3779 inexact Newton method", SISC 17 (1), pp.16-32, 1996. 3780 3781 .keywords: SNES, KSP, Eisenstat, Walker, convergence, test, inexact, Newton 3782 3783 .seealso: SNESKSPGetUseEW(), SNESKSPGetParametersEW(), SNESKSPSetParametersEW() 3784 @*/ 3785 PetscErrorCode SNESKSPSetUseEW(SNES snes,PetscBool flag) 3786 { 3787 PetscFunctionBegin; 3788 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3789 PetscValidLogicalCollectiveBool(snes,flag,2); 3790 snes->ksp_ewconv = flag; 3791 PetscFunctionReturn(0); 3792 } 3793 3794 #undef __FUNCT__ 3795 #define __FUNCT__ "SNESKSPGetUseEW" 3796 /*@ 3797 SNESKSPGetUseEW - Gets if SNES is using Eisenstat-Walker method 3798 for computing relative tolerance for linear solvers within an 3799 inexact Newton method. 3800 3801 Not Collective 3802 3803 Input Parameter: 3804 . snes - SNES context 3805 3806 Output Parameter: 3807 . flag - PETSC_TRUE or PETSC_FALSE 3808 3809 Notes: 3810 Currently, the default is to use a constant relative tolerance for 3811 the inner linear solvers. Alternatively, one can use the 3812 Eisenstat-Walker method, where the relative convergence tolerance 3813 is reset at each Newton iteration according progress of the nonlinear 3814 solver. 3815 3816 Level: advanced 3817 3818 Reference: 3819 S. C. Eisenstat and H. F. Walker, "Choosing the forcing terms in an 3820 inexact Newton method", SISC 17 (1), pp.16-32, 1996. 3821 3822 .keywords: SNES, KSP, Eisenstat, Walker, convergence, test, inexact, Newton 3823 3824 .seealso: SNESKSPSetUseEW(), SNESKSPGetParametersEW(), SNESKSPSetParametersEW() 3825 @*/ 3826 PetscErrorCode SNESKSPGetUseEW(SNES snes, PetscBool *flag) 3827 { 3828 PetscFunctionBegin; 3829 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3830 PetscValidPointer(flag,2); 3831 *flag = snes->ksp_ewconv; 3832 PetscFunctionReturn(0); 3833 } 3834 3835 #undef __FUNCT__ 3836 #define __FUNCT__ "SNESKSPSetParametersEW" 3837 /*@ 3838 SNESKSPSetParametersEW - Sets parameters for Eisenstat-Walker 3839 convergence criteria for the linear solvers within an inexact 3840 Newton method. 3841 3842 Logically Collective on SNES 3843 3844 Input Parameters: 3845 + snes - SNES context 3846 . version - version 1, 2 (default is 2) or 3 3847 . rtol_0 - initial relative tolerance (0 <= rtol_0 < 1) 3848 . rtol_max - maximum relative tolerance (0 <= rtol_max < 1) 3849 . gamma - multiplicative factor for version 2 rtol computation 3850 (0 <= gamma2 <= 1) 3851 . alpha - power for version 2 rtol computation (1 < alpha <= 2) 3852 . alpha2 - power for safeguard 3853 - threshold - threshold for imposing safeguard (0 < threshold < 1) 3854 3855 Note: 3856 Version 3 was contributed by Luis Chacon, June 2006. 3857 3858 Use PETSC_DEFAULT to retain the default for any of the parameters. 3859 3860 Level: advanced 3861 3862 Reference: 3863 S. C. Eisenstat and H. F. Walker, "Choosing the forcing terms in an 3864 inexact Newton method", Utah State University Math. Stat. Dept. Res. 3865 Report 6/94/75, June, 1994, to appear in SIAM J. Sci. Comput. 3866 3867 .keywords: SNES, KSP, Eisenstat, Walker, set, parameters 3868 3869 .seealso: SNESKSPSetUseEW(), SNESKSPGetUseEW(), SNESKSPGetParametersEW() 3870 @*/ 3871 PetscErrorCode SNESKSPSetParametersEW(SNES snes,PetscInt version,PetscReal rtol_0,PetscReal rtol_max, 3872 PetscReal gamma,PetscReal alpha,PetscReal alpha2,PetscReal threshold) 3873 { 3874 SNESKSPEW *kctx; 3875 PetscFunctionBegin; 3876 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3877 kctx = (SNESKSPEW*)snes->kspconvctx; 3878 if (!kctx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"No Eisenstat-Walker context existing"); 3879 PetscValidLogicalCollectiveInt(snes,version,2); 3880 PetscValidLogicalCollectiveReal(snes,rtol_0,3); 3881 PetscValidLogicalCollectiveReal(snes,rtol_max,4); 3882 PetscValidLogicalCollectiveReal(snes,gamma,5); 3883 PetscValidLogicalCollectiveReal(snes,alpha,6); 3884 PetscValidLogicalCollectiveReal(snes,alpha2,7); 3885 PetscValidLogicalCollectiveReal(snes,threshold,8); 3886 3887 if (version != PETSC_DEFAULT) kctx->version = version; 3888 if (rtol_0 != PETSC_DEFAULT) kctx->rtol_0 = rtol_0; 3889 if (rtol_max != PETSC_DEFAULT) kctx->rtol_max = rtol_max; 3890 if (gamma != PETSC_DEFAULT) kctx->gamma = gamma; 3891 if (alpha != PETSC_DEFAULT) kctx->alpha = alpha; 3892 if (alpha2 != PETSC_DEFAULT) kctx->alpha2 = alpha2; 3893 if (threshold != PETSC_DEFAULT) kctx->threshold = threshold; 3894 3895 if (kctx->version < 1 || kctx->version > 3) { 3896 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Only versions 1, 2 and 3 are supported: %D",kctx->version); 3897 } 3898 if (kctx->rtol_0 < 0.0 || kctx->rtol_0 >= 1.0) { 3899 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"0.0 <= rtol_0 < 1.0: %G",kctx->rtol_0); 3900 } 3901 if (kctx->rtol_max < 0.0 || kctx->rtol_max >= 1.0) { 3902 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"0.0 <= rtol_max (%G) < 1.0\n",kctx->rtol_max); 3903 } 3904 if (kctx->gamma < 0.0 || kctx->gamma > 1.0) { 3905 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"0.0 <= gamma (%G) <= 1.0\n",kctx->gamma); 3906 } 3907 if (kctx->alpha <= 1.0 || kctx->alpha > 2.0) { 3908 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"1.0 < alpha (%G) <= 2.0\n",kctx->alpha); 3909 } 3910 if (kctx->threshold <= 0.0 || kctx->threshold >= 1.0) { 3911 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"0.0 < threshold (%G) < 1.0\n",kctx->threshold); 3912 } 3913 PetscFunctionReturn(0); 3914 } 3915 3916 #undef __FUNCT__ 3917 #define __FUNCT__ "SNESKSPGetParametersEW" 3918 /*@ 3919 SNESKSPGetParametersEW - Gets parameters for Eisenstat-Walker 3920 convergence criteria for the linear solvers within an inexact 3921 Newton method. 3922 3923 Not Collective 3924 3925 Input Parameters: 3926 snes - SNES context 3927 3928 Output Parameters: 3929 + version - version 1, 2 (default is 2) or 3 3930 . rtol_0 - initial relative tolerance (0 <= rtol_0 < 1) 3931 . rtol_max - maximum relative tolerance (0 <= rtol_max < 1) 3932 . gamma - multiplicative factor for version 2 rtol computation 3933 (0 <= gamma2 <= 1) 3934 . alpha - power for version 2 rtol computation (1 < alpha <= 2) 3935 . alpha2 - power for safeguard 3936 - threshold - threshold for imposing safeguard (0 < threshold < 1) 3937 3938 Level: advanced 3939 3940 .keywords: SNES, KSP, Eisenstat, Walker, get, parameters 3941 3942 .seealso: SNESKSPSetUseEW(), SNESKSPGetUseEW(), SNESKSPSetParametersEW() 3943 @*/ 3944 PetscErrorCode SNESKSPGetParametersEW(SNES snes,PetscInt *version,PetscReal *rtol_0,PetscReal *rtol_max, 3945 PetscReal *gamma,PetscReal *alpha,PetscReal *alpha2,PetscReal *threshold) 3946 { 3947 SNESKSPEW *kctx; 3948 PetscFunctionBegin; 3949 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3950 kctx = (SNESKSPEW*)snes->kspconvctx; 3951 if (!kctx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"No Eisenstat-Walker context existing"); 3952 if(version) *version = kctx->version; 3953 if(rtol_0) *rtol_0 = kctx->rtol_0; 3954 if(rtol_max) *rtol_max = kctx->rtol_max; 3955 if(gamma) *gamma = kctx->gamma; 3956 if(alpha) *alpha = kctx->alpha; 3957 if(alpha2) *alpha2 = kctx->alpha2; 3958 if(threshold) *threshold = kctx->threshold; 3959 PetscFunctionReturn(0); 3960 } 3961 3962 #undef __FUNCT__ 3963 #define __FUNCT__ "SNESKSPEW_PreSolve" 3964 static PetscErrorCode SNESKSPEW_PreSolve(SNES snes, KSP ksp, Vec b, Vec x) 3965 { 3966 PetscErrorCode ierr; 3967 SNESKSPEW *kctx = (SNESKSPEW*)snes->kspconvctx; 3968 PetscReal rtol=PETSC_DEFAULT,stol; 3969 3970 PetscFunctionBegin; 3971 if (!kctx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"No Eisenstat-Walker context exists"); 3972 if (!snes->iter) { /* first time in, so use the original user rtol */ 3973 rtol = kctx->rtol_0; 3974 } else { 3975 if (kctx->version == 1) { 3976 rtol = (snes->norm - kctx->lresid_last)/kctx->norm_last; 3977 if (rtol < 0.0) rtol = -rtol; 3978 stol = pow(kctx->rtol_last,kctx->alpha2); 3979 if (stol > kctx->threshold) rtol = PetscMax(rtol,stol); 3980 } else if (kctx->version == 2) { 3981 rtol = kctx->gamma * pow(snes->norm/kctx->norm_last,kctx->alpha); 3982 stol = kctx->gamma * pow(kctx->rtol_last,kctx->alpha); 3983 if (stol > kctx->threshold) rtol = PetscMax(rtol,stol); 3984 } else if (kctx->version == 3) {/* contributed by Luis Chacon, June 2006. */ 3985 rtol = kctx->gamma * pow(snes->norm/kctx->norm_last,kctx->alpha); 3986 /* safeguard: avoid sharp decrease of rtol */ 3987 stol = kctx->gamma*pow(kctx->rtol_last,kctx->alpha); 3988 stol = PetscMax(rtol,stol); 3989 rtol = PetscMin(kctx->rtol_0,stol); 3990 /* safeguard: avoid oversolving */ 3991 stol = kctx->gamma*(snes->ttol)/snes->norm; 3992 stol = PetscMax(rtol,stol); 3993 rtol = PetscMin(kctx->rtol_0,stol); 3994 } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Only versions 1, 2 or 3 are supported: %D",kctx->version); 3995 } 3996 /* safeguard: avoid rtol greater than one */ 3997 rtol = PetscMin(rtol,kctx->rtol_max); 3998 ierr = KSPSetTolerances(ksp,rtol,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);CHKERRQ(ierr); 3999 ierr = PetscInfo3(snes,"iter %D, Eisenstat-Walker (version %D) KSP rtol=%G\n",snes->iter,kctx->version,rtol);CHKERRQ(ierr); 4000 PetscFunctionReturn(0); 4001 } 4002 4003 #undef __FUNCT__ 4004 #define __FUNCT__ "SNESKSPEW_PostSolve" 4005 static PetscErrorCode SNESKSPEW_PostSolve(SNES snes, KSP ksp, Vec b, Vec x) 4006 { 4007 PetscErrorCode ierr; 4008 SNESKSPEW *kctx = (SNESKSPEW*)snes->kspconvctx; 4009 PCSide pcside; 4010 Vec lres; 4011 4012 PetscFunctionBegin; 4013 if (!kctx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"No Eisenstat-Walker context exists"); 4014 ierr = KSPGetTolerances(ksp,&kctx->rtol_last,0,0,0);CHKERRQ(ierr); 4015 ierr = SNESGetFunctionNorm(snes,&kctx->norm_last);CHKERRQ(ierr); 4016 if (kctx->version == 1) { 4017 ierr = KSPGetPCSide(ksp,&pcside);CHKERRQ(ierr); 4018 if (pcside == PC_RIGHT) { /* XXX Should we also test KSP_UNPRECONDITIONED_NORM ? */ 4019 /* KSP residual is true linear residual */ 4020 ierr = KSPGetResidualNorm(ksp,&kctx->lresid_last);CHKERRQ(ierr); 4021 } else { 4022 /* KSP residual is preconditioned residual */ 4023 /* compute true linear residual norm */ 4024 ierr = VecDuplicate(b,&lres);CHKERRQ(ierr); 4025 ierr = MatMult(snes->jacobian,x,lres);CHKERRQ(ierr); 4026 ierr = VecAYPX(lres,-1.0,b);CHKERRQ(ierr); 4027 ierr = VecNorm(lres,NORM_2,&kctx->lresid_last);CHKERRQ(ierr); 4028 ierr = VecDestroy(&lres);CHKERRQ(ierr); 4029 } 4030 } 4031 PetscFunctionReturn(0); 4032 } 4033 4034 #undef __FUNCT__ 4035 #define __FUNCT__ "SNES_KSPSolve" 4036 PetscErrorCode SNES_KSPSolve(SNES snes, KSP ksp, Vec b, Vec x) 4037 { 4038 PetscErrorCode ierr; 4039 4040 PetscFunctionBegin; 4041 if (snes->ksp_ewconv) { ierr = SNESKSPEW_PreSolve(snes,ksp,b,x);CHKERRQ(ierr); } 4042 ierr = KSPSolve(ksp,b,x);CHKERRQ(ierr); 4043 if (snes->ksp_ewconv) { ierr = SNESKSPEW_PostSolve(snes,ksp,b,x);CHKERRQ(ierr); } 4044 PetscFunctionReturn(0); 4045 } 4046 4047 #undef __FUNCT__ 4048 #define __FUNCT__ "SNESSetDM" 4049 /*@ 4050 SNESSetDM - Sets the DM that may be used by some preconditioners 4051 4052 Logically Collective on SNES 4053 4054 Input Parameters: 4055 + snes - the preconditioner context 4056 - dm - the dm 4057 4058 Level: intermediate 4059 4060 4061 .seealso: SNESGetDM(), KSPSetDM(), KSPGetDM() 4062 @*/ 4063 PetscErrorCode SNESSetDM(SNES snes,DM dm) 4064 { 4065 PetscErrorCode ierr; 4066 KSP ksp; 4067 SNESDM sdm; 4068 4069 PetscFunctionBegin; 4070 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4071 if (dm) {ierr = PetscObjectReference((PetscObject)dm);CHKERRQ(ierr);} 4072 if (snes->dm) { /* Move the SNESDM context over to the new DM unless the new DM already has one */ 4073 PetscContainer oldcontainer,container; 4074 ierr = PetscObjectQuery((PetscObject)snes->dm,"SNESDM",(PetscObject*)&oldcontainer);CHKERRQ(ierr); 4075 ierr = PetscObjectQuery((PetscObject)dm,"SNESDM",(PetscObject*)&container);CHKERRQ(ierr); 4076 if (oldcontainer && !container) { 4077 ierr = DMSNESCopyContext(snes->dm,dm);CHKERRQ(ierr); 4078 ierr = DMSNESGetContext(snes->dm,&sdm);CHKERRQ(ierr); 4079 if (sdm->originaldm == snes->dm) { /* Grant write privileges to the replacement DM */ 4080 sdm->originaldm = dm; 4081 } 4082 } 4083 ierr = DMDestroy(&snes->dm);CHKERRQ(ierr); 4084 } 4085 snes->dm = dm; 4086 ierr = SNESGetKSP(snes,&ksp);CHKERRQ(ierr); 4087 ierr = KSPSetDM(ksp,dm);CHKERRQ(ierr); 4088 ierr = KSPSetDMActive(ksp,PETSC_FALSE);CHKERRQ(ierr); 4089 if (snes->pc) { 4090 ierr = SNESSetDM(snes->pc, snes->dm);CHKERRQ(ierr); 4091 } 4092 PetscFunctionReturn(0); 4093 } 4094 4095 #undef __FUNCT__ 4096 #define __FUNCT__ "SNESGetDM" 4097 /*@ 4098 SNESGetDM - Gets the DM that may be used by some preconditioners 4099 4100 Not Collective but DM obtained is parallel on SNES 4101 4102 Input Parameter: 4103 . snes - the preconditioner context 4104 4105 Output Parameter: 4106 . dm - the dm 4107 4108 Level: intermediate 4109 4110 4111 .seealso: SNESSetDM(), KSPSetDM(), KSPGetDM() 4112 @*/ 4113 PetscErrorCode SNESGetDM(SNES snes,DM *dm) 4114 { 4115 PetscErrorCode ierr; 4116 4117 PetscFunctionBegin; 4118 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4119 if (!snes->dm) { 4120 ierr = DMShellCreate(((PetscObject)snes)->comm,&snes->dm);CHKERRQ(ierr); 4121 } 4122 *dm = snes->dm; 4123 PetscFunctionReturn(0); 4124 } 4125 4126 #undef __FUNCT__ 4127 #define __FUNCT__ "SNESSetPC" 4128 /*@ 4129 SNESSetPC - Sets the nonlinear preconditioner to be used. 4130 4131 Collective on SNES 4132 4133 Input Parameters: 4134 + snes - iterative context obtained from SNESCreate() 4135 - pc - the preconditioner object 4136 4137 Notes: 4138 Use SNESGetPC() to retrieve the preconditioner context (for example, 4139 to configure it using the API). 4140 4141 Level: developer 4142 4143 .keywords: SNES, set, precondition 4144 .seealso: SNESGetPC() 4145 @*/ 4146 PetscErrorCode SNESSetPC(SNES snes, SNES pc) 4147 { 4148 PetscErrorCode ierr; 4149 4150 PetscFunctionBegin; 4151 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); 4152 PetscValidHeaderSpecific(pc, SNES_CLASSID, 2); 4153 PetscCheckSameComm(snes, 1, pc, 2); 4154 ierr = PetscObjectReference((PetscObject) pc);CHKERRQ(ierr); 4155 ierr = SNESDestroy(&snes->pc);CHKERRQ(ierr); 4156 snes->pc = pc; 4157 ierr = PetscLogObjectParent(snes, snes->pc);CHKERRQ(ierr); 4158 PetscFunctionReturn(0); 4159 } 4160 4161 #undef __FUNCT__ 4162 #define __FUNCT__ "SNESGetPC" 4163 /*@ 4164 SNESGetPC - Returns a pointer to the nonlinear preconditioning context set with SNESSetPC(). 4165 4166 Not Collective 4167 4168 Input Parameter: 4169 . snes - iterative context obtained from SNESCreate() 4170 4171 Output Parameter: 4172 . pc - preconditioner context 4173 4174 Level: developer 4175 4176 .keywords: SNES, get, preconditioner 4177 .seealso: SNESSetPC() 4178 @*/ 4179 PetscErrorCode SNESGetPC(SNES snes, SNES *pc) 4180 { 4181 PetscErrorCode ierr; 4182 4183 PetscFunctionBegin; 4184 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); 4185 PetscValidPointer(pc, 2); 4186 if (!snes->pc) { 4187 ierr = SNESCreate(((PetscObject) snes)->comm, &snes->pc);CHKERRQ(ierr); 4188 ierr = PetscObjectIncrementTabLevel((PetscObject) snes->pc, (PetscObject) snes, 1);CHKERRQ(ierr); 4189 ierr = PetscLogObjectParent(snes, snes->pc);CHKERRQ(ierr); 4190 } 4191 *pc = snes->pc; 4192 PetscFunctionReturn(0); 4193 } 4194 4195 #if defined(PETSC_HAVE_MATLAB_ENGINE) 4196 #include <mex.h> 4197 4198 typedef struct {char *funcname; mxArray *ctx;} SNESMatlabContext; 4199 4200 #undef __FUNCT__ 4201 #define __FUNCT__ "SNESComputeFunction_Matlab" 4202 /* 4203 SNESComputeFunction_Matlab - Calls the function that has been set with 4204 SNESSetFunctionMatlab(). 4205 4206 Collective on SNES 4207 4208 Input Parameters: 4209 + snes - the SNES context 4210 - x - input vector 4211 4212 Output Parameter: 4213 . y - function vector, as set by SNESSetFunction() 4214 4215 Notes: 4216 SNESComputeFunction() is typically used within nonlinear solvers 4217 implementations, so most users would not generally call this routine 4218 themselves. 4219 4220 Level: developer 4221 4222 .keywords: SNES, nonlinear, compute, function 4223 4224 .seealso: SNESSetFunction(), SNESGetFunction() 4225 */ 4226 PetscErrorCode SNESComputeFunction_Matlab(SNES snes,Vec x,Vec y, void *ctx) 4227 { 4228 PetscErrorCode ierr; 4229 SNESMatlabContext *sctx = (SNESMatlabContext *)ctx; 4230 int nlhs = 1,nrhs = 5; 4231 mxArray *plhs[1],*prhs[5]; 4232 long long int lx = 0,ly = 0,ls = 0; 4233 4234 PetscFunctionBegin; 4235 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4236 PetscValidHeaderSpecific(x,VEC_CLASSID,2); 4237 PetscValidHeaderSpecific(y,VEC_CLASSID,3); 4238 PetscCheckSameComm(snes,1,x,2); 4239 PetscCheckSameComm(snes,1,y,3); 4240 4241 /* call Matlab function in ctx with arguments x and y */ 4242 4243 ierr = PetscMemcpy(&ls,&snes,sizeof(snes));CHKERRQ(ierr); 4244 ierr = PetscMemcpy(&lx,&x,sizeof(x));CHKERRQ(ierr); 4245 ierr = PetscMemcpy(&ly,&y,sizeof(x));CHKERRQ(ierr); 4246 prhs[0] = mxCreateDoubleScalar((double)ls); 4247 prhs[1] = mxCreateDoubleScalar((double)lx); 4248 prhs[2] = mxCreateDoubleScalar((double)ly); 4249 prhs[3] = mxCreateString(sctx->funcname); 4250 prhs[4] = sctx->ctx; 4251 ierr = mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscSNESComputeFunctionInternal");CHKERRQ(ierr); 4252 ierr = mxGetScalar(plhs[0]);CHKERRQ(ierr); 4253 mxDestroyArray(prhs[0]); 4254 mxDestroyArray(prhs[1]); 4255 mxDestroyArray(prhs[2]); 4256 mxDestroyArray(prhs[3]); 4257 mxDestroyArray(plhs[0]); 4258 PetscFunctionReturn(0); 4259 } 4260 4261 4262 #undef __FUNCT__ 4263 #define __FUNCT__ "SNESSetFunctionMatlab" 4264 /* 4265 SNESSetFunctionMatlab - Sets the function evaluation routine and function 4266 vector for use by the SNES routines in solving systems of nonlinear 4267 equations from MATLAB. Here the function is a string containing the name of a MATLAB function 4268 4269 Logically Collective on SNES 4270 4271 Input Parameters: 4272 + snes - the SNES context 4273 . r - vector to store function value 4274 - func - function evaluation routine 4275 4276 Calling sequence of func: 4277 $ func (SNES snes,Vec x,Vec f,void *ctx); 4278 4279 4280 Notes: 4281 The Newton-like methods typically solve linear systems of the form 4282 $ f'(x) x = -f(x), 4283 where f'(x) denotes the Jacobian matrix and f(x) is the function. 4284 4285 Level: beginner 4286 4287 .keywords: SNES, nonlinear, set, function 4288 4289 .seealso: SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESSetFunction() 4290 */ 4291 PetscErrorCode SNESSetFunctionMatlab(SNES snes,Vec r,const char *func,mxArray *ctx) 4292 { 4293 PetscErrorCode ierr; 4294 SNESMatlabContext *sctx; 4295 4296 PetscFunctionBegin; 4297 /* currently sctx is memory bleed */ 4298 ierr = PetscMalloc(sizeof(SNESMatlabContext),&sctx);CHKERRQ(ierr); 4299 ierr = PetscStrallocpy(func,&sctx->funcname);CHKERRQ(ierr); 4300 /* 4301 This should work, but it doesn't 4302 sctx->ctx = ctx; 4303 mexMakeArrayPersistent(sctx->ctx); 4304 */ 4305 sctx->ctx = mxDuplicateArray(ctx); 4306 ierr = SNESSetFunction(snes,r,SNESComputeFunction_Matlab,sctx);CHKERRQ(ierr); 4307 PetscFunctionReturn(0); 4308 } 4309 4310 #undef __FUNCT__ 4311 #define __FUNCT__ "SNESComputeJacobian_Matlab" 4312 /* 4313 SNESComputeJacobian_Matlab - Calls the function that has been set with 4314 SNESSetJacobianMatlab(). 4315 4316 Collective on SNES 4317 4318 Input Parameters: 4319 + snes - the SNES context 4320 . x - input vector 4321 . A, B - the matrices 4322 - ctx - user context 4323 4324 Output Parameter: 4325 . flag - structure of the matrix 4326 4327 Level: developer 4328 4329 .keywords: SNES, nonlinear, compute, function 4330 4331 .seealso: SNESSetFunction(), SNESGetFunction() 4332 @*/ 4333 PetscErrorCode SNESComputeJacobian_Matlab(SNES snes,Vec x,Mat *A,Mat *B,MatStructure *flag, void *ctx) 4334 { 4335 PetscErrorCode ierr; 4336 SNESMatlabContext *sctx = (SNESMatlabContext *)ctx; 4337 int nlhs = 2,nrhs = 6; 4338 mxArray *plhs[2],*prhs[6]; 4339 long long int lx = 0,lA = 0,ls = 0, lB = 0; 4340 4341 PetscFunctionBegin; 4342 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4343 PetscValidHeaderSpecific(x,VEC_CLASSID,2); 4344 4345 /* call Matlab function in ctx with arguments x and y */ 4346 4347 ierr = PetscMemcpy(&ls,&snes,sizeof(snes));CHKERRQ(ierr); 4348 ierr = PetscMemcpy(&lx,&x,sizeof(x));CHKERRQ(ierr); 4349 ierr = PetscMemcpy(&lA,A,sizeof(x));CHKERRQ(ierr); 4350 ierr = PetscMemcpy(&lB,B,sizeof(x));CHKERRQ(ierr); 4351 prhs[0] = mxCreateDoubleScalar((double)ls); 4352 prhs[1] = mxCreateDoubleScalar((double)lx); 4353 prhs[2] = mxCreateDoubleScalar((double)lA); 4354 prhs[3] = mxCreateDoubleScalar((double)lB); 4355 prhs[4] = mxCreateString(sctx->funcname); 4356 prhs[5] = sctx->ctx; 4357 ierr = mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscSNESComputeJacobianInternal");CHKERRQ(ierr); 4358 ierr = mxGetScalar(plhs[0]);CHKERRQ(ierr); 4359 *flag = (MatStructure) mxGetScalar(plhs[1]);CHKERRQ(ierr); 4360 mxDestroyArray(prhs[0]); 4361 mxDestroyArray(prhs[1]); 4362 mxDestroyArray(prhs[2]); 4363 mxDestroyArray(prhs[3]); 4364 mxDestroyArray(prhs[4]); 4365 mxDestroyArray(plhs[0]); 4366 mxDestroyArray(plhs[1]); 4367 PetscFunctionReturn(0); 4368 } 4369 4370 4371 #undef __FUNCT__ 4372 #define __FUNCT__ "SNESSetJacobianMatlab" 4373 /* 4374 SNESSetJacobianMatlab - Sets the Jacobian function evaluation routine and two empty Jacobian matrices 4375 vector for use by the SNES routines in solving systems of nonlinear 4376 equations from MATLAB. Here the function is a string containing the name of a MATLAB function 4377 4378 Logically Collective on SNES 4379 4380 Input Parameters: 4381 + snes - the SNES context 4382 . A,B - Jacobian matrices 4383 . func - function evaluation routine 4384 - ctx - user context 4385 4386 Calling sequence of func: 4387 $ flag = func (SNES snes,Vec x,Mat A,Mat B,void *ctx); 4388 4389 4390 Level: developer 4391 4392 .keywords: SNES, nonlinear, set, function 4393 4394 .seealso: SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESSetFunction() 4395 */ 4396 PetscErrorCode SNESSetJacobianMatlab(SNES snes,Mat A,Mat B,const char *func,mxArray *ctx) 4397 { 4398 PetscErrorCode ierr; 4399 SNESMatlabContext *sctx; 4400 4401 PetscFunctionBegin; 4402 /* currently sctx is memory bleed */ 4403 ierr = PetscMalloc(sizeof(SNESMatlabContext),&sctx);CHKERRQ(ierr); 4404 ierr = PetscStrallocpy(func,&sctx->funcname);CHKERRQ(ierr); 4405 /* 4406 This should work, but it doesn't 4407 sctx->ctx = ctx; 4408 mexMakeArrayPersistent(sctx->ctx); 4409 */ 4410 sctx->ctx = mxDuplicateArray(ctx); 4411 ierr = SNESSetJacobian(snes,A,B,SNESComputeJacobian_Matlab,sctx);CHKERRQ(ierr); 4412 PetscFunctionReturn(0); 4413 } 4414 4415 #undef __FUNCT__ 4416 #define __FUNCT__ "SNESMonitor_Matlab" 4417 /* 4418 SNESMonitor_Matlab - Calls the function that has been set with SNESMonitorSetMatlab(). 4419 4420 Collective on SNES 4421 4422 .seealso: SNESSetFunction(), SNESGetFunction() 4423 @*/ 4424 PetscErrorCode SNESMonitor_Matlab(SNES snes,PetscInt it, PetscReal fnorm, void *ctx) 4425 { 4426 PetscErrorCode ierr; 4427 SNESMatlabContext *sctx = (SNESMatlabContext *)ctx; 4428 int nlhs = 1,nrhs = 6; 4429 mxArray *plhs[1],*prhs[6]; 4430 long long int lx = 0,ls = 0; 4431 Vec x=snes->vec_sol; 4432 4433 PetscFunctionBegin; 4434 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4435 4436 ierr = PetscMemcpy(&ls,&snes,sizeof(snes));CHKERRQ(ierr); 4437 ierr = PetscMemcpy(&lx,&x,sizeof(x));CHKERRQ(ierr); 4438 prhs[0] = mxCreateDoubleScalar((double)ls); 4439 prhs[1] = mxCreateDoubleScalar((double)it); 4440 prhs[2] = mxCreateDoubleScalar((double)fnorm); 4441 prhs[3] = mxCreateDoubleScalar((double)lx); 4442 prhs[4] = mxCreateString(sctx->funcname); 4443 prhs[5] = sctx->ctx; 4444 ierr = mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscSNESMonitorInternal");CHKERRQ(ierr); 4445 ierr = mxGetScalar(plhs[0]);CHKERRQ(ierr); 4446 mxDestroyArray(prhs[0]); 4447 mxDestroyArray(prhs[1]); 4448 mxDestroyArray(prhs[2]); 4449 mxDestroyArray(prhs[3]); 4450 mxDestroyArray(prhs[4]); 4451 mxDestroyArray(plhs[0]); 4452 PetscFunctionReturn(0); 4453 } 4454 4455 4456 #undef __FUNCT__ 4457 #define __FUNCT__ "SNESMonitorSetMatlab" 4458 /* 4459 SNESMonitorSetMatlab - Sets the monitor function from MATLAB 4460 4461 Level: developer 4462 4463 .keywords: SNES, nonlinear, set, function 4464 4465 .seealso: SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESSetFunction() 4466 */ 4467 PetscErrorCode SNESMonitorSetMatlab(SNES snes,const char *func,mxArray *ctx) 4468 { 4469 PetscErrorCode ierr; 4470 SNESMatlabContext *sctx; 4471 4472 PetscFunctionBegin; 4473 /* currently sctx is memory bleed */ 4474 ierr = PetscMalloc(sizeof(SNESMatlabContext),&sctx);CHKERRQ(ierr); 4475 ierr = PetscStrallocpy(func,&sctx->funcname);CHKERRQ(ierr); 4476 /* 4477 This should work, but it doesn't 4478 sctx->ctx = ctx; 4479 mexMakeArrayPersistent(sctx->ctx); 4480 */ 4481 sctx->ctx = mxDuplicateArray(ctx); 4482 ierr = SNESMonitorSet(snes,SNESMonitor_Matlab,sctx,PETSC_NULL);CHKERRQ(ierr); 4483 PetscFunctionReturn(0); 4484 } 4485 4486 #endif 4487