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