Actual source code: mprk.c

      1: /*
      2:   Code for time stepping with the Multirate Partitioned Runge-Kutta method

      4:   Notes:
      5:   1) The general system is written as
      6:      Udot = F(t,U)
      7:      if one does not split the RHS function, but gives the indexes for both slow and fast components;
      8:   2) The general system is written as
      9:      Usdot = Fs(t,Us,Uf)
     10:      Ufdot = Ff(t,Us,Uf)
     11:      for component-wise partitioned system,
     12:      users should split the RHS function themselves and also provide the indexes for both slow and fast components.
     13:   3) To correct The confusing terminology in the paper, we use 'slow method', 'slow buffer method' and 'fast method' to denote the methods applied to 'slow region', 'slow buffer region' and 'fast region' respectively. The 'slow method' in the original paper actually means the 'slow buffer method'.
     14:   4) Why does the buffer region have to be inside the slow region? The buffer region is treated with a slow method essentially. Applying the slow method to a region with a fast characteristic time scale is apparently not a good choice.

     16:   Reference:
     17:   Emil M. Constantinescu, Adrian Sandu, Multirate Timestepping Methods for Hyperbolic Conservation Laws, Journal of Scientific Computing 2007
     18: */

     20: #include <petsc/private/tsimpl.h>
     21: #include <petscdm.h>

     23: static TSMPRKType TSMPRKDefault = TSMPRK2A22;
     24: static PetscBool  TSMPRKRegisterAllCalled;
     25: static PetscBool  TSMPRKPackageInitialized;

     27: typedef struct _MPRKTableau *MPRKTableau;
     28: struct _MPRKTableau {
     29:   char      *name;
     30:   PetscInt   order;           /* Classical approximation order of the method i */
     31:   PetscInt   sbase;           /* Number of stages in the base method*/
     32:   PetscInt   s;               /* Number of stages */
     33:   PetscInt   np;              /* Number of partitions */
     34:   PetscReal *Af, *bf, *cf;    /* Tableau for fast components */
     35:   PetscReal *Amb, *bmb, *cmb; /* Tableau for medium components */
     36:   PetscInt  *rmb;             /* Array of flags for repeated stages in medium method */
     37:   PetscReal *Asb, *bsb, *csb; /* Tableau for slow components */
     38:   PetscInt  *rsb;             /* Array of flags for repeated staged in slow method*/
     39: };
     40: typedef struct _MPRKTableauLink *MPRKTableauLink;
     41: struct _MPRKTableauLink {
     42:   struct _MPRKTableau tab;
     43:   MPRKTableauLink     next;
     44: };
     45: static MPRKTableauLink MPRKTableauList;

     47: typedef struct {
     48:   MPRKTableau  tableau;
     49:   Vec         *Y; /* States computed during the step                           */
     50:   Vec         *YdotRHS;
     51:   Vec         *YdotRHS_slow;         /* Function evaluations by slow tableau for slow components  */
     52:   Vec         *YdotRHS_slowbuffer;   /* Function evaluations by medium tableau for slow components  */
     53:   Vec         *YdotRHS_medium;       /* Function evaluations by medium tableau for medium components*/
     54:   Vec         *YdotRHS_mediumbuffer; /* Function evaluations by fast tableau for medium components  */
     55:   Vec         *YdotRHS_fast;         /* Function evaluations by fast tableau for fast components  */
     56:   PetscScalar *work_slow;            /* Scalar work_slow by slow tableau                          */
     57:   PetscScalar *work_slowbuffer;      /* Scalar work_slow by slow tableau                          */
     58:   PetscScalar *work_medium;          /* Scalar work_slow by medium tableau                        */
     59:   PetscScalar *work_mediumbuffer;    /* Scalar work_mediumbuffer by fast tableau                          */
     60:   PetscScalar *work_fast;            /* Scalar work_fast by fast tableau                          */
     61:   PetscReal    stage_time;
     62:   TSStepStatus status;
     63:   PetscReal    ptime;
     64:   PetscReal    time_step;
     65:   IS           is_slow, is_slowbuffer, is_medium, is_mediumbuffer, is_fast;
     66:   TS           subts_slow, subts_slowbuffer, subts_medium, subts_mediumbuffer, subts_fast;
     67: } TS_MPRK;

     69: static PetscErrorCode TSMPRKGenerateTableau2(PetscInt ratio, PetscInt s, const PetscReal Abase[], const PetscReal bbase[], PetscReal A1[], PetscReal b1[], PetscReal A2[], PetscReal b2[])
     70: {
     71:   PetscInt i, j, k, l;

     73:   PetscFunctionBegin;
     74:   for (k = 0; k < ratio; k++) {
     75:     /* diagonal blocks */
     76:     for (i = 0; i < s; i++)
     77:       for (j = 0; j < s; j++) {
     78:         A1[(k * s + i) * ratio * s + k * s + j] = Abase[i * s + j];
     79:         A2[(k * s + i) * ratio * s + k * s + j] = Abase[i * s + j] / ratio;
     80:       }
     81:     /* off-diagonal blocks */
     82:     for (l = 0; l < k; l++)
     83:       for (i = 0; i < s; i++)
     84:         for (j = 0; j < s; j++) A2[(k * s + i) * ratio * s + l * s + j] = bbase[j] / ratio;
     85:     for (j = 0; j < s; j++) {
     86:       b1[k * s + j] = bbase[j] / ratio;
     87:       b2[k * s + j] = bbase[j] / ratio;
     88:     }
     89:   }
     90:   PetscFunctionReturn(PETSC_SUCCESS);
     91: }

     93: static PetscErrorCode TSMPRKGenerateTableau3(PetscInt ratio, PetscInt s, const PetscReal Abase[], const PetscReal bbase[], PetscReal A1[], PetscReal b1[], PetscReal A2[], PetscReal b2[], PetscReal A3[], PetscReal b3[])
     94: {
     95:   PetscInt i, j, k, l, m, n;

     97:   PetscFunctionBegin;
     98:   for (k = 0; k < ratio; k++) { /* diagonal blocks of size ratio*s by ratio*s */
     99:     for (l = 0; l < ratio; l++) /* diagonal sub-blocks of size s by s */
    100:       for (i = 0; i < s; i++)
    101:         for (j = 0; j < s; j++) {
    102:           A1[((k * ratio + l) * s + i) * ratio * ratio * s + (k * ratio + l) * s + j] = Abase[i * s + j];
    103:           A2[((k * ratio + l) * s + i) * ratio * ratio * s + (k * ratio + l) * s + j] = Abase[i * s + j] / ratio;
    104:           A3[((k * ratio + l) * s + i) * ratio * ratio * s + (k * ratio + l) * s + j] = Abase[i * s + j] / ratio / ratio;
    105:         }
    106:     for (l = 0; l < k; l++) /* off-diagonal blocks of size ratio*s by ratio*s */
    107:       for (m = 0; m < ratio; m++)
    108:         for (n = 0; n < ratio; n++)
    109:           for (i = 0; i < s; i++)
    110:             for (j = 0; j < s; j++) {
    111:               A2[((k * ratio + m) * s + i) * ratio * ratio * s + (l * ratio + n) * s + j] = bbase[j] / ratio / ratio;
    112:               A3[((k * ratio + m) * s + i) * ratio * ratio * s + (l * ratio + n) * s + j] = bbase[j] / ratio / ratio;
    113:             }
    114:     for (m = 0; m < ratio; m++)
    115:       for (n = 0; n < m; n++) /* off-diagonal sub-blocks of size s by s in the diagonal blocks */
    116:         for (i = 0; i < s; i++)
    117:           for (j = 0; j < s; j++) A3[((k * ratio + m) * s + i) * ratio * ratio * s + (k * ratio + n) * s + j] = bbase[j] / ratio / ratio;
    118:     for (n = 0; n < ratio; n++)
    119:       for (j = 0; j < s; j++) {
    120:         b1[(k * ratio + n) * s + j] = bbase[j] / ratio / ratio;
    121:         b2[(k * ratio + n) * s + j] = bbase[j] / ratio / ratio;
    122:         b3[(k * ratio + n) * s + j] = bbase[j] / ratio / ratio;
    123:       }
    124:   }
    125:   PetscFunctionReturn(PETSC_SUCCESS);
    126: }

    128: /*MC
    129:      TSMPRK2A22 - Second Order Multirate Partitioned Runge Kutta scheme based on RK2A.

    131:      This method has four stages for slow and fast parts. The refinement factor of the stepsize is 2.
    132:      r = 2, np = 2

    134:      Options Database Key:
    135: .     -ts_mprk_type 2a22 - select this scheme

    137:      Level: advanced

    139: .seealso: [](ch_ts), `TSMPRK`, `TSMPRKType`, `TSMPRKSetType()`
    140: M*/
    141: /*MC
    142:      TSMPRK2A23 - Second Order Multirate Partitioned Runge-Kutta scheme based on RK2A.

    144:      This method has eight stages for slow and medium and fast parts. The refinement factor of the stepsize is 2.
    145:      r = 2, np = 3

    147:      Options Database Key:
    148: .     -ts_mprk_type 2a23 - select this scheme

    150:      Level: advanced

    152: .seealso: [](ch_ts), `TSMPRK`, `TSMPRKType`, `TSMPRKSetType()`
    153: M*/
    154: /*MC
    155:      TSMPRK2A32 - Second Order Multirate Partitioned Runge-Kutta scheme based on RK2A.

    157:      This method has four stages for slow and fast parts. The refinement factor of the stepsize is 3.
    158:      r = 3, np = 2

    160:      Options Database Key:
    161: .     -ts_mprk_type 2a32 - select this scheme

    163:      Level: advanced

    165: .seealso: [](ch_ts), `TSMPRK`, `TSMPRKType`, `TSMPRKSetType()`
    166: M*/
    167: /*MC
    168:      TSMPRK2A33 - Second Order Multirate Partitioned Runge-Kutta scheme based on RK2A.

    170:      This method has eight stages for slow and medium and fast parts. The refinement factor of the stepsize is 3.
    171:      r = 3, np = 3

    173:      Options Database Key:
    174: .     -ts_mprk_type 2a33- select this scheme

    176:      Level: advanced

    178: .seealso: [](ch_ts), `TSMPRK`, `TSMPRKType`, `TSMPRKSetType()`
    179: M*/
    180: /*MC
    181:      TSMPRK3P2M - Third Order Multirate Partitioned Runge-Kutta scheme.

    183:      This method has eight stages for both slow and fast parts.

    185:      Options Database Key:
    186: .     -ts_mprk_type pm3 - select this scheme

    188:      Level: advanced

    190: .seealso: [](ch_ts), `TSMPRK`, `TSMPRKType`, `TSMPRKSetType()`
    191: M*/
    192: /*MC
    193:      TSMPRKP2 - Second Order Multirate Partitioned Runge-Kutta scheme.

    195:      This method has five stages for both slow and fast parts.

    197:      Options Database Key:
    198: .     -ts_mprk_type p2 - select this scheme

    200:      Level: advanced

    202: .seealso: [](ch_ts), `TSMPRK`, `TSMPRKType`, `TSMPRKSetType()`
    203: M*/
    204: /*MC
    205:      TSMPRKP3 - Third Order Multirate Partitioned Runge-Kutta scheme.

    207:      This method has ten stages for both slow and fast parts.

    209:      Options Database Key:
    210: .     -ts_mprk_type p3 - select this scheme

    212:      Level: advanced

    214: .seealso: [](ch_ts), `TSMPRK`, `TSMPRKType`, `TSMPRKSetType()`
    215: M*/

    217: /*@C
    218:   TSMPRKRegisterAll - Registers all of the Partitioned Runge-Kutta explicit methods in `TSMPRK`

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

    222:   Level: advanced

    224: .seealso: [](ch_ts), `TSMPRK`, `TSMPRKRegisterDestroy()`
    225: @*/
    226: PetscErrorCode TSMPRKRegisterAll(void)
    227: {
    228:   PetscFunctionBegin;
    229:   if (TSMPRKRegisterAllCalled) PetscFunctionReturn(PETSC_SUCCESS);
    230:   TSMPRKRegisterAllCalled = PETSC_TRUE;

    232: #define RC PetscRealConstant
    233:   {
    234:     const PetscReal Abase[2][2] =
    235:       {
    236:         {0,       0},
    237:         {RC(1.0), 0}
    238:     },
    239:                     bbase[2] = {RC(0.5), RC(0.5)};
    240:     PetscReal Asb[4][4] = {{0}}, Af[4][4] = {{0}}, bsb[4] = {0}, bf[4] = {0};
    241:     PetscInt  rsb[4] = {0, 0, 1, 2};
    242:     PetscCall(TSMPRKGenerateTableau2(2, 2, &Abase[0][0], bbase, &Asb[0][0], bsb, &Af[0][0], bf));
    243:     PetscCall(TSMPRKRegister(TSMPRK2A22, 2, 2, 2, 1, &Asb[0][0], bsb, NULL, rsb, NULL, NULL, NULL, NULL, &Af[0][0], bf, NULL));
    244:   }
    245:   {
    246:     const PetscReal Abase[2][2] =
    247:       {
    248:         {0,       0},
    249:         {RC(1.0), 0}
    250:     },
    251:                     bbase[2] = {RC(0.5), RC(0.5)};
    252:     PetscReal Asb[8][8] = {{0}}, Amb[8][8] = {{0}}, Af[8][8] = {{0}}, bsb[8] = {0}, bmb[8] = {0}, bf[8] = {0};
    253:     PetscInt  rsb[8] = {0, 0, 1, 2, 1, 2, 1, 2}, rmb[8] = {0, 0, 1, 2, 0, 0, 5, 6};
    254:     PetscCall(TSMPRKGenerateTableau3(2, 2, &Abase[0][0], bbase, &Asb[0][0], bsb, &Amb[0][0], bmb, &Af[0][0], bf));
    255:     PetscCall(TSMPRKRegister(TSMPRK2A23, 2, 2, 2, 2, &Asb[0][0], bsb, NULL, rsb, &Amb[0][0], bmb, NULL, rmb, &Af[0][0], bf, NULL));
    256:   }
    257:   {
    258:     const PetscReal Abase[2][2] =
    259:       {
    260:         {0,       0},
    261:         {RC(1.0), 0}
    262:     },
    263:                     bbase[2] = {RC(0.5), RC(0.5)};
    264:     PetscReal Asb[6][6] = {{0}}, Af[6][6] = {{0}}, bsb[6] = {0}, bf[6] = {0};
    265:     PetscInt  rsb[6] = {0, 0, 1, 2, 1, 2};
    266:     PetscCall(TSMPRKGenerateTableau2(3, 2, &Abase[0][0], bbase, &Asb[0][0], bsb, &Af[0][0], bf));
    267:     PetscCall(TSMPRKRegister(TSMPRK2A32, 2, 2, 3, 1, &Asb[0][0], bsb, NULL, rsb, NULL, NULL, NULL, NULL, &Af[0][0], bf, NULL));
    268:   }
    269:   {
    270:     const PetscReal Abase[2][2] =
    271:       {
    272:         {0,       0},
    273:         {RC(1.0), 0}
    274:     },
    275:                     bbase[2] = {RC(0.5), RC(0.5)};
    276:     PetscReal Asb[18][18] = {{0}}, Amb[18][18] = {{0}}, Af[18][18] = {{0}}, bsb[18] = {0}, bmb[18] = {0}, bf[18] = {0};
    277:     PetscInt  rsb[18] = {0, 0, 1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1, 2}, rmb[18] = {0, 0, 1, 2, 1, 2, 0, 0, 7, 8, 7, 8, 0, 0, 13, 14, 13, 14};
    278:     PetscCall(TSMPRKGenerateTableau3(3, 2, &Abase[0][0], bbase, &Asb[0][0], bsb, &Amb[0][0], bmb, &Af[0][0], bf));
    279:     PetscCall(TSMPRKRegister(TSMPRK2A33, 2, 2, 3, 3, &Asb[0][0], bsb, NULL, rsb, &Amb[0][0], bmb, NULL, rmb, &Af[0][0], bf, NULL));
    280:   }
    281:   /*
    282:     PetscReal
    283:       Asb[8][8] = {{Abase[0][0],Abase[0][1],0,0,0,0,0,0},
    284:                    {Abase[1][0],Abase[1][1],0,0,0,0,0,0},
    285:                    {0,0,Abase[0][0],Abase[0][1],0,0,0,0},
    286:                    {0,0,Abase[1][0],Abase[1][1],0,0,0,0},
    287:                    {0,0,0,0,Abase[0][0],Abase[0][1],0,0},
    288:                    {0,0,0,0,Abase[1][0],Abase[1][1],0,0},
    289:                    {0,0,0,0,0,0,Abase[0][0],Abase[0][1]},
    290:                    {0,0,0,0,0,0,Abase[1][0],Abase[1][1]}},
    291:       Amb[8][8] = {{Abase[0][0]/m,Abase[0][1]/m,0,0,0,0,0,0},
    292:                    {Abase[1][0]/m,Abase[1][1]/m,0,0,0,0,0,0},
    293:                    {0,0,Abase[0][0]/m,Abase[0][1]/m,0,0,0,0},
    294:                    {0,0,Abase[1][0]/m,Abase[1][1]/m,0,0,0,0},
    295:                    {bbase[0]/m,bbase[1]/m,bbase[0]/m,bbase[1]/m,Abase[0][0]/m,Abase[0][1]/m,0,0},
    296:                    {bbase[0]/m,bbase[1]/m,bbase[0]/m,bbase[1]/m,Abase[1][0]/m,Abase[1][1]/m,0,0},
    297:                    {bbase[0]/m,bbase[1]/m,bbase[0]/m,bbase[1]/m,0,0,Abase[0][0]/m,Abase[0][1]/m},
    298:                    {bbase[0]/m,bbase[1]/m,bbase[0]/m,bbase[1]/m,0,0,Abase[1][0]/m,Abase[1][1]/m}},
    299:       Af[8][8] = {{Abase[0][0]/m/m,Abase[0][1]/m/m,0,0,0,0,0,0},
    300:                    {Abase[1][0]/m/m,Abase[1][1]/m/m,0,0,0,0,0,0},
    301:                    {0,0,Abase[0][0]/m/m,Abase[0][1]/m/m,0,0,0,0},
    302:                    {0,0,Abase[1][0]/m/m,Abase[1][1]/m/m,0,0,0,0},
    303:                    {bbase[0]/m/m,bbase[1]/m/m,bbase[0]/m/m,bbase[1]/m/m,Abase[0][0]/m/m,Abase[0][1]/m/m,0,0},
    304:                    {bbase[0]/m/m,bbase[1]/m/m,bbase[0]/m/m,bbase[1]/m/m,Abase[1][0]/m/m,Abase[1][1]/m/m,0,0},
    305:                    {bbase[0]/m/m,bbase[1]/m/m,bbase[0]/m/m,bbase[1]/m/m,bbase[0]/m/m,bbase[1]/m/m,Abase[0][0]/m,Abase[0][1]/m},
    306:                    {bbase[0]/m/m,bbase[1]/m/m,bbase[0]/m/m,bbase[1]/m/m,bbase[0]/m/m,bbase[1]/m/m,Abase[1][0]/m,Abase[1][1]/m}},
    307:       bsb[8]    = {bbase[0]/m/m,bbase[1]/m/m,bbase[0]/m/m,bbase[1]/m/m,bbase[0]/m/m,bbase[1]/m/m,bbase[0]/m/m,bbase[1]/m/m},
    308:       bmb[8]    = {bbase[0]/m/m,bbase[1]/m/m,bbase[0]/m/m,bbase[1]/m/m,bbase[0]/m/m,bbase[1]/m/m,bbase[0]/m,bbase[1]/m/m},
    309:       bf[8]     = {bbase[0]/m/m,bbase[1]/m/m,bbase[0]/m/m,bbase[1]/m/m,bbase[0]/m/m,bbase[1]/m,bbase[0]/m/m,bbase[1]/m/m},
    310: */
    311:   /*{
    312:       const PetscReal
    313:         As[8][8] = {{0,0,0,0,0,0,0,0},
    314:                     {RC(1.0)/RC(2.0),0,0,0,0,0,0,0},
    315:                     {RC(-1.0)/RC(6.0),RC(2.0)/RC(3.0),0,0,0,0,0,0},
    316:                     {RC(1.0)/RC(3.0),RC(-1.0)/RC(3.0),RC(1.0),0,0,0,0,0},
    317:                     {0,0,0,0,0,0,0,0},
    318:                     {0,0,0,0,RC(1.0)/RC(2.0),0,0,0},
    319:                     {0,0,0,0,RC(-1.0)/RC(6.0),RC(2.0)/RC(3.0),0,0},
    320:                     {0,0,0,0,RC(1.0)/RC(3.0),RC(-1.0)/RC(3.0),RC(1.0),0}},
    321:          A[8][8] = {{0,0,0,0,0,0,0,0},
    322:                     {RC(1.0)/RC(4.0),0,0,0,0,0,0,0},
    323:                     {RC(-1.0)/RC(12.0),RC(1.0)/RC(3.0),0,0,0,0,0,0},
    324:                     {RC(1.0)/RC(6.0),RC(-1.0)/RC(6.0),RC(1.0)/RC(2.0),0,0,0,0,0},
    325:                     {RC(1.0)/RC(12.0),RC(1.0)/RC(6.0),RC(1.0)/RC(6.0),RC(1.0)/RC(12.0),0,0,0,0},
    326:                     {RC(1.0)/RC(12.0),RC(1.0)/RC(6.0),RC(1.0)/RC(6.0),RC(1.0)/RC(12.0),RC(1.0)/RC(4.0),0,0,0},
    327:                     {RC(1.0)/RC(12.0),RC(1.0)/RC(6.0),RC(1.0)/RC(6.0),RC(1.0)/RC(12.0),RC(-1.0)/RC(12.0),RC(1.0)/RC(3.0),0,0},
    328:                     {RC(1.0)/RC(12.0),RC(1.0)/RC(6.0),RC(1.0)/RC(6.0),RC(1.0)/RC(12.0),RC(1.0)/RC(6.0),RC(-1.0)/RC(6.0),RC(1.0)/RC(2.0),0}},
    329:           bs[8] = {RC(1.0)/RC(12.0),RC(1.0)/RC(6.0),RC(1.0)/RC(6.0),RC(1.0)/RC(12.0),RC(1.0)/RC(12.0),RC(1.0)/RC(6.0),RC(1.0)/RC(6.0),RC(1.0)/RC(12.0)},
    330:            b[8] = {RC(1.0)/RC(12.0),RC(1.0)/RC(6.0),RC(1.0)/RC(6.0),RC(1.0)/RC(12.0),RC(1.0)/RC(12.0),RC(1.0)/RC(6.0),RC(1.0)/RC(6.0),RC(1.0)/RC(12.0)};
    331:            PetscCall(TSMPRKRegister(TSMPRKPM3,3,8,&As[0][0],bs,NULL,&A[0][0],b,NULL));
    332:   }*/

    334:   {
    335:     const PetscReal Asb[5][5] =
    336:       {
    337:         {0,                 0, 0, 0, 0},
    338:         {RC(1.0) / RC(2.0), 0, 0, 0, 0},
    339:         {RC(1.0) / RC(2.0), 0, 0, 0, 0},
    340:         {RC(1.0),           0, 0, 0, 0},
    341:         {RC(1.0),           0, 0, 0, 0}
    342:     },
    343:                     Af[5][5] = {{0, 0, 0, 0, 0}, {RC(1.0) / RC(2.0), 0, 0, 0, 0}, {RC(1.0) / RC(4.0), RC(1.0) / RC(4.0), 0, 0, 0}, {RC(1.0) / RC(4.0), RC(1.0) / RC(4.0), RC(1.0) / RC(2.0), 0, 0}, {RC(1.0) / RC(4.0), RC(1.0) / RC(4.0), RC(1.0) / RC(4.0), RC(1.0) / RC(4.0), 0}}, bsb[5] = {RC(1.0) / RC(2.0), 0, 0, 0, RC(1.0) / RC(2.0)}, bf[5] = {RC(1.0) / RC(4.0), RC(1.0) / RC(4.0), RC(1.0) / RC(4.0), RC(1.0) / RC(4.0), 0};
    344:     const PetscInt rsb[5] = {0, 0, 2, 0, 4};
    345:     PetscCall(TSMPRKRegister(TSMPRKP2, 2, 5, 1, 1, &Asb[0][0], bsb, NULL, rsb, NULL, NULL, NULL, NULL, &Af[0][0], bf, NULL));
    346:   }

    348:   {
    349:     const PetscReal Asb[10][10] =
    350:       {
    351:         {0,                  0, 0, 0, 0,                  0,                 0, 0, 0, 0},
    352:         {RC(1.0) / RC(4.0),  0, 0, 0, 0,                  0,                 0, 0, 0, 0},
    353:         {RC(1.0) / RC(4.0),  0, 0, 0, 0,                  0,                 0, 0, 0, 0},
    354:         {RC(1.0) / RC(2.0),  0, 0, 0, 0,                  0,                 0, 0, 0, 0},
    355:         {RC(1.0) / RC(2.0),  0, 0, 0, 0,                  0,                 0, 0, 0, 0},
    356:         {RC(-1.0) / RC(6.0), 0, 0, 0, RC(2.0) / RC(3.0),  0,                 0, 0, 0, 0},
    357:         {RC(1.0) / RC(12.0), 0, 0, 0, RC(1.0) / RC(6.0),  RC(1.0) / RC(2.0), 0, 0, 0, 0},
    358:         {RC(1.0) / RC(12.0), 0, 0, 0, RC(1.0) / RC(6.0),  RC(1.0) / RC(2.0), 0, 0, 0, 0},
    359:         {RC(1.0) / RC(3.0),  0, 0, 0, RC(-1.0) / RC(3.0), RC(1.0),           0, 0, 0, 0},
    360:         {RC(1.0) / RC(3.0),  0, 0, 0, RC(-1.0) / RC(3.0), RC(1.0),           0, 0, 0, 0}
    361:     },
    362:                     Af[10][10] = {{0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {RC(1.0) / RC(4.0), 0, 0, 0, 0, 0, 0, 0, 0, 0}, {RC(-1.0) / RC(12.0), RC(1.0) / RC(3.0), 0, 0, 0, 0, 0, 0, 0, 0}, {RC(1.0) / RC(6.0), RC(-1.0) / RC(6.0), RC(1.0) / RC(2.0), 0, 0, 0, 0, 0, 0, 0}, {RC(1.0) / RC(12.0), RC(1.0) / RC(6.0), RC(1.0) / RC(6.0), RC(1.0) / RC(12.0), 0, 0, 0, 0, 0, 0}, {RC(1.0) / RC(12.0), RC(1.0) / RC(6.0), RC(1.0) / RC(6.0), RC(1.0) / RC(12.0), 0, 0, 0, 0, 0, 0}, {RC(1.0) / RC(12.0), RC(1.0) / RC(6.0), RC(1.0) / RC(6.0), RC(1.0) / RC(12.0), 0, RC(1.0) / RC(4.0), 0, 0, 0, 0}, {RC(1.0) / RC(12.0), RC(1.0) / RC(6.0), RC(1.0) / RC(6.0), RC(1.0) / RC(12.0), 0, RC(-1.0) / RC(12.0), RC(1.0) / RC(3.0), 0, 0, 0}, {RC(1.0) / RC(12.0), RC(1.0) / RC(6.0), RC(1.0) / RC(6.0), RC(1.0) / RC(12.0), 0, RC(1.0) / RC(6.0), RC(-1.0) / RC(6.0), RC(1.0) / RC(2.0), 0, 0}, {RC(1.0) / RC(12.0), RC(1.0) / RC(6.0), RC(1.0) / RC(6.0), RC(1.0) / RC(12.0), 0, RC(1.0) / RC(12.0), RC(1.0) / RC(6.0), RC(1.0) / RC(6.0), RC(1.0) / RC(12.0), 0}}, bsb[10] = {RC(1.0) / RC(6.0), 0, 0, 0, RC(1.0) / RC(3.0), RC(1.0) / RC(3.0), 0, 0, 0, RC(1.0) / RC(6.0)}, bf[10] = {RC(1.0) / RC(12.0), RC(1.0) / RC(6.0), RC(1.0) / RC(6.0), RC(1.0) / RC(12.0), 0, RC(1.0) / RC(12.0), RC(1.0) / RC(6.0), RC(1.0) / RC(6.0), RC(1.0) / RC(12.0), 0};
    363:     const PetscInt rsb[10] = {0, 0, 2, 0, 4, 0, 0, 7, 0, 9};
    364:     PetscCall(TSMPRKRegister(TSMPRKP3, 3, 5, 2, 1, &Asb[0][0], bsb, NULL, rsb, NULL, NULL, NULL, NULL, &Af[0][0], bf, NULL));
    365:   }
    366: #undef RC
    367:   PetscFunctionReturn(PETSC_SUCCESS);
    368: }

    370: /*@C
    371:   TSMPRKRegisterDestroy - Frees the list of schemes that were registered by `TSMPRKRegister()`.

    373:   Not Collective

    375:   Level: advanced

    377: .seealso: [](ch_ts), `TSMPRK`, `TSMPRKRegister()`, `TSMPRKRegisterAll()`
    378: @*/
    379: PetscErrorCode TSMPRKRegisterDestroy(void)
    380: {
    381:   MPRKTableauLink link;

    383:   PetscFunctionBegin;
    384:   while ((link = MPRKTableauList)) {
    385:     MPRKTableau t   = &link->tab;
    386:     MPRKTableauList = link->next;
    387:     PetscCall(PetscFree3(t->Asb, t->bsb, t->csb));
    388:     PetscCall(PetscFree3(t->Amb, t->bmb, t->cmb));
    389:     PetscCall(PetscFree3(t->Af, t->bf, t->cf));
    390:     PetscCall(PetscFree(t->rsb));
    391:     PetscCall(PetscFree(t->rmb));
    392:     PetscCall(PetscFree(t->name));
    393:     PetscCall(PetscFree(link));
    394:   }
    395:   TSMPRKRegisterAllCalled = PETSC_FALSE;
    396:   PetscFunctionReturn(PETSC_SUCCESS);
    397: }

    399: /*@C
    400:   TSMPRKInitializePackage - This function initializes everything in the `TSMPRK` package. It is called
    401:   from `PetscDLLibraryRegister()` when using dynamic libraries, and on the first call to `TSCreate_MPRK()`
    402:   when using static libraries.

    404:   Level: developer

    406: .seealso: [](ch_ts), `TSMPRK`, `PetscInitialize()`
    407: @*/
    408: PetscErrorCode TSMPRKInitializePackage(void)
    409: {
    410:   PetscFunctionBegin;
    411:   if (TSMPRKPackageInitialized) PetscFunctionReturn(PETSC_SUCCESS);
    412:   TSMPRKPackageInitialized = PETSC_TRUE;
    413:   PetscCall(TSMPRKRegisterAll());
    414:   PetscCall(PetscRegisterFinalize(TSMPRKFinalizePackage));
    415:   PetscFunctionReturn(PETSC_SUCCESS);
    416: }

    418: /*@C
    419:   TSMPRKFinalizePackage - This function destroys everything in the `TSMPRK` package. It is
    420:   called from `PetscFinalize()`.

    422:   Level: developer

    424: .seealso: [](ch_ts), `TSMPRK`, `PetscFinalize()`
    425: @*/
    426: PetscErrorCode TSMPRKFinalizePackage(void)
    427: {
    428:   PetscFunctionBegin;
    429:   TSMPRKPackageInitialized = PETSC_FALSE;
    430:   PetscCall(TSMPRKRegisterDestroy());
    431:   PetscFunctionReturn(PETSC_SUCCESS);
    432: }

    434: /*@C
    435:   TSMPRKRegister - register a `TSMPRK` scheme by providing the entries in the Butcher tableau

    437:   Not Collective, but the same schemes should be registered on all processes on which they will be used, No Fortran Support

    439:   Input Parameters:
    440: + name   - identifier for method
    441: . order  - approximation order of method
    442: . sbase  - number of stages in the base methods
    443: . ratio1 - stepsize ratio at 1st level (e.g. slow/medium)
    444: . ratio2 - stepsize ratio at 2nd level (e.g. medium/fast)
    445: . Asb    - stage coefficients for slow components(dimension s*s, row-major)
    446: . bsb    - step completion table for slow components(dimension s)
    447: . csb    - abscissa for slow components(dimension s)
    448: . rsb    - array of flags for repeated stages for slow components (dimension s)
    449: . Amb    - stage coefficients for medium components(dimension s*s, row-major)
    450: . bmb    - step completion table for medium components(dimension s)
    451: . cmb    - abscissa for medium components(dimension s)
    452: . rmb    - array of flags for repeated stages for medium components (dimension s)
    453: . Af     - stage coefficients for fast components(dimension s*s, row-major)
    454: . bf     - step completion table for fast components(dimension s)
    455: - cf     - abscissa for fast components(dimension s)

    457:   Level: advanced

    459:   Note:
    460:   Several `TSMPRK` methods are provided, this function is only needed to create new methods.

    462: .seealso: [](ch_ts), `TSMPRK`
    463: @*/
    464: PetscErrorCode TSMPRKRegister(TSMPRKType name, PetscInt order, PetscInt sbase, PetscInt ratio1, PetscInt ratio2, const PetscReal Asb[], const PetscReal bsb[], const PetscReal csb[], const PetscInt rsb[], const PetscReal Amb[], const PetscReal bmb[], const PetscReal cmb[], const PetscInt rmb[], const PetscReal Af[], const PetscReal bf[], const PetscReal cf[])
    465: {
    466:   MPRKTableauLink link;
    467:   MPRKTableau     t;
    468:   PetscInt        s, i, j;

    470:   PetscFunctionBegin;
    471:   PetscAssertPointer(name, 1);
    472:   PetscAssertPointer(Asb, 6);
    473:   if (bsb) PetscAssertPointer(bsb, 7);
    474:   if (csb) PetscAssertPointer(csb, 8);
    475:   if (rsb) PetscAssertPointer(rsb, 9);
    476:   if (Amb) PetscAssertPointer(Amb, 10);
    477:   if (bmb) PetscAssertPointer(bmb, 11);
    478:   if (cmb) PetscAssertPointer(cmb, 12);
    479:   if (rmb) PetscAssertPointer(rmb, 13);
    480:   PetscAssertPointer(Af, 14);
    481:   if (bf) PetscAssertPointer(bf, 15);
    482:   if (cf) PetscAssertPointer(cf, 16);

    484:   PetscCall(PetscNew(&link));
    485:   t = &link->tab;

    487:   PetscCall(PetscStrallocpy(name, &t->name));
    488:   s        = sbase * ratio1 * ratio2; /*  this is the dimension of the matrices below */
    489:   t->order = order;
    490:   t->sbase = sbase;
    491:   t->s     = s;
    492:   t->np    = 2;

    494:   PetscCall(PetscMalloc3(s * s, &t->Af, s, &t->bf, s, &t->cf));
    495:   PetscCall(PetscArraycpy(t->Af, Af, s * s));
    496:   if (bf) {
    497:     PetscCall(PetscArraycpy(t->bf, bf, s));
    498:   } else
    499:     for (i = 0; i < s; i++) t->bf[i] = Af[(s - 1) * s + i];
    500:   if (cf) {
    501:     PetscCall(PetscArraycpy(t->cf, cf, s));
    502:   } else {
    503:     for (i = 0; i < s; i++)
    504:       for (j = 0, t->cf[i] = 0; j < s; j++) t->cf[i] += Af[i * s + j];
    505:   }

    507:   if (Amb) {
    508:     t->np = 3;
    509:     PetscCall(PetscMalloc3(s * s, &t->Amb, s, &t->bmb, s, &t->cmb));
    510:     PetscCall(PetscArraycpy(t->Amb, Amb, s * s));
    511:     if (bmb) {
    512:       PetscCall(PetscArraycpy(t->bmb, bmb, s));
    513:     } else {
    514:       for (i = 0; i < s; i++) t->bmb[i] = Amb[(s - 1) * s + i];
    515:     }
    516:     if (cmb) {
    517:       PetscCall(PetscArraycpy(t->cmb, cmb, s));
    518:     } else {
    519:       for (i = 0; i < s; i++)
    520:         for (j = 0, t->cmb[i] = 0; j < s; j++) t->cmb[i] += Amb[i * s + j];
    521:     }
    522:     if (rmb) {
    523:       PetscCall(PetscMalloc1(s, &t->rmb));
    524:       PetscCall(PetscArraycpy(t->rmb, rmb, s));
    525:     } else {
    526:       PetscCall(PetscCalloc1(s, &t->rmb));
    527:     }
    528:   }

    530:   PetscCall(PetscMalloc3(s * s, &t->Asb, s, &t->bsb, s, &t->csb));
    531:   PetscCall(PetscArraycpy(t->Asb, Asb, s * s));
    532:   if (bsb) {
    533:     PetscCall(PetscArraycpy(t->bsb, bsb, s));
    534:   } else
    535:     for (i = 0; i < s; i++) t->bsb[i] = Asb[(s - 1) * s + i];
    536:   if (csb) {
    537:     PetscCall(PetscArraycpy(t->csb, csb, s));
    538:   } else {
    539:     for (i = 0; i < s; i++)
    540:       for (j = 0, t->csb[i] = 0; j < s; j++) t->csb[i] += Asb[i * s + j];
    541:   }
    542:   if (rsb) {
    543:     PetscCall(PetscMalloc1(s, &t->rsb));
    544:     PetscCall(PetscArraycpy(t->rsb, rsb, s));
    545:   } else {
    546:     PetscCall(PetscCalloc1(s, &t->rsb));
    547:   }
    548:   link->next      = MPRKTableauList;
    549:   MPRKTableauList = link;
    550:   PetscFunctionReturn(PETSC_SUCCESS);
    551: }

    553: static PetscErrorCode TSMPRKSetSplits(TS ts)
    554: {
    555:   TS_MPRK    *mprk = (TS_MPRK *)ts->data;
    556:   MPRKTableau tab  = mprk->tableau;
    557:   DM          dm, subdm, newdm;

    559:   PetscFunctionBegin;
    560:   PetscCall(TSRHSSplitGetSubTS(ts, "slow", &mprk->subts_slow));
    561:   PetscCall(TSRHSSplitGetSubTS(ts, "fast", &mprk->subts_fast));
    562:   PetscCheck(mprk->subts_slow && mprk->subts_fast, PetscObjectComm((PetscObject)ts), PETSC_ERR_USER, "Must set up the RHSFunctions for 'slow' and 'fast' components using TSRHSSplitSetRHSFunction() or calling TSSetRHSFunction() for each sub-TS");

    564:   /* Only copy the DM */
    565:   PetscCall(TSGetDM(ts, &dm));

    567:   PetscCall(TSRHSSplitGetSubTS(ts, "slowbuffer", &mprk->subts_slowbuffer));
    568:   if (!mprk->subts_slowbuffer) {
    569:     mprk->subts_slowbuffer = mprk->subts_slow;
    570:     mprk->subts_slow       = NULL;
    571:   }
    572:   if (mprk->subts_slow) {
    573:     PetscCall(DMClone(dm, &newdm));
    574:     PetscCall(TSGetDM(mprk->subts_slow, &subdm));
    575:     PetscCall(DMCopyDMTS(subdm, newdm));
    576:     PetscCall(DMCopyDMSNES(subdm, newdm));
    577:     PetscCall(TSSetDM(mprk->subts_slow, newdm));
    578:     PetscCall(DMDestroy(&newdm));
    579:   }
    580:   PetscCall(DMClone(dm, &newdm));
    581:   PetscCall(TSGetDM(mprk->subts_slowbuffer, &subdm));
    582:   PetscCall(DMCopyDMTS(subdm, newdm));
    583:   PetscCall(DMCopyDMSNES(subdm, newdm));
    584:   PetscCall(TSSetDM(mprk->subts_slowbuffer, newdm));
    585:   PetscCall(DMDestroy(&newdm));

    587:   PetscCall(DMClone(dm, &newdm));
    588:   PetscCall(TSGetDM(mprk->subts_fast, &subdm));
    589:   PetscCall(DMCopyDMTS(subdm, newdm));
    590:   PetscCall(DMCopyDMSNES(subdm, newdm));
    591:   PetscCall(TSSetDM(mprk->subts_fast, newdm));
    592:   PetscCall(DMDestroy(&newdm));

    594:   if (tab->np == 3) {
    595:     PetscCall(TSRHSSplitGetSubTS(ts, "medium", &mprk->subts_medium));
    596:     PetscCall(TSRHSSplitGetSubTS(ts, "mediumbuffer", &mprk->subts_mediumbuffer));
    597:     if (mprk->subts_medium && !mprk->subts_mediumbuffer) {
    598:       mprk->subts_mediumbuffer = mprk->subts_medium;
    599:       mprk->subts_medium       = NULL;
    600:     }
    601:     if (mprk->subts_medium) {
    602:       PetscCall(DMClone(dm, &newdm));
    603:       PetscCall(TSGetDM(mprk->subts_medium, &subdm));
    604:       PetscCall(DMCopyDMTS(subdm, newdm));
    605:       PetscCall(DMCopyDMSNES(subdm, newdm));
    606:       PetscCall(TSSetDM(mprk->subts_medium, newdm));
    607:       PetscCall(DMDestroy(&newdm));
    608:     }
    609:     PetscCall(DMClone(dm, &newdm));
    610:     PetscCall(TSGetDM(mprk->subts_mediumbuffer, &subdm));
    611:     PetscCall(DMCopyDMTS(subdm, newdm));
    612:     PetscCall(DMCopyDMSNES(subdm, newdm));
    613:     PetscCall(TSSetDM(mprk->subts_mediumbuffer, newdm));
    614:     PetscCall(DMDestroy(&newdm));
    615:   }
    616:   PetscFunctionReturn(PETSC_SUCCESS);
    617: }

    619: /*
    620:  This if for nonsplit RHS MPRK
    621:  The step completion formula is

    623:  x1 = x0 + h b^T YdotRHS

    625: */
    626: static PetscErrorCode TSEvaluateStep_MPRK(TS ts, PetscInt order, Vec X, PetscBool *done)
    627: {
    628:   TS_MPRK     *mprk = (TS_MPRK *)ts->data;
    629:   MPRKTableau  tab  = mprk->tableau;
    630:   PetscScalar *wf   = mprk->work_fast;
    631:   PetscReal    h    = ts->time_step;
    632:   PetscInt     s    = tab->s, j;

    634:   PetscFunctionBegin;
    635:   for (j = 0; j < s; j++) wf[j] = h * tab->bf[j];
    636:   PetscCall(VecCopy(ts->vec_sol, X));
    637:   PetscCall(VecMAXPY(X, s, wf, mprk->YdotRHS));
    638:   PetscFunctionReturn(PETSC_SUCCESS);
    639: }

    641: static PetscErrorCode TSStep_MPRK(TS ts)
    642: {
    643:   TS_MPRK         *mprk = (TS_MPRK *)ts->data;
    644:   Vec             *Y = mprk->Y, *YdotRHS = mprk->YdotRHS, *YdotRHS_fast = mprk->YdotRHS_fast, *YdotRHS_slow = mprk->YdotRHS_slow, *YdotRHS_slowbuffer = mprk->YdotRHS_slowbuffer;
    645:   Vec              Yslow, Yslowbuffer, Yfast;
    646:   MPRKTableau      tab = mprk->tableau;
    647:   const PetscInt   s   = tab->s;
    648:   const PetscReal *Af = tab->Af, *cf = tab->cf, *Asb = tab->Asb, *csb = tab->csb;
    649:   PetscScalar     *wf = mprk->work_fast, *wsb = mprk->work_slowbuffer;
    650:   PetscInt         i, j;
    651:   PetscReal        next_time_step = ts->time_step, t = ts->ptime, h = ts->time_step;

    653:   PetscFunctionBegin;
    654:   for (i = 0; i < s; i++) {
    655:     mprk->stage_time = t + h * cf[i];
    656:     PetscCall(TSPreStage(ts, mprk->stage_time));
    657:     PetscCall(VecCopy(ts->vec_sol, Y[i]));

    659:     /* slow buffer region */
    660:     for (j = 0; j < i; j++) wsb[j] = h * Asb[i * s + j];
    661:     for (j = 0; j < i; j++) PetscCall(VecGetSubVector(YdotRHS[j], mprk->is_slowbuffer, &YdotRHS_slowbuffer[j]));
    662:     PetscCall(VecGetSubVector(Y[i], mprk->is_slowbuffer, &Yslowbuffer));
    663:     PetscCall(VecMAXPY(Yslowbuffer, i, wsb, mprk->YdotRHS_slowbuffer));
    664:     PetscCall(VecRestoreSubVector(Y[i], mprk->is_slowbuffer, &Yslowbuffer));
    665:     for (j = 0; j < i; j++) PetscCall(VecRestoreSubVector(YdotRHS[j], mprk->is_slowbuffer, &YdotRHS_slowbuffer[j]));
    666:     /* slow region */
    667:     if (mprk->is_slow) {
    668:       for (j = 0; j < i; j++) PetscCall(VecGetSubVector(YdotRHS[j], mprk->is_slow, &YdotRHS_slow[j]));
    669:       PetscCall(VecGetSubVector(Y[i], mprk->is_slow, &Yslow));
    670:       PetscCall(VecMAXPY(Yslow, i, wsb, mprk->YdotRHS_slow));
    671:       PetscCall(VecRestoreSubVector(Y[i], mprk->is_slow, &Yslow));
    672:       for (j = 0; j < i; j++) PetscCall(VecRestoreSubVector(YdotRHS[j], mprk->is_slow, &YdotRHS_slow[j]));
    673:     }

    675:     /* fast region */
    676:     for (j = 0; j < i; j++) wf[j] = h * Af[i * s + j];
    677:     for (j = 0; j < i; j++) PetscCall(VecGetSubVector(YdotRHS[j], mprk->is_fast, &YdotRHS_fast[j]));
    678:     PetscCall(VecGetSubVector(Y[i], mprk->is_fast, &Yfast));
    679:     PetscCall(VecMAXPY(Yfast, i, wf, mprk->YdotRHS_fast));
    680:     PetscCall(VecRestoreSubVector(Y[i], mprk->is_fast, &Yfast));
    681:     for (j = 0; j < i; j++) PetscCall(VecRestoreSubVector(YdotRHS[j], mprk->is_fast, &YdotRHS_fast[j]));
    682:     if (tab->np == 3) {
    683:       Vec         *YdotRHS_medium = mprk->YdotRHS_medium, *YdotRHS_mediumbuffer = mprk->YdotRHS_mediumbuffer;
    684:       Vec          Ymedium, Ymediumbuffer;
    685:       PetscScalar *wmb = mprk->work_mediumbuffer;

    687:       for (j = 0; j < i; j++) wmb[j] = h * tab->Amb[i * s + j];
    688:       /* medium region */
    689:       if (mprk->is_medium) {
    690:         for (j = 0; j < i; j++) PetscCall(VecGetSubVector(YdotRHS[j], mprk->is_medium, &YdotRHS_medium[j]));
    691:         PetscCall(VecGetSubVector(Y[i], mprk->is_medium, &Ymedium));
    692:         PetscCall(VecMAXPY(Ymedium, i, wmb, mprk->YdotRHS_medium));
    693:         PetscCall(VecRestoreSubVector(Y[i], mprk->is_medium, &Ymedium));
    694:         for (j = 0; j < i; j++) PetscCall(VecRestoreSubVector(YdotRHS[j], mprk->is_medium, &YdotRHS_medium[j]));
    695:       }
    696:       /* medium buffer region */
    697:       for (j = 0; j < i; j++) PetscCall(VecGetSubVector(YdotRHS[j], mprk->is_mediumbuffer, &YdotRHS_mediumbuffer[j]));
    698:       PetscCall(VecGetSubVector(Y[i], mprk->is_mediumbuffer, &Ymediumbuffer));
    699:       PetscCall(VecMAXPY(Ymediumbuffer, i, wmb, mprk->YdotRHS_mediumbuffer));
    700:       PetscCall(VecRestoreSubVector(Y[i], mprk->is_mediumbuffer, &Ymediumbuffer));
    701:       for (j = 0; j < i; j++) PetscCall(VecRestoreSubVector(YdotRHS[j], mprk->is_mediumbuffer, &YdotRHS_mediumbuffer[j]));
    702:     }
    703:     PetscCall(TSPostStage(ts, mprk->stage_time, i, Y));
    704:     /* compute the stage RHS by fast and slow tableau respectively */
    705:     PetscCall(TSComputeRHSFunction(ts, t + h * csb[i], Y[i], YdotRHS[i]));
    706:   }
    707:   PetscCall(TSEvaluateStep(ts, tab->order, ts->vec_sol, NULL));
    708:   ts->ptime += ts->time_step;
    709:   ts->time_step = next_time_step;
    710:   PetscFunctionReturn(PETSC_SUCCESS);
    711: }

    713: /*
    714:  This if for the case when split RHS is used
    715:  The step completion formula is
    716:  x1 = x0 + h b^T YdotRHS
    717: */
    718: static PetscErrorCode TSEvaluateStep_MPRKSPLIT(TS ts, PetscInt order, Vec X, PetscBool *done)
    719: {
    720:   TS_MPRK     *mprk = (TS_MPRK *)ts->data;
    721:   MPRKTableau  tab  = mprk->tableau;
    722:   Vec          Xslow, Xfast, Xslowbuffer; /* subvectors for slow and fast components in X respectively */
    723:   PetscScalar *wf = mprk->work_fast, *ws = mprk->work_slow, *wsb = mprk->work_slowbuffer;
    724:   PetscReal    h = ts->time_step;
    725:   PetscInt     s = tab->s, j, computedstages;

    727:   PetscFunctionBegin;
    728:   PetscCall(VecCopy(ts->vec_sol, X));

    730:   /* slow region */
    731:   if (mprk->is_slow) {
    732:     computedstages = 0;
    733:     for (j = 0; j < s; j++) {
    734:       if (tab->rsb[j]) ws[tab->rsb[j] - 1] += h * tab->bsb[j];
    735:       else ws[computedstages++] = h * tab->bsb[j];
    736:     }
    737:     PetscCall(VecGetSubVector(X, mprk->is_slow, &Xslow));
    738:     PetscCall(VecMAXPY(Xslow, computedstages, ws, mprk->YdotRHS_slow));
    739:     PetscCall(VecRestoreSubVector(X, mprk->is_slow, &Xslow));
    740:   }

    742:   if (tab->np == 3 && mprk->is_medium) {
    743:     computedstages = 0;
    744:     for (j = 0; j < s; j++) {
    745:       if (tab->rmb[j]) wsb[computedstages - tab->sbase + (tab->rmb[j] - 1) % tab->sbase] += h * tab->bsb[j];
    746:       else wsb[computedstages++] = h * tab->bsb[j];
    747:     }
    748:     PetscCall(VecGetSubVector(X, mprk->is_slowbuffer, &Xslowbuffer));
    749:     PetscCall(VecMAXPY(Xslowbuffer, computedstages, wsb, mprk->YdotRHS_slowbuffer));
    750:     PetscCall(VecRestoreSubVector(X, mprk->is_slowbuffer, &Xslowbuffer));
    751:   } else {
    752:     /* slow buffer region */
    753:     for (j = 0; j < s; j++) wsb[j] = h * tab->bsb[j];
    754:     PetscCall(VecGetSubVector(X, mprk->is_slowbuffer, &Xslowbuffer));
    755:     PetscCall(VecMAXPY(Xslowbuffer, s, wsb, mprk->YdotRHS_slowbuffer));
    756:     PetscCall(VecRestoreSubVector(X, mprk->is_slowbuffer, &Xslowbuffer));
    757:   }
    758:   if (tab->np == 3) {
    759:     Vec          Xmedium, Xmediumbuffer;
    760:     PetscScalar *wm = mprk->work_medium, *wmb = mprk->work_mediumbuffer;
    761:     /* medium region and slow buffer region */
    762:     if (mprk->is_medium) {
    763:       computedstages = 0;
    764:       for (j = 0; j < s; j++) {
    765:         if (tab->rmb[j]) wm[computedstages - tab->sbase + (tab->rmb[j] - 1) % tab->sbase] += h * tab->bmb[j];
    766:         else wm[computedstages++] = h * tab->bmb[j];
    767:       }
    768:       PetscCall(VecGetSubVector(X, mprk->is_medium, &Xmedium));
    769:       PetscCall(VecMAXPY(Xmedium, computedstages, wm, mprk->YdotRHS_medium));
    770:       PetscCall(VecRestoreSubVector(X, mprk->is_medium, &Xmedium));
    771:     }
    772:     /* medium buffer region */
    773:     for (j = 0; j < s; j++) wmb[j] = h * tab->bmb[j];
    774:     PetscCall(VecGetSubVector(X, mprk->is_mediumbuffer, &Xmediumbuffer));
    775:     PetscCall(VecMAXPY(Xmediumbuffer, s, wmb, mprk->YdotRHS_mediumbuffer));
    776:     PetscCall(VecRestoreSubVector(X, mprk->is_mediumbuffer, &Xmediumbuffer));
    777:   }
    778:   /* fast region */
    779:   for (j = 0; j < s; j++) wf[j] = h * tab->bf[j];
    780:   PetscCall(VecGetSubVector(X, mprk->is_fast, &Xfast));
    781:   PetscCall(VecMAXPY(Xfast, s, wf, mprk->YdotRHS_fast));
    782:   PetscCall(VecRestoreSubVector(X, mprk->is_fast, &Xfast));
    783:   PetscFunctionReturn(PETSC_SUCCESS);
    784: }

    786: static PetscErrorCode TSStep_MPRKSPLIT(TS ts)
    787: {
    788:   TS_MPRK         *mprk = (TS_MPRK *)ts->data;
    789:   MPRKTableau      tab  = mprk->tableau;
    790:   Vec             *Y = mprk->Y, *YdotRHS_fast = mprk->YdotRHS_fast, *YdotRHS_slow = mprk->YdotRHS_slow, *YdotRHS_slowbuffer = mprk->YdotRHS_slowbuffer;
    791:   Vec              Yslow, Yslowbuffer, Yfast; /* subvectors for slow and fast components in Y[i] respectively */
    792:   PetscInt         s  = tab->s;
    793:   const PetscReal *Af = tab->Af, *cf = tab->cf, *Asb = tab->Asb, *csb = tab->csb;
    794:   PetscScalar     *wf = mprk->work_fast, *ws = mprk->work_slow, *wsb = mprk->work_slowbuffer;
    795:   PetscInt         i, j, computedstages;
    796:   PetscReal        next_time_step = ts->time_step, t = ts->ptime, h = ts->time_step;

    798:   PetscFunctionBegin;
    799:   for (i = 0; i < s; i++) {
    800:     mprk->stage_time = t + h * cf[i];
    801:     PetscCall(TSPreStage(ts, mprk->stage_time));
    802:     /* calculate the stage value for fast and slow components respectively */
    803:     PetscCall(VecCopy(ts->vec_sol, Y[i]));
    804:     for (j = 0; j < i; j++) wsb[j] = h * Asb[i * s + j];

    806:     /* slow buffer region */
    807:     if (tab->np == 3 && mprk->is_medium) {
    808:       if (tab->rmb[i]) {
    809:         PetscCall(VecGetSubVector(Y[i], mprk->is_slowbuffer, &Yslowbuffer));
    810:         PetscCall(VecISCopy(Y[tab->rmb[i] - 1], mprk->is_slowbuffer, SCATTER_REVERSE, Yslowbuffer));
    811:         PetscCall(VecRestoreSubVector(Y[i], mprk->is_slowbuffer, &Yslowbuffer));
    812:       } else {
    813:         PetscScalar *wm = mprk->work_medium;
    814:         computedstages  = 0;
    815:         for (j = 0; j < i; j++) {
    816:           if (tab->rmb[j]) wm[computedstages - tab->sbase + (tab->rmb[j] - 1) % tab->sbase] += wsb[j];
    817:           else wm[computedstages++] = wsb[j];
    818:         }
    819:         PetscCall(VecGetSubVector(Y[i], mprk->is_slowbuffer, &Yslowbuffer));
    820:         PetscCall(VecMAXPY(Yslowbuffer, computedstages, wm, YdotRHS_slowbuffer));
    821:         PetscCall(VecRestoreSubVector(Y[i], mprk->is_slowbuffer, &Yslowbuffer));
    822:       }
    823:     } else {
    824:       PetscCall(VecGetSubVector(Y[i], mprk->is_slowbuffer, &Yslowbuffer));
    825:       PetscCall(VecMAXPY(Yslowbuffer, i, wsb, YdotRHS_slowbuffer));
    826:       PetscCall(VecRestoreSubVector(Y[i], mprk->is_slowbuffer, &Yslowbuffer));
    827:     }

    829:     /* slow region */
    830:     if (mprk->is_slow) {
    831:       if (tab->rsb[i]) { /* repeat previous stage */
    832:         PetscCall(VecGetSubVector(Y[i], mprk->is_slow, &Yslow));
    833:         PetscCall(VecISCopy(Y[tab->rsb[i] - 1], mprk->is_slow, SCATTER_REVERSE, Yslow));
    834:         PetscCall(VecRestoreSubVector(Y[i], mprk->is_slow, &Yslow));
    835:       } else {
    836:         computedstages = 0;
    837:         for (j = 0; j < i; j++) {
    838:           if (tab->rsb[j]) ws[tab->rsb[j] - 1] += wsb[j];
    839:           else ws[computedstages++] = wsb[j];
    840:         }
    841:         PetscCall(VecGetSubVector(Y[i], mprk->is_slow, &Yslow));
    842:         PetscCall(VecMAXPY(Yslow, computedstages, ws, YdotRHS_slow));
    843:         PetscCall(VecRestoreSubVector(Y[i], mprk->is_slow, &Yslow));
    844:         /* only depends on the slow buffer region */
    845:         PetscCall(TSComputeRHSFunction(mprk->subts_slow, t + h * csb[i], Y[i], YdotRHS_slow[computedstages]));
    846:       }
    847:     }

    849:     /* fast region */
    850:     for (j = 0; j < i; j++) wf[j] = h * Af[i * s + j];
    851:     PetscCall(VecGetSubVector(Y[i], mprk->is_fast, &Yfast));
    852:     PetscCall(VecMAXPY(Yfast, i, wf, YdotRHS_fast));
    853:     PetscCall(VecRestoreSubVector(Y[i], mprk->is_fast, &Yfast));

    855:     if (tab->np == 3) {
    856:       Vec             *YdotRHS_medium = mprk->YdotRHS_medium, *YdotRHS_mediumbuffer = mprk->YdotRHS_mediumbuffer;
    857:       Vec              Ymedium, Ymediumbuffer;
    858:       const PetscReal *Amb = tab->Amb, *cmb = tab->cmb;
    859:       PetscScalar     *wm = mprk->work_medium, *wmb = mprk->work_mediumbuffer;

    861:       for (j = 0; j < i; j++) wmb[j] = h * Amb[i * s + j];
    862:       /* medium buffer region */
    863:       PetscCall(VecGetSubVector(Y[i], mprk->is_mediumbuffer, &Ymediumbuffer));
    864:       PetscCall(VecMAXPY(Ymediumbuffer, i, wmb, YdotRHS_mediumbuffer));
    865:       PetscCall(VecRestoreSubVector(Y[i], mprk->is_mediumbuffer, &Ymediumbuffer));
    866:       /* medium region */
    867:       if (mprk->is_medium) {
    868:         if (tab->rmb[i]) { /* repeat previous stage */
    869:           PetscCall(VecGetSubVector(Y[i], mprk->is_medium, &Ymedium));
    870:           PetscCall(VecISCopy(Y[tab->rmb[i] - 1], mprk->is_medium, SCATTER_REVERSE, Ymedium));
    871:           PetscCall(VecRestoreSubVector(Y[i], mprk->is_medium, &Ymedium));
    872:         } else {
    873:           computedstages = 0;
    874:           for (j = 0; j < i; j++) {
    875:             if (tab->rmb[j]) wm[computedstages - tab->sbase + (tab->rmb[j] - 1) % tab->sbase] += wmb[j];
    876:             else wm[computedstages++] = wmb[j];
    877:           }
    878:           PetscCall(VecGetSubVector(Y[i], mprk->is_medium, &Ymedium));
    879:           PetscCall(VecMAXPY(Ymedium, computedstages, wm, YdotRHS_medium));
    880:           PetscCall(VecRestoreSubVector(Y[i], mprk->is_medium, &Ymedium));
    881:           /* only depends on the medium buffer region */
    882:           PetscCall(TSComputeRHSFunction(mprk->subts_medium, t + h * cmb[i], Y[i], YdotRHS_medium[computedstages]));
    883:           /* must be computed after all regions are updated in Y */
    884:           PetscCall(TSComputeRHSFunction(mprk->subts_slowbuffer, t + h * csb[i], Y[i], YdotRHS_slowbuffer[computedstages]));
    885:         }
    886:       }
    887:       /* must be computed after fast region and slow region are updated in Y */
    888:       PetscCall(TSComputeRHSFunction(mprk->subts_mediumbuffer, t + h * cmb[i], Y[i], YdotRHS_mediumbuffer[i]));
    889:     }
    890:     if (!(tab->np == 3 && mprk->is_medium)) PetscCall(TSComputeRHSFunction(mprk->subts_slowbuffer, t + h * csb[i], Y[i], YdotRHS_slowbuffer[i]));
    891:     PetscCall(TSComputeRHSFunction(mprk->subts_fast, t + h * cf[i], Y[i], YdotRHS_fast[i]));
    892:   }

    894:   PetscCall(TSEvaluateStep(ts, tab->order, ts->vec_sol, NULL));
    895:   ts->ptime += ts->time_step;
    896:   ts->time_step = next_time_step;
    897:   PetscFunctionReturn(PETSC_SUCCESS);
    898: }

    900: static PetscErrorCode TSMPRKTableauReset(TS ts)
    901: {
    902:   TS_MPRK    *mprk = (TS_MPRK *)ts->data;
    903:   MPRKTableau tab  = mprk->tableau;

    905:   PetscFunctionBegin;
    906:   if (!tab) PetscFunctionReturn(PETSC_SUCCESS);
    907:   PetscCall(PetscFree(mprk->work_fast));
    908:   PetscCall(PetscFree(mprk->work_slow));
    909:   PetscCall(PetscFree(mprk->work_slowbuffer));
    910:   PetscCall(PetscFree(mprk->work_medium));
    911:   PetscCall(PetscFree(mprk->work_mediumbuffer));
    912:   PetscCall(VecDestroyVecs(tab->s, &mprk->Y));
    913:   if (ts->use_splitrhsfunction) {
    914:     PetscCall(VecDestroyVecs(tab->s, &mprk->YdotRHS_fast));
    915:     PetscCall(VecDestroyVecs(tab->s, &mprk->YdotRHS_slow));
    916:     PetscCall(VecDestroyVecs(tab->s, &mprk->YdotRHS_slowbuffer));
    917:     PetscCall(VecDestroyVecs(tab->s, &mprk->YdotRHS_medium));
    918:     PetscCall(VecDestroyVecs(tab->s, &mprk->YdotRHS_mediumbuffer));
    919:   } else {
    920:     PetscCall(VecDestroyVecs(tab->s, &mprk->YdotRHS));
    921:     if (mprk->is_slow) PetscCall(PetscFree(mprk->YdotRHS_slow));
    922:     PetscCall(PetscFree(mprk->YdotRHS_slowbuffer));
    923:     if (tab->np == 3) {
    924:       if (mprk->is_medium) PetscCall(PetscFree(mprk->YdotRHS_medium));
    925:       PetscCall(PetscFree(mprk->YdotRHS_mediumbuffer));
    926:     }
    927:     PetscCall(PetscFree(mprk->YdotRHS_fast));
    928:   }
    929:   PetscFunctionReturn(PETSC_SUCCESS);
    930: }

    932: static PetscErrorCode TSReset_MPRK(TS ts)
    933: {
    934:   PetscFunctionBegin;
    935:   PetscCall(TSMPRKTableauReset(ts));
    936:   PetscFunctionReturn(PETSC_SUCCESS);
    937: }

    939: static PetscErrorCode DMCoarsenHook_TSMPRK(DM fine, DM coarse, void *ctx)
    940: {
    941:   PetscFunctionBegin;
    942:   PetscFunctionReturn(PETSC_SUCCESS);
    943: }

    945: static PetscErrorCode DMRestrictHook_TSMPRK(DM fine, Mat restrct, Vec rscale, Mat inject, DM coarse, void *ctx)
    946: {
    947:   PetscFunctionBegin;
    948:   PetscFunctionReturn(PETSC_SUCCESS);
    949: }

    951: static PetscErrorCode DMSubDomainHook_TSMPRK(DM dm, DM subdm, void *ctx)
    952: {
    953:   PetscFunctionBegin;
    954:   PetscFunctionReturn(PETSC_SUCCESS);
    955: }

    957: static PetscErrorCode DMSubDomainRestrictHook_TSMPRK(DM dm, VecScatter gscat, VecScatter lscat, DM subdm, void *ctx)
    958: {
    959:   PetscFunctionBegin;
    960:   PetscFunctionReturn(PETSC_SUCCESS);
    961: }

    963: static PetscErrorCode TSMPRKTableauSetUp(TS ts)
    964: {
    965:   TS_MPRK    *mprk = (TS_MPRK *)ts->data;
    966:   MPRKTableau tab  = mprk->tableau;
    967:   Vec         YdotRHS_slow, YdotRHS_slowbuffer, YdotRHS_medium, YdotRHS_mediumbuffer, YdotRHS_fast;

    969:   PetscFunctionBegin;
    970:   PetscCall(VecDuplicateVecs(ts->vec_sol, tab->s, &mprk->Y));
    971:   if (mprk->is_slow) PetscCall(PetscMalloc1(tab->s, &mprk->work_slow));
    972:   PetscCall(PetscMalloc1(tab->s, &mprk->work_slowbuffer));
    973:   if (tab->np == 3) {
    974:     if (mprk->is_medium) PetscCall(PetscMalloc1(tab->s, &mprk->work_medium));
    975:     PetscCall(PetscMalloc1(tab->s, &mprk->work_mediumbuffer));
    976:   }
    977:   PetscCall(PetscMalloc1(tab->s, &mprk->work_fast));

    979:   if (ts->use_splitrhsfunction) {
    980:     if (mprk->is_slow) {
    981:       PetscCall(VecGetSubVector(ts->vec_sol, mprk->is_slow, &YdotRHS_slow));
    982:       PetscCall(VecDuplicateVecs(YdotRHS_slow, tab->s, &mprk->YdotRHS_slow));
    983:       PetscCall(VecRestoreSubVector(ts->vec_sol, mprk->is_slow, &YdotRHS_slow));
    984:     }
    985:     PetscCall(VecGetSubVector(ts->vec_sol, mprk->is_slowbuffer, &YdotRHS_slowbuffer));
    986:     PetscCall(VecDuplicateVecs(YdotRHS_slowbuffer, tab->s, &mprk->YdotRHS_slowbuffer));
    987:     PetscCall(VecRestoreSubVector(ts->vec_sol, mprk->is_slowbuffer, &YdotRHS_slowbuffer));
    988:     if (tab->np == 3) {
    989:       if (mprk->is_medium) {
    990:         PetscCall(VecGetSubVector(ts->vec_sol, mprk->is_medium, &YdotRHS_medium));
    991:         PetscCall(VecDuplicateVecs(YdotRHS_medium, tab->s, &mprk->YdotRHS_medium));
    992:         PetscCall(VecRestoreSubVector(ts->vec_sol, mprk->is_medium, &YdotRHS_medium));
    993:       }
    994:       PetscCall(VecGetSubVector(ts->vec_sol, mprk->is_mediumbuffer, &YdotRHS_mediumbuffer));
    995:       PetscCall(VecDuplicateVecs(YdotRHS_mediumbuffer, tab->s, &mprk->YdotRHS_mediumbuffer));
    996:       PetscCall(VecRestoreSubVector(ts->vec_sol, mprk->is_mediumbuffer, &YdotRHS_mediumbuffer));
    997:     }
    998:     PetscCall(VecGetSubVector(ts->vec_sol, mprk->is_fast, &YdotRHS_fast));
    999:     PetscCall(VecDuplicateVecs(YdotRHS_fast, tab->s, &mprk->YdotRHS_fast));
    1000:     PetscCall(VecRestoreSubVector(ts->vec_sol, mprk->is_fast, &YdotRHS_fast));
    1001:   } else {
    1002:     PetscCall(VecDuplicateVecs(ts->vec_sol, tab->s, &mprk->YdotRHS));
    1003:     if (mprk->is_slow) PetscCall(PetscMalloc1(tab->s, &mprk->YdotRHS_slow));
    1004:     PetscCall(PetscMalloc1(tab->s, &mprk->YdotRHS_slowbuffer));
    1005:     if (tab->np == 3) {
    1006:       if (mprk->is_medium) PetscCall(PetscMalloc1(tab->s, &mprk->YdotRHS_medium));
    1007:       PetscCall(PetscMalloc1(tab->s, &mprk->YdotRHS_mediumbuffer));
    1008:     }
    1009:     PetscCall(PetscMalloc1(tab->s, &mprk->YdotRHS_fast));
    1010:   }
    1011:   PetscFunctionReturn(PETSC_SUCCESS);
    1012: }

    1014: static PetscErrorCode TSSetUp_MPRK(TS ts)
    1015: {
    1016:   TS_MPRK    *mprk = (TS_MPRK *)ts->data;
    1017:   MPRKTableau tab  = mprk->tableau;
    1018:   DM          dm;

    1020:   PetscFunctionBegin;
    1021:   PetscCall(TSRHSSplitGetIS(ts, "slow", &mprk->is_slow));
    1022:   PetscCall(TSRHSSplitGetIS(ts, "fast", &mprk->is_fast));
    1023:   PetscCheck(mprk->is_slow && mprk->is_fast, PetscObjectComm((PetscObject)ts), PETSC_ERR_USER, "Must set up RHSSplits with TSRHSSplitSetIS() using split names 'slow' and 'fast' respectively in order to use the method '%s'", tab->name);

    1025:   if (tab->np == 3) {
    1026:     PetscCall(TSRHSSplitGetIS(ts, "medium", &mprk->is_medium));
    1027:     PetscCheck(mprk->is_medium, PetscObjectComm((PetscObject)ts), PETSC_ERR_USER, "Must set up RHSSplits with TSRHSSplitSetIS() using split names 'slow' and 'medium' and 'fast' respectively in order to use the method '%s'", tab->name);
    1028:     PetscCall(TSRHSSplitGetIS(ts, "mediumbuffer", &mprk->is_mediumbuffer));
    1029:     if (!mprk->is_mediumbuffer) { /* let medium buffer cover whole medium region */
    1030:       mprk->is_mediumbuffer = mprk->is_medium;
    1031:       mprk->is_medium       = NULL;
    1032:     }
    1033:   }

    1035:   /* If users do not provide buffer region settings, the solver will do them automatically, but with a performance penalty */
    1036:   PetscCall(TSRHSSplitGetIS(ts, "slowbuffer", &mprk->is_slowbuffer));
    1037:   if (!mprk->is_slowbuffer) { /* let slow buffer cover whole slow region */
    1038:     mprk->is_slowbuffer = mprk->is_slow;
    1039:     mprk->is_slow       = NULL;
    1040:   }
    1041:   PetscCall(TSCheckImplicitTerm(ts));
    1042:   PetscCall(TSMPRKTableauSetUp(ts));
    1043:   PetscCall(TSGetDM(ts, &dm));
    1044:   PetscCall(DMCoarsenHookAdd(dm, DMCoarsenHook_TSMPRK, DMRestrictHook_TSMPRK, ts));
    1045:   PetscCall(DMSubDomainHookAdd(dm, DMSubDomainHook_TSMPRK, DMSubDomainRestrictHook_TSMPRK, ts));
    1046:   if (ts->use_splitrhsfunction) {
    1047:     ts->ops->step         = TSStep_MPRKSPLIT;
    1048:     ts->ops->evaluatestep = TSEvaluateStep_MPRKSPLIT;
    1049:     PetscCall(TSMPRKSetSplits(ts));
    1050:   } else {
    1051:     ts->ops->step         = TSStep_MPRK;
    1052:     ts->ops->evaluatestep = TSEvaluateStep_MPRK;
    1053:   }
    1054:   PetscFunctionReturn(PETSC_SUCCESS);
    1055: }

    1057: static PetscErrorCode TSSetFromOptions_MPRK(TS ts, PetscOptionItems PetscOptionsObject)
    1058: {
    1059:   TS_MPRK *mprk = (TS_MPRK *)ts->data;

    1061:   PetscFunctionBegin;
    1062:   PetscOptionsHeadBegin(PetscOptionsObject, "PRK ODE solver options");
    1063:   {
    1064:     MPRKTableauLink link;
    1065:     PetscInt        count, choice;
    1066:     PetscBool       flg;
    1067:     const char    **namelist;
    1068:     for (link = MPRKTableauList, count = 0; link; link = link->next, count++);
    1069:     PetscCall(PetscMalloc1(count, (char ***)&namelist));
    1070:     for (link = MPRKTableauList, count = 0; link; link = link->next, count++) namelist[count] = link->tab.name;
    1071:     PetscCall(PetscOptionsEList("-ts_mprk_type", "Family of MPRK method", "TSMPRKSetType", (const char *const *)namelist, count, mprk->tableau->name, &choice, &flg));
    1072:     if (flg) PetscCall(TSMPRKSetType(ts, namelist[choice]));
    1073:     PetscCall(PetscFree(namelist));
    1074:   }
    1075:   PetscOptionsHeadEnd();
    1076:   PetscFunctionReturn(PETSC_SUCCESS);
    1077: }

    1079: static PetscErrorCode TSView_MPRK(TS ts, PetscViewer viewer)
    1080: {
    1081:   TS_MPRK  *mprk = (TS_MPRK *)ts->data;
    1082:   PetscBool isascii;

    1084:   PetscFunctionBegin;
    1085:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
    1086:   if (isascii) {
    1087:     MPRKTableau tab = mprk->tableau;
    1088:     TSMPRKType  mprktype;
    1089:     char        fbuf[512];
    1090:     char        sbuf[512];
    1091:     PetscInt    i;
    1092:     PetscCall(TSMPRKGetType(ts, &mprktype));
    1093:     PetscCall(PetscViewerASCIIPrintf(viewer, "  MPRK type %s\n", mprktype));
    1094:     PetscCall(PetscViewerASCIIPrintf(viewer, "  Order: %" PetscInt_FMT "\n", tab->order));

    1096:     PetscCall(PetscFormatRealArray(fbuf, sizeof(fbuf), "% 8.6f", tab->s, tab->cf));
    1097:     PetscCall(PetscViewerASCIIPrintf(viewer, "  Abscissa cf = %s\n", fbuf));
    1098:     PetscCall(PetscViewerASCIIPrintf(viewer, "  Af = \n"));
    1099:     for (i = 0; i < tab->s; i++) {
    1100:       PetscCall(PetscFormatRealArray(fbuf, sizeof(fbuf), "% 8.6f", tab->s, &tab->Af[i * tab->s]));
    1101:       PetscCall(PetscViewerASCIIPrintf(viewer, "    %s\n", fbuf));
    1102:     }
    1103:     PetscCall(PetscFormatRealArray(fbuf, sizeof(fbuf), "% 8.6f", tab->s, tab->bf));
    1104:     PetscCall(PetscViewerASCIIPrintf(viewer, "  bf = %s\n", fbuf));

    1106:     PetscCall(PetscFormatRealArray(sbuf, sizeof(sbuf), "% 8.6f", tab->s, tab->csb));
    1107:     PetscCall(PetscViewerASCIIPrintf(viewer, "  Abscissa csb = %s\n", sbuf));
    1108:     PetscCall(PetscViewerASCIIPrintf(viewer, "  Asb = \n"));
    1109:     for (i = 0; i < tab->s; i++) {
    1110:       PetscCall(PetscFormatRealArray(sbuf, sizeof(sbuf), "% 8.6f", tab->s, &tab->Asb[i * tab->s]));
    1111:       PetscCall(PetscViewerASCIIPrintf(viewer, "    %s\n", sbuf));
    1112:     }
    1113:     PetscCall(PetscFormatRealArray(sbuf, sizeof(sbuf), "% 8.6f", tab->s, tab->bsb));
    1114:     PetscCall(PetscViewerASCIIPrintf(viewer, "  bsb = %s\n", sbuf));

    1116:     if (tab->np == 3) {
    1117:       char mbuf[512];
    1118:       PetscCall(PetscFormatRealArray(mbuf, sizeof(mbuf), "% 8.6f", tab->s, tab->cmb));
    1119:       PetscCall(PetscViewerASCIIPrintf(viewer, "  Abscissa cmb = %s\n", mbuf));
    1120:       PetscCall(PetscViewerASCIIPrintf(viewer, "  Amb = \n"));
    1121:       for (i = 0; i < tab->s; i++) {
    1122:         PetscCall(PetscFormatRealArray(mbuf, sizeof(mbuf), "% 8.6f", tab->s, &tab->Amb[i * tab->s]));
    1123:         PetscCall(PetscViewerASCIIPrintf(viewer, "    %s\n", mbuf));
    1124:       }
    1125:       PetscCall(PetscFormatRealArray(mbuf, sizeof(mbuf), "% 8.6f", tab->s, tab->bmb));
    1126:       PetscCall(PetscViewerASCIIPrintf(viewer, "  bmb = %s\n", mbuf));
    1127:     }
    1128:   }
    1129:   PetscFunctionReturn(PETSC_SUCCESS);
    1130: }

    1132: static PetscErrorCode TSLoad_MPRK(TS ts, PetscViewer viewer)
    1133: {
    1134:   TSAdapt adapt;

    1136:   PetscFunctionBegin;
    1137:   PetscCall(TSGetAdapt(ts, &adapt));
    1138:   PetscCall(TSAdaptLoad(adapt, viewer));
    1139:   PetscFunctionReturn(PETSC_SUCCESS);
    1140: }

    1142: /*@
    1143:   TSMPRKSetType - Set the type of `TSMPRK` scheme

    1145:   Not Collective

    1147:   Input Parameters:
    1148: + ts       - timestepping context
    1149: - mprktype - type of `TSMPRK` scheme

    1151:   Options Database Key:
    1152: . -ts_mprk_type - <pm2,p2,p3> - select the specific scheme

    1154:   Level: intermediate

    1156: .seealso: [](ch_ts), `TSMPRKGetType()`, `TSMPRK`, `TSMPRKType`
    1157: @*/
    1158: PetscErrorCode TSMPRKSetType(TS ts, TSMPRKType mprktype)
    1159: {
    1160:   PetscFunctionBegin;
    1162:   PetscAssertPointer(mprktype, 2);
    1163:   PetscTryMethod(ts, "TSMPRKSetType_C", (TS, TSMPRKType), (ts, mprktype));
    1164:   PetscFunctionReturn(PETSC_SUCCESS);
    1165: }

    1167: /*@
    1168:   TSMPRKGetType - Get the type of `TSMPRK` scheme

    1170:   Not Collective

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

    1175:   Output Parameter:
    1176: . mprktype - type of `TSMPRK` scheme

    1178:   Level: intermediate

    1180: .seealso: [](ch_ts), `TSMPRK`
    1181: @*/
    1182: PetscErrorCode TSMPRKGetType(TS ts, TSMPRKType *mprktype)
    1183: {
    1184:   PetscFunctionBegin;
    1186:   PetscUseMethod(ts, "TSMPRKGetType_C", (TS, TSMPRKType *), (ts, mprktype));
    1187:   PetscFunctionReturn(PETSC_SUCCESS);
    1188: }

    1190: static PetscErrorCode TSMPRKGetType_MPRK(TS ts, TSMPRKType *mprktype)
    1191: {
    1192:   TS_MPRK *mprk = (TS_MPRK *)ts->data;

    1194:   PetscFunctionBegin;
    1195:   *mprktype = mprk->tableau->name;
    1196:   PetscFunctionReturn(PETSC_SUCCESS);
    1197: }

    1199: static PetscErrorCode TSMPRKSetType_MPRK(TS ts, TSMPRKType mprktype)
    1200: {
    1201:   TS_MPRK        *mprk = (TS_MPRK *)ts->data;
    1202:   PetscBool       match;
    1203:   MPRKTableauLink link;

    1205:   PetscFunctionBegin;
    1206:   if (mprk->tableau) {
    1207:     PetscCall(PetscStrcmp(mprk->tableau->name, mprktype, &match));
    1208:     if (match) PetscFunctionReturn(PETSC_SUCCESS);
    1209:   }
    1210:   for (link = MPRKTableauList; link; link = link->next) {
    1211:     PetscCall(PetscStrcmp(link->tab.name, mprktype, &match));
    1212:     if (match) {
    1213:       if (ts->setupcalled) PetscCall(TSMPRKTableauReset(ts));
    1214:       mprk->tableau = &link->tab;
    1215:       if (ts->setupcalled) PetscCall(TSMPRKTableauSetUp(ts));
    1216:       PetscFunctionReturn(PETSC_SUCCESS);
    1217:     }
    1218:   }
    1219:   SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_UNKNOWN_TYPE, "Could not find '%s'", mprktype);
    1220: }

    1222: static PetscErrorCode TSGetStages_MPRK(TS ts, PetscInt *ns, Vec **Y)
    1223: {
    1224:   TS_MPRK *mprk = (TS_MPRK *)ts->data;

    1226:   PetscFunctionBegin;
    1227:   *ns = mprk->tableau->s;
    1228:   if (Y) *Y = mprk->Y;
    1229:   PetscFunctionReturn(PETSC_SUCCESS);
    1230: }

    1232: static PetscErrorCode TSDestroy_MPRK(TS ts)
    1233: {
    1234:   PetscFunctionBegin;
    1235:   PetscCall(TSReset_MPRK(ts));
    1236:   if (ts->dm) {
    1237:     PetscCall(DMCoarsenHookRemove(ts->dm, DMCoarsenHook_TSMPRK, DMRestrictHook_TSMPRK, ts));
    1238:     PetscCall(DMSubDomainHookRemove(ts->dm, DMSubDomainHook_TSMPRK, DMSubDomainRestrictHook_TSMPRK, ts));
    1239:   }
    1240:   PetscCall(PetscFree(ts->data));
    1241:   PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSMPRKGetType_C", NULL));
    1242:   PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSMPRKSetType_C", NULL));
    1243:   PetscFunctionReturn(PETSC_SUCCESS);
    1244: }

    1246: /*MC
    1247:       TSMPRK - ODE solver using Multirate Partitioned Runge-Kutta schemes

    1249:   The user should provide the right-hand side of the equation using `TSSetRHSFunction()`.

    1251:   Level: beginner

    1253:   Note:
    1254:   The default is `TSMPRKPM2`, it can be changed with `TSMPRKSetType()` or -ts_mprk_type

    1256: .seealso: [](ch_ts), `TSCreate()`, `TS`, `TSSetType()`, `TSMPRKSetType()`, `TSMPRKGetType()`, `TSMPRKType`, `TSMPRKRegister()`, `TSMPRKSetMultirateType()`
    1257:           `TSMPRKM2`, `TSMPRKM3`, `TSMPRKRFSMR3`, `TSMPRKRFSMR2`, `TSType`
    1258: M*/
    1259: PETSC_EXTERN PetscErrorCode TSCreate_MPRK(TS ts)
    1260: {
    1261:   TS_MPRK *mprk;

    1263:   PetscFunctionBegin;
    1264:   PetscCall(TSMPRKInitializePackage());

    1266:   ts->ops->reset          = TSReset_MPRK;
    1267:   ts->ops->destroy        = TSDestroy_MPRK;
    1268:   ts->ops->view           = TSView_MPRK;
    1269:   ts->ops->load           = TSLoad_MPRK;
    1270:   ts->ops->setup          = TSSetUp_MPRK;
    1271:   ts->ops->step           = TSStep_MPRK;
    1272:   ts->ops->evaluatestep   = TSEvaluateStep_MPRK;
    1273:   ts->ops->setfromoptions = TSSetFromOptions_MPRK;
    1274:   ts->ops->getstages      = TSGetStages_MPRK;

    1276:   PetscCall(PetscNew(&mprk));
    1277:   ts->data = (void *)mprk;

    1279:   PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSMPRKGetType_C", TSMPRKGetType_MPRK));
    1280:   PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSMPRKSetType_C", TSMPRKSetType_MPRK));

    1282:   PetscCall(TSMPRKSetType(ts, TSMPRKDefault));
    1283:   PetscFunctionReturn(PETSC_SUCCESS);
    1284: }