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