Actual source code: rk.c

petsc-3.6.1 2015-07-22
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:   Vec          *VecDeltaLam;     /* Increment of the adjoint sensitivity w.r.t IC at stage*/
 42:   Vec          *VecDeltaMu;      /* Increment of the adjoint sensitivity w.r.t P at stage*/
 43:   Vec          *VecSensiTemp;    /* Vector to be timed with Jacobian transpose*/
 44:   PetscScalar  *work;            /* Scalar work */
 45:   PetscReal    stage_time;
 46:   TSStepStatus status;
 47: } TS_RK;

 49: /*MC
 50:      TSRK1 - First order forward Euler scheme.

 52:      This method has one stage.

 54:      Level: advanced

 56: .seealso: TSRK
 57: M*/
 58: /*MC
 59:      TSRK2A - Second order RK scheme.

 61:      This method has two stages.

 63:      Level: advanced

 65: .seealso: TSRK
 66: M*/
 67: /*MC
 68:      TSRK3 - Third order RK scheme.

 70:      This method has three stages.

 72:      Level: advanced

 74: .seealso: TSRK
 75: M*/
 76: /*MC
 77:      TSRK3BS - Third order RK scheme of Bogacki-Shampine with 2nd order embedded method.

 79:      This method has four stages.

 81:      Level: advanced

 83: .seealso: TSRK
 84: M*/
 85: /*MC
 86:      TSRK4 - Fourth order RK scheme.

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

 90:      Level: advanced

 92: .seealso: TSRK
 93: M*/
 94: /*MC
 95:      TSRK5F - Fifth order Fehlberg RK scheme with a 4th order embedded method.

 97:      This method has six stages.

 99:      Level: advanced

101: .seealso: TSRK
102: M*/
103: /*MC
104:      TSRK5DP - Fifth order Dormand-Prince RK scheme with the 4th order embedded method.

106:      This method has seven stages.

108:      Level: advanced

110: .seealso: TSRK
111: M*/

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

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

120:   Level: advanced

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

124: .seealso:  TSRKRegisterDestroy()
125: @*/
126: PetscErrorCode TSRKRegisterAll(void)
127: {

131:   if (TSRKRegisterAllCalled) return(0);
132:   TSRKRegisterAllCalled = PETSC_TRUE;

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

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

208:    Not Collective

210:    Level: advanced

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

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

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

241:   Level: developer

243: .keywords: TS, TSRK, initialize, package
244: .seealso: PetscInitialize()
245: @*/
246: PetscErrorCode TSRKInitializePackage(void)
247: {

251:   if (TSRKPackageInitialized) return(0);
252:   TSRKPackageInitialized = PETSC_TRUE;
253:   TSRKRegisterAll();
254:   PetscObjectComposedDataRegister(&explicit_stage_time_id);
255:   PetscRegisterFinalize(TSRKFinalizePackage);
256:   return(0);
257: }

261: /*@C
262:   TSRKFinalizePackage - This function destroys everything in the TSRK package. It is
263:   called from PetscFinalize().

265:   Level: developer

267: .keywords: Petsc, destroy, package
268: .seealso: PetscFinalize()
269: @*/
270: PetscErrorCode TSRKFinalizePackage(void)
271: {

275:   TSRKPackageInitialized = PETSC_FALSE;
276:   TSRKRegisterDestroy();
277:   return(0);
278: }

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

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

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

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

301:    Level: advanced

303: .keywords: TS, register

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

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

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

347: /*
348:  The step completion formula is

350:  x1 = x0 + h b^T YdotRHS

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

356:  x1e = x0 + h be^T YdotRHS
357:      = x1 - h b^T YdotRHS + h be^T YdotRHS
358:      = x1 + h (be - b)^T YdotRHS

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

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

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


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


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

450:     /* Register only the current method as a candidate because we're not supporting multiple candidates yet. */
451:     TSGetAdapt(ts,&adapt);
452:     TSAdaptCandidatesClear(adapt);
453:     TSAdaptCandidateAdd(adapt,tab->name,tab->order,1,tab->ccfl,1.*tab->s,PETSC_TRUE);
454:     TSAdaptChoose(adapt,ts,ts->time_step,&next_scheme,&next_time_step,&accept);
455:     if (accept) {
456:       if (ts->costintegralfwd) {
457:         /* Evolve ts->vec_costintegral to compute integrals */
458:         for (i=0; i<s; i++) {
459:           TSAdjointComputeCostIntegrand(ts,t+h*c[i],Y[i],ts->vec_costintegrand);
460:           VecAXPY(ts->vec_costintegral,h*b[i],ts->vec_costintegrand);
461:         }
462:       }

464:       /* ignore next_scheme for now */
465:       ts->ptime    += ts->time_step;
466:       ts->time_step = next_time_step;
467:       ts->steps++;
468:       rk->status = TS_STEP_COMPLETE;
469:       PetscObjectComposedDataSetReal((PetscObject)ts->vec_sol,explicit_stage_time_id,ts->ptime);
470:       break;
471:     } else {                    /* Roll back the current step */
472:       for (j=0; j<s; j++) w[j] = -h*b[j];
473:       VecMAXPY(ts->vec_sol,s,w,rk->YdotRHS);
474:       ts->time_step = next_time_step;
475:       rk->status   = TS_STEP_INCOMPLETE;
476:     }
477: reject_step: continue;
478:   }
479:   if (rk->status != TS_STEP_COMPLETE && !ts->reason) ts->reason = TS_DIVERGED_STEP_REJECTED;
480:   return(0);
481: }

485: static PetscErrorCode TSAdjointSetUp_RK(TS ts)
486: {
487:   TS_RK         *rk = (TS_RK*)ts->data;
488:   RKTableau      tab;
489:   PetscInt       s;

493:   if (ts->adjointsetupcalled++) return(0);
494:   tab  = rk->tableau;
495:   s    = tab->s;
496:   VecDuplicateVecs(ts->vecs_sensi[0],s*ts->numcost,&rk->VecDeltaLam);
497:   if(ts->vecs_sensip) {
498:     VecDuplicateVecs(ts->vecs_sensip[0],s*ts->numcost,&rk->VecDeltaMu);
499:   }
500:   VecDuplicateVecs(ts->vecs_sensi[0],ts->numcost,&rk->VecSensiTemp);
501:   return(0);
502: }

506: static PetscErrorCode TSAdjointStep_RK(TS ts)
507: {
508:   TS_RK           *rk   = (TS_RK*)ts->data;
509:   RKTableau        tab  = rk->tableau;
510:   const PetscInt   s    = tab->s;
511:   const PetscReal *A = tab->A,*b = tab->b,*c = tab->c;
512:   PetscScalar     *w    = rk->work;
513:   Vec             *Y    = rk->Y,*VecDeltaLam = rk->VecDeltaLam,*VecDeltaMu = rk->VecDeltaMu,*VecSensiTemp = rk->VecSensiTemp;
514:   PetscInt         i,j,nadj;
515:   PetscReal        t;
516:   PetscErrorCode   ierr;
517:   PetscReal        h = ts->time_step;
518:   Mat              J,Jp;

521:   t          = ts->ptime;
522:   rk->status = TS_STEP_INCOMPLETE;
523:   h = ts->time_step;
524:   TSPreStep(ts);
525:   for (i=s-1; i>=0; i--) {
526:     rk->stage_time = t + h*(1.0-c[i]);
527:     for (nadj=0; nadj<ts->numcost; nadj++) {
528:       VecCopy(ts->vecs_sensi[nadj],VecSensiTemp[nadj]);
529:       VecScale(VecSensiTemp[nadj],-h*b[i]);
530:       for (j=i+1; j<s; j++) {
531:         VecAXPY(VecSensiTemp[nadj],-h*A[j*s+i],VecDeltaLam[nadj*s+j]);
532:       }
533:     }
534:     /* Stage values of lambda */
535:     TSGetRHSJacobian(ts,&J,&Jp,NULL,NULL);
536:     TSComputeRHSJacobian(ts,rk->stage_time,Y[i],J,Jp);
537:     if (ts->vec_costintegral) {
538:       TSAdjointComputeDRDYFunction(ts,rk->stage_time,Y[i],ts->vecs_drdy);
539:       if (!ts->costintegralfwd) {
540:         /* Evolve ts->vec_costintegral to compute integrals */
541:         TSAdjointComputeCostIntegrand(ts,rk->stage_time,Y[i],ts->vec_costintegrand);
542:         VecAXPY(ts->vec_costintegral,-h*b[i],ts->vec_costintegrand);
543:       }
544:     }
545:     for (nadj=0; nadj<ts->numcost; nadj++) {
546:       MatMultTranspose(J,VecSensiTemp[nadj],VecDeltaLam[nadj*s+i]);
547:       if (ts->vec_costintegral) {
548:         VecAXPY(VecDeltaLam[nadj*s+i],-h*b[i],ts->vecs_drdy[nadj]);
549:       }
550:     }

552:     /* Stage values of mu */
553:     if(ts->vecs_sensip) {
554:       TSAdjointComputeRHSJacobian(ts,rk->stage_time,Y[i],ts->Jacp);
555:       if (ts->vec_costintegral) {
556:         TSAdjointComputeDRDPFunction(ts,rk->stage_time,Y[i],ts->vecs_drdp);
557:       }

559:       for (nadj=0; nadj<ts->numcost; nadj++) {
560:         MatMultTranspose(ts->Jacp,VecSensiTemp[nadj],VecDeltaMu[nadj*s+i]);
561:         if (ts->vec_costintegral) {
562:           VecAXPY(VecDeltaMu[nadj*s+i],-h*b[i],ts->vecs_drdp[nadj]);
563:         }
564:       }
565:     }
566:   }

568:   for (j=0; j<s; j++) w[j] = 1.0;
569:   for (nadj=0; nadj<ts->numcost; nadj++) {
570:     VecMAXPY(ts->vecs_sensi[nadj],s,w,&VecDeltaLam[nadj*s]);
571:     if(ts->vecs_sensip) {
572:       VecMAXPY(ts->vecs_sensip[nadj],s,w,&VecDeltaMu[nadj*s]);
573:     }
574:   }
575:   ts->ptime += ts->time_step;
576:   ts->steps++;
577:   rk->status = TS_STEP_COMPLETE;
578:   return(0);
579: }

583: static PetscErrorCode TSInterpolate_RK(TS ts,PetscReal itime,Vec X)
584: {
585:   TS_RK           *rk = (TS_RK*)ts->data;
586:   PetscInt         s  = rk->tableau->s,pinterp = rk->tableau->pinterp,i,j;
587:   PetscReal        h;
588:   PetscReal        tt,t;
589:   PetscScalar     *b;
590:   const PetscReal *B = rk->tableau->binterp;
591:   PetscErrorCode   ierr;

594:   if (!B) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"TSRK %s does not have an interpolation formula",rk->tableau->name);
595:   switch (rk->status) {
596:   case TS_STEP_INCOMPLETE:
597:   case TS_STEP_PENDING:
598:     h = ts->time_step;
599:     t = (itime - ts->ptime)/h;
600:     break;
601:   case TS_STEP_COMPLETE:
602:     h = ts->time_step_prev;
603:     t = (itime - ts->ptime)/h + 1; /* In the interval [0,1] */
604:     break;
605:   default: SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_PLIB,"Invalid TSStepStatus");
606:   }
607:   PetscMalloc1(s,&b);
608:   for (i=0; i<s; i++) b[i] = 0;
609:   for (j=0,tt=t; j<pinterp; j++,tt*=t) {
610:     for (i=0; i<s; i++) {
611:       b[i]  += h * B[i*pinterp+j] * tt;
612:     }
613:   }
614:   VecCopy(rk->Y[0],X);
615:   VecMAXPY(X,s,b,rk->YdotRHS);
616:   PetscFree(b);
617:   return(0);
618: }

620: /*------------------------------------------------------------*/
623: static PetscErrorCode TSReset_RK(TS ts)
624: {
625:   TS_RK         *rk = (TS_RK*)ts->data;
626:   PetscInt       s;

630:   if (!rk->tableau) return(0);
631:   s    = rk->tableau->s;
632:   VecDestroyVecs(s,&rk->Y);
633:   VecDestroyVecs(s,&rk->YdotRHS);
634:   VecDestroyVecs(s*ts->numcost,&rk->VecDeltaLam);
635:   VecDestroyVecs(s*ts->numcost,&rk->VecDeltaMu);
636:   VecDestroyVecs(ts->numcost,&rk->VecSensiTemp);
637:   PetscFree(rk->work);
638:   return(0);
639: }

643: static PetscErrorCode TSDestroy_RK(TS ts)
644: {

648:   TSReset_RK(ts);
649:   PetscFree(ts->data);
650:   PetscObjectComposeFunction((PetscObject)ts,"TSRKGetType_C",NULL);
651:   PetscObjectComposeFunction((PetscObject)ts,"TSRKSetType_C",NULL);
652:   return(0);
653: }


658: static PetscErrorCode DMCoarsenHook_TSRK(DM fine,DM coarse,void *ctx)
659: {
661:   return(0);
662: }

666: static PetscErrorCode DMRestrictHook_TSRK(DM fine,Mat restrct,Vec rscale,Mat inject,DM coarse,void *ctx)
667: {
669:   return(0);
670: }


675: static PetscErrorCode DMSubDomainHook_TSRK(DM dm,DM subdm,void *ctx)
676: {
678:   return(0);
679: }

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

687:   return(0);
688: }
689: /*
692: static PetscErrorCode RKSetAdjCoe(RKTableau tab)
693: {
694:   PetscReal *A,*b;
695:   PetscInt        s,i,j;
696:   PetscErrorCode  ierr;

699:   s    = tab->s;
700:   PetscMalloc2(s*s,&A,s,&b);

702:   for (i=0; i<s; i++)
703:     for (j=0; j<s; j++) {
704:       A[i*s+j] = (tab->b[s-1-i]==0) ? 0: -tab->A[s-1-i+(s-1-j)*s] * tab->b[s-1-j] / tab->b[s-1-i];
705:       PetscPrintf(PETSC_COMM_WORLD,"Coefficients: A[%D][%D]=%.6f\n",i,j,A[i*s+j]);
706:     }

708:   for (i=0; i<s; i++) b[i] = (tab->b[s-1-i]==0)? 0: -tab->b[s-1-i];

710:   PetscMemcpy(tab->A,A,s*s*sizeof(A[0]));
711:   PetscMemcpy(tab->b,b,s*sizeof(b[0]));
712:   PetscFree2(A,b);
713:   return(0);
714: }
715: */
718: static PetscErrorCode TSSetUp_RK(TS ts)
719: {
720:   TS_RK         *rk = (TS_RK*)ts->data;
721:   RKTableau      tab;
722:   PetscInt       s;
724:   DM             dm;

727:   if (!rk->tableau) {
728:     TSRKSetType(ts,TSRKDefault);
729:   }
730:   tab  = rk->tableau;
731:   s    = tab->s;
732:   VecDuplicateVecs(ts->vec_sol,s,&rk->Y);
733:   VecDuplicateVecs(ts->vec_sol,s,&rk->YdotRHS);
734:   PetscMalloc1(s,&rk->work);
735:   TSGetDM(ts,&dm);
736:   if (dm) {
737:     DMCoarsenHookAdd(dm,DMCoarsenHook_TSRK,DMRestrictHook_TSRK,ts);
738:     DMSubDomainHookAdd(dm,DMSubDomainHook_TSRK,DMSubDomainRestrictHook_TSRK,ts);
739:   }
740:   return(0);
741: }


744: /*------------------------------------------------------------*/

748: static PetscErrorCode TSSetFromOptions_RK(PetscOptions *PetscOptionsObject,TS ts)
749: {
751:   char           rktype[256];

754:   PetscOptionsHead(PetscOptionsObject,"RK ODE solver options");
755:   {
756:     RKTableauLink  link;
757:     PetscInt       count,choice;
758:     PetscBool      flg;
759:     const char   **namelist;
760:     PetscStrncpy(rktype,TSRKDefault,sizeof(rktype));
761:     for (link=RKTableauList,count=0; link; link=link->next,count++) ;
762:     PetscMalloc1(count,&namelist);
763:     for (link=RKTableauList,count=0; link; link=link->next,count++) namelist[count] = link->tab.name;
764:     PetscOptionsEList("-ts_rk_type","Family of RK method","TSRKSetType",(const char*const*)namelist,count,rktype,&choice,&flg);
765:     TSRKSetType(ts,flg ? namelist[choice] : rktype);
766:     PetscFree(namelist);
767:   }
768:   PetscOptionsTail();
769:   return(0);
770: }

774: static PetscErrorCode PetscFormatRealArray(char buf[],size_t len,const char *fmt,PetscInt n,const PetscReal x[])
775: {
777:   PetscInt       i;
778:   size_t         left,count;
779:   char           *p;

782:   for (i=0,p=buf,left=len; i<n; i++) {
783:     PetscSNPrintfCount(p,left,fmt,&count,x[i]);
784:     if (count >= left) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Insufficient space in buffer");
785:     left -= count;
786:     p    += count;
787:     *p++  = ' ';
788:   }
789:   p[i ? 0 : -1] = 0;
790:   return(0);
791: }

795: static PetscErrorCode TSView_RK(TS ts,PetscViewer viewer)
796: {
797:   TS_RK         *rk   = (TS_RK*)ts->data;
798:   RKTableau      tab  = rk->tableau;
799:   PetscBool      iascii;
801:   TSAdapt        adapt;

804:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
805:   if (iascii) {
806:     TSRKType rktype;
807:     char     buf[512];
808:     TSRKGetType(ts,&rktype);
809:     PetscViewerASCIIPrintf(viewer,"  RK %s\n",rktype);
810:     PetscFormatRealArray(buf,sizeof(buf),"% 8.6f",tab->s,tab->c);
811:     PetscViewerASCIIPrintf(viewer,"  Abscissa     c = %s\n",buf);
812:     PetscViewerASCIIPrintf(viewer,"FSAL: %s\n",tab->FSAL ? "yes" : "no");
813:   }
814:   TSGetAdapt(ts,&adapt);
815:   TSAdaptView(adapt,viewer);
816:   return(0);
817: }

821: static PetscErrorCode TSLoad_RK(TS ts,PetscViewer viewer)
822: {
824:   TSAdapt        tsadapt;

827:   TSGetAdapt(ts,&tsadapt);
828:   TSAdaptLoad(tsadapt,viewer);
829:   return(0);
830: }

834: /*@C
835:   TSRKSetType - Set the type of RK scheme

837:   Logically collective

839:   Input Parameter:
840: +  ts - timestepping context
841: -  rktype - type of RK-scheme

843:   Level: intermediate

845: .seealso: TSRKGetType(), TSRK, TSRK2, TSRK3, TSRKPRSSP2, TSRK4, TSRK5
846: @*/
847: PetscErrorCode TSRKSetType(TS ts,TSRKType rktype)
848: {

853:   PetscTryMethod(ts,"TSRKSetType_C",(TS,TSRKType),(ts,rktype));
854:   return(0);
855: }

859: /*@C
860:   TSRKGetType - Get the type of RK scheme

862:   Logically collective

864:   Input Parameter:
865: .  ts - timestepping context

867:   Output Parameter:
868: .  rktype - type of RK-scheme

870:   Level: intermediate

872: .seealso: TSRKGetType()
873: @*/
874: PetscErrorCode TSRKGetType(TS ts,TSRKType *rktype)
875: {

880:   PetscUseMethod(ts,"TSRKGetType_C",(TS,TSRKType*),(ts,rktype));
881:   return(0);
882: }

886: PetscErrorCode  TSRKGetType_RK(TS ts,TSRKType *rktype)
887: {
888:   TS_RK     *rk = (TS_RK*)ts->data;

892:   if (!rk->tableau) {
893:     TSRKSetType(ts,TSRKDefault);
894:   }
895:   *rktype = rk->tableau->name;
896:   return(0);
897: }
900: PetscErrorCode  TSRKSetType_RK(TS ts,TSRKType rktype)
901: {
902:   TS_RK     *rk = (TS_RK*)ts->data;
904:   PetscBool      match;
905:   RKTableauLink link;

908:   if (rk->tableau) {
909:     PetscStrcmp(rk->tableau->name,rktype,&match);
910:     if (match) return(0);
911:   }
912:   for (link = RKTableauList; link; link=link->next) {
913:     PetscStrcmp(link->tab.name,rktype,&match);
914:     if (match) {
915:       TSReset_RK(ts);
916:       rk->tableau = &link->tab;
917:       return(0);
918:     }
919:   }
920:   SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_UNKNOWN_TYPE,"Could not find '%s'",rktype);
921:   return(0);
922: }

926: static PetscErrorCode  TSGetStages_RK(TS ts,PetscInt *ns,Vec **Y)
927: {
928:   TS_RK          *rk = (TS_RK*)ts->data;

931:   *ns = rk->tableau->s;
932:   if(Y) *Y  = rk->Y;
933:   return(0);
934: }


937: /* ------------------------------------------------------------ */
938: /*MC
939:       TSRK - ODE and DAE solver using Runge-Kutta schemes

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

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

947:   Level: beginner

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

952: M*/
955: PETSC_EXTERN PetscErrorCode TSCreate_RK(TS ts)
956: {
957:   TS_RK     *th;

961: #if !defined(PETSC_USE_DYNAMIC_LIBRARIES)
962:   TSRKInitializePackage();
963: #endif

965:   ts->ops->reset          = TSReset_RK;
966:   ts->ops->destroy        = TSDestroy_RK;
967:   ts->ops->view           = TSView_RK;
968:   ts->ops->load           = TSLoad_RK;
969:   ts->ops->setup          = TSSetUp_RK;
970:   ts->ops->adjointsetup   = TSAdjointSetUp_RK;
971:   ts->ops->step           = TSStep_RK;
972:   ts->ops->interpolate    = TSInterpolate_RK;
973:   ts->ops->evaluatestep   = TSEvaluateStep_RK;
974:   ts->ops->setfromoptions = TSSetFromOptions_RK;
975:   ts->ops->getstages      = TSGetStages_RK;
976:   ts->ops->adjointstep    = TSAdjointStep_RK;

978:   PetscNewLog(ts,&th);
979:   ts->data = (void*)th;

981:   PetscObjectComposeFunction((PetscObject)ts,"TSRKGetType_C",TSRKGetType_RK);
982:   PetscObjectComposeFunction((PetscObject)ts,"TSRKSetType_C",TSRKSetType_RK);
983:   return(0);
984: }