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