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