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