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