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 Level: intermediate 3729 3730 .keywords: SNES, set, type 3731 3732 .seealso: SNESType, SNESCreate() 3733 3734 @*/ 3735 PetscErrorCode SNESSetType(SNES snes,SNESType type) 3736 { 3737 PetscErrorCode ierr,(*r)(SNES); 3738 PetscBool match; 3739 3740 PetscFunctionBegin; 3741 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3742 PetscValidCharPointer(type,2); 3743 3744 ierr = PetscObjectTypeCompare((PetscObject)snes,type,&match);CHKERRQ(ierr); 3745 if (match) PetscFunctionReturn(0); 3746 3747 ierr = PetscFunctionListFind(SNESList,type,&r);CHKERRQ(ierr); 3748 if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unable to find requested SNES type %s",type); 3749 /* Destroy the previous private SNES context */ 3750 if (snes->ops->destroy) { 3751 ierr = (*(snes)->ops->destroy)(snes);CHKERRQ(ierr); 3752 snes->ops->destroy = NULL; 3753 } 3754 /* Reinitialize function pointers in SNESOps structure */ 3755 snes->ops->setup = 0; 3756 snes->ops->solve = 0; 3757 snes->ops->view = 0; 3758 snes->ops->setfromoptions = 0; 3759 snes->ops->destroy = 0; 3760 /* Call the SNESCreate_XXX routine for this particular Nonlinear solver */ 3761 snes->setupcalled = PETSC_FALSE; 3762 3763 ierr = PetscObjectChangeTypeName((PetscObject)snes,type);CHKERRQ(ierr); 3764 ierr = (*r)(snes);CHKERRQ(ierr); 3765 PetscFunctionReturn(0); 3766 } 3767 3768 #undef __FUNCT__ 3769 #define __FUNCT__ "SNESGetType" 3770 /*@C 3771 SNESGetType - Gets the SNES method type and name (as a string). 3772 3773 Not Collective 3774 3775 Input Parameter: 3776 . snes - nonlinear solver context 3777 3778 Output Parameter: 3779 . type - SNES method (a character string) 3780 3781 Level: intermediate 3782 3783 .keywords: SNES, nonlinear, get, type, name 3784 @*/ 3785 PetscErrorCode SNESGetType(SNES snes,SNESType *type) 3786 { 3787 PetscFunctionBegin; 3788 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3789 PetscValidPointer(type,2); 3790 *type = ((PetscObject)snes)->type_name; 3791 PetscFunctionReturn(0); 3792 } 3793 3794 #undef __FUNCT__ 3795 #define __FUNCT__ "SNESGetSolution" 3796 /*@ 3797 SNESGetSolution - Returns the vector where the approximate solution is 3798 stored. This is the fine grid solution when using SNESSetGridSequence(). 3799 3800 Not Collective, but Vec is parallel if SNES is parallel 3801 3802 Input Parameter: 3803 . snes - the SNES context 3804 3805 Output Parameter: 3806 . x - the solution 3807 3808 Level: intermediate 3809 3810 .keywords: SNES, nonlinear, get, solution 3811 3812 .seealso: SNESGetSolutionUpdate(), SNESGetFunction() 3813 @*/ 3814 PetscErrorCode SNESGetSolution(SNES snes,Vec *x) 3815 { 3816 PetscFunctionBegin; 3817 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3818 PetscValidPointer(x,2); 3819 *x = snes->vec_sol; 3820 PetscFunctionReturn(0); 3821 } 3822 3823 #undef __FUNCT__ 3824 #define __FUNCT__ "SNESGetSolutionUpdate" 3825 /*@ 3826 SNESGetSolutionUpdate - Returns the vector where the solution update is 3827 stored. 3828 3829 Not Collective, but Vec is parallel if SNES is parallel 3830 3831 Input Parameter: 3832 . snes - the SNES context 3833 3834 Output Parameter: 3835 . x - the solution update 3836 3837 Level: advanced 3838 3839 .keywords: SNES, nonlinear, get, solution, update 3840 3841 .seealso: SNESGetSolution(), SNESGetFunction() 3842 @*/ 3843 PetscErrorCode SNESGetSolutionUpdate(SNES snes,Vec *x) 3844 { 3845 PetscFunctionBegin; 3846 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3847 PetscValidPointer(x,2); 3848 *x = snes->vec_sol_update; 3849 PetscFunctionReturn(0); 3850 } 3851 3852 #undef __FUNCT__ 3853 #define __FUNCT__ "SNESGetFunction" 3854 /*@C 3855 SNESGetFunction - Returns the vector where the function is stored. 3856 3857 Not Collective, but Vec is parallel if SNES is parallel. Collective if Vec is requested, but has not been created yet. 3858 3859 Input Parameter: 3860 . snes - the SNES context 3861 3862 Output Parameter: 3863 + r - the vector that is used to store residuals (or NULL if you don't want it) 3864 . SNESFunction- the function (or NULL if you don't want it) 3865 - ctx - the function context (or NULL if you don't want it) 3866 3867 Level: advanced 3868 3869 .keywords: SNES, nonlinear, get, function 3870 3871 .seealso: SNESSetFunction(), SNESGetSolution(), SNESFunction 3872 @*/ 3873 PetscErrorCode SNESGetFunction(SNES snes,Vec *r,PetscErrorCode (**SNESFunction)(SNES,Vec,Vec,void*),void **ctx) 3874 { 3875 PetscErrorCode ierr; 3876 DM dm; 3877 3878 PetscFunctionBegin; 3879 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3880 if (r) { 3881 if (!snes->vec_func) { 3882 if (snes->vec_rhs) { 3883 ierr = VecDuplicate(snes->vec_rhs,&snes->vec_func);CHKERRQ(ierr); 3884 } else if (snes->vec_sol) { 3885 ierr = VecDuplicate(snes->vec_sol,&snes->vec_func);CHKERRQ(ierr); 3886 } else if (snes->dm) { 3887 ierr = DMCreateGlobalVector(snes->dm,&snes->vec_func);CHKERRQ(ierr); 3888 } 3889 } 3890 *r = snes->vec_func; 3891 } 3892 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 3893 ierr = DMSNESGetFunction(dm,SNESFunction,ctx);CHKERRQ(ierr); 3894 PetscFunctionReturn(0); 3895 } 3896 3897 /*@C 3898 SNESGetGS - Returns the GS function and context. 3899 3900 Input Parameter: 3901 . snes - the SNES context 3902 3903 Output Parameter: 3904 + SNESGSFunction - the function (or NULL) 3905 - ctx - the function context (or NULL) 3906 3907 Level: advanced 3908 3909 .keywords: SNES, nonlinear, get, function 3910 3911 .seealso: SNESSetGS(), SNESGetFunction() 3912 @*/ 3913 3914 #undef __FUNCT__ 3915 #define __FUNCT__ "SNESGetGS" 3916 PetscErrorCode SNESGetGS (SNES snes, PetscErrorCode (**SNESGSFunction)(SNES, Vec, Vec, void*), void ** ctx) 3917 { 3918 PetscErrorCode ierr; 3919 DM dm; 3920 3921 PetscFunctionBegin; 3922 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3923 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 3924 ierr = DMSNESGetGS(dm,SNESGSFunction,ctx);CHKERRQ(ierr); 3925 PetscFunctionReturn(0); 3926 } 3927 3928 #undef __FUNCT__ 3929 #define __FUNCT__ "SNESSetOptionsPrefix" 3930 /*@C 3931 SNESSetOptionsPrefix - Sets the prefix used for searching for all 3932 SNES options in the database. 3933 3934 Logically Collective on SNES 3935 3936 Input Parameter: 3937 + snes - the SNES context 3938 - prefix - the prefix to prepend to all option names 3939 3940 Notes: 3941 A hyphen (-) must NOT be given at the beginning of the prefix name. 3942 The first character of all runtime options is AUTOMATICALLY the hyphen. 3943 3944 Level: advanced 3945 3946 .keywords: SNES, set, options, prefix, database 3947 3948 .seealso: SNESSetFromOptions() 3949 @*/ 3950 PetscErrorCode SNESSetOptionsPrefix(SNES snes,const char prefix[]) 3951 { 3952 PetscErrorCode ierr; 3953 3954 PetscFunctionBegin; 3955 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3956 ierr = PetscObjectSetOptionsPrefix((PetscObject)snes,prefix);CHKERRQ(ierr); 3957 if (!snes->ksp) {ierr = SNESGetKSP(snes,&snes->ksp);CHKERRQ(ierr);} 3958 if (snes->linesearch) { 3959 ierr = SNESGetLineSearch(snes,&snes->linesearch);CHKERRQ(ierr); 3960 ierr = PetscObjectSetOptionsPrefix((PetscObject)snes->linesearch,prefix);CHKERRQ(ierr); 3961 } 3962 ierr = KSPSetOptionsPrefix(snes->ksp,prefix);CHKERRQ(ierr); 3963 PetscFunctionReturn(0); 3964 } 3965 3966 #undef __FUNCT__ 3967 #define __FUNCT__ "SNESAppendOptionsPrefix" 3968 /*@C 3969 SNESAppendOptionsPrefix - Appends to the prefix used for searching for all 3970 SNES options in the database. 3971 3972 Logically Collective on SNES 3973 3974 Input Parameters: 3975 + snes - the SNES context 3976 - prefix - the prefix to prepend to all option names 3977 3978 Notes: 3979 A hyphen (-) must NOT be given at the beginning of the prefix name. 3980 The first character of all runtime options is AUTOMATICALLY the hyphen. 3981 3982 Level: advanced 3983 3984 .keywords: SNES, append, options, prefix, database 3985 3986 .seealso: SNESGetOptionsPrefix() 3987 @*/ 3988 PetscErrorCode SNESAppendOptionsPrefix(SNES snes,const char prefix[]) 3989 { 3990 PetscErrorCode ierr; 3991 3992 PetscFunctionBegin; 3993 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3994 ierr = PetscObjectAppendOptionsPrefix((PetscObject)snes,prefix);CHKERRQ(ierr); 3995 if (!snes->ksp) {ierr = SNESGetKSP(snes,&snes->ksp);CHKERRQ(ierr);} 3996 if (snes->linesearch) { 3997 ierr = SNESGetLineSearch(snes,&snes->linesearch);CHKERRQ(ierr); 3998 ierr = PetscObjectAppendOptionsPrefix((PetscObject)snes->linesearch,prefix);CHKERRQ(ierr); 3999 } 4000 ierr = KSPAppendOptionsPrefix(snes->ksp,prefix);CHKERRQ(ierr); 4001 PetscFunctionReturn(0); 4002 } 4003 4004 #undef __FUNCT__ 4005 #define __FUNCT__ "SNESGetOptionsPrefix" 4006 /*@C 4007 SNESGetOptionsPrefix - Sets the prefix used for searching for all 4008 SNES options in the database. 4009 4010 Not Collective 4011 4012 Input Parameter: 4013 . snes - the SNES context 4014 4015 Output Parameter: 4016 . prefix - pointer to the prefix string used 4017 4018 Notes: On the fortran side, the user should pass in a string 'prefix' of 4019 sufficient length to hold the prefix. 4020 4021 Level: advanced 4022 4023 .keywords: SNES, get, options, prefix, database 4024 4025 .seealso: SNESAppendOptionsPrefix() 4026 @*/ 4027 PetscErrorCode SNESGetOptionsPrefix(SNES snes,const char *prefix[]) 4028 { 4029 PetscErrorCode ierr; 4030 4031 PetscFunctionBegin; 4032 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4033 ierr = PetscObjectGetOptionsPrefix((PetscObject)snes,prefix);CHKERRQ(ierr); 4034 PetscFunctionReturn(0); 4035 } 4036 4037 4038 #undef __FUNCT__ 4039 #define __FUNCT__ "SNESRegister" 4040 /*@C 4041 SNESRegister - Adds a method to the nonlinear solver package. 4042 4043 Not collective 4044 4045 Input Parameters: 4046 + name_solver - name of a new user-defined solver 4047 - routine_create - routine to create method context 4048 4049 Notes: 4050 SNESRegister() may be called multiple times to add several user-defined solvers. 4051 4052 Sample usage: 4053 .vb 4054 SNESRegister("my_solver",MySolverCreate); 4055 .ve 4056 4057 Then, your solver can be chosen with the procedural interface via 4058 $ SNESSetType(snes,"my_solver") 4059 or at runtime via the option 4060 $ -snes_type my_solver 4061 4062 Level: advanced 4063 4064 Note: If your function is not being put into a shared library then use SNESRegister() instead 4065 4066 .keywords: SNES, nonlinear, register 4067 4068 .seealso: SNESRegisterAll(), SNESRegisterDestroy() 4069 4070 Level: advanced 4071 @*/ 4072 PetscErrorCode SNESRegister(const char sname[],PetscErrorCode (*function)(SNES)) 4073 { 4074 PetscErrorCode ierr; 4075 4076 PetscFunctionBegin; 4077 ierr = PetscFunctionListAdd(&SNESList,sname,function);CHKERRQ(ierr); 4078 PetscFunctionReturn(0); 4079 } 4080 4081 #undef __FUNCT__ 4082 #define __FUNCT__ "SNESTestLocalMin" 4083 PetscErrorCode SNESTestLocalMin(SNES snes) 4084 { 4085 PetscErrorCode ierr; 4086 PetscInt N,i,j; 4087 Vec u,uh,fh; 4088 PetscScalar value; 4089 PetscReal norm; 4090 4091 PetscFunctionBegin; 4092 ierr = SNESGetSolution(snes,&u);CHKERRQ(ierr); 4093 ierr = VecDuplicate(u,&uh);CHKERRQ(ierr); 4094 ierr = VecDuplicate(u,&fh);CHKERRQ(ierr); 4095 4096 /* currently only works for sequential */ 4097 ierr = PetscPrintf(PETSC_COMM_WORLD,"Testing FormFunction() for local min\n");CHKERRQ(ierr); 4098 ierr = VecGetSize(u,&N);CHKERRQ(ierr); 4099 for (i=0; i<N; i++) { 4100 ierr = VecCopy(u,uh);CHKERRQ(ierr); 4101 ierr = PetscPrintf(PETSC_COMM_WORLD,"i = %D\n",i);CHKERRQ(ierr); 4102 for (j=-10; j<11; j++) { 4103 value = PetscSign(j)*exp(PetscAbs(j)-10.0); 4104 ierr = VecSetValue(uh,i,value,ADD_VALUES);CHKERRQ(ierr); 4105 ierr = SNESComputeFunction(snes,uh,fh);CHKERRQ(ierr); 4106 ierr = VecNorm(fh,NORM_2,&norm);CHKERRQ(ierr); 4107 ierr = PetscPrintf(PETSC_COMM_WORLD," j norm %D %18.16e\n",j,norm);CHKERRQ(ierr); 4108 value = -value; 4109 ierr = VecSetValue(uh,i,value,ADD_VALUES);CHKERRQ(ierr); 4110 } 4111 } 4112 ierr = VecDestroy(&uh);CHKERRQ(ierr); 4113 ierr = VecDestroy(&fh);CHKERRQ(ierr); 4114 PetscFunctionReturn(0); 4115 } 4116 4117 #undef __FUNCT__ 4118 #define __FUNCT__ "SNESKSPSetUseEW" 4119 /*@ 4120 SNESKSPSetUseEW - Sets SNES use Eisenstat-Walker method for 4121 computing relative tolerance for linear solvers within an inexact 4122 Newton method. 4123 4124 Logically Collective on SNES 4125 4126 Input Parameters: 4127 + snes - SNES context 4128 - flag - PETSC_TRUE or PETSC_FALSE 4129 4130 Options Database: 4131 + -snes_ksp_ew - use Eisenstat-Walker method for determining linear system convergence 4132 . -snes_ksp_ew_version ver - version of Eisenstat-Walker method 4133 . -snes_ksp_ew_rtol0 <rtol0> - Sets rtol0 4134 . -snes_ksp_ew_rtolmax <rtolmax> - Sets rtolmax 4135 . -snes_ksp_ew_gamma <gamma> - Sets gamma 4136 . -snes_ksp_ew_alpha <alpha> - Sets alpha 4137 . -snes_ksp_ew_alpha2 <alpha2> - Sets alpha2 4138 - -snes_ksp_ew_threshold <threshold> - Sets threshold 4139 4140 Notes: 4141 Currently, the default is to use a constant relative tolerance for 4142 the inner linear solvers. Alternatively, one can use the 4143 Eisenstat-Walker method, where the relative convergence tolerance 4144 is reset at each Newton iteration according progress of the nonlinear 4145 solver. 4146 4147 Level: advanced 4148 4149 Reference: 4150 S. C. Eisenstat and H. F. Walker, "Choosing the forcing terms in an 4151 inexact Newton method", SISC 17 (1), pp.16-32, 1996. 4152 4153 .keywords: SNES, KSP, Eisenstat, Walker, convergence, test, inexact, Newton 4154 4155 .seealso: SNESKSPGetUseEW(), SNESKSPGetParametersEW(), SNESKSPSetParametersEW() 4156 @*/ 4157 PetscErrorCode SNESKSPSetUseEW(SNES snes,PetscBool flag) 4158 { 4159 PetscFunctionBegin; 4160 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4161 PetscValidLogicalCollectiveBool(snes,flag,2); 4162 snes->ksp_ewconv = flag; 4163 PetscFunctionReturn(0); 4164 } 4165 4166 #undef __FUNCT__ 4167 #define __FUNCT__ "SNESKSPGetUseEW" 4168 /*@ 4169 SNESKSPGetUseEW - Gets if SNES is using Eisenstat-Walker method 4170 for computing relative tolerance for linear solvers within an 4171 inexact Newton method. 4172 4173 Not Collective 4174 4175 Input Parameter: 4176 . snes - SNES context 4177 4178 Output Parameter: 4179 . flag - PETSC_TRUE or PETSC_FALSE 4180 4181 Notes: 4182 Currently, the default is to use a constant relative tolerance for 4183 the inner linear solvers. Alternatively, one can use the 4184 Eisenstat-Walker method, where the relative convergence tolerance 4185 is reset at each Newton iteration according progress of the nonlinear 4186 solver. 4187 4188 Level: advanced 4189 4190 Reference: 4191 S. C. Eisenstat and H. F. Walker, "Choosing the forcing terms in an 4192 inexact Newton method", SISC 17 (1), pp.16-32, 1996. 4193 4194 .keywords: SNES, KSP, Eisenstat, Walker, convergence, test, inexact, Newton 4195 4196 .seealso: SNESKSPSetUseEW(), SNESKSPGetParametersEW(), SNESKSPSetParametersEW() 4197 @*/ 4198 PetscErrorCode SNESKSPGetUseEW(SNES snes, PetscBool *flag) 4199 { 4200 PetscFunctionBegin; 4201 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4202 PetscValidPointer(flag,2); 4203 *flag = snes->ksp_ewconv; 4204 PetscFunctionReturn(0); 4205 } 4206 4207 #undef __FUNCT__ 4208 #define __FUNCT__ "SNESKSPSetParametersEW" 4209 /*@ 4210 SNESKSPSetParametersEW - Sets parameters for Eisenstat-Walker 4211 convergence criteria for the linear solvers within an inexact 4212 Newton method. 4213 4214 Logically Collective on SNES 4215 4216 Input Parameters: 4217 + snes - SNES context 4218 . version - version 1, 2 (default is 2) or 3 4219 . rtol_0 - initial relative tolerance (0 <= rtol_0 < 1) 4220 . rtol_max - maximum relative tolerance (0 <= rtol_max < 1) 4221 . gamma - multiplicative factor for version 2 rtol computation 4222 (0 <= gamma2 <= 1) 4223 . alpha - power for version 2 rtol computation (1 < alpha <= 2) 4224 . alpha2 - power for safeguard 4225 - threshold - threshold for imposing safeguard (0 < threshold < 1) 4226 4227 Note: 4228 Version 3 was contributed by Luis Chacon, June 2006. 4229 4230 Use PETSC_DEFAULT to retain the default for any of the parameters. 4231 4232 Level: advanced 4233 4234 Reference: 4235 S. C. Eisenstat and H. F. Walker, "Choosing the forcing terms in an 4236 inexact Newton method", Utah State University Math. Stat. Dept. Res. 4237 Report 6/94/75, June, 1994, to appear in SIAM J. Sci. Comput. 4238 4239 .keywords: SNES, KSP, Eisenstat, Walker, set, parameters 4240 4241 .seealso: SNESKSPSetUseEW(), SNESKSPGetUseEW(), SNESKSPGetParametersEW() 4242 @*/ 4243 PetscErrorCode SNESKSPSetParametersEW(SNES snes,PetscInt version,PetscReal rtol_0,PetscReal rtol_max,PetscReal gamma,PetscReal alpha,PetscReal alpha2,PetscReal threshold) 4244 { 4245 SNESKSPEW *kctx; 4246 4247 PetscFunctionBegin; 4248 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4249 kctx = (SNESKSPEW*)snes->kspconvctx; 4250 if (!kctx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"No Eisenstat-Walker context existing"); 4251 PetscValidLogicalCollectiveInt(snes,version,2); 4252 PetscValidLogicalCollectiveReal(snes,rtol_0,3); 4253 PetscValidLogicalCollectiveReal(snes,rtol_max,4); 4254 PetscValidLogicalCollectiveReal(snes,gamma,5); 4255 PetscValidLogicalCollectiveReal(snes,alpha,6); 4256 PetscValidLogicalCollectiveReal(snes,alpha2,7); 4257 PetscValidLogicalCollectiveReal(snes,threshold,8); 4258 4259 if (version != PETSC_DEFAULT) kctx->version = version; 4260 if (rtol_0 != PETSC_DEFAULT) kctx->rtol_0 = rtol_0; 4261 if (rtol_max != PETSC_DEFAULT) kctx->rtol_max = rtol_max; 4262 if (gamma != PETSC_DEFAULT) kctx->gamma = gamma; 4263 if (alpha != PETSC_DEFAULT) kctx->alpha = alpha; 4264 if (alpha2 != PETSC_DEFAULT) kctx->alpha2 = alpha2; 4265 if (threshold != PETSC_DEFAULT) kctx->threshold = threshold; 4266 4267 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); 4268 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); 4269 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); 4270 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); 4271 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); 4272 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); 4273 PetscFunctionReturn(0); 4274 } 4275 4276 #undef __FUNCT__ 4277 #define __FUNCT__ "SNESKSPGetParametersEW" 4278 /*@ 4279 SNESKSPGetParametersEW - Gets parameters for Eisenstat-Walker 4280 convergence criteria for the linear solvers within an inexact 4281 Newton method. 4282 4283 Not Collective 4284 4285 Input Parameters: 4286 snes - SNES context 4287 4288 Output Parameters: 4289 + version - version 1, 2 (default is 2) or 3 4290 . rtol_0 - initial relative tolerance (0 <= rtol_0 < 1) 4291 . rtol_max - maximum relative tolerance (0 <= rtol_max < 1) 4292 . gamma - multiplicative factor for version 2 rtol computation (0 <= gamma2 <= 1) 4293 . alpha - power for version 2 rtol computation (1 < alpha <= 2) 4294 . alpha2 - power for safeguard 4295 - threshold - threshold for imposing safeguard (0 < threshold < 1) 4296 4297 Level: advanced 4298 4299 .keywords: SNES, KSP, Eisenstat, Walker, get, parameters 4300 4301 .seealso: SNESKSPSetUseEW(), SNESKSPGetUseEW(), SNESKSPSetParametersEW() 4302 @*/ 4303 PetscErrorCode SNESKSPGetParametersEW(SNES snes,PetscInt *version,PetscReal *rtol_0,PetscReal *rtol_max,PetscReal *gamma,PetscReal *alpha,PetscReal *alpha2,PetscReal *threshold) 4304 { 4305 SNESKSPEW *kctx; 4306 4307 PetscFunctionBegin; 4308 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4309 kctx = (SNESKSPEW*)snes->kspconvctx; 4310 if (!kctx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"No Eisenstat-Walker context existing"); 4311 if (version) *version = kctx->version; 4312 if (rtol_0) *rtol_0 = kctx->rtol_0; 4313 if (rtol_max) *rtol_max = kctx->rtol_max; 4314 if (gamma) *gamma = kctx->gamma; 4315 if (alpha) *alpha = kctx->alpha; 4316 if (alpha2) *alpha2 = kctx->alpha2; 4317 if (threshold) *threshold = kctx->threshold; 4318 PetscFunctionReturn(0); 4319 } 4320 4321 #undef __FUNCT__ 4322 #define __FUNCT__ "SNESKSPEW_PreSolve" 4323 PetscErrorCode SNESKSPEW_PreSolve(KSP ksp, Vec b, Vec x, SNES snes) 4324 { 4325 PetscErrorCode ierr; 4326 SNESKSPEW *kctx = (SNESKSPEW*)snes->kspconvctx; 4327 PetscReal rtol = PETSC_DEFAULT,stol; 4328 4329 PetscFunctionBegin; 4330 if (!snes->ksp_ewconv) PetscFunctionReturn(0); 4331 if (!snes->iter) rtol = kctx->rtol_0; /* first time in, so use the original user rtol */ 4332 else { 4333 if (kctx->version == 1) { 4334 rtol = (snes->norm - kctx->lresid_last)/kctx->norm_last; 4335 if (rtol < 0.0) rtol = -rtol; 4336 stol = PetscPowReal(kctx->rtol_last,kctx->alpha2); 4337 if (stol > kctx->threshold) rtol = PetscMax(rtol,stol); 4338 } else if (kctx->version == 2) { 4339 rtol = kctx->gamma * PetscPowReal(snes->norm/kctx->norm_last,kctx->alpha); 4340 stol = kctx->gamma * PetscPowReal(kctx->rtol_last,kctx->alpha); 4341 if (stol > kctx->threshold) rtol = PetscMax(rtol,stol); 4342 } else if (kctx->version == 3) { /* contributed by Luis Chacon, June 2006. */ 4343 rtol = kctx->gamma * PetscPowReal(snes->norm/kctx->norm_last,kctx->alpha); 4344 /* safeguard: avoid sharp decrease of rtol */ 4345 stol = kctx->gamma*PetscPowReal(kctx->rtol_last,kctx->alpha); 4346 stol = PetscMax(rtol,stol); 4347 rtol = PetscMin(kctx->rtol_0,stol); 4348 /* safeguard: avoid oversolving */ 4349 stol = kctx->gamma*(snes->ttol)/snes->norm; 4350 stol = PetscMax(rtol,stol); 4351 rtol = PetscMin(kctx->rtol_0,stol); 4352 } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Only versions 1, 2 or 3 are supported: %D",kctx->version); 4353 } 4354 /* safeguard: avoid rtol greater than one */ 4355 rtol = PetscMin(rtol,kctx->rtol_max); 4356 ierr = KSPSetTolerances(ksp,rtol,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);CHKERRQ(ierr); 4357 ierr = PetscInfo3(snes,"iter %D, Eisenstat-Walker (version %D) KSP rtol=%G\n",snes->iter,kctx->version,rtol);CHKERRQ(ierr); 4358 PetscFunctionReturn(0); 4359 } 4360 4361 #undef __FUNCT__ 4362 #define __FUNCT__ "SNESKSPEW_PostSolve" 4363 PetscErrorCode SNESKSPEW_PostSolve(KSP ksp, Vec b, Vec x, SNES snes) 4364 { 4365 PetscErrorCode ierr; 4366 SNESKSPEW *kctx = (SNESKSPEW*)snes->kspconvctx; 4367 PCSide pcside; 4368 Vec lres; 4369 4370 PetscFunctionBegin; 4371 if (!snes->ksp_ewconv) PetscFunctionReturn(0); 4372 ierr = KSPGetTolerances(ksp,&kctx->rtol_last,0,0,0);CHKERRQ(ierr); 4373 ierr = SNESGetFunctionNorm(snes,&kctx->norm_last);CHKERRQ(ierr); 4374 if (kctx->version == 1) { 4375 ierr = KSPGetPCSide(ksp,&pcside);CHKERRQ(ierr); 4376 if (pcside == PC_RIGHT) { /* XXX Should we also test KSP_UNPRECONDITIONED_NORM ? */ 4377 /* KSP residual is true linear residual */ 4378 ierr = KSPGetResidualNorm(ksp,&kctx->lresid_last);CHKERRQ(ierr); 4379 } else { 4380 /* KSP residual is preconditioned residual */ 4381 /* compute true linear residual norm */ 4382 ierr = VecDuplicate(b,&lres);CHKERRQ(ierr); 4383 ierr = MatMult(snes->jacobian,x,lres);CHKERRQ(ierr); 4384 ierr = VecAYPX(lres,-1.0,b);CHKERRQ(ierr); 4385 ierr = VecNorm(lres,NORM_2,&kctx->lresid_last);CHKERRQ(ierr); 4386 ierr = VecDestroy(&lres);CHKERRQ(ierr); 4387 } 4388 } 4389 PetscFunctionReturn(0); 4390 } 4391 4392 #undef __FUNCT__ 4393 #define __FUNCT__ "SNESGetKSP" 4394 /*@ 4395 SNESGetKSP - Returns the KSP context for a SNES solver. 4396 4397 Not Collective, but if SNES object is parallel, then KSP object is parallel 4398 4399 Input Parameter: 4400 . snes - the SNES context 4401 4402 Output Parameter: 4403 . ksp - the KSP context 4404 4405 Notes: 4406 The user can then directly manipulate the KSP context to set various 4407 options, etc. Likewise, the user can then extract and manipulate the 4408 PC contexts as well. 4409 4410 Level: beginner 4411 4412 .keywords: SNES, nonlinear, get, KSP, context 4413 4414 .seealso: KSPGetPC(), SNESCreate(), KSPCreate(), SNESSetKSP() 4415 @*/ 4416 PetscErrorCode SNESGetKSP(SNES snes,KSP *ksp) 4417 { 4418 PetscErrorCode ierr; 4419 4420 PetscFunctionBegin; 4421 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4422 PetscValidPointer(ksp,2); 4423 4424 if (!snes->ksp) { 4425 ierr = KSPCreate(PetscObjectComm((PetscObject)snes),&snes->ksp);CHKERRQ(ierr); 4426 ierr = PetscObjectIncrementTabLevel((PetscObject)snes->ksp,(PetscObject)snes,1);CHKERRQ(ierr); 4427 ierr = PetscLogObjectParent(snes,snes->ksp);CHKERRQ(ierr); 4428 4429 ierr = KSPSetPreSolve(snes->ksp,(PetscErrorCode (*)(KSP,Vec,Vec,void*))SNESKSPEW_PreSolve,snes);CHKERRQ(ierr); 4430 ierr = KSPSetPostSolve(snes->ksp,(PetscErrorCode (*)(KSP,Vec,Vec,void*))SNESKSPEW_PostSolve,snes);CHKERRQ(ierr); 4431 } 4432 *ksp = snes->ksp; 4433 PetscFunctionReturn(0); 4434 } 4435 4436 4437 #include <petsc-private/dmimpl.h> 4438 #undef __FUNCT__ 4439 #define __FUNCT__ "SNESSetDM" 4440 /*@ 4441 SNESSetDM - Sets the DM that may be used by some preconditioners 4442 4443 Logically Collective on SNES 4444 4445 Input Parameters: 4446 + snes - the preconditioner context 4447 - dm - the dm 4448 4449 Level: intermediate 4450 4451 .seealso: SNESGetDM(), KSPSetDM(), KSPGetDM() 4452 @*/ 4453 PetscErrorCode SNESSetDM(SNES snes,DM dm) 4454 { 4455 PetscErrorCode ierr; 4456 KSP ksp; 4457 DMSNES sdm; 4458 4459 PetscFunctionBegin; 4460 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4461 if (dm) {ierr = PetscObjectReference((PetscObject)dm);CHKERRQ(ierr);} 4462 if (snes->dm) { /* Move the DMSNES context over to the new DM unless the new DM already has one */ 4463 if (snes->dm->dmsnes && snes->dmAuto && !dm->dmsnes) { 4464 ierr = DMCopyDMSNES(snes->dm,dm);CHKERRQ(ierr); 4465 ierr = DMGetDMSNES(snes->dm,&sdm);CHKERRQ(ierr); 4466 if (sdm->originaldm == snes->dm) sdm->originaldm = dm; /* Grant write privileges to the replacement DM */ 4467 } 4468 ierr = DMDestroy(&snes->dm);CHKERRQ(ierr); 4469 } 4470 snes->dm = dm; 4471 snes->dmAuto = PETSC_FALSE; 4472 4473 ierr = SNESGetKSP(snes,&ksp);CHKERRQ(ierr); 4474 ierr = KSPSetDM(ksp,dm);CHKERRQ(ierr); 4475 ierr = KSPSetDMActive(ksp,PETSC_FALSE);CHKERRQ(ierr); 4476 if (snes->pc) { 4477 ierr = SNESSetDM(snes->pc, snes->dm);CHKERRQ(ierr); 4478 ierr = SNESSetPCSide(snes,snes->pcside);CHKERRQ(ierr); 4479 } 4480 PetscFunctionReturn(0); 4481 } 4482 4483 #undef __FUNCT__ 4484 #define __FUNCT__ "SNESGetDM" 4485 /*@ 4486 SNESGetDM - Gets the DM that may be used by some preconditioners 4487 4488 Not Collective but DM obtained is parallel on SNES 4489 4490 Input Parameter: 4491 . snes - the preconditioner context 4492 4493 Output Parameter: 4494 . dm - the dm 4495 4496 Level: intermediate 4497 4498 .seealso: SNESSetDM(), KSPSetDM(), KSPGetDM() 4499 @*/ 4500 PetscErrorCode SNESGetDM(SNES snes,DM *dm) 4501 { 4502 PetscErrorCode ierr; 4503 4504 PetscFunctionBegin; 4505 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4506 if (!snes->dm) { 4507 ierr = DMShellCreate(PetscObjectComm((PetscObject)snes),&snes->dm);CHKERRQ(ierr); 4508 snes->dmAuto = PETSC_TRUE; 4509 } 4510 *dm = snes->dm; 4511 PetscFunctionReturn(0); 4512 } 4513 4514 #undef __FUNCT__ 4515 #define __FUNCT__ "SNESSetPC" 4516 /*@ 4517 SNESSetPC - Sets the nonlinear preconditioner to be used. 4518 4519 Collective on SNES 4520 4521 Input Parameters: 4522 + snes - iterative context obtained from SNESCreate() 4523 - pc - the preconditioner object 4524 4525 Notes: 4526 Use SNESGetPC() to retrieve the preconditioner context (for example, 4527 to configure it using the API). 4528 4529 Level: developer 4530 4531 .keywords: SNES, set, precondition 4532 .seealso: SNESGetPC() 4533 @*/ 4534 PetscErrorCode SNESSetPC(SNES snes, SNES pc) 4535 { 4536 PetscErrorCode ierr; 4537 4538 PetscFunctionBegin; 4539 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); 4540 PetscValidHeaderSpecific(pc, SNES_CLASSID, 2); 4541 PetscCheckSameComm(snes, 1, pc, 2); 4542 ierr = PetscObjectReference((PetscObject) pc);CHKERRQ(ierr); 4543 ierr = SNESDestroy(&snes->pc);CHKERRQ(ierr); 4544 snes->pc = pc; 4545 ierr = PetscLogObjectParent(snes, snes->pc);CHKERRQ(ierr); 4546 PetscFunctionReturn(0); 4547 } 4548 4549 #undef __FUNCT__ 4550 #define __FUNCT__ "SNESGetPC" 4551 /*@ 4552 SNESGetPC - Returns a pointer to the nonlinear preconditioning context set with SNESSetPC(). 4553 4554 Not Collective 4555 4556 Input Parameter: 4557 . snes - iterative context obtained from SNESCreate() 4558 4559 Output Parameter: 4560 . pc - preconditioner context 4561 4562 Level: developer 4563 4564 .keywords: SNES, get, preconditioner 4565 .seealso: SNESSetPC() 4566 @*/ 4567 PetscErrorCode SNESGetPC(SNES snes, SNES *pc) 4568 { 4569 PetscErrorCode ierr; 4570 const char *optionsprefix; 4571 4572 PetscFunctionBegin; 4573 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); 4574 PetscValidPointer(pc, 2); 4575 if (!snes->pc) { 4576 ierr = SNESCreate(PetscObjectComm((PetscObject)snes),&snes->pc);CHKERRQ(ierr); 4577 ierr = PetscObjectIncrementTabLevel((PetscObject)snes->pc,(PetscObject)snes,1);CHKERRQ(ierr); 4578 ierr = PetscLogObjectParent(snes,snes->pc);CHKERRQ(ierr); 4579 ierr = SNESGetOptionsPrefix(snes,&optionsprefix);CHKERRQ(ierr); 4580 ierr = SNESSetOptionsPrefix(snes->pc,optionsprefix);CHKERRQ(ierr); 4581 ierr = SNESAppendOptionsPrefix(snes->pc,"npc_");CHKERRQ(ierr); 4582 } 4583 *pc = snes->pc; 4584 PetscFunctionReturn(0); 4585 } 4586 4587 #undef __FUNCT__ 4588 #define __FUNCT__ "SNESSetPCSide" 4589 /*@ 4590 SNESSetPCSide - Sets the preconditioning side. 4591 4592 Logically Collective on SNES 4593 4594 Input Parameter: 4595 . snes - iterative context obtained from SNESCreate() 4596 4597 Output Parameter: 4598 . side - the preconditioning side, where side is one of 4599 .vb 4600 PC_LEFT - left preconditioning (default) 4601 PC_RIGHT - right preconditioning 4602 .ve 4603 4604 Options Database Keys: 4605 . -snes_pc_side <right,left> 4606 4607 Level: intermediate 4608 4609 .keywords: SNES, set, right, left, side, preconditioner, flag 4610 4611 .seealso: SNESGetPCSide(), KSPSetPCSide() 4612 @*/ 4613 PetscErrorCode SNESSetPCSide(SNES snes,PCSide side) 4614 { 4615 PetscFunctionBegin; 4616 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4617 PetscValidLogicalCollectiveEnum(snes,side,2); 4618 snes->pcside = side; 4619 PetscFunctionReturn(0); 4620 } 4621 4622 #undef __FUNCT__ 4623 #define __FUNCT__ "SNESGetPCSide" 4624 /*@ 4625 SNESGetPCSide - Gets the preconditioning side. 4626 4627 Not Collective 4628 4629 Input Parameter: 4630 . snes - iterative context obtained from SNESCreate() 4631 4632 Output Parameter: 4633 . side - the preconditioning side, where side is one of 4634 .vb 4635 PC_LEFT - left preconditioning (default) 4636 PC_RIGHT - right preconditioning 4637 .ve 4638 4639 Level: intermediate 4640 4641 .keywords: SNES, get, right, left, side, preconditioner, flag 4642 4643 .seealso: SNESSetPCSide(), KSPGetPCSide() 4644 @*/ 4645 PetscErrorCode SNESGetPCSide(SNES snes,PCSide *side) 4646 { 4647 PetscFunctionBegin; 4648 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4649 PetscValidPointer(side,2); 4650 *side = snes->pcside; 4651 PetscFunctionReturn(0); 4652 } 4653 4654 #undef __FUNCT__ 4655 #define __FUNCT__ "SNESSetLineSearch" 4656 /*@ 4657 SNESSetLineSearch - Sets the linesearch on the SNES instance. 4658 4659 Collective on SNES 4660 4661 Input Parameters: 4662 + snes - iterative context obtained from SNESCreate() 4663 - linesearch - the linesearch object 4664 4665 Notes: 4666 Use SNESGetLineSearch() to retrieve the preconditioner context (for example, 4667 to configure it using the API). 4668 4669 Level: developer 4670 4671 .keywords: SNES, set, linesearch 4672 .seealso: SNESGetLineSearch() 4673 @*/ 4674 PetscErrorCode SNESSetLineSearch(SNES snes, SNESLineSearch linesearch) 4675 { 4676 PetscErrorCode ierr; 4677 4678 PetscFunctionBegin; 4679 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); 4680 PetscValidHeaderSpecific(linesearch, SNESLINESEARCH_CLASSID, 2); 4681 PetscCheckSameComm(snes, 1, linesearch, 2); 4682 ierr = PetscObjectReference((PetscObject) linesearch);CHKERRQ(ierr); 4683 ierr = SNESLineSearchDestroy(&snes->linesearch);CHKERRQ(ierr); 4684 4685 snes->linesearch = linesearch; 4686 4687 ierr = PetscLogObjectParent(snes, snes->linesearch);CHKERRQ(ierr); 4688 PetscFunctionReturn(0); 4689 } 4690 4691 #undef __FUNCT__ 4692 #define __FUNCT__ "SNESGetLineSearch" 4693 /*@ 4694 SNESGetLineSearch - Returns a pointer to the line search context set with SNESSetLineSearch() 4695 or creates a default line search instance associated with the SNES and returns it. 4696 4697 Not Collective 4698 4699 Input Parameter: 4700 . snes - iterative context obtained from SNESCreate() 4701 4702 Output Parameter: 4703 . linesearch - linesearch context 4704 4705 Level: developer 4706 4707 .keywords: SNES, get, linesearch 4708 .seealso: SNESSetLineSearch() 4709 @*/ 4710 PetscErrorCode SNESGetLineSearch(SNES snes, SNESLineSearch *linesearch) 4711 { 4712 PetscErrorCode ierr; 4713 const char *optionsprefix; 4714 4715 PetscFunctionBegin; 4716 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); 4717 PetscValidPointer(linesearch, 2); 4718 if (!snes->linesearch) { 4719 ierr = SNESGetOptionsPrefix(snes, &optionsprefix);CHKERRQ(ierr); 4720 ierr = SNESLineSearchCreate(PetscObjectComm((PetscObject)snes), &snes->linesearch);CHKERRQ(ierr); 4721 ierr = SNESLineSearchSetSNES(snes->linesearch, snes);CHKERRQ(ierr); 4722 ierr = SNESLineSearchAppendOptionsPrefix(snes->linesearch, optionsprefix);CHKERRQ(ierr); 4723 ierr = PetscObjectIncrementTabLevel((PetscObject) snes->linesearch, (PetscObject) snes, 1);CHKERRQ(ierr); 4724 ierr = PetscLogObjectParent(snes, snes->linesearch);CHKERRQ(ierr); 4725 } 4726 *linesearch = snes->linesearch; 4727 PetscFunctionReturn(0); 4728 } 4729 4730 #if defined(PETSC_HAVE_MATLAB_ENGINE) 4731 #include <mex.h> 4732 4733 typedef struct {char *funcname; mxArray *ctx;} SNESMatlabContext; 4734 4735 #undef __FUNCT__ 4736 #define __FUNCT__ "SNESComputeFunction_Matlab" 4737 /* 4738 SNESComputeFunction_Matlab - Calls the function that has been set with SNESSetFunctionMatlab(). 4739 4740 Collective on SNES 4741 4742 Input Parameters: 4743 + snes - the SNES context 4744 - x - input vector 4745 4746 Output Parameter: 4747 . y - function vector, as set by SNESSetFunction() 4748 4749 Notes: 4750 SNESComputeFunction() is typically used within nonlinear solvers 4751 implementations, so most users would not generally call this routine 4752 themselves. 4753 4754 Level: developer 4755 4756 .keywords: SNES, nonlinear, compute, function 4757 4758 .seealso: SNESSetFunction(), SNESGetFunction() 4759 */ 4760 PetscErrorCode SNESComputeFunction_Matlab(SNES snes,Vec x,Vec y, void *ctx) 4761 { 4762 PetscErrorCode ierr; 4763 SNESMatlabContext *sctx = (SNESMatlabContext*)ctx; 4764 int nlhs = 1,nrhs = 5; 4765 mxArray *plhs[1],*prhs[5]; 4766 long long int lx = 0,ly = 0,ls = 0; 4767 4768 PetscFunctionBegin; 4769 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4770 PetscValidHeaderSpecific(x,VEC_CLASSID,2); 4771 PetscValidHeaderSpecific(y,VEC_CLASSID,3); 4772 PetscCheckSameComm(snes,1,x,2); 4773 PetscCheckSameComm(snes,1,y,3); 4774 4775 /* call Matlab function in ctx with arguments x and y */ 4776 4777 ierr = PetscMemcpy(&ls,&snes,sizeof(snes));CHKERRQ(ierr); 4778 ierr = PetscMemcpy(&lx,&x,sizeof(x));CHKERRQ(ierr); 4779 ierr = PetscMemcpy(&ly,&y,sizeof(x));CHKERRQ(ierr); 4780 prhs[0] = mxCreateDoubleScalar((double)ls); 4781 prhs[1] = mxCreateDoubleScalar((double)lx); 4782 prhs[2] = mxCreateDoubleScalar((double)ly); 4783 prhs[3] = mxCreateString(sctx->funcname); 4784 prhs[4] = sctx->ctx; 4785 ierr = mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscSNESComputeFunctionInternal");CHKERRQ(ierr); 4786 ierr = mxGetScalar(plhs[0]);CHKERRQ(ierr); 4787 mxDestroyArray(prhs[0]); 4788 mxDestroyArray(prhs[1]); 4789 mxDestroyArray(prhs[2]); 4790 mxDestroyArray(prhs[3]); 4791 mxDestroyArray(plhs[0]); 4792 PetscFunctionReturn(0); 4793 } 4794 4795 #undef __FUNCT__ 4796 #define __FUNCT__ "SNESSetFunctionMatlab" 4797 /* 4798 SNESSetFunctionMatlab - Sets the function evaluation routine and function 4799 vector for use by the SNES routines in solving systems of nonlinear 4800 equations from MATLAB. Here the function is a string containing the name of a MATLAB function 4801 4802 Logically Collective on SNES 4803 4804 Input Parameters: 4805 + snes - the SNES context 4806 . r - vector to store function value 4807 - SNESFunction - function evaluation routine 4808 4809 Notes: 4810 The Newton-like methods typically solve linear systems of the form 4811 $ f'(x) x = -f(x), 4812 where f'(x) denotes the Jacobian matrix and f(x) is the function. 4813 4814 Level: beginner 4815 4816 Developer Note: This bleeds the allocated memory SNESMatlabContext *sctx; 4817 4818 .keywords: SNES, nonlinear, set, function 4819 4820 .seealso: SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESSetFunction() 4821 */ 4822 PetscErrorCode SNESSetFunctionMatlab(SNES snes,Vec r,const char *SNESFunction,mxArray *ctx) 4823 { 4824 PetscErrorCode ierr; 4825 SNESMatlabContext *sctx; 4826 4827 PetscFunctionBegin; 4828 /* currently sctx is memory bleed */ 4829 ierr = PetscMalloc(sizeof(SNESMatlabContext),&sctx);CHKERRQ(ierr); 4830 ierr = PetscStrallocpy(SNESFunction,&sctx->funcname);CHKERRQ(ierr); 4831 /* 4832 This should work, but it doesn't 4833 sctx->ctx = ctx; 4834 mexMakeArrayPersistent(sctx->ctx); 4835 */ 4836 sctx->ctx = mxDuplicateArray(ctx); 4837 ierr = SNESSetFunction(snes,r,SNESComputeFunction_Matlab,sctx);CHKERRQ(ierr); 4838 PetscFunctionReturn(0); 4839 } 4840 4841 #undef __FUNCT__ 4842 #define __FUNCT__ "SNESComputeJacobian_Matlab" 4843 /* 4844 SNESComputeJacobian_Matlab - Calls the function that has been set with SNESSetJacobianMatlab(). 4845 4846 Collective on SNES 4847 4848 Input Parameters: 4849 + snes - the SNES context 4850 . x - input vector 4851 . A, B - the matrices 4852 - ctx - user context 4853 4854 Output Parameter: 4855 . flag - structure of the matrix 4856 4857 Level: developer 4858 4859 .keywords: SNES, nonlinear, compute, function 4860 4861 .seealso: SNESSetFunction(), SNESGetFunction() 4862 @*/ 4863 PetscErrorCode SNESComputeJacobian_Matlab(SNES snes,Vec x,Mat *A,Mat *B,MatStructure *flag, void *ctx) 4864 { 4865 PetscErrorCode ierr; 4866 SNESMatlabContext *sctx = (SNESMatlabContext*)ctx; 4867 int nlhs = 2,nrhs = 6; 4868 mxArray *plhs[2],*prhs[6]; 4869 long long int lx = 0,lA = 0,ls = 0, lB = 0; 4870 4871 PetscFunctionBegin; 4872 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4873 PetscValidHeaderSpecific(x,VEC_CLASSID,2); 4874 4875 /* call Matlab function in ctx with arguments x and y */ 4876 4877 ierr = PetscMemcpy(&ls,&snes,sizeof(snes));CHKERRQ(ierr); 4878 ierr = PetscMemcpy(&lx,&x,sizeof(x));CHKERRQ(ierr); 4879 ierr = PetscMemcpy(&lA,A,sizeof(x));CHKERRQ(ierr); 4880 ierr = PetscMemcpy(&lB,B,sizeof(x));CHKERRQ(ierr); 4881 prhs[0] = mxCreateDoubleScalar((double)ls); 4882 prhs[1] = mxCreateDoubleScalar((double)lx); 4883 prhs[2] = mxCreateDoubleScalar((double)lA); 4884 prhs[3] = mxCreateDoubleScalar((double)lB); 4885 prhs[4] = mxCreateString(sctx->funcname); 4886 prhs[5] = sctx->ctx; 4887 ierr = mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscSNESComputeJacobianInternal");CHKERRQ(ierr); 4888 ierr = mxGetScalar(plhs[0]);CHKERRQ(ierr); 4889 *flag = (MatStructure) mxGetScalar(plhs[1]);CHKERRQ(ierr); 4890 mxDestroyArray(prhs[0]); 4891 mxDestroyArray(prhs[1]); 4892 mxDestroyArray(prhs[2]); 4893 mxDestroyArray(prhs[3]); 4894 mxDestroyArray(prhs[4]); 4895 mxDestroyArray(plhs[0]); 4896 mxDestroyArray(plhs[1]); 4897 PetscFunctionReturn(0); 4898 } 4899 4900 #undef __FUNCT__ 4901 #define __FUNCT__ "SNESSetJacobianMatlab" 4902 /* 4903 SNESSetJacobianMatlab - Sets the Jacobian function evaluation routine and two empty Jacobian matrices 4904 vector for use by the SNES routines in solving systems of nonlinear 4905 equations from MATLAB. Here the function is a string containing the name of a MATLAB function 4906 4907 Logically Collective on SNES 4908 4909 Input Parameters: 4910 + snes - the SNES context 4911 . A,B - Jacobian matrices 4912 . SNESJacobianFunction - function evaluation routine 4913 - ctx - user context 4914 4915 Level: developer 4916 4917 Developer Note: This bleeds the allocated memory SNESMatlabContext *sctx; 4918 4919 .keywords: SNES, nonlinear, set, function 4920 4921 .seealso: SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESSetFunction(), SNESJacobianFunction 4922 */ 4923 PetscErrorCode SNESSetJacobianMatlab(SNES snes,Mat A,Mat B,const char *SNESJacobianFunction,mxArray *ctx) 4924 { 4925 PetscErrorCode ierr; 4926 SNESMatlabContext *sctx; 4927 4928 PetscFunctionBegin; 4929 /* currently sctx is memory bleed */ 4930 ierr = PetscMalloc(sizeof(SNESMatlabContext),&sctx);CHKERRQ(ierr); 4931 ierr = PetscStrallocpy(SNESJacobianFunction,&sctx->funcname);CHKERRQ(ierr); 4932 /* 4933 This should work, but it doesn't 4934 sctx->ctx = ctx; 4935 mexMakeArrayPersistent(sctx->ctx); 4936 */ 4937 sctx->ctx = mxDuplicateArray(ctx); 4938 ierr = SNESSetJacobian(snes,A,B,SNESComputeJacobian_Matlab,sctx);CHKERRQ(ierr); 4939 PetscFunctionReturn(0); 4940 } 4941 4942 #undef __FUNCT__ 4943 #define __FUNCT__ "SNESMonitor_Matlab" 4944 /* 4945 SNESMonitor_Matlab - Calls the function that has been set with SNESMonitorSetMatlab(). 4946 4947 Collective on SNES 4948 4949 .seealso: SNESSetFunction(), SNESGetFunction() 4950 @*/ 4951 PetscErrorCode SNESMonitor_Matlab(SNES snes,PetscInt it, PetscReal fnorm, void *ctx) 4952 { 4953 PetscErrorCode ierr; 4954 SNESMatlabContext *sctx = (SNESMatlabContext*)ctx; 4955 int nlhs = 1,nrhs = 6; 4956 mxArray *plhs[1],*prhs[6]; 4957 long long int lx = 0,ls = 0; 4958 Vec x = snes->vec_sol; 4959 4960 PetscFunctionBegin; 4961 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4962 4963 ierr = PetscMemcpy(&ls,&snes,sizeof(snes));CHKERRQ(ierr); 4964 ierr = PetscMemcpy(&lx,&x,sizeof(x));CHKERRQ(ierr); 4965 prhs[0] = mxCreateDoubleScalar((double)ls); 4966 prhs[1] = mxCreateDoubleScalar((double)it); 4967 prhs[2] = mxCreateDoubleScalar((double)fnorm); 4968 prhs[3] = mxCreateDoubleScalar((double)lx); 4969 prhs[4] = mxCreateString(sctx->funcname); 4970 prhs[5] = sctx->ctx; 4971 ierr = mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscSNESMonitorInternal");CHKERRQ(ierr); 4972 ierr = mxGetScalar(plhs[0]);CHKERRQ(ierr); 4973 mxDestroyArray(prhs[0]); 4974 mxDestroyArray(prhs[1]); 4975 mxDestroyArray(prhs[2]); 4976 mxDestroyArray(prhs[3]); 4977 mxDestroyArray(prhs[4]); 4978 mxDestroyArray(plhs[0]); 4979 PetscFunctionReturn(0); 4980 } 4981 4982 #undef __FUNCT__ 4983 #define __FUNCT__ "SNESMonitorSetMatlab" 4984 /* 4985 SNESMonitorSetMatlab - Sets the monitor function from MATLAB 4986 4987 Level: developer 4988 4989 Developer Note: This bleeds the allocated memory SNESMatlabContext *sctx; 4990 4991 .keywords: SNES, nonlinear, set, function 4992 4993 .seealso: SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESSetFunction() 4994 */ 4995 PetscErrorCode SNESMonitorSetMatlab(SNES snes,const char *SNESMonitorFunction,mxArray *ctx) 4996 { 4997 PetscErrorCode ierr; 4998 SNESMatlabContext *sctx; 4999 5000 PetscFunctionBegin; 5001 /* currently sctx is memory bleed */ 5002 ierr = PetscMalloc(sizeof(SNESMatlabContext),&sctx);CHKERRQ(ierr); 5003 ierr = PetscStrallocpy(SNESMonitorFunction,&sctx->funcname);CHKERRQ(ierr); 5004 /* 5005 This should work, but it doesn't 5006 sctx->ctx = ctx; 5007 mexMakeArrayPersistent(sctx->ctx); 5008 */ 5009 sctx->ctx = mxDuplicateArray(ctx); 5010 ierr = SNESMonitorSet(snes,SNESMonitor_Matlab,sctx,NULL);CHKERRQ(ierr); 5011 PetscFunctionReturn(0); 5012 } 5013 5014 #endif 5015