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