Actual source code: mrk.c
petsc-3.14.6 2021-03-30
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,¤tlevelrk->YdotRHS_fast);
459: PetscMalloc1(rk->tableau->s,¤tlevelrk->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: }