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