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