Actual source code: mrk.c

petsc-master 2019-08-22
Report Typos and Errors
  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;

 26:   VecDestroy(&rk->X0);
 27:   VecDestroyVecs(tab->s,&rk->YdotRHS_slow);
 28:   return(0);
 29: }

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

 42:   if (!B) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"TSRK %s does not have an interpolation formula",rk->tableau->name);
 43:   t = (itime - rk->ptime)/h;
 44:   PetscMalloc1(s,&b);
 45:   for (i=0; i<s; i++) b[i] = 0;
 46:   for (j=0,tt=t; j<p; j++,tt*=t) {
 47:     for (i=0; i<s; i++) {
 48:       b[i]  += h * B[i*p+j] * tt;
 49:     }
 50:   }
 51:   VecCopy(rk->X0,X);
 52:   VecMAXPY(X,s,b,rk->YdotRHS_slow);
 53:   PetscFree(b);
 54:   return(0);
 55: }

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

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

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

108:       TSStepRefine_RK_MultirateNonsplit(ts);

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

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

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

165:   rk->subts_current = ts;
166:   rk->ptime = t;
167:   rk->time_step = h;
168:   TSStepRefine_RK_MultirateNonsplit(ts);

170:   ts->ptime = t + ts->time_step;
171:   ts->time_step = next_time_step;
172:   rk->status = TS_STEP_COMPLETE;

174:   /* free memory of work vectors */
175:   VecDestroy(&stage_slow);
176:   VecDestroy(&sol_slow);
177:   return(0);
178: }

180: static PetscErrorCode TSSetUp_RK_MultirateNonsplit(TS ts)
181: {
182:   TS_RK          *rk = (TS_RK*)ts->data;
183:   RKTableau      tab = rk->tableau;

187:   TSRHSSplitGetIS(ts,"slow",&rk->is_slow);
188:   TSRHSSplitGetIS(ts,"fast",&rk->is_fast);
189:   if (!rk->is_slow || !rk->is_fast) SETERRQ(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");
190:   TSRHSSplitGetSubTS(ts,"slow",&rk->subts_slow);
191:   TSRHSSplitGetSubTS(ts,"fast",&rk->subts_fast);
192:   if (!rk->subts_slow || !rk->subts_fast) SETERRQ(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");
193:   VecDuplicate(ts->vec_sol,&rk->X0);
194:   VecDuplicateVecs(ts->vec_sol,tab->s,&rk->YdotRHS_slow);
195:   rk->subts_current = rk->subts_fast;

197:   ts->ops->step        = TSStep_RK_MultirateNonsplit;
198:   ts->ops->interpolate = TSInterpolate_RK_MultirateNonsplit;
199:   return(0);
200: }

202: /*
203:   Copy DM from tssrc to tsdest, while keeping the original DMTS and DMSNES in tsdest.
204: */
205: static PetscErrorCode TSCopyDM(TS tssrc,TS tsdest)
206: {
207:   DM             newdm,dmsrc,dmdest;

211:   TSGetDM(tssrc,&dmsrc);
212:   DMClone(dmsrc,&newdm);
213:   TSGetDM(tsdest,&dmdest);
214:   DMCopyDMTS(dmdest,newdm);
215:   DMCopyDMSNES(dmdest,newdm);
216:   TSSetDM(tsdest,newdm);
217:   DMDestroy(&newdm);
218:   return(0);
219: }

221: static PetscErrorCode TSReset_RK_MultirateSplit(TS ts)
222: {
223:   TS_RK          *rk = (TS_RK*)ts->data;

227:   if (rk->subts_slow) {
228:     PetscFree(rk->subts_slow->data);
229:     rk->subts_slow = NULL;
230:   }
231:   if (rk->subts_fast) {
232:     PetscFree(rk->YdotRHS_fast);
233:     PetscFree(rk->YdotRHS_slow);
234:     VecDestroy(&rk->X0);
235:     TSReset_RK_MultirateSplit(rk->subts_fast);
236:     PetscFree(rk->subts_fast->data);
237:     rk->subts_fast = NULL;
238:   }
239:   return(0);
240: }

242: static PetscErrorCode TSInterpolate_RK_MultirateSplit(TS ts,PetscReal itime,Vec X)
243: {
244:   TS_RK           *rk = (TS_RK*)ts->data;
245:   Vec             Xslow;
246:   PetscInt        s = rk->tableau->s,p = rk->tableau->p,i,j;
247:   PetscReal       h;
248:   PetscReal       tt,t;
249:   PetscScalar     *b;
250:   const PetscReal *B = rk->tableau->binterp;
251:   PetscErrorCode  ierr;

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

256:   switch (rk->status) {
257:     case TS_STEP_INCOMPLETE:
258:     case TS_STEP_PENDING:
259:       h = ts->time_step;
260:       t = (itime - ts->ptime)/h;
261:       break;
262:     case TS_STEP_COMPLETE:
263:       h = ts->ptime - ts->ptime_prev;
264:       t = (itime - ts->ptime)/h + 1; /* In the interval [0,1] */
265:       break;
266:     default: SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_PLIB,"Invalid TSStepStatus");
267:   }
268:   PetscMalloc1(s,&b);
269:   for (i=0; i<s; i++) b[i] = 0;
270:   for (j=0,tt=t; j<p; j++,tt*=t) {
271:     for (i=0; i<s; i++) {
272:       b[i]  += h * B[i*p+j] * tt;
273:     }
274:   }
275:   for (i=0; i<s; i++) {
276:     VecGetSubVector(rk->YdotRHS[i],rk->is_slow,&rk->YdotRHS_slow[i]);
277:   }
278:   VecGetSubVector(X,rk->is_slow,&Xslow);
279:   VecISCopy(rk->X0,rk->is_slow,SCATTER_REVERSE,Xslow);
280:   VecMAXPY(Xslow,s,b,rk->YdotRHS_slow);
281:   VecRestoreSubVector(X,rk->is_slow,&Xslow);
282:   for (i=0; i<s; i++) {
283:     VecRestoreSubVector(rk->YdotRHS[i],rk->is_slow,&rk->YdotRHS_slow[i]);
284:   }
285:   PetscFree(b);
286:   return(0);
287: }

289: /*
290:  This is for partitioned RHS multirate RK method
291:  The step completion formula is

293:  x1 = x0 + h b^T YdotRHS

295: */
296: static PetscErrorCode TSEvaluateStep_RK_MultirateSplit(TS ts,PetscInt order,Vec X,PetscBool *done)
297: {
298:   TS_RK          *rk = (TS_RK*)ts->data;
299:   RKTableau      tab = rk->tableau;
300:   Vec            Xslow,Xfast;                  /* subvectors of X which store slow components and fast components respectively */
301:   PetscScalar    *w = rk->work;
302:   PetscReal      h = ts->time_step;
303:   PetscInt       s = tab->s,j;

307:   VecCopy(ts->vec_sol,X);
308:   if (rk->slow) {
309:     for (j=0; j<s; j++) w[j] = h*tab->b[j];
310:     VecGetSubVector(ts->vec_sol,rk->is_slow,&Xslow);
311:     VecMAXPY(Xslow,s,w,rk->YdotRHS_slow);
312:     VecRestoreSubVector(ts->vec_sol,rk->is_slow,&Xslow);;
313:   } else {
314:     for (j=0; j<s; j++) w[j] = h/rk->dtratio*tab->b[j];
315:     VecGetSubVector(X,rk->is_fast,&Xfast);
316:     VecMAXPY(Xfast,s,w,rk->YdotRHS_fast);
317:     VecRestoreSubVector(X,rk->is_fast,&Xfast);
318:   }
319:   return(0);
320: }

322: static PetscErrorCode TSStepRefine_RK_MultirateSplit(TS ts)
323: {
324:   TS_RK           *rk = (TS_RK*)ts->data;
325:   TS              subts_fast = rk->subts_fast,currentlevelts;
326:   TS_RK           *subrk_fast = (TS_RK*)subts_fast->data;
327:   RKTableau       tab = rk->tableau;
328:   Vec             *Y = rk->Y;
329:   Vec             *YdotRHS = rk->YdotRHS,*YdotRHS_fast = rk->YdotRHS_fast;
330:   Vec             Yfast,Xfast;
331:   const PetscInt  s = tab->s;
332:   const PetscReal *A = tab->A,*c = tab->c;
333:   PetscScalar     *w = rk->work;
334:   PetscInt        i,j,k;
335:   PetscReal       t = ts->ptime,h = ts->time_step;
336:   PetscErrorCode  ierr;

338:   for (k=0; k<rk->dtratio; k++) {
339:     VecGetSubVector(ts->vec_sol,rk->is_fast,&Xfast);
340:     for (i=0; i<s; i++) {
341:       VecGetSubVector(YdotRHS[i],rk->is_fast,&YdotRHS_fast[i]);
342:     }
343:     /* propogate fast component using small time steps */
344:     for (i=0; i<s; i++) {
345:       /* stage value for slow components */
346:       TSInterpolate_RK_MultirateSplit(rk->ts_root,t+k*h/rk->dtratio+h/rk->dtratio*c[i],Y[i]);
347:       currentlevelts = rk->ts_root;
348:       while (currentlevelts != ts) { /* all the slow parts need to be interpolated separated */
349:         currentlevelts = ((TS_RK*)currentlevelts->data)->subts_fast;
350:         TSInterpolate_RK_MultirateSplit(currentlevelts,t+k*h/rk->dtratio+h/rk->dtratio*c[i],Y[i]);
351:       }
352:       for (j=0; j<i; j++) w[j] = h/rk->dtratio*A[i*s+j];
353:       subrk_fast->stage_time = t + h/rk->dtratio*c[i];
354:       TSPreStage(subts_fast,subrk_fast->stage_time);
355:       /* stage value for fast components */
356:       VecGetSubVector(Y[i],rk->is_fast,&Yfast);
357:       VecCopy(Xfast,Yfast);
358:       VecMAXPY(Yfast,i,w,YdotRHS_fast);
359:       VecRestoreSubVector(Y[i],rk->is_fast,&Yfast);
360:       TSPostStage(subts_fast,subrk_fast->stage_time,i,Y);
361:       /* compute the stage RHS for fast components */
362:       TSComputeRHSFunction(subts_fast,t+k*h*rk->dtratio+h/rk->dtratio*c[i],Y[i],YdotRHS_fast[i]);
363:     }
364:     VecRestoreSubVector(ts->vec_sol,rk->is_fast,&Xfast);
365:     /* update the value of fast components using fast time step */
366:     rk->slow = PETSC_FALSE;
367:     TSEvaluateStep_RK_MultirateSplit(ts,tab->order,ts->vec_sol,NULL);
368:     for (i=0; i<s; i++) {
369:       VecRestoreSubVector(YdotRHS[i],rk->is_fast,&YdotRHS_fast[i]);
370:     }

372:     if (subrk_fast->subts_fast) {
373:       subts_fast->ptime = t+k*h/rk->dtratio;
374:       subts_fast->time_step = h/rk->dtratio;
375:       TSStepRefine_RK_MultirateSplit(subts_fast);
376:     }
377:     /* update the fast components of the solution */
378:     VecGetSubVector(ts->vec_sol,rk->is_fast,&Xfast);
379:     VecISCopy(rk->X0,rk->is_fast,SCATTER_FORWARD,Xfast);
380:     VecRestoreSubVector(ts->vec_sol,rk->is_fast,&Xfast);
381:   }
382:   return(0);
383: }

385: static PetscErrorCode TSStep_RK_MultirateSplit(TS ts)
386: {
387:   TS_RK           *rk = (TS_RK*)ts->data;
388:   RKTableau       tab = rk->tableau;
389:   Vec             *Y = rk->Y,*YdotRHS = rk->YdotRHS;
390:   Vec             *YdotRHS_fast = rk->YdotRHS_fast,*YdotRHS_slow = rk->YdotRHS_slow;
391:   Vec             Yslow,Yfast; /* subvectors store the stges of slow components and fast components respectively                           */
392:   const PetscInt  s = tab->s;
393:   const PetscReal *A = tab->A,*c = tab->c;
394:   PetscScalar     *w = rk->work;
395:   PetscInt        i,j;
396:   PetscReal       next_time_step = ts->time_step,t = ts->ptime,h = ts->time_step;
397:   PetscErrorCode  ierr;

400:   rk->status = TS_STEP_INCOMPLETE;
401:   for (i=0; i<s; i++) {
402:     VecGetSubVector(YdotRHS[i],rk->is_slow,&YdotRHS_slow[i]);
403:     VecGetSubVector(YdotRHS[i],rk->is_fast,&YdotRHS_fast[i]);
404:   }
405:   VecCopy(ts->vec_sol,rk->X0);
406:   /* propogate both slow and fast components using large time steps */
407:   for (i=0; i<s; i++) {
408:     rk->stage_time = t + h*c[i];
409:     TSPreStage(ts,rk->stage_time);
410:     VecCopy(ts->vec_sol,Y[i]);
411:     VecGetSubVector(Y[i],rk->is_fast,&Yfast);
412:     VecGetSubVector(Y[i],rk->is_slow,&Yslow);
413:     for (j=0; j<i; j++) w[j] = h*A[i*s+j];
414:     VecMAXPY(Yfast,i,w,YdotRHS_fast);
415:     VecMAXPY(Yslow,i,w,YdotRHS_slow);
416:     VecRestoreSubVector(Y[i],rk->is_fast,&Yfast);
417:     VecRestoreSubVector(Y[i],rk->is_slow,&Yslow);
418:     TSPostStage(ts,rk->stage_time,i,Y);
419:     TSComputeRHSFunction(rk->subts_slow,t+h*c[i],Y[i],YdotRHS_slow[i]);
420:     TSComputeRHSFunction(rk->subts_fast,t+h*c[i],Y[i],YdotRHS_fast[i]);
421:   }
422:   rk->slow = PETSC_TRUE;
423:   /* update the slow components of the solution using slow time step */
424:   TSEvaluateStep_RK_MultirateSplit(ts,tab->order,ts->vec_sol,NULL);
425:   /* YdotRHS will be used for interpolation during refinement */
426:   for (i=0; i<s; i++) {
427:     VecRestoreSubVector(YdotRHS[i],rk->is_slow,&YdotRHS_slow[i]);
428:     VecRestoreSubVector(YdotRHS[i],rk->is_fast,&YdotRHS_fast[i]);
429:   }

431:   TSStepRefine_RK_MultirateSplit(ts);

433:   ts->ptime = t + ts->time_step;
434:   ts->time_step = next_time_step;
435:   rk->status = TS_STEP_COMPLETE;
436:   return(0);
437: }

439: static PetscErrorCode TSSetUp_RK_MultirateSplit(TS ts)
440: {
441:   TS_RK          *rk = (TS_RK*)ts->data,*nextlevelrk,*currentlevelrk;
442:   TS             nextlevelts;
443:   Vec            X0;

447:   TSRHSSplitGetIS(ts,"slow",&rk->is_slow);
448:   TSRHSSplitGetIS(ts,"fast",&rk->is_fast);
449:   if (!rk->is_slow || !rk->is_fast) SETERRQ(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");
450:   TSRHSSplitGetSubTS(ts,"slow",&rk->subts_slow);
451:   TSRHSSplitGetSubTS(ts,"fast",&rk->subts_fast);
452:   if (!rk->subts_slow || !rk->subts_fast) SETERRQ(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");

454:   VecDuplicate(ts->vec_sol,&X0);
455:   /* The TS at each level share the same tableau, work array, solution vector, stage values and stage derivatives */
456:   currentlevelrk = rk;
457:   while (currentlevelrk->subts_fast) {
458:     PetscMalloc1(rk->tableau->s,&currentlevelrk->YdotRHS_fast);
459:     PetscMalloc1(rk->tableau->s,&currentlevelrk->YdotRHS_slow);
460:     PetscObjectReference((PetscObject)X0);
461:     currentlevelrk->X0 = X0;
462:     currentlevelrk->ts_root = ts;

464:     /* set up the ts for the slow part */
465:     nextlevelts = currentlevelrk->subts_slow;
466:     PetscNewLog(nextlevelts,&nextlevelrk);
467:     nextlevelrk->tableau = rk->tableau;
468:     nextlevelrk->work = rk->work;
469:     nextlevelrk->Y = rk->Y;
470:     nextlevelrk->YdotRHS = rk->YdotRHS;
471:     nextlevelts->data = (void*)nextlevelrk;
472:     TSCopyDM(ts,nextlevelts);
473:     TSSetSolution(nextlevelts,ts->vec_sol);

475:     /* set up the ts for the fast part */
476:     nextlevelts = currentlevelrk->subts_fast;
477:     PetscNewLog(nextlevelts,&nextlevelrk);
478:     nextlevelrk->tableau = rk->tableau;
479:     nextlevelrk->work = rk->work;
480:     nextlevelrk->Y = rk->Y;
481:     nextlevelrk->YdotRHS = rk->YdotRHS;
482:     nextlevelrk->dtratio = rk->dtratio;
483:     TSRHSSplitGetIS(nextlevelts,"slow",&nextlevelrk->is_slow);
484:     TSRHSSplitGetSubTS(nextlevelts,"slow",&nextlevelrk->subts_slow);
485:     TSRHSSplitGetIS(nextlevelts,"fast",&nextlevelrk->is_fast);
486:     TSRHSSplitGetSubTS(nextlevelts,"fast",&nextlevelrk->subts_fast);
487:     nextlevelts->data = (void*)nextlevelrk;
488:     TSCopyDM(ts,nextlevelts);
489:     TSSetSolution(nextlevelts,ts->vec_sol);

491:     currentlevelrk = nextlevelrk;
492:   }
493:   VecDestroy(&X0);

495:   ts->ops->step         = TSStep_RK_MultirateSplit;
496:   ts->ops->evaluatestep = TSEvaluateStep_RK_MultirateSplit;
497:   ts->ops->interpolate  = TSInterpolate_RK_MultirateSplit;
498:   return(0);
499: }

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

508:   rk->use_multirate = use_multirate;
509:   if (use_multirate) {
510:     rk->dtratio = 2;
511:     PetscObjectComposeFunction((PetscObject)ts,"TSSetUp_RK_MultirateSplit_C",TSSetUp_RK_MultirateSplit);
512:     PetscObjectComposeFunction((PetscObject)ts,"TSReset_RK_MultirateSplit_C",TSReset_RK_MultirateSplit);
513:     PetscObjectComposeFunction((PetscObject)ts,"TSSetUp_RK_MultirateNonsplit_C",TSSetUp_RK_MultirateNonsplit);
514:     PetscObjectComposeFunction((PetscObject)ts,"TSReset_RK_MultirateNonsplit_C",TSReset_RK_MultirateNonsplit);
515:   } else {
516:     rk->dtratio = 0;
517:     PetscObjectComposeFunction((PetscObject)ts,"TSSetUp_RK_MultirateSplit_C",NULL);
518:     PetscObjectComposeFunction((PetscObject)ts,"TSReset_RK_MultirateSplit_C",NULL);
519:     PetscObjectComposeFunction((PetscObject)ts,"TSSetUp_RK_MultirateNonsplit_C",NULL);
520:     PetscObjectComposeFunction((PetscObject)ts,"TSReset_RK_MultirateNonsplit_C",NULL);
521:   }
522:   return(0);
523: }

525: PetscErrorCode TSRKGetMultirate_RK(TS ts,PetscBool *use_multirate)
526: {
527:   TS_RK *rk = (TS_RK*)ts->data;

531:   *use_multirate = rk->use_multirate;
532:   return(0);
533: }