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