1 2 #include <private/snesimpl.h> /*I "petscsnes.h" I*/ 3 #include <petscdmshell.h> /*I "petscdmshell.h" I*/ 4 5 PetscBool SNESRegisterAllCalled = PETSC_FALSE; 6 PetscFList SNESList = PETSC_NULL; 7 8 /* Logging support */ 9 PetscClassId SNES_CLASSID; 10 PetscLogEvent SNES_Solve, SNES_FunctionEval, SNES_JacobianEval, SNES_GSEval; 11 12 #undef __FUNCT__ 13 #define __FUNCT__ "SNESDMComputeJacobian" 14 /* 15 Translates from a SNES call to a DM call in computing a Jacobian 16 17 This is a legacy calling sequence, should transition to dispatching through the SNESDM. 18 */ 19 PetscErrorCode SNESDMComputeJacobian(SNES snes,Vec X,Mat *J,Mat *B,MatStructure *flag,void *ptr) 20 { 21 PetscErrorCode ierr; 22 DM dm; 23 24 PetscFunctionBegin; 25 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 26 ierr = DMComputeJacobian(dm,X,*J,*B,flag);CHKERRQ(ierr); 27 PetscFunctionReturn(0); 28 } 29 30 #undef __FUNCT__ 31 #define __FUNCT__ "SNESSetErrorIfNotConverged" 32 /*@ 33 SNESSetErrorIfNotConverged - Causes SNESSolve() to generate an error if the solver has not converged. 34 35 Logically Collective on SNES 36 37 Input Parameters: 38 + snes - iterative context obtained from SNESCreate() 39 - flg - PETSC_TRUE indicates you want the error generated 40 41 Options database keys: 42 . -snes_error_if_not_converged : this takes an optional truth value (0/1/no/yes/true/false) 43 44 Level: intermediate 45 46 Notes: 47 Normally PETSc continues if a linear solver fails to converge, you can call SNESGetConvergedReason() after a SNESSolve() 48 to determine if it has converged. 49 50 .keywords: SNES, set, initial guess, nonzero 51 52 .seealso: SNESGetErrorIfNotConverged(), KSPGetErrorIfNotConverged(), KSPSetErrorIFNotConverged() 53 @*/ 54 PetscErrorCode SNESSetErrorIfNotConverged(SNES snes,PetscBool flg) 55 { 56 PetscFunctionBegin; 57 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 58 PetscValidLogicalCollectiveBool(snes,flg,2); 59 snes->errorifnotconverged = flg; 60 61 PetscFunctionReturn(0); 62 } 63 64 #undef __FUNCT__ 65 #define __FUNCT__ "SNESGetErrorIfNotConverged" 66 /*@ 67 SNESGetErrorIfNotConverged - Will SNESSolve() generate an error if the solver does not converge? 68 69 Not Collective 70 71 Input Parameter: 72 . snes - iterative context obtained from SNESCreate() 73 74 Output Parameter: 75 . flag - PETSC_TRUE if it will generate an error, else PETSC_FALSE 76 77 Level: intermediate 78 79 .keywords: SNES, set, initial guess, nonzero 80 81 .seealso: SNESSetErrorIfNotConverged(), KSPGetErrorIfNotConverged(), KSPSetErrorIFNotConverged() 82 @*/ 83 PetscErrorCode SNESGetErrorIfNotConverged(SNES snes,PetscBool *flag) 84 { 85 PetscFunctionBegin; 86 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 87 PetscValidPointer(flag,2); 88 *flag = snes->errorifnotconverged; 89 PetscFunctionReturn(0); 90 } 91 92 #undef __FUNCT__ 93 #define __FUNCT__ "SNESSetFunctionDomainError" 94 /*@ 95 SNESSetFunctionDomainError - tells SNES that the input vector to your FormFunction is not 96 in the functions domain. For example, negative pressure. 97 98 Logically Collective on SNES 99 100 Input Parameters: 101 . snes - the SNES context 102 103 Level: advanced 104 105 .keywords: SNES, view 106 107 .seealso: SNESCreate(), SNESSetFunction() 108 @*/ 109 PetscErrorCode SNESSetFunctionDomainError(SNES snes) 110 { 111 PetscFunctionBegin; 112 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 113 snes->domainerror = PETSC_TRUE; 114 PetscFunctionReturn(0); 115 } 116 117 118 #undef __FUNCT__ 119 #define __FUNCT__ "SNESGetFunctionDomainError" 120 /*@ 121 SNESGetFunctionDomainError - Gets the status of the domain error after a call to SNESComputeFunction; 122 123 Logically Collective on SNES 124 125 Input Parameters: 126 . snes - the SNES context 127 128 Output Parameters: 129 . domainerror Set to PETSC_TRUE if there's a domain error; PETSC_FALSE otherwise. 130 131 Level: advanced 132 133 .keywords: SNES, view 134 135 .seealso: SNESSetFunctionDomainError, SNESComputeFunction() 136 @*/ 137 PetscErrorCode SNESGetFunctionDomainError(SNES snes, PetscBool *domainerror) 138 { 139 PetscFunctionBegin; 140 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 141 PetscValidPointer(domainerror, 2); 142 *domainerror = snes->domainerror; 143 PetscFunctionReturn(0); 144 } 145 146 147 #undef __FUNCT__ 148 #define __FUNCT__ "SNESView" 149 /*@C 150 SNESView - Prints the SNES data structure. 151 152 Collective on SNES 153 154 Input Parameters: 155 + SNES - the SNES context 156 - viewer - visualization context 157 158 Options Database Key: 159 . -snes_view - Calls SNESView() at end of SNESSolve() 160 161 Notes: 162 The available visualization contexts include 163 + PETSC_VIEWER_STDOUT_SELF - standard output (default) 164 - PETSC_VIEWER_STDOUT_WORLD - synchronized standard 165 output where only the first processor opens 166 the file. All other processors send their 167 data to the first processor to print. 168 169 The user can open an alternative visualization context with 170 PetscViewerASCIIOpen() - output to a specified file. 171 172 Level: beginner 173 174 .keywords: SNES, view 175 176 .seealso: PetscViewerASCIIOpen() 177 @*/ 178 PetscErrorCode SNESView(SNES snes,PetscViewer viewer) 179 { 180 SNESKSPEW *kctx; 181 PetscErrorCode ierr; 182 KSP ksp; 183 SNESLineSearch linesearch; 184 PetscBool iascii,isstring; 185 186 PetscFunctionBegin; 187 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 188 if (!viewer) { 189 ierr = PetscViewerASCIIGetStdout(((PetscObject)snes)->comm,&viewer);CHKERRQ(ierr); 190 } 191 PetscValidHeaderSpecific(viewer,PETSC_VIEWER_CLASSID,2); 192 PetscCheckSameComm(snes,1,viewer,2); 193 194 ierr = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 195 ierr = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERSTRING,&isstring);CHKERRQ(ierr); 196 if (iascii) { 197 ierr = PetscObjectPrintClassNamePrefixType((PetscObject)snes,viewer,"SNES Object");CHKERRQ(ierr); 198 if (snes->ops->view) { 199 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 200 ierr = (*snes->ops->view)(snes,viewer);CHKERRQ(ierr); 201 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 202 } 203 ierr = PetscViewerASCIIPrintf(viewer," maximum iterations=%D, maximum function evaluations=%D\n",snes->max_its,snes->max_funcs);CHKERRQ(ierr); 204 ierr = PetscViewerASCIIPrintf(viewer," tolerances: relative=%G, absolute=%G, solution=%G\n", 205 snes->rtol,snes->abstol,snes->stol);CHKERRQ(ierr); 206 ierr = PetscViewerASCIIPrintf(viewer," total number of linear solver iterations=%D\n",snes->linear_its);CHKERRQ(ierr); 207 ierr = PetscViewerASCIIPrintf(viewer," total number of function evaluations=%D\n",snes->nfuncs);CHKERRQ(ierr); 208 if (snes->ksp_ewconv) { 209 kctx = (SNESKSPEW *)snes->kspconvctx; 210 if (kctx) { 211 ierr = PetscViewerASCIIPrintf(viewer," Eisenstat-Walker computation of KSP relative tolerance (version %D)\n",kctx->version);CHKERRQ(ierr); 212 ierr = PetscViewerASCIIPrintf(viewer," rtol_0=%G, rtol_max=%G, threshold=%G\n",kctx->rtol_0,kctx->rtol_max,kctx->threshold);CHKERRQ(ierr); 213 ierr = PetscViewerASCIIPrintf(viewer," gamma=%G, alpha=%G, alpha2=%G\n",kctx->gamma,kctx->alpha,kctx->alpha2);CHKERRQ(ierr); 214 } 215 } 216 if (snes->lagpreconditioner == -1) { 217 ierr = PetscViewerASCIIPrintf(viewer," Preconditioned is never rebuilt\n");CHKERRQ(ierr); 218 } else if (snes->lagpreconditioner > 1) { 219 ierr = PetscViewerASCIIPrintf(viewer," Preconditioned is rebuilt every %D new Jacobians\n",snes->lagpreconditioner);CHKERRQ(ierr); 220 } 221 if (snes->lagjacobian == -1) { 222 ierr = PetscViewerASCIIPrintf(viewer," Jacobian is never rebuilt\n");CHKERRQ(ierr); 223 } else if (snes->lagjacobian > 1) { 224 ierr = PetscViewerASCIIPrintf(viewer," Jacobian is rebuilt every %D SNES iterations\n",snes->lagjacobian);CHKERRQ(ierr); 225 } 226 } else if (isstring) { 227 const char *type; 228 ierr = SNESGetType(snes,&type);CHKERRQ(ierr); 229 ierr = PetscViewerStringSPrintf(viewer," %-3.3s",type);CHKERRQ(ierr); 230 } 231 if (snes->pc && snes->usespc) { 232 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 233 ierr = SNESView(snes->pc, viewer);CHKERRQ(ierr); 234 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 235 } 236 if (snes->usesksp) { 237 ierr = SNESGetKSP(snes,&ksp);CHKERRQ(ierr); 238 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 239 ierr = KSPView(ksp,viewer);CHKERRQ(ierr); 240 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 241 } 242 if (snes->linesearch) { 243 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 244 ierr = SNESGetSNESLineSearch(snes, &linesearch);CHKERRQ(ierr); 245 ierr = SNESLineSearchView(linesearch, viewer);CHKERRQ(ierr); 246 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 247 } 248 PetscFunctionReturn(0); 249 } 250 251 /* 252 We retain a list of functions that also take SNES command 253 line options. These are called at the end SNESSetFromOptions() 254 */ 255 #define MAXSETFROMOPTIONS 5 256 static PetscInt numberofsetfromoptions; 257 static PetscErrorCode (*othersetfromoptions[MAXSETFROMOPTIONS])(SNES); 258 259 #undef __FUNCT__ 260 #define __FUNCT__ "SNESAddOptionsChecker" 261 /*@C 262 SNESAddOptionsChecker - Adds an additional function to check for SNES options. 263 264 Not Collective 265 266 Input Parameter: 267 . snescheck - function that checks for options 268 269 Level: developer 270 271 .seealso: SNESSetFromOptions() 272 @*/ 273 PetscErrorCode SNESAddOptionsChecker(PetscErrorCode (*snescheck)(SNES)) 274 { 275 PetscFunctionBegin; 276 if (numberofsetfromoptions >= MAXSETFROMOPTIONS) { 277 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE, "Too many options checkers, only %D allowed", MAXSETFROMOPTIONS); 278 } 279 othersetfromoptions[numberofsetfromoptions++] = snescheck; 280 PetscFunctionReturn(0); 281 } 282 283 extern PetscErrorCode SNESDefaultMatrixFreeCreate2(SNES,Vec,Mat*); 284 285 #undef __FUNCT__ 286 #define __FUNCT__ "SNESSetUpMatrixFree_Private" 287 static PetscErrorCode SNESSetUpMatrixFree_Private(SNES snes, PetscBool hasOperator, PetscInt version) 288 { 289 Mat J; 290 KSP ksp; 291 PC pc; 292 PetscBool match; 293 PetscErrorCode ierr; 294 295 PetscFunctionBegin; 296 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 297 298 if(!snes->vec_func && (snes->jacobian || snes->jacobian_pre)) { 299 Mat A = snes->jacobian, B = snes->jacobian_pre; 300 ierr = MatGetVecs(A ? A : B, PETSC_NULL,&snes->vec_func);CHKERRQ(ierr); 301 } 302 303 if (version == 1) { 304 ierr = MatCreateSNESMF(snes,&J);CHKERRQ(ierr); 305 ierr = MatMFFDSetOptionsPrefix(J,((PetscObject)snes)->prefix);CHKERRQ(ierr); 306 ierr = MatSetFromOptions(J);CHKERRQ(ierr); 307 } else if (version == 2) { 308 if (!snes->vec_func) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"SNESSetFunction() must be called first"); 309 #if !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_REAL_SINGLE) && !defined(PETSC_USE_REAL___FLOAT128) 310 ierr = SNESDefaultMatrixFreeCreate2(snes,snes->vec_func,&J);CHKERRQ(ierr); 311 #else 312 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP, "matrix-free operator rutines (version 2)"); 313 #endif 314 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE, "matrix-free operator rutines, only version 1 and 2"); 315 316 ierr = PetscInfo1(snes,"Setting default matrix-free operator routines (version %D)\n", version);CHKERRQ(ierr); 317 if (hasOperator) { 318 /* This version replaces the user provided Jacobian matrix with a 319 matrix-free version but still employs the user-provided preconditioner matrix. */ 320 ierr = SNESSetJacobian(snes,J,0,0,0);CHKERRQ(ierr); 321 } else { 322 /* This version replaces both the user-provided Jacobian and the user- 323 provided preconditioner matrix with the default matrix free version. */ 324 void *functx; 325 ierr = SNESGetFunction(snes,PETSC_NULL,PETSC_NULL,&functx);CHKERRQ(ierr); 326 ierr = SNESSetJacobian(snes,J,J,MatMFFDComputeJacobian,functx);CHKERRQ(ierr); 327 /* Force no preconditioner */ 328 ierr = SNESGetKSP(snes,&ksp);CHKERRQ(ierr); 329 ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr); 330 ierr = PetscTypeCompare((PetscObject)pc,PCSHELL,&match);CHKERRQ(ierr); 331 if (!match) { 332 ierr = PetscInfo(snes,"Setting default matrix-free preconditioner routines\nThat is no preconditioner is being used\n");CHKERRQ(ierr); 333 ierr = PCSetType(pc,PCNONE);CHKERRQ(ierr); 334 } 335 } 336 ierr = MatDestroy(&J);CHKERRQ(ierr); 337 PetscFunctionReturn(0); 338 } 339 340 #undef __FUNCT__ 341 #define __FUNCT__ "KSPComputeOperators_SNES" 342 static PetscErrorCode KSPComputeOperators_SNES(KSP ksp,Mat A,Mat B,MatStructure *mstruct,void *ctx) 343 { 344 SNES snes = (SNES)ctx; 345 PetscErrorCode ierr; 346 Mat Asave = A,Bsave = B; 347 348 PetscFunctionBegin; 349 ierr = SNESComputeJacobian(snes,snes->vec_sol,&A,&B,mstruct);CHKERRQ(ierr); 350 if (A != Asave || B != Bsave) SETERRQ(((PetscObject)snes)->comm,PETSC_ERR_SUP,"No support for changing matrices at this time"); 351 PetscFunctionReturn(0); 352 } 353 354 #undef __FUNCT__ 355 #define __FUNCT__ "SNESSetUpMatrices" 356 /*@ 357 SNESSetUpMatrices - ensures that matrices are available for SNES, to be called by SNESSetUp_XXX() 358 359 Collective 360 361 Input Arguments: 362 . snes - snes to configure 363 364 Level: developer 365 366 .seealso: SNESSetUp() 367 @*/ 368 PetscErrorCode SNESSetUpMatrices(SNES snes) 369 { 370 PetscErrorCode ierr; 371 DM dm; 372 SNESDM sdm; 373 374 PetscFunctionBegin; 375 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 376 ierr = DMSNESGetContext(dm,&sdm);CHKERRQ(ierr); 377 if (!sdm->computejacobian) { 378 Mat J,B; 379 ierr = DMCreateMatrix(snes->dm,MATAIJ,&B);CHKERRQ(ierr); 380 if (snes->mf_operator) { 381 ierr = MatCreateSNESMF(snes,&J);CHKERRQ(ierr); 382 ierr = MatMFFDSetOptionsPrefix(J,((PetscObject)snes)->prefix);CHKERRQ(ierr); 383 ierr = MatSetFromOptions(J);CHKERRQ(ierr); 384 } else { 385 J = B; 386 ierr = PetscObjectReference((PetscObject)J);CHKERRQ(ierr); 387 } 388 ierr = SNESSetJacobian(snes,J,B,SNESDMComputeJacobian,PETSC_NULL);CHKERRQ(ierr); 389 ierr = MatDestroy(&J);CHKERRQ(ierr); 390 ierr = MatDestroy(&B);CHKERRQ(ierr); 391 } else if (!snes->jacobian && sdm->computejacobian == MatMFFDComputeJacobian) { 392 Mat J; 393 void *functx; 394 ierr = MatCreateSNESMF(snes,&J);CHKERRQ(ierr); 395 ierr = MatMFFDSetOptionsPrefix(J,((PetscObject)snes)->prefix);CHKERRQ(ierr); 396 ierr = MatSetFromOptions(J);CHKERRQ(ierr); 397 ierr = SNESGetFunction(snes,PETSC_NULL,PETSC_NULL,&functx);CHKERRQ(ierr); 398 ierr = SNESSetJacobian(snes,J,J,MatMFFDComputeJacobian,functx);CHKERRQ(ierr); 399 ierr = MatDestroy(&J);CHKERRQ(ierr); 400 } else if (snes->mf_operator && !snes->jacobian_pre && !snes->jacobian) { 401 Mat J,B; 402 void *functx; 403 ierr = MatCreateSNESMF(snes,&J);CHKERRQ(ierr); 404 ierr = MatMFFDSetOptionsPrefix(J,((PetscObject)snes)->prefix);CHKERRQ(ierr); 405 ierr = MatSetFromOptions(J);CHKERRQ(ierr); 406 ierr = DMCreateMatrix(snes->dm,MATAIJ,&B);CHKERRQ(ierr); 407 ierr = SNESGetFunction(snes,PETSC_NULL,PETSC_NULL,&functx);CHKERRQ(ierr); 408 ierr = SNESSetJacobian(snes,J,B,SNESDMComputeJacobian,functx);CHKERRQ(ierr); 409 ierr = MatDestroy(&J);CHKERRQ(ierr); 410 ierr = MatDestroy(&B);CHKERRQ(ierr); 411 } else if (!snes->jacobian_pre) { 412 Mat J,B; 413 J = snes->jacobian; 414 ierr = DMCreateMatrix(snes->dm,MATAIJ,&B);CHKERRQ(ierr); 415 ierr = SNESSetJacobian(snes,J?J:B,B,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr); 416 ierr = MatDestroy(&B);CHKERRQ(ierr); 417 } 418 { 419 PetscBool flg = PETSC_FALSE; 420 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_kspcompute",&flg,PETSC_NULL);CHKERRQ(ierr); 421 if (flg) { /* Plan to transition to this model */ 422 KSP ksp; 423 ierr = SNESGetKSP(snes,&ksp);CHKERRQ(ierr); 424 ierr = KSPSetComputeOperators(ksp,KSPComputeOperators_SNES,snes);CHKERRQ(ierr); 425 } 426 } 427 PetscFunctionReturn(0); 428 } 429 430 #undef __FUNCT__ 431 #define __FUNCT__ "SNESSetFromOptions" 432 /*@ 433 SNESSetFromOptions - Sets various SNES and KSP parameters from user options. 434 435 Collective on SNES 436 437 Input Parameter: 438 . snes - the SNES context 439 440 Options Database Keys: 441 + -snes_type <type> - ls, tr, ngmres, ncg, richardson, qn, vi, fas 442 . -snes_stol - convergence tolerance in terms of the norm 443 of the change in the solution between steps 444 . -snes_atol <abstol> - absolute tolerance of residual norm 445 . -snes_rtol <rtol> - relative decrease in tolerance norm from initial 446 . -snes_max_it <max_it> - maximum number of iterations 447 . -snes_max_funcs <max_funcs> - maximum number of function evaluations 448 . -snes_max_fail <max_fail> - maximum number of line search failures allowed before stopping, default is none 449 . -snes_max_linear_solve_fail - number of linear solver failures before SNESSolve() stops 450 . -snes_lag_preconditioner <lag> - how often preconditioner is rebuilt (use -1 to never rebuild) 451 . -snes_lag_jacobian <lag> - how often Jacobian is rebuilt (use -1 to never rebuild) 452 . -snes_trtol <trtol> - trust region tolerance 453 . -snes_no_convergence_test - skip convergence test in nonlinear 454 solver; hence iterations will continue until max_it 455 or some other criterion is reached. Saves expense 456 of convergence test 457 . -snes_monitor <optional filename> - prints residual norm at each iteration. if no 458 filename given prints to stdout 459 . -snes_monitor_solution - plots solution at each iteration 460 . -snes_monitor_residual - plots residual (not its norm) at each iteration 461 . -snes_monitor_solution_update - plots update to solution at each iteration 462 . -snes_monitor_draw - plots residual norm at each iteration 463 . -snes_fd - use finite differences to compute Jacobian; very slow, only for testing 464 . -snes_mf_ksp_monitor - if using matrix-free multiply then print h at each KSP iteration 465 - -snes_converged_reason - print the reason for convergence/divergence after each solve 466 467 Options Database for Eisenstat-Walker method: 468 + -snes_ksp_ew - use Eisenstat-Walker method for determining linear system convergence 469 . -snes_ksp_ew_version ver - version of Eisenstat-Walker method 470 . -snes_ksp_ew_rtol0 <rtol0> - Sets rtol0 471 . -snes_ksp_ew_rtolmax <rtolmax> - Sets rtolmax 472 . -snes_ksp_ew_gamma <gamma> - Sets gamma 473 . -snes_ksp_ew_alpha <alpha> - Sets alpha 474 . -snes_ksp_ew_alpha2 <alpha2> - Sets alpha2 475 - -snes_ksp_ew_threshold <threshold> - Sets threshold 476 477 Notes: 478 To see all options, run your program with the -help option or consult 479 the <A href="../../docs/manual.pdf#nameddest=ch_snes">SNES chapter of the users manual</A>. 480 481 Level: beginner 482 483 .keywords: SNES, nonlinear, set, options, database 484 485 .seealso: SNESSetOptionsPrefix() 486 @*/ 487 PetscErrorCode SNESSetFromOptions(SNES snes) 488 { 489 PetscBool flg,mf,mf_operator,pcset; 490 PetscInt i,indx,lag,mf_version,grids; 491 MatStructure matflag; 492 const char *deft = SNESLS; 493 const char *convtests[] = {"default","skip"}; 494 SNESKSPEW *kctx = NULL; 495 char type[256], monfilename[PETSC_MAX_PATH_LEN]; 496 const char *optionsprefix; 497 PetscViewer monviewer; 498 PetscErrorCode ierr; 499 500 PetscFunctionBegin; 501 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 502 503 if (!SNESRegisterAllCalled) {ierr = SNESRegisterAll(PETSC_NULL);CHKERRQ(ierr);} 504 ierr = PetscObjectOptionsBegin((PetscObject)snes);CHKERRQ(ierr); 505 if (((PetscObject)snes)->type_name) { deft = ((PetscObject)snes)->type_name; } 506 ierr = PetscOptionsList("-snes_type","Nonlinear solver method","SNESSetType",SNESList,deft,type,256,&flg);CHKERRQ(ierr); 507 if (flg) { 508 ierr = SNESSetType(snes,type);CHKERRQ(ierr); 509 } else if (!((PetscObject)snes)->type_name) { 510 ierr = SNESSetType(snes,deft);CHKERRQ(ierr); 511 } 512 /* not used here, but called so will go into help messaage */ 513 ierr = PetscOptionsName("-snes_view","Print detailed information on solver used","SNESView",0);CHKERRQ(ierr); 514 515 ierr = PetscOptionsReal("-snes_stol","Stop if step length less than","SNESSetTolerances",snes->stol,&snes->stol,0);CHKERRQ(ierr); 516 ierr = PetscOptionsReal("-snes_atol","Stop if function norm less than","SNESSetTolerances",snes->abstol,&snes->abstol,0);CHKERRQ(ierr); 517 518 ierr = PetscOptionsReal("-snes_rtol","Stop if decrease in function norm less than","SNESSetTolerances",snes->rtol,&snes->rtol,0);CHKERRQ(ierr); 519 ierr = PetscOptionsInt("-snes_max_it","Maximum iterations","SNESSetTolerances",snes->max_its,&snes->max_its,PETSC_NULL);CHKERRQ(ierr); 520 ierr = PetscOptionsInt("-snes_max_funcs","Maximum function evaluations","SNESSetTolerances",snes->max_funcs,&snes->max_funcs,PETSC_NULL);CHKERRQ(ierr); 521 ierr = PetscOptionsInt("-snes_max_fail","Maximum failures","SNESSetTolerances",snes->maxFailures,&snes->maxFailures,PETSC_NULL);CHKERRQ(ierr); 522 ierr = PetscOptionsInt("-snes_max_linear_solve_fail","Maximum failures in linear solves allowed","SNESSetMaxLinearSolveFailures",snes->maxLinearSolveFailures,&snes->maxLinearSolveFailures,PETSC_NULL);CHKERRQ(ierr); 523 ierr = PetscOptionsBool("-snes_error_if_not_converged","Generate error if solver does not converge","SNESSetErrorIfNotConverged",snes->errorifnotconverged,&snes->errorifnotconverged,PETSC_NULL);CHKERRQ(ierr); 524 525 ierr = PetscOptionsInt("-snes_lag_preconditioner","How often to rebuild preconditioner","SNESSetLagPreconditioner",snes->lagpreconditioner,&lag,&flg);CHKERRQ(ierr); 526 if (flg) { 527 ierr = SNESSetLagPreconditioner(snes,lag);CHKERRQ(ierr); 528 } 529 ierr = PetscOptionsInt("-snes_lag_jacobian","How often to rebuild Jacobian","SNESSetLagJacobian",snes->lagjacobian,&lag,&flg);CHKERRQ(ierr); 530 if (flg) { 531 ierr = SNESSetLagJacobian(snes,lag);CHKERRQ(ierr); 532 } 533 ierr = PetscOptionsInt("-snes_grid_sequence","Use grid sequencing to generate initial guess","SNESSetGridSequence",snes->gridsequence,&grids,&flg);CHKERRQ(ierr); 534 if (flg) { 535 ierr = SNESSetGridSequence(snes,grids);CHKERRQ(ierr); 536 } 537 538 ierr = PetscOptionsEList("-snes_convergence_test","Convergence test","SNESSetConvergenceTest",convtests,2,"default",&indx,&flg);CHKERRQ(ierr); 539 if (flg) { 540 switch (indx) { 541 case 0: ierr = SNESSetConvergenceTest(snes,SNESDefaultConverged,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr); break; 542 case 1: ierr = SNESSetConvergenceTest(snes,SNESSkipConverged,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr); break; 543 } 544 } 545 546 ierr = PetscOptionsBool("-snes_converged_reason","Print reason for converged or diverged","SNESSolve",snes->printreason,&snes->printreason,PETSC_NULL);CHKERRQ(ierr); 547 548 kctx = (SNESKSPEW *)snes->kspconvctx; 549 550 ierr = PetscOptionsBool("-snes_ksp_ew","Use Eisentat-Walker linear system convergence test","SNESKSPSetUseEW",snes->ksp_ewconv,&snes->ksp_ewconv,PETSC_NULL);CHKERRQ(ierr); 551 552 ierr = PetscOptionsInt("-snes_ksp_ew_version","Version 1, 2 or 3","SNESKSPSetParametersEW",kctx->version,&kctx->version,0);CHKERRQ(ierr); 553 ierr = PetscOptionsReal("-snes_ksp_ew_rtol0","0 <= rtol0 < 1","SNESKSPSetParametersEW",kctx->rtol_0,&kctx->rtol_0,0);CHKERRQ(ierr); 554 ierr = PetscOptionsReal("-snes_ksp_ew_rtolmax","0 <= rtolmax < 1","SNESKSPSetParametersEW",kctx->rtol_max,&kctx->rtol_max,0);CHKERRQ(ierr); 555 ierr = PetscOptionsReal("-snes_ksp_ew_gamma","0 <= gamma <= 1","SNESKSPSetParametersEW",kctx->gamma,&kctx->gamma,0);CHKERRQ(ierr); 556 ierr = PetscOptionsReal("-snes_ksp_ew_alpha","1 < alpha <= 2","SNESKSPSetParametersEW",kctx->alpha,&kctx->alpha,0);CHKERRQ(ierr); 557 ierr = PetscOptionsReal("-snes_ksp_ew_alpha2","alpha2","SNESKSPSetParametersEW",kctx->alpha2,&kctx->alpha2,0);CHKERRQ(ierr); 558 ierr = PetscOptionsReal("-snes_ksp_ew_threshold","0 < threshold < 1","SNESKSPSetParametersEW",kctx->threshold,&kctx->threshold,0);CHKERRQ(ierr); 559 560 flg = PETSC_FALSE; 561 ierr = PetscOptionsBool("-snes_monitor_cancel","Remove all monitors","SNESMonitorCancel",flg,&flg,PETSC_NULL);CHKERRQ(ierr); 562 if (flg) {ierr = SNESMonitorCancel(snes);CHKERRQ(ierr);} 563 564 ierr = PetscOptionsString("-snes_monitor","Monitor norm of function","SNESMonitorSet","stdout",monfilename,PETSC_MAX_PATH_LEN,&flg);CHKERRQ(ierr); 565 if (flg) { 566 ierr = PetscViewerASCIIOpen(((PetscObject)snes)->comm,monfilename,&monviewer);CHKERRQ(ierr); 567 ierr = SNESMonitorSet(snes,SNESMonitorDefault,monviewer,(PetscErrorCode (*)(void**))PetscViewerDestroy);CHKERRQ(ierr); 568 } 569 570 ierr = PetscOptionsString("-snes_monitor_range","Monitor range of elements of function","SNESMonitorSet","stdout",monfilename,PETSC_MAX_PATH_LEN,&flg);CHKERRQ(ierr); 571 if (flg) { 572 ierr = SNESMonitorSet(snes,SNESMonitorRange,0,0);CHKERRQ(ierr); 573 } 574 575 ierr = PetscOptionsString("-snes_ratiomonitor","Monitor ratios of norms of function","SNESMonitorSetRatio","stdout",monfilename,PETSC_MAX_PATH_LEN,&flg);CHKERRQ(ierr); 576 if (flg) { 577 ierr = PetscViewerASCIIOpen(((PetscObject)snes)->comm,monfilename,&monviewer);CHKERRQ(ierr); 578 ierr = SNESMonitorSetRatio(snes,monviewer);CHKERRQ(ierr); 579 } 580 581 ierr = PetscOptionsString("-snes_monitor_short","Monitor norm of function (fewer digits)","SNESMonitorSet","stdout",monfilename,PETSC_MAX_PATH_LEN,&flg);CHKERRQ(ierr); 582 if (flg) { 583 ierr = PetscViewerASCIIOpen(((PetscObject)snes)->comm,monfilename,&monviewer);CHKERRQ(ierr); 584 ierr = SNESMonitorSet(snes,SNESMonitorDefaultShort,monviewer,(PetscErrorCode (*)(void**))PetscViewerDestroy);CHKERRQ(ierr); 585 } 586 587 ierr = PetscOptionsString("-snes_monitor_python","Use Python function","SNESMonitorSet",0,monfilename,PETSC_MAX_PATH_LEN,&flg);CHKERRQ(ierr); 588 if (flg) {ierr = PetscPythonMonitorSet((PetscObject)snes,monfilename);CHKERRQ(ierr);} 589 590 flg = PETSC_FALSE; 591 ierr = PetscOptionsBool("-snes_monitor_solution","Plot solution at each iteration","SNESMonitorSolution",flg,&flg,PETSC_NULL);CHKERRQ(ierr); 592 if (flg) {ierr = SNESMonitorSet(snes,SNESMonitorSolution,0,0);CHKERRQ(ierr);} 593 flg = PETSC_FALSE; 594 ierr = PetscOptionsBool("-snes_monitor_solution_update","Plot correction at each iteration","SNESMonitorSolutionUpdate",flg,&flg,PETSC_NULL);CHKERRQ(ierr); 595 if (flg) {ierr = SNESMonitorSet(snes,SNESMonitorSolutionUpdate,0,0);CHKERRQ(ierr);} 596 flg = PETSC_FALSE; 597 ierr = PetscOptionsBool("-snes_monitor_residual","Plot residual at each iteration","SNESMonitorResidual",flg,&flg,PETSC_NULL);CHKERRQ(ierr); 598 if (flg) {ierr = SNESMonitorSet(snes,SNESMonitorResidual,0,0);CHKERRQ(ierr);} 599 flg = PETSC_FALSE; 600 ierr = PetscOptionsBool("-snes_monitor_draw","Plot function norm at each iteration","SNESMonitorLG",flg,&flg,PETSC_NULL);CHKERRQ(ierr); 601 if (flg) {ierr = SNESMonitorSet(snes,SNESMonitorLG,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr);} 602 flg = PETSC_FALSE; 603 ierr = PetscOptionsBool("-snes_monitor_range_draw","Plot function range at each iteration","SNESMonitorLG",flg,&flg,PETSC_NULL);CHKERRQ(ierr); 604 if (flg) {ierr = SNESMonitorSet(snes,SNESMonitorLGRange,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr);} 605 606 flg = PETSC_FALSE; 607 ierr = PetscOptionsBool("-snes_fd","Use finite differences (slow) to compute Jacobian","SNESDefaultComputeJacobian",flg,&flg,PETSC_NULL);CHKERRQ(ierr); 608 if (flg) { 609 void *functx; 610 ierr = SNESGetFunction(snes,PETSC_NULL,PETSC_NULL,&functx);CHKERRQ(ierr); 611 ierr = SNESSetJacobian(snes,snes->jacobian,snes->jacobian_pre,SNESDefaultComputeJacobian,functx);CHKERRQ(ierr); 612 ierr = PetscInfo(snes,"Setting default finite difference Jacobian matrix\n");CHKERRQ(ierr); 613 } 614 615 mf = mf_operator = PETSC_FALSE; 616 flg = PETSC_FALSE; 617 ierr = PetscOptionsBool("-snes_mf_operator","Use a Matrix-Free Jacobian with user-provided preconditioner matrix","MatCreateSNESMF",PETSC_FALSE,&mf_operator,&flg);CHKERRQ(ierr); 618 if (flg && mf_operator) { 619 snes->mf_operator = PETSC_TRUE; 620 mf = PETSC_TRUE; 621 } 622 flg = PETSC_FALSE; 623 ierr = PetscOptionsBool("-snes_mf","Use a Matrix-Free Jacobian with no preconditioner matrix","MatCreateSNESMF",PETSC_FALSE,&mf,&flg);CHKERRQ(ierr); 624 if (!flg && mf_operator) mf = PETSC_TRUE; 625 mf_version = 1; 626 ierr = PetscOptionsInt("-snes_mf_version","Matrix-Free routines version 1 or 2","None",mf_version,&mf_version,0);CHKERRQ(ierr); 627 628 629 /* GS Options */ 630 ierr = PetscOptionsInt("-snes_gs_sweeps","Number of sweeps of GS to apply","SNESComputeGS",snes->gssweeps,&snes->gssweeps,PETSC_NULL);CHKERRQ(ierr); 631 632 for(i = 0; i < numberofsetfromoptions; i++) { 633 ierr = (*othersetfromoptions[i])(snes);CHKERRQ(ierr); 634 } 635 636 if (snes->ops->setfromoptions) { 637 ierr = (*snes->ops->setfromoptions)(snes);CHKERRQ(ierr); 638 } 639 640 /* process any options handlers added with PetscObjectAddOptionsHandler() */ 641 ierr = PetscObjectProcessOptionsHandlers((PetscObject)snes);CHKERRQ(ierr); 642 ierr = PetscOptionsEnd();CHKERRQ(ierr); 643 644 if (mf) { ierr = SNESSetUpMatrixFree_Private(snes, mf_operator, mf_version);CHKERRQ(ierr); } 645 646 if (!snes->ksp) {ierr = SNESGetKSP(snes,&snes->ksp);CHKERRQ(ierr);} 647 ierr = KSPGetOperators(snes->ksp,PETSC_NULL,PETSC_NULL,&matflag); 648 ierr = KSPSetOperators(snes->ksp,snes->jacobian,snes->jacobian_pre,matflag);CHKERRQ(ierr); 649 ierr = KSPSetFromOptions(snes->ksp);CHKERRQ(ierr); 650 651 if (!snes->linesearch) { 652 ierr = SNESGetSNESLineSearch(snes, &snes->linesearch);CHKERRQ(ierr); 653 } 654 ierr = SNESLineSearchSetFromOptions(snes->linesearch);CHKERRQ(ierr); 655 656 /* if someone has set the SNES PC type, create it. */ 657 ierr = SNESGetOptionsPrefix(snes, &optionsprefix);CHKERRQ(ierr); 658 ierr = PetscOptionsHasName(optionsprefix, "-npc_snes_type", &pcset);CHKERRQ(ierr); 659 if (pcset && (!snes->pc)) { 660 ierr = SNESGetPC(snes, &snes->pc);CHKERRQ(ierr); 661 } 662 if (snes->pc) { 663 ierr = SNESSetOptionsPrefix(snes->pc, optionsprefix);CHKERRQ(ierr); 664 ierr = SNESAppendOptionsPrefix(snes->pc, "npc_");CHKERRQ(ierr); 665 ierr = SNESSetDM(snes->pc, snes->dm);CHKERRQ(ierr); 666 /* 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 ierr = DMCreateInterpolation(snes->dm,fine,&interp,PETSC_NULL);CHKERRQ(ierr); 3326 ierr = DMCreateGlobalVector(fine,&xnew);CHKERRQ(ierr); 3327 ierr = MatInterpolate(interp,x,xnew);CHKERRQ(ierr); 3328 ierr = MatDestroy(&interp);CHKERRQ(ierr); 3329 x = xnew; 3330 3331 ierr = SNESReset(snes);CHKERRQ(ierr); 3332 ierr = SNESSetDM(snes,fine);CHKERRQ(ierr); 3333 ierr = DMDestroy(&fine);CHKERRQ(ierr); 3334 ierr = PetscViewerASCIIPopTab(PETSC_VIEWER_STDOUT_(((PetscObject)snes)->comm));CHKERRQ(ierr); 3335 } 3336 } 3337 ierr = VecDestroy(&xcreated);CHKERRQ(ierr); 3338 PetscFunctionReturn(0); 3339 } 3340 3341 /* --------- Internal routines for SNES Package --------- */ 3342 3343 #undef __FUNCT__ 3344 #define __FUNCT__ "SNESSetType" 3345 /*@C 3346 SNESSetType - Sets the method for the nonlinear solver. 3347 3348 Collective on SNES 3349 3350 Input Parameters: 3351 + snes - the SNES context 3352 - type - a known method 3353 3354 Options Database Key: 3355 . -snes_type <type> - Sets the method; use -help for a list 3356 of available methods (for instance, ls or tr) 3357 3358 Notes: 3359 See "petsc/include/petscsnes.h" for available methods (for instance) 3360 + SNESLS - Newton's method with line search 3361 (systems of nonlinear equations) 3362 . SNESTR - Newton's method with trust region 3363 (systems of nonlinear equations) 3364 3365 Normally, it is best to use the SNESSetFromOptions() command and then 3366 set the SNES solver type from the options database rather than by using 3367 this routine. Using the options database provides the user with 3368 maximum flexibility in evaluating the many nonlinear solvers. 3369 The SNESSetType() routine is provided for those situations where it 3370 is necessary to set the nonlinear solver independently of the command 3371 line or options database. This might be the case, for example, when 3372 the choice of solver changes during the execution of the program, 3373 and the user's application is taking responsibility for choosing the 3374 appropriate method. 3375 3376 Level: intermediate 3377 3378 .keywords: SNES, set, type 3379 3380 .seealso: SNESType, SNESCreate() 3381 3382 @*/ 3383 PetscErrorCode SNESSetType(SNES snes,const SNESType type) 3384 { 3385 PetscErrorCode ierr,(*r)(SNES); 3386 PetscBool match; 3387 3388 PetscFunctionBegin; 3389 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3390 PetscValidCharPointer(type,2); 3391 3392 ierr = PetscTypeCompare((PetscObject)snes,type,&match);CHKERRQ(ierr); 3393 if (match) PetscFunctionReturn(0); 3394 3395 ierr = PetscFListFind(SNESList,((PetscObject)snes)->comm,type,PETSC_TRUE,(void (**)(void)) &r);CHKERRQ(ierr); 3396 if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unable to find requested SNES type %s",type); 3397 /* Destroy the previous private SNES context */ 3398 if (snes->ops->destroy) { 3399 ierr = (*(snes)->ops->destroy)(snes);CHKERRQ(ierr); 3400 snes->ops->destroy = PETSC_NULL; 3401 } 3402 /* Reinitialize function pointers in SNESOps structure */ 3403 snes->ops->setup = 0; 3404 snes->ops->solve = 0; 3405 snes->ops->view = 0; 3406 snes->ops->setfromoptions = 0; 3407 snes->ops->destroy = 0; 3408 /* Call the SNESCreate_XXX routine for this particular Nonlinear solver */ 3409 snes->setupcalled = PETSC_FALSE; 3410 ierr = PetscObjectChangeTypeName((PetscObject)snes,type);CHKERRQ(ierr); 3411 ierr = (*r)(snes);CHKERRQ(ierr); 3412 #if defined(PETSC_HAVE_AMS) 3413 if (PetscAMSPublishAll) { 3414 ierr = PetscObjectAMSPublish((PetscObject)snes);CHKERRQ(ierr); 3415 } 3416 #endif 3417 PetscFunctionReturn(0); 3418 } 3419 3420 3421 /* --------------------------------------------------------------------- */ 3422 #undef __FUNCT__ 3423 #define __FUNCT__ "SNESRegisterDestroy" 3424 /*@ 3425 SNESRegisterDestroy - Frees the list of nonlinear solvers that were 3426 registered by SNESRegisterDynamic(). 3427 3428 Not Collective 3429 3430 Level: advanced 3431 3432 .keywords: SNES, nonlinear, register, destroy 3433 3434 .seealso: SNESRegisterAll(), SNESRegisterAll() 3435 @*/ 3436 PetscErrorCode SNESRegisterDestroy(void) 3437 { 3438 PetscErrorCode ierr; 3439 3440 PetscFunctionBegin; 3441 ierr = PetscFListDestroy(&SNESList);CHKERRQ(ierr); 3442 SNESRegisterAllCalled = PETSC_FALSE; 3443 PetscFunctionReturn(0); 3444 } 3445 3446 #undef __FUNCT__ 3447 #define __FUNCT__ "SNESGetType" 3448 /*@C 3449 SNESGetType - Gets the SNES method type and name (as a string). 3450 3451 Not Collective 3452 3453 Input Parameter: 3454 . snes - nonlinear solver context 3455 3456 Output Parameter: 3457 . type - SNES method (a character string) 3458 3459 Level: intermediate 3460 3461 .keywords: SNES, nonlinear, get, type, name 3462 @*/ 3463 PetscErrorCode SNESGetType(SNES snes,const SNESType *type) 3464 { 3465 PetscFunctionBegin; 3466 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3467 PetscValidPointer(type,2); 3468 *type = ((PetscObject)snes)->type_name; 3469 PetscFunctionReturn(0); 3470 } 3471 3472 #undef __FUNCT__ 3473 #define __FUNCT__ "SNESGetSolution" 3474 /*@ 3475 SNESGetSolution - Returns the vector where the approximate solution is 3476 stored. This is the fine grid solution when using SNESSetGridSequence(). 3477 3478 Not Collective, but Vec is parallel if SNES is parallel 3479 3480 Input Parameter: 3481 . snes - the SNES context 3482 3483 Output Parameter: 3484 . x - the solution 3485 3486 Level: intermediate 3487 3488 .keywords: SNES, nonlinear, get, solution 3489 3490 .seealso: SNESGetSolutionUpdate(), SNESGetFunction() 3491 @*/ 3492 PetscErrorCode SNESGetSolution(SNES snes,Vec *x) 3493 { 3494 PetscFunctionBegin; 3495 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3496 PetscValidPointer(x,2); 3497 *x = snes->vec_sol; 3498 PetscFunctionReturn(0); 3499 } 3500 3501 #undef __FUNCT__ 3502 #define __FUNCT__ "SNESGetSolutionUpdate" 3503 /*@ 3504 SNESGetSolutionUpdate - Returns the vector where the solution update is 3505 stored. 3506 3507 Not Collective, but Vec is parallel if SNES is parallel 3508 3509 Input Parameter: 3510 . snes - the SNES context 3511 3512 Output Parameter: 3513 . x - the solution update 3514 3515 Level: advanced 3516 3517 .keywords: SNES, nonlinear, get, solution, update 3518 3519 .seealso: SNESGetSolution(), SNESGetFunction() 3520 @*/ 3521 PetscErrorCode SNESGetSolutionUpdate(SNES snes,Vec *x) 3522 { 3523 PetscFunctionBegin; 3524 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3525 PetscValidPointer(x,2); 3526 *x = snes->vec_sol_update; 3527 PetscFunctionReturn(0); 3528 } 3529 3530 #undef __FUNCT__ 3531 #define __FUNCT__ "SNESGetFunction" 3532 /*@C 3533 SNESGetFunction - Returns the vector where the function is stored. 3534 3535 Not Collective, but Vec is parallel if SNES is parallel. Collective if Vec is requested, but has not been created yet. 3536 3537 Input Parameter: 3538 . snes - the SNES context 3539 3540 Output Parameter: 3541 + r - the function (or PETSC_NULL) 3542 . func - the function (or PETSC_NULL) 3543 - ctx - the function context (or PETSC_NULL) 3544 3545 Level: advanced 3546 3547 .keywords: SNES, nonlinear, get, function 3548 3549 .seealso: SNESSetFunction(), SNESGetSolution() 3550 @*/ 3551 PetscErrorCode SNESGetFunction(SNES snes,Vec *r,PetscErrorCode (**func)(SNES,Vec,Vec,void*),void **ctx) 3552 { 3553 PetscErrorCode ierr; 3554 DM dm; 3555 3556 PetscFunctionBegin; 3557 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3558 if (r) { 3559 if (!snes->vec_func) { 3560 if (snes->vec_rhs) { 3561 ierr = VecDuplicate(snes->vec_rhs,&snes->vec_func);CHKERRQ(ierr); 3562 } else if (snes->vec_sol) { 3563 ierr = VecDuplicate(snes->vec_sol,&snes->vec_func);CHKERRQ(ierr); 3564 } else if (snes->dm) { 3565 ierr = DMCreateGlobalVector(snes->dm,&snes->vec_func);CHKERRQ(ierr); 3566 } 3567 } 3568 *r = snes->vec_func; 3569 } 3570 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 3571 ierr = DMSNESGetFunction(dm,func,ctx);CHKERRQ(ierr); 3572 PetscFunctionReturn(0); 3573 } 3574 3575 /*@C 3576 SNESGetGS - Returns the GS function and context. 3577 3578 Input Parameter: 3579 . snes - the SNES context 3580 3581 Output Parameter: 3582 + gsfunc - the function (or PETSC_NULL) 3583 - ctx - the function context (or PETSC_NULL) 3584 3585 Level: advanced 3586 3587 .keywords: SNES, nonlinear, get, function 3588 3589 .seealso: SNESSetGS(), SNESGetFunction() 3590 @*/ 3591 3592 #undef __FUNCT__ 3593 #define __FUNCT__ "SNESGetGS" 3594 PetscErrorCode SNESGetGS (SNES snes, PetscErrorCode(**func)(SNES, Vec, Vec, void*), void ** ctx) 3595 { 3596 PetscErrorCode ierr; 3597 DM dm; 3598 3599 PetscFunctionBegin; 3600 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3601 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 3602 ierr = DMSNESGetGS(dm,func,ctx);CHKERRQ(ierr); 3603 PetscFunctionReturn(0); 3604 } 3605 3606 #undef __FUNCT__ 3607 #define __FUNCT__ "SNESSetOptionsPrefix" 3608 /*@C 3609 SNESSetOptionsPrefix - Sets the prefix used for searching for all 3610 SNES options in the database. 3611 3612 Logically Collective on SNES 3613 3614 Input Parameter: 3615 + snes - the SNES context 3616 - prefix - the prefix to prepend to all option names 3617 3618 Notes: 3619 A hyphen (-) must NOT be given at the beginning of the prefix name. 3620 The first character of all runtime options is AUTOMATICALLY the hyphen. 3621 3622 Level: advanced 3623 3624 .keywords: SNES, set, options, prefix, database 3625 3626 .seealso: SNESSetFromOptions() 3627 @*/ 3628 PetscErrorCode SNESSetOptionsPrefix(SNES snes,const char prefix[]) 3629 { 3630 PetscErrorCode ierr; 3631 3632 PetscFunctionBegin; 3633 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3634 ierr = PetscObjectSetOptionsPrefix((PetscObject)snes,prefix);CHKERRQ(ierr); 3635 if (!snes->ksp) {ierr = SNESGetKSP(snes,&snes->ksp);CHKERRQ(ierr);} 3636 ierr = KSPSetOptionsPrefix(snes->ksp,prefix);CHKERRQ(ierr); 3637 PetscFunctionReturn(0); 3638 } 3639 3640 #undef __FUNCT__ 3641 #define __FUNCT__ "SNESAppendOptionsPrefix" 3642 /*@C 3643 SNESAppendOptionsPrefix - Appends to the prefix used for searching for all 3644 SNES options in the database. 3645 3646 Logically Collective on SNES 3647 3648 Input Parameters: 3649 + snes - the SNES context 3650 - prefix - the prefix to prepend to all option names 3651 3652 Notes: 3653 A hyphen (-) must NOT be given at the beginning of the prefix name. 3654 The first character of all runtime options is AUTOMATICALLY the hyphen. 3655 3656 Level: advanced 3657 3658 .keywords: SNES, append, options, prefix, database 3659 3660 .seealso: SNESGetOptionsPrefix() 3661 @*/ 3662 PetscErrorCode SNESAppendOptionsPrefix(SNES snes,const char prefix[]) 3663 { 3664 PetscErrorCode ierr; 3665 3666 PetscFunctionBegin; 3667 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3668 ierr = PetscObjectAppendOptionsPrefix((PetscObject)snes,prefix);CHKERRQ(ierr); 3669 if (!snes->ksp) {ierr = SNESGetKSP(snes,&snes->ksp);CHKERRQ(ierr);} 3670 ierr = KSPAppendOptionsPrefix(snes->ksp,prefix);CHKERRQ(ierr); 3671 PetscFunctionReturn(0); 3672 } 3673 3674 #undef __FUNCT__ 3675 #define __FUNCT__ "SNESGetOptionsPrefix" 3676 /*@C 3677 SNESGetOptionsPrefix - Sets the prefix used for searching for all 3678 SNES options in the database. 3679 3680 Not Collective 3681 3682 Input Parameter: 3683 . snes - the SNES context 3684 3685 Output Parameter: 3686 . prefix - pointer to the prefix string used 3687 3688 Notes: On the fortran side, the user should pass in a string 'prefix' of 3689 sufficient length to hold the prefix. 3690 3691 Level: advanced 3692 3693 .keywords: SNES, get, options, prefix, database 3694 3695 .seealso: SNESAppendOptionsPrefix() 3696 @*/ 3697 PetscErrorCode SNESGetOptionsPrefix(SNES snes,const char *prefix[]) 3698 { 3699 PetscErrorCode ierr; 3700 3701 PetscFunctionBegin; 3702 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3703 ierr = PetscObjectGetOptionsPrefix((PetscObject)snes,prefix);CHKERRQ(ierr); 3704 PetscFunctionReturn(0); 3705 } 3706 3707 3708 #undef __FUNCT__ 3709 #define __FUNCT__ "SNESRegister" 3710 /*@C 3711 SNESRegister - See SNESRegisterDynamic() 3712 3713 Level: advanced 3714 @*/ 3715 PetscErrorCode SNESRegister(const char sname[],const char path[],const char name[],PetscErrorCode (*function)(SNES)) 3716 { 3717 char fullname[PETSC_MAX_PATH_LEN]; 3718 PetscErrorCode ierr; 3719 3720 PetscFunctionBegin; 3721 ierr = PetscFListConcat(path,name,fullname);CHKERRQ(ierr); 3722 ierr = PetscFListAdd(&SNESList,sname,fullname,(void (*)(void))function);CHKERRQ(ierr); 3723 PetscFunctionReturn(0); 3724 } 3725 3726 #undef __FUNCT__ 3727 #define __FUNCT__ "SNESTestLocalMin" 3728 PetscErrorCode SNESTestLocalMin(SNES snes) 3729 { 3730 PetscErrorCode ierr; 3731 PetscInt N,i,j; 3732 Vec u,uh,fh; 3733 PetscScalar value; 3734 PetscReal norm; 3735 3736 PetscFunctionBegin; 3737 ierr = SNESGetSolution(snes,&u);CHKERRQ(ierr); 3738 ierr = VecDuplicate(u,&uh);CHKERRQ(ierr); 3739 ierr = VecDuplicate(u,&fh);CHKERRQ(ierr); 3740 3741 /* currently only works for sequential */ 3742 ierr = PetscPrintf(PETSC_COMM_WORLD,"Testing FormFunction() for local min\n"); 3743 ierr = VecGetSize(u,&N);CHKERRQ(ierr); 3744 for (i=0; i<N; i++) { 3745 ierr = VecCopy(u,uh);CHKERRQ(ierr); 3746 ierr = PetscPrintf(PETSC_COMM_WORLD,"i = %D\n",i);CHKERRQ(ierr); 3747 for (j=-10; j<11; j++) { 3748 value = PetscSign(j)*exp(PetscAbs(j)-10.0); 3749 ierr = VecSetValue(uh,i,value,ADD_VALUES);CHKERRQ(ierr); 3750 ierr = SNESComputeFunction(snes,uh,fh);CHKERRQ(ierr); 3751 ierr = VecNorm(fh,NORM_2,&norm);CHKERRQ(ierr); 3752 ierr = PetscPrintf(PETSC_COMM_WORLD," j norm %D %18.16e\n",j,norm);CHKERRQ(ierr); 3753 value = -value; 3754 ierr = VecSetValue(uh,i,value,ADD_VALUES);CHKERRQ(ierr); 3755 } 3756 } 3757 ierr = VecDestroy(&uh);CHKERRQ(ierr); 3758 ierr = VecDestroy(&fh);CHKERRQ(ierr); 3759 PetscFunctionReturn(0); 3760 } 3761 3762 #undef __FUNCT__ 3763 #define __FUNCT__ "SNESKSPSetUseEW" 3764 /*@ 3765 SNESKSPSetUseEW - Sets SNES use Eisenstat-Walker method for 3766 computing relative tolerance for linear solvers within an inexact 3767 Newton method. 3768 3769 Logically Collective on SNES 3770 3771 Input Parameters: 3772 + snes - SNES context 3773 - flag - PETSC_TRUE or PETSC_FALSE 3774 3775 Options Database: 3776 + -snes_ksp_ew - use Eisenstat-Walker method for determining linear system convergence 3777 . -snes_ksp_ew_version ver - version of Eisenstat-Walker method 3778 . -snes_ksp_ew_rtol0 <rtol0> - Sets rtol0 3779 . -snes_ksp_ew_rtolmax <rtolmax> - Sets rtolmax 3780 . -snes_ksp_ew_gamma <gamma> - Sets gamma 3781 . -snes_ksp_ew_alpha <alpha> - Sets alpha 3782 . -snes_ksp_ew_alpha2 <alpha2> - Sets alpha2 3783 - -snes_ksp_ew_threshold <threshold> - Sets threshold 3784 3785 Notes: 3786 Currently, the default is to use a constant relative tolerance for 3787 the inner linear solvers. Alternatively, one can use the 3788 Eisenstat-Walker method, where the relative convergence tolerance 3789 is reset at each Newton iteration according progress of the nonlinear 3790 solver. 3791 3792 Level: advanced 3793 3794 Reference: 3795 S. C. Eisenstat and H. F. Walker, "Choosing the forcing terms in an 3796 inexact Newton method", SISC 17 (1), pp.16-32, 1996. 3797 3798 .keywords: SNES, KSP, Eisenstat, Walker, convergence, test, inexact, Newton 3799 3800 .seealso: SNESKSPGetUseEW(), SNESKSPGetParametersEW(), SNESKSPSetParametersEW() 3801 @*/ 3802 PetscErrorCode SNESKSPSetUseEW(SNES snes,PetscBool flag) 3803 { 3804 PetscFunctionBegin; 3805 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3806 PetscValidLogicalCollectiveBool(snes,flag,2); 3807 snes->ksp_ewconv = flag; 3808 PetscFunctionReturn(0); 3809 } 3810 3811 #undef __FUNCT__ 3812 #define __FUNCT__ "SNESKSPGetUseEW" 3813 /*@ 3814 SNESKSPGetUseEW - Gets if SNES is using Eisenstat-Walker method 3815 for computing relative tolerance for linear solvers within an 3816 inexact Newton method. 3817 3818 Not Collective 3819 3820 Input Parameter: 3821 . snes - SNES context 3822 3823 Output Parameter: 3824 . flag - PETSC_TRUE or PETSC_FALSE 3825 3826 Notes: 3827 Currently, the default is to use a constant relative tolerance for 3828 the inner linear solvers. Alternatively, one can use the 3829 Eisenstat-Walker method, where the relative convergence tolerance 3830 is reset at each Newton iteration according progress of the nonlinear 3831 solver. 3832 3833 Level: advanced 3834 3835 Reference: 3836 S. C. Eisenstat and H. F. Walker, "Choosing the forcing terms in an 3837 inexact Newton method", SISC 17 (1), pp.16-32, 1996. 3838 3839 .keywords: SNES, KSP, Eisenstat, Walker, convergence, test, inexact, Newton 3840 3841 .seealso: SNESKSPSetUseEW(), SNESKSPGetParametersEW(), SNESKSPSetParametersEW() 3842 @*/ 3843 PetscErrorCode SNESKSPGetUseEW(SNES snes, PetscBool *flag) 3844 { 3845 PetscFunctionBegin; 3846 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3847 PetscValidPointer(flag,2); 3848 *flag = snes->ksp_ewconv; 3849 PetscFunctionReturn(0); 3850 } 3851 3852 #undef __FUNCT__ 3853 #define __FUNCT__ "SNESKSPSetParametersEW" 3854 /*@ 3855 SNESKSPSetParametersEW - Sets parameters for Eisenstat-Walker 3856 convergence criteria for the linear solvers within an inexact 3857 Newton method. 3858 3859 Logically Collective on SNES 3860 3861 Input Parameters: 3862 + snes - SNES context 3863 . version - version 1, 2 (default is 2) or 3 3864 . rtol_0 - initial relative tolerance (0 <= rtol_0 < 1) 3865 . rtol_max - maximum relative tolerance (0 <= rtol_max < 1) 3866 . gamma - multiplicative factor for version 2 rtol computation 3867 (0 <= gamma2 <= 1) 3868 . alpha - power for version 2 rtol computation (1 < alpha <= 2) 3869 . alpha2 - power for safeguard 3870 - threshold - threshold for imposing safeguard (0 < threshold < 1) 3871 3872 Note: 3873 Version 3 was contributed by Luis Chacon, June 2006. 3874 3875 Use PETSC_DEFAULT to retain the default for any of the parameters. 3876 3877 Level: advanced 3878 3879 Reference: 3880 S. C. Eisenstat and H. F. Walker, "Choosing the forcing terms in an 3881 inexact Newton method", Utah State University Math. Stat. Dept. Res. 3882 Report 6/94/75, June, 1994, to appear in SIAM J. Sci. Comput. 3883 3884 .keywords: SNES, KSP, Eisenstat, Walker, set, parameters 3885 3886 .seealso: SNESKSPSetUseEW(), SNESKSPGetUseEW(), SNESKSPGetParametersEW() 3887 @*/ 3888 PetscErrorCode SNESKSPSetParametersEW(SNES snes,PetscInt version,PetscReal rtol_0,PetscReal rtol_max, 3889 PetscReal gamma,PetscReal alpha,PetscReal alpha2,PetscReal threshold) 3890 { 3891 SNESKSPEW *kctx; 3892 PetscFunctionBegin; 3893 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3894 kctx = (SNESKSPEW*)snes->kspconvctx; 3895 if (!kctx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"No Eisenstat-Walker context existing"); 3896 PetscValidLogicalCollectiveInt(snes,version,2); 3897 PetscValidLogicalCollectiveReal(snes,rtol_0,3); 3898 PetscValidLogicalCollectiveReal(snes,rtol_max,4); 3899 PetscValidLogicalCollectiveReal(snes,gamma,5); 3900 PetscValidLogicalCollectiveReal(snes,alpha,6); 3901 PetscValidLogicalCollectiveReal(snes,alpha2,7); 3902 PetscValidLogicalCollectiveReal(snes,threshold,8); 3903 3904 if (version != PETSC_DEFAULT) kctx->version = version; 3905 if (rtol_0 != PETSC_DEFAULT) kctx->rtol_0 = rtol_0; 3906 if (rtol_max != PETSC_DEFAULT) kctx->rtol_max = rtol_max; 3907 if (gamma != PETSC_DEFAULT) kctx->gamma = gamma; 3908 if (alpha != PETSC_DEFAULT) kctx->alpha = alpha; 3909 if (alpha2 != PETSC_DEFAULT) kctx->alpha2 = alpha2; 3910 if (threshold != PETSC_DEFAULT) kctx->threshold = threshold; 3911 3912 if (kctx->version < 1 || kctx->version > 3) { 3913 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Only versions 1, 2 and 3 are supported: %D",kctx->version); 3914 } 3915 if (kctx->rtol_0 < 0.0 || kctx->rtol_0 >= 1.0) { 3916 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"0.0 <= rtol_0 < 1.0: %G",kctx->rtol_0); 3917 } 3918 if (kctx->rtol_max < 0.0 || kctx->rtol_max >= 1.0) { 3919 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"0.0 <= rtol_max (%G) < 1.0\n",kctx->rtol_max); 3920 } 3921 if (kctx->gamma < 0.0 || kctx->gamma > 1.0) { 3922 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"0.0 <= gamma (%G) <= 1.0\n",kctx->gamma); 3923 } 3924 if (kctx->alpha <= 1.0 || kctx->alpha > 2.0) { 3925 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"1.0 < alpha (%G) <= 2.0\n",kctx->alpha); 3926 } 3927 if (kctx->threshold <= 0.0 || kctx->threshold >= 1.0) { 3928 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"0.0 < threshold (%G) < 1.0\n",kctx->threshold); 3929 } 3930 PetscFunctionReturn(0); 3931 } 3932 3933 #undef __FUNCT__ 3934 #define __FUNCT__ "SNESKSPGetParametersEW" 3935 /*@ 3936 SNESKSPGetParametersEW - Gets parameters for Eisenstat-Walker 3937 convergence criteria for the linear solvers within an inexact 3938 Newton method. 3939 3940 Not Collective 3941 3942 Input Parameters: 3943 snes - SNES context 3944 3945 Output Parameters: 3946 + version - version 1, 2 (default is 2) or 3 3947 . rtol_0 - initial relative tolerance (0 <= rtol_0 < 1) 3948 . rtol_max - maximum relative tolerance (0 <= rtol_max < 1) 3949 . gamma - multiplicative factor for version 2 rtol computation 3950 (0 <= gamma2 <= 1) 3951 . alpha - power for version 2 rtol computation (1 < alpha <= 2) 3952 . alpha2 - power for safeguard 3953 - threshold - threshold for imposing safeguard (0 < threshold < 1) 3954 3955 Level: advanced 3956 3957 .keywords: SNES, KSP, Eisenstat, Walker, get, parameters 3958 3959 .seealso: SNESKSPSetUseEW(), SNESKSPGetUseEW(), SNESKSPSetParametersEW() 3960 @*/ 3961 PetscErrorCode SNESKSPGetParametersEW(SNES snes,PetscInt *version,PetscReal *rtol_0,PetscReal *rtol_max, 3962 PetscReal *gamma,PetscReal *alpha,PetscReal *alpha2,PetscReal *threshold) 3963 { 3964 SNESKSPEW *kctx; 3965 PetscFunctionBegin; 3966 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3967 kctx = (SNESKSPEW*)snes->kspconvctx; 3968 if (!kctx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"No Eisenstat-Walker context existing"); 3969 if(version) *version = kctx->version; 3970 if(rtol_0) *rtol_0 = kctx->rtol_0; 3971 if(rtol_max) *rtol_max = kctx->rtol_max; 3972 if(gamma) *gamma = kctx->gamma; 3973 if(alpha) *alpha = kctx->alpha; 3974 if(alpha2) *alpha2 = kctx->alpha2; 3975 if(threshold) *threshold = kctx->threshold; 3976 PetscFunctionReturn(0); 3977 } 3978 3979 #undef __FUNCT__ 3980 #define __FUNCT__ "SNESKSPEW_PreSolve" 3981 static PetscErrorCode SNESKSPEW_PreSolve(SNES snes, KSP ksp, Vec b, Vec x) 3982 { 3983 PetscErrorCode ierr; 3984 SNESKSPEW *kctx = (SNESKSPEW*)snes->kspconvctx; 3985 PetscReal rtol=PETSC_DEFAULT,stol; 3986 3987 PetscFunctionBegin; 3988 if (!kctx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"No Eisenstat-Walker context exists"); 3989 if (!snes->iter) { /* first time in, so use the original user rtol */ 3990 rtol = kctx->rtol_0; 3991 } else { 3992 if (kctx->version == 1) { 3993 rtol = (snes->norm - kctx->lresid_last)/kctx->norm_last; 3994 if (rtol < 0.0) rtol = -rtol; 3995 stol = pow(kctx->rtol_last,kctx->alpha2); 3996 if (stol > kctx->threshold) rtol = PetscMax(rtol,stol); 3997 } else if (kctx->version == 2) { 3998 rtol = kctx->gamma * pow(snes->norm/kctx->norm_last,kctx->alpha); 3999 stol = kctx->gamma * pow(kctx->rtol_last,kctx->alpha); 4000 if (stol > kctx->threshold) rtol = PetscMax(rtol,stol); 4001 } else if (kctx->version == 3) {/* contributed by Luis Chacon, June 2006. */ 4002 rtol = kctx->gamma * pow(snes->norm/kctx->norm_last,kctx->alpha); 4003 /* safeguard: avoid sharp decrease of rtol */ 4004 stol = kctx->gamma*pow(kctx->rtol_last,kctx->alpha); 4005 stol = PetscMax(rtol,stol); 4006 rtol = PetscMin(kctx->rtol_0,stol); 4007 /* safeguard: avoid oversolving */ 4008 stol = kctx->gamma*(snes->ttol)/snes->norm; 4009 stol = PetscMax(rtol,stol); 4010 rtol = PetscMin(kctx->rtol_0,stol); 4011 } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Only versions 1, 2 or 3 are supported: %D",kctx->version); 4012 } 4013 /* safeguard: avoid rtol greater than one */ 4014 rtol = PetscMin(rtol,kctx->rtol_max); 4015 ierr = KSPSetTolerances(ksp,rtol,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);CHKERRQ(ierr); 4016 ierr = PetscInfo3(snes,"iter %D, Eisenstat-Walker (version %D) KSP rtol=%G\n",snes->iter,kctx->version,rtol);CHKERRQ(ierr); 4017 PetscFunctionReturn(0); 4018 } 4019 4020 #undef __FUNCT__ 4021 #define __FUNCT__ "SNESKSPEW_PostSolve" 4022 static PetscErrorCode SNESKSPEW_PostSolve(SNES snes, KSP ksp, Vec b, Vec x) 4023 { 4024 PetscErrorCode ierr; 4025 SNESKSPEW *kctx = (SNESKSPEW*)snes->kspconvctx; 4026 PCSide pcside; 4027 Vec lres; 4028 4029 PetscFunctionBegin; 4030 if (!kctx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"No Eisenstat-Walker context exists"); 4031 ierr = KSPGetTolerances(ksp,&kctx->rtol_last,0,0,0);CHKERRQ(ierr); 4032 ierr = SNESGetFunctionNorm(snes,&kctx->norm_last);CHKERRQ(ierr); 4033 if (kctx->version == 1) { 4034 ierr = KSPGetPCSide(ksp,&pcside);CHKERRQ(ierr); 4035 if (pcside == PC_RIGHT) { /* XXX Should we also test KSP_UNPRECONDITIONED_NORM ? */ 4036 /* KSP residual is true linear residual */ 4037 ierr = KSPGetResidualNorm(ksp,&kctx->lresid_last);CHKERRQ(ierr); 4038 } else { 4039 /* KSP residual is preconditioned residual */ 4040 /* compute true linear residual norm */ 4041 ierr = VecDuplicate(b,&lres);CHKERRQ(ierr); 4042 ierr = MatMult(snes->jacobian,x,lres);CHKERRQ(ierr); 4043 ierr = VecAYPX(lres,-1.0,b);CHKERRQ(ierr); 4044 ierr = VecNorm(lres,NORM_2,&kctx->lresid_last);CHKERRQ(ierr); 4045 ierr = VecDestroy(&lres);CHKERRQ(ierr); 4046 } 4047 } 4048 PetscFunctionReturn(0); 4049 } 4050 4051 #undef __FUNCT__ 4052 #define __FUNCT__ "SNES_KSPSolve" 4053 PetscErrorCode SNES_KSPSolve(SNES snes, KSP ksp, Vec b, Vec x) 4054 { 4055 PetscErrorCode ierr; 4056 4057 PetscFunctionBegin; 4058 if (snes->ksp_ewconv) { ierr = SNESKSPEW_PreSolve(snes,ksp,b,x);CHKERRQ(ierr); } 4059 ierr = KSPSolve(ksp,b,x);CHKERRQ(ierr); 4060 if (snes->ksp_ewconv) { ierr = SNESKSPEW_PostSolve(snes,ksp,b,x);CHKERRQ(ierr); } 4061 PetscFunctionReturn(0); 4062 } 4063 4064 #undef __FUNCT__ 4065 #define __FUNCT__ "SNESSetDM" 4066 /*@ 4067 SNESSetDM - Sets the DM that may be used by some preconditioners 4068 4069 Logically Collective on SNES 4070 4071 Input Parameters: 4072 + snes - the preconditioner context 4073 - dm - the dm 4074 4075 Level: intermediate 4076 4077 4078 .seealso: SNESGetDM(), KSPSetDM(), KSPGetDM() 4079 @*/ 4080 PetscErrorCode SNESSetDM(SNES snes,DM dm) 4081 { 4082 PetscErrorCode ierr; 4083 KSP ksp; 4084 SNESDM sdm; 4085 4086 PetscFunctionBegin; 4087 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4088 if (dm) {ierr = PetscObjectReference((PetscObject)dm);CHKERRQ(ierr);} 4089 if (snes->dm) { /* Move the SNESDM context over to the new DM unless the new DM already has one */ 4090 PetscContainer oldcontainer,container; 4091 ierr = PetscObjectQuery((PetscObject)snes->dm,"SNESDM",(PetscObject*)&oldcontainer);CHKERRQ(ierr); 4092 ierr = PetscObjectQuery((PetscObject)dm,"SNESDM",(PetscObject*)&container);CHKERRQ(ierr); 4093 if (oldcontainer && !container) { 4094 ierr = DMSNESCopyContext(snes->dm,dm);CHKERRQ(ierr); 4095 ierr = DMSNESGetContext(snes->dm,&sdm);CHKERRQ(ierr); 4096 if (sdm->originaldm == snes->dm) { /* Grant write privileges to the replacement DM */ 4097 sdm->originaldm = dm; 4098 } 4099 } 4100 ierr = DMDestroy(&snes->dm);CHKERRQ(ierr); 4101 } 4102 snes->dm = dm; 4103 ierr = SNESGetKSP(snes,&ksp);CHKERRQ(ierr); 4104 ierr = KSPSetDM(ksp,dm);CHKERRQ(ierr); 4105 ierr = KSPSetDMActive(ksp,PETSC_FALSE);CHKERRQ(ierr); 4106 if (snes->pc) { 4107 ierr = SNESSetDM(snes->pc, snes->dm);CHKERRQ(ierr); 4108 } 4109 PetscFunctionReturn(0); 4110 } 4111 4112 #undef __FUNCT__ 4113 #define __FUNCT__ "SNESGetDM" 4114 /*@ 4115 SNESGetDM - Gets the DM that may be used by some preconditioners 4116 4117 Not Collective but DM obtained is parallel on SNES 4118 4119 Input Parameter: 4120 . snes - the preconditioner context 4121 4122 Output Parameter: 4123 . dm - the dm 4124 4125 Level: intermediate 4126 4127 4128 .seealso: SNESSetDM(), KSPSetDM(), KSPGetDM() 4129 @*/ 4130 PetscErrorCode SNESGetDM(SNES snes,DM *dm) 4131 { 4132 PetscErrorCode ierr; 4133 4134 PetscFunctionBegin; 4135 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4136 if (!snes->dm) { 4137 ierr = DMShellCreate(((PetscObject)snes)->comm,&snes->dm);CHKERRQ(ierr); 4138 } 4139 *dm = snes->dm; 4140 PetscFunctionReturn(0); 4141 } 4142 4143 #undef __FUNCT__ 4144 #define __FUNCT__ "SNESSetPC" 4145 /*@ 4146 SNESSetPC - Sets the nonlinear preconditioner to be used. 4147 4148 Collective on SNES 4149 4150 Input Parameters: 4151 + snes - iterative context obtained from SNESCreate() 4152 - pc - the preconditioner object 4153 4154 Notes: 4155 Use SNESGetPC() to retrieve the preconditioner context (for example, 4156 to configure it using the API). 4157 4158 Level: developer 4159 4160 .keywords: SNES, set, precondition 4161 .seealso: SNESGetPC() 4162 @*/ 4163 PetscErrorCode SNESSetPC(SNES snes, SNES pc) 4164 { 4165 PetscErrorCode ierr; 4166 4167 PetscFunctionBegin; 4168 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); 4169 PetscValidHeaderSpecific(pc, SNES_CLASSID, 2); 4170 PetscCheckSameComm(snes, 1, pc, 2); 4171 ierr = PetscObjectReference((PetscObject) pc);CHKERRQ(ierr); 4172 ierr = SNESDestroy(&snes->pc);CHKERRQ(ierr); 4173 snes->pc = pc; 4174 ierr = PetscLogObjectParent(snes, snes->pc);CHKERRQ(ierr); 4175 PetscFunctionReturn(0); 4176 } 4177 4178 #undef __FUNCT__ 4179 #define __FUNCT__ "SNESGetPC" 4180 /*@ 4181 SNESGetPC - Returns a pointer to the nonlinear preconditioning context set with SNESSetPC(). 4182 4183 Not Collective 4184 4185 Input Parameter: 4186 . snes - iterative context obtained from SNESCreate() 4187 4188 Output Parameter: 4189 . pc - preconditioner context 4190 4191 Level: developer 4192 4193 .keywords: SNES, get, preconditioner 4194 .seealso: SNESSetPC() 4195 @*/ 4196 PetscErrorCode SNESGetPC(SNES snes, SNES *pc) 4197 { 4198 PetscErrorCode ierr; 4199 4200 PetscFunctionBegin; 4201 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); 4202 PetscValidPointer(pc, 2); 4203 if (!snes->pc) { 4204 ierr = SNESCreate(((PetscObject) snes)->comm, &snes->pc);CHKERRQ(ierr); 4205 ierr = PetscObjectIncrementTabLevel((PetscObject) snes->pc, (PetscObject) snes, 1);CHKERRQ(ierr); 4206 ierr = PetscLogObjectParent(snes, snes->pc);CHKERRQ(ierr); 4207 } 4208 *pc = snes->pc; 4209 PetscFunctionReturn(0); 4210 } 4211 4212 #undef __FUNCT__ 4213 #define __FUNCT__ "SNESSetSNESLineSearch" 4214 /*@ 4215 SNESSetSNESLineSearch - Sets the linesearch. 4216 4217 Collective on SNES 4218 4219 Input Parameters: 4220 + snes - iterative context obtained from SNESCreate() 4221 - linesearch - the linesearch object 4222 4223 Notes: 4224 Use SNESGetSNESLineSearch() to retrieve the preconditioner context (for example, 4225 to configure it using the API). 4226 4227 Level: developer 4228 4229 .keywords: SNES, set, linesearch 4230 .seealso: SNESGetSNESLineSearch() 4231 @*/ 4232 PetscErrorCode SNESSetSNESLineSearch(SNES snes, SNESLineSearch linesearch) 4233 { 4234 PetscErrorCode ierr; 4235 4236 PetscFunctionBegin; 4237 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); 4238 PetscValidHeaderSpecific(linesearch, SNESLINESEARCH_CLASSID, 2); 4239 PetscCheckSameComm(snes, 1, linesearch, 2); 4240 ierr = PetscObjectReference((PetscObject) linesearch);CHKERRQ(ierr); 4241 ierr = SNESLineSearchDestroy(&snes->linesearch);CHKERRQ(ierr); 4242 snes->linesearch = linesearch; 4243 ierr = PetscLogObjectParent(snes, snes->linesearch);CHKERRQ(ierr); 4244 PetscFunctionReturn(0); 4245 } 4246 4247 #undef __FUNCT__ 4248 #define __FUNCT__ "SNESGetSNESLineSearch" 4249 /*@C 4250 SNESGetSNESLineSearch - Returns a pointer to the line search context set with SNESSetLineSearch(). 4251 4252 Not Collective 4253 4254 Input Parameter: 4255 . snes - iterative context obtained from SNESCreate() 4256 4257 Output Parameter: 4258 . linesearch - linesearch context 4259 4260 Level: developer 4261 4262 .keywords: SNES, get, linesearch 4263 .seealso: SNESSetSNESLineSearch() 4264 @*/ 4265 PetscErrorCode SNESGetSNESLineSearch(SNES snes, SNESLineSearch *linesearch) 4266 { 4267 PetscErrorCode ierr; 4268 const char *optionsprefix; 4269 4270 PetscFunctionBegin; 4271 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); 4272 PetscValidPointer(linesearch, 2); 4273 if (!snes->linesearch) { 4274 ierr = SNESGetOptionsPrefix(snes, &optionsprefix);CHKERRQ(ierr); 4275 ierr = SNESLineSearchCreate(((PetscObject) snes)->comm, &snes->linesearch);CHKERRQ(ierr); 4276 ierr = SNESLineSearchSetSNES(snes->linesearch, snes);CHKERRQ(ierr); 4277 ierr = PetscObjectIncrementTabLevel((PetscObject) snes->linesearch, (PetscObject) snes, 1);CHKERRQ(ierr); 4278 ierr = PetscLogObjectParent(snes, snes->linesearch);CHKERRQ(ierr); 4279 } 4280 *linesearch = snes->linesearch; 4281 PetscFunctionReturn(0); 4282 } 4283 4284 #if defined(PETSC_HAVE_MATLAB_ENGINE) 4285 #include <mex.h> 4286 4287 typedef struct {char *funcname; mxArray *ctx;} SNESMatlabContext; 4288 4289 #undef __FUNCT__ 4290 #define __FUNCT__ "SNESComputeFunction_Matlab" 4291 /* 4292 SNESComputeFunction_Matlab - Calls the function that has been set with 4293 SNESSetFunctionMatlab(). 4294 4295 Collective on SNES 4296 4297 Input Parameters: 4298 + snes - the SNES context 4299 - x - input vector 4300 4301 Output Parameter: 4302 . y - function vector, as set by SNESSetFunction() 4303 4304 Notes: 4305 SNESComputeFunction() is typically used within nonlinear solvers 4306 implementations, so most users would not generally call this routine 4307 themselves. 4308 4309 Level: developer 4310 4311 .keywords: SNES, nonlinear, compute, function 4312 4313 .seealso: SNESSetFunction(), SNESGetFunction() 4314 */ 4315 PetscErrorCode SNESComputeFunction_Matlab(SNES snes,Vec x,Vec y, void *ctx) 4316 { 4317 PetscErrorCode ierr; 4318 SNESMatlabContext *sctx = (SNESMatlabContext *)ctx; 4319 int nlhs = 1,nrhs = 5; 4320 mxArray *plhs[1],*prhs[5]; 4321 long long int lx = 0,ly = 0,ls = 0; 4322 4323 PetscFunctionBegin; 4324 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4325 PetscValidHeaderSpecific(x,VEC_CLASSID,2); 4326 PetscValidHeaderSpecific(y,VEC_CLASSID,3); 4327 PetscCheckSameComm(snes,1,x,2); 4328 PetscCheckSameComm(snes,1,y,3); 4329 4330 /* call Matlab function in ctx with arguments x and y */ 4331 4332 ierr = PetscMemcpy(&ls,&snes,sizeof(snes));CHKERRQ(ierr); 4333 ierr = PetscMemcpy(&lx,&x,sizeof(x));CHKERRQ(ierr); 4334 ierr = PetscMemcpy(&ly,&y,sizeof(x));CHKERRQ(ierr); 4335 prhs[0] = mxCreateDoubleScalar((double)ls); 4336 prhs[1] = mxCreateDoubleScalar((double)lx); 4337 prhs[2] = mxCreateDoubleScalar((double)ly); 4338 prhs[3] = mxCreateString(sctx->funcname); 4339 prhs[4] = sctx->ctx; 4340 ierr = mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscSNESComputeFunctionInternal");CHKERRQ(ierr); 4341 ierr = mxGetScalar(plhs[0]);CHKERRQ(ierr); 4342 mxDestroyArray(prhs[0]); 4343 mxDestroyArray(prhs[1]); 4344 mxDestroyArray(prhs[2]); 4345 mxDestroyArray(prhs[3]); 4346 mxDestroyArray(plhs[0]); 4347 PetscFunctionReturn(0); 4348 } 4349 4350 4351 #undef __FUNCT__ 4352 #define __FUNCT__ "SNESSetFunctionMatlab" 4353 /* 4354 SNESSetFunctionMatlab - Sets the function evaluation routine and function 4355 vector for use by the SNES routines in solving systems of nonlinear 4356 equations from MATLAB. Here the function is a string containing the name of a MATLAB function 4357 4358 Logically Collective on SNES 4359 4360 Input Parameters: 4361 + snes - the SNES context 4362 . r - vector to store function value 4363 - func - function evaluation routine 4364 4365 Calling sequence of func: 4366 $ func (SNES snes,Vec x,Vec f,void *ctx); 4367 4368 4369 Notes: 4370 The Newton-like methods typically solve linear systems of the form 4371 $ f'(x) x = -f(x), 4372 where f'(x) denotes the Jacobian matrix and f(x) is the function. 4373 4374 Level: beginner 4375 4376 .keywords: SNES, nonlinear, set, function 4377 4378 .seealso: SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESSetFunction() 4379 */ 4380 PetscErrorCode SNESSetFunctionMatlab(SNES snes,Vec r,const char *func,mxArray *ctx) 4381 { 4382 PetscErrorCode ierr; 4383 SNESMatlabContext *sctx; 4384 4385 PetscFunctionBegin; 4386 /* currently sctx is memory bleed */ 4387 ierr = PetscMalloc(sizeof(SNESMatlabContext),&sctx);CHKERRQ(ierr); 4388 ierr = PetscStrallocpy(func,&sctx->funcname);CHKERRQ(ierr); 4389 /* 4390 This should work, but it doesn't 4391 sctx->ctx = ctx; 4392 mexMakeArrayPersistent(sctx->ctx); 4393 */ 4394 sctx->ctx = mxDuplicateArray(ctx); 4395 ierr = SNESSetFunction(snes,r,SNESComputeFunction_Matlab,sctx);CHKERRQ(ierr); 4396 PetscFunctionReturn(0); 4397 } 4398 4399 #undef __FUNCT__ 4400 #define __FUNCT__ "SNESComputeJacobian_Matlab" 4401 /* 4402 SNESComputeJacobian_Matlab - Calls the function that has been set with 4403 SNESSetJacobianMatlab(). 4404 4405 Collective on SNES 4406 4407 Input Parameters: 4408 + snes - the SNES context 4409 . x - input vector 4410 . A, B - the matrices 4411 - ctx - user context 4412 4413 Output Parameter: 4414 . flag - structure of the matrix 4415 4416 Level: developer 4417 4418 .keywords: SNES, nonlinear, compute, function 4419 4420 .seealso: SNESSetFunction(), SNESGetFunction() 4421 @*/ 4422 PetscErrorCode SNESComputeJacobian_Matlab(SNES snes,Vec x,Mat *A,Mat *B,MatStructure *flag, void *ctx) 4423 { 4424 PetscErrorCode ierr; 4425 SNESMatlabContext *sctx = (SNESMatlabContext *)ctx; 4426 int nlhs = 2,nrhs = 6; 4427 mxArray *plhs[2],*prhs[6]; 4428 long long int lx = 0,lA = 0,ls = 0, lB = 0; 4429 4430 PetscFunctionBegin; 4431 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4432 PetscValidHeaderSpecific(x,VEC_CLASSID,2); 4433 4434 /* call Matlab function in ctx with arguments x and y */ 4435 4436 ierr = PetscMemcpy(&ls,&snes,sizeof(snes));CHKERRQ(ierr); 4437 ierr = PetscMemcpy(&lx,&x,sizeof(x));CHKERRQ(ierr); 4438 ierr = PetscMemcpy(&lA,A,sizeof(x));CHKERRQ(ierr); 4439 ierr = PetscMemcpy(&lB,B,sizeof(x));CHKERRQ(ierr); 4440 prhs[0] = mxCreateDoubleScalar((double)ls); 4441 prhs[1] = mxCreateDoubleScalar((double)lx); 4442 prhs[2] = mxCreateDoubleScalar((double)lA); 4443 prhs[3] = mxCreateDoubleScalar((double)lB); 4444 prhs[4] = mxCreateString(sctx->funcname); 4445 prhs[5] = sctx->ctx; 4446 ierr = mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscSNESComputeJacobianInternal");CHKERRQ(ierr); 4447 ierr = mxGetScalar(plhs[0]);CHKERRQ(ierr); 4448 *flag = (MatStructure) mxGetScalar(plhs[1]);CHKERRQ(ierr); 4449 mxDestroyArray(prhs[0]); 4450 mxDestroyArray(prhs[1]); 4451 mxDestroyArray(prhs[2]); 4452 mxDestroyArray(prhs[3]); 4453 mxDestroyArray(prhs[4]); 4454 mxDestroyArray(plhs[0]); 4455 mxDestroyArray(plhs[1]); 4456 PetscFunctionReturn(0); 4457 } 4458 4459 4460 #undef __FUNCT__ 4461 #define __FUNCT__ "SNESSetJacobianMatlab" 4462 /* 4463 SNESSetJacobianMatlab - Sets the Jacobian function evaluation routine and two empty Jacobian matrices 4464 vector for use by the SNES routines in solving systems of nonlinear 4465 equations from MATLAB. Here the function is a string containing the name of a MATLAB function 4466 4467 Logically Collective on SNES 4468 4469 Input Parameters: 4470 + snes - the SNES context 4471 . A,B - Jacobian matrices 4472 . func - function evaluation routine 4473 - ctx - user context 4474 4475 Calling sequence of func: 4476 $ flag = func (SNES snes,Vec x,Mat A,Mat B,void *ctx); 4477 4478 4479 Level: developer 4480 4481 .keywords: SNES, nonlinear, set, function 4482 4483 .seealso: SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESSetFunction() 4484 */ 4485 PetscErrorCode SNESSetJacobianMatlab(SNES snes,Mat A,Mat B,const char *func,mxArray *ctx) 4486 { 4487 PetscErrorCode ierr; 4488 SNESMatlabContext *sctx; 4489 4490 PetscFunctionBegin; 4491 /* currently sctx is memory bleed */ 4492 ierr = PetscMalloc(sizeof(SNESMatlabContext),&sctx);CHKERRQ(ierr); 4493 ierr = PetscStrallocpy(func,&sctx->funcname);CHKERRQ(ierr); 4494 /* 4495 This should work, but it doesn't 4496 sctx->ctx = ctx; 4497 mexMakeArrayPersistent(sctx->ctx); 4498 */ 4499 sctx->ctx = mxDuplicateArray(ctx); 4500 ierr = SNESSetJacobian(snes,A,B,SNESComputeJacobian_Matlab,sctx);CHKERRQ(ierr); 4501 PetscFunctionReturn(0); 4502 } 4503 4504 #undef __FUNCT__ 4505 #define __FUNCT__ "SNESMonitor_Matlab" 4506 /* 4507 SNESMonitor_Matlab - Calls the function that has been set with SNESMonitorSetMatlab(). 4508 4509 Collective on SNES 4510 4511 .seealso: SNESSetFunction(), SNESGetFunction() 4512 @*/ 4513 PetscErrorCode SNESMonitor_Matlab(SNES snes,PetscInt it, PetscReal fnorm, void *ctx) 4514 { 4515 PetscErrorCode ierr; 4516 SNESMatlabContext *sctx = (SNESMatlabContext *)ctx; 4517 int nlhs = 1,nrhs = 6; 4518 mxArray *plhs[1],*prhs[6]; 4519 long long int lx = 0,ls = 0; 4520 Vec x=snes->vec_sol; 4521 4522 PetscFunctionBegin; 4523 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4524 4525 ierr = PetscMemcpy(&ls,&snes,sizeof(snes));CHKERRQ(ierr); 4526 ierr = PetscMemcpy(&lx,&x,sizeof(x));CHKERRQ(ierr); 4527 prhs[0] = mxCreateDoubleScalar((double)ls); 4528 prhs[1] = mxCreateDoubleScalar((double)it); 4529 prhs[2] = mxCreateDoubleScalar((double)fnorm); 4530 prhs[3] = mxCreateDoubleScalar((double)lx); 4531 prhs[4] = mxCreateString(sctx->funcname); 4532 prhs[5] = sctx->ctx; 4533 ierr = mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscSNESMonitorInternal");CHKERRQ(ierr); 4534 ierr = mxGetScalar(plhs[0]);CHKERRQ(ierr); 4535 mxDestroyArray(prhs[0]); 4536 mxDestroyArray(prhs[1]); 4537 mxDestroyArray(prhs[2]); 4538 mxDestroyArray(prhs[3]); 4539 mxDestroyArray(prhs[4]); 4540 mxDestroyArray(plhs[0]); 4541 PetscFunctionReturn(0); 4542 } 4543 4544 4545 #undef __FUNCT__ 4546 #define __FUNCT__ "SNESMonitorSetMatlab" 4547 /* 4548 SNESMonitorSetMatlab - Sets the monitor function from MATLAB 4549 4550 Level: developer 4551 4552 .keywords: SNES, nonlinear, set, function 4553 4554 .seealso: SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESSetFunction() 4555 */ 4556 PetscErrorCode SNESMonitorSetMatlab(SNES snes,const char *func,mxArray *ctx) 4557 { 4558 PetscErrorCode ierr; 4559 SNESMatlabContext *sctx; 4560 4561 PetscFunctionBegin; 4562 /* currently sctx is memory bleed */ 4563 ierr = PetscMalloc(sizeof(SNESMatlabContext),&sctx);CHKERRQ(ierr); 4564 ierr = PetscStrallocpy(func,&sctx->funcname);CHKERRQ(ierr); 4565 /* 4566 This should work, but it doesn't 4567 sctx->ctx = ctx; 4568 mexMakeArrayPersistent(sctx->ctx); 4569 */ 4570 sctx->ctx = mxDuplicateArray(ctx); 4571 ierr = SNESMonitorSet(snes,SNESMonitor_Matlab,sctx,PETSC_NULL);CHKERRQ(ierr); 4572 PetscFunctionReturn(0); 4573 } 4574 4575 #endif 4576