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