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