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