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