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