Actual source code: rk.c

petsc-master 2019-05-21
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: */

 11:  #include <petsc/private/tsimpl.h>
 12:  #include <petscdm.h>
 13:  #include <../src/ts/impls/explicit/rk/rk.h>
 14:  #include <../src/ts/impls/explicit/rk/mrk.h>

 16: static TSRKType  TSRKDefault = TSRK3BS;
 17: static PetscBool TSRKRegisterAllCalled;
 18: static PetscBool TSRKPackageInitialized;

 20: static RKTableauLink RKTableauList;

 22: /*MC
 23:      TSRK1FE - First order forward Euler scheme.

 25:      This method has one stage.

 27:      Options database:
 28: .     -ts_rk_type 1fe

 30:      Level: advanced

 32: .seealso: TSRK, TSRKType, TSRKSetType()
 33: M*/
 34: /*MC
 35:      TSRK2A - Second order RK scheme.

 37:      This method has two stages.

 39:      Options database:
 40: .     -ts_rk_type 2a

 42:      Level: advanced

 44: .seealso: TSRK, TSRKType, TSRKSetType()
 45: M*/
 46: /*MC
 47:      TSRK3 - Third order RK scheme.

 49:      This method has three stages.

 51:      Options database:
 52: .     -ts_rk_type 3

 54:      Level: advanced

 56: .seealso: TSRK, TSRKType, TSRKSetType()
 57: M*/
 58: /*MC
 59:      TSRK3BS - Third order RK scheme of Bogacki-Shampine with 2nd order embedded method.

 61:      This method has four stages with the First Same As Last (FSAL) property.

 63:      Options database:
 64: .     -ts_rk_type 3bs

 66:      Level: advanced

 68:      References: https://doi.org/10.1016/0893-9659(89)90079-7

 70: .seealso: TSRK, TSRKType, TSRKSetType()
 71: M*/
 72: /*MC
 73:      TSRK4 - Fourth order RK scheme.

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

 77:      Options database:
 78: .     -ts_rk_type 4

 80:      Level: advanced

 82: .seealso: TSRK, TSRKType, TSRKSetType()
 83: M*/
 84: /*MC
 85:      TSRK5F - Fifth order Fehlberg RK scheme with a 4th order embedded method.

 87:      This method has six stages.

 89:      Options database:
 90: .     -ts_rk_type 5f

 92:      Level: advanced

 94: .seealso: TSRK, TSRKType, TSRKSetType()
 95: M*/
 96: /*MC
 97:      TSRK5DP - Fifth order Dormand-Prince RK scheme with the 4th order embedded method.

 99:      This method has seven stages with the First Same As Last (FSAL) property.

101:      Options database:
102: .     -ts_rk_type 5dp

104:      Level: advanced

106:      References: https://doi.org/10.1016/0771-050X(80)90013-3

108: .seealso: TSRK, TSRKType, TSRKSetType()
109: M*/
110: /*MC
111:      TSRK5BS - Fifth order Bogacki-Shampine RK scheme with 4th order embedded method.

113:      This method has eight stages with the First Same As Last (FSAL) property.

115:      Options database:
116: .     -ts_rk_type 5bs

118:      Level: advanced

120:      References: https://doi.org/10.1016/0898-1221(96)00141-1

122: .seealso: TSRK, TSRKType, TSRKSetType()
123: M*/

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

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

130:   Level: advanced

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

134: .seealso:  TSRKRegisterDestroy()
135: @*/
136: PetscErrorCode TSRKRegisterAll(void)
137: {

141:   if (TSRKRegisterAllCalled) return(0);
142:   TSRKRegisterAllCalled = PETSC_TRUE;

144: #define RC PetscRealConstant
145:   {
146:     const PetscReal
147:       A[1][1] = {{0}},
148:       b[1]    = {RC(1.0)};
149:     TSRKRegister(TSRK1FE,1,1,&A[0][0],b,NULL,NULL,0,NULL);
150:   }
151:   {
152:     const PetscReal
153:       A[2][2]   = {{0,0},
154:                    {RC(1.0),0}},
155:       b[2]      =  {RC(0.5),RC(0.5)},
156:       bembed[2] =  {RC(1.0),0};
157:     TSRKRegister(TSRK2A,2,2,&A[0][0],b,NULL,bembed,0,NULL);
158:   }
159:   {
160:     const PetscReal
161:       A[3][3] = {{0,0,0},
162:                  {RC(2.0)/RC(3.0),0,0},
163:                  {RC(-1.0)/RC(3.0),RC(1.0),0}},
164:       b[3]    =  {RC(0.25),RC(0.5),RC(0.25)};
165:     TSRKRegister(TSRK3,3,3,&A[0][0],b,NULL,NULL,0,NULL);
166:   }
167:   {
168:     const PetscReal
169:       A[4][4]   = {{0,0,0,0},
170:                    {RC(1.0)/RC(2.0),0,0,0},
171:                    {0,RC(3.0)/RC(4.0),0,0},
172:                    {RC(2.0)/RC(9.0),RC(1.0)/RC(3.0),RC(4.0)/RC(9.0),0}},
173:       b[4]      =  {RC(2.0)/RC(9.0),RC(1.0)/RC(3.0),RC(4.0)/RC(9.0),0},
174:       bembed[4] =  {RC(7.0)/RC(24.0),RC(1.0)/RC(4.0),RC(1.0)/RC(3.0),RC(1.0)/RC(8.0)};
175:     TSRKRegister(TSRK3BS,3,4,&A[0][0],b,NULL,bembed,0,NULL);
176:   }
177:   {
178:     const PetscReal
179:       A[4][4] = {{0,0,0,0},
180:                  {RC(0.5),0,0,0},
181:                  {0,RC(0.5),0,0},
182:                  {0,0,RC(1.0),0}},
183:       b[4]    =  {RC(1.0)/RC(6.0),RC(1.0)/RC(3.0),RC(1.0)/RC(3.0),RC(1.0)/RC(6.0)};
184:     TSRKRegister(TSRK4,4,4,&A[0][0],b,NULL,NULL,0,NULL);
185:   }
186:   {
187:     const PetscReal
188:       A[6][6]   = {{0,0,0,0,0,0},
189:                    {RC(0.25),0,0,0,0,0},
190:                    {RC(3.0)/RC(32.0),RC(9.0)/RC(32.0),0,0,0,0},
191:                    {RC(1932.0)/RC(2197.0),RC(-7200.0)/RC(2197.0),RC(7296.0)/RC(2197.0),0,0,0},
192:                    {RC(439.0)/RC(216.0),RC(-8.0),RC(3680.0)/RC(513.0),RC(-845.0)/RC(4104.0),0,0},
193:                    {RC(-8.0)/RC(27.0),RC(2.0),RC(-3544.0)/RC(2565.0),RC(1859.0)/RC(4104.0),RC(-11.0)/RC(40.0),0}},
194:       b[6]      =  {RC(16.0)/RC(135.0),0,RC(6656.0)/RC(12825.0),RC(28561.0)/RC(56430.0),RC(-9.0)/RC(50.0),RC(2.0)/RC(55.0)},
195:       bembed[6] =  {RC(25.0)/RC(216.0),0,RC(1408.0)/RC(2565.0),RC(2197.0)/RC(4104.0),RC(-1.0)/RC(5.0),0};
196:     TSRKRegister(TSRK5F,5,6,&A[0][0],b,NULL,bembed,0,NULL);
197:   }
198:   {
199:     const PetscReal
200:       A[7][7]   = {{0,0,0,0,0,0,0},
201:                    {RC(1.0)/RC(5.0),0,0,0,0,0,0},
202:                    {RC(3.0)/RC(40.0),RC(9.0)/RC(40.0),0,0,0,0,0},
203:                    {RC(44.0)/RC(45.0),RC(-56.0)/RC(15.0),RC(32.0)/RC(9.0),0,0,0,0},
204:                    {RC(19372.0)/RC(6561.0),RC(-25360.0)/RC(2187.0),RC(64448.0)/RC(6561.0),RC(-212.0)/RC(729.0),0,0,0},
205:                    {RC(9017.0)/RC(3168.0),RC(-355.0)/RC(33.0),RC(46732.0)/RC(5247.0),RC(49.0)/RC(176.0),RC(-5103.0)/RC(18656.0),0,0},
206:                    {RC(35.0)/RC(384.0),0,RC(500.0)/RC(1113.0),RC(125.0)/RC(192.0),RC(-2187.0)/RC(6784.0),RC(11.0)/RC(84.0),0}},
207:       b[7]      =  {RC(35.0)/RC(384.0),0,RC(500.0)/RC(1113.0),RC(125.0)/RC(192.0),RC(-2187.0)/RC(6784.0),RC(11.0)/RC(84.0),0},
208:         bembed[7] =  {RC(5179.0)/RC(57600.0),0,RC(7571.0)/RC(16695.0),RC(393.0)/RC(640.0),RC(-92097.0)/RC(339200.0),RC(187.0)/RC(2100.0),RC(1.0)/RC(40.0)},
209:         binterp[7][5] =  {{RC(1.0),RC(-4034104133.0)/RC(1410260304.0),RC(105330401.0)/RC(33982176.0),RC(-13107642775.0)/RC(11282082432.0),RC(6542295.0)/RC(470086768.0)},
210:                     {0,0,0,0,0},
211:                     {0,RC(132343189600.0)/RC(32700410799.0),RC(-833316000.0)/RC(131326951.0),RC(91412856700.0)/RC(32700410799.0),RC(-523383600.0)/RC(10900136933.0)},
212:                     {0,RC(-115792950.0)/RC(29380423.0),RC(185270875.0)/RC(16991088.0),RC(-12653452475.0)/RC(1880347072.0),RC(98134425.0)/RC(235043384.0)},
213:                     {0,RC(70805911779.0)/RC(24914598704.0),RC(-4531260609.0)/RC(600351776.0),RC(988140236175.0)/RC(199316789632.0),RC(-14307999165.0)/RC(24914598704.0)},
214:                     {0,RC(-331320693.0)/RC(205662961.0),RC(31361737.0)/RC(7433601.0),RC(-2426908385.0)/RC(822651844.0),RC(97305120.0)/RC(205662961.0)},
215:                     {0,RC(44764047.0)/RC(29380423.0),RC(-1532549.0)/RC(353981.0),RC(90730570.0)/RC(29380423.0),RC(-8293050.0)/RC(29380423.0)}};

217:         TSRKRegister(TSRK5DP,5,7,&A[0][0],b,NULL,bembed,5,binterp[0]);
218:   }
219:   {
220:     const PetscReal
221:       A[8][8]   = {{0,0,0,0,0,0,0,0},
222:                    {RC(1.0)/RC(6.0),0,0,0,0,0,0,0},
223:                    {RC(2.0)/RC(27.0),RC(4.0)/RC(27.0),0,0,0,0,0,0},
224:                    {RC(183.0)/RC(1372.0),RC(-162.0)/RC(343.0),RC(1053.0)/RC(1372.0),0,0,0,0,0},
225:                    {RC(68.0)/RC(297.0),RC(-4.0)/RC(11.0),RC(42.0)/RC(143.0),RC(1960.0)/RC(3861.0),0,0,0,0},
226:                    {RC(597.0)/RC(22528.0),RC(81.0)/RC(352.0),RC(63099.0)/RC(585728.0),RC(58653.0)/RC(366080.0),RC(4617.0)/RC(20480.0),0,0,0},
227:                    {RC(174197.0)/RC(959244.0),RC(-30942.0)/RC(79937.0),RC(8152137.0)/RC(19744439.0),RC(666106.0)/RC(1039181.0),RC(-29421.0)/RC(29068.0),RC(482048.0)/RC(414219.0),0,0},
228:                    {RC(587.0)/RC(8064.0),0,RC(4440339.0)/RC(15491840.0),RC(24353.0)/RC(124800.0),RC(387.0)/RC(44800.0),RC(2152.0)/RC(5985.0),RC(7267.0)/RC(94080.0),0}},
229:       b[8]      =  {RC(587.0)/RC(8064.0),0,RC(4440339.0)/RC(15491840.0),RC(24353.0)/RC(124800.0),RC(387.0)/RC(44800.0),RC(2152.0)/RC(5985.0),RC(7267.0)/RC(94080.0),0},
230:       bembed[8] =  {RC(2479.0)/RC(34992.0),0,RC(123.0)/RC(416.0),RC(612941.0)/RC(3411720.0),RC(43.0)/RC(1440.0),RC(2272.0)/RC(6561.0),RC(79937.0)/RC(1113912.0),RC(3293.0)/RC(556956.0)};
231:     TSRKRegister(TSRK5BS,5,8,&A[0][0],b,NULL,bembed,0,NULL);
232:   }
233: #undef RC
234:   return(0);
235: }

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

240:    Not Collective

242:    Level: advanced

244: .keywords: TSRK, register, destroy
245: .seealso: TSRKRegister(), TSRKRegisterAll()
246: @*/
247: PetscErrorCode TSRKRegisterDestroy(void)
248: {
250:   RKTableauLink  link;

253:   while ((link = RKTableauList)) {
254:     RKTableau t = &link->tab;
255:     RKTableauList = link->next;
256:     PetscFree3(t->A,t->b,t->c);
257:     PetscFree (t->bembed);
258:     PetscFree (t->binterp);
259:     PetscFree (t->name);
260:     PetscFree (link);
261:   }
262:   TSRKRegisterAllCalled = PETSC_FALSE;
263:   return(0);
264: }

266: /*@C
267:   TSRKInitializePackage - This function initializes everything in the TSRK package. It is called
268:   from TSInitializePackage().

270:   Level: developer

272: .keywords: TS, TSRK, initialize, package
273: .seealso: PetscInitialize()
274: @*/
275: PetscErrorCode TSRKInitializePackage(void)
276: {

280:   if (TSRKPackageInitialized) return(0);
281:   TSRKPackageInitialized = PETSC_TRUE;
282:   TSRKRegisterAll();
283:   PetscRegisterFinalize(TSRKFinalizePackage);
284:   return(0);
285: }

287: /*@C
288:   TSRKFinalizePackage - This function destroys everything in the TSRK package. It is
289:   called from PetscFinalize().

291:   Level: developer

293: .keywords: Petsc, destroy, package
294: .seealso: PetscFinalize()
295: @*/
296: PetscErrorCode TSRKFinalizePackage(void)
297: {

301:   TSRKPackageInitialized = PETSC_FALSE;
302:   TSRKRegisterDestroy();
303:   return(0);
304: }

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

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

311:    Input Parameters:
312: +  name - identifier for method
313: .  order - approximation order of method
314: .  s - number of stages, this is the dimension of the matrices below
315: .  A - stage coefficients (dimension s*s, row-major)
316: .  b - step completion table (dimension s; NULL to use last row of A)
317: .  c - abscissa (dimension s; NULL to use row sums of A)
318: .  bembed - completion table for embedded method (dimension s; NULL if not available)
319: .  p - Order of the interpolation scheme, equal to the number of columns of binterp
320: -  binterp - Coefficients of the interpolation formula (dimension s*p; NULL to reuse b with p=1)

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

325:    Level: advanced

327: .keywords: TS, register

329: .seealso: TSRK
330: @*/
331: PetscErrorCode TSRKRegister(TSRKType name,PetscInt order,PetscInt s,
332:                             const PetscReal A[],const PetscReal b[],const PetscReal c[],
333:                             const PetscReal bembed[],PetscInt p,const PetscReal binterp[])
334: {
335:   PetscErrorCode  ierr;
336:   RKTableauLink   link;
337:   RKTableau       t;
338:   PetscInt        i,j;


348:   TSRKInitializePackage();
349:   PetscNew(&link);
350:   t = &link->tab;

352:   PetscStrallocpy(name,&t->name);
353:   t->order = order;
354:   t->s = s;
355:   PetscMalloc3(s*s,&t->A,s,&t->b,s,&t->c);
356:   PetscMemcpy(t->A,A,s*s*sizeof(A[0]));
357:   if (b)  { PetscMemcpy(t->b,b,s*sizeof(b[0])); }
358:   else for (i=0; i<s; i++) t->b[i] = A[(s-1)*s+i];
359:   if (c)  { PetscMemcpy(t->c,c,s*sizeof(c[0])); }
360:   else for (i=0; i<s; i++) for (j=0,t->c[i]=0; j<s; j++) t->c[i] += A[i*s+j];
361:   t->FSAL = PETSC_TRUE;
362:   for (i=0; i<s; i++) if (t->A[(s-1)*s+i] != t->b[i]) t->FSAL = PETSC_FALSE;

364:   if (bembed) {
365:     PetscMalloc1(s,&t->bembed);
366:     PetscMemcpy(t->bembed,bembed,s*sizeof(bembed[0]));
367:   }

369:   if (!binterp) { p = 1; binterp = t->b; }
370:   t->p = p;
371:   PetscMalloc1(s*p,&t->binterp);
372:   PetscMemcpy(t->binterp,binterp,s*p*sizeof(binterp[0]));

374:   link->next = RKTableauList;
375:   RKTableauList = link;
376:   return(0);
377: }

379: /*
380:  This is for single-step RK method
381:  The step completion formula is

383:  x1 = x0 + h b^T YdotRHS

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

389:  x1e = x0 + h be^T YdotRHS
390:      = x1 - h b^T YdotRHS + h be^T YdotRHS
391:      = x1 + h (be - b)^T YdotRHS

393:  so we can evaluate the method with different order even after the step has been optimistically completed.
394: */
395: static PetscErrorCode TSEvaluateStep_RK(TS ts,PetscInt order,Vec X,PetscBool *done)
396: {
397:   TS_RK          *rk   = (TS_RK*)ts->data;
398:   RKTableau      tab  = rk->tableau;
399:   PetscScalar    *w    = rk->work;
400:   PetscReal      h;
401:   PetscInt       s    = tab->s,j;

405:   switch (rk->status) {
406:   case TS_STEP_INCOMPLETE:
407:   case TS_STEP_PENDING:
408:     h = ts->time_step; break;
409:   case TS_STEP_COMPLETE:
410:     h = ts->ptime - ts->ptime_prev; break;
411:   default: SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_PLIB,"Invalid TSStepStatus");
412:   }
413:   if (order == tab->order) {
414:     if (rk->status == TS_STEP_INCOMPLETE) {
415:       VecCopy(ts->vec_sol,X);
416:       for (j=0; j<s; j++) w[j] = h*tab->b[j]/rk->dtratio;
417:       VecMAXPY(X,s,w,rk->YdotRHS);
418:     } else {VecCopy(ts->vec_sol,X);}
419:     return(0);
420:   } else if (order == tab->order-1) {
421:     if (!tab->bembed) goto unavailable;
422:     if (rk->status == TS_STEP_INCOMPLETE) { /*Complete with the embedded method (be)*/
423:       VecCopy(ts->vec_sol,X);
424:       for (j=0; j<s; j++) w[j] = h*tab->bembed[j];
425:       VecMAXPY(X,s,w,rk->YdotRHS);
426:     } else {  /*Rollback and re-complete using (be-b) */
427:       VecCopy(ts->vec_sol,X);
428:       for (j=0; j<s; j++) w[j] = h*(tab->bembed[j] - tab->b[j]);
429:       VecMAXPY(X,s,w,rk->YdotRHS);
430:     }
431:     if (done) *done = PETSC_TRUE;
432:     return(0);
433:   }
434: unavailable:
435:   if (done) *done = PETSC_FALSE;
436:   else SETERRQ3(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"RK '%s' of order %D cannot evaluate step at order %D. Consider using -ts_adapt_type none or a different method that has an embedded estimate.",tab->name,tab->order,order);
437:   return(0);
438: }

440: static PetscErrorCode TSForwardCostIntegral_RK(TS ts)
441: {
442:   TS_RK           *rk = (TS_RK*)ts->data;
443:   TS              quadts = ts->quadraturets;
444:   RKTableau       tab = rk->tableau;
445:   const PetscInt  s = tab->s;
446:   const PetscReal *b = tab->b,*c = tab->c;
447:   Vec             *Y = rk->Y;
448:   PetscInt        i;
449:   PetscErrorCode  ierr;

452:   /* No need to backup quadts->vec_sol since it can be reverted in TSRollBack_RK */
453:   for (i=s-1; i>=0; i--) {
454:     /* Evolve quadrature TS solution to compute integrals */
455:     TSComputeRHSFunction(quadts,rk->ptime+rk->time_step*c[i],Y[i],ts->vec_costintegrand);
456:     VecAXPY(quadts->vec_sol,rk->time_step*b[i],ts->vec_costintegrand);
457:   }
458:   return(0);
459: }

461: static PetscErrorCode TSAdjointCostIntegral_RK(TS ts)
462: {
463:   TS_RK           *rk = (TS_RK*)ts->data;
464:   RKTableau       tab = rk->tableau;
465:   TS              quadts = ts->quadraturets;
466:   const PetscInt  s = tab->s;
467:   const PetscReal *b = tab->b,*c = tab->c;
468:   Vec             *Y = rk->Y;
469:   PetscInt        i;
470:   PetscErrorCode  ierr;

473:   for (i=s-1; i>=0; i--) {
474:     /* Evolve quadrature TS solution to compute integrals */
475:     TSComputeRHSFunction(quadts,ts->ptime+ts->time_step*(1.0-c[i]),Y[i],ts->vec_costintegrand);
476:     VecAXPY(quadts->vec_sol,-ts->time_step*b[i],ts->vec_costintegrand);
477:   }
478:   return(0);
479: }

481: static PetscErrorCode TSRollBack_RK(TS ts)
482: {
483:   TS_RK           *rk = (TS_RK*)ts->data;
484:   TS              quadts = ts->quadraturets;
485:   RKTableau       tab = rk->tableau;
486:   const PetscInt  s  = tab->s;
487:   const PetscReal *b = tab->b,*c = tab->c;
488:   PetscScalar     *w = rk->work;
489:   Vec             *Y = rk->Y,*YdotRHS = rk->YdotRHS;
490:   PetscInt        j;
491:   PetscReal       h;
492:   PetscErrorCode  ierr;

495:   switch (rk->status) {
496:   case TS_STEP_INCOMPLETE:
497:   case TS_STEP_PENDING:
498:     h = ts->time_step; break;
499:   case TS_STEP_COMPLETE:
500:     h = ts->ptime - ts->ptime_prev; break;
501:   default: SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_PLIB,"Invalid TSStepStatus");
502:   }
503:   for (j=0; j<s; j++) w[j] = -h*b[j];
504:   VecMAXPY(ts->vec_sol,s,w,YdotRHS);
505:   if (quadts && ts->costintegralfwd) {
506:     for (j=0; j<s; j++) {
507:       /* Revert the quadrature TS solution */
508:       TSComputeRHSFunction(quadts,rk->ptime+h*c[j],Y[j],ts->vec_costintegrand);
509:       VecAXPY(quadts->vec_sol,-h*b[j],ts->vec_costintegrand);
510:     }
511:   }
512:   return(0);
513: }

515: static PetscErrorCode TSForwardStep_RK(TS ts)
516: {
517:   TS_RK           *rk = (TS_RK*)ts->data;
518:   RKTableau       tab = rk->tableau;
519:   Mat             J,*MatsFwdSensipTemp = rk->MatsFwdSensipTemp;
520:   const PetscInt  s = tab->s;
521:   const PetscReal *A = tab->A,*c = tab->c,*b = tab->b;
522:   Vec             *Y = rk->Y;
523:   PetscInt        i,j;
524:   PetscReal       stage_time,h = ts->time_step;
525:   PetscBool       zero;
526:   PetscErrorCode  ierr;

529:   MatCopy(ts->mat_sensip,rk->MatFwdSensip0,SAME_NONZERO_PATTERN);
530:   TSGetRHSJacobian(ts,&J,NULL,NULL,NULL);

532:   for (i=0; i<s; i++) {
533:     stage_time = ts->ptime + h*c[i];
534:     zero = PETSC_FALSE;
535:     if (b[i] == 0 && i == s-1) zero = PETSC_TRUE;
536:     /* TLM Stage values */
537:     if(!i) {
538:       MatCopy(ts->mat_sensip,rk->MatsFwdStageSensip[i],SAME_NONZERO_PATTERN);
539:     } else if (!zero) {
540:       MatZeroEntries(rk->MatsFwdStageSensip[i]);
541:       for (j=0; j<i; j++) {
542:         MatAXPY(rk->MatsFwdStageSensip[i],h*A[i*s+j],MatsFwdSensipTemp[j],SAME_NONZERO_PATTERN);
543:       }
544:       MatAXPY(rk->MatsFwdStageSensip[i],1.,ts->mat_sensip,SAME_NONZERO_PATTERN);
545:     } else {
546:       MatZeroEntries(rk->MatsFwdStageSensip[i]);
547:     }

549:     TSComputeRHSJacobian(ts,stage_time,Y[i],J,J);
550:     MatMatMult(J,rk->MatsFwdStageSensip[i],MAT_REUSE_MATRIX,PETSC_DEFAULT,&MatsFwdSensipTemp[i]);
551:     if (ts->Jacprhs) {
552:       TSComputeRHSJacobianP(ts,stage_time,Y[i],ts->Jacprhs); /* get f_p */
553:       if (ts->vecs_sensi2p) { /* TLM used for 2nd-order adjoint */
554:         PetscScalar *xarr;
555:         MatDenseGetColumn(MatsFwdSensipTemp[i],0,&xarr);
556:         VecPlaceArray(rk->VecDeltaFwdSensipCol,xarr);
557:         MatMultAdd(ts->Jacprhs,ts->vec_dir,rk->VecDeltaFwdSensipCol,rk->VecDeltaFwdSensipCol);
558:         VecResetArray(rk->VecDeltaFwdSensipCol);
559:         MatDenseRestoreColumn(MatsFwdSensipTemp[i],&xarr);
560:       } else {
561:         MatAXPY(MatsFwdSensipTemp[i],1.,ts->Jacprhs,SUBSET_NONZERO_PATTERN);
562:       }
563:     }
564:   }

566:   for (i=0; i<s; i++) {
567:     MatAXPY(ts->mat_sensip,h*b[i],rk->MatsFwdSensipTemp[i],SAME_NONZERO_PATTERN);
568:   }
569:   rk->status = TS_STEP_COMPLETE;
570:   return(0);
571: }

573: static PetscErrorCode TSForwardGetStages_RK(TS ts,PetscInt *ns,Mat **stagesensip)
574: {
575:   TS_RK     *rk = (TS_RK*)ts->data;
576:   RKTableau tab  = rk->tableau;

579:   if (ns) *ns = tab->s;
580:   if (stagesensip) *stagesensip = rk->MatsFwdStageSensip;
581:   return(0);
582: }

584: static PetscErrorCode TSForwardSetUp_RK(TS ts)
585: {
586:   TS_RK          *rk = (TS_RK*)ts->data;
587:   RKTableau      tab  = rk->tableau;
588:   PetscInt       i;

592:   /* backup sensitivity results for roll-backs */
593:   MatDuplicate(ts->mat_sensip,MAT_DO_NOT_COPY_VALUES,&rk->MatFwdSensip0);

595:   PetscMalloc1(tab->s,&rk->MatsFwdStageSensip);
596:   PetscMalloc1(tab->s,&rk->MatsFwdSensipTemp);
597:   for(i=0; i<tab->s; i++) {
598:     MatDuplicate(ts->mat_sensip,MAT_DO_NOT_COPY_VALUES,&rk->MatsFwdStageSensip[i]);
599:     MatDuplicate(ts->mat_sensip,MAT_DO_NOT_COPY_VALUES,&rk->MatsFwdSensipTemp[i]);
600:   }
601:   VecDuplicate(ts->vec_sol,&rk->VecDeltaFwdSensipCol);
602:   return(0);
603: }

605: static PetscErrorCode TSForwardReset_RK(TS ts)
606: {
607:   TS_RK          *rk = (TS_RK*)ts->data;
608:   RKTableau      tab  = rk->tableau;
609:   PetscInt       i;

613:   MatDestroy(&rk->MatFwdSensip0);
614:   if (rk->MatsFwdStageSensip) {
615:     for (i=0; i<tab->s; i++) {
616:       MatDestroy(&rk->MatsFwdStageSensip[i]);
617:     }
618:     PetscFree(rk->MatsFwdStageSensip);
619:   }
620:   if (rk->MatsFwdSensipTemp) {
621:     for (i=0; i<tab->s; i++) {
622:       MatDestroy(&rk->MatsFwdSensipTemp[i]);
623:     }
624:     PetscFree(rk->MatsFwdSensipTemp);
625:   }
626:   VecDestroy(&rk->VecDeltaFwdSensipCol);
627:   return(0);
628: }

630: static PetscErrorCode TSStep_RK(TS ts)
631: {
632:   TS_RK           *rk  = (TS_RK*)ts->data;
633:   RKTableau       tab  = rk->tableau;
634:   const PetscInt  s = tab->s;
635:   const PetscReal *A = tab->A,*c = tab->c;
636:   PetscScalar     *w = rk->work;
637:   Vec             *Y = rk->Y,*YdotRHS = rk->YdotRHS;
638:   PetscBool       FSAL = tab->FSAL;
639:   TSAdapt         adapt;
640:   PetscInt        i,j;
641:   PetscInt        rejections = 0;
642:   PetscBool       stageok,accept = PETSC_TRUE;
643:   PetscReal       next_time_step = ts->time_step;
644:   PetscErrorCode  ierr;

647:   if (ts->steprollback || ts->steprestart) FSAL = PETSC_FALSE;
648:   if (FSAL) { VecCopy(YdotRHS[s-1],YdotRHS[0]); }

650:   rk->status = TS_STEP_INCOMPLETE;
651:   while (!ts->reason && rk->status != TS_STEP_COMPLETE) {
652:     PetscReal t = ts->ptime;
653:     PetscReal h = ts->time_step;
654:     for (i=0; i<s; i++) {
655:       rk->stage_time = t + h*c[i];
656:       TSPreStage(ts,rk->stage_time);
657:       VecCopy(ts->vec_sol,Y[i]);
658:       for (j=0; j<i; j++) w[j] = h*A[i*s+j];
659:       VecMAXPY(Y[i],i,w,YdotRHS);
660:       TSPostStage(ts,rk->stage_time,i,Y);
661:       TSGetAdapt(ts,&adapt);
662:       TSAdaptCheckStage(adapt,ts,rk->stage_time,Y[i],&stageok);
663:       if (!stageok) goto reject_step;
664:       if (FSAL && !i) continue;
665:       TSComputeRHSFunction(ts,t+h*c[i],Y[i],YdotRHS[i]);
666:     }

668:     rk->status = TS_STEP_INCOMPLETE;
669:     TSEvaluateStep(ts,tab->order,ts->vec_sol,NULL);
670:     rk->status = TS_STEP_PENDING;
671:     TSGetAdapt(ts,&adapt);
672:     TSAdaptCandidatesClear(adapt);
673:     TSAdaptCandidateAdd(adapt,tab->name,tab->order,1,tab->ccfl,(PetscReal)tab->s,PETSC_TRUE);
674:     TSAdaptChoose(adapt,ts,ts->time_step,NULL,&next_time_step,&accept);
675:     rk->status = accept ? TS_STEP_COMPLETE : TS_STEP_INCOMPLETE;
676:     if (!accept) { /* Roll back the current step */
677:       TSRollBack_RK(ts);
678:       ts->time_step = next_time_step;
679:       goto reject_step;
680:     }

682:     if (ts->costintegralfwd) { /* Save the info for the later use in cost integral evaluation */
683:       rk->ptime     = ts->ptime;
684:       rk->time_step = ts->time_step;
685:     }

687:     ts->ptime += ts->time_step;
688:     ts->time_step = next_time_step;
689:     break;

691:     reject_step:
692:     ts->reject++; accept = PETSC_FALSE;
693:     if (!ts->reason && ++rejections > ts->max_reject && ts->max_reject >= 0) {
694:       ts->reason = TS_DIVERGED_STEP_REJECTED;
695:       PetscInfo2(ts,"Step=%D, step rejections %D greater than current TS allowed, stopping solve\n",ts->steps,rejections);
696:     }
697:   }
698:   return(0);
699: }

701: static PetscErrorCode TSAdjointSetUp_RK(TS ts)
702: {
703:   TS_RK          *rk  = (TS_RK*)ts->data;
704:   RKTableau      tab = rk->tableau;
705:   PetscInt       s   = tab->s;

709:   if (ts->adjointsetupcalled++) return(0);
710:   VecDuplicateVecs(ts->vecs_sensi[0],s*ts->numcost,&rk->VecsDeltaLam);
711:   VecDuplicateVecs(ts->vecs_sensi[0],ts->numcost,&rk->VecsSensiTemp);
712:   if(ts->vecs_sensip) {
713:     VecDuplicate(ts->vecs_sensip[0],&rk->VecDeltaMu);
714:   }
715:   if (ts->vecs_sensi2) {
716:     VecDuplicateVecs(ts->vecs_sensi[0],s*ts->numcost,&rk->VecsDeltaLam2);
717:     VecDuplicateVecs(ts->vecs_sensi2[0],ts->numcost,&rk->VecsSensi2Temp);
718:   }
719:   if (ts->vecs_sensi2p) {
720:     VecDuplicate(ts->vecs_sensi2p[0],&rk->VecDeltaMu2);
721:   }
722:   return(0);
723: }

725: /*
726:   Assumptions:
727:     - TSStep_RK() always evaluates the step with b, not bembed.
728: */
729: static PetscErrorCode TSAdjointStep_RK(TS ts)
730: {
731:   TS_RK            *rk = (TS_RK*)ts->data;
732:   TS               quadts = ts->quadraturets;
733:   RKTableau        tab = rk->tableau;
734:   Mat              J,Jquad;
735:   const PetscInt   s = tab->s;
736:   const PetscReal  *A = tab->A,*b = tab->b,*c = tab->c;
737:   PetscScalar      *w = rk->work,*xarr;
738:   Vec              *Y = rk->Y,*VecsDeltaLam = rk->VecsDeltaLam,VecDeltaMu = rk->VecDeltaMu,*VecsSensiTemp = rk->VecsSensiTemp;
739:   Vec              *VecsDeltaLam2 = rk->VecsDeltaLam2,VecDeltaMu2 = rk->VecDeltaMu2,*VecsSensi2Temp = rk->VecsSensi2Temp;
740:   Vec              VecDRDUTransCol = ts->vec_drdu_col,VecDRDPTransCol = ts->vec_drdp_col;
741:   PetscInt         i,j,nadj;
742:   PetscReal        t = ts->ptime;
743:   PetscReal        h = ts->time_step,stage_time;
744:   PetscErrorCode   ierr;

747:   rk->status = TS_STEP_INCOMPLETE;

749:   TSGetRHSJacobian(ts,&J,NULL,NULL,NULL);
750:   if (quadts) {
751:     TSGetRHSJacobian(quadts,&Jquad,NULL,NULL,NULL);
752:   }
753:   for (i=s-1; i>=0; i--) {
754:     if (tab->FSAL && i == s-1) {
755:       /* VecsDeltaLam[nadj*s+s-1] are initialized with zeros and the values never change.*/
756:       continue;
757:     }
758:     stage_time = t + h*(1.0-c[i]);
759:     TSComputeRHSJacobian(ts,stage_time,Y[i],J,J);
760:     if (quadts) {
761:       TSComputeRHSJacobian(quadts,stage_time,Y[i],Jquad,Jquad); /* get r_u^T */
762:     }
763:     if (ts->vecs_sensip) {
764:       TSComputeRHSJacobianP(ts,stage_time,Y[i],ts->Jacprhs); /* get f_p */
765:       if (quadts) {
766:         TSComputeRHSJacobianP(quadts,stage_time,Y[i],quadts->Jacprhs); /* get f_p for the quadrature */
767:       }
768:     }

770:     if (b[i]) {
771:       for (j=i+1; j<s; j++) w[j-i-1] = A[j*s+i]/b[i]; /* coefficients for computing VecsSensiTemp */
772:     } else {
773:       for (j=i+1; j<s; j++) w[j-i-1] = A[j*s+i]; /* coefficients for computing VecsSensiTemp */
774:     }

776:     for (nadj=0; nadj<ts->numcost; nadj++) {
777:       /* Stage values of lambda */
778:       if (b[i]) {
779:         /* lambda_{n+1} + \sum_{j=i+1}^s a_{ji}/b[i]*lambda_{s,j} */
780:         VecCopy(ts->vecs_sensi[nadj],VecsSensiTemp[nadj]); /* VecDeltaLam is an vec array of size s by numcost */
781:         VecMAXPY(VecsSensiTemp[nadj],s-i-1,w,&VecsDeltaLam[nadj*s+i+1]);
782:         MatMultTranspose(J,VecsSensiTemp[nadj],VecsDeltaLam[nadj*s+i]); /* VecsSensiTemp will be reused by 2nd-order adjoint */
783:         VecScale(VecsDeltaLam[nadj*s+i],-h*b[i]);
784:         if (quadts) {
785:           MatDenseGetColumn(Jquad,nadj,&xarr);
786:           VecPlaceArray(VecDRDUTransCol,xarr);
787:           VecAXPY(VecsDeltaLam[nadj*s+i],-h*b[i],VecDRDUTransCol);
788:           VecResetArray(VecDRDUTransCol);
789:           MatDenseRestoreColumn(Jquad,&xarr);
790:         }
791:       } else {
792:         /* \sum_{j=i+1}^s a_{ji}*lambda_{s,j} */
793:         VecSet(VecsSensiTemp[nadj],0);
794:         VecMAXPY(VecsSensiTemp[nadj],s-i-1,w,&VecsDeltaLam[nadj*s+i+1]);
795:         MatMultTranspose(J,VecsSensiTemp[nadj],VecsDeltaLam[nadj*s+i]);
796:         VecScale(VecsDeltaLam[nadj*s+i],-h);
797:       }

799:       /* Stage values of mu */
800:       if (ts->vecs_sensip) {
801:         MatMultTranspose(ts->Jacprhs,VecsSensiTemp[nadj],VecDeltaMu);
802:         if (b[i]) {
803:           VecScale(VecDeltaMu,-h*b[i]);
804:           if (quadts) {
805:             MatDenseGetColumn(quadts->Jacprhs,nadj,&xarr);
806:             VecPlaceArray(VecDRDPTransCol,xarr);
807:             VecAXPY(VecDeltaMu,-h*b[i],VecDRDPTransCol);
808:             VecResetArray(VecDRDPTransCol);
809:             MatDenseRestoreColumn(quadts->Jacprhs,&xarr);
810:           }
811:         } else {
812:           VecScale(VecDeltaMu,-h);
813:         }
814:         VecAXPY(ts->vecs_sensip[nadj],1.,VecDeltaMu); /* update sensip for each stage */
815:       }
816:     }

818:     if (ts->vecs_sensi2 && ts->forward_solve) { /* 2nd-order adjoint, TLM mode has to be turned on */
819:       /* Get w1 at t_{n+1} from TLM matrix */
820:       MatDenseGetColumn(rk->MatsFwdStageSensip[i],0,&xarr);
821:       VecPlaceArray(ts->vec_sensip_col,xarr);
822:       /* lambda_s^T F_UU w_1 */
823:       TSComputeRHSHessianProductFunctionUU(ts,stage_time,Y[i],VecsSensiTemp,ts->vec_sensip_col,ts->vecs_guu);
824:       if (quadts)  {
825:         /* R_UU w_1 */
826:         TSComputeRHSHessianProductFunctionUU(quadts,stage_time,Y[i],NULL,ts->vec_sensip_col,ts->vecs_guu);
827:       }
828:       if (ts->vecs_sensip) {
829:         /* lambda_s^T F_UP w_2 */
830:         TSComputeRHSHessianProductFunctionUP(ts,stage_time,Y[i],VecsSensiTemp,ts->vec_dir,ts->vecs_gup);
831:         if (quadts)  {
832:           /* R_UP w_2 */
833:           TSComputeRHSHessianProductFunctionUP(quadts,stage_time,Y[i],NULL,ts->vec_sensip_col,ts->vecs_gup);
834:         }
835:       }
836:       if (ts->vecs_sensi2p) {
837:         /* lambda_s^T F_PU w_1 */
838:         TSComputeRHSHessianProductFunctionPU(ts,stage_time,Y[i],VecsSensiTemp,ts->vec_sensip_col,ts->vecs_gpu);
839:         /* lambda_s^T F_PP w_2 */
840:         TSComputeRHSHessianProductFunctionPP(ts,stage_time,Y[i],VecsSensiTemp,ts->vec_dir,ts->vecs_gpp);
841:         if (b[i] && quadts) {
842:           /* R_PU w_1 */
843:           TSComputeRHSHessianProductFunctionPU(quadts,stage_time,Y[i],NULL,ts->vec_sensip_col,ts->vecs_gpu);
844:           /* R_PP w_2 */
845:           TSComputeRHSHessianProductFunctionPP(quadts,stage_time,Y[i],NULL,ts->vec_dir,ts->vecs_gpp);
846:         }
847:       }
848:       VecResetArray(ts->vec_sensip_col);
849:       MatDenseRestoreColumn(rk->MatsFwdStageSensip[i],&xarr);

851:       for (nadj=0; nadj<ts->numcost; nadj++) {
852:         /* Stage values of lambda */
853:         if (b[i]) {
854:           /* J_i^T*(Lambda_{n+1}+\sum_{j=i+1}^s a_{ji}/b_i*Lambda_{s,j} */
855:           VecCopy(ts->vecs_sensi2[nadj],VecsSensi2Temp[nadj]);
856:           VecMAXPY(VecsSensi2Temp[nadj],s-i-1,w,&VecsDeltaLam2[nadj*s+i+1]);
857:           MatMultTranspose(J,VecsSensi2Temp[nadj],VecsDeltaLam2[nadj*s+i]);
858:           VecScale(VecsDeltaLam2[nadj*s+i],-h*b[i]);
859:           VecAXPY(VecsDeltaLam2[nadj*s+i],-h*b[i],ts->vecs_guu[nadj]);
860:           if (ts->vecs_sensip) {
861:             VecAXPY(VecsDeltaLam2[nadj*s+i],-h*b[i],ts->vecs_gup[nadj]);
862:           }
863:         } else {
864:           /* \sum_{j=i+1}^s a_{ji}*Lambda_{s,j} */
865:           VecSet(VecsDeltaLam2[nadj*s+i],0);
866:           VecMAXPY(VecsSensi2Temp[nadj],s-i-1,w,&VecsDeltaLam2[nadj*s+i+1]);
867:           MatMultTranspose(J,VecsSensi2Temp[nadj],VecsDeltaLam2[nadj*s+i]);
868:           VecScale(VecsDeltaLam2[nadj*s+i],-h);
869:           VecAXPY(VecsDeltaLam2[nadj*s+i],-h,ts->vecs_guu[nadj]);
870:           if (ts->vecs_sensip) {
871:             VecAXPY(VecsDeltaLam2[nadj*s+i],-h,ts->vecs_gup[nadj]);
872:           }
873:         }
874:         if (ts->vecs_sensi2p) { /* 2nd-order adjoint for parameters */
875:           MatMultTranspose(ts->Jacprhs,VecsSensi2Temp[nadj],VecDeltaMu2);
876:           if (b[i]) {
877:             VecScale(VecDeltaMu2,-h*b[i]);
878:             VecAXPY(VecDeltaMu2,-h*b[i],ts->vecs_gpu[nadj]);
879:             VecAXPY(VecDeltaMu2,-h*b[i],ts->vecs_gpp[nadj]);
880:           } else {
881:             VecScale(VecDeltaMu2,-h);
882:             VecAXPY(VecDeltaMu2,-h,ts->vecs_gpu[nadj]);
883:             VecAXPY(VecDeltaMu2,-h,ts->vecs_gpp[nadj]);
884:           }
885:           VecAXPY(ts->vecs_sensi2p[nadj],1,VecDeltaMu2); /* update sensi2p for each stage */
886:         }
887:       }
888:     }
889:   }

891:   for (j=0; j<s; j++) w[j] = 1.0;
892:   for (nadj=0; nadj<ts->numcost; nadj++) { /* no need to do this for mu's */
893:     VecMAXPY(ts->vecs_sensi[nadj],s,w,&VecsDeltaLam[nadj*s]);
894:     if (ts->vecs_sensi2) {
895:       VecMAXPY(ts->vecs_sensi2[nadj],s,w,&VecsDeltaLam2[nadj*s]);
896:     }
897:   }
898:   rk->status = TS_STEP_COMPLETE;
899:   return(0);
900: }

902: static PetscErrorCode TSAdjointReset_RK(TS ts)
903: {
904:   TS_RK          *rk = (TS_RK*)ts->data;
905:   RKTableau      tab = rk->tableau;

909:   VecDestroyVecs(tab->s*ts->numcost,&rk->VecsDeltaLam);
910:   VecDestroyVecs(ts->numcost,&rk->VecsSensiTemp);
911:   VecDestroy(&rk->VecDeltaMu);
912:   VecDestroyVecs(tab->s*ts->numcost,&rk->VecsDeltaLam2);
913:   VecDestroy(&rk->VecDeltaMu2);
914:   VecDestroyVecs(ts->numcost,&rk->VecsSensi2Temp);
915:   return(0);
916: }

918: static PetscErrorCode TSInterpolate_RK(TS ts,PetscReal itime,Vec X)
919: {
920:   TS_RK            *rk = (TS_RK*)ts->data;
921:   PetscInt         s  = rk->tableau->s,p = rk->tableau->p,i,j;
922:   PetscReal        h;
923:   PetscReal        tt,t;
924:   PetscScalar      *b;
925:   const PetscReal  *B = rk->tableau->binterp;
926:   PetscErrorCode   ierr;

929:   if (!B) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"TSRK %s does not have an interpolation formula",rk->tableau->name);

931:   switch (rk->status) {
932:     case TS_STEP_INCOMPLETE:
933:     case TS_STEP_PENDING:
934:       h = ts->time_step;
935:       t = (itime - ts->ptime)/h;
936:       break;
937:     case TS_STEP_COMPLETE:
938:       h = ts->ptime - ts->ptime_prev;
939:       t = (itime - ts->ptime)/h + 1; /* In the interval [0,1] */
940:       break;
941:     default: SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_PLIB,"Invalid TSStepStatus");
942:   }
943:   PetscMalloc1(s,&b);
944:   for (i=0; i<s; i++) b[i] = 0;
945:   for (j=0,tt=t; j<p; j++,tt*=t) {
946:     for (i=0; i<s; i++) {
947:       b[i]  += h * B[i*p+j] * tt;
948:     }
949:   }
950:   VecCopy(rk->Y[0],X);
951:   VecMAXPY(X,s,b,rk->YdotRHS);
952:   PetscFree(b);
953:   return(0);
954: }

956: /*------------------------------------------------------------*/

958: static PetscErrorCode TSRKTableauReset(TS ts)
959: {
960:   TS_RK          *rk = (TS_RK*)ts->data;
961:   RKTableau      tab = rk->tableau;

965:   if (!tab) return(0);
966:   PetscFree(rk->work);
967:   VecDestroyVecs(tab->s,&rk->Y);
968:   VecDestroyVecs(tab->s,&rk->YdotRHS);
969:   VecDestroyVecs(tab->s*ts->numcost,&rk->VecsDeltaLam);
970:   VecDestroyVecs(ts->numcost,&rk->VecsSensiTemp);
971:   VecDestroy(&rk->VecDeltaMu);
972:   return(0);
973: }

975: static PetscErrorCode TSReset_RK(TS ts)
976: {

980:   TSRKTableauReset(ts);
981:   if (ts->use_splitrhsfunction) {
982:     PetscTryMethod(ts,"TSReset_RK_MultirateSplit_C",(TS),(ts));
983:   } else {
984:     PetscTryMethod(ts,"TSReset_RK_MultirateNonsplit_C",(TS),(ts));
985:   }
986:   return(0);
987: }

989: static PetscErrorCode DMCoarsenHook_TSRK(DM fine,DM coarse,void *ctx)
990: {
992:   return(0);
993: }

995: static PetscErrorCode DMRestrictHook_TSRK(DM fine,Mat restrct,Vec rscale,Mat inject,DM coarse,void *ctx)
996: {
998:   return(0);
999: }


1002: static PetscErrorCode DMSubDomainHook_TSRK(DM dm,DM subdm,void *ctx)
1003: {
1005:   return(0);
1006: }

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

1012:   return(0);
1013: }
1014: /*
1015: static PetscErrorCode RKSetAdjCoe(RKTableau tab)
1016: {
1017:   PetscReal *A,*b;
1018:   PetscInt        s,i,j;
1019:   PetscErrorCode  ierr;

1022:   s    = tab->s;
1023:   PetscMalloc2(s*s,&A,s,&b);

1025:   for (i=0; i<s; i++)
1026:     for (j=0; j<s; j++) {
1027:       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];
1028:       PetscPrintf(PETSC_COMM_WORLD,"Coefficients: A[%D][%D]=%.6f\n",i,j,A[i*s+j]);
1029:     }

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

1033:   PetscMemcpy(tab->A,A,s*s*sizeof(A[0]));
1034:   PetscMemcpy(tab->b,b,s*sizeof(b[0]));
1035:   PetscFree2(A,b);
1036:   return(0);
1037: }
1038: */

1040: static PetscErrorCode TSRKTableauSetUp(TS ts)
1041: {
1042:   TS_RK          *rk  = (TS_RK*)ts->data;
1043:   RKTableau      tab = rk->tableau;

1047:   PetscMalloc1(tab->s,&rk->work);
1048:   VecDuplicateVecs(ts->vec_sol,tab->s,&rk->Y);
1049:   VecDuplicateVecs(ts->vec_sol,tab->s,&rk->YdotRHS);
1050:   return(0);
1051: }

1053: static PetscErrorCode TSSetUp_RK(TS ts)
1054: {
1055:   TS             quadts = ts->quadraturets;
1057:   DM             dm;

1060:   TSCheckImplicitTerm(ts);
1061:   TSRKTableauSetUp(ts);
1062:   if (quadts && ts->costintegralfwd) {
1063:     Mat Jquad;
1064:     TSGetRHSJacobian(quadts,&Jquad,NULL,NULL,NULL);
1065:   }
1066:   TSGetDM(ts,&dm);
1067:   DMCoarsenHookAdd(dm,DMCoarsenHook_TSRK,DMRestrictHook_TSRK,ts);
1068:   DMSubDomainHookAdd(dm,DMSubDomainHook_TSRK,DMSubDomainRestrictHook_TSRK,ts);
1069:   if (ts->use_splitrhsfunction) {
1070:     PetscTryMethod(ts,"TSSetUp_RK_MultirateSplit_C",(TS),(ts));
1071:   } else {
1072:     PetscTryMethod(ts,"TSSetUp_RK_MultirateNonsplit_C",(TS),(ts));
1073:   }
1074:   return(0);
1075: }

1077: static PetscErrorCode TSSetFromOptions_RK(PetscOptionItems *PetscOptionsObject,TS ts)
1078: {
1079:   TS_RK          *rk = (TS_RK*)ts->data;

1083:   PetscOptionsHead(PetscOptionsObject,"RK ODE solver options");
1084:   {
1085:     RKTableauLink link;
1086:     PetscInt      count,choice;
1087:     PetscBool     flg,use_multirate = PETSC_FALSE;
1088:     const char    **namelist;

1090:     for (link=RKTableauList,count=0; link; link=link->next,count++) ;
1091:     PetscMalloc1(count,(char***)&namelist);
1092:     for (link=RKTableauList,count=0; link; link=link->next,count++) namelist[count] = link->tab.name;
1093:     PetscOptionsBool("-ts_rk_multirate","Use interpolation-based multirate RK method","TSRKSetMultirate",rk->use_multirate,&use_multirate,&flg);
1094:     if (flg) {
1095:       TSRKSetMultirate(ts,use_multirate);
1096:     }
1097:     PetscOptionsEList("-ts_rk_type","Family of RK method","TSRKSetType",(const char*const*)namelist,count,rk->tableau->name,&choice,&flg);
1098:     if (flg) {TSRKSetType(ts,namelist[choice]);}
1099:     PetscFree(namelist);
1100:   }
1101:   PetscOptionsTail();
1102:   PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"Multirate methods options","");
1103:   PetscOptionsInt("-ts_rk_dtratio","time step ratio between slow and fast","",rk->dtratio,&rk->dtratio,NULL);
1104:   PetscOptionsEnd();
1105:   return(0);
1106: }

1108: static PetscErrorCode TSView_RK(TS ts,PetscViewer viewer)
1109: {
1110:   TS_RK          *rk = (TS_RK*)ts->data;
1111:   PetscBool      iascii;

1115:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
1116:   if (iascii) {
1117:     RKTableau tab  = rk->tableau;
1118:     TSRKType  rktype;
1119:     char      buf[512];
1120:     TSRKGetType(ts,&rktype);
1121:     PetscViewerASCIIPrintf(viewer,"  RK type %s\n",rktype);
1122:     PetscViewerASCIIPrintf(viewer,"  Order: %D\n",tab->order);
1123:     PetscViewerASCIIPrintf(viewer,"  FSAL property: %s\n",tab->FSAL ? "yes" : "no");
1124:     PetscFormatRealArray(buf,sizeof(buf),"% 8.6f",tab->s,tab->c);
1125:     PetscViewerASCIIPrintf(viewer,"  Abscissa c = %s\n",buf);
1126:   }
1127:   return(0);
1128: }

1130: static PetscErrorCode TSLoad_RK(TS ts,PetscViewer viewer)
1131: {
1133:   TSAdapt        adapt;

1136:   TSGetAdapt(ts,&adapt);
1137:   TSAdaptLoad(adapt,viewer);
1138:   return(0);
1139: }

1141: /*@C
1142:   TSRKSetType - Set the type of RK scheme

1144:   Logically collective

1146:   Input Parameter:
1147: +  ts - timestepping context
1148: -  rktype - type of RK-scheme

1150:   Options Database:
1151: .   -ts_rk_type - <1fe,2a,3,3bs,4,5f,5dp,5bs>

1153:   Level: intermediate

1155: .seealso: TSRKGetType(), TSRK, TSRKType, TSRK1FE, TSRK2A, TSRK3, TSRK3BS, TSRK4, TSRK5F, TSRK5DP, TSRK5BS
1156: @*/
1157: PetscErrorCode TSRKSetType(TS ts,TSRKType rktype)
1158: {

1164:   PetscTryMethod(ts,"TSRKSetType_C",(TS,TSRKType),(ts,rktype));
1165:   return(0);
1166: }

1168: /*@C
1169:   TSRKGetType - Get the type of RK scheme

1171:   Logically collective

1173:   Input Parameter:
1174: .  ts - timestepping context

1176:   Output Parameter:
1177: .  rktype - type of RK-scheme

1179:   Level: intermediate

1181: .seealso: TSRKGetType()
1182: @*/
1183: PetscErrorCode TSRKGetType(TS ts,TSRKType *rktype)
1184: {

1189:   PetscUseMethod(ts,"TSRKGetType_C",(TS,TSRKType*),(ts,rktype));
1190:   return(0);
1191: }

1193: static PetscErrorCode TSRKGetType_RK(TS ts,TSRKType *rktype)
1194: {
1195:   TS_RK *rk = (TS_RK*)ts->data;

1198:   *rktype = rk->tableau->name;
1199:   return(0);
1200: }

1202: static PetscErrorCode TSRKSetType_RK(TS ts,TSRKType rktype)
1203: {
1204:   TS_RK          *rk = (TS_RK*)ts->data;
1206:   PetscBool      match;
1207:   RKTableauLink  link;

1210:   if (rk->tableau) {
1211:     PetscStrcmp(rk->tableau->name,rktype,&match);
1212:     if (match) return(0);
1213:   }
1214:   for (link = RKTableauList; link; link=link->next) {
1215:     PetscStrcmp(link->tab.name,rktype,&match);
1216:     if (match) {
1217:       if (ts->setupcalled) {TSRKTableauReset(ts);}
1218:       rk->tableau = &link->tab;
1219:       if (ts->setupcalled) {TSRKTableauSetUp(ts);}
1220:       ts->default_adapt_type = rk->tableau->bembed ? TSADAPTBASIC : TSADAPTNONE;
1221:       return(0);
1222:     }
1223:   }
1224:   SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_UNKNOWN_TYPE,"Could not find '%s'",rktype);
1225:   return(0);
1226: }

1228: static PetscErrorCode  TSGetStages_RK(TS ts,PetscInt *ns,Vec **Y)
1229: {
1230:   TS_RK *rk = (TS_RK*)ts->data;

1233:   if (ns) *ns = rk->tableau->s;
1234:   if (Y)   *Y = rk->Y;
1235:   return(0);
1236: }

1238: static PetscErrorCode TSDestroy_RK(TS ts)
1239: {

1243:   TSReset_RK(ts);
1244:   if (ts->dm) {
1245:     DMCoarsenHookRemove(ts->dm,DMCoarsenHook_TSRK,DMRestrictHook_TSRK,ts);
1246:     DMSubDomainHookRemove(ts->dm,DMSubDomainHook_TSRK,DMSubDomainRestrictHook_TSRK,ts);
1247:   }
1248:   PetscFree(ts->data);
1249:   PetscObjectComposeFunction((PetscObject)ts,"TSRKGetType_C",NULL);
1250:   PetscObjectComposeFunction((PetscObject)ts,"TSRKSetType_C",NULL);
1251:   PetscObjectComposeFunction((PetscObject)ts,"TSRKSetMultirate_C",NULL);
1252:   PetscObjectComposeFunction((PetscObject)ts,"TSRKGetMultirate_C",NULL);
1253:   return(0);
1254: }

1256: /*@C
1257:   TSRKSetMultirate - Use the interpolation-based multirate RK method

1259:   Logically collective

1261:   Input Parameter:
1262: +  ts - timestepping context
1263: -  use_multirate - PETSC_TRUE enables the multirate RK method, sets the basic method to be RK2A and sets the ratio between slow stepsize and fast stepsize to be 2

1265:   Options Database:
1266: .   -ts_rk_multirate - <true,false>

1268:   Notes:
1269:   The multirate method requires interpolation. The default interpolation works for 1st- and 2nd- order RK, but not for high-order RKs except TSRK5DP which comes with the interpolation coeffcients (binterp).

1271:   Level: intermediate

1273: .seealso: TSRKGetMultirate()
1274: @*/
1275: PetscErrorCode TSRKSetMultirate(TS ts,PetscBool use_multirate)
1276: {

1280:   PetscTryMethod(ts,"TSRKSetMultirate_C",(TS,PetscBool),(ts,use_multirate));
1281:   return(0);
1282: }

1284: /*@C
1285:   TSRKGetMultirate - Gets whether to Use the interpolation-based multirate RK method

1287:   Not collective

1289:   Input Parameter:
1290: .  ts - timestepping context

1292:   Output Parameter:
1293: .  use_multirate - PETSC_TRUE if the multirate RK method is enabled, PETSC_FALSE otherwise

1295:   Level: intermediate

1297: .seealso: TSRKSetMultirate()
1298: @*/
1299: PetscErrorCode TSRKGetMultirate(TS ts,PetscBool *use_multirate)
1300: {

1304:   PetscUseMethod(ts,"TSRKGetMultirate_C",(TS,PetscBool*),(ts,use_multirate));
1305:   return(0);
1306: }

1308: /*MC
1309:       TSRK - ODE and DAE solver using Runge-Kutta schemes

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

1314:   Notes:
1315:   The default is TSRK3BS, it can be changed with TSRKSetType() or -ts_rk_type

1317:   Level: beginner

1319: .seealso:  TSCreate(), TS, TSSetType(), TSRKSetType(), TSRKGetType(), TSRKSetFullyImplicit(), TSRK2D, TTSRK2E, TSRK3,
1320:            TSRK4, TSRK5, TSRKPRSSP2, TSRKBPR3, TSRKType, TSRKRegister(), TSRKSetMultirate(), TSRKGetMultirate()

1322: M*/
1323: PETSC_EXTERN PetscErrorCode TSCreate_RK(TS ts)
1324: {
1325:   TS_RK          *rk;

1329:   TSRKInitializePackage();

1331:   ts->ops->reset          = TSReset_RK;
1332:   ts->ops->destroy        = TSDestroy_RK;
1333:   ts->ops->view           = TSView_RK;
1334:   ts->ops->load           = TSLoad_RK;
1335:   ts->ops->setup          = TSSetUp_RK;
1336:   ts->ops->interpolate    = TSInterpolate_RK;
1337:   ts->ops->step           = TSStep_RK;
1338:   ts->ops->evaluatestep   = TSEvaluateStep_RK;
1339:   ts->ops->rollback       = TSRollBack_RK;
1340:   ts->ops->setfromoptions = TSSetFromOptions_RK;
1341:   ts->ops->getstages      = TSGetStages_RK;

1343:   ts->ops->adjointintegral = TSAdjointCostIntegral_RK;
1344:   ts->ops->adjointsetup    = TSAdjointSetUp_RK;
1345:   ts->ops->adjointstep     = TSAdjointStep_RK;
1346:   ts->ops->adjointreset    = TSAdjointReset_RK;

1348:   ts->ops->forwardintegral = TSForwardCostIntegral_RK;
1349:   ts->ops->forwardsetup    = TSForwardSetUp_RK;
1350:   ts->ops->forwardreset    = TSForwardReset_RK;
1351:   ts->ops->forwardstep     = TSForwardStep_RK;
1352:   ts->ops->forwardgetstages= TSForwardGetStages_RK;

1354:   PetscNewLog(ts,&rk);
1355:   ts->data = (void*)rk;

1357:   PetscObjectComposeFunction((PetscObject)ts,"TSRKGetType_C",TSRKGetType_RK);
1358:   PetscObjectComposeFunction((PetscObject)ts,"TSRKSetType_C",TSRKSetType_RK);
1359:   PetscObjectComposeFunction((PetscObject)ts,"TSRKSetMultirate_C",TSRKSetMultirate_RK);
1360:   PetscObjectComposeFunction((PetscObject)ts,"TSRKGetMultirate_C",TSRKGetMultirate_RK);

1362:   TSRKSetType(ts,TSRKDefault);
1363:   rk->dtratio = 1;
1364:   return(0);
1365: }