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