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 /* Create context to compute Eisenstat-Walker relative tolerance for KSP */ 1510 ierr = PetscNewLog(snes,&kctx);CHKERRQ(ierr); 1511 1512 snes->kspconvctx = (void*)kctx; 1513 kctx->version = 2; 1514 kctx->rtol_0 = .3; /* Eisenstat and Walker suggest rtol_0=.5, but 1515 this was too large for some test cases */ 1516 kctx->rtol_last = 0.0; 1517 kctx->rtol_max = .9; 1518 kctx->gamma = 1.0; 1519 kctx->alpha = .5*(1.0 + PetscSqrtReal(5.0)); 1520 kctx->alpha2 = kctx->alpha; 1521 kctx->threshold = .1; 1522 kctx->lresid_last = 0.0; 1523 kctx->norm_last = 0.0; 1524 1525 *outsnes = snes; 1526 PetscFunctionReturn(0); 1527 } 1528 1529 /*MC 1530 SNESFunction - Functional form used to convey the nonlinear function to be solved by SNES 1531 1532 Synopsis: 1533 #include "petscsnes.h" 1534 PetscErrorCode SNESFunction(SNES snes,Vec x,Vec f,void *ctx); 1535 1536 Input Parameters: 1537 + snes - the SNES context 1538 . x - state at which to evaluate residual 1539 - ctx - optional user-defined function context, passed in with SNESSetFunction() 1540 1541 Output Parameter: 1542 . f - vector to put residual (function value) 1543 1544 Level: intermediate 1545 1546 .seealso: SNESSetFunction(), SNESGetFunction() 1547 M*/ 1548 1549 #undef __FUNCT__ 1550 #define __FUNCT__ "SNESSetFunction" 1551 /*@C 1552 SNESSetFunction - Sets the function evaluation routine and function 1553 vector for use by the SNES routines in solving systems of nonlinear 1554 equations. 1555 1556 Logically Collective on SNES 1557 1558 Input Parameters: 1559 + snes - the SNES context 1560 . r - vector to store function value 1561 . f - function evaluation routine; see SNESFunction for calling sequence details 1562 - ctx - [optional] user-defined context for private data for the 1563 function evaluation routine (may be NULL) 1564 1565 Notes: 1566 The Newton-like methods typically solve linear systems of the form 1567 $ f'(x) x = -f(x), 1568 where f'(x) denotes the Jacobian matrix and f(x) is the function. 1569 1570 Level: beginner 1571 1572 .keywords: SNES, nonlinear, set, function 1573 1574 .seealso: SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESSetPicard(), SNESFunction 1575 @*/ 1576 PetscErrorCode SNESSetFunction(SNES snes,Vec r,PetscErrorCode (*f)(SNES,Vec,Vec,void*),void *ctx) 1577 { 1578 PetscErrorCode ierr; 1579 DM dm; 1580 1581 PetscFunctionBegin; 1582 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1583 if (r) { 1584 PetscValidHeaderSpecific(r,VEC_CLASSID,2); 1585 PetscCheckSameComm(snes,1,r,2); 1586 ierr = PetscObjectReference((PetscObject)r);CHKERRQ(ierr); 1587 ierr = VecDestroy(&snes->vec_func);CHKERRQ(ierr); 1588 1589 snes->vec_func = r; 1590 } 1591 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 1592 ierr = DMSNESSetFunction(dm,f,ctx);CHKERRQ(ierr); 1593 PetscFunctionReturn(0); 1594 } 1595 1596 1597 #undef __FUNCT__ 1598 #define __FUNCT__ "SNESSetInitialFunction" 1599 /*@C 1600 SNESSetInitialFunction - Sets the function vector to be used as the 1601 function norm at the initialization of the method. In some 1602 instances, the user has precomputed the function before calling 1603 SNESSolve. This function allows one to avoid a redundant call 1604 to SNESComputeFunction in that case. 1605 1606 Logically Collective on SNES 1607 1608 Input Parameters: 1609 + snes - the SNES context 1610 - f - vector to store function value 1611 1612 Notes: 1613 This should not be modified during the solution procedure. 1614 1615 This is used extensively in the SNESFAS hierarchy and in nonlinear preconditioning. 1616 1617 Level: developer 1618 1619 .keywords: SNES, nonlinear, set, function 1620 1621 .seealso: SNESSetFunction(), SNESComputeFunction(), SNESSetInitialFunctionNorm() 1622 @*/ 1623 PetscErrorCode SNESSetInitialFunction(SNES snes, Vec f) 1624 { 1625 PetscErrorCode ierr; 1626 Vec vec_func; 1627 1628 PetscFunctionBegin; 1629 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1630 PetscValidHeaderSpecific(f,VEC_CLASSID,2); 1631 PetscCheckSameComm(snes,1,f,2); 1632 if (snes->pcside == PC_LEFT && snes->functype == SNES_FUNCTION_PRECONDITIONED) { 1633 snes->vec_func_init_set = PETSC_FALSE; 1634 PetscFunctionReturn(0); 1635 } 1636 ierr = SNESGetFunction(snes,&vec_func,NULL,NULL);CHKERRQ(ierr); 1637 ierr = VecCopy(f, vec_func);CHKERRQ(ierr); 1638 1639 snes->vec_func_init_set = PETSC_TRUE; 1640 PetscFunctionReturn(0); 1641 } 1642 1643 #undef __FUNCT__ 1644 #define __FUNCT__ "SNESSetNormSchedule" 1645 /*@ 1646 SNESSetNormSchedule - Sets the SNESNormSchedule used in covergence and monitoring 1647 of the SNES method. 1648 1649 Logically Collective on SNES 1650 1651 Input Parameters: 1652 + snes - the SNES context 1653 - normschedule - the frequency of norm computation 1654 1655 Options Database Key: 1656 . -snes_norm_schedule <none, always, initialonly, finalonly, initalfinalonly> 1657 1658 Notes: 1659 Only certain SNES methods support certain SNESNormSchedules. Most require evaluation 1660 of the nonlinear function and the taking of its norm at every iteration to 1661 even ensure convergence at all. However, methods such as custom Gauss-Seidel methods 1662 (SNESNGS) and the like do not require the norm of the function to be computed, and therfore 1663 may either be monitored for convergence or not. As these are often used as nonlinear 1664 preconditioners, monitoring the norm of their error is not a useful enterprise within 1665 their solution. 1666 1667 Level: developer 1668 1669 .keywords: SNES, nonlinear, set, function, norm, type 1670 1671 .seealso: SNESGetNormSchedule(), SNESComputeFunction(), VecNorm(), SNESSetFunction(), SNESSetInitialFunction(), SNESNormSchedule 1672 @*/ 1673 PetscErrorCode SNESSetNormSchedule(SNES snes, SNESNormSchedule normschedule) 1674 { 1675 PetscFunctionBegin; 1676 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1677 snes->normschedule = normschedule; 1678 PetscFunctionReturn(0); 1679 } 1680 1681 1682 #undef __FUNCT__ 1683 #define __FUNCT__ "SNESGetNormSchedule" 1684 /*@ 1685 SNESGetNormSchedule - Gets the SNESNormSchedule used in covergence and monitoring 1686 of the SNES method. 1687 1688 Logically Collective on SNES 1689 1690 Input Parameters: 1691 + snes - the SNES context 1692 - normschedule - the type of the norm used 1693 1694 Level: advanced 1695 1696 .keywords: SNES, nonlinear, set, function, norm, type 1697 1698 .seealso: SNESSetNormSchedule(), SNESComputeFunction(), VecNorm(), SNESSetFunction(), SNESSetInitialFunction(), SNESNormSchedule 1699 @*/ 1700 PetscErrorCode SNESGetNormSchedule(SNES snes, SNESNormSchedule *normschedule) 1701 { 1702 PetscFunctionBegin; 1703 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1704 *normschedule = snes->normschedule; 1705 PetscFunctionReturn(0); 1706 } 1707 1708 1709 #undef __FUNCT__ 1710 #define __FUNCT__ "SNESSetFunctionType" 1711 /*@C 1712 SNESSetFunctionType - Sets the SNESNormSchedule used in covergence and monitoring 1713 of the SNES method. 1714 1715 Logically Collective on SNES 1716 1717 Input Parameters: 1718 + snes - the SNES context 1719 - normschedule - the frequency of norm computation 1720 1721 Notes: 1722 Only certain SNES methods support certain SNESNormSchedules. Most require evaluation 1723 of the nonlinear function and the taking of its norm at every iteration to 1724 even ensure convergence at all. However, methods such as custom Gauss-Seidel methods 1725 (SNESNGS) and the like do not require the norm of the function to be computed, and therfore 1726 may either be monitored for convergence or not. As these are often used as nonlinear 1727 preconditioners, monitoring the norm of their error is not a useful enterprise within 1728 their solution. 1729 1730 Level: developer 1731 1732 .keywords: SNES, nonlinear, set, function, norm, type 1733 1734 .seealso: SNESGetNormSchedule(), SNESComputeFunction(), VecNorm(), SNESSetFunction(), SNESSetInitialFunction(), SNESNormSchedule 1735 @*/ 1736 PetscErrorCode SNESSetFunctionType(SNES snes, SNESFunctionType type) 1737 { 1738 PetscFunctionBegin; 1739 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1740 snes->functype = type; 1741 PetscFunctionReturn(0); 1742 } 1743 1744 1745 #undef __FUNCT__ 1746 #define __FUNCT__ "SNESGetFunctionType" 1747 /*@C 1748 SNESGetFunctionType - Gets the SNESNormSchedule used in covergence and monitoring 1749 of the SNES method. 1750 1751 Logically Collective on SNES 1752 1753 Input Parameters: 1754 + snes - the SNES context 1755 - normschedule - the type of the norm used 1756 1757 Level: advanced 1758 1759 .keywords: SNES, nonlinear, set, function, norm, type 1760 1761 .seealso: SNESSetNormSchedule(), SNESComputeFunction(), VecNorm(), SNESSetFunction(), SNESSetInitialFunction(), SNESNormSchedule 1762 @*/ 1763 PetscErrorCode SNESGetFunctionType(SNES snes, SNESFunctionType *type) 1764 { 1765 PetscFunctionBegin; 1766 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1767 *type = snes->functype; 1768 PetscFunctionReturn(0); 1769 } 1770 1771 /*MC 1772 SNESNGSFunction - function used to convey a Gauss-Seidel sweep on the nonlinear function 1773 1774 Synopsis: 1775 #include <petscsnes.h> 1776 $ SNESNGSFunction(SNES snes,Vec x,Vec b,void *ctx); 1777 1778 + X - solution vector 1779 . B - RHS vector 1780 - ctx - optional user-defined Gauss-Seidel context 1781 1782 Level: intermediate 1783 1784 .seealso: SNESSetNGS(), SNESGetNGS() 1785 M*/ 1786 1787 #undef __FUNCT__ 1788 #define __FUNCT__ "SNESSetNGS" 1789 /*@C 1790 SNESSetNGS - Sets the user nonlinear Gauss-Seidel routine for 1791 use with composed nonlinear solvers. 1792 1793 Input Parameters: 1794 + snes - the SNES context 1795 . f - function evaluation routine to apply Gauss-Seidel see SNESNGSFunction 1796 - ctx - [optional] user-defined context for private data for the 1797 smoother evaluation routine (may be NULL) 1798 1799 Notes: 1800 The NGS routines are used by the composed nonlinear solver to generate 1801 a problem appropriate update to the solution, particularly FAS. 1802 1803 Level: intermediate 1804 1805 .keywords: SNES, nonlinear, set, Gauss-Seidel 1806 1807 .seealso: SNESGetFunction(), SNESComputeNGS() 1808 @*/ 1809 PetscErrorCode SNESSetNGS(SNES snes,PetscErrorCode (*f)(SNES,Vec,Vec,void*),void *ctx) 1810 { 1811 PetscErrorCode ierr; 1812 DM dm; 1813 1814 PetscFunctionBegin; 1815 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1816 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 1817 ierr = DMSNESSetNGS(dm,f,ctx);CHKERRQ(ierr); 1818 PetscFunctionReturn(0); 1819 } 1820 1821 #undef __FUNCT__ 1822 #define __FUNCT__ "SNESPicardComputeFunction" 1823 PETSC_EXTERN PetscErrorCode SNESPicardComputeFunction(SNES snes,Vec x,Vec f,void *ctx) 1824 { 1825 PetscErrorCode ierr; 1826 DM dm; 1827 DMSNES sdm; 1828 1829 PetscFunctionBegin; 1830 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 1831 ierr = DMGetDMSNES(dm,&sdm);CHKERRQ(ierr); 1832 /* A(x)*x - b(x) */ 1833 if (sdm->ops->computepfunction) { 1834 ierr = (*sdm->ops->computepfunction)(snes,x,f,sdm->pctx);CHKERRQ(ierr); 1835 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE, "Must call SNESSetPicard() to provide Picard function."); 1836 1837 if (sdm->ops->computepjacobian) { 1838 ierr = (*sdm->ops->computepjacobian)(snes,x,snes->jacobian,snes->jacobian_pre,sdm->pctx);CHKERRQ(ierr); 1839 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE, "Must call SNESSetPicard() to provide Picard matrix."); 1840 ierr = VecScale(f,-1.0);CHKERRQ(ierr); 1841 ierr = MatMultAdd(snes->jacobian,x,f,f);CHKERRQ(ierr); 1842 PetscFunctionReturn(0); 1843 } 1844 1845 #undef __FUNCT__ 1846 #define __FUNCT__ "SNESPicardComputeJacobian" 1847 PETSC_EXTERN PetscErrorCode SNESPicardComputeJacobian(SNES snes,Vec x1,Mat J,Mat B,void *ctx) 1848 { 1849 PetscFunctionBegin; 1850 /* the jacobian matrix should be pre-filled in SNESPicardComputeFunction */ 1851 PetscFunctionReturn(0); 1852 } 1853 1854 #undef __FUNCT__ 1855 #define __FUNCT__ "SNESSetPicard" 1856 /*@C 1857 SNESSetPicard - Use SNES to solve the semilinear-system A(x) x = b(x) via a Picard type iteration (Picard linearization) 1858 1859 Logically Collective on SNES 1860 1861 Input Parameters: 1862 + snes - the SNES context 1863 . r - vector to store function value 1864 . b - function evaluation routine 1865 . Amat - matrix with which A(x) x - b(x) is to be computed 1866 . Pmat - matrix from which preconditioner is computed (usually the same as Amat) 1867 . J - function to compute matrix value, see SNESJacobianFunction for details on its calling sequence 1868 - ctx - [optional] user-defined context for private data for the 1869 function evaluation routine (may be NULL) 1870 1871 Notes: 1872 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 1873 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. 1874 1875 One can call SNESSetPicard() or SNESSetFunction() (and possibly SNESSetJacobian()) but cannot call both 1876 1877 $ 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} 1878 $ Note that when an exact solver is used this corresponds to the "classic" Picard A(x^{n}) x^{n+1} = b(x^{n}) iteration. 1879 1880 Run with -snes_mf_operator to solve the system with Newton's method using A(x^{n}) to construct the preconditioner. 1881 1882 We implement the defect correction form of the Picard iteration because it converges much more generally when inexact linear solvers are used then 1883 the direct Picard iteration A(x^n) x^{n+1} = b(x^n) 1884 1885 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 1886 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 1887 different please contact us at petsc-dev@mcs.anl.gov and we'll have an entirely new argument :-). 1888 1889 Level: intermediate 1890 1891 .keywords: SNES, nonlinear, set, function 1892 1893 .seealso: SNESGetFunction(), SNESSetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESGetPicard(), SNESLineSearchPreCheckPicard(), SNESJacobianFunction 1894 @*/ 1895 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) 1896 { 1897 PetscErrorCode ierr; 1898 DM dm; 1899 1900 PetscFunctionBegin; 1901 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1902 ierr = SNESGetDM(snes, &dm);CHKERRQ(ierr); 1903 ierr = DMSNESSetPicard(dm,b,J,ctx);CHKERRQ(ierr); 1904 ierr = SNESSetFunction(snes,r,SNESPicardComputeFunction,ctx);CHKERRQ(ierr); 1905 ierr = SNESSetJacobian(snes,Amat,Pmat,SNESPicardComputeJacobian,ctx);CHKERRQ(ierr); 1906 PetscFunctionReturn(0); 1907 } 1908 1909 #undef __FUNCT__ 1910 #define __FUNCT__ "SNESGetPicard" 1911 /*@C 1912 SNESGetPicard - Returns the context for the Picard iteration 1913 1914 Not Collective, but Vec is parallel if SNES is parallel. Collective if Vec is requested, but has not been created yet. 1915 1916 Input Parameter: 1917 . snes - the SNES context 1918 1919 Output Parameter: 1920 + r - the function (or NULL) 1921 . f - the function (or NULL); see SNESFunction for calling sequence details 1922 . Amat - the matrix used to defined the operation A(x) x - b(x) (or NULL) 1923 . Pmat - the matrix from which the preconditioner will be constructed (or NULL) 1924 . J - the function for matrix evaluation (or NULL); see SNESJacobianFunction for calling sequence details 1925 - ctx - the function context (or NULL) 1926 1927 Level: advanced 1928 1929 .keywords: SNES, nonlinear, get, function 1930 1931 .seealso: SNESSetPicard(), SNESGetFunction(), SNESGetJacobian(), SNESGetDM(), SNESFunction, SNESJacobianFunction 1932 @*/ 1933 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) 1934 { 1935 PetscErrorCode ierr; 1936 DM dm; 1937 1938 PetscFunctionBegin; 1939 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1940 ierr = SNESGetFunction(snes,r,NULL,NULL);CHKERRQ(ierr); 1941 ierr = SNESGetJacobian(snes,Amat,Pmat,NULL,NULL);CHKERRQ(ierr); 1942 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 1943 ierr = DMSNESGetPicard(dm,f,J,ctx);CHKERRQ(ierr); 1944 PetscFunctionReturn(0); 1945 } 1946 1947 #undef __FUNCT__ 1948 #define __FUNCT__ "SNESSetComputeInitialGuess" 1949 /*@C 1950 SNESSetComputeInitialGuess - Sets a routine used to compute an initial guess for the problem 1951 1952 Logically Collective on SNES 1953 1954 Input Parameters: 1955 + snes - the SNES context 1956 . func - function evaluation routine 1957 - ctx - [optional] user-defined context for private data for the 1958 function evaluation routine (may be NULL) 1959 1960 Calling sequence of func: 1961 $ func (SNES snes,Vec x,void *ctx); 1962 1963 . f - function vector 1964 - ctx - optional user-defined function context 1965 1966 Level: intermediate 1967 1968 .keywords: SNES, nonlinear, set, function 1969 1970 .seealso: SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian() 1971 @*/ 1972 PetscErrorCode SNESSetComputeInitialGuess(SNES snes,PetscErrorCode (*func)(SNES,Vec,void*),void *ctx) 1973 { 1974 PetscFunctionBegin; 1975 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1976 if (func) snes->ops->computeinitialguess = func; 1977 if (ctx) snes->initialguessP = ctx; 1978 PetscFunctionReturn(0); 1979 } 1980 1981 /* --------------------------------------------------------------- */ 1982 #undef __FUNCT__ 1983 #define __FUNCT__ "SNESGetRhs" 1984 /*@C 1985 SNESGetRhs - Gets the vector for solving F(x) = rhs. If rhs is not set 1986 it assumes a zero right hand side. 1987 1988 Logically Collective on SNES 1989 1990 Input Parameter: 1991 . snes - the SNES context 1992 1993 Output Parameter: 1994 . rhs - the right hand side vector or NULL if the right hand side vector is null 1995 1996 Level: intermediate 1997 1998 .keywords: SNES, nonlinear, get, function, right hand side 1999 2000 .seealso: SNESGetSolution(), SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESSetFunction() 2001 @*/ 2002 PetscErrorCode SNESGetRhs(SNES snes,Vec *rhs) 2003 { 2004 PetscFunctionBegin; 2005 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2006 PetscValidPointer(rhs,2); 2007 *rhs = snes->vec_rhs; 2008 PetscFunctionReturn(0); 2009 } 2010 2011 #undef __FUNCT__ 2012 #define __FUNCT__ "SNESComputeFunction" 2013 /*@ 2014 SNESComputeFunction - Calls the function that has been set with SNESSetFunction(). 2015 2016 Collective on SNES 2017 2018 Input Parameters: 2019 + snes - the SNES context 2020 - x - input vector 2021 2022 Output Parameter: 2023 . y - function vector, as set by SNESSetFunction() 2024 2025 Notes: 2026 SNESComputeFunction() is typically used within nonlinear solvers 2027 implementations, so most users would not generally call this routine 2028 themselves. 2029 2030 Level: developer 2031 2032 .keywords: SNES, nonlinear, compute, function 2033 2034 .seealso: SNESSetFunction(), SNESGetFunction() 2035 @*/ 2036 PetscErrorCode SNESComputeFunction(SNES snes,Vec x,Vec y) 2037 { 2038 PetscErrorCode ierr; 2039 DM dm; 2040 DMSNES sdm; 2041 2042 PetscFunctionBegin; 2043 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2044 PetscValidHeaderSpecific(x,VEC_CLASSID,2); 2045 PetscValidHeaderSpecific(y,VEC_CLASSID,3); 2046 PetscCheckSameComm(snes,1,x,2); 2047 PetscCheckSameComm(snes,1,y,3); 2048 ierr = VecValidValues(x,2,PETSC_TRUE);CHKERRQ(ierr); 2049 2050 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 2051 ierr = DMGetDMSNES(dm,&sdm);CHKERRQ(ierr); 2052 if (sdm->ops->computefunction) { 2053 ierr = PetscLogEventBegin(SNES_FunctionEval,snes,x,y,0);CHKERRQ(ierr); 2054 ierr = VecLockPush(x);CHKERRQ(ierr); 2055 PetscStackPush("SNES user function"); 2056 ierr = (*sdm->ops->computefunction)(snes,x,y,sdm->functionctx);CHKERRQ(ierr); 2057 PetscStackPop; 2058 ierr = VecLockPop(x);CHKERRQ(ierr); 2059 ierr = PetscLogEventEnd(SNES_FunctionEval,snes,x,y,0);CHKERRQ(ierr); 2060 } else if (snes->vec_rhs) { 2061 ierr = MatMult(snes->jacobian, x, y);CHKERRQ(ierr); 2062 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE, "Must call SNESSetFunction() or SNESSetDM() before SNESComputeFunction(), likely called from SNESSolve()."); 2063 if (snes->vec_rhs) { 2064 ierr = VecAXPY(y,-1.0,snes->vec_rhs);CHKERRQ(ierr); 2065 } 2066 snes->nfuncs++; 2067 /* 2068 domainerror might not be set on all processes; so we tag vector locally with Inf and the next inner product or norm will 2069 propagate the value to all processes 2070 */ 2071 if (snes->domainerror) { 2072 ierr = VecSetInf(y);CHKERRQ(ierr); 2073 } 2074 PetscFunctionReturn(0); 2075 } 2076 2077 #undef __FUNCT__ 2078 #define __FUNCT__ "SNESComputeNGS" 2079 /*@ 2080 SNESComputeNGS - Calls the Gauss-Seidel function that has been set with SNESSetNGS(). 2081 2082 Collective on SNES 2083 2084 Input Parameters: 2085 + snes - the SNES context 2086 . x - input vector 2087 - b - rhs vector 2088 2089 Output Parameter: 2090 . x - new solution vector 2091 2092 Notes: 2093 SNESComputeNGS() is typically used within composed nonlinear solver 2094 implementations, so most users would not generally call this routine 2095 themselves. 2096 2097 Level: developer 2098 2099 .keywords: SNES, nonlinear, compute, function 2100 2101 .seealso: SNESSetNGS(), SNESComputeFunction() 2102 @*/ 2103 PetscErrorCode SNESComputeNGS(SNES snes,Vec b,Vec x) 2104 { 2105 PetscErrorCode ierr; 2106 DM dm; 2107 DMSNES sdm; 2108 2109 PetscFunctionBegin; 2110 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2111 PetscValidHeaderSpecific(x,VEC_CLASSID,2); 2112 if (b) PetscValidHeaderSpecific(b,VEC_CLASSID,3); 2113 PetscCheckSameComm(snes,1,x,2); 2114 if (b) PetscCheckSameComm(snes,1,b,3); 2115 if (b) {ierr = VecValidValues(b,2,PETSC_TRUE);CHKERRQ(ierr);} 2116 ierr = PetscLogEventBegin(SNES_NGSEval,snes,x,b,0);CHKERRQ(ierr); 2117 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 2118 ierr = DMGetDMSNES(dm,&sdm);CHKERRQ(ierr); 2119 if (sdm->ops->computegs) { 2120 if (b) {ierr = VecLockPush(b);CHKERRQ(ierr);} 2121 PetscStackPush("SNES user NGS"); 2122 ierr = (*sdm->ops->computegs)(snes,x,b,sdm->gsctx);CHKERRQ(ierr); 2123 PetscStackPop; 2124 if (b) {ierr = VecLockPop(b);CHKERRQ(ierr);} 2125 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE, "Must call SNESSetNGS() before SNESComputeNGS(), likely called from SNESSolve()."); 2126 ierr = PetscLogEventEnd(SNES_NGSEval,snes,x,b,0);CHKERRQ(ierr); 2127 PetscFunctionReturn(0); 2128 } 2129 2130 #undef __FUNCT__ 2131 #define __FUNCT__ "SNESComputeJacobian" 2132 /*@ 2133 SNESComputeJacobian - Computes the Jacobian matrix that has been set with SNESSetJacobian(). 2134 2135 Collective on SNES and Mat 2136 2137 Input Parameters: 2138 + snes - the SNES context 2139 - x - input vector 2140 2141 Output Parameters: 2142 + A - Jacobian matrix 2143 - B - optional preconditioning matrix 2144 2145 Options Database Keys: 2146 + -snes_lag_preconditioner <lag> 2147 . -snes_lag_jacobian <lag> 2148 . -snes_compare_explicit - Compare the computed Jacobian to the finite difference Jacobian and output the differences 2149 . -snes_compare_explicit_draw - Compare the computed Jacobian to the finite difference Jacobian and draw the result 2150 . -snes_compare_explicit_contour - Compare the computed Jacobian to the finite difference Jacobian and draw a contour plot with the result 2151 . -snes_compare_operator - Make the comparison options above use the operator instead of the preconditioning matrix 2152 . -snes_compare_coloring - Compute the finite differece Jacobian using coloring and display norms of difference 2153 . -snes_compare_coloring_display - Compute the finite differece Jacobian using coloring and display verbose differences 2154 . -snes_compare_coloring_threshold - Display only those matrix entries that differ by more than a given threshold 2155 . -snes_compare_coloring_threshold_atol - Absolute tolerance for difference in matrix entries to be displayed by -snes_compare_coloring_threshold 2156 . -snes_compare_coloring_threshold_rtol - Relative tolerance for difference in matrix entries to be displayed by -snes_compare_coloring_threshold 2157 . -snes_compare_coloring_draw - Compute the finite differece Jacobian using coloring and draw differences 2158 - -snes_compare_coloring_draw_contour - Compute the finite differece Jacobian using coloring and show contours of matrices and differences 2159 2160 2161 Notes: 2162 Most users should not need to explicitly call this routine, as it 2163 is used internally within the nonlinear solvers. 2164 2165 Level: developer 2166 2167 .keywords: SNES, compute, Jacobian, matrix 2168 2169 .seealso: SNESSetJacobian(), KSPSetOperators(), MatStructure, SNESSetLagPreconditioner(), SNESSetLagJacobian() 2170 @*/ 2171 PetscErrorCode SNESComputeJacobian(SNES snes,Vec X,Mat A,Mat B) 2172 { 2173 PetscErrorCode ierr; 2174 PetscBool flag; 2175 DM dm; 2176 DMSNES sdm; 2177 KSP ksp; 2178 2179 PetscFunctionBegin; 2180 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2181 PetscValidHeaderSpecific(X,VEC_CLASSID,2); 2182 PetscCheckSameComm(snes,1,X,2); 2183 ierr = VecValidValues(X,2,PETSC_TRUE);CHKERRQ(ierr); 2184 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 2185 ierr = DMGetDMSNES(dm,&sdm);CHKERRQ(ierr); 2186 2187 if (!sdm->ops->computejacobian) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_USER,"Must call SNESSetJacobian(), DMSNESSetJacobian(), DMDASNESSetJacobianLocal(), etc"); 2188 2189 /* make sure that MatAssemblyBegin/End() is called on A matrix if it is matrix free */ 2190 2191 if (snes->lagjacobian == -2) { 2192 snes->lagjacobian = -1; 2193 2194 ierr = PetscInfo(snes,"Recomputing Jacobian/preconditioner because lag is -2 (means compute Jacobian, but then never again) \n");CHKERRQ(ierr); 2195 } else if (snes->lagjacobian == -1) { 2196 ierr = PetscInfo(snes,"Reusing Jacobian/preconditioner because lag is -1\n");CHKERRQ(ierr); 2197 ierr = PetscObjectTypeCompare((PetscObject)A,MATMFFD,&flag);CHKERRQ(ierr); 2198 if (flag) { 2199 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2200 ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2201 } 2202 PetscFunctionReturn(0); 2203 } else if (snes->lagjacobian > 1 && (snes->iter + snes->jac_iter) % snes->lagjacobian) { 2204 ierr = PetscInfo2(snes,"Reusing Jacobian/preconditioner because lag is %D and SNES iteration is %D\n",snes->lagjacobian,snes->iter);CHKERRQ(ierr); 2205 ierr = PetscObjectTypeCompare((PetscObject)A,MATMFFD,&flag);CHKERRQ(ierr); 2206 if (flag) { 2207 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2208 ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2209 } 2210 PetscFunctionReturn(0); 2211 } 2212 if (snes->pc && snes->pcside == PC_LEFT) { 2213 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2214 ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2215 PetscFunctionReturn(0); 2216 } 2217 2218 ierr = PetscLogEventBegin(SNES_JacobianEval,snes,X,A,B);CHKERRQ(ierr); 2219 ierr = VecLockPush(X);CHKERRQ(ierr); 2220 PetscStackPush("SNES user Jacobian function"); 2221 ierr = (*sdm->ops->computejacobian)(snes,X,A,B,sdm->jacobianctx);CHKERRQ(ierr); 2222 PetscStackPop; 2223 ierr = VecLockPop(X);CHKERRQ(ierr); 2224 ierr = PetscLogEventEnd(SNES_JacobianEval,snes,X,A,B);CHKERRQ(ierr); 2225 2226 /* the next line ensures that snes->ksp exists */ 2227 ierr = SNESGetKSP(snes,&ksp);CHKERRQ(ierr); 2228 if (snes->lagpreconditioner == -2) { 2229 ierr = PetscInfo(snes,"Rebuilding preconditioner exactly once since lag is -2\n");CHKERRQ(ierr); 2230 ierr = KSPSetReusePreconditioner(snes->ksp,PETSC_FALSE);CHKERRQ(ierr); 2231 snes->lagpreconditioner = -1; 2232 } else if (snes->lagpreconditioner == -1) { 2233 ierr = PetscInfo(snes,"Reusing preconditioner because lag is -1\n");CHKERRQ(ierr); 2234 ierr = KSPSetReusePreconditioner(snes->ksp,PETSC_TRUE);CHKERRQ(ierr); 2235 } else if (snes->lagpreconditioner > 1 && (snes->iter + snes->pre_iter) % snes->lagpreconditioner) { 2236 ierr = PetscInfo2(snes,"Reusing preconditioner because lag is %D and SNES iteration is %D\n",snes->lagpreconditioner,snes->iter);CHKERRQ(ierr); 2237 ierr = KSPSetReusePreconditioner(snes->ksp,PETSC_TRUE);CHKERRQ(ierr); 2238 } else { 2239 ierr = PetscInfo(snes,"Rebuilding preconditioner\n");CHKERRQ(ierr); 2240 ierr = KSPSetReusePreconditioner(snes->ksp,PETSC_FALSE);CHKERRQ(ierr); 2241 } 2242 2243 /* make sure user returned a correct Jacobian and preconditioner */ 2244 /* PetscValidHeaderSpecific(A,MAT_CLASSID,3); 2245 PetscValidHeaderSpecific(B,MAT_CLASSID,4); */ 2246 { 2247 PetscBool flag = PETSC_FALSE,flag_draw = PETSC_FALSE,flag_contour = PETSC_FALSE,flag_operator = PETSC_FALSE; 2248 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_explicit",&flag,NULL);CHKERRQ(ierr); 2249 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_explicit_draw",&flag_draw,NULL);CHKERRQ(ierr); 2250 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_explicit_draw_contour",&flag_contour,NULL);CHKERRQ(ierr); 2251 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_operator",&flag_operator,NULL);CHKERRQ(ierr); 2252 if (flag || flag_draw || flag_contour) { 2253 Mat Bexp_mine = NULL,Bexp,FDexp; 2254 PetscViewer vdraw,vstdout; 2255 PetscBool flg; 2256 if (flag_operator) { 2257 ierr = MatComputeExplicitOperator(A,&Bexp_mine);CHKERRQ(ierr); 2258 Bexp = Bexp_mine; 2259 } else { 2260 /* See if the preconditioning matrix can be viewed and added directly */ 2261 ierr = PetscObjectTypeCompareAny((PetscObject)B,&flg,MATSEQAIJ,MATMPIAIJ,MATSEQDENSE,MATMPIDENSE,MATSEQBAIJ,MATMPIBAIJ,MATSEQSBAIJ,MATMPIBAIJ,"");CHKERRQ(ierr); 2262 if (flg) Bexp = B; 2263 else { 2264 /* If the "preconditioning" matrix is itself MATSHELL or some other type without direct support */ 2265 ierr = MatComputeExplicitOperator(B,&Bexp_mine);CHKERRQ(ierr); 2266 Bexp = Bexp_mine; 2267 } 2268 } 2269 ierr = MatConvert(Bexp,MATSAME,MAT_INITIAL_MATRIX,&FDexp);CHKERRQ(ierr); 2270 ierr = SNESComputeJacobianDefault(snes,X,FDexp,FDexp,NULL);CHKERRQ(ierr); 2271 ierr = PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)snes),&vstdout);CHKERRQ(ierr); 2272 if (flag_draw || flag_contour) { 2273 ierr = PetscViewerDrawOpen(PetscObjectComm((PetscObject)snes),0,"Explicit Jacobians",PETSC_DECIDE,PETSC_DECIDE,300,300,&vdraw);CHKERRQ(ierr); 2274 if (flag_contour) {ierr = PetscViewerPushFormat(vdraw,PETSC_VIEWER_DRAW_CONTOUR);CHKERRQ(ierr);} 2275 } else vdraw = NULL; 2276 ierr = PetscViewerASCIIPrintf(vstdout,"Explicit %s\n",flag_operator ? "Jacobian" : "preconditioning Jacobian");CHKERRQ(ierr); 2277 if (flag) {ierr = MatView(Bexp,vstdout);CHKERRQ(ierr);} 2278 if (vdraw) {ierr = MatView(Bexp,vdraw);CHKERRQ(ierr);} 2279 ierr = PetscViewerASCIIPrintf(vstdout,"Finite difference Jacobian\n");CHKERRQ(ierr); 2280 if (flag) {ierr = MatView(FDexp,vstdout);CHKERRQ(ierr);} 2281 if (vdraw) {ierr = MatView(FDexp,vdraw);CHKERRQ(ierr);} 2282 ierr = MatAYPX(FDexp,-1.0,Bexp,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 2283 ierr = PetscViewerASCIIPrintf(vstdout,"User-provided matrix minus finite difference Jacobian\n");CHKERRQ(ierr); 2284 if (flag) {ierr = MatView(FDexp,vstdout);CHKERRQ(ierr);} 2285 if (vdraw) { /* Always use contour for the difference */ 2286 ierr = PetscViewerPushFormat(vdraw,PETSC_VIEWER_DRAW_CONTOUR);CHKERRQ(ierr); 2287 ierr = MatView(FDexp,vdraw);CHKERRQ(ierr); 2288 ierr = PetscViewerPopFormat(vdraw);CHKERRQ(ierr); 2289 } 2290 if (flag_contour) {ierr = PetscViewerPopFormat(vdraw);CHKERRQ(ierr);} 2291 ierr = PetscViewerDestroy(&vdraw);CHKERRQ(ierr); 2292 ierr = MatDestroy(&Bexp_mine);CHKERRQ(ierr); 2293 ierr = MatDestroy(&FDexp);CHKERRQ(ierr); 2294 } 2295 } 2296 { 2297 PetscBool flag = PETSC_FALSE,flag_display = PETSC_FALSE,flag_draw = PETSC_FALSE,flag_contour = PETSC_FALSE,flag_threshold = PETSC_FALSE; 2298 PetscReal threshold_atol = PETSC_SQRT_MACHINE_EPSILON,threshold_rtol = 10*PETSC_SQRT_MACHINE_EPSILON; 2299 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_coloring",&flag,NULL);CHKERRQ(ierr); 2300 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_coloring_display",&flag_display,NULL);CHKERRQ(ierr); 2301 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_coloring_draw",&flag_draw,NULL);CHKERRQ(ierr); 2302 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_coloring_draw_contour",&flag_contour,NULL);CHKERRQ(ierr); 2303 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_coloring_threshold",&flag_threshold,NULL);CHKERRQ(ierr); 2304 ierr = PetscOptionsGetReal(((PetscObject)snes)->prefix,"-snes_compare_coloring_threshold_rtol",&threshold_rtol,NULL);CHKERRQ(ierr); 2305 ierr = PetscOptionsGetReal(((PetscObject)snes)->prefix,"-snes_compare_coloring_threshold_atol",&threshold_atol,NULL);CHKERRQ(ierr); 2306 if (flag || flag_display || flag_draw || flag_contour || flag_threshold) { 2307 Mat Bfd; 2308 PetscViewer vdraw,vstdout; 2309 MatColoring coloring; 2310 ISColoring iscoloring; 2311 MatFDColoring matfdcoloring; 2312 PetscErrorCode (*func)(SNES,Vec,Vec,void*); 2313 void *funcctx; 2314 PetscReal norm1,norm2,normmax; 2315 2316 ierr = MatDuplicate(B,MAT_DO_NOT_COPY_VALUES,&Bfd);CHKERRQ(ierr); 2317 ierr = MatColoringCreate(Bfd,&coloring);CHKERRQ(ierr); 2318 ierr = MatColoringSetType(coloring,MATCOLORINGSL);CHKERRQ(ierr); 2319 ierr = MatColoringSetFromOptions(coloring);CHKERRQ(ierr); 2320 ierr = MatColoringApply(coloring,&iscoloring);CHKERRQ(ierr); 2321 ierr = MatColoringDestroy(&coloring);CHKERRQ(ierr); 2322 ierr = MatFDColoringCreate(Bfd,iscoloring,&matfdcoloring);CHKERRQ(ierr); 2323 ierr = MatFDColoringSetFromOptions(matfdcoloring);CHKERRQ(ierr); 2324 ierr = MatFDColoringSetUp(Bfd,iscoloring,matfdcoloring);CHKERRQ(ierr); 2325 ierr = ISColoringDestroy(&iscoloring);CHKERRQ(ierr); 2326 2327 /* This method of getting the function is currently unreliable since it doesn't work for DM local functions. */ 2328 ierr = SNESGetFunction(snes,NULL,&func,&funcctx);CHKERRQ(ierr); 2329 ierr = MatFDColoringSetFunction(matfdcoloring,(PetscErrorCode (*)(void))func,funcctx);CHKERRQ(ierr); 2330 ierr = PetscObjectSetOptionsPrefix((PetscObject)matfdcoloring,((PetscObject)snes)->prefix);CHKERRQ(ierr); 2331 ierr = PetscObjectAppendOptionsPrefix((PetscObject)matfdcoloring,"coloring_");CHKERRQ(ierr); 2332 ierr = MatFDColoringSetFromOptions(matfdcoloring);CHKERRQ(ierr); 2333 ierr = MatFDColoringApply(Bfd,matfdcoloring,X,snes);CHKERRQ(ierr); 2334 ierr = MatFDColoringDestroy(&matfdcoloring);CHKERRQ(ierr); 2335 2336 ierr = PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)snes),&vstdout);CHKERRQ(ierr); 2337 if (flag_draw || flag_contour) { 2338 ierr = PetscViewerDrawOpen(PetscObjectComm((PetscObject)snes),0,"Colored Jacobians",PETSC_DECIDE,PETSC_DECIDE,300,300,&vdraw);CHKERRQ(ierr); 2339 if (flag_contour) {ierr = PetscViewerPushFormat(vdraw,PETSC_VIEWER_DRAW_CONTOUR);CHKERRQ(ierr);} 2340 } else vdraw = NULL; 2341 ierr = PetscViewerASCIIPrintf(vstdout,"Explicit preconditioning Jacobian\n");CHKERRQ(ierr); 2342 if (flag_display) {ierr = MatView(B,vstdout);CHKERRQ(ierr);} 2343 if (vdraw) {ierr = MatView(B,vdraw);CHKERRQ(ierr);} 2344 ierr = PetscViewerASCIIPrintf(vstdout,"Colored Finite difference Jacobian\n");CHKERRQ(ierr); 2345 if (flag_display) {ierr = MatView(Bfd,vstdout);CHKERRQ(ierr);} 2346 if (vdraw) {ierr = MatView(Bfd,vdraw);CHKERRQ(ierr);} 2347 ierr = MatAYPX(Bfd,-1.0,B,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 2348 ierr = MatNorm(Bfd,NORM_1,&norm1);CHKERRQ(ierr); 2349 ierr = MatNorm(Bfd,NORM_FROBENIUS,&norm2);CHKERRQ(ierr); 2350 ierr = MatNorm(Bfd,NORM_MAX,&normmax);CHKERRQ(ierr); 2351 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); 2352 if (flag_display) {ierr = MatView(Bfd,vstdout);CHKERRQ(ierr);} 2353 if (vdraw) { /* Always use contour for the difference */ 2354 ierr = PetscViewerPushFormat(vdraw,PETSC_VIEWER_DRAW_CONTOUR);CHKERRQ(ierr); 2355 ierr = MatView(Bfd,vdraw);CHKERRQ(ierr); 2356 ierr = PetscViewerPopFormat(vdraw);CHKERRQ(ierr); 2357 } 2358 if (flag_contour) {ierr = PetscViewerPopFormat(vdraw);CHKERRQ(ierr);} 2359 2360 if (flag_threshold) { 2361 PetscInt bs,rstart,rend,i; 2362 ierr = MatGetBlockSize(B,&bs);CHKERRQ(ierr); 2363 ierr = MatGetOwnershipRange(B,&rstart,&rend);CHKERRQ(ierr); 2364 for (i=rstart; i<rend; i++) { 2365 const PetscScalar *ba,*ca; 2366 const PetscInt *bj,*cj; 2367 PetscInt bn,cn,j,maxentrycol = -1,maxdiffcol = -1,maxrdiffcol = -1; 2368 PetscReal maxentry = 0,maxdiff = 0,maxrdiff = 0; 2369 ierr = MatGetRow(B,i,&bn,&bj,&ba);CHKERRQ(ierr); 2370 ierr = MatGetRow(Bfd,i,&cn,&cj,&ca);CHKERRQ(ierr); 2371 if (bn != cn) SETERRQ(((PetscObject)A)->comm,PETSC_ERR_PLIB,"Unexpected different nonzero pattern in -snes_compare_coloring_threshold"); 2372 for (j=0; j<bn; j++) { 2373 PetscReal rdiff = PetscAbsScalar(ca[j]) / (threshold_atol + threshold_rtol*PetscAbsScalar(ba[j])); 2374 if (PetscAbsScalar(ba[j]) > PetscAbs(maxentry)) { 2375 maxentrycol = bj[j]; 2376 maxentry = PetscRealPart(ba[j]); 2377 } 2378 if (PetscAbsScalar(ca[j]) > PetscAbs(maxdiff)) { 2379 maxdiffcol = bj[j]; 2380 maxdiff = PetscRealPart(ca[j]); 2381 } 2382 if (rdiff > maxrdiff) { 2383 maxrdiffcol = bj[j]; 2384 maxrdiff = rdiff; 2385 } 2386 } 2387 if (maxrdiff > 1) { 2388 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); 2389 for (j=0; j<bn; j++) { 2390 PetscReal rdiff; 2391 rdiff = PetscAbsScalar(ca[j]) / (threshold_atol + threshold_rtol*PetscAbsScalar(ba[j])); 2392 if (rdiff > 1) { 2393 ierr = PetscViewerASCIIPrintf(vstdout," (%D,%g:%g)",bj[j],(double)PetscRealPart(ba[j]),(double)PetscRealPart(ca[j]));CHKERRQ(ierr); 2394 } 2395 } 2396 ierr = PetscViewerASCIIPrintf(vstdout,"\n",i,maxentry,maxdiff,maxrdiff);CHKERRQ(ierr); 2397 } 2398 ierr = MatRestoreRow(B,i,&bn,&bj,&ba);CHKERRQ(ierr); 2399 ierr = MatRestoreRow(Bfd,i,&cn,&cj,&ca);CHKERRQ(ierr); 2400 } 2401 } 2402 ierr = PetscViewerDestroy(&vdraw);CHKERRQ(ierr); 2403 ierr = MatDestroy(&Bfd);CHKERRQ(ierr); 2404 } 2405 } 2406 PetscFunctionReturn(0); 2407 } 2408 2409 /*MC 2410 SNESJacobianFunction - Function used to convey the nonlinear Jacobian of the function to be solved by SNES 2411 2412 Synopsis: 2413 #include "petscsnes.h" 2414 PetscErrorCode SNESJacobianFunction(SNES snes,Vec x,Mat Amat,Mat Pmat,void *ctx); 2415 2416 + x - input vector 2417 . Amat - the matrix that defines the (approximate) Jacobian 2418 . Pmat - the matrix to be used in constructing the preconditioner, usually the same as Amat. 2419 - ctx - [optional] user-defined Jacobian context 2420 2421 Level: intermediate 2422 2423 .seealso: SNESSetFunction(), SNESGetFunction(), SNESSetJacobian(), SNESGetJacobian() 2424 M*/ 2425 2426 #undef __FUNCT__ 2427 #define __FUNCT__ "SNESSetJacobian" 2428 /*@C 2429 SNESSetJacobian - Sets the function to compute Jacobian as well as the 2430 location to store the matrix. 2431 2432 Logically Collective on SNES and Mat 2433 2434 Input Parameters: 2435 + snes - the SNES context 2436 . Amat - the matrix that defines the (approximate) Jacobian 2437 . Pmat - the matrix to be used in constructing the preconditioner, usually the same as Amat. 2438 . J - Jacobian evaluation routine (if NULL then SNES retains any previously set value), see SNESJacobianFunction for details 2439 - ctx - [optional] user-defined context for private data for the 2440 Jacobian evaluation routine (may be NULL) (if NULL then SNES retains any previously set value) 2441 2442 Notes: 2443 If the Amat matrix and Pmat matrix are different you must call MatAssemblyBegin/End() on 2444 each matrix. 2445 2446 If using SNESComputeJacobianDefaultColor() to assemble a Jacobian, the ctx argument 2447 must be a MatFDColoring. 2448 2449 Other defect-correction schemes can be used by computing a different matrix in place of the Jacobian. One common 2450 example is to use the "Picard linearization" which only differentiates through the highest order parts of each term. 2451 2452 Level: beginner 2453 2454 .keywords: SNES, nonlinear, set, Jacobian, matrix 2455 2456 .seealso: KSPSetOperators(), SNESSetFunction(), MatMFFDComputeJacobian(), SNESComputeJacobianDefaultColor(), MatStructure, J, 2457 SNESSetPicard(), SNESJacobianFunction 2458 @*/ 2459 PetscErrorCode SNESSetJacobian(SNES snes,Mat Amat,Mat Pmat,PetscErrorCode (*J)(SNES,Vec,Mat,Mat,void*),void *ctx) 2460 { 2461 PetscErrorCode ierr; 2462 DM dm; 2463 2464 PetscFunctionBegin; 2465 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2466 if (Amat) PetscValidHeaderSpecific(Amat,MAT_CLASSID,2); 2467 if (Pmat) PetscValidHeaderSpecific(Pmat,MAT_CLASSID,3); 2468 if (Amat) PetscCheckSameComm(snes,1,Amat,2); 2469 if (Pmat) PetscCheckSameComm(snes,1,Pmat,3); 2470 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 2471 ierr = DMSNESSetJacobian(dm,J,ctx);CHKERRQ(ierr); 2472 if (Amat) { 2473 ierr = PetscObjectReference((PetscObject)Amat);CHKERRQ(ierr); 2474 ierr = MatDestroy(&snes->jacobian);CHKERRQ(ierr); 2475 2476 snes->jacobian = Amat; 2477 } 2478 if (Pmat) { 2479 ierr = PetscObjectReference((PetscObject)Pmat);CHKERRQ(ierr); 2480 ierr = MatDestroy(&snes->jacobian_pre);CHKERRQ(ierr); 2481 2482 snes->jacobian_pre = Pmat; 2483 } 2484 PetscFunctionReturn(0); 2485 } 2486 2487 #undef __FUNCT__ 2488 #define __FUNCT__ "SNESGetJacobian" 2489 /*@C 2490 SNESGetJacobian - Returns the Jacobian matrix and optionally the user 2491 provided context for evaluating the Jacobian. 2492 2493 Not Collective, but Mat object will be parallel if SNES object is 2494 2495 Input Parameter: 2496 . snes - the nonlinear solver context 2497 2498 Output Parameters: 2499 + Amat - location to stash (approximate) Jacobian matrix (or NULL) 2500 . Pmat - location to stash matrix used to compute the preconditioner (or NULL) 2501 . J - location to put Jacobian function (or NULL), see SNESJacobianFunction for details on its calling sequence 2502 - ctx - location to stash Jacobian ctx (or NULL) 2503 2504 Level: advanced 2505 2506 .seealso: SNESSetJacobian(), SNESComputeJacobian(), SNESJacobianFunction, SNESGetFunction() 2507 @*/ 2508 PetscErrorCode SNESGetJacobian(SNES snes,Mat *Amat,Mat *Pmat,PetscErrorCode (**J)(SNES,Vec,Mat,Mat,void*),void **ctx) 2509 { 2510 PetscErrorCode ierr; 2511 DM dm; 2512 DMSNES sdm; 2513 2514 PetscFunctionBegin; 2515 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2516 if (Amat) *Amat = snes->jacobian; 2517 if (Pmat) *Pmat = snes->jacobian_pre; 2518 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 2519 ierr = DMGetDMSNES(dm,&sdm);CHKERRQ(ierr); 2520 if (J) *J = sdm->ops->computejacobian; 2521 if (ctx) *ctx = sdm->jacobianctx; 2522 PetscFunctionReturn(0); 2523 } 2524 2525 #undef __FUNCT__ 2526 #define __FUNCT__ "SNESSetUp" 2527 /*@ 2528 SNESSetUp - Sets up the internal data structures for the later use 2529 of a nonlinear solver. 2530 2531 Collective on SNES 2532 2533 Input Parameters: 2534 . snes - the SNES context 2535 2536 Notes: 2537 For basic use of the SNES solvers the user need not explicitly call 2538 SNESSetUp(), since these actions will automatically occur during 2539 the call to SNESSolve(). However, if one wishes to control this 2540 phase separately, SNESSetUp() should be called after SNESCreate() 2541 and optional routines of the form SNESSetXXX(), but before SNESSolve(). 2542 2543 Level: advanced 2544 2545 .keywords: SNES, nonlinear, setup 2546 2547 .seealso: SNESCreate(), SNESSolve(), SNESDestroy() 2548 @*/ 2549 PetscErrorCode SNESSetUp(SNES snes) 2550 { 2551 PetscErrorCode ierr; 2552 DM dm; 2553 DMSNES sdm; 2554 SNESLineSearch linesearch, pclinesearch; 2555 void *lsprectx,*lspostctx; 2556 PetscErrorCode (*precheck)(SNESLineSearch,Vec,Vec,PetscBool*,void*); 2557 PetscErrorCode (*postcheck)(SNESLineSearch,Vec,Vec,Vec,PetscBool*,PetscBool*,void*); 2558 PetscErrorCode (*func)(SNES,Vec,Vec,void*); 2559 Vec f,fpc; 2560 void *funcctx; 2561 PetscErrorCode (*jac)(SNES,Vec,Mat,Mat,void*); 2562 void *jacctx,*appctx; 2563 Mat j,jpre; 2564 2565 PetscFunctionBegin; 2566 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2567 if (snes->setupcalled) PetscFunctionReturn(0); 2568 2569 if (!((PetscObject)snes)->type_name) { 2570 ierr = SNESSetType(snes,SNESNEWTONLS);CHKERRQ(ierr); 2571 } 2572 2573 ierr = SNESGetFunction(snes,&snes->vec_func,NULL,NULL);CHKERRQ(ierr); 2574 2575 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 2576 ierr = DMGetDMSNES(dm,&sdm);CHKERRQ(ierr); 2577 if (!sdm->ops->computefunction) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"Function never provided to SNES object"); 2578 if (!sdm->ops->computejacobian) { 2579 ierr = DMSNESSetJacobian(dm,SNESComputeJacobianDefaultColor,NULL);CHKERRQ(ierr); 2580 } 2581 if (!snes->vec_func) { 2582 ierr = DMCreateGlobalVector(dm,&snes->vec_func);CHKERRQ(ierr); 2583 } 2584 2585 if (!snes->ksp) { 2586 ierr = SNESGetKSP(snes, &snes->ksp);CHKERRQ(ierr); 2587 } 2588 2589 if (!snes->linesearch) { 2590 ierr = SNESGetLineSearch(snes, &snes->linesearch);CHKERRQ(ierr); 2591 } 2592 ierr = SNESLineSearchSetFunction(snes->linesearch,SNESComputeFunction);CHKERRQ(ierr); 2593 2594 if (snes->pc && (snes->pcside == PC_LEFT)) { 2595 snes->mf = PETSC_TRUE; 2596 snes->mf_operator = PETSC_FALSE; 2597 } 2598 2599 if (snes->pc) { 2600 /* copy the DM over */ 2601 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 2602 ierr = SNESSetDM(snes->pc,dm);CHKERRQ(ierr); 2603 2604 ierr = SNESGetFunction(snes,&f,&func,&funcctx);CHKERRQ(ierr); 2605 ierr = VecDuplicate(f,&fpc);CHKERRQ(ierr); 2606 ierr = SNESSetFunction(snes->pc,fpc,func,funcctx);CHKERRQ(ierr); 2607 ierr = SNESGetJacobian(snes,&j,&jpre,&jac,&jacctx);CHKERRQ(ierr); 2608 ierr = SNESSetJacobian(snes->pc,j,jpre,jac,jacctx);CHKERRQ(ierr); 2609 ierr = SNESGetApplicationContext(snes,&appctx);CHKERRQ(ierr); 2610 ierr = SNESSetApplicationContext(snes->pc,appctx);CHKERRQ(ierr); 2611 ierr = VecDestroy(&fpc);CHKERRQ(ierr); 2612 2613 /* copy the function pointers over */ 2614 ierr = PetscObjectCopyFortranFunctionPointers((PetscObject)snes,(PetscObject)snes->pc);CHKERRQ(ierr); 2615 2616 /* default to 1 iteration */ 2617 ierr = SNESSetTolerances(snes->pc,0.0,0.0,0.0,1,snes->pc->max_funcs);CHKERRQ(ierr); 2618 if (snes->pcside==PC_RIGHT) { 2619 ierr = SNESSetNormSchedule(snes->pc,SNES_NORM_FINAL_ONLY);CHKERRQ(ierr); 2620 } else { 2621 ierr = SNESSetNormSchedule(snes->pc,SNES_NORM_NONE);CHKERRQ(ierr); 2622 } 2623 ierr = SNESSetFromOptions(snes->pc);CHKERRQ(ierr); 2624 2625 /* copy the line search context over */ 2626 ierr = SNESGetLineSearch(snes,&linesearch);CHKERRQ(ierr); 2627 ierr = SNESGetLineSearch(snes->pc,&pclinesearch);CHKERRQ(ierr); 2628 ierr = SNESLineSearchGetPreCheck(linesearch,&precheck,&lsprectx);CHKERRQ(ierr); 2629 ierr = SNESLineSearchGetPostCheck(linesearch,&postcheck,&lspostctx);CHKERRQ(ierr); 2630 ierr = SNESLineSearchSetPreCheck(pclinesearch,precheck,lsprectx);CHKERRQ(ierr); 2631 ierr = SNESLineSearchSetPostCheck(pclinesearch,postcheck,lspostctx);CHKERRQ(ierr); 2632 ierr = PetscObjectCopyFortranFunctionPointers((PetscObject)linesearch, (PetscObject)pclinesearch);CHKERRQ(ierr); 2633 } 2634 if (snes->mf) { 2635 ierr = SNESSetUpMatrixFree_Private(snes, snes->mf_operator, snes->mf_version);CHKERRQ(ierr); 2636 } 2637 if (snes->ops->usercompute && !snes->user) { 2638 ierr = (*snes->ops->usercompute)(snes,(void**)&snes->user);CHKERRQ(ierr); 2639 } 2640 2641 snes->jac_iter = 0; 2642 snes->pre_iter = 0; 2643 2644 if (snes->ops->setup) { 2645 ierr = (*snes->ops->setup)(snes);CHKERRQ(ierr); 2646 } 2647 2648 if (snes->pc && (snes->pcside == PC_LEFT)) { 2649 if (snes->functype == SNES_FUNCTION_PRECONDITIONED) { 2650 ierr = SNESGetLineSearch(snes,&linesearch);CHKERRQ(ierr); 2651 ierr = SNESLineSearchSetFunction(linesearch,SNESComputeFunctionDefaultNPC);CHKERRQ(ierr); 2652 } 2653 } 2654 2655 snes->setupcalled = PETSC_TRUE; 2656 PetscFunctionReturn(0); 2657 } 2658 2659 #undef __FUNCT__ 2660 #define __FUNCT__ "SNESReset" 2661 /*@ 2662 SNESReset - Resets a SNES context to the snessetupcalled = 0 state and removes any allocated Vecs and Mats 2663 2664 Collective on SNES 2665 2666 Input Parameter: 2667 . snes - iterative context obtained from SNESCreate() 2668 2669 Level: intermediate 2670 2671 Notes: Also calls the application context destroy routine set with SNESSetComputeApplicationContext() 2672 2673 .keywords: SNES, destroy 2674 2675 .seealso: SNESCreate(), SNESSetUp(), SNESSolve() 2676 @*/ 2677 PetscErrorCode SNESReset(SNES snes) 2678 { 2679 PetscErrorCode ierr; 2680 2681 PetscFunctionBegin; 2682 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2683 if (snes->ops->userdestroy && snes->user) { 2684 ierr = (*snes->ops->userdestroy)((void**)&snes->user);CHKERRQ(ierr); 2685 snes->user = NULL; 2686 } 2687 if (snes->pc) { 2688 ierr = SNESReset(snes->pc);CHKERRQ(ierr); 2689 } 2690 2691 if (snes->ops->reset) { 2692 ierr = (*snes->ops->reset)(snes);CHKERRQ(ierr); 2693 } 2694 if (snes->ksp) { 2695 ierr = KSPReset(snes->ksp);CHKERRQ(ierr); 2696 } 2697 2698 if (snes->linesearch) { 2699 ierr = SNESLineSearchReset(snes->linesearch);CHKERRQ(ierr); 2700 } 2701 2702 ierr = VecDestroy(&snes->vec_rhs);CHKERRQ(ierr); 2703 ierr = VecDestroy(&snes->vec_sol);CHKERRQ(ierr); 2704 ierr = VecDestroy(&snes->vec_sol_update);CHKERRQ(ierr); 2705 ierr = VecDestroy(&snes->vec_func);CHKERRQ(ierr); 2706 ierr = MatDestroy(&snes->jacobian);CHKERRQ(ierr); 2707 ierr = MatDestroy(&snes->jacobian_pre);CHKERRQ(ierr); 2708 ierr = VecDestroyVecs(snes->nwork,&snes->work);CHKERRQ(ierr); 2709 ierr = VecDestroyVecs(snes->nvwork,&snes->vwork);CHKERRQ(ierr); 2710 2711 snes->nwork = snes->nvwork = 0; 2712 snes->setupcalled = PETSC_FALSE; 2713 PetscFunctionReturn(0); 2714 } 2715 2716 #undef __FUNCT__ 2717 #define __FUNCT__ "SNESDestroy" 2718 /*@ 2719 SNESDestroy - Destroys the nonlinear solver context that was created 2720 with SNESCreate(). 2721 2722 Collective on SNES 2723 2724 Input Parameter: 2725 . snes - the SNES context 2726 2727 Level: beginner 2728 2729 .keywords: SNES, nonlinear, destroy 2730 2731 .seealso: SNESCreate(), SNESSolve() 2732 @*/ 2733 PetscErrorCode SNESDestroy(SNES *snes) 2734 { 2735 PetscErrorCode ierr; 2736 2737 PetscFunctionBegin; 2738 if (!*snes) PetscFunctionReturn(0); 2739 PetscValidHeaderSpecific((*snes),SNES_CLASSID,1); 2740 if (--((PetscObject)(*snes))->refct > 0) {*snes = 0; PetscFunctionReturn(0);} 2741 2742 ierr = SNESReset((*snes));CHKERRQ(ierr); 2743 ierr = SNESDestroy(&(*snes)->pc);CHKERRQ(ierr); 2744 2745 /* if memory was published with SAWs then destroy it */ 2746 ierr = PetscObjectSAWsViewOff((PetscObject)*snes);CHKERRQ(ierr); 2747 if ((*snes)->ops->destroy) {ierr = (*((*snes))->ops->destroy)((*snes));CHKERRQ(ierr);} 2748 2749 ierr = DMDestroy(&(*snes)->dm);CHKERRQ(ierr); 2750 ierr = KSPDestroy(&(*snes)->ksp);CHKERRQ(ierr); 2751 ierr = SNESLineSearchDestroy(&(*snes)->linesearch);CHKERRQ(ierr); 2752 2753 ierr = PetscFree((*snes)->kspconvctx);CHKERRQ(ierr); 2754 if ((*snes)->ops->convergeddestroy) { 2755 ierr = (*(*snes)->ops->convergeddestroy)((*snes)->cnvP);CHKERRQ(ierr); 2756 } 2757 if ((*snes)->conv_malloc) { 2758 ierr = PetscFree((*snes)->conv_hist);CHKERRQ(ierr); 2759 ierr = PetscFree((*snes)->conv_hist_its);CHKERRQ(ierr); 2760 } 2761 ierr = SNESMonitorCancel((*snes));CHKERRQ(ierr); 2762 ierr = PetscHeaderDestroy(snes);CHKERRQ(ierr); 2763 PetscFunctionReturn(0); 2764 } 2765 2766 /* ----------- Routines to set solver parameters ---------- */ 2767 2768 #undef __FUNCT__ 2769 #define __FUNCT__ "SNESSetLagPreconditioner" 2770 /*@ 2771 SNESSetLagPreconditioner - Determines when the preconditioner is rebuilt in the nonlinear solve. 2772 2773 Logically Collective on SNES 2774 2775 Input Parameters: 2776 + snes - the SNES context 2777 - lag - -1 indicates NEVER rebuild, 1 means rebuild every time the Jacobian is computed within a single nonlinear solve, 2 means every second time 2778 the Jacobian is built etc. -2 indicates rebuild preconditioner at next chance but then never rebuild after that 2779 2780 Options Database Keys: 2781 . -snes_lag_preconditioner <lag> 2782 2783 Notes: 2784 The default is 1 2785 The preconditioner is ALWAYS built in the first iteration of a nonlinear solve unless lag is -1 2786 If -1 is used before the very first nonlinear solve the preconditioner is still built because there is no previous preconditioner to use 2787 2788 Level: intermediate 2789 2790 .keywords: SNES, nonlinear, set, convergence, tolerances 2791 2792 .seealso: SNESSetTrustRegionTolerance(), SNESGetLagPreconditioner(), SNESSetLagJacobian(), SNESGetLagJacobian() 2793 2794 @*/ 2795 PetscErrorCode SNESSetLagPreconditioner(SNES snes,PetscInt lag) 2796 { 2797 PetscFunctionBegin; 2798 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2799 if (lag < -2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Lag must be -2, -1, 1 or greater"); 2800 if (!lag) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Lag cannot be 0"); 2801 PetscValidLogicalCollectiveInt(snes,lag,2); 2802 snes->lagpreconditioner = lag; 2803 PetscFunctionReturn(0); 2804 } 2805 2806 #undef __FUNCT__ 2807 #define __FUNCT__ "SNESSetGridSequence" 2808 /*@ 2809 SNESSetGridSequence - sets the number of steps of grid sequencing that SNES does 2810 2811 Logically Collective on SNES 2812 2813 Input Parameters: 2814 + snes - the SNES context 2815 - steps - the number of refinements to do, defaults to 0 2816 2817 Options Database Keys: 2818 . -snes_grid_sequence <steps> 2819 2820 Level: intermediate 2821 2822 Notes: 2823 Use SNESGetSolution() to extract the fine grid solution after grid sequencing. 2824 2825 .keywords: SNES, nonlinear, set, convergence, tolerances 2826 2827 .seealso: SNESSetTrustRegionTolerance(), SNESGetLagPreconditioner(), SNESSetLagJacobian(), SNESGetLagJacobian(), SNESGetGridSequence() 2828 2829 @*/ 2830 PetscErrorCode SNESSetGridSequence(SNES snes,PetscInt steps) 2831 { 2832 PetscFunctionBegin; 2833 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2834 PetscValidLogicalCollectiveInt(snes,steps,2); 2835 snes->gridsequence = steps; 2836 PetscFunctionReturn(0); 2837 } 2838 2839 #undef __FUNCT__ 2840 #define __FUNCT__ "SNESGetGridSequence" 2841 /*@ 2842 SNESGetGridSequence - gets the number of steps of grid sequencing that SNES does 2843 2844 Logically Collective on SNES 2845 2846 Input Parameter: 2847 . snes - the SNES context 2848 2849 Output Parameter: 2850 . steps - the number of refinements to do, defaults to 0 2851 2852 Options Database Keys: 2853 . -snes_grid_sequence <steps> 2854 2855 Level: intermediate 2856 2857 Notes: 2858 Use SNESGetSolution() to extract the fine grid solution after grid sequencing. 2859 2860 .keywords: SNES, nonlinear, set, convergence, tolerances 2861 2862 .seealso: SNESSetTrustRegionTolerance(), SNESGetLagPreconditioner(), SNESSetLagJacobian(), SNESGetLagJacobian(), SNESSetGridSequence() 2863 2864 @*/ 2865 PetscErrorCode SNESGetGridSequence(SNES snes,PetscInt *steps) 2866 { 2867 PetscFunctionBegin; 2868 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2869 *steps = snes->gridsequence; 2870 PetscFunctionReturn(0); 2871 } 2872 2873 #undef __FUNCT__ 2874 #define __FUNCT__ "SNESGetLagPreconditioner" 2875 /*@ 2876 SNESGetLagPreconditioner - Indicates how often the preconditioner is rebuilt 2877 2878 Not Collective 2879 2880 Input Parameter: 2881 . snes - the SNES context 2882 2883 Output Parameter: 2884 . lag - -1 indicates NEVER rebuild, 1 means rebuild every time the Jacobian is computed within a single nonlinear solve, 2 means every second time 2885 the Jacobian is built etc. -2 indicates rebuild preconditioner at next chance but then never rebuild after that 2886 2887 Options Database Keys: 2888 . -snes_lag_preconditioner <lag> 2889 2890 Notes: 2891 The default is 1 2892 The preconditioner is ALWAYS built in the first iteration of a nonlinear solve unless lag is -1 2893 2894 Level: intermediate 2895 2896 .keywords: SNES, nonlinear, set, convergence, tolerances 2897 2898 .seealso: SNESSetTrustRegionTolerance(), SNESSetLagPreconditioner() 2899 2900 @*/ 2901 PetscErrorCode SNESGetLagPreconditioner(SNES snes,PetscInt *lag) 2902 { 2903 PetscFunctionBegin; 2904 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2905 *lag = snes->lagpreconditioner; 2906 PetscFunctionReturn(0); 2907 } 2908 2909 #undef __FUNCT__ 2910 #define __FUNCT__ "SNESSetLagJacobian" 2911 /*@ 2912 SNESSetLagJacobian - Determines when the Jacobian is rebuilt in the nonlinear solve. See SNESSetLagPreconditioner() for determining how 2913 often the preconditioner is rebuilt. 2914 2915 Logically Collective on SNES 2916 2917 Input Parameters: 2918 + snes - the SNES context 2919 - lag - -1 indicates NEVER rebuild, 1 means rebuild every time the Jacobian is computed within a single nonlinear solve, 2 means every second time 2920 the Jacobian is built etc. -2 means rebuild at next chance but then never again 2921 2922 Options Database Keys: 2923 . -snes_lag_jacobian <lag> 2924 2925 Notes: 2926 The default is 1 2927 The Jacobian is ALWAYS built in the first iteration of a nonlinear solve unless lag is -1 2928 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 2929 at the next Newton step but never again (unless it is reset to another value) 2930 2931 Level: intermediate 2932 2933 .keywords: SNES, nonlinear, set, convergence, tolerances 2934 2935 .seealso: SNESSetTrustRegionTolerance(), SNESGetLagPreconditioner(), SNESSetLagPreconditioner(), SNESGetLagJacobian() 2936 2937 @*/ 2938 PetscErrorCode SNESSetLagJacobian(SNES snes,PetscInt lag) 2939 { 2940 PetscFunctionBegin; 2941 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2942 if (lag < -2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Lag must be -2, -1, 1 or greater"); 2943 if (!lag) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Lag cannot be 0"); 2944 PetscValidLogicalCollectiveInt(snes,lag,2); 2945 snes->lagjacobian = lag; 2946 PetscFunctionReturn(0); 2947 } 2948 2949 #undef __FUNCT__ 2950 #define __FUNCT__ "SNESGetLagJacobian" 2951 /*@ 2952 SNESGetLagJacobian - Indicates how often the Jacobian is rebuilt. See SNESGetLagPreconditioner() to determine when the preconditioner is rebuilt 2953 2954 Not Collective 2955 2956 Input Parameter: 2957 . snes - the SNES context 2958 2959 Output Parameter: 2960 . lag - -1 indicates NEVER rebuild, 1 means rebuild every time the Jacobian is computed within a single nonlinear solve, 2 means every second time 2961 the Jacobian is built etc. 2962 2963 Options Database Keys: 2964 . -snes_lag_jacobian <lag> 2965 2966 Notes: 2967 The default is 1 2968 The jacobian is ALWAYS built in the first iteration of a nonlinear solve unless lag is -1 2969 2970 Level: intermediate 2971 2972 .keywords: SNES, nonlinear, set, convergence, tolerances 2973 2974 .seealso: SNESSetTrustRegionTolerance(), SNESSetLagJacobian(), SNESSetLagPreconditioner(), SNESGetLagPreconditioner() 2975 2976 @*/ 2977 PetscErrorCode SNESGetLagJacobian(SNES snes,PetscInt *lag) 2978 { 2979 PetscFunctionBegin; 2980 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2981 *lag = snes->lagjacobian; 2982 PetscFunctionReturn(0); 2983 } 2984 2985 #undef __FUNCT__ 2986 #define __FUNCT__ "SNESSetLagJacobianPersists" 2987 /*@ 2988 SNESSetLagJacobianPersists - Set whether or not the Jacobian lagging persists through multiple solves 2989 2990 Logically collective on SNES 2991 2992 Input Parameter: 2993 + snes - the SNES context 2994 - flg - jacobian lagging persists if true 2995 2996 Options Database Keys: 2997 . -snes_lag_jacobian_persists <flg> 2998 2999 Notes: This is useful both for nonlinear preconditioning, where it's appropriate to have the Jacobian be stale by 3000 several solves, and for implicit time-stepping, where Jacobian lagging in the inner nonlinear solve over several 3001 timesteps may present huge efficiency gains. 3002 3003 Level: developer 3004 3005 .keywords: SNES, nonlinear, lag 3006 3007 .seealso: SNESSetLagPreconditionerPersists(), SNESSetLagJacobian(), SNESGetLagJacobian(), SNESGetNPC() 3008 3009 @*/ 3010 PetscErrorCode SNESSetLagJacobianPersists(SNES snes,PetscBool flg) 3011 { 3012 PetscFunctionBegin; 3013 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3014 PetscValidLogicalCollectiveBool(snes,flg,2); 3015 snes->lagjac_persist = flg; 3016 PetscFunctionReturn(0); 3017 } 3018 3019 #undef __FUNCT__ 3020 #define __FUNCT__ "SNESSetLagPreconditionerPersists" 3021 /*@ 3022 SNESSetLagPreconditionerPersists - Set whether or not the preconditioner lagging persists through multiple solves 3023 3024 Logically Collective on SNES 3025 3026 Input Parameter: 3027 + snes - the SNES context 3028 - flg - preconditioner lagging persists if true 3029 3030 Options Database Keys: 3031 . -snes_lag_jacobian_persists <flg> 3032 3033 Notes: This is useful both for nonlinear preconditioning, where it's appropriate to have the preconditioner be stale 3034 by several solves, and for implicit time-stepping, where preconditioner lagging in the inner nonlinear solve over 3035 several timesteps may present huge efficiency gains. 3036 3037 Level: developer 3038 3039 .keywords: SNES, nonlinear, lag 3040 3041 .seealso: SNESSetLagJacobianPersists(), SNESSetLagJacobian(), SNESGetLagJacobian(), SNESGetNPC() 3042 3043 @*/ 3044 PetscErrorCode SNESSetLagPreconditionerPersists(SNES snes,PetscBool flg) 3045 { 3046 PetscFunctionBegin; 3047 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3048 PetscValidLogicalCollectiveBool(snes,flg,2); 3049 snes->lagpre_persist = flg; 3050 PetscFunctionReturn(0); 3051 } 3052 3053 #undef __FUNCT__ 3054 #define __FUNCT__ "SNESSetTolerances" 3055 /*@ 3056 SNESSetTolerances - Sets various parameters used in convergence tests. 3057 3058 Logically Collective on SNES 3059 3060 Input Parameters: 3061 + snes - the SNES context 3062 . abstol - absolute convergence tolerance 3063 . rtol - relative convergence tolerance 3064 . stol - convergence tolerance in terms of the norm of the change in the solution between steps, || delta x || < stol*|| x || 3065 . maxit - maximum number of iterations 3066 - maxf - maximum number of function evaluations 3067 3068 Options Database Keys: 3069 + -snes_atol <abstol> - Sets abstol 3070 . -snes_rtol <rtol> - Sets rtol 3071 . -snes_stol <stol> - Sets stol 3072 . -snes_max_it <maxit> - Sets maxit 3073 - -snes_max_funcs <maxf> - Sets maxf 3074 3075 Notes: 3076 The default maximum number of iterations is 50. 3077 The default maximum number of function evaluations is 1000. 3078 3079 Level: intermediate 3080 3081 .keywords: SNES, nonlinear, set, convergence, tolerances 3082 3083 .seealso: SNESSetTrustRegionTolerance() 3084 @*/ 3085 PetscErrorCode SNESSetTolerances(SNES snes,PetscReal abstol,PetscReal rtol,PetscReal stol,PetscInt maxit,PetscInt maxf) 3086 { 3087 PetscFunctionBegin; 3088 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3089 PetscValidLogicalCollectiveReal(snes,abstol,2); 3090 PetscValidLogicalCollectiveReal(snes,rtol,3); 3091 PetscValidLogicalCollectiveReal(snes,stol,4); 3092 PetscValidLogicalCollectiveInt(snes,maxit,5); 3093 PetscValidLogicalCollectiveInt(snes,maxf,6); 3094 3095 if (abstol != PETSC_DEFAULT) { 3096 if (abstol < 0.0) SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_OUTOFRANGE,"Absolute tolerance %g must be non-negative",(double)abstol); 3097 snes->abstol = abstol; 3098 } 3099 if (rtol != PETSC_DEFAULT) { 3100 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); 3101 snes->rtol = rtol; 3102 } 3103 if (stol != PETSC_DEFAULT) { 3104 if (stol < 0.0) SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_OUTOFRANGE,"Step tolerance %g must be non-negative",(double)stol); 3105 snes->stol = stol; 3106 } 3107 if (maxit != PETSC_DEFAULT) { 3108 if (maxit < 0) SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_OUTOFRANGE,"Maximum number of iterations %D must be non-negative",maxit); 3109 snes->max_its = maxit; 3110 } 3111 if (maxf != PETSC_DEFAULT) { 3112 if (maxf < 0) SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_OUTOFRANGE,"Maximum number of function evaluations %D must be non-negative",maxf); 3113 snes->max_funcs = maxf; 3114 } 3115 snes->tolerancesset = PETSC_TRUE; 3116 PetscFunctionReturn(0); 3117 } 3118 3119 #undef __FUNCT__ 3120 #define __FUNCT__ "SNESGetTolerances" 3121 /*@ 3122 SNESGetTolerances - Gets various parameters used in convergence tests. 3123 3124 Not Collective 3125 3126 Input Parameters: 3127 + snes - the SNES context 3128 . atol - absolute convergence tolerance 3129 . rtol - relative convergence tolerance 3130 . stol - convergence tolerance in terms of the norm 3131 of the change in the solution between steps 3132 . maxit - maximum number of iterations 3133 - maxf - maximum number of function evaluations 3134 3135 Notes: 3136 The user can specify NULL for any parameter that is not needed. 3137 3138 Level: intermediate 3139 3140 .keywords: SNES, nonlinear, get, convergence, tolerances 3141 3142 .seealso: SNESSetTolerances() 3143 @*/ 3144 PetscErrorCode SNESGetTolerances(SNES snes,PetscReal *atol,PetscReal *rtol,PetscReal *stol,PetscInt *maxit,PetscInt *maxf) 3145 { 3146 PetscFunctionBegin; 3147 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3148 if (atol) *atol = snes->abstol; 3149 if (rtol) *rtol = snes->rtol; 3150 if (stol) *stol = snes->stol; 3151 if (maxit) *maxit = snes->max_its; 3152 if (maxf) *maxf = snes->max_funcs; 3153 PetscFunctionReturn(0); 3154 } 3155 3156 #undef __FUNCT__ 3157 #define __FUNCT__ "SNESSetTrustRegionTolerance" 3158 /*@ 3159 SNESSetTrustRegionTolerance - Sets the trust region parameter tolerance. 3160 3161 Logically Collective on SNES 3162 3163 Input Parameters: 3164 + snes - the SNES context 3165 - tol - tolerance 3166 3167 Options Database Key: 3168 . -snes_trtol <tol> - Sets tol 3169 3170 Level: intermediate 3171 3172 .keywords: SNES, nonlinear, set, trust region, tolerance 3173 3174 .seealso: SNESSetTolerances() 3175 @*/ 3176 PetscErrorCode SNESSetTrustRegionTolerance(SNES snes,PetscReal tol) 3177 { 3178 PetscFunctionBegin; 3179 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3180 PetscValidLogicalCollectiveReal(snes,tol,2); 3181 snes->deltatol = tol; 3182 PetscFunctionReturn(0); 3183 } 3184 3185 /* 3186 Duplicate the lg monitors for SNES from KSP; for some reason with 3187 dynamic libraries things don't work under Sun4 if we just use 3188 macros instead of functions 3189 */ 3190 #undef __FUNCT__ 3191 #define __FUNCT__ "SNESMonitorLGResidualNorm" 3192 PetscErrorCode SNESMonitorLGResidualNorm(SNES snes,PetscInt it,PetscReal norm,PetscObject *objs) 3193 { 3194 PetscErrorCode ierr; 3195 3196 PetscFunctionBegin; 3197 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3198 ierr = KSPMonitorLGResidualNorm((KSP)snes,it,norm,objs);CHKERRQ(ierr); 3199 PetscFunctionReturn(0); 3200 } 3201 3202 #undef __FUNCT__ 3203 #define __FUNCT__ "SNESMonitorLGCreate" 3204 PetscErrorCode SNESMonitorLGCreate(const char host[],const char label[],int x,int y,int m,int n,PetscObject **draw) 3205 { 3206 PetscErrorCode ierr; 3207 3208 PetscFunctionBegin; 3209 ierr = KSPMonitorLGResidualNormCreate(host,label,x,y,m,n,draw);CHKERRQ(ierr); 3210 PetscFunctionReturn(0); 3211 } 3212 3213 #undef __FUNCT__ 3214 #define __FUNCT__ "SNESMonitorLGDestroy" 3215 PetscErrorCode SNESMonitorLGDestroy(PetscObject **objs) 3216 { 3217 PetscErrorCode ierr; 3218 3219 PetscFunctionBegin; 3220 ierr = KSPMonitorLGResidualNormDestroy(objs);CHKERRQ(ierr); 3221 PetscFunctionReturn(0); 3222 } 3223 3224 extern PetscErrorCode SNESMonitorRange_Private(SNES,PetscInt,PetscReal*); 3225 #undef __FUNCT__ 3226 #define __FUNCT__ "SNESMonitorLGRange" 3227 PetscErrorCode SNESMonitorLGRange(SNES snes,PetscInt n,PetscReal rnorm,void *monctx) 3228 { 3229 PetscDrawLG lg; 3230 PetscErrorCode ierr; 3231 PetscReal x,y,per; 3232 PetscViewer v = (PetscViewer)monctx; 3233 static PetscReal prev; /* should be in the context */ 3234 PetscDraw draw; 3235 3236 PetscFunctionBegin; 3237 ierr = PetscViewerDrawGetDrawLG(v,0,&lg);CHKERRQ(ierr); 3238 if (!n) {ierr = PetscDrawLGReset(lg);CHKERRQ(ierr);} 3239 ierr = PetscDrawLGGetDraw(lg,&draw);CHKERRQ(ierr); 3240 ierr = PetscDrawSetTitle(draw,"Residual norm");CHKERRQ(ierr); 3241 x = (PetscReal)n; 3242 if (rnorm > 0.0) y = PetscLog10Real(rnorm); 3243 else y = -15.0; 3244 ierr = PetscDrawLGAddPoint(lg,&x,&y);CHKERRQ(ierr); 3245 if (n < 20 || !(n % 5)) { 3246 ierr = PetscDrawLGDraw(lg);CHKERRQ(ierr); 3247 } 3248 3249 ierr = PetscViewerDrawGetDrawLG(v,1,&lg);CHKERRQ(ierr); 3250 if (!n) {ierr = PetscDrawLGReset(lg);CHKERRQ(ierr);} 3251 ierr = PetscDrawLGGetDraw(lg,&draw);CHKERRQ(ierr); 3252 ierr = PetscDrawSetTitle(draw,"% elemts > .2*max elemt");CHKERRQ(ierr); 3253 ierr = SNESMonitorRange_Private(snes,n,&per);CHKERRQ(ierr); 3254 x = (PetscReal)n; 3255 y = 100.0*per; 3256 ierr = PetscDrawLGAddPoint(lg,&x,&y);CHKERRQ(ierr); 3257 if (n < 20 || !(n % 5)) { 3258 ierr = PetscDrawLGDraw(lg);CHKERRQ(ierr); 3259 } 3260 3261 ierr = PetscViewerDrawGetDrawLG(v,2,&lg);CHKERRQ(ierr); 3262 if (!n) {prev = rnorm;ierr = PetscDrawLGReset(lg);CHKERRQ(ierr);} 3263 ierr = PetscDrawLGGetDraw(lg,&draw);CHKERRQ(ierr); 3264 ierr = PetscDrawSetTitle(draw,"(norm -oldnorm)/oldnorm");CHKERRQ(ierr); 3265 x = (PetscReal)n; 3266 y = (prev - rnorm)/prev; 3267 ierr = PetscDrawLGAddPoint(lg,&x,&y);CHKERRQ(ierr); 3268 if (n < 20 || !(n % 5)) { 3269 ierr = PetscDrawLGDraw(lg);CHKERRQ(ierr); 3270 } 3271 3272 ierr = PetscViewerDrawGetDrawLG(v,3,&lg);CHKERRQ(ierr); 3273 if (!n) {ierr = PetscDrawLGReset(lg);CHKERRQ(ierr);} 3274 ierr = PetscDrawLGGetDraw(lg,&draw);CHKERRQ(ierr); 3275 ierr = PetscDrawSetTitle(draw,"(norm -oldnorm)/oldnorm*(% > .2 max)");CHKERRQ(ierr); 3276 x = (PetscReal)n; 3277 y = (prev - rnorm)/(prev*per); 3278 if (n > 2) { /*skip initial crazy value */ 3279 ierr = PetscDrawLGAddPoint(lg,&x,&y);CHKERRQ(ierr); 3280 } 3281 if (n < 20 || !(n % 5)) { 3282 ierr = PetscDrawLGDraw(lg);CHKERRQ(ierr); 3283 } 3284 prev = rnorm; 3285 PetscFunctionReturn(0); 3286 } 3287 3288 #undef __FUNCT__ 3289 #define __FUNCT__ "SNESMonitor" 3290 /*@ 3291 SNESMonitor - runs the user provided monitor routines, if they exist 3292 3293 Collective on SNES 3294 3295 Input Parameters: 3296 + snes - nonlinear solver context obtained from SNESCreate() 3297 . iter - iteration number 3298 - rnorm - relative norm of the residual 3299 3300 Notes: 3301 This routine is called by the SNES implementations. 3302 It does not typically need to be called by the user. 3303 3304 Level: developer 3305 3306 .seealso: SNESMonitorSet() 3307 @*/ 3308 PetscErrorCode SNESMonitor(SNES snes,PetscInt iter,PetscReal rnorm) 3309 { 3310 PetscErrorCode ierr; 3311 PetscInt i,n = snes->numbermonitors; 3312 3313 PetscFunctionBegin; 3314 ierr = VecLockPush(snes->vec_sol);CHKERRQ(ierr); 3315 for (i=0; i<n; i++) { 3316 ierr = (*snes->monitor[i])(snes,iter,rnorm,snes->monitorcontext[i]);CHKERRQ(ierr); 3317 } 3318 ierr = VecLockPop(snes->vec_sol);CHKERRQ(ierr); 3319 PetscFunctionReturn(0); 3320 } 3321 3322 /* ------------ Routines to set performance monitoring options ----------- */ 3323 3324 /*MC 3325 SNESMonitorFunction - functional form passed to SNESMonitorSet() to monitor convergence of nonlinear solver 3326 3327 Synopsis: 3328 #include <petscsnes.h> 3329 $ PetscErrorCode SNESMonitorFunction(SNES snes,PetscInt its, PetscReal norm,void *mctx) 3330 3331 + snes - the SNES context 3332 . its - iteration number 3333 . norm - 2-norm function value (may be estimated) 3334 - mctx - [optional] monitoring context 3335 3336 Level: advanced 3337 3338 .seealso: SNESMonitorSet(), SNESMonitorGet() 3339 M*/ 3340 3341 #undef __FUNCT__ 3342 #define __FUNCT__ "SNESMonitorSet" 3343 /*@C 3344 SNESMonitorSet - Sets an ADDITIONAL function that is to be used at every 3345 iteration of the nonlinear solver to display the iteration's 3346 progress. 3347 3348 Logically Collective on SNES 3349 3350 Input Parameters: 3351 + snes - the SNES context 3352 . f - the monitor function, see SNESMonitorFunction for the calling sequence 3353 . mctx - [optional] user-defined context for private data for the 3354 monitor routine (use NULL if no context is desired) 3355 - monitordestroy - [optional] routine that frees monitor context 3356 (may be NULL) 3357 3358 Options Database Keys: 3359 + -snes_monitor - sets SNESMonitorDefault() 3360 . -snes_monitor_lg_residualnorm - sets line graph monitor, 3361 uses SNESMonitorLGCreate() 3362 - -snes_monitor_cancel - cancels all monitors that have 3363 been hardwired into a code by 3364 calls to SNESMonitorSet(), but 3365 does not cancel those set via 3366 the options database. 3367 3368 Notes: 3369 Several different monitoring routines may be set by calling 3370 SNESMonitorSet() multiple times; all will be called in the 3371 order in which they were set. 3372 3373 Fortran notes: Only a single monitor function can be set for each SNES object 3374 3375 Level: intermediate 3376 3377 .keywords: SNES, nonlinear, set, monitor 3378 3379 .seealso: SNESMonitorDefault(), SNESMonitorCancel(), SNESMonitorFunction 3380 @*/ 3381 PetscErrorCode SNESMonitorSet(SNES snes,PetscErrorCode (*f)(SNES,PetscInt,PetscReal,void*),void *mctx,PetscErrorCode (*monitordestroy)(void**)) 3382 { 3383 PetscInt i; 3384 PetscErrorCode ierr; 3385 3386 PetscFunctionBegin; 3387 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3388 if (snes->numbermonitors >= MAXSNESMONITORS) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Too many monitors set"); 3389 for (i=0; i<snes->numbermonitors;i++) { 3390 if (f == snes->monitor[i] && monitordestroy == snes->monitordestroy[i] && mctx == snes->monitorcontext[i]) { 3391 if (monitordestroy) { 3392 ierr = (*monitordestroy)(&mctx);CHKERRQ(ierr); 3393 } 3394 PetscFunctionReturn(0); 3395 } 3396 } 3397 snes->monitor[snes->numbermonitors] = f; 3398 snes->monitordestroy[snes->numbermonitors] = monitordestroy; 3399 snes->monitorcontext[snes->numbermonitors++] = (void*)mctx; 3400 PetscFunctionReturn(0); 3401 } 3402 3403 #undef __FUNCT__ 3404 #define __FUNCT__ "SNESMonitorCancel" 3405 /*@ 3406 SNESMonitorCancel - Clears all the monitor functions for a SNES object. 3407 3408 Logically Collective on SNES 3409 3410 Input Parameters: 3411 . snes - the SNES context 3412 3413 Options Database Key: 3414 . -snes_monitor_cancel - cancels all monitors that have been hardwired 3415 into a code by calls to SNESMonitorSet(), but does not cancel those 3416 set via the options database 3417 3418 Notes: 3419 There is no way to clear one specific monitor from a SNES object. 3420 3421 Level: intermediate 3422 3423 .keywords: SNES, nonlinear, set, monitor 3424 3425 .seealso: SNESMonitorDefault(), SNESMonitorSet() 3426 @*/ 3427 PetscErrorCode SNESMonitorCancel(SNES snes) 3428 { 3429 PetscErrorCode ierr; 3430 PetscInt i; 3431 3432 PetscFunctionBegin; 3433 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3434 for (i=0; i<snes->numbermonitors; i++) { 3435 if (snes->monitordestroy[i]) { 3436 ierr = (*snes->monitordestroy[i])(&snes->monitorcontext[i]);CHKERRQ(ierr); 3437 } 3438 } 3439 snes->numbermonitors = 0; 3440 PetscFunctionReturn(0); 3441 } 3442 3443 /*MC 3444 SNESConvergenceTestFunction - functional form used for testing of convergence of nonlinear solver 3445 3446 Synopsis: 3447 #include <petscsnes.h> 3448 $ PetscErrorCode SNESConvergenceTest(SNES snes,PetscInt it,PetscReal xnorm,PetscReal gnorm,PetscReal f,SNESConvergedReason *reason,void *cctx) 3449 3450 + snes - the SNES context 3451 . it - current iteration (0 is the first and is before any Newton step) 3452 . cctx - [optional] convergence context 3453 . reason - reason for convergence/divergence 3454 . xnorm - 2-norm of current iterate 3455 . gnorm - 2-norm of current step 3456 - f - 2-norm of function 3457 3458 Level: intermediate 3459 3460 .seealso: SNESSetConvergenceTest(), SNESGetConvergenceTest() 3461 M*/ 3462 3463 #undef __FUNCT__ 3464 #define __FUNCT__ "SNESSetConvergenceTest" 3465 /*@C 3466 SNESSetConvergenceTest - Sets the function that is to be used 3467 to test for convergence of the nonlinear iterative solution. 3468 3469 Logically Collective on SNES 3470 3471 Input Parameters: 3472 + snes - the SNES context 3473 . SNESConvergenceTestFunction - routine to test for convergence 3474 . cctx - [optional] context for private data for the convergence routine (may be NULL) 3475 - destroy - [optional] destructor for the context (may be NULL; NULL_FUNCTION in Fortran) 3476 3477 Level: advanced 3478 3479 .keywords: SNES, nonlinear, set, convergence, test 3480 3481 .seealso: SNESConvergedDefault(), SNESConvergedSkip(), SNESConvergenceTestFunction 3482 @*/ 3483 PetscErrorCode SNESSetConvergenceTest(SNES snes,PetscErrorCode (*SNESConvergenceTestFunction)(SNES,PetscInt,PetscReal,PetscReal,PetscReal,SNESConvergedReason*,void*),void *cctx,PetscErrorCode (*destroy)(void*)) 3484 { 3485 PetscErrorCode ierr; 3486 3487 PetscFunctionBegin; 3488 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3489 if (!SNESConvergenceTestFunction) SNESConvergenceTestFunction = SNESConvergedSkip; 3490 if (snes->ops->convergeddestroy) { 3491 ierr = (*snes->ops->convergeddestroy)(snes->cnvP);CHKERRQ(ierr); 3492 } 3493 snes->ops->converged = SNESConvergenceTestFunction; 3494 snes->ops->convergeddestroy = destroy; 3495 snes->cnvP = cctx; 3496 PetscFunctionReturn(0); 3497 } 3498 3499 #undef __FUNCT__ 3500 #define __FUNCT__ "SNESGetConvergedReason" 3501 /*@ 3502 SNESGetConvergedReason - Gets the reason the SNES iteration was stopped. 3503 3504 Not Collective 3505 3506 Input Parameter: 3507 . snes - the SNES context 3508 3509 Output Parameter: 3510 . reason - negative value indicates diverged, positive value converged, see SNESConvergedReason or the 3511 manual pages for the individual convergence tests for complete lists 3512 3513 Level: intermediate 3514 3515 Notes: Can only be called after the call the SNESSolve() is complete. 3516 3517 .keywords: SNES, nonlinear, set, convergence, test 3518 3519 .seealso: SNESSetConvergenceTest(), SNESConvergedReason 3520 @*/ 3521 PetscErrorCode SNESGetConvergedReason(SNES snes,SNESConvergedReason *reason) 3522 { 3523 PetscFunctionBegin; 3524 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3525 PetscValidPointer(reason,2); 3526 *reason = snes->reason; 3527 PetscFunctionReturn(0); 3528 } 3529 3530 #undef __FUNCT__ 3531 #define __FUNCT__ "SNESSetConvergenceHistory" 3532 /*@ 3533 SNESSetConvergenceHistory - Sets the array used to hold the convergence history. 3534 3535 Logically Collective on SNES 3536 3537 Input Parameters: 3538 + snes - iterative context obtained from SNESCreate() 3539 . a - array to hold history, this array will contain the function norms computed at each step 3540 . its - integer array holds the number of linear iterations for each solve. 3541 . na - size of a and its 3542 - reset - PETSC_TRUE indicates each new nonlinear solve resets the history counter to zero, 3543 else it continues storing new values for new nonlinear solves after the old ones 3544 3545 Notes: 3546 If 'a' and 'its' are NULL then space is allocated for the history. If 'na' PETSC_DECIDE or PETSC_DEFAULT then a 3547 default array of length 10000 is allocated. 3548 3549 This routine is useful, e.g., when running a code for purposes 3550 of accurate performance monitoring, when no I/O should be done 3551 during the section of code that is being timed. 3552 3553 Level: intermediate 3554 3555 .keywords: SNES, set, convergence, history 3556 3557 .seealso: SNESGetConvergenceHistory() 3558 3559 @*/ 3560 PetscErrorCode SNESSetConvergenceHistory(SNES snes,PetscReal a[],PetscInt its[],PetscInt na,PetscBool reset) 3561 { 3562 PetscErrorCode ierr; 3563 3564 PetscFunctionBegin; 3565 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3566 if (a) PetscValidScalarPointer(a,2); 3567 if (its) PetscValidIntPointer(its,3); 3568 if (!a) { 3569 if (na == PETSC_DECIDE || na == PETSC_DEFAULT) na = 1000; 3570 ierr = PetscCalloc1(na,&a);CHKERRQ(ierr); 3571 ierr = PetscCalloc1(na,&its);CHKERRQ(ierr); 3572 3573 snes->conv_malloc = PETSC_TRUE; 3574 } 3575 snes->conv_hist = a; 3576 snes->conv_hist_its = its; 3577 snes->conv_hist_max = na; 3578 snes->conv_hist_len = 0; 3579 snes->conv_hist_reset = reset; 3580 PetscFunctionReturn(0); 3581 } 3582 3583 #if defined(PETSC_HAVE_MATLAB_ENGINE) 3584 #include <engine.h> /* MATLAB include file */ 3585 #include <mex.h> /* MATLAB include file */ 3586 3587 #undef __FUNCT__ 3588 #define __FUNCT__ "SNESGetConvergenceHistoryMatlab" 3589 PETSC_EXTERN mxArray *SNESGetConvergenceHistoryMatlab(SNES snes) 3590 { 3591 mxArray *mat; 3592 PetscInt i; 3593 PetscReal *ar; 3594 3595 PetscFunctionBegin; 3596 mat = mxCreateDoubleMatrix(snes->conv_hist_len,1,mxREAL); 3597 ar = (PetscReal*) mxGetData(mat); 3598 for (i=0; i<snes->conv_hist_len; i++) ar[i] = snes->conv_hist[i]; 3599 PetscFunctionReturn(mat); 3600 } 3601 #endif 3602 3603 #undef __FUNCT__ 3604 #define __FUNCT__ "SNESGetConvergenceHistory" 3605 /*@C 3606 SNESGetConvergenceHistory - Gets the array used to hold the convergence history. 3607 3608 Not Collective 3609 3610 Input Parameter: 3611 . snes - iterative context obtained from SNESCreate() 3612 3613 Output Parameters: 3614 . a - array to hold history 3615 . its - integer array holds the number of linear iterations (or 3616 negative if not converged) for each solve. 3617 - na - size of a and its 3618 3619 Notes: 3620 The calling sequence for this routine in Fortran is 3621 $ call SNESGetConvergenceHistory(SNES snes, integer na, integer ierr) 3622 3623 This routine is useful, e.g., when running a code for purposes 3624 of accurate performance monitoring, when no I/O should be done 3625 during the section of code that is being timed. 3626 3627 Level: intermediate 3628 3629 .keywords: SNES, get, convergence, history 3630 3631 .seealso: SNESSetConvergencHistory() 3632 3633 @*/ 3634 PetscErrorCode SNESGetConvergenceHistory(SNES snes,PetscReal *a[],PetscInt *its[],PetscInt *na) 3635 { 3636 PetscFunctionBegin; 3637 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3638 if (a) *a = snes->conv_hist; 3639 if (its) *its = snes->conv_hist_its; 3640 if (na) *na = snes->conv_hist_len; 3641 PetscFunctionReturn(0); 3642 } 3643 3644 #undef __FUNCT__ 3645 #define __FUNCT__ "SNESSetUpdate" 3646 /*@C 3647 SNESSetUpdate - Sets the general-purpose update function called 3648 at the beginning of every iteration of the nonlinear solve. Specifically 3649 it is called just before the Jacobian is "evaluated". 3650 3651 Logically Collective on SNES 3652 3653 Input Parameters: 3654 . snes - The nonlinear solver context 3655 . func - The function 3656 3657 Calling sequence of func: 3658 . func (SNES snes, PetscInt step); 3659 3660 . step - The current step of the iteration 3661 3662 Level: advanced 3663 3664 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() 3665 This is not used by most users. 3666 3667 .keywords: SNES, update 3668 3669 .seealso SNESSetJacobian(), SNESSolve() 3670 @*/ 3671 PetscErrorCode SNESSetUpdate(SNES snes, PetscErrorCode (*func)(SNES, PetscInt)) 3672 { 3673 PetscFunctionBegin; 3674 PetscValidHeaderSpecific(snes, SNES_CLASSID,1); 3675 snes->ops->update = func; 3676 PetscFunctionReturn(0); 3677 } 3678 3679 #undef __FUNCT__ 3680 #define __FUNCT__ "SNESScaleStep_Private" 3681 /* 3682 SNESScaleStep_Private - Scales a step so that its length is less than the 3683 positive parameter delta. 3684 3685 Input Parameters: 3686 + snes - the SNES context 3687 . y - approximate solution of linear system 3688 . fnorm - 2-norm of current function 3689 - delta - trust region size 3690 3691 Output Parameters: 3692 + gpnorm - predicted function norm at the new point, assuming local 3693 linearization. The value is zero if the step lies within the trust 3694 region, and exceeds zero otherwise. 3695 - ynorm - 2-norm of the step 3696 3697 Note: 3698 For non-trust region methods such as SNESNEWTONLS, the parameter delta 3699 is set to be the maximum allowable step size. 3700 3701 .keywords: SNES, nonlinear, scale, step 3702 */ 3703 PetscErrorCode SNESScaleStep_Private(SNES snes,Vec y,PetscReal *fnorm,PetscReal *delta,PetscReal *gpnorm,PetscReal *ynorm) 3704 { 3705 PetscReal nrm; 3706 PetscScalar cnorm; 3707 PetscErrorCode ierr; 3708 3709 PetscFunctionBegin; 3710 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3711 PetscValidHeaderSpecific(y,VEC_CLASSID,2); 3712 PetscCheckSameComm(snes,1,y,2); 3713 3714 ierr = VecNorm(y,NORM_2,&nrm);CHKERRQ(ierr); 3715 if (nrm > *delta) { 3716 nrm = *delta/nrm; 3717 *gpnorm = (1.0 - nrm)*(*fnorm); 3718 cnorm = nrm; 3719 ierr = VecScale(y,cnorm);CHKERRQ(ierr); 3720 *ynorm = *delta; 3721 } else { 3722 *gpnorm = 0.0; 3723 *ynorm = nrm; 3724 } 3725 PetscFunctionReturn(0); 3726 } 3727 3728 #undef __FUNCT__ 3729 #define __FUNCT__ "SNESReasonView" 3730 /*@ 3731 SNESReasonView - Displays the reason a SNES solve converged or diverged to a viewer 3732 3733 Collective on SNES 3734 3735 Parameter: 3736 + snes - iterative context obtained from SNESCreate() 3737 - viewer - the viewer to display the reason 3738 3739 3740 Options Database Keys: 3741 . -snes_converged_reason - print reason for converged or diverged, also prints number of iterations 3742 3743 Level: beginner 3744 3745 .keywords: SNES, solve, linear system 3746 3747 .seealso: SNESCreate(), SNESSetUp(), SNESDestroy(), SNESSetTolerances(), SNESConvergedDefault() 3748 3749 @*/ 3750 PetscErrorCode SNESReasonView(SNES snes,PetscViewer viewer) 3751 { 3752 PetscErrorCode ierr; 3753 PetscBool isAscii; 3754 3755 PetscFunctionBegin; 3756 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isAscii);CHKERRQ(ierr); 3757 if (isAscii) { 3758 ierr = PetscViewerASCIIAddTab(viewer,((PetscObject)snes)->tablevel);CHKERRQ(ierr); 3759 if (snes->reason > 0) { 3760 if (((PetscObject) snes)->prefix) { 3761 ierr = PetscViewerASCIIPrintf(viewer,"Nonlinear %s solve converged due to %s iterations %D\n",((PetscObject) snes)->prefix,SNESConvergedReasons[snes->reason],snes->iter);CHKERRQ(ierr); 3762 } else { 3763 ierr = PetscViewerASCIIPrintf(viewer,"Nonlinear solve converged due to %s iterations %D\n",SNESConvergedReasons[snes->reason],snes->iter);CHKERRQ(ierr); 3764 } 3765 } else { 3766 if (((PetscObject) snes)->prefix) { 3767 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); 3768 } else { 3769 ierr = PetscViewerASCIIPrintf(viewer,"Nonlinear solve did not converge due to %s iterations %D\n",SNESConvergedReasons[snes->reason],snes->iter);CHKERRQ(ierr); 3770 } 3771 } 3772 ierr = PetscViewerASCIISubtractTab(viewer,((PetscObject)snes)->tablevel);CHKERRQ(ierr); 3773 } 3774 PetscFunctionReturn(0); 3775 } 3776 3777 #undef __FUNCT__ 3778 #define __FUNCT__ "SNESReasonViewFromOptions" 3779 /*@C 3780 SNESReasonViewFromOptions - Processes command line options to determine if/how a SNESReason is to be viewed. 3781 3782 Collective on SNES 3783 3784 Input Parameters: 3785 . snes - the SNES object 3786 3787 Level: intermediate 3788 3789 @*/ 3790 PetscErrorCode SNESReasonViewFromOptions(SNES snes) 3791 { 3792 PetscErrorCode ierr; 3793 PetscViewer viewer; 3794 PetscBool flg; 3795 static PetscBool incall = PETSC_FALSE; 3796 PetscViewerFormat format; 3797 3798 PetscFunctionBegin; 3799 if (incall) PetscFunctionReturn(0); 3800 incall = PETSC_TRUE; 3801 ierr = PetscOptionsGetViewer(PetscObjectComm((PetscObject)snes),((PetscObject)snes)->prefix,"-snes_converged_reason",&viewer,&format,&flg);CHKERRQ(ierr); 3802 if (flg) { 3803 ierr = PetscViewerPushFormat(viewer,format);CHKERRQ(ierr); 3804 ierr = SNESReasonView(snes,viewer);CHKERRQ(ierr); 3805 ierr = PetscViewerPopFormat(viewer);CHKERRQ(ierr); 3806 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 3807 } 3808 incall = PETSC_FALSE; 3809 PetscFunctionReturn(0); 3810 } 3811 3812 #undef __FUNCT__ 3813 #define __FUNCT__ "SNESSolve" 3814 /*@C 3815 SNESSolve - Solves a nonlinear system F(x) = b. 3816 Call SNESSolve() after calling SNESCreate() and optional routines of the form SNESSetXXX(). 3817 3818 Collective on SNES 3819 3820 Input Parameters: 3821 + snes - the SNES context 3822 . b - the constant part of the equation F(x) = b, or NULL to use zero. 3823 - x - the solution vector. 3824 3825 Notes: 3826 The user should initialize the vector,x, with the initial guess 3827 for the nonlinear solve prior to calling SNESSolve. In particular, 3828 to employ an initial guess of zero, the user should explicitly set 3829 this vector to zero by calling VecSet(). 3830 3831 Level: beginner 3832 3833 .keywords: SNES, nonlinear, solve 3834 3835 .seealso: SNESCreate(), SNESDestroy(), SNESSetFunction(), SNESSetJacobian(), SNESSetGridSequence(), SNESGetSolution() 3836 @*/ 3837 PetscErrorCode SNESSolve(SNES snes,Vec b,Vec x) 3838 { 3839 PetscErrorCode ierr; 3840 PetscBool flg; 3841 PetscInt grid; 3842 Vec xcreated = NULL; 3843 DM dm; 3844 3845 PetscFunctionBegin; 3846 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3847 if (x) PetscValidHeaderSpecific(x,VEC_CLASSID,3); 3848 if (x) PetscCheckSameComm(snes,1,x,3); 3849 if (b) PetscValidHeaderSpecific(b,VEC_CLASSID,2); 3850 if (b) PetscCheckSameComm(snes,1,b,2); 3851 3852 if (!x) { 3853 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 3854 ierr = DMCreateGlobalVector(dm,&xcreated);CHKERRQ(ierr); 3855 x = xcreated; 3856 } 3857 ierr = SNESViewFromOptions(snes,NULL,"-snes_view_pre");CHKERRQ(ierr); 3858 3859 for (grid=0; grid<snes->gridsequence; grid++) {ierr = PetscViewerASCIIPushTab(PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)snes)));CHKERRQ(ierr);} 3860 for (grid=0; grid<snes->gridsequence+1; grid++) { 3861 3862 /* set solution vector */ 3863 if (!grid) {ierr = PetscObjectReference((PetscObject)x);CHKERRQ(ierr);} 3864 ierr = VecDestroy(&snes->vec_sol);CHKERRQ(ierr); 3865 snes->vec_sol = x; 3866 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 3867 3868 /* set affine vector if provided */ 3869 if (b) { ierr = PetscObjectReference((PetscObject)b);CHKERRQ(ierr); } 3870 ierr = VecDestroy(&snes->vec_rhs);CHKERRQ(ierr); 3871 snes->vec_rhs = b; 3872 3873 if (snes->vec_func == snes->vec_sol) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_IDN,"Solution vector cannot be function vector"); 3874 if (snes->vec_rhs == snes->vec_sol) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_IDN,"Solution vector cannot be right hand side vector"); 3875 if (!snes->vec_sol_update /* && snes->vec_sol */) { 3876 ierr = VecDuplicate(snes->vec_sol,&snes->vec_sol_update);CHKERRQ(ierr); 3877 ierr = PetscLogObjectParent((PetscObject)snes,(PetscObject)snes->vec_sol_update);CHKERRQ(ierr); 3878 } 3879 ierr = DMShellSetGlobalVector(dm,snes->vec_sol);CHKERRQ(ierr); 3880 ierr = SNESSetUp(snes);CHKERRQ(ierr); 3881 3882 if (!grid) { 3883 if (snes->ops->computeinitialguess) { 3884 ierr = (*snes->ops->computeinitialguess)(snes,snes->vec_sol,snes->initialguessP);CHKERRQ(ierr); 3885 } 3886 } 3887 3888 if (snes->conv_hist_reset) snes->conv_hist_len = 0; 3889 if (snes->counters_reset) {snes->nfuncs = 0; snes->linear_its = 0; snes->numFailures = 0;} 3890 3891 ierr = PetscLogEventBegin(SNES_Solve,snes,0,0,0);CHKERRQ(ierr); 3892 ierr = (*snes->ops->solve)(snes);CHKERRQ(ierr); 3893 ierr = PetscLogEventEnd(SNES_Solve,snes,0,0,0);CHKERRQ(ierr); 3894 if (!snes->reason) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Internal error, solver returned without setting converged reason"); 3895 snes->domainerror = PETSC_FALSE; /* clear the flag if it has been set */ 3896 3897 if (snes->lagjac_persist) snes->jac_iter += snes->iter; 3898 if (snes->lagpre_persist) snes->pre_iter += snes->iter; 3899 3900 flg = PETSC_FALSE; 3901 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_test_local_min",&flg,NULL);CHKERRQ(ierr); 3902 if (flg && !PetscPreLoadingOn) { ierr = SNESTestLocalMin(snes);CHKERRQ(ierr); } 3903 ierr = SNESReasonViewFromOptions(snes);CHKERRQ(ierr); 3904 3905 if (snes->errorifnotconverged && snes->reason < 0) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_NOT_CONVERGED,"SNESSolve has not converged"); 3906 if (grid < snes->gridsequence) { 3907 DM fine; 3908 Vec xnew; 3909 Mat interp; 3910 3911 ierr = DMRefine(snes->dm,PetscObjectComm((PetscObject)snes),&fine);CHKERRQ(ierr); 3912 if (!fine) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_INCOMP,"DMRefine() did not perform any refinement, cannot continue grid sequencing"); 3913 ierr = DMCreateInterpolation(snes->dm,fine,&interp,NULL);CHKERRQ(ierr); 3914 ierr = DMCreateGlobalVector(fine,&xnew);CHKERRQ(ierr); 3915 ierr = MatInterpolate(interp,x,xnew);CHKERRQ(ierr); 3916 ierr = DMInterpolate(snes->dm,interp,fine);CHKERRQ(ierr); 3917 ierr = MatDestroy(&interp);CHKERRQ(ierr); 3918 x = xnew; 3919 3920 ierr = SNESReset(snes);CHKERRQ(ierr); 3921 ierr = SNESSetDM(snes,fine);CHKERRQ(ierr); 3922 ierr = DMDestroy(&fine);CHKERRQ(ierr); 3923 ierr = PetscViewerASCIIPopTab(PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)snes)));CHKERRQ(ierr); 3924 } 3925 } 3926 ierr = SNESViewFromOptions(snes,NULL,"-snes_view");CHKERRQ(ierr); 3927 ierr = VecViewFromOptions(snes->vec_sol,((PetscObject)snes)->prefix,"-snes_view_solution");CHKERRQ(ierr); 3928 3929 ierr = VecDestroy(&xcreated);CHKERRQ(ierr); 3930 ierr = PetscObjectSAWsBlock((PetscObject)snes);CHKERRQ(ierr); 3931 PetscFunctionReturn(0); 3932 } 3933 3934 /* --------- Internal routines for SNES Package --------- */ 3935 3936 #undef __FUNCT__ 3937 #define __FUNCT__ "SNESSetType" 3938 /*@C 3939 SNESSetType - Sets the method for the nonlinear solver. 3940 3941 Collective on SNES 3942 3943 Input Parameters: 3944 + snes - the SNES context 3945 - type - a known method 3946 3947 Options Database Key: 3948 . -snes_type <type> - Sets the method; use -help for a list 3949 of available methods (for instance, newtonls or newtontr) 3950 3951 Notes: 3952 See "petsc/include/petscsnes.h" for available methods (for instance) 3953 + SNESNEWTONLS - Newton's method with line search 3954 (systems of nonlinear equations) 3955 . SNESNEWTONTR - Newton's method with trust region 3956 (systems of nonlinear equations) 3957 3958 Normally, it is best to use the SNESSetFromOptions() command and then 3959 set the SNES solver type from the options database rather than by using 3960 this routine. Using the options database provides the user with 3961 maximum flexibility in evaluating the many nonlinear solvers. 3962 The SNESSetType() routine is provided for those situations where it 3963 is necessary to set the nonlinear solver independently of the command 3964 line or options database. This might be the case, for example, when 3965 the choice of solver changes during the execution of the program, 3966 and the user's application is taking responsibility for choosing the 3967 appropriate method. 3968 3969 Developer Notes: SNESRegister() adds a constructor for a new SNESType to SNESList, SNESSetType() locates 3970 the constructor in that list and calls it to create the spexific object. 3971 3972 Level: intermediate 3973 3974 .keywords: SNES, set, type 3975 3976 .seealso: SNESType, SNESCreate(), SNESDestroy(), SNESGetType(), SNESSetFromOptions() 3977 3978 @*/ 3979 PetscErrorCode SNESSetType(SNES snes,SNESType type) 3980 { 3981 PetscErrorCode ierr,(*r)(SNES); 3982 PetscBool match; 3983 3984 PetscFunctionBegin; 3985 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3986 PetscValidCharPointer(type,2); 3987 3988 ierr = PetscObjectTypeCompare((PetscObject)snes,type,&match);CHKERRQ(ierr); 3989 if (match) PetscFunctionReturn(0); 3990 3991 ierr = PetscFunctionListFind(SNESList,type,&r);CHKERRQ(ierr); 3992 if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unable to find requested SNES type %s",type); 3993 /* Destroy the previous private SNES context */ 3994 if (snes->ops->destroy) { 3995 ierr = (*(snes)->ops->destroy)(snes);CHKERRQ(ierr); 3996 snes->ops->destroy = NULL; 3997 } 3998 /* Reinitialize function pointers in SNESOps structure */ 3999 snes->ops->setup = 0; 4000 snes->ops->solve = 0; 4001 snes->ops->view = 0; 4002 snes->ops->setfromoptions = 0; 4003 snes->ops->destroy = 0; 4004 /* Call the SNESCreate_XXX routine for this particular Nonlinear solver */ 4005 snes->setupcalled = PETSC_FALSE; 4006 4007 ierr = PetscObjectChangeTypeName((PetscObject)snes,type);CHKERRQ(ierr); 4008 ierr = (*r)(snes);CHKERRQ(ierr); 4009 PetscFunctionReturn(0); 4010 } 4011 4012 #undef __FUNCT__ 4013 #define __FUNCT__ "SNESGetType" 4014 /*@C 4015 SNESGetType - Gets the SNES method type and name (as a string). 4016 4017 Not Collective 4018 4019 Input Parameter: 4020 . snes - nonlinear solver context 4021 4022 Output Parameter: 4023 . type - SNES method (a character string) 4024 4025 Level: intermediate 4026 4027 .keywords: SNES, nonlinear, get, type, name 4028 @*/ 4029 PetscErrorCode SNESGetType(SNES snes,SNESType *type) 4030 { 4031 PetscFunctionBegin; 4032 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4033 PetscValidPointer(type,2); 4034 *type = ((PetscObject)snes)->type_name; 4035 PetscFunctionReturn(0); 4036 } 4037 4038 #undef __FUNCT__ 4039 #define __FUNCT__ "SNESSetSolution" 4040 /*@ 4041 SNESSetSolution - Sets the solution vector for use by the SNES routines. 4042 4043 Logically Collective on SNES and Vec 4044 4045 Input Parameters: 4046 + snes - the SNES context obtained from SNESCreate() 4047 - u - the solution vector 4048 4049 Level: beginner 4050 4051 .keywords: SNES, set, solution 4052 @*/ 4053 PetscErrorCode SNESSetSolution(SNES snes, Vec u) 4054 { 4055 DM dm; 4056 PetscErrorCode ierr; 4057 4058 PetscFunctionBegin; 4059 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); 4060 PetscValidHeaderSpecific(u, VEC_CLASSID, 2); 4061 ierr = PetscObjectReference((PetscObject) u);CHKERRQ(ierr); 4062 ierr = VecDestroy(&snes->vec_sol);CHKERRQ(ierr); 4063 4064 snes->vec_sol = u; 4065 4066 ierr = SNESGetDM(snes, &dm);CHKERRQ(ierr); 4067 ierr = DMShellSetGlobalVector(dm, u);CHKERRQ(ierr); 4068 PetscFunctionReturn(0); 4069 } 4070 4071 #undef __FUNCT__ 4072 #define __FUNCT__ "SNESGetSolution" 4073 /*@ 4074 SNESGetSolution - Returns the vector where the approximate solution is 4075 stored. This is the fine grid solution when using SNESSetGridSequence(). 4076 4077 Not Collective, but Vec is parallel if SNES is parallel 4078 4079 Input Parameter: 4080 . snes - the SNES context 4081 4082 Output Parameter: 4083 . x - the solution 4084 4085 Level: intermediate 4086 4087 .keywords: SNES, nonlinear, get, solution 4088 4089 .seealso: SNESGetSolutionUpdate(), SNESGetFunction() 4090 @*/ 4091 PetscErrorCode SNESGetSolution(SNES snes,Vec *x) 4092 { 4093 PetscFunctionBegin; 4094 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4095 PetscValidPointer(x,2); 4096 *x = snes->vec_sol; 4097 PetscFunctionReturn(0); 4098 } 4099 4100 #undef __FUNCT__ 4101 #define __FUNCT__ "SNESGetSolutionUpdate" 4102 /*@ 4103 SNESGetSolutionUpdate - Returns the vector where the solution update is 4104 stored. 4105 4106 Not Collective, but Vec is parallel if SNES is parallel 4107 4108 Input Parameter: 4109 . snes - the SNES context 4110 4111 Output Parameter: 4112 . x - the solution update 4113 4114 Level: advanced 4115 4116 .keywords: SNES, nonlinear, get, solution, update 4117 4118 .seealso: SNESGetSolution(), SNESGetFunction() 4119 @*/ 4120 PetscErrorCode SNESGetSolutionUpdate(SNES snes,Vec *x) 4121 { 4122 PetscFunctionBegin; 4123 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4124 PetscValidPointer(x,2); 4125 *x = snes->vec_sol_update; 4126 PetscFunctionReturn(0); 4127 } 4128 4129 #undef __FUNCT__ 4130 #define __FUNCT__ "SNESGetFunction" 4131 /*@C 4132 SNESGetFunction - Returns the vector where the function is stored. 4133 4134 Not Collective, but Vec is parallel if SNES is parallel. Collective if Vec is requested, but has not been created yet. 4135 4136 Input Parameter: 4137 . snes - the SNES context 4138 4139 Output Parameter: 4140 + r - the vector that is used to store residuals (or NULL if you don't want it) 4141 . f - the function (or NULL if you don't want it); see SNESFunction for calling sequence details 4142 - ctx - the function context (or NULL if you don't want it) 4143 4144 Level: advanced 4145 4146 .keywords: SNES, nonlinear, get, function 4147 4148 .seealso: SNESSetFunction(), SNESGetSolution(), SNESFunction 4149 @*/ 4150 PetscErrorCode SNESGetFunction(SNES snes,Vec *r,PetscErrorCode (**f)(SNES,Vec,Vec,void*),void **ctx) 4151 { 4152 PetscErrorCode ierr; 4153 DM dm; 4154 4155 PetscFunctionBegin; 4156 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4157 if (r) { 4158 if (!snes->vec_func) { 4159 if (snes->vec_rhs) { 4160 ierr = VecDuplicate(snes->vec_rhs,&snes->vec_func);CHKERRQ(ierr); 4161 } else if (snes->vec_sol) { 4162 ierr = VecDuplicate(snes->vec_sol,&snes->vec_func);CHKERRQ(ierr); 4163 } else if (snes->dm) { 4164 ierr = DMCreateGlobalVector(snes->dm,&snes->vec_func);CHKERRQ(ierr); 4165 } 4166 } 4167 *r = snes->vec_func; 4168 } 4169 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 4170 ierr = DMSNESGetFunction(dm,f,ctx);CHKERRQ(ierr); 4171 PetscFunctionReturn(0); 4172 } 4173 4174 /*@C 4175 SNESGetNGS - Returns the NGS function and context. 4176 4177 Input Parameter: 4178 . snes - the SNES context 4179 4180 Output Parameter: 4181 + f - the function (or NULL) see SNESNGSFunction for details 4182 - ctx - the function context (or NULL) 4183 4184 Level: advanced 4185 4186 .keywords: SNES, nonlinear, get, function 4187 4188 .seealso: SNESSetNGS(), SNESGetFunction() 4189 @*/ 4190 4191 #undef __FUNCT__ 4192 #define __FUNCT__ "SNESGetNGS" 4193 PetscErrorCode SNESGetNGS (SNES snes, PetscErrorCode (**f)(SNES, Vec, Vec, void*), void ** ctx) 4194 { 4195 PetscErrorCode ierr; 4196 DM dm; 4197 4198 PetscFunctionBegin; 4199 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4200 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 4201 ierr = DMSNESGetNGS(dm,f,ctx);CHKERRQ(ierr); 4202 PetscFunctionReturn(0); 4203 } 4204 4205 #undef __FUNCT__ 4206 #define __FUNCT__ "SNESSetOptionsPrefix" 4207 /*@C 4208 SNESSetOptionsPrefix - Sets the prefix used for searching for all 4209 SNES options in the database. 4210 4211 Logically Collective on SNES 4212 4213 Input Parameter: 4214 + snes - the SNES context 4215 - prefix - the prefix to prepend to all option names 4216 4217 Notes: 4218 A hyphen (-) must NOT be given at the beginning of the prefix name. 4219 The first character of all runtime options is AUTOMATICALLY the hyphen. 4220 4221 Level: advanced 4222 4223 .keywords: SNES, set, options, prefix, database 4224 4225 .seealso: SNESSetFromOptions() 4226 @*/ 4227 PetscErrorCode SNESSetOptionsPrefix(SNES snes,const char prefix[]) 4228 { 4229 PetscErrorCode ierr; 4230 4231 PetscFunctionBegin; 4232 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4233 ierr = PetscObjectSetOptionsPrefix((PetscObject)snes,prefix);CHKERRQ(ierr); 4234 if (!snes->ksp) {ierr = SNESGetKSP(snes,&snes->ksp);CHKERRQ(ierr);} 4235 if (snes->linesearch) { 4236 ierr = SNESGetLineSearch(snes,&snes->linesearch);CHKERRQ(ierr); 4237 ierr = PetscObjectSetOptionsPrefix((PetscObject)snes->linesearch,prefix);CHKERRQ(ierr); 4238 } 4239 ierr = KSPSetOptionsPrefix(snes->ksp,prefix);CHKERRQ(ierr); 4240 PetscFunctionReturn(0); 4241 } 4242 4243 #undef __FUNCT__ 4244 #define __FUNCT__ "SNESAppendOptionsPrefix" 4245 /*@C 4246 SNESAppendOptionsPrefix - Appends to the prefix used for searching for all 4247 SNES options in the database. 4248 4249 Logically Collective on SNES 4250 4251 Input Parameters: 4252 + snes - the SNES context 4253 - prefix - the prefix to prepend to all option names 4254 4255 Notes: 4256 A hyphen (-) must NOT be given at the beginning of the prefix name. 4257 The first character of all runtime options is AUTOMATICALLY the hyphen. 4258 4259 Level: advanced 4260 4261 .keywords: SNES, append, options, prefix, database 4262 4263 .seealso: SNESGetOptionsPrefix() 4264 @*/ 4265 PetscErrorCode SNESAppendOptionsPrefix(SNES snes,const char prefix[]) 4266 { 4267 PetscErrorCode ierr; 4268 4269 PetscFunctionBegin; 4270 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4271 ierr = PetscObjectAppendOptionsPrefix((PetscObject)snes,prefix);CHKERRQ(ierr); 4272 if (!snes->ksp) {ierr = SNESGetKSP(snes,&snes->ksp);CHKERRQ(ierr);} 4273 if (snes->linesearch) { 4274 ierr = SNESGetLineSearch(snes,&snes->linesearch);CHKERRQ(ierr); 4275 ierr = PetscObjectAppendOptionsPrefix((PetscObject)snes->linesearch,prefix);CHKERRQ(ierr); 4276 } 4277 ierr = KSPAppendOptionsPrefix(snes->ksp,prefix);CHKERRQ(ierr); 4278 PetscFunctionReturn(0); 4279 } 4280 4281 #undef __FUNCT__ 4282 #define __FUNCT__ "SNESGetOptionsPrefix" 4283 /*@C 4284 SNESGetOptionsPrefix - Sets the prefix used for searching for all 4285 SNES options in the database. 4286 4287 Not Collective 4288 4289 Input Parameter: 4290 . snes - the SNES context 4291 4292 Output Parameter: 4293 . prefix - pointer to the prefix string used 4294 4295 Notes: On the fortran side, the user should pass in a string 'prefix' of 4296 sufficient length to hold the prefix. 4297 4298 Level: advanced 4299 4300 .keywords: SNES, get, options, prefix, database 4301 4302 .seealso: SNESAppendOptionsPrefix() 4303 @*/ 4304 PetscErrorCode SNESGetOptionsPrefix(SNES snes,const char *prefix[]) 4305 { 4306 PetscErrorCode ierr; 4307 4308 PetscFunctionBegin; 4309 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4310 ierr = PetscObjectGetOptionsPrefix((PetscObject)snes,prefix);CHKERRQ(ierr); 4311 PetscFunctionReturn(0); 4312 } 4313 4314 4315 #undef __FUNCT__ 4316 #define __FUNCT__ "SNESRegister" 4317 /*@C 4318 SNESRegister - Adds a method to the nonlinear solver package. 4319 4320 Not collective 4321 4322 Input Parameters: 4323 + name_solver - name of a new user-defined solver 4324 - routine_create - routine to create method context 4325 4326 Notes: 4327 SNESRegister() may be called multiple times to add several user-defined solvers. 4328 4329 Sample usage: 4330 .vb 4331 SNESRegister("my_solver",MySolverCreate); 4332 .ve 4333 4334 Then, your solver can be chosen with the procedural interface via 4335 $ SNESSetType(snes,"my_solver") 4336 or at runtime via the option 4337 $ -snes_type my_solver 4338 4339 Level: advanced 4340 4341 Note: If your function is not being put into a shared library then use SNESRegister() instead 4342 4343 .keywords: SNES, nonlinear, register 4344 4345 .seealso: SNESRegisterAll(), SNESRegisterDestroy() 4346 4347 Level: advanced 4348 @*/ 4349 PetscErrorCode SNESRegister(const char sname[],PetscErrorCode (*function)(SNES)) 4350 { 4351 PetscErrorCode ierr; 4352 4353 PetscFunctionBegin; 4354 ierr = PetscFunctionListAdd(&SNESList,sname,function);CHKERRQ(ierr); 4355 PetscFunctionReturn(0); 4356 } 4357 4358 #undef __FUNCT__ 4359 #define __FUNCT__ "SNESTestLocalMin" 4360 PetscErrorCode SNESTestLocalMin(SNES snes) 4361 { 4362 PetscErrorCode ierr; 4363 PetscInt N,i,j; 4364 Vec u,uh,fh; 4365 PetscScalar value; 4366 PetscReal norm; 4367 4368 PetscFunctionBegin; 4369 ierr = SNESGetSolution(snes,&u);CHKERRQ(ierr); 4370 ierr = VecDuplicate(u,&uh);CHKERRQ(ierr); 4371 ierr = VecDuplicate(u,&fh);CHKERRQ(ierr); 4372 4373 /* currently only works for sequential */ 4374 ierr = PetscPrintf(PETSC_COMM_WORLD,"Testing FormFunction() for local min\n");CHKERRQ(ierr); 4375 ierr = VecGetSize(u,&N);CHKERRQ(ierr); 4376 for (i=0; i<N; i++) { 4377 ierr = VecCopy(u,uh);CHKERRQ(ierr); 4378 ierr = PetscPrintf(PETSC_COMM_WORLD,"i = %D\n",i);CHKERRQ(ierr); 4379 for (j=-10; j<11; j++) { 4380 value = PetscSign(j)*PetscExpReal(PetscAbs(j)-10.0); 4381 ierr = VecSetValue(uh,i,value,ADD_VALUES);CHKERRQ(ierr); 4382 ierr = SNESComputeFunction(snes,uh,fh);CHKERRQ(ierr); 4383 ierr = VecNorm(fh,NORM_2,&norm);CHKERRQ(ierr); 4384 ierr = PetscPrintf(PETSC_COMM_WORLD," j norm %D %18.16e\n",j,norm);CHKERRQ(ierr); 4385 value = -value; 4386 ierr = VecSetValue(uh,i,value,ADD_VALUES);CHKERRQ(ierr); 4387 } 4388 } 4389 ierr = VecDestroy(&uh);CHKERRQ(ierr); 4390 ierr = VecDestroy(&fh);CHKERRQ(ierr); 4391 PetscFunctionReturn(0); 4392 } 4393 4394 #undef __FUNCT__ 4395 #define __FUNCT__ "SNESKSPSetUseEW" 4396 /*@ 4397 SNESKSPSetUseEW - Sets SNES use Eisenstat-Walker method for 4398 computing relative tolerance for linear solvers within an inexact 4399 Newton method. 4400 4401 Logically Collective on SNES 4402 4403 Input Parameters: 4404 + snes - SNES context 4405 - flag - PETSC_TRUE or PETSC_FALSE 4406 4407 Options Database: 4408 + -snes_ksp_ew - use Eisenstat-Walker method for determining linear system convergence 4409 . -snes_ksp_ew_version ver - version of Eisenstat-Walker method 4410 . -snes_ksp_ew_rtol0 <rtol0> - Sets rtol0 4411 . -snes_ksp_ew_rtolmax <rtolmax> - Sets rtolmax 4412 . -snes_ksp_ew_gamma <gamma> - Sets gamma 4413 . -snes_ksp_ew_alpha <alpha> - Sets alpha 4414 . -snes_ksp_ew_alpha2 <alpha2> - Sets alpha2 4415 - -snes_ksp_ew_threshold <threshold> - Sets threshold 4416 4417 Notes: 4418 Currently, the default is to use a constant relative tolerance for 4419 the inner linear solvers. Alternatively, one can use the 4420 Eisenstat-Walker method, where the relative convergence tolerance 4421 is reset at each Newton iteration according progress of the nonlinear 4422 solver. 4423 4424 Level: advanced 4425 4426 Reference: 4427 S. C. Eisenstat and H. F. Walker, "Choosing the forcing terms in an 4428 inexact Newton method", SISC 17 (1), pp.16-32, 1996. 4429 4430 .keywords: SNES, KSP, Eisenstat, Walker, convergence, test, inexact, Newton 4431 4432 .seealso: SNESKSPGetUseEW(), SNESKSPGetParametersEW(), SNESKSPSetParametersEW() 4433 @*/ 4434 PetscErrorCode SNESKSPSetUseEW(SNES snes,PetscBool flag) 4435 { 4436 PetscFunctionBegin; 4437 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4438 PetscValidLogicalCollectiveBool(snes,flag,2); 4439 snes->ksp_ewconv = flag; 4440 PetscFunctionReturn(0); 4441 } 4442 4443 #undef __FUNCT__ 4444 #define __FUNCT__ "SNESKSPGetUseEW" 4445 /*@ 4446 SNESKSPGetUseEW - Gets if SNES is using Eisenstat-Walker method 4447 for computing relative tolerance for linear solvers within an 4448 inexact Newton method. 4449 4450 Not Collective 4451 4452 Input Parameter: 4453 . snes - SNES context 4454 4455 Output Parameter: 4456 . flag - PETSC_TRUE or PETSC_FALSE 4457 4458 Notes: 4459 Currently, the default is to use a constant relative tolerance for 4460 the inner linear solvers. Alternatively, one can use the 4461 Eisenstat-Walker method, where the relative convergence tolerance 4462 is reset at each Newton iteration according progress of the nonlinear 4463 solver. 4464 4465 Level: advanced 4466 4467 Reference: 4468 S. C. Eisenstat and H. F. Walker, "Choosing the forcing terms in an 4469 inexact Newton method", SISC 17 (1), pp.16-32, 1996. 4470 4471 .keywords: SNES, KSP, Eisenstat, Walker, convergence, test, inexact, Newton 4472 4473 .seealso: SNESKSPSetUseEW(), SNESKSPGetParametersEW(), SNESKSPSetParametersEW() 4474 @*/ 4475 PetscErrorCode SNESKSPGetUseEW(SNES snes, PetscBool *flag) 4476 { 4477 PetscFunctionBegin; 4478 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4479 PetscValidPointer(flag,2); 4480 *flag = snes->ksp_ewconv; 4481 PetscFunctionReturn(0); 4482 } 4483 4484 #undef __FUNCT__ 4485 #define __FUNCT__ "SNESKSPSetParametersEW" 4486 /*@ 4487 SNESKSPSetParametersEW - Sets parameters for Eisenstat-Walker 4488 convergence criteria for the linear solvers within an inexact 4489 Newton method. 4490 4491 Logically Collective on SNES 4492 4493 Input Parameters: 4494 + snes - SNES context 4495 . version - version 1, 2 (default is 2) or 3 4496 . rtol_0 - initial relative tolerance (0 <= rtol_0 < 1) 4497 . rtol_max - maximum relative tolerance (0 <= rtol_max < 1) 4498 . gamma - multiplicative factor for version 2 rtol computation 4499 (0 <= gamma2 <= 1) 4500 . alpha - power for version 2 rtol computation (1 < alpha <= 2) 4501 . alpha2 - power for safeguard 4502 - threshold - threshold for imposing safeguard (0 < threshold < 1) 4503 4504 Note: 4505 Version 3 was contributed by Luis Chacon, June 2006. 4506 4507 Use PETSC_DEFAULT to retain the default for any of the parameters. 4508 4509 Level: advanced 4510 4511 Reference: 4512 S. C. Eisenstat and H. F. Walker, "Choosing the forcing terms in an 4513 inexact Newton method", Utah State University Math. Stat. Dept. Res. 4514 Report 6/94/75, June, 1994, to appear in SIAM J. Sci. Comput. 4515 4516 .keywords: SNES, KSP, Eisenstat, Walker, set, parameters 4517 4518 .seealso: SNESKSPSetUseEW(), SNESKSPGetUseEW(), SNESKSPGetParametersEW() 4519 @*/ 4520 PetscErrorCode SNESKSPSetParametersEW(SNES snes,PetscInt version,PetscReal rtol_0,PetscReal rtol_max,PetscReal gamma,PetscReal alpha,PetscReal alpha2,PetscReal threshold) 4521 { 4522 SNESKSPEW *kctx; 4523 4524 PetscFunctionBegin; 4525 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4526 kctx = (SNESKSPEW*)snes->kspconvctx; 4527 if (!kctx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"No Eisenstat-Walker context existing"); 4528 PetscValidLogicalCollectiveInt(snes,version,2); 4529 PetscValidLogicalCollectiveReal(snes,rtol_0,3); 4530 PetscValidLogicalCollectiveReal(snes,rtol_max,4); 4531 PetscValidLogicalCollectiveReal(snes,gamma,5); 4532 PetscValidLogicalCollectiveReal(snes,alpha,6); 4533 PetscValidLogicalCollectiveReal(snes,alpha2,7); 4534 PetscValidLogicalCollectiveReal(snes,threshold,8); 4535 4536 if (version != PETSC_DEFAULT) kctx->version = version; 4537 if (rtol_0 != PETSC_DEFAULT) kctx->rtol_0 = rtol_0; 4538 if (rtol_max != PETSC_DEFAULT) kctx->rtol_max = rtol_max; 4539 if (gamma != PETSC_DEFAULT) kctx->gamma = gamma; 4540 if (alpha != PETSC_DEFAULT) kctx->alpha = alpha; 4541 if (alpha2 != PETSC_DEFAULT) kctx->alpha2 = alpha2; 4542 if (threshold != PETSC_DEFAULT) kctx->threshold = threshold; 4543 4544 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); 4545 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); 4546 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); 4547 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); 4548 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); 4549 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); 4550 PetscFunctionReturn(0); 4551 } 4552 4553 #undef __FUNCT__ 4554 #define __FUNCT__ "SNESKSPGetParametersEW" 4555 /*@ 4556 SNESKSPGetParametersEW - Gets parameters for Eisenstat-Walker 4557 convergence criteria for the linear solvers within an inexact 4558 Newton method. 4559 4560 Not Collective 4561 4562 Input Parameters: 4563 snes - SNES context 4564 4565 Output Parameters: 4566 + version - version 1, 2 (default is 2) or 3 4567 . rtol_0 - initial relative tolerance (0 <= rtol_0 < 1) 4568 . rtol_max - maximum relative tolerance (0 <= rtol_max < 1) 4569 . gamma - multiplicative factor for version 2 rtol computation (0 <= gamma2 <= 1) 4570 . alpha - power for version 2 rtol computation (1 < alpha <= 2) 4571 . alpha2 - power for safeguard 4572 - threshold - threshold for imposing safeguard (0 < threshold < 1) 4573 4574 Level: advanced 4575 4576 .keywords: SNES, KSP, Eisenstat, Walker, get, parameters 4577 4578 .seealso: SNESKSPSetUseEW(), SNESKSPGetUseEW(), SNESKSPSetParametersEW() 4579 @*/ 4580 PetscErrorCode SNESKSPGetParametersEW(SNES snes,PetscInt *version,PetscReal *rtol_0,PetscReal *rtol_max,PetscReal *gamma,PetscReal *alpha,PetscReal *alpha2,PetscReal *threshold) 4581 { 4582 SNESKSPEW *kctx; 4583 4584 PetscFunctionBegin; 4585 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4586 kctx = (SNESKSPEW*)snes->kspconvctx; 4587 if (!kctx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"No Eisenstat-Walker context existing"); 4588 if (version) *version = kctx->version; 4589 if (rtol_0) *rtol_0 = kctx->rtol_0; 4590 if (rtol_max) *rtol_max = kctx->rtol_max; 4591 if (gamma) *gamma = kctx->gamma; 4592 if (alpha) *alpha = kctx->alpha; 4593 if (alpha2) *alpha2 = kctx->alpha2; 4594 if (threshold) *threshold = kctx->threshold; 4595 PetscFunctionReturn(0); 4596 } 4597 4598 #undef __FUNCT__ 4599 #define __FUNCT__ "KSPPreSolve_SNESEW" 4600 PetscErrorCode KSPPreSolve_SNESEW(KSP ksp, Vec b, Vec x, SNES snes) 4601 { 4602 PetscErrorCode ierr; 4603 SNESKSPEW *kctx = (SNESKSPEW*)snes->kspconvctx; 4604 PetscReal rtol = PETSC_DEFAULT,stol; 4605 4606 PetscFunctionBegin; 4607 if (!snes->ksp_ewconv) PetscFunctionReturn(0); 4608 if (!snes->iter) { 4609 rtol = kctx->rtol_0; /* first time in, so use the original user rtol */ 4610 ierr = VecNorm(snes->vec_func,NORM_2,&kctx->norm_first);CHKERRQ(ierr); 4611 } 4612 else { 4613 if (kctx->version == 1) { 4614 rtol = (snes->norm - kctx->lresid_last)/kctx->norm_last; 4615 if (rtol < 0.0) rtol = -rtol; 4616 stol = PetscPowReal(kctx->rtol_last,kctx->alpha2); 4617 if (stol > kctx->threshold) rtol = PetscMax(rtol,stol); 4618 } else if (kctx->version == 2) { 4619 rtol = kctx->gamma * PetscPowReal(snes->norm/kctx->norm_last,kctx->alpha); 4620 stol = kctx->gamma * PetscPowReal(kctx->rtol_last,kctx->alpha); 4621 if (stol > kctx->threshold) rtol = PetscMax(rtol,stol); 4622 } else if (kctx->version == 3) { /* contributed by Luis Chacon, June 2006. */ 4623 rtol = kctx->gamma * PetscPowReal(snes->norm/kctx->norm_last,kctx->alpha); 4624 /* safeguard: avoid sharp decrease of rtol */ 4625 stol = kctx->gamma*PetscPowReal(kctx->rtol_last,kctx->alpha); 4626 stol = PetscMax(rtol,stol); 4627 rtol = PetscMin(kctx->rtol_0,stol); 4628 /* safeguard: avoid oversolving */ 4629 stol = kctx->gamma*(kctx->norm_first*snes->rtol)/snes->norm; 4630 stol = PetscMax(rtol,stol); 4631 rtol = PetscMin(kctx->rtol_0,stol); 4632 } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Only versions 1, 2 or 3 are supported: %D",kctx->version); 4633 } 4634 /* safeguard: avoid rtol greater than one */ 4635 rtol = PetscMin(rtol,kctx->rtol_max); 4636 ierr = KSPSetTolerances(ksp,rtol,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);CHKERRQ(ierr); 4637 ierr = PetscInfo3(snes,"iter %D, Eisenstat-Walker (version %D) KSP rtol=%g\n",snes->iter,kctx->version,(double)rtol);CHKERRQ(ierr); 4638 PetscFunctionReturn(0); 4639 } 4640 4641 #undef __FUNCT__ 4642 #define __FUNCT__ "KSPPostSolve_SNESEW" 4643 PetscErrorCode KSPPostSolve_SNESEW(KSP ksp, Vec b, Vec x, SNES snes) 4644 { 4645 PetscErrorCode ierr; 4646 SNESKSPEW *kctx = (SNESKSPEW*)snes->kspconvctx; 4647 PCSide pcside; 4648 Vec lres; 4649 4650 PetscFunctionBegin; 4651 if (!snes->ksp_ewconv) PetscFunctionReturn(0); 4652 ierr = KSPGetTolerances(ksp,&kctx->rtol_last,0,0,0);CHKERRQ(ierr); 4653 kctx->norm_last = snes->norm; 4654 if (kctx->version == 1) { 4655 PC pc; 4656 PetscBool isNone; 4657 4658 ierr = KSPGetPC(ksp, &pc);CHKERRQ(ierr); 4659 ierr = PetscObjectTypeCompare((PetscObject) pc, PCNONE, &isNone);CHKERRQ(ierr); 4660 ierr = KSPGetPCSide(ksp,&pcside);CHKERRQ(ierr); 4661 if (pcside == PC_RIGHT || isNone) { /* XXX Should we also test KSP_UNPRECONDITIONED_NORM ? */ 4662 /* KSP residual is true linear residual */ 4663 ierr = KSPGetResidualNorm(ksp,&kctx->lresid_last);CHKERRQ(ierr); 4664 } else { 4665 /* KSP residual is preconditioned residual */ 4666 /* compute true linear residual norm */ 4667 ierr = VecDuplicate(b,&lres);CHKERRQ(ierr); 4668 ierr = MatMult(snes->jacobian,x,lres);CHKERRQ(ierr); 4669 ierr = VecAYPX(lres,-1.0,b);CHKERRQ(ierr); 4670 ierr = VecNorm(lres,NORM_2,&kctx->lresid_last);CHKERRQ(ierr); 4671 ierr = VecDestroy(&lres);CHKERRQ(ierr); 4672 } 4673 } 4674 PetscFunctionReturn(0); 4675 } 4676 4677 #undef __FUNCT__ 4678 #define __FUNCT__ "SNESGetKSP" 4679 /*@ 4680 SNESGetKSP - Returns the KSP context for a SNES solver. 4681 4682 Not Collective, but if SNES object is parallel, then KSP object is parallel 4683 4684 Input Parameter: 4685 . snes - the SNES context 4686 4687 Output Parameter: 4688 . ksp - the KSP context 4689 4690 Notes: 4691 The user can then directly manipulate the KSP context to set various 4692 options, etc. Likewise, the user can then extract and manipulate the 4693 PC contexts as well. 4694 4695 Level: beginner 4696 4697 .keywords: SNES, nonlinear, get, KSP, context 4698 4699 .seealso: KSPGetPC(), SNESCreate(), KSPCreate(), SNESSetKSP() 4700 @*/ 4701 PetscErrorCode SNESGetKSP(SNES snes,KSP *ksp) 4702 { 4703 PetscErrorCode ierr; 4704 4705 PetscFunctionBegin; 4706 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4707 PetscValidPointer(ksp,2); 4708 4709 if (!snes->ksp) { 4710 PetscBool monitor = PETSC_FALSE; 4711 4712 ierr = KSPCreate(PetscObjectComm((PetscObject)snes),&snes->ksp);CHKERRQ(ierr); 4713 ierr = PetscObjectIncrementTabLevel((PetscObject)snes->ksp,(PetscObject)snes,1);CHKERRQ(ierr); 4714 ierr = PetscLogObjectParent((PetscObject)snes,(PetscObject)snes->ksp);CHKERRQ(ierr); 4715 4716 ierr = KSPSetPreSolve(snes->ksp,(PetscErrorCode (*)(KSP,Vec,Vec,void*))KSPPreSolve_SNESEW,snes);CHKERRQ(ierr); 4717 ierr = KSPSetPostSolve(snes->ksp,(PetscErrorCode (*)(KSP,Vec,Vec,void*))KSPPostSolve_SNESEW,snes);CHKERRQ(ierr); 4718 4719 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-ksp_monitor_snes",&monitor,NULL);CHKERRQ(ierr); 4720 if (monitor) { 4721 ierr = KSPMonitorSet(snes->ksp,KSPMonitorSNES,snes,NULL);CHKERRQ(ierr); 4722 } 4723 monitor = PETSC_FALSE; 4724 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-ksp_monitor_snes_lg",&monitor,NULL);CHKERRQ(ierr); 4725 if (monitor) { 4726 PetscObject *objs; 4727 ierr = KSPMonitorSNESLGResidualNormCreate(0,0,PETSC_DECIDE,PETSC_DECIDE,600,600,&objs);CHKERRQ(ierr); 4728 objs[0] = (PetscObject) snes; 4729 ierr = KSPMonitorSet(snes->ksp,(PetscErrorCode (*)(KSP,PetscInt,PetscReal,void*))KSPMonitorSNESLGResidualNorm,objs,(PetscErrorCode (*)(void**))KSPMonitorSNESLGResidualNormDestroy);CHKERRQ(ierr); 4730 } 4731 } 4732 *ksp = snes->ksp; 4733 PetscFunctionReturn(0); 4734 } 4735 4736 4737 #include <petsc/private/dmimpl.h> 4738 #undef __FUNCT__ 4739 #define __FUNCT__ "SNESSetDM" 4740 /*@ 4741 SNESSetDM - Sets the DM that may be used by some preconditioners 4742 4743 Logically Collective on SNES 4744 4745 Input Parameters: 4746 + snes - the preconditioner context 4747 - dm - the dm 4748 4749 Level: intermediate 4750 4751 .seealso: SNESGetDM(), KSPSetDM(), KSPGetDM() 4752 @*/ 4753 PetscErrorCode SNESSetDM(SNES snes,DM dm) 4754 { 4755 PetscErrorCode ierr; 4756 KSP ksp; 4757 DMSNES sdm; 4758 4759 PetscFunctionBegin; 4760 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4761 if (dm) {ierr = PetscObjectReference((PetscObject)dm);CHKERRQ(ierr);} 4762 if (snes->dm) { /* Move the DMSNES context over to the new DM unless the new DM already has one */ 4763 if (snes->dm->dmsnes && snes->dmAuto && !dm->dmsnes) { 4764 ierr = DMCopyDMSNES(snes->dm,dm);CHKERRQ(ierr); 4765 ierr = DMGetDMSNES(snes->dm,&sdm);CHKERRQ(ierr); 4766 if (sdm->originaldm == snes->dm) sdm->originaldm = dm; /* Grant write privileges to the replacement DM */ 4767 } 4768 ierr = DMDestroy(&snes->dm);CHKERRQ(ierr); 4769 } 4770 snes->dm = dm; 4771 snes->dmAuto = PETSC_FALSE; 4772 4773 ierr = SNESGetKSP(snes,&ksp);CHKERRQ(ierr); 4774 ierr = KSPSetDM(ksp,dm);CHKERRQ(ierr); 4775 ierr = KSPSetDMActive(ksp,PETSC_FALSE);CHKERRQ(ierr); 4776 if (snes->pc) { 4777 ierr = SNESSetDM(snes->pc, snes->dm);CHKERRQ(ierr); 4778 ierr = SNESSetNPCSide(snes,snes->pcside);CHKERRQ(ierr); 4779 } 4780 PetscFunctionReturn(0); 4781 } 4782 4783 #undef __FUNCT__ 4784 #define __FUNCT__ "SNESGetDM" 4785 /*@ 4786 SNESGetDM - Gets the DM that may be used by some preconditioners 4787 4788 Not Collective but DM obtained is parallel on SNES 4789 4790 Input Parameter: 4791 . snes - the preconditioner context 4792 4793 Output Parameter: 4794 . dm - the dm 4795 4796 Level: intermediate 4797 4798 .seealso: SNESSetDM(), KSPSetDM(), KSPGetDM() 4799 @*/ 4800 PetscErrorCode SNESGetDM(SNES snes,DM *dm) 4801 { 4802 PetscErrorCode ierr; 4803 4804 PetscFunctionBegin; 4805 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4806 if (!snes->dm) { 4807 ierr = DMShellCreate(PetscObjectComm((PetscObject)snes),&snes->dm);CHKERRQ(ierr); 4808 snes->dmAuto = PETSC_TRUE; 4809 } 4810 *dm = snes->dm; 4811 PetscFunctionReturn(0); 4812 } 4813 4814 #undef __FUNCT__ 4815 #define __FUNCT__ "SNESSetNPC" 4816 /*@ 4817 SNESSetNPC - Sets the nonlinear preconditioner to be used. 4818 4819 Collective on SNES 4820 4821 Input Parameters: 4822 + snes - iterative context obtained from SNESCreate() 4823 - pc - the preconditioner object 4824 4825 Notes: 4826 Use SNESGetNPC() to retrieve the preconditioner context (for example, 4827 to configure it using the API). 4828 4829 Level: developer 4830 4831 .keywords: SNES, set, precondition 4832 .seealso: SNESGetNPC(), SNESHasNPC() 4833 @*/ 4834 PetscErrorCode SNESSetNPC(SNES snes, SNES pc) 4835 { 4836 PetscErrorCode ierr; 4837 4838 PetscFunctionBegin; 4839 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); 4840 PetscValidHeaderSpecific(pc, SNES_CLASSID, 2); 4841 PetscCheckSameComm(snes, 1, pc, 2); 4842 ierr = PetscObjectReference((PetscObject) pc);CHKERRQ(ierr); 4843 ierr = SNESDestroy(&snes->pc);CHKERRQ(ierr); 4844 snes->pc = pc; 4845 ierr = PetscLogObjectParent((PetscObject)snes, (PetscObject)snes->pc);CHKERRQ(ierr); 4846 PetscFunctionReturn(0); 4847 } 4848 4849 #undef __FUNCT__ 4850 #define __FUNCT__ "SNESGetNPC" 4851 /*@ 4852 SNESGetNPC - Creates a nonlinear preconditioning solver (SNES) to be used to precondition the nonlinear solver. 4853 4854 Not Collective 4855 4856 Input Parameter: 4857 . snes - iterative context obtained from SNESCreate() 4858 4859 Output Parameter: 4860 . pc - preconditioner context 4861 4862 Notes: If a SNES was previously set with SNESSetNPC() then that SNES is returned. 4863 4864 Level: developer 4865 4866 .keywords: SNES, get, preconditioner 4867 .seealso: SNESSetNPC(), SNESHasNPC() 4868 @*/ 4869 PetscErrorCode SNESGetNPC(SNES snes, SNES *pc) 4870 { 4871 PetscErrorCode ierr; 4872 const char *optionsprefix; 4873 4874 PetscFunctionBegin; 4875 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); 4876 PetscValidPointer(pc, 2); 4877 if (!snes->pc) { 4878 ierr = SNESCreate(PetscObjectComm((PetscObject)snes),&snes->pc);CHKERRQ(ierr); 4879 ierr = PetscObjectIncrementTabLevel((PetscObject)snes->pc,(PetscObject)snes,1);CHKERRQ(ierr); 4880 ierr = PetscLogObjectParent((PetscObject)snes,(PetscObject)snes->pc);CHKERRQ(ierr); 4881 ierr = SNESGetOptionsPrefix(snes,&optionsprefix);CHKERRQ(ierr); 4882 ierr = SNESSetOptionsPrefix(snes->pc,optionsprefix);CHKERRQ(ierr); 4883 ierr = SNESAppendOptionsPrefix(snes->pc,"npc_");CHKERRQ(ierr); 4884 ierr = SNESSetCountersReset(snes->pc,PETSC_FALSE);CHKERRQ(ierr); 4885 } 4886 *pc = snes->pc; 4887 PetscFunctionReturn(0); 4888 } 4889 4890 #undef __FUNCT__ 4891 #define __FUNCT__ "SNESHasNPC" 4892 /*@ 4893 SNESHasNPC - Returns whether a nonlinear preconditioner exists 4894 4895 Not Collective 4896 4897 Input Parameter: 4898 . snes - iterative context obtained from SNESCreate() 4899 4900 Output Parameter: 4901 . has_npc - whether the SNES has an NPC or not 4902 4903 Level: developer 4904 4905 .keywords: SNES, has, preconditioner 4906 .seealso: SNESSetNPC(), SNESGetNPC() 4907 @*/ 4908 PetscErrorCode SNESHasNPC(SNES snes, PetscBool *has_npc) 4909 { 4910 PetscFunctionBegin; 4911 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); 4912 *has_npc = (PetscBool) (snes->pc != NULL); 4913 PetscFunctionReturn(0); 4914 } 4915 4916 #undef __FUNCT__ 4917 #define __FUNCT__ "SNESSetNPCSide" 4918 /*@ 4919 SNESSetNPCSide - Sets the preconditioning side. 4920 4921 Logically Collective on SNES 4922 4923 Input Parameter: 4924 . snes - iterative context obtained from SNESCreate() 4925 4926 Output Parameter: 4927 . side - the preconditioning side, where side is one of 4928 .vb 4929 PC_LEFT - left preconditioning (default) 4930 PC_RIGHT - right preconditioning 4931 .ve 4932 4933 Options Database Keys: 4934 . -snes_pc_side <right,left> 4935 4936 Level: intermediate 4937 4938 .keywords: SNES, set, right, left, side, preconditioner, flag 4939 4940 .seealso: SNESGetNPCSide(), KSPSetPCSide() 4941 @*/ 4942 PetscErrorCode SNESSetNPCSide(SNES snes,PCSide side) 4943 { 4944 PetscFunctionBegin; 4945 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4946 PetscValidLogicalCollectiveEnum(snes,side,2); 4947 snes->pcside = side; 4948 PetscFunctionReturn(0); 4949 } 4950 4951 #undef __FUNCT__ 4952 #define __FUNCT__ "SNESGetNPCSide" 4953 /*@ 4954 SNESGetNPCSide - Gets the preconditioning side. 4955 4956 Not Collective 4957 4958 Input Parameter: 4959 . snes - iterative context obtained from SNESCreate() 4960 4961 Output Parameter: 4962 . side - the preconditioning side, where side is one of 4963 .vb 4964 PC_LEFT - left preconditioning (default) 4965 PC_RIGHT - right preconditioning 4966 .ve 4967 4968 Level: intermediate 4969 4970 .keywords: SNES, get, right, left, side, preconditioner, flag 4971 4972 .seealso: SNESSetNPCSide(), KSPGetPCSide() 4973 @*/ 4974 PetscErrorCode SNESGetNPCSide(SNES snes,PCSide *side) 4975 { 4976 PetscFunctionBegin; 4977 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4978 PetscValidPointer(side,2); 4979 *side = snes->pcside; 4980 PetscFunctionReturn(0); 4981 } 4982 4983 #undef __FUNCT__ 4984 #define __FUNCT__ "SNESSetLineSearch" 4985 /*@ 4986 SNESSetLineSearch - Sets the linesearch on the SNES instance. 4987 4988 Collective on SNES 4989 4990 Input Parameters: 4991 + snes - iterative context obtained from SNESCreate() 4992 - linesearch - the linesearch object 4993 4994 Notes: 4995 Use SNESGetLineSearch() to retrieve the preconditioner context (for example, 4996 to configure it using the API). 4997 4998 Level: developer 4999 5000 .keywords: SNES, set, linesearch 5001 .seealso: SNESGetLineSearch() 5002 @*/ 5003 PetscErrorCode SNESSetLineSearch(SNES snes, SNESLineSearch linesearch) 5004 { 5005 PetscErrorCode ierr; 5006 5007 PetscFunctionBegin; 5008 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); 5009 PetscValidHeaderSpecific(linesearch, SNESLINESEARCH_CLASSID, 2); 5010 PetscCheckSameComm(snes, 1, linesearch, 2); 5011 ierr = PetscObjectReference((PetscObject) linesearch);CHKERRQ(ierr); 5012 ierr = SNESLineSearchDestroy(&snes->linesearch);CHKERRQ(ierr); 5013 5014 snes->linesearch = linesearch; 5015 5016 ierr = PetscLogObjectParent((PetscObject)snes, (PetscObject)snes->linesearch);CHKERRQ(ierr); 5017 PetscFunctionReturn(0); 5018 } 5019 5020 #undef __FUNCT__ 5021 #define __FUNCT__ "SNESGetLineSearch" 5022 /*@ 5023 SNESGetLineSearch - Returns a pointer to the line search context set with SNESSetLineSearch() 5024 or creates a default line search instance associated with the SNES and returns it. 5025 5026 Not Collective 5027 5028 Input Parameter: 5029 . snes - iterative context obtained from SNESCreate() 5030 5031 Output Parameter: 5032 . linesearch - linesearch context 5033 5034 Level: beginner 5035 5036 .keywords: SNES, get, linesearch 5037 .seealso: SNESSetLineSearch(), SNESLineSearchCreate() 5038 @*/ 5039 PetscErrorCode SNESGetLineSearch(SNES snes, SNESLineSearch *linesearch) 5040 { 5041 PetscErrorCode ierr; 5042 const char *optionsprefix; 5043 5044 PetscFunctionBegin; 5045 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); 5046 PetscValidPointer(linesearch, 2); 5047 if (!snes->linesearch) { 5048 ierr = SNESGetOptionsPrefix(snes, &optionsprefix);CHKERRQ(ierr); 5049 ierr = SNESLineSearchCreate(PetscObjectComm((PetscObject)snes), &snes->linesearch);CHKERRQ(ierr); 5050 ierr = SNESLineSearchSetSNES(snes->linesearch, snes);CHKERRQ(ierr); 5051 ierr = SNESLineSearchAppendOptionsPrefix(snes->linesearch, optionsprefix);CHKERRQ(ierr); 5052 ierr = PetscObjectIncrementTabLevel((PetscObject) snes->linesearch, (PetscObject) snes, 1);CHKERRQ(ierr); 5053 ierr = PetscLogObjectParent((PetscObject)snes, (PetscObject)snes->linesearch);CHKERRQ(ierr); 5054 } 5055 *linesearch = snes->linesearch; 5056 PetscFunctionReturn(0); 5057 } 5058 5059 #if defined(PETSC_HAVE_MATLAB_ENGINE) 5060 #include <mex.h> 5061 5062 typedef struct {char *funcname; mxArray *ctx;} SNESMatlabContext; 5063 5064 #undef __FUNCT__ 5065 #define __FUNCT__ "SNESComputeFunction_Matlab" 5066 /* 5067 SNESComputeFunction_Matlab - Calls the function that has been set with SNESSetFunctionMatlab(). 5068 5069 Collective on SNES 5070 5071 Input Parameters: 5072 + snes - the SNES context 5073 - x - input vector 5074 5075 Output Parameter: 5076 . y - function vector, as set by SNESSetFunction() 5077 5078 Notes: 5079 SNESComputeFunction() is typically used within nonlinear solvers 5080 implementations, so most users would not generally call this routine 5081 themselves. 5082 5083 Level: developer 5084 5085 .keywords: SNES, nonlinear, compute, function 5086 5087 .seealso: SNESSetFunction(), SNESGetFunction() 5088 */ 5089 PetscErrorCode SNESComputeFunction_Matlab(SNES snes,Vec x,Vec y, void *ctx) 5090 { 5091 PetscErrorCode ierr; 5092 SNESMatlabContext *sctx = (SNESMatlabContext*)ctx; 5093 int nlhs = 1,nrhs = 5; 5094 mxArray *plhs[1],*prhs[5]; 5095 long long int lx = 0,ly = 0,ls = 0; 5096 5097 PetscFunctionBegin; 5098 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 5099 PetscValidHeaderSpecific(x,VEC_CLASSID,2); 5100 PetscValidHeaderSpecific(y,VEC_CLASSID,3); 5101 PetscCheckSameComm(snes,1,x,2); 5102 PetscCheckSameComm(snes,1,y,3); 5103 5104 /* call Matlab function in ctx with arguments x and y */ 5105 5106 ierr = PetscMemcpy(&ls,&snes,sizeof(snes));CHKERRQ(ierr); 5107 ierr = PetscMemcpy(&lx,&x,sizeof(x));CHKERRQ(ierr); 5108 ierr = PetscMemcpy(&ly,&y,sizeof(x));CHKERRQ(ierr); 5109 prhs[0] = mxCreateDoubleScalar((double)ls); 5110 prhs[1] = mxCreateDoubleScalar((double)lx); 5111 prhs[2] = mxCreateDoubleScalar((double)ly); 5112 prhs[3] = mxCreateString(sctx->funcname); 5113 prhs[4] = sctx->ctx; 5114 ierr = mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscSNESComputeFunctionInternal");CHKERRQ(ierr); 5115 ierr = mxGetScalar(plhs[0]);CHKERRQ(ierr); 5116 mxDestroyArray(prhs[0]); 5117 mxDestroyArray(prhs[1]); 5118 mxDestroyArray(prhs[2]); 5119 mxDestroyArray(prhs[3]); 5120 mxDestroyArray(plhs[0]); 5121 PetscFunctionReturn(0); 5122 } 5123 5124 #undef __FUNCT__ 5125 #define __FUNCT__ "SNESSetFunctionMatlab" 5126 /* 5127 SNESSetFunctionMatlab - Sets the function evaluation routine and function 5128 vector for use by the SNES routines in solving systems of nonlinear 5129 equations from MATLAB. Here the function is a string containing the name of a MATLAB function 5130 5131 Logically Collective on SNES 5132 5133 Input Parameters: 5134 + snes - the SNES context 5135 . r - vector to store function value 5136 - f - function evaluation routine 5137 5138 Notes: 5139 The Newton-like methods typically solve linear systems of the form 5140 $ f'(x) x = -f(x), 5141 where f'(x) denotes the Jacobian matrix and f(x) is the function. 5142 5143 Level: beginner 5144 5145 Developer Note: This bleeds the allocated memory SNESMatlabContext *sctx; 5146 5147 .keywords: SNES, nonlinear, set, function 5148 5149 .seealso: SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESSetFunction() 5150 */ 5151 PetscErrorCode SNESSetFunctionMatlab(SNES snes,Vec r,const char *f,mxArray *ctx) 5152 { 5153 PetscErrorCode ierr; 5154 SNESMatlabContext *sctx; 5155 5156 PetscFunctionBegin; 5157 /* currently sctx is memory bleed */ 5158 ierr = PetscNew(&sctx);CHKERRQ(ierr); 5159 ierr = PetscStrallocpy(f,&sctx->funcname);CHKERRQ(ierr); 5160 /* 5161 This should work, but it doesn't 5162 sctx->ctx = ctx; 5163 mexMakeArrayPersistent(sctx->ctx); 5164 */ 5165 sctx->ctx = mxDuplicateArray(ctx); 5166 ierr = SNESSetFunction(snes,r,SNESComputeFunction_Matlab,sctx);CHKERRQ(ierr); 5167 PetscFunctionReturn(0); 5168 } 5169 5170 #undef __FUNCT__ 5171 #define __FUNCT__ "SNESComputeJacobian_Matlab" 5172 /* 5173 SNESComputeJacobian_Matlab - Calls the function that has been set with SNESSetJacobianMatlab(). 5174 5175 Collective on SNES 5176 5177 Input Parameters: 5178 + snes - the SNES context 5179 . x - input vector 5180 . A, B - the matrices 5181 - ctx - user context 5182 5183 Level: developer 5184 5185 .keywords: SNES, nonlinear, compute, function 5186 5187 .seealso: SNESSetFunction(), SNESGetFunction() 5188 @*/ 5189 PetscErrorCode SNESComputeJacobian_Matlab(SNES snes,Vec x,Mat A,Mat B,void *ctx) 5190 { 5191 PetscErrorCode ierr; 5192 SNESMatlabContext *sctx = (SNESMatlabContext*)ctx; 5193 int nlhs = 2,nrhs = 6; 5194 mxArray *plhs[2],*prhs[6]; 5195 long long int lx = 0,lA = 0,ls = 0, lB = 0; 5196 5197 PetscFunctionBegin; 5198 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 5199 PetscValidHeaderSpecific(x,VEC_CLASSID,2); 5200 5201 /* call Matlab function in ctx with arguments x and y */ 5202 5203 ierr = PetscMemcpy(&ls,&snes,sizeof(snes));CHKERRQ(ierr); 5204 ierr = PetscMemcpy(&lx,&x,sizeof(x));CHKERRQ(ierr); 5205 ierr = PetscMemcpy(&lA,A,sizeof(x));CHKERRQ(ierr); 5206 ierr = PetscMemcpy(&lB,B,sizeof(x));CHKERRQ(ierr); 5207 prhs[0] = mxCreateDoubleScalar((double)ls); 5208 prhs[1] = mxCreateDoubleScalar((double)lx); 5209 prhs[2] = mxCreateDoubleScalar((double)lA); 5210 prhs[3] = mxCreateDoubleScalar((double)lB); 5211 prhs[4] = mxCreateString(sctx->funcname); 5212 prhs[5] = sctx->ctx; 5213 ierr = mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscSNESComputeJacobianInternal");CHKERRQ(ierr); 5214 ierr = mxGetScalar(plhs[0]);CHKERRQ(ierr); 5215 mxDestroyArray(prhs[0]); 5216 mxDestroyArray(prhs[1]); 5217 mxDestroyArray(prhs[2]); 5218 mxDestroyArray(prhs[3]); 5219 mxDestroyArray(prhs[4]); 5220 mxDestroyArray(plhs[0]); 5221 mxDestroyArray(plhs[1]); 5222 PetscFunctionReturn(0); 5223 } 5224 5225 #undef __FUNCT__ 5226 #define __FUNCT__ "SNESSetJacobianMatlab" 5227 /* 5228 SNESSetJacobianMatlab - Sets the Jacobian function evaluation routine and two empty Jacobian matrices 5229 vector for use by the SNES routines in solving systems of nonlinear 5230 equations from MATLAB. Here the function is a string containing the name of a MATLAB function 5231 5232 Logically Collective on SNES 5233 5234 Input Parameters: 5235 + snes - the SNES context 5236 . A,B - Jacobian matrices 5237 . J - function evaluation routine 5238 - ctx - user context 5239 5240 Level: developer 5241 5242 Developer Note: This bleeds the allocated memory SNESMatlabContext *sctx; 5243 5244 .keywords: SNES, nonlinear, set, function 5245 5246 .seealso: SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESSetFunction(), J 5247 */ 5248 PetscErrorCode SNESSetJacobianMatlab(SNES snes,Mat A,Mat B,const char *J,mxArray *ctx) 5249 { 5250 PetscErrorCode ierr; 5251 SNESMatlabContext *sctx; 5252 5253 PetscFunctionBegin; 5254 /* currently sctx is memory bleed */ 5255 ierr = PetscNew(&sctx);CHKERRQ(ierr); 5256 ierr = PetscStrallocpy(J,&sctx->funcname);CHKERRQ(ierr); 5257 /* 5258 This should work, but it doesn't 5259 sctx->ctx = ctx; 5260 mexMakeArrayPersistent(sctx->ctx); 5261 */ 5262 sctx->ctx = mxDuplicateArray(ctx); 5263 ierr = SNESSetJacobian(snes,A,B,SNESComputeJacobian_Matlab,sctx);CHKERRQ(ierr); 5264 PetscFunctionReturn(0); 5265 } 5266 5267 #undef __FUNCT__ 5268 #define __FUNCT__ "SNESMonitor_Matlab" 5269 /* 5270 SNESMonitor_Matlab - Calls the function that has been set with SNESMonitorSetMatlab(). 5271 5272 Collective on SNES 5273 5274 .seealso: SNESSetFunction(), SNESGetFunction() 5275 @*/ 5276 PetscErrorCode SNESMonitor_Matlab(SNES snes,PetscInt it, PetscReal fnorm, void *ctx) 5277 { 5278 PetscErrorCode ierr; 5279 SNESMatlabContext *sctx = (SNESMatlabContext*)ctx; 5280 int nlhs = 1,nrhs = 6; 5281 mxArray *plhs[1],*prhs[6]; 5282 long long int lx = 0,ls = 0; 5283 Vec x = snes->vec_sol; 5284 5285 PetscFunctionBegin; 5286 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 5287 5288 ierr = PetscMemcpy(&ls,&snes,sizeof(snes));CHKERRQ(ierr); 5289 ierr = PetscMemcpy(&lx,&x,sizeof(x));CHKERRQ(ierr); 5290 prhs[0] = mxCreateDoubleScalar((double)ls); 5291 prhs[1] = mxCreateDoubleScalar((double)it); 5292 prhs[2] = mxCreateDoubleScalar((double)fnorm); 5293 prhs[3] = mxCreateDoubleScalar((double)lx); 5294 prhs[4] = mxCreateString(sctx->funcname); 5295 prhs[5] = sctx->ctx; 5296 ierr = mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscSNESMonitorInternal");CHKERRQ(ierr); 5297 ierr = mxGetScalar(plhs[0]);CHKERRQ(ierr); 5298 mxDestroyArray(prhs[0]); 5299 mxDestroyArray(prhs[1]); 5300 mxDestroyArray(prhs[2]); 5301 mxDestroyArray(prhs[3]); 5302 mxDestroyArray(prhs[4]); 5303 mxDestroyArray(plhs[0]); 5304 PetscFunctionReturn(0); 5305 } 5306 5307 #undef __FUNCT__ 5308 #define __FUNCT__ "SNESMonitorSetMatlab" 5309 /* 5310 SNESMonitorSetMatlab - Sets the monitor function from MATLAB 5311 5312 Level: developer 5313 5314 Developer Note: This bleeds the allocated memory SNESMatlabContext *sctx; 5315 5316 .keywords: SNES, nonlinear, set, function 5317 5318 .seealso: SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESSetFunction() 5319 */ 5320 PetscErrorCode SNESMonitorSetMatlab(SNES snes,const char *f,mxArray *ctx) 5321 { 5322 PetscErrorCode ierr; 5323 SNESMatlabContext *sctx; 5324 5325 PetscFunctionBegin; 5326 /* currently sctx is memory bleed */ 5327 ierr = PetscNew(&sctx);CHKERRQ(ierr); 5328 ierr = PetscStrallocpy(f,&sctx->funcname);CHKERRQ(ierr); 5329 /* 5330 This should work, but it doesn't 5331 sctx->ctx = ctx; 5332 mexMakeArrayPersistent(sctx->ctx); 5333 */ 5334 sctx->ctx = mxDuplicateArray(ctx); 5335 ierr = SNESMonitorSet(snes,SNESMonitor_Matlab,sctx,NULL);CHKERRQ(ierr); 5336 PetscFunctionReturn(0); 5337 } 5338 5339 #endif 5340