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