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