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