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: TSGLEEi1 - Second order three stage implicit GLEE method
67: This method has two stages.
68: s = 3, r = 2
70: Level: advanced
72: .seealso: [](ch_ts), `TSGLEE`
73: M*/
74: /*MC
75: TSGLEE23 - Second order three stage explicit GLEE method
77: This method has three stages.
78: s = 3, r = 2
80: Level: advanced
82: .seealso: [](ch_ts), `TSGLEE`
83: M*/
84: /*MC
85: TSGLEE24 - Second order four stage explicit GLEE method
87: This method has four stages.
88: s = 4, r = 2
90: Level: advanced
92: .seealso: [](ch_ts), `TSGLEE`
93: M*/
94: /*MC
95: TSGLEE25i - Second order five stage explict GLEE method
97: This method has five stages.
98: s = 5, r = 2
100: Level: advanced
102: .seealso: [](ch_ts), `TSGLEE`
103: M*/
104: /*MC
105: TSGLEE35 - Third order five stage explicit GLEE method
107: This method has five stages.
108: s = 5, r = 2
110: Level: advanced
112: .seealso: [](ch_ts), `TSGLEE`
113: M*/
114: /*MC
115: TSGLEEEXRK2A - Second order six stage explict GLEE method
117: This method has six stages.
118: s = 6, r = 2
120: Level: advanced
122: .seealso: [](ch_ts), `TSGLEE`
123: M*/
124: /*MC
125: TSGLEERK32G1 - Third order eight stage explicit GLEE method
127: This method has eight stages.
128: s = 8, r = 2
130: Level: advanced
132: .seealso: [](ch_ts), `TSGLEE`
133: M*/
134: /*MC
135: TSGLEERK285EX - Second order nine stage explicit GLEE method
137: This method has nine stages.
138: s = 9, r = 2
140: Level: advanced
142: .seealso: [](ch_ts), `TSGLEE`
143: M*/
145: /*@C
146: TSGLEERegisterAll - Registers all of the General Linear with Error Estimation methods in `TSGLEE`
148: Not Collective, but should be called by all processes which will need the schemes to be registered
150: Level: advanced
152: .seealso: [](ch_ts), `TSGLEERegisterDestroy()`
153: @*/
154: PetscErrorCode TSGLEERegisterAll(void)
155: {
156: PetscFunctionBegin;
157: if (TSGLEERegisterAllCalled) PetscFunctionReturn(PETSC_SUCCESS);
158: TSGLEERegisterAllCalled = PETSC_TRUE;
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: }
293: /*@C
294: TSGLEERegisterDestroy - Frees the list of schemes that were registered by `TSGLEERegister()`.
296: Not Collective
298: Level: advanced
300: .seealso: [](ch_ts), `TSGLEERegister()`, `TSGLEERegisterAll()`
301: @*/
302: PetscErrorCode TSGLEERegisterDestroy(void)
303: {
304: GLEETableauLink link;
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: }
323: /*@C
324: TSGLEEInitializePackage - This function initializes everything in the `TSGLEE` package. It is called
325: from `TSInitializePackage()`.
327: Level: developer
329: .seealso: [](ch_ts), `PetscInitialize()`
330: @*/
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: }
342: /*@C
343: TSGLEEFinalizePackage - This function destroys everything in the `TSGLEE` package. It is
344: called from `PetscFinalize()`.
346: Level: developer
348: .seealso: [](ch_ts), `PetscFinalize()`
349: @*/
350: PetscErrorCode TSGLEEFinalizePackage(void)
351: {
352: PetscFunctionBegin;
353: TSGLEEPackageInitialized = PETSC_FALSE;
354: PetscCall(TSGLEERegisterDestroy());
355: PetscFunctionReturn(PETSC_SUCCESS);
356: }
358: /*@C
359: TSGLEERegister - register a new `TSGLEE` scheme by providing the entries in the Butcher tableau
361: Not Collective, but the same schemes should be registered on all processes on which they will be used, No Fortran Support
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)
382: Level: advanced
384: Note:
385: Several `TSGLEE` methods are provided, this function is only needed to create new methods.
387: .seealso: [](ch_ts), `TSGLEE`
388: @*/
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;
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));
428: link->next = GLEETableauList;
429: GLEETableauList = link;
430: PetscFunctionReturn(PETSC_SUCCESS);
431: }
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;
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: }
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);
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));
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: }
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;
510: PetscFunctionBegin;
511: PetscCall(PetscCitationsRegister(citation, &cited));
513: for (i = 0; i < r; i++) PetscCall(VecCopy(Y[i], X[i]));
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;
521: for (reject = 0; reject < ts->max_reject && !ts->reason; reject++, ts->reject++) {
522: PetscReal h = ts->time_step;
523: PetscCall(TSPreStep(ts));
525: for (i = 0; i < s; i++) {
526: glee->stage_time = t + h * c[i];
527: PetscCall(TSPreStage(ts, glee->stage_time));
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;
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: }
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;
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: }
624: static PetscErrorCode TSReset_GLEE(TS ts)
625: {
626: TS_GLEE *glee = (TS_GLEE *)ts->data;
627: PetscInt s, r;
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: }
644: static PetscErrorCode TSGLEEGetVecs(TS ts, DM dm, Vec *Ydot)
645: {
646: TS_GLEE *glee = (TS_GLEE *)ts->data;
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: }
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: }
665: /*
666: This defines the nonlinear equation that is to be solved with SNES
667: */
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;
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;
684: PetscCall(TSComputeIFunction(ts, glee->stage_time, X, Ydot, F, PETSC_FALSE));
686: ts->dm = dmsave;
687: PetscCall(TSGLEERestoreVecs(ts, dm, &Ydot));
688: PetscFunctionReturn(PETSC_SUCCESS);
689: }
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;
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;
705: PetscCall(TSComputeIJacobian(ts, glee->stage_time, X, Ydot, shift, A, B, PETSC_FALSE));
707: ts->dm = dmsave;
708: PetscCall(TSGLEERestoreVecs(ts, dm, &Ydot));
709: PetscFunctionReturn(PETSC_SUCCESS);
710: }
712: static PetscErrorCode DMCoarsenHook_TSGLEE(DM fine, DM coarse, void *ctx)
713: {
714: PetscFunctionBegin;
715: PetscFunctionReturn(PETSC_SUCCESS);
716: }
718: static PetscErrorCode DMRestrictHook_TSGLEE(DM fine, Mat restrct, Vec rscale, Mat inject, DM coarse, void *ctx)
719: {
720: PetscFunctionBegin;
721: PetscFunctionReturn(PETSC_SUCCESS);
722: }
724: static PetscErrorCode DMSubDomainHook_TSGLEE(DM dm, DM subdm, void *ctx)
725: {
726: PetscFunctionBegin;
727: PetscFunctionReturn(PETSC_SUCCESS);
728: }
730: static PetscErrorCode DMSubDomainRestrictHook_TSGLEE(DM dm, VecScatter gscat, VecScatter lscat, DM subdm, void *ctx)
731: {
732: PetscFunctionBegin;
733: PetscFunctionReturn(PETSC_SUCCESS);
734: }
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;
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: }
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;
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: }
778: static PetscErrorCode TSSetFromOptions_GLEE(TS ts, PetscOptionItems PetscOptionsObject)
779: {
780: char gleetype[256];
782: PetscFunctionBegin;
783: PetscOptionsHeadBegin(PetscOptionsObject, "GLEE ODE solver options");
784: {
785: GLEETableauLink link;
786: PetscInt count, choice;
787: PetscBool flg;
788: const char **namelist;
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: }
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;
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: }
822: static PetscErrorCode TSLoad_GLEE(TS ts, PetscViewer viewer)
823: {
824: SNES snes;
825: TSAdapt tsadapt;
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: }
838: /*@
839: TSGLEESetType - Set the type of `TSGLEE` scheme
841: Logically Collective
843: Input Parameters:
844: + ts - timestepping context
845: - gleetype - type of `TSGLEE` scheme
847: Level: intermediate
849: .seealso: [](ch_ts), `TSGLEEGetType()`, `TSGLEE`
850: @*/
851: PetscErrorCode TSGLEESetType(TS ts, TSGLEEType gleetype)
852: {
853: PetscFunctionBegin;
855: PetscAssertPointer(gleetype, 2);
856: PetscTryMethod(ts, "TSGLEESetType_C", (TS, TSGLEEType), (ts, gleetype));
857: PetscFunctionReturn(PETSC_SUCCESS);
858: }
860: /*@
861: TSGLEEGetType - Get the type of `TSGLEE` scheme
863: Logically Collective
865: Input Parameter:
866: . ts - timestepping context
868: Output Parameter:
869: . gleetype - type of `TSGLEE` scheme
871: Level: intermediate
873: .seealso: [](ch_ts), `TSGLEE`, `TSGLEESetType()`
874: @*/
875: PetscErrorCode TSGLEEGetType(TS ts, TSGLEEType *gleetype)
876: {
877: PetscFunctionBegin;
879: PetscUseMethod(ts, "TSGLEEGetType_C", (TS, TSGLEEType *), (ts, gleetype));
880: PetscFunctionReturn(PETSC_SUCCESS);
881: }
883: static PetscErrorCode TSGLEEGetType_GLEE(TS ts, TSGLEEType *gleetype)
884: {
885: TS_GLEE *glee = (TS_GLEE *)ts->data;
887: PetscFunctionBegin;
888: if (!glee->tableau) PetscCall(TSGLEESetType(ts, TSGLEEDefaultType));
889: *gleetype = glee->tableau->name;
890: PetscFunctionReturn(PETSC_SUCCESS);
891: }
892: static PetscErrorCode TSGLEESetType_GLEE(TS ts, TSGLEEType gleetype)
893: {
894: TS_GLEE *glee = (TS_GLEE *)ts->data;
895: PetscBool match;
896: GLEETableauLink link;
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: }
914: static PetscErrorCode TSGetStages_GLEE(TS ts, PetscInt *ns, Vec **Y)
915: {
916: TS_GLEE *glee = (TS_GLEE *)ts->data;
918: PetscFunctionBegin;
919: if (ns) *ns = glee->tableau->s;
920: if (Y) *Y = glee->YStage;
921: PetscFunctionReturn(PETSC_SUCCESS);
922: }
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;
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: }
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;
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: }
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;
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: }
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;
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: }
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: }
1008: /*MC
1009: TSGLEE - ODE and DAE solver using General Linear with Error Estimation schemes
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()`
1014: Level: beginner
1016: Notes:
1017: The default is `TSGLEE35`, it can be changed with `TSGLEESetType()` or `-ts_glee_type type`
1019: The only implicit scheme is `TSGLEEi1`
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*/
1025: PETSC_EXTERN PetscErrorCode TSCreate_GLEE(TS ts)
1026: {
1027: TS_GLEE *th;
1029: PetscFunctionBegin;
1030: PetscCall(TSGLEEInitializePackage());
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;
1051: ts->usessnes = PETSC_TRUE;
1053: PetscCall(PetscNew(&th));
1054: ts->data = (void *)th;
1056: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSGLEEGetType_C", TSGLEEGetType_GLEE));
1057: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSGLEESetType_C", TSGLEESetType_GLEE));
1058: PetscFunctionReturn(PETSC_SUCCESS);
1059: }