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