Actual source code: rk.c

petsc-3.5.1 2014-07-24
Report Typos and Errors
  1: /*
  2:   Code for time stepping with the Runge-Kutta method

  4:   Notes:
  5:   The general system is written as

  7:   Udot = F(t,U)

  9: */
 10: #include <petsc-private/tsimpl.h>                /*I   "petscts.h"   I*/
 11: #include <petscdm.h>

 13: static TSRKType  TSRKDefault = TSRK3BS;
 14: static PetscBool TSRKRegisterAllCalled;
 15: static PetscBool TSRKPackageInitialized;
 16: static PetscInt  explicit_stage_time_id;

 18: typedef struct _RKTableau *RKTableau;
 19: struct _RKTableau {
 20:   char      *name;
 21:   PetscInt   order;               /* Classical approximation order of the method i              */
 22:   PetscInt   s;                   /* Number of stages                                           */
 23:   PetscBool  FSAL;                /* flag to indicate if tableau is FSAL                        */
 24:   PetscInt   pinterp;             /* Interpolation order                                        */
 25:   PetscReal *A,*b,*c;             /* Tableau                                                    */
 26:   PetscReal *bembed;              /* Embedded formula of order one less (order-1)               */
 27:   PetscReal *binterp;             /* Dense output formula                                       */
 28:   PetscReal  ccfl;                /* Placeholder for CFL coefficient relative to forward Euler  */
 29: };
 30: typedef struct _RKTableauLink *RKTableauLink;
 31: struct _RKTableauLink {
 32:   struct _RKTableau tab;
 33:   RKTableauLink     next;
 34: };
 35: static RKTableauLink RKTableauList;

 37: typedef struct {
 38:   RKTableau   tableau;
 39:   Vec          *Y;               /* States computed during the step */
 40:   Vec          *YdotRHS;         /* Function evaluations for the non-stiff part */
 41:   PetscScalar  *work;            /* Scalar work */
 42:   PetscReal    stage_time;
 43:   TSStepStatus status;
 44: } TS_RK;

 46: /*MC
 47:      TSRK1 - First order forward Euler scheme.

 49:      This method has one stage.

 51:      Level: advanced

 53: .seealso: TSRK
 54: M*/
 55: /*MC
 56:      TSRK2A - Second order RK scheme.

 58:      This method has two stages.

 60:      Level: advanced

 62: .seealso: TSRK
 63: M*/
 64: /*MC
 65:      TSRK3 - Third order RK scheme.

 67:      This method has three stages.

 69:      Level: advanced

 71: .seealso: TSRK
 72: M*/
 73: /*MC
 74:      TSRK3BS - Third order RK scheme of Bogacki-Shampine with 2nd order embedded method.

 76:      This method has four stages.

 78:      Level: advanced

 80: .seealso: TSRK
 81: M*/
 82: /*MC
 83:      TSRK4 - Fourth order RK scheme.

 85:      This is the classical Runge-Kutta method with four stages.

 87:      Level: advanced

 89: .seealso: TSRK
 90: M*/
 91: /*MC
 92:      TSRK5F - Fifth order Fehlberg RK scheme with a 4th order embedded method.

 94:      This method has six stages.

 96:      Level: advanced

 98: .seealso: TSRK
 99: M*/
100: /*MC
101:      TSRK5DP - Fifth order Dormand-Prince RK scheme with the 4th order embedded method.

103:      This method has seven stages.

105:      Level: advanced

107: .seealso: TSRK
108: M*/

112: /*@C
113:   TSRKRegisterAll - Registers all of the Runge-Kutta explicit methods in TSRK

115:   Not Collective, but should be called by all processes which will need the schemes to be registered

117:   Level: advanced

119: .keywords: TS, TSRK, register, all

121: .seealso:  TSRKRegisterDestroy()
122: @*/
123: PetscErrorCode TSRKRegisterAll(void)
124: {

128:   if (TSRKRegisterAllCalled) return(0);
129:   TSRKRegisterAllCalled = PETSC_TRUE;

131:   {
132:     const PetscReal
133:       A[1][1] = {{0.0}},
134:       b[1]    = {1.0};
135:     TSRKRegister(TSRK1FE,1,1,&A[0][0],b,NULL,NULL,1,b);
136:   }
137:   {
138:     const PetscReal
139:       A[2][2]     = {{0.0,0.0},
140:                     {1.0,0.0}},
141:       b[2]        = {0.5,0.5},
142:       bembed[2]   = {1.0,0};
143:     TSRKRegister(TSRK2A,2,2,&A[0][0],b,NULL,bembed,2,b);
144:   }
145:   {
146:     const PetscReal
147:       A[3][3] = {{0,0,0},
148:                  {2.0/3.0,0,0},
149:                  {-1.0/3.0,1.0,0}},
150:       b[3]    = {0.25,0.5,0.25};
151:     TSRKRegister(TSRK3,3,3,&A[0][0],b,NULL,NULL,3,b);
152:   }
153:   {
154:     const PetscReal
155:       A[4][4] = {{0,0,0,0},
156:                  {1.0/2.0,0,0,0},
157:                  {0,3.0/4.0,0,0},
158:                  {2.0/9.0,1.0/3.0,4.0/9.0,0}},
159:       b[4]    = {2.0/9.0,1.0/3.0,4.0/9.0,0},
160:       bembed[4] = {7.0/24.0,1.0/4.0,1.0/3.0,1.0/8.0};
161:     TSRKRegister(TSRK3BS,3,4,&A[0][0],b,NULL,bembed,3,b);
162:   }
163:   {
164:     const PetscReal
165:       A[4][4] = {{0,0,0,0},
166:                  {0.5,0,0,0},
167:                  {0,0.5,0,0},
168:                  {0,0,1.0,0}},
169:       b[4]    = {1.0/6.0,1.0/3.0,1.0/3.0,1.0/6.0};
170:     TSRKRegister(TSRK4,4,4,&A[0][0],b,NULL,NULL,4,b);
171:   }
172:   {
173:     const PetscReal
174:       A[6][6]   = {{0,0,0,0,0,0},
175:                    {0.25,0,0,0,0,0},
176:                    {3.0/32.0,9.0/32.0,0,0,0,0},
177:                    {1932.0/2197.0,-7200.0/2197.0,7296.0/2197.0,0,0,0},
178:                    {439.0/216.0,-8.0,3680.0/513.0,-845.0/4104.0,0,0},
179:                    {-8.0/27.0,2.0,-3544.0/2565.0,1859.0/4104.0,-11.0/40.0,0}},
180:       b[6]      = {16.0/135.0,0,6656.0/12825.0,28561.0/56430.0,-9.0/50.0,2.0/55.0},
181:       bembed[6] = {25.0/216.0,0,1408.0/2565.0,2197.0/4104.0,-1.0/5.0,0};
182:     TSRKRegister(TSRK5F,5,6,&A[0][0],b,NULL,bembed,5,b);
183:   }
184:   {
185:     const PetscReal
186:       A[7][7]   = {{0,0,0,0,0,0,0},
187:                    {1.0/5.0,0,0,0,0,0,0},
188:                    {3.0/40.0,9.0/40.0,0,0,0,0,0},
189:                    {44.0/45.0,-56.0/15.0,32.0/9.0,0,0,0,0},
190:                    {19372.0/6561.0,-25360.0/2187.0,64448.0/6561.0,-212.0/729.0,0,0,0},
191:                    {9017.0/3168.0,-355.0/33.0,46732.0/5247.0,49.0/176.0,-5103.0/18656.0,0,0},
192:                    {35.0/384.0,0,500.0/1113.0,125.0/192.0,-2187.0/6784.0,11.0/84.0,0}},
193:       b[7]      = {35.0/384.0,0,500.0/1113.0,125.0/192.0,-2187.0/6784.0,11.0/84.0,0},
194:       bembed[7] = {5179.0/57600.0,0,7571.0/16695.0,393.0/640.0,-92097.0/339200.0,187.0/2100.0,1.0/40.0};
195:     TSRKRegister(TSRK5DP,5,7,&A[0][0],b,NULL,bembed,5,b);
196:   }
197:   return(0);
198: }

202: /*@C
203:    TSRKRegisterDestroy - Frees the list of schemes that were registered by TSRKRegister().

205:    Not Collective

207:    Level: advanced

209: .keywords: TSRK, register, destroy
210: .seealso: TSRKRegister(), TSRKRegisterAll()
211: @*/
212: PetscErrorCode TSRKRegisterDestroy(void)
213: {
215:   RKTableauLink link;

218:   while ((link = RKTableauList)) {
219:     RKTableau t = &link->tab;
220:     RKTableauList = link->next;
221:     PetscFree3(t->A,t->b,t->c);
222:     PetscFree (t->bembed);
223:     PetscFree (t->binterp);
224:     PetscFree (t->name);
225:     PetscFree (link);
226:   }
227:   TSRKRegisterAllCalled = PETSC_FALSE;
228:   return(0);
229: }

233: /*@C
234:   TSRKInitializePackage - This function initializes everything in the TSRK package. It is called
235:   from PetscDLLibraryRegister() when using dynamic libraries, and on the first call to TSCreate_RK()
236:   when using static libraries.

238:   Level: developer

240: .keywords: TS, TSRK, initialize, package
241: .seealso: PetscInitialize()
242: @*/
243: PetscErrorCode TSRKInitializePackage(void)
244: {

248:   if (TSRKPackageInitialized) return(0);
249:   TSRKPackageInitialized = PETSC_TRUE;
250:   TSRKRegisterAll();
251:   PetscObjectComposedDataRegister(&explicit_stage_time_id);
252:   PetscRegisterFinalize(TSRKFinalizePackage);
253:   return(0);
254: }

258: /*@C
259:   TSRKFinalizePackage - This function destroys everything in the TSRK package. It is
260:   called from PetscFinalize().

262:   Level: developer

264: .keywords: Petsc, destroy, package
265: .seealso: PetscFinalize()
266: @*/
267: PetscErrorCode TSRKFinalizePackage(void)
268: {

272:   TSRKPackageInitialized = PETSC_FALSE;
273:   TSRKRegisterDestroy();
274:   return(0);
275: }

279: /*@C
280:    TSRKRegister - register an RK scheme by providing the entries in the Butcher tableau and optionally embedded approximations and interpolation

282:    Not Collective, but the same schemes should be registered on all processes on which they will be used

284:    Input Parameters:
285: +  name - identifier for method
286: .  order - approximation order of method
287: .  s - number of stages, this is the dimension of the matrices below
288: .  A - stage coefficients (dimension s*s, row-major)
289: .  b - step completion table (dimension s; NULL to use last row of A)
290: .  c - abscissa (dimension s; NULL to use row sums of A)
291: .  bembed - completion table for embedded method (dimension s; NULL if not available)
292: .  pinterp - Order of the interpolation scheme, equal to the number of columns of binterp
293: -  binterp - Coefficients of the interpolation formula (dimension s*pinterp; NULL to reuse binterpt)

295:    Notes:
296:    Several RK methods are provided, this function is only needed to create new methods.

298:    Level: advanced

300: .keywords: TS, register

302: .seealso: TSRK
303: @*/
304: PetscErrorCode TSRKRegister(TSRKType name,PetscInt order,PetscInt s,
305:                                  const PetscReal A[],const PetscReal b[],const PetscReal c[],
306:                                  const PetscReal bembed[],
307:                                  PetscInt pinterp,const PetscReal binterp[])
308: {
309:   PetscErrorCode  ierr;
310:   RKTableauLink   link;
311:   RKTableau       t;
312:   PetscInt        i,j;

315:   PetscMalloc(sizeof(*link),&link);
316:   PetscMemzero(link,sizeof(*link));
317:   t        = &link->tab;
318:   PetscStrallocpy(name,&t->name);
319:   t->order = order;
320:   t->s     = s;
321:   PetscMalloc3(s*s,&t->A,s,&t->b,s,&t->c);
322:   PetscMemcpy(t->A,A,s*s*sizeof(A[0]));
323:   if (b)  { PetscMemcpy(t->b,b,s*sizeof(b[0])); }
324:   else for (i=0; i<s; i++) t->b[i] = A[(s-1)*s+i];
325:   if (c)  { PetscMemcpy(t->c,c,s*sizeof(c[0])); }
326:   else for (i=0; i<s; i++) for (j=0,t->c[i]=0; j<s; j++) t->c[i] += A[i*s+j];
327:   t->FSAL = PETSC_TRUE;
328:   for (i=0; i<s; i++) if (t->A[(s-1)*s+i] != t->b[i]) t->FSAL = PETSC_FALSE;
329:   if (bembed) {
330:     PetscMalloc1(s,&t->bembed);
331:     PetscMemcpy(t->bembed,bembed,s*sizeof(bembed[0]));
332:   }

334:   t->pinterp     = pinterp;
335:   PetscMalloc1(s*pinterp,&t->binterp);
336:   PetscMemcpy(t->binterp,binterp,s*pinterp*sizeof(binterp[0]));
337:   link->next     = RKTableauList;
338:   RKTableauList = link;
339:   return(0);
340: }

344: /*
345:  The step completion formula is

347:  x1 = x0 + h b^T YdotRHS

349:  This function can be called before or after ts->vec_sol has been updated.
350:  Suppose we have a completion formula (b) and an embedded formula (be) of different order.
351:  We can write

353:  x1e = x0 + h be^T YdotRHS
354:      = x1 - h b^T YdotRHS + h be^T YdotRHS
355:      = x1 + h (be - b)^T YdotRHS

357:  so we can evaluate the method with different order even after the step has been optimistically completed.
358: */
359: static PetscErrorCode TSEvaluateStep_RK(TS ts,PetscInt order,Vec X,PetscBool *done)
360: {
361:   TS_RK         *rk   = (TS_RK*)ts->data;
362:   RKTableau      tab  = rk->tableau;
363:   PetscScalar   *w    = rk->work;
364:   PetscReal      h;
365:   PetscInt       s    = tab->s,j;

369:   switch (rk->status) {
370:   case TS_STEP_INCOMPLETE:
371:   case TS_STEP_PENDING:
372:     h = ts->time_step; break;
373:   case TS_STEP_COMPLETE:
374:     h = ts->time_step_prev; break;
375:   default: SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_PLIB,"Invalid TSStepStatus");
376:   }
377:   if (order == tab->order) {
378:     if (rk->status == TS_STEP_INCOMPLETE) {
379:       VecCopy(ts->vec_sol,X);
380:       for (j=0; j<s; j++) w[j] = h*tab->b[j];
381:       VecMAXPY(X,s,w,rk->YdotRHS);
382:     } else {VecCopy(ts->vec_sol,X);}
383:     if (done) *done = PETSC_TRUE;
384:     return(0);
385:   } else if (order == tab->order-1) {
386:     if (!tab->bembed) goto unavailable;
387:     if (rk->status == TS_STEP_INCOMPLETE) { /* Complete with the embedded method (be) */
388:       VecCopy(ts->vec_sol,X);
389:       for (j=0; j<s; j++) w[j] = h*tab->bembed[j];
390:       VecMAXPY(X,s,w,rk->YdotRHS);
391:     } else {                    /* Rollback and re-complete using (be-b) */
392:       VecCopy(ts->vec_sol,X);
393:       for (j=0; j<s; j++) w[j] = h*(tab->bembed[j] - tab->b[j]);
394:       VecMAXPY(X,s,w,rk->YdotRHS);
395:     }
396:     if (done) *done = PETSC_TRUE;
397:     return(0);
398:   }
399: unavailable:
400:   if (done) *done = PETSC_FALSE;
401:   else SETERRQ3(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"RK '%s' of order %D cannot evaluate step at order %D",tab->name,tab->order,order);
402:   return(0);
403: }

407: static PetscErrorCode TSStep_RK(TS ts)
408: {
409:   TS_RK           *rk   = (TS_RK*)ts->data;
410:   RKTableau        tab  = rk->tableau;
411:   const PetscInt   s    = tab->s;
412:   const PetscReal *A = tab->A,*b = tab->b,*c = tab->c;
413:   PetscScalar     *w    = rk->work;
414:   Vec             *Y    = rk->Y,*YdotRHS = rk->YdotRHS;
415:   TSAdapt          adapt;
416:   PetscInt         i,j,reject,next_scheme;
417:   PetscReal        next_time_step;
418:   PetscReal        t;
419:   PetscBool        accept;
420:   PetscErrorCode   ierr;


424:   next_time_step = ts->time_step;
425:   t              = ts->ptime;
426:   accept         = PETSC_TRUE;
427:   rk->status     = TS_STEP_INCOMPLETE;


430:   for (reject=0; reject<ts->max_reject && !ts->reason; reject++,ts->reject++) {
431:     PetscReal h = ts->time_step;
432:     TSPreStep(ts);
433:     for (i=0; i<s; i++) {
434:       rk->stage_time = t + h*c[i];
435:       TSPreStage(ts,rk->stage_time);
436:       VecCopy(ts->vec_sol,Y[i]);
437:       for (j=0; j<i; j++) w[j] = h*A[i*s+j];
438:       VecMAXPY(Y[i],i,w,YdotRHS);
439:       TSPostStage(ts,rk->stage_time,i,Y);
440:       TSGetAdapt(ts,&adapt);
441:       TSAdaptCheckStage(adapt,ts,&accept);
442:       if (!accept) goto reject_step;
443:       TSComputeRHSFunction(ts,t+h*c[i],Y[i],YdotRHS[i]);
444:     }
445:     TSEvaluateStep(ts,tab->order,ts->vec_sol,NULL);
446:     rk->status = TS_STEP_PENDING;

448:     /* Register only the current method as a candidate because we're not supporting multiple candidates yet. */
449:     TSGetAdapt(ts,&adapt);
450:     TSAdaptCandidatesClear(adapt);
451:     TSAdaptCandidateAdd(adapt,tab->name,tab->order,1,tab->ccfl,1.*tab->s,PETSC_TRUE);
452:     TSAdaptChoose(adapt,ts,ts->time_step,&next_scheme,&next_time_step,&accept);
453:     if (accept) {
454:       /* ignore next_scheme for now */
455:       ts->ptime    += ts->time_step;
456:       ts->time_step = next_time_step;
457:       ts->steps++;
458:       rk->status = TS_STEP_COMPLETE;
459:       PetscObjectComposedDataSetReal((PetscObject)ts->vec_sol,explicit_stage_time_id,ts->ptime);
460:       break;
461:     } else {                    /* Roll back the current step */
462:       for (j=0; j<s; j++) w[j] = -h*b[j];
463:       VecMAXPY(ts->vec_sol,s,w,rk->YdotRHS);
464:       ts->time_step = next_time_step;
465:       rk->status   = TS_STEP_INCOMPLETE;
466:     }
467: reject_step: continue;
468:   }
469:   if (rk->status != TS_STEP_COMPLETE && !ts->reason) ts->reason = TS_DIVERGED_STEP_REJECTED;
470:   return(0);
471: }

475: static PetscErrorCode TSInterpolate_RK(TS ts,PetscReal itime,Vec X)
476: {
477:   TS_RK           *rk = (TS_RK*)ts->data;
478:   PetscInt         s  = rk->tableau->s,pinterp = rk->tableau->pinterp,i,j;
479:   PetscReal        h;
480:   PetscReal        tt,t;
481:   PetscScalar     *b;
482:   const PetscReal *B = rk->tableau->binterp;
483:   PetscErrorCode   ierr;

486:   if (!B) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"TSRK %s does not have an interpolation formula",rk->tableau->name);
487:   switch (rk->status) {
488:   case TS_STEP_INCOMPLETE:
489:   case TS_STEP_PENDING:
490:     h = ts->time_step;
491:     t = (itime - ts->ptime)/h;
492:     break;
493:   case TS_STEP_COMPLETE:
494:     h = ts->time_step_prev;
495:     t = (itime - ts->ptime)/h + 1; /* In the interval [0,1] */
496:     break;
497:   default: SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_PLIB,"Invalid TSStepStatus");
498:   }
499:   PetscMalloc1(s,&b);
500:   for (i=0; i<s; i++) b[i] = 0;
501:   for (j=0,tt=t; j<pinterp; j++,tt*=t) {
502:     for (i=0; i<s; i++) {
503:       b[i]  += h * B[i*pinterp+j] * tt;
504:     }
505:   }
506:   VecCopy(rk->Y[0],X);
507:   VecMAXPY(X,s,b,rk->YdotRHS);
508:   PetscFree(b);
509:   return(0);
510: }

512: /*------------------------------------------------------------*/
515: static PetscErrorCode TSReset_RK(TS ts)
516: {
517:   TS_RK         *rk = (TS_RK*)ts->data;
518:   PetscInt       s;

522:   if (!rk->tableau) return(0);
523:   s    = rk->tableau->s;
524:   VecDestroyVecs(s,&rk->Y);
525:   VecDestroyVecs(s,&rk->YdotRHS);
526:   PetscFree(rk->work);
527:   return(0);
528: }

532: static PetscErrorCode TSDestroy_RK(TS ts)
533: {

537:   TSReset_RK(ts);
538:   PetscFree(ts->data);
539:   PetscObjectComposeFunction((PetscObject)ts,"TSRKGetType_C",NULL);
540:   PetscObjectComposeFunction((PetscObject)ts,"TSRKSetType_C",NULL);
541:   return(0);
542: }


547: static PetscErrorCode DMCoarsenHook_TSRK(DM fine,DM coarse,void *ctx)
548: {
550:   return(0);
551: }

555: static PetscErrorCode DMRestrictHook_TSRK(DM fine,Mat restrct,Vec rscale,Mat inject,DM coarse,void *ctx)
556: {
558:   return(0);
559: }


564: static PetscErrorCode DMSubDomainHook_TSRK(DM dm,DM subdm,void *ctx)
565: {
567:   return(0);
568: }

572: static PetscErrorCode DMSubDomainRestrictHook_TSRK(DM dm,VecScatter gscat,VecScatter lscat,DM subdm,void *ctx)
573: {

576:   return(0);
577: }

581: static PetscErrorCode TSSetUp_RK(TS ts)
582: {
583:   TS_RK         *rk = (TS_RK*)ts->data;
584:   RKTableau      tab;
585:   PetscInt       s;
587:   DM             dm;

590:   if (!rk->tableau) {
591:     TSRKSetType(ts,TSRKDefault);
592:   }
593:   tab  = rk->tableau;
594:   s    = tab->s;
595:   VecDuplicateVecs(ts->vec_sol,s,&rk->Y);
596:   VecDuplicateVecs(ts->vec_sol,s,&rk->YdotRHS);
597:   PetscMalloc1(s,&rk->work);
598:   TSGetDM(ts,&dm);
599:   if (dm) {
600:     DMCoarsenHookAdd(dm,DMCoarsenHook_TSRK,DMRestrictHook_TSRK,ts);
601:     DMSubDomainHookAdd(dm,DMSubDomainHook_TSRK,DMSubDomainRestrictHook_TSRK,ts);
602:   }
603:   return(0);
604: }
605: /*------------------------------------------------------------*/

609: static PetscErrorCode TSSetFromOptions_RK(TS ts)
610: {
612:   char           rktype[256];

615:   PetscOptionsHead("RK ODE solver options");
616:   {
617:     RKTableauLink  link;
618:     PetscInt       count,choice;
619:     PetscBool      flg;
620:     const char   **namelist;
621:     PetscStrncpy(rktype,TSRKDefault,sizeof(rktype));
622:     for (link=RKTableauList,count=0; link; link=link->next,count++) ;
623:     PetscMalloc1(count,&namelist);
624:     for (link=RKTableauList,count=0; link; link=link->next,count++) namelist[count] = link->tab.name;
625:     PetscOptionsEList("-ts_rk_type","Family of RK method","TSRKSetType",(const char*const*)namelist,count,rktype,&choice,&flg);
626:     TSRKSetType(ts,flg ? namelist[choice] : rktype);
627:     PetscFree(namelist);
628:   }
629:   PetscOptionsTail();
630:   return(0);
631: }

635: static PetscErrorCode PetscFormatRealArray(char buf[],size_t len,const char *fmt,PetscInt n,const PetscReal x[])
636: {
638:   PetscInt       i;
639:   size_t         left,count;
640:   char           *p;

643:   for (i=0,p=buf,left=len; i<n; i++) {
644:     PetscSNPrintfCount(p,left,fmt,&count,x[i]);
645:     if (count >= left) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Insufficient space in buffer");
646:     left -= count;
647:     p    += count;
648:     *p++  = ' ';
649:   }
650:   p[i ? 0 : -1] = 0;
651:   return(0);
652: }

656: static PetscErrorCode TSView_RK(TS ts,PetscViewer viewer)
657: {
658:   TS_RK         *rk   = (TS_RK*)ts->data;
659:   RKTableau      tab  = rk->tableau;
660:   PetscBool      iascii;
662:   TSAdapt        adapt;

665:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
666:   if (iascii) {
667:     TSRKType rktype;
668:     char     buf[512];
669:     TSRKGetType(ts,&rktype);
670:     PetscViewerASCIIPrintf(viewer,"  RK %s\n",rktype);
671:     PetscFormatRealArray(buf,sizeof(buf),"% 8.6f",tab->s,tab->c);
672:     PetscViewerASCIIPrintf(viewer,"  Abscissa     c = %s\n",buf);
673:     PetscViewerASCIIPrintf(viewer,"FSAL: %s\n",tab->FSAL ? "yes" : "no");
674:   }
675:   TSGetAdapt(ts,&adapt);
676:   TSAdaptView(adapt,viewer);
677:   return(0);
678: }

682: static PetscErrorCode TSLoad_RK(TS ts,PetscViewer viewer)
683: {
685:   TSAdapt        tsadapt;

688:   TSGetAdapt(ts,&tsadapt);
689:   TSAdaptLoad(tsadapt,viewer);
690:   return(0);
691: }

695: /*@C
696:   TSRKSetType - Set the type of RK scheme

698:   Logically collective

700:   Input Parameter:
701: +  ts - timestepping context
702: -  rktype - type of RK-scheme

704:   Level: intermediate

706: .seealso: TSRKGetType(), TSRK, TSRK2, TSRK3, TSRKPRSSP2, TSRK4, TSRK5
707: @*/
708: PetscErrorCode TSRKSetType(TS ts,TSRKType rktype)
709: {

714:   PetscTryMethod(ts,"TSRKSetType_C",(TS,TSRKType),(ts,rktype));
715:   return(0);
716: }

720: /*@C
721:   TSRKGetType - Get the type of RK scheme

723:   Logically collective

725:   Input Parameter:
726: .  ts - timestepping context

728:   Output Parameter:
729: .  rktype - type of RK-scheme

731:   Level: intermediate

733: .seealso: TSRKGetType()
734: @*/
735: PetscErrorCode TSRKGetType(TS ts,TSRKType *rktype)
736: {

741:   PetscUseMethod(ts,"TSRKGetType_C",(TS,TSRKType*),(ts,rktype));
742:   return(0);
743: }

747: PetscErrorCode  TSRKGetType_RK(TS ts,TSRKType *rktype)
748: {
749:   TS_RK     *rk = (TS_RK*)ts->data;

753:   if (!rk->tableau) {
754:     TSRKSetType(ts,TSRKDefault);
755:   }
756:   *rktype = rk->tableau->name;
757:   return(0);
758: }
761: PetscErrorCode  TSRKSetType_RK(TS ts,TSRKType rktype)
762: {
763:   TS_RK     *rk = (TS_RK*)ts->data;
765:   PetscBool      match;
766:   RKTableauLink link;

769:   if (rk->tableau) {
770:     PetscStrcmp(rk->tableau->name,rktype,&match);
771:     if (match) return(0);
772:   }
773:   for (link = RKTableauList; link; link=link->next) {
774:     PetscStrcmp(link->tab.name,rktype,&match);
775:     if (match) {
776:       TSReset_RK(ts);
777:       rk->tableau = &link->tab;
778:       return(0);
779:     }
780:   }
781:   SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_UNKNOWN_TYPE,"Could not find '%s'",rktype);
782:   return(0);
783: }

785: /* ------------------------------------------------------------ */
786: /*MC
787:       TSRK - ODE and DAE solver using Runge-Kutta schemes

789:   The user should provide the right hand side of the equation 
790:   using TSSetRHSFunction().

792:   Notes:
793:   The default is TSRK3, it can be changed with TSRKSetType() or -ts_rk_type

795:   Level: beginner

797: .seealso:  TSCreate(), TS, TSSetType(), TSRKSetType(), TSRKGetType(), TSRKSetFullyImplicit(), TSRK2D, TTSRK2E, TSRK3,
798:            TSRK4, TSRK5, TSRKPRSSP2, TSRKBPR3, TSRKType, TSRKRegister()

800: M*/
803: PETSC_EXTERN PetscErrorCode TSCreate_RK(TS ts)
804: {
805:   TS_RK     *th;

809: #if !defined(PETSC_USE_DYNAMIC_LIBRARIES)
810:   TSRKInitializePackage();
811: #endif

813:   ts->ops->reset          = TSReset_RK;
814:   ts->ops->destroy        = TSDestroy_RK;
815:   ts->ops->view           = TSView_RK;
816:   ts->ops->load           = TSLoad_RK;
817:   ts->ops->setup          = TSSetUp_RK;
818:   ts->ops->step           = TSStep_RK;
819:   ts->ops->interpolate    = TSInterpolate_RK;
820:   ts->ops->evaluatestep   = TSEvaluateStep_RK;
821:   ts->ops->setfromoptions = TSSetFromOptions_RK;

823:   PetscNewLog(ts,&th);
824:   ts->data = (void*)th;

826:   PetscObjectComposeFunction((PetscObject)ts,"TSRKGetType_C",TSRKGetType_RK);
827:   PetscObjectComposeFunction((PetscObject)ts,"TSRKSetType_C",TSRKSetType_RK);
828:   return(0);
829: }