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