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