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. Pseudocode 2 of {cite}`ketcheson_2008`
49: Level: beginner
51: .seealso: [](ch_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: PetscFunctionBegin;
60: s = ssp->nstages;
61: PetscCall(TSSSPGetWorkVectors(ts, 2, &work));
62: F = work[1];
63: PetscCall(VecCopy(sol, work[0]));
64: for (i = 0; i < s - 1; i++) {
65: PetscReal stage_time = t0 + dt * (i / (s - 1.));
66: PetscCall(TSPreStage(ts, stage_time));
67: PetscCall(TSComputeRHSFunction(ts, stage_time, work[0], F));
68: PetscCall(VecAXPY(work[0], dt / (s - 1.), F));
69: }
70: PetscCall(TSComputeRHSFunction(ts, t0 + dt, work[0], F));
71: PetscCall(VecAXPBYPCZ(sol, (s - 1.) / s, dt / s, 1. / s, work[0], F));
72: PetscCall(TSSSPRestoreWorkVectors(ts, 2, &work));
73: PetscFunctionReturn(PETSC_SUCCESS);
74: }
76: /*MC
77: TSSSPRKS3 - Optimal third order SSP Runge-Kutta, low-storage, $c_eff=(PetscSqrtReal(s)-1)/PetscSqrtReal(s)$, where `PetscSqrtReal`(s) is an integer
79: Pseudocode 2 of {cite}`ketcheson_2008`
81: Level: beginner
83: .seealso: [](ch_ts), `TSSSP`, `TSSSPSetType()`, `TSSSPSetNumStages()`
84: M*/
85: static PetscErrorCode TSSSPStep_RK_3(TS ts, PetscReal t0, PetscReal dt, Vec sol)
86: {
87: TS_SSP *ssp = (TS_SSP *)ts->data;
88: Vec *work, F;
89: PetscInt i, s, n, r;
90: PetscReal c, stage_time;
92: PetscFunctionBegin;
93: s = ssp->nstages;
94: n = (PetscInt)(PetscSqrtReal((PetscReal)s) + 0.001);
95: r = s - n;
96: 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);
97: PetscCall(TSSSPGetWorkVectors(ts, 3, &work));
98: F = work[2];
99: PetscCall(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: PetscCall(TSPreStage(ts, stage_time));
104: PetscCall(TSComputeRHSFunction(ts, stage_time, work[0], F));
105: PetscCall(VecAXPY(work[0], dt / r, F));
106: }
107: PetscCall(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: PetscCall(TSPreStage(ts, stage_time));
112: PetscCall(TSComputeRHSFunction(ts, stage_time, work[0], F));
113: PetscCall(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: PetscCall(TSPreStage(ts, stage_time));
119: PetscCall(TSComputeRHSFunction(ts, stage_time, work[0], F));
120: PetscCall(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: PetscCall(TSPreStage(ts, stage_time));
127: PetscCall(TSComputeRHSFunction(ts, stage_time, work[0], F));
128: PetscCall(VecAXPY(work[0], dt / r, F));
129: }
130: PetscCall(VecCopy(work[0], sol));
131: PetscCall(TSSSPRestoreWorkVectors(ts, 3, &work));
132: PetscFunctionReturn(PETSC_SUCCESS);
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 {cite}`ketcheson_2008`
140: Level: beginner
142: .seealso: [](ch_ts), `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: PetscFunctionBegin;
152: PetscCall(TSSSPGetWorkVectors(ts, 3, &work));
153: F = work[2];
154: PetscCall(VecCopy(sol, work[0]));
155: for (i = 0; i < 5; i++) {
156: stage_time = t0 + c[i] * dt;
157: PetscCall(TSPreStage(ts, stage_time));
158: PetscCall(TSComputeRHSFunction(ts, stage_time, work[0], F));
159: PetscCall(VecAXPY(work[0], dt / 6, F));
160: }
161: PetscCall(VecAXPBYPCZ(work[1], 1. / 25, 9. / 25, 0, sol, work[0]));
162: PetscCall(VecAXPBY(work[0], 15, -5, work[1]));
163: for (; i < 9; i++) {
164: stage_time = t0 + c[i] * dt;
165: PetscCall(TSPreStage(ts, stage_time));
166: PetscCall(TSComputeRHSFunction(ts, stage_time, work[0], F));
167: PetscCall(VecAXPY(work[0], dt / 6, F));
168: }
169: stage_time = t0 + dt;
170: PetscCall(TSPreStage(ts, stage_time));
171: PetscCall(TSComputeRHSFunction(ts, stage_time, work[0], F));
172: PetscCall(VecAXPBYPCZ(work[1], 3. / 5, dt / 10, 1, work[0], F));
173: PetscCall(VecCopy(work[1], sol));
174: PetscCall(TSSSPRestoreWorkVectors(ts, 3, &work));
175: PetscFunctionReturn(PETSC_SUCCESS);
176: }
178: static PetscErrorCode TSSetUp_SSP(TS ts)
179: {
180: PetscFunctionBegin;
181: PetscCall(TSCheckImplicitTerm(ts));
182: PetscCall(TSGetAdapt(ts, &ts->adapt));
183: PetscCall(TSAdaptCandidatesClear(ts->adapt));
184: PetscFunctionReturn(PETSC_SUCCESS);
185: }
187: static PetscErrorCode TSStep_SSP(TS ts)
188: {
189: TS_SSP *ssp = (TS_SSP *)ts->data;
190: Vec sol = ts->vec_sol;
191: PetscBool stageok, accept = PETSC_TRUE;
192: PetscReal next_time_step = ts->time_step;
194: PetscFunctionBegin;
195: PetscCall((*ssp->onestep)(ts, ts->ptime, ts->time_step, sol));
196: PetscCall(TSPostStage(ts, ts->ptime, 0, &sol));
197: PetscCall(TSAdaptCheckStage(ts->adapt, ts, ts->ptime + ts->time_step, sol, &stageok));
198: if (!stageok) {
199: ts->reason = TS_DIVERGED_STEP_REJECTED;
200: PetscFunctionReturn(PETSC_SUCCESS);
201: }
203: PetscCall(TSAdaptChoose(ts->adapt, ts, ts->time_step, NULL, &next_time_step, &accept));
204: if (!accept) {
205: ts->reason = TS_DIVERGED_STEP_REJECTED;
206: PetscFunctionReturn(PETSC_SUCCESS);
207: }
209: ts->ptime += ts->time_step;
210: ts->time_step = next_time_step;
211: PetscFunctionReturn(PETSC_SUCCESS);
212: }
213: /*------------------------------------------------------------*/
215: static PetscErrorCode TSReset_SSP(TS ts)
216: {
217: TS_SSP *ssp = (TS_SSP *)ts->data;
219: PetscFunctionBegin;
220: if (ssp->work) PetscCall(VecDestroyVecs(ssp->nwork, &ssp->work));
221: ssp->nwork = 0;
222: ssp->workout = PETSC_FALSE;
223: PetscFunctionReturn(PETSC_SUCCESS);
224: }
226: static PetscErrorCode TSDestroy_SSP(TS ts)
227: {
228: TS_SSP *ssp = (TS_SSP *)ts->data;
230: PetscFunctionBegin;
231: PetscCall(TSReset_SSP(ts));
232: PetscCall(PetscFree(ssp->type_name));
233: PetscCall(PetscFree(ts->data));
234: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSSSPGetType_C", NULL));
235: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSSSPSetType_C", NULL));
236: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSSSPGetNumStages_C", NULL));
237: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSSSPSetNumStages_C", NULL));
238: PetscFunctionReturn(PETSC_SUCCESS);
239: }
240: /*------------------------------------------------------------*/
242: /*@
243: TSSSPSetType - set the `TSSSP` time integration scheme to use
245: Logically Collective
247: Input Parameters:
248: + ts - time stepping object
249: - ssptype - type of scheme to use
251: Options Database Keys:
252: + -ts_ssp_type <rks2> - Type of `TSSSP` method (one of) rks2 rks3 rk104
253: - -ts_ssp_nstages<rks2: 5, rks3: 9> - Number of stages
255: Level: beginner
257: .seealso: [](ch_ts), `TSSSP`, `TSSSPGetType()`, `TSSSPSetNumStages()`, `TSSSPRKS2`, `TSSSPRKS3`, `TSSSPRK104`
258: @*/
259: PetscErrorCode TSSSPSetType(TS ts, TSSSPType ssptype)
260: {
261: PetscFunctionBegin;
263: PetscAssertPointer(ssptype, 2);
264: PetscTryMethod(ts, "TSSSPSetType_C", (TS, TSSSPType), (ts, ssptype));
265: PetscFunctionReturn(PETSC_SUCCESS);
266: }
268: /*@
269: TSSSPGetType - get the `TSSSP` time integration scheme
271: Logically Collective
273: Input Parameter:
274: . ts - time stepping object
276: Output Parameter:
277: . type - type of scheme being used
279: Level: beginner
281: .seealso: [](ch_ts), `TSSSP`, `TSSSPSettype()`, `TSSSPSetNumStages()`, `TSSSPRKS2`, `TSSSPRKS3`, `TSSSPRK104`
282: @*/
283: PetscErrorCode TSSSPGetType(TS ts, TSSSPType *type)
284: {
285: PetscFunctionBegin;
287: PetscUseMethod(ts, "TSSSPGetType_C", (TS, TSSSPType *), (ts, type));
288: PetscFunctionReturn(PETSC_SUCCESS);
289: }
291: /*@
292: TSSSPSetNumStages - set the number of stages to use with the `TSSSP` method. Must be called after
293: `TSSSPSetType()`.
295: Logically Collective
297: Input Parameters:
298: + ts - time stepping object
299: - nstages - number of stages
301: Options Database Keys:
302: + -ts_ssp_type <rks2> - Type of `TSSSP` method (one of) rks2 rks3 rk104
303: - -ts_ssp_nstages<rks2: 5, rks3: 9> - Number of stages
305: Level: beginner
307: .seealso: [](ch_ts), `TSSSP`, `TSSSPGetNumStages()`, `TSSSPRKS2`, `TSSSPRKS3`, `TSSSPRK104`
308: @*/
309: PetscErrorCode TSSSPSetNumStages(TS ts, PetscInt nstages)
310: {
311: PetscFunctionBegin;
313: PetscTryMethod(ts, "TSSSPSetNumStages_C", (TS, PetscInt), (ts, nstages));
314: PetscFunctionReturn(PETSC_SUCCESS);
315: }
317: /*@
318: TSSSPGetNumStages - get the number of stages in the `TSSSP` time integration scheme
320: Logically Collective
322: Input Parameter:
323: . ts - time stepping object
325: Output Parameter:
326: . nstages - number of stages
328: Level: beginner
330: .seealso: [](ch_ts), `TSSSP`, `TSSSPGetType()`, `TSSSPSetNumStages()`, `TSSSPRKS2`, `TSSSPRKS3`, `TSSSPRK104`
331: @*/
332: PetscErrorCode TSSSPGetNumStages(TS ts, PetscInt *nstages)
333: {
334: PetscFunctionBegin;
336: PetscUseMethod(ts, "TSSSPGetNumStages_C", (TS, PetscInt *), (ts, nstages));
337: PetscFunctionReturn(PETSC_SUCCESS);
338: }
340: static PetscErrorCode TSSSPSetType_SSP(TS ts, TSSSPType type)
341: {
342: TS_SSP *ssp = (TS_SSP *)ts->data;
343: PetscErrorCode (*r)(TS, PetscReal, PetscReal, Vec);
344: PetscBool flag;
346: PetscFunctionBegin;
347: PetscCall(PetscFunctionListFind(TSSSPList, type, &r));
348: PetscCheck(r, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown TS_SSP type %s given", type);
349: ssp->onestep = r;
350: PetscCall(PetscFree(ssp->type_name));
351: PetscCall(PetscStrallocpy(type, &ssp->type_name));
352: ts->default_adapt_type = TSADAPTNONE;
353: PetscCall(PetscStrcmp(type, TSSSPRKS2, &flag));
354: if (flag) {
355: ssp->nstages = 5;
356: } else {
357: PetscCall(PetscStrcmp(type, TSSSPRKS3, &flag));
358: if (flag) {
359: ssp->nstages = 9;
360: } else {
361: ssp->nstages = 5;
362: }
363: }
364: PetscFunctionReturn(PETSC_SUCCESS);
365: }
366: static PetscErrorCode TSSSPGetType_SSP(TS ts, TSSSPType *type)
367: {
368: TS_SSP *ssp = (TS_SSP *)ts->data;
370: PetscFunctionBegin;
371: *type = ssp->type_name;
372: PetscFunctionReturn(PETSC_SUCCESS);
373: }
374: static PetscErrorCode TSSSPSetNumStages_SSP(TS ts, PetscInt nstages)
375: {
376: TS_SSP *ssp = (TS_SSP *)ts->data;
378: PetscFunctionBegin;
379: ssp->nstages = nstages;
380: PetscFunctionReturn(PETSC_SUCCESS);
381: }
382: static PetscErrorCode TSSSPGetNumStages_SSP(TS ts, PetscInt *nstages)
383: {
384: TS_SSP *ssp = (TS_SSP *)ts->data;
386: PetscFunctionBegin;
387: *nstages = ssp->nstages;
388: PetscFunctionReturn(PETSC_SUCCESS);
389: }
391: static PetscErrorCode TSSetFromOptions_SSP(TS ts, PetscOptionItems *PetscOptionsObject)
392: {
393: char tname[256] = TSSSPRKS2;
394: TS_SSP *ssp = (TS_SSP *)ts->data;
395: PetscBool flg;
397: PetscFunctionBegin;
398: PetscOptionsHeadBegin(PetscOptionsObject, "SSP ODE solver options");
399: {
400: PetscCall(PetscOptionsFList("-ts_ssp_type", "Type of SSP method", "TSSSPSetType", TSSSPList, tname, tname, sizeof(tname), &flg));
401: if (flg) PetscCall(TSSSPSetType(ts, tname));
402: PetscCall(PetscOptionsInt("-ts_ssp_nstages", "Number of stages", "TSSSPSetNumStages", ssp->nstages, &ssp->nstages, NULL));
403: }
404: PetscOptionsHeadEnd();
405: PetscFunctionReturn(PETSC_SUCCESS);
406: }
408: static PetscErrorCode TSView_SSP(TS ts, PetscViewer viewer)
409: {
410: TS_SSP *ssp = (TS_SSP *)ts->data;
411: PetscBool ascii;
413: PetscFunctionBegin;
414: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &ascii));
415: if (ascii) PetscCall(PetscViewerASCIIPrintf(viewer, " Scheme: %s\n", ssp->type_name));
416: PetscFunctionReturn(PETSC_SUCCESS);
417: }
419: /* ------------------------------------------------------------ */
421: /*MC
422: TSSSP - Explicit strong stability preserving ODE solver {cite}`ketcheson_2008` {cite}`gottliebketchesonshu2009`
424: Most hyperbolic conservation laws have exact solutions that are total variation diminishing (TVD) or total variation
425: bounded (TVB) although these solutions often contain discontinuities. Spatial discretizations such as Godunov's
426: scheme and high-resolution finite volume methods (TVD limiters, ENO/WENO) are designed to preserve these properties,
427: but they are usually formulated using a forward Euler time discretization or by coupling the space and time
428: discretization as in the classical Lax-Wendroff scheme. When the space and time discretization is coupled, it is very
429: difficult to produce schemes with high temporal accuracy while preserving TVD properties. An alternative is the
430: semidiscrete formulation where we choose a spatial discretization that is TVD with forward Euler and then choose a
431: time discretization that preserves the TVD property. Such integrators are called strong stability preserving (SSP).
433: Let c_eff be the minimum number of function evaluations required to step as far as one step of forward Euler while
434: still being SSP. Some theoretical bounds
436: 1. There are no explicit methods with c_eff > 1.
438: 2. There are no explicit methods beyond order 4 (for nonlinear problems) and c_eff > 0.
440: 3. There are no implicit methods with order greater than 1 and c_eff > 2.
442: This integrator provides Runge-Kutta methods of order 2, 3, and 4 with maximal values of c_eff. More stages allows
443: for larger values of c_eff which improves efficiency. These implementations are low-memory and only use 2 or 3 work
444: vectors regardless of the total number of stages, so e.g. 25-stage 3rd order methods may be an excellent choice.
446: Methods can be chosen with -ts_ssp_type {rks2,rks3,rk104}
448: rks2: Second order methods with any number s>1 of stages. c_eff = (s-1)/s
450: rks3: Third order methods with s=n^2 stages, n>1. c_eff = (s-n)/s
452: rk104: A 10-stage fourth order method. c_eff = 0.6
454: Level: beginner
456: .seealso: [](ch_ts), `TSCreate()`, `TS`, `TSSetType()`
457: M*/
458: PETSC_EXTERN PetscErrorCode TSCreate_SSP(TS ts)
459: {
460: TS_SSP *ssp;
462: PetscFunctionBegin;
463: PetscCall(TSSSPInitializePackage());
465: ts->ops->setup = TSSetUp_SSP;
466: ts->ops->step = TSStep_SSP;
467: ts->ops->reset = TSReset_SSP;
468: ts->ops->destroy = TSDestroy_SSP;
469: ts->ops->setfromoptions = TSSetFromOptions_SSP;
470: ts->ops->view = TSView_SSP;
472: PetscCall(PetscNew(&ssp));
473: ts->data = (void *)ssp;
475: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSSSPGetType_C", TSSSPGetType_SSP));
476: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSSSPSetType_C", TSSSPSetType_SSP));
477: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSSSPGetNumStages_C", TSSSPGetNumStages_SSP));
478: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSSSPSetNumStages_C", TSSSPSetNumStages_SSP));
480: PetscCall(TSSSPSetType(ts, TSSSPRKS2));
481: PetscFunctionReturn(PETSC_SUCCESS);
482: }
484: /*@C
485: TSSSPInitializePackage - This function initializes everything in the `TSSSP` package. It is called
486: from `TSInitializePackage()`.
488: Level: developer
490: .seealso: [](ch_ts), `PetscInitialize()`, `TSSSPFinalizePackage()`, `TSInitializePackage()`
491: @*/
492: PetscErrorCode TSSSPInitializePackage(void)
493: {
494: PetscFunctionBegin;
495: if (TSSSPPackageInitialized) PetscFunctionReturn(PETSC_SUCCESS);
496: TSSSPPackageInitialized = PETSC_TRUE;
497: PetscCall(PetscFunctionListAdd(&TSSSPList, TSSSPRKS2, TSSSPStep_RK_2));
498: PetscCall(PetscFunctionListAdd(&TSSSPList, TSSSPRKS3, TSSSPStep_RK_3));
499: PetscCall(PetscFunctionListAdd(&TSSSPList, TSSSPRK104, TSSSPStep_RK_10_4));
500: PetscCall(PetscRegisterFinalize(TSSSPFinalizePackage));
501: PetscFunctionReturn(PETSC_SUCCESS);
502: }
504: /*@C
505: TSSSPFinalizePackage - This function destroys everything in the `TSSSP` package. It is
506: called from `PetscFinalize()`.
508: Level: developer
510: .seealso: [](ch_ts), `PetscFinalize()`, `TSSSPInitiallizePackage()`, `TSInitializePackage()`
511: @*/
512: PetscErrorCode TSSSPFinalizePackage(void)
513: {
514: PetscFunctionBegin;
515: TSSSPPackageInitialized = PETSC_FALSE;
516: PetscCall(PetscFunctionListDestroy(&TSSSPList));
517: PetscFunctionReturn(PETSC_SUCCESS);
518: }