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