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