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