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