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