1 #include <private/snesimpl.h> /*I "petscsnes.h" I*/ 2 3 static const SNESMSType SNESMSDefault = SNESMSM62; 4 static PetscBool SNESMSRegisterAllCalled; 5 static PetscBool SNESMSPackageInitialized; 6 7 typedef struct _SNESMSTableau *SNESMSTableau; 8 struct _SNESMSTableau { 9 char *name; 10 PetscInt nstages; /* Number of stages */ 11 PetscInt nregisters; /* Number of registers */ 12 PetscReal stability; /* Scaled stability region */ 13 PetscReal *gamma; /* Coefficients of 3S* method */ 14 PetscReal *delta; /* Coefficients of 3S* method */ 15 PetscReal *betasub; /* Subdiagonal of beta in Shu-Osher form */ 16 }; 17 18 typedef struct _SNESMSTableauLink *SNESMSTableauLink; 19 struct _SNESMSTableauLink { 20 struct _SNESMSTableau tab; 21 SNESMSTableauLink next; 22 }; 23 static SNESMSTableauLink SNESMSTableauList; 24 25 typedef struct { 26 SNESMSTableau tableau; /* Tableau in low-storage form */ 27 PetscReal damping; /* Damping parameter, like length of (pseudo) time step */ 28 PetscBool norms; /* Compute norms, usually only for monitoring purposes */ 29 } SNES_MS; 30 31 #undef __FUNCT__ 32 #define __FUNCT__ "SNESMSRegisterAll" 33 /*@C 34 SNESMSRegisterAll - Registers all of the multi-stage methods in SNESMS 35 36 Not Collective, but should be called by all processes which will need the schemes to be registered 37 38 Level: advanced 39 40 .keywords: SNES, SNESMS, register, all 41 42 .seealso: SNESMSRegisterDestroy() 43 @*/ 44 PetscErrorCode SNESMSRegisterAll(void) 45 { 46 PetscErrorCode ierr; 47 48 PetscFunctionBegin; 49 if (SNESMSRegisterAllCalled) PetscFunctionReturn(0); 50 SNESMSRegisterAllCalled = PETSC_TRUE; 51 52 { 53 PetscReal gamma[3][1] = {{1.0},{0.0},{0.0}}; 54 PetscReal delta[1] = {0.0}; 55 PetscReal betasub[1] = {1.0}; 56 ierr = SNESMSRegister(SNESMSEULER,1,3,1.0,&gamma[0][0],delta,betasub);CHKERRQ(ierr); 57 } 58 59 { 60 PetscReal gamma[3][6] = { 61 {0.0000000000000000E+00, -7.0304722367110606E-01, -1.9836719667506464E-01, -1.6023843981863788E+00, 9.4483822882855284E-02, -1.4204296130641869E-01}, 62 {1.0000000000000000E+00, 1.1111025767083920E+00, 5.6150921583923230E-01, 7.4151723494934041E-01, 3.1714538168600587E-01, 4.6479276238548706E-01}, 63 {0.0000000000000000E+00, 0.0000000000000000E+00, 0.0000000000000000E+00, 6.7968174970583317E-01, -4.1755042846051737E-03, -1.9115668129923846E-01} 64 }; 65 PetscReal delta[6] = {1.0000000000000000E+00, 5.3275427433201750E-01, 6.0143544663985238E-01, 4.5874077053842177E-01, 2.7544386906104651E-01, 0.0000000000000000E+00}; 66 PetscReal betasub[6] = {8.4753115429481929E-01, 7.4018896368655618E-01, 6.5963574086583309E-03, 4.6747795645517759E-01, 1.3314545813643919E-01, 5.3260800028018784E-01}; 67 ierr = SNESMSRegister(SNESMSM62,6,3,1.0,&gamma[0][0],delta,betasub);CHKERRQ(ierr); 68 } 69 PetscFunctionReturn(0); 70 } 71 72 #undef __FUNCT__ 73 #define __FUNCT__ "SNESMSRegisterDestroy" 74 /*@C 75 SNESMSRegisterDestroy - Frees the list of schemes that were registered by TSRosWRegister(). 76 77 Not Collective 78 79 Level: advanced 80 81 .keywords: TSRosW, register, destroy 82 .seealso: TSRosWRegister(), TSRosWRegisterAll(), TSRosWRegisterDynamic() 83 @*/ 84 PetscErrorCode SNESMSRegisterDestroy(void) 85 { 86 PetscErrorCode ierr; 87 SNESMSTableauLink link; 88 89 PetscFunctionBegin; 90 while ((link = SNESMSTableauList)) { 91 SNESMSTableau t = &link->tab; 92 SNESMSTableauList = link->next; 93 ierr = PetscFree3(t->gamma,t->delta,t->betasub);CHKERRQ(ierr); 94 ierr = PetscFree(t->name);CHKERRQ(ierr); 95 ierr = PetscFree(link);CHKERRQ(ierr); 96 } 97 SNESMSRegisterAllCalled = PETSC_FALSE; 98 PetscFunctionReturn(0); 99 } 100 101 #undef __FUNCT__ 102 #define __FUNCT__ "SNESMSInitializePackage" 103 /*@C 104 SNESMSInitializePackage - This function initializes everything in the SNESMS package. It is called 105 from PetscDLLibraryRegister() when using dynamic libraries, and on the first call to SNESCreate_MS() 106 when using static libraries. 107 108 Input Parameter: 109 path - The dynamic library path, or PETSC_NULL 110 111 Level: developer 112 113 .keywords: SNES, SNESMS, initialize, package 114 .seealso: PetscInitialize() 115 @*/ 116 PetscErrorCode SNESMSInitializePackage(const char path[]) 117 { 118 PetscErrorCode ierr; 119 120 PetscFunctionBegin; 121 if (SNESMSPackageInitialized) PetscFunctionReturn(0); 122 SNESMSPackageInitialized = PETSC_TRUE; 123 ierr = SNESMSRegisterAll();CHKERRQ(ierr); 124 ierr = PetscRegisterFinalize(SNESMSFinalizePackage);CHKERRQ(ierr); 125 PetscFunctionReturn(0); 126 } 127 128 #undef __FUNCT__ 129 #define __FUNCT__ "SNESMSFinalizePackage" 130 /*@C 131 SNESMSFinalizePackage - This function destroys everything in the SNESMS package. It is 132 called from PetscFinalize(). 133 134 Level: developer 135 136 .keywords: Petsc, destroy, package 137 .seealso: PetscFinalize() 138 @*/ 139 PetscErrorCode SNESMSFinalizePackage(void) 140 { 141 PetscErrorCode ierr; 142 143 PetscFunctionBegin; 144 SNESMSPackageInitialized = PETSC_FALSE; 145 ierr = SNESMSRegisterDestroy();CHKERRQ(ierr); 146 PetscFunctionReturn(0); 147 } 148 149 #undef __FUNCT__ 150 #define __FUNCT__ "SNESMSRegister" 151 /*@C 152 SNESMSRegister - register a multistage scheme 153 154 Not Collective, but the same schemes should be registered on all processes on which they will be used 155 156 Input Parameters: 157 + name - identifier for method 158 . nstages - number of stages 159 . nregisters - number of registers used by low-storage implementation 160 . gamma - coefficients, see Ketcheson's paper 161 . delta - coefficients, see Ketcheson's paper 162 - betasub - subdiagonal of Shu-Osher form 163 164 Notes: 165 The notation is described in Ketcheson (2010) Runge-Kutta methods with minimum storage implementations. 166 167 Level: advanced 168 169 .keywords: SNES, register 170 171 .seealso: SNESMS 172 @*/ 173 PetscErrorCode SNESMSRegister(const SNESMSType name,PetscInt nstages,PetscInt nregisters,PetscReal stability,const PetscReal gamma[],const PetscReal delta[],const PetscReal betasub[]) 174 { 175 PetscErrorCode ierr; 176 SNESMSTableauLink link; 177 SNESMSTableau t; 178 179 PetscFunctionBegin; 180 PetscValidCharPointer(name,1); 181 if (nstages < 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Must have at least one stage"); 182 if (nregisters != 3) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only support for methods written in 3-register form"); 183 PetscValidPointer(gamma,4); 184 PetscValidPointer(delta,5); 185 PetscValidPointer(betasub,6); 186 187 ierr = PetscMalloc(sizeof(*link),&link);CHKERRQ(ierr); 188 ierr = PetscMemzero(link,sizeof(*link));CHKERRQ(ierr); 189 t = &link->tab; 190 ierr = PetscStrallocpy(name,&t->name);CHKERRQ(ierr); 191 t->nstages = nstages; 192 t->nregisters = nregisters; 193 t->stability = stability; 194 ierr = PetscMalloc3(nstages*nregisters,PetscReal,&t->gamma,nstages,PetscReal,&t->delta,nstages,PetscReal,&t->betasub);CHKERRQ(ierr); 195 ierr = PetscMemcpy(t->gamma,gamma,nstages*nregisters*sizeof(PetscReal));CHKERRQ(ierr); 196 ierr = PetscMemcpy(t->delta,delta,nstages*sizeof(PetscReal));CHKERRQ(ierr); 197 ierr = PetscMemcpy(t->betasub,betasub,nstages*sizeof(PetscReal));CHKERRQ(ierr); 198 199 link->next = SNESMSTableauList; 200 SNESMSTableauList = link; 201 PetscFunctionReturn(0); 202 } 203 204 #undef __FUNCT__ 205 #define __FUNCT__ "SNESMSStep_3Sstar" 206 /* 207 X - initial state, updated in-place. 208 F - residual, computed at the initial X on input 209 */ 210 static PetscErrorCode SNESMSStep_3Sstar(SNES snes,Vec X,Vec F) 211 { 212 PetscErrorCode ierr; 213 SNES_MS *ms = (SNES_MS*)snes->data; 214 SNESMSTableau t = ms->tableau; 215 const PetscReal *gamma = t->gamma,*delta = t->delta,*betasub = t->betasub; 216 Vec S1,S2,S3,Y; 217 PetscInt i,nstages = t->nstages;; 218 219 220 PetscFunctionBegin; 221 Y = snes->work[0]; 222 S1 = X; 223 S2 = snes->work[1]; 224 S3 = snes->work[2]; 225 ierr = VecZeroEntries(S2);CHKERRQ(ierr); 226 ierr = VecCopy(X,S3);CHKERRQ(ierr); 227 for (i=0; i<nstages; i++) { 228 Vec Ss[4] = {S1,S2,S3,Y}; 229 PetscReal scoeff[4] = {gamma[0*nstages+i]-1.0,gamma[1*nstages+i],gamma[2*nstages+i],-betasub[i]*ms->damping}; 230 ierr = VecAXPY(S2,delta[i],S1);CHKERRQ(ierr); 231 if (i>0) { 232 ierr = SNESComputeFunction(snes,S1,F);CHKERRQ(ierr); 233 if (snes->domainerror) { 234 snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; 235 PetscFunctionReturn(0); 236 } 237 } 238 ierr = SNES_KSPSolve(snes,snes->ksp,F,Y);CHKERRQ(ierr); 239 ierr = VecMAXPY(S1,4,scoeff,Ss);CHKERRQ(ierr); 240 } 241 PetscFunctionReturn(0); 242 } 243 244 #undef __FUNCT__ 245 #define __FUNCT__ "SNESSolve_MS" 246 static PetscErrorCode SNESSolve_MS(SNES snes) 247 { 248 SNES_MS *ms = (SNES_MS*)snes->data; 249 Vec X = snes->vec_sol,F = snes->vec_func; 250 PetscReal fnorm; 251 MatStructure mstruct; 252 PetscInt i; 253 PetscErrorCode ierr; 254 255 PetscFunctionBegin; 256 snes->reason = SNES_CONVERGED_ITERATING; 257 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 258 snes->iter = 0; 259 snes->norm = 0.; 260 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 261 262 ierr = SNESComputeFunction(snes,X,F);CHKERRQ(ierr); 263 if (snes->domainerror) { 264 snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; 265 PetscFunctionReturn(0); 266 } 267 if (snes->jacobian) { /* This method does not require a Jacobian, but it is usually preconditioned by PBJacobi */ 268 ierr = SNESComputeJacobian(snes,snes->vec_sol,&snes->jacobian,&snes->jacobian_pre,&mstruct);CHKERRQ(ierr); 269 ierr = KSPSetOperators(snes->ksp,snes->jacobian,snes->jacobian_pre,mstruct);CHKERRQ(ierr); 270 } 271 if (ms->norms) { 272 ierr = VecNorm(F,NORM_2,&fnorm);CHKERRQ(ierr); /* fnorm <- ||F|| */ 273 if (PetscIsInfOrNanReal(fnorm)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"Infinite or not-a-number generated in norm"); 274 /* Monitor convergence */ 275 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 276 snes->iter = 0; 277 snes->norm = fnorm; 278 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 279 SNESLogConvHistory(snes,snes->norm,0); 280 ierr = SNESMonitor(snes,snes->iter,snes->norm);CHKERRQ(ierr); 281 282 /* set parameter for default relative tolerance convergence test */ 283 snes->ttol = fnorm*snes->rtol; 284 /* Test for convergence */ 285 ierr = (*snes->ops->converged)(snes,snes->iter,0.0,0.0,fnorm,&snes->reason,snes->cnvP);CHKERRQ(ierr); 286 if (snes->reason) PetscFunctionReturn(0); 287 } 288 289 /* Call general purpose update function */ 290 if (snes->ops->update) { 291 ierr = (*snes->ops->update)(snes,snes->iter);CHKERRQ(ierr); 292 } 293 for (i = 0; i < snes->max_its; i++) { 294 ierr = SNESMSStep_3Sstar(snes,X,F);CHKERRQ(ierr); 295 296 if (i+1 < snes->max_its || ms->norms) { 297 ierr = SNESComputeFunction(snes,X,F);CHKERRQ(ierr); 298 if (snes->domainerror) { 299 snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; 300 PetscFunctionReturn(0); 301 } 302 } 303 304 if (ms->norms) { 305 ierr = VecNorm(F,NORM_2,&fnorm);CHKERRQ(ierr); /* fnorm <- ||F|| */ 306 if (PetscIsInfOrNanReal(fnorm)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"Infinite or not-a-number generated in norm"); 307 /* Monitor convergence */ 308 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 309 snes->iter = i+1; 310 snes->norm = fnorm; 311 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 312 SNESLogConvHistory(snes,snes->norm,0); 313 ierr = SNESMonitor(snes,snes->iter,snes->norm);CHKERRQ(ierr); 314 315 /* Test for convergence */ 316 ierr = (*snes->ops->converged)(snes,snes->iter,0.0,0.0,fnorm,&snes->reason,snes->cnvP);CHKERRQ(ierr); 317 if (snes->reason) PetscFunctionReturn(0); 318 } 319 320 /* Call general purpose update function */ 321 if (snes->ops->update) { 322 ierr = (*snes->ops->update)(snes, snes->iter);CHKERRQ(ierr); 323 } 324 } 325 if (!snes->reason) snes->reason = SNES_CONVERGED_ITS; 326 PetscFunctionReturn(0); 327 } 328 329 #undef __FUNCT__ 330 #define __FUNCT__ "SNESSetUp_MS" 331 static PetscErrorCode SNESSetUp_MS(SNES snes) 332 { 333 SNES_MS * ms = (SNES_MS *)snes->data; 334 PetscErrorCode ierr; 335 336 PetscFunctionBegin; 337 if (!ms->tableau) {ierr = SNESMSSetType(snes,SNESMSDefault);CHKERRQ(ierr);} 338 ierr = SNESDefaultGetWork(snes,3);CHKERRQ(ierr); 339 PetscFunctionReturn(0); 340 } 341 342 #undef __FUNCT__ 343 #define __FUNCT__ "SNESReset_MS" 344 static PetscErrorCode SNESReset_MS(SNES snes) 345 { 346 347 PetscFunctionBegin; 348 PetscFunctionReturn(0); 349 } 350 351 #undef __FUNCT__ 352 #define __FUNCT__ "SNESDestroy_MS" 353 static PetscErrorCode SNESDestroy_MS(SNES snes) 354 { 355 PetscErrorCode ierr; 356 357 PetscFunctionBegin; 358 ierr = PetscFree(snes->data);CHKERRQ(ierr); 359 ierr = PetscObjectComposeFunctionDynamic((PetscObject)snes,"","",PETSC_NULL);CHKERRQ(ierr); 360 PetscFunctionReturn(0); 361 } 362 363 #undef __FUNCT__ 364 #define __FUNCT__ "SNESView_MS" 365 static PetscErrorCode SNESView_MS(SNES snes,PetscViewer viewer) 366 { 367 PetscBool iascii; 368 PetscErrorCode ierr; 369 SNES_MS *ms = (SNES_MS*)snes->data; 370 371 PetscFunctionBegin; 372 ierr = PetscTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);CHKERRQ(ierr); 373 if (iascii) { 374 SNESMSTableau tab = ms->tableau; 375 ierr = PetscViewerASCIIPrintf(viewer," multi-stage method type: %s\n",tab?tab->name:"not yet set");CHKERRQ(ierr); 376 } 377 PetscFunctionReturn(0); 378 } 379 380 #undef __FUNCT__ 381 #define __FUNCT__ "SNESSetFromOptions_MS" 382 static PetscErrorCode SNESSetFromOptions_MS(SNES snes) 383 { 384 SNES_MS *ms = (SNES_MS*)snes->data;; 385 PetscErrorCode ierr; 386 387 PetscFunctionBegin; 388 ierr = PetscOptionsHead("SNES MS options");CHKERRQ(ierr); 389 { 390 SNESMSTableauLink link; 391 PetscInt count,choice; 392 PetscBool flg; 393 const char **namelist; 394 char mstype[256]; 395 396 ierr = PetscStrncpy(mstype,SNESMSDefault,sizeof mstype);CHKERRQ(ierr); 397 for (link=SNESMSTableauList,count=0; link; link=link->next,count++) ; 398 ierr = PetscMalloc(count*sizeof(char*),&namelist);CHKERRQ(ierr); 399 for (link=SNESMSTableauList,count=0; link; link=link->next,count++) namelist[count] = link->tab.name; 400 ierr = PetscOptionsEList("-snes_ms_type","Multistage smoother type","SNESMSSetType",(const char*const*)namelist,count,mstype,&choice,&flg);CHKERRQ(ierr); 401 ierr = SNESMSSetType(snes,flg ? namelist[choice] : mstype);CHKERRQ(ierr); 402 ierr = PetscFree(namelist);CHKERRQ(ierr); 403 ierr = PetscOptionsReal("-snes_ms_damping","Damping for multistage method","SNESMSSetDamping",ms->damping,&ms->damping,PETSC_NULL);CHKERRQ(ierr); 404 ierr = PetscOptionsBool("-snes_ms_norms","Compute norms for monitoring","none",ms->norms,&ms->norms,PETSC_NULL);CHKERRQ(ierr); 405 } 406 ierr = PetscOptionsTail();CHKERRQ(ierr); 407 PetscFunctionReturn(0); 408 } 409 410 EXTERN_C_BEGIN 411 #undef __FUNCT__ 412 #define __FUNCT__ "SNESMSSetType_MS" 413 PetscErrorCode SNESMSSetType_MS(SNES snes,SNESMSType mstype) 414 { 415 PetscErrorCode ierr; 416 SNES_MS *ms = (SNES_MS*)snes->data; 417 SNESMSTableauLink link; 418 PetscBool match; 419 420 PetscFunctionBegin; 421 if (ms->tableau) { 422 ierr = PetscStrcmp(ms->tableau->name,mstype,&match);CHKERRQ(ierr); 423 if (match) PetscFunctionReturn(0); 424 } 425 for (link = SNESMSTableauList; link; link=link->next) { 426 ierr = PetscStrcmp(link->tab.name,mstype,&match);CHKERRQ(ierr); 427 if (match) { 428 ierr = SNESReset_MS(snes);CHKERRQ(ierr); 429 ms->tableau = &link->tab; 430 PetscFunctionReturn(0); 431 } 432 } 433 SETERRQ1(((PetscObject)snes)->comm,PETSC_ERR_ARG_UNKNOWN_TYPE,"Could not find '%s'",mstype); 434 PetscFunctionReturn(0); 435 } 436 EXTERN_C_END 437 438 439 #undef __FUNCT__ 440 #define __FUNCT__ "SNESMSSetType" 441 /*@C 442 SNESMSSetType - Set the type of multistage smoother 443 444 Logically collective 445 446 Input Parameter: 447 + snes - nonlinear solver context 448 - mstype - type of multistage method 449 450 Level: beginner 451 452 .seealso: SNESMSGetType(), SNESMS 453 @*/ 454 PetscErrorCode SNESMSSetType(SNES snes,const SNESMSType rostype) 455 { 456 PetscErrorCode ierr; 457 458 PetscFunctionBegin; 459 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 460 ierr = PetscTryMethod(snes,"SNESMSSetType_C",(SNES,const SNESMSType),(snes,rostype));CHKERRQ(ierr); 461 PetscFunctionReturn(0); 462 } 463 464 /* -------------------------------------------------------------------------- */ 465 /*MC 466 SNESMS - multi-stage smoothers 467 468 Options Database: 469 470 + -snes_ms_type - type of multi-stage smoother 471 - -snes_ms_damping - damping for multi-stage method 472 473 Notes: 474 These multistage methods are explicit Runge-Kutta methods that are often used as smoothers for 475 FAS multigrid for transport problems. In the linear case, these are equivalent to polynomial smoothers (such as Chebychev). 476 477 Multi-stage smoothers should usually be preconditioned by point-block Jacobi to ensure proper scaling and to normalize the wave speeds. 478 479 The methods are specified in low storage form (Ketcheson 2010). New methods can be registered with SNESMSRegister(). 480 481 References: 482 483 Ketcheson (2010) Runge-Kutta methods with minimum storage implementations. 484 485 Jameson (1983) Solution of the Euler equations for two dimensional transonic flow by a multigrid method. 486 487 Pierce and Giles (1997) Preconditioned multigrid methods for compressible flow calculations on stretched meshes. 488 489 Level: beginner 490 491 .seealso: SNESCreate(), SNES, SNESSetType(), SNESMS, SNESFAS, KSPCHEBYCHEV 492 493 M*/ 494 EXTERN_C_BEGIN 495 #undef __FUNCT__ 496 #define __FUNCT__ "SNESCreate_MS" 497 PetscErrorCode SNESCreate_MS(SNES snes) 498 { 499 PetscErrorCode ierr; 500 SNES_MS *ms; 501 502 PetscFunctionBegin; 503 #if !defined(PETSC_USE_DYNAMIC_LIBRARIES) 504 ierr = SNESMSInitializePackage(PETSC_NULL);CHKERRQ(ierr); 505 #endif 506 507 snes->ops->setup = SNESSetUp_MS; 508 snes->ops->solve = SNESSolve_MS; 509 snes->ops->destroy = SNESDestroy_MS; 510 snes->ops->setfromoptions = SNESSetFromOptions_MS; 511 snes->ops->view = SNESView_MS; 512 snes->ops->reset = SNESReset_MS; 513 514 snes->usespc = PETSC_FALSE; 515 snes->usesksp = PETSC_TRUE; 516 517 ierr = PetscNewLog(snes,SNES_MS,&ms);CHKERRQ(ierr); 518 snes->data = (void*)ms; 519 ms->damping = 0.9; 520 ms->norms = PETSC_FALSE; 521 522 ierr = PetscObjectComposeFunctionDynamic((PetscObject)snes,"SNESMSSetType_C","SNESMSSetType_MS",SNESMSSetType_MS);CHKERRQ(ierr); 523 PetscFunctionReturn(0); 524 } 525 EXTERN_C_END 526