1 /*
2 Code for time stepping with the General Linear with Error Estimation method
3
4 Notes:
5 The general system is written as
6
7 Udot = F(t,U)
8
9 */
10 #include <petsc/private/tsimpl.h> /*I "petscts.h" I*/
11 #include <petscdm.h>
12
13 static PetscBool cited = PETSC_FALSE;
14 static const char citation[] = "@ARTICLE{Constantinescu_TR2016b,\n"
15 " author = {Constantinescu, E.M.},\n"
16 " title = {Estimating Global Errors in Time Stepping},\n"
17 " journal = {ArXiv e-prints},\n"
18 " year = 2016,\n"
19 " adsurl = {http://adsabs.harvard.edu/abs/2015arXiv150305166C}\n}\n";
20
21 static TSGLEEType TSGLEEDefaultType = TSGLEE35;
22 static PetscBool TSGLEERegisterAllCalled;
23 static PetscBool TSGLEEPackageInitialized;
24 static PetscInt explicit_stage_time_id;
25
26 typedef struct _GLEETableau *GLEETableau;
27 struct _GLEETableau {
28 char *name;
29 PetscInt order; /* Classical approximation order of the method i*/
30 PetscInt s; /* Number of stages */
31 PetscInt r; /* Number of steps */
32 PetscReal gamma; /* LTE ratio */
33 PetscReal *A, *B, *U, *V, *S, *F, *c; /* Tableau */
34 PetscReal *Fembed; /* Embedded final method coefficients */
35 PetscReal *Ferror; /* Coefficients for computing error */
36 PetscReal *Serror; /* Coefficients for initializing the error */
37 PetscInt pinterp; /* Interpolation order */
38 PetscReal *binterp; /* Interpolation coefficients */
39 PetscReal ccfl; /* Placeholder for CFL coefficient relative to forward Euler */
40 };
41 typedef struct _GLEETableauLink *GLEETableauLink;
42 struct _GLEETableauLink {
43 struct _GLEETableau tab;
44 GLEETableauLink next;
45 };
46 static GLEETableauLink GLEETableauList;
47
48 typedef struct {
49 GLEETableau tableau;
50 Vec *Y; /* Solution vector (along with auxiliary solution y~ or eps) */
51 Vec *X; /* Temporary solution vector */
52 Vec *YStage; /* Stage values */
53 Vec *YdotStage; /* Stage right-hand side */
54 Vec W; /* Right-hand-side for implicit stage solve */
55 Vec Ydot; /* Work vector holding Ydot during residual evaluation */
56 Vec yGErr; /* Vector holding the global error after a step is completed */
57 PetscScalar *swork; /* Scalar work (size of the number of stages)*/
58 PetscScalar *rwork; /* Scalar work (size of the number of steps)*/
59 PetscReal scoeff; /* shift = scoeff/dt */
60 PetscReal stage_time;
61 TSStepStatus status;
62 } TS_GLEE;
63
64 /*MC
65 TSGLEEi1 - Second order three stage implicit GLEE method
66
67 This method has two stages.
68 s = 3, r = 2
69
70 Level: advanced
71
72 .seealso: [](ch_ts), `TSGLEE`
73 M*/
74 /*MC
75 TSGLEE23 - Second order three stage explicit GLEE method
76
77 This method has three stages.
78 s = 3, r = 2
79
80 Level: advanced
81
82 .seealso: [](ch_ts), `TSGLEE`
83 M*/
84 /*MC
85 TSGLEE24 - Second order four stage explicit GLEE method
86
87 This method has four stages.
88 s = 4, r = 2
89
90 Level: advanced
91
92 .seealso: [](ch_ts), `TSGLEE`
93 M*/
94 /*MC
95 TSGLEE25i - Second order five stage explicit GLEE method
96
97 This method has five stages.
98 s = 5, r = 2
99
100 Level: advanced
101
102 .seealso: [](ch_ts), `TSGLEE`
103 M*/
104 /*MC
105 TSGLEE35 - Third order five stage explicit GLEE method
106
107 This method has five stages.
108 s = 5, r = 2
109
110 Level: advanced
111
112 .seealso: [](ch_ts), `TSGLEE`
113 M*/
114 /*MC
115 TSGLEEEXRK2A - Second order six stage explicit GLEE method
116
117 This method has six stages.
118 s = 6, r = 2
119
120 Level: advanced
121
122 .seealso: [](ch_ts), `TSGLEE`
123 M*/
124 /*MC
125 TSGLEERK32G1 - Third order eight stage explicit GLEE method
126
127 This method has eight stages.
128 s = 8, r = 2
129
130 Level: advanced
131
132 .seealso: [](ch_ts), `TSGLEE`
133 M*/
134 /*MC
135 TSGLEERK285EX - Second order nine stage explicit GLEE method
136
137 This method has nine stages.
138 s = 9, r = 2
139
140 Level: advanced
141
142 .seealso: [](ch_ts), `TSGLEE`
143 M*/
144
145 /*@C
146 TSGLEERegisterAll - Registers all of the General Linear with Error Estimation methods in `TSGLEE`
147
148 Not Collective, but should be called by all processes which will need the schemes to be registered
149
150 Level: advanced
151
152 .seealso: [](ch_ts), `TSGLEERegisterDestroy()`
153 @*/
TSGLEERegisterAll(void)154 PetscErrorCode TSGLEERegisterAll(void)
155 {
156 PetscFunctionBegin;
157 if (TSGLEERegisterAllCalled) PetscFunctionReturn(PETSC_SUCCESS);
158 TSGLEERegisterAllCalled = PETSC_TRUE;
159
160 {
161 #define GAMMA 0.5
162 /* y-eps form */
163 const PetscInt p = 1, s = 3, r = 2;
164 const PetscReal A[3][3] =
165 {
166 {1.0, 0, 0 },
167 {0, 0.5, 0 },
168 {0, 0.5, 0.5}
169 },
170 B[2][3] = {{1.0, 0, 0}, {-2.0, 1.0, 1.0}}, U[3][2] = {{1.0, 0}, {1.0, 0.5}, {1.0, 0.5}}, V[2][2] = {{1, 0}, {0, 1}}, S[2] = {1, 0}, F[2] = {1, 0}, Fembed[2] = {1, 1 - GAMMA}, Ferror[2] = {0, 1}, Serror[2] = {1, 0};
171 PetscCall(TSGLEERegister(TSGLEEi1, p, s, r, GAMMA, &A[0][0], &B[0][0], &U[0][0], &V[0][0], S, F, NULL, Fembed, Ferror, Serror, 0, NULL));
172 }
173 {
174 #undef GAMMA
175 #define GAMMA 0.0
176 /* y-eps form */
177 const PetscInt p = 2, s = 3, r = 2;
178 const PetscReal A[3][3] =
179 {
180 {0, 0, 0},
181 {1, 0, 0},
182 {0.25, 0.25, 0}
183 },
184 B[2][3] = {{1.0 / 12.0, 1.0 / 12.0, 5.0 / 6.0}, {1.0 / 12.0, 1.0 / 12.0, -1.0 / 6.0}}, U[3][2] = {{1, 0}, {1, 10}, {1, -1}}, V[2][2] = {{1, 0}, {0, 1}}, S[2] = {1, 0}, F[2] = {1, 0}, Fembed[2] = {1, 1 - GAMMA}, Ferror[2] = {0, 1}, Serror[2] = {1, 0};
185 PetscCall(TSGLEERegister(TSGLEE23, p, s, r, GAMMA, &A[0][0], &B[0][0], &U[0][0], &V[0][0], S, F, NULL, Fembed, Ferror, Serror, 0, NULL));
186 }
187 {
188 #undef GAMMA
189 #define GAMMA 0.0
190 /* y-y~ form */
191 const PetscInt p = 2, s = 4, r = 2;
192 const PetscReal A[4][4] =
193 {
194 {0, 0, 0, 0},
195 {0.75, 0, 0, 0},
196 {0.25, 29.0 / 60.0, 0, 0},
197 {-21.0 / 44.0, 145.0 / 44.0, -20.0 / 11.0, 0}
198 },
199 B[2][4] = {{109.0 / 275.0, 58.0 / 75.0, -37.0 / 110.0, 1.0 / 6.0}, {3.0 / 11.0, 0, 75.0 / 88.0, -1.0 / 8.0}}, U[4][2] = {{0, 1}, {75.0 / 58.0, -17.0 / 58.0}, {0, 1}, {0, 1}}, V[2][2] = {{1, 0}, {0, 1}}, S[2] = {1, 1}, F[2] = {1, 0}, Fembed[2] = {0, 1}, Ferror[2] = {-1.0 / (1.0 - GAMMA), 1.0 / (1.0 - GAMMA)}, Serror[2] = {1.0 - GAMMA, 1.0};
200 PetscCall(TSGLEERegister(TSGLEE24, p, s, r, GAMMA, &A[0][0], &B[0][0], &U[0][0], &V[0][0], S, F, NULL, Fembed, Ferror, Serror, 0, NULL));
201 }
202 {
203 #undef GAMMA
204 #define GAMMA 0.0
205 /* y-y~ form */
206 const PetscInt p = 2, s = 5, r = 2;
207 const PetscReal A[5][5] =
208 {
209 {0, 0, 0, 0, 0},
210 {-0.94079244066783383269, 0, 0, 0, 0},
211 {0.64228187778301907108, 0.10915356933958500042, 0, 0, 0},
212 {-0.51764297742287450812, 0.74414270351096040738, -0.71404164927824538121, 0, 0},
213 {-0.44696561556825969206, -0.76768425657590196518, 0.20111608138142987881, 0.93828186737840469796, 0}
214 },
215 B[2][5] = {{-0.029309178948150356153, -0.49671981884013874923, 0.34275801517650053274, 0.32941112623949194988, 0.85385985637229662276}, {0.78133219686062535272, 0.074238691892675897635, 0.57957363498384957966, -0.24638502829674959968, -0.18875949544040123033}}, U[5][2] = {{0.16911424754448327735, 0.83088575245551672265}, {0.53638465733199574340, 0.46361534266800425660}, {0.39901579167169582526, 0.60098420832830417474}, {0.87689005530618575480, 0.12310994469381424520}, {0.99056100455550913009, 0.0094389954444908699092}}, V[2][2] = {{1, 0}, {0, 1}}, S[2] = {1, 1}, F[2] = {1, 0}, Fembed[2] = {0, 1}, Ferror[2] = {-1.0 / (1.0 - GAMMA), 1.0 / (1.0 - GAMMA)}, Serror[2] = {1.0 - GAMMA, 1.0};
216 PetscCall(TSGLEERegister(TSGLEE25I, p, s, r, GAMMA, &A[0][0], &B[0][0], &U[0][0], &V[0][0], S, F, NULL, Fembed, Ferror, Serror, 0, NULL));
217 }
218 {
219 #undef GAMMA
220 #define GAMMA 0.0
221 /* y-y~ form */
222 const PetscInt p = 3, s = 5, r = 2;
223 const PetscReal A[5][5] =
224 {
225 {0, 0, 0, 0, 0},
226 {-2169604947363702313.0 / 24313474998937147335.0, 0, 0, 0, 0},
227 {46526746497697123895.0 / 94116917485856474137.0, -10297879244026594958.0 / 49199457603717988219.0, 0, 0, 0},
228 {23364788935845982499.0 / 87425311444725389446.0, -79205144337496116638.0 / 148994349441340815519.0, 40051189859317443782.0 / 36487615018004984309.0, 0, 0},
229 {42089522664062539205.0 / 124911313006412840286.0, -15074384760342762939.0 / 137927286865289746282.0, -62274678522253371016.0 / 125918573676298591413.0, 13755475729852471739.0 / 79257927066651693390.0, 0}
230 },
231 B[2][5] = {{61546696837458703723.0 / 56982519523786160813.0, -55810892792806293355.0 / 206957624151308356511.0, 24061048952676379087.0 / 158739347956038723465.0, 3577972206874351339.0 / 7599733370677197135.0, -59449832954780563947.0 / 137360038685338563670.0}, {-9738262186984159168.0 / 99299082461487742983.0, -32797097931948613195.0 / 61521565616362163366.0, 42895514606418420631.0 / 71714201188501437336.0, 22608567633166065068.0 / 55371917805607957003.0, 94655809487476459565.0 / 151517167160302729021.0}}, U[5][2] = {{70820309139834661559.0 / 80863923579509469826.0, 10043614439674808267.0 / 80863923579509469826.0}, {161694774978034105510.0 / 106187653640211060371.0, -55507121337823045139.0 / 106187653640211060371.0}, {78486094644566264568.0 / 88171030896733822981.0, 9684936252167558413.0 / 88171030896733822981.0}, {65394922146334854435.0 / 84570853840405479554.0, 19175931694070625119.0 / 84570853840405479554.0}, {8607282770183754108.0 / 108658046436496925911.0, 100050763666313171803.0 / 108658046436496925911.0}}, V[2][2] = {{1, 0}, {0, 1}}, S[2] = {1, 1}, F[2] = {1, 0}, Fembed[2] = {0, 1}, Ferror[2] = {-1.0 / (1.0 - GAMMA), 1.0 / (1.0 - GAMMA)}, Serror[2] = {1.0 - GAMMA, 1.0};
232 PetscCall(TSGLEERegister(TSGLEE35, p, s, r, GAMMA, &A[0][0], &B[0][0], &U[0][0], &V[0][0], S, F, NULL, Fembed, Ferror, Serror, 0, NULL));
233 }
234 {
235 #undef GAMMA
236 #define GAMMA 0.25
237 /* y-eps form */
238 const PetscInt p = 2, s = 6, r = 2;
239 const PetscReal A[6][6] =
240 {
241 {0, 0, 0, 0, 0, 0},
242 {1, 0, 0, 0, 0, 0},
243 {0, 0, 0, 0, 0, 0},
244 {0, 0, 0.5, 0, 0, 0},
245 {0, 0, 0.25, 0.25, 0, 0},
246 {0, 0, 0.25, 0.25, 0.5, 0}
247 },
248 B[2][6] = {{0.5, 0.5, 0, 0, 0, 0}, {-2.0 / 3.0, -2.0 / 3.0, 1.0 / 3.0, 1.0 / 3.0, 1.0 / 3.0, 1.0 / 3.0}}, U[6][2] = {{1, 0}, {1, 0}, {1, 0.75}, {1, 0.75}, {1, 0.75}, {1, 0.75}}, V[2][2] = {{1, 0}, {0, 1}}, S[2] = {1, 0}, F[2] = {1, 0}, Fembed[2] = {1, 1 - GAMMA}, Ferror[2] = {0, 1}, Serror[2] = {1, 0};
249 PetscCall(TSGLEERegister(TSGLEEEXRK2A, p, s, r, GAMMA, &A[0][0], &B[0][0], &U[0][0], &V[0][0], S, F, NULL, Fembed, Ferror, Serror, 0, NULL));
250 }
251 {
252 #undef GAMMA
253 #define GAMMA 0.0
254 /* y-eps form */
255 const PetscInt p = 3, s = 8, r = 2;
256 const PetscReal A[8][8] =
257 {
258 {0, 0, 0, 0, 0, 0, 0, 0},
259 {0.5, 0, 0, 0, 0, 0, 0, 0},
260 {-1, 2, 0, 0, 0, 0, 0, 0},
261 {1.0 / 6.0, 2.0 / 3.0, 1.0 / 6.0, 0, 0, 0, 0, 0},
262 {0, 0, 0, 0, 0, 0, 0, 0},
263 {-7.0 / 24.0, 1.0 / 3.0, 1.0 / 12.0, -1.0 / 8.0, 0.5, 0, 0, 0},
264 {7.0 / 6.0, -4.0 / 3.0, -1.0 / 3.0, 0.5, -1.0, 2.0, 0, 0},
265 {0, 0, 0, 0, 1.0 / 6.0, 2.0 / 3.0, 1.0 / 6.0, 0}
266 },
267 B[2][8] = {{1.0 / 6.0, 2.0 / 3.0, 1.0 / 6.0, 0, 0, 0, 0, 0}, {-1.0 / 6.0, -2.0 / 3.0, -1.0 / 6.0, 0, 1.0 / 6.0, 2.0 / 3.0, 1.0 / 6.0, 0}}, U[8][2] = {{1, 0}, {1, 0}, {1, 0}, {1, 0}, {1, 1}, {1, 1}, {1, 1}, {1, 1}}, V[2][2] = {{1, 0}, {0, 1}}, S[2] = {1, 0}, F[2] = {1, 0}, Fembed[2] = {1, 1 - GAMMA}, Ferror[2] = {0, 1}, Serror[2] = {1, 0};
268 PetscCall(TSGLEERegister(TSGLEERK32G1, p, s, r, GAMMA, &A[0][0], &B[0][0], &U[0][0], &V[0][0], S, F, NULL, Fembed, Ferror, Serror, 0, NULL));
269 }
270 {
271 #undef GAMMA
272 #define GAMMA 0.25
273 /* y-eps form */
274 const PetscInt p = 2, s = 9, r = 2;
275 const PetscReal A[9][9] =
276 {
277 {0, 0, 0, 0, 0, 0, 0, 0, 0},
278 {0.585786437626904966, 0, 0, 0, 0, 0, 0, 0, 0},
279 {0.149999999999999994, 0.849999999999999978, 0, 0, 0, 0, 0, 0, 0},
280 {0, 0, 0, 0, 0, 0, 0, 0, 0},
281 {0, 0, 0, 0.292893218813452483, 0, 0, 0, 0, 0},
282 {0, 0, 0, 0.0749999999999999972, 0.424999999999999989, 0, 0, 0, 0},
283 {0, 0, 0, 0.176776695296636893, 0.176776695296636893, 0.146446609406726241, 0, 0, 0},
284 {0, 0, 0, 0.176776695296636893, 0.176776695296636893, 0.146446609406726241, 0.292893218813452483, 0, 0},
285 {0, 0, 0, 0.176776695296636893, 0.176776695296636893, 0.146446609406726241, 0.0749999999999999972, 0.424999999999999989, 0}
286 },
287 B[2][9] = {{0.353553390593273786, 0.353553390593273786, 0.292893218813452483, 0, 0, 0, 0, 0, 0}, {-0.471404520791031678, -0.471404520791031678, -0.390524291751269959, 0.235702260395515839, 0.235702260395515839, 0.195262145875634979, 0.235702260395515839, 0.235702260395515839, 0.195262145875634979}}, U[9][2] = {{1, 0}, {1, 0}, {1, 0}, {1, 0.75}, {1, 0.75}, {1, 0.75}, {1, 0.75}, {1, 0.75}, {1, 0.75}}, V[2][2] = {{1, 0}, {0, 1}}, S[2] = {1, 0}, F[2] = {1, 0}, Fembed[2] = {1, 1 - GAMMA}, Ferror[2] = {0, 1}, Serror[2] = {1, 0};
288 PetscCall(TSGLEERegister(TSGLEERK285EX, p, s, r, GAMMA, &A[0][0], &B[0][0], &U[0][0], &V[0][0], S, F, NULL, Fembed, Ferror, Serror, 0, NULL));
289 }
290 PetscFunctionReturn(PETSC_SUCCESS);
291 }
292
293 /*@C
294 TSGLEERegisterDestroy - Frees the list of schemes that were registered by `TSGLEERegister()`.
295
296 Not Collective
297
298 Level: advanced
299
300 .seealso: [](ch_ts), `TSGLEERegister()`, `TSGLEERegisterAll()`
301 @*/
TSGLEERegisterDestroy(void)302 PetscErrorCode TSGLEERegisterDestroy(void)
303 {
304 GLEETableauLink link;
305
306 PetscFunctionBegin;
307 while ((link = GLEETableauList)) {
308 GLEETableau t = &link->tab;
309 GLEETableauList = link->next;
310 PetscCall(PetscFree5(t->A, t->B, t->U, t->V, t->c));
311 PetscCall(PetscFree2(t->S, t->F));
312 PetscCall(PetscFree(t->Fembed));
313 PetscCall(PetscFree(t->Ferror));
314 PetscCall(PetscFree(t->Serror));
315 PetscCall(PetscFree(t->binterp));
316 PetscCall(PetscFree(t->name));
317 PetscCall(PetscFree(link));
318 }
319 TSGLEERegisterAllCalled = PETSC_FALSE;
320 PetscFunctionReturn(PETSC_SUCCESS);
321 }
322
323 /*@C
324 TSGLEEInitializePackage - This function initializes everything in the `TSGLEE` package. It is called
325 from `TSInitializePackage()`.
326
327 Level: developer
328
329 .seealso: [](ch_ts), `PetscInitialize()`
330 @*/
TSGLEEInitializePackage(void)331 PetscErrorCode TSGLEEInitializePackage(void)
332 {
333 PetscFunctionBegin;
334 if (TSGLEEPackageInitialized) PetscFunctionReturn(PETSC_SUCCESS);
335 TSGLEEPackageInitialized = PETSC_TRUE;
336 PetscCall(TSGLEERegisterAll());
337 PetscCall(PetscObjectComposedDataRegister(&explicit_stage_time_id));
338 PetscCall(PetscRegisterFinalize(TSGLEEFinalizePackage));
339 PetscFunctionReturn(PETSC_SUCCESS);
340 }
341
342 /*@C
343 TSGLEEFinalizePackage - This function destroys everything in the `TSGLEE` package. It is
344 called from `PetscFinalize()`.
345
346 Level: developer
347
348 .seealso: [](ch_ts), `PetscFinalize()`
349 @*/
TSGLEEFinalizePackage(void)350 PetscErrorCode TSGLEEFinalizePackage(void)
351 {
352 PetscFunctionBegin;
353 TSGLEEPackageInitialized = PETSC_FALSE;
354 PetscCall(TSGLEERegisterDestroy());
355 PetscFunctionReturn(PETSC_SUCCESS);
356 }
357
358 /*@C
359 TSGLEERegister - register a new `TSGLEE` scheme by providing the entries in the Butcher tableau
360
361 Not Collective, but the same schemes should be registered on all processes on which they will be used, No Fortran Support
362
363 Input Parameters:
364 + name - identifier for method
365 . order - order of method
366 . s - number of stages
367 . r - number of steps
368 . gamma - LTE ratio
369 . A - stage coefficients (dimension s*s, row-major)
370 . B - step completion coefficients (dimension r*s, row-major)
371 . U - method coefficients (dimension s*r, row-major)
372 . V - method coefficients (dimension r*r, row-major)
373 . S - starting coefficients
374 . F - finishing coefficients
375 . c - abscissa (dimension s; NULL to use row sums of A)
376 . Fembed - step completion coefficients for embedded method
377 . Ferror - error computation coefficients
378 . Serror - error initialization coefficients
379 . pinterp - order of interpolation (0 if unavailable)
380 - binterp - array of interpolation coefficients (NULL if unavailable)
381
382 Level: advanced
383
384 Note:
385 Several `TSGLEE` methods are provided, this function is only needed to create new methods.
386
387 .seealso: [](ch_ts), `TSGLEE`
388 @*/
TSGLEERegister(TSGLEEType name,PetscInt order,PetscInt s,PetscInt r,PetscReal gamma,const PetscReal A[],const PetscReal B[],const PetscReal U[],const PetscReal V[],const PetscReal S[],const PetscReal F[],const PetscReal c[],const PetscReal Fembed[],const PetscReal Ferror[],const PetscReal Serror[],PetscInt pinterp,const PetscReal binterp[])389 PetscErrorCode TSGLEERegister(TSGLEEType name, PetscInt order, PetscInt s, PetscInt r, PetscReal gamma, const PetscReal A[], const PetscReal B[], const PetscReal U[], const PetscReal V[], const PetscReal S[], const PetscReal F[], const PetscReal c[], const PetscReal Fembed[], const PetscReal Ferror[], const PetscReal Serror[], PetscInt pinterp, const PetscReal binterp[])
390 {
391 GLEETableauLink link;
392 GLEETableau t;
393 PetscInt i, j;
394
395 PetscFunctionBegin;
396 PetscCall(TSGLEEInitializePackage());
397 PetscCall(PetscNew(&link));
398 t = &link->tab;
399 PetscCall(PetscStrallocpy(name, &t->name));
400 t->order = order;
401 t->s = s;
402 t->r = r;
403 t->gamma = gamma;
404 PetscCall(PetscMalloc5(s * s, &t->A, r * r, &t->V, s, &t->c, r * s, &t->B, s * r, &t->U));
405 PetscCall(PetscMalloc2(r, &t->S, r, &t->F));
406 PetscCall(PetscArraycpy(t->A, A, s * s));
407 PetscCall(PetscArraycpy(t->B, B, r * s));
408 PetscCall(PetscArraycpy(t->U, U, s * r));
409 PetscCall(PetscArraycpy(t->V, V, r * r));
410 PetscCall(PetscArraycpy(t->S, S, r));
411 PetscCall(PetscArraycpy(t->F, F, r));
412 if (c) {
413 PetscCall(PetscArraycpy(t->c, c, s));
414 } else {
415 for (i = 0; i < s; i++)
416 for (j = 0, t->c[i] = 0; j < s; j++) t->c[i] += A[i * s + j];
417 }
418 PetscCall(PetscMalloc1(r, &t->Fembed));
419 PetscCall(PetscMalloc1(r, &t->Ferror));
420 PetscCall(PetscMalloc1(r, &t->Serror));
421 PetscCall(PetscArraycpy(t->Fembed, Fembed, r));
422 PetscCall(PetscArraycpy(t->Ferror, Ferror, r));
423 PetscCall(PetscArraycpy(t->Serror, Serror, r));
424 t->pinterp = pinterp;
425 PetscCall(PetscMalloc1(s * pinterp, &t->binterp));
426 PetscCall(PetscArraycpy(t->binterp, binterp, s * pinterp));
427
428 link->next = GLEETableauList;
429 GLEETableauList = link;
430 PetscFunctionReturn(PETSC_SUCCESS);
431 }
432
TSEvaluateStep_GLEE(TS ts,PetscInt order,Vec X,PetscBool * done)433 static PetscErrorCode TSEvaluateStep_GLEE(TS ts, PetscInt order, Vec X, PetscBool *done)
434 {
435 TS_GLEE *glee = (TS_GLEE *)ts->data;
436 GLEETableau tab = glee->tableau;
437 PetscReal h, *B = tab->B, *V = tab->V, *F = tab->F, *Fembed = tab->Fembed;
438 PetscInt s = tab->s, r = tab->r, i, j;
439 Vec *Y = glee->Y, *YdotStage = glee->YdotStage;
440 PetscScalar *ws = glee->swork, *wr = glee->rwork;
441
442 PetscFunctionBegin;
443 switch (glee->status) {
444 case TS_STEP_INCOMPLETE:
445 case TS_STEP_PENDING:
446 h = ts->time_step;
447 break;
448 case TS_STEP_COMPLETE:
449 h = ts->ptime - ts->ptime_prev;
450 break;
451 default:
452 SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_PLIB, "Invalid TSStepStatus");
453 }
454
455 if (order == tab->order) {
456 /* Note: Irrespective of whether status is TS_STEP_INCOMPLETE
457 or TS_STEP_COMPLETE, glee->X has the solution at the
458 beginning of the time step. So no need to roll-back.
459 */
460 if (glee->status == TS_STEP_INCOMPLETE) {
461 for (i = 0; i < r; i++) {
462 PetscCall(VecZeroEntries(Y[i]));
463 for (j = 0; j < r; j++) wr[j] = V[i * r + j];
464 PetscCall(VecMAXPY(Y[i], r, wr, glee->X));
465 for (j = 0; j < s; j++) ws[j] = h * B[i * s + j];
466 PetscCall(VecMAXPY(Y[i], s, ws, YdotStage));
467 }
468 PetscCall(VecZeroEntries(X));
469 for (j = 0; j < r; j++) wr[j] = F[j];
470 PetscCall(VecMAXPY(X, r, wr, Y));
471 } else PetscCall(VecCopy(ts->vec_sol, X));
472 PetscFunctionReturn(PETSC_SUCCESS);
473
474 } else if (order == tab->order - 1) {
475 /* Complete with the embedded method (Fembed) */
476 for (i = 0; i < r; i++) {
477 PetscCall(VecZeroEntries(Y[i]));
478 for (j = 0; j < r; j++) wr[j] = V[i * r + j];
479 PetscCall(VecMAXPY(Y[i], r, wr, glee->X));
480 for (j = 0; j < s; j++) ws[j] = h * B[i * s + j];
481 PetscCall(VecMAXPY(Y[i], s, ws, YdotStage));
482 }
483 PetscCall(VecZeroEntries(X));
484 for (j = 0; j < r; j++) wr[j] = Fembed[j];
485 PetscCall(VecMAXPY(X, r, wr, Y));
486
487 if (done) *done = PETSC_TRUE;
488 PetscFunctionReturn(PETSC_SUCCESS);
489 }
490 PetscCheck(done, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "GLEE '%s' of order %" PetscInt_FMT " cannot evaluate step at order %" PetscInt_FMT, tab->name, tab->order, order);
491 *done = PETSC_FALSE;
492 PetscFunctionReturn(PETSC_SUCCESS);
493 }
494
TSStep_GLEE(TS ts)495 static PetscErrorCode TSStep_GLEE(TS ts)
496 {
497 TS_GLEE *glee = (TS_GLEE *)ts->data;
498 GLEETableau tab = glee->tableau;
499 const PetscInt s = tab->s, r = tab->r;
500 PetscReal *A = tab->A, *U = tab->U, *F = tab->F, *c = tab->c;
501 Vec *Y = glee->Y, *X = glee->X, *YStage = glee->YStage, *YdotStage = glee->YdotStage, W = glee->W;
502 SNES snes;
503 PetscScalar *ws = glee->swork, *wr = glee->rwork;
504 TSAdapt adapt;
505 PetscInt i, j, reject, next_scheme, its, lits;
506 PetscReal next_time_step;
507 PetscReal t;
508 PetscBool accept;
509
510 PetscFunctionBegin;
511 PetscCall(PetscCitationsRegister(citation, &cited));
512
513 for (i = 0; i < r; i++) PetscCall(VecCopy(Y[i], X[i]));
514
515 PetscCall(TSGetSNES(ts, &snes));
516 next_time_step = ts->time_step;
517 t = ts->ptime;
518 accept = PETSC_TRUE;
519 glee->status = TS_STEP_INCOMPLETE;
520
521 for (reject = 0; reject < ts->max_reject && !ts->reason; reject++, ts->reject++) {
522 PetscReal h = ts->time_step;
523 PetscCall(TSPreStep(ts));
524
525 for (i = 0; i < s; i++) {
526 glee->stage_time = t + h * c[i];
527 PetscCall(TSPreStage(ts, glee->stage_time));
528
529 if (A[i * s + i] == 0) { /* Explicit stage */
530 PetscCall(VecZeroEntries(YStage[i]));
531 for (j = 0; j < r; j++) wr[j] = U[i * r + j];
532 PetscCall(VecMAXPY(YStage[i], r, wr, X));
533 for (j = 0; j < i; j++) ws[j] = h * A[i * s + j];
534 PetscCall(VecMAXPY(YStage[i], i, ws, YdotStage));
535 } else { /* Implicit stage */
536 glee->scoeff = 1.0 / A[i * s + i];
537 /* compute right-hand side */
538 PetscCall(VecZeroEntries(W));
539 for (j = 0; j < r; j++) wr[j] = U[i * r + j];
540 PetscCall(VecMAXPY(W, r, wr, X));
541 for (j = 0; j < i; j++) ws[j] = h * A[i * s + j];
542 PetscCall(VecMAXPY(W, i, ws, YdotStage));
543 PetscCall(VecScale(W, glee->scoeff / h));
544 /* set initial guess */
545 PetscCall(VecCopy(i > 0 ? YStage[i - 1] : ts->vec_sol, YStage[i]));
546 /* solve for this stage */
547 PetscCall(SNESSolve(snes, W, YStage[i]));
548 PetscCall(SNESGetIterationNumber(snes, &its));
549 PetscCall(SNESGetLinearSolveIterations(snes, &lits));
550 ts->snes_its += its;
551 ts->ksp_its += lits;
552 }
553 PetscCall(TSGetAdapt(ts, &adapt));
554 PetscCall(TSAdaptCheckStage(adapt, ts, glee->stage_time, YStage[i], &accept));
555 if (!accept) goto reject_step;
556 PetscCall(TSPostStage(ts, glee->stage_time, i, YStage));
557 PetscCall(TSComputeRHSFunction(ts, t + h * c[i], YStage[i], YdotStage[i]));
558 }
559 PetscCall(TSEvaluateStep(ts, tab->order, ts->vec_sol, NULL));
560 glee->status = TS_STEP_PENDING;
561
562 /* Register only the current method as a candidate because we're not supporting multiple candidates yet. */
563 PetscCall(TSGetAdapt(ts, &adapt));
564 PetscCall(TSAdaptCandidatesClear(adapt));
565 PetscCall(TSAdaptCandidateAdd(adapt, tab->name, tab->order, 1, tab->ccfl, (PetscReal)tab->s, PETSC_TRUE));
566 PetscCall(TSAdaptChoose(adapt, ts, ts->time_step, &next_scheme, &next_time_step, &accept));
567 if (accept) {
568 /* ignore next_scheme for now */
569 ts->ptime += ts->time_step;
570 ts->time_step = next_time_step;
571 glee->status = TS_STEP_COMPLETE;
572 /* compute and store the global error */
573 /* Note: this is not needed if TSAdaptGLEE is not used */
574 PetscCall(TSGetTimeError(ts, 0, &glee->yGErr));
575 PetscCall(PetscObjectComposedDataSetReal((PetscObject)ts->vec_sol, explicit_stage_time_id, ts->ptime));
576 break;
577 } else { /* Roll back the current step */
578 for (j = 0; j < r; j++) wr[j] = F[j];
579 PetscCall(VecMAXPY(ts->vec_sol, r, wr, X));
580 ts->time_step = next_time_step;
581 glee->status = TS_STEP_INCOMPLETE;
582 }
583 reject_step:
584 continue;
585 }
586 if (glee->status != TS_STEP_COMPLETE && !ts->reason) ts->reason = TS_DIVERGED_STEP_REJECTED;
587 PetscFunctionReturn(PETSC_SUCCESS);
588 }
589
TSInterpolate_GLEE(TS ts,PetscReal itime,Vec X)590 static PetscErrorCode TSInterpolate_GLEE(TS ts, PetscReal itime, Vec X)
591 {
592 TS_GLEE *glee = (TS_GLEE *)ts->data;
593 PetscInt s = glee->tableau->s, pinterp = glee->tableau->pinterp, i, j;
594 PetscReal h, tt, t;
595 PetscScalar *b;
596 const PetscReal *B = glee->tableau->binterp;
597
598 PetscFunctionBegin;
599 PetscCheck(B, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "TSGLEE %s does not have an interpolation formula", glee->tableau->name);
600 switch (glee->status) {
601 case TS_STEP_INCOMPLETE:
602 case TS_STEP_PENDING:
603 h = ts->time_step;
604 t = (itime - ts->ptime) / h;
605 break;
606 case TS_STEP_COMPLETE:
607 h = ts->ptime - ts->ptime_prev;
608 t = (itime - ts->ptime) / h + 1; /* In the interval [0,1] */
609 break;
610 default:
611 SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_PLIB, "Invalid TSStepStatus");
612 }
613 PetscCall(PetscMalloc1(s, &b));
614 for (i = 0; i < s; i++) b[i] = 0;
615 for (j = 0, tt = t; j < pinterp; j++, tt *= t) {
616 for (i = 0; i < s; i++) b[i] += h * B[i * pinterp + j] * tt;
617 }
618 PetscCall(VecCopy(glee->YStage[0], X));
619 PetscCall(VecMAXPY(X, s, b, glee->YdotStage));
620 PetscCall(PetscFree(b));
621 PetscFunctionReturn(PETSC_SUCCESS);
622 }
623
TSReset_GLEE(TS ts)624 static PetscErrorCode TSReset_GLEE(TS ts)
625 {
626 TS_GLEE *glee = (TS_GLEE *)ts->data;
627 PetscInt s, r;
628
629 PetscFunctionBegin;
630 if (!glee->tableau) PetscFunctionReturn(PETSC_SUCCESS);
631 s = glee->tableau->s;
632 r = glee->tableau->r;
633 PetscCall(VecDestroyVecs(r, &glee->Y));
634 PetscCall(VecDestroyVecs(r, &glee->X));
635 PetscCall(VecDestroyVecs(s, &glee->YStage));
636 PetscCall(VecDestroyVecs(s, &glee->YdotStage));
637 PetscCall(VecDestroy(&glee->Ydot));
638 PetscCall(VecDestroy(&glee->yGErr));
639 PetscCall(VecDestroy(&glee->W));
640 PetscCall(PetscFree2(glee->swork, glee->rwork));
641 PetscFunctionReturn(PETSC_SUCCESS);
642 }
643
TSGLEEGetVecs(TS ts,DM dm,Vec * Ydot)644 static PetscErrorCode TSGLEEGetVecs(TS ts, DM dm, Vec *Ydot)
645 {
646 TS_GLEE *glee = (TS_GLEE *)ts->data;
647
648 PetscFunctionBegin;
649 if (Ydot) {
650 if (dm && dm != ts->dm) PetscCall(DMGetNamedGlobalVector(dm, "TSGLEE_Ydot", Ydot));
651 else *Ydot = glee->Ydot;
652 }
653 PetscFunctionReturn(PETSC_SUCCESS);
654 }
655
TSGLEERestoreVecs(TS ts,DM dm,Vec * Ydot)656 static PetscErrorCode TSGLEERestoreVecs(TS ts, DM dm, Vec *Ydot)
657 {
658 PetscFunctionBegin;
659 if (Ydot) {
660 if (dm && dm != ts->dm) PetscCall(DMRestoreNamedGlobalVector(dm, "TSGLEE_Ydot", Ydot));
661 }
662 PetscFunctionReturn(PETSC_SUCCESS);
663 }
664
665 /*
666 This defines the nonlinear equation that is to be solved with SNES
667 */
SNESTSFormFunction_GLEE(SNES snes,Vec X,Vec F,TS ts)668 static PetscErrorCode SNESTSFormFunction_GLEE(SNES snes, Vec X, Vec F, TS ts)
669 {
670 TS_GLEE *glee = (TS_GLEE *)ts->data;
671 DM dm, dmsave;
672 Vec Ydot;
673 PetscReal shift = glee->scoeff / ts->time_step;
674
675 PetscFunctionBegin;
676 PetscCall(SNESGetDM(snes, &dm));
677 PetscCall(TSGLEEGetVecs(ts, dm, &Ydot));
678 /* Set Ydot = shift*X */
679 PetscCall(VecCopy(X, Ydot));
680 PetscCall(VecScale(Ydot, shift));
681 dmsave = ts->dm;
682 ts->dm = dm;
683
684 PetscCall(TSComputeIFunction(ts, glee->stage_time, X, Ydot, F, PETSC_FALSE));
685
686 ts->dm = dmsave;
687 PetscCall(TSGLEERestoreVecs(ts, dm, &Ydot));
688 PetscFunctionReturn(PETSC_SUCCESS);
689 }
690
SNESTSFormJacobian_GLEE(SNES snes,Vec X,Mat A,Mat B,TS ts)691 static PetscErrorCode SNESTSFormJacobian_GLEE(SNES snes, Vec X, Mat A, Mat B, TS ts)
692 {
693 TS_GLEE *glee = (TS_GLEE *)ts->data;
694 DM dm, dmsave;
695 Vec Ydot;
696 PetscReal shift = glee->scoeff / ts->time_step;
697
698 PetscFunctionBegin;
699 PetscCall(SNESGetDM(snes, &dm));
700 PetscCall(TSGLEEGetVecs(ts, dm, &Ydot));
701 /* glee->Ydot has already been computed in SNESTSFormFunction_GLEE (SNES guarantees this) */
702 dmsave = ts->dm;
703 ts->dm = dm;
704
705 PetscCall(TSComputeIJacobian(ts, glee->stage_time, X, Ydot, shift, A, B, PETSC_FALSE));
706
707 ts->dm = dmsave;
708 PetscCall(TSGLEERestoreVecs(ts, dm, &Ydot));
709 PetscFunctionReturn(PETSC_SUCCESS);
710 }
711
DMCoarsenHook_TSGLEE(DM fine,DM coarse,PetscCtx ctx)712 static PetscErrorCode DMCoarsenHook_TSGLEE(DM fine, DM coarse, PetscCtx ctx)
713 {
714 PetscFunctionBegin;
715 PetscFunctionReturn(PETSC_SUCCESS);
716 }
717
DMRestrictHook_TSGLEE(DM fine,Mat restrct,Vec rscale,Mat inject,DM coarse,PetscCtx ctx)718 static PetscErrorCode DMRestrictHook_TSGLEE(DM fine, Mat restrct, Vec rscale, Mat inject, DM coarse, PetscCtx ctx)
719 {
720 PetscFunctionBegin;
721 PetscFunctionReturn(PETSC_SUCCESS);
722 }
723
DMSubDomainHook_TSGLEE(DM dm,DM subdm,PetscCtx ctx)724 static PetscErrorCode DMSubDomainHook_TSGLEE(DM dm, DM subdm, PetscCtx ctx)
725 {
726 PetscFunctionBegin;
727 PetscFunctionReturn(PETSC_SUCCESS);
728 }
729
DMSubDomainRestrictHook_TSGLEE(DM dm,VecScatter gscat,VecScatter lscat,DM subdm,PetscCtx ctx)730 static PetscErrorCode DMSubDomainRestrictHook_TSGLEE(DM dm, VecScatter gscat, VecScatter lscat, DM subdm, PetscCtx ctx)
731 {
732 PetscFunctionBegin;
733 PetscFunctionReturn(PETSC_SUCCESS);
734 }
735
TSSetUp_GLEE(TS ts)736 static PetscErrorCode TSSetUp_GLEE(TS ts)
737 {
738 TS_GLEE *glee = (TS_GLEE *)ts->data;
739 GLEETableau tab;
740 PetscInt s, r;
741 DM dm;
742
743 PetscFunctionBegin;
744 if (!glee->tableau) PetscCall(TSGLEESetType(ts, TSGLEEDefaultType));
745 tab = glee->tableau;
746 s = tab->s;
747 r = tab->r;
748 PetscCall(VecDuplicateVecs(ts->vec_sol, r, &glee->Y));
749 PetscCall(VecDuplicateVecs(ts->vec_sol, r, &glee->X));
750 PetscCall(VecDuplicateVecs(ts->vec_sol, s, &glee->YStage));
751 PetscCall(VecDuplicateVecs(ts->vec_sol, s, &glee->YdotStage));
752 PetscCall(VecDuplicate(ts->vec_sol, &glee->Ydot));
753 PetscCall(VecDuplicate(ts->vec_sol, &glee->yGErr));
754 PetscCall(VecZeroEntries(glee->yGErr));
755 PetscCall(VecDuplicate(ts->vec_sol, &glee->W));
756 PetscCall(PetscMalloc2(s, &glee->swork, r, &glee->rwork));
757 PetscCall(TSGetDM(ts, &dm));
758 PetscCall(DMCoarsenHookAdd(dm, DMCoarsenHook_TSGLEE, DMRestrictHook_TSGLEE, ts));
759 PetscCall(DMSubDomainHookAdd(dm, DMSubDomainHook_TSGLEE, DMSubDomainRestrictHook_TSGLEE, ts));
760 PetscFunctionReturn(PETSC_SUCCESS);
761 }
762
TSStartingMethod_GLEE(TS ts)763 static PetscErrorCode TSStartingMethod_GLEE(TS ts)
764 {
765 TS_GLEE *glee = (TS_GLEE *)ts->data;
766 GLEETableau tab = glee->tableau;
767 PetscInt r = tab->r, i;
768 PetscReal *S = tab->S;
769
770 PetscFunctionBegin;
771 for (i = 0; i < r; i++) {
772 PetscCall(VecZeroEntries(glee->Y[i]));
773 PetscCall(VecAXPY(glee->Y[i], S[i], ts->vec_sol));
774 }
775 PetscFunctionReturn(PETSC_SUCCESS);
776 }
777
TSSetFromOptions_GLEE(TS ts,PetscOptionItems PetscOptionsObject)778 static PetscErrorCode TSSetFromOptions_GLEE(TS ts, PetscOptionItems PetscOptionsObject)
779 {
780 char gleetype[256];
781
782 PetscFunctionBegin;
783 PetscOptionsHeadBegin(PetscOptionsObject, "GLEE ODE solver options");
784 {
785 GLEETableauLink link;
786 PetscInt count, choice;
787 PetscBool flg;
788 const char **namelist;
789
790 PetscCall(PetscStrncpy(gleetype, TSGLEEDefaultType, sizeof(gleetype)));
791 for (link = GLEETableauList, count = 0; link; link = link->next, count++);
792 PetscCall(PetscMalloc1(count, (char ***)&namelist));
793 for (link = GLEETableauList, count = 0; link; link = link->next, count++) namelist[count] = link->tab.name;
794 PetscCall(PetscOptionsEList("-ts_glee_type", "Family of GLEE method", "TSGLEESetType", (const char *const *)namelist, count, gleetype, &choice, &flg));
795 PetscCall(TSGLEESetType(ts, flg ? namelist[choice] : gleetype));
796 PetscCall(PetscFree(namelist));
797 }
798 PetscOptionsHeadEnd();
799 PetscFunctionReturn(PETSC_SUCCESS);
800 }
801
TSView_GLEE(TS ts,PetscViewer viewer)802 static PetscErrorCode TSView_GLEE(TS ts, PetscViewer viewer)
803 {
804 TS_GLEE *glee = (TS_GLEE *)ts->data;
805 GLEETableau tab = glee->tableau;
806 PetscBool isascii;
807
808 PetscFunctionBegin;
809 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
810 if (isascii) {
811 TSGLEEType gleetype;
812 char buf[512];
813 PetscCall(TSGLEEGetType(ts, &gleetype));
814 PetscCall(PetscViewerASCIIPrintf(viewer, " GLEE type %s\n", gleetype));
815 PetscCall(PetscFormatRealArray(buf, sizeof(buf), "% 8.6f", tab->s, tab->c));
816 PetscCall(PetscViewerASCIIPrintf(viewer, " Abscissa c = %s\n", buf));
817 /* Note: print out r as well */
818 }
819 PetscFunctionReturn(PETSC_SUCCESS);
820 }
821
TSLoad_GLEE(TS ts,PetscViewer viewer)822 static PetscErrorCode TSLoad_GLEE(TS ts, PetscViewer viewer)
823 {
824 SNES snes;
825 TSAdapt tsadapt;
826
827 PetscFunctionBegin;
828 PetscCall(TSGetAdapt(ts, &tsadapt));
829 PetscCall(TSAdaptLoad(tsadapt, viewer));
830 PetscCall(TSGetSNES(ts, &snes));
831 PetscCall(SNESLoad(snes, viewer));
832 /* function and Jacobian context for SNES when used with TS is always ts object */
833 PetscCall(SNESSetFunction(snes, NULL, NULL, ts));
834 PetscCall(SNESSetJacobian(snes, NULL, NULL, NULL, ts));
835 PetscFunctionReturn(PETSC_SUCCESS);
836 }
837
838 /*@
839 TSGLEESetType - Set the type of `TSGLEE` scheme
840
841 Logically Collective
842
843 Input Parameters:
844 + ts - timestepping context
845 - gleetype - type of `TSGLEE` scheme
846
847 Level: intermediate
848
849 .seealso: [](ch_ts), `TSGLEEGetType()`, `TSGLEE`
850 @*/
TSGLEESetType(TS ts,TSGLEEType gleetype)851 PetscErrorCode TSGLEESetType(TS ts, TSGLEEType gleetype)
852 {
853 PetscFunctionBegin;
854 PetscValidHeaderSpecific(ts, TS_CLASSID, 1);
855 PetscAssertPointer(gleetype, 2);
856 PetscTryMethod(ts, "TSGLEESetType_C", (TS, TSGLEEType), (ts, gleetype));
857 PetscFunctionReturn(PETSC_SUCCESS);
858 }
859
860 /*@
861 TSGLEEGetType - Get the type of `TSGLEE` scheme
862
863 Logically Collective
864
865 Input Parameter:
866 . ts - timestepping context
867
868 Output Parameter:
869 . gleetype - type of `TSGLEE` scheme
870
871 Level: intermediate
872
873 .seealso: [](ch_ts), `TSGLEE`, `TSGLEESetType()`
874 @*/
TSGLEEGetType(TS ts,TSGLEEType * gleetype)875 PetscErrorCode TSGLEEGetType(TS ts, TSGLEEType *gleetype)
876 {
877 PetscFunctionBegin;
878 PetscValidHeaderSpecific(ts, TS_CLASSID, 1);
879 PetscUseMethod(ts, "TSGLEEGetType_C", (TS, TSGLEEType *), (ts, gleetype));
880 PetscFunctionReturn(PETSC_SUCCESS);
881 }
882
TSGLEEGetType_GLEE(TS ts,TSGLEEType * gleetype)883 static PetscErrorCode TSGLEEGetType_GLEE(TS ts, TSGLEEType *gleetype)
884 {
885 TS_GLEE *glee = (TS_GLEE *)ts->data;
886
887 PetscFunctionBegin;
888 if (!glee->tableau) PetscCall(TSGLEESetType(ts, TSGLEEDefaultType));
889 *gleetype = glee->tableau->name;
890 PetscFunctionReturn(PETSC_SUCCESS);
891 }
TSGLEESetType_GLEE(TS ts,TSGLEEType gleetype)892 static PetscErrorCode TSGLEESetType_GLEE(TS ts, TSGLEEType gleetype)
893 {
894 TS_GLEE *glee = (TS_GLEE *)ts->data;
895 PetscBool match;
896 GLEETableauLink link;
897
898 PetscFunctionBegin;
899 if (glee->tableau) {
900 PetscCall(PetscStrcmp(glee->tableau->name, gleetype, &match));
901 if (match) PetscFunctionReturn(PETSC_SUCCESS);
902 }
903 for (link = GLEETableauList; link; link = link->next) {
904 PetscCall(PetscStrcmp(link->tab.name, gleetype, &match));
905 if (match) {
906 PetscCall(TSReset_GLEE(ts));
907 glee->tableau = &link->tab;
908 PetscFunctionReturn(PETSC_SUCCESS);
909 }
910 }
911 SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_UNKNOWN_TYPE, "Could not find '%s'", gleetype);
912 }
913
TSGetStages_GLEE(TS ts,PetscInt * ns,Vec ** Y)914 static PetscErrorCode TSGetStages_GLEE(TS ts, PetscInt *ns, Vec **Y)
915 {
916 TS_GLEE *glee = (TS_GLEE *)ts->data;
917
918 PetscFunctionBegin;
919 if (ns) *ns = glee->tableau->s;
920 if (Y) *Y = glee->YStage;
921 PetscFunctionReturn(PETSC_SUCCESS);
922 }
923
TSGetSolutionComponents_GLEE(TS ts,PetscInt * n,Vec * Y)924 static PetscErrorCode TSGetSolutionComponents_GLEE(TS ts, PetscInt *n, Vec *Y)
925 {
926 TS_GLEE *glee = (TS_GLEE *)ts->data;
927 GLEETableau tab = glee->tableau;
928
929 PetscFunctionBegin;
930 if (!Y) *n = tab->r;
931 else {
932 PetscCheck(*n >= 0 && *n < tab->r, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_OUTOFRANGE, "Second argument (%" PetscInt_FMT ") out of range[0,%" PetscInt_FMT "].", *n, tab->r - 1);
933 PetscCall(VecCopy(glee->Y[*n], *Y));
934 }
935 PetscFunctionReturn(PETSC_SUCCESS);
936 }
937
TSGetAuxSolution_GLEE(TS ts,Vec * X)938 static PetscErrorCode TSGetAuxSolution_GLEE(TS ts, Vec *X)
939 {
940 TS_GLEE *glee = (TS_GLEE *)ts->data;
941 GLEETableau tab = glee->tableau;
942 PetscReal *F = tab->Fembed;
943 PetscInt r = tab->r;
944 Vec *Y = glee->Y;
945 PetscScalar *wr = glee->rwork;
946 PetscInt i;
947
948 PetscFunctionBegin;
949 PetscCall(VecZeroEntries(*X));
950 for (i = 0; i < r; i++) wr[i] = F[i];
951 PetscCall(VecMAXPY(*X, r, wr, Y));
952 PetscFunctionReturn(PETSC_SUCCESS);
953 }
954
TSGetTimeError_GLEE(TS ts,PetscInt n,Vec * X)955 static PetscErrorCode TSGetTimeError_GLEE(TS ts, PetscInt n, Vec *X)
956 {
957 TS_GLEE *glee = (TS_GLEE *)ts->data;
958 GLEETableau tab = glee->tableau;
959 PetscReal *F = tab->Ferror;
960 PetscInt r = tab->r;
961 Vec *Y = glee->Y;
962 PetscScalar *wr = glee->rwork;
963 PetscInt i;
964
965 PetscFunctionBegin;
966 PetscCall(VecZeroEntries(*X));
967 if (n == 0) {
968 for (i = 0; i < r; i++) wr[i] = F[i];
969 PetscCall(VecMAXPY(*X, r, wr, Y));
970 } else if (n == -1) {
971 *X = glee->yGErr;
972 }
973 PetscFunctionReturn(PETSC_SUCCESS);
974 }
975
TSSetTimeError_GLEE(TS ts,Vec X)976 static PetscErrorCode TSSetTimeError_GLEE(TS ts, Vec X)
977 {
978 TS_GLEE *glee = (TS_GLEE *)ts->data;
979 GLEETableau tab = glee->tableau;
980 PetscReal *S = tab->Serror;
981 PetscInt r = tab->r, i;
982 Vec *Y = glee->Y;
983
984 PetscFunctionBegin;
985 PetscCheck(r == 2, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "TSSetTimeError_GLEE not supported for '%s' with r=%" PetscInt_FMT ".", tab->name, tab->r);
986 for (i = 1; i < r; i++) {
987 PetscCall(VecCopy(ts->vec_sol, Y[i]));
988 PetscCall(VecAXPBY(Y[i], S[0], S[1], X));
989 PetscCall(VecCopy(X, glee->yGErr));
990 }
991 PetscFunctionReturn(PETSC_SUCCESS);
992 }
993
TSDestroy_GLEE(TS ts)994 static PetscErrorCode TSDestroy_GLEE(TS ts)
995 {
996 PetscFunctionBegin;
997 PetscCall(TSReset_GLEE(ts));
998 if (ts->dm) {
999 PetscCall(DMCoarsenHookRemove(ts->dm, DMCoarsenHook_TSGLEE, DMRestrictHook_TSGLEE, ts));
1000 PetscCall(DMSubDomainHookRemove(ts->dm, DMSubDomainHook_TSGLEE, DMSubDomainRestrictHook_TSGLEE, ts));
1001 }
1002 PetscCall(PetscFree(ts->data));
1003 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSGLEEGetType_C", NULL));
1004 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSGLEESetType_C", NULL));
1005 PetscFunctionReturn(PETSC_SUCCESS);
1006 }
1007
1008 /*MC
1009 TSGLEE - ODE and DAE solver using General Linear with Error Estimation schemes
1010
1011 The user should provide the right-hand side of the equation using `TSSetRHSFunction()`
1012 and for `TSGLEEi1` the Jacobian of the right-hand side using `TSSetRHSJacobian()`
1013
1014 Level: beginner
1015
1016 Notes:
1017 The default is `TSGLEE35`, it can be changed with `TSGLEESetType()` or `-ts_glee_type type`
1018
1019 The only implicit scheme is `TSGLEEi1`
1020
1021 .seealso: [](ch_ts), [](sec_ts_glee), `TSCreate()`, `TS`, `TSSetType()`, `TSGLEESetType()`, `TSGLEEGetType()`,
1022 `TSGLEE23`, `TSGLEE24`, `TSGLEE35`, `TSGLEE25I`, `TSGLEEEXRK2A`,
1023 `TSGLEERK32G1`, `TSGLEERK285EX`, `TSGLEEType`, `TSGLEERegister()`, `TSType`
1024 M*/
TSCreate_GLEE(TS ts)1025 PETSC_EXTERN PetscErrorCode TSCreate_GLEE(TS ts)
1026 {
1027 TS_GLEE *th;
1028
1029 PetscFunctionBegin;
1030 PetscCall(TSGLEEInitializePackage());
1031
1032 ts->ops->reset = TSReset_GLEE;
1033 ts->ops->destroy = TSDestroy_GLEE;
1034 ts->ops->view = TSView_GLEE;
1035 ts->ops->load = TSLoad_GLEE;
1036 ts->ops->setup = TSSetUp_GLEE;
1037 ts->ops->step = TSStep_GLEE;
1038 ts->ops->interpolate = TSInterpolate_GLEE;
1039 ts->ops->evaluatestep = TSEvaluateStep_GLEE;
1040 ts->ops->setfromoptions = TSSetFromOptions_GLEE;
1041 ts->ops->getstages = TSGetStages_GLEE;
1042 ts->ops->snesfunction = SNESTSFormFunction_GLEE;
1043 ts->ops->snesjacobian = SNESTSFormJacobian_GLEE;
1044 ts->ops->getsolutioncomponents = TSGetSolutionComponents_GLEE;
1045 ts->ops->getauxsolution = TSGetAuxSolution_GLEE;
1046 ts->ops->gettimeerror = TSGetTimeError_GLEE;
1047 ts->ops->settimeerror = TSSetTimeError_GLEE;
1048 ts->ops->startingmethod = TSStartingMethod_GLEE;
1049 ts->default_adapt_type = TSADAPTGLEE;
1050
1051 ts->usessnes = PETSC_TRUE;
1052
1053 PetscCall(PetscNew(&th));
1054 ts->data = (void *)th;
1055
1056 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSGLEEGetType_C", TSGLEEGetType_GLEE));
1057 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSGLEESetType_C", TSGLEESetType_GLEE));
1058 PetscFunctionReturn(PETSC_SUCCESS);
1059 }
1060