Actual source code: ssp.c

petsc-3.8.3 2017-12-09
Report Typos and Errors
  1: /*
  2:        Code for Timestepping with explicit SSP.
  3: */
  4:  #include <petsc/private/tsimpl.h>

  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;


 19: static PetscErrorCode TSSSPGetWorkVectors(TS ts,PetscInt n,Vec **work)
 20: {
 21:   TS_SSP         *ssp = (TS_SSP*)ts->data;

 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:       VecDestroyVecs(ssp->nwork,&ssp->work);
 29:     }
 30:     VecDuplicateVecs(ts->vec_sol,n,&ssp->work);
 31:     ssp->nwork = n;
 32:   }
 33:   *work = ssp->work;
 34:   ssp->workout = PETSC_TRUE;
 35:   return(0);
 36: }

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

 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:   return(0);
 48: }

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

 53:    Pseudocode 2 of Ketcheson 2008

 55:    Level: beginner

 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;

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

 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

 86:    Pseudocode 2 of Ketcheson 2008

 88:    Level: beginner

 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;

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:   TSSSPGetWorkVectors(ts,3,&work);
106:   F    = work[2];
107:   VecCopy(sol,work[0]);
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:     TSPreStage(ts,stage_time);
112:     TSComputeRHSFunction(ts,stage_time,work[0],F);
113:     VecAXPY(work[0],dt/r,F);
114:   }
115:   VecCopy(work[0],work[1]);
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:     TSPreStage(ts,stage_time);
120:     TSComputeRHSFunction(ts,stage_time,work[0],F);
121:     VecAXPY(work[0],dt/r,F);
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:     TSPreStage(ts,stage_time);
127:     TSComputeRHSFunction(ts,stage_time,work[0],F);
128:     VecAXPBYPCZ(work[0],1.*n/(2*n-1.),(n-1.)*dt/(r*(2*n-1)),(n-1.)/(2*n-1.),work[1],F);
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:     TSPreStage(ts,stage_time);
135:     TSComputeRHSFunction(ts,stage_time,work[0],F);
136:     VecAXPY(work[0],dt/r,F);
137:   }
138:   VecCopy(work[0],sol);
139:   TSSSPRestoreWorkVectors(ts,3,&work);
140:   return(0);
141: }

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

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

148:    Level: beginner

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;

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


188: static PetscErrorCode TSSetUp_SSP(TS ts)
189: {

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

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;

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

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

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

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

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

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

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

251: /*@C
252:    TSSSPSetType - set the SSP time integration scheme to use

254:    Logically Collective

256:    Input Arguments:
257:    ts - time stepping object
258:    ssptype - type of scheme to use

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

264:    Level: beginner

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

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

279: /*@C
280:    TSSSPGetType - get the SSP time integration scheme

282:    Logically Collective

284:    Input Argument:
285:    ts - time stepping object

287:    Output Argument:
288:    type - type of scheme being used

290:    Level: beginner

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

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

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

307:    Logically Collective

309:    Input Arguments:
310:    ts - time stepping object
311:    nstages - number of stages

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

317:    Level: beginner

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

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

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

334:    Logically Collective

336:    Input Argument:
337:    ts - time stepping object

339:    Output Argument:
340:    nstages - number of stages

342:    Level: beginner

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

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

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

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

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

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

391:   *nstages = ssp->nstages;
392:   return(0);
393: }

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

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

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

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

427: /* ------------------------------------------------------------ */

429: /*MC
430:       TSSSP - Explicit strong stability preserving ODE solver

432:   Most hyperbolic conservation laws have exact solutions that are total variation diminishing (TVD) or total variation
433:   bounded (TVB) although these solutions often contain discontinuities.  Spatial discretizations such as Godunov's
434:   scheme and high-resolution finite volume methods (TVD limiters, ENO/WENO) are designed to preserve these properties,
435:   but they are usually formulated using a forward Euler time discretization or by coupling the space and time
436:   discretization as in the classical Lax-Wendroff scheme.  When the space and time discretization is coupled, it is very
437:   difficult to produce schemes with high temporal accuracy while preserving TVD properties.  An alternative is the
438:   semidiscrete formulation where we choose a spatial discretization that is TVD with forward Euler and then choose a
439:   time discretization that preserves the TVD property.  Such integrators are called strong stability preserving (SSP).

441:   Let c_eff be the minimum number of function evaluations required to step as far as one step of forward Euler while
442:   still being SSP.  Some theoretical bounds

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

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

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

450:   This integrator provides Runge-Kutta methods of order 2, 3, and 4 with maximal values of c_eff.  More stages allows
451:   for larger values of c_eff which improves efficiency.  These implementations are low-memory and only use 2 or 3 work
452:   vectors regardless of the total number of stages, so e.g. 25-stage 3rd order methods may be an excellent choice.

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

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

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

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

462:   Level: beginner

464:   References:
465: +  1. - Ketcheson, Highly efficient strong stability preserving Runge Kutta methods with low storage implementations, SISC, 2008.
466: -  2. - Gottlieb, Ketcheson, and Shu, High order strong stability preserving time discretizations, J Scientific Computing, 2009.

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

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

477:   TSSSPInitializePackage();

479:   ts->ops->setup          = TSSetUp_SSP;
480:   ts->ops->step           = TSStep_SSP;
481:   ts->ops->reset          = TSReset_SSP;
482:   ts->ops->destroy        = TSDestroy_SSP;
483:   ts->ops->setfromoptions = TSSetFromOptions_SSP;
484:   ts->ops->view           = TSView_SSP;

486:   PetscNewLog(ts,&ssp);
487:   ts->data = (void*)ssp;

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

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

499: /*@C
500:   TSSSPInitializePackage - This function initializes everything in the TSSSP package. It is called
501:   from PetscDLLibraryRegister() when using dynamic libraries, and on the first call to TSCreate_SSP()
502:   when using static libraries.

504:   Level: developer

506: .keywords: TS, TSSSP, initialize, package
507: .seealso: PetscInitialize()
508: @*/
509: PetscErrorCode TSSSPInitializePackage(void)
510: {

514:   if (TSSSPPackageInitialized) return(0);
515:   TSSSPPackageInitialized = PETSC_TRUE;
516:   PetscFunctionListAdd(&TSSSPList,TSSSPRKS2, TSSSPStep_RK_2);
517:   PetscFunctionListAdd(&TSSSPList,TSSSPRKS3, TSSSPStep_RK_3);
518:   PetscFunctionListAdd(&TSSSPList,TSSSPRK104,TSSSPStep_RK_10_4);
519:   PetscRegisterFinalize(TSSSPFinalizePackage);
520:   return(0);
521: }

523: /*@C
524:   TSSSPFinalizePackage - This function destroys everything in the TSSSP package. It is
525:   called from PetscFinalize().

527:   Level: developer

529: .keywords: Petsc, destroy, package
530: .seealso: PetscFinalize()
531: @*/
532: PetscErrorCode TSSSPFinalizePackage(void)
533: {

537:   TSSSPPackageInitialized = PETSC_FALSE;
538:   PetscFunctionListDestroy(&TSSSPList);
539:   return(0);
540: }