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