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