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