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