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