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