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 MatColoring coloring; 2327 ISColoring iscoloring; 2328 MatFDColoring matfdcoloring; 2329 PetscErrorCode (*func)(SNES,Vec,Vec,void*); 2330 void *funcctx; 2331 PetscReal norm1,norm2,normmax; 2332 2333 ierr = MatDuplicate(*B,MAT_DO_NOT_COPY_VALUES,&Bfd);CHKERRQ(ierr); 2334 ierr = MatColoringCreate(Bfd,&coloring);CHKERRQ(ierr); 2335 ierr = MatColoringSetType(coloring,MATCOLORINGSL);CHKERRQ(ierr); 2336 ierr = MatColoringSetFromOptions(coloring);CHKERRQ(ierr); 2337 ierr = MatColoringApply(coloring,&iscoloring);CHKERRQ(ierr); 2338 ierr = MatColoringDestroy(&coloring);CHKERRQ(ierr); 2339 ierr = MatFDColoringCreate(Bfd,iscoloring,&matfdcoloring);CHKERRQ(ierr); 2340 ierr = MatFDColoringSetFromOptions(matfdcoloring);CHKERRQ(ierr); 2341 ierr = MatFDColoringSetUp(Bfd,iscoloring,matfdcoloring);CHKERRQ(ierr); 2342 ierr = ISColoringDestroy(&iscoloring);CHKERRQ(ierr); 2343 2344 /* This method of getting the function is currently unreliable since it doesn't work for DM local functions. */ 2345 ierr = SNESGetFunction(snes,NULL,&func,&funcctx);CHKERRQ(ierr); 2346 ierr = MatFDColoringSetFunction(matfdcoloring,(PetscErrorCode (*)(void))func,funcctx);CHKERRQ(ierr); 2347 ierr = PetscObjectSetOptionsPrefix((PetscObject)matfdcoloring,((PetscObject)snes)->prefix);CHKERRQ(ierr); 2348 ierr = PetscObjectAppendOptionsPrefix((PetscObject)matfdcoloring,"coloring_");CHKERRQ(ierr); 2349 ierr = MatFDColoringSetFromOptions(matfdcoloring);CHKERRQ(ierr); 2350 ierr = MatFDColoringApply(Bfd,matfdcoloring,X,&mstruct,snes);CHKERRQ(ierr); 2351 ierr = MatFDColoringDestroy(&matfdcoloring);CHKERRQ(ierr); 2352 2353 ierr = PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)snes),&vstdout);CHKERRQ(ierr); 2354 if (flag_draw || flag_contour) { 2355 ierr = PetscViewerDrawOpen(PetscObjectComm((PetscObject)snes),0,"Colored Jacobians",PETSC_DECIDE,PETSC_DECIDE,300,300,&vdraw);CHKERRQ(ierr); 2356 if (flag_contour) {ierr = PetscViewerPushFormat(vdraw,PETSC_VIEWER_DRAW_CONTOUR);CHKERRQ(ierr);} 2357 } else vdraw = NULL; 2358 ierr = PetscViewerASCIIPrintf(vstdout,"Explicit preconditioning Jacobian\n");CHKERRQ(ierr); 2359 if (flag_display) {ierr = MatView(*B,vstdout);CHKERRQ(ierr);} 2360 if (vdraw) {ierr = MatView(*B,vdraw);CHKERRQ(ierr);} 2361 ierr = PetscViewerASCIIPrintf(vstdout,"Colored Finite difference Jacobian\n");CHKERRQ(ierr); 2362 if (flag_display) {ierr = MatView(Bfd,vstdout);CHKERRQ(ierr);} 2363 if (vdraw) {ierr = MatView(Bfd,vdraw);CHKERRQ(ierr);} 2364 ierr = MatAYPX(Bfd,-1.0,*B,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 2365 ierr = MatNorm(Bfd,NORM_1,&norm1);CHKERRQ(ierr); 2366 ierr = MatNorm(Bfd,NORM_FROBENIUS,&norm2);CHKERRQ(ierr); 2367 ierr = MatNorm(Bfd,NORM_MAX,&normmax);CHKERRQ(ierr); 2368 ierr = PetscViewerASCIIPrintf(vstdout,"User-provided matrix minus finite difference Jacobian, norm1=%G normFrob=%G normmax=%G\n",norm1,norm2,normmax);CHKERRQ(ierr); 2369 if (flag_display) {ierr = MatView(Bfd,vstdout);CHKERRQ(ierr);} 2370 if (vdraw) { /* Always use contour for the difference */ 2371 ierr = PetscViewerPushFormat(vdraw,PETSC_VIEWER_DRAW_CONTOUR);CHKERRQ(ierr); 2372 ierr = MatView(Bfd,vdraw);CHKERRQ(ierr); 2373 ierr = PetscViewerPopFormat(vdraw);CHKERRQ(ierr); 2374 } 2375 if (flag_contour) {ierr = PetscViewerPopFormat(vdraw);CHKERRQ(ierr);} 2376 2377 if (flag_threshold) { 2378 PetscInt bs,rstart,rend,i; 2379 ierr = MatGetBlockSize(*B,&bs);CHKERRQ(ierr); 2380 ierr = MatGetOwnershipRange(*B,&rstart,&rend);CHKERRQ(ierr); 2381 for (i=rstart; i<rend; i++) { 2382 const PetscScalar *ba,*ca; 2383 const PetscInt *bj,*cj; 2384 PetscInt bn,cn,j,maxentrycol = -1,maxdiffcol = -1,maxrdiffcol = -1; 2385 PetscReal maxentry = 0,maxdiff = 0,maxrdiff = 0; 2386 ierr = MatGetRow(*B,i,&bn,&bj,&ba);CHKERRQ(ierr); 2387 ierr = MatGetRow(Bfd,i,&cn,&cj,&ca);CHKERRQ(ierr); 2388 if (bn != cn) SETERRQ(((PetscObject)*A)->comm,PETSC_ERR_PLIB,"Unexpected different nonzero pattern in -snes_compare_coloring_threshold"); 2389 for (j=0; j<bn; j++) { 2390 PetscReal rdiff = PetscAbsScalar(ca[j]) / (threshold_atol + threshold_rtol*PetscAbsScalar(ba[j])); 2391 if (PetscAbsScalar(ba[j]) > PetscAbs(maxentry)) { 2392 maxentrycol = bj[j]; 2393 maxentry = PetscRealPart(ba[j]); 2394 } 2395 if (PetscAbsScalar(ca[j]) > PetscAbs(maxdiff)) { 2396 maxdiffcol = bj[j]; 2397 maxdiff = PetscRealPart(ca[j]); 2398 } 2399 if (rdiff > maxrdiff) { 2400 maxrdiffcol = bj[j]; 2401 maxrdiff = rdiff; 2402 } 2403 } 2404 if (maxrdiff > 1) { 2405 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); 2406 for (j=0; j<bn; j++) { 2407 PetscReal rdiff; 2408 rdiff = PetscAbsScalar(ca[j]) / (threshold_atol + threshold_rtol*PetscAbsScalar(ba[j])); 2409 if (rdiff > 1) { 2410 ierr = PetscViewerASCIIPrintf(vstdout," (%D,%G:%G)",bj[j],PetscRealPart(ba[j]),PetscRealPart(ca[j]));CHKERRQ(ierr); 2411 } 2412 } 2413 ierr = PetscViewerASCIIPrintf(vstdout,"\n",i,maxentry,maxdiff,maxrdiff);CHKERRQ(ierr); 2414 } 2415 ierr = MatRestoreRow(*B,i,&bn,&bj,&ba);CHKERRQ(ierr); 2416 ierr = MatRestoreRow(Bfd,i,&cn,&cj,&ca);CHKERRQ(ierr); 2417 } 2418 } 2419 ierr = PetscViewerDestroy(&vdraw);CHKERRQ(ierr); 2420 ierr = MatDestroy(&Bfd);CHKERRQ(ierr); 2421 } 2422 } 2423 PetscFunctionReturn(0); 2424 } 2425 2426 /*MC 2427 SNESJacobianFunction - function used to convey the nonlinear Jacobian of the function to be solved by SNES 2428 2429 Synopsis: 2430 #include "petscsnes.h" 2431 $ SNESJacobianFunction(SNES snes,Vec x,Mat *Amat,Mat *Pmat,int *flag,void *ctx); 2432 2433 + x - input vector 2434 . Amat - the matrix that defines the (approximate) Jacobian 2435 . Pmat - the matrix to be used in constructing the preconditioner, usually the same as Amat. 2436 . flag - flag indicating information about the preconditioner matrix 2437 structure (same as flag in KSPSetOperators()), one of SAME_NONZERO_PATTERN,DIFFERENT_NONZERO_PATTERN,SAME_PRECONDITIONER 2438 - ctx - [optional] user-defined Jacobian context 2439 2440 Level: intermediate 2441 2442 .seealso: SNESSetFunction(), SNESGetFunction(), SNESSetJacobian(), SNESGetJacobian() 2443 M*/ 2444 2445 #undef __FUNCT__ 2446 #define __FUNCT__ "SNESSetJacobian" 2447 /*@C 2448 SNESSetJacobian - Sets the function to compute Jacobian as well as the 2449 location to store the matrix. 2450 2451 Logically Collective on SNES and Mat 2452 2453 Input Parameters: 2454 + snes - the SNES context 2455 . Amat - the matrix that defines the (approximate) Jacobian 2456 . Pmat - the matrix to be used in constructing the preconditioner, usually the same as Amat. 2457 . SNESJacobianFunction - Jacobian evaluation routine (if NULL then SNES retains any previously set value) 2458 - ctx - [optional] user-defined context for private data for the 2459 Jacobian evaluation routine (may be NULL) (if NULL then SNES retains any previously set value) 2460 2461 Notes: 2462 See KSPSetOperators() for important information about setting the flag 2463 output parameter in the routine func(). Be sure to read this information! 2464 2465 The routine func() takes Mat * as the matrix arguments rather than Mat. 2466 This allows the Jacobian evaluation routine to replace A and/or B with a 2467 completely new new matrix structure (not just different matrix elements) 2468 when appropriate, for instance, if the nonzero structure is changing 2469 throughout the global iterations. 2470 2471 If the Amat matrix and Pmat matrix are different you must call MatAssemblyBegin/End() on 2472 each matrix. 2473 2474 If using SNESComputeJacobianDefaultColor() to assemble a Jacobian, the ctx argument 2475 must be a MatFDColoring. 2476 2477 Other defect-correction schemes can be used by computing a different matrix in place of the Jacobian. One common 2478 example is to use the "Picard linearization" which only differentiates through the highest order parts of each term. 2479 2480 Level: beginner 2481 2482 .keywords: SNES, nonlinear, set, Jacobian, matrix 2483 2484 .seealso: KSPSetOperators(), SNESSetFunction(), MatMFFDComputeJacobian(), SNESComputeJacobianDefaultColor(), MatStructure, SNESJacobianFunction, SNESSetPicard() 2485 @*/ 2486 PetscErrorCode SNESSetJacobian(SNES snes,Mat Amat,Mat Pmat,PetscErrorCode (*SNESJacobianFunction)(SNES,Vec,Mat*,Mat*,MatStructure*,void*),void *ctx) 2487 { 2488 PetscErrorCode ierr; 2489 DM dm; 2490 2491 PetscFunctionBegin; 2492 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2493 if (Amat) PetscValidHeaderSpecific(Amat,MAT_CLASSID,2); 2494 if (Pmat) PetscValidHeaderSpecific(Pmat,MAT_CLASSID,3); 2495 if (Amat) PetscCheckSameComm(snes,1,Amat,2); 2496 if (Pmat) PetscCheckSameComm(snes,1,Pmat,3); 2497 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 2498 ierr = DMSNESSetJacobian(dm,SNESJacobianFunction,ctx);CHKERRQ(ierr); 2499 if (Amat) { 2500 ierr = PetscObjectReference((PetscObject)Amat);CHKERRQ(ierr); 2501 ierr = MatDestroy(&snes->jacobian);CHKERRQ(ierr); 2502 2503 snes->jacobian = Amat; 2504 } 2505 if (Pmat) { 2506 ierr = PetscObjectReference((PetscObject)Pmat);CHKERRQ(ierr); 2507 ierr = MatDestroy(&snes->jacobian_pre);CHKERRQ(ierr); 2508 2509 snes->jacobian_pre = Pmat; 2510 } 2511 PetscFunctionReturn(0); 2512 } 2513 2514 #undef __FUNCT__ 2515 #define __FUNCT__ "SNESGetJacobian" 2516 /*@C 2517 SNESGetJacobian - Returns the Jacobian matrix and optionally the user 2518 provided context for evaluating the Jacobian. 2519 2520 Not Collective, but Mat object will be parallel if SNES object is 2521 2522 Input Parameter: 2523 . snes - the nonlinear solver context 2524 2525 Output Parameters: 2526 + Amat - location to stash (approximate) Jacobian matrix (or NULL) 2527 . Pmat - location to stash matrix used to compute the preconditioner (or NULL) 2528 . SNESJacobianFunction - location to put Jacobian function (or NULL) 2529 - ctx - location to stash Jacobian ctx (or NULL) 2530 2531 Level: advanced 2532 2533 .seealso: SNESSetJacobian(), SNESComputeJacobian() 2534 @*/ 2535 PetscErrorCode SNESGetJacobian(SNES snes,Mat *Amat,Mat *Pmat,PetscErrorCode (**SNESJacobianFunction)(SNES,Vec,Mat*,Mat*,MatStructure*,void*),void **ctx) 2536 { 2537 PetscErrorCode ierr; 2538 DM dm; 2539 DMSNES sdm; 2540 2541 PetscFunctionBegin; 2542 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2543 if (Amat) *Amat = snes->jacobian; 2544 if (Pmat) *Pmat = snes->jacobian_pre; 2545 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 2546 ierr = DMGetDMSNES(dm,&sdm);CHKERRQ(ierr); 2547 if (SNESJacobianFunction) *SNESJacobianFunction = sdm->ops->computejacobian; 2548 if (ctx) *ctx = sdm->jacobianctx; 2549 PetscFunctionReturn(0); 2550 } 2551 2552 #undef __FUNCT__ 2553 #define __FUNCT__ "SNESSetUp" 2554 /*@ 2555 SNESSetUp - Sets up the internal data structures for the later use 2556 of a nonlinear solver. 2557 2558 Collective on SNES 2559 2560 Input Parameters: 2561 . snes - the SNES context 2562 2563 Notes: 2564 For basic use of the SNES solvers the user need not explicitly call 2565 SNESSetUp(), since these actions will automatically occur during 2566 the call to SNESSolve(). However, if one wishes to control this 2567 phase separately, SNESSetUp() should be called after SNESCreate() 2568 and optional routines of the form SNESSetXXX(), but before SNESSolve(). 2569 2570 Level: advanced 2571 2572 .keywords: SNES, nonlinear, setup 2573 2574 .seealso: SNESCreate(), SNESSolve(), SNESDestroy() 2575 @*/ 2576 PetscErrorCode SNESSetUp(SNES snes) 2577 { 2578 PetscErrorCode ierr; 2579 DM dm; 2580 DMSNES sdm; 2581 SNESLineSearch linesearch, pclinesearch; 2582 void *lsprectx,*lspostctx; 2583 PetscErrorCode (*precheck)(SNESLineSearch,Vec,Vec,PetscBool*,void*); 2584 PetscErrorCode (*postcheck)(SNESLineSearch,Vec,Vec,Vec,PetscBool*,PetscBool*,void*); 2585 PetscErrorCode (*func)(SNES,Vec,Vec,void*); 2586 Vec f,fpc; 2587 void *funcctx; 2588 PetscErrorCode (*jac)(SNES,Vec,Mat*,Mat*,MatStructure*,void*); 2589 void *jacctx,*appctx; 2590 2591 PetscFunctionBegin; 2592 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2593 if (snes->setupcalled) PetscFunctionReturn(0); 2594 2595 if (!((PetscObject)snes)->type_name) { 2596 ierr = SNESSetType(snes,SNESNEWTONLS);CHKERRQ(ierr); 2597 } 2598 2599 ierr = SNESGetFunction(snes,&snes->vec_func,NULL,NULL);CHKERRQ(ierr); 2600 2601 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 2602 ierr = DMGetDMSNES(dm,&sdm);CHKERRQ(ierr); 2603 if (!sdm->ops->computefunction) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"Function never provided to SNES object"); 2604 if (!sdm->ops->computejacobian) { 2605 ierr = DMSNESSetJacobian(dm,SNESComputeJacobianDefaultColor,NULL);CHKERRQ(ierr); 2606 } 2607 if (!snes->vec_func) { 2608 ierr = DMCreateGlobalVector(dm,&snes->vec_func);CHKERRQ(ierr); 2609 } 2610 2611 if (!snes->ksp) { 2612 ierr = SNESGetKSP(snes, &snes->ksp);CHKERRQ(ierr); 2613 } 2614 2615 if (!snes->linesearch) { 2616 ierr = SNESGetLineSearch(snes, &snes->linesearch);CHKERRQ(ierr); 2617 } 2618 ierr = SNESLineSearchSetFunction(snes->linesearch,SNESComputeFunction);CHKERRQ(ierr); 2619 2620 if (snes->pc && (snes->pcside == PC_LEFT)) { 2621 snes->mf = PETSC_TRUE; 2622 snes->mf_operator = PETSC_FALSE; 2623 } 2624 2625 if (snes->mf) { 2626 ierr = SNESSetUpMatrixFree_Private(snes, snes->mf_operator, snes->mf_version);CHKERRQ(ierr); 2627 } 2628 2629 if (snes->ops->usercompute && !snes->user) { 2630 ierr = (*snes->ops->usercompute)(snes,(void**)&snes->user);CHKERRQ(ierr); 2631 } 2632 2633 if (snes->pc) { 2634 /* copy the DM over */ 2635 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 2636 ierr = SNESSetDM(snes->pc,dm);CHKERRQ(ierr); 2637 2638 ierr = SNESGetFunction(snes,&f,&func,&funcctx);CHKERRQ(ierr); 2639 ierr = VecDuplicate(f,&fpc);CHKERRQ(ierr); 2640 ierr = SNESSetFunction(snes->pc,fpc,func,funcctx);CHKERRQ(ierr); 2641 ierr = SNESGetJacobian(snes,NULL,NULL,&jac,&jacctx);CHKERRQ(ierr); 2642 ierr = SNESSetJacobian(snes->pc,NULL,NULL,jac,jacctx);CHKERRQ(ierr); 2643 ierr = SNESGetApplicationContext(snes,&appctx);CHKERRQ(ierr); 2644 ierr = SNESSetApplicationContext(snes->pc,appctx);CHKERRQ(ierr); 2645 ierr = VecDestroy(&fpc);CHKERRQ(ierr); 2646 2647 /* copy the function pointers over */ 2648 ierr = PetscObjectCopyFortranFunctionPointers((PetscObject)snes,(PetscObject)snes->pc);CHKERRQ(ierr); 2649 2650 /* default to 1 iteration */ 2651 ierr = SNESSetTolerances(snes->pc,0.0,0.0,0.0,1,snes->pc->max_funcs);CHKERRQ(ierr); 2652 if (snes->pcside==PC_RIGHT) { 2653 ierr = SNESSetNormSchedule(snes->pc,SNES_NORM_FINAL_ONLY);CHKERRQ(ierr); 2654 } else { 2655 ierr = SNESSetNormSchedule(snes->pc,SNES_NORM_NONE);CHKERRQ(ierr); 2656 } 2657 ierr = SNESSetFromOptions(snes->pc);CHKERRQ(ierr); 2658 2659 /* copy the line search context over */ 2660 ierr = SNESGetLineSearch(snes,&linesearch);CHKERRQ(ierr); 2661 ierr = SNESGetLineSearch(snes->pc,&pclinesearch);CHKERRQ(ierr); 2662 ierr = SNESLineSearchGetPreCheck(linesearch,&precheck,&lsprectx);CHKERRQ(ierr); 2663 ierr = SNESLineSearchGetPostCheck(linesearch,&postcheck,&lspostctx);CHKERRQ(ierr); 2664 ierr = SNESLineSearchSetPreCheck(pclinesearch,precheck,lsprectx);CHKERRQ(ierr); 2665 ierr = SNESLineSearchSetPostCheck(pclinesearch,postcheck,lspostctx);CHKERRQ(ierr); 2666 ierr = PetscObjectCopyFortranFunctionPointers((PetscObject)linesearch, (PetscObject)pclinesearch);CHKERRQ(ierr); 2667 } 2668 2669 snes->jac_iter = 0; 2670 snes->pre_iter = 0; 2671 2672 if (snes->ops->setup) { 2673 ierr = (*snes->ops->setup)(snes);CHKERRQ(ierr); 2674 } 2675 2676 if (snes->pc && (snes->pcside == PC_LEFT)) { 2677 if (snes->functype == SNES_FUNCTION_PRECONDITIONED) { 2678 ierr = SNESGetLineSearch(snes,&linesearch);CHKERRQ(ierr); 2679 ierr = SNESLineSearchSetFunction(linesearch,SNESComputeFunctionDefaultPC);CHKERRQ(ierr); 2680 } 2681 } 2682 2683 snes->setupcalled = PETSC_TRUE; 2684 PetscFunctionReturn(0); 2685 } 2686 2687 #undef __FUNCT__ 2688 #define __FUNCT__ "SNESReset" 2689 /*@ 2690 SNESReset - Resets a SNES context to the snessetupcalled = 0 state and removes any allocated Vecs and Mats 2691 2692 Collective on SNES 2693 2694 Input Parameter: 2695 . snes - iterative context obtained from SNESCreate() 2696 2697 Level: intermediate 2698 2699 Notes: Also calls the application context destroy routine set with SNESSetComputeApplicationContext() 2700 2701 .keywords: SNES, destroy 2702 2703 .seealso: SNESCreate(), SNESSetUp(), SNESSolve() 2704 @*/ 2705 PetscErrorCode SNESReset(SNES snes) 2706 { 2707 PetscErrorCode ierr; 2708 2709 PetscFunctionBegin; 2710 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2711 if (snes->ops->userdestroy && snes->user) { 2712 ierr = (*snes->ops->userdestroy)((void**)&snes->user);CHKERRQ(ierr); 2713 snes->user = NULL; 2714 } 2715 if (snes->pc) { 2716 ierr = SNESReset(snes->pc);CHKERRQ(ierr); 2717 } 2718 2719 if (snes->ops->reset) { 2720 ierr = (*snes->ops->reset)(snes);CHKERRQ(ierr); 2721 } 2722 if (snes->ksp) { 2723 ierr = KSPReset(snes->ksp);CHKERRQ(ierr); 2724 } 2725 2726 if (snes->linesearch) { 2727 ierr = SNESLineSearchReset(snes->linesearch);CHKERRQ(ierr); 2728 } 2729 2730 ierr = VecDestroy(&snes->vec_rhs);CHKERRQ(ierr); 2731 ierr = VecDestroy(&snes->vec_sol);CHKERRQ(ierr); 2732 ierr = VecDestroy(&snes->vec_sol_update);CHKERRQ(ierr); 2733 ierr = VecDestroy(&snes->vec_func);CHKERRQ(ierr); 2734 ierr = MatDestroy(&snes->jacobian);CHKERRQ(ierr); 2735 ierr = MatDestroy(&snes->jacobian_pre);CHKERRQ(ierr); 2736 ierr = VecDestroyVecs(snes->nwork,&snes->work);CHKERRQ(ierr); 2737 ierr = VecDestroyVecs(snes->nvwork,&snes->vwork);CHKERRQ(ierr); 2738 2739 snes->nwork = snes->nvwork = 0; 2740 snes->setupcalled = PETSC_FALSE; 2741 PetscFunctionReturn(0); 2742 } 2743 2744 #undef __FUNCT__ 2745 #define __FUNCT__ "SNESDestroy" 2746 /*@ 2747 SNESDestroy - Destroys the nonlinear solver context that was created 2748 with SNESCreate(). 2749 2750 Collective on SNES 2751 2752 Input Parameter: 2753 . snes - the SNES context 2754 2755 Level: beginner 2756 2757 .keywords: SNES, nonlinear, destroy 2758 2759 .seealso: SNESCreate(), SNESSolve() 2760 @*/ 2761 PetscErrorCode SNESDestroy(SNES *snes) 2762 { 2763 PetscErrorCode ierr; 2764 2765 PetscFunctionBegin; 2766 if (!*snes) PetscFunctionReturn(0); 2767 PetscValidHeaderSpecific((*snes),SNES_CLASSID,1); 2768 if (--((PetscObject)(*snes))->refct > 0) {*snes = 0; PetscFunctionReturn(0);} 2769 2770 ierr = SNESReset((*snes));CHKERRQ(ierr); 2771 ierr = SNESDestroy(&(*snes)->pc);CHKERRQ(ierr); 2772 2773 /* if memory was published with SAWs then destroy it */ 2774 ierr = PetscObjectSAWsViewOff((PetscObject)*snes);CHKERRQ(ierr); 2775 if ((*snes)->ops->destroy) {ierr = (*((*snes))->ops->destroy)((*snes));CHKERRQ(ierr);} 2776 2777 ierr = DMDestroy(&(*snes)->dm);CHKERRQ(ierr); 2778 ierr = KSPDestroy(&(*snes)->ksp);CHKERRQ(ierr); 2779 ierr = SNESLineSearchDestroy(&(*snes)->linesearch);CHKERRQ(ierr); 2780 2781 ierr = PetscFree((*snes)->kspconvctx);CHKERRQ(ierr); 2782 if ((*snes)->ops->convergeddestroy) { 2783 ierr = (*(*snes)->ops->convergeddestroy)((*snes)->cnvP);CHKERRQ(ierr); 2784 } 2785 if ((*snes)->conv_malloc) { 2786 ierr = PetscFree((*snes)->conv_hist);CHKERRQ(ierr); 2787 ierr = PetscFree((*snes)->conv_hist_its);CHKERRQ(ierr); 2788 } 2789 ierr = SNESMonitorCancel((*snes));CHKERRQ(ierr); 2790 ierr = PetscHeaderDestroy(snes);CHKERRQ(ierr); 2791 PetscFunctionReturn(0); 2792 } 2793 2794 /* ----------- Routines to set solver parameters ---------- */ 2795 2796 #undef __FUNCT__ 2797 #define __FUNCT__ "SNESSetLagPreconditioner" 2798 /*@ 2799 SNESSetLagPreconditioner - Determines when the preconditioner is rebuilt in the nonlinear solve. 2800 2801 Logically Collective on SNES 2802 2803 Input Parameters: 2804 + snes - the SNES context 2805 - lag - -1 indicates NEVER rebuild, 1 means rebuild every time the Jacobian is computed within a single nonlinear solve, 2 means every second time 2806 the Jacobian is built etc. -2 indicates rebuild preconditioner at next chance but then never rebuild after that 2807 2808 Options Database Keys: 2809 . -snes_lag_preconditioner <lag> 2810 2811 Notes: 2812 The default is 1 2813 The preconditioner is ALWAYS built in the first iteration of a nonlinear solve unless lag is -1 2814 If -1 is used before the very first nonlinear solve the preconditioner is still built because there is no previous preconditioner to use 2815 2816 Level: intermediate 2817 2818 .keywords: SNES, nonlinear, set, convergence, tolerances 2819 2820 .seealso: SNESSetTrustRegionTolerance(), SNESGetLagPreconditioner(), SNESSetLagJacobian(), SNESGetLagJacobian() 2821 2822 @*/ 2823 PetscErrorCode SNESSetLagPreconditioner(SNES snes,PetscInt lag) 2824 { 2825 PetscFunctionBegin; 2826 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2827 if (lag < -2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Lag must be -2, -1, 1 or greater"); 2828 if (!lag) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Lag cannot be 0"); 2829 PetscValidLogicalCollectiveInt(snes,lag,2); 2830 snes->lagpreconditioner = lag; 2831 PetscFunctionReturn(0); 2832 } 2833 2834 #undef __FUNCT__ 2835 #define __FUNCT__ "SNESSetGridSequence" 2836 /*@ 2837 SNESSetGridSequence - sets the number of steps of grid sequencing that SNES does 2838 2839 Logically Collective on SNES 2840 2841 Input Parameters: 2842 + snes - the SNES context 2843 - steps - the number of refinements to do, defaults to 0 2844 2845 Options Database Keys: 2846 . -snes_grid_sequence <steps> 2847 2848 Level: intermediate 2849 2850 Notes: 2851 Use SNESGetSolution() to extract the fine grid solution after grid sequencing. 2852 2853 .keywords: SNES, nonlinear, set, convergence, tolerances 2854 2855 .seealso: SNESSetTrustRegionTolerance(), SNESGetLagPreconditioner(), SNESSetLagJacobian(), SNESGetLagJacobian() 2856 2857 @*/ 2858 PetscErrorCode SNESSetGridSequence(SNES snes,PetscInt steps) 2859 { 2860 PetscFunctionBegin; 2861 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2862 PetscValidLogicalCollectiveInt(snes,steps,2); 2863 snes->gridsequence = steps; 2864 PetscFunctionReturn(0); 2865 } 2866 2867 #undef __FUNCT__ 2868 #define __FUNCT__ "SNESGetLagPreconditioner" 2869 /*@ 2870 SNESGetLagPreconditioner - Indicates how often the preconditioner is rebuilt 2871 2872 Not Collective 2873 2874 Input Parameter: 2875 . snes - the SNES context 2876 2877 Output Parameter: 2878 . lag - -1 indicates NEVER rebuild, 1 means rebuild every time the Jacobian is computed within a single nonlinear solve, 2 means every second time 2879 the Jacobian is built etc. -2 indicates rebuild preconditioner at next chance but then never rebuild after that 2880 2881 Options Database Keys: 2882 . -snes_lag_preconditioner <lag> 2883 2884 Notes: 2885 The default is 1 2886 The preconditioner is ALWAYS built in the first iteration of a nonlinear solve unless lag is -1 2887 2888 Level: intermediate 2889 2890 .keywords: SNES, nonlinear, set, convergence, tolerances 2891 2892 .seealso: SNESSetTrustRegionTolerance(), SNESSetLagPreconditioner() 2893 2894 @*/ 2895 PetscErrorCode SNESGetLagPreconditioner(SNES snes,PetscInt *lag) 2896 { 2897 PetscFunctionBegin; 2898 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2899 *lag = snes->lagpreconditioner; 2900 PetscFunctionReturn(0); 2901 } 2902 2903 #undef __FUNCT__ 2904 #define __FUNCT__ "SNESSetLagJacobian" 2905 /*@ 2906 SNESSetLagJacobian - Determines when the Jacobian is rebuilt in the nonlinear solve. See SNESSetLagPreconditioner() for determining how 2907 often the preconditioner is rebuilt. 2908 2909 Logically Collective on SNES 2910 2911 Input Parameters: 2912 + snes - the SNES context 2913 - lag - -1 indicates NEVER rebuild, 1 means rebuild every time the Jacobian is computed within a single nonlinear solve, 2 means every second time 2914 the Jacobian is built etc. -2 means rebuild at next chance but then never again 2915 2916 Options Database Keys: 2917 . -snes_lag_jacobian <lag> 2918 2919 Notes: 2920 The default is 1 2921 The Jacobian is ALWAYS built in the first iteration of a nonlinear solve unless lag is -1 2922 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 2923 at the next Newton step but never again (unless it is reset to another value) 2924 2925 Level: intermediate 2926 2927 .keywords: SNES, nonlinear, set, convergence, tolerances 2928 2929 .seealso: SNESSetTrustRegionTolerance(), SNESGetLagPreconditioner(), SNESSetLagPreconditioner(), SNESGetLagJacobian() 2930 2931 @*/ 2932 PetscErrorCode SNESSetLagJacobian(SNES snes,PetscInt lag) 2933 { 2934 PetscFunctionBegin; 2935 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2936 if (lag < -2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Lag must be -2, -1, 1 or greater"); 2937 if (!lag) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Lag cannot be 0"); 2938 PetscValidLogicalCollectiveInt(snes,lag,2); 2939 snes->lagjacobian = lag; 2940 PetscFunctionReturn(0); 2941 } 2942 2943 #undef __FUNCT__ 2944 #define __FUNCT__ "SNESGetLagJacobian" 2945 /*@ 2946 SNESGetLagJacobian - Indicates how often the Jacobian is rebuilt. See SNESGetLagPreconditioner() to determine when the preconditioner is rebuilt 2947 2948 Not Collective 2949 2950 Input Parameter: 2951 . snes - the SNES context 2952 2953 Output Parameter: 2954 . lag - -1 indicates NEVER rebuild, 1 means rebuild every time the Jacobian is computed within a single nonlinear solve, 2 means every second time 2955 the Jacobian is built etc. 2956 2957 Options Database Keys: 2958 . -snes_lag_jacobian <lag> 2959 2960 Notes: 2961 The default is 1 2962 The jacobian is ALWAYS built in the first iteration of a nonlinear solve unless lag is -1 2963 2964 Level: intermediate 2965 2966 .keywords: SNES, nonlinear, set, convergence, tolerances 2967 2968 .seealso: SNESSetTrustRegionTolerance(), SNESSetLagJacobian(), SNESSetLagPreconditioner(), SNESGetLagPreconditioner() 2969 2970 @*/ 2971 PetscErrorCode SNESGetLagJacobian(SNES snes,PetscInt *lag) 2972 { 2973 PetscFunctionBegin; 2974 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2975 *lag = snes->lagjacobian; 2976 PetscFunctionReturn(0); 2977 } 2978 2979 #undef __FUNCT__ 2980 #define __FUNCT__ "SNESSetLagJacobianPersists" 2981 /*@ 2982 SNESSetLagJacobianPersists - Set whether or not the Jacobian lagging persists through multiple solves 2983 2984 Logically collective on SNES 2985 2986 Input Parameter: 2987 + snes - the SNES context 2988 - flg - jacobian lagging persists if true 2989 2990 Options Database Keys: 2991 . -snes_lag_jacobian_persists <flg> 2992 2993 Notes: This is useful both for nonlinear preconditioning, where it's appropriate to have the Jacobian be stale by 2994 several solves, and for implicit time-stepping, where Jacobian lagging in the inner nonlinear solve over several 2995 timesteps may present huge efficiency gains. 2996 2997 Level: developer 2998 2999 .keywords: SNES, nonlinear, lag, NPC 3000 3001 .seealso: SNESSetLagPreconditionerPersists(), SNESSetLagJacobian(), SNESGetLagJacobian(), SNESGetPC() 3002 3003 @*/ 3004 PetscErrorCode SNESSetLagJacobianPersists(SNES snes,PetscBool flg) 3005 { 3006 PetscFunctionBegin; 3007 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3008 PetscValidLogicalCollectiveBool(snes,flg,2); 3009 snes->lagjac_persist = flg; 3010 PetscFunctionReturn(0); 3011 } 3012 3013 #undef __FUNCT__ 3014 #define __FUNCT__ "SNESSetLagPreconditionerPersists" 3015 /*@ 3016 SNESSetLagPreconditionerPersists - Set whether or not the preconditioner lagging persists through multiple solves 3017 3018 Logically Collective on SNES 3019 3020 Input Parameter: 3021 + snes - the SNES context 3022 - flg - preconditioner lagging persists if true 3023 3024 Options Database Keys: 3025 . -snes_lag_jacobian_persists <flg> 3026 3027 Notes: This is useful both for nonlinear preconditioning, where it's appropriate to have the preconditioner be stale 3028 by several solves, and for implicit time-stepping, where preconditioner lagging in the inner nonlinear solve over 3029 several timesteps may present huge efficiency gains. 3030 3031 Level: developer 3032 3033 .keywords: SNES, nonlinear, lag, NPC 3034 3035 .seealso: SNESSetLagJacobianPersists(), SNESSetLagJacobian(), SNESGetLagJacobian(), SNESGetPC() 3036 3037 @*/ 3038 PetscErrorCode SNESSetLagPreconditionerPersists(SNES snes,PetscBool flg) 3039 { 3040 PetscFunctionBegin; 3041 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3042 PetscValidLogicalCollectiveBool(snes,flg,2); 3043 snes->lagpre_persist = flg; 3044 PetscFunctionReturn(0); 3045 } 3046 3047 #undef __FUNCT__ 3048 #define __FUNCT__ "SNESSetTolerances" 3049 /*@ 3050 SNESSetTolerances - Sets various parameters used in convergence tests. 3051 3052 Logically Collective on SNES 3053 3054 Input Parameters: 3055 + snes - the SNES context 3056 . abstol - absolute convergence tolerance 3057 . rtol - relative convergence tolerance 3058 . stol - convergence tolerance in terms of the norm of the change in the solution between steps, || delta x || < stol*|| x || 3059 . maxit - maximum number of iterations 3060 - maxf - maximum number of function evaluations 3061 3062 Options Database Keys: 3063 + -snes_atol <abstol> - Sets abstol 3064 . -snes_rtol <rtol> - Sets rtol 3065 . -snes_stol <stol> - Sets stol 3066 . -snes_max_it <maxit> - Sets maxit 3067 - -snes_max_funcs <maxf> - Sets maxf 3068 3069 Notes: 3070 The default maximum number of iterations is 50. 3071 The default maximum number of function evaluations is 1000. 3072 3073 Level: intermediate 3074 3075 .keywords: SNES, nonlinear, set, convergence, tolerances 3076 3077 .seealso: SNESSetTrustRegionTolerance() 3078 @*/ 3079 PetscErrorCode SNESSetTolerances(SNES snes,PetscReal abstol,PetscReal rtol,PetscReal stol,PetscInt maxit,PetscInt maxf) 3080 { 3081 PetscFunctionBegin; 3082 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3083 PetscValidLogicalCollectiveReal(snes,abstol,2); 3084 PetscValidLogicalCollectiveReal(snes,rtol,3); 3085 PetscValidLogicalCollectiveReal(snes,stol,4); 3086 PetscValidLogicalCollectiveInt(snes,maxit,5); 3087 PetscValidLogicalCollectiveInt(snes,maxf,6); 3088 3089 if (abstol != PETSC_DEFAULT) { 3090 if (abstol < 0.0) SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_OUTOFRANGE,"Absolute tolerance %G must be non-negative",abstol); 3091 snes->abstol = abstol; 3092 } 3093 if (rtol != PETSC_DEFAULT) { 3094 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); 3095 snes->rtol = rtol; 3096 } 3097 if (stol != PETSC_DEFAULT) { 3098 if (stol < 0.0) SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_OUTOFRANGE,"Step tolerance %G must be non-negative",stol); 3099 snes->stol = stol; 3100 } 3101 if (maxit != PETSC_DEFAULT) { 3102 if (maxit < 0) SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_OUTOFRANGE,"Maximum number of iterations %D must be non-negative",maxit); 3103 snes->max_its = maxit; 3104 } 3105 if (maxf != PETSC_DEFAULT) { 3106 if (maxf < 0) SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_OUTOFRANGE,"Maximum number of function evaluations %D must be non-negative",maxf); 3107 snes->max_funcs = maxf; 3108 } 3109 snes->tolerancesset = PETSC_TRUE; 3110 PetscFunctionReturn(0); 3111 } 3112 3113 #undef __FUNCT__ 3114 #define __FUNCT__ "SNESGetTolerances" 3115 /*@ 3116 SNESGetTolerances - Gets various parameters used in convergence tests. 3117 3118 Not Collective 3119 3120 Input Parameters: 3121 + snes - the SNES context 3122 . atol - absolute convergence tolerance 3123 . rtol - relative convergence tolerance 3124 . stol - convergence tolerance in terms of the norm 3125 of the change in the solution between steps 3126 . maxit - maximum number of iterations 3127 - maxf - maximum number of function evaluations 3128 3129 Notes: 3130 The user can specify NULL for any parameter that is not needed. 3131 3132 Level: intermediate 3133 3134 .keywords: SNES, nonlinear, get, convergence, tolerances 3135 3136 .seealso: SNESSetTolerances() 3137 @*/ 3138 PetscErrorCode SNESGetTolerances(SNES snes,PetscReal *atol,PetscReal *rtol,PetscReal *stol,PetscInt *maxit,PetscInt *maxf) 3139 { 3140 PetscFunctionBegin; 3141 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3142 if (atol) *atol = snes->abstol; 3143 if (rtol) *rtol = snes->rtol; 3144 if (stol) *stol = snes->stol; 3145 if (maxit) *maxit = snes->max_its; 3146 if (maxf) *maxf = snes->max_funcs; 3147 PetscFunctionReturn(0); 3148 } 3149 3150 #undef __FUNCT__ 3151 #define __FUNCT__ "SNESSetTrustRegionTolerance" 3152 /*@ 3153 SNESSetTrustRegionTolerance - Sets the trust region parameter tolerance. 3154 3155 Logically Collective on SNES 3156 3157 Input Parameters: 3158 + snes - the SNES context 3159 - tol - tolerance 3160 3161 Options Database Key: 3162 . -snes_trtol <tol> - Sets tol 3163 3164 Level: intermediate 3165 3166 .keywords: SNES, nonlinear, set, trust region, tolerance 3167 3168 .seealso: SNESSetTolerances() 3169 @*/ 3170 PetscErrorCode SNESSetTrustRegionTolerance(SNES snes,PetscReal tol) 3171 { 3172 PetscFunctionBegin; 3173 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3174 PetscValidLogicalCollectiveReal(snes,tol,2); 3175 snes->deltatol = tol; 3176 PetscFunctionReturn(0); 3177 } 3178 3179 /* 3180 Duplicate the lg monitors for SNES from KSP; for some reason with 3181 dynamic libraries things don't work under Sun4 if we just use 3182 macros instead of functions 3183 */ 3184 #undef __FUNCT__ 3185 #define __FUNCT__ "SNESMonitorLGResidualNorm" 3186 PetscErrorCode SNESMonitorLGResidualNorm(SNES snes,PetscInt it,PetscReal norm,void *ctx) 3187 { 3188 PetscErrorCode ierr; 3189 3190 PetscFunctionBegin; 3191 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3192 ierr = KSPMonitorLGResidualNorm((KSP)snes,it,norm,ctx);CHKERRQ(ierr); 3193 PetscFunctionReturn(0); 3194 } 3195 3196 #undef __FUNCT__ 3197 #define __FUNCT__ "SNESMonitorLGCreate" 3198 PetscErrorCode SNESMonitorLGCreate(const char host[],const char label[],int x,int y,int m,int n,PetscDrawLG *draw) 3199 { 3200 PetscErrorCode ierr; 3201 3202 PetscFunctionBegin; 3203 ierr = KSPMonitorLGResidualNormCreate(host,label,x,y,m,n,draw);CHKERRQ(ierr); 3204 PetscFunctionReturn(0); 3205 } 3206 3207 #undef __FUNCT__ 3208 #define __FUNCT__ "SNESMonitorLGDestroy" 3209 PetscErrorCode SNESMonitorLGDestroy(PetscDrawLG *draw) 3210 { 3211 PetscErrorCode ierr; 3212 3213 PetscFunctionBegin; 3214 ierr = KSPMonitorLGResidualNormDestroy(draw);CHKERRQ(ierr); 3215 PetscFunctionReturn(0); 3216 } 3217 3218 extern PetscErrorCode SNESMonitorRange_Private(SNES,PetscInt,PetscReal*); 3219 #undef __FUNCT__ 3220 #define __FUNCT__ "SNESMonitorLGRange" 3221 PetscErrorCode SNESMonitorLGRange(SNES snes,PetscInt n,PetscReal rnorm,void *monctx) 3222 { 3223 PetscDrawLG lg; 3224 PetscErrorCode ierr; 3225 PetscReal x,y,per; 3226 PetscViewer v = (PetscViewer)monctx; 3227 static PetscReal prev; /* should be in the context */ 3228 PetscDraw draw; 3229 3230 PetscFunctionBegin; 3231 ierr = PetscViewerDrawGetDrawLG(v,0,&lg);CHKERRQ(ierr); 3232 if (!n) {ierr = PetscDrawLGReset(lg);CHKERRQ(ierr);} 3233 ierr = PetscDrawLGGetDraw(lg,&draw);CHKERRQ(ierr); 3234 ierr = PetscDrawSetTitle(draw,"Residual norm");CHKERRQ(ierr); 3235 x = (PetscReal)n; 3236 if (rnorm > 0.0) y = log10(rnorm); 3237 else y = -15.0; 3238 ierr = PetscDrawLGAddPoint(lg,&x,&y);CHKERRQ(ierr); 3239 if (n < 20 || !(n % 5)) { 3240 ierr = PetscDrawLGDraw(lg);CHKERRQ(ierr); 3241 } 3242 3243 ierr = PetscViewerDrawGetDrawLG(v,1,&lg);CHKERRQ(ierr); 3244 if (!n) {ierr = PetscDrawLGReset(lg);CHKERRQ(ierr);} 3245 ierr = PetscDrawLGGetDraw(lg,&draw);CHKERRQ(ierr); 3246 ierr = PetscDrawSetTitle(draw,"% elemts > .2*max elemt");CHKERRQ(ierr); 3247 ierr = SNESMonitorRange_Private(snes,n,&per);CHKERRQ(ierr); 3248 x = (PetscReal)n; 3249 y = 100.0*per; 3250 ierr = PetscDrawLGAddPoint(lg,&x,&y);CHKERRQ(ierr); 3251 if (n < 20 || !(n % 5)) { 3252 ierr = PetscDrawLGDraw(lg);CHKERRQ(ierr); 3253 } 3254 3255 ierr = PetscViewerDrawGetDrawLG(v,2,&lg);CHKERRQ(ierr); 3256 if (!n) {prev = rnorm;ierr = PetscDrawLGReset(lg);CHKERRQ(ierr);} 3257 ierr = PetscDrawLGGetDraw(lg,&draw);CHKERRQ(ierr); 3258 ierr = PetscDrawSetTitle(draw,"(norm -oldnorm)/oldnorm");CHKERRQ(ierr); 3259 x = (PetscReal)n; 3260 y = (prev - rnorm)/prev; 3261 ierr = PetscDrawLGAddPoint(lg,&x,&y);CHKERRQ(ierr); 3262 if (n < 20 || !(n % 5)) { 3263 ierr = PetscDrawLGDraw(lg);CHKERRQ(ierr); 3264 } 3265 3266 ierr = PetscViewerDrawGetDrawLG(v,3,&lg);CHKERRQ(ierr); 3267 if (!n) {ierr = PetscDrawLGReset(lg);CHKERRQ(ierr);} 3268 ierr = PetscDrawLGGetDraw(lg,&draw);CHKERRQ(ierr); 3269 ierr = PetscDrawSetTitle(draw,"(norm -oldnorm)/oldnorm*(% > .2 max)");CHKERRQ(ierr); 3270 x = (PetscReal)n; 3271 y = (prev - rnorm)/(prev*per); 3272 if (n > 2) { /*skip initial crazy value */ 3273 ierr = PetscDrawLGAddPoint(lg,&x,&y);CHKERRQ(ierr); 3274 } 3275 if (n < 20 || !(n % 5)) { 3276 ierr = PetscDrawLGDraw(lg);CHKERRQ(ierr); 3277 } 3278 prev = rnorm; 3279 PetscFunctionReturn(0); 3280 } 3281 3282 #undef __FUNCT__ 3283 #define __FUNCT__ "SNESMonitor" 3284 /*@ 3285 SNESMonitor - runs the user provided monitor routines, if they exist 3286 3287 Collective on SNES 3288 3289 Input Parameters: 3290 + snes - nonlinear solver context obtained from SNESCreate() 3291 . iter - iteration number 3292 - rnorm - relative norm of the residual 3293 3294 Notes: 3295 This routine is called by the SNES implementations. 3296 It does not typically need to be called by the user. 3297 3298 Level: developer 3299 3300 .seealso: SNESMonitorSet() 3301 @*/ 3302 PetscErrorCode SNESMonitor(SNES snes,PetscInt iter,PetscReal rnorm) 3303 { 3304 PetscErrorCode ierr; 3305 PetscInt i,n = snes->numbermonitors; 3306 3307 PetscFunctionBegin; 3308 for (i=0; i<n; i++) { 3309 ierr = (*snes->monitor[i])(snes,iter,rnorm,snes->monitorcontext[i]);CHKERRQ(ierr); 3310 } 3311 PetscFunctionReturn(0); 3312 } 3313 3314 /* ------------ Routines to set performance monitoring options ----------- */ 3315 3316 /*MC 3317 SNESMonitorFunction - functional form passed to SNESMonitorSet() to monitor convergence of nonlinear solver 3318 3319 Synopsis: 3320 #include "petscsnes.h" 3321 $ PetscErrorCode SNESMonitorFunction(SNES snes,PetscInt its, PetscReal norm,void *mctx) 3322 3323 + snes - the SNES context 3324 . its - iteration number 3325 . norm - 2-norm function value (may be estimated) 3326 - mctx - [optional] monitoring context 3327 3328 Level: advanced 3329 3330 .seealso: SNESMonitorSet(), SNESMonitorGet() 3331 M*/ 3332 3333 #undef __FUNCT__ 3334 #define __FUNCT__ "SNESMonitorSet" 3335 /*@C 3336 SNESMonitorSet - Sets an ADDITIONAL function that is to be used at every 3337 iteration of the nonlinear solver to display the iteration's 3338 progress. 3339 3340 Logically Collective on SNES 3341 3342 Input Parameters: 3343 + snes - the SNES context 3344 . SNESMonitorFunction - monitoring routine 3345 . mctx - [optional] user-defined context for private data for the 3346 monitor routine (use NULL if no context is desired) 3347 - monitordestroy - [optional] routine that frees monitor context 3348 (may be NULL) 3349 3350 Options Database Keys: 3351 + -snes_monitor - sets SNESMonitorDefault() 3352 . -snes_monitor_lg_residualnorm - sets line graph monitor, 3353 uses SNESMonitorLGCreate() 3354 - -snes_monitor_cancel - cancels all monitors that have 3355 been hardwired into a code by 3356 calls to SNESMonitorSet(), but 3357 does not cancel those set via 3358 the options database. 3359 3360 Notes: 3361 Several different monitoring routines may be set by calling 3362 SNESMonitorSet() multiple times; all will be called in the 3363 order in which they were set. 3364 3365 Fortran notes: Only a single monitor function can be set for each SNES object 3366 3367 Level: intermediate 3368 3369 .keywords: SNES, nonlinear, set, monitor 3370 3371 .seealso: SNESMonitorDefault(), SNESMonitorCancel(), SNESMonitorFunction 3372 @*/ 3373 PetscErrorCode SNESMonitorSet(SNES snes,PetscErrorCode (*SNESMonitorFunction)(SNES,PetscInt,PetscReal,void*),void *mctx,PetscErrorCode (*monitordestroy)(void**)) 3374 { 3375 PetscInt i; 3376 PetscErrorCode ierr; 3377 3378 PetscFunctionBegin; 3379 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3380 if (snes->numbermonitors >= MAXSNESMONITORS) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Too many monitors set"); 3381 for (i=0; i<snes->numbermonitors;i++) { 3382 if (SNESMonitorFunction == snes->monitor[i] && monitordestroy == snes->monitordestroy[i] && mctx == snes->monitorcontext[i]) { 3383 if (monitordestroy) { 3384 ierr = (*monitordestroy)(&mctx);CHKERRQ(ierr); 3385 } 3386 PetscFunctionReturn(0); 3387 } 3388 } 3389 snes->monitor[snes->numbermonitors] = SNESMonitorFunction; 3390 snes->monitordestroy[snes->numbermonitors] = monitordestroy; 3391 snes->monitorcontext[snes->numbermonitors++] = (void*)mctx; 3392 PetscFunctionReturn(0); 3393 } 3394 3395 #undef __FUNCT__ 3396 #define __FUNCT__ "SNESMonitorCancel" 3397 /*@ 3398 SNESMonitorCancel - Clears all the monitor functions for a SNES object. 3399 3400 Logically Collective on SNES 3401 3402 Input Parameters: 3403 . snes - the SNES context 3404 3405 Options Database Key: 3406 . -snes_monitor_cancel - cancels all monitors that have been hardwired 3407 into a code by calls to SNESMonitorSet(), but does not cancel those 3408 set via the options database 3409 3410 Notes: 3411 There is no way to clear one specific monitor from a SNES object. 3412 3413 Level: intermediate 3414 3415 .keywords: SNES, nonlinear, set, monitor 3416 3417 .seealso: SNESMonitorDefault(), SNESMonitorSet() 3418 @*/ 3419 PetscErrorCode SNESMonitorCancel(SNES snes) 3420 { 3421 PetscErrorCode ierr; 3422 PetscInt i; 3423 3424 PetscFunctionBegin; 3425 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3426 for (i=0; i<snes->numbermonitors; i++) { 3427 if (snes->monitordestroy[i]) { 3428 ierr = (*snes->monitordestroy[i])(&snes->monitorcontext[i]);CHKERRQ(ierr); 3429 } 3430 } 3431 snes->numbermonitors = 0; 3432 PetscFunctionReturn(0); 3433 } 3434 3435 /*MC 3436 SNESConvergenceTestFunction - functional form used for testing of convergence of nonlinear solver 3437 3438 Synopsis: 3439 #include "petscsnes.h" 3440 $ PetscErrorCode SNESConvergenceTest(SNES snes,PetscInt it,PetscReal xnorm,PetscReal gnorm,PetscReal f,SNESConvergedReason *reason,void *cctx) 3441 3442 + snes - the SNES context 3443 . it - current iteration (0 is the first and is before any Newton step) 3444 . cctx - [optional] convergence context 3445 . reason - reason for convergence/divergence 3446 . xnorm - 2-norm of current iterate 3447 . gnorm - 2-norm of current step 3448 - f - 2-norm of function 3449 3450 Level: intermediate 3451 3452 .seealso: SNESSetConvergenceTest(), SNESGetConvergenceTest() 3453 M*/ 3454 3455 #undef __FUNCT__ 3456 #define __FUNCT__ "SNESSetConvergenceTest" 3457 /*@C 3458 SNESSetConvergenceTest - Sets the function that is to be used 3459 to test for convergence of the nonlinear iterative solution. 3460 3461 Logically Collective on SNES 3462 3463 Input Parameters: 3464 + snes - the SNES context 3465 . SNESConvergenceTestFunction - routine to test for convergence 3466 . cctx - [optional] context for private data for the convergence routine (may be NULL) 3467 - destroy - [optional] destructor for the context (may be NULL; NULL_FUNCTION in Fortran) 3468 3469 Level: advanced 3470 3471 .keywords: SNES, nonlinear, set, convergence, test 3472 3473 .seealso: SNESConvergedDefault(), SNESConvergedSkip(), SNESConvergenceTestFunction 3474 @*/ 3475 PetscErrorCode SNESSetConvergenceTest(SNES snes,PetscErrorCode (*SNESConvergenceTestFunction)(SNES,PetscInt,PetscReal,PetscReal,PetscReal,SNESConvergedReason*,void*),void *cctx,PetscErrorCode (*destroy)(void*)) 3476 { 3477 PetscErrorCode ierr; 3478 3479 PetscFunctionBegin; 3480 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3481 if (!SNESConvergenceTestFunction) SNESConvergenceTestFunction = SNESConvergedSkip; 3482 if (snes->ops->convergeddestroy) { 3483 ierr = (*snes->ops->convergeddestroy)(snes->cnvP);CHKERRQ(ierr); 3484 } 3485 snes->ops->converged = SNESConvergenceTestFunction; 3486 snes->ops->convergeddestroy = destroy; 3487 snes->cnvP = cctx; 3488 PetscFunctionReturn(0); 3489 } 3490 3491 #undef __FUNCT__ 3492 #define __FUNCT__ "SNESGetConvergedReason" 3493 /*@ 3494 SNESGetConvergedReason - Gets the reason the SNES iteration was stopped. 3495 3496 Not Collective 3497 3498 Input Parameter: 3499 . snes - the SNES context 3500 3501 Output Parameter: 3502 . reason - negative value indicates diverged, positive value converged, see SNESConvergedReason or the 3503 manual pages for the individual convergence tests for complete lists 3504 3505 Level: intermediate 3506 3507 Notes: Can only be called after the call the SNESSolve() is complete. 3508 3509 .keywords: SNES, nonlinear, set, convergence, test 3510 3511 .seealso: SNESSetConvergenceTest(), SNESConvergedReason 3512 @*/ 3513 PetscErrorCode SNESGetConvergedReason(SNES snes,SNESConvergedReason *reason) 3514 { 3515 PetscFunctionBegin; 3516 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3517 PetscValidPointer(reason,2); 3518 *reason = snes->reason; 3519 PetscFunctionReturn(0); 3520 } 3521 3522 #undef __FUNCT__ 3523 #define __FUNCT__ "SNESSetConvergenceHistory" 3524 /*@ 3525 SNESSetConvergenceHistory - Sets the array used to hold the convergence history. 3526 3527 Logically Collective on SNES 3528 3529 Input Parameters: 3530 + snes - iterative context obtained from SNESCreate() 3531 . a - array to hold history, this array will contain the function norms computed at each step 3532 . its - integer array holds the number of linear iterations for each solve. 3533 . na - size of a and its 3534 - reset - PETSC_TRUE indicates each new nonlinear solve resets the history counter to zero, 3535 else it continues storing new values for new nonlinear solves after the old ones 3536 3537 Notes: 3538 If 'a' and 'its' are NULL then space is allocated for the history. If 'na' PETSC_DECIDE or PETSC_DEFAULT then a 3539 default array of length 10000 is allocated. 3540 3541 This routine is useful, e.g., when running a code for purposes 3542 of accurate performance monitoring, when no I/O should be done 3543 during the section of code that is being timed. 3544 3545 Level: intermediate 3546 3547 .keywords: SNES, set, convergence, history 3548 3549 .seealso: SNESGetConvergenceHistory() 3550 3551 @*/ 3552 PetscErrorCode SNESSetConvergenceHistory(SNES snes,PetscReal a[],PetscInt its[],PetscInt na,PetscBool reset) 3553 { 3554 PetscErrorCode ierr; 3555 3556 PetscFunctionBegin; 3557 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3558 if (a) PetscValidScalarPointer(a,2); 3559 if (its) PetscValidIntPointer(its,3); 3560 if (!a) { 3561 if (na == PETSC_DECIDE || na == PETSC_DEFAULT) na = 1000; 3562 ierr = PetscMalloc(na*sizeof(PetscReal),&a);CHKERRQ(ierr); 3563 ierr = PetscMalloc(na*sizeof(PetscInt),&its);CHKERRQ(ierr); 3564 3565 snes->conv_malloc = PETSC_TRUE; 3566 } 3567 snes->conv_hist = a; 3568 snes->conv_hist_its = its; 3569 snes->conv_hist_max = na; 3570 snes->conv_hist_len = 0; 3571 snes->conv_hist_reset = reset; 3572 PetscFunctionReturn(0); 3573 } 3574 3575 #if defined(PETSC_HAVE_MATLAB_ENGINE) 3576 #include <engine.h> /* MATLAB include file */ 3577 #include <mex.h> /* MATLAB include file */ 3578 3579 #undef __FUNCT__ 3580 #define __FUNCT__ "SNESGetConvergenceHistoryMatlab" 3581 PETSC_EXTERN mxArray *SNESGetConvergenceHistoryMatlab(SNES snes) 3582 { 3583 mxArray *mat; 3584 PetscInt i; 3585 PetscReal *ar; 3586 3587 PetscFunctionBegin; 3588 mat = mxCreateDoubleMatrix(snes->conv_hist_len,1,mxREAL); 3589 ar = (PetscReal*) mxGetData(mat); 3590 for (i=0; i<snes->conv_hist_len; i++) ar[i] = snes->conv_hist[i]; 3591 PetscFunctionReturn(mat); 3592 } 3593 #endif 3594 3595 #undef __FUNCT__ 3596 #define __FUNCT__ "SNESGetConvergenceHistory" 3597 /*@C 3598 SNESGetConvergenceHistory - Gets the array used to hold the convergence history. 3599 3600 Not Collective 3601 3602 Input Parameter: 3603 . snes - iterative context obtained from SNESCreate() 3604 3605 Output Parameters: 3606 . a - array to hold history 3607 . its - integer array holds the number of linear iterations (or 3608 negative if not converged) for each solve. 3609 - na - size of a and its 3610 3611 Notes: 3612 The calling sequence for this routine in Fortran is 3613 $ call SNESGetConvergenceHistory(SNES snes, integer na, integer ierr) 3614 3615 This routine is useful, e.g., when running a code for purposes 3616 of accurate performance monitoring, when no I/O should be done 3617 during the section of code that is being timed. 3618 3619 Level: intermediate 3620 3621 .keywords: SNES, get, convergence, history 3622 3623 .seealso: SNESSetConvergencHistory() 3624 3625 @*/ 3626 PetscErrorCode SNESGetConvergenceHistory(SNES snes,PetscReal *a[],PetscInt *its[],PetscInt *na) 3627 { 3628 PetscFunctionBegin; 3629 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3630 if (a) *a = snes->conv_hist; 3631 if (its) *its = snes->conv_hist_its; 3632 if (na) *na = snes->conv_hist_len; 3633 PetscFunctionReturn(0); 3634 } 3635 3636 #undef __FUNCT__ 3637 #define __FUNCT__ "SNESSetUpdate" 3638 /*@C 3639 SNESSetUpdate - Sets the general-purpose update function called 3640 at the beginning of every iteration of the nonlinear solve. Specifically 3641 it is called just before the Jacobian is "evaluated". 3642 3643 Logically Collective on SNES 3644 3645 Input Parameters: 3646 . snes - The nonlinear solver context 3647 . func - The function 3648 3649 Calling sequence of func: 3650 . func (SNES snes, PetscInt step); 3651 3652 . step - The current step of the iteration 3653 3654 Level: advanced 3655 3656 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() 3657 This is not used by most users. 3658 3659 .keywords: SNES, update 3660 3661 .seealso SNESSetJacobian(), SNESSolve() 3662 @*/ 3663 PetscErrorCode SNESSetUpdate(SNES snes, PetscErrorCode (*func)(SNES, PetscInt)) 3664 { 3665 PetscFunctionBegin; 3666 PetscValidHeaderSpecific(snes, SNES_CLASSID,1); 3667 snes->ops->update = func; 3668 PetscFunctionReturn(0); 3669 } 3670 3671 #undef __FUNCT__ 3672 #define __FUNCT__ "SNESScaleStep_Private" 3673 /* 3674 SNESScaleStep_Private - Scales a step so that its length is less than the 3675 positive parameter delta. 3676 3677 Input Parameters: 3678 + snes - the SNES context 3679 . y - approximate solution of linear system 3680 . fnorm - 2-norm of current function 3681 - delta - trust region size 3682 3683 Output Parameters: 3684 + gpnorm - predicted function norm at the new point, assuming local 3685 linearization. The value is zero if the step lies within the trust 3686 region, and exceeds zero otherwise. 3687 - ynorm - 2-norm of the step 3688 3689 Note: 3690 For non-trust region methods such as SNESNEWTONLS, the parameter delta 3691 is set to be the maximum allowable step size. 3692 3693 .keywords: SNES, nonlinear, scale, step 3694 */ 3695 PetscErrorCode SNESScaleStep_Private(SNES snes,Vec y,PetscReal *fnorm,PetscReal *delta,PetscReal *gpnorm,PetscReal *ynorm) 3696 { 3697 PetscReal nrm; 3698 PetscScalar cnorm; 3699 PetscErrorCode ierr; 3700 3701 PetscFunctionBegin; 3702 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3703 PetscValidHeaderSpecific(y,VEC_CLASSID,2); 3704 PetscCheckSameComm(snes,1,y,2); 3705 3706 ierr = VecNorm(y,NORM_2,&nrm);CHKERRQ(ierr); 3707 if (nrm > *delta) { 3708 nrm = *delta/nrm; 3709 *gpnorm = (1.0 - nrm)*(*fnorm); 3710 cnorm = nrm; 3711 ierr = VecScale(y,cnorm);CHKERRQ(ierr); 3712 *ynorm = *delta; 3713 } else { 3714 *gpnorm = 0.0; 3715 *ynorm = nrm; 3716 } 3717 PetscFunctionReturn(0); 3718 } 3719 3720 #undef __FUNCT__ 3721 #define __FUNCT__ "SNESSolve" 3722 /*@C 3723 SNESSolve - Solves a nonlinear system F(x) = b. 3724 Call SNESSolve() after calling SNESCreate() and optional routines of the form SNESSetXXX(). 3725 3726 Collective on SNES 3727 3728 Input Parameters: 3729 + snes - the SNES context 3730 . b - the constant part of the equation F(x) = b, or NULL to use zero. 3731 - x - the solution vector. 3732 3733 Notes: 3734 The user should initialize the vector,x, with the initial guess 3735 for the nonlinear solve prior to calling SNESSolve. In particular, 3736 to employ an initial guess of zero, the user should explicitly set 3737 this vector to zero by calling VecSet(). 3738 3739 Level: beginner 3740 3741 .keywords: SNES, nonlinear, solve 3742 3743 .seealso: SNESCreate(), SNESDestroy(), SNESSetFunction(), SNESSetJacobian(), SNESSetGridSequence(), SNESGetSolution() 3744 @*/ 3745 PetscErrorCode SNESSolve(SNES snes,Vec b,Vec x) 3746 { 3747 PetscErrorCode ierr; 3748 PetscBool flg; 3749 PetscViewer viewer; 3750 PetscInt grid; 3751 Vec xcreated = NULL; 3752 DM dm; 3753 PetscViewerFormat format; 3754 3755 PetscFunctionBegin; 3756 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3757 if (x) PetscValidHeaderSpecific(x,VEC_CLASSID,3); 3758 if (x) PetscCheckSameComm(snes,1,x,3); 3759 if (b) PetscValidHeaderSpecific(b,VEC_CLASSID,2); 3760 if (b) PetscCheckSameComm(snes,1,b,2); 3761 3762 if (!x) { 3763 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 3764 ierr = DMCreateGlobalVector(dm,&xcreated);CHKERRQ(ierr); 3765 x = xcreated; 3766 } 3767 3768 ierr = PetscOptionsGetViewer(PetscObjectComm((PetscObject)snes),((PetscObject)snes)->prefix,"-snes_view_pre",&viewer,&format,&flg);CHKERRQ(ierr); 3769 if (flg && !PetscPreLoadingOn) { 3770 ierr = PetscViewerPushFormat(viewer,format);CHKERRQ(ierr); 3771 ierr = SNESView(snes,viewer);CHKERRQ(ierr); 3772 ierr = PetscViewerPopFormat(viewer);CHKERRQ(ierr); 3773 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 3774 } 3775 3776 for (grid=0; grid<snes->gridsequence; grid++) {ierr = PetscViewerASCIIPushTab(PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)snes)));CHKERRQ(ierr);} 3777 for (grid=0; grid<snes->gridsequence+1; grid++) { 3778 3779 /* set solution vector */ 3780 if (!grid) {ierr = PetscObjectReference((PetscObject)x);CHKERRQ(ierr);} 3781 ierr = VecDestroy(&snes->vec_sol);CHKERRQ(ierr); 3782 snes->vec_sol = x; 3783 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 3784 3785 /* set affine vector if provided */ 3786 if (b) { ierr = PetscObjectReference((PetscObject)b);CHKERRQ(ierr); } 3787 ierr = VecDestroy(&snes->vec_rhs);CHKERRQ(ierr); 3788 snes->vec_rhs = b; 3789 3790 if (snes->vec_func == snes->vec_sol) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_IDN,"Solution vector cannot be function vector"); 3791 if (snes->vec_rhs == snes->vec_sol) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_IDN,"Solution vector cannot be right hand side vector"); 3792 if (!snes->vec_sol_update /* && snes->vec_sol */) { 3793 ierr = VecDuplicate(snes->vec_sol,&snes->vec_sol_update);CHKERRQ(ierr); 3794 ierr = PetscLogObjectParent((PetscObject)snes,(PetscObject)snes->vec_sol_update);CHKERRQ(ierr); 3795 } 3796 ierr = DMShellSetGlobalVector(dm,snes->vec_sol);CHKERRQ(ierr); 3797 ierr = SNESSetUp(snes);CHKERRQ(ierr); 3798 3799 if (!grid) { 3800 if (snes->ops->computeinitialguess) { 3801 ierr = (*snes->ops->computeinitialguess)(snes,snes->vec_sol,snes->initialguessP);CHKERRQ(ierr); 3802 } 3803 } 3804 3805 if (snes->conv_hist_reset) snes->conv_hist_len = 0; 3806 if (snes->counters_reset) {snes->nfuncs = 0; snes->linear_its = 0; snes->numFailures = 0;} 3807 3808 ierr = PetscLogEventBegin(SNES_Solve,snes,0,0,0);CHKERRQ(ierr); 3809 ierr = (*snes->ops->solve)(snes);CHKERRQ(ierr); 3810 ierr = PetscLogEventEnd(SNES_Solve,snes,0,0,0);CHKERRQ(ierr); 3811 if (snes->domainerror) { 3812 snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; 3813 snes->domainerror = PETSC_FALSE; 3814 } 3815 if (!snes->reason) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Internal error, solver returned without setting converged reason"); 3816 3817 if (snes->lagjac_persist) snes->jac_iter += snes->iter; 3818 if (snes->lagpre_persist) snes->pre_iter += snes->iter; 3819 3820 flg = PETSC_FALSE; 3821 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_test_local_min",&flg,NULL);CHKERRQ(ierr); 3822 if (flg && !PetscPreLoadingOn) { ierr = SNESTestLocalMin(snes);CHKERRQ(ierr); } 3823 if (snes->printreason) { 3824 ierr = PetscViewerASCIIAddTab(PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)snes)),((PetscObject)snes)->tablevel);CHKERRQ(ierr); 3825 if (snes->reason > 0) { 3826 ierr = PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)snes)),"Nonlinear solve converged due to %s iterations %D\n",SNESConvergedReasons[snes->reason],snes->iter);CHKERRQ(ierr); 3827 } else { 3828 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); 3829 } 3830 ierr = PetscViewerASCIISubtractTab(PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)snes)),((PetscObject)snes)->tablevel);CHKERRQ(ierr); 3831 } 3832 3833 if (snes->errorifnotconverged && snes->reason < 0) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_NOT_CONVERGED,"SNESSolve has not converged"); 3834 if (grid < snes->gridsequence) { 3835 DM fine; 3836 Vec xnew; 3837 Mat interp; 3838 3839 ierr = DMRefine(snes->dm,PetscObjectComm((PetscObject)snes),&fine);CHKERRQ(ierr); 3840 if (!fine) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_INCOMP,"DMRefine() did not perform any refinement, cannot continue grid sequencing"); 3841 ierr = DMCreateInterpolation(snes->dm,fine,&interp,NULL);CHKERRQ(ierr); 3842 ierr = DMCreateGlobalVector(fine,&xnew);CHKERRQ(ierr); 3843 ierr = MatInterpolate(interp,x,xnew);CHKERRQ(ierr); 3844 ierr = DMInterpolate(snes->dm,interp,fine);CHKERRQ(ierr); 3845 ierr = MatDestroy(&interp);CHKERRQ(ierr); 3846 x = xnew; 3847 3848 ierr = SNESReset(snes);CHKERRQ(ierr); 3849 ierr = SNESSetDM(snes,fine);CHKERRQ(ierr); 3850 ierr = DMDestroy(&fine);CHKERRQ(ierr); 3851 ierr = PetscViewerASCIIPopTab(PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)snes)));CHKERRQ(ierr); 3852 } 3853 } 3854 /* monitoring and viewing */ 3855 ierr = PetscOptionsGetViewer(PetscObjectComm((PetscObject)snes),((PetscObject)snes)->prefix,"-snes_view",&viewer,&format,&flg);CHKERRQ(ierr); 3856 if (flg && !PetscPreLoadingOn) { 3857 ierr = PetscViewerPushFormat(viewer,format);CHKERRQ(ierr); 3858 ierr = SNESView(snes,viewer);CHKERRQ(ierr); 3859 ierr = PetscViewerPopFormat(viewer);CHKERRQ(ierr); 3860 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 3861 } 3862 ierr = VecViewFromOptions(snes->vec_sol,((PetscObject)snes)->prefix,"-snes_view_solution");CHKERRQ(ierr); 3863 3864 ierr = VecDestroy(&xcreated);CHKERRQ(ierr); 3865 ierr = PetscObjectSAWsBlock((PetscObject)snes);CHKERRQ(ierr); 3866 PetscFunctionReturn(0); 3867 } 3868 3869 /* --------- Internal routines for SNES Package --------- */ 3870 3871 #undef __FUNCT__ 3872 #define __FUNCT__ "SNESSetType" 3873 /*@C 3874 SNESSetType - Sets the method for the nonlinear solver. 3875 3876 Collective on SNES 3877 3878 Input Parameters: 3879 + snes - the SNES context 3880 - type - a known method 3881 3882 Options Database Key: 3883 . -snes_type <type> - Sets the method; use -help for a list 3884 of available methods (for instance, newtonls or newtontr) 3885 3886 Notes: 3887 See "petsc/include/petscsnes.h" for available methods (for instance) 3888 + SNESNEWTONLS - Newton's method with line search 3889 (systems of nonlinear equations) 3890 . SNESNEWTONTR - Newton's method with trust region 3891 (systems of nonlinear equations) 3892 3893 Normally, it is best to use the SNESSetFromOptions() command and then 3894 set the SNES solver type from the options database rather than by using 3895 this routine. Using the options database provides the user with 3896 maximum flexibility in evaluating the many nonlinear solvers. 3897 The SNESSetType() routine is provided for those situations where it 3898 is necessary to set the nonlinear solver independently of the command 3899 line or options database. This might be the case, for example, when 3900 the choice of solver changes during the execution of the program, 3901 and the user's application is taking responsibility for choosing the 3902 appropriate method. 3903 3904 Developer Notes: SNESRegister() adds a constructor for a new SNESType to SNESList, SNESSetType() locates 3905 the constructor in that list and calls it to create the spexific object. 3906 3907 Level: intermediate 3908 3909 .keywords: SNES, set, type 3910 3911 .seealso: SNESType, SNESCreate(), SNESDestroy(), SNESGetType(), SNESSetFromOptions() 3912 3913 @*/ 3914 PetscErrorCode SNESSetType(SNES snes,SNESType type) 3915 { 3916 PetscErrorCode ierr,(*r)(SNES); 3917 PetscBool match; 3918 3919 PetscFunctionBegin; 3920 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3921 PetscValidCharPointer(type,2); 3922 3923 ierr = PetscObjectTypeCompare((PetscObject)snes,type,&match);CHKERRQ(ierr); 3924 if (match) PetscFunctionReturn(0); 3925 3926 ierr = PetscFunctionListFind(SNESList,type,&r);CHKERRQ(ierr); 3927 if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unable to find requested SNES type %s",type); 3928 /* Destroy the previous private SNES context */ 3929 if (snes->ops->destroy) { 3930 ierr = (*(snes)->ops->destroy)(snes);CHKERRQ(ierr); 3931 snes->ops->destroy = NULL; 3932 } 3933 /* Reinitialize function pointers in SNESOps structure */ 3934 snes->ops->setup = 0; 3935 snes->ops->solve = 0; 3936 snes->ops->view = 0; 3937 snes->ops->setfromoptions = 0; 3938 snes->ops->destroy = 0; 3939 /* Call the SNESCreate_XXX routine for this particular Nonlinear solver */ 3940 snes->setupcalled = PETSC_FALSE; 3941 3942 ierr = PetscObjectChangeTypeName((PetscObject)snes,type);CHKERRQ(ierr); 3943 ierr = (*r)(snes);CHKERRQ(ierr); 3944 PetscFunctionReturn(0); 3945 } 3946 3947 #undef __FUNCT__ 3948 #define __FUNCT__ "SNESGetType" 3949 /*@C 3950 SNESGetType - Gets the SNES method type and name (as a string). 3951 3952 Not Collective 3953 3954 Input Parameter: 3955 . snes - nonlinear solver context 3956 3957 Output Parameter: 3958 . type - SNES method (a character string) 3959 3960 Level: intermediate 3961 3962 .keywords: SNES, nonlinear, get, type, name 3963 @*/ 3964 PetscErrorCode SNESGetType(SNES snes,SNESType *type) 3965 { 3966 PetscFunctionBegin; 3967 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3968 PetscValidPointer(type,2); 3969 *type = ((PetscObject)snes)->type_name; 3970 PetscFunctionReturn(0); 3971 } 3972 3973 #undef __FUNCT__ 3974 #define __FUNCT__ "SNESGetSolution" 3975 /*@ 3976 SNESGetSolution - Returns the vector where the approximate solution is 3977 stored. This is the fine grid solution when using SNESSetGridSequence(). 3978 3979 Not Collective, but Vec is parallel if SNES is parallel 3980 3981 Input Parameter: 3982 . snes - the SNES context 3983 3984 Output Parameter: 3985 . x - the solution 3986 3987 Level: intermediate 3988 3989 .keywords: SNES, nonlinear, get, solution 3990 3991 .seealso: SNESGetSolutionUpdate(), SNESGetFunction() 3992 @*/ 3993 PetscErrorCode SNESGetSolution(SNES snes,Vec *x) 3994 { 3995 PetscFunctionBegin; 3996 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3997 PetscValidPointer(x,2); 3998 *x = snes->vec_sol; 3999 PetscFunctionReturn(0); 4000 } 4001 4002 #undef __FUNCT__ 4003 #define __FUNCT__ "SNESGetSolutionUpdate" 4004 /*@ 4005 SNESGetSolutionUpdate - Returns the vector where the solution update is 4006 stored. 4007 4008 Not Collective, but Vec is parallel if SNES is parallel 4009 4010 Input Parameter: 4011 . snes - the SNES context 4012 4013 Output Parameter: 4014 . x - the solution update 4015 4016 Level: advanced 4017 4018 .keywords: SNES, nonlinear, get, solution, update 4019 4020 .seealso: SNESGetSolution(), SNESGetFunction() 4021 @*/ 4022 PetscErrorCode SNESGetSolutionUpdate(SNES snes,Vec *x) 4023 { 4024 PetscFunctionBegin; 4025 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4026 PetscValidPointer(x,2); 4027 *x = snes->vec_sol_update; 4028 PetscFunctionReturn(0); 4029 } 4030 4031 #undef __FUNCT__ 4032 #define __FUNCT__ "SNESGetFunction" 4033 /*@C 4034 SNESGetFunction - Returns the vector where the function is stored. 4035 4036 Not Collective, but Vec is parallel if SNES is parallel. Collective if Vec is requested, but has not been created yet. 4037 4038 Input Parameter: 4039 . snes - the SNES context 4040 4041 Output Parameter: 4042 + r - the vector that is used to store residuals (or NULL if you don't want it) 4043 . SNESFunction- the function (or NULL if you don't want it) 4044 - ctx - the function context (or NULL if you don't want it) 4045 4046 Level: advanced 4047 4048 .keywords: SNES, nonlinear, get, function 4049 4050 .seealso: SNESSetFunction(), SNESGetSolution(), SNESFunction 4051 @*/ 4052 PetscErrorCode SNESGetFunction(SNES snes,Vec *r,PetscErrorCode (**SNESFunction)(SNES,Vec,Vec,void*),void **ctx) 4053 { 4054 PetscErrorCode ierr; 4055 DM dm; 4056 4057 PetscFunctionBegin; 4058 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4059 if (r) { 4060 if (!snes->vec_func) { 4061 if (snes->vec_rhs) { 4062 ierr = VecDuplicate(snes->vec_rhs,&snes->vec_func);CHKERRQ(ierr); 4063 } else if (snes->vec_sol) { 4064 ierr = VecDuplicate(snes->vec_sol,&snes->vec_func);CHKERRQ(ierr); 4065 } else if (snes->dm) { 4066 ierr = DMCreateGlobalVector(snes->dm,&snes->vec_func);CHKERRQ(ierr); 4067 } 4068 } 4069 *r = snes->vec_func; 4070 } 4071 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 4072 ierr = DMSNESGetFunction(dm,SNESFunction,ctx);CHKERRQ(ierr); 4073 PetscFunctionReturn(0); 4074 } 4075 4076 /*@C 4077 SNESGetGS - Returns the GS function and context. 4078 4079 Input Parameter: 4080 . snes - the SNES context 4081 4082 Output Parameter: 4083 + SNESGSFunction - the function (or NULL) 4084 - ctx - the function context (or NULL) 4085 4086 Level: advanced 4087 4088 .keywords: SNES, nonlinear, get, function 4089 4090 .seealso: SNESSetGS(), SNESGetFunction() 4091 @*/ 4092 4093 #undef __FUNCT__ 4094 #define __FUNCT__ "SNESGetGS" 4095 PetscErrorCode SNESGetGS (SNES snes, PetscErrorCode (**SNESGSFunction)(SNES, Vec, Vec, void*), void ** ctx) 4096 { 4097 PetscErrorCode ierr; 4098 DM dm; 4099 4100 PetscFunctionBegin; 4101 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4102 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 4103 ierr = DMSNESGetGS(dm,SNESGSFunction,ctx);CHKERRQ(ierr); 4104 PetscFunctionReturn(0); 4105 } 4106 4107 #undef __FUNCT__ 4108 #define __FUNCT__ "SNESSetOptionsPrefix" 4109 /*@C 4110 SNESSetOptionsPrefix - Sets the prefix used for searching for all 4111 SNES options in the database. 4112 4113 Logically Collective on SNES 4114 4115 Input Parameter: 4116 + snes - the SNES context 4117 - prefix - the prefix to prepend to all option names 4118 4119 Notes: 4120 A hyphen (-) must NOT be given at the beginning of the prefix name. 4121 The first character of all runtime options is AUTOMATICALLY the hyphen. 4122 4123 Level: advanced 4124 4125 .keywords: SNES, set, options, prefix, database 4126 4127 .seealso: SNESSetFromOptions() 4128 @*/ 4129 PetscErrorCode SNESSetOptionsPrefix(SNES snes,const char prefix[]) 4130 { 4131 PetscErrorCode ierr; 4132 4133 PetscFunctionBegin; 4134 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4135 ierr = PetscObjectSetOptionsPrefix((PetscObject)snes,prefix);CHKERRQ(ierr); 4136 if (!snes->ksp) {ierr = SNESGetKSP(snes,&snes->ksp);CHKERRQ(ierr);} 4137 if (snes->linesearch) { 4138 ierr = SNESGetLineSearch(snes,&snes->linesearch);CHKERRQ(ierr); 4139 ierr = PetscObjectSetOptionsPrefix((PetscObject)snes->linesearch,prefix);CHKERRQ(ierr); 4140 } 4141 ierr = KSPSetOptionsPrefix(snes->ksp,prefix);CHKERRQ(ierr); 4142 PetscFunctionReturn(0); 4143 } 4144 4145 #undef __FUNCT__ 4146 #define __FUNCT__ "SNESAppendOptionsPrefix" 4147 /*@C 4148 SNESAppendOptionsPrefix - Appends to the prefix used for searching for all 4149 SNES options in the database. 4150 4151 Logically Collective on SNES 4152 4153 Input Parameters: 4154 + snes - the SNES context 4155 - prefix - the prefix to prepend to all option names 4156 4157 Notes: 4158 A hyphen (-) must NOT be given at the beginning of the prefix name. 4159 The first character of all runtime options is AUTOMATICALLY the hyphen. 4160 4161 Level: advanced 4162 4163 .keywords: SNES, append, options, prefix, database 4164 4165 .seealso: SNESGetOptionsPrefix() 4166 @*/ 4167 PetscErrorCode SNESAppendOptionsPrefix(SNES snes,const char prefix[]) 4168 { 4169 PetscErrorCode ierr; 4170 4171 PetscFunctionBegin; 4172 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4173 ierr = PetscObjectAppendOptionsPrefix((PetscObject)snes,prefix);CHKERRQ(ierr); 4174 if (!snes->ksp) {ierr = SNESGetKSP(snes,&snes->ksp);CHKERRQ(ierr);} 4175 if (snes->linesearch) { 4176 ierr = SNESGetLineSearch(snes,&snes->linesearch);CHKERRQ(ierr); 4177 ierr = PetscObjectAppendOptionsPrefix((PetscObject)snes->linesearch,prefix);CHKERRQ(ierr); 4178 } 4179 ierr = KSPAppendOptionsPrefix(snes->ksp,prefix);CHKERRQ(ierr); 4180 PetscFunctionReturn(0); 4181 } 4182 4183 #undef __FUNCT__ 4184 #define __FUNCT__ "SNESGetOptionsPrefix" 4185 /*@C 4186 SNESGetOptionsPrefix - Sets the prefix used for searching for all 4187 SNES options in the database. 4188 4189 Not Collective 4190 4191 Input Parameter: 4192 . snes - the SNES context 4193 4194 Output Parameter: 4195 . prefix - pointer to the prefix string used 4196 4197 Notes: On the fortran side, the user should pass in a string 'prefix' of 4198 sufficient length to hold the prefix. 4199 4200 Level: advanced 4201 4202 .keywords: SNES, get, options, prefix, database 4203 4204 .seealso: SNESAppendOptionsPrefix() 4205 @*/ 4206 PetscErrorCode SNESGetOptionsPrefix(SNES snes,const char *prefix[]) 4207 { 4208 PetscErrorCode ierr; 4209 4210 PetscFunctionBegin; 4211 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4212 ierr = PetscObjectGetOptionsPrefix((PetscObject)snes,prefix);CHKERRQ(ierr); 4213 PetscFunctionReturn(0); 4214 } 4215 4216 4217 #undef __FUNCT__ 4218 #define __FUNCT__ "SNESRegister" 4219 /*@C 4220 SNESRegister - Adds a method to the nonlinear solver package. 4221 4222 Not collective 4223 4224 Input Parameters: 4225 + name_solver - name of a new user-defined solver 4226 - routine_create - routine to create method context 4227 4228 Notes: 4229 SNESRegister() may be called multiple times to add several user-defined solvers. 4230 4231 Sample usage: 4232 .vb 4233 SNESRegister("my_solver",MySolverCreate); 4234 .ve 4235 4236 Then, your solver can be chosen with the procedural interface via 4237 $ SNESSetType(snes,"my_solver") 4238 or at runtime via the option 4239 $ -snes_type my_solver 4240 4241 Level: advanced 4242 4243 Note: If your function is not being put into a shared library then use SNESRegister() instead 4244 4245 .keywords: SNES, nonlinear, register 4246 4247 .seealso: SNESRegisterAll(), SNESRegisterDestroy() 4248 4249 Level: advanced 4250 @*/ 4251 PetscErrorCode SNESRegister(const char sname[],PetscErrorCode (*function)(SNES)) 4252 { 4253 PetscErrorCode ierr; 4254 4255 PetscFunctionBegin; 4256 ierr = PetscFunctionListAdd(&SNESList,sname,function);CHKERRQ(ierr); 4257 PetscFunctionReturn(0); 4258 } 4259 4260 #undef __FUNCT__ 4261 #define __FUNCT__ "SNESTestLocalMin" 4262 PetscErrorCode SNESTestLocalMin(SNES snes) 4263 { 4264 PetscErrorCode ierr; 4265 PetscInt N,i,j; 4266 Vec u,uh,fh; 4267 PetscScalar value; 4268 PetscReal norm; 4269 4270 PetscFunctionBegin; 4271 ierr = SNESGetSolution(snes,&u);CHKERRQ(ierr); 4272 ierr = VecDuplicate(u,&uh);CHKERRQ(ierr); 4273 ierr = VecDuplicate(u,&fh);CHKERRQ(ierr); 4274 4275 /* currently only works for sequential */ 4276 ierr = PetscPrintf(PETSC_COMM_WORLD,"Testing FormFunction() for local min\n");CHKERRQ(ierr); 4277 ierr = VecGetSize(u,&N);CHKERRQ(ierr); 4278 for (i=0; i<N; i++) { 4279 ierr = VecCopy(u,uh);CHKERRQ(ierr); 4280 ierr = PetscPrintf(PETSC_COMM_WORLD,"i = %D\n",i);CHKERRQ(ierr); 4281 for (j=-10; j<11; j++) { 4282 value = PetscSign(j)*exp(PetscAbs(j)-10.0); 4283 ierr = VecSetValue(uh,i,value,ADD_VALUES);CHKERRQ(ierr); 4284 ierr = SNESComputeFunction(snes,uh,fh);CHKERRQ(ierr); 4285 ierr = VecNorm(fh,NORM_2,&norm);CHKERRQ(ierr); 4286 ierr = PetscPrintf(PETSC_COMM_WORLD," j norm %D %18.16e\n",j,norm);CHKERRQ(ierr); 4287 value = -value; 4288 ierr = VecSetValue(uh,i,value,ADD_VALUES);CHKERRQ(ierr); 4289 } 4290 } 4291 ierr = VecDestroy(&uh);CHKERRQ(ierr); 4292 ierr = VecDestroy(&fh);CHKERRQ(ierr); 4293 PetscFunctionReturn(0); 4294 } 4295 4296 #undef __FUNCT__ 4297 #define __FUNCT__ "SNESKSPSetUseEW" 4298 /*@ 4299 SNESKSPSetUseEW - Sets SNES use Eisenstat-Walker method for 4300 computing relative tolerance for linear solvers within an inexact 4301 Newton method. 4302 4303 Logically Collective on SNES 4304 4305 Input Parameters: 4306 + snes - SNES context 4307 - flag - PETSC_TRUE or PETSC_FALSE 4308 4309 Options Database: 4310 + -snes_ksp_ew - use Eisenstat-Walker method for determining linear system convergence 4311 . -snes_ksp_ew_version ver - version of Eisenstat-Walker method 4312 . -snes_ksp_ew_rtol0 <rtol0> - Sets rtol0 4313 . -snes_ksp_ew_rtolmax <rtolmax> - Sets rtolmax 4314 . -snes_ksp_ew_gamma <gamma> - Sets gamma 4315 . -snes_ksp_ew_alpha <alpha> - Sets alpha 4316 . -snes_ksp_ew_alpha2 <alpha2> - Sets alpha2 4317 - -snes_ksp_ew_threshold <threshold> - Sets threshold 4318 4319 Notes: 4320 Currently, the default is to use a constant relative tolerance for 4321 the inner linear solvers. Alternatively, one can use the 4322 Eisenstat-Walker method, where the relative convergence tolerance 4323 is reset at each Newton iteration according progress of the nonlinear 4324 solver. 4325 4326 Level: advanced 4327 4328 Reference: 4329 S. C. Eisenstat and H. F. Walker, "Choosing the forcing terms in an 4330 inexact Newton method", SISC 17 (1), pp.16-32, 1996. 4331 4332 .keywords: SNES, KSP, Eisenstat, Walker, convergence, test, inexact, Newton 4333 4334 .seealso: SNESKSPGetUseEW(), SNESKSPGetParametersEW(), SNESKSPSetParametersEW() 4335 @*/ 4336 PetscErrorCode SNESKSPSetUseEW(SNES snes,PetscBool flag) 4337 { 4338 PetscFunctionBegin; 4339 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4340 PetscValidLogicalCollectiveBool(snes,flag,2); 4341 snes->ksp_ewconv = flag; 4342 PetscFunctionReturn(0); 4343 } 4344 4345 #undef __FUNCT__ 4346 #define __FUNCT__ "SNESKSPGetUseEW" 4347 /*@ 4348 SNESKSPGetUseEW - Gets if SNES is using Eisenstat-Walker method 4349 for computing relative tolerance for linear solvers within an 4350 inexact Newton method. 4351 4352 Not Collective 4353 4354 Input Parameter: 4355 . snes - SNES context 4356 4357 Output Parameter: 4358 . flag - PETSC_TRUE or PETSC_FALSE 4359 4360 Notes: 4361 Currently, the default is to use a constant relative tolerance for 4362 the inner linear solvers. Alternatively, one can use the 4363 Eisenstat-Walker method, where the relative convergence tolerance 4364 is reset at each Newton iteration according progress of the nonlinear 4365 solver. 4366 4367 Level: advanced 4368 4369 Reference: 4370 S. C. Eisenstat and H. F. Walker, "Choosing the forcing terms in an 4371 inexact Newton method", SISC 17 (1), pp.16-32, 1996. 4372 4373 .keywords: SNES, KSP, Eisenstat, Walker, convergence, test, inexact, Newton 4374 4375 .seealso: SNESKSPSetUseEW(), SNESKSPGetParametersEW(), SNESKSPSetParametersEW() 4376 @*/ 4377 PetscErrorCode SNESKSPGetUseEW(SNES snes, PetscBool *flag) 4378 { 4379 PetscFunctionBegin; 4380 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4381 PetscValidPointer(flag,2); 4382 *flag = snes->ksp_ewconv; 4383 PetscFunctionReturn(0); 4384 } 4385 4386 #undef __FUNCT__ 4387 #define __FUNCT__ "SNESKSPSetParametersEW" 4388 /*@ 4389 SNESKSPSetParametersEW - Sets parameters for Eisenstat-Walker 4390 convergence criteria for the linear solvers within an inexact 4391 Newton method. 4392 4393 Logically Collective on SNES 4394 4395 Input Parameters: 4396 + snes - SNES context 4397 . version - version 1, 2 (default is 2) or 3 4398 . rtol_0 - initial relative tolerance (0 <= rtol_0 < 1) 4399 . rtol_max - maximum relative tolerance (0 <= rtol_max < 1) 4400 . gamma - multiplicative factor for version 2 rtol computation 4401 (0 <= gamma2 <= 1) 4402 . alpha - power for version 2 rtol computation (1 < alpha <= 2) 4403 . alpha2 - power for safeguard 4404 - threshold - threshold for imposing safeguard (0 < threshold < 1) 4405 4406 Note: 4407 Version 3 was contributed by Luis Chacon, June 2006. 4408 4409 Use PETSC_DEFAULT to retain the default for any of the parameters. 4410 4411 Level: advanced 4412 4413 Reference: 4414 S. C. Eisenstat and H. F. Walker, "Choosing the forcing terms in an 4415 inexact Newton method", Utah State University Math. Stat. Dept. Res. 4416 Report 6/94/75, June, 1994, to appear in SIAM J. Sci. Comput. 4417 4418 .keywords: SNES, KSP, Eisenstat, Walker, set, parameters 4419 4420 .seealso: SNESKSPSetUseEW(), SNESKSPGetUseEW(), SNESKSPGetParametersEW() 4421 @*/ 4422 PetscErrorCode SNESKSPSetParametersEW(SNES snes,PetscInt version,PetscReal rtol_0,PetscReal rtol_max,PetscReal gamma,PetscReal alpha,PetscReal alpha2,PetscReal threshold) 4423 { 4424 SNESKSPEW *kctx; 4425 4426 PetscFunctionBegin; 4427 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4428 kctx = (SNESKSPEW*)snes->kspconvctx; 4429 if (!kctx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"No Eisenstat-Walker context existing"); 4430 PetscValidLogicalCollectiveInt(snes,version,2); 4431 PetscValidLogicalCollectiveReal(snes,rtol_0,3); 4432 PetscValidLogicalCollectiveReal(snes,rtol_max,4); 4433 PetscValidLogicalCollectiveReal(snes,gamma,5); 4434 PetscValidLogicalCollectiveReal(snes,alpha,6); 4435 PetscValidLogicalCollectiveReal(snes,alpha2,7); 4436 PetscValidLogicalCollectiveReal(snes,threshold,8); 4437 4438 if (version != PETSC_DEFAULT) kctx->version = version; 4439 if (rtol_0 != PETSC_DEFAULT) kctx->rtol_0 = rtol_0; 4440 if (rtol_max != PETSC_DEFAULT) kctx->rtol_max = rtol_max; 4441 if (gamma != PETSC_DEFAULT) kctx->gamma = gamma; 4442 if (alpha != PETSC_DEFAULT) kctx->alpha = alpha; 4443 if (alpha2 != PETSC_DEFAULT) kctx->alpha2 = alpha2; 4444 if (threshold != PETSC_DEFAULT) kctx->threshold = threshold; 4445 4446 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); 4447 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); 4448 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); 4449 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); 4450 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); 4451 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); 4452 PetscFunctionReturn(0); 4453 } 4454 4455 #undef __FUNCT__ 4456 #define __FUNCT__ "SNESKSPGetParametersEW" 4457 /*@ 4458 SNESKSPGetParametersEW - Gets parameters for Eisenstat-Walker 4459 convergence criteria for the linear solvers within an inexact 4460 Newton method. 4461 4462 Not Collective 4463 4464 Input Parameters: 4465 snes - SNES context 4466 4467 Output Parameters: 4468 + version - version 1, 2 (default is 2) or 3 4469 . rtol_0 - initial relative tolerance (0 <= rtol_0 < 1) 4470 . rtol_max - maximum relative tolerance (0 <= rtol_max < 1) 4471 . gamma - multiplicative factor for version 2 rtol computation (0 <= gamma2 <= 1) 4472 . alpha - power for version 2 rtol computation (1 < alpha <= 2) 4473 . alpha2 - power for safeguard 4474 - threshold - threshold for imposing safeguard (0 < threshold < 1) 4475 4476 Level: advanced 4477 4478 .keywords: SNES, KSP, Eisenstat, Walker, get, parameters 4479 4480 .seealso: SNESKSPSetUseEW(), SNESKSPGetUseEW(), SNESKSPSetParametersEW() 4481 @*/ 4482 PetscErrorCode SNESKSPGetParametersEW(SNES snes,PetscInt *version,PetscReal *rtol_0,PetscReal *rtol_max,PetscReal *gamma,PetscReal *alpha,PetscReal *alpha2,PetscReal *threshold) 4483 { 4484 SNESKSPEW *kctx; 4485 4486 PetscFunctionBegin; 4487 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4488 kctx = (SNESKSPEW*)snes->kspconvctx; 4489 if (!kctx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"No Eisenstat-Walker context existing"); 4490 if (version) *version = kctx->version; 4491 if (rtol_0) *rtol_0 = kctx->rtol_0; 4492 if (rtol_max) *rtol_max = kctx->rtol_max; 4493 if (gamma) *gamma = kctx->gamma; 4494 if (alpha) *alpha = kctx->alpha; 4495 if (alpha2) *alpha2 = kctx->alpha2; 4496 if (threshold) *threshold = kctx->threshold; 4497 PetscFunctionReturn(0); 4498 } 4499 4500 #undef __FUNCT__ 4501 #define __FUNCT__ "KSPPreSolve_SNESEW" 4502 PetscErrorCode KSPPreSolve_SNESEW(KSP ksp, Vec b, Vec x, SNES snes) 4503 { 4504 PetscErrorCode ierr; 4505 SNESKSPEW *kctx = (SNESKSPEW*)snes->kspconvctx; 4506 PetscReal rtol = PETSC_DEFAULT,stol; 4507 4508 PetscFunctionBegin; 4509 if (!snes->ksp_ewconv) PetscFunctionReturn(0); 4510 if (!snes->iter) { 4511 rtol = kctx->rtol_0; /* first time in, so use the original user rtol */ 4512 ierr = VecNorm(snes->vec_func,NORM_2,&kctx->norm_first);CHKERRQ(ierr); 4513 } 4514 else { 4515 if (kctx->version == 1) { 4516 rtol = (snes->norm - kctx->lresid_last)/kctx->norm_last; 4517 if (rtol < 0.0) rtol = -rtol; 4518 stol = PetscPowReal(kctx->rtol_last,kctx->alpha2); 4519 if (stol > kctx->threshold) rtol = PetscMax(rtol,stol); 4520 } else if (kctx->version == 2) { 4521 rtol = kctx->gamma * PetscPowReal(snes->norm/kctx->norm_last,kctx->alpha); 4522 stol = kctx->gamma * PetscPowReal(kctx->rtol_last,kctx->alpha); 4523 if (stol > kctx->threshold) rtol = PetscMax(rtol,stol); 4524 } else if (kctx->version == 3) { /* contributed by Luis Chacon, June 2006. */ 4525 rtol = kctx->gamma * PetscPowReal(snes->norm/kctx->norm_last,kctx->alpha); 4526 /* safeguard: avoid sharp decrease of rtol */ 4527 stol = kctx->gamma*PetscPowReal(kctx->rtol_last,kctx->alpha); 4528 stol = PetscMax(rtol,stol); 4529 rtol = PetscMin(kctx->rtol_0,stol); 4530 /* safeguard: avoid oversolving */ 4531 stol = kctx->gamma*(kctx->norm_first*snes->rtol)/snes->norm; 4532 stol = PetscMax(rtol,stol); 4533 rtol = PetscMin(kctx->rtol_0,stol); 4534 } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Only versions 1, 2 or 3 are supported: %D",kctx->version); 4535 } 4536 /* safeguard: avoid rtol greater than one */ 4537 rtol = PetscMin(rtol,kctx->rtol_max); 4538 ierr = KSPSetTolerances(ksp,rtol,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);CHKERRQ(ierr); 4539 ierr = PetscInfo3(snes,"iter %D, Eisenstat-Walker (version %D) KSP rtol=%G\n",snes->iter,kctx->version,rtol);CHKERRQ(ierr); 4540 PetscFunctionReturn(0); 4541 } 4542 4543 #undef __FUNCT__ 4544 #define __FUNCT__ "KSPPostSolve_SNESEW" 4545 PetscErrorCode KSPPostSolve_SNESEW(KSP ksp, Vec b, Vec x, SNES snes) 4546 { 4547 PetscErrorCode ierr; 4548 SNESKSPEW *kctx = (SNESKSPEW*)snes->kspconvctx; 4549 PCSide pcside; 4550 Vec lres; 4551 4552 PetscFunctionBegin; 4553 if (!snes->ksp_ewconv) PetscFunctionReturn(0); 4554 ierr = KSPGetTolerances(ksp,&kctx->rtol_last,0,0,0);CHKERRQ(ierr); 4555 ierr = SNESGetFunctionNorm(snes,&kctx->norm_last);CHKERRQ(ierr); 4556 if (kctx->version == 1) { 4557 ierr = KSPGetPCSide(ksp,&pcside);CHKERRQ(ierr); 4558 if (pcside == PC_RIGHT) { /* XXX Should we also test KSP_UNPRECONDITIONED_NORM ? */ 4559 /* KSP residual is true linear residual */ 4560 ierr = KSPGetResidualNorm(ksp,&kctx->lresid_last);CHKERRQ(ierr); 4561 } else { 4562 /* KSP residual is preconditioned residual */ 4563 /* compute true linear residual norm */ 4564 ierr = VecDuplicate(b,&lres);CHKERRQ(ierr); 4565 ierr = MatMult(snes->jacobian,x,lres);CHKERRQ(ierr); 4566 ierr = VecAYPX(lres,-1.0,b);CHKERRQ(ierr); 4567 ierr = VecNorm(lres,NORM_2,&kctx->lresid_last);CHKERRQ(ierr); 4568 ierr = VecDestroy(&lres);CHKERRQ(ierr); 4569 } 4570 } 4571 PetscFunctionReturn(0); 4572 } 4573 4574 #undef __FUNCT__ 4575 #define __FUNCT__ "SNESGetKSP" 4576 /*@ 4577 SNESGetKSP - Returns the KSP context for a SNES solver. 4578 4579 Not Collective, but if SNES object is parallel, then KSP object is parallel 4580 4581 Input Parameter: 4582 . snes - the SNES context 4583 4584 Output Parameter: 4585 . ksp - the KSP context 4586 4587 Notes: 4588 The user can then directly manipulate the KSP context to set various 4589 options, etc. Likewise, the user can then extract and manipulate the 4590 PC contexts as well. 4591 4592 Level: beginner 4593 4594 .keywords: SNES, nonlinear, get, KSP, context 4595 4596 .seealso: KSPGetPC(), SNESCreate(), KSPCreate(), SNESSetKSP() 4597 @*/ 4598 PetscErrorCode SNESGetKSP(SNES snes,KSP *ksp) 4599 { 4600 PetscErrorCode ierr; 4601 4602 PetscFunctionBegin; 4603 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4604 PetscValidPointer(ksp,2); 4605 4606 if (!snes->ksp) { 4607 ierr = KSPCreate(PetscObjectComm((PetscObject)snes),&snes->ksp);CHKERRQ(ierr); 4608 ierr = PetscObjectIncrementTabLevel((PetscObject)snes->ksp,(PetscObject)snes,1);CHKERRQ(ierr); 4609 ierr = PetscLogObjectParent((PetscObject)snes,(PetscObject)snes->ksp);CHKERRQ(ierr); 4610 4611 ierr = KSPSetPreSolve(snes->ksp,(PetscErrorCode (*)(KSP,Vec,Vec,void*))KSPPreSolve_SNESEW,snes);CHKERRQ(ierr); 4612 ierr = KSPSetPostSolve(snes->ksp,(PetscErrorCode (*)(KSP,Vec,Vec,void*))KSPPostSolve_SNESEW,snes);CHKERRQ(ierr); 4613 } 4614 *ksp = snes->ksp; 4615 PetscFunctionReturn(0); 4616 } 4617 4618 4619 #include <petsc-private/dmimpl.h> 4620 #undef __FUNCT__ 4621 #define __FUNCT__ "SNESSetDM" 4622 /*@ 4623 SNESSetDM - Sets the DM that may be used by some preconditioners 4624 4625 Logically Collective on SNES 4626 4627 Input Parameters: 4628 + snes - the preconditioner context 4629 - dm - the dm 4630 4631 Level: intermediate 4632 4633 .seealso: SNESGetDM(), KSPSetDM(), KSPGetDM() 4634 @*/ 4635 PetscErrorCode SNESSetDM(SNES snes,DM dm) 4636 { 4637 PetscErrorCode ierr; 4638 KSP ksp; 4639 DMSNES sdm; 4640 4641 PetscFunctionBegin; 4642 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4643 if (dm) {ierr = PetscObjectReference((PetscObject)dm);CHKERRQ(ierr);} 4644 if (snes->dm) { /* Move the DMSNES context over to the new DM unless the new DM already has one */ 4645 if (snes->dm->dmsnes && snes->dmAuto && !dm->dmsnes) { 4646 ierr = DMCopyDMSNES(snes->dm,dm);CHKERRQ(ierr); 4647 ierr = DMGetDMSNES(snes->dm,&sdm);CHKERRQ(ierr); 4648 if (sdm->originaldm == snes->dm) sdm->originaldm = dm; /* Grant write privileges to the replacement DM */ 4649 } 4650 ierr = DMDestroy(&snes->dm);CHKERRQ(ierr); 4651 } 4652 snes->dm = dm; 4653 snes->dmAuto = PETSC_FALSE; 4654 4655 ierr = SNESGetKSP(snes,&ksp);CHKERRQ(ierr); 4656 ierr = KSPSetDM(ksp,dm);CHKERRQ(ierr); 4657 ierr = KSPSetDMActive(ksp,PETSC_FALSE);CHKERRQ(ierr); 4658 if (snes->pc) { 4659 ierr = SNESSetDM(snes->pc, snes->dm);CHKERRQ(ierr); 4660 ierr = SNESSetPCSide(snes,snes->pcside);CHKERRQ(ierr); 4661 } 4662 PetscFunctionReturn(0); 4663 } 4664 4665 #undef __FUNCT__ 4666 #define __FUNCT__ "SNESGetDM" 4667 /*@ 4668 SNESGetDM - Gets the DM that may be used by some preconditioners 4669 4670 Not Collective but DM obtained is parallel on SNES 4671 4672 Input Parameter: 4673 . snes - the preconditioner context 4674 4675 Output Parameter: 4676 . dm - the dm 4677 4678 Level: intermediate 4679 4680 .seealso: SNESSetDM(), KSPSetDM(), KSPGetDM() 4681 @*/ 4682 PetscErrorCode SNESGetDM(SNES snes,DM *dm) 4683 { 4684 PetscErrorCode ierr; 4685 4686 PetscFunctionBegin; 4687 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4688 if (!snes->dm) { 4689 ierr = DMShellCreate(PetscObjectComm((PetscObject)snes),&snes->dm);CHKERRQ(ierr); 4690 snes->dmAuto = PETSC_TRUE; 4691 } 4692 *dm = snes->dm; 4693 PetscFunctionReturn(0); 4694 } 4695 4696 #undef __FUNCT__ 4697 #define __FUNCT__ "SNESSetPC" 4698 /*@ 4699 SNESSetPC - Sets the nonlinear preconditioner to be used. 4700 4701 Collective on SNES 4702 4703 Input Parameters: 4704 + snes - iterative context obtained from SNESCreate() 4705 - pc - the preconditioner object 4706 4707 Notes: 4708 Use SNESGetPC() to retrieve the preconditioner context (for example, 4709 to configure it using the API). 4710 4711 Level: developer 4712 4713 .keywords: SNES, set, precondition 4714 .seealso: SNESGetPC() 4715 @*/ 4716 PetscErrorCode SNESSetPC(SNES snes, SNES pc) 4717 { 4718 PetscErrorCode ierr; 4719 4720 PetscFunctionBegin; 4721 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); 4722 PetscValidHeaderSpecific(pc, SNES_CLASSID, 2); 4723 PetscCheckSameComm(snes, 1, pc, 2); 4724 ierr = PetscObjectReference((PetscObject) pc);CHKERRQ(ierr); 4725 ierr = SNESDestroy(&snes->pc);CHKERRQ(ierr); 4726 snes->pc = pc; 4727 ierr = PetscLogObjectParent((PetscObject)snes, (PetscObject)snes->pc);CHKERRQ(ierr); 4728 PetscFunctionReturn(0); 4729 } 4730 4731 #undef __FUNCT__ 4732 #define __FUNCT__ "SNESGetPC" 4733 /*@ 4734 SNESGetPC - Returns a pointer to the nonlinear preconditioning context set with SNESSetPC(). 4735 4736 Not Collective 4737 4738 Input Parameter: 4739 . snes - iterative context obtained from SNESCreate() 4740 4741 Output Parameter: 4742 . pc - preconditioner context 4743 4744 Level: developer 4745 4746 .keywords: SNES, get, preconditioner 4747 .seealso: SNESSetPC() 4748 @*/ 4749 PetscErrorCode SNESGetPC(SNES snes, SNES *pc) 4750 { 4751 PetscErrorCode ierr; 4752 const char *optionsprefix; 4753 4754 PetscFunctionBegin; 4755 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); 4756 PetscValidPointer(pc, 2); 4757 if (!snes->pc) { 4758 ierr = SNESCreate(PetscObjectComm((PetscObject)snes),&snes->pc);CHKERRQ(ierr); 4759 ierr = PetscObjectIncrementTabLevel((PetscObject)snes->pc,(PetscObject)snes,1);CHKERRQ(ierr); 4760 ierr = PetscLogObjectParent((PetscObject)snes,(PetscObject)snes->pc);CHKERRQ(ierr); 4761 ierr = SNESGetOptionsPrefix(snes,&optionsprefix);CHKERRQ(ierr); 4762 ierr = SNESSetOptionsPrefix(snes->pc,optionsprefix);CHKERRQ(ierr); 4763 ierr = SNESAppendOptionsPrefix(snes->pc,"npc_");CHKERRQ(ierr); 4764 ierr = SNESSetCountersReset(snes->pc,PETSC_FALSE);CHKERRQ(ierr); 4765 } 4766 *pc = snes->pc; 4767 PetscFunctionReturn(0); 4768 } 4769 4770 #undef __FUNCT__ 4771 #define __FUNCT__ "SNESSetPCSide" 4772 /*@ 4773 SNESSetPCSide - Sets the preconditioning side. 4774 4775 Logically Collective on SNES 4776 4777 Input Parameter: 4778 . snes - iterative context obtained from SNESCreate() 4779 4780 Output Parameter: 4781 . side - the preconditioning side, where side is one of 4782 .vb 4783 PC_LEFT - left preconditioning (default) 4784 PC_RIGHT - right preconditioning 4785 .ve 4786 4787 Options Database Keys: 4788 . -snes_pc_side <right,left> 4789 4790 Level: intermediate 4791 4792 .keywords: SNES, set, right, left, side, preconditioner, flag 4793 4794 .seealso: SNESGetPCSide(), KSPSetPCSide() 4795 @*/ 4796 PetscErrorCode SNESSetPCSide(SNES snes,PCSide side) 4797 { 4798 PetscFunctionBegin; 4799 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4800 PetscValidLogicalCollectiveEnum(snes,side,2); 4801 snes->pcside = side; 4802 PetscFunctionReturn(0); 4803 } 4804 4805 #undef __FUNCT__ 4806 #define __FUNCT__ "SNESGetPCSide" 4807 /*@ 4808 SNESGetPCSide - Gets the preconditioning side. 4809 4810 Not Collective 4811 4812 Input Parameter: 4813 . snes - iterative context obtained from SNESCreate() 4814 4815 Output Parameter: 4816 . side - the preconditioning side, where side is one of 4817 .vb 4818 PC_LEFT - left preconditioning (default) 4819 PC_RIGHT - right preconditioning 4820 .ve 4821 4822 Level: intermediate 4823 4824 .keywords: SNES, get, right, left, side, preconditioner, flag 4825 4826 .seealso: SNESSetPCSide(), KSPGetPCSide() 4827 @*/ 4828 PetscErrorCode SNESGetPCSide(SNES snes,PCSide *side) 4829 { 4830 PetscFunctionBegin; 4831 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4832 PetscValidPointer(side,2); 4833 *side = snes->pcside; 4834 PetscFunctionReturn(0); 4835 } 4836 4837 #undef __FUNCT__ 4838 #define __FUNCT__ "SNESSetLineSearch" 4839 /*@ 4840 SNESSetLineSearch - Sets the linesearch on the SNES instance. 4841 4842 Collective on SNES 4843 4844 Input Parameters: 4845 + snes - iterative context obtained from SNESCreate() 4846 - linesearch - the linesearch object 4847 4848 Notes: 4849 Use SNESGetLineSearch() to retrieve the preconditioner context (for example, 4850 to configure it using the API). 4851 4852 Level: developer 4853 4854 .keywords: SNES, set, linesearch 4855 .seealso: SNESGetLineSearch() 4856 @*/ 4857 PetscErrorCode SNESSetLineSearch(SNES snes, SNESLineSearch linesearch) 4858 { 4859 PetscErrorCode ierr; 4860 4861 PetscFunctionBegin; 4862 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); 4863 PetscValidHeaderSpecific(linesearch, SNESLINESEARCH_CLASSID, 2); 4864 PetscCheckSameComm(snes, 1, linesearch, 2); 4865 ierr = PetscObjectReference((PetscObject) linesearch);CHKERRQ(ierr); 4866 ierr = SNESLineSearchDestroy(&snes->linesearch);CHKERRQ(ierr); 4867 4868 snes->linesearch = linesearch; 4869 4870 ierr = PetscLogObjectParent((PetscObject)snes, (PetscObject)snes->linesearch);CHKERRQ(ierr); 4871 PetscFunctionReturn(0); 4872 } 4873 4874 #undef __FUNCT__ 4875 #define __FUNCT__ "SNESGetLineSearch" 4876 /*@ 4877 SNESGetLineSearch - Returns a pointer to the line search context set with SNESSetLineSearch() 4878 or creates a default line search instance associated with the SNES and returns it. 4879 4880 Not Collective 4881 4882 Input Parameter: 4883 . snes - iterative context obtained from SNESCreate() 4884 4885 Output Parameter: 4886 . linesearch - linesearch context 4887 4888 Level: developer 4889 4890 .keywords: SNES, get, linesearch 4891 .seealso: SNESSetLineSearch() 4892 @*/ 4893 PetscErrorCode SNESGetLineSearch(SNES snes, SNESLineSearch *linesearch) 4894 { 4895 PetscErrorCode ierr; 4896 const char *optionsprefix; 4897 4898 PetscFunctionBegin; 4899 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); 4900 PetscValidPointer(linesearch, 2); 4901 if (!snes->linesearch) { 4902 ierr = SNESGetOptionsPrefix(snes, &optionsprefix);CHKERRQ(ierr); 4903 ierr = SNESLineSearchCreate(PetscObjectComm((PetscObject)snes), &snes->linesearch);CHKERRQ(ierr); 4904 ierr = SNESLineSearchSetSNES(snes->linesearch, snes);CHKERRQ(ierr); 4905 ierr = SNESLineSearchAppendOptionsPrefix(snes->linesearch, optionsprefix);CHKERRQ(ierr); 4906 ierr = PetscObjectIncrementTabLevel((PetscObject) snes->linesearch, (PetscObject) snes, 1);CHKERRQ(ierr); 4907 ierr = PetscLogObjectParent((PetscObject)snes, (PetscObject)snes->linesearch);CHKERRQ(ierr); 4908 } 4909 *linesearch = snes->linesearch; 4910 PetscFunctionReturn(0); 4911 } 4912 4913 #if defined(PETSC_HAVE_MATLAB_ENGINE) 4914 #include <mex.h> 4915 4916 typedef struct {char *funcname; mxArray *ctx;} SNESMatlabContext; 4917 4918 #undef __FUNCT__ 4919 #define __FUNCT__ "SNESComputeFunction_Matlab" 4920 /* 4921 SNESComputeFunction_Matlab - Calls the function that has been set with SNESSetFunctionMatlab(). 4922 4923 Collective on SNES 4924 4925 Input Parameters: 4926 + snes - the SNES context 4927 - x - input vector 4928 4929 Output Parameter: 4930 . y - function vector, as set by SNESSetFunction() 4931 4932 Notes: 4933 SNESComputeFunction() is typically used within nonlinear solvers 4934 implementations, so most users would not generally call this routine 4935 themselves. 4936 4937 Level: developer 4938 4939 .keywords: SNES, nonlinear, compute, function 4940 4941 .seealso: SNESSetFunction(), SNESGetFunction() 4942 */ 4943 PetscErrorCode SNESComputeFunction_Matlab(SNES snes,Vec x,Vec y, void *ctx) 4944 { 4945 PetscErrorCode ierr; 4946 SNESMatlabContext *sctx = (SNESMatlabContext*)ctx; 4947 int nlhs = 1,nrhs = 5; 4948 mxArray *plhs[1],*prhs[5]; 4949 long long int lx = 0,ly = 0,ls = 0; 4950 4951 PetscFunctionBegin; 4952 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4953 PetscValidHeaderSpecific(x,VEC_CLASSID,2); 4954 PetscValidHeaderSpecific(y,VEC_CLASSID,3); 4955 PetscCheckSameComm(snes,1,x,2); 4956 PetscCheckSameComm(snes,1,y,3); 4957 4958 /* call Matlab function in ctx with arguments x and y */ 4959 4960 ierr = PetscMemcpy(&ls,&snes,sizeof(snes));CHKERRQ(ierr); 4961 ierr = PetscMemcpy(&lx,&x,sizeof(x));CHKERRQ(ierr); 4962 ierr = PetscMemcpy(&ly,&y,sizeof(x));CHKERRQ(ierr); 4963 prhs[0] = mxCreateDoubleScalar((double)ls); 4964 prhs[1] = mxCreateDoubleScalar((double)lx); 4965 prhs[2] = mxCreateDoubleScalar((double)ly); 4966 prhs[3] = mxCreateString(sctx->funcname); 4967 prhs[4] = sctx->ctx; 4968 ierr = mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscSNESComputeFunctionInternal");CHKERRQ(ierr); 4969 ierr = mxGetScalar(plhs[0]);CHKERRQ(ierr); 4970 mxDestroyArray(prhs[0]); 4971 mxDestroyArray(prhs[1]); 4972 mxDestroyArray(prhs[2]); 4973 mxDestroyArray(prhs[3]); 4974 mxDestroyArray(plhs[0]); 4975 PetscFunctionReturn(0); 4976 } 4977 4978 #undef __FUNCT__ 4979 #define __FUNCT__ "SNESSetFunctionMatlab" 4980 /* 4981 SNESSetFunctionMatlab - Sets the function evaluation routine and function 4982 vector for use by the SNES routines in solving systems of nonlinear 4983 equations from MATLAB. Here the function is a string containing the name of a MATLAB function 4984 4985 Logically Collective on SNES 4986 4987 Input Parameters: 4988 + snes - the SNES context 4989 . r - vector to store function value 4990 - SNESFunction - function evaluation routine 4991 4992 Notes: 4993 The Newton-like methods typically solve linear systems of the form 4994 $ f'(x) x = -f(x), 4995 where f'(x) denotes the Jacobian matrix and f(x) is the function. 4996 4997 Level: beginner 4998 4999 Developer Note: This bleeds the allocated memory SNESMatlabContext *sctx; 5000 5001 .keywords: SNES, nonlinear, set, function 5002 5003 .seealso: SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESSetFunction() 5004 */ 5005 PetscErrorCode SNESSetFunctionMatlab(SNES snes,Vec r,const char *SNESFunction,mxArray *ctx) 5006 { 5007 PetscErrorCode ierr; 5008 SNESMatlabContext *sctx; 5009 5010 PetscFunctionBegin; 5011 /* currently sctx is memory bleed */ 5012 ierr = PetscMalloc(sizeof(SNESMatlabContext),&sctx);CHKERRQ(ierr); 5013 ierr = PetscStrallocpy(SNESFunction,&sctx->funcname);CHKERRQ(ierr); 5014 /* 5015 This should work, but it doesn't 5016 sctx->ctx = ctx; 5017 mexMakeArrayPersistent(sctx->ctx); 5018 */ 5019 sctx->ctx = mxDuplicateArray(ctx); 5020 ierr = SNESSetFunction(snes,r,SNESComputeFunction_Matlab,sctx);CHKERRQ(ierr); 5021 PetscFunctionReturn(0); 5022 } 5023 5024 #undef __FUNCT__ 5025 #define __FUNCT__ "SNESComputeJacobian_Matlab" 5026 /* 5027 SNESComputeJacobian_Matlab - Calls the function that has been set with SNESSetJacobianMatlab(). 5028 5029 Collective on SNES 5030 5031 Input Parameters: 5032 + snes - the SNES context 5033 . x - input vector 5034 . A, B - the matrices 5035 - ctx - user context 5036 5037 Output Parameter: 5038 . flag - structure of the matrix 5039 5040 Level: developer 5041 5042 .keywords: SNES, nonlinear, compute, function 5043 5044 .seealso: SNESSetFunction(), SNESGetFunction() 5045 @*/ 5046 PetscErrorCode SNESComputeJacobian_Matlab(SNES snes,Vec x,Mat *A,Mat *B,MatStructure *flag, void *ctx) 5047 { 5048 PetscErrorCode ierr; 5049 SNESMatlabContext *sctx = (SNESMatlabContext*)ctx; 5050 int nlhs = 2,nrhs = 6; 5051 mxArray *plhs[2],*prhs[6]; 5052 long long int lx = 0,lA = 0,ls = 0, lB = 0; 5053 5054 PetscFunctionBegin; 5055 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 5056 PetscValidHeaderSpecific(x,VEC_CLASSID,2); 5057 5058 /* call Matlab function in ctx with arguments x and y */ 5059 5060 ierr = PetscMemcpy(&ls,&snes,sizeof(snes));CHKERRQ(ierr); 5061 ierr = PetscMemcpy(&lx,&x,sizeof(x));CHKERRQ(ierr); 5062 ierr = PetscMemcpy(&lA,A,sizeof(x));CHKERRQ(ierr); 5063 ierr = PetscMemcpy(&lB,B,sizeof(x));CHKERRQ(ierr); 5064 prhs[0] = mxCreateDoubleScalar((double)ls); 5065 prhs[1] = mxCreateDoubleScalar((double)lx); 5066 prhs[2] = mxCreateDoubleScalar((double)lA); 5067 prhs[3] = mxCreateDoubleScalar((double)lB); 5068 prhs[4] = mxCreateString(sctx->funcname); 5069 prhs[5] = sctx->ctx; 5070 ierr = mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscSNESComputeJacobianInternal");CHKERRQ(ierr); 5071 ierr = mxGetScalar(plhs[0]);CHKERRQ(ierr); 5072 *flag = (MatStructure) mxGetScalar(plhs[1]);CHKERRQ(ierr); 5073 mxDestroyArray(prhs[0]); 5074 mxDestroyArray(prhs[1]); 5075 mxDestroyArray(prhs[2]); 5076 mxDestroyArray(prhs[3]); 5077 mxDestroyArray(prhs[4]); 5078 mxDestroyArray(plhs[0]); 5079 mxDestroyArray(plhs[1]); 5080 PetscFunctionReturn(0); 5081 } 5082 5083 #undef __FUNCT__ 5084 #define __FUNCT__ "SNESSetJacobianMatlab" 5085 /* 5086 SNESSetJacobianMatlab - Sets the Jacobian function evaluation routine and two empty Jacobian matrices 5087 vector for use by the SNES routines in solving systems of nonlinear 5088 equations from MATLAB. Here the function is a string containing the name of a MATLAB function 5089 5090 Logically Collective on SNES 5091 5092 Input Parameters: 5093 + snes - the SNES context 5094 . A,B - Jacobian matrices 5095 . SNESJacobianFunction - function evaluation routine 5096 - ctx - user context 5097 5098 Level: developer 5099 5100 Developer Note: This bleeds the allocated memory SNESMatlabContext *sctx; 5101 5102 .keywords: SNES, nonlinear, set, function 5103 5104 .seealso: SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESSetFunction(), SNESJacobianFunction 5105 */ 5106 PetscErrorCode SNESSetJacobianMatlab(SNES snes,Mat A,Mat B,const char *SNESJacobianFunction,mxArray *ctx) 5107 { 5108 PetscErrorCode ierr; 5109 SNESMatlabContext *sctx; 5110 5111 PetscFunctionBegin; 5112 /* currently sctx is memory bleed */ 5113 ierr = PetscMalloc(sizeof(SNESMatlabContext),&sctx);CHKERRQ(ierr); 5114 ierr = PetscStrallocpy(SNESJacobianFunction,&sctx->funcname);CHKERRQ(ierr); 5115 /* 5116 This should work, but it doesn't 5117 sctx->ctx = ctx; 5118 mexMakeArrayPersistent(sctx->ctx); 5119 */ 5120 sctx->ctx = mxDuplicateArray(ctx); 5121 ierr = SNESSetJacobian(snes,A,B,SNESComputeJacobian_Matlab,sctx);CHKERRQ(ierr); 5122 PetscFunctionReturn(0); 5123 } 5124 5125 #undef __FUNCT__ 5126 #define __FUNCT__ "SNESMonitor_Matlab" 5127 /* 5128 SNESMonitor_Matlab - Calls the function that has been set with SNESMonitorSetMatlab(). 5129 5130 Collective on SNES 5131 5132 .seealso: SNESSetFunction(), SNESGetFunction() 5133 @*/ 5134 PetscErrorCode SNESMonitor_Matlab(SNES snes,PetscInt it, PetscReal fnorm, void *ctx) 5135 { 5136 PetscErrorCode ierr; 5137 SNESMatlabContext *sctx = (SNESMatlabContext*)ctx; 5138 int nlhs = 1,nrhs = 6; 5139 mxArray *plhs[1],*prhs[6]; 5140 long long int lx = 0,ls = 0; 5141 Vec x = snes->vec_sol; 5142 5143 PetscFunctionBegin; 5144 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 5145 5146 ierr = PetscMemcpy(&ls,&snes,sizeof(snes));CHKERRQ(ierr); 5147 ierr = PetscMemcpy(&lx,&x,sizeof(x));CHKERRQ(ierr); 5148 prhs[0] = mxCreateDoubleScalar((double)ls); 5149 prhs[1] = mxCreateDoubleScalar((double)it); 5150 prhs[2] = mxCreateDoubleScalar((double)fnorm); 5151 prhs[3] = mxCreateDoubleScalar((double)lx); 5152 prhs[4] = mxCreateString(sctx->funcname); 5153 prhs[5] = sctx->ctx; 5154 ierr = mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscSNESMonitorInternal");CHKERRQ(ierr); 5155 ierr = mxGetScalar(plhs[0]);CHKERRQ(ierr); 5156 mxDestroyArray(prhs[0]); 5157 mxDestroyArray(prhs[1]); 5158 mxDestroyArray(prhs[2]); 5159 mxDestroyArray(prhs[3]); 5160 mxDestroyArray(prhs[4]); 5161 mxDestroyArray(plhs[0]); 5162 PetscFunctionReturn(0); 5163 } 5164 5165 #undef __FUNCT__ 5166 #define __FUNCT__ "SNESMonitorSetMatlab" 5167 /* 5168 SNESMonitorSetMatlab - Sets the monitor function from MATLAB 5169 5170 Level: developer 5171 5172 Developer Note: This bleeds the allocated memory SNESMatlabContext *sctx; 5173 5174 .keywords: SNES, nonlinear, set, function 5175 5176 .seealso: SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESSetFunction() 5177 */ 5178 PetscErrorCode SNESMonitorSetMatlab(SNES snes,const char *SNESMonitorFunction,mxArray *ctx) 5179 { 5180 PetscErrorCode ierr; 5181 SNESMatlabContext *sctx; 5182 5183 PetscFunctionBegin; 5184 /* currently sctx is memory bleed */ 5185 ierr = PetscMalloc(sizeof(SNESMatlabContext),&sctx);CHKERRQ(ierr); 5186 ierr = PetscStrallocpy(SNESMonitorFunction,&sctx->funcname);CHKERRQ(ierr); 5187 /* 5188 This should work, but it doesn't 5189 sctx->ctx = ctx; 5190 mexMakeArrayPersistent(sctx->ctx); 5191 */ 5192 sctx->ctx = mxDuplicateArray(ctx); 5193 ierr = SNESMonitorSet(snes,SNESMonitor_Matlab,sctx,NULL);CHKERRQ(ierr); 5194 PetscFunctionReturn(0); 5195 } 5196 5197 #endif 5198