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