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