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