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 (sdm->ops->computefunction) { 2042 ierr = PetscLogEventBegin(SNES_FunctionEval,snes,x,y,0);CHKERRQ(ierr); 2043 PetscStackPush("SNES user function"); 2044 ierr = (*sdm->ops->computefunction)(snes,x,y,sdm->functionctx);CHKERRQ(ierr); 2045 PetscStackPop; 2046 ierr = PetscLogEventEnd(SNES_FunctionEval,snes,x,y,0);CHKERRQ(ierr); 2047 } else if (snes->vec_rhs) { 2048 ierr = MatMult(snes->jacobian, x, y);CHKERRQ(ierr); 2049 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE, "Must call SNESSetFunction() or SNESSetDM() before SNESComputeFunction(), likely called from SNESSolve()."); 2050 if (snes->vec_rhs) { 2051 ierr = VecAXPY(y,-1.0,snes->vec_rhs);CHKERRQ(ierr); 2052 } 2053 snes->nfuncs++; 2054 ierr = VecValidValues(y,3,PETSC_FALSE);CHKERRQ(ierr); 2055 PetscFunctionReturn(0); 2056 } 2057 2058 #undef __FUNCT__ 2059 #define __FUNCT__ "SNESComputeGS" 2060 /*@ 2061 SNESComputeGS - Calls the Gauss-Seidel function that has been set with SNESSetGS(). 2062 2063 Collective on SNES 2064 2065 Input Parameters: 2066 + snes - the SNES context 2067 . x - input vector 2068 - b - rhs vector 2069 2070 Output Parameter: 2071 . x - new solution vector 2072 2073 Notes: 2074 SNESComputeGS() is typically used within composed nonlinear solver 2075 implementations, so most users would not generally call this routine 2076 themselves. 2077 2078 Level: developer 2079 2080 .keywords: SNES, nonlinear, compute, function 2081 2082 .seealso: SNESSetGS(), SNESComputeFunction() 2083 @*/ 2084 PetscErrorCode SNESComputeGS(SNES snes,Vec b,Vec x) 2085 { 2086 PetscErrorCode ierr; 2087 DM dm; 2088 DMSNES sdm; 2089 2090 PetscFunctionBegin; 2091 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2092 PetscValidHeaderSpecific(x,VEC_CLASSID,2); 2093 if (b) PetscValidHeaderSpecific(b,VEC_CLASSID,3); 2094 PetscCheckSameComm(snes,1,x,2); 2095 if (b) PetscCheckSameComm(snes,1,b,3); 2096 if (b) {ierr = VecValidValues(b,2,PETSC_TRUE);CHKERRQ(ierr);} 2097 ierr = PetscLogEventBegin(SNES_GSEval,snes,x,b,0);CHKERRQ(ierr); 2098 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 2099 ierr = DMGetDMSNES(dm,&sdm);CHKERRQ(ierr); 2100 if (sdm->ops->computegs) { 2101 PetscStackPush("SNES user GS"); 2102 ierr = (*sdm->ops->computegs)(snes,x,b,sdm->gsctx);CHKERRQ(ierr); 2103 PetscStackPop; 2104 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE, "Must call SNESSetGS() before SNESComputeGS(), likely called from SNESSolve()."); 2105 ierr = PetscLogEventEnd(SNES_GSEval,snes,x,b,0);CHKERRQ(ierr); 2106 ierr = VecValidValues(x,3,PETSC_FALSE);CHKERRQ(ierr); 2107 PetscFunctionReturn(0); 2108 } 2109 2110 #undef __FUNCT__ 2111 #define __FUNCT__ "SNESComputeJacobian" 2112 /*@ 2113 SNESComputeJacobian - Computes the Jacobian matrix that has been set with SNESSetJacobian(). 2114 2115 Collective on SNES and Mat 2116 2117 Input Parameters: 2118 + snes - the SNES context 2119 - x - input vector 2120 2121 Output Parameters: 2122 + A - Jacobian matrix 2123 . B - optional preconditioning matrix 2124 - flag - flag indicating matrix structure (one of, SAME_NONZERO_PATTERN,DIFFERENT_NONZERO_PATTERN,SAME_PRECONDITIONER) 2125 2126 Options Database Keys: 2127 + -snes_lag_preconditioner <lag> 2128 . -snes_lag_jacobian <lag> 2129 . -snes_compare_explicit - Compare the computed Jacobian to the finite difference Jacobian and output the differences 2130 . -snes_compare_explicit_draw - Compare the computed Jacobian to the finite difference Jacobian and draw the result 2131 . -snes_compare_explicit_contour - Compare the computed Jacobian to the finite difference Jacobian and draw a contour plot with the result 2132 . -snes_compare_operator - Make the comparison options above use the operator instead of the preconditioning matrix 2133 . -snes_compare_coloring - Compute the finite differece Jacobian using coloring and display norms of difference 2134 . -snes_compare_coloring_display - Compute the finite differece Jacobian using coloring and display verbose differences 2135 . -snes_compare_coloring_threshold - Display only those matrix entries that differ by more than a given threshold 2136 . -snes_compare_coloring_threshold_atol - Absolute tolerance for difference in matrix entries to be displayed by -snes_compare_coloring_threshold 2137 . -snes_compare_coloring_threshold_rtol - Relative tolerance for difference in matrix entries to be displayed by -snes_compare_coloring_threshold 2138 . -snes_compare_coloring_draw - Compute the finite differece Jacobian using coloring and draw differences 2139 - -snes_compare_coloring_draw_contour - Compute the finite differece Jacobian using coloring and show contours of matrices and differences 2140 2141 2142 Notes: 2143 Most users should not need to explicitly call this routine, as it 2144 is used internally within the nonlinear solvers. 2145 2146 See KSPSetOperators() for important information about setting the 2147 flag parameter. 2148 2149 Level: developer 2150 2151 .keywords: SNES, compute, Jacobian, matrix 2152 2153 .seealso: SNESSetJacobian(), KSPSetOperators(), MatStructure, SNESSetLagPreconditioner(), SNESSetLagJacobian() 2154 @*/ 2155 PetscErrorCode SNESComputeJacobian(SNES snes,Vec X,Mat *A,Mat *B,MatStructure *flg) 2156 { 2157 PetscErrorCode ierr; 2158 PetscBool flag; 2159 DM dm; 2160 DMSNES sdm; 2161 2162 PetscFunctionBegin; 2163 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2164 PetscValidHeaderSpecific(X,VEC_CLASSID,2); 2165 PetscValidPointer(flg,5); 2166 PetscCheckSameComm(snes,1,X,2); 2167 ierr = VecValidValues(X,2,PETSC_TRUE);CHKERRQ(ierr); 2168 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 2169 ierr = DMGetDMSNES(dm,&sdm);CHKERRQ(ierr); 2170 2171 if (!sdm->ops->computejacobian) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_USER,"Must call SNESSetJacobian(), DMSNESSetJacobian(), DMDASNESSetJacobianLocal(), etc"); 2172 2173 /* make sure that MatAssemblyBegin/End() is called on A matrix if it is matrix free */ 2174 2175 if (snes->lagjacobian == -2) { 2176 snes->lagjacobian = -1; 2177 2178 ierr = PetscInfo(snes,"Recomputing Jacobian/preconditioner because lag is -2 (means compute Jacobian, but then never again) \n");CHKERRQ(ierr); 2179 } else if (snes->lagjacobian == -1) { 2180 *flg = SAME_PRECONDITIONER; 2181 ierr = PetscInfo(snes,"Reusing Jacobian/preconditioner because lag is -1\n");CHKERRQ(ierr); 2182 ierr = PetscObjectTypeCompare((PetscObject)*A,MATMFFD,&flag);CHKERRQ(ierr); 2183 if (flag) { 2184 ierr = MatAssemblyBegin(*A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2185 ierr = MatAssemblyEnd(*A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2186 } 2187 PetscFunctionReturn(0); 2188 } else if (snes->lagjacobian > 1 && snes->iter % snes->lagjacobian) { 2189 *flg = SAME_PRECONDITIONER; 2190 ierr = PetscInfo2(snes,"Reusing Jacobian/preconditioner because lag is %D and SNES iteration is %D\n",snes->lagjacobian,snes->iter);CHKERRQ(ierr); 2191 ierr = PetscObjectTypeCompare((PetscObject)*A,MATMFFD,&flag);CHKERRQ(ierr); 2192 if (flag) { 2193 ierr = MatAssemblyBegin(*A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2194 ierr = MatAssemblyEnd(*A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2195 } 2196 PetscFunctionReturn(0); 2197 } 2198 if (snes->pc && snes->pcside == PC_LEFT) { 2199 ierr = MatAssemblyBegin(*A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2200 ierr = MatAssemblyEnd(*A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2201 PetscFunctionReturn(0); 2202 } 2203 2204 *flg = DIFFERENT_NONZERO_PATTERN; 2205 ierr = PetscLogEventBegin(SNES_JacobianEval,snes,X,*A,*B);CHKERRQ(ierr); 2206 2207 PetscStackPush("SNES user Jacobian function"); 2208 ierr = (*sdm->ops->computejacobian)(snes,X,A,B,flg,sdm->jacobianctx);CHKERRQ(ierr); 2209 PetscStackPop; 2210 2211 ierr = PetscLogEventEnd(SNES_JacobianEval,snes,X,*A,*B);CHKERRQ(ierr); 2212 2213 if (snes->lagpreconditioner == -2) { 2214 ierr = PetscInfo(snes,"Rebuilding preconditioner exactly once since lag is -2\n");CHKERRQ(ierr); 2215 2216 snes->lagpreconditioner = -1; 2217 } else if (snes->lagpreconditioner == -1) { 2218 *flg = SAME_PRECONDITIONER; 2219 ierr = PetscInfo(snes,"Reusing preconditioner because lag is -1\n");CHKERRQ(ierr); 2220 } else if (snes->lagpreconditioner > 1 && snes->iter % snes->lagpreconditioner) { 2221 *flg = SAME_PRECONDITIONER; 2222 ierr = PetscInfo2(snes,"Reusing preconditioner because lag is %D and SNES iteration is %D\n",snes->lagpreconditioner,snes->iter);CHKERRQ(ierr); 2223 } 2224 2225 /* make sure user returned a correct Jacobian and preconditioner */ 2226 /* PetscValidHeaderSpecific(*A,MAT_CLASSID,3); 2227 PetscValidHeaderSpecific(*B,MAT_CLASSID,4); */ 2228 { 2229 PetscBool flag = PETSC_FALSE,flag_draw = PETSC_FALSE,flag_contour = PETSC_FALSE,flag_operator = PETSC_FALSE; 2230 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_explicit",&flag,NULL);CHKERRQ(ierr); 2231 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_explicit_draw",&flag_draw,NULL);CHKERRQ(ierr); 2232 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_explicit_draw_contour",&flag_contour,NULL);CHKERRQ(ierr); 2233 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_operator",&flag_operator,NULL);CHKERRQ(ierr); 2234 if (flag || flag_draw || flag_contour) { 2235 Mat Bexp_mine = NULL,Bexp,FDexp; 2236 MatStructure mstruct; 2237 PetscViewer vdraw,vstdout; 2238 PetscBool flg; 2239 if (flag_operator) { 2240 ierr = MatComputeExplicitOperator(*A,&Bexp_mine);CHKERRQ(ierr); 2241 Bexp = Bexp_mine; 2242 } else { 2243 /* See if the preconditioning matrix can be viewed and added directly */ 2244 ierr = PetscObjectTypeCompareAny((PetscObject)*B,&flg,MATSEQAIJ,MATMPIAIJ,MATSEQDENSE,MATMPIDENSE,MATSEQBAIJ,MATMPIBAIJ,MATSEQSBAIJ,MATMPIBAIJ,"");CHKERRQ(ierr); 2245 if (flg) Bexp = *B; 2246 else { 2247 /* If the "preconditioning" matrix is itself MATSHELL or some other type without direct support */ 2248 ierr = MatComputeExplicitOperator(*B,&Bexp_mine);CHKERRQ(ierr); 2249 Bexp = Bexp_mine; 2250 } 2251 } 2252 ierr = MatConvert(Bexp,MATSAME,MAT_INITIAL_MATRIX,&FDexp);CHKERRQ(ierr); 2253 ierr = SNESComputeJacobianDefault(snes,X,&FDexp,&FDexp,&mstruct,NULL);CHKERRQ(ierr); 2254 ierr = PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)snes),&vstdout);CHKERRQ(ierr); 2255 if (flag_draw || flag_contour) { 2256 ierr = PetscViewerDrawOpen(PetscObjectComm((PetscObject)snes),0,"Explicit Jacobians",PETSC_DECIDE,PETSC_DECIDE,300,300,&vdraw);CHKERRQ(ierr); 2257 if (flag_contour) {ierr = PetscViewerPushFormat(vdraw,PETSC_VIEWER_DRAW_CONTOUR);CHKERRQ(ierr);} 2258 } else vdraw = NULL; 2259 ierr = PetscViewerASCIIPrintf(vstdout,"Explicit %s\n",flag_operator ? "Jacobian" : "preconditioning Jacobian");CHKERRQ(ierr); 2260 if (flag) {ierr = MatView(Bexp,vstdout);CHKERRQ(ierr);} 2261 if (vdraw) {ierr = MatView(Bexp,vdraw);CHKERRQ(ierr);} 2262 ierr = PetscViewerASCIIPrintf(vstdout,"Finite difference Jacobian\n");CHKERRQ(ierr); 2263 if (flag) {ierr = MatView(FDexp,vstdout);CHKERRQ(ierr);} 2264 if (vdraw) {ierr = MatView(FDexp,vdraw);CHKERRQ(ierr);} 2265 ierr = MatAYPX(FDexp,-1.0,Bexp,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 2266 ierr = PetscViewerASCIIPrintf(vstdout,"User-provided matrix minus finite difference Jacobian\n");CHKERRQ(ierr); 2267 if (flag) {ierr = MatView(FDexp,vstdout);CHKERRQ(ierr);} 2268 if (vdraw) { /* Always use contour for the difference */ 2269 ierr = PetscViewerPushFormat(vdraw,PETSC_VIEWER_DRAW_CONTOUR);CHKERRQ(ierr); 2270 ierr = MatView(FDexp,vdraw);CHKERRQ(ierr); 2271 ierr = PetscViewerPopFormat(vdraw);CHKERRQ(ierr); 2272 } 2273 if (flag_contour) {ierr = PetscViewerPopFormat(vdraw);CHKERRQ(ierr);} 2274 ierr = PetscViewerDestroy(&vdraw);CHKERRQ(ierr); 2275 ierr = MatDestroy(&Bexp_mine);CHKERRQ(ierr); 2276 ierr = MatDestroy(&FDexp);CHKERRQ(ierr); 2277 } 2278 } 2279 { 2280 PetscBool flag = PETSC_FALSE,flag_display = PETSC_FALSE,flag_draw = PETSC_FALSE,flag_contour = PETSC_FALSE,flag_threshold = PETSC_FALSE; 2281 PetscReal threshold_atol = PETSC_SQRT_MACHINE_EPSILON,threshold_rtol = 10*PETSC_SQRT_MACHINE_EPSILON; 2282 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_coloring",&flag,NULL);CHKERRQ(ierr); 2283 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_coloring_display",&flag_display,NULL);CHKERRQ(ierr); 2284 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_coloring_draw",&flag_draw,NULL);CHKERRQ(ierr); 2285 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_coloring_draw_contour",&flag_contour,NULL);CHKERRQ(ierr); 2286 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_coloring_threshold",&flag_threshold,NULL);CHKERRQ(ierr); 2287 ierr = PetscOptionsGetReal(((PetscObject)snes)->prefix,"-snes_compare_coloring_threshold_rtol",&threshold_rtol,NULL);CHKERRQ(ierr); 2288 ierr = PetscOptionsGetReal(((PetscObject)snes)->prefix,"-snes_compare_coloring_threshold_atol",&threshold_atol,NULL);CHKERRQ(ierr); 2289 if (flag || flag_display || flag_draw || flag_contour || flag_threshold) { 2290 Mat Bfd; 2291 MatStructure mstruct; 2292 PetscViewer vdraw,vstdout; 2293 ISColoring iscoloring; 2294 MatFDColoring matfdcoloring; 2295 PetscErrorCode (*func)(SNES,Vec,Vec,void*); 2296 void *funcctx; 2297 PetscReal norm1,norm2,normmax; 2298 2299 ierr = MatDuplicate(*B,MAT_DO_NOT_COPY_VALUES,&Bfd);CHKERRQ(ierr); 2300 ierr = MatGetColoring(Bfd,MATCOLORINGSL,&iscoloring);CHKERRQ(ierr); 2301 ierr = MatFDColoringCreate(Bfd,iscoloring,&matfdcoloring);CHKERRQ(ierr); 2302 ierr = ISColoringDestroy(&iscoloring);CHKERRQ(ierr); 2303 2304 /* This method of getting the function is currently unreliable since it doesn't work for DM local functions. */ 2305 ierr = SNESGetFunction(snes,NULL,&func,&funcctx);CHKERRQ(ierr); 2306 ierr = MatFDColoringSetFunction(matfdcoloring,(PetscErrorCode (*)(void))func,funcctx);CHKERRQ(ierr); 2307 ierr = PetscObjectSetOptionsPrefix((PetscObject)matfdcoloring,((PetscObject)snes)->prefix);CHKERRQ(ierr); 2308 ierr = PetscObjectAppendOptionsPrefix((PetscObject)matfdcoloring,"coloring_");CHKERRQ(ierr); 2309 ierr = MatFDColoringSetFromOptions(matfdcoloring);CHKERRQ(ierr); 2310 ierr = MatFDColoringApply(Bfd,matfdcoloring,X,&mstruct,snes);CHKERRQ(ierr); 2311 ierr = MatFDColoringDestroy(&matfdcoloring);CHKERRQ(ierr); 2312 2313 ierr = PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)snes),&vstdout);CHKERRQ(ierr); 2314 if (flag_draw || flag_contour) { 2315 ierr = PetscViewerDrawOpen(PetscObjectComm((PetscObject)snes),0,"Colored Jacobians",PETSC_DECIDE,PETSC_DECIDE,300,300,&vdraw);CHKERRQ(ierr); 2316 if (flag_contour) {ierr = PetscViewerPushFormat(vdraw,PETSC_VIEWER_DRAW_CONTOUR);CHKERRQ(ierr);} 2317 } else vdraw = NULL; 2318 ierr = PetscViewerASCIIPrintf(vstdout,"Explicit preconditioning Jacobian\n");CHKERRQ(ierr); 2319 if (flag_display) {ierr = MatView(*B,vstdout);CHKERRQ(ierr);} 2320 if (vdraw) {ierr = MatView(*B,vdraw);CHKERRQ(ierr);} 2321 ierr = PetscViewerASCIIPrintf(vstdout,"Colored Finite difference Jacobian\n");CHKERRQ(ierr); 2322 if (flag_display) {ierr = MatView(Bfd,vstdout);CHKERRQ(ierr);} 2323 if (vdraw) {ierr = MatView(Bfd,vdraw);CHKERRQ(ierr);} 2324 ierr = MatAYPX(Bfd,-1.0,*B,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 2325 ierr = MatNorm(Bfd,NORM_1,&norm1);CHKERRQ(ierr); 2326 ierr = MatNorm(Bfd,NORM_FROBENIUS,&norm2);CHKERRQ(ierr); 2327 ierr = MatNorm(Bfd,NORM_MAX,&normmax);CHKERRQ(ierr); 2328 ierr = PetscViewerASCIIPrintf(vstdout,"User-provided matrix minus finite difference Jacobian, norm1=%G normFrob=%G normmax=%G\n",norm1,norm2,normmax);CHKERRQ(ierr); 2329 if (flag_display) {ierr = MatView(Bfd,vstdout);CHKERRQ(ierr);} 2330 if (vdraw) { /* Always use contour for the difference */ 2331 ierr = PetscViewerPushFormat(vdraw,PETSC_VIEWER_DRAW_CONTOUR);CHKERRQ(ierr); 2332 ierr = MatView(Bfd,vdraw);CHKERRQ(ierr); 2333 ierr = PetscViewerPopFormat(vdraw);CHKERRQ(ierr); 2334 } 2335 if (flag_contour) {ierr = PetscViewerPopFormat(vdraw);CHKERRQ(ierr);} 2336 2337 if (flag_threshold) { 2338 PetscInt bs,rstart,rend,i; 2339 ierr = MatGetBlockSize(*B,&bs);CHKERRQ(ierr); 2340 ierr = MatGetOwnershipRange(*B,&rstart,&rend);CHKERRQ(ierr); 2341 for (i=rstart; i<rend; i++) { 2342 const PetscScalar *ba,*ca; 2343 const PetscInt *bj,*cj; 2344 PetscInt bn,cn,j,maxentrycol = -1,maxdiffcol = -1,maxrdiffcol = -1; 2345 PetscReal maxentry = 0,maxdiff = 0,maxrdiff = 0; 2346 ierr = MatGetRow(*B,i,&bn,&bj,&ba);CHKERRQ(ierr); 2347 ierr = MatGetRow(Bfd,i,&cn,&cj,&ca);CHKERRQ(ierr); 2348 if (bn != cn) SETERRQ(((PetscObject)*A)->comm,PETSC_ERR_PLIB,"Unexpected different nonzero pattern in -snes_compare_coloring_threshold"); 2349 for (j=0; j<bn; j++) { 2350 PetscReal rdiff = PetscAbsScalar(ca[j]) / (threshold_atol + threshold_rtol*PetscAbsScalar(ba[j])); 2351 if (PetscAbsScalar(ba[j]) > PetscAbs(maxentry)) { 2352 maxentrycol = bj[j]; 2353 maxentry = PetscRealPart(ba[j]); 2354 } 2355 if (PetscAbsScalar(ca[j]) > PetscAbs(maxdiff)) { 2356 maxdiffcol = bj[j]; 2357 maxdiff = PetscRealPart(ca[j]); 2358 } 2359 if (rdiff > maxrdiff) { 2360 maxrdiffcol = bj[j]; 2361 maxrdiff = rdiff; 2362 } 2363 } 2364 if (maxrdiff > 1) { 2365 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); 2366 for (j=0; j<bn; j++) { 2367 PetscReal rdiff; 2368 rdiff = PetscAbsScalar(ca[j]) / (threshold_atol + threshold_rtol*PetscAbsScalar(ba[j])); 2369 if (rdiff > 1) { 2370 ierr = PetscViewerASCIIPrintf(vstdout," (%D,%G:%G)",bj[j],PetscRealPart(ba[j]),PetscRealPart(ca[j]));CHKERRQ(ierr); 2371 } 2372 } 2373 ierr = PetscViewerASCIIPrintf(vstdout,"\n",i,maxentry,maxdiff,maxrdiff);CHKERRQ(ierr); 2374 } 2375 ierr = MatRestoreRow(*B,i,&bn,&bj,&ba);CHKERRQ(ierr); 2376 ierr = MatRestoreRow(Bfd,i,&cn,&cj,&ca);CHKERRQ(ierr); 2377 } 2378 } 2379 ierr = PetscViewerDestroy(&vdraw);CHKERRQ(ierr); 2380 ierr = MatDestroy(&Bfd);CHKERRQ(ierr); 2381 } 2382 } 2383 PetscFunctionReturn(0); 2384 } 2385 2386 /*MC 2387 SNESJacobianFunction - function used to convey the nonlinear Jacobian of the function to be solved by SNES 2388 2389 Synopsis: 2390 #include "petscsnes.h" 2391 $ SNESJacobianFunction(SNES snes,Vec x,Mat *Amat,Mat *Pmat,int *flag,void *ctx); 2392 2393 + x - input vector 2394 . Amat - the matrix that defines the (approximate) Jacobian 2395 . Pmat - the matrix to be used in constructing the preconditioner, usually the same as Amat. 2396 . flag - flag indicating information about the preconditioner matrix 2397 structure (same as flag in KSPSetOperators()), one of SAME_NONZERO_PATTERN,DIFFERENT_NONZERO_PATTERN,SAME_PRECONDITIONER 2398 - ctx - [optional] user-defined Jacobian context 2399 2400 Level: intermediate 2401 2402 .seealso: SNESSetFunction(), SNESGetFunction(), SNESSetJacobian(), SNESGetJacobian() 2403 M*/ 2404 2405 #undef __FUNCT__ 2406 #define __FUNCT__ "SNESSetJacobian" 2407 /*@C 2408 SNESSetJacobian - Sets the function to compute Jacobian as well as the 2409 location to store the matrix. 2410 2411 Logically Collective on SNES and Mat 2412 2413 Input Parameters: 2414 + snes - the SNES context 2415 . Amat - the matrix that defines the (approximate) Jacobian 2416 . Pmat - the matrix to be used in constructing the preconditioner, usually the same as Amat. 2417 . SNESJacobianFunction - Jacobian evaluation routine (if NULL then SNES retains any previously set value) 2418 - ctx - [optional] user-defined context for private data for the 2419 Jacobian evaluation routine (may be NULL) (if NULL then SNES retains any previously set value) 2420 2421 Notes: 2422 See KSPSetOperators() for important information about setting the flag 2423 output parameter in the routine func(). Be sure to read this information! 2424 2425 The routine func() takes Mat * as the matrix arguments rather than Mat. 2426 This allows the Jacobian evaluation routine to replace A and/or B with a 2427 completely new new matrix structure (not just different matrix elements) 2428 when appropriate, for instance, if the nonzero structure is changing 2429 throughout the global iterations. 2430 2431 If the Amat matrix and Pmat matrix are different you must call MatAssemblyBegin/End() on 2432 each matrix. 2433 2434 If using SNESComputeJacobianDefaultColor() to assemble a Jacobian, the ctx argument 2435 must be a MatFDColoring. 2436 2437 Other defect-correction schemes can be used by computing a different matrix in place of the Jacobian. One common 2438 example is to use the "Picard linearization" which only differentiates through the highest order parts of each term. 2439 2440 Level: beginner 2441 2442 .keywords: SNES, nonlinear, set, Jacobian, matrix 2443 2444 .seealso: KSPSetOperators(), SNESSetFunction(), MatMFFDComputeJacobian(), SNESComputeJacobianDefaultColor(), MatStructure, SNESJacobianFunction, SNESSetPicard() 2445 @*/ 2446 PetscErrorCode SNESSetJacobian(SNES snes,Mat Amat,Mat Pmat,PetscErrorCode (*SNESJacobianFunction)(SNES,Vec,Mat*,Mat*,MatStructure*,void*),void *ctx) 2447 { 2448 PetscErrorCode ierr; 2449 DM dm; 2450 2451 PetscFunctionBegin; 2452 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2453 if (Amat) PetscValidHeaderSpecific(Amat,MAT_CLASSID,2); 2454 if (Pmat) PetscValidHeaderSpecific(Pmat,MAT_CLASSID,3); 2455 if (Amat) PetscCheckSameComm(snes,1,Amat,2); 2456 if (Pmat) PetscCheckSameComm(snes,1,Pmat,3); 2457 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 2458 ierr = DMSNESSetJacobian(dm,SNESJacobianFunction,ctx);CHKERRQ(ierr); 2459 if (Amat) { 2460 ierr = PetscObjectReference((PetscObject)Amat);CHKERRQ(ierr); 2461 ierr = MatDestroy(&snes->jacobian);CHKERRQ(ierr); 2462 2463 snes->jacobian = Amat; 2464 } 2465 if (Pmat) { 2466 ierr = PetscObjectReference((PetscObject)Pmat);CHKERRQ(ierr); 2467 ierr = MatDestroy(&snes->jacobian_pre);CHKERRQ(ierr); 2468 2469 snes->jacobian_pre = Pmat; 2470 } 2471 PetscFunctionReturn(0); 2472 } 2473 2474 #undef __FUNCT__ 2475 #define __FUNCT__ "SNESGetJacobian" 2476 /*@C 2477 SNESGetJacobian - Returns the Jacobian matrix and optionally the user 2478 provided context for evaluating the Jacobian. 2479 2480 Not Collective, but Mat object will be parallel if SNES object is 2481 2482 Input Parameter: 2483 . snes - the nonlinear solver context 2484 2485 Output Parameters: 2486 + Amat - location to stash (approximate) Jacobian matrix (or NULL) 2487 . Pmat - location to stash matrix used to compute the preconditioner (or NULL) 2488 . SNESJacobianFunction - location to put Jacobian function (or NULL) 2489 - ctx - location to stash Jacobian ctx (or NULL) 2490 2491 Level: advanced 2492 2493 .seealso: SNESSetJacobian(), SNESComputeJacobian() 2494 @*/ 2495 PetscErrorCode SNESGetJacobian(SNES snes,Mat *Amat,Mat *Pmat,PetscErrorCode (**SNESJacobianFunction)(SNES,Vec,Mat*,Mat*,MatStructure*,void*),void **ctx) 2496 { 2497 PetscErrorCode ierr; 2498 DM dm; 2499 DMSNES sdm; 2500 2501 PetscFunctionBegin; 2502 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2503 if (Amat) *Amat = snes->jacobian; 2504 if (Pmat) *Pmat = snes->jacobian_pre; 2505 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 2506 ierr = DMGetDMSNES(dm,&sdm);CHKERRQ(ierr); 2507 if (SNESJacobianFunction) *SNESJacobianFunction = sdm->ops->computejacobian; 2508 if (ctx) *ctx = sdm->jacobianctx; 2509 PetscFunctionReturn(0); 2510 } 2511 2512 #undef __FUNCT__ 2513 #define __FUNCT__ "SNESSetUp" 2514 /*@ 2515 SNESSetUp - Sets up the internal data structures for the later use 2516 of a nonlinear solver. 2517 2518 Collective on SNES 2519 2520 Input Parameters: 2521 . snes - the SNES context 2522 2523 Notes: 2524 For basic use of the SNES solvers the user need not explicitly call 2525 SNESSetUp(), since these actions will automatically occur during 2526 the call to SNESSolve(). However, if one wishes to control this 2527 phase separately, SNESSetUp() should be called after SNESCreate() 2528 and optional routines of the form SNESSetXXX(), but before SNESSolve(). 2529 2530 Level: advanced 2531 2532 .keywords: SNES, nonlinear, setup 2533 2534 .seealso: SNESCreate(), SNESSolve(), SNESDestroy() 2535 @*/ 2536 PetscErrorCode SNESSetUp(SNES snes) 2537 { 2538 PetscErrorCode ierr; 2539 DM dm; 2540 DMSNES sdm; 2541 SNESLineSearch linesearch, pclinesearch; 2542 void *lsprectx,*lspostctx; 2543 PetscErrorCode (*precheck)(SNESLineSearch,Vec,Vec,PetscBool*,void*); 2544 PetscErrorCode (*postcheck)(SNESLineSearch,Vec,Vec,Vec,PetscBool*,PetscBool*,void*); 2545 PetscErrorCode (*func)(SNES,Vec,Vec,void*); 2546 Vec f,fpc; 2547 void *funcctx; 2548 PetscErrorCode (*jac)(SNES,Vec,Mat*,Mat*,MatStructure*,void*); 2549 void *jacctx,*appctx; 2550 Mat A,B; 2551 2552 PetscFunctionBegin; 2553 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2554 if (snes->setupcalled) PetscFunctionReturn(0); 2555 2556 if (!((PetscObject)snes)->type_name) { 2557 ierr = SNESSetType(snes,SNESNEWTONLS);CHKERRQ(ierr); 2558 } 2559 2560 ierr = SNESGetFunction(snes,&snes->vec_func,NULL,NULL);CHKERRQ(ierr); 2561 if (snes->vec_func == snes->vec_sol) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_IDN,"Solution vector cannot be function vector"); 2562 if (snes->vec_rhs == snes->vec_sol) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_IDN,"Solution vector cannot be right hand side vector"); 2563 2564 if (!snes->vec_sol_update /* && snes->vec_sol */) { 2565 ierr = VecDuplicate(snes->vec_sol,&snes->vec_sol_update);CHKERRQ(ierr); 2566 ierr = PetscLogObjectParent(snes,snes->vec_sol_update);CHKERRQ(ierr); 2567 } 2568 2569 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 2570 ierr = DMShellSetGlobalVector(snes->dm,snes->vec_sol);CHKERRQ(ierr); 2571 ierr = DMGetDMSNES(dm,&sdm);CHKERRQ(ierr); 2572 if (!sdm->ops->computefunction) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"Function never provided to SNES object"); 2573 if (!sdm->ops->computejacobian) { 2574 ierr = DMSNESSetJacobian(dm,SNESComputeJacobianDefaultColor,NULL);CHKERRQ(ierr); 2575 } 2576 if (!snes->vec_func) { 2577 ierr = DMCreateGlobalVector(dm,&snes->vec_func);CHKERRQ(ierr); 2578 } 2579 2580 if (!snes->ksp) { 2581 ierr = SNESGetKSP(snes, &snes->ksp);CHKERRQ(ierr); 2582 } 2583 2584 if (!snes->linesearch) { 2585 ierr = SNESGetLineSearch(snes, &snes->linesearch);CHKERRQ(ierr); 2586 } 2587 2588 if (snes->pc && (snes->pcside == PC_LEFT)) { 2589 snes->mf = PETSC_TRUE; 2590 snes->mf_operator = PETSC_FALSE; 2591 } 2592 2593 if (snes->mf) { 2594 ierr = SNESSetUpMatrixFree_Private(snes, snes->mf_operator, snes->mf_version);CHKERRQ(ierr); 2595 } 2596 2597 if (snes->ops->usercompute && !snes->user) { 2598 ierr = (*snes->ops->usercompute)(snes,(void**)&snes->user);CHKERRQ(ierr); 2599 } 2600 2601 if (snes->pc) { 2602 /* copy the DM over */ 2603 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 2604 ierr = SNESSetDM(snes->pc,dm);CHKERRQ(ierr); 2605 2606 ierr = SNESGetFunction(snes,&f,&func,&funcctx);CHKERRQ(ierr); 2607 ierr = VecDuplicate(f,&fpc);CHKERRQ(ierr); 2608 ierr = SNESSetFunction(snes->pc,fpc,func,funcctx);CHKERRQ(ierr); 2609 ierr = SNESGetJacobian(snes,&A,&B,&jac,&jacctx);CHKERRQ(ierr); 2610 ierr = SNESSetJacobian(snes->pc,A,B,jac,jacctx);CHKERRQ(ierr); 2611 ierr = SNESGetApplicationContext(snes,&appctx);CHKERRQ(ierr); 2612 ierr = SNESSetApplicationContext(snes->pc,appctx);CHKERRQ(ierr); 2613 ierr = VecDestroy(&fpc);CHKERRQ(ierr); 2614 2615 /* copy the function pointers over */ 2616 ierr = PetscObjectCopyFortranFunctionPointers((PetscObject)snes,(PetscObject)snes->pc);CHKERRQ(ierr); 2617 2618 /* default to 1 iteration */ 2619 ierr = SNESSetTolerances(snes->pc,0.0,0.0,0.0,1,snes->pc->max_funcs);CHKERRQ(ierr); 2620 if (snes->pcside==PC_RIGHT) { 2621 ierr = SNESSetNormSchedule(snes->pc,SNES_NORM_FINAL_ONLY);CHKERRQ(ierr); 2622 } else { 2623 ierr = SNESSetNormSchedule(snes->pc,SNES_NORM_NONE);CHKERRQ(ierr); 2624 } 2625 ierr = SNESSetFromOptions(snes->pc);CHKERRQ(ierr); 2626 2627 /* copy the line search context over */ 2628 ierr = SNESGetLineSearch(snes,&linesearch);CHKERRQ(ierr); 2629 ierr = SNESGetLineSearch(snes->pc,&pclinesearch);CHKERRQ(ierr); 2630 ierr = SNESLineSearchGetPreCheck(linesearch,&precheck,&lsprectx);CHKERRQ(ierr); 2631 ierr = SNESLineSearchGetPostCheck(linesearch,&postcheck,&lspostctx);CHKERRQ(ierr); 2632 ierr = SNESLineSearchSetPreCheck(pclinesearch,precheck,lsprectx);CHKERRQ(ierr); 2633 ierr = SNESLineSearchSetPostCheck(pclinesearch,postcheck,lspostctx);CHKERRQ(ierr); 2634 ierr = PetscObjectCopyFortranFunctionPointers((PetscObject)linesearch, (PetscObject)pclinesearch);CHKERRQ(ierr); 2635 } 2636 2637 if (snes->ops->setup) { 2638 ierr = (*snes->ops->setup)(snes);CHKERRQ(ierr); 2639 } 2640 2641 snes->setupcalled = PETSC_TRUE; 2642 PetscFunctionReturn(0); 2643 } 2644 2645 #undef __FUNCT__ 2646 #define __FUNCT__ "SNESReset" 2647 /*@ 2648 SNESReset - Resets a SNES context to the snessetupcalled = 0 state and removes any allocated Vecs and Mats 2649 2650 Collective on SNES 2651 2652 Input Parameter: 2653 . snes - iterative context obtained from SNESCreate() 2654 2655 Level: intermediate 2656 2657 Notes: Also calls the application context destroy routine set with SNESSetComputeApplicationContext() 2658 2659 .keywords: SNES, destroy 2660 2661 .seealso: SNESCreate(), SNESSetUp(), SNESSolve() 2662 @*/ 2663 PetscErrorCode SNESReset(SNES snes) 2664 { 2665 PetscErrorCode ierr; 2666 2667 PetscFunctionBegin; 2668 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2669 if (snes->ops->userdestroy && snes->user) { 2670 ierr = (*snes->ops->userdestroy)((void**)&snes->user);CHKERRQ(ierr); 2671 snes->user = NULL; 2672 } 2673 if (snes->pc) { 2674 ierr = SNESReset(snes->pc);CHKERRQ(ierr); 2675 } 2676 2677 if (snes->ops->reset) { 2678 ierr = (*snes->ops->reset)(snes);CHKERRQ(ierr); 2679 } 2680 if (snes->ksp) { 2681 ierr = KSPReset(snes->ksp);CHKERRQ(ierr); 2682 } 2683 2684 if (snes->linesearch) { 2685 ierr = SNESLineSearchReset(snes->linesearch);CHKERRQ(ierr); 2686 } 2687 2688 ierr = VecDestroy(&snes->vec_rhs);CHKERRQ(ierr); 2689 ierr = VecDestroy(&snes->vec_sol);CHKERRQ(ierr); 2690 ierr = VecDestroy(&snes->vec_sol_update);CHKERRQ(ierr); 2691 ierr = VecDestroy(&snes->vec_func);CHKERRQ(ierr); 2692 ierr = MatDestroy(&snes->jacobian);CHKERRQ(ierr); 2693 ierr = MatDestroy(&snes->jacobian_pre);CHKERRQ(ierr); 2694 ierr = VecDestroyVecs(snes->nwork,&snes->work);CHKERRQ(ierr); 2695 ierr = VecDestroyVecs(snes->nvwork,&snes->vwork);CHKERRQ(ierr); 2696 2697 snes->nwork = snes->nvwork = 0; 2698 snes->setupcalled = PETSC_FALSE; 2699 PetscFunctionReturn(0); 2700 } 2701 2702 #undef __FUNCT__ 2703 #define __FUNCT__ "SNESDestroy" 2704 /*@ 2705 SNESDestroy - Destroys the nonlinear solver context that was created 2706 with SNESCreate(). 2707 2708 Collective on SNES 2709 2710 Input Parameter: 2711 . snes - the SNES context 2712 2713 Level: beginner 2714 2715 .keywords: SNES, nonlinear, destroy 2716 2717 .seealso: SNESCreate(), SNESSolve() 2718 @*/ 2719 PetscErrorCode SNESDestroy(SNES *snes) 2720 { 2721 PetscErrorCode ierr; 2722 2723 PetscFunctionBegin; 2724 if (!*snes) PetscFunctionReturn(0); 2725 PetscValidHeaderSpecific((*snes),SNES_CLASSID,1); 2726 if (--((PetscObject)(*snes))->refct > 0) {*snes = 0; PetscFunctionReturn(0);} 2727 2728 ierr = SNESReset((*snes));CHKERRQ(ierr); 2729 ierr = SNESDestroy(&(*snes)->pc);CHKERRQ(ierr); 2730 2731 /* if memory was published with AMS then destroy it */ 2732 ierr = PetscObjectAMSViewOff((PetscObject)*snes);CHKERRQ(ierr); 2733 if ((*snes)->ops->destroy) {ierr = (*((*snes))->ops->destroy)((*snes));CHKERRQ(ierr);} 2734 2735 ierr = DMDestroy(&(*snes)->dm);CHKERRQ(ierr); 2736 ierr = KSPDestroy(&(*snes)->ksp);CHKERRQ(ierr); 2737 ierr = SNESLineSearchDestroy(&(*snes)->linesearch);CHKERRQ(ierr); 2738 2739 ierr = PetscFree((*snes)->kspconvctx);CHKERRQ(ierr); 2740 if ((*snes)->ops->convergeddestroy) { 2741 ierr = (*(*snes)->ops->convergeddestroy)((*snes)->cnvP);CHKERRQ(ierr); 2742 } 2743 if ((*snes)->conv_malloc) { 2744 ierr = PetscFree((*snes)->conv_hist);CHKERRQ(ierr); 2745 ierr = PetscFree((*snes)->conv_hist_its);CHKERRQ(ierr); 2746 } 2747 ierr = SNESMonitorCancel((*snes));CHKERRQ(ierr); 2748 ierr = PetscHeaderDestroy(snes);CHKERRQ(ierr); 2749 PetscFunctionReturn(0); 2750 } 2751 2752 /* ----------- Routines to set solver parameters ---------- */ 2753 2754 #undef __FUNCT__ 2755 #define __FUNCT__ "SNESSetLagPreconditioner" 2756 /*@ 2757 SNESSetLagPreconditioner - Determines when the preconditioner is rebuilt in the nonlinear solve. 2758 2759 Logically Collective on SNES 2760 2761 Input Parameters: 2762 + snes - the SNES context 2763 - lag - -1 indicates NEVER rebuild, 1 means rebuild every time the Jacobian is computed within a single nonlinear solve, 2 means every second time 2764 the Jacobian is built etc. -2 indicates rebuild preconditioner at next chance but then never rebuild after that 2765 2766 Options Database Keys: 2767 . -snes_lag_preconditioner <lag> 2768 2769 Notes: 2770 The default is 1 2771 The preconditioner is ALWAYS built in the first iteration of a nonlinear solve unless lag is -1 2772 If -1 is used before the very first nonlinear solve the preconditioner is still built because there is no previous preconditioner to use 2773 2774 Level: intermediate 2775 2776 .keywords: SNES, nonlinear, set, convergence, tolerances 2777 2778 .seealso: SNESSetTrustRegionTolerance(), SNESGetLagPreconditioner(), SNESSetLagJacobian(), SNESGetLagJacobian() 2779 2780 @*/ 2781 PetscErrorCode SNESSetLagPreconditioner(SNES snes,PetscInt lag) 2782 { 2783 PetscFunctionBegin; 2784 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2785 if (lag < -2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Lag must be -2, -1, 1 or greater"); 2786 if (!lag) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Lag cannot be 0"); 2787 PetscValidLogicalCollectiveInt(snes,lag,2); 2788 snes->lagpreconditioner = lag; 2789 PetscFunctionReturn(0); 2790 } 2791 2792 #undef __FUNCT__ 2793 #define __FUNCT__ "SNESSetGridSequence" 2794 /*@ 2795 SNESSetGridSequence - sets the number of steps of grid sequencing that SNES does 2796 2797 Logically Collective on SNES 2798 2799 Input Parameters: 2800 + snes - the SNES context 2801 - steps - the number of refinements to do, defaults to 0 2802 2803 Options Database Keys: 2804 . -snes_grid_sequence <steps> 2805 2806 Level: intermediate 2807 2808 Notes: 2809 Use SNESGetSolution() to extract the fine grid solution after grid sequencing. 2810 2811 .keywords: SNES, nonlinear, set, convergence, tolerances 2812 2813 .seealso: SNESSetTrustRegionTolerance(), SNESGetLagPreconditioner(), SNESSetLagJacobian(), SNESGetLagJacobian() 2814 2815 @*/ 2816 PetscErrorCode SNESSetGridSequence(SNES snes,PetscInt steps) 2817 { 2818 PetscFunctionBegin; 2819 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2820 PetscValidLogicalCollectiveInt(snes,steps,2); 2821 snes->gridsequence = steps; 2822 PetscFunctionReturn(0); 2823 } 2824 2825 #undef __FUNCT__ 2826 #define __FUNCT__ "SNESGetLagPreconditioner" 2827 /*@ 2828 SNESGetLagPreconditioner - Indicates how often the preconditioner is rebuilt 2829 2830 Not Collective 2831 2832 Input Parameter: 2833 . snes - the SNES context 2834 2835 Output Parameter: 2836 . lag - -1 indicates NEVER rebuild, 1 means rebuild every time the Jacobian is computed within a single nonlinear solve, 2 means every second time 2837 the Jacobian is built etc. -2 indicates rebuild preconditioner at next chance but then never rebuild after that 2838 2839 Options Database Keys: 2840 . -snes_lag_preconditioner <lag> 2841 2842 Notes: 2843 The default is 1 2844 The preconditioner is ALWAYS built in the first iteration of a nonlinear solve unless lag is -1 2845 2846 Level: intermediate 2847 2848 .keywords: SNES, nonlinear, set, convergence, tolerances 2849 2850 .seealso: SNESSetTrustRegionTolerance(), SNESSetLagPreconditioner() 2851 2852 @*/ 2853 PetscErrorCode SNESGetLagPreconditioner(SNES snes,PetscInt *lag) 2854 { 2855 PetscFunctionBegin; 2856 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2857 *lag = snes->lagpreconditioner; 2858 PetscFunctionReturn(0); 2859 } 2860 2861 #undef __FUNCT__ 2862 #define __FUNCT__ "SNESSetLagJacobian" 2863 /*@ 2864 SNESSetLagJacobian - Determines when the Jacobian is rebuilt in the nonlinear solve. See SNESSetLagPreconditioner() for determining how 2865 often the preconditioner is rebuilt. 2866 2867 Logically Collective on SNES 2868 2869 Input Parameters: 2870 + snes - the SNES context 2871 - lag - -1 indicates NEVER rebuild, 1 means rebuild every time the Jacobian is computed within a single nonlinear solve, 2 means every second time 2872 the Jacobian is built etc. -2 means rebuild at next chance but then never again 2873 2874 Options Database Keys: 2875 . -snes_lag_jacobian <lag> 2876 2877 Notes: 2878 The default is 1 2879 The Jacobian is ALWAYS built in the first iteration of a nonlinear solve unless lag is -1 2880 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 2881 at the next Newton step but never again (unless it is reset to another value) 2882 2883 Level: intermediate 2884 2885 .keywords: SNES, nonlinear, set, convergence, tolerances 2886 2887 .seealso: SNESSetTrustRegionTolerance(), SNESGetLagPreconditioner(), SNESSetLagPreconditioner(), SNESGetLagJacobian() 2888 2889 @*/ 2890 PetscErrorCode SNESSetLagJacobian(SNES snes,PetscInt lag) 2891 { 2892 PetscFunctionBegin; 2893 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2894 if (lag < -2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Lag must be -2, -1, 1 or greater"); 2895 if (!lag) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Lag cannot be 0"); 2896 PetscValidLogicalCollectiveInt(snes,lag,2); 2897 snes->lagjacobian = lag; 2898 PetscFunctionReturn(0); 2899 } 2900 2901 #undef __FUNCT__ 2902 #define __FUNCT__ "SNESGetLagJacobian" 2903 /*@ 2904 SNESGetLagJacobian - Indicates how often the Jacobian is rebuilt. See SNESGetLagPreconditioner() to determine when the preconditioner is rebuilt 2905 2906 Not Collective 2907 2908 Input Parameter: 2909 . snes - the SNES context 2910 2911 Output Parameter: 2912 . lag - -1 indicates NEVER rebuild, 1 means rebuild every time the Jacobian is computed within a single nonlinear solve, 2 means every second time 2913 the Jacobian is built etc. 2914 2915 Options Database Keys: 2916 . -snes_lag_jacobian <lag> 2917 2918 Notes: 2919 The default is 1 2920 The jacobian is ALWAYS built in the first iteration of a nonlinear solve unless lag is -1 2921 2922 Level: intermediate 2923 2924 .keywords: SNES, nonlinear, set, convergence, tolerances 2925 2926 .seealso: SNESSetTrustRegionTolerance(), SNESSetLagJacobian(), SNESSetLagPreconditioner(), SNESGetLagPreconditioner() 2927 2928 @*/ 2929 PetscErrorCode SNESGetLagJacobian(SNES snes,PetscInt *lag) 2930 { 2931 PetscFunctionBegin; 2932 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2933 *lag = snes->lagjacobian; 2934 PetscFunctionReturn(0); 2935 } 2936 2937 #undef __FUNCT__ 2938 #define __FUNCT__ "SNESSetTolerances" 2939 /*@ 2940 SNESSetTolerances - Sets various parameters used in convergence tests. 2941 2942 Logically Collective on SNES 2943 2944 Input Parameters: 2945 + snes - the SNES context 2946 . abstol - absolute convergence tolerance 2947 . rtol - relative convergence tolerance 2948 . stol - convergence tolerance in terms of the norm of the change in the solution between steps, || delta x || < stol*|| x || 2949 . maxit - maximum number of iterations 2950 - maxf - maximum number of function evaluations 2951 2952 Options Database Keys: 2953 + -snes_atol <abstol> - Sets abstol 2954 . -snes_rtol <rtol> - Sets rtol 2955 . -snes_stol <stol> - Sets stol 2956 . -snes_max_it <maxit> - Sets maxit 2957 - -snes_max_funcs <maxf> - Sets maxf 2958 2959 Notes: 2960 The default maximum number of iterations is 50. 2961 The default maximum number of function evaluations is 1000. 2962 2963 Level: intermediate 2964 2965 .keywords: SNES, nonlinear, set, convergence, tolerances 2966 2967 .seealso: SNESSetTrustRegionTolerance() 2968 @*/ 2969 PetscErrorCode SNESSetTolerances(SNES snes,PetscReal abstol,PetscReal rtol,PetscReal stol,PetscInt maxit,PetscInt maxf) 2970 { 2971 PetscFunctionBegin; 2972 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2973 PetscValidLogicalCollectiveReal(snes,abstol,2); 2974 PetscValidLogicalCollectiveReal(snes,rtol,3); 2975 PetscValidLogicalCollectiveReal(snes,stol,4); 2976 PetscValidLogicalCollectiveInt(snes,maxit,5); 2977 PetscValidLogicalCollectiveInt(snes,maxf,6); 2978 2979 if (abstol != PETSC_DEFAULT) { 2980 if (abstol < 0.0) SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_OUTOFRANGE,"Absolute tolerance %G must be non-negative",abstol); 2981 snes->abstol = abstol; 2982 } 2983 if (rtol != PETSC_DEFAULT) { 2984 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); 2985 snes->rtol = rtol; 2986 } 2987 if (stol != PETSC_DEFAULT) { 2988 if (stol < 0.0) SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_OUTOFRANGE,"Step tolerance %G must be non-negative",stol); 2989 snes->stol = stol; 2990 } 2991 if (maxit != PETSC_DEFAULT) { 2992 if (maxit < 0) SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_OUTOFRANGE,"Maximum number of iterations %D must be non-negative",maxit); 2993 snes->max_its = maxit; 2994 } 2995 if (maxf != PETSC_DEFAULT) { 2996 if (maxf < 0) SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_OUTOFRANGE,"Maximum number of function evaluations %D must be non-negative",maxf); 2997 snes->max_funcs = maxf; 2998 } 2999 snes->tolerancesset = PETSC_TRUE; 3000 PetscFunctionReturn(0); 3001 } 3002 3003 #undef __FUNCT__ 3004 #define __FUNCT__ "SNESGetTolerances" 3005 /*@ 3006 SNESGetTolerances - Gets various parameters used in convergence tests. 3007 3008 Not Collective 3009 3010 Input Parameters: 3011 + snes - the SNES context 3012 . atol - absolute convergence tolerance 3013 . rtol - relative convergence tolerance 3014 . stol - convergence tolerance in terms of the norm 3015 of the change in the solution between steps 3016 . maxit - maximum number of iterations 3017 - maxf - maximum number of function evaluations 3018 3019 Notes: 3020 The user can specify NULL for any parameter that is not needed. 3021 3022 Level: intermediate 3023 3024 .keywords: SNES, nonlinear, get, convergence, tolerances 3025 3026 .seealso: SNESSetTolerances() 3027 @*/ 3028 PetscErrorCode SNESGetTolerances(SNES snes,PetscReal *atol,PetscReal *rtol,PetscReal *stol,PetscInt *maxit,PetscInt *maxf) 3029 { 3030 PetscFunctionBegin; 3031 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3032 if (atol) *atol = snes->abstol; 3033 if (rtol) *rtol = snes->rtol; 3034 if (stol) *stol = snes->stol; 3035 if (maxit) *maxit = snes->max_its; 3036 if (maxf) *maxf = snes->max_funcs; 3037 PetscFunctionReturn(0); 3038 } 3039 3040 #undef __FUNCT__ 3041 #define __FUNCT__ "SNESSetTrustRegionTolerance" 3042 /*@ 3043 SNESSetTrustRegionTolerance - Sets the trust region parameter tolerance. 3044 3045 Logically Collective on SNES 3046 3047 Input Parameters: 3048 + snes - the SNES context 3049 - tol - tolerance 3050 3051 Options Database Key: 3052 . -snes_trtol <tol> - Sets tol 3053 3054 Level: intermediate 3055 3056 .keywords: SNES, nonlinear, set, trust region, tolerance 3057 3058 .seealso: SNESSetTolerances() 3059 @*/ 3060 PetscErrorCode SNESSetTrustRegionTolerance(SNES snes,PetscReal tol) 3061 { 3062 PetscFunctionBegin; 3063 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3064 PetscValidLogicalCollectiveReal(snes,tol,2); 3065 snes->deltatol = tol; 3066 PetscFunctionReturn(0); 3067 } 3068 3069 /* 3070 Duplicate the lg monitors for SNES from KSP; for some reason with 3071 dynamic libraries things don't work under Sun4 if we just use 3072 macros instead of functions 3073 */ 3074 #undef __FUNCT__ 3075 #define __FUNCT__ "SNESMonitorLGResidualNorm" 3076 PetscErrorCode SNESMonitorLGResidualNorm(SNES snes,PetscInt it,PetscReal norm,void *ctx) 3077 { 3078 PetscErrorCode ierr; 3079 3080 PetscFunctionBegin; 3081 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3082 ierr = KSPMonitorLGResidualNorm((KSP)snes,it,norm,ctx);CHKERRQ(ierr); 3083 PetscFunctionReturn(0); 3084 } 3085 3086 #undef __FUNCT__ 3087 #define __FUNCT__ "SNESMonitorLGCreate" 3088 PetscErrorCode SNESMonitorLGCreate(const char host[],const char label[],int x,int y,int m,int n,PetscDrawLG *draw) 3089 { 3090 PetscErrorCode ierr; 3091 3092 PetscFunctionBegin; 3093 ierr = KSPMonitorLGResidualNormCreate(host,label,x,y,m,n,draw);CHKERRQ(ierr); 3094 PetscFunctionReturn(0); 3095 } 3096 3097 #undef __FUNCT__ 3098 #define __FUNCT__ "SNESMonitorLGDestroy" 3099 PetscErrorCode SNESMonitorLGDestroy(PetscDrawLG *draw) 3100 { 3101 PetscErrorCode ierr; 3102 3103 PetscFunctionBegin; 3104 ierr = KSPMonitorLGResidualNormDestroy(draw);CHKERRQ(ierr); 3105 PetscFunctionReturn(0); 3106 } 3107 3108 extern PetscErrorCode SNESMonitorRange_Private(SNES,PetscInt,PetscReal*); 3109 #undef __FUNCT__ 3110 #define __FUNCT__ "SNESMonitorLGRange" 3111 PetscErrorCode SNESMonitorLGRange(SNES snes,PetscInt n,PetscReal rnorm,void *monctx) 3112 { 3113 PetscDrawLG lg; 3114 PetscErrorCode ierr; 3115 PetscReal x,y,per; 3116 PetscViewer v = (PetscViewer)monctx; 3117 static PetscReal prev; /* should be in the context */ 3118 PetscDraw draw; 3119 3120 PetscFunctionBegin; 3121 ierr = PetscViewerDrawGetDrawLG(v,0,&lg);CHKERRQ(ierr); 3122 if (!n) {ierr = PetscDrawLGReset(lg);CHKERRQ(ierr);} 3123 ierr = PetscDrawLGGetDraw(lg,&draw);CHKERRQ(ierr); 3124 ierr = PetscDrawSetTitle(draw,"Residual norm");CHKERRQ(ierr); 3125 x = (PetscReal)n; 3126 if (rnorm > 0.0) y = log10(rnorm); 3127 else y = -15.0; 3128 ierr = PetscDrawLGAddPoint(lg,&x,&y);CHKERRQ(ierr); 3129 if (n < 20 || !(n % 5)) { 3130 ierr = PetscDrawLGDraw(lg);CHKERRQ(ierr); 3131 } 3132 3133 ierr = PetscViewerDrawGetDrawLG(v,1,&lg);CHKERRQ(ierr); 3134 if (!n) {ierr = PetscDrawLGReset(lg);CHKERRQ(ierr);} 3135 ierr = PetscDrawLGGetDraw(lg,&draw);CHKERRQ(ierr); 3136 ierr = PetscDrawSetTitle(draw,"% elemts > .2*max elemt");CHKERRQ(ierr); 3137 ierr = SNESMonitorRange_Private(snes,n,&per);CHKERRQ(ierr); 3138 x = (PetscReal)n; 3139 y = 100.0*per; 3140 ierr = PetscDrawLGAddPoint(lg,&x,&y);CHKERRQ(ierr); 3141 if (n < 20 || !(n % 5)) { 3142 ierr = PetscDrawLGDraw(lg);CHKERRQ(ierr); 3143 } 3144 3145 ierr = PetscViewerDrawGetDrawLG(v,2,&lg);CHKERRQ(ierr); 3146 if (!n) {prev = rnorm;ierr = PetscDrawLGReset(lg);CHKERRQ(ierr);} 3147 ierr = PetscDrawLGGetDraw(lg,&draw);CHKERRQ(ierr); 3148 ierr = PetscDrawSetTitle(draw,"(norm -oldnorm)/oldnorm");CHKERRQ(ierr); 3149 x = (PetscReal)n; 3150 y = (prev - rnorm)/prev; 3151 ierr = PetscDrawLGAddPoint(lg,&x,&y);CHKERRQ(ierr); 3152 if (n < 20 || !(n % 5)) { 3153 ierr = PetscDrawLGDraw(lg);CHKERRQ(ierr); 3154 } 3155 3156 ierr = PetscViewerDrawGetDrawLG(v,3,&lg);CHKERRQ(ierr); 3157 if (!n) {ierr = PetscDrawLGReset(lg);CHKERRQ(ierr);} 3158 ierr = PetscDrawLGGetDraw(lg,&draw);CHKERRQ(ierr); 3159 ierr = PetscDrawSetTitle(draw,"(norm -oldnorm)/oldnorm*(% > .2 max)");CHKERRQ(ierr); 3160 x = (PetscReal)n; 3161 y = (prev - rnorm)/(prev*per); 3162 if (n > 2) { /*skip initial crazy value */ 3163 ierr = PetscDrawLGAddPoint(lg,&x,&y);CHKERRQ(ierr); 3164 } 3165 if (n < 20 || !(n % 5)) { 3166 ierr = PetscDrawLGDraw(lg);CHKERRQ(ierr); 3167 } 3168 prev = rnorm; 3169 PetscFunctionReturn(0); 3170 } 3171 3172 #undef __FUNCT__ 3173 #define __FUNCT__ "SNESMonitor" 3174 /*@ 3175 SNESMonitor - runs the user provided monitor routines, if they exist 3176 3177 Collective on SNES 3178 3179 Input Parameters: 3180 + snes - nonlinear solver context obtained from SNESCreate() 3181 . iter - iteration number 3182 - rnorm - relative norm of the residual 3183 3184 Notes: 3185 This routine is called by the SNES implementations. 3186 It does not typically need to be called by the user. 3187 3188 Level: developer 3189 3190 .seealso: SNESMonitorSet() 3191 @*/ 3192 PetscErrorCode SNESMonitor(SNES snes,PetscInt iter,PetscReal rnorm) 3193 { 3194 PetscErrorCode ierr; 3195 PetscInt i,n = snes->numbermonitors; 3196 3197 PetscFunctionBegin; 3198 for (i=0; i<n; i++) { 3199 ierr = (*snes->monitor[i])(snes,iter,rnorm,snes->monitorcontext[i]);CHKERRQ(ierr); 3200 } 3201 PetscFunctionReturn(0); 3202 } 3203 3204 /* ------------ Routines to set performance monitoring options ----------- */ 3205 3206 /*MC 3207 SNESMonitorFunction - functional form passed to SNESMonitorSet() to monitor convergence of nonlinear solver 3208 3209 Synopsis: 3210 #include "petscsnes.h" 3211 $ PetscErrorCode SNESMonitorFunction(SNES snes,PetscInt its, PetscReal norm,void *mctx) 3212 3213 + snes - the SNES context 3214 . its - iteration number 3215 . norm - 2-norm function value (may be estimated) 3216 - mctx - [optional] monitoring context 3217 3218 Level: advanced 3219 3220 .seealso: SNESMonitorSet(), SNESMonitorGet() 3221 M*/ 3222 3223 #undef __FUNCT__ 3224 #define __FUNCT__ "SNESMonitorSet" 3225 /*@C 3226 SNESMonitorSet - Sets an ADDITIONAL function that is to be used at every 3227 iteration of the nonlinear solver to display the iteration's 3228 progress. 3229 3230 Logically Collective on SNES 3231 3232 Input Parameters: 3233 + snes - the SNES context 3234 . SNESMonitorFunction - monitoring routine 3235 . mctx - [optional] user-defined context for private data for the 3236 monitor routine (use NULL if no context is desired) 3237 - monitordestroy - [optional] routine that frees monitor context 3238 (may be NULL) 3239 3240 Options Database Keys: 3241 + -snes_monitor - sets SNESMonitorDefault() 3242 . -snes_monitor_lg_residualnorm - sets line graph monitor, 3243 uses SNESMonitorLGCreate() 3244 - -snes_monitor_cancel - cancels all monitors that have 3245 been hardwired into a code by 3246 calls to SNESMonitorSet(), but 3247 does not cancel those set via 3248 the options database. 3249 3250 Notes: 3251 Several different monitoring routines may be set by calling 3252 SNESMonitorSet() multiple times; all will be called in the 3253 order in which they were set. 3254 3255 Fortran notes: Only a single monitor function can be set for each SNES object 3256 3257 Level: intermediate 3258 3259 .keywords: SNES, nonlinear, set, monitor 3260 3261 .seealso: SNESMonitorDefault(), SNESMonitorCancel(), SNESMonitorFunction 3262 @*/ 3263 PetscErrorCode SNESMonitorSet(SNES snes,PetscErrorCode (*SNESMonitorFunction)(SNES,PetscInt,PetscReal,void*),void *mctx,PetscErrorCode (*monitordestroy)(void**)) 3264 { 3265 PetscInt i; 3266 PetscErrorCode ierr; 3267 3268 PetscFunctionBegin; 3269 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3270 if (snes->numbermonitors >= MAXSNESMONITORS) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Too many monitors set"); 3271 for (i=0; i<snes->numbermonitors;i++) { 3272 if (SNESMonitorFunction == snes->monitor[i] && monitordestroy == snes->monitordestroy[i] && mctx == snes->monitorcontext[i]) { 3273 if (monitordestroy) { 3274 ierr = (*monitordestroy)(&mctx);CHKERRQ(ierr); 3275 } 3276 PetscFunctionReturn(0); 3277 } 3278 } 3279 snes->monitor[snes->numbermonitors] = SNESMonitorFunction; 3280 snes->monitordestroy[snes->numbermonitors] = monitordestroy; 3281 snes->monitorcontext[snes->numbermonitors++] = (void*)mctx; 3282 PetscFunctionReturn(0); 3283 } 3284 3285 #undef __FUNCT__ 3286 #define __FUNCT__ "SNESMonitorCancel" 3287 /*@C 3288 SNESMonitorCancel - Clears all the monitor functions for a SNES object. 3289 3290 Logically Collective on SNES 3291 3292 Input Parameters: 3293 . snes - the SNES context 3294 3295 Options Database Key: 3296 . -snes_monitor_cancel - cancels all monitors that have been hardwired 3297 into a code by calls to SNESMonitorSet(), but does not cancel those 3298 set via the options database 3299 3300 Notes: 3301 There is no way to clear one specific monitor from a SNES object. 3302 3303 Level: intermediate 3304 3305 .keywords: SNES, nonlinear, set, monitor 3306 3307 .seealso: SNESMonitorDefault(), SNESMonitorSet() 3308 @*/ 3309 PetscErrorCode SNESMonitorCancel(SNES snes) 3310 { 3311 PetscErrorCode ierr; 3312 PetscInt i; 3313 3314 PetscFunctionBegin; 3315 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3316 for (i=0; i<snes->numbermonitors; i++) { 3317 if (snes->monitordestroy[i]) { 3318 ierr = (*snes->monitordestroy[i])(&snes->monitorcontext[i]);CHKERRQ(ierr); 3319 } 3320 } 3321 snes->numbermonitors = 0; 3322 PetscFunctionReturn(0); 3323 } 3324 3325 /*MC 3326 SNESConvergenceTestFunction - functional form used for testing of convergence of nonlinear solver 3327 3328 Synopsis: 3329 #include "petscsnes.h" 3330 $ PetscErrorCode SNESConvergenceTest(SNES snes,PetscInt it,PetscReal xnorm,PetscReal gnorm,PetscReal f,SNESConvergedReason *reason,void *cctx) 3331 3332 + snes - the SNES context 3333 . it - current iteration (0 is the first and is before any Newton step) 3334 . cctx - [optional] convergence context 3335 . reason - reason for convergence/divergence 3336 . xnorm - 2-norm of current iterate 3337 . gnorm - 2-norm of current step 3338 - f - 2-norm of function 3339 3340 Level: intermediate 3341 3342 .seealso: SNESSetConvergenceTest(), SNESGetConvergenceTest() 3343 M*/ 3344 3345 #undef __FUNCT__ 3346 #define __FUNCT__ "SNESSetConvergenceTest" 3347 /*@C 3348 SNESSetConvergenceTest - Sets the function that is to be used 3349 to test for convergence of the nonlinear iterative solution. 3350 3351 Logically Collective on SNES 3352 3353 Input Parameters: 3354 + snes - the SNES context 3355 . SNESConvergenceTestFunction - routine to test for convergence 3356 . cctx - [optional] context for private data for the convergence routine (may be NULL) 3357 - destroy - [optional] destructor for the context (may be NULL; NULL_FUNCTION in Fortran) 3358 3359 Level: advanced 3360 3361 .keywords: SNES, nonlinear, set, convergence, test 3362 3363 .seealso: SNESConvergedDefault(), SNESSkipConverged(), SNESConvergenceTestFunction 3364 @*/ 3365 PetscErrorCode SNESSetConvergenceTest(SNES snes,PetscErrorCode (*SNESConvergenceTestFunction)(SNES,PetscInt,PetscReal,PetscReal,PetscReal,SNESConvergedReason*,void*),void *cctx,PetscErrorCode (*destroy)(void*)) 3366 { 3367 PetscErrorCode ierr; 3368 3369 PetscFunctionBegin; 3370 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3371 if (!SNESConvergenceTestFunction) SNESConvergenceTestFunction = SNESSkipConverged; 3372 if (snes->ops->convergeddestroy) { 3373 ierr = (*snes->ops->convergeddestroy)(snes->cnvP);CHKERRQ(ierr); 3374 } 3375 snes->ops->converged = SNESConvergenceTestFunction; 3376 snes->ops->convergeddestroy = destroy; 3377 snes->cnvP = cctx; 3378 PetscFunctionReturn(0); 3379 } 3380 3381 #undef __FUNCT__ 3382 #define __FUNCT__ "SNESGetConvergedReason" 3383 /*@ 3384 SNESGetConvergedReason - Gets the reason the SNES iteration was stopped. 3385 3386 Not Collective 3387 3388 Input Parameter: 3389 . snes - the SNES context 3390 3391 Output Parameter: 3392 . reason - negative value indicates diverged, positive value converged, see SNESConvergedReason or the 3393 manual pages for the individual convergence tests for complete lists 3394 3395 Level: intermediate 3396 3397 Notes: Can only be called after the call the SNESSolve() is complete. 3398 3399 .keywords: SNES, nonlinear, set, convergence, test 3400 3401 .seealso: SNESSetConvergenceTest(), SNESConvergedReason 3402 @*/ 3403 PetscErrorCode SNESGetConvergedReason(SNES snes,SNESConvergedReason *reason) 3404 { 3405 PetscFunctionBegin; 3406 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3407 PetscValidPointer(reason,2); 3408 *reason = snes->reason; 3409 PetscFunctionReturn(0); 3410 } 3411 3412 #undef __FUNCT__ 3413 #define __FUNCT__ "SNESSetConvergenceHistory" 3414 /*@ 3415 SNESSetConvergenceHistory - Sets the array used to hold the convergence history. 3416 3417 Logically Collective on SNES 3418 3419 Input Parameters: 3420 + snes - iterative context obtained from SNESCreate() 3421 . a - array to hold history, this array will contain the function norms computed at each step 3422 . its - integer array holds the number of linear iterations for each solve. 3423 . na - size of a and its 3424 - reset - PETSC_TRUE indicates each new nonlinear solve resets the history counter to zero, 3425 else it continues storing new values for new nonlinear solves after the old ones 3426 3427 Notes: 3428 If 'a' and 'its' are NULL then space is allocated for the history. If 'na' PETSC_DECIDE or PETSC_DEFAULT then a 3429 default array of length 10000 is allocated. 3430 3431 This routine is useful, e.g., when running a code for purposes 3432 of accurate performance monitoring, when no I/O should be done 3433 during the section of code that is being timed. 3434 3435 Level: intermediate 3436 3437 .keywords: SNES, set, convergence, history 3438 3439 .seealso: SNESGetConvergenceHistory() 3440 3441 @*/ 3442 PetscErrorCode SNESSetConvergenceHistory(SNES snes,PetscReal a[],PetscInt its[],PetscInt na,PetscBool reset) 3443 { 3444 PetscErrorCode ierr; 3445 3446 PetscFunctionBegin; 3447 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3448 if (a) PetscValidScalarPointer(a,2); 3449 if (its) PetscValidIntPointer(its,3); 3450 if (!a) { 3451 if (na == PETSC_DECIDE || na == PETSC_DEFAULT) na = 1000; 3452 ierr = PetscMalloc(na*sizeof(PetscReal),&a);CHKERRQ(ierr); 3453 ierr = PetscMalloc(na*sizeof(PetscInt),&its);CHKERRQ(ierr); 3454 3455 snes->conv_malloc = PETSC_TRUE; 3456 } 3457 snes->conv_hist = a; 3458 snes->conv_hist_its = its; 3459 snes->conv_hist_max = na; 3460 snes->conv_hist_len = 0; 3461 snes->conv_hist_reset = reset; 3462 PetscFunctionReturn(0); 3463 } 3464 3465 #if defined(PETSC_HAVE_MATLAB_ENGINE) 3466 #include <engine.h> /* MATLAB include file */ 3467 #include <mex.h> /* MATLAB include file */ 3468 3469 #undef __FUNCT__ 3470 #define __FUNCT__ "SNESGetConvergenceHistoryMatlab" 3471 PETSC_EXTERN mxArray *SNESGetConvergenceHistoryMatlab(SNES snes) 3472 { 3473 mxArray *mat; 3474 PetscInt i; 3475 PetscReal *ar; 3476 3477 PetscFunctionBegin; 3478 mat = mxCreateDoubleMatrix(snes->conv_hist_len,1,mxREAL); 3479 ar = (PetscReal*) mxGetData(mat); 3480 for (i=0; i<snes->conv_hist_len; i++) ar[i] = snes->conv_hist[i]; 3481 PetscFunctionReturn(mat); 3482 } 3483 #endif 3484 3485 #undef __FUNCT__ 3486 #define __FUNCT__ "SNESGetConvergenceHistory" 3487 /*@C 3488 SNESGetConvergenceHistory - Gets the array used to hold the convergence history. 3489 3490 Not Collective 3491 3492 Input Parameter: 3493 . snes - iterative context obtained from SNESCreate() 3494 3495 Output Parameters: 3496 . a - array to hold history 3497 . its - integer array holds the number of linear iterations (or 3498 negative if not converged) for each solve. 3499 - na - size of a and its 3500 3501 Notes: 3502 The calling sequence for this routine in Fortran is 3503 $ call SNESGetConvergenceHistory(SNES snes, integer na, integer ierr) 3504 3505 This routine is useful, e.g., when running a code for purposes 3506 of accurate performance monitoring, when no I/O should be done 3507 during the section of code that is being timed. 3508 3509 Level: intermediate 3510 3511 .keywords: SNES, get, convergence, history 3512 3513 .seealso: SNESSetConvergencHistory() 3514 3515 @*/ 3516 PetscErrorCode SNESGetConvergenceHistory(SNES snes,PetscReal *a[],PetscInt *its[],PetscInt *na) 3517 { 3518 PetscFunctionBegin; 3519 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3520 if (a) *a = snes->conv_hist; 3521 if (its) *its = snes->conv_hist_its; 3522 if (na) *na = snes->conv_hist_len; 3523 PetscFunctionReturn(0); 3524 } 3525 3526 #undef __FUNCT__ 3527 #define __FUNCT__ "SNESSetUpdate" 3528 /*@C 3529 SNESSetUpdate - Sets the general-purpose update function called 3530 at the beginning of every iteration of the nonlinear solve. Specifically 3531 it is called just before the Jacobian is "evaluated". 3532 3533 Logically Collective on SNES 3534 3535 Input Parameters: 3536 . snes - The nonlinear solver context 3537 . func - The function 3538 3539 Calling sequence of func: 3540 . func (SNES snes, PetscInt step); 3541 3542 . step - The current step of the iteration 3543 3544 Level: advanced 3545 3546 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() 3547 This is not used by most users. 3548 3549 .keywords: SNES, update 3550 3551 .seealso SNESSetJacobian(), SNESSolve() 3552 @*/ 3553 PetscErrorCode SNESSetUpdate(SNES snes, PetscErrorCode (*func)(SNES, PetscInt)) 3554 { 3555 PetscFunctionBegin; 3556 PetscValidHeaderSpecific(snes, SNES_CLASSID,1); 3557 snes->ops->update = func; 3558 PetscFunctionReturn(0); 3559 } 3560 3561 #undef __FUNCT__ 3562 #define __FUNCT__ "SNESScaleStep_Private" 3563 /* 3564 SNESScaleStep_Private - Scales a step so that its length is less than the 3565 positive parameter delta. 3566 3567 Input Parameters: 3568 + snes - the SNES context 3569 . y - approximate solution of linear system 3570 . fnorm - 2-norm of current function 3571 - delta - trust region size 3572 3573 Output Parameters: 3574 + gpnorm - predicted function norm at the new point, assuming local 3575 linearization. The value is zero if the step lies within the trust 3576 region, and exceeds zero otherwise. 3577 - ynorm - 2-norm of the step 3578 3579 Note: 3580 For non-trust region methods such as SNESNEWTONLS, the parameter delta 3581 is set to be the maximum allowable step size. 3582 3583 .keywords: SNES, nonlinear, scale, step 3584 */ 3585 PetscErrorCode SNESScaleStep_Private(SNES snes,Vec y,PetscReal *fnorm,PetscReal *delta,PetscReal *gpnorm,PetscReal *ynorm) 3586 { 3587 PetscReal nrm; 3588 PetscScalar cnorm; 3589 PetscErrorCode ierr; 3590 3591 PetscFunctionBegin; 3592 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3593 PetscValidHeaderSpecific(y,VEC_CLASSID,2); 3594 PetscCheckSameComm(snes,1,y,2); 3595 3596 ierr = VecNorm(y,NORM_2,&nrm);CHKERRQ(ierr); 3597 if (nrm > *delta) { 3598 nrm = *delta/nrm; 3599 *gpnorm = (1.0 - nrm)*(*fnorm); 3600 cnorm = nrm; 3601 ierr = VecScale(y,cnorm);CHKERRQ(ierr); 3602 *ynorm = *delta; 3603 } else { 3604 *gpnorm = 0.0; 3605 *ynorm = nrm; 3606 } 3607 PetscFunctionReturn(0); 3608 } 3609 3610 #undef __FUNCT__ 3611 #define __FUNCT__ "SNESSolve" 3612 /*@C 3613 SNESSolve - Solves a nonlinear system F(x) = b. 3614 Call SNESSolve() after calling SNESCreate() and optional routines of the form SNESSetXXX(). 3615 3616 Collective on SNES 3617 3618 Input Parameters: 3619 + snes - the SNES context 3620 . b - the constant part of the equation F(x) = b, or NULL to use zero. 3621 - x - the solution vector. 3622 3623 Notes: 3624 The user should initialize the vector,x, with the initial guess 3625 for the nonlinear solve prior to calling SNESSolve. In particular, 3626 to employ an initial guess of zero, the user should explicitly set 3627 this vector to zero by calling VecSet(). 3628 3629 Level: beginner 3630 3631 .keywords: SNES, nonlinear, solve 3632 3633 .seealso: SNESCreate(), SNESDestroy(), SNESSetFunction(), SNESSetJacobian(), SNESSetGridSequence(), SNESGetSolution() 3634 @*/ 3635 PetscErrorCode SNESSolve(SNES snes,Vec b,Vec x) 3636 { 3637 PetscErrorCode ierr; 3638 PetscBool flg; 3639 PetscViewer viewer; 3640 PetscInt grid; 3641 Vec xcreated = NULL; 3642 DM dm; 3643 PetscViewerFormat format; 3644 3645 PetscFunctionBegin; 3646 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3647 if (x) PetscValidHeaderSpecific(x,VEC_CLASSID,3); 3648 if (x) PetscCheckSameComm(snes,1,x,3); 3649 if (b) PetscValidHeaderSpecific(b,VEC_CLASSID,2); 3650 if (b) PetscCheckSameComm(snes,1,b,2); 3651 3652 if (!x) { 3653 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 3654 ierr = DMCreateGlobalVector(dm,&xcreated);CHKERRQ(ierr); 3655 x = xcreated; 3656 } 3657 3658 ierr = PetscOptionsGetViewer(PetscObjectComm((PetscObject)snes),((PetscObject)snes)->prefix,"-snes_view_pre",&viewer,&format,&flg);CHKERRQ(ierr); 3659 if (flg && !PetscPreLoadingOn) { 3660 ierr = PetscViewerPushFormat(viewer,format);CHKERRQ(ierr); 3661 ierr = SNESView(snes,viewer);CHKERRQ(ierr); 3662 ierr = PetscViewerPopFormat(viewer);CHKERRQ(ierr); 3663 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 3664 } 3665 3666 for (grid=0; grid<snes->gridsequence; grid++) {ierr = PetscViewerASCIIPushTab(PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)snes)));CHKERRQ(ierr);} 3667 for (grid=0; grid<snes->gridsequence+1; grid++) { 3668 3669 /* set solution vector */ 3670 if (!grid) {ierr = PetscObjectReference((PetscObject)x);CHKERRQ(ierr);} 3671 ierr = VecDestroy(&snes->vec_sol);CHKERRQ(ierr); 3672 snes->vec_sol = x; 3673 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 3674 3675 /* set affine vector if provided */ 3676 if (b) { ierr = PetscObjectReference((PetscObject)b);CHKERRQ(ierr); } 3677 ierr = VecDestroy(&snes->vec_rhs);CHKERRQ(ierr); 3678 snes->vec_rhs = b; 3679 3680 ierr = SNESSetUp(snes);CHKERRQ(ierr); 3681 3682 if (!grid) { 3683 if (snes->ops->computeinitialguess) { 3684 ierr = (*snes->ops->computeinitialguess)(snes,snes->vec_sol,snes->initialguessP);CHKERRQ(ierr); 3685 } 3686 } 3687 3688 if (snes->conv_hist_reset) snes->conv_hist_len = 0; 3689 snes->nfuncs = 0; snes->linear_its = 0; snes->numFailures = 0; 3690 3691 ierr = PetscLogEventBegin(SNES_Solve,snes,0,0,0);CHKERRQ(ierr); 3692 ierr = (*snes->ops->solve)(snes);CHKERRQ(ierr); 3693 ierr = PetscLogEventEnd(SNES_Solve,snes,0,0,0);CHKERRQ(ierr); 3694 if (snes->domainerror) { 3695 snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; 3696 snes->domainerror = PETSC_FALSE; 3697 } 3698 if (!snes->reason) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Internal error, solver returned without setting converged reason"); 3699 3700 flg = PETSC_FALSE; 3701 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_test_local_min",&flg,NULL);CHKERRQ(ierr); 3702 if (flg && !PetscPreLoadingOn) { ierr = SNESTestLocalMin(snes);CHKERRQ(ierr); } 3703 if (snes->printreason) { 3704 ierr = PetscViewerASCIIAddTab(PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)snes)),((PetscObject)snes)->tablevel);CHKERRQ(ierr); 3705 if (snes->reason > 0) { 3706 ierr = PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)snes)),"Nonlinear solve converged due to %s iterations %D\n",SNESConvergedReasons[snes->reason],snes->iter);CHKERRQ(ierr); 3707 } else { 3708 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); 3709 } 3710 ierr = PetscViewerASCIISubtractTab(PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)snes)),((PetscObject)snes)->tablevel);CHKERRQ(ierr); 3711 } 3712 3713 if (snes->errorifnotconverged && snes->reason < 0) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_NOT_CONVERGED,"SNESSolve has not converged"); 3714 if (grid < snes->gridsequence) { 3715 DM fine; 3716 Vec xnew; 3717 Mat interp; 3718 3719 ierr = DMRefine(snes->dm,PetscObjectComm((PetscObject)snes),&fine);CHKERRQ(ierr); 3720 if (!fine) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_INCOMP,"DMRefine() did not perform any refinement, cannot continue grid sequencing"); 3721 ierr = DMCreateInterpolation(snes->dm,fine,&interp,NULL);CHKERRQ(ierr); 3722 ierr = DMCreateGlobalVector(fine,&xnew);CHKERRQ(ierr); 3723 ierr = MatInterpolate(interp,x,xnew);CHKERRQ(ierr); 3724 ierr = DMInterpolate(snes->dm,interp,fine);CHKERRQ(ierr); 3725 ierr = MatDestroy(&interp);CHKERRQ(ierr); 3726 x = xnew; 3727 3728 ierr = SNESReset(snes);CHKERRQ(ierr); 3729 ierr = SNESSetDM(snes,fine);CHKERRQ(ierr); 3730 ierr = DMDestroy(&fine);CHKERRQ(ierr); 3731 ierr = PetscViewerASCIIPopTab(PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)snes)));CHKERRQ(ierr); 3732 } 3733 } 3734 /* monitoring and viewing */ 3735 ierr = PetscOptionsGetViewer(PetscObjectComm((PetscObject)snes),((PetscObject)snes)->prefix,"-snes_view",&viewer,&format,&flg);CHKERRQ(ierr); 3736 if (flg && !PetscPreLoadingOn) { 3737 ierr = PetscViewerPushFormat(viewer,format);CHKERRQ(ierr); 3738 ierr = SNESView(snes,viewer);CHKERRQ(ierr); 3739 ierr = PetscViewerPopFormat(viewer);CHKERRQ(ierr); 3740 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 3741 } 3742 ierr = VecViewFromOptions(snes->vec_sol,((PetscObject)snes)->prefix,"-snes_view_solution");CHKERRQ(ierr); 3743 3744 ierr = VecDestroy(&xcreated);CHKERRQ(ierr); 3745 ierr = PetscObjectAMSBlock((PetscObject)snes);CHKERRQ(ierr); 3746 PetscFunctionReturn(0); 3747 } 3748 3749 /* --------- Internal routines for SNES Package --------- */ 3750 3751 #undef __FUNCT__ 3752 #define __FUNCT__ "SNESSetType" 3753 /*@C 3754 SNESSetType - Sets the method for the nonlinear solver. 3755 3756 Collective on SNES 3757 3758 Input Parameters: 3759 + snes - the SNES context 3760 - type - a known method 3761 3762 Options Database Key: 3763 . -snes_type <type> - Sets the method; use -help for a list 3764 of available methods (for instance, newtonls or newtontr) 3765 3766 Notes: 3767 See "petsc/include/petscsnes.h" for available methods (for instance) 3768 + SNESNEWTONLS - Newton's method with line search 3769 (systems of nonlinear equations) 3770 . SNESNEWTONTR - Newton's method with trust region 3771 (systems of nonlinear equations) 3772 3773 Normally, it is best to use the SNESSetFromOptions() command and then 3774 set the SNES solver type from the options database rather than by using 3775 this routine. Using the options database provides the user with 3776 maximum flexibility in evaluating the many nonlinear solvers. 3777 The SNESSetType() routine is provided for those situations where it 3778 is necessary to set the nonlinear solver independently of the command 3779 line or options database. This might be the case, for example, when 3780 the choice of solver changes during the execution of the program, 3781 and the user's application is taking responsibility for choosing the 3782 appropriate method. 3783 3784 Level: intermediate 3785 3786 .keywords: SNES, set, type 3787 3788 .seealso: SNESType, SNESCreate() 3789 3790 @*/ 3791 PetscErrorCode SNESSetType(SNES snes,SNESType type) 3792 { 3793 PetscErrorCode ierr,(*r)(SNES); 3794 PetscBool match; 3795 3796 PetscFunctionBegin; 3797 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3798 PetscValidCharPointer(type,2); 3799 3800 ierr = PetscObjectTypeCompare((PetscObject)snes,type,&match);CHKERRQ(ierr); 3801 if (match) PetscFunctionReturn(0); 3802 3803 ierr = PetscFunctionListFind(SNESList,type,&r);CHKERRQ(ierr); 3804 if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unable to find requested SNES type %s",type); 3805 /* Destroy the previous private SNES context */ 3806 if (snes->ops->destroy) { 3807 ierr = (*(snes)->ops->destroy)(snes);CHKERRQ(ierr); 3808 snes->ops->destroy = NULL; 3809 } 3810 /* Reinitialize function pointers in SNESOps structure */ 3811 snes->ops->setup = 0; 3812 snes->ops->solve = 0; 3813 snes->ops->view = 0; 3814 snes->ops->setfromoptions = 0; 3815 snes->ops->destroy = 0; 3816 /* Call the SNESCreate_XXX routine for this particular Nonlinear solver */ 3817 snes->setupcalled = PETSC_FALSE; 3818 3819 ierr = PetscObjectChangeTypeName((PetscObject)snes,type);CHKERRQ(ierr); 3820 ierr = (*r)(snes);CHKERRQ(ierr); 3821 PetscFunctionReturn(0); 3822 } 3823 3824 #undef __FUNCT__ 3825 #define __FUNCT__ "SNESGetType" 3826 /*@C 3827 SNESGetType - Gets the SNES method type and name (as a string). 3828 3829 Not Collective 3830 3831 Input Parameter: 3832 . snes - nonlinear solver context 3833 3834 Output Parameter: 3835 . type - SNES method (a character string) 3836 3837 Level: intermediate 3838 3839 .keywords: SNES, nonlinear, get, type, name 3840 @*/ 3841 PetscErrorCode SNESGetType(SNES snes,SNESType *type) 3842 { 3843 PetscFunctionBegin; 3844 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3845 PetscValidPointer(type,2); 3846 *type = ((PetscObject)snes)->type_name; 3847 PetscFunctionReturn(0); 3848 } 3849 3850 #undef __FUNCT__ 3851 #define __FUNCT__ "SNESGetSolution" 3852 /*@ 3853 SNESGetSolution - Returns the vector where the approximate solution is 3854 stored. This is the fine grid solution when using SNESSetGridSequence(). 3855 3856 Not Collective, but Vec is parallel if SNES is parallel 3857 3858 Input Parameter: 3859 . snes - the SNES context 3860 3861 Output Parameter: 3862 . x - the solution 3863 3864 Level: intermediate 3865 3866 .keywords: SNES, nonlinear, get, solution 3867 3868 .seealso: SNESGetSolutionUpdate(), SNESGetFunction() 3869 @*/ 3870 PetscErrorCode SNESGetSolution(SNES snes,Vec *x) 3871 { 3872 PetscFunctionBegin; 3873 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3874 PetscValidPointer(x,2); 3875 *x = snes->vec_sol; 3876 PetscFunctionReturn(0); 3877 } 3878 3879 #undef __FUNCT__ 3880 #define __FUNCT__ "SNESGetSolutionUpdate" 3881 /*@ 3882 SNESGetSolutionUpdate - Returns the vector where the solution update is 3883 stored. 3884 3885 Not Collective, but Vec is parallel if SNES is parallel 3886 3887 Input Parameter: 3888 . snes - the SNES context 3889 3890 Output Parameter: 3891 . x - the solution update 3892 3893 Level: advanced 3894 3895 .keywords: SNES, nonlinear, get, solution, update 3896 3897 .seealso: SNESGetSolution(), SNESGetFunction() 3898 @*/ 3899 PetscErrorCode SNESGetSolutionUpdate(SNES snes,Vec *x) 3900 { 3901 PetscFunctionBegin; 3902 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3903 PetscValidPointer(x,2); 3904 *x = snes->vec_sol_update; 3905 PetscFunctionReturn(0); 3906 } 3907 3908 #undef __FUNCT__ 3909 #define __FUNCT__ "SNESGetFunction" 3910 /*@C 3911 SNESGetFunction - Returns the vector where the function is stored. 3912 3913 Not Collective, but Vec is parallel if SNES is parallel. Collective if Vec is requested, but has not been created yet. 3914 3915 Input Parameter: 3916 . snes - the SNES context 3917 3918 Output Parameter: 3919 + r - the vector that is used to store residuals (or NULL if you don't want it) 3920 . SNESFunction- the function (or NULL if you don't want it) 3921 - ctx - the function context (or NULL if you don't want it) 3922 3923 Level: advanced 3924 3925 .keywords: SNES, nonlinear, get, function 3926 3927 .seealso: SNESSetFunction(), SNESGetSolution(), SNESFunction 3928 @*/ 3929 PetscErrorCode SNESGetFunction(SNES snes,Vec *r,PetscErrorCode (**SNESFunction)(SNES,Vec,Vec,void*),void **ctx) 3930 { 3931 PetscErrorCode ierr; 3932 DM dm; 3933 3934 PetscFunctionBegin; 3935 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3936 if (r) { 3937 if (!snes->vec_func) { 3938 if (snes->vec_rhs) { 3939 ierr = VecDuplicate(snes->vec_rhs,&snes->vec_func);CHKERRQ(ierr); 3940 } else if (snes->vec_sol) { 3941 ierr = VecDuplicate(snes->vec_sol,&snes->vec_func);CHKERRQ(ierr); 3942 } else if (snes->dm) { 3943 ierr = DMCreateGlobalVector(snes->dm,&snes->vec_func);CHKERRQ(ierr); 3944 } 3945 } 3946 *r = snes->vec_func; 3947 } 3948 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 3949 ierr = DMSNESGetFunction(dm,SNESFunction,ctx);CHKERRQ(ierr); 3950 PetscFunctionReturn(0); 3951 } 3952 3953 /*@C 3954 SNESGetGS - Returns the GS function and context. 3955 3956 Input Parameter: 3957 . snes - the SNES context 3958 3959 Output Parameter: 3960 + SNESGSFunction - the function (or NULL) 3961 - ctx - the function context (or NULL) 3962 3963 Level: advanced 3964 3965 .keywords: SNES, nonlinear, get, function 3966 3967 .seealso: SNESSetGS(), SNESGetFunction() 3968 @*/ 3969 3970 #undef __FUNCT__ 3971 #define __FUNCT__ "SNESGetGS" 3972 PetscErrorCode SNESGetGS (SNES snes, PetscErrorCode (**SNESGSFunction)(SNES, Vec, Vec, void*), void ** ctx) 3973 { 3974 PetscErrorCode ierr; 3975 DM dm; 3976 3977 PetscFunctionBegin; 3978 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3979 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 3980 ierr = DMSNESGetGS(dm,SNESGSFunction,ctx);CHKERRQ(ierr); 3981 PetscFunctionReturn(0); 3982 } 3983 3984 #undef __FUNCT__ 3985 #define __FUNCT__ "SNESSetOptionsPrefix" 3986 /*@C 3987 SNESSetOptionsPrefix - Sets the prefix used for searching for all 3988 SNES options in the database. 3989 3990 Logically Collective on SNES 3991 3992 Input Parameter: 3993 + snes - the SNES context 3994 - prefix - the prefix to prepend to all option names 3995 3996 Notes: 3997 A hyphen (-) must NOT be given at the beginning of the prefix name. 3998 The first character of all runtime options is AUTOMATICALLY the hyphen. 3999 4000 Level: advanced 4001 4002 .keywords: SNES, set, options, prefix, database 4003 4004 .seealso: SNESSetFromOptions() 4005 @*/ 4006 PetscErrorCode SNESSetOptionsPrefix(SNES snes,const char prefix[]) 4007 { 4008 PetscErrorCode ierr; 4009 4010 PetscFunctionBegin; 4011 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4012 ierr = PetscObjectSetOptionsPrefix((PetscObject)snes,prefix);CHKERRQ(ierr); 4013 if (!snes->ksp) {ierr = SNESGetKSP(snes,&snes->ksp);CHKERRQ(ierr);} 4014 if (snes->linesearch) { 4015 ierr = SNESGetLineSearch(snes,&snes->linesearch);CHKERRQ(ierr); 4016 ierr = PetscObjectSetOptionsPrefix((PetscObject)snes->linesearch,prefix);CHKERRQ(ierr); 4017 } 4018 ierr = KSPSetOptionsPrefix(snes->ksp,prefix);CHKERRQ(ierr); 4019 PetscFunctionReturn(0); 4020 } 4021 4022 #undef __FUNCT__ 4023 #define __FUNCT__ "SNESAppendOptionsPrefix" 4024 /*@C 4025 SNESAppendOptionsPrefix - Appends to the prefix used for searching for all 4026 SNES options in the database. 4027 4028 Logically Collective on SNES 4029 4030 Input Parameters: 4031 + snes - the SNES context 4032 - prefix - the prefix to prepend to all option names 4033 4034 Notes: 4035 A hyphen (-) must NOT be given at the beginning of the prefix name. 4036 The first character of all runtime options is AUTOMATICALLY the hyphen. 4037 4038 Level: advanced 4039 4040 .keywords: SNES, append, options, prefix, database 4041 4042 .seealso: SNESGetOptionsPrefix() 4043 @*/ 4044 PetscErrorCode SNESAppendOptionsPrefix(SNES snes,const char prefix[]) 4045 { 4046 PetscErrorCode ierr; 4047 4048 PetscFunctionBegin; 4049 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4050 ierr = PetscObjectAppendOptionsPrefix((PetscObject)snes,prefix);CHKERRQ(ierr); 4051 if (!snes->ksp) {ierr = SNESGetKSP(snes,&snes->ksp);CHKERRQ(ierr);} 4052 if (snes->linesearch) { 4053 ierr = SNESGetLineSearch(snes,&snes->linesearch);CHKERRQ(ierr); 4054 ierr = PetscObjectAppendOptionsPrefix((PetscObject)snes->linesearch,prefix);CHKERRQ(ierr); 4055 } 4056 ierr = KSPAppendOptionsPrefix(snes->ksp,prefix);CHKERRQ(ierr); 4057 PetscFunctionReturn(0); 4058 } 4059 4060 #undef __FUNCT__ 4061 #define __FUNCT__ "SNESGetOptionsPrefix" 4062 /*@C 4063 SNESGetOptionsPrefix - Sets the prefix used for searching for all 4064 SNES options in the database. 4065 4066 Not Collective 4067 4068 Input Parameter: 4069 . snes - the SNES context 4070 4071 Output Parameter: 4072 . prefix - pointer to the prefix string used 4073 4074 Notes: On the fortran side, the user should pass in a string 'prefix' of 4075 sufficient length to hold the prefix. 4076 4077 Level: advanced 4078 4079 .keywords: SNES, get, options, prefix, database 4080 4081 .seealso: SNESAppendOptionsPrefix() 4082 @*/ 4083 PetscErrorCode SNESGetOptionsPrefix(SNES snes,const char *prefix[]) 4084 { 4085 PetscErrorCode ierr; 4086 4087 PetscFunctionBegin; 4088 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4089 ierr = PetscObjectGetOptionsPrefix((PetscObject)snes,prefix);CHKERRQ(ierr); 4090 PetscFunctionReturn(0); 4091 } 4092 4093 4094 #undef __FUNCT__ 4095 #define __FUNCT__ "SNESRegister" 4096 /*@C 4097 SNESRegister - Adds a method to the nonlinear solver package. 4098 4099 Not collective 4100 4101 Input Parameters: 4102 + name_solver - name of a new user-defined solver 4103 - routine_create - routine to create method context 4104 4105 Notes: 4106 SNESRegister() may be called multiple times to add several user-defined solvers. 4107 4108 Sample usage: 4109 .vb 4110 SNESRegister("my_solver",MySolverCreate); 4111 .ve 4112 4113 Then, your solver can be chosen with the procedural interface via 4114 $ SNESSetType(snes,"my_solver") 4115 or at runtime via the option 4116 $ -snes_type my_solver 4117 4118 Level: advanced 4119 4120 Note: If your function is not being put into a shared library then use SNESRegister() instead 4121 4122 .keywords: SNES, nonlinear, register 4123 4124 .seealso: SNESRegisterAll(), SNESRegisterDestroy() 4125 4126 Level: advanced 4127 @*/ 4128 PetscErrorCode SNESRegister(const char sname[],PetscErrorCode (*function)(SNES)) 4129 { 4130 PetscErrorCode ierr; 4131 4132 PetscFunctionBegin; 4133 ierr = PetscFunctionListAdd(&SNESList,sname,function);CHKERRQ(ierr); 4134 PetscFunctionReturn(0); 4135 } 4136 4137 #undef __FUNCT__ 4138 #define __FUNCT__ "SNESTestLocalMin" 4139 PetscErrorCode SNESTestLocalMin(SNES snes) 4140 { 4141 PetscErrorCode ierr; 4142 PetscInt N,i,j; 4143 Vec u,uh,fh; 4144 PetscScalar value; 4145 PetscReal norm; 4146 4147 PetscFunctionBegin; 4148 ierr = SNESGetSolution(snes,&u);CHKERRQ(ierr); 4149 ierr = VecDuplicate(u,&uh);CHKERRQ(ierr); 4150 ierr = VecDuplicate(u,&fh);CHKERRQ(ierr); 4151 4152 /* currently only works for sequential */ 4153 ierr = PetscPrintf(PETSC_COMM_WORLD,"Testing FormFunction() for local min\n");CHKERRQ(ierr); 4154 ierr = VecGetSize(u,&N);CHKERRQ(ierr); 4155 for (i=0; i<N; i++) { 4156 ierr = VecCopy(u,uh);CHKERRQ(ierr); 4157 ierr = PetscPrintf(PETSC_COMM_WORLD,"i = %D\n",i);CHKERRQ(ierr); 4158 for (j=-10; j<11; j++) { 4159 value = PetscSign(j)*exp(PetscAbs(j)-10.0); 4160 ierr = VecSetValue(uh,i,value,ADD_VALUES);CHKERRQ(ierr); 4161 ierr = SNESComputeFunction(snes,uh,fh);CHKERRQ(ierr); 4162 ierr = VecNorm(fh,NORM_2,&norm);CHKERRQ(ierr); 4163 ierr = PetscPrintf(PETSC_COMM_WORLD," j norm %D %18.16e\n",j,norm);CHKERRQ(ierr); 4164 value = -value; 4165 ierr = VecSetValue(uh,i,value,ADD_VALUES);CHKERRQ(ierr); 4166 } 4167 } 4168 ierr = VecDestroy(&uh);CHKERRQ(ierr); 4169 ierr = VecDestroy(&fh);CHKERRQ(ierr); 4170 PetscFunctionReturn(0); 4171 } 4172 4173 #undef __FUNCT__ 4174 #define __FUNCT__ "SNESKSPSetUseEW" 4175 /*@ 4176 SNESKSPSetUseEW - Sets SNES use Eisenstat-Walker method for 4177 computing relative tolerance for linear solvers within an inexact 4178 Newton method. 4179 4180 Logically Collective on SNES 4181 4182 Input Parameters: 4183 + snes - SNES context 4184 - flag - PETSC_TRUE or PETSC_FALSE 4185 4186 Options Database: 4187 + -snes_ksp_ew - use Eisenstat-Walker method for determining linear system convergence 4188 . -snes_ksp_ew_version ver - version of Eisenstat-Walker method 4189 . -snes_ksp_ew_rtol0 <rtol0> - Sets rtol0 4190 . -snes_ksp_ew_rtolmax <rtolmax> - Sets rtolmax 4191 . -snes_ksp_ew_gamma <gamma> - Sets gamma 4192 . -snes_ksp_ew_alpha <alpha> - Sets alpha 4193 . -snes_ksp_ew_alpha2 <alpha2> - Sets alpha2 4194 - -snes_ksp_ew_threshold <threshold> - Sets threshold 4195 4196 Notes: 4197 Currently, the default is to use a constant relative tolerance for 4198 the inner linear solvers. Alternatively, one can use the 4199 Eisenstat-Walker method, where the relative convergence tolerance 4200 is reset at each Newton iteration according progress of the nonlinear 4201 solver. 4202 4203 Level: advanced 4204 4205 Reference: 4206 S. C. Eisenstat and H. F. Walker, "Choosing the forcing terms in an 4207 inexact Newton method", SISC 17 (1), pp.16-32, 1996. 4208 4209 .keywords: SNES, KSP, Eisenstat, Walker, convergence, test, inexact, Newton 4210 4211 .seealso: SNESKSPGetUseEW(), SNESKSPGetParametersEW(), SNESKSPSetParametersEW() 4212 @*/ 4213 PetscErrorCode SNESKSPSetUseEW(SNES snes,PetscBool flag) 4214 { 4215 PetscFunctionBegin; 4216 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4217 PetscValidLogicalCollectiveBool(snes,flag,2); 4218 snes->ksp_ewconv = flag; 4219 PetscFunctionReturn(0); 4220 } 4221 4222 #undef __FUNCT__ 4223 #define __FUNCT__ "SNESKSPGetUseEW" 4224 /*@ 4225 SNESKSPGetUseEW - Gets if SNES is using Eisenstat-Walker method 4226 for computing relative tolerance for linear solvers within an 4227 inexact Newton method. 4228 4229 Not Collective 4230 4231 Input Parameter: 4232 . snes - SNES context 4233 4234 Output Parameter: 4235 . flag - PETSC_TRUE or PETSC_FALSE 4236 4237 Notes: 4238 Currently, the default is to use a constant relative tolerance for 4239 the inner linear solvers. Alternatively, one can use the 4240 Eisenstat-Walker method, where the relative convergence tolerance 4241 is reset at each Newton iteration according progress of the nonlinear 4242 solver. 4243 4244 Level: advanced 4245 4246 Reference: 4247 S. C. Eisenstat and H. F. Walker, "Choosing the forcing terms in an 4248 inexact Newton method", SISC 17 (1), pp.16-32, 1996. 4249 4250 .keywords: SNES, KSP, Eisenstat, Walker, convergence, test, inexact, Newton 4251 4252 .seealso: SNESKSPSetUseEW(), SNESKSPGetParametersEW(), SNESKSPSetParametersEW() 4253 @*/ 4254 PetscErrorCode SNESKSPGetUseEW(SNES snes, PetscBool *flag) 4255 { 4256 PetscFunctionBegin; 4257 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4258 PetscValidPointer(flag,2); 4259 *flag = snes->ksp_ewconv; 4260 PetscFunctionReturn(0); 4261 } 4262 4263 #undef __FUNCT__ 4264 #define __FUNCT__ "SNESKSPSetParametersEW" 4265 /*@ 4266 SNESKSPSetParametersEW - Sets parameters for Eisenstat-Walker 4267 convergence criteria for the linear solvers within an inexact 4268 Newton method. 4269 4270 Logically Collective on SNES 4271 4272 Input Parameters: 4273 + snes - SNES context 4274 . version - version 1, 2 (default is 2) or 3 4275 . rtol_0 - initial relative tolerance (0 <= rtol_0 < 1) 4276 . rtol_max - maximum relative tolerance (0 <= rtol_max < 1) 4277 . gamma - multiplicative factor for version 2 rtol computation 4278 (0 <= gamma2 <= 1) 4279 . alpha - power for version 2 rtol computation (1 < alpha <= 2) 4280 . alpha2 - power for safeguard 4281 - threshold - threshold for imposing safeguard (0 < threshold < 1) 4282 4283 Note: 4284 Version 3 was contributed by Luis Chacon, June 2006. 4285 4286 Use PETSC_DEFAULT to retain the default for any of the parameters. 4287 4288 Level: advanced 4289 4290 Reference: 4291 S. C. Eisenstat and H. F. Walker, "Choosing the forcing terms in an 4292 inexact Newton method", Utah State University Math. Stat. Dept. Res. 4293 Report 6/94/75, June, 1994, to appear in SIAM J. Sci. Comput. 4294 4295 .keywords: SNES, KSP, Eisenstat, Walker, set, parameters 4296 4297 .seealso: SNESKSPSetUseEW(), SNESKSPGetUseEW(), SNESKSPGetParametersEW() 4298 @*/ 4299 PetscErrorCode SNESKSPSetParametersEW(SNES snes,PetscInt version,PetscReal rtol_0,PetscReal rtol_max,PetscReal gamma,PetscReal alpha,PetscReal alpha2,PetscReal threshold) 4300 { 4301 SNESKSPEW *kctx; 4302 4303 PetscFunctionBegin; 4304 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4305 kctx = (SNESKSPEW*)snes->kspconvctx; 4306 if (!kctx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"No Eisenstat-Walker context existing"); 4307 PetscValidLogicalCollectiveInt(snes,version,2); 4308 PetscValidLogicalCollectiveReal(snes,rtol_0,3); 4309 PetscValidLogicalCollectiveReal(snes,rtol_max,4); 4310 PetscValidLogicalCollectiveReal(snes,gamma,5); 4311 PetscValidLogicalCollectiveReal(snes,alpha,6); 4312 PetscValidLogicalCollectiveReal(snes,alpha2,7); 4313 PetscValidLogicalCollectiveReal(snes,threshold,8); 4314 4315 if (version != PETSC_DEFAULT) kctx->version = version; 4316 if (rtol_0 != PETSC_DEFAULT) kctx->rtol_0 = rtol_0; 4317 if (rtol_max != PETSC_DEFAULT) kctx->rtol_max = rtol_max; 4318 if (gamma != PETSC_DEFAULT) kctx->gamma = gamma; 4319 if (alpha != PETSC_DEFAULT) kctx->alpha = alpha; 4320 if (alpha2 != PETSC_DEFAULT) kctx->alpha2 = alpha2; 4321 if (threshold != PETSC_DEFAULT) kctx->threshold = threshold; 4322 4323 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); 4324 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); 4325 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); 4326 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); 4327 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); 4328 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); 4329 PetscFunctionReturn(0); 4330 } 4331 4332 #undef __FUNCT__ 4333 #define __FUNCT__ "SNESKSPGetParametersEW" 4334 /*@ 4335 SNESKSPGetParametersEW - Gets parameters for Eisenstat-Walker 4336 convergence criteria for the linear solvers within an inexact 4337 Newton method. 4338 4339 Not Collective 4340 4341 Input Parameters: 4342 snes - SNES context 4343 4344 Output Parameters: 4345 + version - version 1, 2 (default is 2) or 3 4346 . rtol_0 - initial relative tolerance (0 <= rtol_0 < 1) 4347 . rtol_max - maximum relative tolerance (0 <= rtol_max < 1) 4348 . gamma - multiplicative factor for version 2 rtol computation (0 <= gamma2 <= 1) 4349 . alpha - power for version 2 rtol computation (1 < alpha <= 2) 4350 . alpha2 - power for safeguard 4351 - threshold - threshold for imposing safeguard (0 < threshold < 1) 4352 4353 Level: advanced 4354 4355 .keywords: SNES, KSP, Eisenstat, Walker, get, parameters 4356 4357 .seealso: SNESKSPSetUseEW(), SNESKSPGetUseEW(), SNESKSPSetParametersEW() 4358 @*/ 4359 PetscErrorCode SNESKSPGetParametersEW(SNES snes,PetscInt *version,PetscReal *rtol_0,PetscReal *rtol_max,PetscReal *gamma,PetscReal *alpha,PetscReal *alpha2,PetscReal *threshold) 4360 { 4361 SNESKSPEW *kctx; 4362 4363 PetscFunctionBegin; 4364 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4365 kctx = (SNESKSPEW*)snes->kspconvctx; 4366 if (!kctx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"No Eisenstat-Walker context existing"); 4367 if (version) *version = kctx->version; 4368 if (rtol_0) *rtol_0 = kctx->rtol_0; 4369 if (rtol_max) *rtol_max = kctx->rtol_max; 4370 if (gamma) *gamma = kctx->gamma; 4371 if (alpha) *alpha = kctx->alpha; 4372 if (alpha2) *alpha2 = kctx->alpha2; 4373 if (threshold) *threshold = kctx->threshold; 4374 PetscFunctionReturn(0); 4375 } 4376 4377 #undef __FUNCT__ 4378 #define __FUNCT__ "SNESKSPEW_PreSolve" 4379 PetscErrorCode SNESKSPEW_PreSolve(KSP ksp, Vec b, Vec x, SNES snes) 4380 { 4381 PetscErrorCode ierr; 4382 SNESKSPEW *kctx = (SNESKSPEW*)snes->kspconvctx; 4383 PetscReal rtol = PETSC_DEFAULT,stol; 4384 4385 PetscFunctionBegin; 4386 if (!snes->ksp_ewconv) PetscFunctionReturn(0); 4387 if (!snes->iter) rtol = kctx->rtol_0; /* first time in, so use the original user rtol */ 4388 else { 4389 if (kctx->version == 1) { 4390 rtol = (snes->norm - kctx->lresid_last)/kctx->norm_last; 4391 if (rtol < 0.0) rtol = -rtol; 4392 stol = PetscPowReal(kctx->rtol_last,kctx->alpha2); 4393 if (stol > kctx->threshold) rtol = PetscMax(rtol,stol); 4394 } else if (kctx->version == 2) { 4395 rtol = kctx->gamma * PetscPowReal(snes->norm/kctx->norm_last,kctx->alpha); 4396 stol = kctx->gamma * PetscPowReal(kctx->rtol_last,kctx->alpha); 4397 if (stol > kctx->threshold) rtol = PetscMax(rtol,stol); 4398 } else if (kctx->version == 3) { /* contributed by Luis Chacon, June 2006. */ 4399 rtol = kctx->gamma * PetscPowReal(snes->norm/kctx->norm_last,kctx->alpha); 4400 /* safeguard: avoid sharp decrease of rtol */ 4401 stol = kctx->gamma*PetscPowReal(kctx->rtol_last,kctx->alpha); 4402 stol = PetscMax(rtol,stol); 4403 rtol = PetscMin(kctx->rtol_0,stol); 4404 /* safeguard: avoid oversolving */ 4405 stol = kctx->gamma*(snes->ttol)/snes->norm; 4406 stol = PetscMax(rtol,stol); 4407 rtol = PetscMin(kctx->rtol_0,stol); 4408 } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Only versions 1, 2 or 3 are supported: %D",kctx->version); 4409 } 4410 /* safeguard: avoid rtol greater than one */ 4411 rtol = PetscMin(rtol,kctx->rtol_max); 4412 ierr = KSPSetTolerances(ksp,rtol,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);CHKERRQ(ierr); 4413 ierr = PetscInfo3(snes,"iter %D, Eisenstat-Walker (version %D) KSP rtol=%G\n",snes->iter,kctx->version,rtol);CHKERRQ(ierr); 4414 PetscFunctionReturn(0); 4415 } 4416 4417 #undef __FUNCT__ 4418 #define __FUNCT__ "SNESKSPEW_PostSolve" 4419 PetscErrorCode SNESKSPEW_PostSolve(KSP ksp, Vec b, Vec x, SNES snes) 4420 { 4421 PetscErrorCode ierr; 4422 SNESKSPEW *kctx = (SNESKSPEW*)snes->kspconvctx; 4423 PCSide pcside; 4424 Vec lres; 4425 4426 PetscFunctionBegin; 4427 if (!snes->ksp_ewconv) PetscFunctionReturn(0); 4428 ierr = KSPGetTolerances(ksp,&kctx->rtol_last,0,0,0);CHKERRQ(ierr); 4429 ierr = SNESGetFunctionNorm(snes,&kctx->norm_last);CHKERRQ(ierr); 4430 if (kctx->version == 1) { 4431 ierr = KSPGetPCSide(ksp,&pcside);CHKERRQ(ierr); 4432 if (pcside == PC_RIGHT) { /* XXX Should we also test KSP_UNPRECONDITIONED_NORM ? */ 4433 /* KSP residual is true linear residual */ 4434 ierr = KSPGetResidualNorm(ksp,&kctx->lresid_last);CHKERRQ(ierr); 4435 } else { 4436 /* KSP residual is preconditioned residual */ 4437 /* compute true linear residual norm */ 4438 ierr = VecDuplicate(b,&lres);CHKERRQ(ierr); 4439 ierr = MatMult(snes->jacobian,x,lres);CHKERRQ(ierr); 4440 ierr = VecAYPX(lres,-1.0,b);CHKERRQ(ierr); 4441 ierr = VecNorm(lres,NORM_2,&kctx->lresid_last);CHKERRQ(ierr); 4442 ierr = VecDestroy(&lres);CHKERRQ(ierr); 4443 } 4444 } 4445 PetscFunctionReturn(0); 4446 } 4447 4448 #undef __FUNCT__ 4449 #define __FUNCT__ "SNESGetKSP" 4450 /*@ 4451 SNESGetKSP - Returns the KSP context for a SNES solver. 4452 4453 Not Collective, but if SNES object is parallel, then KSP object is parallel 4454 4455 Input Parameter: 4456 . snes - the SNES context 4457 4458 Output Parameter: 4459 . ksp - the KSP context 4460 4461 Notes: 4462 The user can then directly manipulate the KSP context to set various 4463 options, etc. Likewise, the user can then extract and manipulate the 4464 PC contexts as well. 4465 4466 Level: beginner 4467 4468 .keywords: SNES, nonlinear, get, KSP, context 4469 4470 .seealso: KSPGetPC(), SNESCreate(), KSPCreate(), SNESSetKSP() 4471 @*/ 4472 PetscErrorCode SNESGetKSP(SNES snes,KSP *ksp) 4473 { 4474 PetscErrorCode ierr; 4475 4476 PetscFunctionBegin; 4477 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4478 PetscValidPointer(ksp,2); 4479 4480 if (!snes->ksp) { 4481 ierr = KSPCreate(PetscObjectComm((PetscObject)snes),&snes->ksp);CHKERRQ(ierr); 4482 ierr = PetscObjectIncrementTabLevel((PetscObject)snes->ksp,(PetscObject)snes,1);CHKERRQ(ierr); 4483 ierr = PetscLogObjectParent(snes,snes->ksp);CHKERRQ(ierr); 4484 4485 ierr = KSPSetPreSolve(snes->ksp,(PetscErrorCode (*)(KSP,Vec,Vec,void*))SNESKSPEW_PreSolve,snes);CHKERRQ(ierr); 4486 ierr = KSPSetPostSolve(snes->ksp,(PetscErrorCode (*)(KSP,Vec,Vec,void*))SNESKSPEW_PostSolve,snes);CHKERRQ(ierr); 4487 } 4488 *ksp = snes->ksp; 4489 PetscFunctionReturn(0); 4490 } 4491 4492 4493 #include <petsc-private/dmimpl.h> 4494 #undef __FUNCT__ 4495 #define __FUNCT__ "SNESSetDM" 4496 /*@ 4497 SNESSetDM - Sets the DM that may be used by some preconditioners 4498 4499 Logically Collective on SNES 4500 4501 Input Parameters: 4502 + snes - the preconditioner context 4503 - dm - the dm 4504 4505 Level: intermediate 4506 4507 .seealso: SNESGetDM(), KSPSetDM(), KSPGetDM() 4508 @*/ 4509 PetscErrorCode SNESSetDM(SNES snes,DM dm) 4510 { 4511 PetscErrorCode ierr; 4512 KSP ksp; 4513 DMSNES sdm; 4514 4515 PetscFunctionBegin; 4516 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4517 if (dm) {ierr = PetscObjectReference((PetscObject)dm);CHKERRQ(ierr);} 4518 if (snes->dm) { /* Move the DMSNES context over to the new DM unless the new DM already has one */ 4519 if (snes->dm->dmsnes && snes->dmAuto && !dm->dmsnes) { 4520 ierr = DMCopyDMSNES(snes->dm,dm);CHKERRQ(ierr); 4521 ierr = DMGetDMSNES(snes->dm,&sdm);CHKERRQ(ierr); 4522 if (sdm->originaldm == snes->dm) sdm->originaldm = dm; /* Grant write privileges to the replacement DM */ 4523 } 4524 ierr = DMDestroy(&snes->dm);CHKERRQ(ierr); 4525 } 4526 snes->dm = dm; 4527 snes->dmAuto = PETSC_FALSE; 4528 4529 ierr = SNESGetKSP(snes,&ksp);CHKERRQ(ierr); 4530 ierr = KSPSetDM(ksp,dm);CHKERRQ(ierr); 4531 ierr = KSPSetDMActive(ksp,PETSC_FALSE);CHKERRQ(ierr); 4532 if (snes->pc) { 4533 ierr = SNESSetDM(snes->pc, snes->dm);CHKERRQ(ierr); 4534 ierr = SNESSetPCSide(snes,snes->pcside);CHKERRQ(ierr); 4535 } 4536 PetscFunctionReturn(0); 4537 } 4538 4539 #undef __FUNCT__ 4540 #define __FUNCT__ "SNESGetDM" 4541 /*@ 4542 SNESGetDM - Gets the DM that may be used by some preconditioners 4543 4544 Not Collective but DM obtained is parallel on SNES 4545 4546 Input Parameter: 4547 . snes - the preconditioner context 4548 4549 Output Parameter: 4550 . dm - the dm 4551 4552 Level: intermediate 4553 4554 .seealso: SNESSetDM(), KSPSetDM(), KSPGetDM() 4555 @*/ 4556 PetscErrorCode SNESGetDM(SNES snes,DM *dm) 4557 { 4558 PetscErrorCode ierr; 4559 4560 PetscFunctionBegin; 4561 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4562 if (!snes->dm) { 4563 ierr = DMShellCreate(PetscObjectComm((PetscObject)snes),&snes->dm);CHKERRQ(ierr); 4564 snes->dmAuto = PETSC_TRUE; 4565 } 4566 *dm = snes->dm; 4567 PetscFunctionReturn(0); 4568 } 4569 4570 #undef __FUNCT__ 4571 #define __FUNCT__ "SNESSetPC" 4572 /*@ 4573 SNESSetPC - Sets the nonlinear preconditioner to be used. 4574 4575 Collective on SNES 4576 4577 Input Parameters: 4578 + snes - iterative context obtained from SNESCreate() 4579 - pc - the preconditioner object 4580 4581 Notes: 4582 Use SNESGetPC() to retrieve the preconditioner context (for example, 4583 to configure it using the API). 4584 4585 Level: developer 4586 4587 .keywords: SNES, set, precondition 4588 .seealso: SNESGetPC() 4589 @*/ 4590 PetscErrorCode SNESSetPC(SNES snes, SNES pc) 4591 { 4592 PetscErrorCode ierr; 4593 4594 PetscFunctionBegin; 4595 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); 4596 PetscValidHeaderSpecific(pc, SNES_CLASSID, 2); 4597 PetscCheckSameComm(snes, 1, pc, 2); 4598 ierr = PetscObjectReference((PetscObject) pc);CHKERRQ(ierr); 4599 ierr = SNESDestroy(&snes->pc);CHKERRQ(ierr); 4600 snes->pc = pc; 4601 ierr = PetscLogObjectParent(snes, snes->pc);CHKERRQ(ierr); 4602 PetscFunctionReturn(0); 4603 } 4604 4605 #undef __FUNCT__ 4606 #define __FUNCT__ "SNESGetPC" 4607 /*@ 4608 SNESGetPC - Returns a pointer to the nonlinear preconditioning context set with SNESSetPC(). 4609 4610 Not Collective 4611 4612 Input Parameter: 4613 . snes - iterative context obtained from SNESCreate() 4614 4615 Output Parameter: 4616 . pc - preconditioner context 4617 4618 Level: developer 4619 4620 .keywords: SNES, get, preconditioner 4621 .seealso: SNESSetPC() 4622 @*/ 4623 PetscErrorCode SNESGetPC(SNES snes, SNES *pc) 4624 { 4625 PetscErrorCode ierr; 4626 const char *optionsprefix; 4627 4628 PetscFunctionBegin; 4629 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); 4630 PetscValidPointer(pc, 2); 4631 if (!snes->pc) { 4632 ierr = SNESCreate(PetscObjectComm((PetscObject)snes),&snes->pc);CHKERRQ(ierr); 4633 ierr = PetscObjectIncrementTabLevel((PetscObject)snes->pc,(PetscObject)snes,1);CHKERRQ(ierr); 4634 ierr = PetscLogObjectParent(snes,snes->pc);CHKERRQ(ierr); 4635 ierr = SNESGetOptionsPrefix(snes,&optionsprefix);CHKERRQ(ierr); 4636 ierr = SNESSetOptionsPrefix(snes->pc,optionsprefix);CHKERRQ(ierr); 4637 ierr = SNESAppendOptionsPrefix(snes->pc,"npc_");CHKERRQ(ierr); 4638 } 4639 *pc = snes->pc; 4640 PetscFunctionReturn(0); 4641 } 4642 4643 #undef __FUNCT__ 4644 #define __FUNCT__ "SNESSetPCSide" 4645 /*@ 4646 SNESSetPCSide - Sets the preconditioning side. 4647 4648 Logically Collective on SNES 4649 4650 Input Parameter: 4651 . snes - iterative context obtained from SNESCreate() 4652 4653 Output Parameter: 4654 . side - the preconditioning side, where side is one of 4655 .vb 4656 PC_LEFT - left preconditioning (default) 4657 PC_RIGHT - right preconditioning 4658 .ve 4659 4660 Options Database Keys: 4661 . -snes_pc_side <right,left> 4662 4663 Level: intermediate 4664 4665 .keywords: SNES, set, right, left, side, preconditioner, flag 4666 4667 .seealso: SNESGetPCSide(), KSPSetPCSide() 4668 @*/ 4669 PetscErrorCode SNESSetPCSide(SNES snes,PCSide side) 4670 { 4671 PetscFunctionBegin; 4672 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4673 PetscValidLogicalCollectiveEnum(snes,side,2); 4674 snes->pcside = side; 4675 PetscFunctionReturn(0); 4676 } 4677 4678 #undef __FUNCT__ 4679 #define __FUNCT__ "SNESGetPCSide" 4680 /*@ 4681 SNESGetPCSide - Gets the preconditioning side. 4682 4683 Not Collective 4684 4685 Input Parameter: 4686 . snes - iterative context obtained from SNESCreate() 4687 4688 Output Parameter: 4689 . side - the preconditioning side, where side is one of 4690 .vb 4691 PC_LEFT - left preconditioning (default) 4692 PC_RIGHT - right preconditioning 4693 .ve 4694 4695 Level: intermediate 4696 4697 .keywords: SNES, get, right, left, side, preconditioner, flag 4698 4699 .seealso: SNESSetPCSide(), KSPGetPCSide() 4700 @*/ 4701 PetscErrorCode SNESGetPCSide(SNES snes,PCSide *side) 4702 { 4703 PetscFunctionBegin; 4704 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4705 PetscValidPointer(side,2); 4706 *side = snes->pcside; 4707 PetscFunctionReturn(0); 4708 } 4709 4710 #undef __FUNCT__ 4711 #define __FUNCT__ "SNESSetLineSearch" 4712 /*@ 4713 SNESSetLineSearch - Sets the linesearch on the SNES instance. 4714 4715 Collective on SNES 4716 4717 Input Parameters: 4718 + snes - iterative context obtained from SNESCreate() 4719 - linesearch - the linesearch object 4720 4721 Notes: 4722 Use SNESGetLineSearch() to retrieve the preconditioner context (for example, 4723 to configure it using the API). 4724 4725 Level: developer 4726 4727 .keywords: SNES, set, linesearch 4728 .seealso: SNESGetLineSearch() 4729 @*/ 4730 PetscErrorCode SNESSetLineSearch(SNES snes, SNESLineSearch linesearch) 4731 { 4732 PetscErrorCode ierr; 4733 4734 PetscFunctionBegin; 4735 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); 4736 PetscValidHeaderSpecific(linesearch, SNESLINESEARCH_CLASSID, 2); 4737 PetscCheckSameComm(snes, 1, linesearch, 2); 4738 ierr = PetscObjectReference((PetscObject) linesearch);CHKERRQ(ierr); 4739 ierr = SNESLineSearchDestroy(&snes->linesearch);CHKERRQ(ierr); 4740 4741 snes->linesearch = linesearch; 4742 4743 ierr = PetscLogObjectParent(snes, snes->linesearch);CHKERRQ(ierr); 4744 PetscFunctionReturn(0); 4745 } 4746 4747 #undef __FUNCT__ 4748 #define __FUNCT__ "SNESGetLineSearch" 4749 /*@ 4750 SNESGetLineSearch - Returns a pointer to the line search context set with SNESSetLineSearch() 4751 or creates a default line search instance associated with the SNES and returns it. 4752 4753 Not Collective 4754 4755 Input Parameter: 4756 . snes - iterative context obtained from SNESCreate() 4757 4758 Output Parameter: 4759 . linesearch - linesearch context 4760 4761 Level: developer 4762 4763 .keywords: SNES, get, linesearch 4764 .seealso: SNESSetLineSearch() 4765 @*/ 4766 PetscErrorCode SNESGetLineSearch(SNES snes, SNESLineSearch *linesearch) 4767 { 4768 PetscErrorCode ierr; 4769 const char *optionsprefix; 4770 4771 PetscFunctionBegin; 4772 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); 4773 PetscValidPointer(linesearch, 2); 4774 if (!snes->linesearch) { 4775 ierr = SNESGetOptionsPrefix(snes, &optionsprefix);CHKERRQ(ierr); 4776 ierr = SNESLineSearchCreate(PetscObjectComm((PetscObject)snes), &snes->linesearch);CHKERRQ(ierr); 4777 ierr = SNESLineSearchSetSNES(snes->linesearch, snes);CHKERRQ(ierr); 4778 ierr = SNESLineSearchAppendOptionsPrefix(snes->linesearch, optionsprefix);CHKERRQ(ierr); 4779 ierr = PetscObjectIncrementTabLevel((PetscObject) snes->linesearch, (PetscObject) snes, 1);CHKERRQ(ierr); 4780 ierr = PetscLogObjectParent(snes, snes->linesearch);CHKERRQ(ierr); 4781 } 4782 *linesearch = snes->linesearch; 4783 PetscFunctionReturn(0); 4784 } 4785 4786 #if defined(PETSC_HAVE_MATLAB_ENGINE) 4787 #include <mex.h> 4788 4789 typedef struct {char *funcname; mxArray *ctx;} SNESMatlabContext; 4790 4791 #undef __FUNCT__ 4792 #define __FUNCT__ "SNESComputeFunction_Matlab" 4793 /* 4794 SNESComputeFunction_Matlab - Calls the function that has been set with SNESSetFunctionMatlab(). 4795 4796 Collective on SNES 4797 4798 Input Parameters: 4799 + snes - the SNES context 4800 - x - input vector 4801 4802 Output Parameter: 4803 . y - function vector, as set by SNESSetFunction() 4804 4805 Notes: 4806 SNESComputeFunction() is typically used within nonlinear solvers 4807 implementations, so most users would not generally call this routine 4808 themselves. 4809 4810 Level: developer 4811 4812 .keywords: SNES, nonlinear, compute, function 4813 4814 .seealso: SNESSetFunction(), SNESGetFunction() 4815 */ 4816 PetscErrorCode SNESComputeFunction_Matlab(SNES snes,Vec x,Vec y, void *ctx) 4817 { 4818 PetscErrorCode ierr; 4819 SNESMatlabContext *sctx = (SNESMatlabContext*)ctx; 4820 int nlhs = 1,nrhs = 5; 4821 mxArray *plhs[1],*prhs[5]; 4822 long long int lx = 0,ly = 0,ls = 0; 4823 4824 PetscFunctionBegin; 4825 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4826 PetscValidHeaderSpecific(x,VEC_CLASSID,2); 4827 PetscValidHeaderSpecific(y,VEC_CLASSID,3); 4828 PetscCheckSameComm(snes,1,x,2); 4829 PetscCheckSameComm(snes,1,y,3); 4830 4831 /* call Matlab function in ctx with arguments x and y */ 4832 4833 ierr = PetscMemcpy(&ls,&snes,sizeof(snes));CHKERRQ(ierr); 4834 ierr = PetscMemcpy(&lx,&x,sizeof(x));CHKERRQ(ierr); 4835 ierr = PetscMemcpy(&ly,&y,sizeof(x));CHKERRQ(ierr); 4836 prhs[0] = mxCreateDoubleScalar((double)ls); 4837 prhs[1] = mxCreateDoubleScalar((double)lx); 4838 prhs[2] = mxCreateDoubleScalar((double)ly); 4839 prhs[3] = mxCreateString(sctx->funcname); 4840 prhs[4] = sctx->ctx; 4841 ierr = mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscSNESComputeFunctionInternal");CHKERRQ(ierr); 4842 ierr = mxGetScalar(plhs[0]);CHKERRQ(ierr); 4843 mxDestroyArray(prhs[0]); 4844 mxDestroyArray(prhs[1]); 4845 mxDestroyArray(prhs[2]); 4846 mxDestroyArray(prhs[3]); 4847 mxDestroyArray(plhs[0]); 4848 PetscFunctionReturn(0); 4849 } 4850 4851 #undef __FUNCT__ 4852 #define __FUNCT__ "SNESSetFunctionMatlab" 4853 /* 4854 SNESSetFunctionMatlab - Sets the function evaluation routine and function 4855 vector for use by the SNES routines in solving systems of nonlinear 4856 equations from MATLAB. Here the function is a string containing the name of a MATLAB function 4857 4858 Logically Collective on SNES 4859 4860 Input Parameters: 4861 + snes - the SNES context 4862 . r - vector to store function value 4863 - SNESFunction - function evaluation routine 4864 4865 Notes: 4866 The Newton-like methods typically solve linear systems of the form 4867 $ f'(x) x = -f(x), 4868 where f'(x) denotes the Jacobian matrix and f(x) is the function. 4869 4870 Level: beginner 4871 4872 Developer Note: This bleeds the allocated memory SNESMatlabContext *sctx; 4873 4874 .keywords: SNES, nonlinear, set, function 4875 4876 .seealso: SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESSetFunction() 4877 */ 4878 PetscErrorCode SNESSetFunctionMatlab(SNES snes,Vec r,const char *SNESFunction,mxArray *ctx) 4879 { 4880 PetscErrorCode ierr; 4881 SNESMatlabContext *sctx; 4882 4883 PetscFunctionBegin; 4884 /* currently sctx is memory bleed */ 4885 ierr = PetscMalloc(sizeof(SNESMatlabContext),&sctx);CHKERRQ(ierr); 4886 ierr = PetscStrallocpy(SNESFunction,&sctx->funcname);CHKERRQ(ierr); 4887 /* 4888 This should work, but it doesn't 4889 sctx->ctx = ctx; 4890 mexMakeArrayPersistent(sctx->ctx); 4891 */ 4892 sctx->ctx = mxDuplicateArray(ctx); 4893 ierr = SNESSetFunction(snes,r,SNESComputeFunction_Matlab,sctx);CHKERRQ(ierr); 4894 PetscFunctionReturn(0); 4895 } 4896 4897 #undef __FUNCT__ 4898 #define __FUNCT__ "SNESComputeJacobian_Matlab" 4899 /* 4900 SNESComputeJacobian_Matlab - Calls the function that has been set with SNESSetJacobianMatlab(). 4901 4902 Collective on SNES 4903 4904 Input Parameters: 4905 + snes - the SNES context 4906 . x - input vector 4907 . A, B - the matrices 4908 - ctx - user context 4909 4910 Output Parameter: 4911 . flag - structure of the matrix 4912 4913 Level: developer 4914 4915 .keywords: SNES, nonlinear, compute, function 4916 4917 .seealso: SNESSetFunction(), SNESGetFunction() 4918 @*/ 4919 PetscErrorCode SNESComputeJacobian_Matlab(SNES snes,Vec x,Mat *A,Mat *B,MatStructure *flag, void *ctx) 4920 { 4921 PetscErrorCode ierr; 4922 SNESMatlabContext *sctx = (SNESMatlabContext*)ctx; 4923 int nlhs = 2,nrhs = 6; 4924 mxArray *plhs[2],*prhs[6]; 4925 long long int lx = 0,lA = 0,ls = 0, lB = 0; 4926 4927 PetscFunctionBegin; 4928 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4929 PetscValidHeaderSpecific(x,VEC_CLASSID,2); 4930 4931 /* call Matlab function in ctx with arguments x and y */ 4932 4933 ierr = PetscMemcpy(&ls,&snes,sizeof(snes));CHKERRQ(ierr); 4934 ierr = PetscMemcpy(&lx,&x,sizeof(x));CHKERRQ(ierr); 4935 ierr = PetscMemcpy(&lA,A,sizeof(x));CHKERRQ(ierr); 4936 ierr = PetscMemcpy(&lB,B,sizeof(x));CHKERRQ(ierr); 4937 prhs[0] = mxCreateDoubleScalar((double)ls); 4938 prhs[1] = mxCreateDoubleScalar((double)lx); 4939 prhs[2] = mxCreateDoubleScalar((double)lA); 4940 prhs[3] = mxCreateDoubleScalar((double)lB); 4941 prhs[4] = mxCreateString(sctx->funcname); 4942 prhs[5] = sctx->ctx; 4943 ierr = mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscSNESComputeJacobianInternal");CHKERRQ(ierr); 4944 ierr = mxGetScalar(plhs[0]);CHKERRQ(ierr); 4945 *flag = (MatStructure) mxGetScalar(plhs[1]);CHKERRQ(ierr); 4946 mxDestroyArray(prhs[0]); 4947 mxDestroyArray(prhs[1]); 4948 mxDestroyArray(prhs[2]); 4949 mxDestroyArray(prhs[3]); 4950 mxDestroyArray(prhs[4]); 4951 mxDestroyArray(plhs[0]); 4952 mxDestroyArray(plhs[1]); 4953 PetscFunctionReturn(0); 4954 } 4955 4956 #undef __FUNCT__ 4957 #define __FUNCT__ "SNESSetJacobianMatlab" 4958 /* 4959 SNESSetJacobianMatlab - Sets the Jacobian function evaluation routine and two empty Jacobian matrices 4960 vector for use by the SNES routines in solving systems of nonlinear 4961 equations from MATLAB. Here the function is a string containing the name of a MATLAB function 4962 4963 Logically Collective on SNES 4964 4965 Input Parameters: 4966 + snes - the SNES context 4967 . A,B - Jacobian matrices 4968 . SNESJacobianFunction - function evaluation routine 4969 - ctx - user context 4970 4971 Level: developer 4972 4973 Developer Note: This bleeds the allocated memory SNESMatlabContext *sctx; 4974 4975 .keywords: SNES, nonlinear, set, function 4976 4977 .seealso: SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESSetFunction(), SNESJacobianFunction 4978 */ 4979 PetscErrorCode SNESSetJacobianMatlab(SNES snes,Mat A,Mat B,const char *SNESJacobianFunction,mxArray *ctx) 4980 { 4981 PetscErrorCode ierr; 4982 SNESMatlabContext *sctx; 4983 4984 PetscFunctionBegin; 4985 /* currently sctx is memory bleed */ 4986 ierr = PetscMalloc(sizeof(SNESMatlabContext),&sctx);CHKERRQ(ierr); 4987 ierr = PetscStrallocpy(SNESJacobianFunction,&sctx->funcname);CHKERRQ(ierr); 4988 /* 4989 This should work, but it doesn't 4990 sctx->ctx = ctx; 4991 mexMakeArrayPersistent(sctx->ctx); 4992 */ 4993 sctx->ctx = mxDuplicateArray(ctx); 4994 ierr = SNESSetJacobian(snes,A,B,SNESComputeJacobian_Matlab,sctx);CHKERRQ(ierr); 4995 PetscFunctionReturn(0); 4996 } 4997 4998 #undef __FUNCT__ 4999 #define __FUNCT__ "SNESMonitor_Matlab" 5000 /* 5001 SNESMonitor_Matlab - Calls the function that has been set with SNESMonitorSetMatlab(). 5002 5003 Collective on SNES 5004 5005 .seealso: SNESSetFunction(), SNESGetFunction() 5006 @*/ 5007 PetscErrorCode SNESMonitor_Matlab(SNES snes,PetscInt it, PetscReal fnorm, void *ctx) 5008 { 5009 PetscErrorCode ierr; 5010 SNESMatlabContext *sctx = (SNESMatlabContext*)ctx; 5011 int nlhs = 1,nrhs = 6; 5012 mxArray *plhs[1],*prhs[6]; 5013 long long int lx = 0,ls = 0; 5014 Vec x = snes->vec_sol; 5015 5016 PetscFunctionBegin; 5017 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 5018 5019 ierr = PetscMemcpy(&ls,&snes,sizeof(snes));CHKERRQ(ierr); 5020 ierr = PetscMemcpy(&lx,&x,sizeof(x));CHKERRQ(ierr); 5021 prhs[0] = mxCreateDoubleScalar((double)ls); 5022 prhs[1] = mxCreateDoubleScalar((double)it); 5023 prhs[2] = mxCreateDoubleScalar((double)fnorm); 5024 prhs[3] = mxCreateDoubleScalar((double)lx); 5025 prhs[4] = mxCreateString(sctx->funcname); 5026 prhs[5] = sctx->ctx; 5027 ierr = mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscSNESMonitorInternal");CHKERRQ(ierr); 5028 ierr = mxGetScalar(plhs[0]);CHKERRQ(ierr); 5029 mxDestroyArray(prhs[0]); 5030 mxDestroyArray(prhs[1]); 5031 mxDestroyArray(prhs[2]); 5032 mxDestroyArray(prhs[3]); 5033 mxDestroyArray(prhs[4]); 5034 mxDestroyArray(plhs[0]); 5035 PetscFunctionReturn(0); 5036 } 5037 5038 #undef __FUNCT__ 5039 #define __FUNCT__ "SNESMonitorSetMatlab" 5040 /* 5041 SNESMonitorSetMatlab - Sets the monitor function from MATLAB 5042 5043 Level: developer 5044 5045 Developer Note: This bleeds the allocated memory SNESMatlabContext *sctx; 5046 5047 .keywords: SNES, nonlinear, set, function 5048 5049 .seealso: SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESSetFunction() 5050 */ 5051 PetscErrorCode SNESMonitorSetMatlab(SNES snes,const char *SNESMonitorFunction,mxArray *ctx) 5052 { 5053 PetscErrorCode ierr; 5054 SNESMatlabContext *sctx; 5055 5056 PetscFunctionBegin; 5057 /* currently sctx is memory bleed */ 5058 ierr = PetscMalloc(sizeof(SNESMatlabContext),&sctx);CHKERRQ(ierr); 5059 ierr = PetscStrallocpy(SNESMonitorFunction,&sctx->funcname);CHKERRQ(ierr); 5060 /* 5061 This should work, but it doesn't 5062 sctx->ctx = ctx; 5063 mexMakeArrayPersistent(sctx->ctx); 5064 */ 5065 sctx->ctx = mxDuplicateArray(ctx); 5066 ierr = SNESMonitorSet(snes,SNESMonitor_Matlab,sctx,NULL);CHKERRQ(ierr); 5067 PetscFunctionReturn(0); 5068 } 5069 5070 #endif 5071