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