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