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: {
147: if (TSGLEERegisterAllCalled) return 0;
148: TSGLEERegisterAllCalled = PETSC_TRUE;
150: {
151: #define GAMMA 0.5
152: /* y-eps form */
153: const PetscInt
154: p = 1,
155: s = 3,
156: r = 2;
157: const PetscReal
158: A[3][3] = {{1.0,0,0},{0,0.5,0},{0,0.5,0.5}},
159: B[2][3] = {{1.0,0,0},{-2.0,1.0,1.0}},
160: U[3][2] = {{1.0,0},{1.0,0.5},{1.0,0.5}},
161: V[2][2] = {{1,0},{0,1}},
162: S[2] = {1,0},
163: F[2] = {1,0},
164: Fembed[2] = {1,1-GAMMA},
165: Ferror[2] = {0,1},
166: Serror[2] = {1,0};
167: 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);
168: }
169: {
170: #undef GAMMA
171: #define GAMMA 0.0
172: /* y-eps form */
173: const PetscInt
174: p = 2,
175: s = 3,
176: r = 2;
177: const PetscReal
178: A[3][3] = {{0,0,0},{1,0,0},{0.25,0.25,0}},
179: 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}},
180: U[3][2] = {{1,0},{1,10},{1,-1}},
181: V[2][2] = {{1,0},{0,1}},
182: S[2] = {1,0},
183: F[2] = {1,0},
184: Fembed[2] = {1,1-GAMMA},
185: Ferror[2] = {0,1},
186: Serror[2] = {1,0};
187: 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);
188: }
189: {
190: #undef GAMMA
191: #define GAMMA 0.0
192: /* y-y~ form */
193: const PetscInt
194: p = 2,
195: s = 4,
196: r = 2;
197: const PetscReal
198: 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}},
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}},
200: U[4][2] = {{0,1},{75.0/58.0,-17.0/58.0},{0,1},{0,1}},
201: V[2][2] = {{1,0},{0,1}},
202: S[2] = {1,1},
203: F[2] = {1,0},
204: Fembed[2] = {0,1},
205: Ferror[2] = {-1.0/(1.0-GAMMA),1.0/(1.0-GAMMA)},
206: Serror[2] = {1.0-GAMMA,1.0};
207: 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);
208: }
209: {
210: #undef GAMMA
211: #define GAMMA 0.0
212: /* y-y~ form */
213: const PetscInt
214: p = 2,
215: s = 5,
216: r = 2;
217: const PetscReal
218: A[5][5] = {{0,0,0,0,0},
219: {-0.94079244066783383269,0,0,0,0},
220: {0.64228187778301907108,0.10915356933958500042,0,0,0},
221: {-0.51764297742287450812,0.74414270351096040738,-0.71404164927824538121,0,0},
222: {-0.44696561556825969206,-0.76768425657590196518,0.20111608138142987881,0.93828186737840469796,0}},
223: B[2][5] = {{-0.029309178948150356153,-0.49671981884013874923,0.34275801517650053274,0.32941112623949194988,0.85385985637229662276},
224: {0.78133219686062535272,0.074238691892675897635,0.57957363498384957966,-0.24638502829674959968,-0.18875949544040123033}},
225: U[5][2] = {{0.16911424754448327735,0.83088575245551672265},
226: {0.53638465733199574340,0.46361534266800425660},
227: {0.39901579167169582526,0.60098420832830417474},
228: {0.87689005530618575480,0.12310994469381424520},
229: {0.99056100455550913009,0.0094389954444908699092}},
230: V[2][2] = {{1,0},{0,1}},
231: S[2] = {1,1},
232: F[2] = {1,0},
233: Fembed[2] = {0,1},
234: Ferror[2] = {-1.0/(1.0-GAMMA),1.0/(1.0-GAMMA)},
235: Serror[2] = {1.0-GAMMA,1.0};
236: 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);
237: }
238: {
239: #undef GAMMA
240: #define GAMMA 0.0
241: /* y-y~ form */
242: const PetscInt
243: p = 3,
244: s = 5,
245: r = 2;
246: const PetscReal
247: A[5][5] = {{0,0,0,0,0},
248: {- 2169604947363702313.0 / 24313474998937147335.0,0,0,0,0},
249: {46526746497697123895.0 / 94116917485856474137.0,-10297879244026594958.0 / 49199457603717988219.0,0,0,0},
250: {23364788935845982499.0 / 87425311444725389446.0,-79205144337496116638.0 / 148994349441340815519.0,40051189859317443782.0 / 36487615018004984309.0,0,0},
251: {42089522664062539205.0 / 124911313006412840286.0,-15074384760342762939.0 / 137927286865289746282.0,-62274678522253371016.0 / 125918573676298591413.0,13755475729852471739.0 / 79257927066651693390.0,0}},
252: B[2][5] = {{61546696837458703723.0 / 56982519523786160813.0,-55810892792806293355.0 / 206957624151308356511.0,24061048952676379087.0 / 158739347956038723465.0,3577972206874351339.0 / 7599733370677197135.0,-59449832954780563947.0 / 137360038685338563670.0},
253: {- 9738262186984159168.0 / 99299082461487742983.0,-32797097931948613195.0 / 61521565616362163366.0,42895514606418420631.0 / 71714201188501437336.0,22608567633166065068.0 / 55371917805607957003.0,94655809487476459565.0 / 151517167160302729021.0}},
254: U[5][2] = {{70820309139834661559.0 / 80863923579509469826.0,10043614439674808267.0 / 80863923579509469826.0},
255: {161694774978034105510.0 / 106187653640211060371.0,-55507121337823045139.0 / 106187653640211060371.0},
256: {78486094644566264568.0 / 88171030896733822981.0,9684936252167558413.0 / 88171030896733822981.0},
257: {65394922146334854435.0 / 84570853840405479554.0,19175931694070625119.0 / 84570853840405479554.0},
258: {8607282770183754108.0 / 108658046436496925911.0,100050763666313171803.0 / 108658046436496925911.0}},
259: V[2][2] = {{1,0},{0,1}},
260: S[2] = {1,1},
261: F[2] = {1,0},
262: Fembed[2] = {0,1},
263: Ferror[2] = {-1.0/(1.0-GAMMA),1.0/(1.0-GAMMA)},
264: Serror[2] = {1.0-GAMMA,1.0};
265: 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);
266: }
267: {
268: #undef GAMMA
269: #define GAMMA 0.25
270: /* y-eps form */
271: const PetscInt
272: p = 2,
273: s = 6,
274: r = 2;
275: const PetscReal
276: A[6][6] = {{0,0,0,0,0,0},
277: {1,0,0,0,0,0},
278: {0,0,0,0,0,0},
279: {0,0,0.5,0,0,0},
280: {0,0,0.25,0.25,0,0},
281: {0,0,0.25,0.25,0.5,0}},
282: B[2][6] = {{0.5,0.5,0,0,0,0},
283: {-2.0/3.0,-2.0/3.0,1.0/3.0,1.0/3.0,1.0/3.0,1.0/3.0}},
284: U[6][2] = {{1,0},{1,0},{1,0.75},{1,0.75},{1,0.75},{1,0.75}},
285: V[2][2] = {{1,0},{0,1}},
286: S[2] = {1,0},
287: F[2] = {1,0},
288: Fembed[2] = {1,1-GAMMA},
289: Ferror[2] = {0,1},
290: Serror[2] = {1,0};
291: 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);
292: }
293: {
294: #undef GAMMA
295: #define GAMMA 0.0
296: /* y-eps form */
297: const PetscInt
298: p = 3,
299: s = 8,
300: r = 2;
301: const PetscReal
302: A[8][8] = {{0,0,0,0,0,0,0,0},
303: {0.5,0,0,0,0,0,0,0},
304: {-1,2,0,0,0,0,0,0},
305: {1.0/6.0,2.0/3.0,1.0/6.0,0,0,0,0,0},
306: {0,0,0,0,0,0,0,0},
307: {-7.0/24.0,1.0/3.0,1.0/12.0,-1.0/8.0,0.5,0,0,0},
308: {7.0/6.0,-4.0/3.0,-1.0/3.0,0.5,-1.0,2.0,0,0},
309: {0,0,0,0,1.0/6.0,2.0/3.0,1.0/6.0,0}},
310: B[2][8] = {{1.0/6.0,2.0/3.0,1.0/6.0,0,0,0,0,0},
311: {-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}},
312: U[8][2] = {{1,0},{1,0},{1,0},{1,0},{1,1},{1,1},{1,1},{1,1}},
313: V[2][2] = {{1,0},{0,1}},
314: S[2] = {1,0},
315: F[2] = {1,0},
316: Fembed[2] = {1,1-GAMMA},
317: Ferror[2] = {0,1},
318: Serror[2] = {1,0};
319: 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);
320: }
321: {
322: #undef GAMMA
323: #define GAMMA 0.25
324: /* y-eps form */
325: const PetscInt
326: p = 2,
327: s = 9,
328: r = 2;
329: const PetscReal
330: A[9][9] = {{0,0,0,0,0,0,0,0,0},
331: {0.585786437626904966,0,0,0,0,0,0,0,0},
332: {0.149999999999999994,0.849999999999999978,0,0,0,0,0,0,0},
333: {0,0,0,0,0,0,0,0,0},
334: {0,0,0,0.292893218813452483,0,0,0,0,0},
335: {0,0,0,0.0749999999999999972,0.424999999999999989,0,0,0,0},
336: {0,0,0,0.176776695296636893,0.176776695296636893,0.146446609406726241,0,0,0},
337: {0,0,0,0.176776695296636893,0.176776695296636893,0.146446609406726241,0.292893218813452483,0,0},
338: {0,0,0,0.176776695296636893,0.176776695296636893,0.146446609406726241,0.0749999999999999972,0.424999999999999989,0}},
339: B[2][9] = {{0.353553390593273786,0.353553390593273786,0.292893218813452483,0,0,0,0,0,0},
340: {-0.471404520791031678,-0.471404520791031678,-0.390524291751269959,0.235702260395515839,0.235702260395515839,0.195262145875634979,0.235702260395515839,0.235702260395515839,0.195262145875634979}},
341: 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}},
342: V[2][2] = {{1,0},{0,1}},
343: S[2] = {1,0},
344: F[2] = {1,0},
345: Fembed[2] = {1,1-GAMMA},
346: Ferror[2] = {0,1},
347: Serror[2] = {1,0};
348: 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);
349: }
351: return 0;
352: }
354: /*@C
355: TSGLEERegisterDestroy - Frees the list of schemes that were registered by TSGLEERegister().
357: Not Collective
359: Level: advanced
361: .seealso: TSGLEERegister(), TSGLEERegisterAll()
362: @*/
363: PetscErrorCode TSGLEERegisterDestroy(void)
364: {
366: GLEETableauLink link;
368: while ((link = GLEETableauList)) {
369: GLEETableau t = &link->tab;
370: GLEETableauList = link->next;
371: PetscFree5(t->A,t->B,t->U,t->V,t->c);
372: PetscFree2(t->S,t->F); ierr;
373: PetscFree (t->Fembed); ierr;
374: PetscFree (t->Ferror); ierr;
375: PetscFree (t->Serror); ierr;
376: PetscFree (t->binterp); ierr;
377: PetscFree (t->name); ierr;
378: PetscFree (link); ierr;
379: }
380: TSGLEERegisterAllCalled = PETSC_FALSE;
381: return 0;
382: }
384: /*@C
385: TSGLEEInitializePackage - This function initializes everything in the TSGLEE package. It is called
386: from TSInitializePackage().
388: Level: developer
390: .seealso: PetscInitialize()
391: @*/
392: PetscErrorCode TSGLEEInitializePackage(void)
393: {
394: if (TSGLEEPackageInitialized) return 0;
395: TSGLEEPackageInitialized = PETSC_TRUE;
396: TSGLEERegisterAll();
397: PetscObjectComposedDataRegister(&explicit_stage_time_id);
398: PetscRegisterFinalize(TSGLEEFinalizePackage);
399: return 0;
400: }
402: /*@C
403: TSGLEEFinalizePackage - This function destroys everything in the TSGLEE package. It is
404: called from PetscFinalize().
406: Level: developer
408: .seealso: PetscFinalize()
409: @*/
410: PetscErrorCode TSGLEEFinalizePackage(void)
411: {
412: TSGLEEPackageInitialized = PETSC_FALSE;
413: TSGLEERegisterDestroy();
414: return 0;
415: }
417: /*@C
418: TSGLEERegister - register an GLEE scheme by providing the entries in the Butcher tableau
420: Not Collective, but the same schemes should be registered on all processes on which they will be used
422: Input Parameters:
423: + name - identifier for method
424: . order - order of method
425: . s - number of stages
426: . r - number of steps
427: . gamma - LTE ratio
428: . A - stage coefficients (dimension s*s, row-major)
429: . B - step completion coefficients (dimension r*s, row-major)
430: . U - method coefficients (dimension s*r, row-major)
431: . V - method coefficients (dimension r*r, row-major)
432: . S - starting coefficients
433: . F - finishing coefficients
434: . c - abscissa (dimension s; NULL to use row sums of A)
435: . Fembed - step completion coefficients for embedded method
436: . Ferror - error computation coefficients
437: . Serror - error initialization coefficients
438: . pinterp - order of interpolation (0 if unavailable)
439: - binterp - array of interpolation coefficients (NULL if unavailable)
441: Notes:
442: Several GLEE methods are provided, this function is only needed to create new methods.
444: Level: advanced
446: .seealso: TSGLEE
447: @*/
448: PetscErrorCode TSGLEERegister(TSGLEEType name,PetscInt order,PetscInt s, PetscInt r,
449: PetscReal gamma,
450: const PetscReal A[],const PetscReal B[],
451: const PetscReal U[],const PetscReal V[],
452: const PetscReal S[],const PetscReal F[],
453: const PetscReal c[],
454: const PetscReal Fembed[],const PetscReal Ferror[],
455: const PetscReal Serror[],
456: PetscInt pinterp, const PetscReal binterp[])
457: {
458: GLEETableauLink link;
459: GLEETableau t;
460: PetscInt i,j;
462: TSGLEEInitializePackage();
463: PetscNew(&link);
464: t = &link->tab;
465: PetscStrallocpy(name,&t->name);
466: t->order = order;
467: t->s = s;
468: t->r = r;
469: t->gamma = gamma;
470: PetscMalloc5(s*s,&t->A,r*r,&t->V,s,&t->c,r*s,&t->B,s*r,&t->U);
471: PetscMalloc2(r,&t->S,r,&t->F);
472: PetscArraycpy(t->A,A,s*s);
473: PetscArraycpy(t->B,B,r*s);
474: PetscArraycpy(t->U,U,s*r);
475: PetscArraycpy(t->V,V,r*r);
476: PetscArraycpy(t->S,S,r);
477: PetscArraycpy(t->F,F,r);
478: if (c) {
479: PetscArraycpy(t->c,c,s);
480: } else {
481: for (i=0; i<s; i++) for (j=0,t->c[i]=0; j<s; j++) t->c[i] += A[i*s+j];
482: }
483: PetscMalloc1(r,&t->Fembed);
484: PetscMalloc1(r,&t->Ferror);
485: PetscMalloc1(r,&t->Serror);
486: PetscArraycpy(t->Fembed,Fembed,r);
487: PetscArraycpy(t->Ferror,Ferror,r);
488: PetscArraycpy(t->Serror,Serror,r);
489: t->pinterp = pinterp;
490: PetscMalloc1(s*pinterp,&t->binterp);
491: PetscArraycpy(t->binterp,binterp,s*pinterp);
493: link->next = GLEETableauList;
494: GLEETableauList = link;
495: return 0;
496: }
498: static PetscErrorCode TSEvaluateStep_GLEE(TS ts,PetscInt order,Vec X,PetscBool *done)
499: {
500: TS_GLEE *glee = (TS_GLEE*) ts->data;
501: GLEETableau tab = glee->tableau;
502: PetscReal h, *B = tab->B, *V = tab->V,
503: *F = tab->F,
504: *Fembed = tab->Fembed;
505: PetscInt s = tab->s, r = tab->r, i, j;
506: Vec *Y = glee->Y, *YdotStage = glee->YdotStage;
507: PetscScalar *ws = glee->swork, *wr = glee->rwork;
510: switch (glee->status) {
511: case TS_STEP_INCOMPLETE:
512: case TS_STEP_PENDING:
513: h = ts->time_step; break;
514: case TS_STEP_COMPLETE:
515: h = ts->ptime - ts->ptime_prev; break;
516: default: SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_PLIB,"Invalid TSStepStatus");
517: }
519: if (order == tab->order) {
521: /* Note: Irrespective of whether status is TS_STEP_INCOMPLETE
522: or TS_STEP_COMPLETE, glee->X has the solution at the
523: beginning of the time step. So no need to roll-back.
524: */
525: if (glee->status == TS_STEP_INCOMPLETE) {
526: for (i=0; i<r; i++) {
527: VecZeroEntries(Y[i]);
528: for (j=0; j<r; j++) wr[j] = V[i*r+j];
529: VecMAXPY(Y[i],r,wr,glee->X);
530: for (j=0; j<s; j++) ws[j] = h*B[i*s+j];
531: VecMAXPY(Y[i],s,ws,YdotStage);
532: }
533: VecZeroEntries(X);
534: for (j=0; j<r; j++) wr[j] = F[j];
535: VecMAXPY(X,r,wr,Y);
536: } else VecCopy(ts->vec_sol,X);
537: return 0;
539: } else if (order == tab->order-1) {
541: /* Complete with the embedded method (Fembed) */
542: for (i=0; i<r; i++) {
543: VecZeroEntries(Y[i]);
544: for (j=0; j<r; j++) wr[j] = V[i*r+j];
545: VecMAXPY(Y[i],r,wr,glee->X);
546: for (j=0; j<s; j++) ws[j] = h*B[i*s+j];
547: VecMAXPY(Y[i],s,ws,YdotStage);
548: }
549: VecZeroEntries(X);
550: for (j=0; j<r; j++) wr[j] = Fembed[j];
551: VecMAXPY(X,r,wr,Y);
553: if (done) *done = PETSC_TRUE;
554: return 0;
555: }
556: if (done) *done = PETSC_FALSE;
557: else SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"GLEE '%s' of order %D cannot evaluate step at order %D",tab->name,tab->order,order);
558: return 0;
559: }
561: static PetscErrorCode TSStep_GLEE(TS ts)
562: {
563: TS_GLEE *glee = (TS_GLEE*)ts->data;
564: GLEETableau tab = glee->tableau;
565: const PetscInt s = tab->s, r = tab->r;
566: PetscReal *A = tab->A, *U = tab->U,
567: *F = tab->F,
568: *c = tab->c;
569: Vec *Y = glee->Y, *X = glee->X,
570: *YStage = glee->YStage,
571: *YdotStage = glee->YdotStage,
572: W = glee->W;
573: SNES snes;
574: PetscScalar *ws = glee->swork, *wr = glee->rwork;
575: TSAdapt adapt;
576: PetscInt i,j,reject,next_scheme,its,lits;
577: PetscReal next_time_step;
578: PetscReal t;
579: PetscBool accept;
581: PetscCitationsRegister(citation,&cited);
583: for (i=0; i<r; i++) VecCopy(Y[i],X[i]);
585: TSGetSNES(ts,&snes);
586: next_time_step = ts->time_step;
587: t = ts->ptime;
588: accept = PETSC_TRUE;
589: glee->status = TS_STEP_INCOMPLETE;
591: for (reject=0; reject<ts->max_reject && !ts->reason; reject++,ts->reject++) {
593: PetscReal h = ts->time_step;
594: TSPreStep(ts);
596: for (i=0; i<s; i++) {
598: glee->stage_time = t + h*c[i];
599: TSPreStage(ts,glee->stage_time);
601: if (A[i*s+i] == 0) { /* Explicit stage */
602: VecZeroEntries(YStage[i]);
603: for (j=0; j<r; j++) wr[j] = U[i*r+j];
604: VecMAXPY(YStage[i],r,wr,X);
605: for (j=0; j<i; j++) ws[j] = h*A[i*s+j];
606: VecMAXPY(YStage[i],i,ws,YdotStage);
607: } else { /* Implicit stage */
608: glee->scoeff = 1.0/A[i*s+i];
609: /* compute right-hand-side */
610: VecZeroEntries(W);
611: for (j=0; j<r; j++) wr[j] = U[i*r+j];
612: VecMAXPY(W,r,wr,X);
613: for (j=0; j<i; j++) ws[j] = h*A[i*s+j];
614: VecMAXPY(W,i,ws,YdotStage);
615: VecScale(W,glee->scoeff/h);
616: /* set initial guess */
617: VecCopy(i>0 ? YStage[i-1] : ts->vec_sol,YStage[i]);
618: /* solve for this stage */
619: SNESSolve(snes,W,YStage[i]);
620: SNESGetIterationNumber(snes,&its);
621: SNESGetLinearSolveIterations(snes,&lits);
622: ts->snes_its += its; ts->ksp_its += lits;
623: }
624: TSGetAdapt(ts,&adapt);
625: TSAdaptCheckStage(adapt,ts,glee->stage_time,YStage[i],&accept);
626: if (!accept) goto reject_step;
627: TSPostStage(ts,glee->stage_time,i,YStage);
628: TSComputeRHSFunction(ts,t+h*c[i],YStage[i],YdotStage[i]);
629: }
630: TSEvaluateStep(ts,tab->order,ts->vec_sol,NULL);
631: glee->status = TS_STEP_PENDING;
633: /* Register only the current method as a candidate because we're not supporting multiple candidates yet. */
634: TSGetAdapt(ts,&adapt);
635: TSAdaptCandidatesClear(adapt);
636: TSAdaptCandidateAdd(adapt,tab->name,tab->order,1,tab->ccfl,(PetscReal)tab->s,PETSC_TRUE);
637: TSAdaptChoose(adapt,ts,ts->time_step,&next_scheme,&next_time_step,&accept);
638: if (accept) {
639: /* ignore next_scheme for now */
640: ts->ptime += ts->time_step;
641: ts->time_step = next_time_step;
642: glee->status = TS_STEP_COMPLETE;
643: /* compute and store the global error */
644: /* Note: this is not needed if TSAdaptGLEE is not used */
645: TSGetTimeError(ts,0,&(glee->yGErr));
646: PetscObjectComposedDataSetReal((PetscObject)ts->vec_sol,explicit_stage_time_id,ts->ptime);
647: break;
648: } else { /* Roll back the current step */
649: for (j=0; j<r; j++) wr[j] = F[j];
650: VecMAXPY(ts->vec_sol,r,wr,X);
651: ts->time_step = next_time_step;
652: glee->status = TS_STEP_INCOMPLETE;
653: }
654: reject_step: continue;
655: }
656: if (glee->status != TS_STEP_COMPLETE && !ts->reason) ts->reason = TS_DIVERGED_STEP_REJECTED;
657: return 0;
658: }
660: static PetscErrorCode TSInterpolate_GLEE(TS ts,PetscReal itime,Vec X)
661: {
662: TS_GLEE *glee = (TS_GLEE*)ts->data;
663: PetscInt s=glee->tableau->s, pinterp=glee->tableau->pinterp,i,j;
664: PetscReal h,tt,t;
665: PetscScalar *b;
666: const PetscReal *B = glee->tableau->binterp;
669: switch (glee->status) {
670: case TS_STEP_INCOMPLETE:
671: case TS_STEP_PENDING:
672: h = ts->time_step;
673: t = (itime - ts->ptime)/h;
674: break;
675: case TS_STEP_COMPLETE:
676: h = ts->ptime - ts->ptime_prev;
677: t = (itime - ts->ptime)/h + 1; /* In the interval [0,1] */
678: break;
679: default: SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_PLIB,"Invalid TSStepStatus");
680: }
681: PetscMalloc1(s,&b);
682: for (i=0; i<s; i++) b[i] = 0;
683: for (j=0,tt=t; j<pinterp; j++,tt*=t) {
684: for (i=0; i<s; i++) {
685: b[i] += h * B[i*pinterp+j] * tt;
686: }
687: }
688: VecCopy(glee->YStage[0],X);
689: VecMAXPY(X,s,b,glee->YdotStage);
690: PetscFree(b);
691: return 0;
692: }
694: /*------------------------------------------------------------*/
695: static PetscErrorCode TSReset_GLEE(TS ts)
696: {
697: TS_GLEE *glee = (TS_GLEE*)ts->data;
698: PetscInt s, r;
700: if (!glee->tableau) return 0;
701: s = glee->tableau->s;
702: r = glee->tableau->r;
703: VecDestroyVecs(r,&glee->Y);
704: VecDestroyVecs(r,&glee->X);
705: VecDestroyVecs(s,&glee->YStage);
706: VecDestroyVecs(s,&glee->YdotStage);
707: VecDestroy(&glee->Ydot);
708: VecDestroy(&glee->yGErr);
709: VecDestroy(&glee->W);
710: PetscFree2(glee->swork,glee->rwork);
711: return 0;
712: }
714: static PetscErrorCode TSGLEEGetVecs(TS ts,DM dm,Vec *Ydot)
715: {
716: TS_GLEE *glee = (TS_GLEE*)ts->data;
718: if (Ydot) {
719: if (dm && dm != ts->dm) {
720: DMGetNamedGlobalVector(dm,"TSGLEE_Ydot",Ydot);
721: } else *Ydot = glee->Ydot;
722: }
723: return 0;
724: }
726: static PetscErrorCode TSGLEERestoreVecs(TS ts,DM dm,Vec *Ydot)
727: {
728: if (Ydot) {
729: if (dm && dm != ts->dm) {
730: DMRestoreNamedGlobalVector(dm,"TSGLEE_Ydot",Ydot);
731: }
732: }
733: return 0;
734: }
736: /*
737: This defines the nonlinear equation that is to be solved with SNES
738: */
739: static PetscErrorCode SNESTSFormFunction_GLEE(SNES snes,Vec X,Vec F,TS ts)
740: {
741: TS_GLEE *glee = (TS_GLEE*)ts->data;
742: DM dm,dmsave;
743: Vec Ydot;
744: PetscReal shift = glee->scoeff / ts->time_step;
746: SNESGetDM(snes,&dm);
747: TSGLEEGetVecs(ts,dm,&Ydot);
748: /* Set Ydot = shift*X */
749: VecCopy(X,Ydot);
750: VecScale(Ydot,shift);
751: dmsave = ts->dm;
752: ts->dm = dm;
754: TSComputeIFunction(ts,glee->stage_time,X,Ydot,F,PETSC_FALSE);
756: ts->dm = dmsave;
757: TSGLEERestoreVecs(ts,dm,&Ydot);
758: return 0;
759: }
761: static PetscErrorCode SNESTSFormJacobian_GLEE(SNES snes,Vec X,Mat A,Mat B,TS ts)
762: {
763: TS_GLEE *glee = (TS_GLEE*)ts->data;
764: DM dm,dmsave;
765: Vec Ydot;
766: PetscReal shift = glee->scoeff / ts->time_step;
768: SNESGetDM(snes,&dm);
769: TSGLEEGetVecs(ts,dm,&Ydot);
770: /* glee->Ydot has already been computed in SNESTSFormFunction_GLEE (SNES guarantees this) */
771: dmsave = ts->dm;
772: ts->dm = dm;
774: TSComputeIJacobian(ts,glee->stage_time,X,Ydot,shift,A,B,PETSC_FALSE);
776: ts->dm = dmsave;
777: TSGLEERestoreVecs(ts,dm,&Ydot);
778: return 0;
779: }
781: static PetscErrorCode DMCoarsenHook_TSGLEE(DM fine,DM coarse,void *ctx)
782: {
783: return 0;
784: }
786: static PetscErrorCode DMRestrictHook_TSGLEE(DM fine,Mat restrct,Vec rscale,Mat inject,DM coarse,void *ctx)
787: {
788: return 0;
789: }
791: static PetscErrorCode DMSubDomainHook_TSGLEE(DM dm,DM subdm,void *ctx)
792: {
793: return 0;
794: }
796: static PetscErrorCode DMSubDomainRestrictHook_TSGLEE(DM dm,VecScatter gscat,VecScatter lscat,DM subdm,void *ctx)
797: {
798: return 0;
799: }
801: static PetscErrorCode TSSetUp_GLEE(TS ts)
802: {
803: TS_GLEE *glee = (TS_GLEE*)ts->data;
804: GLEETableau tab;
805: PetscInt s,r;
806: DM dm;
808: if (!glee->tableau) {
809: TSGLEESetType(ts,TSGLEEDefaultType);
810: }
811: tab = glee->tableau;
812: s = tab->s;
813: r = tab->r;
814: VecDuplicateVecs(ts->vec_sol,r,&glee->Y);
815: VecDuplicateVecs(ts->vec_sol,r,&glee->X);
816: VecDuplicateVecs(ts->vec_sol,s,&glee->YStage);
817: VecDuplicateVecs(ts->vec_sol,s,&glee->YdotStage);
818: VecDuplicate(ts->vec_sol,&glee->Ydot);
819: VecDuplicate(ts->vec_sol,&glee->yGErr);
820: VecZeroEntries(glee->yGErr);
821: VecDuplicate(ts->vec_sol,&glee->W);
822: PetscMalloc2(s,&glee->swork,r,&glee->rwork);
823: TSGetDM(ts,&dm);
824: DMCoarsenHookAdd(dm,DMCoarsenHook_TSGLEE,DMRestrictHook_TSGLEE,ts);
825: DMSubDomainHookAdd(dm,DMSubDomainHook_TSGLEE,DMSubDomainRestrictHook_TSGLEE,ts);
826: return 0;
827: }
829: PetscErrorCode TSStartingMethod_GLEE(TS ts)
830: {
831: TS_GLEE *glee = (TS_GLEE*)ts->data;
832: GLEETableau tab = glee->tableau;
833: PetscInt r=tab->r,i;
834: PetscReal *S=tab->S;
836: for (i=0; i<r; i++) {
837: VecZeroEntries(glee->Y[i]);
838: VecAXPY(glee->Y[i],S[i],ts->vec_sol);
839: }
841: return 0;
842: }
844: /*------------------------------------------------------------*/
846: static PetscErrorCode TSSetFromOptions_GLEE(PetscOptionItems *PetscOptionsObject,TS ts)
847: {
848: char gleetype[256];
850: PetscOptionsHead(PetscOptionsObject,"GLEE ODE solver options");
851: {
852: GLEETableauLink link;
853: PetscInt count,choice;
854: PetscBool flg;
855: const char **namelist;
857: PetscStrncpy(gleetype,TSGLEEDefaultType,sizeof(gleetype));
858: for (link=GLEETableauList,count=0; link; link=link->next,count++) ;
859: PetscMalloc1(count,(char***)&namelist);
860: for (link=GLEETableauList,count=0; link; link=link->next,count++) namelist[count] = link->tab.name;
861: PetscOptionsEList("-ts_glee_type","Family of GLEE method","TSGLEESetType",(const char*const*)namelist,count,gleetype,&choice,&flg);
862: TSGLEESetType(ts,flg ? namelist[choice] : gleetype);
863: PetscFree(namelist);
864: }
865: PetscOptionsTail();
866: return 0;
867: }
869: static PetscErrorCode TSView_GLEE(TS ts,PetscViewer viewer)
870: {
871: TS_GLEE *glee = (TS_GLEE*)ts->data;
872: GLEETableau tab = glee->tableau;
873: PetscBool iascii;
875: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
876: if (iascii) {
877: TSGLEEType gleetype;
878: char buf[512];
879: TSGLEEGetType(ts,&gleetype);
880: PetscViewerASCIIPrintf(viewer," GLEE type %s\n",gleetype);
881: PetscFormatRealArray(buf,sizeof(buf),"% 8.6f",tab->s,tab->c);
882: PetscViewerASCIIPrintf(viewer," Abscissa c = %s\n",buf);
883: /* Note: print out r as well */
884: }
885: return 0;
886: }
888: static PetscErrorCode TSLoad_GLEE(TS ts,PetscViewer viewer)
889: {
890: SNES snes;
891: TSAdapt tsadapt;
893: TSGetAdapt(ts,&tsadapt);
894: TSAdaptLoad(tsadapt,viewer);
895: TSGetSNES(ts,&snes);
896: SNESLoad(snes,viewer);
897: /* function and Jacobian context for SNES when used with TS is always ts object */
898: SNESSetFunction(snes,NULL,NULL,ts);
899: SNESSetJacobian(snes,NULL,NULL,NULL,ts);
900: return 0;
901: }
903: /*@C
904: TSGLEESetType - Set the type of GLEE scheme
906: Logically collective
908: Input Parameters:
909: + ts - timestepping context
910: - gleetype - type of GLEE-scheme
912: Level: intermediate
914: .seealso: TSGLEEGetType(), TSGLEE
915: @*/
916: PetscErrorCode TSGLEESetType(TS ts,TSGLEEType gleetype)
917: {
920: PetscTryMethod(ts,"TSGLEESetType_C",(TS,TSGLEEType),(ts,gleetype));
921: return 0;
922: }
924: /*@C
925: TSGLEEGetType - Get the type of GLEE scheme
927: Logically collective
929: Input Parameter:
930: . ts - timestepping context
932: Output Parameter:
933: . gleetype - type of GLEE-scheme
935: Level: intermediate
937: .seealso: TSGLEESetType()
938: @*/
939: PetscErrorCode TSGLEEGetType(TS ts,TSGLEEType *gleetype)
940: {
942: PetscUseMethod(ts,"TSGLEEGetType_C",(TS,TSGLEEType*),(ts,gleetype));
943: return 0;
944: }
946: PetscErrorCode TSGLEEGetType_GLEE(TS ts,TSGLEEType *gleetype)
947: {
948: TS_GLEE *glee = (TS_GLEE*)ts->data;
950: if (!glee->tableau) {
951: TSGLEESetType(ts,TSGLEEDefaultType);
952: }
953: *gleetype = glee->tableau->name;
954: return 0;
955: }
956: PetscErrorCode TSGLEESetType_GLEE(TS ts,TSGLEEType gleetype)
957: {
958: TS_GLEE *glee = (TS_GLEE*)ts->data;
959: PetscBool match;
960: GLEETableauLink link;
962: if (glee->tableau) {
963: PetscStrcmp(glee->tableau->name,gleetype,&match);
964: if (match) return 0;
965: }
966: for (link = GLEETableauList; link; link=link->next) {
967: PetscStrcmp(link->tab.name,gleetype,&match);
968: if (match) {
969: TSReset_GLEE(ts);
970: glee->tableau = &link->tab;
971: return 0;
972: }
973: }
974: SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_UNKNOWN_TYPE,"Could not find '%s'",gleetype);
975: }
977: static PetscErrorCode TSGetStages_GLEE(TS ts,PetscInt *ns,Vec **Y)
978: {
979: TS_GLEE *glee = (TS_GLEE*)ts->data;
981: if (ns) *ns = glee->tableau->s;
982: if (Y) *Y = glee->YStage;
983: return 0;
984: }
986: PetscErrorCode TSGetSolutionComponents_GLEE(TS ts,PetscInt *n,Vec *Y)
987: {
988: TS_GLEE *glee = (TS_GLEE*)ts->data;
989: GLEETableau tab = glee->tableau;
991: if (!Y) *n = tab->r;
992: else {
993: if ((*n >= 0) && (*n < tab->r)) {
994: VecCopy(glee->Y[*n],*Y);
995: } else SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_OUTOFRANGE,"Second argument (%d) out of range[%d,%d].",*n,0,tab->r-1);
996: }
997: return 0;
998: }
1000: PetscErrorCode TSGetAuxSolution_GLEE(TS ts,Vec *X)
1001: {
1002: TS_GLEE *glee = (TS_GLEE*)ts->data;
1003: GLEETableau tab = glee->tableau;
1004: PetscReal *F = tab->Fembed;
1005: PetscInt r = tab->r;
1006: Vec *Y = glee->Y;
1007: PetscScalar *wr = glee->rwork;
1008: PetscInt i;
1010: VecZeroEntries(*X);
1011: for (i=0; i<r; i++) wr[i] = F[i];
1012: VecMAXPY((*X),r,wr,Y);
1013: return 0;
1014: }
1016: PetscErrorCode TSGetTimeError_GLEE(TS ts,PetscInt n,Vec *X)
1017: {
1018: TS_GLEE *glee = (TS_GLEE*)ts->data;
1019: GLEETableau tab = glee->tableau;
1020: PetscReal *F = tab->Ferror;
1021: PetscInt r = tab->r;
1022: Vec *Y = glee->Y;
1023: PetscScalar *wr = glee->rwork;
1024: PetscInt i;
1026: VecZeroEntries(*X);
1027: if (n==0) {
1028: for (i=0; i<r; i++) wr[i] = F[i];
1029: VecMAXPY((*X),r,wr,Y);
1030: } else if (n==-1) {
1031: *X=glee->yGErr;
1032: }
1033: return 0;
1034: }
1036: PetscErrorCode TSSetTimeError_GLEE(TS ts,Vec X)
1037: {
1038: TS_GLEE *glee = (TS_GLEE*)ts->data;
1039: GLEETableau tab = glee->tableau;
1040: PetscReal *S = tab->Serror;
1041: PetscInt r = tab->r,i;
1042: Vec *Y = glee->Y;
1045: for (i=1; i<r; i++) {
1046: VecCopy(ts->vec_sol,Y[i]);
1047: VecAXPBY(Y[i],S[0],S[1],X);
1048: VecCopy(X,glee->yGErr);
1049: }
1050: return 0;
1051: }
1053: static PetscErrorCode TSDestroy_GLEE(TS ts)
1054: {
1055: TSReset_GLEE(ts);
1056: if (ts->dm) {
1057: DMCoarsenHookRemove(ts->dm,DMCoarsenHook_TSGLEE,DMRestrictHook_TSGLEE,ts);
1058: DMSubDomainHookRemove(ts->dm,DMSubDomainHook_TSGLEE,DMSubDomainRestrictHook_TSGLEE,ts);
1059: }
1060: PetscFree(ts->data);
1061: PetscObjectComposeFunction((PetscObject)ts,"TSGLEEGetType_C",NULL);
1062: PetscObjectComposeFunction((PetscObject)ts,"TSGLEESetType_C",NULL);
1063: return 0;
1064: }
1066: /* ------------------------------------------------------------ */
1067: /*MC
1068: TSGLEE - ODE and DAE solver using General Linear with Error Estimation schemes
1070: The user should provide the right hand side of the equation
1071: using TSSetRHSFunction().
1073: Notes:
1074: The default is TSGLEE35, it can be changed with TSGLEESetType() or -ts_glee_type
1076: Level: beginner
1078: .seealso: TSCreate(), TS, TSSetType(), TSGLEESetType(), TSGLEEGetType(),
1079: TSGLEE23, TTSGLEE24, TSGLEE35, TSGLEE25I, TSGLEEEXRK2A,
1080: TSGLEERK32G1, TSGLEERK285EX, TSGLEEType, TSGLEERegister()
1082: M*/
1083: PETSC_EXTERN PetscErrorCode TSCreate_GLEE(TS ts)
1084: {
1085: TS_GLEE *th;
1087: TSGLEEInitializePackage();
1089: ts->ops->reset = TSReset_GLEE;
1090: ts->ops->destroy = TSDestroy_GLEE;
1091: ts->ops->view = TSView_GLEE;
1092: ts->ops->load = TSLoad_GLEE;
1093: ts->ops->setup = TSSetUp_GLEE;
1094: ts->ops->step = TSStep_GLEE;
1095: ts->ops->interpolate = TSInterpolate_GLEE;
1096: ts->ops->evaluatestep = TSEvaluateStep_GLEE;
1097: ts->ops->setfromoptions = TSSetFromOptions_GLEE;
1098: ts->ops->getstages = TSGetStages_GLEE;
1099: ts->ops->snesfunction = SNESTSFormFunction_GLEE;
1100: ts->ops->snesjacobian = SNESTSFormJacobian_GLEE;
1101: ts->ops->getsolutioncomponents = TSGetSolutionComponents_GLEE;
1102: ts->ops->getauxsolution = TSGetAuxSolution_GLEE;
1103: ts->ops->gettimeerror = TSGetTimeError_GLEE;
1104: ts->ops->settimeerror = TSSetTimeError_GLEE;
1105: ts->ops->startingmethod = TSStartingMethod_GLEE;
1106: ts->default_adapt_type = TSADAPTGLEE;
1108: ts->usessnes = PETSC_TRUE;
1110: PetscNewLog(ts,&th);
1111: ts->data = (void*)th;
1113: PetscObjectComposeFunction((PetscObject)ts,"TSGLEEGetType_C",TSGLEEGetType_GLEE);
1114: PetscObjectComposeFunction((PetscObject)ts,"TSGLEESetType_C",TSGLEESetType_GLEE);
1115: return 0;
1116: }