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