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