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