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