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