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