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) || snes->reason) { 3326 ierr = PetscDrawLGDraw(lg);CHKERRQ(ierr); 3327 ierr = PetscDrawLGSave(lg);CHKERRQ(ierr); 3328 } 3329 3330 ierr = PetscViewerDrawGetDrawLG(v,1,&lg);CHKERRQ(ierr); 3331 if (!n) {ierr = PetscDrawLGReset(lg);CHKERRQ(ierr);} 3332 ierr = PetscDrawLGGetDraw(lg,&draw);CHKERRQ(ierr); 3333 ierr = PetscDrawSetTitle(draw,"% elemts > .2*max elemt");CHKERRQ(ierr); 3334 ierr = SNESMonitorRange_Private(snes,n,&per);CHKERRQ(ierr); 3335 x = (PetscReal)n; 3336 y = 100.0*per; 3337 ierr = PetscDrawLGAddPoint(lg,&x,&y);CHKERRQ(ierr); 3338 if (n < 20 || !(n % 5) || snes->reason) { 3339 ierr = PetscDrawLGDraw(lg);CHKERRQ(ierr); 3340 ierr = PetscDrawLGSave(lg);CHKERRQ(ierr); 3341 } 3342 3343 ierr = PetscViewerDrawGetDrawLG(v,2,&lg);CHKERRQ(ierr); 3344 if (!n) {prev = rnorm;ierr = PetscDrawLGReset(lg);CHKERRQ(ierr);} 3345 ierr = PetscDrawLGGetDraw(lg,&draw);CHKERRQ(ierr); 3346 ierr = PetscDrawSetTitle(draw,"(norm -oldnorm)/oldnorm");CHKERRQ(ierr); 3347 x = (PetscReal)n; 3348 y = (prev - rnorm)/prev; 3349 ierr = PetscDrawLGAddPoint(lg,&x,&y);CHKERRQ(ierr); 3350 if (n < 20 || !(n % 5) || snes->reason) { 3351 ierr = PetscDrawLGDraw(lg);CHKERRQ(ierr); 3352 ierr = PetscDrawLGSave(lg);CHKERRQ(ierr); 3353 } 3354 3355 ierr = PetscViewerDrawGetDrawLG(v,3,&lg);CHKERRQ(ierr); 3356 if (!n) {ierr = PetscDrawLGReset(lg);CHKERRQ(ierr);} 3357 ierr = PetscDrawLGGetDraw(lg,&draw);CHKERRQ(ierr); 3358 ierr = PetscDrawSetTitle(draw,"(norm -oldnorm)/oldnorm*(% > .2 max)");CHKERRQ(ierr); 3359 x = (PetscReal)n; 3360 y = (prev - rnorm)/(prev*per); 3361 if (n > 2) { /*skip initial crazy value */ 3362 ierr = PetscDrawLGAddPoint(lg,&x,&y);CHKERRQ(ierr); 3363 } 3364 if (n < 20 || !(n % 5) || snes->reason) { 3365 ierr = PetscDrawLGDraw(lg);CHKERRQ(ierr); 3366 ierr = PetscDrawLGSave(lg);CHKERRQ(ierr); 3367 } 3368 prev = rnorm; 3369 PetscFunctionReturn(0); 3370 } 3371 3372 #undef __FUNCT__ 3373 #define __FUNCT__ "SNESMonitor" 3374 /*@ 3375 SNESMonitor - runs the user provided monitor routines, if they exist 3376 3377 Collective on SNES 3378 3379 Input Parameters: 3380 + snes - nonlinear solver context obtained from SNESCreate() 3381 . iter - iteration number 3382 - rnorm - relative norm of the residual 3383 3384 Notes: 3385 This routine is called by the SNES implementations. 3386 It does not typically need to be called by the user. 3387 3388 Level: developer 3389 3390 .seealso: SNESMonitorSet() 3391 @*/ 3392 PetscErrorCode SNESMonitor(SNES snes,PetscInt iter,PetscReal rnorm) 3393 { 3394 PetscErrorCode ierr; 3395 PetscInt i,n = snes->numbermonitors; 3396 3397 PetscFunctionBegin; 3398 ierr = VecLockPush(snes->vec_sol);CHKERRQ(ierr); 3399 for (i=0; i<n; i++) { 3400 ierr = (*snes->monitor[i])(snes,iter,rnorm,snes->monitorcontext[i]);CHKERRQ(ierr); 3401 } 3402 ierr = VecLockPop(snes->vec_sol);CHKERRQ(ierr); 3403 PetscFunctionReturn(0); 3404 } 3405 3406 /* ------------ Routines to set performance monitoring options ----------- */ 3407 3408 /*MC 3409 SNESMonitorFunction - functional form passed to SNESMonitorSet() to monitor convergence of nonlinear solver 3410 3411 Synopsis: 3412 #include <petscsnes.h> 3413 $ PetscErrorCode SNESMonitorFunction(SNES snes,PetscInt its, PetscReal norm,void *mctx) 3414 3415 + snes - the SNES context 3416 . its - iteration number 3417 . norm - 2-norm function value (may be estimated) 3418 - mctx - [optional] monitoring context 3419 3420 Level: advanced 3421 3422 .seealso: SNESMonitorSet(), SNESMonitorGet() 3423 M*/ 3424 3425 #undef __FUNCT__ 3426 #define __FUNCT__ "SNESMonitorSet" 3427 /*@C 3428 SNESMonitorSet - Sets an ADDITIONAL function that is to be used at every 3429 iteration of the nonlinear solver to display the iteration's 3430 progress. 3431 3432 Logically Collective on SNES 3433 3434 Input Parameters: 3435 + snes - the SNES context 3436 . f - the monitor function, see SNESMonitorFunction for the calling sequence 3437 . mctx - [optional] user-defined context for private data for the 3438 monitor routine (use NULL if no context is desired) 3439 - monitordestroy - [optional] routine that frees monitor context 3440 (may be NULL) 3441 3442 Options Database Keys: 3443 + -snes_monitor - sets SNESMonitorDefault() 3444 . -snes_monitor_lg_residualnorm - sets line graph monitor, 3445 uses SNESMonitorLGCreate() 3446 - -snes_monitor_cancel - cancels all monitors that have 3447 been hardwired into a code by 3448 calls to SNESMonitorSet(), but 3449 does not cancel those set via 3450 the options database. 3451 3452 Notes: 3453 Several different monitoring routines may be set by calling 3454 SNESMonitorSet() multiple times; all will be called in the 3455 order in which they were set. 3456 3457 Fortran notes: Only a single monitor function can be set for each SNES object 3458 3459 Level: intermediate 3460 3461 .keywords: SNES, nonlinear, set, monitor 3462 3463 .seealso: SNESMonitorDefault(), SNESMonitorCancel(), SNESMonitorFunction 3464 @*/ 3465 PetscErrorCode SNESMonitorSet(SNES snes,PetscErrorCode (*f)(SNES,PetscInt,PetscReal,void*),void *mctx,PetscErrorCode (*monitordestroy)(void**)) 3466 { 3467 PetscInt i; 3468 PetscErrorCode ierr; 3469 3470 PetscFunctionBegin; 3471 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3472 if (snes->numbermonitors >= MAXSNESMONITORS) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Too many monitors set"); 3473 for (i=0; i<snes->numbermonitors;i++) { 3474 if (f == snes->monitor[i] && monitordestroy == snes->monitordestroy[i] && mctx == snes->monitorcontext[i]) { 3475 if (monitordestroy) { 3476 ierr = (*monitordestroy)(&mctx);CHKERRQ(ierr); 3477 } 3478 PetscFunctionReturn(0); 3479 } 3480 } 3481 snes->monitor[snes->numbermonitors] = f; 3482 snes->monitordestroy[snes->numbermonitors] = monitordestroy; 3483 snes->monitorcontext[snes->numbermonitors++] = (void*)mctx; 3484 PetscFunctionReturn(0); 3485 } 3486 3487 #undef __FUNCT__ 3488 #define __FUNCT__ "SNESMonitorCancel" 3489 /*@ 3490 SNESMonitorCancel - Clears all the monitor functions for a SNES object. 3491 3492 Logically Collective on SNES 3493 3494 Input Parameters: 3495 . snes - the SNES context 3496 3497 Options Database Key: 3498 . -snes_monitor_cancel - cancels all monitors that have been hardwired 3499 into a code by calls to SNESMonitorSet(), but does not cancel those 3500 set via the options database 3501 3502 Notes: 3503 There is no way to clear one specific monitor from a SNES object. 3504 3505 Level: intermediate 3506 3507 .keywords: SNES, nonlinear, set, monitor 3508 3509 .seealso: SNESMonitorDefault(), SNESMonitorSet() 3510 @*/ 3511 PetscErrorCode SNESMonitorCancel(SNES snes) 3512 { 3513 PetscErrorCode ierr; 3514 PetscInt i; 3515 3516 PetscFunctionBegin; 3517 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3518 for (i=0; i<snes->numbermonitors; i++) { 3519 if (snes->monitordestroy[i]) { 3520 ierr = (*snes->monitordestroy[i])(&snes->monitorcontext[i]);CHKERRQ(ierr); 3521 } 3522 } 3523 snes->numbermonitors = 0; 3524 PetscFunctionReturn(0); 3525 } 3526 3527 /*MC 3528 SNESConvergenceTestFunction - functional form used for testing of convergence of nonlinear solver 3529 3530 Synopsis: 3531 #include <petscsnes.h> 3532 $ PetscErrorCode SNESConvergenceTest(SNES snes,PetscInt it,PetscReal xnorm,PetscReal gnorm,PetscReal f,SNESConvergedReason *reason,void *cctx) 3533 3534 + snes - the SNES context 3535 . it - current iteration (0 is the first and is before any Newton step) 3536 . cctx - [optional] convergence context 3537 . reason - reason for convergence/divergence 3538 . xnorm - 2-norm of current iterate 3539 . gnorm - 2-norm of current step 3540 - f - 2-norm of function 3541 3542 Level: intermediate 3543 3544 .seealso: SNESSetConvergenceTest(), SNESGetConvergenceTest() 3545 M*/ 3546 3547 #undef __FUNCT__ 3548 #define __FUNCT__ "SNESSetConvergenceTest" 3549 /*@C 3550 SNESSetConvergenceTest - Sets the function that is to be used 3551 to test for convergence of the nonlinear iterative solution. 3552 3553 Logically Collective on SNES 3554 3555 Input Parameters: 3556 + snes - the SNES context 3557 . SNESConvergenceTestFunction - routine to test for convergence 3558 . cctx - [optional] context for private data for the convergence routine (may be NULL) 3559 - destroy - [optional] destructor for the context (may be NULL; NULL_FUNCTION in Fortran) 3560 3561 Level: advanced 3562 3563 .keywords: SNES, nonlinear, set, convergence, test 3564 3565 .seealso: SNESConvergedDefault(), SNESConvergedSkip(), SNESConvergenceTestFunction 3566 @*/ 3567 PetscErrorCode SNESSetConvergenceTest(SNES snes,PetscErrorCode (*SNESConvergenceTestFunction)(SNES,PetscInt,PetscReal,PetscReal,PetscReal,SNESConvergedReason*,void*),void *cctx,PetscErrorCode (*destroy)(void*)) 3568 { 3569 PetscErrorCode ierr; 3570 3571 PetscFunctionBegin; 3572 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3573 if (!SNESConvergenceTestFunction) SNESConvergenceTestFunction = SNESConvergedSkip; 3574 if (snes->ops->convergeddestroy) { 3575 ierr = (*snes->ops->convergeddestroy)(snes->cnvP);CHKERRQ(ierr); 3576 } 3577 snes->ops->converged = SNESConvergenceTestFunction; 3578 snes->ops->convergeddestroy = destroy; 3579 snes->cnvP = cctx; 3580 PetscFunctionReturn(0); 3581 } 3582 3583 #undef __FUNCT__ 3584 #define __FUNCT__ "SNESGetConvergedReason" 3585 /*@ 3586 SNESGetConvergedReason - Gets the reason the SNES iteration was stopped. 3587 3588 Not Collective 3589 3590 Input Parameter: 3591 . snes - the SNES context 3592 3593 Output Parameter: 3594 . reason - negative value indicates diverged, positive value converged, see SNESConvergedReason or the 3595 manual pages for the individual convergence tests for complete lists 3596 3597 Level: intermediate 3598 3599 Notes: Can only be called after the call the SNESSolve() is complete. 3600 3601 .keywords: SNES, nonlinear, set, convergence, test 3602 3603 .seealso: SNESSetConvergenceTest(), SNESSetConvergedReason(), SNESConvergedReason 3604 @*/ 3605 PetscErrorCode SNESGetConvergedReason(SNES snes,SNESConvergedReason *reason) 3606 { 3607 PetscFunctionBegin; 3608 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3609 PetscValidPointer(reason,2); 3610 *reason = snes->reason; 3611 PetscFunctionReturn(0); 3612 } 3613 3614 #undef __FUNCT__ 3615 #define __FUNCT__ "SNESSetConvergedReason" 3616 /*@ 3617 SNESSetConvergedReason - Sets the reason the SNES iteration was stopped. 3618 3619 Not Collective 3620 3621 Input Parameters: 3622 + snes - the SNES context 3623 - reason - negative value indicates diverged, positive value converged, see SNESConvergedReason or the 3624 manual pages for the individual convergence tests for complete lists 3625 3626 Level: intermediate 3627 3628 .keywords: SNES, nonlinear, set, convergence, test 3629 .seealso: SNESGetConvergedReason(), SNESSetConvergenceTest(), SNESConvergedReason 3630 @*/ 3631 PetscErrorCode SNESSetConvergedReason(SNES snes,SNESConvergedReason reason) 3632 { 3633 PetscFunctionBegin; 3634 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3635 snes->reason = reason; 3636 PetscFunctionReturn(0); 3637 } 3638 3639 #undef __FUNCT__ 3640 #define __FUNCT__ "SNESSetConvergenceHistory" 3641 /*@ 3642 SNESSetConvergenceHistory - Sets the array used to hold the convergence history. 3643 3644 Logically Collective on SNES 3645 3646 Input Parameters: 3647 + snes - iterative context obtained from SNESCreate() 3648 . a - array to hold history, this array will contain the function norms computed at each step 3649 . its - integer array holds the number of linear iterations for each solve. 3650 . na - size of a and its 3651 - reset - PETSC_TRUE indicates each new nonlinear solve resets the history counter to zero, 3652 else it continues storing new values for new nonlinear solves after the old ones 3653 3654 Notes: 3655 If 'a' and 'its' are NULL then space is allocated for the history. If 'na' PETSC_DECIDE or PETSC_DEFAULT then a 3656 default array of length 10000 is allocated. 3657 3658 This routine is useful, e.g., when running a code for purposes 3659 of accurate performance monitoring, when no I/O should be done 3660 during the section of code that is being timed. 3661 3662 Level: intermediate 3663 3664 .keywords: SNES, set, convergence, history 3665 3666 .seealso: SNESGetConvergenceHistory() 3667 3668 @*/ 3669 PetscErrorCode SNESSetConvergenceHistory(SNES snes,PetscReal a[],PetscInt its[],PetscInt na,PetscBool reset) 3670 { 3671 PetscErrorCode ierr; 3672 3673 PetscFunctionBegin; 3674 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3675 if (a) PetscValidScalarPointer(a,2); 3676 if (its) PetscValidIntPointer(its,3); 3677 if (!a) { 3678 if (na == PETSC_DECIDE || na == PETSC_DEFAULT) na = 1000; 3679 ierr = PetscCalloc1(na,&a);CHKERRQ(ierr); 3680 ierr = PetscCalloc1(na,&its);CHKERRQ(ierr); 3681 3682 snes->conv_malloc = PETSC_TRUE; 3683 } 3684 snes->conv_hist = a; 3685 snes->conv_hist_its = its; 3686 snes->conv_hist_max = na; 3687 snes->conv_hist_len = 0; 3688 snes->conv_hist_reset = reset; 3689 PetscFunctionReturn(0); 3690 } 3691 3692 #if defined(PETSC_HAVE_MATLAB_ENGINE) 3693 #include <engine.h> /* MATLAB include file */ 3694 #include <mex.h> /* MATLAB include file */ 3695 3696 #undef __FUNCT__ 3697 #define __FUNCT__ "SNESGetConvergenceHistoryMatlab" 3698 PETSC_EXTERN mxArray *SNESGetConvergenceHistoryMatlab(SNES snes) 3699 { 3700 mxArray *mat; 3701 PetscInt i; 3702 PetscReal *ar; 3703 3704 PetscFunctionBegin; 3705 mat = mxCreateDoubleMatrix(snes->conv_hist_len,1,mxREAL); 3706 ar = (PetscReal*) mxGetData(mat); 3707 for (i=0; i<snes->conv_hist_len; i++) ar[i] = snes->conv_hist[i]; 3708 PetscFunctionReturn(mat); 3709 } 3710 #endif 3711 3712 #undef __FUNCT__ 3713 #define __FUNCT__ "SNESGetConvergenceHistory" 3714 /*@C 3715 SNESGetConvergenceHistory - Gets the array used to hold the convergence history. 3716 3717 Not Collective 3718 3719 Input Parameter: 3720 . snes - iterative context obtained from SNESCreate() 3721 3722 Output Parameters: 3723 . a - array to hold history 3724 . its - integer array holds the number of linear iterations (or 3725 negative if not converged) for each solve. 3726 - na - size of a and its 3727 3728 Notes: 3729 The calling sequence for this routine in Fortran is 3730 $ call SNESGetConvergenceHistory(SNES snes, integer na, integer ierr) 3731 3732 This routine is useful, e.g., when running a code for purposes 3733 of accurate performance monitoring, when no I/O should be done 3734 during the section of code that is being timed. 3735 3736 Level: intermediate 3737 3738 .keywords: SNES, get, convergence, history 3739 3740 .seealso: SNESSetConvergencHistory() 3741 3742 @*/ 3743 PetscErrorCode SNESGetConvergenceHistory(SNES snes,PetscReal *a[],PetscInt *its[],PetscInt *na) 3744 { 3745 PetscFunctionBegin; 3746 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3747 if (a) *a = snes->conv_hist; 3748 if (its) *its = snes->conv_hist_its; 3749 if (na) *na = snes->conv_hist_len; 3750 PetscFunctionReturn(0); 3751 } 3752 3753 #undef __FUNCT__ 3754 #define __FUNCT__ "SNESSetUpdate" 3755 /*@C 3756 SNESSetUpdate - Sets the general-purpose update function called 3757 at the beginning of every iteration of the nonlinear solve. Specifically 3758 it is called just before the Jacobian is "evaluated". 3759 3760 Logically Collective on SNES 3761 3762 Input Parameters: 3763 . snes - The nonlinear solver context 3764 . func - The function 3765 3766 Calling sequence of func: 3767 . func (SNES snes, PetscInt step); 3768 3769 . step - The current step of the iteration 3770 3771 Level: advanced 3772 3773 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() 3774 This is not used by most users. 3775 3776 .keywords: SNES, update 3777 3778 .seealso SNESSetJacobian(), SNESSolve() 3779 @*/ 3780 PetscErrorCode SNESSetUpdate(SNES snes, PetscErrorCode (*func)(SNES, PetscInt)) 3781 { 3782 PetscFunctionBegin; 3783 PetscValidHeaderSpecific(snes, SNES_CLASSID,1); 3784 snes->ops->update = func; 3785 PetscFunctionReturn(0); 3786 } 3787 3788 #undef __FUNCT__ 3789 #define __FUNCT__ "SNESScaleStep_Private" 3790 /* 3791 SNESScaleStep_Private - Scales a step so that its length is less than the 3792 positive parameter delta. 3793 3794 Input Parameters: 3795 + snes - the SNES context 3796 . y - approximate solution of linear system 3797 . fnorm - 2-norm of current function 3798 - delta - trust region size 3799 3800 Output Parameters: 3801 + gpnorm - predicted function norm at the new point, assuming local 3802 linearization. The value is zero if the step lies within the trust 3803 region, and exceeds zero otherwise. 3804 - ynorm - 2-norm of the step 3805 3806 Note: 3807 For non-trust region methods such as SNESNEWTONLS, the parameter delta 3808 is set to be the maximum allowable step size. 3809 3810 .keywords: SNES, nonlinear, scale, step 3811 */ 3812 PetscErrorCode SNESScaleStep_Private(SNES snes,Vec y,PetscReal *fnorm,PetscReal *delta,PetscReal *gpnorm,PetscReal *ynorm) 3813 { 3814 PetscReal nrm; 3815 PetscScalar cnorm; 3816 PetscErrorCode ierr; 3817 3818 PetscFunctionBegin; 3819 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3820 PetscValidHeaderSpecific(y,VEC_CLASSID,2); 3821 PetscCheckSameComm(snes,1,y,2); 3822 3823 ierr = VecNorm(y,NORM_2,&nrm);CHKERRQ(ierr); 3824 if (nrm > *delta) { 3825 nrm = *delta/nrm; 3826 *gpnorm = (1.0 - nrm)*(*fnorm); 3827 cnorm = nrm; 3828 ierr = VecScale(y,cnorm);CHKERRQ(ierr); 3829 *ynorm = *delta; 3830 } else { 3831 *gpnorm = 0.0; 3832 *ynorm = nrm; 3833 } 3834 PetscFunctionReturn(0); 3835 } 3836 3837 #undef __FUNCT__ 3838 #define __FUNCT__ "SNESReasonView" 3839 /*@ 3840 SNESReasonView - Displays the reason a SNES solve converged or diverged to a viewer 3841 3842 Collective on SNES 3843 3844 Parameter: 3845 + snes - iterative context obtained from SNESCreate() 3846 - viewer - the viewer to display the reason 3847 3848 3849 Options Database Keys: 3850 . -snes_converged_reason - print reason for converged or diverged, also prints number of iterations 3851 3852 Level: beginner 3853 3854 .keywords: SNES, solve, linear system 3855 3856 .seealso: SNESCreate(), SNESSetUp(), SNESDestroy(), SNESSetTolerances(), SNESConvergedDefault() 3857 3858 @*/ 3859 PetscErrorCode SNESReasonView(SNES snes,PetscViewer viewer) 3860 { 3861 PetscErrorCode ierr; 3862 PetscBool isAscii; 3863 3864 PetscFunctionBegin; 3865 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isAscii);CHKERRQ(ierr); 3866 if (isAscii) { 3867 ierr = PetscViewerASCIIAddTab(viewer,((PetscObject)snes)->tablevel);CHKERRQ(ierr); 3868 if (snes->reason > 0) { 3869 if (((PetscObject) snes)->prefix) { 3870 ierr = PetscViewerASCIIPrintf(viewer,"Nonlinear %s solve converged due to %s iterations %D\n",((PetscObject) snes)->prefix,SNESConvergedReasons[snes->reason],snes->iter);CHKERRQ(ierr); 3871 } else { 3872 ierr = PetscViewerASCIIPrintf(viewer,"Nonlinear solve converged due to %s iterations %D\n",SNESConvergedReasons[snes->reason],snes->iter);CHKERRQ(ierr); 3873 } 3874 } else { 3875 if (((PetscObject) snes)->prefix) { 3876 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); 3877 } else { 3878 ierr = PetscViewerASCIIPrintf(viewer,"Nonlinear solve did not converge due to %s iterations %D\n",SNESConvergedReasons[snes->reason],snes->iter);CHKERRQ(ierr); 3879 } 3880 } 3881 ierr = PetscViewerASCIISubtractTab(viewer,((PetscObject)snes)->tablevel);CHKERRQ(ierr); 3882 } 3883 PetscFunctionReturn(0); 3884 } 3885 3886 #undef __FUNCT__ 3887 #define __FUNCT__ "SNESReasonViewFromOptions" 3888 /*@C 3889 SNESReasonViewFromOptions - Processes command line options to determine if/how a SNESReason is to be viewed. 3890 3891 Collective on SNES 3892 3893 Input Parameters: 3894 . snes - the SNES object 3895 3896 Level: intermediate 3897 3898 @*/ 3899 PetscErrorCode SNESReasonViewFromOptions(SNES snes) 3900 { 3901 PetscErrorCode ierr; 3902 PetscViewer viewer; 3903 PetscBool flg; 3904 static PetscBool incall = PETSC_FALSE; 3905 PetscViewerFormat format; 3906 3907 PetscFunctionBegin; 3908 if (incall) PetscFunctionReturn(0); 3909 incall = PETSC_TRUE; 3910 ierr = PetscOptionsGetViewer(PetscObjectComm((PetscObject)snes),((PetscObject)snes)->prefix,"-snes_converged_reason",&viewer,&format,&flg);CHKERRQ(ierr); 3911 if (flg) { 3912 ierr = PetscViewerPushFormat(viewer,format);CHKERRQ(ierr); 3913 ierr = SNESReasonView(snes,viewer);CHKERRQ(ierr); 3914 ierr = PetscViewerPopFormat(viewer);CHKERRQ(ierr); 3915 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 3916 } 3917 incall = PETSC_FALSE; 3918 PetscFunctionReturn(0); 3919 } 3920 3921 #undef __FUNCT__ 3922 #define __FUNCT__ "SNESSolve" 3923 /*@C 3924 SNESSolve - Solves a nonlinear system F(x) = b. 3925 Call SNESSolve() after calling SNESCreate() and optional routines of the form SNESSetXXX(). 3926 3927 Collective on SNES 3928 3929 Input Parameters: 3930 + snes - the SNES context 3931 . b - the constant part of the equation F(x) = b, or NULL to use zero. 3932 - x - the solution vector. 3933 3934 Notes: 3935 The user should initialize the vector,x, with the initial guess 3936 for the nonlinear solve prior to calling SNESSolve. In particular, 3937 to employ an initial guess of zero, the user should explicitly set 3938 this vector to zero by calling VecSet(). 3939 3940 Level: beginner 3941 3942 .keywords: SNES, nonlinear, solve 3943 3944 .seealso: SNESCreate(), SNESDestroy(), SNESSetFunction(), SNESSetJacobian(), SNESSetGridSequence(), SNESGetSolution() 3945 @*/ 3946 PetscErrorCode SNESSolve(SNES snes,Vec b,Vec x) 3947 { 3948 PetscErrorCode ierr; 3949 PetscBool flg; 3950 PetscInt grid; 3951 Vec xcreated = NULL; 3952 DM dm; 3953 3954 PetscFunctionBegin; 3955 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3956 if (x) PetscValidHeaderSpecific(x,VEC_CLASSID,3); 3957 if (x) PetscCheckSameComm(snes,1,x,3); 3958 if (b) PetscValidHeaderSpecific(b,VEC_CLASSID,2); 3959 if (b) PetscCheckSameComm(snes,1,b,2); 3960 3961 if (!x) { 3962 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 3963 ierr = DMCreateGlobalVector(dm,&xcreated);CHKERRQ(ierr); 3964 x = xcreated; 3965 } 3966 ierr = SNESViewFromOptions(snes,NULL,"-snes_view_pre");CHKERRQ(ierr); 3967 3968 for (grid=0; grid<snes->gridsequence; grid++) {ierr = PetscViewerASCIIPushTab(PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)snes)));CHKERRQ(ierr);} 3969 for (grid=0; grid<snes->gridsequence+1; grid++) { 3970 3971 /* set solution vector */ 3972 if (!grid) {ierr = PetscObjectReference((PetscObject)x);CHKERRQ(ierr);} 3973 ierr = VecDestroy(&snes->vec_sol);CHKERRQ(ierr); 3974 snes->vec_sol = x; 3975 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 3976 3977 /* set affine vector if provided */ 3978 if (b) { ierr = PetscObjectReference((PetscObject)b);CHKERRQ(ierr); } 3979 ierr = VecDestroy(&snes->vec_rhs);CHKERRQ(ierr); 3980 snes->vec_rhs = b; 3981 3982 if (snes->vec_func == snes->vec_sol) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_IDN,"Solution vector cannot be function vector"); 3983 if (snes->vec_rhs == snes->vec_sol) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_IDN,"Solution vector cannot be right hand side vector"); 3984 if (!snes->vec_sol_update /* && snes->vec_sol */) { 3985 ierr = VecDuplicate(snes->vec_sol,&snes->vec_sol_update);CHKERRQ(ierr); 3986 ierr = PetscLogObjectParent((PetscObject)snes,(PetscObject)snes->vec_sol_update);CHKERRQ(ierr); 3987 } 3988 ierr = DMShellSetGlobalVector(dm,snes->vec_sol);CHKERRQ(ierr); 3989 ierr = SNESSetUp(snes);CHKERRQ(ierr); 3990 3991 if (!grid) { 3992 if (snes->ops->computeinitialguess) { 3993 ierr = (*snes->ops->computeinitialguess)(snes,snes->vec_sol,snes->initialguessP);CHKERRQ(ierr); 3994 } 3995 } 3996 3997 if (snes->conv_hist_reset) snes->conv_hist_len = 0; 3998 if (snes->counters_reset) {snes->nfuncs = 0; snes->linear_its = 0; snes->numFailures = 0;} 3999 4000 ierr = PetscLogEventBegin(SNES_Solve,snes,0,0,0);CHKERRQ(ierr); 4001 ierr = (*snes->ops->solve)(snes);CHKERRQ(ierr); 4002 ierr = PetscLogEventEnd(SNES_Solve,snes,0,0,0);CHKERRQ(ierr); 4003 if (!snes->reason) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Internal error, solver returned without setting converged reason"); 4004 snes->domainerror = PETSC_FALSE; /* clear the flag if it has been set */ 4005 4006 if (snes->lagjac_persist) snes->jac_iter += snes->iter; 4007 if (snes->lagpre_persist) snes->pre_iter += snes->iter; 4008 4009 ierr = PetscOptionsGetViewer(PetscObjectComm((PetscObject)snes),((PetscObject)snes)->prefix,"-snes_test_local_min",NULL,NULL,&flg);CHKERRQ(ierr); 4010 if (flg && !PetscPreLoadingOn) { ierr = SNESTestLocalMin(snes);CHKERRQ(ierr); } 4011 ierr = SNESReasonViewFromOptions(snes);CHKERRQ(ierr); 4012 4013 if (snes->errorifnotconverged && snes->reason < 0) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_NOT_CONVERGED,"SNESSolve has not converged"); 4014 if (snes->reason < 0) break; 4015 if (grid < snes->gridsequence) { 4016 DM fine; 4017 Vec xnew; 4018 Mat interp; 4019 4020 ierr = DMRefine(snes->dm,PetscObjectComm((PetscObject)snes),&fine);CHKERRQ(ierr); 4021 if (!fine) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_INCOMP,"DMRefine() did not perform any refinement, cannot continue grid sequencing"); 4022 ierr = DMCreateInterpolation(snes->dm,fine,&interp,NULL);CHKERRQ(ierr); 4023 ierr = DMCreateGlobalVector(fine,&xnew);CHKERRQ(ierr); 4024 ierr = MatInterpolate(interp,x,xnew);CHKERRQ(ierr); 4025 ierr = DMInterpolate(snes->dm,interp,fine);CHKERRQ(ierr); 4026 ierr = MatDestroy(&interp);CHKERRQ(ierr); 4027 x = xnew; 4028 4029 ierr = SNESReset(snes);CHKERRQ(ierr); 4030 ierr = SNESSetDM(snes,fine);CHKERRQ(ierr); 4031 ierr = DMDestroy(&fine);CHKERRQ(ierr); 4032 ierr = PetscViewerASCIIPopTab(PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)snes)));CHKERRQ(ierr); 4033 } 4034 } 4035 ierr = SNESViewFromOptions(snes,NULL,"-snes_view");CHKERRQ(ierr); 4036 ierr = VecViewFromOptions(snes->vec_sol,(PetscObject)snes,"-snes_view_solution");CHKERRQ(ierr); 4037 4038 ierr = VecDestroy(&xcreated);CHKERRQ(ierr); 4039 ierr = PetscObjectSAWsBlock((PetscObject)snes);CHKERRQ(ierr); 4040 PetscFunctionReturn(0); 4041 } 4042 4043 /* --------- Internal routines for SNES Package --------- */ 4044 4045 #undef __FUNCT__ 4046 #define __FUNCT__ "SNESSetType" 4047 /*@C 4048 SNESSetType - Sets the method for the nonlinear solver. 4049 4050 Collective on SNES 4051 4052 Input Parameters: 4053 + snes - the SNES context 4054 - type - a known method 4055 4056 Options Database Key: 4057 . -snes_type <type> - Sets the method; use -help for a list 4058 of available methods (for instance, newtonls or newtontr) 4059 4060 Notes: 4061 See "petsc/include/petscsnes.h" for available methods (for instance) 4062 + SNESNEWTONLS - Newton's method with line search 4063 (systems of nonlinear equations) 4064 . SNESNEWTONTR - Newton's method with trust region 4065 (systems of nonlinear equations) 4066 4067 Normally, it is best to use the SNESSetFromOptions() command and then 4068 set the SNES solver type from the options database rather than by using 4069 this routine. Using the options database provides the user with 4070 maximum flexibility in evaluating the many nonlinear solvers. 4071 The SNESSetType() routine is provided for those situations where it 4072 is necessary to set the nonlinear solver independently of the command 4073 line or options database. This might be the case, for example, when 4074 the choice of solver changes during the execution of the program, 4075 and the user's application is taking responsibility for choosing the 4076 appropriate method. 4077 4078 Developer Notes: SNESRegister() adds a constructor for a new SNESType to SNESList, SNESSetType() locates 4079 the constructor in that list and calls it to create the spexific object. 4080 4081 Level: intermediate 4082 4083 .keywords: SNES, set, type 4084 4085 .seealso: SNESType, SNESCreate(), SNESDestroy(), SNESGetType(), SNESSetFromOptions() 4086 4087 @*/ 4088 PetscErrorCode SNESSetType(SNES snes,SNESType type) 4089 { 4090 PetscErrorCode ierr,(*r)(SNES); 4091 PetscBool match; 4092 4093 PetscFunctionBegin; 4094 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4095 PetscValidCharPointer(type,2); 4096 4097 ierr = PetscObjectTypeCompare((PetscObject)snes,type,&match);CHKERRQ(ierr); 4098 if (match) PetscFunctionReturn(0); 4099 4100 ierr = PetscFunctionListFind(SNESList,type,&r);CHKERRQ(ierr); 4101 if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unable to find requested SNES type %s",type); 4102 /* Destroy the previous private SNES context */ 4103 if (snes->ops->destroy) { 4104 ierr = (*(snes)->ops->destroy)(snes);CHKERRQ(ierr); 4105 snes->ops->destroy = NULL; 4106 } 4107 /* Reinitialize function pointers in SNESOps structure */ 4108 snes->ops->setup = 0; 4109 snes->ops->solve = 0; 4110 snes->ops->view = 0; 4111 snes->ops->setfromoptions = 0; 4112 snes->ops->destroy = 0; 4113 /* Call the SNESCreate_XXX routine for this particular Nonlinear solver */ 4114 snes->setupcalled = PETSC_FALSE; 4115 4116 ierr = PetscObjectChangeTypeName((PetscObject)snes,type);CHKERRQ(ierr); 4117 ierr = (*r)(snes);CHKERRQ(ierr); 4118 PetscFunctionReturn(0); 4119 } 4120 4121 #undef __FUNCT__ 4122 #define __FUNCT__ "SNESGetType" 4123 /*@C 4124 SNESGetType - Gets the SNES method type and name (as a string). 4125 4126 Not Collective 4127 4128 Input Parameter: 4129 . snes - nonlinear solver context 4130 4131 Output Parameter: 4132 . type - SNES method (a character string) 4133 4134 Level: intermediate 4135 4136 .keywords: SNES, nonlinear, get, type, name 4137 @*/ 4138 PetscErrorCode SNESGetType(SNES snes,SNESType *type) 4139 { 4140 PetscFunctionBegin; 4141 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4142 PetscValidPointer(type,2); 4143 *type = ((PetscObject)snes)->type_name; 4144 PetscFunctionReturn(0); 4145 } 4146 4147 #undef __FUNCT__ 4148 #define __FUNCT__ "SNESSetSolution" 4149 /*@ 4150 SNESSetSolution - Sets the solution vector for use by the SNES routines. 4151 4152 Logically Collective on SNES and Vec 4153 4154 Input Parameters: 4155 + snes - the SNES context obtained from SNESCreate() 4156 - u - the solution vector 4157 4158 Level: beginner 4159 4160 .keywords: SNES, set, solution 4161 @*/ 4162 PetscErrorCode SNESSetSolution(SNES snes, Vec u) 4163 { 4164 DM dm; 4165 PetscErrorCode ierr; 4166 4167 PetscFunctionBegin; 4168 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); 4169 PetscValidHeaderSpecific(u, VEC_CLASSID, 2); 4170 ierr = PetscObjectReference((PetscObject) u);CHKERRQ(ierr); 4171 ierr = VecDestroy(&snes->vec_sol);CHKERRQ(ierr); 4172 4173 snes->vec_sol = u; 4174 4175 ierr = SNESGetDM(snes, &dm);CHKERRQ(ierr); 4176 ierr = DMShellSetGlobalVector(dm, u);CHKERRQ(ierr); 4177 PetscFunctionReturn(0); 4178 } 4179 4180 #undef __FUNCT__ 4181 #define __FUNCT__ "SNESGetSolution" 4182 /*@ 4183 SNESGetSolution - Returns the vector where the approximate solution is 4184 stored. This is the fine grid solution when using SNESSetGridSequence(). 4185 4186 Not Collective, but Vec is parallel if SNES is parallel 4187 4188 Input Parameter: 4189 . snes - the SNES context 4190 4191 Output Parameter: 4192 . x - the solution 4193 4194 Level: intermediate 4195 4196 .keywords: SNES, nonlinear, get, solution 4197 4198 .seealso: SNESGetSolutionUpdate(), SNESGetFunction() 4199 @*/ 4200 PetscErrorCode SNESGetSolution(SNES snes,Vec *x) 4201 { 4202 PetscFunctionBegin; 4203 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4204 PetscValidPointer(x,2); 4205 *x = snes->vec_sol; 4206 PetscFunctionReturn(0); 4207 } 4208 4209 #undef __FUNCT__ 4210 #define __FUNCT__ "SNESGetSolutionUpdate" 4211 /*@ 4212 SNESGetSolutionUpdate - Returns the vector where the solution update is 4213 stored. 4214 4215 Not Collective, but Vec is parallel if SNES is parallel 4216 4217 Input Parameter: 4218 . snes - the SNES context 4219 4220 Output Parameter: 4221 . x - the solution update 4222 4223 Level: advanced 4224 4225 .keywords: SNES, nonlinear, get, solution, update 4226 4227 .seealso: SNESGetSolution(), SNESGetFunction() 4228 @*/ 4229 PetscErrorCode SNESGetSolutionUpdate(SNES snes,Vec *x) 4230 { 4231 PetscFunctionBegin; 4232 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4233 PetscValidPointer(x,2); 4234 *x = snes->vec_sol_update; 4235 PetscFunctionReturn(0); 4236 } 4237 4238 #undef __FUNCT__ 4239 #define __FUNCT__ "SNESGetFunction" 4240 /*@C 4241 SNESGetFunction - Returns the vector where the function is stored. 4242 4243 Not Collective, but Vec is parallel if SNES is parallel. Collective if Vec is requested, but has not been created yet. 4244 4245 Input Parameter: 4246 . snes - the SNES context 4247 4248 Output Parameter: 4249 + r - the vector that is used to store residuals (or NULL if you don't want it) 4250 . f - the function (or NULL if you don't want it); see SNESFunction for calling sequence details 4251 - ctx - the function context (or NULL if you don't want it) 4252 4253 Level: advanced 4254 4255 .keywords: SNES, nonlinear, get, function 4256 4257 .seealso: SNESSetFunction(), SNESGetSolution(), SNESFunction 4258 @*/ 4259 PetscErrorCode SNESGetFunction(SNES snes,Vec *r,PetscErrorCode (**f)(SNES,Vec,Vec,void*),void **ctx) 4260 { 4261 PetscErrorCode ierr; 4262 DM dm; 4263 4264 PetscFunctionBegin; 4265 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4266 if (r) { 4267 if (!snes->vec_func) { 4268 if (snes->vec_rhs) { 4269 ierr = VecDuplicate(snes->vec_rhs,&snes->vec_func);CHKERRQ(ierr); 4270 } else if (snes->vec_sol) { 4271 ierr = VecDuplicate(snes->vec_sol,&snes->vec_func);CHKERRQ(ierr); 4272 } else if (snes->dm) { 4273 ierr = DMCreateGlobalVector(snes->dm,&snes->vec_func);CHKERRQ(ierr); 4274 } 4275 } 4276 *r = snes->vec_func; 4277 } 4278 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 4279 ierr = DMSNESGetFunction(dm,f,ctx);CHKERRQ(ierr); 4280 PetscFunctionReturn(0); 4281 } 4282 4283 /*@C 4284 SNESGetNGS - Returns the NGS function and context. 4285 4286 Input Parameter: 4287 . snes - the SNES context 4288 4289 Output Parameter: 4290 + f - the function (or NULL) see SNESNGSFunction for details 4291 - ctx - the function context (or NULL) 4292 4293 Level: advanced 4294 4295 .keywords: SNES, nonlinear, get, function 4296 4297 .seealso: SNESSetNGS(), SNESGetFunction() 4298 @*/ 4299 4300 #undef __FUNCT__ 4301 #define __FUNCT__ "SNESGetNGS" 4302 PetscErrorCode SNESGetNGS (SNES snes, PetscErrorCode (**f)(SNES, Vec, Vec, void*), void ** ctx) 4303 { 4304 PetscErrorCode ierr; 4305 DM dm; 4306 4307 PetscFunctionBegin; 4308 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4309 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 4310 ierr = DMSNESGetNGS(dm,f,ctx);CHKERRQ(ierr); 4311 PetscFunctionReturn(0); 4312 } 4313 4314 #undef __FUNCT__ 4315 #define __FUNCT__ "SNESSetOptionsPrefix" 4316 /*@C 4317 SNESSetOptionsPrefix - Sets the prefix used for searching for all 4318 SNES options in the database. 4319 4320 Logically Collective on SNES 4321 4322 Input Parameter: 4323 + snes - the SNES context 4324 - prefix - the prefix to prepend to all option names 4325 4326 Notes: 4327 A hyphen (-) must NOT be given at the beginning of the prefix name. 4328 The first character of all runtime options is AUTOMATICALLY the hyphen. 4329 4330 Level: advanced 4331 4332 .keywords: SNES, set, options, prefix, database 4333 4334 .seealso: SNESSetFromOptions() 4335 @*/ 4336 PetscErrorCode SNESSetOptionsPrefix(SNES snes,const char prefix[]) 4337 { 4338 PetscErrorCode ierr; 4339 4340 PetscFunctionBegin; 4341 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4342 ierr = PetscObjectSetOptionsPrefix((PetscObject)snes,prefix);CHKERRQ(ierr); 4343 if (!snes->ksp) {ierr = SNESGetKSP(snes,&snes->ksp);CHKERRQ(ierr);} 4344 if (snes->linesearch) { 4345 ierr = SNESGetLineSearch(snes,&snes->linesearch);CHKERRQ(ierr); 4346 ierr = PetscObjectSetOptionsPrefix((PetscObject)snes->linesearch,prefix);CHKERRQ(ierr); 4347 } 4348 ierr = KSPSetOptionsPrefix(snes->ksp,prefix);CHKERRQ(ierr); 4349 PetscFunctionReturn(0); 4350 } 4351 4352 #undef __FUNCT__ 4353 #define __FUNCT__ "SNESAppendOptionsPrefix" 4354 /*@C 4355 SNESAppendOptionsPrefix - Appends to the prefix used for searching for all 4356 SNES options in the database. 4357 4358 Logically Collective on SNES 4359 4360 Input Parameters: 4361 + snes - the SNES context 4362 - prefix - the prefix to prepend to all option names 4363 4364 Notes: 4365 A hyphen (-) must NOT be given at the beginning of the prefix name. 4366 The first character of all runtime options is AUTOMATICALLY the hyphen. 4367 4368 Level: advanced 4369 4370 .keywords: SNES, append, options, prefix, database 4371 4372 .seealso: SNESGetOptionsPrefix() 4373 @*/ 4374 PetscErrorCode SNESAppendOptionsPrefix(SNES snes,const char prefix[]) 4375 { 4376 PetscErrorCode ierr; 4377 4378 PetscFunctionBegin; 4379 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4380 ierr = PetscObjectAppendOptionsPrefix((PetscObject)snes,prefix);CHKERRQ(ierr); 4381 if (!snes->ksp) {ierr = SNESGetKSP(snes,&snes->ksp);CHKERRQ(ierr);} 4382 if (snes->linesearch) { 4383 ierr = SNESGetLineSearch(snes,&snes->linesearch);CHKERRQ(ierr); 4384 ierr = PetscObjectAppendOptionsPrefix((PetscObject)snes->linesearch,prefix);CHKERRQ(ierr); 4385 } 4386 ierr = KSPAppendOptionsPrefix(snes->ksp,prefix);CHKERRQ(ierr); 4387 PetscFunctionReturn(0); 4388 } 4389 4390 #undef __FUNCT__ 4391 #define __FUNCT__ "SNESGetOptionsPrefix" 4392 /*@C 4393 SNESGetOptionsPrefix - Sets the prefix used for searching for all 4394 SNES options in the database. 4395 4396 Not Collective 4397 4398 Input Parameter: 4399 . snes - the SNES context 4400 4401 Output Parameter: 4402 . prefix - pointer to the prefix string used 4403 4404 Notes: On the fortran side, the user should pass in a string 'prefix' of 4405 sufficient length to hold the prefix. 4406 4407 Level: advanced 4408 4409 .keywords: SNES, get, options, prefix, database 4410 4411 .seealso: SNESAppendOptionsPrefix() 4412 @*/ 4413 PetscErrorCode SNESGetOptionsPrefix(SNES snes,const char *prefix[]) 4414 { 4415 PetscErrorCode ierr; 4416 4417 PetscFunctionBegin; 4418 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4419 ierr = PetscObjectGetOptionsPrefix((PetscObject)snes,prefix);CHKERRQ(ierr); 4420 PetscFunctionReturn(0); 4421 } 4422 4423 4424 #undef __FUNCT__ 4425 #define __FUNCT__ "SNESRegister" 4426 /*@C 4427 SNESRegister - Adds a method to the nonlinear solver package. 4428 4429 Not collective 4430 4431 Input Parameters: 4432 + name_solver - name of a new user-defined solver 4433 - routine_create - routine to create method context 4434 4435 Notes: 4436 SNESRegister() may be called multiple times to add several user-defined solvers. 4437 4438 Sample usage: 4439 .vb 4440 SNESRegister("my_solver",MySolverCreate); 4441 .ve 4442 4443 Then, your solver can be chosen with the procedural interface via 4444 $ SNESSetType(snes,"my_solver") 4445 or at runtime via the option 4446 $ -snes_type my_solver 4447 4448 Level: advanced 4449 4450 Note: If your function is not being put into a shared library then use SNESRegister() instead 4451 4452 .keywords: SNES, nonlinear, register 4453 4454 .seealso: SNESRegisterAll(), SNESRegisterDestroy() 4455 4456 Level: advanced 4457 @*/ 4458 PetscErrorCode SNESRegister(const char sname[],PetscErrorCode (*function)(SNES)) 4459 { 4460 PetscErrorCode ierr; 4461 4462 PetscFunctionBegin; 4463 ierr = PetscFunctionListAdd(&SNESList,sname,function);CHKERRQ(ierr); 4464 PetscFunctionReturn(0); 4465 } 4466 4467 #undef __FUNCT__ 4468 #define __FUNCT__ "SNESTestLocalMin" 4469 PetscErrorCode SNESTestLocalMin(SNES snes) 4470 { 4471 PetscErrorCode ierr; 4472 PetscInt N,i,j; 4473 Vec u,uh,fh; 4474 PetscScalar value; 4475 PetscReal norm; 4476 4477 PetscFunctionBegin; 4478 ierr = SNESGetSolution(snes,&u);CHKERRQ(ierr); 4479 ierr = VecDuplicate(u,&uh);CHKERRQ(ierr); 4480 ierr = VecDuplicate(u,&fh);CHKERRQ(ierr); 4481 4482 /* currently only works for sequential */ 4483 ierr = PetscPrintf(PETSC_COMM_WORLD,"Testing FormFunction() for local min\n");CHKERRQ(ierr); 4484 ierr = VecGetSize(u,&N);CHKERRQ(ierr); 4485 for (i=0; i<N; i++) { 4486 ierr = VecCopy(u,uh);CHKERRQ(ierr); 4487 ierr = PetscPrintf(PETSC_COMM_WORLD,"i = %D\n",i);CHKERRQ(ierr); 4488 for (j=-10; j<11; j++) { 4489 value = PetscSign(j)*PetscExpReal(PetscAbs(j)-10.0); 4490 ierr = VecSetValue(uh,i,value,ADD_VALUES);CHKERRQ(ierr); 4491 ierr = SNESComputeFunction(snes,uh,fh);CHKERRQ(ierr); 4492 ierr = VecNorm(fh,NORM_2,&norm);CHKERRQ(ierr); 4493 ierr = PetscPrintf(PETSC_COMM_WORLD," j norm %D %18.16e\n",j,norm);CHKERRQ(ierr); 4494 value = -value; 4495 ierr = VecSetValue(uh,i,value,ADD_VALUES);CHKERRQ(ierr); 4496 } 4497 } 4498 ierr = VecDestroy(&uh);CHKERRQ(ierr); 4499 ierr = VecDestroy(&fh);CHKERRQ(ierr); 4500 PetscFunctionReturn(0); 4501 } 4502 4503 #undef __FUNCT__ 4504 #define __FUNCT__ "SNESKSPSetUseEW" 4505 /*@ 4506 SNESKSPSetUseEW - Sets SNES use Eisenstat-Walker method for 4507 computing relative tolerance for linear solvers within an inexact 4508 Newton method. 4509 4510 Logically Collective on SNES 4511 4512 Input Parameters: 4513 + snes - SNES context 4514 - flag - PETSC_TRUE or PETSC_FALSE 4515 4516 Options Database: 4517 + -snes_ksp_ew - use Eisenstat-Walker method for determining linear system convergence 4518 . -snes_ksp_ew_version ver - version of Eisenstat-Walker method 4519 . -snes_ksp_ew_rtol0 <rtol0> - Sets rtol0 4520 . -snes_ksp_ew_rtolmax <rtolmax> - Sets rtolmax 4521 . -snes_ksp_ew_gamma <gamma> - Sets gamma 4522 . -snes_ksp_ew_alpha <alpha> - Sets alpha 4523 . -snes_ksp_ew_alpha2 <alpha2> - Sets alpha2 4524 - -snes_ksp_ew_threshold <threshold> - Sets threshold 4525 4526 Notes: 4527 Currently, the default is to use a constant relative tolerance for 4528 the inner linear solvers. Alternatively, one can use the 4529 Eisenstat-Walker method, where the relative convergence tolerance 4530 is reset at each Newton iteration according progress of the nonlinear 4531 solver. 4532 4533 Level: advanced 4534 4535 Reference: 4536 S. C. Eisenstat and H. F. Walker, "Choosing the forcing terms in an 4537 inexact Newton method", SISC 17 (1), pp.16-32, 1996. 4538 4539 .keywords: SNES, KSP, Eisenstat, Walker, convergence, test, inexact, Newton 4540 4541 .seealso: SNESKSPGetUseEW(), SNESKSPGetParametersEW(), SNESKSPSetParametersEW() 4542 @*/ 4543 PetscErrorCode SNESKSPSetUseEW(SNES snes,PetscBool flag) 4544 { 4545 PetscFunctionBegin; 4546 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4547 PetscValidLogicalCollectiveBool(snes,flag,2); 4548 snes->ksp_ewconv = flag; 4549 PetscFunctionReturn(0); 4550 } 4551 4552 #undef __FUNCT__ 4553 #define __FUNCT__ "SNESKSPGetUseEW" 4554 /*@ 4555 SNESKSPGetUseEW - Gets if SNES is using Eisenstat-Walker method 4556 for computing relative tolerance for linear solvers within an 4557 inexact Newton method. 4558 4559 Not Collective 4560 4561 Input Parameter: 4562 . snes - SNES context 4563 4564 Output Parameter: 4565 . flag - PETSC_TRUE or PETSC_FALSE 4566 4567 Notes: 4568 Currently, the default is to use a constant relative tolerance for 4569 the inner linear solvers. Alternatively, one can use the 4570 Eisenstat-Walker method, where the relative convergence tolerance 4571 is reset at each Newton iteration according progress of the nonlinear 4572 solver. 4573 4574 Level: advanced 4575 4576 Reference: 4577 S. C. Eisenstat and H. F. Walker, "Choosing the forcing terms in an 4578 inexact Newton method", SISC 17 (1), pp.16-32, 1996. 4579 4580 .keywords: SNES, KSP, Eisenstat, Walker, convergence, test, inexact, Newton 4581 4582 .seealso: SNESKSPSetUseEW(), SNESKSPGetParametersEW(), SNESKSPSetParametersEW() 4583 @*/ 4584 PetscErrorCode SNESKSPGetUseEW(SNES snes, PetscBool *flag) 4585 { 4586 PetscFunctionBegin; 4587 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4588 PetscValidPointer(flag,2); 4589 *flag = snes->ksp_ewconv; 4590 PetscFunctionReturn(0); 4591 } 4592 4593 #undef __FUNCT__ 4594 #define __FUNCT__ "SNESKSPSetParametersEW" 4595 /*@ 4596 SNESKSPSetParametersEW - Sets parameters for Eisenstat-Walker 4597 convergence criteria for the linear solvers within an inexact 4598 Newton method. 4599 4600 Logically Collective on SNES 4601 4602 Input Parameters: 4603 + snes - SNES context 4604 . version - version 1, 2 (default is 2) or 3 4605 . rtol_0 - initial relative tolerance (0 <= rtol_0 < 1) 4606 . rtol_max - maximum relative tolerance (0 <= rtol_max < 1) 4607 . gamma - multiplicative factor for version 2 rtol computation 4608 (0 <= gamma2 <= 1) 4609 . alpha - power for version 2 rtol computation (1 < alpha <= 2) 4610 . alpha2 - power for safeguard 4611 - threshold - threshold for imposing safeguard (0 < threshold < 1) 4612 4613 Note: 4614 Version 3 was contributed by Luis Chacon, June 2006. 4615 4616 Use PETSC_DEFAULT to retain the default for any of the parameters. 4617 4618 Level: advanced 4619 4620 Reference: 4621 S. C. Eisenstat and H. F. Walker, "Choosing the forcing terms in an 4622 inexact Newton method", Utah State University Math. Stat. Dept. Res. 4623 Report 6/94/75, June, 1994, to appear in SIAM J. Sci. Comput. 4624 4625 .keywords: SNES, KSP, Eisenstat, Walker, set, parameters 4626 4627 .seealso: SNESKSPSetUseEW(), SNESKSPGetUseEW(), SNESKSPGetParametersEW() 4628 @*/ 4629 PetscErrorCode SNESKSPSetParametersEW(SNES snes,PetscInt version,PetscReal rtol_0,PetscReal rtol_max,PetscReal gamma,PetscReal alpha,PetscReal alpha2,PetscReal threshold) 4630 { 4631 SNESKSPEW *kctx; 4632 4633 PetscFunctionBegin; 4634 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4635 kctx = (SNESKSPEW*)snes->kspconvctx; 4636 if (!kctx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"No Eisenstat-Walker context existing"); 4637 PetscValidLogicalCollectiveInt(snes,version,2); 4638 PetscValidLogicalCollectiveReal(snes,rtol_0,3); 4639 PetscValidLogicalCollectiveReal(snes,rtol_max,4); 4640 PetscValidLogicalCollectiveReal(snes,gamma,5); 4641 PetscValidLogicalCollectiveReal(snes,alpha,6); 4642 PetscValidLogicalCollectiveReal(snes,alpha2,7); 4643 PetscValidLogicalCollectiveReal(snes,threshold,8); 4644 4645 if (version != PETSC_DEFAULT) kctx->version = version; 4646 if (rtol_0 != PETSC_DEFAULT) kctx->rtol_0 = rtol_0; 4647 if (rtol_max != PETSC_DEFAULT) kctx->rtol_max = rtol_max; 4648 if (gamma != PETSC_DEFAULT) kctx->gamma = gamma; 4649 if (alpha != PETSC_DEFAULT) kctx->alpha = alpha; 4650 if (alpha2 != PETSC_DEFAULT) kctx->alpha2 = alpha2; 4651 if (threshold != PETSC_DEFAULT) kctx->threshold = threshold; 4652 4653 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); 4654 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); 4655 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); 4656 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); 4657 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); 4658 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); 4659 PetscFunctionReturn(0); 4660 } 4661 4662 #undef __FUNCT__ 4663 #define __FUNCT__ "SNESKSPGetParametersEW" 4664 /*@ 4665 SNESKSPGetParametersEW - Gets parameters for Eisenstat-Walker 4666 convergence criteria for the linear solvers within an inexact 4667 Newton method. 4668 4669 Not Collective 4670 4671 Input Parameters: 4672 snes - SNES context 4673 4674 Output Parameters: 4675 + version - version 1, 2 (default is 2) or 3 4676 . rtol_0 - initial relative tolerance (0 <= rtol_0 < 1) 4677 . rtol_max - maximum relative tolerance (0 <= rtol_max < 1) 4678 . gamma - multiplicative factor for version 2 rtol computation (0 <= gamma2 <= 1) 4679 . alpha - power for version 2 rtol computation (1 < alpha <= 2) 4680 . alpha2 - power for safeguard 4681 - threshold - threshold for imposing safeguard (0 < threshold < 1) 4682 4683 Level: advanced 4684 4685 .keywords: SNES, KSP, Eisenstat, Walker, get, parameters 4686 4687 .seealso: SNESKSPSetUseEW(), SNESKSPGetUseEW(), SNESKSPSetParametersEW() 4688 @*/ 4689 PetscErrorCode SNESKSPGetParametersEW(SNES snes,PetscInt *version,PetscReal *rtol_0,PetscReal *rtol_max,PetscReal *gamma,PetscReal *alpha,PetscReal *alpha2,PetscReal *threshold) 4690 { 4691 SNESKSPEW *kctx; 4692 4693 PetscFunctionBegin; 4694 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4695 kctx = (SNESKSPEW*)snes->kspconvctx; 4696 if (!kctx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"No Eisenstat-Walker context existing"); 4697 if (version) *version = kctx->version; 4698 if (rtol_0) *rtol_0 = kctx->rtol_0; 4699 if (rtol_max) *rtol_max = kctx->rtol_max; 4700 if (gamma) *gamma = kctx->gamma; 4701 if (alpha) *alpha = kctx->alpha; 4702 if (alpha2) *alpha2 = kctx->alpha2; 4703 if (threshold) *threshold = kctx->threshold; 4704 PetscFunctionReturn(0); 4705 } 4706 4707 #undef __FUNCT__ 4708 #define __FUNCT__ "KSPPreSolve_SNESEW" 4709 PetscErrorCode KSPPreSolve_SNESEW(KSP ksp, Vec b, Vec x, SNES snes) 4710 { 4711 PetscErrorCode ierr; 4712 SNESKSPEW *kctx = (SNESKSPEW*)snes->kspconvctx; 4713 PetscReal rtol = PETSC_DEFAULT,stol; 4714 4715 PetscFunctionBegin; 4716 if (!snes->ksp_ewconv) PetscFunctionReturn(0); 4717 if (!snes->iter) { 4718 rtol = kctx->rtol_0; /* first time in, so use the original user rtol */ 4719 ierr = VecNorm(snes->vec_func,NORM_2,&kctx->norm_first);CHKERRQ(ierr); 4720 } 4721 else { 4722 if (kctx->version == 1) { 4723 rtol = (snes->norm - kctx->lresid_last)/kctx->norm_last; 4724 if (rtol < 0.0) rtol = -rtol; 4725 stol = PetscPowReal(kctx->rtol_last,kctx->alpha2); 4726 if (stol > kctx->threshold) rtol = PetscMax(rtol,stol); 4727 } else if (kctx->version == 2) { 4728 rtol = kctx->gamma * PetscPowReal(snes->norm/kctx->norm_last,kctx->alpha); 4729 stol = kctx->gamma * PetscPowReal(kctx->rtol_last,kctx->alpha); 4730 if (stol > kctx->threshold) rtol = PetscMax(rtol,stol); 4731 } else if (kctx->version == 3) { /* contributed by Luis Chacon, June 2006. */ 4732 rtol = kctx->gamma * PetscPowReal(snes->norm/kctx->norm_last,kctx->alpha); 4733 /* safeguard: avoid sharp decrease of rtol */ 4734 stol = kctx->gamma*PetscPowReal(kctx->rtol_last,kctx->alpha); 4735 stol = PetscMax(rtol,stol); 4736 rtol = PetscMin(kctx->rtol_0,stol); 4737 /* safeguard: avoid oversolving */ 4738 stol = kctx->gamma*(kctx->norm_first*snes->rtol)/snes->norm; 4739 stol = PetscMax(rtol,stol); 4740 rtol = PetscMin(kctx->rtol_0,stol); 4741 } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Only versions 1, 2 or 3 are supported: %D",kctx->version); 4742 } 4743 /* safeguard: avoid rtol greater than one */ 4744 rtol = PetscMin(rtol,kctx->rtol_max); 4745 ierr = KSPSetTolerances(ksp,rtol,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);CHKERRQ(ierr); 4746 ierr = PetscInfo3(snes,"iter %D, Eisenstat-Walker (version %D) KSP rtol=%g\n",snes->iter,kctx->version,(double)rtol);CHKERRQ(ierr); 4747 PetscFunctionReturn(0); 4748 } 4749 4750 #undef __FUNCT__ 4751 #define __FUNCT__ "KSPPostSolve_SNESEW" 4752 PetscErrorCode KSPPostSolve_SNESEW(KSP ksp, Vec b, Vec x, SNES snes) 4753 { 4754 PetscErrorCode ierr; 4755 SNESKSPEW *kctx = (SNESKSPEW*)snes->kspconvctx; 4756 PCSide pcside; 4757 Vec lres; 4758 4759 PetscFunctionBegin; 4760 if (!snes->ksp_ewconv) PetscFunctionReturn(0); 4761 ierr = KSPGetTolerances(ksp,&kctx->rtol_last,0,0,0);CHKERRQ(ierr); 4762 kctx->norm_last = snes->norm; 4763 if (kctx->version == 1) { 4764 PC pc; 4765 PetscBool isNone; 4766 4767 ierr = KSPGetPC(ksp, &pc);CHKERRQ(ierr); 4768 ierr = PetscObjectTypeCompare((PetscObject) pc, PCNONE, &isNone);CHKERRQ(ierr); 4769 ierr = KSPGetPCSide(ksp,&pcside);CHKERRQ(ierr); 4770 if (pcside == PC_RIGHT || isNone) { /* XXX Should we also test KSP_UNPRECONDITIONED_NORM ? */ 4771 /* KSP residual is true linear residual */ 4772 ierr = KSPGetResidualNorm(ksp,&kctx->lresid_last);CHKERRQ(ierr); 4773 } else { 4774 /* KSP residual is preconditioned residual */ 4775 /* compute true linear residual norm */ 4776 ierr = VecDuplicate(b,&lres);CHKERRQ(ierr); 4777 ierr = MatMult(snes->jacobian,x,lres);CHKERRQ(ierr); 4778 ierr = VecAYPX(lres,-1.0,b);CHKERRQ(ierr); 4779 ierr = VecNorm(lres,NORM_2,&kctx->lresid_last);CHKERRQ(ierr); 4780 ierr = VecDestroy(&lres);CHKERRQ(ierr); 4781 } 4782 } 4783 PetscFunctionReturn(0); 4784 } 4785 4786 #undef __FUNCT__ 4787 #define __FUNCT__ "SNESGetKSP" 4788 /*@ 4789 SNESGetKSP - Returns the KSP context for a SNES solver. 4790 4791 Not Collective, but if SNES object is parallel, then KSP object is parallel 4792 4793 Input Parameter: 4794 . snes - the SNES context 4795 4796 Output Parameter: 4797 . ksp - the KSP context 4798 4799 Notes: 4800 The user can then directly manipulate the KSP context to set various 4801 options, etc. Likewise, the user can then extract and manipulate the 4802 PC contexts as well. 4803 4804 Level: beginner 4805 4806 .keywords: SNES, nonlinear, get, KSP, context 4807 4808 .seealso: KSPGetPC(), SNESCreate(), KSPCreate(), SNESSetKSP() 4809 @*/ 4810 PetscErrorCode SNESGetKSP(SNES snes,KSP *ksp) 4811 { 4812 PetscErrorCode ierr; 4813 4814 PetscFunctionBegin; 4815 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4816 PetscValidPointer(ksp,2); 4817 4818 if (!snes->ksp) { 4819 PetscBool monitor = PETSC_FALSE; 4820 4821 ierr = KSPCreate(PetscObjectComm((PetscObject)snes),&snes->ksp);CHKERRQ(ierr); 4822 ierr = PetscObjectIncrementTabLevel((PetscObject)snes->ksp,(PetscObject)snes,1);CHKERRQ(ierr); 4823 ierr = PetscLogObjectParent((PetscObject)snes,(PetscObject)snes->ksp);CHKERRQ(ierr); 4824 4825 ierr = KSPSetPreSolve(snes->ksp,(PetscErrorCode (*)(KSP,Vec,Vec,void*))KSPPreSolve_SNESEW,snes);CHKERRQ(ierr); 4826 ierr = KSPSetPostSolve(snes->ksp,(PetscErrorCode (*)(KSP,Vec,Vec,void*))KSPPostSolve_SNESEW,snes);CHKERRQ(ierr); 4827 4828 ierr = PetscOptionsGetBool(((PetscObject)snes)->options,((PetscObject)snes)->prefix,"-ksp_monitor_snes",&monitor,NULL);CHKERRQ(ierr); 4829 if (monitor) { 4830 ierr = KSPMonitorSet(snes->ksp,KSPMonitorSNES,snes,NULL);CHKERRQ(ierr); 4831 } 4832 monitor = PETSC_FALSE; 4833 ierr = PetscOptionsGetBool(((PetscObject)snes)->options,((PetscObject)snes)->prefix,"-ksp_monitor_snes_lg",&monitor,NULL);CHKERRQ(ierr); 4834 if (monitor) { 4835 PetscObject *objs; 4836 ierr = KSPMonitorSNESLGResidualNormCreate(PetscObjectComm((PetscObject)snes),NULL,NULL,PETSC_DECIDE,PETSC_DECIDE,600,600,&objs);CHKERRQ(ierr); 4837 objs[0] = (PetscObject) snes; 4838 ierr = KSPMonitorSet(snes->ksp,(PetscErrorCode (*)(KSP,PetscInt,PetscReal,void*))KSPMonitorSNESLGResidualNorm,objs,(PetscErrorCode (*)(void**))KSPMonitorSNESLGResidualNormDestroy);CHKERRQ(ierr); 4839 } 4840 } 4841 *ksp = snes->ksp; 4842 PetscFunctionReturn(0); 4843 } 4844 4845 4846 #include <petsc/private/dmimpl.h> 4847 #undef __FUNCT__ 4848 #define __FUNCT__ "SNESSetDM" 4849 /*@ 4850 SNESSetDM - Sets the DM that may be used by some preconditioners 4851 4852 Logically Collective on SNES 4853 4854 Input Parameters: 4855 + snes - the preconditioner context 4856 - dm - the dm 4857 4858 Level: intermediate 4859 4860 .seealso: SNESGetDM(), KSPSetDM(), KSPGetDM() 4861 @*/ 4862 PetscErrorCode SNESSetDM(SNES snes,DM dm) 4863 { 4864 PetscErrorCode ierr; 4865 KSP ksp; 4866 DMSNES sdm; 4867 4868 PetscFunctionBegin; 4869 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4870 if (dm) {ierr = PetscObjectReference((PetscObject)dm);CHKERRQ(ierr);} 4871 if (snes->dm) { /* Move the DMSNES context over to the new DM unless the new DM already has one */ 4872 if (snes->dm->dmsnes && snes->dmAuto && !dm->dmsnes) { 4873 ierr = DMCopyDMSNES(snes->dm,dm);CHKERRQ(ierr); 4874 ierr = DMGetDMSNES(snes->dm,&sdm);CHKERRQ(ierr); 4875 if (sdm->originaldm == snes->dm) sdm->originaldm = dm; /* Grant write privileges to the replacement DM */ 4876 } 4877 ierr = DMDestroy(&snes->dm);CHKERRQ(ierr); 4878 } 4879 snes->dm = dm; 4880 snes->dmAuto = PETSC_FALSE; 4881 4882 ierr = SNESGetKSP(snes,&ksp);CHKERRQ(ierr); 4883 ierr = KSPSetDM(ksp,dm);CHKERRQ(ierr); 4884 ierr = KSPSetDMActive(ksp,PETSC_FALSE);CHKERRQ(ierr); 4885 if (snes->pc) { 4886 ierr = SNESSetDM(snes->pc, snes->dm);CHKERRQ(ierr); 4887 ierr = SNESSetNPCSide(snes,snes->pcside);CHKERRQ(ierr); 4888 } 4889 PetscFunctionReturn(0); 4890 } 4891 4892 #undef __FUNCT__ 4893 #define __FUNCT__ "SNESGetDM" 4894 /*@ 4895 SNESGetDM - Gets the DM that may be used by some preconditioners 4896 4897 Not Collective but DM obtained is parallel on SNES 4898 4899 Input Parameter: 4900 . snes - the preconditioner context 4901 4902 Output Parameter: 4903 . dm - the dm 4904 4905 Level: intermediate 4906 4907 .seealso: SNESSetDM(), KSPSetDM(), KSPGetDM() 4908 @*/ 4909 PetscErrorCode SNESGetDM(SNES snes,DM *dm) 4910 { 4911 PetscErrorCode ierr; 4912 4913 PetscFunctionBegin; 4914 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4915 if (!snes->dm) { 4916 ierr = DMShellCreate(PetscObjectComm((PetscObject)snes),&snes->dm);CHKERRQ(ierr); 4917 snes->dmAuto = PETSC_TRUE; 4918 } 4919 *dm = snes->dm; 4920 PetscFunctionReturn(0); 4921 } 4922 4923 #undef __FUNCT__ 4924 #define __FUNCT__ "SNESSetNPC" 4925 /*@ 4926 SNESSetNPC - Sets the nonlinear preconditioner to be used. 4927 4928 Collective on SNES 4929 4930 Input Parameters: 4931 + snes - iterative context obtained from SNESCreate() 4932 - pc - the preconditioner object 4933 4934 Notes: 4935 Use SNESGetNPC() to retrieve the preconditioner context (for example, 4936 to configure it using the API). 4937 4938 Level: developer 4939 4940 .keywords: SNES, set, precondition 4941 .seealso: SNESGetNPC(), SNESHasNPC() 4942 @*/ 4943 PetscErrorCode SNESSetNPC(SNES snes, SNES pc) 4944 { 4945 PetscErrorCode ierr; 4946 4947 PetscFunctionBegin; 4948 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); 4949 PetscValidHeaderSpecific(pc, SNES_CLASSID, 2); 4950 PetscCheckSameComm(snes, 1, pc, 2); 4951 ierr = PetscObjectReference((PetscObject) pc);CHKERRQ(ierr); 4952 ierr = SNESDestroy(&snes->pc);CHKERRQ(ierr); 4953 snes->pc = pc; 4954 ierr = PetscLogObjectParent((PetscObject)snes, (PetscObject)snes->pc);CHKERRQ(ierr); 4955 PetscFunctionReturn(0); 4956 } 4957 4958 #undef __FUNCT__ 4959 #define __FUNCT__ "SNESGetNPC" 4960 /*@ 4961 SNESGetNPC - Creates a nonlinear preconditioning solver (SNES) to be used to precondition the nonlinear solver. 4962 4963 Not Collective 4964 4965 Input Parameter: 4966 . snes - iterative context obtained from SNESCreate() 4967 4968 Output Parameter: 4969 . pc - preconditioner context 4970 4971 Notes: If a SNES was previously set with SNESSetNPC() then that SNES is returned. 4972 4973 Level: developer 4974 4975 .keywords: SNES, get, preconditioner 4976 .seealso: SNESSetNPC(), SNESHasNPC() 4977 @*/ 4978 PetscErrorCode SNESGetNPC(SNES snes, SNES *pc) 4979 { 4980 PetscErrorCode ierr; 4981 const char *optionsprefix; 4982 4983 PetscFunctionBegin; 4984 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); 4985 PetscValidPointer(pc, 2); 4986 if (!snes->pc) { 4987 ierr = SNESCreate(PetscObjectComm((PetscObject)snes),&snes->pc);CHKERRQ(ierr); 4988 ierr = PetscObjectIncrementTabLevel((PetscObject)snes->pc,(PetscObject)snes,1);CHKERRQ(ierr); 4989 ierr = PetscLogObjectParent((PetscObject)snes,(PetscObject)snes->pc);CHKERRQ(ierr); 4990 ierr = SNESGetOptionsPrefix(snes,&optionsprefix);CHKERRQ(ierr); 4991 ierr = SNESSetOptionsPrefix(snes->pc,optionsprefix);CHKERRQ(ierr); 4992 ierr = SNESAppendOptionsPrefix(snes->pc,"npc_");CHKERRQ(ierr); 4993 ierr = SNESSetCountersReset(snes->pc,PETSC_FALSE);CHKERRQ(ierr); 4994 } 4995 *pc = snes->pc; 4996 PetscFunctionReturn(0); 4997 } 4998 4999 #undef __FUNCT__ 5000 #define __FUNCT__ "SNESHasNPC" 5001 /*@ 5002 SNESHasNPC - Returns whether a nonlinear preconditioner exists 5003 5004 Not Collective 5005 5006 Input Parameter: 5007 . snes - iterative context obtained from SNESCreate() 5008 5009 Output Parameter: 5010 . has_npc - whether the SNES has an NPC or not 5011 5012 Level: developer 5013 5014 .keywords: SNES, has, preconditioner 5015 .seealso: SNESSetNPC(), SNESGetNPC() 5016 @*/ 5017 PetscErrorCode SNESHasNPC(SNES snes, PetscBool *has_npc) 5018 { 5019 PetscFunctionBegin; 5020 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); 5021 *has_npc = (PetscBool) (snes->pc ? PETSC_TRUE : PETSC_FALSE); 5022 PetscFunctionReturn(0); 5023 } 5024 5025 #undef __FUNCT__ 5026 #define __FUNCT__ "SNESSetNPCSide" 5027 /*@ 5028 SNESSetNPCSide - Sets the preconditioning side. 5029 5030 Logically Collective on SNES 5031 5032 Input Parameter: 5033 . snes - iterative context obtained from SNESCreate() 5034 5035 Output Parameter: 5036 . side - the preconditioning side, where side is one of 5037 .vb 5038 PC_LEFT - left preconditioning (default) 5039 PC_RIGHT - right preconditioning 5040 .ve 5041 5042 Options Database Keys: 5043 . -snes_pc_side <right,left> 5044 5045 Level: intermediate 5046 5047 .keywords: SNES, set, right, left, side, preconditioner, flag 5048 5049 .seealso: SNESGetNPCSide(), KSPSetPCSide() 5050 @*/ 5051 PetscErrorCode SNESSetNPCSide(SNES snes,PCSide side) 5052 { 5053 PetscFunctionBegin; 5054 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 5055 PetscValidLogicalCollectiveEnum(snes,side,2); 5056 snes->pcside = side; 5057 PetscFunctionReturn(0); 5058 } 5059 5060 #undef __FUNCT__ 5061 #define __FUNCT__ "SNESGetNPCSide" 5062 /*@ 5063 SNESGetNPCSide - Gets the preconditioning side. 5064 5065 Not Collective 5066 5067 Input Parameter: 5068 . snes - iterative context obtained from SNESCreate() 5069 5070 Output Parameter: 5071 . side - the preconditioning side, where side is one of 5072 .vb 5073 PC_LEFT - left preconditioning (default) 5074 PC_RIGHT - right preconditioning 5075 .ve 5076 5077 Level: intermediate 5078 5079 .keywords: SNES, get, right, left, side, preconditioner, flag 5080 5081 .seealso: SNESSetNPCSide(), KSPGetPCSide() 5082 @*/ 5083 PetscErrorCode SNESGetNPCSide(SNES snes,PCSide *side) 5084 { 5085 PetscFunctionBegin; 5086 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 5087 PetscValidPointer(side,2); 5088 *side = snes->pcside; 5089 PetscFunctionReturn(0); 5090 } 5091 5092 #undef __FUNCT__ 5093 #define __FUNCT__ "SNESSetLineSearch" 5094 /*@ 5095 SNESSetLineSearch - Sets the linesearch on the SNES instance. 5096 5097 Collective on SNES 5098 5099 Input Parameters: 5100 + snes - iterative context obtained from SNESCreate() 5101 - linesearch - the linesearch object 5102 5103 Notes: 5104 Use SNESGetLineSearch() to retrieve the preconditioner context (for example, 5105 to configure it using the API). 5106 5107 Level: developer 5108 5109 .keywords: SNES, set, linesearch 5110 .seealso: SNESGetLineSearch() 5111 @*/ 5112 PetscErrorCode SNESSetLineSearch(SNES snes, SNESLineSearch linesearch) 5113 { 5114 PetscErrorCode ierr; 5115 5116 PetscFunctionBegin; 5117 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); 5118 PetscValidHeaderSpecific(linesearch, SNESLINESEARCH_CLASSID, 2); 5119 PetscCheckSameComm(snes, 1, linesearch, 2); 5120 ierr = PetscObjectReference((PetscObject) linesearch);CHKERRQ(ierr); 5121 ierr = SNESLineSearchDestroy(&snes->linesearch);CHKERRQ(ierr); 5122 5123 snes->linesearch = linesearch; 5124 5125 ierr = PetscLogObjectParent((PetscObject)snes, (PetscObject)snes->linesearch);CHKERRQ(ierr); 5126 PetscFunctionReturn(0); 5127 } 5128 5129 #undef __FUNCT__ 5130 #define __FUNCT__ "SNESGetLineSearch" 5131 /*@ 5132 SNESGetLineSearch - Returns a pointer to the line search context set with SNESSetLineSearch() 5133 or creates a default line search instance associated with the SNES and returns it. 5134 5135 Not Collective 5136 5137 Input Parameter: 5138 . snes - iterative context obtained from SNESCreate() 5139 5140 Output Parameter: 5141 . linesearch - linesearch context 5142 5143 Level: beginner 5144 5145 .keywords: SNES, get, linesearch 5146 .seealso: SNESSetLineSearch(), SNESLineSearchCreate() 5147 @*/ 5148 PetscErrorCode SNESGetLineSearch(SNES snes, SNESLineSearch *linesearch) 5149 { 5150 PetscErrorCode ierr; 5151 const char *optionsprefix; 5152 5153 PetscFunctionBegin; 5154 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); 5155 PetscValidPointer(linesearch, 2); 5156 if (!snes->linesearch) { 5157 ierr = SNESGetOptionsPrefix(snes, &optionsprefix);CHKERRQ(ierr); 5158 ierr = SNESLineSearchCreate(PetscObjectComm((PetscObject)snes), &snes->linesearch);CHKERRQ(ierr); 5159 ierr = SNESLineSearchSetSNES(snes->linesearch, snes);CHKERRQ(ierr); 5160 ierr = SNESLineSearchAppendOptionsPrefix(snes->linesearch, optionsprefix);CHKERRQ(ierr); 5161 ierr = PetscObjectIncrementTabLevel((PetscObject) snes->linesearch, (PetscObject) snes, 1);CHKERRQ(ierr); 5162 ierr = PetscLogObjectParent((PetscObject)snes, (PetscObject)snes->linesearch);CHKERRQ(ierr); 5163 } 5164 *linesearch = snes->linesearch; 5165 PetscFunctionReturn(0); 5166 } 5167 5168 #if defined(PETSC_HAVE_MATLAB_ENGINE) 5169 #include <mex.h> 5170 5171 typedef struct {char *funcname; mxArray *ctx;} SNESMatlabContext; 5172 5173 #undef __FUNCT__ 5174 #define __FUNCT__ "SNESComputeFunction_Matlab" 5175 /* 5176 SNESComputeFunction_Matlab - Calls the function that has been set with SNESSetFunctionMatlab(). 5177 5178 Collective on SNES 5179 5180 Input Parameters: 5181 + snes - the SNES context 5182 - x - input vector 5183 5184 Output Parameter: 5185 . y - function vector, as set by SNESSetFunction() 5186 5187 Notes: 5188 SNESComputeFunction() is typically used within nonlinear solvers 5189 implementations, so most users would not generally call this routine 5190 themselves. 5191 5192 Level: developer 5193 5194 .keywords: SNES, nonlinear, compute, function 5195 5196 .seealso: SNESSetFunction(), SNESGetFunction() 5197 */ 5198 PetscErrorCode SNESComputeFunction_Matlab(SNES snes,Vec x,Vec y, void *ctx) 5199 { 5200 PetscErrorCode ierr; 5201 SNESMatlabContext *sctx = (SNESMatlabContext*)ctx; 5202 int nlhs = 1,nrhs = 5; 5203 mxArray *plhs[1],*prhs[5]; 5204 long long int lx = 0,ly = 0,ls = 0; 5205 5206 PetscFunctionBegin; 5207 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 5208 PetscValidHeaderSpecific(x,VEC_CLASSID,2); 5209 PetscValidHeaderSpecific(y,VEC_CLASSID,3); 5210 PetscCheckSameComm(snes,1,x,2); 5211 PetscCheckSameComm(snes,1,y,3); 5212 5213 /* call Matlab function in ctx with arguments x and y */ 5214 5215 ierr = PetscMemcpy(&ls,&snes,sizeof(snes));CHKERRQ(ierr); 5216 ierr = PetscMemcpy(&lx,&x,sizeof(x));CHKERRQ(ierr); 5217 ierr = PetscMemcpy(&ly,&y,sizeof(x));CHKERRQ(ierr); 5218 prhs[0] = mxCreateDoubleScalar((double)ls); 5219 prhs[1] = mxCreateDoubleScalar((double)lx); 5220 prhs[2] = mxCreateDoubleScalar((double)ly); 5221 prhs[3] = mxCreateString(sctx->funcname); 5222 prhs[4] = sctx->ctx; 5223 ierr = mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscSNESComputeFunctionInternal");CHKERRQ(ierr); 5224 ierr = mxGetScalar(plhs[0]);CHKERRQ(ierr); 5225 mxDestroyArray(prhs[0]); 5226 mxDestroyArray(prhs[1]); 5227 mxDestroyArray(prhs[2]); 5228 mxDestroyArray(prhs[3]); 5229 mxDestroyArray(plhs[0]); 5230 PetscFunctionReturn(0); 5231 } 5232 5233 #undef __FUNCT__ 5234 #define __FUNCT__ "SNESSetFunctionMatlab" 5235 /* 5236 SNESSetFunctionMatlab - Sets the function evaluation routine and function 5237 vector for use by the SNES routines in solving systems of nonlinear 5238 equations from MATLAB. Here the function is a string containing the name of a MATLAB function 5239 5240 Logically Collective on SNES 5241 5242 Input Parameters: 5243 + snes - the SNES context 5244 . r - vector to store function value 5245 - f - function evaluation routine 5246 5247 Notes: 5248 The Newton-like methods typically solve linear systems of the form 5249 $ f'(x) x = -f(x), 5250 where f'(x) denotes the Jacobian matrix and f(x) is the function. 5251 5252 Level: beginner 5253 5254 Developer Note: This bleeds the allocated memory SNESMatlabContext *sctx; 5255 5256 .keywords: SNES, nonlinear, set, function 5257 5258 .seealso: SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESSetFunction() 5259 */ 5260 PetscErrorCode SNESSetFunctionMatlab(SNES snes,Vec r,const char *f,mxArray *ctx) 5261 { 5262 PetscErrorCode ierr; 5263 SNESMatlabContext *sctx; 5264 5265 PetscFunctionBegin; 5266 /* currently sctx is memory bleed */ 5267 ierr = PetscNew(&sctx);CHKERRQ(ierr); 5268 ierr = PetscStrallocpy(f,&sctx->funcname);CHKERRQ(ierr); 5269 /* 5270 This should work, but it doesn't 5271 sctx->ctx = ctx; 5272 mexMakeArrayPersistent(sctx->ctx); 5273 */ 5274 sctx->ctx = mxDuplicateArray(ctx); 5275 ierr = SNESSetFunction(snes,r,SNESComputeFunction_Matlab,sctx);CHKERRQ(ierr); 5276 PetscFunctionReturn(0); 5277 } 5278 5279 #undef __FUNCT__ 5280 #define __FUNCT__ "SNESComputeJacobian_Matlab" 5281 /* 5282 SNESComputeJacobian_Matlab - Calls the function that has been set with SNESSetJacobianMatlab(). 5283 5284 Collective on SNES 5285 5286 Input Parameters: 5287 + snes - the SNES context 5288 . x - input vector 5289 . A, B - the matrices 5290 - ctx - user context 5291 5292 Level: developer 5293 5294 .keywords: SNES, nonlinear, compute, function 5295 5296 .seealso: SNESSetFunction(), SNESGetFunction() 5297 @*/ 5298 PetscErrorCode SNESComputeJacobian_Matlab(SNES snes,Vec x,Mat A,Mat B,void *ctx) 5299 { 5300 PetscErrorCode ierr; 5301 SNESMatlabContext *sctx = (SNESMatlabContext*)ctx; 5302 int nlhs = 2,nrhs = 6; 5303 mxArray *plhs[2],*prhs[6]; 5304 long long int lx = 0,lA = 0,ls = 0, lB = 0; 5305 5306 PetscFunctionBegin; 5307 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 5308 PetscValidHeaderSpecific(x,VEC_CLASSID,2); 5309 5310 /* call Matlab function in ctx with arguments x and y */ 5311 5312 ierr = PetscMemcpy(&ls,&snes,sizeof(snes));CHKERRQ(ierr); 5313 ierr = PetscMemcpy(&lx,&x,sizeof(x));CHKERRQ(ierr); 5314 ierr = PetscMemcpy(&lA,A,sizeof(x));CHKERRQ(ierr); 5315 ierr = PetscMemcpy(&lB,B,sizeof(x));CHKERRQ(ierr); 5316 prhs[0] = mxCreateDoubleScalar((double)ls); 5317 prhs[1] = mxCreateDoubleScalar((double)lx); 5318 prhs[2] = mxCreateDoubleScalar((double)lA); 5319 prhs[3] = mxCreateDoubleScalar((double)lB); 5320 prhs[4] = mxCreateString(sctx->funcname); 5321 prhs[5] = sctx->ctx; 5322 ierr = mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscSNESComputeJacobianInternal");CHKERRQ(ierr); 5323 ierr = mxGetScalar(plhs[0]);CHKERRQ(ierr); 5324 mxDestroyArray(prhs[0]); 5325 mxDestroyArray(prhs[1]); 5326 mxDestroyArray(prhs[2]); 5327 mxDestroyArray(prhs[3]); 5328 mxDestroyArray(prhs[4]); 5329 mxDestroyArray(plhs[0]); 5330 mxDestroyArray(plhs[1]); 5331 PetscFunctionReturn(0); 5332 } 5333 5334 #undef __FUNCT__ 5335 #define __FUNCT__ "SNESSetJacobianMatlab" 5336 /* 5337 SNESSetJacobianMatlab - Sets the Jacobian function evaluation routine and two empty Jacobian matrices 5338 vector for use by the SNES routines in solving systems of nonlinear 5339 equations from MATLAB. Here the function is a string containing the name of a MATLAB function 5340 5341 Logically Collective on SNES 5342 5343 Input Parameters: 5344 + snes - the SNES context 5345 . A,B - Jacobian matrices 5346 . J - function evaluation routine 5347 - ctx - user context 5348 5349 Level: developer 5350 5351 Developer Note: This bleeds the allocated memory SNESMatlabContext *sctx; 5352 5353 .keywords: SNES, nonlinear, set, function 5354 5355 .seealso: SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESSetFunction(), J 5356 */ 5357 PetscErrorCode SNESSetJacobianMatlab(SNES snes,Mat A,Mat B,const char *J,mxArray *ctx) 5358 { 5359 PetscErrorCode ierr; 5360 SNESMatlabContext *sctx; 5361 5362 PetscFunctionBegin; 5363 /* currently sctx is memory bleed */ 5364 ierr = PetscNew(&sctx);CHKERRQ(ierr); 5365 ierr = PetscStrallocpy(J,&sctx->funcname);CHKERRQ(ierr); 5366 /* 5367 This should work, but it doesn't 5368 sctx->ctx = ctx; 5369 mexMakeArrayPersistent(sctx->ctx); 5370 */ 5371 sctx->ctx = mxDuplicateArray(ctx); 5372 ierr = SNESSetJacobian(snes,A,B,SNESComputeJacobian_Matlab,sctx);CHKERRQ(ierr); 5373 PetscFunctionReturn(0); 5374 } 5375 5376 #undef __FUNCT__ 5377 #define __FUNCT__ "SNESMonitor_Matlab" 5378 /* 5379 SNESMonitor_Matlab - Calls the function that has been set with SNESMonitorSetMatlab(). 5380 5381 Collective on SNES 5382 5383 .seealso: SNESSetFunction(), SNESGetFunction() 5384 @*/ 5385 PetscErrorCode SNESMonitor_Matlab(SNES snes,PetscInt it, PetscReal fnorm, void *ctx) 5386 { 5387 PetscErrorCode ierr; 5388 SNESMatlabContext *sctx = (SNESMatlabContext*)ctx; 5389 int nlhs = 1,nrhs = 6; 5390 mxArray *plhs[1],*prhs[6]; 5391 long long int lx = 0,ls = 0; 5392 Vec x = snes->vec_sol; 5393 5394 PetscFunctionBegin; 5395 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 5396 5397 ierr = PetscMemcpy(&ls,&snes,sizeof(snes));CHKERRQ(ierr); 5398 ierr = PetscMemcpy(&lx,&x,sizeof(x));CHKERRQ(ierr); 5399 prhs[0] = mxCreateDoubleScalar((double)ls); 5400 prhs[1] = mxCreateDoubleScalar((double)it); 5401 prhs[2] = mxCreateDoubleScalar((double)fnorm); 5402 prhs[3] = mxCreateDoubleScalar((double)lx); 5403 prhs[4] = mxCreateString(sctx->funcname); 5404 prhs[5] = sctx->ctx; 5405 ierr = mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscSNESMonitorInternal");CHKERRQ(ierr); 5406 ierr = mxGetScalar(plhs[0]);CHKERRQ(ierr); 5407 mxDestroyArray(prhs[0]); 5408 mxDestroyArray(prhs[1]); 5409 mxDestroyArray(prhs[2]); 5410 mxDestroyArray(prhs[3]); 5411 mxDestroyArray(prhs[4]); 5412 mxDestroyArray(plhs[0]); 5413 PetscFunctionReturn(0); 5414 } 5415 5416 #undef __FUNCT__ 5417 #define __FUNCT__ "SNESMonitorSetMatlab" 5418 /* 5419 SNESMonitorSetMatlab - Sets the monitor function from MATLAB 5420 5421 Level: developer 5422 5423 Developer Note: This bleeds the allocated memory SNESMatlabContext *sctx; 5424 5425 .keywords: SNES, nonlinear, set, function 5426 5427 .seealso: SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESSetFunction() 5428 */ 5429 PetscErrorCode SNESMonitorSetMatlab(SNES snes,const char *f,mxArray *ctx) 5430 { 5431 PetscErrorCode ierr; 5432 SNESMatlabContext *sctx; 5433 5434 PetscFunctionBegin; 5435 /* currently sctx is memory bleed */ 5436 ierr = PetscNew(&sctx);CHKERRQ(ierr); 5437 ierr = PetscStrallocpy(f,&sctx->funcname);CHKERRQ(ierr); 5438 /* 5439 This should work, but it doesn't 5440 sctx->ctx = ctx; 5441 mexMakeArrayPersistent(sctx->ctx); 5442 */ 5443 sctx->ctx = mxDuplicateArray(ctx); 5444 ierr = SNESMonitorSet(snes,SNESMonitor_Matlab,sctx,NULL);CHKERRQ(ierr); 5445 PetscFunctionReturn(0); 5446 } 5447 5448 #endif 5449