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