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