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