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 Mat j,jpre; 2532 2533 PetscFunctionBegin; 2534 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2535 if (snes->setupcalled) PetscFunctionReturn(0); 2536 2537 if (!((PetscObject)snes)->type_name) { 2538 ierr = SNESSetType(snes,SNESNEWTONLS);CHKERRQ(ierr); 2539 } 2540 2541 ierr = SNESGetFunction(snes,&snes->vec_func,NULL,NULL);CHKERRQ(ierr); 2542 2543 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 2544 ierr = DMGetDMSNES(dm,&sdm);CHKERRQ(ierr); 2545 if (!sdm->ops->computefunction) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"Function never provided to SNES object"); 2546 if (!sdm->ops->computejacobian) { 2547 ierr = DMSNESSetJacobian(dm,SNESComputeJacobianDefaultColor,NULL);CHKERRQ(ierr); 2548 } 2549 if (!snes->vec_func) { 2550 ierr = DMCreateGlobalVector(dm,&snes->vec_func);CHKERRQ(ierr); 2551 } 2552 2553 if (!snes->ksp) { 2554 ierr = SNESGetKSP(snes, &snes->ksp);CHKERRQ(ierr); 2555 } 2556 2557 if (!snes->linesearch) { 2558 ierr = SNESGetLineSearch(snes, &snes->linesearch);CHKERRQ(ierr); 2559 } 2560 ierr = SNESLineSearchSetFunction(snes->linesearch,SNESComputeFunction);CHKERRQ(ierr); 2561 2562 if (snes->pc && (snes->pcside == PC_LEFT)) { 2563 snes->mf = PETSC_TRUE; 2564 snes->mf_operator = PETSC_FALSE; 2565 } 2566 2567 if (snes->pc) { 2568 /* copy the DM over */ 2569 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 2570 ierr = SNESSetDM(snes->pc,dm);CHKERRQ(ierr); 2571 2572 ierr = SNESGetFunction(snes,&f,&func,&funcctx);CHKERRQ(ierr); 2573 ierr = VecDuplicate(f,&fpc);CHKERRQ(ierr); 2574 ierr = SNESSetFunction(snes->pc,fpc,func,funcctx);CHKERRQ(ierr); 2575 ierr = SNESGetJacobian(snes,&j,&jpre,&jac,&jacctx);CHKERRQ(ierr); 2576 ierr = SNESSetJacobian(snes->pc,j,jpre,jac,jacctx);CHKERRQ(ierr); 2577 ierr = SNESGetApplicationContext(snes,&appctx);CHKERRQ(ierr); 2578 ierr = SNESSetApplicationContext(snes->pc,appctx);CHKERRQ(ierr); 2579 ierr = VecDestroy(&fpc);CHKERRQ(ierr); 2580 2581 /* copy the function pointers over */ 2582 ierr = PetscObjectCopyFortranFunctionPointers((PetscObject)snes,(PetscObject)snes->pc);CHKERRQ(ierr); 2583 2584 /* default to 1 iteration */ 2585 ierr = SNESSetTolerances(snes->pc,0.0,0.0,0.0,1,snes->pc->max_funcs);CHKERRQ(ierr); 2586 if (snes->pcside==PC_RIGHT) { 2587 ierr = SNESSetNormSchedule(snes->pc,SNES_NORM_FINAL_ONLY);CHKERRQ(ierr); 2588 } else { 2589 ierr = SNESSetNormSchedule(snes->pc,SNES_NORM_NONE);CHKERRQ(ierr); 2590 } 2591 ierr = SNESSetFromOptions(snes->pc);CHKERRQ(ierr); 2592 2593 /* copy the line search context over */ 2594 ierr = SNESGetLineSearch(snes,&linesearch);CHKERRQ(ierr); 2595 ierr = SNESGetLineSearch(snes->pc,&pclinesearch);CHKERRQ(ierr); 2596 ierr = SNESLineSearchGetPreCheck(linesearch,&precheck,&lsprectx);CHKERRQ(ierr); 2597 ierr = SNESLineSearchGetPostCheck(linesearch,&postcheck,&lspostctx);CHKERRQ(ierr); 2598 ierr = SNESLineSearchSetPreCheck(pclinesearch,precheck,lsprectx);CHKERRQ(ierr); 2599 ierr = SNESLineSearchSetPostCheck(pclinesearch,postcheck,lspostctx);CHKERRQ(ierr); 2600 ierr = PetscObjectCopyFortranFunctionPointers((PetscObject)linesearch, (PetscObject)pclinesearch);CHKERRQ(ierr); 2601 } 2602 if (snes->mf) { 2603 ierr = SNESSetUpMatrixFree_Private(snes, snes->mf_operator, snes->mf_version);CHKERRQ(ierr); 2604 } 2605 if (snes->ops->usercompute && !snes->user) { 2606 ierr = (*snes->ops->usercompute)(snes,(void**)&snes->user);CHKERRQ(ierr); 2607 } 2608 2609 snes->jac_iter = 0; 2610 snes->pre_iter = 0; 2611 2612 if (snes->ops->setup) { 2613 ierr = (*snes->ops->setup)(snes);CHKERRQ(ierr); 2614 } 2615 2616 if (snes->pc && (snes->pcside == PC_LEFT)) { 2617 if (snes->functype == SNES_FUNCTION_PRECONDITIONED) { 2618 ierr = SNESGetLineSearch(snes,&linesearch);CHKERRQ(ierr); 2619 ierr = SNESLineSearchSetFunction(linesearch,SNESComputeFunctionDefaultNPC);CHKERRQ(ierr); 2620 } 2621 } 2622 2623 snes->setupcalled = PETSC_TRUE; 2624 PetscFunctionReturn(0); 2625 } 2626 2627 #undef __FUNCT__ 2628 #define __FUNCT__ "SNESReset" 2629 /*@ 2630 SNESReset - Resets a SNES context to the snessetupcalled = 0 state and removes any allocated Vecs and Mats 2631 2632 Collective on SNES 2633 2634 Input Parameter: 2635 . snes - iterative context obtained from SNESCreate() 2636 2637 Level: intermediate 2638 2639 Notes: Also calls the application context destroy routine set with SNESSetComputeApplicationContext() 2640 2641 .keywords: SNES, destroy 2642 2643 .seealso: SNESCreate(), SNESSetUp(), SNESSolve() 2644 @*/ 2645 PetscErrorCode SNESReset(SNES snes) 2646 { 2647 PetscErrorCode ierr; 2648 2649 PetscFunctionBegin; 2650 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2651 if (snes->ops->userdestroy && snes->user) { 2652 ierr = (*snes->ops->userdestroy)((void**)&snes->user);CHKERRQ(ierr); 2653 snes->user = NULL; 2654 } 2655 if (snes->pc) { 2656 ierr = SNESReset(snes->pc);CHKERRQ(ierr); 2657 } 2658 2659 if (snes->ops->reset) { 2660 ierr = (*snes->ops->reset)(snes);CHKERRQ(ierr); 2661 } 2662 if (snes->ksp) { 2663 ierr = KSPReset(snes->ksp);CHKERRQ(ierr); 2664 } 2665 2666 if (snes->linesearch) { 2667 ierr = SNESLineSearchReset(snes->linesearch);CHKERRQ(ierr); 2668 } 2669 2670 ierr = VecDestroy(&snes->vec_rhs);CHKERRQ(ierr); 2671 ierr = VecDestroy(&snes->vec_sol);CHKERRQ(ierr); 2672 ierr = VecDestroy(&snes->vec_sol_update);CHKERRQ(ierr); 2673 ierr = VecDestroy(&snes->vec_func);CHKERRQ(ierr); 2674 ierr = MatDestroy(&snes->jacobian);CHKERRQ(ierr); 2675 ierr = MatDestroy(&snes->jacobian_pre);CHKERRQ(ierr); 2676 ierr = VecDestroyVecs(snes->nwork,&snes->work);CHKERRQ(ierr); 2677 ierr = VecDestroyVecs(snes->nvwork,&snes->vwork);CHKERRQ(ierr); 2678 2679 snes->nwork = snes->nvwork = 0; 2680 snes->setupcalled = PETSC_FALSE; 2681 PetscFunctionReturn(0); 2682 } 2683 2684 #undef __FUNCT__ 2685 #define __FUNCT__ "SNESDestroy" 2686 /*@ 2687 SNESDestroy - Destroys the nonlinear solver context that was created 2688 with SNESCreate(). 2689 2690 Collective on SNES 2691 2692 Input Parameter: 2693 . snes - the SNES context 2694 2695 Level: beginner 2696 2697 .keywords: SNES, nonlinear, destroy 2698 2699 .seealso: SNESCreate(), SNESSolve() 2700 @*/ 2701 PetscErrorCode SNESDestroy(SNES *snes) 2702 { 2703 PetscErrorCode ierr; 2704 2705 PetscFunctionBegin; 2706 if (!*snes) PetscFunctionReturn(0); 2707 PetscValidHeaderSpecific((*snes),SNES_CLASSID,1); 2708 if (--((PetscObject)(*snes))->refct > 0) {*snes = 0; PetscFunctionReturn(0);} 2709 2710 ierr = SNESReset((*snes));CHKERRQ(ierr); 2711 ierr = SNESDestroy(&(*snes)->pc);CHKERRQ(ierr); 2712 2713 /* if memory was published with SAWs then destroy it */ 2714 ierr = PetscObjectSAWsViewOff((PetscObject)*snes);CHKERRQ(ierr); 2715 if ((*snes)->ops->destroy) {ierr = (*((*snes))->ops->destroy)((*snes));CHKERRQ(ierr);} 2716 2717 ierr = DMDestroy(&(*snes)->dm);CHKERRQ(ierr); 2718 ierr = KSPDestroy(&(*snes)->ksp);CHKERRQ(ierr); 2719 ierr = SNESLineSearchDestroy(&(*snes)->linesearch);CHKERRQ(ierr); 2720 2721 ierr = PetscFree((*snes)->kspconvctx);CHKERRQ(ierr); 2722 if ((*snes)->ops->convergeddestroy) { 2723 ierr = (*(*snes)->ops->convergeddestroy)((*snes)->cnvP);CHKERRQ(ierr); 2724 } 2725 if ((*snes)->conv_malloc) { 2726 ierr = PetscFree((*snes)->conv_hist);CHKERRQ(ierr); 2727 ierr = PetscFree((*snes)->conv_hist_its);CHKERRQ(ierr); 2728 } 2729 ierr = SNESMonitorCancel((*snes));CHKERRQ(ierr); 2730 ierr = PetscHeaderDestroy(snes);CHKERRQ(ierr); 2731 PetscFunctionReturn(0); 2732 } 2733 2734 /* ----------- Routines to set solver parameters ---------- */ 2735 2736 #undef __FUNCT__ 2737 #define __FUNCT__ "SNESSetLagPreconditioner" 2738 /*@ 2739 SNESSetLagPreconditioner - Determines when the preconditioner is rebuilt in the nonlinear solve. 2740 2741 Logically Collective on SNES 2742 2743 Input Parameters: 2744 + snes - the SNES context 2745 - lag - -1 indicates NEVER rebuild, 1 means rebuild every time the Jacobian is computed within a single nonlinear solve, 2 means every second time 2746 the Jacobian is built etc. -2 indicates rebuild preconditioner at next chance but then never rebuild after that 2747 2748 Options Database Keys: 2749 . -snes_lag_preconditioner <lag> 2750 2751 Notes: 2752 The default is 1 2753 The preconditioner is ALWAYS built in the first iteration of a nonlinear solve unless lag is -1 2754 If -1 is used before the very first nonlinear solve the preconditioner is still built because there is no previous preconditioner to use 2755 2756 Level: intermediate 2757 2758 .keywords: SNES, nonlinear, set, convergence, tolerances 2759 2760 .seealso: SNESSetTrustRegionTolerance(), SNESGetLagPreconditioner(), SNESSetLagJacobian(), SNESGetLagJacobian() 2761 2762 @*/ 2763 PetscErrorCode SNESSetLagPreconditioner(SNES snes,PetscInt lag) 2764 { 2765 PetscFunctionBegin; 2766 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2767 if (lag < -2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Lag must be -2, -1, 1 or greater"); 2768 if (!lag) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Lag cannot be 0"); 2769 PetscValidLogicalCollectiveInt(snes,lag,2); 2770 snes->lagpreconditioner = lag; 2771 PetscFunctionReturn(0); 2772 } 2773 2774 #undef __FUNCT__ 2775 #define __FUNCT__ "SNESSetGridSequence" 2776 /*@ 2777 SNESSetGridSequence - sets the number of steps of grid sequencing that SNES does 2778 2779 Logically Collective on SNES 2780 2781 Input Parameters: 2782 + snes - the SNES context 2783 - steps - the number of refinements to do, defaults to 0 2784 2785 Options Database Keys: 2786 . -snes_grid_sequence <steps> 2787 2788 Level: intermediate 2789 2790 Notes: 2791 Use SNESGetSolution() to extract the fine grid solution after grid sequencing. 2792 2793 .keywords: SNES, nonlinear, set, convergence, tolerances 2794 2795 .seealso: SNESSetTrustRegionTolerance(), SNESGetLagPreconditioner(), SNESSetLagJacobian(), SNESGetLagJacobian() 2796 2797 @*/ 2798 PetscErrorCode SNESSetGridSequence(SNES snes,PetscInt steps) 2799 { 2800 PetscFunctionBegin; 2801 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2802 PetscValidLogicalCollectiveInt(snes,steps,2); 2803 snes->gridsequence = steps; 2804 PetscFunctionReturn(0); 2805 } 2806 2807 #undef __FUNCT__ 2808 #define __FUNCT__ "SNESGetLagPreconditioner" 2809 /*@ 2810 SNESGetLagPreconditioner - Indicates how often the preconditioner is rebuilt 2811 2812 Not Collective 2813 2814 Input Parameter: 2815 . snes - the SNES context 2816 2817 Output Parameter: 2818 . lag - -1 indicates NEVER rebuild, 1 means rebuild every time the Jacobian is computed within a single nonlinear solve, 2 means every second time 2819 the Jacobian is built etc. -2 indicates rebuild preconditioner at next chance but then never rebuild after that 2820 2821 Options Database Keys: 2822 . -snes_lag_preconditioner <lag> 2823 2824 Notes: 2825 The default is 1 2826 The preconditioner is ALWAYS built in the first iteration of a nonlinear solve unless lag is -1 2827 2828 Level: intermediate 2829 2830 .keywords: SNES, nonlinear, set, convergence, tolerances 2831 2832 .seealso: SNESSetTrustRegionTolerance(), SNESSetLagPreconditioner() 2833 2834 @*/ 2835 PetscErrorCode SNESGetLagPreconditioner(SNES snes,PetscInt *lag) 2836 { 2837 PetscFunctionBegin; 2838 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2839 *lag = snes->lagpreconditioner; 2840 PetscFunctionReturn(0); 2841 } 2842 2843 #undef __FUNCT__ 2844 #define __FUNCT__ "SNESSetLagJacobian" 2845 /*@ 2846 SNESSetLagJacobian - Determines when the Jacobian is rebuilt in the nonlinear solve. See SNESSetLagPreconditioner() for determining how 2847 often the preconditioner is rebuilt. 2848 2849 Logically Collective on SNES 2850 2851 Input Parameters: 2852 + snes - the SNES context 2853 - lag - -1 indicates NEVER rebuild, 1 means rebuild every time the Jacobian is computed within a single nonlinear solve, 2 means every second time 2854 the Jacobian is built etc. -2 means rebuild at next chance but then never again 2855 2856 Options Database Keys: 2857 . -snes_lag_jacobian <lag> 2858 2859 Notes: 2860 The default is 1 2861 The Jacobian is ALWAYS built in the first iteration of a nonlinear solve unless lag is -1 2862 If -1 is used before the very first nonlinear solve the CODE WILL FAIL! because no Jacobian is used, use -2 to indicate you want it recomputed 2863 at the next Newton step but never again (unless it is reset to another value) 2864 2865 Level: intermediate 2866 2867 .keywords: SNES, nonlinear, set, convergence, tolerances 2868 2869 .seealso: SNESSetTrustRegionTolerance(), SNESGetLagPreconditioner(), SNESSetLagPreconditioner(), SNESGetLagJacobian() 2870 2871 @*/ 2872 PetscErrorCode SNESSetLagJacobian(SNES snes,PetscInt lag) 2873 { 2874 PetscFunctionBegin; 2875 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2876 if (lag < -2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Lag must be -2, -1, 1 or greater"); 2877 if (!lag) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Lag cannot be 0"); 2878 PetscValidLogicalCollectiveInt(snes,lag,2); 2879 snes->lagjacobian = lag; 2880 PetscFunctionReturn(0); 2881 } 2882 2883 #undef __FUNCT__ 2884 #define __FUNCT__ "SNESGetLagJacobian" 2885 /*@ 2886 SNESGetLagJacobian - Indicates how often the Jacobian is rebuilt. See SNESGetLagPreconditioner() to determine when the preconditioner is rebuilt 2887 2888 Not Collective 2889 2890 Input Parameter: 2891 . snes - the SNES context 2892 2893 Output Parameter: 2894 . lag - -1 indicates NEVER rebuild, 1 means rebuild every time the Jacobian is computed within a single nonlinear solve, 2 means every second time 2895 the Jacobian is built etc. 2896 2897 Options Database Keys: 2898 . -snes_lag_jacobian <lag> 2899 2900 Notes: 2901 The default is 1 2902 The jacobian is ALWAYS built in the first iteration of a nonlinear solve unless lag is -1 2903 2904 Level: intermediate 2905 2906 .keywords: SNES, nonlinear, set, convergence, tolerances 2907 2908 .seealso: SNESSetTrustRegionTolerance(), SNESSetLagJacobian(), SNESSetLagPreconditioner(), SNESGetLagPreconditioner() 2909 2910 @*/ 2911 PetscErrorCode SNESGetLagJacobian(SNES snes,PetscInt *lag) 2912 { 2913 PetscFunctionBegin; 2914 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2915 *lag = snes->lagjacobian; 2916 PetscFunctionReturn(0); 2917 } 2918 2919 #undef __FUNCT__ 2920 #define __FUNCT__ "SNESSetLagJacobianPersists" 2921 /*@ 2922 SNESSetLagJacobianPersists - Set whether or not the Jacobian lagging persists through multiple solves 2923 2924 Logically collective on SNES 2925 2926 Input Parameter: 2927 + snes - the SNES context 2928 - flg - jacobian lagging persists if true 2929 2930 Options Database Keys: 2931 . -snes_lag_jacobian_persists <flg> 2932 2933 Notes: This is useful both for nonlinear preconditioning, where it's appropriate to have the Jacobian be stale by 2934 several solves, and for implicit time-stepping, where Jacobian lagging in the inner nonlinear solve over several 2935 timesteps may present huge efficiency gains. 2936 2937 Level: developer 2938 2939 .keywords: SNES, nonlinear, lag 2940 2941 .seealso: SNESSetLagPreconditionerPersists(), SNESSetLagJacobian(), SNESGetLagJacobian(), SNESGetNPC() 2942 2943 @*/ 2944 PetscErrorCode SNESSetLagJacobianPersists(SNES snes,PetscBool flg) 2945 { 2946 PetscFunctionBegin; 2947 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2948 PetscValidLogicalCollectiveBool(snes,flg,2); 2949 snes->lagjac_persist = flg; 2950 PetscFunctionReturn(0); 2951 } 2952 2953 #undef __FUNCT__ 2954 #define __FUNCT__ "SNESSetLagPreconditionerPersists" 2955 /*@ 2956 SNESSetLagPreconditionerPersists - Set whether or not the preconditioner lagging persists through multiple solves 2957 2958 Logically Collective on SNES 2959 2960 Input Parameter: 2961 + snes - the SNES context 2962 - flg - preconditioner lagging persists if true 2963 2964 Options Database Keys: 2965 . -snes_lag_jacobian_persists <flg> 2966 2967 Notes: This is useful both for nonlinear preconditioning, where it's appropriate to have the preconditioner be stale 2968 by several solves, and for implicit time-stepping, where preconditioner lagging in the inner nonlinear solve over 2969 several timesteps may present huge efficiency gains. 2970 2971 Level: developer 2972 2973 .keywords: SNES, nonlinear, lag 2974 2975 .seealso: SNESSetLagJacobianPersists(), SNESSetLagJacobian(), SNESGetLagJacobian(), SNESGetNPC() 2976 2977 @*/ 2978 PetscErrorCode SNESSetLagPreconditionerPersists(SNES snes,PetscBool flg) 2979 { 2980 PetscFunctionBegin; 2981 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2982 PetscValidLogicalCollectiveBool(snes,flg,2); 2983 snes->lagpre_persist = flg; 2984 PetscFunctionReturn(0); 2985 } 2986 2987 #undef __FUNCT__ 2988 #define __FUNCT__ "SNESSetTolerances" 2989 /*@ 2990 SNESSetTolerances - Sets various parameters used in convergence tests. 2991 2992 Logically Collective on SNES 2993 2994 Input Parameters: 2995 + snes - the SNES context 2996 . abstol - absolute convergence tolerance 2997 . rtol - relative convergence tolerance 2998 . stol - convergence tolerance in terms of the norm of the change in the solution between steps, || delta x || < stol*|| x || 2999 . maxit - maximum number of iterations 3000 - maxf - maximum number of function evaluations 3001 3002 Options Database Keys: 3003 + -snes_atol <abstol> - Sets abstol 3004 . -snes_rtol <rtol> - Sets rtol 3005 . -snes_stol <stol> - Sets stol 3006 . -snes_max_it <maxit> - Sets maxit 3007 - -snes_max_funcs <maxf> - Sets maxf 3008 3009 Notes: 3010 The default maximum number of iterations is 50. 3011 The default maximum number of function evaluations is 1000. 3012 3013 Level: intermediate 3014 3015 .keywords: SNES, nonlinear, set, convergence, tolerances 3016 3017 .seealso: SNESSetTrustRegionTolerance() 3018 @*/ 3019 PetscErrorCode SNESSetTolerances(SNES snes,PetscReal abstol,PetscReal rtol,PetscReal stol,PetscInt maxit,PetscInt maxf) 3020 { 3021 PetscFunctionBegin; 3022 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3023 PetscValidLogicalCollectiveReal(snes,abstol,2); 3024 PetscValidLogicalCollectiveReal(snes,rtol,3); 3025 PetscValidLogicalCollectiveReal(snes,stol,4); 3026 PetscValidLogicalCollectiveInt(snes,maxit,5); 3027 PetscValidLogicalCollectiveInt(snes,maxf,6); 3028 3029 if (abstol != PETSC_DEFAULT) { 3030 if (abstol < 0.0) SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_OUTOFRANGE,"Absolute tolerance %g must be non-negative",(double)abstol); 3031 snes->abstol = abstol; 3032 } 3033 if (rtol != PETSC_DEFAULT) { 3034 if (rtol < 0.0 || 1.0 <= rtol) SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_OUTOFRANGE,"Relative tolerance %g must be non-negative and less than 1.0",(double)rtol); 3035 snes->rtol = rtol; 3036 } 3037 if (stol != PETSC_DEFAULT) { 3038 if (stol < 0.0) SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_OUTOFRANGE,"Step tolerance %g must be non-negative",(double)stol); 3039 snes->stol = stol; 3040 } 3041 if (maxit != PETSC_DEFAULT) { 3042 if (maxit < 0) SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_OUTOFRANGE,"Maximum number of iterations %D must be non-negative",maxit); 3043 snes->max_its = maxit; 3044 } 3045 if (maxf != PETSC_DEFAULT) { 3046 if (maxf < 0) SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_OUTOFRANGE,"Maximum number of function evaluations %D must be non-negative",maxf); 3047 snes->max_funcs = maxf; 3048 } 3049 snes->tolerancesset = PETSC_TRUE; 3050 PetscFunctionReturn(0); 3051 } 3052 3053 #undef __FUNCT__ 3054 #define __FUNCT__ "SNESGetTolerances" 3055 /*@ 3056 SNESGetTolerances - Gets various parameters used in convergence tests. 3057 3058 Not Collective 3059 3060 Input Parameters: 3061 + snes - the SNES context 3062 . atol - absolute convergence tolerance 3063 . rtol - relative convergence tolerance 3064 . stol - convergence tolerance in terms of the norm 3065 of the change in the solution between steps 3066 . maxit - maximum number of iterations 3067 - maxf - maximum number of function evaluations 3068 3069 Notes: 3070 The user can specify NULL for any parameter that is not needed. 3071 3072 Level: intermediate 3073 3074 .keywords: SNES, nonlinear, get, convergence, tolerances 3075 3076 .seealso: SNESSetTolerances() 3077 @*/ 3078 PetscErrorCode SNESGetTolerances(SNES snes,PetscReal *atol,PetscReal *rtol,PetscReal *stol,PetscInt *maxit,PetscInt *maxf) 3079 { 3080 PetscFunctionBegin; 3081 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3082 if (atol) *atol = snes->abstol; 3083 if (rtol) *rtol = snes->rtol; 3084 if (stol) *stol = snes->stol; 3085 if (maxit) *maxit = snes->max_its; 3086 if (maxf) *maxf = snes->max_funcs; 3087 PetscFunctionReturn(0); 3088 } 3089 3090 #undef __FUNCT__ 3091 #define __FUNCT__ "SNESSetTrustRegionTolerance" 3092 /*@ 3093 SNESSetTrustRegionTolerance - Sets the trust region parameter tolerance. 3094 3095 Logically Collective on SNES 3096 3097 Input Parameters: 3098 + snes - the SNES context 3099 - tol - tolerance 3100 3101 Options Database Key: 3102 . -snes_trtol <tol> - Sets tol 3103 3104 Level: intermediate 3105 3106 .keywords: SNES, nonlinear, set, trust region, tolerance 3107 3108 .seealso: SNESSetTolerances() 3109 @*/ 3110 PetscErrorCode SNESSetTrustRegionTolerance(SNES snes,PetscReal tol) 3111 { 3112 PetscFunctionBegin; 3113 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3114 PetscValidLogicalCollectiveReal(snes,tol,2); 3115 snes->deltatol = tol; 3116 PetscFunctionReturn(0); 3117 } 3118 3119 /* 3120 Duplicate the lg monitors for SNES from KSP; for some reason with 3121 dynamic libraries things don't work under Sun4 if we just use 3122 macros instead of functions 3123 */ 3124 #undef __FUNCT__ 3125 #define __FUNCT__ "SNESMonitorLGResidualNorm" 3126 PetscErrorCode SNESMonitorLGResidualNorm(SNES snes,PetscInt it,PetscReal norm,void *ctx) 3127 { 3128 PetscErrorCode ierr; 3129 3130 PetscFunctionBegin; 3131 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3132 ierr = KSPMonitorLGResidualNorm((KSP)snes,it,norm,ctx);CHKERRQ(ierr); 3133 PetscFunctionReturn(0); 3134 } 3135 3136 #undef __FUNCT__ 3137 #define __FUNCT__ "SNESMonitorLGCreate" 3138 PetscErrorCode SNESMonitorLGCreate(const char host[],const char label[],int x,int y,int m,int n,PetscDrawLG *draw) 3139 { 3140 PetscErrorCode ierr; 3141 3142 PetscFunctionBegin; 3143 ierr = KSPMonitorLGResidualNormCreate(host,label,x,y,m,n,draw);CHKERRQ(ierr); 3144 PetscFunctionReturn(0); 3145 } 3146 3147 #undef __FUNCT__ 3148 #define __FUNCT__ "SNESMonitorLGDestroy" 3149 PetscErrorCode SNESMonitorLGDestroy(PetscDrawLG *draw) 3150 { 3151 PetscErrorCode ierr; 3152 3153 PetscFunctionBegin; 3154 ierr = KSPMonitorLGResidualNormDestroy(draw);CHKERRQ(ierr); 3155 PetscFunctionReturn(0); 3156 } 3157 3158 extern PetscErrorCode SNESMonitorRange_Private(SNES,PetscInt,PetscReal*); 3159 #undef __FUNCT__ 3160 #define __FUNCT__ "SNESMonitorLGRange" 3161 PetscErrorCode SNESMonitorLGRange(SNES snes,PetscInt n,PetscReal rnorm,void *monctx) 3162 { 3163 PetscDrawLG lg; 3164 PetscErrorCode ierr; 3165 PetscReal x,y,per; 3166 PetscViewer v = (PetscViewer)monctx; 3167 static PetscReal prev; /* should be in the context */ 3168 PetscDraw draw; 3169 3170 PetscFunctionBegin; 3171 ierr = PetscViewerDrawGetDrawLG(v,0,&lg);CHKERRQ(ierr); 3172 if (!n) {ierr = PetscDrawLGReset(lg);CHKERRQ(ierr);} 3173 ierr = PetscDrawLGGetDraw(lg,&draw);CHKERRQ(ierr); 3174 ierr = PetscDrawSetTitle(draw,"Residual norm");CHKERRQ(ierr); 3175 x = (PetscReal)n; 3176 if (rnorm > 0.0) y = PetscLog10Real(rnorm); 3177 else y = -15.0; 3178 ierr = PetscDrawLGAddPoint(lg,&x,&y);CHKERRQ(ierr); 3179 if (n < 20 || !(n % 5)) { 3180 ierr = PetscDrawLGDraw(lg);CHKERRQ(ierr); 3181 } 3182 3183 ierr = PetscViewerDrawGetDrawLG(v,1,&lg);CHKERRQ(ierr); 3184 if (!n) {ierr = PetscDrawLGReset(lg);CHKERRQ(ierr);} 3185 ierr = PetscDrawLGGetDraw(lg,&draw);CHKERRQ(ierr); 3186 ierr = PetscDrawSetTitle(draw,"% elemts > .2*max elemt");CHKERRQ(ierr); 3187 ierr = SNESMonitorRange_Private(snes,n,&per);CHKERRQ(ierr); 3188 x = (PetscReal)n; 3189 y = 100.0*per; 3190 ierr = PetscDrawLGAddPoint(lg,&x,&y);CHKERRQ(ierr); 3191 if (n < 20 || !(n % 5)) { 3192 ierr = PetscDrawLGDraw(lg);CHKERRQ(ierr); 3193 } 3194 3195 ierr = PetscViewerDrawGetDrawLG(v,2,&lg);CHKERRQ(ierr); 3196 if (!n) {prev = rnorm;ierr = PetscDrawLGReset(lg);CHKERRQ(ierr);} 3197 ierr = PetscDrawLGGetDraw(lg,&draw);CHKERRQ(ierr); 3198 ierr = PetscDrawSetTitle(draw,"(norm -oldnorm)/oldnorm");CHKERRQ(ierr); 3199 x = (PetscReal)n; 3200 y = (prev - rnorm)/prev; 3201 ierr = PetscDrawLGAddPoint(lg,&x,&y);CHKERRQ(ierr); 3202 if (n < 20 || !(n % 5)) { 3203 ierr = PetscDrawLGDraw(lg);CHKERRQ(ierr); 3204 } 3205 3206 ierr = PetscViewerDrawGetDrawLG(v,3,&lg);CHKERRQ(ierr); 3207 if (!n) {ierr = PetscDrawLGReset(lg);CHKERRQ(ierr);} 3208 ierr = PetscDrawLGGetDraw(lg,&draw);CHKERRQ(ierr); 3209 ierr = PetscDrawSetTitle(draw,"(norm -oldnorm)/oldnorm*(% > .2 max)");CHKERRQ(ierr); 3210 x = (PetscReal)n; 3211 y = (prev - rnorm)/(prev*per); 3212 if (n > 2) { /*skip initial crazy value */ 3213 ierr = PetscDrawLGAddPoint(lg,&x,&y);CHKERRQ(ierr); 3214 } 3215 if (n < 20 || !(n % 5)) { 3216 ierr = PetscDrawLGDraw(lg);CHKERRQ(ierr); 3217 } 3218 prev = rnorm; 3219 PetscFunctionReturn(0); 3220 } 3221 3222 #undef __FUNCT__ 3223 #define __FUNCT__ "SNESMonitor" 3224 /*@ 3225 SNESMonitor - runs the user provided monitor routines, if they exist 3226 3227 Collective on SNES 3228 3229 Input Parameters: 3230 + snes - nonlinear solver context obtained from SNESCreate() 3231 . iter - iteration number 3232 - rnorm - relative norm of the residual 3233 3234 Notes: 3235 This routine is called by the SNES implementations. 3236 It does not typically need to be called by the user. 3237 3238 Level: developer 3239 3240 .seealso: SNESMonitorSet() 3241 @*/ 3242 PetscErrorCode SNESMonitor(SNES snes,PetscInt iter,PetscReal rnorm) 3243 { 3244 PetscErrorCode ierr; 3245 PetscInt i,n = snes->numbermonitors; 3246 3247 PetscFunctionBegin; 3248 for (i=0; i<n; i++) { 3249 ierr = (*snes->monitor[i])(snes,iter,rnorm,snes->monitorcontext[i]);CHKERRQ(ierr); 3250 } 3251 PetscFunctionReturn(0); 3252 } 3253 3254 /* ------------ Routines to set performance monitoring options ----------- */ 3255 3256 /*MC 3257 SNESMonitorFunction - functional form passed to SNESMonitorSet() to monitor convergence of nonlinear solver 3258 3259 Synopsis: 3260 #include <petscsnes.h> 3261 $ PetscErrorCode SNESMonitorFunction(SNES snes,PetscInt its, PetscReal norm,void *mctx) 3262 3263 + snes - the SNES context 3264 . its - iteration number 3265 . norm - 2-norm function value (may be estimated) 3266 - mctx - [optional] monitoring context 3267 3268 Level: advanced 3269 3270 .seealso: SNESMonitorSet(), SNESMonitorGet() 3271 M*/ 3272 3273 #undef __FUNCT__ 3274 #define __FUNCT__ "SNESMonitorSet" 3275 /*@C 3276 SNESMonitorSet - Sets an ADDITIONAL function that is to be used at every 3277 iteration of the nonlinear solver to display the iteration's 3278 progress. 3279 3280 Logically Collective on SNES 3281 3282 Input Parameters: 3283 + snes - the SNES context 3284 . f - the monitor function, see SNESMonitorFunction for the calling sequence 3285 . mctx - [optional] user-defined context for private data for the 3286 monitor routine (use NULL if no context is desired) 3287 - monitordestroy - [optional] routine that frees monitor context 3288 (may be NULL) 3289 3290 Options Database Keys: 3291 + -snes_monitor - sets SNESMonitorDefault() 3292 . -snes_monitor_lg_residualnorm - sets line graph monitor, 3293 uses SNESMonitorLGCreate() 3294 - -snes_monitor_cancel - cancels all monitors that have 3295 been hardwired into a code by 3296 calls to SNESMonitorSet(), but 3297 does not cancel those set via 3298 the options database. 3299 3300 Notes: 3301 Several different monitoring routines may be set by calling 3302 SNESMonitorSet() multiple times; all will be called in the 3303 order in which they were set. 3304 3305 Fortran notes: Only a single monitor function can be set for each SNES object 3306 3307 Level: intermediate 3308 3309 .keywords: SNES, nonlinear, set, monitor 3310 3311 .seealso: SNESMonitorDefault(), SNESMonitorCancel(), SNESMonitorFunction 3312 @*/ 3313 PetscErrorCode SNESMonitorSet(SNES snes,PetscErrorCode (*f)(SNES,PetscInt,PetscReal,void*),void *mctx,PetscErrorCode (*monitordestroy)(void**)) 3314 { 3315 PetscInt i; 3316 PetscErrorCode ierr; 3317 3318 PetscFunctionBegin; 3319 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3320 if (snes->numbermonitors >= MAXSNESMONITORS) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Too many monitors set"); 3321 for (i=0; i<snes->numbermonitors;i++) { 3322 if (f == snes->monitor[i] && monitordestroy == snes->monitordestroy[i] && mctx == snes->monitorcontext[i]) { 3323 if (monitordestroy) { 3324 ierr = (*monitordestroy)(&mctx);CHKERRQ(ierr); 3325 } 3326 PetscFunctionReturn(0); 3327 } 3328 } 3329 snes->monitor[snes->numbermonitors] = f; 3330 snes->monitordestroy[snes->numbermonitors] = monitordestroy; 3331 snes->monitorcontext[snes->numbermonitors++] = (void*)mctx; 3332 PetscFunctionReturn(0); 3333 } 3334 3335 #undef __FUNCT__ 3336 #define __FUNCT__ "SNESMonitorCancel" 3337 /*@ 3338 SNESMonitorCancel - Clears all the monitor functions for a SNES object. 3339 3340 Logically Collective on SNES 3341 3342 Input Parameters: 3343 . snes - the SNES context 3344 3345 Options Database Key: 3346 . -snes_monitor_cancel - cancels all monitors that have been hardwired 3347 into a code by calls to SNESMonitorSet(), but does not cancel those 3348 set via the options database 3349 3350 Notes: 3351 There is no way to clear one specific monitor from a SNES object. 3352 3353 Level: intermediate 3354 3355 .keywords: SNES, nonlinear, set, monitor 3356 3357 .seealso: SNESMonitorDefault(), SNESMonitorSet() 3358 @*/ 3359 PetscErrorCode SNESMonitorCancel(SNES snes) 3360 { 3361 PetscErrorCode ierr; 3362 PetscInt i; 3363 3364 PetscFunctionBegin; 3365 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3366 for (i=0; i<snes->numbermonitors; i++) { 3367 if (snes->monitordestroy[i]) { 3368 ierr = (*snes->monitordestroy[i])(&snes->monitorcontext[i]);CHKERRQ(ierr); 3369 } 3370 } 3371 snes->numbermonitors = 0; 3372 PetscFunctionReturn(0); 3373 } 3374 3375 /*MC 3376 SNESConvergenceTestFunction - functional form used for testing of convergence of nonlinear solver 3377 3378 Synopsis: 3379 #include <petscsnes.h> 3380 $ PetscErrorCode SNESConvergenceTest(SNES snes,PetscInt it,PetscReal xnorm,PetscReal gnorm,PetscReal f,SNESConvergedReason *reason,void *cctx) 3381 3382 + snes - the SNES context 3383 . it - current iteration (0 is the first and is before any Newton step) 3384 . cctx - [optional] convergence context 3385 . reason - reason for convergence/divergence 3386 . xnorm - 2-norm of current iterate 3387 . gnorm - 2-norm of current step 3388 - f - 2-norm of function 3389 3390 Level: intermediate 3391 3392 .seealso: SNESSetConvergenceTest(), SNESGetConvergenceTest() 3393 M*/ 3394 3395 #undef __FUNCT__ 3396 #define __FUNCT__ "SNESSetConvergenceTest" 3397 /*@C 3398 SNESSetConvergenceTest - Sets the function that is to be used 3399 to test for convergence of the nonlinear iterative solution. 3400 3401 Logically Collective on SNES 3402 3403 Input Parameters: 3404 + snes - the SNES context 3405 . SNESConvergenceTestFunction - routine to test for convergence 3406 . cctx - [optional] context for private data for the convergence routine (may be NULL) 3407 - destroy - [optional] destructor for the context (may be NULL; NULL_FUNCTION in Fortran) 3408 3409 Level: advanced 3410 3411 .keywords: SNES, nonlinear, set, convergence, test 3412 3413 .seealso: SNESConvergedDefault(), SNESConvergedSkip(), SNESConvergenceTestFunction 3414 @*/ 3415 PetscErrorCode SNESSetConvergenceTest(SNES snes,PetscErrorCode (*SNESConvergenceTestFunction)(SNES,PetscInt,PetscReal,PetscReal,PetscReal,SNESConvergedReason*,void*),void *cctx,PetscErrorCode (*destroy)(void*)) 3416 { 3417 PetscErrorCode ierr; 3418 3419 PetscFunctionBegin; 3420 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3421 if (!SNESConvergenceTestFunction) SNESConvergenceTestFunction = SNESConvergedSkip; 3422 if (snes->ops->convergeddestroy) { 3423 ierr = (*snes->ops->convergeddestroy)(snes->cnvP);CHKERRQ(ierr); 3424 } 3425 snes->ops->converged = SNESConvergenceTestFunction; 3426 snes->ops->convergeddestroy = destroy; 3427 snes->cnvP = cctx; 3428 PetscFunctionReturn(0); 3429 } 3430 3431 #undef __FUNCT__ 3432 #define __FUNCT__ "SNESGetConvergedReason" 3433 /*@ 3434 SNESGetConvergedReason - Gets the reason the SNES iteration was stopped. 3435 3436 Not Collective 3437 3438 Input Parameter: 3439 . snes - the SNES context 3440 3441 Output Parameter: 3442 . reason - negative value indicates diverged, positive value converged, see SNESConvergedReason or the 3443 manual pages for the individual convergence tests for complete lists 3444 3445 Level: intermediate 3446 3447 Notes: Can only be called after the call the SNESSolve() is complete. 3448 3449 .keywords: SNES, nonlinear, set, convergence, test 3450 3451 .seealso: SNESSetConvergenceTest(), SNESConvergedReason 3452 @*/ 3453 PetscErrorCode SNESGetConvergedReason(SNES snes,SNESConvergedReason *reason) 3454 { 3455 PetscFunctionBegin; 3456 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3457 PetscValidPointer(reason,2); 3458 *reason = snes->reason; 3459 PetscFunctionReturn(0); 3460 } 3461 3462 #undef __FUNCT__ 3463 #define __FUNCT__ "SNESSetConvergenceHistory" 3464 /*@ 3465 SNESSetConvergenceHistory - Sets the array used to hold the convergence history. 3466 3467 Logically Collective on SNES 3468 3469 Input Parameters: 3470 + snes - iterative context obtained from SNESCreate() 3471 . a - array to hold history, this array will contain the function norms computed at each step 3472 . its - integer array holds the number of linear iterations for each solve. 3473 . na - size of a and its 3474 - reset - PETSC_TRUE indicates each new nonlinear solve resets the history counter to zero, 3475 else it continues storing new values for new nonlinear solves after the old ones 3476 3477 Notes: 3478 If 'a' and 'its' are NULL then space is allocated for the history. If 'na' PETSC_DECIDE or PETSC_DEFAULT then a 3479 default array of length 10000 is allocated. 3480 3481 This routine is useful, e.g., when running a code for purposes 3482 of accurate performance monitoring, when no I/O should be done 3483 during the section of code that is being timed. 3484 3485 Level: intermediate 3486 3487 .keywords: SNES, set, convergence, history 3488 3489 .seealso: SNESGetConvergenceHistory() 3490 3491 @*/ 3492 PetscErrorCode SNESSetConvergenceHistory(SNES snes,PetscReal a[],PetscInt its[],PetscInt na,PetscBool reset) 3493 { 3494 PetscErrorCode ierr; 3495 3496 PetscFunctionBegin; 3497 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3498 if (a) PetscValidScalarPointer(a,2); 3499 if (its) PetscValidIntPointer(its,3); 3500 if (!a) { 3501 if (na == PETSC_DECIDE || na == PETSC_DEFAULT) na = 1000; 3502 ierr = PetscCalloc1(na,&a);CHKERRQ(ierr); 3503 ierr = PetscCalloc1(na,&its);CHKERRQ(ierr); 3504 3505 snes->conv_malloc = PETSC_TRUE; 3506 } 3507 snes->conv_hist = a; 3508 snes->conv_hist_its = its; 3509 snes->conv_hist_max = na; 3510 snes->conv_hist_len = 0; 3511 snes->conv_hist_reset = reset; 3512 PetscFunctionReturn(0); 3513 } 3514 3515 #if defined(PETSC_HAVE_MATLAB_ENGINE) 3516 #include <engine.h> /* MATLAB include file */ 3517 #include <mex.h> /* MATLAB include file */ 3518 3519 #undef __FUNCT__ 3520 #define __FUNCT__ "SNESGetConvergenceHistoryMatlab" 3521 PETSC_EXTERN mxArray *SNESGetConvergenceHistoryMatlab(SNES snes) 3522 { 3523 mxArray *mat; 3524 PetscInt i; 3525 PetscReal *ar; 3526 3527 PetscFunctionBegin; 3528 mat = mxCreateDoubleMatrix(snes->conv_hist_len,1,mxREAL); 3529 ar = (PetscReal*) mxGetData(mat); 3530 for (i=0; i<snes->conv_hist_len; i++) ar[i] = snes->conv_hist[i]; 3531 PetscFunctionReturn(mat); 3532 } 3533 #endif 3534 3535 #undef __FUNCT__ 3536 #define __FUNCT__ "SNESGetConvergenceHistory" 3537 /*@C 3538 SNESGetConvergenceHistory - Gets the array used to hold the convergence history. 3539 3540 Not Collective 3541 3542 Input Parameter: 3543 . snes - iterative context obtained from SNESCreate() 3544 3545 Output Parameters: 3546 . a - array to hold history 3547 . its - integer array holds the number of linear iterations (or 3548 negative if not converged) for each solve. 3549 - na - size of a and its 3550 3551 Notes: 3552 The calling sequence for this routine in Fortran is 3553 $ call SNESGetConvergenceHistory(SNES snes, integer na, integer ierr) 3554 3555 This routine is useful, e.g., when running a code for purposes 3556 of accurate performance monitoring, when no I/O should be done 3557 during the section of code that is being timed. 3558 3559 Level: intermediate 3560 3561 .keywords: SNES, get, convergence, history 3562 3563 .seealso: SNESSetConvergencHistory() 3564 3565 @*/ 3566 PetscErrorCode SNESGetConvergenceHistory(SNES snes,PetscReal *a[],PetscInt *its[],PetscInt *na) 3567 { 3568 PetscFunctionBegin; 3569 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3570 if (a) *a = snes->conv_hist; 3571 if (its) *its = snes->conv_hist_its; 3572 if (na) *na = snes->conv_hist_len; 3573 PetscFunctionReturn(0); 3574 } 3575 3576 #undef __FUNCT__ 3577 #define __FUNCT__ "SNESSetUpdate" 3578 /*@C 3579 SNESSetUpdate - Sets the general-purpose update function called 3580 at the beginning of every iteration of the nonlinear solve. Specifically 3581 it is called just before the Jacobian is "evaluated". 3582 3583 Logically Collective on SNES 3584 3585 Input Parameters: 3586 . snes - The nonlinear solver context 3587 . func - The function 3588 3589 Calling sequence of func: 3590 . func (SNES snes, PetscInt step); 3591 3592 . step - The current step of the iteration 3593 3594 Level: advanced 3595 3596 Note: This is NOT what one uses to update the ghost points before a function evaluation, that should be done at the beginning of your FormFunction() 3597 This is not used by most users. 3598 3599 .keywords: SNES, update 3600 3601 .seealso SNESSetJacobian(), SNESSolve() 3602 @*/ 3603 PetscErrorCode SNESSetUpdate(SNES snes, PetscErrorCode (*func)(SNES, PetscInt)) 3604 { 3605 PetscFunctionBegin; 3606 PetscValidHeaderSpecific(snes, SNES_CLASSID,1); 3607 snes->ops->update = func; 3608 PetscFunctionReturn(0); 3609 } 3610 3611 #undef __FUNCT__ 3612 #define __FUNCT__ "SNESScaleStep_Private" 3613 /* 3614 SNESScaleStep_Private - Scales a step so that its length is less than the 3615 positive parameter delta. 3616 3617 Input Parameters: 3618 + snes - the SNES context 3619 . y - approximate solution of linear system 3620 . fnorm - 2-norm of current function 3621 - delta - trust region size 3622 3623 Output Parameters: 3624 + gpnorm - predicted function norm at the new point, assuming local 3625 linearization. The value is zero if the step lies within the trust 3626 region, and exceeds zero otherwise. 3627 - ynorm - 2-norm of the step 3628 3629 Note: 3630 For non-trust region methods such as SNESNEWTONLS, the parameter delta 3631 is set to be the maximum allowable step size. 3632 3633 .keywords: SNES, nonlinear, scale, step 3634 */ 3635 PetscErrorCode SNESScaleStep_Private(SNES snes,Vec y,PetscReal *fnorm,PetscReal *delta,PetscReal *gpnorm,PetscReal *ynorm) 3636 { 3637 PetscReal nrm; 3638 PetscScalar cnorm; 3639 PetscErrorCode ierr; 3640 3641 PetscFunctionBegin; 3642 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3643 PetscValidHeaderSpecific(y,VEC_CLASSID,2); 3644 PetscCheckSameComm(snes,1,y,2); 3645 3646 ierr = VecNorm(y,NORM_2,&nrm);CHKERRQ(ierr); 3647 if (nrm > *delta) { 3648 nrm = *delta/nrm; 3649 *gpnorm = (1.0 - nrm)*(*fnorm); 3650 cnorm = nrm; 3651 ierr = VecScale(y,cnorm);CHKERRQ(ierr); 3652 *ynorm = *delta; 3653 } else { 3654 *gpnorm = 0.0; 3655 *ynorm = nrm; 3656 } 3657 PetscFunctionReturn(0); 3658 } 3659 3660 #undef __FUNCT__ 3661 #define __FUNCT__ "SNESSolve" 3662 /*@C 3663 SNESSolve - Solves a nonlinear system F(x) = b. 3664 Call SNESSolve() after calling SNESCreate() and optional routines of the form SNESSetXXX(). 3665 3666 Collective on SNES 3667 3668 Input Parameters: 3669 + snes - the SNES context 3670 . b - the constant part of the equation F(x) = b, or NULL to use zero. 3671 - x - the solution vector. 3672 3673 Notes: 3674 The user should initialize the vector,x, with the initial guess 3675 for the nonlinear solve prior to calling SNESSolve. In particular, 3676 to employ an initial guess of zero, the user should explicitly set 3677 this vector to zero by calling VecSet(). 3678 3679 Level: beginner 3680 3681 .keywords: SNES, nonlinear, solve 3682 3683 .seealso: SNESCreate(), SNESDestroy(), SNESSetFunction(), SNESSetJacobian(), SNESSetGridSequence(), SNESGetSolution() 3684 @*/ 3685 PetscErrorCode SNESSolve(SNES snes,Vec b,Vec x) 3686 { 3687 PetscErrorCode ierr; 3688 PetscBool flg; 3689 PetscInt grid; 3690 Vec xcreated = NULL; 3691 DM dm; 3692 3693 PetscFunctionBegin; 3694 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3695 if (x) PetscValidHeaderSpecific(x,VEC_CLASSID,3); 3696 if (x) PetscCheckSameComm(snes,1,x,3); 3697 if (b) PetscValidHeaderSpecific(b,VEC_CLASSID,2); 3698 if (b) PetscCheckSameComm(snes,1,b,2); 3699 3700 if (!x) { 3701 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 3702 ierr = DMCreateGlobalVector(dm,&xcreated);CHKERRQ(ierr); 3703 x = xcreated; 3704 } 3705 ierr = SNESViewFromOptions(snes,NULL,"-snes_view_pre");CHKERRQ(ierr); 3706 3707 for (grid=0; grid<snes->gridsequence; grid++) {ierr = PetscViewerASCIIPushTab(PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)snes)));CHKERRQ(ierr);} 3708 for (grid=0; grid<snes->gridsequence+1; grid++) { 3709 3710 /* set solution vector */ 3711 if (!grid) {ierr = PetscObjectReference((PetscObject)x);CHKERRQ(ierr);} 3712 ierr = VecDestroy(&snes->vec_sol);CHKERRQ(ierr); 3713 snes->vec_sol = x; 3714 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 3715 3716 /* set affine vector if provided */ 3717 if (b) { ierr = PetscObjectReference((PetscObject)b);CHKERRQ(ierr); } 3718 ierr = VecDestroy(&snes->vec_rhs);CHKERRQ(ierr); 3719 snes->vec_rhs = b; 3720 3721 if (snes->vec_func == snes->vec_sol) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_IDN,"Solution vector cannot be function vector"); 3722 if (snes->vec_rhs == snes->vec_sol) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_IDN,"Solution vector cannot be right hand side vector"); 3723 if (!snes->vec_sol_update /* && snes->vec_sol */) { 3724 ierr = VecDuplicate(snes->vec_sol,&snes->vec_sol_update);CHKERRQ(ierr); 3725 ierr = PetscLogObjectParent((PetscObject)snes,(PetscObject)snes->vec_sol_update);CHKERRQ(ierr); 3726 } 3727 ierr = DMShellSetGlobalVector(dm,snes->vec_sol);CHKERRQ(ierr); 3728 ierr = SNESSetUp(snes);CHKERRQ(ierr); 3729 3730 if (!grid) { 3731 if (snes->ops->computeinitialguess) { 3732 ierr = (*snes->ops->computeinitialguess)(snes,snes->vec_sol,snes->initialguessP);CHKERRQ(ierr); 3733 } 3734 } 3735 3736 if (snes->conv_hist_reset) snes->conv_hist_len = 0; 3737 if (snes->counters_reset) {snes->nfuncs = 0; snes->linear_its = 0; snes->numFailures = 0;} 3738 3739 ierr = PetscLogEventBegin(SNES_Solve,snes,0,0,0);CHKERRQ(ierr); 3740 ierr = (*snes->ops->solve)(snes);CHKERRQ(ierr); 3741 ierr = PetscLogEventEnd(SNES_Solve,snes,0,0,0);CHKERRQ(ierr); 3742 if (snes->domainerror) { 3743 snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; 3744 snes->domainerror = PETSC_FALSE; 3745 } 3746 if (!snes->reason) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Internal error, solver returned without setting converged reason"); 3747 3748 if (snes->lagjac_persist) snes->jac_iter += snes->iter; 3749 if (snes->lagpre_persist) snes->pre_iter += snes->iter; 3750 3751 flg = PETSC_FALSE; 3752 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_test_local_min",&flg,NULL);CHKERRQ(ierr); 3753 if (flg && !PetscPreLoadingOn) { ierr = SNESTestLocalMin(snes);CHKERRQ(ierr); } 3754 if (snes->printreason) { 3755 ierr = PetscViewerASCIIAddTab(PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)snes)),((PetscObject)snes)->tablevel);CHKERRQ(ierr); 3756 if (snes->reason > 0) { 3757 ierr = PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)snes)),"Nonlinear solve converged due to %s iterations %D\n",SNESConvergedReasons[snes->reason],snes->iter);CHKERRQ(ierr); 3758 } else { 3759 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); 3760 } 3761 ierr = PetscViewerASCIISubtractTab(PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)snes)),((PetscObject)snes)->tablevel);CHKERRQ(ierr); 3762 } 3763 3764 if (snes->errorifnotconverged && snes->reason < 0) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_NOT_CONVERGED,"SNESSolve has not converged"); 3765 if (grid < snes->gridsequence) { 3766 DM fine; 3767 Vec xnew; 3768 Mat interp; 3769 3770 ierr = DMRefine(snes->dm,PetscObjectComm((PetscObject)snes),&fine);CHKERRQ(ierr); 3771 if (!fine) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_INCOMP,"DMRefine() did not perform any refinement, cannot continue grid sequencing"); 3772 ierr = DMCreateInterpolation(snes->dm,fine,&interp,NULL);CHKERRQ(ierr); 3773 ierr = DMCreateGlobalVector(fine,&xnew);CHKERRQ(ierr); 3774 ierr = MatInterpolate(interp,x,xnew);CHKERRQ(ierr); 3775 ierr = DMInterpolate(snes->dm,interp,fine);CHKERRQ(ierr); 3776 ierr = MatDestroy(&interp);CHKERRQ(ierr); 3777 x = xnew; 3778 3779 ierr = SNESReset(snes);CHKERRQ(ierr); 3780 ierr = SNESSetDM(snes,fine);CHKERRQ(ierr); 3781 ierr = DMDestroy(&fine);CHKERRQ(ierr); 3782 ierr = PetscViewerASCIIPopTab(PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)snes)));CHKERRQ(ierr); 3783 } 3784 } 3785 ierr = SNESViewFromOptions(snes,NULL,"-snes_view");CHKERRQ(ierr); 3786 ierr = VecViewFromOptions(snes->vec_sol,((PetscObject)snes)->prefix,"-snes_view_solution");CHKERRQ(ierr); 3787 3788 ierr = VecDestroy(&xcreated);CHKERRQ(ierr); 3789 ierr = PetscObjectSAWsBlock((PetscObject)snes);CHKERRQ(ierr); 3790 PetscFunctionReturn(0); 3791 } 3792 3793 /* --------- Internal routines for SNES Package --------- */ 3794 3795 #undef __FUNCT__ 3796 #define __FUNCT__ "SNESSetType" 3797 /*@C 3798 SNESSetType - Sets the method for the nonlinear solver. 3799 3800 Collective on SNES 3801 3802 Input Parameters: 3803 + snes - the SNES context 3804 - type - a known method 3805 3806 Options Database Key: 3807 . -snes_type <type> - Sets the method; use -help for a list 3808 of available methods (for instance, newtonls or newtontr) 3809 3810 Notes: 3811 See "petsc/include/petscsnes.h" for available methods (for instance) 3812 + SNESNEWTONLS - Newton's method with line search 3813 (systems of nonlinear equations) 3814 . SNESNEWTONTR - Newton's method with trust region 3815 (systems of nonlinear equations) 3816 3817 Normally, it is best to use the SNESSetFromOptions() command and then 3818 set the SNES solver type from the options database rather than by using 3819 this routine. Using the options database provides the user with 3820 maximum flexibility in evaluating the many nonlinear solvers. 3821 The SNESSetType() routine is provided for those situations where it 3822 is necessary to set the nonlinear solver independently of the command 3823 line or options database. This might be the case, for example, when 3824 the choice of solver changes during the execution of the program, 3825 and the user's application is taking responsibility for choosing the 3826 appropriate method. 3827 3828 Developer Notes: SNESRegister() adds a constructor for a new SNESType to SNESList, SNESSetType() locates 3829 the constructor in that list and calls it to create the spexific object. 3830 3831 Level: intermediate 3832 3833 .keywords: SNES, set, type 3834 3835 .seealso: SNESType, SNESCreate(), SNESDestroy(), SNESGetType(), SNESSetFromOptions() 3836 3837 @*/ 3838 PetscErrorCode SNESSetType(SNES snes,SNESType type) 3839 { 3840 PetscErrorCode ierr,(*r)(SNES); 3841 PetscBool match; 3842 3843 PetscFunctionBegin; 3844 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3845 PetscValidCharPointer(type,2); 3846 3847 ierr = PetscObjectTypeCompare((PetscObject)snes,type,&match);CHKERRQ(ierr); 3848 if (match) PetscFunctionReturn(0); 3849 3850 ierr = PetscFunctionListFind(SNESList,type,&r);CHKERRQ(ierr); 3851 if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unable to find requested SNES type %s",type); 3852 /* Destroy the previous private SNES context */ 3853 if (snes->ops->destroy) { 3854 ierr = (*(snes)->ops->destroy)(snes);CHKERRQ(ierr); 3855 snes->ops->destroy = NULL; 3856 } 3857 /* Reinitialize function pointers in SNESOps structure */ 3858 snes->ops->setup = 0; 3859 snes->ops->solve = 0; 3860 snes->ops->view = 0; 3861 snes->ops->setfromoptions = 0; 3862 snes->ops->destroy = 0; 3863 /* Call the SNESCreate_XXX routine for this particular Nonlinear solver */ 3864 snes->setupcalled = PETSC_FALSE; 3865 3866 ierr = PetscObjectChangeTypeName((PetscObject)snes,type);CHKERRQ(ierr); 3867 ierr = (*r)(snes);CHKERRQ(ierr); 3868 PetscFunctionReturn(0); 3869 } 3870 3871 #undef __FUNCT__ 3872 #define __FUNCT__ "SNESGetType" 3873 /*@C 3874 SNESGetType - Gets the SNES method type and name (as a string). 3875 3876 Not Collective 3877 3878 Input Parameter: 3879 . snes - nonlinear solver context 3880 3881 Output Parameter: 3882 . type - SNES method (a character string) 3883 3884 Level: intermediate 3885 3886 .keywords: SNES, nonlinear, get, type, name 3887 @*/ 3888 PetscErrorCode SNESGetType(SNES snes,SNESType *type) 3889 { 3890 PetscFunctionBegin; 3891 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3892 PetscValidPointer(type,2); 3893 *type = ((PetscObject)snes)->type_name; 3894 PetscFunctionReturn(0); 3895 } 3896 3897 #undef __FUNCT__ 3898 #define __FUNCT__ "SNESGetSolution" 3899 /*@ 3900 SNESGetSolution - Returns the vector where the approximate solution is 3901 stored. This is the fine grid solution when using SNESSetGridSequence(). 3902 3903 Not Collective, but Vec is parallel if SNES is parallel 3904 3905 Input Parameter: 3906 . snes - the SNES context 3907 3908 Output Parameter: 3909 . x - the solution 3910 3911 Level: intermediate 3912 3913 .keywords: SNES, nonlinear, get, solution 3914 3915 .seealso: SNESGetSolutionUpdate(), SNESGetFunction() 3916 @*/ 3917 PetscErrorCode SNESGetSolution(SNES snes,Vec *x) 3918 { 3919 PetscFunctionBegin; 3920 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3921 PetscValidPointer(x,2); 3922 *x = snes->vec_sol; 3923 PetscFunctionReturn(0); 3924 } 3925 3926 #undef __FUNCT__ 3927 #define __FUNCT__ "SNESGetSolutionUpdate" 3928 /*@ 3929 SNESGetSolutionUpdate - Returns the vector where the solution update is 3930 stored. 3931 3932 Not Collective, but Vec is parallel if SNES is parallel 3933 3934 Input Parameter: 3935 . snes - the SNES context 3936 3937 Output Parameter: 3938 . x - the solution update 3939 3940 Level: advanced 3941 3942 .keywords: SNES, nonlinear, get, solution, update 3943 3944 .seealso: SNESGetSolution(), SNESGetFunction() 3945 @*/ 3946 PetscErrorCode SNESGetSolutionUpdate(SNES snes,Vec *x) 3947 { 3948 PetscFunctionBegin; 3949 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3950 PetscValidPointer(x,2); 3951 *x = snes->vec_sol_update; 3952 PetscFunctionReturn(0); 3953 } 3954 3955 #undef __FUNCT__ 3956 #define __FUNCT__ "SNESGetFunction" 3957 /*@C 3958 SNESGetFunction - Returns the vector where the function is stored. 3959 3960 Not Collective, but Vec is parallel if SNES is parallel. Collective if Vec is requested, but has not been created yet. 3961 3962 Input Parameter: 3963 . snes - the SNES context 3964 3965 Output Parameter: 3966 + r - the vector that is used to store residuals (or NULL if you don't want it) 3967 . f - the function (or NULL if you don't want it); see SNESFunction for calling sequence details 3968 - ctx - the function context (or NULL if you don't want it) 3969 3970 Level: advanced 3971 3972 .keywords: SNES, nonlinear, get, function 3973 3974 .seealso: SNESSetFunction(), SNESGetSolution(), SNESFunction 3975 @*/ 3976 PetscErrorCode SNESGetFunction(SNES snes,Vec *r,PetscErrorCode (**f)(SNES,Vec,Vec,void*),void **ctx) 3977 { 3978 PetscErrorCode ierr; 3979 DM dm; 3980 3981 PetscFunctionBegin; 3982 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3983 if (r) { 3984 if (!snes->vec_func) { 3985 if (snes->vec_rhs) { 3986 ierr = VecDuplicate(snes->vec_rhs,&snes->vec_func);CHKERRQ(ierr); 3987 } else if (snes->vec_sol) { 3988 ierr = VecDuplicate(snes->vec_sol,&snes->vec_func);CHKERRQ(ierr); 3989 } else if (snes->dm) { 3990 ierr = DMCreateGlobalVector(snes->dm,&snes->vec_func);CHKERRQ(ierr); 3991 } 3992 } 3993 *r = snes->vec_func; 3994 } 3995 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 3996 ierr = DMSNESGetFunction(dm,f,ctx);CHKERRQ(ierr); 3997 PetscFunctionReturn(0); 3998 } 3999 4000 /*@C 4001 SNESGetNGS - Returns the NGS function and context. 4002 4003 Input Parameter: 4004 . snes - the SNES context 4005 4006 Output Parameter: 4007 + f - the function (or NULL) see SNESNGSFunction for details 4008 - ctx - the function context (or NULL) 4009 4010 Level: advanced 4011 4012 .keywords: SNES, nonlinear, get, function 4013 4014 .seealso: SNESSetNGS(), SNESGetFunction() 4015 @*/ 4016 4017 #undef __FUNCT__ 4018 #define __FUNCT__ "SNESGetNGS" 4019 PetscErrorCode SNESGetNGS (SNES snes, PetscErrorCode (**f)(SNES, Vec, Vec, void*), void ** ctx) 4020 { 4021 PetscErrorCode ierr; 4022 DM dm; 4023 4024 PetscFunctionBegin; 4025 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4026 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 4027 ierr = DMSNESGetNGS(dm,f,ctx);CHKERRQ(ierr); 4028 PetscFunctionReturn(0); 4029 } 4030 4031 #undef __FUNCT__ 4032 #define __FUNCT__ "SNESSetOptionsPrefix" 4033 /*@C 4034 SNESSetOptionsPrefix - Sets the prefix used for searching for all 4035 SNES options in the database. 4036 4037 Logically Collective on SNES 4038 4039 Input Parameter: 4040 + snes - the SNES context 4041 - prefix - the prefix to prepend to all option names 4042 4043 Notes: 4044 A hyphen (-) must NOT be given at the beginning of the prefix name. 4045 The first character of all runtime options is AUTOMATICALLY the hyphen. 4046 4047 Level: advanced 4048 4049 .keywords: SNES, set, options, prefix, database 4050 4051 .seealso: SNESSetFromOptions() 4052 @*/ 4053 PetscErrorCode SNESSetOptionsPrefix(SNES snes,const char prefix[]) 4054 { 4055 PetscErrorCode ierr; 4056 4057 PetscFunctionBegin; 4058 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4059 ierr = PetscObjectSetOptionsPrefix((PetscObject)snes,prefix);CHKERRQ(ierr); 4060 if (!snes->ksp) {ierr = SNESGetKSP(snes,&snes->ksp);CHKERRQ(ierr);} 4061 if (snes->linesearch) { 4062 ierr = SNESGetLineSearch(snes,&snes->linesearch);CHKERRQ(ierr); 4063 ierr = PetscObjectSetOptionsPrefix((PetscObject)snes->linesearch,prefix);CHKERRQ(ierr); 4064 } 4065 ierr = KSPSetOptionsPrefix(snes->ksp,prefix);CHKERRQ(ierr); 4066 PetscFunctionReturn(0); 4067 } 4068 4069 #undef __FUNCT__ 4070 #define __FUNCT__ "SNESAppendOptionsPrefix" 4071 /*@C 4072 SNESAppendOptionsPrefix - Appends to the prefix used for searching for all 4073 SNES options in the database. 4074 4075 Logically Collective on SNES 4076 4077 Input Parameters: 4078 + snes - the SNES context 4079 - prefix - the prefix to prepend to all option names 4080 4081 Notes: 4082 A hyphen (-) must NOT be given at the beginning of the prefix name. 4083 The first character of all runtime options is AUTOMATICALLY the hyphen. 4084 4085 Level: advanced 4086 4087 .keywords: SNES, append, options, prefix, database 4088 4089 .seealso: SNESGetOptionsPrefix() 4090 @*/ 4091 PetscErrorCode SNESAppendOptionsPrefix(SNES snes,const char prefix[]) 4092 { 4093 PetscErrorCode ierr; 4094 4095 PetscFunctionBegin; 4096 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4097 ierr = PetscObjectAppendOptionsPrefix((PetscObject)snes,prefix);CHKERRQ(ierr); 4098 if (!snes->ksp) {ierr = SNESGetKSP(snes,&snes->ksp);CHKERRQ(ierr);} 4099 if (snes->linesearch) { 4100 ierr = SNESGetLineSearch(snes,&snes->linesearch);CHKERRQ(ierr); 4101 ierr = PetscObjectAppendOptionsPrefix((PetscObject)snes->linesearch,prefix);CHKERRQ(ierr); 4102 } 4103 ierr = KSPAppendOptionsPrefix(snes->ksp,prefix);CHKERRQ(ierr); 4104 PetscFunctionReturn(0); 4105 } 4106 4107 #undef __FUNCT__ 4108 #define __FUNCT__ "SNESGetOptionsPrefix" 4109 /*@C 4110 SNESGetOptionsPrefix - Sets the prefix used for searching for all 4111 SNES options in the database. 4112 4113 Not Collective 4114 4115 Input Parameter: 4116 . snes - the SNES context 4117 4118 Output Parameter: 4119 . prefix - pointer to the prefix string used 4120 4121 Notes: On the fortran side, the user should pass in a string 'prefix' of 4122 sufficient length to hold the prefix. 4123 4124 Level: advanced 4125 4126 .keywords: SNES, get, options, prefix, database 4127 4128 .seealso: SNESAppendOptionsPrefix() 4129 @*/ 4130 PetscErrorCode SNESGetOptionsPrefix(SNES snes,const char *prefix[]) 4131 { 4132 PetscErrorCode ierr; 4133 4134 PetscFunctionBegin; 4135 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4136 ierr = PetscObjectGetOptionsPrefix((PetscObject)snes,prefix);CHKERRQ(ierr); 4137 PetscFunctionReturn(0); 4138 } 4139 4140 4141 #undef __FUNCT__ 4142 #define __FUNCT__ "SNESRegister" 4143 /*@C 4144 SNESRegister - Adds a method to the nonlinear solver package. 4145 4146 Not collective 4147 4148 Input Parameters: 4149 + name_solver - name of a new user-defined solver 4150 - routine_create - routine to create method context 4151 4152 Notes: 4153 SNESRegister() may be called multiple times to add several user-defined solvers. 4154 4155 Sample usage: 4156 .vb 4157 SNESRegister("my_solver",MySolverCreate); 4158 .ve 4159 4160 Then, your solver can be chosen with the procedural interface via 4161 $ SNESSetType(snes,"my_solver") 4162 or at runtime via the option 4163 $ -snes_type my_solver 4164 4165 Level: advanced 4166 4167 Note: If your function is not being put into a shared library then use SNESRegister() instead 4168 4169 .keywords: SNES, nonlinear, register 4170 4171 .seealso: SNESRegisterAll(), SNESRegisterDestroy() 4172 4173 Level: advanced 4174 @*/ 4175 PetscErrorCode SNESRegister(const char sname[],PetscErrorCode (*function)(SNES)) 4176 { 4177 PetscErrorCode ierr; 4178 4179 PetscFunctionBegin; 4180 ierr = PetscFunctionListAdd(&SNESList,sname,function);CHKERRQ(ierr); 4181 PetscFunctionReturn(0); 4182 } 4183 4184 #undef __FUNCT__ 4185 #define __FUNCT__ "SNESTestLocalMin" 4186 PetscErrorCode SNESTestLocalMin(SNES snes) 4187 { 4188 PetscErrorCode ierr; 4189 PetscInt N,i,j; 4190 Vec u,uh,fh; 4191 PetscScalar value; 4192 PetscReal norm; 4193 4194 PetscFunctionBegin; 4195 ierr = SNESGetSolution(snes,&u);CHKERRQ(ierr); 4196 ierr = VecDuplicate(u,&uh);CHKERRQ(ierr); 4197 ierr = VecDuplicate(u,&fh);CHKERRQ(ierr); 4198 4199 /* currently only works for sequential */ 4200 ierr = PetscPrintf(PETSC_COMM_WORLD,"Testing FormFunction() for local min\n");CHKERRQ(ierr); 4201 ierr = VecGetSize(u,&N);CHKERRQ(ierr); 4202 for (i=0; i<N; i++) { 4203 ierr = VecCopy(u,uh);CHKERRQ(ierr); 4204 ierr = PetscPrintf(PETSC_COMM_WORLD,"i = %D\n",i);CHKERRQ(ierr); 4205 for (j=-10; j<11; j++) { 4206 value = PetscSign(j)*PetscExpReal(PetscAbs(j)-10.0); 4207 ierr = VecSetValue(uh,i,value,ADD_VALUES);CHKERRQ(ierr); 4208 ierr = SNESComputeFunction(snes,uh,fh);CHKERRQ(ierr); 4209 ierr = VecNorm(fh,NORM_2,&norm);CHKERRQ(ierr); 4210 ierr = PetscPrintf(PETSC_COMM_WORLD," j norm %D %18.16e\n",j,norm);CHKERRQ(ierr); 4211 value = -value; 4212 ierr = VecSetValue(uh,i,value,ADD_VALUES);CHKERRQ(ierr); 4213 } 4214 } 4215 ierr = VecDestroy(&uh);CHKERRQ(ierr); 4216 ierr = VecDestroy(&fh);CHKERRQ(ierr); 4217 PetscFunctionReturn(0); 4218 } 4219 4220 #undef __FUNCT__ 4221 #define __FUNCT__ "SNESKSPSetUseEW" 4222 /*@ 4223 SNESKSPSetUseEW - Sets SNES use Eisenstat-Walker method for 4224 computing relative tolerance for linear solvers within an inexact 4225 Newton method. 4226 4227 Logically Collective on SNES 4228 4229 Input Parameters: 4230 + snes - SNES context 4231 - flag - PETSC_TRUE or PETSC_FALSE 4232 4233 Options Database: 4234 + -snes_ksp_ew - use Eisenstat-Walker method for determining linear system convergence 4235 . -snes_ksp_ew_version ver - version of Eisenstat-Walker method 4236 . -snes_ksp_ew_rtol0 <rtol0> - Sets rtol0 4237 . -snes_ksp_ew_rtolmax <rtolmax> - Sets rtolmax 4238 . -snes_ksp_ew_gamma <gamma> - Sets gamma 4239 . -snes_ksp_ew_alpha <alpha> - Sets alpha 4240 . -snes_ksp_ew_alpha2 <alpha2> - Sets alpha2 4241 - -snes_ksp_ew_threshold <threshold> - Sets threshold 4242 4243 Notes: 4244 Currently, the default is to use a constant relative tolerance for 4245 the inner linear solvers. Alternatively, one can use the 4246 Eisenstat-Walker method, where the relative convergence tolerance 4247 is reset at each Newton iteration according progress of the nonlinear 4248 solver. 4249 4250 Level: advanced 4251 4252 Reference: 4253 S. C. Eisenstat and H. F. Walker, "Choosing the forcing terms in an 4254 inexact Newton method", SISC 17 (1), pp.16-32, 1996. 4255 4256 .keywords: SNES, KSP, Eisenstat, Walker, convergence, test, inexact, Newton 4257 4258 .seealso: SNESKSPGetUseEW(), SNESKSPGetParametersEW(), SNESKSPSetParametersEW() 4259 @*/ 4260 PetscErrorCode SNESKSPSetUseEW(SNES snes,PetscBool flag) 4261 { 4262 PetscFunctionBegin; 4263 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4264 PetscValidLogicalCollectiveBool(snes,flag,2); 4265 snes->ksp_ewconv = flag; 4266 PetscFunctionReturn(0); 4267 } 4268 4269 #undef __FUNCT__ 4270 #define __FUNCT__ "SNESKSPGetUseEW" 4271 /*@ 4272 SNESKSPGetUseEW - Gets if SNES is using Eisenstat-Walker method 4273 for computing relative tolerance for linear solvers within an 4274 inexact Newton method. 4275 4276 Not Collective 4277 4278 Input Parameter: 4279 . snes - SNES context 4280 4281 Output Parameter: 4282 . flag - PETSC_TRUE or PETSC_FALSE 4283 4284 Notes: 4285 Currently, the default is to use a constant relative tolerance for 4286 the inner linear solvers. Alternatively, one can use the 4287 Eisenstat-Walker method, where the relative convergence tolerance 4288 is reset at each Newton iteration according progress of the nonlinear 4289 solver. 4290 4291 Level: advanced 4292 4293 Reference: 4294 S. C. Eisenstat and H. F. Walker, "Choosing the forcing terms in an 4295 inexact Newton method", SISC 17 (1), pp.16-32, 1996. 4296 4297 .keywords: SNES, KSP, Eisenstat, Walker, convergence, test, inexact, Newton 4298 4299 .seealso: SNESKSPSetUseEW(), SNESKSPGetParametersEW(), SNESKSPSetParametersEW() 4300 @*/ 4301 PetscErrorCode SNESKSPGetUseEW(SNES snes, PetscBool *flag) 4302 { 4303 PetscFunctionBegin; 4304 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4305 PetscValidPointer(flag,2); 4306 *flag = snes->ksp_ewconv; 4307 PetscFunctionReturn(0); 4308 } 4309 4310 #undef __FUNCT__ 4311 #define __FUNCT__ "SNESKSPSetParametersEW" 4312 /*@ 4313 SNESKSPSetParametersEW - Sets parameters for Eisenstat-Walker 4314 convergence criteria for the linear solvers within an inexact 4315 Newton method. 4316 4317 Logically Collective on SNES 4318 4319 Input Parameters: 4320 + snes - SNES context 4321 . version - version 1, 2 (default is 2) or 3 4322 . rtol_0 - initial relative tolerance (0 <= rtol_0 < 1) 4323 . rtol_max - maximum relative tolerance (0 <= rtol_max < 1) 4324 . gamma - multiplicative factor for version 2 rtol computation 4325 (0 <= gamma2 <= 1) 4326 . alpha - power for version 2 rtol computation (1 < alpha <= 2) 4327 . alpha2 - power for safeguard 4328 - threshold - threshold for imposing safeguard (0 < threshold < 1) 4329 4330 Note: 4331 Version 3 was contributed by Luis Chacon, June 2006. 4332 4333 Use PETSC_DEFAULT to retain the default for any of the parameters. 4334 4335 Level: advanced 4336 4337 Reference: 4338 S. C. Eisenstat and H. F. Walker, "Choosing the forcing terms in an 4339 inexact Newton method", Utah State University Math. Stat. Dept. Res. 4340 Report 6/94/75, June, 1994, to appear in SIAM J. Sci. Comput. 4341 4342 .keywords: SNES, KSP, Eisenstat, Walker, set, parameters 4343 4344 .seealso: SNESKSPSetUseEW(), SNESKSPGetUseEW(), SNESKSPGetParametersEW() 4345 @*/ 4346 PetscErrorCode SNESKSPSetParametersEW(SNES snes,PetscInt version,PetscReal rtol_0,PetscReal rtol_max,PetscReal gamma,PetscReal alpha,PetscReal alpha2,PetscReal threshold) 4347 { 4348 SNESKSPEW *kctx; 4349 4350 PetscFunctionBegin; 4351 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4352 kctx = (SNESKSPEW*)snes->kspconvctx; 4353 if (!kctx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"No Eisenstat-Walker context existing"); 4354 PetscValidLogicalCollectiveInt(snes,version,2); 4355 PetscValidLogicalCollectiveReal(snes,rtol_0,3); 4356 PetscValidLogicalCollectiveReal(snes,rtol_max,4); 4357 PetscValidLogicalCollectiveReal(snes,gamma,5); 4358 PetscValidLogicalCollectiveReal(snes,alpha,6); 4359 PetscValidLogicalCollectiveReal(snes,alpha2,7); 4360 PetscValidLogicalCollectiveReal(snes,threshold,8); 4361 4362 if (version != PETSC_DEFAULT) kctx->version = version; 4363 if (rtol_0 != PETSC_DEFAULT) kctx->rtol_0 = rtol_0; 4364 if (rtol_max != PETSC_DEFAULT) kctx->rtol_max = rtol_max; 4365 if (gamma != PETSC_DEFAULT) kctx->gamma = gamma; 4366 if (alpha != PETSC_DEFAULT) kctx->alpha = alpha; 4367 if (alpha2 != PETSC_DEFAULT) kctx->alpha2 = alpha2; 4368 if (threshold != PETSC_DEFAULT) kctx->threshold = threshold; 4369 4370 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); 4371 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); 4372 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); 4373 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); 4374 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); 4375 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); 4376 PetscFunctionReturn(0); 4377 } 4378 4379 #undef __FUNCT__ 4380 #define __FUNCT__ "SNESKSPGetParametersEW" 4381 /*@ 4382 SNESKSPGetParametersEW - Gets parameters for Eisenstat-Walker 4383 convergence criteria for the linear solvers within an inexact 4384 Newton method. 4385 4386 Not Collective 4387 4388 Input Parameters: 4389 snes - SNES context 4390 4391 Output Parameters: 4392 + version - version 1, 2 (default is 2) or 3 4393 . rtol_0 - initial relative tolerance (0 <= rtol_0 < 1) 4394 . rtol_max - maximum relative tolerance (0 <= rtol_max < 1) 4395 . gamma - multiplicative factor for version 2 rtol computation (0 <= gamma2 <= 1) 4396 . alpha - power for version 2 rtol computation (1 < alpha <= 2) 4397 . alpha2 - power for safeguard 4398 - threshold - threshold for imposing safeguard (0 < threshold < 1) 4399 4400 Level: advanced 4401 4402 .keywords: SNES, KSP, Eisenstat, Walker, get, parameters 4403 4404 .seealso: SNESKSPSetUseEW(), SNESKSPGetUseEW(), SNESKSPSetParametersEW() 4405 @*/ 4406 PetscErrorCode SNESKSPGetParametersEW(SNES snes,PetscInt *version,PetscReal *rtol_0,PetscReal *rtol_max,PetscReal *gamma,PetscReal *alpha,PetscReal *alpha2,PetscReal *threshold) 4407 { 4408 SNESKSPEW *kctx; 4409 4410 PetscFunctionBegin; 4411 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4412 kctx = (SNESKSPEW*)snes->kspconvctx; 4413 if (!kctx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"No Eisenstat-Walker context existing"); 4414 if (version) *version = kctx->version; 4415 if (rtol_0) *rtol_0 = kctx->rtol_0; 4416 if (rtol_max) *rtol_max = kctx->rtol_max; 4417 if (gamma) *gamma = kctx->gamma; 4418 if (alpha) *alpha = kctx->alpha; 4419 if (alpha2) *alpha2 = kctx->alpha2; 4420 if (threshold) *threshold = kctx->threshold; 4421 PetscFunctionReturn(0); 4422 } 4423 4424 #undef __FUNCT__ 4425 #define __FUNCT__ "KSPPreSolve_SNESEW" 4426 PetscErrorCode KSPPreSolve_SNESEW(KSP ksp, Vec b, Vec x, SNES snes) 4427 { 4428 PetscErrorCode ierr; 4429 SNESKSPEW *kctx = (SNESKSPEW*)snes->kspconvctx; 4430 PetscReal rtol = PETSC_DEFAULT,stol; 4431 4432 PetscFunctionBegin; 4433 if (!snes->ksp_ewconv) PetscFunctionReturn(0); 4434 if (!snes->iter) { 4435 rtol = kctx->rtol_0; /* first time in, so use the original user rtol */ 4436 ierr = VecNorm(snes->vec_func,NORM_2,&kctx->norm_first);CHKERRQ(ierr); 4437 } 4438 else { 4439 if (kctx->version == 1) { 4440 rtol = (snes->norm - kctx->lresid_last)/kctx->norm_last; 4441 if (rtol < 0.0) rtol = -rtol; 4442 stol = PetscPowReal(kctx->rtol_last,kctx->alpha2); 4443 if (stol > kctx->threshold) rtol = PetscMax(rtol,stol); 4444 } else if (kctx->version == 2) { 4445 rtol = kctx->gamma * PetscPowReal(snes->norm/kctx->norm_last,kctx->alpha); 4446 stol = kctx->gamma * PetscPowReal(kctx->rtol_last,kctx->alpha); 4447 if (stol > kctx->threshold) rtol = PetscMax(rtol,stol); 4448 } else if (kctx->version == 3) { /* contributed by Luis Chacon, June 2006. */ 4449 rtol = kctx->gamma * PetscPowReal(snes->norm/kctx->norm_last,kctx->alpha); 4450 /* safeguard: avoid sharp decrease of rtol */ 4451 stol = kctx->gamma*PetscPowReal(kctx->rtol_last,kctx->alpha); 4452 stol = PetscMax(rtol,stol); 4453 rtol = PetscMin(kctx->rtol_0,stol); 4454 /* safeguard: avoid oversolving */ 4455 stol = kctx->gamma*(kctx->norm_first*snes->rtol)/snes->norm; 4456 stol = PetscMax(rtol,stol); 4457 rtol = PetscMin(kctx->rtol_0,stol); 4458 } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Only versions 1, 2 or 3 are supported: %D",kctx->version); 4459 } 4460 /* safeguard: avoid rtol greater than one */ 4461 rtol = PetscMin(rtol,kctx->rtol_max); 4462 ierr = KSPSetTolerances(ksp,rtol,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);CHKERRQ(ierr); 4463 ierr = PetscInfo3(snes,"iter %D, Eisenstat-Walker (version %D) KSP rtol=%g\n",snes->iter,kctx->version,(double)rtol);CHKERRQ(ierr); 4464 PetscFunctionReturn(0); 4465 } 4466 4467 #undef __FUNCT__ 4468 #define __FUNCT__ "KSPPostSolve_SNESEW" 4469 PetscErrorCode KSPPostSolve_SNESEW(KSP ksp, Vec b, Vec x, SNES snes) 4470 { 4471 PetscErrorCode ierr; 4472 SNESKSPEW *kctx = (SNESKSPEW*)snes->kspconvctx; 4473 PCSide pcside; 4474 Vec lres; 4475 4476 PetscFunctionBegin; 4477 if (!snes->ksp_ewconv) PetscFunctionReturn(0); 4478 ierr = KSPGetTolerances(ksp,&kctx->rtol_last,0,0,0);CHKERRQ(ierr); 4479 kctx->norm_last = snes->norm; 4480 if (kctx->version == 1) { 4481 ierr = KSPGetPCSide(ksp,&pcside);CHKERRQ(ierr); 4482 if (pcside == PC_RIGHT) { /* XXX Should we also test KSP_UNPRECONDITIONED_NORM ? */ 4483 /* KSP residual is true linear residual */ 4484 ierr = KSPGetResidualNorm(ksp,&kctx->lresid_last);CHKERRQ(ierr); 4485 } else { 4486 /* KSP residual is preconditioned residual */ 4487 /* compute true linear residual norm */ 4488 ierr = VecDuplicate(b,&lres);CHKERRQ(ierr); 4489 ierr = MatMult(snes->jacobian,x,lres);CHKERRQ(ierr); 4490 ierr = VecAYPX(lres,-1.0,b);CHKERRQ(ierr); 4491 ierr = VecNorm(lres,NORM_2,&kctx->lresid_last);CHKERRQ(ierr); 4492 ierr = VecDestroy(&lres);CHKERRQ(ierr); 4493 } 4494 } 4495 PetscFunctionReturn(0); 4496 } 4497 4498 #undef __FUNCT__ 4499 #define __FUNCT__ "SNESGetKSP" 4500 /*@ 4501 SNESGetKSP - Returns the KSP context for a SNES solver. 4502 4503 Not Collective, but if SNES object is parallel, then KSP object is parallel 4504 4505 Input Parameter: 4506 . snes - the SNES context 4507 4508 Output Parameter: 4509 . ksp - the KSP context 4510 4511 Notes: 4512 The user can then directly manipulate the KSP context to set various 4513 options, etc. Likewise, the user can then extract and manipulate the 4514 PC contexts as well. 4515 4516 Level: beginner 4517 4518 .keywords: SNES, nonlinear, get, KSP, context 4519 4520 .seealso: KSPGetPC(), SNESCreate(), KSPCreate(), SNESSetKSP() 4521 @*/ 4522 PetscErrorCode SNESGetKSP(SNES snes,KSP *ksp) 4523 { 4524 PetscErrorCode ierr; 4525 4526 PetscFunctionBegin; 4527 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4528 PetscValidPointer(ksp,2); 4529 4530 if (!snes->ksp) { 4531 ierr = KSPCreate(PetscObjectComm((PetscObject)snes),&snes->ksp);CHKERRQ(ierr); 4532 ierr = PetscObjectIncrementTabLevel((PetscObject)snes->ksp,(PetscObject)snes,1);CHKERRQ(ierr); 4533 ierr = PetscLogObjectParent((PetscObject)snes,(PetscObject)snes->ksp);CHKERRQ(ierr); 4534 4535 ierr = KSPSetPreSolve(snes->ksp,(PetscErrorCode (*)(KSP,Vec,Vec,void*))KSPPreSolve_SNESEW,snes);CHKERRQ(ierr); 4536 ierr = KSPSetPostSolve(snes->ksp,(PetscErrorCode (*)(KSP,Vec,Vec,void*))KSPPostSolve_SNESEW,snes);CHKERRQ(ierr); 4537 } 4538 *ksp = snes->ksp; 4539 PetscFunctionReturn(0); 4540 } 4541 4542 4543 #include <petsc-private/dmimpl.h> 4544 #undef __FUNCT__ 4545 #define __FUNCT__ "SNESSetDM" 4546 /*@ 4547 SNESSetDM - Sets the DM that may be used by some preconditioners 4548 4549 Logically Collective on SNES 4550 4551 Input Parameters: 4552 + snes - the preconditioner context 4553 - dm - the dm 4554 4555 Level: intermediate 4556 4557 .seealso: SNESGetDM(), KSPSetDM(), KSPGetDM() 4558 @*/ 4559 PetscErrorCode SNESSetDM(SNES snes,DM dm) 4560 { 4561 PetscErrorCode ierr; 4562 KSP ksp; 4563 DMSNES sdm; 4564 4565 PetscFunctionBegin; 4566 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4567 if (dm) {ierr = PetscObjectReference((PetscObject)dm);CHKERRQ(ierr);} 4568 if (snes->dm) { /* Move the DMSNES context over to the new DM unless the new DM already has one */ 4569 if (snes->dm->dmsnes && snes->dmAuto && !dm->dmsnes) { 4570 ierr = DMCopyDMSNES(snes->dm,dm);CHKERRQ(ierr); 4571 ierr = DMGetDMSNES(snes->dm,&sdm);CHKERRQ(ierr); 4572 if (sdm->originaldm == snes->dm) sdm->originaldm = dm; /* Grant write privileges to the replacement DM */ 4573 } 4574 ierr = DMDestroy(&snes->dm);CHKERRQ(ierr); 4575 } 4576 snes->dm = dm; 4577 snes->dmAuto = PETSC_FALSE; 4578 4579 ierr = SNESGetKSP(snes,&ksp);CHKERRQ(ierr); 4580 ierr = KSPSetDM(ksp,dm);CHKERRQ(ierr); 4581 ierr = KSPSetDMActive(ksp,PETSC_FALSE);CHKERRQ(ierr); 4582 if (snes->pc) { 4583 ierr = SNESSetDM(snes->pc, snes->dm);CHKERRQ(ierr); 4584 ierr = SNESSetNPCSide(snes,snes->pcside);CHKERRQ(ierr); 4585 } 4586 PetscFunctionReturn(0); 4587 } 4588 4589 #undef __FUNCT__ 4590 #define __FUNCT__ "SNESGetDM" 4591 /*@ 4592 SNESGetDM - Gets the DM that may be used by some preconditioners 4593 4594 Not Collective but DM obtained is parallel on SNES 4595 4596 Input Parameter: 4597 . snes - the preconditioner context 4598 4599 Output Parameter: 4600 . dm - the dm 4601 4602 Level: intermediate 4603 4604 .seealso: SNESSetDM(), KSPSetDM(), KSPGetDM() 4605 @*/ 4606 PetscErrorCode SNESGetDM(SNES snes,DM *dm) 4607 { 4608 PetscErrorCode ierr; 4609 4610 PetscFunctionBegin; 4611 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4612 if (!snes->dm) { 4613 ierr = DMShellCreate(PetscObjectComm((PetscObject)snes),&snes->dm);CHKERRQ(ierr); 4614 snes->dmAuto = PETSC_TRUE; 4615 } 4616 *dm = snes->dm; 4617 PetscFunctionReturn(0); 4618 } 4619 4620 #undef __FUNCT__ 4621 #define __FUNCT__ "SNESSetNPC" 4622 /*@ 4623 SNESSetNPC - Sets the nonlinear preconditioner to be used. 4624 4625 Collective on SNES 4626 4627 Input Parameters: 4628 + snes - iterative context obtained from SNESCreate() 4629 - pc - the preconditioner object 4630 4631 Notes: 4632 Use SNESGetNPC() to retrieve the preconditioner context (for example, 4633 to configure it using the API). 4634 4635 Level: developer 4636 4637 .keywords: SNES, set, precondition 4638 .seealso: SNESGetNPC() 4639 @*/ 4640 PetscErrorCode SNESSetNPC(SNES snes, SNES pc) 4641 { 4642 PetscErrorCode ierr; 4643 4644 PetscFunctionBegin; 4645 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); 4646 PetscValidHeaderSpecific(pc, SNES_CLASSID, 2); 4647 PetscCheckSameComm(snes, 1, pc, 2); 4648 ierr = PetscObjectReference((PetscObject) pc);CHKERRQ(ierr); 4649 ierr = SNESDestroy(&snes->pc);CHKERRQ(ierr); 4650 snes->pc = pc; 4651 ierr = PetscLogObjectParent((PetscObject)snes, (PetscObject)snes->pc);CHKERRQ(ierr); 4652 PetscFunctionReturn(0); 4653 } 4654 4655 #undef __FUNCT__ 4656 #define __FUNCT__ "SNESGetNPC" 4657 /*@ 4658 SNESGetNPC - Creates a nonlinear preconditioning solver (SNES) to be used to precondition the nonlinear solver. 4659 4660 Not Collective 4661 4662 Input Parameter: 4663 . snes - iterative context obtained from SNESCreate() 4664 4665 Output Parameter: 4666 . pc - preconditioner context 4667 4668 Notes: If a SNES was previously set with SNESSetNPC() then that SNES is returned. 4669 4670 Level: developer 4671 4672 .keywords: SNES, get, preconditioner 4673 .seealso: SNESSetNPC() 4674 @*/ 4675 PetscErrorCode SNESGetNPC(SNES snes, SNES *pc) 4676 { 4677 PetscErrorCode ierr; 4678 const char *optionsprefix; 4679 4680 PetscFunctionBegin; 4681 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); 4682 PetscValidPointer(pc, 2); 4683 if (!snes->pc) { 4684 ierr = SNESCreate(PetscObjectComm((PetscObject)snes),&snes->pc);CHKERRQ(ierr); 4685 ierr = PetscObjectIncrementTabLevel((PetscObject)snes->pc,(PetscObject)snes,1);CHKERRQ(ierr); 4686 ierr = PetscLogObjectParent((PetscObject)snes,(PetscObject)snes->pc);CHKERRQ(ierr); 4687 ierr = SNESGetOptionsPrefix(snes,&optionsprefix);CHKERRQ(ierr); 4688 ierr = SNESSetOptionsPrefix(snes->pc,optionsprefix);CHKERRQ(ierr); 4689 ierr = SNESAppendOptionsPrefix(snes->pc,"npc_");CHKERRQ(ierr); 4690 ierr = SNESSetCountersReset(snes->pc,PETSC_FALSE);CHKERRQ(ierr); 4691 } 4692 *pc = snes->pc; 4693 PetscFunctionReturn(0); 4694 } 4695 4696 #undef __FUNCT__ 4697 #define __FUNCT__ "SNESSetNPCSide" 4698 /*@ 4699 SNESSetNPCSide - Sets the preconditioning side. 4700 4701 Logically Collective on SNES 4702 4703 Input Parameter: 4704 . snes - iterative context obtained from SNESCreate() 4705 4706 Output Parameter: 4707 . side - the preconditioning side, where side is one of 4708 .vb 4709 PC_LEFT - left preconditioning (default) 4710 PC_RIGHT - right preconditioning 4711 .ve 4712 4713 Options Database Keys: 4714 . -snes_pc_side <right,left> 4715 4716 Level: intermediate 4717 4718 .keywords: SNES, set, right, left, side, preconditioner, flag 4719 4720 .seealso: SNESGetNPCSide(), KSPSetPCSide() 4721 @*/ 4722 PetscErrorCode SNESSetNPCSide(SNES snes,PCSide side) 4723 { 4724 PetscFunctionBegin; 4725 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4726 PetscValidLogicalCollectiveEnum(snes,side,2); 4727 snes->pcside = side; 4728 PetscFunctionReturn(0); 4729 } 4730 4731 #undef __FUNCT__ 4732 #define __FUNCT__ "SNESGetNPCSide" 4733 /*@ 4734 SNESGetNPCSide - Gets the preconditioning side. 4735 4736 Not Collective 4737 4738 Input Parameter: 4739 . snes - iterative context obtained from SNESCreate() 4740 4741 Output Parameter: 4742 . side - the preconditioning side, where side is one of 4743 .vb 4744 PC_LEFT - left preconditioning (default) 4745 PC_RIGHT - right preconditioning 4746 .ve 4747 4748 Level: intermediate 4749 4750 .keywords: SNES, get, right, left, side, preconditioner, flag 4751 4752 .seealso: SNESSetNPCSide(), KSPGetPCSide() 4753 @*/ 4754 PetscErrorCode SNESGetNPCSide(SNES snes,PCSide *side) 4755 { 4756 PetscFunctionBegin; 4757 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4758 PetscValidPointer(side,2); 4759 *side = snes->pcside; 4760 PetscFunctionReturn(0); 4761 } 4762 4763 #undef __FUNCT__ 4764 #define __FUNCT__ "SNESSetLineSearch" 4765 /*@ 4766 SNESSetLineSearch - Sets the linesearch on the SNES instance. 4767 4768 Collective on SNES 4769 4770 Input Parameters: 4771 + snes - iterative context obtained from SNESCreate() 4772 - linesearch - the linesearch object 4773 4774 Notes: 4775 Use SNESGetLineSearch() to retrieve the preconditioner context (for example, 4776 to configure it using the API). 4777 4778 Level: developer 4779 4780 .keywords: SNES, set, linesearch 4781 .seealso: SNESGetLineSearch() 4782 @*/ 4783 PetscErrorCode SNESSetLineSearch(SNES snes, SNESLineSearch linesearch) 4784 { 4785 PetscErrorCode ierr; 4786 4787 PetscFunctionBegin; 4788 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); 4789 PetscValidHeaderSpecific(linesearch, SNESLINESEARCH_CLASSID, 2); 4790 PetscCheckSameComm(snes, 1, linesearch, 2); 4791 ierr = PetscObjectReference((PetscObject) linesearch);CHKERRQ(ierr); 4792 ierr = SNESLineSearchDestroy(&snes->linesearch);CHKERRQ(ierr); 4793 4794 snes->linesearch = linesearch; 4795 4796 ierr = PetscLogObjectParent((PetscObject)snes, (PetscObject)snes->linesearch);CHKERRQ(ierr); 4797 PetscFunctionReturn(0); 4798 } 4799 4800 #undef __FUNCT__ 4801 #define __FUNCT__ "SNESGetLineSearch" 4802 /*@ 4803 SNESGetLineSearch - Returns a pointer to the line search context set with SNESSetLineSearch() 4804 or creates a default line search instance associated with the SNES and returns it. 4805 4806 Not Collective 4807 4808 Input Parameter: 4809 . snes - iterative context obtained from SNESCreate() 4810 4811 Output Parameter: 4812 . linesearch - linesearch context 4813 4814 Level: beginner 4815 4816 .keywords: SNES, get, linesearch 4817 .seealso: SNESSetLineSearch(), SNESLineSearchCreate() 4818 @*/ 4819 PetscErrorCode SNESGetLineSearch(SNES snes, SNESLineSearch *linesearch) 4820 { 4821 PetscErrorCode ierr; 4822 const char *optionsprefix; 4823 4824 PetscFunctionBegin; 4825 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); 4826 PetscValidPointer(linesearch, 2); 4827 if (!snes->linesearch) { 4828 ierr = SNESGetOptionsPrefix(snes, &optionsprefix);CHKERRQ(ierr); 4829 ierr = SNESLineSearchCreate(PetscObjectComm((PetscObject)snes), &snes->linesearch);CHKERRQ(ierr); 4830 ierr = SNESLineSearchSetSNES(snes->linesearch, snes);CHKERRQ(ierr); 4831 ierr = SNESLineSearchAppendOptionsPrefix(snes->linesearch, optionsprefix);CHKERRQ(ierr); 4832 ierr = PetscObjectIncrementTabLevel((PetscObject) snes->linesearch, (PetscObject) snes, 1);CHKERRQ(ierr); 4833 ierr = PetscLogObjectParent((PetscObject)snes, (PetscObject)snes->linesearch);CHKERRQ(ierr); 4834 } 4835 *linesearch = snes->linesearch; 4836 PetscFunctionReturn(0); 4837 } 4838 4839 #if defined(PETSC_HAVE_MATLAB_ENGINE) 4840 #include <mex.h> 4841 4842 typedef struct {char *funcname; mxArray *ctx;} SNESMatlabContext; 4843 4844 #undef __FUNCT__ 4845 #define __FUNCT__ "SNESComputeFunction_Matlab" 4846 /* 4847 SNESComputeFunction_Matlab - Calls the function that has been set with SNESSetFunctionMatlab(). 4848 4849 Collective on SNES 4850 4851 Input Parameters: 4852 + snes - the SNES context 4853 - x - input vector 4854 4855 Output Parameter: 4856 . y - function vector, as set by SNESSetFunction() 4857 4858 Notes: 4859 SNESComputeFunction() is typically used within nonlinear solvers 4860 implementations, so most users would not generally call this routine 4861 themselves. 4862 4863 Level: developer 4864 4865 .keywords: SNES, nonlinear, compute, function 4866 4867 .seealso: SNESSetFunction(), SNESGetFunction() 4868 */ 4869 PetscErrorCode SNESComputeFunction_Matlab(SNES snes,Vec x,Vec y, void *ctx) 4870 { 4871 PetscErrorCode ierr; 4872 SNESMatlabContext *sctx = (SNESMatlabContext*)ctx; 4873 int nlhs = 1,nrhs = 5; 4874 mxArray *plhs[1],*prhs[5]; 4875 long long int lx = 0,ly = 0,ls = 0; 4876 4877 PetscFunctionBegin; 4878 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4879 PetscValidHeaderSpecific(x,VEC_CLASSID,2); 4880 PetscValidHeaderSpecific(y,VEC_CLASSID,3); 4881 PetscCheckSameComm(snes,1,x,2); 4882 PetscCheckSameComm(snes,1,y,3); 4883 4884 /* call Matlab function in ctx with arguments x and y */ 4885 4886 ierr = PetscMemcpy(&ls,&snes,sizeof(snes));CHKERRQ(ierr); 4887 ierr = PetscMemcpy(&lx,&x,sizeof(x));CHKERRQ(ierr); 4888 ierr = PetscMemcpy(&ly,&y,sizeof(x));CHKERRQ(ierr); 4889 prhs[0] = mxCreateDoubleScalar((double)ls); 4890 prhs[1] = mxCreateDoubleScalar((double)lx); 4891 prhs[2] = mxCreateDoubleScalar((double)ly); 4892 prhs[3] = mxCreateString(sctx->funcname); 4893 prhs[4] = sctx->ctx; 4894 ierr = mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscSNESComputeFunctionInternal");CHKERRQ(ierr); 4895 ierr = mxGetScalar(plhs[0]);CHKERRQ(ierr); 4896 mxDestroyArray(prhs[0]); 4897 mxDestroyArray(prhs[1]); 4898 mxDestroyArray(prhs[2]); 4899 mxDestroyArray(prhs[3]); 4900 mxDestroyArray(plhs[0]); 4901 PetscFunctionReturn(0); 4902 } 4903 4904 #undef __FUNCT__ 4905 #define __FUNCT__ "SNESSetFunctionMatlab" 4906 /* 4907 SNESSetFunctionMatlab - Sets the function evaluation routine and function 4908 vector for use by the SNES routines in solving systems of nonlinear 4909 equations from MATLAB. Here the function is a string containing the name of a MATLAB function 4910 4911 Logically Collective on SNES 4912 4913 Input Parameters: 4914 + snes - the SNES context 4915 . r - vector to store function value 4916 - f - function evaluation routine 4917 4918 Notes: 4919 The Newton-like methods typically solve linear systems of the form 4920 $ f'(x) x = -f(x), 4921 where f'(x) denotes the Jacobian matrix and f(x) is the function. 4922 4923 Level: beginner 4924 4925 Developer Note: This bleeds the allocated memory SNESMatlabContext *sctx; 4926 4927 .keywords: SNES, nonlinear, set, function 4928 4929 .seealso: SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESSetFunction() 4930 */ 4931 PetscErrorCode SNESSetFunctionMatlab(SNES snes,Vec r,const char *f,mxArray *ctx) 4932 { 4933 PetscErrorCode ierr; 4934 SNESMatlabContext *sctx; 4935 4936 PetscFunctionBegin; 4937 /* currently sctx is memory bleed */ 4938 ierr = PetscMalloc(sizeof(SNESMatlabContext),&sctx);CHKERRQ(ierr); 4939 ierr = PetscStrallocpy(f,&sctx->funcname);CHKERRQ(ierr); 4940 /* 4941 This should work, but it doesn't 4942 sctx->ctx = ctx; 4943 mexMakeArrayPersistent(sctx->ctx); 4944 */ 4945 sctx->ctx = mxDuplicateArray(ctx); 4946 ierr = SNESSetFunction(snes,r,SNESComputeFunction_Matlab,sctx);CHKERRQ(ierr); 4947 PetscFunctionReturn(0); 4948 } 4949 4950 #undef __FUNCT__ 4951 #define __FUNCT__ "SNESComputeJacobian_Matlab" 4952 /* 4953 SNESComputeJacobian_Matlab - Calls the function that has been set with SNESSetJacobianMatlab(). 4954 4955 Collective on SNES 4956 4957 Input Parameters: 4958 + snes - the SNES context 4959 . x - input vector 4960 . A, B - the matrices 4961 - ctx - user context 4962 4963 Output Parameter: 4964 . flag - structure of the matrix 4965 4966 Level: developer 4967 4968 .keywords: SNES, nonlinear, compute, function 4969 4970 .seealso: SNESSetFunction(), SNESGetFunction() 4971 @*/ 4972 PetscErrorCode SNESComputeJacobian_Matlab(SNES snes,Vec x,Mat A,Mat B,MatStructure *flag, void *ctx) 4973 { 4974 PetscErrorCode ierr; 4975 SNESMatlabContext *sctx = (SNESMatlabContext*)ctx; 4976 int nlhs = 2,nrhs = 6; 4977 mxArray *plhs[2],*prhs[6]; 4978 long long int lx = 0,lA = 0,ls = 0, lB = 0; 4979 4980 PetscFunctionBegin; 4981 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4982 PetscValidHeaderSpecific(x,VEC_CLASSID,2); 4983 4984 /* call Matlab function in ctx with arguments x and y */ 4985 4986 ierr = PetscMemcpy(&ls,&snes,sizeof(snes));CHKERRQ(ierr); 4987 ierr = PetscMemcpy(&lx,&x,sizeof(x));CHKERRQ(ierr); 4988 ierr = PetscMemcpy(&lA,A,sizeof(x));CHKERRQ(ierr); 4989 ierr = PetscMemcpy(&lB,B,sizeof(x));CHKERRQ(ierr); 4990 prhs[0] = mxCreateDoubleScalar((double)ls); 4991 prhs[1] = mxCreateDoubleScalar((double)lx); 4992 prhs[2] = mxCreateDoubleScalar((double)lA); 4993 prhs[3] = mxCreateDoubleScalar((double)lB); 4994 prhs[4] = mxCreateString(sctx->funcname); 4995 prhs[5] = sctx->ctx; 4996 ierr = mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscSNESComputeJacobianInternal");CHKERRQ(ierr); 4997 ierr = mxGetScalar(plhs[0]);CHKERRQ(ierr); 4998 *flag = (MatStructure) mxGetScalar(plhs[1]);CHKERRQ(ierr); 4999 mxDestroyArray(prhs[0]); 5000 mxDestroyArray(prhs[1]); 5001 mxDestroyArray(prhs[2]); 5002 mxDestroyArray(prhs[3]); 5003 mxDestroyArray(prhs[4]); 5004 mxDestroyArray(plhs[0]); 5005 mxDestroyArray(plhs[1]); 5006 PetscFunctionReturn(0); 5007 } 5008 5009 #undef __FUNCT__ 5010 #define __FUNCT__ "SNESSetJacobianMatlab" 5011 /* 5012 SNESSetJacobianMatlab - Sets the Jacobian function evaluation routine and two empty Jacobian matrices 5013 vector for use by the SNES routines in solving systems of nonlinear 5014 equations from MATLAB. Here the function is a string containing the name of a MATLAB function 5015 5016 Logically Collective on SNES 5017 5018 Input Parameters: 5019 + snes - the SNES context 5020 . A,B - Jacobian matrices 5021 . J - function evaluation routine 5022 - ctx - user context 5023 5024 Level: developer 5025 5026 Developer Note: This bleeds the allocated memory SNESMatlabContext *sctx; 5027 5028 .keywords: SNES, nonlinear, set, function 5029 5030 .seealso: SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESSetFunction(), J 5031 */ 5032 PetscErrorCode SNESSetJacobianMatlab(SNES snes,Mat A,Mat B,const char *J,mxArray *ctx) 5033 { 5034 PetscErrorCode ierr; 5035 SNESMatlabContext *sctx; 5036 5037 PetscFunctionBegin; 5038 /* currently sctx is memory bleed */ 5039 ierr = PetscMalloc(sizeof(SNESMatlabContext),&sctx);CHKERRQ(ierr); 5040 ierr = PetscStrallocpy(J,&sctx->funcname);CHKERRQ(ierr); 5041 /* 5042 This should work, but it doesn't 5043 sctx->ctx = ctx; 5044 mexMakeArrayPersistent(sctx->ctx); 5045 */ 5046 sctx->ctx = mxDuplicateArray(ctx); 5047 ierr = SNESSetJacobian(snes,A,B,SNESComputeJacobian_Matlab,sctx);CHKERRQ(ierr); 5048 PetscFunctionReturn(0); 5049 } 5050 5051 #undef __FUNCT__ 5052 #define __FUNCT__ "SNESMonitor_Matlab" 5053 /* 5054 SNESMonitor_Matlab - Calls the function that has been set with SNESMonitorSetMatlab(). 5055 5056 Collective on SNES 5057 5058 .seealso: SNESSetFunction(), SNESGetFunction() 5059 @*/ 5060 PetscErrorCode SNESMonitor_Matlab(SNES snes,PetscInt it, PetscReal fnorm, void *ctx) 5061 { 5062 PetscErrorCode ierr; 5063 SNESMatlabContext *sctx = (SNESMatlabContext*)ctx; 5064 int nlhs = 1,nrhs = 6; 5065 mxArray *plhs[1],*prhs[6]; 5066 long long int lx = 0,ls = 0; 5067 Vec x = snes->vec_sol; 5068 5069 PetscFunctionBegin; 5070 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 5071 5072 ierr = PetscMemcpy(&ls,&snes,sizeof(snes));CHKERRQ(ierr); 5073 ierr = PetscMemcpy(&lx,&x,sizeof(x));CHKERRQ(ierr); 5074 prhs[0] = mxCreateDoubleScalar((double)ls); 5075 prhs[1] = mxCreateDoubleScalar((double)it); 5076 prhs[2] = mxCreateDoubleScalar((double)fnorm); 5077 prhs[3] = mxCreateDoubleScalar((double)lx); 5078 prhs[4] = mxCreateString(sctx->funcname); 5079 prhs[5] = sctx->ctx; 5080 ierr = mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscSNESMonitorInternal");CHKERRQ(ierr); 5081 ierr = mxGetScalar(plhs[0]);CHKERRQ(ierr); 5082 mxDestroyArray(prhs[0]); 5083 mxDestroyArray(prhs[1]); 5084 mxDestroyArray(prhs[2]); 5085 mxDestroyArray(prhs[3]); 5086 mxDestroyArray(prhs[4]); 5087 mxDestroyArray(plhs[0]); 5088 PetscFunctionReturn(0); 5089 } 5090 5091 #undef __FUNCT__ 5092 #define __FUNCT__ "SNESMonitorSetMatlab" 5093 /* 5094 SNESMonitorSetMatlab - Sets the monitor function from MATLAB 5095 5096 Level: developer 5097 5098 Developer Note: This bleeds the allocated memory SNESMatlabContext *sctx; 5099 5100 .keywords: SNES, nonlinear, set, function 5101 5102 .seealso: SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESSetFunction() 5103 */ 5104 PetscErrorCode SNESMonitorSetMatlab(SNES snes,const char *f,mxArray *ctx) 5105 { 5106 PetscErrorCode ierr; 5107 SNESMatlabContext *sctx; 5108 5109 PetscFunctionBegin; 5110 /* currently sctx is memory bleed */ 5111 ierr = PetscMalloc(sizeof(SNESMatlabContext),&sctx);CHKERRQ(ierr); 5112 ierr = PetscStrallocpy(f,&sctx->funcname);CHKERRQ(ierr); 5113 /* 5114 This should work, but it doesn't 5115 sctx->ctx = ctx; 5116 mexMakeArrayPersistent(sctx->ctx); 5117 */ 5118 sctx->ctx = mxDuplicateArray(ctx); 5119 ierr = SNESMonitorSet(snes,SNESMonitor_Matlab,sctx,NULL);CHKERRQ(ierr); 5120 PetscFunctionReturn(0); 5121 } 5122 5123 #endif 5124