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