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