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