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