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