Actual source code: ssp.c
1: /*
2: Code for Timestepping with explicit SSP.
3: */
4: #include <petsc/private/tsimpl.h>
6: PetscFunctionList TSSSPList = NULL;
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;
18: static PetscErrorCode TSSSPGetWorkVectors(TS ts,PetscInt n,Vec **work)
19: {
20: TS_SSP *ssp = (TS_SSP*)ts->data;
23: if (ssp->nwork < n) {
24: if (ssp->nwork > 0) {
25: VecDestroyVecs(ssp->nwork,&ssp->work);
26: }
27: VecDuplicateVecs(ts->vec_sol,n,&ssp->work);
28: ssp->nwork = n;
29: }
30: *work = ssp->work;
31: ssp->workout = PETSC_TRUE;
32: return 0;
33: }
35: static PetscErrorCode TSSSPRestoreWorkVectors(TS ts,PetscInt n,Vec **work)
36: {
37: TS_SSP *ssp = (TS_SSP*)ts->data;
41: ssp->workout = PETSC_FALSE;
42: *work = NULL;
43: return 0;
44: }
46: /*MC
47: TSSSPRKS2 - Optimal second order SSP Runge-Kutta method, low-storage, c_eff=(s-1)/s
49: Pseudocode 2 of Ketcheson 2008
51: Level: beginner
53: .seealso: TSSSP, TSSSPSetType(), TSSSPSetNumStages()
54: M*/
55: static PetscErrorCode TSSSPStep_RK_2(TS ts,PetscReal t0,PetscReal dt,Vec sol)
56: {
57: TS_SSP *ssp = (TS_SSP*)ts->data;
58: Vec *work,F;
59: PetscInt i,s;
61: s = ssp->nstages;
62: TSSSPGetWorkVectors(ts,2,&work);
63: F = work[1];
64: VecCopy(sol,work[0]);
65: for (i=0; i<s-1; i++) {
66: PetscReal stage_time = t0+dt*(i/(s-1.));
67: TSPreStage(ts,stage_time);
68: TSComputeRHSFunction(ts,stage_time,work[0],F);
69: VecAXPY(work[0],dt/(s-1.),F);
70: }
71: TSComputeRHSFunction(ts,t0+dt,work[0],F);
72: VecAXPBYPCZ(sol,(s-1.)/s,dt/s,1./s,work[0],F);
73: TSSSPRestoreWorkVectors(ts,2,&work);
74: return 0;
75: }
77: /*MC
78: TSSSPRKS3 - Optimal third order SSP Runge-Kutta, low-storage, c_eff=(PetscSqrtReal(s)-1)/PetscSqrtReal(s), where PetscSqrtReal(s) is an integer
80: Pseudocode 2 of Ketcheson 2008
82: Level: beginner
84: .seealso: TSSSP, TSSSPSetType(), TSSSPSetNumStages()
85: M*/
86: static PetscErrorCode TSSSPStep_RK_3(TS ts,PetscReal t0,PetscReal dt,Vec sol)
87: {
88: TS_SSP *ssp = (TS_SSP*)ts->data;
89: Vec *work,F;
90: PetscInt i,s,n,r;
91: PetscReal c,stage_time;
93: s = ssp->nstages;
94: n = (PetscInt)(PetscSqrtReal((PetscReal)s)+0.001);
95: r = s-n;
97: TSSSPGetWorkVectors(ts,3,&work);
98: F = work[2];
99: VecCopy(sol,work[0]);
100: for (i=0; i<(n-1)*(n-2)/2; i++) {
101: c = (i<n*(n+1)/2) ? 1.*i/(s-n) : (1.*i-n)/(s-n);
102: stage_time = t0+c*dt;
103: TSPreStage(ts,stage_time);
104: TSComputeRHSFunction(ts,stage_time,work[0],F);
105: VecAXPY(work[0],dt/r,F);
106: }
107: VecCopy(work[0],work[1]);
108: for (; i<n*(n+1)/2-1; i++) {
109: c = (i<n*(n+1)/2) ? 1.*i/(s-n) : (1.*i-n)/(s-n);
110: stage_time = t0+c*dt;
111: TSPreStage(ts,stage_time);
112: TSComputeRHSFunction(ts,stage_time,work[0],F);
113: VecAXPY(work[0],dt/r,F);
114: }
115: {
116: c = (i<n*(n+1)/2) ? 1.*i/(s-n) : (1.*i-n)/(s-n);
117: stage_time = t0+c*dt;
118: TSPreStage(ts,stage_time);
119: TSComputeRHSFunction(ts,stage_time,work[0],F);
120: VecAXPBYPCZ(work[0],1.*n/(2*n-1.),(n-1.)*dt/(r*(2*n-1)),(n-1.)/(2*n-1.),work[1],F);
121: i++;
122: }
123: for (; i<s; i++) {
124: c = (i<n*(n+1)/2) ? 1.*i/(s-n) : (1.*i-n)/(s-n);
125: stage_time = t0+c*dt;
126: TSPreStage(ts,stage_time);
127: TSComputeRHSFunction(ts,stage_time,work[0],F);
128: VecAXPY(work[0],dt/r,F);
129: }
130: VecCopy(work[0],sol);
131: TSSSPRestoreWorkVectors(ts,3,&work);
132: return 0;
133: }
135: /*MC
136: TSSSPRKS104 - Optimal fourth order SSP Runge-Kutta, low-storage (2N), c_eff=0.6
138: SSPRK(10,4), Pseudocode 3 of Ketcheson 2008
140: Level: beginner
142: .seealso: TSSSP, TSSSPSetType()
143: M*/
144: static PetscErrorCode TSSSPStep_RK_10_4(TS ts,PetscReal t0,PetscReal dt,Vec sol)
145: {
146: const PetscReal c[10] = {0, 1./6, 2./6, 3./6, 4./6, 2./6, 3./6, 4./6, 5./6, 1};
147: Vec *work,F;
148: PetscInt i;
149: PetscReal stage_time;
151: TSSSPGetWorkVectors(ts,3,&work);
152: F = work[2];
153: VecCopy(sol,work[0]);
154: for (i=0; i<5; i++) {
155: stage_time = t0+c[i]*dt;
156: TSPreStage(ts,stage_time);
157: TSComputeRHSFunction(ts,stage_time,work[0],F);
158: VecAXPY(work[0],dt/6,F);
159: }
160: VecAXPBYPCZ(work[1],1./25,9./25,0,sol,work[0]);
161: VecAXPBY(work[0],15,-5,work[1]);
162: for (; i<9; i++) {
163: stage_time = t0+c[i]*dt;
164: TSPreStage(ts,stage_time);
165: TSComputeRHSFunction(ts,stage_time,work[0],F);
166: VecAXPY(work[0],dt/6,F);
167: }
168: stage_time = t0+dt;
169: TSPreStage(ts,stage_time);
170: TSComputeRHSFunction(ts,stage_time,work[0],F);
171: VecAXPBYPCZ(work[1],3./5,dt/10,1,work[0],F);
172: VecCopy(work[1],sol);
173: TSSSPRestoreWorkVectors(ts,3,&work);
174: return 0;
175: }
177: static PetscErrorCode TSSetUp_SSP(TS ts)
178: {
179: TSCheckImplicitTerm(ts);
180: TSGetAdapt(ts,&ts->adapt);
181: TSAdaptCandidatesClear(ts->adapt);
182: return 0;
183: }
185: static PetscErrorCode TSStep_SSP(TS ts)
186: {
187: TS_SSP *ssp = (TS_SSP*)ts->data;
188: Vec sol = ts->vec_sol;
189: PetscBool stageok,accept = PETSC_TRUE;
190: PetscReal next_time_step = ts->time_step;
192: (*ssp->onestep)(ts,ts->ptime,ts->time_step,sol);
193: TSPostStage(ts,ts->ptime,0,&sol);
194: TSAdaptCheckStage(ts->adapt,ts,ts->ptime+ts->time_step,sol,&stageok);
195: if (!stageok) {ts->reason = TS_DIVERGED_STEP_REJECTED; return 0;}
197: TSAdaptChoose(ts->adapt,ts,ts->time_step,NULL,&next_time_step,&accept);
198: if (!accept) {ts->reason = TS_DIVERGED_STEP_REJECTED; return 0;}
200: ts->ptime += ts->time_step;
201: ts->time_step = next_time_step;
202: return 0;
203: }
204: /*------------------------------------------------------------*/
206: static PetscErrorCode TSReset_SSP(TS ts)
207: {
208: TS_SSP *ssp = (TS_SSP*)ts->data;
210: if (ssp->work) VecDestroyVecs(ssp->nwork,&ssp->work);
211: ssp->nwork = 0;
212: ssp->workout = PETSC_FALSE;
213: return 0;
214: }
216: static PetscErrorCode TSDestroy_SSP(TS ts)
217: {
218: TS_SSP *ssp = (TS_SSP*)ts->data;
220: TSReset_SSP(ts);
221: PetscFree(ssp->type_name);
222: PetscFree(ts->data);
223: PetscObjectComposeFunction((PetscObject)ts,"TSSSPGetType_C",NULL);
224: PetscObjectComposeFunction((PetscObject)ts,"TSSSPSetType_C",NULL);
225: PetscObjectComposeFunction((PetscObject)ts,"TSSSPGetNumStages_C",NULL);
226: PetscObjectComposeFunction((PetscObject)ts,"TSSSPSetNumStages_C",NULL);
227: return 0;
228: }
229: /*------------------------------------------------------------*/
231: /*@C
232: TSSSPSetType - set the SSP time integration scheme to use
234: Logically Collective
236: Input Parameters:
237: + ts - time stepping object
238: - ssptype - type of scheme to use
240: Options Database Keys:
241: -ts_ssp_type <rks2>: Type of SSP method (one of) rks2 rks3 rk104
242: -ts_ssp_nstages <5>: Number of stages
244: Level: beginner
246: .seealso: TSSSP, TSSSPGetType(), TSSSPSetNumStages(), TSSSPRKS2, TSSSPRKS3, TSSSPRK104
247: @*/
248: PetscErrorCode TSSSPSetType(TS ts,TSSSPType ssptype)
249: {
252: PetscTryMethod(ts,"TSSSPSetType_C",(TS,TSSSPType),(ts,ssptype));
253: return 0;
254: }
256: /*@C
257: TSSSPGetType - get the SSP time integration scheme
259: Logically Collective
261: Input Parameter:
262: . ts - time stepping object
264: Output Parameter:
265: . type - type of scheme being used
267: Level: beginner
269: .seealso: TSSSP, TSSSPSettype(), TSSSPSetNumStages(), TSSSPRKS2, TSSSPRKS3, TSSSPRK104
270: @*/
271: PetscErrorCode TSSSPGetType(TS ts,TSSSPType *type)
272: {
274: PetscUseMethod(ts,"TSSSPGetType_C",(TS,TSSSPType*),(ts,type));
275: return 0;
276: }
278: /*@
279: TSSSPSetNumStages - set the number of stages to use with the SSP method
281: Logically Collective
283: Input Parameters:
284: + ts - time stepping object
285: - nstages - number of stages
287: Options Database Keys:
288: -ts_ssp_type <rks2>: NumStages of SSP method (one of) rks2 rks3 rk104
289: -ts_ssp_nstages <5>: Number of stages
291: Level: beginner
293: .seealso: TSSSP, TSSSPGetNumStages(), TSSSPSetNumStages(), TSSSPRKS2, TSSSPRKS3, TSSSPRK104
294: @*/
295: PetscErrorCode TSSSPSetNumStages(TS ts,PetscInt nstages)
296: {
298: PetscTryMethod(ts,"TSSSPSetNumStages_C",(TS,PetscInt),(ts,nstages));
299: return 0;
300: }
302: /*@
303: TSSSPGetNumStages - get the number of stages in the SSP time integration scheme
305: Logically Collective
307: Input Parameter:
308: . ts - time stepping object
310: Output Parameter:
311: . nstages - number of stages
313: Level: beginner
315: .seealso: TSSSP, TSSSPGetType(), TSSSPSetNumStages(), TSSSPRKS2, TSSSPRKS3, TSSSPRK104
316: @*/
317: PetscErrorCode TSSSPGetNumStages(TS ts,PetscInt *nstages)
318: {
320: PetscUseMethod(ts,"TSSSPGetNumStages_C",(TS,PetscInt*),(ts,nstages));
321: return 0;
322: }
324: static PetscErrorCode TSSSPSetType_SSP(TS ts,TSSSPType type)
325: {
326: TS_SSP *ssp = (TS_SSP*)ts->data;
327: PetscErrorCode (*r)(TS,PetscReal,PetscReal,Vec);
329: PetscFunctionListFind(TSSSPList,type,&r);
331: ssp->onestep = r;
332: PetscFree(ssp->type_name);
333: PetscStrallocpy(type,&ssp->type_name);
334: ts->default_adapt_type = TSADAPTNONE;
335: return 0;
336: }
337: static PetscErrorCode TSSSPGetType_SSP(TS ts,TSSSPType *type)
338: {
339: TS_SSP *ssp = (TS_SSP*)ts->data;
341: *type = ssp->type_name;
342: return 0;
343: }
344: static PetscErrorCode TSSSPSetNumStages_SSP(TS ts,PetscInt nstages)
345: {
346: TS_SSP *ssp = (TS_SSP*)ts->data;
348: ssp->nstages = nstages;
349: return 0;
350: }
351: static PetscErrorCode TSSSPGetNumStages_SSP(TS ts,PetscInt *nstages)
352: {
353: TS_SSP *ssp = (TS_SSP*)ts->data;
355: *nstages = ssp->nstages;
356: return 0;
357: }
359: static PetscErrorCode TSSetFromOptions_SSP(PetscOptionItems *PetscOptionsObject,TS ts)
360: {
361: char tname[256] = TSSSPRKS2;
362: TS_SSP *ssp = (TS_SSP*)ts->data;
363: PetscBool flg;
365: PetscOptionsHead(PetscOptionsObject,"SSP ODE solver options");
366: {
367: PetscOptionsFList("-ts_ssp_type","Type of SSP method","TSSSPSetType",TSSSPList,tname,tname,sizeof(tname),&flg);
368: if (flg) {
369: TSSSPSetType(ts,tname);
370: }
371: PetscOptionsInt("-ts_ssp_nstages","Number of stages","TSSSPSetNumStages",ssp->nstages,&ssp->nstages,NULL);
372: }
373: PetscOptionsTail();
374: return 0;
375: }
377: static PetscErrorCode TSView_SSP(TS ts,PetscViewer viewer)
378: {
379: TS_SSP *ssp = (TS_SSP*)ts->data;
380: PetscBool ascii;
382: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&ascii);
383: if (ascii) PetscViewerASCIIPrintf(viewer," Scheme: %s\n",ssp->type_name);
384: return 0;
385: }
387: /* ------------------------------------------------------------ */
389: /*MC
390: TSSSP - Explicit strong stability preserving ODE solver
392: Most hyperbolic conservation laws have exact solutions that are total variation diminishing (TVD) or total variation
393: bounded (TVB) although these solutions often contain discontinuities. Spatial discretizations such as Godunov's
394: scheme and high-resolution finite volume methods (TVD limiters, ENO/WENO) are designed to preserve these properties,
395: but they are usually formulated using a forward Euler time discretization or by coupling the space and time
396: discretization as in the classical Lax-Wendroff scheme. When the space and time discretization is coupled, it is very
397: difficult to produce schemes with high temporal accuracy while preserving TVD properties. An alternative is the
398: semidiscrete formulation where we choose a spatial discretization that is TVD with forward Euler and then choose a
399: time discretization that preserves the TVD property. Such integrators are called strong stability preserving (SSP).
401: Let c_eff be the minimum number of function evaluations required to step as far as one step of forward Euler while
402: still being SSP. Some theoretical bounds
404: 1. There are no explicit methods with c_eff > 1.
406: 2. There are no explicit methods beyond order 4 (for nonlinear problems) and c_eff > 0.
408: 3. There are no implicit methods with order greater than 1 and c_eff > 2.
410: This integrator provides Runge-Kutta methods of order 2, 3, and 4 with maximal values of c_eff. More stages allows
411: for larger values of c_eff which improves efficiency. These implementations are low-memory and only use 2 or 3 work
412: vectors regardless of the total number of stages, so e.g. 25-stage 3rd order methods may be an excellent choice.
414: Methods can be chosen with -ts_ssp_type {rks2,rks3,rk104}
416: rks2: Second order methods with any number s>1 of stages. c_eff = (s-1)/s
418: rks3: Third order methods with s=n^2 stages, n>1. c_eff = (s-n)/s
420: rk104: A 10-stage fourth order method. c_eff = 0.6
422: Level: beginner
424: References:
425: + * - Ketcheson, Highly efficient strong stability preserving Runge Kutta methods with low storage implementations, SISC, 2008.
426: - * - Gottlieb, Ketcheson, and Shu, High order strong stability preserving time discretizations, J Scientific Computing, 2009.
428: .seealso: TSCreate(), TS, TSSetType()
430: M*/
431: PETSC_EXTERN PetscErrorCode TSCreate_SSP(TS ts)
432: {
433: TS_SSP *ssp;
435: TSSSPInitializePackage();
437: ts->ops->setup = TSSetUp_SSP;
438: ts->ops->step = TSStep_SSP;
439: ts->ops->reset = TSReset_SSP;
440: ts->ops->destroy = TSDestroy_SSP;
441: ts->ops->setfromoptions = TSSetFromOptions_SSP;
442: ts->ops->view = TSView_SSP;
444: PetscNewLog(ts,&ssp);
445: ts->data = (void*)ssp;
447: PetscObjectComposeFunction((PetscObject)ts,"TSSSPGetType_C",TSSSPGetType_SSP);
448: PetscObjectComposeFunction((PetscObject)ts,"TSSSPSetType_C",TSSSPSetType_SSP);
449: PetscObjectComposeFunction((PetscObject)ts,"TSSSPGetNumStages_C",TSSSPGetNumStages_SSP);
450: PetscObjectComposeFunction((PetscObject)ts,"TSSSPSetNumStages_C",TSSSPSetNumStages_SSP);
452: TSSSPSetType(ts,TSSSPRKS2);
453: ssp->nstages = 5;
454: return 0;
455: }
457: /*@C
458: TSSSPInitializePackage - This function initializes everything in the TSSSP package. It is called
459: from TSInitializePackage().
461: Level: developer
463: .seealso: PetscInitialize()
464: @*/
465: PetscErrorCode TSSSPInitializePackage(void)
466: {
467: if (TSSSPPackageInitialized) return 0;
468: TSSSPPackageInitialized = PETSC_TRUE;
469: PetscFunctionListAdd(&TSSSPList,TSSSPRKS2, TSSSPStep_RK_2);
470: PetscFunctionListAdd(&TSSSPList,TSSSPRKS3, TSSSPStep_RK_3);
471: PetscFunctionListAdd(&TSSSPList,TSSSPRK104,TSSSPStep_RK_10_4);
472: PetscRegisterFinalize(TSSSPFinalizePackage);
473: return 0;
474: }
476: /*@C
477: TSSSPFinalizePackage - This function destroys everything in the TSSSP package. It is
478: called from PetscFinalize().
480: Level: developer
482: .seealso: PetscFinalize()
483: @*/
484: PetscErrorCode TSSSPFinalizePackage(void)
485: {
486: TSSSPPackageInitialized = PETSC_FALSE;
487: PetscFunctionListDestroy(&TSSSPList);
488: return 0;
489: }