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