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