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:   VecDestroy(&rk->X0);
 25:   VecDestroyVecs(tab->s, &rk->YdotRHS_slow);
 26:   return 0;
 27: }

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

 39:   t = (itime - rk->ptime) / h;
 40:   PetscMalloc1(s, &b);
 41:   for (i = 0; i < s; i++) b[i] = 0;
 42:   for (j = 0, tt = t; j < p; j++, tt *= t) {
 43:     for (i = 0; i < s; i++) b[i] += h * B[i * p + j] * tt;
 44:   }
 45:   VecCopy(rk->X0, X);
 46:   VecMAXPY(X, s, b, rk->YdotRHS_slow);
 47:   PetscFree(b);
 48:   return 0;
 49: }

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

 64:   VecDuplicate(ts->vec_sol, &vec_fast);
 65:   previousts = rk->subts_current;
 66:   TSRHSSplitGetSubTS(rk->subts_current, "fast", &subts);
 67:   TSRHSSplitGetSubTS(subts, "fast", &subts);
 68:   for (k = 0; k < rk->dtratio; k++) {
 69:     for (i = 0; i < s; i++) {
 70:       TSInterpolate_RK_MultirateNonsplit(ts, t + k * h / rk->dtratio + h / rk->dtratio * c[i], Y[i]);
 71:       for (j = 0; j < i; j++) w[j] = h / rk->dtratio * A[i * s + j];
 72:       /* 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 */
 73:       VecCopy(ts->vec_sol, vec_fast);
 74:       VecMAXPY(vec_fast, i, w, YdotRHS);
 75:       /* update the fast components in the stage value */
 76:       VecGetSubVector(vec_fast, rk->is_fast, &subvec_fast);
 77:       VecISCopy(Y[i], rk->is_fast, SCATTER_FORWARD, subvec_fast);
 78:       VecRestoreSubVector(vec_fast, rk->is_fast, &subvec_fast);
 79:       /* compute the stage RHS */
 80:       TSComputeRHSFunction(ts, t + k * h / rk->dtratio + h / rk->dtratio * c[i], Y[i], YdotRHS[i]);
 81:     }
 82:     VecCopy(ts->vec_sol, vec_fast);
 83:     TSEvaluateStep(ts, tab->order, vec_fast, NULL);
 84:     /* update the fast components in the solution */
 85:     VecGetSubVector(vec_fast, rk->is_fast, &subvec_fast);
 86:     VecISCopy(ts->vec_sol, rk->is_fast, SCATTER_FORWARD, subvec_fast);
 87:     VecRestoreSubVector(vec_fast, rk->is_fast, &subvec_fast);

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

101:       TSStepRefine_RK_MultirateNonsplit(ts);

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

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

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

154:   rk->subts_current = ts;
155:   rk->ptime         = t;
156:   rk->time_step     = h;
157:   TSStepRefine_RK_MultirateNonsplit(ts);

159:   ts->ptime     = t + ts->time_step;
160:   ts->time_step = next_time_step;
161:   rk->status    = TS_STEP_COMPLETE;

163:   /* free memory of work vectors */
164:   VecDestroy(&stage_slow);
165:   VecDestroy(&sol_slow);
166:   return 0;
167: }

169: static PetscErrorCode TSSetUp_RK_MultirateNonsplit(TS ts)
170: {
171:   TS_RK    *rk  = (TS_RK *)ts->data;
172:   RKTableau tab = rk->tableau;

174:   TSRHSSplitGetIS(ts, "slow", &rk->is_slow);
175:   TSRHSSplitGetIS(ts, "fast", &rk->is_fast);
177:   TSRHSSplitGetSubTS(ts, "slow", &rk->subts_slow);
178:   TSRHSSplitGetSubTS(ts, "fast", &rk->subts_fast);
180:   VecDuplicate(ts->vec_sol, &rk->X0);
181:   VecDuplicateVecs(ts->vec_sol, tab->s, &rk->YdotRHS_slow);
182:   rk->subts_current = rk->subts_fast;

184:   ts->ops->step        = TSStep_RK_MultirateNonsplit;
185:   ts->ops->interpolate = TSInterpolate_RK_MultirateNonsplit;
186:   return 0;
187: }

189: /*
190:   Copy DM from tssrc to tsdest, while keeping the original DMTS and DMSNES in tsdest.
191: */
192: static PetscErrorCode TSCopyDM(TS tssrc, TS tsdest)
193: {
194:   DM newdm, dmsrc, dmdest;

196:   TSGetDM(tssrc, &dmsrc);
197:   DMClone(dmsrc, &newdm);
198:   TSGetDM(tsdest, &dmdest);
199:   DMCopyDMTS(dmdest, newdm);
200:   DMCopyDMSNES(dmdest, newdm);
201:   TSSetDM(tsdest, newdm);
202:   DMDestroy(&newdm);
203:   return 0;
204: }

206: static PetscErrorCode TSReset_RK_MultirateSplit(TS ts)
207: {
208:   TS_RK *rk = (TS_RK *)ts->data;

210:   if (rk->subts_slow) {
211:     PetscFree(rk->subts_slow->data);
212:     rk->subts_slow = NULL;
213:   }
214:   if (rk->subts_fast) {
215:     PetscFree(rk->YdotRHS_fast);
216:     PetscFree(rk->YdotRHS_slow);
217:     VecDestroy(&rk->X0);
218:     TSReset_RK_MultirateSplit(rk->subts_fast);
219:     PetscFree(rk->subts_fast->data);
220:     rk->subts_fast = NULL;
221:   }
222:   return 0;
223: }

225: static PetscErrorCode TSInterpolate_RK_MultirateSplit(TS ts, PetscReal itime, Vec X)
226: {
227:   TS_RK           *rk = (TS_RK *)ts->data;
228:   Vec              Xslow;
229:   PetscInt         s = rk->tableau->s, p = rk->tableau->p, i, j;
230:   PetscReal        h;
231:   PetscReal        tt, t;
232:   PetscScalar     *b;
233:   const PetscReal *B = rk->tableau->binterp;


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

265: /*
266:  This is for partitioned RHS multirate RK method
267:  The step completion formula is

269:  x1 = x0 + h b^T YdotRHS

271: */
272: static PetscErrorCode TSEvaluateStep_RK_MultirateSplit(TS ts, PetscInt order, Vec X, PetscBool *done)
273: {
274:   TS_RK       *rk  = (TS_RK *)ts->data;
275:   RKTableau    tab = rk->tableau;
276:   Vec          Xslow, Xfast; /* subvectors of X which store slow components and fast components respectively */
277:   PetscScalar *w = rk->work;
278:   PetscReal    h = ts->time_step;
279:   PetscInt     s = tab->s, j;

281:   VecCopy(ts->vec_sol, X);
282:   if (rk->slow) {
283:     for (j = 0; j < s; j++) w[j] = h * tab->b[j];
284:     VecGetSubVector(ts->vec_sol, rk->is_slow, &Xslow);
285:     VecMAXPY(Xslow, s, w, rk->YdotRHS_slow);
286:     VecRestoreSubVector(ts->vec_sol, rk->is_slow, &Xslow);
287:   } else {
288:     for (j = 0; j < s; j++) w[j] = h / rk->dtratio * tab->b[j];
289:     VecGetSubVector(X, rk->is_fast, &Xfast);
290:     VecMAXPY(Xfast, s, w, rk->YdotRHS_fast);
291:     VecRestoreSubVector(X, rk->is_fast, &Xfast);
292:   }
293:   return 0;
294: }

296: static PetscErrorCode TSStepRefine_RK_MultirateSplit(TS ts)
297: {
298:   TS_RK           *rk         = (TS_RK *)ts->data;
299:   TS               subts_fast = rk->subts_fast, currentlevelts;
300:   TS_RK           *subrk_fast = (TS_RK *)subts_fast->data;
301:   RKTableau        tab        = rk->tableau;
302:   Vec             *Y          = rk->Y;
303:   Vec             *YdotRHS = rk->YdotRHS, *YdotRHS_fast = rk->YdotRHS_fast;
304:   Vec              Yfast, Xfast;
305:   const PetscInt   s = tab->s;
306:   const PetscReal *A = tab->A, *c = tab->c;
307:   PetscScalar     *w = rk->work;
308:   PetscInt         i, j, k;
309:   PetscReal        t = ts->ptime, h = ts->time_step;

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

341:     if (subrk_fast->subts_fast) {
342:       subts_fast->ptime     = t + k * h / rk->dtratio;
343:       subts_fast->time_step = h / rk->dtratio;
344:       TSStepRefine_RK_MultirateSplit(subts_fast);
345:     }
346:     /* update the fast components of the solution */
347:     VecGetSubVector(ts->vec_sol, rk->is_fast, &Xfast);
348:     VecISCopy(rk->X0, rk->is_fast, SCATTER_FORWARD, Xfast);
349:     VecRestoreSubVector(ts->vec_sol, rk->is_fast, &Xfast);
350:   }
351:   return 0;
352: }

354: static PetscErrorCode TSStep_RK_MultirateSplit(TS ts)
355: {
356:   TS_RK           *rk  = (TS_RK *)ts->data;
357:   RKTableau        tab = rk->tableau;
358:   Vec             *Y = rk->Y, *YdotRHS = rk->YdotRHS;
359:   Vec             *YdotRHS_fast = rk->YdotRHS_fast, *YdotRHS_slow = rk->YdotRHS_slow;
360:   Vec              Yslow, Yfast; /* subvectors store the stges of slow components and fast components respectively                           */
361:   const PetscInt   s = tab->s;
362:   const PetscReal *A = tab->A, *c = tab->c;
363:   PetscScalar     *w = rk->work;
364:   PetscInt         i, j;
365:   PetscReal        next_time_step = ts->time_step, t = ts->ptime, h = ts->time_step;

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

398:   TSStepRefine_RK_MultirateSplit(ts);

400:   ts->ptime     = t + ts->time_step;
401:   ts->time_step = next_time_step;
402:   rk->status    = TS_STEP_COMPLETE;
403:   return 0;
404: }

406: static PetscErrorCode TSSetUp_RK_MultirateSplit(TS ts)
407: {
408:   TS_RK *rk = (TS_RK *)ts->data, *nextlevelrk, *currentlevelrk;
409:   TS     nextlevelts;
410:   Vec    X0;

412:   TSRHSSplitGetIS(ts, "slow", &rk->is_slow);
413:   TSRHSSplitGetIS(ts, "fast", &rk->is_fast);
415:   TSRHSSplitGetSubTS(ts, "slow", &rk->subts_slow);
416:   TSRHSSplitGetSubTS(ts, "fast", &rk->subts_fast);

419:   VecDuplicate(ts->vec_sol, &X0);
420:   /* The TS at each level share the same tableau, work array, solution vector, stage values and stage derivatives */
421:   currentlevelrk = rk;
422:   while (currentlevelrk->subts_fast) {
423:     PetscMalloc1(rk->tableau->s, &currentlevelrk->YdotRHS_fast);
424:     PetscMalloc1(rk->tableau->s, &currentlevelrk->YdotRHS_slow);
425:     PetscObjectReference((PetscObject)X0);
426:     currentlevelrk->X0      = X0;
427:     currentlevelrk->ts_root = ts;

429:     /* set up the ts for the slow part */
430:     nextlevelts = currentlevelrk->subts_slow;
431:     PetscNew(&nextlevelrk);
432:     nextlevelrk->tableau = rk->tableau;
433:     nextlevelrk->work    = rk->work;
434:     nextlevelrk->Y       = rk->Y;
435:     nextlevelrk->YdotRHS = rk->YdotRHS;
436:     nextlevelts->data    = (void *)nextlevelrk;
437:     TSCopyDM(ts, nextlevelts);
438:     TSSetSolution(nextlevelts, ts->vec_sol);

440:     /* set up the ts for the fast part */
441:     nextlevelts = currentlevelrk->subts_fast;
442:     PetscNew(&nextlevelrk);
443:     nextlevelrk->tableau = rk->tableau;
444:     nextlevelrk->work    = rk->work;
445:     nextlevelrk->Y       = rk->Y;
446:     nextlevelrk->YdotRHS = rk->YdotRHS;
447:     nextlevelrk->dtratio = rk->dtratio;
448:     TSRHSSplitGetIS(nextlevelts, "slow", &nextlevelrk->is_slow);
449:     TSRHSSplitGetSubTS(nextlevelts, "slow", &nextlevelrk->subts_slow);
450:     TSRHSSplitGetIS(nextlevelts, "fast", &nextlevelrk->is_fast);
451:     TSRHSSplitGetSubTS(nextlevelts, "fast", &nextlevelrk->subts_fast);
452:     nextlevelts->data = (void *)nextlevelrk;
453:     TSCopyDM(ts, nextlevelts);
454:     TSSetSolution(nextlevelts, ts->vec_sol);

456:     currentlevelrk = nextlevelrk;
457:   }
458:   VecDestroy(&X0);

460:   ts->ops->step         = TSStep_RK_MultirateSplit;
461:   ts->ops->evaluatestep = TSEvaluateStep_RK_MultirateSplit;
462:   ts->ops->interpolate  = TSInterpolate_RK_MultirateSplit;
463:   return 0;
464: }

466: PetscErrorCode TSRKSetMultirate_RK(TS ts, PetscBool use_multirate)
467: {
468:   TS_RK *rk = (TS_RK *)ts->data;

471:   rk->use_multirate = use_multirate;
472:   if (use_multirate) {
473:     rk->dtratio = 2;
474:     PetscObjectComposeFunction((PetscObject)ts, "TSSetUp_RK_MultirateSplit_C", TSSetUp_RK_MultirateSplit);
475:     PetscObjectComposeFunction((PetscObject)ts, "TSReset_RK_MultirateSplit_C", TSReset_RK_MultirateSplit);
476:     PetscObjectComposeFunction((PetscObject)ts, "TSSetUp_RK_MultirateNonsplit_C", TSSetUp_RK_MultirateNonsplit);
477:     PetscObjectComposeFunction((PetscObject)ts, "TSReset_RK_MultirateNonsplit_C", TSReset_RK_MultirateNonsplit);
478:   } else {
479:     rk->dtratio = 0;
480:     PetscObjectComposeFunction((PetscObject)ts, "TSSetUp_RK_MultirateSplit_C", NULL);
481:     PetscObjectComposeFunction((PetscObject)ts, "TSReset_RK_MultirateSplit_C", NULL);
482:     PetscObjectComposeFunction((PetscObject)ts, "TSSetUp_RK_MultirateNonsplit_C", NULL);
483:     PetscObjectComposeFunction((PetscObject)ts, "TSReset_RK_MultirateNonsplit_C", NULL);
484:   }
485:   return 0;
486: }

488: PetscErrorCode TSRKGetMultirate_RK(TS ts, PetscBool *use_multirate)
489: {
490:   TS_RK *rk = (TS_RK *)ts->data;

493:   *use_multirate = rk->use_multirate;
494:   return 0;
495: }