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