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