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