Actual source code: glee.c
1: /*
2: Code for time stepping with the General Linear with Error Estimation method
4: Notes:
5: The general system is written as
7: Udot = F(t,U)
9: */
10: #include <petsc/private/tsimpl.h>
11: #include <petscdm.h>
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";
21: static TSGLEEType TSGLEEDefaultType = TSGLEE35;
22: static PetscBool TSGLEERegisterAllCalled;
23: static PetscBool TSGLEEPackageInitialized;
24: static PetscInt explicit_stage_time_id;
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;
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;
64: /*MC
65: TSGLEE23 - Second order three stage GLEE method
67: This method has three stages.
68: s = 3, r = 2
70: Level: advanced
72: .seealso: `TSGLEE`
73: M*/
74: /*MC
75: TSGLEE24 - Second order four stage GLEE method
77: This method has four stages.
78: s = 4, r = 2
80: Level: advanced
82: .seealso: `TSGLEE`
83: M*/
84: /*MC
85: TSGLEE25i - Second order five stage GLEE method
87: This method has five stages.
88: s = 5, r = 2
90: Level: advanced
92: .seealso: `TSGLEE`
93: M*/
94: /*MC
95: TSGLEE35 - Third order five stage GLEE method
97: This method has five stages.
98: s = 5, r = 2
100: Level: advanced
102: .seealso: `TSGLEE`
103: M*/
104: /*MC
105: TSGLEEEXRK2A - Second order six stage GLEE method
107: This method has six stages.
108: s = 6, r = 2
110: Level: advanced
112: .seealso: `TSGLEE`
113: M*/
114: /*MC
115: TSGLEERK32G1 - Third order eight stage GLEE method
117: This method has eight stages.
118: s = 8, r = 2
120: Level: advanced
122: .seealso: `TSGLEE`
123: M*/
124: /*MC
125: TSGLEERK285EX - Second order nine stage GLEE method
127: This method has nine stages.
128: s = 9, r = 2
130: Level: advanced
132: .seealso: `TSGLEE`
133: M*/
135: /*@C
136: TSGLEERegisterAll - Registers all of the General Linear with Error Estimation methods in TSGLEE
138: Not Collective, but should be called by all processes which will need the schemes to be registered
140: Level: advanced
142: .seealso: `TSGLEERegisterDestroy()`
143: @*/
144: PetscErrorCode TSGLEERegisterAll(void)
145: {
146: if (TSGLEERegisterAllCalled) return 0;
147: TSGLEERegisterAllCalled = PETSC_TRUE;
149: {
150: #define GAMMA 0.5
151: /* y-eps form */
152: const PetscInt p = 1, s = 3, r = 2;
153: const PetscReal A[3][3] =
154: {
155: {1.0, 0, 0 },
156: {0, 0.5, 0 },
157: {0, 0.5, 0.5}
158: },
159: 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};
160: 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);
161: }
162: {
163: #undef GAMMA
164: #define GAMMA 0.0
165: /* y-eps form */
166: const PetscInt p = 2, s = 3, r = 2;
167: const PetscReal A[3][3] =
168: {
169: {0, 0, 0},
170: {1, 0, 0},
171: {0.25, 0.25, 0}
172: },
173: 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};
174: 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);
175: }
176: {
177: #undef GAMMA
178: #define GAMMA 0.0
179: /* y-y~ form */
180: const PetscInt p = 2, s = 4, r = 2;
181: const PetscReal A[4][4] =
182: {
183: {0, 0, 0, 0},
184: {0.75, 0, 0, 0},
185: {0.25, 29.0 / 60.0, 0, 0},
186: {-21.0 / 44.0, 145.0 / 44.0, -20.0 / 11.0, 0}
187: },
188: 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};
189: 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);
190: }
191: {
192: #undef GAMMA
193: #define GAMMA 0.0
194: /* y-y~ form */
195: const PetscInt p = 2, s = 5, r = 2;
196: const PetscReal A[5][5] =
197: {
198: {0, 0, 0, 0, 0},
199: {-0.94079244066783383269, 0, 0, 0, 0},
200: {0.64228187778301907108, 0.10915356933958500042, 0, 0, 0},
201: {-0.51764297742287450812, 0.74414270351096040738, -0.71404164927824538121, 0, 0},
202: {-0.44696561556825969206, -0.76768425657590196518, 0.20111608138142987881, 0.93828186737840469796, 0}
203: },
204: 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};
205: 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);
206: }
207: {
208: #undef GAMMA
209: #define GAMMA 0.0
210: /* y-y~ form */
211: const PetscInt p = 3, s = 5, r = 2;
212: const PetscReal A[5][5] =
213: {
214: {0, 0, 0, 0, 0},
215: {-2169604947363702313.0 / 24313474998937147335.0, 0, 0, 0, 0},
216: {46526746497697123895.0 / 94116917485856474137.0, -10297879244026594958.0 / 49199457603717988219.0, 0, 0, 0},
217: {23364788935845982499.0 / 87425311444725389446.0, -79205144337496116638.0 / 148994349441340815519.0, 40051189859317443782.0 / 36487615018004984309.0, 0, 0},
218: {42089522664062539205.0 / 124911313006412840286.0, -15074384760342762939.0 / 137927286865289746282.0, -62274678522253371016.0 / 125918573676298591413.0, 13755475729852471739.0 / 79257927066651693390.0, 0}
219: },
220: 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};
221: 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);
222: }
223: {
224: #undef GAMMA
225: #define GAMMA 0.25
226: /* y-eps form */
227: const PetscInt p = 2, s = 6, r = 2;
228: const PetscReal A[6][6] =
229: {
230: {0, 0, 0, 0, 0, 0},
231: {1, 0, 0, 0, 0, 0},
232: {0, 0, 0, 0, 0, 0},
233: {0, 0, 0.5, 0, 0, 0},
234: {0, 0, 0.25, 0.25, 0, 0},
235: {0, 0, 0.25, 0.25, 0.5, 0}
236: },
237: 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};
238: 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);
239: }
240: {
241: #undef GAMMA
242: #define GAMMA 0.0
243: /* y-eps form */
244: const PetscInt p = 3, s = 8, r = 2;
245: const PetscReal A[8][8] =
246: {
247: {0, 0, 0, 0, 0, 0, 0, 0},
248: {0.5, 0, 0, 0, 0, 0, 0, 0},
249: {-1, 2, 0, 0, 0, 0, 0, 0},
250: {1.0 / 6.0, 2.0 / 3.0, 1.0 / 6.0, 0, 0, 0, 0, 0},
251: {0, 0, 0, 0, 0, 0, 0, 0},
252: {-7.0 / 24.0, 1.0 / 3.0, 1.0 / 12.0, -1.0 / 8.0, 0.5, 0, 0, 0},
253: {7.0 / 6.0, -4.0 / 3.0, -1.0 / 3.0, 0.5, -1.0, 2.0, 0, 0},
254: {0, 0, 0, 0, 1.0 / 6.0, 2.0 / 3.0, 1.0 / 6.0, 0}
255: },
256: 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};
257: 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);
258: }
259: {
260: #undef GAMMA
261: #define GAMMA 0.25
262: /* y-eps form */
263: const PetscInt p = 2, s = 9, r = 2;
264: const PetscReal A[9][9] =
265: {
266: {0, 0, 0, 0, 0, 0, 0, 0, 0},
267: {0.585786437626904966, 0, 0, 0, 0, 0, 0, 0, 0},
268: {0.149999999999999994, 0.849999999999999978, 0, 0, 0, 0, 0, 0, 0},
269: {0, 0, 0, 0, 0, 0, 0, 0, 0},
270: {0, 0, 0, 0.292893218813452483, 0, 0, 0, 0, 0},
271: {0, 0, 0, 0.0749999999999999972, 0.424999999999999989, 0, 0, 0, 0},
272: {0, 0, 0, 0.176776695296636893, 0.176776695296636893, 0.146446609406726241, 0, 0, 0},
273: {0, 0, 0, 0.176776695296636893, 0.176776695296636893, 0.146446609406726241, 0.292893218813452483, 0, 0},
274: {0, 0, 0, 0.176776695296636893, 0.176776695296636893, 0.146446609406726241, 0.0749999999999999972, 0.424999999999999989, 0}
275: },
276: 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};
277: 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);
278: }
280: return 0;
281: }
283: /*@C
284: TSGLEERegisterDestroy - Frees the list of schemes that were registered by TSGLEERegister().
286: Not Collective
288: Level: advanced
290: .seealso: `TSGLEERegister()`, `TSGLEERegisterAll()`
291: @*/
292: PetscErrorCode TSGLEERegisterDestroy(void)
293: {
294: GLEETableauLink link;
296: while ((link = GLEETableauList)) {
297: GLEETableau t = &link->tab;
298: GLEETableauList = link->next;
299: PetscFree5(t->A, t->B, t->U, t->V, t->c);
300: PetscFree2(t->S, t->F);
301: PetscFree(t->Fembed);
302: PetscFree(t->Ferror);
303: PetscFree(t->Serror);
304: PetscFree(t->binterp);
305: PetscFree(t->name);
306: PetscFree(link);
307: }
308: TSGLEERegisterAllCalled = PETSC_FALSE;
309: return 0;
310: }
312: /*@C
313: TSGLEEInitializePackage - This function initializes everything in the TSGLEE package. It is called
314: from TSInitializePackage().
316: Level: developer
318: .seealso: `PetscInitialize()`
319: @*/
320: PetscErrorCode TSGLEEInitializePackage(void)
321: {
322: if (TSGLEEPackageInitialized) return 0;
323: TSGLEEPackageInitialized = PETSC_TRUE;
324: TSGLEERegisterAll();
325: PetscObjectComposedDataRegister(&explicit_stage_time_id);
326: PetscRegisterFinalize(TSGLEEFinalizePackage);
327: return 0;
328: }
330: /*@C
331: TSGLEEFinalizePackage - This function destroys everything in the TSGLEE package. It is
332: called from PetscFinalize().
334: Level: developer
336: .seealso: `PetscFinalize()`
337: @*/
338: PetscErrorCode TSGLEEFinalizePackage(void)
339: {
340: TSGLEEPackageInitialized = PETSC_FALSE;
341: TSGLEERegisterDestroy();
342: return 0;
343: }
345: /*@C
346: TSGLEERegister - register an GLEE scheme by providing the entries in the Butcher tableau
348: Not Collective, but the same schemes should be registered on all processes on which they will be used
350: Input Parameters:
351: + name - identifier for method
352: . order - order of method
353: . s - number of stages
354: . r - number of steps
355: . gamma - LTE ratio
356: . A - stage coefficients (dimension s*s, row-major)
357: . B - step completion coefficients (dimension r*s, row-major)
358: . U - method coefficients (dimension s*r, row-major)
359: . V - method coefficients (dimension r*r, row-major)
360: . S - starting coefficients
361: . F - finishing coefficients
362: . c - abscissa (dimension s; NULL to use row sums of A)
363: . Fembed - step completion coefficients for embedded method
364: . Ferror - error computation coefficients
365: . Serror - error initialization coefficients
366: . pinterp - order of interpolation (0 if unavailable)
367: - binterp - array of interpolation coefficients (NULL if unavailable)
369: Notes:
370: Several GLEE methods are provided, this function is only needed to create new methods.
372: Level: advanced
374: .seealso: `TSGLEE`
375: @*/
376: 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[])
377: {
378: GLEETableauLink link;
379: GLEETableau t;
380: PetscInt i, j;
382: TSGLEEInitializePackage();
383: PetscNew(&link);
384: t = &link->tab;
385: PetscStrallocpy(name, &t->name);
386: t->order = order;
387: t->s = s;
388: t->r = r;
389: t->gamma = gamma;
390: PetscMalloc5(s * s, &t->A, r * r, &t->V, s, &t->c, r * s, &t->B, s * r, &t->U);
391: PetscMalloc2(r, &t->S, r, &t->F);
392: PetscArraycpy(t->A, A, s * s);
393: PetscArraycpy(t->B, B, r * s);
394: PetscArraycpy(t->U, U, s * r);
395: PetscArraycpy(t->V, V, r * r);
396: PetscArraycpy(t->S, S, r);
397: PetscArraycpy(t->F, F, r);
398: if (c) {
399: PetscArraycpy(t->c, c, s);
400: } else {
401: for (i = 0; i < s; i++)
402: for (j = 0, t->c[i] = 0; j < s; j++) t->c[i] += A[i * s + j];
403: }
404: PetscMalloc1(r, &t->Fembed);
405: PetscMalloc1(r, &t->Ferror);
406: PetscMalloc1(r, &t->Serror);
407: PetscArraycpy(t->Fembed, Fembed, r);
408: PetscArraycpy(t->Ferror, Ferror, r);
409: PetscArraycpy(t->Serror, Serror, r);
410: t->pinterp = pinterp;
411: PetscMalloc1(s * pinterp, &t->binterp);
412: PetscArraycpy(t->binterp, binterp, s * pinterp);
414: link->next = GLEETableauList;
415: GLEETableauList = link;
416: return 0;
417: }
419: static PetscErrorCode TSEvaluateStep_GLEE(TS ts, PetscInt order, Vec X, PetscBool *done)
420: {
421: TS_GLEE *glee = (TS_GLEE *)ts->data;
422: GLEETableau tab = glee->tableau;
423: PetscReal h, *B = tab->B, *V = tab->V, *F = tab->F, *Fembed = tab->Fembed;
424: PetscInt s = tab->s, r = tab->r, i, j;
425: Vec *Y = glee->Y, *YdotStage = glee->YdotStage;
426: PetscScalar *ws = glee->swork, *wr = glee->rwork;
429: switch (glee->status) {
430: case TS_STEP_INCOMPLETE:
431: case TS_STEP_PENDING:
432: h = ts->time_step;
433: break;
434: case TS_STEP_COMPLETE:
435: h = ts->ptime - ts->ptime_prev;
436: break;
437: default:
438: SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_PLIB, "Invalid TSStepStatus");
439: }
441: if (order == tab->order) {
442: /* Note: Irrespective of whether status is TS_STEP_INCOMPLETE
443: or TS_STEP_COMPLETE, glee->X has the solution at the
444: beginning of the time step. So no need to roll-back.
445: */
446: if (glee->status == TS_STEP_INCOMPLETE) {
447: for (i = 0; i < r; i++) {
448: VecZeroEntries(Y[i]);
449: for (j = 0; j < r; j++) wr[j] = V[i * r + j];
450: VecMAXPY(Y[i], r, wr, glee->X);
451: for (j = 0; j < s; j++) ws[j] = h * B[i * s + j];
452: VecMAXPY(Y[i], s, ws, YdotStage);
453: }
454: VecZeroEntries(X);
455: for (j = 0; j < r; j++) wr[j] = F[j];
456: VecMAXPY(X, r, wr, Y);
457: } else VecCopy(ts->vec_sol, X);
458: return 0;
460: } else if (order == tab->order - 1) {
461: /* Complete with the embedded method (Fembed) */
462: for (i = 0; i < r; i++) {
463: VecZeroEntries(Y[i]);
464: for (j = 0; j < r; j++) wr[j] = V[i * r + j];
465: VecMAXPY(Y[i], r, wr, glee->X);
466: for (j = 0; j < s; j++) ws[j] = h * B[i * s + j];
467: VecMAXPY(Y[i], s, ws, YdotStage);
468: }
469: VecZeroEntries(X);
470: for (j = 0; j < r; j++) wr[j] = Fembed[j];
471: VecMAXPY(X, r, wr, Y);
473: if (done) *done = PETSC_TRUE;
474: return 0;
475: }
476: if (done) *done = PETSC_FALSE;
477: else SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "GLEE '%s' of order %" PetscInt_FMT " cannot evaluate step at order %" PetscInt_FMT, tab->name, tab->order, order);
478: return 0;
479: }
481: static PetscErrorCode TSStep_GLEE(TS ts)
482: {
483: TS_GLEE *glee = (TS_GLEE *)ts->data;
484: GLEETableau tab = glee->tableau;
485: const PetscInt s = tab->s, r = tab->r;
486: PetscReal *A = tab->A, *U = tab->U, *F = tab->F, *c = tab->c;
487: Vec *Y = glee->Y, *X = glee->X, *YStage = glee->YStage, *YdotStage = glee->YdotStage, W = glee->W;
488: SNES snes;
489: PetscScalar *ws = glee->swork, *wr = glee->rwork;
490: TSAdapt adapt;
491: PetscInt i, j, reject, next_scheme, its, lits;
492: PetscReal next_time_step;
493: PetscReal t;
494: PetscBool accept;
496: PetscCitationsRegister(citation, &cited);
498: for (i = 0; i < r; i++) VecCopy(Y[i], X[i]);
500: TSGetSNES(ts, &snes);
501: next_time_step = ts->time_step;
502: t = ts->ptime;
503: accept = PETSC_TRUE;
504: glee->status = TS_STEP_INCOMPLETE;
506: for (reject = 0; reject < ts->max_reject && !ts->reason; reject++, ts->reject++) {
507: PetscReal h = ts->time_step;
508: TSPreStep(ts);
510: for (i = 0; i < s; i++) {
511: glee->stage_time = t + h * c[i];
512: TSPreStage(ts, glee->stage_time);
514: if (A[i * s + i] == 0) { /* Explicit stage */
515: VecZeroEntries(YStage[i]);
516: for (j = 0; j < r; j++) wr[j] = U[i * r + j];
517: VecMAXPY(YStage[i], r, wr, X);
518: for (j = 0; j < i; j++) ws[j] = h * A[i * s + j];
519: VecMAXPY(YStage[i], i, ws, YdotStage);
520: } else { /* Implicit stage */
521: glee->scoeff = 1.0 / A[i * s + i];
522: /* compute right-hand-side */
523: VecZeroEntries(W);
524: for (j = 0; j < r; j++) wr[j] = U[i * r + j];
525: VecMAXPY(W, r, wr, X);
526: for (j = 0; j < i; j++) ws[j] = h * A[i * s + j];
527: VecMAXPY(W, i, ws, YdotStage);
528: VecScale(W, glee->scoeff / h);
529: /* set initial guess */
530: VecCopy(i > 0 ? YStage[i - 1] : ts->vec_sol, YStage[i]);
531: /* solve for this stage */
532: SNESSolve(snes, W, YStage[i]);
533: SNESGetIterationNumber(snes, &its);
534: SNESGetLinearSolveIterations(snes, &lits);
535: ts->snes_its += its;
536: ts->ksp_its += lits;
537: }
538: TSGetAdapt(ts, &adapt);
539: TSAdaptCheckStage(adapt, ts, glee->stage_time, YStage[i], &accept);
540: if (!accept) goto reject_step;
541: TSPostStage(ts, glee->stage_time, i, YStage);
542: TSComputeRHSFunction(ts, t + h * c[i], YStage[i], YdotStage[i]);
543: }
544: TSEvaluateStep(ts, tab->order, ts->vec_sol, NULL);
545: glee->status = TS_STEP_PENDING;
547: /* Register only the current method as a candidate because we're not supporting multiple candidates yet. */
548: TSGetAdapt(ts, &adapt);
549: TSAdaptCandidatesClear(adapt);
550: TSAdaptCandidateAdd(adapt, tab->name, tab->order, 1, tab->ccfl, (PetscReal)tab->s, PETSC_TRUE);
551: TSAdaptChoose(adapt, ts, ts->time_step, &next_scheme, &next_time_step, &accept);
552: if (accept) {
553: /* ignore next_scheme for now */
554: ts->ptime += ts->time_step;
555: ts->time_step = next_time_step;
556: glee->status = TS_STEP_COMPLETE;
557: /* compute and store the global error */
558: /* Note: this is not needed if TSAdaptGLEE is not used */
559: TSGetTimeError(ts, 0, &(glee->yGErr));
560: PetscObjectComposedDataSetReal((PetscObject)ts->vec_sol, explicit_stage_time_id, ts->ptime);
561: break;
562: } else { /* Roll back the current step */
563: for (j = 0; j < r; j++) wr[j] = F[j];
564: VecMAXPY(ts->vec_sol, r, wr, X);
565: ts->time_step = next_time_step;
566: glee->status = TS_STEP_INCOMPLETE;
567: }
568: reject_step:
569: continue;
570: }
571: if (glee->status != TS_STEP_COMPLETE && !ts->reason) ts->reason = TS_DIVERGED_STEP_REJECTED;
572: return 0;
573: }
575: static PetscErrorCode TSInterpolate_GLEE(TS ts, PetscReal itime, Vec X)
576: {
577: TS_GLEE *glee = (TS_GLEE *)ts->data;
578: PetscInt s = glee->tableau->s, pinterp = glee->tableau->pinterp, i, j;
579: PetscReal h, tt, t;
580: PetscScalar *b;
581: const PetscReal *B = glee->tableau->binterp;
584: switch (glee->status) {
585: case TS_STEP_INCOMPLETE:
586: case TS_STEP_PENDING:
587: h = ts->time_step;
588: t = (itime - ts->ptime) / h;
589: break;
590: case TS_STEP_COMPLETE:
591: h = ts->ptime - ts->ptime_prev;
592: t = (itime - ts->ptime) / h + 1; /* In the interval [0,1] */
593: break;
594: default:
595: SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_PLIB, "Invalid TSStepStatus");
596: }
597: PetscMalloc1(s, &b);
598: for (i = 0; i < s; i++) b[i] = 0;
599: for (j = 0, tt = t; j < pinterp; j++, tt *= t) {
600: for (i = 0; i < s; i++) b[i] += h * B[i * pinterp + j] * tt;
601: }
602: VecCopy(glee->YStage[0], X);
603: VecMAXPY(X, s, b, glee->YdotStage);
604: PetscFree(b);
605: return 0;
606: }
608: /*------------------------------------------------------------*/
609: static PetscErrorCode TSReset_GLEE(TS ts)
610: {
611: TS_GLEE *glee = (TS_GLEE *)ts->data;
612: PetscInt s, r;
614: if (!glee->tableau) return 0;
615: s = glee->tableau->s;
616: r = glee->tableau->r;
617: VecDestroyVecs(r, &glee->Y);
618: VecDestroyVecs(r, &glee->X);
619: VecDestroyVecs(s, &glee->YStage);
620: VecDestroyVecs(s, &glee->YdotStage);
621: VecDestroy(&glee->Ydot);
622: VecDestroy(&glee->yGErr);
623: VecDestroy(&glee->W);
624: PetscFree2(glee->swork, glee->rwork);
625: return 0;
626: }
628: static PetscErrorCode TSGLEEGetVecs(TS ts, DM dm, Vec *Ydot)
629: {
630: TS_GLEE *glee = (TS_GLEE *)ts->data;
632: if (Ydot) {
633: if (dm && dm != ts->dm) {
634: DMGetNamedGlobalVector(dm, "TSGLEE_Ydot", Ydot);
635: } else *Ydot = glee->Ydot;
636: }
637: return 0;
638: }
640: static PetscErrorCode TSGLEERestoreVecs(TS ts, DM dm, Vec *Ydot)
641: {
642: if (Ydot) {
643: if (dm && dm != ts->dm) DMRestoreNamedGlobalVector(dm, "TSGLEE_Ydot", Ydot);
644: }
645: return 0;
646: }
648: /*
649: This defines the nonlinear equation that is to be solved with SNES
650: */
651: static PetscErrorCode SNESTSFormFunction_GLEE(SNES snes, Vec X, Vec F, TS ts)
652: {
653: TS_GLEE *glee = (TS_GLEE *)ts->data;
654: DM dm, dmsave;
655: Vec Ydot;
656: PetscReal shift = glee->scoeff / ts->time_step;
658: SNESGetDM(snes, &dm);
659: TSGLEEGetVecs(ts, dm, &Ydot);
660: /* Set Ydot = shift*X */
661: VecCopy(X, Ydot);
662: VecScale(Ydot, shift);
663: dmsave = ts->dm;
664: ts->dm = dm;
666: TSComputeIFunction(ts, glee->stage_time, X, Ydot, F, PETSC_FALSE);
668: ts->dm = dmsave;
669: TSGLEERestoreVecs(ts, dm, &Ydot);
670: return 0;
671: }
673: static PetscErrorCode SNESTSFormJacobian_GLEE(SNES snes, Vec X, Mat A, Mat B, TS ts)
674: {
675: TS_GLEE *glee = (TS_GLEE *)ts->data;
676: DM dm, dmsave;
677: Vec Ydot;
678: PetscReal shift = glee->scoeff / ts->time_step;
680: SNESGetDM(snes, &dm);
681: TSGLEEGetVecs(ts, dm, &Ydot);
682: /* glee->Ydot has already been computed in SNESTSFormFunction_GLEE (SNES guarantees this) */
683: dmsave = ts->dm;
684: ts->dm = dm;
686: TSComputeIJacobian(ts, glee->stage_time, X, Ydot, shift, A, B, PETSC_FALSE);
688: ts->dm = dmsave;
689: TSGLEERestoreVecs(ts, dm, &Ydot);
690: return 0;
691: }
693: static PetscErrorCode DMCoarsenHook_TSGLEE(DM fine, DM coarse, void *ctx)
694: {
695: return 0;
696: }
698: static PetscErrorCode DMRestrictHook_TSGLEE(DM fine, Mat restrct, Vec rscale, Mat inject, DM coarse, void *ctx)
699: {
700: return 0;
701: }
703: static PetscErrorCode DMSubDomainHook_TSGLEE(DM dm, DM subdm, void *ctx)
704: {
705: return 0;
706: }
708: static PetscErrorCode DMSubDomainRestrictHook_TSGLEE(DM dm, VecScatter gscat, VecScatter lscat, DM subdm, void *ctx)
709: {
710: return 0;
711: }
713: static PetscErrorCode TSSetUp_GLEE(TS ts)
714: {
715: TS_GLEE *glee = (TS_GLEE *)ts->data;
716: GLEETableau tab;
717: PetscInt s, r;
718: DM dm;
720: if (!glee->tableau) TSGLEESetType(ts, TSGLEEDefaultType);
721: tab = glee->tableau;
722: s = tab->s;
723: r = tab->r;
724: VecDuplicateVecs(ts->vec_sol, r, &glee->Y);
725: VecDuplicateVecs(ts->vec_sol, r, &glee->X);
726: VecDuplicateVecs(ts->vec_sol, s, &glee->YStage);
727: VecDuplicateVecs(ts->vec_sol, s, &glee->YdotStage);
728: VecDuplicate(ts->vec_sol, &glee->Ydot);
729: VecDuplicate(ts->vec_sol, &glee->yGErr);
730: VecZeroEntries(glee->yGErr);
731: VecDuplicate(ts->vec_sol, &glee->W);
732: PetscMalloc2(s, &glee->swork, r, &glee->rwork);
733: TSGetDM(ts, &dm);
734: DMCoarsenHookAdd(dm, DMCoarsenHook_TSGLEE, DMRestrictHook_TSGLEE, ts);
735: DMSubDomainHookAdd(dm, DMSubDomainHook_TSGLEE, DMSubDomainRestrictHook_TSGLEE, ts);
736: return 0;
737: }
739: PetscErrorCode TSStartingMethod_GLEE(TS ts)
740: {
741: TS_GLEE *glee = (TS_GLEE *)ts->data;
742: GLEETableau tab = glee->tableau;
743: PetscInt r = tab->r, i;
744: PetscReal *S = tab->S;
746: for (i = 0; i < r; i++) {
747: VecZeroEntries(glee->Y[i]);
748: VecAXPY(glee->Y[i], S[i], ts->vec_sol);
749: }
751: return 0;
752: }
754: /*------------------------------------------------------------*/
756: static PetscErrorCode TSSetFromOptions_GLEE(TS ts, PetscOptionItems *PetscOptionsObject)
757: {
758: char gleetype[256];
760: PetscOptionsHeadBegin(PetscOptionsObject, "GLEE ODE solver options");
761: {
762: GLEETableauLink link;
763: PetscInt count, choice;
764: PetscBool flg;
765: const char **namelist;
767: PetscStrncpy(gleetype, TSGLEEDefaultType, sizeof(gleetype));
768: for (link = GLEETableauList, count = 0; link; link = link->next, count++)
769: ;
770: PetscMalloc1(count, (char ***)&namelist);
771: for (link = GLEETableauList, count = 0; link; link = link->next, count++) namelist[count] = link->tab.name;
772: PetscOptionsEList("-ts_glee_type", "Family of GLEE method", "TSGLEESetType", (const char *const *)namelist, count, gleetype, &choice, &flg);
773: TSGLEESetType(ts, flg ? namelist[choice] : gleetype);
774: PetscFree(namelist);
775: }
776: PetscOptionsHeadEnd();
777: return 0;
778: }
780: static PetscErrorCode TSView_GLEE(TS ts, PetscViewer viewer)
781: {
782: TS_GLEE *glee = (TS_GLEE *)ts->data;
783: GLEETableau tab = glee->tableau;
784: PetscBool iascii;
786: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
787: if (iascii) {
788: TSGLEEType gleetype;
789: char buf[512];
790: TSGLEEGetType(ts, &gleetype);
791: PetscViewerASCIIPrintf(viewer, " GLEE type %s\n", gleetype);
792: PetscFormatRealArray(buf, sizeof(buf), "% 8.6f", tab->s, tab->c);
793: PetscViewerASCIIPrintf(viewer, " Abscissa c = %s\n", buf);
794: /* Note: print out r as well */
795: }
796: return 0;
797: }
799: static PetscErrorCode TSLoad_GLEE(TS ts, PetscViewer viewer)
800: {
801: SNES snes;
802: TSAdapt tsadapt;
804: TSGetAdapt(ts, &tsadapt);
805: TSAdaptLoad(tsadapt, viewer);
806: TSGetSNES(ts, &snes);
807: SNESLoad(snes, viewer);
808: /* function and Jacobian context for SNES when used with TS is always ts object */
809: SNESSetFunction(snes, NULL, NULL, ts);
810: SNESSetJacobian(snes, NULL, NULL, NULL, ts);
811: return 0;
812: }
814: /*@C
815: TSGLEESetType - Set the type of GLEE scheme
817: Logically collective
819: Input Parameters:
820: + ts - timestepping context
821: - gleetype - type of GLEE-scheme
823: Level: intermediate
825: .seealso: `TSGLEEGetType()`, `TSGLEE`
826: @*/
827: PetscErrorCode TSGLEESetType(TS ts, TSGLEEType gleetype)
828: {
831: PetscTryMethod(ts, "TSGLEESetType_C", (TS, TSGLEEType), (ts, gleetype));
832: return 0;
833: }
835: /*@C
836: TSGLEEGetType - Get the type of GLEE scheme
838: Logically collective
840: Input Parameter:
841: . ts - timestepping context
843: Output Parameter:
844: . gleetype - type of GLEE-scheme
846: Level: intermediate
848: .seealso: `TSGLEESetType()`
849: @*/
850: PetscErrorCode TSGLEEGetType(TS ts, TSGLEEType *gleetype)
851: {
853: PetscUseMethod(ts, "TSGLEEGetType_C", (TS, TSGLEEType *), (ts, gleetype));
854: return 0;
855: }
857: PetscErrorCode TSGLEEGetType_GLEE(TS ts, TSGLEEType *gleetype)
858: {
859: TS_GLEE *glee = (TS_GLEE *)ts->data;
861: if (!glee->tableau) TSGLEESetType(ts, TSGLEEDefaultType);
862: *gleetype = glee->tableau->name;
863: return 0;
864: }
865: PetscErrorCode TSGLEESetType_GLEE(TS ts, TSGLEEType gleetype)
866: {
867: TS_GLEE *glee = (TS_GLEE *)ts->data;
868: PetscBool match;
869: GLEETableauLink link;
871: if (glee->tableau) {
872: PetscStrcmp(glee->tableau->name, gleetype, &match);
873: if (match) return 0;
874: }
875: for (link = GLEETableauList; link; link = link->next) {
876: PetscStrcmp(link->tab.name, gleetype, &match);
877: if (match) {
878: TSReset_GLEE(ts);
879: glee->tableau = &link->tab;
880: return 0;
881: }
882: }
883: SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_UNKNOWN_TYPE, "Could not find '%s'", gleetype);
884: }
886: static PetscErrorCode TSGetStages_GLEE(TS ts, PetscInt *ns, Vec **Y)
887: {
888: TS_GLEE *glee = (TS_GLEE *)ts->data;
890: if (ns) *ns = glee->tableau->s;
891: if (Y) *Y = glee->YStage;
892: return 0;
893: }
895: PetscErrorCode TSGetSolutionComponents_GLEE(TS ts, PetscInt *n, Vec *Y)
896: {
897: TS_GLEE *glee = (TS_GLEE *)ts->data;
898: GLEETableau tab = glee->tableau;
900: if (!Y) *n = tab->r;
901: else {
902: if ((*n >= 0) && (*n < tab->r)) {
903: VecCopy(glee->Y[*n], *Y);
904: } else SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_OUTOFRANGE, "Second argument (%" PetscInt_FMT ") out of range[0,%" PetscInt_FMT "].", *n, tab->r - 1);
905: }
906: return 0;
907: }
909: PetscErrorCode TSGetAuxSolution_GLEE(TS ts, Vec *X)
910: {
911: TS_GLEE *glee = (TS_GLEE *)ts->data;
912: GLEETableau tab = glee->tableau;
913: PetscReal *F = tab->Fembed;
914: PetscInt r = tab->r;
915: Vec *Y = glee->Y;
916: PetscScalar *wr = glee->rwork;
917: PetscInt i;
919: VecZeroEntries(*X);
920: for (i = 0; i < r; i++) wr[i] = F[i];
921: VecMAXPY((*X), r, wr, Y);
922: return 0;
923: }
925: PetscErrorCode TSGetTimeError_GLEE(TS ts, PetscInt n, Vec *X)
926: {
927: TS_GLEE *glee = (TS_GLEE *)ts->data;
928: GLEETableau tab = glee->tableau;
929: PetscReal *F = tab->Ferror;
930: PetscInt r = tab->r;
931: Vec *Y = glee->Y;
932: PetscScalar *wr = glee->rwork;
933: PetscInt i;
935: VecZeroEntries(*X);
936: if (n == 0) {
937: for (i = 0; i < r; i++) wr[i] = F[i];
938: VecMAXPY((*X), r, wr, Y);
939: } else if (n == -1) {
940: *X = glee->yGErr;
941: }
942: return 0;
943: }
945: PetscErrorCode TSSetTimeError_GLEE(TS ts, Vec X)
946: {
947: TS_GLEE *glee = (TS_GLEE *)ts->data;
948: GLEETableau tab = glee->tableau;
949: PetscReal *S = tab->Serror;
950: PetscInt r = tab->r, i;
951: Vec *Y = glee->Y;
954: for (i = 1; i < r; i++) {
955: VecCopy(ts->vec_sol, Y[i]);
956: VecAXPBY(Y[i], S[0], S[1], X);
957: VecCopy(X, glee->yGErr);
958: }
959: return 0;
960: }
962: static PetscErrorCode TSDestroy_GLEE(TS ts)
963: {
964: TSReset_GLEE(ts);
965: if (ts->dm) {
966: DMCoarsenHookRemove(ts->dm, DMCoarsenHook_TSGLEE, DMRestrictHook_TSGLEE, ts);
967: DMSubDomainHookRemove(ts->dm, DMSubDomainHook_TSGLEE, DMSubDomainRestrictHook_TSGLEE, ts);
968: }
969: PetscFree(ts->data);
970: PetscObjectComposeFunction((PetscObject)ts, "TSGLEEGetType_C", NULL);
971: PetscObjectComposeFunction((PetscObject)ts, "TSGLEESetType_C", NULL);
972: return 0;
973: }
975: /* ------------------------------------------------------------ */
976: /*MC
977: TSGLEE - ODE and DAE solver using General Linear with Error Estimation schemes
979: The user should provide the right hand side of the equation
980: using TSSetRHSFunction().
982: Notes:
983: The default is TSGLEE35, it can be changed with TSGLEESetType() or -ts_glee_type
985: Level: beginner
987: .seealso: `TSCreate()`, `TS`, `TSSetType()`, `TSGLEESetType()`, `TSGLEEGetType()`,
988: `TSGLEE23`, `TTSGLEE24`, `TSGLEE35`, `TSGLEE25I`, `TSGLEEEXRK2A`,
989: `TSGLEERK32G1`, `TSGLEERK285EX`, `TSGLEEType`, `TSGLEERegister()`
991: M*/
992: PETSC_EXTERN PetscErrorCode TSCreate_GLEE(TS ts)
993: {
994: TS_GLEE *th;
996: TSGLEEInitializePackage();
998: ts->ops->reset = TSReset_GLEE;
999: ts->ops->destroy = TSDestroy_GLEE;
1000: ts->ops->view = TSView_GLEE;
1001: ts->ops->load = TSLoad_GLEE;
1002: ts->ops->setup = TSSetUp_GLEE;
1003: ts->ops->step = TSStep_GLEE;
1004: ts->ops->interpolate = TSInterpolate_GLEE;
1005: ts->ops->evaluatestep = TSEvaluateStep_GLEE;
1006: ts->ops->setfromoptions = TSSetFromOptions_GLEE;
1007: ts->ops->getstages = TSGetStages_GLEE;
1008: ts->ops->snesfunction = SNESTSFormFunction_GLEE;
1009: ts->ops->snesjacobian = SNESTSFormJacobian_GLEE;
1010: ts->ops->getsolutioncomponents = TSGetSolutionComponents_GLEE;
1011: ts->ops->getauxsolution = TSGetAuxSolution_GLEE;
1012: ts->ops->gettimeerror = TSGetTimeError_GLEE;
1013: ts->ops->settimeerror = TSSetTimeError_GLEE;
1014: ts->ops->startingmethod = TSStartingMethod_GLEE;
1015: ts->default_adapt_type = TSADAPTGLEE;
1017: ts->usessnes = PETSC_TRUE;
1019: PetscNew(&th);
1020: ts->data = (void *)th;
1022: PetscObjectComposeFunction((PetscObject)ts, "TSGLEEGetType_C", TSGLEEGetType_GLEE);
1023: PetscObjectComposeFunction((PetscObject)ts, "TSGLEESetType_C", TSGLEESetType_GLEE);
1024: return 0;
1025: }