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