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