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