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