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