Actual source code: ssp.c
petsc-3.7.7 2017-09-25
1: /*
2: Code for Timestepping with explicit SSP.
3: */
4: #include <petsc/private/tsimpl.h> /*I "petscts.h" I*/
6: PetscFunctionList TSSSPList = 0;
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;
21: static PetscErrorCode TSSSPGetWorkVectors(TS ts,PetscInt n,Vec **work)
22: {
23: TS_SSP *ssp = (TS_SSP*)ts->data;
27: if (ssp->workout) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Work vectors already gotten");
28: if (ssp->nwork < n) {
29: if (ssp->nwork > 0) {
30: VecDestroyVecs(ssp->nwork,&ssp->work);
31: }
32: VecDuplicateVecs(ts->vec_sol,n,&ssp->work);
33: ssp->nwork = n;
34: }
35: *work = ssp->work;
36: ssp->workout = PETSC_TRUE;
37: return(0);
38: }
42: static PetscErrorCode TSSSPRestoreWorkVectors(TS ts,PetscInt n,Vec **work)
43: {
44: TS_SSP *ssp = (TS_SSP*)ts->data;
47: if (!ssp->workout) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Work vectors have not been gotten");
48: if (*work != ssp->work) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Wrong work vectors checked out");
49: ssp->workout = PETSC_FALSE;
50: *work = NULL;
51: return(0);
52: }
56: /*MC
57: TSSSPRKS2 - Optimal second order SSP Runge-Kutta method, low-storage, c_eff=(s-1)/s
59: Pseudocode 2 of Ketcheson 2008
61: Level: beginner
63: .seealso: TSSSP, TSSSPSetType(), TSSSPSetNumStages()
64: M*/
65: static PetscErrorCode TSSSPStep_RK_2(TS ts,PetscReal t0,PetscReal dt,Vec sol)
66: {
67: TS_SSP *ssp = (TS_SSP*)ts->data;
68: Vec *work,F;
69: PetscInt i,s;
73: s = ssp->nstages;
74: TSSSPGetWorkVectors(ts,2,&work);
75: F = work[1];
76: VecCopy(sol,work[0]);
77: for (i=0; i<s-1; i++) {
78: PetscReal stage_time = t0+dt*(i/(s-1.));
79: TSPreStage(ts,stage_time);
80: TSComputeRHSFunction(ts,stage_time,work[0],F);
81: VecAXPY(work[0],dt/(s-1.),F);
82: }
83: TSComputeRHSFunction(ts,t0+dt,work[0],F);
84: VecAXPBYPCZ(sol,(s-1.)/s,dt/s,1./s,work[0],F);
85: TSSSPRestoreWorkVectors(ts,2,&work);
86: return(0);
87: }
91: /*MC
92: TSSSPRKS3 - Optimal third order SSP Runge-Kutta, low-storage, c_eff=(PetscSqrtReal(s)-1)/PetscSqrtReal(s), where PetscSqrtReal(s) is an integer
94: Pseudocode 2 of Ketcheson 2008
96: Level: beginner
98: .seealso: TSSSP, TSSSPSetType(), TSSSPSetNumStages()
99: M*/
100: static PetscErrorCode TSSSPStep_RK_3(TS ts,PetscReal t0,PetscReal dt,Vec sol)
101: {
102: TS_SSP *ssp = (TS_SSP*)ts->data;
103: Vec *work,F;
104: PetscInt i,s,n,r;
105: PetscReal c,stage_time;
109: s = ssp->nstages;
110: n = (PetscInt)(PetscSqrtReal((PetscReal)s)+0.001);
111: r = s-n;
112: 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);
113: TSSSPGetWorkVectors(ts,3,&work);
114: F = work[2];
115: VecCopy(sol,work[0]);
116: for (i=0; i<(n-1)*(n-2)/2; i++) {
117: c = (i<n*(n+1)/2) ? 1.*i/(s-n) : (1.*i-n)/(s-n);
118: stage_time = t0+c*dt;
119: TSPreStage(ts,stage_time);
120: TSComputeRHSFunction(ts,stage_time,work[0],F);
121: VecAXPY(work[0],dt/r,F);
122: }
123: VecCopy(work[0],work[1]);
124: for (; i<n*(n+1)/2-1; i++) {
125: c = (i<n*(n+1)/2) ? 1.*i/(s-n) : (1.*i-n)/(s-n);
126: stage_time = t0+c*dt;
127: TSPreStage(ts,stage_time);
128: TSComputeRHSFunction(ts,stage_time,work[0],F);
129: VecAXPY(work[0],dt/r,F);
130: }
131: {
132: c = (i<n*(n+1)/2) ? 1.*i/(s-n) : (1.*i-n)/(s-n);
133: stage_time = t0+c*dt;
134: TSPreStage(ts,stage_time);
135: TSComputeRHSFunction(ts,stage_time,work[0],F);
136: VecAXPBYPCZ(work[0],1.*n/(2*n-1.),(n-1.)*dt/(r*(2*n-1)),(n-1.)/(2*n-1.),work[1],F);
137: i++;
138: }
139: for (; i<s; i++) {
140: c = (i<n*(n+1)/2) ? 1.*i/(s-n) : (1.*i-n)/(s-n);
141: stage_time = t0+c*dt;
142: TSPreStage(ts,stage_time);
143: TSComputeRHSFunction(ts,stage_time,work[0],F);
144: VecAXPY(work[0],dt/r,F);
145: }
146: VecCopy(work[0],sol);
147: TSSSPRestoreWorkVectors(ts,3,&work);
148: return(0);
149: }
153: /*MC
154: TSSSPRKS104 - Optimal fourth order SSP Runge-Kutta, low-storage (2N), c_eff=0.6
156: SSPRK(10,4), Pseudocode 3 of Ketcheson 2008
158: Level: beginner
160: .seealso: TSSSP, TSSSPSetType()
161: M*/
162: static PetscErrorCode TSSSPStep_RK_10_4(TS ts,PetscReal t0,PetscReal dt,Vec sol)
163: {
164: const PetscReal c[10] = {0, 1./6, 2./6, 3./6, 4./6, 2./6, 3./6, 4./6, 5./6, 1};
165: Vec *work,F;
166: PetscInt i;
167: PetscReal stage_time;
168: PetscErrorCode ierr;
171: TSSSPGetWorkVectors(ts,3,&work);
172: F = work[2];
173: VecCopy(sol,work[0]);
174: for (i=0; i<5; i++) {
175: stage_time = t0+c[i]*dt;
176: TSPreStage(ts,stage_time);
177: TSComputeRHSFunction(ts,stage_time,work[0],F);
178: VecAXPY(work[0],dt/6,F);
179: }
180: VecAXPBYPCZ(work[1],1./25,9./25,0,sol,work[0]);
181: VecAXPBY(work[0],15,-5,work[1]);
182: for (; i<9; i++) {
183: stage_time = t0+c[i]*dt;
184: TSPreStage(ts,stage_time);
185: TSComputeRHSFunction(ts,stage_time,work[0],F);
186: VecAXPY(work[0],dt/6,F);
187: }
188: stage_time = t0+dt;
189: TSPreStage(ts,stage_time);
190: TSComputeRHSFunction(ts,stage_time,work[0],F);
191: VecAXPBYPCZ(work[1],3./5,dt/10,1,work[0],F);
192: VecCopy(work[1],sol);
193: TSSSPRestoreWorkVectors(ts,3,&work);
194: return(0);
195: }
200: static PetscErrorCode TSSetUp_SSP(TS ts)
201: {
205: TSCheckImplicitTerm(ts);
206: TSGetAdapt(ts,&ts->adapt);
207: TSAdaptCandidatesClear(ts->adapt);
208: return(0);
209: }
213: static PetscErrorCode TSStep_SSP(TS ts)
214: {
215: TS_SSP *ssp = (TS_SSP*)ts->data;
216: Vec sol = ts->vec_sol;
217: PetscBool stageok,accept = PETSC_TRUE;
218: PetscReal next_time_step = ts->time_step;
222: (*ssp->onestep)(ts,ts->ptime,ts->time_step,sol);
223: TSPostStage(ts,ts->ptime,0,&sol);
224: TSAdaptCheckStage(ts->adapt,ts,ts->ptime+ts->time_step,sol,&stageok);
225: if(!stageok) {ts->reason = TS_DIVERGED_STEP_REJECTED; return(0);}
227: TSAdaptChoose(ts->adapt,ts,ts->time_step,NULL,&next_time_step,&accept);
228: if (!accept) {ts->reason = TS_DIVERGED_STEP_REJECTED; return(0);}
230: ts->ptime += ts->time_step;
231: ts->time_step = next_time_step;
232: return(0);
233: }
234: /*------------------------------------------------------------*/
238: static PetscErrorCode TSReset_SSP(TS ts)
239: {
240: TS_SSP *ssp = (TS_SSP*)ts->data;
244: if (ssp->work) {VecDestroyVecs(ssp->nwork,&ssp->work);}
245: ssp->nwork = 0;
246: ssp->workout = PETSC_FALSE;
247: return(0);
248: }
252: static PetscErrorCode TSDestroy_SSP(TS ts)
253: {
254: TS_SSP *ssp = (TS_SSP*)ts->data;
258: TSReset_SSP(ts);
259: PetscFree(ssp->type_name);
260: PetscFree(ts->data);
261: PetscObjectComposeFunction((PetscObject)ts,"TSSSPGetType_C",NULL);
262: PetscObjectComposeFunction((PetscObject)ts,"TSSSPSetType_C",NULL);
263: PetscObjectComposeFunction((PetscObject)ts,"TSSSPGetNumStages_C",NULL);
264: PetscObjectComposeFunction((PetscObject)ts,"TSSSPSetNumStages_C",NULL);
265: return(0);
266: }
267: /*------------------------------------------------------------*/
271: /*@C
272: TSSSPSetType - set the SSP time integration scheme to use
274: Logically Collective
276: Input Arguments:
277: ts - time stepping object
278: type - type of scheme to use
280: Options Database Keys:
281: -ts_ssp_type <rks2>: Type of SSP method (one of) rks2 rks3 rk104
282: -ts_ssp_nstages <5>: Number of stages
284: Level: beginner
286: .seealso: TSSSP, TSSSPGetType(), TSSSPSetNumStages(), TSSSPRKS2, TSSSPRKS3, TSSSPRK104
287: @*/
288: PetscErrorCode TSSSPSetType(TS ts,TSSSPType type)
289: {
294: PetscTryMethod(ts,"TSSSPSetType_C",(TS,TSSSPType),(ts,type));
295: return(0);
296: }
300: /*@C
301: TSSSPGetType - get the SSP time integration scheme
303: Logically Collective
305: Input Argument:
306: ts - time stepping object
308: Output Argument:
309: type - type of scheme being used
311: Level: beginner
313: .seealso: TSSSP, TSSSPSettype(), TSSSPSetNumStages(), TSSSPRKS2, TSSSPRKS3, TSSSPRK104
314: @*/
315: PetscErrorCode TSSSPGetType(TS ts,TSSSPType *type)
316: {
321: PetscUseMethod(ts,"TSSSPGetType_C",(TS,TSSSPType*),(ts,type));
322: return(0);
323: }
327: /*@
328: TSSSPSetNumStages - set the number of stages to use with the SSP method
330: Logically Collective
332: Input Arguments:
333: ts - time stepping object
334: nstages - number of stages
336: Options Database Keys:
337: -ts_ssp_type <rks2>: NumStages of SSP method (one of) rks2 rks3 rk104
338: -ts_ssp_nstages <5>: Number of stages
340: Level: beginner
342: .seealso: TSSSP, TSSSPGetNumStages(), TSSSPSetNumStages(), TSSSPRKS2, TSSSPRKS3, TSSSPRK104
343: @*/
344: PetscErrorCode TSSSPSetNumStages(TS ts,PetscInt nstages)
345: {
350: PetscTryMethod(ts,"TSSSPSetNumStages_C",(TS,PetscInt),(ts,nstages));
351: return(0);
352: }
356: /*@
357: TSSSPGetNumStages - get the number of stages in the SSP time integration scheme
359: Logically Collective
361: Input Argument:
362: ts - time stepping object
364: Output Argument:
365: nstages - number of stages
367: Level: beginner
369: .seealso: TSSSP, TSSSPGetType(), TSSSPSetNumStages(), TSSSPRKS2, TSSSPRKS3, TSSSPRK104
370: @*/
371: PetscErrorCode TSSSPGetNumStages(TS ts,PetscInt *nstages)
372: {
377: PetscUseMethod(ts,"TSSSPGetNumStages_C",(TS,PetscInt*),(ts,nstages));
378: return(0);
379: }
383: static PetscErrorCode TSSSPSetType_SSP(TS ts,TSSSPType type)
384: {
385: PetscErrorCode ierr,(*r)(TS,PetscReal,PetscReal,Vec);
386: TS_SSP *ssp = (TS_SSP*)ts->data;
389: PetscFunctionListFind(TSSSPList,type,&r);
390: if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown TS_SSP type %s given",type);
391: ssp->onestep = r;
392: PetscFree(ssp->type_name);
393: PetscStrallocpy(type,&ssp->type_name);
394: return(0);
395: }
398: static PetscErrorCode TSSSPGetType_SSP(TS ts,TSSSPType *type)
399: {
400: TS_SSP *ssp = (TS_SSP*)ts->data;
403: *type = ssp->type_name;
404: return(0);
405: }
408: static PetscErrorCode TSSSPSetNumStages_SSP(TS ts,PetscInt nstages)
409: {
410: TS_SSP *ssp = (TS_SSP*)ts->data;
413: ssp->nstages = nstages;
414: return(0);
415: }
418: static PetscErrorCode TSSSPGetNumStages_SSP(TS ts,PetscInt *nstages)
419: {
420: TS_SSP *ssp = (TS_SSP*)ts->data;
423: *nstages = ssp->nstages;
424: return(0);
425: }
429: static PetscErrorCode TSSetFromOptions_SSP(PetscOptionItems *PetscOptionsObject,TS ts)
430: {
431: char tname[256] = TSSSPRKS2;
432: TS_SSP *ssp = (TS_SSP*)ts->data;
434: PetscBool flg;
437: PetscOptionsHead(PetscOptionsObject,"SSP ODE solver options");
438: {
439: PetscOptionsFList("-ts_ssp_type","Type of SSP method","TSSSPSetType",TSSSPList,tname,tname,sizeof(tname),&flg);
440: if (flg) {
441: TSSSPSetType(ts,tname);
442: }
443: PetscOptionsInt("-ts_ssp_nstages","Number of stages","TSSSPSetNumStages",ssp->nstages,&ssp->nstages,NULL);
444: }
445: PetscOptionsTail();
446: return(0);
447: }
451: static PetscErrorCode TSView_SSP(TS ts,PetscViewer viewer)
452: {
453: TS_SSP *ssp = (TS_SSP*)ts->data;
454: PetscBool ascii;
458: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&ascii);
459: if (ascii) {PetscViewerASCIIPrintf(viewer," Scheme: %s\n",ssp->type_name);}
460: return(0);
461: }
463: /* ------------------------------------------------------------ */
465: /*MC
466: TSSSP - Explicit strong stability preserving ODE solver
468: Most hyperbolic conservation laws have exact solutions that are total variation diminishing (TVD) or total variation
469: bounded (TVB) although these solutions often contain discontinuities. Spatial discretizations such as Godunov's
470: scheme and high-resolution finite volume methods (TVD limiters, ENO/WENO) are designed to preserve these properties,
471: but they are usually formulated using a forward Euler time discretization or by coupling the space and time
472: discretization as in the classical Lax-Wendroff scheme. When the space and time discretization is coupled, it is very
473: difficult to produce schemes with high temporal accuracy while preserving TVD properties. An alternative is the
474: semidiscrete formulation where we choose a spatial discretization that is TVD with forward Euler and then choose a
475: time discretization that preserves the TVD property. Such integrators are called strong stability preserving (SSP).
477: Let c_eff be the minimum number of function evaluations required to step as far as one step of forward Euler while
478: still being SSP. Some theoretical bounds
480: 1. There are no explicit methods with c_eff > 1.
482: 2. There are no explicit methods beyond order 4 (for nonlinear problems) and c_eff > 0.
484: 3. There are no implicit methods with order greater than 1 and c_eff > 2.
486: This integrator provides Runge-Kutta methods of order 2, 3, and 4 with maximal values of c_eff. More stages allows
487: for larger values of c_eff which improves efficiency. These implementations are low-memory and only use 2 or 3 work
488: vectors regardless of the total number of stages, so e.g. 25-stage 3rd order methods may be an excellent choice.
490: Methods can be chosen with -ts_ssp_type {rks2,rks3,rk104}
492: rks2: Second order methods with any number s>1 of stages. c_eff = (s-1)/s
494: rks3: Third order methods with s=n^2 stages, n>1. c_eff = (s-n)/s
496: rk104: A 10-stage fourth order method. c_eff = 0.6
498: Level: beginner
500: References:
501: + 1. - Ketcheson, Highly efficient strong stability preserving Runge Kutta methods with low storage implementations, SISC, 2008.
502: - 2. - Gottlieb, Ketcheson, and Shu, High order strong stability preserving time discretizations, J Scientific Computing, 2009.
504: .seealso: TSCreate(), TS, TSSetType()
506: M*/
509: PETSC_EXTERN PetscErrorCode TSCreate_SSP(TS ts)
510: {
511: TS_SSP *ssp;
515: TSSSPInitializePackage();
517: ts->ops->setup = TSSetUp_SSP;
518: ts->ops->step = TSStep_SSP;
519: ts->ops->reset = TSReset_SSP;
520: ts->ops->destroy = TSDestroy_SSP;
521: ts->ops->setfromoptions = TSSetFromOptions_SSP;
522: ts->ops->view = TSView_SSP;
524: PetscNewLog(ts,&ssp);
525: ts->data = (void*)ssp;
527: PetscObjectComposeFunction((PetscObject)ts,"TSSSPGetType_C",TSSSPGetType_SSP);
528: PetscObjectComposeFunction((PetscObject)ts,"TSSSPSetType_C",TSSSPSetType_SSP);
529: PetscObjectComposeFunction((PetscObject)ts,"TSSSPGetNumStages_C",TSSSPGetNumStages_SSP);
530: PetscObjectComposeFunction((PetscObject)ts,"TSSSPSetNumStages_C",TSSSPSetNumStages_SSP);
532: TSSSPSetType(ts,TSSSPRKS2);
533: ssp->nstages = 5;
534: return(0);
535: }
539: /*@C
540: TSSSPInitializePackage - This function initializes everything in the TSSSP package. It is called
541: from PetscDLLibraryRegister() when using dynamic libraries, and on the first call to TSCreate_SSP()
542: when using static libraries.
544: Level: developer
546: .keywords: TS, TSSSP, initialize, package
547: .seealso: PetscInitialize()
548: @*/
549: PetscErrorCode TSSSPInitializePackage(void)
550: {
554: if (TSSSPPackageInitialized) return(0);
555: TSSSPPackageInitialized = PETSC_TRUE;
556: PetscFunctionListAdd(&TSSSPList,TSSSPRKS2, TSSSPStep_RK_2);
557: PetscFunctionListAdd(&TSSSPList,TSSSPRKS3, TSSSPStep_RK_3);
558: PetscFunctionListAdd(&TSSSPList,TSSSPRK104,TSSSPStep_RK_10_4);
559: PetscRegisterFinalize(TSSSPFinalizePackage);
560: return(0);
561: }
565: /*@C
566: TSSSPFinalizePackage - This function destroys everything in the TSSSP package. It is
567: called from PetscFinalize().
569: Level: developer
571: .keywords: Petsc, destroy, package
572: .seealso: PetscFinalize()
573: @*/
574: PetscErrorCode TSSSPFinalizePackage(void)
575: {
579: TSSSPPackageInitialized = PETSC_FALSE;
580: PetscFunctionListDestroy(&TSSSPList);
581: return(0);
582: }