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 1315 #undef __FUNCT__ 1316 #define __FUNCT__ "SNESSetKSP" 1317 /*@ 1318 SNESSetKSP - Sets a KSP context for the SNES object to use 1319 1320 Not Collective, but the SNES and KSP objects must live on the same MPI_Comm 1321 1322 Input Parameters: 1323 + snes - the SNES context 1324 - ksp - the KSP context 1325 1326 Notes: 1327 The SNES object already has its KSP object, you can obtain with SNESGetKSP() 1328 so this routine is rarely needed. 1329 1330 The KSP object that is already in the SNES object has its reference count 1331 decreased by one. 1332 1333 Level: developer 1334 1335 .keywords: SNES, nonlinear, get, KSP, context 1336 1337 .seealso: KSPGetPC(), SNESCreate(), KSPCreate(), SNESSetKSP() 1338 @*/ 1339 PetscErrorCode SNESSetKSP(SNES snes,KSP ksp) 1340 { 1341 PetscErrorCode ierr; 1342 1343 PetscFunctionBegin; 1344 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1345 PetscValidHeaderSpecific(ksp,KSP_CLASSID,2); 1346 PetscCheckSameComm(snes,1,ksp,2); 1347 ierr = PetscObjectReference((PetscObject)ksp);CHKERRQ(ierr); 1348 if (snes->ksp) {ierr = PetscObjectDereference((PetscObject)snes->ksp);CHKERRQ(ierr);} 1349 snes->ksp = ksp; 1350 PetscFunctionReturn(0); 1351 } 1352 1353 #if defined(PETSC_HAVE_AMS) 1354 #undef __FUNCT__ 1355 #define __FUNCT__ "SNESPublish_Petsc" 1356 static PetscErrorCode SNESPublish_Petsc(PetscObject obj) 1357 { 1358 PetscFunctionBegin; 1359 PetscFunctionReturn(0); 1360 } 1361 #endif 1362 1363 /* -----------------------------------------------------------*/ 1364 #undef __FUNCT__ 1365 #define __FUNCT__ "SNESCreate" 1366 /*@ 1367 SNESCreate - Creates a nonlinear solver context. 1368 1369 Collective on MPI_Comm 1370 1371 Input Parameters: 1372 . comm - MPI communicator 1373 1374 Output Parameter: 1375 . outsnes - the new SNES context 1376 1377 Options Database Keys: 1378 + -snes_mf - Activates default matrix-free Jacobian-vector products, 1379 and no preconditioning matrix 1380 . -snes_mf_operator - Activates default matrix-free Jacobian-vector 1381 products, and a user-provided preconditioning matrix 1382 as set by SNESSetJacobian() 1383 - -snes_fd - Uses (slow!) finite differences to compute Jacobian 1384 1385 Level: beginner 1386 1387 .keywords: SNES, nonlinear, create, context 1388 1389 .seealso: SNESSolve(), SNESDestroy(), SNES, SNESSetLagPreconditioner() 1390 1391 @*/ 1392 PetscErrorCode SNESCreate(MPI_Comm comm,SNES *outsnes) 1393 { 1394 PetscErrorCode ierr; 1395 SNES snes; 1396 SNESKSPEW *kctx; 1397 1398 PetscFunctionBegin; 1399 PetscValidPointer(outsnes,2); 1400 *outsnes = NULL; 1401 #if !defined(PETSC_USE_DYNAMIC_LIBRARIES) 1402 ierr = SNESInitializePackage(NULL);CHKERRQ(ierr); 1403 #endif 1404 1405 ierr = PetscHeaderCreate(snes,_p_SNES,struct _SNESOps,SNES_CLASSID,"SNES","Nonlinear solver","SNES",comm,SNESDestroy,SNESView);CHKERRQ(ierr); 1406 1407 snes->ops->converged = SNESConvergedDefault; 1408 snes->usesksp = PETSC_TRUE; 1409 snes->tolerancesset = PETSC_FALSE; 1410 snes->max_its = 50; 1411 snes->max_funcs = 10000; 1412 snes->norm = 0.0; 1413 snes->normtype = SNES_NORM_FUNCTION; 1414 snes->rtol = 1.e-8; 1415 snes->ttol = 0.0; 1416 snes->abstol = 1.e-50; 1417 snes->stol = 1.e-8; 1418 snes->deltatol = 1.e-12; 1419 snes->nfuncs = 0; 1420 snes->numFailures = 0; 1421 snes->maxFailures = 1; 1422 snes->linear_its = 0; 1423 snes->lagjacobian = 1; 1424 snes->lagpreconditioner = 1; 1425 snes->numbermonitors = 0; 1426 snes->data = 0; 1427 snes->setupcalled = PETSC_FALSE; 1428 snes->ksp_ewconv = PETSC_FALSE; 1429 snes->nwork = 0; 1430 snes->work = 0; 1431 snes->nvwork = 0; 1432 snes->vwork = 0; 1433 snes->conv_hist_len = 0; 1434 snes->conv_hist_max = 0; 1435 snes->conv_hist = NULL; 1436 snes->conv_hist_its = NULL; 1437 snes->conv_hist_reset = PETSC_TRUE; 1438 snes->vec_func_init_set = PETSC_FALSE; 1439 snes->norm_init = 0.; 1440 snes->norm_init_set = PETSC_FALSE; 1441 snes->reason = SNES_CONVERGED_ITERATING; 1442 snes->pcside = PC_RIGHT; 1443 1444 #if defined(PETSC_HAVE_AMS) 1445 ((PetscObject)snes)->bops->publish = SNESPublish_Petsc; 1446 #endif 1447 1448 snes->mf = PETSC_FALSE; 1449 snes->mf_operator = PETSC_FALSE; 1450 snes->mf_version = 1; 1451 1452 snes->numLinearSolveFailures = 0; 1453 snes->maxLinearSolveFailures = 1; 1454 1455 /* Create context to compute Eisenstat-Walker relative tolerance for KSP */ 1456 ierr = PetscNewLog(snes,SNESKSPEW,&kctx);CHKERRQ(ierr); 1457 1458 snes->kspconvctx = (void*)kctx; 1459 kctx->version = 2; 1460 kctx->rtol_0 = .3; /* Eisenstat and Walker suggest rtol_0=.5, but 1461 this was too large for some test cases */ 1462 kctx->rtol_last = 0.0; 1463 kctx->rtol_max = .9; 1464 kctx->gamma = 1.0; 1465 kctx->alpha = .5*(1.0 + PetscSqrtReal(5.0)); 1466 kctx->alpha2 = kctx->alpha; 1467 kctx->threshold = .1; 1468 kctx->lresid_last = 0.0; 1469 kctx->norm_last = 0.0; 1470 1471 *outsnes = snes; 1472 PetscFunctionReturn(0); 1473 } 1474 1475 /*MC 1476 SNESFunction - functional form used to convey the nonlinear function to be solved by SNES 1477 1478 Synopsis: 1479 #include "petscsnes.h" 1480 SNESFunction(SNES snes,Vec x,Vec f,void *ctx); 1481 1482 Input Parameters: 1483 + snes - the SNES context 1484 . x - state at which to evaluate residual 1485 - ctx - optional user-defined function context, passed in with SNESSetFunction() 1486 1487 Output Parameter: 1488 . f - vector to put residual (function value) 1489 1490 Level: intermediate 1491 1492 .seealso: SNESSetFunction(), SNESGetFunction() 1493 M*/ 1494 1495 #undef __FUNCT__ 1496 #define __FUNCT__ "SNESSetFunction" 1497 /*@C 1498 SNESSetFunction - Sets the function evaluation routine and function 1499 vector for use by the SNES routines in solving systems of nonlinear 1500 equations. 1501 1502 Logically Collective on SNES 1503 1504 Input Parameters: 1505 + snes - the SNES context 1506 . r - vector to store function value 1507 . SNESFunction - function evaluation routine 1508 - ctx - [optional] user-defined context for private data for the 1509 function evaluation routine (may be NULL) 1510 1511 Notes: 1512 The Newton-like methods typically solve linear systems of the form 1513 $ f'(x) x = -f(x), 1514 where f'(x) denotes the Jacobian matrix and f(x) is the function. 1515 1516 Level: beginner 1517 1518 .keywords: SNES, nonlinear, set, function 1519 1520 .seealso: SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESSetPicard(), SNESFunction 1521 @*/ 1522 PetscErrorCode SNESSetFunction(SNES snes,Vec r,PetscErrorCode (*SNESFunction)(SNES,Vec,Vec,void*),void *ctx) 1523 { 1524 PetscErrorCode ierr; 1525 DM dm; 1526 1527 PetscFunctionBegin; 1528 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1529 if (r) { 1530 PetscValidHeaderSpecific(r,VEC_CLASSID,2); 1531 PetscCheckSameComm(snes,1,r,2); 1532 ierr = PetscObjectReference((PetscObject)r);CHKERRQ(ierr); 1533 ierr = VecDestroy(&snes->vec_func);CHKERRQ(ierr); 1534 1535 snes->vec_func = r; 1536 } 1537 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 1538 ierr = DMSNESSetFunction(dm,SNESFunction,ctx);CHKERRQ(ierr); 1539 PetscFunctionReturn(0); 1540 } 1541 1542 1543 #undef __FUNCT__ 1544 #define __FUNCT__ "SNESSetInitialFunction" 1545 /*@C 1546 SNESSetInitialFunction - Sets the function vector to be used as the 1547 function norm at the initialization of the method. In some 1548 instances, the user has precomputed the function before calling 1549 SNESSolve. This function allows one to avoid a redundant call 1550 to SNESComputeFunction in that case. 1551 1552 Logically Collective on SNES 1553 1554 Input Parameters: 1555 + snes - the SNES context 1556 - f - vector to store function value 1557 1558 Notes: 1559 This should not be modified during the solution procedure. 1560 1561 This is used extensively in the SNESFAS hierarchy and in nonlinear preconditioning. 1562 1563 Level: developer 1564 1565 .keywords: SNES, nonlinear, set, function 1566 1567 .seealso: SNESSetFunction(), SNESComputeFunction(), SNESSetInitialFunctionNorm() 1568 @*/ 1569 PetscErrorCode SNESSetInitialFunction(SNES snes, Vec f) 1570 { 1571 PetscErrorCode ierr; 1572 Vec vec_func; 1573 1574 PetscFunctionBegin; 1575 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1576 PetscValidHeaderSpecific(f,VEC_CLASSID,2); 1577 PetscCheckSameComm(snes,1,f,2); 1578 ierr = SNESGetFunction(snes,&vec_func,NULL,NULL);CHKERRQ(ierr); 1579 ierr = VecCopy(f, vec_func);CHKERRQ(ierr); 1580 1581 snes->vec_func_init_set = PETSC_TRUE; 1582 PetscFunctionReturn(0); 1583 } 1584 1585 1586 #undef __FUNCT__ 1587 #define __FUNCT__ "SNESSetInitialFunctionNorm" 1588 /*@C 1589 SNESSetInitialFunctionNorm - Sets the function norm to be used as the function norm 1590 at the initialization of the method. In some instances, the user has precomputed 1591 the function and its norm before calling SNESSolve. This function allows one to 1592 avoid a redundant call to SNESComputeFunction() and VecNorm() in that case. 1593 1594 Logically Collective on SNES 1595 1596 Input Parameters: 1597 + snes - the SNES context 1598 - fnorm - the norm of F as set by SNESSetInitialFunction() 1599 1600 This is used extensively in the SNESFAS hierarchy and in nonlinear preconditioning. 1601 1602 Level: developer 1603 1604 .keywords: SNES, nonlinear, set, function, norm 1605 1606 .seealso: SNESSetFunction(), SNESComputeFunction(), SNESSetInitialFunction() 1607 @*/ 1608 PetscErrorCode SNESSetInitialFunctionNorm(SNES snes, PetscReal fnorm) 1609 { 1610 PetscFunctionBegin; 1611 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1612 snes->norm_init = fnorm; 1613 snes->norm_init_set = PETSC_TRUE; 1614 PetscFunctionReturn(0); 1615 } 1616 1617 #undef __FUNCT__ 1618 #define __FUNCT__ "SNESSetNormType" 1619 /*@ 1620 SNESSetNormType - Sets the SNESNormType used in covergence and monitoring 1621 of the SNES method. 1622 1623 Logically Collective on SNES 1624 1625 Input Parameters: 1626 + snes - the SNES context 1627 - normtype - the type of the norm used 1628 1629 Notes: 1630 Only certain SNES methods support certain SNESNormTypes. Most require evaluation 1631 of the nonlinear function and the taking of its norm at every iteration to 1632 even ensure convergence at all. However, methods such as custom Gauss-Seidel methods 1633 (SNESGS) and the like do not require the norm of the function to be computed, and therfore 1634 may either be monitored for convergence or not. As these are often used as nonlinear 1635 preconditioners, monitoring the norm of their error is not a useful enterprise within 1636 their solution. 1637 1638 Level: developer 1639 1640 .keywords: SNES, nonlinear, set, function, norm, type 1641 1642 .seealso: SNESGetNormType(), SNESComputeFunction(), VecNorm(), SNESSetFunction(), SNESSetInitialFunction(), SNESNormType 1643 @*/ 1644 PetscErrorCode SNESSetNormType(SNES snes, SNESNormType normtype) 1645 { 1646 PetscFunctionBegin; 1647 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1648 snes->normtype = normtype; 1649 PetscFunctionReturn(0); 1650 } 1651 1652 1653 #undef __FUNCT__ 1654 #define __FUNCT__ "SNESGetNormType" 1655 /*@ 1656 SNESGetNormType - Gets the SNESNormType used in covergence and monitoring 1657 of the SNES method. 1658 1659 Logically Collective on SNES 1660 1661 Input Parameters: 1662 + snes - the SNES context 1663 - normtype - the type of the norm used 1664 1665 Level: advanced 1666 1667 .keywords: SNES, nonlinear, set, function, norm, type 1668 1669 .seealso: SNESSetNormType(), SNESComputeFunction(), VecNorm(), SNESSetFunction(), SNESSetInitialFunction(), SNESNormType 1670 @*/ 1671 PetscErrorCode SNESGetNormType(SNES snes, SNESNormType *normtype) 1672 { 1673 PetscFunctionBegin; 1674 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1675 *normtype = snes->normtype; 1676 PetscFunctionReturn(0); 1677 } 1678 1679 /*MC 1680 SNESGSFunction - function used to convey a Gauss-Seidel sweep on the nonlinear function 1681 1682 Synopsis: 1683 #include "petscsnes.h" 1684 $ SNESGSFunction(SNES snes,Vec x,Vec b,void *ctx); 1685 1686 + X - solution vector 1687 . B - RHS vector 1688 - ctx - optional user-defined Gauss-Seidel context 1689 1690 Level: intermediate 1691 1692 .seealso: SNESSetGS(), SNESGetGS() 1693 M*/ 1694 1695 #undef __FUNCT__ 1696 #define __FUNCT__ "SNESSetGS" 1697 /*@C 1698 SNESSetGS - Sets the user nonlinear Gauss-Seidel routine for 1699 use with composed nonlinear solvers. 1700 1701 Input Parameters: 1702 + snes - the SNES context 1703 . SNESGSFunction - function evaluation routine 1704 - ctx - [optional] user-defined context for private data for the 1705 smoother evaluation routine (may be NULL) 1706 1707 Notes: 1708 The GS routines are used by the composed nonlinear solver to generate 1709 a problem appropriate update to the solution, particularly FAS. 1710 1711 Level: intermediate 1712 1713 .keywords: SNES, nonlinear, set, Gauss-Seidel 1714 1715 .seealso: SNESGetFunction(), SNESComputeGS() 1716 @*/ 1717 PetscErrorCode SNESSetGS(SNES snes,PetscErrorCode (*SNESGSFunction)(SNES,Vec,Vec,void*),void *ctx) 1718 { 1719 PetscErrorCode ierr; 1720 DM dm; 1721 1722 PetscFunctionBegin; 1723 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1724 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 1725 ierr = DMSNESSetGS(dm,SNESGSFunction,ctx);CHKERRQ(ierr); 1726 PetscFunctionReturn(0); 1727 } 1728 1729 #undef __FUNCT__ 1730 #define __FUNCT__ "SNESPicardComputeFunction" 1731 PetscErrorCode SNESPicardComputeFunction(SNES snes,Vec x,Vec f,void *ctx) 1732 { 1733 PetscErrorCode ierr; 1734 DM dm; 1735 DMSNES sdm; 1736 1737 PetscFunctionBegin; 1738 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 1739 ierr = DMGetDMSNES(dm,&sdm);CHKERRQ(ierr); 1740 /* A(x)*x - b(x) */ 1741 if (sdm->ops->computepfunction) { 1742 ierr = (*sdm->ops->computepfunction)(snes,x,f,sdm->pctx);CHKERRQ(ierr); 1743 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE, "Must call SNESSetPicard() to provide Picard function."); 1744 1745 if (sdm->ops->computepjacobian) { 1746 ierr = (*sdm->ops->computepjacobian)(snes,x,&snes->jacobian,&snes->jacobian_pre,&snes->matstruct,sdm->pctx);CHKERRQ(ierr); 1747 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE, "Must call SNESSetPicard() to provide Picard matrix."); 1748 ierr = VecScale(f,-1.0);CHKERRQ(ierr); 1749 ierr = MatMultAdd(snes->jacobian,x,f,f);CHKERRQ(ierr); 1750 PetscFunctionReturn(0); 1751 } 1752 1753 #undef __FUNCT__ 1754 #define __FUNCT__ "SNESPicardComputeJacobian" 1755 PetscErrorCode SNESPicardComputeJacobian(SNES snes,Vec x1,Mat *J,Mat *B,MatStructure *flag,void *ctx) 1756 { 1757 PetscFunctionBegin; 1758 /* the jacobian matrix should be pre-filled in SNESPicardComputeFunction */ 1759 *flag = snes->matstruct; 1760 PetscFunctionReturn(0); 1761 } 1762 1763 #undef __FUNCT__ 1764 #define __FUNCT__ "SNESSetPicard" 1765 /*@C 1766 SNESSetPicard - Use SNES to solve the semilinear-system A(x) x = b(x) via a Picard type iteration (Picard linearization) 1767 1768 Logically Collective on SNES 1769 1770 Input Parameters: 1771 + snes - the SNES context 1772 . r - vector to store function value 1773 . SNESFunction - function evaluation routine 1774 . Amat - matrix with which A(x) x - b(x) is to be computed 1775 . Pmat - matrix from which preconditioner is computed (usually the same as Amat) 1776 . SNESJacobianFunction - function to compute matrix value 1777 - ctx - [optional] user-defined context for private data for the 1778 function evaluation routine (may be NULL) 1779 1780 Notes: 1781 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 1782 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. 1783 1784 One can call SNESSetPicard() or SNESSetFunction() (and possibly SNESSetJacobian()) but cannot call both 1785 1786 $ 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} 1787 $ Note that when an exact solver is used this corresponds to the "classic" Picard A(x^{n}) x^{n+1} = b(x^{n}) iteration. 1788 1789 Run with -snes_mf_operator to solve the system with Newton's method using A(x^{n}) to construct the preconditioner. 1790 1791 We implement the defect correction form of the Picard iteration because it converges much more generally when inexact linear solvers are used then 1792 the direct Picard iteration A(x^n) x^{n+1} = b(x^n) 1793 1794 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 1795 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 1796 different please contact us at petsc-dev@mcs.anl.gov and we'll have an entirely new argument :-). 1797 1798 Level: intermediate 1799 1800 .keywords: SNES, nonlinear, set, function 1801 1802 .seealso: SNESGetFunction(), SNESSetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESGetPicard(), SNESLineSearchPreCheckPicard() 1803 @*/ 1804 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) 1805 { 1806 PetscErrorCode ierr; 1807 DM dm; 1808 1809 PetscFunctionBegin; 1810 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1811 ierr = SNESGetDM(snes, &dm);CHKERRQ(ierr); 1812 ierr = DMSNESSetPicard(dm,SNESFunction,SNESJacobianFunction,ctx);CHKERRQ(ierr); 1813 ierr = SNESSetFunction(snes,r,SNESPicardComputeFunction,ctx);CHKERRQ(ierr); 1814 ierr = SNESSetJacobian(snes,Amat,Pmat,SNESPicardComputeJacobian,ctx);CHKERRQ(ierr); 1815 PetscFunctionReturn(0); 1816 } 1817 1818 #undef __FUNCT__ 1819 #define __FUNCT__ "SNESGetPicard" 1820 /*@C 1821 SNESGetPicard - Returns the context for the Picard iteration 1822 1823 Not Collective, but Vec is parallel if SNES is parallel. Collective if Vec is requested, but has not been created yet. 1824 1825 Input Parameter: 1826 . snes - the SNES context 1827 1828 Output Parameter: 1829 + r - the function (or NULL) 1830 . SNESFunction - the function (or NULL) 1831 . Amat - the matrix used to defined the operation A(x) x - b(x) (or NULL) 1832 . Pmat - the matrix from which the preconditioner will be constructed (or NULL) 1833 . SNESJacobianFunction - the function for matrix evaluation (or NULL) 1834 - ctx - the function context (or NULL) 1835 1836 Level: advanced 1837 1838 .keywords: SNES, nonlinear, get, function 1839 1840 .seealso: SNESSetPicard(), SNESGetFunction(), SNESGetJacobian(), SNESGetDM(), SNESFunction, SNESJacobianFunction 1841 @*/ 1842 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) 1843 { 1844 PetscErrorCode ierr; 1845 DM dm; 1846 1847 PetscFunctionBegin; 1848 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1849 ierr = SNESGetFunction(snes,r,NULL,NULL);CHKERRQ(ierr); 1850 ierr = SNESGetJacobian(snes,Amat,Pmat,NULL,NULL);CHKERRQ(ierr); 1851 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 1852 ierr = DMSNESGetPicard(dm,SNESFunction,SNESJacobianFunction,ctx);CHKERRQ(ierr); 1853 PetscFunctionReturn(0); 1854 } 1855 1856 #undef __FUNCT__ 1857 #define __FUNCT__ "SNESSetComputeInitialGuess" 1858 /*@C 1859 SNESSetComputeInitialGuess - Sets a routine used to compute an initial guess for the problem 1860 1861 Logically Collective on SNES 1862 1863 Input Parameters: 1864 + snes - the SNES context 1865 . func - function evaluation routine 1866 - ctx - [optional] user-defined context for private data for the 1867 function evaluation routine (may be NULL) 1868 1869 Calling sequence of func: 1870 $ func (SNES snes,Vec x,void *ctx); 1871 1872 . f - function vector 1873 - ctx - optional user-defined function context 1874 1875 Level: intermediate 1876 1877 .keywords: SNES, nonlinear, set, function 1878 1879 .seealso: SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian() 1880 @*/ 1881 PetscErrorCode SNESSetComputeInitialGuess(SNES snes,PetscErrorCode (*func)(SNES,Vec,void*),void *ctx) 1882 { 1883 PetscFunctionBegin; 1884 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1885 if (func) snes->ops->computeinitialguess = func; 1886 if (ctx) snes->initialguessP = ctx; 1887 PetscFunctionReturn(0); 1888 } 1889 1890 /* --------------------------------------------------------------- */ 1891 #undef __FUNCT__ 1892 #define __FUNCT__ "SNESGetRhs" 1893 /*@C 1894 SNESGetRhs - Gets the vector for solving F(x) = rhs. If rhs is not set 1895 it assumes a zero right hand side. 1896 1897 Logically Collective on SNES 1898 1899 Input Parameter: 1900 . snes - the SNES context 1901 1902 Output Parameter: 1903 . rhs - the right hand side vector or NULL if the right hand side vector is null 1904 1905 Level: intermediate 1906 1907 .keywords: SNES, nonlinear, get, function, right hand side 1908 1909 .seealso: SNESGetSolution(), SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESSetFunction() 1910 @*/ 1911 PetscErrorCode SNESGetRhs(SNES snes,Vec *rhs) 1912 { 1913 PetscFunctionBegin; 1914 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1915 PetscValidPointer(rhs,2); 1916 *rhs = snes->vec_rhs; 1917 PetscFunctionReturn(0); 1918 } 1919 1920 #undef __FUNCT__ 1921 #define __FUNCT__ "SNESComputeFunction" 1922 /*@ 1923 SNESComputeFunction - Calls the function that has been set with SNESSetFunction(). 1924 1925 Collective on SNES 1926 1927 Input Parameters: 1928 + snes - the SNES context 1929 - x - input vector 1930 1931 Output Parameter: 1932 . y - function vector, as set by SNESSetFunction() 1933 1934 Notes: 1935 SNESComputeFunction() is typically used within nonlinear solvers 1936 implementations, so most users would not generally call this routine 1937 themselves. 1938 1939 Level: developer 1940 1941 .keywords: SNES, nonlinear, compute, function 1942 1943 .seealso: SNESSetFunction(), SNESGetFunction() 1944 @*/ 1945 PetscErrorCode SNESComputeFunction(SNES snes,Vec x,Vec y) 1946 { 1947 PetscErrorCode ierr; 1948 DM dm; 1949 DMSNES sdm; 1950 1951 PetscFunctionBegin; 1952 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1953 PetscValidHeaderSpecific(x,VEC_CLASSID,2); 1954 PetscValidHeaderSpecific(y,VEC_CLASSID,3); 1955 PetscCheckSameComm(snes,1,x,2); 1956 PetscCheckSameComm(snes,1,y,3); 1957 ierr = VecValidValues(x,2,PETSC_TRUE);CHKERRQ(ierr); 1958 1959 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 1960 ierr = DMGetDMSNES(dm,&sdm);CHKERRQ(ierr); 1961 if (snes->pc && snes->pcside == PC_LEFT) { 1962 ierr = VecCopy(x,y);CHKERRQ(ierr); 1963 ierr = PetscLogEventBegin(SNES_NPCSolve,snes->pc,x,y,0);CHKERRQ(ierr); 1964 ierr = SNESSolve(snes->pc,snes->vec_rhs,y);CHKERRQ(ierr); 1965 ierr = PetscLogEventEnd(SNES_NPCSolve,snes->pc,x,y,0);CHKERRQ(ierr); 1966 ierr = VecAYPX(y,-1.0,x);CHKERRQ(ierr); 1967 } else if (sdm->ops->computefunction) { 1968 ierr = PetscLogEventBegin(SNES_FunctionEval,snes,x,y,0);CHKERRQ(ierr); 1969 PetscStackPush("SNES user function"); 1970 ierr = (*sdm->ops->computefunction)(snes,x,y,sdm->functionctx);CHKERRQ(ierr); 1971 PetscStackPop; 1972 ierr = PetscLogEventEnd(SNES_FunctionEval,snes,x,y,0);CHKERRQ(ierr); 1973 } else if (snes->vec_rhs) { 1974 ierr = MatMult(snes->jacobian, x, y);CHKERRQ(ierr); 1975 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE, "Must call SNESSetFunction() or SNESSetDM() before SNESComputeFunction(), likely called from SNESSolve()."); 1976 if (snes->vec_rhs) { 1977 ierr = VecAXPY(y,-1.0,snes->vec_rhs);CHKERRQ(ierr); 1978 } 1979 snes->nfuncs++; 1980 ierr = VecValidValues(y,3,PETSC_FALSE);CHKERRQ(ierr); 1981 PetscFunctionReturn(0); 1982 } 1983 1984 #undef __FUNCT__ 1985 #define __FUNCT__ "SNESComputeGS" 1986 /*@ 1987 SNESComputeGS - Calls the Gauss-Seidel function that has been set with SNESSetGS(). 1988 1989 Collective on SNES 1990 1991 Input Parameters: 1992 + snes - the SNES context 1993 . x - input vector 1994 - b - rhs vector 1995 1996 Output Parameter: 1997 . x - new solution vector 1998 1999 Notes: 2000 SNESComputeGS() is typically used within composed nonlinear solver 2001 implementations, so most users would not generally call this routine 2002 themselves. 2003 2004 Level: developer 2005 2006 .keywords: SNES, nonlinear, compute, function 2007 2008 .seealso: SNESSetGS(), SNESComputeFunction() 2009 @*/ 2010 PetscErrorCode SNESComputeGS(SNES snes,Vec b,Vec x) 2011 { 2012 PetscErrorCode ierr; 2013 DM dm; 2014 DMSNES sdm; 2015 2016 PetscFunctionBegin; 2017 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2018 PetscValidHeaderSpecific(x,VEC_CLASSID,2); 2019 if (b) PetscValidHeaderSpecific(b,VEC_CLASSID,3); 2020 PetscCheckSameComm(snes,1,x,2); 2021 if (b) PetscCheckSameComm(snes,1,b,3); 2022 if (b) {ierr = VecValidValues(b,2,PETSC_TRUE);CHKERRQ(ierr);} 2023 ierr = PetscLogEventBegin(SNES_GSEval,snes,x,b,0);CHKERRQ(ierr); 2024 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 2025 ierr = DMGetDMSNES(dm,&sdm);CHKERRQ(ierr); 2026 if (sdm->ops->computegs) { 2027 PetscStackPush("SNES user GS"); 2028 ierr = (*sdm->ops->computegs)(snes,x,b,sdm->gsctx);CHKERRQ(ierr); 2029 PetscStackPop; 2030 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE, "Must call SNESSetGS() before SNESComputeGS(), likely called from SNESSolve()."); 2031 ierr = PetscLogEventEnd(SNES_GSEval,snes,x,b,0);CHKERRQ(ierr); 2032 ierr = VecValidValues(x,3,PETSC_FALSE);CHKERRQ(ierr); 2033 PetscFunctionReturn(0); 2034 } 2035 2036 #undef __FUNCT__ 2037 #define __FUNCT__ "SNESComputeJacobian" 2038 /*@ 2039 SNESComputeJacobian - Computes the Jacobian matrix that has been set with SNESSetJacobian(). 2040 2041 Collective on SNES and Mat 2042 2043 Input Parameters: 2044 + snes - the SNES context 2045 - x - input vector 2046 2047 Output Parameters: 2048 + A - Jacobian matrix 2049 . B - optional preconditioning matrix 2050 - flag - flag indicating matrix structure (one of, SAME_NONZERO_PATTERN,DIFFERENT_NONZERO_PATTERN,SAME_PRECONDITIONER) 2051 2052 Options Database Keys: 2053 + -snes_lag_preconditioner <lag> 2054 . -snes_lag_jacobian <lag> 2055 . -snes_compare_explicit - Compare the computed Jacobian to the finite difference Jacobian and output the differences 2056 . -snes_compare_explicit_draw - Compare the computed Jacobian to the finite difference Jacobian and draw the result 2057 . -snes_compare_explicit_contour - Compare the computed Jacobian to the finite difference Jacobian and draw a contour plot with the result 2058 . -snes_compare_operator - Make the comparison options above use the operator instead of the preconditioning matrix 2059 . -snes_compare_coloring - Compute the finite differece Jacobian using coloring and display norms of difference 2060 . -snes_compare_coloring_display - Compute the finite differece Jacobian using coloring and display verbose differences 2061 . -snes_compare_coloring_threshold - Display only those matrix entries that differ by more than a given threshold 2062 . -snes_compare_coloring_threshold_atol - Absolute tolerance for difference in matrix entries to be displayed by -snes_compare_coloring_threshold 2063 . -snes_compare_coloring_threshold_rtol - Relative tolerance for difference in matrix entries to be displayed by -snes_compare_coloring_threshold 2064 . -snes_compare_coloring_draw - Compute the finite differece Jacobian using coloring and draw differences 2065 - -snes_compare_coloring_draw_contour - Compute the finite differece Jacobian using coloring and show contours of matrices and differences 2066 2067 2068 Notes: 2069 Most users should not need to explicitly call this routine, as it 2070 is used internally within the nonlinear solvers. 2071 2072 See KSPSetOperators() for important information about setting the 2073 flag parameter. 2074 2075 Level: developer 2076 2077 .keywords: SNES, compute, Jacobian, matrix 2078 2079 .seealso: SNESSetJacobian(), KSPSetOperators(), MatStructure, SNESSetLagPreconditioner(), SNESSetLagJacobian() 2080 @*/ 2081 PetscErrorCode SNESComputeJacobian(SNES snes,Vec X,Mat *A,Mat *B,MatStructure *flg) 2082 { 2083 PetscErrorCode ierr; 2084 PetscBool flag; 2085 DM dm; 2086 DMSNES sdm; 2087 2088 PetscFunctionBegin; 2089 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2090 PetscValidHeaderSpecific(X,VEC_CLASSID,2); 2091 PetscValidPointer(flg,5); 2092 PetscCheckSameComm(snes,1,X,2); 2093 ierr = VecValidValues(X,2,PETSC_TRUE);CHKERRQ(ierr); 2094 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 2095 ierr = DMGetDMSNES(dm,&sdm);CHKERRQ(ierr); 2096 2097 if (!sdm->ops->computejacobian) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_USER,"Must call SNESSetJacobian(), DMSNESSetJacobian(), DMDASNESSetJacobianLocal(), etc"); 2098 2099 /* make sure that MatAssemblyBegin/End() is called on A matrix if it is matrix free */ 2100 2101 if (snes->lagjacobian == -2) { 2102 snes->lagjacobian = -1; 2103 2104 ierr = PetscInfo(snes,"Recomputing Jacobian/preconditioner because lag is -2 (means compute Jacobian, but then never again) \n");CHKERRQ(ierr); 2105 } else if (snes->lagjacobian == -1) { 2106 *flg = SAME_PRECONDITIONER; 2107 ierr = PetscInfo(snes,"Reusing Jacobian/preconditioner because lag is -1\n");CHKERRQ(ierr); 2108 ierr = PetscObjectTypeCompare((PetscObject)*A,MATMFFD,&flag);CHKERRQ(ierr); 2109 if (flag) { 2110 ierr = MatAssemblyBegin(*A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2111 ierr = MatAssemblyEnd(*A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2112 } 2113 PetscFunctionReturn(0); 2114 } else if (snes->lagjacobian > 1 && snes->iter % snes->lagjacobian) { 2115 *flg = SAME_PRECONDITIONER; 2116 ierr = PetscInfo2(snes,"Reusing Jacobian/preconditioner because lag is %D and SNES iteration is %D\n",snes->lagjacobian,snes->iter);CHKERRQ(ierr); 2117 ierr = PetscObjectTypeCompare((PetscObject)*A,MATMFFD,&flag);CHKERRQ(ierr); 2118 if (flag) { 2119 ierr = MatAssemblyBegin(*A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2120 ierr = MatAssemblyEnd(*A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2121 } 2122 PetscFunctionReturn(0); 2123 } 2124 if (snes->pc && snes->pcside == PC_LEFT) { 2125 ierr = MatAssemblyBegin(*A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2126 ierr = MatAssemblyEnd(*A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2127 PetscFunctionReturn(0); 2128 } 2129 2130 *flg = DIFFERENT_NONZERO_PATTERN; 2131 ierr = PetscLogEventBegin(SNES_JacobianEval,snes,X,*A,*B);CHKERRQ(ierr); 2132 2133 PetscStackPush("SNES user Jacobian function"); 2134 ierr = (*sdm->ops->computejacobian)(snes,X,A,B,flg,sdm->jacobianctx);CHKERRQ(ierr); 2135 PetscStackPop; 2136 2137 ierr = PetscLogEventEnd(SNES_JacobianEval,snes,X,*A,*B);CHKERRQ(ierr); 2138 2139 if (snes->lagpreconditioner == -2) { 2140 ierr = PetscInfo(snes,"Rebuilding preconditioner exactly once since lag is -2\n");CHKERRQ(ierr); 2141 2142 snes->lagpreconditioner = -1; 2143 } else if (snes->lagpreconditioner == -1) { 2144 *flg = SAME_PRECONDITIONER; 2145 ierr = PetscInfo(snes,"Reusing preconditioner because lag is -1\n");CHKERRQ(ierr); 2146 } else if (snes->lagpreconditioner > 1 && snes->iter % snes->lagpreconditioner) { 2147 *flg = SAME_PRECONDITIONER; 2148 ierr = PetscInfo2(snes,"Reusing preconditioner because lag is %D and SNES iteration is %D\n",snes->lagpreconditioner,snes->iter);CHKERRQ(ierr); 2149 } 2150 2151 /* make sure user returned a correct Jacobian and preconditioner */ 2152 /* PetscValidHeaderSpecific(*A,MAT_CLASSID,3); 2153 PetscValidHeaderSpecific(*B,MAT_CLASSID,4); */ 2154 { 2155 PetscBool flag = PETSC_FALSE,flag_draw = PETSC_FALSE,flag_contour = PETSC_FALSE,flag_operator = PETSC_FALSE; 2156 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_explicit",&flag,NULL);CHKERRQ(ierr); 2157 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_explicit_draw",&flag_draw,NULL);CHKERRQ(ierr); 2158 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_explicit_draw_contour",&flag_contour,NULL);CHKERRQ(ierr); 2159 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_operator",&flag_operator,NULL);CHKERRQ(ierr); 2160 if (flag || flag_draw || flag_contour) { 2161 Mat Bexp_mine = NULL,Bexp,FDexp; 2162 MatStructure mstruct; 2163 PetscViewer vdraw,vstdout; 2164 PetscBool flg; 2165 if (flag_operator) { 2166 ierr = MatComputeExplicitOperator(*A,&Bexp_mine);CHKERRQ(ierr); 2167 Bexp = Bexp_mine; 2168 } else { 2169 /* See if the preconditioning matrix can be viewed and added directly */ 2170 ierr = PetscObjectTypeCompareAny((PetscObject)*B,&flg,MATSEQAIJ,MATMPIAIJ,MATSEQDENSE,MATMPIDENSE,MATSEQBAIJ,MATMPIBAIJ,MATSEQSBAIJ,MATMPIBAIJ,"");CHKERRQ(ierr); 2171 if (flg) Bexp = *B; 2172 else { 2173 /* If the "preconditioning" matrix is itself MATSHELL or some other type without direct support */ 2174 ierr = MatComputeExplicitOperator(*B,&Bexp_mine);CHKERRQ(ierr); 2175 Bexp = Bexp_mine; 2176 } 2177 } 2178 ierr = MatConvert(Bexp,MATSAME,MAT_INITIAL_MATRIX,&FDexp);CHKERRQ(ierr); 2179 ierr = SNESComputeJacobianDefault(snes,X,&FDexp,&FDexp,&mstruct,NULL);CHKERRQ(ierr); 2180 ierr = PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)snes),&vstdout);CHKERRQ(ierr); 2181 if (flag_draw || flag_contour) { 2182 ierr = PetscViewerDrawOpen(PetscObjectComm((PetscObject)snes),0,"Explicit Jacobians",PETSC_DECIDE,PETSC_DECIDE,300,300,&vdraw);CHKERRQ(ierr); 2183 if (flag_contour) {ierr = PetscViewerPushFormat(vdraw,PETSC_VIEWER_DRAW_CONTOUR);CHKERRQ(ierr);} 2184 } else vdraw = NULL; 2185 ierr = PetscViewerASCIIPrintf(vstdout,"Explicit %s\n",flag_operator ? "Jacobian" : "preconditioning Jacobian");CHKERRQ(ierr); 2186 if (flag) {ierr = MatView(Bexp,vstdout);CHKERRQ(ierr);} 2187 if (vdraw) {ierr = MatView(Bexp,vdraw);CHKERRQ(ierr);} 2188 ierr = PetscViewerASCIIPrintf(vstdout,"Finite difference Jacobian\n");CHKERRQ(ierr); 2189 if (flag) {ierr = MatView(FDexp,vstdout);CHKERRQ(ierr);} 2190 if (vdraw) {ierr = MatView(FDexp,vdraw);CHKERRQ(ierr);} 2191 ierr = MatAYPX(FDexp,-1.0,Bexp,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 2192 ierr = PetscViewerASCIIPrintf(vstdout,"User-provided matrix minus finite difference Jacobian\n");CHKERRQ(ierr); 2193 if (flag) {ierr = MatView(FDexp,vstdout);CHKERRQ(ierr);} 2194 if (vdraw) { /* Always use contour for the difference */ 2195 ierr = PetscViewerPushFormat(vdraw,PETSC_VIEWER_DRAW_CONTOUR);CHKERRQ(ierr); 2196 ierr = MatView(FDexp,vdraw);CHKERRQ(ierr); 2197 ierr = PetscViewerPopFormat(vdraw);CHKERRQ(ierr); 2198 } 2199 if (flag_contour) {ierr = PetscViewerPopFormat(vdraw);CHKERRQ(ierr);} 2200 ierr = PetscViewerDestroy(&vdraw);CHKERRQ(ierr); 2201 ierr = MatDestroy(&Bexp_mine);CHKERRQ(ierr); 2202 ierr = MatDestroy(&FDexp);CHKERRQ(ierr); 2203 } 2204 } 2205 { 2206 PetscBool flag = PETSC_FALSE,flag_display = PETSC_FALSE,flag_draw = PETSC_FALSE,flag_contour = PETSC_FALSE,flag_threshold = PETSC_FALSE; 2207 PetscReal threshold_atol = PETSC_SQRT_MACHINE_EPSILON,threshold_rtol = 10*PETSC_SQRT_MACHINE_EPSILON; 2208 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_coloring",&flag,NULL);CHKERRQ(ierr); 2209 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_coloring_display",&flag_display,NULL);CHKERRQ(ierr); 2210 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_coloring_draw",&flag_draw,NULL);CHKERRQ(ierr); 2211 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_coloring_draw_contour",&flag_contour,NULL);CHKERRQ(ierr); 2212 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_coloring_threshold",&flag_threshold,NULL);CHKERRQ(ierr); 2213 ierr = PetscOptionsGetReal(((PetscObject)snes)->prefix,"-snes_compare_coloring_threshold_rtol",&threshold_rtol,NULL);CHKERRQ(ierr); 2214 ierr = PetscOptionsGetReal(((PetscObject)snes)->prefix,"-snes_compare_coloring_threshold_atol",&threshold_atol,NULL);CHKERRQ(ierr); 2215 if (flag || flag_display || flag_draw || flag_contour || flag_threshold) { 2216 Mat Bfd; 2217 MatStructure mstruct; 2218 PetscViewer vdraw,vstdout; 2219 ISColoring iscoloring; 2220 MatFDColoring matfdcoloring; 2221 PetscErrorCode (*func)(SNES,Vec,Vec,void*); 2222 void *funcctx; 2223 PetscReal norm1,norm2,normmax; 2224 2225 ierr = MatDuplicate(*B,MAT_DO_NOT_COPY_VALUES,&Bfd);CHKERRQ(ierr); 2226 ierr = MatGetColoring(Bfd,MATCOLORINGSL,&iscoloring);CHKERRQ(ierr); 2227 ierr = MatFDColoringCreate(Bfd,iscoloring,&matfdcoloring);CHKERRQ(ierr); 2228 ierr = ISColoringDestroy(&iscoloring);CHKERRQ(ierr); 2229 2230 /* This method of getting the function is currently unreliable since it doesn't work for DM local functions. */ 2231 ierr = SNESGetFunction(snes,NULL,&func,&funcctx);CHKERRQ(ierr); 2232 ierr = MatFDColoringSetFunction(matfdcoloring,(PetscErrorCode (*)(void))func,funcctx);CHKERRQ(ierr); 2233 ierr = PetscObjectSetOptionsPrefix((PetscObject)matfdcoloring,((PetscObject)snes)->prefix);CHKERRQ(ierr); 2234 ierr = PetscObjectAppendOptionsPrefix((PetscObject)matfdcoloring,"coloring_");CHKERRQ(ierr); 2235 ierr = MatFDColoringSetFromOptions(matfdcoloring);CHKERRQ(ierr); 2236 ierr = MatFDColoringApply(Bfd,matfdcoloring,X,&mstruct,snes);CHKERRQ(ierr); 2237 ierr = MatFDColoringDestroy(&matfdcoloring);CHKERRQ(ierr); 2238 2239 ierr = PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)snes),&vstdout);CHKERRQ(ierr); 2240 if (flag_draw || flag_contour) { 2241 ierr = PetscViewerDrawOpen(PetscObjectComm((PetscObject)snes),0,"Colored Jacobians",PETSC_DECIDE,PETSC_DECIDE,300,300,&vdraw);CHKERRQ(ierr); 2242 if (flag_contour) {ierr = PetscViewerPushFormat(vdraw,PETSC_VIEWER_DRAW_CONTOUR);CHKERRQ(ierr);} 2243 } else vdraw = NULL; 2244 ierr = PetscViewerASCIIPrintf(vstdout,"Explicit preconditioning Jacobian\n");CHKERRQ(ierr); 2245 if (flag_display) {ierr = MatView(*B,vstdout);CHKERRQ(ierr);} 2246 if (vdraw) {ierr = MatView(*B,vdraw);CHKERRQ(ierr);} 2247 ierr = PetscViewerASCIIPrintf(vstdout,"Colored Finite difference Jacobian\n");CHKERRQ(ierr); 2248 if (flag_display) {ierr = MatView(Bfd,vstdout);CHKERRQ(ierr);} 2249 if (vdraw) {ierr = MatView(Bfd,vdraw);CHKERRQ(ierr);} 2250 ierr = MatAYPX(Bfd,-1.0,*B,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 2251 ierr = MatNorm(Bfd,NORM_1,&norm1);CHKERRQ(ierr); 2252 ierr = MatNorm(Bfd,NORM_FROBENIUS,&norm2);CHKERRQ(ierr); 2253 ierr = MatNorm(Bfd,NORM_MAX,&normmax);CHKERRQ(ierr); 2254 ierr = PetscViewerASCIIPrintf(vstdout,"User-provided matrix minus finite difference Jacobian, norm1=%G normFrob=%G normmax=%G\n",norm1,norm2,normmax);CHKERRQ(ierr); 2255 if (flag_display) {ierr = MatView(Bfd,vstdout);CHKERRQ(ierr);} 2256 if (vdraw) { /* Always use contour for the difference */ 2257 ierr = PetscViewerPushFormat(vdraw,PETSC_VIEWER_DRAW_CONTOUR);CHKERRQ(ierr); 2258 ierr = MatView(Bfd,vdraw);CHKERRQ(ierr); 2259 ierr = PetscViewerPopFormat(vdraw);CHKERRQ(ierr); 2260 } 2261 if (flag_contour) {ierr = PetscViewerPopFormat(vdraw);CHKERRQ(ierr);} 2262 2263 if (flag_threshold) { 2264 PetscInt bs,rstart,rend,i; 2265 ierr = MatGetBlockSize(*B,&bs);CHKERRQ(ierr); 2266 ierr = MatGetOwnershipRange(*B,&rstart,&rend);CHKERRQ(ierr); 2267 for (i=rstart; i<rend; i++) { 2268 const PetscScalar *ba,*ca; 2269 const PetscInt *bj,*cj; 2270 PetscInt bn,cn,j,maxentrycol = -1,maxdiffcol = -1,maxrdiffcol = -1; 2271 PetscReal maxentry = 0,maxdiff = 0,maxrdiff = 0; 2272 ierr = MatGetRow(*B,i,&bn,&bj,&ba);CHKERRQ(ierr); 2273 ierr = MatGetRow(Bfd,i,&cn,&cj,&ca);CHKERRQ(ierr); 2274 if (bn != cn) SETERRQ(((PetscObject)*A)->comm,PETSC_ERR_PLIB,"Unexpected different nonzero pattern in -snes_compare_coloring_threshold"); 2275 for (j=0; j<bn; j++) { 2276 PetscReal rdiff = PetscAbsScalar(ca[j]) / (threshold_atol + threshold_rtol*PetscAbsScalar(ba[j])); 2277 if (PetscAbsScalar(ba[j]) > PetscAbs(maxentry)) { 2278 maxentrycol = bj[j]; 2279 maxentry = PetscRealPart(ba[j]); 2280 } 2281 if (PetscAbsScalar(ca[j]) > PetscAbs(maxdiff)) { 2282 maxdiffcol = bj[j]; 2283 maxdiff = PetscRealPart(ca[j]); 2284 } 2285 if (rdiff > maxrdiff) { 2286 maxrdiffcol = bj[j]; 2287 maxrdiff = rdiff; 2288 } 2289 } 2290 if (maxrdiff > 1) { 2291 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); 2292 for (j=0; j<bn; j++) { 2293 PetscReal rdiff; 2294 rdiff = PetscAbsScalar(ca[j]) / (threshold_atol + threshold_rtol*PetscAbsScalar(ba[j])); 2295 if (rdiff > 1) { 2296 ierr = PetscViewerASCIIPrintf(vstdout," (%D,%G:%G)",bj[j],PetscRealPart(ba[j]),PetscRealPart(ca[j]));CHKERRQ(ierr); 2297 } 2298 } 2299 ierr = PetscViewerASCIIPrintf(vstdout,"\n",i,maxentry,maxdiff,maxrdiff);CHKERRQ(ierr); 2300 } 2301 ierr = MatRestoreRow(*B,i,&bn,&bj,&ba);CHKERRQ(ierr); 2302 ierr = MatRestoreRow(Bfd,i,&cn,&cj,&ca);CHKERRQ(ierr); 2303 } 2304 } 2305 ierr = PetscViewerDestroy(&vdraw);CHKERRQ(ierr); 2306 ierr = MatDestroy(&Bfd);CHKERRQ(ierr); 2307 } 2308 } 2309 PetscFunctionReturn(0); 2310 } 2311 2312 /*MC 2313 SNESJacobianFunction - function used to convey the nonlinear Jacobian of the function to be solved by SNES 2314 2315 Synopsis: 2316 #include "petscsnes.h" 2317 $ SNESJacobianFunction(SNES snes,Vec x,Mat *Amat,Mat *Pmat,int *flag,void *ctx); 2318 2319 + x - input vector 2320 . Amat - the matrix that defines the (approximate) Jacobian 2321 . Pmat - the matrix to be used in constructing the preconditioner, usually the same as Amat. 2322 . flag - flag indicating information about the preconditioner matrix 2323 structure (same as flag in KSPSetOperators()), one of SAME_NONZERO_PATTERN,DIFFERENT_NONZERO_PATTERN,SAME_PRECONDITIONER 2324 - ctx - [optional] user-defined Jacobian context 2325 2326 Level: intermediate 2327 2328 .seealso: SNESSetFunction(), SNESGetFunction(), SNESSetJacobian(), SNESGetJacobian() 2329 M*/ 2330 2331 #undef __FUNCT__ 2332 #define __FUNCT__ "SNESSetJacobian" 2333 /*@C 2334 SNESSetJacobian - Sets the function to compute Jacobian as well as the 2335 location to store the matrix. 2336 2337 Logically Collective on SNES and Mat 2338 2339 Input Parameters: 2340 + snes - the SNES context 2341 . Amat - the matrix that defines the (approximate) Jacobian 2342 . Pmat - the matrix to be used in constructing the preconditioner, usually the same as Amat. 2343 . SNESJacobianFunction - Jacobian evaluation routine (if NULL then SNES retains any previously set value) 2344 - ctx - [optional] user-defined context for private data for the 2345 Jacobian evaluation routine (may be NULL) (if NULL then SNES retains any previously set value) 2346 2347 Notes: 2348 See KSPSetOperators() for important information about setting the flag 2349 output parameter in the routine func(). Be sure to read this information! 2350 2351 The routine func() takes Mat * as the matrix arguments rather than Mat. 2352 This allows the Jacobian evaluation routine to replace A and/or B with a 2353 completely new new matrix structure (not just different matrix elements) 2354 when appropriate, for instance, if the nonzero structure is changing 2355 throughout the global iterations. 2356 2357 If the Amat matrix and Pmat matrix are different you must call MatAssemblyBegin/End() on 2358 each matrix. 2359 2360 If using SNESComputeJacobianDefaultColor() to assemble a Jacobian, the ctx argument 2361 must be a MatFDColoring. 2362 2363 Other defect-correction schemes can be used by computing a different matrix in place of the Jacobian. One common 2364 example is to use the "Picard linearization" which only differentiates through the highest order parts of each term. 2365 2366 Level: beginner 2367 2368 .keywords: SNES, nonlinear, set, Jacobian, matrix 2369 2370 .seealso: KSPSetOperators(), SNESSetFunction(), MatMFFDComputeJacobian(), SNESComputeJacobianDefaultColor(), MatStructure, SNESJacobianFunction, SNESSetPicard() 2371 @*/ 2372 PetscErrorCode SNESSetJacobian(SNES snes,Mat Amat,Mat Pmat,PetscErrorCode (*SNESJacobianFunction)(SNES,Vec,Mat*,Mat*,MatStructure*,void*),void *ctx) 2373 { 2374 PetscErrorCode ierr; 2375 DM dm; 2376 2377 PetscFunctionBegin; 2378 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2379 if (Amat) PetscValidHeaderSpecific(Amat,MAT_CLASSID,2); 2380 if (Pmat) PetscValidHeaderSpecific(Pmat,MAT_CLASSID,3); 2381 if (Amat) PetscCheckSameComm(snes,1,Amat,2); 2382 if (Pmat) PetscCheckSameComm(snes,1,Pmat,3); 2383 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 2384 ierr = DMSNESSetJacobian(dm,SNESJacobianFunction,ctx);CHKERRQ(ierr); 2385 if (Amat) { 2386 ierr = PetscObjectReference((PetscObject)Amat);CHKERRQ(ierr); 2387 ierr = MatDestroy(&snes->jacobian);CHKERRQ(ierr); 2388 2389 snes->jacobian = Amat; 2390 } 2391 if (Pmat) { 2392 ierr = PetscObjectReference((PetscObject)Pmat);CHKERRQ(ierr); 2393 ierr = MatDestroy(&snes->jacobian_pre);CHKERRQ(ierr); 2394 2395 snes->jacobian_pre = Pmat; 2396 } 2397 PetscFunctionReturn(0); 2398 } 2399 2400 #undef __FUNCT__ 2401 #define __FUNCT__ "SNESGetJacobian" 2402 /*@C 2403 SNESGetJacobian - Returns the Jacobian matrix and optionally the user 2404 provided context for evaluating the Jacobian. 2405 2406 Not Collective, but Mat object will be parallel if SNES object is 2407 2408 Input Parameter: 2409 . snes - the nonlinear solver context 2410 2411 Output Parameters: 2412 + Amat - location to stash (approximate) Jacobian matrix (or NULL) 2413 . Pmat - location to stash matrix used to compute the preconditioner (or NULL) 2414 . SNESJacobianFunction - location to put Jacobian function (or NULL) 2415 - ctx - location to stash Jacobian ctx (or NULL) 2416 2417 Level: advanced 2418 2419 .seealso: SNESSetJacobian(), SNESComputeJacobian() 2420 @*/ 2421 PetscErrorCode SNESGetJacobian(SNES snes,Mat *Amat,Mat *Pmat,PetscErrorCode (**SNESJacobianFunction)(SNES,Vec,Mat*,Mat*,MatStructure*,void*),void **ctx) 2422 { 2423 PetscErrorCode ierr; 2424 DM dm; 2425 DMSNES sdm; 2426 2427 PetscFunctionBegin; 2428 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2429 if (Amat) *Amat = snes->jacobian; 2430 if (Pmat) *Pmat = snes->jacobian_pre; 2431 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 2432 ierr = DMGetDMSNES(dm,&sdm);CHKERRQ(ierr); 2433 if (SNESJacobianFunction) *SNESJacobianFunction = sdm->ops->computejacobian; 2434 if (ctx) *ctx = sdm->jacobianctx; 2435 PetscFunctionReturn(0); 2436 } 2437 2438 #undef __FUNCT__ 2439 #define __FUNCT__ "SNESSetUp" 2440 /*@ 2441 SNESSetUp - Sets up the internal data structures for the later use 2442 of a nonlinear solver. 2443 2444 Collective on SNES 2445 2446 Input Parameters: 2447 . snes - the SNES context 2448 2449 Notes: 2450 For basic use of the SNES solvers the user need not explicitly call 2451 SNESSetUp(), since these actions will automatically occur during 2452 the call to SNESSolve(). However, if one wishes to control this 2453 phase separately, SNESSetUp() should be called after SNESCreate() 2454 and optional routines of the form SNESSetXXX(), but before SNESSolve(). 2455 2456 Level: advanced 2457 2458 .keywords: SNES, nonlinear, setup 2459 2460 .seealso: SNESCreate(), SNESSolve(), SNESDestroy() 2461 @*/ 2462 PetscErrorCode SNESSetUp(SNES snes) 2463 { 2464 PetscErrorCode ierr; 2465 DM dm; 2466 DMSNES sdm; 2467 SNESLineSearch linesearch, pclinesearch; 2468 void *lsprectx,*lspostctx; 2469 PetscErrorCode (*precheck)(SNESLineSearch,Vec,Vec,PetscBool*,void*); 2470 PetscErrorCode (*postcheck)(SNESLineSearch,Vec,Vec,Vec,PetscBool*,PetscBool*,void*); 2471 PetscErrorCode (*func)(SNES,Vec,Vec,void*); 2472 Vec f,fpc; 2473 void *funcctx; 2474 PetscErrorCode (*jac)(SNES,Vec,Mat*,Mat*,MatStructure*,void*); 2475 void *jacctx,*appctx; 2476 Mat A,B; 2477 2478 PetscFunctionBegin; 2479 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2480 if (snes->setupcalled) PetscFunctionReturn(0); 2481 2482 if (!((PetscObject)snes)->type_name) { 2483 ierr = SNESSetType(snes,SNESNEWTONLS);CHKERRQ(ierr); 2484 } 2485 2486 ierr = SNESGetFunction(snes,&snes->vec_func,NULL,NULL);CHKERRQ(ierr); 2487 if (snes->vec_func == snes->vec_sol) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_IDN,"Solution vector cannot be function vector"); 2488 if (snes->vec_rhs == snes->vec_sol) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_IDN,"Solution vector cannot be right hand side vector"); 2489 2490 if (!snes->vec_sol_update /* && snes->vec_sol */) { 2491 ierr = VecDuplicate(snes->vec_sol,&snes->vec_sol_update);CHKERRQ(ierr); 2492 ierr = PetscLogObjectParent(snes,snes->vec_sol_update);CHKERRQ(ierr); 2493 } 2494 2495 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 2496 ierr = DMShellSetGlobalVector(snes->dm,snes->vec_sol);CHKERRQ(ierr); 2497 ierr = DMGetDMSNES(dm,&sdm);CHKERRQ(ierr); 2498 if (!sdm->ops->computefunction) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"Function never provided to SNES object"); 2499 if (!sdm->ops->computejacobian) { 2500 ierr = DMSNESSetJacobian(dm,SNESComputeJacobianDefaultColor,NULL);CHKERRQ(ierr); 2501 } 2502 if (!snes->vec_func) { 2503 ierr = DMCreateGlobalVector(dm,&snes->vec_func);CHKERRQ(ierr); 2504 } 2505 2506 if (!snes->ksp) { 2507 ierr = SNESGetKSP(snes, &snes->ksp);CHKERRQ(ierr); 2508 } 2509 2510 if (!snes->linesearch) { 2511 ierr = SNESGetSNESLineSearch(snes, &snes->linesearch);CHKERRQ(ierr); 2512 } 2513 2514 if (snes->pc && (snes->pcside == PC_LEFT)) { 2515 snes->mf = PETSC_TRUE; 2516 snes->mf_operator = PETSC_FALSE; 2517 } 2518 2519 if (snes->mf) { 2520 ierr = SNESSetUpMatrixFree_Private(snes, snes->mf_operator, snes->mf_version);CHKERRQ(ierr); 2521 } 2522 2523 if (snes->ops->usercompute && !snes->user) { 2524 ierr = (*snes->ops->usercompute)(snes,(void**)&snes->user);CHKERRQ(ierr); 2525 } 2526 2527 if (snes->pc) { 2528 /* copy the DM over */ 2529 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 2530 ierr = SNESSetDM(snes->pc,dm);CHKERRQ(ierr); 2531 2532 ierr = SNESGetFunction(snes,&f,&func,&funcctx);CHKERRQ(ierr); 2533 ierr = VecDuplicate(f,&fpc);CHKERRQ(ierr); 2534 ierr = SNESSetFunction(snes->pc,fpc,func,funcctx);CHKERRQ(ierr); 2535 ierr = SNESGetJacobian(snes,&A,&B,&jac,&jacctx);CHKERRQ(ierr); 2536 ierr = SNESSetJacobian(snes->pc,A,B,jac,jacctx);CHKERRQ(ierr); 2537 ierr = SNESGetApplicationContext(snes,&appctx);CHKERRQ(ierr); 2538 ierr = SNESSetApplicationContext(snes->pc,appctx);CHKERRQ(ierr); 2539 ierr = VecDestroy(&fpc);CHKERRQ(ierr); 2540 2541 /* copy the function pointers over */ 2542 ierr = PetscObjectCopyFortranFunctionPointers((PetscObject)snes,(PetscObject)snes->pc);CHKERRQ(ierr); 2543 2544 /* default to 1 iteration */ 2545 ierr = SNESSetTolerances(snes->pc,0.0,0.0,0.0,1,snes->pc->max_funcs);CHKERRQ(ierr); 2546 if (snes->pcside==PC_RIGHT) { 2547 ierr = SNESSetNormType(snes->pc,SNES_NORM_FINAL_ONLY);CHKERRQ(ierr); 2548 } else { 2549 ierr = SNESSetNormType(snes->pc,SNES_NORM_NONE);CHKERRQ(ierr); 2550 } 2551 ierr = SNESSetFromOptions(snes->pc);CHKERRQ(ierr); 2552 2553 /* copy the line search context over */ 2554 ierr = SNESGetSNESLineSearch(snes,&linesearch);CHKERRQ(ierr); 2555 ierr = SNESGetSNESLineSearch(snes->pc,&pclinesearch);CHKERRQ(ierr); 2556 ierr = SNESLineSearchGetPreCheck(linesearch,&precheck,&lsprectx);CHKERRQ(ierr); 2557 ierr = SNESLineSearchGetPostCheck(linesearch,&postcheck,&lspostctx);CHKERRQ(ierr); 2558 ierr = SNESLineSearchSetPreCheck(pclinesearch,precheck,lsprectx);CHKERRQ(ierr); 2559 ierr = SNESLineSearchSetPostCheck(pclinesearch,postcheck,lspostctx);CHKERRQ(ierr); 2560 ierr = PetscObjectCopyFortranFunctionPointers((PetscObject)linesearch, (PetscObject)pclinesearch);CHKERRQ(ierr); 2561 } 2562 2563 if (snes->ops->setup) { 2564 ierr = (*snes->ops->setup)(snes);CHKERRQ(ierr); 2565 } 2566 2567 snes->setupcalled = PETSC_TRUE; 2568 PetscFunctionReturn(0); 2569 } 2570 2571 #undef __FUNCT__ 2572 #define __FUNCT__ "SNESReset" 2573 /*@ 2574 SNESReset - Resets a SNES context to the snessetupcalled = 0 state and removes any allocated Vecs and Mats 2575 2576 Collective on SNES 2577 2578 Input Parameter: 2579 . snes - iterative context obtained from SNESCreate() 2580 2581 Level: intermediate 2582 2583 Notes: Also calls the application context destroy routine set with SNESSetComputeApplicationContext() 2584 2585 .keywords: SNES, destroy 2586 2587 .seealso: SNESCreate(), SNESSetUp(), SNESSolve() 2588 @*/ 2589 PetscErrorCode SNESReset(SNES snes) 2590 { 2591 PetscErrorCode ierr; 2592 2593 PetscFunctionBegin; 2594 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2595 if (snes->ops->userdestroy && snes->user) { 2596 ierr = (*snes->ops->userdestroy)((void**)&snes->user);CHKERRQ(ierr); 2597 snes->user = NULL; 2598 } 2599 if (snes->pc) { 2600 ierr = SNESReset(snes->pc);CHKERRQ(ierr); 2601 } 2602 2603 if (snes->ops->reset) { 2604 ierr = (*snes->ops->reset)(snes);CHKERRQ(ierr); 2605 } 2606 if (snes->ksp) { 2607 ierr = KSPReset(snes->ksp);CHKERRQ(ierr); 2608 } 2609 2610 if (snes->linesearch) { 2611 ierr = SNESLineSearchReset(snes->linesearch);CHKERRQ(ierr); 2612 } 2613 2614 ierr = VecDestroy(&snes->vec_rhs);CHKERRQ(ierr); 2615 ierr = VecDestroy(&snes->vec_sol);CHKERRQ(ierr); 2616 ierr = VecDestroy(&snes->vec_sol_update);CHKERRQ(ierr); 2617 ierr = VecDestroy(&snes->vec_func);CHKERRQ(ierr); 2618 ierr = MatDestroy(&snes->jacobian);CHKERRQ(ierr); 2619 ierr = MatDestroy(&snes->jacobian_pre);CHKERRQ(ierr); 2620 ierr = VecDestroyVecs(snes->nwork,&snes->work);CHKERRQ(ierr); 2621 ierr = VecDestroyVecs(snes->nvwork,&snes->vwork);CHKERRQ(ierr); 2622 2623 snes->nwork = snes->nvwork = 0; 2624 snes->setupcalled = PETSC_FALSE; 2625 PetscFunctionReturn(0); 2626 } 2627 2628 #undef __FUNCT__ 2629 #define __FUNCT__ "SNESDestroy" 2630 /*@ 2631 SNESDestroy - Destroys the nonlinear solver context that was created 2632 with SNESCreate(). 2633 2634 Collective on SNES 2635 2636 Input Parameter: 2637 . snes - the SNES context 2638 2639 Level: beginner 2640 2641 .keywords: SNES, nonlinear, destroy 2642 2643 .seealso: SNESCreate(), SNESSolve() 2644 @*/ 2645 PetscErrorCode SNESDestroy(SNES *snes) 2646 { 2647 PetscErrorCode ierr; 2648 2649 PetscFunctionBegin; 2650 if (!*snes) PetscFunctionReturn(0); 2651 PetscValidHeaderSpecific((*snes),SNES_CLASSID,1); 2652 if (--((PetscObject)(*snes))->refct > 0) {*snes = 0; PetscFunctionReturn(0);} 2653 2654 ierr = SNESReset((*snes));CHKERRQ(ierr); 2655 ierr = SNESDestroy(&(*snes)->pc);CHKERRQ(ierr); 2656 2657 /* if memory was published with AMS then destroy it */ 2658 ierr = PetscObjectDepublish((*snes));CHKERRQ(ierr); 2659 if ((*snes)->ops->destroy) {ierr = (*((*snes))->ops->destroy)((*snes));CHKERRQ(ierr);} 2660 2661 ierr = DMDestroy(&(*snes)->dm);CHKERRQ(ierr); 2662 ierr = KSPDestroy(&(*snes)->ksp);CHKERRQ(ierr); 2663 ierr = SNESLineSearchDestroy(&(*snes)->linesearch);CHKERRQ(ierr); 2664 2665 ierr = PetscFree((*snes)->kspconvctx);CHKERRQ(ierr); 2666 if ((*snes)->ops->convergeddestroy) { 2667 ierr = (*(*snes)->ops->convergeddestroy)((*snes)->cnvP);CHKERRQ(ierr); 2668 } 2669 if ((*snes)->conv_malloc) { 2670 ierr = PetscFree((*snes)->conv_hist);CHKERRQ(ierr); 2671 ierr = PetscFree((*snes)->conv_hist_its);CHKERRQ(ierr); 2672 } 2673 ierr = SNESMonitorCancel((*snes));CHKERRQ(ierr); 2674 ierr = PetscHeaderDestroy(snes);CHKERRQ(ierr); 2675 PetscFunctionReturn(0); 2676 } 2677 2678 /* ----------- Routines to set solver parameters ---------- */ 2679 2680 #undef __FUNCT__ 2681 #define __FUNCT__ "SNESSetLagPreconditioner" 2682 /*@ 2683 SNESSetLagPreconditioner - Determines when the preconditioner is rebuilt in the nonlinear solve. 2684 2685 Logically Collective on SNES 2686 2687 Input Parameters: 2688 + snes - the SNES context 2689 - lag - -1 indicates NEVER rebuild, 1 means rebuild every time the Jacobian is computed within a single nonlinear solve, 2 means every second time 2690 the Jacobian is built etc. -2 indicates rebuild preconditioner at next chance but then never rebuild after that 2691 2692 Options Database Keys: 2693 . -snes_lag_preconditioner <lag> 2694 2695 Notes: 2696 The default is 1 2697 The preconditioner is ALWAYS built in the first iteration of a nonlinear solve unless lag is -1 2698 If -1 is used before the very first nonlinear solve the preconditioner is still built because there is no previous preconditioner to use 2699 2700 Level: intermediate 2701 2702 .keywords: SNES, nonlinear, set, convergence, tolerances 2703 2704 .seealso: SNESSetTrustRegionTolerance(), SNESGetLagPreconditioner(), SNESSetLagJacobian(), SNESGetLagJacobian() 2705 2706 @*/ 2707 PetscErrorCode SNESSetLagPreconditioner(SNES snes,PetscInt lag) 2708 { 2709 PetscFunctionBegin; 2710 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2711 if (lag < -2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Lag must be -2, -1, 1 or greater"); 2712 if (!lag) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Lag cannot be 0"); 2713 PetscValidLogicalCollectiveInt(snes,lag,2); 2714 snes->lagpreconditioner = lag; 2715 PetscFunctionReturn(0); 2716 } 2717 2718 #undef __FUNCT__ 2719 #define __FUNCT__ "SNESSetGridSequence" 2720 /*@ 2721 SNESSetGridSequence - sets the number of steps of grid sequencing that SNES does 2722 2723 Logically Collective on SNES 2724 2725 Input Parameters: 2726 + snes - the SNES context 2727 - steps - the number of refinements to do, defaults to 0 2728 2729 Options Database Keys: 2730 . -snes_grid_sequence <steps> 2731 2732 Level: intermediate 2733 2734 Notes: 2735 Use SNESGetSolution() to extract the fine grid solution after grid sequencing. 2736 2737 .keywords: SNES, nonlinear, set, convergence, tolerances 2738 2739 .seealso: SNESSetTrustRegionTolerance(), SNESGetLagPreconditioner(), SNESSetLagJacobian(), SNESGetLagJacobian() 2740 2741 @*/ 2742 PetscErrorCode SNESSetGridSequence(SNES snes,PetscInt steps) 2743 { 2744 PetscFunctionBegin; 2745 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2746 PetscValidLogicalCollectiveInt(snes,steps,2); 2747 snes->gridsequence = steps; 2748 PetscFunctionReturn(0); 2749 } 2750 2751 #undef __FUNCT__ 2752 #define __FUNCT__ "SNESGetLagPreconditioner" 2753 /*@ 2754 SNESGetLagPreconditioner - Indicates how often the preconditioner is rebuilt 2755 2756 Not Collective 2757 2758 Input Parameter: 2759 . snes - the SNES context 2760 2761 Output Parameter: 2762 . lag - -1 indicates NEVER rebuild, 1 means rebuild every time the Jacobian is computed within a single nonlinear solve, 2 means every second time 2763 the Jacobian is built etc. -2 indicates rebuild preconditioner at next chance but then never rebuild after that 2764 2765 Options Database Keys: 2766 . -snes_lag_preconditioner <lag> 2767 2768 Notes: 2769 The default is 1 2770 The preconditioner is ALWAYS built in the first iteration of a nonlinear solve unless lag is -1 2771 2772 Level: intermediate 2773 2774 .keywords: SNES, nonlinear, set, convergence, tolerances 2775 2776 .seealso: SNESSetTrustRegionTolerance(), SNESSetLagPreconditioner() 2777 2778 @*/ 2779 PetscErrorCode SNESGetLagPreconditioner(SNES snes,PetscInt *lag) 2780 { 2781 PetscFunctionBegin; 2782 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2783 *lag = snes->lagpreconditioner; 2784 PetscFunctionReturn(0); 2785 } 2786 2787 #undef __FUNCT__ 2788 #define __FUNCT__ "SNESSetLagJacobian" 2789 /*@ 2790 SNESSetLagJacobian - Determines when the Jacobian is rebuilt in the nonlinear solve. See SNESSetLagPreconditioner() for determining how 2791 often the preconditioner is rebuilt. 2792 2793 Logically Collective on SNES 2794 2795 Input Parameters: 2796 + snes - the SNES context 2797 - lag - -1 indicates NEVER rebuild, 1 means rebuild every time the Jacobian is computed within a single nonlinear solve, 2 means every second time 2798 the Jacobian is built etc. -2 means rebuild at next chance but then never again 2799 2800 Options Database Keys: 2801 . -snes_lag_jacobian <lag> 2802 2803 Notes: 2804 The default is 1 2805 The Jacobian is ALWAYS built in the first iteration of a nonlinear solve unless lag is -1 2806 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 2807 at the next Newton step but never again (unless it is reset to another value) 2808 2809 Level: intermediate 2810 2811 .keywords: SNES, nonlinear, set, convergence, tolerances 2812 2813 .seealso: SNESSetTrustRegionTolerance(), SNESGetLagPreconditioner(), SNESSetLagPreconditioner(), SNESGetLagJacobian() 2814 2815 @*/ 2816 PetscErrorCode SNESSetLagJacobian(SNES snes,PetscInt lag) 2817 { 2818 PetscFunctionBegin; 2819 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2820 if (lag < -2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Lag must be -2, -1, 1 or greater"); 2821 if (!lag) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Lag cannot be 0"); 2822 PetscValidLogicalCollectiveInt(snes,lag,2); 2823 snes->lagjacobian = lag; 2824 PetscFunctionReturn(0); 2825 } 2826 2827 #undef __FUNCT__ 2828 #define __FUNCT__ "SNESGetLagJacobian" 2829 /*@ 2830 SNESGetLagJacobian - Indicates how often the Jacobian is rebuilt. See SNESGetLagPreconditioner() to determine when the preconditioner is rebuilt 2831 2832 Not Collective 2833 2834 Input Parameter: 2835 . snes - the SNES context 2836 2837 Output Parameter: 2838 . lag - -1 indicates NEVER rebuild, 1 means rebuild every time the Jacobian is computed within a single nonlinear solve, 2 means every second time 2839 the Jacobian is built etc. 2840 2841 Options Database Keys: 2842 . -snes_lag_jacobian <lag> 2843 2844 Notes: 2845 The default is 1 2846 The jacobian is ALWAYS built in the first iteration of a nonlinear solve unless lag is -1 2847 2848 Level: intermediate 2849 2850 .keywords: SNES, nonlinear, set, convergence, tolerances 2851 2852 .seealso: SNESSetTrustRegionTolerance(), SNESSetLagJacobian(), SNESSetLagPreconditioner(), SNESGetLagPreconditioner() 2853 2854 @*/ 2855 PetscErrorCode SNESGetLagJacobian(SNES snes,PetscInt *lag) 2856 { 2857 PetscFunctionBegin; 2858 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2859 *lag = snes->lagjacobian; 2860 PetscFunctionReturn(0); 2861 } 2862 2863 #undef __FUNCT__ 2864 #define __FUNCT__ "SNESSetTolerances" 2865 /*@ 2866 SNESSetTolerances - Sets various parameters used in convergence tests. 2867 2868 Logically Collective on SNES 2869 2870 Input Parameters: 2871 + snes - the SNES context 2872 . abstol - absolute convergence tolerance 2873 . rtol - relative convergence tolerance 2874 . stol - convergence tolerance in terms of the norm of the change in the solution between steps, || delta x || < stol*|| x || 2875 . maxit - maximum number of iterations 2876 - maxf - maximum number of function evaluations 2877 2878 Options Database Keys: 2879 + -snes_atol <abstol> - Sets abstol 2880 . -snes_rtol <rtol> - Sets rtol 2881 . -snes_stol <stol> - Sets stol 2882 . -snes_max_it <maxit> - Sets maxit 2883 - -snes_max_funcs <maxf> - Sets maxf 2884 2885 Notes: 2886 The default maximum number of iterations is 50. 2887 The default maximum number of function evaluations is 1000. 2888 2889 Level: intermediate 2890 2891 .keywords: SNES, nonlinear, set, convergence, tolerances 2892 2893 .seealso: SNESSetTrustRegionTolerance() 2894 @*/ 2895 PetscErrorCode SNESSetTolerances(SNES snes,PetscReal abstol,PetscReal rtol,PetscReal stol,PetscInt maxit,PetscInt maxf) 2896 { 2897 PetscFunctionBegin; 2898 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2899 PetscValidLogicalCollectiveReal(snes,abstol,2); 2900 PetscValidLogicalCollectiveReal(snes,rtol,3); 2901 PetscValidLogicalCollectiveReal(snes,stol,4); 2902 PetscValidLogicalCollectiveInt(snes,maxit,5); 2903 PetscValidLogicalCollectiveInt(snes,maxf,6); 2904 2905 if (abstol != PETSC_DEFAULT) { 2906 if (abstol < 0.0) SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_OUTOFRANGE,"Absolute tolerance %G must be non-negative",abstol); 2907 snes->abstol = abstol; 2908 } 2909 if (rtol != PETSC_DEFAULT) { 2910 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); 2911 snes->rtol = rtol; 2912 } 2913 if (stol != PETSC_DEFAULT) { 2914 if (stol < 0.0) SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_OUTOFRANGE,"Step tolerance %G must be non-negative",stol); 2915 snes->stol = stol; 2916 } 2917 if (maxit != PETSC_DEFAULT) { 2918 if (maxit < 0) SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_OUTOFRANGE,"Maximum number of iterations %D must be non-negative",maxit); 2919 snes->max_its = maxit; 2920 } 2921 if (maxf != PETSC_DEFAULT) { 2922 if (maxf < 0) SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_OUTOFRANGE,"Maximum number of function evaluations %D must be non-negative",maxf); 2923 snes->max_funcs = maxf; 2924 } 2925 snes->tolerancesset = PETSC_TRUE; 2926 PetscFunctionReturn(0); 2927 } 2928 2929 #undef __FUNCT__ 2930 #define __FUNCT__ "SNESGetTolerances" 2931 /*@ 2932 SNESGetTolerances - Gets various parameters used in convergence tests. 2933 2934 Not Collective 2935 2936 Input Parameters: 2937 + snes - the SNES context 2938 . atol - absolute convergence tolerance 2939 . rtol - relative convergence tolerance 2940 . stol - convergence tolerance in terms of the norm 2941 of the change in the solution between steps 2942 . maxit - maximum number of iterations 2943 - maxf - maximum number of function evaluations 2944 2945 Notes: 2946 The user can specify NULL for any parameter that is not needed. 2947 2948 Level: intermediate 2949 2950 .keywords: SNES, nonlinear, get, convergence, tolerances 2951 2952 .seealso: SNESSetTolerances() 2953 @*/ 2954 PetscErrorCode SNESGetTolerances(SNES snes,PetscReal *atol,PetscReal *rtol,PetscReal *stol,PetscInt *maxit,PetscInt *maxf) 2955 { 2956 PetscFunctionBegin; 2957 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2958 if (atol) *atol = snes->abstol; 2959 if (rtol) *rtol = snes->rtol; 2960 if (stol) *stol = snes->stol; 2961 if (maxit) *maxit = snes->max_its; 2962 if (maxf) *maxf = snes->max_funcs; 2963 PetscFunctionReturn(0); 2964 } 2965 2966 #undef __FUNCT__ 2967 #define __FUNCT__ "SNESSetTrustRegionTolerance" 2968 /*@ 2969 SNESSetTrustRegionTolerance - Sets the trust region parameter tolerance. 2970 2971 Logically Collective on SNES 2972 2973 Input Parameters: 2974 + snes - the SNES context 2975 - tol - tolerance 2976 2977 Options Database Key: 2978 . -snes_trtol <tol> - Sets tol 2979 2980 Level: intermediate 2981 2982 .keywords: SNES, nonlinear, set, trust region, tolerance 2983 2984 .seealso: SNESSetTolerances() 2985 @*/ 2986 PetscErrorCode SNESSetTrustRegionTolerance(SNES snes,PetscReal tol) 2987 { 2988 PetscFunctionBegin; 2989 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2990 PetscValidLogicalCollectiveReal(snes,tol,2); 2991 snes->deltatol = tol; 2992 PetscFunctionReturn(0); 2993 } 2994 2995 /* 2996 Duplicate the lg monitors for SNES from KSP; for some reason with 2997 dynamic libraries things don't work under Sun4 if we just use 2998 macros instead of functions 2999 */ 3000 #undef __FUNCT__ 3001 #define __FUNCT__ "SNESMonitorLGResidualNorm" 3002 PetscErrorCode SNESMonitorLGResidualNorm(SNES snes,PetscInt it,PetscReal norm,void *ctx) 3003 { 3004 PetscErrorCode ierr; 3005 3006 PetscFunctionBegin; 3007 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3008 ierr = KSPMonitorLGResidualNorm((KSP)snes,it,norm,ctx);CHKERRQ(ierr); 3009 PetscFunctionReturn(0); 3010 } 3011 3012 #undef __FUNCT__ 3013 #define __FUNCT__ "SNESMonitorLGCreate" 3014 PetscErrorCode SNESMonitorLGCreate(const char host[],const char label[],int x,int y,int m,int n,PetscDrawLG *draw) 3015 { 3016 PetscErrorCode ierr; 3017 3018 PetscFunctionBegin; 3019 ierr = KSPMonitorLGResidualNormCreate(host,label,x,y,m,n,draw);CHKERRQ(ierr); 3020 PetscFunctionReturn(0); 3021 } 3022 3023 #undef __FUNCT__ 3024 #define __FUNCT__ "SNESMonitorLGDestroy" 3025 PetscErrorCode SNESMonitorLGDestroy(PetscDrawLG *draw) 3026 { 3027 PetscErrorCode ierr; 3028 3029 PetscFunctionBegin; 3030 ierr = KSPMonitorLGResidualNormDestroy(draw);CHKERRQ(ierr); 3031 PetscFunctionReturn(0); 3032 } 3033 3034 extern PetscErrorCode SNESMonitorRange_Private(SNES,PetscInt,PetscReal*); 3035 #undef __FUNCT__ 3036 #define __FUNCT__ "SNESMonitorLGRange" 3037 PetscErrorCode SNESMonitorLGRange(SNES snes,PetscInt n,PetscReal rnorm,void *monctx) 3038 { 3039 PetscDrawLG lg; 3040 PetscErrorCode ierr; 3041 PetscReal x,y,per; 3042 PetscViewer v = (PetscViewer)monctx; 3043 static PetscReal prev; /* should be in the context */ 3044 PetscDraw draw; 3045 3046 PetscFunctionBegin; 3047 ierr = PetscViewerDrawGetDrawLG(v,0,&lg);CHKERRQ(ierr); 3048 if (!n) {ierr = PetscDrawLGReset(lg);CHKERRQ(ierr);} 3049 ierr = PetscDrawLGGetDraw(lg,&draw);CHKERRQ(ierr); 3050 ierr = PetscDrawSetTitle(draw,"Residual norm");CHKERRQ(ierr); 3051 x = (PetscReal)n; 3052 if (rnorm > 0.0) y = log10(rnorm); 3053 else y = -15.0; 3054 ierr = PetscDrawLGAddPoint(lg,&x,&y);CHKERRQ(ierr); 3055 if (n < 20 || !(n % 5)) { 3056 ierr = PetscDrawLGDraw(lg);CHKERRQ(ierr); 3057 } 3058 3059 ierr = PetscViewerDrawGetDrawLG(v,1,&lg);CHKERRQ(ierr); 3060 if (!n) {ierr = PetscDrawLGReset(lg);CHKERRQ(ierr);} 3061 ierr = PetscDrawLGGetDraw(lg,&draw);CHKERRQ(ierr); 3062 ierr = PetscDrawSetTitle(draw,"% elemts > .2*max elemt");CHKERRQ(ierr); 3063 ierr = SNESMonitorRange_Private(snes,n,&per);CHKERRQ(ierr); 3064 x = (PetscReal)n; 3065 y = 100.0*per; 3066 ierr = PetscDrawLGAddPoint(lg,&x,&y);CHKERRQ(ierr); 3067 if (n < 20 || !(n % 5)) { 3068 ierr = PetscDrawLGDraw(lg);CHKERRQ(ierr); 3069 } 3070 3071 ierr = PetscViewerDrawGetDrawLG(v,2,&lg);CHKERRQ(ierr); 3072 if (!n) {prev = rnorm;ierr = PetscDrawLGReset(lg);CHKERRQ(ierr);} 3073 ierr = PetscDrawLGGetDraw(lg,&draw);CHKERRQ(ierr); 3074 ierr = PetscDrawSetTitle(draw,"(norm -oldnorm)/oldnorm");CHKERRQ(ierr); 3075 x = (PetscReal)n; 3076 y = (prev - rnorm)/prev; 3077 ierr = PetscDrawLGAddPoint(lg,&x,&y);CHKERRQ(ierr); 3078 if (n < 20 || !(n % 5)) { 3079 ierr = PetscDrawLGDraw(lg);CHKERRQ(ierr); 3080 } 3081 3082 ierr = PetscViewerDrawGetDrawLG(v,3,&lg);CHKERRQ(ierr); 3083 if (!n) {ierr = PetscDrawLGReset(lg);CHKERRQ(ierr);} 3084 ierr = PetscDrawLGGetDraw(lg,&draw);CHKERRQ(ierr); 3085 ierr = PetscDrawSetTitle(draw,"(norm -oldnorm)/oldnorm*(% > .2 max)");CHKERRQ(ierr); 3086 x = (PetscReal)n; 3087 y = (prev - rnorm)/(prev*per); 3088 if (n > 2) { /*skip initial crazy value */ 3089 ierr = PetscDrawLGAddPoint(lg,&x,&y);CHKERRQ(ierr); 3090 } 3091 if (n < 20 || !(n % 5)) { 3092 ierr = PetscDrawLGDraw(lg);CHKERRQ(ierr); 3093 } 3094 prev = rnorm; 3095 PetscFunctionReturn(0); 3096 } 3097 3098 #undef __FUNCT__ 3099 #define __FUNCT__ "SNESMonitor" 3100 /*@ 3101 SNESMonitor - runs the user provided monitor routines, if they exist 3102 3103 Collective on SNES 3104 3105 Input Parameters: 3106 + snes - nonlinear solver context obtained from SNESCreate() 3107 . iter - iteration number 3108 - rnorm - relative norm of the residual 3109 3110 Notes: 3111 This routine is called by the SNES implementations. 3112 It does not typically need to be called by the user. 3113 3114 Level: developer 3115 3116 .seealso: SNESMonitorSet() 3117 @*/ 3118 PetscErrorCode SNESMonitor(SNES snes,PetscInt iter,PetscReal rnorm) 3119 { 3120 PetscErrorCode ierr; 3121 PetscInt i,n = snes->numbermonitors; 3122 3123 PetscFunctionBegin; 3124 for (i=0; i<n; i++) { 3125 ierr = (*snes->monitor[i])(snes,iter,rnorm,snes->monitorcontext[i]);CHKERRQ(ierr); 3126 } 3127 PetscFunctionReturn(0); 3128 } 3129 3130 /* ------------ Routines to set performance monitoring options ----------- */ 3131 3132 /*MC 3133 SNESMonitorFunction - functional form passed to SNESMonitorSet() to monitor convergence of nonlinear solver 3134 3135 Synopsis: 3136 #include "petscsnes.h" 3137 $ PetscErrorCode SNESMonitorFunction(SNES snes,PetscInt its, PetscReal norm,void *mctx) 3138 3139 + snes - the SNES context 3140 . its - iteration number 3141 . norm - 2-norm function value (may be estimated) 3142 - mctx - [optional] monitoring context 3143 3144 Level: advanced 3145 3146 .seealso: SNESMonitorSet(), SNESMonitorGet() 3147 M*/ 3148 3149 #undef __FUNCT__ 3150 #define __FUNCT__ "SNESMonitorSet" 3151 /*@C 3152 SNESMonitorSet - Sets an ADDITIONAL function that is to be used at every 3153 iteration of the nonlinear solver to display the iteration's 3154 progress. 3155 3156 Logically Collective on SNES 3157 3158 Input Parameters: 3159 + snes - the SNES context 3160 . SNESMonitorFunction - monitoring routine 3161 . mctx - [optional] user-defined context for private data for the 3162 monitor routine (use NULL if no context is desired) 3163 - monitordestroy - [optional] routine that frees monitor context 3164 (may be NULL) 3165 3166 Options Database Keys: 3167 + -snes_monitor - sets SNESMonitorDefault() 3168 . -snes_monitor_lg_residualnorm - sets line graph monitor, 3169 uses SNESMonitorLGCreate() 3170 - -snes_monitor_cancel - cancels all monitors that have 3171 been hardwired into a code by 3172 calls to SNESMonitorSet(), but 3173 does not cancel those set via 3174 the options database. 3175 3176 Notes: 3177 Several different monitoring routines may be set by calling 3178 SNESMonitorSet() multiple times; all will be called in the 3179 order in which they were set. 3180 3181 Fortran notes: Only a single monitor function can be set for each SNES object 3182 3183 Level: intermediate 3184 3185 .keywords: SNES, nonlinear, set, monitor 3186 3187 .seealso: SNESMonitorDefault(), SNESMonitorCancel(), SNESMonitorFunction 3188 @*/ 3189 PetscErrorCode SNESMonitorSet(SNES snes,PetscErrorCode (*SNESMonitorFunction)(SNES,PetscInt,PetscReal,void*),void *mctx,PetscErrorCode (*monitordestroy)(void**)) 3190 { 3191 PetscInt i; 3192 PetscErrorCode ierr; 3193 3194 PetscFunctionBegin; 3195 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3196 if (snes->numbermonitors >= MAXSNESMONITORS) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Too many monitors set"); 3197 for (i=0; i<snes->numbermonitors;i++) { 3198 if (SNESMonitorFunction == snes->monitor[i] && monitordestroy == snes->monitordestroy[i] && mctx == snes->monitorcontext[i]) { 3199 if (monitordestroy) { 3200 ierr = (*monitordestroy)(&mctx);CHKERRQ(ierr); 3201 } 3202 PetscFunctionReturn(0); 3203 } 3204 } 3205 snes->monitor[snes->numbermonitors] = SNESMonitorFunction; 3206 snes->monitordestroy[snes->numbermonitors] = monitordestroy; 3207 snes->monitorcontext[snes->numbermonitors++] = (void*)mctx; 3208 PetscFunctionReturn(0); 3209 } 3210 3211 #undef __FUNCT__ 3212 #define __FUNCT__ "SNESMonitorCancel" 3213 /*@C 3214 SNESMonitorCancel - Clears all the monitor functions for a SNES object. 3215 3216 Logically Collective on SNES 3217 3218 Input Parameters: 3219 . snes - the SNES context 3220 3221 Options Database Key: 3222 . -snes_monitor_cancel - cancels all monitors that have been hardwired 3223 into a code by calls to SNESMonitorSet(), but does not cancel those 3224 set via the options database 3225 3226 Notes: 3227 There is no way to clear one specific monitor from a SNES object. 3228 3229 Level: intermediate 3230 3231 .keywords: SNES, nonlinear, set, monitor 3232 3233 .seealso: SNESMonitorDefault(), SNESMonitorSet() 3234 @*/ 3235 PetscErrorCode SNESMonitorCancel(SNES snes) 3236 { 3237 PetscErrorCode ierr; 3238 PetscInt i; 3239 3240 PetscFunctionBegin; 3241 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3242 for (i=0; i<snes->numbermonitors; i++) { 3243 if (snes->monitordestroy[i]) { 3244 ierr = (*snes->monitordestroy[i])(&snes->monitorcontext[i]);CHKERRQ(ierr); 3245 } 3246 } 3247 snes->numbermonitors = 0; 3248 PetscFunctionReturn(0); 3249 } 3250 3251 /*MC 3252 SNESConvergenceTestFunction - functional form used for testing of convergence of nonlinear solver 3253 3254 Synopsis: 3255 #include "petscsnes.h" 3256 $ PetscErrorCode SNESConvergenceTest(SNES snes,PetscInt it,PetscReal xnorm,PetscReal gnorm,PetscReal f,SNESConvergedReason *reason,void *cctx) 3257 3258 + snes - the SNES context 3259 . it - current iteration (0 is the first and is before any Newton step) 3260 . cctx - [optional] convergence context 3261 . reason - reason for convergence/divergence 3262 . xnorm - 2-norm of current iterate 3263 . gnorm - 2-norm of current step 3264 - f - 2-norm of function 3265 3266 Level: intermediate 3267 3268 .seealso: SNESSetConvergenceTest(), SNESGetConvergenceTest() 3269 M*/ 3270 3271 #undef __FUNCT__ 3272 #define __FUNCT__ "SNESSetConvergenceTest" 3273 /*@C 3274 SNESSetConvergenceTest - Sets the function that is to be used 3275 to test for convergence of the nonlinear iterative solution. 3276 3277 Logically Collective on SNES 3278 3279 Input Parameters: 3280 + snes - the SNES context 3281 . SNESConvergenceTestFunction - routine to test for convergence 3282 . cctx - [optional] context for private data for the convergence routine (may be NULL) 3283 - destroy - [optional] destructor for the context (may be NULL; NULL_FUNCTION in Fortran) 3284 3285 Level: advanced 3286 3287 .keywords: SNES, nonlinear, set, convergence, test 3288 3289 .seealso: SNESConvergedDefault(), SNESSkipConverged(), SNESConvergenceTestFunction 3290 @*/ 3291 PetscErrorCode SNESSetConvergenceTest(SNES snes,PetscErrorCode (*SNESConvergenceTestFunction)(SNES,PetscInt,PetscReal,PetscReal,PetscReal,SNESConvergedReason*,void*),void *cctx,PetscErrorCode (*destroy)(void*)) 3292 { 3293 PetscErrorCode ierr; 3294 3295 PetscFunctionBegin; 3296 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3297 if (!SNESConvergenceTestFunction) SNESConvergenceTestFunction = SNESSkipConverged; 3298 if (snes->ops->convergeddestroy) { 3299 ierr = (*snes->ops->convergeddestroy)(snes->cnvP);CHKERRQ(ierr); 3300 } 3301 snes->ops->converged = SNESConvergenceTestFunction; 3302 snes->ops->convergeddestroy = destroy; 3303 snes->cnvP = cctx; 3304 PetscFunctionReturn(0); 3305 } 3306 3307 #undef __FUNCT__ 3308 #define __FUNCT__ "SNESGetConvergedReason" 3309 /*@ 3310 SNESGetConvergedReason - Gets the reason the SNES iteration was stopped. 3311 3312 Not Collective 3313 3314 Input Parameter: 3315 . snes - the SNES context 3316 3317 Output Parameter: 3318 . reason - negative value indicates diverged, positive value converged, see SNESConvergedReason or the 3319 manual pages for the individual convergence tests for complete lists 3320 3321 Level: intermediate 3322 3323 Notes: Can only be called after the call the SNESSolve() is complete. 3324 3325 .keywords: SNES, nonlinear, set, convergence, test 3326 3327 .seealso: SNESSetConvergenceTest(), SNESConvergedReason 3328 @*/ 3329 PetscErrorCode SNESGetConvergedReason(SNES snes,SNESConvergedReason *reason) 3330 { 3331 PetscFunctionBegin; 3332 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3333 PetscValidPointer(reason,2); 3334 *reason = snes->reason; 3335 PetscFunctionReturn(0); 3336 } 3337 3338 #undef __FUNCT__ 3339 #define __FUNCT__ "SNESSetConvergenceHistory" 3340 /*@ 3341 SNESSetConvergenceHistory - Sets the array used to hold the convergence history. 3342 3343 Logically Collective on SNES 3344 3345 Input Parameters: 3346 + snes - iterative context obtained from SNESCreate() 3347 . a - array to hold history, this array will contain the function norms computed at each step 3348 . its - integer array holds the number of linear iterations for each solve. 3349 . na - size of a and its 3350 - reset - PETSC_TRUE indicates each new nonlinear solve resets the history counter to zero, 3351 else it continues storing new values for new nonlinear solves after the old ones 3352 3353 Notes: 3354 If 'a' and 'its' are NULL then space is allocated for the history. If 'na' PETSC_DECIDE or PETSC_DEFAULT then a 3355 default array of length 10000 is allocated. 3356 3357 This routine is useful, e.g., when running a code for purposes 3358 of accurate performance monitoring, when no I/O should be done 3359 during the section of code that is being timed. 3360 3361 Level: intermediate 3362 3363 .keywords: SNES, set, convergence, history 3364 3365 .seealso: SNESGetConvergenceHistory() 3366 3367 @*/ 3368 PetscErrorCode SNESSetConvergenceHistory(SNES snes,PetscReal a[],PetscInt its[],PetscInt na,PetscBool reset) 3369 { 3370 PetscErrorCode ierr; 3371 3372 PetscFunctionBegin; 3373 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3374 if (na) PetscValidScalarPointer(a,2); 3375 if (its) PetscValidIntPointer(its,3); 3376 if (na == PETSC_DECIDE || na == PETSC_DEFAULT || !a) { 3377 if (na == PETSC_DECIDE || na == PETSC_DEFAULT) na = 1000; 3378 ierr = PetscMalloc(na*sizeof(PetscReal),&a);CHKERRQ(ierr); 3379 ierr = PetscMalloc(na*sizeof(PetscInt),&its);CHKERRQ(ierr); 3380 3381 snes->conv_malloc = PETSC_TRUE; 3382 } 3383 snes->conv_hist = a; 3384 snes->conv_hist_its = its; 3385 snes->conv_hist_max = na; 3386 snes->conv_hist_len = 0; 3387 snes->conv_hist_reset = reset; 3388 PetscFunctionReturn(0); 3389 } 3390 3391 #if defined(PETSC_HAVE_MATLAB_ENGINE) 3392 #include <engine.h> /* MATLAB include file */ 3393 #include <mex.h> /* MATLAB include file */ 3394 3395 #undef __FUNCT__ 3396 #define __FUNCT__ "SNESGetConvergenceHistoryMatlab" 3397 PETSC_EXTERN mxArray *SNESGetConvergenceHistoryMatlab(SNES snes) 3398 { 3399 mxArray *mat; 3400 PetscInt i; 3401 PetscReal *ar; 3402 3403 PetscFunctionBegin; 3404 mat = mxCreateDoubleMatrix(snes->conv_hist_len,1,mxREAL); 3405 ar = (PetscReal*) mxGetData(mat); 3406 for (i=0; i<snes->conv_hist_len; i++) ar[i] = snes->conv_hist[i]; 3407 PetscFunctionReturn(mat); 3408 } 3409 #endif 3410 3411 #undef __FUNCT__ 3412 #define __FUNCT__ "SNESGetConvergenceHistory" 3413 /*@C 3414 SNESGetConvergenceHistory - Gets the array used to hold the convergence history. 3415 3416 Not Collective 3417 3418 Input Parameter: 3419 . snes - iterative context obtained from SNESCreate() 3420 3421 Output Parameters: 3422 . a - array to hold history 3423 . its - integer array holds the number of linear iterations (or 3424 negative if not converged) for each solve. 3425 - na - size of a and its 3426 3427 Notes: 3428 The calling sequence for this routine in Fortran is 3429 $ call SNESGetConvergenceHistory(SNES snes, integer na, integer ierr) 3430 3431 This routine is useful, e.g., when running a code for purposes 3432 of accurate performance monitoring, when no I/O should be done 3433 during the section of code that is being timed. 3434 3435 Level: intermediate 3436 3437 .keywords: SNES, get, convergence, history 3438 3439 .seealso: SNESSetConvergencHistory() 3440 3441 @*/ 3442 PetscErrorCode SNESGetConvergenceHistory(SNES snes,PetscReal *a[],PetscInt *its[],PetscInt *na) 3443 { 3444 PetscFunctionBegin; 3445 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3446 if (a) *a = snes->conv_hist; 3447 if (its) *its = snes->conv_hist_its; 3448 if (na) *na = snes->conv_hist_len; 3449 PetscFunctionReturn(0); 3450 } 3451 3452 #undef __FUNCT__ 3453 #define __FUNCT__ "SNESSetUpdate" 3454 /*@C 3455 SNESSetUpdate - Sets the general-purpose update function called 3456 at the beginning of every iteration of the nonlinear solve. Specifically 3457 it is called just before the Jacobian is "evaluated". 3458 3459 Logically Collective on SNES 3460 3461 Input Parameters: 3462 . snes - The nonlinear solver context 3463 . func - The function 3464 3465 Calling sequence of func: 3466 . func (SNES snes, PetscInt step); 3467 3468 . step - The current step of the iteration 3469 3470 Level: advanced 3471 3472 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() 3473 This is not used by most users. 3474 3475 .keywords: SNES, update 3476 3477 .seealso SNESSetJacobian(), SNESSolve() 3478 @*/ 3479 PetscErrorCode SNESSetUpdate(SNES snes, PetscErrorCode (*func)(SNES, PetscInt)) 3480 { 3481 PetscFunctionBegin; 3482 PetscValidHeaderSpecific(snes, SNES_CLASSID,1); 3483 snes->ops->update = func; 3484 PetscFunctionReturn(0); 3485 } 3486 3487 #undef __FUNCT__ 3488 #define __FUNCT__ "SNESScaleStep_Private" 3489 /* 3490 SNESScaleStep_Private - Scales a step so that its length is less than the 3491 positive parameter delta. 3492 3493 Input Parameters: 3494 + snes - the SNES context 3495 . y - approximate solution of linear system 3496 . fnorm - 2-norm of current function 3497 - delta - trust region size 3498 3499 Output Parameters: 3500 + gpnorm - predicted function norm at the new point, assuming local 3501 linearization. The value is zero if the step lies within the trust 3502 region, and exceeds zero otherwise. 3503 - ynorm - 2-norm of the step 3504 3505 Note: 3506 For non-trust region methods such as SNESNEWTONLS, the parameter delta 3507 is set to be the maximum allowable step size. 3508 3509 .keywords: SNES, nonlinear, scale, step 3510 */ 3511 PetscErrorCode SNESScaleStep_Private(SNES snes,Vec y,PetscReal *fnorm,PetscReal *delta,PetscReal *gpnorm,PetscReal *ynorm) 3512 { 3513 PetscReal nrm; 3514 PetscScalar cnorm; 3515 PetscErrorCode ierr; 3516 3517 PetscFunctionBegin; 3518 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3519 PetscValidHeaderSpecific(y,VEC_CLASSID,2); 3520 PetscCheckSameComm(snes,1,y,2); 3521 3522 ierr = VecNorm(y,NORM_2,&nrm);CHKERRQ(ierr); 3523 if (nrm > *delta) { 3524 nrm = *delta/nrm; 3525 *gpnorm = (1.0 - nrm)*(*fnorm); 3526 cnorm = nrm; 3527 ierr = VecScale(y,cnorm);CHKERRQ(ierr); 3528 *ynorm = *delta; 3529 } else { 3530 *gpnorm = 0.0; 3531 *ynorm = nrm; 3532 } 3533 PetscFunctionReturn(0); 3534 } 3535 3536 #undef __FUNCT__ 3537 #define __FUNCT__ "SNESSolve" 3538 /*@C 3539 SNESSolve - Solves a nonlinear system F(x) = b. 3540 Call SNESSolve() after calling SNESCreate() and optional routines of the form SNESSetXXX(). 3541 3542 Collective on SNES 3543 3544 Input Parameters: 3545 + snes - the SNES context 3546 . b - the constant part of the equation F(x) = b, or NULL to use zero. 3547 - x - the solution vector. 3548 3549 Notes: 3550 The user should initialize the vector,x, with the initial guess 3551 for the nonlinear solve prior to calling SNESSolve. In particular, 3552 to employ an initial guess of zero, the user should explicitly set 3553 this vector to zero by calling VecSet(). 3554 3555 Level: beginner 3556 3557 .keywords: SNES, nonlinear, solve 3558 3559 .seealso: SNESCreate(), SNESDestroy(), SNESSetFunction(), SNESSetJacobian(), SNESSetGridSequence(), SNESGetSolution() 3560 @*/ 3561 PetscErrorCode SNESSolve(SNES snes,Vec b,Vec x) 3562 { 3563 PetscErrorCode ierr; 3564 PetscBool flg; 3565 PetscViewer viewer; 3566 PetscInt grid; 3567 Vec xcreated = NULL; 3568 DM dm; 3569 PetscViewerFormat format; 3570 3571 PetscFunctionBegin; 3572 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3573 if (x) PetscValidHeaderSpecific(x,VEC_CLASSID,3); 3574 if (x) PetscCheckSameComm(snes,1,x,3); 3575 if (b) PetscValidHeaderSpecific(b,VEC_CLASSID,2); 3576 if (b) PetscCheckSameComm(snes,1,b,2); 3577 3578 if (!x) { 3579 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 3580 ierr = DMCreateGlobalVector(dm,&xcreated);CHKERRQ(ierr); 3581 x = xcreated; 3582 } 3583 3584 for (grid=0; grid<snes->gridsequence; grid++) {ierr = PetscViewerASCIIPushTab(PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)snes)));CHKERRQ(ierr);} 3585 for (grid=0; grid<snes->gridsequence+1; grid++) { 3586 3587 /* set solution vector */ 3588 if (!grid) {ierr = PetscObjectReference((PetscObject)x);CHKERRQ(ierr);} 3589 ierr = VecDestroy(&snes->vec_sol);CHKERRQ(ierr); 3590 snes->vec_sol = x; 3591 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 3592 3593 /* set affine vector if provided */ 3594 if (b) { ierr = PetscObjectReference((PetscObject)b);CHKERRQ(ierr); } 3595 ierr = VecDestroy(&snes->vec_rhs);CHKERRQ(ierr); 3596 snes->vec_rhs = b; 3597 3598 ierr = SNESSetUp(snes);CHKERRQ(ierr); 3599 3600 if (!grid) { 3601 if (snes->ops->computeinitialguess) { 3602 ierr = (*snes->ops->computeinitialguess)(snes,snes->vec_sol,snes->initialguessP);CHKERRQ(ierr); 3603 } 3604 } 3605 3606 if (snes->conv_hist_reset) snes->conv_hist_len = 0; 3607 snes->nfuncs = 0; snes->linear_its = 0; snes->numFailures = 0; 3608 3609 ierr = PetscLogEventBegin(SNES_Solve,snes,0,0,0);CHKERRQ(ierr); 3610 ierr = (*snes->ops->solve)(snes);CHKERRQ(ierr); 3611 ierr = PetscLogEventEnd(SNES_Solve,snes,0,0,0);CHKERRQ(ierr); 3612 if (snes->domainerror) { 3613 snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; 3614 snes->domainerror = PETSC_FALSE; 3615 } 3616 if (!snes->reason) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Internal error, solver returned without setting converged reason"); 3617 3618 flg = PETSC_FALSE; 3619 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_test_local_min",&flg,NULL);CHKERRQ(ierr); 3620 if (flg && !PetscPreLoadingOn) { ierr = SNESTestLocalMin(snes);CHKERRQ(ierr); } 3621 if (snes->printreason) { 3622 ierr = PetscViewerASCIIAddTab(PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)snes)),((PetscObject)snes)->tablevel);CHKERRQ(ierr); 3623 if (snes->reason > 0) { 3624 ierr = PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)snes)),"Nonlinear solve converged due to %s iterations %D\n",SNESConvergedReasons[snes->reason],snes->iter);CHKERRQ(ierr); 3625 } else { 3626 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); 3627 } 3628 ierr = PetscViewerASCIISubtractTab(PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)snes)),((PetscObject)snes)->tablevel);CHKERRQ(ierr); 3629 } 3630 3631 if (snes->errorifnotconverged && snes->reason < 0) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_NOT_CONVERGED,"SNESSolve has not converged"); 3632 if (grid < snes->gridsequence) { 3633 DM fine; 3634 Vec xnew; 3635 Mat interp; 3636 3637 ierr = DMRefine(snes->dm,PetscObjectComm((PetscObject)snes),&fine);CHKERRQ(ierr); 3638 if (!fine) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_INCOMP,"DMRefine() did not perform any refinement, cannot continue grid sequencing"); 3639 ierr = DMCreateInterpolation(snes->dm,fine,&interp,NULL);CHKERRQ(ierr); 3640 ierr = DMCreateGlobalVector(fine,&xnew);CHKERRQ(ierr); 3641 ierr = MatInterpolate(interp,x,xnew);CHKERRQ(ierr); 3642 ierr = DMInterpolate(snes->dm,interp,fine);CHKERRQ(ierr); 3643 ierr = MatDestroy(&interp);CHKERRQ(ierr); 3644 x = xnew; 3645 3646 ierr = SNESReset(snes);CHKERRQ(ierr); 3647 ierr = SNESSetDM(snes,fine);CHKERRQ(ierr); 3648 ierr = DMDestroy(&fine);CHKERRQ(ierr); 3649 ierr = PetscViewerASCIIPopTab(PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)snes)));CHKERRQ(ierr); 3650 } 3651 } 3652 /* monitoring and viewing */ 3653 ierr = PetscOptionsGetViewer(PetscObjectComm((PetscObject)snes),((PetscObject)snes)->prefix,"-snes_view",&viewer,&format,&flg);CHKERRQ(ierr); 3654 if (flg && !PetscPreLoadingOn) { 3655 ierr = PetscViewerPushFormat(viewer,format);CHKERRQ(ierr); 3656 ierr = SNESView(snes,viewer);CHKERRQ(ierr); 3657 ierr = PetscViewerPopFormat(viewer);CHKERRQ(ierr); 3658 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 3659 } 3660 ierr = VecViewFromOptions(snes->vec_sol,"-snes_view_solution");CHKERRQ(ierr); 3661 3662 ierr = VecDestroy(&xcreated);CHKERRQ(ierr); 3663 #if defined(PETSC_HAVE_AMS) 3664 ierr = PetscObjectAMSBlock((PetscObject)snes);CHKERRQ(ierr); 3665 #endif 3666 PetscFunctionReturn(0); 3667 } 3668 3669 /* --------- Internal routines for SNES Package --------- */ 3670 3671 #undef __FUNCT__ 3672 #define __FUNCT__ "SNESSetType" 3673 /*@C 3674 SNESSetType - Sets the method for the nonlinear solver. 3675 3676 Collective on SNES 3677 3678 Input Parameters: 3679 + snes - the SNES context 3680 - type - a known method 3681 3682 Options Database Key: 3683 . -snes_type <type> - Sets the method; use -help for a list 3684 of available methods (for instance, newtonls or newtontr) 3685 3686 Notes: 3687 See "petsc/include/petscsnes.h" for available methods (for instance) 3688 + SNESNEWTONLS - Newton's method with line search 3689 (systems of nonlinear equations) 3690 . SNESNEWTONTR - Newton's method with trust region 3691 (systems of nonlinear equations) 3692 3693 Normally, it is best to use the SNESSetFromOptions() command and then 3694 set the SNES solver type from the options database rather than by using 3695 this routine. Using the options database provides the user with 3696 maximum flexibility in evaluating the many nonlinear solvers. 3697 The SNESSetType() routine is provided for those situations where it 3698 is necessary to set the nonlinear solver independently of the command 3699 line or options database. This might be the case, for example, when 3700 the choice of solver changes during the execution of the program, 3701 and the user's application is taking responsibility for choosing the 3702 appropriate method. 3703 3704 Level: intermediate 3705 3706 .keywords: SNES, set, type 3707 3708 .seealso: SNESType, SNESCreate() 3709 3710 @*/ 3711 PetscErrorCode SNESSetType(SNES snes,SNESType type) 3712 { 3713 PetscErrorCode ierr,(*r)(SNES); 3714 PetscBool match; 3715 3716 PetscFunctionBegin; 3717 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3718 PetscValidCharPointer(type,2); 3719 3720 ierr = PetscObjectTypeCompare((PetscObject)snes,type,&match);CHKERRQ(ierr); 3721 if (match) PetscFunctionReturn(0); 3722 3723 ierr = PetscFunctionListFind(PetscObjectComm((PetscObject)snes),SNESList,type,PETSC_TRUE,(void (**)(void)) &r);CHKERRQ(ierr); 3724 if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unable to find requested SNES type %s",type); 3725 /* Destroy the previous private SNES context */ 3726 if (snes->ops->destroy) { 3727 ierr = (*(snes)->ops->destroy)(snes);CHKERRQ(ierr); 3728 snes->ops->destroy = NULL; 3729 } 3730 /* Reinitialize function pointers in SNESOps structure */ 3731 snes->ops->setup = 0; 3732 snes->ops->solve = 0; 3733 snes->ops->view = 0; 3734 snes->ops->setfromoptions = 0; 3735 snes->ops->destroy = 0; 3736 /* Call the SNESCreate_XXX routine for this particular Nonlinear solver */ 3737 snes->setupcalled = PETSC_FALSE; 3738 3739 ierr = PetscObjectChangeTypeName((PetscObject)snes,type);CHKERRQ(ierr); 3740 ierr = (*r)(snes);CHKERRQ(ierr); 3741 #if defined(PETSC_HAVE_AMS) 3742 if (PetscAMSPublishAll) { 3743 ierr = PetscObjectAMSPublish((PetscObject)snes);CHKERRQ(ierr); 3744 } 3745 #endif 3746 PetscFunctionReturn(0); 3747 } 3748 3749 /* --------------------------------------------------------------------- */ 3750 #undef __FUNCT__ 3751 #define __FUNCT__ "SNESRegisterDestroy" 3752 /*@ 3753 SNESRegisterDestroy - Frees the list of nonlinear solvers that were 3754 registered by SNESRegisterDynamic(). 3755 3756 Not Collective 3757 3758 Level: advanced 3759 3760 .keywords: SNES, nonlinear, register, destroy 3761 3762 .seealso: SNESRegisterAll() 3763 @*/ 3764 PetscErrorCode SNESRegisterDestroy(void) 3765 { 3766 PetscErrorCode ierr; 3767 3768 PetscFunctionBegin; 3769 ierr = PetscFunctionListDestroy(&SNESList);CHKERRQ(ierr); 3770 3771 SNESRegisterAllCalled = PETSC_FALSE; 3772 PetscFunctionReturn(0); 3773 } 3774 3775 #undef __FUNCT__ 3776 #define __FUNCT__ "SNESGetType" 3777 /*@C 3778 SNESGetType - Gets the SNES method type and name (as a string). 3779 3780 Not Collective 3781 3782 Input Parameter: 3783 . snes - nonlinear solver context 3784 3785 Output Parameter: 3786 . type - SNES method (a character string) 3787 3788 Level: intermediate 3789 3790 .keywords: SNES, nonlinear, get, type, name 3791 @*/ 3792 PetscErrorCode SNESGetType(SNES snes,SNESType *type) 3793 { 3794 PetscFunctionBegin; 3795 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3796 PetscValidPointer(type,2); 3797 *type = ((PetscObject)snes)->type_name; 3798 PetscFunctionReturn(0); 3799 } 3800 3801 #undef __FUNCT__ 3802 #define __FUNCT__ "SNESGetSolution" 3803 /*@ 3804 SNESGetSolution - Returns the vector where the approximate solution is 3805 stored. This is the fine grid solution when using SNESSetGridSequence(). 3806 3807 Not Collective, but Vec is parallel if SNES is parallel 3808 3809 Input Parameter: 3810 . snes - the SNES context 3811 3812 Output Parameter: 3813 . x - the solution 3814 3815 Level: intermediate 3816 3817 .keywords: SNES, nonlinear, get, solution 3818 3819 .seealso: SNESGetSolutionUpdate(), SNESGetFunction() 3820 @*/ 3821 PetscErrorCode SNESGetSolution(SNES snes,Vec *x) 3822 { 3823 PetscFunctionBegin; 3824 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3825 PetscValidPointer(x,2); 3826 *x = snes->vec_sol; 3827 PetscFunctionReturn(0); 3828 } 3829 3830 #undef __FUNCT__ 3831 #define __FUNCT__ "SNESGetSolutionUpdate" 3832 /*@ 3833 SNESGetSolutionUpdate - Returns the vector where the solution update is 3834 stored. 3835 3836 Not Collective, but Vec is parallel if SNES is parallel 3837 3838 Input Parameter: 3839 . snes - the SNES context 3840 3841 Output Parameter: 3842 . x - the solution update 3843 3844 Level: advanced 3845 3846 .keywords: SNES, nonlinear, get, solution, update 3847 3848 .seealso: SNESGetSolution(), SNESGetFunction() 3849 @*/ 3850 PetscErrorCode SNESGetSolutionUpdate(SNES snes,Vec *x) 3851 { 3852 PetscFunctionBegin; 3853 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3854 PetscValidPointer(x,2); 3855 *x = snes->vec_sol_update; 3856 PetscFunctionReturn(0); 3857 } 3858 3859 #undef __FUNCT__ 3860 #define __FUNCT__ "SNESGetFunction" 3861 /*@C 3862 SNESGetFunction - Returns the vector where the function is stored. 3863 3864 Not Collective, but Vec is parallel if SNES is parallel. Collective if Vec is requested, but has not been created yet. 3865 3866 Input Parameter: 3867 . snes - the SNES context 3868 3869 Output Parameter: 3870 + r - the vector that is used to store residuals (or NULL if you don't want it) 3871 . SNESFunction- the function (or NULL if you don't want it) 3872 - ctx - the function context (or NULL if you don't want it) 3873 3874 Level: advanced 3875 3876 .keywords: SNES, nonlinear, get, function 3877 3878 .seealso: SNESSetFunction(), SNESGetSolution(), SNESFunction 3879 @*/ 3880 PetscErrorCode SNESGetFunction(SNES snes,Vec *r,PetscErrorCode (**SNESFunction)(SNES,Vec,Vec,void*),void **ctx) 3881 { 3882 PetscErrorCode ierr; 3883 DM dm; 3884 3885 PetscFunctionBegin; 3886 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3887 if (r) { 3888 if (!snes->vec_func) { 3889 if (snes->vec_rhs) { 3890 ierr = VecDuplicate(snes->vec_rhs,&snes->vec_func);CHKERRQ(ierr); 3891 } else if (snes->vec_sol) { 3892 ierr = VecDuplicate(snes->vec_sol,&snes->vec_func);CHKERRQ(ierr); 3893 } else if (snes->dm) { 3894 ierr = DMCreateGlobalVector(snes->dm,&snes->vec_func);CHKERRQ(ierr); 3895 } 3896 } 3897 *r = snes->vec_func; 3898 } 3899 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 3900 ierr = DMSNESGetFunction(dm,SNESFunction,ctx);CHKERRQ(ierr); 3901 PetscFunctionReturn(0); 3902 } 3903 3904 /*@C 3905 SNESGetGS - Returns the GS function and context. 3906 3907 Input Parameter: 3908 . snes - the SNES context 3909 3910 Output Parameter: 3911 + SNESGSFunction - the function (or NULL) 3912 - ctx - the function context (or NULL) 3913 3914 Level: advanced 3915 3916 .keywords: SNES, nonlinear, get, function 3917 3918 .seealso: SNESSetGS(), SNESGetFunction() 3919 @*/ 3920 3921 #undef __FUNCT__ 3922 #define __FUNCT__ "SNESGetGS" 3923 PetscErrorCode SNESGetGS (SNES snes, PetscErrorCode (**SNESGSFunction)(SNES, Vec, Vec, void*), void ** ctx) 3924 { 3925 PetscErrorCode ierr; 3926 DM dm; 3927 3928 PetscFunctionBegin; 3929 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3930 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 3931 ierr = DMSNESGetGS(dm,SNESGSFunction,ctx);CHKERRQ(ierr); 3932 PetscFunctionReturn(0); 3933 } 3934 3935 #undef __FUNCT__ 3936 #define __FUNCT__ "SNESSetOptionsPrefix" 3937 /*@C 3938 SNESSetOptionsPrefix - Sets the prefix used for searching for all 3939 SNES options in the database. 3940 3941 Logically Collective on SNES 3942 3943 Input Parameter: 3944 + snes - the SNES context 3945 - prefix - the prefix to prepend to all option names 3946 3947 Notes: 3948 A hyphen (-) must NOT be given at the beginning of the prefix name. 3949 The first character of all runtime options is AUTOMATICALLY the hyphen. 3950 3951 Level: advanced 3952 3953 .keywords: SNES, set, options, prefix, database 3954 3955 .seealso: SNESSetFromOptions() 3956 @*/ 3957 PetscErrorCode SNESSetOptionsPrefix(SNES snes,const char prefix[]) 3958 { 3959 PetscErrorCode ierr; 3960 3961 PetscFunctionBegin; 3962 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3963 ierr = PetscObjectSetOptionsPrefix((PetscObject)snes,prefix);CHKERRQ(ierr); 3964 if (!snes->ksp) {ierr = SNESGetKSP(snes,&snes->ksp);CHKERRQ(ierr);} 3965 if (snes->linesearch) { 3966 ierr = SNESGetSNESLineSearch(snes,&snes->linesearch);CHKERRQ(ierr); 3967 ierr = PetscObjectSetOptionsPrefix((PetscObject)snes->linesearch,prefix);CHKERRQ(ierr); 3968 } 3969 ierr = KSPSetOptionsPrefix(snes->ksp,prefix);CHKERRQ(ierr); 3970 PetscFunctionReturn(0); 3971 } 3972 3973 #undef __FUNCT__ 3974 #define __FUNCT__ "SNESAppendOptionsPrefix" 3975 /*@C 3976 SNESAppendOptionsPrefix - Appends to the prefix used for searching for all 3977 SNES options in the database. 3978 3979 Logically Collective on SNES 3980 3981 Input Parameters: 3982 + snes - the SNES context 3983 - prefix - the prefix to prepend to all option names 3984 3985 Notes: 3986 A hyphen (-) must NOT be given at the beginning of the prefix name. 3987 The first character of all runtime options is AUTOMATICALLY the hyphen. 3988 3989 Level: advanced 3990 3991 .keywords: SNES, append, options, prefix, database 3992 3993 .seealso: SNESGetOptionsPrefix() 3994 @*/ 3995 PetscErrorCode SNESAppendOptionsPrefix(SNES snes,const char prefix[]) 3996 { 3997 PetscErrorCode ierr; 3998 3999 PetscFunctionBegin; 4000 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4001 ierr = PetscObjectAppendOptionsPrefix((PetscObject)snes,prefix);CHKERRQ(ierr); 4002 if (!snes->ksp) {ierr = SNESGetKSP(snes,&snes->ksp);CHKERRQ(ierr);} 4003 if (snes->linesearch) { 4004 ierr = SNESGetSNESLineSearch(snes,&snes->linesearch);CHKERRQ(ierr); 4005 ierr = PetscObjectAppendOptionsPrefix((PetscObject)snes->linesearch,prefix);CHKERRQ(ierr); 4006 } 4007 ierr = KSPAppendOptionsPrefix(snes->ksp,prefix);CHKERRQ(ierr); 4008 PetscFunctionReturn(0); 4009 } 4010 4011 #undef __FUNCT__ 4012 #define __FUNCT__ "SNESGetOptionsPrefix" 4013 /*@C 4014 SNESGetOptionsPrefix - Sets the prefix used for searching for all 4015 SNES options in the database. 4016 4017 Not Collective 4018 4019 Input Parameter: 4020 . snes - the SNES context 4021 4022 Output Parameter: 4023 . prefix - pointer to the prefix string used 4024 4025 Notes: On the fortran side, the user should pass in a string 'prefix' of 4026 sufficient length to hold the prefix. 4027 4028 Level: advanced 4029 4030 .keywords: SNES, get, options, prefix, database 4031 4032 .seealso: SNESAppendOptionsPrefix() 4033 @*/ 4034 PetscErrorCode SNESGetOptionsPrefix(SNES snes,const char *prefix[]) 4035 { 4036 PetscErrorCode ierr; 4037 4038 PetscFunctionBegin; 4039 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4040 ierr = PetscObjectGetOptionsPrefix((PetscObject)snes,prefix);CHKERRQ(ierr); 4041 PetscFunctionReturn(0); 4042 } 4043 4044 4045 #undef __FUNCT__ 4046 #define __FUNCT__ "SNESRegister" 4047 /*@C 4048 SNESRegister - See SNESRegisterDynamic() 4049 4050 Level: advanced 4051 @*/ 4052 PetscErrorCode SNESRegister(const char sname[],const char path[],const char name[],PetscErrorCode (*function)(SNES)) 4053 { 4054 char fullname[PETSC_MAX_PATH_LEN]; 4055 PetscErrorCode ierr; 4056 4057 PetscFunctionBegin; 4058 ierr = PetscFunctionListConcat(path,name,fullname);CHKERRQ(ierr); 4059 ierr = PetscFunctionListAdd(PETSC_COMM_WORLD,&SNESList,sname,fullname,(void (*)(void))function);CHKERRQ(ierr); 4060 PetscFunctionReturn(0); 4061 } 4062 4063 #undef __FUNCT__ 4064 #define __FUNCT__ "SNESTestLocalMin" 4065 PetscErrorCode SNESTestLocalMin(SNES snes) 4066 { 4067 PetscErrorCode ierr; 4068 PetscInt N,i,j; 4069 Vec u,uh,fh; 4070 PetscScalar value; 4071 PetscReal norm; 4072 4073 PetscFunctionBegin; 4074 ierr = SNESGetSolution(snes,&u);CHKERRQ(ierr); 4075 ierr = VecDuplicate(u,&uh);CHKERRQ(ierr); 4076 ierr = VecDuplicate(u,&fh);CHKERRQ(ierr); 4077 4078 /* currently only works for sequential */ 4079 ierr = PetscPrintf(PETSC_COMM_WORLD,"Testing FormFunction() for local min\n");CHKERRQ(ierr); 4080 ierr = VecGetSize(u,&N);CHKERRQ(ierr); 4081 for (i=0; i<N; i++) { 4082 ierr = VecCopy(u,uh);CHKERRQ(ierr); 4083 ierr = PetscPrintf(PETSC_COMM_WORLD,"i = %D\n",i);CHKERRQ(ierr); 4084 for (j=-10; j<11; j++) { 4085 value = PetscSign(j)*exp(PetscAbs(j)-10.0); 4086 ierr = VecSetValue(uh,i,value,ADD_VALUES);CHKERRQ(ierr); 4087 ierr = SNESComputeFunction(snes,uh,fh);CHKERRQ(ierr); 4088 ierr = VecNorm(fh,NORM_2,&norm);CHKERRQ(ierr); 4089 ierr = PetscPrintf(PETSC_COMM_WORLD," j norm %D %18.16e\n",j,norm);CHKERRQ(ierr); 4090 value = -value; 4091 ierr = VecSetValue(uh,i,value,ADD_VALUES);CHKERRQ(ierr); 4092 } 4093 } 4094 ierr = VecDestroy(&uh);CHKERRQ(ierr); 4095 ierr = VecDestroy(&fh);CHKERRQ(ierr); 4096 PetscFunctionReturn(0); 4097 } 4098 4099 #undef __FUNCT__ 4100 #define __FUNCT__ "SNESKSPSetUseEW" 4101 /*@ 4102 SNESKSPSetUseEW - Sets SNES use Eisenstat-Walker method for 4103 computing relative tolerance for linear solvers within an inexact 4104 Newton method. 4105 4106 Logically Collective on SNES 4107 4108 Input Parameters: 4109 + snes - SNES context 4110 - flag - PETSC_TRUE or PETSC_FALSE 4111 4112 Options Database: 4113 + -snes_ksp_ew - use Eisenstat-Walker method for determining linear system convergence 4114 . -snes_ksp_ew_version ver - version of Eisenstat-Walker method 4115 . -snes_ksp_ew_rtol0 <rtol0> - Sets rtol0 4116 . -snes_ksp_ew_rtolmax <rtolmax> - Sets rtolmax 4117 . -snes_ksp_ew_gamma <gamma> - Sets gamma 4118 . -snes_ksp_ew_alpha <alpha> - Sets alpha 4119 . -snes_ksp_ew_alpha2 <alpha2> - Sets alpha2 4120 - -snes_ksp_ew_threshold <threshold> - Sets threshold 4121 4122 Notes: 4123 Currently, the default is to use a constant relative tolerance for 4124 the inner linear solvers. Alternatively, one can use the 4125 Eisenstat-Walker method, where the relative convergence tolerance 4126 is reset at each Newton iteration according progress of the nonlinear 4127 solver. 4128 4129 Level: advanced 4130 4131 Reference: 4132 S. C. Eisenstat and H. F. Walker, "Choosing the forcing terms in an 4133 inexact Newton method", SISC 17 (1), pp.16-32, 1996. 4134 4135 .keywords: SNES, KSP, Eisenstat, Walker, convergence, test, inexact, Newton 4136 4137 .seealso: SNESKSPGetUseEW(), SNESKSPGetParametersEW(), SNESKSPSetParametersEW() 4138 @*/ 4139 PetscErrorCode SNESKSPSetUseEW(SNES snes,PetscBool flag) 4140 { 4141 PetscFunctionBegin; 4142 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4143 PetscValidLogicalCollectiveBool(snes,flag,2); 4144 snes->ksp_ewconv = flag; 4145 PetscFunctionReturn(0); 4146 } 4147 4148 #undef __FUNCT__ 4149 #define __FUNCT__ "SNESKSPGetUseEW" 4150 /*@ 4151 SNESKSPGetUseEW - Gets if SNES is using Eisenstat-Walker method 4152 for computing relative tolerance for linear solvers within an 4153 inexact Newton method. 4154 4155 Not Collective 4156 4157 Input Parameter: 4158 . snes - SNES context 4159 4160 Output Parameter: 4161 . flag - PETSC_TRUE or PETSC_FALSE 4162 4163 Notes: 4164 Currently, the default is to use a constant relative tolerance for 4165 the inner linear solvers. Alternatively, one can use the 4166 Eisenstat-Walker method, where the relative convergence tolerance 4167 is reset at each Newton iteration according progress of the nonlinear 4168 solver. 4169 4170 Level: advanced 4171 4172 Reference: 4173 S. C. Eisenstat and H. F. Walker, "Choosing the forcing terms in an 4174 inexact Newton method", SISC 17 (1), pp.16-32, 1996. 4175 4176 .keywords: SNES, KSP, Eisenstat, Walker, convergence, test, inexact, Newton 4177 4178 .seealso: SNESKSPSetUseEW(), SNESKSPGetParametersEW(), SNESKSPSetParametersEW() 4179 @*/ 4180 PetscErrorCode SNESKSPGetUseEW(SNES snes, PetscBool *flag) 4181 { 4182 PetscFunctionBegin; 4183 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4184 PetscValidPointer(flag,2); 4185 *flag = snes->ksp_ewconv; 4186 PetscFunctionReturn(0); 4187 } 4188 4189 #undef __FUNCT__ 4190 #define __FUNCT__ "SNESKSPSetParametersEW" 4191 /*@ 4192 SNESKSPSetParametersEW - Sets parameters for Eisenstat-Walker 4193 convergence criteria for the linear solvers within an inexact 4194 Newton method. 4195 4196 Logically Collective on SNES 4197 4198 Input Parameters: 4199 + snes - SNES context 4200 . version - version 1, 2 (default is 2) or 3 4201 . rtol_0 - initial relative tolerance (0 <= rtol_0 < 1) 4202 . rtol_max - maximum relative tolerance (0 <= rtol_max < 1) 4203 . gamma - multiplicative factor for version 2 rtol computation 4204 (0 <= gamma2 <= 1) 4205 . alpha - power for version 2 rtol computation (1 < alpha <= 2) 4206 . alpha2 - power for safeguard 4207 - threshold - threshold for imposing safeguard (0 < threshold < 1) 4208 4209 Note: 4210 Version 3 was contributed by Luis Chacon, June 2006. 4211 4212 Use PETSC_DEFAULT to retain the default for any of the parameters. 4213 4214 Level: advanced 4215 4216 Reference: 4217 S. C. Eisenstat and H. F. Walker, "Choosing the forcing terms in an 4218 inexact Newton method", Utah State University Math. Stat. Dept. Res. 4219 Report 6/94/75, June, 1994, to appear in SIAM J. Sci. Comput. 4220 4221 .keywords: SNES, KSP, Eisenstat, Walker, set, parameters 4222 4223 .seealso: SNESKSPSetUseEW(), SNESKSPGetUseEW(), SNESKSPGetParametersEW() 4224 @*/ 4225 PetscErrorCode SNESKSPSetParametersEW(SNES snes,PetscInt version,PetscReal rtol_0,PetscReal rtol_max,PetscReal gamma,PetscReal alpha,PetscReal alpha2,PetscReal threshold) 4226 { 4227 SNESKSPEW *kctx; 4228 4229 PetscFunctionBegin; 4230 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4231 kctx = (SNESKSPEW*)snes->kspconvctx; 4232 if (!kctx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"No Eisenstat-Walker context existing"); 4233 PetscValidLogicalCollectiveInt(snes,version,2); 4234 PetscValidLogicalCollectiveReal(snes,rtol_0,3); 4235 PetscValidLogicalCollectiveReal(snes,rtol_max,4); 4236 PetscValidLogicalCollectiveReal(snes,gamma,5); 4237 PetscValidLogicalCollectiveReal(snes,alpha,6); 4238 PetscValidLogicalCollectiveReal(snes,alpha2,7); 4239 PetscValidLogicalCollectiveReal(snes,threshold,8); 4240 4241 if (version != PETSC_DEFAULT) kctx->version = version; 4242 if (rtol_0 != PETSC_DEFAULT) kctx->rtol_0 = rtol_0; 4243 if (rtol_max != PETSC_DEFAULT) kctx->rtol_max = rtol_max; 4244 if (gamma != PETSC_DEFAULT) kctx->gamma = gamma; 4245 if (alpha != PETSC_DEFAULT) kctx->alpha = alpha; 4246 if (alpha2 != PETSC_DEFAULT) kctx->alpha2 = alpha2; 4247 if (threshold != PETSC_DEFAULT) kctx->threshold = threshold; 4248 4249 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); 4250 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); 4251 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); 4252 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); 4253 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); 4254 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); 4255 PetscFunctionReturn(0); 4256 } 4257 4258 #undef __FUNCT__ 4259 #define __FUNCT__ "SNESKSPGetParametersEW" 4260 /*@ 4261 SNESKSPGetParametersEW - Gets parameters for Eisenstat-Walker 4262 convergence criteria for the linear solvers within an inexact 4263 Newton method. 4264 4265 Not Collective 4266 4267 Input Parameters: 4268 snes - SNES context 4269 4270 Output Parameters: 4271 + version - version 1, 2 (default is 2) or 3 4272 . rtol_0 - initial relative tolerance (0 <= rtol_0 < 1) 4273 . rtol_max - maximum relative tolerance (0 <= rtol_max < 1) 4274 . gamma - multiplicative factor for version 2 rtol computation (0 <= gamma2 <= 1) 4275 . alpha - power for version 2 rtol computation (1 < alpha <= 2) 4276 . alpha2 - power for safeguard 4277 - threshold - threshold for imposing safeguard (0 < threshold < 1) 4278 4279 Level: advanced 4280 4281 .keywords: SNES, KSP, Eisenstat, Walker, get, parameters 4282 4283 .seealso: SNESKSPSetUseEW(), SNESKSPGetUseEW(), SNESKSPSetParametersEW() 4284 @*/ 4285 PetscErrorCode SNESKSPGetParametersEW(SNES snes,PetscInt *version,PetscReal *rtol_0,PetscReal *rtol_max,PetscReal *gamma,PetscReal *alpha,PetscReal *alpha2,PetscReal *threshold) 4286 { 4287 SNESKSPEW *kctx; 4288 4289 PetscFunctionBegin; 4290 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4291 kctx = (SNESKSPEW*)snes->kspconvctx; 4292 if (!kctx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"No Eisenstat-Walker context existing"); 4293 if (version) *version = kctx->version; 4294 if (rtol_0) *rtol_0 = kctx->rtol_0; 4295 if (rtol_max) *rtol_max = kctx->rtol_max; 4296 if (gamma) *gamma = kctx->gamma; 4297 if (alpha) *alpha = kctx->alpha; 4298 if (alpha2) *alpha2 = kctx->alpha2; 4299 if (threshold) *threshold = kctx->threshold; 4300 PetscFunctionReturn(0); 4301 } 4302 4303 #undef __FUNCT__ 4304 #define __FUNCT__ "SNESKSPEW_PreSolve" 4305 PetscErrorCode SNESKSPEW_PreSolve(KSP ksp, Vec b, Vec x, SNES snes) 4306 { 4307 PetscErrorCode ierr; 4308 SNESKSPEW *kctx = (SNESKSPEW*)snes->kspconvctx; 4309 PetscReal rtol = PETSC_DEFAULT,stol; 4310 4311 PetscFunctionBegin; 4312 if (!snes->ksp_ewconv) PetscFunctionReturn(0); 4313 if (!snes->iter) rtol = kctx->rtol_0; /* first time in, so use the original user rtol */ 4314 else { 4315 if (kctx->version == 1) { 4316 rtol = (snes->norm - kctx->lresid_last)/kctx->norm_last; 4317 if (rtol < 0.0) rtol = -rtol; 4318 stol = PetscPowReal(kctx->rtol_last,kctx->alpha2); 4319 if (stol > kctx->threshold) rtol = PetscMax(rtol,stol); 4320 } else if (kctx->version == 2) { 4321 rtol = kctx->gamma * PetscPowReal(snes->norm/kctx->norm_last,kctx->alpha); 4322 stol = kctx->gamma * PetscPowReal(kctx->rtol_last,kctx->alpha); 4323 if (stol > kctx->threshold) rtol = PetscMax(rtol,stol); 4324 } else if (kctx->version == 3) { /* contributed by Luis Chacon, June 2006. */ 4325 rtol = kctx->gamma * PetscPowReal(snes->norm/kctx->norm_last,kctx->alpha); 4326 /* safeguard: avoid sharp decrease of rtol */ 4327 stol = kctx->gamma*PetscPowReal(kctx->rtol_last,kctx->alpha); 4328 stol = PetscMax(rtol,stol); 4329 rtol = PetscMin(kctx->rtol_0,stol); 4330 /* safeguard: avoid oversolving */ 4331 stol = kctx->gamma*(snes->ttol)/snes->norm; 4332 stol = PetscMax(rtol,stol); 4333 rtol = PetscMin(kctx->rtol_0,stol); 4334 } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Only versions 1, 2 or 3 are supported: %D",kctx->version); 4335 } 4336 /* safeguard: avoid rtol greater than one */ 4337 rtol = PetscMin(rtol,kctx->rtol_max); 4338 ierr = KSPSetTolerances(ksp,rtol,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);CHKERRQ(ierr); 4339 ierr = PetscInfo3(snes,"iter %D, Eisenstat-Walker (version %D) KSP rtol=%G\n",snes->iter,kctx->version,rtol);CHKERRQ(ierr); 4340 PetscFunctionReturn(0); 4341 } 4342 4343 #undef __FUNCT__ 4344 #define __FUNCT__ "SNESKSPEW_PostSolve" 4345 PetscErrorCode SNESKSPEW_PostSolve(KSP ksp, Vec b, Vec x, SNES snes) 4346 { 4347 PetscErrorCode ierr; 4348 SNESKSPEW *kctx = (SNESKSPEW*)snes->kspconvctx; 4349 PCSide pcside; 4350 Vec lres; 4351 4352 PetscFunctionBegin; 4353 if (!snes->ksp_ewconv) PetscFunctionReturn(0); 4354 ierr = KSPGetTolerances(ksp,&kctx->rtol_last,0,0,0);CHKERRQ(ierr); 4355 ierr = SNESGetFunctionNorm(snes,&kctx->norm_last);CHKERRQ(ierr); 4356 if (kctx->version == 1) { 4357 ierr = KSPGetPCSide(ksp,&pcside);CHKERRQ(ierr); 4358 if (pcside == PC_RIGHT) { /* XXX Should we also test KSP_UNPRECONDITIONED_NORM ? */ 4359 /* KSP residual is true linear residual */ 4360 ierr = KSPGetResidualNorm(ksp,&kctx->lresid_last);CHKERRQ(ierr); 4361 } else { 4362 /* KSP residual is preconditioned residual */ 4363 /* compute true linear residual norm */ 4364 ierr = VecDuplicate(b,&lres);CHKERRQ(ierr); 4365 ierr = MatMult(snes->jacobian,x,lres);CHKERRQ(ierr); 4366 ierr = VecAYPX(lres,-1.0,b);CHKERRQ(ierr); 4367 ierr = VecNorm(lres,NORM_2,&kctx->lresid_last);CHKERRQ(ierr); 4368 ierr = VecDestroy(&lres);CHKERRQ(ierr); 4369 } 4370 } 4371 PetscFunctionReturn(0); 4372 } 4373 4374 #undef __FUNCT__ 4375 #define __FUNCT__ "SNESGetKSP" 4376 /*@ 4377 SNESGetKSP - Returns the KSP context for a SNES solver. 4378 4379 Not Collective, but if SNES object is parallel, then KSP object is parallel 4380 4381 Input Parameter: 4382 . snes - the SNES context 4383 4384 Output Parameter: 4385 . ksp - the KSP context 4386 4387 Notes: 4388 The user can then directly manipulate the KSP context to set various 4389 options, etc. Likewise, the user can then extract and manipulate the 4390 PC contexts as well. 4391 4392 Level: beginner 4393 4394 .keywords: SNES, nonlinear, get, KSP, context 4395 4396 .seealso: KSPGetPC(), SNESCreate(), KSPCreate(), SNESSetKSP() 4397 @*/ 4398 PetscErrorCode SNESGetKSP(SNES snes,KSP *ksp) 4399 { 4400 PetscErrorCode ierr; 4401 4402 PetscFunctionBegin; 4403 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4404 PetscValidPointer(ksp,2); 4405 4406 if (!snes->ksp) { 4407 ierr = KSPCreate(PetscObjectComm((PetscObject)snes),&snes->ksp);CHKERRQ(ierr); 4408 ierr = PetscObjectIncrementTabLevel((PetscObject)snes->ksp,(PetscObject)snes,1);CHKERRQ(ierr); 4409 ierr = PetscLogObjectParent(snes,snes->ksp);CHKERRQ(ierr); 4410 4411 ierr = KSPSetPreSolve(snes->ksp,(PetscErrorCode (*)(KSP,Vec,Vec,void*))SNESKSPEW_PreSolve,snes);CHKERRQ(ierr); 4412 ierr = KSPSetPostSolve(snes->ksp,(PetscErrorCode (*)(KSP,Vec,Vec,void*))SNESKSPEW_PostSolve,snes);CHKERRQ(ierr); 4413 } 4414 *ksp = snes->ksp; 4415 PetscFunctionReturn(0); 4416 } 4417 4418 4419 #include <petsc-private/dmimpl.h> 4420 #undef __FUNCT__ 4421 #define __FUNCT__ "SNESSetDM" 4422 /*@ 4423 SNESSetDM - Sets the DM that may be used by some preconditioners 4424 4425 Logically Collective on SNES 4426 4427 Input Parameters: 4428 + snes - the preconditioner context 4429 - dm - the dm 4430 4431 Level: intermediate 4432 4433 .seealso: SNESGetDM(), KSPSetDM(), KSPGetDM() 4434 @*/ 4435 PetscErrorCode SNESSetDM(SNES snes,DM dm) 4436 { 4437 PetscErrorCode ierr; 4438 KSP ksp; 4439 DMSNES sdm; 4440 4441 PetscFunctionBegin; 4442 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4443 if (dm) {ierr = PetscObjectReference((PetscObject)dm);CHKERRQ(ierr);} 4444 if (snes->dm) { /* Move the DMSNES context over to the new DM unless the new DM already has one */ 4445 if (snes->dm->dmsnes && snes->dmAuto && !dm->dmsnes) { 4446 ierr = DMCopyDMSNES(snes->dm,dm);CHKERRQ(ierr); 4447 ierr = DMGetDMSNES(snes->dm,&sdm);CHKERRQ(ierr); 4448 if (sdm->originaldm == snes->dm) sdm->originaldm = dm; /* Grant write privileges to the replacement DM */ 4449 } 4450 ierr = DMDestroy(&snes->dm);CHKERRQ(ierr); 4451 } 4452 snes->dm = dm; 4453 snes->dmAuto = PETSC_FALSE; 4454 4455 ierr = SNESGetKSP(snes,&ksp);CHKERRQ(ierr); 4456 ierr = KSPSetDM(ksp,dm);CHKERRQ(ierr); 4457 ierr = KSPSetDMActive(ksp,PETSC_FALSE);CHKERRQ(ierr); 4458 if (snes->pc) { 4459 ierr = SNESSetDM(snes->pc, snes->dm);CHKERRQ(ierr); 4460 ierr = SNESSetPCSide(snes,snes->pcside);CHKERRQ(ierr); 4461 } 4462 PetscFunctionReturn(0); 4463 } 4464 4465 #undef __FUNCT__ 4466 #define __FUNCT__ "SNESGetDM" 4467 /*@ 4468 SNESGetDM - Gets the DM that may be used by some preconditioners 4469 4470 Not Collective but DM obtained is parallel on SNES 4471 4472 Input Parameter: 4473 . snes - the preconditioner context 4474 4475 Output Parameter: 4476 . dm - the dm 4477 4478 Level: intermediate 4479 4480 .seealso: SNESSetDM(), KSPSetDM(), KSPGetDM() 4481 @*/ 4482 PetscErrorCode SNESGetDM(SNES snes,DM *dm) 4483 { 4484 PetscErrorCode ierr; 4485 4486 PetscFunctionBegin; 4487 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4488 if (!snes->dm) { 4489 ierr = DMShellCreate(PetscObjectComm((PetscObject)snes),&snes->dm);CHKERRQ(ierr); 4490 snes->dmAuto = PETSC_TRUE; 4491 } 4492 *dm = snes->dm; 4493 PetscFunctionReturn(0); 4494 } 4495 4496 #undef __FUNCT__ 4497 #define __FUNCT__ "SNESSetPC" 4498 /*@ 4499 SNESSetPC - Sets the nonlinear preconditioner to be used. 4500 4501 Collective on SNES 4502 4503 Input Parameters: 4504 + snes - iterative context obtained from SNESCreate() 4505 - pc - the preconditioner object 4506 4507 Notes: 4508 Use SNESGetPC() to retrieve the preconditioner context (for example, 4509 to configure it using the API). 4510 4511 Level: developer 4512 4513 .keywords: SNES, set, precondition 4514 .seealso: SNESGetPC() 4515 @*/ 4516 PetscErrorCode SNESSetPC(SNES snes, SNES pc) 4517 { 4518 PetscErrorCode ierr; 4519 4520 PetscFunctionBegin; 4521 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); 4522 PetscValidHeaderSpecific(pc, SNES_CLASSID, 2); 4523 PetscCheckSameComm(snes, 1, pc, 2); 4524 ierr = PetscObjectReference((PetscObject) pc);CHKERRQ(ierr); 4525 ierr = SNESDestroy(&snes->pc);CHKERRQ(ierr); 4526 snes->pc = pc; 4527 ierr = PetscLogObjectParent(snes, snes->pc);CHKERRQ(ierr); 4528 PetscFunctionReturn(0); 4529 } 4530 4531 #undef __FUNCT__ 4532 #define __FUNCT__ "SNESGetPC" 4533 /*@ 4534 SNESGetPC - Returns a pointer to the nonlinear preconditioning context set with SNESSetPC(). 4535 4536 Not Collective 4537 4538 Input Parameter: 4539 . snes - iterative context obtained from SNESCreate() 4540 4541 Output Parameter: 4542 . pc - preconditioner context 4543 4544 Level: developer 4545 4546 .keywords: SNES, get, preconditioner 4547 .seealso: SNESSetPC() 4548 @*/ 4549 PetscErrorCode SNESGetPC(SNES snes, SNES *pc) 4550 { 4551 PetscErrorCode ierr; 4552 const char *optionsprefix; 4553 4554 PetscFunctionBegin; 4555 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); 4556 PetscValidPointer(pc, 2); 4557 if (!snes->pc) { 4558 ierr = SNESCreate(PetscObjectComm((PetscObject)snes),&snes->pc);CHKERRQ(ierr); 4559 ierr = PetscObjectIncrementTabLevel((PetscObject)snes->pc,(PetscObject)snes,1);CHKERRQ(ierr); 4560 ierr = PetscLogObjectParent(snes,snes->pc);CHKERRQ(ierr); 4561 ierr = SNESGetOptionsPrefix(snes,&optionsprefix);CHKERRQ(ierr); 4562 ierr = SNESSetOptionsPrefix(snes->pc,optionsprefix);CHKERRQ(ierr); 4563 ierr = SNESAppendOptionsPrefix(snes->pc,"npc_");CHKERRQ(ierr); 4564 } 4565 *pc = snes->pc; 4566 PetscFunctionReturn(0); 4567 } 4568 4569 #undef __FUNCT__ 4570 #define __FUNCT__ "SNESSetPCSide" 4571 /*@ 4572 SNESSetPCSide - Sets the preconditioning side. 4573 4574 Logically Collective on SNES 4575 4576 Input Parameter: 4577 . snes - iterative context obtained from SNESCreate() 4578 4579 Output Parameter: 4580 . side - the preconditioning side, where side is one of 4581 .vb 4582 PC_LEFT - left preconditioning (default) 4583 PC_RIGHT - right preconditioning 4584 .ve 4585 4586 Options Database Keys: 4587 . -snes_pc_side <right,left> 4588 4589 Level: intermediate 4590 4591 .keywords: SNES, set, right, left, side, preconditioner, flag 4592 4593 .seealso: SNESGetPCSide(), KSPSetPCSide() 4594 @*/ 4595 PetscErrorCode SNESSetPCSide(SNES snes,PCSide side) 4596 { 4597 PetscFunctionBegin; 4598 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4599 PetscValidLogicalCollectiveEnum(snes,side,2); 4600 snes->pcside = side; 4601 PetscFunctionReturn(0); 4602 } 4603 4604 #undef __FUNCT__ 4605 #define __FUNCT__ "SNESGetPCSide" 4606 /*@ 4607 SNESGetPCSide - Gets the preconditioning side. 4608 4609 Not Collective 4610 4611 Input Parameter: 4612 . snes - iterative context obtained from SNESCreate() 4613 4614 Output Parameter: 4615 . side - the preconditioning side, where side is one of 4616 .vb 4617 PC_LEFT - left preconditioning (default) 4618 PC_RIGHT - right preconditioning 4619 .ve 4620 4621 Level: intermediate 4622 4623 .keywords: SNES, get, right, left, side, preconditioner, flag 4624 4625 .seealso: SNESSetPCSide(), KSPGetPCSide() 4626 @*/ 4627 PetscErrorCode SNESGetPCSide(SNES snes,PCSide *side) 4628 { 4629 PetscFunctionBegin; 4630 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4631 PetscValidPointer(side,2); 4632 *side = snes->pcside; 4633 PetscFunctionReturn(0); 4634 } 4635 4636 #undef __FUNCT__ 4637 #define __FUNCT__ "SNESSetSNESLineSearch" 4638 /*@ 4639 SNESSetSNESLineSearch - Sets the linesearch on the SNES instance. 4640 4641 Collective on SNES 4642 4643 Input Parameters: 4644 + snes - iterative context obtained from SNESCreate() 4645 - linesearch - the linesearch object 4646 4647 Notes: 4648 Use SNESGetSNESLineSearch() to retrieve the preconditioner context (for example, 4649 to configure it using the API). 4650 4651 Level: developer 4652 4653 .keywords: SNES, set, linesearch 4654 .seealso: SNESGetSNESLineSearch() 4655 @*/ 4656 PetscErrorCode SNESSetSNESLineSearch(SNES snes, SNESLineSearch linesearch) 4657 { 4658 PetscErrorCode ierr; 4659 4660 PetscFunctionBegin; 4661 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); 4662 PetscValidHeaderSpecific(linesearch, SNESLINESEARCH_CLASSID, 2); 4663 PetscCheckSameComm(snes, 1, linesearch, 2); 4664 ierr = PetscObjectReference((PetscObject) linesearch);CHKERRQ(ierr); 4665 ierr = SNESLineSearchDestroy(&snes->linesearch);CHKERRQ(ierr); 4666 4667 snes->linesearch = linesearch; 4668 4669 ierr = PetscLogObjectParent(snes, snes->linesearch);CHKERRQ(ierr); 4670 PetscFunctionReturn(0); 4671 } 4672 4673 #undef __FUNCT__ 4674 #define __FUNCT__ "SNESGetSNESLineSearch" 4675 /*@C 4676 SNESGetSNESLineSearch - Returns a pointer to the line search context set with SNESSetLineSearch() 4677 or creates a default line search instance associated with the SNES and returns it. 4678 4679 Not Collective 4680 4681 Input Parameter: 4682 . snes - iterative context obtained from SNESCreate() 4683 4684 Output Parameter: 4685 . linesearch - linesearch context 4686 4687 Level: developer 4688 4689 .keywords: SNES, get, linesearch 4690 .seealso: SNESSetSNESLineSearch() 4691 @*/ 4692 PetscErrorCode SNESGetSNESLineSearch(SNES snes, SNESLineSearch *linesearch) 4693 { 4694 PetscErrorCode ierr; 4695 const char *optionsprefix; 4696 4697 PetscFunctionBegin; 4698 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); 4699 PetscValidPointer(linesearch, 2); 4700 if (!snes->linesearch) { 4701 ierr = SNESGetOptionsPrefix(snes, &optionsprefix);CHKERRQ(ierr); 4702 ierr = SNESLineSearchCreate(PetscObjectComm((PetscObject)snes), &snes->linesearch);CHKERRQ(ierr); 4703 ierr = SNESLineSearchSetSNES(snes->linesearch, snes);CHKERRQ(ierr); 4704 ierr = SNESLineSearchAppendOptionsPrefix(snes->linesearch, optionsprefix);CHKERRQ(ierr); 4705 ierr = PetscObjectIncrementTabLevel((PetscObject) snes->linesearch, (PetscObject) snes, 1);CHKERRQ(ierr); 4706 ierr = PetscLogObjectParent(snes, snes->linesearch);CHKERRQ(ierr); 4707 } 4708 *linesearch = snes->linesearch; 4709 PetscFunctionReturn(0); 4710 } 4711 4712 #if defined(PETSC_HAVE_MATLAB_ENGINE) 4713 #include <mex.h> 4714 4715 typedef struct {char *funcname; mxArray *ctx;} SNESMatlabContext; 4716 4717 #undef __FUNCT__ 4718 #define __FUNCT__ "SNESComputeFunction_Matlab" 4719 /* 4720 SNESComputeFunction_Matlab - Calls the function that has been set with SNESSetFunctionMatlab(). 4721 4722 Collective on SNES 4723 4724 Input Parameters: 4725 + snes - the SNES context 4726 - x - input vector 4727 4728 Output Parameter: 4729 . y - function vector, as set by SNESSetFunction() 4730 4731 Notes: 4732 SNESComputeFunction() is typically used within nonlinear solvers 4733 implementations, so most users would not generally call this routine 4734 themselves. 4735 4736 Level: developer 4737 4738 .keywords: SNES, nonlinear, compute, function 4739 4740 .seealso: SNESSetFunction(), SNESGetFunction() 4741 */ 4742 PetscErrorCode SNESComputeFunction_Matlab(SNES snes,Vec x,Vec y, void *ctx) 4743 { 4744 PetscErrorCode ierr; 4745 SNESMatlabContext *sctx = (SNESMatlabContext*)ctx; 4746 int nlhs = 1,nrhs = 5; 4747 mxArray *plhs[1],*prhs[5]; 4748 long long int lx = 0,ly = 0,ls = 0; 4749 4750 PetscFunctionBegin; 4751 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4752 PetscValidHeaderSpecific(x,VEC_CLASSID,2); 4753 PetscValidHeaderSpecific(y,VEC_CLASSID,3); 4754 PetscCheckSameComm(snes,1,x,2); 4755 PetscCheckSameComm(snes,1,y,3); 4756 4757 /* call Matlab function in ctx with arguments x and y */ 4758 4759 ierr = PetscMemcpy(&ls,&snes,sizeof(snes));CHKERRQ(ierr); 4760 ierr = PetscMemcpy(&lx,&x,sizeof(x));CHKERRQ(ierr); 4761 ierr = PetscMemcpy(&ly,&y,sizeof(x));CHKERRQ(ierr); 4762 prhs[0] = mxCreateDoubleScalar((double)ls); 4763 prhs[1] = mxCreateDoubleScalar((double)lx); 4764 prhs[2] = mxCreateDoubleScalar((double)ly); 4765 prhs[3] = mxCreateString(sctx->funcname); 4766 prhs[4] = sctx->ctx; 4767 ierr = mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscSNESComputeFunctionInternal");CHKERRQ(ierr); 4768 ierr = mxGetScalar(plhs[0]);CHKERRQ(ierr); 4769 mxDestroyArray(prhs[0]); 4770 mxDestroyArray(prhs[1]); 4771 mxDestroyArray(prhs[2]); 4772 mxDestroyArray(prhs[3]); 4773 mxDestroyArray(plhs[0]); 4774 PetscFunctionReturn(0); 4775 } 4776 4777 #undef __FUNCT__ 4778 #define __FUNCT__ "SNESSetFunctionMatlab" 4779 /* 4780 SNESSetFunctionMatlab - Sets the function evaluation routine and function 4781 vector for use by the SNES routines in solving systems of nonlinear 4782 equations from MATLAB. Here the function is a string containing the name of a MATLAB function 4783 4784 Logically Collective on SNES 4785 4786 Input Parameters: 4787 + snes - the SNES context 4788 . r - vector to store function value 4789 - SNESFunction - function evaluation routine 4790 4791 Notes: 4792 The Newton-like methods typically solve linear systems of the form 4793 $ f'(x) x = -f(x), 4794 where f'(x) denotes the Jacobian matrix and f(x) is the function. 4795 4796 Level: beginner 4797 4798 .keywords: SNES, nonlinear, set, function 4799 4800 .seealso: SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESSetFunction() 4801 */ 4802 PetscErrorCode SNESSetFunctionMatlab(SNES snes,Vec r,const char *SNESFunction,mxArray *ctx) 4803 { 4804 PetscErrorCode ierr; 4805 SNESMatlabContext *sctx; 4806 4807 PetscFunctionBegin; 4808 /* currently sctx is memory bleed */ 4809 ierr = PetscMalloc(sizeof(SNESMatlabContext),&sctx);CHKERRQ(ierr); 4810 ierr = PetscStrallocpy(SNESFunction,&sctx->funcname);CHKERRQ(ierr); 4811 /* 4812 This should work, but it doesn't 4813 sctx->ctx = ctx; 4814 mexMakeArrayPersistent(sctx->ctx); 4815 */ 4816 sctx->ctx = mxDuplicateArray(ctx); 4817 ierr = SNESSetFunction(snes,r,SNESComputeFunction_Matlab,sctx);CHKERRQ(ierr); 4818 PetscFunctionReturn(0); 4819 } 4820 4821 #undef __FUNCT__ 4822 #define __FUNCT__ "SNESComputeJacobian_Matlab" 4823 /* 4824 SNESComputeJacobian_Matlab - Calls the function that has been set with SNESSetJacobianMatlab(). 4825 4826 Collective on SNES 4827 4828 Input Parameters: 4829 + snes - the SNES context 4830 . x - input vector 4831 . A, B - the matrices 4832 - ctx - user context 4833 4834 Output Parameter: 4835 . flag - structure of the matrix 4836 4837 Level: developer 4838 4839 .keywords: SNES, nonlinear, compute, function 4840 4841 .seealso: SNESSetFunction(), SNESGetFunction() 4842 @*/ 4843 PetscErrorCode SNESComputeJacobian_Matlab(SNES snes,Vec x,Mat *A,Mat *B,MatStructure *flag, void *ctx) 4844 { 4845 PetscErrorCode ierr; 4846 SNESMatlabContext *sctx = (SNESMatlabContext*)ctx; 4847 int nlhs = 2,nrhs = 6; 4848 mxArray *plhs[2],*prhs[6]; 4849 long long int lx = 0,lA = 0,ls = 0, lB = 0; 4850 4851 PetscFunctionBegin; 4852 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4853 PetscValidHeaderSpecific(x,VEC_CLASSID,2); 4854 4855 /* call Matlab function in ctx with arguments x and y */ 4856 4857 ierr = PetscMemcpy(&ls,&snes,sizeof(snes));CHKERRQ(ierr); 4858 ierr = PetscMemcpy(&lx,&x,sizeof(x));CHKERRQ(ierr); 4859 ierr = PetscMemcpy(&lA,A,sizeof(x));CHKERRQ(ierr); 4860 ierr = PetscMemcpy(&lB,B,sizeof(x));CHKERRQ(ierr); 4861 prhs[0] = mxCreateDoubleScalar((double)ls); 4862 prhs[1] = mxCreateDoubleScalar((double)lx); 4863 prhs[2] = mxCreateDoubleScalar((double)lA); 4864 prhs[3] = mxCreateDoubleScalar((double)lB); 4865 prhs[4] = mxCreateString(sctx->funcname); 4866 prhs[5] = sctx->ctx; 4867 ierr = mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscSNESComputeJacobianInternal");CHKERRQ(ierr); 4868 ierr = mxGetScalar(plhs[0]);CHKERRQ(ierr); 4869 *flag = (MatStructure) mxGetScalar(plhs[1]);CHKERRQ(ierr); 4870 mxDestroyArray(prhs[0]); 4871 mxDestroyArray(prhs[1]); 4872 mxDestroyArray(prhs[2]); 4873 mxDestroyArray(prhs[3]); 4874 mxDestroyArray(prhs[4]); 4875 mxDestroyArray(plhs[0]); 4876 mxDestroyArray(plhs[1]); 4877 PetscFunctionReturn(0); 4878 } 4879 4880 #undef __FUNCT__ 4881 #define __FUNCT__ "SNESSetJacobianMatlab" 4882 /* 4883 SNESSetJacobianMatlab - Sets the Jacobian function evaluation routine and two empty Jacobian matrices 4884 vector for use by the SNES routines in solving systems of nonlinear 4885 equations from MATLAB. Here the function is a string containing the name of a MATLAB function 4886 4887 Logically Collective on SNES 4888 4889 Input Parameters: 4890 + snes - the SNES context 4891 . A,B - Jacobian matrices 4892 . SNESJacobianFunction - function evaluation routine 4893 - ctx - user context 4894 4895 Level: developer 4896 4897 .keywords: SNES, nonlinear, set, function 4898 4899 .seealso: SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESSetFunction(), SNESJacobianFunction 4900 */ 4901 PetscErrorCode SNESSetJacobianMatlab(SNES snes,Mat A,Mat B,const char *SNESJacobianFunction,mxArray *ctx) 4902 { 4903 PetscErrorCode ierr; 4904 SNESMatlabContext *sctx; 4905 4906 PetscFunctionBegin; 4907 /* currently sctx is memory bleed */ 4908 ierr = PetscMalloc(sizeof(SNESMatlabContext),&sctx);CHKERRQ(ierr); 4909 ierr = PetscStrallocpy(SNESJacobianFunction,&sctx->funcname);CHKERRQ(ierr); 4910 /* 4911 This should work, but it doesn't 4912 sctx->ctx = ctx; 4913 mexMakeArrayPersistent(sctx->ctx); 4914 */ 4915 sctx->ctx = mxDuplicateArray(ctx); 4916 ierr = SNESSetJacobian(snes,A,B,SNESComputeJacobian_Matlab,sctx);CHKERRQ(ierr); 4917 PetscFunctionReturn(0); 4918 } 4919 4920 #undef __FUNCT__ 4921 #define __FUNCT__ "SNESMonitor_Matlab" 4922 /* 4923 SNESMonitor_Matlab - Calls the function that has been set with SNESMonitorSetMatlab(). 4924 4925 Collective on SNES 4926 4927 .seealso: SNESSetFunction(), SNESGetFunction() 4928 @*/ 4929 PetscErrorCode SNESMonitor_Matlab(SNES snes,PetscInt it, PetscReal fnorm, void *ctx) 4930 { 4931 PetscErrorCode ierr; 4932 SNESMatlabContext *sctx = (SNESMatlabContext*)ctx; 4933 int nlhs = 1,nrhs = 6; 4934 mxArray *plhs[1],*prhs[6]; 4935 long long int lx = 0,ls = 0; 4936 Vec x = snes->vec_sol; 4937 4938 PetscFunctionBegin; 4939 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4940 4941 ierr = PetscMemcpy(&ls,&snes,sizeof(snes));CHKERRQ(ierr); 4942 ierr = PetscMemcpy(&lx,&x,sizeof(x));CHKERRQ(ierr); 4943 prhs[0] = mxCreateDoubleScalar((double)ls); 4944 prhs[1] = mxCreateDoubleScalar((double)it); 4945 prhs[2] = mxCreateDoubleScalar((double)fnorm); 4946 prhs[3] = mxCreateDoubleScalar((double)lx); 4947 prhs[4] = mxCreateString(sctx->funcname); 4948 prhs[5] = sctx->ctx; 4949 ierr = mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscSNESMonitorInternal");CHKERRQ(ierr); 4950 ierr = mxGetScalar(plhs[0]);CHKERRQ(ierr); 4951 mxDestroyArray(prhs[0]); 4952 mxDestroyArray(prhs[1]); 4953 mxDestroyArray(prhs[2]); 4954 mxDestroyArray(prhs[3]); 4955 mxDestroyArray(prhs[4]); 4956 mxDestroyArray(plhs[0]); 4957 PetscFunctionReturn(0); 4958 } 4959 4960 #undef __FUNCT__ 4961 #define __FUNCT__ "SNESMonitorSetMatlab" 4962 /* 4963 SNESMonitorSetMatlab - Sets the monitor function from MATLAB 4964 4965 Level: developer 4966 4967 .keywords: SNES, nonlinear, set, function 4968 4969 .seealso: SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESSetFunction() 4970 */ 4971 PetscErrorCode SNESMonitorSetMatlab(SNES snes,const char *SNESMonitorFunction,mxArray *ctx) 4972 { 4973 PetscErrorCode ierr; 4974 SNESMatlabContext *sctx; 4975 4976 PetscFunctionBegin; 4977 /* currently sctx is memory bleed */ 4978 ierr = PetscMalloc(sizeof(SNESMatlabContext),&sctx);CHKERRQ(ierr); 4979 ierr = PetscStrallocpy(SNESMonitorFunction,&sctx->funcname);CHKERRQ(ierr); 4980 /* 4981 This should work, but it doesn't 4982 sctx->ctx = ctx; 4983 mexMakeArrayPersistent(sctx->ctx); 4984 */ 4985 sctx->ctx = mxDuplicateArray(ctx); 4986 ierr = SNESMonitorSet(snes,SNESMonitor_Matlab,sctx,NULL);CHKERRQ(ierr); 4987 PetscFunctionReturn(0); 4988 } 4989 4990 #endif 4991