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