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