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