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