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

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

103:       TSStepRefine_RK_MultirateNonsplit(ts);

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

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

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

158:   rk->subts_current = ts;
159:   rk->ptime = t;
160:   rk->time_step = h;
161:   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:   VecDestroy(&stage_slow);
169:   VecDestroy(&sol_slow);
170:   return 0;
171: }

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

178:   TSRHSSplitGetIS(ts,"slow",&rk->is_slow);
179:   TSRHSSplitGetIS(ts,"fast",&rk->is_fast);
181:   TSRHSSplitGetSubTS(ts,"slow",&rk->subts_slow);
182:   TSRHSSplitGetSubTS(ts,"fast",&rk->subts_fast);
184:   VecDuplicate(ts->vec_sol,&rk->X0);
185:   VecDuplicateVecs(ts->vec_sol,tab->s,&rk->YdotRHS_slow);
186:   rk->subts_current = rk->subts_fast;

188:   ts->ops->step        = TSStep_RK_MultirateNonsplit;
189:   ts->ops->interpolate = TSInterpolate_RK_MultirateNonsplit;
190:   return 0;
191: }

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

200:   TSGetDM(tssrc,&dmsrc);
201:   DMClone(dmsrc,&newdm);
202:   TSGetDM(tsdest,&dmdest);
203:   DMCopyDMTS(dmdest,newdm);
204:   DMCopyDMSNES(dmdest,newdm);
205:   TSSetDM(tsdest,newdm);
206:   DMDestroy(&newdm);
207:   return 0;
208: }

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

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

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


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

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

278:  x1 = x0 + h b^T YdotRHS

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

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

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

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

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

411:   TSStepRefine_RK_MultirateSplit(ts);

413:   ts->ptime = t + ts->time_step;
414:   ts->time_step = next_time_step;
415:   rk->status = TS_STEP_COMPLETE;
416:   return 0;
417: }

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

425:   TSRHSSplitGetIS(ts,"slow",&rk->is_slow);
426:   TSRHSSplitGetIS(ts,"fast",&rk->is_fast);
428:   TSRHSSplitGetSubTS(ts,"slow",&rk->subts_slow);
429:   TSRHSSplitGetSubTS(ts,"fast",&rk->subts_fast);

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

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

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

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

473:   ts->ops->step         = TSStep_RK_MultirateSplit;
474:   ts->ops->evaluatestep = TSEvaluateStep_RK_MultirateSplit;
475:   ts->ops->interpolate  = TSInterpolate_RK_MultirateSplit;
476:   return 0;
477: }

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

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

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

506:   *use_multirate = rk->use_multirate;
507:   return 0;
508: }