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