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