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,&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,&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,SNESConvergedSkip,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 ierr = SNESInitializePackage();CHKERRQ(ierr); 1466 1467 ierr = PetscHeaderCreate(snes,_p_SNES,struct _SNESOps,SNES_CLASSID,"SNES","Nonlinear solver","SNES",comm,SNESDestroy,SNESView);CHKERRQ(ierr); 1468 1469 snes->ops->converged = SNESConvergedDefault; 1470 snes->usesksp = PETSC_TRUE; 1471 snes->tolerancesset = PETSC_FALSE; 1472 snes->max_its = 50; 1473 snes->max_funcs = 10000; 1474 snes->norm = 0.0; 1475 snes->normschedule = SNES_NORM_ALWAYS; 1476 snes->functype = SNES_FUNCTION_DEFAULT; 1477 snes->rtol = 1.e-8; 1478 snes->ttol = 0.0; 1479 snes->abstol = 1.e-50; 1480 snes->stol = 1.e-8; 1481 snes->deltatol = 1.e-12; 1482 snes->nfuncs = 0; 1483 snes->numFailures = 0; 1484 snes->maxFailures = 1; 1485 snes->linear_its = 0; 1486 snes->lagjacobian = 1; 1487 snes->jac_iter = 0; 1488 snes->lagjac_persist = PETSC_FALSE; 1489 snes->lagpreconditioner = 1; 1490 snes->pre_iter = 0; 1491 snes->lagpre_persist = PETSC_FALSE; 1492 snes->numbermonitors = 0; 1493 snes->data = 0; 1494 snes->setupcalled = PETSC_FALSE; 1495 snes->ksp_ewconv = PETSC_FALSE; 1496 snes->nwork = 0; 1497 snes->work = 0; 1498 snes->nvwork = 0; 1499 snes->vwork = 0; 1500 snes->conv_hist_len = 0; 1501 snes->conv_hist_max = 0; 1502 snes->conv_hist = NULL; 1503 snes->conv_hist_its = NULL; 1504 snes->conv_hist_reset = PETSC_TRUE; 1505 snes->counters_reset = PETSC_TRUE; 1506 snes->vec_func_init_set = PETSC_FALSE; 1507 snes->reason = SNES_CONVERGED_ITERATING; 1508 snes->pcside = PC_RIGHT; 1509 1510 snes->mf = PETSC_FALSE; 1511 snes->mf_operator = PETSC_FALSE; 1512 snes->mf_version = 1; 1513 1514 snes->numLinearSolveFailures = 0; 1515 snes->maxLinearSolveFailures = 1; 1516 1517 /* Create context to compute Eisenstat-Walker relative tolerance for KSP */ 1518 ierr = PetscNewLog(snes,SNESKSPEW,&kctx);CHKERRQ(ierr); 1519 1520 snes->kspconvctx = (void*)kctx; 1521 kctx->version = 2; 1522 kctx->rtol_0 = .3; /* Eisenstat and Walker suggest rtol_0=.5, but 1523 this was too large for some test cases */ 1524 kctx->rtol_last = 0.0; 1525 kctx->rtol_max = .9; 1526 kctx->gamma = 1.0; 1527 kctx->alpha = .5*(1.0 + PetscSqrtReal(5.0)); 1528 kctx->alpha2 = kctx->alpha; 1529 kctx->threshold = .1; 1530 kctx->lresid_last = 0.0; 1531 kctx->norm_last = 0.0; 1532 1533 *outsnes = snes; 1534 PetscFunctionReturn(0); 1535 } 1536 1537 /*MC 1538 SNESFunction - functional form used to convey the nonlinear function to be solved by SNES 1539 1540 Synopsis: 1541 #include "petscsnes.h" 1542 SNESFunction(SNES snes,Vec x,Vec f,void *ctx); 1543 1544 Input Parameters: 1545 + snes - the SNES context 1546 . x - state at which to evaluate residual 1547 - ctx - optional user-defined function context, passed in with SNESSetFunction() 1548 1549 Output Parameter: 1550 . f - vector to put residual (function value) 1551 1552 Level: intermediate 1553 1554 .seealso: SNESSetFunction(), SNESGetFunction() 1555 M*/ 1556 1557 #undef __FUNCT__ 1558 #define __FUNCT__ "SNESSetFunction" 1559 /*@C 1560 SNESSetFunction - Sets the function evaluation routine and function 1561 vector for use by the SNES routines in solving systems of nonlinear 1562 equations. 1563 1564 Logically Collective on SNES 1565 1566 Input Parameters: 1567 + snes - the SNES context 1568 . r - vector to store function value 1569 . SNESFunction - function evaluation routine 1570 - ctx - [optional] user-defined context for private data for the 1571 function evaluation routine (may be NULL) 1572 1573 Notes: 1574 The Newton-like methods typically solve linear systems of the form 1575 $ f'(x) x = -f(x), 1576 where f'(x) denotes the Jacobian matrix and f(x) is the function. 1577 1578 Level: beginner 1579 1580 .keywords: SNES, nonlinear, set, function 1581 1582 .seealso: SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESSetPicard(), SNESFunction 1583 @*/ 1584 PetscErrorCode SNESSetFunction(SNES snes,Vec r,PetscErrorCode (*SNESFunction)(SNES,Vec,Vec,void*),void *ctx) 1585 { 1586 PetscErrorCode ierr; 1587 DM dm; 1588 1589 PetscFunctionBegin; 1590 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1591 if (r) { 1592 PetscValidHeaderSpecific(r,VEC_CLASSID,2); 1593 PetscCheckSameComm(snes,1,r,2); 1594 ierr = PetscObjectReference((PetscObject)r);CHKERRQ(ierr); 1595 ierr = VecDestroy(&snes->vec_func);CHKERRQ(ierr); 1596 1597 snes->vec_func = r; 1598 } 1599 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 1600 ierr = DMSNESSetFunction(dm,SNESFunction,ctx);CHKERRQ(ierr); 1601 PetscFunctionReturn(0); 1602 } 1603 1604 1605 #undef __FUNCT__ 1606 #define __FUNCT__ "SNESSetInitialFunction" 1607 /*@C 1608 SNESSetInitialFunction - Sets the function vector to be used as the 1609 function norm at the initialization of the method. In some 1610 instances, the user has precomputed the function before calling 1611 SNESSolve. This function allows one to avoid a redundant call 1612 to SNESComputeFunction in that case. 1613 1614 Logically Collective on SNES 1615 1616 Input Parameters: 1617 + snes - the SNES context 1618 - f - vector to store function value 1619 1620 Notes: 1621 This should not be modified during the solution procedure. 1622 1623 This is used extensively in the SNESFAS hierarchy and in nonlinear preconditioning. 1624 1625 Level: developer 1626 1627 .keywords: SNES, nonlinear, set, function 1628 1629 .seealso: SNESSetFunction(), SNESComputeFunction(), SNESSetInitialFunctionNorm() 1630 @*/ 1631 PetscErrorCode SNESSetInitialFunction(SNES snes, Vec f) 1632 { 1633 PetscErrorCode ierr; 1634 Vec vec_func; 1635 1636 PetscFunctionBegin; 1637 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1638 PetscValidHeaderSpecific(f,VEC_CLASSID,2); 1639 PetscCheckSameComm(snes,1,f,2); 1640 if (snes->pcside == PC_LEFT && snes->functype == SNES_FUNCTION_PRECONDITIONED) { 1641 snes->vec_func_init_set = PETSC_FALSE; 1642 PetscFunctionReturn(0); 1643 } 1644 ierr = SNESGetFunction(snes,&vec_func,NULL,NULL);CHKERRQ(ierr); 1645 ierr = VecCopy(f, vec_func);CHKERRQ(ierr); 1646 1647 snes->vec_func_init_set = PETSC_TRUE; 1648 PetscFunctionReturn(0); 1649 } 1650 1651 #undef __FUNCT__ 1652 #define __FUNCT__ "SNESSetNormSchedule" 1653 /*@ 1654 SNESSetNormSchedule - Sets the SNESNormSchedule used in covergence and monitoring 1655 of the SNES method. 1656 1657 Logically Collective on SNES 1658 1659 Input Parameters: 1660 + snes - the SNES context 1661 - normschedule - the frequency of norm computation 1662 1663 Notes: 1664 Only certain SNES methods support certain SNESNormSchedules. Most require evaluation 1665 of the nonlinear function and the taking of its norm at every iteration to 1666 even ensure convergence at all. However, methods such as custom Gauss-Seidel methods 1667 (SNESGS) and the like do not require the norm of the function to be computed, and therfore 1668 may either be monitored for convergence or not. As these are often used as nonlinear 1669 preconditioners, monitoring the norm of their error is not a useful enterprise within 1670 their solution. 1671 1672 Level: developer 1673 1674 .keywords: SNES, nonlinear, set, function, norm, type 1675 1676 .seealso: SNESGetNormSchedule(), SNESComputeFunction(), VecNorm(), SNESSetFunction(), SNESSetInitialFunction(), SNESNormSchedule 1677 @*/ 1678 PetscErrorCode SNESSetNormSchedule(SNES snes, SNESNormSchedule normschedule) 1679 { 1680 PetscFunctionBegin; 1681 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1682 snes->normschedule = normschedule; 1683 PetscFunctionReturn(0); 1684 } 1685 1686 1687 #undef __FUNCT__ 1688 #define __FUNCT__ "SNESGetNormSchedule" 1689 /*@ 1690 SNESGetNormSchedule - Gets the SNESNormSchedule used in covergence and monitoring 1691 of the SNES method. 1692 1693 Logically Collective on SNES 1694 1695 Input Parameters: 1696 + snes - the SNES context 1697 - normschedule - the type of the norm used 1698 1699 Level: advanced 1700 1701 .keywords: SNES, nonlinear, set, function, norm, type 1702 1703 .seealso: SNESSetNormSchedule(), SNESComputeFunction(), VecNorm(), SNESSetFunction(), SNESSetInitialFunction(), SNESNormSchedule 1704 @*/ 1705 PetscErrorCode SNESGetNormSchedule(SNES snes, SNESNormSchedule *normschedule) 1706 { 1707 PetscFunctionBegin; 1708 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1709 *normschedule = snes->normschedule; 1710 PetscFunctionReturn(0); 1711 } 1712 1713 1714 #undef __FUNCT__ 1715 #define __FUNCT__ "SNESSetFunctionType" 1716 /*@C 1717 SNESSetFunctionType - Sets the SNESNormSchedule used in covergence and monitoring 1718 of the SNES method. 1719 1720 Logically Collective on SNES 1721 1722 Input Parameters: 1723 + snes - the SNES context 1724 - normschedule - the frequency of norm computation 1725 1726 Notes: 1727 Only certain SNES methods support certain SNESNormSchedules. Most require evaluation 1728 of the nonlinear function and the taking of its norm at every iteration to 1729 even ensure convergence at all. However, methods such as custom Gauss-Seidel methods 1730 (SNESGS) and the like do not require the norm of the function to be computed, and therfore 1731 may either be monitored for convergence or not. As these are often used as nonlinear 1732 preconditioners, monitoring the norm of their error is not a useful enterprise within 1733 their solution. 1734 1735 Level: developer 1736 1737 .keywords: SNES, nonlinear, set, function, norm, type 1738 1739 .seealso: SNESGetNormSchedule(), SNESComputeFunction(), VecNorm(), SNESSetFunction(), SNESSetInitialFunction(), SNESNormSchedule 1740 @*/ 1741 PetscErrorCode SNESSetFunctionType(SNES snes, SNESFunctionType type) 1742 { 1743 PetscFunctionBegin; 1744 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1745 snes->functype = type; 1746 PetscFunctionReturn(0); 1747 } 1748 1749 1750 #undef __FUNCT__ 1751 #define __FUNCT__ "SNESGetFunctionType" 1752 /*@C 1753 SNESGetFunctionType - Gets the SNESNormSchedule used in covergence and monitoring 1754 of the SNES method. 1755 1756 Logically Collective on SNES 1757 1758 Input Parameters: 1759 + snes - the SNES context 1760 - normschedule - the type of the norm used 1761 1762 Level: advanced 1763 1764 .keywords: SNES, nonlinear, set, function, norm, type 1765 1766 .seealso: SNESSetNormSchedule(), SNESComputeFunction(), VecNorm(), SNESSetFunction(), SNESSetInitialFunction(), SNESNormSchedule 1767 @*/ 1768 PetscErrorCode SNESGetFunctionType(SNES snes, SNESFunctionType *type) 1769 { 1770 PetscFunctionBegin; 1771 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1772 *type = snes->functype; 1773 PetscFunctionReturn(0); 1774 } 1775 1776 /*MC 1777 SNESGSFunction - function used to convey a Gauss-Seidel sweep on the nonlinear function 1778 1779 Synopsis: 1780 #include "petscsnes.h" 1781 $ SNESGSFunction(SNES snes,Vec x,Vec b,void *ctx); 1782 1783 + X - solution vector 1784 . B - RHS vector 1785 - ctx - optional user-defined Gauss-Seidel context 1786 1787 Level: intermediate 1788 1789 .seealso: SNESSetGS(), SNESGetGS() 1790 M*/ 1791 1792 #undef __FUNCT__ 1793 #define __FUNCT__ "SNESSetGS" 1794 /*@C 1795 SNESSetGS - Sets the user nonlinear Gauss-Seidel routine for 1796 use with composed nonlinear solvers. 1797 1798 Input Parameters: 1799 + snes - the SNES context 1800 . SNESGSFunction - function evaluation routine 1801 - ctx - [optional] user-defined context for private data for the 1802 smoother evaluation routine (may be NULL) 1803 1804 Notes: 1805 The GS routines are used by the composed nonlinear solver to generate 1806 a problem appropriate update to the solution, particularly FAS. 1807 1808 Level: intermediate 1809 1810 .keywords: SNES, nonlinear, set, Gauss-Seidel 1811 1812 .seealso: SNESGetFunction(), SNESComputeGS() 1813 @*/ 1814 PetscErrorCode SNESSetGS(SNES snes,PetscErrorCode (*SNESGSFunction)(SNES,Vec,Vec,void*),void *ctx) 1815 { 1816 PetscErrorCode ierr; 1817 DM dm; 1818 1819 PetscFunctionBegin; 1820 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1821 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 1822 ierr = DMSNESSetGS(dm,SNESGSFunction,ctx);CHKERRQ(ierr); 1823 PetscFunctionReturn(0); 1824 } 1825 1826 #undef __FUNCT__ 1827 #define __FUNCT__ "SNESPicardComputeFunction" 1828 PETSC_EXTERN PetscErrorCode SNESPicardComputeFunction(SNES snes,Vec x,Vec f,void *ctx) 1829 { 1830 PetscErrorCode ierr; 1831 DM dm; 1832 DMSNES sdm; 1833 1834 PetscFunctionBegin; 1835 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 1836 ierr = DMGetDMSNES(dm,&sdm);CHKERRQ(ierr); 1837 /* A(x)*x - b(x) */ 1838 if (sdm->ops->computepfunction) { 1839 ierr = (*sdm->ops->computepfunction)(snes,x,f,sdm->pctx);CHKERRQ(ierr); 1840 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE, "Must call SNESSetPicard() to provide Picard function."); 1841 1842 if (sdm->ops->computepjacobian) { 1843 ierr = (*sdm->ops->computepjacobian)(snes,x,&snes->jacobian,&snes->jacobian_pre,&snes->matstruct,sdm->pctx);CHKERRQ(ierr); 1844 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE, "Must call SNESSetPicard() to provide Picard matrix."); 1845 ierr = VecScale(f,-1.0);CHKERRQ(ierr); 1846 ierr = MatMultAdd(snes->jacobian,x,f,f);CHKERRQ(ierr); 1847 PetscFunctionReturn(0); 1848 } 1849 1850 #undef __FUNCT__ 1851 #define __FUNCT__ "SNESPicardComputeJacobian" 1852 PETSC_EXTERN PetscErrorCode SNESPicardComputeJacobian(SNES snes,Vec x1,Mat *J,Mat *B,MatStructure *flag,void *ctx) 1853 { 1854 PetscFunctionBegin; 1855 /* the jacobian matrix should be pre-filled in SNESPicardComputeFunction */ 1856 *flag = snes->matstruct; 1857 PetscFunctionReturn(0); 1858 } 1859 1860 #undef __FUNCT__ 1861 #define __FUNCT__ "SNESSetPicard" 1862 /*@C 1863 SNESSetPicard - Use SNES to solve the semilinear-system A(x) x = b(x) via a Picard type iteration (Picard linearization) 1864 1865 Logically Collective on SNES 1866 1867 Input Parameters: 1868 + snes - the SNES context 1869 . r - vector to store function value 1870 . SNESFunction - function evaluation routine 1871 . Amat - matrix with which A(x) x - b(x) is to be computed 1872 . Pmat - matrix from which preconditioner is computed (usually the same as Amat) 1873 . SNESJacobianFunction - function to compute matrix value 1874 - ctx - [optional] user-defined context for private data for the 1875 function evaluation routine (may be NULL) 1876 1877 Notes: 1878 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 1879 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. 1880 1881 One can call SNESSetPicard() or SNESSetFunction() (and possibly SNESSetJacobian()) but cannot call both 1882 1883 $ 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} 1884 $ Note that when an exact solver is used this corresponds to the "classic" Picard A(x^{n}) x^{n+1} = b(x^{n}) iteration. 1885 1886 Run with -snes_mf_operator to solve the system with Newton's method using A(x^{n}) to construct the preconditioner. 1887 1888 We implement the defect correction form of the Picard iteration because it converges much more generally when inexact linear solvers are used then 1889 the direct Picard iteration A(x^n) x^{n+1} = b(x^n) 1890 1891 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 1892 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 1893 different please contact us at petsc-dev@mcs.anl.gov and we'll have an entirely new argument :-). 1894 1895 Level: intermediate 1896 1897 .keywords: SNES, nonlinear, set, function 1898 1899 .seealso: SNESGetFunction(), SNESSetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESGetPicard(), SNESLineSearchPreCheckPicard() 1900 @*/ 1901 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) 1902 { 1903 PetscErrorCode ierr; 1904 DM dm; 1905 1906 PetscFunctionBegin; 1907 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1908 ierr = SNESGetDM(snes, &dm);CHKERRQ(ierr); 1909 ierr = DMSNESSetPicard(dm,SNESFunction,SNESJacobianFunction,ctx);CHKERRQ(ierr); 1910 ierr = SNESSetFunction(snes,r,SNESPicardComputeFunction,ctx);CHKERRQ(ierr); 1911 ierr = SNESSetJacobian(snes,Amat,Pmat,SNESPicardComputeJacobian,ctx);CHKERRQ(ierr); 1912 PetscFunctionReturn(0); 1913 } 1914 1915 #undef __FUNCT__ 1916 #define __FUNCT__ "SNESGetPicard" 1917 /*@C 1918 SNESGetPicard - Returns the context for the Picard iteration 1919 1920 Not Collective, but Vec is parallel if SNES is parallel. Collective if Vec is requested, but has not been created yet. 1921 1922 Input Parameter: 1923 . snes - the SNES context 1924 1925 Output Parameter: 1926 + r - the function (or NULL) 1927 . SNESFunction - the function (or NULL) 1928 . Amat - the matrix used to defined the operation A(x) x - b(x) (or NULL) 1929 . Pmat - the matrix from which the preconditioner will be constructed (or NULL) 1930 . SNESJacobianFunction - the function for matrix evaluation (or NULL) 1931 - ctx - the function context (or NULL) 1932 1933 Level: advanced 1934 1935 .keywords: SNES, nonlinear, get, function 1936 1937 .seealso: SNESSetPicard(), SNESGetFunction(), SNESGetJacobian(), SNESGetDM(), SNESFunction, SNESJacobianFunction 1938 @*/ 1939 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) 1940 { 1941 PetscErrorCode ierr; 1942 DM dm; 1943 1944 PetscFunctionBegin; 1945 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1946 ierr = SNESGetFunction(snes,r,NULL,NULL);CHKERRQ(ierr); 1947 ierr = SNESGetJacobian(snes,Amat,Pmat,NULL,NULL);CHKERRQ(ierr); 1948 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 1949 ierr = DMSNESGetPicard(dm,SNESFunction,SNESJacobianFunction,ctx);CHKERRQ(ierr); 1950 PetscFunctionReturn(0); 1951 } 1952 1953 #undef __FUNCT__ 1954 #define __FUNCT__ "SNESSetComputeInitialGuess" 1955 /*@C 1956 SNESSetComputeInitialGuess - Sets a routine used to compute an initial guess for the problem 1957 1958 Logically Collective on SNES 1959 1960 Input Parameters: 1961 + snes - the SNES context 1962 . func - function evaluation routine 1963 - ctx - [optional] user-defined context for private data for the 1964 function evaluation routine (may be NULL) 1965 1966 Calling sequence of func: 1967 $ func (SNES snes,Vec x,void *ctx); 1968 1969 . f - function vector 1970 - ctx - optional user-defined function context 1971 1972 Level: intermediate 1973 1974 .keywords: SNES, nonlinear, set, function 1975 1976 .seealso: SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian() 1977 @*/ 1978 PetscErrorCode SNESSetComputeInitialGuess(SNES snes,PetscErrorCode (*func)(SNES,Vec,void*),void *ctx) 1979 { 1980 PetscFunctionBegin; 1981 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1982 if (func) snes->ops->computeinitialguess = func; 1983 if (ctx) snes->initialguessP = ctx; 1984 PetscFunctionReturn(0); 1985 } 1986 1987 /* --------------------------------------------------------------- */ 1988 #undef __FUNCT__ 1989 #define __FUNCT__ "SNESGetRhs" 1990 /*@C 1991 SNESGetRhs - Gets the vector for solving F(x) = rhs. If rhs is not set 1992 it assumes a zero right hand side. 1993 1994 Logically Collective on SNES 1995 1996 Input Parameter: 1997 . snes - the SNES context 1998 1999 Output Parameter: 2000 . rhs - the right hand side vector or NULL if the right hand side vector is null 2001 2002 Level: intermediate 2003 2004 .keywords: SNES, nonlinear, get, function, right hand side 2005 2006 .seealso: SNESGetSolution(), SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESSetFunction() 2007 @*/ 2008 PetscErrorCode SNESGetRhs(SNES snes,Vec *rhs) 2009 { 2010 PetscFunctionBegin; 2011 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2012 PetscValidPointer(rhs,2); 2013 *rhs = snes->vec_rhs; 2014 PetscFunctionReturn(0); 2015 } 2016 2017 #undef __FUNCT__ 2018 #define __FUNCT__ "SNESComputeFunction" 2019 /*@ 2020 SNESComputeFunction - Calls the function that has been set with SNESSetFunction(). 2021 2022 Collective on SNES 2023 2024 Input Parameters: 2025 + snes - the SNES context 2026 - x - input vector 2027 2028 Output Parameter: 2029 . y - function vector, as set by SNESSetFunction() 2030 2031 Notes: 2032 SNESComputeFunction() is typically used within nonlinear solvers 2033 implementations, so most users would not generally call this routine 2034 themselves. 2035 2036 Level: developer 2037 2038 .keywords: SNES, nonlinear, compute, function 2039 2040 .seealso: SNESSetFunction(), SNESGetFunction() 2041 @*/ 2042 PetscErrorCode SNESComputeFunction(SNES snes,Vec x,Vec y) 2043 { 2044 PetscErrorCode ierr; 2045 DM dm; 2046 DMSNES sdm; 2047 2048 PetscFunctionBegin; 2049 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2050 PetscValidHeaderSpecific(x,VEC_CLASSID,2); 2051 PetscValidHeaderSpecific(y,VEC_CLASSID,3); 2052 PetscCheckSameComm(snes,1,x,2); 2053 PetscCheckSameComm(snes,1,y,3); 2054 ierr = VecValidValues(x,2,PETSC_TRUE);CHKERRQ(ierr); 2055 2056 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 2057 ierr = DMGetDMSNES(dm,&sdm);CHKERRQ(ierr); 2058 if (sdm->ops->computefunction) { 2059 ierr = PetscLogEventBegin(SNES_FunctionEval,snes,x,y,0);CHKERRQ(ierr); 2060 PetscStackPush("SNES user function"); 2061 ierr = (*sdm->ops->computefunction)(snes,x,y,sdm->functionctx);CHKERRQ(ierr); 2062 PetscStackPop; 2063 ierr = PetscLogEventEnd(SNES_FunctionEval,snes,x,y,0);CHKERRQ(ierr); 2064 } else if (snes->vec_rhs) { 2065 ierr = MatMult(snes->jacobian, x, y);CHKERRQ(ierr); 2066 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE, "Must call SNESSetFunction() or SNESSetDM() before SNESComputeFunction(), likely called from SNESSolve()."); 2067 if (snes->vec_rhs) { 2068 ierr = VecAXPY(y,-1.0,snes->vec_rhs);CHKERRQ(ierr); 2069 } 2070 snes->nfuncs++; 2071 ierr = VecValidValues(y,3,PETSC_FALSE);CHKERRQ(ierr); 2072 PetscFunctionReturn(0); 2073 } 2074 2075 #undef __FUNCT__ 2076 #define __FUNCT__ "SNESComputeGS" 2077 /*@ 2078 SNESComputeGS - Calls the Gauss-Seidel function that has been set with SNESSetGS(). 2079 2080 Collective on SNES 2081 2082 Input Parameters: 2083 + snes - the SNES context 2084 . x - input vector 2085 - b - rhs vector 2086 2087 Output Parameter: 2088 . x - new solution vector 2089 2090 Notes: 2091 SNESComputeGS() is typically used within composed nonlinear solver 2092 implementations, so most users would not generally call this routine 2093 themselves. 2094 2095 Level: developer 2096 2097 .keywords: SNES, nonlinear, compute, function 2098 2099 .seealso: SNESSetGS(), SNESComputeFunction() 2100 @*/ 2101 PetscErrorCode SNESComputeGS(SNES snes,Vec b,Vec x) 2102 { 2103 PetscErrorCode ierr; 2104 DM dm; 2105 DMSNES sdm; 2106 2107 PetscFunctionBegin; 2108 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2109 PetscValidHeaderSpecific(x,VEC_CLASSID,2); 2110 if (b) PetscValidHeaderSpecific(b,VEC_CLASSID,3); 2111 PetscCheckSameComm(snes,1,x,2); 2112 if (b) PetscCheckSameComm(snes,1,b,3); 2113 if (b) {ierr = VecValidValues(b,2,PETSC_TRUE);CHKERRQ(ierr);} 2114 ierr = PetscLogEventBegin(SNES_GSEval,snes,x,b,0);CHKERRQ(ierr); 2115 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 2116 ierr = DMGetDMSNES(dm,&sdm);CHKERRQ(ierr); 2117 if (sdm->ops->computegs) { 2118 PetscStackPush("SNES user GS"); 2119 ierr = (*sdm->ops->computegs)(snes,x,b,sdm->gsctx);CHKERRQ(ierr); 2120 PetscStackPop; 2121 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE, "Must call SNESSetGS() before SNESComputeGS(), likely called from SNESSolve()."); 2122 ierr = PetscLogEventEnd(SNES_GSEval,snes,x,b,0);CHKERRQ(ierr); 2123 ierr = VecValidValues(x,3,PETSC_FALSE);CHKERRQ(ierr); 2124 PetscFunctionReturn(0); 2125 } 2126 2127 #undef __FUNCT__ 2128 #define __FUNCT__ "SNESComputeJacobian" 2129 /*@ 2130 SNESComputeJacobian - Computes the Jacobian matrix that has been set with SNESSetJacobian(). 2131 2132 Collective on SNES and Mat 2133 2134 Input Parameters: 2135 + snes - the SNES context 2136 - x - input vector 2137 2138 Output Parameters: 2139 + A - Jacobian matrix 2140 . B - optional preconditioning matrix 2141 - flag - flag indicating matrix structure (one of, SAME_NONZERO_PATTERN,DIFFERENT_NONZERO_PATTERN,SAME_PRECONDITIONER) 2142 2143 Options Database Keys: 2144 + -snes_lag_preconditioner <lag> 2145 . -snes_lag_jacobian <lag> 2146 . -snes_compare_explicit - Compare the computed Jacobian to the finite difference Jacobian and output the differences 2147 . -snes_compare_explicit_draw - Compare the computed Jacobian to the finite difference Jacobian and draw the result 2148 . -snes_compare_explicit_contour - Compare the computed Jacobian to the finite difference Jacobian and draw a contour plot with the result 2149 . -snes_compare_operator - Make the comparison options above use the operator instead of the preconditioning matrix 2150 . -snes_compare_coloring - Compute the finite differece Jacobian using coloring and display norms of difference 2151 . -snes_compare_coloring_display - Compute the finite differece Jacobian using coloring and display verbose differences 2152 . -snes_compare_coloring_threshold - Display only those matrix entries that differ by more than a given threshold 2153 . -snes_compare_coloring_threshold_atol - Absolute tolerance for difference in matrix entries to be displayed by -snes_compare_coloring_threshold 2154 . -snes_compare_coloring_threshold_rtol - Relative tolerance for difference in matrix entries to be displayed by -snes_compare_coloring_threshold 2155 . -snes_compare_coloring_draw - Compute the finite differece Jacobian using coloring and draw differences 2156 - -snes_compare_coloring_draw_contour - Compute the finite differece Jacobian using coloring and show contours of matrices and differences 2157 2158 2159 Notes: 2160 Most users should not need to explicitly call this routine, as it 2161 is used internally within the nonlinear solvers. 2162 2163 See KSPSetOperators() for important information about setting the 2164 flag parameter. 2165 2166 Level: developer 2167 2168 .keywords: SNES, compute, Jacobian, matrix 2169 2170 .seealso: SNESSetJacobian(), KSPSetOperators(), MatStructure, SNESSetLagPreconditioner(), SNESSetLagJacobian() 2171 @*/ 2172 PetscErrorCode SNESComputeJacobian(SNES snes,Vec X,Mat *A,Mat *B,MatStructure *flg) 2173 { 2174 PetscErrorCode ierr; 2175 PetscBool flag; 2176 DM dm; 2177 DMSNES sdm; 2178 2179 PetscFunctionBegin; 2180 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2181 PetscValidHeaderSpecific(X,VEC_CLASSID,2); 2182 PetscValidPointer(flg,5); 2183 PetscCheckSameComm(snes,1,X,2); 2184 ierr = VecValidValues(X,2,PETSC_TRUE);CHKERRQ(ierr); 2185 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 2186 ierr = DMGetDMSNES(dm,&sdm);CHKERRQ(ierr); 2187 2188 if (!sdm->ops->computejacobian) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_USER,"Must call SNESSetJacobian(), DMSNESSetJacobian(), DMDASNESSetJacobianLocal(), etc"); 2189 2190 /* make sure that MatAssemblyBegin/End() is called on A matrix if it is matrix free */ 2191 2192 if (snes->lagjacobian == -2) { 2193 snes->lagjacobian = -1; 2194 2195 ierr = PetscInfo(snes,"Recomputing Jacobian/preconditioner because lag is -2 (means compute Jacobian, but then never again) \n");CHKERRQ(ierr); 2196 } else if (snes->lagjacobian == -1) { 2197 *flg = SAME_PRECONDITIONER; 2198 ierr = PetscInfo(snes,"Reusing Jacobian/preconditioner because lag is -1\n");CHKERRQ(ierr); 2199 ierr = PetscObjectTypeCompare((PetscObject)*A,MATMFFD,&flag);CHKERRQ(ierr); 2200 if (flag) { 2201 ierr = MatAssemblyBegin(*A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2202 ierr = MatAssemblyEnd(*A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2203 } 2204 PetscFunctionReturn(0); 2205 } else if (snes->lagjacobian > 1 && (snes->iter + snes->jac_iter) % snes->lagjacobian) { 2206 *flg = SAME_PRECONDITIONER; 2207 ierr = PetscInfo2(snes,"Reusing Jacobian/preconditioner because lag is %D and SNES iteration is %D\n",snes->lagjacobian,snes->iter);CHKERRQ(ierr); 2208 ierr = PetscObjectTypeCompare((PetscObject)*A,MATMFFD,&flag);CHKERRQ(ierr); 2209 if (flag) { 2210 ierr = MatAssemblyBegin(*A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2211 ierr = MatAssemblyEnd(*A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2212 } 2213 PetscFunctionReturn(0); 2214 } 2215 if (snes->pc && snes->pcside == PC_LEFT) { 2216 ierr = MatAssemblyBegin(*A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2217 ierr = MatAssemblyEnd(*A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2218 PetscFunctionReturn(0); 2219 } 2220 2221 *flg = DIFFERENT_NONZERO_PATTERN; 2222 ierr = PetscLogEventBegin(SNES_JacobianEval,snes,X,*A,*B);CHKERRQ(ierr); 2223 2224 PetscStackPush("SNES user Jacobian function"); 2225 ierr = (*sdm->ops->computejacobian)(snes,X,A,B,flg,sdm->jacobianctx);CHKERRQ(ierr); 2226 PetscStackPop; 2227 2228 ierr = PetscLogEventEnd(SNES_JacobianEval,snes,X,*A,*B);CHKERRQ(ierr); 2229 2230 if (snes->lagpreconditioner == -2) { 2231 ierr = PetscInfo(snes,"Rebuilding preconditioner exactly once since lag is -2\n");CHKERRQ(ierr); 2232 2233 snes->lagpreconditioner = -1; 2234 } else if (snes->lagpreconditioner == -1) { 2235 *flg = SAME_PRECONDITIONER; 2236 ierr = PetscInfo(snes,"Reusing preconditioner because lag is -1\n");CHKERRQ(ierr); 2237 } else if (snes->lagpreconditioner > 1 && (snes->iter + snes->pre_iter) % snes->lagpreconditioner) { 2238 *flg = SAME_PRECONDITIONER; 2239 ierr = PetscInfo2(snes,"Reusing preconditioner because lag is %D and SNES iteration is %D\n",snes->lagpreconditioner,snes->iter);CHKERRQ(ierr); 2240 } 2241 2242 /* make sure user returned a correct Jacobian and preconditioner */ 2243 /* PetscValidHeaderSpecific(*A,MAT_CLASSID,3); 2244 PetscValidHeaderSpecific(*B,MAT_CLASSID,4); */ 2245 { 2246 PetscBool flag = PETSC_FALSE,flag_draw = PETSC_FALSE,flag_contour = PETSC_FALSE,flag_operator = PETSC_FALSE; 2247 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_explicit",&flag,NULL);CHKERRQ(ierr); 2248 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_explicit_draw",&flag_draw,NULL);CHKERRQ(ierr); 2249 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_explicit_draw_contour",&flag_contour,NULL);CHKERRQ(ierr); 2250 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_operator",&flag_operator,NULL);CHKERRQ(ierr); 2251 if (flag || flag_draw || flag_contour) { 2252 Mat Bexp_mine = NULL,Bexp,FDexp; 2253 MatStructure mstruct; 2254 PetscViewer vdraw,vstdout; 2255 PetscBool flg; 2256 if (flag_operator) { 2257 ierr = MatComputeExplicitOperator(*A,&Bexp_mine);CHKERRQ(ierr); 2258 Bexp = Bexp_mine; 2259 } else { 2260 /* See if the preconditioning matrix can be viewed and added directly */ 2261 ierr = PetscObjectTypeCompareAny((PetscObject)*B,&flg,MATSEQAIJ,MATMPIAIJ,MATSEQDENSE,MATMPIDENSE,MATSEQBAIJ,MATMPIBAIJ,MATSEQSBAIJ,MATMPIBAIJ,"");CHKERRQ(ierr); 2262 if (flg) Bexp = *B; 2263 else { 2264 /* If the "preconditioning" matrix is itself MATSHELL or some other type without direct support */ 2265 ierr = MatComputeExplicitOperator(*B,&Bexp_mine);CHKERRQ(ierr); 2266 Bexp = Bexp_mine; 2267 } 2268 } 2269 ierr = MatConvert(Bexp,MATSAME,MAT_INITIAL_MATRIX,&FDexp);CHKERRQ(ierr); 2270 ierr = SNESComputeJacobianDefault(snes,X,&FDexp,&FDexp,&mstruct,NULL);CHKERRQ(ierr); 2271 ierr = PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)snes),&vstdout);CHKERRQ(ierr); 2272 if (flag_draw || flag_contour) { 2273 ierr = PetscViewerDrawOpen(PetscObjectComm((PetscObject)snes),0,"Explicit Jacobians",PETSC_DECIDE,PETSC_DECIDE,300,300,&vdraw);CHKERRQ(ierr); 2274 if (flag_contour) {ierr = PetscViewerPushFormat(vdraw,PETSC_VIEWER_DRAW_CONTOUR);CHKERRQ(ierr);} 2275 } else vdraw = NULL; 2276 ierr = PetscViewerASCIIPrintf(vstdout,"Explicit %s\n",flag_operator ? "Jacobian" : "preconditioning Jacobian");CHKERRQ(ierr); 2277 if (flag) {ierr = MatView(Bexp,vstdout);CHKERRQ(ierr);} 2278 if (vdraw) {ierr = MatView(Bexp,vdraw);CHKERRQ(ierr);} 2279 ierr = PetscViewerASCIIPrintf(vstdout,"Finite difference Jacobian\n");CHKERRQ(ierr); 2280 if (flag) {ierr = MatView(FDexp,vstdout);CHKERRQ(ierr);} 2281 if (vdraw) {ierr = MatView(FDexp,vdraw);CHKERRQ(ierr);} 2282 ierr = MatAYPX(FDexp,-1.0,Bexp,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 2283 ierr = PetscViewerASCIIPrintf(vstdout,"User-provided matrix minus finite difference Jacobian\n");CHKERRQ(ierr); 2284 if (flag) {ierr = MatView(FDexp,vstdout);CHKERRQ(ierr);} 2285 if (vdraw) { /* Always use contour for the difference */ 2286 ierr = PetscViewerPushFormat(vdraw,PETSC_VIEWER_DRAW_CONTOUR);CHKERRQ(ierr); 2287 ierr = MatView(FDexp,vdraw);CHKERRQ(ierr); 2288 ierr = PetscViewerPopFormat(vdraw);CHKERRQ(ierr); 2289 } 2290 if (flag_contour) {ierr = PetscViewerPopFormat(vdraw);CHKERRQ(ierr);} 2291 ierr = PetscViewerDestroy(&vdraw);CHKERRQ(ierr); 2292 ierr = MatDestroy(&Bexp_mine);CHKERRQ(ierr); 2293 ierr = MatDestroy(&FDexp);CHKERRQ(ierr); 2294 } 2295 } 2296 { 2297 PetscBool flag = PETSC_FALSE,flag_display = PETSC_FALSE,flag_draw = PETSC_FALSE,flag_contour = PETSC_FALSE,flag_threshold = PETSC_FALSE; 2298 PetscReal threshold_atol = PETSC_SQRT_MACHINE_EPSILON,threshold_rtol = 10*PETSC_SQRT_MACHINE_EPSILON; 2299 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_coloring",&flag,NULL);CHKERRQ(ierr); 2300 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_coloring_display",&flag_display,NULL);CHKERRQ(ierr); 2301 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_coloring_draw",&flag_draw,NULL);CHKERRQ(ierr); 2302 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_coloring_draw_contour",&flag_contour,NULL);CHKERRQ(ierr); 2303 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_coloring_threshold",&flag_threshold,NULL);CHKERRQ(ierr); 2304 ierr = PetscOptionsGetReal(((PetscObject)snes)->prefix,"-snes_compare_coloring_threshold_rtol",&threshold_rtol,NULL);CHKERRQ(ierr); 2305 ierr = PetscOptionsGetReal(((PetscObject)snes)->prefix,"-snes_compare_coloring_threshold_atol",&threshold_atol,NULL);CHKERRQ(ierr); 2306 if (flag || flag_display || flag_draw || flag_contour || flag_threshold) { 2307 Mat Bfd; 2308 MatStructure mstruct; 2309 PetscViewer vdraw,vstdout; 2310 MatColoring coloring; 2311 ISColoring iscoloring; 2312 MatFDColoring matfdcoloring; 2313 PetscErrorCode (*func)(SNES,Vec,Vec,void*); 2314 void *funcctx; 2315 PetscReal norm1,norm2,normmax; 2316 2317 ierr = MatDuplicate(*B,MAT_DO_NOT_COPY_VALUES,&Bfd);CHKERRQ(ierr); 2318 ierr = MatColoringCreate(Bfd,&coloring);CHKERRQ(ierr); 2319 ierr = MatColoringSetType(coloring,MATCOLORINGSL);CHKERRQ(ierr); 2320 ierr = MatColoringSetFromOptions(coloring);CHKERRQ(ierr); 2321 ierr = MatColoringApply(coloring,&iscoloring);CHKERRQ(ierr); 2322 ierr = MatColoringDestroy(&coloring);CHKERRQ(ierr); 2323 ierr = MatFDColoringCreate(Bfd,iscoloring,&matfdcoloring);CHKERRQ(ierr); 2324 ierr = ISColoringDestroy(&iscoloring);CHKERRQ(ierr); 2325 2326 /* This method of getting the function is currently unreliable since it doesn't work for DM local functions. */ 2327 ierr = SNESGetFunction(snes,NULL,&func,&funcctx);CHKERRQ(ierr); 2328 ierr = MatFDColoringSetFunction(matfdcoloring,(PetscErrorCode (*)(void))func,funcctx);CHKERRQ(ierr); 2329 ierr = PetscObjectSetOptionsPrefix((PetscObject)matfdcoloring,((PetscObject)snes)->prefix);CHKERRQ(ierr); 2330 ierr = PetscObjectAppendOptionsPrefix((PetscObject)matfdcoloring,"coloring_");CHKERRQ(ierr); 2331 ierr = MatFDColoringSetFromOptions(matfdcoloring);CHKERRQ(ierr); 2332 ierr = MatFDColoringApply(Bfd,matfdcoloring,X,&mstruct,snes);CHKERRQ(ierr); 2333 ierr = MatFDColoringDestroy(&matfdcoloring);CHKERRQ(ierr); 2334 2335 ierr = PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)snes),&vstdout);CHKERRQ(ierr); 2336 if (flag_draw || flag_contour) { 2337 ierr = PetscViewerDrawOpen(PetscObjectComm((PetscObject)snes),0,"Colored Jacobians",PETSC_DECIDE,PETSC_DECIDE,300,300,&vdraw);CHKERRQ(ierr); 2338 if (flag_contour) {ierr = PetscViewerPushFormat(vdraw,PETSC_VIEWER_DRAW_CONTOUR);CHKERRQ(ierr);} 2339 } else vdraw = NULL; 2340 ierr = PetscViewerASCIIPrintf(vstdout,"Explicit preconditioning Jacobian\n");CHKERRQ(ierr); 2341 if (flag_display) {ierr = MatView(*B,vstdout);CHKERRQ(ierr);} 2342 if (vdraw) {ierr = MatView(*B,vdraw);CHKERRQ(ierr);} 2343 ierr = PetscViewerASCIIPrintf(vstdout,"Colored Finite difference Jacobian\n");CHKERRQ(ierr); 2344 if (flag_display) {ierr = MatView(Bfd,vstdout);CHKERRQ(ierr);} 2345 if (vdraw) {ierr = MatView(Bfd,vdraw);CHKERRQ(ierr);} 2346 ierr = MatAYPX(Bfd,-1.0,*B,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 2347 ierr = MatNorm(Bfd,NORM_1,&norm1);CHKERRQ(ierr); 2348 ierr = MatNorm(Bfd,NORM_FROBENIUS,&norm2);CHKERRQ(ierr); 2349 ierr = MatNorm(Bfd,NORM_MAX,&normmax);CHKERRQ(ierr); 2350 ierr = PetscViewerASCIIPrintf(vstdout,"User-provided matrix minus finite difference Jacobian, norm1=%G normFrob=%G normmax=%G\n",norm1,norm2,normmax);CHKERRQ(ierr); 2351 if (flag_display) {ierr = MatView(Bfd,vstdout);CHKERRQ(ierr);} 2352 if (vdraw) { /* Always use contour for the difference */ 2353 ierr = PetscViewerPushFormat(vdraw,PETSC_VIEWER_DRAW_CONTOUR);CHKERRQ(ierr); 2354 ierr = MatView(Bfd,vdraw);CHKERRQ(ierr); 2355 ierr = PetscViewerPopFormat(vdraw);CHKERRQ(ierr); 2356 } 2357 if (flag_contour) {ierr = PetscViewerPopFormat(vdraw);CHKERRQ(ierr);} 2358 2359 if (flag_threshold) { 2360 PetscInt bs,rstart,rend,i; 2361 ierr = MatGetBlockSize(*B,&bs);CHKERRQ(ierr); 2362 ierr = MatGetOwnershipRange(*B,&rstart,&rend);CHKERRQ(ierr); 2363 for (i=rstart; i<rend; i++) { 2364 const PetscScalar *ba,*ca; 2365 const PetscInt *bj,*cj; 2366 PetscInt bn,cn,j,maxentrycol = -1,maxdiffcol = -1,maxrdiffcol = -1; 2367 PetscReal maxentry = 0,maxdiff = 0,maxrdiff = 0; 2368 ierr = MatGetRow(*B,i,&bn,&bj,&ba);CHKERRQ(ierr); 2369 ierr = MatGetRow(Bfd,i,&cn,&cj,&ca);CHKERRQ(ierr); 2370 if (bn != cn) SETERRQ(((PetscObject)*A)->comm,PETSC_ERR_PLIB,"Unexpected different nonzero pattern in -snes_compare_coloring_threshold"); 2371 for (j=0; j<bn; j++) { 2372 PetscReal rdiff = PetscAbsScalar(ca[j]) / (threshold_atol + threshold_rtol*PetscAbsScalar(ba[j])); 2373 if (PetscAbsScalar(ba[j]) > PetscAbs(maxentry)) { 2374 maxentrycol = bj[j]; 2375 maxentry = PetscRealPart(ba[j]); 2376 } 2377 if (PetscAbsScalar(ca[j]) > PetscAbs(maxdiff)) { 2378 maxdiffcol = bj[j]; 2379 maxdiff = PetscRealPart(ca[j]); 2380 } 2381 if (rdiff > maxrdiff) { 2382 maxrdiffcol = bj[j]; 2383 maxrdiff = rdiff; 2384 } 2385 } 2386 if (maxrdiff > 1) { 2387 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); 2388 for (j=0; j<bn; j++) { 2389 PetscReal rdiff; 2390 rdiff = PetscAbsScalar(ca[j]) / (threshold_atol + threshold_rtol*PetscAbsScalar(ba[j])); 2391 if (rdiff > 1) { 2392 ierr = PetscViewerASCIIPrintf(vstdout," (%D,%G:%G)",bj[j],PetscRealPart(ba[j]),PetscRealPart(ca[j]));CHKERRQ(ierr); 2393 } 2394 } 2395 ierr = PetscViewerASCIIPrintf(vstdout,"\n",i,maxentry,maxdiff,maxrdiff);CHKERRQ(ierr); 2396 } 2397 ierr = MatRestoreRow(*B,i,&bn,&bj,&ba);CHKERRQ(ierr); 2398 ierr = MatRestoreRow(Bfd,i,&cn,&cj,&ca);CHKERRQ(ierr); 2399 } 2400 } 2401 ierr = PetscViewerDestroy(&vdraw);CHKERRQ(ierr); 2402 ierr = MatDestroy(&Bfd);CHKERRQ(ierr); 2403 } 2404 } 2405 PetscFunctionReturn(0); 2406 } 2407 2408 /*MC 2409 SNESJacobianFunction - function used to convey the nonlinear Jacobian of the function to be solved by SNES 2410 2411 Synopsis: 2412 #include "petscsnes.h" 2413 $ SNESJacobianFunction(SNES snes,Vec x,Mat *Amat,Mat *Pmat,int *flag,void *ctx); 2414 2415 + x - input vector 2416 . Amat - the matrix that defines the (approximate) Jacobian 2417 . Pmat - the matrix to be used in constructing the preconditioner, usually the same as Amat. 2418 . flag - flag indicating information about the preconditioner matrix 2419 structure (same as flag in KSPSetOperators()), one of SAME_NONZERO_PATTERN,DIFFERENT_NONZERO_PATTERN,SAME_PRECONDITIONER 2420 - ctx - [optional] user-defined Jacobian context 2421 2422 Level: intermediate 2423 2424 .seealso: SNESSetFunction(), SNESGetFunction(), SNESSetJacobian(), SNESGetJacobian() 2425 M*/ 2426 2427 #undef __FUNCT__ 2428 #define __FUNCT__ "SNESSetJacobian" 2429 /*@C 2430 SNESSetJacobian - Sets the function to compute Jacobian as well as the 2431 location to store the matrix. 2432 2433 Logically Collective on SNES and Mat 2434 2435 Input Parameters: 2436 + snes - the SNES context 2437 . Amat - the matrix that defines the (approximate) Jacobian 2438 . Pmat - the matrix to be used in constructing the preconditioner, usually the same as Amat. 2439 . SNESJacobianFunction - Jacobian evaluation routine (if NULL then SNES retains any previously set value) 2440 - ctx - [optional] user-defined context for private data for the 2441 Jacobian evaluation routine (may be NULL) (if NULL then SNES retains any previously set value) 2442 2443 Notes: 2444 See KSPSetOperators() for important information about setting the flag 2445 output parameter in the routine func(). Be sure to read this information! 2446 2447 The routine func() takes Mat * as the matrix arguments rather than Mat. 2448 This allows the Jacobian evaluation routine to replace A and/or B with a 2449 completely new new matrix structure (not just different matrix elements) 2450 when appropriate, for instance, if the nonzero structure is changing 2451 throughout the global iterations. 2452 2453 If the Amat matrix and Pmat matrix are different you must call MatAssemblyBegin/End() on 2454 each matrix. 2455 2456 If using SNESComputeJacobianDefaultColor() to assemble a Jacobian, the ctx argument 2457 must be a MatFDColoring. 2458 2459 Other defect-correction schemes can be used by computing a different matrix in place of the Jacobian. One common 2460 example is to use the "Picard linearization" which only differentiates through the highest order parts of each term. 2461 2462 Level: beginner 2463 2464 .keywords: SNES, nonlinear, set, Jacobian, matrix 2465 2466 .seealso: KSPSetOperators(), SNESSetFunction(), MatMFFDComputeJacobian(), SNESComputeJacobianDefaultColor(), MatStructure, SNESJacobianFunction, SNESSetPicard() 2467 @*/ 2468 PetscErrorCode SNESSetJacobian(SNES snes,Mat Amat,Mat Pmat,PetscErrorCode (*SNESJacobianFunction)(SNES,Vec,Mat*,Mat*,MatStructure*,void*),void *ctx) 2469 { 2470 PetscErrorCode ierr; 2471 DM dm; 2472 2473 PetscFunctionBegin; 2474 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2475 if (Amat) PetscValidHeaderSpecific(Amat,MAT_CLASSID,2); 2476 if (Pmat) PetscValidHeaderSpecific(Pmat,MAT_CLASSID,3); 2477 if (Amat) PetscCheckSameComm(snes,1,Amat,2); 2478 if (Pmat) PetscCheckSameComm(snes,1,Pmat,3); 2479 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 2480 ierr = DMSNESSetJacobian(dm,SNESJacobianFunction,ctx);CHKERRQ(ierr); 2481 if (Amat) { 2482 ierr = PetscObjectReference((PetscObject)Amat);CHKERRQ(ierr); 2483 ierr = MatDestroy(&snes->jacobian);CHKERRQ(ierr); 2484 2485 snes->jacobian = Amat; 2486 } 2487 if (Pmat) { 2488 ierr = PetscObjectReference((PetscObject)Pmat);CHKERRQ(ierr); 2489 ierr = MatDestroy(&snes->jacobian_pre);CHKERRQ(ierr); 2490 2491 snes->jacobian_pre = Pmat; 2492 } 2493 PetscFunctionReturn(0); 2494 } 2495 2496 #undef __FUNCT__ 2497 #define __FUNCT__ "SNESGetJacobian" 2498 /*@C 2499 SNESGetJacobian - Returns the Jacobian matrix and optionally the user 2500 provided context for evaluating the Jacobian. 2501 2502 Not Collective, but Mat object will be parallel if SNES object is 2503 2504 Input Parameter: 2505 . snes - the nonlinear solver context 2506 2507 Output Parameters: 2508 + Amat - location to stash (approximate) Jacobian matrix (or NULL) 2509 . Pmat - location to stash matrix used to compute the preconditioner (or NULL) 2510 . SNESJacobianFunction - location to put Jacobian function (or NULL) 2511 - ctx - location to stash Jacobian ctx (or NULL) 2512 2513 Level: advanced 2514 2515 .seealso: SNESSetJacobian(), SNESComputeJacobian() 2516 @*/ 2517 PetscErrorCode SNESGetJacobian(SNES snes,Mat *Amat,Mat *Pmat,PetscErrorCode (**SNESJacobianFunction)(SNES,Vec,Mat*,Mat*,MatStructure*,void*),void **ctx) 2518 { 2519 PetscErrorCode ierr; 2520 DM dm; 2521 DMSNES sdm; 2522 2523 PetscFunctionBegin; 2524 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2525 if (Amat) *Amat = snes->jacobian; 2526 if (Pmat) *Pmat = snes->jacobian_pre; 2527 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 2528 ierr = DMGetDMSNES(dm,&sdm);CHKERRQ(ierr); 2529 if (SNESJacobianFunction) *SNESJacobianFunction = sdm->ops->computejacobian; 2530 if (ctx) *ctx = sdm->jacobianctx; 2531 PetscFunctionReturn(0); 2532 } 2533 2534 #undef __FUNCT__ 2535 #define __FUNCT__ "SNESSetUp" 2536 /*@ 2537 SNESSetUp - Sets up the internal data structures for the later use 2538 of a nonlinear solver. 2539 2540 Collective on SNES 2541 2542 Input Parameters: 2543 . snes - the SNES context 2544 2545 Notes: 2546 For basic use of the SNES solvers the user need not explicitly call 2547 SNESSetUp(), since these actions will automatically occur during 2548 the call to SNESSolve(). However, if one wishes to control this 2549 phase separately, SNESSetUp() should be called after SNESCreate() 2550 and optional routines of the form SNESSetXXX(), but before SNESSolve(). 2551 2552 Level: advanced 2553 2554 .keywords: SNES, nonlinear, setup 2555 2556 .seealso: SNESCreate(), SNESSolve(), SNESDestroy() 2557 @*/ 2558 PetscErrorCode SNESSetUp(SNES snes) 2559 { 2560 PetscErrorCode ierr; 2561 DM dm; 2562 DMSNES sdm; 2563 SNESLineSearch linesearch, pclinesearch; 2564 void *lsprectx,*lspostctx; 2565 PetscErrorCode (*precheck)(SNESLineSearch,Vec,Vec,PetscBool*,void*); 2566 PetscErrorCode (*postcheck)(SNESLineSearch,Vec,Vec,Vec,PetscBool*,PetscBool*,void*); 2567 PetscErrorCode (*func)(SNES,Vec,Vec,void*); 2568 Vec f,fpc; 2569 void *funcctx; 2570 PetscErrorCode (*jac)(SNES,Vec,Mat*,Mat*,MatStructure*,void*); 2571 void *jacctx,*appctx; 2572 2573 PetscFunctionBegin; 2574 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2575 if (snes->setupcalled) PetscFunctionReturn(0); 2576 2577 if (!((PetscObject)snes)->type_name) { 2578 ierr = SNESSetType(snes,SNESNEWTONLS);CHKERRQ(ierr); 2579 } 2580 2581 ierr = SNESGetFunction(snes,&snes->vec_func,NULL,NULL);CHKERRQ(ierr); 2582 if (snes->vec_func == snes->vec_sol) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_IDN,"Solution vector cannot be function vector"); 2583 if (snes->vec_rhs == snes->vec_sol) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_IDN,"Solution vector cannot be right hand side vector"); 2584 2585 if (!snes->vec_sol_update /* && snes->vec_sol */) { 2586 ierr = VecDuplicate(snes->vec_sol,&snes->vec_sol_update);CHKERRQ(ierr); 2587 ierr = PetscLogObjectParent((PetscObject)snes,(PetscObject)snes->vec_sol_update);CHKERRQ(ierr); 2588 } 2589 2590 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 2591 ierr = DMShellSetGlobalVector(snes->dm,snes->vec_sol);CHKERRQ(ierr); 2592 ierr = DMGetDMSNES(dm,&sdm);CHKERRQ(ierr); 2593 if (!sdm->ops->computefunction) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"Function never provided to SNES object"); 2594 if (!sdm->ops->computejacobian) { 2595 ierr = DMSNESSetJacobian(dm,SNESComputeJacobianDefaultColor,NULL);CHKERRQ(ierr); 2596 } 2597 if (!snes->vec_func) { 2598 ierr = DMCreateGlobalVector(dm,&snes->vec_func);CHKERRQ(ierr); 2599 } 2600 2601 if (!snes->ksp) { 2602 ierr = SNESGetKSP(snes, &snes->ksp);CHKERRQ(ierr); 2603 } 2604 2605 if (!snes->linesearch) { 2606 ierr = SNESGetLineSearch(snes, &snes->linesearch);CHKERRQ(ierr); 2607 } 2608 ierr = SNESLineSearchSetFunction(snes->linesearch,SNESComputeFunction);CHKERRQ(ierr); 2609 2610 if (snes->pc && (snes->pcside == PC_LEFT)) { 2611 snes->mf = PETSC_TRUE; 2612 snes->mf_operator = PETSC_FALSE; 2613 } 2614 2615 if (snes->mf) { 2616 ierr = SNESSetUpMatrixFree_Private(snes, snes->mf_operator, snes->mf_version);CHKERRQ(ierr); 2617 } 2618 2619 if (snes->ops->usercompute && !snes->user) { 2620 ierr = (*snes->ops->usercompute)(snes,(void**)&snes->user);CHKERRQ(ierr); 2621 } 2622 2623 if (snes->pc) { 2624 /* copy the DM over */ 2625 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 2626 ierr = SNESSetDM(snes->pc,dm);CHKERRQ(ierr); 2627 2628 ierr = SNESGetFunction(snes,&f,&func,&funcctx);CHKERRQ(ierr); 2629 ierr = VecDuplicate(f,&fpc);CHKERRQ(ierr); 2630 ierr = SNESSetFunction(snes->pc,fpc,func,funcctx);CHKERRQ(ierr); 2631 ierr = SNESGetJacobian(snes,NULL,NULL,&jac,&jacctx);CHKERRQ(ierr); 2632 ierr = SNESSetJacobian(snes->pc,NULL,NULL,jac,jacctx);CHKERRQ(ierr); 2633 ierr = SNESGetApplicationContext(snes,&appctx);CHKERRQ(ierr); 2634 ierr = SNESSetApplicationContext(snes->pc,appctx);CHKERRQ(ierr); 2635 ierr = VecDestroy(&fpc);CHKERRQ(ierr); 2636 2637 /* copy the function pointers over */ 2638 ierr = PetscObjectCopyFortranFunctionPointers((PetscObject)snes,(PetscObject)snes->pc);CHKERRQ(ierr); 2639 2640 /* default to 1 iteration */ 2641 ierr = SNESSetTolerances(snes->pc,0.0,0.0,0.0,1,snes->pc->max_funcs);CHKERRQ(ierr); 2642 if (snes->pcside==PC_RIGHT) { 2643 ierr = SNESSetNormSchedule(snes->pc,SNES_NORM_FINAL_ONLY);CHKERRQ(ierr); 2644 } else { 2645 ierr = SNESSetNormSchedule(snes->pc,SNES_NORM_NONE);CHKERRQ(ierr); 2646 } 2647 ierr = SNESSetFromOptions(snes->pc);CHKERRQ(ierr); 2648 2649 /* copy the line search context over */ 2650 ierr = SNESGetLineSearch(snes,&linesearch);CHKERRQ(ierr); 2651 ierr = SNESGetLineSearch(snes->pc,&pclinesearch);CHKERRQ(ierr); 2652 ierr = SNESLineSearchGetPreCheck(linesearch,&precheck,&lsprectx);CHKERRQ(ierr); 2653 ierr = SNESLineSearchGetPostCheck(linesearch,&postcheck,&lspostctx);CHKERRQ(ierr); 2654 ierr = SNESLineSearchSetPreCheck(pclinesearch,precheck,lsprectx);CHKERRQ(ierr); 2655 ierr = SNESLineSearchSetPostCheck(pclinesearch,postcheck,lspostctx);CHKERRQ(ierr); 2656 ierr = PetscObjectCopyFortranFunctionPointers((PetscObject)linesearch, (PetscObject)pclinesearch);CHKERRQ(ierr); 2657 } 2658 2659 snes->jac_iter = 0; 2660 snes->pre_iter = 0; 2661 2662 if (snes->ops->setup) { 2663 ierr = (*snes->ops->setup)(snes);CHKERRQ(ierr); 2664 } 2665 2666 if (snes->pc && (snes->pcside == PC_LEFT)) { 2667 if (snes->functype == SNES_FUNCTION_PRECONDITIONED) { 2668 ierr = SNESGetLineSearch(snes,&linesearch);CHKERRQ(ierr); 2669 ierr = SNESLineSearchSetFunction(linesearch,SNESComputeFunctionDefaultPC);CHKERRQ(ierr); 2670 } 2671 } 2672 2673 snes->setupcalled = PETSC_TRUE; 2674 PetscFunctionReturn(0); 2675 } 2676 2677 #undef __FUNCT__ 2678 #define __FUNCT__ "SNESReset" 2679 /*@ 2680 SNESReset - Resets a SNES context to the snessetupcalled = 0 state and removes any allocated Vecs and Mats 2681 2682 Collective on SNES 2683 2684 Input Parameter: 2685 . snes - iterative context obtained from SNESCreate() 2686 2687 Level: intermediate 2688 2689 Notes: Also calls the application context destroy routine set with SNESSetComputeApplicationContext() 2690 2691 .keywords: SNES, destroy 2692 2693 .seealso: SNESCreate(), SNESSetUp(), SNESSolve() 2694 @*/ 2695 PetscErrorCode SNESReset(SNES snes) 2696 { 2697 PetscErrorCode ierr; 2698 2699 PetscFunctionBegin; 2700 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2701 if (snes->ops->userdestroy && snes->user) { 2702 ierr = (*snes->ops->userdestroy)((void**)&snes->user);CHKERRQ(ierr); 2703 snes->user = NULL; 2704 } 2705 if (snes->pc) { 2706 ierr = SNESReset(snes->pc);CHKERRQ(ierr); 2707 } 2708 2709 if (snes->ops->reset) { 2710 ierr = (*snes->ops->reset)(snes);CHKERRQ(ierr); 2711 } 2712 if (snes->ksp) { 2713 ierr = KSPReset(snes->ksp);CHKERRQ(ierr); 2714 } 2715 2716 if (snes->linesearch) { 2717 ierr = SNESLineSearchReset(snes->linesearch);CHKERRQ(ierr); 2718 } 2719 2720 ierr = VecDestroy(&snes->vec_rhs);CHKERRQ(ierr); 2721 ierr = VecDestroy(&snes->vec_sol);CHKERRQ(ierr); 2722 ierr = VecDestroy(&snes->vec_sol_update);CHKERRQ(ierr); 2723 ierr = VecDestroy(&snes->vec_func);CHKERRQ(ierr); 2724 ierr = MatDestroy(&snes->jacobian);CHKERRQ(ierr); 2725 ierr = MatDestroy(&snes->jacobian_pre);CHKERRQ(ierr); 2726 ierr = VecDestroyVecs(snes->nwork,&snes->work);CHKERRQ(ierr); 2727 ierr = VecDestroyVecs(snes->nvwork,&snes->vwork);CHKERRQ(ierr); 2728 2729 snes->nwork = snes->nvwork = 0; 2730 snes->setupcalled = PETSC_FALSE; 2731 PetscFunctionReturn(0); 2732 } 2733 2734 #undef __FUNCT__ 2735 #define __FUNCT__ "SNESDestroy" 2736 /*@ 2737 SNESDestroy - Destroys the nonlinear solver context that was created 2738 with SNESCreate(). 2739 2740 Collective on SNES 2741 2742 Input Parameter: 2743 . snes - the SNES context 2744 2745 Level: beginner 2746 2747 .keywords: SNES, nonlinear, destroy 2748 2749 .seealso: SNESCreate(), SNESSolve() 2750 @*/ 2751 PetscErrorCode SNESDestroy(SNES *snes) 2752 { 2753 PetscErrorCode ierr; 2754 2755 PetscFunctionBegin; 2756 if (!*snes) PetscFunctionReturn(0); 2757 PetscValidHeaderSpecific((*snes),SNES_CLASSID,1); 2758 if (--((PetscObject)(*snes))->refct > 0) {*snes = 0; PetscFunctionReturn(0);} 2759 2760 ierr = SNESReset((*snes));CHKERRQ(ierr); 2761 ierr = SNESDestroy(&(*snes)->pc);CHKERRQ(ierr); 2762 2763 /* if memory was published with AMS then destroy it */ 2764 ierr = PetscObjectAMSViewOff((PetscObject)*snes);CHKERRQ(ierr); 2765 if ((*snes)->ops->destroy) {ierr = (*((*snes))->ops->destroy)((*snes));CHKERRQ(ierr);} 2766 2767 ierr = DMDestroy(&(*snes)->dm);CHKERRQ(ierr); 2768 ierr = KSPDestroy(&(*snes)->ksp);CHKERRQ(ierr); 2769 ierr = SNESLineSearchDestroy(&(*snes)->linesearch);CHKERRQ(ierr); 2770 2771 ierr = PetscFree((*snes)->kspconvctx);CHKERRQ(ierr); 2772 if ((*snes)->ops->convergeddestroy) { 2773 ierr = (*(*snes)->ops->convergeddestroy)((*snes)->cnvP);CHKERRQ(ierr); 2774 } 2775 if ((*snes)->conv_malloc) { 2776 ierr = PetscFree((*snes)->conv_hist);CHKERRQ(ierr); 2777 ierr = PetscFree((*snes)->conv_hist_its);CHKERRQ(ierr); 2778 } 2779 ierr = SNESMonitorCancel((*snes));CHKERRQ(ierr); 2780 ierr = PetscHeaderDestroy(snes);CHKERRQ(ierr); 2781 PetscFunctionReturn(0); 2782 } 2783 2784 /* ----------- Routines to set solver parameters ---------- */ 2785 2786 #undef __FUNCT__ 2787 #define __FUNCT__ "SNESSetLagPreconditioner" 2788 /*@ 2789 SNESSetLagPreconditioner - Determines when the preconditioner is rebuilt in the nonlinear solve. 2790 2791 Logically Collective on SNES 2792 2793 Input Parameters: 2794 + snes - the SNES context 2795 - lag - -1 indicates NEVER rebuild, 1 means rebuild every time the Jacobian is computed within a single nonlinear solve, 2 means every second time 2796 the Jacobian is built etc. -2 indicates rebuild preconditioner at next chance but then never rebuild after that 2797 2798 Options Database Keys: 2799 . -snes_lag_preconditioner <lag> 2800 2801 Notes: 2802 The default is 1 2803 The preconditioner is ALWAYS built in the first iteration of a nonlinear solve unless lag is -1 2804 If -1 is used before the very first nonlinear solve the preconditioner is still built because there is no previous preconditioner to use 2805 2806 Level: intermediate 2807 2808 .keywords: SNES, nonlinear, set, convergence, tolerances 2809 2810 .seealso: SNESSetTrustRegionTolerance(), SNESGetLagPreconditioner(), SNESSetLagJacobian(), SNESGetLagJacobian() 2811 2812 @*/ 2813 PetscErrorCode SNESSetLagPreconditioner(SNES snes,PetscInt lag) 2814 { 2815 PetscFunctionBegin; 2816 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2817 if (lag < -2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Lag must be -2, -1, 1 or greater"); 2818 if (!lag) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Lag cannot be 0"); 2819 PetscValidLogicalCollectiveInt(snes,lag,2); 2820 snes->lagpreconditioner = lag; 2821 PetscFunctionReturn(0); 2822 } 2823 2824 #undef __FUNCT__ 2825 #define __FUNCT__ "SNESSetGridSequence" 2826 /*@ 2827 SNESSetGridSequence - sets the number of steps of grid sequencing that SNES does 2828 2829 Logically Collective on SNES 2830 2831 Input Parameters: 2832 + snes - the SNES context 2833 - steps - the number of refinements to do, defaults to 0 2834 2835 Options Database Keys: 2836 . -snes_grid_sequence <steps> 2837 2838 Level: intermediate 2839 2840 Notes: 2841 Use SNESGetSolution() to extract the fine grid solution after grid sequencing. 2842 2843 .keywords: SNES, nonlinear, set, convergence, tolerances 2844 2845 .seealso: SNESSetTrustRegionTolerance(), SNESGetLagPreconditioner(), SNESSetLagJacobian(), SNESGetLagJacobian() 2846 2847 @*/ 2848 PetscErrorCode SNESSetGridSequence(SNES snes,PetscInt steps) 2849 { 2850 PetscFunctionBegin; 2851 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2852 PetscValidLogicalCollectiveInt(snes,steps,2); 2853 snes->gridsequence = steps; 2854 PetscFunctionReturn(0); 2855 } 2856 2857 #undef __FUNCT__ 2858 #define __FUNCT__ "SNESGetLagPreconditioner" 2859 /*@ 2860 SNESGetLagPreconditioner - Indicates how often the preconditioner is rebuilt 2861 2862 Not Collective 2863 2864 Input Parameter: 2865 . snes - the SNES context 2866 2867 Output Parameter: 2868 . lag - -1 indicates NEVER rebuild, 1 means rebuild every time the Jacobian is computed within a single nonlinear solve, 2 means every second time 2869 the Jacobian is built etc. -2 indicates rebuild preconditioner at next chance but then never rebuild after that 2870 2871 Options Database Keys: 2872 . -snes_lag_preconditioner <lag> 2873 2874 Notes: 2875 The default is 1 2876 The preconditioner is ALWAYS built in the first iteration of a nonlinear solve unless lag is -1 2877 2878 Level: intermediate 2879 2880 .keywords: SNES, nonlinear, set, convergence, tolerances 2881 2882 .seealso: SNESSetTrustRegionTolerance(), SNESSetLagPreconditioner() 2883 2884 @*/ 2885 PetscErrorCode SNESGetLagPreconditioner(SNES snes,PetscInt *lag) 2886 { 2887 PetscFunctionBegin; 2888 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2889 *lag = snes->lagpreconditioner; 2890 PetscFunctionReturn(0); 2891 } 2892 2893 #undef __FUNCT__ 2894 #define __FUNCT__ "SNESSetLagJacobian" 2895 /*@ 2896 SNESSetLagJacobian - Determines when the Jacobian is rebuilt in the nonlinear solve. See SNESSetLagPreconditioner() for determining how 2897 often the preconditioner is rebuilt. 2898 2899 Logically Collective on SNES 2900 2901 Input Parameters: 2902 + snes - the SNES context 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 means rebuild at next chance but then never again 2905 2906 Options Database Keys: 2907 . -snes_lag_jacobian <lag> 2908 2909 Notes: 2910 The default is 1 2911 The Jacobian is ALWAYS built in the first iteration of a nonlinear solve unless lag is -1 2912 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 2913 at the next Newton step but never again (unless it is reset to another value) 2914 2915 Level: intermediate 2916 2917 .keywords: SNES, nonlinear, set, convergence, tolerances 2918 2919 .seealso: SNESSetTrustRegionTolerance(), SNESGetLagPreconditioner(), SNESSetLagPreconditioner(), SNESGetLagJacobian() 2920 2921 @*/ 2922 PetscErrorCode SNESSetLagJacobian(SNES snes,PetscInt lag) 2923 { 2924 PetscFunctionBegin; 2925 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2926 if (lag < -2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Lag must be -2, -1, 1 or greater"); 2927 if (!lag) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Lag cannot be 0"); 2928 PetscValidLogicalCollectiveInt(snes,lag,2); 2929 snes->lagjacobian = lag; 2930 PetscFunctionReturn(0); 2931 } 2932 2933 #undef __FUNCT__ 2934 #define __FUNCT__ "SNESGetLagJacobian" 2935 /*@ 2936 SNESGetLagJacobian - Indicates how often the Jacobian is rebuilt. See SNESGetLagPreconditioner() to determine when the preconditioner is rebuilt 2937 2938 Not Collective 2939 2940 Input Parameter: 2941 . snes - the SNES context 2942 2943 Output Parameter: 2944 . lag - -1 indicates NEVER rebuild, 1 means rebuild every time the Jacobian is computed within a single nonlinear solve, 2 means every second time 2945 the Jacobian is built etc. 2946 2947 Options Database Keys: 2948 . -snes_lag_jacobian <lag> 2949 2950 Notes: 2951 The default is 1 2952 The jacobian is ALWAYS built in the first iteration of a nonlinear solve unless lag is -1 2953 2954 Level: intermediate 2955 2956 .keywords: SNES, nonlinear, set, convergence, tolerances 2957 2958 .seealso: SNESSetTrustRegionTolerance(), SNESSetLagJacobian(), SNESSetLagPreconditioner(), SNESGetLagPreconditioner() 2959 2960 @*/ 2961 PetscErrorCode SNESGetLagJacobian(SNES snes,PetscInt *lag) 2962 { 2963 PetscFunctionBegin; 2964 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2965 *lag = snes->lagjacobian; 2966 PetscFunctionReturn(0); 2967 } 2968 2969 #undef __FUNCT__ 2970 #define __FUNCT__ "SNESSetLagJacobianPersists" 2971 /*@ 2972 SNESSetLagJacobianPersists - Set whether or not the Jacobian lagging persists through multiple solves 2973 2974 Logically collective on SNES 2975 2976 Input Parameter: 2977 + snes - the SNES context 2978 - flg - jacobian lagging persists if true 2979 2980 Options Database Keys: 2981 . -snes_lag_jacobian_persists <flg> 2982 2983 Notes: This is useful both for nonlinear preconditioning, where it's appropriate to have the Jacobian be stale by 2984 several solves, and for implicit time-stepping, where Jacobian lagging in the inner nonlinear solve over several 2985 timesteps may present huge efficiency gains. 2986 2987 Level: developer 2988 2989 .keywords: SNES, nonlinear, lag, NPC 2990 2991 .seealso: SNESSetLagPreconditionerPersists(), SNESSetLagJacobian(), SNESGetLagJacobian(), SNESGetPC() 2992 2993 @*/ 2994 PetscErrorCode SNESSetLagJacobianPersists(SNES snes,PetscBool flg) 2995 { 2996 PetscFunctionBegin; 2997 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2998 PetscValidLogicalCollectiveBool(snes,flg,2); 2999 snes->lagjac_persist = flg; 3000 PetscFunctionReturn(0); 3001 } 3002 3003 #undef __FUNCT__ 3004 #define __FUNCT__ "SNESSetLagPreconditionerPersists" 3005 /*@ 3006 SNESSetLagPreconditionerPersists - Set whether or not the preconditioner lagging persists through multiple solves 3007 3008 Logically Collective on SNES 3009 3010 Input Parameter: 3011 + snes - the SNES context 3012 - flg - preconditioner lagging persists if true 3013 3014 Options Database Keys: 3015 . -snes_lag_jacobian_persists <flg> 3016 3017 Notes: This is useful both for nonlinear preconditioning, where it's appropriate to have the preconditioner be stale 3018 by several solves, and for implicit time-stepping, where preconditioner lagging in the inner nonlinear solve over 3019 several timesteps may present huge efficiency gains. 3020 3021 Level: developer 3022 3023 .keywords: SNES, nonlinear, lag, NPC 3024 3025 .seealso: SNESSetLagJacobianPersists(), SNESSetLagJacobian(), SNESGetLagJacobian(), SNESGetPC() 3026 3027 @*/ 3028 PetscErrorCode SNESSetLagPreconditionerPersists(SNES snes,PetscBool flg) 3029 { 3030 PetscFunctionBegin; 3031 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3032 PetscValidLogicalCollectiveBool(snes,flg,2); 3033 snes->lagpre_persist = flg; 3034 PetscFunctionReturn(0); 3035 } 3036 3037 #undef __FUNCT__ 3038 #define __FUNCT__ "SNESSetTolerances" 3039 /*@ 3040 SNESSetTolerances - Sets various parameters used in convergence tests. 3041 3042 Logically Collective on SNES 3043 3044 Input Parameters: 3045 + snes - the SNES context 3046 . abstol - absolute convergence tolerance 3047 . rtol - relative convergence tolerance 3048 . stol - convergence tolerance in terms of the norm of the change in the solution between steps, || delta x || < stol*|| x || 3049 . maxit - maximum number of iterations 3050 - maxf - maximum number of function evaluations 3051 3052 Options Database Keys: 3053 + -snes_atol <abstol> - Sets abstol 3054 . -snes_rtol <rtol> - Sets rtol 3055 . -snes_stol <stol> - Sets stol 3056 . -snes_max_it <maxit> - Sets maxit 3057 - -snes_max_funcs <maxf> - Sets maxf 3058 3059 Notes: 3060 The default maximum number of iterations is 50. 3061 The default maximum number of function evaluations is 1000. 3062 3063 Level: intermediate 3064 3065 .keywords: SNES, nonlinear, set, convergence, tolerances 3066 3067 .seealso: SNESSetTrustRegionTolerance() 3068 @*/ 3069 PetscErrorCode SNESSetTolerances(SNES snes,PetscReal abstol,PetscReal rtol,PetscReal stol,PetscInt maxit,PetscInt maxf) 3070 { 3071 PetscFunctionBegin; 3072 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3073 PetscValidLogicalCollectiveReal(snes,abstol,2); 3074 PetscValidLogicalCollectiveReal(snes,rtol,3); 3075 PetscValidLogicalCollectiveReal(snes,stol,4); 3076 PetscValidLogicalCollectiveInt(snes,maxit,5); 3077 PetscValidLogicalCollectiveInt(snes,maxf,6); 3078 3079 if (abstol != PETSC_DEFAULT) { 3080 if (abstol < 0.0) SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_OUTOFRANGE,"Absolute tolerance %G must be non-negative",abstol); 3081 snes->abstol = abstol; 3082 } 3083 if (rtol != PETSC_DEFAULT) { 3084 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); 3085 snes->rtol = rtol; 3086 } 3087 if (stol != PETSC_DEFAULT) { 3088 if (stol < 0.0) SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_OUTOFRANGE,"Step tolerance %G must be non-negative",stol); 3089 snes->stol = stol; 3090 } 3091 if (maxit != PETSC_DEFAULT) { 3092 if (maxit < 0) SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_OUTOFRANGE,"Maximum number of iterations %D must be non-negative",maxit); 3093 snes->max_its = maxit; 3094 } 3095 if (maxf != PETSC_DEFAULT) { 3096 if (maxf < 0) SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_OUTOFRANGE,"Maximum number of function evaluations %D must be non-negative",maxf); 3097 snes->max_funcs = maxf; 3098 } 3099 snes->tolerancesset = PETSC_TRUE; 3100 PetscFunctionReturn(0); 3101 } 3102 3103 #undef __FUNCT__ 3104 #define __FUNCT__ "SNESGetTolerances" 3105 /*@ 3106 SNESGetTolerances - Gets various parameters used in convergence tests. 3107 3108 Not Collective 3109 3110 Input Parameters: 3111 + snes - the SNES context 3112 . atol - absolute convergence tolerance 3113 . rtol - relative convergence tolerance 3114 . stol - convergence tolerance in terms of the norm 3115 of the change in the solution between steps 3116 . maxit - maximum number of iterations 3117 - maxf - maximum number of function evaluations 3118 3119 Notes: 3120 The user can specify NULL for any parameter that is not needed. 3121 3122 Level: intermediate 3123 3124 .keywords: SNES, nonlinear, get, convergence, tolerances 3125 3126 .seealso: SNESSetTolerances() 3127 @*/ 3128 PetscErrorCode SNESGetTolerances(SNES snes,PetscReal *atol,PetscReal *rtol,PetscReal *stol,PetscInt *maxit,PetscInt *maxf) 3129 { 3130 PetscFunctionBegin; 3131 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3132 if (atol) *atol = snes->abstol; 3133 if (rtol) *rtol = snes->rtol; 3134 if (stol) *stol = snes->stol; 3135 if (maxit) *maxit = snes->max_its; 3136 if (maxf) *maxf = snes->max_funcs; 3137 PetscFunctionReturn(0); 3138 } 3139 3140 #undef __FUNCT__ 3141 #define __FUNCT__ "SNESSetTrustRegionTolerance" 3142 /*@ 3143 SNESSetTrustRegionTolerance - Sets the trust region parameter tolerance. 3144 3145 Logically Collective on SNES 3146 3147 Input Parameters: 3148 + snes - the SNES context 3149 - tol - tolerance 3150 3151 Options Database Key: 3152 . -snes_trtol <tol> - Sets tol 3153 3154 Level: intermediate 3155 3156 .keywords: SNES, nonlinear, set, trust region, tolerance 3157 3158 .seealso: SNESSetTolerances() 3159 @*/ 3160 PetscErrorCode SNESSetTrustRegionTolerance(SNES snes,PetscReal tol) 3161 { 3162 PetscFunctionBegin; 3163 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3164 PetscValidLogicalCollectiveReal(snes,tol,2); 3165 snes->deltatol = tol; 3166 PetscFunctionReturn(0); 3167 } 3168 3169 /* 3170 Duplicate the lg monitors for SNES from KSP; for some reason with 3171 dynamic libraries things don't work under Sun4 if we just use 3172 macros instead of functions 3173 */ 3174 #undef __FUNCT__ 3175 #define __FUNCT__ "SNESMonitorLGResidualNorm" 3176 PetscErrorCode SNESMonitorLGResidualNorm(SNES snes,PetscInt it,PetscReal norm,void *ctx) 3177 { 3178 PetscErrorCode ierr; 3179 3180 PetscFunctionBegin; 3181 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3182 ierr = KSPMonitorLGResidualNorm((KSP)snes,it,norm,ctx);CHKERRQ(ierr); 3183 PetscFunctionReturn(0); 3184 } 3185 3186 #undef __FUNCT__ 3187 #define __FUNCT__ "SNESMonitorLGCreate" 3188 PetscErrorCode SNESMonitorLGCreate(const char host[],const char label[],int x,int y,int m,int n,PetscDrawLG *draw) 3189 { 3190 PetscErrorCode ierr; 3191 3192 PetscFunctionBegin; 3193 ierr = KSPMonitorLGResidualNormCreate(host,label,x,y,m,n,draw);CHKERRQ(ierr); 3194 PetscFunctionReturn(0); 3195 } 3196 3197 #undef __FUNCT__ 3198 #define __FUNCT__ "SNESMonitorLGDestroy" 3199 PetscErrorCode SNESMonitorLGDestroy(PetscDrawLG *draw) 3200 { 3201 PetscErrorCode ierr; 3202 3203 PetscFunctionBegin; 3204 ierr = KSPMonitorLGResidualNormDestroy(draw);CHKERRQ(ierr); 3205 PetscFunctionReturn(0); 3206 } 3207 3208 extern PetscErrorCode SNESMonitorRange_Private(SNES,PetscInt,PetscReal*); 3209 #undef __FUNCT__ 3210 #define __FUNCT__ "SNESMonitorLGRange" 3211 PetscErrorCode SNESMonitorLGRange(SNES snes,PetscInt n,PetscReal rnorm,void *monctx) 3212 { 3213 PetscDrawLG lg; 3214 PetscErrorCode ierr; 3215 PetscReal x,y,per; 3216 PetscViewer v = (PetscViewer)monctx; 3217 static PetscReal prev; /* should be in the context */ 3218 PetscDraw draw; 3219 3220 PetscFunctionBegin; 3221 ierr = PetscViewerDrawGetDrawLG(v,0,&lg);CHKERRQ(ierr); 3222 if (!n) {ierr = PetscDrawLGReset(lg);CHKERRQ(ierr);} 3223 ierr = PetscDrawLGGetDraw(lg,&draw);CHKERRQ(ierr); 3224 ierr = PetscDrawSetTitle(draw,"Residual norm");CHKERRQ(ierr); 3225 x = (PetscReal)n; 3226 if (rnorm > 0.0) y = log10(rnorm); 3227 else y = -15.0; 3228 ierr = PetscDrawLGAddPoint(lg,&x,&y);CHKERRQ(ierr); 3229 if (n < 20 || !(n % 5)) { 3230 ierr = PetscDrawLGDraw(lg);CHKERRQ(ierr); 3231 } 3232 3233 ierr = PetscViewerDrawGetDrawLG(v,1,&lg);CHKERRQ(ierr); 3234 if (!n) {ierr = PetscDrawLGReset(lg);CHKERRQ(ierr);} 3235 ierr = PetscDrawLGGetDraw(lg,&draw);CHKERRQ(ierr); 3236 ierr = PetscDrawSetTitle(draw,"% elemts > .2*max elemt");CHKERRQ(ierr); 3237 ierr = SNESMonitorRange_Private(snes,n,&per);CHKERRQ(ierr); 3238 x = (PetscReal)n; 3239 y = 100.0*per; 3240 ierr = PetscDrawLGAddPoint(lg,&x,&y);CHKERRQ(ierr); 3241 if (n < 20 || !(n % 5)) { 3242 ierr = PetscDrawLGDraw(lg);CHKERRQ(ierr); 3243 } 3244 3245 ierr = PetscViewerDrawGetDrawLG(v,2,&lg);CHKERRQ(ierr); 3246 if (!n) {prev = rnorm;ierr = PetscDrawLGReset(lg);CHKERRQ(ierr);} 3247 ierr = PetscDrawLGGetDraw(lg,&draw);CHKERRQ(ierr); 3248 ierr = PetscDrawSetTitle(draw,"(norm -oldnorm)/oldnorm");CHKERRQ(ierr); 3249 x = (PetscReal)n; 3250 y = (prev - rnorm)/prev; 3251 ierr = PetscDrawLGAddPoint(lg,&x,&y);CHKERRQ(ierr); 3252 if (n < 20 || !(n % 5)) { 3253 ierr = PetscDrawLGDraw(lg);CHKERRQ(ierr); 3254 } 3255 3256 ierr = PetscViewerDrawGetDrawLG(v,3,&lg);CHKERRQ(ierr); 3257 if (!n) {ierr = PetscDrawLGReset(lg);CHKERRQ(ierr);} 3258 ierr = PetscDrawLGGetDraw(lg,&draw);CHKERRQ(ierr); 3259 ierr = PetscDrawSetTitle(draw,"(norm -oldnorm)/oldnorm*(% > .2 max)");CHKERRQ(ierr); 3260 x = (PetscReal)n; 3261 y = (prev - rnorm)/(prev*per); 3262 if (n > 2) { /*skip initial crazy value */ 3263 ierr = PetscDrawLGAddPoint(lg,&x,&y);CHKERRQ(ierr); 3264 } 3265 if (n < 20 || !(n % 5)) { 3266 ierr = PetscDrawLGDraw(lg);CHKERRQ(ierr); 3267 } 3268 prev = rnorm; 3269 PetscFunctionReturn(0); 3270 } 3271 3272 #undef __FUNCT__ 3273 #define __FUNCT__ "SNESMonitor" 3274 /*@ 3275 SNESMonitor - runs the user provided monitor routines, if they exist 3276 3277 Collective on SNES 3278 3279 Input Parameters: 3280 + snes - nonlinear solver context obtained from SNESCreate() 3281 . iter - iteration number 3282 - rnorm - relative norm of the residual 3283 3284 Notes: 3285 This routine is called by the SNES implementations. 3286 It does not typically need to be called by the user. 3287 3288 Level: developer 3289 3290 .seealso: SNESMonitorSet() 3291 @*/ 3292 PetscErrorCode SNESMonitor(SNES snes,PetscInt iter,PetscReal rnorm) 3293 { 3294 PetscErrorCode ierr; 3295 PetscInt i,n = snes->numbermonitors; 3296 3297 PetscFunctionBegin; 3298 for (i=0; i<n; i++) { 3299 ierr = (*snes->monitor[i])(snes,iter,rnorm,snes->monitorcontext[i]);CHKERRQ(ierr); 3300 } 3301 PetscFunctionReturn(0); 3302 } 3303 3304 /* ------------ Routines to set performance monitoring options ----------- */ 3305 3306 /*MC 3307 SNESMonitorFunction - functional form passed to SNESMonitorSet() to monitor convergence of nonlinear solver 3308 3309 Synopsis: 3310 #include "petscsnes.h" 3311 $ PetscErrorCode SNESMonitorFunction(SNES snes,PetscInt its, PetscReal norm,void *mctx) 3312 3313 + snes - the SNES context 3314 . its - iteration number 3315 . norm - 2-norm function value (may be estimated) 3316 - mctx - [optional] monitoring context 3317 3318 Level: advanced 3319 3320 .seealso: SNESMonitorSet(), SNESMonitorGet() 3321 M*/ 3322 3323 #undef __FUNCT__ 3324 #define __FUNCT__ "SNESMonitorSet" 3325 /*@C 3326 SNESMonitorSet - Sets an ADDITIONAL function that is to be used at every 3327 iteration of the nonlinear solver to display the iteration's 3328 progress. 3329 3330 Logically Collective on SNES 3331 3332 Input Parameters: 3333 + snes - the SNES context 3334 . SNESMonitorFunction - monitoring routine 3335 . mctx - [optional] user-defined context for private data for the 3336 monitor routine (use NULL if no context is desired) 3337 - monitordestroy - [optional] routine that frees monitor context 3338 (may be NULL) 3339 3340 Options Database Keys: 3341 + -snes_monitor - sets SNESMonitorDefault() 3342 . -snes_monitor_lg_residualnorm - sets line graph monitor, 3343 uses SNESMonitorLGCreate() 3344 - -snes_monitor_cancel - cancels all monitors that have 3345 been hardwired into a code by 3346 calls to SNESMonitorSet(), but 3347 does not cancel those set via 3348 the options database. 3349 3350 Notes: 3351 Several different monitoring routines may be set by calling 3352 SNESMonitorSet() multiple times; all will be called in the 3353 order in which they were set. 3354 3355 Fortran notes: Only a single monitor function can be set for each SNES object 3356 3357 Level: intermediate 3358 3359 .keywords: SNES, nonlinear, set, monitor 3360 3361 .seealso: SNESMonitorDefault(), SNESMonitorCancel(), SNESMonitorFunction 3362 @*/ 3363 PetscErrorCode SNESMonitorSet(SNES snes,PetscErrorCode (*SNESMonitorFunction)(SNES,PetscInt,PetscReal,void*),void *mctx,PetscErrorCode (*monitordestroy)(void**)) 3364 { 3365 PetscInt i; 3366 PetscErrorCode ierr; 3367 3368 PetscFunctionBegin; 3369 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3370 if (snes->numbermonitors >= MAXSNESMONITORS) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Too many monitors set"); 3371 for (i=0; i<snes->numbermonitors;i++) { 3372 if (SNESMonitorFunction == snes->monitor[i] && monitordestroy == snes->monitordestroy[i] && mctx == snes->monitorcontext[i]) { 3373 if (monitordestroy) { 3374 ierr = (*monitordestroy)(&mctx);CHKERRQ(ierr); 3375 } 3376 PetscFunctionReturn(0); 3377 } 3378 } 3379 snes->monitor[snes->numbermonitors] = SNESMonitorFunction; 3380 snes->monitordestroy[snes->numbermonitors] = monitordestroy; 3381 snes->monitorcontext[snes->numbermonitors++] = (void*)mctx; 3382 PetscFunctionReturn(0); 3383 } 3384 3385 #undef __FUNCT__ 3386 #define __FUNCT__ "SNESMonitorCancel" 3387 /*@ 3388 SNESMonitorCancel - Clears all the monitor functions for a SNES object. 3389 3390 Logically Collective on SNES 3391 3392 Input Parameters: 3393 . snes - the SNES context 3394 3395 Options Database Key: 3396 . -snes_monitor_cancel - cancels all monitors that have been hardwired 3397 into a code by calls to SNESMonitorSet(), but does not cancel those 3398 set via the options database 3399 3400 Notes: 3401 There is no way to clear one specific monitor from a SNES object. 3402 3403 Level: intermediate 3404 3405 .keywords: SNES, nonlinear, set, monitor 3406 3407 .seealso: SNESMonitorDefault(), SNESMonitorSet() 3408 @*/ 3409 PetscErrorCode SNESMonitorCancel(SNES snes) 3410 { 3411 PetscErrorCode ierr; 3412 PetscInt i; 3413 3414 PetscFunctionBegin; 3415 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3416 for (i=0; i<snes->numbermonitors; i++) { 3417 if (snes->monitordestroy[i]) { 3418 ierr = (*snes->monitordestroy[i])(&snes->monitorcontext[i]);CHKERRQ(ierr); 3419 } 3420 } 3421 snes->numbermonitors = 0; 3422 PetscFunctionReturn(0); 3423 } 3424 3425 /*MC 3426 SNESConvergenceTestFunction - functional form used for testing of convergence of nonlinear solver 3427 3428 Synopsis: 3429 #include "petscsnes.h" 3430 $ PetscErrorCode SNESConvergenceTest(SNES snes,PetscInt it,PetscReal xnorm,PetscReal gnorm,PetscReal f,SNESConvergedReason *reason,void *cctx) 3431 3432 + snes - the SNES context 3433 . it - current iteration (0 is the first and is before any Newton step) 3434 . cctx - [optional] convergence context 3435 . reason - reason for convergence/divergence 3436 . xnorm - 2-norm of current iterate 3437 . gnorm - 2-norm of current step 3438 - f - 2-norm of function 3439 3440 Level: intermediate 3441 3442 .seealso: SNESSetConvergenceTest(), SNESGetConvergenceTest() 3443 M*/ 3444 3445 #undef __FUNCT__ 3446 #define __FUNCT__ "SNESSetConvergenceTest" 3447 /*@C 3448 SNESSetConvergenceTest - Sets the function that is to be used 3449 to test for convergence of the nonlinear iterative solution. 3450 3451 Logically Collective on SNES 3452 3453 Input Parameters: 3454 + snes - the SNES context 3455 . SNESConvergenceTestFunction - routine to test for convergence 3456 . cctx - [optional] context for private data for the convergence routine (may be NULL) 3457 - destroy - [optional] destructor for the context (may be NULL; NULL_FUNCTION in Fortran) 3458 3459 Level: advanced 3460 3461 .keywords: SNES, nonlinear, set, convergence, test 3462 3463 .seealso: SNESConvergedDefault(), SNESConvergedSkip(), SNESConvergenceTestFunction 3464 @*/ 3465 PetscErrorCode SNESSetConvergenceTest(SNES snes,PetscErrorCode (*SNESConvergenceTestFunction)(SNES,PetscInt,PetscReal,PetscReal,PetscReal,SNESConvergedReason*,void*),void *cctx,PetscErrorCode (*destroy)(void*)) 3466 { 3467 PetscErrorCode ierr; 3468 3469 PetscFunctionBegin; 3470 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3471 if (!SNESConvergenceTestFunction) SNESConvergenceTestFunction = SNESConvergedSkip; 3472 if (snes->ops->convergeddestroy) { 3473 ierr = (*snes->ops->convergeddestroy)(snes->cnvP);CHKERRQ(ierr); 3474 } 3475 snes->ops->converged = SNESConvergenceTestFunction; 3476 snes->ops->convergeddestroy = destroy; 3477 snes->cnvP = cctx; 3478 PetscFunctionReturn(0); 3479 } 3480 3481 #undef __FUNCT__ 3482 #define __FUNCT__ "SNESGetConvergedReason" 3483 /*@ 3484 SNESGetConvergedReason - Gets the reason the SNES iteration was stopped. 3485 3486 Not Collective 3487 3488 Input Parameter: 3489 . snes - the SNES context 3490 3491 Output Parameter: 3492 . reason - negative value indicates diverged, positive value converged, see SNESConvergedReason or the 3493 manual pages for the individual convergence tests for complete lists 3494 3495 Level: intermediate 3496 3497 Notes: Can only be called after the call the SNESSolve() is complete. 3498 3499 .keywords: SNES, nonlinear, set, convergence, test 3500 3501 .seealso: SNESSetConvergenceTest(), SNESConvergedReason 3502 @*/ 3503 PetscErrorCode SNESGetConvergedReason(SNES snes,SNESConvergedReason *reason) 3504 { 3505 PetscFunctionBegin; 3506 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3507 PetscValidPointer(reason,2); 3508 *reason = snes->reason; 3509 PetscFunctionReturn(0); 3510 } 3511 3512 #undef __FUNCT__ 3513 #define __FUNCT__ "SNESSetConvergenceHistory" 3514 /*@ 3515 SNESSetConvergenceHistory - Sets the array used to hold the convergence history. 3516 3517 Logically Collective on SNES 3518 3519 Input Parameters: 3520 + snes - iterative context obtained from SNESCreate() 3521 . a - array to hold history, this array will contain the function norms computed at each step 3522 . its - integer array holds the number of linear iterations for each solve. 3523 . na - size of a and its 3524 - reset - PETSC_TRUE indicates each new nonlinear solve resets the history counter to zero, 3525 else it continues storing new values for new nonlinear solves after the old ones 3526 3527 Notes: 3528 If 'a' and 'its' are NULL then space is allocated for the history. If 'na' PETSC_DECIDE or PETSC_DEFAULT then a 3529 default array of length 10000 is allocated. 3530 3531 This routine is useful, e.g., when running a code for purposes 3532 of accurate performance monitoring, when no I/O should be done 3533 during the section of code that is being timed. 3534 3535 Level: intermediate 3536 3537 .keywords: SNES, set, convergence, history 3538 3539 .seealso: SNESGetConvergenceHistory() 3540 3541 @*/ 3542 PetscErrorCode SNESSetConvergenceHistory(SNES snes,PetscReal a[],PetscInt its[],PetscInt na,PetscBool reset) 3543 { 3544 PetscErrorCode ierr; 3545 3546 PetscFunctionBegin; 3547 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3548 if (a) PetscValidScalarPointer(a,2); 3549 if (its) PetscValidIntPointer(its,3); 3550 if (!a) { 3551 if (na == PETSC_DECIDE || na == PETSC_DEFAULT) na = 1000; 3552 ierr = PetscMalloc(na*sizeof(PetscReal),&a);CHKERRQ(ierr); 3553 ierr = PetscMalloc(na*sizeof(PetscInt),&its);CHKERRQ(ierr); 3554 3555 snes->conv_malloc = PETSC_TRUE; 3556 } 3557 snes->conv_hist = a; 3558 snes->conv_hist_its = its; 3559 snes->conv_hist_max = na; 3560 snes->conv_hist_len = 0; 3561 snes->conv_hist_reset = reset; 3562 PetscFunctionReturn(0); 3563 } 3564 3565 #if defined(PETSC_HAVE_MATLAB_ENGINE) 3566 #include <engine.h> /* MATLAB include file */ 3567 #include <mex.h> /* MATLAB include file */ 3568 3569 #undef __FUNCT__ 3570 #define __FUNCT__ "SNESGetConvergenceHistoryMatlab" 3571 PETSC_EXTERN mxArray *SNESGetConvergenceHistoryMatlab(SNES snes) 3572 { 3573 mxArray *mat; 3574 PetscInt i; 3575 PetscReal *ar; 3576 3577 PetscFunctionBegin; 3578 mat = mxCreateDoubleMatrix(snes->conv_hist_len,1,mxREAL); 3579 ar = (PetscReal*) mxGetData(mat); 3580 for (i=0; i<snes->conv_hist_len; i++) ar[i] = snes->conv_hist[i]; 3581 PetscFunctionReturn(mat); 3582 } 3583 #endif 3584 3585 #undef __FUNCT__ 3586 #define __FUNCT__ "SNESGetConvergenceHistory" 3587 /*@C 3588 SNESGetConvergenceHistory - Gets the array used to hold the convergence history. 3589 3590 Not Collective 3591 3592 Input Parameter: 3593 . snes - iterative context obtained from SNESCreate() 3594 3595 Output Parameters: 3596 . a - array to hold history 3597 . its - integer array holds the number of linear iterations (or 3598 negative if not converged) for each solve. 3599 - na - size of a and its 3600 3601 Notes: 3602 The calling sequence for this routine in Fortran is 3603 $ call SNESGetConvergenceHistory(SNES snes, integer na, integer ierr) 3604 3605 This routine is useful, e.g., when running a code for purposes 3606 of accurate performance monitoring, when no I/O should be done 3607 during the section of code that is being timed. 3608 3609 Level: intermediate 3610 3611 .keywords: SNES, get, convergence, history 3612 3613 .seealso: SNESSetConvergencHistory() 3614 3615 @*/ 3616 PetscErrorCode SNESGetConvergenceHistory(SNES snes,PetscReal *a[],PetscInt *its[],PetscInt *na) 3617 { 3618 PetscFunctionBegin; 3619 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3620 if (a) *a = snes->conv_hist; 3621 if (its) *its = snes->conv_hist_its; 3622 if (na) *na = snes->conv_hist_len; 3623 PetscFunctionReturn(0); 3624 } 3625 3626 #undef __FUNCT__ 3627 #define __FUNCT__ "SNESSetUpdate" 3628 /*@C 3629 SNESSetUpdate - Sets the general-purpose update function called 3630 at the beginning of every iteration of the nonlinear solve. Specifically 3631 it is called just before the Jacobian is "evaluated". 3632 3633 Logically Collective on SNES 3634 3635 Input Parameters: 3636 . snes - The nonlinear solver context 3637 . func - The function 3638 3639 Calling sequence of func: 3640 . func (SNES snes, PetscInt step); 3641 3642 . step - The current step of the iteration 3643 3644 Level: advanced 3645 3646 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() 3647 This is not used by most users. 3648 3649 .keywords: SNES, update 3650 3651 .seealso SNESSetJacobian(), SNESSolve() 3652 @*/ 3653 PetscErrorCode SNESSetUpdate(SNES snes, PetscErrorCode (*func)(SNES, PetscInt)) 3654 { 3655 PetscFunctionBegin; 3656 PetscValidHeaderSpecific(snes, SNES_CLASSID,1); 3657 snes->ops->update = func; 3658 PetscFunctionReturn(0); 3659 } 3660 3661 #undef __FUNCT__ 3662 #define __FUNCT__ "SNESScaleStep_Private" 3663 /* 3664 SNESScaleStep_Private - Scales a step so that its length is less than the 3665 positive parameter delta. 3666 3667 Input Parameters: 3668 + snes - the SNES context 3669 . y - approximate solution of linear system 3670 . fnorm - 2-norm of current function 3671 - delta - trust region size 3672 3673 Output Parameters: 3674 + gpnorm - predicted function norm at the new point, assuming local 3675 linearization. The value is zero if the step lies within the trust 3676 region, and exceeds zero otherwise. 3677 - ynorm - 2-norm of the step 3678 3679 Note: 3680 For non-trust region methods such as SNESNEWTONLS, the parameter delta 3681 is set to be the maximum allowable step size. 3682 3683 .keywords: SNES, nonlinear, scale, step 3684 */ 3685 PetscErrorCode SNESScaleStep_Private(SNES snes,Vec y,PetscReal *fnorm,PetscReal *delta,PetscReal *gpnorm,PetscReal *ynorm) 3686 { 3687 PetscReal nrm; 3688 PetscScalar cnorm; 3689 PetscErrorCode ierr; 3690 3691 PetscFunctionBegin; 3692 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3693 PetscValidHeaderSpecific(y,VEC_CLASSID,2); 3694 PetscCheckSameComm(snes,1,y,2); 3695 3696 ierr = VecNorm(y,NORM_2,&nrm);CHKERRQ(ierr); 3697 if (nrm > *delta) { 3698 nrm = *delta/nrm; 3699 *gpnorm = (1.0 - nrm)*(*fnorm); 3700 cnorm = nrm; 3701 ierr = VecScale(y,cnorm);CHKERRQ(ierr); 3702 *ynorm = *delta; 3703 } else { 3704 *gpnorm = 0.0; 3705 *ynorm = nrm; 3706 } 3707 PetscFunctionReturn(0); 3708 } 3709 3710 #undef __FUNCT__ 3711 #define __FUNCT__ "SNESSolve" 3712 /*@C 3713 SNESSolve - Solves a nonlinear system F(x) = b. 3714 Call SNESSolve() after calling SNESCreate() and optional routines of the form SNESSetXXX(). 3715 3716 Collective on SNES 3717 3718 Input Parameters: 3719 + snes - the SNES context 3720 . b - the constant part of the equation F(x) = b, or NULL to use zero. 3721 - x - the solution vector. 3722 3723 Notes: 3724 The user should initialize the vector,x, with the initial guess 3725 for the nonlinear solve prior to calling SNESSolve. In particular, 3726 to employ an initial guess of zero, the user should explicitly set 3727 this vector to zero by calling VecSet(). 3728 3729 Level: beginner 3730 3731 .keywords: SNES, nonlinear, solve 3732 3733 .seealso: SNESCreate(), SNESDestroy(), SNESSetFunction(), SNESSetJacobian(), SNESSetGridSequence(), SNESGetSolution() 3734 @*/ 3735 PetscErrorCode SNESSolve(SNES snes,Vec b,Vec x) 3736 { 3737 PetscErrorCode ierr; 3738 PetscBool flg; 3739 PetscViewer viewer; 3740 PetscInt grid; 3741 Vec xcreated = NULL; 3742 DM dm; 3743 PetscViewerFormat format; 3744 3745 PetscFunctionBegin; 3746 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3747 if (x) PetscValidHeaderSpecific(x,VEC_CLASSID,3); 3748 if (x) PetscCheckSameComm(snes,1,x,3); 3749 if (b) PetscValidHeaderSpecific(b,VEC_CLASSID,2); 3750 if (b) PetscCheckSameComm(snes,1,b,2); 3751 3752 if (!x) { 3753 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 3754 ierr = DMCreateGlobalVector(dm,&xcreated);CHKERRQ(ierr); 3755 x = xcreated; 3756 } 3757 3758 ierr = PetscOptionsGetViewer(PetscObjectComm((PetscObject)snes),((PetscObject)snes)->prefix,"-snes_view_pre",&viewer,&format,&flg);CHKERRQ(ierr); 3759 if (flg && !PetscPreLoadingOn) { 3760 ierr = PetscViewerPushFormat(viewer,format);CHKERRQ(ierr); 3761 ierr = SNESView(snes,viewer);CHKERRQ(ierr); 3762 ierr = PetscViewerPopFormat(viewer);CHKERRQ(ierr); 3763 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 3764 } 3765 3766 for (grid=0; grid<snes->gridsequence; grid++) {ierr = PetscViewerASCIIPushTab(PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)snes)));CHKERRQ(ierr);} 3767 for (grid=0; grid<snes->gridsequence+1; grid++) { 3768 3769 /* set solution vector */ 3770 if (!grid) {ierr = PetscObjectReference((PetscObject)x);CHKERRQ(ierr);} 3771 ierr = VecDestroy(&snes->vec_sol);CHKERRQ(ierr); 3772 snes->vec_sol = x; 3773 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 3774 3775 /* set affine vector if provided */ 3776 if (b) { ierr = PetscObjectReference((PetscObject)b);CHKERRQ(ierr); } 3777 ierr = VecDestroy(&snes->vec_rhs);CHKERRQ(ierr); 3778 snes->vec_rhs = b; 3779 3780 ierr = SNESSetUp(snes);CHKERRQ(ierr); 3781 3782 if (!grid) { 3783 if (snes->ops->computeinitialguess) { 3784 ierr = (*snes->ops->computeinitialguess)(snes,snes->vec_sol,snes->initialguessP);CHKERRQ(ierr); 3785 } 3786 } 3787 3788 if (snes->conv_hist_reset) snes->conv_hist_len = 0; 3789 if (snes->counters_reset) {snes->nfuncs = 0; snes->linear_its = 0; snes->numFailures = 0;} 3790 3791 ierr = PetscLogEventBegin(SNES_Solve,snes,0,0,0);CHKERRQ(ierr); 3792 ierr = (*snes->ops->solve)(snes);CHKERRQ(ierr); 3793 ierr = PetscLogEventEnd(SNES_Solve,snes,0,0,0);CHKERRQ(ierr); 3794 if (snes->domainerror) { 3795 snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; 3796 snes->domainerror = PETSC_FALSE; 3797 } 3798 if (!snes->reason) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Internal error, solver returned without setting converged reason"); 3799 3800 if (snes->lagjac_persist) snes->jac_iter += snes->iter; 3801 if (snes->lagpre_persist) snes->pre_iter += snes->iter; 3802 3803 flg = PETSC_FALSE; 3804 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_test_local_min",&flg,NULL);CHKERRQ(ierr); 3805 if (flg && !PetscPreLoadingOn) { ierr = SNESTestLocalMin(snes);CHKERRQ(ierr); } 3806 if (snes->printreason) { 3807 ierr = PetscViewerASCIIAddTab(PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)snes)),((PetscObject)snes)->tablevel);CHKERRQ(ierr); 3808 if (snes->reason > 0) { 3809 ierr = PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)snes)),"Nonlinear solve converged due to %s iterations %D\n",SNESConvergedReasons[snes->reason],snes->iter);CHKERRQ(ierr); 3810 } else { 3811 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); 3812 } 3813 ierr = PetscViewerASCIISubtractTab(PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)snes)),((PetscObject)snes)->tablevel);CHKERRQ(ierr); 3814 } 3815 3816 if (snes->errorifnotconverged && snes->reason < 0) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_NOT_CONVERGED,"SNESSolve has not converged"); 3817 if (grid < snes->gridsequence) { 3818 DM fine; 3819 Vec xnew; 3820 Mat interp; 3821 3822 ierr = DMRefine(snes->dm,PetscObjectComm((PetscObject)snes),&fine);CHKERRQ(ierr); 3823 if (!fine) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_INCOMP,"DMRefine() did not perform any refinement, cannot continue grid sequencing"); 3824 ierr = DMCreateInterpolation(snes->dm,fine,&interp,NULL);CHKERRQ(ierr); 3825 ierr = DMCreateGlobalVector(fine,&xnew);CHKERRQ(ierr); 3826 ierr = MatInterpolate(interp,x,xnew);CHKERRQ(ierr); 3827 ierr = DMInterpolate(snes->dm,interp,fine);CHKERRQ(ierr); 3828 ierr = MatDestroy(&interp);CHKERRQ(ierr); 3829 x = xnew; 3830 3831 ierr = SNESReset(snes);CHKERRQ(ierr); 3832 ierr = SNESSetDM(snes,fine);CHKERRQ(ierr); 3833 ierr = DMDestroy(&fine);CHKERRQ(ierr); 3834 ierr = PetscViewerASCIIPopTab(PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)snes)));CHKERRQ(ierr); 3835 } 3836 } 3837 /* monitoring and viewing */ 3838 ierr = PetscOptionsGetViewer(PetscObjectComm((PetscObject)snes),((PetscObject)snes)->prefix,"-snes_view",&viewer,&format,&flg);CHKERRQ(ierr); 3839 if (flg && !PetscPreLoadingOn) { 3840 ierr = PetscViewerPushFormat(viewer,format);CHKERRQ(ierr); 3841 ierr = SNESView(snes,viewer);CHKERRQ(ierr); 3842 ierr = PetscViewerPopFormat(viewer);CHKERRQ(ierr); 3843 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 3844 } 3845 ierr = VecViewFromOptions(snes->vec_sol,((PetscObject)snes)->prefix,"-snes_view_solution");CHKERRQ(ierr); 3846 3847 ierr = VecDestroy(&xcreated);CHKERRQ(ierr); 3848 ierr = PetscObjectAMSBlock((PetscObject)snes);CHKERRQ(ierr); 3849 PetscFunctionReturn(0); 3850 } 3851 3852 /* --------- Internal routines for SNES Package --------- */ 3853 3854 #undef __FUNCT__ 3855 #define __FUNCT__ "SNESSetType" 3856 /*@C 3857 SNESSetType - Sets the method for the nonlinear solver. 3858 3859 Collective on SNES 3860 3861 Input Parameters: 3862 + snes - the SNES context 3863 - type - a known method 3864 3865 Options Database Key: 3866 . -snes_type <type> - Sets the method; use -help for a list 3867 of available methods (for instance, newtonls or newtontr) 3868 3869 Notes: 3870 See "petsc/include/petscsnes.h" for available methods (for instance) 3871 + SNESNEWTONLS - Newton's method with line search 3872 (systems of nonlinear equations) 3873 . SNESNEWTONTR - Newton's method with trust region 3874 (systems of nonlinear equations) 3875 3876 Normally, it is best to use the SNESSetFromOptions() command and then 3877 set the SNES solver type from the options database rather than by using 3878 this routine. Using the options database provides the user with 3879 maximum flexibility in evaluating the many nonlinear solvers. 3880 The SNESSetType() routine is provided for those situations where it 3881 is necessary to set the nonlinear solver independently of the command 3882 line or options database. This might be the case, for example, when 3883 the choice of solver changes during the execution of the program, 3884 and the user's application is taking responsibility for choosing the 3885 appropriate method. 3886 3887 Developer Notes: SNESRegister() adds a constructor for a new SNESType to SNESList, SNESSetType() locates 3888 the constructor in that list and calls it to create the spexific object. 3889 3890 Level: intermediate 3891 3892 .keywords: SNES, set, type 3893 3894 .seealso: SNESType, SNESCreate(), SNESDestroy(), SNESGetType(), SNESSetFromOptions() 3895 3896 @*/ 3897 PetscErrorCode SNESSetType(SNES snes,SNESType type) 3898 { 3899 PetscErrorCode ierr,(*r)(SNES); 3900 PetscBool match; 3901 3902 PetscFunctionBegin; 3903 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3904 PetscValidCharPointer(type,2); 3905 3906 ierr = PetscObjectTypeCompare((PetscObject)snes,type,&match);CHKERRQ(ierr); 3907 if (match) PetscFunctionReturn(0); 3908 3909 ierr = PetscFunctionListFind(SNESList,type,&r);CHKERRQ(ierr); 3910 if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unable to find requested SNES type %s",type); 3911 /* Destroy the previous private SNES context */ 3912 if (snes->ops->destroy) { 3913 ierr = (*(snes)->ops->destroy)(snes);CHKERRQ(ierr); 3914 snes->ops->destroy = NULL; 3915 } 3916 /* Reinitialize function pointers in SNESOps structure */ 3917 snes->ops->setup = 0; 3918 snes->ops->solve = 0; 3919 snes->ops->view = 0; 3920 snes->ops->setfromoptions = 0; 3921 snes->ops->destroy = 0; 3922 /* Call the SNESCreate_XXX routine for this particular Nonlinear solver */ 3923 snes->setupcalled = PETSC_FALSE; 3924 3925 ierr = PetscObjectChangeTypeName((PetscObject)snes,type);CHKERRQ(ierr); 3926 ierr = (*r)(snes);CHKERRQ(ierr); 3927 PetscFunctionReturn(0); 3928 } 3929 3930 #undef __FUNCT__ 3931 #define __FUNCT__ "SNESGetType" 3932 /*@C 3933 SNESGetType - Gets the SNES method type and name (as a string). 3934 3935 Not Collective 3936 3937 Input Parameter: 3938 . snes - nonlinear solver context 3939 3940 Output Parameter: 3941 . type - SNES method (a character string) 3942 3943 Level: intermediate 3944 3945 .keywords: SNES, nonlinear, get, type, name 3946 @*/ 3947 PetscErrorCode SNESGetType(SNES snes,SNESType *type) 3948 { 3949 PetscFunctionBegin; 3950 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3951 PetscValidPointer(type,2); 3952 *type = ((PetscObject)snes)->type_name; 3953 PetscFunctionReturn(0); 3954 } 3955 3956 #undef __FUNCT__ 3957 #define __FUNCT__ "SNESGetSolution" 3958 /*@ 3959 SNESGetSolution - Returns the vector where the approximate solution is 3960 stored. This is the fine grid solution when using SNESSetGridSequence(). 3961 3962 Not Collective, but Vec is parallel if SNES is parallel 3963 3964 Input Parameter: 3965 . snes - the SNES context 3966 3967 Output Parameter: 3968 . x - the solution 3969 3970 Level: intermediate 3971 3972 .keywords: SNES, nonlinear, get, solution 3973 3974 .seealso: SNESGetSolutionUpdate(), SNESGetFunction() 3975 @*/ 3976 PetscErrorCode SNESGetSolution(SNES snes,Vec *x) 3977 { 3978 PetscFunctionBegin; 3979 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3980 PetscValidPointer(x,2); 3981 *x = snes->vec_sol; 3982 PetscFunctionReturn(0); 3983 } 3984 3985 #undef __FUNCT__ 3986 #define __FUNCT__ "SNESGetSolutionUpdate" 3987 /*@ 3988 SNESGetSolutionUpdate - Returns the vector where the solution update is 3989 stored. 3990 3991 Not Collective, but Vec is parallel if SNES is parallel 3992 3993 Input Parameter: 3994 . snes - the SNES context 3995 3996 Output Parameter: 3997 . x - the solution update 3998 3999 Level: advanced 4000 4001 .keywords: SNES, nonlinear, get, solution, update 4002 4003 .seealso: SNESGetSolution(), SNESGetFunction() 4004 @*/ 4005 PetscErrorCode SNESGetSolutionUpdate(SNES snes,Vec *x) 4006 { 4007 PetscFunctionBegin; 4008 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4009 PetscValidPointer(x,2); 4010 *x = snes->vec_sol_update; 4011 PetscFunctionReturn(0); 4012 } 4013 4014 #undef __FUNCT__ 4015 #define __FUNCT__ "SNESGetFunction" 4016 /*@C 4017 SNESGetFunction - Returns the vector where the function is stored. 4018 4019 Not Collective, but Vec is parallel if SNES is parallel. Collective if Vec is requested, but has not been created yet. 4020 4021 Input Parameter: 4022 . snes - the SNES context 4023 4024 Output Parameter: 4025 + r - the vector that is used to store residuals (or NULL if you don't want it) 4026 . SNESFunction- the function (or NULL if you don't want it) 4027 - ctx - the function context (or NULL if you don't want it) 4028 4029 Level: advanced 4030 4031 .keywords: SNES, nonlinear, get, function 4032 4033 .seealso: SNESSetFunction(), SNESGetSolution(), SNESFunction 4034 @*/ 4035 PetscErrorCode SNESGetFunction(SNES snes,Vec *r,PetscErrorCode (**SNESFunction)(SNES,Vec,Vec,void*),void **ctx) 4036 { 4037 PetscErrorCode ierr; 4038 DM dm; 4039 4040 PetscFunctionBegin; 4041 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4042 if (r) { 4043 if (!snes->vec_func) { 4044 if (snes->vec_rhs) { 4045 ierr = VecDuplicate(snes->vec_rhs,&snes->vec_func);CHKERRQ(ierr); 4046 } else if (snes->vec_sol) { 4047 ierr = VecDuplicate(snes->vec_sol,&snes->vec_func);CHKERRQ(ierr); 4048 } else if (snes->dm) { 4049 ierr = DMCreateGlobalVector(snes->dm,&snes->vec_func);CHKERRQ(ierr); 4050 } 4051 } 4052 *r = snes->vec_func; 4053 } 4054 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 4055 ierr = DMSNESGetFunction(dm,SNESFunction,ctx);CHKERRQ(ierr); 4056 PetscFunctionReturn(0); 4057 } 4058 4059 /*@C 4060 SNESGetGS - Returns the GS function and context. 4061 4062 Input Parameter: 4063 . snes - the SNES context 4064 4065 Output Parameter: 4066 + SNESGSFunction - the function (or NULL) 4067 - ctx - the function context (or NULL) 4068 4069 Level: advanced 4070 4071 .keywords: SNES, nonlinear, get, function 4072 4073 .seealso: SNESSetGS(), SNESGetFunction() 4074 @*/ 4075 4076 #undef __FUNCT__ 4077 #define __FUNCT__ "SNESGetGS" 4078 PetscErrorCode SNESGetGS (SNES snes, PetscErrorCode (**SNESGSFunction)(SNES, Vec, Vec, void*), void ** ctx) 4079 { 4080 PetscErrorCode ierr; 4081 DM dm; 4082 4083 PetscFunctionBegin; 4084 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4085 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 4086 ierr = DMSNESGetGS(dm,SNESGSFunction,ctx);CHKERRQ(ierr); 4087 PetscFunctionReturn(0); 4088 } 4089 4090 #undef __FUNCT__ 4091 #define __FUNCT__ "SNESSetOptionsPrefix" 4092 /*@C 4093 SNESSetOptionsPrefix - Sets the prefix used for searching for all 4094 SNES options in the database. 4095 4096 Logically Collective on SNES 4097 4098 Input Parameter: 4099 + snes - the SNES context 4100 - prefix - the prefix to prepend to all option names 4101 4102 Notes: 4103 A hyphen (-) must NOT be given at the beginning of the prefix name. 4104 The first character of all runtime options is AUTOMATICALLY the hyphen. 4105 4106 Level: advanced 4107 4108 .keywords: SNES, set, options, prefix, database 4109 4110 .seealso: SNESSetFromOptions() 4111 @*/ 4112 PetscErrorCode SNESSetOptionsPrefix(SNES snes,const char prefix[]) 4113 { 4114 PetscErrorCode ierr; 4115 4116 PetscFunctionBegin; 4117 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4118 ierr = PetscObjectSetOptionsPrefix((PetscObject)snes,prefix);CHKERRQ(ierr); 4119 if (!snes->ksp) {ierr = SNESGetKSP(snes,&snes->ksp);CHKERRQ(ierr);} 4120 if (snes->linesearch) { 4121 ierr = SNESGetLineSearch(snes,&snes->linesearch);CHKERRQ(ierr); 4122 ierr = PetscObjectSetOptionsPrefix((PetscObject)snes->linesearch,prefix);CHKERRQ(ierr); 4123 } 4124 ierr = KSPSetOptionsPrefix(snes->ksp,prefix);CHKERRQ(ierr); 4125 PetscFunctionReturn(0); 4126 } 4127 4128 #undef __FUNCT__ 4129 #define __FUNCT__ "SNESAppendOptionsPrefix" 4130 /*@C 4131 SNESAppendOptionsPrefix - Appends to the prefix used for searching for all 4132 SNES options in the database. 4133 4134 Logically Collective on SNES 4135 4136 Input Parameters: 4137 + snes - the SNES context 4138 - prefix - the prefix to prepend to all option names 4139 4140 Notes: 4141 A hyphen (-) must NOT be given at the beginning of the prefix name. 4142 The first character of all runtime options is AUTOMATICALLY the hyphen. 4143 4144 Level: advanced 4145 4146 .keywords: SNES, append, options, prefix, database 4147 4148 .seealso: SNESGetOptionsPrefix() 4149 @*/ 4150 PetscErrorCode SNESAppendOptionsPrefix(SNES snes,const char prefix[]) 4151 { 4152 PetscErrorCode ierr; 4153 4154 PetscFunctionBegin; 4155 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4156 ierr = PetscObjectAppendOptionsPrefix((PetscObject)snes,prefix);CHKERRQ(ierr); 4157 if (!snes->ksp) {ierr = SNESGetKSP(snes,&snes->ksp);CHKERRQ(ierr);} 4158 if (snes->linesearch) { 4159 ierr = SNESGetLineSearch(snes,&snes->linesearch);CHKERRQ(ierr); 4160 ierr = PetscObjectAppendOptionsPrefix((PetscObject)snes->linesearch,prefix);CHKERRQ(ierr); 4161 } 4162 ierr = KSPAppendOptionsPrefix(snes->ksp,prefix);CHKERRQ(ierr); 4163 PetscFunctionReturn(0); 4164 } 4165 4166 #undef __FUNCT__ 4167 #define __FUNCT__ "SNESGetOptionsPrefix" 4168 /*@C 4169 SNESGetOptionsPrefix - Sets the prefix used for searching for all 4170 SNES options in the database. 4171 4172 Not Collective 4173 4174 Input Parameter: 4175 . snes - the SNES context 4176 4177 Output Parameter: 4178 . prefix - pointer to the prefix string used 4179 4180 Notes: On the fortran side, the user should pass in a string 'prefix' of 4181 sufficient length to hold the prefix. 4182 4183 Level: advanced 4184 4185 .keywords: SNES, get, options, prefix, database 4186 4187 .seealso: SNESAppendOptionsPrefix() 4188 @*/ 4189 PetscErrorCode SNESGetOptionsPrefix(SNES snes,const char *prefix[]) 4190 { 4191 PetscErrorCode ierr; 4192 4193 PetscFunctionBegin; 4194 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4195 ierr = PetscObjectGetOptionsPrefix((PetscObject)snes,prefix);CHKERRQ(ierr); 4196 PetscFunctionReturn(0); 4197 } 4198 4199 4200 #undef __FUNCT__ 4201 #define __FUNCT__ "SNESRegister" 4202 /*@C 4203 SNESRegister - Adds a method to the nonlinear solver package. 4204 4205 Not collective 4206 4207 Input Parameters: 4208 + name_solver - name of a new user-defined solver 4209 - routine_create - routine to create method context 4210 4211 Notes: 4212 SNESRegister() may be called multiple times to add several user-defined solvers. 4213 4214 Sample usage: 4215 .vb 4216 SNESRegister("my_solver",MySolverCreate); 4217 .ve 4218 4219 Then, your solver can be chosen with the procedural interface via 4220 $ SNESSetType(snes,"my_solver") 4221 or at runtime via the option 4222 $ -snes_type my_solver 4223 4224 Level: advanced 4225 4226 Note: If your function is not being put into a shared library then use SNESRegister() instead 4227 4228 .keywords: SNES, nonlinear, register 4229 4230 .seealso: SNESRegisterAll(), SNESRegisterDestroy() 4231 4232 Level: advanced 4233 @*/ 4234 PetscErrorCode SNESRegister(const char sname[],PetscErrorCode (*function)(SNES)) 4235 { 4236 PetscErrorCode ierr; 4237 4238 PetscFunctionBegin; 4239 ierr = PetscFunctionListAdd(&SNESList,sname,function);CHKERRQ(ierr); 4240 PetscFunctionReturn(0); 4241 } 4242 4243 #undef __FUNCT__ 4244 #define __FUNCT__ "SNESTestLocalMin" 4245 PetscErrorCode SNESTestLocalMin(SNES snes) 4246 { 4247 PetscErrorCode ierr; 4248 PetscInt N,i,j; 4249 Vec u,uh,fh; 4250 PetscScalar value; 4251 PetscReal norm; 4252 4253 PetscFunctionBegin; 4254 ierr = SNESGetSolution(snes,&u);CHKERRQ(ierr); 4255 ierr = VecDuplicate(u,&uh);CHKERRQ(ierr); 4256 ierr = VecDuplicate(u,&fh);CHKERRQ(ierr); 4257 4258 /* currently only works for sequential */ 4259 ierr = PetscPrintf(PETSC_COMM_WORLD,"Testing FormFunction() for local min\n");CHKERRQ(ierr); 4260 ierr = VecGetSize(u,&N);CHKERRQ(ierr); 4261 for (i=0; i<N; i++) { 4262 ierr = VecCopy(u,uh);CHKERRQ(ierr); 4263 ierr = PetscPrintf(PETSC_COMM_WORLD,"i = %D\n",i);CHKERRQ(ierr); 4264 for (j=-10; j<11; j++) { 4265 value = PetscSign(j)*exp(PetscAbs(j)-10.0); 4266 ierr = VecSetValue(uh,i,value,ADD_VALUES);CHKERRQ(ierr); 4267 ierr = SNESComputeFunction(snes,uh,fh);CHKERRQ(ierr); 4268 ierr = VecNorm(fh,NORM_2,&norm);CHKERRQ(ierr); 4269 ierr = PetscPrintf(PETSC_COMM_WORLD," j norm %D %18.16e\n",j,norm);CHKERRQ(ierr); 4270 value = -value; 4271 ierr = VecSetValue(uh,i,value,ADD_VALUES);CHKERRQ(ierr); 4272 } 4273 } 4274 ierr = VecDestroy(&uh);CHKERRQ(ierr); 4275 ierr = VecDestroy(&fh);CHKERRQ(ierr); 4276 PetscFunctionReturn(0); 4277 } 4278 4279 #undef __FUNCT__ 4280 #define __FUNCT__ "SNESKSPSetUseEW" 4281 /*@ 4282 SNESKSPSetUseEW - Sets SNES use Eisenstat-Walker method for 4283 computing relative tolerance for linear solvers within an inexact 4284 Newton method. 4285 4286 Logically Collective on SNES 4287 4288 Input Parameters: 4289 + snes - SNES context 4290 - flag - PETSC_TRUE or PETSC_FALSE 4291 4292 Options Database: 4293 + -snes_ksp_ew - use Eisenstat-Walker method for determining linear system convergence 4294 . -snes_ksp_ew_version ver - version of Eisenstat-Walker method 4295 . -snes_ksp_ew_rtol0 <rtol0> - Sets rtol0 4296 . -snes_ksp_ew_rtolmax <rtolmax> - Sets rtolmax 4297 . -snes_ksp_ew_gamma <gamma> - Sets gamma 4298 . -snes_ksp_ew_alpha <alpha> - Sets alpha 4299 . -snes_ksp_ew_alpha2 <alpha2> - Sets alpha2 4300 - -snes_ksp_ew_threshold <threshold> - Sets threshold 4301 4302 Notes: 4303 Currently, the default is to use a constant relative tolerance for 4304 the inner linear solvers. Alternatively, one can use the 4305 Eisenstat-Walker method, where the relative convergence tolerance 4306 is reset at each Newton iteration according progress of the nonlinear 4307 solver. 4308 4309 Level: advanced 4310 4311 Reference: 4312 S. C. Eisenstat and H. F. Walker, "Choosing the forcing terms in an 4313 inexact Newton method", SISC 17 (1), pp.16-32, 1996. 4314 4315 .keywords: SNES, KSP, Eisenstat, Walker, convergence, test, inexact, Newton 4316 4317 .seealso: SNESKSPGetUseEW(), SNESKSPGetParametersEW(), SNESKSPSetParametersEW() 4318 @*/ 4319 PetscErrorCode SNESKSPSetUseEW(SNES snes,PetscBool flag) 4320 { 4321 PetscFunctionBegin; 4322 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4323 PetscValidLogicalCollectiveBool(snes,flag,2); 4324 snes->ksp_ewconv = flag; 4325 PetscFunctionReturn(0); 4326 } 4327 4328 #undef __FUNCT__ 4329 #define __FUNCT__ "SNESKSPGetUseEW" 4330 /*@ 4331 SNESKSPGetUseEW - Gets if SNES is using Eisenstat-Walker method 4332 for computing relative tolerance for linear solvers within an 4333 inexact Newton method. 4334 4335 Not Collective 4336 4337 Input Parameter: 4338 . snes - SNES context 4339 4340 Output Parameter: 4341 . flag - PETSC_TRUE or PETSC_FALSE 4342 4343 Notes: 4344 Currently, the default is to use a constant relative tolerance for 4345 the inner linear solvers. Alternatively, one can use the 4346 Eisenstat-Walker method, where the relative convergence tolerance 4347 is reset at each Newton iteration according progress of the nonlinear 4348 solver. 4349 4350 Level: advanced 4351 4352 Reference: 4353 S. C. Eisenstat and H. F. Walker, "Choosing the forcing terms in an 4354 inexact Newton method", SISC 17 (1), pp.16-32, 1996. 4355 4356 .keywords: SNES, KSP, Eisenstat, Walker, convergence, test, inexact, Newton 4357 4358 .seealso: SNESKSPSetUseEW(), SNESKSPGetParametersEW(), SNESKSPSetParametersEW() 4359 @*/ 4360 PetscErrorCode SNESKSPGetUseEW(SNES snes, PetscBool *flag) 4361 { 4362 PetscFunctionBegin; 4363 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4364 PetscValidPointer(flag,2); 4365 *flag = snes->ksp_ewconv; 4366 PetscFunctionReturn(0); 4367 } 4368 4369 #undef __FUNCT__ 4370 #define __FUNCT__ "SNESKSPSetParametersEW" 4371 /*@ 4372 SNESKSPSetParametersEW - Sets parameters for Eisenstat-Walker 4373 convergence criteria for the linear solvers within an inexact 4374 Newton method. 4375 4376 Logically Collective on SNES 4377 4378 Input Parameters: 4379 + snes - SNES context 4380 . version - version 1, 2 (default is 2) or 3 4381 . rtol_0 - initial relative tolerance (0 <= rtol_0 < 1) 4382 . rtol_max - maximum relative tolerance (0 <= rtol_max < 1) 4383 . gamma - multiplicative factor for version 2 rtol computation 4384 (0 <= gamma2 <= 1) 4385 . alpha - power for version 2 rtol computation (1 < alpha <= 2) 4386 . alpha2 - power for safeguard 4387 - threshold - threshold for imposing safeguard (0 < threshold < 1) 4388 4389 Note: 4390 Version 3 was contributed by Luis Chacon, June 2006. 4391 4392 Use PETSC_DEFAULT to retain the default for any of the parameters. 4393 4394 Level: advanced 4395 4396 Reference: 4397 S. C. Eisenstat and H. F. Walker, "Choosing the forcing terms in an 4398 inexact Newton method", Utah State University Math. Stat. Dept. Res. 4399 Report 6/94/75, June, 1994, to appear in SIAM J. Sci. Comput. 4400 4401 .keywords: SNES, KSP, Eisenstat, Walker, set, parameters 4402 4403 .seealso: SNESKSPSetUseEW(), SNESKSPGetUseEW(), SNESKSPGetParametersEW() 4404 @*/ 4405 PetscErrorCode SNESKSPSetParametersEW(SNES snes,PetscInt version,PetscReal rtol_0,PetscReal rtol_max,PetscReal gamma,PetscReal alpha,PetscReal alpha2,PetscReal threshold) 4406 { 4407 SNESKSPEW *kctx; 4408 4409 PetscFunctionBegin; 4410 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4411 kctx = (SNESKSPEW*)snes->kspconvctx; 4412 if (!kctx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"No Eisenstat-Walker context existing"); 4413 PetscValidLogicalCollectiveInt(snes,version,2); 4414 PetscValidLogicalCollectiveReal(snes,rtol_0,3); 4415 PetscValidLogicalCollectiveReal(snes,rtol_max,4); 4416 PetscValidLogicalCollectiveReal(snes,gamma,5); 4417 PetscValidLogicalCollectiveReal(snes,alpha,6); 4418 PetscValidLogicalCollectiveReal(snes,alpha2,7); 4419 PetscValidLogicalCollectiveReal(snes,threshold,8); 4420 4421 if (version != PETSC_DEFAULT) kctx->version = version; 4422 if (rtol_0 != PETSC_DEFAULT) kctx->rtol_0 = rtol_0; 4423 if (rtol_max != PETSC_DEFAULT) kctx->rtol_max = rtol_max; 4424 if (gamma != PETSC_DEFAULT) kctx->gamma = gamma; 4425 if (alpha != PETSC_DEFAULT) kctx->alpha = alpha; 4426 if (alpha2 != PETSC_DEFAULT) kctx->alpha2 = alpha2; 4427 if (threshold != PETSC_DEFAULT) kctx->threshold = threshold; 4428 4429 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); 4430 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); 4431 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); 4432 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); 4433 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); 4434 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); 4435 PetscFunctionReturn(0); 4436 } 4437 4438 #undef __FUNCT__ 4439 #define __FUNCT__ "SNESKSPGetParametersEW" 4440 /*@ 4441 SNESKSPGetParametersEW - Gets parameters for Eisenstat-Walker 4442 convergence criteria for the linear solvers within an inexact 4443 Newton method. 4444 4445 Not Collective 4446 4447 Input Parameters: 4448 snes - SNES context 4449 4450 Output Parameters: 4451 + version - version 1, 2 (default is 2) or 3 4452 . rtol_0 - initial relative tolerance (0 <= rtol_0 < 1) 4453 . rtol_max - maximum relative tolerance (0 <= rtol_max < 1) 4454 . gamma - multiplicative factor for version 2 rtol computation (0 <= gamma2 <= 1) 4455 . alpha - power for version 2 rtol computation (1 < alpha <= 2) 4456 . alpha2 - power for safeguard 4457 - threshold - threshold for imposing safeguard (0 < threshold < 1) 4458 4459 Level: advanced 4460 4461 .keywords: SNES, KSP, Eisenstat, Walker, get, parameters 4462 4463 .seealso: SNESKSPSetUseEW(), SNESKSPGetUseEW(), SNESKSPSetParametersEW() 4464 @*/ 4465 PetscErrorCode SNESKSPGetParametersEW(SNES snes,PetscInt *version,PetscReal *rtol_0,PetscReal *rtol_max,PetscReal *gamma,PetscReal *alpha,PetscReal *alpha2,PetscReal *threshold) 4466 { 4467 SNESKSPEW *kctx; 4468 4469 PetscFunctionBegin; 4470 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4471 kctx = (SNESKSPEW*)snes->kspconvctx; 4472 if (!kctx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"No Eisenstat-Walker context existing"); 4473 if (version) *version = kctx->version; 4474 if (rtol_0) *rtol_0 = kctx->rtol_0; 4475 if (rtol_max) *rtol_max = kctx->rtol_max; 4476 if (gamma) *gamma = kctx->gamma; 4477 if (alpha) *alpha = kctx->alpha; 4478 if (alpha2) *alpha2 = kctx->alpha2; 4479 if (threshold) *threshold = kctx->threshold; 4480 PetscFunctionReturn(0); 4481 } 4482 4483 #undef __FUNCT__ 4484 #define __FUNCT__ "KSPPreSolve_SNESEW" 4485 PetscErrorCode KSPPreSolve_SNESEW(KSP ksp, Vec b, Vec x, SNES snes) 4486 { 4487 PetscErrorCode ierr; 4488 SNESKSPEW *kctx = (SNESKSPEW*)snes->kspconvctx; 4489 PetscReal rtol = PETSC_DEFAULT,stol; 4490 4491 PetscFunctionBegin; 4492 if (!snes->ksp_ewconv) PetscFunctionReturn(0); 4493 if (!snes->iter) { 4494 rtol = kctx->rtol_0; /* first time in, so use the original user rtol */ 4495 ierr = VecNorm(snes->vec_func,NORM_2,&kctx->norm_first);CHKERRQ(ierr); 4496 } 4497 else { 4498 if (kctx->version == 1) { 4499 rtol = (snes->norm - kctx->lresid_last)/kctx->norm_last; 4500 if (rtol < 0.0) rtol = -rtol; 4501 stol = PetscPowReal(kctx->rtol_last,kctx->alpha2); 4502 if (stol > kctx->threshold) rtol = PetscMax(rtol,stol); 4503 } else if (kctx->version == 2) { 4504 rtol = kctx->gamma * PetscPowReal(snes->norm/kctx->norm_last,kctx->alpha); 4505 stol = kctx->gamma * PetscPowReal(kctx->rtol_last,kctx->alpha); 4506 if (stol > kctx->threshold) rtol = PetscMax(rtol,stol); 4507 } else if (kctx->version == 3) { /* contributed by Luis Chacon, June 2006. */ 4508 rtol = kctx->gamma * PetscPowReal(snes->norm/kctx->norm_last,kctx->alpha); 4509 /* safeguard: avoid sharp decrease of rtol */ 4510 stol = kctx->gamma*PetscPowReal(kctx->rtol_last,kctx->alpha); 4511 stol = PetscMax(rtol,stol); 4512 rtol = PetscMin(kctx->rtol_0,stol); 4513 /* safeguard: avoid oversolving */ 4514 stol = kctx->gamma*(kctx->norm_first*snes->rtol)/snes->norm; 4515 stol = PetscMax(rtol,stol); 4516 rtol = PetscMin(kctx->rtol_0,stol); 4517 } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Only versions 1, 2 or 3 are supported: %D",kctx->version); 4518 } 4519 /* safeguard: avoid rtol greater than one */ 4520 rtol = PetscMin(rtol,kctx->rtol_max); 4521 ierr = KSPSetTolerances(ksp,rtol,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);CHKERRQ(ierr); 4522 ierr = PetscInfo3(snes,"iter %D, Eisenstat-Walker (version %D) KSP rtol=%G\n",snes->iter,kctx->version,rtol);CHKERRQ(ierr); 4523 PetscFunctionReturn(0); 4524 } 4525 4526 #undef __FUNCT__ 4527 #define __FUNCT__ "KSPPostSolve_SNESEW" 4528 PetscErrorCode KSPPostSolve_SNESEW(KSP ksp, Vec b, Vec x, SNES snes) 4529 { 4530 PetscErrorCode ierr; 4531 SNESKSPEW *kctx = (SNESKSPEW*)snes->kspconvctx; 4532 PCSide pcside; 4533 Vec lres; 4534 4535 PetscFunctionBegin; 4536 if (!snes->ksp_ewconv) PetscFunctionReturn(0); 4537 ierr = KSPGetTolerances(ksp,&kctx->rtol_last,0,0,0);CHKERRQ(ierr); 4538 ierr = SNESGetFunctionNorm(snes,&kctx->norm_last);CHKERRQ(ierr); 4539 if (kctx->version == 1) { 4540 ierr = KSPGetPCSide(ksp,&pcside);CHKERRQ(ierr); 4541 if (pcside == PC_RIGHT) { /* XXX Should we also test KSP_UNPRECONDITIONED_NORM ? */ 4542 /* KSP residual is true linear residual */ 4543 ierr = KSPGetResidualNorm(ksp,&kctx->lresid_last);CHKERRQ(ierr); 4544 } else { 4545 /* KSP residual is preconditioned residual */ 4546 /* compute true linear residual norm */ 4547 ierr = VecDuplicate(b,&lres);CHKERRQ(ierr); 4548 ierr = MatMult(snes->jacobian,x,lres);CHKERRQ(ierr); 4549 ierr = VecAYPX(lres,-1.0,b);CHKERRQ(ierr); 4550 ierr = VecNorm(lres,NORM_2,&kctx->lresid_last);CHKERRQ(ierr); 4551 ierr = VecDestroy(&lres);CHKERRQ(ierr); 4552 } 4553 } 4554 PetscFunctionReturn(0); 4555 } 4556 4557 #undef __FUNCT__ 4558 #define __FUNCT__ "SNESGetKSP" 4559 /*@ 4560 SNESGetKSP - Returns the KSP context for a SNES solver. 4561 4562 Not Collective, but if SNES object is parallel, then KSP object is parallel 4563 4564 Input Parameter: 4565 . snes - the SNES context 4566 4567 Output Parameter: 4568 . ksp - the KSP context 4569 4570 Notes: 4571 The user can then directly manipulate the KSP context to set various 4572 options, etc. Likewise, the user can then extract and manipulate the 4573 PC contexts as well. 4574 4575 Level: beginner 4576 4577 .keywords: SNES, nonlinear, get, KSP, context 4578 4579 .seealso: KSPGetPC(), SNESCreate(), KSPCreate(), SNESSetKSP() 4580 @*/ 4581 PetscErrorCode SNESGetKSP(SNES snes,KSP *ksp) 4582 { 4583 PetscErrorCode ierr; 4584 4585 PetscFunctionBegin; 4586 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4587 PetscValidPointer(ksp,2); 4588 4589 if (!snes->ksp) { 4590 ierr = KSPCreate(PetscObjectComm((PetscObject)snes),&snes->ksp);CHKERRQ(ierr); 4591 ierr = PetscObjectIncrementTabLevel((PetscObject)snes->ksp,(PetscObject)snes,1);CHKERRQ(ierr); 4592 ierr = PetscLogObjectParent((PetscObject)snes,(PetscObject)snes->ksp);CHKERRQ(ierr); 4593 4594 ierr = KSPSetPreSolve(snes->ksp,(PetscErrorCode (*)(KSP,Vec,Vec,void*))KSPPreSolve_SNESEW,snes);CHKERRQ(ierr); 4595 ierr = KSPSetPostSolve(snes->ksp,(PetscErrorCode (*)(KSP,Vec,Vec,void*))KSPPostSolve_SNESEW,snes);CHKERRQ(ierr); 4596 } 4597 *ksp = snes->ksp; 4598 PetscFunctionReturn(0); 4599 } 4600 4601 4602 #include <petsc-private/dmimpl.h> 4603 #undef __FUNCT__ 4604 #define __FUNCT__ "SNESSetDM" 4605 /*@ 4606 SNESSetDM - Sets the DM that may be used by some preconditioners 4607 4608 Logically Collective on SNES 4609 4610 Input Parameters: 4611 + snes - the preconditioner context 4612 - dm - the dm 4613 4614 Level: intermediate 4615 4616 .seealso: SNESGetDM(), KSPSetDM(), KSPGetDM() 4617 @*/ 4618 PetscErrorCode SNESSetDM(SNES snes,DM dm) 4619 { 4620 PetscErrorCode ierr; 4621 KSP ksp; 4622 DMSNES sdm; 4623 4624 PetscFunctionBegin; 4625 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4626 if (dm) {ierr = PetscObjectReference((PetscObject)dm);CHKERRQ(ierr);} 4627 if (snes->dm) { /* Move the DMSNES context over to the new DM unless the new DM already has one */ 4628 if (snes->dm->dmsnes && snes->dmAuto && !dm->dmsnes) { 4629 ierr = DMCopyDMSNES(snes->dm,dm);CHKERRQ(ierr); 4630 ierr = DMGetDMSNES(snes->dm,&sdm);CHKERRQ(ierr); 4631 if (sdm->originaldm == snes->dm) sdm->originaldm = dm; /* Grant write privileges to the replacement DM */ 4632 } 4633 ierr = DMDestroy(&snes->dm);CHKERRQ(ierr); 4634 } 4635 snes->dm = dm; 4636 snes->dmAuto = PETSC_FALSE; 4637 4638 ierr = SNESGetKSP(snes,&ksp);CHKERRQ(ierr); 4639 ierr = KSPSetDM(ksp,dm);CHKERRQ(ierr); 4640 ierr = KSPSetDMActive(ksp,PETSC_FALSE);CHKERRQ(ierr); 4641 if (snes->pc) { 4642 ierr = SNESSetDM(snes->pc, snes->dm);CHKERRQ(ierr); 4643 ierr = SNESSetPCSide(snes,snes->pcside);CHKERRQ(ierr); 4644 } 4645 PetscFunctionReturn(0); 4646 } 4647 4648 #undef __FUNCT__ 4649 #define __FUNCT__ "SNESGetDM" 4650 /*@ 4651 SNESGetDM - Gets the DM that may be used by some preconditioners 4652 4653 Not Collective but DM obtained is parallel on SNES 4654 4655 Input Parameter: 4656 . snes - the preconditioner context 4657 4658 Output Parameter: 4659 . dm - the dm 4660 4661 Level: intermediate 4662 4663 .seealso: SNESSetDM(), KSPSetDM(), KSPGetDM() 4664 @*/ 4665 PetscErrorCode SNESGetDM(SNES snes,DM *dm) 4666 { 4667 PetscErrorCode ierr; 4668 4669 PetscFunctionBegin; 4670 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4671 if (!snes->dm) { 4672 ierr = DMShellCreate(PetscObjectComm((PetscObject)snes),&snes->dm);CHKERRQ(ierr); 4673 snes->dmAuto = PETSC_TRUE; 4674 } 4675 *dm = snes->dm; 4676 PetscFunctionReturn(0); 4677 } 4678 4679 #undef __FUNCT__ 4680 #define __FUNCT__ "SNESSetPC" 4681 /*@ 4682 SNESSetPC - Sets the nonlinear preconditioner to be used. 4683 4684 Collective on SNES 4685 4686 Input Parameters: 4687 + snes - iterative context obtained from SNESCreate() 4688 - pc - the preconditioner object 4689 4690 Notes: 4691 Use SNESGetPC() to retrieve the preconditioner context (for example, 4692 to configure it using the API). 4693 4694 Level: developer 4695 4696 .keywords: SNES, set, precondition 4697 .seealso: SNESGetPC() 4698 @*/ 4699 PetscErrorCode SNESSetPC(SNES snes, SNES pc) 4700 { 4701 PetscErrorCode ierr; 4702 4703 PetscFunctionBegin; 4704 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); 4705 PetscValidHeaderSpecific(pc, SNES_CLASSID, 2); 4706 PetscCheckSameComm(snes, 1, pc, 2); 4707 ierr = PetscObjectReference((PetscObject) pc);CHKERRQ(ierr); 4708 ierr = SNESDestroy(&snes->pc);CHKERRQ(ierr); 4709 snes->pc = pc; 4710 ierr = PetscLogObjectParent((PetscObject)snes, (PetscObject)snes->pc);CHKERRQ(ierr); 4711 PetscFunctionReturn(0); 4712 } 4713 4714 #undef __FUNCT__ 4715 #define __FUNCT__ "SNESGetPC" 4716 /*@ 4717 SNESGetPC - Returns a pointer to the nonlinear preconditioning context set with SNESSetPC(). 4718 4719 Not Collective 4720 4721 Input Parameter: 4722 . snes - iterative context obtained from SNESCreate() 4723 4724 Output Parameter: 4725 . pc - preconditioner context 4726 4727 Level: developer 4728 4729 .keywords: SNES, get, preconditioner 4730 .seealso: SNESSetPC() 4731 @*/ 4732 PetscErrorCode SNESGetPC(SNES snes, SNES *pc) 4733 { 4734 PetscErrorCode ierr; 4735 const char *optionsprefix; 4736 4737 PetscFunctionBegin; 4738 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); 4739 PetscValidPointer(pc, 2); 4740 if (!snes->pc) { 4741 ierr = SNESCreate(PetscObjectComm((PetscObject)snes),&snes->pc);CHKERRQ(ierr); 4742 ierr = PetscObjectIncrementTabLevel((PetscObject)snes->pc,(PetscObject)snes,1);CHKERRQ(ierr); 4743 ierr = PetscLogObjectParent((PetscObject)snes,(PetscObject)snes->pc);CHKERRQ(ierr); 4744 ierr = SNESGetOptionsPrefix(snes,&optionsprefix);CHKERRQ(ierr); 4745 ierr = SNESSetOptionsPrefix(snes->pc,optionsprefix);CHKERRQ(ierr); 4746 ierr = SNESAppendOptionsPrefix(snes->pc,"npc_");CHKERRQ(ierr); 4747 ierr = SNESSetCountersReset(snes->pc,PETSC_FALSE);CHKERRQ(ierr); 4748 } 4749 *pc = snes->pc; 4750 PetscFunctionReturn(0); 4751 } 4752 4753 #undef __FUNCT__ 4754 #define __FUNCT__ "SNESSetPCSide" 4755 /*@ 4756 SNESSetPCSide - Sets the preconditioning side. 4757 4758 Logically Collective on SNES 4759 4760 Input Parameter: 4761 . snes - iterative context obtained from SNESCreate() 4762 4763 Output Parameter: 4764 . side - the preconditioning side, where side is one of 4765 .vb 4766 PC_LEFT - left preconditioning (default) 4767 PC_RIGHT - right preconditioning 4768 .ve 4769 4770 Options Database Keys: 4771 . -snes_pc_side <right,left> 4772 4773 Level: intermediate 4774 4775 .keywords: SNES, set, right, left, side, preconditioner, flag 4776 4777 .seealso: SNESGetPCSide(), KSPSetPCSide() 4778 @*/ 4779 PetscErrorCode SNESSetPCSide(SNES snes,PCSide side) 4780 { 4781 PetscFunctionBegin; 4782 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4783 PetscValidLogicalCollectiveEnum(snes,side,2); 4784 snes->pcside = side; 4785 PetscFunctionReturn(0); 4786 } 4787 4788 #undef __FUNCT__ 4789 #define __FUNCT__ "SNESGetPCSide" 4790 /*@ 4791 SNESGetPCSide - Gets the preconditioning side. 4792 4793 Not Collective 4794 4795 Input Parameter: 4796 . snes - iterative context obtained from SNESCreate() 4797 4798 Output Parameter: 4799 . side - the preconditioning side, where side is one of 4800 .vb 4801 PC_LEFT - left preconditioning (default) 4802 PC_RIGHT - right preconditioning 4803 .ve 4804 4805 Level: intermediate 4806 4807 .keywords: SNES, get, right, left, side, preconditioner, flag 4808 4809 .seealso: SNESSetPCSide(), KSPGetPCSide() 4810 @*/ 4811 PetscErrorCode SNESGetPCSide(SNES snes,PCSide *side) 4812 { 4813 PetscFunctionBegin; 4814 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4815 PetscValidPointer(side,2); 4816 *side = snes->pcside; 4817 PetscFunctionReturn(0); 4818 } 4819 4820 #undef __FUNCT__ 4821 #define __FUNCT__ "SNESSetLineSearch" 4822 /*@ 4823 SNESSetLineSearch - Sets the linesearch on the SNES instance. 4824 4825 Collective on SNES 4826 4827 Input Parameters: 4828 + snes - iterative context obtained from SNESCreate() 4829 - linesearch - the linesearch object 4830 4831 Notes: 4832 Use SNESGetLineSearch() to retrieve the preconditioner context (for example, 4833 to configure it using the API). 4834 4835 Level: developer 4836 4837 .keywords: SNES, set, linesearch 4838 .seealso: SNESGetLineSearch() 4839 @*/ 4840 PetscErrorCode SNESSetLineSearch(SNES snes, SNESLineSearch linesearch) 4841 { 4842 PetscErrorCode ierr; 4843 4844 PetscFunctionBegin; 4845 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); 4846 PetscValidHeaderSpecific(linesearch, SNESLINESEARCH_CLASSID, 2); 4847 PetscCheckSameComm(snes, 1, linesearch, 2); 4848 ierr = PetscObjectReference((PetscObject) linesearch);CHKERRQ(ierr); 4849 ierr = SNESLineSearchDestroy(&snes->linesearch);CHKERRQ(ierr); 4850 4851 snes->linesearch = linesearch; 4852 4853 ierr = PetscLogObjectParent((PetscObject)snes, (PetscObject)snes->linesearch);CHKERRQ(ierr); 4854 PetscFunctionReturn(0); 4855 } 4856 4857 #undef __FUNCT__ 4858 #define __FUNCT__ "SNESGetLineSearch" 4859 /*@ 4860 SNESGetLineSearch - Returns a pointer to the line search context set with SNESSetLineSearch() 4861 or creates a default line search instance associated with the SNES and returns it. 4862 4863 Not Collective 4864 4865 Input Parameter: 4866 . snes - iterative context obtained from SNESCreate() 4867 4868 Output Parameter: 4869 . linesearch - linesearch context 4870 4871 Level: developer 4872 4873 .keywords: SNES, get, linesearch 4874 .seealso: SNESSetLineSearch() 4875 @*/ 4876 PetscErrorCode SNESGetLineSearch(SNES snes, SNESLineSearch *linesearch) 4877 { 4878 PetscErrorCode ierr; 4879 const char *optionsprefix; 4880 4881 PetscFunctionBegin; 4882 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); 4883 PetscValidPointer(linesearch, 2); 4884 if (!snes->linesearch) { 4885 ierr = SNESGetOptionsPrefix(snes, &optionsprefix);CHKERRQ(ierr); 4886 ierr = SNESLineSearchCreate(PetscObjectComm((PetscObject)snes), &snes->linesearch);CHKERRQ(ierr); 4887 ierr = SNESLineSearchSetSNES(snes->linesearch, snes);CHKERRQ(ierr); 4888 ierr = SNESLineSearchAppendOptionsPrefix(snes->linesearch, optionsprefix);CHKERRQ(ierr); 4889 ierr = PetscObjectIncrementTabLevel((PetscObject) snes->linesearch, (PetscObject) snes, 1);CHKERRQ(ierr); 4890 ierr = PetscLogObjectParent((PetscObject)snes, (PetscObject)snes->linesearch);CHKERRQ(ierr); 4891 } 4892 *linesearch = snes->linesearch; 4893 PetscFunctionReturn(0); 4894 } 4895 4896 #if defined(PETSC_HAVE_MATLAB_ENGINE) 4897 #include <mex.h> 4898 4899 typedef struct {char *funcname; mxArray *ctx;} SNESMatlabContext; 4900 4901 #undef __FUNCT__ 4902 #define __FUNCT__ "SNESComputeFunction_Matlab" 4903 /* 4904 SNESComputeFunction_Matlab - Calls the function that has been set with SNESSetFunctionMatlab(). 4905 4906 Collective on SNES 4907 4908 Input Parameters: 4909 + snes - the SNES context 4910 - x - input vector 4911 4912 Output Parameter: 4913 . y - function vector, as set by SNESSetFunction() 4914 4915 Notes: 4916 SNESComputeFunction() is typically used within nonlinear solvers 4917 implementations, so most users would not generally call this routine 4918 themselves. 4919 4920 Level: developer 4921 4922 .keywords: SNES, nonlinear, compute, function 4923 4924 .seealso: SNESSetFunction(), SNESGetFunction() 4925 */ 4926 PetscErrorCode SNESComputeFunction_Matlab(SNES snes,Vec x,Vec y, void *ctx) 4927 { 4928 PetscErrorCode ierr; 4929 SNESMatlabContext *sctx = (SNESMatlabContext*)ctx; 4930 int nlhs = 1,nrhs = 5; 4931 mxArray *plhs[1],*prhs[5]; 4932 long long int lx = 0,ly = 0,ls = 0; 4933 4934 PetscFunctionBegin; 4935 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4936 PetscValidHeaderSpecific(x,VEC_CLASSID,2); 4937 PetscValidHeaderSpecific(y,VEC_CLASSID,3); 4938 PetscCheckSameComm(snes,1,x,2); 4939 PetscCheckSameComm(snes,1,y,3); 4940 4941 /* call Matlab function in ctx with arguments x and y */ 4942 4943 ierr = PetscMemcpy(&ls,&snes,sizeof(snes));CHKERRQ(ierr); 4944 ierr = PetscMemcpy(&lx,&x,sizeof(x));CHKERRQ(ierr); 4945 ierr = PetscMemcpy(&ly,&y,sizeof(x));CHKERRQ(ierr); 4946 prhs[0] = mxCreateDoubleScalar((double)ls); 4947 prhs[1] = mxCreateDoubleScalar((double)lx); 4948 prhs[2] = mxCreateDoubleScalar((double)ly); 4949 prhs[3] = mxCreateString(sctx->funcname); 4950 prhs[4] = sctx->ctx; 4951 ierr = mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscSNESComputeFunctionInternal");CHKERRQ(ierr); 4952 ierr = mxGetScalar(plhs[0]);CHKERRQ(ierr); 4953 mxDestroyArray(prhs[0]); 4954 mxDestroyArray(prhs[1]); 4955 mxDestroyArray(prhs[2]); 4956 mxDestroyArray(prhs[3]); 4957 mxDestroyArray(plhs[0]); 4958 PetscFunctionReturn(0); 4959 } 4960 4961 #undef __FUNCT__ 4962 #define __FUNCT__ "SNESSetFunctionMatlab" 4963 /* 4964 SNESSetFunctionMatlab - Sets the function evaluation routine and function 4965 vector for use by the SNES routines in solving systems of nonlinear 4966 equations from MATLAB. Here the function is a string containing the name of a MATLAB function 4967 4968 Logically Collective on SNES 4969 4970 Input Parameters: 4971 + snes - the SNES context 4972 . r - vector to store function value 4973 - SNESFunction - function evaluation routine 4974 4975 Notes: 4976 The Newton-like methods typically solve linear systems of the form 4977 $ f'(x) x = -f(x), 4978 where f'(x) denotes the Jacobian matrix and f(x) is the function. 4979 4980 Level: beginner 4981 4982 Developer Note: This bleeds the allocated memory SNESMatlabContext *sctx; 4983 4984 .keywords: SNES, nonlinear, set, function 4985 4986 .seealso: SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESSetFunction() 4987 */ 4988 PetscErrorCode SNESSetFunctionMatlab(SNES snes,Vec r,const char *SNESFunction,mxArray *ctx) 4989 { 4990 PetscErrorCode ierr; 4991 SNESMatlabContext *sctx; 4992 4993 PetscFunctionBegin; 4994 /* currently sctx is memory bleed */ 4995 ierr = PetscMalloc(sizeof(SNESMatlabContext),&sctx);CHKERRQ(ierr); 4996 ierr = PetscStrallocpy(SNESFunction,&sctx->funcname);CHKERRQ(ierr); 4997 /* 4998 This should work, but it doesn't 4999 sctx->ctx = ctx; 5000 mexMakeArrayPersistent(sctx->ctx); 5001 */ 5002 sctx->ctx = mxDuplicateArray(ctx); 5003 ierr = SNESSetFunction(snes,r,SNESComputeFunction_Matlab,sctx);CHKERRQ(ierr); 5004 PetscFunctionReturn(0); 5005 } 5006 5007 #undef __FUNCT__ 5008 #define __FUNCT__ "SNESComputeJacobian_Matlab" 5009 /* 5010 SNESComputeJacobian_Matlab - Calls the function that has been set with SNESSetJacobianMatlab(). 5011 5012 Collective on SNES 5013 5014 Input Parameters: 5015 + snes - the SNES context 5016 . x - input vector 5017 . A, B - the matrices 5018 - ctx - user context 5019 5020 Output Parameter: 5021 . flag - structure of the matrix 5022 5023 Level: developer 5024 5025 .keywords: SNES, nonlinear, compute, function 5026 5027 .seealso: SNESSetFunction(), SNESGetFunction() 5028 @*/ 5029 PetscErrorCode SNESComputeJacobian_Matlab(SNES snes,Vec x,Mat *A,Mat *B,MatStructure *flag, void *ctx) 5030 { 5031 PetscErrorCode ierr; 5032 SNESMatlabContext *sctx = (SNESMatlabContext*)ctx; 5033 int nlhs = 2,nrhs = 6; 5034 mxArray *plhs[2],*prhs[6]; 5035 long long int lx = 0,lA = 0,ls = 0, lB = 0; 5036 5037 PetscFunctionBegin; 5038 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 5039 PetscValidHeaderSpecific(x,VEC_CLASSID,2); 5040 5041 /* call Matlab function in ctx with arguments x and y */ 5042 5043 ierr = PetscMemcpy(&ls,&snes,sizeof(snes));CHKERRQ(ierr); 5044 ierr = PetscMemcpy(&lx,&x,sizeof(x));CHKERRQ(ierr); 5045 ierr = PetscMemcpy(&lA,A,sizeof(x));CHKERRQ(ierr); 5046 ierr = PetscMemcpy(&lB,B,sizeof(x));CHKERRQ(ierr); 5047 prhs[0] = mxCreateDoubleScalar((double)ls); 5048 prhs[1] = mxCreateDoubleScalar((double)lx); 5049 prhs[2] = mxCreateDoubleScalar((double)lA); 5050 prhs[3] = mxCreateDoubleScalar((double)lB); 5051 prhs[4] = mxCreateString(sctx->funcname); 5052 prhs[5] = sctx->ctx; 5053 ierr = mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscSNESComputeJacobianInternal");CHKERRQ(ierr); 5054 ierr = mxGetScalar(plhs[0]);CHKERRQ(ierr); 5055 *flag = (MatStructure) mxGetScalar(plhs[1]);CHKERRQ(ierr); 5056 mxDestroyArray(prhs[0]); 5057 mxDestroyArray(prhs[1]); 5058 mxDestroyArray(prhs[2]); 5059 mxDestroyArray(prhs[3]); 5060 mxDestroyArray(prhs[4]); 5061 mxDestroyArray(plhs[0]); 5062 mxDestroyArray(plhs[1]); 5063 PetscFunctionReturn(0); 5064 } 5065 5066 #undef __FUNCT__ 5067 #define __FUNCT__ "SNESSetJacobianMatlab" 5068 /* 5069 SNESSetJacobianMatlab - Sets the Jacobian function evaluation routine and two empty Jacobian matrices 5070 vector for use by the SNES routines in solving systems of nonlinear 5071 equations from MATLAB. Here the function is a string containing the name of a MATLAB function 5072 5073 Logically Collective on SNES 5074 5075 Input Parameters: 5076 + snes - the SNES context 5077 . A,B - Jacobian matrices 5078 . SNESJacobianFunction - function evaluation routine 5079 - ctx - user context 5080 5081 Level: developer 5082 5083 Developer Note: This bleeds the allocated memory SNESMatlabContext *sctx; 5084 5085 .keywords: SNES, nonlinear, set, function 5086 5087 .seealso: SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESSetFunction(), SNESJacobianFunction 5088 */ 5089 PetscErrorCode SNESSetJacobianMatlab(SNES snes,Mat A,Mat B,const char *SNESJacobianFunction,mxArray *ctx) 5090 { 5091 PetscErrorCode ierr; 5092 SNESMatlabContext *sctx; 5093 5094 PetscFunctionBegin; 5095 /* currently sctx is memory bleed */ 5096 ierr = PetscMalloc(sizeof(SNESMatlabContext),&sctx);CHKERRQ(ierr); 5097 ierr = PetscStrallocpy(SNESJacobianFunction,&sctx->funcname);CHKERRQ(ierr); 5098 /* 5099 This should work, but it doesn't 5100 sctx->ctx = ctx; 5101 mexMakeArrayPersistent(sctx->ctx); 5102 */ 5103 sctx->ctx = mxDuplicateArray(ctx); 5104 ierr = SNESSetJacobian(snes,A,B,SNESComputeJacobian_Matlab,sctx);CHKERRQ(ierr); 5105 PetscFunctionReturn(0); 5106 } 5107 5108 #undef __FUNCT__ 5109 #define __FUNCT__ "SNESMonitor_Matlab" 5110 /* 5111 SNESMonitor_Matlab - Calls the function that has been set with SNESMonitorSetMatlab(). 5112 5113 Collective on SNES 5114 5115 .seealso: SNESSetFunction(), SNESGetFunction() 5116 @*/ 5117 PetscErrorCode SNESMonitor_Matlab(SNES snes,PetscInt it, PetscReal fnorm, void *ctx) 5118 { 5119 PetscErrorCode ierr; 5120 SNESMatlabContext *sctx = (SNESMatlabContext*)ctx; 5121 int nlhs = 1,nrhs = 6; 5122 mxArray *plhs[1],*prhs[6]; 5123 long long int lx = 0,ls = 0; 5124 Vec x = snes->vec_sol; 5125 5126 PetscFunctionBegin; 5127 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 5128 5129 ierr = PetscMemcpy(&ls,&snes,sizeof(snes));CHKERRQ(ierr); 5130 ierr = PetscMemcpy(&lx,&x,sizeof(x));CHKERRQ(ierr); 5131 prhs[0] = mxCreateDoubleScalar((double)ls); 5132 prhs[1] = mxCreateDoubleScalar((double)it); 5133 prhs[2] = mxCreateDoubleScalar((double)fnorm); 5134 prhs[3] = mxCreateDoubleScalar((double)lx); 5135 prhs[4] = mxCreateString(sctx->funcname); 5136 prhs[5] = sctx->ctx; 5137 ierr = mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscSNESMonitorInternal");CHKERRQ(ierr); 5138 ierr = mxGetScalar(plhs[0]);CHKERRQ(ierr); 5139 mxDestroyArray(prhs[0]); 5140 mxDestroyArray(prhs[1]); 5141 mxDestroyArray(prhs[2]); 5142 mxDestroyArray(prhs[3]); 5143 mxDestroyArray(prhs[4]); 5144 mxDestroyArray(plhs[0]); 5145 PetscFunctionReturn(0); 5146 } 5147 5148 #undef __FUNCT__ 5149 #define __FUNCT__ "SNESMonitorSetMatlab" 5150 /* 5151 SNESMonitorSetMatlab - Sets the monitor function from MATLAB 5152 5153 Level: developer 5154 5155 Developer Note: This bleeds the allocated memory SNESMatlabContext *sctx; 5156 5157 .keywords: SNES, nonlinear, set, function 5158 5159 .seealso: SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESSetFunction() 5160 */ 5161 PetscErrorCode SNESMonitorSetMatlab(SNES snes,const char *SNESMonitorFunction,mxArray *ctx) 5162 { 5163 PetscErrorCode ierr; 5164 SNESMatlabContext *sctx; 5165 5166 PetscFunctionBegin; 5167 /* currently sctx is memory bleed */ 5168 ierr = PetscMalloc(sizeof(SNESMatlabContext),&sctx);CHKERRQ(ierr); 5169 ierr = PetscStrallocpy(SNESMonitorFunction,&sctx->funcname);CHKERRQ(ierr); 5170 /* 5171 This should work, but it doesn't 5172 sctx->ctx = ctx; 5173 mexMakeArrayPersistent(sctx->ctx); 5174 */ 5175 sctx->ctx = mxDuplicateArray(ctx); 5176 ierr = SNESMonitorSet(snes,SNESMonitor_Matlab,sctx,NULL);CHKERRQ(ierr); 5177 PetscFunctionReturn(0); 5178 } 5179 5180 #endif 5181