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