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