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: }
214: static PetscErrorCode TSReset_SSP(TS ts)
215: {
216: TS_SSP *ssp = (TS_SSP *)ts->data;
218: PetscFunctionBegin;
219: if (ssp->work) PetscCall(VecDestroyVecs(ssp->nwork, &ssp->work));
220: ssp->nwork = 0;
221: ssp->workout = PETSC_FALSE;
222: PetscFunctionReturn(PETSC_SUCCESS);
223: }
225: static PetscErrorCode TSDestroy_SSP(TS ts)
226: {
227: TS_SSP *ssp = (TS_SSP *)ts->data;
229: PetscFunctionBegin;
230: PetscCall(TSReset_SSP(ts));
231: PetscCall(PetscFree(ssp->type_name));
232: PetscCall(PetscFree(ts->data));
233: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSSSPGetType_C", NULL));
234: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSSSPSetType_C", NULL));
235: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSSSPGetNumStages_C", NULL));
236: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSSSPSetNumStages_C", NULL));
237: PetscFunctionReturn(PETSC_SUCCESS);
238: }
240: /*@
241: TSSSPSetType - set the `TSSSP` time integration scheme to use
243: Logically Collective
245: Input Parameters:
246: + ts - time stepping object
247: - ssptype - type of scheme to use
249: Options Database Keys:
250: + -ts_ssp_type <rks2> - Type of `TSSSP` method (one of) rks2 rks3 rk104
251: - -ts_ssp_nstages<rks2: 5, rks3: 9> - Number of stages
253: Level: beginner
255: .seealso: [](ch_ts), `TSSSP`, `TSSSPGetType()`, `TSSSPSetNumStages()`, `TSSSPRKS2`, `TSSSPRKS3`, `TSSSPRK104`
256: @*/
257: PetscErrorCode TSSSPSetType(TS ts, TSSSPType ssptype)
258: {
259: PetscFunctionBegin;
261: PetscAssertPointer(ssptype, 2);
262: PetscTryMethod(ts, "TSSSPSetType_C", (TS, TSSSPType), (ts, ssptype));
263: PetscFunctionReturn(PETSC_SUCCESS);
264: }
266: /*@
267: TSSSPGetType - get the `TSSSP` time integration scheme
269: Logically Collective
271: Input Parameter:
272: . ts - time stepping object
274: Output Parameter:
275: . type - type of scheme being used
277: Level: beginner
279: .seealso: [](ch_ts), `TSSSP`, `TSSSPSetType()`, `TSSSPSetNumStages()`, `TSSSPRKS2`, `TSSSPRKS3`, `TSSSPRK104`
280: @*/
281: PetscErrorCode TSSSPGetType(TS ts, TSSSPType *type)
282: {
283: PetscFunctionBegin;
285: PetscUseMethod(ts, "TSSSPGetType_C", (TS, TSSSPType *), (ts, type));
286: PetscFunctionReturn(PETSC_SUCCESS);
287: }
289: /*@
290: TSSSPSetNumStages - set the number of stages to use with the `TSSSP` method. Must be called after
291: `TSSSPSetType()`.
293: Logically Collective
295: Input Parameters:
296: + ts - time stepping object
297: - nstages - number of stages
299: Options Database Keys:
300: + -ts_ssp_type <rks2> - Type of `TSSSP` method (one of) rks2 rks3 rk104
301: - -ts_ssp_nstages<rks2: 5, rks3: 9> - Number of stages
303: Level: beginner
305: .seealso: [](ch_ts), `TSSSP`, `TSSSPGetNumStages()`, `TSSSPRKS2`, `TSSSPRKS3`, `TSSSPRK104`
306: @*/
307: PetscErrorCode TSSSPSetNumStages(TS ts, PetscInt nstages)
308: {
309: PetscFunctionBegin;
311: PetscTryMethod(ts, "TSSSPSetNumStages_C", (TS, PetscInt), (ts, nstages));
312: PetscFunctionReturn(PETSC_SUCCESS);
313: }
315: /*@
316: TSSSPGetNumStages - get the number of stages in the `TSSSP` time integration scheme
318: Logically Collective
320: Input Parameter:
321: . ts - time stepping object
323: Output Parameter:
324: . nstages - number of stages
326: Level: beginner
328: .seealso: [](ch_ts), `TSSSP`, `TSSSPGetType()`, `TSSSPSetNumStages()`, `TSSSPRKS2`, `TSSSPRKS3`, `TSSSPRK104`
329: @*/
330: PetscErrorCode TSSSPGetNumStages(TS ts, PetscInt *nstages)
331: {
332: PetscFunctionBegin;
334: PetscUseMethod(ts, "TSSSPGetNumStages_C", (TS, PetscInt *), (ts, nstages));
335: PetscFunctionReturn(PETSC_SUCCESS);
336: }
338: static PetscErrorCode TSSSPSetType_SSP(TS ts, TSSSPType type)
339: {
340: TS_SSP *ssp = (TS_SSP *)ts->data;
341: PetscErrorCode (*r)(TS, PetscReal, PetscReal, Vec);
342: PetscBool flag;
344: PetscFunctionBegin;
345: PetscCall(PetscFunctionListFind(TSSSPList, type, &r));
346: PetscCheck(r, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown TS_SSP type %s given", type);
347: ssp->onestep = r;
348: PetscCall(PetscFree(ssp->type_name));
349: PetscCall(PetscStrallocpy(type, &ssp->type_name));
350: ts->default_adapt_type = TSADAPTNONE;
351: PetscCall(PetscStrcmp(type, TSSSPRKS2, &flag));
352: if (flag) {
353: ssp->nstages = 5;
354: } else {
355: PetscCall(PetscStrcmp(type, TSSSPRKS3, &flag));
356: if (flag) {
357: ssp->nstages = 9;
358: } else {
359: ssp->nstages = 5;
360: }
361: }
362: PetscFunctionReturn(PETSC_SUCCESS);
363: }
364: static PetscErrorCode TSSSPGetType_SSP(TS ts, TSSSPType *type)
365: {
366: TS_SSP *ssp = (TS_SSP *)ts->data;
368: PetscFunctionBegin;
369: *type = ssp->type_name;
370: PetscFunctionReturn(PETSC_SUCCESS);
371: }
372: static PetscErrorCode TSSSPSetNumStages_SSP(TS ts, PetscInt nstages)
373: {
374: TS_SSP *ssp = (TS_SSP *)ts->data;
376: PetscFunctionBegin;
377: ssp->nstages = nstages;
378: PetscFunctionReturn(PETSC_SUCCESS);
379: }
380: static PetscErrorCode TSSSPGetNumStages_SSP(TS ts, PetscInt *nstages)
381: {
382: TS_SSP *ssp = (TS_SSP *)ts->data;
384: PetscFunctionBegin;
385: *nstages = ssp->nstages;
386: PetscFunctionReturn(PETSC_SUCCESS);
387: }
389: static PetscErrorCode TSSetFromOptions_SSP(TS ts, PetscOptionItems PetscOptionsObject)
390: {
391: char tname[256] = TSSSPRKS2;
392: TS_SSP *ssp = (TS_SSP *)ts->data;
393: PetscBool flg;
395: PetscFunctionBegin;
396: PetscOptionsHeadBegin(PetscOptionsObject, "SSP ODE solver options");
397: {
398: PetscCall(PetscOptionsFList("-ts_ssp_type", "Type of SSP method", "TSSSPSetType", TSSSPList, tname, tname, sizeof(tname), &flg));
399: if (flg) PetscCall(TSSSPSetType(ts, tname));
400: PetscCall(PetscOptionsInt("-ts_ssp_nstages", "Number of stages", "TSSSPSetNumStages", ssp->nstages, &ssp->nstages, NULL));
401: }
402: PetscOptionsHeadEnd();
403: PetscFunctionReturn(PETSC_SUCCESS);
404: }
406: static PetscErrorCode TSView_SSP(TS ts, PetscViewer viewer)
407: {
408: TS_SSP *ssp = (TS_SSP *)ts->data;
409: PetscBool ascii;
411: PetscFunctionBegin;
412: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &ascii));
413: if (ascii) PetscCall(PetscViewerASCIIPrintf(viewer, " Scheme: %s\n", ssp->type_name));
414: PetscFunctionReturn(PETSC_SUCCESS);
415: }
417: /*MC
418: TSSSP - Explicit strong stability preserving ODE solver {cite}`ketcheson_2008` {cite}`gottliebketchesonshu2009`
420: Most hyperbolic conservation laws have exact solutions that are total variation diminishing (TVD) or total variation
421: bounded (TVB) although these solutions often contain discontinuities. Spatial discretizations such as Godunov's
422: scheme and high-resolution finite volume methods (TVD limiters, ENO/WENO) are designed to preserve these properties,
423: but they are usually formulated using a forward Euler time discretization or by coupling the space and time
424: discretization as in the classical Lax-Wendroff scheme. When the space and time discretization is coupled, it is very
425: difficult to produce schemes with high temporal accuracy while preserving TVD properties. An alternative is the
426: semidiscrete formulation where we choose a spatial discretization that is TVD with forward Euler and then choose a
427: time discretization that preserves the TVD property. Such integrators are called strong stability preserving (SSP).
429: Let c_eff be the minimum number of function evaluations required to step as far as one step of forward Euler while
430: still being SSP. Some theoretical bounds
432: 1. There are no explicit methods with c_eff > 1.
434: 2. There are no explicit methods beyond order 4 (for nonlinear problems) and c_eff > 0.
436: 3. There are no implicit methods with order greater than 1 and c_eff > 2.
438: This integrator provides Runge-Kutta methods of order 2, 3, and 4 with maximal values of c_eff. More stages allows
439: for larger values of c_eff which improves efficiency. These implementations are low-memory and only use 2 or 3 work
440: vectors regardless of the total number of stages, so e.g. 25-stage 3rd order methods may be an excellent choice.
442: Methods can be chosen with -ts_ssp_type {rks2,rks3,rk104}
444: rks2: Second order methods with any number s>1 of stages. c_eff = (s-1)/s
446: rks3: Third order methods with s=n^2 stages, n>1. c_eff = (s-n)/s
448: rk104: A 10-stage fourth order method. c_eff = 0.6
450: Level: beginner
452: .seealso: [](ch_ts), `TSCreate()`, `TS`, `TSSetType()`
453: M*/
454: PETSC_EXTERN PetscErrorCode TSCreate_SSP(TS ts)
455: {
456: TS_SSP *ssp;
458: PetscFunctionBegin;
459: PetscCall(TSSSPInitializePackage());
461: ts->ops->setup = TSSetUp_SSP;
462: ts->ops->step = TSStep_SSP;
463: ts->ops->reset = TSReset_SSP;
464: ts->ops->destroy = TSDestroy_SSP;
465: ts->ops->setfromoptions = TSSetFromOptions_SSP;
466: ts->ops->view = TSView_SSP;
468: PetscCall(PetscNew(&ssp));
469: ts->data = (void *)ssp;
471: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSSSPGetType_C", TSSSPGetType_SSP));
472: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSSSPSetType_C", TSSSPSetType_SSP));
473: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSSSPGetNumStages_C", TSSSPGetNumStages_SSP));
474: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSSSPSetNumStages_C", TSSSPSetNumStages_SSP));
476: PetscCall(TSSSPSetType(ts, TSSSPRKS2));
477: PetscFunctionReturn(PETSC_SUCCESS);
478: }
480: /*@C
481: TSSSPInitializePackage - This function initializes everything in the `TSSSP` package. It is called
482: from `TSInitializePackage()`.
484: Level: developer
486: .seealso: [](ch_ts), `PetscInitialize()`, `TSSSPFinalizePackage()`, `TSInitializePackage()`
487: @*/
488: PetscErrorCode TSSSPInitializePackage(void)
489: {
490: PetscFunctionBegin;
491: if (TSSSPPackageInitialized) PetscFunctionReturn(PETSC_SUCCESS);
492: TSSSPPackageInitialized = PETSC_TRUE;
493: PetscCall(PetscFunctionListAdd(&TSSSPList, TSSSPRKS2, TSSSPStep_RK_2));
494: PetscCall(PetscFunctionListAdd(&TSSSPList, TSSSPRKS3, TSSSPStep_RK_3));
495: PetscCall(PetscFunctionListAdd(&TSSSPList, TSSSPRK104, TSSSPStep_RK_10_4));
496: PetscCall(PetscRegisterFinalize(TSSSPFinalizePackage));
497: PetscFunctionReturn(PETSC_SUCCESS);
498: }
500: /*@C
501: TSSSPFinalizePackage - This function destroys everything in the `TSSSP` package. It is
502: called from `PetscFinalize()`.
504: Level: developer
506: .seealso: [](ch_ts), `PetscFinalize()`, `TSSSPInitiallizePackage()`, `TSInitializePackage()`
507: @*/
508: PetscErrorCode TSSSPFinalizePackage(void)
509: {
510: PetscFunctionBegin;
511: TSSSPPackageInitialized = PETSC_FALSE;
512: PetscCall(PetscFunctionListDestroy(&TSSSPList));
513: PetscFunctionReturn(PETSC_SUCCESS);
514: }