1 2 #include <petsc-private/snesimpl.h> /*I "petscsnes.h" I*/ 3 #include <petscdmshell.h> 4 5 PetscBool SNESRegisterAllCalled = PETSC_FALSE; 6 PetscFunctionList SNESList = NULL; 7 8 /* Logging support */ 9 PetscClassId SNES_CLASSID, DMSNES_CLASSID; 10 PetscLogEvent SNES_Solve, SNES_FunctionEval, SNES_JacobianEval, SNES_NGSEval, SNES_NGSFuncEval, SNES_NPCSolve; 11 12 #undef __FUNCT__ 13 #define __FUNCT__ "SNESSetErrorIfNotConverged" 14 /*@ 15 SNESSetErrorIfNotConverged - Causes SNESSolve() to generate an error if the solver has not converged. 16 17 Logically Collective on SNES 18 19 Input Parameters: 20 + snes - iterative context obtained from SNESCreate() 21 - flg - PETSC_TRUE indicates you want the error generated 22 23 Options database keys: 24 . -snes_error_if_not_converged : this takes an optional truth value (0/1/no/yes/true/false) 25 26 Level: intermediate 27 28 Notes: 29 Normally PETSc continues if a linear solver fails to converge, you can call SNESGetConvergedReason() after a SNESSolve() 30 to determine if it has converged. 31 32 .keywords: SNES, set, initial guess, nonzero 33 34 .seealso: SNESGetErrorIfNotConverged(), KSPGetErrorIfNotConverged(), KSPSetErrorIFNotConverged() 35 @*/ 36 PetscErrorCode SNESSetErrorIfNotConverged(SNES snes,PetscBool flg) 37 { 38 PetscFunctionBegin; 39 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 40 PetscValidLogicalCollectiveBool(snes,flg,2); 41 snes->errorifnotconverged = flg; 42 PetscFunctionReturn(0); 43 } 44 45 #undef __FUNCT__ 46 #define __FUNCT__ "SNESGetErrorIfNotConverged" 47 /*@ 48 SNESGetErrorIfNotConverged - Will SNESSolve() generate an error if the solver does not converge? 49 50 Not Collective 51 52 Input Parameter: 53 . snes - iterative context obtained from SNESCreate() 54 55 Output Parameter: 56 . flag - PETSC_TRUE if it will generate an error, else PETSC_FALSE 57 58 Level: intermediate 59 60 .keywords: SNES, set, initial guess, nonzero 61 62 .seealso: SNESSetErrorIfNotConverged(), KSPGetErrorIfNotConverged(), KSPSetErrorIFNotConverged() 63 @*/ 64 PetscErrorCode SNESGetErrorIfNotConverged(SNES snes,PetscBool *flag) 65 { 66 PetscFunctionBegin; 67 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 68 PetscValidPointer(flag,2); 69 *flag = snes->errorifnotconverged; 70 PetscFunctionReturn(0); 71 } 72 73 #undef __FUNCT__ 74 #define __FUNCT__ "SNESSetFunctionDomainError" 75 /*@ 76 SNESSetFunctionDomainError - tells SNES that the input vector to your SNESFunction is not 77 in the functions domain. For example, negative pressure. 78 79 Logically Collective on SNES 80 81 Input Parameters: 82 . snes - the SNES context 83 84 Level: advanced 85 86 .keywords: SNES, view 87 88 .seealso: SNESCreate(), SNESSetFunction(), SNESFunction 89 @*/ 90 PetscErrorCode SNESSetFunctionDomainError(SNES snes) 91 { 92 PetscFunctionBegin; 93 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 94 snes->domainerror = PETSC_TRUE; 95 PetscFunctionReturn(0); 96 } 97 98 #undef __FUNCT__ 99 #define __FUNCT__ "SNESGetFunctionDomainError" 100 /*@ 101 SNESGetFunctionDomainError - Gets the status of the domain error after a call to SNESComputeFunction; 102 103 Logically Collective on SNES 104 105 Input Parameters: 106 . snes - the SNES context 107 108 Output Parameters: 109 . domainerror - Set to PETSC_TRUE if there's a domain error; PETSC_FALSE otherwise. 110 111 Level: advanced 112 113 .keywords: SNES, view 114 115 .seealso: SNESSetFunctionDomainError(), SNESComputeFunction() 116 @*/ 117 PetscErrorCode SNESGetFunctionDomainError(SNES snes, PetscBool *domainerror) 118 { 119 PetscFunctionBegin; 120 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 121 PetscValidPointer(domainerror, 2); 122 *domainerror = snes->domainerror; 123 PetscFunctionReturn(0); 124 } 125 126 #undef __FUNCT__ 127 #define __FUNCT__ "SNESLoad" 128 /*@C 129 SNESLoad - Loads a SNES that has been stored in binary with SNESView(). 130 131 Collective on PetscViewer 132 133 Input Parameters: 134 + newdm - the newly loaded SNES, this needs to have been created with SNESCreate() or 135 some related function before a call to SNESLoad(). 136 - viewer - binary file viewer, obtained from PetscViewerBinaryOpen() 137 138 Level: intermediate 139 140 Notes: 141 The type is determined by the data in the file, any type set into the SNES before this call is ignored. 142 143 Notes for advanced users: 144 Most users should not need to know the details of the binary storage 145 format, since SNESLoad() and TSView() completely hide these details. 146 But for anyone who's interested, the standard binary matrix storage 147 format is 148 .vb 149 has not yet been determined 150 .ve 151 152 .seealso: PetscViewerBinaryOpen(), SNESView(), MatLoad(), VecLoad() 153 @*/ 154 PetscErrorCode SNESLoad(SNES snes, PetscViewer viewer) 155 { 156 PetscErrorCode ierr; 157 PetscBool isbinary; 158 PetscInt classid; 159 char type[256]; 160 KSP ksp; 161 DM dm; 162 DMSNES dmsnes; 163 164 PetscFunctionBegin; 165 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 166 PetscValidHeaderSpecific(viewer,PETSC_VIEWER_CLASSID,2); 167 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr); 168 if (!isbinary) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Invalid viewer; open viewer with PetscViewerBinaryOpen()"); 169 170 ierr = PetscViewerBinaryRead(viewer,&classid,1,PETSC_INT);CHKERRQ(ierr); 171 if (classid != SNES_FILE_CLASSID) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONG,"Not SNES next in file"); 172 ierr = PetscViewerBinaryRead(viewer,type,256,PETSC_CHAR);CHKERRQ(ierr); 173 ierr = SNESSetType(snes, type);CHKERRQ(ierr); 174 if (snes->ops->load) { 175 ierr = (*snes->ops->load)(snes,viewer);CHKERRQ(ierr); 176 } 177 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 178 ierr = DMGetDMSNES(dm,&dmsnes);CHKERRQ(ierr); 179 ierr = DMSNESLoad(dmsnes,viewer);CHKERRQ(ierr); 180 ierr = SNESGetKSP(snes,&ksp);CHKERRQ(ierr); 181 ierr = KSPLoad(ksp,viewer);CHKERRQ(ierr); 182 PetscFunctionReturn(0); 183 } 184 185 #include <petscdraw.h> 186 #if defined(PETSC_HAVE_SAWS) 187 #include <petscviewersaws.h> 188 #endif 189 #undef __FUNCT__ 190 #define __FUNCT__ "SNESView" 191 /*@C 192 SNESView - Prints the SNES data structure. 193 194 Collective on SNES 195 196 Input Parameters: 197 + SNES - the SNES context 198 - viewer - visualization context 199 200 Options Database Key: 201 . -snes_view - Calls SNESView() at end of SNESSolve() 202 203 Notes: 204 The available visualization contexts include 205 + PETSC_VIEWER_STDOUT_SELF - standard output (default) 206 - PETSC_VIEWER_STDOUT_WORLD - synchronized standard 207 output where only the first processor opens 208 the file. All other processors send their 209 data to the first processor to print. 210 211 The user can open an alternative visualization context with 212 PetscViewerASCIIOpen() - output to a specified file. 213 214 Level: beginner 215 216 .keywords: SNES, view 217 218 .seealso: PetscViewerASCIIOpen() 219 @*/ 220 PetscErrorCode SNESView(SNES snes,PetscViewer viewer) 221 { 222 SNESKSPEW *kctx; 223 PetscErrorCode ierr; 224 KSP ksp; 225 SNESLineSearch linesearch; 226 PetscBool iascii,isstring,isbinary,isdraw; 227 DMSNES dmsnes; 228 #if defined(PETSC_HAVE_SAWS) 229 PetscBool issaws; 230 #endif 231 232 PetscFunctionBegin; 233 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 234 if (!viewer) { 235 ierr = PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)snes),&viewer);CHKERRQ(ierr); 236 } 237 PetscValidHeaderSpecific(viewer,PETSC_VIEWER_CLASSID,2); 238 PetscCheckSameComm(snes,1,viewer,2); 239 240 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 241 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSTRING,&isstring);CHKERRQ(ierr); 242 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr); 243 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);CHKERRQ(ierr); 244 #if defined(PETSC_HAVE_SAWS) 245 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSAWS,&issaws);CHKERRQ(ierr); 246 #endif 247 if (iascii) { 248 SNESNormSchedule normschedule; 249 250 ierr = PetscObjectPrintClassNamePrefixType((PetscObject)snes,viewer);CHKERRQ(ierr); 251 if (!snes->setupcalled) { 252 ierr = PetscViewerASCIIPrintf(viewer," SNES has not been set up so information may be incomplete\n");CHKERRQ(ierr); 253 } 254 if (snes->ops->view) { 255 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 256 ierr = (*snes->ops->view)(snes,viewer);CHKERRQ(ierr); 257 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 258 } 259 ierr = PetscViewerASCIIPrintf(viewer," maximum iterations=%D, maximum function evaluations=%D\n",snes->max_its,snes->max_funcs);CHKERRQ(ierr); 260 ierr = PetscViewerASCIIPrintf(viewer," tolerances: relative=%g, absolute=%g, solution=%g\n",(double)snes->rtol,(double)snes->abstol,(double)snes->stol);CHKERRQ(ierr); 261 ierr = PetscViewerASCIIPrintf(viewer," total number of linear solver iterations=%D\n",snes->linear_its);CHKERRQ(ierr); 262 ierr = PetscViewerASCIIPrintf(viewer," total number of function evaluations=%D\n",snes->nfuncs);CHKERRQ(ierr); 263 ierr = SNESGetNormSchedule(snes, &normschedule);CHKERRQ(ierr); 264 if (normschedule > 0) {ierr = PetscViewerASCIIPrintf(viewer," norm schedule %s\n",SNESNormSchedules[normschedule]);CHKERRQ(ierr);} 265 if (snes->gridsequence) { 266 ierr = PetscViewerASCIIPrintf(viewer," total number of grid sequence refinements=%D\n",snes->gridsequence);CHKERRQ(ierr); 267 } 268 if (snes->ksp_ewconv) { 269 kctx = (SNESKSPEW*)snes->kspconvctx; 270 if (kctx) { 271 ierr = PetscViewerASCIIPrintf(viewer," Eisenstat-Walker computation of KSP relative tolerance (version %D)\n",kctx->version);CHKERRQ(ierr); 272 ierr = PetscViewerASCIIPrintf(viewer," rtol_0=%g, rtol_max=%g, threshold=%g\n",(double)kctx->rtol_0,(double)kctx->rtol_max,(double)kctx->threshold);CHKERRQ(ierr); 273 ierr = PetscViewerASCIIPrintf(viewer," gamma=%g, alpha=%g, alpha2=%g\n",(double)kctx->gamma,(double)kctx->alpha,(double)kctx->alpha2);CHKERRQ(ierr); 274 } 275 } 276 if (snes->lagpreconditioner == -1) { 277 ierr = PetscViewerASCIIPrintf(viewer," Preconditioned is never rebuilt\n");CHKERRQ(ierr); 278 } else if (snes->lagpreconditioner > 1) { 279 ierr = PetscViewerASCIIPrintf(viewer," Preconditioned is rebuilt every %D new Jacobians\n",snes->lagpreconditioner);CHKERRQ(ierr); 280 } 281 if (snes->lagjacobian == -1) { 282 ierr = PetscViewerASCIIPrintf(viewer," Jacobian is never rebuilt\n");CHKERRQ(ierr); 283 } else if (snes->lagjacobian > 1) { 284 ierr = PetscViewerASCIIPrintf(viewer," Jacobian is rebuilt every %D SNES iterations\n",snes->lagjacobian);CHKERRQ(ierr); 285 } 286 } else if (isstring) { 287 const char *type; 288 ierr = SNESGetType(snes,&type);CHKERRQ(ierr); 289 ierr = PetscViewerStringSPrintf(viewer," %-3.3s",type);CHKERRQ(ierr); 290 } else if (isbinary) { 291 PetscInt classid = SNES_FILE_CLASSID; 292 MPI_Comm comm; 293 PetscMPIInt rank; 294 char type[256]; 295 296 ierr = PetscObjectGetComm((PetscObject)snes,&comm);CHKERRQ(ierr); 297 ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr); 298 if (!rank) { 299 ierr = PetscViewerBinaryWrite(viewer,&classid,1,PETSC_INT,PETSC_FALSE);CHKERRQ(ierr); 300 ierr = PetscStrncpy(type,((PetscObject)snes)->type_name,sizeof(type));CHKERRQ(ierr); 301 ierr = PetscViewerBinaryWrite(viewer,type,sizeof(type),PETSC_CHAR,PETSC_FALSE);CHKERRQ(ierr); 302 } 303 if (snes->ops->view) { 304 ierr = (*snes->ops->view)(snes,viewer);CHKERRQ(ierr); 305 } 306 } else if (isdraw) { 307 PetscDraw draw; 308 char str[36]; 309 PetscReal x,y,bottom,h; 310 311 ierr = PetscViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr); 312 ierr = PetscDrawGetCurrentPoint(draw,&x,&y);CHKERRQ(ierr); 313 ierr = PetscStrcpy(str,"SNES: ");CHKERRQ(ierr); 314 ierr = PetscStrcat(str,((PetscObject)snes)->type_name);CHKERRQ(ierr); 315 ierr = PetscDrawBoxedString(draw,x,y,PETSC_DRAW_BLUE,PETSC_DRAW_BLACK,str,NULL,&h);CHKERRQ(ierr); 316 bottom = y - h; 317 ierr = PetscDrawPushCurrentPoint(draw,x,bottom);CHKERRQ(ierr); 318 if (snes->ops->view) { 319 ierr = (*snes->ops->view)(snes,viewer);CHKERRQ(ierr); 320 } 321 #if defined(PETSC_HAVE_SAWS) 322 } else if (issaws) { 323 PetscMPIInt rank; 324 const char *name; 325 326 ierr = PetscObjectGetName((PetscObject)snes,&name);CHKERRQ(ierr); 327 ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr); 328 if (!((PetscObject)snes)->amsmem && !rank) { 329 char dir[1024]; 330 331 ierr = PetscObjectViewSAWs((PetscObject)snes,viewer);CHKERRQ(ierr); 332 ierr = PetscSNPrintf(dir,1024,"/PETSc/Objects/%s/its",name);CHKERRQ(ierr); 333 PetscStackCallSAWs(SAWs_Register,(dir,&snes->iter,1,SAWs_READ,SAWs_INT)); 334 if (!snes->conv_hist) { 335 ierr = SNESSetConvergenceHistory(snes,NULL,NULL,PETSC_DECIDE,PETSC_TRUE);CHKERRQ(ierr); 336 } 337 ierr = PetscSNPrintf(dir,1024,"/PETSc/Objects/%s/conv_hist",name);CHKERRQ(ierr); 338 PetscStackCallSAWs(SAWs_Register,(dir,snes->conv_hist,10,SAWs_READ,SAWs_DOUBLE)); 339 } 340 #endif 341 } 342 if (snes->linesearch) { 343 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 344 ierr = SNESGetLineSearch(snes, &linesearch);CHKERRQ(ierr); 345 ierr = SNESLineSearchView(linesearch, viewer);CHKERRQ(ierr); 346 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 347 } 348 if (snes->pc && snes->usespc) { 349 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 350 ierr = SNESView(snes->pc, viewer);CHKERRQ(ierr); 351 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 352 } 353 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 354 ierr = DMGetDMSNES(snes->dm,&dmsnes);CHKERRQ(ierr); 355 ierr = DMSNESView(dmsnes, viewer);CHKERRQ(ierr); 356 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 357 if (snes->usesksp) { 358 ierr = SNESGetKSP(snes,&ksp);CHKERRQ(ierr); 359 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 360 ierr = KSPView(ksp,viewer);CHKERRQ(ierr); 361 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 362 } 363 if (isdraw) { 364 PetscDraw draw; 365 ierr = PetscViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr); 366 ierr = PetscDrawPopCurrentPoint(draw);CHKERRQ(ierr); 367 } 368 PetscFunctionReturn(0); 369 } 370 371 /* 372 We retain a list of functions that also take SNES command 373 line options. These are called at the end SNESSetFromOptions() 374 */ 375 #define MAXSETFROMOPTIONS 5 376 static PetscInt numberofsetfromoptions; 377 static PetscErrorCode (*othersetfromoptions[MAXSETFROMOPTIONS])(SNES); 378 379 #undef __FUNCT__ 380 #define __FUNCT__ "SNESAddOptionsChecker" 381 /*@C 382 SNESAddOptionsChecker - Adds an additional function to check for SNES options. 383 384 Not Collective 385 386 Input Parameter: 387 . snescheck - function that checks for options 388 389 Level: developer 390 391 .seealso: SNESSetFromOptions() 392 @*/ 393 PetscErrorCode SNESAddOptionsChecker(PetscErrorCode (*snescheck)(SNES)) 394 { 395 PetscFunctionBegin; 396 if (numberofsetfromoptions >= MAXSETFROMOPTIONS) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE, "Too many options checkers, only %D allowed", MAXSETFROMOPTIONS); 397 othersetfromoptions[numberofsetfromoptions++] = snescheck; 398 PetscFunctionReturn(0); 399 } 400 401 extern PetscErrorCode SNESDefaultMatrixFreeCreate2(SNES,Vec,Mat*); 402 403 #undef __FUNCT__ 404 #define __FUNCT__ "SNESSetUpMatrixFree_Private" 405 static PetscErrorCode SNESSetUpMatrixFree_Private(SNES snes, PetscBool hasOperator, PetscInt version) 406 { 407 Mat J; 408 KSP ksp; 409 PC pc; 410 PetscBool match; 411 PetscErrorCode ierr; 412 413 PetscFunctionBegin; 414 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 415 416 if (!snes->vec_func && (snes->jacobian || snes->jacobian_pre)) { 417 Mat A = snes->jacobian, B = snes->jacobian_pre; 418 ierr = MatCreateVecs(A ? A : B, NULL,&snes->vec_func);CHKERRQ(ierr); 419 } 420 421 if (version == 1) { 422 ierr = MatCreateSNESMF(snes,&J);CHKERRQ(ierr); 423 ierr = MatMFFDSetOptionsPrefix(J,((PetscObject)snes)->prefix);CHKERRQ(ierr); 424 ierr = MatSetFromOptions(J);CHKERRQ(ierr); 425 } else if (version == 2) { 426 if (!snes->vec_func) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"SNESSetFunction() must be called first"); 427 #if !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_REAL_SINGLE) && !defined(PETSC_USE_REAL___FLOAT128) 428 ierr = SNESDefaultMatrixFreeCreate2(snes,snes->vec_func,&J);CHKERRQ(ierr); 429 #else 430 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP, "matrix-free operator rutines (version 2)"); 431 #endif 432 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE, "matrix-free operator rutines, only version 1 and 2"); 433 434 ierr = PetscInfo1(snes,"Setting default matrix-free operator routines (version %D)\n", version);CHKERRQ(ierr); 435 if (hasOperator) { 436 437 /* This version replaces the user provided Jacobian matrix with a 438 matrix-free version but still employs the user-provided preconditioner matrix. */ 439 ierr = SNESSetJacobian(snes,J,0,0,0);CHKERRQ(ierr); 440 } else { 441 /* This version replaces both the user-provided Jacobian and the user- 442 provided preconditioner Jacobian with the default matrix free version. */ 443 if ((snes->pcside == PC_LEFT) && snes->pc) { 444 if (!snes->jacobian){ierr = SNESSetJacobian(snes,J,0,0,0);CHKERRQ(ierr);} 445 } else { 446 ierr = SNESSetJacobian(snes,J,J,MatMFFDComputeJacobian,0);CHKERRQ(ierr); 447 } 448 /* Force no preconditioner */ 449 ierr = SNESGetKSP(snes,&ksp);CHKERRQ(ierr); 450 ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr); 451 ierr = PetscObjectTypeCompare((PetscObject)pc,PCSHELL,&match);CHKERRQ(ierr); 452 if (!match) { 453 ierr = PetscInfo(snes,"Setting default matrix-free preconditioner routines\nThat is no preconditioner is being used\n");CHKERRQ(ierr); 454 ierr = PCSetType(pc,PCNONE);CHKERRQ(ierr); 455 } 456 } 457 ierr = MatDestroy(&J);CHKERRQ(ierr); 458 PetscFunctionReturn(0); 459 } 460 461 #undef __FUNCT__ 462 #define __FUNCT__ "DMRestrictHook_SNESVecSol" 463 static PetscErrorCode DMRestrictHook_SNESVecSol(DM dmfine,Mat Restrict,Vec Rscale,Mat Inject,DM dmcoarse,void *ctx) 464 { 465 SNES snes = (SNES)ctx; 466 PetscErrorCode ierr; 467 Vec Xfine,Xfine_named = NULL,Xcoarse; 468 469 PetscFunctionBegin; 470 if (PetscLogPrintInfo) { 471 PetscInt finelevel,coarselevel,fineclevel,coarseclevel; 472 ierr = DMGetRefineLevel(dmfine,&finelevel);CHKERRQ(ierr); 473 ierr = DMGetCoarsenLevel(dmfine,&fineclevel);CHKERRQ(ierr); 474 ierr = DMGetRefineLevel(dmcoarse,&coarselevel);CHKERRQ(ierr); 475 ierr = DMGetCoarsenLevel(dmcoarse,&coarseclevel);CHKERRQ(ierr); 476 ierr = PetscInfo4(dmfine,"Restricting SNES solution vector from level %D-%D to level %D-%D\n",finelevel,fineclevel,coarselevel,coarseclevel);CHKERRQ(ierr); 477 } 478 if (dmfine == snes->dm) Xfine = snes->vec_sol; 479 else { 480 ierr = DMGetNamedGlobalVector(dmfine,"SNESVecSol",&Xfine_named);CHKERRQ(ierr); 481 Xfine = Xfine_named; 482 } 483 ierr = DMGetNamedGlobalVector(dmcoarse,"SNESVecSol",&Xcoarse);CHKERRQ(ierr); 484 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 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)(PetscOptionsObject,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->linesearch) { 896 ierr = SNESGetLineSearch(snes, &snes->linesearch);CHKERRQ(ierr); 897 } 898 ierr = SNESLineSearchSetFromOptions(snes->linesearch);CHKERRQ(ierr); 899 900 if (!snes->ksp) {ierr = SNESGetKSP(snes,&snes->ksp);CHKERRQ(ierr);} 901 ierr = KSPSetOperators(snes->ksp,snes->jacobian,snes->jacobian_pre);CHKERRQ(ierr); 902 ierr = KSPSetFromOptions(snes->ksp);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 PetscErrorCode 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, see SNESJacobianFunction for details on its calling sequence 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(), SNESJacobianFunction 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 ierr = VecLockPush(x);CHKERRQ(ierr); 2038 PetscStackPush("SNES user function"); 2039 ierr = (*sdm->ops->computefunction)(snes,x,y,sdm->functionctx);CHKERRQ(ierr); 2040 PetscStackPop; 2041 ierr = VecLockPop(x);CHKERRQ(ierr); 2042 ierr = PetscLogEventEnd(SNES_FunctionEval,snes,x,y,0);CHKERRQ(ierr); 2043 } else if (snes->vec_rhs) { 2044 ierr = MatMult(snes->jacobian, x, y);CHKERRQ(ierr); 2045 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE, "Must call SNESSetFunction() or SNESSetDM() before SNESComputeFunction(), likely called from SNESSolve()."); 2046 if (snes->vec_rhs) { 2047 ierr = VecAXPY(y,-1.0,snes->vec_rhs);CHKERRQ(ierr); 2048 } 2049 snes->nfuncs++; 2050 ierr = VecValidValues(y,3,PETSC_FALSE);CHKERRQ(ierr); 2051 PetscFunctionReturn(0); 2052 } 2053 2054 #undef __FUNCT__ 2055 #define __FUNCT__ "SNESComputeNGS" 2056 /*@ 2057 SNESComputeNGS - Calls the Gauss-Seidel function that has been set with SNESSetNGS(). 2058 2059 Collective on SNES 2060 2061 Input Parameters: 2062 + snes - the SNES context 2063 . x - input vector 2064 - b - rhs vector 2065 2066 Output Parameter: 2067 . x - new solution vector 2068 2069 Notes: 2070 SNESComputeNGS() is typically used within composed nonlinear solver 2071 implementations, so most users would not generally call this routine 2072 themselves. 2073 2074 Level: developer 2075 2076 .keywords: SNES, nonlinear, compute, function 2077 2078 .seealso: SNESSetNGS(), SNESComputeFunction() 2079 @*/ 2080 PetscErrorCode SNESComputeNGS(SNES snes,Vec b,Vec x) 2081 { 2082 PetscErrorCode ierr; 2083 DM dm; 2084 DMSNES sdm; 2085 2086 PetscFunctionBegin; 2087 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2088 PetscValidHeaderSpecific(x,VEC_CLASSID,2); 2089 if (b) PetscValidHeaderSpecific(b,VEC_CLASSID,3); 2090 PetscCheckSameComm(snes,1,x,2); 2091 if (b) PetscCheckSameComm(snes,1,b,3); 2092 if (b) {ierr = VecValidValues(b,2,PETSC_TRUE);CHKERRQ(ierr);} 2093 ierr = PetscLogEventBegin(SNES_NGSEval,snes,x,b,0);CHKERRQ(ierr); 2094 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 2095 ierr = DMGetDMSNES(dm,&sdm);CHKERRQ(ierr); 2096 if (sdm->ops->computegs) { 2097 if (b) {ierr = VecLockPush(b);CHKERRQ(ierr);} 2098 PetscStackPush("SNES user NGS"); 2099 ierr = (*sdm->ops->computegs)(snes,x,b,sdm->gsctx);CHKERRQ(ierr); 2100 PetscStackPop; 2101 if (b) {ierr = VecLockPop(b);CHKERRQ(ierr);} 2102 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE, "Must call SNESSetNGS() before SNESComputeNGS(), likely called from SNESSolve()."); 2103 ierr = PetscLogEventEnd(SNES_NGSEval,snes,x,b,0);CHKERRQ(ierr); 2104 ierr = VecValidValues(x,3,PETSC_FALSE);CHKERRQ(ierr); 2105 PetscFunctionReturn(0); 2106 } 2107 2108 #undef __FUNCT__ 2109 #define __FUNCT__ "SNESComputeJacobian" 2110 /*@ 2111 SNESComputeJacobian - Computes the Jacobian matrix that has been set with SNESSetJacobian(). 2112 2113 Collective on SNES and Mat 2114 2115 Input Parameters: 2116 + snes - the SNES context 2117 - x - input vector 2118 2119 Output Parameters: 2120 + A - Jacobian matrix 2121 - B - optional preconditioning matrix 2122 2123 Options Database Keys: 2124 + -snes_lag_preconditioner <lag> 2125 . -snes_lag_jacobian <lag> 2126 . -snes_compare_explicit - Compare the computed Jacobian to the finite difference Jacobian and output the differences 2127 . -snes_compare_explicit_draw - Compare the computed Jacobian to the finite difference Jacobian and draw the result 2128 . -snes_compare_explicit_contour - Compare the computed Jacobian to the finite difference Jacobian and draw a contour plot with the result 2129 . -snes_compare_operator - Make the comparison options above use the operator instead of the preconditioning matrix 2130 . -snes_compare_coloring - Compute the finite differece Jacobian using coloring and display norms of difference 2131 . -snes_compare_coloring_display - Compute the finite differece Jacobian using coloring and display verbose differences 2132 . -snes_compare_coloring_threshold - Display only those matrix entries that differ by more than a given threshold 2133 . -snes_compare_coloring_threshold_atol - Absolute tolerance for difference in matrix entries to be displayed by -snes_compare_coloring_threshold 2134 . -snes_compare_coloring_threshold_rtol - Relative tolerance for difference in matrix entries to be displayed by -snes_compare_coloring_threshold 2135 . -snes_compare_coloring_draw - Compute the finite differece Jacobian using coloring and draw differences 2136 - -snes_compare_coloring_draw_contour - Compute the finite differece Jacobian using coloring and show contours of matrices and differences 2137 2138 2139 Notes: 2140 Most users should not need to explicitly call this routine, as it 2141 is used internally within the nonlinear solvers. 2142 2143 Level: developer 2144 2145 .keywords: SNES, compute, Jacobian, matrix 2146 2147 .seealso: SNESSetJacobian(), KSPSetOperators(), MatStructure, SNESSetLagPreconditioner(), SNESSetLagJacobian() 2148 @*/ 2149 PetscErrorCode SNESComputeJacobian(SNES snes,Vec X,Mat A,Mat B) 2150 { 2151 PetscErrorCode ierr; 2152 PetscBool flag; 2153 DM dm; 2154 DMSNES sdm; 2155 KSP ksp; 2156 2157 PetscFunctionBegin; 2158 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2159 PetscValidHeaderSpecific(X,VEC_CLASSID,2); 2160 PetscCheckSameComm(snes,1,X,2); 2161 ierr = VecValidValues(X,2,PETSC_TRUE);CHKERRQ(ierr); 2162 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 2163 ierr = DMGetDMSNES(dm,&sdm);CHKERRQ(ierr); 2164 2165 if (!sdm->ops->computejacobian) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_USER,"Must call SNESSetJacobian(), DMSNESSetJacobian(), DMDASNESSetJacobianLocal(), etc"); 2166 2167 /* make sure that MatAssemblyBegin/End() is called on A matrix if it is matrix free */ 2168 2169 if (snes->lagjacobian == -2) { 2170 snes->lagjacobian = -1; 2171 2172 ierr = PetscInfo(snes,"Recomputing Jacobian/preconditioner because lag is -2 (means compute Jacobian, but then never again) \n");CHKERRQ(ierr); 2173 } else if (snes->lagjacobian == -1) { 2174 ierr = PetscInfo(snes,"Reusing Jacobian/preconditioner because lag is -1\n");CHKERRQ(ierr); 2175 ierr = PetscObjectTypeCompare((PetscObject)A,MATMFFD,&flag);CHKERRQ(ierr); 2176 if (flag) { 2177 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2178 ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2179 } 2180 PetscFunctionReturn(0); 2181 } else if (snes->lagjacobian > 1 && (snes->iter + snes->jac_iter) % snes->lagjacobian) { 2182 ierr = PetscInfo2(snes,"Reusing Jacobian/preconditioner because lag is %D and SNES iteration is %D\n",snes->lagjacobian,snes->iter);CHKERRQ(ierr); 2183 ierr = PetscObjectTypeCompare((PetscObject)A,MATMFFD,&flag);CHKERRQ(ierr); 2184 if (flag) { 2185 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2186 ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2187 } 2188 PetscFunctionReturn(0); 2189 } 2190 if (snes->pc && snes->pcside == PC_LEFT) { 2191 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2192 ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2193 PetscFunctionReturn(0); 2194 } 2195 2196 ierr = PetscLogEventBegin(SNES_JacobianEval,snes,X,A,B);CHKERRQ(ierr); 2197 ierr = VecLockPush(X);CHKERRQ(ierr); 2198 PetscStackPush("SNES user Jacobian function"); 2199 ierr = (*sdm->ops->computejacobian)(snes,X,A,B,sdm->jacobianctx);CHKERRQ(ierr); 2200 PetscStackPop; 2201 ierr = VecLockPop(X);CHKERRQ(ierr); 2202 ierr = PetscLogEventEnd(SNES_JacobianEval,snes,X,A,B);CHKERRQ(ierr); 2203 2204 /* the next line ensures that snes->ksp exists */ 2205 ierr = SNESGetKSP(snes,&ksp);CHKERRQ(ierr); 2206 if (snes->lagpreconditioner == -2) { 2207 ierr = PetscInfo(snes,"Rebuilding preconditioner exactly once since lag is -2\n");CHKERRQ(ierr); 2208 ierr = KSPSetReusePreconditioner(snes->ksp,PETSC_FALSE);CHKERRQ(ierr); 2209 snes->lagpreconditioner = -1; 2210 } else if (snes->lagpreconditioner == -1) { 2211 ierr = PetscInfo(snes,"Reusing preconditioner because lag is -1\n");CHKERRQ(ierr); 2212 ierr = KSPSetReusePreconditioner(snes->ksp,PETSC_TRUE);CHKERRQ(ierr); 2213 } else if (snes->lagpreconditioner > 1 && (snes->iter + snes->pre_iter) % snes->lagpreconditioner) { 2214 ierr = PetscInfo2(snes,"Reusing preconditioner because lag is %D and SNES iteration is %D\n",snes->lagpreconditioner,snes->iter);CHKERRQ(ierr); 2215 ierr = KSPSetReusePreconditioner(snes->ksp,PETSC_TRUE);CHKERRQ(ierr); 2216 } else { 2217 ierr = PetscInfo(snes,"Rebuilding preconditioner\n");CHKERRQ(ierr); 2218 ierr = KSPSetReusePreconditioner(snes->ksp,PETSC_FALSE);CHKERRQ(ierr); 2219 } 2220 2221 /* make sure user returned a correct Jacobian and preconditioner */ 2222 /* PetscValidHeaderSpecific(A,MAT_CLASSID,3); 2223 PetscValidHeaderSpecific(B,MAT_CLASSID,4); */ 2224 { 2225 PetscBool flag = PETSC_FALSE,flag_draw = PETSC_FALSE,flag_contour = PETSC_FALSE,flag_operator = PETSC_FALSE; 2226 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_explicit",&flag,NULL);CHKERRQ(ierr); 2227 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_explicit_draw",&flag_draw,NULL);CHKERRQ(ierr); 2228 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_explicit_draw_contour",&flag_contour,NULL);CHKERRQ(ierr); 2229 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_operator",&flag_operator,NULL);CHKERRQ(ierr); 2230 if (flag || flag_draw || flag_contour) { 2231 Mat Bexp_mine = NULL,Bexp,FDexp; 2232 PetscViewer vdraw,vstdout; 2233 PetscBool flg; 2234 if (flag_operator) { 2235 ierr = MatComputeExplicitOperator(A,&Bexp_mine);CHKERRQ(ierr); 2236 Bexp = Bexp_mine; 2237 } else { 2238 /* See if the preconditioning matrix can be viewed and added directly */ 2239 ierr = PetscObjectTypeCompareAny((PetscObject)B,&flg,MATSEQAIJ,MATMPIAIJ,MATSEQDENSE,MATMPIDENSE,MATSEQBAIJ,MATMPIBAIJ,MATSEQSBAIJ,MATMPIBAIJ,"");CHKERRQ(ierr); 2240 if (flg) Bexp = B; 2241 else { 2242 /* If the "preconditioning" matrix is itself MATSHELL or some other type without direct support */ 2243 ierr = MatComputeExplicitOperator(B,&Bexp_mine);CHKERRQ(ierr); 2244 Bexp = Bexp_mine; 2245 } 2246 } 2247 ierr = MatConvert(Bexp,MATSAME,MAT_INITIAL_MATRIX,&FDexp);CHKERRQ(ierr); 2248 ierr = SNESComputeJacobianDefault(snes,X,FDexp,FDexp,NULL);CHKERRQ(ierr); 2249 ierr = PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)snes),&vstdout);CHKERRQ(ierr); 2250 if (flag_draw || flag_contour) { 2251 ierr = PetscViewerDrawOpen(PetscObjectComm((PetscObject)snes),0,"Explicit Jacobians",PETSC_DECIDE,PETSC_DECIDE,300,300,&vdraw);CHKERRQ(ierr); 2252 if (flag_contour) {ierr = PetscViewerPushFormat(vdraw,PETSC_VIEWER_DRAW_CONTOUR);CHKERRQ(ierr);} 2253 } else vdraw = NULL; 2254 ierr = PetscViewerASCIIPrintf(vstdout,"Explicit %s\n",flag_operator ? "Jacobian" : "preconditioning Jacobian");CHKERRQ(ierr); 2255 if (flag) {ierr = MatView(Bexp,vstdout);CHKERRQ(ierr);} 2256 if (vdraw) {ierr = MatView(Bexp,vdraw);CHKERRQ(ierr);} 2257 ierr = PetscViewerASCIIPrintf(vstdout,"Finite difference Jacobian\n");CHKERRQ(ierr); 2258 if (flag) {ierr = MatView(FDexp,vstdout);CHKERRQ(ierr);} 2259 if (vdraw) {ierr = MatView(FDexp,vdraw);CHKERRQ(ierr);} 2260 ierr = MatAYPX(FDexp,-1.0,Bexp,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 2261 ierr = PetscViewerASCIIPrintf(vstdout,"User-provided matrix minus finite difference Jacobian\n");CHKERRQ(ierr); 2262 if (flag) {ierr = MatView(FDexp,vstdout);CHKERRQ(ierr);} 2263 if (vdraw) { /* Always use contour for the difference */ 2264 ierr = PetscViewerPushFormat(vdraw,PETSC_VIEWER_DRAW_CONTOUR);CHKERRQ(ierr); 2265 ierr = MatView(FDexp,vdraw);CHKERRQ(ierr); 2266 ierr = PetscViewerPopFormat(vdraw);CHKERRQ(ierr); 2267 } 2268 if (flag_contour) {ierr = PetscViewerPopFormat(vdraw);CHKERRQ(ierr);} 2269 ierr = PetscViewerDestroy(&vdraw);CHKERRQ(ierr); 2270 ierr = MatDestroy(&Bexp_mine);CHKERRQ(ierr); 2271 ierr = MatDestroy(&FDexp);CHKERRQ(ierr); 2272 } 2273 } 2274 { 2275 PetscBool flag = PETSC_FALSE,flag_display = PETSC_FALSE,flag_draw = PETSC_FALSE,flag_contour = PETSC_FALSE,flag_threshold = PETSC_FALSE; 2276 PetscReal threshold_atol = PETSC_SQRT_MACHINE_EPSILON,threshold_rtol = 10*PETSC_SQRT_MACHINE_EPSILON; 2277 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_coloring",&flag,NULL);CHKERRQ(ierr); 2278 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_coloring_display",&flag_display,NULL);CHKERRQ(ierr); 2279 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_coloring_draw",&flag_draw,NULL);CHKERRQ(ierr); 2280 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_coloring_draw_contour",&flag_contour,NULL);CHKERRQ(ierr); 2281 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_coloring_threshold",&flag_threshold,NULL);CHKERRQ(ierr); 2282 ierr = PetscOptionsGetReal(((PetscObject)snes)->prefix,"-snes_compare_coloring_threshold_rtol",&threshold_rtol,NULL);CHKERRQ(ierr); 2283 ierr = PetscOptionsGetReal(((PetscObject)snes)->prefix,"-snes_compare_coloring_threshold_atol",&threshold_atol,NULL);CHKERRQ(ierr); 2284 if (flag || flag_display || flag_draw || flag_contour || flag_threshold) { 2285 Mat Bfd; 2286 PetscViewer vdraw,vstdout; 2287 MatColoring coloring; 2288 ISColoring iscoloring; 2289 MatFDColoring matfdcoloring; 2290 PetscErrorCode (*func)(SNES,Vec,Vec,void*); 2291 void *funcctx; 2292 PetscReal norm1,norm2,normmax; 2293 2294 ierr = MatDuplicate(B,MAT_DO_NOT_COPY_VALUES,&Bfd);CHKERRQ(ierr); 2295 ierr = MatColoringCreate(Bfd,&coloring);CHKERRQ(ierr); 2296 ierr = MatColoringSetType(coloring,MATCOLORINGSL);CHKERRQ(ierr); 2297 ierr = MatColoringSetFromOptions(coloring);CHKERRQ(ierr); 2298 ierr = MatColoringApply(coloring,&iscoloring);CHKERRQ(ierr); 2299 ierr = MatColoringDestroy(&coloring);CHKERRQ(ierr); 2300 ierr = MatFDColoringCreate(Bfd,iscoloring,&matfdcoloring);CHKERRQ(ierr); 2301 ierr = MatFDColoringSetFromOptions(matfdcoloring);CHKERRQ(ierr); 2302 ierr = MatFDColoringSetUp(Bfd,iscoloring,matfdcoloring);CHKERRQ(ierr); 2303 ierr = ISColoringDestroy(&iscoloring);CHKERRQ(ierr); 2304 2305 /* This method of getting the function is currently unreliable since it doesn't work for DM local functions. */ 2306 ierr = SNESGetFunction(snes,NULL,&func,&funcctx);CHKERRQ(ierr); 2307 ierr = MatFDColoringSetFunction(matfdcoloring,(PetscErrorCode (*)(void))func,funcctx);CHKERRQ(ierr); 2308 ierr = PetscObjectSetOptionsPrefix((PetscObject)matfdcoloring,((PetscObject)snes)->prefix);CHKERRQ(ierr); 2309 ierr = PetscObjectAppendOptionsPrefix((PetscObject)matfdcoloring,"coloring_");CHKERRQ(ierr); 2310 ierr = MatFDColoringSetFromOptions(matfdcoloring);CHKERRQ(ierr); 2311 ierr = MatFDColoringApply(Bfd,matfdcoloring,X,snes);CHKERRQ(ierr); 2312 ierr = MatFDColoringDestroy(&matfdcoloring);CHKERRQ(ierr); 2313 2314 ierr = PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)snes),&vstdout);CHKERRQ(ierr); 2315 if (flag_draw || flag_contour) { 2316 ierr = PetscViewerDrawOpen(PetscObjectComm((PetscObject)snes),0,"Colored Jacobians",PETSC_DECIDE,PETSC_DECIDE,300,300,&vdraw);CHKERRQ(ierr); 2317 if (flag_contour) {ierr = PetscViewerPushFormat(vdraw,PETSC_VIEWER_DRAW_CONTOUR);CHKERRQ(ierr);} 2318 } else vdraw = NULL; 2319 ierr = PetscViewerASCIIPrintf(vstdout,"Explicit preconditioning Jacobian\n");CHKERRQ(ierr); 2320 if (flag_display) {ierr = MatView(B,vstdout);CHKERRQ(ierr);} 2321 if (vdraw) {ierr = MatView(B,vdraw);CHKERRQ(ierr);} 2322 ierr = PetscViewerASCIIPrintf(vstdout,"Colored Finite difference Jacobian\n");CHKERRQ(ierr); 2323 if (flag_display) {ierr = MatView(Bfd,vstdout);CHKERRQ(ierr);} 2324 if (vdraw) {ierr = MatView(Bfd,vdraw);CHKERRQ(ierr);} 2325 ierr = MatAYPX(Bfd,-1.0,B,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 2326 ierr = MatNorm(Bfd,NORM_1,&norm1);CHKERRQ(ierr); 2327 ierr = MatNorm(Bfd,NORM_FROBENIUS,&norm2);CHKERRQ(ierr); 2328 ierr = MatNorm(Bfd,NORM_MAX,&normmax);CHKERRQ(ierr); 2329 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); 2330 if (flag_display) {ierr = MatView(Bfd,vstdout);CHKERRQ(ierr);} 2331 if (vdraw) { /* Always use contour for the difference */ 2332 ierr = PetscViewerPushFormat(vdraw,PETSC_VIEWER_DRAW_CONTOUR);CHKERRQ(ierr); 2333 ierr = MatView(Bfd,vdraw);CHKERRQ(ierr); 2334 ierr = PetscViewerPopFormat(vdraw);CHKERRQ(ierr); 2335 } 2336 if (flag_contour) {ierr = PetscViewerPopFormat(vdraw);CHKERRQ(ierr);} 2337 2338 if (flag_threshold) { 2339 PetscInt bs,rstart,rend,i; 2340 ierr = MatGetBlockSize(B,&bs);CHKERRQ(ierr); 2341 ierr = MatGetOwnershipRange(B,&rstart,&rend);CHKERRQ(ierr); 2342 for (i=rstart; i<rend; i++) { 2343 const PetscScalar *ba,*ca; 2344 const PetscInt *bj,*cj; 2345 PetscInt bn,cn,j,maxentrycol = -1,maxdiffcol = -1,maxrdiffcol = -1; 2346 PetscReal maxentry = 0,maxdiff = 0,maxrdiff = 0; 2347 ierr = MatGetRow(B,i,&bn,&bj,&ba);CHKERRQ(ierr); 2348 ierr = MatGetRow(Bfd,i,&cn,&cj,&ca);CHKERRQ(ierr); 2349 if (bn != cn) SETERRQ(((PetscObject)A)->comm,PETSC_ERR_PLIB,"Unexpected different nonzero pattern in -snes_compare_coloring_threshold"); 2350 for (j=0; j<bn; j++) { 2351 PetscReal rdiff = PetscAbsScalar(ca[j]) / (threshold_atol + threshold_rtol*PetscAbsScalar(ba[j])); 2352 if (PetscAbsScalar(ba[j]) > PetscAbs(maxentry)) { 2353 maxentrycol = bj[j]; 2354 maxentry = PetscRealPart(ba[j]); 2355 } 2356 if (PetscAbsScalar(ca[j]) > PetscAbs(maxdiff)) { 2357 maxdiffcol = bj[j]; 2358 maxdiff = PetscRealPart(ca[j]); 2359 } 2360 if (rdiff > maxrdiff) { 2361 maxrdiffcol = bj[j]; 2362 maxrdiff = rdiff; 2363 } 2364 } 2365 if (maxrdiff > 1) { 2366 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); 2367 for (j=0; j<bn; j++) { 2368 PetscReal rdiff; 2369 rdiff = PetscAbsScalar(ca[j]) / (threshold_atol + threshold_rtol*PetscAbsScalar(ba[j])); 2370 if (rdiff > 1) { 2371 ierr = PetscViewerASCIIPrintf(vstdout," (%D,%g:%g)",bj[j],(double)PetscRealPart(ba[j]),(double)PetscRealPart(ca[j]));CHKERRQ(ierr); 2372 } 2373 } 2374 ierr = PetscViewerASCIIPrintf(vstdout,"\n",i,maxentry,maxdiff,maxrdiff);CHKERRQ(ierr); 2375 } 2376 ierr = MatRestoreRow(B,i,&bn,&bj,&ba);CHKERRQ(ierr); 2377 ierr = MatRestoreRow(Bfd,i,&cn,&cj,&ca);CHKERRQ(ierr); 2378 } 2379 } 2380 ierr = PetscViewerDestroy(&vdraw);CHKERRQ(ierr); 2381 ierr = MatDestroy(&Bfd);CHKERRQ(ierr); 2382 } 2383 } 2384 PetscFunctionReturn(0); 2385 } 2386 2387 /*MC 2388 SNESJacobianFunction - Function used to convey the nonlinear Jacobian of the function to be solved by SNES 2389 2390 Synopsis: 2391 #include "petscsnes.h" 2392 PetscErrorCode SNESJacobianFunction(SNES snes,Vec x,Mat Amat,Mat Pmat,void *ctx); 2393 2394 + x - input vector 2395 . Amat - the matrix that defines the (approximate) Jacobian 2396 . Pmat - the matrix to be used in constructing the preconditioner, usually the same as Amat. 2397 - ctx - [optional] user-defined Jacobian context 2398 2399 Level: intermediate 2400 2401 .seealso: SNESSetFunction(), SNESGetFunction(), SNESSetJacobian(), SNESGetJacobian() 2402 M*/ 2403 2404 #undef __FUNCT__ 2405 #define __FUNCT__ "SNESSetJacobian" 2406 /*@C 2407 SNESSetJacobian - Sets the function to compute Jacobian as well as the 2408 location to store the matrix. 2409 2410 Logically Collective on SNES and Mat 2411 2412 Input Parameters: 2413 + snes - the SNES context 2414 . Amat - the matrix that defines the (approximate) Jacobian 2415 . Pmat - the matrix to be used in constructing the preconditioner, usually the same as Amat. 2416 . J - Jacobian evaluation routine (if NULL then SNES retains any previously set value), see SNESJacobianFunction for details 2417 - ctx - [optional] user-defined context for private data for the 2418 Jacobian evaluation routine (may be NULL) (if NULL then SNES retains any previously set value) 2419 2420 Notes: 2421 If the Amat matrix and Pmat matrix are different you must call MatAssemblyBegin/End() on 2422 each matrix. 2423 2424 If using SNESComputeJacobianDefaultColor() to assemble a Jacobian, the ctx argument 2425 must be a MatFDColoring. 2426 2427 Other defect-correction schemes can be used by computing a different matrix in place of the Jacobian. One common 2428 example is to use the "Picard linearization" which only differentiates through the highest order parts of each term. 2429 2430 Level: beginner 2431 2432 .keywords: SNES, nonlinear, set, Jacobian, matrix 2433 2434 .seealso: KSPSetOperators(), SNESSetFunction(), MatMFFDComputeJacobian(), SNESComputeJacobianDefaultColor(), MatStructure, J, 2435 SNESSetPicard(), SNESJacobianFunction 2436 @*/ 2437 PetscErrorCode SNESSetJacobian(SNES snes,Mat Amat,Mat Pmat,PetscErrorCode (*J)(SNES,Vec,Mat,Mat,void*),void *ctx) 2438 { 2439 PetscErrorCode ierr; 2440 DM dm; 2441 2442 PetscFunctionBegin; 2443 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2444 if (Amat) PetscValidHeaderSpecific(Amat,MAT_CLASSID,2); 2445 if (Pmat) PetscValidHeaderSpecific(Pmat,MAT_CLASSID,3); 2446 if (Amat) PetscCheckSameComm(snes,1,Amat,2); 2447 if (Pmat) PetscCheckSameComm(snes,1,Pmat,3); 2448 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 2449 ierr = DMSNESSetJacobian(dm,J,ctx);CHKERRQ(ierr); 2450 if (Amat) { 2451 ierr = PetscObjectReference((PetscObject)Amat);CHKERRQ(ierr); 2452 ierr = MatDestroy(&snes->jacobian);CHKERRQ(ierr); 2453 2454 snes->jacobian = Amat; 2455 } 2456 if (Pmat) { 2457 ierr = PetscObjectReference((PetscObject)Pmat);CHKERRQ(ierr); 2458 ierr = MatDestroy(&snes->jacobian_pre);CHKERRQ(ierr); 2459 2460 snes->jacobian_pre = Pmat; 2461 } 2462 PetscFunctionReturn(0); 2463 } 2464 2465 #undef __FUNCT__ 2466 #define __FUNCT__ "SNESGetJacobian" 2467 /*@C 2468 SNESGetJacobian - Returns the Jacobian matrix and optionally the user 2469 provided context for evaluating the Jacobian. 2470 2471 Not Collective, but Mat object will be parallel if SNES object is 2472 2473 Input Parameter: 2474 . snes - the nonlinear solver context 2475 2476 Output Parameters: 2477 + Amat - location to stash (approximate) Jacobian matrix (or NULL) 2478 . Pmat - location to stash matrix used to compute the preconditioner (or NULL) 2479 . J - location to put Jacobian function (or NULL), see SNESJacobianFunction for details on its calling sequence 2480 - ctx - location to stash Jacobian ctx (or NULL) 2481 2482 Level: advanced 2483 2484 .seealso: SNESSetJacobian(), SNESComputeJacobian(), SNESJacobianFunction, SNESGetFunction() 2485 @*/ 2486 PetscErrorCode SNESGetJacobian(SNES snes,Mat *Amat,Mat *Pmat,PetscErrorCode (**J)(SNES,Vec,Mat,Mat,void*),void **ctx) 2487 { 2488 PetscErrorCode ierr; 2489 DM dm; 2490 DMSNES sdm; 2491 2492 PetscFunctionBegin; 2493 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2494 if (Amat) *Amat = snes->jacobian; 2495 if (Pmat) *Pmat = snes->jacobian_pre; 2496 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 2497 ierr = DMGetDMSNES(dm,&sdm);CHKERRQ(ierr); 2498 if (J) *J = sdm->ops->computejacobian; 2499 if (ctx) *ctx = sdm->jacobianctx; 2500 PetscFunctionReturn(0); 2501 } 2502 2503 #undef __FUNCT__ 2504 #define __FUNCT__ "SNESSetUp" 2505 /*@ 2506 SNESSetUp - Sets up the internal data structures for the later use 2507 of a nonlinear solver. 2508 2509 Collective on SNES 2510 2511 Input Parameters: 2512 . snes - the SNES context 2513 2514 Notes: 2515 For basic use of the SNES solvers the user need not explicitly call 2516 SNESSetUp(), since these actions will automatically occur during 2517 the call to SNESSolve(). However, if one wishes to control this 2518 phase separately, SNESSetUp() should be called after SNESCreate() 2519 and optional routines of the form SNESSetXXX(), but before SNESSolve(). 2520 2521 Level: advanced 2522 2523 .keywords: SNES, nonlinear, setup 2524 2525 .seealso: SNESCreate(), SNESSolve(), SNESDestroy() 2526 @*/ 2527 PetscErrorCode SNESSetUp(SNES snes) 2528 { 2529 PetscErrorCode ierr; 2530 DM dm; 2531 DMSNES sdm; 2532 SNESLineSearch linesearch, pclinesearch; 2533 void *lsprectx,*lspostctx; 2534 PetscErrorCode (*precheck)(SNESLineSearch,Vec,Vec,PetscBool*,void*); 2535 PetscErrorCode (*postcheck)(SNESLineSearch,Vec,Vec,Vec,PetscBool*,PetscBool*,void*); 2536 PetscErrorCode (*func)(SNES,Vec,Vec,void*); 2537 Vec f,fpc; 2538 void *funcctx; 2539 PetscErrorCode (*jac)(SNES,Vec,Mat,Mat,void*); 2540 void *jacctx,*appctx; 2541 Mat j,jpre; 2542 2543 PetscFunctionBegin; 2544 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2545 if (snes->setupcalled) PetscFunctionReturn(0); 2546 2547 if (!((PetscObject)snes)->type_name) { 2548 ierr = SNESSetType(snes,SNESNEWTONLS);CHKERRQ(ierr); 2549 } 2550 2551 ierr = SNESGetFunction(snes,&snes->vec_func,NULL,NULL);CHKERRQ(ierr); 2552 2553 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 2554 ierr = DMGetDMSNES(dm,&sdm);CHKERRQ(ierr); 2555 if (!sdm->ops->computefunction) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"Function never provided to SNES object"); 2556 if (!sdm->ops->computejacobian) { 2557 ierr = DMSNESSetJacobian(dm,SNESComputeJacobianDefaultColor,NULL);CHKERRQ(ierr); 2558 } 2559 if (!snes->vec_func) { 2560 ierr = DMCreateGlobalVector(dm,&snes->vec_func);CHKERRQ(ierr); 2561 } 2562 2563 if (!snes->ksp) { 2564 ierr = SNESGetKSP(snes, &snes->ksp);CHKERRQ(ierr); 2565 } 2566 2567 if (!snes->linesearch) { 2568 ierr = SNESGetLineSearch(snes, &snes->linesearch);CHKERRQ(ierr); 2569 } 2570 ierr = SNESLineSearchSetFunction(snes->linesearch,SNESComputeFunction);CHKERRQ(ierr); 2571 2572 if (snes->pc && (snes->pcside == PC_LEFT)) { 2573 snes->mf = PETSC_TRUE; 2574 snes->mf_operator = PETSC_FALSE; 2575 } 2576 2577 if (snes->pc) { 2578 /* copy the DM over */ 2579 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 2580 ierr = SNESSetDM(snes->pc,dm);CHKERRQ(ierr); 2581 2582 ierr = SNESGetFunction(snes,&f,&func,&funcctx);CHKERRQ(ierr); 2583 ierr = VecDuplicate(f,&fpc);CHKERRQ(ierr); 2584 ierr = SNESSetFunction(snes->pc,fpc,func,funcctx);CHKERRQ(ierr); 2585 ierr = SNESGetJacobian(snes,&j,&jpre,&jac,&jacctx);CHKERRQ(ierr); 2586 ierr = SNESSetJacobian(snes->pc,j,jpre,jac,jacctx);CHKERRQ(ierr); 2587 ierr = SNESGetApplicationContext(snes,&appctx);CHKERRQ(ierr); 2588 ierr = SNESSetApplicationContext(snes->pc,appctx);CHKERRQ(ierr); 2589 ierr = VecDestroy(&fpc);CHKERRQ(ierr); 2590 2591 /* copy the function pointers over */ 2592 ierr = PetscObjectCopyFortranFunctionPointers((PetscObject)snes,(PetscObject)snes->pc);CHKERRQ(ierr); 2593 2594 /* default to 1 iteration */ 2595 ierr = SNESSetTolerances(snes->pc,0.0,0.0,0.0,1,snes->pc->max_funcs);CHKERRQ(ierr); 2596 if (snes->pcside==PC_RIGHT) { 2597 ierr = SNESSetNormSchedule(snes->pc,SNES_NORM_FINAL_ONLY);CHKERRQ(ierr); 2598 } else { 2599 ierr = SNESSetNormSchedule(snes->pc,SNES_NORM_NONE);CHKERRQ(ierr); 2600 } 2601 ierr = SNESSetFromOptions(snes->pc);CHKERRQ(ierr); 2602 2603 /* copy the line search context over */ 2604 ierr = SNESGetLineSearch(snes,&linesearch);CHKERRQ(ierr); 2605 ierr = SNESGetLineSearch(snes->pc,&pclinesearch);CHKERRQ(ierr); 2606 ierr = SNESLineSearchGetPreCheck(linesearch,&precheck,&lsprectx);CHKERRQ(ierr); 2607 ierr = SNESLineSearchGetPostCheck(linesearch,&postcheck,&lspostctx);CHKERRQ(ierr); 2608 ierr = SNESLineSearchSetPreCheck(pclinesearch,precheck,lsprectx);CHKERRQ(ierr); 2609 ierr = SNESLineSearchSetPostCheck(pclinesearch,postcheck,lspostctx);CHKERRQ(ierr); 2610 ierr = PetscObjectCopyFortranFunctionPointers((PetscObject)linesearch, (PetscObject)pclinesearch);CHKERRQ(ierr); 2611 } 2612 if (snes->mf) { 2613 ierr = SNESSetUpMatrixFree_Private(snes, snes->mf_operator, snes->mf_version);CHKERRQ(ierr); 2614 } 2615 if (snes->ops->usercompute && !snes->user) { 2616 ierr = (*snes->ops->usercompute)(snes,(void**)&snes->user);CHKERRQ(ierr); 2617 } 2618 2619 snes->jac_iter = 0; 2620 snes->pre_iter = 0; 2621 2622 if (snes->ops->setup) { 2623 ierr = (*snes->ops->setup)(snes);CHKERRQ(ierr); 2624 } 2625 2626 if (snes->pc && (snes->pcside == PC_LEFT)) { 2627 if (snes->functype == SNES_FUNCTION_PRECONDITIONED) { 2628 ierr = SNESGetLineSearch(snes,&linesearch);CHKERRQ(ierr); 2629 ierr = SNESLineSearchSetFunction(linesearch,SNESComputeFunctionDefaultNPC);CHKERRQ(ierr); 2630 } 2631 } 2632 2633 snes->setupcalled = PETSC_TRUE; 2634 PetscFunctionReturn(0); 2635 } 2636 2637 #undef __FUNCT__ 2638 #define __FUNCT__ "SNESReset" 2639 /*@ 2640 SNESReset - Resets a SNES context to the snessetupcalled = 0 state and removes any allocated Vecs and Mats 2641 2642 Collective on SNES 2643 2644 Input Parameter: 2645 . snes - iterative context obtained from SNESCreate() 2646 2647 Level: intermediate 2648 2649 Notes: Also calls the application context destroy routine set with SNESSetComputeApplicationContext() 2650 2651 .keywords: SNES, destroy 2652 2653 .seealso: SNESCreate(), SNESSetUp(), SNESSolve() 2654 @*/ 2655 PetscErrorCode SNESReset(SNES snes) 2656 { 2657 PetscErrorCode ierr; 2658 2659 PetscFunctionBegin; 2660 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2661 if (snes->ops->userdestroy && snes->user) { 2662 ierr = (*snes->ops->userdestroy)((void**)&snes->user);CHKERRQ(ierr); 2663 snes->user = NULL; 2664 } 2665 if (snes->pc) { 2666 ierr = SNESReset(snes->pc);CHKERRQ(ierr); 2667 } 2668 2669 if (snes->ops->reset) { 2670 ierr = (*snes->ops->reset)(snes);CHKERRQ(ierr); 2671 } 2672 if (snes->ksp) { 2673 ierr = KSPReset(snes->ksp);CHKERRQ(ierr); 2674 } 2675 2676 if (snes->linesearch) { 2677 ierr = SNESLineSearchReset(snes->linesearch);CHKERRQ(ierr); 2678 } 2679 2680 ierr = VecDestroy(&snes->vec_rhs);CHKERRQ(ierr); 2681 ierr = VecDestroy(&snes->vec_sol);CHKERRQ(ierr); 2682 ierr = VecDestroy(&snes->vec_sol_update);CHKERRQ(ierr); 2683 ierr = VecDestroy(&snes->vec_func);CHKERRQ(ierr); 2684 ierr = MatDestroy(&snes->jacobian);CHKERRQ(ierr); 2685 ierr = MatDestroy(&snes->jacobian_pre);CHKERRQ(ierr); 2686 ierr = VecDestroyVecs(snes->nwork,&snes->work);CHKERRQ(ierr); 2687 ierr = VecDestroyVecs(snes->nvwork,&snes->vwork);CHKERRQ(ierr); 2688 2689 snes->nwork = snes->nvwork = 0; 2690 snes->setupcalled = PETSC_FALSE; 2691 PetscFunctionReturn(0); 2692 } 2693 2694 #undef __FUNCT__ 2695 #define __FUNCT__ "SNESDestroy" 2696 /*@ 2697 SNESDestroy - Destroys the nonlinear solver context that was created 2698 with SNESCreate(). 2699 2700 Collective on SNES 2701 2702 Input Parameter: 2703 . snes - the SNES context 2704 2705 Level: beginner 2706 2707 .keywords: SNES, nonlinear, destroy 2708 2709 .seealso: SNESCreate(), SNESSolve() 2710 @*/ 2711 PetscErrorCode SNESDestroy(SNES *snes) 2712 { 2713 PetscErrorCode ierr; 2714 2715 PetscFunctionBegin; 2716 if (!*snes) PetscFunctionReturn(0); 2717 PetscValidHeaderSpecific((*snes),SNES_CLASSID,1); 2718 if (--((PetscObject)(*snes))->refct > 0) {*snes = 0; PetscFunctionReturn(0);} 2719 2720 ierr = SNESReset((*snes));CHKERRQ(ierr); 2721 ierr = SNESDestroy(&(*snes)->pc);CHKERRQ(ierr); 2722 2723 /* if memory was published with SAWs then destroy it */ 2724 ierr = PetscObjectSAWsViewOff((PetscObject)*snes);CHKERRQ(ierr); 2725 if ((*snes)->ops->destroy) {ierr = (*((*snes))->ops->destroy)((*snes));CHKERRQ(ierr);} 2726 2727 ierr = DMDestroy(&(*snes)->dm);CHKERRQ(ierr); 2728 ierr = KSPDestroy(&(*snes)->ksp);CHKERRQ(ierr); 2729 ierr = SNESLineSearchDestroy(&(*snes)->linesearch);CHKERRQ(ierr); 2730 2731 ierr = PetscFree((*snes)->kspconvctx);CHKERRQ(ierr); 2732 if ((*snes)->ops->convergeddestroy) { 2733 ierr = (*(*snes)->ops->convergeddestroy)((*snes)->cnvP);CHKERRQ(ierr); 2734 } 2735 if ((*snes)->conv_malloc) { 2736 ierr = PetscFree((*snes)->conv_hist);CHKERRQ(ierr); 2737 ierr = PetscFree((*snes)->conv_hist_its);CHKERRQ(ierr); 2738 } 2739 ierr = SNESMonitorCancel((*snes));CHKERRQ(ierr); 2740 ierr = PetscHeaderDestroy(snes);CHKERRQ(ierr); 2741 PetscFunctionReturn(0); 2742 } 2743 2744 /* ----------- Routines to set solver parameters ---------- */ 2745 2746 #undef __FUNCT__ 2747 #define __FUNCT__ "SNESSetLagPreconditioner" 2748 /*@ 2749 SNESSetLagPreconditioner - Determines when the preconditioner is rebuilt in the nonlinear solve. 2750 2751 Logically Collective on SNES 2752 2753 Input Parameters: 2754 + snes - the SNES context 2755 - lag - -1 indicates NEVER rebuild, 1 means rebuild every time the Jacobian is computed within a single nonlinear solve, 2 means every second time 2756 the Jacobian is built etc. -2 indicates rebuild preconditioner at next chance but then never rebuild after that 2757 2758 Options Database Keys: 2759 . -snes_lag_preconditioner <lag> 2760 2761 Notes: 2762 The default is 1 2763 The preconditioner is ALWAYS built in the first iteration of a nonlinear solve unless lag is -1 2764 If -1 is used before the very first nonlinear solve the preconditioner is still built because there is no previous preconditioner to use 2765 2766 Level: intermediate 2767 2768 .keywords: SNES, nonlinear, set, convergence, tolerances 2769 2770 .seealso: SNESSetTrustRegionTolerance(), SNESGetLagPreconditioner(), SNESSetLagJacobian(), SNESGetLagJacobian() 2771 2772 @*/ 2773 PetscErrorCode SNESSetLagPreconditioner(SNES snes,PetscInt lag) 2774 { 2775 PetscFunctionBegin; 2776 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2777 if (lag < -2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Lag must be -2, -1, 1 or greater"); 2778 if (!lag) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Lag cannot be 0"); 2779 PetscValidLogicalCollectiveInt(snes,lag,2); 2780 snes->lagpreconditioner = lag; 2781 PetscFunctionReturn(0); 2782 } 2783 2784 #undef __FUNCT__ 2785 #define __FUNCT__ "SNESSetGridSequence" 2786 /*@ 2787 SNESSetGridSequence - sets the number of steps of grid sequencing that SNES does 2788 2789 Logically Collective on SNES 2790 2791 Input Parameters: 2792 + snes - the SNES context 2793 - steps - the number of refinements to do, defaults to 0 2794 2795 Options Database Keys: 2796 . -snes_grid_sequence <steps> 2797 2798 Level: intermediate 2799 2800 Notes: 2801 Use SNESGetSolution() to extract the fine grid solution after grid sequencing. 2802 2803 .keywords: SNES, nonlinear, set, convergence, tolerances 2804 2805 .seealso: SNESSetTrustRegionTolerance(), SNESGetLagPreconditioner(), SNESSetLagJacobian(), SNESGetLagJacobian() 2806 2807 @*/ 2808 PetscErrorCode SNESSetGridSequence(SNES snes,PetscInt steps) 2809 { 2810 PetscFunctionBegin; 2811 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2812 PetscValidLogicalCollectiveInt(snes,steps,2); 2813 snes->gridsequence = steps; 2814 PetscFunctionReturn(0); 2815 } 2816 2817 #undef __FUNCT__ 2818 #define __FUNCT__ "SNESGetLagPreconditioner" 2819 /*@ 2820 SNESGetLagPreconditioner - Indicates how often the preconditioner is rebuilt 2821 2822 Not Collective 2823 2824 Input Parameter: 2825 . snes - the SNES context 2826 2827 Output Parameter: 2828 . lag - -1 indicates NEVER rebuild, 1 means rebuild every time the Jacobian is computed within a single nonlinear solve, 2 means every second time 2829 the Jacobian is built etc. -2 indicates rebuild preconditioner at next chance but then never rebuild after that 2830 2831 Options Database Keys: 2832 . -snes_lag_preconditioner <lag> 2833 2834 Notes: 2835 The default is 1 2836 The preconditioner is ALWAYS built in the first iteration of a nonlinear solve unless lag is -1 2837 2838 Level: intermediate 2839 2840 .keywords: SNES, nonlinear, set, convergence, tolerances 2841 2842 .seealso: SNESSetTrustRegionTolerance(), SNESSetLagPreconditioner() 2843 2844 @*/ 2845 PetscErrorCode SNESGetLagPreconditioner(SNES snes,PetscInt *lag) 2846 { 2847 PetscFunctionBegin; 2848 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2849 *lag = snes->lagpreconditioner; 2850 PetscFunctionReturn(0); 2851 } 2852 2853 #undef __FUNCT__ 2854 #define __FUNCT__ "SNESSetLagJacobian" 2855 /*@ 2856 SNESSetLagJacobian - Determines when the Jacobian is rebuilt in the nonlinear solve. See SNESSetLagPreconditioner() for determining how 2857 often the preconditioner is rebuilt. 2858 2859 Logically Collective on SNES 2860 2861 Input Parameters: 2862 + snes - the SNES context 2863 - lag - -1 indicates NEVER rebuild, 1 means rebuild every time the Jacobian is computed within a single nonlinear solve, 2 means every second time 2864 the Jacobian is built etc. -2 means rebuild at next chance but then never again 2865 2866 Options Database Keys: 2867 . -snes_lag_jacobian <lag> 2868 2869 Notes: 2870 The default is 1 2871 The Jacobian is ALWAYS built in the first iteration of a nonlinear solve unless lag is -1 2872 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 2873 at the next Newton step but never again (unless it is reset to another value) 2874 2875 Level: intermediate 2876 2877 .keywords: SNES, nonlinear, set, convergence, tolerances 2878 2879 .seealso: SNESSetTrustRegionTolerance(), SNESGetLagPreconditioner(), SNESSetLagPreconditioner(), SNESGetLagJacobian() 2880 2881 @*/ 2882 PetscErrorCode SNESSetLagJacobian(SNES snes,PetscInt lag) 2883 { 2884 PetscFunctionBegin; 2885 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2886 if (lag < -2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Lag must be -2, -1, 1 or greater"); 2887 if (!lag) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Lag cannot be 0"); 2888 PetscValidLogicalCollectiveInt(snes,lag,2); 2889 snes->lagjacobian = lag; 2890 PetscFunctionReturn(0); 2891 } 2892 2893 #undef __FUNCT__ 2894 #define __FUNCT__ "SNESGetLagJacobian" 2895 /*@ 2896 SNESGetLagJacobian - Indicates how often the Jacobian is rebuilt. See SNESGetLagPreconditioner() to determine when the preconditioner is rebuilt 2897 2898 Not Collective 2899 2900 Input Parameter: 2901 . snes - the SNES context 2902 2903 Output Parameter: 2904 . lag - -1 indicates NEVER rebuild, 1 means rebuild every time the Jacobian is computed within a single nonlinear solve, 2 means every second time 2905 the Jacobian is built etc. 2906 2907 Options Database Keys: 2908 . -snes_lag_jacobian <lag> 2909 2910 Notes: 2911 The default is 1 2912 The jacobian is ALWAYS built in the first iteration of a nonlinear solve unless lag is -1 2913 2914 Level: intermediate 2915 2916 .keywords: SNES, nonlinear, set, convergence, tolerances 2917 2918 .seealso: SNESSetTrustRegionTolerance(), SNESSetLagJacobian(), SNESSetLagPreconditioner(), SNESGetLagPreconditioner() 2919 2920 @*/ 2921 PetscErrorCode SNESGetLagJacobian(SNES snes,PetscInt *lag) 2922 { 2923 PetscFunctionBegin; 2924 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2925 *lag = snes->lagjacobian; 2926 PetscFunctionReturn(0); 2927 } 2928 2929 #undef __FUNCT__ 2930 #define __FUNCT__ "SNESSetLagJacobianPersists" 2931 /*@ 2932 SNESSetLagJacobianPersists - Set whether or not the Jacobian lagging persists through multiple solves 2933 2934 Logically collective on SNES 2935 2936 Input Parameter: 2937 + snes - the SNES context 2938 - flg - jacobian lagging persists if true 2939 2940 Options Database Keys: 2941 . -snes_lag_jacobian_persists <flg> 2942 2943 Notes: This is useful both for nonlinear preconditioning, where it's appropriate to have the Jacobian be stale by 2944 several solves, and for implicit time-stepping, where Jacobian lagging in the inner nonlinear solve over several 2945 timesteps may present huge efficiency gains. 2946 2947 Level: developer 2948 2949 .keywords: SNES, nonlinear, lag 2950 2951 .seealso: SNESSetLagPreconditionerPersists(), SNESSetLagJacobian(), SNESGetLagJacobian(), SNESGetNPC() 2952 2953 @*/ 2954 PetscErrorCode SNESSetLagJacobianPersists(SNES snes,PetscBool flg) 2955 { 2956 PetscFunctionBegin; 2957 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2958 PetscValidLogicalCollectiveBool(snes,flg,2); 2959 snes->lagjac_persist = flg; 2960 PetscFunctionReturn(0); 2961 } 2962 2963 #undef __FUNCT__ 2964 #define __FUNCT__ "SNESSetLagPreconditionerPersists" 2965 /*@ 2966 SNESSetLagPreconditionerPersists - Set whether or not the preconditioner lagging persists through multiple solves 2967 2968 Logically Collective on SNES 2969 2970 Input Parameter: 2971 + snes - the SNES context 2972 - flg - preconditioner lagging persists if true 2973 2974 Options Database Keys: 2975 . -snes_lag_jacobian_persists <flg> 2976 2977 Notes: This is useful both for nonlinear preconditioning, where it's appropriate to have the preconditioner be stale 2978 by several solves, and for implicit time-stepping, where preconditioner lagging in the inner nonlinear solve over 2979 several timesteps may present huge efficiency gains. 2980 2981 Level: developer 2982 2983 .keywords: SNES, nonlinear, lag 2984 2985 .seealso: SNESSetLagJacobianPersists(), SNESSetLagJacobian(), SNESGetLagJacobian(), SNESGetNPC() 2986 2987 @*/ 2988 PetscErrorCode SNESSetLagPreconditionerPersists(SNES snes,PetscBool flg) 2989 { 2990 PetscFunctionBegin; 2991 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2992 PetscValidLogicalCollectiveBool(snes,flg,2); 2993 snes->lagpre_persist = flg; 2994 PetscFunctionReturn(0); 2995 } 2996 2997 #undef __FUNCT__ 2998 #define __FUNCT__ "SNESSetTolerances" 2999 /*@ 3000 SNESSetTolerances - Sets various parameters used in convergence tests. 3001 3002 Logically Collective on SNES 3003 3004 Input Parameters: 3005 + snes - the SNES context 3006 . abstol - absolute convergence tolerance 3007 . rtol - relative convergence tolerance 3008 . stol - convergence tolerance in terms of the norm of the change in the solution between steps, || delta x || < stol*|| x || 3009 . maxit - maximum number of iterations 3010 - maxf - maximum number of function evaluations 3011 3012 Options Database Keys: 3013 + -snes_atol <abstol> - Sets abstol 3014 . -snes_rtol <rtol> - Sets rtol 3015 . -snes_stol <stol> - Sets stol 3016 . -snes_max_it <maxit> - Sets maxit 3017 - -snes_max_funcs <maxf> - Sets maxf 3018 3019 Notes: 3020 The default maximum number of iterations is 50. 3021 The default maximum number of function evaluations is 1000. 3022 3023 Level: intermediate 3024 3025 .keywords: SNES, nonlinear, set, convergence, tolerances 3026 3027 .seealso: SNESSetTrustRegionTolerance() 3028 @*/ 3029 PetscErrorCode SNESSetTolerances(SNES snes,PetscReal abstol,PetscReal rtol,PetscReal stol,PetscInt maxit,PetscInt maxf) 3030 { 3031 PetscFunctionBegin; 3032 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3033 PetscValidLogicalCollectiveReal(snes,abstol,2); 3034 PetscValidLogicalCollectiveReal(snes,rtol,3); 3035 PetscValidLogicalCollectiveReal(snes,stol,4); 3036 PetscValidLogicalCollectiveInt(snes,maxit,5); 3037 PetscValidLogicalCollectiveInt(snes,maxf,6); 3038 3039 if (abstol != PETSC_DEFAULT) { 3040 if (abstol < 0.0) SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_OUTOFRANGE,"Absolute tolerance %g must be non-negative",(double)abstol); 3041 snes->abstol = abstol; 3042 } 3043 if (rtol != PETSC_DEFAULT) { 3044 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); 3045 snes->rtol = rtol; 3046 } 3047 if (stol != PETSC_DEFAULT) { 3048 if (stol < 0.0) SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_OUTOFRANGE,"Step tolerance %g must be non-negative",(double)stol); 3049 snes->stol = stol; 3050 } 3051 if (maxit != PETSC_DEFAULT) { 3052 if (maxit < 0) SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_OUTOFRANGE,"Maximum number of iterations %D must be non-negative",maxit); 3053 snes->max_its = maxit; 3054 } 3055 if (maxf != PETSC_DEFAULT) { 3056 if (maxf < 0) SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_OUTOFRANGE,"Maximum number of function evaluations %D must be non-negative",maxf); 3057 snes->max_funcs = maxf; 3058 } 3059 snes->tolerancesset = PETSC_TRUE; 3060 PetscFunctionReturn(0); 3061 } 3062 3063 #undef __FUNCT__ 3064 #define __FUNCT__ "SNESGetTolerances" 3065 /*@ 3066 SNESGetTolerances - Gets various parameters used in convergence tests. 3067 3068 Not Collective 3069 3070 Input Parameters: 3071 + snes - the SNES context 3072 . atol - absolute convergence tolerance 3073 . rtol - relative convergence tolerance 3074 . stol - convergence tolerance in terms of the norm 3075 of the change in the solution between steps 3076 . maxit - maximum number of iterations 3077 - maxf - maximum number of function evaluations 3078 3079 Notes: 3080 The user can specify NULL for any parameter that is not needed. 3081 3082 Level: intermediate 3083 3084 .keywords: SNES, nonlinear, get, convergence, tolerances 3085 3086 .seealso: SNESSetTolerances() 3087 @*/ 3088 PetscErrorCode SNESGetTolerances(SNES snes,PetscReal *atol,PetscReal *rtol,PetscReal *stol,PetscInt *maxit,PetscInt *maxf) 3089 { 3090 PetscFunctionBegin; 3091 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3092 if (atol) *atol = snes->abstol; 3093 if (rtol) *rtol = snes->rtol; 3094 if (stol) *stol = snes->stol; 3095 if (maxit) *maxit = snes->max_its; 3096 if (maxf) *maxf = snes->max_funcs; 3097 PetscFunctionReturn(0); 3098 } 3099 3100 #undef __FUNCT__ 3101 #define __FUNCT__ "SNESSetTrustRegionTolerance" 3102 /*@ 3103 SNESSetTrustRegionTolerance - Sets the trust region parameter tolerance. 3104 3105 Logically Collective on SNES 3106 3107 Input Parameters: 3108 + snes - the SNES context 3109 - tol - tolerance 3110 3111 Options Database Key: 3112 . -snes_trtol <tol> - Sets tol 3113 3114 Level: intermediate 3115 3116 .keywords: SNES, nonlinear, set, trust region, tolerance 3117 3118 .seealso: SNESSetTolerances() 3119 @*/ 3120 PetscErrorCode SNESSetTrustRegionTolerance(SNES snes,PetscReal tol) 3121 { 3122 PetscFunctionBegin; 3123 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3124 PetscValidLogicalCollectiveReal(snes,tol,2); 3125 snes->deltatol = tol; 3126 PetscFunctionReturn(0); 3127 } 3128 3129 /* 3130 Duplicate the lg monitors for SNES from KSP; for some reason with 3131 dynamic libraries things don't work under Sun4 if we just use 3132 macros instead of functions 3133 */ 3134 #undef __FUNCT__ 3135 #define __FUNCT__ "SNESMonitorLGResidualNorm" 3136 PetscErrorCode SNESMonitorLGResidualNorm(SNES snes,PetscInt it,PetscReal norm,PetscObject *objs) 3137 { 3138 PetscErrorCode ierr; 3139 3140 PetscFunctionBegin; 3141 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3142 ierr = KSPMonitorLGResidualNorm((KSP)snes,it,norm,objs);CHKERRQ(ierr); 3143 PetscFunctionReturn(0); 3144 } 3145 3146 #undef __FUNCT__ 3147 #define __FUNCT__ "SNESMonitorLGCreate" 3148 PetscErrorCode SNESMonitorLGCreate(const char host[],const char label[],int x,int y,int m,int n,PetscObject **draw) 3149 { 3150 PetscErrorCode ierr; 3151 3152 PetscFunctionBegin; 3153 ierr = KSPMonitorLGResidualNormCreate(host,label,x,y,m,n,draw);CHKERRQ(ierr); 3154 PetscFunctionReturn(0); 3155 } 3156 3157 #undef __FUNCT__ 3158 #define __FUNCT__ "SNESMonitorLGDestroy" 3159 PetscErrorCode SNESMonitorLGDestroy(PetscObject **objs) 3160 { 3161 PetscErrorCode ierr; 3162 3163 PetscFunctionBegin; 3164 ierr = KSPMonitorLGResidualNormDestroy(objs);CHKERRQ(ierr); 3165 PetscFunctionReturn(0); 3166 } 3167 3168 extern PetscErrorCode SNESMonitorRange_Private(SNES,PetscInt,PetscReal*); 3169 #undef __FUNCT__ 3170 #define __FUNCT__ "SNESMonitorLGRange" 3171 PetscErrorCode SNESMonitorLGRange(SNES snes,PetscInt n,PetscReal rnorm,void *monctx) 3172 { 3173 PetscDrawLG lg; 3174 PetscErrorCode ierr; 3175 PetscReal x,y,per; 3176 PetscViewer v = (PetscViewer)monctx; 3177 static PetscReal prev; /* should be in the context */ 3178 PetscDraw draw; 3179 3180 PetscFunctionBegin; 3181 ierr = PetscViewerDrawGetDrawLG(v,0,&lg);CHKERRQ(ierr); 3182 if (!n) {ierr = PetscDrawLGReset(lg);CHKERRQ(ierr);} 3183 ierr = PetscDrawLGGetDraw(lg,&draw);CHKERRQ(ierr); 3184 ierr = PetscDrawSetTitle(draw,"Residual norm");CHKERRQ(ierr); 3185 x = (PetscReal)n; 3186 if (rnorm > 0.0) y = PetscLog10Real(rnorm); 3187 else y = -15.0; 3188 ierr = PetscDrawLGAddPoint(lg,&x,&y);CHKERRQ(ierr); 3189 if (n < 20 || !(n % 5)) { 3190 ierr = PetscDrawLGDraw(lg);CHKERRQ(ierr); 3191 } 3192 3193 ierr = PetscViewerDrawGetDrawLG(v,1,&lg);CHKERRQ(ierr); 3194 if (!n) {ierr = PetscDrawLGReset(lg);CHKERRQ(ierr);} 3195 ierr = PetscDrawLGGetDraw(lg,&draw);CHKERRQ(ierr); 3196 ierr = PetscDrawSetTitle(draw,"% elemts > .2*max elemt");CHKERRQ(ierr); 3197 ierr = SNESMonitorRange_Private(snes,n,&per);CHKERRQ(ierr); 3198 x = (PetscReal)n; 3199 y = 100.0*per; 3200 ierr = PetscDrawLGAddPoint(lg,&x,&y);CHKERRQ(ierr); 3201 if (n < 20 || !(n % 5)) { 3202 ierr = PetscDrawLGDraw(lg);CHKERRQ(ierr); 3203 } 3204 3205 ierr = PetscViewerDrawGetDrawLG(v,2,&lg);CHKERRQ(ierr); 3206 if (!n) {prev = rnorm;ierr = PetscDrawLGReset(lg);CHKERRQ(ierr);} 3207 ierr = PetscDrawLGGetDraw(lg,&draw);CHKERRQ(ierr); 3208 ierr = PetscDrawSetTitle(draw,"(norm -oldnorm)/oldnorm");CHKERRQ(ierr); 3209 x = (PetscReal)n; 3210 y = (prev - rnorm)/prev; 3211 ierr = PetscDrawLGAddPoint(lg,&x,&y);CHKERRQ(ierr); 3212 if (n < 20 || !(n % 5)) { 3213 ierr = PetscDrawLGDraw(lg);CHKERRQ(ierr); 3214 } 3215 3216 ierr = PetscViewerDrawGetDrawLG(v,3,&lg);CHKERRQ(ierr); 3217 if (!n) {ierr = PetscDrawLGReset(lg);CHKERRQ(ierr);} 3218 ierr = PetscDrawLGGetDraw(lg,&draw);CHKERRQ(ierr); 3219 ierr = PetscDrawSetTitle(draw,"(norm -oldnorm)/oldnorm*(% > .2 max)");CHKERRQ(ierr); 3220 x = (PetscReal)n; 3221 y = (prev - rnorm)/(prev*per); 3222 if (n > 2) { /*skip initial crazy value */ 3223 ierr = PetscDrawLGAddPoint(lg,&x,&y);CHKERRQ(ierr); 3224 } 3225 if (n < 20 || !(n % 5)) { 3226 ierr = PetscDrawLGDraw(lg);CHKERRQ(ierr); 3227 } 3228 prev = rnorm; 3229 PetscFunctionReturn(0); 3230 } 3231 3232 #undef __FUNCT__ 3233 #define __FUNCT__ "SNESMonitor" 3234 /*@ 3235 SNESMonitor - runs the user provided monitor routines, if they exist 3236 3237 Collective on SNES 3238 3239 Input Parameters: 3240 + snes - nonlinear solver context obtained from SNESCreate() 3241 . iter - iteration number 3242 - rnorm - relative norm of the residual 3243 3244 Notes: 3245 This routine is called by the SNES implementations. 3246 It does not typically need to be called by the user. 3247 3248 Level: developer 3249 3250 .seealso: SNESMonitorSet() 3251 @*/ 3252 PetscErrorCode SNESMonitor(SNES snes,PetscInt iter,PetscReal rnorm) 3253 { 3254 PetscErrorCode ierr; 3255 PetscInt i,n = snes->numbermonitors; 3256 3257 PetscFunctionBegin; 3258 ierr = VecLockPush(snes->vec_sol);CHKERRQ(ierr); 3259 for (i=0; i<n; i++) { 3260 ierr = (*snes->monitor[i])(snes,iter,rnorm,snes->monitorcontext[i]);CHKERRQ(ierr); 3261 } 3262 ierr = VecLockPop(snes->vec_sol);CHKERRQ(ierr); 3263 PetscFunctionReturn(0); 3264 } 3265 3266 /* ------------ Routines to set performance monitoring options ----------- */ 3267 3268 /*MC 3269 SNESMonitorFunction - functional form passed to SNESMonitorSet() to monitor convergence of nonlinear solver 3270 3271 Synopsis: 3272 #include <petscsnes.h> 3273 $ PetscErrorCode SNESMonitorFunction(SNES snes,PetscInt its, PetscReal norm,void *mctx) 3274 3275 + snes - the SNES context 3276 . its - iteration number 3277 . norm - 2-norm function value (may be estimated) 3278 - mctx - [optional] monitoring context 3279 3280 Level: advanced 3281 3282 .seealso: SNESMonitorSet(), SNESMonitorGet() 3283 M*/ 3284 3285 #undef __FUNCT__ 3286 #define __FUNCT__ "SNESMonitorSet" 3287 /*@C 3288 SNESMonitorSet - Sets an ADDITIONAL function that is to be used at every 3289 iteration of the nonlinear solver to display the iteration's 3290 progress. 3291 3292 Logically Collective on SNES 3293 3294 Input Parameters: 3295 + snes - the SNES context 3296 . f - the monitor function, see SNESMonitorFunction for the calling sequence 3297 . mctx - [optional] user-defined context for private data for the 3298 monitor routine (use NULL if no context is desired) 3299 - monitordestroy - [optional] routine that frees monitor context 3300 (may be NULL) 3301 3302 Options Database Keys: 3303 + -snes_monitor - sets SNESMonitorDefault() 3304 . -snes_monitor_lg_residualnorm - sets line graph monitor, 3305 uses SNESMonitorLGCreate() 3306 - -snes_monitor_cancel - cancels all monitors that have 3307 been hardwired into a code by 3308 calls to SNESMonitorSet(), but 3309 does not cancel those set via 3310 the options database. 3311 3312 Notes: 3313 Several different monitoring routines may be set by calling 3314 SNESMonitorSet() multiple times; all will be called in the 3315 order in which they were set. 3316 3317 Fortran notes: Only a single monitor function can be set for each SNES object 3318 3319 Level: intermediate 3320 3321 .keywords: SNES, nonlinear, set, monitor 3322 3323 .seealso: SNESMonitorDefault(), SNESMonitorCancel(), SNESMonitorFunction 3324 @*/ 3325 PetscErrorCode SNESMonitorSet(SNES snes,PetscErrorCode (*f)(SNES,PetscInt,PetscReal,void*),void *mctx,PetscErrorCode (*monitordestroy)(void**)) 3326 { 3327 PetscInt i; 3328 PetscErrorCode ierr; 3329 3330 PetscFunctionBegin; 3331 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3332 if (snes->numbermonitors >= MAXSNESMONITORS) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Too many monitors set"); 3333 for (i=0; i<snes->numbermonitors;i++) { 3334 if (f == snes->monitor[i] && monitordestroy == snes->monitordestroy[i] && mctx == snes->monitorcontext[i]) { 3335 if (monitordestroy) { 3336 ierr = (*monitordestroy)(&mctx);CHKERRQ(ierr); 3337 } 3338 PetscFunctionReturn(0); 3339 } 3340 } 3341 snes->monitor[snes->numbermonitors] = f; 3342 snes->monitordestroy[snes->numbermonitors] = monitordestroy; 3343 snes->monitorcontext[snes->numbermonitors++] = (void*)mctx; 3344 PetscFunctionReturn(0); 3345 } 3346 3347 #undef __FUNCT__ 3348 #define __FUNCT__ "SNESMonitorCancel" 3349 /*@ 3350 SNESMonitorCancel - Clears all the monitor functions for a SNES object. 3351 3352 Logically Collective on SNES 3353 3354 Input Parameters: 3355 . snes - the SNES context 3356 3357 Options Database Key: 3358 . -snes_monitor_cancel - cancels all monitors that have been hardwired 3359 into a code by calls to SNESMonitorSet(), but does not cancel those 3360 set via the options database 3361 3362 Notes: 3363 There is no way to clear one specific monitor from a SNES object. 3364 3365 Level: intermediate 3366 3367 .keywords: SNES, nonlinear, set, monitor 3368 3369 .seealso: SNESMonitorDefault(), SNESMonitorSet() 3370 @*/ 3371 PetscErrorCode SNESMonitorCancel(SNES snes) 3372 { 3373 PetscErrorCode ierr; 3374 PetscInt i; 3375 3376 PetscFunctionBegin; 3377 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3378 for (i=0; i<snes->numbermonitors; i++) { 3379 if (snes->monitordestroy[i]) { 3380 ierr = (*snes->monitordestroy[i])(&snes->monitorcontext[i]);CHKERRQ(ierr); 3381 } 3382 } 3383 snes->numbermonitors = 0; 3384 PetscFunctionReturn(0); 3385 } 3386 3387 /*MC 3388 SNESConvergenceTestFunction - functional form used for testing of convergence of nonlinear solver 3389 3390 Synopsis: 3391 #include <petscsnes.h> 3392 $ PetscErrorCode SNESConvergenceTest(SNES snes,PetscInt it,PetscReal xnorm,PetscReal gnorm,PetscReal f,SNESConvergedReason *reason,void *cctx) 3393 3394 + snes - the SNES context 3395 . it - current iteration (0 is the first and is before any Newton step) 3396 . cctx - [optional] convergence context 3397 . reason - reason for convergence/divergence 3398 . xnorm - 2-norm of current iterate 3399 . gnorm - 2-norm of current step 3400 - f - 2-norm of function 3401 3402 Level: intermediate 3403 3404 .seealso: SNESSetConvergenceTest(), SNESGetConvergenceTest() 3405 M*/ 3406 3407 #undef __FUNCT__ 3408 #define __FUNCT__ "SNESSetConvergenceTest" 3409 /*@C 3410 SNESSetConvergenceTest - Sets the function that is to be used 3411 to test for convergence of the nonlinear iterative solution. 3412 3413 Logically Collective on SNES 3414 3415 Input Parameters: 3416 + snes - the SNES context 3417 . SNESConvergenceTestFunction - routine to test for convergence 3418 . cctx - [optional] context for private data for the convergence routine (may be NULL) 3419 - destroy - [optional] destructor for the context (may be NULL; NULL_FUNCTION in Fortran) 3420 3421 Level: advanced 3422 3423 .keywords: SNES, nonlinear, set, convergence, test 3424 3425 .seealso: SNESConvergedDefault(), SNESConvergedSkip(), SNESConvergenceTestFunction 3426 @*/ 3427 PetscErrorCode SNESSetConvergenceTest(SNES snes,PetscErrorCode (*SNESConvergenceTestFunction)(SNES,PetscInt,PetscReal,PetscReal,PetscReal,SNESConvergedReason*,void*),void *cctx,PetscErrorCode (*destroy)(void*)) 3428 { 3429 PetscErrorCode ierr; 3430 3431 PetscFunctionBegin; 3432 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3433 if (!SNESConvergenceTestFunction) SNESConvergenceTestFunction = SNESConvergedSkip; 3434 if (snes->ops->convergeddestroy) { 3435 ierr = (*snes->ops->convergeddestroy)(snes->cnvP);CHKERRQ(ierr); 3436 } 3437 snes->ops->converged = SNESConvergenceTestFunction; 3438 snes->ops->convergeddestroy = destroy; 3439 snes->cnvP = cctx; 3440 PetscFunctionReturn(0); 3441 } 3442 3443 #undef __FUNCT__ 3444 #define __FUNCT__ "SNESGetConvergedReason" 3445 /*@ 3446 SNESGetConvergedReason - Gets the reason the SNES iteration was stopped. 3447 3448 Not Collective 3449 3450 Input Parameter: 3451 . snes - the SNES context 3452 3453 Output Parameter: 3454 . reason - negative value indicates diverged, positive value converged, see SNESConvergedReason or the 3455 manual pages for the individual convergence tests for complete lists 3456 3457 Level: intermediate 3458 3459 Notes: Can only be called after the call the SNESSolve() is complete. 3460 3461 .keywords: SNES, nonlinear, set, convergence, test 3462 3463 .seealso: SNESSetConvergenceTest(), SNESConvergedReason 3464 @*/ 3465 PetscErrorCode SNESGetConvergedReason(SNES snes,SNESConvergedReason *reason) 3466 { 3467 PetscFunctionBegin; 3468 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3469 PetscValidPointer(reason,2); 3470 *reason = snes->reason; 3471 PetscFunctionReturn(0); 3472 } 3473 3474 #undef __FUNCT__ 3475 #define __FUNCT__ "SNESSetConvergenceHistory" 3476 /*@ 3477 SNESSetConvergenceHistory - Sets the array used to hold the convergence history. 3478 3479 Logically Collective on SNES 3480 3481 Input Parameters: 3482 + snes - iterative context obtained from SNESCreate() 3483 . a - array to hold history, this array will contain the function norms computed at each step 3484 . its - integer array holds the number of linear iterations for each solve. 3485 . na - size of a and its 3486 - reset - PETSC_TRUE indicates each new nonlinear solve resets the history counter to zero, 3487 else it continues storing new values for new nonlinear solves after the old ones 3488 3489 Notes: 3490 If 'a' and 'its' are NULL then space is allocated for the history. If 'na' PETSC_DECIDE or PETSC_DEFAULT then a 3491 default array of length 10000 is allocated. 3492 3493 This routine is useful, e.g., when running a code for purposes 3494 of accurate performance monitoring, when no I/O should be done 3495 during the section of code that is being timed. 3496 3497 Level: intermediate 3498 3499 .keywords: SNES, set, convergence, history 3500 3501 .seealso: SNESGetConvergenceHistory() 3502 3503 @*/ 3504 PetscErrorCode SNESSetConvergenceHistory(SNES snes,PetscReal a[],PetscInt its[],PetscInt na,PetscBool reset) 3505 { 3506 PetscErrorCode ierr; 3507 3508 PetscFunctionBegin; 3509 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3510 if (a) PetscValidScalarPointer(a,2); 3511 if (its) PetscValidIntPointer(its,3); 3512 if (!a) { 3513 if (na == PETSC_DECIDE || na == PETSC_DEFAULT) na = 1000; 3514 ierr = PetscCalloc1(na,&a);CHKERRQ(ierr); 3515 ierr = PetscCalloc1(na,&its);CHKERRQ(ierr); 3516 3517 snes->conv_malloc = PETSC_TRUE; 3518 } 3519 snes->conv_hist = a; 3520 snes->conv_hist_its = its; 3521 snes->conv_hist_max = na; 3522 snes->conv_hist_len = 0; 3523 snes->conv_hist_reset = reset; 3524 PetscFunctionReturn(0); 3525 } 3526 3527 #if defined(PETSC_HAVE_MATLAB_ENGINE) 3528 #include <engine.h> /* MATLAB include file */ 3529 #include <mex.h> /* MATLAB include file */ 3530 3531 #undef __FUNCT__ 3532 #define __FUNCT__ "SNESGetConvergenceHistoryMatlab" 3533 PETSC_EXTERN mxArray *SNESGetConvergenceHistoryMatlab(SNES snes) 3534 { 3535 mxArray *mat; 3536 PetscInt i; 3537 PetscReal *ar; 3538 3539 PetscFunctionBegin; 3540 mat = mxCreateDoubleMatrix(snes->conv_hist_len,1,mxREAL); 3541 ar = (PetscReal*) mxGetData(mat); 3542 for (i=0; i<snes->conv_hist_len; i++) ar[i] = snes->conv_hist[i]; 3543 PetscFunctionReturn(mat); 3544 } 3545 #endif 3546 3547 #undef __FUNCT__ 3548 #define __FUNCT__ "SNESGetConvergenceHistory" 3549 /*@C 3550 SNESGetConvergenceHistory - Gets the array used to hold the convergence history. 3551 3552 Not Collective 3553 3554 Input Parameter: 3555 . snes - iterative context obtained from SNESCreate() 3556 3557 Output Parameters: 3558 . a - array to hold history 3559 . its - integer array holds the number of linear iterations (or 3560 negative if not converged) for each solve. 3561 - na - size of a and its 3562 3563 Notes: 3564 The calling sequence for this routine in Fortran is 3565 $ call SNESGetConvergenceHistory(SNES snes, integer na, integer ierr) 3566 3567 This routine is useful, e.g., when running a code for purposes 3568 of accurate performance monitoring, when no I/O should be done 3569 during the section of code that is being timed. 3570 3571 Level: intermediate 3572 3573 .keywords: SNES, get, convergence, history 3574 3575 .seealso: SNESSetConvergencHistory() 3576 3577 @*/ 3578 PetscErrorCode SNESGetConvergenceHistory(SNES snes,PetscReal *a[],PetscInt *its[],PetscInt *na) 3579 { 3580 PetscFunctionBegin; 3581 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3582 if (a) *a = snes->conv_hist; 3583 if (its) *its = snes->conv_hist_its; 3584 if (na) *na = snes->conv_hist_len; 3585 PetscFunctionReturn(0); 3586 } 3587 3588 #undef __FUNCT__ 3589 #define __FUNCT__ "SNESSetUpdate" 3590 /*@C 3591 SNESSetUpdate - Sets the general-purpose update function called 3592 at the beginning of every iteration of the nonlinear solve. Specifically 3593 it is called just before the Jacobian is "evaluated". 3594 3595 Logically Collective on SNES 3596 3597 Input Parameters: 3598 . snes - The nonlinear solver context 3599 . func - The function 3600 3601 Calling sequence of func: 3602 . func (SNES snes, PetscInt step); 3603 3604 . step - The current step of the iteration 3605 3606 Level: advanced 3607 3608 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() 3609 This is not used by most users. 3610 3611 .keywords: SNES, update 3612 3613 .seealso SNESSetJacobian(), SNESSolve() 3614 @*/ 3615 PetscErrorCode SNESSetUpdate(SNES snes, PetscErrorCode (*func)(SNES, PetscInt)) 3616 { 3617 PetscFunctionBegin; 3618 PetscValidHeaderSpecific(snes, SNES_CLASSID,1); 3619 snes->ops->update = func; 3620 PetscFunctionReturn(0); 3621 } 3622 3623 #undef __FUNCT__ 3624 #define __FUNCT__ "SNESScaleStep_Private" 3625 /* 3626 SNESScaleStep_Private - Scales a step so that its length is less than the 3627 positive parameter delta. 3628 3629 Input Parameters: 3630 + snes - the SNES context 3631 . y - approximate solution of linear system 3632 . fnorm - 2-norm of current function 3633 - delta - trust region size 3634 3635 Output Parameters: 3636 + gpnorm - predicted function norm at the new point, assuming local 3637 linearization. The value is zero if the step lies within the trust 3638 region, and exceeds zero otherwise. 3639 - ynorm - 2-norm of the step 3640 3641 Note: 3642 For non-trust region methods such as SNESNEWTONLS, the parameter delta 3643 is set to be the maximum allowable step size. 3644 3645 .keywords: SNES, nonlinear, scale, step 3646 */ 3647 PetscErrorCode SNESScaleStep_Private(SNES snes,Vec y,PetscReal *fnorm,PetscReal *delta,PetscReal *gpnorm,PetscReal *ynorm) 3648 { 3649 PetscReal nrm; 3650 PetscScalar cnorm; 3651 PetscErrorCode ierr; 3652 3653 PetscFunctionBegin; 3654 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3655 PetscValidHeaderSpecific(y,VEC_CLASSID,2); 3656 PetscCheckSameComm(snes,1,y,2); 3657 3658 ierr = VecNorm(y,NORM_2,&nrm);CHKERRQ(ierr); 3659 if (nrm > *delta) { 3660 nrm = *delta/nrm; 3661 *gpnorm = (1.0 - nrm)*(*fnorm); 3662 cnorm = nrm; 3663 ierr = VecScale(y,cnorm);CHKERRQ(ierr); 3664 *ynorm = *delta; 3665 } else { 3666 *gpnorm = 0.0; 3667 *ynorm = nrm; 3668 } 3669 PetscFunctionReturn(0); 3670 } 3671 3672 #undef __FUNCT__ 3673 #define __FUNCT__ "SNESReasonView" 3674 /*@ 3675 SNESReasonView - Displays the reason a SNES solve converged or diverged to a viewer 3676 3677 Collective on SNES 3678 3679 Parameter: 3680 + snes - iterative context obtained from SNESCreate() 3681 - viewer - the viewer to display the reason 3682 3683 3684 Options Database Keys: 3685 . -snes_converged_reason - print reason for converged or diverged, also prints number of iterations 3686 3687 Level: beginner 3688 3689 .keywords: SNES, solve, linear system 3690 3691 .seealso: SNESCreate(), SNESSetUp(), SNESDestroy(), SNESSetTolerances(), SNESConvergedDefault() 3692 3693 @*/ 3694 PetscErrorCode SNESReasonView(SNES snes,PetscViewer viewer) 3695 { 3696 PetscErrorCode ierr; 3697 PetscBool isAscii; 3698 3699 PetscFunctionBegin; 3700 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isAscii);CHKERRQ(ierr); 3701 if (isAscii) { 3702 ierr = PetscViewerASCIIAddTab(viewer,((PetscObject)snes)->tablevel);CHKERRQ(ierr); 3703 if (snes->reason > 0) { 3704 if (((PetscObject) snes)->prefix) { 3705 ierr = PetscViewerASCIIPrintf(viewer,"Nonlinear %s solve converged due to %s iterations %D\n",((PetscObject) snes)->prefix,SNESConvergedReasons[snes->reason],snes->iter);CHKERRQ(ierr); 3706 } else { 3707 ierr = PetscViewerASCIIPrintf(viewer,"Nonlinear solve converged due to %s iterations %D\n",SNESConvergedReasons[snes->reason],snes->iter);CHKERRQ(ierr); 3708 } 3709 } else { 3710 if (((PetscObject) snes)->prefix) { 3711 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); 3712 } else { 3713 ierr = PetscViewerASCIIPrintf(viewer,"Nonlinear solve did not converge due to %s iterations %D\n",SNESConvergedReasons[snes->reason],snes->iter);CHKERRQ(ierr); 3714 } 3715 } 3716 ierr = PetscViewerASCIISubtractTab(viewer,((PetscObject)snes)->tablevel);CHKERRQ(ierr); 3717 } 3718 PetscFunctionReturn(0); 3719 } 3720 3721 #undef __FUNCT__ 3722 #define __FUNCT__ "SNESReasonViewFromOptions" 3723 /*@C 3724 SNESReasonViewFromOptions - Processes command line options to determine if/how a SNESReason is to be viewed. 3725 3726 Collective on SNES 3727 3728 Input Parameters: 3729 . snes - the SNES object 3730 3731 Level: intermediate 3732 3733 @*/ 3734 PetscErrorCode SNESReasonViewFromOptions(SNES snes) 3735 { 3736 PetscErrorCode ierr; 3737 PetscViewer viewer; 3738 PetscBool flg; 3739 static PetscBool incall = PETSC_FALSE; 3740 PetscViewerFormat format; 3741 3742 PetscFunctionBegin; 3743 if (incall) PetscFunctionReturn(0); 3744 incall = PETSC_TRUE; 3745 ierr = PetscOptionsGetViewer(PetscObjectComm((PetscObject)snes),((PetscObject)snes)->prefix,"-snes_converged_reason",&viewer,&format,&flg);CHKERRQ(ierr); 3746 if (flg) { 3747 ierr = PetscViewerPushFormat(viewer,format);CHKERRQ(ierr); 3748 ierr = SNESReasonView(snes,viewer);CHKERRQ(ierr); 3749 ierr = PetscViewerPopFormat(viewer);CHKERRQ(ierr); 3750 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 3751 } 3752 incall = PETSC_FALSE; 3753 PetscFunctionReturn(0); 3754 } 3755 3756 #undef __FUNCT__ 3757 #define __FUNCT__ "SNESSolve" 3758 /*@C 3759 SNESSolve - Solves a nonlinear system F(x) = b. 3760 Call SNESSolve() after calling SNESCreate() and optional routines of the form SNESSetXXX(). 3761 3762 Collective on SNES 3763 3764 Input Parameters: 3765 + snes - the SNES context 3766 . b - the constant part of the equation F(x) = b, or NULL to use zero. 3767 - x - the solution vector. 3768 3769 Notes: 3770 The user should initialize the vector,x, with the initial guess 3771 for the nonlinear solve prior to calling SNESSolve. In particular, 3772 to employ an initial guess of zero, the user should explicitly set 3773 this vector to zero by calling VecSet(). 3774 3775 Level: beginner 3776 3777 .keywords: SNES, nonlinear, solve 3778 3779 .seealso: SNESCreate(), SNESDestroy(), SNESSetFunction(), SNESSetJacobian(), SNESSetGridSequence(), SNESGetSolution() 3780 @*/ 3781 PetscErrorCode SNESSolve(SNES snes,Vec b,Vec x) 3782 { 3783 PetscErrorCode ierr; 3784 PetscBool flg; 3785 PetscInt grid; 3786 Vec xcreated = NULL; 3787 DM dm; 3788 3789 PetscFunctionBegin; 3790 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3791 if (x) PetscValidHeaderSpecific(x,VEC_CLASSID,3); 3792 if (x) PetscCheckSameComm(snes,1,x,3); 3793 if (b) PetscValidHeaderSpecific(b,VEC_CLASSID,2); 3794 if (b) PetscCheckSameComm(snes,1,b,2); 3795 3796 if (!x) { 3797 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 3798 ierr = DMCreateGlobalVector(dm,&xcreated);CHKERRQ(ierr); 3799 x = xcreated; 3800 } 3801 ierr = SNESViewFromOptions(snes,NULL,"-snes_view_pre");CHKERRQ(ierr); 3802 3803 for (grid=0; grid<snes->gridsequence; grid++) {ierr = PetscViewerASCIIPushTab(PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)snes)));CHKERRQ(ierr);} 3804 for (grid=0; grid<snes->gridsequence+1; grid++) { 3805 3806 /* set solution vector */ 3807 if (!grid) {ierr = PetscObjectReference((PetscObject)x);CHKERRQ(ierr);} 3808 ierr = VecDestroy(&snes->vec_sol);CHKERRQ(ierr); 3809 snes->vec_sol = x; 3810 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 3811 3812 /* set affine vector if provided */ 3813 if (b) { ierr = PetscObjectReference((PetscObject)b);CHKERRQ(ierr); } 3814 ierr = VecDestroy(&snes->vec_rhs);CHKERRQ(ierr); 3815 snes->vec_rhs = b; 3816 3817 if (snes->vec_func == snes->vec_sol) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_IDN,"Solution vector cannot be function vector"); 3818 if (snes->vec_rhs == snes->vec_sol) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_IDN,"Solution vector cannot be right hand side vector"); 3819 if (!snes->vec_sol_update /* && snes->vec_sol */) { 3820 ierr = VecDuplicate(snes->vec_sol,&snes->vec_sol_update);CHKERRQ(ierr); 3821 ierr = PetscLogObjectParent((PetscObject)snes,(PetscObject)snes->vec_sol_update);CHKERRQ(ierr); 3822 } 3823 ierr = DMShellSetGlobalVector(dm,snes->vec_sol);CHKERRQ(ierr); 3824 ierr = SNESSetUp(snes);CHKERRQ(ierr); 3825 3826 if (!grid) { 3827 if (snes->ops->computeinitialguess) { 3828 ierr = (*snes->ops->computeinitialguess)(snes,snes->vec_sol,snes->initialguessP);CHKERRQ(ierr); 3829 } 3830 } 3831 3832 if (snes->conv_hist_reset) snes->conv_hist_len = 0; 3833 if (snes->counters_reset) {snes->nfuncs = 0; snes->linear_its = 0; snes->numFailures = 0;} 3834 3835 ierr = PetscLogEventBegin(SNES_Solve,snes,0,0,0);CHKERRQ(ierr); 3836 ierr = (*snes->ops->solve)(snes);CHKERRQ(ierr); 3837 ierr = PetscLogEventEnd(SNES_Solve,snes,0,0,0);CHKERRQ(ierr); 3838 if (snes->domainerror) { 3839 snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; 3840 snes->domainerror = PETSC_FALSE; 3841 } 3842 if (!snes->reason) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Internal error, solver returned without setting converged reason"); 3843 3844 if (snes->lagjac_persist) snes->jac_iter += snes->iter; 3845 if (snes->lagpre_persist) snes->pre_iter += snes->iter; 3846 3847 flg = PETSC_FALSE; 3848 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_test_local_min",&flg,NULL);CHKERRQ(ierr); 3849 if (flg && !PetscPreLoadingOn) { ierr = SNESTestLocalMin(snes);CHKERRQ(ierr); } 3850 ierr = SNESReasonViewFromOptions(snes);CHKERRQ(ierr); 3851 3852 if (snes->errorifnotconverged && snes->reason < 0) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_NOT_CONVERGED,"SNESSolve has not converged"); 3853 if (grid < snes->gridsequence) { 3854 DM fine; 3855 Vec xnew; 3856 Mat interp; 3857 3858 ierr = DMRefine(snes->dm,PetscObjectComm((PetscObject)snes),&fine);CHKERRQ(ierr); 3859 if (!fine) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_INCOMP,"DMRefine() did not perform any refinement, cannot continue grid sequencing"); 3860 ierr = DMCreateInterpolation(snes->dm,fine,&interp,NULL);CHKERRQ(ierr); 3861 ierr = DMCreateGlobalVector(fine,&xnew);CHKERRQ(ierr); 3862 ierr = MatInterpolate(interp,x,xnew);CHKERRQ(ierr); 3863 ierr = DMInterpolate(snes->dm,interp,fine);CHKERRQ(ierr); 3864 ierr = MatDestroy(&interp);CHKERRQ(ierr); 3865 x = xnew; 3866 3867 ierr = SNESReset(snes);CHKERRQ(ierr); 3868 ierr = SNESSetDM(snes,fine);CHKERRQ(ierr); 3869 ierr = DMDestroy(&fine);CHKERRQ(ierr); 3870 ierr = PetscViewerASCIIPopTab(PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)snes)));CHKERRQ(ierr); 3871 } 3872 } 3873 ierr = SNESViewFromOptions(snes,NULL,"-snes_view");CHKERRQ(ierr); 3874 ierr = VecViewFromOptions(snes->vec_sol,((PetscObject)snes)->prefix,"-snes_view_solution");CHKERRQ(ierr); 3875 3876 ierr = VecDestroy(&xcreated);CHKERRQ(ierr); 3877 ierr = PetscObjectSAWsBlock((PetscObject)snes);CHKERRQ(ierr); 3878 PetscFunctionReturn(0); 3879 } 3880 3881 /* --------- Internal routines for SNES Package --------- */ 3882 3883 #undef __FUNCT__ 3884 #define __FUNCT__ "SNESSetType" 3885 /*@C 3886 SNESSetType - Sets the method for the nonlinear solver. 3887 3888 Collective on SNES 3889 3890 Input Parameters: 3891 + snes - the SNES context 3892 - type - a known method 3893 3894 Options Database Key: 3895 . -snes_type <type> - Sets the method; use -help for a list 3896 of available methods (for instance, newtonls or newtontr) 3897 3898 Notes: 3899 See "petsc/include/petscsnes.h" for available methods (for instance) 3900 + SNESNEWTONLS - Newton's method with line search 3901 (systems of nonlinear equations) 3902 . SNESNEWTONTR - Newton's method with trust region 3903 (systems of nonlinear equations) 3904 3905 Normally, it is best to use the SNESSetFromOptions() command and then 3906 set the SNES solver type from the options database rather than by using 3907 this routine. Using the options database provides the user with 3908 maximum flexibility in evaluating the many nonlinear solvers. 3909 The SNESSetType() routine is provided for those situations where it 3910 is necessary to set the nonlinear solver independently of the command 3911 line or options database. This might be the case, for example, when 3912 the choice of solver changes during the execution of the program, 3913 and the user's application is taking responsibility for choosing the 3914 appropriate method. 3915 3916 Developer Notes: SNESRegister() adds a constructor for a new SNESType to SNESList, SNESSetType() locates 3917 the constructor in that list and calls it to create the spexific object. 3918 3919 Level: intermediate 3920 3921 .keywords: SNES, set, type 3922 3923 .seealso: SNESType, SNESCreate(), SNESDestroy(), SNESGetType(), SNESSetFromOptions() 3924 3925 @*/ 3926 PetscErrorCode SNESSetType(SNES snes,SNESType type) 3927 { 3928 PetscErrorCode ierr,(*r)(SNES); 3929 PetscBool match; 3930 3931 PetscFunctionBegin; 3932 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3933 PetscValidCharPointer(type,2); 3934 3935 ierr = PetscObjectTypeCompare((PetscObject)snes,type,&match);CHKERRQ(ierr); 3936 if (match) PetscFunctionReturn(0); 3937 3938 ierr = PetscFunctionListFind(SNESList,type,&r);CHKERRQ(ierr); 3939 if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unable to find requested SNES type %s",type); 3940 /* Destroy the previous private SNES context */ 3941 if (snes->ops->destroy) { 3942 ierr = (*(snes)->ops->destroy)(snes);CHKERRQ(ierr); 3943 snes->ops->destroy = NULL; 3944 } 3945 /* Reinitialize function pointers in SNESOps structure */ 3946 snes->ops->setup = 0; 3947 snes->ops->solve = 0; 3948 snes->ops->view = 0; 3949 snes->ops->setfromoptions = 0; 3950 snes->ops->destroy = 0; 3951 /* Call the SNESCreate_XXX routine for this particular Nonlinear solver */ 3952 snes->setupcalled = PETSC_FALSE; 3953 3954 ierr = PetscObjectChangeTypeName((PetscObject)snes,type);CHKERRQ(ierr); 3955 ierr = (*r)(snes);CHKERRQ(ierr); 3956 PetscFunctionReturn(0); 3957 } 3958 3959 #undef __FUNCT__ 3960 #define __FUNCT__ "SNESGetType" 3961 /*@C 3962 SNESGetType - Gets the SNES method type and name (as a string). 3963 3964 Not Collective 3965 3966 Input Parameter: 3967 . snes - nonlinear solver context 3968 3969 Output Parameter: 3970 . type - SNES method (a character string) 3971 3972 Level: intermediate 3973 3974 .keywords: SNES, nonlinear, get, type, name 3975 @*/ 3976 PetscErrorCode SNESGetType(SNES snes,SNESType *type) 3977 { 3978 PetscFunctionBegin; 3979 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3980 PetscValidPointer(type,2); 3981 *type = ((PetscObject)snes)->type_name; 3982 PetscFunctionReturn(0); 3983 } 3984 3985 #undef __FUNCT__ 3986 #define __FUNCT__ "SNESSetSolution" 3987 /*@ 3988 SNESSetSolution - Sets the solution vector for use by the SNES routines. 3989 3990 Logically Collective on SNES and Vec 3991 3992 Input Parameters: 3993 + snes - the SNES context obtained from SNESCreate() 3994 - u - the solution vector 3995 3996 Level: beginner 3997 3998 .keywords: SNES, set, solution 3999 @*/ 4000 PetscErrorCode SNESSetSolution(SNES snes, Vec u) 4001 { 4002 DM dm; 4003 PetscErrorCode ierr; 4004 4005 PetscFunctionBegin; 4006 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); 4007 PetscValidHeaderSpecific(u, VEC_CLASSID, 2); 4008 ierr = PetscObjectReference((PetscObject) u);CHKERRQ(ierr); 4009 ierr = VecDestroy(&snes->vec_sol);CHKERRQ(ierr); 4010 4011 snes->vec_sol = u; 4012 4013 ierr = SNESGetDM(snes, &dm);CHKERRQ(ierr); 4014 ierr = DMShellSetGlobalVector(dm, u);CHKERRQ(ierr); 4015 PetscFunctionReturn(0); 4016 } 4017 4018 #undef __FUNCT__ 4019 #define __FUNCT__ "SNESGetSolution" 4020 /*@ 4021 SNESGetSolution - Returns the vector where the approximate solution is 4022 stored. This is the fine grid solution when using SNESSetGridSequence(). 4023 4024 Not Collective, but Vec is parallel if SNES is parallel 4025 4026 Input Parameter: 4027 . snes - the SNES context 4028 4029 Output Parameter: 4030 . x - the solution 4031 4032 Level: intermediate 4033 4034 .keywords: SNES, nonlinear, get, solution 4035 4036 .seealso: SNESGetSolutionUpdate(), SNESGetFunction() 4037 @*/ 4038 PetscErrorCode SNESGetSolution(SNES snes,Vec *x) 4039 { 4040 PetscFunctionBegin; 4041 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4042 PetscValidPointer(x,2); 4043 *x = snes->vec_sol; 4044 PetscFunctionReturn(0); 4045 } 4046 4047 #undef __FUNCT__ 4048 #define __FUNCT__ "SNESGetSolutionUpdate" 4049 /*@ 4050 SNESGetSolutionUpdate - Returns the vector where the solution update is 4051 stored. 4052 4053 Not Collective, but Vec is parallel if SNES is parallel 4054 4055 Input Parameter: 4056 . snes - the SNES context 4057 4058 Output Parameter: 4059 . x - the solution update 4060 4061 Level: advanced 4062 4063 .keywords: SNES, nonlinear, get, solution, update 4064 4065 .seealso: SNESGetSolution(), SNESGetFunction() 4066 @*/ 4067 PetscErrorCode SNESGetSolutionUpdate(SNES snes,Vec *x) 4068 { 4069 PetscFunctionBegin; 4070 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4071 PetscValidPointer(x,2); 4072 *x = snes->vec_sol_update; 4073 PetscFunctionReturn(0); 4074 } 4075 4076 #undef __FUNCT__ 4077 #define __FUNCT__ "SNESGetFunction" 4078 /*@C 4079 SNESGetFunction - Returns the vector where the function is stored. 4080 4081 Not Collective, but Vec is parallel if SNES is parallel. Collective if Vec is requested, but has not been created yet. 4082 4083 Input Parameter: 4084 . snes - the SNES context 4085 4086 Output Parameter: 4087 + r - the vector that is used to store residuals (or NULL if you don't want it) 4088 . f - the function (or NULL if you don't want it); see SNESFunction for calling sequence details 4089 - ctx - the function context (or NULL if you don't want it) 4090 4091 Level: advanced 4092 4093 .keywords: SNES, nonlinear, get, function 4094 4095 .seealso: SNESSetFunction(), SNESGetSolution(), SNESFunction 4096 @*/ 4097 PetscErrorCode SNESGetFunction(SNES snes,Vec *r,PetscErrorCode (**f)(SNES,Vec,Vec,void*),void **ctx) 4098 { 4099 PetscErrorCode ierr; 4100 DM dm; 4101 4102 PetscFunctionBegin; 4103 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4104 if (r) { 4105 if (!snes->vec_func) { 4106 if (snes->vec_rhs) { 4107 ierr = VecDuplicate(snes->vec_rhs,&snes->vec_func);CHKERRQ(ierr); 4108 } else if (snes->vec_sol) { 4109 ierr = VecDuplicate(snes->vec_sol,&snes->vec_func);CHKERRQ(ierr); 4110 } else if (snes->dm) { 4111 ierr = DMCreateGlobalVector(snes->dm,&snes->vec_func);CHKERRQ(ierr); 4112 } 4113 } 4114 *r = snes->vec_func; 4115 } 4116 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 4117 ierr = DMSNESGetFunction(dm,f,ctx);CHKERRQ(ierr); 4118 PetscFunctionReturn(0); 4119 } 4120 4121 /*@C 4122 SNESGetNGS - Returns the NGS function and context. 4123 4124 Input Parameter: 4125 . snes - the SNES context 4126 4127 Output Parameter: 4128 + f - the function (or NULL) see SNESNGSFunction for details 4129 - ctx - the function context (or NULL) 4130 4131 Level: advanced 4132 4133 .keywords: SNES, nonlinear, get, function 4134 4135 .seealso: SNESSetNGS(), SNESGetFunction() 4136 @*/ 4137 4138 #undef __FUNCT__ 4139 #define __FUNCT__ "SNESGetNGS" 4140 PetscErrorCode SNESGetNGS (SNES snes, PetscErrorCode (**f)(SNES, Vec, Vec, void*), void ** ctx) 4141 { 4142 PetscErrorCode ierr; 4143 DM dm; 4144 4145 PetscFunctionBegin; 4146 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4147 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 4148 ierr = DMSNESGetNGS(dm,f,ctx);CHKERRQ(ierr); 4149 PetscFunctionReturn(0); 4150 } 4151 4152 #undef __FUNCT__ 4153 #define __FUNCT__ "SNESSetOptionsPrefix" 4154 /*@C 4155 SNESSetOptionsPrefix - Sets the prefix used for searching for all 4156 SNES options in the database. 4157 4158 Logically Collective on SNES 4159 4160 Input Parameter: 4161 + snes - the SNES context 4162 - prefix - the prefix to prepend to all option names 4163 4164 Notes: 4165 A hyphen (-) must NOT be given at the beginning of the prefix name. 4166 The first character of all runtime options is AUTOMATICALLY the hyphen. 4167 4168 Level: advanced 4169 4170 .keywords: SNES, set, options, prefix, database 4171 4172 .seealso: SNESSetFromOptions() 4173 @*/ 4174 PetscErrorCode SNESSetOptionsPrefix(SNES snes,const char prefix[]) 4175 { 4176 PetscErrorCode ierr; 4177 4178 PetscFunctionBegin; 4179 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4180 ierr = PetscObjectSetOptionsPrefix((PetscObject)snes,prefix);CHKERRQ(ierr); 4181 if (!snes->ksp) {ierr = SNESGetKSP(snes,&snes->ksp);CHKERRQ(ierr);} 4182 if (snes->linesearch) { 4183 ierr = SNESGetLineSearch(snes,&snes->linesearch);CHKERRQ(ierr); 4184 ierr = PetscObjectSetOptionsPrefix((PetscObject)snes->linesearch,prefix);CHKERRQ(ierr); 4185 } 4186 ierr = KSPSetOptionsPrefix(snes->ksp,prefix);CHKERRQ(ierr); 4187 PetscFunctionReturn(0); 4188 } 4189 4190 #undef __FUNCT__ 4191 #define __FUNCT__ "SNESAppendOptionsPrefix" 4192 /*@C 4193 SNESAppendOptionsPrefix - Appends to the prefix used for searching for all 4194 SNES options in the database. 4195 4196 Logically Collective on SNES 4197 4198 Input Parameters: 4199 + snes - the SNES context 4200 - prefix - the prefix to prepend to all option names 4201 4202 Notes: 4203 A hyphen (-) must NOT be given at the beginning of the prefix name. 4204 The first character of all runtime options is AUTOMATICALLY the hyphen. 4205 4206 Level: advanced 4207 4208 .keywords: SNES, append, options, prefix, database 4209 4210 .seealso: SNESGetOptionsPrefix() 4211 @*/ 4212 PetscErrorCode SNESAppendOptionsPrefix(SNES snes,const char prefix[]) 4213 { 4214 PetscErrorCode ierr; 4215 4216 PetscFunctionBegin; 4217 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4218 ierr = PetscObjectAppendOptionsPrefix((PetscObject)snes,prefix);CHKERRQ(ierr); 4219 if (!snes->ksp) {ierr = SNESGetKSP(snes,&snes->ksp);CHKERRQ(ierr);} 4220 if (snes->linesearch) { 4221 ierr = SNESGetLineSearch(snes,&snes->linesearch);CHKERRQ(ierr); 4222 ierr = PetscObjectAppendOptionsPrefix((PetscObject)snes->linesearch,prefix);CHKERRQ(ierr); 4223 } 4224 ierr = KSPAppendOptionsPrefix(snes->ksp,prefix);CHKERRQ(ierr); 4225 PetscFunctionReturn(0); 4226 } 4227 4228 #undef __FUNCT__ 4229 #define __FUNCT__ "SNESGetOptionsPrefix" 4230 /*@C 4231 SNESGetOptionsPrefix - Sets the prefix used for searching for all 4232 SNES options in the database. 4233 4234 Not Collective 4235 4236 Input Parameter: 4237 . snes - the SNES context 4238 4239 Output Parameter: 4240 . prefix - pointer to the prefix string used 4241 4242 Notes: On the fortran side, the user should pass in a string 'prefix' of 4243 sufficient length to hold the prefix. 4244 4245 Level: advanced 4246 4247 .keywords: SNES, get, options, prefix, database 4248 4249 .seealso: SNESAppendOptionsPrefix() 4250 @*/ 4251 PetscErrorCode SNESGetOptionsPrefix(SNES snes,const char *prefix[]) 4252 { 4253 PetscErrorCode ierr; 4254 4255 PetscFunctionBegin; 4256 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4257 ierr = PetscObjectGetOptionsPrefix((PetscObject)snes,prefix);CHKERRQ(ierr); 4258 PetscFunctionReturn(0); 4259 } 4260 4261 4262 #undef __FUNCT__ 4263 #define __FUNCT__ "SNESRegister" 4264 /*@C 4265 SNESRegister - Adds a method to the nonlinear solver package. 4266 4267 Not collective 4268 4269 Input Parameters: 4270 + name_solver - name of a new user-defined solver 4271 - routine_create - routine to create method context 4272 4273 Notes: 4274 SNESRegister() may be called multiple times to add several user-defined solvers. 4275 4276 Sample usage: 4277 .vb 4278 SNESRegister("my_solver",MySolverCreate); 4279 .ve 4280 4281 Then, your solver can be chosen with the procedural interface via 4282 $ SNESSetType(snes,"my_solver") 4283 or at runtime via the option 4284 $ -snes_type my_solver 4285 4286 Level: advanced 4287 4288 Note: If your function is not being put into a shared library then use SNESRegister() instead 4289 4290 .keywords: SNES, nonlinear, register 4291 4292 .seealso: SNESRegisterAll(), SNESRegisterDestroy() 4293 4294 Level: advanced 4295 @*/ 4296 PetscErrorCode SNESRegister(const char sname[],PetscErrorCode (*function)(SNES)) 4297 { 4298 PetscErrorCode ierr; 4299 4300 PetscFunctionBegin; 4301 ierr = PetscFunctionListAdd(&SNESList,sname,function);CHKERRQ(ierr); 4302 PetscFunctionReturn(0); 4303 } 4304 4305 #undef __FUNCT__ 4306 #define __FUNCT__ "SNESTestLocalMin" 4307 PetscErrorCode SNESTestLocalMin(SNES snes) 4308 { 4309 PetscErrorCode ierr; 4310 PetscInt N,i,j; 4311 Vec u,uh,fh; 4312 PetscScalar value; 4313 PetscReal norm; 4314 4315 PetscFunctionBegin; 4316 ierr = SNESGetSolution(snes,&u);CHKERRQ(ierr); 4317 ierr = VecDuplicate(u,&uh);CHKERRQ(ierr); 4318 ierr = VecDuplicate(u,&fh);CHKERRQ(ierr); 4319 4320 /* currently only works for sequential */ 4321 ierr = PetscPrintf(PETSC_COMM_WORLD,"Testing FormFunction() for local min\n");CHKERRQ(ierr); 4322 ierr = VecGetSize(u,&N);CHKERRQ(ierr); 4323 for (i=0; i<N; i++) { 4324 ierr = VecCopy(u,uh);CHKERRQ(ierr); 4325 ierr = PetscPrintf(PETSC_COMM_WORLD,"i = %D\n",i);CHKERRQ(ierr); 4326 for (j=-10; j<11; j++) { 4327 value = PetscSign(j)*PetscExpReal(PetscAbs(j)-10.0); 4328 ierr = VecSetValue(uh,i,value,ADD_VALUES);CHKERRQ(ierr); 4329 ierr = SNESComputeFunction(snes,uh,fh);CHKERRQ(ierr); 4330 ierr = VecNorm(fh,NORM_2,&norm);CHKERRQ(ierr); 4331 ierr = PetscPrintf(PETSC_COMM_WORLD," j norm %D %18.16e\n",j,norm);CHKERRQ(ierr); 4332 value = -value; 4333 ierr = VecSetValue(uh,i,value,ADD_VALUES);CHKERRQ(ierr); 4334 } 4335 } 4336 ierr = VecDestroy(&uh);CHKERRQ(ierr); 4337 ierr = VecDestroy(&fh);CHKERRQ(ierr); 4338 PetscFunctionReturn(0); 4339 } 4340 4341 #undef __FUNCT__ 4342 #define __FUNCT__ "SNESKSPSetUseEW" 4343 /*@ 4344 SNESKSPSetUseEW - Sets SNES use Eisenstat-Walker method for 4345 computing relative tolerance for linear solvers within an inexact 4346 Newton method. 4347 4348 Logically Collective on SNES 4349 4350 Input Parameters: 4351 + snes - SNES context 4352 - flag - PETSC_TRUE or PETSC_FALSE 4353 4354 Options Database: 4355 + -snes_ksp_ew - use Eisenstat-Walker method for determining linear system convergence 4356 . -snes_ksp_ew_version ver - version of Eisenstat-Walker method 4357 . -snes_ksp_ew_rtol0 <rtol0> - Sets rtol0 4358 . -snes_ksp_ew_rtolmax <rtolmax> - Sets rtolmax 4359 . -snes_ksp_ew_gamma <gamma> - Sets gamma 4360 . -snes_ksp_ew_alpha <alpha> - Sets alpha 4361 . -snes_ksp_ew_alpha2 <alpha2> - Sets alpha2 4362 - -snes_ksp_ew_threshold <threshold> - Sets threshold 4363 4364 Notes: 4365 Currently, the default is to use a constant relative tolerance for 4366 the inner linear solvers. Alternatively, one can use the 4367 Eisenstat-Walker method, where the relative convergence tolerance 4368 is reset at each Newton iteration according progress of the nonlinear 4369 solver. 4370 4371 Level: advanced 4372 4373 Reference: 4374 S. C. Eisenstat and H. F. Walker, "Choosing the forcing terms in an 4375 inexact Newton method", SISC 17 (1), pp.16-32, 1996. 4376 4377 .keywords: SNES, KSP, Eisenstat, Walker, convergence, test, inexact, Newton 4378 4379 .seealso: SNESKSPGetUseEW(), SNESKSPGetParametersEW(), SNESKSPSetParametersEW() 4380 @*/ 4381 PetscErrorCode SNESKSPSetUseEW(SNES snes,PetscBool flag) 4382 { 4383 PetscFunctionBegin; 4384 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4385 PetscValidLogicalCollectiveBool(snes,flag,2); 4386 snes->ksp_ewconv = flag; 4387 PetscFunctionReturn(0); 4388 } 4389 4390 #undef __FUNCT__ 4391 #define __FUNCT__ "SNESKSPGetUseEW" 4392 /*@ 4393 SNESKSPGetUseEW - Gets if SNES is using Eisenstat-Walker method 4394 for computing relative tolerance for linear solvers within an 4395 inexact Newton method. 4396 4397 Not Collective 4398 4399 Input Parameter: 4400 . snes - SNES context 4401 4402 Output Parameter: 4403 . flag - PETSC_TRUE or PETSC_FALSE 4404 4405 Notes: 4406 Currently, the default is to use a constant relative tolerance for 4407 the inner linear solvers. Alternatively, one can use the 4408 Eisenstat-Walker method, where the relative convergence tolerance 4409 is reset at each Newton iteration according progress of the nonlinear 4410 solver. 4411 4412 Level: advanced 4413 4414 Reference: 4415 S. C. Eisenstat and H. F. Walker, "Choosing the forcing terms in an 4416 inexact Newton method", SISC 17 (1), pp.16-32, 1996. 4417 4418 .keywords: SNES, KSP, Eisenstat, Walker, convergence, test, inexact, Newton 4419 4420 .seealso: SNESKSPSetUseEW(), SNESKSPGetParametersEW(), SNESKSPSetParametersEW() 4421 @*/ 4422 PetscErrorCode SNESKSPGetUseEW(SNES snes, PetscBool *flag) 4423 { 4424 PetscFunctionBegin; 4425 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4426 PetscValidPointer(flag,2); 4427 *flag = snes->ksp_ewconv; 4428 PetscFunctionReturn(0); 4429 } 4430 4431 #undef __FUNCT__ 4432 #define __FUNCT__ "SNESKSPSetParametersEW" 4433 /*@ 4434 SNESKSPSetParametersEW - Sets parameters for Eisenstat-Walker 4435 convergence criteria for the linear solvers within an inexact 4436 Newton method. 4437 4438 Logically Collective on SNES 4439 4440 Input Parameters: 4441 + snes - SNES context 4442 . version - version 1, 2 (default is 2) or 3 4443 . rtol_0 - initial relative tolerance (0 <= rtol_0 < 1) 4444 . rtol_max - maximum relative tolerance (0 <= rtol_max < 1) 4445 . gamma - multiplicative factor for version 2 rtol computation 4446 (0 <= gamma2 <= 1) 4447 . alpha - power for version 2 rtol computation (1 < alpha <= 2) 4448 . alpha2 - power for safeguard 4449 - threshold - threshold for imposing safeguard (0 < threshold < 1) 4450 4451 Note: 4452 Version 3 was contributed by Luis Chacon, June 2006. 4453 4454 Use PETSC_DEFAULT to retain the default for any of the parameters. 4455 4456 Level: advanced 4457 4458 Reference: 4459 S. C. Eisenstat and H. F. Walker, "Choosing the forcing terms in an 4460 inexact Newton method", Utah State University Math. Stat. Dept. Res. 4461 Report 6/94/75, June, 1994, to appear in SIAM J. Sci. Comput. 4462 4463 .keywords: SNES, KSP, Eisenstat, Walker, set, parameters 4464 4465 .seealso: SNESKSPSetUseEW(), SNESKSPGetUseEW(), SNESKSPGetParametersEW() 4466 @*/ 4467 PetscErrorCode SNESKSPSetParametersEW(SNES snes,PetscInt version,PetscReal rtol_0,PetscReal rtol_max,PetscReal gamma,PetscReal alpha,PetscReal alpha2,PetscReal threshold) 4468 { 4469 SNESKSPEW *kctx; 4470 4471 PetscFunctionBegin; 4472 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4473 kctx = (SNESKSPEW*)snes->kspconvctx; 4474 if (!kctx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"No Eisenstat-Walker context existing"); 4475 PetscValidLogicalCollectiveInt(snes,version,2); 4476 PetscValidLogicalCollectiveReal(snes,rtol_0,3); 4477 PetscValidLogicalCollectiveReal(snes,rtol_max,4); 4478 PetscValidLogicalCollectiveReal(snes,gamma,5); 4479 PetscValidLogicalCollectiveReal(snes,alpha,6); 4480 PetscValidLogicalCollectiveReal(snes,alpha2,7); 4481 PetscValidLogicalCollectiveReal(snes,threshold,8); 4482 4483 if (version != PETSC_DEFAULT) kctx->version = version; 4484 if (rtol_0 != PETSC_DEFAULT) kctx->rtol_0 = rtol_0; 4485 if (rtol_max != PETSC_DEFAULT) kctx->rtol_max = rtol_max; 4486 if (gamma != PETSC_DEFAULT) kctx->gamma = gamma; 4487 if (alpha != PETSC_DEFAULT) kctx->alpha = alpha; 4488 if (alpha2 != PETSC_DEFAULT) kctx->alpha2 = alpha2; 4489 if (threshold != PETSC_DEFAULT) kctx->threshold = threshold; 4490 4491 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); 4492 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); 4493 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); 4494 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); 4495 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); 4496 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); 4497 PetscFunctionReturn(0); 4498 } 4499 4500 #undef __FUNCT__ 4501 #define __FUNCT__ "SNESKSPGetParametersEW" 4502 /*@ 4503 SNESKSPGetParametersEW - Gets parameters for Eisenstat-Walker 4504 convergence criteria for the linear solvers within an inexact 4505 Newton method. 4506 4507 Not Collective 4508 4509 Input Parameters: 4510 snes - SNES context 4511 4512 Output Parameters: 4513 + version - version 1, 2 (default is 2) or 3 4514 . rtol_0 - initial relative tolerance (0 <= rtol_0 < 1) 4515 . rtol_max - maximum relative tolerance (0 <= rtol_max < 1) 4516 . gamma - multiplicative factor for version 2 rtol computation (0 <= gamma2 <= 1) 4517 . alpha - power for version 2 rtol computation (1 < alpha <= 2) 4518 . alpha2 - power for safeguard 4519 - threshold - threshold for imposing safeguard (0 < threshold < 1) 4520 4521 Level: advanced 4522 4523 .keywords: SNES, KSP, Eisenstat, Walker, get, parameters 4524 4525 .seealso: SNESKSPSetUseEW(), SNESKSPGetUseEW(), SNESKSPSetParametersEW() 4526 @*/ 4527 PetscErrorCode SNESKSPGetParametersEW(SNES snes,PetscInt *version,PetscReal *rtol_0,PetscReal *rtol_max,PetscReal *gamma,PetscReal *alpha,PetscReal *alpha2,PetscReal *threshold) 4528 { 4529 SNESKSPEW *kctx; 4530 4531 PetscFunctionBegin; 4532 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4533 kctx = (SNESKSPEW*)snes->kspconvctx; 4534 if (!kctx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"No Eisenstat-Walker context existing"); 4535 if (version) *version = kctx->version; 4536 if (rtol_0) *rtol_0 = kctx->rtol_0; 4537 if (rtol_max) *rtol_max = kctx->rtol_max; 4538 if (gamma) *gamma = kctx->gamma; 4539 if (alpha) *alpha = kctx->alpha; 4540 if (alpha2) *alpha2 = kctx->alpha2; 4541 if (threshold) *threshold = kctx->threshold; 4542 PetscFunctionReturn(0); 4543 } 4544 4545 #undef __FUNCT__ 4546 #define __FUNCT__ "KSPPreSolve_SNESEW" 4547 PetscErrorCode KSPPreSolve_SNESEW(KSP ksp, Vec b, Vec x, SNES snes) 4548 { 4549 PetscErrorCode ierr; 4550 SNESKSPEW *kctx = (SNESKSPEW*)snes->kspconvctx; 4551 PetscReal rtol = PETSC_DEFAULT,stol; 4552 4553 PetscFunctionBegin; 4554 if (!snes->ksp_ewconv) PetscFunctionReturn(0); 4555 if (!snes->iter) { 4556 rtol = kctx->rtol_0; /* first time in, so use the original user rtol */ 4557 ierr = VecNorm(snes->vec_func,NORM_2,&kctx->norm_first);CHKERRQ(ierr); 4558 } 4559 else { 4560 if (kctx->version == 1) { 4561 rtol = (snes->norm - kctx->lresid_last)/kctx->norm_last; 4562 if (rtol < 0.0) rtol = -rtol; 4563 stol = PetscPowReal(kctx->rtol_last,kctx->alpha2); 4564 if (stol > kctx->threshold) rtol = PetscMax(rtol,stol); 4565 } else if (kctx->version == 2) { 4566 rtol = kctx->gamma * PetscPowReal(snes->norm/kctx->norm_last,kctx->alpha); 4567 stol = kctx->gamma * PetscPowReal(kctx->rtol_last,kctx->alpha); 4568 if (stol > kctx->threshold) rtol = PetscMax(rtol,stol); 4569 } else if (kctx->version == 3) { /* contributed by Luis Chacon, June 2006. */ 4570 rtol = kctx->gamma * PetscPowReal(snes->norm/kctx->norm_last,kctx->alpha); 4571 /* safeguard: avoid sharp decrease of rtol */ 4572 stol = kctx->gamma*PetscPowReal(kctx->rtol_last,kctx->alpha); 4573 stol = PetscMax(rtol,stol); 4574 rtol = PetscMin(kctx->rtol_0,stol); 4575 /* safeguard: avoid oversolving */ 4576 stol = kctx->gamma*(kctx->norm_first*snes->rtol)/snes->norm; 4577 stol = PetscMax(rtol,stol); 4578 rtol = PetscMin(kctx->rtol_0,stol); 4579 } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Only versions 1, 2 or 3 are supported: %D",kctx->version); 4580 } 4581 /* safeguard: avoid rtol greater than one */ 4582 rtol = PetscMin(rtol,kctx->rtol_max); 4583 ierr = KSPSetTolerances(ksp,rtol,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);CHKERRQ(ierr); 4584 ierr = PetscInfo3(snes,"iter %D, Eisenstat-Walker (version %D) KSP rtol=%g\n",snes->iter,kctx->version,(double)rtol);CHKERRQ(ierr); 4585 PetscFunctionReturn(0); 4586 } 4587 4588 #undef __FUNCT__ 4589 #define __FUNCT__ "KSPPostSolve_SNESEW" 4590 PetscErrorCode KSPPostSolve_SNESEW(KSP ksp, Vec b, Vec x, SNES snes) 4591 { 4592 PetscErrorCode ierr; 4593 SNESKSPEW *kctx = (SNESKSPEW*)snes->kspconvctx; 4594 PCSide pcside; 4595 Vec lres; 4596 4597 PetscFunctionBegin; 4598 if (!snes->ksp_ewconv) PetscFunctionReturn(0); 4599 ierr = KSPGetTolerances(ksp,&kctx->rtol_last,0,0,0);CHKERRQ(ierr); 4600 kctx->norm_last = snes->norm; 4601 if (kctx->version == 1) { 4602 PC pc; 4603 PetscBool isNone; 4604 4605 ierr = KSPGetPC(ksp, &pc);CHKERRQ(ierr); 4606 ierr = PetscObjectTypeCompare((PetscObject) pc, PCNONE, &isNone);CHKERRQ(ierr); 4607 ierr = KSPGetPCSide(ksp,&pcside);CHKERRQ(ierr); 4608 if (pcside == PC_RIGHT || isNone) { /* XXX Should we also test KSP_UNPRECONDITIONED_NORM ? */ 4609 /* KSP residual is true linear residual */ 4610 ierr = KSPGetResidualNorm(ksp,&kctx->lresid_last);CHKERRQ(ierr); 4611 } else { 4612 /* KSP residual is preconditioned residual */ 4613 /* compute true linear residual norm */ 4614 ierr = VecDuplicate(b,&lres);CHKERRQ(ierr); 4615 ierr = MatMult(snes->jacobian,x,lres);CHKERRQ(ierr); 4616 ierr = VecAYPX(lres,-1.0,b);CHKERRQ(ierr); 4617 ierr = VecNorm(lres,NORM_2,&kctx->lresid_last);CHKERRQ(ierr); 4618 ierr = VecDestroy(&lres);CHKERRQ(ierr); 4619 } 4620 } 4621 PetscFunctionReturn(0); 4622 } 4623 4624 #undef __FUNCT__ 4625 #define __FUNCT__ "SNESGetKSP" 4626 /*@ 4627 SNESGetKSP - Returns the KSP context for a SNES solver. 4628 4629 Not Collective, but if SNES object is parallel, then KSP object is parallel 4630 4631 Input Parameter: 4632 . snes - the SNES context 4633 4634 Output Parameter: 4635 . ksp - the KSP context 4636 4637 Notes: 4638 The user can then directly manipulate the KSP context to set various 4639 options, etc. Likewise, the user can then extract and manipulate the 4640 PC contexts as well. 4641 4642 Level: beginner 4643 4644 .keywords: SNES, nonlinear, get, KSP, context 4645 4646 .seealso: KSPGetPC(), SNESCreate(), KSPCreate(), SNESSetKSP() 4647 @*/ 4648 PetscErrorCode SNESGetKSP(SNES snes,KSP *ksp) 4649 { 4650 PetscErrorCode ierr; 4651 4652 PetscFunctionBegin; 4653 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4654 PetscValidPointer(ksp,2); 4655 4656 if (!snes->ksp) { 4657 PetscBool monitor = PETSC_FALSE; 4658 4659 ierr = KSPCreate(PetscObjectComm((PetscObject)snes),&snes->ksp);CHKERRQ(ierr); 4660 ierr = PetscObjectIncrementTabLevel((PetscObject)snes->ksp,(PetscObject)snes,1);CHKERRQ(ierr); 4661 ierr = PetscLogObjectParent((PetscObject)snes,(PetscObject)snes->ksp);CHKERRQ(ierr); 4662 4663 ierr = KSPSetPreSolve(snes->ksp,(PetscErrorCode (*)(KSP,Vec,Vec,void*))KSPPreSolve_SNESEW,snes);CHKERRQ(ierr); 4664 ierr = KSPSetPostSolve(snes->ksp,(PetscErrorCode (*)(KSP,Vec,Vec,void*))KSPPostSolve_SNESEW,snes);CHKERRQ(ierr); 4665 4666 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-ksp_monitor_snes",&monitor,NULL);CHKERRQ(ierr); 4667 if (monitor) { 4668 ierr = KSPMonitorSet(snes->ksp,KSPMonitorSNES,snes,NULL);CHKERRQ(ierr); 4669 } 4670 monitor = PETSC_FALSE; 4671 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-ksp_monitor_snes_lg",&monitor,NULL);CHKERRQ(ierr); 4672 if (monitor) { 4673 PetscObject *objs; 4674 ierr = KSPMonitorSNESLGResidualNormCreate(0,0,PETSC_DECIDE,PETSC_DECIDE,600,600,&objs);CHKERRQ(ierr); 4675 objs[0] = (PetscObject) snes; 4676 ierr = KSPMonitorSet(snes->ksp,(PetscErrorCode (*)(KSP,PetscInt,PetscReal,void*))KSPMonitorSNESLGResidualNorm,objs,(PetscErrorCode (*)(void**))KSPMonitorSNESLGResidualNormDestroy);CHKERRQ(ierr); 4677 } 4678 } 4679 *ksp = snes->ksp; 4680 PetscFunctionReturn(0); 4681 } 4682 4683 4684 #include <petsc-private/dmimpl.h> 4685 #undef __FUNCT__ 4686 #define __FUNCT__ "SNESSetDM" 4687 /*@ 4688 SNESSetDM - Sets the DM that may be used by some preconditioners 4689 4690 Logically Collective on SNES 4691 4692 Input Parameters: 4693 + snes - the preconditioner context 4694 - dm - the dm 4695 4696 Level: intermediate 4697 4698 .seealso: SNESGetDM(), KSPSetDM(), KSPGetDM() 4699 @*/ 4700 PetscErrorCode SNESSetDM(SNES snes,DM dm) 4701 { 4702 PetscErrorCode ierr; 4703 KSP ksp; 4704 DMSNES sdm; 4705 4706 PetscFunctionBegin; 4707 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4708 if (dm) {ierr = PetscObjectReference((PetscObject)dm);CHKERRQ(ierr);} 4709 if (snes->dm) { /* Move the DMSNES context over to the new DM unless the new DM already has one */ 4710 if (snes->dm->dmsnes && snes->dmAuto && !dm->dmsnes) { 4711 ierr = DMCopyDMSNES(snes->dm,dm);CHKERRQ(ierr); 4712 ierr = DMGetDMSNES(snes->dm,&sdm);CHKERRQ(ierr); 4713 if (sdm->originaldm == snes->dm) sdm->originaldm = dm; /* Grant write privileges to the replacement DM */ 4714 } 4715 ierr = DMDestroy(&snes->dm);CHKERRQ(ierr); 4716 } 4717 snes->dm = dm; 4718 snes->dmAuto = PETSC_FALSE; 4719 4720 ierr = SNESGetKSP(snes,&ksp);CHKERRQ(ierr); 4721 ierr = KSPSetDM(ksp,dm);CHKERRQ(ierr); 4722 ierr = KSPSetDMActive(ksp,PETSC_FALSE);CHKERRQ(ierr); 4723 if (snes->pc) { 4724 ierr = SNESSetDM(snes->pc, snes->dm);CHKERRQ(ierr); 4725 ierr = SNESSetNPCSide(snes,snes->pcside);CHKERRQ(ierr); 4726 } 4727 PetscFunctionReturn(0); 4728 } 4729 4730 #undef __FUNCT__ 4731 #define __FUNCT__ "SNESGetDM" 4732 /*@ 4733 SNESGetDM - Gets the DM that may be used by some preconditioners 4734 4735 Not Collective but DM obtained is parallel on SNES 4736 4737 Input Parameter: 4738 . snes - the preconditioner context 4739 4740 Output Parameter: 4741 . dm - the dm 4742 4743 Level: intermediate 4744 4745 .seealso: SNESSetDM(), KSPSetDM(), KSPGetDM() 4746 @*/ 4747 PetscErrorCode SNESGetDM(SNES snes,DM *dm) 4748 { 4749 PetscErrorCode ierr; 4750 4751 PetscFunctionBegin; 4752 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4753 if (!snes->dm) { 4754 ierr = DMShellCreate(PetscObjectComm((PetscObject)snes),&snes->dm);CHKERRQ(ierr); 4755 snes->dmAuto = PETSC_TRUE; 4756 } 4757 *dm = snes->dm; 4758 PetscFunctionReturn(0); 4759 } 4760 4761 #undef __FUNCT__ 4762 #define __FUNCT__ "SNESSetNPC" 4763 /*@ 4764 SNESSetNPC - Sets the nonlinear preconditioner to be used. 4765 4766 Collective on SNES 4767 4768 Input Parameters: 4769 + snes - iterative context obtained from SNESCreate() 4770 - pc - the preconditioner object 4771 4772 Notes: 4773 Use SNESGetNPC() to retrieve the preconditioner context (for example, 4774 to configure it using the API). 4775 4776 Level: developer 4777 4778 .keywords: SNES, set, precondition 4779 .seealso: SNESGetNPC(), SNESHasNPC() 4780 @*/ 4781 PetscErrorCode SNESSetNPC(SNES snes, SNES pc) 4782 { 4783 PetscErrorCode ierr; 4784 4785 PetscFunctionBegin; 4786 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); 4787 PetscValidHeaderSpecific(pc, SNES_CLASSID, 2); 4788 PetscCheckSameComm(snes, 1, pc, 2); 4789 ierr = PetscObjectReference((PetscObject) pc);CHKERRQ(ierr); 4790 ierr = SNESDestroy(&snes->pc);CHKERRQ(ierr); 4791 snes->pc = pc; 4792 ierr = PetscLogObjectParent((PetscObject)snes, (PetscObject)snes->pc);CHKERRQ(ierr); 4793 PetscFunctionReturn(0); 4794 } 4795 4796 #undef __FUNCT__ 4797 #define __FUNCT__ "SNESGetNPC" 4798 /*@ 4799 SNESGetNPC - Creates a nonlinear preconditioning solver (SNES) to be used to precondition the nonlinear solver. 4800 4801 Not Collective 4802 4803 Input Parameter: 4804 . snes - iterative context obtained from SNESCreate() 4805 4806 Output Parameter: 4807 . pc - preconditioner context 4808 4809 Notes: If a SNES was previously set with SNESSetNPC() then that SNES is returned. 4810 4811 Level: developer 4812 4813 .keywords: SNES, get, preconditioner 4814 .seealso: SNESSetNPC(), SNESHasNPC() 4815 @*/ 4816 PetscErrorCode SNESGetNPC(SNES snes, SNES *pc) 4817 { 4818 PetscErrorCode ierr; 4819 const char *optionsprefix; 4820 4821 PetscFunctionBegin; 4822 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); 4823 PetscValidPointer(pc, 2); 4824 if (!snes->pc) { 4825 ierr = SNESCreate(PetscObjectComm((PetscObject)snes),&snes->pc);CHKERRQ(ierr); 4826 ierr = PetscObjectIncrementTabLevel((PetscObject)snes->pc,(PetscObject)snes,1);CHKERRQ(ierr); 4827 ierr = PetscLogObjectParent((PetscObject)snes,(PetscObject)snes->pc);CHKERRQ(ierr); 4828 ierr = SNESGetOptionsPrefix(snes,&optionsprefix);CHKERRQ(ierr); 4829 ierr = SNESSetOptionsPrefix(snes->pc,optionsprefix);CHKERRQ(ierr); 4830 ierr = SNESAppendOptionsPrefix(snes->pc,"npc_");CHKERRQ(ierr); 4831 ierr = SNESSetCountersReset(snes->pc,PETSC_FALSE);CHKERRQ(ierr); 4832 } 4833 *pc = snes->pc; 4834 PetscFunctionReturn(0); 4835 } 4836 4837 #undef __FUNCT__ 4838 #define __FUNCT__ "SNESHasNPC" 4839 /*@ 4840 SNESHasNPC - Returns whether a nonlinear preconditioner exists 4841 4842 Not Collective 4843 4844 Input Parameter: 4845 . snes - iterative context obtained from SNESCreate() 4846 4847 Output Parameter: 4848 . has_npc - whether the SNES has an NPC or not 4849 4850 Level: developer 4851 4852 .keywords: SNES, has, preconditioner 4853 .seealso: SNESSetNPC(), SNESGetNPC() 4854 @*/ 4855 PetscErrorCode SNESHasNPC(SNES snes, PetscBool *has_npc) 4856 { 4857 PetscFunctionBegin; 4858 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); 4859 *has_npc = (PetscBool) (snes->pc != NULL); 4860 PetscFunctionReturn(0); 4861 } 4862 4863 #undef __FUNCT__ 4864 #define __FUNCT__ "SNESSetNPCSide" 4865 /*@ 4866 SNESSetNPCSide - Sets the preconditioning side. 4867 4868 Logically Collective on SNES 4869 4870 Input Parameter: 4871 . snes - iterative context obtained from SNESCreate() 4872 4873 Output Parameter: 4874 . side - the preconditioning side, where side is one of 4875 .vb 4876 PC_LEFT - left preconditioning (default) 4877 PC_RIGHT - right preconditioning 4878 .ve 4879 4880 Options Database Keys: 4881 . -snes_pc_side <right,left> 4882 4883 Level: intermediate 4884 4885 .keywords: SNES, set, right, left, side, preconditioner, flag 4886 4887 .seealso: SNESGetNPCSide(), KSPSetPCSide() 4888 @*/ 4889 PetscErrorCode SNESSetNPCSide(SNES snes,PCSide side) 4890 { 4891 PetscFunctionBegin; 4892 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4893 PetscValidLogicalCollectiveEnum(snes,side,2); 4894 snes->pcside = side; 4895 PetscFunctionReturn(0); 4896 } 4897 4898 #undef __FUNCT__ 4899 #define __FUNCT__ "SNESGetNPCSide" 4900 /*@ 4901 SNESGetNPCSide - Gets the preconditioning side. 4902 4903 Not Collective 4904 4905 Input Parameter: 4906 . snes - iterative context obtained from SNESCreate() 4907 4908 Output Parameter: 4909 . side - the preconditioning side, where side is one of 4910 .vb 4911 PC_LEFT - left preconditioning (default) 4912 PC_RIGHT - right preconditioning 4913 .ve 4914 4915 Level: intermediate 4916 4917 .keywords: SNES, get, right, left, side, preconditioner, flag 4918 4919 .seealso: SNESSetNPCSide(), KSPGetPCSide() 4920 @*/ 4921 PetscErrorCode SNESGetNPCSide(SNES snes,PCSide *side) 4922 { 4923 PetscFunctionBegin; 4924 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4925 PetscValidPointer(side,2); 4926 *side = snes->pcside; 4927 PetscFunctionReturn(0); 4928 } 4929 4930 #undef __FUNCT__ 4931 #define __FUNCT__ "SNESSetLineSearch" 4932 /*@ 4933 SNESSetLineSearch - Sets the linesearch on the SNES instance. 4934 4935 Collective on SNES 4936 4937 Input Parameters: 4938 + snes - iterative context obtained from SNESCreate() 4939 - linesearch - the linesearch object 4940 4941 Notes: 4942 Use SNESGetLineSearch() to retrieve the preconditioner context (for example, 4943 to configure it using the API). 4944 4945 Level: developer 4946 4947 .keywords: SNES, set, linesearch 4948 .seealso: SNESGetLineSearch() 4949 @*/ 4950 PetscErrorCode SNESSetLineSearch(SNES snes, SNESLineSearch linesearch) 4951 { 4952 PetscErrorCode ierr; 4953 4954 PetscFunctionBegin; 4955 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); 4956 PetscValidHeaderSpecific(linesearch, SNESLINESEARCH_CLASSID, 2); 4957 PetscCheckSameComm(snes, 1, linesearch, 2); 4958 ierr = PetscObjectReference((PetscObject) linesearch);CHKERRQ(ierr); 4959 ierr = SNESLineSearchDestroy(&snes->linesearch);CHKERRQ(ierr); 4960 4961 snes->linesearch = linesearch; 4962 4963 ierr = PetscLogObjectParent((PetscObject)snes, (PetscObject)snes->linesearch);CHKERRQ(ierr); 4964 PetscFunctionReturn(0); 4965 } 4966 4967 #undef __FUNCT__ 4968 #define __FUNCT__ "SNESGetLineSearch" 4969 /*@ 4970 SNESGetLineSearch - Returns a pointer to the line search context set with SNESSetLineSearch() 4971 or creates a default line search instance associated with the SNES and returns it. 4972 4973 Not Collective 4974 4975 Input Parameter: 4976 . snes - iterative context obtained from SNESCreate() 4977 4978 Output Parameter: 4979 . linesearch - linesearch context 4980 4981 Level: beginner 4982 4983 .keywords: SNES, get, linesearch 4984 .seealso: SNESSetLineSearch(), SNESLineSearchCreate() 4985 @*/ 4986 PetscErrorCode SNESGetLineSearch(SNES snes, SNESLineSearch *linesearch) 4987 { 4988 PetscErrorCode ierr; 4989 const char *optionsprefix; 4990 4991 PetscFunctionBegin; 4992 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); 4993 PetscValidPointer(linesearch, 2); 4994 if (!snes->linesearch) { 4995 ierr = SNESGetOptionsPrefix(snes, &optionsprefix);CHKERRQ(ierr); 4996 ierr = SNESLineSearchCreate(PetscObjectComm((PetscObject)snes), &snes->linesearch);CHKERRQ(ierr); 4997 ierr = SNESLineSearchSetSNES(snes->linesearch, snes);CHKERRQ(ierr); 4998 ierr = SNESLineSearchAppendOptionsPrefix(snes->linesearch, optionsprefix);CHKERRQ(ierr); 4999 ierr = PetscObjectIncrementTabLevel((PetscObject) snes->linesearch, (PetscObject) snes, 1);CHKERRQ(ierr); 5000 ierr = PetscLogObjectParent((PetscObject)snes, (PetscObject)snes->linesearch);CHKERRQ(ierr); 5001 } 5002 *linesearch = snes->linesearch; 5003 PetscFunctionReturn(0); 5004 } 5005 5006 #if defined(PETSC_HAVE_MATLAB_ENGINE) 5007 #include <mex.h> 5008 5009 typedef struct {char *funcname; mxArray *ctx;} SNESMatlabContext; 5010 5011 #undef __FUNCT__ 5012 #define __FUNCT__ "SNESComputeFunction_Matlab" 5013 /* 5014 SNESComputeFunction_Matlab - Calls the function that has been set with SNESSetFunctionMatlab(). 5015 5016 Collective on SNES 5017 5018 Input Parameters: 5019 + snes - the SNES context 5020 - x - input vector 5021 5022 Output Parameter: 5023 . y - function vector, as set by SNESSetFunction() 5024 5025 Notes: 5026 SNESComputeFunction() is typically used within nonlinear solvers 5027 implementations, so most users would not generally call this routine 5028 themselves. 5029 5030 Level: developer 5031 5032 .keywords: SNES, nonlinear, compute, function 5033 5034 .seealso: SNESSetFunction(), SNESGetFunction() 5035 */ 5036 PetscErrorCode SNESComputeFunction_Matlab(SNES snes,Vec x,Vec y, void *ctx) 5037 { 5038 PetscErrorCode ierr; 5039 SNESMatlabContext *sctx = (SNESMatlabContext*)ctx; 5040 int nlhs = 1,nrhs = 5; 5041 mxArray *plhs[1],*prhs[5]; 5042 long long int lx = 0,ly = 0,ls = 0; 5043 5044 PetscFunctionBegin; 5045 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 5046 PetscValidHeaderSpecific(x,VEC_CLASSID,2); 5047 PetscValidHeaderSpecific(y,VEC_CLASSID,3); 5048 PetscCheckSameComm(snes,1,x,2); 5049 PetscCheckSameComm(snes,1,y,3); 5050 5051 /* call Matlab function in ctx with arguments x and y */ 5052 5053 ierr = PetscMemcpy(&ls,&snes,sizeof(snes));CHKERRQ(ierr); 5054 ierr = PetscMemcpy(&lx,&x,sizeof(x));CHKERRQ(ierr); 5055 ierr = PetscMemcpy(&ly,&y,sizeof(x));CHKERRQ(ierr); 5056 prhs[0] = mxCreateDoubleScalar((double)ls); 5057 prhs[1] = mxCreateDoubleScalar((double)lx); 5058 prhs[2] = mxCreateDoubleScalar((double)ly); 5059 prhs[3] = mxCreateString(sctx->funcname); 5060 prhs[4] = sctx->ctx; 5061 ierr = mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscSNESComputeFunctionInternal");CHKERRQ(ierr); 5062 ierr = mxGetScalar(plhs[0]);CHKERRQ(ierr); 5063 mxDestroyArray(prhs[0]); 5064 mxDestroyArray(prhs[1]); 5065 mxDestroyArray(prhs[2]); 5066 mxDestroyArray(prhs[3]); 5067 mxDestroyArray(plhs[0]); 5068 PetscFunctionReturn(0); 5069 } 5070 5071 #undef __FUNCT__ 5072 #define __FUNCT__ "SNESSetFunctionMatlab" 5073 /* 5074 SNESSetFunctionMatlab - Sets the function evaluation routine and function 5075 vector for use by the SNES routines in solving systems of nonlinear 5076 equations from MATLAB. Here the function is a string containing the name of a MATLAB function 5077 5078 Logically Collective on SNES 5079 5080 Input Parameters: 5081 + snes - the SNES context 5082 . r - vector to store function value 5083 - f - function evaluation routine 5084 5085 Notes: 5086 The Newton-like methods typically solve linear systems of the form 5087 $ f'(x) x = -f(x), 5088 where f'(x) denotes the Jacobian matrix and f(x) is the function. 5089 5090 Level: beginner 5091 5092 Developer Note: This bleeds the allocated memory SNESMatlabContext *sctx; 5093 5094 .keywords: SNES, nonlinear, set, function 5095 5096 .seealso: SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESSetFunction() 5097 */ 5098 PetscErrorCode SNESSetFunctionMatlab(SNES snes,Vec r,const char *f,mxArray *ctx) 5099 { 5100 PetscErrorCode ierr; 5101 SNESMatlabContext *sctx; 5102 5103 PetscFunctionBegin; 5104 /* currently sctx is memory bleed */ 5105 ierr = PetscNew(&sctx);CHKERRQ(ierr); 5106 ierr = PetscStrallocpy(f,&sctx->funcname);CHKERRQ(ierr); 5107 /* 5108 This should work, but it doesn't 5109 sctx->ctx = ctx; 5110 mexMakeArrayPersistent(sctx->ctx); 5111 */ 5112 sctx->ctx = mxDuplicateArray(ctx); 5113 ierr = SNESSetFunction(snes,r,SNESComputeFunction_Matlab,sctx);CHKERRQ(ierr); 5114 PetscFunctionReturn(0); 5115 } 5116 5117 #undef __FUNCT__ 5118 #define __FUNCT__ "SNESComputeJacobian_Matlab" 5119 /* 5120 SNESComputeJacobian_Matlab - Calls the function that has been set with SNESSetJacobianMatlab(). 5121 5122 Collective on SNES 5123 5124 Input Parameters: 5125 + snes - the SNES context 5126 . x - input vector 5127 . A, B - the matrices 5128 - ctx - user context 5129 5130 Level: developer 5131 5132 .keywords: SNES, nonlinear, compute, function 5133 5134 .seealso: SNESSetFunction(), SNESGetFunction() 5135 @*/ 5136 PetscErrorCode SNESComputeJacobian_Matlab(SNES snes,Vec x,Mat A,Mat B,void *ctx) 5137 { 5138 PetscErrorCode ierr; 5139 SNESMatlabContext *sctx = (SNESMatlabContext*)ctx; 5140 int nlhs = 2,nrhs = 6; 5141 mxArray *plhs[2],*prhs[6]; 5142 long long int lx = 0,lA = 0,ls = 0, lB = 0; 5143 5144 PetscFunctionBegin; 5145 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 5146 PetscValidHeaderSpecific(x,VEC_CLASSID,2); 5147 5148 /* call Matlab function in ctx with arguments x and y */ 5149 5150 ierr = PetscMemcpy(&ls,&snes,sizeof(snes));CHKERRQ(ierr); 5151 ierr = PetscMemcpy(&lx,&x,sizeof(x));CHKERRQ(ierr); 5152 ierr = PetscMemcpy(&lA,A,sizeof(x));CHKERRQ(ierr); 5153 ierr = PetscMemcpy(&lB,B,sizeof(x));CHKERRQ(ierr); 5154 prhs[0] = mxCreateDoubleScalar((double)ls); 5155 prhs[1] = mxCreateDoubleScalar((double)lx); 5156 prhs[2] = mxCreateDoubleScalar((double)lA); 5157 prhs[3] = mxCreateDoubleScalar((double)lB); 5158 prhs[4] = mxCreateString(sctx->funcname); 5159 prhs[5] = sctx->ctx; 5160 ierr = mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscSNESComputeJacobianInternal");CHKERRQ(ierr); 5161 ierr = mxGetScalar(plhs[0]);CHKERRQ(ierr); 5162 mxDestroyArray(prhs[0]); 5163 mxDestroyArray(prhs[1]); 5164 mxDestroyArray(prhs[2]); 5165 mxDestroyArray(prhs[3]); 5166 mxDestroyArray(prhs[4]); 5167 mxDestroyArray(plhs[0]); 5168 mxDestroyArray(plhs[1]); 5169 PetscFunctionReturn(0); 5170 } 5171 5172 #undef __FUNCT__ 5173 #define __FUNCT__ "SNESSetJacobianMatlab" 5174 /* 5175 SNESSetJacobianMatlab - Sets the Jacobian function evaluation routine and two empty Jacobian matrices 5176 vector for use by the SNES routines in solving systems of nonlinear 5177 equations from MATLAB. Here the function is a string containing the name of a MATLAB function 5178 5179 Logically Collective on SNES 5180 5181 Input Parameters: 5182 + snes - the SNES context 5183 . A,B - Jacobian matrices 5184 . J - function evaluation routine 5185 - ctx - user context 5186 5187 Level: developer 5188 5189 Developer Note: This bleeds the allocated memory SNESMatlabContext *sctx; 5190 5191 .keywords: SNES, nonlinear, set, function 5192 5193 .seealso: SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESSetFunction(), J 5194 */ 5195 PetscErrorCode SNESSetJacobianMatlab(SNES snes,Mat A,Mat B,const char *J,mxArray *ctx) 5196 { 5197 PetscErrorCode ierr; 5198 SNESMatlabContext *sctx; 5199 5200 PetscFunctionBegin; 5201 /* currently sctx is memory bleed */ 5202 ierr = PetscNew(&sctx);CHKERRQ(ierr); 5203 ierr = PetscStrallocpy(J,&sctx->funcname);CHKERRQ(ierr); 5204 /* 5205 This should work, but it doesn't 5206 sctx->ctx = ctx; 5207 mexMakeArrayPersistent(sctx->ctx); 5208 */ 5209 sctx->ctx = mxDuplicateArray(ctx); 5210 ierr = SNESSetJacobian(snes,A,B,SNESComputeJacobian_Matlab,sctx);CHKERRQ(ierr); 5211 PetscFunctionReturn(0); 5212 } 5213 5214 #undef __FUNCT__ 5215 #define __FUNCT__ "SNESMonitor_Matlab" 5216 /* 5217 SNESMonitor_Matlab - Calls the function that has been set with SNESMonitorSetMatlab(). 5218 5219 Collective on SNES 5220 5221 .seealso: SNESSetFunction(), SNESGetFunction() 5222 @*/ 5223 PetscErrorCode SNESMonitor_Matlab(SNES snes,PetscInt it, PetscReal fnorm, void *ctx) 5224 { 5225 PetscErrorCode ierr; 5226 SNESMatlabContext *sctx = (SNESMatlabContext*)ctx; 5227 int nlhs = 1,nrhs = 6; 5228 mxArray *plhs[1],*prhs[6]; 5229 long long int lx = 0,ls = 0; 5230 Vec x = snes->vec_sol; 5231 5232 PetscFunctionBegin; 5233 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 5234 5235 ierr = PetscMemcpy(&ls,&snes,sizeof(snes));CHKERRQ(ierr); 5236 ierr = PetscMemcpy(&lx,&x,sizeof(x));CHKERRQ(ierr); 5237 prhs[0] = mxCreateDoubleScalar((double)ls); 5238 prhs[1] = mxCreateDoubleScalar((double)it); 5239 prhs[2] = mxCreateDoubleScalar((double)fnorm); 5240 prhs[3] = mxCreateDoubleScalar((double)lx); 5241 prhs[4] = mxCreateString(sctx->funcname); 5242 prhs[5] = sctx->ctx; 5243 ierr = mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscSNESMonitorInternal");CHKERRQ(ierr); 5244 ierr = mxGetScalar(plhs[0]);CHKERRQ(ierr); 5245 mxDestroyArray(prhs[0]); 5246 mxDestroyArray(prhs[1]); 5247 mxDestroyArray(prhs[2]); 5248 mxDestroyArray(prhs[3]); 5249 mxDestroyArray(prhs[4]); 5250 mxDestroyArray(plhs[0]); 5251 PetscFunctionReturn(0); 5252 } 5253 5254 #undef __FUNCT__ 5255 #define __FUNCT__ "SNESMonitorSetMatlab" 5256 /* 5257 SNESMonitorSetMatlab - Sets the monitor function from MATLAB 5258 5259 Level: developer 5260 5261 Developer Note: This bleeds the allocated memory SNESMatlabContext *sctx; 5262 5263 .keywords: SNES, nonlinear, set, function 5264 5265 .seealso: SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESSetFunction() 5266 */ 5267 PetscErrorCode SNESMonitorSetMatlab(SNES snes,const char *f,mxArray *ctx) 5268 { 5269 PetscErrorCode ierr; 5270 SNESMatlabContext *sctx; 5271 5272 PetscFunctionBegin; 5273 /* currently sctx is memory bleed */ 5274 ierr = PetscNew(&sctx);CHKERRQ(ierr); 5275 ierr = PetscStrallocpy(f,&sctx->funcname);CHKERRQ(ierr); 5276 /* 5277 This should work, but it doesn't 5278 sctx->ctx = ctx; 5279 mexMakeArrayPersistent(sctx->ctx); 5280 */ 5281 sctx->ctx = mxDuplicateArray(ctx); 5282 ierr = SNESMonitorSet(snes,SNESMonitor_Matlab,sctx,NULL);CHKERRQ(ierr); 5283 PetscFunctionReturn(0); 5284 } 5285 5286 #endif 5287