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