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