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