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