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;
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;
72: VecDuplicate(ts->vec_sol,&vec_fast);
73: previousts = rk->subts_current;
74: TSRHSSplitGetSubTS(rk->subts_current,"fast",&subts);
75: TSRHSSplitGetSubTS(subts,"fast",&subts);
76: for (k=0; k<rk->dtratio; k++) {
77: for (i=0; i<s; i++) {
78: TSInterpolate_RK_MultirateNonsplit(ts,t+k*h/rk->dtratio+h/rk->dtratio*c[i],Y[i]);
79: for (j=0; j<i; j++) w[j] = h/rk->dtratio*A[i*s+j];
80: /* 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 */
81: VecCopy(ts->vec_sol,vec_fast);
82: VecMAXPY(vec_fast,i,w,YdotRHS);
83: /* update the fast components in the stage value */
84: VecGetSubVector(vec_fast,rk->is_fast,&subvec_fast);
85: VecISCopy(Y[i],rk->is_fast,SCATTER_FORWARD,subvec_fast);
86: VecRestoreSubVector(vec_fast,rk->is_fast,&subvec_fast);
87: /* compute the stage RHS */
88: TSComputeRHSFunction(ts,t+k*h/rk->dtratio+h/rk->dtratio*c[i],Y[i],YdotRHS[i]);
89: }
90: VecCopy(ts->vec_sol,vec_fast);
91: TSEvaluateStep(ts,tab->order,vec_fast,NULL);
92: /* update the fast components in the solution */
93: VecGetSubVector(vec_fast,rk->is_fast,&subvec_fast);
94: VecISCopy(ts->vec_sol,rk->is_fast,SCATTER_FORWARD,subvec_fast);
95: VecRestoreSubVector(vec_fast,rk->is_fast,&subvec_fast);
97: if (subts) {
98: Vec *YdotRHS_copy;
99: VecDuplicateVecs(ts->vec_sol,s,&YdotRHS_copy);
100: rk->subts_current = rk->subts_fast;
101: ts->ptime = t+k*h/rk->dtratio;
102: ts->time_step = h/rk->dtratio;
103: TSRHSSplitGetIS(rk->subts_current,"fast",&rk->is_fast);
104: for (i=0; i<s; i++) {
105: VecCopy(rk->YdotRHS_slow[i],YdotRHS_copy[i]);
106: VecCopy(YdotRHS[i],rk->YdotRHS_slow[i]);
107: }
109: TSStepRefine_RK_MultirateNonsplit(ts);
111: rk->subts_current = previousts;
112: ts->ptime = t;
113: ts->time_step = h;
114: TSRHSSplitGetIS(previousts,"fast",&rk->is_fast);
115: for (i=0; i<s; i++) {
116: VecCopy(YdotRHS_copy[i],rk->YdotRHS_slow[i]);
117: }
118: VecDestroyVecs(s,&YdotRHS_copy);
119: }
120: }
121: VecDestroy(&vec_fast);
122: return(0);
123: }
125: static PetscErrorCode TSStep_RK_MultirateNonsplit(TS ts)
126: {
127: TS_RK *rk = (TS_RK*)ts->data;
128: RKTableau tab = rk->tableau;
129: Vec *Y = rk->Y,*YdotRHS = rk->YdotRHS,*YdotRHS_slow = rk->YdotRHS_slow;
130: Vec stage_slow,sol_slow; /* vectors store the slow components */
131: Vec subvec_slow; /* sub vector to store the slow components */
132: IS is_slow = rk->is_slow;
133: const PetscInt s = tab->s;
134: const PetscReal *A = tab->A,*c = tab->c;
135: PetscScalar *w = rk->work;
136: PetscInt i,j,dtratio = rk->dtratio;
137: PetscReal next_time_step = ts->time_step,t = ts->ptime,h = ts->time_step;
138: PetscErrorCode ierr;
141: rk->status = TS_STEP_INCOMPLETE;
142: VecDuplicate(ts->vec_sol,&stage_slow);
143: VecDuplicate(ts->vec_sol,&sol_slow);
144: VecCopy(ts->vec_sol,rk->X0);
145: for (i=0; i<s; i++) {
146: rk->stage_time = t + h*c[i];
147: TSPreStage(ts,rk->stage_time);
148: VecCopy(ts->vec_sol,Y[i]);
149: for (j=0; j<i; j++) w[j] = h*A[i*s+j];
150: VecMAXPY(Y[i],i,w,YdotRHS_slow);
151: TSPostStage(ts,rk->stage_time,i,Y);
152: /* compute the stage RHS */
153: TSComputeRHSFunction(ts,t+h*c[i],Y[i],YdotRHS_slow[i]);
154: }
155: /* update the slow components in the solution */
156: rk->YdotRHS = YdotRHS_slow;
157: rk->dtratio = 1;
158: TSEvaluateStep(ts,tab->order,sol_slow,NULL);
159: rk->dtratio = dtratio;
160: rk->YdotRHS = YdotRHS;
161: /* update the slow components in the solution */
162: VecGetSubVector(sol_slow,is_slow,&subvec_slow);
163: VecISCopy(ts->vec_sol,is_slow,SCATTER_FORWARD,subvec_slow);
164: VecRestoreSubVector(sol_slow,is_slow,&subvec_slow);
166: rk->subts_current = ts;
167: rk->ptime = t;
168: rk->time_step = h;
169: TSStepRefine_RK_MultirateNonsplit(ts);
171: ts->ptime = t + ts->time_step;
172: ts->time_step = next_time_step;
173: rk->status = TS_STEP_COMPLETE;
175: /* free memory of work vectors */
176: VecDestroy(&stage_slow);
177: VecDestroy(&sol_slow);
178: return(0);
179: }
181: static PetscErrorCode TSSetUp_RK_MultirateNonsplit(TS ts)
182: {
183: TS_RK *rk = (TS_RK*)ts->data;
184: RKTableau tab = rk->tableau;
188: TSRHSSplitGetIS(ts,"slow",&rk->is_slow);
189: TSRHSSplitGetIS(ts,"fast",&rk->is_fast);
190: 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");
191: TSRHSSplitGetSubTS(ts,"slow",&rk->subts_slow);
192: TSRHSSplitGetSubTS(ts,"fast",&rk->subts_fast);
193: 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");
194: VecDuplicate(ts->vec_sol,&rk->X0);
195: VecDuplicateVecs(ts->vec_sol,tab->s,&rk->YdotRHS_slow);
196: rk->subts_current = rk->subts_fast;
198: ts->ops->step = TSStep_RK_MultirateNonsplit;
199: ts->ops->interpolate = TSInterpolate_RK_MultirateNonsplit;
200: return(0);
201: }
203: /*
204: Copy DM from tssrc to tsdest, while keeping the original DMTS and DMSNES in tsdest.
205: */
206: static PetscErrorCode TSCopyDM(TS tssrc,TS tsdest)
207: {
208: DM newdm,dmsrc,dmdest;
212: TSGetDM(tssrc,&dmsrc);
213: DMClone(dmsrc,&newdm);
214: TSGetDM(tsdest,&dmdest);
215: DMCopyDMTS(dmdest,newdm);
216: DMCopyDMSNES(dmdest,newdm);
217: TSSetDM(tsdest,newdm);
218: DMDestroy(&newdm);
219: return(0);
220: }
222: static PetscErrorCode TSReset_RK_MultirateSplit(TS ts)
223: {
224: TS_RK *rk = (TS_RK*)ts->data;
228: if (rk->subts_slow) {
229: PetscFree(rk->subts_slow->data);
230: rk->subts_slow = NULL;
231: }
232: if (rk->subts_fast) {
233: PetscFree(rk->YdotRHS_fast);
234: PetscFree(rk->YdotRHS_slow);
235: VecDestroy(&rk->X0);
236: TSReset_RK_MultirateSplit(rk->subts_fast);
237: PetscFree(rk->subts_fast->data);
238: rk->subts_fast = NULL;
239: }
240: return(0);
241: }
243: static PetscErrorCode TSInterpolate_RK_MultirateSplit(TS ts,PetscReal itime,Vec X)
244: {
245: TS_RK *rk = (TS_RK*)ts->data;
246: Vec Xslow;
247: PetscInt s = rk->tableau->s,p = rk->tableau->p,i,j;
248: PetscReal h;
249: PetscReal tt,t;
250: PetscScalar *b;
251: const PetscReal *B = rk->tableau->binterp;
252: PetscErrorCode ierr;
255: if (!B) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"TSRK %s does not have an interpolation formula",rk->tableau->name);
257: switch (rk->status) {
258: case TS_STEP_INCOMPLETE:
259: case TS_STEP_PENDING:
260: h = ts->time_step;
261: t = (itime - ts->ptime)/h;
262: break;
263: case TS_STEP_COMPLETE:
264: h = ts->ptime - ts->ptime_prev;
265: t = (itime - ts->ptime)/h + 1; /* In the interval [0,1] */
266: break;
267: default: SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_PLIB,"Invalid TSStepStatus");
268: }
269: PetscMalloc1(s,&b);
270: for (i=0; i<s; i++) b[i] = 0;
271: for (j=0,tt=t; j<p; j++,tt*=t) {
272: for (i=0; i<s; i++) {
273: b[i] += h * B[i*p+j] * tt;
274: }
275: }
276: for (i=0; i<s; i++) {
277: VecGetSubVector(rk->YdotRHS[i],rk->is_slow,&rk->YdotRHS_slow[i]);
278: }
279: VecGetSubVector(X,rk->is_slow,&Xslow);
280: VecISCopy(rk->X0,rk->is_slow,SCATTER_REVERSE,Xslow);
281: VecMAXPY(Xslow,s,b,rk->YdotRHS_slow);
282: VecRestoreSubVector(X,rk->is_slow,&Xslow);
283: for (i=0; i<s; i++) {
284: VecRestoreSubVector(rk->YdotRHS[i],rk->is_slow,&rk->YdotRHS_slow[i]);
285: }
286: PetscFree(b);
287: return(0);
288: }
290: /*
291: This is for partitioned RHS multirate RK method
292: The step completion formula is
294: x1 = x0 + h b^T YdotRHS
296: */
297: static PetscErrorCode TSEvaluateStep_RK_MultirateSplit(TS ts,PetscInt order,Vec X,PetscBool *done)
298: {
299: TS_RK *rk = (TS_RK*)ts->data;
300: RKTableau tab = rk->tableau;
301: Vec Xslow,Xfast; /* subvectors of X which store slow components and fast components respectively */
302: PetscScalar *w = rk->work;
303: PetscReal h = ts->time_step;
304: PetscInt s = tab->s,j;
308: VecCopy(ts->vec_sol,X);
309: if (rk->slow) {
310: for (j=0; j<s; j++) w[j] = h*tab->b[j];
311: VecGetSubVector(ts->vec_sol,rk->is_slow,&Xslow);
312: VecMAXPY(Xslow,s,w,rk->YdotRHS_slow);
313: VecRestoreSubVector(ts->vec_sol,rk->is_slow,&Xslow);
314: } else {
315: for (j=0; j<s; j++) w[j] = h/rk->dtratio*tab->b[j];
316: VecGetSubVector(X,rk->is_fast,&Xfast);
317: VecMAXPY(Xfast,s,w,rk->YdotRHS_fast);
318: VecRestoreSubVector(X,rk->is_fast,&Xfast);
319: }
320: return(0);
321: }
323: static PetscErrorCode TSStepRefine_RK_MultirateSplit(TS ts)
324: {
325: TS_RK *rk = (TS_RK*)ts->data;
326: TS subts_fast = rk->subts_fast,currentlevelts;
327: TS_RK *subrk_fast = (TS_RK*)subts_fast->data;
328: RKTableau tab = rk->tableau;
329: Vec *Y = rk->Y;
330: Vec *YdotRHS = rk->YdotRHS,*YdotRHS_fast = rk->YdotRHS_fast;
331: Vec Yfast,Xfast;
332: const PetscInt s = tab->s;
333: const PetscReal *A = tab->A,*c = tab->c;
334: PetscScalar *w = rk->work;
335: PetscInt i,j,k;
336: PetscReal t = ts->ptime,h = ts->time_step;
337: PetscErrorCode ierr;
340: for (k=0; k<rk->dtratio; k++) {
341: VecGetSubVector(ts->vec_sol,rk->is_fast,&Xfast);
342: for (i=0; i<s; i++) {
343: VecGetSubVector(YdotRHS[i],rk->is_fast,&YdotRHS_fast[i]);
344: }
345: /* propagate fast component using small time steps */
346: for (i=0; i<s; i++) {
347: /* stage value for slow components */
348: TSInterpolate_RK_MultirateSplit(rk->ts_root,t+k*h/rk->dtratio+h/rk->dtratio*c[i],Y[i]);
349: currentlevelts = rk->ts_root;
350: while (currentlevelts != ts) { /* all the slow parts need to be interpolated separated */
351: currentlevelts = ((TS_RK*)currentlevelts->data)->subts_fast;
352: TSInterpolate_RK_MultirateSplit(currentlevelts,t+k*h/rk->dtratio+h/rk->dtratio*c[i],Y[i]);
353: }
354: for (j=0; j<i; j++) w[j] = h/rk->dtratio*A[i*s+j];
355: subrk_fast->stage_time = t + h/rk->dtratio*c[i];
356: TSPreStage(subts_fast,subrk_fast->stage_time);
357: /* stage value for fast components */
358: VecGetSubVector(Y[i],rk->is_fast,&Yfast);
359: VecCopy(Xfast,Yfast);
360: VecMAXPY(Yfast,i,w,YdotRHS_fast);
361: VecRestoreSubVector(Y[i],rk->is_fast,&Yfast);
362: TSPostStage(subts_fast,subrk_fast->stage_time,i,Y);
363: /* compute the stage RHS for fast components */
364: TSComputeRHSFunction(subts_fast,t+k*h*rk->dtratio+h/rk->dtratio*c[i],Y[i],YdotRHS_fast[i]);
365: }
366: VecRestoreSubVector(ts->vec_sol,rk->is_fast,&Xfast);
367: /* update the value of fast components using fast time step */
368: rk->slow = PETSC_FALSE;
369: TSEvaluateStep_RK_MultirateSplit(ts,tab->order,ts->vec_sol,NULL);
370: for (i=0; i<s; i++) {
371: VecRestoreSubVector(YdotRHS[i],rk->is_fast,&YdotRHS_fast[i]);
372: }
374: if (subrk_fast->subts_fast) {
375: subts_fast->ptime = t+k*h/rk->dtratio;
376: subts_fast->time_step = h/rk->dtratio;
377: TSStepRefine_RK_MultirateSplit(subts_fast);
378: }
379: /* update the fast components of the solution */
380: VecGetSubVector(ts->vec_sol,rk->is_fast,&Xfast);
381: VecISCopy(rk->X0,rk->is_fast,SCATTER_FORWARD,Xfast);
382: VecRestoreSubVector(ts->vec_sol,rk->is_fast,&Xfast);
383: }
384: return(0);
385: }
387: static PetscErrorCode TSStep_RK_MultirateSplit(TS ts)
388: {
389: TS_RK *rk = (TS_RK*)ts->data;
390: RKTableau tab = rk->tableau;
391: Vec *Y = rk->Y,*YdotRHS = rk->YdotRHS;
392: Vec *YdotRHS_fast = rk->YdotRHS_fast,*YdotRHS_slow = rk->YdotRHS_slow;
393: Vec Yslow,Yfast; /* subvectors store the stges of slow components and fast components respectively */
394: const PetscInt s = tab->s;
395: const PetscReal *A = tab->A,*c = tab->c;
396: PetscScalar *w = rk->work;
397: PetscInt i,j;
398: PetscReal next_time_step = ts->time_step,t = ts->ptime,h = ts->time_step;
399: PetscErrorCode ierr;
402: rk->status = TS_STEP_INCOMPLETE;
403: for (i=0; i<s; i++) {
404: VecGetSubVector(YdotRHS[i],rk->is_slow,&YdotRHS_slow[i]);
405: VecGetSubVector(YdotRHS[i],rk->is_fast,&YdotRHS_fast[i]);
406: }
407: VecCopy(ts->vec_sol,rk->X0);
408: /* propagate both slow and fast components using large time steps */
409: for (i=0; i<s; i++) {
410: rk->stage_time = t + h*c[i];
411: TSPreStage(ts,rk->stage_time);
412: VecCopy(ts->vec_sol,Y[i]);
413: VecGetSubVector(Y[i],rk->is_fast,&Yfast);
414: VecGetSubVector(Y[i],rk->is_slow,&Yslow);
415: for (j=0; j<i; j++) w[j] = h*A[i*s+j];
416: VecMAXPY(Yfast,i,w,YdotRHS_fast);
417: VecMAXPY(Yslow,i,w,YdotRHS_slow);
418: VecRestoreSubVector(Y[i],rk->is_fast,&Yfast);
419: VecRestoreSubVector(Y[i],rk->is_slow,&Yslow);
420: TSPostStage(ts,rk->stage_time,i,Y);
421: TSComputeRHSFunction(rk->subts_slow,t+h*c[i],Y[i],YdotRHS_slow[i]);
422: TSComputeRHSFunction(rk->subts_fast,t+h*c[i],Y[i],YdotRHS_fast[i]);
423: }
424: rk->slow = PETSC_TRUE;
425: /* update the slow components of the solution using slow time step */
426: TSEvaluateStep_RK_MultirateSplit(ts,tab->order,ts->vec_sol,NULL);
427: /* YdotRHS will be used for interpolation during refinement */
428: for (i=0; i<s; i++) {
429: VecRestoreSubVector(YdotRHS[i],rk->is_slow,&YdotRHS_slow[i]);
430: VecRestoreSubVector(YdotRHS[i],rk->is_fast,&YdotRHS_fast[i]);
431: }
433: TSStepRefine_RK_MultirateSplit(ts);
435: ts->ptime = t + ts->time_step;
436: ts->time_step = next_time_step;
437: rk->status = TS_STEP_COMPLETE;
438: return(0);
439: }
441: static PetscErrorCode TSSetUp_RK_MultirateSplit(TS ts)
442: {
443: TS_RK *rk = (TS_RK*)ts->data,*nextlevelrk,*currentlevelrk;
444: TS nextlevelts;
445: Vec X0;
449: TSRHSSplitGetIS(ts,"slow",&rk->is_slow);
450: TSRHSSplitGetIS(ts,"fast",&rk->is_fast);
451: 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");
452: TSRHSSplitGetSubTS(ts,"slow",&rk->subts_slow);
453: TSRHSSplitGetSubTS(ts,"fast",&rk->subts_fast);
454: 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");
456: VecDuplicate(ts->vec_sol,&X0);
457: /* The TS at each level share the same tableau, work array, solution vector, stage values and stage derivatives */
458: currentlevelrk = rk;
459: while (currentlevelrk->subts_fast) {
460: PetscMalloc1(rk->tableau->s,¤tlevelrk->YdotRHS_fast);
461: PetscMalloc1(rk->tableau->s,¤tlevelrk->YdotRHS_slow);
462: PetscObjectReference((PetscObject)X0);
463: currentlevelrk->X0 = X0;
464: currentlevelrk->ts_root = ts;
466: /* set up the ts for the slow part */
467: nextlevelts = currentlevelrk->subts_slow;
468: PetscNewLog(nextlevelts,&nextlevelrk);
469: nextlevelrk->tableau = rk->tableau;
470: nextlevelrk->work = rk->work;
471: nextlevelrk->Y = rk->Y;
472: nextlevelrk->YdotRHS = rk->YdotRHS;
473: nextlevelts->data = (void*)nextlevelrk;
474: TSCopyDM(ts,nextlevelts);
475: TSSetSolution(nextlevelts,ts->vec_sol);
477: /* set up the ts for the fast part */
478: nextlevelts = currentlevelrk->subts_fast;
479: PetscNewLog(nextlevelts,&nextlevelrk);
480: nextlevelrk->tableau = rk->tableau;
481: nextlevelrk->work = rk->work;
482: nextlevelrk->Y = rk->Y;
483: nextlevelrk->YdotRHS = rk->YdotRHS;
484: nextlevelrk->dtratio = rk->dtratio;
485: TSRHSSplitGetIS(nextlevelts,"slow",&nextlevelrk->is_slow);
486: TSRHSSplitGetSubTS(nextlevelts,"slow",&nextlevelrk->subts_slow);
487: TSRHSSplitGetIS(nextlevelts,"fast",&nextlevelrk->is_fast);
488: TSRHSSplitGetSubTS(nextlevelts,"fast",&nextlevelrk->subts_fast);
489: nextlevelts->data = (void*)nextlevelrk;
490: TSCopyDM(ts,nextlevelts);
491: TSSetSolution(nextlevelts,ts->vec_sol);
493: currentlevelrk = nextlevelrk;
494: }
495: VecDestroy(&X0);
497: ts->ops->step = TSStep_RK_MultirateSplit;
498: ts->ops->evaluatestep = TSEvaluateStep_RK_MultirateSplit;
499: ts->ops->interpolate = TSInterpolate_RK_MultirateSplit;
500: return(0);
501: }
503: PetscErrorCode TSRKSetMultirate_RK(TS ts,PetscBool use_multirate)
504: {
505: TS_RK *rk = (TS_RK*)ts->data;
510: rk->use_multirate = use_multirate;
511: if (use_multirate) {
512: rk->dtratio = 2;
513: PetscObjectComposeFunction((PetscObject)ts,"TSSetUp_RK_MultirateSplit_C",TSSetUp_RK_MultirateSplit);
514: PetscObjectComposeFunction((PetscObject)ts,"TSReset_RK_MultirateSplit_C",TSReset_RK_MultirateSplit);
515: PetscObjectComposeFunction((PetscObject)ts,"TSSetUp_RK_MultirateNonsplit_C",TSSetUp_RK_MultirateNonsplit);
516: PetscObjectComposeFunction((PetscObject)ts,"TSReset_RK_MultirateNonsplit_C",TSReset_RK_MultirateNonsplit);
517: } else {
518: rk->dtratio = 0;
519: PetscObjectComposeFunction((PetscObject)ts,"TSSetUp_RK_MultirateSplit_C",NULL);
520: PetscObjectComposeFunction((PetscObject)ts,"TSReset_RK_MultirateSplit_C",NULL);
521: PetscObjectComposeFunction((PetscObject)ts,"TSSetUp_RK_MultirateNonsplit_C",NULL);
522: PetscObjectComposeFunction((PetscObject)ts,"TSReset_RK_MultirateNonsplit_C",NULL);
523: }
524: return(0);
525: }
527: PetscErrorCode TSRKGetMultirate_RK(TS ts,PetscBool *use_multirate)
528: {
529: TS_RK *rk = (TS_RK*)ts->data;
533: *use_multirate = rk->use_multirate;
534: return(0);
535: }