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