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