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,&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 = PetscMalloc1(na,&a);CHKERRQ(ierr); 3563 ierr = PetscMalloc1(na,&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 PetscInt grid; 3750 Vec xcreated = NULL; 3751 DM dm; 3752 3753 PetscFunctionBegin; 3754 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3755 if (x) PetscValidHeaderSpecific(x,VEC_CLASSID,3); 3756 if (x) PetscCheckSameComm(snes,1,x,3); 3757 if (b) PetscValidHeaderSpecific(b,VEC_CLASSID,2); 3758 if (b) PetscCheckSameComm(snes,1,b,2); 3759 3760 if (!x) { 3761 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 3762 ierr = DMCreateGlobalVector(dm,&xcreated);CHKERRQ(ierr); 3763 x = xcreated; 3764 } 3765 ierr = SNESViewFromOptions(snes,NULL,"-snes_view_pre");CHKERRQ(ierr); 3766 3767 for (grid=0; grid<snes->gridsequence; grid++) {ierr = PetscViewerASCIIPushTab(PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)snes)));CHKERRQ(ierr);} 3768 for (grid=0; grid<snes->gridsequence+1; grid++) { 3769 3770 /* set solution vector */ 3771 if (!grid) {ierr = PetscObjectReference((PetscObject)x);CHKERRQ(ierr);} 3772 ierr = VecDestroy(&snes->vec_sol);CHKERRQ(ierr); 3773 snes->vec_sol = x; 3774 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 3775 3776 /* set affine vector if provided */ 3777 if (b) { ierr = PetscObjectReference((PetscObject)b);CHKERRQ(ierr); } 3778 ierr = VecDestroy(&snes->vec_rhs);CHKERRQ(ierr); 3779 snes->vec_rhs = b; 3780 3781 if (snes->vec_func == snes->vec_sol) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_IDN,"Solution vector cannot be function vector"); 3782 if (snes->vec_rhs == snes->vec_sol) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_IDN,"Solution vector cannot be right hand side vector"); 3783 if (!snes->vec_sol_update /* && snes->vec_sol */) { 3784 ierr = VecDuplicate(snes->vec_sol,&snes->vec_sol_update);CHKERRQ(ierr); 3785 ierr = PetscLogObjectParent((PetscObject)snes,(PetscObject)snes->vec_sol_update);CHKERRQ(ierr); 3786 } 3787 ierr = DMShellSetGlobalVector(dm,snes->vec_sol);CHKERRQ(ierr); 3788 ierr = SNESSetUp(snes);CHKERRQ(ierr); 3789 3790 if (!grid) { 3791 if (snes->ops->computeinitialguess) { 3792 ierr = (*snes->ops->computeinitialguess)(snes,snes->vec_sol,snes->initialguessP);CHKERRQ(ierr); 3793 } 3794 } 3795 3796 if (snes->conv_hist_reset) snes->conv_hist_len = 0; 3797 if (snes->counters_reset) {snes->nfuncs = 0; snes->linear_its = 0; snes->numFailures = 0;} 3798 3799 ierr = PetscLogEventBegin(SNES_Solve,snes,0,0,0);CHKERRQ(ierr); 3800 ierr = (*snes->ops->solve)(snes);CHKERRQ(ierr); 3801 ierr = PetscLogEventEnd(SNES_Solve,snes,0,0,0);CHKERRQ(ierr); 3802 if (snes->domainerror) { 3803 snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; 3804 snes->domainerror = PETSC_FALSE; 3805 } 3806 if (!snes->reason) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Internal error, solver returned without setting converged reason"); 3807 3808 if (snes->lagjac_persist) snes->jac_iter += snes->iter; 3809 if (snes->lagpre_persist) snes->pre_iter += snes->iter; 3810 3811 flg = PETSC_FALSE; 3812 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_test_local_min",&flg,NULL);CHKERRQ(ierr); 3813 if (flg && !PetscPreLoadingOn) { ierr = SNESTestLocalMin(snes);CHKERRQ(ierr); } 3814 if (snes->printreason) { 3815 ierr = PetscViewerASCIIAddTab(PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)snes)),((PetscObject)snes)->tablevel);CHKERRQ(ierr); 3816 if (snes->reason > 0) { 3817 ierr = PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)snes)),"Nonlinear solve converged due to %s iterations %D\n",SNESConvergedReasons[snes->reason],snes->iter);CHKERRQ(ierr); 3818 } else { 3819 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); 3820 } 3821 ierr = PetscViewerASCIISubtractTab(PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)snes)),((PetscObject)snes)->tablevel);CHKERRQ(ierr); 3822 } 3823 3824 if (snes->errorifnotconverged && snes->reason < 0) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_NOT_CONVERGED,"SNESSolve has not converged"); 3825 if (grid < snes->gridsequence) { 3826 DM fine; 3827 Vec xnew; 3828 Mat interp; 3829 3830 ierr = DMRefine(snes->dm,PetscObjectComm((PetscObject)snes),&fine);CHKERRQ(ierr); 3831 if (!fine) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_INCOMP,"DMRefine() did not perform any refinement, cannot continue grid sequencing"); 3832 ierr = DMCreateInterpolation(snes->dm,fine,&interp,NULL);CHKERRQ(ierr); 3833 ierr = DMCreateGlobalVector(fine,&xnew);CHKERRQ(ierr); 3834 ierr = MatInterpolate(interp,x,xnew);CHKERRQ(ierr); 3835 ierr = DMInterpolate(snes->dm,interp,fine);CHKERRQ(ierr); 3836 ierr = MatDestroy(&interp);CHKERRQ(ierr); 3837 x = xnew; 3838 3839 ierr = SNESReset(snes);CHKERRQ(ierr); 3840 ierr = SNESSetDM(snes,fine);CHKERRQ(ierr); 3841 ierr = DMDestroy(&fine);CHKERRQ(ierr); 3842 ierr = PetscViewerASCIIPopTab(PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)snes)));CHKERRQ(ierr); 3843 } 3844 } 3845 ierr = SNESViewFromOptions(snes,NULL,"-snes_view");CHKERRQ(ierr); 3846 ierr = VecViewFromOptions(snes->vec_sol,((PetscObject)snes)->prefix,"-snes_view_solution");CHKERRQ(ierr); 3847 3848 ierr = VecDestroy(&xcreated);CHKERRQ(ierr); 3849 ierr = PetscObjectSAWsBlock((PetscObject)snes);CHKERRQ(ierr); 3850 PetscFunctionReturn(0); 3851 } 3852 3853 /* --------- Internal routines for SNES Package --------- */ 3854 3855 #undef __FUNCT__ 3856 #define __FUNCT__ "SNESSetType" 3857 /*@C 3858 SNESSetType - Sets the method for the nonlinear solver. 3859 3860 Collective on SNES 3861 3862 Input Parameters: 3863 + snes - the SNES context 3864 - type - a known method 3865 3866 Options Database Key: 3867 . -snes_type <type> - Sets the method; use -help for a list 3868 of available methods (for instance, newtonls or newtontr) 3869 3870 Notes: 3871 See "petsc/include/petscsnes.h" for available methods (for instance) 3872 + SNESNEWTONLS - Newton's method with line search 3873 (systems of nonlinear equations) 3874 . SNESNEWTONTR - Newton's method with trust region 3875 (systems of nonlinear equations) 3876 3877 Normally, it is best to use the SNESSetFromOptions() command and then 3878 set the SNES solver type from the options database rather than by using 3879 this routine. Using the options database provides the user with 3880 maximum flexibility in evaluating the many nonlinear solvers. 3881 The SNESSetType() routine is provided for those situations where it 3882 is necessary to set the nonlinear solver independently of the command 3883 line or options database. This might be the case, for example, when 3884 the choice of solver changes during the execution of the program, 3885 and the user's application is taking responsibility for choosing the 3886 appropriate method. 3887 3888 Developer Notes: SNESRegister() adds a constructor for a new SNESType to SNESList, SNESSetType() locates 3889 the constructor in that list and calls it to create the spexific object. 3890 3891 Level: intermediate 3892 3893 .keywords: SNES, set, type 3894 3895 .seealso: SNESType, SNESCreate(), SNESDestroy(), SNESGetType(), SNESSetFromOptions() 3896 3897 @*/ 3898 PetscErrorCode SNESSetType(SNES snes,SNESType type) 3899 { 3900 PetscErrorCode ierr,(*r)(SNES); 3901 PetscBool match; 3902 3903 PetscFunctionBegin; 3904 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3905 PetscValidCharPointer(type,2); 3906 3907 ierr = PetscObjectTypeCompare((PetscObject)snes,type,&match);CHKERRQ(ierr); 3908 if (match) PetscFunctionReturn(0); 3909 3910 ierr = PetscFunctionListFind(SNESList,type,&r);CHKERRQ(ierr); 3911 if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unable to find requested SNES type %s",type); 3912 /* Destroy the previous private SNES context */ 3913 if (snes->ops->destroy) { 3914 ierr = (*(snes)->ops->destroy)(snes);CHKERRQ(ierr); 3915 snes->ops->destroy = NULL; 3916 } 3917 /* Reinitialize function pointers in SNESOps structure */ 3918 snes->ops->setup = 0; 3919 snes->ops->solve = 0; 3920 snes->ops->view = 0; 3921 snes->ops->setfromoptions = 0; 3922 snes->ops->destroy = 0; 3923 /* Call the SNESCreate_XXX routine for this particular Nonlinear solver */ 3924 snes->setupcalled = PETSC_FALSE; 3925 3926 ierr = PetscObjectChangeTypeName((PetscObject)snes,type);CHKERRQ(ierr); 3927 ierr = (*r)(snes);CHKERRQ(ierr); 3928 PetscFunctionReturn(0); 3929 } 3930 3931 #undef __FUNCT__ 3932 #define __FUNCT__ "SNESGetType" 3933 /*@C 3934 SNESGetType - Gets the SNES method type and name (as a string). 3935 3936 Not Collective 3937 3938 Input Parameter: 3939 . snes - nonlinear solver context 3940 3941 Output Parameter: 3942 . type - SNES method (a character string) 3943 3944 Level: intermediate 3945 3946 .keywords: SNES, nonlinear, get, type, name 3947 @*/ 3948 PetscErrorCode SNESGetType(SNES snes,SNESType *type) 3949 { 3950 PetscFunctionBegin; 3951 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3952 PetscValidPointer(type,2); 3953 *type = ((PetscObject)snes)->type_name; 3954 PetscFunctionReturn(0); 3955 } 3956 3957 #undef __FUNCT__ 3958 #define __FUNCT__ "SNESGetSolution" 3959 /*@ 3960 SNESGetSolution - Returns the vector where the approximate solution is 3961 stored. This is the fine grid solution when using SNESSetGridSequence(). 3962 3963 Not Collective, but Vec is parallel if SNES is parallel 3964 3965 Input Parameter: 3966 . snes - the SNES context 3967 3968 Output Parameter: 3969 . x - the solution 3970 3971 Level: intermediate 3972 3973 .keywords: SNES, nonlinear, get, solution 3974 3975 .seealso: SNESGetSolutionUpdate(), SNESGetFunction() 3976 @*/ 3977 PetscErrorCode SNESGetSolution(SNES snes,Vec *x) 3978 { 3979 PetscFunctionBegin; 3980 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3981 PetscValidPointer(x,2); 3982 *x = snes->vec_sol; 3983 PetscFunctionReturn(0); 3984 } 3985 3986 #undef __FUNCT__ 3987 #define __FUNCT__ "SNESGetSolutionUpdate" 3988 /*@ 3989 SNESGetSolutionUpdate - Returns the vector where the solution update is 3990 stored. 3991 3992 Not Collective, but Vec is parallel if SNES is parallel 3993 3994 Input Parameter: 3995 . snes - the SNES context 3996 3997 Output Parameter: 3998 . x - the solution update 3999 4000 Level: advanced 4001 4002 .keywords: SNES, nonlinear, get, solution, update 4003 4004 .seealso: SNESGetSolution(), SNESGetFunction() 4005 @*/ 4006 PetscErrorCode SNESGetSolutionUpdate(SNES snes,Vec *x) 4007 { 4008 PetscFunctionBegin; 4009 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4010 PetscValidPointer(x,2); 4011 *x = snes->vec_sol_update; 4012 PetscFunctionReturn(0); 4013 } 4014 4015 #undef __FUNCT__ 4016 #define __FUNCT__ "SNESGetFunction" 4017 /*@C 4018 SNESGetFunction - Returns the vector where the function is stored. 4019 4020 Not Collective, but Vec is parallel if SNES is parallel. Collective if Vec is requested, but has not been created yet. 4021 4022 Input Parameter: 4023 . snes - the SNES context 4024 4025 Output Parameter: 4026 + r - the vector that is used to store residuals (or NULL if you don't want it) 4027 . SNESFunction- the function (or NULL if you don't want it) 4028 - ctx - the function context (or NULL if you don't want it) 4029 4030 Level: advanced 4031 4032 .keywords: SNES, nonlinear, get, function 4033 4034 .seealso: SNESSetFunction(), SNESGetSolution(), SNESFunction 4035 @*/ 4036 PetscErrorCode SNESGetFunction(SNES snes,Vec *r,PetscErrorCode (**SNESFunction)(SNES,Vec,Vec,void*),void **ctx) 4037 { 4038 PetscErrorCode ierr; 4039 DM dm; 4040 4041 PetscFunctionBegin; 4042 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4043 if (r) { 4044 if (!snes->vec_func) { 4045 if (snes->vec_rhs) { 4046 ierr = VecDuplicate(snes->vec_rhs,&snes->vec_func);CHKERRQ(ierr); 4047 } else if (snes->vec_sol) { 4048 ierr = VecDuplicate(snes->vec_sol,&snes->vec_func);CHKERRQ(ierr); 4049 } else if (snes->dm) { 4050 ierr = DMCreateGlobalVector(snes->dm,&snes->vec_func);CHKERRQ(ierr); 4051 } 4052 } 4053 *r = snes->vec_func; 4054 } 4055 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 4056 ierr = DMSNESGetFunction(dm,SNESFunction,ctx);CHKERRQ(ierr); 4057 PetscFunctionReturn(0); 4058 } 4059 4060 /*@C 4061 SNESGetGS - Returns the GS function and context. 4062 4063 Input Parameter: 4064 . snes - the SNES context 4065 4066 Output Parameter: 4067 + SNESGSFunction - the function (or NULL) 4068 - ctx - the function context (or NULL) 4069 4070 Level: advanced 4071 4072 .keywords: SNES, nonlinear, get, function 4073 4074 .seealso: SNESSetGS(), SNESGetFunction() 4075 @*/ 4076 4077 #undef __FUNCT__ 4078 #define __FUNCT__ "SNESGetGS" 4079 PetscErrorCode SNESGetGS (SNES snes, PetscErrorCode (**SNESGSFunction)(SNES, Vec, Vec, void*), void ** ctx) 4080 { 4081 PetscErrorCode ierr; 4082 DM dm; 4083 4084 PetscFunctionBegin; 4085 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4086 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 4087 ierr = DMSNESGetGS(dm,SNESGSFunction,ctx);CHKERRQ(ierr); 4088 PetscFunctionReturn(0); 4089 } 4090 4091 #undef __FUNCT__ 4092 #define __FUNCT__ "SNESSetOptionsPrefix" 4093 /*@C 4094 SNESSetOptionsPrefix - Sets the prefix used for searching for all 4095 SNES options in the database. 4096 4097 Logically Collective on SNES 4098 4099 Input Parameter: 4100 + snes - the SNES context 4101 - prefix - the prefix to prepend to all option names 4102 4103 Notes: 4104 A hyphen (-) must NOT be given at the beginning of the prefix name. 4105 The first character of all runtime options is AUTOMATICALLY the hyphen. 4106 4107 Level: advanced 4108 4109 .keywords: SNES, set, options, prefix, database 4110 4111 .seealso: SNESSetFromOptions() 4112 @*/ 4113 PetscErrorCode SNESSetOptionsPrefix(SNES snes,const char prefix[]) 4114 { 4115 PetscErrorCode ierr; 4116 4117 PetscFunctionBegin; 4118 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4119 ierr = PetscObjectSetOptionsPrefix((PetscObject)snes,prefix);CHKERRQ(ierr); 4120 if (!snes->ksp) {ierr = SNESGetKSP(snes,&snes->ksp);CHKERRQ(ierr);} 4121 if (snes->linesearch) { 4122 ierr = SNESGetLineSearch(snes,&snes->linesearch);CHKERRQ(ierr); 4123 ierr = PetscObjectSetOptionsPrefix((PetscObject)snes->linesearch,prefix);CHKERRQ(ierr); 4124 } 4125 ierr = KSPSetOptionsPrefix(snes->ksp,prefix);CHKERRQ(ierr); 4126 PetscFunctionReturn(0); 4127 } 4128 4129 #undef __FUNCT__ 4130 #define __FUNCT__ "SNESAppendOptionsPrefix" 4131 /*@C 4132 SNESAppendOptionsPrefix - Appends to the prefix used for searching for all 4133 SNES options in the database. 4134 4135 Logically Collective on SNES 4136 4137 Input Parameters: 4138 + snes - the SNES context 4139 - prefix - the prefix to prepend to all option names 4140 4141 Notes: 4142 A hyphen (-) must NOT be given at the beginning of the prefix name. 4143 The first character of all runtime options is AUTOMATICALLY the hyphen. 4144 4145 Level: advanced 4146 4147 .keywords: SNES, append, options, prefix, database 4148 4149 .seealso: SNESGetOptionsPrefix() 4150 @*/ 4151 PetscErrorCode SNESAppendOptionsPrefix(SNES snes,const char prefix[]) 4152 { 4153 PetscErrorCode ierr; 4154 4155 PetscFunctionBegin; 4156 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4157 ierr = PetscObjectAppendOptionsPrefix((PetscObject)snes,prefix);CHKERRQ(ierr); 4158 if (!snes->ksp) {ierr = SNESGetKSP(snes,&snes->ksp);CHKERRQ(ierr);} 4159 if (snes->linesearch) { 4160 ierr = SNESGetLineSearch(snes,&snes->linesearch);CHKERRQ(ierr); 4161 ierr = PetscObjectAppendOptionsPrefix((PetscObject)snes->linesearch,prefix);CHKERRQ(ierr); 4162 } 4163 ierr = KSPAppendOptionsPrefix(snes->ksp,prefix);CHKERRQ(ierr); 4164 PetscFunctionReturn(0); 4165 } 4166 4167 #undef __FUNCT__ 4168 #define __FUNCT__ "SNESGetOptionsPrefix" 4169 /*@C 4170 SNESGetOptionsPrefix - Sets the prefix used for searching for all 4171 SNES options in the database. 4172 4173 Not Collective 4174 4175 Input Parameter: 4176 . snes - the SNES context 4177 4178 Output Parameter: 4179 . prefix - pointer to the prefix string used 4180 4181 Notes: On the fortran side, the user should pass in a string 'prefix' of 4182 sufficient length to hold the prefix. 4183 4184 Level: advanced 4185 4186 .keywords: SNES, get, options, prefix, database 4187 4188 .seealso: SNESAppendOptionsPrefix() 4189 @*/ 4190 PetscErrorCode SNESGetOptionsPrefix(SNES snes,const char *prefix[]) 4191 { 4192 PetscErrorCode ierr; 4193 4194 PetscFunctionBegin; 4195 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4196 ierr = PetscObjectGetOptionsPrefix((PetscObject)snes,prefix);CHKERRQ(ierr); 4197 PetscFunctionReturn(0); 4198 } 4199 4200 4201 #undef __FUNCT__ 4202 #define __FUNCT__ "SNESRegister" 4203 /*@C 4204 SNESRegister - Adds a method to the nonlinear solver package. 4205 4206 Not collective 4207 4208 Input Parameters: 4209 + name_solver - name of a new user-defined solver 4210 - routine_create - routine to create method context 4211 4212 Notes: 4213 SNESRegister() may be called multiple times to add several user-defined solvers. 4214 4215 Sample usage: 4216 .vb 4217 SNESRegister("my_solver",MySolverCreate); 4218 .ve 4219 4220 Then, your solver can be chosen with the procedural interface via 4221 $ SNESSetType(snes,"my_solver") 4222 or at runtime via the option 4223 $ -snes_type my_solver 4224 4225 Level: advanced 4226 4227 Note: If your function is not being put into a shared library then use SNESRegister() instead 4228 4229 .keywords: SNES, nonlinear, register 4230 4231 .seealso: SNESRegisterAll(), SNESRegisterDestroy() 4232 4233 Level: advanced 4234 @*/ 4235 PetscErrorCode SNESRegister(const char sname[],PetscErrorCode (*function)(SNES)) 4236 { 4237 PetscErrorCode ierr; 4238 4239 PetscFunctionBegin; 4240 ierr = PetscFunctionListAdd(&SNESList,sname,function);CHKERRQ(ierr); 4241 PetscFunctionReturn(0); 4242 } 4243 4244 #undef __FUNCT__ 4245 #define __FUNCT__ "SNESTestLocalMin" 4246 PetscErrorCode SNESTestLocalMin(SNES snes) 4247 { 4248 PetscErrorCode ierr; 4249 PetscInt N,i,j; 4250 Vec u,uh,fh; 4251 PetscScalar value; 4252 PetscReal norm; 4253 4254 PetscFunctionBegin; 4255 ierr = SNESGetSolution(snes,&u);CHKERRQ(ierr); 4256 ierr = VecDuplicate(u,&uh);CHKERRQ(ierr); 4257 ierr = VecDuplicate(u,&fh);CHKERRQ(ierr); 4258 4259 /* currently only works for sequential */ 4260 ierr = PetscPrintf(PETSC_COMM_WORLD,"Testing FormFunction() for local min\n");CHKERRQ(ierr); 4261 ierr = VecGetSize(u,&N);CHKERRQ(ierr); 4262 for (i=0; i<N; i++) { 4263 ierr = VecCopy(u,uh);CHKERRQ(ierr); 4264 ierr = PetscPrintf(PETSC_COMM_WORLD,"i = %D\n",i);CHKERRQ(ierr); 4265 for (j=-10; j<11; j++) { 4266 value = PetscSign(j)*exp(PetscAbs(j)-10.0); 4267 ierr = VecSetValue(uh,i,value,ADD_VALUES);CHKERRQ(ierr); 4268 ierr = SNESComputeFunction(snes,uh,fh);CHKERRQ(ierr); 4269 ierr = VecNorm(fh,NORM_2,&norm);CHKERRQ(ierr); 4270 ierr = PetscPrintf(PETSC_COMM_WORLD," j norm %D %18.16e\n",j,norm);CHKERRQ(ierr); 4271 value = -value; 4272 ierr = VecSetValue(uh,i,value,ADD_VALUES);CHKERRQ(ierr); 4273 } 4274 } 4275 ierr = VecDestroy(&uh);CHKERRQ(ierr); 4276 ierr = VecDestroy(&fh);CHKERRQ(ierr); 4277 PetscFunctionReturn(0); 4278 } 4279 4280 #undef __FUNCT__ 4281 #define __FUNCT__ "SNESKSPSetUseEW" 4282 /*@ 4283 SNESKSPSetUseEW - Sets SNES use Eisenstat-Walker method for 4284 computing relative tolerance for linear solvers within an inexact 4285 Newton method. 4286 4287 Logically Collective on SNES 4288 4289 Input Parameters: 4290 + snes - SNES context 4291 - flag - PETSC_TRUE or PETSC_FALSE 4292 4293 Options Database: 4294 + -snes_ksp_ew - use Eisenstat-Walker method for determining linear system convergence 4295 . -snes_ksp_ew_version ver - version of Eisenstat-Walker method 4296 . -snes_ksp_ew_rtol0 <rtol0> - Sets rtol0 4297 . -snes_ksp_ew_rtolmax <rtolmax> - Sets rtolmax 4298 . -snes_ksp_ew_gamma <gamma> - Sets gamma 4299 . -snes_ksp_ew_alpha <alpha> - Sets alpha 4300 . -snes_ksp_ew_alpha2 <alpha2> - Sets alpha2 4301 - -snes_ksp_ew_threshold <threshold> - Sets threshold 4302 4303 Notes: 4304 Currently, the default is to use a constant relative tolerance for 4305 the inner linear solvers. Alternatively, one can use the 4306 Eisenstat-Walker method, where the relative convergence tolerance 4307 is reset at each Newton iteration according progress of the nonlinear 4308 solver. 4309 4310 Level: advanced 4311 4312 Reference: 4313 S. C. Eisenstat and H. F. Walker, "Choosing the forcing terms in an 4314 inexact Newton method", SISC 17 (1), pp.16-32, 1996. 4315 4316 .keywords: SNES, KSP, Eisenstat, Walker, convergence, test, inexact, Newton 4317 4318 .seealso: SNESKSPGetUseEW(), SNESKSPGetParametersEW(), SNESKSPSetParametersEW() 4319 @*/ 4320 PetscErrorCode SNESKSPSetUseEW(SNES snes,PetscBool flag) 4321 { 4322 PetscFunctionBegin; 4323 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4324 PetscValidLogicalCollectiveBool(snes,flag,2); 4325 snes->ksp_ewconv = flag; 4326 PetscFunctionReturn(0); 4327 } 4328 4329 #undef __FUNCT__ 4330 #define __FUNCT__ "SNESKSPGetUseEW" 4331 /*@ 4332 SNESKSPGetUseEW - Gets if SNES is using Eisenstat-Walker method 4333 for computing relative tolerance for linear solvers within an 4334 inexact Newton method. 4335 4336 Not Collective 4337 4338 Input Parameter: 4339 . snes - SNES context 4340 4341 Output Parameter: 4342 . flag - PETSC_TRUE or PETSC_FALSE 4343 4344 Notes: 4345 Currently, the default is to use a constant relative tolerance for 4346 the inner linear solvers. Alternatively, one can use the 4347 Eisenstat-Walker method, where the relative convergence tolerance 4348 is reset at each Newton iteration according progress of the nonlinear 4349 solver. 4350 4351 Level: advanced 4352 4353 Reference: 4354 S. C. Eisenstat and H. F. Walker, "Choosing the forcing terms in an 4355 inexact Newton method", SISC 17 (1), pp.16-32, 1996. 4356 4357 .keywords: SNES, KSP, Eisenstat, Walker, convergence, test, inexact, Newton 4358 4359 .seealso: SNESKSPSetUseEW(), SNESKSPGetParametersEW(), SNESKSPSetParametersEW() 4360 @*/ 4361 PetscErrorCode SNESKSPGetUseEW(SNES snes, PetscBool *flag) 4362 { 4363 PetscFunctionBegin; 4364 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4365 PetscValidPointer(flag,2); 4366 *flag = snes->ksp_ewconv; 4367 PetscFunctionReturn(0); 4368 } 4369 4370 #undef __FUNCT__ 4371 #define __FUNCT__ "SNESKSPSetParametersEW" 4372 /*@ 4373 SNESKSPSetParametersEW - Sets parameters for Eisenstat-Walker 4374 convergence criteria for the linear solvers within an inexact 4375 Newton method. 4376 4377 Logically Collective on SNES 4378 4379 Input Parameters: 4380 + snes - SNES context 4381 . version - version 1, 2 (default is 2) or 3 4382 . rtol_0 - initial relative tolerance (0 <= rtol_0 < 1) 4383 . rtol_max - maximum relative tolerance (0 <= rtol_max < 1) 4384 . gamma - multiplicative factor for version 2 rtol computation 4385 (0 <= gamma2 <= 1) 4386 . alpha - power for version 2 rtol computation (1 < alpha <= 2) 4387 . alpha2 - power for safeguard 4388 - threshold - threshold for imposing safeguard (0 < threshold < 1) 4389 4390 Note: 4391 Version 3 was contributed by Luis Chacon, June 2006. 4392 4393 Use PETSC_DEFAULT to retain the default for any of the parameters. 4394 4395 Level: advanced 4396 4397 Reference: 4398 S. C. Eisenstat and H. F. Walker, "Choosing the forcing terms in an 4399 inexact Newton method", Utah State University Math. Stat. Dept. Res. 4400 Report 6/94/75, June, 1994, to appear in SIAM J. Sci. Comput. 4401 4402 .keywords: SNES, KSP, Eisenstat, Walker, set, parameters 4403 4404 .seealso: SNESKSPSetUseEW(), SNESKSPGetUseEW(), SNESKSPGetParametersEW() 4405 @*/ 4406 PetscErrorCode SNESKSPSetParametersEW(SNES snes,PetscInt version,PetscReal rtol_0,PetscReal rtol_max,PetscReal gamma,PetscReal alpha,PetscReal alpha2,PetscReal threshold) 4407 { 4408 SNESKSPEW *kctx; 4409 4410 PetscFunctionBegin; 4411 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4412 kctx = (SNESKSPEW*)snes->kspconvctx; 4413 if (!kctx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"No Eisenstat-Walker context existing"); 4414 PetscValidLogicalCollectiveInt(snes,version,2); 4415 PetscValidLogicalCollectiveReal(snes,rtol_0,3); 4416 PetscValidLogicalCollectiveReal(snes,rtol_max,4); 4417 PetscValidLogicalCollectiveReal(snes,gamma,5); 4418 PetscValidLogicalCollectiveReal(snes,alpha,6); 4419 PetscValidLogicalCollectiveReal(snes,alpha2,7); 4420 PetscValidLogicalCollectiveReal(snes,threshold,8); 4421 4422 if (version != PETSC_DEFAULT) kctx->version = version; 4423 if (rtol_0 != PETSC_DEFAULT) kctx->rtol_0 = rtol_0; 4424 if (rtol_max != PETSC_DEFAULT) kctx->rtol_max = rtol_max; 4425 if (gamma != PETSC_DEFAULT) kctx->gamma = gamma; 4426 if (alpha != PETSC_DEFAULT) kctx->alpha = alpha; 4427 if (alpha2 != PETSC_DEFAULT) kctx->alpha2 = alpha2; 4428 if (threshold != PETSC_DEFAULT) kctx->threshold = threshold; 4429 4430 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); 4431 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); 4432 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); 4433 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); 4434 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); 4435 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); 4436 PetscFunctionReturn(0); 4437 } 4438 4439 #undef __FUNCT__ 4440 #define __FUNCT__ "SNESKSPGetParametersEW" 4441 /*@ 4442 SNESKSPGetParametersEW - Gets parameters for Eisenstat-Walker 4443 convergence criteria for the linear solvers within an inexact 4444 Newton method. 4445 4446 Not Collective 4447 4448 Input Parameters: 4449 snes - SNES context 4450 4451 Output Parameters: 4452 + version - version 1, 2 (default is 2) or 3 4453 . rtol_0 - initial relative tolerance (0 <= rtol_0 < 1) 4454 . rtol_max - maximum relative tolerance (0 <= rtol_max < 1) 4455 . gamma - multiplicative factor for version 2 rtol computation (0 <= gamma2 <= 1) 4456 . alpha - power for version 2 rtol computation (1 < alpha <= 2) 4457 . alpha2 - power for safeguard 4458 - threshold - threshold for imposing safeguard (0 < threshold < 1) 4459 4460 Level: advanced 4461 4462 .keywords: SNES, KSP, Eisenstat, Walker, get, parameters 4463 4464 .seealso: SNESKSPSetUseEW(), SNESKSPGetUseEW(), SNESKSPSetParametersEW() 4465 @*/ 4466 PetscErrorCode SNESKSPGetParametersEW(SNES snes,PetscInt *version,PetscReal *rtol_0,PetscReal *rtol_max,PetscReal *gamma,PetscReal *alpha,PetscReal *alpha2,PetscReal *threshold) 4467 { 4468 SNESKSPEW *kctx; 4469 4470 PetscFunctionBegin; 4471 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4472 kctx = (SNESKSPEW*)snes->kspconvctx; 4473 if (!kctx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"No Eisenstat-Walker context existing"); 4474 if (version) *version = kctx->version; 4475 if (rtol_0) *rtol_0 = kctx->rtol_0; 4476 if (rtol_max) *rtol_max = kctx->rtol_max; 4477 if (gamma) *gamma = kctx->gamma; 4478 if (alpha) *alpha = kctx->alpha; 4479 if (alpha2) *alpha2 = kctx->alpha2; 4480 if (threshold) *threshold = kctx->threshold; 4481 PetscFunctionReturn(0); 4482 } 4483 4484 #undef __FUNCT__ 4485 #define __FUNCT__ "KSPPreSolve_SNESEW" 4486 PetscErrorCode KSPPreSolve_SNESEW(KSP ksp, Vec b, Vec x, SNES snes) 4487 { 4488 PetscErrorCode ierr; 4489 SNESKSPEW *kctx = (SNESKSPEW*)snes->kspconvctx; 4490 PetscReal rtol = PETSC_DEFAULT,stol; 4491 4492 PetscFunctionBegin; 4493 if (!snes->ksp_ewconv) PetscFunctionReturn(0); 4494 if (!snes->iter) { 4495 rtol = kctx->rtol_0; /* first time in, so use the original user rtol */ 4496 ierr = VecNorm(snes->vec_func,NORM_2,&kctx->norm_first);CHKERRQ(ierr); 4497 } 4498 else { 4499 if (kctx->version == 1) { 4500 rtol = (snes->norm - kctx->lresid_last)/kctx->norm_last; 4501 if (rtol < 0.0) rtol = -rtol; 4502 stol = PetscPowReal(kctx->rtol_last,kctx->alpha2); 4503 if (stol > kctx->threshold) rtol = PetscMax(rtol,stol); 4504 } else if (kctx->version == 2) { 4505 rtol = kctx->gamma * PetscPowReal(snes->norm/kctx->norm_last,kctx->alpha); 4506 stol = kctx->gamma * PetscPowReal(kctx->rtol_last,kctx->alpha); 4507 if (stol > kctx->threshold) rtol = PetscMax(rtol,stol); 4508 } else if (kctx->version == 3) { /* contributed by Luis Chacon, June 2006. */ 4509 rtol = kctx->gamma * PetscPowReal(snes->norm/kctx->norm_last,kctx->alpha); 4510 /* safeguard: avoid sharp decrease of rtol */ 4511 stol = kctx->gamma*PetscPowReal(kctx->rtol_last,kctx->alpha); 4512 stol = PetscMax(rtol,stol); 4513 rtol = PetscMin(kctx->rtol_0,stol); 4514 /* safeguard: avoid oversolving */ 4515 stol = kctx->gamma*(kctx->norm_first*snes->rtol)/snes->norm; 4516 stol = PetscMax(rtol,stol); 4517 rtol = PetscMin(kctx->rtol_0,stol); 4518 } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Only versions 1, 2 or 3 are supported: %D",kctx->version); 4519 } 4520 /* safeguard: avoid rtol greater than one */ 4521 rtol = PetscMin(rtol,kctx->rtol_max); 4522 ierr = KSPSetTolerances(ksp,rtol,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);CHKERRQ(ierr); 4523 ierr = PetscInfo3(snes,"iter %D, Eisenstat-Walker (version %D) KSP rtol=%G\n",snes->iter,kctx->version,rtol);CHKERRQ(ierr); 4524 PetscFunctionReturn(0); 4525 } 4526 4527 #undef __FUNCT__ 4528 #define __FUNCT__ "KSPPostSolve_SNESEW" 4529 PetscErrorCode KSPPostSolve_SNESEW(KSP ksp, Vec b, Vec x, SNES snes) 4530 { 4531 PetscErrorCode ierr; 4532 SNESKSPEW *kctx = (SNESKSPEW*)snes->kspconvctx; 4533 PCSide pcside; 4534 Vec lres; 4535 4536 PetscFunctionBegin; 4537 if (!snes->ksp_ewconv) PetscFunctionReturn(0); 4538 ierr = KSPGetTolerances(ksp,&kctx->rtol_last,0,0,0);CHKERRQ(ierr); 4539 ierr = SNESGetFunctionNorm(snes,&kctx->norm_last);CHKERRQ(ierr); 4540 if (kctx->version == 1) { 4541 ierr = KSPGetPCSide(ksp,&pcside);CHKERRQ(ierr); 4542 if (pcside == PC_RIGHT) { /* XXX Should we also test KSP_UNPRECONDITIONED_NORM ? */ 4543 /* KSP residual is true linear residual */ 4544 ierr = KSPGetResidualNorm(ksp,&kctx->lresid_last);CHKERRQ(ierr); 4545 } else { 4546 /* KSP residual is preconditioned residual */ 4547 /* compute true linear residual norm */ 4548 ierr = VecDuplicate(b,&lres);CHKERRQ(ierr); 4549 ierr = MatMult(snes->jacobian,x,lres);CHKERRQ(ierr); 4550 ierr = VecAYPX(lres,-1.0,b);CHKERRQ(ierr); 4551 ierr = VecNorm(lres,NORM_2,&kctx->lresid_last);CHKERRQ(ierr); 4552 ierr = VecDestroy(&lres);CHKERRQ(ierr); 4553 } 4554 } 4555 PetscFunctionReturn(0); 4556 } 4557 4558 #undef __FUNCT__ 4559 #define __FUNCT__ "SNESGetKSP" 4560 /*@ 4561 SNESGetKSP - Returns the KSP context for a SNES solver. 4562 4563 Not Collective, but if SNES object is parallel, then KSP object is parallel 4564 4565 Input Parameter: 4566 . snes - the SNES context 4567 4568 Output Parameter: 4569 . ksp - the KSP context 4570 4571 Notes: 4572 The user can then directly manipulate the KSP context to set various 4573 options, etc. Likewise, the user can then extract and manipulate the 4574 PC contexts as well. 4575 4576 Level: beginner 4577 4578 .keywords: SNES, nonlinear, get, KSP, context 4579 4580 .seealso: KSPGetPC(), SNESCreate(), KSPCreate(), SNESSetKSP() 4581 @*/ 4582 PetscErrorCode SNESGetKSP(SNES snes,KSP *ksp) 4583 { 4584 PetscErrorCode ierr; 4585 4586 PetscFunctionBegin; 4587 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4588 PetscValidPointer(ksp,2); 4589 4590 if (!snes->ksp) { 4591 ierr = KSPCreate(PetscObjectComm((PetscObject)snes),&snes->ksp);CHKERRQ(ierr); 4592 ierr = PetscObjectIncrementTabLevel((PetscObject)snes->ksp,(PetscObject)snes,1);CHKERRQ(ierr); 4593 ierr = PetscLogObjectParent((PetscObject)snes,(PetscObject)snes->ksp);CHKERRQ(ierr); 4594 4595 ierr = KSPSetPreSolve(snes->ksp,(PetscErrorCode (*)(KSP,Vec,Vec,void*))KSPPreSolve_SNESEW,snes);CHKERRQ(ierr); 4596 ierr = KSPSetPostSolve(snes->ksp,(PetscErrorCode (*)(KSP,Vec,Vec,void*))KSPPostSolve_SNESEW,snes);CHKERRQ(ierr); 4597 } 4598 *ksp = snes->ksp; 4599 PetscFunctionReturn(0); 4600 } 4601 4602 4603 #include <petsc-private/dmimpl.h> 4604 #undef __FUNCT__ 4605 #define __FUNCT__ "SNESSetDM" 4606 /*@ 4607 SNESSetDM - Sets the DM that may be used by some preconditioners 4608 4609 Logically Collective on SNES 4610 4611 Input Parameters: 4612 + snes - the preconditioner context 4613 - dm - the dm 4614 4615 Level: intermediate 4616 4617 .seealso: SNESGetDM(), KSPSetDM(), KSPGetDM() 4618 @*/ 4619 PetscErrorCode SNESSetDM(SNES snes,DM dm) 4620 { 4621 PetscErrorCode ierr; 4622 KSP ksp; 4623 DMSNES sdm; 4624 4625 PetscFunctionBegin; 4626 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4627 if (dm) {ierr = PetscObjectReference((PetscObject)dm);CHKERRQ(ierr);} 4628 if (snes->dm) { /* Move the DMSNES context over to the new DM unless the new DM already has one */ 4629 if (snes->dm->dmsnes && snes->dmAuto && !dm->dmsnes) { 4630 ierr = DMCopyDMSNES(snes->dm,dm);CHKERRQ(ierr); 4631 ierr = DMGetDMSNES(snes->dm,&sdm);CHKERRQ(ierr); 4632 if (sdm->originaldm == snes->dm) sdm->originaldm = dm; /* Grant write privileges to the replacement DM */ 4633 } 4634 ierr = DMDestroy(&snes->dm);CHKERRQ(ierr); 4635 } 4636 snes->dm = dm; 4637 snes->dmAuto = PETSC_FALSE; 4638 4639 ierr = SNESGetKSP(snes,&ksp);CHKERRQ(ierr); 4640 ierr = KSPSetDM(ksp,dm);CHKERRQ(ierr); 4641 ierr = KSPSetDMActive(ksp,PETSC_FALSE);CHKERRQ(ierr); 4642 if (snes->pc) { 4643 ierr = SNESSetDM(snes->pc, snes->dm);CHKERRQ(ierr); 4644 ierr = SNESSetPCSide(snes,snes->pcside);CHKERRQ(ierr); 4645 } 4646 PetscFunctionReturn(0); 4647 } 4648 4649 #undef __FUNCT__ 4650 #define __FUNCT__ "SNESGetDM" 4651 /*@ 4652 SNESGetDM - Gets the DM that may be used by some preconditioners 4653 4654 Not Collective but DM obtained is parallel on SNES 4655 4656 Input Parameter: 4657 . snes - the preconditioner context 4658 4659 Output Parameter: 4660 . dm - the dm 4661 4662 Level: intermediate 4663 4664 .seealso: SNESSetDM(), KSPSetDM(), KSPGetDM() 4665 @*/ 4666 PetscErrorCode SNESGetDM(SNES snes,DM *dm) 4667 { 4668 PetscErrorCode ierr; 4669 4670 PetscFunctionBegin; 4671 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4672 if (!snes->dm) { 4673 ierr = DMShellCreate(PetscObjectComm((PetscObject)snes),&snes->dm);CHKERRQ(ierr); 4674 snes->dmAuto = PETSC_TRUE; 4675 } 4676 *dm = snes->dm; 4677 PetscFunctionReturn(0); 4678 } 4679 4680 #undef __FUNCT__ 4681 #define __FUNCT__ "SNESSetPC" 4682 /*@ 4683 SNESSetPC - Sets the nonlinear preconditioner to be used. 4684 4685 Collective on SNES 4686 4687 Input Parameters: 4688 + snes - iterative context obtained from SNESCreate() 4689 - pc - the preconditioner object 4690 4691 Notes: 4692 Use SNESGetPC() to retrieve the preconditioner context (for example, 4693 to configure it using the API). 4694 4695 Level: developer 4696 4697 .keywords: SNES, set, precondition 4698 .seealso: SNESGetPC() 4699 @*/ 4700 PetscErrorCode SNESSetPC(SNES snes, SNES pc) 4701 { 4702 PetscErrorCode ierr; 4703 4704 PetscFunctionBegin; 4705 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); 4706 PetscValidHeaderSpecific(pc, SNES_CLASSID, 2); 4707 PetscCheckSameComm(snes, 1, pc, 2); 4708 ierr = PetscObjectReference((PetscObject) pc);CHKERRQ(ierr); 4709 ierr = SNESDestroy(&snes->pc);CHKERRQ(ierr); 4710 snes->pc = pc; 4711 ierr = PetscLogObjectParent((PetscObject)snes, (PetscObject)snes->pc);CHKERRQ(ierr); 4712 PetscFunctionReturn(0); 4713 } 4714 4715 #undef __FUNCT__ 4716 #define __FUNCT__ "SNESGetPC" 4717 /*@ 4718 SNESGetPC - Returns a pointer to the nonlinear preconditioning context set with SNESSetPC(). 4719 4720 Not Collective 4721 4722 Input Parameter: 4723 . snes - iterative context obtained from SNESCreate() 4724 4725 Output Parameter: 4726 . pc - preconditioner context 4727 4728 Level: developer 4729 4730 .keywords: SNES, get, preconditioner 4731 .seealso: SNESSetPC() 4732 @*/ 4733 PetscErrorCode SNESGetPC(SNES snes, SNES *pc) 4734 { 4735 PetscErrorCode ierr; 4736 const char *optionsprefix; 4737 4738 PetscFunctionBegin; 4739 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); 4740 PetscValidPointer(pc, 2); 4741 if (!snes->pc) { 4742 ierr = SNESCreate(PetscObjectComm((PetscObject)snes),&snes->pc);CHKERRQ(ierr); 4743 ierr = PetscObjectIncrementTabLevel((PetscObject)snes->pc,(PetscObject)snes,1);CHKERRQ(ierr); 4744 ierr = PetscLogObjectParent((PetscObject)snes,(PetscObject)snes->pc);CHKERRQ(ierr); 4745 ierr = SNESGetOptionsPrefix(snes,&optionsprefix);CHKERRQ(ierr); 4746 ierr = SNESSetOptionsPrefix(snes->pc,optionsprefix);CHKERRQ(ierr); 4747 ierr = SNESAppendOptionsPrefix(snes->pc,"npc_");CHKERRQ(ierr); 4748 ierr = SNESSetCountersReset(snes->pc,PETSC_FALSE);CHKERRQ(ierr); 4749 } 4750 *pc = snes->pc; 4751 PetscFunctionReturn(0); 4752 } 4753 4754 #undef __FUNCT__ 4755 #define __FUNCT__ "SNESSetPCSide" 4756 /*@ 4757 SNESSetPCSide - Sets the preconditioning side. 4758 4759 Logically Collective on SNES 4760 4761 Input Parameter: 4762 . snes - iterative context obtained from SNESCreate() 4763 4764 Output Parameter: 4765 . side - the preconditioning side, where side is one of 4766 .vb 4767 PC_LEFT - left preconditioning (default) 4768 PC_RIGHT - right preconditioning 4769 .ve 4770 4771 Options Database Keys: 4772 . -snes_pc_side <right,left> 4773 4774 Level: intermediate 4775 4776 .keywords: SNES, set, right, left, side, preconditioner, flag 4777 4778 .seealso: SNESGetPCSide(), KSPSetPCSide() 4779 @*/ 4780 PetscErrorCode SNESSetPCSide(SNES snes,PCSide side) 4781 { 4782 PetscFunctionBegin; 4783 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4784 PetscValidLogicalCollectiveEnum(snes,side,2); 4785 snes->pcside = side; 4786 PetscFunctionReturn(0); 4787 } 4788 4789 #undef __FUNCT__ 4790 #define __FUNCT__ "SNESGetPCSide" 4791 /*@ 4792 SNESGetPCSide - Gets the preconditioning side. 4793 4794 Not Collective 4795 4796 Input Parameter: 4797 . snes - iterative context obtained from SNESCreate() 4798 4799 Output Parameter: 4800 . side - the preconditioning side, where side is one of 4801 .vb 4802 PC_LEFT - left preconditioning (default) 4803 PC_RIGHT - right preconditioning 4804 .ve 4805 4806 Level: intermediate 4807 4808 .keywords: SNES, get, right, left, side, preconditioner, flag 4809 4810 .seealso: SNESSetPCSide(), KSPGetPCSide() 4811 @*/ 4812 PetscErrorCode SNESGetPCSide(SNES snes,PCSide *side) 4813 { 4814 PetscFunctionBegin; 4815 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4816 PetscValidPointer(side,2); 4817 *side = snes->pcside; 4818 PetscFunctionReturn(0); 4819 } 4820 4821 #undef __FUNCT__ 4822 #define __FUNCT__ "SNESSetLineSearch" 4823 /*@ 4824 SNESSetLineSearch - Sets the linesearch on the SNES instance. 4825 4826 Collective on SNES 4827 4828 Input Parameters: 4829 + snes - iterative context obtained from SNESCreate() 4830 - linesearch - the linesearch object 4831 4832 Notes: 4833 Use SNESGetLineSearch() to retrieve the preconditioner context (for example, 4834 to configure it using the API). 4835 4836 Level: developer 4837 4838 .keywords: SNES, set, linesearch 4839 .seealso: SNESGetLineSearch() 4840 @*/ 4841 PetscErrorCode SNESSetLineSearch(SNES snes, SNESLineSearch linesearch) 4842 { 4843 PetscErrorCode ierr; 4844 4845 PetscFunctionBegin; 4846 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); 4847 PetscValidHeaderSpecific(linesearch, SNESLINESEARCH_CLASSID, 2); 4848 PetscCheckSameComm(snes, 1, linesearch, 2); 4849 ierr = PetscObjectReference((PetscObject) linesearch);CHKERRQ(ierr); 4850 ierr = SNESLineSearchDestroy(&snes->linesearch);CHKERRQ(ierr); 4851 4852 snes->linesearch = linesearch; 4853 4854 ierr = PetscLogObjectParent((PetscObject)snes, (PetscObject)snes->linesearch);CHKERRQ(ierr); 4855 PetscFunctionReturn(0); 4856 } 4857 4858 #undef __FUNCT__ 4859 #define __FUNCT__ "SNESGetLineSearch" 4860 /*@ 4861 SNESGetLineSearch - Returns a pointer to the line search context set with SNESSetLineSearch() 4862 or creates a default line search instance associated with the SNES and returns it. 4863 4864 Not Collective 4865 4866 Input Parameter: 4867 . snes - iterative context obtained from SNESCreate() 4868 4869 Output Parameter: 4870 . linesearch - linesearch context 4871 4872 Level: developer 4873 4874 .keywords: SNES, get, linesearch 4875 .seealso: SNESSetLineSearch() 4876 @*/ 4877 PetscErrorCode SNESGetLineSearch(SNES snes, SNESLineSearch *linesearch) 4878 { 4879 PetscErrorCode ierr; 4880 const char *optionsprefix; 4881 4882 PetscFunctionBegin; 4883 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); 4884 PetscValidPointer(linesearch, 2); 4885 if (!snes->linesearch) { 4886 ierr = SNESGetOptionsPrefix(snes, &optionsprefix);CHKERRQ(ierr); 4887 ierr = SNESLineSearchCreate(PetscObjectComm((PetscObject)snes), &snes->linesearch);CHKERRQ(ierr); 4888 ierr = SNESLineSearchSetSNES(snes->linesearch, snes);CHKERRQ(ierr); 4889 ierr = SNESLineSearchAppendOptionsPrefix(snes->linesearch, optionsprefix);CHKERRQ(ierr); 4890 ierr = PetscObjectIncrementTabLevel((PetscObject) snes->linesearch, (PetscObject) snes, 1);CHKERRQ(ierr); 4891 ierr = PetscLogObjectParent((PetscObject)snes, (PetscObject)snes->linesearch);CHKERRQ(ierr); 4892 } 4893 *linesearch = snes->linesearch; 4894 PetscFunctionReturn(0); 4895 } 4896 4897 #if defined(PETSC_HAVE_MATLAB_ENGINE) 4898 #include <mex.h> 4899 4900 typedef struct {char *funcname; mxArray *ctx;} SNESMatlabContext; 4901 4902 #undef __FUNCT__ 4903 #define __FUNCT__ "SNESComputeFunction_Matlab" 4904 /* 4905 SNESComputeFunction_Matlab - Calls the function that has been set with SNESSetFunctionMatlab(). 4906 4907 Collective on SNES 4908 4909 Input Parameters: 4910 + snes - the SNES context 4911 - x - input vector 4912 4913 Output Parameter: 4914 . y - function vector, as set by SNESSetFunction() 4915 4916 Notes: 4917 SNESComputeFunction() is typically used within nonlinear solvers 4918 implementations, so most users would not generally call this routine 4919 themselves. 4920 4921 Level: developer 4922 4923 .keywords: SNES, nonlinear, compute, function 4924 4925 .seealso: SNESSetFunction(), SNESGetFunction() 4926 */ 4927 PetscErrorCode SNESComputeFunction_Matlab(SNES snes,Vec x,Vec y, void *ctx) 4928 { 4929 PetscErrorCode ierr; 4930 SNESMatlabContext *sctx = (SNESMatlabContext*)ctx; 4931 int nlhs = 1,nrhs = 5; 4932 mxArray *plhs[1],*prhs[5]; 4933 long long int lx = 0,ly = 0,ls = 0; 4934 4935 PetscFunctionBegin; 4936 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4937 PetscValidHeaderSpecific(x,VEC_CLASSID,2); 4938 PetscValidHeaderSpecific(y,VEC_CLASSID,3); 4939 PetscCheckSameComm(snes,1,x,2); 4940 PetscCheckSameComm(snes,1,y,3); 4941 4942 /* call Matlab function in ctx with arguments x and y */ 4943 4944 ierr = PetscMemcpy(&ls,&snes,sizeof(snes));CHKERRQ(ierr); 4945 ierr = PetscMemcpy(&lx,&x,sizeof(x));CHKERRQ(ierr); 4946 ierr = PetscMemcpy(&ly,&y,sizeof(x));CHKERRQ(ierr); 4947 prhs[0] = mxCreateDoubleScalar((double)ls); 4948 prhs[1] = mxCreateDoubleScalar((double)lx); 4949 prhs[2] = mxCreateDoubleScalar((double)ly); 4950 prhs[3] = mxCreateString(sctx->funcname); 4951 prhs[4] = sctx->ctx; 4952 ierr = mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscSNESComputeFunctionInternal");CHKERRQ(ierr); 4953 ierr = mxGetScalar(plhs[0]);CHKERRQ(ierr); 4954 mxDestroyArray(prhs[0]); 4955 mxDestroyArray(prhs[1]); 4956 mxDestroyArray(prhs[2]); 4957 mxDestroyArray(prhs[3]); 4958 mxDestroyArray(plhs[0]); 4959 PetscFunctionReturn(0); 4960 } 4961 4962 #undef __FUNCT__ 4963 #define __FUNCT__ "SNESSetFunctionMatlab" 4964 /* 4965 SNESSetFunctionMatlab - Sets the function evaluation routine and function 4966 vector for use by the SNES routines in solving systems of nonlinear 4967 equations from MATLAB. Here the function is a string containing the name of a MATLAB function 4968 4969 Logically Collective on SNES 4970 4971 Input Parameters: 4972 + snes - the SNES context 4973 . r - vector to store function value 4974 - SNESFunction - function evaluation routine 4975 4976 Notes: 4977 The Newton-like methods typically solve linear systems of the form 4978 $ f'(x) x = -f(x), 4979 where f'(x) denotes the Jacobian matrix and f(x) is the function. 4980 4981 Level: beginner 4982 4983 Developer Note: This bleeds the allocated memory SNESMatlabContext *sctx; 4984 4985 .keywords: SNES, nonlinear, set, function 4986 4987 .seealso: SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESSetFunction() 4988 */ 4989 PetscErrorCode SNESSetFunctionMatlab(SNES snes,Vec r,const char *SNESFunction,mxArray *ctx) 4990 { 4991 PetscErrorCode ierr; 4992 SNESMatlabContext *sctx; 4993 4994 PetscFunctionBegin; 4995 /* currently sctx is memory bleed */ 4996 ierr = PetscMalloc(sizeof(SNESMatlabContext),&sctx);CHKERRQ(ierr); 4997 ierr = PetscStrallocpy(SNESFunction,&sctx->funcname);CHKERRQ(ierr); 4998 /* 4999 This should work, but it doesn't 5000 sctx->ctx = ctx; 5001 mexMakeArrayPersistent(sctx->ctx); 5002 */ 5003 sctx->ctx = mxDuplicateArray(ctx); 5004 ierr = SNESSetFunction(snes,r,SNESComputeFunction_Matlab,sctx);CHKERRQ(ierr); 5005 PetscFunctionReturn(0); 5006 } 5007 5008 #undef __FUNCT__ 5009 #define __FUNCT__ "SNESComputeJacobian_Matlab" 5010 /* 5011 SNESComputeJacobian_Matlab - Calls the function that has been set with SNESSetJacobianMatlab(). 5012 5013 Collective on SNES 5014 5015 Input Parameters: 5016 + snes - the SNES context 5017 . x - input vector 5018 . A, B - the matrices 5019 - ctx - user context 5020 5021 Output Parameter: 5022 . flag - structure of the matrix 5023 5024 Level: developer 5025 5026 .keywords: SNES, nonlinear, compute, function 5027 5028 .seealso: SNESSetFunction(), SNESGetFunction() 5029 @*/ 5030 PetscErrorCode SNESComputeJacobian_Matlab(SNES snes,Vec x,Mat *A,Mat *B,MatStructure *flag, void *ctx) 5031 { 5032 PetscErrorCode ierr; 5033 SNESMatlabContext *sctx = (SNESMatlabContext*)ctx; 5034 int nlhs = 2,nrhs = 6; 5035 mxArray *plhs[2],*prhs[6]; 5036 long long int lx = 0,lA = 0,ls = 0, lB = 0; 5037 5038 PetscFunctionBegin; 5039 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 5040 PetscValidHeaderSpecific(x,VEC_CLASSID,2); 5041 5042 /* call Matlab function in ctx with arguments x and y */ 5043 5044 ierr = PetscMemcpy(&ls,&snes,sizeof(snes));CHKERRQ(ierr); 5045 ierr = PetscMemcpy(&lx,&x,sizeof(x));CHKERRQ(ierr); 5046 ierr = PetscMemcpy(&lA,A,sizeof(x));CHKERRQ(ierr); 5047 ierr = PetscMemcpy(&lB,B,sizeof(x));CHKERRQ(ierr); 5048 prhs[0] = mxCreateDoubleScalar((double)ls); 5049 prhs[1] = mxCreateDoubleScalar((double)lx); 5050 prhs[2] = mxCreateDoubleScalar((double)lA); 5051 prhs[3] = mxCreateDoubleScalar((double)lB); 5052 prhs[4] = mxCreateString(sctx->funcname); 5053 prhs[5] = sctx->ctx; 5054 ierr = mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscSNESComputeJacobianInternal");CHKERRQ(ierr); 5055 ierr = mxGetScalar(plhs[0]);CHKERRQ(ierr); 5056 *flag = (MatStructure) mxGetScalar(plhs[1]);CHKERRQ(ierr); 5057 mxDestroyArray(prhs[0]); 5058 mxDestroyArray(prhs[1]); 5059 mxDestroyArray(prhs[2]); 5060 mxDestroyArray(prhs[3]); 5061 mxDestroyArray(prhs[4]); 5062 mxDestroyArray(plhs[0]); 5063 mxDestroyArray(plhs[1]); 5064 PetscFunctionReturn(0); 5065 } 5066 5067 #undef __FUNCT__ 5068 #define __FUNCT__ "SNESSetJacobianMatlab" 5069 /* 5070 SNESSetJacobianMatlab - Sets the Jacobian function evaluation routine and two empty Jacobian matrices 5071 vector for use by the SNES routines in solving systems of nonlinear 5072 equations from MATLAB. Here the function is a string containing the name of a MATLAB function 5073 5074 Logically Collective on SNES 5075 5076 Input Parameters: 5077 + snes - the SNES context 5078 . A,B - Jacobian matrices 5079 . SNESJacobianFunction - function evaluation routine 5080 - ctx - user context 5081 5082 Level: developer 5083 5084 Developer Note: This bleeds the allocated memory SNESMatlabContext *sctx; 5085 5086 .keywords: SNES, nonlinear, set, function 5087 5088 .seealso: SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESSetFunction(), SNESJacobianFunction 5089 */ 5090 PetscErrorCode SNESSetJacobianMatlab(SNES snes,Mat A,Mat B,const char *SNESJacobianFunction,mxArray *ctx) 5091 { 5092 PetscErrorCode ierr; 5093 SNESMatlabContext *sctx; 5094 5095 PetscFunctionBegin; 5096 /* currently sctx is memory bleed */ 5097 ierr = PetscMalloc(sizeof(SNESMatlabContext),&sctx);CHKERRQ(ierr); 5098 ierr = PetscStrallocpy(SNESJacobianFunction,&sctx->funcname);CHKERRQ(ierr); 5099 /* 5100 This should work, but it doesn't 5101 sctx->ctx = ctx; 5102 mexMakeArrayPersistent(sctx->ctx); 5103 */ 5104 sctx->ctx = mxDuplicateArray(ctx); 5105 ierr = SNESSetJacobian(snes,A,B,SNESComputeJacobian_Matlab,sctx);CHKERRQ(ierr); 5106 PetscFunctionReturn(0); 5107 } 5108 5109 #undef __FUNCT__ 5110 #define __FUNCT__ "SNESMonitor_Matlab" 5111 /* 5112 SNESMonitor_Matlab - Calls the function that has been set with SNESMonitorSetMatlab(). 5113 5114 Collective on SNES 5115 5116 .seealso: SNESSetFunction(), SNESGetFunction() 5117 @*/ 5118 PetscErrorCode SNESMonitor_Matlab(SNES snes,PetscInt it, PetscReal fnorm, void *ctx) 5119 { 5120 PetscErrorCode ierr; 5121 SNESMatlabContext *sctx = (SNESMatlabContext*)ctx; 5122 int nlhs = 1,nrhs = 6; 5123 mxArray *plhs[1],*prhs[6]; 5124 long long int lx = 0,ls = 0; 5125 Vec x = snes->vec_sol; 5126 5127 PetscFunctionBegin; 5128 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 5129 5130 ierr = PetscMemcpy(&ls,&snes,sizeof(snes));CHKERRQ(ierr); 5131 ierr = PetscMemcpy(&lx,&x,sizeof(x));CHKERRQ(ierr); 5132 prhs[0] = mxCreateDoubleScalar((double)ls); 5133 prhs[1] = mxCreateDoubleScalar((double)it); 5134 prhs[2] = mxCreateDoubleScalar((double)fnorm); 5135 prhs[3] = mxCreateDoubleScalar((double)lx); 5136 prhs[4] = mxCreateString(sctx->funcname); 5137 prhs[5] = sctx->ctx; 5138 ierr = mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscSNESMonitorInternal");CHKERRQ(ierr); 5139 ierr = mxGetScalar(plhs[0]);CHKERRQ(ierr); 5140 mxDestroyArray(prhs[0]); 5141 mxDestroyArray(prhs[1]); 5142 mxDestroyArray(prhs[2]); 5143 mxDestroyArray(prhs[3]); 5144 mxDestroyArray(prhs[4]); 5145 mxDestroyArray(plhs[0]); 5146 PetscFunctionReturn(0); 5147 } 5148 5149 #undef __FUNCT__ 5150 #define __FUNCT__ "SNESMonitorSetMatlab" 5151 /* 5152 SNESMonitorSetMatlab - Sets the monitor function from MATLAB 5153 5154 Level: developer 5155 5156 Developer Note: This bleeds the allocated memory SNESMatlabContext *sctx; 5157 5158 .keywords: SNES, nonlinear, set, function 5159 5160 .seealso: SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESSetFunction() 5161 */ 5162 PetscErrorCode SNESMonitorSetMatlab(SNES snes,const char *SNESMonitorFunction,mxArray *ctx) 5163 { 5164 PetscErrorCode ierr; 5165 SNESMatlabContext *sctx; 5166 5167 PetscFunctionBegin; 5168 /* currently sctx is memory bleed */ 5169 ierr = PetscMalloc(sizeof(SNESMatlabContext),&sctx);CHKERRQ(ierr); 5170 ierr = PetscStrallocpy(SNESMonitorFunction,&sctx->funcname);CHKERRQ(ierr); 5171 /* 5172 This should work, but it doesn't 5173 sctx->ctx = ctx; 5174 mexMakeArrayPersistent(sctx->ctx); 5175 */ 5176 sctx->ctx = mxDuplicateArray(ctx); 5177 ierr = SNESMonitorSet(snes,SNESMonitor_Matlab,sctx,NULL);CHKERRQ(ierr); 5178 PetscFunctionReturn(0); 5179 } 5180 5181 #endif 5182