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