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