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