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