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