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