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