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;
24: if (ssp->workout) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Work vectors already gotten");
25: if (ssp->nwork < n) {
26: if (ssp->nwork > 0) {
27: VecDestroyVecs(ssp->nwork,&ssp->work);
28: }
29: VecDuplicateVecs(ts->vec_sol,n,&ssp->work);
30: ssp->nwork = n;
31: }
32: *work = ssp->work;
33: ssp->workout = PETSC_TRUE;
34: return(0);
35: }
37: static PetscErrorCode TSSSPRestoreWorkVectors(TS ts,PetscInt n,Vec **work)
38: {
39: TS_SSP *ssp = (TS_SSP*)ts->data;
42: if (!ssp->workout) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Work vectors have not been gotten");
43: if (*work != ssp->work) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Wrong work vectors checked out");
44: ssp->workout = PETSC_FALSE;
45: *work = NULL;
46: return(0);
47: }
49: /*MC
50: TSSSPRKS2 - Optimal second order SSP Runge-Kutta method, low-storage, c_eff=(s-1)/s
52: Pseudocode 2 of Ketcheson 2008
54: Level: beginner
56: .seealso: TSSSP, TSSSPSetType(), TSSSPSetNumStages()
57: M*/
58: static PetscErrorCode TSSSPStep_RK_2(TS ts,PetscReal t0,PetscReal dt,Vec sol)
59: {
60: TS_SSP *ssp = (TS_SSP*)ts->data;
61: Vec *work,F;
62: PetscInt i,s;
66: s = ssp->nstages;
67: TSSSPGetWorkVectors(ts,2,&work);
68: F = work[1];
69: VecCopy(sol,work[0]);
70: for (i=0; i<s-1; i++) {
71: PetscReal stage_time = t0+dt*(i/(s-1.));
72: TSPreStage(ts,stage_time);
73: TSComputeRHSFunction(ts,stage_time,work[0],F);
74: VecAXPY(work[0],dt/(s-1.),F);
75: }
76: TSComputeRHSFunction(ts,t0+dt,work[0],F);
77: VecAXPBYPCZ(sol,(s-1.)/s,dt/s,1./s,work[0],F);
78: TSSSPRestoreWorkVectors(ts,2,&work);
79: return(0);
80: }
82: /*MC
83: TSSSPRKS3 - Optimal third order SSP Runge-Kutta, low-storage, c_eff=(PetscSqrtReal(s)-1)/PetscSqrtReal(s), where PetscSqrtReal(s) is an integer
85: Pseudocode 2 of Ketcheson 2008
87: Level: beginner
89: .seealso: TSSSP, TSSSPSetType(), TSSSPSetNumStages()
90: M*/
91: static PetscErrorCode TSSSPStep_RK_3(TS ts,PetscReal t0,PetscReal dt,Vec sol)
92: {
93: TS_SSP *ssp = (TS_SSP*)ts->data;
94: Vec *work,F;
95: PetscInt i,s,n,r;
96: PetscReal c,stage_time;
100: s = ssp->nstages;
101: n = (PetscInt)(PetscSqrtReal((PetscReal)s)+0.001);
102: r = s-n;
103: 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);
104: TSSSPGetWorkVectors(ts,3,&work);
105: F = work[2];
106: VecCopy(sol,work[0]);
107: for (i=0; i<(n-1)*(n-2)/2; i++) {
108: c = (i<n*(n+1)/2) ? 1.*i/(s-n) : (1.*i-n)/(s-n);
109: stage_time = t0+c*dt;
110: TSPreStage(ts,stage_time);
111: TSComputeRHSFunction(ts,stage_time,work[0],F);
112: VecAXPY(work[0],dt/r,F);
113: }
114: VecCopy(work[0],work[1]);
115: for (; i<n*(n+1)/2-1; i++) {
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: VecAXPY(work[0],dt/r,F);
121: }
122: {
123: c = (i<n*(n+1)/2) ? 1.*i/(s-n) : (1.*i-n)/(s-n);
124: stage_time = t0+c*dt;
125: TSPreStage(ts,stage_time);
126: TSComputeRHSFunction(ts,stage_time,work[0],F);
127: VecAXPBYPCZ(work[0],1.*n/(2*n-1.),(n-1.)*dt/(r*(2*n-1)),(n-1.)/(2*n-1.),work[1],F);
128: i++;
129: }
130: for (; i<s; i++) {
131: c = (i<n*(n+1)/2) ? 1.*i/(s-n) : (1.*i-n)/(s-n);
132: stage_time = t0+c*dt;
133: TSPreStage(ts,stage_time);
134: TSComputeRHSFunction(ts,stage_time,work[0],F);
135: VecAXPY(work[0],dt/r,F);
136: }
137: VecCopy(work[0],sol);
138: TSSSPRestoreWorkVectors(ts,3,&work);
139: return(0);
140: }
142: /*MC
143: TSSSPRKS104 - Optimal fourth order SSP Runge-Kutta, low-storage (2N), c_eff=0.6
145: SSPRK(10,4), Pseudocode 3 of Ketcheson 2008
147: Level: beginner
149: .seealso: TSSSP, TSSSPSetType()
150: M*/
151: static PetscErrorCode TSSSPStep_RK_10_4(TS ts,PetscReal t0,PetscReal dt,Vec sol)
152: {
153: const PetscReal c[10] = {0, 1./6, 2./6, 3./6, 4./6, 2./6, 3./6, 4./6, 5./6, 1};
154: Vec *work,F;
155: PetscInt i;
156: PetscReal stage_time;
157: PetscErrorCode ierr;
160: TSSSPGetWorkVectors(ts,3,&work);
161: F = work[2];
162: VecCopy(sol,work[0]);
163: for (i=0; i<5; i++) {
164: stage_time = t0+c[i]*dt;
165: TSPreStage(ts,stage_time);
166: TSComputeRHSFunction(ts,stage_time,work[0],F);
167: VecAXPY(work[0],dt/6,F);
168: }
169: VecAXPBYPCZ(work[1],1./25,9./25,0,sol,work[0]);
170: VecAXPBY(work[0],15,-5,work[1]);
171: for (; i<9; i++) {
172: stage_time = t0+c[i]*dt;
173: TSPreStage(ts,stage_time);
174: TSComputeRHSFunction(ts,stage_time,work[0],F);
175: VecAXPY(work[0],dt/6,F);
176: }
177: stage_time = t0+dt;
178: TSPreStage(ts,stage_time);
179: TSComputeRHSFunction(ts,stage_time,work[0],F);
180: VecAXPBYPCZ(work[1],3./5,dt/10,1,work[0],F);
181: VecCopy(work[1],sol);
182: TSSSPRestoreWorkVectors(ts,3,&work);
183: return(0);
184: }
186: static PetscErrorCode TSSetUp_SSP(TS ts)
187: {
191: TSCheckImplicitTerm(ts);
192: TSGetAdapt(ts,&ts->adapt);
193: TSAdaptCandidatesClear(ts->adapt);
194: return(0);
195: }
197: static PetscErrorCode TSStep_SSP(TS ts)
198: {
199: TS_SSP *ssp = (TS_SSP*)ts->data;
200: Vec sol = ts->vec_sol;
201: PetscBool stageok,accept = PETSC_TRUE;
202: PetscReal next_time_step = ts->time_step;
206: (*ssp->onestep)(ts,ts->ptime,ts->time_step,sol);
207: TSPostStage(ts,ts->ptime,0,&sol);
208: TSAdaptCheckStage(ts->adapt,ts,ts->ptime+ts->time_step,sol,&stageok);
209: if (!stageok) {ts->reason = TS_DIVERGED_STEP_REJECTED; return(0);}
211: TSAdaptChoose(ts->adapt,ts,ts->time_step,NULL,&next_time_step,&accept);
212: if (!accept) {ts->reason = TS_DIVERGED_STEP_REJECTED; return(0);}
214: ts->ptime += ts->time_step;
215: ts->time_step = next_time_step;
216: return(0);
217: }
218: /*------------------------------------------------------------*/
220: static PetscErrorCode TSReset_SSP(TS ts)
221: {
222: TS_SSP *ssp = (TS_SSP*)ts->data;
226: if (ssp->work) {VecDestroyVecs(ssp->nwork,&ssp->work);}
227: ssp->nwork = 0;
228: ssp->workout = PETSC_FALSE;
229: return(0);
230: }
232: static PetscErrorCode TSDestroy_SSP(TS ts)
233: {
234: TS_SSP *ssp = (TS_SSP*)ts->data;
238: TSReset_SSP(ts);
239: PetscFree(ssp->type_name);
240: PetscFree(ts->data);
241: PetscObjectComposeFunction((PetscObject)ts,"TSSSPGetType_C",NULL);
242: PetscObjectComposeFunction((PetscObject)ts,"TSSSPSetType_C",NULL);
243: PetscObjectComposeFunction((PetscObject)ts,"TSSSPGetNumStages_C",NULL);
244: PetscObjectComposeFunction((PetscObject)ts,"TSSSPSetNumStages_C",NULL);
245: return(0);
246: }
247: /*------------------------------------------------------------*/
249: /*@C
250: TSSSPSetType - set the SSP time integration scheme to use
252: Logically Collective
254: Input Parameters:
255: + ts - time stepping object
256: - ssptype - type of scheme to use
258: Options Database Keys:
259: -ts_ssp_type <rks2>: Type of SSP method (one of) rks2 rks3 rk104
260: -ts_ssp_nstages <5>: Number of stages
262: Level: beginner
264: .seealso: TSSSP, TSSSPGetType(), TSSSPSetNumStages(), TSSSPRKS2, TSSSPRKS3, TSSSPRK104
265: @*/
266: PetscErrorCode TSSSPSetType(TS ts,TSSSPType ssptype)
267: {
273: PetscTryMethod(ts,"TSSSPSetType_C",(TS,TSSSPType),(ts,ssptype));
274: return(0);
275: }
277: /*@C
278: TSSSPGetType - get the SSP time integration scheme
280: Logically Collective
282: Input Parameter:
283: . ts - time stepping object
285: Output Parameter:
286: . type - type of scheme being used
288: Level: beginner
290: .seealso: TSSSP, TSSSPSettype(), TSSSPSetNumStages(), TSSSPRKS2, TSSSPRKS3, TSSSPRK104
291: @*/
292: PetscErrorCode TSSSPGetType(TS ts,TSSSPType *type)
293: {
298: PetscUseMethod(ts,"TSSSPGetType_C",(TS,TSSSPType*),(ts,type));
299: return(0);
300: }
302: /*@
303: TSSSPSetNumStages - set the number of stages to use with the SSP method
305: Logically Collective
307: Input Parameters:
308: + ts - time stepping object
309: - nstages - number of stages
311: Options Database Keys:
312: -ts_ssp_type <rks2>: NumStages of SSP method (one of) rks2 rks3 rk104
313: -ts_ssp_nstages <5>: Number of stages
315: Level: beginner
317: .seealso: TSSSP, TSSSPGetNumStages(), TSSSPSetNumStages(), TSSSPRKS2, TSSSPRKS3, TSSSPRK104
318: @*/
319: PetscErrorCode TSSSPSetNumStages(TS ts,PetscInt nstages)
320: {
325: PetscTryMethod(ts,"TSSSPSetNumStages_C",(TS,PetscInt),(ts,nstages));
326: return(0);
327: }
329: /*@
330: TSSSPGetNumStages - get the number of stages in the SSP time integration scheme
332: Logically Collective
334: Input Parameter:
335: . ts - time stepping object
337: Output Parameter:
338: . nstages - number of stages
340: Level: beginner
342: .seealso: TSSSP, TSSSPGetType(), TSSSPSetNumStages(), TSSSPRKS2, TSSSPRKS3, TSSSPRK104
343: @*/
344: PetscErrorCode TSSSPGetNumStages(TS ts,PetscInt *nstages)
345: {
350: PetscUseMethod(ts,"TSSSPGetNumStages_C",(TS,PetscInt*),(ts,nstages));
351: return(0);
352: }
354: static PetscErrorCode TSSSPSetType_SSP(TS ts,TSSSPType type)
355: {
356: PetscErrorCode ierr,(*r)(TS,PetscReal,PetscReal,Vec);
357: TS_SSP *ssp = (TS_SSP*)ts->data;
360: PetscFunctionListFind(TSSSPList,type,&r);
361: if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown TS_SSP type %s given",type);
362: ssp->onestep = r;
363: PetscFree(ssp->type_name);
364: PetscStrallocpy(type,&ssp->type_name);
365: ts->default_adapt_type = TSADAPTNONE;
366: return(0);
367: }
368: static PetscErrorCode TSSSPGetType_SSP(TS ts,TSSSPType *type)
369: {
370: TS_SSP *ssp = (TS_SSP*)ts->data;
373: *type = ssp->type_name;
374: return(0);
375: }
376: static PetscErrorCode TSSSPSetNumStages_SSP(TS ts,PetscInt nstages)
377: {
378: TS_SSP *ssp = (TS_SSP*)ts->data;
381: ssp->nstages = nstages;
382: return(0);
383: }
384: static PetscErrorCode TSSSPGetNumStages_SSP(TS ts,PetscInt *nstages)
385: {
386: TS_SSP *ssp = (TS_SSP*)ts->data;
389: *nstages = ssp->nstages;
390: return(0);
391: }
393: static PetscErrorCode TSSetFromOptions_SSP(PetscOptionItems *PetscOptionsObject,TS ts)
394: {
395: char tname[256] = TSSSPRKS2;
396: TS_SSP *ssp = (TS_SSP*)ts->data;
398: PetscBool flg;
401: PetscOptionsHead(PetscOptionsObject,"SSP ODE solver options");
402: {
403: PetscOptionsFList("-ts_ssp_type","Type of SSP method","TSSSPSetType",TSSSPList,tname,tname,sizeof(tname),&flg);
404: if (flg) {
405: TSSSPSetType(ts,tname);
406: }
407: PetscOptionsInt("-ts_ssp_nstages","Number of stages","TSSSPSetNumStages",ssp->nstages,&ssp->nstages,NULL);
408: }
409: PetscOptionsTail();
410: return(0);
411: }
413: static PetscErrorCode TSView_SSP(TS ts,PetscViewer viewer)
414: {
415: TS_SSP *ssp = (TS_SSP*)ts->data;
416: PetscBool ascii;
420: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&ascii);
421: if (ascii) {PetscViewerASCIIPrintf(viewer," Scheme: %s\n",ssp->type_name);}
422: return(0);
423: }
425: /* ------------------------------------------------------------ */
427: /*MC
428: TSSSP - Explicit strong stability preserving ODE solver
430: Most hyperbolic conservation laws have exact solutions that are total variation diminishing (TVD) or total variation
431: bounded (TVB) although these solutions often contain discontinuities. Spatial discretizations such as Godunov's
432: scheme and high-resolution finite volume methods (TVD limiters, ENO/WENO) are designed to preserve these properties,
433: but they are usually formulated using a forward Euler time discretization or by coupling the space and time
434: discretization as in the classical Lax-Wendroff scheme. When the space and time discretization is coupled, it is very
435: difficult to produce schemes with high temporal accuracy while preserving TVD properties. An alternative is the
436: semidiscrete formulation where we choose a spatial discretization that is TVD with forward Euler and then choose a
437: time discretization that preserves the TVD property. Such integrators are called strong stability preserving (SSP).
439: Let c_eff be the minimum number of function evaluations required to step as far as one step of forward Euler while
440: still being SSP. Some theoretical bounds
442: 1. There are no explicit methods with c_eff > 1.
444: 2. There are no explicit methods beyond order 4 (for nonlinear problems) and c_eff > 0.
446: 3. There are no implicit methods with order greater than 1 and c_eff > 2.
448: This integrator provides Runge-Kutta methods of order 2, 3, and 4 with maximal values of c_eff. More stages allows
449: for larger values of c_eff which improves efficiency. These implementations are low-memory and only use 2 or 3 work
450: vectors regardless of the total number of stages, so e.g. 25-stage 3rd order methods may be an excellent choice.
452: Methods can be chosen with -ts_ssp_type {rks2,rks3,rk104}
454: rks2: Second order methods with any number s>1 of stages. c_eff = (s-1)/s
456: rks3: Third order methods with s=n^2 stages, n>1. c_eff = (s-n)/s
458: rk104: A 10-stage fourth order method. c_eff = 0.6
460: Level: beginner
462: References:
463: + 1. - Ketcheson, Highly efficient strong stability preserving Runge Kutta methods with low storage implementations, SISC, 2008.
464: - 2. - Gottlieb, Ketcheson, and Shu, High order strong stability preserving time discretizations, J Scientific Computing, 2009.
466: .seealso: TSCreate(), TS, TSSetType()
468: M*/
469: PETSC_EXTERN PetscErrorCode TSCreate_SSP(TS ts)
470: {
471: TS_SSP *ssp;
475: TSSSPInitializePackage();
477: ts->ops->setup = TSSetUp_SSP;
478: ts->ops->step = TSStep_SSP;
479: ts->ops->reset = TSReset_SSP;
480: ts->ops->destroy = TSDestroy_SSP;
481: ts->ops->setfromoptions = TSSetFromOptions_SSP;
482: ts->ops->view = TSView_SSP;
484: PetscNewLog(ts,&ssp);
485: ts->data = (void*)ssp;
487: PetscObjectComposeFunction((PetscObject)ts,"TSSSPGetType_C",TSSSPGetType_SSP);
488: PetscObjectComposeFunction((PetscObject)ts,"TSSSPSetType_C",TSSSPSetType_SSP);
489: PetscObjectComposeFunction((PetscObject)ts,"TSSSPGetNumStages_C",TSSSPGetNumStages_SSP);
490: PetscObjectComposeFunction((PetscObject)ts,"TSSSPSetNumStages_C",TSSSPSetNumStages_SSP);
492: TSSSPSetType(ts,TSSSPRKS2);
493: ssp->nstages = 5;
494: return(0);
495: }
497: /*@C
498: TSSSPInitializePackage - This function initializes everything in the TSSSP package. It is called
499: from TSInitializePackage().
501: Level: developer
503: .seealso: PetscInitialize()
504: @*/
505: PetscErrorCode TSSSPInitializePackage(void)
506: {
510: if (TSSSPPackageInitialized) return(0);
511: TSSSPPackageInitialized = PETSC_TRUE;
512: PetscFunctionListAdd(&TSSSPList,TSSSPRKS2, TSSSPStep_RK_2);
513: PetscFunctionListAdd(&TSSSPList,TSSSPRKS3, TSSSPStep_RK_3);
514: PetscFunctionListAdd(&TSSSPList,TSSSPRK104,TSSSPStep_RK_10_4);
515: PetscRegisterFinalize(TSSSPFinalizePackage);
516: return(0);
517: }
519: /*@C
520: TSSSPFinalizePackage - This function destroys everything in the TSSSP package. It is
521: called from PetscFinalize().
523: Level: developer
525: .seealso: PetscFinalize()
526: @*/
527: PetscErrorCode TSSSPFinalizePackage(void)
528: {
532: TSSSPPackageInitialized = PETSC_FALSE;
533: PetscFunctionListDestroy(&TSSSPList);
534: return(0);
535: }