Actual source code: ssp.c
petsc-3.5.4 2015-05-23
1: /*
2: Code for Timestepping with explicit SSP.
3: */
4: #include <petsc-private/tsimpl.h> /*I "petscts.h" I*/
6: PetscFunctionList TSSSPList = 0;
7: static PetscBool TSSSPPackageInitialized;
9: typedef struct {
10: PetscErrorCode (*onestep)(TS,PetscReal,PetscReal,Vec);
11: char *type_name;
12: PetscInt nstages;
13: Vec *work;
14: PetscInt nwork;
15: PetscBool workout;
16: } TS_SSP;
21: static PetscErrorCode TSSSPGetWorkVectors(TS ts,PetscInt n,Vec **work)
22: {
23: TS_SSP *ssp = (TS_SSP*)ts->data;
27: if (ssp->workout) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Work vectors already gotten");
28: if (ssp->nwork < n) {
29: if (ssp->nwork > 0) {
30: VecDestroyVecs(ssp->nwork,&ssp->work);
31: }
32: VecDuplicateVecs(ts->vec_sol,n,&ssp->work);
33: ssp->nwork = n;
34: }
35: *work = ssp->work;
36: ssp->workout = PETSC_TRUE;
37: return(0);
38: }
42: static PetscErrorCode TSSSPRestoreWorkVectors(TS ts,PetscInt n,Vec **work)
43: {
44: TS_SSP *ssp = (TS_SSP*)ts->data;
47: if (!ssp->workout) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Work vectors have not been gotten");
48: if (*work != ssp->work) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Wrong work vectors checked out");
49: ssp->workout = PETSC_FALSE;
50: *work = NULL;
51: return(0);
52: }
56: /*MC
57: TSSSPRKS2 - Optimal second order SSP Runge-Kutta method, low-storage, c_eff=(s-1)/s
59: Pseudocode 2 of Ketcheson 2008
61: Level: beginner
63: .seealso: TSSSP, TSSSPSetType(), TSSSPSetNumStages()
64: M*/
65: static PetscErrorCode TSSSPStep_RK_2(TS ts,PetscReal t0,PetscReal dt,Vec sol)
66: {
67: TS_SSP *ssp = (TS_SSP*)ts->data;
68: Vec *work,F;
69: PetscInt i,s;
73: s = ssp->nstages;
74: TSSSPGetWorkVectors(ts,2,&work);
75: F = work[1];
76: VecCopy(sol,work[0]);
77: for (i=0; i<s-1; i++) {
78: PetscReal stage_time = t0+dt*(i/(s-1.));
79: TSPreStage(ts,stage_time);
80: TSComputeRHSFunction(ts,stage_time,work[0],F);
81: VecAXPY(work[0],dt/(s-1.),F);
82: }
83: TSComputeRHSFunction(ts,t0+dt,work[0],F);
84: VecAXPBYPCZ(sol,(s-1.)/s,dt/s,1./s,work[0],F);
85: TSSSPRestoreWorkVectors(ts,2,&work);
86: return(0);
87: }
91: /*MC
92: TSSSPRKS3 - Optimal third order SSP Runge-Kutta, low-storage, c_eff=(PetscSqrtReal(s)-1)/PetscSqrtReal(s), where PetscSqrtReal(s) is an integer
94: Pseudocode 2 of Ketcheson 2008
96: Level: beginner
98: .seealso: TSSSP, TSSSPSetType(), TSSSPSetNumStages()
99: M*/
100: static PetscErrorCode TSSSPStep_RK_3(TS ts,PetscReal t0,PetscReal dt,Vec sol)
101: {
102: TS_SSP *ssp = (TS_SSP*)ts->data;
103: Vec *work,F;
104: PetscInt i,s,n,r;
105: PetscReal c,stage_time;
109: s = ssp->nstages;
110: n = (PetscInt)(PetscSqrtReal((PetscReal)s)+0.001);
111: r = s-n;
112: if (n*n != s) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support for optimal third order schemes with %d stages, must be a square number at least 4",s);
113: TSSSPGetWorkVectors(ts,3,&work);
114: F = work[2];
115: VecCopy(sol,work[0]);
116: for (i=0; i<(n-1)*(n-2)/2; i++) {
117: c = (i<n*(n+1)/2) ? 1.*i/(s-n) : (1.*i-n)/(s-n);
118: stage_time = t0+c*dt;
119: TSPreStage(ts,stage_time);
120: TSComputeRHSFunction(ts,stage_time,work[0],F);
121: VecAXPY(work[0],dt/r,F);
122: }
123: VecCopy(work[0],work[1]);
124: for (; i<n*(n+1)/2-1; i++) {
125: c = (i<n*(n+1)/2) ? 1.*i/(s-n) : (1.*i-n)/(s-n);
126: stage_time = t0+c*dt;
127: TSPreStage(ts,stage_time);
128: TSComputeRHSFunction(ts,stage_time,work[0],F);
129: VecAXPY(work[0],dt/r,F);
130: }
131: {
132: c = (i<n*(n+1)/2) ? 1.*i/(s-n) : (1.*i-n)/(s-n);
133: stage_time = t0+c*dt;
134: TSPreStage(ts,stage_time);
135: TSComputeRHSFunction(ts,stage_time,work[0],F);
136: VecAXPBYPCZ(work[0],1.*n/(2*n-1.),(n-1.)*dt/(r*(2*n-1)),(n-1.)/(2*n-1.),work[1],F);
137: i++;
138: }
139: for (; i<s; i++) {
140: c = (i<n*(n+1)/2) ? 1.*i/(s-n) : (1.*i-n)/(s-n);
141: stage_time = t0+c*dt;
142: TSPreStage(ts,stage_time);
143: TSComputeRHSFunction(ts,stage_time,work[0],F);
144: VecAXPY(work[0],dt/r,F);
145: }
146: VecCopy(work[0],sol);
147: TSSSPRestoreWorkVectors(ts,3,&work);
148: return(0);
149: }
153: /*MC
154: TSSSPRKS104 - Optimal fourth order SSP Runge-Kutta, low-storage (2N), c_eff=0.6
156: SSPRK(10,4), Pseudocode 3 of Ketcheson 2008
158: Level: beginner
160: .seealso: TSSSP, TSSSPSetType()
161: M*/
162: static PetscErrorCode TSSSPStep_RK_10_4(TS ts,PetscReal t0,PetscReal dt,Vec sol)
163: {
164: const PetscReal c[10] = {0, 1./6, 2./6, 3./6, 4./6, 2./6, 3./6, 4./6, 5./6, 1};
165: Vec *work,F;
166: PetscInt i;
167: PetscReal stage_time;
168: PetscErrorCode ierr;
171: TSSSPGetWorkVectors(ts,3,&work);
172: F = work[2];
173: VecCopy(sol,work[0]);
174: for (i=0; i<5; i++) {
175: stage_time = t0+c[i]*dt;
176: TSPreStage(ts,stage_time);
177: TSComputeRHSFunction(ts,stage_time,work[0],F);
178: VecAXPY(work[0],dt/6,F);
179: }
180: VecAXPBYPCZ(work[1],1./25,9./25,0,sol,work[0]);
181: VecAXPBY(work[0],15,-5,work[1]);
182: for (; i<9; i++) {
183: stage_time = t0+c[i]*dt;
184: TSPreStage(ts,stage_time);
185: TSComputeRHSFunction(ts,stage_time,work[0],F);
186: VecAXPY(work[0],dt/6,F);
187: }
188: stage_time = t0+dt;
189: TSPreStage(ts,stage_time);
190: TSComputeRHSFunction(ts,stage_time,work[0],F);
191: VecAXPBYPCZ(work[1],3./5,dt/10,1,work[0],F);
192: VecCopy(work[1],sol);
193: TSSSPRestoreWorkVectors(ts,3,&work);
194: return(0);
195: }
200: static PetscErrorCode TSSetUp_SSP(TS ts)
201: {
204: return(0);
205: }
209: static PetscErrorCode TSStep_SSP(TS ts)
210: {
211: TS_SSP *ssp = (TS_SSP*)ts->data;
212: Vec sol = ts->vec_sol;
216: TSPreStep(ts);
217: (*ssp->onestep)(ts,ts->ptime,ts->time_step,sol);
218: ts->ptime += ts->time_step;
219: ts->steps++;
220: return(0);
221: }
222: /*------------------------------------------------------------*/
225: static PetscErrorCode TSReset_SSP(TS ts)
226: {
227: TS_SSP *ssp = (TS_SSP*)ts->data;
231: if (ssp->work) {VecDestroyVecs(ssp->nwork,&ssp->work);}
232: ssp->nwork = 0;
233: ssp->workout = PETSC_FALSE;
234: return(0);
235: }
239: static PetscErrorCode TSDestroy_SSP(TS ts)
240: {
241: TS_SSP *ssp = (TS_SSP*)ts->data;
245: TSReset_SSP(ts);
246: PetscFree(ssp->type_name);
247: PetscFree(ts->data);
248: PetscObjectComposeFunction((PetscObject)ts,"TSSSPGetType_C",NULL);
249: PetscObjectComposeFunction((PetscObject)ts,"TSSSPSetType_C",NULL);
250: PetscObjectComposeFunction((PetscObject)ts,"TSSSPGetNumStages_C",NULL);
251: PetscObjectComposeFunction((PetscObject)ts,"TSSSPSetNumStages_C",NULL);
252: return(0);
253: }
254: /*------------------------------------------------------------*/
258: /*@C
259: TSSSPSetType - set the SSP time integration scheme to use
261: Logically Collective
263: Input Arguments:
264: ts - time stepping object
265: type - type of scheme to use
267: Options Database Keys:
268: -ts_ssp_type <rks2>: Type of SSP method (one of) rks2 rks3 rk104
269: -ts_ssp_nstages <5>: Number of stages
271: Level: beginner
273: .seealso: TSSSP, TSSSPGetType(), TSSSPSetNumStages(), TSSSPRKS2, TSSSPRKS3, TSSSPRK104
274: @*/
275: PetscErrorCode TSSSPSetType(TS ts,TSSSPType type)
276: {
281: PetscTryMethod(ts,"TSSSPSetType_C",(TS,TSSSPType),(ts,type));
282: return(0);
283: }
287: /*@C
288: TSSSPGetType - get the SSP time integration scheme
290: Logically Collective
292: Input Argument:
293: ts - time stepping object
295: Output Argument:
296: type - type of scheme being used
298: Level: beginner
300: .seealso: TSSSP, TSSSPSettype(), TSSSPSetNumStages(), TSSSPRKS2, TSSSPRKS3, TSSSPRK104
301: @*/
302: PetscErrorCode TSSSPGetType(TS ts,TSSSPType *type)
303: {
308: PetscTryMethod(ts,"TSSSPGetType_C",(TS,TSSSPType*),(ts,type));
309: return(0);
310: }
314: /*@
315: TSSSPSetNumStages - set the number of stages to use with the SSP method
317: Logically Collective
319: Input Arguments:
320: ts - time stepping object
321: nstages - number of stages
323: Options Database Keys:
324: -ts_ssp_type <rks2>: NumStages of SSP method (one of) rks2 rks3 rk104
325: -ts_ssp_nstages <5>: Number of stages
327: Level: beginner
329: .seealso: TSSSP, TSSSPGetNumStages(), TSSSPSetNumStages(), TSSSPRKS2, TSSSPRKS3, TSSSPRK104
330: @*/
331: PetscErrorCode TSSSPSetNumStages(TS ts,PetscInt nstages)
332: {
337: PetscTryMethod(ts,"TSSSPSetNumStages_C",(TS,PetscInt),(ts,nstages));
338: return(0);
339: }
343: /*@
344: TSSSPGetNumStages - get the number of stages in the SSP time integration scheme
346: Logically Collective
348: Input Argument:
349: ts - time stepping object
351: Output Argument:
352: nstages - number of stages
354: Level: beginner
356: .seealso: TSSSP, TSSSPGetType(), TSSSPSetNumStages(), TSSSPRKS2, TSSSPRKS3, TSSSPRK104
357: @*/
358: PetscErrorCode TSSSPGetNumStages(TS ts,PetscInt *nstages)
359: {
364: PetscTryMethod(ts,"TSSSPGetNumStages_C",(TS,PetscInt*),(ts,nstages));
365: return(0);
366: }
370: PETSC_EXTERN PetscErrorCode TSSSPSetType_SSP(TS ts,TSSSPType type)
371: {
372: PetscErrorCode ierr,(*r)(TS,PetscReal,PetscReal,Vec);
373: TS_SSP *ssp = (TS_SSP*)ts->data;
376: PetscFunctionListFind(TSSSPList,type,&r);
377: if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown TS_SSP type %s given",type);
378: ssp->onestep = r;
379: PetscFree(ssp->type_name);
380: PetscStrallocpy(type,&ssp->type_name);
381: return(0);
382: }
385: PetscErrorCode TSSSPGetType_SSP(TS ts,TSSSPType *type)
386: {
387: TS_SSP *ssp = (TS_SSP*)ts->data;
390: *type = ssp->type_name;
391: return(0);
392: }
395: PetscErrorCode TSSSPSetNumStages_SSP(TS ts,PetscInt nstages)
396: {
397: TS_SSP *ssp = (TS_SSP*)ts->data;
400: ssp->nstages = nstages;
401: return(0);
402: }
405: PetscErrorCode TSSSPGetNumStages_SSP(TS ts,PetscInt *nstages)
406: {
407: TS_SSP *ssp = (TS_SSP*)ts->data;
410: *nstages = ssp->nstages;
411: return(0);
412: }
416: static PetscErrorCode TSSetFromOptions_SSP(TS ts)
417: {
418: char tname[256] = TSSSPRKS2;
419: TS_SSP *ssp = (TS_SSP*)ts->data;
421: PetscBool flg;
424: PetscOptionsHead("SSP ODE solver options");
425: {
426: PetscOptionsFList("-ts_ssp_type","Type of SSP method","TSSSPSetType",TSSSPList,tname,tname,sizeof(tname),&flg);
427: if (flg) {
428: TSSSPSetType(ts,tname);
429: }
430: PetscOptionsInt("-ts_ssp_nstages","Number of stages","TSSSPSetNumStages",ssp->nstages,&ssp->nstages,NULL);
431: }
432: PetscOptionsTail();
433: return(0);
434: }
438: static PetscErrorCode TSView_SSP(TS ts,PetscViewer viewer)
439: {
441: return(0);
442: }
444: /* ------------------------------------------------------------ */
446: /*MC
447: TSSSP - Explicit strong stability preserving ODE solver
449: Most hyperbolic conservation laws have exact solutions that are total variation diminishing (TVD) or total variation
450: bounded (TVB) although these solutions often contain discontinuities. Spatial discretizations such as Godunov's
451: scheme and high-resolution finite volume methods (TVD limiters, ENO/WENO) are designed to preserve these properties,
452: but they are usually formulated using a forward Euler time discretization or by coupling the space and time
453: discretization as in the classical Lax-Wendroff scheme. When the space and time discretization is coupled, it is very
454: difficult to produce schemes with high temporal accuracy while preserving TVD properties. An alternative is the
455: semidiscrete formulation where we choose a spatial discretization that is TVD with forward Euler and then choose a
456: time discretization that preserves the TVD property. Such integrators are called strong stability preserving (SSP).
458: Let c_eff be the minimum number of function evaluations required to step as far as one step of forward Euler while
459: still being SSP. Some theoretical bounds
461: 1. There are no explicit methods with c_eff > 1.
463: 2. There are no explicit methods beyond order 4 (for nonlinear problems) and c_eff > 0.
465: 3. There are no implicit methods with order greater than 1 and c_eff > 2.
467: This integrator provides Runge-Kutta methods of order 2, 3, and 4 with maximal values of c_eff. More stages allows
468: for larger values of c_eff which improves efficiency. These implementations are low-memory and only use 2 or 3 work
469: vectors regardless of the total number of stages, so e.g. 25-stage 3rd order methods may be an excellent choice.
471: Methods can be chosen with -ts_ssp_type {rks2,rks3,rk104}
473: rks2: Second order methods with any number s>1 of stages. c_eff = (s-1)/s
475: rks3: Third order methods with s=n^2 stages, n>1. c_eff = (s-n)/s
477: rk104: A 10-stage fourth order method. c_eff = 0.6
479: Level: beginner
481: References:
482: Ketcheson, Highly efficient strong stability preserving Runge-Kutta methods with low-storage implementations, SISC, 2008.
484: Gottlieb, Ketcheson, and Shu, High order strong stability preserving time discretizations, J Scientific Computing, 2009.
486: .seealso: TSCreate(), TS, TSSetType()
488: M*/
491: PETSC_EXTERN PetscErrorCode TSCreate_SSP(TS ts)
492: {
493: TS_SSP *ssp;
497: TSSSPInitializePackage();
499: ts->ops->setup = TSSetUp_SSP;
500: ts->ops->step = TSStep_SSP;
501: ts->ops->reset = TSReset_SSP;
502: ts->ops->destroy = TSDestroy_SSP;
503: ts->ops->setfromoptions = TSSetFromOptions_SSP;
504: ts->ops->view = TSView_SSP;
506: PetscNewLog(ts,&ssp);
507: ts->data = (void*)ssp;
509: PetscObjectComposeFunction((PetscObject)ts,"TSSSPGetType_C",TSSSPGetType_SSP);
510: PetscObjectComposeFunction((PetscObject)ts,"TSSSPSetType_C",TSSSPSetType_SSP);
511: PetscObjectComposeFunction((PetscObject)ts,"TSSSPGetNumStages_C",TSSSPGetNumStages_SSP);
512: PetscObjectComposeFunction((PetscObject)ts,"TSSSPSetNumStages_C",TSSSPSetNumStages_SSP);
514: TSSSPSetType(ts,TSSSPRKS2);
515: ssp->nstages = 5;
516: return(0);
517: }
521: /*@C
522: TSSSPInitializePackage - This function initializes everything in the TSSSP package. It is called
523: from PetscDLLibraryRegister() when using dynamic libraries, and on the first call to TSCreate_SSP()
524: when using static libraries.
526: Level: developer
528: .keywords: TS, TSSSP, initialize, package
529: .seealso: PetscInitialize()
530: @*/
531: PetscErrorCode TSSSPInitializePackage(void)
532: {
536: if (TSSSPPackageInitialized) return(0);
537: TSSSPPackageInitialized = PETSC_TRUE;
538: PetscFunctionListAdd(&TSSSPList,TSSSPRKS2, TSSSPStep_RK_2);
539: PetscFunctionListAdd(&TSSSPList,TSSSPRKS3, TSSSPStep_RK_3);
540: PetscFunctionListAdd(&TSSSPList,TSSSPRK104,TSSSPStep_RK_10_4);
541: PetscRegisterFinalize(TSSSPFinalizePackage);
542: return(0);
543: }
547: /*@C
548: TSSSPFinalizePackage - This function destroys everything in the TSSSP package. It is
549: called from PetscFinalize().
551: Level: developer
553: .keywords: Petsc, destroy, package
554: .seealso: PetscFinalize()
555: @*/
556: PetscErrorCode TSSSPFinalizePackage(void)
557: {
561: TSSSPPackageInitialized = PETSC_FALSE;
562: PetscFunctionListDestroy(&TSSSPList);
563: return(0);
564: }