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