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