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;
22: PetscFunctionBegin;
23: PetscCheck(!ssp->workout, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Work vectors already gotten");
24: if (ssp->nwork < n) {
25: if (ssp->nwork > 0) PetscCall(VecDestroyVecs(ssp->nwork, &ssp->work));
26: PetscCall(VecDuplicateVecs(ts->vec_sol, n, &ssp->work));
27: ssp->nwork = n;
28: }
29: *work = ssp->work;
30: ssp->workout = PETSC_TRUE;
31: PetscFunctionReturn(PETSC_SUCCESS);
32: }
34: static PetscErrorCode TSSSPRestoreWorkVectors(TS ts, PetscInt n, Vec **work)
35: {
36: TS_SSP *ssp = (TS_SSP *)ts->data;
38: PetscFunctionBegin;
39: PetscCheck(ssp->workout, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Work vectors have not been gotten");
40: PetscCheck(*work == ssp->work, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Wrong work vectors checked out");
41: ssp->workout = PETSC_FALSE;
42: *work = NULL;
43: PetscFunctionReturn(PETSC_SUCCESS);
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: [](ch_ts), `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: PetscFunctionBegin;
62: s = ssp->nstages;
63: PetscCall(TSSSPGetWorkVectors(ts, 2, &work));
64: F = work[1];
65: PetscCall(VecCopy(sol, work[0]));
66: for (i = 0; i < s - 1; i++) {
67: PetscReal stage_time = t0 + dt * (i / (s - 1.));
68: PetscCall(TSPreStage(ts, stage_time));
69: PetscCall(TSComputeRHSFunction(ts, stage_time, work[0], F));
70: PetscCall(VecAXPY(work[0], dt / (s - 1.), F));
71: }
72: PetscCall(TSComputeRHSFunction(ts, t0 + dt, work[0], F));
73: PetscCall(VecAXPBYPCZ(sol, (s - 1.) / s, dt / s, 1. / s, work[0], F));
74: PetscCall(TSSSPRestoreWorkVectors(ts, 2, &work));
75: PetscFunctionReturn(PETSC_SUCCESS);
76: }
78: /*MC
79: TSSSPRKS3 - Optimal third order SSP Runge-Kutta, low-storage, c_eff=(PetscSqrtReal(s)-1)/PetscSqrtReal(s), where PetscSqrtReal(s) is an integer
81: Pseudocode 2 of Ketcheson 2008
83: Level: beginner
85: .seealso: [](ch_ts), `TSSSP`, `TSSSPSetType()`, `TSSSPSetNumStages()`
86: M*/
87: static PetscErrorCode TSSSPStep_RK_3(TS ts, PetscReal t0, PetscReal dt, Vec sol)
88: {
89: TS_SSP *ssp = (TS_SSP *)ts->data;
90: Vec *work, F;
91: PetscInt i, s, n, r;
92: PetscReal c, stage_time;
94: PetscFunctionBegin;
95: s = ssp->nstages;
96: n = (PetscInt)(PetscSqrtReal((PetscReal)s) + 0.001);
97: r = s - n;
98: PetscCheck(n * n == s, PETSC_COMM_SELF, PETSC_ERR_SUP, "No support for optimal third order schemes with %" PetscInt_FMT " stages, must be a square number at least 4", s);
99: PetscCall(TSSSPGetWorkVectors(ts, 3, &work));
100: F = work[2];
101: PetscCall(VecCopy(sol, work[0]));
102: for (i = 0; i < (n - 1) * (n - 2) / 2; i++) {
103: c = (i < n * (n + 1) / 2) ? 1. * i / (s - n) : (1. * i - n) / (s - n);
104: stage_time = t0 + c * dt;
105: PetscCall(TSPreStage(ts, stage_time));
106: PetscCall(TSComputeRHSFunction(ts, stage_time, work[0], F));
107: PetscCall(VecAXPY(work[0], dt / r, F));
108: }
109: PetscCall(VecCopy(work[0], work[1]));
110: for (; i < n * (n + 1) / 2 - 1; i++) {
111: c = (i < n * (n + 1) / 2) ? 1. * i / (s - n) : (1. * i - n) / (s - n);
112: stage_time = t0 + c * dt;
113: PetscCall(TSPreStage(ts, stage_time));
114: PetscCall(TSComputeRHSFunction(ts, stage_time, work[0], F));
115: PetscCall(VecAXPY(work[0], dt / r, F));
116: }
117: {
118: c = (i < n * (n + 1) / 2) ? 1. * i / (s - n) : (1. * i - n) / (s - n);
119: stage_time = t0 + c * dt;
120: PetscCall(TSPreStage(ts, stage_time));
121: PetscCall(TSComputeRHSFunction(ts, stage_time, work[0], F));
122: PetscCall(VecAXPBYPCZ(work[0], 1. * n / (2 * n - 1.), (n - 1.) * dt / (r * (2 * n - 1)), (n - 1.) / (2 * n - 1.), work[1], F));
123: i++;
124: }
125: for (; i < s; i++) {
126: c = (i < n * (n + 1) / 2) ? 1. * i / (s - n) : (1. * i - n) / (s - n);
127: stage_time = t0 + c * dt;
128: PetscCall(TSPreStage(ts, stage_time));
129: PetscCall(TSComputeRHSFunction(ts, stage_time, work[0], F));
130: PetscCall(VecAXPY(work[0], dt / r, F));
131: }
132: PetscCall(VecCopy(work[0], sol));
133: PetscCall(TSSSPRestoreWorkVectors(ts, 3, &work));
134: PetscFunctionReturn(PETSC_SUCCESS);
135: }
137: /*MC
138: TSSSPRKS104 - Optimal fourth order SSP Runge-Kutta, low-storage (2N), c_eff=0.6
140: SSPRK(10,4), Pseudocode 3 of Ketcheson 2008
142: Level: beginner
144: .seealso: [](ch_ts), `TSSSP`, `TSSSPSetType()`
145: M*/
146: static PetscErrorCode TSSSPStep_RK_10_4(TS ts, PetscReal t0, PetscReal dt, Vec sol)
147: {
148: const PetscReal c[10] = {0, 1. / 6, 2. / 6, 3. / 6, 4. / 6, 2. / 6, 3. / 6, 4. / 6, 5. / 6, 1};
149: Vec *work, F;
150: PetscInt i;
151: PetscReal stage_time;
153: PetscFunctionBegin;
154: PetscCall(TSSSPGetWorkVectors(ts, 3, &work));
155: F = work[2];
156: PetscCall(VecCopy(sol, work[0]));
157: for (i = 0; i < 5; i++) {
158: stage_time = t0 + c[i] * dt;
159: PetscCall(TSPreStage(ts, stage_time));
160: PetscCall(TSComputeRHSFunction(ts, stage_time, work[0], F));
161: PetscCall(VecAXPY(work[0], dt / 6, F));
162: }
163: PetscCall(VecAXPBYPCZ(work[1], 1. / 25, 9. / 25, 0, sol, work[0]));
164: PetscCall(VecAXPBY(work[0], 15, -5, work[1]));
165: for (; i < 9; i++) {
166: stage_time = t0 + c[i] * dt;
167: PetscCall(TSPreStage(ts, stage_time));
168: PetscCall(TSComputeRHSFunction(ts, stage_time, work[0], F));
169: PetscCall(VecAXPY(work[0], dt / 6, F));
170: }
171: stage_time = t0 + dt;
172: PetscCall(TSPreStage(ts, stage_time));
173: PetscCall(TSComputeRHSFunction(ts, stage_time, work[0], F));
174: PetscCall(VecAXPBYPCZ(work[1], 3. / 5, dt / 10, 1, work[0], F));
175: PetscCall(VecCopy(work[1], sol));
176: PetscCall(TSSSPRestoreWorkVectors(ts, 3, &work));
177: PetscFunctionReturn(PETSC_SUCCESS);
178: }
180: static PetscErrorCode TSSetUp_SSP(TS ts)
181: {
182: PetscFunctionBegin;
183: PetscCall(TSCheckImplicitTerm(ts));
184: PetscCall(TSGetAdapt(ts, &ts->adapt));
185: PetscCall(TSAdaptCandidatesClear(ts->adapt));
186: PetscFunctionReturn(PETSC_SUCCESS);
187: }
189: static PetscErrorCode TSStep_SSP(TS ts)
190: {
191: TS_SSP *ssp = (TS_SSP *)ts->data;
192: Vec sol = ts->vec_sol;
193: PetscBool stageok, accept = PETSC_TRUE;
194: PetscReal next_time_step = ts->time_step;
196: PetscFunctionBegin;
197: PetscCall((*ssp->onestep)(ts, ts->ptime, ts->time_step, sol));
198: PetscCall(TSPostStage(ts, ts->ptime, 0, &sol));
199: PetscCall(TSAdaptCheckStage(ts->adapt, ts, ts->ptime + ts->time_step, sol, &stageok));
200: if (!stageok) {
201: ts->reason = TS_DIVERGED_STEP_REJECTED;
202: PetscFunctionReturn(PETSC_SUCCESS);
203: }
205: PetscCall(TSAdaptChoose(ts->adapt, ts, ts->time_step, NULL, &next_time_step, &accept));
206: if (!accept) {
207: ts->reason = TS_DIVERGED_STEP_REJECTED;
208: PetscFunctionReturn(PETSC_SUCCESS);
209: }
211: ts->ptime += ts->time_step;
212: ts->time_step = next_time_step;
213: PetscFunctionReturn(PETSC_SUCCESS);
214: }
215: /*------------------------------------------------------------*/
217: static PetscErrorCode TSReset_SSP(TS ts)
218: {
219: TS_SSP *ssp = (TS_SSP *)ts->data;
221: PetscFunctionBegin;
222: if (ssp->work) PetscCall(VecDestroyVecs(ssp->nwork, &ssp->work));
223: ssp->nwork = 0;
224: ssp->workout = PETSC_FALSE;
225: PetscFunctionReturn(PETSC_SUCCESS);
226: }
228: static PetscErrorCode TSDestroy_SSP(TS ts)
229: {
230: TS_SSP *ssp = (TS_SSP *)ts->data;
232: PetscFunctionBegin;
233: PetscCall(TSReset_SSP(ts));
234: PetscCall(PetscFree(ssp->type_name));
235: PetscCall(PetscFree(ts->data));
236: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSSSPGetType_C", NULL));
237: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSSSPSetType_C", NULL));
238: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSSSPGetNumStages_C", NULL));
239: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSSSPSetNumStages_C", NULL));
240: PetscFunctionReturn(PETSC_SUCCESS);
241: }
242: /*------------------------------------------------------------*/
244: /*@C
245: TSSSPSetType - set the `TSSSP` time integration scheme to use
247: Logically Collective
249: Input Parameters:
250: + ts - time stepping object
251: - ssptype - type of scheme to use
253: Options Database Keys:
254: + -ts_ssp_type <rks2> - Type of `TSSSP` method (one of) rks2 rks3 rk104
255: - -ts_ssp_nstages<rks2: 5, rks3: 9> - Number of stages
257: Level: beginner
259: .seealso: [](ch_ts), `TSSSP`, `TSSSPGetType()`, `TSSSPSetNumStages()`, `TSSSPRKS2`, `TSSSPRKS3`, `TSSSPRK104`
260: @*/
261: PetscErrorCode TSSSPSetType(TS ts, TSSSPType ssptype)
262: {
263: PetscFunctionBegin;
266: PetscTryMethod(ts, "TSSSPSetType_C", (TS, TSSSPType), (ts, ssptype));
267: PetscFunctionReturn(PETSC_SUCCESS);
268: }
270: /*@C
271: TSSSPGetType - get the `TSSSP` time integration scheme
273: Logically Collective
275: Input Parameter:
276: . ts - time stepping object
278: Output Parameter:
279: . type - type of scheme being used
281: Level: beginner
283: .seealso: [](ch_ts), `TSSSP`, `TSSSPSettype()`, `TSSSPSetNumStages()`, `TSSSPRKS2`, `TSSSPRKS3`, `TSSSPRK104`
284: @*/
285: PetscErrorCode TSSSPGetType(TS ts, TSSSPType *type)
286: {
287: PetscFunctionBegin;
289: PetscUseMethod(ts, "TSSSPGetType_C", (TS, TSSSPType *), (ts, type));
290: PetscFunctionReturn(PETSC_SUCCESS);
291: }
293: /*@
294: TSSSPSetNumStages - set the number of stages to use with the `TSSSP` method. Must be called after
295: `TSSSPSetType()`.
297: Logically Collective
299: Input Parameters:
300: + ts - time stepping object
301: - nstages - number of stages
303: Options Database Keys:
304: + -ts_ssp_type <rks2> - Type of `TSSSP` method (one of) rks2 rks3 rk104
305: - -ts_ssp_nstages<rks2: 5, rks3: 9> - Number of stages
307: Level: beginner
309: .seealso: [](ch_ts), `TSSSP`, `TSSSPGetNumStages()`, `TSSSPSetNumStages()`, `TSSSPRKS2`, `TSSSPRKS3`, `TSSSPRK104`
310: @*/
311: PetscErrorCode TSSSPSetNumStages(TS ts, PetscInt nstages)
312: {
313: PetscFunctionBegin;
315: PetscTryMethod(ts, "TSSSPSetNumStages_C", (TS, PetscInt), (ts, nstages));
316: PetscFunctionReturn(PETSC_SUCCESS);
317: }
319: /*@
320: TSSSPGetNumStages - get the number of stages in the `TSSSP` time integration scheme
322: Logically Collective
324: Input Parameter:
325: . ts - time stepping object
327: Output Parameter:
328: . nstages - number of stages
330: Level: beginner
332: .seealso: [](ch_ts), `TSSSP`, `TSSSPGetType()`, `TSSSPSetNumStages()`, `TSSSPRKS2`, `TSSSPRKS3`, `TSSSPRK104`
333: @*/
334: PetscErrorCode TSSSPGetNumStages(TS ts, PetscInt *nstages)
335: {
336: PetscFunctionBegin;
338: PetscUseMethod(ts, "TSSSPGetNumStages_C", (TS, PetscInt *), (ts, nstages));
339: PetscFunctionReturn(PETSC_SUCCESS);
340: }
342: static PetscErrorCode TSSSPSetType_SSP(TS ts, TSSSPType type)
343: {
344: TS_SSP *ssp = (TS_SSP *)ts->data;
345: PetscErrorCode (*r)(TS, PetscReal, PetscReal, Vec);
346: PetscBool flag;
348: PetscFunctionBegin;
349: PetscCall(PetscFunctionListFind(TSSSPList, type, &r));
350: PetscCheck(r, PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown TS_SSP type %s given", type);
351: ssp->onestep = r;
352: PetscCall(PetscFree(ssp->type_name));
353: PetscCall(PetscStrallocpy(type, &ssp->type_name));
354: ts->default_adapt_type = TSADAPTNONE;
355: PetscCall(PetscStrcmp(type, TSSSPRKS2, &flag));
356: if (flag) {
357: ssp->nstages = 5;
358: } else {
359: PetscCall(PetscStrcmp(type, TSSSPRKS3, &flag));
360: if (flag) {
361: ssp->nstages = 9;
362: } else {
363: ssp->nstages = 5;
364: }
365: }
366: PetscFunctionReturn(PETSC_SUCCESS);
367: }
368: static PetscErrorCode TSSSPGetType_SSP(TS ts, TSSSPType *type)
369: {
370: TS_SSP *ssp = (TS_SSP *)ts->data;
372: PetscFunctionBegin;
373: *type = ssp->type_name;
374: PetscFunctionReturn(PETSC_SUCCESS);
375: }
376: static PetscErrorCode TSSSPSetNumStages_SSP(TS ts, PetscInt nstages)
377: {
378: TS_SSP *ssp = (TS_SSP *)ts->data;
380: PetscFunctionBegin;
381: ssp->nstages = nstages;
382: PetscFunctionReturn(PETSC_SUCCESS);
383: }
384: static PetscErrorCode TSSSPGetNumStages_SSP(TS ts, PetscInt *nstages)
385: {
386: TS_SSP *ssp = (TS_SSP *)ts->data;
388: PetscFunctionBegin;
389: *nstages = ssp->nstages;
390: PetscFunctionReturn(PETSC_SUCCESS);
391: }
393: static PetscErrorCode TSSetFromOptions_SSP(TS ts, PetscOptionItems *PetscOptionsObject)
394: {
395: char tname[256] = TSSSPRKS2;
396: TS_SSP *ssp = (TS_SSP *)ts->data;
397: PetscBool flg;
399: PetscFunctionBegin;
400: PetscOptionsHeadBegin(PetscOptionsObject, "SSP ODE solver options");
401: {
402: PetscCall(PetscOptionsFList("-ts_ssp_type", "Type of SSP method", "TSSSPSetType", TSSSPList, tname, tname, sizeof(tname), &flg));
403: if (flg) PetscCall(TSSSPSetType(ts, tname));
404: PetscCall(PetscOptionsInt("-ts_ssp_nstages", "Number of stages", "TSSSPSetNumStages", ssp->nstages, &ssp->nstages, NULL));
405: }
406: PetscOptionsHeadEnd();
407: PetscFunctionReturn(PETSC_SUCCESS);
408: }
410: static PetscErrorCode TSView_SSP(TS ts, PetscViewer viewer)
411: {
412: TS_SSP *ssp = (TS_SSP *)ts->data;
413: PetscBool ascii;
415: PetscFunctionBegin;
416: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &ascii));
417: if (ascii) PetscCall(PetscViewerASCIIPrintf(viewer, " Scheme: %s\n", ssp->type_name));
418: PetscFunctionReturn(PETSC_SUCCESS);
419: }
421: /* ------------------------------------------------------------ */
423: /*MC
424: TSSSP - Explicit strong stability preserving ODE solver
426: Most hyperbolic conservation laws have exact solutions that are total variation diminishing (TVD) or total variation
427: bounded (TVB) although these solutions often contain discontinuities. Spatial discretizations such as Godunov's
428: scheme and high-resolution finite volume methods (TVD limiters, ENO/WENO) are designed to preserve these properties,
429: but they are usually formulated using a forward Euler time discretization or by coupling the space and time
430: discretization as in the classical Lax-Wendroff scheme. When the space and time discretization is coupled, it is very
431: difficult to produce schemes with high temporal accuracy while preserving TVD properties. An alternative is the
432: semidiscrete formulation where we choose a spatial discretization that is TVD with forward Euler and then choose a
433: time discretization that preserves the TVD property. Such integrators are called strong stability preserving (SSP).
435: Let c_eff be the minimum number of function evaluations required to step as far as one step of forward Euler while
436: still being SSP. Some theoretical bounds
438: 1. There are no explicit methods with c_eff > 1.
440: 2. There are no explicit methods beyond order 4 (for nonlinear problems) and c_eff > 0.
442: 3. There are no implicit methods with order greater than 1 and c_eff > 2.
444: This integrator provides Runge-Kutta methods of order 2, 3, and 4 with maximal values of c_eff. More stages allows
445: for larger values of c_eff which improves efficiency. These implementations are low-memory and only use 2 or 3 work
446: vectors regardless of the total number of stages, so e.g. 25-stage 3rd order methods may be an excellent choice.
448: Methods can be chosen with -ts_ssp_type {rks2,rks3,rk104}
450: rks2: Second order methods with any number s>1 of stages. c_eff = (s-1)/s
452: rks3: Third order methods with s=n^2 stages, n>1. c_eff = (s-n)/s
454: rk104: A 10-stage fourth order method. c_eff = 0.6
456: Level: beginner
458: References:
459: + * - Ketcheson, Highly efficient strong stability preserving Runge Kutta methods with low storage implementations, SISC, 2008.
460: - * - Gottlieb, Ketcheson, and Shu, High order strong stability preserving time discretizations, J Scientific Computing, 2009.
462: .seealso: [](ch_ts), `TSCreate()`, `TS`, `TSSetType()`
464: M*/
465: PETSC_EXTERN PetscErrorCode TSCreate_SSP(TS ts)
466: {
467: TS_SSP *ssp;
469: PetscFunctionBegin;
470: PetscCall(TSSSPInitializePackage());
472: ts->ops->setup = TSSetUp_SSP;
473: ts->ops->step = TSStep_SSP;
474: ts->ops->reset = TSReset_SSP;
475: ts->ops->destroy = TSDestroy_SSP;
476: ts->ops->setfromoptions = TSSetFromOptions_SSP;
477: ts->ops->view = TSView_SSP;
479: PetscCall(PetscNew(&ssp));
480: ts->data = (void *)ssp;
482: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSSSPGetType_C", TSSSPGetType_SSP));
483: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSSSPSetType_C", TSSSPSetType_SSP));
484: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSSSPGetNumStages_C", TSSSPGetNumStages_SSP));
485: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSSSPSetNumStages_C", TSSSPSetNumStages_SSP));
487: PetscCall(TSSSPSetType(ts, TSSSPRKS2));
488: PetscFunctionReturn(PETSC_SUCCESS);
489: }
491: /*@C
492: TSSSPInitializePackage - This function initializes everything in the `TSSSP` package. It is called
493: from `TSInitializePackage()`.
495: Level: developer
497: .seealso: [](ch_ts), `PetscInitialize()`, `TSSSPFinalizePackage()`, `TSInitializePackage()`
498: @*/
499: PetscErrorCode TSSSPInitializePackage(void)
500: {
501: PetscFunctionBegin;
502: if (TSSSPPackageInitialized) PetscFunctionReturn(PETSC_SUCCESS);
503: TSSSPPackageInitialized = PETSC_TRUE;
504: PetscCall(PetscFunctionListAdd(&TSSSPList, TSSSPRKS2, TSSSPStep_RK_2));
505: PetscCall(PetscFunctionListAdd(&TSSSPList, TSSSPRKS3, TSSSPStep_RK_3));
506: PetscCall(PetscFunctionListAdd(&TSSSPList, TSSSPRK104, TSSSPStep_RK_10_4));
507: PetscCall(PetscRegisterFinalize(TSSSPFinalizePackage));
508: PetscFunctionReturn(PETSC_SUCCESS);
509: }
511: /*@C
512: TSSSPFinalizePackage - This function destroys everything in the `TSSSP` package. It is
513: called from `PetscFinalize()`.
515: Level: developer
517: .seealso: [](ch_ts), `PetscFinalize()`, `TSSSPInitiallizePackage()`, `TSInitializePackage()`
518: @*/
519: PetscErrorCode TSSSPFinalizePackage(void)
520: {
521: PetscFunctionBegin;
522: TSSSPPackageInitialized = PETSC_FALSE;
523: PetscCall(PetscFunctionListDestroy(&TSSSPList));
524: PetscFunctionReturn(PETSC_SUCCESS);
525: }