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