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