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