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