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: PetscFunctionBegin;
25: PetscCall(VecDestroy(&rk->X0));
26: PetscCall(VecDestroyVecs(tab->s, &rk->YdotRHS_slow));
27: PetscFunctionReturn(PETSC_SUCCESS);
28: }
30: static PetscErrorCode TSInterpolate_RK_MultirateNonsplit(TS ts, PetscReal itime, Vec X)
31: {
32: TS_RK *rk = (TS_RK *)ts->data;
33: PetscInt s = rk->tableau->s, p = rk->tableau->p, i, j;
34: PetscReal h = ts->time_step;
35: PetscReal tt, t;
36: PetscScalar *b;
37: const PetscReal *B = rk->tableau->binterp;
39: PetscFunctionBegin;
40: PetscCheck(B, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "TSRK %s does not have an interpolation formula", rk->tableau->name);
41: t = (itime - rk->ptime) / h;
42: PetscCall(PetscMalloc1(s, &b));
43: for (i = 0; i < s; i++) b[i] = 0;
44: for (j = 0, tt = t; j < p; j++, tt *= t) {
45: for (i = 0; i < s; i++) b[i] += h * B[i * p + j] * tt;
46: }
47: PetscCall(VecCopy(rk->X0, X));
48: PetscCall(VecMAXPY(X, s, b, rk->YdotRHS_slow));
49: PetscCall(PetscFree(b));
50: PetscFunctionReturn(PETSC_SUCCESS);
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: PetscFunctionBegin;
67: PetscCall(VecDuplicate(ts->vec_sol, &vec_fast));
68: previousts = rk->subts_current;
69: PetscCall(TSRHSSplitGetSubTS(rk->subts_current, "fast", &subts));
70: PetscCall(TSRHSSplitGetSubTS(subts, "fast", &subts));
71: for (k = 0; k < rk->dtratio; k++) {
72: for (i = 0; i < s; i++) {
73: PetscCall(TSInterpolate_RK_MultirateNonsplit(ts, t + k * h / rk->dtratio + h / rk->dtratio * c[i], Y[i]));
74: for (j = 0; j < i; j++) w[j] = h / rk->dtratio * A[i * s + j];
75: /* 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 */
76: PetscCall(VecCopy(ts->vec_sol, vec_fast));
77: PetscCall(VecMAXPY(vec_fast, i, w, YdotRHS));
78: /* update the fast components in the stage value */
79: PetscCall(VecGetSubVector(vec_fast, rk->is_fast, &subvec_fast));
80: PetscCall(VecISCopy(Y[i], rk->is_fast, SCATTER_FORWARD, subvec_fast));
81: PetscCall(VecRestoreSubVector(vec_fast, rk->is_fast, &subvec_fast));
82: /* compute the stage RHS */
83: PetscCall(TSComputeRHSFunction(ts, t + k * h / rk->dtratio + h / rk->dtratio * c[i], Y[i], YdotRHS[i]));
84: }
85: PetscCall(VecCopy(ts->vec_sol, vec_fast));
86: PetscCall(TSEvaluateStep(ts, tab->order, vec_fast, NULL));
87: /* update the fast components in the solution */
88: PetscCall(VecGetSubVector(vec_fast, rk->is_fast, &subvec_fast));
89: PetscCall(VecISCopy(ts->vec_sol, rk->is_fast, SCATTER_FORWARD, subvec_fast));
90: PetscCall(VecRestoreSubVector(vec_fast, rk->is_fast, &subvec_fast));
92: if (subts) {
93: Vec *YdotRHS_copy;
94: PetscCall(VecDuplicateVecs(ts->vec_sol, s, &YdotRHS_copy));
95: rk->subts_current = rk->subts_fast;
96: ts->ptime = t + k * h / rk->dtratio;
97: ts->time_step = h / rk->dtratio;
98: PetscCall(TSRHSSplitGetIS(rk->subts_current, "fast", &rk->is_fast));
99: for (i = 0; i < s; i++) {
100: PetscCall(VecCopy(rk->YdotRHS_slow[i], YdotRHS_copy[i]));
101: PetscCall(VecCopy(YdotRHS[i], rk->YdotRHS_slow[i]));
102: }
104: PetscCall(TSStepRefine_RK_MultirateNonsplit(ts));
106: rk->subts_current = previousts;
107: ts->ptime = t;
108: ts->time_step = h;
109: PetscCall(TSRHSSplitGetIS(previousts, "fast", &rk->is_fast));
110: for (i = 0; i < s; i++) PetscCall(VecCopy(YdotRHS_copy[i], rk->YdotRHS_slow[i]));
111: PetscCall(VecDestroyVecs(s, &YdotRHS_copy));
112: }
113: }
114: PetscCall(VecDestroy(&vec_fast));
115: PetscFunctionReturn(PETSC_SUCCESS);
116: }
118: static PetscErrorCode TSStep_RK_MultirateNonsplit(TS ts)
119: {
120: TS_RK *rk = (TS_RK *)ts->data;
121: RKTableau tab = rk->tableau;
122: Vec *Y = rk->Y, *YdotRHS = rk->YdotRHS, *YdotRHS_slow = rk->YdotRHS_slow;
123: Vec stage_slow, sol_slow; /* vectors store the slow components */
124: Vec subvec_slow; /* sub vector to store the slow components */
125: IS is_slow = rk->is_slow;
126: const PetscInt s = tab->s;
127: const PetscReal *A = tab->A, *c = tab->c;
128: PetscScalar *w = rk->work;
129: PetscInt i, j, dtratio = rk->dtratio;
130: PetscReal next_time_step = ts->time_step, t = ts->ptime, h = ts->time_step;
132: PetscFunctionBegin;
133: rk->status = TS_STEP_INCOMPLETE;
134: PetscCall(VecDuplicate(ts->vec_sol, &stage_slow));
135: PetscCall(VecDuplicate(ts->vec_sol, &sol_slow));
136: PetscCall(VecCopy(ts->vec_sol, rk->X0));
137: for (i = 0; i < s; i++) {
138: rk->stage_time = t + h * c[i];
139: PetscCall(TSPreStage(ts, rk->stage_time));
140: PetscCall(VecCopy(ts->vec_sol, Y[i]));
141: for (j = 0; j < i; j++) w[j] = h * A[i * s + j];
142: PetscCall(VecMAXPY(Y[i], i, w, YdotRHS_slow));
143: PetscCall(TSPostStage(ts, rk->stage_time, i, Y));
144: /* compute the stage RHS */
145: PetscCall(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: PetscCall(TSEvaluateStep(ts, tab->order, sol_slow, NULL));
151: rk->dtratio = dtratio;
152: rk->YdotRHS = YdotRHS;
153: /* update the slow components in the solution */
154: PetscCall(VecGetSubVector(sol_slow, is_slow, &subvec_slow));
155: PetscCall(VecISCopy(ts->vec_sol, is_slow, SCATTER_FORWARD, subvec_slow));
156: PetscCall(VecRestoreSubVector(sol_slow, is_slow, &subvec_slow));
158: rk->subts_current = ts;
159: rk->ptime = t;
160: rk->time_step = h;
161: PetscCall(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: PetscCall(VecDestroy(&stage_slow));
169: PetscCall(VecDestroy(&sol_slow));
170: PetscFunctionReturn(PETSC_SUCCESS);
171: }
173: static PetscErrorCode TSSetUp_RK_MultirateNonsplit(TS ts)
174: {
175: TS_RK *rk = (TS_RK *)ts->data;
176: RKTableau tab = rk->tableau;
178: PetscFunctionBegin;
179: PetscCall(TSRHSSplitGetIS(ts, "slow", &rk->is_slow));
180: PetscCall(TSRHSSplitGetIS(ts, "fast", &rk->is_fast));
181: PetscCheck(rk->is_slow && rk->is_fast, 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");
182: PetscCall(TSRHSSplitGetSubTS(ts, "slow", &rk->subts_slow));
183: PetscCall(TSRHSSplitGetSubTS(ts, "fast", &rk->subts_fast));
184: PetscCheck(rk->subts_slow && rk->subts_fast, 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");
185: PetscCall(VecDuplicate(ts->vec_sol, &rk->X0));
186: PetscCall(VecDuplicateVecs(ts->vec_sol, tab->s, &rk->YdotRHS_slow));
187: rk->subts_current = rk->subts_fast;
189: ts->ops->step = TSStep_RK_MultirateNonsplit;
190: ts->ops->interpolate = TSInterpolate_RK_MultirateNonsplit;
191: PetscFunctionReturn(PETSC_SUCCESS);
192: }
194: /*
195: Copy DM from tssrc to tsdest, while keeping the original DMTS and DMSNES in tsdest.
196: */
197: static PetscErrorCode TSCopyDM(TS tssrc, TS tsdest)
198: {
199: DM newdm, dmsrc, dmdest;
201: PetscFunctionBegin;
202: PetscCall(TSGetDM(tssrc, &dmsrc));
203: PetscCall(DMClone(dmsrc, &newdm));
204: PetscCall(TSGetDM(tsdest, &dmdest));
205: PetscCall(DMCopyDMTS(dmdest, newdm));
206: PetscCall(DMCopyDMSNES(dmdest, newdm));
207: PetscCall(TSSetDM(tsdest, newdm));
208: PetscCall(DMDestroy(&newdm));
209: PetscFunctionReturn(PETSC_SUCCESS);
210: }
212: static PetscErrorCode TSReset_RK_MultirateSplit(TS ts)
213: {
214: TS_RK *rk = (TS_RK *)ts->data;
216: PetscFunctionBegin;
217: if (rk->subts_slow) {
218: PetscCall(PetscFree(rk->subts_slow->data));
219: rk->subts_slow = NULL;
220: }
221: if (rk->subts_fast) {
222: PetscCall(PetscFree(rk->YdotRHS_fast));
223: PetscCall(PetscFree(rk->YdotRHS_slow));
224: PetscCall(VecDestroy(&rk->X0));
225: PetscCall(TSReset_RK_MultirateSplit(rk->subts_fast));
226: PetscCall(PetscFree(rk->subts_fast->data));
227: rk->subts_fast = NULL;
228: }
229: PetscFunctionReturn(PETSC_SUCCESS);
230: }
232: static PetscErrorCode TSInterpolate_RK_MultirateSplit(TS ts, PetscReal itime, Vec X)
233: {
234: TS_RK *rk = (TS_RK *)ts->data;
235: Vec Xslow;
236: PetscInt s = rk->tableau->s, p = rk->tableau->p, i, j;
237: PetscReal h;
238: PetscReal tt, t;
239: PetscScalar *b;
240: const PetscReal *B = rk->tableau->binterp;
242: PetscFunctionBegin;
243: PetscCheck(B, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "TSRK %s does not have an interpolation formula", rk->tableau->name);
245: switch (rk->status) {
246: case TS_STEP_INCOMPLETE:
247: case TS_STEP_PENDING:
248: h = ts->time_step;
249: t = (itime - ts->ptime) / h;
250: break;
251: case TS_STEP_COMPLETE:
252: h = ts->ptime - ts->ptime_prev;
253: t = (itime - ts->ptime) / h + 1; /* In the interval [0,1] */
254: break;
255: default:
256: SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_PLIB, "Invalid TSStepStatus");
257: }
258: PetscCall(PetscMalloc1(s, &b));
259: for (i = 0; i < s; i++) b[i] = 0;
260: for (j = 0, tt = t; j < p; j++, tt *= t) {
261: for (i = 0; i < s; i++) b[i] += h * B[i * p + j] * tt;
262: }
263: for (i = 0; i < s; i++) PetscCall(VecGetSubVector(rk->YdotRHS[i], rk->is_slow, &rk->YdotRHS_slow[i]));
264: PetscCall(VecGetSubVector(X, rk->is_slow, &Xslow));
265: PetscCall(VecISCopy(rk->X0, rk->is_slow, SCATTER_REVERSE, Xslow));
266: PetscCall(VecMAXPY(Xslow, s, b, rk->YdotRHS_slow));
267: PetscCall(VecRestoreSubVector(X, rk->is_slow, &Xslow));
268: for (i = 0; i < s; i++) PetscCall(VecRestoreSubVector(rk->YdotRHS[i], rk->is_slow, &rk->YdotRHS_slow[i]));
269: PetscCall(PetscFree(b));
270: PetscFunctionReturn(PETSC_SUCCESS);
271: }
273: /*
274: This is for partitioned RHS multirate RK method
275: The step completion formula is
277: x1 = x0 + h b^T YdotRHS
279: */
280: static PetscErrorCode TSEvaluateStep_RK_MultirateSplit(TS ts, PetscInt order, Vec X, PetscBool *done)
281: {
282: TS_RK *rk = (TS_RK *)ts->data;
283: RKTableau tab = rk->tableau;
284: Vec Xslow, Xfast; /* subvectors of X which store slow components and fast components respectively */
285: PetscScalar *w = rk->work;
286: PetscReal h = ts->time_step;
287: PetscInt s = tab->s, j;
289: PetscFunctionBegin;
290: PetscCall(VecCopy(ts->vec_sol, X));
291: if (rk->slow) {
292: for (j = 0; j < s; j++) w[j] = h * tab->b[j];
293: PetscCall(VecGetSubVector(ts->vec_sol, rk->is_slow, &Xslow));
294: PetscCall(VecMAXPY(Xslow, s, w, rk->YdotRHS_slow));
295: PetscCall(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: PetscCall(VecGetSubVector(X, rk->is_fast, &Xfast));
299: PetscCall(VecMAXPY(Xfast, s, w, rk->YdotRHS_fast));
300: PetscCall(VecRestoreSubVector(X, rk->is_fast, &Xfast));
301: }
302: PetscFunctionReturn(PETSC_SUCCESS);
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: PetscFunctionBegin;
321: for (k = 0; k < rk->dtratio; k++) {
322: PetscCall(VecGetSubVector(ts->vec_sol, rk->is_fast, &Xfast));
323: for (i = 0; i < s; i++) PetscCall(VecGetSubVector(YdotRHS[i], rk->is_fast, &YdotRHS_fast[i]));
324: /* propagate fast component using small time steps */
325: for (i = 0; i < s; i++) {
326: /* stage value for slow components */
327: PetscCall(TSInterpolate_RK_MultirateSplit(rk->ts_root, t + k * h / rk->dtratio + h / rk->dtratio * c[i], Y[i]));
328: currentlevelts = rk->ts_root;
329: while (currentlevelts != ts) { /* all the slow parts need to be interpolated separated */
330: currentlevelts = ((TS_RK *)currentlevelts->data)->subts_fast;
331: PetscCall(TSInterpolate_RK_MultirateSplit(currentlevelts, t + k * h / rk->dtratio + h / rk->dtratio * c[i], Y[i]));
332: }
333: for (j = 0; j < i; j++) w[j] = h / rk->dtratio * A[i * s + j];
334: subrk_fast->stage_time = t + h / rk->dtratio * c[i];
335: PetscCall(TSPreStage(subts_fast, subrk_fast->stage_time));
336: /* stage value for fast components */
337: PetscCall(VecGetSubVector(Y[i], rk->is_fast, &Yfast));
338: PetscCall(VecCopy(Xfast, Yfast));
339: PetscCall(VecMAXPY(Yfast, i, w, YdotRHS_fast));
340: PetscCall(VecRestoreSubVector(Y[i], rk->is_fast, &Yfast));
341: PetscCall(TSPostStage(subts_fast, subrk_fast->stage_time, i, Y));
342: /* compute the stage RHS for fast components */
343: PetscCall(TSComputeRHSFunction(subts_fast, t + k * h * rk->dtratio + h / rk->dtratio * c[i], Y[i], YdotRHS_fast[i]));
344: }
345: PetscCall(VecRestoreSubVector(ts->vec_sol, rk->is_fast, &Xfast));
346: /* update the value of fast components using fast time step */
347: rk->slow = PETSC_FALSE;
348: PetscCall(TSEvaluateStep_RK_MultirateSplit(ts, tab->order, ts->vec_sol, NULL));
349: for (i = 0; i < s; i++) PetscCall(VecRestoreSubVector(YdotRHS[i], rk->is_fast, &YdotRHS_fast[i]));
351: if (subrk_fast->subts_fast) {
352: subts_fast->ptime = t + k * h / rk->dtratio;
353: subts_fast->time_step = h / rk->dtratio;
354: PetscCall(TSStepRefine_RK_MultirateSplit(subts_fast));
355: }
356: /* update the fast components of the solution */
357: PetscCall(VecGetSubVector(ts->vec_sol, rk->is_fast, &Xfast));
358: PetscCall(VecISCopy(rk->X0, rk->is_fast, SCATTER_FORWARD, Xfast));
359: PetscCall(VecRestoreSubVector(ts->vec_sol, rk->is_fast, &Xfast));
360: }
361: PetscFunctionReturn(PETSC_SUCCESS);
362: }
364: static PetscErrorCode TSStep_RK_MultirateSplit(TS ts)
365: {
366: TS_RK *rk = (TS_RK *)ts->data;
367: RKTableau tab = rk->tableau;
368: Vec *Y = rk->Y, *YdotRHS = rk->YdotRHS;
369: Vec *YdotRHS_fast = rk->YdotRHS_fast, *YdotRHS_slow = rk->YdotRHS_slow;
370: Vec Yslow, Yfast; /* subvectors store the stges of slow components and fast components respectively */
371: const PetscInt s = tab->s;
372: const PetscReal *A = tab->A, *c = tab->c;
373: PetscScalar *w = rk->work;
374: PetscInt i, j;
375: PetscReal next_time_step = ts->time_step, t = ts->ptime, h = ts->time_step;
377: PetscFunctionBegin;
378: rk->status = TS_STEP_INCOMPLETE;
379: for (i = 0; i < s; i++) {
380: PetscCall(VecGetSubVector(YdotRHS[i], rk->is_slow, &YdotRHS_slow[i]));
381: PetscCall(VecGetSubVector(YdotRHS[i], rk->is_fast, &YdotRHS_fast[i]));
382: }
383: PetscCall(VecCopy(ts->vec_sol, rk->X0));
384: /* propagate both slow and fast components using large time steps */
385: for (i = 0; i < s; i++) {
386: rk->stage_time = t + h * c[i];
387: PetscCall(TSPreStage(ts, rk->stage_time));
388: PetscCall(VecCopy(ts->vec_sol, Y[i]));
389: PetscCall(VecGetSubVector(Y[i], rk->is_fast, &Yfast));
390: PetscCall(VecGetSubVector(Y[i], rk->is_slow, &Yslow));
391: for (j = 0; j < i; j++) w[j] = h * A[i * s + j];
392: PetscCall(VecMAXPY(Yfast, i, w, YdotRHS_fast));
393: PetscCall(VecMAXPY(Yslow, i, w, YdotRHS_slow));
394: PetscCall(VecRestoreSubVector(Y[i], rk->is_fast, &Yfast));
395: PetscCall(VecRestoreSubVector(Y[i], rk->is_slow, &Yslow));
396: PetscCall(TSPostStage(ts, rk->stage_time, i, Y));
397: PetscCall(TSComputeRHSFunction(rk->subts_slow, t + h * c[i], Y[i], YdotRHS_slow[i]));
398: PetscCall(TSComputeRHSFunction(rk->subts_fast, t + h * c[i], Y[i], YdotRHS_fast[i]));
399: }
400: rk->slow = PETSC_TRUE;
401: /* update the slow components of the solution using slow time step */
402: PetscCall(TSEvaluateStep_RK_MultirateSplit(ts, tab->order, ts->vec_sol, NULL));
403: /* YdotRHS will be used for interpolation during refinement */
404: for (i = 0; i < s; i++) {
405: PetscCall(VecRestoreSubVector(YdotRHS[i], rk->is_slow, &YdotRHS_slow[i]));
406: PetscCall(VecRestoreSubVector(YdotRHS[i], rk->is_fast, &YdotRHS_fast[i]));
407: }
409: PetscCall(TSStepRefine_RK_MultirateSplit(ts));
411: ts->ptime = t + ts->time_step;
412: ts->time_step = next_time_step;
413: rk->status = TS_STEP_COMPLETE;
414: PetscFunctionReturn(PETSC_SUCCESS);
415: }
417: static PetscErrorCode TSSetUp_RK_MultirateSplit(TS ts)
418: {
419: TS_RK *rk = (TS_RK *)ts->data, *nextlevelrk, *currentlevelrk;
420: TS nextlevelts;
421: Vec X0;
423: PetscFunctionBegin;
424: PetscCall(TSRHSSplitGetIS(ts, "slow", &rk->is_slow));
425: PetscCall(TSRHSSplitGetIS(ts, "fast", &rk->is_fast));
426: PetscCheck(rk->is_slow && rk->is_fast, 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");
427: PetscCall(TSRHSSplitGetSubTS(ts, "slow", &rk->subts_slow));
428: PetscCall(TSRHSSplitGetSubTS(ts, "fast", &rk->subts_fast));
429: PetscCheck(rk->subts_slow && rk->subts_fast, 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");
431: PetscCall(VecDuplicate(ts->vec_sol, &X0));
432: /* The TS at each level share the same tableau, work array, solution vector, stage values and stage derivatives */
433: currentlevelrk = rk;
434: while (currentlevelrk->subts_fast) {
435: PetscCall(PetscMalloc1(rk->tableau->s, ¤tlevelrk->YdotRHS_fast));
436: PetscCall(PetscMalloc1(rk->tableau->s, ¤tlevelrk->YdotRHS_slow));
437: PetscCall(PetscObjectReference((PetscObject)X0));
438: currentlevelrk->X0 = X0;
439: currentlevelrk->ts_root = ts;
441: /* set up the ts for the slow part */
442: nextlevelts = currentlevelrk->subts_slow;
443: PetscCall(PetscNew(&nextlevelrk));
444: nextlevelrk->tableau = rk->tableau;
445: nextlevelrk->work = rk->work;
446: nextlevelrk->Y = rk->Y;
447: nextlevelrk->YdotRHS = rk->YdotRHS;
448: nextlevelts->data = (void *)nextlevelrk;
449: PetscCall(TSCopyDM(ts, nextlevelts));
450: PetscCall(TSSetSolution(nextlevelts, ts->vec_sol));
452: /* set up the ts for the fast part */
453: nextlevelts = currentlevelrk->subts_fast;
454: PetscCall(PetscNew(&nextlevelrk));
455: nextlevelrk->tableau = rk->tableau;
456: nextlevelrk->work = rk->work;
457: nextlevelrk->Y = rk->Y;
458: nextlevelrk->YdotRHS = rk->YdotRHS;
459: nextlevelrk->dtratio = rk->dtratio;
460: PetscCall(TSRHSSplitGetIS(nextlevelts, "slow", &nextlevelrk->is_slow));
461: PetscCall(TSRHSSplitGetSubTS(nextlevelts, "slow", &nextlevelrk->subts_slow));
462: PetscCall(TSRHSSplitGetIS(nextlevelts, "fast", &nextlevelrk->is_fast));
463: PetscCall(TSRHSSplitGetSubTS(nextlevelts, "fast", &nextlevelrk->subts_fast));
464: nextlevelts->data = (void *)nextlevelrk;
465: PetscCall(TSCopyDM(ts, nextlevelts));
466: PetscCall(TSSetSolution(nextlevelts, ts->vec_sol));
468: currentlevelrk = nextlevelrk;
469: }
470: PetscCall(VecDestroy(&X0));
472: ts->ops->step = TSStep_RK_MultirateSplit;
473: ts->ops->evaluatestep = TSEvaluateStep_RK_MultirateSplit;
474: ts->ops->interpolate = TSInterpolate_RK_MultirateSplit;
475: PetscFunctionReturn(PETSC_SUCCESS);
476: }
478: PetscErrorCode TSRKSetMultirate_RK(TS ts, PetscBool use_multirate)
479: {
480: TS_RK *rk = (TS_RK *)ts->data;
482: PetscFunctionBegin;
484: rk->use_multirate = use_multirate;
485: if (use_multirate) {
486: rk->dtratio = 2;
487: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSSetUp_RK_MultirateSplit_C", TSSetUp_RK_MultirateSplit));
488: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSReset_RK_MultirateSplit_C", TSReset_RK_MultirateSplit));
489: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSSetUp_RK_MultirateNonsplit_C", TSSetUp_RK_MultirateNonsplit));
490: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSReset_RK_MultirateNonsplit_C", TSReset_RK_MultirateNonsplit));
491: } else {
492: rk->dtratio = 0;
493: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSSetUp_RK_MultirateSplit_C", NULL));
494: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSReset_RK_MultirateSplit_C", NULL));
495: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSSetUp_RK_MultirateNonsplit_C", NULL));
496: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSReset_RK_MultirateNonsplit_C", NULL));
497: }
498: PetscFunctionReturn(PETSC_SUCCESS);
499: }
501: PetscErrorCode TSRKGetMultirate_RK(TS ts, PetscBool *use_multirate)
502: {
503: TS_RK *rk = (TS_RK *)ts->data;
505: PetscFunctionBegin;
507: *use_multirate = rk->use_multirate;
508: PetscFunctionReturn(PETSC_SUCCESS);
509: }