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