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;
19: static PetscErrorCode TSSSPGetWorkVectors(TS ts,PetscInt n,Vec **work)
20: {
21: TS_SSP *ssp = (TS_SSP*)ts->data;
25: if (ssp->workout) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Work vectors already gotten");
26: if (ssp->nwork < n) {
27: if (ssp->nwork > 0) {
28: VecDestroyVecs(ssp->nwork,&ssp->work);
29: }
30: VecDuplicateVecs(ts->vec_sol,n,&ssp->work);
31: ssp->nwork = n;
32: }
33: *work = ssp->work;
34: ssp->workout = PETSC_TRUE;
35: return(0);
36: }
38: static PetscErrorCode TSSSPRestoreWorkVectors(TS ts,PetscInt n,Vec **work)
39: {
40: TS_SSP *ssp = (TS_SSP*)ts->data;
43: if (!ssp->workout) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Work vectors have not been gotten");
44: if (*work != ssp->work) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Wrong work vectors checked out");
45: ssp->workout = PETSC_FALSE;
46: *work = NULL;
47: return(0);
48: }
50: /*MC
51: TSSSPRKS2 - Optimal second order SSP Runge-Kutta method, low-storage, c_eff=(s-1)/s
53: Pseudocode 2 of Ketcheson 2008
55: Level: beginner
57: .seealso: TSSSP, TSSSPSetType(), TSSSPSetNumStages()
58: M*/
59: static PetscErrorCode TSSSPStep_RK_2(TS ts,PetscReal t0,PetscReal dt,Vec sol)
60: {
61: TS_SSP *ssp = (TS_SSP*)ts->data;
62: Vec *work,F;
63: PetscInt i,s;
67: s = ssp->nstages;
68: TSSSPGetWorkVectors(ts,2,&work);
69: F = work[1];
70: VecCopy(sol,work[0]);
71: for (i=0; i<s-1; i++) {
72: PetscReal stage_time = t0+dt*(i/(s-1.));
73: TSPreStage(ts,stage_time);
74: TSComputeRHSFunction(ts,stage_time,work[0],F);
75: VecAXPY(work[0],dt/(s-1.),F);
76: }
77: TSComputeRHSFunction(ts,t0+dt,work[0],F);
78: VecAXPBYPCZ(sol,(s-1.)/s,dt/s,1./s,work[0],F);
79: TSSSPRestoreWorkVectors(ts,2,&work);
80: return(0);
81: }
83: /*MC
84: TSSSPRKS3 - Optimal third order SSP Runge-Kutta, low-storage, c_eff=(PetscSqrtReal(s)-1)/PetscSqrtReal(s), where PetscSqrtReal(s) is an integer
86: Pseudocode 2 of Ketcheson 2008
88: Level: beginner
90: .seealso: TSSSP, TSSSPSetType(), TSSSPSetNumStages()
91: M*/
92: static PetscErrorCode TSSSPStep_RK_3(TS ts,PetscReal t0,PetscReal dt,Vec sol)
93: {
94: TS_SSP *ssp = (TS_SSP*)ts->data;
95: Vec *work,F;
96: PetscInt i,s,n,r;
97: PetscReal c,stage_time;
101: s = ssp->nstages;
102: n = (PetscInt)(PetscSqrtReal((PetscReal)s)+0.001);
103: r = s-n;
104: 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);
105: TSSSPGetWorkVectors(ts,3,&work);
106: F = work[2];
107: VecCopy(sol,work[0]);
108: for (i=0; i<(n-1)*(n-2)/2; 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: VecCopy(work[0],work[1]);
116: for (; i<n*(n+1)/2-1; 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: {
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: VecAXPBYPCZ(work[0],1.*n/(2*n-1.),(n-1.)*dt/(r*(2*n-1)),(n-1.)/(2*n-1.),work[1],F);
129: i++;
130: }
131: for (; i<s; i++) {
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: VecAXPY(work[0],dt/r,F);
137: }
138: VecCopy(work[0],sol);
139: TSSSPRestoreWorkVectors(ts,3,&work);
140: return(0);
141: }
143: /*MC
144: TSSSPRKS104 - Optimal fourth order SSP Runge-Kutta, low-storage (2N), c_eff=0.6
146: SSPRK(10,4), Pseudocode 3 of Ketcheson 2008
148: Level: beginner
150: .seealso: TSSSP, TSSSPSetType()
151: M*/
152: static PetscErrorCode TSSSPStep_RK_10_4(TS ts,PetscReal t0,PetscReal dt,Vec sol)
153: {
154: const PetscReal c[10] = {0, 1./6, 2./6, 3./6, 4./6, 2./6, 3./6, 4./6, 5./6, 1};
155: Vec *work,F;
156: PetscInt i;
157: PetscReal stage_time;
158: PetscErrorCode ierr;
161: TSSSPGetWorkVectors(ts,3,&work);
162: F = work[2];
163: VecCopy(sol,work[0]);
164: for (i=0; i<5; i++) {
165: stage_time = t0+c[i]*dt;
166: TSPreStage(ts,stage_time);
167: TSComputeRHSFunction(ts,stage_time,work[0],F);
168: VecAXPY(work[0],dt/6,F);
169: }
170: VecAXPBYPCZ(work[1],1./25,9./25,0,sol,work[0]);
171: VecAXPBY(work[0],15,-5,work[1]);
172: for (; i<9; i++) {
173: stage_time = t0+c[i]*dt;
174: TSPreStage(ts,stage_time);
175: TSComputeRHSFunction(ts,stage_time,work[0],F);
176: VecAXPY(work[0],dt/6,F);
177: }
178: stage_time = t0+dt;
179: TSPreStage(ts,stage_time);
180: TSComputeRHSFunction(ts,stage_time,work[0],F);
181: VecAXPBYPCZ(work[1],3./5,dt/10,1,work[0],F);
182: VecCopy(work[1],sol);
183: TSSSPRestoreWorkVectors(ts,3,&work);
184: return(0);
185: }
188: static PetscErrorCode TSSetUp_SSP(TS ts)
189: {
193: TSCheckImplicitTerm(ts);
194: TSGetAdapt(ts,&ts->adapt);
195: TSAdaptCandidatesClear(ts->adapt);
196: return(0);
197: }
199: static PetscErrorCode TSStep_SSP(TS ts)
200: {
201: TS_SSP *ssp = (TS_SSP*)ts->data;
202: Vec sol = ts->vec_sol;
203: PetscBool stageok,accept = PETSC_TRUE;
204: PetscReal next_time_step = ts->time_step;
208: (*ssp->onestep)(ts,ts->ptime,ts->time_step,sol);
209: TSPostStage(ts,ts->ptime,0,&sol);
210: TSAdaptCheckStage(ts->adapt,ts,ts->ptime+ts->time_step,sol,&stageok);
211: if (!stageok) {ts->reason = TS_DIVERGED_STEP_REJECTED; return(0);}
213: TSAdaptChoose(ts->adapt,ts,ts->time_step,NULL,&next_time_step,&accept);
214: if (!accept) {ts->reason = TS_DIVERGED_STEP_REJECTED; return(0);}
216: ts->ptime += ts->time_step;
217: ts->time_step = next_time_step;
218: return(0);
219: }
220: /*------------------------------------------------------------*/
222: static PetscErrorCode TSReset_SSP(TS ts)
223: {
224: TS_SSP *ssp = (TS_SSP*)ts->data;
228: if (ssp->work) {VecDestroyVecs(ssp->nwork,&ssp->work);}
229: ssp->nwork = 0;
230: ssp->workout = PETSC_FALSE;
231: return(0);
232: }
234: static PetscErrorCode TSDestroy_SSP(TS ts)
235: {
236: TS_SSP *ssp = (TS_SSP*)ts->data;
240: TSReset_SSP(ts);
241: PetscFree(ssp->type_name);
242: PetscFree(ts->data);
243: PetscObjectComposeFunction((PetscObject)ts,"TSSSPGetType_C",NULL);
244: PetscObjectComposeFunction((PetscObject)ts,"TSSSPSetType_C",NULL);
245: PetscObjectComposeFunction((PetscObject)ts,"TSSSPGetNumStages_C",NULL);
246: PetscObjectComposeFunction((PetscObject)ts,"TSSSPSetNumStages_C",NULL);
247: return(0);
248: }
249: /*------------------------------------------------------------*/
251: /*@C
252: TSSSPSetType - set the SSP time integration scheme to use
254: Logically Collective
256: Input Arguments:
257: ts - time stepping object
258: ssptype - type of scheme to use
260: Options Database Keys:
261: -ts_ssp_type <rks2>: Type of SSP method (one of) rks2 rks3 rk104
262: -ts_ssp_nstages <5>: Number of stages
264: Level: beginner
266: .seealso: TSSSP, TSSSPGetType(), TSSSPSetNumStages(), TSSSPRKS2, TSSSPRKS3, TSSSPRK104
267: @*/
268: PetscErrorCode TSSSPSetType(TS ts,TSSSPType ssptype)
269: {
275: PetscTryMethod(ts,"TSSSPSetType_C",(TS,TSSSPType),(ts,ssptype));
276: return(0);
277: }
279: /*@C
280: TSSSPGetType - get the SSP time integration scheme
282: Logically Collective
284: Input Argument:
285: ts - time stepping object
287: Output Argument:
288: type - type of scheme being used
290: Level: beginner
292: .seealso: TSSSP, TSSSPSettype(), TSSSPSetNumStages(), TSSSPRKS2, TSSSPRKS3, TSSSPRK104
293: @*/
294: PetscErrorCode TSSSPGetType(TS ts,TSSSPType *type)
295: {
300: PetscUseMethod(ts,"TSSSPGetType_C",(TS,TSSSPType*),(ts,type));
301: return(0);
302: }
304: /*@
305: TSSSPSetNumStages - set the number of stages to use with the SSP method
307: Logically Collective
309: Input Arguments:
310: ts - time stepping object
311: nstages - number of stages
313: Options Database Keys:
314: -ts_ssp_type <rks2>: NumStages of SSP method (one of) rks2 rks3 rk104
315: -ts_ssp_nstages <5>: Number of stages
317: Level: beginner
319: .seealso: TSSSP, TSSSPGetNumStages(), TSSSPSetNumStages(), TSSSPRKS2, TSSSPRKS3, TSSSPRK104
320: @*/
321: PetscErrorCode TSSSPSetNumStages(TS ts,PetscInt nstages)
322: {
327: PetscTryMethod(ts,"TSSSPSetNumStages_C",(TS,PetscInt),(ts,nstages));
328: return(0);
329: }
331: /*@
332: TSSSPGetNumStages - get the number of stages in the SSP time integration scheme
334: Logically Collective
336: Input Argument:
337: ts - time stepping object
339: Output Argument:
340: nstages - number of stages
342: Level: beginner
344: .seealso: TSSSP, TSSSPGetType(), TSSSPSetNumStages(), TSSSPRKS2, TSSSPRKS3, TSSSPRK104
345: @*/
346: PetscErrorCode TSSSPGetNumStages(TS ts,PetscInt *nstages)
347: {
352: PetscUseMethod(ts,"TSSSPGetNumStages_C",(TS,PetscInt*),(ts,nstages));
353: return(0);
354: }
356: static PetscErrorCode TSSSPSetType_SSP(TS ts,TSSSPType type)
357: {
358: PetscErrorCode ierr,(*r)(TS,PetscReal,PetscReal,Vec);
359: TS_SSP *ssp = (TS_SSP*)ts->data;
362: PetscFunctionListFind(TSSSPList,type,&r);
363: if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown TS_SSP type %s given",type);
364: ssp->onestep = r;
365: PetscFree(ssp->type_name);
366: PetscStrallocpy(type,&ssp->type_name);
367: ts->default_adapt_type = TSADAPTNONE;
368: return(0);
369: }
370: static PetscErrorCode TSSSPGetType_SSP(TS ts,TSSSPType *type)
371: {
372: TS_SSP *ssp = (TS_SSP*)ts->data;
375: *type = ssp->type_name;
376: return(0);
377: }
378: static PetscErrorCode TSSSPSetNumStages_SSP(TS ts,PetscInt nstages)
379: {
380: TS_SSP *ssp = (TS_SSP*)ts->data;
383: ssp->nstages = nstages;
384: return(0);
385: }
386: static PetscErrorCode TSSSPGetNumStages_SSP(TS ts,PetscInt *nstages)
387: {
388: TS_SSP *ssp = (TS_SSP*)ts->data;
391: *nstages = ssp->nstages;
392: return(0);
393: }
395: static PetscErrorCode TSSetFromOptions_SSP(PetscOptionItems *PetscOptionsObject,TS ts)
396: {
397: char tname[256] = TSSSPRKS2;
398: TS_SSP *ssp = (TS_SSP*)ts->data;
400: PetscBool flg;
403: PetscOptionsHead(PetscOptionsObject,"SSP ODE solver options");
404: {
405: PetscOptionsFList("-ts_ssp_type","Type of SSP method","TSSSPSetType",TSSSPList,tname,tname,sizeof(tname),&flg);
406: if (flg) {
407: TSSSPSetType(ts,tname);
408: }
409: PetscOptionsInt("-ts_ssp_nstages","Number of stages","TSSSPSetNumStages",ssp->nstages,&ssp->nstages,NULL);
410: }
411: PetscOptionsTail();
412: return(0);
413: }
415: static PetscErrorCode TSView_SSP(TS ts,PetscViewer viewer)
416: {
417: TS_SSP *ssp = (TS_SSP*)ts->data;
418: PetscBool ascii;
422: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&ascii);
423: if (ascii) {PetscViewerASCIIPrintf(viewer," Scheme: %s\n",ssp->type_name);}
424: return(0);
425: }
427: /* ------------------------------------------------------------ */
429: /*MC
430: TSSSP - Explicit strong stability preserving ODE solver
432: Most hyperbolic conservation laws have exact solutions that are total variation diminishing (TVD) or total variation
433: bounded (TVB) although these solutions often contain discontinuities. Spatial discretizations such as Godunov's
434: scheme and high-resolution finite volume methods (TVD limiters, ENO/WENO) are designed to preserve these properties,
435: but they are usually formulated using a forward Euler time discretization or by coupling the space and time
436: discretization as in the classical Lax-Wendroff scheme. When the space and time discretization is coupled, it is very
437: difficult to produce schemes with high temporal accuracy while preserving TVD properties. An alternative is the
438: semidiscrete formulation where we choose a spatial discretization that is TVD with forward Euler and then choose a
439: time discretization that preserves the TVD property. Such integrators are called strong stability preserving (SSP).
441: Let c_eff be the minimum number of function evaluations required to step as far as one step of forward Euler while
442: still being SSP. Some theoretical bounds
444: 1. There are no explicit methods with c_eff > 1.
446: 2. There are no explicit methods beyond order 4 (for nonlinear problems) and c_eff > 0.
448: 3. There are no implicit methods with order greater than 1 and c_eff > 2.
450: This integrator provides Runge-Kutta methods of order 2, 3, and 4 with maximal values of c_eff. More stages allows
451: for larger values of c_eff which improves efficiency. These implementations are low-memory and only use 2 or 3 work
452: vectors regardless of the total number of stages, so e.g. 25-stage 3rd order methods may be an excellent choice.
454: Methods can be chosen with -ts_ssp_type {rks2,rks3,rk104}
456: rks2: Second order methods with any number s>1 of stages. c_eff = (s-1)/s
458: rks3: Third order methods with s=n^2 stages, n>1. c_eff = (s-n)/s
460: rk104: A 10-stage fourth order method. c_eff = 0.6
462: Level: beginner
464: References:
465: + 1. - Ketcheson, Highly efficient strong stability preserving Runge Kutta methods with low storage implementations, SISC, 2008.
466: - 2. - Gottlieb, Ketcheson, and Shu, High order strong stability preserving time discretizations, J Scientific Computing, 2009.
468: .seealso: TSCreate(), TS, TSSetType()
470: M*/
471: PETSC_EXTERN PetscErrorCode TSCreate_SSP(TS ts)
472: {
473: TS_SSP *ssp;
477: TSSSPInitializePackage();
479: ts->ops->setup = TSSetUp_SSP;
480: ts->ops->step = TSStep_SSP;
481: ts->ops->reset = TSReset_SSP;
482: ts->ops->destroy = TSDestroy_SSP;
483: ts->ops->setfromoptions = TSSetFromOptions_SSP;
484: ts->ops->view = TSView_SSP;
486: PetscNewLog(ts,&ssp);
487: ts->data = (void*)ssp;
489: PetscObjectComposeFunction((PetscObject)ts,"TSSSPGetType_C",TSSSPGetType_SSP);
490: PetscObjectComposeFunction((PetscObject)ts,"TSSSPSetType_C",TSSSPSetType_SSP);
491: PetscObjectComposeFunction((PetscObject)ts,"TSSSPGetNumStages_C",TSSSPGetNumStages_SSP);
492: PetscObjectComposeFunction((PetscObject)ts,"TSSSPSetNumStages_C",TSSSPSetNumStages_SSP);
494: TSSSPSetType(ts,TSSSPRKS2);
495: ssp->nstages = 5;
496: return(0);
497: }
499: /*@C
500: TSSSPInitializePackage - This function initializes everything in the TSSSP package. It is called
501: from TSInitializePackage().
503: Level: developer
505: .seealso: PetscInitialize()
506: @*/
507: PetscErrorCode TSSSPInitializePackage(void)
508: {
512: if (TSSSPPackageInitialized) return(0);
513: TSSSPPackageInitialized = PETSC_TRUE;
514: PetscFunctionListAdd(&TSSSPList,TSSSPRKS2, TSSSPStep_RK_2);
515: PetscFunctionListAdd(&TSSSPList,TSSSPRKS3, TSSSPStep_RK_3);
516: PetscFunctionListAdd(&TSSSPList,TSSSPRK104,TSSSPStep_RK_10_4);
517: PetscRegisterFinalize(TSSSPFinalizePackage);
518: return(0);
519: }
521: /*@C
522: TSSSPFinalizePackage - This function destroys everything in the TSSSP package. It is
523: called from PetscFinalize().
525: Level: developer
527: .seealso: PetscFinalize()
528: @*/
529: PetscErrorCode TSSSPFinalizePackage(void)
530: {
534: TSSSPPackageInitialized = PETSC_FALSE;
535: PetscFunctionListDestroy(&TSSSPList);
536: return(0);
537: }