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) VecDestroyVecs(ssp->nwork, &ssp->work);
25: VecDuplicateVecs(ts->vec_sol, n, &ssp->work);
26: ssp->nwork = n;
27: }
28: *work = ssp->work;
29: ssp->workout = PETSC_TRUE;
30: return 0;
31: }
33: static PetscErrorCode TSSSPRestoreWorkVectors(TS ts, PetscInt n, Vec **work)
34: {
35: TS_SSP *ssp = (TS_SSP *)ts->data;
39: ssp->workout = PETSC_FALSE;
40: *work = NULL;
41: return 0;
42: }
44: /*MC
45: TSSSPRKS2 - Optimal second order SSP Runge-Kutta method, low-storage, c_eff=(s-1)/s
47: Pseudocode 2 of Ketcheson 2008
49: Level: beginner
51: .seealso: [](chapter_ts), `TSSSP`, `TSSSPSetType()`, `TSSSPSetNumStages()`
52: M*/
53: static PetscErrorCode TSSSPStep_RK_2(TS ts, PetscReal t0, PetscReal dt, Vec sol)
54: {
55: TS_SSP *ssp = (TS_SSP *)ts->data;
56: Vec *work, F;
57: PetscInt i, s;
59: s = ssp->nstages;
60: TSSSPGetWorkVectors(ts, 2, &work);
61: F = work[1];
62: VecCopy(sol, work[0]);
63: for (i = 0; i < s - 1; i++) {
64: PetscReal stage_time = t0 + dt * (i / (s - 1.));
65: TSPreStage(ts, stage_time);
66: TSComputeRHSFunction(ts, stage_time, work[0], F);
67: VecAXPY(work[0], dt / (s - 1.), F);
68: }
69: TSComputeRHSFunction(ts, t0 + dt, work[0], F);
70: VecAXPBYPCZ(sol, (s - 1.) / s, dt / s, 1. / s, work[0], F);
71: TSSSPRestoreWorkVectors(ts, 2, &work);
72: return 0;
73: }
75: /*MC
76: TSSSPRKS3 - Optimal third order SSP Runge-Kutta, low-storage, c_eff=(PetscSqrtReal(s)-1)/PetscSqrtReal(s), where PetscSqrtReal(s) is an integer
78: Pseudocode 2 of Ketcheson 2008
80: Level: beginner
82: .seealso: [](chapter_ts), `TSSSP`, `TSSSPSetType()`, `TSSSPSetNumStages()`
83: M*/
84: static PetscErrorCode TSSSPStep_RK_3(TS ts, PetscReal t0, PetscReal dt, Vec sol)
85: {
86: TS_SSP *ssp = (TS_SSP *)ts->data;
87: Vec *work, F;
88: PetscInt i, s, n, r;
89: PetscReal c, stage_time;
91: s = ssp->nstages;
92: n = (PetscInt)(PetscSqrtReal((PetscReal)s) + 0.001);
93: r = s - n;
95: TSSSPGetWorkVectors(ts, 3, &work);
96: F = work[2];
97: VecCopy(sol, work[0]);
98: for (i = 0; i < (n - 1) * (n - 2) / 2; i++) {
99: c = (i < n * (n + 1) / 2) ? 1. * i / (s - n) : (1. * i - n) / (s - n);
100: stage_time = t0 + c * dt;
101: TSPreStage(ts, stage_time);
102: TSComputeRHSFunction(ts, stage_time, work[0], F);
103: VecAXPY(work[0], dt / r, F);
104: }
105: VecCopy(work[0], work[1]);
106: for (; i < n * (n + 1) / 2 - 1; i++) {
107: c = (i < n * (n + 1) / 2) ? 1. * i / (s - n) : (1. * i - n) / (s - n);
108: stage_time = t0 + c * dt;
109: TSPreStage(ts, stage_time);
110: TSComputeRHSFunction(ts, stage_time, work[0], F);
111: VecAXPY(work[0], dt / r, F);
112: }
113: {
114: c = (i < n * (n + 1) / 2) ? 1. * i / (s - n) : (1. * i - n) / (s - n);
115: stage_time = t0 + c * dt;
116: TSPreStage(ts, stage_time);
117: TSComputeRHSFunction(ts, stage_time, work[0], F);
118: VecAXPBYPCZ(work[0], 1. * n / (2 * n - 1.), (n - 1.) * dt / (r * (2 * n - 1)), (n - 1.) / (2 * n - 1.), work[1], F);
119: i++;
120: }
121: for (; i < s; i++) {
122: c = (i < n * (n + 1) / 2) ? 1. * i / (s - n) : (1. * i - n) / (s - n);
123: stage_time = t0 + c * dt;
124: TSPreStage(ts, stage_time);
125: TSComputeRHSFunction(ts, stage_time, work[0], F);
126: VecAXPY(work[0], dt / r, F);
127: }
128: VecCopy(work[0], sol);
129: TSSSPRestoreWorkVectors(ts, 3, &work);
130: return 0;
131: }
133: /*MC
134: TSSSPRKS104 - Optimal fourth order SSP Runge-Kutta, low-storage (2N), c_eff=0.6
136: SSPRK(10,4), Pseudocode 3 of Ketcheson 2008
138: Level: beginner
140: .seealso: [](chapter_ts), `TSSSP`, `TSSSPSetType()`
141: M*/
142: static PetscErrorCode TSSSPStep_RK_10_4(TS ts, PetscReal t0, PetscReal dt, Vec sol)
143: {
144: const PetscReal c[10] = {0, 1. / 6, 2. / 6, 3. / 6, 4. / 6, 2. / 6, 3. / 6, 4. / 6, 5. / 6, 1};
145: Vec *work, F;
146: PetscInt i;
147: PetscReal stage_time;
149: TSSSPGetWorkVectors(ts, 3, &work);
150: F = work[2];
151: VecCopy(sol, work[0]);
152: for (i = 0; i < 5; i++) {
153: stage_time = t0 + c[i] * dt;
154: TSPreStage(ts, stage_time);
155: TSComputeRHSFunction(ts, stage_time, work[0], F);
156: VecAXPY(work[0], dt / 6, F);
157: }
158: VecAXPBYPCZ(work[1], 1. / 25, 9. / 25, 0, sol, work[0]);
159: VecAXPBY(work[0], 15, -5, work[1]);
160: for (; i < 9; i++) {
161: stage_time = t0 + c[i] * dt;
162: TSPreStage(ts, stage_time);
163: TSComputeRHSFunction(ts, stage_time, work[0], F);
164: VecAXPY(work[0], dt / 6, F);
165: }
166: stage_time = t0 + dt;
167: TSPreStage(ts, stage_time);
168: TSComputeRHSFunction(ts, stage_time, work[0], F);
169: VecAXPBYPCZ(work[1], 3. / 5, dt / 10, 1, work[0], F);
170: VecCopy(work[1], sol);
171: TSSSPRestoreWorkVectors(ts, 3, &work);
172: return 0;
173: }
175: static PetscErrorCode TSSetUp_SSP(TS ts)
176: {
177: TSCheckImplicitTerm(ts);
178: TSGetAdapt(ts, &ts->adapt);
179: TSAdaptCandidatesClear(ts->adapt);
180: return 0;
181: }
183: static PetscErrorCode TSStep_SSP(TS ts)
184: {
185: TS_SSP *ssp = (TS_SSP *)ts->data;
186: Vec sol = ts->vec_sol;
187: PetscBool stageok, accept = PETSC_TRUE;
188: PetscReal next_time_step = ts->time_step;
190: (*ssp->onestep)(ts, ts->ptime, ts->time_step, sol);
191: TSPostStage(ts, ts->ptime, 0, &sol);
192: TSAdaptCheckStage(ts->adapt, ts, ts->ptime + ts->time_step, sol, &stageok);
193: if (!stageok) {
194: ts->reason = TS_DIVERGED_STEP_REJECTED;
195: return 0;
196: }
198: TSAdaptChoose(ts->adapt, ts, ts->time_step, NULL, &next_time_step, &accept);
199: if (!accept) {
200: ts->reason = TS_DIVERGED_STEP_REJECTED;
201: return 0;
202: }
204: ts->ptime += ts->time_step;
205: ts->time_step = next_time_step;
206: return 0;
207: }
208: /*------------------------------------------------------------*/
210: static PetscErrorCode TSReset_SSP(TS ts)
211: {
212: TS_SSP *ssp = (TS_SSP *)ts->data;
214: if (ssp->work) VecDestroyVecs(ssp->nwork, &ssp->work);
215: ssp->nwork = 0;
216: ssp->workout = PETSC_FALSE;
217: return 0;
218: }
220: static PetscErrorCode TSDestroy_SSP(TS ts)
221: {
222: TS_SSP *ssp = (TS_SSP *)ts->data;
224: TSReset_SSP(ts);
225: PetscFree(ssp->type_name);
226: PetscFree(ts->data);
227: PetscObjectComposeFunction((PetscObject)ts, "TSSSPGetType_C", NULL);
228: PetscObjectComposeFunction((PetscObject)ts, "TSSSPSetType_C", NULL);
229: PetscObjectComposeFunction((PetscObject)ts, "TSSSPGetNumStages_C", NULL);
230: PetscObjectComposeFunction((PetscObject)ts, "TSSSPSetNumStages_C", NULL);
231: return 0;
232: }
233: /*------------------------------------------------------------*/
235: /*@C
236: TSSSPSetType - set the `TSSSP` time integration scheme to use
238: Logically Collective
240: Input Parameters:
241: + ts - time stepping object
242: - ssptype - type of scheme to use
244: Options Database Keys:
245: + -ts_ssp_type <rks2> - Type of `TSSSP` method (one of) rks2 rks3 rk104
246: - -ts_ssp_nstages<rks2: 5, rks3: 9> - Number of stages
248: Level: beginner
250: .seealso: [](chapter_ts), `TSSSP`, `TSSSPGetType()`, `TSSSPSetNumStages()`, `TSSSPRKS2`, `TSSSPRKS3`, `TSSSPRK104`
251: @*/
252: PetscErrorCode TSSSPSetType(TS ts, TSSSPType ssptype)
253: {
256: PetscTryMethod(ts, "TSSSPSetType_C", (TS, TSSSPType), (ts, ssptype));
257: return 0;
258: }
260: /*@C
261: TSSSPGetType - get the `TSSSP` time integration scheme
263: Logically Collective
265: Input Parameter:
266: . ts - time stepping object
268: Output Parameter:
269: . type - type of scheme being used
271: Level: beginner
273: .seealso: [](chapter_ts), `TSSSP`, `TSSSPSettype()`, `TSSSPSetNumStages()`, `TSSSPRKS2`, `TSSSPRKS3`, `TSSSPRK104`
274: @*/
275: PetscErrorCode TSSSPGetType(TS ts, TSSSPType *type)
276: {
278: PetscUseMethod(ts, "TSSSPGetType_C", (TS, TSSSPType *), (ts, type));
279: return 0;
280: }
282: /*@
283: TSSSPSetNumStages - set the number of stages to use with the `TSSSP` method. Must be called after
284: `TSSSPSetType()`.
286: Logically Collective
288: Input Parameters:
289: + ts - time stepping object
290: - nstages - number of stages
292: Options Database Keys:
293: + -ts_ssp_type <rks2> - Type of `TSSSP` method (one of) rks2 rks3 rk104
294: - -ts_ssp_nstages<rks2: 5, rks3: 9> - Number of stages
296: Level: beginner
298: .seealso: [](chapter_ts), `TSSSP`, `TSSSPGetNumStages()`, `TSSSPSetNumStages()`, `TSSSPRKS2`, `TSSSPRKS3`, `TSSSPRK104`
299: @*/
300: PetscErrorCode TSSSPSetNumStages(TS ts, PetscInt nstages)
301: {
303: PetscTryMethod(ts, "TSSSPSetNumStages_C", (TS, PetscInt), (ts, nstages));
304: return 0;
305: }
307: /*@
308: TSSSPGetNumStages - get the number of stages in the `TSSSP` time integration scheme
310: Logically Collective
312: Input Parameter:
313: . ts - time stepping object
315: Output Parameter:
316: . nstages - number of stages
318: Level: beginner
320: .seealso: [](chapter_ts), `TSSSP`, `TSSSPGetType()`, `TSSSPSetNumStages()`, `TSSSPRKS2`, `TSSSPRKS3`, `TSSSPRK104`
321: @*/
322: PetscErrorCode TSSSPGetNumStages(TS ts, PetscInt *nstages)
323: {
325: PetscUseMethod(ts, "TSSSPGetNumStages_C", (TS, PetscInt *), (ts, nstages));
326: return 0;
327: }
329: static PetscErrorCode TSSSPSetType_SSP(TS ts, TSSSPType type)
330: {
331: TS_SSP *ssp = (TS_SSP *)ts->data;
332: PetscErrorCode (*r)(TS, PetscReal, PetscReal, Vec);
333: PetscBool flag;
335: PetscFunctionListFind(TSSSPList, type, &r);
337: ssp->onestep = r;
338: PetscFree(ssp->type_name);
339: PetscStrallocpy(type, &ssp->type_name);
340: ts->default_adapt_type = TSADAPTNONE;
341: PetscStrcmp(type, TSSSPRKS2, &flag);
342: if (flag) {
343: ssp->nstages = 5;
344: } else {
345: PetscStrcmp(type, TSSSPRKS3, &flag);
346: if (flag) {
347: ssp->nstages = 9;
348: } else {
349: ssp->nstages = 5;
350: }
351: }
352: return 0;
353: }
354: static PetscErrorCode TSSSPGetType_SSP(TS ts, TSSSPType *type)
355: {
356: TS_SSP *ssp = (TS_SSP *)ts->data;
358: *type = ssp->type_name;
359: return 0;
360: }
361: static PetscErrorCode TSSSPSetNumStages_SSP(TS ts, PetscInt nstages)
362: {
363: TS_SSP *ssp = (TS_SSP *)ts->data;
365: ssp->nstages = nstages;
366: return 0;
367: }
368: static PetscErrorCode TSSSPGetNumStages_SSP(TS ts, PetscInt *nstages)
369: {
370: TS_SSP *ssp = (TS_SSP *)ts->data;
372: *nstages = ssp->nstages;
373: return 0;
374: }
376: static PetscErrorCode TSSetFromOptions_SSP(TS ts, PetscOptionItems *PetscOptionsObject)
377: {
378: char tname[256] = TSSSPRKS2;
379: TS_SSP *ssp = (TS_SSP *)ts->data;
380: PetscBool flg;
382: PetscOptionsHeadBegin(PetscOptionsObject, "SSP ODE solver options");
383: {
384: PetscOptionsFList("-ts_ssp_type", "Type of SSP method", "TSSSPSetType", TSSSPList, tname, tname, sizeof(tname), &flg);
385: if (flg) TSSSPSetType(ts, tname);
386: PetscOptionsInt("-ts_ssp_nstages", "Number of stages", "TSSSPSetNumStages", ssp->nstages, &ssp->nstages, NULL);
387: }
388: PetscOptionsHeadEnd();
389: return 0;
390: }
392: static PetscErrorCode TSView_SSP(TS ts, PetscViewer viewer)
393: {
394: TS_SSP *ssp = (TS_SSP *)ts->data;
395: PetscBool ascii;
397: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &ascii);
398: if (ascii) PetscViewerASCIIPrintf(viewer, " Scheme: %s\n", ssp->type_name);
399: return 0;
400: }
402: /* ------------------------------------------------------------ */
404: /*MC
405: TSSSP - Explicit strong stability preserving ODE solver
407: Most hyperbolic conservation laws have exact solutions that are total variation diminishing (TVD) or total variation
408: bounded (TVB) although these solutions often contain discontinuities. Spatial discretizations such as Godunov's
409: scheme and high-resolution finite volume methods (TVD limiters, ENO/WENO) are designed to preserve these properties,
410: but they are usually formulated using a forward Euler time discretization or by coupling the space and time
411: discretization as in the classical Lax-Wendroff scheme. When the space and time discretization is coupled, it is very
412: difficult to produce schemes with high temporal accuracy while preserving TVD properties. An alternative is the
413: semidiscrete formulation where we choose a spatial discretization that is TVD with forward Euler and then choose a
414: time discretization that preserves the TVD property. Such integrators are called strong stability preserving (SSP).
416: Let c_eff be the minimum number of function evaluations required to step as far as one step of forward Euler while
417: still being SSP. Some theoretical bounds
419: 1. There are no explicit methods with c_eff > 1.
421: 2. There are no explicit methods beyond order 4 (for nonlinear problems) and c_eff > 0.
423: 3. There are no implicit methods with order greater than 1 and c_eff > 2.
425: This integrator provides Runge-Kutta methods of order 2, 3, and 4 with maximal values of c_eff. More stages allows
426: for larger values of c_eff which improves efficiency. These implementations are low-memory and only use 2 or 3 work
427: vectors regardless of the total number of stages, so e.g. 25-stage 3rd order methods may be an excellent choice.
429: Methods can be chosen with -ts_ssp_type {rks2,rks3,rk104}
431: rks2: Second order methods with any number s>1 of stages. c_eff = (s-1)/s
433: rks3: Third order methods with s=n^2 stages, n>1. c_eff = (s-n)/s
435: rk104: A 10-stage fourth order method. c_eff = 0.6
437: Level: beginner
439: References:
440: + * - Ketcheson, Highly efficient strong stability preserving Runge Kutta methods with low storage implementations, SISC, 2008.
441: - * - Gottlieb, Ketcheson, and Shu, High order strong stability preserving time discretizations, J Scientific Computing, 2009.
443: .seealso: [](chapter_ts), `TSCreate()`, `TS`, `TSSetType()`
445: M*/
446: PETSC_EXTERN PetscErrorCode TSCreate_SSP(TS ts)
447: {
448: TS_SSP *ssp;
450: TSSSPInitializePackage();
452: ts->ops->setup = TSSetUp_SSP;
453: ts->ops->step = TSStep_SSP;
454: ts->ops->reset = TSReset_SSP;
455: ts->ops->destroy = TSDestroy_SSP;
456: ts->ops->setfromoptions = TSSetFromOptions_SSP;
457: ts->ops->view = TSView_SSP;
459: PetscNew(&ssp);
460: ts->data = (void *)ssp;
462: PetscObjectComposeFunction((PetscObject)ts, "TSSSPGetType_C", TSSSPGetType_SSP);
463: PetscObjectComposeFunction((PetscObject)ts, "TSSSPSetType_C", TSSSPSetType_SSP);
464: PetscObjectComposeFunction((PetscObject)ts, "TSSSPGetNumStages_C", TSSSPGetNumStages_SSP);
465: PetscObjectComposeFunction((PetscObject)ts, "TSSSPSetNumStages_C", TSSSPSetNumStages_SSP);
467: TSSSPSetType(ts, TSSSPRKS2);
468: return 0;
469: }
471: /*@C
472: TSSSPInitializePackage - This function initializes everything in the `TSSSP` package. It is called
473: from `TSInitializePackage()`.
475: Level: developer
477: .seealso: [](chapter_ts), `PetscInitialize()`, `TSSSPFinalizePackage()`, `TSInitializePackage()`
478: @*/
479: PetscErrorCode TSSSPInitializePackage(void)
480: {
481: if (TSSSPPackageInitialized) return 0;
482: TSSSPPackageInitialized = PETSC_TRUE;
483: PetscFunctionListAdd(&TSSSPList, TSSSPRKS2, TSSSPStep_RK_2);
484: PetscFunctionListAdd(&TSSSPList, TSSSPRKS3, TSSSPStep_RK_3);
485: PetscFunctionListAdd(&TSSSPList, TSSSPRK104, TSSSPStep_RK_10_4);
486: PetscRegisterFinalize(TSSSPFinalizePackage);
487: return 0;
488: }
490: /*@C
491: TSSSPFinalizePackage - This function destroys everything in the `TSSSP` package. It is
492: called from `PetscFinalize()`.
494: Level: developer
496: .seealso: [](chapter_ts), `PetscFinalize()`, `TSSSPInitiallizePackage()`, `TSInitializePackage()`
497: @*/
498: PetscErrorCode TSSSPFinalizePackage(void)
499: {
500: TSSSPPackageInitialized = PETSC_FALSE;
501: PetscFunctionListDestroy(&TSSSPList);
502: return 0;
503: }