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