Actual source code: ssp.c

petsc-3.7.6 2017-04-24
Report Typos and Errors
  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: }