Actual source code: ssp.c
petsc-3.7.5 2017-01-01
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: TSGetAdapt(ts,&ts->adapt);
206: TSAdaptCandidatesClear(ts->adapt);
207: return(0);
208: }
212: static PetscErrorCode TSStep_SSP(TS ts)
213: {
214: TS_SSP *ssp = (TS_SSP*)ts->data;
215: Vec sol = ts->vec_sol;
216: PetscBool stageok,accept = PETSC_TRUE;
217: PetscReal next_time_step = ts->time_step;
221: (*ssp->onestep)(ts,ts->ptime,ts->time_step,sol);
222: TSPostStage(ts,ts->ptime,0,&sol);
223: TSAdaptCheckStage(ts->adapt,ts,ts->ptime+ts->time_step,sol,&stageok);
224: if(!stageok) {ts->reason = TS_DIVERGED_STEP_REJECTED; return(0);}
226: TSAdaptChoose(ts->adapt,ts,ts->time_step,NULL,&next_time_step,&accept);
227: if (!accept) {ts->reason = TS_DIVERGED_STEP_REJECTED; return(0);}
229: ts->ptime += ts->time_step;
230: ts->time_step = next_time_step;
231: return(0);
232: }
233: /*------------------------------------------------------------*/
237: static PetscErrorCode TSReset_SSP(TS ts)
238: {
239: TS_SSP *ssp = (TS_SSP*)ts->data;
243: if (ssp->work) {VecDestroyVecs(ssp->nwork,&ssp->work);}
244: ssp->nwork = 0;
245: ssp->workout = PETSC_FALSE;
246: return(0);
247: }
251: static PetscErrorCode TSDestroy_SSP(TS ts)
252: {
253: TS_SSP *ssp = (TS_SSP*)ts->data;
257: TSReset_SSP(ts);
258: PetscFree(ssp->type_name);
259: PetscFree(ts->data);
260: PetscObjectComposeFunction((PetscObject)ts,"TSSSPGetType_C",NULL);
261: PetscObjectComposeFunction((PetscObject)ts,"TSSSPSetType_C",NULL);
262: PetscObjectComposeFunction((PetscObject)ts,"TSSSPGetNumStages_C",NULL);
263: PetscObjectComposeFunction((PetscObject)ts,"TSSSPSetNumStages_C",NULL);
264: return(0);
265: }
266: /*------------------------------------------------------------*/
270: /*@C
271: TSSSPSetType - set the SSP time integration scheme to use
273: Logically Collective
275: Input Arguments:
276: ts - time stepping object
277: type - type of scheme to use
279: Options Database Keys:
280: -ts_ssp_type <rks2>: Type of SSP method (one of) rks2 rks3 rk104
281: -ts_ssp_nstages <5>: Number of stages
283: Level: beginner
285: .seealso: TSSSP, TSSSPGetType(), TSSSPSetNumStages(), TSSSPRKS2, TSSSPRKS3, TSSSPRK104
286: @*/
287: PetscErrorCode TSSSPSetType(TS ts,TSSSPType type)
288: {
293: PetscTryMethod(ts,"TSSSPSetType_C",(TS,TSSSPType),(ts,type));
294: return(0);
295: }
299: /*@C
300: TSSSPGetType - get the SSP time integration scheme
302: Logically Collective
304: Input Argument:
305: ts - time stepping object
307: Output Argument:
308: type - type of scheme being used
310: Level: beginner
312: .seealso: TSSSP, TSSSPSettype(), TSSSPSetNumStages(), TSSSPRKS2, TSSSPRKS3, TSSSPRK104
313: @*/
314: PetscErrorCode TSSSPGetType(TS ts,TSSSPType *type)
315: {
320: PetscUseMethod(ts,"TSSSPGetType_C",(TS,TSSSPType*),(ts,type));
321: return(0);
322: }
326: /*@
327: TSSSPSetNumStages - set the number of stages to use with the SSP method
329: Logically Collective
331: Input Arguments:
332: ts - time stepping object
333: nstages - number of stages
335: Options Database Keys:
336: -ts_ssp_type <rks2>: NumStages of SSP method (one of) rks2 rks3 rk104
337: -ts_ssp_nstages <5>: Number of stages
339: Level: beginner
341: .seealso: TSSSP, TSSSPGetNumStages(), TSSSPSetNumStages(), TSSSPRKS2, TSSSPRKS3, TSSSPRK104
342: @*/
343: PetscErrorCode TSSSPSetNumStages(TS ts,PetscInt nstages)
344: {
349: PetscTryMethod(ts,"TSSSPSetNumStages_C",(TS,PetscInt),(ts,nstages));
350: return(0);
351: }
355: /*@
356: TSSSPGetNumStages - get the number of stages in the SSP time integration scheme
358: Logically Collective
360: Input Argument:
361: ts - time stepping object
363: Output Argument:
364: nstages - number of stages
366: Level: beginner
368: .seealso: TSSSP, TSSSPGetType(), TSSSPSetNumStages(), TSSSPRKS2, TSSSPRKS3, TSSSPRK104
369: @*/
370: PetscErrorCode TSSSPGetNumStages(TS ts,PetscInt *nstages)
371: {
376: PetscUseMethod(ts,"TSSSPGetNumStages_C",(TS,PetscInt*),(ts,nstages));
377: return(0);
378: }
382: static PetscErrorCode TSSSPSetType_SSP(TS ts,TSSSPType type)
383: {
384: PetscErrorCode ierr,(*r)(TS,PetscReal,PetscReal,Vec);
385: TS_SSP *ssp = (TS_SSP*)ts->data;
388: PetscFunctionListFind(TSSSPList,type,&r);
389: if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown TS_SSP type %s given",type);
390: ssp->onestep = r;
391: PetscFree(ssp->type_name);
392: PetscStrallocpy(type,&ssp->type_name);
393: return(0);
394: }
397: static PetscErrorCode TSSSPGetType_SSP(TS ts,TSSSPType *type)
398: {
399: TS_SSP *ssp = (TS_SSP*)ts->data;
402: *type = ssp->type_name;
403: return(0);
404: }
407: static PetscErrorCode TSSSPSetNumStages_SSP(TS ts,PetscInt nstages)
408: {
409: TS_SSP *ssp = (TS_SSP*)ts->data;
412: ssp->nstages = nstages;
413: return(0);
414: }
417: static PetscErrorCode TSSSPGetNumStages_SSP(TS ts,PetscInt *nstages)
418: {
419: TS_SSP *ssp = (TS_SSP*)ts->data;
422: *nstages = ssp->nstages;
423: return(0);
424: }
428: static PetscErrorCode TSSetFromOptions_SSP(PetscOptionItems *PetscOptionsObject,TS ts)
429: {
430: char tname[256] = TSSSPRKS2;
431: TS_SSP *ssp = (TS_SSP*)ts->data;
433: PetscBool flg;
436: PetscOptionsHead(PetscOptionsObject,"SSP ODE solver options");
437: {
438: PetscOptionsFList("-ts_ssp_type","Type of SSP method","TSSSPSetType",TSSSPList,tname,tname,sizeof(tname),&flg);
439: if (flg) {
440: TSSSPSetType(ts,tname);
441: }
442: PetscOptionsInt("-ts_ssp_nstages","Number of stages","TSSSPSetNumStages",ssp->nstages,&ssp->nstages,NULL);
443: }
444: PetscOptionsTail();
445: return(0);
446: }
450: static PetscErrorCode TSView_SSP(TS ts,PetscViewer viewer)
451: {
452: TS_SSP *ssp = (TS_SSP*)ts->data;
453: PetscBool ascii;
457: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&ascii);
458: if (ascii) {PetscViewerASCIIPrintf(viewer," Scheme: %s\n",ssp->type_name);}
459: return(0);
460: }
462: /* ------------------------------------------------------------ */
464: /*MC
465: TSSSP - Explicit strong stability preserving ODE solver
467: Most hyperbolic conservation laws have exact solutions that are total variation diminishing (TVD) or total variation
468: bounded (TVB) although these solutions often contain discontinuities. Spatial discretizations such as Godunov's
469: scheme and high-resolution finite volume methods (TVD limiters, ENO/WENO) are designed to preserve these properties,
470: but they are usually formulated using a forward Euler time discretization or by coupling the space and time
471: discretization as in the classical Lax-Wendroff scheme. When the space and time discretization is coupled, it is very
472: difficult to produce schemes with high temporal accuracy while preserving TVD properties. An alternative is the
473: semidiscrete formulation where we choose a spatial discretization that is TVD with forward Euler and then choose a
474: time discretization that preserves the TVD property. Such integrators are called strong stability preserving (SSP).
476: Let c_eff be the minimum number of function evaluations required to step as far as one step of forward Euler while
477: still being SSP. Some theoretical bounds
479: 1. There are no explicit methods with c_eff > 1.
481: 2. There are no explicit methods beyond order 4 (for nonlinear problems) and c_eff > 0.
483: 3. There are no implicit methods with order greater than 1 and c_eff > 2.
485: This integrator provides Runge-Kutta methods of order 2, 3, and 4 with maximal values of c_eff. More stages allows
486: for larger values of c_eff which improves efficiency. These implementations are low-memory and only use 2 or 3 work
487: vectors regardless of the total number of stages, so e.g. 25-stage 3rd order methods may be an excellent choice.
489: Methods can be chosen with -ts_ssp_type {rks2,rks3,rk104}
491: rks2: Second order methods with any number s>1 of stages. c_eff = (s-1)/s
493: rks3: Third order methods with s=n^2 stages, n>1. c_eff = (s-n)/s
495: rk104: A 10-stage fourth order method. c_eff = 0.6
497: Level: beginner
499: References:
500: + 1. - Ketcheson, Highly efficient strong stability preserving Runge Kutta methods with low storage implementations, SISC, 2008.
501: - 2. - Gottlieb, Ketcheson, and Shu, High order strong stability preserving time discretizations, J Scientific Computing, 2009.
503: .seealso: TSCreate(), TS, TSSetType()
505: M*/
508: PETSC_EXTERN PetscErrorCode TSCreate_SSP(TS ts)
509: {
510: TS_SSP *ssp;
514: TSSSPInitializePackage();
516: ts->ops->setup = TSSetUp_SSP;
517: ts->ops->step = TSStep_SSP;
518: ts->ops->reset = TSReset_SSP;
519: ts->ops->destroy = TSDestroy_SSP;
520: ts->ops->setfromoptions = TSSetFromOptions_SSP;
521: ts->ops->view = TSView_SSP;
523: PetscNewLog(ts,&ssp);
524: ts->data = (void*)ssp;
526: PetscObjectComposeFunction((PetscObject)ts,"TSSSPGetType_C",TSSSPGetType_SSP);
527: PetscObjectComposeFunction((PetscObject)ts,"TSSSPSetType_C",TSSSPSetType_SSP);
528: PetscObjectComposeFunction((PetscObject)ts,"TSSSPGetNumStages_C",TSSSPGetNumStages_SSP);
529: PetscObjectComposeFunction((PetscObject)ts,"TSSSPSetNumStages_C",TSSSPSetNumStages_SSP);
531: TSSSPSetType(ts,TSSSPRKS2);
532: ssp->nstages = 5;
533: return(0);
534: }
538: /*@C
539: TSSSPInitializePackage - This function initializes everything in the TSSSP package. It is called
540: from PetscDLLibraryRegister() when using dynamic libraries, and on the first call to TSCreate_SSP()
541: when using static libraries.
543: Level: developer
545: .keywords: TS, TSSSP, initialize, package
546: .seealso: PetscInitialize()
547: @*/
548: PetscErrorCode TSSSPInitializePackage(void)
549: {
553: if (TSSSPPackageInitialized) return(0);
554: TSSSPPackageInitialized = PETSC_TRUE;
555: PetscFunctionListAdd(&TSSSPList,TSSSPRKS2, TSSSPStep_RK_2);
556: PetscFunctionListAdd(&TSSSPList,TSSSPRKS3, TSSSPStep_RK_3);
557: PetscFunctionListAdd(&TSSSPList,TSSSPRK104,TSSSPStep_RK_10_4);
558: PetscRegisterFinalize(TSSSPFinalizePackage);
559: return(0);
560: }
564: /*@C
565: TSSSPFinalizePackage - This function destroys everything in the TSSSP package. It is
566: called from PetscFinalize().
568: Level: developer
570: .keywords: Petsc, destroy, package
571: .seealso: PetscFinalize()
572: @*/
573: PetscErrorCode TSSSPFinalizePackage(void)
574: {
578: TSSSPPackageInitialized = PETSC_FALSE;
579: PetscFunctionListDestroy(&TSSSPList);
580: return(0);
581: }