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;

 23:   if (ssp->nwork < n) {
 24:     if (ssp->nwork > 0) VecDestroyVecs(ssp->nwork, &ssp->work);
 25:     VecDuplicateVecs(ts->vec_sol, n, &ssp->work);
 26:     ssp->nwork = n;
 27:   }
 28:   *work        = ssp->work;
 29:   ssp->workout = PETSC_TRUE;
 30:   return 0;
 31: }

 33: static PetscErrorCode TSSSPRestoreWorkVectors(TS ts, PetscInt n, Vec **work)
 34: {
 35:   TS_SSP *ssp = (TS_SSP *)ts->data;

 39:   ssp->workout = PETSC_FALSE;
 40:   *work        = NULL;
 41:   return 0;
 42: }

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

 47:    Pseudocode 2 of Ketcheson 2008

 49:    Level: beginner

 51: .seealso: `TSSSP`, `TSSSPSetType()`, `TSSSPSetNumStages()`
 52: M*/
 53: static PetscErrorCode TSSSPStep_RK_2(TS ts, PetscReal t0, PetscReal dt, Vec sol)
 54: {
 55:   TS_SSP  *ssp = (TS_SSP *)ts->data;
 56:   Vec     *work, F;
 57:   PetscInt i, s;

 59:   s = ssp->nstages;
 60:   TSSSPGetWorkVectors(ts, 2, &work);
 61:   F = work[1];
 62:   VecCopy(sol, work[0]);
 63:   for (i = 0; i < s - 1; i++) {
 64:     PetscReal stage_time = t0 + dt * (i / (s - 1.));
 65:     TSPreStage(ts, stage_time);
 66:     TSComputeRHSFunction(ts, stage_time, work[0], F);
 67:     VecAXPY(work[0], dt / (s - 1.), F);
 68:   }
 69:   TSComputeRHSFunction(ts, t0 + dt, work[0], F);
 70:   VecAXPBYPCZ(sol, (s - 1.) / s, dt / s, 1. / s, work[0], F);
 71:   TSSSPRestoreWorkVectors(ts, 2, &work);
 72:   return 0;
 73: }

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

 78:    Pseudocode 2 of Ketcheson 2008

 80:    Level: beginner

 82: .seealso: `TSSSP`, `TSSSPSetType()`, `TSSSPSetNumStages()`
 83: M*/
 84: static PetscErrorCode TSSSPStep_RK_3(TS ts, PetscReal t0, PetscReal dt, Vec sol)
 85: {
 86:   TS_SSP   *ssp = (TS_SSP *)ts->data;
 87:   Vec      *work, F;
 88:   PetscInt  i, s, n, r;
 89:   PetscReal c, stage_time;

 91:   s = ssp->nstages;
 92:   n = (PetscInt)(PetscSqrtReal((PetscReal)s) + 0.001);
 93:   r = s - n;
 95:   TSSSPGetWorkVectors(ts, 3, &work);
 96:   F = work[2];
 97:   VecCopy(sol, work[0]);
 98:   for (i = 0; i < (n - 1) * (n - 2) / 2; i++) {
 99:     c          = (i < n * (n + 1) / 2) ? 1. * i / (s - n) : (1. * i - n) / (s - n);
100:     stage_time = t0 + c * dt;
101:     TSPreStage(ts, stage_time);
102:     TSComputeRHSFunction(ts, stage_time, work[0], F);
103:     VecAXPY(work[0], dt / r, F);
104:   }
105:   VecCopy(work[0], work[1]);
106:   for (; i < n * (n + 1) / 2 - 1; i++) {
107:     c          = (i < n * (n + 1) / 2) ? 1. * i / (s - n) : (1. * i - n) / (s - n);
108:     stage_time = t0 + c * dt;
109:     TSPreStage(ts, stage_time);
110:     TSComputeRHSFunction(ts, stage_time, work[0], F);
111:     VecAXPY(work[0], dt / r, F);
112:   }
113:   {
114:     c          = (i < n * (n + 1) / 2) ? 1. * i / (s - n) : (1. * i - n) / (s - n);
115:     stage_time = t0 + c * dt;
116:     TSPreStage(ts, stage_time);
117:     TSComputeRHSFunction(ts, stage_time, work[0], F);
118:     VecAXPBYPCZ(work[0], 1. * n / (2 * n - 1.), (n - 1.) * dt / (r * (2 * n - 1)), (n - 1.) / (2 * n - 1.), work[1], F);
119:     i++;
120:   }
121:   for (; i < s; i++) {
122:     c          = (i < n * (n + 1) / 2) ? 1. * i / (s - n) : (1. * i - n) / (s - n);
123:     stage_time = t0 + c * dt;
124:     TSPreStage(ts, stage_time);
125:     TSComputeRHSFunction(ts, stage_time, work[0], F);
126:     VecAXPY(work[0], dt / r, F);
127:   }
128:   VecCopy(work[0], sol);
129:   TSSSPRestoreWorkVectors(ts, 3, &work);
130:   return 0;
131: }

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

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

138:    Level: beginner

140: .seealso: `TSSSP`, `TSSSPSetType()`
141: M*/
142: static PetscErrorCode TSSSPStep_RK_10_4(TS ts, PetscReal t0, PetscReal dt, Vec sol)
143: {
144:   const PetscReal c[10] = {0, 1. / 6, 2. / 6, 3. / 6, 4. / 6, 2. / 6, 3. / 6, 4. / 6, 5. / 6, 1};
145:   Vec            *work, F;
146:   PetscInt        i;
147:   PetscReal       stage_time;

149:   TSSSPGetWorkVectors(ts, 3, &work);
150:   F = work[2];
151:   VecCopy(sol, work[0]);
152:   for (i = 0; i < 5; i++) {
153:     stage_time = t0 + c[i] * dt;
154:     TSPreStage(ts, stage_time);
155:     TSComputeRHSFunction(ts, stage_time, work[0], F);
156:     VecAXPY(work[0], dt / 6, F);
157:   }
158:   VecAXPBYPCZ(work[1], 1. / 25, 9. / 25, 0, sol, work[0]);
159:   VecAXPBY(work[0], 15, -5, work[1]);
160:   for (; i < 9; i++) {
161:     stage_time = t0 + c[i] * dt;
162:     TSPreStage(ts, stage_time);
163:     TSComputeRHSFunction(ts, stage_time, work[0], F);
164:     VecAXPY(work[0], dt / 6, F);
165:   }
166:   stage_time = t0 + dt;
167:   TSPreStage(ts, stage_time);
168:   TSComputeRHSFunction(ts, stage_time, work[0], F);
169:   VecAXPBYPCZ(work[1], 3. / 5, dt / 10, 1, work[0], F);
170:   VecCopy(work[1], sol);
171:   TSSSPRestoreWorkVectors(ts, 3, &work);
172:   return 0;
173: }

175: static PetscErrorCode TSSetUp_SSP(TS ts)
176: {
177:   TSCheckImplicitTerm(ts);
178:   TSGetAdapt(ts, &ts->adapt);
179:   TSAdaptCandidatesClear(ts->adapt);
180:   return 0;
181: }

183: static PetscErrorCode TSStep_SSP(TS ts)
184: {
185:   TS_SSP   *ssp = (TS_SSP *)ts->data;
186:   Vec       sol = ts->vec_sol;
187:   PetscBool stageok, accept = PETSC_TRUE;
188:   PetscReal next_time_step = ts->time_step;

190:   (*ssp->onestep)(ts, ts->ptime, ts->time_step, sol);
191:   TSPostStage(ts, ts->ptime, 0, &sol);
192:   TSAdaptCheckStage(ts->adapt, ts, ts->ptime + ts->time_step, sol, &stageok);
193:   if (!stageok) {
194:     ts->reason = TS_DIVERGED_STEP_REJECTED;
195:     return 0;
196:   }

198:   TSAdaptChoose(ts->adapt, ts, ts->time_step, NULL, &next_time_step, &accept);
199:   if (!accept) {
200:     ts->reason = TS_DIVERGED_STEP_REJECTED;
201:     return 0;
202:   }

204:   ts->ptime += ts->time_step;
205:   ts->time_step = next_time_step;
206:   return 0;
207: }
208: /*------------------------------------------------------------*/

210: static PetscErrorCode TSReset_SSP(TS ts)
211: {
212:   TS_SSP *ssp = (TS_SSP *)ts->data;

214:   if (ssp->work) VecDestroyVecs(ssp->nwork, &ssp->work);
215:   ssp->nwork   = 0;
216:   ssp->workout = PETSC_FALSE;
217:   return 0;
218: }

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

224:   TSReset_SSP(ts);
225:   PetscFree(ssp->type_name);
226:   PetscFree(ts->data);
227:   PetscObjectComposeFunction((PetscObject)ts, "TSSSPGetType_C", NULL);
228:   PetscObjectComposeFunction((PetscObject)ts, "TSSSPSetType_C", NULL);
229:   PetscObjectComposeFunction((PetscObject)ts, "TSSSPGetNumStages_C", NULL);
230:   PetscObjectComposeFunction((PetscObject)ts, "TSSSPSetNumStages_C", NULL);
231:   return 0;
232: }
233: /*------------------------------------------------------------*/

235: /*@C
236:    TSSSPSetType - set the SSP time integration scheme to use

238:    Logically Collective

240:    Input Parameters:
241: +  ts - time stepping object
242: -  ssptype - type of scheme to use

244:    Options Database Keys:
245:    -ts_ssp_type <rks2>               : Type of SSP method (one of) rks2 rks3 rk104
246:    -ts_ssp_nstages<rks2: 5, rks3: 9> : Number of stages

248:    Level: beginner

250: .seealso: `TSSSP`, `TSSSPGetType()`, `TSSSPSetNumStages()`, `TSSSPRKS2`, `TSSSPRKS3`, `TSSSPRK104`
251: @*/
252: PetscErrorCode TSSSPSetType(TS ts, TSSSPType ssptype)
253: {
256:   PetscTryMethod(ts, "TSSSPSetType_C", (TS, TSSSPType), (ts, ssptype));
257:   return 0;
258: }

260: /*@C
261:    TSSSPGetType - get the SSP time integration scheme

263:    Logically Collective

265:    Input Parameter:
266: .  ts - time stepping object

268:    Output Parameter:
269: .  type - type of scheme being used

271:    Level: beginner

273: .seealso: `TSSSP`, `TSSSPSettype()`, `TSSSPSetNumStages()`, `TSSSPRKS2`, `TSSSPRKS3`, `TSSSPRK104`
274: @*/
275: PetscErrorCode TSSSPGetType(TS ts, TSSSPType *type)
276: {
278:   PetscUseMethod(ts, "TSSSPGetType_C", (TS, TSSSPType *), (ts, type));
279:   return 0;
280: }

282: /*@
283:    TSSSPSetNumStages - set the number of stages to use with the SSP method. Must be called after
284:    TSSSPSetType().

286:    Logically Collective

288:    Input Parameters:
289: +  ts - time stepping object
290: -  nstages - number of stages

292:    Options Database Keys:
293:    -ts_ssp_type <rks2>               : Type of SSP method (one of) rks2 rks3 rk104
294:    -ts_ssp_nstages<rks2: 5, rks3: 9> : Number of stages

296:    Level: beginner

298: .seealso: `TSSSP`, `TSSSPGetNumStages()`, `TSSSPSetNumStages()`, `TSSSPRKS2`, `TSSSPRKS3`, `TSSSPRK104`
299: @*/
300: PetscErrorCode TSSSPSetNumStages(TS ts, PetscInt nstages)
301: {
303:   PetscTryMethod(ts, "TSSSPSetNumStages_C", (TS, PetscInt), (ts, nstages));
304:   return 0;
305: }

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

310:    Logically Collective

312:    Input Parameter:
313: .  ts - time stepping object

315:    Output Parameter:
316: .  nstages - number of stages

318:    Level: beginner

320: .seealso: `TSSSP`, `TSSSPGetType()`, `TSSSPSetNumStages()`, `TSSSPRKS2`, `TSSSPRKS3`, `TSSSPRK104`
321: @*/
322: PetscErrorCode TSSSPGetNumStages(TS ts, PetscInt *nstages)
323: {
325:   PetscUseMethod(ts, "TSSSPGetNumStages_C", (TS, PetscInt *), (ts, nstages));
326:   return 0;
327: }

329: static PetscErrorCode TSSSPSetType_SSP(TS ts, TSSSPType type)
330: {
331:   TS_SSP *ssp = (TS_SSP *)ts->data;
332:   PetscErrorCode (*r)(TS, PetscReal, PetscReal, Vec);
333:   PetscBool flag;

335:   PetscFunctionListFind(TSSSPList, type, &r);
337:   ssp->onestep = r;
338:   PetscFree(ssp->type_name);
339:   PetscStrallocpy(type, &ssp->type_name);
340:   ts->default_adapt_type = TSADAPTNONE;
341:   PetscStrcmp(type, TSSSPRKS2, &flag);
342:   if (flag) {
343:     ssp->nstages = 5;
344:   } else {
345:     PetscStrcmp(type, TSSSPRKS3, &flag);
346:     if (flag) {
347:       ssp->nstages = 9;
348:     } else {
349:       ssp->nstages = 5;
350:     }
351:   }
352:   return 0;
353: }
354: static PetscErrorCode TSSSPGetType_SSP(TS ts, TSSSPType *type)
355: {
356:   TS_SSP *ssp = (TS_SSP *)ts->data;

358:   *type = ssp->type_name;
359:   return 0;
360: }
361: static PetscErrorCode TSSSPSetNumStages_SSP(TS ts, PetscInt nstages)
362: {
363:   TS_SSP *ssp = (TS_SSP *)ts->data;

365:   ssp->nstages = nstages;
366:   return 0;
367: }
368: static PetscErrorCode TSSSPGetNumStages_SSP(TS ts, PetscInt *nstages)
369: {
370:   TS_SSP *ssp = (TS_SSP *)ts->data;

372:   *nstages = ssp->nstages;
373:   return 0;
374: }

376: static PetscErrorCode TSSetFromOptions_SSP(TS ts, PetscOptionItems *PetscOptionsObject)
377: {
378:   char      tname[256] = TSSSPRKS2;
379:   TS_SSP   *ssp        = (TS_SSP *)ts->data;
380:   PetscBool flg;

382:   PetscOptionsHeadBegin(PetscOptionsObject, "SSP ODE solver options");
383:   {
384:     PetscOptionsFList("-ts_ssp_type", "Type of SSP method", "TSSSPSetType", TSSSPList, tname, tname, sizeof(tname), &flg);
385:     if (flg) TSSSPSetType(ts, tname);
386:     PetscOptionsInt("-ts_ssp_nstages", "Number of stages", "TSSSPSetNumStages", ssp->nstages, &ssp->nstages, NULL);
387:   }
388:   PetscOptionsHeadEnd();
389:   return 0;
390: }

392: static PetscErrorCode TSView_SSP(TS ts, PetscViewer viewer)
393: {
394:   TS_SSP   *ssp = (TS_SSP *)ts->data;
395:   PetscBool ascii;

397:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &ascii);
398:   if (ascii) PetscViewerASCIIPrintf(viewer, "  Scheme: %s\n", ssp->type_name);
399:   return 0;
400: }

402: /* ------------------------------------------------------------ */

404: /*MC
405:       TSSSP - Explicit strong stability preserving ODE solver

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

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

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

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

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

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

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

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

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

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

437:   Level: beginner

439:   References:
440: +  * - Ketcheson, Highly efficient strong stability preserving Runge Kutta methods with low storage implementations, SISC, 2008.
441: -  * - Gottlieb, Ketcheson, and Shu, High order strong stability preserving time discretizations, J Scientific Computing, 2009.

443: .seealso: `TSCreate()`, `TS`, `TSSetType()`

445: M*/
446: PETSC_EXTERN PetscErrorCode TSCreate_SSP(TS ts)
447: {
448:   TS_SSP *ssp;

450:   TSSSPInitializePackage();

452:   ts->ops->setup          = TSSetUp_SSP;
453:   ts->ops->step           = TSStep_SSP;
454:   ts->ops->reset          = TSReset_SSP;
455:   ts->ops->destroy        = TSDestroy_SSP;
456:   ts->ops->setfromoptions = TSSetFromOptions_SSP;
457:   ts->ops->view           = TSView_SSP;

459:   PetscNew(&ssp);
460:   ts->data = (void *)ssp;

462:   PetscObjectComposeFunction((PetscObject)ts, "TSSSPGetType_C", TSSSPGetType_SSP);
463:   PetscObjectComposeFunction((PetscObject)ts, "TSSSPSetType_C", TSSSPSetType_SSP);
464:   PetscObjectComposeFunction((PetscObject)ts, "TSSSPGetNumStages_C", TSSSPGetNumStages_SSP);
465:   PetscObjectComposeFunction((PetscObject)ts, "TSSSPSetNumStages_C", TSSSPSetNumStages_SSP);

467:   TSSSPSetType(ts, TSSSPRKS2);
468:   return 0;
469: }

471: /*@C
472:   TSSSPInitializePackage - This function initializes everything in the TSSSP package. It is called
473:   from TSInitializePackage().

475:   Level: developer

477: .seealso: `PetscInitialize()`
478: @*/
479: PetscErrorCode TSSSPInitializePackage(void)
480: {
481:   if (TSSSPPackageInitialized) return 0;
482:   TSSSPPackageInitialized = PETSC_TRUE;
483:   PetscFunctionListAdd(&TSSSPList, TSSSPRKS2, TSSSPStep_RK_2);
484:   PetscFunctionListAdd(&TSSSPList, TSSSPRKS3, TSSSPStep_RK_3);
485:   PetscFunctionListAdd(&TSSSPList, TSSSPRK104, TSSSPStep_RK_10_4);
486:   PetscRegisterFinalize(TSSSPFinalizePackage);
487:   return 0;
488: }

490: /*@C
491:   TSSSPFinalizePackage - This function destroys everything in the TSSSP package. It is
492:   called from PetscFinalize().

494:   Level: developer

496: .seealso: `PetscFinalize()`
497: @*/
498: PetscErrorCode TSSSPFinalizePackage(void)
499: {
500:   TSSSPPackageInitialized = PETSC_FALSE;
501:   PetscFunctionListDestroy(&TSSSPList);
502:   return 0;
503: }