1 /*
2 Code for time stepping with the Multirate Partitioned Runge-Kutta method
3
4 Notes:
5 1) The general system is written as
6 Udot = F(t,U)
7 if one does not split the RHS function, but gives 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)
11 for component-wise partitioned system,
12 users should split the RHS function themselves and also provide the indexes for both slow and fast components.
13 3) To correct The confusing terminology in the paper, we use 'slow method', 'slow buffer method' and 'fast method' to denote the methods applied to 'slow region', 'slow buffer region' and 'fast region' respectively. The 'slow method' in the original paper actually means the 'slow buffer method'.
14 4) Why does the buffer region have to be inside the slow region? The buffer region is treated with a slow method essentially. Applying the slow method to a region with a fast characteristic time scale is apparently not a good choice.
15
16 Reference:
17 Emil M. Constantinescu, Adrian Sandu, Multirate Timestepping Methods for Hyperbolic Conservation Laws, Journal of Scientific Computing 2007
18 */
19
20 #include <petsc/private/tsimpl.h> /*I "petscts.h" I*/
21 #include <petscdm.h>
22
23 static TSMPRKType TSMPRKDefault = TSMPRK2A22;
24 static PetscBool TSMPRKRegisterAllCalled;
25 static PetscBool TSMPRKPackageInitialized;
26
27 typedef struct _MPRKTableau *MPRKTableau;
28 struct _MPRKTableau {
29 char *name;
30 PetscInt order; /* Classical approximation order of the method i */
31 PetscInt sbase; /* Number of stages in the base method*/
32 PetscInt s; /* Number of stages */
33 PetscInt np; /* Number of partitions */
34 PetscReal *Af, *bf, *cf; /* Tableau for fast components */
35 PetscReal *Amb, *bmb, *cmb; /* Tableau for medium components */
36 PetscInt *rmb; /* Array of flags for repeated stages in medium method */
37 PetscReal *Asb, *bsb, *csb; /* Tableau for slow components */
38 PetscInt *rsb; /* Array of flags for repeated staged in slow method*/
39 };
40 typedef struct _MPRKTableauLink *MPRKTableauLink;
41 struct _MPRKTableauLink {
42 struct _MPRKTableau tab;
43 MPRKTableauLink next;
44 };
45 static MPRKTableauLink MPRKTableauList;
46
47 typedef struct {
48 MPRKTableau tableau;
49 Vec *Y; /* States computed during the step */
50 Vec *YdotRHS;
51 Vec *YdotRHS_slow; /* Function evaluations by slow tableau for slow components */
52 Vec *YdotRHS_slowbuffer; /* Function evaluations by medium tableau for slow components */
53 Vec *YdotRHS_medium; /* Function evaluations by medium tableau for medium components*/
54 Vec *YdotRHS_mediumbuffer; /* Function evaluations by fast tableau for medium components */
55 Vec *YdotRHS_fast; /* Function evaluations by fast tableau for fast components */
56 PetscScalar *work_slow; /* Scalar work_slow by slow tableau */
57 PetscScalar *work_slowbuffer; /* Scalar work_slow by slow tableau */
58 PetscScalar *work_medium; /* Scalar work_slow by medium tableau */
59 PetscScalar *work_mediumbuffer; /* Scalar work_mediumbuffer by fast tableau */
60 PetscScalar *work_fast; /* Scalar work_fast by fast tableau */
61 PetscReal stage_time;
62 TSStepStatus status;
63 PetscReal ptime;
64 PetscReal time_step;
65 IS is_slow, is_slowbuffer, is_medium, is_mediumbuffer, is_fast;
66 TS subts_slow, subts_slowbuffer, subts_medium, subts_mediumbuffer, subts_fast;
67 } TS_MPRK;
68
TSMPRKGenerateTableau2(PetscInt ratio,PetscInt s,const PetscReal Abase[],const PetscReal bbase[],PetscReal A1[],PetscReal b1[],PetscReal A2[],PetscReal b2[])69 static PetscErrorCode TSMPRKGenerateTableau2(PetscInt ratio, PetscInt s, const PetscReal Abase[], const PetscReal bbase[], PetscReal A1[], PetscReal b1[], PetscReal A2[], PetscReal b2[])
70 {
71 PetscInt i, j, k, l;
72
73 PetscFunctionBegin;
74 for (k = 0; k < ratio; k++) {
75 /* diagonal blocks */
76 for (i = 0; i < s; i++)
77 for (j = 0; j < s; j++) {
78 A1[(k * s + i) * ratio * s + k * s + j] = Abase[i * s + j];
79 A2[(k * s + i) * ratio * s + k * s + j] = Abase[i * s + j] / ratio;
80 }
81 /* off-diagonal blocks */
82 for (l = 0; l < k; l++)
83 for (i = 0; i < s; i++)
84 for (j = 0; j < s; j++) A2[(k * s + i) * ratio * s + l * s + j] = bbase[j] / ratio;
85 for (j = 0; j < s; j++) {
86 b1[k * s + j] = bbase[j] / ratio;
87 b2[k * s + j] = bbase[j] / ratio;
88 }
89 }
90 PetscFunctionReturn(PETSC_SUCCESS);
91 }
92
TSMPRKGenerateTableau3(PetscInt ratio,PetscInt s,const PetscReal Abase[],const PetscReal bbase[],PetscReal A1[],PetscReal b1[],PetscReal A2[],PetscReal b2[],PetscReal A3[],PetscReal b3[])93 static PetscErrorCode TSMPRKGenerateTableau3(PetscInt ratio, PetscInt s, const PetscReal Abase[], const PetscReal bbase[], PetscReal A1[], PetscReal b1[], PetscReal A2[], PetscReal b2[], PetscReal A3[], PetscReal b3[])
94 {
95 PetscInt i, j, k, l, m, n;
96
97 PetscFunctionBegin;
98 for (k = 0; k < ratio; k++) { /* diagonal blocks of size ratio*s by ratio*s */
99 for (l = 0; l < ratio; l++) /* diagonal sub-blocks of size s by s */
100 for (i = 0; i < s; i++)
101 for (j = 0; j < s; j++) {
102 A1[((k * ratio + l) * s + i) * ratio * ratio * s + (k * ratio + l) * s + j] = Abase[i * s + j];
103 A2[((k * ratio + l) * s + i) * ratio * ratio * s + (k * ratio + l) * s + j] = Abase[i * s + j] / ratio;
104 A3[((k * ratio + l) * s + i) * ratio * ratio * s + (k * ratio + l) * s + j] = Abase[i * s + j] / ratio / ratio;
105 }
106 for (l = 0; l < k; l++) /* off-diagonal blocks of size ratio*s by ratio*s */
107 for (m = 0; m < ratio; m++)
108 for (n = 0; n < ratio; n++)
109 for (i = 0; i < s; i++)
110 for (j = 0; j < s; j++) {
111 A2[((k * ratio + m) * s + i) * ratio * ratio * s + (l * ratio + n) * s + j] = bbase[j] / ratio / ratio;
112 A3[((k * ratio + m) * s + i) * ratio * ratio * s + (l * ratio + n) * s + j] = bbase[j] / ratio / ratio;
113 }
114 for (m = 0; m < ratio; m++)
115 for (n = 0; n < m; n++) /* off-diagonal sub-blocks of size s by s in the diagonal blocks */
116 for (i = 0; i < s; i++)
117 for (j = 0; j < s; j++) A3[((k * ratio + m) * s + i) * ratio * ratio * s + (k * ratio + n) * s + j] = bbase[j] / ratio / ratio;
118 for (n = 0; n < ratio; n++)
119 for (j = 0; j < s; j++) {
120 b1[(k * ratio + n) * s + j] = bbase[j] / ratio / ratio;
121 b2[(k * ratio + n) * s + j] = bbase[j] / ratio / ratio;
122 b3[(k * ratio + n) * s + j] = bbase[j] / ratio / ratio;
123 }
124 }
125 PetscFunctionReturn(PETSC_SUCCESS);
126 }
127
128 /*MC
129 TSMPRK2A22 - Second Order Multirate Partitioned Runge Kutta scheme based on RK2A.
130
131 This method has four stages for slow and fast parts. The refinement factor of the stepsize is 2.
132 r = 2, np = 2
133
134 Options Database Key:
135 . -ts_mprk_type 2a22 - select this scheme
136
137 Level: advanced
138
139 .seealso: [](ch_ts), `TSMPRK`, `TSMPRKType`, `TSMPRKSetType()`
140 M*/
141 /*MC
142 TSMPRK2A23 - Second Order Multirate Partitioned Runge-Kutta scheme based on RK2A.
143
144 This method has eight stages for slow and medium and fast parts. The refinement factor of the stepsize is 2.
145 r = 2, np = 3
146
147 Options Database Key:
148 . -ts_mprk_type 2a23 - select this scheme
149
150 Level: advanced
151
152 .seealso: [](ch_ts), `TSMPRK`, `TSMPRKType`, `TSMPRKSetType()`
153 M*/
154 /*MC
155 TSMPRK2A32 - Second Order Multirate Partitioned Runge-Kutta scheme based on RK2A.
156
157 This method has four stages for slow and fast parts. The refinement factor of the stepsize is 3.
158 r = 3, np = 2
159
160 Options Database Key:
161 . -ts_mprk_type 2a32 - select this scheme
162
163 Level: advanced
164
165 .seealso: [](ch_ts), `TSMPRK`, `TSMPRKType`, `TSMPRKSetType()`
166 M*/
167 /*MC
168 TSMPRK2A33 - Second Order Multirate Partitioned Runge-Kutta scheme based on RK2A.
169
170 This method has eight stages for slow and medium and fast parts. The refinement factor of the stepsize is 3.
171 r = 3, np = 3
172
173 Options Database Key:
174 . -ts_mprk_type 2a33- select this scheme
175
176 Level: advanced
177
178 .seealso: [](ch_ts), `TSMPRK`, `TSMPRKType`, `TSMPRKSetType()`
179 M*/
180 /*MC
181 TSMPRK3P2M - Third Order Multirate Partitioned Runge-Kutta scheme.
182
183 This method has eight stages for both slow and fast parts.
184
185 Options Database Key:
186 . -ts_mprk_type pm3 - select this scheme
187
188 Level: advanced
189
190 .seealso: [](ch_ts), `TSMPRK`, `TSMPRKType`, `TSMPRKSetType()`
191 M*/
192 /*MC
193 TSMPRKP2 - Second Order Multirate Partitioned Runge-Kutta scheme.
194
195 This method has five stages for both slow and fast parts.
196
197 Options Database Key:
198 . -ts_mprk_type p2 - select this scheme
199
200 Level: advanced
201
202 .seealso: [](ch_ts), `TSMPRK`, `TSMPRKType`, `TSMPRKSetType()`
203 M*/
204 /*MC
205 TSMPRKP3 - Third Order Multirate Partitioned Runge-Kutta scheme.
206
207 This method has ten stages for both slow and fast parts.
208
209 Options Database Key:
210 . -ts_mprk_type p3 - select this scheme
211
212 Level: advanced
213
214 .seealso: [](ch_ts), `TSMPRK`, `TSMPRKType`, `TSMPRKSetType()`
215 M*/
216
217 /*@C
218 TSMPRKRegisterAll - Registers all of the Partitioned Runge-Kutta explicit methods in `TSMPRK`
219
220 Not Collective, but should be called by all processes which will need the schemes to be registered
221
222 Level: advanced
223
224 .seealso: [](ch_ts), `TSMPRK`, `TSMPRKRegisterDestroy()`
225 @*/
TSMPRKRegisterAll(void)226 PetscErrorCode TSMPRKRegisterAll(void)
227 {
228 PetscFunctionBegin;
229 if (TSMPRKRegisterAllCalled) PetscFunctionReturn(PETSC_SUCCESS);
230 TSMPRKRegisterAllCalled = PETSC_TRUE;
231
232 #define RC PetscRealConstant
233 {
234 const PetscReal Abase[2][2] =
235 {
236 {0, 0},
237 {RC(1.0), 0}
238 },
239 bbase[2] = {RC(0.5), RC(0.5)};
240 PetscReal Asb[4][4] = {{0}}, Af[4][4] = {{0}}, bsb[4] = {0}, bf[4] = {0};
241 PetscInt rsb[4] = {0, 0, 1, 2};
242 PetscCall(TSMPRKGenerateTableau2(2, 2, &Abase[0][0], bbase, &Asb[0][0], bsb, &Af[0][0], bf));
243 PetscCall(TSMPRKRegister(TSMPRK2A22, 2, 2, 2, 1, &Asb[0][0], bsb, NULL, rsb, NULL, NULL, NULL, NULL, &Af[0][0], bf, NULL));
244 }
245 {
246 const PetscReal Abase[2][2] =
247 {
248 {0, 0},
249 {RC(1.0), 0}
250 },
251 bbase[2] = {RC(0.5), RC(0.5)};
252 PetscReal Asb[8][8] = {{0}}, Amb[8][8] = {{0}}, Af[8][8] = {{0}}, bsb[8] = {0}, bmb[8] = {0}, bf[8] = {0};
253 PetscInt rsb[8] = {0, 0, 1, 2, 1, 2, 1, 2}, rmb[8] = {0, 0, 1, 2, 0, 0, 5, 6};
254 PetscCall(TSMPRKGenerateTableau3(2, 2, &Abase[0][0], bbase, &Asb[0][0], bsb, &Amb[0][0], bmb, &Af[0][0], bf));
255 PetscCall(TSMPRKRegister(TSMPRK2A23, 2, 2, 2, 2, &Asb[0][0], bsb, NULL, rsb, &Amb[0][0], bmb, NULL, rmb, &Af[0][0], bf, NULL));
256 }
257 {
258 const PetscReal Abase[2][2] =
259 {
260 {0, 0},
261 {RC(1.0), 0}
262 },
263 bbase[2] = {RC(0.5), RC(0.5)};
264 PetscReal Asb[6][6] = {{0}}, Af[6][6] = {{0}}, bsb[6] = {0}, bf[6] = {0};
265 PetscInt rsb[6] = {0, 0, 1, 2, 1, 2};
266 PetscCall(TSMPRKGenerateTableau2(3, 2, &Abase[0][0], bbase, &Asb[0][0], bsb, &Af[0][0], bf));
267 PetscCall(TSMPRKRegister(TSMPRK2A32, 2, 2, 3, 1, &Asb[0][0], bsb, NULL, rsb, NULL, NULL, NULL, NULL, &Af[0][0], bf, NULL));
268 }
269 {
270 const PetscReal Abase[2][2] =
271 {
272 {0, 0},
273 {RC(1.0), 0}
274 },
275 bbase[2] = {RC(0.5), RC(0.5)};
276 PetscReal Asb[18][18] = {{0}}, Amb[18][18] = {{0}}, Af[18][18] = {{0}}, bsb[18] = {0}, bmb[18] = {0}, bf[18] = {0};
277 PetscInt rsb[18] = {0, 0, 1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1, 2}, rmb[18] = {0, 0, 1, 2, 1, 2, 0, 0, 7, 8, 7, 8, 0, 0, 13, 14, 13, 14};
278 PetscCall(TSMPRKGenerateTableau3(3, 2, &Abase[0][0], bbase, &Asb[0][0], bsb, &Amb[0][0], bmb, &Af[0][0], bf));
279 PetscCall(TSMPRKRegister(TSMPRK2A33, 2, 2, 3, 3, &Asb[0][0], bsb, NULL, rsb, &Amb[0][0], bmb, NULL, rmb, &Af[0][0], bf, NULL));
280 }
281 /*
282 PetscReal
283 Asb[8][8] = {{Abase[0][0],Abase[0][1],0,0,0,0,0,0},
284 {Abase[1][0],Abase[1][1],0,0,0,0,0,0},
285 {0,0,Abase[0][0],Abase[0][1],0,0,0,0},
286 {0,0,Abase[1][0],Abase[1][1],0,0,0,0},
287 {0,0,0,0,Abase[0][0],Abase[0][1],0,0},
288 {0,0,0,0,Abase[1][0],Abase[1][1],0,0},
289 {0,0,0,0,0,0,Abase[0][0],Abase[0][1]},
290 {0,0,0,0,0,0,Abase[1][0],Abase[1][1]}},
291 Amb[8][8] = {{Abase[0][0]/m,Abase[0][1]/m,0,0,0,0,0,0},
292 {Abase[1][0]/m,Abase[1][1]/m,0,0,0,0,0,0},
293 {0,0,Abase[0][0]/m,Abase[0][1]/m,0,0,0,0},
294 {0,0,Abase[1][0]/m,Abase[1][1]/m,0,0,0,0},
295 {bbase[0]/m,bbase[1]/m,bbase[0]/m,bbase[1]/m,Abase[0][0]/m,Abase[0][1]/m,0,0},
296 {bbase[0]/m,bbase[1]/m,bbase[0]/m,bbase[1]/m,Abase[1][0]/m,Abase[1][1]/m,0,0},
297 {bbase[0]/m,bbase[1]/m,bbase[0]/m,bbase[1]/m,0,0,Abase[0][0]/m,Abase[0][1]/m},
298 {bbase[0]/m,bbase[1]/m,bbase[0]/m,bbase[1]/m,0,0,Abase[1][0]/m,Abase[1][1]/m}},
299 Af[8][8] = {{Abase[0][0]/m/m,Abase[0][1]/m/m,0,0,0,0,0,0},
300 {Abase[1][0]/m/m,Abase[1][1]/m/m,0,0,0,0,0,0},
301 {0,0,Abase[0][0]/m/m,Abase[0][1]/m/m,0,0,0,0},
302 {0,0,Abase[1][0]/m/m,Abase[1][1]/m/m,0,0,0,0},
303 {bbase[0]/m/m,bbase[1]/m/m,bbase[0]/m/m,bbase[1]/m/m,Abase[0][0]/m/m,Abase[0][1]/m/m,0,0},
304 {bbase[0]/m/m,bbase[1]/m/m,bbase[0]/m/m,bbase[1]/m/m,Abase[1][0]/m/m,Abase[1][1]/m/m,0,0},
305 {bbase[0]/m/m,bbase[1]/m/m,bbase[0]/m/m,bbase[1]/m/m,bbase[0]/m/m,bbase[1]/m/m,Abase[0][0]/m,Abase[0][1]/m},
306 {bbase[0]/m/m,bbase[1]/m/m,bbase[0]/m/m,bbase[1]/m/m,bbase[0]/m/m,bbase[1]/m/m,Abase[1][0]/m,Abase[1][1]/m}},
307 bsb[8] = {bbase[0]/m/m,bbase[1]/m/m,bbase[0]/m/m,bbase[1]/m/m,bbase[0]/m/m,bbase[1]/m/m,bbase[0]/m/m,bbase[1]/m/m},
308 bmb[8] = {bbase[0]/m/m,bbase[1]/m/m,bbase[0]/m/m,bbase[1]/m/m,bbase[0]/m/m,bbase[1]/m/m,bbase[0]/m,bbase[1]/m/m},
309 bf[8] = {bbase[0]/m/m,bbase[1]/m/m,bbase[0]/m/m,bbase[1]/m/m,bbase[0]/m/m,bbase[1]/m,bbase[0]/m/m,bbase[1]/m/m},
310 */
311 /*{
312 const PetscReal
313 As[8][8] = {{0,0,0,0,0,0,0,0},
314 {RC(1.0)/RC(2.0),0,0,0,0,0,0,0},
315 {RC(-1.0)/RC(6.0),RC(2.0)/RC(3.0),0,0,0,0,0,0},
316 {RC(1.0)/RC(3.0),RC(-1.0)/RC(3.0),RC(1.0),0,0,0,0,0},
317 {0,0,0,0,0,0,0,0},
318 {0,0,0,0,RC(1.0)/RC(2.0),0,0,0},
319 {0,0,0,0,RC(-1.0)/RC(6.0),RC(2.0)/RC(3.0),0,0},
320 {0,0,0,0,RC(1.0)/RC(3.0),RC(-1.0)/RC(3.0),RC(1.0),0}},
321 A[8][8] = {{0,0,0,0,0,0,0,0},
322 {RC(1.0)/RC(4.0),0,0,0,0,0,0,0},
323 {RC(-1.0)/RC(12.0),RC(1.0)/RC(3.0),0,0,0,0,0,0},
324 {RC(1.0)/RC(6.0),RC(-1.0)/RC(6.0),RC(1.0)/RC(2.0),0,0,0,0,0},
325 {RC(1.0)/RC(12.0),RC(1.0)/RC(6.0),RC(1.0)/RC(6.0),RC(1.0)/RC(12.0),0,0,0,0},
326 {RC(1.0)/RC(12.0),RC(1.0)/RC(6.0),RC(1.0)/RC(6.0),RC(1.0)/RC(12.0),RC(1.0)/RC(4.0),0,0,0},
327 {RC(1.0)/RC(12.0),RC(1.0)/RC(6.0),RC(1.0)/RC(6.0),RC(1.0)/RC(12.0),RC(-1.0)/RC(12.0),RC(1.0)/RC(3.0),0,0},
328 {RC(1.0)/RC(12.0),RC(1.0)/RC(6.0),RC(1.0)/RC(6.0),RC(1.0)/RC(12.0),RC(1.0)/RC(6.0),RC(-1.0)/RC(6.0),RC(1.0)/RC(2.0),0}},
329 bs[8] = {RC(1.0)/RC(12.0),RC(1.0)/RC(6.0),RC(1.0)/RC(6.0),RC(1.0)/RC(12.0),RC(1.0)/RC(12.0),RC(1.0)/RC(6.0),RC(1.0)/RC(6.0),RC(1.0)/RC(12.0)},
330 b[8] = {RC(1.0)/RC(12.0),RC(1.0)/RC(6.0),RC(1.0)/RC(6.0),RC(1.0)/RC(12.0),RC(1.0)/RC(12.0),RC(1.0)/RC(6.0),RC(1.0)/RC(6.0),RC(1.0)/RC(12.0)};
331 PetscCall(TSMPRKRegister(TSMPRKPM3,3,8,&As[0][0],bs,NULL,&A[0][0],b,NULL));
332 }*/
333
334 {
335 const PetscReal Asb[5][5] =
336 {
337 {0, 0, 0, 0, 0},
338 {RC(1.0) / RC(2.0), 0, 0, 0, 0},
339 {RC(1.0) / RC(2.0), 0, 0, 0, 0},
340 {RC(1.0), 0, 0, 0, 0},
341 {RC(1.0), 0, 0, 0, 0}
342 },
343 Af[5][5] = {{0, 0, 0, 0, 0}, {RC(1.0) / RC(2.0), 0, 0, 0, 0}, {RC(1.0) / RC(4.0), RC(1.0) / RC(4.0), 0, 0, 0}, {RC(1.0) / RC(4.0), RC(1.0) / RC(4.0), RC(1.0) / RC(2.0), 0, 0}, {RC(1.0) / RC(4.0), RC(1.0) / RC(4.0), RC(1.0) / RC(4.0), RC(1.0) / RC(4.0), 0}}, bsb[5] = {RC(1.0) / RC(2.0), 0, 0, 0, RC(1.0) / RC(2.0)}, bf[5] = {RC(1.0) / RC(4.0), RC(1.0) / RC(4.0), RC(1.0) / RC(4.0), RC(1.0) / RC(4.0), 0};
344 const PetscInt rsb[5] = {0, 0, 2, 0, 4};
345 PetscCall(TSMPRKRegister(TSMPRKP2, 2, 5, 1, 1, &Asb[0][0], bsb, NULL, rsb, NULL, NULL, NULL, NULL, &Af[0][0], bf, NULL));
346 }
347
348 {
349 const PetscReal Asb[10][10] =
350 {
351 {0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
352 {RC(1.0) / RC(4.0), 0, 0, 0, 0, 0, 0, 0, 0, 0},
353 {RC(1.0) / RC(4.0), 0, 0, 0, 0, 0, 0, 0, 0, 0},
354 {RC(1.0) / RC(2.0), 0, 0, 0, 0, 0, 0, 0, 0, 0},
355 {RC(1.0) / RC(2.0), 0, 0, 0, 0, 0, 0, 0, 0, 0},
356 {RC(-1.0) / RC(6.0), 0, 0, 0, RC(2.0) / RC(3.0), 0, 0, 0, 0, 0},
357 {RC(1.0) / RC(12.0), 0, 0, 0, RC(1.0) / RC(6.0), RC(1.0) / RC(2.0), 0, 0, 0, 0},
358 {RC(1.0) / RC(12.0), 0, 0, 0, RC(1.0) / RC(6.0), RC(1.0) / RC(2.0), 0, 0, 0, 0},
359 {RC(1.0) / RC(3.0), 0, 0, 0, RC(-1.0) / RC(3.0), RC(1.0), 0, 0, 0, 0},
360 {RC(1.0) / RC(3.0), 0, 0, 0, RC(-1.0) / RC(3.0), RC(1.0), 0, 0, 0, 0}
361 },
362 Af[10][10] = {{0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {RC(1.0) / RC(4.0), 0, 0, 0, 0, 0, 0, 0, 0, 0}, {RC(-1.0) / RC(12.0), RC(1.0) / RC(3.0), 0, 0, 0, 0, 0, 0, 0, 0}, {RC(1.0) / RC(6.0), RC(-1.0) / RC(6.0), RC(1.0) / RC(2.0), 0, 0, 0, 0, 0, 0, 0}, {RC(1.0) / RC(12.0), RC(1.0) / RC(6.0), RC(1.0) / RC(6.0), RC(1.0) / RC(12.0), 0, 0, 0, 0, 0, 0}, {RC(1.0) / RC(12.0), RC(1.0) / RC(6.0), RC(1.0) / RC(6.0), RC(1.0) / RC(12.0), 0, 0, 0, 0, 0, 0}, {RC(1.0) / RC(12.0), RC(1.0) / RC(6.0), RC(1.0) / RC(6.0), RC(1.0) / RC(12.0), 0, RC(1.0) / RC(4.0), 0, 0, 0, 0}, {RC(1.0) / RC(12.0), RC(1.0) / RC(6.0), RC(1.0) / RC(6.0), RC(1.0) / RC(12.0), 0, RC(-1.0) / RC(12.0), RC(1.0) / RC(3.0), 0, 0, 0}, {RC(1.0) / RC(12.0), RC(1.0) / RC(6.0), RC(1.0) / RC(6.0), RC(1.0) / RC(12.0), 0, RC(1.0) / RC(6.0), RC(-1.0) / RC(6.0), RC(1.0) / RC(2.0), 0, 0}, {RC(1.0) / RC(12.0), RC(1.0) / RC(6.0), RC(1.0) / RC(6.0), RC(1.0) / RC(12.0), 0, RC(1.0) / RC(12.0), RC(1.0) / RC(6.0), RC(1.0) / RC(6.0), RC(1.0) / RC(12.0), 0}}, bsb[10] = {RC(1.0) / RC(6.0), 0, 0, 0, RC(1.0) / RC(3.0), RC(1.0) / RC(3.0), 0, 0, 0, RC(1.0) / RC(6.0)}, bf[10] = {RC(1.0) / RC(12.0), RC(1.0) / RC(6.0), RC(1.0) / RC(6.0), RC(1.0) / RC(12.0), 0, RC(1.0) / RC(12.0), RC(1.0) / RC(6.0), RC(1.0) / RC(6.0), RC(1.0) / RC(12.0), 0};
363 const PetscInt rsb[10] = {0, 0, 2, 0, 4, 0, 0, 7, 0, 9};
364 PetscCall(TSMPRKRegister(TSMPRKP3, 3, 5, 2, 1, &Asb[0][0], bsb, NULL, rsb, NULL, NULL, NULL, NULL, &Af[0][0], bf, NULL));
365 }
366 #undef RC
367 PetscFunctionReturn(PETSC_SUCCESS);
368 }
369
370 /*@C
371 TSMPRKRegisterDestroy - Frees the list of schemes that were registered by `TSMPRKRegister()`.
372
373 Not Collective
374
375 Level: advanced
376
377 .seealso: [](ch_ts), `TSMPRK`, `TSMPRKRegister()`, `TSMPRKRegisterAll()`
378 @*/
TSMPRKRegisterDestroy(void)379 PetscErrorCode TSMPRKRegisterDestroy(void)
380 {
381 MPRKTableauLink link;
382
383 PetscFunctionBegin;
384 while ((link = MPRKTableauList)) {
385 MPRKTableau t = &link->tab;
386 MPRKTableauList = link->next;
387 PetscCall(PetscFree3(t->Asb, t->bsb, t->csb));
388 PetscCall(PetscFree3(t->Amb, t->bmb, t->cmb));
389 PetscCall(PetscFree3(t->Af, t->bf, t->cf));
390 PetscCall(PetscFree(t->rsb));
391 PetscCall(PetscFree(t->rmb));
392 PetscCall(PetscFree(t->name));
393 PetscCall(PetscFree(link));
394 }
395 TSMPRKRegisterAllCalled = PETSC_FALSE;
396 PetscFunctionReturn(PETSC_SUCCESS);
397 }
398
399 /*@C
400 TSMPRKInitializePackage - This function initializes everything in the `TSMPRK` package. It is called
401 from `PetscDLLibraryRegister()` when using dynamic libraries, and on the first call to `TSCreate_MPRK()`
402 when using static libraries.
403
404 Level: developer
405
406 .seealso: [](ch_ts), `TSMPRK`, `PetscInitialize()`
407 @*/
TSMPRKInitializePackage(void)408 PetscErrorCode TSMPRKInitializePackage(void)
409 {
410 PetscFunctionBegin;
411 if (TSMPRKPackageInitialized) PetscFunctionReturn(PETSC_SUCCESS);
412 TSMPRKPackageInitialized = PETSC_TRUE;
413 PetscCall(TSMPRKRegisterAll());
414 PetscCall(PetscRegisterFinalize(TSMPRKFinalizePackage));
415 PetscFunctionReturn(PETSC_SUCCESS);
416 }
417
418 /*@C
419 TSMPRKFinalizePackage - This function destroys everything in the `TSMPRK` package. It is
420 called from `PetscFinalize()`.
421
422 Level: developer
423
424 .seealso: [](ch_ts), `TSMPRK`, `PetscFinalize()`
425 @*/
TSMPRKFinalizePackage(void)426 PetscErrorCode TSMPRKFinalizePackage(void)
427 {
428 PetscFunctionBegin;
429 TSMPRKPackageInitialized = PETSC_FALSE;
430 PetscCall(TSMPRKRegisterDestroy());
431 PetscFunctionReturn(PETSC_SUCCESS);
432 }
433
434 /*@C
435 TSMPRKRegister - register a `TSMPRK` scheme by providing the entries in the Butcher tableau
436
437 Not Collective, but the same schemes should be registered on all processes on which they will be used, No Fortran Support
438
439 Input Parameters:
440 + name - identifier for method
441 . order - approximation order of method
442 . sbase - number of stages in the base methods
443 . ratio1 - stepsize ratio at 1st level (e.g. slow/medium)
444 . ratio2 - stepsize ratio at 2nd level (e.g. medium/fast)
445 . Asb - stage coefficients for slow components(dimension s*s, row-major)
446 . bsb - step completion table for slow components(dimension s)
447 . csb - abscissa for slow components(dimension s)
448 . rsb - array of flags for repeated stages for slow components (dimension s)
449 . Amb - stage coefficients for medium components(dimension s*s, row-major)
450 . bmb - step completion table for medium components(dimension s)
451 . cmb - abscissa for medium components(dimension s)
452 . rmb - array of flags for repeated stages for medium components (dimension s)
453 . Af - stage coefficients for fast components(dimension s*s, row-major)
454 . bf - step completion table for fast components(dimension s)
455 - cf - abscissa for fast components(dimension s)
456
457 Level: advanced
458
459 Note:
460 Several `TSMPRK` methods are provided, this function is only needed to create new methods.
461
462 .seealso: [](ch_ts), `TSMPRK`
463 @*/
TSMPRKRegister(TSMPRKType name,PetscInt order,PetscInt sbase,PetscInt ratio1,PetscInt ratio2,const PetscReal Asb[],const PetscReal bsb[],const PetscReal csb[],const PetscInt rsb[],const PetscReal Amb[],const PetscReal bmb[],const PetscReal cmb[],const PetscInt rmb[],const PetscReal Af[],const PetscReal bf[],const PetscReal cf[])464 PetscErrorCode TSMPRKRegister(TSMPRKType name, PetscInt order, PetscInt sbase, PetscInt ratio1, PetscInt ratio2, const PetscReal Asb[], const PetscReal bsb[], const PetscReal csb[], const PetscInt rsb[], const PetscReal Amb[], const PetscReal bmb[], const PetscReal cmb[], const PetscInt rmb[], const PetscReal Af[], const PetscReal bf[], const PetscReal cf[])
465 {
466 MPRKTableauLink link;
467 MPRKTableau t;
468 PetscInt s, i, j;
469
470 PetscFunctionBegin;
471 PetscAssertPointer(name, 1);
472 PetscAssertPointer(Asb, 6);
473 if (bsb) PetscAssertPointer(bsb, 7);
474 if (csb) PetscAssertPointer(csb, 8);
475 if (rsb) PetscAssertPointer(rsb, 9);
476 if (Amb) PetscAssertPointer(Amb, 10);
477 if (bmb) PetscAssertPointer(bmb, 11);
478 if (cmb) PetscAssertPointer(cmb, 12);
479 if (rmb) PetscAssertPointer(rmb, 13);
480 PetscAssertPointer(Af, 14);
481 if (bf) PetscAssertPointer(bf, 15);
482 if (cf) PetscAssertPointer(cf, 16);
483
484 PetscCall(PetscNew(&link));
485 t = &link->tab;
486
487 PetscCall(PetscStrallocpy(name, &t->name));
488 s = sbase * ratio1 * ratio2; /* this is the dimension of the matrices below */
489 t->order = order;
490 t->sbase = sbase;
491 t->s = s;
492 t->np = 2;
493
494 PetscCall(PetscMalloc3(s * s, &t->Af, s, &t->bf, s, &t->cf));
495 PetscCall(PetscArraycpy(t->Af, Af, s * s));
496 if (bf) {
497 PetscCall(PetscArraycpy(t->bf, bf, s));
498 } else
499 for (i = 0; i < s; i++) t->bf[i] = Af[(s - 1) * s + i];
500 if (cf) {
501 PetscCall(PetscArraycpy(t->cf, cf, s));
502 } else {
503 for (i = 0; i < s; i++)
504 for (j = 0, t->cf[i] = 0; j < s; j++) t->cf[i] += Af[i * s + j];
505 }
506
507 if (Amb) {
508 t->np = 3;
509 PetscCall(PetscMalloc3(s * s, &t->Amb, s, &t->bmb, s, &t->cmb));
510 PetscCall(PetscArraycpy(t->Amb, Amb, s * s));
511 if (bmb) {
512 PetscCall(PetscArraycpy(t->bmb, bmb, s));
513 } else {
514 for (i = 0; i < s; i++) t->bmb[i] = Amb[(s - 1) * s + i];
515 }
516 if (cmb) {
517 PetscCall(PetscArraycpy(t->cmb, cmb, s));
518 } else {
519 for (i = 0; i < s; i++)
520 for (j = 0, t->cmb[i] = 0; j < s; j++) t->cmb[i] += Amb[i * s + j];
521 }
522 if (rmb) {
523 PetscCall(PetscMalloc1(s, &t->rmb));
524 PetscCall(PetscArraycpy(t->rmb, rmb, s));
525 } else {
526 PetscCall(PetscCalloc1(s, &t->rmb));
527 }
528 }
529
530 PetscCall(PetscMalloc3(s * s, &t->Asb, s, &t->bsb, s, &t->csb));
531 PetscCall(PetscArraycpy(t->Asb, Asb, s * s));
532 if (bsb) {
533 PetscCall(PetscArraycpy(t->bsb, bsb, s));
534 } else
535 for (i = 0; i < s; i++) t->bsb[i] = Asb[(s - 1) * s + i];
536 if (csb) {
537 PetscCall(PetscArraycpy(t->csb, csb, s));
538 } else {
539 for (i = 0; i < s; i++)
540 for (j = 0, t->csb[i] = 0; j < s; j++) t->csb[i] += Asb[i * s + j];
541 }
542 if (rsb) {
543 PetscCall(PetscMalloc1(s, &t->rsb));
544 PetscCall(PetscArraycpy(t->rsb, rsb, s));
545 } else {
546 PetscCall(PetscCalloc1(s, &t->rsb));
547 }
548 link->next = MPRKTableauList;
549 MPRKTableauList = link;
550 PetscFunctionReturn(PETSC_SUCCESS);
551 }
552
TSMPRKSetSplits(TS ts)553 static PetscErrorCode TSMPRKSetSplits(TS ts)
554 {
555 TS_MPRK *mprk = (TS_MPRK *)ts->data;
556 MPRKTableau tab = mprk->tableau;
557 DM dm, subdm, newdm;
558
559 PetscFunctionBegin;
560 PetscCall(TSRHSSplitGetSubTS(ts, "slow", &mprk->subts_slow));
561 PetscCall(TSRHSSplitGetSubTS(ts, "fast", &mprk->subts_fast));
562 PetscCheck(mprk->subts_slow && mprk->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");
563
564 /* Only copy the DM */
565 PetscCall(TSGetDM(ts, &dm));
566
567 PetscCall(TSRHSSplitGetSubTS(ts, "slowbuffer", &mprk->subts_slowbuffer));
568 if (!mprk->subts_slowbuffer) {
569 mprk->subts_slowbuffer = mprk->subts_slow;
570 mprk->subts_slow = NULL;
571 }
572 if (mprk->subts_slow) {
573 PetscCall(DMClone(dm, &newdm));
574 PetscCall(TSGetDM(mprk->subts_slow, &subdm));
575 PetscCall(DMCopyDMTS(subdm, newdm));
576 PetscCall(DMCopyDMSNES(subdm, newdm));
577 PetscCall(TSSetDM(mprk->subts_slow, newdm));
578 PetscCall(DMDestroy(&newdm));
579 }
580 PetscCall(DMClone(dm, &newdm));
581 PetscCall(TSGetDM(mprk->subts_slowbuffer, &subdm));
582 PetscCall(DMCopyDMTS(subdm, newdm));
583 PetscCall(DMCopyDMSNES(subdm, newdm));
584 PetscCall(TSSetDM(mprk->subts_slowbuffer, newdm));
585 PetscCall(DMDestroy(&newdm));
586
587 PetscCall(DMClone(dm, &newdm));
588 PetscCall(TSGetDM(mprk->subts_fast, &subdm));
589 PetscCall(DMCopyDMTS(subdm, newdm));
590 PetscCall(DMCopyDMSNES(subdm, newdm));
591 PetscCall(TSSetDM(mprk->subts_fast, newdm));
592 PetscCall(DMDestroy(&newdm));
593
594 if (tab->np == 3) {
595 PetscCall(TSRHSSplitGetSubTS(ts, "medium", &mprk->subts_medium));
596 PetscCall(TSRHSSplitGetSubTS(ts, "mediumbuffer", &mprk->subts_mediumbuffer));
597 if (mprk->subts_medium && !mprk->subts_mediumbuffer) {
598 mprk->subts_mediumbuffer = mprk->subts_medium;
599 mprk->subts_medium = NULL;
600 }
601 if (mprk->subts_medium) {
602 PetscCall(DMClone(dm, &newdm));
603 PetscCall(TSGetDM(mprk->subts_medium, &subdm));
604 PetscCall(DMCopyDMTS(subdm, newdm));
605 PetscCall(DMCopyDMSNES(subdm, newdm));
606 PetscCall(TSSetDM(mprk->subts_medium, newdm));
607 PetscCall(DMDestroy(&newdm));
608 }
609 PetscCall(DMClone(dm, &newdm));
610 PetscCall(TSGetDM(mprk->subts_mediumbuffer, &subdm));
611 PetscCall(DMCopyDMTS(subdm, newdm));
612 PetscCall(DMCopyDMSNES(subdm, newdm));
613 PetscCall(TSSetDM(mprk->subts_mediumbuffer, newdm));
614 PetscCall(DMDestroy(&newdm));
615 }
616 PetscFunctionReturn(PETSC_SUCCESS);
617 }
618
619 /*
620 This if for nonsplit RHS MPRK
621 The step completion formula is
622
623 x1 = x0 + h b^T YdotRHS
624
625 */
TSEvaluateStep_MPRK(TS ts,PetscInt order,Vec X,PetscBool * done)626 static PetscErrorCode TSEvaluateStep_MPRK(TS ts, PetscInt order, Vec X, PetscBool *done)
627 {
628 TS_MPRK *mprk = (TS_MPRK *)ts->data;
629 MPRKTableau tab = mprk->tableau;
630 PetscScalar *wf = mprk->work_fast;
631 PetscReal h = ts->time_step;
632 PetscInt s = tab->s, j;
633
634 PetscFunctionBegin;
635 for (j = 0; j < s; j++) wf[j] = h * tab->bf[j];
636 PetscCall(VecCopy(ts->vec_sol, X));
637 PetscCall(VecMAXPY(X, s, wf, mprk->YdotRHS));
638 PetscFunctionReturn(PETSC_SUCCESS);
639 }
640
TSStep_MPRK(TS ts)641 static PetscErrorCode TSStep_MPRK(TS ts)
642 {
643 TS_MPRK *mprk = (TS_MPRK *)ts->data;
644 Vec *Y = mprk->Y, *YdotRHS = mprk->YdotRHS, *YdotRHS_fast = mprk->YdotRHS_fast, *YdotRHS_slow = mprk->YdotRHS_slow, *YdotRHS_slowbuffer = mprk->YdotRHS_slowbuffer;
645 Vec Yslow, Yslowbuffer, Yfast;
646 MPRKTableau tab = mprk->tableau;
647 const PetscInt s = tab->s;
648 const PetscReal *Af = tab->Af, *cf = tab->cf, *Asb = tab->Asb, *csb = tab->csb;
649 PetscScalar *wf = mprk->work_fast, *wsb = mprk->work_slowbuffer;
650 PetscInt i, j;
651 PetscReal next_time_step = ts->time_step, t = ts->ptime, h = ts->time_step;
652
653 PetscFunctionBegin;
654 for (i = 0; i < s; i++) {
655 mprk->stage_time = t + h * cf[i];
656 PetscCall(TSPreStage(ts, mprk->stage_time));
657 PetscCall(VecCopy(ts->vec_sol, Y[i]));
658
659 /* slow buffer region */
660 for (j = 0; j < i; j++) wsb[j] = h * Asb[i * s + j];
661 for (j = 0; j < i; j++) PetscCall(VecGetSubVector(YdotRHS[j], mprk->is_slowbuffer, &YdotRHS_slowbuffer[j]));
662 PetscCall(VecGetSubVector(Y[i], mprk->is_slowbuffer, &Yslowbuffer));
663 PetscCall(VecMAXPY(Yslowbuffer, i, wsb, mprk->YdotRHS_slowbuffer));
664 PetscCall(VecRestoreSubVector(Y[i], mprk->is_slowbuffer, &Yslowbuffer));
665 for (j = 0; j < i; j++) PetscCall(VecRestoreSubVector(YdotRHS[j], mprk->is_slowbuffer, &YdotRHS_slowbuffer[j]));
666 /* slow region */
667 if (mprk->is_slow) {
668 for (j = 0; j < i; j++) PetscCall(VecGetSubVector(YdotRHS[j], mprk->is_slow, &YdotRHS_slow[j]));
669 PetscCall(VecGetSubVector(Y[i], mprk->is_slow, &Yslow));
670 PetscCall(VecMAXPY(Yslow, i, wsb, mprk->YdotRHS_slow));
671 PetscCall(VecRestoreSubVector(Y[i], mprk->is_slow, &Yslow));
672 for (j = 0; j < i; j++) PetscCall(VecRestoreSubVector(YdotRHS[j], mprk->is_slow, &YdotRHS_slow[j]));
673 }
674
675 /* fast region */
676 for (j = 0; j < i; j++) wf[j] = h * Af[i * s + j];
677 for (j = 0; j < i; j++) PetscCall(VecGetSubVector(YdotRHS[j], mprk->is_fast, &YdotRHS_fast[j]));
678 PetscCall(VecGetSubVector(Y[i], mprk->is_fast, &Yfast));
679 PetscCall(VecMAXPY(Yfast, i, wf, mprk->YdotRHS_fast));
680 PetscCall(VecRestoreSubVector(Y[i], mprk->is_fast, &Yfast));
681 for (j = 0; j < i; j++) PetscCall(VecRestoreSubVector(YdotRHS[j], mprk->is_fast, &YdotRHS_fast[j]));
682 if (tab->np == 3) {
683 Vec *YdotRHS_medium = mprk->YdotRHS_medium, *YdotRHS_mediumbuffer = mprk->YdotRHS_mediumbuffer;
684 Vec Ymedium, Ymediumbuffer;
685 PetscScalar *wmb = mprk->work_mediumbuffer;
686
687 for (j = 0; j < i; j++) wmb[j] = h * tab->Amb[i * s + j];
688 /* medium region */
689 if (mprk->is_medium) {
690 for (j = 0; j < i; j++) PetscCall(VecGetSubVector(YdotRHS[j], mprk->is_medium, &YdotRHS_medium[j]));
691 PetscCall(VecGetSubVector(Y[i], mprk->is_medium, &Ymedium));
692 PetscCall(VecMAXPY(Ymedium, i, wmb, mprk->YdotRHS_medium));
693 PetscCall(VecRestoreSubVector(Y[i], mprk->is_medium, &Ymedium));
694 for (j = 0; j < i; j++) PetscCall(VecRestoreSubVector(YdotRHS[j], mprk->is_medium, &YdotRHS_medium[j]));
695 }
696 /* medium buffer region */
697 for (j = 0; j < i; j++) PetscCall(VecGetSubVector(YdotRHS[j], mprk->is_mediumbuffer, &YdotRHS_mediumbuffer[j]));
698 PetscCall(VecGetSubVector(Y[i], mprk->is_mediumbuffer, &Ymediumbuffer));
699 PetscCall(VecMAXPY(Ymediumbuffer, i, wmb, mprk->YdotRHS_mediumbuffer));
700 PetscCall(VecRestoreSubVector(Y[i], mprk->is_mediumbuffer, &Ymediumbuffer));
701 for (j = 0; j < i; j++) PetscCall(VecRestoreSubVector(YdotRHS[j], mprk->is_mediumbuffer, &YdotRHS_mediumbuffer[j]));
702 }
703 PetscCall(TSPostStage(ts, mprk->stage_time, i, Y));
704 /* compute the stage RHS by fast and slow tableau respectively */
705 PetscCall(TSComputeRHSFunction(ts, t + h * csb[i], Y[i], YdotRHS[i]));
706 }
707 PetscCall(TSEvaluateStep(ts, tab->order, ts->vec_sol, NULL));
708 ts->ptime += ts->time_step;
709 ts->time_step = next_time_step;
710 PetscFunctionReturn(PETSC_SUCCESS);
711 }
712
713 /*
714 This if for the case when split RHS is used
715 The step completion formula is
716 x1 = x0 + h b^T YdotRHS
717 */
TSEvaluateStep_MPRKSPLIT(TS ts,PetscInt order,Vec X,PetscBool * done)718 static PetscErrorCode TSEvaluateStep_MPRKSPLIT(TS ts, PetscInt order, Vec X, PetscBool *done)
719 {
720 TS_MPRK *mprk = (TS_MPRK *)ts->data;
721 MPRKTableau tab = mprk->tableau;
722 Vec Xslow, Xfast, Xslowbuffer; /* subvectors for slow and fast components in X respectively */
723 PetscScalar *wf = mprk->work_fast, *ws = mprk->work_slow, *wsb = mprk->work_slowbuffer;
724 PetscReal h = ts->time_step;
725 PetscInt s = tab->s, j, computedstages;
726
727 PetscFunctionBegin;
728 PetscCall(VecCopy(ts->vec_sol, X));
729
730 /* slow region */
731 if (mprk->is_slow) {
732 computedstages = 0;
733 for (j = 0; j < s; j++) {
734 if (tab->rsb[j]) ws[tab->rsb[j] - 1] += h * tab->bsb[j];
735 else ws[computedstages++] = h * tab->bsb[j];
736 }
737 PetscCall(VecGetSubVector(X, mprk->is_slow, &Xslow));
738 PetscCall(VecMAXPY(Xslow, computedstages, ws, mprk->YdotRHS_slow));
739 PetscCall(VecRestoreSubVector(X, mprk->is_slow, &Xslow));
740 }
741
742 if (tab->np == 3 && mprk->is_medium) {
743 computedstages = 0;
744 for (j = 0; j < s; j++) {
745 if (tab->rmb[j]) wsb[computedstages - tab->sbase + (tab->rmb[j] - 1) % tab->sbase] += h * tab->bsb[j];
746 else wsb[computedstages++] = h * tab->bsb[j];
747 }
748 PetscCall(VecGetSubVector(X, mprk->is_slowbuffer, &Xslowbuffer));
749 PetscCall(VecMAXPY(Xslowbuffer, computedstages, wsb, mprk->YdotRHS_slowbuffer));
750 PetscCall(VecRestoreSubVector(X, mprk->is_slowbuffer, &Xslowbuffer));
751 } else {
752 /* slow buffer region */
753 for (j = 0; j < s; j++) wsb[j] = h * tab->bsb[j];
754 PetscCall(VecGetSubVector(X, mprk->is_slowbuffer, &Xslowbuffer));
755 PetscCall(VecMAXPY(Xslowbuffer, s, wsb, mprk->YdotRHS_slowbuffer));
756 PetscCall(VecRestoreSubVector(X, mprk->is_slowbuffer, &Xslowbuffer));
757 }
758 if (tab->np == 3) {
759 Vec Xmedium, Xmediumbuffer;
760 PetscScalar *wm = mprk->work_medium, *wmb = mprk->work_mediumbuffer;
761 /* medium region and slow buffer region */
762 if (mprk->is_medium) {
763 computedstages = 0;
764 for (j = 0; j < s; j++) {
765 if (tab->rmb[j]) wm[computedstages - tab->sbase + (tab->rmb[j] - 1) % tab->sbase] += h * tab->bmb[j];
766 else wm[computedstages++] = h * tab->bmb[j];
767 }
768 PetscCall(VecGetSubVector(X, mprk->is_medium, &Xmedium));
769 PetscCall(VecMAXPY(Xmedium, computedstages, wm, mprk->YdotRHS_medium));
770 PetscCall(VecRestoreSubVector(X, mprk->is_medium, &Xmedium));
771 }
772 /* medium buffer region */
773 for (j = 0; j < s; j++) wmb[j] = h * tab->bmb[j];
774 PetscCall(VecGetSubVector(X, mprk->is_mediumbuffer, &Xmediumbuffer));
775 PetscCall(VecMAXPY(Xmediumbuffer, s, wmb, mprk->YdotRHS_mediumbuffer));
776 PetscCall(VecRestoreSubVector(X, mprk->is_mediumbuffer, &Xmediumbuffer));
777 }
778 /* fast region */
779 for (j = 0; j < s; j++) wf[j] = h * tab->bf[j];
780 PetscCall(VecGetSubVector(X, mprk->is_fast, &Xfast));
781 PetscCall(VecMAXPY(Xfast, s, wf, mprk->YdotRHS_fast));
782 PetscCall(VecRestoreSubVector(X, mprk->is_fast, &Xfast));
783 PetscFunctionReturn(PETSC_SUCCESS);
784 }
785
TSStep_MPRKSPLIT(TS ts)786 static PetscErrorCode TSStep_MPRKSPLIT(TS ts)
787 {
788 TS_MPRK *mprk = (TS_MPRK *)ts->data;
789 MPRKTableau tab = mprk->tableau;
790 Vec *Y = mprk->Y, *YdotRHS_fast = mprk->YdotRHS_fast, *YdotRHS_slow = mprk->YdotRHS_slow, *YdotRHS_slowbuffer = mprk->YdotRHS_slowbuffer;
791 Vec Yslow, Yslowbuffer, Yfast; /* subvectors for slow and fast components in Y[i] respectively */
792 PetscInt s = tab->s;
793 const PetscReal *Af = tab->Af, *cf = tab->cf, *Asb = tab->Asb, *csb = tab->csb;
794 PetscScalar *wf = mprk->work_fast, *ws = mprk->work_slow, *wsb = mprk->work_slowbuffer;
795 PetscInt i, j, computedstages;
796 PetscReal next_time_step = ts->time_step, t = ts->ptime, h = ts->time_step;
797
798 PetscFunctionBegin;
799 for (i = 0; i < s; i++) {
800 mprk->stage_time = t + h * cf[i];
801 PetscCall(TSPreStage(ts, mprk->stage_time));
802 /* calculate the stage value for fast and slow components respectively */
803 PetscCall(VecCopy(ts->vec_sol, Y[i]));
804 for (j = 0; j < i; j++) wsb[j] = h * Asb[i * s + j];
805
806 /* slow buffer region */
807 if (tab->np == 3 && mprk->is_medium) {
808 if (tab->rmb[i]) {
809 PetscCall(VecGetSubVector(Y[i], mprk->is_slowbuffer, &Yslowbuffer));
810 PetscCall(VecISCopy(Y[tab->rmb[i] - 1], mprk->is_slowbuffer, SCATTER_REVERSE, Yslowbuffer));
811 PetscCall(VecRestoreSubVector(Y[i], mprk->is_slowbuffer, &Yslowbuffer));
812 } else {
813 PetscScalar *wm = mprk->work_medium;
814 computedstages = 0;
815 for (j = 0; j < i; j++) {
816 if (tab->rmb[j]) wm[computedstages - tab->sbase + (tab->rmb[j] - 1) % tab->sbase] += wsb[j];
817 else wm[computedstages++] = wsb[j];
818 }
819 PetscCall(VecGetSubVector(Y[i], mprk->is_slowbuffer, &Yslowbuffer));
820 PetscCall(VecMAXPY(Yslowbuffer, computedstages, wm, YdotRHS_slowbuffer));
821 PetscCall(VecRestoreSubVector(Y[i], mprk->is_slowbuffer, &Yslowbuffer));
822 }
823 } else {
824 PetscCall(VecGetSubVector(Y[i], mprk->is_slowbuffer, &Yslowbuffer));
825 PetscCall(VecMAXPY(Yslowbuffer, i, wsb, YdotRHS_slowbuffer));
826 PetscCall(VecRestoreSubVector(Y[i], mprk->is_slowbuffer, &Yslowbuffer));
827 }
828
829 /* slow region */
830 if (mprk->is_slow) {
831 if (tab->rsb[i]) { /* repeat previous stage */
832 PetscCall(VecGetSubVector(Y[i], mprk->is_slow, &Yslow));
833 PetscCall(VecISCopy(Y[tab->rsb[i] - 1], mprk->is_slow, SCATTER_REVERSE, Yslow));
834 PetscCall(VecRestoreSubVector(Y[i], mprk->is_slow, &Yslow));
835 } else {
836 computedstages = 0;
837 for (j = 0; j < i; j++) {
838 if (tab->rsb[j]) ws[tab->rsb[j] - 1] += wsb[j];
839 else ws[computedstages++] = wsb[j];
840 }
841 PetscCall(VecGetSubVector(Y[i], mprk->is_slow, &Yslow));
842 PetscCall(VecMAXPY(Yslow, computedstages, ws, YdotRHS_slow));
843 PetscCall(VecRestoreSubVector(Y[i], mprk->is_slow, &Yslow));
844 /* only depends on the slow buffer region */
845 PetscCall(TSComputeRHSFunction(mprk->subts_slow, t + h * csb[i], Y[i], YdotRHS_slow[computedstages]));
846 }
847 }
848
849 /* fast region */
850 for (j = 0; j < i; j++) wf[j] = h * Af[i * s + j];
851 PetscCall(VecGetSubVector(Y[i], mprk->is_fast, &Yfast));
852 PetscCall(VecMAXPY(Yfast, i, wf, YdotRHS_fast));
853 PetscCall(VecRestoreSubVector(Y[i], mprk->is_fast, &Yfast));
854
855 if (tab->np == 3) {
856 Vec *YdotRHS_medium = mprk->YdotRHS_medium, *YdotRHS_mediumbuffer = mprk->YdotRHS_mediumbuffer;
857 Vec Ymedium, Ymediumbuffer;
858 const PetscReal *Amb = tab->Amb, *cmb = tab->cmb;
859 PetscScalar *wm = mprk->work_medium, *wmb = mprk->work_mediumbuffer;
860
861 for (j = 0; j < i; j++) wmb[j] = h * Amb[i * s + j];
862 /* medium buffer region */
863 PetscCall(VecGetSubVector(Y[i], mprk->is_mediumbuffer, &Ymediumbuffer));
864 PetscCall(VecMAXPY(Ymediumbuffer, i, wmb, YdotRHS_mediumbuffer));
865 PetscCall(VecRestoreSubVector(Y[i], mprk->is_mediumbuffer, &Ymediumbuffer));
866 /* medium region */
867 if (mprk->is_medium) {
868 if (tab->rmb[i]) { /* repeat previous stage */
869 PetscCall(VecGetSubVector(Y[i], mprk->is_medium, &Ymedium));
870 PetscCall(VecISCopy(Y[tab->rmb[i] - 1], mprk->is_medium, SCATTER_REVERSE, Ymedium));
871 PetscCall(VecRestoreSubVector(Y[i], mprk->is_medium, &Ymedium));
872 } else {
873 computedstages = 0;
874 for (j = 0; j < i; j++) {
875 if (tab->rmb[j]) wm[computedstages - tab->sbase + (tab->rmb[j] - 1) % tab->sbase] += wmb[j];
876 else wm[computedstages++] = wmb[j];
877 }
878 PetscCall(VecGetSubVector(Y[i], mprk->is_medium, &Ymedium));
879 PetscCall(VecMAXPY(Ymedium, computedstages, wm, YdotRHS_medium));
880 PetscCall(VecRestoreSubVector(Y[i], mprk->is_medium, &Ymedium));
881 /* only depends on the medium buffer region */
882 PetscCall(TSComputeRHSFunction(mprk->subts_medium, t + h * cmb[i], Y[i], YdotRHS_medium[computedstages]));
883 /* must be computed after all regions are updated in Y */
884 PetscCall(TSComputeRHSFunction(mprk->subts_slowbuffer, t + h * csb[i], Y[i], YdotRHS_slowbuffer[computedstages]));
885 }
886 }
887 /* must be computed after fast region and slow region are updated in Y */
888 PetscCall(TSComputeRHSFunction(mprk->subts_mediumbuffer, t + h * cmb[i], Y[i], YdotRHS_mediumbuffer[i]));
889 }
890 if (!(tab->np == 3 && mprk->is_medium)) PetscCall(TSComputeRHSFunction(mprk->subts_slowbuffer, t + h * csb[i], Y[i], YdotRHS_slowbuffer[i]));
891 PetscCall(TSComputeRHSFunction(mprk->subts_fast, t + h * cf[i], Y[i], YdotRHS_fast[i]));
892 }
893
894 PetscCall(TSEvaluateStep(ts, tab->order, ts->vec_sol, NULL));
895 ts->ptime += ts->time_step;
896 ts->time_step = next_time_step;
897 PetscFunctionReturn(PETSC_SUCCESS);
898 }
899
TSMPRKTableauReset(TS ts)900 static PetscErrorCode TSMPRKTableauReset(TS ts)
901 {
902 TS_MPRK *mprk = (TS_MPRK *)ts->data;
903 MPRKTableau tab = mprk->tableau;
904
905 PetscFunctionBegin;
906 if (!tab) PetscFunctionReturn(PETSC_SUCCESS);
907 PetscCall(PetscFree(mprk->work_fast));
908 PetscCall(PetscFree(mprk->work_slow));
909 PetscCall(PetscFree(mprk->work_slowbuffer));
910 PetscCall(PetscFree(mprk->work_medium));
911 PetscCall(PetscFree(mprk->work_mediumbuffer));
912 PetscCall(VecDestroyVecs(tab->s, &mprk->Y));
913 if (ts->use_splitrhsfunction) {
914 PetscCall(VecDestroyVecs(tab->s, &mprk->YdotRHS_fast));
915 PetscCall(VecDestroyVecs(tab->s, &mprk->YdotRHS_slow));
916 PetscCall(VecDestroyVecs(tab->s, &mprk->YdotRHS_slowbuffer));
917 PetscCall(VecDestroyVecs(tab->s, &mprk->YdotRHS_medium));
918 PetscCall(VecDestroyVecs(tab->s, &mprk->YdotRHS_mediumbuffer));
919 } else {
920 PetscCall(VecDestroyVecs(tab->s, &mprk->YdotRHS));
921 if (mprk->is_slow) PetscCall(PetscFree(mprk->YdotRHS_slow));
922 PetscCall(PetscFree(mprk->YdotRHS_slowbuffer));
923 if (tab->np == 3) {
924 if (mprk->is_medium) PetscCall(PetscFree(mprk->YdotRHS_medium));
925 PetscCall(PetscFree(mprk->YdotRHS_mediumbuffer));
926 }
927 PetscCall(PetscFree(mprk->YdotRHS_fast));
928 }
929 PetscFunctionReturn(PETSC_SUCCESS);
930 }
931
TSReset_MPRK(TS ts)932 static PetscErrorCode TSReset_MPRK(TS ts)
933 {
934 PetscFunctionBegin;
935 PetscCall(TSMPRKTableauReset(ts));
936 PetscFunctionReturn(PETSC_SUCCESS);
937 }
938
DMCoarsenHook_TSMPRK(DM fine,DM coarse,PetscCtx ctx)939 static PetscErrorCode DMCoarsenHook_TSMPRK(DM fine, DM coarse, PetscCtx ctx)
940 {
941 PetscFunctionBegin;
942 PetscFunctionReturn(PETSC_SUCCESS);
943 }
944
DMRestrictHook_TSMPRK(DM fine,Mat restrct,Vec rscale,Mat inject,DM coarse,PetscCtx ctx)945 static PetscErrorCode DMRestrictHook_TSMPRK(DM fine, Mat restrct, Vec rscale, Mat inject, DM coarse, PetscCtx ctx)
946 {
947 PetscFunctionBegin;
948 PetscFunctionReturn(PETSC_SUCCESS);
949 }
950
DMSubDomainHook_TSMPRK(DM dm,DM subdm,PetscCtx ctx)951 static PetscErrorCode DMSubDomainHook_TSMPRK(DM dm, DM subdm, PetscCtx ctx)
952 {
953 PetscFunctionBegin;
954 PetscFunctionReturn(PETSC_SUCCESS);
955 }
956
DMSubDomainRestrictHook_TSMPRK(DM dm,VecScatter gscat,VecScatter lscat,DM subdm,PetscCtx ctx)957 static PetscErrorCode DMSubDomainRestrictHook_TSMPRK(DM dm, VecScatter gscat, VecScatter lscat, DM subdm, PetscCtx ctx)
958 {
959 PetscFunctionBegin;
960 PetscFunctionReturn(PETSC_SUCCESS);
961 }
962
TSMPRKTableauSetUp(TS ts)963 static PetscErrorCode TSMPRKTableauSetUp(TS ts)
964 {
965 TS_MPRK *mprk = (TS_MPRK *)ts->data;
966 MPRKTableau tab = mprk->tableau;
967 Vec YdotRHS_slow, YdotRHS_slowbuffer, YdotRHS_medium, YdotRHS_mediumbuffer, YdotRHS_fast;
968
969 PetscFunctionBegin;
970 PetscCall(VecDuplicateVecs(ts->vec_sol, tab->s, &mprk->Y));
971 if (mprk->is_slow) PetscCall(PetscMalloc1(tab->s, &mprk->work_slow));
972 PetscCall(PetscMalloc1(tab->s, &mprk->work_slowbuffer));
973 if (tab->np == 3) {
974 if (mprk->is_medium) PetscCall(PetscMalloc1(tab->s, &mprk->work_medium));
975 PetscCall(PetscMalloc1(tab->s, &mprk->work_mediumbuffer));
976 }
977 PetscCall(PetscMalloc1(tab->s, &mprk->work_fast));
978
979 if (ts->use_splitrhsfunction) {
980 if (mprk->is_slow) {
981 PetscCall(VecGetSubVector(ts->vec_sol, mprk->is_slow, &YdotRHS_slow));
982 PetscCall(VecDuplicateVecs(YdotRHS_slow, tab->s, &mprk->YdotRHS_slow));
983 PetscCall(VecRestoreSubVector(ts->vec_sol, mprk->is_slow, &YdotRHS_slow));
984 }
985 PetscCall(VecGetSubVector(ts->vec_sol, mprk->is_slowbuffer, &YdotRHS_slowbuffer));
986 PetscCall(VecDuplicateVecs(YdotRHS_slowbuffer, tab->s, &mprk->YdotRHS_slowbuffer));
987 PetscCall(VecRestoreSubVector(ts->vec_sol, mprk->is_slowbuffer, &YdotRHS_slowbuffer));
988 if (tab->np == 3) {
989 if (mprk->is_medium) {
990 PetscCall(VecGetSubVector(ts->vec_sol, mprk->is_medium, &YdotRHS_medium));
991 PetscCall(VecDuplicateVecs(YdotRHS_medium, tab->s, &mprk->YdotRHS_medium));
992 PetscCall(VecRestoreSubVector(ts->vec_sol, mprk->is_medium, &YdotRHS_medium));
993 }
994 PetscCall(VecGetSubVector(ts->vec_sol, mprk->is_mediumbuffer, &YdotRHS_mediumbuffer));
995 PetscCall(VecDuplicateVecs(YdotRHS_mediumbuffer, tab->s, &mprk->YdotRHS_mediumbuffer));
996 PetscCall(VecRestoreSubVector(ts->vec_sol, mprk->is_mediumbuffer, &YdotRHS_mediumbuffer));
997 }
998 PetscCall(VecGetSubVector(ts->vec_sol, mprk->is_fast, &YdotRHS_fast));
999 PetscCall(VecDuplicateVecs(YdotRHS_fast, tab->s, &mprk->YdotRHS_fast));
1000 PetscCall(VecRestoreSubVector(ts->vec_sol, mprk->is_fast, &YdotRHS_fast));
1001 } else {
1002 PetscCall(VecDuplicateVecs(ts->vec_sol, tab->s, &mprk->YdotRHS));
1003 if (mprk->is_slow) PetscCall(PetscMalloc1(tab->s, &mprk->YdotRHS_slow));
1004 PetscCall(PetscMalloc1(tab->s, &mprk->YdotRHS_slowbuffer));
1005 if (tab->np == 3) {
1006 if (mprk->is_medium) PetscCall(PetscMalloc1(tab->s, &mprk->YdotRHS_medium));
1007 PetscCall(PetscMalloc1(tab->s, &mprk->YdotRHS_mediumbuffer));
1008 }
1009 PetscCall(PetscMalloc1(tab->s, &mprk->YdotRHS_fast));
1010 }
1011 PetscFunctionReturn(PETSC_SUCCESS);
1012 }
1013
TSSetUp_MPRK(TS ts)1014 static PetscErrorCode TSSetUp_MPRK(TS ts)
1015 {
1016 TS_MPRK *mprk = (TS_MPRK *)ts->data;
1017 MPRKTableau tab = mprk->tableau;
1018 DM dm;
1019
1020 PetscFunctionBegin;
1021 PetscCall(TSRHSSplitGetIS(ts, "slow", &mprk->is_slow));
1022 PetscCall(TSRHSSplitGetIS(ts, "fast", &mprk->is_fast));
1023 PetscCheck(mprk->is_slow && mprk->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 the method '%s'", tab->name);
1024
1025 if (tab->np == 3) {
1026 PetscCall(TSRHSSplitGetIS(ts, "medium", &mprk->is_medium));
1027 PetscCheck(mprk->is_medium, PetscObjectComm((PetscObject)ts), PETSC_ERR_USER, "Must set up RHSSplits with TSRHSSplitSetIS() using split names 'slow' and 'medium' and 'fast' respectively in order to use the method '%s'", tab->name);
1028 PetscCall(TSRHSSplitGetIS(ts, "mediumbuffer", &mprk->is_mediumbuffer));
1029 if (!mprk->is_mediumbuffer) { /* let medium buffer cover whole medium region */
1030 mprk->is_mediumbuffer = mprk->is_medium;
1031 mprk->is_medium = NULL;
1032 }
1033 }
1034
1035 /* If users do not provide buffer region settings, the solver will do them automatically, but with a performance penalty */
1036 PetscCall(TSRHSSplitGetIS(ts, "slowbuffer", &mprk->is_slowbuffer));
1037 if (!mprk->is_slowbuffer) { /* let slow buffer cover whole slow region */
1038 mprk->is_slowbuffer = mprk->is_slow;
1039 mprk->is_slow = NULL;
1040 }
1041 PetscCall(TSCheckImplicitTerm(ts));
1042 PetscCall(TSMPRKTableauSetUp(ts));
1043 PetscCall(TSGetDM(ts, &dm));
1044 PetscCall(DMCoarsenHookAdd(dm, DMCoarsenHook_TSMPRK, DMRestrictHook_TSMPRK, ts));
1045 PetscCall(DMSubDomainHookAdd(dm, DMSubDomainHook_TSMPRK, DMSubDomainRestrictHook_TSMPRK, ts));
1046 if (ts->use_splitrhsfunction) {
1047 ts->ops->step = TSStep_MPRKSPLIT;
1048 ts->ops->evaluatestep = TSEvaluateStep_MPRKSPLIT;
1049 PetscCall(TSMPRKSetSplits(ts));
1050 } else {
1051 ts->ops->step = TSStep_MPRK;
1052 ts->ops->evaluatestep = TSEvaluateStep_MPRK;
1053 }
1054 PetscFunctionReturn(PETSC_SUCCESS);
1055 }
1056
TSSetFromOptions_MPRK(TS ts,PetscOptionItems PetscOptionsObject)1057 static PetscErrorCode TSSetFromOptions_MPRK(TS ts, PetscOptionItems PetscOptionsObject)
1058 {
1059 TS_MPRK *mprk = (TS_MPRK *)ts->data;
1060
1061 PetscFunctionBegin;
1062 PetscOptionsHeadBegin(PetscOptionsObject, "PRK ODE solver options");
1063 {
1064 MPRKTableauLink link;
1065 PetscInt count, choice;
1066 PetscBool flg;
1067 const char **namelist;
1068 for (link = MPRKTableauList, count = 0; link; link = link->next, count++);
1069 PetscCall(PetscMalloc1(count, (char ***)&namelist));
1070 for (link = MPRKTableauList, count = 0; link; link = link->next, count++) namelist[count] = link->tab.name;
1071 PetscCall(PetscOptionsEList("-ts_mprk_type", "Family of MPRK method", "TSMPRKSetType", (const char *const *)namelist, count, mprk->tableau->name, &choice, &flg));
1072 if (flg) PetscCall(TSMPRKSetType(ts, namelist[choice]));
1073 PetscCall(PetscFree(namelist));
1074 }
1075 PetscOptionsHeadEnd();
1076 PetscFunctionReturn(PETSC_SUCCESS);
1077 }
1078
TSView_MPRK(TS ts,PetscViewer viewer)1079 static PetscErrorCode TSView_MPRK(TS ts, PetscViewer viewer)
1080 {
1081 TS_MPRK *mprk = (TS_MPRK *)ts->data;
1082 PetscBool isascii;
1083
1084 PetscFunctionBegin;
1085 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
1086 if (isascii) {
1087 MPRKTableau tab = mprk->tableau;
1088 TSMPRKType mprktype;
1089 char fbuf[512];
1090 char sbuf[512];
1091 PetscInt i;
1092 PetscCall(TSMPRKGetType(ts, &mprktype));
1093 PetscCall(PetscViewerASCIIPrintf(viewer, " MPRK type %s\n", mprktype));
1094 PetscCall(PetscViewerASCIIPrintf(viewer, " Order: %" PetscInt_FMT "\n", tab->order));
1095
1096 PetscCall(PetscFormatRealArray(fbuf, sizeof(fbuf), "% 8.6f", tab->s, tab->cf));
1097 PetscCall(PetscViewerASCIIPrintf(viewer, " Abscissa cf = %s\n", fbuf));
1098 PetscCall(PetscViewerASCIIPrintf(viewer, " Af = \n"));
1099 for (i = 0; i < tab->s; i++) {
1100 PetscCall(PetscFormatRealArray(fbuf, sizeof(fbuf), "% 8.6f", tab->s, &tab->Af[i * tab->s]));
1101 PetscCall(PetscViewerASCIIPrintf(viewer, " %s\n", fbuf));
1102 }
1103 PetscCall(PetscFormatRealArray(fbuf, sizeof(fbuf), "% 8.6f", tab->s, tab->bf));
1104 PetscCall(PetscViewerASCIIPrintf(viewer, " bf = %s\n", fbuf));
1105
1106 PetscCall(PetscFormatRealArray(sbuf, sizeof(sbuf), "% 8.6f", tab->s, tab->csb));
1107 PetscCall(PetscViewerASCIIPrintf(viewer, " Abscissa csb = %s\n", sbuf));
1108 PetscCall(PetscViewerASCIIPrintf(viewer, " Asb = \n"));
1109 for (i = 0; i < tab->s; i++) {
1110 PetscCall(PetscFormatRealArray(sbuf, sizeof(sbuf), "% 8.6f", tab->s, &tab->Asb[i * tab->s]));
1111 PetscCall(PetscViewerASCIIPrintf(viewer, " %s\n", sbuf));
1112 }
1113 PetscCall(PetscFormatRealArray(sbuf, sizeof(sbuf), "% 8.6f", tab->s, tab->bsb));
1114 PetscCall(PetscViewerASCIIPrintf(viewer, " bsb = %s\n", sbuf));
1115
1116 if (tab->np == 3) {
1117 char mbuf[512];
1118 PetscCall(PetscFormatRealArray(mbuf, sizeof(mbuf), "% 8.6f", tab->s, tab->cmb));
1119 PetscCall(PetscViewerASCIIPrintf(viewer, " Abscissa cmb = %s\n", mbuf));
1120 PetscCall(PetscViewerASCIIPrintf(viewer, " Amb = \n"));
1121 for (i = 0; i < tab->s; i++) {
1122 PetscCall(PetscFormatRealArray(mbuf, sizeof(mbuf), "% 8.6f", tab->s, &tab->Amb[i * tab->s]));
1123 PetscCall(PetscViewerASCIIPrintf(viewer, " %s\n", mbuf));
1124 }
1125 PetscCall(PetscFormatRealArray(mbuf, sizeof(mbuf), "% 8.6f", tab->s, tab->bmb));
1126 PetscCall(PetscViewerASCIIPrintf(viewer, " bmb = %s\n", mbuf));
1127 }
1128 }
1129 PetscFunctionReturn(PETSC_SUCCESS);
1130 }
1131
TSLoad_MPRK(TS ts,PetscViewer viewer)1132 static PetscErrorCode TSLoad_MPRK(TS ts, PetscViewer viewer)
1133 {
1134 TSAdapt adapt;
1135
1136 PetscFunctionBegin;
1137 PetscCall(TSGetAdapt(ts, &adapt));
1138 PetscCall(TSAdaptLoad(adapt, viewer));
1139 PetscFunctionReturn(PETSC_SUCCESS);
1140 }
1141
1142 /*@
1143 TSMPRKSetType - Set the type of `TSMPRK` scheme
1144
1145 Not Collective
1146
1147 Input Parameters:
1148 + ts - timestepping context
1149 - mprktype - type of `TSMPRK` scheme
1150
1151 Options Database Key:
1152 . -ts_mprk_type - <pm2,p2,p3> - select the specific scheme
1153
1154 Level: intermediate
1155
1156 .seealso: [](ch_ts), `TSMPRKGetType()`, `TSMPRK`, `TSMPRKType`
1157 @*/
TSMPRKSetType(TS ts,TSMPRKType mprktype)1158 PetscErrorCode TSMPRKSetType(TS ts, TSMPRKType mprktype)
1159 {
1160 PetscFunctionBegin;
1161 PetscValidHeaderSpecific(ts, TS_CLASSID, 1);
1162 PetscAssertPointer(mprktype, 2);
1163 PetscTryMethod(ts, "TSMPRKSetType_C", (TS, TSMPRKType), (ts, mprktype));
1164 PetscFunctionReturn(PETSC_SUCCESS);
1165 }
1166
1167 /*@
1168 TSMPRKGetType - Get the type of `TSMPRK` scheme
1169
1170 Not Collective
1171
1172 Input Parameter:
1173 . ts - timestepping context
1174
1175 Output Parameter:
1176 . mprktype - type of `TSMPRK` scheme
1177
1178 Level: intermediate
1179
1180 .seealso: [](ch_ts), `TSMPRK`
1181 @*/
TSMPRKGetType(TS ts,TSMPRKType * mprktype)1182 PetscErrorCode TSMPRKGetType(TS ts, TSMPRKType *mprktype)
1183 {
1184 PetscFunctionBegin;
1185 PetscValidHeaderSpecific(ts, TS_CLASSID, 1);
1186 PetscUseMethod(ts, "TSMPRKGetType_C", (TS, TSMPRKType *), (ts, mprktype));
1187 PetscFunctionReturn(PETSC_SUCCESS);
1188 }
1189
TSMPRKGetType_MPRK(TS ts,TSMPRKType * mprktype)1190 static PetscErrorCode TSMPRKGetType_MPRK(TS ts, TSMPRKType *mprktype)
1191 {
1192 TS_MPRK *mprk = (TS_MPRK *)ts->data;
1193
1194 PetscFunctionBegin;
1195 *mprktype = mprk->tableau->name;
1196 PetscFunctionReturn(PETSC_SUCCESS);
1197 }
1198
TSMPRKSetType_MPRK(TS ts,TSMPRKType mprktype)1199 static PetscErrorCode TSMPRKSetType_MPRK(TS ts, TSMPRKType mprktype)
1200 {
1201 TS_MPRK *mprk = (TS_MPRK *)ts->data;
1202 PetscBool match;
1203 MPRKTableauLink link;
1204
1205 PetscFunctionBegin;
1206 if (mprk->tableau) {
1207 PetscCall(PetscStrcmp(mprk->tableau->name, mprktype, &match));
1208 if (match) PetscFunctionReturn(PETSC_SUCCESS);
1209 }
1210 for (link = MPRKTableauList; link; link = link->next) {
1211 PetscCall(PetscStrcmp(link->tab.name, mprktype, &match));
1212 if (match) {
1213 if (ts->setupcalled) PetscCall(TSMPRKTableauReset(ts));
1214 mprk->tableau = &link->tab;
1215 if (ts->setupcalled) PetscCall(TSMPRKTableauSetUp(ts));
1216 PetscFunctionReturn(PETSC_SUCCESS);
1217 }
1218 }
1219 SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_UNKNOWN_TYPE, "Could not find '%s'", mprktype);
1220 }
1221
TSGetStages_MPRK(TS ts,PetscInt * ns,Vec ** Y)1222 static PetscErrorCode TSGetStages_MPRK(TS ts, PetscInt *ns, Vec **Y)
1223 {
1224 TS_MPRK *mprk = (TS_MPRK *)ts->data;
1225
1226 PetscFunctionBegin;
1227 *ns = mprk->tableau->s;
1228 if (Y) *Y = mprk->Y;
1229 PetscFunctionReturn(PETSC_SUCCESS);
1230 }
1231
TSDestroy_MPRK(TS ts)1232 static PetscErrorCode TSDestroy_MPRK(TS ts)
1233 {
1234 PetscFunctionBegin;
1235 PetscCall(TSReset_MPRK(ts));
1236 if (ts->dm) {
1237 PetscCall(DMCoarsenHookRemove(ts->dm, DMCoarsenHook_TSMPRK, DMRestrictHook_TSMPRK, ts));
1238 PetscCall(DMSubDomainHookRemove(ts->dm, DMSubDomainHook_TSMPRK, DMSubDomainRestrictHook_TSMPRK, ts));
1239 }
1240 PetscCall(PetscFree(ts->data));
1241 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSMPRKGetType_C", NULL));
1242 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSMPRKSetType_C", NULL));
1243 PetscFunctionReturn(PETSC_SUCCESS);
1244 }
1245
1246 /*MC
1247 TSMPRK - ODE solver using Multirate Partitioned Runge-Kutta schemes
1248
1249 The user should provide the right-hand side of the equation using `TSSetRHSFunction()`.
1250
1251 Level: beginner
1252
1253 Note:
1254 The default is `TSMPRKPM2`, it can be changed with `TSMPRKSetType()` or -ts_mprk_type
1255
1256 .seealso: [](ch_ts), `TSCreate()`, `TS`, `TSSetType()`, `TSMPRKSetType()`, `TSMPRKGetType()`, `TSMPRKType`, `TSMPRKRegister()`, `TSMPRKSetMultirateType()`
1257 `TSMPRKM2`, `TSMPRKM3`, `TSMPRKRFSMR3`, `TSMPRKRFSMR2`, `TSType`
1258 M*/
TSCreate_MPRK(TS ts)1259 PETSC_EXTERN PetscErrorCode TSCreate_MPRK(TS ts)
1260 {
1261 TS_MPRK *mprk;
1262
1263 PetscFunctionBegin;
1264 PetscCall(TSMPRKInitializePackage());
1265
1266 ts->ops->reset = TSReset_MPRK;
1267 ts->ops->destroy = TSDestroy_MPRK;
1268 ts->ops->view = TSView_MPRK;
1269 ts->ops->load = TSLoad_MPRK;
1270 ts->ops->setup = TSSetUp_MPRK;
1271 ts->ops->step = TSStep_MPRK;
1272 ts->ops->evaluatestep = TSEvaluateStep_MPRK;
1273 ts->ops->setfromoptions = TSSetFromOptions_MPRK;
1274 ts->ops->getstages = TSGetStages_MPRK;
1275
1276 PetscCall(PetscNew(&mprk));
1277 ts->data = (void *)mprk;
1278
1279 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSMPRKGetType_C", TSMPRKGetType_MPRK));
1280 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSMPRKSetType_C", TSMPRKSetType_MPRK));
1281
1282 PetscCall(TSMPRKSetType(ts, TSMPRKDefault));
1283 PetscFunctionReturn(PETSC_SUCCESS);
1284 }
1285