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