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