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