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,B,C,D,jacobian; 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,test = PETSC_FALSE,flg; 2289 PetscViewer viewer,mviewer; 2290 MPI_Comm comm; 2291 PetscInt tabs; 2292 static PetscBool directionsprinted = PETSC_FALSE; 2293 PetscViewerFormat format; 2294 2295 PetscFunctionBegin; 2296 ierr = PetscObjectOptionsBegin((PetscObject)snes); 2297 ierr = PetscOptionsName("-snes_test_jacobian","Compare hand-coded and finite difference Jacobians","None",&test);CHKERRQ(ierr); 2298 ierr = PetscOptionsReal("-snes_test_jacobian", "Threshold for element difference between hand-coded and finite difference being meaningful", "None", threshold, &threshold,NULL);CHKERRQ(ierr); 2299 ierr = PetscOptionsViewer("-snes_test_jacobian_view","View difference between hand-coded and finite difference Jacobians element entries","None",&mviewer,&format,&complete_print);CHKERRQ(ierr); 2300 ierr = PetscOptionsEnd();CHKERRQ(ierr); 2301 if (!test) PetscFunctionReturn(0); 2302 2303 ierr = PetscObjectGetComm((PetscObject)snes,&comm);CHKERRQ(ierr); 2304 ierr = PetscViewerASCIIGetStdout(comm,&viewer);CHKERRQ(ierr); 2305 ierr = PetscViewerASCIIGetTab(viewer, &tabs);CHKERRQ(ierr); 2306 ierr = PetscViewerASCIISetTab(viewer, ((PetscObject)snes)->tablevel);CHKERRQ(ierr); 2307 ierr = PetscViewerASCIIPrintf(viewer," ---------- Testing Jacobian -------------\n");CHKERRQ(ierr); 2308 if (!complete_print && !directionsprinted) { 2309 ierr = PetscViewerASCIIPrintf(viewer," Run with -snes_test_jacobian_view and optionally -snes_test_jacobian <threshold> to show difference\n");CHKERRQ(ierr); 2310 ierr = PetscViewerASCIIPrintf(viewer," of hand-coded and finite difference Jacobian entries great than <threshold>.\n");CHKERRQ(ierr); 2311 } 2312 if (!directionsprinted) { 2313 ierr = PetscViewerASCIIPrintf(viewer," Testing hand-coded Jacobian, if (for double precision runs) ||J - Jfd||_F/||J||_F is\n");CHKERRQ(ierr); 2314 ierr = PetscViewerASCIIPrintf(viewer," O(1.e-8), the hand-coded Jacobian is probably correct.\n");CHKERRQ(ierr); 2315 directionsprinted = PETSC_TRUE; 2316 } 2317 if (complete_print) { 2318 ierr = PetscViewerPushFormat(mviewer,format);CHKERRQ(ierr); 2319 } 2320 2321 /* evaluate the function at this point because SNESComputeJacobianDefault() assumes that the function has been evaluated and put into snes->vec_func */ 2322 ierr = SNESComputeFunction(snes,x,f);CHKERRQ(ierr); 2323 2324 ierr = PetscObjectTypeCompare((PetscObject)snes->jacobian,MATMFFD,&flg);CHKERRQ(ierr); 2325 if (!flg) jacobian = snes->jacobian; 2326 else jacobian = snes->jacobian_pre; 2327 2328 while (jacobian) { 2329 ierr = PetscObjectTypeCompareAny((PetscObject)jacobian,&flg,MATSEQAIJ,MATMPIAIJ,MATSEQDENSE,MATMPIDENSE,MATSEQBAIJ,MATMPIBAIJ,MATSEQSBAIJ,MATMPIBAIJ,"");CHKERRQ(ierr); 2330 if (flg) { 2331 A = jacobian; 2332 ierr = PetscObjectReference((PetscObject)A);CHKERRQ(ierr); 2333 } else { 2334 ierr = MatComputeExplicitOperator(jacobian,&A);CHKERRQ(ierr); 2335 } 2336 2337 ierr = MatCreate(PetscObjectComm((PetscObject)A),&B);CHKERRQ(ierr); 2338 ierr = MatGetSize(A,&M,&N);CHKERRQ(ierr); 2339 ierr = MatGetLocalSize(A,&m,&n);CHKERRQ(ierr); 2340 ierr = MatSetSizes(B,m,n,M,N);CHKERRQ(ierr); 2341 ierr = MatSetType(B,((PetscObject)A)->type_name);CHKERRQ(ierr); 2342 ierr = MatSetUp(B);CHKERRQ(ierr); 2343 ierr = MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);CHKERRQ(ierr); 2344 2345 ierr = SNESGetFunction(snes,NULL,NULL,&functx);CHKERRQ(ierr); 2346 ierr = SNESComputeJacobianDefault(snes,x,B,B,functx);CHKERRQ(ierr); 2347 2348 ierr = MatDuplicate(B,MAT_COPY_VALUES,&D);CHKERRQ(ierr); 2349 ierr = MatAYPX(D,-1.0,A,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr); 2350 ierr = MatNorm(D,NORM_FROBENIUS,&nrm);CHKERRQ(ierr); 2351 ierr = MatNorm(A,NORM_FROBENIUS,&gnorm);CHKERRQ(ierr); 2352 ierr = MatDestroy(&D);CHKERRQ(ierr); 2353 if (!gnorm) gnorm = 1; /* just in case */ 2354 ierr = PetscViewerASCIIPrintf(viewer,"||J - Jfd||_F/||J||_F = %g, ||J - Jfd||_F = %g\n",(double)(nrm/gnorm),(double)nrm);CHKERRQ(ierr); 2355 2356 if (complete_print) { 2357 ierr = PetscViewerASCIIPrintf(viewer," Hand-coded Jacobian ----------\n");CHKERRQ(ierr); 2358 ierr = MatView(jacobian,mviewer);CHKERRQ(ierr); 2359 ierr = PetscViewerASCIIPrintf(viewer," Finite difference Jacobian ----------\n");CHKERRQ(ierr); 2360 ierr = MatView(B,mviewer);CHKERRQ(ierr); 2361 } 2362 2363 if (complete_print) { 2364 PetscInt Istart, Iend, *ccols, bncols, cncols, j, row; 2365 PetscScalar *cvals; 2366 const PetscInt *bcols; 2367 const PetscScalar *bvals; 2368 2369 ierr = MatAYPX(B,-1.0,A,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr); 2370 ierr = MatCreate(PetscObjectComm((PetscObject)A),&C);CHKERRQ(ierr); 2371 ierr = MatSetSizes(C,m,n,M,N);CHKERRQ(ierr); 2372 ierr = MatSetType(C,((PetscObject)A)->type_name);CHKERRQ(ierr); 2373 ierr = MatSetUp(C);CHKERRQ(ierr); 2374 ierr = MatSetOption(C,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);CHKERRQ(ierr); 2375 ierr = MatGetOwnershipRange(B,&Istart,&Iend);CHKERRQ(ierr); 2376 2377 for (row = Istart; row < Iend; row++) { 2378 ierr = MatGetRow(B,row,&bncols,&bcols,&bvals);CHKERRQ(ierr); 2379 ierr = PetscMalloc2(bncols,&ccols,bncols,&cvals);CHKERRQ(ierr); 2380 for (j = 0, cncols = 0; j < bncols; j++) { 2381 if (PetscAbsScalar(bvals[j]) > threshold) { 2382 ccols[cncols] = bcols[j]; 2383 cvals[cncols] = bvals[j]; 2384 cncols += 1; 2385 } 2386 } 2387 if (cncols) { 2388 ierr = MatSetValues(C,1,&row,cncols,ccols,cvals,INSERT_VALUES);CHKERRQ(ierr); 2389 } 2390 ierr = MatRestoreRow(B,row,&bncols,&bcols,&bvals);CHKERRQ(ierr); 2391 ierr = PetscFree2(ccols,cvals);CHKERRQ(ierr); 2392 } 2393 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2394 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2395 ierr = PetscViewerASCIIPrintf(viewer," Hand-coded minus finite-difference Jacobian with tolerance %g ----------\n",(double)threshold);CHKERRQ(ierr); 2396 ierr = MatView(C,mviewer);CHKERRQ(ierr); 2397 ierr = MatDestroy(&C);CHKERRQ(ierr); 2398 } 2399 ierr = MatDestroy(&A);CHKERRQ(ierr); 2400 ierr = MatDestroy(&B);CHKERRQ(ierr); 2401 2402 if (jacobian != snes->jacobian_pre) { 2403 jacobian = snes->jacobian_pre; 2404 ierr = PetscViewerASCIIPrintf(viewer," ---------- Testing Jacobian for preconditioner -------------\n");CHKERRQ(ierr); 2405 } 2406 else jacobian = NULL; 2407 } 2408 if (complete_print) { 2409 ierr = PetscViewerPopFormat(mviewer);CHKERRQ(ierr); 2410 } 2411 ierr = PetscViewerASCIISetTab(viewer,tabs);CHKERRQ(ierr); 2412 PetscFunctionReturn(0); 2413 } 2414 2415 /*@ 2416 SNESComputeJacobian - Computes the Jacobian matrix that has been set with SNESSetJacobian(). 2417 2418 Collective on SNES and Mat 2419 2420 Input Parameters: 2421 + snes - the SNES context 2422 - x - input vector 2423 2424 Output Parameters: 2425 + A - Jacobian matrix 2426 - B - optional preconditioning matrix 2427 2428 Options Database Keys: 2429 + -snes_lag_preconditioner <lag> 2430 . -snes_lag_jacobian <lag> 2431 . -snes_test_jacobian - compare the user provided Jacobian with one compute via finite differences to check for errors 2432 . -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 2433 . -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 2434 . -snes_compare_explicit - Compare the computed Jacobian to the finite difference Jacobian and output the differences 2435 . -snes_compare_explicit_draw - Compare the computed Jacobian to the finite difference Jacobian and draw the result 2436 . -snes_compare_explicit_contour - Compare the computed Jacobian to the finite difference Jacobian and draw a contour plot with the result 2437 . -snes_compare_operator - Make the comparison options above use the operator instead of the preconditioning matrix 2438 . -snes_compare_coloring - Compute the finite difference Jacobian using coloring and display norms of difference 2439 . -snes_compare_coloring_display - Compute the finite differece Jacobian using coloring and display verbose differences 2440 . -snes_compare_coloring_threshold - Display only those matrix entries that differ by more than a given threshold 2441 . -snes_compare_coloring_threshold_atol - Absolute tolerance for difference in matrix entries to be displayed by -snes_compare_coloring_threshold 2442 . -snes_compare_coloring_threshold_rtol - Relative tolerance for difference in matrix entries to be displayed by -snes_compare_coloring_threshold 2443 . -snes_compare_coloring_draw - Compute the finite differece Jacobian using coloring and draw differences 2444 - -snes_compare_coloring_draw_contour - Compute the finite differece Jacobian using coloring and show contours of matrices and differences 2445 2446 2447 Notes: 2448 Most users should not need to explicitly call this routine, as it 2449 is used internally within the nonlinear solvers. 2450 2451 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 2452 for with the SNESType of test that has been removed. 2453 2454 Level: developer 2455 2456 .keywords: SNES, compute, Jacobian, matrix 2457 2458 .seealso: SNESSetJacobian(), KSPSetOperators(), MatStructure, SNESSetLagPreconditioner(), SNESSetLagJacobian() 2459 @*/ 2460 PetscErrorCode SNESComputeJacobian(SNES snes,Vec X,Mat A,Mat B) 2461 { 2462 PetscErrorCode ierr; 2463 PetscBool flag; 2464 DM dm; 2465 DMSNES sdm; 2466 KSP ksp; 2467 2468 PetscFunctionBegin; 2469 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2470 PetscValidHeaderSpecific(X,VEC_CLASSID,2); 2471 PetscCheckSameComm(snes,1,X,2); 2472 ierr = VecValidValues(X,2,PETSC_TRUE);CHKERRQ(ierr); 2473 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 2474 ierr = DMGetDMSNES(dm,&sdm);CHKERRQ(ierr); 2475 2476 if (!sdm->ops->computejacobian) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_USER,"Must call SNESSetJacobian(), DMSNESSetJacobian(), DMDASNESSetJacobianLocal(), etc"); 2477 2478 /* make sure that MatAssemblyBegin/End() is called on A matrix if it is matrix free */ 2479 2480 if (snes->lagjacobian == -2) { 2481 snes->lagjacobian = -1; 2482 2483 ierr = PetscInfo(snes,"Recomputing Jacobian/preconditioner because lag is -2 (means compute Jacobian, but then never again) \n");CHKERRQ(ierr); 2484 } else if (snes->lagjacobian == -1) { 2485 ierr = PetscInfo(snes,"Reusing Jacobian/preconditioner because lag is -1\n");CHKERRQ(ierr); 2486 ierr = PetscObjectTypeCompare((PetscObject)A,MATMFFD,&flag);CHKERRQ(ierr); 2487 if (flag) { 2488 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2489 ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2490 } 2491 PetscFunctionReturn(0); 2492 } else if (snes->lagjacobian > 1 && (snes->iter + snes->jac_iter) % snes->lagjacobian) { 2493 ierr = PetscInfo2(snes,"Reusing Jacobian/preconditioner because lag is %D and SNES iteration is %D\n",snes->lagjacobian,snes->iter);CHKERRQ(ierr); 2494 ierr = PetscObjectTypeCompare((PetscObject)A,MATMFFD,&flag);CHKERRQ(ierr); 2495 if (flag) { 2496 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2497 ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2498 } 2499 PetscFunctionReturn(0); 2500 } 2501 if (snes->npc && snes->npcside== PC_LEFT) { 2502 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2503 ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2504 PetscFunctionReturn(0); 2505 } 2506 2507 ierr = PetscLogEventBegin(SNES_JacobianEval,snes,X,A,B);CHKERRQ(ierr); 2508 ierr = VecLockPush(X);CHKERRQ(ierr); 2509 PetscStackPush("SNES user Jacobian function"); 2510 ierr = (*sdm->ops->computejacobian)(snes,X,A,B,sdm->jacobianctx);CHKERRQ(ierr); 2511 PetscStackPop; 2512 ierr = VecLockPop(X);CHKERRQ(ierr); 2513 ierr = PetscLogEventEnd(SNES_JacobianEval,snes,X,A,B);CHKERRQ(ierr); 2514 2515 /* the next line ensures that snes->ksp exists */ 2516 ierr = SNESGetKSP(snes,&ksp);CHKERRQ(ierr); 2517 if (snes->lagpreconditioner == -2) { 2518 ierr = PetscInfo(snes,"Rebuilding preconditioner exactly once since lag is -2\n");CHKERRQ(ierr); 2519 ierr = KSPSetReusePreconditioner(snes->ksp,PETSC_FALSE);CHKERRQ(ierr); 2520 snes->lagpreconditioner = -1; 2521 } else if (snes->lagpreconditioner == -1) { 2522 ierr = PetscInfo(snes,"Reusing preconditioner because lag is -1\n");CHKERRQ(ierr); 2523 ierr = KSPSetReusePreconditioner(snes->ksp,PETSC_TRUE);CHKERRQ(ierr); 2524 } else if (snes->lagpreconditioner > 1 && (snes->iter + snes->pre_iter) % snes->lagpreconditioner) { 2525 ierr = PetscInfo2(snes,"Reusing preconditioner because lag is %D and SNES iteration is %D\n",snes->lagpreconditioner,snes->iter);CHKERRQ(ierr); 2526 ierr = KSPSetReusePreconditioner(snes->ksp,PETSC_TRUE);CHKERRQ(ierr); 2527 } else { 2528 ierr = PetscInfo(snes,"Rebuilding preconditioner\n");CHKERRQ(ierr); 2529 ierr = KSPSetReusePreconditioner(snes->ksp,PETSC_FALSE);CHKERRQ(ierr); 2530 } 2531 2532 ierr = SNESTestJacobian(snes);CHKERRQ(ierr); 2533 /* make sure user returned a correct Jacobian and preconditioner */ 2534 /* PetscValidHeaderSpecific(A,MAT_CLASSID,3); 2535 PetscValidHeaderSpecific(B,MAT_CLASSID,4); */ 2536 { 2537 PetscBool flag = PETSC_FALSE,flag_draw = PETSC_FALSE,flag_contour = PETSC_FALSE,flag_operator = PETSC_FALSE; 2538 ierr = PetscOptionsGetViewer(PetscObjectComm((PetscObject)snes),((PetscObject)snes)->prefix,"-snes_compare_explicit",NULL,NULL,&flag);CHKERRQ(ierr); 2539 ierr = PetscOptionsGetViewer(PetscObjectComm((PetscObject)snes),((PetscObject)snes)->prefix,"-snes_compare_explicit_draw",NULL,NULL,&flag_draw);CHKERRQ(ierr); 2540 ierr = PetscOptionsGetViewer(PetscObjectComm((PetscObject)snes),((PetscObject)snes)->prefix,"-snes_compare_explicit_draw_contour",NULL,NULL,&flag_contour);CHKERRQ(ierr); 2541 ierr = PetscOptionsGetViewer(PetscObjectComm((PetscObject)snes),((PetscObject)snes)->prefix,"-snes_compare_operator",NULL,NULL,&flag_operator);CHKERRQ(ierr); 2542 if (flag || flag_draw || flag_contour) { 2543 Mat Bexp_mine = NULL,Bexp,FDexp; 2544 PetscViewer vdraw,vstdout; 2545 PetscBool flg; 2546 if (flag_operator) { 2547 ierr = MatComputeExplicitOperator(A,&Bexp_mine);CHKERRQ(ierr); 2548 Bexp = Bexp_mine; 2549 } else { 2550 /* See if the preconditioning matrix can be viewed and added directly */ 2551 ierr = PetscObjectTypeCompareAny((PetscObject)B,&flg,MATSEQAIJ,MATMPIAIJ,MATSEQDENSE,MATMPIDENSE,MATSEQBAIJ,MATMPIBAIJ,MATSEQSBAIJ,MATMPIBAIJ,"");CHKERRQ(ierr); 2552 if (flg) Bexp = B; 2553 else { 2554 /* If the "preconditioning" matrix is itself MATSHELL or some other type without direct support */ 2555 ierr = MatComputeExplicitOperator(B,&Bexp_mine);CHKERRQ(ierr); 2556 Bexp = Bexp_mine; 2557 } 2558 } 2559 ierr = MatConvert(Bexp,MATSAME,MAT_INITIAL_MATRIX,&FDexp);CHKERRQ(ierr); 2560 ierr = SNESComputeJacobianDefault(snes,X,FDexp,FDexp,NULL);CHKERRQ(ierr); 2561 ierr = PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)snes),&vstdout);CHKERRQ(ierr); 2562 if (flag_draw || flag_contour) { 2563 ierr = PetscViewerDrawOpen(PetscObjectComm((PetscObject)snes),0,"Explicit Jacobians",PETSC_DECIDE,PETSC_DECIDE,300,300,&vdraw);CHKERRQ(ierr); 2564 if (flag_contour) {ierr = PetscViewerPushFormat(vdraw,PETSC_VIEWER_DRAW_CONTOUR);CHKERRQ(ierr);} 2565 } else vdraw = NULL; 2566 ierr = PetscViewerASCIIPrintf(vstdout,"Explicit %s\n",flag_operator ? "Jacobian" : "preconditioning Jacobian");CHKERRQ(ierr); 2567 if (flag) {ierr = MatView(Bexp,vstdout);CHKERRQ(ierr);} 2568 if (vdraw) {ierr = MatView(Bexp,vdraw);CHKERRQ(ierr);} 2569 ierr = PetscViewerASCIIPrintf(vstdout,"Finite difference Jacobian\n");CHKERRQ(ierr); 2570 if (flag) {ierr = MatView(FDexp,vstdout);CHKERRQ(ierr);} 2571 if (vdraw) {ierr = MatView(FDexp,vdraw);CHKERRQ(ierr);} 2572 ierr = MatAYPX(FDexp,-1.0,Bexp,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 2573 ierr = PetscViewerASCIIPrintf(vstdout,"User-provided matrix minus finite difference Jacobian\n");CHKERRQ(ierr); 2574 if (flag) {ierr = MatView(FDexp,vstdout);CHKERRQ(ierr);} 2575 if (vdraw) { /* Always use contour for the difference */ 2576 ierr = PetscViewerPushFormat(vdraw,PETSC_VIEWER_DRAW_CONTOUR);CHKERRQ(ierr); 2577 ierr = MatView(FDexp,vdraw);CHKERRQ(ierr); 2578 ierr = PetscViewerPopFormat(vdraw);CHKERRQ(ierr); 2579 } 2580 if (flag_contour) {ierr = PetscViewerPopFormat(vdraw);CHKERRQ(ierr);} 2581 ierr = PetscViewerDestroy(&vdraw);CHKERRQ(ierr); 2582 ierr = MatDestroy(&Bexp_mine);CHKERRQ(ierr); 2583 ierr = MatDestroy(&FDexp);CHKERRQ(ierr); 2584 } 2585 } 2586 { 2587 PetscBool flag = PETSC_FALSE,flag_display = PETSC_FALSE,flag_draw = PETSC_FALSE,flag_contour = PETSC_FALSE,flag_threshold = PETSC_FALSE; 2588 PetscReal threshold_atol = PETSC_SQRT_MACHINE_EPSILON,threshold_rtol = 10*PETSC_SQRT_MACHINE_EPSILON; 2589 ierr = PetscOptionsGetViewer(PetscObjectComm((PetscObject)snes),((PetscObject)snes)->prefix,"-snes_compare_coloring",NULL,NULL,&flag);CHKERRQ(ierr); 2590 ierr = PetscOptionsGetViewer(PetscObjectComm((PetscObject)snes),((PetscObject)snes)->prefix,"-snes_compare_coloring_display",NULL,NULL,&flag_display);CHKERRQ(ierr); 2591 ierr = PetscOptionsGetViewer(PetscObjectComm((PetscObject)snes),((PetscObject)snes)->prefix,"-snes_compare_coloring_draw",NULL,NULL,&flag_draw);CHKERRQ(ierr); 2592 ierr = PetscOptionsGetViewer(PetscObjectComm((PetscObject)snes),((PetscObject)snes)->prefix,"-snes_compare_coloring_draw_contour",NULL,NULL,&flag_contour);CHKERRQ(ierr); 2593 ierr = PetscOptionsGetViewer(PetscObjectComm((PetscObject)snes),((PetscObject)snes)->prefix,"-snes_compare_coloring_threshold",NULL,NULL,&flag_threshold);CHKERRQ(ierr); 2594 if (flag_threshold) { 2595 ierr = PetscOptionsGetReal(((PetscObject)snes)->options,((PetscObject)snes)->prefix,"-snes_compare_coloring_threshold_rtol",&threshold_rtol,NULL);CHKERRQ(ierr); 2596 ierr = PetscOptionsGetReal(((PetscObject)snes)->options,((PetscObject)snes)->prefix,"-snes_compare_coloring_threshold_atol",&threshold_atol,NULL);CHKERRQ(ierr); 2597 } 2598 if (flag || flag_display || flag_draw || flag_contour || flag_threshold) { 2599 Mat Bfd; 2600 PetscViewer vdraw,vstdout; 2601 MatColoring coloring; 2602 ISColoring iscoloring; 2603 MatFDColoring matfdcoloring; 2604 PetscErrorCode (*func)(SNES,Vec,Vec,void*); 2605 void *funcctx; 2606 PetscReal norm1,norm2,normmax; 2607 2608 ierr = MatDuplicate(B,MAT_DO_NOT_COPY_VALUES,&Bfd);CHKERRQ(ierr); 2609 ierr = MatColoringCreate(Bfd,&coloring);CHKERRQ(ierr); 2610 ierr = MatColoringSetType(coloring,MATCOLORINGSL);CHKERRQ(ierr); 2611 ierr = MatColoringSetFromOptions(coloring);CHKERRQ(ierr); 2612 ierr = MatColoringApply(coloring,&iscoloring);CHKERRQ(ierr); 2613 ierr = MatColoringDestroy(&coloring);CHKERRQ(ierr); 2614 ierr = MatFDColoringCreate(Bfd,iscoloring,&matfdcoloring);CHKERRQ(ierr); 2615 ierr = MatFDColoringSetFromOptions(matfdcoloring);CHKERRQ(ierr); 2616 ierr = MatFDColoringSetUp(Bfd,iscoloring,matfdcoloring);CHKERRQ(ierr); 2617 ierr = ISColoringDestroy(&iscoloring);CHKERRQ(ierr); 2618 2619 /* This method of getting the function is currently unreliable since it doesn't work for DM local functions. */ 2620 ierr = SNESGetFunction(snes,NULL,&func,&funcctx);CHKERRQ(ierr); 2621 ierr = MatFDColoringSetFunction(matfdcoloring,(PetscErrorCode (*)(void))func,funcctx);CHKERRQ(ierr); 2622 ierr = PetscObjectSetOptionsPrefix((PetscObject)matfdcoloring,((PetscObject)snes)->prefix);CHKERRQ(ierr); 2623 ierr = PetscObjectAppendOptionsPrefix((PetscObject)matfdcoloring,"coloring_");CHKERRQ(ierr); 2624 ierr = MatFDColoringSetFromOptions(matfdcoloring);CHKERRQ(ierr); 2625 ierr = MatFDColoringApply(Bfd,matfdcoloring,X,snes);CHKERRQ(ierr); 2626 ierr = MatFDColoringDestroy(&matfdcoloring);CHKERRQ(ierr); 2627 2628 ierr = PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)snes),&vstdout);CHKERRQ(ierr); 2629 if (flag_draw || flag_contour) { 2630 ierr = PetscViewerDrawOpen(PetscObjectComm((PetscObject)snes),0,"Colored Jacobians",PETSC_DECIDE,PETSC_DECIDE,300,300,&vdraw);CHKERRQ(ierr); 2631 if (flag_contour) {ierr = PetscViewerPushFormat(vdraw,PETSC_VIEWER_DRAW_CONTOUR);CHKERRQ(ierr);} 2632 } else vdraw = NULL; 2633 ierr = PetscViewerASCIIPrintf(vstdout,"Explicit preconditioning Jacobian\n");CHKERRQ(ierr); 2634 if (flag_display) {ierr = MatView(B,vstdout);CHKERRQ(ierr);} 2635 if (vdraw) {ierr = MatView(B,vdraw);CHKERRQ(ierr);} 2636 ierr = PetscViewerASCIIPrintf(vstdout,"Colored Finite difference Jacobian\n");CHKERRQ(ierr); 2637 if (flag_display) {ierr = MatView(Bfd,vstdout);CHKERRQ(ierr);} 2638 if (vdraw) {ierr = MatView(Bfd,vdraw);CHKERRQ(ierr);} 2639 ierr = MatAYPX(Bfd,-1.0,B,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 2640 ierr = MatNorm(Bfd,NORM_1,&norm1);CHKERRQ(ierr); 2641 ierr = MatNorm(Bfd,NORM_FROBENIUS,&norm2);CHKERRQ(ierr); 2642 ierr = MatNorm(Bfd,NORM_MAX,&normmax);CHKERRQ(ierr); 2643 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); 2644 if (flag_display) {ierr = MatView(Bfd,vstdout);CHKERRQ(ierr);} 2645 if (vdraw) { /* Always use contour for the difference */ 2646 ierr = PetscViewerPushFormat(vdraw,PETSC_VIEWER_DRAW_CONTOUR);CHKERRQ(ierr); 2647 ierr = MatView(Bfd,vdraw);CHKERRQ(ierr); 2648 ierr = PetscViewerPopFormat(vdraw);CHKERRQ(ierr); 2649 } 2650 if (flag_contour) {ierr = PetscViewerPopFormat(vdraw);CHKERRQ(ierr);} 2651 2652 if (flag_threshold) { 2653 PetscInt bs,rstart,rend,i; 2654 ierr = MatGetBlockSize(B,&bs);CHKERRQ(ierr); 2655 ierr = MatGetOwnershipRange(B,&rstart,&rend);CHKERRQ(ierr); 2656 for (i=rstart; i<rend; i++) { 2657 const PetscScalar *ba,*ca; 2658 const PetscInt *bj,*cj; 2659 PetscInt bn,cn,j,maxentrycol = -1,maxdiffcol = -1,maxrdiffcol = -1; 2660 PetscReal maxentry = 0,maxdiff = 0,maxrdiff = 0; 2661 ierr = MatGetRow(B,i,&bn,&bj,&ba);CHKERRQ(ierr); 2662 ierr = MatGetRow(Bfd,i,&cn,&cj,&ca);CHKERRQ(ierr); 2663 if (bn != cn) SETERRQ(((PetscObject)A)->comm,PETSC_ERR_PLIB,"Unexpected different nonzero pattern in -snes_compare_coloring_threshold"); 2664 for (j=0; j<bn; j++) { 2665 PetscReal rdiff = PetscAbsScalar(ca[j]) / (threshold_atol + threshold_rtol*PetscAbsScalar(ba[j])); 2666 if (PetscAbsScalar(ba[j]) > PetscAbs(maxentry)) { 2667 maxentrycol = bj[j]; 2668 maxentry = PetscRealPart(ba[j]); 2669 } 2670 if (PetscAbsScalar(ca[j]) > PetscAbs(maxdiff)) { 2671 maxdiffcol = bj[j]; 2672 maxdiff = PetscRealPart(ca[j]); 2673 } 2674 if (rdiff > maxrdiff) { 2675 maxrdiffcol = bj[j]; 2676 maxrdiff = rdiff; 2677 } 2678 } 2679 if (maxrdiff > 1) { 2680 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); 2681 for (j=0; j<bn; j++) { 2682 PetscReal rdiff; 2683 rdiff = PetscAbsScalar(ca[j]) / (threshold_atol + threshold_rtol*PetscAbsScalar(ba[j])); 2684 if (rdiff > 1) { 2685 ierr = PetscViewerASCIIPrintf(vstdout," (%D,%g:%g)",bj[j],(double)PetscRealPart(ba[j]),(double)PetscRealPart(ca[j]));CHKERRQ(ierr); 2686 } 2687 } 2688 ierr = PetscViewerASCIIPrintf(vstdout,"\n",i,maxentry,maxdiff,maxrdiff);CHKERRQ(ierr); 2689 } 2690 ierr = MatRestoreRow(B,i,&bn,&bj,&ba);CHKERRQ(ierr); 2691 ierr = MatRestoreRow(Bfd,i,&cn,&cj,&ca);CHKERRQ(ierr); 2692 } 2693 } 2694 ierr = PetscViewerDestroy(&vdraw);CHKERRQ(ierr); 2695 ierr = MatDestroy(&Bfd);CHKERRQ(ierr); 2696 } 2697 } 2698 PetscFunctionReturn(0); 2699 } 2700 2701 /*MC 2702 SNESJacobianFunction - Function used to convey the nonlinear Jacobian of the function to be solved by SNES 2703 2704 Synopsis: 2705 #include "petscsnes.h" 2706 PetscErrorCode SNESJacobianFunction(SNES snes,Vec x,Mat Amat,Mat Pmat,void *ctx); 2707 2708 + x - input vector 2709 . Amat - the matrix that defines the (approximate) Jacobian 2710 . Pmat - the matrix to be used in constructing the preconditioner, usually the same as Amat. 2711 - ctx - [optional] user-defined Jacobian context 2712 2713 Level: intermediate 2714 2715 .seealso: SNESSetFunction(), SNESGetFunction(), SNESSetJacobian(), SNESGetJacobian() 2716 M*/ 2717 2718 /*@C 2719 SNESSetJacobian - Sets the function to compute Jacobian as well as the 2720 location to store the matrix. 2721 2722 Logically Collective on SNES and Mat 2723 2724 Input Parameters: 2725 + snes - the SNES context 2726 . Amat - the matrix that defines the (approximate) Jacobian 2727 . Pmat - the matrix to be used in constructing the preconditioner, usually the same as Amat. 2728 . J - Jacobian evaluation routine (if NULL then SNES retains any previously set value), see SNESJacobianFunction for details 2729 - ctx - [optional] user-defined context for private data for the 2730 Jacobian evaluation routine (may be NULL) (if NULL then SNES retains any previously set value) 2731 2732 Notes: 2733 If the Amat matrix and Pmat matrix are different you must call MatAssemblyBegin/End() on 2734 each matrix. 2735 2736 If you know the operator Amat has a null space you can use MatSetNullSpace() and MatSetTransposeNullSpace() to supply the null 2737 space to Amat and the KSP solvers will automatically use that null space as needed during the solution process. 2738 2739 If using SNESComputeJacobianDefaultColor() to assemble a Jacobian, the ctx argument 2740 must be a MatFDColoring. 2741 2742 Other defect-correction schemes can be used by computing a different matrix in place of the Jacobian. One common 2743 example is to use the "Picard linearization" which only differentiates through the highest order parts of each term. 2744 2745 Level: beginner 2746 2747 .keywords: SNES, nonlinear, set, Jacobian, matrix 2748 2749 .seealso: KSPSetOperators(), SNESSetFunction(), MatMFFDComputeJacobian(), SNESComputeJacobianDefaultColor(), MatStructure, J, 2750 SNESSetPicard(), SNESJacobianFunction 2751 @*/ 2752 PetscErrorCode SNESSetJacobian(SNES snes,Mat Amat,Mat Pmat,PetscErrorCode (*J)(SNES,Vec,Mat,Mat,void*),void *ctx) 2753 { 2754 PetscErrorCode ierr; 2755 DM dm; 2756 2757 PetscFunctionBegin; 2758 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2759 if (Amat) PetscValidHeaderSpecific(Amat,MAT_CLASSID,2); 2760 if (Pmat) PetscValidHeaderSpecific(Pmat,MAT_CLASSID,3); 2761 if (Amat) PetscCheckSameComm(snes,1,Amat,2); 2762 if (Pmat) PetscCheckSameComm(snes,1,Pmat,3); 2763 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 2764 ierr = DMSNESSetJacobian(dm,J,ctx);CHKERRQ(ierr); 2765 if (Amat) { 2766 ierr = PetscObjectReference((PetscObject)Amat);CHKERRQ(ierr); 2767 ierr = MatDestroy(&snes->jacobian);CHKERRQ(ierr); 2768 2769 snes->jacobian = Amat; 2770 } 2771 if (Pmat) { 2772 ierr = PetscObjectReference((PetscObject)Pmat);CHKERRQ(ierr); 2773 ierr = MatDestroy(&snes->jacobian_pre);CHKERRQ(ierr); 2774 2775 snes->jacobian_pre = Pmat; 2776 } 2777 PetscFunctionReturn(0); 2778 } 2779 2780 /*@C 2781 SNESGetJacobian - Returns the Jacobian matrix and optionally the user 2782 provided context for evaluating the Jacobian. 2783 2784 Not Collective, but Mat object will be parallel if SNES object is 2785 2786 Input Parameter: 2787 . snes - the nonlinear solver context 2788 2789 Output Parameters: 2790 + Amat - location to stash (approximate) Jacobian matrix (or NULL) 2791 . Pmat - location to stash matrix used to compute the preconditioner (or NULL) 2792 . J - location to put Jacobian function (or NULL), see SNESJacobianFunction for details on its calling sequence 2793 - ctx - location to stash Jacobian ctx (or NULL) 2794 2795 Level: advanced 2796 2797 .seealso: SNESSetJacobian(), SNESComputeJacobian(), SNESJacobianFunction, SNESGetFunction() 2798 @*/ 2799 PetscErrorCode SNESGetJacobian(SNES snes,Mat *Amat,Mat *Pmat,PetscErrorCode (**J)(SNES,Vec,Mat,Mat,void*),void **ctx) 2800 { 2801 PetscErrorCode ierr; 2802 DM dm; 2803 DMSNES sdm; 2804 2805 PetscFunctionBegin; 2806 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2807 if (Amat) *Amat = snes->jacobian; 2808 if (Pmat) *Pmat = snes->jacobian_pre; 2809 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 2810 ierr = DMGetDMSNES(dm,&sdm);CHKERRQ(ierr); 2811 if (J) *J = sdm->ops->computejacobian; 2812 if (ctx) *ctx = sdm->jacobianctx; 2813 PetscFunctionReturn(0); 2814 } 2815 2816 /*@ 2817 SNESSetUp - Sets up the internal data structures for the later use 2818 of a nonlinear solver. 2819 2820 Collective on SNES 2821 2822 Input Parameters: 2823 . snes - the SNES context 2824 2825 Notes: 2826 For basic use of the SNES solvers the user need not explicitly call 2827 SNESSetUp(), since these actions will automatically occur during 2828 the call to SNESSolve(). However, if one wishes to control this 2829 phase separately, SNESSetUp() should be called after SNESCreate() 2830 and optional routines of the form SNESSetXXX(), but before SNESSolve(). 2831 2832 Level: advanced 2833 2834 .keywords: SNES, nonlinear, setup 2835 2836 .seealso: SNESCreate(), SNESSolve(), SNESDestroy() 2837 @*/ 2838 PetscErrorCode SNESSetUp(SNES snes) 2839 { 2840 PetscErrorCode ierr; 2841 DM dm; 2842 DMSNES sdm; 2843 SNESLineSearch linesearch, pclinesearch; 2844 void *lsprectx,*lspostctx; 2845 PetscErrorCode (*precheck)(SNESLineSearch,Vec,Vec,PetscBool*,void*); 2846 PetscErrorCode (*postcheck)(SNESLineSearch,Vec,Vec,Vec,PetscBool*,PetscBool*,void*); 2847 PetscErrorCode (*func)(SNES,Vec,Vec,void*); 2848 Vec f,fpc; 2849 void *funcctx; 2850 PetscErrorCode (*jac)(SNES,Vec,Mat,Mat,void*); 2851 void *jacctx,*appctx; 2852 Mat j,jpre; 2853 2854 PetscFunctionBegin; 2855 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2856 if (snes->setupcalled) PetscFunctionReturn(0); 2857 2858 if (!((PetscObject)snes)->type_name) { 2859 ierr = SNESSetType(snes,SNESNEWTONLS);CHKERRQ(ierr); 2860 } 2861 2862 ierr = SNESGetFunction(snes,&snes->vec_func,NULL,NULL);CHKERRQ(ierr); 2863 2864 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 2865 ierr = DMGetDMSNES(dm,&sdm);CHKERRQ(ierr); 2866 if (!sdm->ops->computefunction) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"Function never provided to SNES object"); 2867 if (!sdm->ops->computejacobian) { 2868 ierr = DMSNESSetJacobian(dm,SNESComputeJacobianDefaultColor,NULL);CHKERRQ(ierr); 2869 } 2870 if (!snes->vec_func) { 2871 ierr = DMCreateGlobalVector(dm,&snes->vec_func);CHKERRQ(ierr); 2872 } 2873 2874 if (!snes->ksp) { 2875 ierr = SNESGetKSP(snes, &snes->ksp);CHKERRQ(ierr); 2876 } 2877 2878 if (!snes->linesearch) { 2879 ierr = SNESGetLineSearch(snes, &snes->linesearch);CHKERRQ(ierr); 2880 } 2881 ierr = SNESLineSearchSetFunction(snes->linesearch,SNESComputeFunction);CHKERRQ(ierr); 2882 2883 if (snes->npc && (snes->npcside== PC_LEFT)) { 2884 snes->mf = PETSC_TRUE; 2885 snes->mf_operator = PETSC_FALSE; 2886 } 2887 2888 if (snes->npc) { 2889 /* copy the DM over */ 2890 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 2891 ierr = SNESSetDM(snes->npc,dm);CHKERRQ(ierr); 2892 2893 ierr = SNESGetFunction(snes,&f,&func,&funcctx);CHKERRQ(ierr); 2894 ierr = VecDuplicate(f,&fpc);CHKERRQ(ierr); 2895 ierr = SNESSetFunction(snes->npc,fpc,func,funcctx);CHKERRQ(ierr); 2896 ierr = SNESGetJacobian(snes,&j,&jpre,&jac,&jacctx);CHKERRQ(ierr); 2897 ierr = SNESSetJacobian(snes->npc,j,jpre,jac,jacctx);CHKERRQ(ierr); 2898 ierr = SNESGetApplicationContext(snes,&appctx);CHKERRQ(ierr); 2899 ierr = SNESSetApplicationContext(snes->npc,appctx);CHKERRQ(ierr); 2900 ierr = VecDestroy(&fpc);CHKERRQ(ierr); 2901 2902 /* copy the function pointers over */ 2903 ierr = PetscObjectCopyFortranFunctionPointers((PetscObject)snes,(PetscObject)snes->npc);CHKERRQ(ierr); 2904 2905 /* default to 1 iteration */ 2906 ierr = SNESSetTolerances(snes->npc,0.0,0.0,0.0,1,snes->npc->max_funcs);CHKERRQ(ierr); 2907 if (snes->npcside==PC_RIGHT) { 2908 ierr = SNESSetNormSchedule(snes->npc,SNES_NORM_FINAL_ONLY);CHKERRQ(ierr); 2909 } else { 2910 ierr = SNESSetNormSchedule(snes->npc,SNES_NORM_NONE);CHKERRQ(ierr); 2911 } 2912 ierr = SNESSetFromOptions(snes->npc);CHKERRQ(ierr); 2913 2914 /* copy the line search context over */ 2915 ierr = SNESGetLineSearch(snes,&linesearch);CHKERRQ(ierr); 2916 ierr = SNESGetLineSearch(snes->npc,&pclinesearch);CHKERRQ(ierr); 2917 ierr = SNESLineSearchGetPreCheck(linesearch,&precheck,&lsprectx);CHKERRQ(ierr); 2918 ierr = SNESLineSearchGetPostCheck(linesearch,&postcheck,&lspostctx);CHKERRQ(ierr); 2919 ierr = SNESLineSearchSetPreCheck(pclinesearch,precheck,lsprectx);CHKERRQ(ierr); 2920 ierr = SNESLineSearchSetPostCheck(pclinesearch,postcheck,lspostctx);CHKERRQ(ierr); 2921 ierr = PetscObjectCopyFortranFunctionPointers((PetscObject)linesearch, (PetscObject)pclinesearch);CHKERRQ(ierr); 2922 } 2923 if (snes->mf) { 2924 ierr = SNESSetUpMatrixFree_Private(snes, snes->mf_operator, snes->mf_version);CHKERRQ(ierr); 2925 } 2926 if (snes->ops->usercompute && !snes->user) { 2927 ierr = (*snes->ops->usercompute)(snes,(void**)&snes->user);CHKERRQ(ierr); 2928 } 2929 2930 snes->jac_iter = 0; 2931 snes->pre_iter = 0; 2932 2933 if (snes->ops->setup) { 2934 ierr = (*snes->ops->setup)(snes);CHKERRQ(ierr); 2935 } 2936 2937 if (snes->npc && (snes->npcside== PC_LEFT)) { 2938 if (snes->functype == SNES_FUNCTION_PRECONDITIONED) { 2939 ierr = SNESGetLineSearch(snes,&linesearch);CHKERRQ(ierr); 2940 ierr = SNESLineSearchSetFunction(linesearch,SNESComputeFunctionDefaultNPC);CHKERRQ(ierr); 2941 } 2942 } 2943 2944 snes->setupcalled = PETSC_TRUE; 2945 PetscFunctionReturn(0); 2946 } 2947 2948 /*@ 2949 SNESReset - Resets a SNES context to the snessetupcalled = 0 state and removes any allocated Vecs and Mats 2950 2951 Collective on SNES 2952 2953 Input Parameter: 2954 . snes - iterative context obtained from SNESCreate() 2955 2956 Level: intermediate 2957 2958 Notes: Also calls the application context destroy routine set with SNESSetComputeApplicationContext() 2959 2960 .keywords: SNES, destroy 2961 2962 .seealso: SNESCreate(), SNESSetUp(), SNESSolve() 2963 @*/ 2964 PetscErrorCode SNESReset(SNES snes) 2965 { 2966 PetscErrorCode ierr; 2967 2968 PetscFunctionBegin; 2969 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2970 if (snes->ops->userdestroy && snes->user) { 2971 ierr = (*snes->ops->userdestroy)((void**)&snes->user);CHKERRQ(ierr); 2972 snes->user = NULL; 2973 } 2974 if (snes->npc) { 2975 ierr = SNESReset(snes->npc);CHKERRQ(ierr); 2976 } 2977 2978 if (snes->ops->reset) { 2979 ierr = (*snes->ops->reset)(snes);CHKERRQ(ierr); 2980 } 2981 if (snes->ksp) { 2982 ierr = KSPReset(snes->ksp);CHKERRQ(ierr); 2983 } 2984 2985 if (snes->linesearch) { 2986 ierr = SNESLineSearchReset(snes->linesearch);CHKERRQ(ierr); 2987 } 2988 2989 ierr = VecDestroy(&snes->vec_rhs);CHKERRQ(ierr); 2990 ierr = VecDestroy(&snes->vec_sol);CHKERRQ(ierr); 2991 ierr = VecDestroy(&snes->vec_sol_update);CHKERRQ(ierr); 2992 ierr = VecDestroy(&snes->vec_func);CHKERRQ(ierr); 2993 ierr = MatDestroy(&snes->jacobian);CHKERRQ(ierr); 2994 ierr = MatDestroy(&snes->jacobian_pre);CHKERRQ(ierr); 2995 ierr = VecDestroyVecs(snes->nwork,&snes->work);CHKERRQ(ierr); 2996 ierr = VecDestroyVecs(snes->nvwork,&snes->vwork);CHKERRQ(ierr); 2997 2998 snes->alwayscomputesfinalresidual = PETSC_FALSE; 2999 3000 snes->nwork = snes->nvwork = 0; 3001 snes->setupcalled = PETSC_FALSE; 3002 PetscFunctionReturn(0); 3003 } 3004 3005 /*@ 3006 SNESDestroy - Destroys the nonlinear solver context that was created 3007 with SNESCreate(). 3008 3009 Collective on SNES 3010 3011 Input Parameter: 3012 . snes - the SNES context 3013 3014 Level: beginner 3015 3016 .keywords: SNES, nonlinear, destroy 3017 3018 .seealso: SNESCreate(), SNESSolve() 3019 @*/ 3020 PetscErrorCode SNESDestroy(SNES *snes) 3021 { 3022 PetscErrorCode ierr; 3023 3024 PetscFunctionBegin; 3025 if (!*snes) PetscFunctionReturn(0); 3026 PetscValidHeaderSpecific((*snes),SNES_CLASSID,1); 3027 if (--((PetscObject)(*snes))->refct > 0) {*snes = 0; PetscFunctionReturn(0);} 3028 3029 ierr = SNESReset((*snes));CHKERRQ(ierr); 3030 ierr = SNESDestroy(&(*snes)->npc);CHKERRQ(ierr); 3031 3032 /* if memory was published with SAWs then destroy it */ 3033 ierr = PetscObjectSAWsViewOff((PetscObject)*snes);CHKERRQ(ierr); 3034 if ((*snes)->ops->destroy) {ierr = (*((*snes))->ops->destroy)((*snes));CHKERRQ(ierr);} 3035 3036 if ((*snes)->dm) {ierr = DMCoarsenHookRemove((*snes)->dm,DMCoarsenHook_SNESVecSol,DMRestrictHook_SNESVecSol,*snes);CHKERRQ(ierr);} 3037 ierr = DMDestroy(&(*snes)->dm);CHKERRQ(ierr); 3038 ierr = KSPDestroy(&(*snes)->ksp);CHKERRQ(ierr); 3039 ierr = SNESLineSearchDestroy(&(*snes)->linesearch);CHKERRQ(ierr); 3040 3041 ierr = PetscFree((*snes)->kspconvctx);CHKERRQ(ierr); 3042 if ((*snes)->ops->convergeddestroy) { 3043 ierr = (*(*snes)->ops->convergeddestroy)((*snes)->cnvP);CHKERRQ(ierr); 3044 } 3045 if ((*snes)->conv_malloc) { 3046 ierr = PetscFree((*snes)->conv_hist);CHKERRQ(ierr); 3047 ierr = PetscFree((*snes)->conv_hist_its);CHKERRQ(ierr); 3048 } 3049 ierr = SNESMonitorCancel((*snes));CHKERRQ(ierr); 3050 ierr = PetscHeaderDestroy(snes);CHKERRQ(ierr); 3051 PetscFunctionReturn(0); 3052 } 3053 3054 /* ----------- Routines to set solver parameters ---------- */ 3055 3056 /*@ 3057 SNESSetLagPreconditioner - Determines when the preconditioner is rebuilt in the nonlinear solve. 3058 3059 Logically Collective on SNES 3060 3061 Input Parameters: 3062 + snes - the SNES context 3063 - lag - -1 indicates NEVER rebuild, 1 means rebuild every time the Jacobian is computed within a single nonlinear solve, 2 means every second time 3064 the Jacobian is built etc. -2 indicates rebuild preconditioner at next chance but then never rebuild after that 3065 3066 Options Database Keys: 3067 . -snes_lag_preconditioner <lag> 3068 3069 Notes: 3070 The default is 1 3071 The preconditioner is ALWAYS built in the first iteration of a nonlinear solve unless lag is -1 3072 If -1 is used before the very first nonlinear solve the preconditioner is still built because there is no previous preconditioner to use 3073 3074 Level: intermediate 3075 3076 .keywords: SNES, nonlinear, set, convergence, tolerances 3077 3078 .seealso: SNESSetTrustRegionTolerance(), SNESGetLagPreconditioner(), SNESSetLagJacobian(), SNESGetLagJacobian() 3079 3080 @*/ 3081 PetscErrorCode SNESSetLagPreconditioner(SNES snes,PetscInt lag) 3082 { 3083 PetscFunctionBegin; 3084 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3085 if (lag < -2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Lag must be -2, -1, 1 or greater"); 3086 if (!lag) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Lag cannot be 0"); 3087 PetscValidLogicalCollectiveInt(snes,lag,2); 3088 snes->lagpreconditioner = lag; 3089 PetscFunctionReturn(0); 3090 } 3091 3092 /*@ 3093 SNESSetGridSequence - sets the number of steps of grid sequencing that SNES does 3094 3095 Logically Collective on SNES 3096 3097 Input Parameters: 3098 + snes - the SNES context 3099 - steps - the number of refinements to do, defaults to 0 3100 3101 Options Database Keys: 3102 . -snes_grid_sequence <steps> 3103 3104 Level: intermediate 3105 3106 Notes: 3107 Use SNESGetSolution() to extract the fine grid solution after grid sequencing. 3108 3109 .keywords: SNES, nonlinear, set, convergence, tolerances 3110 3111 .seealso: SNESSetTrustRegionTolerance(), SNESGetLagPreconditioner(), SNESSetLagJacobian(), SNESGetLagJacobian(), SNESGetGridSequence() 3112 3113 @*/ 3114 PetscErrorCode SNESSetGridSequence(SNES snes,PetscInt steps) 3115 { 3116 PetscFunctionBegin; 3117 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3118 PetscValidLogicalCollectiveInt(snes,steps,2); 3119 snes->gridsequence = steps; 3120 PetscFunctionReturn(0); 3121 } 3122 3123 /*@ 3124 SNESGetGridSequence - gets the number of steps of grid sequencing that SNES does 3125 3126 Logically Collective on SNES 3127 3128 Input Parameter: 3129 . snes - the SNES context 3130 3131 Output Parameter: 3132 . steps - the number of refinements to do, defaults to 0 3133 3134 Options Database Keys: 3135 . -snes_grid_sequence <steps> 3136 3137 Level: intermediate 3138 3139 Notes: 3140 Use SNESGetSolution() to extract the fine grid solution after grid sequencing. 3141 3142 .keywords: SNES, nonlinear, set, convergence, tolerances 3143 3144 .seealso: SNESSetTrustRegionTolerance(), SNESGetLagPreconditioner(), SNESSetLagJacobian(), SNESGetLagJacobian(), SNESSetGridSequence() 3145 3146 @*/ 3147 PetscErrorCode SNESGetGridSequence(SNES snes,PetscInt *steps) 3148 { 3149 PetscFunctionBegin; 3150 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3151 *steps = snes->gridsequence; 3152 PetscFunctionReturn(0); 3153 } 3154 3155 /*@ 3156 SNESGetLagPreconditioner - Indicates how often the preconditioner is rebuilt 3157 3158 Not Collective 3159 3160 Input Parameter: 3161 . snes - the SNES context 3162 3163 Output Parameter: 3164 . lag - -1 indicates NEVER rebuild, 1 means rebuild every time the Jacobian is computed within a single nonlinear solve, 2 means every second time 3165 the Jacobian is built etc. -2 indicates rebuild preconditioner at next chance but then never rebuild after that 3166 3167 Options Database Keys: 3168 . -snes_lag_preconditioner <lag> 3169 3170 Notes: 3171 The default is 1 3172 The preconditioner is ALWAYS built in the first iteration of a nonlinear solve unless lag is -1 3173 3174 Level: intermediate 3175 3176 .keywords: SNES, nonlinear, set, convergence, tolerances 3177 3178 .seealso: SNESSetTrustRegionTolerance(), SNESSetLagPreconditioner() 3179 3180 @*/ 3181 PetscErrorCode SNESGetLagPreconditioner(SNES snes,PetscInt *lag) 3182 { 3183 PetscFunctionBegin; 3184 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3185 *lag = snes->lagpreconditioner; 3186 PetscFunctionReturn(0); 3187 } 3188 3189 /*@ 3190 SNESSetLagJacobian - Determines when the Jacobian is rebuilt in the nonlinear solve. See SNESSetLagPreconditioner() for determining how 3191 often the preconditioner is rebuilt. 3192 3193 Logically Collective on SNES 3194 3195 Input Parameters: 3196 + snes - the SNES context 3197 - lag - -1 indicates NEVER rebuild, 1 means rebuild every time the Jacobian is computed within a single nonlinear solve, 2 means every second time 3198 the Jacobian is built etc. -2 means rebuild at next chance but then never again 3199 3200 Options Database Keys: 3201 . -snes_lag_jacobian <lag> 3202 3203 Notes: 3204 The default is 1 3205 The Jacobian is ALWAYS built in the first iteration of a nonlinear solve unless lag is -1 3206 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 3207 at the next Newton step but never again (unless it is reset to another value) 3208 3209 Level: intermediate 3210 3211 .keywords: SNES, nonlinear, set, convergence, tolerances 3212 3213 .seealso: SNESSetTrustRegionTolerance(), SNESGetLagPreconditioner(), SNESSetLagPreconditioner(), SNESGetLagJacobian() 3214 3215 @*/ 3216 PetscErrorCode SNESSetLagJacobian(SNES snes,PetscInt lag) 3217 { 3218 PetscFunctionBegin; 3219 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3220 if (lag < -2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Lag must be -2, -1, 1 or greater"); 3221 if (!lag) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Lag cannot be 0"); 3222 PetscValidLogicalCollectiveInt(snes,lag,2); 3223 snes->lagjacobian = lag; 3224 PetscFunctionReturn(0); 3225 } 3226 3227 /*@ 3228 SNESGetLagJacobian - Indicates how often the Jacobian is rebuilt. See SNESGetLagPreconditioner() to determine when the preconditioner is rebuilt 3229 3230 Not Collective 3231 3232 Input Parameter: 3233 . snes - the SNES context 3234 3235 Output Parameter: 3236 . lag - -1 indicates NEVER rebuild, 1 means rebuild every time the Jacobian is computed within a single nonlinear solve, 2 means every second time 3237 the Jacobian is built etc. 3238 3239 Options Database Keys: 3240 . -snes_lag_jacobian <lag> 3241 3242 Notes: 3243 The default is 1 3244 The jacobian is ALWAYS built in the first iteration of a nonlinear solve unless lag is -1 3245 3246 Level: intermediate 3247 3248 .keywords: SNES, nonlinear, set, convergence, tolerances 3249 3250 .seealso: SNESSetTrustRegionTolerance(), SNESSetLagJacobian(), SNESSetLagPreconditioner(), SNESGetLagPreconditioner() 3251 3252 @*/ 3253 PetscErrorCode SNESGetLagJacobian(SNES snes,PetscInt *lag) 3254 { 3255 PetscFunctionBegin; 3256 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3257 *lag = snes->lagjacobian; 3258 PetscFunctionReturn(0); 3259 } 3260 3261 /*@ 3262 SNESSetLagJacobianPersists - Set whether or not the Jacobian lagging persists through multiple solves 3263 3264 Logically collective on SNES 3265 3266 Input Parameter: 3267 + snes - the SNES context 3268 - flg - jacobian lagging persists if true 3269 3270 Options Database Keys: 3271 . -snes_lag_jacobian_persists <flg> 3272 3273 Notes: This is useful both for nonlinear preconditioning, where it's appropriate to have the Jacobian be stale by 3274 several solves, and for implicit time-stepping, where Jacobian lagging in the inner nonlinear solve over several 3275 timesteps may present huge efficiency gains. 3276 3277 Level: developer 3278 3279 .keywords: SNES, nonlinear, lag 3280 3281 .seealso: SNESSetLagPreconditionerPersists(), SNESSetLagJacobian(), SNESGetLagJacobian(), SNESGetNPC() 3282 3283 @*/ 3284 PetscErrorCode SNESSetLagJacobianPersists(SNES snes,PetscBool flg) 3285 { 3286 PetscFunctionBegin; 3287 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3288 PetscValidLogicalCollectiveBool(snes,flg,2); 3289 snes->lagjac_persist = flg; 3290 PetscFunctionReturn(0); 3291 } 3292 3293 /*@ 3294 SNESSetLagPreconditionerPersists - Set whether or not the preconditioner lagging persists through multiple solves 3295 3296 Logically Collective on SNES 3297 3298 Input Parameter: 3299 + snes - the SNES context 3300 - flg - preconditioner lagging persists if true 3301 3302 Options Database Keys: 3303 . -snes_lag_jacobian_persists <flg> 3304 3305 Notes: This is useful both for nonlinear preconditioning, where it's appropriate to have the preconditioner be stale 3306 by several solves, and for implicit time-stepping, where preconditioner lagging in the inner nonlinear solve over 3307 several timesteps may present huge efficiency gains. 3308 3309 Level: developer 3310 3311 .keywords: SNES, nonlinear, lag 3312 3313 .seealso: SNESSetLagJacobianPersists(), SNESSetLagJacobian(), SNESGetLagJacobian(), SNESGetNPC() 3314 3315 @*/ 3316 PetscErrorCode SNESSetLagPreconditionerPersists(SNES snes,PetscBool flg) 3317 { 3318 PetscFunctionBegin; 3319 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3320 PetscValidLogicalCollectiveBool(snes,flg,2); 3321 snes->lagpre_persist = flg; 3322 PetscFunctionReturn(0); 3323 } 3324 3325 /*@ 3326 SNESSetForceIteration - force SNESSolve() to take at least one iteration regardless of the initial residual norm 3327 3328 Logically Collective on SNES 3329 3330 Input Parameters: 3331 + snes - the SNES context 3332 - force - PETSC_TRUE require at least one iteration 3333 3334 Options Database Keys: 3335 . -snes_force_iteration <force> - Sets forcing an iteration 3336 3337 Notes: 3338 This is used sometimes with TS to prevent TS from detecting a false steady state solution 3339 3340 Level: intermediate 3341 3342 .keywords: SNES, nonlinear, set, convergence, tolerances 3343 3344 .seealso: SNESSetTrustRegionTolerance(), SNESSetDivergenceTolerance() 3345 @*/ 3346 PetscErrorCode SNESSetForceIteration(SNES snes,PetscBool force) 3347 { 3348 PetscFunctionBegin; 3349 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3350 snes->forceiteration = force; 3351 PetscFunctionReturn(0); 3352 } 3353 3354 /*@ 3355 SNESGetForceIteration - Whether or not to force SNESSolve() take at least one iteration regardless of the initial residual norm 3356 3357 Logically Collective on SNES 3358 3359 Input Parameters: 3360 . snes - the SNES context 3361 3362 Output Parameter: 3363 . force - PETSC_TRUE requires at least one iteration. 3364 3365 .keywords: SNES, nonlinear, set, convergence, tolerances 3366 3367 .seealso: SNESSetForceIteration(), SNESSetTrustRegionTolerance(), SNESSetDivergenceTolerance() 3368 @*/ 3369 PetscErrorCode SNESGetForceIteration(SNES snes,PetscBool *force) 3370 { 3371 PetscFunctionBegin; 3372 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3373 *force = snes->forceiteration; 3374 PetscFunctionReturn(0); 3375 } 3376 3377 /*@ 3378 SNESSetTolerances - Sets various parameters used in convergence tests. 3379 3380 Logically Collective on SNES 3381 3382 Input Parameters: 3383 + snes - the SNES context 3384 . abstol - absolute convergence tolerance 3385 . rtol - relative convergence tolerance 3386 . stol - convergence tolerance in terms of the norm of the change in the solution between steps, || delta x || < stol*|| x || 3387 . maxit - maximum number of iterations 3388 - maxf - maximum number of function evaluations 3389 3390 Options Database Keys: 3391 + -snes_atol <abstol> - Sets abstol 3392 . -snes_rtol <rtol> - Sets rtol 3393 . -snes_stol <stol> - Sets stol 3394 . -snes_max_it <maxit> - Sets maxit 3395 - -snes_max_funcs <maxf> - Sets maxf 3396 3397 Notes: 3398 The default maximum number of iterations is 50. 3399 The default maximum number of function evaluations is 1000. 3400 3401 Level: intermediate 3402 3403 .keywords: SNES, nonlinear, set, convergence, tolerances 3404 3405 .seealso: SNESSetTrustRegionTolerance(), SNESSetDivergenceTolerance(), SNESSetForceIteration() 3406 @*/ 3407 PetscErrorCode SNESSetTolerances(SNES snes,PetscReal abstol,PetscReal rtol,PetscReal stol,PetscInt maxit,PetscInt maxf) 3408 { 3409 PetscFunctionBegin; 3410 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3411 PetscValidLogicalCollectiveReal(snes,abstol,2); 3412 PetscValidLogicalCollectiveReal(snes,rtol,3); 3413 PetscValidLogicalCollectiveReal(snes,stol,4); 3414 PetscValidLogicalCollectiveInt(snes,maxit,5); 3415 PetscValidLogicalCollectiveInt(snes,maxf,6); 3416 3417 if (abstol != PETSC_DEFAULT) { 3418 if (abstol < 0.0) SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_OUTOFRANGE,"Absolute tolerance %g must be non-negative",(double)abstol); 3419 snes->abstol = abstol; 3420 } 3421 if (rtol != PETSC_DEFAULT) { 3422 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); 3423 snes->rtol = rtol; 3424 } 3425 if (stol != PETSC_DEFAULT) { 3426 if (stol < 0.0) SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_OUTOFRANGE,"Step tolerance %g must be non-negative",(double)stol); 3427 snes->stol = stol; 3428 } 3429 if (maxit != PETSC_DEFAULT) { 3430 if (maxit < 0) SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_OUTOFRANGE,"Maximum number of iterations %D must be non-negative",maxit); 3431 snes->max_its = maxit; 3432 } 3433 if (maxf != PETSC_DEFAULT) { 3434 if (maxf < 0) SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_OUTOFRANGE,"Maximum number of function evaluations %D must be non-negative",maxf); 3435 snes->max_funcs = maxf; 3436 } 3437 snes->tolerancesset = PETSC_TRUE; 3438 PetscFunctionReturn(0); 3439 } 3440 3441 /*@ 3442 SNESSetDivergenceTolerance - Sets the divergence tolerance used for the SNES divergence test. 3443 3444 Logically Collective on SNES 3445 3446 Input Parameters: 3447 + snes - the SNES context 3448 - divtol - the divergence tolerance. Use -1 to deactivate the test. 3449 3450 Options Database Keys: 3451 + -snes_divergence_tolerance <divtol> - Sets divtol 3452 3453 Notes: 3454 The default divergence tolerance is 1e4. 3455 3456 Level: intermediate 3457 3458 .keywords: SNES, nonlinear, set, divergence, tolerance 3459 3460 .seealso: SNESSetTolerances(), SNESGetDivergenceTolerance 3461 @*/ 3462 PetscErrorCode SNESSetDivergenceTolerance(SNES snes,PetscReal divtol) 3463 { 3464 PetscFunctionBegin; 3465 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3466 PetscValidLogicalCollectiveReal(snes,divtol,2); 3467 3468 if (divtol != PETSC_DEFAULT) { 3469 snes->divtol = divtol; 3470 } 3471 else { 3472 snes->divtol = 1.0e4; 3473 } 3474 PetscFunctionReturn(0); 3475 } 3476 3477 /*@ 3478 SNESGetTolerances - Gets various parameters used in convergence tests. 3479 3480 Not Collective 3481 3482 Input Parameters: 3483 + snes - the SNES context 3484 . atol - absolute convergence tolerance 3485 . rtol - relative convergence tolerance 3486 . stol - convergence tolerance in terms of the norm 3487 of the change in the solution between steps 3488 . maxit - maximum number of iterations 3489 - maxf - maximum number of function evaluations 3490 3491 Notes: 3492 The user can specify NULL for any parameter that is not needed. 3493 3494 Level: intermediate 3495 3496 .keywords: SNES, nonlinear, get, convergence, tolerances 3497 3498 .seealso: SNESSetTolerances() 3499 @*/ 3500 PetscErrorCode SNESGetTolerances(SNES snes,PetscReal *atol,PetscReal *rtol,PetscReal *stol,PetscInt *maxit,PetscInt *maxf) 3501 { 3502 PetscFunctionBegin; 3503 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3504 if (atol) *atol = snes->abstol; 3505 if (rtol) *rtol = snes->rtol; 3506 if (stol) *stol = snes->stol; 3507 if (maxit) *maxit = snes->max_its; 3508 if (maxf) *maxf = snes->max_funcs; 3509 PetscFunctionReturn(0); 3510 } 3511 3512 /*@ 3513 SNESGetDivergenceTolerance - Gets divergence tolerance used in divergence test. 3514 3515 Not Collective 3516 3517 Input Parameters: 3518 + snes - the SNES context 3519 - divtol - divergence tolerance 3520 3521 Level: intermediate 3522 3523 .keywords: SNES, nonlinear, get, divergence, tolerance 3524 3525 .seealso: SNESSetDivergenceTolerance() 3526 @*/ 3527 PetscErrorCode SNESGetDivergenceTolerance(SNES snes,PetscReal *divtol) 3528 { 3529 PetscFunctionBegin; 3530 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3531 if (divtol) *divtol = snes->divtol; 3532 PetscFunctionReturn(0); 3533 } 3534 3535 /*@ 3536 SNESSetTrustRegionTolerance - Sets the trust region parameter tolerance. 3537 3538 Logically Collective on SNES 3539 3540 Input Parameters: 3541 + snes - the SNES context 3542 - tol - tolerance 3543 3544 Options Database Key: 3545 . -snes_trtol <tol> - Sets tol 3546 3547 Level: intermediate 3548 3549 .keywords: SNES, nonlinear, set, trust region, tolerance 3550 3551 .seealso: SNESSetTolerances() 3552 @*/ 3553 PetscErrorCode SNESSetTrustRegionTolerance(SNES snes,PetscReal tol) 3554 { 3555 PetscFunctionBegin; 3556 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3557 PetscValidLogicalCollectiveReal(snes,tol,2); 3558 snes->deltatol = tol; 3559 PetscFunctionReturn(0); 3560 } 3561 3562 /* 3563 Duplicate the lg monitors for SNES from KSP; for some reason with 3564 dynamic libraries things don't work under Sun4 if we just use 3565 macros instead of functions 3566 */ 3567 PetscErrorCode SNESMonitorLGResidualNorm(SNES snes,PetscInt it,PetscReal norm,void *ctx) 3568 { 3569 PetscErrorCode ierr; 3570 3571 PetscFunctionBegin; 3572 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3573 ierr = KSPMonitorLGResidualNorm((KSP)snes,it,norm,ctx);CHKERRQ(ierr); 3574 PetscFunctionReturn(0); 3575 } 3576 3577 PetscErrorCode SNESMonitorLGCreate(MPI_Comm comm,const char host[],const char label[],int x,int y,int m,int n,PetscDrawLG *lgctx) 3578 { 3579 PetscErrorCode ierr; 3580 3581 PetscFunctionBegin; 3582 ierr = KSPMonitorLGResidualNormCreate(comm,host,label,x,y,m,n,lgctx);CHKERRQ(ierr); 3583 PetscFunctionReturn(0); 3584 } 3585 3586 PETSC_INTERN PetscErrorCode SNESMonitorRange_Private(SNES,PetscInt,PetscReal*); 3587 3588 PetscErrorCode SNESMonitorLGRange(SNES snes,PetscInt n,PetscReal rnorm,void *monctx) 3589 { 3590 PetscDrawLG lg; 3591 PetscErrorCode ierr; 3592 PetscReal x,y,per; 3593 PetscViewer v = (PetscViewer)monctx; 3594 static PetscReal prev; /* should be in the context */ 3595 PetscDraw draw; 3596 3597 PetscFunctionBegin; 3598 PetscValidHeaderSpecific(v,PETSC_VIEWER_CLASSID,4); 3599 ierr = PetscViewerDrawGetDrawLG(v,0,&lg);CHKERRQ(ierr); 3600 if (!n) {ierr = PetscDrawLGReset(lg);CHKERRQ(ierr);} 3601 ierr = PetscDrawLGGetDraw(lg,&draw);CHKERRQ(ierr); 3602 ierr = PetscDrawSetTitle(draw,"Residual norm");CHKERRQ(ierr); 3603 x = (PetscReal)n; 3604 if (rnorm > 0.0) y = PetscLog10Real(rnorm); 3605 else y = -15.0; 3606 ierr = PetscDrawLGAddPoint(lg,&x,&y);CHKERRQ(ierr); 3607 if (n < 20 || !(n % 5) || snes->reason) { 3608 ierr = PetscDrawLGDraw(lg);CHKERRQ(ierr); 3609 ierr = PetscDrawLGSave(lg);CHKERRQ(ierr); 3610 } 3611 3612 ierr = PetscViewerDrawGetDrawLG(v,1,&lg);CHKERRQ(ierr); 3613 if (!n) {ierr = PetscDrawLGReset(lg);CHKERRQ(ierr);} 3614 ierr = PetscDrawLGGetDraw(lg,&draw);CHKERRQ(ierr); 3615 ierr = PetscDrawSetTitle(draw,"% elemts > .2*max elemt");CHKERRQ(ierr); 3616 ierr = SNESMonitorRange_Private(snes,n,&per);CHKERRQ(ierr); 3617 x = (PetscReal)n; 3618 y = 100.0*per; 3619 ierr = PetscDrawLGAddPoint(lg,&x,&y);CHKERRQ(ierr); 3620 if (n < 20 || !(n % 5) || snes->reason) { 3621 ierr = PetscDrawLGDraw(lg);CHKERRQ(ierr); 3622 ierr = PetscDrawLGSave(lg);CHKERRQ(ierr); 3623 } 3624 3625 ierr = PetscViewerDrawGetDrawLG(v,2,&lg);CHKERRQ(ierr); 3626 if (!n) {prev = rnorm;ierr = PetscDrawLGReset(lg);CHKERRQ(ierr);} 3627 ierr = PetscDrawLGGetDraw(lg,&draw);CHKERRQ(ierr); 3628 ierr = PetscDrawSetTitle(draw,"(norm -oldnorm)/oldnorm");CHKERRQ(ierr); 3629 x = (PetscReal)n; 3630 y = (prev - rnorm)/prev; 3631 ierr = PetscDrawLGAddPoint(lg,&x,&y);CHKERRQ(ierr); 3632 if (n < 20 || !(n % 5) || snes->reason) { 3633 ierr = PetscDrawLGDraw(lg);CHKERRQ(ierr); 3634 ierr = PetscDrawLGSave(lg);CHKERRQ(ierr); 3635 } 3636 3637 ierr = PetscViewerDrawGetDrawLG(v,3,&lg);CHKERRQ(ierr); 3638 if (!n) {ierr = PetscDrawLGReset(lg);CHKERRQ(ierr);} 3639 ierr = PetscDrawLGGetDraw(lg,&draw);CHKERRQ(ierr); 3640 ierr = PetscDrawSetTitle(draw,"(norm -oldnorm)/oldnorm*(% > .2 max)");CHKERRQ(ierr); 3641 x = (PetscReal)n; 3642 y = (prev - rnorm)/(prev*per); 3643 if (n > 2) { /*skip initial crazy value */ 3644 ierr = PetscDrawLGAddPoint(lg,&x,&y);CHKERRQ(ierr); 3645 } 3646 if (n < 20 || !(n % 5) || snes->reason) { 3647 ierr = PetscDrawLGDraw(lg);CHKERRQ(ierr); 3648 ierr = PetscDrawLGSave(lg);CHKERRQ(ierr); 3649 } 3650 prev = rnorm; 3651 PetscFunctionReturn(0); 3652 } 3653 3654 /*@ 3655 SNESMonitor - runs the user provided monitor routines, if they exist 3656 3657 Collective on SNES 3658 3659 Input Parameters: 3660 + snes - nonlinear solver context obtained from SNESCreate() 3661 . iter - iteration number 3662 - rnorm - relative norm of the residual 3663 3664 Notes: 3665 This routine is called by the SNES implementations. 3666 It does not typically need to be called by the user. 3667 3668 Level: developer 3669 3670 .seealso: SNESMonitorSet() 3671 @*/ 3672 PetscErrorCode SNESMonitor(SNES snes,PetscInt iter,PetscReal rnorm) 3673 { 3674 PetscErrorCode ierr; 3675 PetscInt i,n = snes->numbermonitors; 3676 3677 PetscFunctionBegin; 3678 ierr = VecLockPush(snes->vec_sol);CHKERRQ(ierr); 3679 for (i=0; i<n; i++) { 3680 ierr = (*snes->monitor[i])(snes,iter,rnorm,snes->monitorcontext[i]);CHKERRQ(ierr); 3681 } 3682 ierr = VecLockPop(snes->vec_sol);CHKERRQ(ierr); 3683 PetscFunctionReturn(0); 3684 } 3685 3686 /* ------------ Routines to set performance monitoring options ----------- */ 3687 3688 /*MC 3689 SNESMonitorFunction - functional form passed to SNESMonitorSet() to monitor convergence of nonlinear solver 3690 3691 Synopsis: 3692 #include <petscsnes.h> 3693 $ PetscErrorCode SNESMonitorFunction(SNES snes,PetscInt its, PetscReal norm,void *mctx) 3694 3695 + snes - the SNES context 3696 . its - iteration number 3697 . norm - 2-norm function value (may be estimated) 3698 - mctx - [optional] monitoring context 3699 3700 Level: advanced 3701 3702 .seealso: SNESMonitorSet(), SNESMonitorGet() 3703 M*/ 3704 3705 /*@C 3706 SNESMonitorSet - Sets an ADDITIONAL function that is to be used at every 3707 iteration of the nonlinear solver to display the iteration's 3708 progress. 3709 3710 Logically Collective on SNES 3711 3712 Input Parameters: 3713 + snes - the SNES context 3714 . f - the monitor function, see SNESMonitorFunction for the calling sequence 3715 . mctx - [optional] user-defined context for private data for the 3716 monitor routine (use NULL if no context is desired) 3717 - monitordestroy - [optional] routine that frees monitor context 3718 (may be NULL) 3719 3720 Options Database Keys: 3721 + -snes_monitor - sets SNESMonitorDefault() 3722 . -snes_monitor_lg_residualnorm - sets line graph monitor, 3723 uses SNESMonitorLGCreate() 3724 - -snes_monitor_cancel - cancels all monitors that have 3725 been hardwired into a code by 3726 calls to SNESMonitorSet(), but 3727 does not cancel those set via 3728 the options database. 3729 3730 Notes: 3731 Several different monitoring routines may be set by calling 3732 SNESMonitorSet() multiple times; all will be called in the 3733 order in which they were set. 3734 3735 Fortran notes: Only a single monitor function can be set for each SNES object 3736 3737 Level: intermediate 3738 3739 .keywords: SNES, nonlinear, set, monitor 3740 3741 .seealso: SNESMonitorDefault(), SNESMonitorCancel(), SNESMonitorFunction 3742 @*/ 3743 PetscErrorCode SNESMonitorSet(SNES snes,PetscErrorCode (*f)(SNES,PetscInt,PetscReal,void*),void *mctx,PetscErrorCode (*monitordestroy)(void**)) 3744 { 3745 PetscInt i; 3746 PetscErrorCode ierr; 3747 PetscBool identical; 3748 3749 PetscFunctionBegin; 3750 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3751 for (i=0; i<snes->numbermonitors;i++) { 3752 ierr = PetscMonitorCompare((PetscErrorCode (*)(void))f,mctx,monitordestroy,(PetscErrorCode (*)(void))snes->monitor[i],snes->monitorcontext[i],snes->monitordestroy[i],&identical);CHKERRQ(ierr); 3753 if (identical) PetscFunctionReturn(0); 3754 } 3755 if (snes->numbermonitors >= MAXSNESMONITORS) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Too many monitors set"); 3756 snes->monitor[snes->numbermonitors] = f; 3757 snes->monitordestroy[snes->numbermonitors] = monitordestroy; 3758 snes->monitorcontext[snes->numbermonitors++] = (void*)mctx; 3759 PetscFunctionReturn(0); 3760 } 3761 3762 /*@ 3763 SNESMonitorCancel - Clears all the monitor functions for a SNES object. 3764 3765 Logically Collective on SNES 3766 3767 Input Parameters: 3768 . snes - the SNES context 3769 3770 Options Database Key: 3771 . -snes_monitor_cancel - cancels all monitors that have been hardwired 3772 into a code by calls to SNESMonitorSet(), but does not cancel those 3773 set via the options database 3774 3775 Notes: 3776 There is no way to clear one specific monitor from a SNES object. 3777 3778 Level: intermediate 3779 3780 .keywords: SNES, nonlinear, set, monitor 3781 3782 .seealso: SNESMonitorDefault(), SNESMonitorSet() 3783 @*/ 3784 PetscErrorCode SNESMonitorCancel(SNES snes) 3785 { 3786 PetscErrorCode ierr; 3787 PetscInt i; 3788 3789 PetscFunctionBegin; 3790 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3791 for (i=0; i<snes->numbermonitors; i++) { 3792 if (snes->monitordestroy[i]) { 3793 ierr = (*snes->monitordestroy[i])(&snes->monitorcontext[i]);CHKERRQ(ierr); 3794 } 3795 } 3796 snes->numbermonitors = 0; 3797 PetscFunctionReturn(0); 3798 } 3799 3800 /*MC 3801 SNESConvergenceTestFunction - functional form used for testing of convergence of nonlinear solver 3802 3803 Synopsis: 3804 #include <petscsnes.h> 3805 $ PetscErrorCode SNESConvergenceTest(SNES snes,PetscInt it,PetscReal xnorm,PetscReal gnorm,PetscReal f,SNESConvergedReason *reason,void *cctx) 3806 3807 + snes - the SNES context 3808 . it - current iteration (0 is the first and is before any Newton step) 3809 . cctx - [optional] convergence context 3810 . reason - reason for convergence/divergence 3811 . xnorm - 2-norm of current iterate 3812 . gnorm - 2-norm of current step 3813 - f - 2-norm of function 3814 3815 Level: intermediate 3816 3817 .seealso: SNESSetConvergenceTest(), SNESGetConvergenceTest() 3818 M*/ 3819 3820 /*@C 3821 SNESSetConvergenceTest - Sets the function that is to be used 3822 to test for convergence of the nonlinear iterative solution. 3823 3824 Logically Collective on SNES 3825 3826 Input Parameters: 3827 + snes - the SNES context 3828 . SNESConvergenceTestFunction - routine to test for convergence 3829 . cctx - [optional] context for private data for the convergence routine (may be NULL) 3830 - destroy - [optional] destructor for the context (may be NULL; NULL_FUNCTION in Fortran) 3831 3832 Level: advanced 3833 3834 .keywords: SNES, nonlinear, set, convergence, test 3835 3836 .seealso: SNESConvergedDefault(), SNESConvergedSkip(), SNESConvergenceTestFunction 3837 @*/ 3838 PetscErrorCode SNESSetConvergenceTest(SNES snes,PetscErrorCode (*SNESConvergenceTestFunction)(SNES,PetscInt,PetscReal,PetscReal,PetscReal,SNESConvergedReason*,void*),void *cctx,PetscErrorCode (*destroy)(void*)) 3839 { 3840 PetscErrorCode ierr; 3841 3842 PetscFunctionBegin; 3843 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3844 if (!SNESConvergenceTestFunction) SNESConvergenceTestFunction = SNESConvergedSkip; 3845 if (snes->ops->convergeddestroy) { 3846 ierr = (*snes->ops->convergeddestroy)(snes->cnvP);CHKERRQ(ierr); 3847 } 3848 snes->ops->converged = SNESConvergenceTestFunction; 3849 snes->ops->convergeddestroy = destroy; 3850 snes->cnvP = cctx; 3851 PetscFunctionReturn(0); 3852 } 3853 3854 /*@ 3855 SNESGetConvergedReason - Gets the reason the SNES iteration was stopped. 3856 3857 Not Collective 3858 3859 Input Parameter: 3860 . snes - the SNES context 3861 3862 Output Parameter: 3863 . reason - negative value indicates diverged, positive value converged, see SNESConvergedReason or the 3864 manual pages for the individual convergence tests for complete lists 3865 3866 Options Database: 3867 . -snes_converged_reason - prints the reason to standard out 3868 3869 Level: intermediate 3870 3871 Notes: Should only be called after the call the SNESSolve() is complete, if it is called earlier it returns the value SNES__CONVERGED_ITERATING. 3872 3873 .keywords: SNES, nonlinear, set, convergence, test 3874 3875 .seealso: SNESSetConvergenceTest(), SNESSetConvergedReason(), SNESConvergedReason 3876 @*/ 3877 PetscErrorCode SNESGetConvergedReason(SNES snes,SNESConvergedReason *reason) 3878 { 3879 PetscFunctionBegin; 3880 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3881 PetscValidPointer(reason,2); 3882 *reason = snes->reason; 3883 PetscFunctionReturn(0); 3884 } 3885 3886 /*@ 3887 SNESSetConvergedReason - Sets the reason the SNES iteration was stopped. 3888 3889 Not Collective 3890 3891 Input Parameters: 3892 + snes - the SNES context 3893 - reason - negative value indicates diverged, positive value converged, see SNESConvergedReason or the 3894 manual pages for the individual convergence tests for complete lists 3895 3896 Level: intermediate 3897 3898 .keywords: SNES, nonlinear, set, convergence, test 3899 .seealso: SNESGetConvergedReason(), SNESSetConvergenceTest(), SNESConvergedReason 3900 @*/ 3901 PetscErrorCode SNESSetConvergedReason(SNES snes,SNESConvergedReason reason) 3902 { 3903 PetscFunctionBegin; 3904 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3905 snes->reason = reason; 3906 PetscFunctionReturn(0); 3907 } 3908 3909 /*@ 3910 SNESSetConvergenceHistory - Sets the array used to hold the convergence history. 3911 3912 Logically Collective on SNES 3913 3914 Input Parameters: 3915 + snes - iterative context obtained from SNESCreate() 3916 . a - array to hold history, this array will contain the function norms computed at each step 3917 . its - integer array holds the number of linear iterations for each solve. 3918 . na - size of a and its 3919 - reset - PETSC_TRUE indicates each new nonlinear solve resets the history counter to zero, 3920 else it continues storing new values for new nonlinear solves after the old ones 3921 3922 Notes: 3923 If 'a' and 'its' are NULL then space is allocated for the history. If 'na' PETSC_DECIDE or PETSC_DEFAULT then a 3924 default array of length 10000 is allocated. 3925 3926 This routine is useful, e.g., when running a code for purposes 3927 of accurate performance monitoring, when no I/O should be done 3928 during the section of code that is being timed. 3929 3930 Level: intermediate 3931 3932 .keywords: SNES, set, convergence, history 3933 3934 .seealso: SNESGetConvergenceHistory() 3935 3936 @*/ 3937 PetscErrorCode SNESSetConvergenceHistory(SNES snes,PetscReal a[],PetscInt its[],PetscInt na,PetscBool reset) 3938 { 3939 PetscErrorCode ierr; 3940 3941 PetscFunctionBegin; 3942 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3943 if (a) PetscValidScalarPointer(a,2); 3944 if (its) PetscValidIntPointer(its,3); 3945 if (!a) { 3946 if (na == PETSC_DECIDE || na == PETSC_DEFAULT) na = 1000; 3947 ierr = PetscCalloc1(na,&a);CHKERRQ(ierr); 3948 ierr = PetscCalloc1(na,&its);CHKERRQ(ierr); 3949 3950 snes->conv_malloc = PETSC_TRUE; 3951 } 3952 snes->conv_hist = a; 3953 snes->conv_hist_its = its; 3954 snes->conv_hist_max = na; 3955 snes->conv_hist_len = 0; 3956 snes->conv_hist_reset = reset; 3957 PetscFunctionReturn(0); 3958 } 3959 3960 #if defined(PETSC_HAVE_MATLAB_ENGINE) 3961 #include <engine.h> /* MATLAB include file */ 3962 #include <mex.h> /* MATLAB include file */ 3963 3964 PETSC_EXTERN mxArray *SNESGetConvergenceHistoryMatlab(SNES snes) 3965 { 3966 mxArray *mat; 3967 PetscInt i; 3968 PetscReal *ar; 3969 3970 PetscFunctionBegin; 3971 mat = mxCreateDoubleMatrix(snes->conv_hist_len,1,mxREAL); 3972 ar = (PetscReal*) mxGetData(mat); 3973 for (i=0; i<snes->conv_hist_len; i++) ar[i] = snes->conv_hist[i]; 3974 PetscFunctionReturn(mat); 3975 } 3976 #endif 3977 3978 /*@C 3979 SNESGetConvergenceHistory - Gets the array used to hold the convergence history. 3980 3981 Not Collective 3982 3983 Input Parameter: 3984 . snes - iterative context obtained from SNESCreate() 3985 3986 Output Parameters: 3987 . a - array to hold history 3988 . its - integer array holds the number of linear iterations (or 3989 negative if not converged) for each solve. 3990 - na - size of a and its 3991 3992 Notes: 3993 The calling sequence for this routine in Fortran is 3994 $ call SNESGetConvergenceHistory(SNES snes, integer na, integer ierr) 3995 3996 This routine is useful, e.g., when running a code for purposes 3997 of accurate performance monitoring, when no I/O should be done 3998 during the section of code that is being timed. 3999 4000 Level: intermediate 4001 4002 .keywords: SNES, get, convergence, history 4003 4004 .seealso: SNESSetConvergencHistory() 4005 4006 @*/ 4007 PetscErrorCode SNESGetConvergenceHistory(SNES snes,PetscReal *a[],PetscInt *its[],PetscInt *na) 4008 { 4009 PetscFunctionBegin; 4010 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4011 if (a) *a = snes->conv_hist; 4012 if (its) *its = snes->conv_hist_its; 4013 if (na) *na = snes->conv_hist_len; 4014 PetscFunctionReturn(0); 4015 } 4016 4017 /*@C 4018 SNESSetUpdate - Sets the general-purpose update function called 4019 at the beginning of every iteration of the nonlinear solve. Specifically 4020 it is called just before the Jacobian is "evaluated". 4021 4022 Logically Collective on SNES 4023 4024 Input Parameters: 4025 . snes - The nonlinear solver context 4026 . func - The function 4027 4028 Calling sequence of func: 4029 . func (SNES snes, PetscInt step); 4030 4031 . step - The current step of the iteration 4032 4033 Level: advanced 4034 4035 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() 4036 This is not used by most users. 4037 4038 .keywords: SNES, update 4039 4040 .seealso SNESSetJacobian(), SNESSolve() 4041 @*/ 4042 PetscErrorCode SNESSetUpdate(SNES snes, PetscErrorCode (*func)(SNES, PetscInt)) 4043 { 4044 PetscFunctionBegin; 4045 PetscValidHeaderSpecific(snes, SNES_CLASSID,1); 4046 snes->ops->update = func; 4047 PetscFunctionReturn(0); 4048 } 4049 4050 /* 4051 SNESScaleStep_Private - Scales a step so that its length is less than the 4052 positive parameter delta. 4053 4054 Input Parameters: 4055 + snes - the SNES context 4056 . y - approximate solution of linear system 4057 . fnorm - 2-norm of current function 4058 - delta - trust region size 4059 4060 Output Parameters: 4061 + gpnorm - predicted function norm at the new point, assuming local 4062 linearization. The value is zero if the step lies within the trust 4063 region, and exceeds zero otherwise. 4064 - ynorm - 2-norm of the step 4065 4066 Note: 4067 For non-trust region methods such as SNESNEWTONLS, the parameter delta 4068 is set to be the maximum allowable step size. 4069 4070 .keywords: SNES, nonlinear, scale, step 4071 */ 4072 PetscErrorCode SNESScaleStep_Private(SNES snes,Vec y,PetscReal *fnorm,PetscReal *delta,PetscReal *gpnorm,PetscReal *ynorm) 4073 { 4074 PetscReal nrm; 4075 PetscScalar cnorm; 4076 PetscErrorCode ierr; 4077 4078 PetscFunctionBegin; 4079 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4080 PetscValidHeaderSpecific(y,VEC_CLASSID,2); 4081 PetscCheckSameComm(snes,1,y,2); 4082 4083 ierr = VecNorm(y,NORM_2,&nrm);CHKERRQ(ierr); 4084 if (nrm > *delta) { 4085 nrm = *delta/nrm; 4086 *gpnorm = (1.0 - nrm)*(*fnorm); 4087 cnorm = nrm; 4088 ierr = VecScale(y,cnorm);CHKERRQ(ierr); 4089 *ynorm = *delta; 4090 } else { 4091 *gpnorm = 0.0; 4092 *ynorm = nrm; 4093 } 4094 PetscFunctionReturn(0); 4095 } 4096 4097 /*@ 4098 SNESReasonView - Displays the reason a SNES solve converged or diverged to a viewer 4099 4100 Collective on SNES 4101 4102 Parameter: 4103 + snes - iterative context obtained from SNESCreate() 4104 - viewer - the viewer to display the reason 4105 4106 4107 Options Database Keys: 4108 . -snes_converged_reason - print reason for converged or diverged, also prints number of iterations 4109 4110 Level: beginner 4111 4112 .keywords: SNES, solve, linear system 4113 4114 .seealso: SNESCreate(), SNESSetUp(), SNESDestroy(), SNESSetTolerances(), SNESConvergedDefault() 4115 4116 @*/ 4117 PetscErrorCode SNESReasonView(SNES snes,PetscViewer viewer) 4118 { 4119 PetscViewerFormat format; 4120 PetscBool isAscii; 4121 PetscErrorCode ierr; 4122 4123 PetscFunctionBegin; 4124 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isAscii);CHKERRQ(ierr); 4125 if (isAscii) { 4126 ierr = PetscViewerGetFormat(viewer, &format);CHKERRQ(ierr); 4127 ierr = PetscViewerASCIIAddTab(viewer,((PetscObject)snes)->tablevel);CHKERRQ(ierr); 4128 if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) { 4129 DM dm; 4130 Vec u; 4131 PetscDS prob; 4132 PetscInt Nf, f; 4133 PetscErrorCode (**exactFuncs)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar[], void *); 4134 PetscReal error; 4135 4136 ierr = SNESGetDM(snes, &dm);CHKERRQ(ierr); 4137 ierr = SNESGetSolution(snes, &u);CHKERRQ(ierr); 4138 ierr = DMGetDS(dm, &prob);CHKERRQ(ierr); 4139 ierr = PetscDSGetNumFields(prob, &Nf);CHKERRQ(ierr); 4140 ierr = PetscMalloc1(Nf, &exactFuncs);CHKERRQ(ierr); 4141 for (f = 0; f < Nf; ++f) {ierr = PetscDSGetExactSolution(prob, f, &exactFuncs[f]);CHKERRQ(ierr);} 4142 ierr = DMComputeL2Diff(dm, 0.0, exactFuncs, NULL, u, &error);CHKERRQ(ierr); 4143 ierr = PetscFree(exactFuncs);CHKERRQ(ierr); 4144 if (error < 1.0e-11) {ierr = PetscViewerASCIIPrintf(viewer, "L_2 Error: < 1.0e-11\n");CHKERRQ(ierr);} 4145 else {ierr = PetscViewerASCIIPrintf(viewer, "L_2 Error: %g\n", error);CHKERRQ(ierr);} 4146 } 4147 if (snes->reason > 0) { 4148 if (((PetscObject) snes)->prefix) { 4149 ierr = PetscViewerASCIIPrintf(viewer,"Nonlinear %s solve converged due to %s iterations %D\n",((PetscObject) snes)->prefix,SNESConvergedReasons[snes->reason],snes->iter);CHKERRQ(ierr); 4150 } else { 4151 ierr = PetscViewerASCIIPrintf(viewer,"Nonlinear solve converged due to %s iterations %D\n",SNESConvergedReasons[snes->reason],snes->iter);CHKERRQ(ierr); 4152 } 4153 } else { 4154 if (((PetscObject) snes)->prefix) { 4155 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); 4156 } else { 4157 ierr = PetscViewerASCIIPrintf(viewer,"Nonlinear solve did not converge due to %s iterations %D\n",SNESConvergedReasons[snes->reason],snes->iter);CHKERRQ(ierr); 4158 } 4159 } 4160 ierr = PetscViewerASCIISubtractTab(viewer,((PetscObject)snes)->tablevel);CHKERRQ(ierr); 4161 } 4162 PetscFunctionReturn(0); 4163 } 4164 4165 /*@C 4166 SNESReasonViewFromOptions - Processes command line options to determine if/how a SNESReason is to be viewed. 4167 4168 Collective on SNES 4169 4170 Input Parameters: 4171 . snes - the SNES object 4172 4173 Level: intermediate 4174 4175 @*/ 4176 PetscErrorCode SNESReasonViewFromOptions(SNES snes) 4177 { 4178 PetscErrorCode ierr; 4179 PetscViewer viewer; 4180 PetscBool flg; 4181 static PetscBool incall = PETSC_FALSE; 4182 PetscViewerFormat format; 4183 4184 PetscFunctionBegin; 4185 if (incall) PetscFunctionReturn(0); 4186 incall = PETSC_TRUE; 4187 ierr = PetscOptionsGetViewer(PetscObjectComm((PetscObject)snes),((PetscObject)snes)->prefix,"-snes_converged_reason",&viewer,&format,&flg);CHKERRQ(ierr); 4188 if (flg) { 4189 ierr = PetscViewerPushFormat(viewer,format);CHKERRQ(ierr); 4190 ierr = SNESReasonView(snes,viewer);CHKERRQ(ierr); 4191 ierr = PetscViewerPopFormat(viewer);CHKERRQ(ierr); 4192 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 4193 } 4194 incall = PETSC_FALSE; 4195 PetscFunctionReturn(0); 4196 } 4197 4198 /*@C 4199 SNESSolve - Solves a nonlinear system F(x) = b. 4200 Call SNESSolve() after calling SNESCreate() and optional routines of the form SNESSetXXX(). 4201 4202 Collective on SNES 4203 4204 Input Parameters: 4205 + snes - the SNES context 4206 . b - the constant part of the equation F(x) = b, or NULL to use zero. 4207 - x - the solution vector. 4208 4209 Notes: 4210 The user should initialize the vector,x, with the initial guess 4211 for the nonlinear solve prior to calling SNESSolve. In particular, 4212 to employ an initial guess of zero, the user should explicitly set 4213 this vector to zero by calling VecSet(). 4214 4215 Level: beginner 4216 4217 .keywords: SNES, nonlinear, solve 4218 4219 .seealso: SNESCreate(), SNESDestroy(), SNESSetFunction(), SNESSetJacobian(), SNESSetGridSequence(), SNESGetSolution() 4220 @*/ 4221 PetscErrorCode SNESSolve(SNES snes,Vec b,Vec x) 4222 { 4223 PetscErrorCode ierr; 4224 PetscBool flg; 4225 PetscInt grid; 4226 Vec xcreated = NULL; 4227 DM dm; 4228 4229 PetscFunctionBegin; 4230 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4231 if (x) PetscValidHeaderSpecific(x,VEC_CLASSID,3); 4232 if (x) PetscCheckSameComm(snes,1,x,3); 4233 if (b) PetscValidHeaderSpecific(b,VEC_CLASSID,2); 4234 if (b) PetscCheckSameComm(snes,1,b,2); 4235 4236 /* High level operations using the nonlinear solver */ 4237 { 4238 PetscViewer viewer; 4239 PetscViewerFormat format; 4240 PetscInt num; 4241 PetscBool flg; 4242 static PetscBool incall = PETSC_FALSE; 4243 4244 if (!incall) { 4245 /* Estimate the convergence rate of the discretization */ 4246 ierr = PetscOptionsGetViewer(PetscObjectComm((PetscObject) snes), ((PetscObject) snes)->prefix, "-snes_convergence_estimate", &viewer, &format, &flg);CHKERRQ(ierr); 4247 if (flg) { 4248 PetscConvEst conv; 4249 PetscReal alpha; /* Convergence rate of the solution error in the L_2 norm */ 4250 4251 incall = PETSC_TRUE; 4252 ierr = PetscConvEstCreate(PetscObjectComm((PetscObject) snes), &conv);CHKERRQ(ierr); 4253 ierr = PetscConvEstSetSolver(conv, snes);CHKERRQ(ierr); 4254 ierr = PetscConvEstSetFromOptions(conv);CHKERRQ(ierr); 4255 ierr = PetscConvEstSetUp(conv);CHKERRQ(ierr); 4256 ierr = PetscConvEstGetConvRate(conv, &alpha);CHKERRQ(ierr); 4257 ierr = PetscViewerPushFormat(viewer, format);CHKERRQ(ierr); 4258 ierr = PetscConvEstRateView(conv, alpha, viewer);CHKERRQ(ierr); 4259 ierr = PetscViewerPopFormat(viewer);CHKERRQ(ierr); 4260 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 4261 ierr = PetscConvEstDestroy(&conv);CHKERRQ(ierr); 4262 incall = PETSC_FALSE; 4263 } 4264 /* Adaptively refine the initial grid */ 4265 num = 1; 4266 ierr = PetscOptionsGetInt(NULL, ((PetscObject) snes)->prefix, "-snes_adapt_initial", &num, &flg);CHKERRQ(ierr); 4267 if (flg) { 4268 DMAdaptor adaptor; 4269 4270 incall = PETSC_TRUE; 4271 ierr = DMAdaptorCreate(PETSC_COMM_WORLD, &adaptor);CHKERRQ(ierr); 4272 ierr = DMAdaptorSetSolver(adaptor, snes);CHKERRQ(ierr); 4273 ierr = DMAdaptorSetSequenceLength(adaptor, num);CHKERRQ(ierr); 4274 ierr = DMAdaptorSetFromOptions(adaptor);CHKERRQ(ierr); 4275 ierr = DMAdaptorSetUp(adaptor);CHKERRQ(ierr); 4276 ierr = DMAdaptorAdapt(adaptor, x, DM_ADAPTATION_INITIAL, &dm, &x);CHKERRQ(ierr); 4277 ierr = DMAdaptorDestroy(&adaptor);CHKERRQ(ierr); 4278 incall = PETSC_FALSE; 4279 } 4280 /* Use grid sequencing to adapt */ 4281 num = 0; 4282 ierr = PetscOptionsGetInt(NULL, ((PetscObject) snes)->prefix, "-snes_adapt_sequence", &num, NULL);CHKERRQ(ierr); 4283 if (num) { 4284 DMAdaptor adaptor; 4285 4286 incall = PETSC_TRUE; 4287 ierr = DMAdaptorCreate(PETSC_COMM_WORLD, &adaptor);CHKERRQ(ierr); 4288 ierr = DMAdaptorSetSolver(adaptor, snes);CHKERRQ(ierr); 4289 ierr = DMAdaptorSetSequenceLength(adaptor, num);CHKERRQ(ierr); 4290 ierr = DMAdaptorSetFromOptions(adaptor);CHKERRQ(ierr); 4291 ierr = DMAdaptorSetUp(adaptor);CHKERRQ(ierr); 4292 ierr = DMAdaptorAdapt(adaptor, x, DM_ADAPTATION_SEQUENTIAL, &dm, &x);CHKERRQ(ierr); 4293 ierr = DMAdaptorDestroy(&adaptor);CHKERRQ(ierr); 4294 incall = PETSC_FALSE; 4295 } 4296 } 4297 } 4298 if (!x) { 4299 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 4300 ierr = DMCreateGlobalVector(dm,&xcreated);CHKERRQ(ierr); 4301 x = xcreated; 4302 } 4303 ierr = SNESViewFromOptions(snes,NULL,"-snes_view_pre");CHKERRQ(ierr); 4304 4305 for (grid=0; grid<snes->gridsequence; grid++) {ierr = PetscViewerASCIIPushTab(PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)snes)));CHKERRQ(ierr);} 4306 for (grid=0; grid<snes->gridsequence+1; grid++) { 4307 4308 /* set solution vector */ 4309 if (!grid) {ierr = PetscObjectReference((PetscObject)x);CHKERRQ(ierr);} 4310 ierr = VecDestroy(&snes->vec_sol);CHKERRQ(ierr); 4311 snes->vec_sol = x; 4312 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 4313 4314 /* set affine vector if provided */ 4315 if (b) { ierr = PetscObjectReference((PetscObject)b);CHKERRQ(ierr); } 4316 ierr = VecDestroy(&snes->vec_rhs);CHKERRQ(ierr); 4317 snes->vec_rhs = b; 4318 4319 if (snes->vec_func == snes->vec_sol) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_IDN,"Solution vector cannot be function vector"); 4320 if (snes->vec_rhs == snes->vec_sol) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_IDN,"Solution vector cannot be right hand side vector"); 4321 if (!snes->vec_sol_update /* && snes->vec_sol */) { 4322 ierr = VecDuplicate(snes->vec_sol,&snes->vec_sol_update);CHKERRQ(ierr); 4323 ierr = PetscLogObjectParent((PetscObject)snes,(PetscObject)snes->vec_sol_update);CHKERRQ(ierr); 4324 } 4325 ierr = DMShellSetGlobalVector(dm,snes->vec_sol);CHKERRQ(ierr); 4326 ierr = SNESSetUp(snes);CHKERRQ(ierr); 4327 4328 if (!grid) { 4329 if (snes->ops->computeinitialguess) { 4330 ierr = (*snes->ops->computeinitialguess)(snes,snes->vec_sol,snes->initialguessP);CHKERRQ(ierr); 4331 } 4332 } 4333 4334 if (snes->conv_hist_reset) snes->conv_hist_len = 0; 4335 if (snes->counters_reset) {snes->nfuncs = 0; snes->linear_its = 0; snes->numFailures = 0;} 4336 4337 ierr = PetscLogEventBegin(SNES_Solve,snes,0,0,0);CHKERRQ(ierr); 4338 ierr = (*snes->ops->solve)(snes);CHKERRQ(ierr); 4339 ierr = PetscLogEventEnd(SNES_Solve,snes,0,0,0);CHKERRQ(ierr); 4340 if (!snes->reason) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Internal error, solver returned without setting converged reason"); 4341 snes->domainerror = PETSC_FALSE; /* clear the flag if it has been set */ 4342 4343 if (snes->lagjac_persist) snes->jac_iter += snes->iter; 4344 if (snes->lagpre_persist) snes->pre_iter += snes->iter; 4345 4346 ierr = PetscOptionsGetViewer(PetscObjectComm((PetscObject)snes),((PetscObject)snes)->prefix,"-snes_test_local_min",NULL,NULL,&flg);CHKERRQ(ierr); 4347 if (flg && !PetscPreLoadingOn) { ierr = SNESTestLocalMin(snes);CHKERRQ(ierr); } 4348 ierr = SNESReasonViewFromOptions(snes);CHKERRQ(ierr); 4349 4350 if (snes->errorifnotconverged && snes->reason < 0) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_NOT_CONVERGED,"SNESSolve has not converged"); 4351 if (snes->reason < 0) break; 4352 if (grid < snes->gridsequence) { 4353 DM fine; 4354 Vec xnew; 4355 Mat interp; 4356 4357 ierr = DMRefine(snes->dm,PetscObjectComm((PetscObject)snes),&fine);CHKERRQ(ierr); 4358 if (!fine) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_INCOMP,"DMRefine() did not perform any refinement, cannot continue grid sequencing"); 4359 ierr = DMCreateInterpolation(snes->dm,fine,&interp,NULL);CHKERRQ(ierr); 4360 ierr = DMCreateGlobalVector(fine,&xnew);CHKERRQ(ierr); 4361 ierr = MatInterpolate(interp,x,xnew);CHKERRQ(ierr); 4362 ierr = DMInterpolate(snes->dm,interp,fine);CHKERRQ(ierr); 4363 ierr = MatDestroy(&interp);CHKERRQ(ierr); 4364 x = xnew; 4365 4366 ierr = SNESReset(snes);CHKERRQ(ierr); 4367 ierr = SNESSetDM(snes,fine);CHKERRQ(ierr); 4368 ierr = DMDestroy(&fine);CHKERRQ(ierr); 4369 ierr = PetscViewerASCIIPopTab(PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)snes)));CHKERRQ(ierr); 4370 } 4371 } 4372 ierr = SNESViewFromOptions(snes,NULL,"-snes_view");CHKERRQ(ierr); 4373 ierr = VecViewFromOptions(snes->vec_sol,(PetscObject)snes,"-snes_view_solution");CHKERRQ(ierr); 4374 4375 ierr = VecDestroy(&xcreated);CHKERRQ(ierr); 4376 ierr = PetscObjectSAWsBlock((PetscObject)snes);CHKERRQ(ierr); 4377 PetscFunctionReturn(0); 4378 } 4379 4380 /* --------- Internal routines for SNES Package --------- */ 4381 4382 /*@C 4383 SNESSetType - Sets the method for the nonlinear solver. 4384 4385 Collective on SNES 4386 4387 Input Parameters: 4388 + snes - the SNES context 4389 - type - a known method 4390 4391 Options Database Key: 4392 . -snes_type <type> - Sets the method; use -help for a list 4393 of available methods (for instance, newtonls or newtontr) 4394 4395 Notes: 4396 See "petsc/include/petscsnes.h" for available methods (for instance) 4397 + SNESNEWTONLS - Newton's method with line search 4398 (systems of nonlinear equations) 4399 . SNESNEWTONTR - Newton's method with trust region 4400 (systems of nonlinear equations) 4401 4402 Normally, it is best to use the SNESSetFromOptions() command and then 4403 set the SNES solver type from the options database rather than by using 4404 this routine. Using the options database provides the user with 4405 maximum flexibility in evaluating the many nonlinear solvers. 4406 The SNESSetType() routine is provided for those situations where it 4407 is necessary to set the nonlinear solver independently of the command 4408 line or options database. This might be the case, for example, when 4409 the choice of solver changes during the execution of the program, 4410 and the user's application is taking responsibility for choosing the 4411 appropriate method. 4412 4413 Developer Notes: SNESRegister() adds a constructor for a new SNESType to SNESList, SNESSetType() locates 4414 the constructor in that list and calls it to create the spexific object. 4415 4416 Level: intermediate 4417 4418 .keywords: SNES, set, type 4419 4420 .seealso: SNESType, SNESCreate(), SNESDestroy(), SNESGetType(), SNESSetFromOptions() 4421 4422 @*/ 4423 PetscErrorCode SNESSetType(SNES snes,SNESType type) 4424 { 4425 PetscErrorCode ierr,(*r)(SNES); 4426 PetscBool match; 4427 4428 PetscFunctionBegin; 4429 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4430 PetscValidCharPointer(type,2); 4431 4432 ierr = PetscObjectTypeCompare((PetscObject)snes,type,&match);CHKERRQ(ierr); 4433 if (match) PetscFunctionReturn(0); 4434 4435 ierr = PetscFunctionListFind(SNESList,type,&r);CHKERRQ(ierr); 4436 if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unable to find requested SNES type %s",type); 4437 /* Destroy the previous private SNES context */ 4438 if (snes->ops->destroy) { 4439 ierr = (*(snes)->ops->destroy)(snes);CHKERRQ(ierr); 4440 snes->ops->destroy = NULL; 4441 } 4442 /* Reinitialize function pointers in SNESOps structure */ 4443 snes->ops->setup = 0; 4444 snes->ops->solve = 0; 4445 snes->ops->view = 0; 4446 snes->ops->setfromoptions = 0; 4447 snes->ops->destroy = 0; 4448 /* Call the SNESCreate_XXX routine for this particular Nonlinear solver */ 4449 snes->setupcalled = PETSC_FALSE; 4450 4451 ierr = PetscObjectChangeTypeName((PetscObject)snes,type);CHKERRQ(ierr); 4452 ierr = (*r)(snes);CHKERRQ(ierr); 4453 PetscFunctionReturn(0); 4454 } 4455 4456 /*@C 4457 SNESGetType - Gets the SNES method type and name (as a string). 4458 4459 Not Collective 4460 4461 Input Parameter: 4462 . snes - nonlinear solver context 4463 4464 Output Parameter: 4465 . type - SNES method (a character string) 4466 4467 Level: intermediate 4468 4469 .keywords: SNES, nonlinear, get, type, name 4470 @*/ 4471 PetscErrorCode SNESGetType(SNES snes,SNESType *type) 4472 { 4473 PetscFunctionBegin; 4474 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4475 PetscValidPointer(type,2); 4476 *type = ((PetscObject)snes)->type_name; 4477 PetscFunctionReturn(0); 4478 } 4479 4480 /*@ 4481 SNESSetSolution - Sets the solution vector for use by the SNES routines. 4482 4483 Logically Collective on SNES and Vec 4484 4485 Input Parameters: 4486 + snes - the SNES context obtained from SNESCreate() 4487 - u - the solution vector 4488 4489 Level: beginner 4490 4491 .keywords: SNES, set, solution 4492 @*/ 4493 PetscErrorCode SNESSetSolution(SNES snes, Vec u) 4494 { 4495 DM dm; 4496 PetscErrorCode ierr; 4497 4498 PetscFunctionBegin; 4499 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); 4500 PetscValidHeaderSpecific(u, VEC_CLASSID, 2); 4501 ierr = PetscObjectReference((PetscObject) u);CHKERRQ(ierr); 4502 ierr = VecDestroy(&snes->vec_sol);CHKERRQ(ierr); 4503 4504 snes->vec_sol = u; 4505 4506 ierr = SNESGetDM(snes, &dm);CHKERRQ(ierr); 4507 ierr = DMShellSetGlobalVector(dm, u);CHKERRQ(ierr); 4508 PetscFunctionReturn(0); 4509 } 4510 4511 /*@ 4512 SNESGetSolution - Returns the vector where the approximate solution is 4513 stored. This is the fine grid solution when using SNESSetGridSequence(). 4514 4515 Not Collective, but Vec is parallel if SNES is parallel 4516 4517 Input Parameter: 4518 . snes - the SNES context 4519 4520 Output Parameter: 4521 . x - the solution 4522 4523 Level: intermediate 4524 4525 .keywords: SNES, nonlinear, get, solution 4526 4527 .seealso: SNESGetSolutionUpdate(), SNESGetFunction() 4528 @*/ 4529 PetscErrorCode SNESGetSolution(SNES snes,Vec *x) 4530 { 4531 PetscFunctionBegin; 4532 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4533 PetscValidPointer(x,2); 4534 *x = snes->vec_sol; 4535 PetscFunctionReturn(0); 4536 } 4537 4538 /*@ 4539 SNESGetSolutionUpdate - Returns the vector where the solution update is 4540 stored. 4541 4542 Not Collective, but Vec is parallel if SNES is parallel 4543 4544 Input Parameter: 4545 . snes - the SNES context 4546 4547 Output Parameter: 4548 . x - the solution update 4549 4550 Level: advanced 4551 4552 .keywords: SNES, nonlinear, get, solution, update 4553 4554 .seealso: SNESGetSolution(), SNESGetFunction() 4555 @*/ 4556 PetscErrorCode SNESGetSolutionUpdate(SNES snes,Vec *x) 4557 { 4558 PetscFunctionBegin; 4559 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4560 PetscValidPointer(x,2); 4561 *x = snes->vec_sol_update; 4562 PetscFunctionReturn(0); 4563 } 4564 4565 /*@C 4566 SNESGetFunction - Returns the vector where the function is stored. 4567 4568 Not Collective, but Vec is parallel if SNES is parallel. Collective if Vec is requested, but has not been created yet. 4569 4570 Input Parameter: 4571 . snes - the SNES context 4572 4573 Output Parameter: 4574 + r - the vector that is used to store residuals (or NULL if you don't want it) 4575 . f - the function (or NULL if you don't want it); see SNESFunction for calling sequence details 4576 - ctx - the function context (or NULL if you don't want it) 4577 4578 Level: advanced 4579 4580 Notes: The vector r DOES NOT, in general contain the current value of the SNES nonlinear function 4581 4582 .keywords: SNES, nonlinear, get, function 4583 4584 .seealso: SNESSetFunction(), SNESGetSolution(), SNESFunction 4585 @*/ 4586 PetscErrorCode SNESGetFunction(SNES snes,Vec *r,PetscErrorCode (**f)(SNES,Vec,Vec,void*),void **ctx) 4587 { 4588 PetscErrorCode ierr; 4589 DM dm; 4590 4591 PetscFunctionBegin; 4592 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4593 if (r) { 4594 if (!snes->vec_func) { 4595 if (snes->vec_rhs) { 4596 ierr = VecDuplicate(snes->vec_rhs,&snes->vec_func);CHKERRQ(ierr); 4597 } else if (snes->vec_sol) { 4598 ierr = VecDuplicate(snes->vec_sol,&snes->vec_func);CHKERRQ(ierr); 4599 } else if (snes->dm) { 4600 ierr = DMCreateGlobalVector(snes->dm,&snes->vec_func);CHKERRQ(ierr); 4601 } 4602 } 4603 *r = snes->vec_func; 4604 } 4605 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 4606 ierr = DMSNESGetFunction(dm,f,ctx);CHKERRQ(ierr); 4607 PetscFunctionReturn(0); 4608 } 4609 4610 /*@C 4611 SNESGetNGS - Returns the NGS function and context. 4612 4613 Input Parameter: 4614 . snes - the SNES context 4615 4616 Output Parameter: 4617 + f - the function (or NULL) see SNESNGSFunction for details 4618 - ctx - the function context (or NULL) 4619 4620 Level: advanced 4621 4622 .keywords: SNES, nonlinear, get, function 4623 4624 .seealso: SNESSetNGS(), SNESGetFunction() 4625 @*/ 4626 4627 PetscErrorCode SNESGetNGS (SNES snes, PetscErrorCode (**f)(SNES, Vec, Vec, void*), void ** ctx) 4628 { 4629 PetscErrorCode ierr; 4630 DM dm; 4631 4632 PetscFunctionBegin; 4633 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4634 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 4635 ierr = DMSNESGetNGS(dm,f,ctx);CHKERRQ(ierr); 4636 PetscFunctionReturn(0); 4637 } 4638 4639 /*@C 4640 SNESSetOptionsPrefix - Sets the prefix used for searching for all 4641 SNES options in the database. 4642 4643 Logically Collective on SNES 4644 4645 Input Parameter: 4646 + snes - the SNES context 4647 - prefix - the prefix to prepend to all option names 4648 4649 Notes: 4650 A hyphen (-) must NOT be given at the beginning of the prefix name. 4651 The first character of all runtime options is AUTOMATICALLY the hyphen. 4652 4653 Level: advanced 4654 4655 .keywords: SNES, set, options, prefix, database 4656 4657 .seealso: SNESSetFromOptions() 4658 @*/ 4659 PetscErrorCode SNESSetOptionsPrefix(SNES snes,const char prefix[]) 4660 { 4661 PetscErrorCode ierr; 4662 4663 PetscFunctionBegin; 4664 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4665 ierr = PetscObjectSetOptionsPrefix((PetscObject)snes,prefix);CHKERRQ(ierr); 4666 if (!snes->ksp) {ierr = SNESGetKSP(snes,&snes->ksp);CHKERRQ(ierr);} 4667 if (snes->linesearch) { 4668 ierr = SNESGetLineSearch(snes,&snes->linesearch);CHKERRQ(ierr); 4669 ierr = PetscObjectSetOptionsPrefix((PetscObject)snes->linesearch,prefix);CHKERRQ(ierr); 4670 } 4671 ierr = KSPSetOptionsPrefix(snes->ksp,prefix);CHKERRQ(ierr); 4672 PetscFunctionReturn(0); 4673 } 4674 4675 /*@C 4676 SNESAppendOptionsPrefix - Appends to the prefix used for searching for all 4677 SNES options in the database. 4678 4679 Logically Collective on SNES 4680 4681 Input Parameters: 4682 + snes - the SNES context 4683 - prefix - the prefix to prepend to all option names 4684 4685 Notes: 4686 A hyphen (-) must NOT be given at the beginning of the prefix name. 4687 The first character of all runtime options is AUTOMATICALLY the hyphen. 4688 4689 Level: advanced 4690 4691 .keywords: SNES, append, options, prefix, database 4692 4693 .seealso: SNESGetOptionsPrefix() 4694 @*/ 4695 PetscErrorCode SNESAppendOptionsPrefix(SNES snes,const char prefix[]) 4696 { 4697 PetscErrorCode ierr; 4698 4699 PetscFunctionBegin; 4700 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4701 ierr = PetscObjectAppendOptionsPrefix((PetscObject)snes,prefix);CHKERRQ(ierr); 4702 if (!snes->ksp) {ierr = SNESGetKSP(snes,&snes->ksp);CHKERRQ(ierr);} 4703 if (snes->linesearch) { 4704 ierr = SNESGetLineSearch(snes,&snes->linesearch);CHKERRQ(ierr); 4705 ierr = PetscObjectAppendOptionsPrefix((PetscObject)snes->linesearch,prefix);CHKERRQ(ierr); 4706 } 4707 ierr = KSPAppendOptionsPrefix(snes->ksp,prefix);CHKERRQ(ierr); 4708 PetscFunctionReturn(0); 4709 } 4710 4711 /*@C 4712 SNESGetOptionsPrefix - Sets the prefix used for searching for all 4713 SNES options in the database. 4714 4715 Not Collective 4716 4717 Input Parameter: 4718 . snes - the SNES context 4719 4720 Output Parameter: 4721 . prefix - pointer to the prefix string used 4722 4723 Notes: On the fortran side, the user should pass in a string 'prefix' of 4724 sufficient length to hold the prefix. 4725 4726 Level: advanced 4727 4728 .keywords: SNES, get, options, prefix, database 4729 4730 .seealso: SNESAppendOptionsPrefix() 4731 @*/ 4732 PetscErrorCode SNESGetOptionsPrefix(SNES snes,const char *prefix[]) 4733 { 4734 PetscErrorCode ierr; 4735 4736 PetscFunctionBegin; 4737 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4738 ierr = PetscObjectGetOptionsPrefix((PetscObject)snes,prefix);CHKERRQ(ierr); 4739 PetscFunctionReturn(0); 4740 } 4741 4742 4743 /*@C 4744 SNESRegister - Adds a method to the nonlinear solver package. 4745 4746 Not collective 4747 4748 Input Parameters: 4749 + name_solver - name of a new user-defined solver 4750 - routine_create - routine to create method context 4751 4752 Notes: 4753 SNESRegister() may be called multiple times to add several user-defined solvers. 4754 4755 Sample usage: 4756 .vb 4757 SNESRegister("my_solver",MySolverCreate); 4758 .ve 4759 4760 Then, your solver can be chosen with the procedural interface via 4761 $ SNESSetType(snes,"my_solver") 4762 or at runtime via the option 4763 $ -snes_type my_solver 4764 4765 Level: advanced 4766 4767 Note: If your function is not being put into a shared library then use SNESRegister() instead 4768 4769 .keywords: SNES, nonlinear, register 4770 4771 .seealso: SNESRegisterAll(), SNESRegisterDestroy() 4772 4773 Level: advanced 4774 @*/ 4775 PetscErrorCode SNESRegister(const char sname[],PetscErrorCode (*function)(SNES)) 4776 { 4777 PetscErrorCode ierr; 4778 4779 PetscFunctionBegin; 4780 ierr = PetscFunctionListAdd(&SNESList,sname,function);CHKERRQ(ierr); 4781 PetscFunctionReturn(0); 4782 } 4783 4784 PetscErrorCode SNESTestLocalMin(SNES snes) 4785 { 4786 PetscErrorCode ierr; 4787 PetscInt N,i,j; 4788 Vec u,uh,fh; 4789 PetscScalar value; 4790 PetscReal norm; 4791 4792 PetscFunctionBegin; 4793 ierr = SNESGetSolution(snes,&u);CHKERRQ(ierr); 4794 ierr = VecDuplicate(u,&uh);CHKERRQ(ierr); 4795 ierr = VecDuplicate(u,&fh);CHKERRQ(ierr); 4796 4797 /* currently only works for sequential */ 4798 ierr = PetscPrintf(PETSC_COMM_WORLD,"Testing FormFunction() for local min\n");CHKERRQ(ierr); 4799 ierr = VecGetSize(u,&N);CHKERRQ(ierr); 4800 for (i=0; i<N; i++) { 4801 ierr = VecCopy(u,uh);CHKERRQ(ierr); 4802 ierr = PetscPrintf(PETSC_COMM_WORLD,"i = %D\n",i);CHKERRQ(ierr); 4803 for (j=-10; j<11; j++) { 4804 value = PetscSign(j)*PetscExpReal(PetscAbs(j)-10.0); 4805 ierr = VecSetValue(uh,i,value,ADD_VALUES);CHKERRQ(ierr); 4806 ierr = SNESComputeFunction(snes,uh,fh);CHKERRQ(ierr); 4807 ierr = VecNorm(fh,NORM_2,&norm);CHKERRQ(ierr); 4808 ierr = PetscPrintf(PETSC_COMM_WORLD," j norm %D %18.16e\n",j,norm);CHKERRQ(ierr); 4809 value = -value; 4810 ierr = VecSetValue(uh,i,value,ADD_VALUES);CHKERRQ(ierr); 4811 } 4812 } 4813 ierr = VecDestroy(&uh);CHKERRQ(ierr); 4814 ierr = VecDestroy(&fh);CHKERRQ(ierr); 4815 PetscFunctionReturn(0); 4816 } 4817 4818 /*@ 4819 SNESKSPSetUseEW - Sets SNES use Eisenstat-Walker method for 4820 computing relative tolerance for linear solvers within an inexact 4821 Newton method. 4822 4823 Logically Collective on SNES 4824 4825 Input Parameters: 4826 + snes - SNES context 4827 - flag - PETSC_TRUE or PETSC_FALSE 4828 4829 Options Database: 4830 + -snes_ksp_ew - use Eisenstat-Walker method for determining linear system convergence 4831 . -snes_ksp_ew_version ver - version of Eisenstat-Walker method 4832 . -snes_ksp_ew_rtol0 <rtol0> - Sets rtol0 4833 . -snes_ksp_ew_rtolmax <rtolmax> - Sets rtolmax 4834 . -snes_ksp_ew_gamma <gamma> - Sets gamma 4835 . -snes_ksp_ew_alpha <alpha> - Sets alpha 4836 . -snes_ksp_ew_alpha2 <alpha2> - Sets alpha2 4837 - -snes_ksp_ew_threshold <threshold> - Sets threshold 4838 4839 Notes: 4840 Currently, the default is to use a constant relative tolerance for 4841 the inner linear solvers. Alternatively, one can use the 4842 Eisenstat-Walker method, where the relative convergence tolerance 4843 is reset at each Newton iteration according progress of the nonlinear 4844 solver. 4845 4846 Level: advanced 4847 4848 Reference: 4849 S. C. Eisenstat and H. F. Walker, "Choosing the forcing terms in an 4850 inexact Newton method", SISC 17 (1), pp.16-32, 1996. 4851 4852 .keywords: SNES, KSP, Eisenstat, Walker, convergence, test, inexact, Newton 4853 4854 .seealso: SNESKSPGetUseEW(), SNESKSPGetParametersEW(), SNESKSPSetParametersEW() 4855 @*/ 4856 PetscErrorCode SNESKSPSetUseEW(SNES snes,PetscBool flag) 4857 { 4858 PetscFunctionBegin; 4859 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4860 PetscValidLogicalCollectiveBool(snes,flag,2); 4861 snes->ksp_ewconv = flag; 4862 PetscFunctionReturn(0); 4863 } 4864 4865 /*@ 4866 SNESKSPGetUseEW - Gets if SNES is using Eisenstat-Walker method 4867 for computing relative tolerance for linear solvers within an 4868 inexact Newton method. 4869 4870 Not Collective 4871 4872 Input Parameter: 4873 . snes - SNES context 4874 4875 Output Parameter: 4876 . flag - PETSC_TRUE or PETSC_FALSE 4877 4878 Notes: 4879 Currently, the default is to use a constant relative tolerance for 4880 the inner linear solvers. Alternatively, one can use the 4881 Eisenstat-Walker method, where the relative convergence tolerance 4882 is reset at each Newton iteration according progress of the nonlinear 4883 solver. 4884 4885 Level: advanced 4886 4887 Reference: 4888 S. C. Eisenstat and H. F. Walker, "Choosing the forcing terms in an 4889 inexact Newton method", SISC 17 (1), pp.16-32, 1996. 4890 4891 .keywords: SNES, KSP, Eisenstat, Walker, convergence, test, inexact, Newton 4892 4893 .seealso: SNESKSPSetUseEW(), SNESKSPGetParametersEW(), SNESKSPSetParametersEW() 4894 @*/ 4895 PetscErrorCode SNESKSPGetUseEW(SNES snes, PetscBool *flag) 4896 { 4897 PetscFunctionBegin; 4898 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4899 PetscValidPointer(flag,2); 4900 *flag = snes->ksp_ewconv; 4901 PetscFunctionReturn(0); 4902 } 4903 4904 /*@ 4905 SNESKSPSetParametersEW - Sets parameters for Eisenstat-Walker 4906 convergence criteria for the linear solvers within an inexact 4907 Newton method. 4908 4909 Logically Collective on SNES 4910 4911 Input Parameters: 4912 + snes - SNES context 4913 . version - version 1, 2 (default is 2) or 3 4914 . rtol_0 - initial relative tolerance (0 <= rtol_0 < 1) 4915 . rtol_max - maximum relative tolerance (0 <= rtol_max < 1) 4916 . gamma - multiplicative factor for version 2 rtol computation 4917 (0 <= gamma2 <= 1) 4918 . alpha - power for version 2 rtol computation (1 < alpha <= 2) 4919 . alpha2 - power for safeguard 4920 - threshold - threshold for imposing safeguard (0 < threshold < 1) 4921 4922 Note: 4923 Version 3 was contributed by Luis Chacon, June 2006. 4924 4925 Use PETSC_DEFAULT to retain the default for any of the parameters. 4926 4927 Level: advanced 4928 4929 Reference: 4930 S. C. Eisenstat and H. F. Walker, "Choosing the forcing terms in an 4931 inexact Newton method", Utah State University Math. Stat. Dept. Res. 4932 Report 6/94/75, June, 1994, to appear in SIAM J. Sci. Comput. 4933 4934 .keywords: SNES, KSP, Eisenstat, Walker, set, parameters 4935 4936 .seealso: SNESKSPSetUseEW(), SNESKSPGetUseEW(), SNESKSPGetParametersEW() 4937 @*/ 4938 PetscErrorCode SNESKSPSetParametersEW(SNES snes,PetscInt version,PetscReal rtol_0,PetscReal rtol_max,PetscReal gamma,PetscReal alpha,PetscReal alpha2,PetscReal threshold) 4939 { 4940 SNESKSPEW *kctx; 4941 4942 PetscFunctionBegin; 4943 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4944 kctx = (SNESKSPEW*)snes->kspconvctx; 4945 if (!kctx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"No Eisenstat-Walker context existing"); 4946 PetscValidLogicalCollectiveInt(snes,version,2); 4947 PetscValidLogicalCollectiveReal(snes,rtol_0,3); 4948 PetscValidLogicalCollectiveReal(snes,rtol_max,4); 4949 PetscValidLogicalCollectiveReal(snes,gamma,5); 4950 PetscValidLogicalCollectiveReal(snes,alpha,6); 4951 PetscValidLogicalCollectiveReal(snes,alpha2,7); 4952 PetscValidLogicalCollectiveReal(snes,threshold,8); 4953 4954 if (version != PETSC_DEFAULT) kctx->version = version; 4955 if (rtol_0 != PETSC_DEFAULT) kctx->rtol_0 = rtol_0; 4956 if (rtol_max != PETSC_DEFAULT) kctx->rtol_max = rtol_max; 4957 if (gamma != PETSC_DEFAULT) kctx->gamma = gamma; 4958 if (alpha != PETSC_DEFAULT) kctx->alpha = alpha; 4959 if (alpha2 != PETSC_DEFAULT) kctx->alpha2 = alpha2; 4960 if (threshold != PETSC_DEFAULT) kctx->threshold = threshold; 4961 4962 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); 4963 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); 4964 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); 4965 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); 4966 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); 4967 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); 4968 PetscFunctionReturn(0); 4969 } 4970 4971 /*@ 4972 SNESKSPGetParametersEW - Gets parameters for Eisenstat-Walker 4973 convergence criteria for the linear solvers within an inexact 4974 Newton method. 4975 4976 Not Collective 4977 4978 Input Parameters: 4979 snes - SNES context 4980 4981 Output Parameters: 4982 + version - version 1, 2 (default is 2) or 3 4983 . rtol_0 - initial relative tolerance (0 <= rtol_0 < 1) 4984 . rtol_max - maximum relative tolerance (0 <= rtol_max < 1) 4985 . gamma - multiplicative factor for version 2 rtol computation (0 <= gamma2 <= 1) 4986 . alpha - power for version 2 rtol computation (1 < alpha <= 2) 4987 . alpha2 - power for safeguard 4988 - threshold - threshold for imposing safeguard (0 < threshold < 1) 4989 4990 Level: advanced 4991 4992 .keywords: SNES, KSP, Eisenstat, Walker, get, parameters 4993 4994 .seealso: SNESKSPSetUseEW(), SNESKSPGetUseEW(), SNESKSPSetParametersEW() 4995 @*/ 4996 PetscErrorCode SNESKSPGetParametersEW(SNES snes,PetscInt *version,PetscReal *rtol_0,PetscReal *rtol_max,PetscReal *gamma,PetscReal *alpha,PetscReal *alpha2,PetscReal *threshold) 4997 { 4998 SNESKSPEW *kctx; 4999 5000 PetscFunctionBegin; 5001 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 5002 kctx = (SNESKSPEW*)snes->kspconvctx; 5003 if (!kctx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"No Eisenstat-Walker context existing"); 5004 if (version) *version = kctx->version; 5005 if (rtol_0) *rtol_0 = kctx->rtol_0; 5006 if (rtol_max) *rtol_max = kctx->rtol_max; 5007 if (gamma) *gamma = kctx->gamma; 5008 if (alpha) *alpha = kctx->alpha; 5009 if (alpha2) *alpha2 = kctx->alpha2; 5010 if (threshold) *threshold = kctx->threshold; 5011 PetscFunctionReturn(0); 5012 } 5013 5014 PetscErrorCode KSPPreSolve_SNESEW(KSP ksp, Vec b, Vec x, SNES snes) 5015 { 5016 PetscErrorCode ierr; 5017 SNESKSPEW *kctx = (SNESKSPEW*)snes->kspconvctx; 5018 PetscReal rtol = PETSC_DEFAULT,stol; 5019 5020 PetscFunctionBegin; 5021 if (!snes->ksp_ewconv) PetscFunctionReturn(0); 5022 if (!snes->iter) { 5023 rtol = kctx->rtol_0; /* first time in, so use the original user rtol */ 5024 ierr = VecNorm(snes->vec_func,NORM_2,&kctx->norm_first);CHKERRQ(ierr); 5025 } 5026 else { 5027 if (kctx->version == 1) { 5028 rtol = (snes->norm - kctx->lresid_last)/kctx->norm_last; 5029 if (rtol < 0.0) rtol = -rtol; 5030 stol = PetscPowReal(kctx->rtol_last,kctx->alpha2); 5031 if (stol > kctx->threshold) rtol = PetscMax(rtol,stol); 5032 } else if (kctx->version == 2) { 5033 rtol = kctx->gamma * PetscPowReal(snes->norm/kctx->norm_last,kctx->alpha); 5034 stol = kctx->gamma * PetscPowReal(kctx->rtol_last,kctx->alpha); 5035 if (stol > kctx->threshold) rtol = PetscMax(rtol,stol); 5036 } else if (kctx->version == 3) { /* contributed by Luis Chacon, June 2006. */ 5037 rtol = kctx->gamma * PetscPowReal(snes->norm/kctx->norm_last,kctx->alpha); 5038 /* safeguard: avoid sharp decrease of rtol */ 5039 stol = kctx->gamma*PetscPowReal(kctx->rtol_last,kctx->alpha); 5040 stol = PetscMax(rtol,stol); 5041 rtol = PetscMin(kctx->rtol_0,stol); 5042 /* safeguard: avoid oversolving */ 5043 stol = kctx->gamma*(kctx->norm_first*snes->rtol)/snes->norm; 5044 stol = PetscMax(rtol,stol); 5045 rtol = PetscMin(kctx->rtol_0,stol); 5046 } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Only versions 1, 2 or 3 are supported: %D",kctx->version); 5047 } 5048 /* safeguard: avoid rtol greater than one */ 5049 rtol = PetscMin(rtol,kctx->rtol_max); 5050 ierr = KSPSetTolerances(ksp,rtol,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);CHKERRQ(ierr); 5051 ierr = PetscInfo3(snes,"iter %D, Eisenstat-Walker (version %D) KSP rtol=%g\n",snes->iter,kctx->version,(double)rtol);CHKERRQ(ierr); 5052 PetscFunctionReturn(0); 5053 } 5054 5055 PetscErrorCode KSPPostSolve_SNESEW(KSP ksp, Vec b, Vec x, SNES snes) 5056 { 5057 PetscErrorCode ierr; 5058 SNESKSPEW *kctx = (SNESKSPEW*)snes->kspconvctx; 5059 PCSide pcside; 5060 Vec lres; 5061 5062 PetscFunctionBegin; 5063 if (!snes->ksp_ewconv) PetscFunctionReturn(0); 5064 ierr = KSPGetTolerances(ksp,&kctx->rtol_last,0,0,0);CHKERRQ(ierr); 5065 kctx->norm_last = snes->norm; 5066 if (kctx->version == 1) { 5067 PC pc; 5068 PetscBool isNone; 5069 5070 ierr = KSPGetPC(ksp, &pc);CHKERRQ(ierr); 5071 ierr = PetscObjectTypeCompare((PetscObject) pc, PCNONE, &isNone);CHKERRQ(ierr); 5072 ierr = KSPGetPCSide(ksp,&pcside);CHKERRQ(ierr); 5073 if (pcside == PC_RIGHT || isNone) { /* XXX Should we also test KSP_UNPRECONDITIONED_NORM ? */ 5074 /* KSP residual is true linear residual */ 5075 ierr = KSPGetResidualNorm(ksp,&kctx->lresid_last);CHKERRQ(ierr); 5076 } else { 5077 /* KSP residual is preconditioned residual */ 5078 /* compute true linear residual norm */ 5079 ierr = VecDuplicate(b,&lres);CHKERRQ(ierr); 5080 ierr = MatMult(snes->jacobian,x,lres);CHKERRQ(ierr); 5081 ierr = VecAYPX(lres,-1.0,b);CHKERRQ(ierr); 5082 ierr = VecNorm(lres,NORM_2,&kctx->lresid_last);CHKERRQ(ierr); 5083 ierr = VecDestroy(&lres);CHKERRQ(ierr); 5084 } 5085 } 5086 PetscFunctionReturn(0); 5087 } 5088 5089 /*@ 5090 SNESGetKSP - Returns the KSP context for a SNES solver. 5091 5092 Not Collective, but if SNES object is parallel, then KSP object is parallel 5093 5094 Input Parameter: 5095 . snes - the SNES context 5096 5097 Output Parameter: 5098 . ksp - the KSP context 5099 5100 Notes: 5101 The user can then directly manipulate the KSP context to set various 5102 options, etc. Likewise, the user can then extract and manipulate the 5103 PC contexts as well. 5104 5105 Level: beginner 5106 5107 .keywords: SNES, nonlinear, get, KSP, context 5108 5109 .seealso: KSPGetPC(), SNESCreate(), KSPCreate(), SNESSetKSP() 5110 @*/ 5111 PetscErrorCode SNESGetKSP(SNES snes,KSP *ksp) 5112 { 5113 PetscErrorCode ierr; 5114 5115 PetscFunctionBegin; 5116 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 5117 PetscValidPointer(ksp,2); 5118 5119 if (!snes->ksp) { 5120 PetscBool monitor = PETSC_FALSE; 5121 5122 ierr = KSPCreate(PetscObjectComm((PetscObject)snes),&snes->ksp);CHKERRQ(ierr); 5123 ierr = PetscObjectIncrementTabLevel((PetscObject)snes->ksp,(PetscObject)snes,1);CHKERRQ(ierr); 5124 ierr = PetscLogObjectParent((PetscObject)snes,(PetscObject)snes->ksp);CHKERRQ(ierr); 5125 5126 ierr = KSPSetPreSolve(snes->ksp,(PetscErrorCode (*)(KSP,Vec,Vec,void*))KSPPreSolve_SNESEW,snes);CHKERRQ(ierr); 5127 ierr = KSPSetPostSolve(snes->ksp,(PetscErrorCode (*)(KSP,Vec,Vec,void*))KSPPostSolve_SNESEW,snes);CHKERRQ(ierr); 5128 5129 ierr = PetscOptionsGetBool(((PetscObject)snes)->options,((PetscObject)snes)->prefix,"-ksp_monitor_snes",&monitor,NULL);CHKERRQ(ierr); 5130 if (monitor) { 5131 ierr = KSPMonitorSet(snes->ksp,KSPMonitorSNES,snes,NULL);CHKERRQ(ierr); 5132 } 5133 monitor = PETSC_FALSE; 5134 ierr = PetscOptionsGetBool(((PetscObject)snes)->options,((PetscObject)snes)->prefix,"-ksp_monitor_snes_lg",&monitor,NULL);CHKERRQ(ierr); 5135 if (monitor) { 5136 PetscObject *objs; 5137 ierr = KSPMonitorSNESLGResidualNormCreate(PetscObjectComm((PetscObject)snes),NULL,NULL,PETSC_DECIDE,PETSC_DECIDE,600,600,&objs);CHKERRQ(ierr); 5138 objs[0] = (PetscObject) snes; 5139 ierr = KSPMonitorSet(snes->ksp,(PetscErrorCode (*)(KSP,PetscInt,PetscReal,void*))KSPMonitorSNESLGResidualNorm,objs,(PetscErrorCode (*)(void**))KSPMonitorSNESLGResidualNormDestroy);CHKERRQ(ierr); 5140 } 5141 } 5142 *ksp = snes->ksp; 5143 PetscFunctionReturn(0); 5144 } 5145 5146 5147 #include <petsc/private/dmimpl.h> 5148 /*@ 5149 SNESSetDM - Sets the DM that may be used by some nonlinear solvers or their underlying preconditioners 5150 5151 Logically Collective on SNES 5152 5153 Input Parameters: 5154 + snes - the nonlinear solver context 5155 - dm - the dm, cannot be NULL 5156 5157 Level: intermediate 5158 5159 .seealso: SNESGetDM(), KSPSetDM(), KSPGetDM() 5160 @*/ 5161 PetscErrorCode SNESSetDM(SNES snes,DM dm) 5162 { 5163 PetscErrorCode ierr; 5164 KSP ksp; 5165 DMSNES sdm; 5166 5167 PetscFunctionBegin; 5168 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 5169 PetscValidHeaderSpecific(dm,DM_CLASSID,2); 5170 ierr = PetscObjectReference((PetscObject)dm);CHKERRQ(ierr); 5171 if (snes->dm) { /* Move the DMSNES context over to the new DM unless the new DM already has one */ 5172 if (snes->dm->dmsnes && !dm->dmsnes) { 5173 ierr = DMCopyDMSNES(snes->dm,dm);CHKERRQ(ierr); 5174 ierr = DMGetDMSNES(snes->dm,&sdm);CHKERRQ(ierr); 5175 if (sdm->originaldm == snes->dm) sdm->originaldm = dm; /* Grant write privileges to the replacement DM */ 5176 } 5177 ierr = DMCoarsenHookRemove(snes->dm,DMCoarsenHook_SNESVecSol,DMRestrictHook_SNESVecSol,snes);CHKERRQ(ierr); 5178 ierr = DMDestroy(&snes->dm);CHKERRQ(ierr); 5179 } 5180 snes->dm = dm; 5181 snes->dmAuto = PETSC_FALSE; 5182 5183 ierr = SNESGetKSP(snes,&ksp);CHKERRQ(ierr); 5184 ierr = KSPSetDM(ksp,dm);CHKERRQ(ierr); 5185 ierr = KSPSetDMActive(ksp,PETSC_FALSE);CHKERRQ(ierr); 5186 if (snes->npc) { 5187 ierr = SNESSetDM(snes->npc, snes->dm);CHKERRQ(ierr); 5188 ierr = SNESSetNPCSide(snes,snes->npcside);CHKERRQ(ierr); 5189 } 5190 PetscFunctionReturn(0); 5191 } 5192 5193 /*@ 5194 SNESGetDM - Gets the DM that may be used by some preconditioners 5195 5196 Not Collective but DM obtained is parallel on SNES 5197 5198 Input Parameter: 5199 . snes - the preconditioner context 5200 5201 Output Parameter: 5202 . dm - the dm 5203 5204 Level: intermediate 5205 5206 .seealso: SNESSetDM(), KSPSetDM(), KSPGetDM() 5207 @*/ 5208 PetscErrorCode SNESGetDM(SNES snes,DM *dm) 5209 { 5210 PetscErrorCode ierr; 5211 5212 PetscFunctionBegin; 5213 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 5214 if (!snes->dm) { 5215 ierr = DMShellCreate(PetscObjectComm((PetscObject)snes),&snes->dm);CHKERRQ(ierr); 5216 snes->dmAuto = PETSC_TRUE; 5217 } 5218 *dm = snes->dm; 5219 PetscFunctionReturn(0); 5220 } 5221 5222 /*@ 5223 SNESSetNPC - Sets the nonlinear preconditioner to be used. 5224 5225 Collective on SNES 5226 5227 Input Parameters: 5228 + snes - iterative context obtained from SNESCreate() 5229 - pc - the preconditioner object 5230 5231 Notes: 5232 Use SNESGetNPC() to retrieve the preconditioner context (for example, 5233 to configure it using the API). 5234 5235 Level: developer 5236 5237 .keywords: SNES, set, precondition 5238 .seealso: SNESGetNPC(), SNESHasNPC() 5239 @*/ 5240 PetscErrorCode SNESSetNPC(SNES snes, SNES pc) 5241 { 5242 PetscErrorCode ierr; 5243 5244 PetscFunctionBegin; 5245 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); 5246 PetscValidHeaderSpecific(pc, SNES_CLASSID, 2); 5247 PetscCheckSameComm(snes, 1, pc, 2); 5248 ierr = PetscObjectReference((PetscObject) pc);CHKERRQ(ierr); 5249 ierr = SNESDestroy(&snes->npc);CHKERRQ(ierr); 5250 snes->npc = pc; 5251 ierr = PetscLogObjectParent((PetscObject)snes, (PetscObject)snes->npc);CHKERRQ(ierr); 5252 PetscFunctionReturn(0); 5253 } 5254 5255 /*@ 5256 SNESGetNPC - Creates a nonlinear preconditioning solver (SNES) to be used to precondition the nonlinear solver. 5257 5258 Not Collective 5259 5260 Input Parameter: 5261 . snes - iterative context obtained from SNESCreate() 5262 5263 Output Parameter: 5264 . pc - preconditioner context 5265 5266 Notes: If a SNES was previously set with SNESSetNPC() then that SNES is returned. 5267 5268 Level: developer 5269 5270 .keywords: SNES, get, preconditioner 5271 .seealso: SNESSetNPC(), SNESHasNPC() 5272 @*/ 5273 PetscErrorCode SNESGetNPC(SNES snes, SNES *pc) 5274 { 5275 PetscErrorCode ierr; 5276 const char *optionsprefix; 5277 5278 PetscFunctionBegin; 5279 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); 5280 PetscValidPointer(pc, 2); 5281 if (!snes->npc) { 5282 ierr = SNESCreate(PetscObjectComm((PetscObject)snes),&snes->npc);CHKERRQ(ierr); 5283 ierr = PetscObjectIncrementTabLevel((PetscObject)snes->npc,(PetscObject)snes,1);CHKERRQ(ierr); 5284 ierr = PetscLogObjectParent((PetscObject)snes,(PetscObject)snes->npc);CHKERRQ(ierr); 5285 ierr = SNESGetOptionsPrefix(snes,&optionsprefix);CHKERRQ(ierr); 5286 ierr = SNESSetOptionsPrefix(snes->npc,optionsprefix);CHKERRQ(ierr); 5287 ierr = SNESAppendOptionsPrefix(snes->npc,"npc_");CHKERRQ(ierr); 5288 ierr = SNESSetCountersReset(snes->npc,PETSC_FALSE);CHKERRQ(ierr); 5289 } 5290 *pc = snes->npc; 5291 PetscFunctionReturn(0); 5292 } 5293 5294 /*@ 5295 SNESHasNPC - Returns whether a nonlinear preconditioner exists 5296 5297 Not Collective 5298 5299 Input Parameter: 5300 . snes - iterative context obtained from SNESCreate() 5301 5302 Output Parameter: 5303 . has_npc - whether the SNES has an NPC or not 5304 5305 Level: developer 5306 5307 .keywords: SNES, has, preconditioner 5308 .seealso: SNESSetNPC(), SNESGetNPC() 5309 @*/ 5310 PetscErrorCode SNESHasNPC(SNES snes, PetscBool *has_npc) 5311 { 5312 PetscFunctionBegin; 5313 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); 5314 *has_npc = (PetscBool) (snes->npc ? PETSC_TRUE : PETSC_FALSE); 5315 PetscFunctionReturn(0); 5316 } 5317 5318 /*@ 5319 SNESSetNPCSide - Sets the preconditioning side. 5320 5321 Logically Collective on SNES 5322 5323 Input Parameter: 5324 . snes - iterative context obtained from SNESCreate() 5325 5326 Output Parameter: 5327 . side - the preconditioning side, where side is one of 5328 .vb 5329 PC_LEFT - left preconditioning 5330 PC_RIGHT - right preconditioning (default for most nonlinear solvers) 5331 .ve 5332 5333 Options Database Keys: 5334 . -snes_pc_side <right,left> 5335 5336 Notes: SNESNRICHARDSON and SNESNCG only support left preconditioning. 5337 5338 Level: intermediate 5339 5340 .keywords: SNES, set, right, left, side, preconditioner, flag 5341 5342 .seealso: SNESGetNPCSide(), KSPSetPCSide() 5343 @*/ 5344 PetscErrorCode SNESSetNPCSide(SNES snes,PCSide side) 5345 { 5346 PetscFunctionBegin; 5347 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 5348 PetscValidLogicalCollectiveEnum(snes,side,2); 5349 snes->npcside= side; 5350 PetscFunctionReturn(0); 5351 } 5352 5353 /*@ 5354 SNESGetNPCSide - Gets the preconditioning side. 5355 5356 Not Collective 5357 5358 Input Parameter: 5359 . snes - iterative context obtained from SNESCreate() 5360 5361 Output Parameter: 5362 . side - the preconditioning side, where side is one of 5363 .vb 5364 PC_LEFT - left preconditioning 5365 PC_RIGHT - right preconditioning (default for most nonlinear solvers) 5366 .ve 5367 5368 Level: intermediate 5369 5370 .keywords: SNES, get, right, left, side, preconditioner, flag 5371 5372 .seealso: SNESSetNPCSide(), KSPGetPCSide() 5373 @*/ 5374 PetscErrorCode SNESGetNPCSide(SNES snes,PCSide *side) 5375 { 5376 PetscFunctionBegin; 5377 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 5378 PetscValidPointer(side,2); 5379 *side = snes->npcside; 5380 PetscFunctionReturn(0); 5381 } 5382 5383 /*@ 5384 SNESSetLineSearch - Sets the linesearch on the SNES instance. 5385 5386 Collective on SNES 5387 5388 Input Parameters: 5389 + snes - iterative context obtained from SNESCreate() 5390 - linesearch - the linesearch object 5391 5392 Notes: 5393 Use SNESGetLineSearch() to retrieve the preconditioner context (for example, 5394 to configure it using the API). 5395 5396 Level: developer 5397 5398 .keywords: SNES, set, linesearch 5399 .seealso: SNESGetLineSearch() 5400 @*/ 5401 PetscErrorCode SNESSetLineSearch(SNES snes, SNESLineSearch linesearch) 5402 { 5403 PetscErrorCode ierr; 5404 5405 PetscFunctionBegin; 5406 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); 5407 PetscValidHeaderSpecific(linesearch, SNESLINESEARCH_CLASSID, 2); 5408 PetscCheckSameComm(snes, 1, linesearch, 2); 5409 ierr = PetscObjectReference((PetscObject) linesearch);CHKERRQ(ierr); 5410 ierr = SNESLineSearchDestroy(&snes->linesearch);CHKERRQ(ierr); 5411 5412 snes->linesearch = linesearch; 5413 5414 ierr = PetscLogObjectParent((PetscObject)snes, (PetscObject)snes->linesearch);CHKERRQ(ierr); 5415 PetscFunctionReturn(0); 5416 } 5417 5418 /*@ 5419 SNESGetLineSearch - Returns a pointer to the line search context set with SNESSetLineSearch() 5420 or creates a default line search instance associated with the SNES and returns it. 5421 5422 Not Collective 5423 5424 Input Parameter: 5425 . snes - iterative context obtained from SNESCreate() 5426 5427 Output Parameter: 5428 . linesearch - linesearch context 5429 5430 Level: beginner 5431 5432 .keywords: SNES, get, linesearch 5433 .seealso: SNESSetLineSearch(), SNESLineSearchCreate() 5434 @*/ 5435 PetscErrorCode SNESGetLineSearch(SNES snes, SNESLineSearch *linesearch) 5436 { 5437 PetscErrorCode ierr; 5438 const char *optionsprefix; 5439 5440 PetscFunctionBegin; 5441 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); 5442 PetscValidPointer(linesearch, 2); 5443 if (!snes->linesearch) { 5444 ierr = SNESGetOptionsPrefix(snes, &optionsprefix);CHKERRQ(ierr); 5445 ierr = SNESLineSearchCreate(PetscObjectComm((PetscObject)snes), &snes->linesearch);CHKERRQ(ierr); 5446 ierr = SNESLineSearchSetSNES(snes->linesearch, snes);CHKERRQ(ierr); 5447 ierr = SNESLineSearchAppendOptionsPrefix(snes->linesearch, optionsprefix);CHKERRQ(ierr); 5448 ierr = PetscObjectIncrementTabLevel((PetscObject) snes->linesearch, (PetscObject) snes, 1);CHKERRQ(ierr); 5449 ierr = PetscLogObjectParent((PetscObject)snes, (PetscObject)snes->linesearch);CHKERRQ(ierr); 5450 } 5451 *linesearch = snes->linesearch; 5452 PetscFunctionReturn(0); 5453 } 5454 5455 #if defined(PETSC_HAVE_MATLAB_ENGINE) 5456 #include <mex.h> 5457 5458 typedef struct {char *funcname; mxArray *ctx;} SNESMatlabContext; 5459 5460 /* 5461 SNESComputeFunction_Matlab - Calls the function that has been set with SNESSetFunctionMatlab(). 5462 5463 Collective on SNES 5464 5465 Input Parameters: 5466 + snes - the SNES context 5467 - x - input vector 5468 5469 Output Parameter: 5470 . y - function vector, as set by SNESSetFunction() 5471 5472 Notes: 5473 SNESComputeFunction() is typically used within nonlinear solvers 5474 implementations, so most users would not generally call this routine 5475 themselves. 5476 5477 Level: developer 5478 5479 .keywords: SNES, nonlinear, compute, function 5480 5481 .seealso: SNESSetFunction(), SNESGetFunction() 5482 */ 5483 PetscErrorCode SNESComputeFunction_Matlab(SNES snes,Vec x,Vec y, void *ctx) 5484 { 5485 PetscErrorCode ierr; 5486 SNESMatlabContext *sctx = (SNESMatlabContext*)ctx; 5487 int nlhs = 1,nrhs = 5; 5488 mxArray *plhs[1],*prhs[5]; 5489 long long int lx = 0,ly = 0,ls = 0; 5490 5491 PetscFunctionBegin; 5492 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 5493 PetscValidHeaderSpecific(x,VEC_CLASSID,2); 5494 PetscValidHeaderSpecific(y,VEC_CLASSID,3); 5495 PetscCheckSameComm(snes,1,x,2); 5496 PetscCheckSameComm(snes,1,y,3); 5497 5498 /* call Matlab function in ctx with arguments x and y */ 5499 5500 ierr = PetscMemcpy(&ls,&snes,sizeof(snes));CHKERRQ(ierr); 5501 ierr = PetscMemcpy(&lx,&x,sizeof(x));CHKERRQ(ierr); 5502 ierr = PetscMemcpy(&ly,&y,sizeof(x));CHKERRQ(ierr); 5503 prhs[0] = mxCreateDoubleScalar((double)ls); 5504 prhs[1] = mxCreateDoubleScalar((double)lx); 5505 prhs[2] = mxCreateDoubleScalar((double)ly); 5506 prhs[3] = mxCreateString(sctx->funcname); 5507 prhs[4] = sctx->ctx; 5508 ierr = mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscSNESComputeFunctionInternal");CHKERRQ(ierr); 5509 ierr = mxGetScalar(plhs[0]);CHKERRQ(ierr); 5510 mxDestroyArray(prhs[0]); 5511 mxDestroyArray(prhs[1]); 5512 mxDestroyArray(prhs[2]); 5513 mxDestroyArray(prhs[3]); 5514 mxDestroyArray(plhs[0]); 5515 PetscFunctionReturn(0); 5516 } 5517 5518 /* 5519 SNESSetFunctionMatlab - Sets the function evaluation routine and function 5520 vector for use by the SNES routines in solving systems of nonlinear 5521 equations from MATLAB. Here the function is a string containing the name of a MATLAB function 5522 5523 Logically Collective on SNES 5524 5525 Input Parameters: 5526 + snes - the SNES context 5527 . r - vector to store function value 5528 - f - function evaluation routine 5529 5530 Notes: 5531 The Newton-like methods typically solve linear systems of the form 5532 $ f'(x) x = -f(x), 5533 where f'(x) denotes the Jacobian matrix and f(x) is the function. 5534 5535 Level: beginner 5536 5537 Developer Note: This bleeds the allocated memory SNESMatlabContext *sctx; 5538 5539 .keywords: SNES, nonlinear, set, function 5540 5541 .seealso: SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESSetFunction() 5542 */ 5543 PetscErrorCode SNESSetFunctionMatlab(SNES snes,Vec r,const char *f,mxArray *ctx) 5544 { 5545 PetscErrorCode ierr; 5546 SNESMatlabContext *sctx; 5547 5548 PetscFunctionBegin; 5549 /* currently sctx is memory bleed */ 5550 ierr = PetscNew(&sctx);CHKERRQ(ierr); 5551 ierr = PetscStrallocpy(f,&sctx->funcname);CHKERRQ(ierr); 5552 /* 5553 This should work, but it doesn't 5554 sctx->ctx = ctx; 5555 mexMakeArrayPersistent(sctx->ctx); 5556 */ 5557 sctx->ctx = mxDuplicateArray(ctx); 5558 ierr = SNESSetFunction(snes,r,SNESComputeFunction_Matlab,sctx);CHKERRQ(ierr); 5559 PetscFunctionReturn(0); 5560 } 5561 5562 /* 5563 SNESComputeJacobian_Matlab - Calls the function that has been set with SNESSetJacobianMatlab(). 5564 5565 Collective on SNES 5566 5567 Input Parameters: 5568 + snes - the SNES context 5569 . x - input vector 5570 . A, B - the matrices 5571 - ctx - user context 5572 5573 Level: developer 5574 5575 .keywords: SNES, nonlinear, compute, function 5576 5577 .seealso: SNESSetFunction(), SNESGetFunction() 5578 @*/ 5579 PetscErrorCode SNESComputeJacobian_Matlab(SNES snes,Vec x,Mat A,Mat B,void *ctx) 5580 { 5581 PetscErrorCode ierr; 5582 SNESMatlabContext *sctx = (SNESMatlabContext*)ctx; 5583 int nlhs = 2,nrhs = 6; 5584 mxArray *plhs[2],*prhs[6]; 5585 long long int lx = 0,lA = 0,ls = 0, lB = 0; 5586 5587 PetscFunctionBegin; 5588 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 5589 PetscValidHeaderSpecific(x,VEC_CLASSID,2); 5590 5591 /* call Matlab function in ctx with arguments x and y */ 5592 5593 ierr = PetscMemcpy(&ls,&snes,sizeof(snes));CHKERRQ(ierr); 5594 ierr = PetscMemcpy(&lx,&x,sizeof(x));CHKERRQ(ierr); 5595 ierr = PetscMemcpy(&lA,A,sizeof(x));CHKERRQ(ierr); 5596 ierr = PetscMemcpy(&lB,B,sizeof(x));CHKERRQ(ierr); 5597 prhs[0] = mxCreateDoubleScalar((double)ls); 5598 prhs[1] = mxCreateDoubleScalar((double)lx); 5599 prhs[2] = mxCreateDoubleScalar((double)lA); 5600 prhs[3] = mxCreateDoubleScalar((double)lB); 5601 prhs[4] = mxCreateString(sctx->funcname); 5602 prhs[5] = sctx->ctx; 5603 ierr = mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscSNESComputeJacobianInternal");CHKERRQ(ierr); 5604 ierr = mxGetScalar(plhs[0]);CHKERRQ(ierr); 5605 mxDestroyArray(prhs[0]); 5606 mxDestroyArray(prhs[1]); 5607 mxDestroyArray(prhs[2]); 5608 mxDestroyArray(prhs[3]); 5609 mxDestroyArray(prhs[4]); 5610 mxDestroyArray(plhs[0]); 5611 mxDestroyArray(plhs[1]); 5612 PetscFunctionReturn(0); 5613 } 5614 5615 /* 5616 SNESSetJacobianMatlab - Sets the Jacobian function evaluation routine and two empty Jacobian matrices 5617 vector for use by the SNES routines in solving systems of nonlinear 5618 equations from MATLAB. Here the function is a string containing the name of a MATLAB function 5619 5620 Logically Collective on SNES 5621 5622 Input Parameters: 5623 + snes - the SNES context 5624 . A,B - Jacobian matrices 5625 . J - function evaluation routine 5626 - ctx - user context 5627 5628 Level: developer 5629 5630 Developer Note: This bleeds the allocated memory SNESMatlabContext *sctx; 5631 5632 .keywords: SNES, nonlinear, set, function 5633 5634 .seealso: SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESSetFunction(), J 5635 */ 5636 PetscErrorCode SNESSetJacobianMatlab(SNES snes,Mat A,Mat B,const char *J,mxArray *ctx) 5637 { 5638 PetscErrorCode ierr; 5639 SNESMatlabContext *sctx; 5640 5641 PetscFunctionBegin; 5642 /* currently sctx is memory bleed */ 5643 ierr = PetscNew(&sctx);CHKERRQ(ierr); 5644 ierr = PetscStrallocpy(J,&sctx->funcname);CHKERRQ(ierr); 5645 /* 5646 This should work, but it doesn't 5647 sctx->ctx = ctx; 5648 mexMakeArrayPersistent(sctx->ctx); 5649 */ 5650 sctx->ctx = mxDuplicateArray(ctx); 5651 ierr = SNESSetJacobian(snes,A,B,SNESComputeJacobian_Matlab,sctx);CHKERRQ(ierr); 5652 PetscFunctionReturn(0); 5653 } 5654 5655 /* 5656 SNESMonitor_Matlab - Calls the function that has been set with SNESMonitorSetMatlab(). 5657 5658 Collective on SNES 5659 5660 .seealso: SNESSetFunction(), SNESGetFunction() 5661 @*/ 5662 PetscErrorCode SNESMonitor_Matlab(SNES snes,PetscInt it, PetscReal fnorm, void *ctx) 5663 { 5664 PetscErrorCode ierr; 5665 SNESMatlabContext *sctx = (SNESMatlabContext*)ctx; 5666 int nlhs = 1,nrhs = 6; 5667 mxArray *plhs[1],*prhs[6]; 5668 long long int lx = 0,ls = 0; 5669 Vec x = snes->vec_sol; 5670 5671 PetscFunctionBegin; 5672 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 5673 5674 ierr = PetscMemcpy(&ls,&snes,sizeof(snes));CHKERRQ(ierr); 5675 ierr = PetscMemcpy(&lx,&x,sizeof(x));CHKERRQ(ierr); 5676 prhs[0] = mxCreateDoubleScalar((double)ls); 5677 prhs[1] = mxCreateDoubleScalar((double)it); 5678 prhs[2] = mxCreateDoubleScalar((double)fnorm); 5679 prhs[3] = mxCreateDoubleScalar((double)lx); 5680 prhs[4] = mxCreateString(sctx->funcname); 5681 prhs[5] = sctx->ctx; 5682 ierr = mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscSNESMonitorInternal");CHKERRQ(ierr); 5683 ierr = mxGetScalar(plhs[0]);CHKERRQ(ierr); 5684 mxDestroyArray(prhs[0]); 5685 mxDestroyArray(prhs[1]); 5686 mxDestroyArray(prhs[2]); 5687 mxDestroyArray(prhs[3]); 5688 mxDestroyArray(prhs[4]); 5689 mxDestroyArray(plhs[0]); 5690 PetscFunctionReturn(0); 5691 } 5692 5693 /* 5694 SNESMonitorSetMatlab - Sets the monitor function from MATLAB 5695 5696 Level: developer 5697 5698 Developer Note: This bleeds the allocated memory SNESMatlabContext *sctx; 5699 5700 .keywords: SNES, nonlinear, set, function 5701 5702 .seealso: SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESSetFunction() 5703 */ 5704 PetscErrorCode SNESMonitorSetMatlab(SNES snes,const char *f,mxArray *ctx) 5705 { 5706 PetscErrorCode ierr; 5707 SNESMatlabContext *sctx; 5708 5709 PetscFunctionBegin; 5710 /* currently sctx is memory bleed */ 5711 ierr = PetscNew(&sctx);CHKERRQ(ierr); 5712 ierr = PetscStrallocpy(f,&sctx->funcname);CHKERRQ(ierr); 5713 /* 5714 This should work, but it doesn't 5715 sctx->ctx = ctx; 5716 mexMakeArrayPersistent(sctx->ctx); 5717 */ 5718 sctx->ctx = mxDuplicateArray(ctx); 5719 ierr = SNESMonitorSet(snes,SNESMonitor_Matlab,sctx,NULL);CHKERRQ(ierr); 5720 PetscFunctionReturn(0); 5721 } 5722 5723 #endif 5724