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;

 24:   if (ssp->workout) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Work vectors already gotten");
 25:   if (ssp->nwork < n) {
 26:     if (ssp->nwork > 0) {
 27:       VecDestroyVecs(ssp->nwork,&ssp->work);
 28:     }
 29:     VecDuplicateVecs(ts->vec_sol,n,&ssp->work);
 30:     ssp->nwork = n;
 31:   }
 32:   *work = ssp->work;
 33:   ssp->workout = PETSC_TRUE;
 34:   return(0);
 35: }

 37: static PetscErrorCode TSSSPRestoreWorkVectors(TS ts,PetscInt n,Vec **work)
 38: {
 39:   TS_SSP *ssp = (TS_SSP*)ts->data;

 42:   if (!ssp->workout) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Work vectors have not been gotten");
 43:   if (*work != ssp->work) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Wrong work vectors checked out");
 44:   ssp->workout = PETSC_FALSE;
 45:   *work = NULL;
 46:   return(0);
 47: }

 49: /*MC
 50:    TSSSPRKS2 - Optimal second order SSP Runge-Kutta method, low-storage, c_eff=(s-1)/s

 52:    Pseudocode 2 of Ketcheson 2008

 54:    Level: beginner

 56: .seealso: TSSSP, TSSSPSetType(), TSSSPSetNumStages()
 57: M*/
 58: static PetscErrorCode TSSSPStep_RK_2(TS ts,PetscReal t0,PetscReal dt,Vec sol)
 59: {
 60:   TS_SSP         *ssp = (TS_SSP*)ts->data;
 61:   Vec            *work,F;
 62:   PetscInt       i,s;

 66:   s    = ssp->nstages;
 67:   TSSSPGetWorkVectors(ts,2,&work);
 68:   F    = work[1];
 69:   VecCopy(sol,work[0]);
 70:   for (i=0; i<s-1; i++) {
 71:     PetscReal stage_time = t0+dt*(i/(s-1.));
 72:     TSPreStage(ts,stage_time);
 73:     TSComputeRHSFunction(ts,stage_time,work[0],F);
 74:     VecAXPY(work[0],dt/(s-1.),F);
 75:   }
 76:   TSComputeRHSFunction(ts,t0+dt,work[0],F);
 77:   VecAXPBYPCZ(sol,(s-1.)/s,dt/s,1./s,work[0],F);
 78:   TSSSPRestoreWorkVectors(ts,2,&work);
 79:   return(0);
 80: }

 82: /*MC
 83:    TSSSPRKS3 - Optimal third order SSP Runge-Kutta, low-storage, c_eff=(PetscSqrtReal(s)-1)/PetscSqrtReal(s), where PetscSqrtReal(s) is an integer

 85:    Pseudocode 2 of Ketcheson 2008

 87:    Level: beginner

 89: .seealso: TSSSP, TSSSPSetType(), TSSSPSetNumStages()
 90: M*/
 91: static PetscErrorCode TSSSPStep_RK_3(TS ts,PetscReal t0,PetscReal dt,Vec sol)
 92: {
 93:   TS_SSP         *ssp = (TS_SSP*)ts->data;
 94:   Vec            *work,F;
 95:   PetscInt       i,s,n,r;
 96:   PetscReal      c,stage_time;

100:   s = ssp->nstages;
101:   n = (PetscInt)(PetscSqrtReal((PetscReal)s)+0.001);
102:   r = s-n;
103:   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);
104:   TSSSPGetWorkVectors(ts,3,&work);
105:   F    = work[2];
106:   VecCopy(sol,work[0]);
107:   for (i=0; i<(n-1)*(n-2)/2; i++) {
108:     c          = (i<n*(n+1)/2) ? 1.*i/(s-n) : (1.*i-n)/(s-n);
109:     stage_time = t0+c*dt;
110:     TSPreStage(ts,stage_time);
111:     TSComputeRHSFunction(ts,stage_time,work[0],F);
112:     VecAXPY(work[0],dt/r,F);
113:   }
114:   VecCopy(work[0],work[1]);
115:   for (; i<n*(n+1)/2-1; i++) {
116:     c          = (i<n*(n+1)/2) ? 1.*i/(s-n) : (1.*i-n)/(s-n);
117:     stage_time = t0+c*dt;
118:     TSPreStage(ts,stage_time);
119:     TSComputeRHSFunction(ts,stage_time,work[0],F);
120:     VecAXPY(work[0],dt/r,F);
121:   }
122:   {
123:     c          = (i<n*(n+1)/2) ? 1.*i/(s-n) : (1.*i-n)/(s-n);
124:     stage_time = t0+c*dt;
125:     TSPreStage(ts,stage_time);
126:     TSComputeRHSFunction(ts,stage_time,work[0],F);
127:     VecAXPBYPCZ(work[0],1.*n/(2*n-1.),(n-1.)*dt/(r*(2*n-1)),(n-1.)/(2*n-1.),work[1],F);
128:     i++;
129:   }
130:   for (; i<s; i++) {
131:     c          = (i<n*(n+1)/2) ? 1.*i/(s-n) : (1.*i-n)/(s-n);
132:     stage_time = t0+c*dt;
133:     TSPreStage(ts,stage_time);
134:     TSComputeRHSFunction(ts,stage_time,work[0],F);
135:     VecAXPY(work[0],dt/r,F);
136:   }
137:   VecCopy(work[0],sol);
138:   TSSSPRestoreWorkVectors(ts,3,&work);
139:   return(0);
140: }

142: /*MC
143:    TSSSPRKS104 - Optimal fourth order SSP Runge-Kutta, low-storage (2N), c_eff=0.6

145:    SSPRK(10,4), Pseudocode 3 of Ketcheson 2008

147:    Level: beginner

149: .seealso: TSSSP, TSSSPSetType()
150: M*/
151: static PetscErrorCode TSSSPStep_RK_10_4(TS ts,PetscReal t0,PetscReal dt,Vec sol)
152: {
153:   const PetscReal c[10] = {0, 1./6, 2./6, 3./6, 4./6, 2./6, 3./6, 4./6, 5./6, 1};
154:   Vec             *work,F;
155:   PetscInt        i;
156:   PetscReal       stage_time;
157:   PetscErrorCode  ierr;

160:   TSSSPGetWorkVectors(ts,3,&work);
161:   F    = work[2];
162:   VecCopy(sol,work[0]);
163:   for (i=0; i<5; i++) {
164:     stage_time = t0+c[i]*dt;
165:     TSPreStage(ts,stage_time);
166:     TSComputeRHSFunction(ts,stage_time,work[0],F);
167:     VecAXPY(work[0],dt/6,F);
168:   }
169:   VecAXPBYPCZ(work[1],1./25,9./25,0,sol,work[0]);
170:   VecAXPBY(work[0],15,-5,work[1]);
171:   for (; i<9; i++) {
172:     stage_time = t0+c[i]*dt;
173:     TSPreStage(ts,stage_time);
174:     TSComputeRHSFunction(ts,stage_time,work[0],F);
175:     VecAXPY(work[0],dt/6,F);
176:   }
177:   stage_time = t0+dt;
178:   TSPreStage(ts,stage_time);
179:   TSComputeRHSFunction(ts,stage_time,work[0],F);
180:   VecAXPBYPCZ(work[1],3./5,dt/10,1,work[0],F);
181:   VecCopy(work[1],sol);
182:   TSSSPRestoreWorkVectors(ts,3,&work);
183:   return(0);
184: }

186: static PetscErrorCode TSSetUp_SSP(TS ts)
187: {

191:   TSCheckImplicitTerm(ts);
192:   TSGetAdapt(ts,&ts->adapt);
193:   TSAdaptCandidatesClear(ts->adapt);
194:   return(0);
195: }

197: static PetscErrorCode TSStep_SSP(TS ts)
198: {
199:   TS_SSP         *ssp = (TS_SSP*)ts->data;
200:   Vec            sol  = ts->vec_sol;
201:   PetscBool      stageok,accept = PETSC_TRUE;
202:   PetscReal      next_time_step = ts->time_step;

206:   (*ssp->onestep)(ts,ts->ptime,ts->time_step,sol);
207:   TSPostStage(ts,ts->ptime,0,&sol);
208:   TSAdaptCheckStage(ts->adapt,ts,ts->ptime+ts->time_step,sol,&stageok);
209:   if (!stageok) {ts->reason = TS_DIVERGED_STEP_REJECTED; return(0);}

211:   TSAdaptChoose(ts->adapt,ts,ts->time_step,NULL,&next_time_step,&accept);
212:   if (!accept) {ts->reason = TS_DIVERGED_STEP_REJECTED; return(0);}

214:   ts->ptime += ts->time_step;
215:   ts->time_step = next_time_step;
216:   return(0);
217: }
218: /*------------------------------------------------------------*/

220: static PetscErrorCode TSReset_SSP(TS ts)
221: {
222:   TS_SSP         *ssp = (TS_SSP*)ts->data;

226:   if (ssp->work) {VecDestroyVecs(ssp->nwork,&ssp->work);}
227:   ssp->nwork   = 0;
228:   ssp->workout = PETSC_FALSE;
229:   return(0);
230: }

232: static PetscErrorCode TSDestroy_SSP(TS ts)
233: {
234:   TS_SSP         *ssp = (TS_SSP*)ts->data;

238:   TSReset_SSP(ts);
239:   PetscFree(ssp->type_name);
240:   PetscFree(ts->data);
241:   PetscObjectComposeFunction((PetscObject)ts,"TSSSPGetType_C",NULL);
242:   PetscObjectComposeFunction((PetscObject)ts,"TSSSPSetType_C",NULL);
243:   PetscObjectComposeFunction((PetscObject)ts,"TSSSPGetNumStages_C",NULL);
244:   PetscObjectComposeFunction((PetscObject)ts,"TSSSPSetNumStages_C",NULL);
245:   return(0);
246: }
247: /*------------------------------------------------------------*/

249: /*@C
250:    TSSSPSetType - set the SSP time integration scheme to use

252:    Logically Collective

254:    Input Parameters:
255: +  ts - time stepping object
256: -  ssptype - type of scheme to use

258:    Options Database Keys:
259:    -ts_ssp_type <rks2>: Type of SSP method (one of) rks2 rks3 rk104
260:    -ts_ssp_nstages <5>: Number of stages

262:    Level: beginner

264: .seealso: TSSSP, TSSSPGetType(), TSSSPSetNumStages(), TSSSPRKS2, TSSSPRKS3, TSSSPRK104
265: @*/
266: PetscErrorCode TSSSPSetType(TS ts,TSSSPType ssptype)
267: {

273:   PetscTryMethod(ts,"TSSSPSetType_C",(TS,TSSSPType),(ts,ssptype));
274:   return(0);
275: }

277: /*@C
278:    TSSSPGetType - get the SSP time integration scheme

280:    Logically Collective

282:    Input Parameter:
283: .  ts - time stepping object

285:    Output Parameter:
286: .  type - type of scheme being used

288:    Level: beginner

290: .seealso: TSSSP, TSSSPSettype(), TSSSPSetNumStages(), TSSSPRKS2, TSSSPRKS3, TSSSPRK104
291: @*/
292: PetscErrorCode TSSSPGetType(TS ts,TSSSPType *type)
293: {

298:   PetscUseMethod(ts,"TSSSPGetType_C",(TS,TSSSPType*),(ts,type));
299:   return(0);
300: }

302: /*@
303:    TSSSPSetNumStages - set the number of stages to use with the SSP method

305:    Logically Collective

307:    Input Parameters:
308: +  ts - time stepping object
309: -  nstages - number of stages

311:    Options Database Keys:
312:    -ts_ssp_type <rks2>: NumStages of SSP method (one of) rks2 rks3 rk104
313:    -ts_ssp_nstages <5>: Number of stages

315:    Level: beginner

317: .seealso: TSSSP, TSSSPGetNumStages(), TSSSPSetNumStages(), TSSSPRKS2, TSSSPRKS3, TSSSPRK104
318: @*/
319: PetscErrorCode TSSSPSetNumStages(TS ts,PetscInt nstages)
320: {

325:   PetscTryMethod(ts,"TSSSPSetNumStages_C",(TS,PetscInt),(ts,nstages));
326:   return(0);
327: }

329: /*@
330:    TSSSPGetNumStages - get the number of stages in the SSP time integration scheme

332:    Logically Collective

334:    Input Parameter:
335: .  ts - time stepping object

337:    Output Parameter:
338: .  nstages - number of stages

340:    Level: beginner

342: .seealso: TSSSP, TSSSPGetType(), TSSSPSetNumStages(), TSSSPRKS2, TSSSPRKS3, TSSSPRK104
343: @*/
344: PetscErrorCode TSSSPGetNumStages(TS ts,PetscInt *nstages)
345: {

350:   PetscUseMethod(ts,"TSSSPGetNumStages_C",(TS,PetscInt*),(ts,nstages));
351:   return(0);
352: }

354: static PetscErrorCode TSSSPSetType_SSP(TS ts,TSSSPType type)
355: {
356:   PetscErrorCode ierr,(*r)(TS,PetscReal,PetscReal,Vec);
357:   TS_SSP         *ssp = (TS_SSP*)ts->data;

360:   PetscFunctionListFind(TSSSPList,type,&r);
361:   if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown TS_SSP type %s given",type);
362:   ssp->onestep = r;
363:   PetscFree(ssp->type_name);
364:   PetscStrallocpy(type,&ssp->type_name);
365:   ts->default_adapt_type = TSADAPTNONE;
366:   return(0);
367: }
368: static PetscErrorCode TSSSPGetType_SSP(TS ts,TSSSPType *type)
369: {
370:   TS_SSP *ssp = (TS_SSP*)ts->data;

373:   *type = ssp->type_name;
374:   return(0);
375: }
376: static PetscErrorCode TSSSPSetNumStages_SSP(TS ts,PetscInt nstages)
377: {
378:   TS_SSP *ssp = (TS_SSP*)ts->data;

381:   ssp->nstages = nstages;
382:   return(0);
383: }
384: static PetscErrorCode TSSSPGetNumStages_SSP(TS ts,PetscInt *nstages)
385: {
386:   TS_SSP *ssp = (TS_SSP*)ts->data;

389:   *nstages = ssp->nstages;
390:   return(0);
391: }

393: static PetscErrorCode TSSetFromOptions_SSP(PetscOptionItems *PetscOptionsObject,TS ts)
394: {
395:   char           tname[256] = TSSSPRKS2;
396:   TS_SSP         *ssp       = (TS_SSP*)ts->data;
398:   PetscBool      flg;

401:   PetscOptionsHead(PetscOptionsObject,"SSP ODE solver options");
402:   {
403:     PetscOptionsFList("-ts_ssp_type","Type of SSP method","TSSSPSetType",TSSSPList,tname,tname,sizeof(tname),&flg);
404:     if (flg) {
405:       TSSSPSetType(ts,tname);
406:     }
407:     PetscOptionsInt("-ts_ssp_nstages","Number of stages","TSSSPSetNumStages",ssp->nstages,&ssp->nstages,NULL);
408:   }
409:   PetscOptionsTail();
410:   return(0);
411: }

413: static PetscErrorCode TSView_SSP(TS ts,PetscViewer viewer)
414: {
415:   TS_SSP         *ssp = (TS_SSP*)ts->data;
416:   PetscBool      ascii;

420:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&ascii);
421:   if (ascii) {PetscViewerASCIIPrintf(viewer,"  Scheme: %s\n",ssp->type_name);}
422:   return(0);
423: }

425: /* ------------------------------------------------------------ */

427: /*MC
428:       TSSSP - Explicit strong stability preserving ODE solver

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).

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

442:   1. There are no explicit methods with c_eff > 1.

444:   2. There are no explicit methods beyond order 4 (for nonlinear problems) and c_eff > 0.

446:   3. There are no implicit methods with order greater than 1 and c_eff > 2.

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.

452:   Methods can be chosen with -ts_ssp_type {rks2,rks3,rk104}

454:   rks2: Second order methods with any number s>1 of stages.  c_eff = (s-1)/s

456:   rks3: Third order methods with s=n^2 stages, n>1.  c_eff = (s-n)/s

458:   rk104: A 10-stage fourth order method.  c_eff = 0.6

460:   Level: beginner

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.

466: .seealso:  TSCreate(), TS, TSSetType()

468: M*/
469: PETSC_EXTERN PetscErrorCode TSCreate_SSP(TS ts)
470: {
471:   TS_SSP         *ssp;

475:   TSSSPInitializePackage();

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;

484:   PetscNewLog(ts,&ssp);
485:   ts->data = (void*)ssp;

487:   PetscObjectComposeFunction((PetscObject)ts,"TSSSPGetType_C",TSSSPGetType_SSP);
488:   PetscObjectComposeFunction((PetscObject)ts,"TSSSPSetType_C",TSSSPSetType_SSP);
489:   PetscObjectComposeFunction((PetscObject)ts,"TSSSPGetNumStages_C",TSSSPGetNumStages_SSP);
490:   PetscObjectComposeFunction((PetscObject)ts,"TSSSPSetNumStages_C",TSSSPSetNumStages_SSP);

492:   TSSSPSetType(ts,TSSSPRKS2);
493:   ssp->nstages = 5;
494:   return(0);
495: }

497: /*@C
498:   TSSSPInitializePackage - This function initializes everything in the TSSSP package. It is called
499:   from TSInitializePackage().

501:   Level: developer

503: .seealso: PetscInitialize()
504: @*/
505: PetscErrorCode TSSSPInitializePackage(void)
506: {

510:   if (TSSSPPackageInitialized) return(0);
511:   TSSSPPackageInitialized = PETSC_TRUE;
512:   PetscFunctionListAdd(&TSSSPList,TSSSPRKS2, TSSSPStep_RK_2);
513:   PetscFunctionListAdd(&TSSSPList,TSSSPRKS3, TSSSPStep_RK_3);
514:   PetscFunctionListAdd(&TSSSPList,TSSSPRK104,TSSSPStep_RK_10_4);
515:   PetscRegisterFinalize(TSSSPFinalizePackage);
516:   return(0);
517: }

519: /*@C
520:   TSSSPFinalizePackage - This function destroys everything in the TSSSP package. It is
521:   called from PetscFinalize().

523:   Level: developer

525: .seealso: PetscFinalize()
526: @*/
527: PetscErrorCode TSSSPFinalizePackage(void)
528: {

532:   TSSSPPackageInitialized = PETSC_FALSE;
533:   PetscFunctionListDestroy(&TSSSPList);
534:   return(0);
535: }