Actual source code: mrk.c

      1: /*
      2:   Code for time stepping with the multi-rate Runge-Kutta method

      4:   Notes:
      5:   1) The general system is written as
      6:      Udot = F(t,U) for the nonsplit version of multi-rate RK,
      7:      user should give 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) for multi-rate RK with RHS splits,
     11:      user should partition RHS by themselves and also provide the indexes for both slow and fast components.
     12: */

     14: #include <petsc/private/tsimpl.h>
     15: #include <petscdm.h>
     16: #include <../src/ts/impls/explicit/rk/rk.h>
     17: #include <../src/ts/impls/explicit/rk/mrk.h>

     19: static PetscErrorCode TSReset_RK_MultirateNonsplit(TS ts)
     20: {
     21:   TS_RK    *rk  = (TS_RK *)ts->data;
     22:   RKTableau tab = rk->tableau;

     24:   PetscFunctionBegin;
     25:   PetscCall(VecDestroy(&rk->X0));
     26:   PetscCall(VecDestroyVecs(tab->s, &rk->YdotRHS_slow));
     27:   PetscFunctionReturn(PETSC_SUCCESS);
     28: }

     30: static PetscErrorCode TSInterpolate_RK_MultirateNonsplit(TS ts, PetscReal itime, Vec X)
     31: {
     32:   TS_RK           *rk = (TS_RK *)ts->data;
     33:   PetscInt         s = rk->tableau->s, p = rk->tableau->p, i, j;
     34:   PetscReal        h = ts->time_step;
     35:   PetscReal        tt, t;
     36:   PetscScalar     *b;
     37:   const PetscReal *B = rk->tableau->binterp;

     39:   PetscFunctionBegin;
     40:   PetscCheck(B, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "TSRK %s does not have an interpolation formula", rk->tableau->name);
     41:   t = (itime - rk->ptime) / h;
     42:   PetscCall(PetscMalloc1(s, &b));
     43:   for (i = 0; i < s; i++) b[i] = 0;
     44:   for (j = 0, tt = t; j < p; j++, tt *= t) {
     45:     for (i = 0; i < s; i++) b[i] += h * B[i * p + j] * tt;
     46:   }
     47:   PetscCall(VecCopy(rk->X0, X));
     48:   PetscCall(VecMAXPY(X, s, b, rk->YdotRHS_slow));
     49:   PetscCall(PetscFree(b));
     50:   PetscFunctionReturn(PETSC_SUCCESS);
     51: }

     53: static PetscErrorCode TSStepRefine_RK_MultirateNonsplit(TS ts)
     54: {
     55:   TS               previousts, subts;
     56:   TS_RK           *rk  = (TS_RK *)ts->data;
     57:   RKTableau        tab = rk->tableau;
     58:   Vec             *Y = rk->Y, *YdotRHS = rk->YdotRHS;
     59:   Vec              vec_fast, subvec_fast;
     60:   const PetscInt   s = tab->s;
     61:   const PetscReal *A = tab->A, *c = tab->c;
     62:   PetscScalar     *w = rk->work;
     63:   PetscInt         i, j, k;
     64:   PetscReal        t = ts->ptime, h = ts->time_step;

     66:   PetscFunctionBegin;
     67:   PetscCall(VecDuplicate(ts->vec_sol, &vec_fast));
     68:   previousts = rk->subts_current;
     69:   PetscCall(TSRHSSplitGetSubTS(rk->subts_current, "fast", &subts));
     70:   PetscCall(TSRHSSplitGetSubTS(subts, "fast", &subts));
     71:   for (k = 0; k < rk->dtratio; k++) {
     72:     for (i = 0; i < s; i++) {
     73:       PetscCall(TSInterpolate_RK_MultirateNonsplit(ts, t + k * h / rk->dtratio + h / rk->dtratio * c[i], Y[i]));
     74:       for (j = 0; j < i; j++) w[j] = h / rk->dtratio * A[i * s + j];
     75:       /* update the fast components in the stage value, the slow components will be ignored, so we do not care the slow part in vec_fast */
     76:       PetscCall(VecCopy(ts->vec_sol, vec_fast));
     77:       PetscCall(VecMAXPY(vec_fast, i, w, YdotRHS));
     78:       /* update the fast components in the stage value */
     79:       PetscCall(VecGetSubVector(vec_fast, rk->is_fast, &subvec_fast));
     80:       PetscCall(VecISCopy(Y[i], rk->is_fast, SCATTER_FORWARD, subvec_fast));
     81:       PetscCall(VecRestoreSubVector(vec_fast, rk->is_fast, &subvec_fast));
     82:       /* compute the stage RHS */
     83:       PetscCall(TSComputeRHSFunction(ts, t + k * h / rk->dtratio + h / rk->dtratio * c[i], Y[i], YdotRHS[i]));
     84:     }
     85:     PetscCall(VecCopy(ts->vec_sol, vec_fast));
     86:     PetscCall(TSEvaluateStep(ts, tab->order, vec_fast, NULL));
     87:     /* update the fast components in the solution */
     88:     PetscCall(VecGetSubVector(vec_fast, rk->is_fast, &subvec_fast));
     89:     PetscCall(VecISCopy(ts->vec_sol, rk->is_fast, SCATTER_FORWARD, subvec_fast));
     90:     PetscCall(VecRestoreSubVector(vec_fast, rk->is_fast, &subvec_fast));

     92:     if (subts) {
     93:       Vec *YdotRHS_copy;
     94:       PetscCall(VecDuplicateVecs(ts->vec_sol, s, &YdotRHS_copy));
     95:       rk->subts_current = rk->subts_fast;
     96:       ts->ptime         = t + k * h / rk->dtratio;
     97:       ts->time_step     = h / rk->dtratio;
     98:       PetscCall(TSRHSSplitGetIS(rk->subts_current, "fast", &rk->is_fast));
     99:       for (i = 0; i < s; i++) {
    100:         PetscCall(VecCopy(rk->YdotRHS_slow[i], YdotRHS_copy[i]));
    101:         PetscCall(VecCopy(YdotRHS[i], rk->YdotRHS_slow[i]));
    102:       }

    104:       PetscCall(TSStepRefine_RK_MultirateNonsplit(ts));

    106:       rk->subts_current = previousts;
    107:       ts->ptime         = t;
    108:       ts->time_step     = h;
    109:       PetscCall(TSRHSSplitGetIS(previousts, "fast", &rk->is_fast));
    110:       for (i = 0; i < s; i++) PetscCall(VecCopy(YdotRHS_copy[i], rk->YdotRHS_slow[i]));
    111:       PetscCall(VecDestroyVecs(s, &YdotRHS_copy));
    112:     }
    113:   }
    114:   PetscCall(VecDestroy(&vec_fast));
    115:   PetscFunctionReturn(PETSC_SUCCESS);
    116: }

    118: static PetscErrorCode TSStep_RK_MultirateNonsplit(TS ts)
    119: {
    120:   TS_RK           *rk  = (TS_RK *)ts->data;
    121:   RKTableau        tab = rk->tableau;
    122:   Vec             *Y = rk->Y, *YdotRHS = rk->YdotRHS, *YdotRHS_slow = rk->YdotRHS_slow;
    123:   Vec              stage_slow, sol_slow; /* vectors store the slow components */
    124:   Vec              subvec_slow;          /* sub vector to store the slow components */
    125:   IS               is_slow = rk->is_slow;
    126:   const PetscInt   s       = tab->s;
    127:   const PetscReal *A = tab->A, *c = tab->c;
    128:   PetscScalar     *w = rk->work;
    129:   PetscInt         i, j, dtratio = rk->dtratio;
    130:   PetscReal        next_time_step = ts->time_step, t = ts->ptime, h = ts->time_step;

    132:   PetscFunctionBegin;
    133:   rk->status = TS_STEP_INCOMPLETE;
    134:   PetscCall(VecDuplicate(ts->vec_sol, &stage_slow));
    135:   PetscCall(VecDuplicate(ts->vec_sol, &sol_slow));
    136:   PetscCall(VecCopy(ts->vec_sol, rk->X0));
    137:   for (i = 0; i < s; i++) {
    138:     rk->stage_time = t + h * c[i];
    139:     PetscCall(TSPreStage(ts, rk->stage_time));
    140:     PetscCall(VecCopy(ts->vec_sol, Y[i]));
    141:     for (j = 0; j < i; j++) w[j] = h * A[i * s + j];
    142:     PetscCall(VecMAXPY(Y[i], i, w, YdotRHS_slow));
    143:     PetscCall(TSPostStage(ts, rk->stage_time, i, Y));
    144:     /* compute the stage RHS */
    145:     PetscCall(TSComputeRHSFunction(ts, t + h * c[i], Y[i], YdotRHS_slow[i]));
    146:   }
    147:   /* update the slow components in the solution */
    148:   rk->YdotRHS = YdotRHS_slow;
    149:   rk->dtratio = 1;
    150:   PetscCall(TSEvaluateStep(ts, tab->order, sol_slow, NULL));
    151:   rk->dtratio = dtratio;
    152:   rk->YdotRHS = YdotRHS;
    153:   /* update the slow components in the solution */
    154:   PetscCall(VecGetSubVector(sol_slow, is_slow, &subvec_slow));
    155:   PetscCall(VecISCopy(ts->vec_sol, is_slow, SCATTER_FORWARD, subvec_slow));
    156:   PetscCall(VecRestoreSubVector(sol_slow, is_slow, &subvec_slow));

    158:   rk->subts_current = ts;
    159:   rk->ptime         = t;
    160:   rk->time_step     = h;
    161:   PetscCall(TSStepRefine_RK_MultirateNonsplit(ts));

    163:   ts->ptime     = t + ts->time_step;
    164:   ts->time_step = next_time_step;
    165:   rk->status    = TS_STEP_COMPLETE;

    167:   /* free memory of work vectors */
    168:   PetscCall(VecDestroy(&stage_slow));
    169:   PetscCall(VecDestroy(&sol_slow));
    170:   PetscFunctionReturn(PETSC_SUCCESS);
    171: }

    173: static PetscErrorCode TSSetUp_RK_MultirateNonsplit(TS ts)
    174: {
    175:   TS_RK    *rk  = (TS_RK *)ts->data;
    176:   RKTableau tab = rk->tableau;

    178:   PetscFunctionBegin;
    179:   PetscCall(TSRHSSplitGetIS(ts, "slow", &rk->is_slow));
    180:   PetscCall(TSRHSSplitGetIS(ts, "fast", &rk->is_fast));
    181:   PetscCheck(rk->is_slow && rk->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 multirate RK");
    182:   PetscCall(TSRHSSplitGetSubTS(ts, "slow", &rk->subts_slow));
    183:   PetscCall(TSRHSSplitGetSubTS(ts, "fast", &rk->subts_fast));
    184:   PetscCheck(rk->subts_slow && rk->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");
    185:   PetscCall(VecDuplicate(ts->vec_sol, &rk->X0));
    186:   PetscCall(VecDuplicateVecs(ts->vec_sol, tab->s, &rk->YdotRHS_slow));
    187:   rk->subts_current = rk->subts_fast;

    189:   ts->ops->step        = TSStep_RK_MultirateNonsplit;
    190:   ts->ops->interpolate = TSInterpolate_RK_MultirateNonsplit;
    191:   PetscFunctionReturn(PETSC_SUCCESS);
    192: }

    194: /*
    195:   Copy DM from tssrc to tsdest, while keeping the original DMTS and DMSNES in tsdest.
    196: */
    197: static PetscErrorCode TSCopyDM(TS tssrc, TS tsdest)
    198: {
    199:   DM newdm, dmsrc, dmdest;

    201:   PetscFunctionBegin;
    202:   PetscCall(TSGetDM(tssrc, &dmsrc));
    203:   PetscCall(DMClone(dmsrc, &newdm));
    204:   PetscCall(TSGetDM(tsdest, &dmdest));
    205:   PetscCall(DMCopyDMTS(dmdest, newdm));
    206:   PetscCall(DMCopyDMSNES(dmdest, newdm));
    207:   PetscCall(TSSetDM(tsdest, newdm));
    208:   PetscCall(DMDestroy(&newdm));
    209:   PetscFunctionReturn(PETSC_SUCCESS);
    210: }

    212: static PetscErrorCode TSReset_RK_MultirateSplit(TS ts)
    213: {
    214:   TS_RK *rk = (TS_RK *)ts->data;

    216:   PetscFunctionBegin;
    217:   if (rk->subts_slow) {
    218:     PetscCall(PetscFree(rk->subts_slow->data));
    219:     rk->subts_slow = NULL;
    220:   }
    221:   if (rk->subts_fast) {
    222:     PetscCall(PetscFree(rk->YdotRHS_fast));
    223:     PetscCall(PetscFree(rk->YdotRHS_slow));
    224:     PetscCall(VecDestroy(&rk->X0));
    225:     PetscCall(TSReset_RK_MultirateSplit(rk->subts_fast));
    226:     PetscCall(PetscFree(rk->subts_fast->data));
    227:     rk->subts_fast = NULL;
    228:   }
    229:   PetscFunctionReturn(PETSC_SUCCESS);
    230: }

    232: static PetscErrorCode TSInterpolate_RK_MultirateSplit(TS ts, PetscReal itime, Vec X)
    233: {
    234:   TS_RK           *rk = (TS_RK *)ts->data;
    235:   Vec              Xslow;
    236:   PetscInt         s = rk->tableau->s, p = rk->tableau->p, i, j;
    237:   PetscReal        h;
    238:   PetscReal        tt, t;
    239:   PetscScalar     *b;
    240:   const PetscReal *B = rk->tableau->binterp;

    242:   PetscFunctionBegin;
    243:   PetscCheck(B, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "TSRK %s does not have an interpolation formula", rk->tableau->name);

    245:   switch (rk->status) {
    246:   case TS_STEP_INCOMPLETE:
    247:   case TS_STEP_PENDING:
    248:     h = ts->time_step;
    249:     t = (itime - ts->ptime) / h;
    250:     break;
    251:   case TS_STEP_COMPLETE:
    252:     h = ts->ptime - ts->ptime_prev;
    253:     t = (itime - ts->ptime) / h + 1; /* In the interval [0,1] */
    254:     break;
    255:   default:
    256:     SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_PLIB, "Invalid TSStepStatus");
    257:   }
    258:   PetscCall(PetscMalloc1(s, &b));
    259:   for (i = 0; i < s; i++) b[i] = 0;
    260:   for (j = 0, tt = t; j < p; j++, tt *= t) {
    261:     for (i = 0; i < s; i++) b[i] += h * B[i * p + j] * tt;
    262:   }
    263:   for (i = 0; i < s; i++) PetscCall(VecGetSubVector(rk->YdotRHS[i], rk->is_slow, &rk->YdotRHS_slow[i]));
    264:   PetscCall(VecGetSubVector(X, rk->is_slow, &Xslow));
    265:   PetscCall(VecISCopy(rk->X0, rk->is_slow, SCATTER_REVERSE, Xslow));
    266:   PetscCall(VecMAXPY(Xslow, s, b, rk->YdotRHS_slow));
    267:   PetscCall(VecRestoreSubVector(X, rk->is_slow, &Xslow));
    268:   for (i = 0; i < s; i++) PetscCall(VecRestoreSubVector(rk->YdotRHS[i], rk->is_slow, &rk->YdotRHS_slow[i]));
    269:   PetscCall(PetscFree(b));
    270:   PetscFunctionReturn(PETSC_SUCCESS);
    271: }

    273: /*
    274:  This is for partitioned RHS multirate RK method
    275:  The step completion formula is

    277:  x1 = x0 + h b^T YdotRHS

    279: */
    280: static PetscErrorCode TSEvaluateStep_RK_MultirateSplit(TS ts, PetscInt order, Vec X, PetscBool *done)
    281: {
    282:   TS_RK       *rk  = (TS_RK *)ts->data;
    283:   RKTableau    tab = rk->tableau;
    284:   Vec          Xslow, Xfast; /* subvectors of X which store slow components and fast components respectively */
    285:   PetscScalar *w = rk->work;
    286:   PetscReal    h = ts->time_step;
    287:   PetscInt     s = tab->s, j;

    289:   PetscFunctionBegin;
    290:   PetscCall(VecCopy(ts->vec_sol, X));
    291:   if (rk->slow) {
    292:     for (j = 0; j < s; j++) w[j] = h * tab->b[j];
    293:     PetscCall(VecGetSubVector(ts->vec_sol, rk->is_slow, &Xslow));
    294:     PetscCall(VecMAXPY(Xslow, s, w, rk->YdotRHS_slow));
    295:     PetscCall(VecRestoreSubVector(ts->vec_sol, rk->is_slow, &Xslow));
    296:   } else {
    297:     for (j = 0; j < s; j++) w[j] = h / rk->dtratio * tab->b[j];
    298:     PetscCall(VecGetSubVector(X, rk->is_fast, &Xfast));
    299:     PetscCall(VecMAXPY(Xfast, s, w, rk->YdotRHS_fast));
    300:     PetscCall(VecRestoreSubVector(X, rk->is_fast, &Xfast));
    301:   }
    302:   PetscFunctionReturn(PETSC_SUCCESS);
    303: }

    305: static PetscErrorCode TSStepRefine_RK_MultirateSplit(TS ts)
    306: {
    307:   TS_RK           *rk         = (TS_RK *)ts->data;
    308:   TS               subts_fast = rk->subts_fast, currentlevelts;
    309:   TS_RK           *subrk_fast = (TS_RK *)subts_fast->data;
    310:   RKTableau        tab        = rk->tableau;
    311:   Vec             *Y          = rk->Y;
    312:   Vec             *YdotRHS = rk->YdotRHS, *YdotRHS_fast = rk->YdotRHS_fast;
    313:   Vec              Yfast, Xfast;
    314:   const PetscInt   s = tab->s;
    315:   const PetscReal *A = tab->A, *c = tab->c;
    316:   PetscScalar     *w = rk->work;
    317:   PetscInt         i, j, k;
    318:   PetscReal        t = ts->ptime, h = ts->time_step;

    320:   PetscFunctionBegin;
    321:   for (k = 0; k < rk->dtratio; k++) {
    322:     PetscCall(VecGetSubVector(ts->vec_sol, rk->is_fast, &Xfast));
    323:     for (i = 0; i < s; i++) PetscCall(VecGetSubVector(YdotRHS[i], rk->is_fast, &YdotRHS_fast[i]));
    324:     /* propagate fast component using small time steps */
    325:     for (i = 0; i < s; i++) {
    326:       /* stage value for slow components */
    327:       PetscCall(TSInterpolate_RK_MultirateSplit(rk->ts_root, t + k * h / rk->dtratio + h / rk->dtratio * c[i], Y[i]));
    328:       currentlevelts = rk->ts_root;
    329:       while (currentlevelts != ts) { /* all the slow parts need to be interpolated separated */
    330:         currentlevelts = ((TS_RK *)currentlevelts->data)->subts_fast;
    331:         PetscCall(TSInterpolate_RK_MultirateSplit(currentlevelts, t + k * h / rk->dtratio + h / rk->dtratio * c[i], Y[i]));
    332:       }
    333:       for (j = 0; j < i; j++) w[j] = h / rk->dtratio * A[i * s + j];
    334:       subrk_fast->stage_time = t + h / rk->dtratio * c[i];
    335:       PetscCall(TSPreStage(subts_fast, subrk_fast->stage_time));
    336:       /* stage value for fast components */
    337:       PetscCall(VecGetSubVector(Y[i], rk->is_fast, &Yfast));
    338:       PetscCall(VecCopy(Xfast, Yfast));
    339:       PetscCall(VecMAXPY(Yfast, i, w, YdotRHS_fast));
    340:       PetscCall(VecRestoreSubVector(Y[i], rk->is_fast, &Yfast));
    341:       PetscCall(TSPostStage(subts_fast, subrk_fast->stage_time, i, Y));
    342:       /* compute the stage RHS for fast components */
    343:       PetscCall(TSComputeRHSFunction(subts_fast, t + k * h * rk->dtratio + h / rk->dtratio * c[i], Y[i], YdotRHS_fast[i]));
    344:     }
    345:     PetscCall(VecRestoreSubVector(ts->vec_sol, rk->is_fast, &Xfast));
    346:     /* update the value of fast components using fast time step */
    347:     rk->slow = PETSC_FALSE;
    348:     PetscCall(TSEvaluateStep_RK_MultirateSplit(ts, tab->order, ts->vec_sol, NULL));
    349:     for (i = 0; i < s; i++) PetscCall(VecRestoreSubVector(YdotRHS[i], rk->is_fast, &YdotRHS_fast[i]));

    351:     if (subrk_fast->subts_fast) {
    352:       subts_fast->ptime     = t + k * h / rk->dtratio;
    353:       subts_fast->time_step = h / rk->dtratio;
    354:       PetscCall(TSStepRefine_RK_MultirateSplit(subts_fast));
    355:     }
    356:     /* update the fast components of the solution */
    357:     PetscCall(VecGetSubVector(ts->vec_sol, rk->is_fast, &Xfast));
    358:     PetscCall(VecISCopy(rk->X0, rk->is_fast, SCATTER_FORWARD, Xfast));
    359:     PetscCall(VecRestoreSubVector(ts->vec_sol, rk->is_fast, &Xfast));
    360:   }
    361:   PetscFunctionReturn(PETSC_SUCCESS);
    362: }

    364: static PetscErrorCode TSStep_RK_MultirateSplit(TS ts)
    365: {
    366:   TS_RK           *rk  = (TS_RK *)ts->data;
    367:   RKTableau        tab = rk->tableau;
    368:   Vec             *Y = rk->Y, *YdotRHS = rk->YdotRHS;
    369:   Vec             *YdotRHS_fast = rk->YdotRHS_fast, *YdotRHS_slow = rk->YdotRHS_slow;
    370:   Vec              Yslow, Yfast; /* subvectors store the stges of slow components and fast components respectively                           */
    371:   const PetscInt   s = tab->s;
    372:   const PetscReal *A = tab->A, *c = tab->c;
    373:   PetscScalar     *w = rk->work;
    374:   PetscInt         i, j;
    375:   PetscReal        next_time_step = ts->time_step, t = ts->ptime, h = ts->time_step;

    377:   PetscFunctionBegin;
    378:   rk->status = TS_STEP_INCOMPLETE;
    379:   for (i = 0; i < s; i++) {
    380:     PetscCall(VecGetSubVector(YdotRHS[i], rk->is_slow, &YdotRHS_slow[i]));
    381:     PetscCall(VecGetSubVector(YdotRHS[i], rk->is_fast, &YdotRHS_fast[i]));
    382:   }
    383:   PetscCall(VecCopy(ts->vec_sol, rk->X0));
    384:   /* propagate both slow and fast components using large time steps */
    385:   for (i = 0; i < s; i++) {
    386:     rk->stage_time = t + h * c[i];
    387:     PetscCall(TSPreStage(ts, rk->stage_time));
    388:     PetscCall(VecCopy(ts->vec_sol, Y[i]));
    389:     PetscCall(VecGetSubVector(Y[i], rk->is_fast, &Yfast));
    390:     PetscCall(VecGetSubVector(Y[i], rk->is_slow, &Yslow));
    391:     for (j = 0; j < i; j++) w[j] = h * A[i * s + j];
    392:     PetscCall(VecMAXPY(Yfast, i, w, YdotRHS_fast));
    393:     PetscCall(VecMAXPY(Yslow, i, w, YdotRHS_slow));
    394:     PetscCall(VecRestoreSubVector(Y[i], rk->is_fast, &Yfast));
    395:     PetscCall(VecRestoreSubVector(Y[i], rk->is_slow, &Yslow));
    396:     PetscCall(TSPostStage(ts, rk->stage_time, i, Y));
    397:     PetscCall(TSComputeRHSFunction(rk->subts_slow, t + h * c[i], Y[i], YdotRHS_slow[i]));
    398:     PetscCall(TSComputeRHSFunction(rk->subts_fast, t + h * c[i], Y[i], YdotRHS_fast[i]));
    399:   }
    400:   rk->slow = PETSC_TRUE;
    401:   /* update the slow components of the solution using slow time step */
    402:   PetscCall(TSEvaluateStep_RK_MultirateSplit(ts, tab->order, ts->vec_sol, NULL));
    403:   /* YdotRHS will be used for interpolation during refinement */
    404:   for (i = 0; i < s; i++) {
    405:     PetscCall(VecRestoreSubVector(YdotRHS[i], rk->is_slow, &YdotRHS_slow[i]));
    406:     PetscCall(VecRestoreSubVector(YdotRHS[i], rk->is_fast, &YdotRHS_fast[i]));
    407:   }

    409:   PetscCall(TSStepRefine_RK_MultirateSplit(ts));

    411:   ts->ptime     = t + ts->time_step;
    412:   ts->time_step = next_time_step;
    413:   rk->status    = TS_STEP_COMPLETE;
    414:   PetscFunctionReturn(PETSC_SUCCESS);
    415: }

    417: static PetscErrorCode TSSetUp_RK_MultirateSplit(TS ts)
    418: {
    419:   TS_RK *rk = (TS_RK *)ts->data, *nextlevelrk, *currentlevelrk;
    420:   TS     nextlevelts;
    421:   Vec    X0;

    423:   PetscFunctionBegin;
    424:   PetscCall(TSRHSSplitGetIS(ts, "slow", &rk->is_slow));
    425:   PetscCall(TSRHSSplitGetIS(ts, "fast", &rk->is_fast));
    426:   PetscCheck(rk->is_slow && rk->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 -ts_type bsi");
    427:   PetscCall(TSRHSSplitGetSubTS(ts, "slow", &rk->subts_slow));
    428:   PetscCall(TSRHSSplitGetSubTS(ts, "fast", &rk->subts_fast));
    429:   PetscCheck(rk->subts_slow && rk->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");

    431:   PetscCall(VecDuplicate(ts->vec_sol, &X0));
    432:   /* The TS at each level share the same tableau, work array, solution vector, stage values and stage derivatives */
    433:   currentlevelrk = rk;
    434:   while (currentlevelrk->subts_fast) {
    435:     PetscCall(PetscMalloc1(rk->tableau->s, &currentlevelrk->YdotRHS_fast));
    436:     PetscCall(PetscMalloc1(rk->tableau->s, &currentlevelrk->YdotRHS_slow));
    437:     PetscCall(PetscObjectReference((PetscObject)X0));
    438:     currentlevelrk->X0      = X0;
    439:     currentlevelrk->ts_root = ts;

    441:     /* set up the ts for the slow part */
    442:     nextlevelts = currentlevelrk->subts_slow;
    443:     PetscCall(PetscNew(&nextlevelrk));
    444:     nextlevelrk->tableau = rk->tableau;
    445:     nextlevelrk->work    = rk->work;
    446:     nextlevelrk->Y       = rk->Y;
    447:     nextlevelrk->YdotRHS = rk->YdotRHS;
    448:     nextlevelts->data    = (void *)nextlevelrk;
    449:     PetscCall(TSCopyDM(ts, nextlevelts));
    450:     PetscCall(TSSetSolution(nextlevelts, ts->vec_sol));

    452:     /* set up the ts for the fast part */
    453:     nextlevelts = currentlevelrk->subts_fast;
    454:     PetscCall(PetscNew(&nextlevelrk));
    455:     nextlevelrk->tableau = rk->tableau;
    456:     nextlevelrk->work    = rk->work;
    457:     nextlevelrk->Y       = rk->Y;
    458:     nextlevelrk->YdotRHS = rk->YdotRHS;
    459:     nextlevelrk->dtratio = rk->dtratio;
    460:     PetscCall(TSRHSSplitGetIS(nextlevelts, "slow", &nextlevelrk->is_slow));
    461:     PetscCall(TSRHSSplitGetSubTS(nextlevelts, "slow", &nextlevelrk->subts_slow));
    462:     PetscCall(TSRHSSplitGetIS(nextlevelts, "fast", &nextlevelrk->is_fast));
    463:     PetscCall(TSRHSSplitGetSubTS(nextlevelts, "fast", &nextlevelrk->subts_fast));
    464:     nextlevelts->data = (void *)nextlevelrk;
    465:     PetscCall(TSCopyDM(ts, nextlevelts));
    466:     PetscCall(TSSetSolution(nextlevelts, ts->vec_sol));

    468:     currentlevelrk = nextlevelrk;
    469:   }
    470:   PetscCall(VecDestroy(&X0));

    472:   ts->ops->step         = TSStep_RK_MultirateSplit;
    473:   ts->ops->evaluatestep = TSEvaluateStep_RK_MultirateSplit;
    474:   ts->ops->interpolate  = TSInterpolate_RK_MultirateSplit;
    475:   PetscFunctionReturn(PETSC_SUCCESS);
    476: }

    478: PetscErrorCode TSRKSetMultirate_RK(TS ts, PetscBool use_multirate)
    479: {
    480:   TS_RK *rk = (TS_RK *)ts->data;

    482:   PetscFunctionBegin;
    484:   rk->use_multirate = use_multirate;
    485:   if (use_multirate) {
    486:     rk->dtratio = 2;
    487:     PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSSetUp_RK_MultirateSplit_C", TSSetUp_RK_MultirateSplit));
    488:     PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSReset_RK_MultirateSplit_C", TSReset_RK_MultirateSplit));
    489:     PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSSetUp_RK_MultirateNonsplit_C", TSSetUp_RK_MultirateNonsplit));
    490:     PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSReset_RK_MultirateNonsplit_C", TSReset_RK_MultirateNonsplit));
    491:   } else {
    492:     rk->dtratio = 0;
    493:     PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSSetUp_RK_MultirateSplit_C", NULL));
    494:     PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSReset_RK_MultirateSplit_C", NULL));
    495:     PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSSetUp_RK_MultirateNonsplit_C", NULL));
    496:     PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSReset_RK_MultirateNonsplit_C", NULL));
    497:   }
    498:   PetscFunctionReturn(PETSC_SUCCESS);
    499: }

    501: PetscErrorCode TSRKGetMultirate_RK(TS ts, PetscBool *use_multirate)
    502: {
    503:   TS_RK *rk = (TS_RK *)ts->data;

    505:   PetscFunctionBegin;
    507:   *use_multirate = rk->use_multirate;
    508:   PetscFunctionReturn(PETSC_SUCCESS);
    509: }