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