Actual source code: ssp.c
petsc-3.3-p7 2013-05-11
1: /*
2: Code for Timestepping with explicit SSP.
3: */
4: #include <petsc-private/tsimpl.h> /*I "petscts.h" I*/
6: PetscFList TSSSPList = 0;
8: typedef struct {
9: PetscErrorCode (*onestep)(TS,PetscReal,PetscReal,Vec);
10: char *type_name;
11: PetscInt nstages;
12: Vec *work;
13: PetscInt nwork;
14: PetscBool workout;
15: } TS_SSP;
20: static PetscErrorCode TSSSPGetWorkVectors(TS ts,PetscInt n,Vec **work)
21: {
22: TS_SSP *ssp = (TS_SSP*)ts->data;
26: if (ssp->workout) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Work vectors already gotten");
27: if (ssp->nwork < n) {
28: if (ssp->nwork > 0) {
29: VecDestroyVecs(ssp->nwork,&ssp->work);
30: }
31: VecDuplicateVecs(ts->vec_sol,n,&ssp->work);
32: ssp->nwork = n;
33: }
34: *work = ssp->work;
35: ssp->workout = PETSC_TRUE;
36: return(0);
37: }
41: static PetscErrorCode TSSSPRestoreWorkVectors(TS ts,PetscInt n,Vec **work)
42: {
43: TS_SSP *ssp = (TS_SSP*)ts->data;
46: if (!ssp->workout) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Work vectors have not been gotten");
47: if (*work != ssp->work) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Wrong work vectors checked out");
48: ssp->workout = PETSC_FALSE;
49: *work = PETSC_NULL;
50: return(0);
51: }
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;
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: PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSSSPGetType_C","",PETSC_NULL);
249: PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSSSPSetType_C","",PETSC_NULL);
250: PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSSSPGetNumStages_C","",PETSC_NULL);
251: PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSSSPSetNumStages_C","",PETSC_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,const TSSSPType type)
276: {
281: PetscTryMethod(ts,"TSSSPSetType_C",(TS,const 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,const TSSSPType *type)
303: {
308: PetscTryMethod(ts,"TSSSPGetType_C",(TS,const 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: }
368: EXTERN_C_BEGIN
371: PetscErrorCode TSSSPSetType_SSP(TS ts,const TSSSPType type)
372: {
373: PetscErrorCode ierr,(*r)(TS,PetscReal,PetscReal,Vec);
374: TS_SSP *ssp = (TS_SSP*)ts->data;
377: PetscFListFind(TSSSPList,((PetscObject)ts)->comm,type,PETSC_TRUE,(PetscVoidStarFunction)&r);
378: if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown TS_SSP type %s given",type);
379: ssp->onestep = r;
380: PetscFree(ssp->type_name);
381: PetscStrallocpy(type,&ssp->type_name);
382: return(0);
383: }
386: PetscErrorCode TSSSPGetType_SSP(TS ts,const TSSSPType *type)
387: {
388: TS_SSP *ssp = (TS_SSP*)ts->data;
391: *type = ssp->type_name;
392: return(0);
393: }
396: PetscErrorCode TSSSPSetNumStages_SSP(TS ts,PetscInt nstages)
397: {
398: TS_SSP *ssp = (TS_SSP*)ts->data;
401: ssp->nstages = nstages;
402: return(0);
403: }
406: PetscErrorCode TSSSPGetNumStages_SSP(TS ts,PetscInt *nstages)
407: {
408: TS_SSP *ssp = (TS_SSP*)ts->data;
411: *nstages = ssp->nstages;
412: return(0);
413: }
414: EXTERN_C_END
418: static PetscErrorCode TSSetFromOptions_SSP(TS ts)
419: {
420: char tname[256] = TSSSPRKS2;
421: TS_SSP *ssp = (TS_SSP*)ts->data;
423: PetscBool flg;
426: PetscOptionsHead("SSP ODE solver options");
427: {
428: PetscOptionsList("-ts_ssp_type","Type of SSP method","TSSSPSetType",TSSSPList,tname,tname,sizeof(tname),&flg);
429: if (flg) {
430: TSSSPSetType(ts,tname);
431: }
432: PetscOptionsInt("-ts_ssp_nstages","Number of stages","TSSSPSetNumStages",ssp->nstages,&ssp->nstages,PETSC_NULL);
433: }
434: PetscOptionsTail();
435: return(0);
436: }
440: static PetscErrorCode TSView_SSP(TS ts,PetscViewer viewer)
441: {
443: return(0);
444: }
446: /* ------------------------------------------------------------ */
448: /*MC
449: TSSSP - Explicit strong stability preserving ODE solver
451: Most hyperbolic conservation laws have exact solutions that are total variation diminishing (TVD) or total variation
452: bounded (TVB) although these solutions often contain discontinuities. Spatial discretizations such as Godunov's
453: scheme and high-resolution finite volume methods (TVD limiters, ENO/WENO) are designed to preserve these properties,
454: but they are usually formulated using a forward Euler time discretization or by coupling the space and time
455: discretization as in the classical Lax-Wendroff scheme. When the space and time discretization is coupled, it is very
456: difficult to produce schemes with high temporal accuracy while preserving TVD properties. An alternative is the
457: semidiscrete formulation where we choose a spatial discretization that is TVD with forward Euler and then choose a
458: time discretization that preserves the TVD property. Such integrators are called strong stability preserving (SSP).
460: Let c_eff be the minimum number of function evaluations required to step as far as one step of forward Euler while
461: still being SSP. Some theoretical bounds
463: 1. There are no explicit methods with c_eff > 1.
465: 2. There are no explicit methods beyond order 4 (for nonlinear problems) and c_eff > 0.
467: 3. There are no implicit methods with order greater than 1 and c_eff > 2.
469: This integrator provides Runge-Kutta methods of order 2, 3, and 4 with maximal values of c_eff. More stages allows
470: for larger values of c_eff which improves efficiency. These implementations are low-memory and only use 2 or 3 work
471: vectors regardless of the total number of stages, so e.g. 25-stage 3rd order methods may be an excellent choice.
473: Methods can be chosen with -ts_ssp_type {rks2,rks3,rk104}
475: rks2: Second order methods with any number s>1 of stages. c_eff = (s-1)/s
477: rks3: Third order methods with s=n^2 stages, n>1. c_eff = (s-n)/s
479: rk104: A 10-stage fourth order method. c_eff = 0.6
481: Level: beginner
483: References:
484: Ketcheson, Highly efficient strong stability preserving Runge-Kutta methods with low-storage implementations, SISC, 2008.
486: Gottlieb, Ketcheson, and Shu, High order strong stability preserving time discretizations, J Scientific Computing, 2009.
488: .seealso: TSCreate(), TS, TSSetType()
490: M*/
491: EXTERN_C_BEGIN
494: PetscErrorCode TSCreate_SSP(TS ts)
495: {
496: TS_SSP *ssp;
500: if (!TSSSPList) {
501: PetscFListAdd(&TSSSPList,TSSSPRKS2, "TSSSPStep_RK_2", (void(*)(void))TSSSPStep_RK_2);
502: PetscFListAdd(&TSSSPList,TSSSPRKS3, "TSSSPStep_RK_3", (void(*)(void))TSSSPStep_RK_3);
503: PetscFListAdd(&TSSSPList,TSSSPRK104, "TSSSPStep_RK_10_4",(void(*)(void))TSSSPStep_RK_10_4);
504: }
506: ts->ops->setup = TSSetUp_SSP;
507: ts->ops->step = TSStep_SSP;
508: ts->ops->reset = TSReset_SSP;
509: ts->ops->destroy = TSDestroy_SSP;
510: ts->ops->setfromoptions = TSSetFromOptions_SSP;
511: ts->ops->view = TSView_SSP;
513: PetscNewLog(ts,TS_SSP,&ssp);
514: ts->data = (void*)ssp;
516: PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSSSPGetType_C","TSSSPGetType_SSP",TSSSPGetType_SSP);
517: PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSSSPSetType_C","TSSSPSetType_SSP",TSSSPSetType_SSP);
518: PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSSSPGetNumStages_C","TSSSPGetNumStages_SSP",TSSSPGetNumStages_SSP);
519: PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSSSPSetNumStages_C","TSSSPSetNumStages_SSP",TSSSPSetNumStages_SSP);
521: TSSSPSetType(ts,TSSSPRKS2);
522: ssp->nstages = 5;
523: return(0);
524: }
525: EXTERN_C_END